OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
splissv.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "vect01_c.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com06_c.inc"
#include "com08_c.inc"
#include "sphcom.inc"
#include "param_c.inc"
#include "scr17_c.inc"
#include "task_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine splissv (x, v, ms, a, spbuf, wa, itab, kxsp, ixsp, nod2sp, d, ispsym, xspsym, vspsym, bufmat, bufgeo, npc, pld, pm, geo, ispcond, xframe, waspsym, ipartsp, partsav, wacomp, wsmcomp, waspact, ipart, itask, sph2sol, sol2sph, irst, ixs, iparg, ngrounc, igrounc, elbuf_tab, iad_elem, fr_elem, igeo, sol2sph_typ, sph_work)

Function/Subroutine Documentation

◆ splissv()

subroutine splissv ( x,
v,
ms,
a,
spbuf,
wa,
integer, dimension(*) itab,
integer, dimension(nisp,*) kxsp,
integer, dimension(kvoisph,*) ixsp,
integer, dimension(*) nod2sp,
d,
integer, dimension(nspcond,*) ispsym,
xspsym,
vspsym,
bufmat,
bufgeo,
integer, dimension(*) npc,
pld,
pm,
geo,
integer, dimension(nispcond,*) ispcond,
xframe,
waspsym,
integer, dimension(*) ipartsp,
partsav,
wacomp,
wsmcomp,
integer, dimension(*) waspact,
integer, dimension(lipart1,*) ipart,
integer itask,
integer, dimension(*) sph2sol,
integer, dimension(2,*) sol2sph,
integer, dimension(3,*) irst,
integer, dimension(nixs,*) ixs,
integer, dimension(nparg,*) iparg,
integer ngrounc,
integer, dimension(*) igrounc,
type (elbuf_struct_), dimension(ngroup), target elbuf_tab,
integer, dimension(2,*) iad_elem,
integer, dimension(*) fr_elem,
integer, dimension(npropgi,*) igeo,
integer, dimension(*) sol2sph_typ,
type (sph_work_) sph_work )

Definition at line 43 of file splissv.F.

53C-----------------------------------------------
54C M o d u l e s
55C-----------------------------------------------
56 USE my_alloc_mod
57 USE elbufdef_mod
58 USE message_mod
59 USE sphbox
60 USE sph_work_mod
61 use element_mod , only : nixs
62C-----------------------------------------------
63C I m p l i c i t T y p e s
64C-----------------------------------------------
65#include "implicit_f.inc"
66#include "comlock.inc"
67C-----------------------------------------------
68C C o m m o n B l o c k s
69C-----------------------------------------------
70#include "vect01_c.inc"
71#include "com01_c.inc"
72#include "com04_c.inc"
73#include "com06_c.inc"
74#include "com08_c.inc"
75#include "sphcom.inc"
76#include "param_c.inc"
77#include "scr17_c.inc"
78#include "task_c.inc"
79C-----------------------------------------------
80C D u m m y A r g u m e n t s
81C-----------------------------------------------
82 INTEGER KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),ITAB(*),
83 . ISPSYM(NSPCOND,*),NPC(*),ISPCOND(NISPCOND,*),
84 . IPARTSP(*),WASPACT(*),IPART(LIPART1,*), ITASK,
85 . SPH2SOL(*),IXS(NIXS,*),IRST(3,*),SOL2SPH(2,*),
86 . IPARG(NPARG,*), NGROUNC, IGROUNC(*),
87 . IAD_ELEM(2,*),FR_ELEM(*),IGEO(NPROPGI,*),SOL2SPH_TYP(*)
88C REAL
89 my_real
90 . x(3,*) ,v(3,*) ,ms(*) ,
91 . a(3,*) ,spbuf(nspbuf,*) ,wa(*),
92 . d(3,*) ,xspsym(3,*) ,vspsym(3,*),
93 . pm(npropm,*), geo(npropg,*),bufmat(*),bufgeo(*),pld(*),
94 . xframe(nxframe,*), waspsym(3,*), partsav(npsav,*),
95 . wacomp(16,*), wsmcomp(6,*)
96 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
97 TYPE (SPH_WORK_) :: SPH_WORK
98C-----------------------------------------------
99C L o c a l V a r i a b l e s
100C-----------------------------------------------
101 INTEGER N,INOD,JNOD,J,NVOIS,M,NN,
102 . IPROP,IMAT,I,
103 . NVOISS,SM,NC,JS,
104 . IS,IC,ISLIDE,IPRT,NS,MYADRN,
105 . NSOL, N1, N2, N3, N4, N5, N6, N7, N8,
106 . IR, IT, NSPHDIR, KP, NP,NELEM, OFFSET, NEL, IG, NG,
107 . II, II1, SZ, LENR
108 my_real
109 . xi,yi,zi,di,rhoi,xj,yj,zj,dj,rhoj,dij,
110 . vxi,vyi,vzi,vxj,vyj,vzj,
111 .
112 . wght,
113 .
114 . alpci,fact,faci,
115 . wax,way,waz,axi,ayi,azi,axj,ayj,azj,an,
116 . vv,kv,ehourt,dtinv,
117 . ox,oy,oz,nx,ny,nz,axs,ays,azs,
118 . alphai,betaxi,betayi,betazi,betai,
119 . alphaj,betaxj,betayj,betazj,betaj,unm,
120 . vx1,vx2,vx3,vx4,vx5,vx6,vx7,vx8,
121 . vy1,vy2,vy3,vy4,vy5,vy6,vy7,vy8,
122 . vz1,vz2,vz3,vz4,vz5,vz6,vz7,vz8,usdt,
123 . phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8,
124 . ksi, eta, zeta
125C-----------------------------------------------
126 my_real get_u_geo
127 EXTERNAL get_u_geo
128C-----
129 TYPE(G_BUFEL_) ,POINTER :: GBUF
130 TYPE(L_BUFEL_) ,POINTER :: LBUF
131 TYPE(BUF_MAT_) ,POINTER :: MBUF
132C-----------------------------------------------
133 my_real
134 . a_gauss(9,9),a_gauss_tetra(9,9)
135 DATA a_gauss /
136 1 0. ,0. ,0. ,
137 1 0. ,0. ,0. ,
138 1 0. ,0. ,0. ,
139 2 -.5 ,0.5 ,0. ,
140 2 0. ,0. ,0. ,
141 2 0. ,0. ,0. ,
142 3 -.666666666666666,0. ,0.666666666666666,
143 3 0. ,0. ,0. ,
144 3 0. ,0. ,0. ,
145 4 -.75 ,-.25 ,0.25 ,
146 4 0.75 ,0. ,0. ,
147 4 0. ,0. ,0. ,
148 5 -.8 ,-.4 ,0. ,
149 5 0.4 ,0.8 ,0. ,
150 5 0. ,0. ,0. ,
151 6 -.833333333333333,-.5 ,-.166666666666666,
152 6 0.166666666666666,0.5 ,0.833333333333333,
153 6 0. ,0. ,0. ,
154 7 -.857142857142857,-.571428571428571,-.285714285714285,
155 7 0. ,0.285714285714285,0.571428571428571,
156 7 0.857142857142857,0. ,0. ,
157 8 -.875 ,-.625 ,-.375 ,
158 8 -.125 ,0.125 ,0.375,
159 8 0.625 ,0.875 ,0. ,
160 9 -.888888888888888,-.666666666666666,-.444444444444444,
161 9 -.222222222222222,0. ,0.222222222222222,
162 9 0.444444444444444,0.666666666666666,0.888888888888888/
163C-----------------------------------------------
164 DATA a_gauss_tetra /
165 1 0.250000000000000,0.000000000000000,0.000000000000000,
166 1 0.000000000000000,0.000000000000000,0.000000000000000,
167 1 0.000000000000000,0.000000000000000,0.000000000000000,
168 2 0.166666666666667,0.500000000000000,0.000000000000000,
169 2 0.000000000000000,0.000000000000000,0.000000000000000,
170 2 0.000000000000000,0.000000000000000,0.000000000000000,
171 3 0.125000000000000,0.375000000000000,0.625000000000000,
172 3 0.000000000000000,0.000000000000000,0.000000000000000,
173 3 0.000000000000000,0.000000000000000,0.000000000000000,
174 4 0.100000000000000,0.300000000000000,0.500000000000000,
175 4 0.700000000000000,0.000000000000000,0.000000000000000,
176 4 0.000000000000000,0.000000000000000,0.000000000000000,
177 5 0.083333333333333,0.250000000000000,0.416666666666667,
178 5 0.583333333333333,0.750000000000000,0.000000000000000,
179 5 0.000000000000000,0.000000000000000,0.000000000000000,
180 6 0.071428571428571,0.214285714285714,0.357142857142857,
181 6 0.500000000000000,0.642857142857143,0.785714285714286,
182 6 0.000000000000000,0.000000000000000,0.000000000000000,
183 7 0.062500000000000,0.187500000000000,0.312500000000000,
184 7 0.437500000000000,0.562500000000000,0.687500000000000,
185 7 0.812500000000000,0.000000000000000,0.000000000000000,
186 8 0.055555555555556,0.166666666666667,0.277777777777778,
187 8 0.388888888888889,0.500000000000000,0.611111111111111,
188 8 0.722222222222222,0.833333333333333,0.000000000000000,
189 9 0.050000000000000,0.150000000000000,0.250000000000000,
190 9 0.350000000000000,0.450000000000000,0.550000000000000,
191 9 0.650000000000000,0.750000000000000,0.850000000000000/
192C-----------------------------------------------
193C
194 IF(sol2sph_flag/=0)THEN
195 IF(itask==0)THEN
196C
197 sph_work%A6(1:6,1:3,1:numnod) = zero
198 IF (ALLOCATED(sph_work%AS))DEALLOCATE(sph_work%AS)
199 CALL my_alloc(sph_work%AS,3,8*nsphact)
200 sph_work%AS= zero
201
202 IF (ALLOCATED(sph_work%AS))DEALLOCATE(sph_work%AS6)
203 CALL my_alloc(sph_work%AS6,6,3,8*nsphact)
204 sph_work%AS6= zero
205
206 sph_work%ITAG= 0
207 END IF
208 ENDIF
209C
210 IF (numsph > 0) THEN
211C
212C-----------------------------------------------
213C mean accelerations on cloud active particles
214C a optimiser : traiter NSPHACT only !!!
215C-----------------------------------------------
216
217 IF (nsphsol > 0) THEN
218 usdt=one/dt12
219!$OMP DO SCHEDULE(DYNAMIC,1)
220 DO ig = 1, ngrounc
221 ng = igrounc(ig)
222 IF(iparg(8,ng)==1)GOTO 250
223 IF (iddw>0) CALL startimeg(ng)
224 DO nelem = 1,iparg(2,ng),nvsiz
225 offset = nelem - 1
226 nel =iparg(2,ng)
227 nft =iparg(3,ng) + offset
228 iad =iparg(4,ng)
229 ity =iparg(5,ng)
230 ipartsph=iparg(69,ng)
231 lft=1
232 llt=min(nvsiz,nel-nelem+1)
233 IF(ity==1.AND.ipartsph/=0) THEN
234C-----------
235 gbuf => elbuf_tab(ng)%GBUF
236 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
237 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
238C-----
239 IF (iparg(28,ng)==4) THEN
240C-----
241C----TETRA---
242C-----
243 DO i=lft,llt
244 n=nft+i
245 IF(gbuf%OFF(i)==zero) cycle
246C
247C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
248 n1=ixs(2,n)
249 vx1=v(1,n1)+dt12*a(1,n1)/ms(n1)
250 vy1=v(2,n1)+dt12*a(2,n1)/ms(n1)
251 vz1=v(3,n1)+dt12*a(3,n1)/ms(n1)
252 n2=ixs(4,n)
253 vx2=v(1,n2)+dt12*a(1,n2)/ms(n2)
254 vy2=v(2,n2)+dt12*a(2,n2)/ms(n2)
255 vz2=v(3,n2)+dt12*a(3,n2)/ms(n2)
256 n3=ixs(7,n)
257 vx3=v(1,n3)+dt12*a(1,n3)/ms(n3)
258 vy3=v(2,n3)+dt12*a(2,n3)/ms(n3)
259 vz3=v(3,n3)+dt12*a(3,n3)/ms(n3)
260 n4=ixs(6,n)
261 vx4=v(1,n4)+dt12*a(1,n4)/ms(n4)
262 vy4=v(2,n4)+dt12*a(2,n4)/ms(n4)
263 vz4=v(3,n4)+dt12*a(3,n4)/ms(n4)
264C
265 nsphdir=igeo(37,ixs(10,n))
266C-----------------------------------------------
267 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
268 np =sol2sph(1,n)+kp
269 ir=irst(1,np-first_sphsol+1)
270 is=irst(2,np-first_sphsol+1)
271 it=irst(3,np-first_sphsol+1)
272C
273 ksi = a_gauss_tetra(ir,nsphdir)
274 eta = a_gauss_tetra(is,nsphdir)
275 zeta = a_gauss_tetra(it,nsphdir)
276C
277 phi1=ksi
278 phi2=eta
279 phi3=zeta
280 phi4=1-ksi-eta-zeta
281C
282 vxi=phi1*vx1+phi2*vx2+phi3*vx3+phi4*vx4
283 vyi=phi1*vy1+phi2*vy2+phi3*vy3+phi4*vy4
284 vzi=phi1*vz1+phi2*vz2+phi3*vz3+phi4*vz4
285C
286 inod=kxsp(3,np)
287 a(1,inod)=ms(inod)*(vxi-v(1,inod))*usdt
288 a(2,inod)=ms(inod)*(vyi-v(2,inod))*usdt
289 a(3,inod)=ms(inod)*(vzi-v(3,inod))*usdt
290 ENDDO
291 ENDDO
292C-----
293 ELSE
294C-----
295C----HEXA---
296C-----
297 DO i=lft,llt
298 n=nft+i
299 IF(gbuf%OFF(i)==zero) cycle
300C
301C SOL2SPH(1,N)+1<=I<=SOLSPH(2,N) <=> N==SPH2SOL(I)
302 n1=ixs(2,n)
303 vx1=v(1,n1)+dt12*a(1,n1)/ms(n1)
304 vy1=v(2,n1)+dt12*a(2,n1)/ms(n1)
305 vz1=v(3,n1)+dt12*a(3,n1)/ms(n1)
306 n2=ixs(3,n)
307 vx2=v(1,n2)+dt12*a(1,n2)/ms(n2)
308 vy2=v(2,n2)+dt12*a(2,n2)/ms(n2)
309 vz2=v(3,n2)+dt12*a(3,n2)/ms(n2)
310 n3=ixs(4,n)
311 vx3=v(1,n3)+dt12*a(1,n3)/ms(n3)
312 vy3=v(2,n3)+dt12*a(2,n3)/ms(n3)
313 vz3=v(3,n3)+dt12*a(3,n3)/ms(n3)
314 n4=ixs(5,n)
315 vx4=v(1,n4)+dt12*a(1,n4)/ms(n4)
316 vy4=v(2,n4)+dt12*a(2,n4)/ms(n4)
317 vz4=v(3,n4)+dt12*a(3,n4)/ms(n4)
318 n5=ixs(6,n)
319 vx5=v(1,n5)+dt12*a(1,n5)/ms(n5)
320 vy5=v(2,n5)+dt12*a(2,n5)/ms(n5)
321 vz5=v(3,n5)+dt12*a(3,n5)/ms(n5)
322 n6=ixs(7,n)
323 vx6=v(1,n6)+dt12*a(1,n6)/ms(n6)
324 vy6=v(2,n6)+dt12*a(2,n6)/ms(n6)
325 vz6=v(3,n6)+dt12*a(3,n6)/ms(n6)
326 n7=ixs(8,n)
327 vx7=v(1,n7)+dt12*a(1,n7)/ms(n7)
328 vy7=v(2,n7)+dt12*a(2,n7)/ms(n7)
329 vz7=v(3,n7)+dt12*a(3,n7)/ms(n7)
330 n8=ixs(9,n)
331 vx8=v(1,n8)+dt12*a(1,n8)/ms(n8)
332 vy8=v(2,n8)+dt12*a(2,n8)/ms(n8)
333 vz8=v(3,n8)+dt12*a(3,n8)/ms(n8)
334C
335 nsphdir=nint((sol2sph(2,n)-sol2sph(1,n))**third)
336C-----------------------------------------------
337 DO kp=1,sol2sph(2,n)-sol2sph(1,n)
338 np =sol2sph(1,n)+kp
339 ir=irst(1,np-first_sphsol+1)
340 is=irst(2,np-first_sphsol+1)
341 it=irst(3,np-first_sphsol+1)
342 ksi = a_gauss(ir,nsphdir)
343 eta = a_gauss(is,nsphdir)
344 zeta = a_gauss(it,nsphdir)
345C
346 phi1=(one-ksi)*(one-eta)*(one-zeta)
347 phi2=(one-ksi)*(one-eta)*(one+zeta)
348 phi3=(one+ksi)*(one-eta)*(one+zeta)
349 phi4=(one+ksi)*(one-eta)*(one-zeta)
350 phi5=(one-ksi)*(one+eta)*(one-zeta)
351 phi6=(one-ksi)*(one+eta)*(one+zeta)
352 phi7=(one+ksi)*(one+eta)*(one+zeta)
353 phi8=(one+ksi)*(one+eta)*(one-zeta)
354 vxi=one_over_8*(phi1*vx1+phi2*vx2+phi3*vx3+phi4*vx4+
355 . phi5*vx5+phi6*vx6+phi7*vx7+phi8*vx8)
356 vyi=one_over_8*(phi1*vy1+phi2*vy2+phi3*vy3+phi4*vy4+
357 . phi5*vy5+phi6*vy6+phi7*vy7+phi8*vy8)
358 vzi=one_over_8*(phi1*vz1+phi2*vz2+phi3*vz3+phi4*vz4+
359 . phi5*vz5+phi6*vz6+phi7*vz7+phi8*vz8)
360 inod=kxsp(3,np)
361 a(1,inod)=ms(inod)*(vxi-v(1,inod))*usdt
362 a(2,inod)=ms(inod)*(vyi-v(2,inod))*usdt
363 a(3,inod)=ms(inod)*(vzi-v(3,inod))*usdt
364 ENDDO
365 ENDDO
366C-----
367 ENDIF
368C-----
369 ENDIF
370 IF (iddw>0) CALL stoptimeg(ng)
371 END DO
372C--------
373 250 CONTINUE
374 END DO
375!$OMP END DO
376 END IF
377
378C-----------------------------------------------
379C Communication H and A on remote cells
380C-----------------------------------------------
381 IF(nspmd>1)THEN
382 IF(itask==0) THEN
383 ALLOCATE(sph_work%ASPHR(4,nsphr))
384 CALL spmd_sphgetf(kxsp,spbuf,a,ms,sph_work%ASPHR)
385c IF(NSPHIO/=0.AND.NSPHACT/=0)THEN
386c CALL SPMD_SPHGETV(KXSP,SPBUF,V,MS)
387c ENDIF
388 END IF
389 CALL my_barrier()
390 END IF
391C-----------------------------------------------
392C compute accelerations of ghost particles (temporary storage).
393C-----------------------------------------------
394 DO nc=1,nspcond
395 is=ispcond(3,nc)
396 ic=ispcond(2,nc)
397 islide=ispcond(5,nc)
398 ox=xframe(10,is)
399 oy=xframe(11,is)
400 oz=xframe(12,is)
401 nx=xframe(3*(ic-1)+1,is)
402 ny=xframe(3*(ic-1)+2,is)
403 nz=xframe(3*(ic-1)+3,is)
404 DO ns =itask+1,nsphact,nthread
405 n=waspact(ns)
406 js=ispsym(nc,n)
407 IF(js>0)THEN
408 inod=kxsp(3,n)
409 xi =x(1,inod)
410 yi =x(2,inod)
411 zi =x(3,inod)
412 axi=a(1,inod)
413 ayi=a(2,inod)
414 azi=a(3,inod)
415 IF(islide==0)THEN
416 axs=-axi
417 ays=-ayi
418 azs=-azi
419 ELSE
420 an=axi*nx+ayi*ny+azi*nz
421 axs=axi-two*an*nx
422 ays=ayi-two*an*ny
423 azs=azi-two*an*nz
424 ENDIF
425 waspsym(1,js)=axs
426 waspsym(2,js)=ays
427 waspsym(3,js)=azs
428 ENDIF
429 ENDDO
430C
431C Symmetrical particles of remote particles
432C
433 DO ns = itask+1,nsphr,nthread
434 js=ispsymr(nc,ns)
435 IF(js>0)THEN
436 xi =xsphr(3,ns)
437 yi =xsphr(4,ns)
438 zi =xsphr(5,ns)
439 axi=sph_work%ASPHR(1,ns)
440 ayi=sph_work%ASPHR(2,ns)
441 azi=sph_work%ASPHR(3,ns)
442 IF(islide==0)THEN
443 axs=-axi
444 ays=-ayi
445 azs=-azi
446 ELSE
447 an=axi*nx+ayi*ny+azi*nz
448 axs=axi-two*an*nx
449 ays=ayi-two*an*ny
450 azs=azi-two*an*nz
451 ENDIF
452 waspsym(1,js)=axs
453 waspsym(2,js)=ays
454 waspsym(3,js)=azs
455 ENDIF
456 END DO
457 ENDDO
458C
459C /---------------/
460 CALL my_barrier
461C /---------------/
462C-------------------------------------------
463 ehourt=zero
464 dtinv=one/dt12
465 DO ns =itask+1,nsphact,nthread
466 n=waspact(ns)
467 inod =kxsp(3,n)
468 unm=one/max(em30,ms(inod))
469 vxi=v(1,inod)+dt12*a(1,inod)*unm
470 vyi=v(2,inod)+dt12*a(2,inod)*unm
471 vzi=v(3,inod)+dt12*a(3,inod)*unm
472 vv=vxi*vxi+vyi*vyi+vzi*vzi
473 kv=half*ms(inod)*vv
474 ehourt=ehourt+kv
475 iprt=ipartsp(n)
476 partsav(8,iprt)=partsav(8,iprt)+kv
477 ENDDO
478C-----------------------------------------------
479C Conservative smoothing of velocities.
480C-----------------------------------------------
481C
482 DO ns =itask+1,nsphact,nthread
483 n =waspact(ns)
484 myadrn=3*(n-1)
485 wa(myadrn+1)=zero
486 wa(myadrn+2)=zero
487 wa(myadrn+3)=zero
488 ENDDO
489C
490 DO ns =itask+1,nsphact,nthread
491 n=waspact(ns)
492 myadrn=3*(n-1)
493 iprt =ipartsp(n)
494 imat =ipart(1,iprt)
495 iprop=ipart(2,iprt)
496 alpci=get_u_geo(4,iprop)
497C------
498 inod =kxsp(3,n)
499 nvois=kxsp(4,n)
500 xi=x(1,inod)
501 yi=x(2,inod)
502 zi=x(3,inod)
503 di=spbuf(1,n)
504 vxi=v(1,inod)
505 vyi=v(2,inod)
506 vzi=v(3,inod)
507 axi=a(1,inod)
508 ayi=a(2,inod)
509 azi=a(3,inod)
510 rhoi=spbuf(2,n)
511C------
512 alphai=wacomp(1,n)
513 betaxi=wacomp(2,n)
514 betayi=wacomp(3,n)
515 betazi=wacomp(4,n)
516C-----
517 DO j=1,nvois
518 jnod=ixsp(j,n)
519 IF(jnod>0)THEN
520 m=nod2sp(jnod)
521C
522C Solids to SPH, no interaction if both particles are inactive
523 IF(kxsp(2,n)<0.AND.kxsp(2,m)<0)cycle
524 xj=x(1,jnod)
525 yj=x(2,jnod)
526 zj=x(3,jnod)
527 rhoj=spbuf(2,m)
528 vxj=v(1,jnod)
529 vyj=v(2,jnod)
530 vzj=v(3,jnod)
531 axj=a(1,jnod)
532 ayj=a(2,jnod)
533 azj=a(3,jnod)
534 dj=spbuf(1,m)
535 dij=half*(di+dj)
536 CALL weight0(xi,yi,zi,xj,yj,zj,dij,wght)
537 betai=one+betaxi*(xi-xj)+betayi*(yi-yj)+betazi*(zi-zj)
538 alphaj=wacomp(1,m)
539 betaxj=wacomp(2,m)
540 betayj=wacomp(3,m)
541 betazj=wacomp(4,m)
542 betaj=one+betaxj*(xj-xi)+betayj*(yj-yi)+betazj*(zj-zi)
543 wght=wght*(alphai*betai+alphaj*betaj)*half
544 fact=two*wght/(rhoi+rhoj)
545 wax=axj-axi+ms(inod)*(vxj-vxi)*dtinv
546 way=ayj-ayi+ms(inod)*(vyj-vyi)*dtinv
547 waz=azj-azi+ms(inod)*(vzj-vzi)*dtinv
548 faci= alpci*spbuf(12,m)*fact
549 ELSE ! cellule remote
550 nn = -jnod
551C
552C Solids to SPH, no interaction if both particles are inactive
553 IF(kxsp(2,n)<=0.AND.xsphr(13,nn)<=0)cycle
554 xj=xsphr(3,nn)
555 yj=xsphr(4,nn)
556 zj=xsphr(5,nn)
557 rhoj=xsphr(7,nn)
558 vxj=xsphr(9,nn)
559 vyj=xsphr(10,nn)
560 vzj=xsphr(11,nn)
561 axj=sph_work%ASPHR(1,nn)
562 ayj=sph_work%ASPHR(2,nn)
563 azj=sph_work%ASPHR(3,nn)
564 dj=xsphr(2,nn)
565 dij=half*(di+dj)
566 CALL weight0(xi,yi,zi,xj,yj,zj,dij,wght)
567 betai=one+betaxi*(xi-xj)+betayi*(yi-yj)+betazi*(zi-zj)
568 alphaj=wacompr(1,nn)
569 betaxj=wacompr(2,nn)
570 betayj=wacompr(3,nn)
571 betazj=wacompr(4,nn)
572 betaj=one+betaxj*(xj-xi)+betayj*(yj-yi)+betazj*(zj-zi)
573 wght=wght*(alphai*betai+alphaj*betaj)*half
574 fact=two*wght/(rhoi+rhoj)
575 wax=axj-axi+ms(inod)*(vxj-vxi)*dtinv
576 way=ayj-ayi+ms(inod)*(vyj-vyi)*dtinv
577 waz=azj-azi+ms(inod)*(vzj-vzi)*dtinv
578 faci= alpci*xsphr(8,nn)*fact
579 END IF
580 wa(myadrn+1)=wa(myadrn+1)+faci*wax
581 wa(myadrn+2)=wa(myadrn+2)+faci*way
582 wa(myadrn+3)=wa(myadrn+3)+faci*waz
583 END DO
584C-----
585C symmetrical part.
586 nvoiss=kxsp(6,n)
587 DO j=kxsp(5,n)+1,kxsp(5,n)+nvoiss
588 js=ixsp(j,n)
589 IF(js>0)THEN
590 sm=js/(nspcond+1)
591C
592C Solids to SPH, no interaction if both particles are inactive
593 IF(kxsp(2,n)<=0.AND.kxsp(2,sm)<=0)cycle
594 nc=mod(js,nspcond+1)
595 js=ispsym(nc,sm)
596 xj=xspsym(1,js)
597 yj=xspsym(2,js)
598 zj=xspsym(3,js)
599 vxj=vspsym(1,js)
600 vyj=vspsym(2,js)
601 vzj=vspsym(3,js)
602 axj=waspsym(1,js)
603 ayj=waspsym(2,js)
604 azj=waspsym(3,js)
605 rhoj=spbuf(2,sm)
606 dj =spbuf(1,sm)
607 dij =half*(di+dj)
608 CALL weight0(xi,yi,zi,xj,yj,zj,dij,wght)
609 jnod=kxsp(3,sm)
610 betai=one +betaxi*(xi-xj)+betayi*(yi-yj)+betazi*(zi-zj)
611 alphaj=wacomp(1,sm)
612C BETAXJ=WACOMP(2,SM)
613C BETAYJ=WACOMP(3,SM)
614C BETAZJ=WACOMP(4,SM)
615 betaxj=wsmcomp(1,js)
616 betayj=wsmcomp(2,js)
617 betazj=wsmcomp(3,js)
618 betaj=one+betaxj*(xj-xi)+betayj*(yj-yi)+betazj*(zj-zi)
619 wght=wght*(alphai*betai+alphaj*betaj)*half
620 fact=alpci*two*spbuf(12,sm)*wght/(rhoi+rhoj)
621 wax=axj-axi+ms(inod)*(vxj-vxi)*dtinv
622 way=ayj-ayi+ms(inod)*(vyj-vyi)*dtinv
623 waz=azj-azi+ms(inod)*(vzj-vzi)*dtinv
624 ELSE ! Symmetrical particle particle Remote
625 sm=-js/(nspcond+1)
626C
627C Solids to SPH, no interaction if both particles are inactive
628 IF(kxsp(2,n)<=0.AND.xsphr(13,sm)<=0)cycle
629 nc=mod(-js,nspcond+1)
630 js=ispsymr(nc,sm)
631 xj=xspsym(1,js)
632 yj=xspsym(2,js)
633 zj=xspsym(3,js)
634 vxj=vspsym(1,js)
635 vyj=vspsym(2,js)
636 vzj=vspsym(3,js)
637 axj=waspsym(1,js)
638 ayj=waspsym(2,js)
639 azj=waspsym(3,js)
640 rhoj=xsphr(7,sm)
641 dj =xsphr(2,sm)
642 dij =half*(di+dj)
643 CALL weight0(xi,yi,zi,xj,yj,zj,dij,wght)
644 betai=one +betaxi*(xi-xj)+betayi*(yi-yj)+betazi*(zi-zj)
645 alphaj=wacompr(1,sm)
646C BETAXJ=WACOMPR(2,SM)
647C BETAYJ=WACOMPR(3,SM)
648C BETAZJ=WACOMPR(4,SM)
649 betaxj=wsmcomp(1,js)
650 betayj=wsmcomp(2,js)
651 betazj=wsmcomp(3,js)
652 betaj=one+betaxj*(xj-xi)+betayj*(yj-yi)+betazj*(zj-zi)
653 wght=wght*(alphai*betai+alphaj*betaj)*half
654 fact=alpci*two*xsphr(8,sm)*wght/(rhoi+rhoj)
655 wax=axj-axi+ms(inod)*(vxj-vxi)*dtinv
656 way=ayj-ayi+ms(inod)*(vyj-vyi)*dtinv
657 waz=azj-azi+ms(inod)*(vzj-vzi)*dtinv
658 END IF
659 wa(myadrn+1)=wa(myadrn+1)+fact*wax
660 wa(myadrn+2)=wa(myadrn+2)+fact*way
661 wa(myadrn+3)=wa(myadrn+3)+fact*waz
662 END DO
663 END DO
664C-----
665
666C barrier on A & ASPHR
667C /---------------/
668 CALL my_barrier
669C /---------------/
670
671C Force assembly in A
672 IF(nsphsol==0)THEN
673 DO ns=itask+1,nsphact,nthread
674 n=waspact(ns)
675 myadrn=3*(n-1)
676 inod=kxsp(3,n)
677 a(1,inod)=a(1,inod)+wa(myadrn+1)
678 a(2,inod)=a(2,inod)+wa(myadrn+2)
679 a(3,inod)=a(3,inod)+wa(myadrn+3)
680 END DO
681 ELSEIF(itask==0)THEN
682C
683C a paralleliser smp
684 ii=0
685 DO ns=1,nsphact
686 n=waspact(ns)
687 myadrn=3*(n-1)
688 IF(sph2sol(n)==0)THEN
689 inod=kxsp(3,n)
690 a(1,inod)=a(1,inod)+wa(myadrn+1)
691 a(2,inod)=a(2,inod)+wa(myadrn+2)
692 a(3,inod)=a(3,inod)+wa(myadrn+3)
693C
694 ELSEIF (sol2sph_typ(sph2sol(n))==4) THEN
695C---------------
696C------ TETRA --
697C---------------
698C for computing Ehour only
699 inod=kxsp(3,n)
700 a(1,inod)=a(1,inod)+wa(myadrn+1)
701 a(2,inod)=a(2,inod)+wa(myadrn+2)
702 a(3,inod)=a(3,inod)+wa(myadrn+3)
703C----------
704 nsol=sph2sol(n)
705C
706 n1=ixs(2,nsol)
707 n2=ixs(4,nsol)
708 n3=ixs(7,nsol)
709 n4=ixs(6,nsol)
710C
711 ir=irst(1,n-first_sphsol+1)
712 is=irst(2,n-first_sphsol+1)
713 it=irst(3,n-first_sphsol+1)
714 nsphdir=igeo(37,ixs(10,nsol))
715C
716 ksi = a_gauss(ir,nsphdir)
717 eta = a_gauss(is,nsphdir)
718 zeta = a_gauss(it,nsphdir)
719C
720 phi1=one_over_8*(one-ksi)*(one-eta)*(one-zeta)
721 phi2=one_over_8*(one-ksi)*(one-eta)*(one+zeta)
722 phi3=one_over_8*(one+ksi)*(one-eta)*(one+zeta)
723 phi4=one_over_8*(one+ksi)*(one-eta)*(one-zeta)
724C
725 ii=ii+1
726 sph_work%AS(1,ii)=phi1*wa(myadrn+1)
727 sph_work%AS(2,ii)=phi1*wa(myadrn+2)
728 sph_work%AS(3,ii)=phi1*wa(myadrn+3)
729
730 ii=ii+1
731 sph_work%AS(1,ii)=phi2*wa(myadrn+1)
732 sph_work%AS(2,ii)=phi2*wa(myadrn+2)
733 sph_work%AS(3,ii)=phi2*wa(myadrn+3)
734
735 ii=ii+1
736 sph_work%AS(1,ii)=phi3*wa(myadrn+1)
737 sph_work%AS(2,ii)=phi3*wa(myadrn+2)
738 sph_work%AS(3,ii)=phi3*wa(myadrn+3)
739
740 ii=ii+1
741 sph_work%AS(1,ii)=phi4*wa(myadrn+1)
742 sph_work%AS(2,ii)=phi4*wa(myadrn+2)
743 sph_work%AS(3,ii)=phi4*wa(myadrn+3)
744C
745 ELSE
746C---------------
747C------ HEXA- --
748C---------------
749C for computing Ehour only
750 inod=kxsp(3,n)
751 a(1,inod)=a(1,inod)+wa(myadrn+1)
752 a(2,inod)=a(2,inod)+wa(myadrn+2)
753 a(3,inod)=a(3,inod)+wa(myadrn+3)
754C----------
755 nsol=sph2sol(n)
756C
757 n1=ixs(2,nsol)
758 n2=ixs(3,nsol)
759 n3=ixs(4,nsol)
760 n4=ixs(5,nsol)
761 n5=ixs(6,nsol)
762 n6=ixs(7,nsol)
763 n7=ixs(8,nsol)
764 n8=ixs(9,nsol)
765C
766 ir=irst(1,n-first_sphsol+1)
767 is=irst(2,n-first_sphsol+1)
768 it=irst(3,n-first_sphsol+1)
769 nsphdir=nint((sol2sph(2,nsol)-sol2sph(1,nsol))**third)
770C
771 ksi = a_gauss(ir,nsphdir)
772 eta = a_gauss(is,nsphdir)
773 zeta = a_gauss(it,nsphdir)
774C
775 phi1=one_over_8*(one-ksi)*(one-eta)*(one-zeta)
776 phi2=one_over_8*(one-ksi)*(one-eta)*(one+zeta)
777 phi3=one_over_8*(one+ksi)*(one-eta)*(one+zeta)
778 phi4=one_over_8*(one+ksi)*(one-eta)*(one-zeta)
779 phi5=one_over_8*(one-ksi)*(one+eta)*(one-zeta)
780 phi6=one_over_8*(one-ksi)*(one+eta)*(one+zeta)
781 phi7=one_over_8*(one+ksi)*(one+eta)*(one+zeta)
782 phi8=one_over_8*(one+ksi)*(one+eta)*(one-zeta)
783C
784 ii=ii+1
785 sph_work%AS(1,ii)=phi1*wa(myadrn+1)
786 sph_work%AS(2,ii)=phi1*wa(myadrn+2)
787 sph_work%AS(3,ii)=phi1*wa(myadrn+3)
788
789 ii=ii+1
790 sph_work%AS(1,ii)=phi2*wa(myadrn+1)
791 sph_work%AS(2,ii)=phi2*wa(myadrn+2)
792 sph_work%AS(3,ii)=phi2*wa(myadrn+3)
793
794 ii=ii+1
795 sph_work%AS(1,ii)=phi3*wa(myadrn+1)
796 sph_work%AS(2,ii)=phi3*wa(myadrn+2)
797 sph_work%AS(3,ii)=phi3*wa(myadrn+3)
798
799 ii=ii+1
800 sph_work%AS(1,ii)=phi4*wa(myadrn+1)
801 sph_work%AS(2,ii)=phi4*wa(myadrn+2)
802 sph_work%AS(3,ii)=phi4*wa(myadrn+3)
803
804 ii=ii+1
805 sph_work%AS(1,ii)=phi5*wa(myadrn+1)
806 sph_work%AS(2,ii)=phi5*wa(myadrn+2)
807 sph_work%AS(3,ii)=phi5*wa(myadrn+3)
808
809 ii=ii+1
810 sph_work%AS(1,ii)=phi6*wa(myadrn+1)
811 sph_work%AS(2,ii)=phi6*wa(myadrn+2)
812 sph_work%AS(3,ii)=phi6*wa(myadrn+3)
813
814 ii=ii+1
815 sph_work%AS(1,ii)=phi7*wa(myadrn+1)
816 sph_work%AS(2,ii)=phi7*wa(myadrn+2)
817 sph_work%AS(3,ii)=phi7*wa(myadrn+3)
818
819 ii=ii+1
820 sph_work%AS(1,ii)=phi8*wa(myadrn+1)
821 sph_work%AS(2,ii)=phi8*wa(myadrn+2)
822 sph_work%AS(3,ii)=phi8*wa(myadrn+3)
823C
824 END IF
825 ENDDO
826C------
827 CALL foat_to_6_float(1,3*ii,sph_work%AS,sph_work%AS6)
828C------
829 ii=0
830 DO ns=1,nsphact
831 n=waspact(ns)
832 IF(sph2sol(n)/=0)THEN
833 IF (sol2sph_typ(sph2sol(n))==4) THEN
834C---------------
835C------ TETRA --
836C---------------
837C
838 nsol=sph2sol(n)
839C
840 n1=ixs(2,nsol)
841 n2=ixs(3,nsol)
842 n3=ixs(4,nsol)
843 n4=ixs(5,nsol)
844C
845 sph_work%ITAG(n1)=1
846 sph_work%ITAG(n2)=1
847 sph_work%ITAG(n3)=1
848 sph_work%ITAG(n4)=1
849C
850 ii1=ii
851 DO j=1,6
852 ii=ii1+1
853 sph_work%A6(j,1,n1)=sph_work%A6(j,1,n1)+sph_work%AS6(j,1,ii)
854 sph_work%A6(j,2,n1)=sph_work%A6(j,2,n1)+sph_work%AS6(j,2,ii)
855 sph_work%A6(j,3,n1)=sph_work%A6(j,3,n1)+sph_work%AS6(j,3,ii)
856C
857 ii=ii+1
858 sph_work%A6(j,1,n2)=sph_work%A6(j,1,n2)+sph_work%AS6(j,1,ii)
859 sph_work%A6(j,2,n2)=sph_work%A6(j,2,n2)+sph_work%AS6(j,2,ii)
860 sph_work%A6(j,3,n2)=sph_work%A6(j,3,n2)+sph_work%AS6(j,3,ii)
861C
862 ii=ii+1
863 sph_work%A6(j,1,n3)=sph_work%A6(j,1,n3)+sph_work%AS6(j,1,ii)
864 sph_work%A6(j,2,n3)=sph_work%A6(j,2,n3)+sph_work%AS6(j,2,ii)
865 sph_work%A6(j,3,n3)=sph_work%A6(j,3,n3)+sph_work%AS6(j,3,ii)
866C
867 ii=ii+1
868 sph_work%A6(j,1,n4)=sph_work%A6(j,1,n4)+sph_work%AS6(j,1,ii)
869 sph_work%A6(j,2,n4)=sph_work%A6(j,2,n4)+sph_work%AS6(j,2,ii)
870 sph_work%A6(j,3,n4)=sph_work%A6(j,3,n4)+sph_work%AS6(j,3,ii)
871C
872 END DO
873C
874 ELSE
875C---------------
876C------ HEXA ---
877C---------------
878C
879 nsol=sph2sol(n)
880C
881 n1=ixs(2,nsol)
882 n2=ixs(3,nsol)
883 n3=ixs(4,nsol)
884 n4=ixs(5,nsol)
885 n5=ixs(6,nsol)
886 n6=ixs(7,nsol)
887 n7=ixs(8,nsol)
888 n8=ixs(9,nsol)
889C
890 sph_work%ITAG(n1)=1
891 sph_work%ITAG(n2)=1
892 sph_work%ITAG(n3)=1
893 sph_work%ITAG(n4)=1
894 sph_work%ITAG(n5)=1
895 sph_work%ITAG(n6)=1
896 sph_work%ITAG(n7)=1
897 sph_work%ITAG(n8)=1
898C
899 ii1=ii
900 DO j=1,6
901 ii=ii1+1
902 sph_work%A6(j,1,n1)=sph_work%A6(j,1,n1)+sph_work%AS6(j,1,ii)
903 sph_work%A6(j,2,n1)=sph_work%A6(j,2,n1)+sph_work%AS6(j,2,ii)
904 sph_work%A6(j,3,n1)=sph_work%A6(j,3,n1)+sph_work%AS6(j,3,ii)
905C
906 ii=ii+1
907 sph_work%A6(j,1,n2)=sph_work%A6(j,1,n2)+sph_work%AS6(j,1,ii)
908 sph_work%A6(j,2,n2)=sph_work%A6(j,2,n2)+sph_work%AS6(j,2,ii)
909 sph_work%A6(j,3,n2)=sph_work%A6(j,3,n2)+sph_work%AS6(j,3,ii)
910C
911 ii=ii+1
912 sph_work%A6(j,1,n3)=sph_work%A6(j,1,n3)+sph_work%AS6(j,1,ii)
913 sph_work%A6(j,2,n3)=sph_work%A6(j,2,n3)+sph_work%AS6(j,2,ii)
914 sph_work%A6(j,3,n3)=sph_work%A6(j,3,n3)+sph_work%AS6(j,3,ii)
915C
916 ii=ii+1
917 sph_work%A6(j,1,n4)=sph_work%A6(j,1,n4)+sph_work%AS6(j,1,ii)
918 sph_work%A6(j,2,n4)=sph_work%A6(j,2,n4)+sph_work%AS6(j,2,ii)
919 sph_work%A6(j,3,n4)=sph_work%A6(j,3,n4)+sph_work%AS6(j,3,ii)
920C
921 ii=ii+1
922 sph_work%A6(j,1,n5)=sph_work%A6(j,1,n5)+sph_work%AS6(j,1,ii)
923 sph_work%A6(j,2,n5)=sph_work%A6(j,2,n5)+sph_work%AS6(j,2,ii)
924 sph_work%A6(j,3,n5)=sph_work%A6(j,3,n5)+sph_work%AS6(j,3,ii)
925C
926 ii=ii+1
927 sph_work%A6(j,1,n6)=sph_work%A6(j,1,n6)+sph_work%AS6(j,1,ii)
928 sph_work%A6(j,2,n6)=sph_work%A6(j,2,n6)+sph_work%AS6(j,2,ii)
929 sph_work%A6(j,3,n6)=sph_work%A6(j,3,n6)+sph_work%AS6(j,3,ii)
930C
931 ii=ii+1
932 sph_work%A6(j,1,n7)=sph_work%A6(j,1,n7)+sph_work%AS6(j,1,ii)
933 sph_work%A6(j,2,n7)=sph_work%A6(j,2,n7)+sph_work%AS6(j,2,ii)
934 sph_work%A6(j,3,n7)=sph_work%A6(j,3,n7)+sph_work%AS6(j,3,ii)
935C
936 ii=ii+1
937 sph_work%A6(j,1,n8)=sph_work%A6(j,1,n8)+sph_work%AS6(j,1,ii)
938 sph_work%A6(j,2,n8)=sph_work%A6(j,2,n8)+sph_work%AS6(j,2,ii)
939 sph_work%A6(j,3,n8)=sph_work%A6(j,3,n8)+sph_work%AS6(j,3,ii)
940C
941 END DO
942C
943 ENDIF
944C
945 END IF
946 END DO
947 END IF
948C
949 ENDIF !(IF (NUMSPH > 0)
950C-----
951 IF ((sol2sph_flag > 0).AND.(itask==0)) THEN
952 IF(nspmd > 1)THEN
953 sz=19
954 lenr = iad_elem(1,nspmd+1)-iad_elem(1,1)
955 CALL spmd_exch_a_sol2sph(
956 1 sph_work%A6 ,sph_work%ITAG ,iad_elem ,fr_elem,sz,
957 2 lenr )
958 END IF
959C-----
960 DO n=1,numnod
961 IF(sph_work%ITAG(n)/=0)THEN
962 a(1,n)=a(1,n)+sph_work%A6(1,1,n)+sph_work%A6(2,1,n)+sph_work%A6(3,1,n)
963 . +sph_work%A6(4,1,n)+sph_work%A6(5,1,n)+sph_work%A6(6,1,n)
964 a(2,n)=a(2,n)+sph_work%A6(1,2,n)+sph_work%A6(2,2,n)+sph_work%A6(3,2,n)
965 . +sph_work%A6(4,2,n)+sph_work%A6(5,2,n)+sph_work%A6(6,2,n)
966 a(3,n)=a(3,n)+sph_work%A6(1,3,n)+sph_work%A6(2,3,n)+sph_work%A6(3,3,n)
967 . +sph_work%A6(4,3,n)+sph_work%A6(5,3,n)+sph_work%A6(6,3,n)
968 END IF
969 END DO
970 ENDIF
971
972C-------------------------------------------
973C /---------------/
974 CALL my_barrier
975C /---------------/
976C
977 IF (numsph > 0) THEN
978 DO ns =itask+1,nsphact,nthread
979 n =waspact(ns)
980 inod =kxsp(3,n)
981 unm=one/max(em30,ms(inod))
982 vxi=v(1,inod)+dt12*a(1,inod)*unm
983 vyi=v(2,inod)+dt12*a(2,inod)*unm
984 vzi=v(3,inod)+dt12*a(3,inod)*unm
985 vv=vxi*vxi+vyi*vyi+vzi*vzi
986 kv=half*ms(inod)*vv
987 ehourt=ehourt-kv
988 iprt=ipartsp(n)
989 partsav(8,iprt)=partsav(8,iprt)-kv
990 ENDDO
991!$OMP ATOMIC
992 ehour=ehour+ehourt
993 ENDIF
994C
995 IF(nspmd>1 .AND. itask==0.AND.ALLOCATED(sph_work%ASPHR)) DEALLOCATE(sph_work%ASPHR)
996C-------------------------------------------
997 RETURN
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
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1257
sphbox::ispsymr
integer, dimension(:,:), allocatable ispsymr
Definition sphbox.F:93
sphbox::nsphr
integer nsphr
Definition sphbox.F:83
foat_to_6_float
subroutine foat_to_6_float(jft, jlt, f, f6)
Definition parit.F:226
spmd_exch_a_sol2sph
subroutine spmd_exch_a_sol2sph(a6, itag, iad_elem, fr_elem, size, lenr)
Definition spmd_sph.F:35
spmd_sphgetf
subroutine spmd_sphgetf(kxsp, spbuf, a, ms, asphr)
Definition spmd_sph.F:935
my_barrier
subroutine my_barrier
Definition machine.F:31
weight0
subroutine weight0(xi, yi, zi, xj, yj, zj, h, w)
Definition weight.F:34