flagshyp.f90
来自「大应变超弹性平面分析程序」· F90 代码 · 共 1,766 行 · 第 1/5 页
F90
1,766 行
READ (5, '(a)') name
IF ( .not. rest ) then
OPEN (2, file = name, status = 'unknown', form = 'formatted', &
err = 20)
ELSE
OPEN (2, file = name, status = 'old', form = 'formatted', &
err = 20, action = 'write', position = 'append')
!b access = 'append')
END IF
!
! Opens the restart file
!
30 WRITE (6, 104)
104 FORMAT(' Enter the restart file name: ')
READ (5, '(a)') name
OPEN (3, file = name, status = 'unknown', form = 'unformatted', &
err = 30)
!
END SUBROUTINE welcome
!-----------------------------------------------------------------------
SUBROUTINE elinfo (mnode, mgaus, ndime, nnode, ngaus, nnodb, &
ngaub, eledb, elbdb, eltyp)
!-----------------------------------------------------------------------
!
! Obtains element information
!
! mnode --> maximum number of nodes per element
! mgaus --> maximum number of gauss points per element
! nnode --> number of nodes per element
! ngaus --> number of gauss points per element
! ndime --> number of dimensions for the given element
! nnodb --> number of nodes per boundary element
! ngaub --> number of gauss points per boundary element
! eledb --> element data matrix of dimensions
! (ndime+1,nnode+1,ngaus):
!
! | N1, N2,.., Nnnode, wgt |
! | |
! | dN1 dN2 dNnnode |
! | ---,---,...,-------, xi |
! | dxi dxi dxi | X number of Gauss points
! | |
! | dN1 dN2 dNnnode |
! | ---,---,...,-------, eta |
! | det det det |
!
! elbdb --> boundary elements data, format as above
! eltyp --> element type
!
!-----------------------------------------------------------------------
!
IMPLICIT none
INTEGER, INTENT (IN ) :: mnode, mgaus
INTEGER, INTENT (OUT) :: ndime, nnode, ngaus, nnodb, ngaub
DOUBLE PRECISION, INTENT (OUT) :: eledb (4, mnode+1, mgaus), &
elbdb (3, mnode+1, mgaus)
CHARACTER (LEN=10), INTENT (OUT) :: eltyp
!
! First reads the element type the finds element information
!
READ (1, '(a10)', err = 10) eltyp
!
! Three noded linear triangle
!
IF ( eltyp (1:5) == 'tria3' ) then
ndime = 2 ; nnode = 3 ; ngaus = 1 ; nnodb = 2 ; ngaub = 1
CALL tria3db (eledb) ; CALL lin2db (elbdb)
!
! Six noded quadratic triangle
!
ELSEIF ( eltyp (1:5) == 'tria6' ) then
ndime = 2 ; nnode = 6 ; ngaus = 3 ; nnodb = 3 ; ngaub = 2
CALL tria6db (eledb) ; CALL qua3db (elbdb)
!
! Four noded bi-linear quadrilateral
!
ELSEIF ( eltyp (1:5) == 'quad4' ) then
ndime = 2 ; nnode = 4 ; ngaus = 4 ; nnodb = 2 ; ngaub = 1
CALL quad4db (eledb) ; CALL lin2db (elbdb)
!
! Four noded linear tetrahedron
!
ELSEIF ( eltyp (1:5) == 'tetr4' ) then
ndime = 3 ; nnode = 4 ; ngaus = 1 ; nnodb = 3 ; ngaub = 1
CALL tetr4db (eledb) ; CALL tria3db (elbdb)
!
! Ten noded quadratic tetrahedron
!
ELSEIF ( eltyp (1:6) == 'tetr10' ) then
ndime = 3 ; nnode = 10 ; ngaus = 5 ; nnodb = 6 ; ngaub = 3
CALL tetr10db (eledb) ; CALL tria6db (elbdb)
!
! Tri-linear brick element
!
ELSEIF ( eltyp (1:5) == 'hexa8' ) then
ndime = 3 ; nnode = 8 ; ngaus = 8 ; nnodb = 4 ; ngaub = 4
CALL hexa8db (eledb) ; CALL quad4db (elbdb)
!
! No more elements implemented yet
!
ELSE
WRITE (6, '(a)') ' Unknown element type'
STOP 'elinfo: Unknown element type'
END IF
!
! Checks that the maximum dimensions are OK
!
IF ( ngaus > mgaus ) then
WRITE (6, 100) ngaus
100 FORMAT(' problem dimensions exceed maximum, ', &
& 'set mgaus to:',i10)
STOP 'elinfo: dimensions exceed mgaus'
ELSEIF ( nnode > mnode ) then
WRITE (6, 101) nnode
101 FORMAT(' problem dimensions exceed maximum, ', &
& 'set mnode to:',i10)
STOP 'elinfo: dimensions exceed mnode'
END IF
RETURN
!
10 WRITE (6, '(a)') ' Error reading element type'
STOP 'elinfo: Error reading element type'
END SUBROUTINE elinfo
!-----------------------------------------------------------------------
SUBROUTINE lin2db (eledb)
!-----------------------------------------------------------------------
!
! Two noded linear element
!
!-----------------------------------------------------------------------
!
IMPLICIT DOUBLE PRECISION (a - h, o - z)
DIMENSION eledb (2, 3, 1)
!
eledb (1, 1, 1) = 0.5d0
eledb (2, 1, 1) = - 0.5d0
eledb (1, 2, 1) = 0.5d0
eledb (2, 2, 1) = 0.5d0
eledb (1, 3, 1) = 2.d0
eledb (2, 3, 1) = 0.d0
END SUBROUTINE lin2db
!-----------------------------------------------------------------------
SUBROUTINE qua3db (eledb)
!-----------------------------------------------------------------------
!
! Three noded quadratic element
!
!-----------------------------------------------------------------------
!
IMPLICIT DOUBLE PRECISION (a - h, o - z)
DIMENSION eledb (2, 4, 2), gauss (2)
DATA gauss / - 0.577350269189626d0, 0.577350269189626d0 /
!
shap1 (xi) = xi * (xi - 1) / 2.d0 ! obsolete under f95
shap2 (xi) = 1 - xi * xi
shap3 (xi) = xi * (xi + 1) / 2.d0
dshp1 (xi) = xi - 0.5d0
dshp2 (xi) = - 2 * xi
dshp3 (xi) = xi + 0.5d0
!
DO ig = 1, 2
xi = gauss (ig)
eledb (1, 1, ig) = shap1 (xi)
eledb (2, 1, ig) = dshp1 (xi)
eledb (1, 2, ig) = shap2 (xi)
eledb (2, 2, ig) = dshp2 (xi)
eledb (1, 3, ig) = shap3 (xi)
eledb (2, 3, ig) = dshp3 (xi)
eledb (1, 4, ig) = 1.d0
eledb (2, 4, ig) = xi
END DO
END SUBROUTINE qua3db
!-----------------------------------------------------------------------
SUBROUTINE tria3db (eledb)
!-----------------------------------------------------------------------
!
! Three noded linear triangle
!
!-----------------------------------------------------------------------
!
IMPLICIT DOUBLE PRECISION (a - h, o - z)
DIMENSION eledb (3, 4, 1)
!
! Gauss point and weight
!
eledb (1, 4, 1) = 1.d0 / 2.d0
eledb (2, 4, 1) = 1.d0 / 3.d0
eledb (3, 4, 1) = 1.d0 / 3.d0
!
! shape function and derivatives at gauss point
!
eledb (1, 1, 1) = 1.d0 / 3.d0
eledb (2, 1, 1) = - 1.d0
eledb (3, 1, 1) = - 1.d0
eledb (1, 2, 1) = 1.d0 / 3.d0
eledb (2, 2, 1) = 1.d0
eledb (3, 2, 1) = 0.d0
eledb (1, 3, 1) = 1.d0 / 3.d0
eledb (2, 3, 1) = 0.d0
eledb (3, 3, 1) = 1.d0
END SUBROUTINE tria3db
!-----------------------------------------------------------------------
SUBROUTINE tria6db (eledb)
!-----------------------------------------------------------------------
!
! Six noded linear triangle with midside integration
!
!-----------------------------------------------------------------------
!
IMPLICIT DOUBLE PRECISION (a - h, o - z)
DIMENSION eledb (3, 7, 3), gauss (2, 3)
DATA (gauss (i, 1), i = 1, 2) / 0.5d0, 0.0d0 /
DATA (gauss (i, 2), i = 1, 2) / 0.5d0, 0.5d0 /
DATA (gauss (i, 3), i = 1, 2) / 0.0d0, 0.5d0 /
!
! Shape functions and derivatives ! obsolete under f95
!
shap1 (xi, et) = (xi + et - 1) * (2 * xi + 2 * et - 1)
shap2 (xi, et) = 4 * xi * (1 - xi - et)
shap3 (xi, et) = xi * (2 * xi - 1)
shap4 (xi, et) = 4 * xi * et
shap5 (xi, et) = et * (2 * et - 1)
shap6 (xi, et) = 4 * et * (1 - xi - et)
dsh1x (xi, et) = 4 * xi + 4 * et - 3
dsh2x (xi, et) = 4 - 8 * xi - 4 * et
dsh3x (xi, et) = 4 * xi - 1
dsh4x (xi, et) = 4 * et
dsh5x (xi, et) = 0
dsh6x (xi, et) = - 4 * et
dsh1e (xi, et) = 4 * xi + 4 * et - 3
dsh2e (xi, et) = - 4 * xi
dsh3e (xi, et) = 0
dsh4e (xi, et) = 4 * xi
dsh5e (xi, et) = 4 * et - 1
dsh6e (xi, et) = 4 - 4 * xi - 8 * et
!
! Constructs the data base. First Gauss point position and weight
!
DO ig = 1, 3
eledb (1, 7, ig) = 1.d0 / 6.d0
xi = gauss (1, ig)
et = gauss (2, ig)
eledb (2, 7, ig) = xi
eledb (3, 7, ig) = et
!
! Shape functions and derivatives
!
eledb (1, 1, ig) = shap1 (xi, et)
eledb (2, 1, ig) = dsh1x (xi, et)
eledb (3, 1, ig) = dsh1e (xi, et)
eledb (1, 2, ig) = shap2 (xi, et)
eledb (2, 2, ig) = dsh2x (xi, et)
eledb (3, 2, ig) = dsh2e (xi, et)
eledb (1, 3, ig) = shap3 (xi, et)
eledb (2, 3, ig) = dsh3x (xi, et)
eledb (3, 3, ig) = dsh3e (xi, et)
eledb (1, 4, ig) = shap4 (xi, et)
eledb (2, 4, ig) = dsh4x (xi, et)
eledb (3, 4, ig) = dsh4e (xi, et)
eledb (1, 5, ig) = shap5 (xi, et)
eledb (2, 5, ig) = dsh5x (xi, et)
eledb (3, 5, ig) = dsh5e (xi, et)
eledb (1, 6, ig) = shap6 (xi, et)
eledb (2, 6, ig) = dsh6x (xi, et)
eledb (3, 6, ig) = dsh6e (xi, et)
END DO
END SUBROUTINE tria6db
!-----------------------------------------------------------------------
SUBROUTINE tetr4db (eledb)
!-----------------------------------------------------------------------
!
! Four noded linear tetrahedron
!
!-----------------------------------------------------------------------
!
IMPLICIT DOUBLE PRECISION (a - h, o - z)
DIMENSION eledb (4, 5, 1)
!
! Gauss point and weight
!
eledb (1, 5, 1) = 1.d0 / 6.d0
eledb (2, 5, 1) = 1.d0 / 4.d0
eledb (3, 5, 1) = 1.d0 / 4.d0
eledb (4, 5, 1) = 1.d0 / 4.d0
!
! shape function and derivatives at gauss point
!
eledb (1, 1, 1) = 1.d0 / 4.d0
eledb (2, 1, 1) = - 1.d0
eledb (3, 1, 1) = - 1.d0
eledb (4, 1, 1) = - 1.d0
eledb (1, 2, 1) = 1.d0 / 4.d0
eledb (2, 2, 1) = 1.d0
eledb (3, 2, 1) = 0.d0
eledb (4, 2, 1) = 0.d0
eledb (1, 3, 1) = 1.d0 / 4.d0
eledb (2, 3, 1) = 0.d0
eledb (3, 3, 1) = 1.d0
eledb (4, 3, 1) = 0.d0
eledb (1, 4, 1) = 1.d0 / 4.d0
eledb (2, 4, 1) = 0.d0
eledb (3, 4, 1) = 0.d0
eledb (4, 4, 1) = 1.d0
END SUBROUTINE tetr4db
!-----------------------------------------------------------------------
SUBROUTINE tetr10db (eledb)
!-----------------------------------------------------------------------
!
! Ten noded quadratic tetrahedron
!
!-----------------------------------------------------------------------
!
IMPLICIT DOUBLE PRECISION (a - h, o - z)
PARAMETER (sixth = 0.166666666666667d0)
DIMENSION eledb (4, 11, 5), gauss (3, 5)
!
! Gauss point positions
!
DATA (gauss (i, 1), i = 1, 3) / 0.25d0, 0.25d0, 0.25d0 /
DATA (gauss (i, 2), i = 1, 3) / sixth, sixth, sixth /
DATA (gauss (i, 3), i = 1, 3) / 0.5d0, sixth, sixth /
DATA (gauss (i, 4), i = 1, 3) / sixth, 0.5d0, sixth /
DATA (gauss (i, 5), i = 1, 3) / sixth, sixth, 0.5d0 /
!
! Shape functions and derivatives ! obsolete under f95
!
shap1 (xi, et, ze) = (1 - 2 * xi - 2 * et - 2 * ze) &
* (1 - xi - et - ze)
shap2 (xi, et, ze) = xi * (2 * xi - 1)
shap3 (xi, et, ze) = et * (2 * et - 1)
shap4 (xi, et, ze) = ze * (2 * ze - 1)
shap5 (xi, et, ze) = 4 * xi * (1 - xi - et - ze)
shap6 (xi, et, ze) = 4 * xi * et
shap7 (xi, et, ze) = 4 * et * (1 - xi - et - ze)
shap8 (xi, et, ze) = 4 * ze * (1 - xi - et - ze)
shap9 (xi, et, ze) = 4 * xi * ze
shap0 (xi, et, ze) = 4 * ze * et
!
dsh1x (xi, et, ze) = 4 * (xi + et + ze) - 3
dsh2x (xi, et, ze) = 4 * xi - 1
dsh3x (xi, et, ze) = 0
dsh4x (xi, et, ze) = 0
dsh5x (xi, et, ze) = 4 * (1 - 2 * xi - et - ze)
dsh6x (xi, et, ze) = 4 * et
dsh7x (xi, et, ze) = - 4 * et
dsh8x (xi, et, ze) = - 4 * ze
dsh9x (xi, et, ze) = 4 * ze
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?