scderi3.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/.
!||====================================================================
!||    scderi3       ../engine/source/elements/thickshell/solidec/scderi3.F
!||--- called by ------------------------------------------------------
!||    scforc3       ../engine/source/elements/thickshell/solidec/scforc3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg        ../engine/source/output/message/message.F
!||--- uses       -----------------------------------------------------
!||    message_mod   ../engine/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE scderi3(
     1   OFF,     DET,     NGL,     X1,
     2   X2,      X3,      X4,      X5,
     3   X6,      X7,      X8,      Y1,
     4   Y2,      Y3,      Y4,      Y5,
     5   Y6,      Y7,      Y8,      Z1,
     6   Z2,      Z3,      Z4,      Z5,
     7   Z6,      Z7,      Z8,      PX1,
     8   PX2,     PX3,     PX4,     PY1,
     9   PY2,     PY3,     PY4,     PZ1,
     A   PZ2,     PZ3,     PZ4,     PX1H1,
     B   PX1H2,   PX1H3,   PX1H4,   PX2H1,
     C   PX2H2,   PX2H3,   PX2H4,   PX3H1,
     D   PX3H2,   PX3H3,   PX3H4,   PX4H1,
     E   PX4H2,   PX4H3,   PX4H4,   JAC1,
     F   JAC2,    JAC3,    JAC4,    JAC5,
     G   JAC6,    RX0,     RY0,     SX0,
     H   SY0,     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
      DOUBLE PRECISION
     .   X1(*), X2(*), X3(*), X4(*), X5(*), X6(*), X7(*), X8(*),
     .   Y1(*), Y2(*), Y3(*), Y4(*), Y5(*), Y6(*), Y7(*), Y8(*),  
     .   Z1(*), Z2(*), Z3(*), Z4(*), Z5(*), Z6(*), Z7(*), Z8(*),
     .   SAV(NEL,21)
     
      my_real
     .   OFF(*),DET(*),
     .   PX1(*), PX2(*), PX3(*), PX4(*),  
     .   PY1(*), PY2(*), PY3(*), PY4(*),  
     .   PZ1(*), PZ2(*), PZ3(*), PZ4(*),  
     .   PX1H1(*), PX1H2(*), PX1H3(*),PX1H4(*),  
     .   PX2H1(*), PX2H2(*), PX2H3(*),PX2H4(*),  
     .   PX3H1(*), PX3H2(*), PX3H3(*),PX3H4(*),  
     .   PX4H1(*), PX4H2(*), PX4H3(*),PX4H4(*),  
     .   jac1(*),jac2(*),jac3(*),
     .   jac4(*),jac5(*),jac6(*),
     .   rx0(*),ry0(*),sx0(*),sy0(*),vzl(*),volg(*),offg(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I,J,NNEGA,INDEX(MVSIZ)
C                                                                     12
      my_real
     .   dett ,  jac7(mvsiz), jac8(mvsiz) ,jac9(mvsiz), 
     .   jaci1, jaci2, jaci3,
     .   jaci4, jaci5, jaci6,
     .   jaci7, jaci8, jaci9,
     .   x17(mvsiz) , x28(mvsiz) , x35(mvsiz) , x46(mvsiz),
     .   y17(mvsiz) , y28(mvsiz) , y35(mvsiz) , y46(mvsiz),
     .   z17(mvsiz) , z28(mvsiz) , z35(mvsiz) , z46(mvsiz),
     .   jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz),jac_19_37(mvsiz),
     .   jaci12, jaci45, jaci78,
     .   x_17_46 , x_28_35 ,
     .   y_17_46 , y_28_35 ,
     .   z_17_46 , z_28_35 ,area(mvsiz) ,
     .   hx,hy,hz,a_i,
     .   h1x(mvsiz),h1y(mvsiz),h2x(mvsiz),h2y(mvsiz),
     .   h1z(mvsiz),h2z(mvsiz),icor
C=======================================================================
C
      nnega=0
      DO i=1,nel
      x17(i)=x7(i)-x1(i)
      x28(i)=x8(i)-x2(i)
      x35(i)=x5(i)-x3(i)
      x46(i)=x6(i)-x4(i)
      y17(i)=y7(i)-y1(i)
      y28(i)=y8(i)-y2(i)
      y35(i)=y5(i)-y3(i)
      y46(i)=y6(i)-y4(i)
      z17(i)=z7(i)-z1(i)
      z28(i)=z8(i)-z2(i)
      z35(i)=z5(i)-z3(i)
      z46(i)=z6(i)-z4(i)
      END DO
C
C Jacobian matrix
      DO i=1,nel
C  -------ri.xi-----------
        jac3(i)=x17(i)+x28(i)-x35(i)-x46(i)
        jac1(i)=y17(i)+y28(i)-y35(i)-y46(i)
        jac2(i)=z17(i)+z28(i)-z35(i)-z46(i)
C-------
        x_17_46=x17(i)+x46(i)
        x_28_35=x28(i)+x35(i)
        y_17_46=y17(i)+y46(i)
        y_28_35=y28(i)+y35(i)
        z_17_46=z17(i)+z46(i)
        z_28_35=z28(i)+z35(i)
C  -------si.xi-----------
        jac6(i)=x_17_46+x_28_35
        jac4(i)=y_17_46+y_28_35
        jac5(i)=z_17_46+z_28_35
C  -------ti.xi-----------
        jac9(i)=x_17_46-x_28_35
        jac7(i)=y_17_46-y_28_35
        jac8(i)=z_17_46-z_28_35
      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_19_37(i)=jac1(i)*jac9(i)-jac3(i)*jac7(i)
        jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
      ENDDO
C
      DO i=1,nel
       area(i) = one_over_64*jac_19_37(i)
       det(i)=one_over_64*(jac1(i)*jac_59_68(i)+jac2(i)*jac_67_49(i)+jac3(i)*jac_48_57(i))
      ENDDO
C
      IF(idtmin(1)==1)THEN
        icor = 0
        DO i=1,nel
          IF(off(i) ==zero)THEN
            area(i) = one
            det(i)=one
          ELSEIF((det(i)<=volmin).OR.(det(i)<=zero)
     .           .OR.(area(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
              area(i) = one
              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.OR.(area(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
            area(i) = one
            det(i)=one
          ELSEIF((det(i)<=volmin).OR.(det(i)<=zero)
     .           .OR.(area(i)<=zero))THEN
            icor=1
          ENDIF
        ENDDO
        IF (icor>0) THEN
          DO i=1,nel
            IF((off(i)/=0.).AND.
     .         (det(i)<=volmin.OR.det(i)<=zero)
     .           .OR.(area(i)<=zero))THEN
              area(i) = one
              det(i)=one
              off(i)=zero
#include "lockon.inc"
                WRITE(istdo,2000) ngl(i)
                WRITE(iout ,2000) ngl(i)
#include "lockoff.inc"
            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.(area(i)<=zero))THEN
            icor = 1
          ENDIF
        ENDDO
        IF (icor>0) THEN
        DO i=1,nel
          IF(off(i) == zero)THEN
            det(i)=one
            area(i) = one
          ELSEIF(offg(i) > one)THEN
            
          ELSEIF((det(i)<=volmin).OR.(area(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
        icor = 0
        DO i=1,nel
          IF(off(i) ==zero)THEN
            det(i)=one
            area(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      
      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)=sav(i,16)
            y6(i)=sav(i,17)
            z6(i)=sav(i,18)
            x7(i)=sav(i,19)
            y7(i)=sav(i,20)
            z7(i)=sav(i,21)
            x8(i)=zero
            y8(i)=zero
            z8(i)=zero
C
C---------
            jac3(i)=x17(i)+x28(i)-x35(i)-x46(i)
            jac1(i)=y17(i)+y28(i)-y35(i)-y46(i)
            jac2(i)=z17(i)+z28(i)-z35(i)-z46(i)
C-------
            x_17_46=x17(i)+x46(i)
            x_28_35=x28(i)+x35(i)
            y_17_46=y17(i)+y46(i)
            y_28_35=y28(i)+y35(i)
            z_17_46=z17(i)+z46(i)
            z_28_35=z28(i)+z35(i)
C----------
            jac6(i)=x_17_46+x_28_35
            jac4(i)=y_17_46+y_28_35
            jac5(i)=z_17_46+z_28_35
            jac9(i)=x_17_46-x_28_35
            jac7(i)=y_17_46-y_28_35
            jac8(i)=z_17_46-z_28_35
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_19_37(i)=jac1(i)*jac9(i)-jac3(i)*jac7(i)
            jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
C
            area(i) = one_over_64*jac_19_37(i)
            det(i)=one_over_64*(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=one_over_64/det(i)
        jaci1=dett*jac_59_68(i)
        jaci4=dett*jac_67_49(i)
        jaci7=dett*jac_48_57(i)
        jaci2=dett*(-jac2(i)*jac9(i)+jac3(i)*jac8(i))
        jaci5=dett*jac_19_37(i)
        jaci8=dett*(-jac1(i)*jac8(i)+jac2(i)*jac7(i))
        jaci3=dett*( jac2(i)*jac6(i)-jac3(i)*jac5(i))
        jaci6=dett*(-jac1(i)*jac6(i)+jac3(i)*jac4(i))
        jaci9=dett*( jac1(i)*jac5(i)-jac2(i)*jac4(i))
C
        jaci12=jaci1-jaci2
        jaci45=jaci4-jaci5
        jaci78=jaci7-jaci8
C
        py3(i)= jaci12+jaci3
        pz3(i)= jaci45+jaci6
        px3(i)= jaci78+jaci9
        py4(i)= jaci12-jaci3
        pz4(i)= jaci45-jaci6
        px4(i)= jaci78-jaci9
 
        jaci12=jaci1+jaci2
        jaci45=jaci4+jaci5
        jaci78=jaci7+jaci8
 
        py1(i)=-jaci12-jaci3
        pz1(i)=-jaci45-jaci6
        px1(i)=-jaci78-jaci9
        py2(i)=-jaci12+jaci3
        pz2(i)=-jaci45+jaci6
        px2(i)=-jaci78+jaci9
      ENDDO
C
       DO i=1,nel
C   h3
C 1 -1 1 -1 1 -1 1 -1
         hx=one_over_8*(x1(i)-x2(i)+x3(i)-x4(i)+x5(i)-x6(i)+x7(i)-x8(i))
         hy=one_over_8*(y1(i)-y2(i)+y3(i)-y4(i)+y5(i)-y6(i)+y7(i)-y8(i))
         hz=one_over_8*(z1(i)-z2(i)+z3(i)-z4(i)+z5(i)-z6(i)+z7(i)-z8(i))
         px1h3(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
         px2h3(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
         px3h3(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
         px4h3(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
       END DO
C   h1
C 1 1 -1 -1 -1 -1 1 1
       DO i=1,nel
         h1x(i)=x1(i)+x2(i)-x3(i)-x4(i)-x5(i)-x6(i)+x7(i)+x8(i)
         h1y(i)=y1(i)+y2(i)-y3(i)-y4(i)-y5(i)-y6(i)+y7(i)+y8(i)
         h1z(i)=z1(i)+z2(i)-z3(i)-z4(i)-z5(i)-z6(i)+z7(i)+z8(i)
         hx=one_over_8*h1x(i)
         hy=one_over_8*h1y(i)
         hz=one_over_8*h1z(i)
         px1h1(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
         px2h1(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
         px3h1(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
         px4h1(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
       END DO
C   h2
C 1 -1 -1 1 -1 1 1 -1
       DO i=1,nel
         h2x(i)=x1(i)-x2(i)-x3(i)+x4(i)-x5(i)+x6(i)+x7(i)-x8(i)
         h2y(i)=y1(i)-y2(i)-y3(i)+y4(i)-y5(i)+y6(i)+y7(i)-y8(i)
         h2z(i)=z1(i)-z2(i)-z3(i)+z4(i)-z5(i)+z6(i)+z7(i)-z8(i)
         hx=one_over_8*h2x(i)
         hy=one_over_8*h2y(i)
         hz=one_over_8*h2z(i)   
         px1h2(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
         px2h2(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
         px3h2(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
         px4h2(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
       END DO
C   h4
C -1 1 -1 1 1 -1 1 -1
       DO i=1,nel
         hx=one_over_8*(-x1(i)+x2(i)-x3(i)+x4(i)+x5(i)-x6(i)+x7(i)-x8(i))
         hy=one_over_8*(-y1(i)+y2(i)-y3(i)+y4(i)+y5(i)-y6(i)+y7(i)-y8(i))
         hz=one_over_8*(-z1(i)+z2(i)-z3(i)+z4(i)+z5(i)-z6(i)+z7(i)-z8(i))   
         px1h4(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
         px2h4(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
         px3h4(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
         px4h4(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
       END DO
C
       DO i=1,nel
        a_i = one_over_8/area(i)
        rx0(i) = jac9(i)*a_i
        ry0(i) = jac7(i)*a_i
        sx0(i) = jac3(i)*a_i
        sy0(i) = jac1(i)*a_i
        vzl(i) = one_over_64*jac5(i)*(
     .          jac9(i)*h1y(i)+jac1(i)*h2x(i)-jac3(i)*h2y(i)-jac7(i)*h1x(i))
     .         + one_over_64*jac4(i)*(
     .          jac3(i)*h2z(i)+jac8(i)*h1x(i)-jac9(i)*h1z(i)-jac2(i)*h2x(i))
     .         + one_over_64*jac6(i)*(
     .          jac7(i)*h1z(i)+jac2(i)*h2y(i)-jac1(i)*h2z(i)-jac8(i)*h1y(i))
C
        volg(i)= det(i)
       ENDDO
C-----------
      RETURN
 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'/) 
C-----------
      END