amomtn3.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/.
!||====================================================================
!||    amomtn3    ../engine/source/elements/solid/solide/amomtn3.F
!||--- called by ------------------------------------------------------
!||    amomt3     ../engine/source/elements/solid/solide/amomt3.f
!||--- calls      -----------------------------------------------------
!||    upwind_v   ../engine/source/elements/solid/solide/upwind_v.F
!||--- uses       -----------------------------------------------------
!||    ale_mod    ../common_source/modules/ale/ale_mod.f
!||====================================================================
      SUBROUTINE amomtn3(
     1   RHO,     DELTAX,  VIS,     VOL,
     2   A0,      A1,      A2,      A3,
     3   AJV,     R,       S,       T,
     4   VDX0,    VDY0,    VDZ0,    VX0,
     5   VY0,     VZ0,     VX1,     VY1,
     6   VZ1,     VX2,     VY2,     VZ2,
     7   VX3,     VY3,     VZ3,     F11,
     8   F21,     F31,     F12,     F22,
     9   F32,     F13,     F23,     F33,
     A   F14,     F24,     F34,     F15,
     B   F25,     F35,     F16,     F26,
     C   F36,     F17,     F27,     F37,
     D   F18,     F28,     F38,     IAD22,
     E   NC1,     NC2,     NC3,     NC4,
     F   NC5,     NC6,     NC7,     NC8,
     G   NALE,    NEL)
C-----------------------------------------------
C   M o d u l e  s
C-----------------------------------------------
      USE ale_mod
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   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "inter22.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: NEL
      my_real
     .   A0(MVSIZ,3),A1(MVSIZ,3),A2(MVSIZ,3),A3(MVSIZ,3),
     .   AJV(*),R(MVSIZ,3),S(MVSIZ,3),T(MVSIZ,3),
     .   VX1(*),VX2(*),VX3(*),VX0(*),VDX0(*),
     .   VY1(*),VY2(*),VY3(*),VY0(*),VDY0(*),
     .   VZ1(*),VZ2(*),VZ3(*),VZ0(*),VDZ0(*),
     .   F11(*),F21(*),F31(*),F12(*),F22(*),F32(*),
     .   F13(*),F23(*),F33(*),F14(*),F24(*),F34(*),
     .   F15(*),F25(*),F35(*),F16(*),F26(*),F36(*),
     .   F17(*),F27(*),F37(*),F18(*),F28(*),F38(*),
     .   RHO(*),DELTAX(*),VIS(*),VOL(*),IAD22(*)
      INTEGER :: NC1(*), NC2(*), NC3(*), NC4(*), 
     .     NC5(*), NC6(*), NC7(*), NC8(*), NALE(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      my_real
     .   gam(mvsiz),
     .   dxxp(mvsiz),dxyp(mvsiz),dxzp(mvsiz),
     .   dyxp(mvsiz),dyyp(mvsiz),dyzp(mvsiz),
     .   dzxp(mvsiz),dzyp(mvsiz),dzzp(mvsiz),
     .   up0(mvsiz),up1(mvsiz),up2(mvsiz),up3(mvsiz),
     .   f1p(mvsiz),f2p(mvsiz),f3p(mvsiz),
     .   f1(mvsiz),f2(mvsiz),f3(mvsiz),fac,
     .   off1(mvsiz), off2(mvsiz), off3(mvsiz), off4(mvsiz), 
     .   off5(mvsiz), off6(mvsiz), off7(mvsiz), off8(mvsiz)  
      INTEGER I
C-----------------------------------------------
      CALL upwind_v(
     1   rho,     vis,     vdx0,    vdy0,
     2   vdz0,    r,       s,       t,
     3   deltax,  gam,     nel)
      off1(1:mvsiz) = one
      off2(1:mvsiz) = one
      off3(1:mvsiz) = one
      off4(1:mvsiz) = one
      off5(1:mvsiz) = one
      off6(1:mvsiz) = one
      off7(1:mvsiz) = one
      off8(1:mvsiz) = one
      DO i=1,nel
         IF (nale(nc1(i)) == 0) off1(i) = zero
         IF (nale(nc2(i)) == 0) off2(i) = zero
         IF (nale(nc3(i)) == 0) off3(i) = zero
         IF (nale(nc4(i)) == 0) off4(i) = zero
         IF (nale(nc5(i)) == 0) off5(i) = zero
         IF (nale(nc6(i)) == 0) off6(i) = zero
         IF (nale(nc7(i)) == 0) off7(i) = zero
         IF (nale(nc8(i)) == 0) off8(i) = zero
      ENDDO
      IF (ale%UPWIND%UPWM==2) THEN
 
      !-------------------------------------------!
      ! MOMENTUM ADVECTION                        !
      ! TG : Taylor Galerkine                     !
      !-------------------------------------------!        
       DO i=1,nel
        gam(i)=gam(i)/ajv(i)
 
        dxxp(i)=a0(i,1)*vx0(i)+a1(i,1)*vx1(i)+a2(i,1)*vx2(i)+a3(i,1)*vx3(i)  
        dxyp(i)=a0(i,2)*vx0(i)+a1(i,2)*vx1(i)+a2(i,2)*vx2(i)+a3(i,2)*vx3(i)  
        dxzp(i)=a0(i,3)*vx0(i)+a1(i,3)*vx1(i)+a2(i,3)*vx2(i)+a3(i,3)*vx3(i)  
        dyxp(i)=a0(i,1)*vy0(i)+a1(i,1)*vy1(i)+a2(i,1)*vy2(i)+a3(i,1)*vy3(i)  
        dyyp(i)=a0(i,2)*vy0(i)+a1(i,2)*vy1(i)+a2(i,2)*vy2(i)+a3(i,2)*vy3(i)  
        dyzp(i)=a0(i,3)*vy0(i)+a1(i,3)*vy1(i)+a2(i,3)*vy2(i)+a3(i,3)*vy3(i)  
        dzxp(i)=a0(i,1)*vz0(i)+a1(i,1)*vz1(i)+a2(i,1)*vz2(i)+a3(i,1)*vz3(i)  
        dzyp(i)=a0(i,2)*vz0(i)+a1(i,2)*vz1(i)+a2(i,2)*vz2(i)+a3(i,2)*vz3(i)  
        dzzp(i)=a0(i,3)*vz0(i)+a1(i,3)*vz1(i)+a2(i,3)*vz2(i)+a3(i,3)*vz3(i)  
 
        fac=eight*gam(i)
        up0(i)=one +fac*(a0(i,1)*vdx0(i)+a0(i,2)*vdy0(i)+a0(i,3)*vdz0(i))
        up1(i)=one +fac*(a1(i,1)*vdx0(i)+a1(i,2)*vdy0(i)+a1(i,3)*vdz0(i))
        up2(i)=one +fac*(a2(i,1)*vdx0(i)+a2(i,2)*vdy0(i)+a2(i,3)*vdz0(i))
        up3(i)=one +fac*(a3(i,1)*vdx0(i)+a3(i,2)*vdy0(i)+a3(i,3)*vdz0(i))
 
        fac=rho(i)*vol(i)/ajv(i)/sixty4
        f1p(i)= fac*(vdx0(i)*dxxp(i)+vdy0(i)*dxyp(i)+vdz0(i)*dxzp(i))
        f2p(i)= fac*(vdx0(i)*dyxp(i)+vdy0(i)*dyyp(i)+vdz0(i)*dyzp(i))
        f3p(i)= fac*(vdx0(i)*dzxp(i)+vdy0(i)*dzyp(i)+vdz0(i)*dzzp(i))
 
        f11(i) = f11(i) - up0(i)*f1p(i)*off1(i)
        f21(i) = f21(i) - up0(i)*f2p(i)*off1(i)
        f31(i) = f31(i) - up0(i)*f3p(i)*off1(i)
        f12(i) = f12(i) - up1(i)*f1p(i)*off2(i)
        f22(i) = f22(i) - up1(i)*f2p(i)*off2(i)
        f32(i) = f32(i) - up1(i)*f3p(i)*off2(i)
        f13(i) = f13(i) - up2(i)*f1p(i)*off3(i)
        f23(i) = f23(i) - up2(i)*f2p(i)*off3(i)
        f33(i) = f33(i) - up2(i)*f3p(i)*off3(i)
        f14(i) = f14(i) - up3(i)*f1p(i)*off4(i)
        f24(i) = f24(i) - up3(i)*f2p(i)*off4(i)
        f34(i) = f34(i) - up3(i)*f3p(i)*off4(i)
        f15(i) = f15(i) - f1p(i)*off5(i)
        f25(i) = f25(i) - f2p(i)*off5(i)
        f35(i) = f35(i) - f3p(i)*off5(i)
        f16(i) = f16(i) - f1p(i)*off6(i)
        f26(i) = f26(i) - f2p(i)*off6(i)
        f36(i) = f36(i) - f3p(i)*off6(i)
        f17(i) = f17(i) - f1p(i)*off7(i)
        f27(i) = f27(i) - f2p(i)*off7(i)
        f37(i) = f37(i) - f3p(i)*off7(i)
        f18(i) = f18(i) - f1p(i)*off8(i)
        f28(i) = f28(i) - f2p(i)*off8(i)
        f38(i) = f38(i) - f3p(i)*off8(i)
       ENDDO
       
       !inter22
       IF(int22>0)THEN
         DO i=1,nel
           IF(nint(iad22(i))==0)THEN
             !currently unplugged only for cut cell buffer      
             f1p(i) = zero
             f2p(i) = zero
             f3p(i) = zero
           ENDIF
         ENDDO
       ENDIF
 
       DO i=1,nel
         f1(i)=-gam(i)*(-f1p(i)*(dyyp(i)+dzzp(i))+f2p(i)*dxyp(i)+f3p(i)*dxzp(i))
         f2(i)=-gam(i)*( f1p(i)*dyxp(i)-f2p(i)*(dxxp(i)+dzzp(i))+f3p(i)*dyzp(i))
         f3(i)=-gam(i)*( f1p(i)*dzxp(i)+f2p(i)*dzyp(i)-f3p(i)*(dyyp(i)+dxxp(i)))
         f11(i) = f11(i) - f1(i)*off1(i)
         f21(i) = f21(i) - f2(i)*off1(i)
         f31(i) = f31(i) - f3(i)*off1(i)
         f12(i) = f12(i) - f1(i)*off2(i)
         f22(i) = f22(i) - f2(i)*off2(i)
         f32(i) = f32(i) - f3(i)*off2(i)
         f13(i) = f13(i) - f1(i)*off3(i)
         f23(i) = f23(i) - f2(i)*off3(i)
         f33(i) = f33(i) - f3(i)*off3(i)
         f14(i) = f14(i) - f1(i)*off4(i)
         f24(i) = f24(i) - f2(i)*off4(i)
         f34(i) = f34(i) - f3(i)*off4(i)
         f15(i) = f15(i) - f1(i)*off5(i)
         f25(i) = f25(i) - f2(i)*off5(i)
         f35(i) = f35(i) - f3(i)*off5(i)
         f16(i) = f16(i) - f1(i)*off6(i)
         f26(i) = f26(i) - f2(i)*off6(i)
         f36(i) = f36(i) - f3(i)*off6(i)
         f17(i) = f17(i) - f1(i)*off7(i)
         f27(i) = f27(i) - f2(i)*off7(i)
         f37(i) = f37(i) - f3(i)*off7(i)
         f18(i) = f18(i) - f1(i)*off8(i)
         f28(i) = f28(i) - f2(i)*off8(i)
         f38(i) = f38(i) - f3(i)*off8(i)
       ENDDO
 
      ELSE
 
      !-------------------------------------------!
      ! MOMENTUM ADVECTION                        !
      ! SUPG : Streamline Upwind Petrov Galerkine !
      !-------------------------------------------! 
       DO i=1,nel
        gam(i)=gam(i)/ajv(i)
        dxxp(i)=a0(i,1)*vx0(i)+a1(i,1)*vx1(i)+a2(i,1)*vx2(i)+a3(i,1)*vx3(i)  
        dxyp(i)=a0(i,2)*vx0(i)+a1(i,2)*vx1(i)+a2(i,2)*vx2(i)+a3(i,2)*vx3(i)  
        dxzp(i)=a0(i,3)*vx0(i)+a1(i,3)*vx1(i)+a2(i,3)*vx2(i)+a3(i,3)*vx3(i)  
        dyxp(i)=a0(i,1)*vy0(i)+a1(i,1)*vy1(i)+a2(i,1)*vy2(i)+a3(i,1)*vy3(i)  
        dyyp(i)=a0(i,2)*vy0(i)+a1(i,2)*vy1(i)+a2(i,2)*vy2(i)+a3(i,2)*vy3(i)  
        dyzp(i)=a0(i,3)*vy0(i)+a1(i,3)*vy1(i)+a2(i,3)*vy2(i)+a3(i,3)*vy3(i)  
        dzxp(i)=a0(i,1)*vz0(i)+a1(i,1)*vz1(i)+a2(i,1)*vz2(i)+a3(i,1)*vz3(i)  
        dzyp(i)=a0(i,2)*vz0(i)+a1(i,2)*vz1(i)+a2(i,2)*vz2(i)+a3(i,2)*vz3(i)  
        dzzp(i)=a0(i,3)*vz0(i)+a1(i,3)*vz1(i)+a2(i,3)*vz2(i)+a3(i,3)*vz3(i)  
 
        fac=eight*gam(i)
        up0(i)=one+fac*(a0(i,1)*vdx0(i)+a0(i,2)*vdy0(i)+a0(i,3)*vdz0(i))
        up1(i)=one+fac*(a1(i,1)*vdx0(i)+a1(i,2)*vdy0(i)+a1(i,3)*vdz0(i))
        up2(i)=one+fac*(a2(i,1)*vdx0(i)+a2(i,2)*vdy0(i)+a2(i,3)*vdz0(i))
        up3(i)=one+fac*(a3(i,1)*vdx0(i)+a3(i,2)*vdy0(i)+a3(i,3)*vdz0(i))
 
        fac=rho(i)*vol(i)/ajv(i)/sixty4
        f1p(i)= fac*(vdx0(i)*dxxp(i)+vdy0(i)*dxyp(i)+vdz0(i)*dxzp(i))
        f2p(i)= fac*(vdx0(i)*dyxp(i)+vdy0(i)*dyyp(i)+vdz0(i)*dyzp(i))
        f3p(i)= fac*(vdx0(i)*dzxp(i)+vdy0(i)*dzyp(i)+vdz0(i)*dzzp(i))
       ENDDO
 
        !inter22
        IF(int22>0)THEN 
          DO i=1,nel       
            IF(nint(iad22(i))==0)THEN
              !currently unplugged only for cut cell buffer
              f1p(i) = zero
              f2p(i) = zero
              f3p(i) = zero
            ENDIF        
          ENDDO
        ENDIF
 
       DO i=1,nel
        f11(i) = f11(i) - up0(i)*f1p(i)*off1(i)
        f21(i) = f21(i) - up0(i)*f2p(i)*off1(i)
        f31(i) = f31(i) - up0(i)*f3p(i)*off1(i)
        f12(i) = f12(i) - up1(i)*f1p(i)*off2(i)
        f22(i) = f22(i) - up1(i)*f2p(i)*off2(i)
        f32(i) = f32(i) - up1(i)*f3p(i)*off2(i)
        f13(i) = f13(i) - up2(i)*f1p(i)*off3(i)
        f23(i) = f23(i) - up2(i)*f2p(i)*off3(i)
        f33(i) = f33(i) - up2(i)*f3p(i)*off3(i)
        f14(i) = f14(i) - up3(i)*f1p(i)*off4(i)
        f24(i) = f24(i) - up3(i)*f2p(i)*off4(i)
        f34(i) = f34(i) - up3(i)*f3p(i)*off4(i)
        f15(i) = f15(i) - f1p(i)*off5(i)
        f25(i) = f25(i) - f2p(i)*off5(i)
        f35(i) = f35(i) - f3p(i)*off5(i)
        f16(i) = f16(i) - f1p(i)*off6(i)
        f26(i) = f26(i) - f2p(i)*off6(i)
        f36(i) = f36(i) - f3p(i)*off6(i)
        f17(i) = f17(i) - f1p(i)*off7(i)
        f27(i) = f27(i) - f2p(i)*off7(i)
        f37(i) = f37(i) - f3p(i)*off7(i)
        f18(i) = f18(i) - f1p(i)*off8(i)
        f28(i) = f28(i) - f2p(i)*off8(i)
        f38(i) = f38(i) - f3p(i)*off8(i)
       ENDDO
 
      ENDIF
 
C-----------------------------------------------
      RETURN
      END