agauge.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/.
!||====================================================================
!||    agauge                 ../engine/source/ale/agauge.F
!||--- called by ------------------------------------------------------
!||    alemain                ../engine/source/ale/alemain.F
!||--- calls      -----------------------------------------------------
!||    agaug3                 ../engine/source/ale/agauge.F
!||    agaug30                ../engine/source/ale/agauge.F
!||    agaug3q                ../engine/source/ale/agaug3q.F
!||    agaug3t                ../engine/source/ale/agaug3t.F
!||    initbuf                ../engine/share/resol/initbuf.F
!||    my_barrier             ../engine/source/system/machine.F
!||    spmd_e1vois            ../engine/source/mpi/fluid/spmd_cfd.F
!||    spmd_rbcast            ../engine/source/mpi/generic/spmd_rbcast.F
!||    spmd_sd_gaug           ../engine/source/mpi/output/spmd_sd_gaug.F
!||--- uses       -----------------------------------------------------
!||    ale_connectivity_mod   ../common_source/modules/ale/ale_connectivity_mod.F
!||    elbufdef_mod           ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    element_mod            ../common_source/modules/elements/element_mod.F90
!||    initbuf_mod            ../engine/share/resol/initbuf.F
!||====================================================================
      SUBROUTINE agauge(
     1  IPARG    ,ELBUF_STR   ,PHI     ,IXS      ,IXQ     ,
     2  X        ,ALE_CONNECT ,ITASK   ,NERCVOIS ,NESDVOIS,
     3  LERCVOIS ,LESDVOIS    ,LENCOM  ,LGAUGE   ,
     4  GAUGE    ,V           ,IGAUP   ,NGAUP    ,IXTG)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE initbuf_mod
      USE elbufdef_mod            
      USE ale_connectivity_mod
      use element_mod , only : nixs,nixq,nixtg
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.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      "com08_c.inc"
#include      "vect01_c.inc"
#include      "param_c.inc"
#include      "task_c.inc"
#include      "tabsiz_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IPARG(NPARG,NGROUP), ITASK, LENCOM,IXTG(NIXTG,NUMELTG),
     .        NERCVOIS(*),NESDVOIS(*),LERCVOIS(*),LESDVOIS(*),
     .        IXS(NIXS,NUMELS),IXQ(NIXQ,NUMELQ),LGAUGE(3,NBGAUGE),IGAUP(NBGAUGE),NGAUP(NSPMD)   
      my_real PHI(SPHI),GAUGE(LLGAUGE,NBGAUGE),X(3,NUMNOD),V(3,NUMNOD)
      TYPE(elbuf_struct_), TARGET, DIMENSION(NGROUP) :: ELBUF_STR
      TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NG, JMUL, IADR, I, II,J,JJ(6),N,IG,IS,IGAUGE,IG0,ITAG(NBGAUGE),NEL,NUMEL,NCONNECT
      my_real P,RHO,E,PA,U2,ALPHA(NBGAUGE)
      TYPE(G_BUFEL_) ,POINTER :: GBUF
C-----------------------------------------------
c     LGAUGE(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:              intantaneous PA
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                     filtered PA
c 32:  Rho                    filtered Rho
c 33:  E                      filtered E       
c 34:              ! Xpoint      !
c 35:              ! Ypoint      !
c 36:              ! Zpoint      !
c 37:              ! Butterworth !
C-----------------------------------------------
 
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      igauge=0
      DO ig=1,nbgauge
        IF(lgauge(1,ig) <= 0 .AND. lgauge(1,ig) >= -(numels+numelq+numeltg))
     .     igauge=1
      END DO
C      
      CALL my_barrier
C      
      IF(igauge == 0)RETURN     
C 
C-----------------------------
C       Recherche des elements T=0
C-----------------------------
      IF(tt==zero)THEN
       DO ng=itask+1,ngroup,nthread
        IF( iparg(5,ng) /= 1 .AND. iparg(5,ng) /= 2 .AND.
     .     (iparg(5,ng) /= 7 .AND. n2d == 0) ) cycle
        CALL initbuf(iparg    ,ng      ,
     2        mtn     ,llt     ,nft     ,iadr    ,ity     ,
     3        npt     ,jale    ,ismstr  ,jeul    ,jtur    ,
     4        jthe    ,jlag    ,jmul    ,jhbe    ,jivf    ,
     5        nvaux   ,jpor    ,jcvt    ,jclose  ,jplasol ,
     6        irep    ,iint    ,igtyp   ,israt   ,isrot   ,
     7        icsen   ,isorth  ,isorthg ,ifailure,jsms    )
        IF (iparg(8,ng) == 1)       cycle
        lft=1
        IF (iparg(5,ng) == 1) THEN
          ! solid (8N)
          numel = numels
          nconnect = 8
          CALL agaug30(lgauge,gauge,ixs ,x  ,nixs,numel,nconnect)
        ELSEIF (iparg(5,ng) == 2) THEN
          ! quad 2D
          numel = numelq
          nconnect = 4
          CALL agaug30(lgauge,gauge,ixq ,x  ,nixq,numel,nconnect)
        ELSEIF (iparg(5,ng) == 7) THEN
          ! tria 2D
          numel = numeltg
          nconnect = 3
          CALL agaug30(lgauge,gauge,ixtg ,x  ,nixtg,numel,nconnect)
        ENDIF
       ENDDO
       CALL my_barrier 
      ENDIF
C
      IF(itask==0)THEN
        DO i=1,max(numels,numelq,numeltg)
          phi(i)=zero
        END DO
C
        DO ig=1,nbgauge
          is = -lgauge(1,ig)
          IF(is > 0 .AND. is <= numels+numelq+numeltg)THEN
               phi(is)= ig 
C
C              PHI(IS)=I
c penetration ratio within element (initialisation)
            gauge(5,ig)=zero
          ENDIF
        END DO
      ENDIF
C
      CALL my_barrier
C-----------------------------
C       SPMD EXCHANGE
C-----------------------------
      IF (nspmd > 1) THEN
!$OMP SINGLE
        CALL spmd_e1vois(phi,nercvois,nesdvois,lercvois,lesdvois,lencom)
!$OMP END SINGLE
      ENDIF
C-----------------------------
C       Searching for new elements
C-----------------------------
      DO ng=itask+1,ngroup,nthread
        IF( iparg(5,ng) /= 1 .AND. iparg(5,ng) /= 2 .AND.
     .     (iparg(5,ng) /= 7 .AND. n2d == 0) ) cycle
        CALL initbuf(iparg    ,ng      ,
     2        mtn     ,llt     ,nft     ,iadr    ,ity     ,
     3        npt     ,jale    ,ismstr  ,jeul    ,jtur    ,
     4        jthe    ,jlag    ,jmul    ,jhbe    ,jivf    ,
     5        nvaux   ,jpor    ,jcvt    ,jclose  ,jplasol ,
     6        irep    ,iint    ,igtyp   ,israt   ,isrot   ,
     7        icsen   ,isorth  ,isorthg ,ifailure,jsms    )
        IF(iparg(8,ng) == 1)       cycle
         
        lft=1    
        IF (iparg(5,ng) == 1) THEN
          ! solid (8N)
          CALL agaug3(lgauge,gauge,phi,ixs  ,x  ,ale_connect )
        ELSEIF (iparg(5,ng) == 2) THEN
          ! quad 2D
          CALL agaug3q(lgauge,gauge,phi,ixq  ,x  ,ale_connect )
        ELSEIF (iparg(5,ng) == 7) THEN
          ! tria 2D
          CALL agaug3t(lgauge,gauge,phi,ixtg ,x  ,ale_connect )
        ENDIF
      END DO
C      
       CALL my_barrier
       IF(itask==0)THEN
        DO i=1,max(numels,numelq,numeltg)
          phi(i)=zero
        END DO
C
C  This array is used when several gauges are on the same element
c
        DO ig= 1,nbgauge
          itag(ig) = 0
        ENDDO
C
        DO ig=1,nbgauge
          is = -lgauge(1,ig)
          IF(is > 0 .AND. is <= (numels+numelq+numeltg))THEN
            ig0 = nint(phi(is))
            IF(ig0 > 0) THEN
              itag(ig) = ig0
            ELSE
              phi(is) = ig  
            ENDIF
           ENDIF 
        END DO
      ENDIF
C      
      CALL my_barrier
C-----------------------------
C       SPMD EXCHANGE
C-----------------------------
      IF (nspmd > 1) THEN
!$OMP SINGLE
        CALL spmd_e1vois(phi,nercvois,nesdvois,lercvois,lesdvois,lencom)
!$OMP END SINGLE
      ENDIF  
C-----------------------------
C     CALCULATION OF GAUGE MEASURES
C-----------------------------
      DO ng=itask+1,ngroup,nthread
C-----------------------------
        IF( iparg(5,ng) /= 1 .AND. iparg(5,ng) /= 2 .AND.
     .     (iparg(5,ng) /= 7 .AND. n2d == 0) ) cycle
        CALL initbuf(iparg    ,ng      ,
     2        mtn     ,llt     ,nft     ,iadr    ,ity     ,
     3        npt     ,jale    ,ismstr  ,jeul    ,jtur    ,
     4        jthe    ,jlag    ,jmul    ,jhbe    ,jivf    ,
     5        nvaux   ,jpor    ,jcvt    ,jclose  ,jplasol ,
     6        irep    ,iint    ,igtyp   ,israt   ,isrot   ,
     7        icsen   ,isorth  ,isorthg ,ifailure,jsms    )
         IF (iparg(8,ng) == 1) cycle
         gbuf => elbuf_str(ng)%GBUF
         lft=1
!
         nel = llt
         DO i=1,6
           jj(i) = nel*(i-1)
         ENDDO
!
         DO i=lft,llt
          ii=i+nft
          ig = nint(phi(ii))
          IF (ig <= 0) cycle
c         
          p = (gbuf%SIG(jj(1)+i) + gbuf%SIG(jj(2)+i) + gbuf%SIG(jj(3)+i))/three
          rho = gbuf%RHO(i)
          e   = gbuf%EINT(i)
          u2 = zero
c
          IF (iparg(5,ng) == 1) THEN
            ! Solid (8N)
            DO j=2,9
              n = ixs(j,ii)
              u2 = u2+ v(1,n)*v(1,n)
              u2 = u2+ v(2,n)*v(2,n)
              u2 = u2+ v(3,n)*v(3,n)
            ENDDO
          ELSEIF (iparg(5,ng) == 2) THEN
            ! Quad 2D
            DO j=2,5
              n = ixq(j,ii)
              u2 = u2+ v(1,n)*v(1,n)
              u2 = u2+ v(2,n)*v(2,n)
              u2 = u2+ v(3,n)*v(3,n)
            ENDDO
          ELSEIF (iparg(5,ng) == 7) THEN
            ! Tria 2D
            DO j=2,4
              n = ixtg(j,ii)
              u2 = u2+ v(1,n)*v(1,n)
              u2 = u2+ v(2,n)*v(2,n)
              u2 = u2+ v(3,n)*v(3,n)
            ENDDO
          ENDIF
          pa = p - rho*u2/sixteen
#include "lockon.inc"
          gauge(30,ig)= -p
          gauge(31,ig)= -pa
          gauge(32,ig)= rho
          gauge(33,ig)= e
#include "lockoff.inc"
          END DO ! i 
      ENDDO ! NG=ITASK+1,NGROUP,NTASK
 
      IF(itask == 0) THEN
       DO ig=1,nbgauge
        ig0= itag(ig)
        IF(ig0 > 0) THEN
#include "lockon.inc"
          gauge(30,ig)= gauge(30,ig0)
          gauge(31,ig)= gauge(31,ig0)
          gauge(32,ig)= gauge(32,ig0)
          gauge(33,ig)= gauge(33,ig0)
#include "lockoff.inc" 
        ENDIF
       ENDDO 
      ENDIF
C       
      CALL my_barrier 
C-------------
      IF(nspmd > 1) THEN        
        IF(itask == 0) THEN
          DO ig=1,nbgauge
           alpha(ig) = gauge(5,ig)
          ENDDO
          CALL spmd_sd_gaug(lgauge,gauge,igaup,ngaup)
          CALL spmd_rbcast(gauge,gauge,llgauge,nbgauge,0,2)
C        
          DO ig=1,nbgauge
           IF(gauge(5,ig) /= alpha(ig))lgauge(1,ig) = 0  
          ENDDO
        ENDIF
      ENDIF 
      RETURN
      END
!||====================================================================
!||    agaug30   ../engine/source/ale/agauge.F
!||--- called by ------------------------------------------------------
!||    agauge    ../engine/source/ale/agauge.F
!||====================================================================
      SUBROUTINE agaug30(LGAUGE ,GAUGE ,IX  ,X ,NIX, NUMEL,NCONNECT)
C-----------------------------------------------
C   Description
C-----------------------------------------------
c  Searching for element associated to the gauge
c     algorithmic complexity:quadratic (numels*nbgauge) at time 0
c     can be improved
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
#include      "vect01_c.inc"
#include      "param_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NIX,NUMEL,NCONNECT
      INTEGER IX(NIX,NUMEL),LGAUGE(3,NBGAUGE)
      my_real X(3,NUMNOD),GAUGE(LLGAUGE,NBGAUGE)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,II,J,IG,IGAUGE
      my_real XX1,YY1,ZZ1,XX2,YY2,ZZ2,XG,YG,ZG
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      DO i=1,llt
          ii=i+nft
          xx1 = ep30
          yy1 = ep30
          zz1 = ep30
          xx2 = -ep30
          yy2 = -ep30
          zz2 = -ep30
          DO j=2,nconnect+1
            xx1 = min(xx1,x(1,ix(j,ii)))
            yy1 = min(yy1,x(2,ix(j,ii)))
            zz1 = min(zz1,x(3,ix(j,ii)))
            xx2 = max(xx2,x(1,ix(j,ii)))
            yy2 = max(yy2,x(2,ix(j,ii)))
            zz2 = max(zz2,x(3,ix(j,ii)))
          ENDDO
C          
          igauge = 0 
          DO ig=1,nbgauge
            IF(lgauge(1,ig) > 0 .OR. 
     .         lgauge(1,ig) < -(numels+numelq+numeltg)) cycle
            xg  = gauge(2,ig)
            yg  = gauge(3,ig)
            zg  = gauge(4,ig)
            IF (nconnect == 8) THEN
              ! solid (8N)
              IF(xg < xx1)cycle
              IF(xg > xx2)cycle
 
              IF(yg < yy1)cycle
              IF(yg > yy2)cycle
 
              IF(zg < zz1)cycle
              IF(zg > zz2)cycle
            ELSEIF (nconnect == 3 . or. nconnect == 4) THEN
              ! tria 2D + quad 2D
              IF(yg < yy1)cycle
              IF(yg > yy2)cycle
 
              IF(zg < zz1)cycle
              IF(zg > zz2)cycle
            ENDIF
 
            igauge=ig
            IF(igauge == 0)cycle
c element associated to the gauge
c in case of multiple solutions do not need to choose the best one : AGAUG3 will update it.           
#include "lockon.inc"
            lgauge(1,igauge)=-ii
#include "lockoff.inc"
         ENDDO
      ENDDO
 
      RETURN
      END
!||====================================================================
!||    agaug3                 ../engine/source/ale/agauge.F
!||--- called by ------------------------------------------------------
!||    agauge                 ../engine/source/ale/agauge.F
!||--- uses       -----------------------------------------------------
!||    ale_connectivity_mod   ../common_source/modules/ale/ale_connectivity_mod.F
!||    element_mod            ../common_source/modules/elements/element_mod.F90
!||====================================================================
      SUBROUTINE agaug3(LGAUGE,GAUGE,PHI,IXS  ,X    ,ALE_CONNECT )
      USE ale_connectivity_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"
#include      "comlock.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
#include      "vect01_c.inc"
#include      "param_c.inc"
#include      "tabsiz_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IXS(NIXS,NUMELS),LGAUGE(3,NBGAUGE)
      my_real PHI(SPHI),X(3,NUMNOD),GAUGE(LLGAUGE,NBGAUGE)
      TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,II,J,N,N1,N2,N3,N4,IG,IAD2
      INTEGER IFACE(4,6)
      my_real ALPHA,XX0,YY0,ZZ0,XX1,YY1,ZZ1,XX2,YY2,ZZ2,A1,A2,A3,
     .        VOL,AREAP32,XG,YG,ZG
      DATA iface/ 1, 2, 3, 4,
     2            4, 3, 7, 8,
     3            8, 7, 6, 5,
     4            5, 6, 2, 1,
     5            2, 6, 7, 3,
     5            1, 4, 8, 5/
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------     
c--------------------------------------------------------------
c  Searching for the new element associated to the gauge
c--------------------------------------------------------------
      DO i=1,llt
        ii=i+nft  
        iad2 = ale_connect%ee_connect%iad_connect(ii)
        DO j=1,6
          n= ale_connect%ee_connect%connected(iad2 + j - 1)
          IF(n<=0)cycle
          ig=nint(phi(n))
          IF(ig==0)cycle
C            
            xg     = gauge(2,ig)
            yg     = gauge(3,ig)
            zg     = gauge(4,ig)
C                    
            n1 = ixs(iface(1,j)+1,ii)
            n2 = ixs(iface(2,j)+1,ii)
            n3 = ixs(iface(3,j)+1,ii)
            n4 = ixs(iface(4,j)+1,ii)
C            
            xx0 = (x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4))*fourth
            yy0 = (x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4))*fourth
            zz0 = (x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4))*fourth
            xx1 = x(1,n3)-x(1,n1)
            yy1 = x(2,n3)-x(2,n1)
            zz1 = x(3,n3)-x(3,n1)
            xx2 = x(1,n4)-x(1,n2)
            yy2 = x(2,n4)-x(2,n2)
            zz2 = x(3,n4)-x(3,n2)
c incoming vector surface
            a1  = yy1*zz2 - yy2*zz1
            a2  = xx2*zz1 - xx1*zz2
            a3  = xx1*yy2 - xx2*yy1
            vol = a1*(xg-xx0) + a2*(yg-yy0) + a3*(zg-zz0) 
            areap32 = (a1*a1+a2*a2+a3*a3)**three_over_4 
            alpha = vol/max(areap32,em20)
#include "lockon.inc"            
            IF(alpha >= zero .AND. alpha >= gauge(5,ig))THEN
c the gauge is changing of element
c element associated to the gauge
              lgauge(1,ig)=-ii
              gauge(5,ig)=alpha 
            ENDIF
#include "lockoff.inc"
          ENDDO
      ENDDO
 
 
      RETURN
      END
!||====================================================================
!||    agauge0       ../engine/source/ale/agauge.F
!||--- called by ------------------------------------------------------
!||    resol         ../engine/source/engine/resol.F
!||--- calls      -----------------------------------------------------
!||    spmd_rbcast   ../engine/source/mpi/generic/spmd_rbcast.F
!||    spmd_sd_gau   ../engine/source/mpi/output/spmd_sd_gau.F
!||--- uses       -----------------------------------------------------
!||    element_mod   ../common_source/modules/elements/element_mod.F90
!||====================================================================
      SUBROUTINE agauge0(LGAUGE ,GAUGE,X ,IXC,IGAUP,NGAUP)
      use element_mod , only : nixc
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
c     LGAUGE(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:              intantaneous PA
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                     filtered PA
c 32:  Rho                    filtered Rho
c 33:  E                      filtered E       
c 34:              ! Xpoint      !
c 35:              ! Ypoint      !
c 36:              ! Zpoint      !
c 37:              ! Butterworth !
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"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IXC(NIXC,NUMELC),LGAUGE(3,NBGAUGE),IGAUP(*),NGAUP(*)
      my_real X(3,NUMNOD),GAUGE(LLGAUGE,NBGAUGE)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER IG,IS,IGAUGE,N,N1,N2,N3,N4
      my_real XX0,YY0,ZZ0,XX1,YY1,ZZ1,XX2,YY2,ZZ2,A1,A2,A3,
     .        AA,DIST
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      igauge=0
      DO ig=1,nbgauge
        IF(lgauge(1,ig) <= 0) igauge=1
      END DO
      IF(igauge == 0)RETURN      
C      
      DO ig=1,nbgauge
        is = -lgauge(1,ig)
C
C ALL GAUGES (ALE, SPH, ..)
        IF(is >= 0)THEN
          gauge(5,ig) = -ep20
          n = lgauge(3,ig)
          IF(n > 0)THEN
            gauge(2,ig) = x(1,n)
            gauge(3,ig) = x(2,n)
            gauge(4,ig) = x(3,n)
          ELSEIF(n < 0)THEN
            n=-n
            dist = gauge(1,ig)
            n1 = ixc(2,n)
            n2 = ixc(3,n)
            n3 = ixc(4,n)
            n4 = ixc(5,n)
            xx0 = (x(1,n1)+x(1,n2)+x(1,n3)+x(1,n4))*fourth
            yy0 = (x(2,n1)+x(2,n2)+x(2,n3)+x(2,n4))*fourth
            zz0 = (x(3,n1)+x(3,n2)+x(3,n3)+x(3,n4))*fourth
            xx1 = x(1,n3)-x(1,n1)
            yy1 = x(2,n3)-x(2,n1)
            zz1 = x(3,n3)-x(3,n1)
            xx2 = x(1,n4)-x(1,n2)
            yy2 = x(2,n4)-x(2,n2)
            zz2 = x(3,n4)-x(3,n2)
cc            A1 = YY1*YY2 - ZZ1*ZZ2            
cc            A2 = ZZ1*ZZ2 - XX1*XX2            
cc            A3 = XX1*XX2 - YY1*YY2
            a1  = yy1*zz2 - yy2*zz1
            a2  = xx2*zz1 - xx1*zz2
            a3  = xx1*yy2 - xx2*yy1            
            aa = dist/sqrt(max(em20,a1*a1+a2*a2+a3*a3))
            gauge(2,ig) = xx0 + aa*a1
            gauge(3,ig) = yy0 + aa*a2
            gauge(4,ig) = zz0 + aa*a3
          ELSE
            ! point coordinates
            gauge(2,ig) = gauge(34,ig)
            gauge(3,ig) = gauge(35,ig)
            gauge(4,ig) = gauge(36,ig)
          ENDIF
        ENDIF
      END DO
C
       IF(nspmd > 1 ) THEN 
          CALL spmd_sd_gau(gauge,igaup,ngaup)
          CALL spmd_rbcast(gauge,gauge,llgauge,nbgauge,0,2)
       ENDIF 
C
      RETURN
      END