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!|| element_mod ../common_source/modules/elements/element_mod.F90
36!|| message_mod ../engine/share/message_module/message_mod.F
37!|| soltosph_mod ../engine/share/modules/soltosph_mod.F
38!||====================================================================
39 SUBROUTINE soltospha(
40 1 ITASK ,V ,A ,MS ,PM ,
41 2 IPART ,IXS ,IPARTS ,KXSP ,IPARTSP ,
42 3 IRST ,SPBUF ,PARTSAV ,SOL2SPH ,IPARG ,
43 4 NGROUNC ,IGROUNC ,ELBUF_TAB,IGEO)
44C-----------------------------------------------
45 USE elbufdef_mod
46 USE soltosph_mod
47 USE message_mod
48 use element_mod , only : nixs
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53#include "comlock.inc"
54C-----------------------------------------------
55C C o m m o n B l o c k s
56C-----------------------------------------------
57#include "com01_c.inc"
58#include "com04_c.inc"
59#include "com08_c.inc"
60#include "param_c.inc"
61#include "sphcom.inc"
62#include "scr11_c.inc"
63#include "scr17_c.inc"
64#include "task_c.inc"
65#include "vect01_c.inc"
66C-----------------------------------------------
67C D u m m y A r g u m e n t s
68C-----------------------------------------------
69 INTEGER IXS(NIXS,*), KXSP(NISP,*), ITASK,
70 . IPARTSP(*), IRST(3,*), NGROUNC, IGROUNC(*),
71 . IPARTS(*), IPARG(NPARG,*), SOL2SPH(2,*), IPART(LIPART1,*),
72 . IGEO(NPROPGI,*)
73 my_real
74 . v(3,*), a(3,*), spbuf(nspbuf,*), ms(*), partsav(npsav,*),
75 . pm(npropm,*)
76 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
77C-----------------------------------------------
78C L o c a l V a r i a b l e s
79C-----------------------------------------------
80 INTEGER I, J, IR, IS, IT, KP, NP, N, INOD, IG, NG, MG, NEL,
81 . N1, N2, N3, N4, N5, N6, N7, N8, IPRT,
82 . NELEM, OFFSET, NSPHDIR, KFT, IERROR
83C
84 my_real
85 . vx1,vx2,vx3,vx4,vx5,vx6,vx7,vx8,
86 . vy1,vy2,vy3,vy4,vy5,vy6,vy7,vy8,
87 . vz1,vz2,vz3,vz4,vz5,vz6,vz7,vz8,
88 . phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8,
89 . ksi, eta, zeta,
90 . vxi, vyi, vzi, usdt, dt05,
91 . xmass_t, xmomt_t, ymomt_t, zmomt_t, encin_t
92C-----
93 TYPE(g_bufel_) ,POINTER :: GBUF, GBUFSP
94 TYPE(L_BUFEL_) ,POINTER :: LBUF
95 TYPE(BUF_MAT_) ,POINTER :: MBUF
96C-----------------------------------------------
97 my_real
98 . a_gauss(9,9),a_gauss_tetra(9,9)
99 DATA a_gauss /
100 1 0. ,0. ,0. ,
101 1 0. ,0. ,0. ,
102 1 0. ,0. ,0. ,
103 2 -.5 ,0.5 ,0. ,
104 2 0. ,0. ,0. ,
105 2 0. ,0. ,0. ,
106 3 -.666666666666666,0. ,0.666666666666666,
107 3 0. ,0. ,0. ,
108 3 0. ,0. ,0. ,
109 4 -.75 ,-.25 ,0.25 ,
110 4 0.75 ,0. ,0. ,
111 4 0. ,0. ,0. ,
112 5 -.8 ,-.4 ,0. ,
113 5 0.4 ,0.8 ,0. ,
114 5 0. ,0. ,0. ,
115 6 -.833333333333333,-.5 ,-.166666666666666,
116 6 0.166666666666666,0.5 ,0.833333333333333,
117 6 0. ,0. ,0. ,
118 7 -.857142857142857,-.571428571428571,-.285714285714285,
119 7 0. ,0.285714285714285,0.571428571428571,
120 7 0.857142857142857,0. ,0. ,
121 8 -.875 ,-.625 ,-.375 ,
122 8 -.125 ,0.125 ,0.375,
123 8 0.625 ,0.875 ,0. ,
124 9 -.888888888888888,-.666666666666666,-.444444444444444,
125 9 -.222222222222222,0. ,0.222222222222222,
126 9 0.444444444444444,0.666666666666666,0.888888888888888/
127C-----------------------------------------------
128 DATA a_gauss_tetra /
129 1 0.250000000000000,0.000000000000000,0.000000000000000,
130 1 0.000000000000000,0.000000000000000,0.000000000000000,
131 1 0.000000000000000,0.000000000000000,0.000000000000000,
132 2 0.166666666666667,0.500000000000000,0.000000000000000,
133 2 0.000000000000000,0.000000000000000,0.000000000000000,
134 2 0.000000000000000,0.000000000000000,0.000000000000000,
135 3 0.125000000000000,0.375000000000000,0.625000000000000,
136 3 0.000000000000000,0.000000000000000,0.000000000000000,
137 3 0.000000000000000,0.000000000000000,0.000000000000000,
138 4 0.100000000000000,0.300000000000000,0.500000000000000,
139 4 0.700000000000000,0.000000000000000,0.000000000000000,
140 4 0.000000000000000,0.000000000000000,0.000000000000000,
141 5 0.083333333333333,0.250000000000000,0.416666666666667,
142 5 0.583333333333333,0.750000000000000,0.000000000000000,
143 5 0.000000000000000,0.000000000000000,0.000000000000000,
144 6 0.071428571428571,0.214285714285714,0.357142857142857,
145 6 0.500000000000000,0.642857142857143,0.785714285714286,
146 6 0.000000000000000,0.000000000000000,0.000000000000000,
147 7 0.062500000000000,0.187500000000000,0.312500000000000,
148 7 0.437500000000000,0.562500000000000,0.687500000000000,
149 7 0.812500000000000,0.000000000000000,0.000000000000000,
150 8 0.055555555555556,0.166666666666667,0.277777777777778,
151 8 0.388888888888889,0.500000000000000,0.611111111111111,
152 8 0.722222222222222,0.833333333333333,0.000000000000000,
153 9 0.050000000000000,0.150000000000000,0.250000000000000,
154 9 0.350000000000000,0.450000000000000,0.550000000000000,
155 9 0.650000000000000,0.750000000000000,0.850000000000000/
156C-----------------------------------------------
157 IF(itask==0)THEN
158 ALLOCATE(wpartsav(nthread,npsav,npart),stat=ierror)
159 IF(ierror/=0) THEN
160 CALL ancmsg(msgid=19,anmode=aninfo,
161 . c1='(SOLIDS to SPH)')
162 CALL arret(2)
163 ENDIF
164 END IF
165C
166 CALL my_barrier
167C
168 usdt=one/dt12
169 dt05=half*dt1
170 xmass_t=zero
171 xmomt_t=zero
172 ymomt_t=zero
173 zmomt_t=zero
174 encin_t=zero
175 DO iprt=itask+1,npart,nthread
176 DO j=1,npsav
177 DO i=1,nthread
178 wpartsav(i,j,iprt)=zero
179 END DO
180 END DO
181 END DO
182C
183!$OMP DO SCHEDULE(DYNAMIC,1)
184 DO ig = 1, ngrounc
185 ng = igrounc(ig)
186 IF(iparg(8,ng)==1)GOTO 250
187 IF (iddw>0) CALL startimeg(ng)
188 DO nelem = 1,iparg(2,ng),nvsiz
189 offset = nelem - 1
190 nel =iparg(2,ng)
191 nft =iparg(3,ng) + offset
192 iad =iparg(4,ng)
193 ity =iparg(5,ng)
194 ipartsph=iparg(69,ng)
195 lft=1
196 llt=min(nvsiz,nel-nelem+1)
197 IF(ity==1.AND.ipartsph/=0) THEN
198C-----------
199 gbuf => elbuf_tab(ng)%GBUF
200 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
201 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
202C-----
203C----TETRA---
204C-----
205 IF (iparg(28,ng)==4) THEN
206C-----
207 DO i=lft,llt
208 n=nft+i
209 IF(gbuf%OFF(i)==zero) cycle
210C
211C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
212 n1=ixs(2,n)
213 vx1=v(1,n1)+dt12*a(1,n1)
214 vy1=v(2,n1)+dt12*a(2,n1)
215 vz1=v(3,n1)+dt12*a(3,n1)
216 n2=ixs(4,n)
217 vx2=v(1,n2)+dt12*a(1,n2)
218 vy2=v(2,n2)+dt12*a(2,n2)
219 vz2=v(3,n2)+dt12*a(3,n2)
220 n3=ixs(7,n)
221 vx3=v(1,n3)+dt12*a(1,n3)
222 vy3=v(2,n3)+dt12*a(2,n3)
223 vz3=v(3,n3)+dt12*a(3,n3)
224 n4=ixs(6,n)
225 vx4=v(1,n4)+dt12*a(1,n4)
226 vy4=v(2,n4)+dt12*a(2,n4)
227 vz4=v(3,n4)+dt12*a(3,n4)
228C
229 nsphdir=igeo(37,ixs(10,n))
230C-----------------------------------------------
231 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
232 np =sol2sph(1,n)+kp
233 ir=irst(1,np-first_sphsol+1)
234 is=irst(2,np-first_sphsol+1)
235 it=irst(3,np-first_sphsol+1)
236 ksi = a_gauss_tetra(ir,nsphdir)
237 eta = a_gauss_tetra(is,nsphdir)
238 zeta = a_gauss_tetra(it,nsphdir)
239C
240 phi1=ksi
241 phi2=eta
242 phi3=zeta
243 phi4=1-ksi-eta-zeta
244C
245 vxi=phi1*vx1+phi2*vx2+phi3*vx3+phi4*vx4
246 vyi=phi1*vy1+phi2*vy2+phi3*vy3+phi4*vy4
247 vzi=phi1*vz1+phi2*vz2+phi3*vz3+phi4*vz4
248C
249 inod=kxsp(3,np)
250 a(1,inod)=(vxi-v(1,inod))*usdt
251 a(2,inod)=(vyi-v(2,inod))*usdt
252 a(3,inod)=(vzi-v(3,inod))*usdt
253 ENDDO
254 ENDDO
255C-----
256 ELSE
257C-----
258C----HEXA---
259C-----
260 DO i=lft,llt
261 n=nft+i
262 IF(gbuf%OFF(i)==zero) cycle
263C
264C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
265 n1=ixs(2,n)
266 vx1=v(1,n1)+dt12*a(1,n1)
267 vy1=v(2,n1)+dt12*a(2,n1)
268 vz1=v(3,n1)+dt12*a(3,n1)
269 n2=ixs(3,n)
270 vx2=v(1,n2)+dt12*a(1,n2)
271 vy2=v(2,n2)+dt12*a(2,n2)
272 vz2=v(3,n2)+dt12*a(3,n2)
273 n3=ixs(4,n)
274 vx3=v(1,n3)+dt12*a(1,n3)
275 vy3=v(2,n3)+dt12*a(2,n3)
276 vz3=v(3,n3)+dt12*a(3,n3)
277 n4=ixs(5,n)
278 vx4=v(1,n4)+dt12*a(1,n4)
279 vy4=v(2,n4)+dt12*a(2,n4)
280 vz4=v(3,n4)+dt12*a(3,n4)
281 n5=ixs(6,n)
282 vx5=v(1,n5)+dt12*a(1,n5)
283 vy5=v(2,n5)+dt12*a(2,n5)
284 vz5=v(3,n5)+dt12*a(3,n5)
285 n6=ixs(7,n)
286 vx6=v(1,n6)+dt12*a(1,n6)
287 vy6=v(2,n6)+dt12*a(2,n6)
288 vz6=v(3,n6)+dt12*a(3,n6)
289 n7=ixs(8,n)
290 vx7=v(1,n7)+dt12*a(1,n7)
291 vy7=v(2,n7)+dt12*a(2,n7)
292 vz7=v(3,n7)+dt12*a(3,n7)
293 n8=ixs(9,n)
294 vx8=v(1,n8)+dt12*a(1,n8)
295 vy8=v(2,n8)+dt12*a(2,n8)
296 vz8=v(3,n8)+dt12*a(3,n8)
297C
298C
299 nsphdir=nint((sol2sph(2,n)-sol2sph(1,n))**third)
300C-----------------------------------------------
301 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
302 np =sol2sph(1,n)+kp
303 ir=irst(1,np-first_sphsol+1)
304 is=irst(2,np-first_sphsol+1)
305 it=irst(3,np-first_sphsol+1)
306 ksi = a_gauss(ir,nsphdir)
307 eta = a_gauss(is,nsphdir)
308 zeta = a_gauss(it,nsphdir)
309C
310 phi1=(one-ksi)*(one-eta)*(one-zeta)
311 phi2=(one-ksi)*(one-eta)*(one+zeta)
312 phi3=(one+ksi)*(one-eta)*(one+zeta)
313 phi4=(one+ksi)*(one-eta)*(one-zeta)
314 phi5=(one-ksi)*(one+eta)*(one-zeta)
315 phi6=(one-ksi)*(one+eta)*(one+zeta)
316 phi7=(one+ksi)*(one+eta)*(one+zeta)
317 phi8=(one+ksi)*(one+eta)*(one-zeta)
318 vxi=one_over_8*(phi1*vx1+phi2*vx2+phi3*vx3+phi4*vx4+
319 . phi5*vx5+phi6*vx6+phi7*vx7+phi8*vx8)
320 vyi=one_over_8*(phi1*vy1+phi2*vy2+phi3*vy3+phi4*vy4+
321 . phi5*vy5+phi6*vy6+phi7*vy7+phi8*vy8)
322 vzi=one_over_8*(phi1*vz1+phi2*vz2+phi3*vz3+phi4*vz4+
323 . phi5*vz5+phi6*vz6+phi7*vz7+phi8*vz8)
324C
325 inod=kxsp(3,np)
326 a(1,inod)=(vxi-v(1,inod))*usdt
327 a(2,inod)=(vyi-v(2,inod))*usdt
328 a(3,inod)=(vzi-v(3,inod))*usdt
329 ENDDO
330 ENDDO
331C-----
332 ENDIF
333C-----
334 ENDIF
335 IF (iddw>0) CALL stoptimeg(ng)
336 END DO
337C--------
338 250 CONTINUE
339 END DO
340!$OMP END DO
341C
342!$OMP DO SCHEDULE(DYNAMIC,1)
343 DO ig = 1, ngrounc
344 ng = igrounc(ig)
345 IF(iparg(8,ng)==1)GOTO 350
346 IF (iddw>0) CALL startimeg(ng)
347 DO nelem = 1,iparg(2,ng),nvsiz
348 offset = nelem - 1
349 nel =iparg(2,ng)
350 nft =iparg(3,ng) + offset
351 iad =iparg(4,ng)
352 ity =iparg(5,ng)
353 ipartsph=iparg(69,ng)
354 lft=1
355 llt=min(nvsiz,nel-nelem+1)
356 IF(ity==1.AND.ipartsph/=0) THEN
357C-----------
358 gbuf => elbuf_tab(ng)%GBUF
359 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
360 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
361C-----
362 DO i=lft,llt
363 n=nft+i
364 IF(gbuf%OFF(i)/=zero) THEN
365C
366C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
367C
368 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
369C
370 np=sol2sph(1,n)+kp
371 IF (kxsp(2,np)/=0) THEN
372C
373 mg =mod(-kxsp(2,np),ngroup+1)
374C
375 inod=kxsp(3,np)
376 vxi=v(1,inod)+dt05*a(1,inod)
377 vyi=v(2,inod)+dt05*a(2,inod)
378 vzi=v(3,inod)+dt05*a(3,inod)
379 xmass_t=xmass_t-ms(inod)
380 xmomt_t=xmomt_t-ms(inod)*vxi
381 ymomt_t=ymomt_t-ms(inod)*vyi
382 zmomt_t=zmomt_t-ms(inod)*vzi
383 encin_t=encin_t
384 . -half*ms(inod)*(vxi*vxi+vyi*vyi+vzi*vzi)
385C
386C IPARG(8,MG)==1 <=> group was not computed
387C (there is no cloud active particle within the group)
388 IF(iparg(8,mg)==1)cycle
389C
390 kft=iparg(3,mg)
391 gbufsp => elbuf_tab(mg)%GBUF
392 iprt=ipartsp(np)
393 wpartsav(itask+1,1,iprt)=wpartsav(itask+1,1,iprt)
394 . -gbufsp%VOL(np-kft)*gbufsp%EINT(np-kft)
395 wpartsav(itask+1,2,iprt)=wpartsav(itask+1,2,iprt)
396 . -half*ms(inod)*(vxi*vxi+vyi*vyi+vzi*vzi)
397 wpartsav(itask+1,3,iprt)=wpartsav(itask+1,3,iprt)
398 . -ms(inod)*vxi
399 wpartsav(itask+1,4,iprt)=wpartsav(itask+1,4,iprt)
400 . -ms(inod)*vyi
401 wpartsav(itask+1,5,iprt)=wpartsav(itask+1,5,iprt)
402 . -ms(inod)*vzi
403 wpartsav(itask+1,6,iprt)=wpartsav(itask+1,6,iprt)
404 . -ms(inod)
405C
406 ENDIF
407C
408 ENDDO
409 END IF
410 END DO
411 END IF
412 END DO
413 IF (iddw>0) CALL stoptimeg(ng)
414 350 CONTINUE
415 END DO
416!$OMP END DO
417C
418#include "lockon.inc"
419 encin=encin + encin_t
420 xmass=xmass + xmass_t
421 xmomt=xmomt + xmomt_t
422 ymomt=ymomt + ymomt_t
423 zmomt=zmomt + zmomt_t
424#include "lockoff.inc"
425 DO iprt=itask+1,npart,nthread
426 DO j=1,npsav
427 DO i=1,nthread
428 partsav(j,iprt)=partsav(j,iprt)+wpartsav(i,j,iprt)
429 END DO
430 END DO
431 END DO
432C
433 CALL my_barrier
434C
435 IF(itask==0) DEALLOCATE(wpartsav)
436C-----------------------------------------------
437 RETURN
438 END SUBROUTINE soltospha
my_real
#define my_real
Definition cppsort.cpp:32
startimeg
subroutine startimeg(ng)
Definition timer.F:1371
stoptimeg
subroutine stoptimeg(ng)
Definition timer.F:1419
min
#define min(a, b)
Definition macros.h:20
message_mod
Definition message_mod.F:1257
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:44
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:895
arret
subroutine arret(nn)
Definition arret.F:86
my_barrier
subroutine my_barrier
Definition machine.F:31