Maximum Program Size

JohnNichols · ‎09-10-2014

Steve:

I can now nicely run C# programs with the MKL libraries. Thanks for the help. The PARDISO routine is 100 times faster than the other inversion program I was using. It is a very nice addition to the FORTRAN tool kit.

Even though it is listed as a sparse matrix solver, it will actually solve dense matrices. It will do a 12500 by 12500 array in 167 seconds. But as soon as I go over the 12,500 I run into the 2GB limit for Windows. I could move to our supercomputer, but that is a pain, can I use allocate to break the 2GB limit?

JMN

andrew_4619 · ‎09-10-2014

You need to build a 64bit application for > 2GB

JohnNichols · ‎09-10-2014

https://software.intel.com/en-us/articles/memory-limits-applications-windows

Found this - which explains it - Steve wrote it

Thanks

JMN

John_Campbell · ‎09-11-2014

As well as the size limitations that Steve outlined for a 64-bit application, there are other limitations that need to be recognised.

Virtual memory limit: The maximum memory that can be allocated and used is limited by the size limits of the pagefile.sys. Although you can allocate a very large array, the limit of physical or virtual memory is not applied until the allocated array is being used. This can result in an error later in the program's execution.

Installed memory limit: The amount of installed memory also applies a practical limitation, as random access of virtual memory can significantly slow down the program run time. While this is not a limit on the addressable memory, it can be a definite practical limit. Use of programming techniques to localise memory usage can mitigate this problem. The aim is to sequentially address memory, which can be readily seen when using multi-dimension arrays. Simply re-ordering multiple DO loop sequences can result in a dramatic improvement of run time, much more so than the gains from other useful approaches like vector instructions or OpenMP approaches. The preferred solution is a combination of all these.

JohnNichols · ‎09-12-2014

Interesting, when I first started on computing we had strict size limits, and we copied clever ways to minimize the pain, for a brief while in the 1990's and 2000's it was no longer a major issue for most programmers, but now it comes back. How big is big enough?

I used to be able to design 4000 tonne bins with this program, I did not write it; it comes from Powell's group at UCB. I think I had the only electronic copy of the code for a while, it was on punch cards, then magnetic tape, then floppy drives and then CD's and now USB. it took me 30 years to find the author - who works at Halliburton in Houston.

1EDGE BEAM FORCES FOR LOADING CASE NUMBER 1 - the 1 on the format line reminds me of the old punch cards at UNI -

!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------    
    
      PROGRAM SHELL 

!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
        COMMON/ALL/NELEMS,R(100),Z(100),NODI(100),NODJ(100),NEQ,MBAND 
        COMMON/MISD/NUMJTS,NUMLDS,KODE(100) 
        DIMENSION TSTIFF(3000)
        CHARACTER*72 TITLE
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------   
100     CALL INPUT
! 
!     ESTABLISH BAND WIDTH
! 
      M=(MBAND-3)/3 
      DO 110 I =1,NELEMS
      J=IABS(NODJ(I)-NODI(I)) 
      IF(J.GT.M)M=J 
  110 CONTINUE
      MBAND=3*M+3 
      NEQ=3*NUMJTS
      NA=NEQ*MBAND
      IF(NA-3000)130,130,120
  120 WRITE(*,2000)
      WRITE(4,2000)
      GOTO 100 
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------        
130 CALL STIFF(TSTIFF,NEQ,MBAND)
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
      CALL DISPL(TSTIFF,NEQ,MBAND)
77    WRITE(*,78)
78    FORMAT(/2X,'DO YOU WISH TO CONTINUE (1=Yes/2=No) ? '\)
      READ(*,88)N
88    FORMAT(I1)
      IF(N.LT.1.OR.N.GT.2)GOTO 77
      CLOSE(3,STATUS='KEEP')
      CLOSE(4,STATUS='KEEP')
      IF(N.EQ.1)GOTO 100
      IF(N.EQ.1)GOTO 100
 
2000  FORMAT(1H1,'NUMBER OF NODAL POINTS AND/OR NODAL POINT','EXCEEDS MAXIMUM ALLOWABLE.')
 
135    END 

!------------------------------------------------------------------------------------------------
!
!         Subroutine Input for Data Files
!
!------------------------------------------------------------------------------------------------

     
        SUBROUTINE INPUT
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------   
        COMMON/ALL/NELEMS,R(100),Z(100),NODI(100),NODJ(100),NEQ,MBAND 
        COMMON/MISD/NUMJTS,NUMLDS,KODE(100) 
        COMMON/ISS/C11I(100),C12I(100),C22I(100),C33I(100),C34I(100), C44I(100),&
        C11J(100),C12J(100),C22J(100),C33J(100),C34J(100), C44J(100)
        COMMON/ISDS/STATIC(600),SPRING(300),TLOAD(300),S(10),NFORM
        COMMON/BEAM/LOCJ(10,6),LOCC(10),BSTRES(30),BMOD(10),BAREA(10),&
                 YY(10),ZZ(10),PHI(10),TDIM(10),NBEAMS 
        CHARACTER*64 FILEI,FILEO
        CHARACTER*72 TITLE
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
        WRITE(*,76)
76      FORMAT(//2X,'------------------------------------------------'/&
       ,/'          SHELL ANALYSIS PROGRAM',//,&
         '   Program date    : 7 July 2013'/'   Code Amendments : J.M. Nichols'/'   Version        : 1.1'//2X,&
                    '------------------------------------------------'//)
      WRITE(*,79)
79    FORMAT(2X,'Input Filename (e.g. B:file.ext) ? '\)
      READ(*,77)FILEI
77    FORMAT(A)
      WRITE(*,86)
86    FORMAT(1X,'Output Filename (e.g. B:file.ext) ? '\)
      READ(*,87)FILEO
87    FORMAT(A)
!     
!     OPEN EXISTING ( INPUT #3 ) AND NEW ( INPUT #4 ) FILES
!
      OPEN(3,FILE=FILEI,STATUS='OLD')
      OPEN(4,FILE=FILEO,STATUS='UNKNOWN')
!
      WRITE(4,111)FILEI,FILEO
111   FORMAT(//2X,'Input File : ',A,//1X,'Output File : ',A//)
!----------------------------------------------------------------------------------------
!
!       Wriie Out Header Data
!
!----------------------------------------------------------------------------------------  
        WRITE(4,76)
! 
!       CONTROL CARDS 
!
!----------------------------------------------------------------------------------------
!
!       Read the Title Cards
!
!----------------------------------------------------------------------------------------  
        READ(3,1000)TITLE
        WRITE(4,2001)TITLE
        Write(*,213)
213     Format('Reading Joint Cards')
        !Pause
        READ(3,*)NUMJTS,NELEMS,NUMLDS,NBEAMS,NFORM 
        WRITE(4,2000)NUMJTS,NELEMS,NUMLDS,NBEAMS,NFORM
        WRITE(*,2000)NUMJTS,NELEMS,NUMLDS,NBEAMS,NFORM
! 
!       JOINT CARDS 
! 
        WRITE(4,2010)
        WRITE(*,2010)
        NJ=1
        MJ=1
        DO 100 N=1,NUMJTS 
            IF(MJ-NJ.GE.1) GOTO 90 
            READ(3,*, END=4000)NJ, R(NJ),Z(NJ),KODE(NJ),S1,S2,S3
            IF(NJ-MJ.LE.1) GOTO 95 
            DIV=NJ-MJ 
            DR=(R(NJ)-R(MJ))/DIV
            DZ=(Z(NJ)-Z(MJ))/DIV
            NNJ=MJ
            MJ=NJ 
            NJ=NNJ
90          NJ=NJ+1 
            R(NJ)=R(NJ-1)+DR
            Z(NJ)=Z(NJ-1)+DZ
            KODE(NJ)=KODE(MJ) 
            SPRING(3*NJ-2)=R(NJ)*S1 
            SPRING(3*NJ-1)=R(NJ)*S2 
            SPRING(3*NJ)=R(NJ)*S3 
            GOTO 1011 
95          MJ=NJ 
            SPRING(3*NJ-2)=R(NJ)*S1 
            SPRING(3*NJ-1)=R(NJ)*S2 
            SPRING(3*NJ)=R(NJ)*S3 
1011        WRITE(*,2020)NJ, R(NJ),Z(NJ),KODE(NJ),S1,S2,S3
            WRITE(4,2020)NJ,R(NJ),Z(NJ),KODE(NJ),S1,S2,S3 
100     End Do      
! 
!       ELEMENT CARDS 
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
      
        Write(*,214)
214     Format('Reading Element Cards')
        WRITE(4,2030)
        Write(*,2030)
        !Pause
        M=0 
        NE=1
        ME=1
        DO 130 N=1,NELEMS 
            IF(ME-NE.GE.1) GO TO 101 
            Write(*,*)NE
            Read(3,*)NE
            Write(*,*)NE
            READ(3,*)NODI(N),NODJ(N),MTYPE,INC
            Write(*,*)NODI(N),NODJ(N),MTYPE,INC
            !Pause
            IF(INC.EQ.0) INC=1
            IF(NE-ME.LE.1)GOTO 105 
!
!           GENERATE INPUT DATA 
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
            TL=0.0
            NN=NE-ME-1
            DO 75 L=1,NN
                NODI(ME+L)=NODI(ME+L-1)+INC 
                NODJ(ME+L)=NODJ(ME+L-1)+INC 
                I=NODI(ME+L)
                J=NODJ(ME+L)
                DR=R(J)-R(I)
                DZ=Z(J)-Z(I)
                DL=SQRT(DR**2+DZ**2)
                TL=TL+DL
75          End Do                
                IF(MTYPE.GT.1) GOTO 80 
            READ(3,*)TMI,E1MI,E2MI,P12MI,TBI,E1BI,E2BI,P12BI 
            Write(3,*)TMI,E1MI,E2MI,P12MI,TBI,E1BI,E2BI,P12BI
            !Pause
            D1=(TMI-TMJ)/TL 
            D2=(E1MI-E1MJ)/TL 
            D3=(E2MI-E2MJ)/TL 
            D4=(P12MI-P12MJ)/TL 
            D5=(TBI-TBJ)/TL 
            D6=(E1BI-E1BJ)/TL 
            D7=(E2BI-E2BJ)/TL 
            D8=(P12BI-P12BJ)/TL
            GOTO 85
80          READ(3,1040)YMODI,AREA1I,AREA2I,BEND1I,BEND2I 
            D1=(YMODI-YMODJ)/TL 
            D2=(AREA1I-AREA1J)/TL 
            D3=(AREA2I-AREA2J)/TL 
            D4=(BEND1I-BEND1J)/TL 
            D5=(BEND2I-BEND2J)/TL 
85          NNE=ME
            ME=NE 
            NE=NNE
101         NE=NE+1 
            IF(NE.EQ.ME) GO TO 102
            I=NODI(NE)
            J=NODJ(NE)
            DR=R(J)-R(I)
            DZ=Z(J)-Z(I)
            DL=SQRT(DR**2+DZ**2)
! 
!           MATERIAL PROPERTIES 
! 
102         IF(MTYPE.GT.1)GOTO 103 
! 
!     MATERIAL ORTHOTROPY 
! 
      TMI=TMJ 
      TMJ=TMI+D1*DL 
      E1MI=E1MJ 
      E1MJ=E1MI+D2*DL 
      E2MI=E2MJ 
      E2MJ=E2MI+D3*DL 
      P12MI=P12MJ 
      P12MJ=P12MI+D4*DL 
      TBI=TBJ 
      TBJ=TBI+D5*DL 
      E1BI=E1BJ 
      E1BJ=E1BI+D6*DL 
      E2BI=E2BJ 
      E2BJ=E2BI+D7*DL 
      P12BI=P12BJ 
      P12BJ=P12BI+D8*DL 
      IF(NE.EQ.ME)READ(3,1030)TMJ,E1MJ,E2MJ,P12MJ,TBJ,E1BJ,E2BJ,P12BJ   
      WRITE(4,2050)NE,NODI(NE),TMI,E1MI,E2MI,P12MI,TBI,E1BI,E2BI,P12BI 
      WRITE(4,2060)NODJ(NE),TMJ,E1MJ,E2MJ,P12MJ,TBJ,E1BJ,E2BJ,P12BJ 
      GOTO 115 
! 
!     GEOMETRIC ORTHOTROPY
! 
  103 YMODI=YMODJ 
      YMODJ=YMODI+D1*DL 
      AREA1I=AREA1J 
      AREA1J=AREA1I+D2*DL 
      AREA2I=AREA2J 
      AREA2J=AREA2I+D3*DL 
      BEND1I=BEND1J 
      BEND1J=BEND1I+D4*DL 
      BEND2I=BEND2J 
      BEND2J=BEND2I+D5*DL 
      IF(NE.EQ.ME)READ(3,1040)YMODJ,AREA1J,AREA2J,BEND1J,BEND2J 
      WRITE(4,2080)NE,NODI(NE),YMODI,AREA1I,AREA2I,BEND1I,BEND2I
      WRITE(4,2090)NODJ(NE),YMODJ,AREA1J,AREA2J,BEND1J,BEND2J 
      GOTO 125 
! 
!     READ INPUT DATA 
! 
!     MATERIAL PROPERTIES 
! 
  105 ME=NE 
      IF(MTYPE.GT.1)GOTO 120 
! 
!     MATERIAL ORTHOTROPY 
! 
      IF(M.EQ.1)GOTO 110 
      WRITE(4,2040)
      WRITE(*,2040)
  110 READ(3,*)TMI,E1MI,E2MI,P12MI,TBI,E1BI,E2BI,P12BI 
      WRITE(4,2050)NE,NODI(NE),TMI,E1MI,E2MI,P12MI,TBI,E1BI,E2BI,P12BI 
      WRITE(*,2050)NE,NODI(NE),TMI,E1MI,E2MI,P12MI,TBI,E1BI,E2BI,P12BI 
      READ(3,*)TMJ,E1MJ,E2MJ,P12MJ,TBJ,E1BJ,E2BJ,P12BJ 
      WRITE(4,2060)NODJ(NE),TMJ,E1MJ,E2MJ,P12MJ,TBJ,E1BJ,E2BJ,P12BJ 
      WRITE(*,2060)NODJ(NE),TMJ,E1MJ,E2MJ,P12MJ,TBJ,E1BJ,E2BJ,P12BJ 
      M=1 
! 
!     FORM MATERIAL STIFFNESS COEFFICIENTS
! 
  115 C11I(NE)=TMI*E1MI**2/(E1MI-E2MI*P12MI**2) 
      C22I(NE)=C11I(NE)*E2MI/E1MI 
      C12I(NE)=P12MI*C22I(NE) 
      C33I(NE)=TBI**3*E1BI**2/(12.0*(E1BI-E2BI*P12BI**2)) 
      C44I(NE)=C33I(NE)*E2BI/E1BI 
      C34I(NE)=-P12BI*C44I(NE)
      C11J(NE)=TMJ*E1MJ**2/(E1MJ-E2MJ*P12MJ**2) 
      C22J(NE)=C11J(NE)*E2MJ/E1MJ 
      C12J(NE)=P12MJ*C22J(NE) 
      C33J(NE)=TBJ**3*E1BJ**2/(12.0*(E1BJ-E2BJ*P12BJ**2)) 
      C44J(NE)=C33J(NE)*E2BJ/E1BJ 
      C34J(NE)=-P12BJ*C44J(NE)
      GOTO 130 
! 
!     GEOMETRIC ORTHOTROPY
! 
  120 IF(M.EQ.2)GOTO 121 
      WRITE(4,2070)
  121 READ(3,*)YMODI,AREA1I,AREA2I,BEND1I,BEND2I 
      WRITE(4,2080)NE,NODI(NE),YMODI,AREA1I,AREA2I,BEND1I,BEND2I
      READ(3,*)YMODJ,AREA1J,AREA2J,BEND1J,BEND2J 
      WRITE(4,2090)NODJ(NE),YMODJ,AREA1J,AREA2J,BEND1J,BEND2J 
      M=2 
! 
!     FORM MATERIAL STIFFNESS COEFFICIENTS
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
  125 C11I(NE)=YMODI*AREA1I 
      C22I(NE)=YMODI*AREA2I 
      C12I(NE)=0.0
      C33I(NE)=YMODI*BEND1I 
      C44I(NE)=YMODI*BEND2I 
      C34I(NE)=0.0
      C11J(NE)=YMODJ*AREA1J 
      C22J(NE)=YMODJ*AREA2J 
      C12J(NE)=0.0
      C33J(NE)=YMODJ*BEND1J 
      C44J(NE)=YMODJ*BEND2J 
      C34J(NE)=0.0
  130 CONTINUE
! 
!     EDGE BEAM CARDS 
! 
      M=0 
      IF(NBEAMS.EQ.0) GOTO 160 
      WRITE(4,3000)
      Write(*,2140)
2140     format('Edge beams')
      !Pause
      WRITE(*,3000)
      DO 150 L=1,NBEAMS 
      READ(3,*)N,LOCC(N),(LOCJ(N,I),I=1,6) 
      READ(3,*)BMOD(N),BAREA(N),YY(N),ZZ(N),PHI(N),TDIM(N) 
      WRITE(4,3010)N,LOCC(N),LOCJ(N,1),BMOD(N),BAREA(N),YY(N),ZZ(N),PHI(N),TDIM(N)
      WRITE(*,3010)N,LOCC(N),LOCJ(N,1),BMOD(N),BAREA(N),YY(N),ZZ(N),PHI(N),TDIM(N)
      KK=0
      DO 140 K=1,6
      IF(LOCJ(N,K).LE.0) GOTO 145
      J=IABS(LOCC(N)-LOCJ(N,K)) 
      IF(J.GT.M) M=J
      KK=KK+1 
  140 CONTINUE
  145 WRITE(4,3020)(LOCJ(N,K),K=2,KK) 
  150 CONTINUE
  160 MBAND=3*M+3 
 
      RETURN
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------   
 1000 FORMAT( A72)
 1010 FORMAT(I3,2F10.1,I3,3F10.1) 
 1020 FORMAT(5I3) 
 1030 FORMAT(3F10.1,F7.4,3F10.1,F7.4) 
 1040 FORMAT(5F10.1)
 1050 FORMAT(8I3) 
 1060 FORMAT(6F10.1)
 2001 FORMAT( A72)
 2000 FORMAT(/////,&
      35H NUMBER OF NODAL POINTS ---------- ,I3// &
      35H NUMBER OF ELEMENTS -------------- ,I3// &
      35H NUMBER OF LOADING CONDITIONS ---- ,I3// &
      35H NUMBER OF EDGE BEAMS ------------ ,I3//&
      35H Stress Condition Output---------- ,I3)
 2010 FORMAT(11H JOINT DATA//5H NODE,5X,6HRADIUS,5X,8HORDINATE,5X,&
      4HCODE,6X,35H********** JOINT SPRINGS **********/44X, &
      35HVERTICAL    HORIZONTAL   ROTATIONAL/)
 2020 FORMAT(1H ,I3,5X,F7.2,5X,F7.2,6X,I3,6X,1PE10.3,3X,E10.3,3X,E10.3) 
 2030 FORMAT('ELEMENT DATA'/) 
 2040 FORMAT(14H0ELEMENT  NODE,5X,47H************* MEMBRANE PROPERTIES *&
      ************,5X,47H************** BENDING PROPERTIES *************&
      //19X,9HTHICKNESS,6X,14HYOUNGS MODULUS,4X,14HPOISSONS RATIO,5X, &
      9HTHICKNESS,6X,14HYOUNGS MODULUS,4X,14HPOISSONS RATIO/33X, &
      5H(ONE),7X,5H(TWO),5X,9H(ONE-TWO),21X,&
      5H(ONE),7X,5H(TWO),5X,9H(ONE-TWO)//)
 2050 FORMAT(2X,I3,5X,I3,7X,F7.4,3X,1PE10.3,2X,E10.3,5X,0PF6.3,10X,F7.4,&
      3X,1PE10.3,2X,E10.3,5X,0PF6.3)
 2060 FORMAT(10X,I3,7X,F7.4,3X,1PE10.3,2X,E10.3,5X,0PF6.3,10X,F7.4,3X,1PE10.3,2X,E10.3,5X,0PF6.3/) 
 2070 FORMAT(14H0ELEMENT  NODE,4X,14HYOUNGS MODULUS,6X,14HEFFECTIVE AREA&
      ,6X,27HEFFECTIVE MOMENT OF INERTIA/37X,5H(ONE),7X,5H(TWO),9X, &
      5H(ONE),7X,5H(TWO)//) 
 2080 FORMAT(2X,I3,5X,I3,7X,1PE10.3,4X,1PE10.3,2X,1PE10.3,4X,1PE10.3,2X,&
      1PE10.3)
 2090 FORMAT(10X,I3,7X,1PE10.3,4X,1PE10.3,2X,1PE10.3,4X,1PE10.3,2X,1PE10.3/) 
 3000 FORMAT(15H EDGE BEAM DATA///10H EDGE BEAM,4X,13HJOINT INDICES,4X,& 
      14HYOUNGS MODULUS,5X,4HAREA,10X,17HMOMENT OF INERTIA,5X,&
      14HPRINCIPAL AXIS,3X,17HLARGEST DIMENSION/&
      13X,15HCENTROID  NODES,36X,5H(IYY),7X,5H(IZZ)/) 
 3010 FORMAT(/4X,I3,8X,I3,5X,I3,7X,1PE10.3,4X,E10.3,4X,E10.3,2X,E10.3,&
      5X,E10.3,8X,E10.3)
3020 FORMAT(23X,I3)
4000 Write(*,4001)
4001 Format('Read Error')
! 
3021 END 
     
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
      SUBROUTINE STIFF(TSTIFF,NQQ,MBND) 
! 
      COMMON/ALL/NELEMS,R(100),Z(100),NODI(100),NODJ(100),NEQ,MBAND 
      COMMON/MISD/NUMJTS,NUMLDS,KODE(100) 
      COMMON/ISS/C11I(100),C12I(100),C22I(100),C33I(100),C34I(100),C44I(100),& 
                 C11J(100),C12J(100),C22J(100),C33J(100),C34J(100), &
                 C44J(100)
      COMMON/ISDS/STATIC(600),SPRING(300),TLOAD(300),S(10),NFORM
      COMMON/BEAM/LOCJ(10,6),LOCC(10),BSTRES(30),BMOD(10),BAREA(10),&
                  YY(10),ZZ(10),PHI(10),TDIM(10),NBEAMS 
      DIMENSION H(10),ESTIFF(7,7),F(7,7),ALAM(7,7),A(6,6),LM(6) 
      DIMENSION TSTIFF(NQQ,MBND)
! 
!     WEIGHTING COEFFICIENTS FOR NUMERICAL INTERGRATION 
! 
      S(1)=0.046910077030668
      S(2)=0.2307653449471585 
      S(3)=0.5
      S(4)=1.0-S(2) 
      S(5)=1.0-S(1) 
      H(1)=0.236926885056189
      H(2)=0.478628670499366
      H(3)=0.568888888888889
      H(4)=H(2) 
      H(5)=H(1) 
! 
!     CLEAR TOTAL STIFFNESS MATRIX
! 
      DO 100 I =1,NEQ 
      DO 100 J =1,MBAND 
  100 TSTIFF(I,J)=0.0 
! 
!     FORM ELEMENT STIFFNESS
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
      DO 200 L =1,NELEMS
      I=NODI(L) 
      J=NODJ(L) 
      DR=R(J)-R(I)
      DZ=Z(J)-Z(I)
      DL=SQRT(DR**2+DZ**2)
      COSA=DR/DL
      SINA=DZ/DL
      DO 110 M=1,7
      DO 110 N=M,7
  110 ESTIFF(M,N)=0.0 
      DO 140 K=1,5
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
! 
!     FORM MATERIAL STIFFNESS MATRIX
! 
      X=S(K)
      C11=C11I(L)+(C11J(L)-C11I(L))*X 
      C12=C12I(L)+(C12J(L)-C12I(L))*X 
      C22=C22I(L)+(C22J(L)-C22I(L))*X 
      C33=C33I(L)+(C33J(L)-C33I(L))*X 
      C34=C34I(L)+(C34J(L)-C34I(L))*X 
      C44=C44I(L)+(C44J(L)-C44I(L))*X 
! 
!     GENERALIZED DISPLACEMENT TRANSFORMATION MATRIX
! 
      IF(R(I).EQ.0.0)R(I)=0.000001
  120 RAD=R(I)+DL*COSA*X
      B1=(-3.0+4.0*X)/DL
      B2=(-1.0+4.0*X)/DL
      B3=(1.0-3.0*X+2.0*X*X)*COSA/RAD 
      B4=(1.0-(3.0-2.0*X)*X**2)*SINA/RAD
      B5=X*(1.0-2.0*X+X**2)*DL*SINA/RAD 
      B6=(-X+2.0*X*X)*COSA/RAD
      B7=X**2*(3.0-2.0*X)*SINA/RAD
      B8=X**2*(-1.0+X)*DL*SINA/RAD
      B9=(-6.0+12.0*X)/DL**2
      B10=(-4.0+6.0*X)/DL 
      B11=-B9 
      B12=(-2.0+6.0*X)/DL 
      B13=(6.0*X-6.0*X*X)*COSA/(DL*RAD) 
      B14=(-1.0+4.0*X-3.0*X*X)*COSA/RAD 
      B15=-B13
      B16=(2.0*X-3.0*X*X)*COSA/RAD
      B17=(4.0-8.0*X)/DL
      B18=(4.0*X-4.0*X*X)*COSA/RAD
! 
!     EVALUATE F(X) 
! 
      F(1,1)=B1*B1*C11+2.0*B1*B3*C12+B3*B3*C22
      F(1,2)=B1*B4*C12+B3*B4*C22
      F(1,3)=B1*B5*C12+B3*B5*C22
      F(1,4)=B1*B2*C11+B2*B3*C12+B1*B6*C12+B3*B6*C22
      F(1,5)=B1*B7*C12+B3*B7*C22
      F(1,6)=B1*B8*C12+B3*B8*C22
      F(1,7)=B1*B17*C11+B3*B17*C12+B1*B18*C12+B3*B18*C22
      F(2,2)=B4*B4*C22+B9*B9*C33+2.0*B9*B13*C34+B13*B13*C44 
      F(2,3)=B4*B5*C22+B9*B10*C33+(B10*B13+B9*B14)*C34+B13*B14*C44
      F(2,4)=B2*B4*C12+B4*B6*C22
      F(2,5)=B4*B7*C22+B9*B11*C33+(B11*B13+B9*B15)*C34+B13*B15*C44
      F(2,6)=B4*B8*C22+B9*B12*C33+(B12*B13+B9*B16)*C34+B13*B16*C44
      F(2,7)=B4*B17*C12+B4*B18*C22
      F(3,3)=B5*B5*C22+B10*B10*C33+2.0*B10*B14*C34+B14*B14*C44
      F(3,4)=B2*B5*C12+B5*B6*C22
      F(3,5)=B5*B7*C22+B10*B11*C33+(B11*B14+B10*B15)*C34+B14*B15*C44
      F(3,7)=B5*B17*C12+B5*B18*C22
      F(4,4)=B2*B2*C11+2.0*B2*B6*C12+B6*B6*C22
      F(4,5)=B2*B7*C12+B6*B7*C22
      F(4,6)=B2*B8*C12+B6*B8*C22
      F(4,7)=B2*B17*C11+(B6*B17+B2*B18)*C12+B6*B18*C22
      F(5,5)=B7*B7*C22+B11*B11*C33+2.0*B11*B15*C34+B15*B15*C44
      F(5,6)=B7*B8*C22+B11*B12*C33+(B12*B15+B11*B16)*C34+B15*B16*C44
      F(5,7)=B7*B17*C12+B7*B18*C22
      F(6,6)=B8*B8*C22+B12*B12*C33+2.0*B12*B16*C34+B16*B16*C44
      F(6,7)=B8*B17*C12+B8*B18*C22
      F(7,7)=B17*B17*C11+2.0*B17*B18*C12+B18*B18*C22
! 
!     ADD PARTIAL AREA TO ELEMENT STIFFNESS MATRIX
! 
      DO 130 M=1,7
      DO 130 N=M,7
  130 ESTIFF(M,N)=ESTIFF(M,N)+H(K)*F(M,N)*RAD*DL/2.0
  140 CONTINUE
! 
!     COMPLETE ELEMENT STIFFNESS
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
      DO 150 M=1,7
      DO 150 N=M,7
  150 ESTIFF(N,M)=ESTIFF(M,N) 
!
!     TRANSFORMATION MATRIX (GLOBAL COORDINATES)
! 
      DO 160 M=1,6
      DO 160 N=1,6
  160 ALAM(M,N)=0.0 
      ALAM(1,1)=SINA
      ALAM(1,2)=COSA
      ALAM(2,1)=-COSA 
      ALAM(2,2)=SINA
      ALAM(3,3)=-1.0
      ALAM(4,4)=SINA
      ALAM(4,5)=COSA
      ALAM(5,4)=-COSA 
      ALAM(5,5)=SINA
      ALAM(6,6)=-1.0
! 
!     STATIC CONDENSATION OF ELEMENT STIFFNESS MATRIX 
! 
      DO 151 M=1,6
  151 ESTIFF(7,M)=-(1.0/ESTIFF(7,7))*ESTIFF(7,M)
      DO 152 M=1,6
      N=6*L+M-6 
      STATIC(N)=0.0 
      DO 152 K=1,6
  152 STATIC(N)=STATIC(N)+ESTIFF(7,K)*ALAM(K,M) 
      DO 153 M=1,6
      DO 153 N=1,6
  153 A(M,N)=ESTIFF(M,7)*ESTIFF(7,N)
      DO 154 M=1,6
      DO 154 N=1,6
  154 ESTIFF(M,N)=ESTIFF(M,N)+A(M,N)
! 
!     MULTIPLY ELEMENT STIFFNESS AND TRANSFORMATION MATRICES
! 
      DO 170 M=1,6
      DO 170 N=1,6
      A(M,N)=0.0
      DO 170 K=1,6
  170 A(M,N)=A(M,N)+ESTIFF(M,K)*ALAM(K,N) 
! 
!     ELEMENT STIFFNESS IN GLOBAL COORDINATES 
! 
      DO 180 M=1,6
      DO 180 N=1,6
      ESTIFF(M,N)=0.0 
      DO 180 K=1,6
  180 ESTIFF(M,N)=ESTIFF(M,N)+ALAM(K,M)*A(K,N)
! 
!     ADD ELEMENT STIFFNESS TO TOTAL STIFFNESS
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
      LM(1)=3*I-2 
      LM(2)=LM(1)+1 
      LM(3)=LM(1)+2 
      LM(4)=3*J-2 
      LM(5)=LM(4)+1 
      LM(6)=LM(4)+2 
      DO 190 M=1,6
      II=LM(M)
      DO 190 N=M,6
      JJ=LM(N)-II+1 
  190 TSTIFF(II,JJ)=TSTIFF(II,JJ)+ESTIFF(M,N) 
  200 CONTINUE
! 
!     ESTABLISH LINKING BEAMS FOR EDGE BEAMS
! 
      IF(NBEAMS.EQ.0) GOTO 260 
      DO 250 L=1,NBEAMS 
      EI=BMOD(L)*TDIM(L)**3/12.0
      EA=BMOD(L)*TDIM(L)
      ANG=PHI(L)*0.01745329 
      DO 230 K=1,6
      IF(LOCJ(L,K).LE.0) GOTO 240
      IF(LOCJ(L,K).GT.LOCC(L)) GOTO 205
      I=LOCJ(L,K)
      J=LOCC(L) 
      GOTO 206 
  205 I=LOCC(L) 
      J=LOCJ(L,K) 
  206 DR=R(J)-R(I)
      DZ=Z(J)-Z(I)
      DL=SQRT(DR**2+DZ**2)
      COSA=DR/DL
      SINA=DZ/DL
      DO 210 M=1,6
      DO 210 N=M,6
  210 ESTIFF(M,N)=0.0 
      ESTIFF(1,1)=12.0*EI*COSA**2/DL**3+EA*SINA**2/DL 
      ESTIFF(1,2)=-12.0*EI*SINA*COSA/DL**3+EA*SINA*COSA/DL
      ESTIFF(1,3)=6.0*EI*COSA/DL**2 
      ESTIFF(1,4)=-ESTIFF(1,1)
      ESTIFF(1,5)=-ESTIFF(1,2)
      ESTIFF(1,6)=ESTIFF(1,3) 
      ESTIFF(2,2)=12.0*EI*SINA**2/DL**3+EA*COSA**2/DL 
      ESTIFF(2,3)=-6.0*EI*SINA/DL**2
      ESTIFF(2,4)=ESTIFF(1,5) 
      ESTIFF(2,5)=-ESTIFF(2,2)
      ESTIFF(2,6)=ESTIFF(2,3) 
      ESTIFF(3,3)=4.0*EI/DL 
      ESTIFF(3,4)=-ESTIFF(1,3)
      ESTIFF(3,5)=-ESTIFF(2,3)
      ESTIFF(3,6)=2.0*EI/DL 
      ESTIFF(4,4)=ESTIFF(1,1) 
      ESTIFF(4,5)=ESTIFF(1,2) 
      ESTIFF(4,6)=-ESTIFF(1,3)
      ESTIFF(5,5)=ESTIFF(2,2) 
      ESTIFF(5,6)=-ESTIFF(2,3)
      ESTIFF(6,6)=ESTIFF(3,3) 
! 
!     ADD LINK BEAM STIFFNESS TO TOTAL STIFFNESS
! 
      LM(1)=3*I-2 
      LM(2)=LM(1)+1 
      LM(3)=LM(1)+2 
      LM(4)=3*J-2 
      LM(5)=LM(4)+1 
      LM(6)=LM(4)+2 
      DO 220 M=1,6
      II=LM(M)
      DO 220 N=M,6
      JJ=LM(N)-II+1 
  220 TSTIFF(II,JJ)=TSTIFF(II,JJ)+ESTIFF(M,N)*(R(I)+R(J))/2.0 
  230 CONTINUE
! 
!     ADD EDGE BEAM STIFFNESS TO TOTAL STIFFNESS
! 
  240 COSA=COS(ANG) 
      SINA=SIN(ANG) 
      I=LOCC(L) 
      II=3*I-1
      TSTIFF(II,1)=TSTIFF(II,1)+BMOD(L)*BAREA(L)/R(I) 
      TSTIFF(II+1,1)=TSTIFF(II+1,1)+(ZZ(L)*SINA**2+YY(L)*COSA**2)*BMOD(L)/R(I)
! 
!     STORE EDGE BEAM STRESS MATRIX 
! 
      BSTRES(3*L-2)=BMOD(L)*BAREA(L)/R(I) 
      BSTRES(3*L-1)=SINA*BMOD(L)*ZZ(L)/R(I) 
  250 BSTRES(3*L)=-COSA*BMOD(L)*YY(L)/R(I)
  260 CONTINUE
! 
!     ADD SPRING STIFFNESS TO TOTAL STIFFNESS 
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
      DO 270 L=1,NUMJTS 
      L4 = (3*L)-2    
      Write(*,*)L4,NUMJTS
      TSTIFF(3*L-2,1)=TSTIFF(3*L-2,1)+SPRING(3*L-2) 
      TSTIFF(3*L-1,1)=TSTIFF(3*L-1,1)+SPRING(3*L-1) 
  270 TSTIFF(3*L,1)=TSTIFF(3*L,1)+SPRING(3*L) 
! 
      RETURN
! 
    END 
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
      SUBROUTINE DISPL(TSTIFF,NQQ,MBND) 
! 
      COMMON/ALL/NELEMS,R(100),Z(100),NODI(100),NODJ(100),NEQ,MBAND 
      COMMON/MISD/NUMJTS,NUMLDS,KODE(100) 
      COMMON/ISDS/STATIC(600),SPRING(300),TLOAD(300),S(10),NFORM
      DIMENSION TSTIFF(NQQ,MBND)
! 
!     MODIFY TOTAL STIFFNESS MATRIX FOR SPECIFIED JOINT DISPLACEMENTS 
! 
      DO 50 N=1,NUMJTS
      KK=KODE(N)
      KD=100
      DO 40 M=1,3 
      IF(KK.LT.KD) GOTO 30 
      NR=3*N-3+M
      DO 10 J=2,MBAND 
      TSTIFF(NR,J)=0.0
      NNR=NR-J+1
      IF(NNR.LT.1) GOTO 10 
      TSTIFF(NNR,J)=0.0 
   10 CONTINUE
   20 TSTIFF(NR,1)=1.0
      KK=KK-KD
   30 KD=KD/10
   40 CONTINUE
   50 CONTINUE
      DO 60 I =1,NEQ
   60 TLOAD(I)=0.0
! 
!     REDUCE THE MODIFIED TOTAL STIFFNESS MATRIX
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
      KKK=1 
      CALL SYMSOL(TSTIFF,TLOAD,NEQ,MBAND,KKK) 
      DO 150 L=1,NUMLDS 
! 
!     READ IN NODAL LOADS FOR ONE LOADING CASE
! 
      DO 70 I=1,NEQ 
   70 TLOAD(I)=0.0
      READ(3,*)LOAD,NUMNLD,NUMELD
      WRITE(4,2000)LOAD,TITLE 
      IF(NUMNLD.LT.1) GOTO 90
      WRITE(4,2010)
      DO 80 I=1,NUMNLD
      READ(3,*)NJ,TLOAD(3*NJ-2),TLOAD(3*NJ-1),TLOAD(3*NJ)
   80 WRITE(4,2020)NJ,TLOAD(3*NJ-2),TLOAD(3*NJ-1),TLOAD(3*NJ) 
! 
!     READ ELEMENT LOADS FOR THIS CASE AND FORM EQUIVALENT NODAL LOADS
! 
   85 IF(NUMELD.LT.1) GOTO 110 
      WRITE(4,2000)LOAD,TITLE 
   90 IF(NUMELD.LT.1) GOTO 110 
      WRITE(4,2030)
      DO 100 K=1,NUMELD 
      READ(3,1020)NE,QNORM,QVERT,QHORIZ 
      WRITE(4,2040)NE,QNORM,QVERT,QHORIZ
      I=NODI(NE)
      J=NODJ(NE)
      DR=R(J)-R(I)
      DZ=Z(J)-Z(I)
      DL=SQRT(DR**2+DZ**2)
      COSA=DR/DL
      SINA=DZ/DL
      QU=QVERT*SINA+QHORIZ*COSA 
      QW=QNORM-QVERT*COSA+QHORIZ*SINA 
! 
!     EQUIVALENT LOADS IN LOCAL COORDINATES 
! 
      R1=(R(I)*DL/2.0+DL*DL*COSA/6.0)*QU
      R2=(R(I)*DL/2.0+3.0*DL*DL*COSA/20.0)*QW 
      R3=(R(I)*DL*DL/12.0+(DL**3)*COSA/30.0)*QW 
      R4=(R(I)*DL/2.0+DL*DL*COSA/3.0)*QU
      R5=(R(I)*DL/2.0+7.0*DL*DL*COSA/20.0)*QW 
      R6=(-R(I)*DL*DL/12.0-(DL**3)*COSA/20.0)*QW
! 
!     TRANSFORM EQUIVALENT LOADS TO GLOBAL COORDINATES
! 
      TLOAD(3*I-2)=TLOAD(3*I-2)+R1*SINA-R2*COSA 
      TLOAD(3*I-1)=TLOAD(3*I-1)+R1*COSA+R2*SINA 
      TLOAD(3*I)=TLOAD(3*I)-R3
      TLOAD(3*J-2)=TLOAD(3*J-2)+R4*SINA-R5*COSA 
      TLOAD(3*J-1)=TLOAD(3*J-1)+R4*COSA+R5*SINA 
  100 TLOAD(3*J)=TLOAD(3*J)-R6
! 
!     MODIFY LOAD VECTOR FOR DISPLACEMENT BOUNDARY CONDITIONS 
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------  
  110 DO 140 N=1,NUMJTS 
      KK=KODE(N)
      KD=100
      DO 130 M=1,3
      IF(KK.LT.KD) GOTO 120
      NR=3*N-3+M
      TLOAD(NR)=0.0 
      KK=KK-KD
  120 KD=KD/10
  130 CONTINUE
  140 CONTINUE
! 
!     SOLVE THIS LOAD CASE FOR JOINT DISPLACEMENTS
! 
      KKK=2 
      CALL SYMSOL(TSTIFF,TLOAD,NEQ,MBAND,KKK) 
! 
!     PRINT JOINT DISPLACEMENTS 
! 
      WRITE(4,2050)LOAD
      WRITE(4,2060)(NJ,TLOAD(3*NJ-2),TLOAD(3*NJ-1),TLOAD(3*NJ),NJ=1,NUMJTS) 
!
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------        
!     SOLVE FOR THE ELEMENT STRESSES
! 
      CALL STRESS(LOAD) 
! 
  150 CONTINUE
! 
      RETURN
! 
 1000 FORMAT(3I3) 
 1010 FORMAT(I3,3F10.2) 
 1020 FORMAT(I3,3F10.2) 
 2000 FORMAT(19HLOADING CASE NUMBER,I3,4H FOR,/1X,12A6) 
 2010 FORMAT(1H0,19HAPPLIED NODAL LOADS//5H NODE,5X,36H VERTICAL    HORIZONTAL   ROTATIONAL//)
 2020 FORMAT(1X,I3,6X,1PE10.3,3X,E10.3,3X,E10.3)
 2030 FORMAT(1H0,21HAPPLIED ELEMENT LOADS//8H ELEMENT,7X,'NORMAL','VERTICAL    HORIZONTAL'//) 
 2040 FORMAT(3X,I3,7X,1PE10.3,3X,E10.3,3X,E10.3)
 2050 FORMAT(1H1,43HJOINT DISPLACEMENTS FOR LOADING CASE NUMBER,I3,//15H NODE          ,5X,' VERTICAL    HORIZONTAL   ROTATIONAL'/)
 2060 FORMAT(1X,I3,6X,1PE10.3,3X,E10.3,3X,E10.3)
! 
      END 
      SUBROUTINE SYMSOL (A,B,NN,MM,KKK) 
! 
      DIMENSION A(NN,MM),B(1) 
! 
!     A = COEFFICIENT MATRIX. 
!     B = RIGHT HAND SIDE MATRIX. 
!     NN = NUMBER OF EQUATIONS (MAX 800 HERE).
!     MM = HALF BAND WIDTH (MAX 20 HERE). 
! 
      GOTO (1000,2000),KKK 
! 
!     REDUCE COEFFICIENT MATRIX.
! 
 1000 DO 280 N=1,NN 
      DO 260 L=2,MM 
      C=A(N,L)/A(N,1) 
      I = N+L-1 
      IF(NN-I) 260,240,240
  240 J=0 
      DO 250 K=L,MM 
      J=J+1 
  250 A(I,J)=A(I,J)-C*A(N,K)
  260 A(N,L)=C
  280 CONTINUE
      GOTO 500 
! 
!     REDUCE RIGHT HAND SIDE MATRIX.
! 
 2000 DO 290 N=1,NN 
      DO 285 L=2,MM 
      I=N+L-1 
      IF(NN-I) 290,285,285
  285 B(I)=B(I)-A(N,L)*B(N) 
  290 B(N)=B(N)/A(N,1)
! 
!     CARRY OUT BACK SUBSTITUTION.
! 
      N=NN
  300 N = N-1 
      IF(N) 350,500,350 
  350 DO 400 K=2,MM 
      L = N+K-1 
      IF(NN-L) 400,370,370
  370 B(N) = B(N) - A(N,K) * B(L) 
  400 CONTINUE
      GOTO 300 
 
  500 RETURN

    END 
    
!----------------------------------------------------------------------------------------
!
!
!
!----------------------------------------------------------------------------------------      
      SUBROUTINE STRESS(LOAD) 
! 
      COMMON/ALL/NELEMS,R(100),Z(100),NODI(100),NODJ(100),NEQ,MBAND 
      COMMON/ISS/C11I(100),C12I(100),C22I(100),C33I(100),C34I(100),C44I(100), &
                C11J(100),C12J(100),C22J(100),C33J(100),C34J(100),C44J(100)
      COMMON/ISDS/STATIC(600),SPRING(300),TLOAD(300),S(10),NFORM
      COMMON/BEAM/LOCJ(10,6),LOCC(10),BSTRES(30),BMOD(10),BAREA(10),&
                 YY(10),ZZ(10),PHI(10),TDIM(10),NBEAMS 
      DIMENSION EJDIS(7),SIGMA(6),STRES(4,7)
! 
      WRITE(4,2000)LOAD 
      DO 120 L=1,NELEMS 
! 
!     FORM MATRIX OF ELEMENT JOINT DISPLACEMENTS
! 
      WRITE(4,2010)L
      I=NODI(L) 
      J=NODJ(L) 
      EJDIS(1)=TLOAD(3*I-2) 
      EJDIS(2)=TLOAD(3*I-1) 
      EJDIS(3)=TLOAD(3*I) 
      EJDIS(4)=TLOAD(3*J-2) 
      EJDIS(5)=TLOAD(3*J-1) 
      EJDIS(6)=TLOAD(3*J) 
! 
!     FIND MID-SIDE DISPLACEMENT
! 
      EJDIS(7)=0.0
      DO 105 M=1,6
      N=6*L+M-6 
  105 EJDIS(7)=EJDIS(7)+STATIC(N)*EJDIS(M)
! 
!     CONSTANT TERMS FOR ELEMENT STRESS MATRIX
! 
      S(1)=0.0
      S(2)=0.5
      S(3)=1.0
      DR=R(J)-R(I)
      DZ=Z(J)-Z(I)
      DL=SQRT(DR**2+DZ**2)
      COSA=DR/DL
      SINA=DZ/DL
      M1=2
      M2=2
      IF(NFORM.GE.1)M1=1
      IF(NFORM.GE.1)M2=3
      DO 110 K=M1,M2
! 
!     FORM MATERIAL STIFFNESS MATRIX
! 
      X=S(K)
      C11=C11I(L)+(C11J(L)-C11I(L))*X 
      C12=C12I(L)+(C12J(L)-C12I(L))*X 
      C22=C22I(L)+(C22J(L)-C22I(L))*X 
      C33=C33I(L)+(C33J(L)-C33I(L))*X 
      C34=C34I(L)+(C34J(L)-C34I(L))*X 
      C44=C44I(L)+(C44J(L)-C44I(L))*X 
! 
!     GENERALIZED DISPLACEMENT TRANSFORMATION MATRIX
! 
      RAD=R(I)+DL*COSA*X
      B1=(-3.0+4.0*X)/DL
      B2=(-1.0+4.0*X)/DL
      B3=(1.0-3.0*X+2.0*X*X)*COSA/RAD 
      B4=(1.0-(3.0-2.0*X)*X**2)*SINA/RAD
      B5=X*(1.0-2.0*X+X**2)*DL*SINA/RAD 
      B6=(-X+2.0*X*X)*COSA/RAD
      B7=X**2*(3.0-2.0*X)*SINA/RAD
      B8=X**2*(-1.0+X)*DL*SINA/RAD
      B9=(-6.0+12.0*X)/DL**2
      B10=(-4.0+6.0*X)/DL 
      B11=-B9 
      B12=(-2.0+6.0*X)/DL 
      B13=(6.0*X-6.0*X*X)*COSA/(DL*RAD) 
      B14=(-1.0+4.0*X-3.0*X*X)*COSA/RAD 
      B15=-B13
      B16=(2.0*X-3.0*X*X)*COSA/RAD
      B17=(4.0-8.0*X)/DL
      B18=(4.0*X-4.0*X*X)*COSA/RAD
! 
!     EVALUATE STRESS MATRIX AT POINT S 
! 
      STRES(1,1)=SINA*B1*C11+(SINA*B3-COSA*B4)*C12
      STRES(1,2)=COSA*B1*C11+(COSA*B3+SINA*B4)*C12
      STRES(1,3)=-B5*C12
      STRES(1,4)=SINA*B2*C11+(SINA*B6-COSA*B7)*C12
      STRES(1,5)=COSA*B2*C11+(COSA*B6+SINA*B7)*C12
      STRES(1,6)=-B8*C12
      STRES(1,7)=B17*C11+B18*C12
      STRES(2,1)=SINA*B1*C12+(SINA*B3-COSA*B4)*C22
      STRES(2,2)=COSA*B1*C12+(COSA*B3+SINA*B4)*C22
      STRES(2,3)=-B5*C22
      STRES(2,4)=SINA*B2*C12+(SINA*B6-COSA*B7)*C22
      STRES(2,5)=COSA*B2*C12+(COSA*B6+SINA*B7)*C22
      STRES(2,6)=-B8*C22
      STRES(2,7)=B17*C12+B18*C22
      STRES(3,1)=-COSA*(B9*C33+B13*C34) 
      STRES(3,2)=SINA*(B9*C33+B13*C34)
      STRES(3,3)=-B10*C33-B14*C34 
      STRES(3,4)=-COSA*(B11*C33+B15*C34)
      STRES(3,5)=SINA*(B11*C33+B15*C34) 
      STRES(3,6)=-B12*C33-B16*C34 
      STRES(3,7)=0.0
      STRES(4,1)=-COSA*(B9*C34+B13*C44) 
      STRES(4,2)=SINA*(B9*C34+B13*C44)
      STRES(4,3)=-B10*C34-B14*C44 
      STRES(4,4)=-COSA*(B11*C34+B15*C44)
      STRES(4,5)=SINA*(B11*C34+B15*C44) 
      STRES(4,6)=-B12*C34-B16*C44 
      STRES(4,7)=0.0
! 
!     MULTIPLY ELEMENT STRESS MATRIX AND JOINT DISPLACEMENT MATRIX
! 
      DO 100 M=1,4
      SIGMA(M)=0.0
      DO 100 N=1,7
  100 SIGMA(M)=SIGMA(M)+STRES(M,N)*EJDIS(N) 
! 
!     PRINT ELEMENT STRESSES AT POINT S 
! 
      DX=X*DL 
      WRITE(4,2020)DX,(SIGMA(M),M=1,4)
! 
  110 CONTINUE
  120 CONTINUE
! 
!     EDGE BEAM FORCES
! 
      IF(NBEAMS.LT.1) GOTO 140 
      WRITE(4,2030)LOAD 
      DO 130 K=1,NBEAMS 
      L=LOCC(K) 
      FXX=BSTRES(3*K-2)*TLOAD(3*L-1)
      ZZM=BSTRES(3*K-1)*TLOAD(3*L)
      YYM=BSTRES(3*K)*TLOAD(3*L)
  130 WRITE(4,2040)K,FXX,YYM,ZZM
  140 CONTINUE
! 
      WRITE(*,141)
141   FORMAT(//2X,'ANALYSIS COMPLETED.'//)
      RETURN
! 
 2000 FORMAT( 'ELEMENT STRESS RESULTANTS FOR LOADING CASE NUMBER :',I4////15X,8HLOCATION,4X,27H****** NORMAL FORCES ******,3X,27H***** BENDING MOMENTS *****//27X,27HMERIDIONAL  CIRCUMFERENTIAL,3X,27HMERIDIONAL  CIRCUMFERENTIAL/)
 2010 FORMAT(8H ELEMENT,I3) 
 2020 FORMAT(15X,F6.2,6X,1PE10.3,4X,E10.3,6X,E10.3,4X,E10.3)
 2030 FORMAT(1H1,40HEDGE BEAM FORCES FOR LOADING CASE NUMBER,I3///10H EDGE BEAM,5X,11HAXIAL FORCE,8X,15HBENDING MOMENTS/18X,5H(FXX),10X,5H(MYY),7X,5H(MZZ)/)
 2040 FORMAT(3X,I3,9X,1PE10.3,5X,E10.3,3X,E10.3)
! 
      END


  Input File : sample.in                                                       

 Output File : s.out                                                           




  ------------------------------------------------

          SHELL ANALYSIS PROGRAM

   Program date    : 7 July 2013
   Code Amendments : J.M. Nichols
   Version        : 1.1

  ------------------------------------------------


This is a Sample input file.                                            





 NUMBER OF NODAL POINTS ----------  11

 NUMBER OF ELEMENTS --------------  10

 NUMBER OF LOADING CONDITIONS ----   1

 NUMBER OF EDGE BEAMS ------------   2

 Stress Condition Output----------   1
 JOINT DATA

 NODE     RADIUS     ORDINATE     CODE      ********** JOINT SPRINGS **********
                                            VERTICAL    HORIZONTAL   ROTATIONAL

   1        1.00        1.00        0       0.000E+00    0.000E+00    0.000E+00
   2        1.00        0.95        0       0.000E+00    0.000E+00    0.000E+00
   3        1.00        0.90        0       0.000E+00    0.000E+00    0.000E+00
   4        1.00        0.85        0       0.000E+00    0.000E+00    0.000E+00
   5        1.00        0.80        0       0.000E+00    0.000E+00    0.000E+00
   6        1.00        0.75        0       0.000E+00    0.000E+00    0.000E+00
   7        1.00        0.70        0       0.000E+00    0.000E+00    0.000E+00
   8        1.00        0.65        0       0.000E+00    0.000E+00    0.000E+00
   9        1.00        0.60        0       0.000E+00    0.000E+00    0.000E+00
  10        1.00        0.55        0       0.000E+00    0.000E+00    0.000E+00
  11        1.00        0.50      111       0.000E+00    0.000E+00    0.000E+00
ELEMENT DATA

0ELEMENT  NODE     ************* MEMBRANE PROPERTIES *************     ************** BENDING PROPERTIES *************

                   THICKNESS      YOUNGS MODULUS    POISSONS RATIO     THICKNESS      YOUNGS MODULUS    POISSONS RATIO
                                 (ONE)       (TWO)     (ONE-TWO)                     (ONE)       (TWO)     (ONE-TWO)


    1       1        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
            2        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

    2       2        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
            3        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

    3       3        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
            4        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

    4       4        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
            5        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

    5       5        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
            6        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

    6       6        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
            7        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

    7       7        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
            8        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

    8       8        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
            9        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

    9       9        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
           10        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

   10      10        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200
           11        0.2000    1.000E+01   1.000E+01      0.200           0.1000    1.100E+01   1.100E+01      0.200

 EDGE BEAM DATA


 EDGE BEAM    JOINT INDICES    YOUNGS MODULUS     AREA          MOMENT OF INERTIA     PRINCIPAL AXIS   LARGEST DIMENSION
             CENTROID  NODES                                    (IYY)       (IZZ)


      1          3       0        2.000E+01     2.000E-01     1.200E+01   1.200E+01      0.000E+00         3.000E+00


      2          5       0        2.000E+01     2.000E-01     1.200E+01   1.200E+01      0.000E+00         3.000E+00

LOADING CASE NUMBER  1 FOR
     pB
0APPLIED NODAL LOADS

 NODE      VERTICAL    HORIZONTAL   ROTATIONAL


   4       0.000E+00    2.000E+00    0.000E+00
1JOINT DISPLACEMENTS FOR LOADING CASE NUMBER  1

 NODE                VERTICAL    HORIZONTAL   ROTATIONAL

   1       4.585E-02   -4.966E-01    6.330E-02
   2       4.090E-02   -4.946E-01    2.831E-02
   3       3.595E-02   -4.947E-01   -1.200E-03
   4       3.055E-02   -6.204E-01   -4.107E+00
   5       2.343E-02   -7.686E-01   -2.099E-03
   6       1.537E-02   -8.404E-01    3.232E-01
   7       7.175E-03   -7.470E-01    6.639E+00
   8       1.286E-03   -3.977E-01    1.059E+01
   9      -7.543E-04   -4.065E-02    6.966E+00
  10      -3.590E-04    6.672E-02    6.097E-01
  11       0.000E+00    0.000E+00    0.000E+00
ELEMENT STRESS RESULTANTS FOR LOADING CASE NUMBER :   1



               LOCATION    ****** NORMAL FORCES ******   ***** BENDING MOMENTS *****

                           MERIDIONAL  CIRCUMFERENTIAL   MERIDIONAL  CIRCUMFERENTIAL

 ELEMENT  1
                 0.00      -7.248E-05    -9.933E-01       1.288E-03    -2.576E-04
                 0.03       2.987E-05    -9.908E-01       6.682E-04    -1.336E-04
                 0.05      -5.007E-05    -9.892E-01       4.840E-05    -9.683E-06
 ELEMENT  2
                 0.00      -8.583E-05    -9.892E-01       2.420E-03    -4.840E-04
                 0.03       2.458E-05    -9.890E-01       5.636E-04    -1.127E-04
                 0.05      -1.907E-05    -9.895E-01      -1.293E-03     2.586E-04
 ELEMENT  3
                 0.00      -8.086E-03    -9.911E-01       1.310E-01    -2.620E-02
                 0.02       3.563E-03    -1.063E+00       7.841E-02    -1.568E-02
                 0.05      -6.173E-03    -1.242E+00       2.582E-02    -5.164E-03
 ELEMENT  4
                 0.00       5.231E-03    -1.240E+00       2.576E-02    -5.151E-03
                 0.03      -3.565E-03    -1.441E+00      -7.839E-02     1.568E-02
                 0.05       9.019E-03    -1.535E+00      -1.825E-01     3.651E-02
 ELEMENT  5
                 0.00      -2.763E-03    -1.538E+00       1.767E-01    -3.534E-02
                 0.03      -2.841E-04    -1.613E+00      -6.212E-03     1.242E-03
                 0.05       3.889E-03    -1.680E+00      -1.891E-01     3.783E-02
 ELEMENT  6
                 0.00       7.600E-03    -1.679E+00       6.437E-02    -1.287E-02
                 0.03      -5.485E-03    -1.667E+00      -1.206E-01     2.412E-02
                 0.05       1.433E-02    -1.491E+00      -3.056E-01     6.112E-02
 ELEMENT  7
                 0.00       3.468E-03    -1.493E+00       1.113E-01    -2.226E-02
                 0.03      -3.436E-03    -1.195E+00      -7.554E-02     1.511E-02
                 0.05       1.026E-02    -7.934E-01      -2.624E-01     5.247E-02
 ELEMENT  8
                 0.00      -9.718E-03    -7.973E-01       2.571E-01    -5.142E-02
                 0.02       3.148E-03    -3.924E-01       6.931E-02    -1.386E-02
                 0.05      -2.888E-03    -8.187E-02      -1.185E-01     2.370E-02
 ELEMENT  9
                 0.00      -1.445E-02    -8.419E-02       3.093E-01    -6.187E-02
                 0.03       5.515E-03     1.066E-01       1.214E-01    -2.428E-02
                 0.05      -7.619E-03     1.319E-01      -6.658E-02     1.332E-02
 ELEMENT 10
                 0.00      -4.475E-03     1.325E-01       1.995E-01    -3.989E-02
                 0.03       5.270E-04     7.444E-02       1.164E-02    -2.329E-03
                 0.05       2.354E-03     4.709E-04      -1.762E-01     3.524E-02
1EDGE BEAM FORCES FOR LOADING CASE NUMBER  1


 EDGE BEAM     AXIAL FORCE        BENDING MOMENTS
                  (FXX)          (MYY)       (MZZ)

     1         -1.979E+00      2.881E-01    0.000E+00
     2         -3.074E+00      5.037E-01    0.000E+00

John_Campbell · ‎09-13-2014

John,

Thanks for the example. This is a very early version of code from Graham Powell as he published a banded solver by 1974. Is this an axisymmetric shell element ? I have not used this analysis. If it is axisymmetric, then it would be a challenge to exceed the memory requirement of even 32-bit.

Do you recall the real and integer precision and computer that was used when this code was developed ?

You will notice that this code was developed before virtual memory was a consideration, as it would now run faster if

      SUBROUTINE SYMSOL (A,B,NN,MM,KKK)  
!  
      DIMENSION A(NN,MM),B(1)    !   A(neq,mband)

was changed to

      SUBROUTINE SYMSOL (A,B,NN,MM,KKK)  
!  
      DIMENSION A(MM,NN),B(1)    !   A(mband,neq)

I have changed the program you have provided to:

Declare variables as either integer*4 or real*8, to cope with the expected precision,
include parameters for key dimensions
Introduced some fortran 90 syntax for the do loops in SYMSOL
Changed the subscript order in TSTIFF

I guessed the sample data file and produced similar results. Hopefully I have included minimal bugs.