s4defo3.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/.
!||====================================================================
!||    s4defo3   ../engine/source/elements/solid/solide4/s4defo3.F
!||--- called by ------------------------------------------------------
!||    s4forc3   ../engine/source/elements/solid/solide4/s4forc3.F
!||====================================================================
      SUBROUTINE s4defo3(
     1   PX1,     PX2,     PX3,     PX4,
     2   PY1,     PY2,     PY3,     PY4,
     3   PZ1,     PZ2,     PZ3,     PZ4,
     4   VX1,     VX2,     VX3,     VX4,
     5   VY1,     VY2,     VY3,     VY4,
     6   VZ1,     VZ2,     VZ3,     VZ4,
     7   DXX,     DXY,     DXZ,     DYX,
     8   DYY,     DYZ,     DZX,     DZY,
     9   DZZ,     D4,      D5,      D6,
     A   WXX,     WYY,     WZZ,     NEL,
     B   ISMSTR)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.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) :: ISMSTR
      my_real
     .   VX1(*), VX2(*), VX3(*), VX4(*),
     .   VY1(*), VY2(*), VY3(*), VY4(*),
     .   VZ1(*), VZ2(*), VZ3(*), VZ4(*),
     .   PX1(*), PX2(*), PX3(*), PX4(*),
     .   PY1(*), PY2(*), PY3(*), PY4(*),
     .   PZ1(*), PZ2(*), PZ3(*), PZ4(*), 
     .   DXX(*), DXY(*), DXZ(*),
     .   DYX(*), DYY(*), DYZ(*),
     .   dzx(*), dzy(*), dzz(*), d4(*), d5(*), d6(*),
     .   wxx(*), wyy(*), wzz(*)
CMasParINCLUDE 'sdefo3.intmap.inc'
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com08_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER  I
C     REAL
      my_real
     .   DT1D2
      my_real 
     .   PXX2,PYY2,PZZ2,PXX2P,PYY2P,PZZ2P,AAA,BBB,CCC,
     .   EXX,EXY,EXZ,EYX,EYY,EYZ,EZX,EZY,EZZ
C-----------------------------------------------
      DO I=1,nel
        dxx(i)=px1(i)*vx1(i)+px2(i)*vx2(i)+
     .         px3(i)*vx3(i)+px4(i)*vx4(i)
        dyy(i)=py1(i)*vy1(i)+py2(i)*vy2(i)+
     .         py3(i)*vy3(i)+py4(i)*vy4(i)
        dzz(i)=pz1(i)*vz1(i)+pz2(i)*vz2(i)+
     .         pz3(i)*vz3(i)+pz4(i)*vz4(i)
        dxy(i)=py1(i)*vx1(i)+py2(i)*vx2(i)+
     .         py3(i)*vx3(i)+py4(i)*vx4(i)
        dxz(i)=pz1(i)*vx1(i)+pz2(i)*vx2(i)+
     .         pz3(i)*vx3(i)+pz4(i)*vx4(i)
        dyx(i)=px1(i)*vy1(i)+px2(i)*vy2(i)+
     .         px3(i)*vy3(i)+px4(i)*vy4(i)
        dyz(i)=pz1(i)*vy1(i)+pz2(i)*vy2(i)+
     .         pz3(i)*vy3(i)+pz4(i)*vy4(i)
        dzx(i)=px1(i)*vz1(i)+px2(i)*vz2(i)+
     .         px3(i)*vz3(i)+px4(i)*vz4(i)
        dzy(i)=py1(i)*vz1(i)+py2(i)*vz2(i)+
     .         py3(i)*vz3(i)+py4(i)*vz4(i)
      ENDDO
C
      dt1d2=half*dt1
C
      IF(ismstr==2.OR.ismstr==4)THEN
C
C       2nd order correction is needed only w/large strain 
C       since small strain is not objective anyway
        DO i=1,nel
C
         exx=dxx(i)
         eyy=dyy(i)
         ezz=dzz(i)
         exy=dxy(i)
         eyx=dyx(i)
         exz=dxz(i)
         ezx=dzx(i)
         eyz=dyz(i)
         ezy=dzy(i)
C
C        Euler-Almansi A=1/2 [I-Invers(B)], with Left Cauchy-Green tensor B = F.Transp(F) 
C                        Invers(B)=Invers(Transp(F)).Invers(F)
C
C        small perturbation  A ~ (I-Transp(gradU)).(I-gradU)
C                            A ~ 1/2 (gradU + Transp(gradU) - Transp(gradU).gradU)
C
         dxx(i)  = dxx(i)-dt1d2*(exx*exx+eyx*eyx+ezx*ezx)
         dyy(i)  = dyy(i)-dt1d2*(eyy*eyy+ezy*ezy+exy*exy)
         dzz(i)  = dzz(i)-dt1d2*(ezz*ezz+exz*exz+eyz*eyz)
         aaa = dt1d2*(exx*exy+eyx*eyy+ezx*ezy)
         dxy(i)  = dxy(i) -aaa
         dyx(i)  = dyx(i) -aaa
         d4(i)   = dxy(i)+dyx(i)
         bbb = dt1d2*(eyy*eyz+ezy*ezz+exy*exz)
         dyz(i)  = dyz(i) -bbb
         dzy(i)  = dzy(i) -bbb
         d5(i)   = dyz(i)+dzy(i)
         ccc = dt1d2*(ezz*ezx+exz*exx+eyz*eyx)
         dxz(i)  = dxz(i) -ccc
         dzx(i)  = dzx(i) -ccc
         d6(i)   = dxz(i)+dzx(i)
C
C        WZZ(I)=DT1D2*(EYX-EXY)
C        WYY(I)=DT1D2*(EXZ-EZX)
C        WXX(I)=DT1D2*(EZY-EYZ)
C
C        B = F Transp(F) = V^2, F = V.R 
C        => R = [Invers(Transp(F)).F]^1/2
C        small perturbation  => R ~ I + 1/2 ( gradU - Transp(gradU) - Transp(gradU).gradU )
C
C                [ 0     wzz    -wyy ]
C        and R = [ wzz     0    -wxx ]
C                [ wyy   wxx       0 ]
C
         wzz(i)=dt1*(half*(eyx-exy)-aaa)
         wyy(i)=dt1*(half*(exz-ezx)-bbb)
         wxx(i)=dt1*(half*(ezy-eyz)-ccc)
        ENDDO
      ELSE ! ISMSTR==1 or Total Strain Formulation (Approximated Strain Rate)
        DO i=1,nel
         d4(i)   = dxy(i)+dyx(i)
         d5(i)   = dyz(i)+dzy(i)
         d6(i)   = dxz(i)+dzx(i)
         wzz(i)=dt1d2*(dyx(i)-dxy(i))
         wyy(i)=dt1d2*(dxz(i)-dzx(i))
         wxx(i)=dt1d2*(dzy(i)-dyz(i))
        ENDDO
      END IF
C
      RETURN
C
      END