scoor3_fvm.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/.
!||====================================================================
!||    scoor3_fvm       ../engine/source/ale/alefvm/scoor3_fvm.F
!||--- called by ------------------------------------------------------
!||    sforc3           ../engine/source/elements/solid/solide/sforc3.F
!||--- calls      -----------------------------------------------------
!||    srepiso3         ../engine/source/elements/solid/solide/srepiso3.F
!||    sreploc3         ../engine/source/elements/solid/solide/sreploc3.F
!||--- uses       -----------------------------------------------------
!||    element_mod      ../common_source/modules/elements/element_mod.F90
!||    i22bufbric_mod   ../common_source/modules/interfaces/cut-cell-search_mod.F
!||    i22tri_mod       ../common_source/modules/interfaces/cut-cell-search_mod.F
!||====================================================================
      SUBROUTINE scoor3_fvm(X,IXS,V,W,GAMA,
     .   X1, X2, X3, X4, X5, X6, X7, X8,
     .   Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8,
     .   Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8,
     .   VX1, VX2, VX3, VX4, VX5, VX6, VX7, VX8,
     .   VY1, VY2, VY3, VY4, VY5, VY6, VY7, VY8,
     .   VZ1, VZ2, VZ3, VZ4, VZ5, VZ6, VZ7, VZ8,
     .   VDX1, VDX2, VDX3, VDX4, VDX5, VDX6, VDX7, VDX8,
     .   VDY1, VDY2, VDY3, VDY4, VDY5, VDY6, VDY7, VDY8,
     .   VDZ1, VDZ2, VDZ3, VDZ4, VDZ5, VDZ6, VDZ7, VDZ8,
     .   VDX,VDY,VDZ,VD2,VIS,OFFG,OFF,RHO,
     .   RHOO,NC1,NC2,NC3,NC4,NC5,NC6,NC7,NC8,NGL,MXT,NGEO,
     .   XD1, XD2  , XD3, XD4, XD5, XD6,  XD7, XD8,
     .   YD1, YD2  , YD3, YD4, YD5, YD6,  YD7, YD8,
     .   ZD1, ZD2  , ZD3, ZD4, ZD5, ZD6,  ZD7, ZD8,
     .   XDP, IPARG, NG , NEL, MOM, TAG22)
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
C 'alefvm' is related to a collocated scheme (built from FVM and based on Godunov scheme)
C  which was temporarily introduced for experimental option /INTER/TYPE22 (FSI coupling with cut cell method)
C This cut cell method is not completed, abandoned, and is not an official option.
C There is no other use for this scheme which is automatically enabled when /INTER/TYPE22 is defined (INT22>0 => IALEFVM=1).
C
C This subroutine is treating an uncut cell.
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE i22bufbric_mod
      USE i22tri_mod
      use element_mod , only : nixs
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      "vect01_c.inc"
#include      "scr05_c.inc"
#include      "param_c.inc"
#include      "inter22.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
! mom(k \in [1:3]) => gbuf%MOM((k-1)*nel + i)
 
      INTEGER IXS(NIXS,*),IPARG(NPARG,*),NG,NEL
      my_real, INTENT(IN) :: MOM(NEL,*), TAG22(NEL)
      my_real
     .  X(3,*),V(3,*),W(3,*), VIS(*),
     .  X1(*), X2(*), X3(*), X4(*), X5(*), X6(*), X7(*), X8(*),
     .  Y1(*), Y2(*), Y3(*), Y4(*), Y5(*), Y6(*), Y7(*), Y8(*),
     .  Z1(*), Z2(*), Z3(*), Z4(*), Z5(*), Z6(*), Z7(*), Z8(*),
     .  VX1(*), VX2(*), VX3(*), VX4(*), VX5(*), VX6(*), VX7(*), VX8(*),
     .  VY1(*), VY2(*), VY3(*), VY4(*), VY5(*), VY6(*), VY7(*), VY8(*),
     .  VZ1(*), VZ2(*), VZ3(*), VZ4(*), VZ5(*), VZ6(*), VZ7(*), VZ8(*),
     .  VDX1(*),VDX2(*),VDX3(*),VDX4(*),VDX5(*),VDX6(*),VDX7(*),VDX8(*),
     .  VDY1(*),VDY2(*),VDY3(*),VDY4(*),VDY5(*),VDY6(*),VDY7(*),VDY8(*),
     .  VDZ1(*),VDZ2(*),VDZ3(*),VDZ4(*),VDZ5(*),VDZ6(*),VDZ7(*),VDZ8(*),
     .  VDX(*), VDY(*), VDZ(*),VD2(*),OFFG(*),OFF(*),RHO(*),
     .  RHOO(*),GAMA(MVSIZ,6)
 
      INTEGER NC1(*), NC2(*), NC3(*), NC4(*),
     .        nc5(*), nc6(*), nc7(*), nc8(*), mxt(*), ngl(*),ngeo(*)
 
      double precision
     .  xdp(3,*),
     .  xd1(*), xd2(*), xd3(*), xd4(*), xd5(*), xd6(*), xd7(*), xd8(*),
     .  yd1(*), yd2(*), yd3(*), yd4(*), yd5(*), yd6(*), yd7(*), yd8(*),
     .  zd1(*), zd2(*), zd3(*), zd4(*), zd5(*), zd6(*), zd7(*), zd8(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, J,MXT_1,IB,NBCUT,NIN,MCELL
      my_real, DIMENSION(MVSIZ) ::
     .   RX, RY, RZ, SX, SY, SZ, TX, TY, TZ,
     .   E1X, E1Y, E1Z, E2X, E2Y, E2Z, E3X, E3Y, E3Z
      my_real
     .   off_l,wdx(mvsiz),wdy(mvsiz),wdz(mvsiz),vw,vwn,vw1,vw2,vw3,surf,sum_surf,
     .   ww1,ww2,ww3
 
      LOGICAL IsCutCell
C=======================================================================
 
      NIN       = 1
      mxt_1     = ixs(1,1)
      DO i=1,nel
        vis(i)  = zero
        ngeo(i) = ixs(10,i)
        ngl(i)  = ixs(11,i)
        mxt(i)  = mxt_1
        nc1(i)  = ixs(2,i)
        nc2(i)  = ixs(3,i)
        nc3(i)  = ixs(4,i)
        nc4(i)  = ixs(5,i)
        nc5(i)  = ixs(6,i)
        nc6(i)  = ixs(7,i)
        nc7(i)  = ixs(8,i)
        nc8(i)  = ixs(9,i)
        rhoo(i) = rho(i)
      ENDDO
 
 
        DO i=1,nel
         vdx1(i)=zero
         vdx2(i)=zero
         vdx3(i)=zero
         vdx4(i)=zero
         vdx5(i)=zero
         vdx6(i)=zero
         vdx7(i)=zero
         vdx8(i)=zero
         vdy1(i)=zero
         vdy2(i)=zero
         vdy3(i)=zero
         vdy4(i)=zero
         vdy5(i)=zero
         vdy6(i)=zero
         vdy7(i)=zero
         vdy8(i)=zero
         vdz1(i)=zero
         vdz2(i)=zero
         vdz3(i)=zero
         vdz4(i)=zero
         vdz5(i)=zero
         vdz6(i)=zero
         vdz7(i)=zero
         vdz8(i)=zero
        ENDDO
 
      off_l  = zero
      !----------------------------
      !     NODAL COORDINATES     |
      !----------------------------
      IF(iresp==1)THEN
        DO i=1,nel
          xd1(i)=xdp(1,nc1(i))
          yd1(i)=xdp(2,nc1(i))
          zd1(i)=xdp(3,nc1(i))
          xd2(i)=xdp(1,nc2(i))
          yd2(i)=xdp(2,nc2(i))
          zd2(i)=xdp(3,nc2(i))
          xd3(i)=xdp(1,nc3(i))
          yd3(i)=xdp(2,nc3(i))
          zd3(i)=xdp(3,nc3(i))
          xd4(i)=xdp(1,nc4(i))
          yd4(i)=xdp(2,nc4(i))
          zd4(i)=xdp(3,nc4(i))
          xd5(i)=xdp(1,nc5(i))
          yd5(i)=xdp(2,nc5(i))
          zd5(i)=xdp(3,nc5(i))
          xd6(i)=xdp(1,nc6(i))
          yd6(i)=xdp(2,nc6(i))
          zd6(i)=xdp(3,nc6(i))
          xd7(i)=xdp(1,nc7(i))
          yd7(i)=xdp(2,nc7(i))
          zd7(i)=xdp(3,nc7(i))
          xd8(i)=xdp(1,nc8(i))
          yd8(i)=xdp(2,nc8(i))
          zd8(i)=xdp(3,nc8(i))
        ENDDO
      ELSE
        DO i=1,nel
          xd1(i)=x(1,nc1(i))
          yd1(i)=x(2,nc1(i))
          zd1(i)=x(3,nc1(i))
          xd2(i)=x(1,nc2(i))
          yd2(i)=x(2,nc2(i))
          zd2(i)=x(3,nc2(i))
          xd3(i)=x(1,nc3(i))
          yd3(i)=x(2,nc3(i))
          zd3(i)=x(3,nc3(i))
          xd4(i)=x(1,nc4(i))
          yd4(i)=x(2,nc4(i))
          zd4(i)=x(3,nc4(i))
          xd5(i)=x(1,nc5(i))
          yd5(i)=x(2,nc5(i))
          zd5(i)=x(3,nc5(i))
          xd6(i)=x(1,nc6(i))
          yd6(i)=x(2,nc6(i))
          zd6(i)=x(3,nc6(i))
          xd7(i)=x(1,nc7(i))
          yd7(i)=x(2,nc7(i))
          zd7(i)=x(3,nc7(i))
          xd8(i)=x(1,nc8(i))
          yd8(i)=x(2,nc8(i))
          zd8(i)=x(3,nc8(i))
        ENDDO
      ENDIF
 
      CALL srepiso3(
     1   xd1,     xd2,     xd3,     xd4,
     2   xd5,     xd6,     xd7,     xd8,
     3   yd1,     yd2,     yd3,     yd4,
     4   yd5,     yd6,     yd7,     yd8,
     5   zd1,     zd2,     zd3,     zd4,
     6   zd5,     zd6,     zd7,     zd8,
     7   rx,      ry,      rz,      sx,
     8   sy,      sz,      tx,      ty,
     9   tz,      nel)
      CALL sreploc3(
     1   rx,      ry,      rz,      sx,
     2   sy,      sz,      tx,      ty,
     3   tz,      e1x,     e2x,     e3x,
     4   e1y,     e2y,     e3y,     e1z,
     5   e2z,     e3z,     llt)
 
        gama(i,1) = one
        gama(i,2) = zero
        gama(i,3) = zero
        gama(i,4) = zero
        gama(i,5) = one
        gama(i,6) = zero
 
      DO i=1,nel
       off(i) = abs(offg(i))
       off_l  = min(off_l,offg(i))
      ENDDO
 
      !copy and cast XD (DP) to X (SP) to assure coherence between XD and X
      DO i=1,nel
        x1(i)= xd1(i)
        y1(i)= yd1(i)
        z1(i)= zd1(i)
        x2(i)= xd2(i)
        y2(i)= yd2(i)
        z2(i)= zd2(i)
        x3(i)= xd3(i)
        y3(i)= yd3(i)
        z3(i)= zd3(i)
        x4(i)= xd4(i)
        y4(i)= yd4(i)
        z4(i)= zd4(i)
        x5(i)= xd5(i)
        y5(i)= yd5(i)
        z5(i)= zd5(i)
        x6(i)= xd6(i)
        y6(i)= yd6(i)
        z6(i)= zd6(i)
        x7(i)= xd7(i)
        y7(i)= yd7(i)
        z7(i)= zd7(i)
        x8(i)= xd8(i)
        y8(i)= yd8(i)
        z8(i)= zd8(i)
      ENDDO
 
      DO i=1,nel
       vx1(i)=v(1,nc1(i))
       vy1(i)=v(2,nc1(i))
       vz1(i)=v(3,nc1(i))
       vx2(i)=v(1,nc2(i))
       vy2(i)=v(2,nc2(i))
       vz2(i)=v(3,nc2(i))
       vx3(i)=v(1,nc3(i))
       vy3(i)=v(2,nc3(i))
       vz3(i)=v(3,nc3(i))
       vx4(i)=v(1,nc4(i))
       vy4(i)=v(2,nc4(i))
       vz4(i)=v(3,nc4(i))
       vx5(i)=v(1,nc5(i))
       vy5(i)=v(2,nc5(i))
       vz5(i)=v(3,nc5(i))
       vx6(i)=v(1,nc6(i))
       vy6(i)=v(2,nc6(i))
       vz6(i)=v(3,nc6(i))
       vx7(i)=v(1,nc7(i))
       vy7(i)=v(2,nc7(i))
       vz7(i)=v(3,nc7(i))
       vx8(i)=v(1,nc8(i))
       vy8(i)=v(2,nc8(i))
       vz8(i)=v(3,nc8(i))
      ENDDO
 
 
      IF(off_l<zero)THEN
        DO i=1,nel
          IF(offg(i)<zero)THEN
            vx1(i)=zero
            vy1(i)=zero
            vz1(i)=zero
            vx2(i)=zero
            vy2(i)=zero
            vz2(i)=zero
            vx3(i)=zero
            vy3(i)=zero
            vz3(i)=zero
            vx4(i)=zero
            vy4(i)=zero
            vz4(i)=zero
            vx5(i)=zero
            vy5(i)=zero
            vz5(i)=zero
            vx6(i)=zero
            vy6(i)=zero
            vz6(i)=zero
            vx7(i)=zero
            vy7(i)=zero
            vz7(i)=zero
            vx8(i)=zero
            vy8(i)=zero
            vz8(i)=zero
          ENDIF
        ENDDO
      ENDIF
 
      IF(iparg(64,ng)==1)THEN
        !boundary material
        vd2(1:nel)=zero
        RETURN
      ENDIF
 
 
       !inter22 euler : w = taken from cut surface only (moving bcs)
       !        ale   : must also take w on polyedra brick nodes which are moving too
 
 
      IF(int22==0)THEN
        IF(jale/=0)THEN
          DO i=1,nel
            wdx(i) = one_over_8*( w(1,nc1(i))+w(1,nc2(i))+w(1,nc3(i))+w(1,nc4(i))+w(1,nc5(i))+w(1,nc6(i))+w(1,nc7(i))+w(1,nc8(i)))
            wdy(i) = one_over_8*( w(2,nc1(i))+w(2,nc2(i))+w(2,nc3(i))+w(2,nc4(i))+w(2,nc5(i))+w(2,nc6(i))+w(2,nc7(i))+w(2,nc8(i)))
            wdz(i) = one_over_8*( w(3,nc1(i))+w(3,nc2(i))+w(3,nc3(i))+w(3,nc4(i))+w(3,nc5(i))+w(3,nc6(i))+w(3,nc7(i))+w(3,nc8(i)))
            vdx(i) = mom(i,1)/rho(i) - wdx(i)
            vdy(i) = mom(i,2)/rho(i) - wdy(i)
            vdz(i) = mom(i,3)/rho(i) - wdz(i)
            vd2(i) = vdx(i)*vdx(i)+vdy(i)*vdy(i)+vdz(i)*vdz(i)
          ENDDO
        ELSEIF(jeul/=0)THEN
          DO i=1,nel
            vdx(i) = mom(i,1)/rho(i)
            vdy(i) = mom(i,2)/rho(i)
            vdz(i) = mom(i,3)/rho(i)
            vd2(i) = vdx(i)*vdx(i)+vdy(i)*vdy(i)+vdz(i)*vdz(i)
          ENDDO
        ENDIF
      ELSEIF(int22/=0)THEN
       IF(jeul/=0)THEN
         DO i=1,nel
           wdx(i) = zero
           wdy(i) = zero
           wdz(i) = zero
           vdx(i) = mom(i,1)/rho(i) - wdx(i)
           vdy(i) = mom(i,2)/rho(i) - wdy(i)
           vdz(i) = mom(i,3)/rho(i) - wdz(i)
           vd2(i) = vdx(i)*vdx(i)+vdy(i)*vdy(i)+vdz(i)*vdz(i)
          ENDDO !next I
       ELSEIF(jale/=0)THEN
         DO i=1,nel
          ib = 0
          IF(int22/=0)ib = nint(tag22(i))
          iscutcell = .false.
          IF(ib/=0)THEN
            nbcut=brick_list(nin,ib)%NBCUT
            IF(nbcut>0)iscutcell=.true.
          ENDIF
          IF(iscutcell .EQV. .false.)THEN
            wdx(i) = one_over_8*( w(1,nc1(i))+w(1,nc2(i))+w(1,nc3(i))+w(1,nc4(i))+w(1,nc5(i))+w(1,nc6(i))+w(1,nc7(i))+w(1,nc8(i)))
            wdy(i) = one_over_8*( w(2,nc1(i))+w(2,nc2(i))+w(2,nc3(i))+w(2,nc4(i))+w(2,nc5(i))+w(2,nc6(i))+w(2,nc7(i))+w(2,nc8(i)))
            wdz(i) = one_over_8*( w(3,nc1(i))+w(3,nc2(i))+w(3,nc3(i))+w(3,nc4(i))+w(3,nc5(i))+w(3,nc6(i))+w(3,nc7(i))+w(3,nc8(i)))
            vdx(i) = mom(i,1)/rho(i) - wdx(i)
            vdy(i) = mom(i,2)/rho(i) - wdy(i)
            vdz(i) = mom(i,3)/rho(i) - wdz(i)
            vd2(i) = vdx(i)*vdx(i)+vdy(i)*vdy(i)+vdz(i)*vdz(i)
          ELSE
            mcell = brick_list(nin,ib)%MainID
            IF (mcell==0)THEN
              wdx(i) = zero
              wdy(i) = zero
              wdz(i) = zero
              vdx(i) = zero
              vdy(i) = zero
              vdz(i) = zero
              vd2(i) = zero
            ELSE
              vw     = zero
              vdx(i) = zero
              vdy(i) = zero
              vdz(i) = zero
              vd2(i) = zero
              j=maxloc(brick_list(nin,ib)%POLY(mcell)%FACE(1:6)%SURF,1)
              surf     = brick_list(nin,ib)%POLY(mcell)%FACE(j)%SURF
              sum_surf = surf
              vw1 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Vel(1)
              vw2 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Vel(2)
              vw3 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Vel(3)
              ww1 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%W(1)
              ww2 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%W(2)
              ww3 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%W(3)
              vw1 = vw1-ww1
              vw2 = vw2-ww2
              vw3 = vw3-ww3
              vwn = vw1*vw1 + vw2*vw2 + vw3*vw3
              vdx(i) = surf*vw1
              vdy(i) = surf*vw2
              vdz(i) = surf*vw3
              vd2(i) = vwn
              DO j=1,6
                surf   = brick_list(nin,ib)%POLY(mcell)%FACE(j)%SURF
                IF(surf>zero)THEN
                  sum_surf = sum_surf + surf
                  vw1 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Vel(1)
                  vw2 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Vel(2)
                  vw3 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%Vel(3)
                  ww1 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%W(1)
                  ww2 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%W(2)
                  ww3 = brick_list(nin,ib)%POLY(mcell)%FACE(j)%W(3)
                  vw1 = vw1-ww1
                  vw2 = vw2-ww2
                  vw3 = vw3-ww3
                  vwn = vw1*vw1 + vw2*vw2 + vw3*vw3
                  vdx(i) = vdx(i)+surf*vw1
                  vdy(i) = vdy(i)+surf*vw2
                  vdz(i) = vdz(i)+surf*vw3
                  vd2(i) = max(vd2(i),vwn)
                ENDIF
              ENDDO
              IF(sum_surf>zero)THEN
                 vdx(i) = vdx(i) / sum_surf
                 vdy(i) = vdy(i) / sum_surf
                 vdz(i) = vdz(i) / sum_surf
               ENDIF
            ENDIF   !(MCELL==0)
          ENDIF  !(IsCutCell .EQV. .FALSE.)
         ENDDO !next I
        ENDIF
      ENDIF
 
 
      if(ibug22_vd2 == -1)then
        do i=lft,llt
          write (*,fmt='(A,I10,4F30.16)') "VD x,y,z,2,   ID=",ixs(11,i), vdx(i),vdy(i),vdz(i),vd2(i)
        enddo
      endif
 
      ! VD2(1:NEL) = ZERO
      ! VDX(1:NEL) = ZERO
      ! VDY(1:NEL) = ZERO
      ! VDZ(1:NEL) = ZERO
 
 
      !Here due to convection subroutine aconve.f
      !        RHO = rho.vol
      !        MOM = rho.U.vol
      ! => MOM/RHO = U
 
 
 
C-----------
      RETURN
      END