📄 front.f90
字号:
SUBROUTINE WavOpt
!.......
! 模块功能
! 波前法求解方程时优化波带宽
!.......
! lnods 单元定义(I) mpntw 最大波宽(O)
! lopts 优化后单元顺序(O) nelem 单元总数(I)
! lbcks 回代信息(O) mnode 单元最大节点数(I)
! lbcks( 1, ip ) 主消节点波宽
! lbcks( 2, ip ) 主消节点总体节点号
! lbcks( 3, ip ) 主消节点波阵序号
! lpnte 节点作为主消节点时须入波阵的单元数(O)
!.......
USE CtrlData
USE MeshData
USE FrontData
IMPLICIT DOUBLE PRECISION( a-h, o-z )
DIMENSION vwork( 3 ), vtemp( 3 )
ALLOCATABLE lwork(:)
IF( nerrc .GT. 0 ) RETURN
ALLOCATE( lwork( nelem ), STAT = ierro )
IF( ierro .NE. 0 ) nerrc = 4
IF( nerrc .GT. 0 ) RETURN
CALL ShowMessage( '带宽优化...' )
!.......初始化
mpntw = 0
kpoin = 0
kelem = 0
npntw = 0
DO ielem = 1, nelem
lwork( ielem ) = 0
END DO
DO ipoin = 1, npoin
lposw( ipoin ) = 0
lpntw( ipoin ) = 0
lpnte( ipoin ) = 0
lbcks( 1, ipoin ) = 0
lbcks( 2, ipoin ) = 0
lbcks( 3, ipoin ) = 0
lbcks( 4, ipoin ) = 0
END DO
!.......设定绕节点的单元数.
DO ielem = 1, nelem
DO inode = 1, mnode
ipoin = IABS( lnods( inode, ielem ) )
IF( ipoin .GT. 0 ) THEN
ipoin = llnks( ipoin )
lpnte( ipoin ) = lpnte( ipoin ) + 1
END IF
END DO
END DO
!.......没有使用的节点首先处理
DO ipoin = 1, npoin
IF( lpnte( ipoin ) .EQ. 0 ) THEN
kpoin = kpoin + 1
lbcks( 1, kpoin ) = 0
lbcks( 3, kpoin ) = 0
lbcks( 2, kpoin ) = ipoin
END IF
END DO
!.......寻找绕节点单元数最少的节点作为第一个主消节点.
jpoin = 0
ipnte = nelem + 1
DO WHILE( jpoin .LT. npoin )
jpoin = jpoin + 1
jpnte = lpnte( jpoin )
IF( jpnte .GT. 0 ) THEN
IF( jpnte .LT. ipnte ) THEN
ipoin = jpoin
ipnte = jpnte
DO idime = 1, ndime
vwork( idime ) = 0.0D0
DO jdime = 1, ndime
vwork( idime ) = vwork( idime ) - axisw( idime, jdime ) * ( &
basew( jdime ) - coord( jdime, ipoin ) )
END DO
END DO
ELSE IF( jpnte .EQ. ipnte ) THEN
DO idime = 1, ndime
vtemp( idime ) = 0.0D0
DO jdime = 1, ndime
vtemp( idime ) = vtemp( idime ) - axisw( idime, jdime ) * ( &
basew( jdime ) - coord( jdime, jpoin ) )
END DO
END DO
iflag = 1
DO idime = 1, ndime
IF( iflag .EQ. 1 ) THEN
IF( vtemp( idime ) .LT. vwork( idime ) ) THEN
iflag = 0
ipoin = jpoin
DO jdime = 1, ndime
vwork( jdime ) = vtemp( jdime )
END DO
ELSE IF( vtemp( idime ) .GT. vwork( idime ) ) THEN
iflag = 0
END IF
END IF
END DO
END IF
END IF
END DO
!.......设定进入波阵的单元顺序.
DO WHILE( kpoin .LT. npoin )
kpoin = kpoin + 1
DO ielem = 1, nelem
IF( lwork( ielem ) .EQ. 0 ) THEN
!..........若当前单元不在波阵中，检查当前单元是否有主消节点或当前单元的节点是否全在波阵中.
iflag = 0
jflag = 0
DO inode = 1, mnode
jpoin = IABS( lnods( inode, ielem ) )
IF( jpoin .GT. 0 ) THEN
jpoin = llnks( jpoin )
IF( ipoin .EQ. jpoin ) iflag = iflag + 1
IF( lposw( jpoin ) .EQ. 0 ) jflag = jflag + 1
END IF
END DO
IF( iflag .GT. 0 .OR. jflag .EQ. 0 ) THEN
!...........若当前单元有主消节点或当前单元的节点全在波阵中，则送入波阵并作相应记录.
kelem = kelem + 1
lwork( ielem ) = kelem
lopts( kelem ) = ielem
ipnte = ipnte - iflag + 1
DO inode = 1, mnode
jpoin = IABS( lnods( inode, ielem ) )
IF( jpoin .GT. 0 ) THEN
!.............不为零的节点进入波阵, 查是否已在波阵中.
jpoin = llnks( jpoin )
ipntw = lposw( jpoin )
lpnte( jpoin ) = lpnte( jpoin ) - 1
IF( ipntw .EQ. 0 ) THEN
!..............如果节点不在波阵中, 进入波阵.
iflag = 0
ipntw = 1
npntw = npntw + 1
DO WHILE( iflag .EQ. 0 )
IF( lpntw( ipntw ) .EQ. 0 ) THEN
lpntw( ipntw ) = jpoin
iflag = 1
ELSE
ipntw = ipntw + 1
END IF
END DO
lposw( jpoin ) = ipntw
lbcks( 4, jpoin ) = kpoin
IF( npntw .GT. mpntw ) mpntw = npntw
END IF
END IF
END DO
END IF
END IF
END DO
!........设定主消节点消元时需要进入波阵的单元数.
lpnte( ipoin ) = ipnte
!........设定主消节点回代时的信息.
ipntw = 1
DO WHILE( ipntw .LE. npoin .AND. lpntw(ipntw) .NE. ipoin )
ipntw = ipntw + 1
END DO
lbcks( 1, kpoin ) = npntw
lbcks( 2, kpoin ) = ipoin
lbcks( 3, kpoin ) = ipntw
!........主消节点消元后调整波阵.
npntw = npntw - 1
lpntw( ipntw ) = 0
!........查找是否有可同时消元的节点, 有则进行消元处理.
ipntw = 1
DO ipntw = 1, npoin
ipoin = lpntw( ipntw )
IF( ipoin .GT. 0 ) THEN
IF( lpnte(ipoin) .EQ. 0 ) THEN
kpoin = kpoin + 1
lbcks( 1, kpoin ) = npntw
lbcks( 2, kpoin ) = ipoin
lbcks( 3, kpoin ) = ipntw
lpntw( ipntw ) = 0
npntw = npntw - 1
END IF
END IF
END DO
!........在波阵中查找环绕单元最小的节点作为下一主消节点.
ipoin = 0
ipnte = nelem + 1
DO ipntw = 1, npoin
jpoin = lpntw( ipntw )
IF( jpoin .GT. 0 ) THEN
IF( lpnte( jpoin ) .LT. ipnte ) THEN
ipoin = jpoin
ipnte = lpnte( jpoin )
DO idime = 1, ndime
vwork( idime ) = 0.0D0
DO jdime = 1, ndime
vwork( idime ) = vwork( idime ) - axisw( idime, jdime )*( &
basew( jdime ) - coord( jdime, ipoin ) )
END DO
END DO
ELSE IF( lpnte( jpoin ) .EQ. ipnte ) THEN
DO idime = 1, ndime
vtemp( idime ) = 0.0D0
DO jdime = 1, ndime
vtemp( idime ) = vtemp( idime ) - axisw( idime, jdime )*( &
basew( jdime ) - coord( jdime, jpoin ) )
END DO
END DO
iflag = 1
DO idime = 1, ndime
IF( iflag .EQ. 1 ) THEN
IF( vtemp( idime ) .LT. vwork( idime ) ) THEN
iflag = 0
ipoin = jpoin
DO jdime = 1, ndime
vwork( jdime ) = vtemp( jdime )
END DO
ELSE IF( vtemp( idime ) .GT. vwork( idime ) ) THEN
iflag = 0
END IF
END IF
END DO
END IF
END IF
END DO
!........如果波阵中已没有节点，重新寻找绕节点单元数最少的节点作为第一个主消节点.
IF( ipoin .EQ. 0 ) THEN
jpoin = 0
ipnte = nelem + 1
DO WHILE( jpoin .LT. npoin )
jpoin = jpoin + 1
IF( lposw( jpoin ) .EQ. 0 ) THEN
jpnte = lpnte( jpoin )
IF( jpnte .GT. 0 ) THEN
IF( jpnte .LT. ipnte ) THEN
ipoin = jpoin
ipnte = jpnte
DO idime = 1, ndime
vwork( idime ) = 0.0D0
DO jdime = 1, ndime
vwork(idime) = vwork(idime) - axisw(idime, jdime) * &
( basew(jdime) - coord(jdime, ipoin) )
END DO
END DO
ELSE IF( jpnte .EQ. ipnte ) THEN
DO idime = 1, ndime
vtemp( idime ) = 0.0D0
DO jdime = 1, ndime
vtemp( idime ) = vtemp( idime ) - axisw( idime, jdime ) * ( &
basew( jdime ) - coord( jdime, jpoin ) )
END DO
END DO
iflag = 1
DO idime = 1, ndime
IF( iflag .EQ. 1 ) THEN
IF( vtemp( idime ) .LT. vwork( idime ) ) THEN
iflag = 0
ipoin = jpoin
DO jdime = 1, ndime
vwork( jdime ) = vtemp( jdime )
END DO
ELSE IF( vtemp( idime ) .GT. vwork( idime ) ) THEN
iflag = 0
END IF
END IF
END DO
END IF
END IF
END IF
END DO
END IF
END DO
mdofw = mpntw * mdofn
nbcks = mdofw * mdofn
mstfw = mdofw * ( mdofw + 1 ) / 2
!.......打印优化后的波阵信息.
WRITE( 12, 2000 )
WRITE( 12, 2200 ) mpntw
WRITE( 12, 2400 )
WRITE( 12, 2600 ) lopts
WRITE( 12, 2800 )
WRITE( 12, 3000 ) lbcks
!.......释放临时工作数组.
DEALLOCATE( lwork )
2000 FORMAT( //20x, '波阵优化信息' )
2200 FORMAT( 10x, '波阵最大带宽', i5, '(节点)' )
2400 FORMAT( //20x, '优化单元顺序' )
2600 FORMAT( 10( 2x, i5 ) )
2800 FORMAT( //20x, '节点在波阵中的信息' / &
10x, ' 当步波宽 总体节点号 波阵位置' / )
3000 FORMAT( 10x, 2x, I6, 7x, I6, 7x, I6, 7x, I6 )
RETURN
END
SUBROUTINE SlvOpt
!........
! 模块功能
! 规范化缓冲区,总体矩阵在内存中定位,一个总体矩阵需要的外存
! 量以及记录数目.
!........
! lsavc 总体矩阵存储信息 lbcks 回代信息(I)
! lsavc(1)刚度矩阵 mdofn 最大节点自由度数(I)
! lsavc(2)质量矩阵 mpntw 最大波阵宽(I)
! lsavc(3)阻尼矩阵 nbufs 缓冲区尺度(I/O)
! (I) == 0 无此矩阵 nerrc 出错信息代码(O)
! == 1 一致矩阵 npoin 节点总数(I)
! == 2 集中矩阵 nrcrd 记录总数(O)
! (O) < 0 内存首址 nexch 数据交换区长度(记录长度)(O)
! == 1 外存存储 ndsks 外存尺度(I)
! == 2 集中矩阵 ndata 一致矩阵尺度(O)
!........
USE CtrlData
USE FrontData
IF( nerrc .GT. 0 ) RETURN
!.......缓冲区太小, 无法求解.
mexch = 0
mdofw = mpntw * mdofn
DO idofn = 1, mdofn
mexch = mexch + mdofw - idofn
END DO
IF( mexch .GT. nbufs ) THEN
WRITE( 12, 2000 )
nerrc = 1562
RETURN
END IF
!.......计算总体矩阵总数.
nmats = 0
DO isavc = 1, 3
IF( lsavc( isavc ) .EQ. 1 ) nmats = nmats + 1
END DO
!.......计算一致矩阵数据量.
ndata = 0
DO ipoin = 1, npoin
npntw = lbcks( 1, ipoin )
IF( npntw .GT. 0 ) THEN
ndofw = npntw * mdofn
DO idofn = 1, mdofn
ndata = ndata + ndofw - idofn
END DO
END IF
END DO
IF( ndata .LE. 0 ) ndata = 1
!.......判断存储方式.
nmatb = nbufs / ndata
IF( nmatb .GE. nmats ) THEN
!........内存足以存储所有的一致矩阵.
nexch = 0
nrcrd = 1
kdata = ndata
ELSE IF( nmatb .GE. 1 ) THEN
!........内存足以存储一个一致矩阵.
nrcrd = 1
nexch = ndata
kdata = ndata
nmatb = nmatb - 1
ELSE
!........内存不足以存储一个一致矩阵.
nrcrd = 0
kdata = 0
nexch = nbufs
DO ipoin = npoin, 1,-1
npntw = lbcks( 1, ipoin )
ndofw = npntw * mdofn
idata = ndofw * mdofn - ( mdofn + 1 ) * mdofn / 2
IF( kdata + idata .GT. nexch ) THEN
nrcrd = nrcrd + 1
kdata = idata
ELSE
kdata = kdata + idata
END IF
END DO
IF( kdata .GT. 0 ) nrcrd = nrcrd + 1
kdata = nrcrd * nexch
END IF
!.......外存可存一致矩阵数目.
nmatd = ndsks / kdata
IF( nmatd + nmatb .LT. nmats ) THEN
!........若内外存存储空间不够,给出出错信息.
WRITE( 12, 2200 ) nmats * kdata
nerrc = 1662
RETURN
END IF
!.......确定内存位置.
nuseb = nexch + 1
DO isavc = 1, 3
IF( lsavc( isavc ) .EQ. 1 ) THEN
IF( nmatb .GT. 0 ) THEN
lsavc( isavc ) =-nuseb
nuseb = nuseb + ndata
nmatb = nmatb - 1
END IF
END IF
END DO
WRITE( 12, 2400 )
IF( lsavc( 1 ) .LT. 0 ) WRITE( 12, 2600 ) -lsavc( 1 )
IF( lsavc( 1 ) .EQ. 1 ) WRITE( 12, 2800 )
IF( lsavc( 1 ) .EQ. 2 ) WRITE( 12, 3000 )
IF( lsavc( 2 ) .LT. 0 ) WRITE( 12, 3200 ) -lsavc( 2 )
IF( lsavc( 2 ) .EQ. 1 ) WRITE( 12, 3400 )
IF( lsavc( 2 ) .EQ. 2 ) WRITE( 12, 3600 )
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -