sptrivox.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/.
 
!||====================================================================
!||    sptrivox   ../starter/source/elements/sph/sptrivox.F
!||--- called by ------------------------------------------------------
!||    spbuc31    ../starter/source/elements/sph/spbuc31.F
!||--- calls      -----------------------------------------------------
!||    sppro31    ../starter/source/elements/sph/spbuc31.F
!||--- uses       -----------------------------------------------------
!||    tri7box    ../starter/share/modules1/tri7box.F
!||====================================================================
      SUBROUTINE sptrivox(
     1      NSN    ,X       ,BMINMA  ,NOD2SP  ,
     2      NBX    ,NBY     ,NBZ     ,
     3      NLIST  ,SPBUF   ,JVOIS   ,JSTOR   ,JPERM   ,
     4      DVOIS  ,IREDUCE ,NSPHACTF,NSPHACTL,VOXEL   ,
     5      KXSP   ,IXSP    ,KREDUCE ,IPARTSP ,SZ_INTP_DIST,
     6      MAX_INTP_DIST_PART,PRE_SEARCH)
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-----------------------------------------------
      INTEGER NSPHACTF, NSPHACTL
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "sphcom.inc"
C-----------------------------------------------
C   PRE_SEARCH = 0 -> full search of neigbours 
C   PRE_SEARCH = 1 -> pre-search of neigbours for computation of max interparticle dist 
C-----------------------------------------------
C   ROLE DE LA ROUTINE:
C   ===================
C   CLASSE LES NOEUDS DANS DES BOITES
C   RECHERCHE POUR CHAQUE FACETTE DES BOITES CONCERNES
C   RECHERCHE DES CANDIDATS
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C
C     NOM          DESCRIPTION                       E/S
C
C     X(3,*)       COORDONNEES NODALES               E
C     XMAX         plus grande abcisse existante     E
C     XMAX         plus grande ordonn. existante     E
C     XMAX         plus grande cote    existante     E
C     VOXEL(ix,iy,iz) contient le numero local du premier noeud de
C                  la boite
C     NEXT_NOD(i)  noeud suivant dans la meme boite (si /= 0)
C     LAST_NOD(i)  dernier noeud dans la meme boite (si /= 0)
C                  utilise uniquement pour aller directement du premier
C                       noeud au dernier
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NSN,NBX,NBY,NBZ,
     .        NLIST(*),NOD2SP(*) ,
     .        VOXEL(NBX+2,NBY+2,NBZ+2),JVOIS(*)  ,JSTOR(*), JPERM(*) ,
     .        IREDUCE,KXSP(NISP,*), IXSP(KVOISPH,*), KREDUCE(*)
      INTEGER ,INTENT(IN) :: IPARTSP(NUMSPH),PRE_SEARCH,SZ_INTP_DIST
C     REAL
      my_real
     .   x(3,*),bminma(6),
     .   spbuf(nspbuf,*),dvois(*)
      my_real ,INTENT(INOUT) :: max_intp_dist_part(sz_intp_dist)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NB_NCN,NB_NCN1,NB_ECN,I,J,DIR,NB_NC,NB_EC,
     .        N1,N2,N3,N4,NN,NE,K,L,II,JJ,JS,NS,N,
     .        nsnf, nsnl,nvois, ig, il,nvoimax
C     REAL
      my_real
     .   dx,dy,dz,xs,ys,zs,xx,sx,sy,sz,s2,xn,yn,zn,
     .   xmin, xmax,ymin, ymax,zmin, zmax, tz,
     .   d1x,d1y,d1z,d2,a2,alpha_marge,distmax
c provisoire
      INTEGER  LAST_NOD(NSN)
      INTEGER  IX,IY,IZ,NEXT,
     .         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,bbb,
     .   aaa2,min_dist
      INTEGER FIRST,NEW,LAST
      INTEGER IPART_I,IPART_JS
      SAVE iix,iiy,iiz,distmax
C-----------------------------------------------
C      IF(ITASK == 0)THEN
      ALLOCATE(next_nod(nsn))
      ALLOCATE(iix(nsn))
      ALLOCATE(iiy(nsn))
      ALLOCATE(iiz(nsn))
C      END IF
C Barrier to wait init voxel and allocation NEX_NOD
C      CALL MY_BARRIER
C Phase initiale de construction de BPE et BPN deplacee de I7BUCE => I7TRI
C
      alpha_marge = sqrt(one +spasort)
 
      xmax = bminma(1)
      ymax = bminma(2)
      zmax = bminma(3)
      xmin = bminma(4)
      ymin = bminma(5)
      zmin = bminma(6)
 
c     dev future: xminb plus grand que xmin...
      xminb = xmin
      yminb = ymin
      zminb = zmin
      xmaxb = xmax
      ymaxb = ymax
      zmaxb = zmax
 
C=======================================================================
C 1   mise des noeuds dans les boites
C=======================================================================
C     IF(ITASK == 0)THEN
 
      distmax = zero
 
      DO i=1,nsn
        iix(i)=0
        iiy(i)=0
        iiz(i)=0
 
        j=nlist(i)
        nn = kxsp(3,j)
 
        distmax = max(distmax,spbuf(1,j))
C Optimisation // recherche les noeuds compris dans xmin xmax des
C elements du processeur
        IF(x(1,nn) < xmin)  cycle
        IF(x(1,nn) > xmax)  cycle
        IF(x(2,nn) < ymin)  cycle
        IF(x(2,nn) > ymax)  cycle
        IF(x(3,nn) < zmin)  cycle
        IF(x(3,nn) > zmax)  cycle
 
        iix(i)=int(nbx*(x(1,nn)-xminb)/(xmaxb-xminb))
        iiy(i)=int(nby*(x(2,nn)-yminb)/(ymaxb-yminb))
        iiz(i)=int(nbz*(x(3,nn)-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     END IF
C Barrier to wait task0 treatment
C     CALL MY_BARRIER
C=======================================================================
C 3   recherche des boites concernant chaque facette
C     et creation des candidats
C=======================================================================
      nvoimax = 0
      DO ne = nsphactf,nsphactl
C on ne retient pas les facettes detruites
c        IF(STF(NE) == ZERO)CYCLE
 
        j=nlist(ne)
        nn = kxsp(3,j)
c        
        min_dist = ep30
        ipart_i=ipartsp(j)
 
c a revoir !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        aaa = (spbuf(1,j)+distmax)* alpha_marge
 
c        indice des voxels occupes par la facette
 
        ix1=int(nbx*(x(1,nn)-aaa-xminb)/(xmaxb-xminb))
        iy1=int(nby*(x(2,nn)-aaa-yminb)/(ymaxb-yminb))
        iz1=int(nbz*(x(3,nn)-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*(x(1,nn)+aaa-xminb)/(xmaxb-xminb))
        iy2=int(nby*(x(2,nn)+aaa-yminb)/(ymaxb-yminb))
        iz2=int(nbz*(x(3,nn)+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))
 
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
 
                js=nlist(jj)
                ns=kxsp(3,js)
                ipart_js=ipartsp(js)
                IF(jj == ne)GOTO 200
                xs = x(1,ns)
                ys = x(2,ns)
                zs = x(3,ns)
 
                aaa = spbuf(1,j)+spbuf(1,js)
 
                bbb = aaa * alpha_marge
 
                IF(xs<=x(1,nn)-bbb)GOTO 200
                IF(xs>=x(1,nn)+bbb)GOTO 200
                IF(ys<=x(2,nn)-bbb)GOTO 200
                IF(ys>=x(2,nn)+bbb)GOTO 200
                IF(zs<=x(3,nn)-bbb)GOTO 200
                IF(zs>=x(3,nn)+bbb)GOTO 200
 
cc             nnr0pelem = nnr0pelem + 1
 
                d1x = xs - x(1,nn)
                d1y = ys - x(2,nn)
                d1z = zs - x(3,nn)
                d2 = d1x*d1x+d1y*d1y+d1z*d1z
                a2 = bbb*bbb
                IF(js==j.or.d2 > a2)GOTO 200
 
cc             nnrpelem = nnrpelem + 1
                aaa2 = aaa*aaa
                nvois=kxsp(5,j)+1
                jvois(nvois)=js
                dvois(nvois)=d2/aaa2
 
cc              Distance to closest neigbhour of the same part 
                IF (ipart_i==ipart_js) min_dist = min(min_dist,sqrt(d2))
 
                kxsp(5,j)  =nvois
 
  200           CONTINUE
 
                jj = next_nod(jj)
 
              ENDDO ! WHILE(JJ /= 0)
 
            ENDDO
          ENDDO
        ENDDO
 
        IF (pre_search==0) THEN 
cc        nbpelg = nbpelg + nbpelem
cc        nnpelg = nnpelg + nnpelem
cc        nnrpelg = nnrpelg + nnrpelem
cc        nnr0pelg = nnr0pelg + nnr0pelem
          CALL sppro31(j    ,kxsp  ,ixsp ,nod2sp,jvois,
     .               jstor,jperm ,dvois,ireduce,kreduce)
        ELSE
cc        Pre-search - computation of max interpaticle distance per part (maximum of dist to closest neigbhour)
          max_intp_dist_part(ipart_i) = max(max_intp_dist_part(ipart_i),min_dist)
          nvoimax = max(nvoimax,nvois)
        ENDIF
 
      ENDDO
 
C-------------------------------------------------------------------------
C     FIN DU TRI
C-------------------------------------------------------------------------
C=======================================================================
C 4   remise a zero des noeuds dans les boites
C=======================================================================
  100 CONTINUE
 
C Barrier to avoid reinitialization before end of sorting
C     CALL MY_BARRIER
 
      DO i=nsphactf,nsphactl
        IF(iix(i)/=0)THEN
          voxel(iix(i),iiy(i),iiz(i))=0
        ENDIF
      ENDDO
C=======================================================================
C
C     CALL MY_BARRIER()
C     IF(ITASK == 0)THEN
      DEALLOCATE(next_nod)
      DEALLOCATE(iix)
      DEALLOCATE(iiy)
      DEALLOCATE(iiz)
C     ENDIF
 
      RETURN
      END