📄 elmt300.f90
字号:
wzgas( izgas ) * xjaco
CALL Dmat300( dmatx, imats, corep, xlocd, ylocd, &
ielem, igeob, imatb, 0 )
IF( nerrc .GT. 0 ) RETURN
CALL Strn300( corgp, strng, shape, nlinc, 3 )
CALL Strs300( corep, strng, strgp, dmatx, xlocd, &
ylocd, imtyp, kgaus, kelem, imatb, &
ietyp, 0 )
IF( iswth .EQ. 7 ) THEN
WRITE( 12, 2400 ) corgp(1), corgp(2), strgp(1), &
strgp(2), strgp(3), corgp(3), &
strgp(4), strgp(5), strgp(6)
ELSE
iloca = 0
jflag = 1
xloca = pxgas( ixgas )
yloca = pygas( iygas )
zloca = pzgas( izgas )
CALL Shap3D( xloca, yloca, zloca, coren, shape, &
xjaco, ndime, nnode, lnode, ielem, &
jflag, nerrc )
DO inode = 1, nnode
IF( lnode( inode ) .NE. 0 ) THEN
iloca = iloca + 1
jloca = 0
DO jnode = 1, nnode
IF( lnode( jnode ) .NE. 0 ) THEN
jloca = jloca + 1
stife(iloca, jloca) = stife(iloca, jloca) + &
shape( 4, inode) * shape( 4, jnode)
END IF
END DO
DO istre = 1, 6
snode(istre, iloca) = snode(istre, iloca) + &
shape( 4, inode) * strgp(istre )
smatx(istre, iloca) = smatx(istre, iloca) + &
shape( 4, inode) * strng(istre )
END DO
END IF
END DO
END IF
END DO
END DO
END DO
IF( iswth .EQ. mswth + 7 ) THEN
knode = 0
DO inode = 1, nnode
IF( lnode( inode ) .NE. 0 ) knode = knode + 1
END DO
CALL Invers( stife, mdofe, knode, nerrc )
IF( nerrc .NE. 0 ) RETURN
DO istre = 1, 6
iloca = 0
DO inode = 1, nnode
IF( lnode( inode ) .NE. 0 ) THEN
force( istre ) = 0.0D0
strng( istre ) = 0.0D0
iloca = iloca + 1
jloca = 0
DO jnode = 1, nnode
IF( lnode( jnode ) .NE. 0 ) THEN
jloca = jloca + 1
force( istre ) = force( istre ) + &
stife( iloca, jloca ) * snode( istre, jloca )
strng( istre ) = strng( istre ) + &
stife( iloca, jloca ) * smatx( istre, jloca )
END IF
END DO
ipoin = lnodx( inode, kelem )
IF( ipoin .LT. 0 ) THEN
ipoin =-ipoin
strss( istre, ipoin ) = strss( istre, ipoin ) + &
force( istre ) * volum
strns( istre, ipoin ) = strns( istre, ipoin ) + &
strng( istre ) * volum
IF( istre .EQ. 1 ) THEN
strss( 7, ipoin ) = strss( 7, ipoin ) + volum
strns( 7, ipoin ) = strns( 7, ipoin ) + volum
END IF
ELSE
strsx( istre, ipoin ) = strsx( istre, ipoin ) + &
force( istre ) * volum
strnx( istre, ipoin ) = strnx( istre, ipoin ) + &
strng( istre ) * volum
IF( istre .EQ. 1 ) THEN
strsx( 7, ipoin ) = strsx( 7, ipoin ) + volum
strnx( 7, ipoin ) = strnx( 7, ipoin ) + volum
END IF
END IF
END IF
END DO
END DO
END IF
kelem = kelem + 1
END DO
ELSE IF( iswth .EQ. 6 ) THEN
DO ibloc = 1, nbloc
kgaus = 0
igeob = props( ibloc * 2 + 2, imats ) + 0.5D0
imatb = props( ibloc * 2 + 3, imats ) + 0.5D0
ietyp = props( 2, igeob ) + 0.5D0
nxgas = props( 5, igeob ) + 0.5D0
nygas = props( 6, igeob ) + 0.5D0
nzgas = props( 7, igeob ) + 0.5D0
nbody = props( 15, imatb ) + 0.5D0
imtyp = props( 16, imatb ) + 0.5D0
CALL Factor( nbody )
IF( nerrc .GT. 0 ) RETURN
fbody( 1 ) = props( 12, imatb ) * fctor
fbody( 2 ) = props( 13, imatb ) * fctor
fbody( 3 ) = props( 14, imatb ) * fctor
summy = DABS( fbody( 1 ) ) + DABS( fbody( 2 ) ) + &
DABS( fbody( 3 ) )
IF( summy .GT. 1.0D-10 ) THEN
CALL Bloc300( coreb, lnodb, igeob, nnodb )
CALL Gauss( pxgas, wxgas, nxgas, nerrc )
CALL Gauss( pygas, wygas, nygas, nerrc )
CALL Gauss( pzgas, wzgas, nzgas, nerrc )
IF( nerrc .GT. 0 ) RETURN
DO ixgas = 1, nxgas
DO iygas = 1, nygas
DO izgas = 1, nzgas
kgaus = kgaus + 1
xloca = pxgas( ixgas )
yloca = pygas( iygas )
zloca = pzgas( izgas )
CALL Shape300( shape, coreb, lnodb, xloca, yloca, &
zloca, xjaco, nnodb, ielem, iflag )
IF( nerrc .GT. 0 ) RETURN
dvolu = wxgas( ixgas ) * wygas( iygas ) * &
wzgas( izgas ) * xjaco
DO inode = 1, nnode
DO idofn = 1, ndofn
idofe = ( inode - 1 ) * ndofn + idofn
force( idofe ) = force( idofe ) + dvolu * &
fbody( idofn ) * shape( 4, inode )
END DO
END DO
END DO
END DO
END DO
END IF
END DO
ELSE IF( iswth .EQ. 11 ) THEN
!.........计算单元历史量
kelem = lprtx( ielem )
DO ibloc = 1, nbloc
kgaus = 0
igeob = props( ibloc * 2 + 2, imats ) + 0.5D0
imatb = props( ibloc * 2 + 3, imats ) + 0.5D0
ietyp = props( 2, igeob ) + 0.5D0
nxgas = props( 5, igeob ) + 0.5D0
nygas = props( 6, igeob ) + 0.5D0
nzgas = props( 7, igeob ) + 0.5D0
imtyp = props( 16, imatb ) + 0.5D0
IF( imtyp .EQ. 1 ) THEN
CALL Bloc300( coreb, lnodb, igeob, nnodb )
CALL MidPln300( coreb, corep, lnodb, ielem, nnodb )
CALL Gauss( pxgas, wxgas, nxgas, nerrc )
CALL Gauss( pygas, wygas, nygas, nerrc )
CALL Gauss( pzgas, wzgas, nzgas, nerrc )
IF( nerrc .NE. 0 ) RETURN
DO ixgas = 1, nxgas
DO iygas = 1, nygas
DO izgas = 1, nzgas
kgaus = kgaus + 1
xloca = pxgas( ixgas )
yloca = pygas( iygas )
zloca = pzgas( izgas )
xlocd = pxgas( ixgas )
ylocd = pygas( iygas )
CALL Shape300( shape, coreb, lnodb, xloca, yloca, &
zloca, xjaco, nnodb, ielem, iflag )
IF( nerrc .GT. 0 ) RETURN
CALL Hstr300( corep, xlocd, ylocd, shape, imats, &
imatb, kgaus, kelem )
END DO
END DO
END DO
END IF
kelem = kelem + 1
END DO
END IF
2200 FORMAT( /20x, ' 单元 = ', i5 / 7x, 'x', 10x, 'y', &
10x, 'S-xx', 10x, 'S-yy', 10x, 'S-zz'/ 5x, 'z(h)', &
20x, 'S-yz', 10x, 'S-zx', 10x, 'S-xy ' )
2400 FORMAT( 2X, 2F10.3, 3E14.6 / 2X, F10.3, 10X, 3E14.6 )
RETURN
END
SUBROUTINE RdMat300( imats )
USE CtrlData
USE MeshData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION lwork( 12 ), lnodb( 20 ), fwork( 3 )
iswth = lprps( mdofn + 3, imats )
IF( iswth .EQ. 0 .OR. iswth .EQ. 1 ) THEN
READ( 11, 2000 ) nbloc, nlinc
IF( nlinc .NE. 0 ) nlinc = 1
props( 2, imats ) = nbloc
props( 3, imats ) = nlinc
IF( nlinc .EQ. 1 ) THEN
IF( nincc .EQ. 1 .OR. nincc .EQ. 3 ) nerrc = 38556
IF( nerrc .NE. 0 ) RETURN
nincc = 2
END IF
IF( nprnc .NE. 0 ) THEN
WRITE( 12, 2200 ) nbloc
IF( nlinc .EQ. 0 ) WRITE( 12, 2400 )
IF( nlinc .GT. 0 ) WRITE( 12, 2401 )
END IF
IF( nbloc * 2 + 3 .GT. nprop ) THEN
WRITE( 12, 2600 ) nbloc * 2 + 3
nerrc = 49456
RETURN
END IF
!.........读入各子块的几何特性序号和材料特性序号
IF( nprnc .NE. 0 ) WRITE( 12, 2800 )
DO iline = 1, ( nbloc + 5 ) / 6
READ( 11, 2000 ) lwork
jbloc = iline * 6
ibloc = jbloc - 5
IF( jbloc .GT. nbloc ) jbloc = nbloc
DO kbloc = ibloc, jbloc
iloca = ( kbloc - ibloc ) * 2
igeom = lwork( iloca + 1 )
imate = lwork( iloca + 2 )
props( kbloc * 2 + 2, imats ) = igeom
props( kbloc * 2 + 3, imats ) = imate
IF( nprnc .NE. 0 ) WRITE( 12, 3000 ) kbloc, igeom, imate
IF( igeom .LT. 1 .OR. igeom .GT. nmats ) THEN
WRITE( 12, 3200 ) igeom
nerrc = 3454
RETURN
END IF
IF( imate .LT. 1 .OR. imate .GT. nmats ) THEN
WRITE( 12, 3400 ) imate
nerrc = 3455
RETURN
END IF
lprps( 1, igeom ) = 300
lprps( mdofn + 3, igeom ) = 2
lprps( 1, imate ) = 300
lprps( mdofn + 3, imate ) = 3
END DO
END DO
ELSE IF( iswth .EQ. 2 ) THEN
READ( 11, 2000 ) ietyp, npntb, nnodb, nxgas, nygas, nzgas
IF( nxgas .LE. 0 ) nxgas = 3
IF( nygas .LE. 0 ) nygas = nxgas
IF( nzgas .LE. 0 ) nzgas = nxgas
nstrh = MAX( nstrh, nxgas * nygas * nzgas * 6 )
props( 2, imats ) = ietyp
props( 3, imats ) = npntb
props( 4, imats ) = nnodb
props( 5, imats ) = nxgas
props( 6, imats ) = nygas
props( 7, imats ) = nzgas
IF( nprnc .NE. 0 ) THEN
WRITE( 12, 3600 )
IF( ietyp .EQ. 0 ) WRITE( 12, 3800 )
IF( ietyp .EQ. 1 ) WRITE( 12, 3801 )
IF( ietyp .EQ. 2 ) WRITE( 12, 3802 )
IF( ietyp .EQ. 3 ) WRITE( 12, 3803 )
IF( ietyp .EQ. 4 ) WRITE( 12, 3804 )
IF( ietyp .EQ. 5 ) WRITE( 12, 3805 )
WRITE( 12, 4000 ) nxgas, nygas, nzgas
END IF
IF( npntb * 3 + nnodb + 3 .GT. nprop ) THEN
WRITE( 12, 2600 ) npntb * 3 + nnodb + 7
nerrc = 49456
RETURN
END IF
!.........读入定义子块的节点局部坐标
DO ipntb = 1, npntb
READ( 11, 4200 ) fwork
iloca = ipntb * 3 + 4
props( iloca + 1, imats ) = fwork( 1 )
props( iloca + 2, imats ) = fwork( 2 )
props( iloca + 3, imats ) = fwork( 3 )
IF( nprnc .NE. 0 ) WRITE( 12, 4400 ) ipntb, fwork
END DO
READ( 11, 4600 ) ( lnodb( inodb ), inodb = 1, nnodb )
IF( nprnc .NE. 0 ) WRITE( 12, 4800 ) ( lnodb( inodb ), &
inodb = 1, nnodb )
iloca = npntb * 3 + 7
DO inodb = 1, nnodb
props( iloca + inodb, imats ) = lnodb( inodb )
IF( lnodb( inodb ) .GT. npntb ) THEN
WRITE( 12, 5000 ) lnodb( inodb )
nerrc = 65967
RETURN
END IF
END DO
ELSE IF( iswth .EQ. 3 ) THEN
READ( 11, 5200 ) qmt11, qmt22, qmt33, qmt12, qmt23, qmt13, &
qmt44, qmt55, qmt66, denst, xbody, ybody, &
zbody, nbody, imtyp, ngaus
IF( nprop .LT. 16 ) WRITE( 12, 2600 ) 16
props( 2, imats ) = qmt11
props( 3, imats ) = qmt22
props( 4, imats ) = qmt33
props( 5, imats ) = qmt12
props( 6, imats ) = qmt23
props( 7, imats ) = qmt13
props( 8, imats ) = qmt44
props( 9, imats ) = qmt55
props( 10, imats ) = qmt66
props( 11, imats ) = denst
props( 12, imats ) = xbody
props( 13, imats ) = ybody
props( 14, imats ) = zbody
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -