OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
soltospha.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| soltospha ../engine/source/elements/sph/soltospha.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../engine/source/output/message/message.F
29!|| arret ../engine/source/system/arret.F
30!|| my_barrier ../engine/source/system/machine.F
31!|| startimeg ../engine/source/system/timer.F
32!|| stoptimeg ../engine/source/system/timer.F
33!||--- uses -----------------------------------------------------
34!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
35!|| message_mod ../engine/share/message_module/message_mod.F
36!|| soltosph_mod ../engine/share/modules/soltosph_mod.f
37!||====================================================================
38 SUBROUTINE soltospha(
39 1 ITASK ,V ,A ,MS ,PM ,
40 2 IPART ,IXS ,IPARTS ,KXSP ,IPARTSP ,
41 3 IRST ,SPBUF ,PARTSAV ,SOL2SPH ,IPARG ,
42 4 NGROUNC ,IGROUNC ,ELBUF_TAB,IGEO)
43C-----------------------------------------------
44 USE elbufdef_mod
45 USE soltosph_mod
46 USE message_mod
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51#include "comlock.inc"
52C-----------------------------------------------
53C C o m m o n B l o c k s
54C-----------------------------------------------
55#include "com01_c.inc"
56#include "com04_c.inc"
57#include "com08_c.inc"
58#include "param_c.inc"
59#include "sphcom.inc"
60#include "scr11_c.inc"
61#include "scr17_c.inc"
62#include "task_c.inc"
63#include "vect01_c.inc"
64C-----------------------------------------------
65C D u m m y A r g u m e n t s
66C-----------------------------------------------
67 INTEGER IXS(NIXS,*), KXSP(NISP,*), ITASK,
68 . IPARTSP(*), IRST(3,*), NGROUNC, IGROUNC(*),
69 . IPARTS(*), IPARG(NPARG,*), SOL2SPH(2,*), IPART(LIPART1,*),
70 . IGEO(NPROPGI,*)
71 my_real
72 . v(3,*), a(3,*), spbuf(nspbuf,*), ms(*), partsav(npsav,*),
73 . pm(npropm,*)
74 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER I, J, IR, IS, IT, KP, NP, N, INOD, IG, NG, MG, NEL,
79 . N1, N2, N3, N4, N5, N6, N7, N8, IPRT, IMAT,
80 . NELEM, OFFSET, NSPHDIR, KFT, IERROR
81C
82 my_real
83 . vx1,vx2,vx3,vx4,vx5,vx6,vx7,vx8,
84 . vy1,vy2,vy3,vy4,vy5,vy6,vy7,vy8,
85 . vz1,vz2,vz3,vz4,vz5,vz6,vz7,vz8,
86 . phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8,
87 . ksi, eta, zeta, wi,
88 . vxi, vyi, vzi, usdt, dt05,
89 . xmass_t, xmomt_t, ymomt_t, zmomt_t, encin_t,
90 . rho0, mass, vi2
91C-----
92 TYPE(g_bufel_) ,POINTER :: GBUF, GBUFSP
93 TYPE(L_BUFEL_) ,POINTER :: LBUF
94 TYPE(BUF_MAT_) ,POINTER :: MBUF
95C-----------------------------------------------
96 my_real
97 . a_gauss(9,9),a_gauss_tetra(9,9)
98 DATA a_gauss /
99 1 0. ,0. ,0. ,
100 1 0. ,0. ,0. ,
101 1 0. ,0. ,0. ,
102 2 -.5 ,0.5 ,0. ,
103 2 0. ,0. ,0. ,
104 2 0. ,0. ,0. ,
105 3 -.666666666666666,0. ,0.666666666666666,
106 3 0. ,0. ,0. ,
107 3 0. ,0. ,0. ,
108 4 -.75 ,-.25 ,0.25 ,
109 4 0.75 ,0. ,0. ,
110 4 0. ,0. ,0. ,
111 5 -.8 ,-.4 ,0. ,
112 5 0.4 ,0.8 ,0. ,
113 5 0. ,0. ,0. ,
114 6 -.833333333333333,-.5 ,-.166666666666666,
115 6 0.166666666666666,0.5 ,0.833333333333333,
116 6 0. ,0. ,0. ,
117 7 -.857142857142857,-.571428571428571,-.285714285714285,
118 7 0. ,0.285714285714285,0.571428571428571,
119 7 0.857142857142857,0. ,0. ,
120 8 -.875 ,-.625 ,-.375 ,
121 8 -.125 ,0.125 ,0.375,
122 8 0.625 ,0.875 ,0. ,
123 9 -.888888888888888,-.666666666666666,-.444444444444444,
124 9 -.222222222222222,0. ,0.222222222222222,
125 9 0.444444444444444,0.666666666666666,0.888888888888888/
126C-----------------------------------------------
127 DATA a_gauss_tetra /
128 1 0.250000000000000,0.000000000000000,0.000000000000000,
129 1 0.000000000000000,0.000000000000000,0.000000000000000,
130 1 0.000000000000000,0.000000000000000,0.000000000000000,
131 2 0.166666666666667,0.500000000000000,0.000000000000000,
132 2 0.000000000000000,0.000000000000000,0.000000000000000,
133 2 0.000000000000000,0.000000000000000,0.000000000000000,
134 3 0.125000000000000,0.375000000000000,0.625000000000000,
135 3 0.000000000000000,0.000000000000000,0.000000000000000,
136 3 0.000000000000000,0.000000000000000,0.000000000000000,
137 4 0.100000000000000,0.300000000000000,0.500000000000000,
138 4 0.700000000000000,0.000000000000000,0.000000000000000,
139 4 0.000000000000000,0.000000000000000,0.000000000000000,
140 5 0.083333333333333,0.250000000000000,0.416666666666667,
141 5 0.583333333333333,0.750000000000000,0.000000000000000,
142 5 0.000000000000000,0.000000000000000,0.000000000000000,
143 6 0.071428571428571,0.214285714285714,0.357142857142857,
144 6 0.500000000000000,0.642857142857143,0.785714285714286,
145 6 0.000000000000000,0.000000000000000,0.000000000000000,
146 7 0.062500000000000,0.187500000000000,0.312500000000000,
147 7 0.437500000000000,0.562500000000000,0.687500000000000,
148 7 0.812500000000000,0.000000000000000,0.000000000000000,
149 8 0.055555555555556,0.166666666666667,0.277777777777778,
150 8 0.388888888888889,0.500000000000000,0.611111111111111,
151 8 0.722222222222222,0.833333333333333,0.000000000000000,
152 9 0.050000000000000,0.150000000000000,0.250000000000000,
153 9 0.350000000000000,0.450000000000000,0.550000000000000,
154 9 0.650000000000000,0.750000000000000,0.850000000000000/
155C-----------------------------------------------
156 IF(itask==0)THEN
157 ALLOCATE(wpartsav(nthread,npsav,npart),stat=ierror)
158 IF(ierror/=0) THEN
159 CALL ancmsg(msgid=19,anmode=aninfo,
160 . c1='(SOLIDS to SPH)')
161 CALL arret(2)
162 ENDIF
163 END IF
164C
165 CALL my_barrier
166C
167 usdt=one/dt12
168 dt05=half*dt1
169 xmass_t=zero
170 xmomt_t=zero
171 ymomt_t=zero
172 zmomt_t=zero
173 encin_t=zero
174 DO iprt=itask+1,npart,nthread
175 DO j=1,npsav
176 DO i=1,nthread
177 wpartsav(i,j,iprt)=zero
178 END DO
179 END DO
180 END DO
181C
182!$OMP DO SCHEDULE(DYNAMIC,1)
183 DO ig = 1, ngrounc
184 ng = igrounc(ig)
185 IF(iparg(8,ng)==1)GOTO 250
186 IF (iddw>0) CALL startimeg(ng)
187 DO nelem = 1,iparg(2,ng),nvsiz
188 offset = nelem - 1
189 nel =iparg(2,ng)
190 nft =iparg(3,ng) + offset
191 iad =iparg(4,ng)
192 ity =iparg(5,ng)
193 ipartsph=iparg(69,ng)
194 lft=1
195 llt=min(nvsiz,nel-nelem+1)
196 IF(ity==1.AND.ipartsph/=0) THEN
197C-----------
198 gbuf => elbuf_tab(ng)%GBUF
199 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
200 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
201C-----
202C----TETRA---
203C-----
204 IF (iparg(28,ng)==4) THEN
205C-----
206 DO i=lft,llt
207 n=nft+i
208 IF(gbuf%OFF(i)==zero) cycle
209C
210C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
211 n1=ixs(2,n)
212 vx1=v(1,n1)+dt12*a(1,n1)
213 vy1=v(2,n1)+dt12*a(2,n1)
214 vz1=v(3,n1)+dt12*a(3,n1)
215 n2=ixs(4,n)
216 vx2=v(1,n2)+dt12*a(1,n2)
217 vy2=v(2,n2)+dt12*a(2,n2)
218 vz2=v(3,n2)+dt12*a(3,n2)
219 n3=ixs(7,n)
220 vx3=v(1,n3)+dt12*a(1,n3)
221 vy3=v(2,n3)+dt12*a(2,n3)
222 vz3=v(3,n3)+dt12*a(3,n3)
223 n4=ixs(6,n)
224 vx4=v(1,n4)+dt12*a(1,n4)
225 vy4=v(2,n4)+dt12*a(2,n4)
226 vz4=v(3,n4)+dt12*a(3,n4)
227C
228 nsphdir=igeo(37,ixs(10,n))
229C-----------------------------------------------
230 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
231 np =sol2sph(1,n)+kp
232 ir=irst(1,np-first_sphsol+1)
233 is=irst(2,np-first_sphsol+1)
234 it=irst(3,np-first_sphsol+1)
235 ksi = a_gauss_tetra(ir,nsphdir)
236 eta = a_gauss_tetra(is,nsphdir)
237 zeta = a_gauss_tetra(it,nsphdir)
238C
239 phi1=ksi
240 phi2=eta
241 phi3=zeta
242 phi4=1-ksi-eta-zeta
243C
244 vxi=phi1*vx1+phi2*vx2+phi3*vx3+phi4*vx4
245 vyi=phi1*vy1+phi2*vy2+phi3*vy3+phi4*vy4
246 vzi=phi1*vz1+phi2*vz2+phi3*vz3+phi4*vz4
247C
248 inod=kxsp(3,np)
249 a(1,inod)=(vxi-v(1,inod))*usdt
250 a(2,inod)=(vyi-v(2,inod))*usdt
251 a(3,inod)=(vzi-v(3,inod))*usdt
252 ENDDO
253 ENDDO
254C-----
255 ELSE
256C-----
257C----HEXA---
258C-----
259 DO i=lft,llt
260 n=nft+i
261 IF(gbuf%OFF(i)==zero) cycle
262C
263C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
264 n1=ixs(2,n)
265 vx1=v(1,n1)+dt12*a(1,n1)
266 vy1=v(2,n1)+dt12*a(2,n1)
267 vz1=v(3,n1)+dt12*a(3,n1)
268 n2=ixs(3,n)
269 vx2=v(1,n2)+dt12*a(1,n2)
270 vy2=v(2,n2)+dt12*a(2,n2)
271 vz2=v(3,n2)+dt12*a(3,n2)
272 n3=ixs(4,n)
273 vx3=v(1,n3)+dt12*a(1,n3)
274 vy3=v(2,n3)+dt12*a(2,n3)
275 vz3=v(3,n3)+dt12*a(3,n3)
276 n4=ixs(5,n)
277 vx4=v(1,n4)+dt12*a(1,n4)
278 vy4=v(2,n4)+dt12*a(2,n4)
279 vz4=v(3,n4)+dt12*a(3,n4)
280 n5=ixs(6,n)
281 vx5=v(1,n5)+dt12*a(1,n5)
282 vy5=v(2,n5)+dt12*a(2,n5)
283 vz5=v(3,n5)+dt12*a(3,n5)
284 n6=ixs(7,n)
285 vx6=v(1,n6)+dt12*a(1,n6)
286 vy6=v(2,n6)+dt12*a(2,n6)
287 vz6=v(3,n6)+dt12*a(3,n6)
288 n7=ixs(8,n)
289 vx7=v(1,n7)+dt12*a(1,n7)
290 vy7=v(2,n7)+dt12*a(2,n7)
291 vz7=v(3,n7)+dt12*a(3,n7)
292 n8=ixs(9,n)
293 vx8=v(1,n8)+dt12*a(1,n8)
294 vy8=v(2,n8)+dt12*a(2,n8)
295 vz8=v(3,n8)+dt12*a(3,n8)
296C
297C
298 nsphdir=nint((sol2sph(2,n)-sol2sph(1,n))**third)
299C-----------------------------------------------
300 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
301 np =sol2sph(1,n)+kp
302 ir=irst(1,np-first_sphsol+1)
303 is=irst(2,np-first_sphsol+1)
304 it=irst(3,np-first_sphsol+1)
305 ksi = a_gauss(ir,nsphdir)
306 eta = a_gauss(is,nsphdir)
307 zeta = a_gauss(it,nsphdir)
308C
309 phi1=(one-ksi)*(one-eta)*(one-zeta)
310 phi2=(one-ksi)*(one-eta)*(one+zeta)
311 phi3=(one+ksi)*(one-eta)*(one+zeta)
312 phi4=(one+ksi)*(one-eta)*(one-zeta)
313 phi5=(one-ksi)*(one+eta)*(one-zeta)
314 phi6=(one-ksi)*(one+eta)*(one+zeta)
315 phi7=(one+ksi)*(one+eta)*(one+zeta)
316 phi8=(one+ksi)*(one+eta)*(one-zeta)
317 vxi=one_over_8*(phi1*vx1+phi2*vx2+phi3*vx3+phi4*vx4+
318 . phi5*vx5+phi6*vx6+phi7*vx7+phi8*vx8)
319 vyi=one_over_8*(phi1*vy1+phi2*vy2+phi3*vy3+phi4*vy4+
320 . phi5*vy5+phi6*vy6+phi7*vy7+phi8*vy8)
321 vzi=one_over_8*(phi1*vz1+phi2*vz2+phi3*vz3+phi4*vz4+
322 . phi5*vz5+phi6*vz6+phi7*vz7+phi8*vz8)
323C
324 inod=kxsp(3,np)
325 a(1,inod)=(vxi-v(1,inod))*usdt
326 a(2,inod)=(vyi-v(2,inod))*usdt
327 a(3,inod)=(vzi-v(3,inod))*usdt
328 ENDDO
329 ENDDO
330C-----
331 ENDIF
332C-----
333 ENDIF
334 IF (iddw>0) CALL stoptimeg(ng)
335 END DO
336C--------
337 250 CONTINUE
338 END DO
339!$omp END DO
340C
341!$OMP DO SCHEDULE(DYNAMIC,1)
342 DO ig = 1, ngrounc
343 ng = igrounc(ig)
344 IF(iparg(8,ng)==1)GOTO 350
345 IF (iddw>0) CALL startimeg(ng)
346 DO nelem = 1,iparg(2,ng),nvsiz
347 offset = nelem - 1
348 nel =iparg(2,ng)
349 nft =iparg(3,ng) + offset
350 iad =iparg(4,ng)
351 ity =iparg(5,ng)
352 ipartsph=iparg(69,ng)
353 lft=1
354 llt=min(nvsiz,nel-nelem+1)
355 IF(ity==1.AND.ipartsph/=0) THEN
356C-----------
357 gbuf => elbuf_tab(ng)%GBUF
358 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
359 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
360C-----
361 DO i=lft,llt
362 n=nft+i
363 IF(gbuf%OFF(i)/=zero) THEN
364C
365C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
366C
367 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
368C
369 np=sol2sph(1,n)+kp
370 mg =mod(-kxsp(2,np),ngroup+1)
371C
372 inod=kxsp(3,np)
373 vxi=v(1,inod)+dt05*a(1,inod)
374 vyi=v(2,inod)+dt05*a(2,inod)
375 vzi=v(3,inod)+dt05*a(3,inod)
376 xmass_t=xmass_t-ms(inod)
377 xmomt_t=xmomt_t-ms(inod)*vxi
378 ymomt_t=ymomt_t-ms(inod)*vyi
379 zmomt_t=zmomt_t-ms(inod)*vzi
380 encin_t=encin_t
381 . -half*ms(inod)*(vxi*vxi+vyi*vyi+vzi*vzi)
382C
383C IPARG(8,MG)==1 <=> group was not computed
384C (there is no cloud active particle within the group)
385 IF(iparg(8,mg)==1)cycle
386C
387 kft=iparg(3,mg)
388 gbufsp => elbuf_tab(mg)%GBUF
389 iprt=ipartsp(np)
390 wpartsav(itask+1,1,iprt)=wpartsav(itask+1,1,iprt)
391 . -gbufsp%VOL(np-kft)*gbufsp%EINT(np-kft)
392 wpartsav(itask+1,2,iprt)=wpartsav(itask+1,2,iprt)
393 . -half*ms(inod)*(vxi*vxi+vyi*vyi+vzi*vzi)
394 wpartsav(itask+1,3,iprt)=wpartsav(itask+1,3,iprt)
395 . -ms(inod)*vxi
396 wpartsav(itask+1,4,iprt)=wpartsav(itask+1,4,iprt)
397 . -ms(inod)*vyi
398 wpartsav(itask+1,5,iprt)=wpartsav(itask+1,5,iprt)
399 . -ms(inod)*vzi
400 wpartsav(itask+1,6,iprt)=wpartsav(itask+1,6,iprt)
401 . -ms(inod)
402C
403 ENDDO
404 END IF
405 END DO
406 END IF
407 END DO
408 IF (iddw>0) CALL stoptimeg(ng)
409 350 CONTINUE
410 END DO
411!$OMP END DO
412C
413#include "lockon.inc"
414 encin=encin + encin_t
415 xmass=xmass + xmass_t
416 xmomt=xmomt + xmomt_t
417 ymomt=ymomt + ymomt_t
418 zmomt=zmomt + zmomt_t
419#include "lockoff.inc"
420 DO iprt=itask+1,npart,nthread
421 DO j=1,npsav
422 DO i=1,nthread
423 partsav(j,iprt)=partsav(j,iprt)+wpartsav(i,j,iprt)
424 END DO
425 END DO
426 END DO
427C
428 CALL my_barrier
429C
430 IF(itask==0) DEALLOCATE(wpartsav)
431C-----------------------------------------------
432 RETURN
433 END SUBROUTINE soltospha
my_real
#define my_real
Definition cppsort.cpp:32
startimeg
subroutine startimeg(ng)
Definition timer.F:1487
stoptimeg
subroutine stoptimeg(ng)
Definition timer.F:1535
min
#define min(a, b)
Definition macros.h:20
message_mod
Definition message_mod.F:1249
soltosph_mod
Definition soltosph_mod.F:32
soltospha
subroutine soltospha(itask, v, a, ms, pm, ipart, ixs, iparts, kxsp, ipartsp, irst, spbuf, partsav, sol2sph, iparg, ngrounc, igrounc, elbuf_tab, igeo)
Definition soltospha.F:43
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87
my_barrier
subroutine my_barrier
Definition machine.F:31