i7trivox.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/.
!||====================================================================
!||    i7trivox        ../engine/source/interfaces/intsort/i7trivox.F
!||--- called by ------------------------------------------------------
!||    i7buce_vox      ../engine/source/interfaces/intsort/i7buce.F
!||--- calls      -----------------------------------------------------
!||    ancmsg          ../engine/source/output/message/message.F
!||    arret           ../engine/source/system/arret.F
!||    i7sto           ../engine/source/interfaces/intsort/i7sto.F
!||    my_barrier      ../engine/source/system/machine.F
!||    spmd_oldnumcd   ../engine/source/mpi/interfaces/spmd_i7tool.F
!||--- uses       -----------------------------------------------------
!||    message_mod     ../engine/share/message_module/message_mod.F
!||    tri7box         ../engine/share/modules/tri7box.F
!||====================================================================
      SUBROUTINE i7trivox(
     1      NSN    ,RENUM       ,NSNR    ,ISZNSNR  ,I_MEM    ,
     2      IRECT  ,X           ,STF     ,STFN     ,XYZM     ,
     3      NSV    ,II_STOK     ,CAND_N  ,ESHIFT   ,CAND_E   ,
     4      MULNSN ,NOINT       ,TZINF   ,GAP_S_L  ,GAP_M_L  ,
     5      VOXEL  ,NBX         ,NBY     ,NBZ      ,INTTH    ,
     6      INACTI ,IFQ         ,CAND_A  ,CAND_P   ,IFPEN    ,
     7      NRTM   ,NSNROLD     ,IGAP    ,GAP      ,GAP_S    ,
     8      GAP_M  ,GAPMIN      ,GAPMAX  ,MARGE    ,CURV_MAX ,
     9      NIN    ,ITASK       ,BGAPSMX ,KREMNOD  ,REMNOD   ,
     A      ITAB   ,FLAGREMNODE ,DRAD    ,ITIED    ,CAND_F   ,
     B      DGAPLOAD,REMOTE_S_NODE,LIST_REMOTE_S_NODE,
     C      TOTAL_NB_NRTM,INTHEAT,IDT_THERM,NODADT_THERM)
C============================================================================
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
      USE tri7box
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     parameter setting the size for the vector (orig version is 128)
      INTEGER NVECSZ
      parameter(nvecsz = mvsiz)
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      "task_c.inc"
#include      "ige3d_c.inc"
C-----------------------------------------------
C   ROUTINE PURPOSE:
C   ================
C   CLASSIFIES NODES INTO BOXES
C   SEARCHES FOR BOXES CONCERNED FOR EACH FACET
C   SEARCHES FOR CANDIDATES
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C
C     NAME         DESCRIPTION                       I/O
C
C     IRECT(4,*)   FACET CONNECTIVITY ARRAY          I
C     X(3,*)       NODAL COORDINATES                 I
C     NSV          NODE NUMBERS                      I
C     XMAX         largest existing abscissa         I
C     XMAX         largest existing ordinate         I
C     XMAX         largest existing coordinate       I
C     I_STOK       storage level of candidate
C                                I/O impact couples
C     CAND_N       node result boxes
C     CAND_E       element result box addresses
C                  MULNSN = MULTIMP*NSN MAX SIZE ALLOWED NOW FOR
C                  NODE,ELEMENT CANDIDATE COUPLES
C     NOINT        USER INTERFACE NUMBER
C     TZINF        INFLUENCE ZONE SIZE
C
C     PROV_N       provisional CAND_N (static variable in i7tri)
C     PROV_E       provisional CAND_E (static variable in i7tri)
 
C     VOXEL(ix,iy,iz) contains the local number of the first node in
C                  the box
C     NEXT_NOD(i)  next node in the same box (if /= 0)
C     LAST_NOD(i)  last node in the same box (if /= 0)
C                  used only to go directly from first
C                       node to last
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER I_MEM,ESHIFT,NSN,ISZNSNR,NSNROLD,NIN,ITASK,
     .        MULNSN,NOINT,INACTI,IFQ,NSNR,IGAP,NBX,NBY,NBZ,
     .        NSV(*),CAND_N(*),CAND_E(*),CAND_A(*),IFPEN(*),RENUM(*),
     .        INTTH,IRECT(4,*), VOXEL(NBX+2,NBY+2,NBZ+2),II_STOK,
     .        KREMNOD(*),REMNOD(*),ITAB(*),FLAGREMNODE,ITIED
      INTEGER, INTENT(in) :: NRTM !< number of segments per threads
      INTEGER, INTENT(in) :: TOTAL_NB_NRTM !< total number of segments
      INTEGER, INTENT(IN) :: INTHEAT
      INTEGER, INTENT(IN) :: IDT_THERM
      INTEGER, INTENT(IN) :: NODADT_THERM
      my_real
     .   X(3,*),XYZM(12),CAND_P(*),STF(*),STFN(*),GAP_S(*),GAP_M(*),
     .   tzinf,marge,gap,gapmin,gapmax,bgapsmx,
     .   curv_max(*),gap_s_l(*),gap_m_l(*),cand_f(*)
      my_real , INTENT(IN) :: drad,dgapload
      INTEGER, INTENT(inout) :: REMOTE_S_NODE
      INTEGER, DIMENSION(NSNR), INTENT(inout) :: LIST_REMOTE_S_NODE
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,
     .        nn,ne,k,l,j_stok,jj,
     .        prov_n(mvsiz),prov_e(mvsiz),
     .        oldnum(isznsnr), nsnf, nsnl,delnod,m
      INTEGER, DIMENSION(:), ALLOCATABLE :: TAGREMNODE
      my_real
     .   xs,ys,zs,sx,sy,sz,s2,
     .   xmin, xmax,ymin, ymax,zmin, zmax,
     .   xx1,xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4,
     .   d1x,d1y,d1z,d2x,d2y,d2z,dd1,dd2,d2,a2
c temporary
      INTEGER  LAST_NOD(NSN+NSNR)
      INTEGER  IX,IY,IZ,M1,M2,M3,M4,
     .         IX1,IY1,IZ1,IX2,IY2,IZ2
      INTEGER,  DIMENSION(:),ALLOCATABLE :: IIX,IIY,IIZ
      my_real
     .   xminb,yminb,zminb,xmaxb,ymaxb,zmaxb,
     .   xmine,ymine,zmine,xmaxe,ymaxe,zmaxe,aaa
      INTEGER CPT_VOX0
      LOGICAL TEST_V0
      INTEGER FIRST,LAST,IERROR
      LOGICAL DBG_type18_fvm
      SAVE IIX,IIY,IIZ
C-----------------------------------------------
      IF(itask == 0)THEN
        remote_s_node = 0
        cpt_vox0 = 0
        ALLOCATE(next_nod(nsn+nsnr),stat=ierror)
        IF(ierror/=0) THEN
          CALL ancmsg(msgid=19,anmode=aninfo,
     .        c1='(/INTER/TYPE7)')
          CALL arret(2)
        ENDIF
        ALLOCATE(iix(nsn+nsnr),iiy(nsn+nsnr),iiz(nsn+nsnr),stat=ierror)
        IF(ierror/=0) THEN
          CALL ancmsg(msgid=19,anmode=aninfo,
     .        c1='(/INTER/TYPE7)')
          CALL arret(2)
        ENDIF
      END IF
 
      CALL my_barrier !Barrier to wait init voxel and allocation NEXT_NOD
 
! Initial phase of BPE and BPN construction moved from I7BUCE => I7TRI
 
      xmin = xyzm(1)
      ymin = xyzm(2)
      zmin = xyzm(3)
      xmax = xyzm(4)
      ymax = xyzm(5)
      zmax = xyzm(6)
 
      !dev future: xminb larger than xmin...
      IF( to_trim(nin) ) THEN
        xminb = xyzm(7)
        yminb = xyzm(8)
        zminb = xyzm(9)
        xmaxb = xyzm(10)
        ymaxb = xyzm(11)
        zmaxb = xyzm(12)
      ELSE
! we use the global bounding box
! the box reduction seemed inefficient
! at the first sorting
! because too many SECONDARY nodes were outside
! the reduced box
        xminb = xmin
        yminb = ymin
        zminb = zmin
        xmaxb = xmax
        ymaxb = ymax
        zmaxb = zmax
      ENDIF
 
      ! In SPMD, for inacti or IFQ, retrieves old numbering of non-local candidates
      IF(nspmd>1.AND.(inacti==5.OR.inacti==6.OR.inacti==7.OR.ifq>0.OR.itied/=0)) THEN
        CALL spmd_oldnumcd(renum,oldnum,isznsnr,nsnrold,intheat,idt_therm,nodadt_therm)
      END IF
 
C=======================================================================
C 1   NODE CLUSTERING
C=======================================================================
      IF(itask==0.AND.total_nb_nrtm>0)THEN
        !==========LAGRANGIAN and STAGGERED ALE==========!
        test_v0 = first_test(nin)
        DO i=1,nsn
          iix(i)=0
          iiy(i)=0
          iiz(i)=0
          IF(stfn(i) == zero)cycle
          j=nsv(i)
          !c optimization // search for nodes within xmin xmax of
          !C processor elements
          IF(x(1,j) < xmin)  cycle
          IF(x(1,j) > xmax)  cycle
          IF(x(2,j) < ymin)  cycle
          IF(x(2,j) > ymax)  cycle
          IF(x(3,j) < zmin)  cycle
          IF(x(3,j) > zmax)  cycle
          iix(i)=int(nbx*(x(1,j)-xminb)/(xmaxb-xminb))
          iiy(i)=int(nby*(x(2,j)-yminb)/(ymaxb-yminb))
          iiz(i)=int(nbz*(x(3,j)-zminb)/(zmaxb-zminb))
          iix(i)=max(1,2+min(nbx,iix(i)))
          iiy(i)=max(1,2+min(nby,iiy(i)))
          iiz(i)=max(1,2+min(nbz,iiz(i)))
          first = voxel(iix(i),iiy(i),iiz(i))
          IF(test_v0) THEN
! count the number of SECONDARY nodes outiside the reduced box
            IF(iix(i) ==  1    .OR.  iiy(i) ==  1    .OR. iiz(i) == 1 .AND.
     .         iix(i) ==  nbx+2 .OR. iiy(i) == nby+2 .OR. iiz(i) == nbz+2) THEN
              cpt_vox0 = cpt_vox0 +1
            ENDIF
          ENDIF
          IF(first == 0)THEN
            !Empty Cell
            voxel(iix(i),iiy(i),iiz(i)) = i ! first
            next_nod(i)                 = 0 ! last one
            last_nod(i)                 = 0 ! no last
          ELSEIF(last_nod(first) == 0)THEN
            !cell containing one node
            !add as next node
            next_nod(first) = i ! next
            last_nod(first) = i ! last
            next_nod(i)     = 0 ! last one
          ELSE
            !
            !jump to the last node of the cell
            last = last_nod(first) ! last node in this voxel
            next_nod(last)  = i ! next
            last_nod(first) = i ! last
            next_nod(i)     = 0 ! last one
          ENDIF
        ENDDO
 
C=======================================================================
C 2   REMOTE NODE CLUSTERING
C=======================================================================
        DO j = 1, nsnr
 
          IF(xrem(1,j) < xmin)  cycle
          IF(xrem(1,j) > xmax)  cycle
          IF(xrem(2,j) < ymin)  cycle
          IF(xrem(2,j) > ymax)  cycle
          IF(xrem(3,j) < zmin)  cycle
          IF(xrem(3,j) > zmax)  cycle
 
          remote_s_node = remote_s_node + 1
          list_remote_s_node( remote_s_node ) = j
          iix(nsn+j)=int(nbx*(xrem(1,j)-xminb)/(xmaxb-xminb))
          iiy(nsn+j)=int(nby*(xrem(2,j)-yminb)/(ymaxb-yminb))
          iiz(nsn+j)=int(nbz*(xrem(3,j)-zminb)/(zmaxb-zminb))
          iix(nsn+j)=max(1,2+min(nbx,iix(nsn+j)))
          iiy(nsn+j)=max(1,2+min(nby,iiy(nsn+j)))
          iiz(nsn+j)=max(1,2+min(nbz,iiz(nsn+j)))
 
          first = voxel(iix(nsn+j),iiy(nsn+j),iiz(nsn+j))
          IF(test_v0) THEN
! count the number of SECONDARY nodes outiside the reduced box
            IF(iix(j+nsn) ==  1     .OR. iiy(j+nsn) ==  1     .OR. iiz(j+nsn) == 1 .AND.
     .         iix(j+nsn) ==  nbx+2 .OR. iiy(j+nsn) ==  nby+2 .OR. iiz(j+nsn) == nbz+2) THEN
              cpt_vox0 = cpt_vox0 +1
            ENDIF
          ENDIF
 
          IF(first == 0)THEN
            ! Empty Cell
            voxel(iix(nsn+j),iiy(nsn+j),iiz(nsn+j)) = nsn+j ! first
            next_nod(nsn+j)     = 0 ! last one
            last_nod(nsn+j)     = 0 ! no last
          ELSEIF(last_nod(first) == 0)THEN
            ! cell containing one node, add it as next node
            next_nod(first) = nsn+j  ! next
            last_nod(first) = nsn+j  ! last
            next_nod(nsn+j)  = 0     ! last one
          ELSE
            ! , jump to the last node of the cell
            last = last_nod(first)  ! last node in this voxel
            next_nod(last)  = nsn+j ! next
            last_nod(first) = nsn+j ! last
            next_nod(nsn+j)     = 0 ! last one
          ENDIF
        ENDDO
      END IF !ITASK == 0
 
      CALL my_barrier ! Barrier to wait task0 treatment
 
      IF(first_test(nin)) THEN
!During the first sorting, we use the reduced box
!If too many (>5 %) nodes are outside the reduced box
!then the box is not reduced anymore
        CALL my_barrier
        IF(itask == 0) THEN
          IF(cpt_vox0 > 5*(remote_s_node + nsn)/100) to_trim(nin) = .false.
          first_test(nin) = .false.
        ENDIF
        CALL my_barrier
      ENDIF
 
 
C=======================================================================
C 3   FACE RECOVERY AND ENUMERATION OF CANDIDATE COUPLES
C=======================================================================
      j_stok = 0
      IF(flagremnode == 2) THEN
        ALLOCATE(tagremnode(numnod+numfakenodigeo))
        DO i=1,numnod+numfakenodigeo
          tagremnode(i) = 0
        ENDDO
      ENDIF
      DO ne=1,nrtm
        IF(stf(ne) == zero)cycle ! do not keep destroyed facets
        IF(flagremnode == 2) THEN
          k = kremnod(2*(ne-1)+1)+1
          l = kremnod(2*(ne-1)+2)
          DO i=k,l
            tagremnode(remnod(i)) = 1
          ENDDO
        ENDIF
        IF(igap == 0)THEN
          aaa = tzinf+curv_max(ne)
        ELSE
          aaa = marge+curv_max(ne)+max(min(gapmax,max(gapmin,bgapsmx+gap_m(ne)))+dgapload,drad)
        ENDIF
        !it is possible to improve the algorithm by subdividing the facet
        !into 2(4,3,6,9...) if the facet is large compared to AAA and inclined
        m1 = irect(1,ne)
        m2 = irect(2,ne)
        m3 = irect(3,ne)
        m4 = irect(4,ne)
 
        xx1=x(1,m1)
        xx2=x(1,m2)
        xx3=x(1,m3)
        xx4=x(1,m4)
        xmaxe=max(xx1,xx2,xx3,xx4)
        xmine=min(xx1,xx2,xx3,xx4)
 
        yy1=x(2,m1)
        yy2=x(2,m2)
        yy3=x(2,m3)
        yy4=x(2,m4)
        ymaxe=max(yy1,yy2,yy3,yy4)
        ymine=min(yy1,yy2,yy3,yy4)
 
        zz1=x(3,m1)
        zz2=x(3,m2)
        zz3=x(3,m3)
        zz4=x(3,m4)
        zmaxe=max(zz1,zz2,zz3,zz4)
        zmine=min(zz1,zz2,zz3,zz4)
 
        !surface calculation (for future candidate elimination)
        sx = (yy3-yy1)*(zz4-zz2) - (zz3-zz1)*(yy4-yy2)
        sy = (zz3-zz1)*(xx4-xx2) - (xx3-xx1)*(zz4-zz2)
        sz = (xx3-xx1)*(yy4-yy2) - (yy3-yy1)*(xx4-xx2)
        s2 = sx*sx + sy*sy + sz*sz
 
        !index of voxels occupied by the facet
        IF(nbx>1) THEN
          ix1=int(nbx*(xmine-aaa-xminb)/(xmaxb-xminb))
          ix2=int(nbx*(xmaxe+aaa-xminb)/(xmaxb-xminb))
        ELSE
          ix1=-2
          ix2=1
        ENDIF
 
        IF(nby>1) THEN
          iy1=int(nby*(ymine-aaa-yminb)/(ymaxb-yminb))
          iy2=int(nby*(ymaxe+aaa-yminb)/(ymaxb-yminb))
        ELSE
          iy1=-2
          iy2=1
        ENDIF
 
        IF(nbz>1) THEN
          iz1=int(nbz*(zmine-aaa-zminb)/(zmaxb-zminb))
          iz2=int(nbz*(zmaxe+aaa-zminb)/(zmaxb-zminb))
        ELSE
          iz1=-2
          iz2=1
        ENDIF
 
        ix1=max(1,2+min(nbx,ix1))
        iy1=max(1,2+min(nby,iy1))
        iz1=max(1,2+min(nbz,iz1))
 
        ix2=max(1,2+min(nbx,ix2))
        iy2=max(1,2+min(nby,iy2))
        iz2=max(1,2+min(nbz,iz2))
 
        !nbpelem = 0
        !nnpelem = 0
        !nnr0pelem = 0
        !nnrpelem = 0
 
        DO iz = iz1,iz2
          DO iy = iy1,iy2
            DO ix = ix1,ix2
 
              ! nbpelem = nbpelem + 1
 
              jj = voxel(ix,iy,iz)
              DO WHILE(jj /= 0)
                delnod = 0
 
                !nnpelem = nnpelem + 1
                IF(jj<=nsn)THEN
                  nn=nsv(jj)
 
                  IF(nn == m1)GOTO 200
                  IF(nn == m2)GOTO 200
                  IF(nn == m3)GOTO 200
                  IF(nn == m4)GOTO 200
 
                  IF(flagremnode == 2) THEN
                    IF( tagremnode(nsv(jj)) == 1)GOTO 200
                  ENDIF
                  xs = x(1,nn)
                  ys = x(2,nn)
                  zs = x(3,nn)
                  IF(igap /= 0)THEN
                    aaa = marge+curv_max(ne)+max(min(gapmax,max(gapmin,gap_s(jj)+gap_m(ne)))+dgapload,drad)
                  ENDIF
                ELSE
                  j=jj-nsn
                  IF(flagremnode == 2) THEN
                    nn = irem(1,j)
                    k = kremnod(2*(ne-1)+2) + 1
                    l = kremnod(2*(ne-1)+3)
                    DO m=k,l
                      IF(remnod(m) == -irem(2,j) ) THEN
                        delnod = delnod + 1
                        EXIT
                      ENDIF
                    ENDDO
                    IF(delnod /= 0)GOTO 200
                  ENDIF
 
                  xs = xrem(1,j)
                  ys = xrem(2,j)
                  zs = xrem(3,j)
                  IF(igap /= 0)THEN
                    aaa = marge+curv_max(ne)+max(min(gapmax,max(gapmin,xrem(9,j)+gap_m(ne)))+dgapload,drad)
                  ENDIF
                ENDIF
 
                IF(xs<=xmine-aaa)GOTO 200
                IF(xs>=xmaxe+aaa)GOTO 200
                IF(ys<=ymine-aaa)GOTO 200
                IF(ys>=ymaxe+aaa)GOTO 200
                IF(zs<=zmine-aaa)GOTO 200
                IF(zs>=zmaxe+aaa)GOTO 200
 
                !underestimation of distance**2 for candidate elimination
                !nnr0pelem = nnr0pelem + 1
 
                d1x = xs - xx1
                d1y = ys - yy1
                d1z = zs - zz1
                d2x = xs - xx2
                d2y = ys - yy2
                d2z = zs - zz2
                dd1 = d1x*sx+d1y*sy+d1z*sz
                dd2 = d2x*sx+d2y*sy+d2z*sz
                IF(dd1*dd2 > zero)THEN
                  d2 = min(dd1*dd1,dd2*dd2)
                  a2 = aaa*aaa*s2
                  IF(d2 > a2)GOTO 200
                ENDIF
 
                !nnrpelem = nnrpelem + 1
 
                j_stok = j_stok + 1
                prov_n(j_stok) = jj
                prov_e(j_stok) = ne
                IF(j_stok == nvsiz)THEN
                  CALL i7sto(
     1               nvsiz ,irect  ,x     ,nsv   ,ii_stok,
     2               cand_n,cand_e ,mulnsn,noint ,marge  ,
     3               i_mem ,prov_n ,prov_e,eshift,inacti ,
     4               ifq   ,cand_a ,cand_p,ifpen ,nsn    ,
     5               oldnum,nsnrold,igap  ,gap   ,gap_s  ,
     6               gap_m ,gapmin ,gapmax,curv_max,nin  ,
     7               gap_s_l,gap_m_l,intth,drad,itied    ,
     8               cand_f ,dgapload)
                  IF(i_mem==2)GOTO 100
                  j_stok = 0
                ENDIF
  200           CONTINUE
                jj = next_nod(jj)
              ENDDO ! WHILE(JJ /= 0)
            ENDDO ! X
          ENDDO  ! Y
        ENDDO   ! Z
        IF(flagremnode == 2) THEN
          k = kremnod(2*(ne-1)+1)+1
          l = kremnod(2*(ne-1)+2)
          DO i=k,l
            tagremnode(remnod(i)) = 0
          ENDDO
        ENDIF
      ENDDO
C=======================================================================
C 4   STORAGE
C=======================================================================
      IF(j_stok/=0)CALL i7sto(
     1              j_stok,irect  ,x     ,nsv   ,ii_stok,
     2              cand_n,cand_e ,mulnsn,noint ,marge  ,
     3              i_mem ,prov_n ,prov_e,eshift,inacti ,
     4              ifq   ,cand_a ,cand_p,ifpen ,nsn    ,
     5              oldnum,nsnrold,igap  ,gap   ,gap_s  ,
     6              gap_m ,gapmin ,gapmax,curv_max,nin  ,
     7              gap_s_l,gap_m_l,intth,drad  ,itied    ,
     8              cand_f ,dgapload)
 
C=======================================================================
C 5   VOXEL CLUSTERING  RESET
C=======================================================================
  100 CONTINUE
      ! Barrier to avoid reinitialization before end of sorting
      CALL my_barrier
      IF(total_nb_nrtm>0) THEN
        nsnf = 1 + itask*nsn / nthread
        nsnl = (itask+1)*nsn / nthread
        DO i=nsnf,nsnl
            IF(iix(i)/=0)THEN
              voxel(iix(i),iiy(i),iiz(i))=0
            ENDIF
        ENDDO
        nsnf = 1 + itask*remote_s_node / nthread
        nsnl = (itask+1)*remote_s_node / nthread
        IF(itask+1==nthread) nsnl=remote_s_node
        DO jj = nsnf, nsnl
            j = list_remote_s_node(jj)
            voxel(iix(nsn+j),iiy(nsn+j),iiz(nsn+j))=0
        ENDDO
      ENDIF
      CALL my_barrier()
 
C=======================================================================
C 6   DEBUG OUTPUT
C=======================================================================
 
      dbg_type18_fvm=.false.
      if(inacti==7 .AND. dbg_type18_fvm)then
        write(*,fmt='(A)')"------------------------------------------"
        write(*,*)"RESULT : Search Algorithm with VOXEL partitioning"
        write(*,*)"  Number of couples =", ii_stok
        if(ii_stok>0)then
          write(*,fmt='(A,(I10))')"  --> SECONDARY Node ids: ", cand_n(1:ii_stok)
          write(*,fmt='(A,(I10))')"  --> Local Face ids: ", cand_e(1:ii_stok)
        endif
        write(*,*)" Structure domain :"
        write(*,fmt='(A,F30.16,A,F30.16)')" Xmin=",xmin,"  Xmax=",xmax
        write(*,fmt='(A,F30.16,A,F30.16)')" Ymin=",ymin,"  Ymax=",ymax
        write(*,fmt='(A,F30.16,A,F30.16)')" Zmin=",zmin,"  Zmax=",zmax
        write(*,*)" Partitioning domain :"
        write(*,*)"   TZINF,AAA=",tzinf,aaa
        write(*,fmt='(A,F30.16,A,F30.16)')" Xmin=",xmin-aaa,"  Xmax=",xmax+aaa
        write(*,fmt='(A,F30.16,A,F30.16)')" Ymin=",ymin-aaa,"  Ymax=",ymax+aaa
        write(*,fmt='(A,F30.16,A,F30.16)')" Zmin=",zmin-aaa,"  Zmax=",zmax+aaa
        write(*,fmt='(a)')"------------------------------------------"
      endif
 
C=======================================================================
C 7   DEALLOCATE
C=======================================================================
      IF(ITASK == 0)THEN
        DEALLOCATE(NEXT_NOD)
        DEALLOCATE(IIX)
        DEALLOCATE(IIY)
        DEALLOCATE(IIZ)
      ENDIF
      IF(FLAGREMNODE == 2) THEN
        IF(ALLOCATED(TAGREMNODE)) DEALLOCATE(TAGREMNODE)
      ENDIF
 
 
      RETURN
      END