s4derit3.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/.
!||====================================================================
!||    s4derit3   ../engine/source/elements/solid/solide4/s4derit3.F
!||--- called by ------------------------------------------------------
!||    s4forc3    ../engine/source/elements/solid/solide4/s4forc3.F
!||    s4ke3      ../engine/source/elements/solid/solide4/s4ke3.f
!||--- calls      -----------------------------------------------------
!||    schkjab3   ../engine/source/elements/solid/solide/schkjab3.F
!||====================================================================
      SUBROUTINE s4derit3(
     1   OFF,     DET,     NGL,     DELTAX,
     2   MXT,     X1,      X2,      X3,
     3   X4,      Y1,      Y2,      Y3,
     4   Y4,      Z1,      Z2,      Z3,
     5   Z4,      PX1,     PX2,     PX3,
     6   PX4,     PY1,     PY2,     PY3,
     7   PY4,     PZ1,     PZ2,     PZ3,
     8   PZ4,     RX,      RY,      RZ,
     9   SX,      SY,      SZ,      TX,
     A   TY,      TZ,      PM,      VOLDP,
     B   NEL,     IFORMDT)
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"
#include      "param_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: NEL
      INTEGER, INTENT(IN) :: IFORMDT
      INTEGER MXT(*)
      DOUBLE PRECISION
     .   X1(*), X2(*), X3(*), X4(*),
     .   Y1(*), Y2(*), Y3(*), Y4(*),
     .   Z1(*), Z2(*), Z3(*), Z4(*),VOLDP(*)
      my_real
     .   OFF(*),DET(*),DELTAX(*),
     .   PX1(*), PX2(*), PX3(*), PX4(*),  
     .   PY1(*), PY2(*), PY3(*), PY4(*),  
     .   pz1(*), pz2(*), pz3(*), pz4(*), 
     .   rx(*), ry(*), rz(*), sx(*), sy(*), sz(*),tx(*), ty(*), tz(*),
     .   pm(npropm,*)
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "scr17_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NGL(*), I
      DOUBLE PRECISION
     .   X41, Y41, Z41, X42, Y42, Z42, X43, Y43, Z43,B1DP,C1DP,D1DP
      my_real
     .   B1(MVSIZ), B2(MVSIZ), B3(MVSIZ), B4(MVSIZ), 
     .   C1(MVSIZ), C2(MVSIZ), C3(MVSIZ), C4(MVSIZ),
     .   D1(MVSIZ), D2(MVSIZ), D3(MVSIZ), D4(MVSIZ),
     .   A1, A2, A3, A4, D, AREAMX2,
     .   PXX, PYY, PZZ, PXY, PYZ, PXZ, GFAC, AA, BB, P, LD
C-----------------------------------------------
C
      DO i=1,nel
       x43 = x4(i) - x3(i)
       y43 = y4(i) - y3(i)
       z43 = z4(i) - z3(i)
       x41 = x4(i) - x1(i)
       y41 = y4(i) - y1(i)
       z41 = z4(i) - z1(i)
       x42 = x4(i) - x2(i)
       y42 = y4(i) - y2(i)
       z42 = z4(i) - z2(i)
C
       rx(i) =  -x41
       ry(i) =  -y41
       rz(i) =  -z41
       sx(i) =  -x42
       sy(i) =  -y42
       sz(i) =  -z42
C
       tx(i) =  -x43
       ty(i) =  -y43
       tz(i) =  -z43
C
       b1dp  =  y43*z42 - y42*z43
       b1(i) =  b1dp
       b2(i) =  y41*z43 - y43*z41
       b3(i) =  y42*z41 - y41*z42
       b4(i) =  -(b1(i) + b2(i) + b3(i))
C
       c1dp  =  z43*x42 - z42*x43
       c1(i) =  c1dp
       c2(i) =  z41*x43 - z43*x41
       c3(i) =  z42*x41 - z41*x42
       c4(i) =  -(c1(i) + c2(i) + c3(i))
C
       d1dp  =  x43*y42 - x42*y43
       d1(i) =  d1dp
       d2(i) =  x41*y43 - x43*y41
       d3(i) =  x42*y41 - x41*y42
       d4(i) =  -(d1(i) + d2(i) + d3(i))
C
       voldp(i) = (x41*b1dp + y41*c1dp + z41*d1dp)*one_over_6
       det(i) = voldp(i)
C
      ENDDO
C
      CALL schkjab3(
     1   off,     det,     ngl,     nel)
C
      DO i=1,nel
        d = one/det(i)/six
        px1(i)=-b1(i)*d
        py1(i)=-c1(i)*d
        pz1(i)=-d1(i)*d
        px2(i)=-b2(i)*d
        py2(i)=-c2(i)*d
        pz2(i)=-d2(i)*d
        px3(i)=-b3(i)*d
        py3(i)=-c3(i)*d
        pz3(i)=-d3(i)*d
        px4(i)=-b4(i)*d
        py4(i)=-c4(i)*d
        pz4(i)=-d4(i)*d
      END DO
 
      IF(idt1sol==0)THEN
 
        DO i=1,nel
          d = max(px1(i)*px1(i)+py1(i)*py1(i)+pz1(i)*pz1(i),
     .      px2(i)*px2(i)+py2(i)*py2(i)+pz2(i)*pz2(i),
     .      px3(i)*px3(i)+py3(i)*py3(i)+pz3(i)*pz3(i),
     .      px4(i)*px4(i)+py4(i)*py4(i)+pz4(i)*pz4(i))
          deltax(i) = one / sqrt(d)
        END DO
 
 
       ELSEIF(iformdt==1)THEN
 
         gfac=pm(105,mxt(1))
         ld  =two*sqrt(max(one-gfac,zero))+one
         DO i=1,nel
           pxx=px1(i)*px1(i)+px2(i)*px2(i)+px3(i)*px3(i)+px4(i)*px4(i)
           pyy=py1(i)*py1(i)+py2(i)*py2(i)+py3(i)*py3(i)+py4(i)*py4(i)
           pzz=pz1(i)*pz1(i)+pz2(i)*pz2(i)+pz3(i)*pz3(i)+pz4(i)*pz4(i)
           pxy=px1(i)*py1(i)+px2(i)*py2(i)+px3(i)*py3(i)+px4(i)*py4(i)
           pxz=px1(i)*pz1(i)+px2(i)*pz2(i)+px3(i)*pz3(i)+px4(i)*pz4(i)
           pyz=py1(i)*pz1(i)+py2(i)*pz2(i)+py3(i)*pz3(i)+py4(i)*pz4(i)
C
           aa = -(pxx+pyy+pzz)
           bb =  (pxx*pyy+pxx*pzz+pyy*pzz-pxy**2-pxz**2-pyz**2) 
           p  = bb-third*aa*aa
           d=two*sqrt(third*max(-p,zero))-third*aa
C
             d=ld*d
C
           deltax(i) = one / sqrt(d)
         END DO
 
       ELSEIF(iformdt==2)THEN
 
         gfac=pm(105,mxt(1))
         DO i=1,nel
           pxx=px1(i)*px1(i)+px2(i)*px2(i)+px3(i)*px3(i)+px4(i)*px4(i)
           pyy=py1(i)*py1(i)+py2(i)*py2(i)+py3(i)*py3(i)+py4(i)*py4(i)
           pzz=pz1(i)*pz1(i)+pz2(i)*pz2(i)+pz3(i)*pz3(i)+pz4(i)*pz4(i)
           pxy=px1(i)*py1(i)+px2(i)*py2(i)+px3(i)*py3(i)+px4(i)*py4(i)
           pxz=px1(i)*pz1(i)+px2(i)*pz2(i)+px3(i)*pz3(i)+px4(i)*pz4(i)
           pyz=py1(i)*pz1(i)+py2(i)*pz2(i)+py3(i)*pz3(i)+py4(i)*pz4(i)
C
           aa = -(pxx+pyy+pzz)
           bb =  gfac*(pxx*pyy+pxx*pzz+pyy*pzz-pxy**2-pxz**2-pyz**2) 
           p  = bb-third*aa*aa
           d  = two*sqrt(third*max(-p,zero))-third*aa
C
           deltax(i) = one / sqrt(d)
         END DO
 
      END IF
C
      RETURN
C
 1000 FORMAT(/' ZERO OR NEGATIVE VOLUME : 3D-ELEMENT NB',i10/)
 2000 FORMAT(/' ZERO OR NEGATIVE VOLUME : DELETE 3D-ELEMENT NB',i10/)
      END