sponof1.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/.
!||====================================================================
!||    sponof1        ../engine/source/elements/sph/sponof1.F
!||--- called by ------------------------------------------------------
!||    sphprep        ../engine/source/elements/sph/sphprep.F
!||--- calls      -----------------------------------------------------
!||    ancmsg         ../engine/source/output/message/message.F
!||    arret          ../engine/source/system/arret.F
!||    initbuf        ../engine/share/resol/initbuf.F
!||    sphreqs        ../engine/source/elements/sph/sphreq.F
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod   ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    groupdef_mod   ../common_source/modules/groupdef_mod.F
!||    initbuf_mod    ../engine/share/resol/initbuf.F
!||    message_mod    ../engine/share/message_module/message_mod.F
!||    outputs_mod    ../common_source/modules/outputs_mod.F
!||====================================================================
      SUBROUTINE sponof1(X ,V       ,D       ,MS      ,SPBUF   ,
     2             ITAB    ,KXSP    ,IXSP    ,NOD2SP  ,NPC     ,
     3             PLD     ,IPARG   ,ELBUF_TAB,ISPHIO  ,VSPHIO  ,
     4             PM      ,IPART   ,IPARTSP ,IGRSURF ,
     5             LPRTSPH ,LONFSPH ,MWA     ,WA      ,VNORMAL ,
     6             XDP, IBUFSSG_IO  ,OFF_SPH_R2R, WFEXT)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE initbuf_mod
      USE elbufdef_mod            
      USE message_mod
      USE groupdef_mod
      USE outputs_mod
C----6------------------------------------------
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "vect01_c.inc"
#include      "com01_c.inc"
#include      "com04_c.inc"
#include      "com06_c.inc"
#include      "com08_c.inc"
#include      "param_c.inc"
#include      "sphcom.inc"
#include      "scr05_c.inc"
#include      "tabsiz_c.inc"
#include      "scr17_c.inc"
#include      "rad2r_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),ITAB(*),NPC(*),
     .        IPARG(NPARG,*),ISPHIO(NISPHIO,*),IPART(LIPART1,*),
     .        LPRTSPH(2,0:NPART),LONFSPH(*),
     .        MWA(*), IBUFSSG_IO(SIBUFSSG_IO),IPARTSP(*),OFF_SPH_R2R(*)
      my_real
     .   x(3,*) ,v(3,*) ,d(3,*) ,ms(*) ,spbuf(nspbuf,*) ,
     .   pld(*) ,vsphio(*) ,pm(npropm,*) ,wa(*),
     .   vnormal(3,*)
      double precision
     .   xdp(3,*)
      TYPE (ELBUF_STRUCT_), TARGET, DIMENSION (NGROUP) :: ELBUF_TAB
      TYPE (SURF_)   , DIMENSION(NSURF)   :: IGRSURF
      DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I, ITYPE, IVAD,IX, N2ENTER,
     .        ii,ipt,jj,npf,
     .        ifvits,ifdens,ifpres,ifener,
     .        inactiv,j,n,inod,iactive,
     .        isu,nseg,in1,in2,in3,in4,
     .        imat,iprop,iprt,nel,kad,ng,k,
     .        kp,np2sort,nband,
     .        knsx,kisx,knsy,kisy,knsz,kisz,
     .        knpx,kipx,knpy,kipy,knpz,kipz,kvnorm,kxproj,kf,
     .        iun,iad2,nsegb
      my_real
     .       pentv,vt,vn,vv,
     .       pentr,rhon,rho0,pentp,pn,pshft,pente,en,ei,
     .       t05,sseg,mp,wfextt,volo,
     .       x1,x2,x3,x4,xp,y1,y2,y3,y4,yp,z1,z2,z3,z4,zp,
     .       x12,y12,z12,x23,y23,z23,x31,y31,z31,nx,ny,nz,nn,
     .       xi,yi,zi,dd,diam,dmax,dps,vol,bid
      TYPE(g_bufel_)  ,POINTER :: GBUF     
C-----------------------------------------------
      my_real
     .         GET_U_GEO
      EXTERNAL GET_U_GEO
      DATA iun/1/
C-----------------------------------------------
      wfextt=zero
      bid=zero
      en = 0
      rhon = zero
C-------
C     reset active particles states.
      DO i=1,nsphio
       iprt   =isphio(2,i)
C      loop over active particles related to IPRT.
       DO iactive=lprtsph(2,iprt-1)+1,lprtsph(1,iprt)
        n=lonfsph(iactive)
        kxsp(2,n)=mod(kxsp(2,n),ngroup+1)
       ENDDO
      ENDDO
C-----------------------------------------------
      DO i=1,nsphio
      itype=isphio(1,i)
      IF(itype==1)THEN
       iprt   =isphio(2,i)
C------
C      general inlet.
       ivad=isphio(4,i)
       t05=tt - half*dt1
C
       ifvits=isphio(8,i)
       vn    =zero
       IF(ifvits/=0)THEN
        npf = (npc(ifvits+1)-npc(ifvits))/2
        ii  = npc(ifvits)
        IF (t05<=pld(ii)) THEN
          pentv=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
          vn   =pld(ii+1)+pentv*(t05-pld(ii))
        ELSEIF (t05>=pld(ii+2*(npf-1))) THEN
          jj=ii+2*(npf-1)
          pentv=(pld(jj+1)-pld(jj-1))/(pld(jj)-pld(jj-2))
          vn   =pld(jj+1)+max(-pld(jj+1),pentv*(t05-pld(jj)))
        ELSE
          DO ipt=1,npf-1
           IF (pld(ii)<=t05.AND.t05<=pld(ii+2)) THEN
            pentv=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
            vn   =pld(ii+1)+pentv*(t05-pld(ii))
            GOTO 20
           ENDIF
           ii=ii+2
          ENDDO
 20       CONTINUE
        ENDIF
       ENDIF
C------
       vv =abs(vn)
       vt =vv*vv
C------
       ifdens=isphio(5,i)
       IF(ifdens==0)THEN
        rhon=vsphio(ivad)
       ELSE
        npf = (npc(ifdens+1)-npc(ifdens))/2
        ii  = npc(ifdens)
        IF (tt<=pld(ii)) THEN
          pentr=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
          rhon =pld(ii+1)+pentr*(tt-pld(ii))
        ELSEIF (tt>=pld(ii+2*(npf-1))) THEN
          jj=ii+2*(npf-1)
          pentr=(pld(jj+1)-pld(jj-1))/(pld(jj)-pld(jj-2))
          rhon =pld(jj+1)+max(-pld(jj+1),pentr*(tt-pld(jj)))
        ELSE
          DO ipt=1,npf-1
           IF (pld(ii)<=tt.AND.tt<=pld(ii+2)) THEN
            pentr=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
            rhon =pld(ii+1)+pentr*(tt-pld(ii))
            GOTO 30
           ENDIF
           ii=ii+2
          ENDDO
 30       CONTINUE
        ENDIF
        rhon=rhon*vsphio(ivad)
       ENDIF
C------
       ifener=isphio(7,i)
       IF(ifener==0)THEN
        en=vsphio(ivad+2)
       ELSE
        npf = (npc(ifener+1)-npc(ifener))/2
        ii  = npc(ifener)
        IF (tt<=pld(ii)) THEN
          pente=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
          en   =pld(ii+1)+pente*(tt-pld(ii))
        ELSEIF (tt>=pld(ii+2*(npf-1))) THEN
          jj=ii+2*(npf-1)
          pente=(pld(jj+1)-pld(jj-1))/(pld(jj)-pld(jj-2))
          en   =pld(jj+1)+max(-pld(jj+1),pente*(tt-pld(jj)))
        ELSE
          DO ipt=1,npf-1
           IF (pld(ii)<=tt.AND.tt<=pld(ii+2)) THEN
            pente=(pld(ii+3)-pld(ii+1))/(pld(ii+2)-pld(ii))
            en   =pld(ii+1)+pente*(tt-pld(ii))
            GOTO 50
           ENDIF
           ii=ii+2
          ENDDO
 50       CONTINUE
        ENDIF
        en=en*vsphio(ivad+2)
       ENDIF
C------
       iprop=ipart(2,iprt)
       diam =get_u_geo(6,iprop)
       imat =ipart(1,iprt)
C      specific energy related to V0=MP/RHO0.
       rho0 =pm(1 ,imat)
       ei   =en*(rho0/rhon)
       
c default value        
       IF(diam==zero) THEN
         mp  = get_u_geo(1,iprop)
         vol = mp/rho0
         diam= (sqr2*vol)**third
       END IF             
C------
C      boite englobante / surface entiere.
       xbmin =ep20
       ybmin =ep20
       zbmin =ep20
       xbmax=-ep20
       ybmax=-ep20
       zbmax=-ep20       
C------
       isu =isphio(3,i)
c NSEGB on split surface (IBUFSSG)
       nsegb=igrsurf(isu)%NSEG
c NSEG on complete surface (IBUFSSG_IO)
       nseg=isphio(10,i)
c IAD on split surface (IBUFSSG)
c IAD2 on complete surface (IBUFSSG_IO)
       iad2 =isphio(11,i)
       dmax=vsphio(ivad+3)
 
c complete surface is used (IBUFSSG_IO) as in sponof2
c call of SPHREQS with IBUFSSG_IO
       DO j=0,nseg-1
          in1=ibufssg_io(iad2+nibsph*j  )
          in2=ibufssg_io(iad2+nibsph*j+1)
          in3=ibufssg_io(iad2+nibsph*j+2)
          in4=ibufssg_io(iad2+nibsph*j+3)
          xp=x(1,in1)
          yp=x(2,in1)
          zp=x(3,in1)
          xbmin=min(xbmin,xp)
          ybmin=min(ybmin,yp)
          zbmin=min(zbmin,zp)
          xbmax=max(xbmax,xp)
          ybmax=max(ybmax,yp)
          zbmax=max(zbmax,zp)
          xp=x(1,in2)
          yp=x(2,in2)
          zp=x(3,in2)
          xbmin=min(xbmin,xp)
          ybmin=min(ybmin,yp)
          zbmin=min(zbmin,zp)
          xbmax=max(xbmax,xp)
          ybmax=max(ybmax,yp)
          zbmax=max(zbmax,zp)
          xp=x(1,in3)
          yp=x(2,in3)
          zp=x(3,in3)
          xbmin=min(xbmin,xp)
          ybmin=min(ybmin,yp)
          zbmin=min(zbmin,zp)
          xbmax=max(xbmax,xp)
          ybmax=max(ybmax,yp)
          zbmax=max(zbmax,zp)
          IF(in4/=in3)THEN
           xp=x(1,in4)
           yp=x(2,in4)
           zp=x(3,in4)
           xbmin=min(xbmin,xp)
           ybmin=min(ybmin,yp)
           zbmin=min(zbmin,zp)
           xbmax=max(xbmax,xp)
           ybmax=max(ybmax,yp)
           zbmax=max(zbmax,zp)
          ENDIF
       ENDDO
       xbmin=xbmin-dmax
       ybmin=ybmin-dmax
       zbmin=zbmin-dmax
       xbmax=xbmax+dmax
       ybmax=ybmax+dmax
       zbmax=zbmax+dmax
C------
       np2sort=0
       DO iactive=lprtsph(2,iprt-1)+1,lprtsph(1,iprt)
        n=lonfsph(iactive)
        inod=kxsp(3,n)
        xi=x(1,inod)
        yi=x(2,inod)
        zi=x(3,inod)
        IF(     xi>xbmin.AND.yi>ybmin.AND.zi>zbmin
     .     .AND.xi<xbmax.AND.yi<ybmax.AND.zi<zbmax)THEN
         np2sort=np2sort+1
         mwa(np2sort)=iactive
        ENDIF
       ENDDO
C------
C      bucket.
       dbucs=dmax
       DO j=1,nseg
          in1=ibufssg_io(iad2+nibsph*(j-1)  )
          in2=ibufssg_io(iad2+nibsph*(j-1)+1)
          in3=ibufssg_io(iad2+nibsph*(j-1)+2)
          in4=ibufssg_io(iad2+nibsph*(j-1)+3)
          x1=x(1,in1)
          y1=x(2,in1)
          z1=x(3,in1)
          x2=x(1,in2)
          y2=x(2,in2)
          z2=x(3,in2)
          x3=x(1,in3)
          y3=x(2,in3)
          z3=x(3,in3)
          dbucs=max(dbucs,abs(x1-x2))
          dbucs=max(dbucs,abs(y1-y2))
          dbucs=max(dbucs,abs(z1-z2))
          dbucs=max(dbucs,abs(x2-x3))
          dbucs=max(dbucs,abs(y2-y3))
          dbucs=max(dbucs,abs(z2-z3))
          dbucs=max(dbucs,abs(x3-x1))
          dbucs=max(dbucs,abs(y3-y1))
          dbucs=max(dbucs,abs(z3-z1))
          IF(in4/=in3)THEN
           x4=x(1,in4)
           y4=x(2,in4)
           z4=x(3,in4)
          ELSE
           x4=x3
           y4=y3
           z4=z3
          ENDIF
          dbucs=max(dbucs,abs(x1-x4))
          dbucs=max(dbucs,abs(y1-y4))
          dbucs=max(dbucs,abs(z1-z4))
          dbucs=max(dbucs,abs(x2-x4))
          dbucs=max(dbucs,abs(y2-y4))
          dbucs=max(dbucs,abs(z2-z4))
          dbucs=max(dbucs,abs(x3-x4))
          dbucs=max(dbucs,abs(y3-y4))
          dbucs=max(dbucs,abs(z3-z4))
       ENDDO
C      if more than 1 box, each box size >= DBUCS
       nbox =max(iun,int((xbmax-xbmin)/dbucs))
       nboy =max(iun,int((ybmax-ybmin)/dbucs))
       nboz =max(iun,int((zbmax-zbmin)/dbucs))
       nband=max(nbox,nboy,nboz)+1
       knsx=numsph+1
       kisx=knsx+nband+1
       knsy=kisx+4*nseg
       kisy=knsy+nband+1
       knsz=kisy+4*nseg
       kisz=knsz+nband+1
       knpx=kisz+4*nseg
       kipx=knpx+nband+1
       knpy=kipx+3*np2sort
       kipy=knpy+nband+1
       knpz=kipy+3*np2sort
       kipz=knpz+nband+1
       kvnorm=kipz+3*np2sort
C      KF    =KVNORM+3*NUMSPH
       CALL sphreqs(nseg ,ibufssg_io(iad2) ,x ,np2sort ,mwa ,
     2    lonfsph  ,kxsp     ,wa(knsx) ,wa(kisx) ,wa(knsy) ,
     3    wa(kisy) ,wa(knsz) ,wa(kisz) ,wa(knpx) ,wa(kipx) ,
     4    wa(knpy) ,wa(kipy) ,wa(knpz) ,wa(kipz) ,wa       ,
     5    wa(kvnorm),bid     ,bid      ,v        ,spbuf    ,
     6    isphio(1,i))
C------
C      check active particles within distance / re-impose kinematic condition.
       DO kp=1,np2sort
        dps=wa(kp)
        IF(dps>=0..AND.dps<dmax)THEN
         iactive=mwa(kp)
         n=lonfsph(iactive)
         inod=kxsp(3,n)
         wfextt=wfextt + half*ms(inod)*(vt-
     . (v(1,inod)*v(1,inod)+v(2,inod)*v(2,inod)+v(3,inod)*v(3,inod)))
         vnormal(1,n)=-wa(kvnorm+3*(n-1))
         vnormal(2,n)=-wa(kvnorm+3*(n-1)+1)
         vnormal(3,n)=-wa(kvnorm+3*(n-1)+2)
         v(1,inod)=vn*vnormal(1,n)
         v(2,inod)=vn*vnormal(2,n)
         v(3,inod)=vn*vnormal(3,n)
         ng=mod(kxsp(2,n),ngroup+1)
         kxsp(2,n) =ng+i*(ngroup+1)
        ENDIF
       ENDDO
C------
C      particles to enter.
       mp   =get_u_geo(1,iprop)
       n2enter=0
       ix=ivad+4
c loop on NSEGB (segments of split surface IBUFSSG) to avoid
c multiple count of cells to enter and to activate in SPMD
       DO j=0,(nsegb-1)
        sseg =vsphio(ix+1)
        vsphio(ix)=vsphio(ix)+rhon*sseg*vv*dt1
        IF(vsphio(ix)>=mp)THEN
         n2enter=n2enter+1
        ENDIF
        ix=ix+2
       ENDDO
C------
       IF(n2enter>0)THEN
        isphbuc=1
C       find N2ENTER inactive particles related to same part.
        inactiv=lprtsph(2,iprt)-lprtsph(1,iprt)
        IF(inactiv<n2enter)THEN
          CALL ancmsg(msgid=175,anmode=aninfo,
     .            i1=isphio(nisphio,i))
          CALL arret(2)
          RETURN
        ENDIF
C------
C       enter N2ENTER particles.
        ix=ivad+4
        iactive=lprtsph(1,iprt)
        DO j=0,(nsegb-1)
        IF(vsphio(ix)>=mp)THEN
C--------
         iactive=iactive+1
          n     =lonfsph(iactive)
C         chainage.
          lprtsph(1,iprt)=iactive
C--------
          in1=igrsurf(isu)%NODES(j+1,1)
          in2=igrsurf(isu)%NODES(j+1,2)
          in3=igrsurf(isu)%NODES(j+1,3)
          in4=igrsurf(isu)%NODES(j+1,4)
c          IN1=IBUFSSG_IO(IAD2+NIBSPH*(J-1)  )
c          IN2=IBUFSSG_IO(IAD2+NIBSPH*(J-1)+1)
c          IN3=IBUFSSG_IO(IAD2+NIBSPH*(J-1)+2)
c          IN4=IBUFSSG_IO(IAD2+NIBSPH*(J-1)+3)
          x1=x(1,in1)
          y1=x(2,in1)
          z1=x(3,in1)
          x2=x(1,in2)
          y2=x(2,in2)
          z2=x(3,in2)
          x3=x(1,in3)
          y3=x(2,in3)
          z3=x(3,in3)
         IF(in4/=in3)THEN
          x4=x(1,in4)
          y4=x(2,in4)
          z4=x(3,in4)
          xp=fourth*(x1+x2+x3+x4)
          yp=fourth*(y1+y2+y3+y4)
          zp=fourth*(z1+z2+z3+z4)    
         ELSE
          x4=x3
          y4=y3
          z4=z3
          xp=third*(x1+x2+x3)
          yp=third*(y1+y2+y3)
          zp=third*(z1+z2+z3)
         ENDIF
C--------
         x12=half*(x2+x3-x1-x4)
         y12=half*(y2+y3-y1-y4)
         z12=half*(z2+z3-z1-z4)
         x23=half*(x3+x4-x1-x2)
         y23=half*(y3+y4-y1-y2)
         z23=half*(z3+z4-z1-z2)   
         nx = y12*z23-z12*y23
         ny =-x12*z23+z12*x23
         nz = x12*y23-y12*x23
         nn =one/max(em20,sqrt(nx*nx+ny*ny+nz*nz))
         nx=nn*nx
         ny=nn*ny
         nz=nn*nz
C--------
         inod=kxsp(3,n)
         IF (itab(inod) >= 1000000000) THEN
           d(1,inod)=xp-xi_res
           d(2,inod)=yp-yi_res
           d(3,inod)=zp-zi_res
         ELSE
           d(1,inod)=xp-x(1,inod)+d(1,inod)
           d(2,inod)=yp-x(2,inod)+d(2,inod)
           d(3,inod)=zp-x(3,inod)+d(3,inod)
         ENDIF
         x(1,inod)=xp
         x(2,inod)=yp
         x(3,inod)=zp
         IF (iresp==1) THEN
           xdp(1,inod)=xp
           xdp(2,inod)=yp
           xdp(3,inod)=zp
         END IF
         wfextt=wfextt+half*mp*vt
         vnormal(1,n)=nx
         vnormal(2,n)=ny
         vnormal(3,n)=nz
         v(1,inod)=vn*nx
         v(2,inod)=vn*ny
         v(3,inod)=vn*nz
         IF (r2r_siu/=0)  off_sph_r2r(inod) = 2
C--------
          ng=-kxsp(2,n)
          iparg(8,ng) =0
          iparg(10,ng)=iparg(10,ng)+1
          CALL initbuf(iparg    ,ng      ,                   
     2       mtn     ,nel     ,nft     ,kad     ,ity     ,    
     3       npt     ,jale    ,ismstr  ,jeul    ,jtur    ,    
     4       jthe    ,jlag    ,jmult   ,jhbe    ,jivf    ,    
     5       nvaux   ,jpor    ,jcvt    ,jclose  ,jplasol ,    
     6       irep    ,iint    ,igtyp   ,israt   ,isrot   ,    
     7       icsen   ,isorth  ,isorthg ,ifailure,jsms    )                      
          gbuf => elbuf_tab(ng)%GBUF
          k=n-nft
          volo=gbuf%VOL(k)
          wfextt=wfextt+ei*volo
          gbuf%EINT(k) =ei
          gbuf%OFF(k)  =one
          gbuf%RHO(k)  =rhon
          spbuf(1,n)=diam
          spbuf(2,n) =rhon
          spbuf(12,n)=ms(inod)
          kxsp(2,n) =ng+i*(ngroup+1)
C-------
C        reset MASS FLOW for this segment.
         vsphio(ix)=vsphio(ix)-mp
        ENDIF
        ix=ix+2
        ENDDO
       ENDIF
      ENDIF
      ENDDO
C-------------------------------------------
!$OMP ATOMIC
       wfext=wfext+wfextt
C-------------------------------------------
      RETURN
      END