spclasv.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/.
 
!||====================================================================
!||    spclasv   ../starter/source/elements/sph/spclasv.F
!||--- called by ------------------------------------------------------
!||    sptri     ../starter/source/elements/sph/sptri.F
!||--- calls      -----------------------------------------------------
!||====================================================================
      SUBROUTINE spclasv(X      ,SPBUF ,KXSP   ,IXSP   ,NOD2SP    ,
     1                   WASPACT,MYSPATRUE,IREDUCE,KREDUCE)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
C----6---------------------------------------------------------------7---------8
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      "sphcom.inc"
C-----------------------------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),
     .        waspact(*), ireduce, kreduce(*)
      my_real
     .   x(3,*),spbuf(nspbuf,*), myspatrue
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      integer
     .   n,inod,jnod,j,nvois,m,ncand,k1,k2,nvois1,nvois2,
     .   nvoiss,nvoiss1,nvoiss2, iaux, ierror,
     .   k, l, jk, nc, js, ns, nn, nb,jj1,jj2, jj, jjj,
     .   mwa(2*kvoisph),jstor(kvoisph), jperm(kvoisph),
     .   lvoired, ig
      my_real
     .       dms,dms2,dk,
     .       xi,yi,zi,di,xj,yj,zj,dj,dd,dvois(kvoisph),
     .       dwa(numsph)
      SAVE lvoired
C-----------------------------------------------
C       force sorting for domain decomposition
        DO ns=1,nsphact
         n=waspact(ns)
         inod=kxsp(3,n)
         xi=x(1,inod)
         yi=x(2,inod)
         zi=x(3,inod)
         di=spbuf(1,n)
         nvois=kxsp(5,n)
         DO k=1,nvois
           jnod = ixsp(k,n)
           m =nod2sp(jnod)
           xj=x(1,jnod)
           yj=x(2,jnod)
           zj=x(3,jnod)
           dj=spbuf(1,m)
           dms =di+dj
           dms2=dms*dms
           dvois(k)=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
           dvois(k)=dvois(k)/dms2
         END DO
C
         CALL myqsort(nvois,dvois,jperm,ierror)
         DO k=1,nvois
           jstor(k) = ixsp(k,n)
         END DO
C
         k1=0
         DO k=1,nvois
           jk=jstor(jperm(k))
           k1=k1+1
           ixsp(k1,n) = jk
         END DO
C
C--- Closest neighbor is stored for tensile instability treatment --
C
         IF (nspbuf==15) THEN
           spbuf(15,n)=two*sqrt(dvois(1))
         ENDIF
C
        END DO
C-----------------------------------------------
        lvoired = 0
        IF(ireduce==0)GO TO 100
C-------------------------------------------
C       tri voisins / ne garder que LVOISPH voisins effectifs
C
        DO ns=1,nsphact
         n=waspact(ns)
         dwa(n)=one
         nvois1 =kxsp(4,n)
         IF(kreduce(n)/=0.OR.nvois1>lvoisph)THEN
C
           IF(nvois1>lvoisph)THEN
             kreduce(n)=kreduce(n)+10
             lvoired = 1
           END IF
C
           inod=kxsp(3,n)
           xi=x(1,inod)
           yi=x(2,inod)
           zi=x(3,inod)
           di=spbuf(1,n)
           nvois=kxsp(5,n)
           DO k=1,nvois
             jnod = ixsp(k,n)
             m =nod2sp(jnod)
             xj=x(1,jnod)
             yj=x(2,jnod)
             zj=x(3,jnod)
             dj=spbuf(1,m)
             dms =di+dj
             dms2=dms*dms
             dvois(k)=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
             dvois(k)=dvois(k)/dms2
           END DO
C
           CALL myqsort(nvois,dvois,jperm,ierror)
           DO k=1,nvois
             jstor(k) = ixsp(k,n)
           END DO
C
           IF(kreduce(n) >= 10)dwa(n)=sqrt(dvois(lvoisph))
C
           k1=0
           DO k=1,nvois
             jk=jstor(jperm(k))
             k1=k1+1
             ixsp(k1,n) = jk
           END DO
C
         END IF
        END DO
C-------------------------------------------
C adapte diametre (reduction only) 
C
        IF(lvoired /= 0)THEN
C
          DO ns=1,nsphact
           n=waspact(ns)
           spbuf(1,n)=min(spbuf(1,n),dwa(n)*spbuf(1,n))
           spbuf(8,n)=spbuf(1,n)
          END DO
        END IF
C
        DO ns=1,nsphact
         n=waspact(ns)
C
         IF(mod(kreduce(n),10)/=0)THEN
C
           nvois1 =kxsp(4,n)
           nvois  =kxsp(5,n)
           inod=kxsp(3,n)
           xi=x(1,inod)
           yi=x(2,inod)
           zi=x(3,inod)
           di=spbuf(1,n)
C
           jnod = ixsp(nvois,n)
           m =nod2sp(jnod)
           xj=x(1,jnod)
           yj=x(2,jnod)
           zj=x(3,jnod)
           dj=spbuf(1,m)
           dd  =(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj)
           dms =di+dj
           dms2=dms*dms
           dk=dd/dms2
           myspatrue=max(zero,min(myspatrue,dk-one))
         END IF
C
       END DO
C-------------------------------------------
 100   CONTINUE
C-----------------------------------------------
      RETURN
      END