Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
29376 Discussions

Error Passing Variable

rdellsy
Beginner
‎11-20-2009 09:39 AM
596 Views
Hey all. In the following code, there appears to be an issue passing the variable O from the main subroutine to the various failure subroutines. Code follows:
[plain]      SUBROUTINE ABQMAIN
C=======================================================================
C     OUTPUT DATA FROM SHOCK ANALYSIS ON PLATE
C=======================================================================
      INCLUDE 'aba_param.inc'
      INCLUDE 'shock_out.inc'
      DIMENSION ARRAY(513),JRRAY(NPRECD,513),LRUNIT(2,1)
      EQUIVALENCE (ARRAY(1),JRRAY(1,1))
C=======================================================================
C     File initialization
C=======================================================================
      WRITE(6,*) ' Enter the name of the results file (w/o .fil):'
      READ(5,'(A)', IOSTAT = J) FNAME
      WRITE(6,'(A)') FNAME
      WRITE(6,*) ' Composite plate (1) or NAB plate (2)'
      READ(5,*) MATL
      IF(MATL.EQ.1) THEN
          WRITE(6,*) 'COMPOSITE'
          WRITE(6,*) 'XT,XC,YT,YC,SL,ST?'
          READ(5,*) XT,XC,YT,YC,SL,ST
          PRINT *, XT,XC,YT,YC,SL,ST
      ELSE IF(MATL.EQ.2) THEN
          WRITE(6,*) 'NAB'
      END IF
      WRITE(6,*) ' Enter the necessary output coordinates: '
      READ(5,*) ROOT,QTR,CTR,TIP,BCNR,TCNR
      WRITE(6,*) ROOT,QTR,CTR,TIP,BCNR,TCNR
      print *, 'Read-in successful.'
      NELEM = 0
      NNODE = 0
      NRU=1
      LRUNIT(1,1)=8
      LRUNIT(2,1)=2
      LOUTF=0
      CALL INITPF(FNAME,NRU,LRUNIT,LOUTF)
      JUNIT=8
      CALL DBRNU(JUNIT)
      OHASH(:,:)=0
      OCHAN(:,:)=0
      OCRSO(:,:)=0
      OCRSN(:,:)=0
      OPUCK(:,:)=0
      OKROL(:,:)=0
      OLACR(:,:)=0
C=======================================================================
C     Initialize the output files
C=======================================================================
      OPEN(UNIT=15,FILE='stress.plt',STATUS='UNKNOWN')
      WRITE(15,*) 'VARIABLES="X","Y","t","S11 (MPa)",
     *"S22 (MPa)","S12 (MPa)"'
      IF(MATL.EQ.1) THEN
          OPEN(UNIT=101,FILE='Hashin.plt',STATUS='UNKNOWN')
          WRITE(101,*) 'VARIABLES="X","Y","t","FIBER","MATRIX"'
          OPEN(UNIT=102,FILE='Chang.plt',STATUS='UNKNOWN')
          WRITE(102,*) 'VARIABLES="X","Y","t","FIBER","MATRIX"'
          OPEN(UNIT=103,FILE='ChristensenOld.plt',STATUS='UNKNOWN')
          WRITE(103,*) 'VARIABLES="X","Y","t","FIBER","MATRIX"'
          OPEN(UNIT=104,FILE='ChristensenNew.plt',STATUS='UNKNOWN')
          WRITE(104,*) 'VARIABLES="X","Y","t","FIBER","MATRIX"'
          OPEN(UNIT=105,FILE='Puck.plt',STATUS='UNKNOWN')
          WRITE(105,*) 'VARIABLES="X","Y","t","FIBER","MATRIX"'
          OPEN(UNIT=106,FILE='Kroll.plt',STATUS='UNKNOWN')
          WRITE(106,*) 'VARIABLES="X","Y","t","FIBER","MATRIX"'
          OPEN(UNIT=107,FILE='LaCR03.plt',STATUS='UNKNOWN')
          WRITE(107,*) 'VARIABLES="X","Y","t","FIBER","MATRIX"'
      ELSE IF(MATL.EQ.2) THEN
          OPEN(UNIT=101,FILE='PEEQ.plt',STATUS='UNKNOWN')
          WRITE(101,*) 'VARIABLES="X","Y","t","PEEQ"'
      END IF
      OPEN(UNIT=16,FILE='reactions.plt',STATUS='UNKNOWN')
      WRITE(16,*) 'VARIABLES="X","Y","t","RF1 (N)","RF2 (N)","RF3 (N)"'
      OPEN(UNIT=17,FILE='pressures.plt',STATUS='UNKNOWN')
      WRITE(17,*) 'VARIABLES="X","Y","Z","P(MPA)"'
      OPEN(UNIT=18,FILE='displacements.plt',STATUS='UNKNOWN')
      WRITE(18,*) 'VARIABLES="X","Y","Z"'
      OPEN(UNIT=112,FILE='root.plt',STATUS='UNKNOWN')
      WRITE(112,*) 'VARIABLES="t","P(MPa)",
     *"U3(m)","U(m)","V3(m/s)",
     *"V(m/s)","A3(m/s2)","A(m/s2)"'
      OPEN(UNIT=113,FILE='qtr_pt.plt',STATUS='UNKNOWN')
      WRITE(113,*) 'VARIABLES="t","P(MPa)",
     *"U3(m)","U(m)","V3(m/s)",
     *"V(m/s)","A3(m/s2)","A(m/s2)"'
      OPEN(UNIT=114,FILE='mid_pt.plt',STATUS='UNKNOWN')
      WRITE(114,*) 'VARIABLES="t","P(MPa)",
     *"U3(m)","U(m)","V3(m/s)",
     *"V(m/s)","A3(m/s2)","A(m/s2)"'
      OPEN(UNIT=115,FILE='tip.plt',STATUS='UNKNOWN')
      WRITE(115,*) 'VARIABLES="t","P(MPa)",
     *"U3(m)","U(m)","V3(m/s)",
     *"V(m/s)","A3(m/s2)","A(m/s2)"'
      OPEN(UNIT=116,FILE='twist.plt',STATUS='UNKNOWN')
      WRITE(116,*) 'VARIABLES="TIME","Phi"'
      print *, 'Outut initialization successful.'
C=======================================================================
C     Loop on all records in results file
C=======================================================================
      CALL DBFILE(0,ARRAY,JRCD)
      DO WHILE (JRCD.EQ.0)
        KEY = JRRAY(1,2)
C=======================================================================
C     READ RECORDS RELATING TO THE OVERALL MODEL - PRINTED ONE TIME
C=======================================================================
C     Read Element Definitions
C=======================================================================
        IF(KEY.EQ.1900) THEN
            IF(JRRAY(1,3).GT.2000000) THEN
                NELEM=NELEM+1
                IELEM(NELEM)=JRRAY(1,3)-2000000
                IREFE(IELEM(NELEM))=NELEM
            END IF                
C=======================================================================
C       Read Nodal Coordinates
C=======================================================================
        ELSE IF(KEY.EQ.1901) THEN
            NNODE=NNODE+1
C            print *, 'Node ',NNODE
            IF(JRRAY(1,3).GT.2000000) THEN
                INODE(NNODE)=JRRAY(1,3)-2000000
                IREFS(INODE(NNODE))=NNODE
            ELSE
                INODE(NNODE)=JRRAY(1,3)
                IREFF(INODE(NNODE))=NNODE
            END IF
            COORD(1,NNODE)=ARRAY(4)
            COORD(2,NNODE)=ARRAY(5)
            COORD(3,NNODE)=ARRAY(6)
C=======================================================================
C     READ RECORDS FOR EACH *STEP - PRINTED AT EACH TIME INCREMENT
C=======================================================================
C     Read the current step (NSTEP), increment number (NINC1), 
C     total time (TTIME), and step time (STIME)
C=======================================================================
        ELSE IF(KEY.EQ.2000) THEN
C        print *, 'Continuing...'
            NSTEP=JRRAY(1,8)
            NINC1=JRRAY(1,9)
            TTIME=ARRAY(3)*1000
            STIME=ARRAY(4)*1000
            PRINT *, NSTEP,NINC1,TTIME,STIME
            WRITE(15,*) 'ZONE T = "',TTIME,'",I=48,J=10,K=1,
     *DATAPACKING=POINT'
            IF(MATL.EQ.1) THEN
                WRITE(101,*) 'ZONE T = "',TTIME,'",I=48,J=10,K=1,
     *DATAPACKING=POINT'
            ELSE IF(MATL.EQ.2) THEN
                WRITE(101,*) 'ZONE T = "',TTIME,'",I=48,J=10,K=1,
     *DATAPACKING=POINT'
            END IF
            WRITE(16,*) 'ZONE T = "',TTIME,'",I=49,J=11,K=1,
     *DATAPACKING=POINT'
            WRITE(17,*) 'ZONE T = "',TTIME,'",I=49,J=49,K=2,
     *DATAPACKING=POINT'
            WRITE(18,*) 'ZONE T = "',TTIME,'",I=49,J=49,K=2,
     *DATAPACKING=POINT'
C=======================================================================
C     Read the element and integration point numbers, JELNUM AND INTPT,
C     and LOCATE, the location of INTPT (0--at int.pt., 1--at centroid,
C     2--at nodes, 4--nodal average) and the number of direct (NDI) and
C     shear (NSHR) components in the analysis
C=======================================================================
        ELSE IF(KEY.EQ.1) THEN
C        print *, 'Reading integration points.'
            JELNUM=JRRAY(1,2)-2000000
            INTPT=JRRAY(1,5)
            LOCATE=JRRAY(1,6)
            NDI=JRRAY(1,8)
            NSHR=JRRAY(1,9)
C=======================================================================
C     Read the direct and shear stresses and write to file
C=======================================================================
        ELSE IF(KEY.EQ.11) THEN
C        print *, 'Reading stresses. INTPT =',INTPT
            I11D=IREFS(JELNUM)
                DO IXX=1,NDI
                      STR(IXX) = ARRAY(IXX+2)
                END DO
                DO IXX=NDI+1,NDI+NSHR
                      STR(IXX) = ARRAY(IXX+2)
                END DO
                WRITE(15,'(7(1X,E14.6))')(COORD(J1+1,I11D),J1=1,2),TTIME
     *             ,(STR(J1)/1000000,J1=1,3)
C=======================================================================
C     Write the Hashin failure criteria
C=======================================================================
                IF(MATL.EQ.1) THEN
                print *, 'Checking failure criteria.'
C                    print *, 'Calling Hashin', OHASH(1,I11D), 
C     1                    OHASH(2,I11D)
C                    CALL HASHIN(STR,OHASH(:,I11D),HASH(:,I11D))
                    CALL CHANG(STR,OCHAN(:,I11D),CHAN(:,I11D))
C                    CALL CHRISTENOLD(STR,OCRSO(:,I11D),CRSO(:,I11D))
C                    CALL CHRISTENNEW(STR,OCRSN(:,I11D),CRSN(:,I11D))
C                    CALL PUCKSH(STR,OPUCK(:,I11D),PUCK(:,I11D))
C                    CALL KROLL(STR,OKROL(:,I11D),KROL(:,I11D))
C                    CALL LACR03(STR,OLACR(:,I11D),LACR(:,I11D))
C                    OHASH(:,I11D)=HASH(:,I11D)
                    OCHAN(:,I11D)=CHAN(:,I11D)
C                    OCRSO(:,I11D)=CRSO(:,I11D)
C                    OCRSN(:,I11D)=CRSN(:,I11D)
C                    OPUCK(:,I11D)=PUCK(:,I11D)
C                    OKROL(:,I11D)=KROL(:,I11D)
C                    OLACR(:,I11D)=LACR(:,I11D)
                    WRITE(101,'(6(1X,E14.3))')(COORD(J1+1,I11D),
     *              J1=1,2),TTIME,(HASH(J1,I11D),J1=1,2)
                    print *, 'Hashin',(COORD(J1+1,I11D), J1=1,2),TTIME,
     1              (HASH(J1,I11D),J1=1,2)
                END IF
C=======================================================================
C     Read the equivalent plastic strain values and write to file
C=======================================================================
        ELSE IF(KEY.EQ.73) THEN
            IF(MATL.EQ.2) THEN
                I73D=IREFS(JELNUM)
                PEEQ = ARRAY(3)
                WRITE(101,'(5(1X,E14.6))')(COORD(J1+1,I73D),
     *          J1=1,2),TTIME,PEEQ
            END IF
C=======================================================================
C     Read the nodal displacements and write to file
C=======================================================================
        ELSE IF(KEY.EQ.101) THEN
            IN=JRRAY(1,3)-1000000
            I101D=IREFS(IN)
            DO IXX=1,3
               DIS(IXX)=ARRAY(3+IXX)
            END DO
            WRITE(18,'(3(1X,E14.6))') (COORD(J1,I101D)+10*DIS(J1),
     *      J1=1,3)
            IF(JRRAY(1,3).EQ.ROOT) THEN
                U1R = DIS(1)
                U2R = DIS(2)
                U3R = DIS(3)
                UR = (DIS(1)**2+DIS(2)**2+DIS(3)**2)**0.5
            ELSE IF(JRRAY(1,3).EQ.QTR) THEN
                U1Q = DIS(1)
                U2Q = DIS(2)
                U3Q = DIS(3)
                UQ = (DIS(1)**2+DIS(2)**2+DIS(3)**2)**0.5
            ELSE IF(JRRAY(1,3).EQ.CTR) THEN
                U1C = DIS(1)
                U2C = DIS(2)
                U3C = DIS(3)
                UC = (DIS(1)**2+DIS(2)**2+DIS(3)**2)**0.5
            ELSE IF(JRRAY(1,3).EQ.TIP) THEN
                U1T = DIS(1)
                U2T = DIS(2)
                U3T = DIS(3)
                UT = (DIS(1)**2+DIS(2)**2+DIS(3)**2)**0.5
            ELSE IF(JRRAY(1,3).EQ.BCNR) THEN
                U3B = DIS(3)
                YBOT = COORD(2,I101D)
            ELSE IF(JRRAY(1,3).EQ.TCNR) THEN
                U3T = DIS(3)
                YTOP = COORD(2,I101D)
                YL = ABS(YTOP-YBOT)
                PHI = (ASIN((U3B-U3T)/YL))*(180/3.1415)
                WRITE(116,'(2(1X,E14.6))') TTIME,PHI
            END IF
C=======================================================================
C     Read the nodal velocities and write to file
C========================================================================
C        ELSE IF(KEY.EQ.102) THEN
C            IN=JRRAY(1,3)-1000000
C            I102D=IREFS(IN)
C            DO IXX=1,3
C               VEL(IXX)=ARRAY(3+IXX)
C            END DO
C            IF(JRRAY(1,3).EQ.ROOT) THEN
C                V1R = VEH(1)
C                V2R = VEH(2)
C                V3R = VEL(3)
C               VR = (VEH(1)**2+VEH(2)**2+VEL(3)**2)**0.5
C           ELSE IF(JRRAY(1,3).EQ.QTR) THEN
C                V1Q = VEH(1)
C                V2Q = VEH(2)
C                V3Q = VEL(3)
C                VQ = (VEH(1)**2+VEH(2)**2+VEL(3)**2)**0.5
C            ELSE IF(JRRAY(1,3).EQ.CTR) THEN
C                V1C = VEH(1)
C                V2C = VEH(2)
C                V3C = VEL(3)
C                VC = (VEH(1)**2+VEH(2)**2+VEL(3)**2)**0.5
C            ELSE IF(JRRAY(1,3).EQ.TIP) THEN
C                V1T = VEH(1)
C                V2T = VEH(2)
C                V3T = VEL(3)
C                VT = (VEH(1)**2+VEH(2)**2+VEL(3)**2)**0.5
C            END IF
C=======================================================================
C     Read the nodal accelerations and write to file
C=======================================================================
C        ELSE IF(KEY.EQ.103) THEN
C            IN=JRRAY(1,3)-1000000
C            I103D=IREFS(IN)
C            DO IXX=1,3
C               ACC(IXX)=ARRAY(3+IXX)
C            END DO
C            IF(JRRAY(1,3).EQ.ROOT) THEN
C                A1R = ACH(1)
C                A2R = ACH(2)
C                A3R = ACH(3)
C                AR = (ACH(1)**2+ACH(2)**2+ACH(3)**2)**0.5
C            ELSE IF(JRRAY(1,3).EQ.QTR) THEN
C                A1Q = ACH(1)
C                A2Q = ACH(2)
C                A3Q = ACH(3)
C                AQ = (ACH(1)**2+ACH(2)**2+ACH(3)**2)**0.5
C            ELSE IF(JRRAY(1,3).EQ.CTR) THEN
C                A1C = ACH(1)
C                A2C = ACH(2)
C                A3C = ACH(3)
C                AC = (ACH(1)**2+ACH(2)**2+ACH(3)**2)**0.5
c            ELSE IF(JRRAY(1,3).EQ.TIP) THEN
C                A1T = ACH(1)
C                A2T = ACH(2)
C                A3T = ACH(3)
C                AT = (ACH(1)**2+ACH(2)**2+ACH(3)**2)**0.5
C            END IF
C=======================================================================
C     Read the reaction forces at the root and write to file
C=======================================================================
        ELSE IF(KEY.EQ.104) THEN
            IN=JRRAY(1,3)-1000000
            I104D=IREFS(IN)
            DO IXX=1,6
               RXN(IXX)=ARRAY(3+IXX)
            END DO
            WRITE(16,'(7(1X,E14.6))')(COORD(J1+1,I104D),J1=1,2),TTIME,
     *             (RXN(J1),J1=1,3)
C=======================================================================
C     Read the nodal fluid pressures at the interface and write to file
C=======================================================================
        ELSE IF(KEY.EQ.108) THEN
            IN=JRRAY(1,3)
            I108D=IREFF(IN)
            PRS=ARRAY(4)/1000000
            WRITE(17,'(5(1X,E14.6))') (COORD(J1,I108D),J1=1,3),
     *             (PRS)
            IF((JRRAY(1,3)+1000000).EQ.ROOT) THEN
                PR = PRS
                WRITE(112,'(8(1X,E14.6))')TTIME,PR,U3R,UR,V3R,VR,A3R,AR
            ELSE IF(JRRAY(1,3)+1000000.EQ.QTR) THEN
                PQ = PRS
                WRITE(113,'(8(1X,E14.6))')TTIME,PQ,U3Q,UQ,V3Q,VQ,A3Q,AQ
            ELSE IF(JRRAY(1,3)+1000000.EQ.CTR) THEN
                PC = PRS
                WRITE(114,'(8(1X,E14.6))')TTIME,PC,U3C,UC,V3C,VC,A3C,AC
            ELSE IF(JRRAY(1,3)+1000000.EQ.TIP) THEN
                PT = PRS
                WRITE(115,'(8(1X,E14.6))')TTIME,PT,U3T,UT,V3T,VT,A3T,AT
            END IF
        END IF
C=======================================================================
C     Loop back to the next data set
C=======================================================================
        CALL DBFILE(0,ARRAY,JRCD)
C        print *,'Next time step.'
      END DO
C=======================================================================
C     Close the output files
C=======================================================================
      CLOSE(15)
      CLOSE(16)
      CLOSE(17)
      CLOSE(18)
      CLOSE(101)
      CLOSE(102)
      CLOSE(103)
      CLOSE(104)
      CLOSE(105)
      CLOSE(106)
C=======================================================================
C     Finish
C=======================================================================
      STOP
      END
C=======================================================================
C    ADDITIONAL SUBROUTINES
C=======================================================================
C
      SUBROUTINE HASHIN(S,O,H)
C
C=======================================================================
C     Subroutine for calculating the Hashin failure criteria
C
C     Fiber direction = 1 ; Transverse direction = 2
C
C     XT - Longitudinal tensile strength
C     XC - Longitudinal compressive strength
C     YT - Transverse tensile strength
C     YC - Transverse compressive strength
C     SL - Longitudinal shear strength
C     ST - Transverse shear strength
C
C     H(1) = Tensile(+1)/Compressive(-1) Fiber Failure
C     H(2) = Tensile(+1)/Compressive(-1) Matrix Failure
C=======================================================================
      INCLUDE 'aba_param.inc'
      INCLUDE 'shock_out.inc'
C..Initialize the Hashin vector
      print *, 'Checking Hashin Criteria'
C..Define the stresses
      S11 = S(1)
      S22 = S(2)
      S12 = S(3)
C..Fiber tensile failure
      IF(S11.GE.0) THEN
          print *, 'Fibre tensile failure'
          IF(O(1).EQ.1.0) THEN
             print *, '1'
             H(1)=1.0
             print *, '2'
          ELSE
             print *, '3'
             H(1)=((S11/XT)**2)+((S12/SL)**2)
             print *, '4'
             IF(H(1).GE.1.0) THEN
                print *, '5'
                H(1)=1.0
                print *, '6'
             END IF
             print *, '7'
          END IF
C..Fiber compressive failure
      ELSE
          print *, 'Fibre compressive failure'
          IF(O(1).EQ.-1.0) THEN
             H(1)=-1.0
          ELSE
             H(1)=(S11/XC)**2
             IF(H(1).LE.-1.0) THEN
                H(1)=-1.0
             END IF
          END IF
      END IF
C..Matrix tensile failure
      IF(S22.GE.0) THEN
          print *, 'Matrix tensile failure'
          IF(O(2).EQ.1.0) THEN
             H(2)=1.0
          ELSE
             H(2)=((S22/YT)**2)+((S12/SL)**2)
             IF(H(2).GE.1.0) THEN
                H(2)=1.0
             END IF
          END IF
C..Matrix compressive failure
      ELSE
          print *, 'Matrix compressive failure'
          IF(O(2).EQ.-1.0) THEN
             H(2)=-1.0
          ELSE
             H(2)=((S22/(2*ST))**2)+((((YC/(2*ST))**2)-1)*(S22/YC))+
     *            ((S12/SL)**2)
             IF(H(2).LE.-1.0) THEN
                H(2)=-1.0
             END IF
          END IF
      END IF
      print *, 'Hashin: Next criteria'
C      
      RETURN
      END

C=======================================================================
C
      SUBROUTINE CHANG(S,O,H)
C
C=======================================================================
C     Subroutine for calculating the Chang failure criteria
C
C     Fiber direction = 1 ; Transverse direction = 2
C
C     XT - Longitudinal tensile strength
C     XC - Longitudinal compressive strength
C     YT - Transverse tensile strength
C     YC - Transverse compressive strength
C     SL - Longitudinal shear strength
C     ST - Transverse shear strength
C
C     Note: H is used because C would lead to a comment.
C     H(1) = Tensile(+1)/Compressive(-1) Fiber Failure
C     H(2) = Tensile(+1)/Compressive(-1) Matrix Failure
C     H(3) = Delamination
C
C     A - Alpha: Linearity constant.
C=======================================================================
      INCLUDE 'aba_param.inc'
      INCLUDE 'shock_out.inc'
C..Initialize the Chang vector
C..Define the stresses
      S11 = S(1)
      S22 = S(2)
      S12 = S(3)
C..Fiber tensile failure
      IF(S11.GE.0) THEN
          IF(O(1).EQ.1.0) THEN
             H(1)=1.0
          ELSE
             H(1)=((S11/XT)**2)+(((S12**2/G12)+(.75*A*S12**4))/
     *            ((ST**2/G12)+(.75*A*ST)))
             IF(H(1).GE.1.0) THEN
                H(1)=1.0
             END IF
          END IF
C..Fiber compressive failure
      ELSE
          IF(O(1).EQ.-1.0) THEN
             H(1)=-1.0
          ELSE
             H(1)=((S11/XC)**2)+(((S12**2/G12)+(.75*A*S12**4))/
     *            ((ST**2/G12)+(.75*A*ST)))
             IF(H(1).LE.-1.0) THEN
                H(1)=-1.0
             END IF
          END IF
      END IF
C..Matrix tensile failure
      IF(S22.GE.0) THEN
          IF(O(2).EQ.1.0) THEN
             S(2)=1.0
          ELSE
             H(2)=((S22/YT)**2)+(((S12**2/G12)+(.75*A*S12**4))/
     *            ((ST**2/G12)+(.75*A*ST)))
             IF(H(2).GE.1.0) THEN
                H(2)=1.0
             END IF
          END IF
C..Matrix compressive failure
      ELSE
          IF(O(2).EQ.-1.0) THEN
             H(2)=-1.0
          ELSE
             H(2)=((S22/YT)**2)+(((S12**2/G12)+(.75*A*S12**4))/
     *            ((ST**2/G12)+(.75*A*ST)))
             IF(H(2).LE.-1.0) THEN
                H(2)=-1.0
             END IF
          END IF
      END IF
C      
      RETURN
      END

C=======================================================================
C
      SUBROUTINE CHRISTENOLD(S,O,H)
C
C=======================================================================
C     Subroutine for calculating the Christensen '87 failure criteria
C
C     Fiber direction = 1 ; Transverse direction = 2
C
C     XT - Longitudinal tensile strength
C     XC - Longitudinal compressive strength
C     YT - Transverse tensile strength
C     YC - Transverse compressive strength
C     SL - Longitudinal shear strength
C     ST - Transverse shear strength
C
C     Note: H is used because C would lead to a comment.
C     H(1) = Tensile(+1)/Compressive(-1) Fiber Failure
C     H(2) = Tensile(+1)/Compressive(-1) Matrix Failure
C     H(3) = Delamination
C=======================================================================
      INCLUDE 'aba_param.inc'
      INCLUDE 'shock_out.inc'
C..Initialize the Chang vector
C..Define the stresses
      S11 = S(1)
      S22 = S(2)
      S12 = S(3)
C..Fiber tensile failure
      IF(S11.GE.0) THEN
          IF(O(1).EQ.1.0) THEN
             H(1)=1.0
          ELSE
             H(1)=((S22/XT)**2)+((S12/SL)**2)
             IF(H(1).GE.1.0) THEN
                H(1)=1.0
             END IF
          END IF
C..Fiber compressive failure
      ELSE
          IF(O(1).EQ.-1.0) THEN
             H(1)=-1.0
          ELSE
             H(1)=((S22/XC)**2)+((S12/SL)**2)
             IF(H(1).LE.-1.0) THEN
                H(1)=-1.0
             END IF
          END IF
      END IF
C..Matrix tensile failure
      IF(S22.GE.0) THEN
          IF(O(2).EQ.1.0) THEN
             H(2)=1.0
          ELSE
             H(2)=((S22/YT)**2)+((S12/SL)**2)
             IF(H(2).GE.1.0) THEN
                H(2)=1.0
             END IF
          END IF
C..Matrix compressive failure
      ELSE
          IF(O(2).EQ.-1.0) THEN
             H(2)=-1.0
          ELSE
             H(2)=((S22/(2*ST))**2)+((((YC/(2*ST))**2)-1)*(S22/YC))+
     *            ((S12/SL)**2)
             IF(H(2).LE.-1.0) THEN
                H(2)=-1.0
             END IF
          END IF
      END IF
C      
      RETURN
      END

C=======================================================================
C
      SUBROUTINE CHRISTENNEW(S,O,H)
C
C=======================================================================
C     Subroutine for calculating the Chang failure criteria
C
C     Fiber direction = 1 ; Transverse direction = 2
C
C     XT - Longitudinal tensile strength
C     XC - Longitudinal compressive strength
C     YT - Transverse tensile strength
C     YC - Transverse compressive strength
C     SL - Longitudinal shear strength
C     ST - Transverse shear strength
C
C     Note: H is used because C would lead to a comment.
C     H(1) = Tensile(+1)/Compressive(-1) Fiber Failure
C     H(2) = Tensile(+1)/Compressive(-1) Matrix Failure
C     H(3) = Delamination
C=======================================================================
      INCLUDE 'aba_param.inc'
      INCLUDE 'shock_out.inc'
C..Initialize the Christensen New vector
C..Define the stresses
      S11 = S(1)
      S22 = S(2)
      S12 = S(3)
C..Fiber tensile failure
      IF(S11.GE.0) THEN
          IF(O(1).EQ.1.0) THEN
             H(1)=1.0
          ELSE
             H(1)=((S22/XT)**2)+((S12/SL)**2)
             IF(H(1).GE.1.0) THEN
                H(1)=1.0
             END IF
          END IF
C..Fiber compressive failure
      ELSE
          IF(O(1).EQ.-1.0) THEN
             H(1)=-1.0
          ELSE
             H(1)=((S22/XC)**2)+((S12/SL)**2)
             IF(H(1).LE.-1.0) THEN
                H(1)=-1.0
             END IF
          END IF
      END IF
C..Matrix tensile failure
      IF(S22.GE.0) THEN
          IF(O(2).EQ.1.0) THEN
             H(2)=1.0
          ELSE
             H(2)=((S22/YT)**2)+((S12/SL)**2)
             IF(H(2).GE.1.0) THEN
                H(2)=1.0
             END IF
          END IF
C..Matrix compressive failure
      ELSE
          IF(O(2).EQ.-1.0) THEN
             H(2)=-1.0
          ELSE
             H(2)=((S22/(2*ST))**2)+((((YC/(2*ST))**2)-1)*(S22/YC))+
     *            ((S12/SL)**2)
             IF(H(2).LE.-1.0) THEN
                H(2)=-1.0
             END IF
          END IF
      END IF
C      
      RETURN
      END

C=======================================================================
C
      SUBROUTINE PUCKSH(S,O,H)
C
C=======================================================================
C     Subroutine for calculating the Puck-Shurmann failure criteria
C
C     Fiber direction = 1 ; Transverse direction = 2
C
C     XT - Longitudinal tensile strength
C     XC - Longitudinal compressive strength
C     YT - Transverse tensile strength
C     YC - Transverse compressive strength
C     SL - Longitudinal shear strength
C     ST - Transverse shear strength
C
C     H(1) = Tensile(+1)/Compressive(-1) Fiber Failure
C     H(2) = Tensile(+1)/Compressive(-1) Matrix Failure
C     H(3) = Delamination
C=======================================================================
      INCLUDE 'aba_param.inc'
      INCLUDE 'shock_out.inc'
C..Initialize the Puck vector
C..Define the stresses
      S11 = S(1)
      S22 = S(2)
      S12 = S(3)
C..Fiber tensile failure
      IF(S11.GE.0) THEN
          IF(O(1).EQ.1.0) THEN
             H(1)=1.0
          ELSE
             H(1)=((S22/XT)**2)+((S12/SL)**2)
             IF(H(1).GE.1.0) THEN
                H(1)=1.0
             END IF
          END IF
C..Fiber compressive failure
      ELSE
          IF(O(1).EQ.-1.0) THEN
             H(1)=-1.0
          ELSE
             H(1)=((S22/XC)**2)+((S12/SL)**2)
             IF(H(1).LE.-1.0) THEN
                H(1)=-1.0
             END IF
          END IF
      END IF
C..Matrix tensile failure
      IF(S22.GE.0) THEN
          IF(O(2).EQ.1.0) THEN
             H(2)=1.0
          ELSE
             H(2)=((S22/YT)**2)+((S12/SL)**2)
             IF(H(2).GE.1.0) THEN
                H(2)=1.0
             END IF
          END IF
C..Matrix compressive failure
      ELSE
          IF(O(2).EQ.-1.0) THEN
             H(2)=-1.0
          ELSE
             H(2)=((S22/(2*ST))**2)+((((YC/(2*ST))**2)-1)*(S22/YC))+
     *            ((S12/SL)**2)
             IF(H(2).LE.-1.0) THEN
                H(2)=-1.0
             END IF
          END IF
      END IF
      
      RETURN
      END

C=======================================================================
C
      SUBROUTINE KROLL(S,O,H)
C
C=======================================================================
C     Subroutine for calculating the Kroll failure criteria
C
C     Fiber direction = 1 ; Transverse direction = 2
C
C     XT - Longitudinal tensile strength
C     XC - Longitudinal compressive strength
C     YT - Transverse tensile strength
C     YC - Transverse compressive strength
C     SL - Longitudinal shear strength
C     ST - Transverse shear strength
C
C     H(1) = Tensile(+1)/Compressive(-1) Fiber Failure
C     H(2) = Tensile(+1)/Compressive(-1) Matrix Failure
C     H(3) = Delamination
C=======================================================================
      INCLUDE 'aba_param.inc'
      INCLUDE 'shock_out.inc'
C..Initialize the Chang vector
C..Define the stresses
      S11 = S(1)
      S22 = S(2)
      S12 = S(3)
C..Fiber tensile failure
      IF(S11.GE.0) THEN
          IF(O(1).EQ.1.0) THEN
             H(1)=1.0
          ELSE
             H(1)=((S22/XT)**2)+((S12/SL)**2)
             IF(H(1).GE.1.0) THEN
                H(1)=1.0
             END IF
          END IF
C..Fiber compressive failure
      ELSE
          IF(O(1).EQ.-1.0) THEN
             H(1)=-1.0
          ELSE
             H(1)=((S22/XC)**2)+((S12/SL)**2)
             IF(H(1).LE.-1.0) THEN
                H(1)=-1.0
             END IF
          END IF
      END IF
C..Matrix tensile failure
      IF(S22.GE.0) THEN
          IF(O(2).EQ.1.0) THEN
             H(2)=1.0
          ELSE
             H(2)=((S22/YT)**2)+((S12/SL)**2)
             IF(H(2).GE.1.0) THEN
                H(2)=1.0
             END IF
          END IF
C..Matrix compressive failure
      ELSE
          IF(O(2).EQ.-1.0) THEN
             H(2)=-1.0
          ELSE
             H(2)=((S22/(2*ST))**2)+((((YC/(2*ST))**2)-1)*(S22/YC))+
     *            ((S12/SL)**2)
             IF(H(2).LE.-1.0) THEN
                H(2)=-1.0
             END IF
          END IF
      END IF
      
      RETURN
      END

C=======================================================================
C
      SUBROUTINE LACR03(S,O,H)
C
C=======================================================================
C     Subroutine for calculating the Chang failure criteria
C
C     Fiber direction = 1 ; Transverse direction = 2
C
C     XT - Longitudinal tensile strength
C     XC - Longitudinal compressive strength
C     YT - Transverse tensile strength
C     YC - Transverse compressive strength
C     SL - Longitudinal shear strength
C     ST - Transverse shear strength
C     
C     EPxx - Strain in the xx dirrection
C     EPxT/L - Transverse or longitudinal strain along the x direction.
C
C     M - Mixture constant (percent fibre)
C     A - Alpha - Fracture Plane Angle (Assume 52 deg, bu 40, 45, 0 ok)
C     A0 - Max Fracture Plane Angle (For Carbon Fibre: 53 deg)
C
C     H(1) = Tensile(+1)/Compressive(-2)/Kinking(-1) Fiber Failure
C     H(2) = Tensile(+1)/Compressive(-1) Matrix Failure
C     LK = Kinking failure test.
C=======================================================================
      INCLUDE 'aba_param.inc'
      INCLUDE 'shock_out.inc'
C..Initialize the LaCR03 vector
C..Define the stresses
      S11 = S(1)
      S22 = S(2)
      S12 = S(3)
C..Specify Constants:
      A=3.1415926*52/180
      A0=3.1316926*53/180
      M=0.6
C..Fiber tensile failure
      IF(S11.GE.0) THEN
          IF(O(1).EQ.1.0) THEN
             H(1)=1.0
          ELSE
C             H(1)=(EP11/EP1T)
             IF(H(1).GE.1.0) THEN
                H(1)=1.0
             END IF
          END IF
C..Fiber compressive failure
      ELSE
          IF(O(1).EQ.-1.0) THEN
             H(1)=-1.0
          ELSE
             LK=.5*(((ABS(S12)+(SL*COS(2*A0))*S22)/SL)+ABS((ABS(S12)+
     *            (SL*COS(2*A0))*S22)/SL))
             IF(LK.LE.-1.0) THEN
                H(1)=-1.0
             END IF
             H(1)=(1-M)*(S22/YT)+M*(S22/YT)**2+(S12/SL)**2
             IF(H(1).LE.-1.0) THEN
                H(1)=-2.0
             END IF
          END IF
      END IF
C..Matrix tensile failure
      IF(S22.GE.0) THEN
          IF(O(2).EQ.1.0) THEN
             H(2)=1.0
          ELSE
             H(2)=(1-M)*(S22/YT)+M*(S22/YT)**2+(S12/SL)**2
             IF(H(2).GE.1.0) THEN
                H(2)=1.0
             END IF
          END IF
C..Matrix compressive failure
      ELSE
          IF(O(2).EQ.-1.0) THEN
             H(2)=-1.0
          ELSE
             H(2)=(.5*(-S22*SIN(A)*COS(A)+ABS(-S22*SIN(A)*COS(A)))/ST)**
     *            2+(.5*(-S12*COS(A)+ABS(-S12*COS(A)))/SL)**2
             IF(H(2).LE.-1.0) THEN
                H(2)=-1.0
             END IF
          END IF
      END IF
C      
      RETURN
      END
[/plain]
And then the include is:
[plain]C===================================================================
C     NNMAX=max. node number in the model
C     NNTOT=max. no. of nodes in the model
C
C     NEMAX=max. element number in the model
C     NETOT=max. no. of elments in the model
C     NNELEMAX=max. no. nodes per element
C
C     MSECPT=max. no. of section points per element
C
C     NMODEMAX=max. no. of output modes
C
C     NCOMPMAX=max. no. of composite sections in the model
C     NLAYERMAX=max. no. of layer in each composite section
C
C     NTMAX=max. no. of time steps
C     NHARMAX=max. no. of harmonics
C
C     NSETMAX=max. no. of output requests per step
C===================================================================

      PARAMETER(NNMAX=2000000,NEMAX=2000000)
      PARAMETER(NNTOT=2000000,NETOT=2000000, NNELEMAX=8)
      PARAMETER(MSECPT=16)
      PARAMETER(NCOMPMAX=10,NLAYERMAX=2)
      PARAMETER(NTMAX=251,NHARMAX=30)
      PARAMETER(NMODEMAX=NTMAX*5,NSDVMAX=10)
      PARAMETER(NSETMAX=300)

      CHARACTER*80 FNAME, DUM

      COMMON/IALLELE/  NELEM,IREFE(NEMAX),IELEM(NETOT),
     1	               NDI,NSHR,NTYPE,INTPT,JELNUM,LOCATE
 
      COMMON/IALLNODE/ NNODE,IREFF(NNMAX),IREFS(NNMAX),INODE(NNTOT)
      COMMON/RALLNODE/ COORD(3,NNTOT)

      COMMON/IALLMODEL/NSTEP,NINC1,ISTEP

      COMMON/IPASS2/   NNDSET,NODEMAX,IELMAX,IN,
     1                 INODENUM,I101D,I102D,I103D,I104D,I108D

      COMMON/RALLMODEL/TTIME,STIME,TSET

      COMMON/RMAT/     XT,XC,YT,YC,SL,ST,MATL,IFAIL

      COMMON/RPTS/     ROOT,QTR,CTR,TIP,BCNR,TCNR
      
      DIMENSION        S(3),H(2)

      DATA ZERO,ONE,TWO,HALF,PI,HALFPI,TWOPI,HUNTPI,RAD/0.0,1.0,2.0,
     1     0.5,3.1415927,1.5707963,6.2831853,15.9154943,1.7453293E-02/

      COMMON/DYN1/     ACC(3),VEL(3),DIS(3),RXN(6)
      
      COMMON/FAIL/     HASH(2,NNMAX),OHASH(2,NNMAX),
     1                 CHAN(2,NNMAX),OCHAN(2,NNMAX),
     2                 CRSO(2,NNMAX),OCRSO(2,NNMAX),
     3                 CRSN(2,NNMAX),OCRSN(2,NNMAX),
     4                 PUCK(2,NNMAX),OPUCK(2,NNMAX),
     5                 KROL(2,NNMAX),OKROL(2,NNMAX),
     6                 LACR(2,NNMAX),OLACR(2,NNMAX),
     7                 FAIL(2,NNMAX),
     8                 PF(NNMAX),PB(NNMAX),
     9                 DCTR(NNMAX),VCTR(NNMAX),STR(6),LE(3)

      COMMON/RALLMODEL2/PRS,PEEQ,QTR1,CTR1
     1                  U1R,U2R,U3R,UR,U1Q,U2Q,U3Q,UQ,
     2                  U1C,U2C,U3C,UC,U1T,U2T,U3T,UT,
     3                  V1R,V2R,V3R,VR,V1Q,V2Q,V3Q,VQ,
     4                  V1C,V2C,V3C,VC,V1T,V2T,V3T,VT,
     5                  A1R,A2R,A3R,AR,A1Q,A2Q,A3Q,AQ,
     6                  A1C,A2C,A3C,AC,A1T,A2T,A3T,AT,
     7                  PR,PQ,PC,PT,PHI

[/plain]
Thanks in advance.
- R
0 Kudos

Link Copied

2 Replies
Lorri_M_Intel
Employee
‎11-20-2009 02:23 PM
596 Views

Well ... first, there are two include files in the original source, so there's not enough here to compile.

It also would make it much easier on us if you could point out which are the "failure" routines, and maybe give a clue as to exactly what doesn't work.

Sometimes I can polish up my crystal ball, but it seems to be in need of an oil change or something ...

thanks -

0 Kudos
Copy link
dannycat
New Contributor I
‎11-21-2009 04:05 AM
596 Views
Is this routineto be used withthe ABAQUS software?

All variables use default types which I would avoid in newly written code where possible. The trend these days is to ensure that arguments match up by type and the default compiler settings expect this to be the case. The arguments have not been dimensioned in the subroutines. You might get this to work by turning off argument checking in compiler diagnostics.

Does the same problem occur witheach failure criteria as only one is activated (uncommented)?

0 Kudos
Copy link
Reply

Community support is provided Monday to Friday. Other contact methods are available here.

Intel does not verify all solutions, including but not limited to any file transfers that may appear in this community. Accordingly, Intel disclaims all express and implied warranties, including without limitation, the implied warranties of merchantability, fitness for a particular purpose, and non-infringement, as well as any warranty arising from course of performance, course of dealing, or usage in trade.

For more complete information about compiler optimizations, see our Optimization Notice.