qdefo2.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/.
!||====================================================================
!||    qdefo2   ../engine/source/elements/solid_2d/quad/qdefo2.F
!||--- called by ------------------------------------------------------
!||    bforc2   ../engine/source/ale/bimat/bforc2.F
!||    qforc2   ../engine/source/elements/solid_2d/quad/qforc2.F
!||====================================================================
      SUBROUTINE qdefo2(
     1   V,       W,       Y1,      Y2,
     2   Y3,      Y4,      Z1,      Z2,
     3   Z3,      Z4,      VY1,     VY2,
     4   VY3,     VY4,     VZ1,     VZ2,
     5   VZ3,     VZ4,     PY1,     PY2,
     6   PZ1,     PZ2,     WYZ,     DYZ,
     7   DZY,     EYY,     EZZ,     ETT,
     8   EYZ,     EYT,     EZT,     RX,
     9   RY,      RZ,      SX,      SY,
     A   SZ,      TX,      TY,      TZ,
     B   VOLN,    AIRE,    AIREM,   NC1,
     C   NC2,     NC3,     NC4,     NEL)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: NEL
      my_real
     .   V(3,*), W(3,*),
     .   Y1(*),Y2(*),Y3(*),Y4(*),Z1(*),Z2(*),Z3(*),Z4(*),
     .   VY1(*),VY2(*),VY3(*),VY4(*),VZ1(*),VZ2(*),VZ3(*),VZ4(*),
     .   PY1(*), PY2(*), PZ1(*), PZ2(*),
     .   WYZ(*), DYZ(*), DZY(*),
     .   EYY(*),EZZ(*), ETT(*), EYZ(*), EYT(*), EZT(*),
     .   RX(*),RY(*), RZ(*), SX(*), SY(*), SZ(*), TX(*), TY(*), TZ(*),
     .   VOLN(*),AIRE(*),AIREM(*)
C
      INTEGER NC1(*), NC2(*), NC3(*), NC4(*)
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "com08_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .   VY13(MVSIZ), VY24(MVSIZ), VZ13(MVSIZ),VZ24(MVSIZ),
     .   ym1(mvsiz), ym2(mvsiz), ym3(mvsiz), ym4(mvsiz),
     .   zm1(mvsiz), zm2(mvsiz), zm3(mvsiz), zm4(mvsiz), dth,
     .   yavg(mvsiz),a1(mvsiz) , a2(mvsiz)
C-----------------------------------------------
      dth=half*dt1
      
      DO i=1,nel
        !from qforc2 subroutine : &W  = &V  
        !from bforc2 subroutine : &W /= &V  
        vy1(i)=w(2,nc1(i))
        vz1(i)=w(3,nc1(i))
        vy2(i)=w(2,nc2(i))
        vz2(i)=w(3,nc2(i))
        vy3(i)=w(2,nc3(i))
        vz3(i)=w(3,nc3(i))
        vy4(i)=w(2,nc4(i))
        vz4(i)=w(3,nc4(i))
      ENDDO
 
      DO i=1,nel
        ym1(i)=y1(i)-dth*vy1(i)
        ym2(i)=y2(i)-dth*vy2(i)
        ym3(i)=y3(i)-dth*vy3(i)
        ym4(i)=y4(i)-dth*vy4(i)
        zm1(i)=z1(i)-dth*vz1(i)
        zm2(i)=z2(i)-dth*vz2(i)
        zm3(i)=z3(i)-dth*vz3(i)
        zm4(i)=z4(i)-dth*vz4(i)
      ENDDO
 
      DO i=1,nel
        vy1(i)=v(2,nc1(i))
        vz1(i)=v(3,nc1(i))
        vy2(i)=v(2,nc2(i))
        vz2(i)=v(3,nc2(i))
        vy3(i)=v(2,nc3(i))
        vz3(i)=v(3,nc3(i))
        vy4(i)=v(2,nc4(i))
        vz4(i)=v(3,nc4(i))
      ENDDO
 
      DO i=1,nel
        py1(i)=half*(zm2(i)-zm4(i))
        py2(i)=half*(zm3(i)-zm1(i))
        pz1(i)=half*(ym4(i)-ym2(i))
        pz2(i)=half*(ym1(i)-ym3(i))      
      ENDDO
 
      DO i=1,nel
        a1(i) =ym2(i)*(zm3(i)-zm4(i))+ym3(i)*(zm4(i)-zm2(i))+ym4(i)*(zm2(i)-zm3(i))
        a2(i) =ym2(i)*(zm4(i)-zm1(i))+ym4(i)*(zm1(i)-zm2(i))+ym1(i)*(zm2(i)-zm4(i))
        airem(i)=(a1(i)+a2(i))*half
        !AUTRE FORMULE AIREM(I)=2.*(PZ2(I)*PY1(I)-PZ1(I)*PY2(I))
        yavg(i) =ym1(i)+ym2(i)+ym3(i)+ym4(i)
      ENDDO
 
      DO i=1,nel
        vy13(i)=vy1(i)-vy3(i)
        vy24(i)=vy2(i)-vy4(i)
        vz13(i)=vz1(i)-vz3(i)
        vz24(i)=vz2(i)-vz4(i)
       ENDDO
 
      DO i=1,nel
        IF(airem(i)>zero) THEN
          eyy(i)=(py1(i)*vy13(i)+py2(i)*vy24(i))/airem(i)
          ezz(i)=(pz1(i)*vz13(i)+pz2(i)*vz24(i))/airem(i)
          ett(i)=zero
          dzy(i)=(py1(i)*vz13(i)+py2(i)*vz24(i))/airem(i)
          dyz(i)=(pz1(i)*vy13(i)+pz2(i)*vy24(i))/airem(i)
          eyt(i)=zero
          ezt(i)=zero   
        ELSE
          eyy(i)=zero
          ezz(i)=zero
          ett(i)=zero
          dzy(i)=zero
          dyz(i)=zero
          eyt(i)=zero
          ezt(i)=zero
        END IF
      ENDDO
 
      IF(n2d==1) THEN
        DO i=1,nel
          IF (yavg(i)<=zero) cycle
          ett(i)=(vy1(i)+vy2(i)+vy3(i)+vy4(i))/yavg(i)
        ENDDO
      ENDIF
 
      DO i=1,nel
        eyz(i)= dzy(i)+dyz(i)
        wyz(i)=half*dt1*(dzy(i)-dyz(i))
      ENDDO
C---------------------------------------
C     DERIVEES POUR LE CALCUL DES FORCES
C---------------------------------------
      DO i=1,nel
        py1(i)=half*(z2(i)-z4(i))
        py2(i)=half*(z3(i)-z1(i))
        pz1(i)=half*(y4(i)-y2(i))
        pz2(i)=half*(y1(i)-y3(i))      
      ENDDO
C---------------------------------------
C     REPERE LOCAL POUR MATERIAUX NON ISOTROPES
C---------------------------------------
      DO i=1,nel
        rx(i)=one
        ry(i)=zero
        rz(i)=zero
        sx(i)=zero    
        sy(i)=-pz1(i)-pz2(i)
        sz(i)=+py1(i)+py2(i)
        tx(i)=zero
        ty(i)=+pz1(i)-pz2(i)
        tz(i)=-py1(i)+py2(i)
      ENDDO
 
      RETURN
      END