This example illustrates how
FMS is used in
finite element programs which use a wavefront solution
algorithm. First you call FMSOS once to create a single
submatrix file on which you store your element data in the
wavefront order. The order that you write your records to
this file will determine the wavefront numbering scheme.
Then you call subroutine FMSWF to generate global equation
numbers and update the submatrix integer file. Subroutine
FMSWF also computes the profile vector LOWEQ and the NEWNUM
vector which relates the equation labels you wrote to the
submatrix integer file to the new global equation numbers.
This example also illustrates how you can write a user
supplied subroutine to modify the equations after assembly
before the factoring process begins. In this example, RSUBLK
is used to modify diagonal terms to impose specified solution
values.
C E X A M P L E 6
C
C Program name:
CHARACTER*9 MYNAME
PARAMETER (MYNAME='EXAMPLE_6')
C
C Problem size parameters:
C
C Number of elements wide and long:
PARAMETER (NEWIDE = 20, NELONG = 100)
C
C Number of nodes wide and long:
PARAMETER (NPWIDE = NEWIDE+1, NPLONG = NELONG+1)
C
C Number of elements and equations:
PARAMETER (NUMEL = NEWIDE*NELONG, MAXEQ = NPWIDE*NPLONG)
C
C Data type = real:
PARAMETER (IDTYPE = 1)
C
C Number of RHS vectors:
PARAMETER (NUMRHS = 1)
C
C Number of element files:
PARAMETER (NUMSF = 1)
C
C Length of element real record:
PARAMETER (LENR = 4*4)
C
C Length of element integer record:
PARAMETER (LENI = 2 + 4)
C
C Length of element vector record:
PARAMETER (LENV = 2)
C
C Number of element vectors:
PARAMETER (NUMELV = 0)
C
C Skip forward reduction during call to RSDS:
C (Performed during call to RSDF)
PARAMETER (ISKIP = 1)
C
C FMS matrix and vector file attributes:
INTEGER LUA(25)
INTEGER LUX(25), LUAX(25)
INTEGER LUS(25)
C
C Solution vector {X}
C RHS vector {B}
C Product [A]{X}, {AX}
REAL*8 RDUM(1), X(MAXEQ), B(MAXEQ), AX(MAXEQ)
COMMON /DATA/ RDUM, X, B, AX
C
C Big diagonal value to impose constraints:
REAL*8 BIGNUM
C
C Local variables:
REAL*8 EI, ERROR
C
C Common block to communicate with RSUBLK:
COMMON /EXAMP6/ BIGNUM, IBC
C
C FMS memory management requires the following arrays:
POINTER (CMD_PTR, CMD)
POINTER (RMD_PTR, RMD)
POINTER (IMD_PTR, IMD)
COMPLEX*16 CMD(0:1)
REAL*8 RMD(0:1)
INTEGER IMD(0:1)
DATA LUA0(1)/0/
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSPSH(MYNAME)
CALL FMSIST ('MPOSDF',1)
CALL FMSIST ('INCORE',2)
CALL FMSIST ('MDATAU',2)
CALL FMSIGT ('MEMPTR', IMD_PTR)
CALL FMSIGT ('MEMPTR', RMD_PTR)
CALL FMSIGT ('MEMPTR', CMD_PTR)
WRITE(6,*) 'NUMBER OF ELEMENTS = ', NUMEL
WRITE(6,*) 'INITIAL NUMBER OF EQUATIONS = ', MAXEQ
BIGNUM = 1.0D30
C
C (2) Open FMS files:
CALL FMSOS (LENR, LENI, LENV, NUMELV, NUMEL, 'LUS', LUS)
C Generate element data:
CALL FMSIMG (IMD, NEWNUM, MAXEQ)
CALL ELGEN (NEWIDE,NELONG,NUMEL,IMD(NEWNUM),LUS,B,MAXEQ)
C Compute profile vector; eliminate zero constraints:
CALL FMSIMG (IMD, LOWEQ, MAXEQ)
CALL FMSWF (LUS, IMD(NEWNUM), MAXEQ, IMD(LOWEQ), NUMEQ)
WRITE(6,*) 'ACTUAL NUMBER OF EQUATIONS = ', NUMEQ
CALL RSDI (IMD(LOWEQ), NUMEQ, 'LUA', LUA)
CALL FMSIMR (IMD, LOWEQ, MAXEQ)
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, X, NUMEQ, LUX)
C
C (3) Write data to FMS files:
C Elements written during ELGEN.
C Assembly during factoring.
C Impose boundary conditions on RHS vector:
C Flag IBC array for constrained diagonals:
C Reorder RHS terms for new numbering system:
CALL FMSIMG (IMD, IBC, NUMEQ)
DO 31 I = 1,MAXEQ
INEW = IMD(NEWNUM-1+I)
IF(INEW .EQ. 0) GO TO 31
IF(INEW .GT. 0) THEN
IMD(IBC-1+INEW) = 0
X(INEW) = B(I)
ELSE
INEW = -INEW
IMD(IBC-1+INEW) = 1
X(INEW) = BIGNUM*B(I)
END IF
31 CONTINUE
C
C (4) Perform matrix algebra:
CALL RSDAF (LUA0, RDUM, 0, LUS, NUMSF, LUA0, LUA,
1 LUX, LUX, NUMRHS)
CALL FMSIMR (IMD, IBC, NUMEQ)
CALL RSDS (LUA, LUX, LUX, NUMRHS, ISKIP)
C
C (5) Read data from FMS files:
C Not required.
C
C Compute [A]{X} from SUM([Si]{X}) to check error:
CALL FMSRMG (RMD, LAX, NUMEQ)
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, RMD(LAX), NUMEQ, LUAX)
CALL RSDSVM (LUS, NUMSF, LUX, LUAX, NUMRHS)
ERROR = 0.0D0
C Expand solution vector to include zero constrained values:
DO 50 I = 1,MAXEQ
INEW = IMD(NEWNUM-1+I)
IF(INEW .EQ. 0) THEN
C Constrained zero value:
B(I) = 0.0D0
ELSE IF(INEW .LT. 0) THEN
C Constrained nonzero value:
INEW = -INEW
B(I) = X(INEW)
ELSE
C Value which was determined:
EI = ABS(RMD(LAX-1+INEW) - B(I))
IF(EI .GT. ERROR) ERROR = EI
B(I) = X(INEW)
END IF
50 CONTINUE
WRITE(6,*) 'MAXIMUM ERROR =',ERROR
CALL FMSIMR (IMD, NEWNUM, MAXEQ)
C
C (6) Close FMS files:
CALL FMSCV (LUAX)
CALL FMSRMR (RMD, LAX, NUMEQ)
CALL FMSCM (LUA)
CALL FMSCV (LUX)
CALL FMSCS (LUS)
CALL FMSPOP(MYNAME)
CALL FMSEND
END
C=======================================================================
SUBROUTINE RSUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2,
1 JCOL1, JCOL2, IJSTEP)
C=======================================================================
INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP
INTEGER LOWEQ(*), LOCEQ(*)
REAL*8 A(0:*), D(*), BIGNUM
COMMON /EXAMP6/ BIGNUM, IBC
C
POINTER (CMD_PTR, CMD)
POINTER (RMD_PTR, RMD)
POINTER (IMD_PTR, IMD)
COMPLEX*16 CMD(0:1)
REAL*8 RMD(0:1)
INTEGER IMD(0:1)
C
CALL FMSIGT ('MEMPTR', IMD_PTR)
CALL FMSIGT ('MEMPTR', RMD_PTR)
CALL FMSIGT ('MEMPTR', CMD_PTR)
IF(IROW2 .EQ. JCOL2) THEN
C This is a diagonal block:
DO 10 IROW = IROW1,IROW2
IF(IMD(IBC-1+IROW) .NE. 0) D(IROW) = BIGNUM
10 CONTINUE
END IF
RETURN
END
C=======================================================================
SUBROUTINE ELGEN (NEWIDE,NELONG,NUMEL,NEWNUM,LUS,B,NUMEQ)
C=======================================================================
INTEGER NEWIDE, NELONG, NUMEL, NUMEQ
INTEGER LUS(25), NEWNUM(NUMEQ)
INTEGER IREC(6), NSUB, ITYPE, ID(4)
CHARACTER*5 MYNAME
PARAMETER (MYNAME='ELGEN')
REAL*8 S(4,4), B(NUMEQ)
EQUIVALENCE (IREC(1),NSUB), (IREC(2),ITYPE), (IREC(3),ID)
C
CALL FMSPSH(MYNAME)
C
C Generate a 4-node quadrilateral finite element mesh that
C is NEWIDE elements wide by NELONG elements long. Use the
C same matrix S(4,4) for all elements. Impose the following
C solution values: X=0 along the top, X=1 along the bottom.
C
C Generate element connectivity array ID(4):
C Write element integer records:
CALL FMSEEK (LUS(2), 0)
NSUB = 4
ITYPE = 1
NPWIDE = NEWIDE + 1
DO 20 I = 1,NELONG
DO 10 J = 1,NEWIDE
ID(1) = I*NPWIDE + J
ID(2) = ID(1) + 1
ID(3) = ID(2) - NPWIDE
ID(4) = ID(1) - NPWIDE
CALL FMSSWR (LUS(2), IREC, 6)
10 CONTINUE
20 CONTINUE
C
C Generate element data array S(4,4):
DO 40 I = 1,4
DO 30 J = 1,4
S(I,J) = -1.0D0
30 CONTINUE
S(I,I) = 3.0D0
40 CONTINUE
C
C Write element real records:
CALL FMSEEK (LUS(1), 0)
DO 50 NE = 1,NUMEL
CALL FMSSWR (LUS(1), S, 16)
50 CONTINUE
C
C Generate right-hand side vector B(NUMEQ)
C Generate boundary condition code vector NEWNUM(NUMEQ)
DO 60 I = 1,NUMEQ
B(I) = 0.0D0
NEWNUM(I) = 1
60 CONTINUE
C Specify X(1:NPWIDE) = 0, X(NUMEQ-NPWIDE+1:NUMEQ) = 1
I2 = NUMEQ - NPWIDE
DO 70 I1 = 1,NPWIDE
NEWNUM(I1) = 0
B(I1) = 0.0D0
I2 = I2 + 1
NEWNUM(I2) = -1
B(I2) = 1.0D0
70 CONTINUE
CALL FMSPOP(MYNAME)
RETURN
END