C E X A M P L E 4
C
C Program name:
CHARACTER*9 MYNAME
PARAMETER (MYNAME='EXAMPLE_4')
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, NUMEQ = 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 file attributes:
INTEGER LUA(25)
INTEGER LUX(25), LUB(25), LUAX(25)
INTEGER LUS(25)
C
C Profile vector and equation numbers:
INTEGER LOWEQ(NUMEQ), NEWNUM(NUMEQ)
INTEGER IDUM(1)
C
C Solution vector {X}
C RHS vector {B}
C Product [A]{X}, {AX}
C Diagonals {D}
REAL*8 X(NUMEQ), B(NUMEQ), AX(NUMEQ), D(NUMEQ)
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, NUMEQ)
WRITE(6,*) 'NUMBER OF ELEMENTS = ', NUMEL
WRITE(6,*) 'NUMBER OF EQUATIONS = ', NUMEQ
C Compute profile vector:
CALL FMSPF (LUS, NUMSF, IDUM, 0, LOWEQ, 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:
CALL FMSRED (LUA(2), 1, D, NUMEQ)
DO 31 I = 1,NUMEQ
INEW = NEWNUM(I)
IF(INEW .LT. 0) THEN
B(I) = BIGNUM*B(I)
D(I) = 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
DO 50 I = 1,NUMEQ
INEW = NEWNUM(I)
IF(INEW .GT. 0) THEN
C Value which was determined:
EI = ABS(AX(I)-B(I))
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