PROGRAM DRIVER C * * * * * * * * * * * * * * * * * * * * * * * * * C DUMMY MAIN TO SET PROGRAM CAPACITY LIMITS: C C TO INCREASE MEMORY STORAGE SIZE CHANGE PARAMETER C MAXR FOR REAL ARRAYS AND/OR C MAXI FOR INTEGER ARRAYS C C TO INCREASE THE NUMBER OF ARRAY NAMES & SIZES CHANGE C NUMR FOR REAL ARRAYS AND/OR C NUMI FOR INTEGER ARRAYS C C ERROR MESSAGES WILL STATE WHEN THESE ARE NECESSARY C * * * * * * * * * * * * * * * * * * * * * * * * * CHARACTER*8 RN, IN PARAMETER (NUMR= 55, MAXR= 55000) PARAMETER (NUMI= 25, MAXI= 5000) DIMENSION R(MAXR), RN(NUMR), I(MAXI), IN(NUMI) DIMENSION J(NUMR), K(NUMI) C C USER ACCESS VIA COMMON IS ALLOWED. STANDARD ACCESS IS C ALWAYS VIA DUMMY DIMENSION ARGUMENTS C COMMON / REAL / R, RN, MMAXR, NNUMR, LLASTR, J COMMON / INTEGER / MMAXI, NNUMI, LLASTI, I, IN, K C C I = VECTOR HOLDING ALL INTEGER ARRAYS (APPENDABLE) C IN = NAMES OF EACH INTEGER ARRAY (APPENDABLE) C J = POINTER TO BEGINNING OF EACH REAL ARRAY (APPENDABLE) C K = POINTER TO BEGINNING OF EACH INTEGER ARRAY (APPENDABLE) C LASTI = NUMBER OF THE LAST ASSIGNED INTEGER ARRAY C LASTR = NUMBER OF THE LAST ASSIGNED REAL ARRAY C LLASTI = NUMBER OF THE LAST ASSIGNED INTEGER ARRAY (COMMON) C LLASTR = NUMBER OF THE LAST ASSIGNED REAL ARRAY (COMMON) C MAXI = STORAGE ALLOWED FOR ALL INTEGER ARRAYS C MAXR = STORAGE ALLOWED FOR ALL REAL ARRAYS C MMAXI = STORAGE ALLOWED FOR ALL INTEGER ARRAYS (COMMON) C MMAXR = STORAGE ALLOWED FOR ALL REAL ARRAYS (COMMON) C NNUMI = ALLOWED NUMBER OF INTEGER ARRAYS (VIA COMMON) C NNUMR = ALLOWED NUMBER OF REAL ARRAYS (VIA COMMON) C NUMI = ALLOWED NUMBER OF INTEGER ARRAYS C NUMR = ALLOWED NUMBER OF REAL ARRAYS C R = VECTOR HOLDING ALL REAL ARRAYS (APPENDABLE) C RN = NAMES OF EACH REAL ARRAY (APPENDABLE) C C ALLOW ACCESS VIA COMMON FOR USER APPLICATIONS (NOT USED) MMAXR = MAXR MMAXI = MAXI NNUMR = NUMR NNUMI = NUMI LLASTI = 0 LLASTR = 0 C C BUILD DYNAMIC DIMENSION STORAGES AND EXECUTE MODEL CALL DYNDIM (NUMR, MAXR, NUMI, MAXI, R, RN, I, IN, 1 J, K ) C NOTE: DYNDIM IS ESSENTIALLY WRITTEN BY PROGRAM DIMMAK.F C USING ARRAYS "REALS" AND "INTEGERS", THEN EDITED C TO REMOVE BLANKS, ETC. IT ALSO WRITES MODEL CALL C AND DIMENSION STATEMENTS STOP 'NORMAL END' END Subroutine notation C ...................... NOTATION ...................... C C AD = VECTOR CONTAINING FLOATING POINT VARIABLES C AJ = JACOBIAN MATRIX C AJINV = INVERSE JACOBIAN MATRIX C C B = STRAIN-DISPLACEMENT (GRADIENT) MATRIX C BODY = BODY FORCE VECTOR C C C = ELEMENT COLUMN MATRIX C CB = BOUNDARY SEGMENT COLUMN MATRIX C CC = COLUMN MATRIX OF SYSTEM EQUATIONS C CEQ = CONSTRAINT EQS COEFFS ARRAY C COORD = SPATIAL COORDINATES OF A SELECTED SET OF NODES C CP = PENALTY CONSTRAINT COLUMN MATRIX C CUTOFF = NUMBER FOR CUTTING OFF ITERATIONS C C D = NODAL PARAMETERS ASSOCIATED WITH A GIVEN ELEMENT C DD = SYSTEM LIST OF NODAL PARAMETERS C DGH = GLOBAL DERIV.S OF INTERPOLATION FUNCTIONS H C DDOLD = SYSTEM LIST OF NODAL DOF FROM LAST ITERATION C DLG = LOCAL DERIVATIVES OF GEOMETRY FUNCTIONS G C DLH = LOCAL DERIVATIVES OF INTERPOLATION FUNCTIONS H C C E = CONSTITUTIVE MATRIX C EB = PRODUCT OF E AND B C ELPROP = ELEMENT ARRAY OF FLOATING POINT PROPERTIES C C FLTNP = REAL PROPERTIES OF SYSTEM NODES C FLTEL = SYSTEM STORAGE OF FLOATING PT ELEMENT PROP C FLTMIS = SYSTEM STORAGE OF FLOATING PT MISC. PROP C FLTNP = SYSTEM STORAGE OF FLOATING PT NODAL PROP C FLUX = SPATIAL COMPONENTS OF SPECIFIED BOUNDARY FLUX C C G = INTERPOLATION FUNCTIONS FOR GEOMETRY C GLOBAL = GLOBAL DERIV.S OF INTERPOLATION FUNCTIONS H C C H = INTERPOLATION FUNCTIONS FOR AN ELEMENT SOLUTION C HINTG = INTEGRAL OF INTERPOLATION FUNCTIONS C C IBC = NODAL POINT BOUNDARY RESTRAINT INDICATOR ARRAY C ID = VECTOR CONTAINING FIXED POINT ARRAYS C IDIAG = DIAGONAL LOCATION IN SKYLINE VECTOR C INDEX = SYSTEM DEGREE OF FREEDOM NUMBERS ARRAY C INRHS > 0, IF INITIAL VALUES OF CC ARE INPUT C IPTEST > 0, IF SOME PROPERTIES ARE DEFINED C ISAY = NO. OF USER REMARKS TO BE I/O C C KFIXED = ALLOCATED SIZE OF ARRAY ID C KFLOAT = ALLOCATED SIZE OF ARRAY AD C KODES = LIST OF DOF RESTRAINT INDICATORS AT A NODE C K1-K5 = NO. OF COLUMNS OF FLOATING PT CONSTRAINT DATA C C LBN = NUMBER OF NODES ON AN ELEMENT BOUNDARY SEGMENT C LEMWRT = 0, IF LIST NODAL PARAMETERS BY ELEMENTS C LHOMO = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS C LNODE = THE N ELEMENT INCIDENCES OF THE ELEMENT C LPFIX = SYSTEM STORAGE ARRAY FOR FIXED PT ELEMENT PROP C LPROP = ARRAY OF FIXED POINT ELEMENT PROPERTIES C LPTEST > 0, IF ELEMENT PROPERTIES ARE DEFINED C C LPPROP = INTEGER PROPERTIES AT EACH ELEMENT NODE C LPROP = ARRAY INTEGER ELEMENT PROPERTIES C LPTEST > 0, IF ELEMENT PROPERTIES HAVE BEEN DEFINED C M = NUMBER OF SYSTEM NODES C MAXACT = NO ACTIVE CONSTRAINT TYPES (<=MAXTYP) C MAXBAN = MAX. HALF BANDWIDTH OF SYSTEM EQUATIONS C MAXTIM > 0, CALCULATE CPU TIMES OF MAJOR SEGMENTS C MAXTYP = MAX NODAL CONSTRAINT TYPE (=3 NOW) C MISCFL = NO. MISC. FLOATING POINT SYSTEM PROPERTIES C MISCFX = NO. MISC. FIXED POINT SYSTEM PROPERTIES C MISFIX = SYSTEM ARRAY OF MISC. FIXED POINT PROPERTIES C MODE = MODE OF STORAGE, 0-SKYLINE, 1-BANDED C MTOTAL = REQUIRED SIZE OF ARRAY AD C M1 TO MNEXT = POINTERS FOR FLOATING POINT ARRAYS C C N = NUMBER OF NODES PER ELEMENT C NCURVE = NO. CONTOUR CURVES CALCULATED PER PARAMETER C NDFREE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM C NDXC = CONSTRAINT EQS DOF NUMBERS ARRAY C NE = NUMBER OF ELEMENTS IN SYSTEM C NELFRE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT C NG = NUMBER OF NODAL PARAMETERS (DOF) PER NODE C NGEOM = NUMBER OF GEOMETRY NODES C NHOMO > 0, IF NODAL SYSTEM PROPERTIES ARE HOMOGENEOUS C NITER = NO. OF ITERATIONS TO BE RUN C NLPFIX = NO. FIXED POINT ELEMENT PROPERTIES C NLPFLO = NO. FLOATING POINT ELEMENT PROPERTIES C NMAT = NUMBER OF MATERIAL TYPES C NNPFIX = NO. FIXED POINT NODAL PROPERTIES C NNPFLO = NO. FLOATING POINT NODAL PROPERTIES C NOCOEF = NO COEFF IN SYSTEM SQ MATRIX C NODES = ELEMENT INCIDENCES OF ALL ELEMENTS C NOTHER = TOTAL NO. OF BOUNDARY RESTRAINTS .GT. TYPE1 C NPARM = DIMENSION OF PARAMWETRIC SPACE C NPFIX = INTEGER PROPERTIES AT ALL NODES C NPROP = NODAL ARRAY OF FIXED POINT PROPERTIES C NPTWRT = 0, LIST NODAL PARAMETERS BY NODES C NQP = NUMBER OF QUADRATURE POINTS C NRANGE = ARRAY CONTAINING NODE NO.S OF EXTREME VALUES C NRB = NUMBER OF ROWS IN B AND E MATRICES C NREQ = NO. OF CONSTRAINT EQS. OF EACH TYPE C NRES = NO. OF CONSTRAINT FLAGS OF EACH TYPE C NSEG = NO OF ELEM BOUNDARY SEGMENTS WITH GIVEN FLUX C NSPACE = DIMENSION OF SPACE C NTAPE1 = UNIT FOR POST SOLUTION MATRICES STORAGE C NTAPE2,3,4 = OPTIONAL UNITS FOR USER (USED WHEN > 0) C NTOTAL = REQUIRED SIZE OF ARRAY ID C NULCOL > 0, IF ELEMENT COLUMN MATRIX IS ALWAYS ZERO C NUMCE = NUMBER OF CONSTRAINT EQS C N1 TO NNEXT = POINTERS FOR FIXED POINT ARRAYS C C PTPROP = NODAL ARRAY OF FLOATING PT PROPERTIES C PRTMAT = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER C PT = QUADRATURE COORDINATES C PRTLPT = FLOATING PT PROP ARRAY OF ELEMENT'S NODES C C RANGE: 1-MAXIMUM VALUE, 2-MINIMUM VALUE OF DOF C C S = ELEMENT SQUARE MATRIX C SB = BOUNDARY SEGMENT SQUARE MATRIX C SS = 'SQUARE' MATRIX OF SYSTEM EQUATIONS C STRAIN = STRAIN OR GRADIENT VECTOR C STRAN0 = INITIAL STRAIN OR GRADIENT VECTOR C STRESS = STRESS VECTOR C C TIME = ARRAY STORING CPU TIMES FOR VARIOUS SEGMENTS C TITLE = PROBLEM TITLE C C USEREL = (USER CHOICE) ELEMENT APPLICATION RESULT C USERPT = (USER CHOICE) NODAL APPLICATION RESULT C C WT = QUADRATURE WEIGHTS C C X = COORDINATES OF SYSTEM NODES C XYZ = SPACE COORDINATES AT A POINT C X = SPATIAL COORDINATES OF ALL NODES IN THE SYSTEM C XPT = SPATIAL COORDINATES OF A CONTOUR POINT C C ...................................................... return end SUBROUTINE APLYBC (MAXACT, NUMCE, NREQ, CEQ, NDXC, 1 NDFREE, NCOEFF, SS, CC, IBW, 2 IDIAG, MODE ) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C APPLY BOUNDARY CONSTRAINT EQUATIONS C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION SS(NCOEFF), CC(NDFREE), NREQ(MAXACT), 1 CEQ(MAXACT,NUMCE), NDXC(MAXACT,NUMCE), 2 IDIAG(NDFREE) C CC = SYS. EQ. COL. MATRIX C CEQ(I,J) = CONSTRAINT EQS COEFF I FOR EQ J C IBW = CURRENT BAND, GROWS WITH CONSTRAINTS C IDIAG = SKY DIAGONAL POINTER FOR EACH DOF C MAXBAN = MAX. HALF-BANDWIDTH OF SYSTEM EQUATIONS C MODE = STORAGE MODE, 0-SKY, 1-BANDED C NCOEFF = NUMBER OF COEFFICIENTS IN SS C NDFREE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM C NDXC(I,J) = CONSTRAINT EQS DOF NO I FOR EQ J C NG = NO. PARAMETERS PER NODE C NREQ = NO OF CONSTRAINT EQS OF EACH TYPE C NUMCE = NUMBER OF CONSTRAINT EQUATIONS C SS = SYS. EQ. SQ. MATRIX IF ( MODE .EQ. 0 ) THEN C SKYLINE MODE c CALL SKYBC (MAXACT, NUMCE, NREQ, CEQ, NDXC, c 1 NDFREE, NCOEFF, SS, CC, IDIAG ) ELSE C BANDED MODE MAXBAN = NCOEFF/NDFREE CALL BANDBC (MAXACT, NUMCE, NREQ, CEQ, NDXC, 1 NDFREE, MAXBAN, SS, CC, IBW ) ENDIF RETURN END SUBROUTINE ASYMBL (NG, NCOEFF, MODE, IDIAG, NODES, SS, CC, 1 M, NE, NDFREE, NITER, LPTEST, LHOMO, NHOMO, NULCOL, N, 2 NSPACE, NELFRE, NRB, NQP, NGEOM, NPARM, NNPFIX, NNPFLO, 3 MISCFX, MISCFL, NLPFIX, NLPFLO, LNODE, INDEX, X, DDOLD, 4 COORD, S, C, H, DGH, B, E, EB, STRAIN, STRAN0, STRESS, 5 BODY, PT, WT, XYZ, DLH, G, DLG, AJ, AJINV, HINTG, D, 6 PRTLPT, FLTNP, FLTEL, FLTMIS, ELPROP, PRTMAT, 7 MISFIX, NPFIX, LPFIX, LPROP, LPPROP, NTAPE1, NTAPE2, 8 NTAPE3, NTAPE4, NTAPE5, LTYPE, NLTYPE, LTDATA, LSHAPE, 9 GPT, GWT ) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C ASSEMBLE SYSTEM EQUATIONS AND STORE POST C SOLUTION ELEMENT DATA C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DATA LASTLT, LTFREE / 0, 0 / DIMENSION CC(NDFREE), SS(NCOEFF), IDIAG(NDFREE), NODES(NE,N), 1 LTYPE(NE), LTDATA(6,NLTYPE) C JUST PASSING THROUGH: SYSTEM DATA DIMENSION X(M,NSPACE), DDOLD(NDFREE), LNODE(N), INDEX(NELFRE) C SYSTEM PROPERTIES DIMENSION PRTLPT(N,0:NNPFLO), FLTNP(M,0:NNPFLO), 1 FLTEL(NE,0:NLPFLO), NPFIX(M,0:NNPFIX), 2 LPFIX(NE,0:NLPFIX) C FOR USE IN ELSQ, ELCOL, OR ELPOST: cb 1 H(N,0:NQP), DGH(NSPACE,N), B(NRB,NELFRE), cb 4 PT(NPARM,0:NQP), WT(0:NQP), DLH(NSPACE,N,0:NQP), cb 7 HINTG(N,0:NQP+1), GPT(0:NQP), GWT(0:NQP), DIMENSION COORD(N,NSPACE), S(NELFRE,NELFRE), C(NELFRE), 1 H(nelfre,0:NQP), DGH(NSPACE,nelfre), B(NRB,NELFRE), 2 E(NRB,NRB), EB(NRB,NELFRE), STRAIN(NRB+2), 3 STRAN0(NRB), STRESS(NRB+2), BODY(NSPACE), 4 PT(NPARM,0:NQP), WT(0:NQP), DLH(NSPACE,nelfre,0:NQP), 5 G(NGEOM,0:NQP), DLG(NPARM,NGEOM,0:NQP), 6 AJ(NSPACE,NSPACE), AJINV(NSPACE,NSPACE), 7 HINTG(nelfre,0:NQP+1), GPT(0:NQP), GWT(0:NQP), 8 XYZ(NSPACE), D(NELFRE), FLTMIS(0:MISCFL), 9 ELPROP(0:NLPFLO), PRTMAT(0:NLPFLO), 1 MISFIX(0:MISCFX), LPROP(0:NLPFIX), 2 LPPROP(0:NNPFIX) C VARIABLES: C AJ = JACOBIAN C AJINV = JACOBIAN INVERSE C B = STRAIN-DISPLACEMENT (GRADIENT) MATRIX C BODY = BODY FORCE VECTOR C CC = SYSTEM EQUATIONS COLUMN MATRIX C COORD = SPATIAL COORDINATES OF ELEMENT'S NODES C D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT C DDOLD = SYSTEM NODAL PARAMETERS FROM LAST ITERATION C DGH = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS C DLG = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION C DLH = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS C E = CONSTITUTIVE MATRIX C EB = PRODUCT OF E*B C ELPROP = ELEMENT ARRAY OF REAL PROPERTIES C FLTEL = REAL PROPERTIES OF ELEMENTS C FLTMIS = MISCELLANEOUS REAL PROPERTIES OF SYSTEM C FLTNP = REAL PROPERTIES OF SYSTEM NODES C G = GEOMETRIC INTERPOLATION FUNCTIONS C H = SOLUTION INTERPOLATION FUNCTIONS C HINTG = INTEGRAL OF INTERPOLATION FUNCTIONS C IDIAG = DIAGONAL LOCATION IN SKYLINE VECTOR C INDEX = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT C LHOMO = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS C LNODE = THE N ELEMENT INCIDENCES OF THE ELEMENT C LPFIX = SYSTEM ARRAY OF INTEGER ELEM PROPERTIES C LPPROP = INTEGER PROPERTIES AT EACH ELEMENT NODE C LPROP = ARRAY INTEGER ELEMENT PROPERTIES C LPTEST > 0, IF ELEMENT PROPERTIES HAVE BEEN DEFINED C LSHAPE = SHAPE FLAG FOR QUADRATURE RULE SELECTION C LTQP = NUMBER OF QUADRATURE PTS FOR ELEMENT TYPE C LTYPE = ELEMENT TYPE NUMBER C M = NUMBER OF SYSTEM NODES C MODE = MODE OF STORAGE, 0-SKYLINE, 1-BANDED C MISFIX = MISCELLANEOUS INTEGER SYSTEM PROPERTIES C N = NUMBER OF NODES PER ELEMENT C NCOEFF = TOTAL NUMBER OF TERMS IN SS C NDFREE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM C NE = NUMBER OF ELEMENTS C NELFRE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT C NG = NUMBER OF PARAMETERS PER NODE C NGEOM = NUMBER OF GEOMETRY NODES C NHOMO = 1, IF NODAL PROPERTIES ARE HOMOGENEOUS C NITER = NO. OF ITERATIONS TO BE RUN (USUALLY 1) C NMAT = NUMBER OF MATERIAL TYPES C NODES = NODAL INCIDENCES OF ALL ELEMENTS C NPARM = DIMENSION OF PARAMWETRIC SPACE C NPFIX = INTEGER PROPERTIES AT ALL NODES C NQP = NUMBER OF QUADRATURE POINTS, >= LTQP C NRB = NUMBER OF ROWS IN B AND E MATRICES C NSPACE = DIMENSION OF SPACE C NTAPE1 = UNIT FOR POST SOLUTION MATRICES STORAGE C NTAPE2,3,4 = OPTIONAL UNITS FOR USER (USED WHEN > 0) C NULCOL > 0, IF ELEMENT COLUMN MATRIX IS ALWAYS ZERO C PRTLPT = REAL PROPERTIES AT ELEMENT NODES C PRTMAT = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER C PT = QUADRATURE COORDINATES C S = ELEMENT SQUARE MATRIX C SS = SYSTEM EQUATIONS SQUARE MATRIX C STRAIN = STRAIN OR GRADIENT VECTOR C STRAN0 = INITIAL STRAIN OR GRADIENT VECTOR C STRESS = STRESS VECTOR C WT = QUADRATURE WEIGHTS C X = COORDINATES OF SYSTEM NODES C XYZ = SPACE COORDINATES AT A POINT C c call at (44) c write(6,*) ntape1,ntape2,ntape3,ntape4,ntape5 C GENERATE ELEMENT EQUATIONS & POST SOLUTION MATRICES DO 10 IE = 1, NE C--> GET ELEMENT TYPE NUMBER LT = 1 IF ( NLTYPE .GT. 1 ) LT = LTYPE(IE) C SAME AS LAST TYPE ? IF ( LT .NE. LASTLT ) THEN LASTLT = LT C GET CONTROLS FOR THIS TYPE CALL GETLT (LT, NLTYPE, LTDATA, LTN, LTQP, LTGEOM, 1 LTPARM, LTSHAP, LTUSER ) LTFREE = LTN*NG C--> GET QUADRATURE RULE FOR ELEMENT TYPE AND SHAPE IF ( LTQP .GT. 0 ) THEN IF ( LTQP .GT. NQP ) STOP 'LTQP > NQP IN ASYMBL' CALL GETQD (LTSHAP, LTQP, NSPACE, GPT, GWT, PT, WT) ENDIF ENDIF C--> EXTRACT ELEMENT NODE NUMBERS CALL LNODES (IE, NE, LTN, NODES, LNODE) C--> CALCULATE DEGREE OF FREEDOM NUMBERS CALL INDXEL (LTN, LTFREE, NG, LNODE, INDEX) c print *, 'ie,ltn,ltfree,ng', ie,ltn,ltfree,ng c call iprint(lnode,1,ltn) c call iprint(index,1,ltfree) c call rprint(d,1,ltfree,1) C--> GENERATE ELEMENT PROBLEM DEPENDENT MATRICES CALL GENELM ( IE, M, NE, NDFREE, NITER, LPTEST, LHOMO, 1 NHOMO, NULCOL, LTN, NSPACE, LTFREE, NRB, LTQP, 2 LTGEOM, LTPARM, NNPFIX, NNPFLO, MISCFX, MISCFL, 3 NLPFIX, NLPFLO, LNODE, INDEX, X, DDOLD, COORD, S, 4 C, H, DGH, B, E, EB, STRAIN, STRAN0, STRESS, BODY, 5 PT, WT, XYZ, DLH, G, DLG, AJ, AJINV, HINTG, D, 6 PRTLPT, FLTNP, FLTEL, FLTMIS, ELPROP, 7 PRTMAT, MISFIX, NPFIX, LPFIX, LPROP, LPPROP, 8 NTAPE1, NTAPE2, NTAPE3, NTAPE4, NTAPE5, LT, 9 LTSHAP, LTUSER, NG ) c print *, 'ie,ltn,ltfree,ng', ie,ltn,ltfree,ng c call iprint(lnode,1,ltn) c call iprint(index,1,ltfree) c call rprint(d,1,ltfree,1) C--> STORE THE MATRICES IN SYSTEM EQUATIONS IF ( MODE .EQ. 0 ) THEN C SKYLINE VECTOR STORAGE MODE c CALL SKYSTR (NCOEFF, NDFREE, LTFREE, INDEX, IDIAG, c 1 S, SS) ELSE C BANDED STORAGE MAXBAN = NCOEFF/NDFREE CALL STORSQ (NDFREE, MAXBAN, LTFREE, INDEX, S, SS) ENDIF IF ( NULCOL .EQ. 0 ) 1 CALL STORCL (NDFREE, LTFREE, INDEX, C, CC) 10 CONTINUE C ASSEMBLY COMPLETED RETURN END SUBROUTINE AT (N) WRITE (6,10) N 10 FORMAT(' -------->>> HERE AT ',I8) RETURN END SUBROUTINE BANCHK (NDFREE, MAXBAN, M, NG, S, C) C * * * * * * * * * * * * * * * * * * * * * * * * * C CHECK BANDED SYSTEM FOR INVALID EQUATIONS & WARN C * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H, O-Z) PARAMETER ( NPRT = 6, ZERO = 0.0 ) DIMENSION S(NDFREE,MAXBAN), C(NDFREE) C C = SYSTEM COLUMN MATRIX C S = SYSTEM SQUARE MATRIX IN BANDED MODE C NDFREE = NUMBER OF EQUATIONS C MAXBAN = HALF BANDWIDTH INCLUDING DIAGONAL SMAX = ZERO DO 10 I = 1, NDFREE TEST = ABS( S(I,1) ) IF ( TEST .GT. SMAX ) SMAX = TEST 10 CONTINUE IF ( SMAX .LE. ZERO ) STOP 1 'ALL ELEMENT STIFFENESSES ZERO, BANCHK' K = 0 DO 20 I = 1, M DO 30 J = 1, NG K = K + 1 TEST = S(K,1) IF ( TEST .LE. ZERO ) THEN IF ( TEST .EQ. ZERO ) WRITE (NPRT,200) I, J 200 FORMAT ('WARNING, NODE ',I5,' DOF',I3,' WAS RESTRAINED') IF ( TEST .LT. ZERO ) WRITE (NPRT,300) I, J 300 FORMAT ('ERROR, NODE ',I5,' DOF',I3,' WAS RESTRAINED') CALL MODFY1 (NDFREE, MAXBAN, K, ZERO, S, C) ENDIF 30 CONTINUE 20 CONTINUE RETURN END SUBROUTINE BANDBC (MAXACT, NUMCE, NREQ, CEQ, NDXC, 1 NDFREE, MAXBAN, SS, CC, IBW ) C * * * * * * * * * * * * * * * * * * * * * * * * * * C APPLY BOUNDARY CONSTRAINT EQUATIONS TO BANDED EQS C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION SS(NDFREE,MAXBAN), CC(NDFREE), NREQ(MAXACT), 1 CEQ(MAXACT,NUMCE), NDXC(MAXACT,NUMCE) C IBW = CURRENT BAND, GROWS WITH CONSTRAINTS C NDFREE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM C MAXBAN = MAX. HALF-BANDWIDTH OF SYSTEM EQUATIONS C SS = SYS. EQ. SQ. MATRIX C CC = SYS. EQ. COL. MATRIX C NG = NO. PARAMETERS PER NODE C NREQ = NO OF CONSTRAINT EQS OF EACH TYPE C CEQ(I,J) = CONSTRAINT EQS COEFF I FOR EQ J C NDXC(I,J) = CONSTRAINT EQS DOF NO I FOR EQ J C NUMCE = NUMBER OF CONSTRAINT EQUATIONS IEQ = 0 c call iprint(nreq,1,maxact) c call iprint(ndxc,maxact,numce) c call rprint(ceq,maxact,numce,0) C DO TYPE ONE LAST c DO 40 IC = MAXACT,1,-1 DO 40 IC = 1, MAXACT NTEST = NREQ(IC) IF ( NTEST .GT. 0 ) THEN IF ( IC .EQ. 1 ) THEN C--> TYPE 1 D(L1) = C1 DO 10 NEQ = 1, NTEST IEQ = IEQ + 1 L1 = NDXC(1,IEQ) C1 = CEQ(1,IEQ) c write(6,*) NDFREE, MAXBAN, IBW, L1, C1 10 CALL MODFY1 (NDFREE, MAXBAN, L1, C1, SS, CC) ELSEIF ( IC .EQ. 2 ) THEN C--> TYPE 2 D(L1)+C1*D(L2)=C2 DO 20 NEQ = 1, NTEST IEQ = IEQ + 1 L1 = NDXC(1,IEQ) L2 = NDXC(2,IEQ) C1 = CEQ(1,IEQ) C2 = CEQ(2,IEQ) c write(6,*) NDFREE, MAXBAN, IBW, L1, L2, C1, C2 20 CALL MODFY2 (NDFREE, MAXBAN, IBW, L1, L2, C1, C2, 1 SS, CC) ELSEIF ( IC .EQ. 3 ) THEN C--> TYPE 3 D(L1)+C1*D(L2)+C2*D(L3)=C3 DO 30 NEQ = 1,NTEST IEQ = IEQ + 1 L1 = NDXC(1,IEQ) L2 = NDXC(2,IEQ) L3 = NDXC(3,IEQ) C1 = CEQ(1,IEQ) C2 = CEQ(2,IEQ) C3 = CEQ(3,IEQ) 30 CALL MODFY3 (NDFREE, MAXBAN, IBW, L1, L2, L3, C1, 1 C2, C3, SS, CC) ELSEIF ( IC .GT. 3 ) THEN C OTHER TYPES NOT DEFINED STOP 'BANLCE NOT INSTALLED, BANDBC' ENDIF ENDIF 40 CONTINUE RETURN END SUBROUTINE BANMLT (NDFREE, MAXBAN, SS, DD, CC, IOPT) C * * * * * * * * * * * * * * * * * * * * * * * * * * C MULTIPLY PACKED SQUARE MATRIX, SS, BY MATRIX DD C IF IOPT = 0 STORE RESULT IN MATRIX CC C OTHERWISE ADD RESULT TO MATRIX CC C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION SS(NDFREE,MAXBAN), DD(NDFREE), CC(NDFREE) C NDFREE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM C MAXBAN = SYSTEM HALF BANDWIDTH MBM1 = MAXBAN - 1 DO 70 I = 1, NDFREE SUM = 0.0 J1 = I - MBM1 J2 = I + MBM1 J1 = MAX0 (J1,1) J2 = MIN0(J2,NDFREE) DO 50 J = J1,J2 IF ( J - I ) 10,20,30 10 JJ = I-J+1 II = J GO TO 40 20 JJ = 1 II = I GO TO 40 30 JJ = J-I+1 II = I 40 SUM = SUM + SS(II,JJ)*DD(J) 50 CONTINUE IF ( IOPT .EQ. 0 ) THEN CC(I) = SUM ELSE CC(I) = CC(I) + SUM ENDIF 70 CONTINUE RETURN END SUBROUTINE BANSUB (I, J, K, L) C * * * * * * * * * * * * * * * * * * * * * * C CONVERT SUBSCRIPTS (I,J) OF SYMMETRIC MATRIX C TO SUBSCRIPTS (K,L) IN UPPER HALF BANDWIDTH C * * * * * * * * * * * * * * * * * * * * * * ITEST = I - J IF ( ITEST ) 10, 20, 30 C BELOW DIAGONAL 10 K = I L = 1 - ITEST RETURN C ON DIAGONAL 20 K = I L = 1 RETURN C ABOVE DIAGONAL 30 K = J L = 1 + ITEST RETURN END SUBROUTINE BAR6T (R, S, T, H, DBH) C * * * * * * * * * * * * * * * * * * * * * * * * * * C BARACENTRIC SHAPES AND DERIVATIVES FOR A T6 ELEMENT C * * * * * * * * * * * * * * * * * * * * * * * * * * PARAMETER ( ZERO = 0.D0, ONE = 1.D0, FOUR = 4.D0 ) DIMENSION H(6), DBH(3,6) C H = SHAPE FUNCTIONS C DBH = BARACENTRIC DERIVATIVES OF H C R,S,T = BARACENTRIC COORDINATES, R+S+T=1 C NODE R S T NODE R S T C 1 1 0 0 4 1/2 1/2 0 3 C 2 0 1 0 5 0 1/2 1/2 6 5 C 3 0 0 1 6 1/2 0 1/2 1 4 2 C SHAPE FUNCTIONS H(1) = R*(R + R - ONE) H(2) = S*(S + S - ONE) H(3) = T*(T + T - ONE) H(4) = FOUR*S*T H(5) = FOUR*R*T H(6) = FOUR*R*S C BARACENTRIC DERIVATIVES DBH(1,1) = FOUR*R - ONE DBH(2,1) = ZERO DBH(3,1) = ZERO DBH(1,2) = ZERO DBH(2,2) = FOUR*S - ONE DBH(3,2) = ZERO DBH(1,3) = ZERO DBH(2,3) = ZERO DBH(3,3) = FOUR*T - ONE DBH(1,4) = ZERO DBH(2,4) = FOUR*T DBH(3,4) = FOUR*S DBH(1,5) = FOUR*T DBH(2,5) = ZERO DBH(3,5) = FOUR*R DBH(1,6) = FOUR*S DBH(2,6) = FOUR*R DBH(3,6) = ZERO RETURN END SUBROUTINE BARPRT (M,NDFREE,NG,NSPACE,IBAR,IPARM, 1 NODIST,X,D,NODBAR) C * * * * * * * * * * * * * * * * * * * * * * * * * * C PRINT-PLOT BAR CHARTS OF NODAL PARAMETER IPARM AT C THE NODES IN ARRAY NODBAR AND SCALE THE RELATIVE C DISTANCE BETWEEN THE POINTS C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) PARAMETER ( NPRT = 6, MID = 25, LINE = 2*MID ) DIMENSION X(M,NSPACE), D(NDFREE), NODBAR(IBAR) CHARACTER ALINE(LINE+1), SKIP(LINE+1) cb DIMENSION ALINE(LINE+1), SKIP(LINE+1) CHARACTER*1 BLANK,DOT,DASH,AX,PLUS,A0 DATA BLANK, DOT, DASH, AX, PLUS, A0 1 / ' ', '+', '-', 'X', '+', '0' / DATA NOPLOT / 0 / C M = TOTAL NUMBER OF NODES IN SYSTEM C IPARM = NODAL PARAMETER TO BE GRAPHED,1<=IPARM<=NG C NODBAR = LIST OF NODES TO BE USED (WHEN IBAR.GT.1) C NODIST = 0 OMIT DISTANCES BETWEEN NODAL BAR LINES C X = ARRAY OF GLOBAL COORDINATES OF ALL NODES C NDFREE = TOTAL NO. OF DEGREES OF FREEDOM IN SYSTEM C D = ARRAY OF ALL NODAL PARAMETERS IN THE SYSTEM C NG = NUMBER OF PARAMETERS PER NODE C NSPACE = DIMENSION OF SOLUTION SPACE C IBAR = NUMBER OF NODES TO BE INCLUDED IN BAR CHART C IF IBAR=1 USE ALL NODES , NODBAR NOT USED CDP SQRT(Z) = DSQRT(Z) CDP ABS(Z) = DABS(Z) NOPLOT = NOPLOT + 1 LIMIT = IBAR IF ( IBAR .EQ. 1 ) LIMIT = M WRITE (NPRT,5000) NOPLOT,IPARM,LIMIT 5000 FORMAT (/, '*** PRINT PLOT NUMBER',I3,' ***',/, 1 'NODAL PARAMETER', I3, 2 ', EVALUATED AT ',I5,' NODE POINTS',/) IPMNG = IPARM - NG IF ( NODIST .EQ. 0 ) GO TO 40 C FIND MINIMUM DISTANCE BETWEEN POINTS DMINSQ = 0.0 DO 10 IS = 1,NSPACE 10 DMINSQ = DMINSQ + ( X(2,IS) - X(1,IS) )**2 DO 30 J = 2,LIMIT I = J IF ( IBAR .GT. 1 ) I = NODBAR(J) DTEST = 0.0 DO 20 IS = 1,NSPACE 20 DTEST = DTEST + ( X(I,IS) - X(I-1,IS) )**2 IF ( DTEST .LT. DMINSQ ) DMINSQ = DTEST 30 CONTINUE DMIN = SQRT(DMINSQ) C--> ESTABLISH GRAPH LIMITING VALUES 40 CONTINUE c YMAX = 0.0 C INITALIZE MAX MIN VALUES OF PARAMETER N = 1 IF ( IBAR .NE. 1 ) N = NODBAR(1) INDEX = NG*N + IPMNG YMAX = D(INDEX) YMIN = YMAX DO 50 I = 1,LIMIT N = I IF ( IBAR .NE. 1 ) N = NODBAR(I) c IF ( IBAR .GT. 1 ) N = NODBAR(I) INDEX = NG*N + IPMNG DTEST = D(INDEX) IF ( DTEST .GT. YMAX ) YMAX = DTEST IF ( DTEST .LT. YMIN ) YMIN = DTEST 50 CONTINUE KOUNT = 1 c 60 IF ( YMAX .LT. 10.0 ) GO TO 70 c KOUNT = KOUNT + 1 c YMAX = YMAX*0.1 c GO TO 60 c 70 YSCALE = 10.**KOUNT c IF ( YMAX .LT. 5.0 ) YSCALE = YSCALE*0.5 c IF ( YMAX .LT. 2.0 ) YSCALE = YSCALE*0.4 c IF ( YMAX .LT. 1.0 ) YSCALE = YSCALE*0.05 WRITE (NPRT,5010) YMIN, YMAX 5010 FORMAT (' RANGE ON GRAPH IS ',1PE12.5,' TO ',1PE12.5/) IF ( YMIN .EQ. YMAX ) THEN WRITE (NPRT,5011) 5011 FORMAT ('CONSTANT VALUE, PLOT SKIPPED') RETURN ENDIF RANGE = YMAX - YMIN C SCALING COMPLETE c WRITE (NPRT,5010) YSCALE c5010 FORMAT (' RANGE ON GRAPH IS +/- ',1PE12.5,/) c CONST = FLOAT(LINE)/(YSCALE + YSCALE) DO 80 I = 2,LINE SKIP(I) = BLANK 80 ALINE(I) = DASH SKIP(1) = PLUS SKIP(MID+1) = DOT SKIP(LINE+1) = PLUS DO 90 I = 1,LINE+1,5 90 ALINE(I) = PLUS ALINE(MID+1) = DOT WRITE (NPRT,5020) ALINE 5020 FORMAT(' NODE VALUE ', (101A1) ) c DO 100 I = 2,LINE c 100 ALINE(I) = BLANK N = 1 IF ( IBAR .GT. 1 ) N = NODBAR(1) DO 210 K = 1,LIMIT NLAST = N N = K IF ( IBAR .GT. 1 ) N = NODBAR(K) INDEX = NG*N + IPMNG DTEST = D(INDEX) c JY = ( D(INDEX) + YSCALE )*CONST + 1.4 DO 100 I = 2,LINE 100 ALINE(I) = BLANK ALINE(1) = PLUS c ALINE(MID+1) = DOT ALINE(LINE+1) = PLUS IF ( JZ .GT. 0 ) ALINE(JZ) = A0 ISPACE = 1 IF ( NODIST .EQ. 0 ) GO TO 120 C--> FIND DISTANCE BETWEEN TWO POINTS DIST = 0.0 DO 110 IS = 1,NSPACE 110 DIST = DIST + ( X(N,IS)-X(NLAST,IS) )**2 C MINIMUM DISTANCE IS ONE SPACE DIST = SQRT(DIST) IF ( K .NE. 1 ) THEN IF ( DMIN .GT. 0.0 ) THEN ISPACE = DIST/DMIN + 0.5 ELSE ISPACE = 1 ENDIF ENDIF IF ( ISPACE .GT. 10 ) ISPACE = 5 120 CONTINUE IF ( NODIST .GT. 0 .AND. K .GT. 1 ) THEN DO 130 I = 1,ISPACE 130 WRITE (NPRT,5030) SKIP 5030 FORMAT (17X, (101A1) ) ENDIF C LINEAR INTERPOLATION FOR COLUMN NUMBER JY = 1*(YMAX - DTEST )/RANGE 1 + (LINE+1)*(DTEST - YMIN)/RANGE + 0.1 C FLAG ZERO VALUE IF ( YMIN .LT. 0.0 .AND. YMAX .GT. 0.0 ) THEN JZ = 1*(YMAX - 0.0 )/RANGE 1 + (LINE+1)*(0.0 - YMIN)/RANGE + 0.1 ELSE JZ = 0 ENDIF c IF ( JY .GT. MID ) GO TO 150 c DO 140 I = JY,MID+1 c 140 ALINE(I) = AX c GO TO 170 c 150 DO 160 I = MID+1,JY c 160 ALINE(I) = AX c 170 CONTINUE DO 165 I = 1, JY 165 ALINE(I) = AX WRITE (NPRT,5040) N, D(INDEX), ALINE 5040 FORMAT (I4,2X,1PE10.3,1X,101A1) c IF ( JY .GT. MID ) GO TO 190 c DO 180 I = JY,MID+1 c 180 ALINE(I) = BLANK c GO TO 210 c 190 DO 200 I = MID+1,JY c 200 ALINE(I) = BLANK 210 CONTINUE DO 220 I = 2,LINE 220 ALINE(I) = DASH DO 230 I = 1,LINE+1,5 230 ALINE(I) = PLUS ALINE(MID+1) = DOT WRITE (NPRT,5020) ALINE RETURN END SUBROUTINE BC2UC (N, NSPACE, DLB, DLU) C * * * * * * * * * * * * * * * * * * * * * * * * * C CONVERT BARACENTRIC LOCAL DERIVATIVES TO C UNIT COORDINATE LOCAL DERIVATIVES C * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION DLU(NSPACE,N), DLB(NSPACE+1,N) C DLB = LOCAL DERIVATIVES IN BARACENTRIC COORDINATES C DLU = LOCAL DERIVATIVES IN UNIT COORDINATES C N = NUMBER OF NODES C NSPACE = DIMENSION OF UNIT SIMPLEX SPACE C LOOP OVER THE INTERPOLATION FUNCTIONS DO 20 K = 1, N DLBLST = DLB(NSPACE+1,K) C SUBTRACT OFF THE LAST BARACENTRIC VALUE DO 10 J = 1, NSPACE 10 DLU(J,K) = DLB(J,K) - DLBLST 20 CONTINUE RETURN END SUBROUTINE BELAST (IOPT, N, NSPACE, NG, GDH, H, 1 R, NS, B) C * * * * * * * * * * * * * * * * * * * * * * * * C ELASTICITY STRAIN-DISPLACEMENT RELATIONS (B) C * * * * * * * * * * * * * * * * * * * * * * * * cb DIMENSION GDH(NSPACE,N), B(NS,N*NG), H(N) DIMENSION GDH(NSPACE,N*NG), B(NS,N*NG), H(N*NG) C B = STRAIN-DISPLACEMENT MATRIX (RETURNED) C GDH = GLOBAL DERIVATIVES OF H C H = ELEMENT INTERPOLATION FUNCTIONS C IOPT = ELASTICITY CLASS C = 1, AXIAL BAR, NG = 1 C = 2, PLANE STRESS, NG = 2 C = 3, PLANE STRAIN, NG = 2 C = 4, AXISYMMETRIC, NG = 2, R = RADIUS C = 5, 3-D SOLID, NG = 3 C N = NUMBER OF NODES PER ELEMENT C NG = NUMBER OF PARAMETERS PER NODE C NS = NUMBER OF STRAINS (ROWS IN B) C NSPACE = DIMENSION OF SPACE IF ( IOPT.LT.1 .OR. IOPT.GT.5 ) STOP 'BELAST' DO 70 J = 1, N K = NG*(J - 1) + 1 L = K + 1 M = L + 1 C--> ONE-DIMENSIONAL AXIAL BAR, IOPT = 1 B(1,K) = GDH(1,J) IF ( IOPT .EQ. 1 ) GO TO 70 C-> PLANE STRESS, PLANE STRAIN, AXISYMMETRIC, 3D 20 B(2,K) = 0.0 B(3,K) = GDH(2,J) B(1,L) = 0.0 B(2,L) = GDH(2,J) B(3,L) = GDH(1,J) IF ( IOPT.EQ.2 .OR. IOPT.EQ.3 ) GO TO 70 C-> AXISYMMETRIC ONLY 30 IF ( IOPT .NE. 4 ) GO TO 40 IF ( R .LE. 0.0 ) STOP 'R=0, BELAST' B(4,K) = H(J)/R B(4,L) = 0.0 GO TO 70 C-> 3D SOLID ONLY 40 B(4,K) = 0.0 B(5,K) = 0.0 B(6,K) = GDH(3,J) B(4,L) = 0.0 B(5,L) = GDH(3,J) B(6,L) = 0.0 B(1,M) = 0.0 B(2,M) = 0.0 B(3,M) = 0.0 B(4,M) = GDH(3,J) B(5,M) = GDH(2,J) B(6,M) = GDH(1,J) 70 CONTINUE RETURN END SUBROUTINE BFLUX (FLUX, COORD, LBN, N, NSPACE, NFLUX, 1 NG, C, S, IOPT, NQP, NPARM, H, DGH, PT, WT, 2 XYZ, DLH, G, DLG, AJ, AJINV, LHOMO, ISEG, 3 NSEG, NBSFIX, NBSFLO, NBSPFX, FLTBS ) C * * * * * * * * * * * * * * * * * * * * * * * * * * C PROBLEM DEPENDENT BOUNDARY FLUX CONTRIBUTIONS C * * * * * * * * * * * * * * * * * * * * * * * * * * C ALWAYS USED DIMENSION COORD(LBN,NSPACE), FLUX(LBN,NG), C(NFLUX), 1 S(NFLUX,NFLUX) C OPTIONAL FOR NUMERICAL INTEGRATION cb DIMENSION H(N), DGH(NSPACE,N), PT(NPARM,0:NQP), cb 1 WT(0:NQP), XYZ(NSPACE), DLH(NSPACE,N), DIMENSION H(N*NG), DGH(NSPACE,N*NG), PT(NPARM,0:NQP), 1 WT(0:NQP), XYZ(NSPACE), DLH(NSPACE,N*NG), 2 G(LBN), DLG(NPARM,LBN), AJ(NSPACE,NSPACE), 3 AJINV(NSPACE,NSPACE) C OPTIONAL SEGMENT PROPERTIES DIMENSION FLTBS(0:NSEG,NBSFLO), NBSPFX(0:NSEG,NBSFIX) C C AJ = JACOBIAN C AJINV = JACOBIAN INVERSE C C = BOUNDARY FLUX COLUMN MATRIX CONTRIBUTIONS C COORD = SPATIAL COORDINATES OF SEGMENT NODES C DGH = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS C DLG = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION C DLH = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS C FLTBS = REAL PROPERTIES ON THE SEGMENTS C FLUX = SPECIFIED BOUNDARY FLUX COMPONENTS C G = GEOMETRIC INTERPOLATION FUNCTIONS C H = SOLUTION INTERPOLATION FUNCTIONS C IOPT = PROBLEM MATRIX REQUIREMENTS, MUST BE SET. C = 1, CALCULATE C ONLY C = 2, CALCULATE S ONLY C = 3, CALCULATE BOTH C AND S C ISEG = SEGMENT NUMBER C LBN = NO. OF NODES ON AN ELEMENT BOUNDARY SEGMENT C LHOMO = 1 IF SEGMENT PROPERTIES ARE HOMOGENEOUS C N = NUMBER OF SOLUTION NODES, N = LBN USUALLY C NBSFIX = NUMBER OF INTEGER PROPERTIES PER SEGMENT C NBSFLO = NUMBER OF REAL PROPERTIES PER SEGMENT C NBSPFX = INTEGER PROPERTIES ON THE SEGMENTS C NFLUX = N*NG = MAXIMUM NUMBER OF FLUX CONTRIBUTIONS C NG = NUMBER OF PARAMETERS PER NODE POINT C NPARM = PARAMETRIC GEOMETRY NODES, = NSPACE USUALLY C NQP = NUMBER OF QUADRATURE POINTS C NSEG = MAX NUMBER OF SEGMENTS C NSPACE = DIMENSION OF SOLUTION SPACE C PT = QUADRATURE COORDINATES C S = BOUNDARY FLUX SQUARE MATRIX C WT = QUADRATURE WEIGHTS C XYZ = SPACE COORDINATES AT A POINT c IOPT = 0 C .................................................... C ** BFLUX PROBLEM DEPENDENT STATEMENTS FOLLOW ** C .................................................... DIMENSION GPT(0:3), GWT(0:3) C T6 OR Q8 OR Q9, LBN = 3, NG = 1 write(6,*) fltbs call rprint (fltbs,nseg+1,nbsflo,0) IF ( LHOMO .EQ. 1) THEN THICK = FLTBS(1,1) ELSE THICK = FLTBS(ISEG,1) ENDIF THICK = 1.0 C GET GAUSS DATA (IN SPACE NPARM - 1) NQ = LBN NQ = MIN0 ( NQ, 3 ) CALL GAUSCO (NQ, GPT, GWT) c write(6,*) 'in bflux' c write(6,*) nparm, nqp c write(6,*) PT c call rprint(pt,nqp+1,nparm,0) c write(6,*) wT c CALL GAUSCO (NQ, PT(1,0), WT) c write(6,*) 'in bflux' c write(6,*) nparm, nqp c write(6,*) PT c call rprint(pt,nqp+1,nparm,0) c write(6,*) wT C ZERO WORKSPACE ARRAY CALL ZEROA (NBFREE,C) C COMPUTE NORMAL FLUX ARRAY BY NUMERICAL INTERGRATION DO 900 IQ = 1, NQ WIQ = GWT(IQ)*THICK PTIQ = GPT(IQ) c write(6,*) IQ, WIQ, ptiq, PT(nparm,IQ) C VARIABLE JACOBIAN ON CURVED LINE CALL DER3L (PTIQ, DLH) DXDR = DOT (LBN, DLH, COORD) DYDR = DOT (LBN, DLH, COORD(1,2)) IF ( DXDR .LT. 0.0 .OR. DYDR .LT. 0.0 ) 1 WRITE (6,*)'WARNING, DISTORTED SEGMENT', ISEG DLDR = SQRT (DXDR*DXDR + DYDR*DYDR) C GET SHAPE, INTERPOLATE NORMAL FLUX, 1---2---3 CALL SHP3L (PTIQ, H) FNORM = DOT (LBN, FLUX, H) C GIVEN HEAT FLUX INWARD DO 54 I =1, LBN 54 C(I) = C(I) + H(I)*WIQ*DLDR*FNORM 900 CONTINUE IOPT = 1 c write(6,*) thick,fnorm,c RETURN END SUBROUTINE BTDB (D, B, S, M, N, IOPT, COEFF) C * * * * * * * * * * * * * * * * * * * * * * C SPECIAL MATRIX MULTIPLICATION OPERATION C IF IOPT=0, S = (B)T*D*B*COEFF PRODUCT C IF IOPT=1, S = (B)T*D*B*COEFF + S NUM. INTG. C * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION D(M,M), B(M,N), S(N,N) C D(M,M) = SYMMETRIC SQUARE MATRIX C B(M,N) = RECTANGLUAR ARRAY C S(N,N) = RETURNED SYMMETRIC SQ MATRIX C COEFF = SCALAR COEFFICIENT DO 40 L = 1,N DO 30 K = 1,N SUM = 0.0 DO 20 I = 1,M DBIK = 0.0 DO 10 J = 1,M C USE SYMMETRY OF D 10 DBIK = DBIK + D(J,I)*B(J,K) SUM = SUM + B(I,L)*DBIK 20 CONTINUE IF ( IOPT .EQ. 0 ) THEN S(L,K) = SUM*COEFF ELSE S(L,K) = S(L,K) + SUM*COEFF ENDIF 30 CONTINUE c 30 S(K,L) = S(L,K) 40 CONTINUE RETURN END SUBROUTINE BTDIAB (DIA, B, S, M, N, IOPT, COEFF) C * * * * * * * * * * * * * * * * * * * * * * * * * * C SPECIAL DIAGONIAL MATRIX MULTIPLICATION OPERATION C IF IOPT=0, S = (B)T*DIA*B*COEFF PRODUCT ONLY C IF IOPT=1, S = (B)T*DIA*B*COEFF + S NUM. INTEG. C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DIA(M), B(M,N), S(N,N) C DIA(M) = DIAGONAL MATRIX C B(M,N) = RECTANGLUAR ARRAY C S(N,N) = RETURNED SYMMETRIC SQ MATRIX C COEFF = SCALAR COEFFICIENT DO 30 L = 1,N DO 20 K = 1,N SUM = 0.0 DO 10 I = 1,M 10 SUM = SUM + B(I,L)*DIA(I)*B(I,K) IF ( IOPT .EQ. 0 ) THEN S(L,K) = SUM*COEFF ELSE S(L,K) = S(L,K) + SUM*COEFF ENDIF 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE CALPRT (N, NNPFLO, H, PRTLPT, VALUES) C * * * * * * * * * * * * * * * * * * * * * * * * * * C CALCULATE NNPFLO PROPERTIES AT A LOCAL PT USING C ELEMENT'S NODAL PROPERTIES, PRTLPT, AND THE N C INTERPOLATION FUNCTIONS, H, AT THE POINT C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION H(N), PRTLPT(N,0:NNPFLO), VALUES(0:NNPFLO) C N = NUMBER OF NODES PER ELEMENT C NNPFLO = NO. OF FLOATING POINT NODAL PROPERTIES C H = C^0 INTERPOLATION FUNCTIONS FOR AN ELEMENT C PRTLPT = FLOATING PT PROPS OF ELEMENT'S NODES C VALUES = LOCAL VALUES OF PROPERTIES IF ( NNPFLO .LT. 1 ) STOP 'NNPFLO = 0, CALPRT' DO 20 I = 1, NNPFLO SUM = 0.0 DO 10 J = 1, N SUM = SUM + H(J)*PRTLPT(J,I) 10 CONTINUE VALUES(I) = SUM 20 CONTINUE RETURN END SUBROUTINE CCOUNT (M, NG, NRES, IBC, KODES, MAXACT, 1 NUMCE, MAXTYP, NREQ) C * * * * * * * * * * * * * * * * * * * * * * * * * * C CALCULATE NUMBER OF CONSTRAINT FLAGS OF EACH TYPE C * * * * * * * * * * * * * * * * * * * * * * * * * * PARAMETER ( NPRT = 6, NBUG = 6) DIMENSION IBC(M), NRES(MAXTYP), KODES(NG), 1 NREQ(MAXTYP) C M = TOTAL NUMBER OF SYSTEM NODES C NG = NO. OF PARAMETERS (DOF) PER NODE C IBC = NODAL POINT BOUNDARY RESTRAINT INDICATOR C KODES = LIST OF RESTRAINT INDICATORS AT A NODE C NRES = LIST OF NUMBER OF FLAGS OF EACH TYPE C = NUMBER OF CONSTR EQS ON EXIT, NREQ C MAXTYP = MAX NO OF DIFFERENT CONSTRAINT TYPES C MAXACT = ACTIVE NO OF TYPES C INITIALIZATION DO 10 I = 1, MAXTYP 10 NRES(I) = 0 DO 30 I = 1, M C DOES NODE I HAVE A NODAL PARAMETER CONSTRAINT ITEST = IABS( IBC(I) ) IF ( ITEST .GT. 0 ) THEN C EXTRACT PARAMETER CODES CALL PTCODE (I,NG,ITEST,KODES) DO 20 J = 1, NG K = KODES(J) C UPDATE CONSTRAINT COUNTERS IF ( K .GT. 0 ) NRES(K) = NRES(K) + 1 20 CONTINUE ENDIF 30 CONTINUE C CONVERT TO EQUATION COUNTERS NUMCE = 0 MAXACT = 1 WRITE (NPRT,5000) 5000 FORMAT ( /, 1 '*** NODAL PARAMETER CONSTRAINT LIST ***', /, 2 'TYPE EQUATIONS') DO 40 I = 1, MAXTYP K = NRES(I) IF ( K .GT. 0 ) MAXACT = I IF ( ((K/I)*I) .LT. K ) WRITE (NBUG,*) 1 'INVALID DATA FOR TYPE', I NREQ(I) = NRES(I)/I IF ( NREQ(I) .GT. 0 ) WRITE (NPRT,5020) I, NREQ(I) 5020 FORMAT ( I4, I10 ) 40 NUMCE = NUMCE + NREQ(I) RETURN END SUBROUTINE CEQBAN (JBW, NREQ, MAXACT, NUMCE, 1 NDXC, NDFREE) C * * * * * * * * * * * * * * * * * * * * * * * * C FIND MAXIMUM HALF BANDWIDTH REQUIRED BY C CONSTRAINT EQUATION MODIFICATION PROCEDURES C * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION NDXC(MAXACT,NUMCE), NREQ(MAXACT) C JBW = MAX HALF BAND FROM CONSTRAINT EQUATIONS C MAXACT = NUMBER ACTIVE CONSTRAINT TYPES C NUMCE = TOTAL NUMBER CONSTRAINT EQS C NDXC(I,J) = CONSTR DOF NO I OF EQ J C NDFREE = TOTAL NO OF SYSTEM DEGREES OF FREEDOM JBW = 1 IEQ = NREQ(1) C--> LOOP OVER NON DIAGONAL CONSTRAINTS DO 30 IC = 2, MAXACT NTEST = NREQ(IC) IF ( NTEST .GT. 0 ) THEN C--> LOOP OVER TYPE IC EQUATIONS DO 20 J = 1, NTEST IEQ = IEQ + 1 IMIN = NDXC(1,IEQ) IMAX = IMIN C--> FIND EQUATION BANDWIDTH DO 10 I = 1, IC INDEX = NDXC(I,IEQ) IF ( INDEX .LT. IMIN ) IMIN = INDEX 10 IF ( INDEX .GT. IMAX ) IMAX = INDEX LBW = IMAX - IMIN + 1 C UPDATE MAXIMUM IF ( LBW .GT. JBW ) JBW = LBW 20 CONTINUE ENDIF 30 CONTINUE RETURN END SUBROUTINE CHANGE (NDFREE, DD, DDOLD, TOTAL, DIFF, 1 RATIO, IPRINT) C * * * * * * * * * * * * * * * * * * * * * * * * * * C CALCULATE THE MEAN CHANGE IN NODAL PARAMETERS FROM C THE LAST ITERATION C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) PARAMETER ( NPRT = 6 ) DIMENSION DD(NDFREE), DDOLD(NDFREE) C * CHANGE SHOULD BE CALLED BEFORE CORECT * C RATIO = DIFF/TOTAL C DIFF = SQRT(SUM OF (DD(I)-DDOLD(I))**2) C TOTAL = SQRT(SUM OF DDOLD(I)**2) C DDOLD = NODAL PARAMETER LIST FROM LAST ITERATION C DD = NODAL PARAMETERS FROM CURRENT ITERATION C NDFREE = TOTAL NO OF DEGREES OF FREEDOM IN SYS C IPRINT > 0, PRINT DIFF, TOTAL, AND RATIO CDP SQRT(Z) = DSQRT(Z) DIFF = 0.0 TOTAL = 0.1E-10 DO 10 I = 1, NDFREE TOTAL = TOTAL + DDOLD(I)*DDOLD(I) 10 DIFF = DIFF + (DD(I)-DDOLD(I))**2 TOTAL = SQRT(TOTAL) DIFF = SQRT(DIFF) RATIO = DIFF/TOTAL IF ( IPRINT .EQ. 0 ) RETURN WRITE (NPRT,5000) DIFF, TOTAL, RATIO 5000 FORMAT ( 1 '*** NODAL DOF FOR CURRENT AND PREVIOUS ITERATIONS ***',/, 2 'ROOT MEAN SQ OF DIFFERENCES . . . . . ',1PE13.5,/, 3 'ROOT MEAN SQ OF PREVIOUS VALUES . . . ',1PE13.5,/, 4 'RATIO OF ABOVE QUANTITIES . . . . . . ',1PE13.5) RETURN END SUBROUTINE CHKSHP (N, NSPACE, LSHAPE, LBN) C * * * * * * * * * * * * * * * * * * * * * * * * * * * C CHECK LSHAPE DATA FOR FREQUENT USER ERRORS C * * * * * * * * * * * * * * * * * * * * * * * * * * * PARAMETER (NPRT = 6 ) C FORCE THE INPUT VALUE IF NEGATIVE IF ( LSHAPE .LT. 0 ) THEN LSHAPE = IABS( LSHAPE ) RETURN ENDIF C LINE IF ( N .EQ. 2 .AND. LSHAPE .NE. 1) THEN LSHAPE = 1 WRITE(NPRT,1000) LSHAPE 1000 FORMAT ('********************************',/, 1 '* WARNING LSHAPE CHANGED TO', I2, ' *', /, 2 '********************************') ENDIF C TRIANGLE IF ( NSPACE .GT. 1 .AND. N .EQ. 3 .AND. LSHAPE .NE. 2) THEN LSHAPE = 2 WRITE(NPRT,1000) LSHAPE ENDIF C QUADRILATERAL IF ( NSPACE .GT. 1 .AND. N .EQ. 4 .AND. LSHAPE .NE. 3) THEN LSHAPE = 3 WRITE(NPRT,1000) LSHAPE ENDIF C TRI OR WEDGE IF ( N .EQ. 6 .AND. LSHAPE .NE. 2) THEN IF ( NSPACE .EQ. 2 ) THEN LSHAPE = 2 WRITE(NPRT,1000) LSHAPE ENDIF ENDIF C QUAD OR HEX IF ( N .EQ. 8 .AND. LSHAPE .NE. 3) THEN IF ( NSPACE .EQ. 3 ) THEN LSHAPE = 3 WRITE(NPRT,1000) LSHAPE ENDIF ENDIF RETURN END SUBROUTINE CONDSE (NTOTAL, NELFRE, S, C) C * * * * * * * * * * * * * * * * * * * * * * * C CONDENSATION OF ELEMENT MATRICES TO REMOVE C INTERNAL DEGREES OF FREEDOM C * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8(A-H,O-Z) PARAMETER ( ZERO = 0.0 ) DIMENSION S(NTOTAL,NTOTAL), C(NTOTAL) C INTERNAL DEGREES OF FREEDOM *MUST* COME LAST C : SAA : SAB : : DA : : CA : C :......:......: :....: = :....: C : SBA : SBB : : DB : : CB : C ENTER FULL S ; RETURN CONDENSED IN SAA AND CA C DIMENSION SAA(NELFRE,NELFRE), CA(NELFRE) C SAA* = (SAA) - (SAB)*(SBB)I*(SAB)T C CA* = (CA) - (SAB)*(SBB)I*(CB) C NTOTAL = ORIG. NO. OF D.O.F. OF ELEMENT C NELFRE = FINAL NO. OF D.O.F. OF ELEMENT C NELIM = NUMBER OF DOF TO ELIMINATE C S = SQUARE ELEMENT MATRIX C C = ELEMENT COLUMN MATRIX NELIM = NTOTAL - NELFRE DO 30 I = 1, NELIM J = NTOTAL - I K = J + 1 SKK = S(K,K) CK = C(K) IF ( SKK .NE. ZERO ) THEN C(K) = CK/SKK DO 20 L = 1, J SLKSKK = S(L,K)/SKK S(L,K) = SLKSKK DO 10 M = L, J S(L,M) = S(L,M) - S(K,M)*SLKSKK 10 S(M,L) = S(L,M) C(L) = C(L) - CK*SLKSKK 20 CONTINUE ENDIF 30 CONTINUE RETURN END SUBROUTINE CONTROL (TITLE, M, NE, NG, N, NSPACE, NSEG, LBN, 1 NITER, NCURVE, INRHS, ISAY, NNPFIX, NNPFLO, NLPFIX, 2 NLPFLO, MISCFX, MISCFL, NHOMO, LHOMO, NPTWRT, LEMWRT, 3 NTAPE1, NTAPE2, NTAPE3, NTAPE4, NTAPE5, NULCOL, 4 NDFREE, NELFRE, NFLUX, IPTEST, LPTEST, NRB, NQP, 5 LSHAPE, NLTYPE, MODE, IBUG, NBSFIX, NBSFLO, NGF) C 1 2 3 4 5 6 712 C23456789012345678901234567890123456789012345678901234567890-----------X cb DIMENSION TITLE(15) CHARACTER*4 TITLE(15) C PRINT AUTHOR CREDITS C CALL TOOT C--> ** READ AND PRINT TITLE AND CONTROL DATA ** NTAPE3 = 0 NTAPE4 = 0 NTAPE5 = 0 READ (5,1234) TITLE WRITE (6,1234) TITLE 1234 FORMAT ( 15A4 ) READ (5,5020) M, NE, NG, N, 1 NSPACE, NSEG, LBN, NITER, 2 NCURVE, INRHS, ISAY, NRB, 3 NQP, LSHAPE, NLTYPE, MODE 5020 FORMAT ( 16I5 ) C CHECK FOR EXTRA OUTPUT REQUESTS IBUG = 0 IF ( NITER .LT. 0 ) THEN IBUG = 1 NITER = IABS ( NITER ) ENDIF IF ( LSHAPE .LT. 1 ) LSHAPE = 1 IF ( N .LT. 1 ) N = 2 IF ( NG .LT. 1 ) NG = 1 IF ( NITER .LT. 1 ) NITER = 1 IF ( NLTYPE .LT. 1 ) NLTYPE = 1 IF ( NRB .LT. 1 ) NRB = 1 IF ( NSPACE .LT. 1 ) NSPACE = 1 C CHECK ELEMENT SHAPE DATA CALL CHKSHP (N, NSPACE, LSHAPE, LBN) MODE = 1 IF ( NITER .LT. 1 ) NITER = 1 WRITE (6,5030) M, NE, NG, N, NSPACE, NSEG, LBN, 1 NITER, NCURVE, INRHS, ISAY, NRB, NQP, 2 LSHAPE, NLTYPE, MODE 5030 FORMAT (/,'**** PROBLEM PARAMETERS (DEFAULT) ****',/, 1 'NUMBER OF NODAL POINTS IN SYSTEM ............',I5,/, 2 'NUMBER OF ELEMENTS IN SYSTEM ................',I5,/, 4 'NUMBER OF PARAMETERS PER NODE (1)............',I5,/, 3 'NUMBER OF NODES PER ELEMENT (2)..............',I5,/, 5 'DIMENSION OF SPACE (1).......................',I5,/, 6 'NUMBER OF BOUNDARIES WITH GIVEN FLUX ........',I5,/, 7 'NUMBER OF NODES ON BOUNDARY SEGMENT .........',I5,/, 8 'NUMBER OF ITERATIONS TO BE RUN (1)...........',I5,/, 9 'NUMBER OF CONTOURS BETWEEN 5 & 95% ..........',I5,/, + 'INITIAL FORCING VECTOR (0-OMIT, 1-READ)......',I5,/, 1 'NUMBER OF USER REMARKS LINES ................',I5,/, 2 'NUMBER OF ROWS IN B MATRIX (1)...............',I5,/, 3 'NUMBER OF QUADRATURE POINTS .................',I5,/, 4 'SHAPE 1-LINE 2-TRI 3-QUAD 4-HEX 5-TET (1)....',I5,/, 5 'NUMBER OF DIFFERENT ELEMENT TYPES (1)........',I5,/, 6 'STIFFNESS STORAGE MODE: (0-SKY, 1-BAND)......',I5) READ (5,5070) NNPFIX, NNPFLO, NLPFIX, NLPFLO, 1 MISCFX, MISCFL, NHOMO, LHOMO, 2 NPTWRT, LEMWRT, NTAPE1, NTAPE2, 3 NGF, NULCOL, NBSFIX, NBSFLO 5070 FORMAT ( 16I5 ) IF ( NSEG .GT. 0 .AND. NGF .LT. 1 ) THEN WRITE (6,*) 'WARNING, NGF < 1, SET TO NGF = NG' NGF = NG ENDIF WRITE (6,5080) NGF, NNPFIX, NNPFLO, NLPFIX, NLPFLO, 1 NBSFIX, NBSFLO, MISCFX, MISCFL 5080 FORMAT ( 8 'NUMBER OF FLUX COMPONENTS PER NODE (NG)......',I5,/, 1 'NUMBER OF INTEGER PROPERTIES PER NODE .......',I5,/, 2 'NUMBER OF REAL PROPERTIES PER NODE ..........',I5,/, 3 'NUMBER OF INTEGER PROPERTIES PER ELEMENT ....',I5,/, 4 'NUMBER OF REAL PROPERTIES PER ELEMENT .......',I5,/, 3 'NUMBER OF INTEGER PROPERTIES PER SEGMENT ....',I5,/, 4 'NUMBER OF REAL PROPERTIES PER SEGMENT .......',I5,/, 5 'NUMBER OF INTEGER MISCELLANEOUS PROPERTIES .',I5,/, 6 'NUMBER OF REAL MISCELLANEOUS PROPERTIES ....',I5) IF ( LBN .GT. N ) 1 STOP 'INCONSISTANT VALUES OF LBN AND N.' NELFRE = N*NG NDFREE = M*NG NFLUX = LBN*NG WRITE (6,5081) NELFRE, NFLUX, NDFREE 5081 FORMAT ( 1 'NUMBER OF D.O.F. FOR ELEMENT ................',I5,/, 2 'NUMBER OF D.O.F. ON FLUX SEGMENT ............',I5,/, 3 'NUMBER OF D.O.F. IN TOTAL SYSTEM ............',I5) IF ( NHOMO .EQ. 1 ) WRITE (6,*) 1 'NODAL POINT PROPERTIES ARE HOMOGENEOUS.' IF ( LHOMO .EQ. 1 ) WRITE (6,*) 1 'ELEMENT PROPERTIES ARE HOMOGENEOUS.' NSUM = NTAPE1 + NTAPE2 + NTAPE3 + NTAPE4 + NTAPE5 IF ( NSUM .GT. 0 ) WRITE (6,5180) 1 NTAPE1, NTAPE2, NTAPE3, NTAPE4, NTAPE5 5180 FORMAT ( /, 'OPTIONAL UNIT NUMBERS (UTILIZED IF > 0)',/, 1 'NTAPE1 = ',I2,', NTAPE2 = ',I2,/,'NTAPE3 = ',I2, 2 ', NTAPE4 = ',I2,', NTAPE5 = ',I2) IF ( NPTWRT .EQ. 0 ) WRITE (6,*) 1 'NODAL PARAMETERS TO BE LISTED BY NODES' IF ( LEMWRT .EQ. 0 ) WRITE (6,*) 1 'NODAL PARAMETERS TO BE LISTED BY ELEMENTS' IF ( NULCOL .NE. 0 ) WRITE (6,*) 1 'ALL ELEMENT COLUMN MATRICES ARE ZERO.' cb IF ( NTAPE1 .GT. 0 ) REWIND NTAPE1 cb IF ( NTAPE2 .GT. 0 ) REWIND NTAPE2 cb IF ( NTAPE3 .GT. 0 ) REWIND NTAPE3 cb IF ( NTAPE4 .GT. 0 ) REWIND NTAPE4 cb IF ( NTAPE5 .GT. 0 ) REWIND NTAPE5 IF ( ISAY .GT. 0 ) CALL IREMRK (ISAY) C SET INITIAL CONSTANTS C LPTEST > 0, ELEMENT PROPERTIES ARE DEFINED C IPTEST > 0, SOME PROPERTIES ARE DEFINED C NBSFIX = NUMBER OF FIXED PT SEGMENT PROP C NBSFLO = NUMBER OF FLOATING PT SEGMENT PROP C NLPFIX = NUMBER OF FIXED PT ELEMENT PROP C NLPFLO = NUMBER OF FLOATING PT ELEMENT PROP C NNPFIX = NUMBER OF FIXED PT NUMBER PROP C NNPFLO = NUMBER OF FLOATING PT NUMBER PROP C MISCFL = NUMBER OF MISC FLOATING PT SYSTEM PROP C MISCFX = NUMBER OF MISC FIXED PT SYSTEM PROP C MAXTYP = MAX ALLOWED CONSTRAINT TYPE C RATIO = CONSTANT FOR ITER CONTROL, SEE MODEL RATIO = 1.0 C MAXTYP = 5 C IF ( NFLUX .LT. 1 ) NFLUX = 1 IPTEST = NNPFIX + NNPFLO + NLPFIX + NLPFLO 1 + NBSFIX + NBSFLO + MISCFX + MISCFL LPTEST = NLPFIX + NLPFLO RETURN END SUBROUTINE CORECT (NDFREE, DD, DDOLD) C * * * * * * * * * * * * * * * * * * * * * * * * * * C CALCULATE NEW STARTING VALUES FOR NEXT ITERATION C * * * * * * * * * * * * * * * * * * * * * * * * * * C OVER RELAXATION METHOD CDP IMPLICIT REAL*8(A-H,O-Z) DIMENSION DD(NDFREE), DDOLD(NDFREE) PARAMETER ( OMEGA = 1.25 ) C DD = CALCULATED DOF FROM LAST ITERATION C DDOLD = DOF TO BE USED TO START NEXT ITERATION C NDFREE = TOTAL NO OF SYS DEGREES OF FREEDOM DO 10 I = 1,NDFREE 10 DDOLD(I) = DDOLD(I) + OMEGA*(DD(I)-DDOLD(I)) RETURN END SUBROUTINE DCHECK (DELTA,N,NSPACE) C * * * * * * * * * * * * * * * * * * * * * * * * * * C CHECKING OF THE LOCAL COORDINATE DERIVATIVES OF THE C N SHAPE FUNCTIONS AT A LOCAL POINT FOR A C0 ELEMENT C * * * * * * * * * * * * * * * * * * * * * * * * * * DOUBLE PRECISION ONE, SUM, TOL PARAMETER ( ONE = 1.0D0, TOL = 1.0D-7, NPRT = 6 ) DIMENSION DELTA(NSPACE,N) C DELTA = LOCAL DERIVATIVES OF SHAPE FUNCTIONS C N = NUMBER OF SHAPE FUNCTIONS C NSPACE = DIMENSION OF LOCAL SPACE IERR = 0 DO 20 J = 1,NSPACE SUM = 0.D0 DO 10 I = 1,N 10 SUM = SUM + DELTA(J,I) IF ( ABS(SUM) .GT. TOL ) THEN IF ( IERR .EQ. 0 ) WRITE (NPRT,*) & 'SUPPLIED DERIVATIVES ARE INCORRECT' IERR = 1 WRITE (NPRT,*) 'J, SUM', J, SUM ENDIF 20 CONTINUE IF ( IERR .NE. 0 ) THEN CALL RPRINT (DELTA,N,NSPACE,1) WRITE (NPRT,*) 'END OF WARNING FROM DCHECK' ENDIF RETURN END SUBROUTINE DEGPAR (IPT, JPARM, NG, INDEX) C * * * * * * * * * * * * * * * * * * * * * * * C DETERMINE THE DEGREE OF FREEDOM NUMBER C OF NODAL PARAMETER JPARM AT NODE POINT IPT C * * * * * * * * * * * * * * * * * * * * * * * C NG = NUMBER OF PARAMETERS PER NODE INDEX = NG*(IPT-1) + JPARM RETURN END SUBROUTINE DELAST (IOPT, E, PR, T, D, NS) C * * * * * * * * * * * * * * * * * * * * * * * C CONSTITUTIVE MATRIX, ELASTICITY (D) C * * * * * * * * * * * * * * * * * * * * * * * DIMENSION D(NS,NS) C D = CONSTITUTIVE MATRIX C E = MODULUS OF ELASTICITY C IOPT = ELASTICITY CLASS C = 1, AXIAL BAR, T = AREA C = 2, PLANE STRESS, T = THICKNESS C = 3, PLANE STRAIN, T = THICKNESS C = 4, AXISYMMETRIC C = 5, 3-D SOLID C NS = NUMBER OF STRAINS (ROWS IN B-MATRIX) C PR = POISSON'S RATIO C T = AREA, OR THICKNESS IF ( T.LE.0.0 ) T = 1.0 IF ( IOPT.LT.1 .OR. IOPT.GT.5 ) STOP 'DELAST' IF ( IOPT.NE.1 ) GO TO 20 C 1-D, SXX D(1,1) = E*T RETURN 20 IF ( IOPT.NE.2 ) GO TO 30 C PLANE STRESS ONLY, SXX, SYY, SXY C = T*E/(1.-PR*PR) D(1,1) = C D(2,1) = C*PR D(3,1) = 0.0 D(1,2) = C*PR D(2,2) = C D(3,2) = 0.0 D(1,3) = 0.0 D(2,3) = 0.0 D(3,3) = 0.5*T*E/(1.+PR) RETURN 30 CONTINUE C PLANE STRAIN OR 3-D, SXX, SYY, SXY C = E*(1.-PR)/(1.+PR)/(1.-PR-PR) C2 = C*PR/(1.-PR) G = 0.5*E/(1.+PR) D(1,1) = C D(2,1) = C2 D(3,1) = 0.0 D(1,2) = C2 D(2,2) = C D(3,2) = 0.0 D(1,3) = 0.0 D(2,3) = 0.0 D(3,3) = G IF ( IOPT.EQ.3 ) RETURN C AXISYMMETRIC ONLY, SXX, SYY, SXY, STT D(4,1) = C2 D(4,2) = C2 D(4,3) = 0.0 D(1,4) = C2 D(2,4) = C2 D(3,4) = 0.0 D(4,4) = C IF ( IOPT.NE.5 ) RETURN C 3-D SOLID ONLY, SXX, SYY, SXY, SZZ, SXZ, SYZ D(5,1) = 0.0 D(6,1) = 0.0 D(5,2) = 0.0 D(6,2) = 0.0 D(5,3) = 0.0 D(6,3) = 0.0 D(5,4) = 0.0 D(6,4) = 0.0 D(1,5) = 0.0 D(2,5) = 0.0 D(3,5) = 0.0 D(4,5) = 0.0 D(5,5) = G D(6,5) = 0.0 D(1,6) = 0.0 D(2,6) = 0.0 D(3,6) = 0.0 D(4,6) = 0.0 D(5,6) = 0.0 D(6,6) = G RETURN END SUBROUTINE DER16QS (R,S,DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C LOCAL DERIVATIVES FOR A SERENDIPITY 16 NODE QUAD C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * PARAMETER ( PT667 = 0.66666666666667 ,PT1D12 = 0.083333333333333 1 ,HALF = 0.5) C SEE SHP16QS FOR NODE LOCATIONS DIMENSION DH(2,16) RR = R*R SS = S*S RRR = R*R*R SSS = S*S*S RP = 1. + R RM = 1. - R R1 = R - 1. SP = 1. + S SM = 1. - S DH(1,1) = PT1D12*SM*(16.*RRR - 12.*RR - 2.*R + 4.*SSS - S + 4.) DH(1,2) = PT1D12*SM*(16.*RRR + 12.*RR - 2.*R - 4.*SSS + S - 4.) DH(1,3) = PT1D12*SP*(16.*RRR + 12.*RR - 2.*R + 4.*SSS - S - 4.) DH(1,4) = PT1D12*SP*(16.*RRR - 12.*RR - 2.*R - 4.*SSS + S + 4.) DH(1,5) = PT667*(4.*R - 1. + 3.*RR - 8.*RRR)*SM DH(1,6) = PT667*S*SM*SP*( -1. + 2.*S) DH(1,7) = PT667*(1. + 4.*R - 3.*RR - 8.*RRR)*SP DH(1,8) = -PT667*S*SM*SP*(1. + 2.*S) DH(1,9) = R*(-5. + 8.*RR)*SM DH(1,10) = HALF*SM*SP*(1. - 2.*S)*(1. + 2.*S) DH(1,11) = R*( -5. + 8.*RR)*SP DH(1,12) = HALF*SM*SP*( -1. + 2.*S)*(1. + 2.*S) DH(1,13) = PT667*(1. + 4.*R - 3.*RR - 8.*RRR)*SM DH(1,14) = PT667*S*SM*SP*(1. + 2.*S) DH(1,15) = PT667*(4.*R - 1. + 3.*RR - 8.*RRR)*SP DH(1,16) = PT667*S*SM*SP*(1. - 2.*S) DH(2,1) = PT1D12*RM*(16.*SSS - 12.*SS - 2.*S + 4.*RRR - R + 4.) DH(2,2) = PT1D12*RP*(16.*SSS - 12.*SS - 2.*S - 4.*RRR + R + 4.) DH(2,3) = PT1D12*RP*(16.*SSS + 12.*SS - 2.*S + 4.*RRR - R - 4.) DH(2,4) = PT1D12*R1*(-16.*SSS - 12.*SS + 2.*S + 4.*RRR - R + 4.) DH(2,5) = PT667*R*R1*RP*(2.*R - 1.) DH(2,6) = PT667*( -1. + 4.*S + 3.*SS - 8.*SSS)*RP DH(2,7) = PT667*R*RM*RP*(1. + 2.*R) DH(2,8) = PT667*( -1. - 4.*S + 3.*SS + 8.*SSS)*R1 DH(2,9) = HALF*RM*RP*(2.*R - 1.)*(1. + 2.*R) DH(2,10) = S*( -5. + 8.*SS)*RP DH(2,11) = HALF*R1*RP*(2.*R - 1.)*(1. + 2.*R) DH(2,12) = S*(5. - 8.*SS)*R1 DH(2,13) = PT667*R*R1*RP*(1. + 2.*R) DH(2,14) = PT667*(1. + 4.*S - 3.*SS - 8.*SSS)*RP DH(2,15) = PT667*R*RM*RP*(2.*R - 1.) DH(2,16) = PT667*(1. - 4.*S - 3.*SS + 8.*SSS)*R1 RETURN END SUBROUTINE DER16R (R,S,A,B,DH) C * * * * * * * * * * * * * * * * * * * * * * * * * C WARNING, this code is not fully checked. No known bugs. C FIRST DERIVATIVES OF A C C1 RECTANGULAR ELEMENT IN UNIT COORDINATES C USING TENSOR PRODUCTS OF 1D BASIS C * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION DH(2,16), HR(4), DHR(4), HS(4), DHS(4) C DOF ARE W W,X W,Y W,XY AT EACH NODE (NG=4) C X // R, Y // S. S C A = PHYSICAL LENGTH IN X 4 -------- 3 C B = PHYSICAL LENGTH IN Y I I C R,S = LOCAL UNIT COORDS I I C 1@(0,0), 3@(1,1) 1 -------- 2 ->R C C Evaluate the 1D interpolations CALL SHPC1L (R,A,HR) CALL SHPC1L (S,B,HS) CALL DERC1L (R,A,DHR) CALL DERC1L (S,B,DHS) C Form tensor products for R direction DH(1,1) = DHR(1)*HS(1) DH(1,2) = DHR(2)*HS(1) DH(1,3) = DHR(1)*HS(2) DH(1,4) = DHR(2)*HS(2) DH(1,5) = DHR(3)*HS(1) DH(1,6) = DHR(4)*HS(1) DH(1,7) = DHR(3)*HS(2) DH(1,8) = DHR(4)*HS(2) DH(1,9) = DHR(3)*HS(3) DH(1,10) = DHR(4)*HS(3) DH(1,11) = DHR(3)*HS(4) DH(1,12) = DHR(4)*HS(4) DH(1,13) = DHR(1)*HS(3) DH(1,14) = DHR(2)*HS(3) DH(1,15) = DHR(1)*HS(4) DH(1,16) = DHR(2)*HS(4) C Form tensor products for S direction DH(2,1) = HR(1)*DHS(1) DH(2,2) = HR(2)*DHS(1) DH(2,3) = HR(1)*DHS(2) DH(2,4) = HR(2)*DHS(2) DH(2,5) = HR(3)*DHS(1) DH(2,6) = HR(4)*DHS(1) DH(2,7) = HR(3)*DHS(2) DH(2,8) = HR(4)*DHS(2) DH(2,9) = HR(3)*DHS(3) DH(2,10) = HR(4)*DHS(3) DH(2,11) = HR(3)*DHS(4) DH(2,12) = HR(4)*DHS(4) DH(2,13) = HR(1)*DHS(3) DH(2,14) = HR(2)*DHS(3) DH(2,15) = HR(1)*DHS(4) DH(2,16) = HR(2)*DHS(4) RETURN END SUBROUTINE DER17Q (R,S,DH) C ****************************************************************** C LOCAL DERIVATIVES OF SERENDIPITY QUAD WITH 17 NODES C ****************************************************************** PARAMETER ( PT667 = 0.66666666666667 ,PT1D12 = 0.083333333333333 1 ,HALF = 0.5) C SEE SHP17Q FOR NODE LOCATIONS DIMENSION DH(2,17) RP = 1. + R RM = 1. - R SP = 1. + S SM = 1. - S R1 = R - 1. RR = R*R SS = S*S RRR = R*R*R SSS = S*S*S DH(1,1) = PT1D12*SM*(16.*RRR-12.*RR-6.*R*S-8.*R+4.*SSS-S+4.) DH(1,2) = PT1D12*SM*(16.*RRR+12.*RR-6.*R*S-8.*R-4.*SSS+S-4.) DH(1,3) = PT1D12*SP*(16.*RRR+12.*RR+6.*R*S-8.*R+4.*SSS-S-4.) DH(1,4) = PT1D12*SP*(16.*RRR-12.*RR+6.*R*S-8.*R-4.*SSS+S+4.) DH(1,5) = -PT667*(1. - 4.*R - 3.*RR + 8.*RRR)*SM DH(1,6) = PT667*S*SM*SP*( -1. + 2.*S) DH(1,7) = -PT667*(-1. - 4.*R + 3.*RR + 8.*RRR)*SP DH(1,8) = -PT667*S*SM*SP*(1. + 2.*S) DH(1,9) = R*SM*(8.*RR + S - 4.) DH(1,10) = HALF*SM*SP*(2.*R - 4.*SS + 1.) DH(1,11) = R*SP*(8.*RR - S - 4.) DH(1,12) = HALF*SM*SP*(2.*R - 1. + 4.*SS) DH(1,13) = PT667*(1. + 4.*R - 3.*RR - 8.*RRR)*SM DH(1,14) = PT667*S*SM*SP*(1. + 2.*S) DH(1,15) = -PT667*(1. - 4.*R - 3.*RR + 8.*RRR)*SP DH(1,16) = PT667*S*SM*SP*(1. - 2.*S) DH(1,17) = -2.*R*SM*SP DH(2,1) = PT1D12*RM*(16.*SSS-12.*SS-6.*R*S-8.*S+4.*RRR-R+4.) DH(2,2) = -PT1D12*RP*(-16.*SSS+12.*SS-6.*R*S+8.*S+4.*RRR-R-4.) DH(2,3) = PT1D12*RP*(16.*SSS+12.*SS+6.*R*S-8.*S+4.*RRR-R-4.) DH(2,4) = PT1D12*R1*(-16.*SSS-12.*SS+6.*R*S+8.*S+4.*RRR-R+4.) DH(2,5) = PT667*R*R1*RP*(2.*R - 1.) DH(2,6) = -PT667*(1. - 4.*S - 3.*SS + 8.*SSS)*RP DH(2,7) = -PT667*R*R1*RP*(1. + 2.*R) DH(2,8) = PT667*(-1. - 4.*S + 3.*SS + 8.*SSS)*R1 DH(2,9) = -HALF*R1*RP*(2.*S - 1. + 4.*RR) DH(2,10) = -S*RP*(-8.*SS + R + 4.) DH(2,11) = HALF*R1*RP*( -2.*S + 4.*RR - 1.) DH(2,12) = -S*R1*(8.*SS + R - 4.) DH(2,13) = PT667*R*R1*RP*(1. + 2.*R) DH(2,14) = -PT667*( -1. - 4.*S + 3.*SS + 8.*SSS)*RP DH(2,15) = -PT667*R*R1*RP*(2.*R - 1.) DH(2,16) = PT667*(1. - 4.*S - 3.*SS + 8.*SSS)*R1 DH(2,17) = 2.*S*R1*RP RETURN END SUBROUTINE DER208 (R, S, T, DH, LNODE) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C LOCAL DERIVATIVES OF INTERPOLATION FUNCTIONS FOR AN C 8 TO 20 NODE HEXAHEDRON, SEE SHP208 FOR TOPOLOGY FIGURE C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * CIBM IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DH(3,20), LNODE(20), I1(20), I2(20) DATA I1 /8*0, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4/ DATA I2 /8*0, 2, 3, 4, 1, 6, 7, 8, 5, 5, 6, 7, 8/ DATA FP, FM / 0.5, -0.5 / C FOR PARAMETER DEFINITIONS SEE SUBROUTINE SHP208 C R,S,T = LOCAL COORDINATES OF THE POINT -1 LE (R,S,T) LE +1 C DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS, 0 IF LNODE(I) = 0 C LNODE = ARRAY OF ELEMENT INCIDENCES, C IF LNODE(I)=0 THEN LOCAL NODE I IS NOT CONSIDERED IN ANALYSIS C I1, I2 = CORNER NODES OF TWELVE EDGES RP = 0.5*(1. + R) SP = 0.5*(1. + S) TP = 0.5*(1. + T) RM = 0.5*(1. - R) SM = 0.5*(1. - S) TM = 0.5*(1. - T) RZ = 1. - R*R SZ = 1. - S*S TZ = 1. - T*T FR = -2.0*R FS = -2.0*S FT = -2.0*T C DERIVATIVES OF TRI-LINEAR CORNERS DH(1,1 ) = TP*SP*FP DH(2,1 ) = TP*FP*RP DH(3,1 ) = FP*SP*RP DH(1,2 ) = TP*SP*FM DH(2,2 ) = TP*FP*RM DH(3,2 ) = FP*SP*RM DH(1,3 ) = TP*SM*FM DH(2,3 ) = TP*FM*RM DH(3,3 ) = FP*SM*RM DH(1,4 ) = TP*SM*FP DH(2,4 ) = TP*FM*RP DH(3,4 ) = FP*SM*RP DH(1,5 ) = TM*SP*FP DH(2,5 ) = TM*FP*RP DH(3,5 ) = FM*SP*RP DH(1,6 ) = TM*SP*FM DH(2,6 ) = TM*FP*RM DH(3,6 ) = FM*SP*RM DH(1,7 ) = TM*SM*FM DH(2,7 ) = TM*FM*RM DH(3,7 ) = FM*SM*RM DH(1,8 ) = TM*SM*FP DH(2,8 ) = TM*FM*RP DH(3,8 ) = FM*SM*RP C DERIVATIVES OF EDGE BUBBLES DH(1,9 ) = TP*SP*FR*0.5 DH(2,9 ) = TP*FP*RZ*0.5 DH(3,9 ) = FP*SP*RZ*0.5 DH(1,10) = TP*SZ*FM*0.5 DH(2,10) = TP*FS*RM*0.5 DH(3,10) = FP*SZ*RM*0.5 DH(1,11) = TP*SM*FR*0.5 DH(2,11) = TP*FM*RZ*0.5 DH(3,11) = FP*SM*RZ*0.5 DH(1,12) = TP*SZ*FP*0.5 DH(2,12) = TP*FS*RP*0.5 DH(3,12) = FP*SZ*RP*0.5 DH(1,13) = TM*SP*FR*0.5 DH(2,13) = TM*FP*RZ*0.5 DH(3,13) = FM*SP*RZ*0.5 DH(1,14) = TM*SZ*FM*0.5 DH(2,14) = TM*FS*RM*0.5 DH(3,14) = FM*SZ*RM*0.5 DH(1,15) = TM*SM*FR*0.5 DH(2,15) = TM*FM*RZ*0.5 DH(3,15) = FM*SM*RZ*0.5 DH(1,16) = TM*SZ*FP*0.5 DH(2,16) = TM*FS*RP*0.5 DH(3,16) = FM*SZ*RP*0.5 DH(1,17) = TZ*SP*FP*0.5 DH(2,17) = TZ*FP*RP*0.5 DH(3,17) = FT*SP*RP*0.5 DH(1,18) = TZ*SP*FM*0.5 DH(2,18) = TZ*FP*RM*0.5 DH(3,18) = FT*SP*RM*0.5 DH(1,19) = TZ*SM*FM*0.5 DH(2,19) = TZ*FM*RM*0.5 DH(3,19) = FT*SM*RM*0.5 DH(1,20) = TZ*SM*FP*0.5 DH(2,20) = TZ*FM*RP*0.5 DH(3,20) = FT*SM*RP*0.5 C LOOP OVER TWELVE EDGES DO 20 K = 9,20 IF ( LNODE(K) .EQ. 0 ) THEN C ZERO EDGE BUBBLE DERIVATIVES DH(1,K) = 0.0 DH(2,K) = 0.0 DH(3,K) = 0.0 ELSE C ENRICH DERIVATIVES AT TWO ENDS OF EDGE DH1 = DH(1,K) DH2 = DH(2,K) DH3 = DH(3,K) K1 = I1(K) K2 = I2(K) DH(1,K1) = DH(1,K1) - DH1 DH(2,K1) = DH(2,K1) - DH2 DH(3,K1) = DH(3,K1) - DH3 DH(1,K2) = DH(1,K2) - DH1 DH(2,K2) = DH(2,K2) - DH2 DH(3,K2) = DH(3,K2) - DH3 DH(1,K) = DH1 + DH1 DH(2,K) = DH2 + DH2 DH(3,K) = DH3 + DH3 ENDIF 20 CONTINUE RETURN END SUBROUTINE DER2C1L (R,A,D2H) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * C 2ND DERIVATIVES OF CUBIC HERMITE IN UNIT COORDINATES C (A C1 ELEMENT, SEE SHPC1L & DERC1L) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION D2H(4) C A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R C R = LOCAL COORDINATE OF POINT R=0 R=1 C D2H = SECOND DERIVATIVES OF H C D()/DX = D()/DR DR/DX = 1/A * D()/DR C DOF ARE VALUE, SLOPE AT EACH END (WRT X) D2H(1) = 6.*(R+R - 1.)/A**2 D2H(2) = (-4. + 6.*R)/A D2H(3) = 6.*(1. - R - R)/A**2 D2H(4) = (6.*R - 2.)/A RETURN END SUBROUTINE DER2C2L (R,A,D2H) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C 2ND DERIVATIVES OF FIFTH ORDER HERMITE ELEMENT IN UNIT C COORDINATES ( A C2 ELEMENT, SEE SUBR SHPC2L & DERC2L ) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION D2H(6) C A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R C R = LOCAL COORDINATE OF POINT R=0 R=1 C D2H = SECOND DERIVATIVES OF H C D()/DX = D()/DR DR/DX = 1/A * D()/DR C DOF ARE VALUE, SLOPE, 2ND DERIV AT EACH END (WRT X) R3 = R*R*R D2H(1) = 30.*(6.*R**2 - 2.*R - 4.*R3)/A**2 D2H(2) = (-36.*R + 96.*R*R - 60.*R3)/A D2H(3) = 0.5*(2. - 18.*R + 36.*R*R - 20.*R3) D2H(4) = 30.*(2.*R - 6.*R*R + 4.*R3)/A**2 D2H(5) = (84.*R*R - 60.*R3 - 24.*R)/A D2H(6) = 0.5*(6.*R - 24.*R*R + 20.*R3) RETURN END SUBROUTINE DER2C3L (R,A,D2H) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * C 2ND DERIV FOR SEVENTH ORDER HERMITE LINE ELEMENT C ( A C3 ELEMENT, SEE SHPC3L & DERC3L ) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION D2H(8) C A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R C R = LOCAL COORDINATE OF POINT R=0 R=1 C D()/DX = D()/DR DR/DX = 1/A * D()/DR C D2H = SECOND DERIVATIVES OF SHAPE FUNCTIONS ARRAY C DOF ARE VALUE, SLOPE, 2ND, 3RD DERIV AT EACH END (WRT X) R2 = R*R R3 = R2*R R4 = R2*R2 D2H(1) = 140.*(6.*R2*R3 - 15.*R4 + 12.*R3 - 3.*R2)/A**2 D2H(2) = -(240.*R2 - 900.*R3 + 1080.*R4 - 420.*R3*R2)/A D2H(3) = (1. - 60.*R2 + 200.*R3 - 225.*R4 + 84.*R2*R3) D2H(4) = (R - 8.*R2 + 20.*R3 - 20.*R4 + 7.*R2*R3)*A D2H(5) = 140.*(3.*R2 - 12.*R3 + 15.*R4 - 6.*R2*R3)/A**2 D2H(6) = (420.*R3*R2 - 1020.*R4 + 780.*R3 - 180.*R2)/A D2H(7) = (30.*R2 - 140.*R3 + 195.*R4 - 84.*R3*R2) D2H(8) = (7.*R3*R2 - 15.*R4 + 10.*R3 - 2.*R2)*A RETURN END SUBROUTINE DER2CU (B, A, D2H) C * * * * * * * * * * * * * * * * * * * * * * * * C SECOND DERIVATIVES OF SHAPE FUNCTIONS FOR 1-D C CUBIC HERMITE ELEMENT (A C1 ELEMENT) C * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION D2H(4) C A = LENGTH OF ELEMENT (SEE SUBR SHPCU) C B = COORDINATE OF POINT C D2H = SECOND DERIVATIVES OF H D2H(1) = (12.*B - 6.)/A/A D2H(2) = (6.*B - 4.)/A D2H(3) = (6. - 12.*B)/A/A D2H(4) = (6.*B - 2.)/A RETURN END SUBROUTINE DER2L (R,DH) C * * * * * * * * * * * * * * * * * * * * * * * * * C DERIVATIVES OF A 2 NODE LINE ELEMENT C * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION DH(2) C R IS UNIT COORD. R=-1 1------------2 R=1 DH(1) = -0.5 DH(2) = 0.5 RETURN END SUBROUTINE DER3L (X, DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * C FIND LOCAL DERIVATIVES FOR A 3 NODE LINE ELEMENT C IN NATURAL COORDINATES C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DH(3) C DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS (SHP3L) C X = LOCAL COORDINATE OF POINT, -1 TO +1 C LOCAL NODE COORD. ARE -1,0,+1. 1----2----3 DH(1) = X - 0.5 DH(2) = -2.*X DH(3) = X + 0.5 RETURN END SUBROUTINE DER3T (S, T, DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * C LOCAL DERIVATIVES OF A THREE NODE UNIT TRIANGLE C SEE SUBROUTINE SHP3T C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DH(2,3) C S,T = LOCAL COORDINATES OF THE POINT C DH(1,K) = DH(K)/DS C DH(2,K) = DH(K)/DT C NODAL COORDS ARE : 1-(0,0) 2-(1,0) 3-(0,1) DH(1,1) = -1. DH(1,2) = 1. DH(1,3) = 0.0 DH(2,1) = -1. DH(2,2) = 0.0 DH(2,3) = 1. RETURN END SUBROUTINE DER4Q (R, S, DELTA) C * * * * * * * * * * * * * * * * * * * * * * * * * * C LOCAL DERIVATIVES OF THE SHAPE FUNCTIONS FOR AN C ISOPARAMETRIC QUADRILATERAL WITH FOUR NODES C SEE SHP4Q C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DELTA(2,4) C DELTA(1,I) = DH/DR C DELTA(2,I) = DH/DS C H = LOCAL INTERPOLATION FUNCTIONS C (R,S) = A POINT IN THE LOCAL COORDINATES C HERE D(H(I))/DR = 0.25*R(I)*(1+S*S(I)), ETC. RP = 1. + R RM = 1. - R SP = 1. + S SM = 1. - S DELTA(1,1) = -0.25*SM DELTA(1,2) = 0.25*SM DELTA(1,3) = 0.25*SP DELTA(1,4) = -0.25*SP DELTA(2,1) = -0.25*RM DELTA(2,2) = -0.25*RP DELTA(2,3) = 0.25*RP DELTA(2,4) = 0.25*RM RETURN END SUBROUTINE DER6T (S,T,DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * C LOCAL DERIVATIVES FOR A SIX NODE UNIT TRIANGLE C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DH(2,6) C S,T = LOCAL COORDINATES, SEE SHP6T C DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS C DH(1,K) = DH(K)/DS, DH(2,K)=DH(K)/DT C NODAL COORDS : 1-(0,0) 2-(1,0) 3-(0,1) C 4-(0.5,0) 5-(0.5,0.5) 6-(0,0.5) DH(1,1) = -3. + 4.*S + 4.*T DH(1,2) = -1. + 4.*S DH(1,3) = 0.0 DH(1,4) = 4. - 8.*S - 4.*T DH(1,5) = 4.*T DH(1,6) = -4.*T DH(2,1) = -3. + 4.*S + 4.*T DH(2,2) = 0.0 DH(2,3) = -1. + 4.*T DH(2,4) = -4.*S DH(2,5) = 4.*S DH(2,6) = 4. -4.*S - 8.*T RETURN END SUBROUTINE DER8H (R,S,T,DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * C LOCAL DERIVATIVES FOR EIGHT NODE HEXAHEDRON C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DH(3,8) C R,S,T = LOCAL COORDINATES OF THE POINT C DH(1,K)=DH/DR, DH(2,K)=DH/DS, DH(3,K)=DH/DT C H = ELEMENT SHAPE FUNCTIONS, SEE SHP8H RP = 1. + R RM = 1. - R SP = 1. + S SM = 1. - S TP = 1. + T TM = 1. - T DH(1,1) = 0.125*SP*TP DH(2,1) = 0.125*RP*TP DH(3,1) = 0.125*RP*SP DH(1,2) = 0.125*SM*TP DH(2,2) = -0.125*RP*TP DH(3,2) = 0.125*RP*SM DH(1,3) = 0.125*SM*TM DH(2,3) = -0.125*RP*TM DH(3,3) = -0.125*RP*SM DH(1,4) = 0.125*SP*TM DH(2,4) = 0.125*RP*TM DH(3,4) = -0.125*RP*SP DH(1,5) = -0.125*SP*TP DH(2,5) = 0.125*RM*TP DH(3,5) = 0.125*RM*SP DH(1,6) = -0.125*SM*TP DH(2,6) = -0.125*RM*TP DH(3,6) = 0.125*RM*SM DH(1,7) = -0.125*SM*TM DH(2,7) = -0.125*RM*TM DH(3,7) = -0.125*RM*SM DH(1,8) = -0.125*SP*TM DH(2,8) = 0.125*RM*TM DH(3,8) = -0.125*RM*SP RETURN END SUBROUTINE DER8Q (S,T,DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * C FIND LOCAL DERIVATIVES OF SHAPE FUNCTIONS FOR AN C EIGHT NODE ISOPARAMETRIC QUADRILATERAL ELEMENT C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DH(2,8) C S,T = LOCAL COORDINATES OF POINT, SEE SHP8Q C DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS, H C DH(1,J) = DH(J)/DS, DH(2,J) = DH(J)/DT C H = SHAPE FUNCTION ARRAY SP = 1. + S SM = 1. - S TP = 1. + T TM = 1. - T DH(1,1) = -0.25*TM*( SM + SM + TM - 3. ) DH(2,1) = -0.25*SM*( TM + SM + TM - 3. ) DH(1,2) = 0.25*TM*( SP + SP + TM - 3. ) DH(2,2) = -0.25*SP*( TM + SP + TM - 3. ) DH(1,3) = 0.25*TP*( SP + SP + TP - 3. ) DH(2,3) = 0.25*SP*( TP + SP + TP - 3. ) DH(1,4) = -0.25*TP*( SM + SM + TP - 3. ) DH(2,4) = 0.25*SM*( TP + SM + TP - 3. ) DH(1,5) = -S*TM DH(2,5) = -0.5*( 1. - S*S ) DH(1,6) = 0.5*( 1. - T*T ) DH(2,6) = -T*SP DH(1,7) = -S*TP DH(2,7) = 0.5*( 1. - S*S ) DH(1,8) = -0.5*( 1. - T*T ) DH(2,8) = -T*SM RETURN END SUBROUTINE DER9Q ( R, S, DH ) C * * * * * * * * * * * * * * * * * * * * * * * C LOCAL DERIVATIVES FOR 9-NODED QUAD C * * * * * * * * * * * * * * * * * * * * * * * C SEE SHP9Q FOR TOPOLOGY CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DH(2,9) RM = R - 1.D0 SM = S - 1.D0 RP = R + 1.D0 SP = S + 1.D0 S2P1 = S + S + 1.D0 S2M1 = S + S - 1.D0 R2P1 = R + R + 1.D0 R2M1 = R + R - 1.D0 DH(1,1) = 0.25D0 * S * SM * R2M1 DH(1,2) = 0.25D0 * S * SM * R2P1 DH(1,3) = 0.25D0 * S * SP * R2P1 DH(1,4) = 0.25D0 * S * SP * R2M1 DH(1,5) = -S * SM * R DH(1,6) = -0.5D0 * SP * SM * R2P1 DH(1,7) = -S * SP * R DH(1,8) = -0.5D0 * SP * SM * R2M1 DH(1,9) = 2.D0 * SP * SM * R DH(2,1) = 0.25D0 * S2M1 * R * RM DH(2,2) = 0.25D0 * S2M1 * R * RP DH(2,3) = 0.25D0 * S2P1 * R * RP DH(2,4) = 0.25D0 * S2P1 * R * RM DH(2,5) = -0.5D0 * S2M1 * RP * RM DH(2,6) = -S * R * RP DH(2,7) = -0.5D0 * S2P1 * RP * RM DH(2,8) = -S * R * RM DH(2,9) = 2.D0 * S * RP * RM RETURN END SUBROUTINE DERC0L (R, A, DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * C FIRST DERIVATIVE OF C0 LINE ELEMENT IN UNIT COORDINATES C (SEE SHPC0L) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION DH(2) C A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R C R = LOCAL COORDINATE OF POINT R=0 R=1 C DH = FIRST PHYSICAL DERIVATIVES OF H C D()/DX = D()/DR DR/DX = 1/A * D()/DR C DOF ARE NODAL VALUES ONLY DH(1) = -1./A DH(2) = 1./A RETURN END SUBROUTINE DERC1L (R,A,DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * C FIRST DERIVATIVES OF CUBIC HERMITE IN UNIT COORDINATES C (A C1 ELEMENT, SEE SHPC1L) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION DH(4) C A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R C R = LOCAL COORDINATE OF POINT R=0 R=1 C DH = FIRST PHYSICAL DERIVATIVES OF H C D()/DX = D()/DR DR/DX = 1/A * D()/DR C DOF ARE FUNCTION AND SLOPE, D/DX, AT EACH NODE DH(1) = 6.0*(R*R - R)/A DH(2) = 1.0 - 4.0*R + 3.0*R*R DH(3) = 6.0*(R - R*R)/A DH(4) = 3.0*R*R - 2.0*R RETURN END SUBROUTINE DERC2L (R,A,DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C FIRST DERIVATIVES OF FIFTH ORDER HERMITE ELEMENT IN UNIT C COORDINATES ( A C2 ELEMENT, SEE SURR SHPC2L ) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION DH(6) C A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R C R = LOCAL COORDINATE OF POINT R=0 R=1 C DH = FIRST PHYSICAL DERIVATIVES OF H C D()/DX = D()/DR DR/DX = 1/A * D()/DR C DOF ARE VALUE, SLOPE, CURVATURE AT EACH NODE (WRT X) R3 = R*R*R DH(1) = 30.*(2.*R3 - R*R - R3*R)/A DH(2) = 1. - 18.*R*R + 32.*R3 - 15.*R3*R DH(3) = 0.5*(2.*R - 9.*R*R + 12.*R3 - 5.*R3*R)*A DH(4) = 30.*R*R*(1. - 2.*R + R*R)/A DH(5) = 28.*R3 - 15.*R3*R - 12.*R*R DH(6) = 0.5*R*R*(3. - 8.*R + 5.*R*R)*A RETURN END SUBROUTINE DERC3L (R,A,DH) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * C FIRST DERIVATIVES FOR 7TH ORDER HERMITE LINE ELEMENT C ( A C3 ELEMENT, SEE SHPC3L ) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * DIMENSION DH(8) C A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R C R = LOCAL COORDINATE OF POINT R=0 R=1 C DH = FIRST PHYSICAL DERIVATIVES OF H C D()/DX = D()/DR DR/DX = 1/A * D()/DR C DOF ARE VALUE, SLOPE, 2ND, 3RD DERIV AT EACH END (WRT X) R2 = R*R R3 = R2*R R4 = R2*R2 DH(1) = 140.*R3*(R3 - 3.*R2 + 3.*R - 1.)/A DH(2) = 1. - R3*(80. - 225.*R + 216.*R2 - 70.*R3) DH(3) = R*(1. - 20.*R2 + 50.*R3 - 45.*R4 + 14.*R2*R3)*A DH(4) = R2*(3. - 16.*R + 30.*R2 - 24.*R3 + 7.*R4)*A*A/6. DH(5) = 140.*R3*(1. - 3.*R + 3.*R2 - R3)/A DH(6) = R3*(70.*R3 - 204.*R2 + 195.*R - 60.) DH(7) = R3*(10. - 35.*R + 39.*R2 - 14.*R3)*A DH(8) = R3*(7.*R3 - 18.*R2 + 15.*R - 4.)*A*A/6. RETURN END SUBROUTINE DERCU (B, A, DH) C * * * * * * * * * * * * * * * * * * * * * * * * C FIRST DERIVATIVES OF SHAPE FUNCTIONS FOR 1-D C CUBIC HERMITE ELEMENT (A C1 ELEMENT) C * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DH(4) C A = LENGTH OF ELEMENT (SEE SUBR SHPCU) C B = COORDINATE OF POINT C DH = FIRST DERIVATIVES OF H DH(1) = 6.*(B*B - B)/A DH(2) = 1. - 4.*B + 3.*B*B DH(3) = 6.*(B - B*B)/A DH(4) = 3.*B*B - 2.*B RETURN END SUBROUTINE DERHQL (NODEDG, LOCATE, NEDGE, LEDGES, NSPACE, & RST, DERIV) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SHAPE FUNCTION DERIVATIVES FOR GENERAL SERENDIPITY C LINE, QUAD, OR OR HEXAHEDRON WITH AN C ARBITRARY NUMBER OF NODES ON EACH EDGE C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) PARAMETER ( MAXDEG = 20 ) DIMENSION RST(3), BLKCRD(3,8), POLI2(3), CDRFN(3), & FARSID(3), DERIV(3), CRDEDG(3,MAXDEG+1), & NODEDG(12), NEATC(3,8), NODEOP(2,12), & NODATC(3), LOCAL(12) DATA BLKCRD &/ -1.,-1.,-1., 1.,-1.,-1., 1.,1.,-1., -1.,1.,-1., & -1.,-1., 1., 1.,-1., 1., 1.,1., 1., -1.,1., 1./ DATA NEATC / 1,4,9, 1,2,10, 3,2,11, 3,4,12, & 5,8,9, 5,6,10, 7,6,11, 7,8,12 / DATA NODEOP / 7,8, 8,5, 5,6, 6,7, & 3,4, 4,1, 1,2, 2,3, & 3,7, 4,8, 1,5, 2,6 / DATA LOCAL / -1, -2, 1, 2, -1, -2, 1, 2, -3, -3, -3, -3 / C BLKCRD = BLOCK CORNER LOCAL COORDINATES C CRDEDG = LOCAL COORDINATES OF SIDE NODES JOINING CORNER C FARSID = FAR SIDE LOCAL COORDINATE C LEDGES = NUMBER OF ELEMENT EDGES, 1, 4, OR 12 C LOCAL = LOCAL COORDINATE PARALLEL TO EACH EDGE C LOCATE = POSITION NUMBER ON EDGE, 0 IF CORNER C MAXDEG = MAXIMUM PLOYNOMIAL DEGREE ON ANY SIDE C NEATC = THE 1, 2, OR 3 EDGES AT A CORNER C NEDGE = EDGE NUMBER OR CORNER NUMBER OF THE NODE COMPUTED C NODATC = NUMBER OF SIDE NODES JOINING A CORNER C NODEDG = NUMBER ON NODES ON 1,4, OR 12 EDGES C NODEOP = 2 DIAGONALLY OPPOSITE NODES FOR EACH EDGE C NSPACE = NUMBER OF SPATIAL DIMENSIONS C RST = LOCAL COORDINATES FOR EVALUATION C VALUE = SHAPE FUNCTION VALUE (RETURNED) C C VALUE = A(R,S,T)*( P1(R) + P2(S) + P3(T) + CONSTANT ) C DERIV = DA(R,S,T)*( P1(R) + P2(S) + P3(T) + CONSTANT ) C + A(R,S,T)*( DP1(R) + DP2(S) + DP3(T) ) C C REF: G. ZAVARISE, ET AL, "AN ALGORITHM FOR GENERATION OF SHAPE C FUNCTIONS IN SERENDIPITY ELEMENTS", ENG COMP,8,19-31,1991 C C T: S C8 *---E7----* C7 T: S 8---15----7 C : / /. /: : / /. /: C :/ / . / : :/ / 22 / : C *---R / E12 / E11 *---R / . / 20 C E8 . E6 : 16 21 / : C / . / : / . / : C / C4*.../.E3..* C3 / 4.13/.12..3 C / . / / / . / / C C5 *--E5-----* C6 / 5---------6 / C : . : / : . : 11 C : E4 : E2 : 14 19 / C E9 . E10 / 17 . : 10 C : . : / : . 18 / C :. :/ :. :/ C C1 *---E1----* C2 1----9----2 C CORNER NODE & EDGE NUMBERS. 22 NODES: CORNERS, THEN BY EDGES. C CCW IF |T|=1, ELSE IN POSITIVE T. C === 3-D FORM === C C C4 *---E3----* C3 4----8----3 C : . : :S : : C : : : : : C E4 E2 : 9 : C : : *---R : : C : : : : C C1 *---E1----* C2 1--5-6-7--2 C CORNER NODE & EDGE NUMBERS. 9 NODES: CORNERS, THEN BY EDGE ORDER. C === 2-D FORM === C C C1 *---E1----* C2 1--2-3-4--5 C CORNER NODE & EDGE NUMBERS. 9 NODES NUMBERED BY EDGE ORDER. C === 1-D FORM === POLI1 = 1. IF ( LOCATE .EQ. 0 ) THEN C C SHAPE FUNCTION FOR CORNER NODES C DO 100 ICORD = 1,NSPACE POLI1 = POLI1*(RST(ICORD) + BLKCRD(ICORD,NEDGE)) & /(2*BLKCRD(ICORD,NEDGE)) 100 CONTINUE CPNUL = 1. POLI2(1) = 0. POLI2(2) = 0. POLI2(3) = 0. DO 200 ICORD = 1,NSPACE NSIDE = NEATC(ICORD,NEDGE) NODATC(ICORD) = NODEDG(NSIDE) - 2 IF ( NODATC(ICORD) .GT. 0 ) THEN IF ( NODATC(ICORD) .GT. MAXDEG ) STOP 'MAXDEG, DERSHAFN' CPNUL = CPNUL - 1. POLI2(ICORD) = 1. FARSID(ICORD) = 2./(NODEDG(NSIDE) - 1) DO 300 INODE = 1,NODATC(ICORD) CRDEDG(ICORD,INODE) = -1. + FARSID(ICORD)*INODE POLI2(ICORD) = POLI2(ICORD)*(RST(ICORD) & - CRDEDG(ICORD,INODE))/(BLKCRD(ICORD,NEDGE) & - CRDEDG(ICORD,INODE)) 300 CONTINUE ENDIF 200 CONTINUE C VALUE = POLI1*(POLI2(1) + POLI2(2) + POLI2(3) + CPNUL) ELSE C C SHAPE FUNCTION FOR EDGE NODES C NOPV1 = NODEOP(1,NEDGE) NOPV2 = NODEOP(2,NEDGE) ISRFN = ABS(LOCAL(NEDGE)) FARSID(1) = 2./(NODEDG(NEDGE) - 1) CDRFN(1) = -BLKCRD(1,NOPV1) CDRFN(2) = -BLKCRD(2,NOPV1) CDRFN(3) = -BLKCRD(3,NOPV1) CDRFN(ISRFN) = (1. - FARSID(1)*LOCATE)*LOCAL(NEDGE)/ISRFN DO 400 ICORD = 1,NSPACE POLI1 = POLI1*(RST(ICORD) - BLKCRD(ICORD,NOPV1)) & /(CDRFN(ICORD) - BLKCRD(ICORD,NOPV1)) 400 CONTINUE PLAN2 = (RST(ISRFN) - BLKCRD(ISRFN,NOPV2)) & /(CDRFN(ISRFN) - BLKCRD(ISRFN,NOPV2)) POLI3 = 1. NODATC(1) = NODEDG(NEDGE) - 2 IF ( NODATC(1) .GT. 0 ) THEN IF ( NODATC(1) .GT. MAXDEG ) STOP 'MAXDEG, DERSHAFN' DO 500 INODE = 1,NODATC(1) CRDEDG(1,INODE) = -1. + FARSID(1)*INODE IF ( ABS(CRDEDG(1,INODE) - CDRFN(ISRFN)) .GT. 0.0001) & THEN POLI3 = POLI3*(RST(ISRFN) - CRDEDG(1,INODE)) & /(CDRFN(ISRFN) - CRDEDG(1,INODE)) ENDIF 500 CONTINUE ENDIF C VALUE = POLI1*PLAN2*POLI3 ENDIF C C DERIVATIVES OF SHAPE FUNCTIONS C DO 600 ICOR1 = 1,NSPACE IF ( LOCATE .EQ. 0 ) THEN C C DERIVATIVES FOR CORNER NODES C DPOL1 = POLI2(1) + POLI2(2) + POLI2(3) + CPNUL DO 700 ICOR2 = 1,NSPACE IF ( ICOR2 .NE. ICOR1 ) THEN DPOL1 = DPOL1*(RST(ICOR2) + BLKCRD(ICOR2,NEDGE)) & /(2*BLKCRD(ICOR2,NEDGE)) ELSE DPOL1 = DPOL1/(2*BLKCRD(ICOR2,NEDGE)) ENDIF 700 CONTINUE DPOL2 = 0. DO 800 INOD1 = 1,NODATC(ICOR1) DETP2 = 1. DO 900 INOD2 = 1,NODATC(ICOR1) IF ( INOD2 .NE. INOD1 ) THEN DETP2 = DETP2*(RST(ICOR1) - CRDEDG(ICOR1,INOD2)) & /(BLKCRD(ICOR1,NEDGE) - CRDEDG(ICOR1,INOD2)) ELSE DETP2 = DETP2/(BLKCRD(ICOR1,NEDGE) & - CRDEDG(ICOR1,INOD2)) ENDIF 900 CONTINUE DPOL2 = DPOL2 + DETP2 800 CONTINUE DPOL2 = DPOL2*POLI1 DERIV(ICOR1) = DPOL1 + DPOL2 ELSE C C DERIVATIVES FOR EDGE NODES C DPOL1 = POLI3*PLAN2 DO 1000 ICOR2 = 1,NSPACE IF ( ICOR2 .NE. ICOR1 ) THEN DPOL1 = DPOL1*(RST(ICOR2) - BLKCRD(ICOR2,NOPV1)) & /(CDRFN(ICOR2) - BLKCRD(ICOR2,NOPV1)) ELSE DPOL1 = DPOL1/(CDRFN(ICOR2) - BLKCRD(ICOR2,NOPV1)) ENDIF 1000 CONTINUE DPLA2 = 0. DPOL3 = 0. IF ( ICOR1 .EQ. ISRFN ) THEN DPLA2 = POLI1*POLI3/(CDRFN(ISRFN) - BLKCRD(ISRFN,NOPV2)) DO 1100 INOD1 = 1,NODATC(1) IF ( ABS(CRDEDG(1,INOD1) - CDRFN(ISRFN)) .GT. 0.0001) & THEN DETP3 = 1. DO 1200 INOD2 = 1,NODATC(1) IF ( ABS(CRDEDG(1,INOD2) - CDRFN(ISRFN)) .GT. & 0.0001 ) THEN IF ( INOD2 .NE. INOD1 ) THEN DETP3 = DETP3*(RST(ISRFN) - CRDEDG(1,INOD2)) & /(CDRFN(ISRFN) - CRDEDG(1,INOD2)) ELSE DETP3 = DETP3/(CDRFN(ISRFN) - CRDEDG(1,INOD2)) ENDIF ENDIF 1200 CONTINUE DPOL3 = DPOL3 + DETP3 ENDIF 1100 CONTINUE DPOL3 = DPOL3*POLI1*PLAN2 ENDIF DERIV(ICOR1) = DPOL1 + DPLA2 + DPOL3 ENDIF 600 CONTINUE RETURN END SUBROUTINE DERIV (PT, DLH, N, NSPACE, LSHAPE, NG, LNODE) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C EVALUATE C0 INTERPOLATION LOCAL DERIVATIVES C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) cb DIMENSION DLH(NSPACE,N), PT(NSPACE), LNODE(N) DIMENSION DLH(NSPACE,N*NG), PT(NSPACE), LNODE(N) C DLH = LOCAL DERIVATIVES AT PT C LNODE = TOPOLOGY LIST, IF VARIABLE C LSHAPE = 1-LINE, 2-TRI, 3-QUAD, 4-HEX, 5-TET, 6-WEDGE C N = NUMBER OF NODES PER ELEMENT C NG = NUMBER OF DEGREES OF FREEDOM PER NODE C NSPACE = NO OF SPATIAL DIMENSIONS C PT = LOCAL COORD OF A POINT C C BRANCH ON SHAPE, THEN NUMBER OF NODES IF ( LSHAPE .LE. 1 ) THEN C--> 1-D ELEMENTS IF ( N .EQ. 2 ) CALL DER2L (PT(1),DLH) c IF ( N .EQ. 3 ) CALL DER3L (PT(1),DLH) RETURN ELSEIF ( LSHAPE .EQ. 2 ) THEN C--> TRIANGULAR 2-D ELEMENTS IF ( N .EQ. 3 ) CALL DER3T (PT(1),PT(2),DLH) C IF ( N .EQ. 4 ) CALL DER4T (PT(1),PT(2),DLH) IF ( N .EQ. 6 ) CALL DER6T (PT(1),PT(2),DLH) C IF ( N .EQ. 7 ) CALL DER7T (PT(1),PT(2),DLH) C IF ( N .EQ. 10 ) CALL DER10T (PT(1),PT(2),DLH) C IF ( N .EQ. 15 ) CALL DER15T (PT(1),PT(2),DLH) RETURN ELSEIF ( LSHAPE .EQ. 3 ) THEN C--> QUADRILATERAL 2-D ELEMENTS IF ( N .EQ. 4 ) CALL DER4Q (PT(1),PT(2),DLH) IF ( N .EQ. 8 ) CALL DER8Q (PT(1),PT(2),DLH) IF ( N .EQ. 9 ) CALL DER9Q (PT(1),PT(2),DLH) C IF ( N .EQ. 12 ) CALL DER412 (PT(1),PT(2),DLH,LNODE) C IF ( N .EQ. 16 ) CALL DER16Q (PT(1),PT(2),DLH) C IF ( N .EQ. 17 ) CALL DER17Q (PT(1),PT(2),DLH) C IF ( N .EQ. 25 ) CALL DER25Q (PT(1),PT(2),DLH) RETURN ELSEIF ( LSHAPE .EQ. 4 ) THEN C--> HEXAHEDRA 3-D ELEMENTS IF ( N .EQ. 8 ) CALL DER8H (PT(1),PT(2),PT(3),DLH) c IF ( N .EQ. 20 ) CALL DER208 (PT(1),PT(2),PT(3),DLH,LNODE) C IF ( N .EQ. 27 ) CALL DER27H (PT(1),PT(2),PT(3),DLH) C IF ( N .EQ. 32 ) CALL DER32H (PT(1),PT(2),PT(3),DLH) RETURN ELSEIF ( LSHAPE .EQ. 5 ) THEN C--> TETRAHEDRA 3-D ELEMENTS (PYRAMIDS) c IF ( N .EQ. 4 ) CALL DER4P (PT(1),PT(2),PT(3),DLH) C IF ( N .EQ. 10 ) CALL DER10P (PT(1),PT(2),PT(3),DLH) C IF ( N .EQ. 21 ) CALL DER21P (PT(1),PT(2),PT(3),DLH) RETURN ELSEIF ( LSHAPE .EQ. 6 ) THEN C--> WEDGE 3-D ELEMENTS STOP 'NO WEDGE IN SHAPE' C IF ( N .EQ. 6 ) CALL DER6W (PT(1),PT(2),PT(3),DLH) C IF ( N .EQ. 15 ) CALL DER15W (PT(1),PT(2),PT(3),DLH) C RETURN ELSEIF ( LSHAPE .EQ. 7 ) THEN C--> USER SUPPLIED ELEMENT C CALL DERUSR (PT(1),PT(2),PT(3),DLH,LNODE) STOP 'NO USER ELEMENT IN DERIV' ELSEIF ( LSHAPE .GT. 7 ) THEN C--> UNSUPPORTED OPTION STOP 'UNSUPPORTED ELEMENT IN DERIV' ENDIF RETURN END SUBROUTINE DERSHPH (NODEDG, LOCATE, NEDGE, LEDGES, NSPACE, & RST, VALUE, DERIV) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C SHAPE FUNCTIONS AND DERIVATIVES FOR GENERAL SERENDIPITY C LINE, QUAD, OR OR HEXAHEDRON WITH AN C ARBITRARY NUMBER OF NODES ON EACH EDGE C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) PARAMETER ( MAXDEG = 20 ) DIMENSION RST(3), BLKCRD(3,8), POLI2(3), CDRFN(3), & FARSID(3), DERIV(3), CRDEDG(3,MAXDEG+1), & NODEDG(12), NEATC(3,8), NODEOP(2,12), & NODATC(3), LOCAL(12) DATA BLKCRD &/ -1.,-1.,-1., 1.,-1.,-1., 1.,1.,-1., -1.,1.,-1., & -1.,-1., 1., 1.,-1., 1., 1.,1., 1., -1.,1., 1./ DATA NEATC / 1,4,9, 1,2,10, 3,2,11, 3,4,12, & 5,8,9, 5,6,10, 7,6,11, 7,8,12 / DATA NODEOP / 7,8, 8,5, 5,6, 6,7, & 3,4, 4,1, 1,2, 2,3, & 3,7, 4,8, 1,5, 2,6 / DATA LOCAL / -1, -2, 1, 2, -1, -2, 1, 2, -3, -3, -3, -3 / C BLKCRD = BLOCK CORNER LOCAL COORDINATES C CRDEDG = LOCAL COORDINATES OF SIDE NODES JOINING CORNER C DERIV = SHAPE FUNCTION DERIVATIVES (RETURNED) C FARSID = FAR SIDE LOCAL COORDINATE C LEDGES = NUMBER OF ELEMENT EDGES, 1, 4, OR 12 C LOCAL = LOCAL COORDINATE PARALLEL TO EACH EDGE C LOCATE = POSITION NUMBER ON EDGE, 0 IF CORNER C MAXDEG = MAXIMUM PLOYNOMIAL DEGREE ON ANY SIDE C NEATC = THE 1, 2, OR 3 EDGES AT A CORNER C NEDGE = EDGE NUMBER OR CORNER NUMBER OF THE NODE COMPUTED C NODATC = NUMBER OF SIDE NODES JOINING A CORNER C NODEDG = NUMBER ON NODES ON 1,4, OR 12 EDGES C NODEOP = 2 DIAGONALLY OPPOSITE NODES FOR EACH EDGE C NSPACE = NUMBER OF SPATIAL DIMENSIONS C RST = LOCAL COORDINATES FOR EVALUATION C VALUE = SHAPE FUNCTION VALUE (RETURNED) C C VALUE = A(R,S,T)*( P1(R) + P2(S) + P3(T) + CONSTANT ) C DERIV = DA(R,S,T)*( P1(R) + P2(S) + P3(T) + CONSTANT ) C + A(R,S,T)*( DP1(R) + DP2(S) + DP3(T) ) C C REF: G. ZAVARISE, ET AL, "AN ALGORITHM FOR GENERATION OF SHAPE C FUNCTIONS IN SERENDIPITY ELEMENTS, I.J.N.M.E." C C T: S C8 *---E7----* C7 T: S 8---15----7 C : / /. /: : / /. /: C :/ / . / : :/ / 22 / : C *---R / E12 / E11 *---R / . / 20 C E8 . E6 : 16 21 / : C / . / : / . / : C / C4*.../.E3..* C3 / 4.13/.12..3 C / . / / / . / / C C5 *--E5-----* C6 / 5---------6 / C : . : / : . : 11 C : E4 : E2 : 14 19 / C E9 . E10 / 17 . : 10 C : . : / : . 18 / C :. :/ :. :/ C C1 *---E1----* C2 1----9----2 C CORNER NODE & EDGE NUMBERS. 22 NODES: CORNERS, THEN BY EDGES. C CCW IF |T|=1, ELSE IN POSITIVE T. C === 3-D FORM === C C C4 *---E3----* C3 4----8----3 C : . : :S : : C : : : : : C E4 E2 : 9 : C : : *---R : : C : : : : C C1 *---E1----* C2 1--5-6-7--2 C CORNER NODE & EDGE NUMBERS. 9 NODES: CORNERS, THEN BY EDGE ORDER. C === 2-D FORM === C C C1 *---E1----* C2 1--2-3-4--5 C CORNER NODE & EDGE NUMBERS. 9 NODES NUMBERED BY EDGE ORDER. C === 1-D FORM === POLI1 = 1. IF ( LOCATE .EQ. 0 ) THEN C C SHAPE FUNCTION FOR CORNER NODES C DO 100 ICORD = 1,NSPACE POLI1 = POLI1*(RST(ICORD) + BLKCRD(ICORD,NEDGE)) & /(2*BLKCRD(ICORD,NEDGE)) 100 CONTINUE CPNUL = 1. POLI2(1) = 0. POLI2(2) = 0. POLI2(3) = 0. DO 200 ICORD = 1,NSPACE NSIDE = NEATC(ICORD,NEDGE) NODATC(ICORD) = NODEDG(NSIDE) - 2 IF ( NODATC(ICORD) .GT. 0 ) THEN IF ( NODATC(ICORD) .GT. MAXDEG ) STOP 'MAXDEG, DERSHAFN' CPNUL = CPNUL - 1. POLI2(ICORD) = 1. FARSID(ICORD) = 2./(NODEDG(NSIDE) - 1) DO 300 INODE = 1,NODATC(ICORD) CRDEDG(ICORD,INODE) = -1. + FARSID(ICORD)*INODE POLI2(ICORD) = POLI2(ICORD)*(RST(ICORD) & - CRDEDG(ICORD,INODE))/(BLKCRD(ICORD,NEDGE) & - CRDEDG(ICORD,INODE)) 300 CONTINUE ENDIF 200 CONTINUE VALUE = POLI1*(POLI2(1) + POLI2(2) + POLI2(3) + CPNUL) ELSE C C SHAPE FUNCTION FOR EDGE NODES C NOPV1 = NODEOP(1,NEDGE) NOPV2 = NODEOP(2,NEDGE) ISRFN = ABS(LOCAL(NEDGE)) FARSID(1) = 2./(NODEDG(NEDGE) - 1) CDRFN(1) = -BLKCRD(1,NOPV1) CDRFN(2) = -BLKCRD(2,NOPV1) CDRFN(3) = -BLKCRD(3,NOPV1) CDRFN(ISRFN) = (1. - FARSID(1)*LOCATE)*LOCAL(NEDGE)/ISRFN DO 400 ICORD = 1,NSPACE POLI1 = POLI1*(RST(ICORD) - BLKCRD(ICORD,NOPV1)) & /(CDRFN(ICORD) - BLKCRD(ICORD,NOPV1)) 400 CONTINUE PLAN2 = (RST(ISRFN) - BLKCRD(ISRFN,NOPV2)) & /(CDRFN(ISRFN) - BLKCRD(ISRFN,NOPV2)) POLI3 = 1. NODATC(1) = NODEDG(NEDGE) - 2 IF ( NODATC(1) .GT. 0 ) THEN IF ( NODATC(1) .GT. MAXDEG ) STOP 'MAXDEG, DERSHAFN' DO 500 INODE = 1,NODATC(1) CRDEDG(1,INODE) = -1. + FARSID(1)*INODE IF ( ABS(CRDEDG(1,INODE) - CDRFN(ISRFN)) .GT. 0.0001) & THEN POLI3 = POLI3*(RST(ISRFN) - CRDEDG(1,INODE)) & /(CDRFN(ISRFN) - CRDEDG(1,INODE)) ENDIF 500 CONTINUE ENDIF VALUE = POLI1*PLAN2*POLI3 ENDIF C C DERIVATIVES OF SHAPE FUNCTIONS C DO 600 ICOR1 = 1,NSPACE IF ( LOCATE .EQ. 0 ) THEN C C DERIVATIVES FOR CORNER NODES C DPOL1 = POLI2(1) + POLI2(2) + POLI2(3) + CPNUL DO 700 ICOR2 = 1,NSPACE IF ( ICOR2 .NE. ICOR1 ) THEN DPOL1 = DPOL1*(RST(ICOR2) + BLKCRD(ICOR2,NEDGE)) & /(2*BLKCRD(ICOR2,NEDGE)) ELSE DPOL1 = DPOL1/(2*BLKCRD(ICOR2,NEDGE)) ENDIF 700 CONTINUE DPOL2 = 0. DO 800 INOD1 = 1,NODATC(ICOR1) DETP2 = 1. DO 900 INOD2 = 1,NODATC(ICOR1) IF ( INOD2 .NE. INOD1 ) THEN DETP2 = DETP2*(RST(ICOR1) - CRDEDG(ICOR1,INOD2)) & /(BLKCRD(ICOR1,NEDGE) - CRDEDG(ICOR1,INOD2)) ELSE DETP2 = DETP2/(BLKCRD(ICOR1,NEDGE) & - CRDEDG(ICOR1,INOD2)) ENDIF 900 CONTINUE DPOL2 = DPOL2 + DETP2 800 CONTINUE DPOL2 = DPOL2*POLI1 DERIV(ICOR1) = DPOL1 + DPOL2 ELSE C C DERIVATIVES FOR EDGE NODES C DPOL1 = POLI3*PLAN2 DO 1000 ICOR2 = 1,NSPACE IF ( ICOR2 .NE. ICOR1 ) THEN DPOL1 = DPOL1*(RST(ICOR2) - BLKCRD(ICOR2,NOPV1)) & /(CDRFN(ICOR2) - BLKCRD(ICOR2,NOPV1)) ELSE DPOL1 = DPOL1/(CDRFN(ICOR2) - BLKCRD(ICOR2,NOPV1)) ENDIF 1000 CONTINUE DPLA2 = 0. DPOL3 = 0. IF ( ICOR1 .EQ. ISRFN ) THEN DPLA2 = POLI1*POLI3/(CDRFN(ISRFN) - BLKCRD(ISRFN,NOPV2)) DO 1100 INOD1 = 1,NODATC(1) IF ( ABS(CRDEDG(1,INOD1) - CDRFN(ISRFN)) .GT. 0.0001) & THEN DETP3 = 1. DO 1200 INOD2 = 1,NODATC(1) IF ( ABS(CRDEDG(1,INOD2) - CDRFN(ISRFN)) .GT. & 0.0001 ) THEN IF ( INOD2 .NE. INOD1 ) THEN DETP3 = DETP3*(RST(ISRFN) - CRDEDG(1,INOD2)) & /(CDRFN(ISRFN) - CRDEDG(1,INOD2)) ELSE DETP3 = DETP3/(CDRFN(ISRFN) - CRDEDG(1,INOD2)) ENDIF ENDIF 1200 CONTINUE DPOL3 = DPOL3 + DETP3 ENDIF 1100 CONTINUE DPOL3 = DPOL3*POLI1*PLAN2 ENDIF DERIV(ICOR1) = DPOL1 + DPLA2 + DPOL3 ENDIF 600 CONTINUE RETURN END FUNCTION DOT (N, A, B) C * * * * * * * * * * * * * * * * * C DOT PRODUCT OF VECTORS A(N)*B(N) C * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(N), B(N) DOT = 0.0 DO 10 I = 1,N 10 DOT = DOT + A(I)*B(I) RETURN END SUBROUTINE DQRULE (IDEG, NQP, NCORD, PT, WT) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DUNAVANT QUADRATURE RULE FOR TRIANGLES, TO DEGREE = 17 C IN AREA COORDINATES (NCORD=3), OR UNIT COORDINATES (NCORD=2) C I.J.N.M.E. VOL. 21, PP.1129-1148, 1985 C INPUT IDEG=0,1,2,3,4,5, 6, 7, 8, 9,10,11,12,13,14,15,16,17,NQP=0 C OR NQP=1,1,3,4,6,7,12,13,16,19,25,27,33,37,42,48,52,61 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C WARNING: REQUIRES COMPILER FLAG FOR 40 CONTINUATION LINES,-Nl40 <===== CDP IMPLICIT REAL*8 (A-H,O,R-V,X-Z) C PT & WT ARE SINGLE PRECISION PARAMETER ( MAXDEG = 17, MAXDAT = 107 ) DIMENSION PT(NCORD,0:*), WT(0:*) DIMENSION AW(MAXDAT), A1(MAXDAT), A2(MAXDAT), A3(MAXDAT), 1 NQPDEG(MAXDEG), ISTART(MAXDEG), LINES(MAXDEG), 2 KOUNTS(MAXDAT) C IDEG = DEGREE OF POLYNOMIAL TO BE INTEGRATED, 0 TO MAXDEG C NQP = NUMBER OF QUADRATURE POINTS, USE 0 IF IDEG GOVERNS C NCORD = NUMBER OF PARAMETRIC DIMENSIONS: 3-AREA, 2-UNIT CORD C PT = RETURNED QUADRATURE COORDINATES, PT(NCORD,NQP) C WT = RETURNED QUADRATURE WEIGHTS, WT(NQP) C C NQPDEG = NUMBER QUADRATURE PTS FOR POLYNOMIAL DEGREE C ISTART = WHERE IDEG RULE DATA STARTS IN DATA TABLES C LINES = NUMBER OF LINES OF DATA FOR EACH RULE C KOUNTS = NUMBER OF TIMES THAT A RULE LINE IS USED C A1,A2,A3 = AREA COORDINATES OF TABLE POINT C AW = AREA WEIGHT OF TABLE POINT DATA NQPDEG /1,3,4,6,7,12,13,16,19,25,27,33,37,42,48,52,61/ DATA LINES /1,1,2,2,3, 3, 4, 5, 6, 6, 7, 8,10,10,11,13,15/ DATA ISTART /1,2,3,5,7,10,13,17,22,28,34,41,49,59,69,80,93/ DATA KOUNTS /1,3,1,3,3,3,1,3,3,3,3,6,1,3,3,6,1,3,3,3,6,1,3,3,3, 1 3,6,1,3,3,6,6,6,3,3,3,3,3,6,6,3,3,3,3,3,6,6,6,1,3, 2 3,3,3,3,3,6,6,6,3,3,3,3,3,3,6,6,6,6,3,3,3,3,3,3,6, 3 6,6,6,6,1,3,3,3,3,3,3,3,6,6,6,6,6,1,3,3,3,3,3,3,3, 4 3,6,6,6,6,6,6 / DATA AW / + 1.000000000000000, 0.333333333333333, -0.562500000000000, + 0.520833333333333, 0.223381589678011, 0.109951743655322, + 0.225000000000000, 0.132394152788506, 0.125939180544827, + 0.116786275726379, 0.050844906370207, 0.082851075618374, + -0.149570044467682, 0.175615257433208, 0.053347235608838, + 0.077113760890257, 0.144315607677787, 0.095091634267285, + 0.103217370534718, 0.032458497623198, 0.027230314174435, + 0.097135796282799, 0.031334700227139, 0.077827541004774, + 0.079647738927210, 0.025577675658698, 0.043283539377289, + 0.090817990382754, 0.036725957756467, 0.045321059435528, + 0.072757916845420, 0.028327242531057, 0.009421666963733, + 0.000927006328961, 0.077149534914813, 0.059322977380774, + 0.036184540503418, 0.013659731002678, 0.052337111962204, + 0.020707659639141, 0.025731066440455, 0.043692544538038, + 0.062858224217885, 0.034796112930709, 0.006166261051559, + 0.040371557766381, 0.022356773202303, 0.017316231108659, + 0.052520923400802, 0.011280145209330, 0.031423518362454, + 0.047072502504194, 0.047363586536355, 0.031167529045794, + 0.007975771465074, 0.036848402728732, 0.017401463303822, + 0.015521786839045, 0.021883581369429, 0.032788353544125, + 0.051774104507292, 0.042162588736993, 0.014433699669777, + 0.004923403602400, 0.024665753212564, 0.038571510787061, + 0.014436308113534, 0.005010228838501, 0.001916875642849, + 0.044249027271145, 0.051186548718852, 0.023687735870688, + 0.013289775690021, 0.004748916608192, 0.038550072599593, + 0.027215814320624, 0.002182077366797, 0.021505319847731, + 0.007673942631049, 0.046875697427642, 0.006405878578585, + 0.041710296739387, 0.026891484250064, 0.042132522761650, + 0.030000266842773, 0.014200098925024, 0.003582462351273, + 0.032773147460627, 0.015298306248441, 0.002386244192839, + 0.019084792755899, 0.006850054546542, 0.033437199290803, + 0.005093415440507, 0.014670864527638, 0.024350878353672, + 0.031107550868969, 0.031257111218620, 0.024815654339665, + 0.014056073070557, 0.003194676173779, 0.008119655318993, + 0.026805742283163, 0.018459993210822, 0.008476868534328, + 0.018292796770025, 0.006665632004165 / DATA A1 / + 0.333333333333333, 0.666666666666667, 0.333333333333333, + 0.600000000000000, 0.108103018168070, 0.816847572980459, + 0.333333333333333, 0.059715871789770, 0.797426985353087, + 0.501426509658179, 0.873821971016996, 0.053145049844817, + 0.333333333333333, 0.479308067841920, 0.869739794195568, + 0.048690315425316, 0.333333333333333, 0.081414823414554, + 0.658861384496480, 0.898905543365938, 0.008394777409958, + 0.333333333333333, 0.020634961602525, 0.125820817014127, + 0.623592928761935, 0.910540973211095, 0.036838412054736, + 0.333333333333333, 0.028844733232685, 0.781036849029926, + 0.141707219414880, 0.025003534762686, 0.009540815400299, + -0.069222096541517, 0.202061394068290, 0.593380199137435, + 0.761298175434837, 0.935270103777448, 0.050178138310495, + 0.021022016536166, 0.023565220452390, 0.120551215411079, + 0.457579229975768, 0.744847708916828, 0.957365299093579, + 0.115343494534698, 0.022838332222257, 0.025734050548330, + 0.333333333333333, 0.009903630120591, 0.062566729780852, + 0.170957326397447, 0.541200855914337, 0.771151009607340, + 0.950377217273082, 0.094853828379579, 0.018100773278807, + 0.022233076674090, 0.022072179275643, 0.164710561319092, + 0.453044943382323, 0.645588935174913, 0.876400233818255, + 0.961218077502598, 0.057124757403648, 0.092916249356972, + 0.014646950055654, 0.001268330932872, -0.013945833716486, + 0.137187291433955, 0.444612710305711, 0.747070217917492, + 0.858383228050628, 0.962069659517853, 0.133734161966621, + 0.036366677396917, -0.010174883126571, 0.036843869875878, + 0.012459809331199, 0.333333333333333, 0.005238916103123, + 0.173061122901295, 0.059082801866017, 0.518892500060958, + 0.704068411554854, 0.849069624685052, 0.966807194753950, + 0.103575692245252, 0.020083411655416, -0.004341002614139, + 0.041941786468010, 0.014317320230681, 0.333333333333333, + 0.005658918886452, 0.035647354750751, 0.099520061958437, + 0.199467521245206, 0.495717464058095, 0.675905990683077, + 0.848248235478508, 0.968690546064356, 0.010186928826919, + 0.135440871671036, 0.054423924290583, 0.012868560833637, + 0.067165782413524, 0.014663182224828 / DATA A2 / + 0.333333333333333, 0.166666666666667, 0.333333333333333, + 0.200000000000000, 0.445948490915965, 0.091576213509771, + 0.333333333333333, 0.470142064105115, 0.101286507323456, + 0.249286745170910, 0.063089014491502, 0.310352451033784, + 0.333333333333333, 0.260345966079040, 0.065130102902216, + 0.312865496004874, 0.333333333333333, 0.459292588292723, + 0.170569307751760, 0.050547228317031, 0.263112829634638, + 0.333333333333333, 0.489682519198738, 0.437089591492937, + 0.188203535619033, 0.044729513394453, 0.221962989160766, + 0.333333333333333, 0.485577633383657, 0.109481575485037, + 0.307939838764121, 0.246672560639903, 0.066803251012200, + 0.534611048270758, 0.398969302965855, 0.203309900431282, + 0.119350912282581, 0.032364948111276, 0.356620648261293, + 0.171488980304042, 0.488217389773805, 0.439724392294460, + 0.271210385012116, 0.127576145541586, 0.021317350453210, + 0.275713269685514, 0.281325580989940, 0.116251915907597, + 0.333333333333333, 0.495048184939705, 0.468716635109574, + 0.414521336801277, 0.229399572042831, 0.114424495196330, + 0.024811391363459, 0.268794997058761, 0.291730066734288, + 0.126357385491669, 0.488963910362179, 0.417644719340454, + 0.273477528308839, 0.177205532412543, 0.061799883090873, + 0.019390961248701, 0.172266687821356, 0.336861459796345, + 0.298372882136258, 0.118974497696957, 0.506972916858243, + 0.431406354283023, 0.277693644847144, 0.126464891041254, + 0.070808385974686, 0.018965170241073, 0.261311371140087, + 0.388046767090269, 0.285712220049916, 0.215599664072284, + 0.103575616576386, 0.333333333333333, 0.497380541948438, + 0.413469438549352, 0.470458599066991, 0.240553749969521, + 0.147965794222573, 0.075465187657474, 0.016596402623025, + 0.296555596579887, 0.337723063403079, 0.204748281642812, + 0.189358492130623, 0.085283615682657, 0.333333333333333, + 0.497170540556774, 0.482176322624625, 0.450239969020782, + 0.400266239377397, 0.252141267970953, 0.162047004658461, + 0.075875882260746, 0.015654726967822, 0.334319867363658, + 0.292221537796944, 0.319574885423190, 0.190704224192292, + 0.180483211648746, 0.080711313679564 / DATA A3 / + 0.333333333333333, 0.166666666666667, 0.333333333333333, + 0.200000000000000, 0.445948490915965, 0.091576213509771, + 0.333333333333333, 0.470142064105115, 0.101286507323456, + 0.249286745170910, 0.063089014491502, 0.636502499121399, + 0.333333333333333, 0.260345966079040, 0.065130102902216, + 0.638444188569810, 0.333333333333333, 0.459292588292723, + 0.170569307751760, 0.050547228317031, 0.728492392955404, + 0.333333333333333, 0.489682519198738, 0.437089591492937, + 0.188203535619033, 0.044729513394453, 0.741198598784498, + 0.333333333333333, 0.485577633383657, 0.109481575485037, + 0.550352941820999, 0.728323904597411, 0.923655933587500, + 0.534611048270758, 0.398969302965855, 0.203309900431282, + 0.119350912282581, 0.032364948111276, 0.593201213428213, + 0.807489003159792, 0.488217389773805, 0.439724392294460, + 0.271210385012116, 0.127576145541586, 0.021317350453210, + 0.608943235779788, 0.695836086787803, 0.858014033544073, + 0.333333333333333, 0.495048184939705, 0.468716635109574, + 0.414521336801277, 0.229399572042831, 0.114424495196330, + 0.024811391363459, 0.636351174561660, 0.690169159986905, + 0.851409537834241, 0.488963910362179, 0.417644719340454, + 0.273477528308839, 0.177205532412543, 0.061799883090873, + 0.019390961248701, 0.770608554774996, 0.570222290846683, + 0.686980167808088, 0.879757171370171, 0.506972916858243, + 0.431406354283023, 0.277693644847144, 0.126464891041254, + 0.070808385974686, 0.018965170241073, 0.604954466893291, + 0.575586555512814, 0.724462663076655, 0.747556466051838, + 0.883964574092416, 0.333333333333333, 0.497380541948438, + 0.413469438549352, 0.470458599066991, 0.240553749969521, + 0.147965794222573, 0.075465187657474, 0.016596402623025, + 0.599868711174861, 0.642193524941505, 0.799592720971327, + 0.768699721401368, 0.900399064086661, 0.333333333333333, + 0.497170540556774, 0.482176322624625, 0.450239969020782, + 0.400266239377397, 0.252141267970953, 0.162047004658461, + 0.075875882260746, 0.015654726967822, 0.655493203809423, + 0.572337590532020, 0.626001190286228, 0.796427214974071, + 0.752351005937729, 0.904625504095608 / C CHECK SPACE TYPE OPTIONS IF ( NCORD .LT. 2 .OR. NCORD .GT. 3 ) THEN STOP 'INVALID COORDINATE TYPE, DQRULE' ENDIF C CHECK FOR IDEG OR NQP CONTROL: LDEG = IDEG IF ( NQP .EQ. 0 ) THEN C USE DEGREE CONTROL IF ( IDEG .EQ. 0 ) NQP = 1 IF ( IDEG .LT. 0 .OR. IDEG .GT. MAXDEG ) THEN STOP 'INVALID IDEG ARGUMENT, DQRULE' ELSE NQP = NQPDEG(IDEG) ENDIF ELSE C USE NQP CONTROL LDEG = 0 DO 10 I = 1, MAXDEG IF ( NQP .EQ. NQPDEG(I) ) LDEG = I 10 CONTINUE IF ( LDEG .EQ. 0 ) STOP 'INVALID NQP ARGUMENT, DQRULE' ENDIF C FOUND VALID RULE, NOW EXPAND TABLE TO FULL RULE IPT = ISTART(LDEG) - 1 IRULE = 0 SUM = 0.D0 DO 20 I = 1, LINES(LDEG) J = IPT + I KOUNT = KOUNTS(J) IRULE = IRULE + 1 SUM = SUM + AW(J)*KOUNT WT(IRULE) = AW(J) PT(1,IRULE) = A1(J) PT(2,IRULE) = A2(J) IF ( NCORD .EQ. 3 ) PT(3,IRULE) = A3(J) IF ( KOUNT .GE. 3 ) THEN IRULE = IRULE + 1 WT(IRULE) = AW(J) PT(1,IRULE) = A3(J) PT(2,IRULE) = A1(J) IF ( NCORD .EQ. 3 ) PT(3,IRULE) = A2(J) IRULE = IRULE + 1 WT(IRULE) = AW(J) PT(1,IRULE) = A2(J) PT(2,IRULE) = A3(J) IF ( NCORD .EQ. 3 ) PT(3,IRULE) = A1(J) ENDIF IF ( KOUNT .EQ. 6 ) THEN IRULE = IRULE + 1 WT(IRULE) = AW(J) PT(1,IRULE) = A1(J) PT(2,IRULE) = A3(J) IF ( NCORD .EQ. 3 ) PT(3,IRULE) = A2(J) IRULE = IRULE + 1 WT(IRULE) = AW(J) PT(1,IRULE) = A3(J) PT(2,IRULE) = A2(J) IF ( NCORD .EQ. 3 ) PT(3,IRULE) = A1(J) IRULE = IRULE + 1 WT(IRULE) = AW(J) PT(1,IRULE) = A2(J) PT(2,IRULE) = A1(J) IF ( NCORD .EQ. 3 ) PT(3,IRULE) = A3(J) ENDIF 20 CONTINUE C CHECK VALIDITY OF RESULTS C IF ( SUM .NE. 1.D0 ) WRITE (6,*) C 1 'WARNING, UNITY NOT', SUM,', DQRULE' IF ( NCORD .EQ. 2 ) THEN DO 30 I = 1, NQP WT(I) = WT(I)*0.5D0 30 CONTINUE ENDIF RETURN END SUBROUTINE DSTART (IPRINT, M, NG, NSPACE, NDFREE, 1 INDEX, X, COORD, DD) C * * * * * * * * * * * * * * * * * * * * * * * * * * C INITIALIZE SYSTEM DOF FOR ITERATIVE SOLUTION C * * * * * * * * * * * * * * * * * * * * * * * * * * CDP IMPLICIT REAL*8(A-H,O-Z) PARAMETER ( NPRT = 6 ) DIMENSION DD(NDFREE), X(M,NSPACE), 1 COORD(1,NSPACE), INDEX(NG) C M = NUMBER OF NODES IN SYSTEM C NG = NUMBER OF PARAMETERS (DOF) PER NODE C NSPACE = DIMENSION OF SPACE C NDFREE = TOTAL NUMBER OF SYSTEM DOF C INDEX = SYSTEM DOF NUMBERS FOR DOF AT NODE C X = COORDINATES OF SYSTEM NODES C COORD = SPATIAL COORDINATE ARRAY OF A NODE C DD = SYSTEM ARRAY OF DEGREES OF FREEDOM C IPRINT > 0, PRINT THE STARTING VALUES IF ( IPRINT .GT. 0 ) WRITE (NPRT,5000) 5000 FORMAT ( /, 1 '** STARTING VALUES FOR ITERATIVE SOLUTION **',/, 2 'NODE PARAMETER VALUE') DO 20 I = 1, M C FIND PT COORDS AND DOF NOS CALL INDXPT (I,NG,INDEX) CALL PTCORD (I,M,NSPACE,X,COORD) DO 10 J = 1, NG INDX = INDEX(J) DD(INDX) = START(J,NSPACE,COORD) C START IS A FUNCTION TO DEFINE INITIAL VALUES C OF THE SYSTEM DEGREES OF FREEDOM IF ( IPRINT .GT. 0 ) WRITE (NPRT,5010) I,J,DD(INDX) 5010 FORMAT ( I5, I10, 2X, 1PE13.5 ) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE DYNDIM (NUMR, MAXR, NUMI, MAXI, R, RN, I, IN, 1 J, K ) C WARNING: COMPILE WITH EXTRA CONTINUE CARD OPTION -NL40 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DYNAMIC DIMENSION CONTROL FOR MODEL C C NOTES: C IF USING FORTRAN 90, THIS CODE IS REALLY NOT NEEDED C SINCE THE "AUTOMATIC ARRAY" FEATURE OF F90 WOULD C LET THE ARRAY SIZES BE DEFINED IN MODEL WHERE THEY C WOULD FIRST APPEAR. IN OTHER WORDS THEY WOULD BE C AUTOMATICALLY "ALLOCATED". THUS, YOU WOULD ONLY C NEED TO "ALLOCATE" THOSE ARRAYS THAT YOU LATER C WANT TO "DEALLOCATE" C C THIS FORTRAN 77 APPROACH TO DYNAMIC DIMENSIONS CALLED C DYNDIM IS ESSENTIALLY WRITTEN BY PROGRAM DIMMAK.F C USING ARRAYS "REALS" AND "INTEGERS", THEN EDITED C TO REMOVE BLANKS, ETC. IT ALSO WRITES MODEL CALL C AND DIMENSION STATEMENTS C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * PARAMETER ( MAXTYP=3) CHARACTER*8 RN, IN cb DIMENSION TITLE(15), R(MAXR), RN(NUMR), I(MAXI), IN(NUMI), CHARACTER*4 TITLE(15) DIMENSION R(MAXR), RN(NUMR), I(MAXI), IN(NUMI), 1 J(NUMR), K(NUMI) C C THE COMMONS BELOW ARE NOT USED, BUT CAN BE ACTIVATED C BY USERS IF THEY SEE A NEED FOR THEM. NOTE THAT THE C DUPLICATE ARRAY SIZE DATA ARE INCLUDED THERE ALSO. C C COMMON / REAL / R, RN, MMAXR, NNUMR, LLASTR, J C COMMON / INTEGER / MMAXI, NNUMI, LLASTI, IN, K C C I = VECTOR HOLDING ALL INTEGER ARRAYS (APPENDABLE) C IN = NAMES OF EACH INTEGER ARRAY (APPENDABLE) C J = POINTER TO BEGINNING OF EACH REAL ARRAY (APPENDABLE) C K = POINTER TO BEGINNING OF EACH INTEGER ARRAY (APPENDABLE) C LASTI = NUMBER OF THE LAST ASSIGNED INTEGER ARRAY C LASTR = NUMBER OF THE LAST ASSIGNED REAL ARRAY C LLASTI = NUMBER OF THE LAST ASSIGNED INTEGER ARRAY (COMMON) C LLASTR = NUMBER OF THE LAST ASSIGNED REAL ARRAY (COMMON) C MAXI = STORAGE ALLOWED FOR ALL INTEGER ARRAYS C MAXR = STORAGE ALLOWED FOR ALL REAL ARRAYS C MMAXI = STORAGE ALLOWED FOR ALL INTEGER ARRAYS (COMMON) C MMAXR = STORAGE ALLOWED FOR ALL REAL ARRAYS (COMMON) C NNUMI = ALLOWED NUMBER OF INTEGER ARRAYS (VIA COMMON) C NNUMR = ALLOWED NUMBER OF REAL ARRAYS (VIA COMMON) C NUMI = ALLOWED NUMBER OF INTEGER ARRAYS C NUMR = ALLOWED NUMBER OF REAL ARRAYS C R = VECTOR HOLDING ALL REAL ARRAYS (APPENDABLE) C RN = NAMES OF EACH REAL ARRAY (APPENDABLE) C C ** READ APPLICATION CONTROL DATA ** C 1 2 3 4 5 6 712 C23456789012345678901234567890123456789012345678901234567890-----------X CALL CONTROL (TITLE, M, NE, NG, N, NSPACE, NSEG, LBN, NITER, 1 NCURVE, INRHS, ISAY, NNPFIX, NNPFLO, NLPFIX, NLPFLO, 2 MISCFX, MISCFL, NHOMO, LHOMO, NPTWRT, LEMWRT, NTAPE1, 3 NTAPE2, NTAPE3, NTAPE4, NTAPE5, NULCOL, NDFREE, NELFRE, 4 NFLUX, IPTEST, LPTEST, NRB, NQP, LSHAPE, NLTYPE, 5 MODE, IBUG, NBSFIX, NBSFLO, NGF) C C SET ADVANCED USER DEFAULTS NC = N NCOEFF = NDFREE*2 NF = NGF NGEOM = N NPARM = NSPACE NPLT = 0 NTMP = 0 NOMAT = NLTYPE C jea 1/2/96 NSYS = 0 C --- Initialize real pointers --- J(1) = 1 K(1) = 1 c write(6,*) k C --- Calculate real array pointers --- C REAL 1, sys pt coords RN( 1) = 'X ' J( 2) = J( 1) + (M )*(NSPACE )*(1 )+0 C --- Calculate integer array pointers --- C INTEGER 1, all pts bc code IN( 1) = 'IBC ' K( 2) = K( 1) + (M )*(1 )*(1 )+0 C INTEGER 2, codes at a poin IN( 2) = 'KODES ' K( 3) = K( 2) + (NG )*(1 )*(1 )+0 C INTEGER 3, system topology IN( 3) = 'NODES ' K( 4) = K( 3) + (NE )*(N )*(1 )+0 C INTEGER 4, no constrain ty IN( 4) = 'NRES ' K( 5) = K( 4) + (MAXTYP )*(1 )*(1 )+0 C INTEGER 5, el type flag IN( 5) = 'LTYPE ' K( 6) = K( 5) + (NE )*(1 )*(1 )+0 C C INPUT NODES, BC FLAGS, ELEMENTS C C CALL INPUT (M, N, NE, NG, NSPACE, X, IBC, NODES, C LTYPE, NLTYPE) CALL INPUT (M, N, NE, NG, NSPACE, R(J(1)), I(K(1)), I(K(3)), 1 I(K(5)), NLTYPE) C C COUNT BC AND CONSTRAINTS, CONVERT NRES TO NREQ C C CALL CCOUNT (M, NG, NRES, IBC, KODES, MAXACT, NUMCE, C 1 MAXTYP, NREQ ) c fix NREQ below CALL CCOUNT (M, NG, I(K(4)), I(K(1)), I(K(2)), MAXACT, NUMCE, 1 MAXTYP, I(K(4)) ) IF ( IBUG .GT. 0 ) THEN C LIST THE ONLY ARRAYS KNOWN AT THIS POINT LOC1 = 1 LOC2 = 5 CALL LISTI (LOC1, LOC2, NEXTI, IN, K, I) LOC1 = 1 LOC2 = 1 CALL LISTR (LOC1, LOC2, NEXTR, RN, J, R) ENDIF C C --- CALCULATE REAL ARRAY POINTERS --- C C REAL 2, jacobian RN( 2) = 'AJ ' J( 3) = J( 2) + (NSPACE )*(NSPACE )*(1 )+0 C REAL 3, inverse jacobia RN( 3) = 'AJINV ' J( 4) = J( 3) + (NSPACE )*(NSPACE )*(1 )+0 C REAL 4, real pt average RN( 4) = 'AVE ' J( 5) = J( 4) + (M +1 )*(NRB+2 )*(1 )+1 C REAL 5, b matrix RN( 5) = 'B ' J( 6) = J( 5) + (NRB )*(NELFRE )*(1 )+0 C REAL 6, body force RN( 6) = 'BODY ' J( 7) = J( 6) + (NSPACE )*(1 )*(1 )+0 C REAL 7, el load vector RN( 7) = 'C ' J( 8) = J( 7) + (NELFRE )*(1 )*(1 )+0 C REAL 8, sys load vector RN( 8) = 'CC ' J( 9) = J( 8) + (NDFREE )*(1 )*(1 )+0 C REAL 9, constrain coeff RN( 9) = 'CEQ ' J( 10) = J( 9) + (MAXACT )*(NUMCE )*(1 )+0 C REAL 10, el or pt coord RN( 10) = 'COORD ' J( 11) = J( 10) + (N )*(NSPACE )*(1 )+0 C REAL 11, el or pt dof RN( 11) = 'D ' J( 12) = J( 11) + (NELFRE )*(1 )*(1 )+0 C REAL 12, old sys dof RN( 12) = 'DDOLD ' J( 13) = J( 12) + (NDFREE )*(1 )*(1 )+0 IF ( NITER .LT. 2 ) J(13) = J(12) + 1 C REAL 13, global deriv h RN( 13) = 'DGH ' cb J( 14) = J( 13) + (NSPACE )*(N )*(NQP +2 )+1 J( 14) = J( 13) + (NSPACE )*(nelfre )*(NQP +2 )+1 C REAL 14, local deriv g RN( 14) = 'DLG ' J( 15) = J( 14) + (NPARM )*(NGEOM )*(NQP +2 )+1 C REAL 15, local deriv h RN( 15) = 'DLH ' cb J( 16) = J( 15) + (NSPACE )*(N )*(NQP +2 )+1 J( 16) = J( 15) + (NSPACE )*(nelfre )*(NQP +2 )+1 C REAL 16, constitutive ma RN( 16) = 'E ' J( 17) = J( 16) + (NRB )*(NRB )*(1 )+0 C REAL 17, matrix product RN( 17) = 'EB ' J( 18) = J( 17) + (NRB )*(NELFRE )*(1 )+0 C REAL 18, real el propert RN( 18) = 'ELPROP ' J( 19) = J( 18) + (NLPFLO+1 )*(1 )*(1 )+1 C REAL 19, real prop all e RN( 19) = 'FLTEL ' J( 20) = J( 19) + (NE )*(NLPFLO+1 )*(1 )+1 C REAL 20, real prop segme RN( 20) = 'FLTBS ' J( 21) = J( 20) + (NSEG +1 )*(NBSFLO+1 )*(1 )+1 C REAL 21, real misc prop RN( 21) = 'FLTMIS ' J( 22) = J( 21) + (MISCFL+1 )*(1 )*(1 )+1 C REAL 22, real prop all p RN( 22) = 'FLTNP ' J( 23) = J( 22) + (M )*(NNPFLO+1 )*(1 )+1 C REAL 23, flux comps on e RN( 23) = 'FLUX ' J( 24) = J( 23) + (NF +1 )*(1 )*(1 )+1 C REAL 24, all flux on nod RN( 24) = 'FLUXBS ' J( 25) = J( 24) + (NSEG +1 )*(NFLUX +1 )*(1 )+1 C REAL 25, geom interpolat RN( 25) = 'G ' J( 26) = J( 25) + (NGEOM )*(NQP +2 )*(1 )+1 C REAL 26, gauss pts 1-d RN( 26) = 'GPT ' J( 27) = J( 26) + (NQP +1 )*(1 )*(1 )+1 C REAL 27, gauss wts 1-d RN( 27) = 'GWT ' J( 28) = J( 27) + (NQP +1 )*(1 )*(1 )+1 C REAL 28, solution interp RN( 28) = 'H ' cb J( 29) = J( 28) + (N )*(NQP +2 )*(1 )+1 J( 29) = J( 28) + (nelfre )*(NQP +2 )*(1 )+1 C REAL 29, integral of h RN( 29) = 'HINTG ' cb J( 30) = J( 29) + (N )*(NQP+3 )*(1 )+0 J( 30) = J( 29) + (nelfre )*(NQP+3 )*(1 )+0 C REAL 30, plotter data RN( 30) = 'PLTSET ' J( 31) = J( 30) + (NPLT +1 )*(1 )*(1 )+1 C REAL 31, real prop el no RN( 31) = 'PRTLPT ' J( 32) = J( 31) + (N )*(NNPFLO+1 )*(1 )+1 C REAL 32, real mat numb p RN( 32) = 'PRTMAT ' J( 33) = J( 32) + (NLPFLO+1 )*(NOMAT +1 )*(1 )+1 C REAL 33, quadrature coor RN( 33) = 'PT ' J( 34) = J( 33) + (NPARM )*(NQP +1 )*(1 )+0 C REAL 34, dof max min val RN( 34) = 'RANGE ' J( 35) = J( 34) + (NG )*(2 )*(1 )+0 C REAL 35, el or edge sq m RN( 35) = 'S ' J( 36) = J( 35) + (NELFRE )*(NELFRE )*(1 )+0 C REAL 36, stress at el po RN( 36) = 'SATPT ' J( 37) = J( 36) + (NRB+2 )*(N )*(1 )+0 C REAL 37, strain or grad RN( 37) = 'STRAIN ' J( 38) = J( 37) + (NRB+2 )*(1 )*(1 )+0 C REAL 38, initial strain RN( 38) = 'STRAN0 ' J( 39) = J( 38) + (NRB )*(1 )*(1 )+0 C REAL 39, stress + rms or RN( 39) = 'STRESS ' J( 40) = J( 39) + (NRB+2 )*(1 )*(1 )+0 C REAL 40, sys control dat RN( 40) = 'SYSDAT ' J( 41) = J( 40) + (NSYS +1 )*(1 )*(1 )+1 C REAL 41, temporary work RN( 41) = 'TMP ' J( 42) = J( 41) + (NTMP +1 )*(1 )*(1 )+1 C REAL 42, values at corne RN( 42) = 'VALC ' J( 43) = J( 42) + (NRB )*(NC +1 )*(1 )+1 C REAL 43, values on edge RN( 43) = 'VALE ' J( 44) = J( 43) + (NRB )*(NC +1 )*(1 )+1 C REAL 44, quadrature weig RN( 44) = 'WT ' J( 45) = J( 44) + (NQP +1 )*(1 )*(1 )+0 C REAL 45, xy ends of a li RN( 45) = 'XPT ' J( 46) = J( 45)