s6cderi3.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/.
!||====================================================================
!||    s6cderi3      ../engine/source/elements/thickshell/solide6c/s6cderi3.f
!||--- called by ------------------------------------------------------
!||    s6cforc3      ../engine/source/elements/thickshell/solide6c/s6cforc3.F
!||    s6cke3        ../engine/source/elements/thickshell/solide6c/s6cke3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg        ../engine/source/output/message/message.F
!||--- uses       -----------------------------------------------------
!||    message_mod   ../engine/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE s6cderi3(
     1   OFF,     DET,     NGL,     X1,
     2   X2,      X3,      X4,      X5,
     3   X6,      Y1,      Y2,      Y3,
     4   Y4,      Y5,      Y6,      Z1,
     5   Z2,      Z3,      Z4,      Z5,
     6   Z6,      PX1,     PX2,     PX3,
     7   PX4,     PY1,     PY2,     PY3,
     8   PY4,     PZ1,     PZ2,     PZ3,
     9   PZ4,     PX1H,    PX2H,    PX3H,
     A   PY1H,    PY2H,    PY3H,    PZ1H,
     B   PZ2H,    PZ3H,    JACOB5,  JACOB6,
     C   JACOB4,  JACOB8,  JACOB9,  JACOB7,
     D   JACI33,  B1X,     B1Y,     B2Y,
     E   B2X,     B1122,   B1221,   B2212,
     F   B1121,   B1XH,    B1YH,    B2XH,
     G   B2YH,    B1122H,  B1221H,  B2212H,
     H   B1121H,  VZL,     VOLG,    SAV,
     I   OFFG,    NEL,     ISMSTR)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com06_c.inc"
#include      "units_c.inc"
#include      "scr07_c.inc"
#include      "scr17_c.inc"
#include      "scr18_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: ISMSTR
      INTEGER :: NEL, NGL(*)
C
      my_real
     .   OFF(*),DET(*),
     .   X1(*), X2(*), X3(*), X4(*), X5(*), X6(*), 
     .   Y1(*), Y2(*), Y3(*), Y4(*), Y5(*), Y6(*),  
     .   Z1(*), Z2(*), Z3(*), Z4(*), Z5(*), Z6(*),   
     .   PX1(*), PX2(*), PX3(*), PX4(*),  
     .   PY1(*), PY2(*), PY3(*), PY4(*),  
     .   PZ1(*), PZ2(*), PZ3(*), PZ4(*),  
     .   PX1H(*), PX2H(*), PX3H(*),  
     .   PY1H(*), PY2H(*), PY3H(*),   
     .   PZ1H(*), PZ2H(*), PZ3H(*),   
     .   JACOB7(*),JACOB8(*),JACOB9(*),
     .   JACOB4(*),JACOB5(*),JACOB6(*),
     .   JACI33(*),B1X(MVSIZ,2),B1Y(MVSIZ,2),B2X(MVSIZ,2),B2Y(MVSIZ,2),
     .   B1XH(MVSIZ,2),B1YH(MVSIZ,2),B2XH(MVSIZ,2),B2YH(MVSIZ,2),
     .   b1122(*),b1221(*),b2212(*),b1121(*),
     .   b1122h(*),b1221h(*),b2212h(*),b1121h(*),
     .   vzl(*),volg(*),offg(*)
      double precision
     .   sav(nel,15)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I, J,ICOR,NNEGA,INDEX(MVSIZ)
C                                                                     12
      my_real
     .   DETT(MVSIZ), JAC1(MVSIZ), JAC2(MVSIZ), JAC3(MVSIZ), 
     .   JAC4(MVSIZ), JAC5(MVSIZ), JAC6(MVSIZ),
     .   JAC7(MVSIZ), JAC8(MVSIZ), JAC9(MVSIZ),
     .   JACI1, JACI2, JACI3,
     .   JACI4, JACI5, JACI6,
     .   JACI7, JACI8, JACI9,
     .   X21(MVSIZ) , X31(MVSIZ) , X54(MVSIZ) , X64(MVSIZ),
     .   Y21(MVSIZ) , Y31(MVSIZ) , Y54(MVSIZ) , Y64(MVSIZ),
     .   Z21(MVSIZ) , Z31(MVSIZ) , Z54(MVSIZ) , Z64(MVSIZ),
     .   X41(MVSIZ) , Y41(MVSIZ) , Z41(MVSIZ) , 
     .   jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz),
     .   jaci12, jaci45, jaci78,
     .   fac
C-----------------------------------------------
      nnega=0
C
      DO i=1,nel
      x21(i)=x2(i)-x1(i)
      x31(i)=x3(i)-x1(i)
      x41(i)=x4(i)-x1(i)
      x54(i)=x5(i)-x4(i)
      x64(i)=x6(i)-x4(i)
      y21(i)=y2(i)-y1(i)
      y31(i)=y3(i)-y1(i)
      y41(i)=y4(i)-y1(i)
      y54(i)=y5(i)-y4(i)
      y64(i)=y6(i)-y4(i)
      z21(i)=z2(i)-z1(i)
      z31(i)=z3(i)-z1(i)
      z41(i)=z4(i)-z1(i)
      z54(i)=z5(i)-z4(i)
      z64(i)=z6(i)-z4(i)
      END DO
C
C Jacobian matrix
      DO i=1,nel
C  -------ri.xi---->ksi--------
       jac1(i)=x21(i)+x54(i)
       jac2(i)=y21(i)+y54(i)
       jac3(i)=z21(i)+z54(i)
      END DO
      DO i=1,nel
C  -------si.xi--->eta--------
       jac4(i)=x31(i)+x64(i)
       jac5(i)=y31(i)+y64(i)
       jac6(i)=z31(i)+z64(i)
C  -------ti.xi----zeta-------
       jac7(i)=third*(x41(i)+x5(i)-x2(i)+x6(i)-x3(i))
       jac8(i)=third*(y41(i)+y5(i)-y2(i)+y6(i)-y3(i))
       jac9(i)=third*(z41(i)+z5(i)-z2(i)+z6(i)-z3(i))
      END DO
C
      DO i=1,nel
       jacob4(i)=jac4(i)
       jacob5(i)=jac5(i)
       jacob6(i)=jac6(i)
       jacob7(i)=jac7(i)
       jacob8(i)=jac8(i)
       jacob9(i)=jac9(i)
      ENDDO 
C
      DO i=1,nel
      jac_59_68(i)=jac5(i)*jac9(i)-jac6(i)*jac8(i)
      jac_67_49(i)=jac6(i)*jac7(i)-jac4(i)*jac9(i)
      jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
      END DO
C
      DO i=1,nel
      det(i)=one_over_8*(jac1(i)*jac_59_68(i)+jac2(i)*jac_67_49(i)+jac3(i)*jac_48_57(i))
      END DO
C
      IF(idtmin(1)==1)THEN
        icor = 0
        DO i=1,nel
           IF(off(i) ==zero)THEN
            det(i)=one
          ELSEIF((det(i)<=volmin).OR.(det(i)<=zero))THEN
            icor = 1
          ENDIF
        ENDDO
        IF (icor>0) THEN
          DO i=1,nel
           IF(off(i)/=zero)THEN
             IF(det(i)<=volmin)THEN
              det(i)=one
              off(i)=zero
#include "lockon.inc"
                WRITE(istdo,2000) ngl(i)
                WRITE(iout ,2000) ngl(i)
#include "lockoff.inc"
             ELSEIF(det(i)<=zero)THEN
               CALL ancmsg(msgid=166,anmode=aninfo,
     .                     i1=ngl(i))
               mstop = 1
             ENDIF
           ENDIF
          ENDDO
        ENDIF
      ELSEIF(idtmin(1)==2)THEN
        icor = 0
        DO i=1,nel
          IF(off(i) ==zero)THEN
            det(i)=one
           ELSEIF((det(i)<=volmin).OR.(det(i)<=zero))THEN
            icor=1
          ENDIF
        ENDDO
        IF (icor>0) THEN
          DO i=1,nel
            IF((off(i)/=zero).AND.
     .         (det(i)<=volmin.OR.det(i)<=zero))THEN
              det(i)=one
              off(i)=zero
#include "lockon.inc"
                WRITE(istdo,2000) ngl(i)
                WRITE(iout ,2000) ngl(i)
#include "lockoff.inc"
              idel7nok = 1
            ENDIF
          ENDDO
        ENDIF
      ELSEIF (ismstr /=4 ) THEN
        icor = 0
        DO i=1,nel
          IF(off(i) ==zero)THEN
            det(i)=one
          ELSEIF((det(i)<=volmin).OR.(det(i)<=zero))THEN
            icor = 1
          ENDIF
        ENDDO
        IF (icor>0) THEN
        DO i=1,nel
          IF(off(i) == zero)THEN
            det(i)=one
          ELSEIF(offg(i) > one)THEN
            
          ELSEIF((det(i)<=volmin).OR.(det(i)<=zero))THEN
            nnega=nnega+1
              index(nnega)=i
#include "lockon.inc"
            WRITE(istdo,3000) ngl(i)
            WRITE(iout ,3000) ngl(i)
#include "lockoff.inc"
          ENDIF
        ENDDO
          IF (ineg_v==0) THEN
           CALL ancmsg(msgid=280,anmode=aninfo)
           mstop = 1
          ENDIF
        END IF !(ICOR>0) THEN
      ELSE
C
        icor = 0
        DO i=1,nel
          IF(off(i) ==zero)THEN
            det(i)=one
          ELSEIF(det(i)<=zero)THEN
            icor=1
          ENDIF
        ENDDO
        IF (icor>0) THEN
          DO i=1,nel
           IF(off(i)/=zero)THEN
            IF(det(i)<=zero)THEN
              CALL ancmsg(msgid=166,anmode=aninfo,
     .                    i1=ngl(i))
              mstop = 1
            ENDIF
           ENDIF
          ENDDO
        ENDIF
      ENDIF
C
C Projection is not changed---
      IF (nnega>0) THEN
#include "vectorize.inc"
         DO j=1,nnega
          i = index(j)
            x1(i)=sav(i,1)
            y1(i)=sav(i,2)
            z1(i)=sav(i,3)
            x2(i)=sav(i,4)
            y2(i)=sav(i,5)
            z2(i)=sav(i,6)
            x3(i)=sav(i,7)
            y3(i)=sav(i,8)
            z3(i)=sav(i,9)
            x4(i)=sav(i,10)
            y4(i)=sav(i,11)
            z4(i)=sav(i,12)
            x5(i)=sav(i,13)
            y5(i)=sav(i,14)
            z5(i)=sav(i,15)
            x6(i)=zero
            y6(i)=zero
            z6(i)=zero
C
           x21(i)=x2(i)-x1(i)
           x31(i)=x3(i)-x1(i)
           x41(i)=x4(i)-x1(i)
           x54(i)=x5(i)-x4(i)
           x64(i)=x6(i)-x4(i)
           y21(i)=y2(i)-y1(i)
           y31(i)=y3(i)-y1(i)
           y41(i)=y4(i)-y1(i)
           y54(i)=y5(i)-y4(i)
           y64(i)=y6(i)-y4(i)
           z21(i)=z2(i)-z1(i)
           z31(i)=z3(i)-z1(i)
           z41(i)=z4(i)-z1(i)
           z54(i)=z5(i)-z4(i)
           z64(i)=z6(i)-z4(i)
C
           jac1(i)=x21(i)+x54(i)
           jac2(i)=y21(i)+y54(i)
           jac3(i)=z21(i)+z54(i)
C----
           jac4(i)=x31(i)+x64(i)
           jac5(i)=y31(i)+y64(i)
           jac6(i)=z31(i)+z64(i)
           jac7(i)=third*(x41(i)+x5(i)-x2(i)+x6(i)-x3(i))
           jac8(i)=third*(y41(i)+y5(i)-y2(i)+y6(i)-y3(i))
           jac9(i)=third*(z41(i)+z5(i)-z2(i)+z6(i)-z3(i))
C
           jacob4(i)=jac4(i)
           jacob5(i)=jac5(i)
           jacob6(i)=jac6(i)
           jacob7(i)=jac7(i)
           jacob8(i)=jac8(i)
           jacob9(i)=jac9(i)
C
           jac_59_68(i)=jac5(i)*jac9(i)-jac6(i)*jac8(i)
           jac_67_49(i)=jac6(i)*jac7(i)-jac4(i)*jac9(i)
           jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
C
           det(i)=one_over_8*(jac1(i)*jac_59_68(i)+jac2(i)*jac_67_49(i)+jac3(i)*jac_48_57(i))
           offg(i) = two
         ENDDO
      END IF
C
C Jacobian matrix inverse
      DO i=1,nel
        dett(i)=one_over_8/det(i)
C
        jaci1=dett(i)*jac_59_68(i)
        jaci4=dett(i)*jac_67_49(i)
        jaci7=dett(i)*jac_48_57(i)
        jaci2=dett(i)*(-jac2(i)*jac9(i)+jac3(i)*jac8(i))
        jaci5=dett(i)*( jac1(i)*jac9(i)-jac3(i)*jac7(i))
        jaci8=dett(i)*(-jac1(i)*jac8(i)+jac2(i)*jac7(i))
        jaci3=dett(i)*( jac2(i)*jac6(i)-jac3(i)*jac5(i))
        jaci6=dett(i)*(-jac1(i)*jac6(i)+jac3(i)*jac4(i))
        jaci9=dett(i)*( jac1(i)*jac5(i)-jac2(i)*jac4(i))
C
        jaci12=jaci1+jaci2
        jaci45=jaci4+jaci5
        jaci78=jaci7+jaci8
C
C Symmetry (a b c a b c)->P1-P3, anti-symmetry(-1 -1 -1 1 1 1)->P4
        px1(i)=-jaci12
        py1(i)=-jaci45
        pz1(i)=-jaci78
        px2(i)=jaci1
        py2(i)=jaci4
        pz2(i)=jaci7
        px3(i)=jaci2
        py3(i)=jaci5
        pz3(i)=jaci8
        px4(i)=third*jaci3
        py4(i)=third*jaci6
        pz4(i)=third*jaci9
        jaci33(i) = jaci9*one_over_12
      ENDDO
C
C For shear traitement----------
      DO i=1,nel
       fac = dett(i)*one_over_12
       b1x(i,1)=-fac*jac1(i)*jac2(i)
       b1x(i,2)=-fac*jac4(i)*jac5(i)
       b1y(i,1)=-fac*jac2(i)*jac2(i)
       b1y(i,2)=-fac*jac5(i)*jac5(i)
       b2x(i,1)=-fac*jac1(i)*jac1(i)
       b2x(i,2)=-fac*jac4(i)*jac4(i)
       b2y(i,1)=b1x(i,1)
       b2y(i,2)=b1x(i,2)
       fac = fac*2.0
       b1122(i)=fac*jac1(i)*jac5(i)
       b1221(i)=fac*jac2(i)*jac4(i)
       b2212(i)=fac*jac5(i)*jac2(i)
       b1121(i)=b2212(i)
C
       b1xh(i,1)=-fac*(x54(i)*y54(i)-x21(i)*y21(i))
       b1xh(i,2)=-fac*(x64(i)*y64(i)-x31(i)*y31(i))
       b1yh(i,1)=-fac*(y54(i)*y54(i)-y21(i)*y21(i))
       b1yh(i,2)=-fac*(y64(i)*y64(i)-y31(i)*y31(i))
       b2xh(i,1)=-fac*(x54(i)*x54(i)-x21(i)*x21(i))
       b2xh(i,2)=-fac*(x64(i)*x64(i)-x31(i)*x31(i))
       b2yh(i,1)=b1xh(i,1)
       b2yh(i,2)=b1xh(i,2)
       fac = fac*two
       b1122h(i)=fac*(x54(i)*y64(i)-x21(i)*y31(i))
       b1221h(i)=fac*(x64(i)*y54(i)-x31(i)*y21(i))
       b2212h(i)=fac*(y54(i)*y64(i)-y21(i)*y31(i))
       b1121h(i)=b2212h(i)
      ENDDO
C Non constant part
      DO i=1,nel
       jac1(i)=-x21(i)+x54(i)
       jac2(i)=-y21(i)+y54(i)
       jac3(i)=-z21(i)+z54(i)
       jac4(i)=-x31(i)+x64(i)
       jac5(i)=-y31(i)+y64(i)
       jac6(i)=-z31(i)+z64(i)
      ENDDO
C
      DO i=1,nel
       jac_59_68(i)=jac5(i)*jac9(i)-jac6(i)*jac8(i)
       jac_67_49(i)=jac6(i)*jac7(i)-jac4(i)*jac9(i)
       jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
      ENDDO
C
      DO i=1,nel
       jaci1=dett(i)*jac_59_68(i)
       jaci4=dett(i)*jac_67_49(i)
       jaci7=dett(i)*jac_48_57(i)
       jaci2=dett(i)*(-jac2(i)*jac9(i)+jac3(i)*jac8(i))
       jaci5=dett(i)*( jac1(i)*jac9(i)-jac3(i)*jac7(i))
       jaci8=dett(i)*(-jac1(i)*jac8(i)+jac2(i)*jac7(i))
C
       jaci12=jaci1+jaci2
       jaci45=jaci4+jaci5
       jaci78=jaci7+jaci8
C
C Symmetry(a b c a b c)->P1-P3
       px1h(i)=-jaci12
       py1h(i)=-jaci45
       pz1h(i)=-jaci78
       px2h(i)=jaci1
       py2h(i)=jaci4
       pz2h(i)=jaci7
       px3h(i)=jaci2
       py3h(i)=jaci5
       pz3h(i)=jaci8
      ENDDO
C                                                                     12
       DO i=1,nel
        vzl(i) = fourth*(jacob9(i)*(
     .          x54(i)*y64(i)-x21(i)*y31(i)-x64(i)*y54(i)+x31(i)*y21(i))
     .                  -jacob8(i)*(
     .          x54(i)*z64(i)+x31(i)*z21(i)-x21(i)*z31(i)-x64(i)*z54(i))
     .                  +jacob7(i)*(
     .          y54(i)*z64(i)+y31(i)*z21(i)-y21(i)*z31(i)-y64(i)*z54(i))
     .                )
        volg(i)= det(i)
       ENDDO
      RETURN
C
 1000 FORMAT(/' ZERO OR NEGATIVE VOLUME : 3D-ELEMENT NB',i10/)
 2000 FORMAT(/' ZERO OR NEGATIVE VOLUME : DELETE 3D-ELEMENT NB',i10/)
 3000 FORMAT(/' ZERO OR NEGATIVE VOLUME : 3D-ELEMENT NB:',i10/,
     +    ' SOLID-SHELL ELEMENT IS SWITCHED TO SMALL STRAIN OPTION'/) 
      END