📄 elmlib.f90
字号:
SUBROUTINE ElmLib( diags, vforc, iswth )
!........
! 模块功能
! 形成总体矩阵: 刚度阵, 阻尼阵, 质量阵.
!........
! iswth 控制开关
! == 0 读入材料信息 == 1 单元刚度阵
! == 2 单元质量阵 == 3 单元阻尼阵
! == 4 单元几何阵 == 5 内力等效节点力
! == 6 体积力等效节点力 == 7 单元高斯点应力
! ==-7 单元节点应力 == 8 单元信息
! == 9 网格输出 ==10 节点反力
!........
USE CtrlData
USE MeshData
USE GlobData
USE SolvData
USE FrontData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION diags( ndofs ), vforc( ndofs )
kcoun = 0 !标识进展
IF( iswth .EQ. 4 ) isavc = lsavc( 2 )
IF( iswth .LT. 4 ) isavc = lsavc( iswth )
IF( isavc .EQ. 0 .AND. iswth .LE. 4 ) RETURN
!.......初始化工作变量和波阵.
knume = 0
npntw = 0
smaxv = 0.0D0
CALL InitInteger( lpntw, mpntw, 0 )
CALL InitFloat( stifw, mstfw, 0.0D0 )
CALL InitFloat( vforc, ndofs, 0.0D0 )
!.......根据记录信息组装矩阵.
DO ircrd = 1, nrcrd
inump = lrcrd( 1, ircrd )
jnump = lrcrd( 2, ircrd )
IF( isavc .EQ. 1 ) kdata = 0
IF( isavc .LT. 0 ) kdata =-isavc - 1
DO knump = inump, jnump
kpoin = lbcks( 2, knump )
kpntw = lbcks( 3, knump )
kpnte = lpnte( kpoin )
DO ipnte = 1, kpnte
knume = knume + 1
kcoun = kcoun + 1
kelem = lopts( knume )
IF( necho .NE. 0 ) CALL ShowProcess( kcoun, nelem )
!.............形成当步波阵节点序号与节点总体编号对照表.
!.............形成元素节点序号与波阵节点序号对照表.
DO inode = 1, mnode
lwave( inode ) = 0
jpoin = IABS( lnods( inode, kelem ) )
IF( jpoin .GT. 0 ) THEN
!.................当前节点为非零节点,查是否在波阵中.
jpoin = llnks( jpoin )
DO jpntw = 1, mpntw
IF( lwave( inode ) .EQ. 0 ) THEN
IF( jpoin .EQ. lpntw( jpntw ) ) &
lwave( inode ) = jpntw
END IF
END DO
!.................如果当前节点不在波阵中,添加该节点.
IF( lwave( inode ) .EQ. 0 ) THEN
npntw = npntw + 1
jpntw = 1
DO WHILE( lwave( inode ) .EQ. 0 )
IF( lpntw( jpntw ) .EQ. 0 ) THEN
lpntw( jpntw ) = jpoin
lwave( inode ) = jpntw
ELSE
jpntw = jpntw + 1
END IF
END DO
END IF
END IF
END DO
!.............形成元素矩阵并组装成总体矩阵.
CALL ElmOpt( kelem, iswth )
CALL slantb
CALL Assemb( vforc, iswth, kelem )
IF( nerrc .GT. 0 ) RETURN
END DO
IF( npntw .GT. 0 .AND. iswth .LE. 4 ) THEN
!.............存储总体矩阵对角元素.
ndofw = npntw * mdofn
DO idofn = 1, mdofn
idofw = ( kpntw - 1 ) * mdofn + idofn
idofs = ( kpoin - 1 ) * mdofn + idofn
istfw = idofw * ( idofw + 1 ) / 2
diags( idofs ) = stifw( istfw )
IF( stifw( istfw ) .GT. smaxv ) &
smaxv = stifw( istfw )
stifw( istfw ) = 0.0D0
END DO
!.............存储总体矩阵元素并收缩波阵
IF( isavc .NE. 2 ) THEN
DO idofn = 1, mdofn
locat = kdata
kdofw = ( kpntw - 1 ) * mdofn
idofw = kdofw + idofn
DO jpntw = 1, mpntw
IF( lpntw( jpntw ) .NE. 0 ) THEN
DO jdofn = 1, mdofn
jdofw = ( jpntw - 1 ) * mdofn + jdofn
IF( jdofw .LE. kdofw ) THEN
!.........................当前行元素存储.
locat = locat + 1
istfw = idofw * ( idofw - 1 ) / 2 + jdofw
buffs( locat ) = stifw( istfw )
stifw( istfw ) = 0.0D0
ELSE IF( jdofw .GT. idofw ) THEN
!.........................当前列元素存储.
locat = locat + 1
istfw = jdofw * ( jdofw - 1 ) / 2 + idofw
buffs( locat ) = stifw( istfw )
stifw( istfw ) = 0.0D0
END IF
END DO
END IF
END DO
kdata = kdata + ndofw - idofn
END DO
END IF
END IF
!...........波阵.
IF( npntw .GT. 0 ) THEN
npntw = npntw - 1
lpntw( kpntw ) = 0
END IF
END DO
IF( isavc .EQ. 1 ) &
CALL OptRec( buffs, ircrd, 2, iswth, nexch )
END DO
RETURN
END
SUBROUTINE ElmOpt( ielem, iswth )
!........
! 模块功能
! 选择单元类型, 计算单元矩阵或应力.
!........
USE CtrlData
USE MeshData
USE ElmtData
USE GlobData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
IF( iswth .EQ. 0 ) THEN
imats = ielem
ELSE IF( iswth .EQ. 8 ) THEN
imats = ielem
nnode = mnode
ELSE IF( iswth .EQ. 9 ) THEN
nnode = mnode
iuses = luses( ielem )
imats = lmats( ielem )
DO inode = 1, mnode
lnode( inode ) = lnods( inode, ielem )
END DO
ELSE
!.........单元材料号和单元类型.
CALL ElmPar( ielem )
iuses = luses( ielem )
IF( nerrc .GT. 0 ) RETURN
!.........设定单元信息.
CALL ElmInf( ielem, iswth )
IF( nerrc .GT. 0 ) RETURN
imats = lmats( ielem )
END IF
!.......调用相应的单元.
itype = IABS( lprps( 1, imats ) )
SELECT CASE( itype )
CASE( 1 )
CALL Elmt001( ielem, imats, iswth, iuses )
CASE( 2 )
CALL Elmt002( ielem, imats, iswth, iuses )
CASE( 3 )
CALL Elmt003( ielem, imats, iswth, iuses )
CASE( 4 )
CALL Elmt004( ielem, imats, iswth, iuses )
CASE( 200 )
CALL Elmt200( ielem, imats, iswth, iuses )
CASE( 300 )
CALL Elmt300( ielem, imats, iswth, iuses )
CASE( 400 )
CALL Elmt400( ielem, imats, iswth, iuses )
CASE( 401 )
CALL Elmt401( ielem, imats, iswth, iuses )
CASE( 501 )
CALL Elmt501( smaxv, ielem, imats, iswth, iuses )
CASE DEFAULT
nerrc = 24245
WRITE( 12, 2000 ) itype
END SELECT
2000 FORMAT( //2x, '错误: 无效的单元类型', I3// )
RETURN
END
SUBROUTINE SlantB
!........
! 模块功能
! 对单元刚度矩阵进行斜约束处理.
!........
USE CtrlData
USE MeshData
USE ElmtData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
ALLOCATABLE tmatx(:,:), vmatx(:,:)
IF( nslnt .LE. 0 ) RETURN
ALLOCATE( tmatx( mdofn, mdofn ), vmatx( mdofn, mdofn ), &
STAT = nerrc )
IF( nerrc .NE. 0 ) nerrc = 24435
IF( nerrc .GT. 0 ) RETURN
DO inode = 1, nnode
ipoin = IABS( lnode( inode ) )
IF( ipoin .GT. 0 ) THEN
label = 0
DO idofn = 1, mdofn
DO jdofn = 1, mdofn
tmatx( idofn, jdofn ) = 0.0D0
END DO
tmatx( idofn, idofn ) = 1.0D0
END DO
DO islnt = 1, nslnt
jpoin = lslnt( 1, islnt )
jdofn = lslnt( 2, islnt )
IF( ipoin .EQ. jpoin ) THEN
label = label + 1
DO idofn = 1, mdofn
tmatx( jdofn, idofn ) = vslnt( idofn, islnt )
END DO
ELSE IF( jpoin .LT. 0 ) THEN
idime =-jpoin
IF( DABS( coord( idime, ipoin ) ) .LT. 1.0D-10 ) THEN
label = label + 1
DO idofn = 1, mdofn
tmatx( jdofn, idofn ) = vslnt( idofn, islnt )
END DO
END IF
END IF
END DO
IF( label .GT. 0 ) THEN
CALL Invers( tmatx, mdofn, mdofn, nerrc )
IF( nerrc .GT. 0 ) RETURN
DO jnode = 1, nnode
DO idofn = 1, ndofn
DO jdofn = 1, ndofn
vmatx( idofn, jdofn ) = 0.0D0
DO kdofn = 1, ndofn
iloca = ( jnode - 1 ) * ndofn + idofn
jloca = ( inode - 1 ) * ndofn + kdofn
vmatx( idofn, jdofn ) = vmatx( idofn, jdofn ) + &
stife( iloca, jloca ) * tmatx( kdofn, jdofn )
END DO
END DO
END DO
DO idofn = 1, ndofn
DO jdofn = 1, ndofn
iloca = ( jnode - 1 ) * ndofn + idofn
jloca = ( inode - 1 ) * ndofn + jdofn
stife( iloca, jloca ) = vmatx( idofn, jdofn )
END DO
END DO
END DO
DO jnode = 1, nnode
DO idofn = 1, ndofn
DO jdofn = 1, ndofn
vmatx( idofn, jdofn ) = 0.0D0
DO kdofn = 1, ndofn
iloca = ( inode - 1 ) * ndofn + kdofn
jloca = ( jnode - 1 ) * ndofn + jdofn
vmatx( idofn, jdofn ) = vmatx( idofn, jdofn ) + &
tmatx( kdofn, idofn ) * stife( iloca, jloca )
END DO
END DO
END DO
DO idofn = 1, ndofn
DO jdofn = 1, ndofn
iloca = ( inode - 1 ) * ndofn + idofn
jloca = ( jnode - 1 ) * ndofn + jdofn
stife( iloca, jloca ) = vmatx( idofn, jdofn )
END DO
END DO
END DO
DO idofn = 1, ndofn
vmatx( idofn, 1 ) = 0.0D0
DO jdofn = 1, ndofn
iloca = ( inode - 1 ) * ndofn + jdofn
vmatx( idofn, 1 ) = vmatx( idofn, 1 ) + &
tmatx( jdofn, idofn ) * force( iloca )
END DO
END DO
DO idofn = 1, ndofn
iloca = ( inode - 1 ) * ndofn + idofn
force( iloca ) = vmatx( idofn, 1 )
END DO
END IF
END IF
END DO
DEALLOCATE( vmatx, tmatx )
RETURN
END
SUBROUTINE Assemb( vforc, iswth, ielem )
USE CtrlData
USE MeshData
USE ElmtData
USE FrontData
!........
! 模块功能
! 对当步元素刚度矩阵和元素分布载荷的等价节点力组装到波阵下三
! 角区和总载荷向量.
!........
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION vforc( ndofs )
imats = lmats( ielem )
DO inode = 1, nnode
ipntw = lwave( inode )
IF( ipntw .NE. 0 ) THEN
ipoin = lpntw( ipntw )
DO idofn = 1, ndofn
kdofn = lprps( idofn + 2, imats )
IF( kdofn .GT. 0 ) THEN
idofs = ( ipoin - 1 ) * mdofn + kdofn
idofw = ( ipntw - 1 ) * mdofn + kdofn
idofe = ( inode - 1 ) * ndofn + idofn
vforc( idofs ) = vforc( idofs ) + force( idofe )
!...............迭加节点载荷.
IF( iswth .LE. 4 ) THEN
DO jnode = 1, nnode
jpntw = lwave( jnode )
IF( jpntw .NE. 0 ) THEN
jpoin = lpntw( jpntw )
DO jdofn = 1, ndofn
kdofn = lprps( jdofn + 2, imats )
IF( kdofn .GT. 0 ) THEN
jdofw = ( jpntw - 1 ) * mdofn + kdofn
jdofe = ( jnode - 1 ) * ndofn + jdofn
IF( jdofw .LE. idofw ) THEN
!...........................迭加单元刚度.
istfw = ( idofw - 1 ) * idofw / 2 + jdofw
stifw( istfw ) = stifw( istfw ) + &
stife( idofe, jdofe )
END IF
END IF
END DO
END IF
END DO
END IF
END IF
END DO
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -