soltosph.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/.
!||====================================================================
!||    soltosphx8       ../starter/source/elements/sph/soltosph.F
!||--- called by ------------------------------------------------------
!||    hm_read_sphcel   ../starter/source/elements/reader/hm_read_sphcel.F
!||--- uses       -----------------------------------------------------
!||    message_mod      ../starter/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE soltosphx8(
     .   NSPHDIR ,NCELL   ,INOD    ,IDS     ,IDMAX  ,
     .   X       ,IXS     ,KXSP    ,IPARTSP ,NOD2SP ,
     .   IRST    )
C-----------------------------------------------
C    M o d u l e s
C-----------------------------------------------
      USE message_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-----------------------------------------------
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 NSPHDIR, NCELL, INOD, IDS, IDMAX, IXS(NIXS),
     .        KXSP(NISP,*), IPARTSP(*), NOD2SP(*), IRST(3,*)
       my_real
     .        x(3,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,  J, IR, IS, IT, 
     .        N1, N2, N3, N4, N5, N6, N7, N8
C                                                                     12
      my_real
     .   x1,x2,x3,x4,x5,x6,x7,x8,
     .   y1,y2,y3,y4,y5,y6,y7,y8,
     .   z1,z2,z3,z4,z5,z6,z7,z8,
     .   phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8,
     .   ksi, eta, zeta, wi,
     .   xi, yi, zi
C-----------------------------------------------
      my_real
     .  a_gauss(9,9)
      DATA a_gauss /
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     2 -.5              ,0.5              ,0.               ,
     2 0.               ,0.               ,0.               ,
     2 0.               ,0.               ,0.               ,
     3 -.666666666666666,0.               ,0.666666666666666,
     3 0.               ,0.               ,0.               ,
     3 0.               ,0.               ,0.               ,
     4 -.75             ,-.25             ,0.25             ,
     4 0.75             ,0.               ,0.               ,
     4 0.               ,0.               ,0.               ,
     5 -.8              ,-.4              ,0.               ,
     5 0.4              ,0.8              ,0.               ,
     5 0.               ,0.               ,0.               ,
     6 -.833333333333333,-.5              ,-.166666666666666,
     6 0.166666666666666,0.5              ,0.833333333333333,
     6 0.               ,0.               ,0.               ,
     7 -.857142857142857,-.571428571428571,-.285714285714285,
     7 0.               ,0.285714285714285,0.571428571428571,
     7 0.857142857142857,0.               ,0.               ,
     8 -.875            ,-.625            ,-.375            ,
     8 -.125            ,0.125            ,0.375,
     8 0.625            ,0.875            ,0.               ,
     9 -.888888888888888,-.666666666666666,-.444444444444444,
     9 -.222222222222222,0.               ,0.222222222222222,
     9 0.444444444444444,0.666666666666666,0.888888888888888/
C-----------------------------------------------
      n1=ixs(2)
      x1=x(1,n1)
      y1=x(2,n1)
      z1=x(3,n1)
      n2=ixs(3)
      x2=x(1,n2)
      y2=x(2,n2)
      z2=x(3,n2)
      n3=ixs(4)
      x3=x(1,n3)
      y3=x(2,n3)
      z3=x(3,n3)
      n4=ixs(5)
      x4=x(1,n4)
      y4=x(2,n4)
      z4=x(3,n4)
      n5=ixs(6)
      x5=x(1,n5)
      y5=x(2,n5)
      z5=x(3,n5)
      n6=ixs(7)
      x6=x(1,n6)
      y6=x(2,n6)
      z6=x(3,n6)
      n7=ixs(8)
      x7=x(1,n7)
      y7=x(2,n7)
      z7=x(3,n7)
      n8=ixs(9)
      x8=x(1,n8)
      y8=x(2,n8)
      z8=x(3,n8)
C-----------------------------------------------
      DO ir=1,nsphdir
      DO is=1,nsphdir
      DO it=1,nsphdir
       ksi  = a_gauss(ir,nsphdir)
       eta  = a_gauss(is,nsphdir)
       zeta = a_gauss(it,nsphdir)
C
       phi1=(one-ksi)*(one-eta)*(one-zeta)
       phi2=(one-ksi)*(one-eta)*(one+zeta)
       phi3=(one+ksi)*(one-eta)*(one+zeta)
       phi4=(one+ksi)*(one-eta)*(one-zeta)
       phi5=(one-ksi)*(one+eta)*(one-zeta)
       phi6=(one-ksi)*(one+eta)*(one+zeta)
       phi7=(one+ksi)*(one+eta)*(one+zeta)
       phi8=(one+ksi)*(one+eta)*(one-zeta)
       xi=one_over_8*(phi1*x1+phi2*x2+phi3*x3+phi4*x4+
     .            phi5*x5+phi6*x6+phi7*x7+phi8*x8)
       yi=one_over_8*(phi1*y1+phi2*y2+phi3*y3+phi4*y4+
     .            phi5*y5+phi6*y6+phi7*y7+phi8*y8)
       zi=one_over_8*(phi1*z1+phi2*z2+phi3*z3+phi4*z4+
     .            phi5*z5+phi6*z6+phi7*z7+phi8*z8)
C
       ncell=ncell+1
       ipartsp(ncell)=ids
       inod =inod+1
       kxsp(3,ncell) =inod
       x(1,inod)=xi
       x(2,inod)=yi
       x(3,inod)=zi
       nod2sp(inod)  =ncell
       kxsp(2,ncell)=-1
       idmax=idmax+1
       kxsp(nisp,ncell)=idmax
       irst(1,ncell-first_sphsol+1)=ir
       irst(2,ncell-first_sphsol+1)=is
       irst(3,ncell-first_sphsol+1)=it
C
      ENDDO
      ENDDO
      ENDDO
C-----------------------------------------------
      RETURN
      END SUBROUTINE soltosphx8
!||====================================================================
!||    soltosphx4       ../starter/source/elements/sph/soltosph.F
!||--- called by ------------------------------------------------------
!||    hm_read_sphcel   ../starter/source/elements/reader/hm_read_sphcel.F
!||--- calls      -----------------------------------------------------
!||    checkvolume_4n   ../starter/source/elements/solid/solide/checksvolume.F
!||--- uses       -----------------------------------------------------
!||    message_mod      ../starter/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE soltosphx4(
     .   NSPHDIR ,NCELL   ,INOD    ,IDS     ,IDMAX  ,
     .   X       ,IXS     ,KXSP    ,IPARTSP ,NOD2SP ,
     .   IRST    )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_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-----------------------------------------------
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 NSPHDIR, NCELL, INOD, IDS, IDMAX, IXS(NIXS),
     .        KXSP(NISP,*), IPARTSP(*), NOD2SP(*), IRST(3,*)
       my_real
     .        X(3,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER IR, IS, IT,N1,N2,N3,N4
C
      my_real
     .   x1,x2,x3,x4,y1,y2,y3,y4,
     .   z1,z2,z3,z4,phi1,phi2,phi3,phi4,ksi,
     .   eta,zeta,xi,yi,zi,a_gauss_tetra(9,9)
C-----------------------------------------------
C   E x t e r n a l  F u n c t i o n s
C-----------------------------------------------
      my_real
     .   checkvolume_4n
C-----------------------------------------------
C     A_GAUSS Generated with (2*IR-1)/(2(NSPHDIR+1))
C-----------------------------------------------
      DATA a_gauss_tetra /
     1 0.250000000000000,0.000000000000000,0.000000000000000,
     1 0.000000000000000,0.000000000000000,0.000000000000000,
     1 0.000000000000000,0.000000000000000,0.000000000000000,
     2 0.166666666666667,0.500000000000000,0.000000000000000,
     2 0.000000000000000,0.000000000000000,0.000000000000000,
     2 0.000000000000000,0.000000000000000,0.000000000000000,
     3 0.125000000000000,0.375000000000000,0.625000000000000,
     3 0.000000000000000,0.000000000000000,0.000000000000000,
     3 0.000000000000000,0.000000000000000,0.000000000000000,
     4 0.100000000000000,0.300000000000000,0.500000000000000,
     4 0.700000000000000,0.000000000000000,0.000000000000000,
     4 0.000000000000000,0.000000000000000,0.000000000000000,
     5 0.083333333333333,0.250000000000000,0.416666666666667,
     5 0.583333333333333,0.750000000000000,0.000000000000000,
     5 0.000000000000000,0.000000000000000,0.000000000000000,
     6 0.071428571428571,0.214285714285714,0.357142857142857,
     6 0.500000000000000,0.642857142857143,0.785714285714286,
     6 0.000000000000000,0.000000000000000,0.000000000000000,
     7 0.062500000000000,0.187500000000000,0.312500000000000,
     7 0.437500000000000,0.562500000000000,0.687500000000000,
     7 0.812500000000000,0.000000000000000,0.000000000000000,
     8 0.055555555555556,0.166666666666667,0.277777777777778,
     8 0.388888888888889,0.500000000000000,0.611111111111111,
     8 0.722222222222222,0.833333333333333,0.000000000000000,
     9 0.050000000000000,0.150000000000000,0.250000000000000,
     9 0.350000000000000,0.450000000000000,0.550000000000000,
     9 0.650000000000000,0.750000000000000,0.850000000000000/
C-----------------------------------------------
C---- KSI - R - N4->N1
C---- ETA - S - N4->N2
C---- ZETA- T - N4->N3  
C-----------------------------------------------
      n1=ixs(2)
      x1=x(1,n1)
      y1=x(2,n1)
      z1=x(3,n1)
      n2=ixs(4)
      x2=x(1,n2)
      y2=x(2,n2)
      z2=x(3,n2)
      n3=ixs(7)
      x3=x(1,n3)
      y3=x(2,n3)
      z3=x(3,n3)
      n4=ixs(6)
      x4=x(1,n4)
      y4=x(2,n4)
      z4=x(3,n4)
 
C------------------------------------------------
C Renumber tetrahedron in case of negative volume
      IF (checkvolume_4n(x ,ixs(1)) < zero) THEN
         n2=ixs(6)
         n4=ixs(4)
      ENDIF
C-----------------------------------------------
      DO ir=1,nsphdir
      DO is=1,nsphdir-ir+1
      DO it=1,nsphdir-is-ir+2
C         
       ksi  = a_gauss_tetra(ir,nsphdir)
       eta  = a_gauss_tetra(is,nsphdir)
       zeta = a_gauss_tetra(it,nsphdir)
C
       phi1=ksi
       phi2=eta
       phi3=zeta
       phi4=one-ksi-eta-zeta
C
       xi=phi1*x1+phi2*x2+phi3*x3+phi4*x4
       yi=phi1*y1+phi2*y2+phi3*y3+phi4*y4
       zi=phi1*z1+phi2*z2+phi3*z3+phi4*z4
C
       ncell=ncell+1
       ipartsp(ncell)=ids
       inod =inod+1
       kxsp(3,ncell) =inod
       x(1,inod)=xi
       x(2,inod)=yi
       x(3,inod)=zi
       nod2sp(inod)  =ncell
       kxsp(2,ncell)=-1
       idmax=idmax+1
       kxsp(nisp,ncell)=idmax
       irst(1,ncell-first_sphsol+1)=ir
       irst(2,ncell-first_sphsol+1)=is
       irst(3,ncell-first_sphsol+1)=it
C
      ENDDO
      ENDDO
      ENDDO
C
C-----------------------------------------------
      RETURN
      END SUBROUTINE soltosphx4
!||====================================================================
!||    soltosphv8    ../starter/source/elements/sph/soltosph.f
!||--- called by ------------------------------------------------------
!||    s8zinit3      ../starter/source/elements/solid/solide8z/s8zinit3.f
!||    sinit3        ../starter/source/elements/solid/solide/sinit3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg        ../starter/source/output/message/message.F
!||--- uses       -----------------------------------------------------
!||    message_mod   ../starter/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE soltosphv8(
     .   NSPHDIR ,RHO     ,NCELL   ,X      ,SPBUF   ,
     .   IXS     ,KXSP    ,IPARTSP ,IRST   )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_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-----------------------------------------------
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 NSPHDIR, NCELL, IXS(NIXS), KXSP(NISP,*),
     .        IPARTSP(*), IRST(3,*)
      my_real
     .        RHO, X(3,*), SPBUF(NSPBUF,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,  J, IR, IS, IT, IP,
     .        n1, n2, n3, n4, n5, n6, n7, n8, np
C                                                                    
      my_real
     .   x1,x2,x3,x4,x5,x6,x7,x8,
     .   y1,y2,y3,y4,y5,y6,y7,y8,
     .   z1,z2,z3,z4,z5,z6,z7,z8,
     .   x17 , x28 , x35 , x46 ,
     .   y17 , y28 , y35 , y46 ,
     .   z17 , z28 , z35 , z46 ,
     .   vol, hx(4), hy(4), hz(4), det, 
     .   jac1 ,jac2 ,jac3 ,
     .   jac4 ,jac5 ,jac6 ,
     .   jac7 ,jac8 ,jac9 ,
     .   cj1 ,cj2 ,cj3 ,
     .   cj4 ,cj5 ,cj6 ,
     .   cj7 ,cj8 ,cj9 ,
     .   jac_59_68, jac_67_49, jac_48_57,
     .   jac_38_29, jac_19_37, jac_27_18,
     .   jac_26_35, jac_34_16, jac_15_24,
     .   x_17_46 , x_28_35 ,
     .   y_17_46 , y_28_35 ,
     .   z_17_46 , z_28_35 ,
     .   ksi, eta, zeta, wi,
     .   xi, yi, zi
C-----------------------------------------------
      my_real
     .  w_gauss(9,9),a_gauss(9,9)
      DATA w_gauss /
     1 2.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     2 1.               ,1.               ,0.               ,
     2 0.               ,0.               ,0.               ,
     2 0.               ,0.               ,0.               ,
     3 0.555555555555556,0.888888888888889,0.555555555555556,
     3 0.               ,0.               ,0.               ,
     3 0.               ,0.               ,0.               ,
     4 0.347854845137454,0.652145154862546,0.652145154862546,
     4 0.347854845137454,0.               ,0.               ,
     4 0.               ,0.               ,0.               ,
     5 0.236926885056189,0.478628670499366,0.568888888888889,
     5 0.478628670499366,0.236926885056189,0.               ,
     5 0.               ,0.               ,0.               ,
     6 0.171324492379170,0.360761573048139,0.467913934572691,
     6 0.467913934572691,0.360761573048139,0.171324492379170,
     6 0.               ,0.               ,0.               ,
     7 0.129484966168870,0.279705391489277,0.381830050505119,
     7 0.417959183673469,0.381830050505119,0.279705391489277,
     7 0.129484966168870,0.               ,0.               ,
     8 0.101228536290376,0.222381034453374,0.313706645877887,
     8 0.362683783378362,0.362683783378362,0.313706645877887,
     8 0.222381034453374,0.101228536290376,0.               ,
     9 0.081274388361574,0.180648160694857,0.260610696402935,
     9 0.312347077040003,0.330239355001260,0.312347077040003,
     9 0.260610696402935,0.180648160694857,0.081274388361574/
      DATA a_gauss /
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     1 0.               ,0.               ,0.               ,
     2 -.5              ,0.5              ,0.               ,
     2 0.               ,0.               ,0.               ,
     2 0.               ,0.               ,0.               ,
     3 -.666666666666666,0.               ,0.666666666666666,
     3 0.               ,0.               ,0.               ,
     3 0.               ,0.               ,0.               ,
     4 -.75             ,-.25             ,0.25             ,
     4 0.75             ,0.               ,0.               ,
     4 0.               ,0.               ,0.               ,
     5 -.8              ,-.4              ,0.               ,
     5 0.4              ,0.8              ,0.               ,
     5 0.               ,0.               ,0.               ,
     6 -.833333333333333,-.5              ,-.166666666666666,
     6 0.166666666666666,0.5              ,0.833333333333333,
     6 0.               ,0.               ,0.               ,
     7 -.857142857142857,-.571428571428571,-.285714285714285,
     7 0.               ,0.285714285714285,0.571428571428571,
     7 0.857142857142857,0.               ,0.               ,
     8 -.875            ,-.625            ,-.375            ,
     8 -.125            ,0.125            ,0.375,
     8 0.625            ,0.875            ,0.               ,
     9 -.888888888888888,-.666666666666666,-.444444444444444,
     9 -.222222222222222,0.               ,0.222222222222222,
     9 0.444444444444444,0.666666666666666,0.888888888888888/
C-----------------------------------------------
      np = nsphdir*nsphdir*nsphdir
C
      n1=ixs(2)
      x1=x(1,n1)
      y1=x(2,n1)
      z1=x(3,n1)
      n2=ixs(3)
      x2=x(1,n2)
      y2=x(2,n2)
      z2=x(3,n2)
      n3=ixs(4)
      x3=x(1,n3)
      y3=x(2,n3)
      z3=x(3,n3)
      n4=ixs(5)
      x4=x(1,n4)
      y4=x(2,n4)
      z4=x(3,n4)
      n5=ixs(6)
      x5=x(1,n5)
      y5=x(2,n5)
      z5=x(3,n5)
      n6=ixs(7)
      x6=x(1,n6)
      y6=x(2,n6)
      z6=x(3,n6)
      n7=ixs(8)
      x7=x(1,n7)
      y7=x(2,n7)
      z7=x(3,n7)
      n8=ixs(9)
      x8=x(1,n8)
      y8=x(2,n8)
      z8=x(3,n8)
C
      x17=x7-x1
      x28=x8-x2
      x35=x5-x3
      x46=x6-x4
      y17=y7-y1
      y28=y8-y2
      y35=y5-y3
      y46=y6-y4
      z17=z7-z1
      z28=z8-z2
      z35=z5-z3
      z46=z6-z4
C
C Jacobian matrix
      cj4=x17+x28-x35-x46
      cj5=y17+y28-y35-y46
      cj6=z17+z28-z35-z46
      x_17_46=x17+x46
      x_28_35=x28+x35
      y_17_46=y17+y46
      y_28_35=y28+y35
      z_17_46=z17+z46
      z_28_35=z28+z35
C
      cj7=x_17_46+x_28_35
      cj8=y_17_46+y_28_35
      cj9=z_17_46+z_28_35
      cj1=x_17_46-x_28_35
      cj2=y_17_46-y_28_35
      cj3=z_17_46-z_28_35
C Hourglass
C mode 1
C 1 1 -1 -1 -1 -1 1 1
      hx(1)=(x1+x2-x3-x4-x5-x6+x7+x8)
      hy(1)=(y1+y2-y3-y4-y5-y6+y7+y8)
      hz(1)=(z1+z2-z3-z4-z5-z6+z7+z8)
C mode 2
C 1 -1 -1 1 -1 1 1 -1
      hx(2)=(x1-x2-x3+x4-x5+x6+x7-x8)
      hy(2)=(y1-y2-y3+y4-y5+y6+y7-y8)
      hz(2)=(z1-z2-z3+z4-z5+z6+z7-z8)
C mode 3
C 1 -1 1 -1 1 -1 1 -1
      hx(3)=(x1-x2+x3-x4+x5-x6+x7-x8)
      hy(3)=(y1-y2+y3-y4+y5-y6+y7-y8)
      hz(3)=(z1-z2+z3-z4+z5-z6+z7-z8)
C mode 4
C -1 1 -1 1 1 -1 1 -1
      hx(4)=(-x1+x2-x3+x4+x5-x6+x7-x8)
      hy(4)=(-y1+y2-y3+y4+y5-y6+y7-y8)
      hz(4)=(-z1+z2-z3+z4+z5-z6+z7-z8)
C-----------------------------------------------
      DO ip=1,ncell
       ir=irst(1,ip)
       is=irst(2,ip)
       it=irst(3,ip)
C
       ksi  = a_gauss(it,nsphdir)
       eta  = a_gauss(ir,nsphdir)
       zeta = a_gauss(is,nsphdir)
C
       wi = eight/np 
C
C Jacobian matrix
       jac1=cj1+hx(3)*eta+(hx(2)+hx(4)*eta)*zeta
       jac2=cj2+hy(3)*eta+(hy(2)+hy(4)*eta)*zeta
       jac3=cj3+hz(3)*eta+(hz(2)+hz(4)*eta)*zeta
C
       jac4=cj4+hx(1)*zeta+(hx(3)+hx(4)*zeta)*ksi
       jac5=cj5+hy(1)*zeta+(hy(3)+hy(4)*zeta)*ksi
       jac6=cj6+hz(1)*zeta+(hz(3)+hz(4)*zeta)*ksi
C
       jac7=cj7+hx(2)*ksi+(hx(1)+hx(4)*ksi)*eta
       jac8=cj8+hy(2)*ksi+(hy(1)+hy(4)*ksi)*eta
       jac9=cj9+hz(2)*ksi+(hz(1)+hz(4)*ksi)*eta
C
       jac_59_68=jac5*jac9-jac6*jac8
       jac_67_49=jac6*jac7-jac4*jac9
       jac_48_57=jac4*jac8-jac5*jac7
C
       det=one_over_512*(jac1*jac_59_68+jac2*jac_67_49+jac3*jac_48_57)
       vol= wi*det
       IF(det<zero)
     . CALL ancmsg(msgid=1038,
     .             msgtype=msgerror,
     .             anmode=aninfo,
     .             i1=ixs(nixs))
C
C Particle volume will be replaced by particle mass later (spinit3)
       spbuf(2,ip) =vol*rho
       spbuf(12,ip)=vol
C
      ENDDO
C-----------------------------------------------
      RETURN
      END SUBROUTINE soltosphv8
!||====================================================================
!||    soltosphv4    ../starter/source/elements/sph/soltosph.F
!||--- called by ------------------------------------------------------
!||    s4init3       ../starter/source/elements/solid/solide4/s4init3.F
!||--- calls      -----------------------------------------------------
!||    ancmsg        ../starter/source/output/message/message.f
!||--- uses       -----------------------------------------------------
!||    message_mod   ../starter/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE soltosphv4(
     .   NSPHDIR ,RHO     ,NCELL   ,X      ,SPBUF   ,
     .   IXS     )
C-----------------------------------------------
      USE message_mod
C-----------------------------------------------
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 NSPHDIR, NCELL, IXS(NIXS)
      my_real
     .        RHO, X(3,*), SPBUF(NSPBUF,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER IT,IP,N1, N2, N3, N4, NT, NP
C                                                                    
      my_real
     .   X1,X2,X3,X4,
     .   y1,y2,y3,y4,
     .   z1,z2,z3,z4,
     .   vol, det, wi,
     .   jac1 ,jac2 ,jac3 ,
     .   jac4 ,jac5 ,jac6 ,
     .   jac7 ,jac8 ,jac9 ,
     .   jac_59_68, jac_67_49, jac_48_57
C-----------------------------------------------
      np = 0
      nt = 0
      DO it=1,nsphdir
        nt=nt+it
        np=np+nt
      END DO
C-----------------------------------------------
      n1=ixs(2)
      x1=x(1,n1)
      y1=x(2,n1)
      z1=x(3,n1)
      n2=ixs(4)
      x2=x(1,n2)
      y2=x(2,n2)
      z2=x(3,n2)
      n3=ixs(7)
      x3=x(1,n3)
      y3=x(2,n3)
      z3=x(3,n3)
      n4=ixs(6)
      x4=x(1,n4)
      y4=x(2,n4)
      z4=x(3,n4)
C
C Jacobian matrix
      jac1=x1-x4
      jac2=y1-y4
      jac3=z1-z4
      jac4=x2-x4
      jac5=y2-y4
      jac6=z2-z4
      jac7=x3-x4
      jac8=y3-y4
      jac9=z3-z4
C
      jac_59_68=jac5*jac9-jac6*jac8
      jac_67_49=jac6*jac7-jac4*jac9
      jac_48_57=jac4*jac8-jac5*jac7
C
      det=jac1*jac_59_68+jac2*jac_67_49+jac3*jac_48_57
C-----------------------------------------------
      DO ip=1,ncell
C
       wi = one/(six*np)
C
       vol= wi*det
       IF(det<zero)
     . CALL ancmsg(msgid=1038,
     .             msgtype=msgerror,
     .             anmode=aninfo,
     .             i1=ixs(nixs))
C
C Particle volume will be replaced by particle mass later (spinit3)
       spbuf(2,ip) =vol*rho
       spbuf(12,ip)=vol
C
      ENDDO
C-----------------------------------------------
      RETURN
      END SUBROUTINE soltosphv4