atf1.f90

薄厚通用的板单元

SUBROUTINE LOADD_2D_Composite_Plate(eload)
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 USE M5_Para_Gauss_Integral
 USE M6_Matrix_Element
 USE M8_Shape_Function
 REAL*8, DIMENSION(NELEM,MEVAB):: eload
 REAL*8, DIMENSION(NDIME,MNODE):: coren
 REAL*8, DIMENSION(:,:), ALLOCATABLE:: cores,press
 INTEGER,DIMENSION(:),   ALLOCATABLE:: noprs 
 INTEGER::idime,idist,inods,ipart,ielem
 INTEGER::nel,nic,nconp,nedge,nodps ! nodps:单元每边的结点数
 INTEGER::etp,nelef,nelel,nqpart
 REAL*8 ::q


 coren=0.0
 READ(1,*) nedge
 WRITE(3,*) ' *There are ',nedge,'   sides under in-plane distributed loadings.'
 IF(nedge>0) THEN
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'
  DO idist=1,nedge 
   READ(1,*) nel
   WRITE(3,*) '  *The No.', idist, '   distributed loading on element', nel
   SELECT CASE(ELMOD(IELES(nel,MNODE+1)))
    CASE(1:2,101)
	 nodps=2
!   CASE(2:4,62)
!    nodps=3
   END SELECT
   ALLOCATE(cores(NDIME,nodps),noprs(nodps),press(nodps,2))
   DO inods=1,nodps
    READ(1,*) noprs(inods), (press(inods,idofn),idofn=1,2)
    WRITE(3,*) noprs(inods), (press(inods,idofn),idofn=1,2)
    DO idime=1,NDIME 
     cores(idime,inods)=COORD(noprs(inods),idime)
    END DO  
    
	IF(inods==1) THEN
     DO inode=1,MNODE 
      nic=IELES(nel,inode)
      DO idime=1,NDIME
       IF(nic/=0) coren(idime,inode)=COORD(nic,idime)
      END DO  
     END DO       
     DO inode=1,MNODE   
	  IF(noprs(inods)==IELES(nel,inode)) nconp=inode
     END DO	     
  
    END IF
   END DO       

   SELECT CASE(ELMOD(IELES(nel,MNODE+1)))
    CASE(1:2,101)
	 CALL LOADD_2D_PLANE_METHOD_1(nel,nodps,cores,noprs,press,eload)
!	CASE(3:4)
!	 CALL LOADD_2D_PLANE_METHOD_2(nel,nconp,nodps,coren,noprs,press,eload)
!	CASE(81:82)
!	 CALL LOADD_2D_PLANE_METHOD_3(nel,nconp,nodps,coren,noprs,press,eload)
   END SELECT
   DEALLOCATE(cores,noprs,press)
   
  END DO
 END IF
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'

 READ(1,*) nqpart
 WRITE(3,*) ' *There are ',nqpart,'   element groups under transverse uniform loadings.'
 IF(nqpart>0) THEN
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'
  DO ipart=1,nqpart 
   READ(1,*) nelef,nelel,q
   WRITE(3,*) '  *The No.', ipart, '   transverse uniform loading from element',nelef,'   to',nelel
   WRITE(3,*) '    and the loading q=',q
   DO ielem=nelef, nelel
    etp=IELES(ielem,MNODE+1)
     SELECT CASE(ELESP(etp))
      CASE(3)
       CALL LOADD_2D_PLATE_T(ielem,q,eload)
	  CASE(4)
       CALL LOADD_2D_PLATE_Q(ielem,q,eload)
	 END SELECT
   END DO
  END DO
 END IF
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'

END SUBROUTINE LOADD_2D_Composite_Plate

!*============================================================================*!
!*S21. LOADD_2D_SINUSOIDAL_RECTANGULAR_PLATE_Q                                *!
!*         功能: 二维四边形层合板单元受双SIN分布荷载（矩形板情形时）          *!
!*============================================================================*!

SUBROUTINE LOADD_2D_SINUSOIDAL_RECTANGULAR_PLATE_Q(eload)
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 USE M4_Support_information
 USE M5_Para_Gauss_Integral
 USE M6_Matrix_Element
 USE M8_Shape_Function
 REAL*8, DIMENSION(NELEM,MEVAB):: eload
 REAL*8, DIMENSION(:),ALLOCATABLE:: rvece   
 REAL*8, DIMENSION(:,:),ALLOCATABLE:: coren   
 REAL*8, DIMENSION(NDIME,MNODE):: coren_1
 REAL*8,DIMENSION(3)::L !b,c,dl,DT
 REAL*8::DJACB,XX,YY,S,T,Pi=3.141592653589793 !,Xmin,Xmax,Ymin,Ymax
 REAL*8::q0, LA,LB,LAX,LBY,GS1,GS2,GS3,GS4,A1,B3,C3,DV
 INTEGER::kg,nic,etp
 INTEGER::ielem,idime,ievab,igaus,inode,jgaus

 ALLOCATE(rvece(MEVAB))
 rvece=0.0
 WRITE(3,*) ' *There are double sinusoidal transverse load on the rectangular plate.'
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'
 READ(1,*) q0,LA,LB,LAX,LBY
 WRITE(3,*)'** q=q0*sin[Pi*(XX+LAX)/LA]*SIN[Pi*(YY+LBY)/LB]'
 WRITE(3,*)'***q0=',q0,'    LA=',LA,'    LB=',LB,'    LAX=',LAX,'    LBY=',LBY

 DO ielem=1,NELEM
  etp=IELES(ielem,MNODE+1) 
  ALLOCATE(coren(NDIME,NNODE(etp)))
  coren=0.0
  DO inode=1,NNODE(etp)
   nic=IELES(ielem,inode)
   DO idime=1,NDIME
    IF(nic/=0) coren(idime,inode)=COORD(nic,idime)
   END DO
  END DO
  rvece=0.0 
  kg=0
   
  SELECT CASE(ELESP(etp))
   CASE(3)
!!@@@
    CALL GAUSS_triangular(NGAUS(etp))
    GS1=POSGP(1);GS2=POSGP(2);GS3=POSGP(3);GS4=POSGP(4)
    A1=WEIGP(1) ; B3=WEIGP(2) ; C3=WEIGP(3)

    IF(A1/=0.0) THEN
     L=1.0/3.0
     kg=kg+1
     CALL SHAPE_N_T3(L)
     CALL JACOBI_2or3_DIMENSION(ielem,djacb,kg,coren)
     DV=0.5*djacb*A1
     DEALLOCATE(CARTD,SHAPE_N,DERIV_N)
     XX=CORGP(1,kg); YY=CORGP(2,kg)
     SELECT CASE(ELMOD(etp))
	  CASE(101)
        CALL SHAPE_N_TLSL_T9(L,coren,ielem)
     END SELECT
     rvece=rvece+DV*q0*SIN(Pi*(XX+LAX)/LA)*SIN(Pi*(YY+LBY)/LB)*SHAPE_N
     DEALLOCATE(CORGP,SHAPE_N,DERIV_N)
    END IF

    IF(B3/=0.0) THEN
     DO igaus=1, 3
      kg=kg+1
      SELECT CASE(igaus)
       CASE(1)
	    L(1)=GS1; L(2)=GS1; L(3)=GS2
       CASE(2)
	    L(1)=GS1; L(2)=GS2; L(3)=GS1
	   CASE(3)
	    L(1)=GS2; L(2)=GS1; L(3)=GS1
      END SELECT
      CALL SHAPE_N_T3(L)
      CALL JACOBI_2or3_DIMENSION(ielem,djacb,kg,coren)
      DV=0.5*djacb*B3
      DEALLOCATE(CARTD,SHAPE_N,DERIV_N)
      XX=CORGP(1,kg); YY=CORGP(2,kg)
      SELECT CASE(ELMOD(etp))
	   CASE(101)
        CALL SHAPE_N_TLSL_T9(L,coren,ielem)
      END SELECT
      rvece=rvece+DV*q0*SIN(Pi*(XX+LAX)/LA)*SIN(Pi*(YY+LBY)/LB)*SHAPE_N
      DEALLOCATE(CORGP,SHAPE_N,DERIV_N)
     END DO 
    END IF

    IF(C3/=0.0) THEN
     DO igaus=1, 3
      kg=kg+1
      SELECT CASE(igaus)
       CASE(1)
        L(1)=GS3; L(2)=GS3; L(3)=GS4
	   CASE(2)
        L(1)=GS3; L(2)=GS4; L(3)=GS3
	   CASE(3)
        L(1)=GS4; L(2)=GS3; L(3)=GS3
      END SELECT
      CALL SHAPE_N_T3(L)
      CALL JACOBI_2or3_DIMENSION(ielem,djacb,kg,coren_1)
      DV=0.5*djacb*C3
      DEALLOCATE(CARTD,SHAPE_N,DERIV_N)
      XX=CORGP(1,kg); YY=CORGP(2,kg)
      SELECT CASE(ELMOD(etp))
	   CASE(101)
        CALL SHAPE_N_TLSL_T9(L,coren,ielem)
      END SELECT
      rvece=rvece+DV*q0*SIN(Pi*(XX+LAX)/LA)*SIN(Pi*(YY+LBY)/LB)*SHAPE_N
      DEALLOCATE(CORGP,SHAPE_N,DERIV_N)
     END DO 
    END IF
!!@@@
   CASE(4)
    CALL GAUSS_quadrilateral(NGAUS(IELES(ielem,MNODE+1)))
    DO igaus=1,NGAUS(IELES(ielem,MNODE+1))
     DO jgaus=1,NGAUS(IELES(ielem,MNODE+1))
      kg=kg+1
	  S=POSGP(igaus);T=POSGP(jgaus)
      CALL SHAPE_N_Q4(S,T)
      CALL JACOBI_2or3_DIMENSION(ielem,DJACB,kg,coren)
      DEALLOCATE(CARTD,SHAPE_N,DERIV_N)
      XX=CORGP(1,kg); YY=CORGP(2,kg)
      SELECT CASE(ELMOD(etp))
	   CASE(1)
        CALL SHAPE_N_TACQ(S,T,coren,ielem)
       CASE(2)
        CALL SHAPE_N_TMQ(S,T)
      END SELECT
      rvece=rvece+DJACB*WEIGP(igaus)*WEIGP(jgaus)        &
	   *q0*SIN(Pi*(XX+LAX)/LA)*SIN(Pi*(YY+LBY)/LB)*SHAPE_N
      DEALLOCATE(CORGP,SHAPE_N,DERIV_N)
     END DO
    END DO
  END SELECT
  DEALLOCATE(POSGP,WEIGP)

  DO ievab=1,MEVAB
   eload(ielem,ievab)=eload(ielem,ievab)+rvece(ievab)
  END DO
  DEALLOCATE(coren)
 END DO
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'

 DEALLOCATE(rvece)

END SUBROUTINE LOADD_2D_SINUSOIDAL_RECTANGULAR_PLATE_Q

!*============================================================================*!
!*S22. Solve_Equation_by_Front   功能: 用波前法解方程，求整体位移             *!
!*============================================================================*!

SUBROUTINE Solve_Equation_by_Front(ASDIS)
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 USE M4_Support_information
 USE M5_Para_Gauss_Integral
 USE M6_Matrix_Element
 USE M8_Shape_Function
 USE M21_INFORMATION_COMPOSITE_Layer
 REAL*8, DIMENSION(NTOTV)          :: ASDIS
 REAL*8,DIMENSION(:,:),ALLOCATABLE :: equat,ESTIF
 INTEGER,DIMENSION(:,:),ALLOCATABLE:: lnods
 REAL*8,DIMENSION(:),ALLOCATABLE   :: eqrhs,force,gload,gstif,vecrv
 INTEGER,DIMENSION(:),ALLOCATABLE  :: locel,namev,nacva,ndest,ndfro,npivo
 INTEGER::idest,idime,idofn,ielem,ielva,ievab,ifron,inode,ipoin,        &
          jdest,jevab,jfron,kelva,kevab,kexis,klast,kstar,lfron,locno,  &
		  mfron,mstif,nbufa,nfron,ngash,ngish,nikno,nlast,nposi
 REAL*8 ::cureq,pivot

 ALLOCATE(ndfro(NELEM),lnods(NELEM,MNODE))
 DO ielem=1,NELEM
  DO inode=1,MNODE
   lnods(ielem,inode)=IELES(ielem,inode)
  END DO
 END DO
 
 ndfro=0
 DO ipoin=1,NPOIN
  kstar=0
  DO ielem=1,NELEM
   DO inode=1,MNODE
    IF(lnods(ielem,inode)/=ipoin) CYCLE
	IF(kstar==0) THEN
	 kstar=ielem
	 ndfro(ielem)=ndfro(ielem)+NDOFN
	END IF
	klast=ielem
	nlast=inode
   END DO
  END DO
  IF(klast<NELEM) ndfro(klast+1)=ndfro(klast+1)-NDOFN
  lnods(klast,nlast)=-ipoin
 END DO   

 nfron=0;mfron=0
 DO ielem=1,NELEM
  nfron=nfron+ndfro(ielem)
  IF(nfron>mfron) mfron=nfron
 END DO
 WRITE(3,*)
 WRITE(3,'(A,I8)') ' *<7>.Maximum frontwidth encountered =',mfron
 mstif=(mfron*mfron-mfron)/2.0+mfron
 DEALLOCATE(ndfro)
 ALLOCATE(eqrhs(NTOTV),equat(mfron,NTOTV),&
  vecrv(mfron),gload(mfron),gstif(mstif),locel(MEVAB),&
  npivo(NTOTV),nacva(mfron),namev(NTOTV),&
  ndest(MEVAB))

 eqrhs=0.0; gstif=0.0; gload=0.0
 vecrv=0.0; nacva=0.0
 equat=0.0

 nbufa=0
 REWIND 7
 REWIND 9

 nfron=0
 kelva=0
 ALLOCATE(ESTIF(MEVAB,MEVAB),force(MEVAB))
 DO ielem=1,NELEM
  READ(9,*) (force(ievab),ievab=1,MEVAB)
  kevab=0
  CALL STIFFNESS_MATRIX_ELE(ielem,ESTIF)
  locel=0  !!! new added
  DO inode=1,MNODE
   locno=lnods(ielem,inode)
   DO idofn=1,NDOFN
    nposi=(inode-1)*NDOFN+idofn
	IF(locno>0) locel(nposi)=(locno-1)*NDOFN+idofn
	IF(locno<0) locel(nposi)=(locno+1)*NDOFN-idofn
   END DO
  END DO
  DO ievab=1,MEVAB
   nikno=IABS(locel(ievab))
   IF(nikno==0) CYCLE  !! NEW ADDED
   kexis=0
   DO ifron=1,nfron
    IF(nikno==nacva(ifron)) THEN
	 kevab=kevab+1
	 kexis=1
	 ndest(kevab)=ifron
	END IF
   END DO
   IF(kexis/=0) CYCLE
   DO ifron=1,mfron
    IF(nacva(ifron)/=0) CYCLE
    nacva(ifron)=nikno
    kevab=kevab+1
    ndest(kevab)=ifron
    EXIT
   END DO
   IF(ndest(kevab)>nfron) nfron=ndest(kevab)
  END DO

  DO ievab=1,MEVAB
   idest=ndest(ievab)
   IF(idest==0) CYCLE !!!
   gload(idest)=gload(idest)+force(ievab)
   DO jevab=1,ievab
    jdest=ndest(jevab)
	ngash=NFUNC(idest,jdest)
	ngish=NFUNC(jdest,idest)
	IF(jdest>=idest) gstif(ngash)=gstif(ngash)+ESTIF(ievab,jevab)
	IF(jdest<idest) gstif(ngish)=gstif(ngish)+ESTIF(ievab,jevab)
   END DO
  END DO

  DO ievab=1,MEVAB
   nikno=-locel(ievab)
   IF(nikno<=0) CYCLE
   DO ifron=1,nfron
    IF(nacva(ifron)/=nikno) CYCLE
	nbufa=nbufa+1
	IF(nbufa>NTOTV) THEN
	 nbufa=1