spinih.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/.
!||====================================================================
!||    spinih         ../starter/source/elements/sph/spinih.F
!||--- called by ------------------------------------------------------
!||    lectur         ../starter/source/starter/lectur.F
!||--- calls      -----------------------------------------------------
!||    get_u_geo      ../starter/source/user_interface/uaccess.F
!||    reset_u_geo    ../starter/source/user_interface/uaccess.F
!||    sptri          ../starter/source/elements/sph/sptri.F
!||--- uses       -----------------------------------------------------
!||====================================================================
      SUBROUTINE spinih(KXSP,IPART,IPARTSP,SPBUF,PM,
     .                  IXSP    ,NOD2SP ,X       ,LPRTSPH,LONFSPH,
     .                  SNOD2SP ,SLONFSPH,NUMNOD,NPART,ITAB)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE my_alloc_mod
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      "param_c.inc"
#include      "sphcom.inc"
#include      "scr17_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER KXSP(NISP,*), IPART(LIPART1,*), IPARTSP(*),ITAB(*)
      INTEGER ,INTENT(IN) :: SNOD2SP,SLONFSPH,NUMNOD,NPART
      INTEGER ,INTENT(INOUT) :: IXSP(KVOISPH,NUMSPH),NOD2SP(SNOD2SP),LPRTSPH(2,0:NPART),
     .                          lonfsph(slonfsph)
      my_real
     .        spbuf(nspbuf,*),pm(npropm,*)
      my_real ,INTENT(IN) :: x(3,numnod)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER N,IPRT,I0,IMAT,PRE_SEARCH,FLG_H,J,NS,NBNOD(NPART),SZ_INTP_DIST
      my_real, DIMENSION(:), ALLOCATABLE :: MAX_INTP_DIST_PART   ! Max of interparticle distance per part
      my_real
     .        dist, mp, rho, vol, sq2
C-----------------------------------------------
      my_real
     .        get_u_geo,h_scal,vol_tot,reset_u_geo
      my_real
     .        xmin(npart),xmax(npart),ymin(npart),ymax(npart),zmin(npart),
     .        zmax(npart),volp_prt(npart)  
C-----------------------------------------------
      pre_search = 0
      sq2 = sqrt(two)
C
      volp_prt(1:npart)=zero
      nbnod(1:npart)=0
      xmin(1:npart)=ep30
      xmax(1:npart)=-ep30
      ymin(1:npart)=ep30
      ymax(1:npart)=-ep30
      zmin(1:npart)=ep30
      zmax(1:npart)=-ep30 
C      
      DO n = 1, numsph
        iprt=ipartsp(n)
        imat=ipart(1,iprt)
        i0  =ipart(2,iprt)
        dist=get_u_geo(6,i0)           
        IF (nint(spbuf(13,n))==1) THEN
C--       type = 1 - particle with mass input   
          mp  = spbuf(12,n)
          rho=pm(1,imat)
          vol=mp/rho
          dist=(sq2*vol)**third    
        ELSEIF (nint(spbuf(13,n))==2) THEN
C--       type = 2 - particle with volume input          
          vol=spbuf(12,n)
          dist=(sq2*vol)**third             
        ELSE           
C--       type = 0 - particle without mass/vol input            
          IF(dist==zero) THEN
            mp  =get_u_geo(1,i0)
            rho=pm(1,imat)
            vol=mp/rho
            dist=(sq2*vol)**third     
          END IF
        ENDIF  
        spbuf(1,n)=dist
C        
        flg_h=nint(get_u_geo(9,i0))           
        IF (flg_h==3) THEN           
C--       FLG_H=3 - computation of H from interparticle distance - presearch of neigbours needed        
          pre_search = 1  
C--       Approximation of H         
          j=kxsp(3,n) 
          nbnod(iprt)=nbnod(iprt)+1
          xmin(iprt)= min(xmin(iprt),x(1,j))
          ymin(iprt)= min(ymin(iprt),x(2,j))
          zmin(iprt)= min(zmin(iprt),x(3,j))
          xmax(iprt)= max(xmax(iprt),x(1,j))
          ymax(iprt)= max(ymax(iprt),x(2,j))
          zmax(iprt)= max(zmax(iprt),x(3,j))          
        ENDIF  
C        
      END DO
C      
      IF (pre_search == 1) THEN
C       Rough approxiamtion of H per part - only used for presearch (important for SOL2SPH because MP not defined)
        DO iprt = 1, npart
          IF (nbnod(iprt) > 0) THEN  
            vol_tot = abs(xmin(iprt)-xmax(iprt))*abs(ymin(iprt)-ymax(iprt))*abs(zmin(iprt)-zmax(iprt))
            volp_prt(iprt) = vol_tot/nbnod(iprt)
          ENDIF
        END DO 
        DO n=1,numsph
          iprt=ipartsp(n)
          i0  =ipart(2,iprt)
          flg_h=nint(get_u_geo(9,i0))
          IF (flg_h==3) THEN
            dist = volp_prt(iprt)**third
            spbuf(1,n) = min(dist,spbuf(1,n)) 
          ENDIF  
        END DO
C             
        CALL my_alloc(max_intp_dist_part,npart)
        max_intp_dist_part(1:npart) = zero
        sz_intp_dist = npart
c
C       Pre-search of neigbours to define interparticular distance  
        CALL sptri(kxsp    ,ixsp    ,nod2sp  ,x      ,spbuf  ,
     .             lprtsph ,lonfsph ,ipartsp ,sz_intp_dist,max_intp_dist_part,
     .             pre_search)                   
C     
        DO n = 1, numsph
          iprt=ipartsp(n)
          i0  =ipart(2,iprt)    
          flg_h=nint(get_u_geo(9,i0))
C          KXSP(9,N) = FLG_H
          IF (flg_h==3) THEN
            h_scal=get_u_geo(12,i0)
            spbuf(1,n)=h_scal*max_intp_dist_part(iprt)    ! DIamter of the particle 
            spbuf(14,n)=h_scal*max_intp_dist_part(iprt)   ! Diameter of the particle at time zero
          ENDIF  
        END DO     
C
        DEALLOCATE(max_intp_dist_part)
C
      ENDIF
C
      RETURN
      END