spgauge.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/.
!||====================================================================
!||    spgauge         ../engine/source/elements/sph/spgauge.F
!||--- called by ------------------------------------------------------
!||    forintp         ../engine/source/elements/forintp.F
!||--- calls      -----------------------------------------------------
!||    spmd_exch_fr6   ../engine/source/mpi/kinematic_conditions/spmd_exch_fr6.F
!||    sum_6_float     ../engine/source/system/parit.F
!||    weight0         ../engine/source/elements/sph/weight.F
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod    ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    initbuf_mod     ../engine/share/resol/initbuf.F
!||    sphbox          ../engine/share/modules/sphbox.F
!||====================================================================
      SUBROUTINE spgauge(LGAUGE  ,GAUGE ,KXSP   ,IXSP  ,
     1                   SPBUF   ,IPARG ,ELBUF_TAB,ISPSYM ,XSPSYM,
     2                   NOD2SP  ,X       ,ITASK  ,WA     ,WASIGSM,
     3                   WAR     ,SPHG_F6)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE initbuf_mod
      USE sphbox
      USE elbufdef_mod            
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      "com04_c.inc"
#include      "param_c.inc"
#include      "sphcom.inc"
#include      "parit_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LGAUGE(3,*), KXSP(NISP,*), IXSP(KVOISPH,*),
     .        IPARG(NPARG,*), ISPSYM(NSPCOND,*), NOD2SP(*), ITASK
C     REAL
      my_real
     .   gauge(llgauge,*), spbuf(nspbuf,*), x(3,*),xspsym(3,*),
     .   wa(kwasph,*), wasigsm(6,*), war(10,*)
      TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
      DOUBLE PRECISION SPHG_F6(4,6,NBGAUGE)
C-----------------------------------------------
c     GAUGE(3,*)
c 1:  -Isolid           -(NUMELS_G+1) if SPH gauge
c 2:  GaugeId
c 3:  +Node or -Shell
c
c     => GAUGE(LLGAUGE,*), LLGAUGE = 37
c 1:  Dist (distance from Shell)     Dist (distance from Shell)
c 2:  XG           XG
c 3:  YG           YG
c 4:  ZG           ZG
c 5:  Alpha (Solid penetration ratio)     not yet used
c 6:               XSAV (SPH sorting)
c 7:               YSAV (SPH sorting)
c 8:               ZSAV (SPH sorting)
c 9:               FF (sph only)
c 10:              intantaneous Pressure
c 11:              not yet available
c 12:              intantaneous Rho
c 13:              intantaneous E
c 14:              ! Butterworth !
c 15:              ! Butterworth !
c 16:              ! Butterworth !
c 17:              ! Butterworth !
c 18:              ! Butterworth !
c 19:              ! Butterworth !
c 20:                    ! Butterworth !
c 21:                    ! Butterworth !
c 22:                    ! Butterworth !
c 23:                    ! Butterworth !  
c 24:              ! Butterworth !
c 25:              ! Butterworth !
c 26:              ! Butterworth !
c 27:              ! Butterworth !
c 28:              ! Butterworth !
c 29:              ! Butterworth !
c 30:  Pressure                 filtered Pressure
c 31:  PA                     not yet available
c 32:  Rho                    filtered Rho
c 33:  E                      filtered E       
c 34:              ! Butterworth !
c 35:              ! Butterworth !
c 36:              ! Butterworth !
c 37:              ! Butterworth !
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,N,INOD,JNOD,J,NVOIS,M,
     .        NVOISS,SM,JS,NC,NS,NN,
     .        IG,NEL,KAD,NG,K,KK,SFGAUGE,FR_RL(NSPMD+2)
      my_real
     .       xi,yi,zi,di,rhoi,xj,yj,zj,dj,rhoj,vj,
     .       wght,alphai,wcompi,
     .       pp,ee,press,ener,rho
      TYPE(g_bufel_)  ,POINTER :: GBUF
      my_real
     .       , DIMENSION(:), ALLOCATABLE :: FGAUGE1,FGAUGE2,FGAUGE3,FGAUGE4        
C-----------------------------------------------
      fr_rl(:) = 1
!$OMP DO SCHEDULE(DYNAMIC,1)
      DO ig=1,nbgauge
        IF(lgauge(1,ig) > -(numels+1)) cycle
C------
        n=numsph+ig
        xi =gauge(2,ig)
        yi =gauge(3,ig)
        zi =gauge(4,ig)
        wcompi =zero
        press  =zero
        ener   =zero
        rho    =zero
C------
          nvois=kxsp(4,n)
C------
        sfgauge = kxsp(4,n) + kxsp(6,n)
        ALLOCATE(fgauge1(sfgauge),fgauge2(sfgauge),fgauge3(sfgauge),fgauge4(sfgauge))
C------
          DO j=1,nvois
            jnod=ixsp(j,n)
            IF(jnod>0)THEN
              m=nod2sp(jnod)
            IF(kxsp(2,m)<=0) cycle
              xj=x(1,jnod)
              yj=x(2,jnod)
              zj=x(3,jnod)
              dj  =spbuf(1,m)
              rhoj=spbuf(2,m)
              vj=spbuf(12,m)/max(em20,rhoj)
              CALL weight0(xi,yi,zi,xj,yj,zj,dj,wght)
            wcompi=vj*wght
            pp=-third*(wa(1,m)+wa(2,m)+wa(3,m))
            fgauge1(j) = wcompi
            fgauge2(j) = vj*wght*pp
            fgauge3(j) = vj*wght*rhoj
            fgauge4(j) = ener
            END IF
          ENDDO
C------
C       partie symetrique.
        nvoiss=kxsp(6,n)
        k = nvois
        DO j=kxsp(5,n)+1,kxsp(5,n)+nvoiss
          js=ixsp(j,n)
          IF(js>0)THEN
            sm=js/(nspcond+1)
            nc=mod(js,nspcond+1)
            js=ispsym(nc,sm)
            xj =xspsym(1,js)
            yj =xspsym(2,js)
            zj =xspsym(3,js)
            dj  =spbuf(1,sm)
            rhoj=spbuf(2,sm)
            vj=spbuf(12,sm)/max(em20,rhoj)
            CALL weight0(xi,yi,zi,xj,yj,zj,dj,wght)
            wcompi=vj*wght
            pp=-third*(wasigsm(1,js)+wasigsm(2,js)+wasigsm(3,js))
            k = k+1
            fgauge1(k) = wcompi
            fgauge2(k) = vj*wght*pp
            fgauge3(k) = vj*wght*rhoj
            fgauge4(k) = ener
          END IF
        ENDDO
C
C Traitement Parith/ON avant echange
C
        DO k = 1, 6
          sphg_f6(1,k,ig) = zero
          sphg_f6(2,k,ig) = zero
          sphg_f6(3,k,ig) = zero
          sphg_f6(4,k,ig) = zero
        END DO
C
        IF(iparit > 0)THEN
          CALL sum_6_float(1  ,sfgauge  ,fgauge1, sphg_f6(1,1,ig),4)
          CALL sum_6_float(1  ,sfgauge  ,fgauge2, sphg_f6(2,1,ig),4)
          CALL sum_6_float(1  ,sfgauge  ,fgauge3, sphg_f6(3,1,ig),4)
          CALL sum_6_float(1  ,sfgauge  ,fgauge4, sphg_f6(4,1,ig),4)
        ELSE
          DO i=1,sfgauge
            k = 1
            sphg_f6(1,k,ig) = sphg_f6(1,k,ig) + fgauge1(i)
            sphg_f6(2,k,ig) = sphg_f6(2,k,ig) + fgauge2(i)
            sphg_f6(3,k,ig) = sphg_f6(3,k,ig) + fgauge3(i)
            sphg_f6(4,k,ig) = sphg_f6(4,k,ig) + fgauge4(i)
          END DO
        END IF
C------
        DEALLOCATE(fgauge1,fgauge2,fgauge3,fgauge4)
C-----------------
      END DO
!$OMP END DO
C
C-------------
      IF(nspmd > 1) THEN    
!$OMP SINGLE
           CALL spmd_exch_fr6(fr_rl,sphg_f6,4*nbgauge*6)
!$OMP END SINGLE
      ENDIF
C-------------
C
C Traitement Parith/ON apres echange
C
!$OMP DO SCHEDULE(DYNAMIC,1)
      DO ig=1,nbgauge
        IF(lgauge(1,ig) > -(numels+1)) cycle
        wcompi =  sphg_f6(1,1,ig)+sphg_f6(1,2,ig)+sphg_f6(1,3,ig)+
     +            sphg_f6(1,4,ig)+sphg_f6(1,5,ig)+sphg_f6(1,6,ig)
        press =   sphg_f6(2,1,ig)+sphg_f6(2,2,ig)+sphg_f6(2,3,ig)+
     +            sphg_f6(2,4,ig)+sphg_f6(2,5,ig)+sphg_f6(2,6,ig)
        rho  =    sphg_f6(3,1,ig)+sphg_f6(3,2,ig)+sphg_f6(3,3,ig)+
     +            sphg_f6(3,4,ig)+sphg_f6(3,5,ig)+sphg_f6(3,6,ig)      
        alphai=one/max(em20,wcompi)
        press=press*alphai
        rho  =rho*alphai
        gauge(10,ig)=press
        gauge(12,ig)=rho
C energy not yet available
        gauge(13,ig)=zero    
      END DO
!$OMP END DO      
C
      RETURN
      END
!||====================================================================
!||    spgauge_f   ../engine/source/elements/sph/spgauge.f
!||--- called by ------------------------------------------------------
!||    resol       ../engine/source/engine/resol.F
!||====================================================================
      SUBROUTINE spgauge_f(P,FF,P0,P1,P2,N)
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      "com08_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER N
C     REAL
      my_real
     . p(n),p0(n,2),p1(n,2),p2(n,2),ff
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER J,IFIRST
      my_real
     . PI1,PI8,PI38,SPI8,SPI38,C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,
     , x1,x2,x3,y1,y2,y3,z1,z2,z3,d,dd,d2,dp,e,g,f
C
      IF(ff==zero)THEN
        DO j=1,n
          p2(j,1) = p(j)
        ENDDO
        RETURN
      END IF
C
      f = min(ff,zep4/dt2)
C-----------------------------------------------
C INITIALIALISATION DES COEFFICIENTS DU FILTRE
C-----------------------------------------------
      pi1   = two*atan2(one,zero)
      pi8  = pi1*one_over_8
      pi38 = three*pi8
      spi8  = sin(pi8)
      spi38 = sin(pi38)
C
      d  = tan(pi1*f*dt2)
      dd = d*d
      d2 = two*d
      dp = one + dd
      e  = d2*spi8
      g  = e + dp
      g  = one/g
C
      c0 = dd * g
      c1 = two* c0
      c2 = c0
      c3 = two * g - c1
      c4 = (e - dp) * g
C
      e  = d2*spi38
      g  = e + dp
      g  = one/g
C
      c5 = dd * g
      c6 = two * c5
      c7 = c5
      c8 = two * g - c6
      c9 = (e - dp) * g
C-----------------------------------------------
C FILTRAGE
C-----------------------------------------------
      DO j=1,n
        x1 = p0(j,2)
        x2 = p0(j,1)
        x3 = p(j)
        y1 = p1(j,2)
        y2 = p1(j,1)
        y3 = c0 * x3 + c1 * x2 + c2 * x1
     .     + c3 * y2 + c4 * y1
        z1 = p2(j,2)
        z2 = p2(j,1)
        z3 = c5 * y3 + c6 * y2 + c7 * y1
     .     + c8 * z2 + c9 * z1
C
        p0(j,2) = x2
        p0(j,1) = x3
        p1(j,2) = y2
        p1(j,1) = y3
        p2(j,2) = z2
        p2(j,1) = z3
      ENDDO
C
      RETURN
      END