C MAIN ROUTINE OF THE PROGRAM PASSFEM C IMPLICIT REAL8(A-H,O-Z) REAL4 S1,S2,S3,S4,S5,S6,S7,S8,SS1,SS2,SS3,SS4 INTEGER CHT COMMON/DIM/N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13,N14 COMMON/PAR/IND,NET,NSN,NMP,NEQ,NSKY,NEQ1,LCOUNT COMMON/TAPES/ISTRES,NDARAY,IPR COMMON/PRECI/ITWO COMMON/MULT/ELMN COMMON/ELEM/NDEL(100) COMMON/ELOAD/M1,M2,AA1,AA2,AA3,AA4,NSLC,NUM DIMENSION SS1(10),SS2(10),SS3(10),SS4(10) DIMENSION LLIB(7,3),NLN(10),TITLE(20),ELM(10) CHARACTER*15 INFILE,OUTFILE C C DIMENSION A(30000) MTOT=30000 CALL SECOND(S1,0.0) C ITWO=1 C C ISTRES=1 NDARAY=2 C
C
C READ & PRINT STRUCTURE DATA
100 READ(5,55)(TITLE(I),I=1,20)
READ(5,11)NSN,NET,NMP,NLC,NEDMAX,MODEX
IF(NSN.EQ.0)STOP
READ(5,12)IPR
WRITE(6,66)(TITLE(I),I=1,20)
WRITE(6,22)NSN,NET,NMP,NLC
C READ,GENERATE & PRINT NODAL DATA & MATERIAL PROPERTIES
N1=1
N2=N1+6NSN
N3=N2+NSNITWO
N4=N3+NSNITWO
N5=N4+NSNITWO
IF(N5.GT.MTOT)CALL ERROR(N5-MTOT,1)
IF(N5.GT.MTOT)STOP
N6=N5+NMPITWO
N7=N6+NMPITWO
N8=N7+NMPITWO
CALL PASSIN (A(N2),A(N3),A(N4),A(N5),A(N6),A(N7),A(N1),
NSN,NMP,NEQ)
C INITIALIZE COLUMN HEIGHT TO '0'
N9=N8+NEQ
DO 10 I=N8,N9
10 A(I)=0.0
C READ AND GENERATE ELEMENT DATA
WRITE(6,33)
NEQ1=NEQ+1
REWIND NDARAY
IND=1
LCOUNT=0
CALL FELIB(A,LLIB,MTOT)
CALL SECOND(S2,S1)
C
C ADDRESS DIAGONAL ELEMENTS OF GSM FOR ASSEMBLY
N11=N10+NEQ1
CALL CADNUM (A(N8),A(N10),NEQ,NEQ1,NSKY,MBAND)
C WRITE DATA FOR SOLUTION PHASE
WRITE(6,88)NEQ,NSKY
C ALLOT SPACE FOR STIFFNESS MATRIX & LOAD VECTOR
N12=N11+NSKYITWO
N13=N12+NEQITWO
N14=N13+NEQ*ITWO
IF(N14.GT.MTOT)CALL ERROR(N14-MTOT,3)
IF(N14.GT.MTOT)STOP
C INITIALIZE GSM,LOAD VECTOR & DISPALCEMENT VECTOR
NSIZE=NSKY+NEQ+NEQ
CALL INISKP(A(N11),NSIZE)
IF (MODEX.NE.0)GO TO 222
C
C CAL LOAD VECTOR & ESM & ASSEMBLE STR STIFFNESS MATRIX
REWIND NDARAY
REWIND ISTRES
IND=2
LCOUNT=0
CALL FELIB(A,LLIB,MTOT)
CALL SECOND(S3,S1)
C
C TRIANGULARIZE STIFFNESS MATRIX
KTR=1
CALL PASOLV(A(N11),A(N13),A(N10),NEQ,NEQ1,NSKY,KTR)
222 CONTINUE
CALL SECOND(S4,S1)
C
C READ & PRINT ELEMENT LOAD MULTIPLIERS
C
READ (5,77)(ELM(J),J=1,NLC)
WRITE(6,78)
WRITE(6,79)(L,L=1,NLC)
WRITE(6,81)(ELM(J),J=1,NLC)
READ(5,44)(NLN(J),J=1,NLC)
NN1=N13-1
CALL SECOND(S5,S1)
S6=0.0
S7=0.0
S8=0.0
DO 80 L=1,NLC
CALL SECOND(SS1(L),0.0)
ELMN=ELM(L)
DO 144 I=N12,NN1
144 A(I+NEQ)=A(I)*ELMN
WRITE(6,99) L
WRITE(6,122)
NC=NLN(L)
IF(NC.EQ.0) GOTO 30
C
C ADD CONC LOAD TO CAL LOAD VECTOR
C
CALL PASLOD(A(N13),A(N1),NC,NSN,NEQ)
C
30 CONTINUE
CALL SECOND(SS3(L),SS1(L))
IF(MODEX.NE.0)GO TO 80
C CAL OF DISPLACEMENTS
KTR=2
CALL PASOLV(A(N11),A(N13),A(N10),NEQ,NEQ1,NSKY,KTR)
C
C PRINT NODAL DISPLACEMENTS
WRITE(6,111)
CALL DISP(A(N13),A(N1),NSN,NEQ)
CALL SECOND(SS3(L),SS1(L))
C
C CAL STRESSES AT SELECTED POINTS OF EACH ELEMENT
REWIND NDARAY
REWIND ISTRES
IND=3
LCOUNT=0
CALL FELIB(A,LLIB,MTOT)
CALL SECOND(SS4(L),SS1(L))
S6=S6+SS2(L)
S7=S7+SS3(L)-SS2(L)
S8=S8+SS4(L)-SS3(L)
80 CONTINUE
IF (MODEX.NE.0) GO TO 200
S1=S2+(S5-S4)+S6
S2=S3-S2
S3=(S4-S3)+S7
S4=S8
S5=S1+S2+S3+S4
WRITE(6,130)
WRITE(6,135)S1,S2,S3,S4,S5
200 CONTINUE
C
C READ NEXT ANALYSIS CASE
C
GO TO 100
1 FORMAT(A15)
11 FORMAT(6I5)
12 FORMAT(I5)
22 FORMAT(53X,'NUMBER OF STRUCTURE NODES....',
*I4/53X,'NUMBER OF ELEMENT TYPES....',I4/53X,
*'NUMBER OF MATERIAL GROUPS....',I4/53X,
'NUMBER OF LOADING CONDITIONS.',I3//)
33 FORMAT(//61X,'ELEMENT DATA'/)
44 FORMAT(10I5)
55 FORMAT(20A4)
66 FORMAT(130('')//27X,20A4//59X,'STRUCTURE DATA',/)
77 FORMAT(10F5.2)
78 FORMAT(//5X,'ELEMENT LOAD MULTIPLIERS'//)
79 FORMAT(5X,'LOAD CASE NO',I4,9I8)
81 FORMAT(15X,10(2X,F6.2))
88 FORMAT(/51X,'DATA OF SOLUTION PHASE'
*/48X,'NUMBER OF EQUATIONS (NEQ)=',
*I7//48X,'NUMBER OF ELEMENT IN SK (NSKY)=',I7//)
99 FORMAT(//5X,'LOAD CASE NO',I4/)
111 FORMAT(///56X,'NODAL DISPLACEMENTS'//5X,'NODE',8X,
*'X-DISPALCEMENTS',6X,'Y-DISPLACEMENTS',4X,'Z-DISPALCEMENTS'
*,8X,'ROTATION-X',9X,'ROTATION Y',9X,'ROTATION-Z')
122 FORMAT(/58X,'NODAL LOAD DATA'/)
130 FORMAT(//5X,'OVER ALL TIMELOG'//)
135 FORMAT(10X,'MODEL INPUT',F8.2/
*10X,'ESM FORMULATION'/10X,'AND ASSEMBLAGE',F8.2/
*10X,'EQUATION SOLUTION AND'/10X,'DISPLACEMENT OUTPUT',F8.2/
*10X,'STRESS OUTPUT',F8.2//10X,'TOTAL SOLUTION TIME',F8.2//)
END
C
C SUBROUTINE TO DETERMINE THE TIME TAKEN FOR EXECUTION CALLED BY MAIN
C
SUBROUTINE SECOND(T2,T1)
T2=SECNDS(T1)
RETURN
END
C
C ROUTINE TO PRINT MSGS WITH SPECIFIED MAIN STORAGE CALLED BY MAIN
C
SUBROUTINE ERROR(N,I)
GO TO(1,2,3),I
1 WRITE(6,11)
GO TO 10
2 WRITE(6,22)
GO TO 10
3 WRITE(6,33)
10 WRITE(6,44)N
RETURN
11 FORMAT(7X,'NOT ENOUGH STORAGE FOR READ-IN OF NDF ARRAY & NODAL
*POINT CO ORDINATES')
22 FORMAT(//////5X,'NOT ENOUGH STORAGE FOR ELEMENTAL NODE
*COORDINATES')
33 FORMAT(//////'NOT ENOUGH STORAGE FOR ASSEMBLAGE OF GLOBAL
STRUCTURE STIFFNESS & DISPLACEMENT IN STRESS SOLUTION PHASE')
44 FORMAT(////'ERRORSTORAGE EXCEEDED BY',I7)
END
C
C
C ROUTINE READS,GENERATES & PRINTS NOAL DATA & MATERIAL PROPERTIES CALLED BY MAIN
C
C
SUBROUTINE PASSIN(X,Y,Z,E,PR,WD,NDF,NSN,NMP,NEQ)
IMPLICIT REAL8(A-H,O-Z)
COMMON/PRECI/ITWO
COMMON/TAPES/ISTRES,NDARAY,IPR
DIMENSION X(NSN),Y(NSN),Z(NSN),E(NMP),PR(NMP),WD(NMP)
DIMENSION NDF(6,NSN),JF(6)
DIMENSION NCYLT(1000)
FLOAT(I)=DBLE(I)
MN=0
100 READ(5,33)NN,(JF(I),I=1,6),X(NN),Y(NN),Z(NN),NI,NCYL
NCYLT(NN)=NCYL
N=MN+NI
MN=MN+1
110 DO 120 I=1,6
NDF(I,NN)=JF(I)
120 CONTINUE
IF (NI.EQ.0) GO TO 130
IF (NN-MN)130,125,140
125 CONTINUE
IF(NSN-NN)170,170,100
130 MN=NN
GO TO 125
C
C AUTOMATIC GENERATION OF NODAL DATA
C
140 NX=(NN-N+NI)/NI
XD=(X(NN)-X(N-NI))/FLOAT(NX)
YD=(Y(NN)-Y(N-NI))/FLOAT(NX)
ZD=(Z(NN)-Z(N-NI))/FLOAT(NX)
MN=NN
150 X(N)=X(N-NI)+XD
Y(N)=Y(N-NI)+YD
Z(N)=Z(N-NI)+ZD
NCYLT(N)=NCYLT(NN)
DO 160 I=1,6
NDF(I,N)=JF(I)
160 CONTINUE
N=N+NI
IF(N.LT.NN) GO TO 150
IF(NSN-NN)170,170,100
170 CONTINUE
DO 180 I=1,NSN
IF(NCYLT(I).EQ.0) GO TO 180
THETA=Y(I)*3.14159/180.0
Y(I)=X(I)*DSIN(THETA)
Y(I)=X(I)*DCOS(THETA)
180 CONTINUE
IF(IPR.EQ.0)WRITE(6,44)(I,(NDF(J,I),J=1,6),X(I),Y(I),Z(I),I=1,NSN)
C CONVERT '0' & '1' OF 'NDF' ARRAY TO EQUATION NUMBERS & '0'S.
NEQ=0
DO 30 N=1,NSN
DO 30 I=1,6
IF(NDF(I,N))10,20,10
20 NEQ=NEQ+1
NDF(I,N)=NEQ
GO TO 30
10 NDF(I,N)=0
30 CONTINUE
IF(IPR.EQ.0)WRITE(6,77)(I,(NDF(J,I),J=1,6),I=1,NSN)
C
C READ & PRINT MATERIAL PROPERTIES
DO 1 J=1,NMP
1 READ(5,55)I,E(I),PR(I),WD(I)
WRITE(6,66)(I,E(I),PR(I),WD(I),I=1,NMP)
C
33 FORMAT(7I5,3F10.3,2I5)
44 FORMAT(//59X,'NODAL POIINT DATA'//5X,'NODE',
*3X,'NODAL D.O.F.',5X,'X-COORD.',5X,'Y-COORD.',5X,'Z-COORD.',
*6X,'NODE',3X,'NODAL DOF',5X,'X-COORD.',5X,'Y-COORD.',5X,
*'Z-COORD.',//(5X,I4,3X,6I2,3X,F10.4,3X,F10.4,3X,F10.4,3X,3X,
*I4,3X,6I2,3X,F10.4,3X,F10.4,3X,F10.4/))
55 FORMAT(I10,E10.3,2F10.4)
66 FORMAT(//57X,'MATERIAL PROPERTIES'//44X,'GROUP',7X,'YOUNGS',
*7X,'POISSON',7X,'WEIGHT',/45X,'NO.',7X,'MODULUS',8X,'RATIO',
7X,'DENSITY'//(45X,I3,5X,F10.2,1X,F10.2,6X,F10.5))
77 FORMAT(/24X,'NODE',11X,'EQUATION NUMBERS',25X,'NODE',12X,
'EQUATION NUMBERS'//(24X,I4,3X,6I5,20X,I4,3X,6I5/))
RETURN
END
C
C
SUBROUTINE FELIB(A,LLIB,MTOT)
IMPLICIT REAL8(A-H,O-Z)
COMMON/DIM/N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13,N14
COMMON/PAR/IND,NET,NSN,NMP,NEQ,NSKY,NEQ1,LCOUNT
COMMON/TAPES/ISTRES,NDARAY,IPR
COMMON/PRECI/ITWO
COMMON/MULT/ELMN
DIMENSION A(MTOT)
DIMENSION LLIB(7,3),LT(7)
DO 100 I=1,NET
IF(IND.NE.1)GO TO 5
READ(5,11)LTYPE,NSHAPE,(LLIB(LTYPE,J),J=1,NSHAPE)
LT(I)=LTYPE
5 LA=LT(I)
GO TO (50),LA
C 10 CALL THREDT(A,LLIB,NSHAPE,MTOT)
C GO TO 100
C 20 CALL THREDB(A,LLIB,NSHAPE,MTOT)
C GO TO 100
C 30 CALL PLANE(A,LLIB,NSHAPE,MTOT)
C GO TO 100
C 40 CALL THREDS(A,LLIB,NSHAPE,MTOT)
C GO TO 100
50 CALL PLATE(A,LLIB,NSHAPE,MTOT)
C GO TO 100
C 60 CALL SHELL(A,LLIB,NSHAPE,MTOT)
C GO TO 100
C 70 CALL BOUND(A,LLIB,NSHAPE,MTOT)
100 CONTINUE
11 FORMAT(5I5)
RETURN
END
SUBROUTINE INISKP(SK, NSIZE)
IMPLICIT REAL8(A-H,O-Z)
DIMENSION SK(NSIZE)
DO 10 I=1,NSIZE
10 SK(I)=0.0
RETURN
END
SUBROUTINE COLUMH(CHT,ND,NED,NEQ)
IMPLICIT REAL*8(A-H,O-Z)
INTEGER CHT(NEQ),ND(NED)
C CALCULATES THE COLUMN HEIGHT S OF EACH COLUMN IN THE GLOBAL STIFFNESS MATRIX LS=100000 DO 30 K=1,NED IF(ND(K)) 10,30,10 10 IF(ND(K)-LS)20,30,30 20 LS=ND(K) 30 CONTINUE DO 40 K=1,NED II=ND(K) IF(II.EQ.0)GO TO 40 ME=II-LS IF(ME.GT.CHT(II))CHT(II)=ME 40 CONTINUE RETURN END
SUBROUTINE CADNUM(CHT,NDS,NEQ,NEQ1,NSKY,MBAND)
IMPLICIT REAL*8(A-H,O-Z)
INTEGER CHT(NEQ),NDS(NEQ1)
C CALCULATES ADDRESSES OF DIAGONAL ELEMENTS IN BANDED MATRIX WHOSE COLOUMN HEIGHTS ARE KNOWN; C CALCULATES THE NO OF ELEMENTS IN THE GLOBAL STIFFNESS MATRIX BELOW THE SKYLINE C DO 10 I =1,NEQ1 10 NDS(I)=0 NDS(1)=1 NDS(2)=2 MBAND=0 IF(NEQ .EQ.1) GO TO 30 DO 20 I=2,NEQ IF(CHT(I).GT.MBAND) MBAND=CHT(I) 20 NDS(I+1)=NDS(I)+CHT(I)+1 30 MBAND=MBAND+1 NSKY=NDS(NEQ1)-1 RETURN END
SUBROUTINE PASSEM(SK,EK,NDS,ND,NED,NEQ1,NSKY,NUED)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION SK(NSKY),NDS(NEQ1), ND(NED),EK(NUED,NUED)
C ASSEMBLE ELEMENT STIFFNESS INTO COMPACTED GLOBAL STIFFNESS DO 70 I=1,NED II=ND(I) IF(II)70,70, 30 30 CONTINUE DO 60 J=1, NED JJ=ND(J) IF (JJ)60,60,40 40 CONTINUE MI=NDS(JJ) IJ=JJ-II IF(JJ)60,50,50 50 KK=MI+IJ SK(KK)=SK(KK)+EK(I,J) 60 CONTINUE 70 CONTINUE RETURN END
SUBROUTINE PASOLV(SK,P,NDS,NN,NEQ1,NSKY,KKK)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION SK(NSKY),P(NN),NDS(NEQ1)
IF(KKK-2) 40,150,150
40 DO 140 N=1,NN KN=NDS(N) KL=KN+1 KU=NDS(N+1)-1 KH=KU-KL IF(KH) 110,90,50 50 K=N-KH IC=0 KLT=KU DO 80 J=1,KH IC=IC+1 KLT=KLT-1 KI=NDS(K) ND=NDS(K+1)-KI-1 IF(ND) 80, 80, 60 60 KK=MIN0(IC,ND) C=0.0 DO 70 L=1, KK 70 C=C+SK(KI+L)SK(KLT+L) SK(KLT)=SK(KLT)-C 80 K=K+1 90 K=N B=0.0 DO 100 KK=KL,KU K=K-1 KI=NDS(K) C=SK(KK)/SK(KI) B=B+CSK(KK) 100 SK(KK)=C SK(KN)=SK(KN)-B 110 IF(SK(KN)) 120,120, 140 120 WRITE(6,222) N, SK(KN) STOP 140 CONTINUE RETURN C C REDUCE RIGHT HAND SIDE LOAD VECTOR 150 DO 180 N=1,NN KL=NDS(N)+1 KU=NDS(N+1)-1 IF(KU-KL) 180,160,160 160 K=N C=0.0 DO 170 KK=KL,KU K=K-1 170 C=C+SK(KK)*P(K) P(N)=P(N)-C 180 CONTINUE C BACK SUBSTITUTION DO 200 N=1,NN K=NDS(N) 200 P(N)=P(N)/SK(K) IF(NN.EQ.1) RETURN N=NN DO 230 L=2,NN KL=NDS(N)+1 KU=NDS(N+1)-1 IF(KU-KL) 230,210,210 210 K=N DO 220 KK=KL,KU K=K-1 220 P(K)=P(K)-SK(KK)*P(N) 230 N=N-1 RETURN 222 FORMAT(//20X,'STOP-STIFFNESS MATRIX NOT POSITIVE DEFINITE', *'NON POSITIVE PIVOT FOR EQUATION',I4,//10X,'PIVOT=',E20.12) END
SUBROUTINE PASLOD(P,NDF,NC,NSN,NEQ)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/TAPES/ISTRES,NDARAY,IPR
DIMENSION P(NEQ),NDF(6,NSN)
DIMENSION CNL(6)
IF(IPR.EQ.0) WRITE (6,30)
DO 20 J=1,NC
READ(5,11) NODE,(CNL(I),I=1,6)
IF (IPR.EQ.0) WRITE (6,40) NODE,(CNL(I),I=1,6)
DO 20 I=1,6
II=NDF(I,NODE)
IF(II) 20,20,10
10 P(II)=P(II)+CNL(I) 20 CONTINUE RETURN 11 FORMAT(I10,6F10.4) 30 FORMAT(10X,'NODE',10X,'X-AXIS',10X,'Y-AXIS',10X,'Z-AXIS', *10X,'X-AXIS',10X,'Y-AXIS',10X,'Z-AXIS'/25X,'FORCE',10X,'FORCE', *10X,'FORCE',10X,'MOMENT',10X,'MOMENT',10X,'MOMENT'/) 40 FORMAT(10X,I4,6F16.3) END
SUBROUTINE DISP(D,NDF,NSN,NEQ)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION D(NEQ),NDF(6,NSN),DISPV(6)
DO 30 J=1,NSN
DO 10 I=1,6
10 DISPV(I)=0.0 DO 20 I=1,6 KK=NDF(I,J) 20 IF(KK.NE.0) DISPV(I)=D(KK) 30 WRITE(6,22) J,(DISPV(I),I=1,6) 22 FORMAT(5X,I3,4X,E19.7,2E19.7,1X,E20.6,2E18.6/) RETURN END C INTEFRACE ROUTINE FROM DEC TOPC FOR GETTING TIME FUNCTION SECNDS(X) C RETURN CURRENT TIME MIDNIGHT- X REAL*4 SECNDS,X
INTEGER*2 IHOUR, IMINUT, ISECON, IHUND
REAL*4 RHOUR,RMINUT, RSECON, RHUND
REAL*4 X1
C CALL GETTIM(IHOUR, IMINUT , ISECON, IHUND) RHOUR = FLOAT( IHOUR ) RMINUT =FLOAT ( IMINUT ) RSECON = FLOAT ( ISECON) RHUND = FLOAT( IHUND ) X1 = RHOUR3600.0+ RMINUT60.0 + RSECON + RHUND/100.0 SECNDS = X1-X RETURN END C CONTROL ROUTINE FOR PLATE BENDING ELEMENTS
SUBROUTINE PLATE(A,LLIB,NSHAPE,MTOT)
IMPLICIT REAL*8(A-H,O-Z)
COMMON/DIM/N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N13,N14
COMMON/PAR/IND,NET,NSN,NMP,NEQ,NSKY,NEQ1,LCOUNT
COMMON/TAPES/ISTRES,NDARAY,IPR
COMMON/PRECI/ITWO
DIMENSION A(MTOT)
DIMENSION LLIB(7,3)
DO 100 J=1,NSHAPE
LB=LLIB(5,J)
GO TO (40),LB
40 NED=12 NSP=4 NS=5 IF(IND.NE.1) GO TO 50 N10=N9+NED N11=N10+1ITWO N12=N11+1ITWO N13=N12+1*ITWO NSKY=1 50 CALL PLATE4(A(N2),A(N3),A(N4),A(N5),A(N6),A(N7),A(N11),A(N12), A(N1),A(N8),A(N9),A(N10),A(N13),NS,NED,NSP) GO TO 100 100 CONTINUE RETURN END C*********************************************************************************** SUBROUTINE PLATE4(X,Y,Z,E,PR,WD,SK,PQ,NDF,CHT,ND,NDS,P,NS, *NED,NSP) C A BILINEAR ISOPARAMETRIC PLATE BENDING ELEMENT-A 4 NODED QUADRILATERAL C FINITE ELEMENT FORMULATION FOR THE SOLUTION OF PLATE BENDING PROBLEMS. C ISOPARAMETRIC CONCEPT AND SELECTIV INTEGRATION TECHNIQUES ARE ADOPTED
IMPLICIT REAL*8(A-H,O-Z)
COMMON/DIM/N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N13,N14
COMMON/PAR/IND,NET,NSN,NMP,NEQ,NSKY,NEQ1,LCOUNT
COMMON/TAPES/ISTRES,NDARAY,IPR
DIMENSION X(NSN),Y(NSN),Z(NSN),E(NMP),PR(NMP),WD(NMP)
DIMENSION SK(NSKY),P(NEQ),NDF(6,NSN),ND(NED),NDS(NEQ1)
DIMENSION ELC(100),STR(100,3)
INTEGER CHT(NEQ)
DIMENSION Q(12),XL(4),YL(4),SF(3,4),B(5,12),WG(4),R(4),S(4)
DIMENSION U(4),V(4),C(5,5),LNC(4),CBS(5,12),CBB(5,12)
DIMENSION XJACI(2,2),XJAC(2,2),SFG(3,4),CB(5,12),MNC(4)
DIMENSION EK(12,12),SIG(6),TMS(4,4),NDEL(21),EKB(12,12)
DIMENSION EKS(12,12),PQ(NEQ)
C SQRT(X) = DSQRT(X) GO TO (310,350),IND 310 READ(5,11) NEL LCOUNT=LCOUNT+1 NDEL(LCOUNT)=NEL IF(NSN.GT.100) WRITE(6,133) IF(NSN.GT.100) STOP DO 314 I=1,NSN 314 ELC(I)=0.0 WRITE(6,33) NEL WRITE(6,44) KK =0 DO 340 I=1,NEL READ(5,22) LNUM,MG,(MNC(J),J=1,4),TH,UDL,KI IF(KI.GT.0) GO TO 316 DO 315 J=1,4 315 LNC(J)=MNC(J) 316 KK =KK+1 IF(LNUM-KK) 2,2,3 3 CONTINUE DO 5 J=1,4 5 LNC(J)=LNC(J)+KI 2 JJ=0 DO 330 J=1,4 K = LNC(J) ELC(K)=ELC(K)+1.0 XL(J)=X(K) YL(J)=Y(K) DO 320 L=3,5 JJ=JJ+1 320 ND(JJ)=NDF(L,K) 330 CONTINUE WRITE(NDARAY)(ND(LL),LL=1,12),(XL(NN),YL(NN),LNC(NN),NN=1,4), TH,MG,UDL CALL COLUMH(CHT,ND,NED,NEQ) IF(IPR.EQ.0)WRITE(6,55)KK,MG,(LNC(J),J=1,4),TH,UDL KK=KK+1 IF(LNUM.GE.KK) GO TO 3 KK=KK-1 IF(LNUM.EQ.NEL) GO TO 341 340 CONTINUE 341 CONTINUE RETURN C****************************************************************************************** C CALCULATE ELEMENT STIFFNESS MATRIX EK 350 CONTINUE DATA R/-1.D0,1.D0,1.D0,-1.D0/ DATA S/-1.D0,-1.D0,1.D0,1.D0/ DATA U/-1.D0,1.D0,1.D0,-1.D0/ DATA V/-1.D0,-1.D0,1.D0,1.D0/ DATA WG/1.D0,1.D0,1.D0,1.D0/ DO 10 I=1,5 DO 10 J=1,5 C(I,J)=0.0 10 CONTINUE C ENTER LOOP OVER EACH ELEMENT LCOUNT=LCOUNT+1 NEL=NDEL(LCOUNT) DO 200 LNUM=1,NEL READ(NDARAY)(ND(LL),LL=1,12),(XL(NN),YL(NN),LNC(NN),NN=1,4), TH,MG,UDL IF(NMP.EQ.1.AND.LNUM.GT.1)GO TO 15 C COMPUTE CONSTITUTIVE MATRIX C C(1,1)=E(MG)TH**3/(12. 0(1.0-PR(MG)**2)) C(1,2)=C(1,1)PR(MG) C(2, 1)=C(1,2) C(2,2)=C(1,1) C(3,3)=C(1,1)(1.0-PR(MG))/2.0 C(4,4)=E(MG)TH/(2.4(1.0+PR(MG))) C(5,5)=C(4,4) 15 CONTINUE C SHEAR ENERGY TERMS EKS (11 ITEGRATION IS USED) C DO 175 I=1,12 DO 175 J=1,12 175 EKS(I,J)=0.0 GPSX=0.0 GPSY=0.0
C EVALUATE SHAPE FUNCTIONS DO 182 I=1,4 AA=(1.0+R(I)GPSX) BB=(1.0+S(I)GPSY) SF(1,I)=0.25R(I)BB SF(2,I)=0.25S(I)AA SF(3,I)=0.25AABB 182 CONTINUE C JACOBIAN MATRIX XJAC DO 184 I=1,2 DO 184 J=1,2 184 XJAC(I,J)=0.0 DO 186 I=1,4 DO 186 K=1,2 XJAC(K,1)=XJAC(K,1)+XL(I)*SF(K,I) XJAC(K,2)=XJAC(K,2)+YL(I)*SF(K,I) 186 CONTINUE C CALCULATE DETERMINANT JACOBIAN DJAC DJAC=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)XJAC(2,1) C COMPUTE JACOBIAN INVERSE XJACI XJACI(1,1)=XJAC(2,2)/DJAC XJACI(1,2)=-XJAC(1,2)/DJAC XJACI(2,1)=-XJAC(2,1)/DJAC XJACI(2,2)=XJAC(1,1)/DJAC DB = 4.0DJAC C FORM GLOBAL DERIVATIVES SFG DO 188 I=1,2 DO 188 J=1,4 188 SFG(I,J)=0.0 DO 190 I=1,2 DO 190 K=1,4 DO 190 J=1,2 SFG(I,K)=SFG(I,K)+XJACI(I,J)*SF(J, K) 190 CONTINUE DO 192 I=1,4 SFG(3,I)=SF(3,I) 192 CONTINUE C COMPUTE CBS-THE SHEAR CONTRIBUTION TO STRESSMATRIX CB
DO 197 I=1,NS
DO 197 J=1,12
197 B(I,J)=0.0 DO 198 I=1,NS K1=3*(I-1)+1 K2=K1+1 K3=K2+1 B(4,K1)=SFG(1,I) B(4,K3)=SFG(3,I) B(5,K1)=SFG(2,I) B(5,K2)=-SFG(3,I) 198 CONTINUE DO 205 I=1,3 DO 205 J=1,12 205 CBS(I,J)=0.0 DO 199 I=4,5 DO 199 J=1,12 CBS(I,J)=C(4,4)*B(I,J) 199 CONTINUE II =0 DO 196 IA=1,4 JJ=0 DO 195 IB =1,4 DO 193 I=1,2 EKS(II+1,JJ+1)=EKS(II+1,JJ+1)+C(4,4)*SFG(I,IA)*SFG(I,IB)*DB EKS(II+I,JJ+3-1)=EKS(II+1,JJ+3-I)-C(4,4)SFG(I+1,IA) *SFG(4-I,IB)*DB 193 EKS(II+I+1,JJ+I+1)=EKS(II+I+1,JJ+I+1)+C(4,4)SFG(3,IA) *SFG(3,IB)DB DO 194 I=1,3,2 194 EKS(II+I,JJ+4-I)=EKS(II+I,JJ+4-I)+C(4,4)SFG(I,IA) SFG(4-I,IB)DB 195 JJ= JJ+3 196 II =II+3 C C INITIALIZE BENDING STIFFNES MATRIX EKB(22 INTEGRATIOB USED) DO 30 I=1,12 DO 30 J=1,12 30 EKB(I,J)=0.0 C ENTER LOOPS FOR NUMERICAL INTEGRATION NSP=4 DO 170 IX=1,NSP C EVALUATE SHAPE FUNCTIONS ITS DERIVATIVES,JACOBIAN C AND DETERMINANT JACOBIAN GP=0.5773502691896 DO 50 I=1,4 AA=(1.0+R(I)(U(IX)GP)) BB=(1.0+S(I)(V(IX)GP)) SF(1,I)=0.25R(I)BB SF(2,I)=0.25S(I)AA SF(3,I)=0.25AABB 50 CONTINUE C COMPUTE JACOBIAN MATRIX XJAC DO 60 I=1,2 DO 60 J=1,2 60 XJAC (I,J)=0.0 DO 70 I=1,4 DO 70 K=1,2 XJAC(K,1)=XJAC(K,1)+XL(I)*SF(K,I) XJAC(K,2)=XJAC(K,2)+YL(I)*SF(K,I) 70 CONTINUE C CALCULATE DETERMINANT JACOBIAN DJAC DJAC =XJAC(1,1)*XJAC(2,2)-XJAC(1,2)XJAC(2,1) C COMPUTE JACOBIAN INVERSE 'XJACI' XJACI(1,1)=XJAC(2,2)/DJAC XJACI(1,2)=-XJAC(1,2)/DJAC XJACI(2,1)=-XJAC(2,1)/DJAC XJACI(2,2)=XJAC(1,1)/DJAC DA=DJACWG(IX) C FROM GLOBAL DERIVATIVES SFG DO 80 I=1,2 DO 80 J=1,4 80 SFG(I,J) =0.0 DO 90 I=1,2 DO 90 K=1,4 DO 90 J=1,2 SFG(I,K)=SFG(I,K)+XJACI(I,J)SF(J,K) 90 CONTINUE DO 92 I=1,4 92 SFG(3,I)=SF(3,I) C COMPUTE CBB BENDING CONTRIBUTION TO STRESS MATRIX CB C COMPUTE STRAIN DISPLACEMENT MATRIX B DO 100 I=1,5 DO 100 J=1,12 100 B(I,J)=0.0 DO 110 I=1,4 K1=3(I-1)+1 K2=K1+1 K3=K2+1 B(1,K3)=SFG(1,I) B(2,K2)=-SFG(2,I) B(3,K3)=SFG(2,I) B(3,K2)=-SFG(1,I) 110 CONTINUE
DO 120 I=1,3
DO 120 J=1,12
CBB(I,J)=C(I,I)*B(I,J)
IF(I.GT.2) GO TO 120
II =3-I
CBB(I,J)=C(I,II)*B(II,J)+CBB(I,J)
120 CONTINUE DO 113 I=4,5 DO 113 J=1,12 113 CBB(I,J)=0.0 C COMPUTE THE STRESS MATRIX 'CB' DO 123 I=1,NS DO 123 J=1,12 123 CB(I,J)=0.0 DO 124 I=1,NS DO 124 J=1,12 IF (IX.EQ.1) GO TO 126 CB(I,J)=CB(I,J)+CBB(I,J) GO TO 124 126 CB(I,J)=CB(I,J)+CBB(I,J)+CBS(I,J) 124 CONTINUE WRITE(ISTRES)((CB(I,J),I=1,NS),J=1,NED) II=0 DO 150 IA=1,4 JJ=0 DO 140 IB =1,4 DO 130 I=2,3 EKB(II+I,JJ+I)=EKB(II+I,JJ+I)+(C(1,1)*SFG(4-I,IA)*SFG(4-I,IB) *+C(3,3)*SFG(I-1,IA)*SFG(I-1,IB))*DA EKB(II+I,JJ+5-I)=EKB(II+I,JJ+5-I)-(C(1,2)SFG(4-I,IA) *SFG(I-1,IB)+C(3,3)*SFG(I-1,IA)*SFG(4-1,IB))DA 130 CONTINUE 140 JJ=JJ+3 150 II=II+3 C C CALCULATE ELEMENT LOAD VECTOR DUE TO UDL DO 810 I=1,4 NN=3I-2 Q(NN)=Q(NN)+SFG(3,I)UDLDA 810 CONTINUE 170 CONTINUE
C COMPUTE ELEMENT STIFFNESS MATRIX EK DO 800 I=1,12 DO 800 J=1,12 EK(I,J)=EK(I,J)+EKB(I,J)+EKS(I,J) 800 CONTINUE DO 820 I=2,12 K= I-1 DO 820 J=1,K 820 EK(I,J)=EK(J,I) DO 860 J=1,NED JJ=ND(J) IF(JJ) 860,860,840 840 PQ(JJ)=PQ(JJ)+Q(J) 860 CONTINUE NUED=NED CALL PASSEM(SK,EK,NDS,ND,NED,NEQ1,NSKY,NUED) 200 CONTINUE 11 FORMAT(I5) 22 FORMAT(6I5,2F10.4,I5) 33 FORMAT(//35X,'NO. OF 4 NODED PLATE ELEMENTS IN THE' *'STRUCTURE....=',I5//) 44 FORMAT(//25X,'ELEMENT',4X,'MATERIAL',5X,'NODAL CONNECTIVITY', *4X,'THICKNESS',9X,'UDL'/28X,'NO.',5X,'GROUP NO.'//) 55 FORMAT(25X,I5,5X,I5,5X,4I5,2X,F10.4,5X,F10.4/) 133 FORMAT(//5X,'NO. OF STRUCTURE NODES EXCEED 100'/5X, *'CHANGE NOS. IN CORRESPONDING DIMENSION STATEMENTS'/5X, *'IN THE SUBROUTINE PLATE4'//) RETURN END