C ******************************************************************** C * * C * FEAST * C * * C * Finite Element Analysis of Spliced Timber * C * * C * * C * Written by Dave Bohnhoff Spring 1986 * C * * C ******************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER(I-N) COMMON IA(2000),A(9000) COMMON /MPOINT/ MPNOD,MPFEXT,MPDISP,MPSTIF,MPEND,MAXREL, * IPEQN,IPLINK,IPMAT,IPTYPE,IPMREL,IPEND,MAXINT, * MPSTEP,MPLOAD,MPMDIS,MPLDIS,MPDIST COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA COMMON /FRMPRP/ NFRMM,FRMP(4,8) COMMON /FFCPRP/ NFFCP,NFFCC,FFCP(6,3),NFFC(2,8),CFFC(3,8) COMMON /GAPPRP/ NGAPS,GAPP(6,8) COMMON /PLTPRP/ NPLTP,NPLTC,NP(20),PLTP(7,4),PLTC(5,20) COMMON /FPCPRP/ NFPCP,NFPCC,NC(2,40),FPCC(9,40),FPCP(13,4) COMMON /GAUSS4/ SF(4,4,4),ETA(4,4,4),XXI(4,4,4),DLOC(4),WT(4) COMMON /GAUSS8/ WGT(8),XII(8) COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT CHARACTER CARD1*80,CARD2*80,DATAIN*30,OUTPUT*30 C C *** Read name of data file, create output file, open files. C IIN=5 IOUT=6 WRITE(*,'(5X,18HINPUT FILE NAME ? ,\)') READ(*,'(A30)')DATAIN OPEN(IIN,FILE=DATAIN) DO 2 I=1,30 IF(DATAIN(I:I).EQ.'.'.OR.DATAIN(I:I).EQ.' ') GOTO 3 2 CONTINUE 3 OUTPUT=DATAIN OUTPUT(I:I+3) = '.FST' WRITE(*,4)OUTPUT 4 FORMAT(5X,'OUTPUT FILE NAME = ',A30) OPEN (IOUT,FILE=OUTPUT,STATUS='NEW') C C *** Initialize variables. C MAXREL=9000 MAXINT=2000 IBUG=1 NDOF=3 NNPE=2 NSD=3 PIR=4.0D0*DATAN(1.0D0)/180.0D0 CALL GAUSSA CALL GAUSSB C C *** Read and write title cards and control data. C READ(IIN,1000)CARD1 READ(IIN,1000)CARD2 WRITE(IOUT,2000)CARD1,CARD2 READ(IIN,*)NUMNOD,NELE,NNL,NDL,NFRMM,NFFCC,NFFCP,NGAPS,NPLTC, * NPLTP,NFPCC,NFPCP WRITE(IOUT,2001) WRITE(IOUT,2002)NUMNOD,NELE,NNL,NDL,NFRMM,NFFCC,NFFCP, * NGAPS,NPLTC,NPLTP,NFPCC,NFPCP NEQ=NUMNOD*NDOF C C *** Initialize memory pointers for arrays IA and A. C MPNOD=1 MPFEXT=MPNOD+NUMNOD*NSD MPSTEP=MPFEXT+NUMNOD*NDOF MPLOAD=MPSTEP+NUMNOD*NDOF MPDIST=MPLOAD+NUMNOD*NDOF MPDISP=MPDIST+NELE*2 MPLDIS=MPDISP+NUMNOD*NDOF MPMDIS=MPLDIS+NUMNOD*NDOF MPSTIF=MPMDIS+NUMNOD*NDOF MPEND=MPSTIF C IPEQN=1 IPLINK=IPEQN+NUMNOD*NDOF IPMAT=IPLINK+NELE*NNPE IPTYPE=IPMAT+NELE IPMREL=IPTYPE+NELE IPNOUT=IPMREL+NELE*2 IPEOUT=IPNOUT+NUMNOD IPEND=IPEOUT+NELE C C *** Check available memory. C CALL MCHECK(MPEND,IPEND,MAXREL,MAXINT) C C *** Complete data input. C CALL INPUT(A(MPNOD),A(MPSTEP),A(MPDIST),IA(IPEQN),IA(IPLINK), * IA(IPMAT),IA(IPTYPE),IA(IPMREL),IA(IPNOUT),IA(IPEOUT),NNL,NDL) C C *** Determine equation bandwidth. C CALL EQBAND(IA(IPEQN),IA(IPLINK)) C C *** Allocate storage for stiffness matrix. C LENGTH=NEQ+(IBAND-1)*(2*NEQ-IBAND)/2 MPEND=MPSTIF+LENGTH WRITE(IOUT,2003)IPEND WRITE(IOUT,2004)MPEND C C *** Check available memory. C CALL MCHECK(MPEND,IPEND,MAXREL,MAXINT) C C *** Clear stiffness. C DO 10 I=1,LENGTH 10 A(MPSTIF+I-1)=0. C C *** Form stiffness matrix. C IF(NELAS.EQ.1) CONFAC=1.0D0 CALL FRMSTF(A(MPNOD),IA(IPEQN),IA(IPLINK),IA(IPMAT),IA(IPTYPE), * IA(IPMREL),A(MPSTIF),A(MPSTEP),A(MPDIST)) CALL FFCSTF(A(MPNOD),IA(IPEQN),IA(IPLINK),IA(IPMAT),IA(IPTYPE), * A(MPSTIF)) CALL FPCSTF(A(MPNOD),IA(IPEQN),IA(IPLINK),IA(IPMAT),IA(IPTYPE), * A(MPSTIF)) CALL PLTSTF(A(MPNOD),IA(IPEQN),IA(IPLINK),IA(IPMAT),IA(IPTYPE), * A(MPSTIF)) C C *** Modify stiffness and load vector to account for prescribed displ. C CALL MODIFY(IA(IPEQN),A(MPSTEP),A(MPSTIF)) C C *** Factorize stiffness. C IOP=0 EPS=1.E-06 MUD=IBAND-1 ZERO=0. CALL MCHB(ZERO,A(MPSTIF),NEQ,0,MUD,IOP,EPS,IER) IF(IER.EQ.0) GO TO 20 C *** Error in factorization. CALL WI('IER-FACT',IER,1) 20 CONTINUE C C *** Transfer external force to displacement array. C DO 30 I=1,NEQ 30 A(MPLDIS+I-1)=A(MPSTEP+I-1) C C *** Solve simultaneous equations for linear-elastic "step" C *** displacements. C IOP=-1 CALL MCHB(A(MPLDIS),A(MPSTIF),NEQ,1,MUD,IOP,EPS,IER) IF(IER.EQ.0) GO TO 40 C *** Error in back substitution. CALL WI('IER-BACK',IER,1) 40 CONTINUE C C *** Loop over the number of load steps. Calculate current linear C *** displacements and loads from step disp. and step load vectors. C DO 50 LS=1,NSTEP WRITE(IOUT,2005)LS DO 60 I=1,NEQ A(MPDISP+I-1)=A(MPLDIS+I-1)*(ILOAD+LS-1) 60 A(MPFEXT+I-1)=A(MPSTEP+I-1)*(ILOAD+LS-1) C C *** Output displacement if only a linear-elastic analysis. C IF(NELAS.EQ.1) GOTO 130 IBETA = 0 C C *** Begin non-linear analysis. C DO 80 ITER=1,NITER WRITE(*,'(3X,I2,2X,1H-,I4)')LS,ITER C C *** Zero and then assemble load modifying vector. C DO 90 I=1,NEQ 90 A(MPLOAD+I-1)=0.D0 IF(NFFCP.GT.0) CALL FFCNLB(A(MPNOD),IA(IPLINK),IA(IPMAT), * IA(IPTYPE),A(MPDISP),A(MPLOAD)) IF(NGAPS.GT.0) CALL GAPNLB(A(MPNOD),IA(IPLINK),IA(IPMAT), * IA(IPTYPE),A(MPDISP),A(MPLOAD),IBETA) IF(NPLTP.GT.0) CALL PLTNLB(A(MPNOD),IA(IPLINK),IA(IPMAT), * IA(IPTYPE),A(MPDISP),A(MPLOAD)) IF(NFPCP.GT.0) CALL FPCNLB(A(MPNOD),IA(IPLINK),IA(IPMAT), * IA(IPTYPE),A(MPDISP),A(MPLOAD)) C C *** Zero loads associated with prescribed DOF's. Add modifying load C *** vector to external load vector and transfer sum to disp. vector. C DO 100 I=1,NEQ IF(IA(IPEQN+I-1).LE.0) A(MPLOAD+I-1)=0.D0 100 A(MPMDIS+I-1)=A(MPLOAD+I-1)+A(MPFEXT+I-1) C C *** Solve simultaneous equations. C IOP=-1 CALL MCHB(A(MPMDIS),A(MPSTIF),NEQ,1,MUD,IOP,EPS,IER) IF(IER.EQ.0) GOTO 110 C *** Error in back substitution CALL WI('IER-BACK',IER,1) 110 CONTINUE C C *** Check convergence. C CALL CONVRG(A(MPMDIS),A(MPDISP),CHECK) IF(CHECK.LE.0.D0) THEN DO 125 I=1,NEQ 125 A(MPDISP+I-1)=A(MPMDIS+I-1) WRITE(IOUT,2030)ITER GOTO 130 END IF IF(ITER.EQ.NITER) THEN WRITE(IOUT,2035)ITER GOTO 200 END IF C C *** Update displacements. Use underrelaxation if one or more gaps C *** have closed (IBETA = 1). C IF(IBETA.EQ.1) THEN DO 128 I=1,NEQ 128 A(MPDISP+I-1) = (1.-BETA)*A(MPDISP+I-1)+BETA*A(MPMDIS+I-1) ELSE DO 129 I=1,NEQ 129 A(MPDISP+I-1) = A(MPMDIS+I-1) END IF C 80 CONTINUE C C *** Output nodal displacements. C 130 IF(NNOUT.EQ.0) GOTO 145 WRITE(IOUT,2040) DO 140 I=1,NUMNOD IF(IA(IPNOUT+I-1).NE.1) GOTO 140 WRITE(IOUT,2050)I,A(MPDISP+I*3-3),A(MPDISP+I*3-2),A(MPDISP+I*3-1) 140 CONTINUE C C *** Post-process data. C 145 IF(NEOUT.EQ.0) GOTO 50 CALL FRMFOR(A(MPNOD),IA(IPLINK),IA(IPMAT),IA(IPTYPE),IA(IPMREL), * IA(IPEOUT),A(MPDISP),A(MPDIST),ILOAD,LS) IF(NFFCP.EQ.0) GOTO 50 CALL FFCFOR(A(MPNOD),IA(IPLINK),IA(IPMAT),IA(IPTYPE),IA(IPEOUT), * A(MPDISP)) C 50 CONTINUE 200 STOP 1000 FORMAT(A80) 1001 FORMAT(6I5) 2000 FORMAT(/,5X,1H ,A80,/,5X,1H ,A80) 2001 FORMAT(//,5X,21H PROGRAM CONTROL DATA,/,6X,64(1H-)) 2002 FORMAT(6X,15HNumber of nodes,45X,1H=,I3,/,6X, * 18HNumber of elements,42X,1H=,I3,/,6X, * 39HNumber of nodes with either an applied ,/,24X, * 39Hload or prescribed nonzero displacement,3X,1H=,I3,/,6X, * 57HNumber of frame elements with uniformly distributed loads, * 3X,1H=,I3,/,6X, * 43HNumber of different frame element materials, * 17X,1H=,I3,/,6X, * 57HNumber of frame-to-frame connector element configurations, * 3X,1H=,I3,/,6X, * 43HNumber of sets of load-slip parameters for ,/,38X, * 25Hframe-to-frame connectors,3X,1H=,I3,/,6X * 36HNumber of gap element specifications, * 24X,1H=,I3,/,6X, * 38HNumber of plate element configurations,22X,1H=,I3,/,6X, * 40HNumber of sets of PLT element properties,20X,1H=,I3,/,6X, * 57HNumber of frame-to-plate connector element configurations, * 3X,1H=,I3,/,6X, * 55HNumber of sets of load-slip parameters for FPC elements, * 5X,1H=,I3) 2003 FORMAT(6X,36HLength of integer storage array {IA},14X,1H=,I6) 2004 FORMAT(6X,32HLength of real storage array {A},18X,1H=,I6) 2005 FORMAT(///,6X,10HLOAD STEP ,I2,/,6X,12(1H-)) 2030 FORMAT(/,6X,35HNUMBER OF ITERATIONS TO CONVERGE = ,I3) 2035 FORMAT(/,6X,28HCONVERGENCE NOT ACHIEVED IN ,I3,11H ITERATIONS) 2040 FORMAT(/,6X,49HNODAL DISPLACEMENT VALUES (IN GLOBAL COORDINATES), * /,6X,49(1H-),/,6X,4HNode,5X,14HX Displacement,4X, * 14HY Displacement,7X,10HZ Rotation,/,6X,58(1H-)) 2050 FORMAT(7X,I3,3X,F15.7,3X,F15.7,3X,F15.7) END SUBROUTINE MCHECK(MPEND,IPEND,MAXREL,MAXINT) IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER(I-N) COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT IF(MPEND.GT.MAXREL .OR. IPEND.GT.MAXINT) WRITE(IOUT,100) RETURN 100 FORMAT(//,5X,34H MEMORY ERROR: INSUFFICIENT MEMORY /) END SUBROUTINE EQBAND(KFIX,LINK) IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER(I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR DIMENSION KFIX(NDOF,1),LINK(NNPE,1) C C *** Scan elements to determine the maximum eqn difference of C *** unconstrained nodes - then compute half bandwidth. C IBAND=0 DO 20 I=1,NELE MEQMIN=10000 MEQMAX=0 DO 10 J=1,NNPE NODE=LINK(J,I) DO 10 K=1,NDOF IEQNUM=(NODE-1)*NDOF+K IF(IEQNUM.LT.MEQMIN .AND. KFIX(K,NODE).NE.0) MEQMIN=IEQNUM IF(IEQNUM.GT.MEQMAX .AND. KFIX(K,NODE).NE.0) MEQMAX=IEQNUM 10 CONTINUE IF(MEQMIN.EQ.10000) GO TO 20 IF(MEQMAX-MEQMIN+1.GT.IBAND) IBAND=MEQMAX-MEQMIN+1 20 CONTINUE RETURN END SUBROUTINE INPUT(X,F,DF,KFIX,LINK,MAT,NTYPE,MREL,NOUT,LOUT,NL,ND) IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER(I-N) CHARACTER TYP(5)*11,CHK(2)*3 COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA COMMON /FRMPRP/ NFRMM,FRMP(4,8) COMMON /FFCPRP/ NFFCP,NFFCC,FFCP(6,3),NFFC(2,8),CFFC(3,8) COMMON /GAPPRP/ NGAPS,GAPP(6,8) COMMON /PLTPRP/ NPLTP,NPLTC,NP(20),PLTP(7,4),PLTC(5,20) COMMON /FPCPRP/ NFPCP,NFPCC,NC(2,40),FPCC(9,40),FPCP(13,4) COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT DIMENSION X(NSD,1),F(NDOF,1),KFIX(NDOF,1),LINK(NNPE,1),MAT(1), * Z1(3),Z2(3),IZ(3),NTYPE(1),DF(2,1),MREL(2,1),NOUT(1),LOUT(1) DATA TYP/'FRM - Matl#','FFC - Conf#','GAP - Spec#','FPC - Conf#', * 'PLT - Conf#'/,CHK/'No ','Yes'/ C C *** Read and write nodal data. C WRITE(IOUT,2000) DO 25 I=1,NUMNOD READ(IIN,*)N,(Z1(J),J=1,3),(IZ(J),J=1,3) DO 10 J=1,NSD 10 X(J,I)=Z1(J) DO 20 J=1,NDOF IF(IZ(J).EQ.0) KFIX(J,I)=1 IF(IZ(J).EQ.1) KFIX(J,I)=0 20 IF(IZ(J).EQ.2) KFIX(J,I)=-1 25 WRITE(IOUT,2010)N,(Z1(J),J=1,3),(IZ(J),J=1,3) C C *** Read and write element data. C WRITE(IOUT,2020) DO 30 I=1,NELE READ(IIN,*)N,NTYPE(I),MAT(I),(LINK(J,I),J=1,NNPE) WRITE(IOUT,2030)I,TYP(NTYPE(I)),MAT(I),(LINK(J,I),J=1,NNPE) IF(NTYPE(I).NE.1) GOTO 30 READ(IIN,*)MREL(1,I),MREL(2,I) IF(MREL(1,I).EQ.0.AND.MREL(2,I).EQ.0)THEN GOTO 30 ELSE IF(MREL(1,I)+MREL(2,I).GT.1.9)THEN WRITE(IOUT,2040)I ELSE IF(MREL(1,I).EQ.1)THEN WRITE(IOUT,2050)I,LINK(1,I) ELSE WRITE(IOUT,2050)I,LINK(2,I) END IF 30 CONTINUE C C *** Read and write nodal loads and/or prescribed displacements. C DO 31 I=1,NUMNOD DO 31 J=1,NDOF 31 F(J,I)=0.D0 IF(NL.EQ.0) GOTO 36 WRITE(IOUT,2060) DO 35 I=1,NL READ(IIN,*)N,(Z2(J),J=1,3) DO 34 J=1,NDOF 34 F(J,N)=Z2(J) 35 WRITE(IOUT,2070)N,(Z2(J),J=1,3) C C *** Read and write data on uniformly distributed loads. C 36 DO 37 I=1,NELE DF(1,I)=0.D0 37 DF(2,I)=0.D0 IF(ND.EQ.0) GOTO 39 WRITE(IOUT,2080) DO 38 I=1,ND READ(IIN,*)N,DFX,DFY DF(1,N)=DFX DF(2,N)=DFY 38 WRITE(IOUT,2090)N,DFX,DFY C C *** Read and write data on frame element materials. C 39 WRITE(IOUT,2100) DO 50 I=1,NFRMM READ(IIN,*)N,(FRMP(J,I),J=1,4) 50 WRITE(IOUT,2110)I,(FRMP(J,I),J=1,4) C C *** Read/write data on frame-to-frame connector elmt configurations. C IF(NFFCP.EQ.0) GOTO 75 WRITE(IOUT,2120) DO 60 I=1,NFFCC READ(IIN,*)N,NFFC(1,I),NFFC(2,I),(CFFC(J,I),J=1,3) 60 WRITE(IOUT,2130)N,NFFC(1,I),NFFC(2,I),(CFFC(J,I),J=1,NFFC(2,I)) C C *** Read/write load-slip parameters for frame-to-frame connectors. C WRITE(IOUT,2140) DO 70 I=1,NFFCP READ(IIN,*)N,(FFCP(J,I),J=1,6) 70 WRITE(IOUT,2150)I,(FFCP(J,I),J=1,6) C C *** Read/write data on gaps. C 75 IF(NGAPS.EQ.0)GOTO 90 WRITE(IOUT,2160) DO 80 I=1,NGAPS READ(IIN,*)N,(GAPP(J,I),J=1,6) 80 WRITE(IOUT,2170)I,(GAPP(J,I),J=1,6) C C *** Read and write data on plate element configurations. C 90 IF(NPLTP.EQ.0)GOTO 115 WRITE(IOUT,2180) DO 100 I=1,NPLTC READ(IIN,*)N,NP(I),(PLTC(J,I),J=1,5) 100 WRITE(IOUT,2190)I,NP(I),(PLTC(J,I),J=1,5) C C *** Read and write data on plate element properties. C WRITE(IOUT,2200) DO 110 I=1,NPLTP READ(IIN,*)N,(PLTP(J,I),J=1,7) 110 WRITE(IOUT,2210)I,(PLTP(J,I),J=1,7) C C *** Read/write frame-to-plate connector element configurations. C 115 IF(NFPCP.EQ.0)GOTO 150 WRITE(IOUT,2220) DO 120 I=1,NFPCC READ(IIN,*)N,NC(1,I),NC(2,I),(FPCC(J,I),J=1,9) WRITE(IOUT,2230)I,NC(1,I),NC(2,I),(FPCC(J,I),J=1,5) 120 WRITE(IOUT,2240)(FPCC(J,I),J=6,9) C C *** Read/write load-slip parameters for frame-to-plate elements. C WRITE(IOUT,2250) DO 130 I=1,NFPCP READ(IIN,*)N,(FPCP(J,I),J=1,7) READ(IIN,*)(FPCP(J,I),J=8,13) 130 WRITE(IOUT,2260)I,(FPCP(J,I),J=1,13) C C *** Read information on desired output. C 150 READ(IIN,*)NNOUT IF(NNOUT.EQ.NUMNOD) THEN DO 160 I=1,NUMNOD 160 NOUT(I) = 1 ELSE DO 170 I=1,NUMNOD 170 NOUT(I) = 0 DO 180 I=1,NNOUT READ(IIN,*) NN 180 NOUT(NN)=1 END IF READ(IIN,*)NEOUT IF(NEOUT.EQ.NELE) THEN DO 190 I=1,NELE 190 LOUT(I) = 1 ELSE DO 200 I=1,NELE 200 LOUT(I) = 0 DO 210 I=1,NEOUT READ(IIN,*) NN 210 LOUT(NN) = 1 END IF C C *** Read and write additional control data. C WRITE(IOUT,2270) READ(IIN,*)ILOAD,NSTEP,NELAS IF(NELAS.NE.1) NELAS=2 WRITE(IOUT,2280)ILOAD,NSTEP,CHK(NELAS) IF(NELAS.EQ.1) GOTO 220 READ(IIN,*)NITER,TOL,CONFAC,BETA WRITE(IOUT,2290)NITER,TOL,CONFAC,BETA C 220 CONTINUE RETURN C 2000 FORMAT(//,5X,' NODAL DATA',/,6X,10(1H-),/,5X,' Node # ', * 'Coordinates (X,Y,Z) Fixity (X,Y,Mz)',/,6X,55(1H-)) 2010 FORMAT(6X,I4,3X,3F10.4,3X,3I4) 2020 FORMAT(//,5X,' ELEMENT DATA',/,6X,12(1H-),/,5X,' Element# ', * ' Type Node 1 Node 2',/,6X,47(1H-)) 2030 FORMAT(2X,I9,7X,A11,I3,2I9) 2040 FORMAT(35X,'(both ends of member ',I2,' are pinned)') 2050 FORMAT(35X,'(member ',I2,' pinned at node ',I2,')') 2060 FORMAT(//,5X,' NODAL LOADS AND/OR PRESCRIBED DISPLACEMENTS',/,6X, * 43(1H-),/,5X,' Node# Forces/Displacements (X,Y,Mz)', * /,6X,49(1H-)) 2070 FORMAT(6X,I4,3F15.4) 2080 FORMAT(//,5X,' UNIFORM DISTRIBUTED LOADS',/,6X,25(1H-),/,5X, * ' Element # Force/Unit Length (X,Y) ',/,6X,38(1H-)) 2090 FORMAT(6X,I5,3X,2F15.4) 2100 FORMAT(//,6X,37HFRAME ELEMENT MATERIAL SPECIFICATIONS,/,6X, * 37(1H-),/,6X,5HMatl#, * 5X,4HArea,11X,1HI,12X,3HMOE,10X,5HDepth,/,6X,57(1H-)) 2110 FORMAT(5X,I4,2X,F11.4,2X,F11.4,2X,F15.4,2X,F9.4) 2120 FORMAT(//,6X,47HFRAME-TO-FRAME CONNECTOR ELEMENT CONFIGURATIONS, * ,/,6X,47(1H-), * /,5X,' Conf# LSP# # and location of connectors',/, * 6X,46(1H-)) 2130 FORMAT(1X,3I8,6F10.5) 2140 FORMAT(//,5X,' LOAD-SLIP PARAMETERS FOR FRAME-TO-FRAME CONNECTOR', * ' ELEMENTS',/,6X,58(1H-),/,5X,' LSP# Parallel to', * ' Grain Perpendicular to Grain',/,16X,' P0 ', * ' P1 K P0 P1 K',/,6X,65(1H-)) 2150 FORMAT(5X,I4,2X,1P6E10.3) 2160 FORMAT(//,6X,26HGAP ELEMENT SPECIFICATIONS,/,6X,26(1H-),/,5X, * ' Spec# fc Width DT(1) DB(1) DT(2) DB(2)', * /,6X,55(1H-)) 2170 FORMAT(4X,I5,3X,F8.2,1X,5(2X,F6.3)) 2180 FORMAT(//,6X,28HPLATE ELEMENT CONFIGURATIONS,/,6X,28(1H-),/,5X, * ' Conf# Prp# Angle',12X,'Coordinates of the Corners', * /,33X,'X(1) Y(1) X(2) Y(2)',/,6X,65(1H-)) 2190 FORMAT(6X,I3,4X,I3,3X,F7.2,3X,4(F9.3,2X)) 2200 FORMAT(//,5X,' PLATE ELEMENT PROPERTIES',/,6X,24(1H-),/,5X, * ' Prp# Density K,axial K,shear Pc Pt ', * ' Ps Length',/,6X,67(1H-)) 2210 FORMAT(I8,3X,F7.4,2E10.3,3(2X,F7.1),2X,F6.3) 2220 FORMAT(//,6X,47HFRAME-TO-PLATE CONNECTOR ELEMENT CONFIGURATIONS, * /,6X,47(1H-),/,6X,26HConf# LSP# Node 3 Angle,10X, * 'Coordinates of the Corners',/,6X,69(1H-)) 2230 FORMAT(6X,I3,4X,I3,4X,I3,2X,F7.2,3X,4HX(I),4(1X,F8.2)) 2240 FORMAT(35X,4HY(I),4(1X,F8.2)) 2250 FORMAT(//,6X,40HLOAD-SLIP PARAMETERS FOR FRAME-TO-PLATE , * 18HCONNECTOR ELEMENTS,/,6X,58(1H-),/,5X,' LSP# Density', * ' Plate Grain',8X,1HK,12X,2HP0,12X,2HP1,/,6X,69(1H-)) 2260 FORMAT(6X,I3,4X,F5.2,4X,'Par Par ',3(2X,1PE12.5),/,22X, * 'Per Par ',3(2X,1PE12.5),/,22X,'Par Per ',3(2X, * 1PE12.5),/,22X,'Per Per ',3(2X,1PE12.5)) 2270 FORMAT(//,5X,21H PROGRAM CONTROL DATA,/,6X,57(1H-)) 2280 FORMAT(6X,48HNumber of load steps comprising the initial load,2X, * 1H=,I6,/,6X,20HNumber of load steps,30X,1H=,I6,/,6X, * 19HNon-linear analysis,31X,1H=,3X,A3) 2290 FORMAT(6X,38HMaximum number of iterations/load step,12X,1H=,I6,/, * 6X,46HTolerance for convergence of iterative process,4X,1H=, * 1X,1PE8.2,/,6X,25HFactor for increasing the,/,20X, * 34Hinitial stiffness of the structure,2X,1H=,1X,1PE8.2,/,6X, * 22HUnderrelaxation factor,28X,1H=,1X,1PE8.2) END SUBROUTINE GAUSSA C C *** This subroutine calculates values for shape functions (of a C *** linear element) and their derivatives at each of 16 gauss C *** points (4-by-4 gauss rule). C C *** Description of subroutine variables. C *** DLOC(L) Location of sampling points. C *** WT(L) Weights for gauss quadrature. C *** SF(N,I,J) Value of shape function N at Xi=DLOC(I) and C *** Eta=DLOC(J). C *** ETA(N,I,J) Value of the derivative of shape function N with C *** respect to Eta at Xi=DLOC(I) and Eta=DLOC(J). C *** XXI(N,I,J) Value of the derivative of shape function N with C *** respect to Xi at Xi=DLOC(I) and Eta=DLOC(J). C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER(I-N) COMMON /GAUSS4/ SF(4,4,4),ETA(4,4,4),XXI(4,4,4),DLOC(4),WT(4) C DLOC(1) = -0.861136311594053 DLOC(2) = -0.339981043584856 DLOC(3) = -DLOC(2) DLOC(4) = -DLOC(1) WT(1) = 0.347854845137454 WT(2) = 0.652145154862546 WT(3) = WT(2) WT(4) = WT(1) DO 10 I=1,4 DUM1 = (1.0D0+DLOC(I))/4.0D0 DUM2 = (1.0D0-DLOC(I))/4.0D0 DO 10 J=1,4 DUM3 = 1.0D0+DLOC(J) DUM4 = 1.0D0-DLOC(J) ETA(1,I,J) =-DUM4/4.0D0 ETA(2,I,J) =-ETA(1,I,J) ETA(3,I,J) = DUM3/4.0D0 ETA(4,I,J) =-ETA(3,I,J) XXI(1,I,J) =-DUM2 XXI(2,I,J) =-DUM1 XXI(3,I,J) = DUM1 XXI(4,I,J) = DUM2 SF(1,I,J) = DUM2*DUM4 SF(2,I,J) = DUM1*DUM4 SF(3,I,J) = DUM1*DUM3 10 SF(4,I,J) = DUM2*DUM3 RETURN END SUBROUTINE GAUSSB C C *** Location and weights for an 8th order Gauss quadrature. C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER(I-N) COMMON /GAUSS8/ WGT(8),XII(8) C XII(1) = -0.960289856497536 XII(2) = -0.796666477413627 XII(3) = -0.525532409916329 XII(4) = -0.183434642495650 WGT(1) = 0.101228536290376 WGT(2) = 0.222381034453374 WGT(3) = 0.313706645877887 WGT(4) = 0.362683783378362 DO 10 I=1,4 XII(4+I) = -XII(5-I) 10 WGT(4+I) = WGT(5-I) RETURN END SUBROUTINE ADSTIF(KFIX,S,LINK,T,ITROW,KELE) IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR DIMENSION KFIX(1),S(1),T(ITROW,1),LINK(NNPE,1) C C *** Assemble element stiffness into upper triangle of global stiff. C DO 100 IGEN=1,NNPE DO 100 JGEN=1,NNPE IGLB=NDOF*(LINK(IGEN,KELE)-1) JGLB=NDOF*(LINK(JGEN,KELE)-1) C DO 100 IDOF=1,NDOF DO 100 JDOF=1,NDOF C C *** If lower triangular element, skip assembly C IF(IGLB+IDOF.GT.JGLB+JDOF) GO TO 100 C C *** Enforce boundary condition: C *** If constrained DOF and nondiagonal term, skip assembly. C IF((KFIX(IGLB+IDOF).EQ.0 .OR. KFIX(JGLB+JDOF).EQ.0) .AND. * IGLB+IDOF.NE.JGLB+JDOF) GO TO 100 C JADD=LLOC(IGLB+IDOF,JGLB+JDOF) S(JADD)=S(JADD)+T(NDOF*(IGEN-1)+IDOF,NDOF*(JGEN-1)+JDOF) 100 CONTINUE RETURN END SUBROUTINE MODIFY(KFIX,F,S) IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER(I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT DIMENSION KFIX(1),F(1),S(1) C C *** Subroutine to modify stiffness matrix and load vector to C *** account for prescribed displacements. C C *** Make sure force is zero for each prescribed zero D.O.F. C *** Setting zero diagonal values to one. DO 5 MEQ=1,NEQ IF (KFIX(MEQ) .EQ. 0) F(MEQ) = 0.D0 IF(S(LLOC(MEQ,MEQ)).GT.0.00001) GOTO 5 S(LLOC(MEQ,MEQ))=1.0 WRITE(IOUT,1000)MEQ 1000 FORMAT(/,5X,' Diagonal value associated with equation ',I3, *' set equal to one.') 5 CONTINUE C DO 100 MEQ=1,NEQ IF(KFIX(MEQ).GE.0) GO TO 100 C C *** First: Modify load vector for equation MEQ. C MSTART=MEQ-IBAND+1 IF(MSTART.LT.1) MSTART=1 MSTOP=MEQ+IBAND-1 IF(MSTOP.GT.NEQ) MSTOP=NEQ C DO 10 K=MSTART,MSTOP IF(K.EQ.MEQ) GO TO 10 IF(K.LT.MEQ) KLOC=LLOC(K,MEQ) IF(K.GT.MEQ) KLOC=LLOC(MEQ,K) C C *** Skip modification if equation K is also a prescribed displ. C IF(KFIX(K).GE.0) F(K)=F(K)-S(KLOC)*F(MEQ) 10 CONTINUE C C *** Second: Modify stiffness - zero row and col and leave diag. C DO 30 K=MSTART,MSTOP IF(K.NE.MEQ) GO TO 20 C C *** Diagonal term- modify prescribed displ and leave diag stiff. C F(MEQ)=F(MEQ)*S(LLOC(MEQ,MEQ)) GO TO 30 C C *** Off diagonal term - zero stiffness. C 20 IF(K.GT.MEQ) S(LLOC(MEQ,K))=0. IF(K.LT.MEQ) S(LLOC(K,MEQ))=0. 30 CONTINUE C 100 CONTINUE RETURN END FUNCTION LLOC(I,J) IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR C C *** Find location in 1-D array of I,J address in 2-D array C *** (location is based on constant bandwidth upper triangular C *** column-wise storage). C LLOC=IBAND*(I-1)+J-I+1 C C *** Adjust address if in last IBAND-2 rows. C IADJST=NEQ-IBAND+2 IF(I.LE.IADJST) RETURN NHT=I-IADJST LLOC=LLOC-NHT*(NHT+1)/2 RETURN END SUBROUTINE CONVRG(DISM,DISP,CHECK) C C *** Convergence can be considered achieved when CHECK is a C *** negative number. C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER(I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA DIMENSION DISM(1),DISP(1) C EI=0.D0 SD=0.D0 DO 10 I=1,NEQ EI=EI+(DISM(I)-DISP(I))**2 10 SD=SD+DISP(I)**2 EI=SQRT(EI) SD=SQRT(SD)*TOL CHECK=EI-SD RETURN END SUBROUTINE WI(ID,I,IOP) IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT CHARACTER*8 ID IF(IOP.GT.IBUG) RETURN WRITE(IOUT,100)ID,I 100 FORMAT(6X,A8,1H=,I10) RETURN END SUBROUTINE WR(ID,R,IOP) IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT CHARACTER*8 ID IF(IOP.GT.IBUG) RETURN WRITE(IOUT,100)ID,R 100 FORMAT(6X,A8,1H=,G15.5) RETURN END SUBROUTINE MCHB (R,A,M,N,MUD,IOP,EPS,IER) C C *** For a given positive definite M by M matrix A with symmetric band C *** structure and (if necessary) a given general M by N matrix R, the C *** following calculations (dependent on the value of the decision C *** parameter IOP) are performed. C *** (1) Matrix A is factorized (If IOP is not negative), that means C *** band matrix TU with upper codiagonals only is generated on the C *** locations of A such that TRANSPOSE(TU)*TU = A. C *** (2) Matrix R is multiplied on the left by INVERSE(TU) and/or C *** INVERSE(TRANSPOSE(TU)) and the result is stored in the C *** location of R. C *** This subroutine especially can be used to solve the system of C *** simultaneous linear equations A*X=R with positive-definite C *** coefficient matrix A of symmetric band structure. C C *** Description of Parameters. C *** R Input IOP = -3,-2,-1,1,2,3 M by N right hand matrix, C *** IOP = 0 irrelevant. C *** Output IOP = 1,-1 INVERSE(A)*R, C *** IOP = 2,-2 INVERSE(TU)*R, C *** IOP = 3,-3 INVERSE(TRANSPOSE(TU))*R, C *** IOP = 0 unchanged. C *** A Input IOP = 0,1,2,3 M x M positive-definite coefficient C *** matrix of symmetric band structure stored in C *** compressed form (see remarks), C *** IOP = -1,-2,-3 M x M band matrix TU with upper C *** codiagonals only, stored in compressed form. C *** Output in all cases band matrix TU with upper codiagonals C *** only, stored in compressed form (that means C *** unchanged if IOP = -1,-2,-3). C *** M Input value specifying the number of rows and columns of A C *** and number of rows of R. C *** N Input value specifying the number of columns of R C *** (irrelevant in case IOP = 0). C *** MUD Input value specifying the no. of upper codiagonals of A. C *** IOP Decision parameter (either -3,-2,-1,0,1,2,3). C *** EPS Input value used as relative tolerance for test on loss of C *** significant digits. C *** IER Resulting error parameter coded as follows, C *** IER = 0 - No error, C *** IER = -1 - No result because of wrong input parameters M, C *** MUD, IOP, or because of nonpositive radicand at C *** some factorization step, or because of zero C *** diagonal element at some division step. C *** IER = K - Warning due to possible loss of significance C *** indicated at factorization step K+1 where radi- C *** cand was no longer greater than EPS*A(K+1,K+1). C C *** Remarks. C *** Upper part of symmetric band matrix A consisting of main diagonal C *** and MUD upper codiagonals (resp. band matrix TU consisting of main C *** diagonal and MUD upper codiagonals) is assumed to be stored in C *** compressed form, i.e. rowwise in totally needed M+MUD*(2M-MUD-1)/2 C *** successive storage locations. On return upper band factor TU (on C *** the locations of A) is stored in the same way. Right hand side C *** matrix R is assumed to be stored columnwise in N*M successive C *** storage locations. On return, INVERSE(A)*R or INVERSE(TU)*R or C *** INVERSE(TRANSPOSE(TU))*R will be stored in the locations of R. C *** Input parameters M, MUD, IOP should satisfy the following C *** restrictions, 1) MUD not less than zero, 2) 1+MUD not greater than C *** M, 3) ABS(IOP) not greater than 3. No action besides error mes- C *** sage IER=-1 takes place if these restrictions are not satisfied. C *** The procedure gives results if the restrictions on input para- C *** meters are satisfied, if radicands at all factorization steps are C *** positive and/or if all diagonal elements of upper band factor TU C *** are nonzero. Factorization is done using Cholesky-S square-root C *** method, which generates the upper band matrix TU such that C *** TRANSPOSE(TU)*TU=A. C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) DIMENSION R(1),A(1) C C *** Test on wrong input parameters. C IF(IABS(IOP)-3)100,100,520 100 IF(MUD)520,110,110 110 MC=MUD+1 IF(M-MC)520,120,120 120 MR=M-MUD IER=0 C C *** MC is the maximum number of elements in the rows of array A. C *** MR is the index of the last row in array A with MC elements. C C C *** Start factorization of matrix A. C IF(IOP)330,130,130 130 IEND=0 LLDST=MUD DO 320 K=1,M IST=IEND+1 IEND=IST+MUD J=K-MR IF(J)150,150,140 140 IEND=IEND-J 150 IF(J-1)170,170,160 160 LLDST=LLDST-1 170 LMAX=MUD J=MC-K IF(J)190,190,180 180 LMAX=LMAX-J 190 ID=0 TOL=A(IST)*EPS C C *** Start factorization-loop over K-th row. C DO 320 I=IST,IEND SUM=0.000 IF(LMAX)230,230,200 C C *** Prepare inner loop. C 200 LL=IST LLD=LLDST C C *** Start inner loop. C DO 220 L=1,LMAX LL=LL-LLD LLL=LL+ID SUM=SUM+A(LL)*A(LLL) IF(LLD-MUD)210,220,220 210 LLD=LLD+1 220 CONTINUE C C *** End of inner loop. C C C *** Transform element A(I). C 230 SUM=A(I)-SUM IF(I-IST)240,240,290 C C *** A(I) is diagonal element. Error test. C 240 IF(SUM)540,540,250 C C *** Test on loss of significant digits and warnings. C 250 IF(SUM-TOL)260,260,280 260 IF(IER)270,270,280 270 IER=K-1 C C *** Computation of pivot element. C 280 PIV= SQRT(SUM) A(I)=PIV PIV= 1.0/PIV GO TO 300 C C *** A(I) is not diagonal element. C 290 A(I)=SUM*PIV C C *** Update ID and LMAX. C 300 ID=ID+1 IF(ID-J)320,320,310 310 LMAX=LMAX-1 320 CONTINUE C C *** End of factorization-loop over K-th row. C *** End of factorization of matrix A. C C *** Prepare matrix divisions. C IF(IOP)330,530,330 330 ID=N*M IEND=IABS(IOP)-2 IF(IEND)340,440,340 C C *** Start division by transpose of matrix TU (TU is stored in C *** locations of A). C 340 IST=1 LMAX =0 J =-MR LLDST =MUD DO 430 K=1,M PIV = A(IST) IF(PIV)350,560,350 350 PIV= 1.0/PIV C C *** Start backsubstitution-loop for K-th row of matrix R. C DO 390 I=K,ID,M SUM=0.000 IF(LMAX)390,390,360 C C *** Prepare inner loop. C 360 LL=IST LLL=I LLD=LLDST C C *** Start inner loop. C DO 380 L=1,LMAX LL=LL-LLD LLL=LLL-1 SUM=SUM+A(LL)*R(LLL) IF(LLD-MUD)370,380,380 370 LLD=LLD+1 380 CONTINUE C C *** End of inner loop. C C *** Transform element R(I). C 390 R(I)=PIV*(R(I)-SUM) C C *** End of backsubstitution-loop for K-th row of matrix R. C C *** Update parameters LMAX, IST AND LLDST. C IF(MC-K)410,410,400 400 LMAX=K 410 IST=IST+MC J=J+1 IF(J)430,430,420 420 IST=IST-J LLDST=LLDST-1 430 CONTINUE C C *** End of division by transpose of matrix TU. C *** Start division by matrix TU (TU is stored on locations of A). C IF(IEND)440,440,530 440 IST=M+(MUD*(M+M-MC))/2+1 LMAX=0 K=M 450 IEND=IST-1 IST=IEND-LMAX PIV=A(IST) IF(PIV)460,520,460 460 PIV= 1.0/PIV L=IST+1 C C *** Start backsubstitution-loop for K-th row of matrix R. C DO 490 I=K,ID,M SUM=0.000 IF(LMAX)490,490,470 470 LLL=I C C *** Start inner loop. C DO 480 LL=L,IEND LLL=LLL+1 480 SUM=SUM+A(LL)*R(LLL) C C *** End of inner loop. C C *** Transform element R(I). C 490 R(I)=PIV*(R(I)-SUM) C C *** End of backsubstitution-loop for K-th row of matrix R. C C *** Update parameters LMAX and K. C IF(K-MR)510,510,500 500 LMAX=LMAX+1 510 K=K-1 IF(K)530,530,450 C C *** End of division by matrix TU. C C *** Error exit in case of wrong input parameters or pivot element C *** less than or equal to zero. C 520 IER=-1 530 RETURN 540 IER=-2 GO TO 530 560 IER=-3 GO TO 530 END SUBROUTINE PLTSTF(XX,KFIX,LINK,MATNUM,NTYPE,S) IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /PLTPRP/ NPLTP,NPLTC,NP(20),PLTP(7,4),PLTC(5,20) COMMON /GAUSS8/ WGT(8),XII(8) COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA DIMENSION XX(NSD,1),KFIX(1),LINK(NNPE,1),MATNUM(1),NTYPE(1), * S(1),X(2),Y(2),ESM(6,6) C C *** Loop over elements. Check to see if plate type element. C DO 100 NN=1,NELE IF(NTYPE(NN).NE.5) GOTO 100 MAT = MATNUM(NN) NT = NP(MAT) C C *** Set nodal coordinates and calculate plate angle relationships. C DO 20 L=1,2 J = LINK(L,NN) X(L) = XX(1,J) 20 Y(L) = XX(2,J) CA = COS(PIR*PLTC(1,MAT)) SA = SIN(PIR*PLTC(1,MAT)) CC = CA*CA SS = SA*SA CS = CA*SA C DL = SQRT((PLTC(4,MAT)-PLTC(2,MAT))**2 + * (PLTC(5,MAT)-PLTC(3,MAT))**2) FAC = PLTP(1,NT)*DL/2.0D0 AK = PLTP(2,NT) SK = PLTP(3,NT) DO 30 I=1,6 DO 30 J=1,6 30 ESM(I,J) = 0.D0 C C *** Loop over Gauss pts to assemble ESM. C DO 50 N=1,8 XL = (1.0D0 + XII(N))/2.0D0 C C *** Locate the X & Y coordinates of the plate element sides at the C *** location of the Gauss point. Calculate the X & Y distance from C *** these side points to their corresponding nodes. C XIP = PLTC(2,MAT)*(1.0D0-XL) + PLTC(4,MAT)*XL YIP = PLTC(3,MAT)*(1.0D0-XL) + PLTC(5,MAT)*XL XJP = XIP + PLTP(7,NT)*CA YJP = YIP + PLTP(7,NT)*SA DXI = XIP - X(1) DYI = YIP - Y(1) DXJ = XJP - X(2) DYJ = YJP - Y(2) C C *** Assemble element stiffness matrix. C ZZ = FAC*WGT(N) ESM(1,1) = ESM(1,1) + ZZ*(AK*CC+SK*SS) ESM(2,1) = ESM(2,1) + ZZ*CS*(AK-SK) ESM(2,2) = ESM(2,2) + ZZ*(AK*SS+SK*CC) ESM(3,1) = ESM(3,1) + ZZ*((DXI*CS-DYI*CC)*AK-(DYI*SS+DXI*CS)*SK) ESM(3,2) = ESM(3,2) + ZZ*((DXI*SS-DYI*CS)*AK+(DYI*CS+DXI*CC)*SK) ESM(3,3) = ESM(3,3) + ZZ*((DYI*CA-DXI*SA)**2*AK + * (DYI*SA+DXI*CA)**2*SK) ESM(6,1) = ESM(6,1) + ZZ*((DYJ*CC-DXJ*CS)*AK+(DYJ*SS+DXJ*CS)*SK) ESM(6,2) = ESM(6,2) + ZZ*((DYJ*CS-DXJ*SS)*AK-(DYJ*CS+DXJ*CC)*SK) ESM(6,3) = ESM(6,3) + ZZ*((DYI*CA-DXI*SA)*(DXJ*SA-DYJ*CA)*AK + * (DYI*SA+DXI*CA)*(-DYJ*SA-DXJ*CA)*SK) 50 ESM(6,6) = ESM(6,6) + ZZ*((DXJ*SA-DYJ*CA)**2*AK + * (DYJ*SA+DXJ*CA)**2*SK) ESM(4,1) =-ESM(1,1) ESM(4,2) =-ESM(2,1) ESM(4,3) =-ESM(3,1) ESM(4,4) = ESM(1,1) ESM(5,1) =-ESM(2,1) ESM(5,2) =-ESM(2,2) ESM(5,3) =-ESM(3,2) ESM(5,4) = ESM(2,1) ESM(5,5) = ESM(2,2) ESM(6,4) =-ESM(6,1) ESM(6,5) =-ESM(6,2) DO 60 I=1,5 K = I+1 DO 60 J=K,6 60 ESM(I,J) = ESM(J,I) C C *** Assemble the stiffness matrix into the global stiffness matrix. C CALL ADSTIF(KFIX,S,LINK,ESM,6,NN) C 100 CONTINUE RETURN END SUBROUTINE PLTNLB(XX,LINK,MATNUM,NTYPE,DISP,W) IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /PLTPRP/ NPLTP,NPLTC,NP(20),PLTP(7,4),PLTC(5,20) COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA COMMON /GAUSS8/ WGT(8),XII(8) DIMENSION XX(NSD,1),LINK(NNPE,1),MATNUM(1),NTYPE(1), * DISP(NDOF,1),W(NDOF,1),DIS(6),QU(6),QV(6) C C *** Loop over elements. Check to see if plate element. C DO 100 NN=1,NELE IF(NTYPE(NN).NE.5) GOTO 100 MAT = MATNUM(NN) NT = NP(MAT) C C *** Set nodal coordinates and displacements. C N1 = LINK(1,NN) N2 = LINK(2,NN) X1 = XX(1,N1) X2 = XX(1,N2) Y1 = XX(2,N1) Y2 = XX(2,N2) DIS(1) = DISP(1,N1) DIS(2) = DISP(2,N1) DIS(3) = DISP(3,N1) DIS(4) = DISP(1,N2) DIS(5) = DISP(2,N2) DIS(6) = DISP(3,N2) C C *** Calculate and/or set plate properties. C CA = COS(PIR*PLTC(1,MAT)) SA = SIN(PIR*PLTC(1,MAT)) QU(1) =-CA QU(2) =-SA QU(4) = CA QU(5) = SA QV(1) =-SA QV(2) = CA QV(4) = SA QV(5) =-CA DL = SQRT((PLTC(4,MAT)-PLTC(2,MAT))**2 + * (PLTC(5,MAT)-PLTC(3,MAT))**2) FAC = PLTP(1,NT)*DL/2.0D0 DT = PLTP(5,NT)/PLTP(2,NT) DC = PLTP(4,NT)/PLTP(2,NT) DS = PLTP(6,NT)/PLTP(3,NT) C C *** Loop over Gauss points to assembly load-modifying vector. C DO 50 N=1,8 C C *** Locate the X & Y coordinates of the plate element sides at the C *** location of the Gauss point. Calculate elongation and horizontal C *** offset of a plate strand at location of the Gauss point. C XL = (1.0D0+XII(N))/2.0D0 XIP = PLTC(2,MAT)*(1.0D0-XL) + PLTC(4,MAT)*XL YIP = PLTC(3,MAT)*(1.0D0-XL) + PLTC(5,MAT)*XL XJP = XIP + PLTP(7,NT)*CA YJP = YIP + PLTP(7,NT)*SA QU(3) = (YIP-Y1)*CA - (XIP-X1)*SA QU(6) =-(YJP-Y2)*CA + (XJP-X2)*SA QV(3) = (YIP-Y1)*SA + (XIP-X1)*CA QV(6) =-(YJP-Y2)*SA - (XJP-X2)*CA DU = 0.0D0 DV = 0.0D0 DO 30 K=1,6 DU = DU + QU(K)*DIS(K) 30 DV = DV + QV(K)*DIS(K) C C *** Checks for plate yielding and add modifying loads to global C *** load-modifying vector. C ADU = ABS(DU) ADV = ABS(DV) ZZ = FAC*WGT(N) C IF(DU.GE.0.D0.AND.DU.LE.DT) THEN GOTO 40 ELSE IF(DU.GE.0.D0.AND.DU.GT.DT) THEN FAU = PLTP(2,NT)*DU - PLTP(5,NT) ELSE IF(ADU.LE.DC) THEN GOTO 40 ELSE FAU = PLTP(2,NT)*DU END IF DO 35 I=1,3 W(I,N1) = W(I,N1) + ZZ*QU(I)*FAU 35 W(I,N2) = W(I,N2) + ZZ*QU(I+3)*FAU C 40 IF(ADV.LE.DS) GOTO 50 IF(DV.GE.0.D0) THEN FAV = PLTP(3,NT)*DV - PLTP(6,NT) ELSE FAV = PLTP(3,NT)*DV + PLTP(6,NT) END IF DO 45 I=1,3 W(I,N1) = W(I,N1) + ZZ*QV(I)*FAV 45 W(I,N2) = W(I,N2) + ZZ*QV(I+3)*FAV C 50 CONTINUE 100 CONTINUE RETURN END SUBROUTINE FPCSTF(XX,KFIX,LINK,MATNUM,NTYPE,S) IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) DOUBLE PRECISION JAC(2,2) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA COMMON /GAUSS4/ SF(4,4,4),ETA(4,4,4),XXI(4,4,4),DLOC(4),WT(4) COMMON /FPCPRP/ NFPCP,NFPCC,NC(2,40),FPCC(9,40),FPCP(13,4) DIMENSION XX(NSD,1),KFIX(1),LINK(NNPE,1),MATNUM(1),NTYPE(1), * S(1),X(3),Y(3),ESM(6,6) C C *** Loop over elements. Check to see if element is a frame-plate C *** connector element. C DO 100 NN=1,NELE IF(NTYPE(NN).NE.4) GOTO 100 C C *** Set coordinates of the nodes. C MAT = MATNUM(NN) NT = NC(1,MAT) DO 20 L=1,2 J= LINK(L,NN) X(L) = XX(1,J) 20 Y(L) = XX(2,J) X(3) = XX(1,NC(2,MAT)) Y(3) = XX(2,NC(2,MAT)) C C *** Calculate connector stiffness with respect to displacements in C *** the global X and Y directions. C DX = X(3) - X(2) DY = Y(3) - Y(2) DL = SQRT(DX**2 + DY**2) SA = DY/DL CA = DX/DL CP = COS(PIR*FPCC(1,MAT)) SP = SIN(PIR*FPCC(1,MAT)) CE = CP*CA + SP*SA SE = CP*SA - SP*CA S1 = FPCP(2,NT)*FPCP(8,NT)/(FPCP(2,NT)*SE**2+FPCP(8,NT)*CE**2) S2 = FPCP(11,NT)*FPCP(5,NT)/(FPCP(11,NT)*SE**2+FPCP(5,NT)*CE**2) UK = CONFAC*S1*S2/(S1*SP**2+S2*CP**2) VK = CONFAC*S1*S2/(S2*SP**2+S1*CP**2) C C *** Zero ESM. Loop over gauss points to assemble ESM. C DO 30 I=1,6 DO 30 J=1,6 30 ESM(I,J) = 0.D0 DO 50 I=1,4 DO 50 J=1,4 C C *** Zero and then calculate the X and Y coordinates, coefficients of C *** the Jacobian matrix, and the determinant of the Jacobian at C *** Gauss point (I,J). C XCOR = 0.D0 YCOR = 0.D0 JAC(1,1) = 0.D0 JAC(2,1) = 0.D0 JAC(1,2) = 0.D0 JAC(2,2) = 0.D0 DO 40 K=1,4 M = K+1 N = K+5 XCOR = XCOR + FPCC(M,MAT)*SF(K,I,J) YCOR = YCOR + FPCC(N,MAT)*SF(K,I,J) JAC(1,1) = JAC(1,1) + FPCC(M,MAT)*ETA(K,I,J) JAC(2,1) = JAC(2,1) + FPCC(M,MAT)*XXI(K,I,J) JAC(1,2) = JAC(1,2) + FPCC(N,MAT)*ETA(K,I,J) 40 JAC(2,2) = JAC(2,2) + FPCC(N,MAT)*XXI(K,I,J) DET = JAC(1,1)*JAC(2,2) - JAC(2,1)*JAC(1,2) C C *** Assembly of the element stiffness matrix. C XI = XCOR - X(1) XJ = XCOR - X(2) YI = YCOR - Y(1) YJ = YCOR - Y(2) CON = DET*FPCP(1,NT)*WT(I)*WT(J) ESM(1,1) = ESM(1,1) + CON*UK ESM(2,2) = ESM(2,2) + CON*VK ESM(3,1) = ESM(3,1) - CON*UK*YI ESM(3,2) = ESM(3,2) + CON*VK*XI ESM(3,3) = ESM(3,3) + CON*(UK*YI*YI + VK*XI*XI) ESM(6,1) = ESM(6,1) + CON*UK*YJ ESM(6,2) = ESM(6,2) - CON*VK*XJ ESM(6,3) = ESM(6,3) - CON*(UK*YI*YJ + VK*XI*XJ) 50 ESM(6,6) = ESM(6,6) + CON*(UK*YJ*YJ + VK*XJ*XJ) ESM(4,1) =-ESM(1,1) ESM(4,3) =-ESM(3,1) ESM(4,4) = ESM(1,1) ESM(5,2) =-ESM(2,2) ESM(5,3) =-ESM(3,2) ESM(5,5) = ESM(2,2) ESM(6,4) =-ESM(6,1) ESM(6,5) =-ESM(6,2) DO 60 I=1,5 K = I+1 DO 60 J=K,6 60 ESM(I,J) = ESM(J,I) C C *** Assemble the stiffness matrix into the global stiffness matrix. C CALL ADSTIF(KFIX,S,LINK,ESM,6,NN) C 100 CONTINUE RETURN END SUBROUTINE FPCNLB(XX,LINK,MATNUM,NTYPE,DISP,W) IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) DOUBLE PRECISION JAC(2,2) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA COMMON /FPCPRP/ NFPCP,NFPCC,NC(2,40),FPCC(9,40),FPCP(13,4) COMMON /GAUSS4/ SF(4,4,4),ETA(4,4,4),XXI(4,4,4),DLOC(4),WT(4) DIMENSION XX(NSD,1),LINK(NNPE,1),MATNUM(1),NTYPE(1), * DISP(NDOF,1),W(NDOF,1),X(3),Y(3),U(2),V(2),R(2) C C *** Loop over elements. Check to see if element is a frame-plate C *** connector element. C DO 100 NN=1,NELE IF(NTYPE(NN).NE.4) GOTO 100 C C *** Set coordinates and displacements of the nodes. C MAT = MATNUM(NN) NT = NC(1,MAT) DO 20 L=1,2 J = LINK(L,NN) X(L) = XX(1,J) Y(L) = XX(2,J) U(L) = DISP(1,J) V(L) = DISP(2,J) 20 R(L) = DISP(3,J) X(3) = XX(1,NC(2,MAT)) Y(3) = XX(2,NC(2,MAT)) N1 = LINK(1,NN) N2 = LINK(2,NN) C C *** Calculate connector stiffness properties with respect to C *** displacements in the global X and Y directions. C DX = X(3) - X(2) DY = Y(3) - Y(2) DL = SQRT(DX**2 + DY**2) SA = DY/DL CA = DX/DL CP = COS(PIR*FPCC(1,MAT)) SP = SIN(PIR*FPCC(1,MAT)) CE = CP*CA + SP*SA SE = CP*SA - SP*CA S1 = FPCP(2,NT)*FPCP(8,NT)/(FPCP(2,NT)*SE**2+FPCP(8,NT)*CE**2) S2 = FPCP(11,NT)*FPCP(5,NT)/(FPCP(11,NT)*SE**2+FPCP(5,NT)*CE**2) UK = S1*S2/(S1*SP**2+S2*CP**2) VK = S1*S2/(S2*SP**2+S1*CP**2) S1 = FPCP(3,NT)*FPCP(9,NT)/(FPCP(3,NT)*SE**2+FPCP(9,NT)*CE**2) S2 = FPCP(12,NT)*FPCP(6,NT)/(FPCP(12,NT)*SE**2+FPCP(6,NT)*CE**2) UP0 = S1*S2/(S1*SP**2+S2*CP**2) VP0 = S1*S2/(S2*SP**2+S1*CP**2) IF(FPCP(4,NT).EQ.0.D0.OR.FPCP(10,NT).EQ.0.D0) THEN S1=0.D0 ELSE S1=FPCP(4,NT)*FPCP(10,NT)/(FPCP(4,NT)*SE**2+FPCP(10,NT)*CE**2) END IF IF(FPCP(13,NT).EQ.0.D0.OR.FPCP(7,NT).EQ.0.D0) THEN S2=0.D0 ELSE S2=FPCP(13,NT)*FPCP(7,NT)/(FPCP(13,NT)*SE**2+FPCP(7,NT)*CE**2) END IF IF(S1.EQ.0.D0.OR.S2.EQ.0.D0) THEN UP1=0.D0 VP1=0.D0 ELSE UP1=S1*S2/(S1*SP**2+S2*CP**2) VP1=S1*S2/(S2*SP**2+S1*CP**2) END IF C C *** Loop over Gauss pts. to assemble element load-modifying vector. C DO 60 I=1,4 DO 60 J=1,4 C C *** Zero and then calculate the X and Y coordinates, coefficients of C *** the Jacobian matrix, and the determinant of the Jacobian at C *** Gauss point (I,J). C XCOR = 0.D0 YCOR = 0.D0 JAC(1,1) = 0.D0 JAC(2,1) = 0.D0 JAC(1,2) = 0.D0 JAC(2,2) = 0.D0 DO 40 K=1,4 M = K+1 N = K+5 XCOR = XCOR + FPCC(M,MAT)*SF(K,I,J) YCOR = YCOR + FPCC(N,MAT)*SF(K,I,J) JAC(1,1) = JAC(1,1) + FPCC(M,MAT)*ETA(K,I,J) JAC(2,1) = JAC(2,1) + FPCC(M,MAT)*XXI(K,I,J) JAC(1,2) = JAC(1,2) + FPCC(N,MAT)*ETA(K,I,J) 40 JAC(2,2) = JAC(2,2) + FPCC(N,MAT)*XXI(K,I,J) DET = JAC(1,1)*JAC(2,2) - JAC(2,1)*JAC(1,2) C C *** Calculate slip at the Gauss point and the load-slip (stiffness) C *** properties in the direction of slip at the point. C DU = U(1) - R(1)*(YCOR-Y(1)) - U(2) + R(2)*(YCOR-Y(2)) DV = V(1) + R(1)*(XCOR-X(1)) - V(2) - R(2)*(XCOR-X(2)) DD = SQRT(DU**2+DV**2) CD = DU/DD SD = DV/DD DK = UK*VK/(UK*SD**2+VK*CD**2) DP0 = UP0*VP0/(UP0*SD**2+VP0*CD**2) IF(UP1.EQ.0.D0.OR.VP1.EQ.0.D0) THEN DP1 = 0.D0 ELSE DP1 = UP1*VP1/(UP1*SD**2+VP1*CD**2) END IF C C *** Calculate modifying loads. Assemble into global load-modifying C *** vector. C CON = DET*FPCP(1,NT)*WT(I)*WT(J) FAC = (DP0+DP1*DD)*(1.0D0-EXP(-DK*DD/DP0)) FX = CON*(UK*CONFAC*DU - FAC*DU/DD) FY = CON*(VK*CONFAC*DV - FAC*DV/DD) W(1,N1) = W(1,N1) + FX W(2,N1) = W(2,N1) + FY W(3,N1) = W(3,N1) - FX*(YCOR-Y(1)) + FY*(XCOR-X(1)) W(1,N2) = W(1,N2) - FX W(2,N2) = W(2,N2) - FY W(3,N2) = W(3,N2) + FX*(YCOR-Y(2)) - FY*(XCOR-X(2)) C 60 CONTINUE 100 CONTINUE RETURN END SUBROUTINE FRMSTF(XX,KFIX,LINK,MATNUM,NTYPE,MREL,S,F,DF) C C *** Subroutine locates each frame element and its corresponding C *** properties and nodal coordinates. The routine then calls C *** subroutine FESM for assembly of the stiffness matrix and sends C *** resulting element matrix to subroutine ADSTIF. Subroutine also C *** loads fixed end forces due to distributed load into load vector. C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /FRMPRP/ NFRMM,FRMP(4,8) COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT COMMON /FRAME / ESM(6,6),AREA,XI,EM,CS,SN,MR(2),FEF(6),DX,DY DIMENSION XX(NSD,1),KFIX(1),LINK(NNPE,1),MATNUM(1),NTYPE(1), * S(1),X(2),Y(2),ELSTF(6,6),F(NDOF,1),DF(2,1),N(2),MREL(2,1) C C *** Loop over elements -- generate [k] and assemble. C DO 100 NN=1,NELE C C *** Check to see if element is a frame type element. C IF(NTYPE(NN).NE.1) GOTO 100 C C *** Set frame element properties. C MAT = MATNUM(NN) AREA= FRMP(1,MAT) XI = FRMP(2,MAT) EM = FRMP(3,MAT) MR(1) = MREL(1,NN) MR(2) = MREL(2,NN) C C *** Set distributed loads on frame element. C NOLOAD = 0 DX = DF(1,NN) DY = DF(2,NN) IF(ABS(DX).LT.0.0001.AND.ABS(DY).LT.0.0001) NOLOAD=1 C C *** Set coordinates of the nodes. C DO 20 L=1,2 N(L) = LINK(L,NN) X(L) = XX(1,N(L)) 20 Y(L) = XX(2,N(L)) C C *** Assemble frame element stiffness matrix. C CALL FESM(X,Y,NOLOAD) C C *** Assemble the stiffness matrix into the structure. C DO 40 K=1,6 DO 40 L=1,6 40 ELSTF(K,L)=ESM(K,L) CALL ADSTIF(KFIX,S,LINK,ELSTF,6,NN) C C *** Assemble vector of fixed end forces into global load vector. C IF(NOLOAD.EQ.1) GOTO 100 DO 50 K=1,2 DO 50 I=1,3 50 F(I,N(K)) = F(I,N(K)) + FEF(K*3-3+I) C 100 CONTINUE C RETURN END SUBROUTINE FESM(XC,YC,NOLD) C C *** Subroutine to assemble the element stiffness matrix and fixed C *** end force vector for a frame element. C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /FRAME / ESM(6,6),AREA,XI,EM,CS,SN,MR(2),FEF(6),DX,DY DIMENSION XC(2),YC(2),DUM(6,6),DUMM(6) C C *** Calculate geometric and stiffness properties of the element. C V = YC(2)-YC(1) H = XC(2)-XC(1) EL = SQRT(H*H+V*V) CS = H/EL SN = V/EL EAOL = EM*AREA/EL EIOL = EM*XI/EL C C *** Calculation of the frame element stiffness matrix ESM. C ESM(1,1) = EAOL*CS*CS+12.*EIOL*SN*SN/EL/EL ESM(2,1) = (EAOL-12.*EIOL/EL/EL)*SN*CS ESM(2,2) = EAOL*SN*SN+12.*EIOL*CS*CS/EL/EL ESM(3,1) =-EIOL*6.*SN/EL ESM(3,2) = EIOL*6.*CS/EL ESM(3,3) = 4.*EIOL ESM(4,1) =-ESM(1,1) ESM(4,2) =-ESM(2,1) ESM(4,3) =-ESM(3,1) ESM(4,4) = ESM(1,1) ESM(5,1) =-ESM(2,1) ESM(5,2) =-ESM(2,2) ESM(5,3) =-ESM(3,2) ESM(5,4) = ESM(2,1) ESM(5,5) = ESM(2,2) ESM(6,1) = ESM(3,1) ESM(6,2) = ESM(3,2) ESM(6,3) = ESM(3,3)/2. ESM(6,4) =-ESM(3,1) ESM(6,5) =-ESM(3,2) ESM(6,6) = ESM(3,3) DO 40 I=1,5 K=I+1 DO 40 J=K,6 40 ESM(I,J)=ESM(J,I) C C *** Formation of fixed end force vector if there is a distributed C *** load on the element. C IF(NOLD.EQ.1) GOTO 45 FEF(1) = DX*ABS(V)/2.0 FEF(2) = DY*ABS(H)/2.0 FEF(3) = (-V*ABS(V)*DX+H*ABS(H)*DY)/12.0 FEF(4) = FEF(1) FEF(5) = FEF(2) FEF(6) =-FEF(3) C C *** Static condensation of fixed end force vector and element C *** stiffness matrix. (Release of moment for hinged frame elements). C 45 DO 70 K=1,2 IF(MR(K).NE.1) GOTO 70 KK = K*3 IF(NOLD.EQ.1) GOTO 55 DO 50 I=1,6 50 DUMM(I) = ESM(I,KK)*FEF(KK)/ESM(KK,KK) DO 52 I=1,6 52 FEF(I) = FEF(I) - DUMM(I) 55 DO 60 I=1,6 DO 60 J=1,6 60 DUM(I,J) = ESM(KK,J)*ESM(I,KK)/ESM(KK,KK) DO 65 I=1,6 DO 65 J=1,6 65 ESM(I,J) = ESM(I,J) - DUM(I,J) 70 CONTINUE C RETURN END SUBROUTINE FRMFOR(XX,LINK,MATNUM,NTYPE,MREL,LOUT,DISP,DF,IL,LSTEP) C C *** Subroutine calculates and outputs nodal forces and stresses C *** for the frame elements. C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /FRMPRP/ NFRMM,FRMP(4,8) COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT COMMON /FRAME / ESM(6,6),AREA,XI,EM,CS,SN,MR(2),FEF(6),DX,DY DIMENSION XX(NSD,1),LINK(NNPE,1),MATNUM(1),NTYPE(1),DISP(1), * X(2),Y(2),N(2),DI(6),EF(6),F(6),DF(2,1),MREL(2,1), * NFRM(100),S(100,6),LOUT(1) C WRITE(IOUT,2000) NFRME=0 C C *** Loop over elements. C DO 100 NN=1,NELE C C *** Check to see if element is a frame type element. C IF(NTYPE(NN).NE.1) GOTO 100 IF(LOUT(NN).NE.1) GOTO 100 C C *** Set frame element properties. C MAT = MATNUM(NN) AREA= FRMP(1,MAT) XI = FRMP(2,MAT) EM = FRMP(3,MAT) DEP = FRMP(4,MAT) MR(1) = MREL(1,NN) MR(2) = MREL(2,NN) C C *** Set distributed loads on frame element. C NOLOAD = 0 DX = DF(1,NN)*(IL+LSTEP-1) DY = DF(2,NN)*(IL+LSTEP-1) IF(ABS(DX).LT.0.00001.AND.ABS(DY).LT.0.00001) NOLOAD=1 C C *** Set coordinates and displacements of the nodes. C DO 20 L=1,2 N(L) = LINK(L,NN) X(L) = XX(1,N(L)) Y(L) = XX(2,N(L)) DO 20 LL=1,3 K = L*3-3+LL 20 DI(K) = DISP(N(L)*3-3+LL) C C *** Assemble frame element stiffness matrix. C CALL FESM(X,Y,NOLOAD) C C *** Calculation of element forces in the global coordinate system. C DO 40 I=1,6 IF(NOLOAD.EQ.1) FEF(I)=0.0D0 EF(I) = -FEF(I) DO 40 J=1,6 40 EF(I) = EF(I) + ESM(I,J)*DI(J) C C *** Transformation of the forces to the member coordinate system. C F(1)=CS*EF(1)+SN*EF(2) F(2)=-SN*EF(1)+CS*EF(2) F(3)=EF(3) F(4)=CS*EF(4)+SN*EF(5) F(5)=-SN*EF(4)+CS*EF(5) F(6)=EF(6) C C *** Output of frame element nodal forces. C WRITE(IOUT,2010)NN,N(1),F(1),F(2),F(3),N(2),F(4),F(5),F(6) C C *** Stress calculations. C NFRME = NFRME + 1 NFRM(NFRME) = NN S(NFRME,1) =-F(1)/AREA S(NFRME,4) =+F(4)/AREA S(NFRME,2) = ABS(1.5*F(2)/AREA) S(NFRME,5) = ABS(1.5*F(5)/AREA) S(NFRME,3) = ABS(F(3)*DEP/XI/2.0) S(NFRME,6) = ABS(F(6)*DEP/XI/2.0) C 100 CONTINUE C C *** Output of stresses. C WRITE(IOUT,2020) DO 110 I=1,NFRME NN = NFRM(I) N(1) = LINK(1,NN) N(2) = LINK(2,NN) APB1 = S(I,1)+S(I,3) AMB1 = S(I,1)-S(I,3) APB2 = S(I,4)+S(I,6) AMB2 = S(I,4)-S(I,6) 110 WRITE(IOUT,2030)NN,N(1),S(I,1),S(I,2),S(I,3),APB1,AMB1, * N(2),S(I,4),S(I,5),S(I,6),APB2,AMB2 C 2000 FORMAT(//,6X,49HFRAME ELEMENT NODAL FORCES (IN LOCAL COORDINATES), * /,6X,49(1H-),/,6X,7HElement,3X,4HNode,4X,11HAxial Force,4X, * 11HShear Force,3X,14HBending Moment,/,6X,61(1H-)) 2010 FORMAT(9X,I2,5X,I3,1X,3F15.5,/,16X,I3,1X,3F15.5) 2020 FORMAT(//,6X,31HSTRESSES AT FRAME ELEMENT NODES,/,6X,31(1H-),/, * 21X,5HAxial,6X,7HMaximum,5X,7HMaximum,7X,5HAxial,7X,5HAxial, * /,6X,45HElmt. Node Comp.(-) Shear Bending,9X,1H+, * 11X,1H-,/,19X,38HTension(+) (Absolute) (Absolute) , * 19HBending Bending,/,6X,70(1H-)) 2030 FORMAT(7X,I3,3X,I3,5(2X,F10.1),/,13X,I3,5(2X,F10.1)) C RETURN END SUBROUTINE FFCSTF(XX,KFIX,LINK,MATNUM,NTYPE,S) C C *** Subroutine to assembly the element stiffness matrix for the C *** frame-to-frame connector element. C C *** Description of some subroutine variables: C *** NP Configuration number associated with the element. C *** NC Connector property set associated with the element. C *** NCN Number of connectors in the configuration. C *** UK Initial stiffness of connector parallel to X axis. C *** VK Initial stiffness of connector parallel to Y axis. C *** D1(I) Vertical distance of connector I from element node 1. C *** D2(I) Vertical distance of connector I from element node 2. C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA COMMON /FFCPRP/ NFFCP,NFFCC,FFCP(6,3),NFFC(2,8),CFFC(3,8) COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT DIMENSION XX(NSD,1),KFIX(1),LINK(NNPE,1),MATNUM(1),NTYPE(1), * D1(3),D2(3),S(1),ESM(6,6) C C *** Loop over elements. Check to see if element is connector type. C DO 100 NN=1,NELE IF (NTYPE(NN).NE.2) GOTO 100 C C *** Set connector element properties. C NP = MATNUM(NN) NC = NFFC(1,NP) NCN = NFFC(2,NP) UK = FFCP(3,NC)*CONFAC VK = FFCP(6,NC)*CONFAC C C *** Set vertical distance between connectors and nodes 1 and 2. C N1 = LINK(1,NN) N2 = LINK(2,NN) YDIS = XX(2,N1)-XX(2,N2) DO 20 I=1,NCN D1(I) = CFFC(I,NP) 20 D2(I) = D1(I) + YDIS C C *** Zero and then assemble connector element stiffness matrix. C DO 30 I=1,6 DO 30 J=1,6 30 ESM(I,J)=0.D0 DO 40 I=1,NCN ESM(1,1) = ESM(1,1)+UK ESM(2,2) = ESM(2,2)+VK ESM(3,1) = ESM(3,1)-D1(I)*UK ESM(3,3) = ESM(3,3)+D1(I)**2*UK ESM(6,1) = ESM(6,1)+D2(I)*UK ESM(6,3) = ESM(6,3)-D1(I)*D2(I)*UK 40 ESM(6,6) = ESM(6,6)+D2(I)**2*UK ESM(4,1) =-ESM(1,1) ESM(4,3) =-ESM(3,1) ESM(4,4) = ESM(1,1) ESM(5,2) =-ESM(2,2) ESM(5,5) = ESM(2,2) ESM(6,4) =-ESM(6,1) DO 50 I=1,5 K=I+1 DO 50 J=K,6 50 ESM(I,J)=ESM(J,I) C C *** Assemble the stiffness matrix into the global stiffness matrix. C CALL ADSTIF(KFIX,S,LINK,ESM,6,NN) C 100 CONTINUE RETURN END SUBROUTINE FFCNLB(XX,LINK,MATNUM,NTYPE,DISP,W) C C *** Subroutine to assemble load-modifying vector so that nonlinear C *** behavior of frame-to-frame connectors is accounted for. C C *** Description of some subroutine variables. C *** NP Configuration number associated with the element. C *** NC Connector property set associated with the element. C *** NCN Number of connectors in the configuration. C *** N1,N2 Local node numbers 1 and 2. C *** D1,D2 Distance between the connector and nodes 1 and 2. C *** DV,DU Vertical and horizontal slip of connector. C *** FY,FX Vertical and horizontal shear force of connector. C *** Load-slip parameters - disp. parallel to the global X axis. C *** FFCP(1,NC)=P0 FFCP(2,NC)=P1 FFCP(3,NC)=K C *** Load-slip parameters - disp. parallel to the global Y axis. C *** FFCP(4,NC)=P0 FFCP(5,NC)=P1 FFCP(6,NC)=K C *** Load-slip parameters in the direction of connector slip. C *** AP0=P0 AP1=P1 AK=K C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA COMMON /FFCPRP/ NFFCP,NFFCC,FFCP(6,3),NFFC(2,8),CFFC(3,8) DIMENSION XX(NSD,1),LINK(NNPE,1),MATNUM(1),NTYPE(1), * DISP(NDOF,1),W(NDOF,1) C C *** Loop over elements. Check to see if element is connector type. C DO 100 NN=1,NELE IF(NTYPE(NN).NE.2) GOTO 100 C C *** Set connector properties and node numbers. C NP = MATNUM(NN) NC = NFFC(1,NP) NCN = NFFC(2,NP) N1 = LINK(1,NN) N2 = LINK(2,NN) YDIS = XX(2,N1)-XX(2,N2) DV = DISP(2,N1)-DISP(2,N2) C C *** Loop over each connector in element. C DO 30 I = 1,NCN C D1 = CFFC(I,NP) D2 = D1+YDIS DU = DISP(1,N1)-D1*DISP(3,N1)-DISP(1,N2)+D2*DISP(3,N2) DD = SQRT(DU**2+DV**2) IF(DD.EQ.0.D0) GOTO 30 CD = (DU/DD)**2 SD = (DV/DD)**2 C C *** Calculate load-slip parameters in the direction of connector disp. C AK = FFCP(3,NC)*FFCP(6,NC)/(FFCP(3,NC)*SD+FFCP(6,NC)*CD) AP0 = FFCP(1,NC)*FFCP(4,NC)/(FFCP(1,NC)*SD+FFCP(4,NC)*CD) IF(FFCP(2,NC).EQ.0.0D0.OR.FFCP(5,NC).EQ.0.D0) THEN AP1 = 0.0D0 ELSE AP1 = FFCP(2,NC)*FFCP(5,NC)/(FFCP(2,NC)*SD+FFCP(5,NC)*CD) END IF C C *** Calculate forces at the nodes and load into global load modifying C *** vector. C FF = (AP0+AP1*DD)*(1.0D0-EXP(-AK*DD/AP0)) FX = FFCP(3,NC)*DU*CONFAC - FF*DU/DD FY = FFCP(6,NC)*DV*CONFAC - FF*DV/DD W(1,N1) = W(1,N1) + FX W(2,N1) = W(2,N1) + FY W(3,N1) = W(3,N1) - FX*D1 W(1,N2) = W(1,N2) - FX W(2,N2) = W(2,N2) - FY W(3,N2) = W(3,N2) + FX*D2 30 CONTINUE C 100 CONTINUE RETURN END SUBROUTINE FFCFOR(XX,LINK,MATNUM,NTYPE,LOUT,DISP) C C *** Subroutine calculates and outputs shear forces the total slip C *** of each connector in a frame-to-frame connector element. C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /LIMITS/ ILOAD,NSTEP,NITER,NELAS,TOL,CONFAC,BETA COMMON /FFCPRP/ NFFCP,NFFCC,FFCP(6,3),NFFC(2,8),CFFC(3,8) COMMON /DEVICE/ IIN,IOUT,IBUG,NNOUT,NEOUT DIMENSION XX(NSD,1),LINK(NNPE,1),MATNUM(1),NTYPE(1), * DISP(NDOF,1),LOUT(1) C WRITE(IOUT,2000) C C *** Loop over elements. Check to see if element is connector type. C DO 100 NN=1,NELE IF(NTYPE(NN).NE.2) GOTO 100 IF(LOUT(NN).NE.1) GOTO 100 C C *** Set connector properties and node numbers. C NP = MATNUM(NN) NC = NFFC(1,NP) NCN = NFFC(2,NP) N1 = LINK(1,NN) N2 = LINK(2,NN) YDIS = XX(2,N1)-XX(2,N2) DV = DISP(2,N1)-DISP(2,N2) C C *** Loop over each connector in element. C DO 30 I = 1,NCN C D1 = CFFC(I,NP) D2 = D1+YDIS DU = DISP(1,N1)-D1*DISP(3,N1)-DISP(1,N2)+D2*DISP(3,N2) DD = SQRT(DU**2+DV**2) IF(DD.EQ.0.D0) THEN DA = 0.D0 FX = 0.D0 FY = 0.D0 FF = 0.D0 GOTO 20 END IF CD = (DU/DD)**2 SD = (DV/DD)**2 DA = DACOS(DU/DD)/PIR C C *** Calculate load-slip parameters in the direction of connector slip. C AK = FFCP(3,NC)*FFCP(6,NC)/(FFCP(3,NC)*SD+FFCP(6,NC)*CD) AP0 = FFCP(1,NC)*FFCP(4,NC)/(FFCP(1,NC)*SD+FFCP(4,NC)*CD) IF(FFCP(2,NC).EQ.0.0D0.OR.FFCP(5,NC).EQ.0.D0) THEN AP1 = 0.0D0 ELSE AP1 = FFCP(2,NC)*FFCP(5,NC)/(FFCP(2,NC)*SD+FFCP(5,NC)*CD) END IF C C *** Calculate and output forces at the connector. C FF = (AP0+AP1*DD)*(1.0D0-EXP(-AK*DD/AP0)) FX = ABS(FF*DU/DD) FY = ABS(FF*DV/DD) 20 IF(I.EQ.1) THEN WRITE(IOUT,2010)NN,DD,DA,FY,FX,FF ELSE WRITE(IOUT,2020)I,DD,DA,FY,FX,FF END IF C 30 CONTINUE 100 CONTINUE C 2000 FORMAT(//,6X,24HCONNECTOR ELEMENT FORCES,/,6X,24(1H-),/,6X, * 5HElmt.,2X,5HConn.,4X,5HTotal,4X,9HDirection,16X, * 11HShear Force,/,23X,4HSlip,5X,7Hof Slip,6X,31(1H-),/,31X, * '(Degrees) Vertical Horizontal Total',/,6X,70(1H-)) 2010 FORMAT(6X,I3,6X,1H1,4X,F8.6,3X,F9.4,3X,3F11.3) 2020 FORMAT(15X,I1,4X,F8.6,3X,F9.4,3X,3F11.3) RETURN END SUBROUTINE GAPNLB(XX,LINK,MATNUM,NTYPE,DISP,W,IBETA) C C *** Subroutine calculates the forces transferred between two frame C *** members due to the closure of a gap between the members. C *** Calculated forces are added to load-modifying vector. C IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) COMMON /KONTRL/ NUMNOD,NDOF,NELE,NNPE,NSD,NEQ,IBAND,PIR COMMON /GAPPRP/ NGAPS,GAPP(6,8) DIMENSION XX(NSD,1),LINK(NNPE,1),MATNUM(1),NTYPE(1), * DISP(NDOF,1),W(NDOF,1),N(2),U(2),V(2),R(2),X(2), * Y(2),XT(2),YT(2),XB(2),YB(2),XXT(2),XXB(2),DT(2),DB(2) C C *** Loop over elements. Check to see if element is a gap. C DO 100 NN=1,NELE IF(NTYPE(NN).NE.3) GOTO 100 C C *** Set element properties. C MAT = MATNUM(NN) STRESS = GAPP(1,MAT) WIDTH = GAPP(2,MAT) DT(1) = GAPP(3,MAT) DB(1) = GAPP(4,MAT) DT(2) = GAPP(5,MAT) DB(2) = GAPP(6,MAT) DO 20 I=1,2 N(I) = LINK(I,NN) X(I) = XX(1,N(I)) Y(I) = XX(2,N(I)) U(I) = DISP(1,N(I)) V(I) = DISP(2,N(I)) R(I) = DISP(3,N(I)) C C *** Calculate corner coordinates for displaced frame members. C XT(I) = X(I)+U(I)-R(I)*DT(I) XB(I) = X(I)+U(I)+R(I)*DB(I) YT(I) = Y(I)+DT(I)+V(I) 20 YB(I) = Y(I)-DB(I)+V(I) C C *** Set highest and lowest point of possible contact between members. C IF(YT(1).GT.YT(2)) THEN YTOP = YT(2) ELSE YTOP = YT(1) END IF IF(YB(1).LT.YB(2)) THEN YBOT = YB(2) ELSE YBOT = YB(1) END IF IF(YBOT.GE.YTOP) GOTO 100 C C *** Check to see if ends of displaced members are still parallel. If C *** they are and gap has closed, calculate crushing forces. C IF(ABS(R(1)-R(2)).GT.0.00001) GOTO 25 GAPN = X(2)+U(2)-X(1)-U(1) IF(GAPN.GE.0.D0) GOTO 100 22 P2 = STRESS*WIDTH*(YTOP-YBOT) YBAR = (YTOP+YBOT)/2. RM1 = P2*(YBAR-Y(1)-V(1)) RM2 =-P2*(YBAR-Y(2)-V(2)) GOTO 50 C C *** Calculate the intersection of the lines defined by the ends of C *** the displaced members (XC,YC). Calculate the slope of the C *** line which bisects the angle formed by the two lines (SINT). C 25 IF(ABS(R(1)).LE.0.00001) THEN XC = XT(1) YC = (XT(2)-XC)/R(2)+YT(2) SINT = -1.0D0/DTAN(DATAN(R(2))/2.0D0) ELSE IF(ABS(R(2)).LE.0.00001) THEN XC = XT(2) YC = (XT(1)-XC)/R(1)+YT(1) SINT = -1.0D0/DTAN(DATAN(R(1))/2.0D0) ELSE IF(ABS(R(1)+R(2)).LE.0.00001) THEN XXT(2) = (XT(1)+XT(2))/2 XXB(2) = XXT(2) YC =(XT(1)-XXT(2))/R(1)+YT(1) GOTO 45 ELSE YC = (R(2)*YT(2)+XT(2)-R(1)*YT(1)-XT(1))/(R(2)-R(1)) XC = R(1)*(YT(1)-YC)+XT(1) SINT= -1.0D0/DTAN((DATAN(R(1))+DATAN(R(2)))/2.0D0) END IF C C *** Determine amount of gap closure and associated crushing forces. C XXT(2) = (YT(2)-YC)/SINT+XC XXB(2) = (YB(2)-YC)/SINT+XC 45 IF(XXT(2).GT.XT(2).AND.XXB(2).GT.XB(2)) THEN GOTO 22 ELSE IF(XXT(2).LT.XT(2).AND.XXB(2).LT.XB(2)) THEN GOTO 100 ELSE IF(XXT(2).GT.XT(2).AND.XXB(2).LT.XB(2)) THEN IF(YC.GE.YTOP) GOTO 100 IF(YC.LE.YBOT) THEN YLIM = YBOT ELSE YLIM = YC END IF P2 = STRESS*WIDTH*(YTOP-YLIM) YBAR = (YTOP+YLIM)/2. RM1 = P2*(YBAR-Y(1)-V(1)) RM2 =-P2*(YBAR-Y(2)-V(2)) GOTO 50 ELSE IF(YC.LE.YBOT) GOTO 100 IF(YC.GE.YTOP) THEN YLIM = YTOP ELSE YLIM = YC END IF P2 = STRESS*WIDTH*(YLIM-YBOT) YBAR = (YLIM+YBOT)/2. RM1 = P2*(YBAR-Y(1)-V(1)) RM2 =-P2*(YBAR-Y(2)-V(2)) END IF C C *** Load forces into load-modifying vector. C 50 IBETA = 1 W(1,N(1)) = W(1,N(1))-P2 W(3,N(1)) = W(3,N(1))+RM1 W(1,N(2)) = W(1,N(2))+P2 W(3,N(2)) = W(3,N(2))+RM2 WRITE(*,2000)NN,P2 2000 FORMAT(15X,I3,F15.5) C 100 CONTINUE RETURN END