This example illustrates how subroutine FMSPF may be used to
eliminate equations whose solution value is constrained to
zero and to modify the profile of nonzero constrained
equations to a diagonal only term. The vector NEWNUM is used
to relate new equation numbers in the condensed system to old
equation numbers of the original system.
C E X A M P L E 5
C
C Program name:
CHARACTER*9 MYNAME
PARAMETER (MYNAME='EXAMPLE_5')
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), LUB(25), LUAX(25)
INTEGER LUS(25)
C
C Profile vector and equation numbers:
INTEGER LOWEQ(MAXEQ), NEWNUM(MAXEQ)
C
C Solution vector {X}
C RHS vector {B}
C Product [A]{X}, {AX}
C Diagonals {D}
REAL*8 X(MAXEQ), B(MAXEQ), AX(MAXEQ), D(MAXEQ)
COMMON X, B, AX, D
C
C Big diagonal value to impose constraints:
REAL*8 BIGNUM
PARAMETER (BIGNUM = 1.0D30)
C
C Local variables:
REAL*8 EI, ERROR
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSPSH(MYNAME)
CALL FMSIST ('MPOSDF',1)
CALL FMSIST ('INCORE',2)
C
C (2) Open FMS files:
CALL FMSOS (LENR, LENI, LENV, NUMELV, NUMEL, 'LUS', LUS)
C Generate test data:
CALL ELGEN (NEWIDE, NELONG, NUMEL, NEWNUM, LUS, B, MAXEQ)
WRITE(6,*) 'NUMBER OF ELEMENTS = ', NUMEL
WRITE(6,*) 'INITIAL NUMBER OF EQUATIONS = ', MAXEQ
C Compute profile vector; eliminate zero constraints:
CALL FMSPF (LUS, NUMSF, NEWNUM, MAXEQ, LOWEQ, NUMEQ)
WRITE(6,*) 'ACTUAL NUMBER OF EQUATIONS = ', NUMEQ
CALL RSDI (LOWEQ, NUMEQ, 'LUA', LUA)
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, X, NUMEQ, LUX)
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, B, NUMEQ, LUB)
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, AX, NUMEQ, LUAX)
C
C (3) Write data to FMS files:
CALL RSDA (LUS, NUMSF, LUA)
C Impose boundary conditions on diagonal and RHS vector:
C Reorder RHS terms for new numbering system:
CALL FMSRED (LUA(2), 1, D, NUMEQ)
DO 31 I = 1,MAXEQ
INEW = NEWNUM(I)
IF(INEW .EQ. 0) GO TO 31
IF(INEW .GT. 0) THEN
B(INEW) = B(I)
ELSE
INEW = -INEW
B(INEW) = BIGNUM*B(I)
D(INEW) = BIGNUM
END IF
31 CONTINUE
CALL FMSWRT (LUA(2), 1, D, NUMEQ)
C
C (4) Perform matrix algebra:
CALL RSDF (LUA, LUA, LUB, LUX, NUMRHS)
CALL RSDS (LUA, LUB, 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 RSDSVM (LUS, NUMSF, LUX, LUAX, NUMRHS)
ERROR = 0.0D0
C Expand solution vector to include zero constrained values:
DO 50 I = MAXEQ,1,-1
INEW = NEWNUM(I)
IF(INEW .EQ. 0) THEN
C Constrained zero value:
X(I) = 0.0D0
ELSE IF(INEW .LT. 0) THEN
C Constrained nonzero value:
INEW = -INEW
X(I) = X(INEW)
ELSE
C Value which was determined:
X(I) = X(INEW)
EI = ABS(AX(INEW)-B(INEW))
IF(EI .GT. ERROR) ERROR = EI
END IF
50 CONTINUE
WRITE(6,*) 'MAXIMUM ERROR =',ERROR
C
C (6) Close FMS files:
CALL FMSCM (LUA)
CALL FMSCV (LUAX)
CALL FMSCV (LUX)
CALL FMSCV (LUB)
CALL FMSCS (LUS)
CALL FMSPOP(MYNAME)
CALL FMSEND
END
C=======================================================================
SUBROUTINE ELGEN (NEWIDE,NELONG,NUMEL,NEWNUM,LUS,B,NUMEQ)
C=======================================================================
CHARACTER*5 MYNAME
PARAMETER (MYNAME='ELGEN')
INTEGER NEWIDE, NELONG, NUMEL, NUMEQ
INTEGER LUS(25), NEWNUM(NUMEQ)
INTEGER IREC(6), NSUB, ITYPE, ID(4)
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