atf1.f90

薄厚通用的板单元

!*============================================================================*!
MODULE M0_PROBLEM_TYPE         ! 问题类型参数
   IMPLICIT NONE
   INTEGER:: TYPE1,TYPE2
END MODULE M0_PROBLEM_TYPE

MODULE M1_Gobal_information    ! 整体控制信息参数
   IMPLICIT NONE
   INTEGER:: NPOIN,NELEM,NVFIX,NMATS,NELTP,NOUTP
   INTEGER:: NTOTV
END MODULE M1_Gobal_information

MODULE M2_Local_information    ! 单元局部信息参数
   IMPLICIT NONE
   INTEGER,DIMENSION(:), ALLOCATABLE :: ELMOD,ELESP
   INTEGER,DIMENSION(:), ALLOCATABLE :: NNODE,NGAUS,NEVAB
   INTEGER:: NDIME, NPROP
   INTEGER:: NDOFN, NSTRE
   INTEGER:: MNODE, MEVAB
END MODULE M2_Local_information

MODULE M3_Stru_information     ! 结构信息
   IMPLICIT NONE
   INTEGER,DIMENSION(:,:),ALLOCATABLE :: IELES
   REAL*8, DIMENSION(:,:),ALLOCATABLE :: COORD
   REAL*8, DIMENSION(:,:),ALLOCATABLE :: PROPS
END MODULE M3_Stru_information

MODULE M4_Support_information  ! 约束信息
  IMPLICIT NONE
  REAL*8,  DIMENSION(:),ALLOCATABLE :: FIXED ! FIXED(NTOTV)
  INTEGER, DIMENSION(:),ALLOCATABLE :: IFFIX ! IFFIX(NTOTV)
  INTEGER, DIMENSION(:),ALLOCATABLE :: NOFIX ! NOFIX(NVFIX)
END MODULE M4_Support_information

MODULE M5_Para_Gauss_Integral  ! 数值积分常数 
   IMPLICIT NONE
   REAL*8,DIMENSION(:),ALLOCATABLE :: POSGP,WEIGP
END MODULE M5_Para_Gauss_Integral

MODULE M6_Matrix_Element       ! 单元计算矩阵 
   IMPLICIT NONE
   REAL*8,DIMENSION(:,:),ALLOCATABLE :: BMATX,BSMATX,DMATX,SMATX,QMATX,CARTD
END MODULE M6_Matrix_Element

MODULE M8_Shape_Function       ! 单元形函数及其导数以及单元高斯积分点坐标
   IMPLICIT NONE
   REAL*8,DIMENSION(:),  ALLOCATABLE :: SHAPE_N
   REAL*8,DIMENSION(:,:),ALLOCATABLE :: DERIV_N,CORGP
   REAL*8,DIMENSION(:,:),ALLOCATABLE :: SHAPE_N_matrix
END MODULE M8_Shape_Function

MODULE M21_INFORMATION_COMPOSITE_Layer  ! 层合系统的材料信息
  IMPLICIT NONE
   INTEGER:: Num_Layer
   REAL*8 :: K1,K2
   REAL*8,DIMENSION(:,:),ALLOCATABLE :: Inform_Layer
   INTEGER,DIMENSION(:), ALLOCATABLE :: Layer_type
   REAL*8,DIMENSION(:,:,:),ALLOCATABLE :: Q_MATRIX_PZ,Q_MATRIX_ZZ
   REAL*8,DIMENSION(:,:,:),ALLOCATABLE :: C_MATRIX_PZ
   REAL*8,DIMENSION(6,6) :: Cp
   REAL*8,DIMENSION(2,2) :: Cs
END MODULE M21_INFORMATION_COMPOSITE_Layer
!*============================================================================*!
!*===================================GCFEM====================================*!
!*2.主程序                                                                    *!
!*===================================GCFEM====================================*!
!*============================================================================*!
PROGRAM GCFEM2000_V1	  
 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(:), ALLOCATABLE :: ASDIS
 CHARACTER(LEN=65):: OUTPUT_FILENAME         ! 结果输出文件名   
 REAL time_begin, time_end                   ! 定义程序始终时间
 CHARACTER(LEN=12):: Real_clock(3)
 INTEGER Date_time(8) 
 					  
 WRITE(*,'(A)') '  ===========================================================' 
 WRITE(*,'(A)') '  *        Welcome to use the FEM-program GCFEM2000         *'
 WRITE(*,'(A)') '  *    This program can solve such problems at present:     *'
 WRITE(*,'(A)') '  *       <1>  SOLID-LINEAR-ELASTIC-STATICS                 *'
 WRITE(*,'(A)') '  *           1. Plane stress problems                      *'            
 WRITE(*,'(A)') '  *           2. Plane strain problems                      *'            
 WRITE(*,'(A)') '  *           3. Plate bending problems                     *'
 WRITE(*,'(A)') '  *           4. Composite plate bending problems           *'
 WRITE(*,'(A)') '  *=========================================================*'
 WRITE(*,*) 

 CALL date_and_time(Real_clock(1),Real_clock(2),Real_clock(3),Date_time)
 CALL OPENFILE(OUTPUT_FILENAME,Date_time)              ! 打开备用文件
 CALL CPU_TIME ( time_begin )

! 读取问题类型, 调用相应求解程序 

 SELECT CASE(TYPE1)
   CASE(1)
     write(*,'(A)') ' Reading and writting input datas  ......'
    CALL Input_Data
     write(*,'(A)') ' Generating efficient nodal Load   ......'
    CALL Generate_Nodal_Load
	 write(*,'(A)') ' done!'
	 write(*,'(A)') ' Solving Equation by Frontal method......'
     write(3,*)
    ALLOCATE(ASDIS(NTOTV))
	ASDIS=0.0
	CALL Solve_Equation_by_Front(ASDIS)
     write(*,'(A)') ' done!'
     write(*,'(A)') ' Outputing results ......'
    CALL Output(ASDIS)
     write(*,'(A)') ' done!'
     write(3,'(A)')'*===========================End of data_out============&
	 =================*'
   CASE DEFAULT
    WRITE(*,'(A)') 'Error: the Program can not run because no such TYPE1 exists'
    WRITE(3,'(A)') 'Error: the Program can not run because no such TYPE1 exists'
    WRITE(5,'(A)') 'Error: the Program can not run because no such TYPE1 exists'
	STOP
 END SELECT

 CALL CPU_TIME ( time_end )
 PRINT*,  'Time of operation was ', time_end-time_begin , ' seconds'
 WRITE(5,*)
 WRITE(5,*)'Time of operation was ', time_end-time_begin , ' seconds'
 WRITE(5,'(A)') '*========================End of the outputfile=========&
 ===============*'
 WRITE(*,'(A)') '============Running  Completed============='
 WRITE(*,'(A,A)') '===Please open and see output file:',OUTPUT_FILENAME

!-------------------------------------------------------------------!
   CLOSE(1)
   CLOSE(3)
   CLOSE(5)
   CLOSE(11)

END PROGRAM GCFEM2000_V1
!*============================================================================*!
!*===================================GCFEM====================================*!
!*3.子程序清单                                                                *!
!*===================================GCFEM====================================*!
!*============================================================================*!
!*S1.     功能: 建立输入输出文件                                     *!
!*============================================================================*!

SUBROUTINE OPENFILE(OUTPUT_FILENAME,Date_time)
 USE M0_PROBLEM_TYPE
 CHARACTER(LEN=65):: INPUT_FILENAME,OUTPUT_FILENAME   
 LOGICAL IF_FILE_EXITS
 CHARACTER(LEN=150):: Title    
 INTEGER Date_time(8) 

   WRITE(*,'(A,\)')'    Please type INPUT-FILE name:?' 
   READ(*,'(A)')INPUT_FILENAME  					   
   INQUIRE(FILE=INPUT_FILENAME,EXIST=IF_FILE_EXITS)	   
     IF(.NOT.IF_FILE_EXITS) THEN					   
     PRINT*,'   Error: This file does not exits.'	   
	 STOP											   
   ENDIF											   
   OPEN(1,FILE=INPUT_FILENAME,STATUS='OLD')            
   WRITE(*,'(A,\)')'    Please type OUTPUT-FILE name:?'
   READ(*,'(A)') OUTPUT_FILENAME					   
   OPEN(5,FILE=OUTPUT_FILENAME,STATUS='UNKNOWN')	   

   WRITE(5,'(A)') '*==========================================================&
===========*'
   WRITE(5,'(A,A)') '*               This file is the result of GCFEMP2000   '&
    ,'              *'
   WRITE(5,'(A)') '*==========================================================&
===========*'
   WRITE(5,*)
   WRITE(5,'(A,I2,A,I2,A,I4)') '*Today is ',Date_time(3),'-',Date_time(2), &
    '-',Date_time(1)
   WRITE(5,'(3(A,I2),A,I3)') '*The beginning time is ',Date_time(5),':',Date_time(6), &
    ':',Date_time(7),'.',Date_time(8)
   WRITE(5,*)

   OPEN(3,FILE='DATA_OUT',STATUS='UNKNOWN')        !建立数据返回文件
   WRITE(3,'(A)') '*==========================================================&
===========*'
   WRITE(3,'(A,A)') '*             This file is the inputdata of GCFEMP2000  '&
     ,'              *'
   WRITE(3,'(A)') '*==========================================================&
===========*'
   WRITE(3,*)
   WRITE(3,'(A,I2,A,I2,A,I4)') '*Today is ',Date_time(3),'-',Date_time(2), &
    '-',Date_time(1)
   WRITE(3,'(3(A,I2),A,I3)') '*The beginning time is ',Date_time(5),':',Date_time(6), &
    ':',Date_time(7),'.',Date_time(8)
   WRITE(3,*)

   READ (1,*) Title
   WRITE(3,'(A,A)') '*',Title
   WRITE(5,'(A,A)') '*',Title

   READ(1,*) TYPE1, TYPE2
   
   SELECT CASE(TYPE1)
    CASE(1)
     write(5,*)
     WRITE(5,'(A)') 'FOR SOLID LINEAR ELASTIC STATIC'
     SELECT CASE(TYPE2)
	  CASE(1)
	   WRITE(5,'(A)')'*   1-1. Plane stress problems                *'
	  CASE(2)
	   WRITE(5,'(A)')'*   1-2. Plane strain problems                *'
      CASE(3)
       WRITE(5,'(A)')'*   1-3. Plate bending problems               *'
      CASE(4)
       WRITE(5,'(A)')'*   1-4. Composite Plate bending problems     *'
      CASE DEFAULT
    WRITE(*,'(A)') 'Error: the Program can not run because no such TYPE2 exists'
    WRITE(3,'(A)') 'Error: the Program can not run because no such TYPE2 exists'
    WRITE(5,'(A)') 'Error: the Program can not run because no such TYPE2 exists'
	  STOP
     END SELECT
    CASE DEFAULT
    WRITE(*,'(A)') 'Error: the Program can not run because no such TYPE1 exists'
    WRITE(3,'(A)') 'Error: the Program can not run because no such TYPE1 exists'
    WRITE(5,'(A)') 'Error: the Program can not run because no such TYPE1 exists'
	 STOP
   END SELECT

   WRITE(5,*)
   WRITE(5,'(A)') '! You are solving equations by FRONTAL method'
   OPEN(7,FILE='FILE_7.DAT',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
   OPEN(9,FILE='FILE_9.DAT',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')

   WRITE(5,*)
   OPEN(11,FILE='SMATX.DAT',STATUS='UNKNOWN')

END SUBROUTINE OPENFILE

!*============================================================================*!
!*S2. Input_Data    功能: 读取并返回输入数据, 建立单元常数                    *!
!*============================================================================*!

SUBROUTINE Input_Data
 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 M21_INFORMATION_COMPOSITE_Layer
 INTEGER::idofn,idist,ielem,ielet,ifdof,ifpre,imats,ilayer,inode,ivfix
 INTEGER::jdime,jdime1,jdofn,jdofn1,jnmat,jnode,jnode1,ngash,nloca
 INTEGER::k,l,nel,nod,nty,etp,n_node,n_gaus,n_evab,n_elsp,n_dofn
 REAL*8,DIMENSION(:,:),ALLOCATABLE :: PRESC
 
!(1).读取返回控制信息

 READ(1,*) NPOIN, NELEM, NVFIX, NMATS, NELTP, NOUTP
 
 ALLOCATE(NNODE(NELTP),NGAUS(NELTP),NEVAB(NELTP),ELMOD(NELTP),ELESP(NELTP))
 SELECT CASE(TYPE2)
  CASE(1:3)
   NDIME=2; NPROP=4; NSTRE=3
  CASE(4)
   NDIME=2; NPROP=7; NSTRE=6
 END SELECT
 READ(1,*) (ELMOD(ielet),ielet=1,NELTP)
 NNODE=0; NGAUS=0; NEVAB=0; ELESP=0; MNODE=0; MEVAB=0; NDOFN=0
 DO ielet=1, NELTP
  etp=ELMOD(ielet)
  SELECT CASE(TYPE2)         	   
   CASE(1:2)
    CALL GCVAR_A_1(etp,n_node,n_gaus,n_evab,n_elsp,n_dofn)  
   CASE(3)	                       
    CALL GCVAR_A_2(etp,n_node,n_gaus,n_evab,n_elsp,n_dofn)
   CASE(4)	                       
    CALL GCVAR_A_3(etp,n_node,n_gaus,n_evab,n_elsp,n_dofn)
  END SELECT
  NNODE(ielet)=n_node; NGAUS(ielet)=n_gaus
  NEVAB(ielet)=n_evab; ELESP(ielet)=n_elsp
  IF(n_node>MNODE) MNODE=n_node
  IF(n_evab>MEVAB) MEVAB=n_evab
  IF(n_dofn>NDOFN) NDOFN=n_dofn
 END DO

 NTOTV=NPOIN*NDOFN  

 WRITE(3,*)
 WRITE(3,*) '*<1>.CONTRAL INFORMATION'
 WRITE(3,FMT='(A,I9)') '   Max   NUMber   of   NODEs   in   STRUture    =',NPOIN
 WRITE(3,FMT='(A,I9)') '   Max   NUMber   of   ELEMENTs   in  STRUture  =',NELEM
 WRITE(3,FMT='(A,I9)') '   Max   NUMber   of   SUPPOPT   POINTs         =',NVFIX
 WRITE(3,FMT='(A,I9)') '   Max           NUMber         of      MATERAL =',NMATS
 WRITE(3,FMT='(A,I9)') '   The Contral Parameter for outputing  results =',NOUTP
 WRITE(3,*) '----------------------------------------------------------------------'
  DO ielet=1,NELTP
  WRITE(3,FMT='(A,I9)') '   Element type sequence No.',ielet
  WRITE(3,FMT='(A,I9)') '   TYPE  NUMber     of      element       MODEL =', &
   ELMOD(ielet)
  END DO
 WRITE(3,*) '----------------------------------------------------------------------'
  DO ielet=1,NELTP
  WRITE(3,FMT='(A,I9)') '   Element type sequence No.',ielet
  WRITE(3,FMT='(A,I9)') '   NUMber    of    NODEs    in   each   ELEment =', &
   NNODE(ielet)
 END DO
 WRITE(3,*) '----------------------------------------------------------------------'
 WRITE(3,FMT='(A,I9)') '   Max   NUMber   of   NODEs   per   element    =',MNODE
 WRITE(3,FMT='(A,I9)') '   NUMber    of    DIMENSION     per    element =',NDIME
 WRITE(3,FMT='(A,I9)') '   NUMber    of    MATERAL       per    ELEment =',NPROP
 WRITE(3,FMT='(A,I9)') '   NUMber  of  DEGREE  of   freedom   per  NODE =',NDOFN
 WRITE(3,FMT='(A,I9)') '   NUMber    of    STRESS   in   the     POINTs =',NSTRE
 WRITE(3,*) '----------------------------------------------------------------------'
 DO ielet=1,NELTP
  WRITE(3,FMT='(A,I9)') '   Element type sequence No.',ielet
  WRITE(3,FMT='(A,I9)') '   NUMber    of    GAUSS  POINTs  per   element =', &
   NGAUS(ielet)
 END DO
 WRITE(3,*) '----------------------------------------------------------------------'
 WRITE(3,FMT='(A,I9)') '   Max   NUMber  of  VARiable  in the STRUcture =',NTOTV
 WRITE(3,*) '----------------------------------------------------------------------'
  DO ielet=1,NELTP
  WRITE(3,FMT='(A,I9)') '   Element type sequence No.',ielet
  WRITE(3,FMT='(A,I9)') '   NUMber     of     VARiable    per    ELEment =', &
   NEVAB(ielet)
  END DO
 WRITE(3,*) '----------------------------------------------------------------------'
 WRITE(3,FMT='(A,I9)') '   Max   NUMber  of  VARiable    per    Element =',MEVAB
 WRITE(3,*)

!(2).分配数组空间,将所有数组填零
 ALLOCATE(IELES(NELEM,MNODE+2), COORD(NPOIN,NDIME))
 IELES=0; COORD=0.0

!(3).读写单元信息
 WRITE(3,*) '*<2>.ELEMENT INFORMATION'
  DO ielem=1,NELEM 
    READ(1,*) nel, (IELES(ielem,jnode1),jnode1=1,MNODE+2)
    WRITE(3,FMT='(2X,I5,10I6)') nel, (IELES(ielem,jnode1),jnode1=1,MNODE+2)
  END DO
 WRITE(3,*)

!(4).读写结点坐标,生成边中结点 
 WRITE(3,*) '*<3>.NODES COODINATE'
 DO inode=1,NPOIN 
   READ(1,*) nod,(COORD(inode,jdime),jdime=1,NDIME)
   IF(nod/=inode) THEN
	DO jdime=1,NDIME 
  	  COORD(inode,jdime)=0.0
 	END DO 
   BACKSPACE 1
   END IF
 END DO           
 IF(TYPE2<3) CALL CORMN	  ! 如有必要计算边中结点坐标
 DO inode=1,NPOIN 
   WRITE(3,FMT='(6X,I5,3F10.3)') inode, (COORD(inode,jdime),jdime=1,NDIME)
 ENDDO
 WRITE(3,*)

!(5).读写结构约束信息
 WRITE(3,*) '*<4>.SUPPORT INFORMATION'
 ALLOCATE(IFFIX(NTOTV), NOFIX(NVFIX), PRESC(NVFIX,NDOFN))
 IFFIX=0

 DO ivfix=1,NVFIX 
  READ(1,*) NOFIX(ivfix),ifpre,(PRESC(ivfix,jdofn),jdofn=1,NDOFN)
  WRITE(3,'(1X,I4,5X,I8,5X,5F10.6)') NOFIX(ivfix),ifpre, &
    (PRESC(ivfix,jdofn),jdofn=1,NDOFN)
  nloca=(NOFIX(ivfix)-1)*NDOFN
  ifdof=10**(NDOFN-1)
 
  DO idofn=1,NDOFN 
   ngash=nloca+idofn
   IF(ifpre>=ifdof) THEN
    IFFIX(ngash)=1
    ifpre=ifpre-ifdof
   END IF
   ifdof=ifdof/10
  END DO
 
 END DO