hm_read_fail_orthbiquad.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    hm_read_fail_orthbiquad    ../starter/source/materials/fail/orthbiquad/hm_read_fail_orthbiquad.f
!||--- called by ------------------------------------------------------
!||    hm_read_fail               ../starter/source/materials/fail/hm_read_fail.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                     ../starter/source/output/message/message.F
!||    hm_get_float_array_index   ../starter/source/devtools/hm_reader/hm_get_float_array_index.F
!||    hm_get_floatv              ../starter/source/devtools/hm_reader/hm_get_floatv.F
!||    hm_get_floatv_dim          ../starter/source/devtools/hm_reader/hm_get_floatv_dim.F
!||    hm_get_intv                ../starter/source/devtools/hm_reader/hm_get_intv.F
!||    hm_option_is_encrypted     ../starter/source/devtools/hm_reader/hm_option_is_encrypted.F
!||--- uses       -----------------------------------------------------
!||    hm_option_read_mod         ../starter/share/modules1/hm_option_read_mod.f
!||    message_mod                ../starter/share/message_module/message_mod.F
!||    submodel_mod               ../starter/share/modules1/submodel_mod.f
!||====================================================================
      SUBROUTINE hm_read_fail_orthbiquad(
     .           FAIL    ,MAT_ID   ,FAIL_ID  ,IRUPT    ,
     .           TITR    ,LSUBMODEL,UNITAB   )
C-----------------------------------------------
c    ROUTINE DESCRIPTION :
c    Orthotropic strain  failure model
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE fail_param_mod
      USE unitab_mod
      USE message_mod
      USE submodel_mod
      USE hm_option_read_mod
      USE names_and_titles_mod , ONLY : nchartitle
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C----------+---------+---+---+--------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "units_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER            ,INTENT(IN) :: FAIL_ID         ! failure model ID
      INTEGER            ,INTENT(IN) :: MAT_ID          ! material law ID
      INTEGER            ,INTENT(IN) :: IRUPT           ! failure model type number
      CHARACTER(LEN=NCHARTITLE) ,INTENT(IN) :: TITR     ! material model title
      TYPE(unit_type_)   ,INTENT(IN) :: UNITAB          ! table of input units
      TYPE(submodel_data),INTENT(IN) :: LSUBMODEL(*)    ! submodel table
      TYPE(fail_param_)  ,INTENT(INOUT) :: FAIL         ! failure model data structure
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: NANGLE,I,J,K,INFO,REG_FUNC,MFLAG,SFLAG,RATE_FUNC,NFUNC,NUPARAM,NUVAR
      INTEGER :: IPIV2(2),IPIV3(3)
      INTEGER ,PARAMETER :: NSIZE = 2
      INTEGER ,DIMENSION(NSIZE) :: IFUNC
      my_real :: pthk,ref_siz,ref_siz_unit,epsd0,cjc,rate_scale,ref_rate_unit,
     .    r1,r2,r4,r5,c5,c5_min,theta_myreal
      my_real, DIMENSION(:), ALLOCATABLE :: c1,c2,c3,c4,inst
      DOUBLE PRECISION A_1(2,2),B_1(2),A_2(3,3),B_2(3),
     .                 triax_1_lin,triax_2_lin,triax_3_lin,
     .                 triax_4_lin,triax_5_lin,triax_1_quad,
     .                 triax_2_quad,triax_3_quad,triax_4_quad,
     .                 triax_5_quad,cos2(10,10),xmin,ymin
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: THETA,THETA_RAD,Q_X11,Q_X12,Q_X13,
     .                                               Q_X21,Q_X22,Q_X23,Q_INST
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: X_1,X_2,AMAT,BVEC
      INTEGER, DIMENSION(:), ALLOCATABLE :: IPIV
      LOGICAL :: IS_AVAILABLE,IS_ENCRYPTED
      DATA triax_1_lin, triax_2_lin, triax_3_lin, triax_4_lin,
     .     triax_5_lin 
     .      / -0.33333333, 0.0, 0.33333333, 0.577350269, 0.66666667 /
      DATA triax_1_quad, triax_2_quad, triax_3_quad, 
     .     triax_4_quad, triax_5_quad 
     .      / 0.11111111, 0.0, 0.11111111, 0.33333333, 0.44444444 /
C
      DATA  cos2/
     1 1.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,   
     2 0.   ,1.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     3 -1.  ,0.   ,2.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     4 0.   ,-3.  ,0.   ,4.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     5 1.   ,0.   ,-8.  ,0.   ,8.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     6 0.   ,5.   ,0.   ,-20. ,0.   ,16.  ,0.   ,0.   ,0.   ,0.   ,     
     7 -1.  ,0.   ,18.  ,0.   ,-48. ,0.   ,32.  ,0.   ,0.   ,0.   ,     
     8 0.   ,-7.  ,0.   ,56.  ,0.   ,-112.,0.   ,64.  ,0.   ,0.   ,     
     9 1.   ,0.   ,-32. ,0.   ,160. ,0.   ,-256.,0.   ,128. ,0.   ,     
     a 0.   ,9.   ,0.   ,-120.,0.   ,432. ,0.   ,-576 ,0.   ,256. /
C=======================================================================
      is_encrypted = .false.
      is_available = .false.
C--------------------------------------------------
C     (IS OPTION CRYPTED)
C--------------------------------------------------
 
      CALL hm_option_is_encrypted(is_encrypted)
 
C======================================================================================
C     EXTRACT DATA
!---------------
! -> Card1
!---------------
! Percentage of integration failure, flags and size regularization
      CALL hm_get_floatv ('Pthk'           ,pthk   ,is_available,lsubmodel,unitab)
      CALL hm_get_intv   ('MAT_MFLAG'      ,mflag  ,is_available,lsubmodel)
      CALL hm_get_intv   ('MAT_SFLAG'      ,sflag  ,is_available,lsubmodel)
      CALL hm_get_intv   ('MAT_refanglemax',nangle ,is_available,lsubmodel)
      ! Error message
      IF (nangle > 10) THEN 
        CALL ancmsg(msgid=2015,msgtype=msgerror,
     .              anmode=aninfo_blind_1,i1=mat_id,c1=titr)      
      ENDIF 
      CALL hm_get_intv   ('fct_IDel' ,reg_func    ,is_available,lsubmodel)
      CALL hm_get_floatv ('EI_ref'   ,ref_siz     ,is_available,lsubmodel,unitab)
! Default values      
      IF (pthk == zero) pthk = one - em06
      pthk = min(pthk, one)
      pthk = max(pthk,-one)
      IF (sflag == 0)    sflag = 2
! Units
      IF ((ref_siz == zero).AND.(reg_func > 0)) THEN
        CALL hm_get_floatv_dim('EI_ref' ,ref_siz_unit,is_available, lsubmodel, unitab)
        ref_siz = one*ref_siz_unit
      ENDIF
!---------------
! -> Card2
!---------------
! Size regularization and percentage of integration failure
      CALL hm_get_floatv ('MAT_C5'    ,c5          ,is_available,lsubmodel,unitab)
      CALL hm_get_floatv ('MAT_EPSD0' ,epsd0       ,is_available,lsubmodel,unitab)
      CALL hm_get_floatv ('MAT_CJC'   ,cjc         ,is_available,lsubmodel,unitab) 
      CALL hm_get_intv   ('fct_IDrate',rate_func   ,is_available,lsubmodel)
      CALL hm_get_floatv ('RATE_scale',rate_scale  ,is_available,lsubmodel,unitab)
! Units and default values
      IF ((rate_scale == zero).AND.(rate_func > 0)) THEN 
        CALL hm_get_floatv_dim('RATE_scale' ,ref_rate_unit  ,is_available, lsubmodel, unitab)
        rate_scale = ref_rate_unit*one
      ENDIF
      IF (rate_func > 0) THEN 
        cjc   = zero
        epsd0 = zero
      ELSE
        rate_scale = zero
      ENDIF
      IF (cjc == zero .OR. epsd0 == zero) THEN 
        cjc   = zero
        epsd0 = zero
      ENDIF
!---------------
! -> Card3
!---------------
! Biquad parameter by angle
      IF (.NOT.ALLOCATED(c1))   ALLOCATE(c1(nangle))
      IF (.NOT.ALLOCATED(c2))   ALLOCATE(c2(nangle))
      IF (.NOT.ALLOCATED(c3))   ALLOCATE(c3(nangle)) 
      IF (.NOT.ALLOCATED(c4))   ALLOCATE(c4(nangle))
      IF (sflag == 3) THEN 
        IF (.NOT.ALLOCATED(inst)) ALLOCATE(inst(nangle))
        inst = zero
      ENDIF
      c5_min = infinity
      ! Material selector
      IF (mflag == 0) THEN 
        ! Loop over angles (must be equally distributed between 0 and pi/2)
        DO j = 1, nangle
          CALL hm_get_float_array_index('MAT_C1',c1(j),j,is_available,lsubmodel,unitab)
          CALL hm_get_float_array_index('mat_c2',C2(J),J,IS_AVAILABLE,LSUBMODEL,UNITAB)
          CALL HM_GET_FLOAT_ARRAY_INDEX('mat_c3',C3(J),J,IS_AVAILABLE,LSUBMODEL,UNITAB)
          CALL HM_GET_FLOAT_ARRAY_INDEX('mat_c4',C4(J),J,IS_AVAILABLE,LSUBMODEL,UNITAB)
          ! Default values
          IF (C3(J) == ZERO) C3(J) = 0.6D0
.AND..AND..AND.          IF (C1(J) == ZERO  C2(J) == ZERO  C4(J) == ZERO  C5 == ZERO) THEN
            C1(J)  = 3.5D0*C3(J)
            C2(J)  = 1.6D0*C3(J)
            C4(J)  = 0.6D0*C3(J)
            C5_MIN = MIN(C5_MIN,1.5D0*C3(J))
          ENDIF
          ! If necking instability is activated
          IF (SFLAG == 3) THEN
            CALL HM_GET_FLOAT_ARRAY_INDEX('mat_inst',INST(J),J,IS_AVAILABLE,LSUBMODEL,UNITAB)
            IF (INST(J) <= ZERO)   THEN
              CALL ANCMSG(MSGID=2016,MSGTYPE=MSGWARNING,
     .                    ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR)    
              SFLAG = 2
            ELSEIF (INST(J) >= C4(J)) THEN
              CALL ANCMSG(MSGID=2017,MSGTYPE=MSGWARNING,
     .                    ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR)
              SFLAG = 2
            ENDIF           
          ENDIF
        ENDDO
      ELSE
        ! RX parameters
        IF (MFLAG == 1) THEN ! Mild Seel (c3 = 0.6)
          R1 = 3.5D0 
          R2 = 1.6D0 
          R4 = 0.6D0 
          R5 = 1.5D0
        ELSEIF (MFLAG == 2) THEN ! DP600 (c3 = 0.5)
          R1 = 4.3D0 
          R2 = 1.4D0
          R4 = 0.6D0
          R5 = 1.6D0
        ELSEIF (MFLAG == 3) THEN ! Boron (c3 = 0.12)
          R1 = 5.2D0 
          R2 = 3.1D0 
          R4 = 0.8D0 
          R5 = 3.5D0 
        ELSEIF (MFLAG == 4) THEN ! Aluminium AA5182 (c3 = 0.3)
          R1 = 5.0D0
          R2 = 1.0D0 
          R4 = 0.4D0
          R5 = 0.8D0
        ELSEIF (MFLAG == 5) THEN ! Aluminium AA6082-T6 (c3 = 0.17)
          R1 = 7.8D0
          R2 = 3.5D0   
          R4 = 0.6D0
          R5 = 2.8D0
        ELSEIF (MFLAG == 6) THEN ! Plastic light_eBody PA6GF30 (c3 = 0.1)
          R1 = 3.6D0
          R2 = 0.6D0
          R4 = 0.5D0
          R5 = 0.6D0
        ELSEIF (MFLAG == 7) THEN ! Plastic light_eBody PP T40 ( c3=0.11 )
          R1 = 10.0D0
          R2 = 2.7D0
          R4 = 0.6D0
          R5 = 0.7D0
        ELSEIF (MFLAG == 99) THEN ! user scalling factors
          CALL HM_GET_FLOATV ('mat_r1' ,R1    ,IS_AVAILABLE,LSUBMODEL,UNITAB)
          CALL HM_GET_FLOATV ('mat_r2' ,R2    ,IS_AVAILABLE,LSUBMODEL,UNITAB)
          CALL HM_GET_FLOATV ('mat_r4' ,R4    ,IS_AVAILABLE,LSUBMODEL,UNITAB)
          CALL HM_GET_FLOATV ('mat_r5' ,R5    ,IS_AVAILABLE,LSUBMODEL,UNITAB)
        ELSE ! ELSE --> Mild Seel
          R1 = 3.5D0
          R2 = 1.6D0
          R4 = 0.6D0
          R5 = 1.5D0
        ENDIF   
        ! Loop over angles (must be equally distributed between 0 and pi/2)
        DO J = 1, NANGLE
          CALL HM_GET_FLOAT_ARRAY_INDEX('mat_c3'  ,C3(J)  ,J,IS_AVAILABLE,LSUBMODEL,UNITAB)
          ! Default values 
          IF (C3(J) == ZERO) THEN 
            IF (MFLAG == 1) THEN 
              C3(J) = 0.6D0
            ELSEIF (MFLAG == 2) THEN 
              C3(J) = 0.5D0
            ELSEIF (MFLAG == 3) THEN
              C3(J) = 0.12D0
            ELSEIF (MFLAG == 4) THEN
              C3(J) = 0.3D0
            ELSEIF (MFLAG == 5) THEN
              C3(J) = 0.17D0
            ELSEIF (MFLAG == 6) THEN
              C3(J) = 0.1D0
            ELSEIF (MFLAG == 7) THEN
              C3(J) = 0.11D0
            ENDIF
          ENDIF
          ! Computation of C1,C2,C4,C5
          C1(J)  = R1*C3(J)
          C2(J)  = R2*C3(J)
          C4(J)  = R4*C3(J)
          C5_MIN = MIN(C5_MIN,R5*C3(J))
          ! If necking instability is activated
          IF (SFLAG == 3) THEN
            CALL HM_GET_FLOAT_ARRAY_INDEX('mat_inst',INST(J),J,IS_AVAILABLE,LSUBMODEL,UNITAB)
            IF (INST(J) <= ZERO)   THEN
              CALL ANCMSG(MSGID=2016,MSGTYPE=MSGWARNING,
     .                    ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR)    
              SFLAG = 2
            ELSEIF (INST(J) >= C4(J)) THEN
              CALL ANCMSG(MSGID=2017,MSGTYPE=MSGWARNING,
     .                    ANMODE=ANINFO_BLIND_1,I1=MAT_ID,C1=TITR)
              SFLAG = 2
            ENDIF           
          ENDIF
        ENDDO
      ENDIF
      ! Default value
      IF (C5 == ZERO) C5 = C5_MIN
C======================================================================================
C    COMPUTING FITTING PARAMETERS
C======================================================================================
.NOT.      IF (ALLOCATED(X_1))  ALLOCATE(X_1(NANGLE,3))
.NOT.      IF (ALLOCATED(X_2))  ALLOCATE(X_2(NANGLE,3))
      !================================================
      ! Loop over the angle
      !================================================
      DO J = 1,NANGLE
C
        ! Coefficient for the first parabole
        ! ---------------------------------------
        A_1(1,1) = TRIAX_1_LIN
        A_1(1,2) = TRIAX_1_QUAD
        A_1(2,1) = TRIAX_3_LIN
        A_1(2,2) = TRIAX_3_QUAD
        B_1(1)   = C1(J) - C2(J)
        B_1(2)   = C3(J) - C2(J)
C
        ! Fitting the first quadratic function
#ifndef WITHOUT_LINALG
        CALL DGESV(2, 1, A_1, 2, IPIV2, B_1, 2, INFO)  
#else
        WRITE(6,*) "Error: Blas/Lapack required"
#endif
        X_1(J,1)   = C2(J)
        X_1(J,2:3) = B_1(1:2)
C      
        ! Coefficient for the second parabole
        !----------------------------------------
        A_2(1,1) = ONE
        A_2(1,2) = TRIAX_3_LIN
        A_2(1,3) = TRIAX_3_QUAD
        A_2(2,1) = ONE
        A_2(2,2) = TRIAX_4_LIN
        A_2(2,3) = TRIAX_4_QUAD
        A_2(3,1) = ONE
        A_2(3,2) = TRIAX_5_LIN
        A_2(3,3) = TRIAX_5_QUAD
        B_2(1)   = C3(J)
        B_2(2)   = C4(J)
        B_2(3)   = C5
C
        ! Fitting the second quadratic function
#ifndef WITHOUT_LINALG
        CALL DGESV(3, 1, A_2, 3, IPIV3, B_2, 3, INFO)
#endif
        X_2(J,1:3) = B_2(1:3)
C
      ENDDO
C
C======================================================================================
C    COMPUTING COSINE INTERPOLATION
C======================================================================================
c      
.NOT.      IF (ALLOCATED(THETA))     ALLOCATE(THETA(NANGLE))
.NOT.      IF (ALLOCATED(THETA_RAD)) ALLOCATE(THETA_RAD(NANGLE))
c
      ! Computation of angles and check curves
      DO J = 1, NANGLE
        IF (NANGLE > 1) THEN 
          THETA(J)     = (J-1)*(90.0D0/(NANGLE-1))
          THETA_RAD(J) = THETA(J)*(PI/180.0D0)
        ELSE
          THETA(J)     = ZERO
          THETA_RAD(J) = ZERO
        ENDIF
c
        ! Check if minimum of first parabola is negative 
        XMIN = -X_1(J,2)/(TWO*X_1(J,3))
        YMIN = X_1(J,3)*(XMIN**2) + X_1(J,2)*XMIN + X_1(J,1)
        IF (YMIN < ZERO) THEN 
          THETA_MYREAL = THETA(J)
          CALL ANCMSG(MSGID=3002,
     .                MSGTYPE=MSGWARNING,
     .                ANMODE=ANINFO_BLIND_1,
     .                I1=MAT_ID,
     .                C1=TITR,
     .                I2=J,
     .                R1=THETA_MYREAL)          
        ENDIF
c
        ! Check if minimum of second parabola is negative 
        IF (SFLAG == 1) THEN 
          XMIN = -X_2(J,2)/(TWO*X_2(J,3))
          YMIN = X_2(J,3)*(XMIN**2) + X_2(J,2)*XMIN + X_2(J,1)
          IF (YMIN < ZERO) THEN 
            THETA_MYREAL = THETA(J)
            CALL ANCMSG(MSGID=3003,
     .                  MSGTYPE=MSGWARNING,
     .                  ANMODE=ANINFO_BLIND_1,
     .                  I1=MAT_ID,
     .                  C1=TITR,
     .                  I2=J,
     .                  R1=THETA_MYREAL)              
          ENDIF
        ENDIF
c
      ENDDO
c
      ! Allocation of the A-MATRIX and the Pivot vector
      ALLOCATE (AMAT(NANGLE,NANGLE),IPIV(NANGLE))
c
      ! Filling the A-MATRIX
      DO J = 1,NANGLE
        DO I = 1,NANGLE
          AMAT(J,I) = ZERO
          DO K = 1,I
            AMAT(J,I) = AMAT(J,I) + COS2(K,I)*(COS(TWO*THETA_RAD(J)))**(K-1)
          ENDDO
        ENDDO
      ENDDO
c
      ! Allocation of factors
      ALLOCATE(Q_X11(NANGLE),Q_X12(NANGLE),Q_X13(NANGLE),
     .         Q_X21(NANGLE),Q_X22(NANGLE),Q_X23(NANGLE))
c
      ! Initialization of tables
      Q_X11(1:NANGLE) = ZERO
      Q_X12(1:NANGLE) = ZERO
      Q_X13(1:NANGLE) = ZERO
      Q_X21(1:NANGLE) = ZERO
      Q_X22(1:NANGLE) = ZERO
      Q_X23(1:NANGLE) = ZERO
      IF (SFLAG == 3) THEN 
        ALLOCATE(Q_INST(NANGLE))
        Q_INST(1:NANGLE) = ZERO
      ENDIF
c
      ! Filling the B vector with experimental points
      IF (SFLAG == 3) THEN 
        ALLOCATE (BVEC(NANGLE,7))
      ELSE
        ALLOCATE (BVEC(NANGLE,6))
      ENDIF
      BVEC(1:NANGLE,1) = X_1(1:NANGLE,1)
      BVEC(1:NANGLE,2) = X_1(1:NANGLE,2)
      BVEC(1:NANGLE,3) = X_1(1:NANGLE,3)
      BVEC(1:NANGLE,4) = X_2(1:NANGLE,1)
      BVEC(1:NANGLE,5) = X_2(1:NANGLE,2)
      BVEC(1:NANGLE,6) = X_2(1:NANGLE,3)
      IF (SFLAG == 3) BVEC(1:NANGLE,7) = INST(1:NANGLE)
c        
      ! Initialization of the Pivot vector
      IPIV(1:NANGLE) = 0
c
      ! Solving the A-MATRIX * x = B vector system
#ifndef WITHOUT_LINALG
      IF (SFLAG == 3) THEN 
        CALL DGESV(NANGLE, 7, AMAT, NANGLE, IPIV, BVEC, NANGLE, INFO)
      ELSE
        CALL DGESV(NANGLE, 6, AMAT, NANGLE, IPIV, BVEC, NANGLE, INFO)
      ENDIF
#else
      WRITE(6,*) "Error: Blas/Lapack required"
#endif
c        
      ! Recovering the solution
      Q_X11(1:NANGLE)  = BVEC(1:NANGLE,1)
      Q_X12(1:NANGLE)  = BVEC(1:NANGLE,2)
      Q_X13(1:NANGLE)  = BVEC(1:NANGLE,3)
      Q_X21(1:NANGLE)  = BVEC(1:NANGLE,4)
      Q_X22(1:NANGLE)  = BVEC(1:NANGLE,5)
      Q_X23(1:NANGLE)  = BVEC(1:NANGLE,6)
      IF (SFLAG == 3) Q_INST(1:NANGLE) = BVEC(1:NANGLE,7)
c
c----------------------------------
      ! -> Number of parameters
      NUPARAM = 7
      IF (SFLAG == 3) THEN 
        NUPARAM = NUPARAM + 7*NANGLE
      ELSE
        NUPARAM = NUPARAM + 6*NANGLE
      ENDIF
      ! -> Number of functions
      NFUNC   = 0
      IF (RATE_FUNC /= 0) THEN 
        NFUNC = NFUNC + 1
        IFUNC(NFUNC) = RATE_FUNC
      ENDIF
      IF (REG_FUNC /= 0) THEN 
        NFUNC = NFUNC + 1 
        IFUNC(NFUNC) = REG_FUNC
      ENDIF
      ! -> Number of user variables
      NUVAR = 3   
C======================================================================================
c     Filling buffer tables
C======================================================================================
      FAIL%KEYWORD = 'orth-biquad' 
      FAIL%IRUPT   = IRUPT 
      FAIL%FAIL_ID = FAIL_ID 
      FAIL%NUPARAM = NUPARAM
      FAIL%NIPARAM = 0
      FAIL%NUVAR   = NUVAR
      FAIL%NFUNC   = NFUNC
      FAIL%NTABLE  = 0
      FAIL%NMOD    = 0
      FAIL%PTHK    = PTHK
c            
      ALLOCATE (FAIL%UPARAM(FAIL%NUPARAM))
      ALLOCATE (FAIL%IPARAM(FAIL%NIPARAM))
      ALLOCATE (FAIL%IFUNC (FAIL%NFUNC))
      ALLOCATE (FAIL%TABLE (FAIL%NTABLE))
c
      FAIL%IFUNC(1:NFUNC) = IFUNC(1:NFUNC)
c
      FAIL%UPARAM(1) = PTHK
      FAIL%UPARAM(2) = SFLAG
      FAIL%UPARAM(3) = REF_SIZ
      FAIL%UPARAM(4) = EPSD0
      FAIL%UPARAM(5) = CJC
      FAIL%UPARAM(6) = RATE_SCALE
      FAIL%UPARAM(7) = NANGLE
      IF (SFLAG == 3) THEN 
        DO J = 1,NANGLE
          FAIL%UPARAM(8  + 7*(J-1)) = Q_X11(J)
          FAIL%UPARAM(9  + 7*(J-1)) = Q_X12(J)
          FAIL%UPARAM(10 + 7*(J-1)) = Q_X13(J)
          FAIL%UPARAM(11 + 7*(J-1)) = Q_X21(J)
          FAIL%UPARAM(12 + 7*(J-1)) = Q_X22(J)
          FAIL%UPARAM(13 + 7*(J-1)) = Q_X23(J)
          FAIL%UPARAM(14 + 7*(J-1)) = Q_INST(J)
        ENDDO    
      ELSE
        DO J = 1,NANGLE
          FAIL%UPARAM(8  + 6*(J-1)) = Q_X11(J)
          FAIL%UPARAM(9  + 6*(J-1)) = Q_X12(J)
          FAIL%UPARAM(10 + 6*(J-1)) = Q_X13(J)
          FAIL%UPARAM(11 + 6*(J-1)) = Q_X21(J)
          FAIL%UPARAM(12 + 6*(J-1)) = Q_X22(J)
          FAIL%UPARAM(13 + 6*(J-1)) = Q_X23(J)
        ENDDO          
      ENDIF
c--------------------------
c     Printout data
c-------------------------- 
      IF (IS_ENCRYPTED) THEN
        WRITE(IOUT,'(5x,a,//)')'confidential data'
      ELSE
        WRITE(IOUT,1000)
        IF (MFLAG /= 0) WRITE(IOUT, 1100) MFLAG
        DO J = 1,NANGLE
          WRITE(IOUT,1200) THETA(J),C1(J),C2(J),C3(J),C4(J),C5,
     &    X_1(J,3),X_1(J,2),X_1(J,1),X_2(J,3),X_2(J,2),X_2(J,1)           
          IF (SFLAG == 3) WRITE(IOUT, 1900) INST(J)
        ENDDO
        WRITE(IOUT,1300) PTHK,SFLAG
        IF (REG_FUNC > 0) WRITE(IOUT,1400) REG_FUNC,REF_SIZ
        IF (EPSD0 > ZERO) THEN 
          WRITE(IOUT,1500) EPSD0,CJC
        ELSEIF (RATE_FUNC > 0) THEN 
          WRITE(IOUT,1600) RATE_FUNC,RATE_SCALE
        ENDIF
        WRITE(IOUT,2000)
      ENDIF
c--------------------------
c     Deallocation
c--------------------------  
      IF (ALLOCATED(C1))        DEALLOCATE(C1)
      IF (ALLOCATED(C2))        DEALLOCATE(C2)
      IF (ALLOCATED(C3))        DEALLOCATE(C3)
      IF (ALLOCATED(C4))        DEALLOCATE(C4)
      IF (ALLOCATED(INST))      DEALLOCATE(INST)
      IF (ALLOCATED(X_1))       DEALLOCATE(X_1)
      IF (ALLOCATED(X_2))       DEALLOCATE(X_2)
      IF (ALLOCATED(THETA))     DEALLOCATE(THETA)
      IF (ALLOCATED(THETA_RAD)) DEALLOCATE(THETA_RAD)
      IF (ALLOCATED(Q_X11))     DEALLOCATE(Q_X11)
      IF (ALLOCATED(Q_X12))     DEALLOCATE(Q_X12)
      IF (ALLOCATED(Q_X13))     DEALLOCATE(Q_X13)
      IF (ALLOCATED(Q_X21))     DEALLOCATE(Q_X21)
      IF (ALLOCATED(Q_X22))     DEALLOCATE(Q_X22)
      IF (ALLOCATED(Q_X23))     DEALLOCATE(Q_X23)
      IF (ALLOCATED(Q_INST))    DEALLOCATE(Q_INST)
      IF (ALLOCATED(AMAT))      DEALLOCATE(AMAT)
      IF (ALLOCATED(IPIV))      DEALLOCATE(IPIV)
c-----------------------------------------------------
 1000 FORMAT(
     & 5X,'  ------------------------------------------   ',/
     & 5X,'    failure criterion : orthotropic biquad     ',/,
     & 5X,'  ------------------------------------------   ',/)
 1100 FORMAT(
     & 5X,'material PARAMETER selector  . . . . . . . . . . .=',I10)
 1200 FORMAT(     
     & 5X,'|| failure strains for angle',F5.1,' deg',/,
     & 5X,'  -------------------------------------------------',/,
     & 5X,'  simple compression c1  . . . . . . . . . . . . .=',1PG20.13/
     & 5X,'  shear c2 . . . . . . . . . . . . . . . . . . . .=',1PG20.13/
     & 5X,'  simple tension c3  . . . . . . . . . . . . . . .=',1PG20.13/
     & 5X,'  plane strain c4  . . . . . . . . . . . . . . . .=',1PG20.13/
     & 5X,'  biaxial tension c5 . . . . . . . . . . . . . . .=',1PG20.13/
     & 5X,'                                                   ',/
     & 5X,'  low stress triaxiality parabola PARAMETER a. . .=',1PG20.13/ 
     & 5X,'  low stress triaxiality parabola PARAMETER b. . .=',1PG20.13/ 
     & 5X,'  low stress triaxiality parabola PARAMETER c. . .=',1PG20.13/ 
     & 5X,'                                                   ',/
     & 5X,'  high stress triaxiality parabola PARAMETER d . .=',1PG20.13/ 
     & 5X,'  high stress triaxiality parabola PARAMETER e . .=',1PG20.13/ 
     & 5X,'  high stress triaxiality parabola PARAMETER f . .=',1PG20.13/)
 1300 FORMAT(
     & 5X,'element deletion :',/,
     & 5X,'shell element deletion PARAMETER pthickfail. . . .=',1PG20.13,/,
     & 5X,'  > 0.0 : fraction of failed thickness             ',/,
     & 5X,'  < 0.0 : fraction of failed intg. points or layers',/,
     & 5X,'s-flag . . . . . . . . . . . . . . . . . . . . . .=',I10,/)
 1400 FORMAT(
     & 5X,'element length regularization used:',/,
     & 5X,'regularization FUNCTION id . . . . . . . . . . . .=',I10,/,
     & 5X,'reference element length . . . . . . . . . . . . .=',1PG20.13,/)
 1500 FORMAT(
     & 5X,'johnson-cook strain-rate dependency:',/,
     & 5X,'reference strain-rate  . . . . . . . . . . . . . .=',1PG20.13,/,
     & 5X,'johnson-cook parameter . . . . . . . . . . . . . .=',1PG20.13,/)
 1600 FORMAT(
     & 5X,'tabulated strain-rate dependency:',/,
     & 5X,'strain-rate dependency function id . . . . . . . .=',I10,/,
     & 5X,'strain-rate scale factor . . . . . . . . . . . . .=',1PG20.13,/)
 1900 FORMAT(
     & 5X,'  instability strain . . . . . . . . . . . . . . .=',1PG20.13,//)
 2000 FORMAT(
     & 5X,' --------------------------------------------------',//)
c----------- 
      END