atf1.f90

薄厚通用的板单元

   POSGP(2)=0.7974269853
   POSGP(3)=0.4701420641
   POSGP(4)=0.0597158717
   WEIGP(1)=0.2250000000
   WEIGP(2)=0.1259391805
   WEIGP(3)=0.1323941527   
 END SELECT
   
END SUBROUTINE GAUSS_triangular

!*============================================================================*!
!*S10. Generate_Nodal_Load   功能: 生成单元等效结点荷载向量                    *!
!*============================================================================*!

SUBROUTINE Generate_Nodal_Load
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 USE M8_Shape_Function
 REAL*8, DIMENSION(NELEM,MEVAB) :: eload
 INTEGER:: igrav,iplod,isufc,idsin
 INTEGER:: ielem,jevab

 WRITE(3,*) '*<6>.LOAD INFORMATION'
 eload=0.0
 
 SELECT CASE(TYPE2)
  CASE(1:2)
   READ(1,*) iplod, isufc, igrav
   WRITE(3, *) '*iplod=',iplod,'      *isufc=',isufc,'      *igrav=',igrav
   IF(iplod/=0) CALL LOADC(eload)
   IF(isufc/=0) CALL LOADD_2D_Plane(eload)
   IF(igrav/=0) CALL LOADG(eload)
  CASE(3)
   READ(1,*) iplod, isufc
   WRITE(3, *) '*iplod=',iplod,'      *isufc=',isufc 
   IF(iplod/=0) CALL LOADC(eload)
   IF(isufc/=0) CALL LOADD_2D_Plate(eload)
  CASE(4)
   READ(1,*) iplod, isufc, idsin
   WRITE(3, *) '*iplod=',iplod,'      *isufc=',isufc,'      *idsin=',idsin 
   IF(iplod/=0) CALL LOADC(eload)
   IF(isufc/=0) CALL LOADD_2D_Composite_Plate(eload)
   IF(idsin/=0) CALL LOADD_2D_SINUSOIDAL_RECTANGULAR_PLATE_Q(eload)
 END SELECT

 DO ielem=1,NELEM
  WRITE(9,*) (eload(ielem,jevab),jevab=1,MEVAB)
 END DO

END SUBROUTINE Generate_Nodal_Load

!*============================================================================*!
!*S11. LOADC   功能: 二维或三维单元结点集中荷载                                *!
!*============================================================================*!

SUBROUTINE LOADC(eload)
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 REAL*8, DIMENSION(NELEM,MEVAB) :: eload
 INTEGER:: idofn,ielem,inode,lodpt,ngash,nplod,nloca,k
 REAL*8, DIMENSION(NDOFN) :: point

 READ(1,*) nplod
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'
 WRITE(3,*) ' *There are ',nplod,'   nodes under concentrated forces.'
 IF(nplod>0) THEN
  k=0
  DO
   k=k+1
   IF(k>nplod) EXIT
   READ(1,*) lodpt, (point(idofn), idofn=1,NDOFN)
   WRITE(3,'(I9,6F10.3)') lodpt, (point(idofn), idofn=1,NDOFN)
   DO ielem=1, NELEM
    DO inode=1, MNODE
	 nloca=IABS(IELES(ielem, inode))
	 IF(lodpt==nloca) EXIT
	END DO
    IF(lodpt==nloca) EXIT
   END DO
   DO idofn=1, NDOFN
	ngash=(inode-1)*NDOFN+idofn
	eload(ielem,ngash)=eload(ielem,ngash)+point(idofn)
   END DO
  END DO
 END IF
 WRITE(3,'(A)') ' -----------------------------------------------------------&
-----------'

END SUBROUTINE LOADC

!*============================================================================*!
!*S12. LOADD_2D_Plane   功能: 二维平面膜元边界分布荷载                        *!
!*============================================================================*!

SUBROUTINE LOADD_2D_Plane(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
 INTEGER::nel,nic,nconp,nedge,nodps ! nodps:单元每边的结点数

 coren=0.0
 READ(1,*) nedge
 WRITE(3,*) ' *There are ',nedge,'   sides under 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,61,81:82)
	 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,61:62)
	 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)') ' -----------------------------------------------------------&
-----------'

END SUBROUTINE LOADD_2D_Plane

!*============================================================================*!
!*S13. LOADD_2D_PLANE_METHOD_1   功能: 二维平面等参单元分布荷载               *!
!*============================================================================*!

SUBROUTINE LOADD_2D_PLANE_METHOD_1(nel,nodps,cores,noprs,press,eload)
 USE M0_PROBLEM_TYPE
 USE M1_Gobal_information
 USE M2_Local_information
 USE M3_Stru_information
 USE M5_Para_Gauss_Integral
 USE M8_Shape_Function
 REAL*8, DIMENSION(NDOFN)         :: dgash,pgash
 REAL*8, DIMENSION(NELEM,MEVAB)   :: eload
 INTEGER, DIMENSION(nodps)        :: noprs
 REAL*8, DIMENSION(nodps,2)       :: press
 REAL*8, DIMENSION(NDIME,nodps)   :: cores
 REAL*8, DIMENSION(MEVAB)         :: rvece
 REAL*8, DIMENSION(3)             :: L
 REAL*8, DIMENSION(:), ALLOCATABLE:: SHAPE_N_SIDE,DERIV_N_SIDE
 REAL*8 :: S, DV, Px, Py
 INTEGER:: nel,nodps,ngaus_s 
 INTEGER:: idofn,igaus,inode,inods 

 rvece=0.0

 SELECT CASE(nodps)
  CASE(2)
   ALLOCATE(SHAPE_N_SIDE(2),DERIV_N_SIDE(2))
   ngaus_s=2
  CASE(3)
   ALLOCATE(SHAPE_N_SIDE(3),DERIV_N_SIDE(3))
   ngaus_s=3
 END SELECT

 CALL GAUSS_quadrilateral(ngaus_s)

 DO igaus=1, ngaus_s
  S=POSGP(igaus)
  SELECT CASE(nodps)
   CASE(2)
    SHAPE_N_SIDE(1)=(1.0-S)/2.0; SHAPE_N_SIDE(2)=(1.0+S)/2.0
    DERIV_N_SIDE(1)=-0.5       ; DERIV_N_SIDE(2)=0.5
   CASE(3)
    SHAPE_N_SIDE(1)=S*(S-1.0)/2.0; SHAPE_N_SIDE(2)=1.0-S*S
    SHAPE_N_SIDE(3)=S*(S+1.0)/2.0; DERIV_N_SIDE(1)=(2*S-1.0)/2.0
    DERIV_N_SIDE(2)=-2.0*S       ; DERIV_N_SIDE(3)=(2*S+1.0)/2.0
  END SELECT
  DO idofn=1,NDOFN     ! Cyc 2s
   pgash(idofn)=0.0; dgash(idofn)=0.0
   DO inods=1,nodps   
    pgash(idofn)=pgash(idofn)+press(inods,idofn)*SHAPE_N_SIDE(inods)
    dgash(idofn)=dgash(idofn)+cores(idofn,inods)*DERIV_N_SIDE(inods)
   END DO      
  END DO      
  DV=WEIGP(igaus)
  Px=dgash(1)*pgash(2)-dgash(2)*pgash(1)
  Py=dgash(1)*pgash(1)+dgash(2)*pgash(2)
  DO inods=1,nodps   
     rvece(NDOFN*(inods-1)+1)=rvece(NDOFN*(inods-1)+1)+Px*SHAPE_N_SIDE(inods)*DV
     rvece(NDOFN*(inods-1)+2)=rvece(NDOFN*(inods-1)+2)+Py*SHAPE_N_SIDE(inods)*DV
  END DO         
 END DO       

 DEALLOCATE(SHAPE_N_SIDE, DERIV_N_SIDE)
 DEALLOCATE(POSGP,WEIGP)
   
 DO inods=1,nodps
  DO inode=1,MNODE
   IF(IELES(nel,inode)==noprs(inods)) THEN
    eload(nel,NDOFN*(inode-1)+1)=eload(nel,NDOFN*(inode-1)+1)+rvece(NDOFN*(inods-1)+1)
    eload(nel,NDOFN*(inode-1)+2)=eload(nel,NDOFN*(inode-1)+2)+rvece(NDOFN*(inods-1)+2)
   END IF
  END DO
 END DO

END SUBROUTINE LOADD_2D_PLANE_METHOD_1

!*============================================================================*!
!*S14. LOADD_2D_PLANE_METHOD_2   功能: 二维平面面积坐标四边形单元分布荷载     *!
!*============================================================================*!

SUBROUTINE LOADD_2D_PLANE_METHOD_2(nel,nconp,nodps,coren,noprs,press,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(nodps)         :: Rx, Ry
 REAL*8, DIMENSION(NELEM,MEVAB)   :: eload
 REAL*8, DIMENSION(NDIME,MNODE)   :: coren
 REAL*8, DIMENSION(nodps,2)       :: press
 REAL*8, DIMENSION(MEVAB)         :: rvece
 REAL*8 :: S,T,DV,dl,sinsita,cossita,a1,a2,a3,b1,b2,b3,RXF,RXY
 INTEGER, DIMENSION(nodps)        :: noprs
 INTEGER:: nel,nconp,nodps
 INTEGER:: ievab,igaus,inode,inods

 rvece=0.0
 dl=DSQRT((COORD(noprs(nodps),1)-COORD(noprs(1),1))**2  &
           +(COORD(noprs(nodps),2)-COORD(noprs(1),2))**2)
 sinsita=(COORD(noprs(nodps),2)-COORD(noprs(1),2))/dl
 cossita=(COORD(noprs(nodps),1)-COORD(noprs(1),1))/dl
 DO inods=1,nodps   
  Rx(inods)=-press(inods,1)*sinsita+press(inods,2)*cossita
  Ry(inods)= press(inods,1)*cossita+press(inods,2)*sinsita
 END DO          
 a1=0.5*Rx(3)+0.5*Rx(1)-Rx(2); a2 =0.5*Rx(3)-0.5*Rx(1); a3=Rx(2)
 b1=0.5*Ry(3)+0.5*Ry(1)-Ry(2); b2 =0.5*Ry(3)-0.5*Ry(1); b3=Ry(2)
 
 CALL GAUSS_quadrilateral(NGAUS(IELES(nel,MNODE+1)))
 DO igaus=1,NGAUS(IELES(nel,MNODE+1))
  SELECT CASE(nconp)
   CASE(1)
    S= POSGP(igaus); T=-1.0
   CASE(2)
    S= 1.0; T= POSGP(igaus)
   CASE(3)
    S= POSGP(igaus); T= 1.0
   CASE(4)
    S=-1.0; T= POSGP(igaus)
  END SELECT
  CALL Shape_Function_Q(nel,S,T,coren)
  DV=WEIGP(igaus)*dl*0.5
  DO inode=1,NNODE(IELES(nel,MNODE+1))
   SELECT CASE(nconp)
	CASE(1)  
	 rvece(2*inode-1)=rvece(2*inode-1)+RXF(a1,a2,a3,S)*SHAPE_N(inode)*DV
     rvece(2*inode)=rvece(2*inode)+RYF(b1,b2,b3,S)*SHAPE_N(inode)*DV
    CASE(2)
     rvece(2*inode-1)=rvece(2*inode-1)+RXF(a1,a2,a3,T)*SHAPE_N(inode)*DV
     rvece(2*inode)=rvece(2*inode)+RYF(b1,b2,b3,T)*SHAPE_N(inode)*DV
    CASE(3)  
     rvece(2*inode-1)=rvece(2*inode-1)-RXF(a1,a2,a3,S)*SHAPE_N(inode)*DV
     rvece(2*inode)=rvece(2*inode)-RYF(b1,b2,b3,S)*SHAPE_N(inode)*DV
    CASE(4)
	 rvece(2*inode-1)=rvece(2*inode-1)-RXF(a1,a2,a3,T)*SHAPE_N(inode)*DV
     rvece(2*inode)=rvece(2*inode)-RYF(b1,b2,b3,T)*SHAPE_N(inode)*DV
   END SELECT
  END DO  
  SELECT CASE(ELMOD(IELES(nel,MNODE+1)))
   CASE(2)
    DEALLOCATE(SHAPE_N,DERIV_N) 
   CASE(3:4)
    DEALLOCATE(SHAPE_N,DERIV_N,CARTD) 
  END SELECT
 END DO

 DO ievab=1,MEVAB 
  eload(nel,ievab)=eload(nel,ievab)+rvece(ievab)
 END DO    

 DEALLOCATE(POSGP,WEIGP)

END SUBROUTINE LOADD_2D_PLANE_METHOD_2
!--------------------------------------------------------------------
FUNCTION RXF(a1,a2,a3,S) RESULT(RXF_RESULT)
   IMPLICIT NONE
   REAL*8::a1,a2,a3,S, RXF_RESULT
   RXF_RESULT=a1*S*S+a2*S+a3
END FUNCTION RXF
FUNCTION RYF(b1,b2,b3,S) RESULT(RYF_RESULT)
   IMPLICIT NONE
   REAL*8::b1,b2,b3,S, RYF_RESULT
   RYF_RESULT=b1*S*S+b2*S+b3
END FUNCTION RYF

!*============================================================================*!
!*S15. LOADD_2D_PLANE_METHOD_3   功能: 二维平面带旋转自由度三角形单元分布荷载 *!
!*============================================================================*!

SUBROUTINE LOADD_2D_PLANE_METHOD_3(nel,nconp,nodps,coren,noprs,press,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(nodps,2)     :: press
 REAL*8, DIMENSION(nodps)       :: Rx,Ry 
 INTEGER, DIMENSION(nodps)      :: noprs 
 REAL*8, DIMENSION(MEVAB)       :: rvece
 REAL*8, DIMENSION(4)           :: P,W 
 INTEGER:: nodps,nconp,nel
 REAL*8 :: dl,sinsita,cossita,DV
 INTEGER:: inods,inode,igaus,ievab
 REAL*8,DIMENSION(3) :: L
      
 dl=DSQRT((COORD(noprs(NODPS),1)-COORD(noprs(1),1))**2  &
         +(COORD(noprs(NODPS),2)-COORD(noprs(1),2))**2)
 sinsita=(COORD(noprs(NODPS),2)-COORD(noprs(1),2))/dl
 cossita=(COORD(noprs(NODPS),1)-COORD(noprs(1),1))/dl
 DO inods=1,nodps 
  Rx(inods)=-press(inods,1)*sinsita+press(inods,2)*cossita
  Ry(inods)= press(inods,1)*cossita+press(inods,2)*sinsita
 END DO

 rvece=0.0
 CALL GAUSS_quadrilateral(4)
 P=POSGP; W=WEIGP
 DEALLOCATE(POSGP,WEIGP)

 DO igaus=1,4
  SELECT CASE(nconp)
   CASE(1)
    L(1)=(P(igaus)+1.0)/2.0; L(2)=1.0-L(1); L(3)=0.0
   CASE(2)
    L(1)=0.0; L(2)=(P(igaus)+1.0)/2.0; L(3)=1.0-L(2)
   CASE(3)
    L(1)=1.0-L(3); L(2)=0.0; L(3)=(P(igaus)+1.0)/2.0
  END SELECT