i24trivox.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/.
!||====================================================================
!||    i24trivox       ../engine/source/interfaces/intsort/i24trivox.F
!||--- called by ------------------------------------------------------
!||    i24buce         ../engine/source/interfaces/intsort/i24buce.F
!||--- calls      -----------------------------------------------------
!||    i24fic_getn     ../engine/source/interfaces/int24/i24for3e.F
!||    i24sto          ../engine/source/interfaces/intsort/i24sto.F
!||    my_barrier      ../engine/source/system/machine.F
!||    spmd_oldnumcd   ../engine/source/mpi/interfaces/spmd_i7tool.F
!||--- uses       -----------------------------------------------------
!||    tri7box         ../engine/share/modules/tri7box.F
!||====================================================================
      SUBROUTINE i24trivox(
     1      NSN    ,NSNR    ,ISZNSNR ,I_MEM   ,VMAXDT   ,
     2      IRECT  ,X       ,STF     ,STFN    ,XYZM     ,
     3      NSV    ,II_STOK ,CAND_N  ,ESHIFT  ,CAND_E   ,
     4      MULNSN ,NOINT   ,V       ,BGAPSMX  ,
     5      VOXEL  ,NBX     ,NBY     ,NBZ     ,PMAX_GAP ,
     6      NRTM   ,GAP_S   ,GAP_M   ,MARGE   ,CURV_MAX ,
     7      NIN    ,ITASK   ,PENE_OLD,ITAB    ,NBINFLG  ,
     8      MBINFLG,ILEV    ,MSEGTYP ,EDGE_L2 ,IEDGE    ,
     9      ISEADD ,ISEDGE  ,CAND_T  ,FLAGREMNODE,KREMNOD,
     A      REMNOD ,CAND_A  ,RENUM   ,NSNROLD ,IRTSE    ,
     B      IS2SE  ,NSNE    ,DGAPLOAD,INTHEAT,IDT_THERM,NODADT_THERM)
C============================================================================
C   M o d u l e s
C-----------------------------------------------
      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"
C-----------------------------------------------
C   role of the routine:
C   ===================
C   classify nodes into boxes
C   search for each facet the concerned boxes
C   search for candidates
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C
C     NOM          DESCRIPTION                       E/S
C
C     IRECT(4,*)   ARRAY OF CONEC FACETTES        E
C     X(3,*)       COORDONNEES NODALES               E
C     NSV          NOS SYSTEMES DES NODES E
C     Xmax larger abcisse existing e
C     XMAX largest order.existing E
C     Xmax larger existing side E
C     I_STOK       storage level of pairs
C                                CANDIDATES impact    E/S
C     CAND_N       boites resultats nodes C     CAND_E       adresses des boites resultat elements
C                  MULNSN = MULTIMP*NSN maximum size now allowed for
C                  COUPLES NODES,ELT CANDIDATES
C     NOINT        INTERFACE USER NUMBER
C
C     Prov_n Provisional Cand_n (static variable in i7tri)
C     PROV_E       CAND_E provisoire (variable static in i7tri)
 
C     VOXEL(ix,iy,iz) contains the local number of the first node of
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 the first
C                       node at the last
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER I_MEM,ESHIFT,NSN,ISZNSNR,NRTM,NIN,ITASK,
     .        MULNSN,NOINT,NSNR,NBX,NBY,NBZ,IEDGE,NSNE,
     .        NSV(*),CAND_N(*),CAND_E(*),
     .        IRECT(4,*), VOXEL(NBX+2,NBY+2,NBZ+2),II_STOK,ITAB(*),
     .        NBINFLG(*),MBINFLG(*),ILEV,MSEGTYP(*),CAND_T(*),
     .        ISEADD(*) ,ISEDGE(*),FLAGREMNODE,KREMNOD(*),REMNOD(*),CAND_A(*),
     .        RENUM(*),NSNROLD,IRTSE(5,*),IS2SE(2,*)
      INTEGER, INTENT(IN) :: INTHEAT
      INTEGER, INTENT(IN) :: IDT_THERM
      INTEGER, INTENT(IN) :: NODADT_THERM
C     REAL
      my_real
     .   x(3,*),v(3,*),xyzm(6),stf(*),stfn(*),gap_s(*),
     .   gap_m(*),curv_max(*),pene_old(5,nsn),edge_l2(*),
     .   marge,bgapsmx,pmax_gap,vmaxdt
      my_real , INTENT(IN) :: dgapload
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,m,nse,ns
C     REAL
      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 provisional
      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 FIRST,LAST
      SAVE IIX,IIY,IIZ
      INTEGER, DIMENSION(NUMNOD+NSNE) :: TAG
C --------------------------------
C TYPE24 E2E - I24FIC_GETN method
C --------------------------------
      INTEGER IK1(4),IK2(4),IED,NS1,NS2,NS1ID,NS2ID
      DATA ik1 /1,2,3,4/
      DATA ik2 /2,3,4,1/
C-----------------------------------------------
      IF(itask == 0)THEN
        ALLOCATE(next_nod(nsn+nsnr))
        ALLOCATE(iix(nsn+nsnr))
        ALLOCATE(iiy(nsn+nsnr))
        ALLOCATE(iiz(nsn+nsnr))
      END IF
C Barrier to wait init voxel and allocation NEX_NOD
      CALL my_barrier
C initial construction phase of BPE and BPN moved from I7BUCE => I7TRI
C
      xmin = xyzm(1)
      ymin = xyzm(2)
      zmin = xyzm(3)
      xmax = xyzm(4)
      ymax = xyzm(5)
      zmax = xyzm(6)
 
c     Dev Future: xminb larger than Xmin ...
      xminb = xmin
      yminb = ymin
      zminb = zmin
      xmaxb = xmax
      ymaxb = ymax
      zmaxb = zmax
 
c!!!!!!!!!!!!!!! A VERIFIER !!!!!!!!!!!!!!!
C In SPMD, for IFQ, finds former number of non -local candidates
c useless for INT 24 !!!!!!!!!!!!!!!!!
      IF(nspmd>1) THEN
        CALL spmd_oldnumcd(renum,oldnum,isznsnr,nsnrold,intheat,idt_therm,nodadt_therm)
      END IF
 
C=======================================================================
C 1   putting nodes in the boxes
C=======================================================================
C Note for Edge2Edge : X is no more the Radioss X Array but an extension
C                      NUMNOD+SNE
C                      It is updated at any cycle
      IF(itask == 0)THEN
        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 included in 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(first == 0)THEN
c         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
c         cell containing one node
c         add as next node
            next_nod(first) = i ! next
            last_nod(first) = i ! last
            next_nod(i)     = 0 ! last one
          ELSE
c
c         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   putting nodes in the boxes
C     non -local candidates in SPMD
C=======================================================================
        DO j = 1, nsnr
 
          IF(irem(8,j)==-1) cycle    ! case IREM / ISEDGE_FI==-1  : Node was added due to Fictive Remote Node only
          ! Do not retain in sorting, otherwise node can be candidate twice
 
          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(first == 0)THEN
c         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
c         cell containing one node
c         add as next node
            next_nod(first) = nsn+j  ! next
            last_nod(first) = nsn+j  ! last
            next_nod(nsn+j)  = 0     ! last one
          ELSE
c
c         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
C Barrier to wait task0 treatment
      CALL my_barrier
C=======================================================================
C 3 Searching boxes concerning each facet
C     and creation of candidates
C=======================================================================
      j_stok = 0
      IF(flagremnode == 2)THEN
        DO i=1,numnod+nsne
          tag(i) = 0
        ENDDO
      END IF
 
      DO ne=1,nrtm
C We do not retain the Destruit facets
        IF(stf(ne) == zero)cycle
 
        aaa = marge+curv_max(ne)+bgapsmx+pmax_gap+vmaxdt
     +      + gap_m(ne)+dgapload
 
 
c     It is possible to improve the algo by cutting the facet
c     in 2 (4,3,6,9 ...) if the facet is large in front of AAA and inclinee
 
        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)
 
 
c        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
 
c        index of voxels occupied by the facet
 
        ix1=int(nbx*(xmine-aaa-xminb)/(xmaxb-xminb))
        iy1=int(nby*(ymine-aaa-yminb)/(ymaxb-yminb))
        iz1=int(nbz*(zmine-aaa-zminb)/(zmaxb-zminb))
 
        ix1=max(1,2+min(nbx,ix1))
        iy1=max(1,2+min(nby,iy1))
        iz1=max(1,2+min(nbz,iz1))
 
        ix2=int(nbx*(xmaxe+aaa-xminb)/(xmaxb-xminb))
        iy2=int(nby*(ymaxe+aaa-yminb)/(ymaxb-yminb))
        iz2=int(nbz*(zmaxe+aaa-zminb)/(zmaxb-zminb))
 
        ix2=max(1,2+min(nbx,ix2))
        iy2=max(1,2+min(nby,iy2))
        iz2=max(1,2+min(nbz,iz2))
 
        IF(flagremnode == 2)THEN
          k = kremnod(2*(ne-1)+1)+1
          l = kremnod(2*(ne-1)+2)
          DO i=k,l
            tag(remnod(i)) = 1
          ENDDO
        END if!(FLAGREMNODE == 2)THEN
cc         nbpelem = 0
cc'nnpelem = 0
cc'nnr0pelem = 0
cc'nnrpelem = 0
 
        DO iz = iz1,iz2
          DO iy = iy1,iy2
            DO ix = ix1,ix2
 
cc             nbpelem = nbpelem + 1
 
              jj = voxel(ix,iy,iz)
 
              DO WHILE(jj /= 0)
 
cc             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(tag(nn) == 1)GOTO 200
                  END IF
C----- fictitious nodes on edges-: auto-empact excluded -------
                  IF (nn >numnod) THEN
                    ns = nn-numnod
                    CALL i24fic_getn(ns      ,irtse   ,is2se   ,nse    ,
     +                               ns1     ,ns2     )
                    IF(ns1 == m1 .OR. ns2 == m1) GOTO 200
                    IF(ns1 == m2 .OR. ns2 == m2) GOTO 200
                    IF(ns1 == m3 .OR. ns2 == m3) GOTO 200
                    IF(ns1 == m4 .OR. ns2 == m4) GOTO 200
                  END IF
                  xs = x(1,nn)
                  ys = x(2,nn)
                  zs = x(3,nn)
c PMAX_GAP is a global overestimate penetration
c NEED to communicate in SPMD
c VMAXDT is a local overestimate of relative incremental displacement
c NO need to communicate in SPMD
 
                  IF (iedge > 0) THEN
                    aaa = marge + curv_max(ne)
     +                  + max(gap_s(jj)+gap_m(ne)+edge_l2(jj)+dgapload
     +                       ,pene_old(3,jj))+vmaxdt
                  ELSE
                    aaa = marge + curv_max(ne)
     +                  + max(gap_s(jj)+gap_m(ne)+dgapload
     +                       ,pene_old(3,jj))+vmaxdt
                  END IF
                ELSE
                  j=jj-nsn
                  IF(flagremnode == 2)THEN
                    k = kremnod(2*(ne-1)+2) + 1
                    l = kremnod(2*(ne-1)+3)
                    IF(irem(8,j)==1) THEN
                      DO m=k,l
                        IF(remnod(m) == -irem(2,j) ) GOTO 200
                      ENDDO
                    ELSE
                      DO m=k,l
                        IF(remnod(m) == -irem(2,j) ) GOTO 200
                      ENDDO
                    ENDIF
                  END if!(FLAGREMNODE == 2)THEN
C
C Auto impact between main surface and second Edge
C can happen with Remote nodes when Second nodes are on border between 2 domains
C and Fictive node is remote
                  IF(irem(8,j)==1) THEN
                    ! Same than in I24FIC_GETN but for Remote Node
                    i24irempnsne=irem(7,j)                     ! in IREM IRTSE is located in IREM(I24IREMPNSNE,J) to IREM(I24IREMPNSNE+4,J)
                    ied = irem(i24irempnsne+4,j)               ! IED = IRTSE(5,xx)
                    ns1 = irem(i24irempnsne-1+ik1(ied),j)      ! NS1 = IRTSE(IK1(IED))
                    ns2 = irem(i24irempnsne-1+ik2(ied),j)      ! NS2 = IRTSE(IK2(IED))
                    ns1id = irem(2,ns1)                        ! ITAB Remote NS1
                    ns2id = irem(2,ns2)                        ! ITAB Remote NS2
                    IF (ns1id == itab(m1) .OR. ns2id == itab(m1)) GOTO 200
                    IF (ns1id == itab(m2) .OR. ns2id == itab(m2)) GOTO 200
                    IF (ns1id == itab(m3) .OR. ns2id == itab(m3)) GOTO 200
                    IF (ns1id == itab(m4) .OR. ns2id == itab(m4)) GOTO 200
                  ENDIF
                  xs = xrem(1,j)
                  ys = xrem(2,j)
                  zs = xrem(3,j)
                  aaa = marge+curv_max(ne)
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c                                      +EDGE_L2(JJ) remote
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     +            + max(xrem(igapxremp,j)+gap_m(ne)+dgapload,xrem(i24xremp+6,j))
     +            + vmaxdt
                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
 
c    underestimation of distance**2 for candidate elimination
 
cc'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
 
cc             nnrpelem = nnrpelem + 1
 
                j_stok = j_stok + 1
                prov_n(j_stok) = jj
                prov_e(j_stok) = ne
                IF(j_stok == nvsiz)THEN
 
                  CALL i24sto(
     1               nvsiz ,irect  ,x     ,nsv   ,ii_stok,
     2               cand_n,cand_e ,mulnsn,noint ,marge  ,
     3               i_mem ,prov_n ,prov_e,eshift,v      ,
     4               nsn   ,gap_s  ,gap_m ,curv_max,nin  ,
     5               pene_old,nbinflg  ,mbinflg,ilev,msegtyp,
     6               edge_l2,iedge,iseadd ,isedge ,cand_t,itab,
     7               cand_a,oldnum,nsnrold,dgapload)
                  IF(i_mem==2)GOTO 100
                  j_stok = 0
                ENDIF
 
  200           CONTINUE
 
                jj = next_nod(jj)
 
              ENDDO ! WHILE(JJ /= 0)
 
            ENDDO
          ENDDO
        ENDDO
cc             nbpelg = nbpelg + nbpelem
cc             nnpelg = nnpelg + nnpelem
cc             nnrpelg = nnrpelg + nnrpelem
cc             nnr0pelg = nnr0pelg + nnr0pelem
        IF(flagremnode == 2)THEN
          k = kremnod(2*(ne-1)+1)+1
          l = kremnod(2*(ne-1)+2)
          DO i=k,l
            tag(remnod(i)) = 0
          ENDDO
        END IF
      ENDDO
 
C-------------------------------------------------------------------------
C     end of sorting
C-------------------------------------------------------------------------
      IF(j_stok/=0)CALL i24sto(
     1              j_stok,irect  ,x     ,nsv   ,ii_stok,
     2              cand_n,cand_e ,mulnsn,noint ,marge  ,
     3              i_mem ,prov_n ,prov_e,eshift,v      ,
     4              nsn   ,gap_s  ,gap_m ,curv_max,nin  ,
     5              pene_old,nbinflg,mbinflg,ilev ,msegtyp,
     6              edge_l2,iedge,iseadd ,isedge ,cand_t,itab,
     7              cand_a,oldnum,nsnrold,dgapload)
 
C=======================================================================
C 4   reset nodes to zero in the boxes
C=======================================================================
  100 CONTINUE
 
C Barrier to avoid reinitialization before end of sorting
      CALL my_barrier
      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
C=======================================================================
C 5   reset nodes to zero in the boxes
C     non -local candidates in SPMD
C=======================================================================
      nsnf = 1 + itask*nsnr / nthread
      nsnl = (itask+1)*nsnr / nthread
      DO j = nsnf, nsnl
        IF(irem(8,j)==-1)cycle
        voxel(iix(nsn+j),iiy(nsn+j),iiz(nsn+j))=0
      ENDDO
 
C
      CALL my_barrier()
      IF(itask == 0)THEN
        DEALLOCATE(next_nod)
        DEALLOCATE(iix)
        DEALLOCATE(iiy)
        DEALLOCATE(iiz)
      ENDIF
 
      RETURN
      END