This example generates a finite element mesh using 8-node
bricks. The mesh is NE_THICK by NE_WIDE by NE_LONG.
FMSlib is then used to solve the problem for a specified
unit force on one end.
C E X A M P L E 2 4
C
C Program name:
CHARACTER*10 MYNAME
PARAMETER (MYNAME='EXAMPLE_24')
C
C Input data functions:
LOGICAL ASK
INTEGER ASK_I
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 Number of element vectors:
PARAMETER (NUMELV = 0)
C
C Skip forward reduction during call to RSDS:
PARAMETER (ISKIP = 0)
C
C FMS matrix and vector file attributes:
INTEGER LUA(25)
INTEGER LUX(25), LUB(25), LUAX(25)
INTEGER LUS(25)
C
C FMS memory management requires the following arrays:
POINTER (IMD_PTR, IMD)
POINTER (RMD_PTR, RMD)
POINTER (CMD_PTR, CMD)
INTEGER IMD(0:1)
REAL*8 RMD(0:1)
COMPLEX*16 CMD(0:1)
C
C Allocated array starting values:
INTEGER NEWNUM
INTEGER LB
INTEGER LX
INTEGER IREC
INTEGER LS
INTEGER LOWEQ
INTEGER LAX
INTEGER LIFIX
C
C FMS Parameter values;
INTEGER MPOSDF
DATA MPOSDF/1/
INTEGER IPRI
DATA IPRI/1/
INTEGER IPRF
DATA IPRF/1026/
INTEGER MDATAU
DATA MDATAU/2/
INTEGER LOGTIM
DATA LOGTIM/3/
INTEGER INSIDE
DATA INSIDE/-1/
INTEGER NEQBLK
DATA NEQBLK/96/
INTEGER INCORE
DATA INCORE/2/
C
C Big diagonal value to impose constraints:
REAL*8 BIGNUM
DATA BIGNUM/1.0D30/
C
C Read in problem data:
INTEGER NETHICK
INTEGER NEWIDE
INTEGER NELONG
INTEGER NDOF
C
C Local variables:
INTEGER MAXEQ, NELDOF, NUMEQ
INTEGER LENR, LENI, LENV
INTEGER INEW, I
INTEGER LENX, LDISK, NV
REAL*8 EI, ERROR
C
C Dummy arguments:
INTEGER LUDUMMY(25)
DATA LUDUMMY(1)/0/
REAL*8 ALPHA(1)
DATA ALPHA(1)/0.0D0/
C
C Common block to communicate with user supplied subroutine.
COMMON/MYDATA/NPTHICK, NPWIDE, NPLONG, NDOF,
1 BIGNUM, LIFIX
INTEGER NPTHICK, NPWIDE, NPLONG
C
C (1) Initialize FMS:
CALL FMSINI
CALL FMSPSH(MYNAME)
C
C Require matrix to be positive definite:
CALL FMSIST ('MPOSDF',MPOSDF)
C
C Call RSUBLK during factoring:
CALL FMSIST ('MDATAU',MDATAU)
C
C Get the memory pointers:
C
C Turn on timings:
CALL FMSIGT ('LOGTIM', LOGTIM)
IF(LOGTIM .LT. 3) CALL FMSIST ('LOGTIM', 3)
C
C Set matrix initialization print level:
CALL FMSIST ('IPRI' , IPRI)
C
C Select outside blocking of matrix profile:
CALL FMSIST ('INSIDE' , INSIDE)
C
C Select block size:
CALL FMSIST ('NEQBLK' , NEQBLK)
C
C Set matrix factoring print level:
CALL FMSIGT ('IPRF', IPRF)
IF(IAND(IPRF, 2) .EQ. 0) IPRF = IPRF + 2
IF(IAND(IPRF,1024) .EQ. 0) IPRF = IPRF + 1024
CALL FMSIST ('IPRF' , IPRF)
C
C Store data in memory if it will fit:
CALL FMSIST ('INCORE' , INCORE)
C
C Read in problem data:
1 CONTINUE
NETHICK = ASK_I('Enter the number of elements thick')
NEWIDE = ASK_I('Enter the number of elements wide')
NELONG = ASK_I('Enter the number of elements long')
NUMEL = NETHICK * NEWIDE * NELONG
IF(NUMEL .LE. 0) GO TO 1
2 CONTINUE
NDOF = ASK_I('Enter the number of degrees of freedom per node')
IF(NDOF .LE. 0) GO TO 2
WRITE(6,*) 'You may now alter any FMS parameter.'
WRITE(6,*) 'When you are finished, type the letters RETURN'
CALL FMSSET
CALL FMSIGT ('MEMPTR',IMD_PTR)
CALL FMSIGT ('MEMPTR',RMD_PTR)
CALL FMSIGT ('MEMPTR',CMD_PTR)
NPTHICK = NETHICK + 1
NPWIDE = NEWIDE + 1
NPLONG = NELONG + 1
MAXEQ = NPTHICK * NPWIDE * NPLONG * NDOF
NELDOF = 8 * NDOF
C
C (2) Open FMS files:
C
C Submatrix File:
LENR = (NELDOF * (NELDOF+1))/2
LENI = 2 + NELDOF
LENV = NELDOF
CALL FMSOS (LENR, LENI, LENV, NUMELV, NUMEL, 'LUS', LUS)
C
C Allocate storage for element integer record:
LENI = LUS(4)
CALL FMSIMG (IMD, IREC, LENI)
C
C Allocate storage for element real record:
LENR = LUS(3)
CALL FMSRMG (RMD, LS, LENR)
C
C Generate element data:
CALL ELGEN (NETHICK, NEWIDE, NELONG, NUMEL,
1 LUS, IMD(IREC), RMD(LS), NELDOF)
WRITE(6,*) 'NUMBER OF ELEMENTS = ', NUMEL
WRITE(6,*) 'INITIAL NUMBER OF EQUATIONS = ', MAXEQ
C
C Free storage for element real record:
CALL FMSRMR (RMD, LS, LENR)
C
C Free storage for element integer record:
CALL FMSIMR (IMD, IREC, LENI)
C
C Allocate storage for NEWNUM vector:
CALL FMSIMG (IMD, NEWNUM, MAXEQ)
C
C Allocate storage for solution vector {X}:
CALL FMSRMG (RMD, LX, MAXEQ)
C
C Impose displacement boundary conditions:
CALL XBC (RMD(LX), IMD(NEWNUM), MAXEQ)
C
C Allocate storage for profile vector:
CALL FMSIMG (IMD, LOWEQ, MAXEQ)
C
C Compute profile vector; eliminate zero constraints:
CALL FMSPF (LUS, NUMSF, NEWNUM, MAXEQ, LOWEQ, NUMEQ)
WRITE(6,*) 'ACTUAL NUMBER OF EQUATIONS = ', NUMEQ
C
C Open the global matrix file:
CALL RSDI (LOWEQ, NUMEQ, 'LUA', LUA)
C
C Free storage for profile vector:
CALL FMSIMR (IMD, LOWEQ, MAXEQ)
C
C (3) Write data to FMS files:
C
C Allocate storage for IFIX vector to modify diagonals for
C nonzero constrained displacements:
CALL FMSIMG (IMD, LIFIX, NUMEQ)
C
C Flag diagonals with nonzero specified displacements:
DO I = 1,MAXEQ
INEW = IMD(NEWNUM+I-1)
IF(INEW .LT. 0) THEN
C Value constrained to a nonzero value:
INEW = -INEW
IMD(LIFIX+INEW-1) = 1
ELSE
IMD(LIFIX+INEW-1) = 0
END IF
END DO
C
C Allocate storage for RHS vector {B}:
CALL FMSRMG (RMD, LB, MAXEQ)
C
C Open an incore file for {X}:
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, RMD(LX), NUMEQ, LUX)
C
C Open an incore file for {B}:
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, RMD(LB), NUMEQ, LUB)
C
C (4) Perform matrix algebra:
C Build and factor the matrix:
CALL RSDAF (LUDUMMY, ALPHA, 0, LUS, NUMSF, LUDUMMY,
1 LUA, LUB, LUX, 0)
C
C Build {B} and solve for {X}:
CALL BBC(RMD(LB), IMD(NEWNUM), MAXEQ)
C
C Solve for {X}:
CALL RSDS (LUA, LUB, LUX, NUMRHS, 0)
C
C (5) Read data from FMS files:
C Not required.
C
C Compute [A]{X} from SUM([Si]{X}) to check error:
C
C Allocate storage for [A]{X} vector:
CALL FMSRMG (RMD, LAX, MAXEQ)
C
C Open an incore file for [A]{X} vector:
CALL FMSOV2 (NUMEQ, IDTYPE, NUMRHS, RMD(LAX), NUMEQ, LUAX)
C
C Compute {AX} = [Si]{X}:
CALL RSDSVM (LUS, NUMSF, LUX, LUAX, NUMRHS)
C
C Compute the error {AX} - {B}:
ERROR = 0.0D0
C Expand solution vector to include zero constrained values:
DO 50 I = MAXEQ,1,-1
INEW = IMD(NEWNUM+I-1)
IF(INEW .EQ. 0) THEN
C Constrained zero value:
RMD(LB+I-1) = 0.0D0
ELSE IF(INEW .LT. 0) THEN
C Constrained nonzero value:
INEW = -INEW
RMD(LB+I-1) = RMD(LB+INEW-1)
ELSE
C Value which was determined:
RMD(LB+1-I) = RMD(LB+INEW-1)
EI = ABS(RMD(LAX+INEW-1)-RMD(LB+INEW-1))
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 FMSRMR (RMD, LAX, MAXEQ)
CALL FMSCV (LUX)
CALL FMSRMR (RMD, LX, MAXEQ)
CALL FMSCV (LUB)
CALL FMSRMR (RMD, LB, MAXEQ)
CALL FMSIMR (IMD, LIFIX, NUMEQ)
CALL FMSCS (LUS)
CALL FMSIMR (IMD, NEWNUM, MAXEQ)
IF(ASK('Do you want another solution?')) GO TO 1
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(*)
POINTER (IMD_PTR, IMD)
INTEGER IMD(0:1)
COMMON/MYDATA/NPTHICK, NPWIDE, NPLONG, NDOF,
1 BIGNUM, LIFIX
INTEGER LIFIX
REAL*8 BIGNUM
C
CALL FMSIGT ('MEMPTR', IMD_PTR)
C Modify diagonal for nonzero constraint equations:
IF(IROW2 .EQ. JCOL2) THEN
C This is a diagonal block:
DO I = IROW1,IROW2
IF(IMD(LIFIX-1+I) .NE. 0) D(I) = BIGNUM
END DO
END IF
RETURN
END
C=======================================================================
SUBROUTINE ELGEN (
1 NETHICK, NEWIDE, NELONG, NUMEL,
2 LUS, IREC, S, NELDOF)
C=======================================================================
INTEGER NETHICK, NEWIDE, NELONG, NUMEL
INTEGER LUS(25)
INTEGER IREC(*)
INTEGER NELDOF
REAL*8 S(*)
COMMON/MYDATA/NPTHICK, NPWIDE, NPLONG, NDOF,
1 BIGNUM, LIFIX
INTEGER NPTHICK, NPWIDE, NPLONG, NDOF
INTEGER LENI, LENR, LENV
REAL*8 DIAG
CHARACTER*5 MYNAME
PARAMETER (MYNAME='ELGEN')
C
CALL FMSPSH(MYNAME)
LENR = LUS(3)
LENI = LUS(4)
LENV = LUS(10)
DIAG = DFLOAT(NELDOF-1)
D WRITE(LUDBUG,2000),
D 1 NETHICK, NEWIDE, NELONG, NUMEL,
D 2 LOC(LUS), LOC(IREC), LOC(S), NELDOF
C
C Generate an 8-node finite element mesh that is NETHICK elements
C thick, NEWIDE elements wide by NELONG elements long.
C Number the nodes first in the thick direction, then in the
C wide direction and then in the long direction.
C Use the same matrix S(NELDOF,NELDOF) for all elements.
C
C ========================
C Element integer records:
C ========================
C Write the element integer file.
C Each record contains the following:
C IREC(1) = Number of element degrees of freedom.
C IREC(2) = Element format type (1:5)
C IREC(3:LENI) = Element global equation numbers
CALL FMSEEK (LUS(2), 0)
IREC(1) = NELDOF
IREC(2) = 4
C DOF increments:
INCTHICK = NDOF
INCWIDE = INCTHICK * NPTHICK
INCLONG = INCWIDE * NPWIDE
D WRITE(LUDBUG,2005) INCTHICK, INCWIDE, INCLONG
C Element numbering:
NP1 = 2
NP2 = NP1 + NDOF
NP3 = NP2 + NDOF
NP4 = NP3 + NDOF
NP5 = NP4 + NDOF
NP6 = NP5 + NDOF
NP7 = NP6 + NDOF
NP8 = NP7 + NDOF
D WRITE(LUDBUG,2001) NP1, NP2, NP3, NP4, NP5, NP6, NP7, NP8
C Global numbering:
NP1G = 0
NP2G = NP1G + NDOF
NP3G = NP1G + INCWIDE
NP4G = NP3G + NDOF
NP5G = NP1G + INCLONG
NP6G = NP5G + NDOF
NP7G = NP5G + INCWIDE
NP8G = NP7G + NDOF
D WRITE(LUDBUG,2002) NP1G, NP2G, NP3G, NP4G, NP5G, NP6G, NP7G,
D 1 NP8G
NDOFLONG1 = 0
NE = 0
D WRITE(LUDBUG,2003)
DO I = 1,NELONG
NDOFWIDE1 = NDOFLONG1
DO J = 1,NEWIDE
NDOFTHICK1 = NDOFWIDE1
DO K = 1,NETHICK
NDOF1 = NDOFTHICK1
NE = NE + 1
DO IDOF = 1,NDOF
NDOF1 = NDOF1 + 1
IREC(NP1+IDOF) = NP1G + NDOF1
IREC(NP2+IDOF) = NP2G + NDOF1
IREC(NP3+IDOF) = NP3G + NDOF1
IREC(NP4+IDOF) = NP4G + NDOF1
IREC(NP5+IDOF) = NP5G + NDOF1
IREC(NP6+IDOF) = NP6G + NDOF1
IREC(NP7+IDOF) = NP7G + NDOF1
IREC(NP8+IDOF) = NP8G + NDOF1
END DO
D WRITE(LUDBUG,2004) NE, (IREC(II),II=1,LENI)
CALL FMSSWR (LUS(2), IREC, LENI)
NDOFTHICK1 = NDOFTHICK1 + INCTHICK
END DO
NDOFWIDE1 = NDOFWIDE1 + INCWIDE
END DO
NDOFLONG1 = NDOFLONG1 + INCLONG
END DO
C =====================
C Element real records:
C =====================
C Use lower triangle stored by row format:
L = 0
DO I = 1,NELDOF
DO J = 1,(I-1)
L = L + 1
S(L) = -1.0D0
END DO
L = L + 1
S(L) = DIAG
END DO
C
C Write element real records:
CALL FMSEEK (LUS(1), 0)
DO NE = 1,NUMEL
CALL FMSSWR (LUS(1), S, LENR)
END DO
CALL FMSPOP(MYNAME)
RETURN
D2000 FORMAT(
D 1 'SUBROUTINE ELGEN('/
D 2 ' NETHICK=',I5,', NEWIDE=',I5,', NELONG=',I5,', NUMEL=',I5/
D 5 ' &LUS =',I15,', &IREC=',I15,', &S=',I15,' NELDOF=',I5,'0'/)
D2001 FORMAT(' NP1 -NP8 :',I4,7I5)
D2002 FORMAT(' NP1G-NP8G:',I4,7I5)
D2003 FORMAT(/
D 1 ' NE SIZE TYPE EQ1'/
D 2 '===== ==== ==== ====')
D2004 FORMAT(15I5/(15X,12I5))
D2005 FORMAT(' INCTHICK, INCWIDE, INCLONG: ',3I5)
END
C=======================================================================
SUBROUTINE XBC (X, NEWNUM, MAXEQ)
C=======================================================================
C Impose displacement boundary conditions:
REAL*8 X(MAXEQ)
INTEGER NEWNUM(MAXEQ)
CHARACTER*3 MYNAME
COMMON/MYDATA/NPTHICK, NPWIDE, NPLONG, NDOF,
1 BIGNUM, LIFIX
INTEGER NPTHICK, NPWIDE, NPLONG
PARAMETER (MYNAME='XBC')
CALL FMSPSH(MYNAME)
C
C Initialize all values to free:
DO I = 1,MAXEQ
NEWNUM(I) = 1
X(I) = 0.0D0
END DO
C
C Constrain all nodes at Z=0 to z=0.
IDOF = 0
DO I=1,NPTHICK
DO J=1,NPWIDE
IDOF = IDOF + NDOF
NEWNUM(IDOF) = 0
END DO
END DO
C
C Constrain node 1 to (0,0,0):
DO I=1,NDOF
NEWNUM(I) = 0
END DO
C
C Constrain node at (1,NEWIDE,1) to x=0:
IDOF = NDOF*(NPTHICK * (NPWIDE-1)) + 1
NEWNUM(IDOF) = 0
D WRITE(6,2000)
D DO I=1,MAXEQ
D WRITE(6,2001) I, NEWNUM(I), X(I)
D END DO
CALL FMSPOP(MYNAME)
RETURN
D2000 FORMAT(/
D 1 ' I NEWNUM X'/
D 2 '===== ====== ======')
D2001 FORMAT(I5, 2I9)
END
C=======================================================================
SUBROUTINE BBC (B, NEWNUM, MAXEQ)
C=======================================================================
C Impose force boundary conditions:
REAL*8 B(MAXEQ)
INTEGER NEWNUM(MAXEQ)
COMMON/MYDATA/NPTHICK, NPWIDE, NPLONG, NDOF,
1 BIGNUM, LIFIX
INTEGER NPTHICK, NPWIDE, NPLONG
CHARACTER*3 MYNAME
PARAMETER (MYNAME='BBC')
CALL FMSPSH(MYNAME)
C
C Apply a force of 1 on all Z nodes at Z=NPLONG:
NFDOF = NDOF * NPTHICK * NPWIDE
DO I = (MAXEQ-NFDOF+NDOF),MAXEQ,NDOF
IF(NEWNUM(I) .LE. 0) THEN
C DOF already has a specified displacement:
PRINT *,'BBC: DOF',I, 'has conflicted boundary conditions'
STOP 'FATAL ERROR IN SUBROUTINE BBC'
END IF
B(I) = 1.0D0
END DO
C
C Reorder {B} for new numbering system:
DO I = 1,MAXEQ
INEW = NEWNUM(I)
IF(INEW .GT. 0) THEN
C Value to be determined:
B(INEW) = B(I)
ELSE IF(INEQ .EQ. 0) THEN
C Value constrained to 0:
B(INEW) = 0.0D0
ELSE
C Value constrained to a nonzero value:
INEW = -INEW
B(INEW) = BIGNUM*B(I)
END IF
END DO
CALL FMSPOP(MYNAME)
RETURN
END
C=======================================================================
LOGICAL FUNCTION ASK(QUESTION)
C=======================================================================
CHARACTER* (*) QUESTION
CHARACTER*1 IYN
WRITE(6,2000) QUESTION
READ (5,1000) IYN
IF( (IYN .EQ. 'Y') .OR. (IYN .EQ. 'y') ) THEN
ASK = .TRUE.
ELSE
ASK = .FALSE.
END IF
RETURN
1000 FORMAT (A)
2000 FORMAT (1X,A,' (y,n)>')
END
C=======================================================================
INTEGER FUNCTION ASK_I(STRING)
C=======================================================================
CHARACTER* (*) STRING
WRITE(6,2000) STRING
READ (5,*) ASK_I
RETURN
2000 FORMAT (1X,A,'>')
END