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 42 of file splissv.F.

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