i9wal3.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/.
!||====================================================================
!||    i9wal3            ../engine/source/interfaces/int09/i9wal3.F
!||--- called by ------------------------------------------------------
!||    i9wale            ../engine/source/interfaces/int09/i9wale.F
!||--- calls      -----------------------------------------------------
!||    bcs2              ../engine/source/constraints/general/bcs/bcs2.F
!||    i9grd3            ../engine/source/interfaces/int09/i9grd3.F
!||    shapeh            ../engine/source/ale/inter/shapeh.F
!||    spmd_glob_dsum9   ../engine/source/mpi/interfaces/spmd_th.F
!||    spmd_glob_isum9   ../engine/source/mpi/interfaces/spmd_th.F
!||    spmd_ibcast       ../engine/source/mpi/generic/spmd_ibcast.F
!||    spmd_rbcast       ../engine/source/mpi/generic/spmd_rbcast.F
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod      ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    element_mod       ../common_source/modules/elements/element_mod.F90
!||====================================================================
      SUBROUTINE i9wal3(X    ,V         ,W         ,A     ,CRST     ,
     2          NSV      ,ILOC      ,IRTL      ,ICODE     ,ISKEW    ,
     3          SKEW     ,MSR       ,LMSR      ,NSEG      ,IRECTS   ,
     4          IRECT    ,UPW       ,IXS       ,ELBUF_TAB,
     5          IPARG    ,PM        ,NALE      ,EE        ,IELES    ,
     6          IELEM    ,TSTIF     ,INTTH     ,IEULT     ,STENS    ,
     7          NOR      ,ISIZES    ,ISIZEM    ,NRTS, NRTM, NSN,NMN  )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE elbufdef_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   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "com08_c.inc"
#include      "scr08_a_c.inc"
#include      "param_c.inc"
#include      "task_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: NRTS, NRTM, NSN,NMN 
      INTEGER IRECT(4,*), NSV(NSN), ILOC(*), IRTL(NSN), ICODE(*), ISKEW(*),
     .   MSR(*), IRECTS(4,*), LMSR(*), NSEG(*),IXS(NIXS,*),
     .   IPARG(NPARG,*), IELES(*), IELEM(*),NALE(*) ,
     .   INTTH, IEULT, ISIZES, ISIZEM
C     REAL
      my_real
     .   upw, tstif,ttt, stens,
     .   x(3,*), v(3,*), w(3,*), a(3,*), crst(2,*), skew(lskew,*),
     .   pm(npropm,*),ee(*),nor(3,*)
      TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER II, N, L, JJ, NN, LL, LL1, LL2, LG, KK, KKK, K1, K2,
     .   I1, I2, IERR, IGROU, IELN,
     .   ixx(4), iperm(0:5),
     .   itemp(2), is, im, ilen,
     .   tags(isizes),tagm(isizem), lists(isizes),listm(isizem),
     .   icomerr(isizem+isizes),icomngr(isizem+isizes),
     .   icomnel(isizem+isizes)
C     REAL
      my_real
     .   h(4), vmx, vmy, vmz, vx, vy, vz, vv, nx, ny, nz, vt,
     .   nnx, nny, nnz, fac, p, x1, y1, z1,x2, y2, z2, tx, ty, tz,
     .   efric, vols, volm, ts, tm ,tstifm, tstifs, dvn, tt2, tt3,
     .   tstift, phi, areas, aream, vn, wn, t2x, t2y, t2z, t2t,
     .   t3x, t3y, t3z, stensx, stensy, stensz,
     .   comarea(isizem+isizes),comstf(isizem+isizes),
     .   comt(isizem+isizes),comvol(isizem+isizes),combuf(isizem+isizes)
 
c
      DATA iperm/ 4, 1, 2, 3, 4, 1/
C-----------------------------------------------
C
C preparation phase for spmd
C
      IF(intth/=zero) THEN
       IF(ispmd==0) THEN
        DO ii = 1, nrts
          tags(ii) = 0
        ENDDO
        DO ii = 1, nrtm
          tagm(ii) = 0
        ENDDO
        is = 0
        im = 0
        DO ii = 1, nsn
          l = irtl(ii)
          IF(iloc(ii)>0.AND.nmn>0.AND.tagm(l)==0)THEN
            im = im + 1
            listm(im) = l
            tagm(l) = im
            ll1=nseg(ii)
            ll2=nseg(ii+1)-1
            DO ll=ll1,ll2
              lg  = lmsr(ll)
              IF(tags(lg)==0) THEN
                is = is + 1
                lists(is) = lg
                tags(lg) = is
              ENDIF
            ENDDO
          ENDIF
        ENDDO
C
C Elem lists compacting
C
        itemp(1) = is
        itemp(2) = im
       ENDIF
C
C Sending second/mains facet list in contact
C
       IF(nspmd > 1) THEN
         CALL spmd_ibcast(itemp,itemp,1,2,0,2)
         is = itemp(1)
         im = itemp(2)
         ilen = im+is
         CALL spmd_ibcast(listm,listm,1,im,0,2)
         CALL spmd_ibcast(lists,lists,1,is,0,2)
       END IF
       DO ii = 1, im
         l = listm(ii)
         ix(1) = msr(irect(1,l))
         ix(2) = msr(irect(2,l))
         ix(3) = msr(irect(3,l))
         ix(4) = msr(irect(4,l))
         IF(ielem(l)>0) THEN
           CALL i9grd3(
     1       ierr ,comarea(ii),comstf(ii),comt(ii),comvol(ii),
     2       ielem(l)   ,x         ,ixs(1,ielem(l)), ix,
     3       iparg,pm         ,elbuf_tab ,igrou   ,ieln    )
           icomerr(ii) = ierr
           icomngr(ii) = igrou
           icomnel(ii) = ieln
         ELSE
           comarea(ii) = zero
           comstf(ii)  = zero
           comt(ii)    = zero
           comvol(ii)  = zero
           icomerr(ii) = 0
           icomngr(ii) = 0
           icomnel(ii) = 0
         ENDIF
         combuf(ii) = zero
       ENDDO
C
       DO ii = 1, is
         l = lists(ii)
         ixx(1)=nsv(irects(1,l))
         ixx(2)=nsv(irects(2,l))
         ixx(3)=nsv(irects(3,l))
         ixx(4)=nsv(irects(4,l))
         IF(ieles(l)>0) THEN
          CALL i9grd3(
     1     ierr ,comarea(im+ii),comstf(im+ii),comt(im+ii),comvol(im+ii),
     2     ieles(l)      ,x,ixs(1,ieles(l))        ,ixx          ,
     3     iparg,pm            ,elbuf_tab    ,igrou   ,ieln    )
           icomerr(im+ii) = ierr
           icomngr(im+ii) = igrou
           icomnel(im+ii) = ieln
         ELSE
           comarea(im+ii) = zero
           comstf(im+ii)  = zero
           comt(im+ii)    = zero
           comvol(im+ii)  = zero
           icomerr(im+ii) = 0
           icomngr(im+ii) = 0
           icomnel(im+ii) = 0
         ENDIF
         combuf(im+ii) = zero
       ENDDO
C
       IF (nspmd > 1) THEN
C
C gather values
C
         CALL spmd_glob_dsum9(comarea,ilen)
         CALL spmd_glob_dsum9(comstf,ilen)
         CALL spmd_glob_dsum9(comt,ilen)
         CALL spmd_glob_dsum9(comvol,ilen)
         CALL spmd_glob_isum9(icomerr,ilen)
                 CALL spmd_glob_isum9(icomngr,ilen)
                 CALL spmd_glob_isum9(icomnel,ilen)
C   node part on P0 only
         IF(ispmd/=0) GOTO 900
       END IF
C   interface processed by p0
      ELSE
       IF(ispmd/=0) RETURN
      ENDIF
C
      DO 800 ii=1,nsn
      ll1=nseg(ii)
      ll2=nseg(ii+1)-1
      n=nsv(ii)
      IF(iloc(ii)>0.AND.nmn>0)THEN
C---------------------------------
C       CONTACT
C---------------------------------
        l=irtl(ii)
        DO 10 jj=1,4
        nn=irect(jj,l)
 10     ix(jj)=msr(nn)
C
        CALL shapeh(h,crst(1,ii),crst(2,ii))
C---------------------------------
C       mesh velocity
C---------------------------------
        vmx=zero
        vmy=zero
        vmz=zero
C
        DO 30 jj=1,4
        vmx=vmx+w(1,ix(jj))*h(jj)
        vmy=vmy+w(2,ix(jj))*h(jj)
 30     vmz=vmz+w(3,ix(jj))*h(jj)
C
        dvn = (vmx-w(1,n)) * nor(1,ii)
     .        + (vmy-w(2,n)) * nor(2,ii)
     .        + (vmz-w(3,n)) * nor(3,ii)
        w(1,n) = w(1,n) + dvn * nor(1,ii)
        w(2,n) = w(2,n) + dvn * nor(2,ii)
        w(3,n) = w(3,n) + dvn * nor(3,ii)
C
C---------------------------------
C       PONT THERMIQUE
C---------------------------------
        IF(intth/=zero)THEN
           kk = tagm(l)
           efric = half * ee(ii) / (ll2-ll1+1)
           ierr = icomerr(kk)
           aream = comarea(kk)
           tstifm = comstf(kk)
           tm = comt(kk)
           volm = comvol(kk)
           IF(ierr==0) THEN
             DO ll = ll1,ll2
               lg = lmsr(ll)
               jj = tags(lg) + im
               ierr = icomerr(jj)
               areas = comarea(jj)
               tstifs = comstf(jj)
               ts = comt(jj)
               vols = comvol(jj)
               IF(ierr==0) THEN
                 tstift = tstifm + tstifs + tstif
                 phi = areas * dt1 * (tm-ts) / tstift
                 combuf(jj) = combuf(jj)
     +                      + (efric+phi)/vols
                 combuf(kk) = combuf(kk)
     +                      + (efric-phi)/volm
               ENDIF
             ENDDO
           ENDIF
c         ENDIF
        ENDIF
C
      ELSEIF(iloc(ii)<0.OR.nmn==0)THEN
C------------------------------------------------------
C       loop over facets connected to node ii
C------------------------------------------------------
        iloc(ii) = -iloc(ii)
C
        vx = v(1,n) - w(1,n)
        vy = v(2,n) - w(2,n)
        vz = v(3,n) - w(3,n)
        vv = max(em30,sqrt(vx**2+vy**2+vz**2))
        nnx = zero
        nny = zero
        nnz = zero
C------------------------------------------------------
C       loop over facets connected to node ii
C------------------------------------------------------
        DO 300 ll=ll1,ll2
         lg=lmsr(ll)
         DO 200 kkk=1,4
          kk=kkk
 200     IF(irects(kk,lg)==ii) GO TO 250
 250     CONTINUE
C------------------------------------------------------
C         calculation of the normal with upwind on the upstream side
C------------------------------------------------------
          k1 = iperm(kk-1)
          k2 = iperm(kk+1)
          i1 = nsv(irects(k1,lg))
          i2 = nsv(irects(k2,lg))
          x1 = x(1,i1) - x(1,n)
          y1 = x(2,i1) - x(2,n)
          z1 = x(3,i1) - x(3,n)
          x2 = x(1,i2) - x(1,n)
          y2 = x(2,i2) - x(2,n)
          z2 = x(3,i2) - x(3,n)
          tx = x1 + x2
          ty = y1 + y2
          tz = z1 + z2
          ttt = max(em30,sqrt(tx**2+ty**2+tz**2))
          vt = v(1,n)*tx + v(2,n)*ty + v(3,n)*tz
          p = onep0001 - upw*(half + sign(half,vt))
          nx = y1 * z2 - z1 * y2
          ny = z1 * x2 - x1 * z2
          nz = x1 * y2 - y1 * x2
C          FAC = P / MAX(EM30,SQRT(NX**2+NY**2+NZ**2))
          fac = p
          nnx = nnx + nx*fac
          nny = nny + ny*fac
          nnz = nnz + nz*fac
C-------------------------------------
C       surface tension
C-------------------------------------
          IF(stens>zero)THEN
            t2x = -x1 + x2
            t2y = -y1 + y2
            t2z = -z1 + z2
            tt2 = max(em30,t2x**2+t2y**2+t2z**2)
            t2t = (t2x*tx + t2y*ty +t2z*tz) / tt2
            t3x = tx - t2x * t2t
            t3y = ty - t2y * t2t
            t3z = tz - t2z * t2t
            tt3 = stens * sqrt(tt2/max(em30,t3x**2+t3y**2+t3z**2))
            stensx = t3x * tt3
            stensy = t3y * tt3
            stensz = t3z * tt3
            a(1,n)  = a(1,n)  + stensx
            a(2,n)  = a(2,n)  + stensy
            a(3,n)  = a(3,n)  + stensz
          ENDIF
 300    CONTINUE
        fac = max(em30,sqrt(nnx**2+nny**2+nnz**2))
        nnx = nnx/fac
        nny = nny/fac
        nnz = nnz/fac
C
C---------------------------------
C       W LAGRANGIEN ACCORDING TO N
C---------------------------------
C---------------------------------
C       grid boundary conditions
C---------------------------------
        IF(icode(n)/=0)THEN
          dvn = vx * nnx + vy * nny + vz * nnz
          w(1,n) = w(1,n) + dvn * nnx
          w(2,n) = w(2,n) + dvn * nny
          w(3,n) = w(3,n) + dvn * nnz
          CALL bcs2(w(1,n),skew(1,iskew(n)),iskew(n),icode(n))
          vn = v(1,n)*nnx + v(2,n)*nny + v(3,n)*nnz
          wn = w(1,n)*nnx + w(2,n)*nny + w(3,n)*nnz
C-------------------------------------
C         W LAGRANGIEN ACCORDING TO N + BCS
C-------------------------------------
          IF(abs(wn)>em30)THEN
            fac = vn / wn
            w(1,n) = w(1,n) * fac
            w(2,n) = w(2,n) * fac
            w(3,n) = w(3,n) * fac
          ENDIF
        ELSEIF(ieult/=0)THEN
C-------------------------------------
C         W LAGRANGIEN ACCORDING TO N
C         W EULERIEN ACCORDING TO T
C---------------------------------
          vn = v(1,n)*nnx + v(2,n)*nny + v(3,n)*nnz
          w(1,n) = vn * nnx
          w(2,n) = vn * nny
          w(3,n) = vn * nnz
        ELSE
C-------------------------------------
C         W LAGRANGIEN ACCORDING TO N
C         LIBRE ACCORDING TO T
C---------------------------------
          dvn = vx * nnx + vy * nny + vz * nnz
          w(1,n) = w(1,n) + dvn * nnx
          w(2,n) = w(2,n) + dvn * nny
          w(3,n) = w(3,n) + dvn * nnz
        ENDIF
      ENDIF
C
  800 CONTINUE
C
C finalization phase for spmd
C
  900 CONTINUE
      IF(intth/=zero) THEN
       IF(nspmd > 1) THEN
C
C Envoi buffer elements updates
C
         CALL spmd_rbcast(combuf,combuf,1,ilen,0,2)
       END IF
C
C Local Elbuf update
C
       DO ii = 1, im
         l = listm(ii)
         IF(ielem(l)>0) THEN
           igrou = icomngr(ii)
           ieln  = icomnel(ii)
           elbuf_tab(igrou)%GBUF%EINT(ieln) = 
     .     elbuf_tab(igrou)%GBUF%EINT(ieln) + combuf(ii)
         ENDIF
       ENDDO
C
       DO ii = 1, is
         l =  lists(ii)
         IF(ieles(l)>0) THEN
           igrou = icomngr(im+ii)
           ieln  = icomnel(im+ii)
           elbuf_tab(igrou)%GBUF%EINT(ieln) = 
     .     elbuf_tab(igrou)%GBUF%EINT(ieln) + combuf(im+ii)
         ENDIF
       ENDDO
      ENDIF
C
      RETURN
      END