OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
forintp.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!|| forintp ../engine/source/elements/forintp.F
25!||--- called by ------------------------------------------------------
26!|| resol ../engine/source/engine/resol.F
27!||--- calls -----------------------------------------------------
28!|| initbuf ../engine/share/resol/initbuf.F
29!|| my_barrier ../engine/source/system/machine.F
30!|| soltosphp ../engine/source/elements/sph/soltosph.F
31!|| spadah ../engine/source/elements/sph/spadah.F
32!|| spdens ../engine/source/elements/sph/spdens.F
33!|| spforcp ../engine/source/elements/sph/spforcp.F
34!|| spgauge ../engine/source/elements/sph/spgauge.F
35!|| spgradt ../engine/source/elements/sph/sptemp.F
36!|| spgtsym ../engine/source/elements/sph/sptemp.F
37!|| splaplt ../engine/source/elements/sph/sptemp.F
38!|| spmd_sphgetg ../engine/source/mpi/elements/spmd_sph.F
39!|| spmd_sphgetstb ../engine/source/mpi/elements/spmd_sph.F
40!|| spmd_sphgett ../engine/source/mpi/elements/spmd_sph.F
41!|| spmd_sphgetw ../engine/source/mpi/elements/spmd_sph.F
42!|| spmd_sphgetwa ../engine/source/mpi/elements/spmd_sph.f
43!|| sponfprs ../engine/source/elements/sph/sponfprs.F
44!|| sponfro ../engine/source/elements/sph/sponfro.F
45!|| spscomp ../engine/source/elements/sph/spcompl.F
46!|| spsgsym ../engine/source/elements/sph/spsgsym.F
47!|| spstabs ../engine/source/elements/sph/spstab.f
48!|| spstabw ../engine/source/elements/sph/spstab.F
49!|| spstres ../engine/source/elements/sph/spstres.F
50!|| startime ../engine/source/system/timer_mod.F90
51!|| startimeg ../engine/source/system/timer.F
52!|| stoptime ../engine/source/system/timer_mod.F90
53!|| stoptimeg ../engine/source/system/timer.F
54!||--- uses -----------------------------------------------------
55!|| dt_mod ../engine/source/modules/dt_mod.F
56!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
57!|| element_mod ../common_source/modules/elements/element_mod.F90
58!|| glob_therm_mod ../common_source/modules/mat_elem/glob_therm_mod.F90
59!|| initbuf_mod ../engine/share/resol/initbuf.F
60!|| mat_elem_mod ../common_source/modules/mat_elem/mat_elem_mod.F90
61!|| nlocal_reg_mod ../common_source/modules/nlocal_reg_mod.F
62!|| output_mod ../common_source/modules/output/output_mod.F90
63!|| sensor_mod ../common_source/modules/sensor_mod.F90
64!|| sph_work_mod ../common_source/modules/mat_elem/sph_work.F90
65!|| sphbox ../engine/share/modules/sphbox.F
66!|| table_mod ../engine/share/modules/table_mod.F
67!|| timer_mod ../engine/source/system/timer_mod.F90
68!||====================================================================
69 SUBROUTINE forintp(TIMERS,
70 1 PM ,GEO ,X ,A ,V ,
71 2 MS ,W ,ELBUF_TAB ,WA ,FV ,
72 3 STIFN ,PLD ,BUFMAT ,PARTSAV ,NLOC_DMG ,
73 4 FSAV ,DT2T ,IADS ,IPARG ,NPC ,
74 5 NELTST ,ITYPTST ,IPART ,ITAB ,ISKY ,
75 6 BUFGEO ,FSKYI ,XFRAME ,KXSP ,IXSP ,
76 7 NOD2SP ,IPARTSP ,SPBUF ,ISPCOND ,ISPSYM ,
77 8 XSPSYM ,VSPSYM ,
78 9 WASPH ,LPRTSPH ,LONFSPH ,WASPACT ,ISPHIO ,
79 A VSPHIO ,SPHVELN ,ITASK ,IPM ,GRESAV ,
80 B GRTH ,IGRTH ,TABLE ,LGAUGE ,GAUGE ,
81 C NGROUNC ,IGROUNC ,IXS ,IRST ,SOL2SPH ,
82 D SPH2SOL ,FSKYV ,FSKY ,IGEO ,TEMP ,
83 E FTHE ,FTHESKYI ,SPHG_F6 ,WSMCOMP ,SOL2SPH_TYP,
84 F MAT_ELEM ,OUTPUT ,SPH_IORD1 ,SNPC ,STF,
85 G SBUFMAT ,IDTMINS ,NSVOIS ,IRESP ,MAXFUNC ,
86 . IMON_MAT ,USERL_AVAIL,impl_s ,idyna ,
87 . DT ,GLOB_THERM,SPH_WORK ,WFEXT ,sensors )
88C-----------------------------------------------
89C M o d u l e s
90C-----------------------------------------------
91 USE timer_mod
92 USE initbuf_mod
93 USE sphbox
94 USE table_mod
95 USE mat_elem_mod
96 USE elbufdef_mod
97 USE nlocal_reg_mod
98 USE output_mod , ONLY : output_
99 USE dt_mod
100 use glob_therm_mod
101 use sensor_mod
102 USE sph_work_mod
103 use element_mod , only : nixs
104C----6---------------------------------------------------------------7---------8
105C I m p l i c i t T y p e s
106C-----------------------------------------------
107#include "implicit_f.inc"
108#include "comlock.inc"
109C-----------------------------------------------
110C G l o b a l P a r a m e t e r s
111C-----------------------------------------------
112#include "mvsiz_p.inc"
113C-----------------------------------------------
114C C o m m o n B l o c k s
115C-----------------------------------------------
116#include "com01_c.inc"
117#include "com04_c.inc"
118#include "com08_c.inc"
119#include "sphcom.inc"
120#include "param_c.inc"
121#include "parit_c.inc"
122#include "vect01_c.inc"
123#include "scr07_c.inc"
124#include "scr17_c.inc"
125#include "task_c.inc"
126#include "units_c.inc"
127#include "scr02_c.inc"
128#include "scr18_c.inc"
129C-----------------------------------------------------------------
130C D u m m y A r g u m e n t s
131C-----------------------------------------------
132 TYPE(timer_), INTENT(INOUT) :: TIMERS
133 INTEGER, INTENT(IN) :: SNPC
134 INTEGER, INTENT(IN) :: STF
135 INTEGER, INTENT(IN) :: SBUFMAT
136 INTEGER, INTENT(IN) :: NSVOIS
137 INTEGER, INTENT(IN) :: IDTMINS
138 INTEGER ,INTENT(IN) :: IRESP
139 INTEGER ,INTENT(IN) :: MAXFUNC
140 INTEGER, INTENT(IN) :: IMPL_S
141 INTEGER, INTENT(IN) :: IDYNA
142 INTEGER, INTENT(IN) :: USERL_AVAIL
143 INTEGER, INTENT(IN) :: IMON_MAT
144 INTEGER IPART(LIPART1,*) ,NPC(*), IPARG(NPARG,*),IADS(8,*),
145 . NELTST, ITYPTST, IPARTSP(*), ISKY(*), ITAB(*),IPM(*),
146 . KXSP(NISP,*),IXSP(KVOISPH,*),NOD2SP(*),
147 . ISPCOND(NISPCOND,*),ISPSYM(NSPCOND,*),
148 . IGEO(NPROPGI,*),
149 . LPRTSPH(2,0:NPART),LONFSPH(*),WASPACT(*),ISPHIO(NISPHIO,*),
150 . itask,grth(*),igrth(*), lgauge(3,*), ngrounc, igrounc(*),
151 . ixs(nixs,*), irst(3,*), sol2sph(2,*), sph2sol(*), sol2sph_typ(*)
152 INTEGER, INTENT(IN) :: SPH_IORD1
153 my_real
154 . x(3,*), v(3,*), ms(*), w(*), pm(npropm,*), geo(npropg,*),
155 . bufmat(*), bufgeo(*), pld(*) ,
156 . fsav(nthvki,*), wa(*), fv(*), a(3,*),
157 . partsav(*), stifn(*), fskyi(lskyi,4) ,
158 . xframe(nxframe,*), spbuf(nspbuf,*), xspsym(3,*), vspsym(3,*),
159 . dt2t, wasph(*), vsphio(*),
160 . sphveln(*),gresav(*), gauge(llgauge,*),
161 . fskyv(lsky,8),fsky(8,lsky),temp(*),fthe(*),ftheskyi(*),wsmcomp(*)
162 TYPE(ttable) TABLE(*)
163 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
164 TYPE (NLOCAL_STR_) , TARGET :: NLOC_DMG
165 DOUBLE PRECISION SPHG_F6(4,6,NBGAUGE)
166
167 TYPE(OUTPUT_), INTENT(INOUT) :: OUTPUT !< output structure
168 TYPE (MAT_ELEM_) ,INTENT(INOUT) :: MAT_ELEM
169 TYPE (DT_), INTENT(IN) :: DT
170 type(glob_therm_) ,intent(inout) :: glob_therm
171 TYPE (SPH_WORK_),INTENT(INOUT) :: SPH_WORK
172 type (sensors_) ,intent(in) :: sensors !< sensor structure
173 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
174C-----------------------------------------------
175C L o c a l V a r i a b l e s
176C-----------------------------------------------
177 INTEGER NDSVOID
178 INTEGER I, N, NG, NVC, K, ISTRA,
179 . KAD, NELEM, NEL, OFFSET, NSG,
180 . inod,mx,ns,kvnorm,myadrn,adrn, nisky_l,
181 . nsol, nski, n1, n2, n3, n4, n5, n6, n7, n8,
182 . k1, k2, k3, k4, k5, k6, k7, k8, ir, is, it, nsphdir, stat,
183 . iexpan,ibid
184 my_real
185 . dtx,dt05,rhoi,vi,
186 . phi1,phi2,phi3,phi4,phi5,phi6,phi7,phi8,
187 . ksi, eta, zeta,
188 . voln(mvsiz)
189 my_real, DIMENSION(MVSIZ,6) :: svis
190C
191 TYPE(g_bufel_) ,POINTER :: GBUF
192C-----------------------------------------------
193 my_real
194 . A_GAUSS(9,9),A_GAUSS_TETRA(9,9)
195 DATA a_gauss /
196 1 0. ,0. ,0. ,
197 1 0. ,0. ,0. ,
198 1 0. ,0. ,0. ,
199 2 -.5 ,0.5 ,0. ,
200 2 0. ,0. ,0. ,
201 2 0. ,0. ,0. ,
202 3 -.666666666666666,0. ,0.666666666666666,
203 3 0. ,0. ,0. ,
204 3 0. ,0. ,0. ,
205 4 -.75 ,-.25 ,0.25 ,
206 4 0.75 ,0. ,0. ,
207 4 0. ,0. ,0. ,
208 5 -.8 ,-.4 ,0. ,
209 5 0.4 ,0.8 ,0. ,
210 5 0. ,0. ,0. ,
211 6 -.833333333333333,-.5 ,-.166666666666666,
212 6 0.166666666666666,0.5 ,0.833333333333333,
213 6 0. ,0. ,0. ,
214 7 -.857142857142857,-.571428571428571,-.285714285714285,
215 7 0. ,0.285714285714285,0.571428571428571,
216 7 0.857142857142857,0. ,0. ,
217 8 -.875 ,-.625 ,-.375 ,
218 8 -.125 ,0.125 ,0.375,
219 8 0.625 ,0.875 ,0. ,
220 9 -.888888888888888,-.666666666666666,-.444444444444444,
221 9 -.222222222222222,0. ,0.222222222222222,
222 9 0.444444444444444,0.666666666666666,0.888888888888888/
223C-----------------------------------------------
224 DATA a_gauss_tetra /
225 1 0.250000000000000,0.000000000000000,0.000000000000000,
226 1 0.000000000000000,0.000000000000000,0.000000000000000,
227 1 0.000000000000000,0.000000000000000,0.000000000000000,
228 2 0.166666666666667,0.500000000000000,0.000000000000000,
229 2 0.000000000000000,0.000000000000000,0.000000000000000,
230 2 0.000000000000000,0.000000000000000,0.000000000000000,
231 3 0.125000000000000,0.375000000000000,0.625000000000000,
232 3 0.000000000000000,0.000000000000000,0.000000000000000,
233 3 0.000000000000000,0.000000000000000,0.000000000000000,
234 4 0.100000000000000,0.300000000000000,0.500000000000000,
235 4 0.700000000000000,0.000000000000000,0.000000000000000,
236 4 0.000000000000000,0.000000000000000,0.000000000000000,
237 5 0.083333333333333,0.250000000000000,0.416666666666667,
238 5 0.583333333333333,0.750000000000000,0.000000000000000,
239 5 0.000000000000000,0.000000000000000,0.000000000000000,
240 6 0.071428571428571,0.214285714285714,0.357142857142857,
241 6 0.500000000000000,0.642857142857143,0.785714285714286,
242 6 0.000000000000000,0.000000000000000,0.000000000000000,
243 7 0.062500000000000,0.187500000000000,0.312500000000000,
244 7 0.437500000000000,0.562500000000000,0.687500000000000,
245 7 0.812500000000000,0.000000000000000,0.000000000000000,
246 8 0.055555555555556,0.166666666666667,0.277777777777778,
247 8 0.388888888888889,0.500000000000000,0.611111111111111,
248 8 0.722222222222222,0.833333333333333,0.000000000000000,
249 9 0.050000000000000,0.150000000000000,0.250000000000000,
250 9 0.350000000000000,0.450000000000000,0.550000000000000,
251 9 0.650000000000000,0.750000000000000,0.850000000000000/
252C-----------------------------------------------
253C calculation of densities and strain rates
254C load wa (1: 14*numsph), which must not be crushed
255C before SPSTRES (for transfer to element buffers)
256C
257C Divv = Wa (13.1: Numsph) and Rotv = Wa (14.1: Numsph) must not be
258C overwritten before the calculation of forces (SPFORC)
259C
260C Divv = WA (13.1: NUMSPH) must not be crushed before spadah
261C (adaptation of the particle diameter)
262C-----------------------------------------------
263 ibid = 0
264C
265C SPMD allowance to do in memory shared
266C
267 IF(itask==0) THEN
268 ALLOCATE(sph_work%WASIGSM(6*nsphsym))
269 sph_work%WASIGSM = zero
270 ENDIF
271 IF(itask==0 .AND. nspmd > 1)THEN
272 ALLOCATE(sph_work%WAR(10,nsphr))
273 ALLOCATE(sph_work%WTR(nsphr))
274 ALLOCATE(sph_work%WGR(3,nsphr))
275 ALLOCATE(sph_work%LAMBDR(nsphr))
276 ALLOCATE(sph_work%WAR2(9,nsphr))
277 END IF
278C-------
279 kvnorm =16*numsph+1
280C-----------------------------------------------
281C old density storage.
282 DO n=itask+1,numsph,nthread
283 wa(kwasph*(n-1)+10)=spbuf(2,n)
284 ENDDO
285
286 IF( (glob_therm%ITHERM/=0) .OR. (glob_therm%ITHERM_FE/=0)) THEN
287 IF(itask==0)THEN
288 ALLOCATE(sph_work%WT(numsph))
289 ALLOCATE(sph_work%WGRADT(3*numsph))
290 ALLOCATE(sph_work%WLAPLT(numsph))
291 ALLOCATE(sph_work%LAMBDA(numsph))
292 ALLOCATE(sph_work%WGRADTSM(3*nsphsym))
293 END IF
294 ngdone = 1
295C /---------------/
296 CALL my_barrier
297C /---------------/
298C
299 50 CONTINUE
300#include "lockon.inc"
301 IF(ngdone>ngroup) THEN
302#include "lockoff.inc"
303 GOTO 51
304 ENDIF
305 ng=ngdone
306 ngdone = ng + 1
307#include "lockoff.inc"
308C--------
309 IF(iparg(8,ng)==1)GOTO 50
310 IF (iddw>0) CALL startimeg(ng)
311 DO nelem = 1,iparg(2,ng),nvsiz
312 offset = nelem - 1
313 nsg =iparg(10,ng)
314 nvc =iparg(19,ng)
315 CALL initbuf(iparg ,ng ,
316 2 mtn ,nel ,nft ,iad ,ity ,
317 3 npt ,jale ,ismstr ,jeul ,jtur ,
318 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
319 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
320 6 irep ,iint ,igtyp ,israt ,isrot ,
321 7 icsen ,isorth ,isorthg ,ifailure,jsms )
322 lft=1
323 llt=min(nvsiz,nel-nelem+1)
324 IF(ity==51) THEN
325C-----------------------------------------------
326 gbuf => elbuf_tab(ng)%GBUF
327 IF(jthe > 0)THEN
328 DO i=lft,llt
329 n=nft+i
330 IF(kxsp(2,n)>0)THEN
331 sph_work%WT(n)=gbuf%TEMP(i)
332 mx =ipart(1,ipartsp(n))
333 IF(sph_work%WT(n)<=pm(80,mx))THEN
334 sph_work%LAMBDA(n)=pm(75,mx)+pm(76,mx)*sph_work%WT(n)
335 ELSE
336 sph_work%LAMBDA(n)=pm(77,mx)+pm(78,mx)*sph_work%WT(n)
337 END IF
338 END IF
339 END DO
340 ELSEIF (jthe < 0) THEN
341 DO i=lft,llt
342 n=nft+i
343 IF(kxsp(2,n)>0)THEN
344 inod = kxsp(3,n)
345 sph_work%WT(n)=temp(inod)
346 mx =ipart(1,ipartsp(n))
347 IF(sph_work%WT(n)<=pm(80,mx))THEN
348 sph_work%LAMBDA(n)=pm(75,mx)+pm(76,mx)*sph_work%WT(n)
349 ELSE
350 sph_work%LAMBDA(n)=pm(77,mx)+pm(78,mx)*sph_work%WT(n)
351 END IF
352 sph_work%LAMBDA(n)=sph_work%LAMBDA(n)*glob_therm%THEACCFACT
353 END IF
354 END DO
355 ELSE
356 DO i=lft,llt
357 n=nft+i
358 sph_work%WT(n) =zero
359 sph_work%LAMBDA(n)=zero
360 END DO
361 END IF
362C-----------------------------------------------
363 ENDIF
364 IF (iddw>0) CALL stoptimeg(ng)
365 END DO
366 GOTO 50
367C-------
368 51 CONTINUE
369C
370 IF(nspmd>1) THEN
371C /---------------/
372 CALL my_barrier
373C /---------------/
374 IF(itask==0) THEN
375 CALL startime(timers,93)
376 CALL spmd_sphgett(sph_work%WT,sph_work%WTR,sph_work%LAMBDA,sph_work%LAMBDR)
377 CALL stoptime(timers,93)
378 END IF
379 END IF
380C
381C /---------------/
382 CALL my_barrier
383C /---------------/
384C
385 ngdone = 1
386C
387C /---------------/
388 CALL my_barrier
389C /---------------/
390C
391 60 CONTINUE
392#include "lockon.inc"
393 IF(ngdone>ngroup) THEN
394#include "lockoff.inc"
395 GOTO 61
396 ENDIF
397 ng=ngdone
398 ngdone = ng + 1
399#include "lockoff.inc"
400C--------
401 IF(iparg(8,ng)==1)GOTO 60
402 IF (iddw>0) CALL startimeg(ng)
403 DO nelem = 1,iparg(2,ng),nvsiz
404 offset = nelem - 1
405 nsg =iparg(10,ng)
406 nvc =iparg(19,ng)
407 CALL initbuf(iparg ,ng ,
408 2 mtn ,nel ,nft ,iad ,ity ,
409 3 npt ,jale ,ismstr ,jeul ,jtur ,
410 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
411 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
412 6 irep ,iint ,igtyp ,israt ,isrot ,
413 7 icsen ,isorth ,isorthg ,ifailure,jsms )
414 lft=1
415 llt=min(nvsiz,nel-nelem+1)
416 IF(ity==51.AND.jthe/=0) THEN
417C-----------------------------------------------
418 CALL spgradt(
419 1 x ,ms ,spbuf ,kxsp ,ixsp ,
420 2 nod2sp ,ispsym ,xspsym ,wa ,wasph ,
421 3 sph_work%WT,sph_work%WTR,sph_work%WGRADT , lft, llt, nft)
422C-----------------------------------------------
423 ENDIF
424 IF (iddw>0) CALL stoptimeg(ng)
425 END DO
426 GOTO 60
427C-------
428 61 CONTINUE
429C
430C /---------------/
431 CALL my_barrier
432C /---------------/
433C
434 IF(nspmd>1) THEN
435 IF(itask==0) THEN
436 CALL startime(timers,93)
437 CALL spmd_sphgetg(sph_work%WGRADT,wasph,sph_work%WGR,sph_iord1)
438 CALL stoptime(timers,93)
439 END IF
440C /---------------/
441 CALL my_barrier
442C /---------------/
443 END IF
444C
445 ngdone = 1
446C
447C-----------------------------------------------
448C PREPARE KERNEL CORRECTIONS FOR SYMMETRIC PARTICLES.
449C-----------------------------------------------
450 CALL spscomp(
451 1 ispsym ,wasph ,ispcond ,xframe ,wsmcomp,
452 2 geo ,ipart ,ipartsp ,waspact ,itask )
453C /---------------/
454 CALL my_barrier
455C /---------------/
456 IF(itask==0)
457 1 CALL spgtsym(
458 1 ispcond, xframe, ispsym, xspsym,
459 2 sph_work%WGRADT, sph_work%WGRADTSM,waspact, sph_work%WGR,
460 3 lft, llt, nft)
461C /---------------/
462 CALL my_barrier
463C /---------------/
464C
465 70 CONTINUE
466#include "lockon.inc"
467 IF(ngdone>ngroup) THEN
468#include "lockoff.inc"
469 GOTO 71
470 ENDIF
471 ng=ngdone
472 ngdone = ng + 1
473#include "lockoff.inc"
474C--------
475 IF(iparg(8,ng)==1)GOTO 70
476 IF (iddw>0) CALL startimeg(ng)
477 DO nelem = 1,iparg(2,ng),nvsiz
478 offset = nelem - 1
479 nsg =iparg(10,ng)
480 nvc =iparg(19,ng)
481 CALL initbuf(iparg ,ng ,
482 2 mtn ,nel ,nft ,iad ,ity ,
483 3 npt ,jale ,ismstr ,jeul ,jtur ,
484 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
485 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
486 6 irep ,iint ,igtyp ,israt ,isrot ,
487 7 icsen ,isorth ,isorthg ,ifailure,jsms )
488 lft=1
489 llt=min(nvsiz,nel-nelem+1)
490 IF(ity==51.AND.jthe==1) THEN
491C-----------------------------------------------
492 CALL splaplt(
493 1 x ,ms ,spbuf ,kxsp ,ixsp ,
494 2 nod2sp ,ispsym ,xspsym ,wa ,wasph ,
495 3 sph_work%WGRADT ,sph_work%WGR ,sph_work%WGRADTSM ,sph_work%WLAPLT ,wsmcomp,
496 4 sph_work%LAMBDA ,sph_work%LAMBDR, lft, llt, nft )
497C-----------------------------------------------
498 gbuf => elbuf_tab(ng)%GBUF
499 DO i=lft,llt
500 n=nft+i
501 IF(kxsp(2,n)>0)THEN
502 inod =kxsp(3,n)
503 rhoi =spbuf(2,n)
504 vi =spbuf(12,n)/max(em20,rhoi)
505 gbuf%EINT(i) = gbuf%EINT(i)
506 . + vi*sph_work%WLAPLT(n)*dt1/max(em20,gbuf%VOL(i))
507 END IF
508 END DO
509 ELSEIF(ity==51.AND.jthe==-1)THEN
510C-----------------------------------------------
511 CALL splaplt(
512 1 x ,ms ,spbuf ,kxsp ,ixsp ,
513 2 nod2sp ,ispsym ,xspsym ,wa ,wasph ,
514 3 sph_work%WGRADT ,sph_work%WGR ,sph_work%WGRADTSM ,sph_work%WLAPLT ,wsmcomp,
515 4 sph_work%LAMBDA ,sph_work%LAMBDR ,lft,llt,nft )
516C-----------------------------------------------
517 gbuf => elbuf_tab(ng)%GBUF
518 DO i=lft,llt
519 n=nft+i
520 IF(kxsp(2,n)>0)THEN
521 myadrn =kwasph*(n-1)
522 inod =kxsp(3,n)
523 rhoi =spbuf(2,n)
524 vi =spbuf(12,n)/max(em20,rhoi)
525 wa(myadrn+15) = vi*sph_work%WLAPLT(n)*dt1
526 END IF
527 END DO
528C-----------------------------------------------
529 ENDIF
530 IF (iddw>0) CALL stoptimeg(ng)
531 END DO
532 GOTO 70
533C-------
534 71 CONTINUE
535C
536C /---------------/
537 CALL my_barrier
538C /---------------/
539C
540 IF(itask==0) DEALLOCATE(sph_work%WT, sph_work%WGRADT, sph_work%WLAPLT, sph_work%LAMBDA, sph_work%WGRADTSM)
541
542 END IF
543C-----------------------------------------------
544C
545 ngdone = 1
546C
547C /---------------/
548 CALL my_barrier
549C /---------------/
550C
551C-------
552100 CONTINUE
553#include "lockon.inc"
554 IF(ngdone>ngroup) THEN
555#include "lockoff.inc"
556 GOTO 101
557 ENDIF
558 ng=ngdone
559 ngdone = ng + 1
560#include "lockoff.inc"
561C--------
562 IF(iparg(8,ng)==1)GOTO 100
563 IF (iddw>0) CALL startimeg(ng)
564 DO nelem = 1,iparg(2,ng),nvsiz
565 offset = nelem - 1
566 nel =iparg(2,ng)
567 nft =iparg(3,ng) + offset
568 iad =iparg(4,ng)
569 ity =iparg(5,ng)
570 lft=1
571 llt=min(nvsiz,nel-nelem+1)
572 isph2sol=iparg(69,ng)
573 IF(ity==51) THEN
574 CALL spdens(
575 1 x ,v ,ms ,spbuf ,itab ,
576 2 kxsp ,ixsp ,nod2sp ,ispsym ,xspsym ,
577 3 vspsym ,iparg ,wa ,wasph )
578 ENDIF
579 IF (iddw>0) CALL stoptimeg(ng)
580 END DO
581 GOTO 100
582 101 CONTINUE
583C
584C /---------------/
585 CALL my_barrier
586C /---------------/
587C
588 IF(itask==0)THEN
589C-------
590C-----------------------------------------------
591C Inlets/outlets : Re-impose rho.
592C-----------------------------------------------
593 IF(nsphio/=0)THEN
594C Comm wa and wa comp on remote cells before treatment with.Sym.
595C we potentially need wa remote
596 IF(nspmd>1)THEN
597 CALL startime(timers,93)
598 CALL spmd_sphgetwa(wa,sph_work%WAR2,kxsp)
599 CALL stoptime(timers,93)
600 ENDIF
601
602 CALL sponfro(x ,v ,a ,ms ,
603 2 spbuf ,itab ,kxsp ,ixsp ,nod2sp ,
604 3 isphio ,ipart ,ipartsp ,waspact ,wa ,
605 4 wasph(kvnorm), sph_work%WAR2 )
606
607 ENDIF
608 ENDIF
609C
610 ngdone = 1
611C
612C /---------------/
613 CALL my_barrier
614C /---------------/
615C
616C-----------------------------------------------
617C Pressure = F (Densite): Material laws.
618C-----------------------------------------------
619250 CONTINUE
620#include "lockon.inc"
621 IF(ngdone>ngroup) THEN
622#include "lockoff.inc"
623 GOTO 251
624 ENDIF
625 ng=ngdone
626 ngdone = ng + 1
627#include "lockoff.inc"
628C--------
629 IF(iparg(8,ng)==1)GOTO 250
630 IF (iddw>0) CALL startimeg(ng)
631 DO nelem = 1,iparg(2,ng),nvsiz
632 offset = nelem - 1
633 nsg =iparg(10,ng)
634 nvc =iparg(19,ng)
635 CALL initbuf(iparg ,ng ,
636 2 mtn ,nel ,nft ,iad ,ity ,
637 3 npt ,jale ,ismstr ,jeul ,jtur ,
638 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
639 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
640 6 irep ,iint ,igtyp ,israt ,isrot ,
641 7 icsen ,isorth ,isorthg ,ifailure,jsms )
642 lft=1
643 llt=min(nvsiz,nel-nelem+1)
644 IF(ity==51) THEN
645 jsph=1
646 jcvt=0
647C full plasticity correction (radial return) by default.
648 jplasol=1
649 istra = iparg(44,ng)
650 isph2sol=iparg(69,ng)
651 iexpan = iparg(49,ng)
652 ipartsph=0
653C-----------------------------------------------
654C Returns GIS, STIF and SSP in WA (1: 8.1: NUMSPH).
655C a ne pas ecraser avt SPFORC.
656C-----------------------------------------------
657 ndsvoid=0
658 CALL spstres(timers,elbuf_tab,ng ,
659 1 pm ,geo ,x ,v ,ms ,
660 2 w ,spbuf ,wa ,nloc_dmg ,
661 3 itab ,pld ,bufmat ,bufgeo ,partsav ,
662 4 fsav ,dt2t ,iparg ,npc ,kxsp ,
663 5 ixsp ,nod2sp ,neltst ,ityptst ,ipart ,
664 6 ipartsp ,fv ,nel ,ipm ,gresav ,
665 7 grth ,igrth ,table ,istra ,voln ,
666 8 igeo ,iexpan ,temp ,itask ,sph2sol ,
667 9 mat_elem ,ibid ,output ,snpc ,stf ,
668 a sbufmat, svis ,nsvois ,idtmins ,iresp,
669 . idel7ng, idel7nok ,idtmin ,maxfunc ,lipart1,
670 . imon_mat, userl_avail,impl_s,
671 v idyna, dt ,glob_therm,sensors)
672 ENDIF
673 IF (iddw>0) CALL stoptimeg(ng)
674 END DO
675 GOTO 250
676C-------
677 251 CONTINUE
678C
679C /---------------/
680 CALL my_barrier
681C /---------------/
682C
683C-----------------------------------------------
684C Download RHO & P from solids to SPH
685C-----------------------------------------------
686 IF(nsphsol/=0)THEN
687C /---------------/
688 CALL my_barrier
689C /---------------/
690 CALL soltosphp(
691 . x ,spbuf ,ixs ,kxsp ,ipartsp ,
692 . irst ,elbuf_tab,iparg ,ngrounc ,igrounc ,
693 . sol2sph,wa ,pm)
694C barrier inside SOLTOSPHP
695 END IF
696C-----------------------------------------------
697C Stabilization (SPH without a tensile instability) :: W(Dp)
698C-----------------------------------------------
699 IF(itask==0)THEN
700 ALLOCATE(sph_work%STAB(7,numsph+nsphr+nsphsym+1),stat=stat)
701 IF (stat/=0)THEN
702 END IF
703 sph_work%STAB=zero
704 END IF
705C /---------------/
706 CALL my_barrier
707C /---------------/
708 CALL spstabw(
709 1 itask ,iparg ,ngrounc ,igrounc ,kxsp ,
710 2 ispcond ,ispsym ,waspact ,sph2sol ,wa ,
711 3 sph_work%WASIGSM,sph_work%WAR ,sph_work%STAB ,ixsp ,nod2sp ,
712 4 spbuf ,x ,ipart ,ipartsp ,xspsym )
713C
714C-------
715C
716C Comm wa and wa comp on remote cells before treatment with.Sym.
717C
718 IF(nspmd>1)THEN
719C /---------------/
720 CALL my_barrier
721C /---------------/
722 IF(itask==0)THEN
723 CALL startime(timers,93)
724 CALL spmd_sphgetw(spbuf,wasph,wa,sph_work%WAR,sph_iord1)
725C To optimize (1 comm)
726 CALL spmd_sphgetstb(sph_work%STAB,sph_work%STAB(1,numsph+1))
727 CALL stoptime(timers,93)
728 ENDIF
729 END IF
730C-----------------------------------------------
731C Outlets : Re-impose P,E.
732C-----------------------------------------------
733 IF(itask==0)THEN
734 IF(nsphio/=0)THEN
735 CALL sponfprs(x ,v ,a ,ms ,
736 2 spbuf ,itab ,kxsp ,ixsp ,nod2sp ,
737 3 isphio ,vsphio ,npc ,pld ,pm ,
738 4 iparg ,elbuf_tab,ipart ,ipartsp ,waspact ,
739 5 wasph(kvnorm),wa ,sphveln ,sph_work%WAR, wfext)
740
741c on doit mettre a jour wa(1),wa(2),wa(3) suite a sponfprs
742 IF(nspmd>1) THEN
743 CALL startime(timers,93)
744 CALL spmd_sphgetw(spbuf,wasph,wa,sph_work%WAR,sph_iord1)
745 CALL stoptime(timers,93)
746 ENDIF
747 ENDIF
748 END IF
749C-----------------------------------------------
750C PREPARE SYMMETRIC PARTICLES STRESSES.
751C-----------------------------------------------
752 IF(itask==0)THEN
753 CALL spsgsym(ispcond ,xframe ,ispsym ,xspsym ,vspsym ,
754 2 wa ,sph_work%WASIGSM,waspact,sph_work%WAR )
755 ENDIF
756C------
757C
758
759C /---------------/
760 CALL my_barrier
761C /---------------/
762
763 DO ns=itask+1,nsphact,nthread
764 n=waspact(ns)
765 spbuf(11,n)=zero
766 ENDDO
767C
768C-----------------------------------------------
769C PREPARE KERNEL CORRECTIONS FOR SYMMETRIC PARTICLES.
770C-----------------------------------------------
771 IF (glob_therm%ITHERM==0)
772 . CALL spscomp(
773 1 ispsym ,wasph ,ispcond ,xframe ,wsmcomp,
774 2 geo ,ipart ,ipartsp ,waspact ,itask )
775C
776C-----------------------------------------------
777C assembly of forces on the particles
778C-----------------------------------------------
779C /---------------/
780 CALL my_barrier
781C /---------------/
782C initialization on task i
783 DO ns=itask+1,nsphact,nthread
784 n =waspact(ns)
785C (re)used for storage of FX, FY, FZ, STIF :
786 myadrn =kwasph*(n-1)
787 wa(myadrn+10)=zero
788 wa(myadrn+11)=zero
789 wa(myadrn+12)=zero
790 wa(myadrn+ 7)=zero
791 ENDDO
792C-----------------------------------------------
793C Stabilization (SPH without a tensile instability) :: Artificial Stress
794C-----------------------------------------------
795 CALL spstabs(
796 1 itask ,iparg ,ngrounc ,igrounc ,kxsp ,
797 2 ispcond ,ispsym ,waspact ,sph2sol ,wa ,
798 3 sph_work%WASIGSM,sph_work%WAR ,sph_work%STAB ,ixsp ,nod2sp ,
799 4 spbuf ,x )
800C
801C Barrier on WA & STAB compulsory, WA & STAB SHARING
802 ngdone = 1
803C /---------------/
804 CALL my_barrier
805C /---------------/
806C-------
807 350 CONTINUE
808#include "lockon.inc"
809 IF(ngdone>ngroup) THEN
810#include "lockoff.inc"
811 GOTO 351
812 ENDIF
813 ng=ngdone
814 ngdone = ng + 1
815#include "lockoff.inc"
816C--------
817C Solids to SPH, particles must be computed when cloud active
818 IF(iparg(8,ng)==1)GOTO 350
819 IF (iddw>0) CALL startimeg(ng)
820 DO nelem = 1,iparg(2,ng),nvsiz
821 offset = nelem - 1
822 nel =iparg(2,ng)
823 nft =iparg(3,ng) + offset
824 iad =iparg(4,ng)
825 ity =iparg(5,ng)
826 isph2sol=iparg(69,ng)
827 ipartsph=0
828 lft=1
829 llt=min(nvsiz,nel-nelem+1)
830 IF(ity==51) THEN
831C-----------
832 CALL spforcp(
833 1 pm ,geo ,x ,v ,ms ,
834 2 spbuf ,itab ,pld ,bufmat ,bufgeo ,
835 3 partsav ,fsav ,dt2t ,iparg ,npc ,
836 4 kxsp ,ixsp ,nod2sp ,neltst ,ityptst ,
837 5 ipart ,ipartsp ,ispcond ,xframe ,ispsym ,
838 6 xspsym ,vspsym ,wa ,sph_work%WASIGSM,wasph ,
839 7 wsmcomp,waspact,sph_work%WAR ,sph_work%STAB, wfext)
840 ENDIF
841 IF (iddw>0) CALL stoptimeg(ng)
842 END DO
843 GOTO 350
844C--------
845 351 CONTINUE
846C
847 nisky_l = nisky
848C /---------------/
849 CALL my_barrier
850C /---------------/
851C--------
852 IF(nsphsol==0)THEN
853 IF (glob_therm%ITHERM_FE > 0)THEN
854 IF(iparit==0)THEN
855 DO ns=itask+1,nsphact,nthread
856 n=waspact(ns)
857 myadrn =kwasph*(n-1)
858 inod=kxsp(3,n)
859 a(1,inod)=a(1,inod)+wa(myadrn+10)
860 a(2,inod)=a(2,inod)+wa(myadrn+11)
861 a(3,inod)=a(3,inod)+wa(myadrn+12)
862 stifn(inod)=stifn(inod)+wa(myadrn+7)
863 fthe(inod)=fthe(inod)+wa(myadrn+15)
864 ENDDO
865 ELSE
866 DO ns=itask+1,nsphact,nthread
867 n=waspact(ns)
868 myadrn =kwasph*(n-1)
869 inod=kxsp(3,n)
870 fskyi(nisky_l+ns,1)=wa(myadrn+10)
871 fskyi(nisky_l+ns,2)=wa(myadrn+11)
872 fskyi(nisky_l+ns,3)=wa(myadrn+12)
873 fskyi(nisky_l+ns,4)=wa(myadrn+7)
874 ftheskyi(nisky_l+ns)=wa(myadrn+15)
875 isky(nisky_l+ns) =inod
876 ENDDO
877 IF(itask==0) nisky = nisky + nsphact
878 ENDIF
879 ELSE
880 IF(iparit==0)THEN
881 DO ns=itask+1,nsphact,nthread
882 n=waspact(ns)
883 myadrn =kwasph*(n-1)
884 inod=kxsp(3,n)
885 a(1,inod)=a(1,inod)+wa(myadrn+10)
886 a(2,inod)=a(2,inod)+wa(myadrn+11)
887 a(3,inod)=a(3,inod)+wa(myadrn+12)
888 stifn(inod)=stifn(inod)+wa(myadrn+7)
889 ENDDO
890 ELSE
891 DO ns=itask+1,nsphact,nthread
892 n=waspact(ns)
893 myadrn =kwasph*(n-1)
894 inod=kxsp(3,n)
895 fskyi(nisky_l+ns,1)=wa(myadrn+10)
896 fskyi(nisky_l+ns,2)=wa(myadrn+11)
897 fskyi(nisky_l+ns,3)=wa(myadrn+12)
898 fskyi(nisky_l+ns,4)=wa(myadrn+7)
899 isky(nisky_l+ns) =inod
900 ENDDO
901 IF(itask==0) nisky = nisky + nsphact
902 ENDIF
903 ENDIF
904 ELSE
905 IF(iparit==0)THEN
906 DO ns=itask+1,nsphact,nthread
907 n=waspact(ns)
908 myadrn =kwasph*(n-1)
909 IF(sph2sol(n)==0)THEN
910 inod=kxsp(3,n)
911 a(1,inod)=a(1,inod)+wa(myadrn+10)
912 a(2,inod)=a(2,inod)+wa(myadrn+11)
913 a(3,inod)=a(3,inod)+wa(myadrn+12)
914 stifn(inod)=stifn(inod)+wa(myadrn+7)
915 ELSEIF (sol2sph_typ(sph2sol(n))==4) THEN
916C---------------
917C------ Tetra --
918C---------------
919 nsol=sph2sol(n)
920C
921 n1=ixs(2,nsol)
922 n2=ixs(4,nsol)
923 n3=ixs(7,nsol)
924 n4=ixs(6,nsol)
925C
926 ir=irst(1,n-first_sphsol+1)
927 is=irst(2,n-first_sphsol+1)
928 it=irst(3,n-first_sphsol+1)
929 nsphdir=igeo(37,ixs(10,nsol))
930C
931 ksi = a_gauss_tetra(ir,nsphdir)
932 eta = a_gauss_tetra(is,nsphdir)
933 zeta = a_gauss_tetra(it,nsphdir)
934C
935 phi1=ksi
936 phi2=eta
937 phi3=zeta
938 phi4=1-ksi-eta-zeta
939C
940 a(1,n1)=a(1,n1)+phi1*wa(myadrn+10)
941 a(2,n1)=a(2,n1)+phi1*wa(myadrn+11)
942 a(3,n1)=a(3,n1)+phi1*wa(myadrn+12)
943 stifn(n1)=stifn(n1)+phi1*wa(myadrn+7)
944
945 a(1,n2)=a(1,n2)+phi2*wa(myadrn+10)
946 a(2,n2)=a(2,n2)+phi2*wa(myadrn+11)
947 a(3,n2)=a(3,n2)+phi2*wa(myadrn+12)
948 stifn(n2)=stifn(n2)+phi2*wa(myadrn+7)
949
950 a(1,n3)=a(1,n3)+phi3*wa(myadrn+10)
951 a(2,n3)=a(2,n3)+phi3*wa(myadrn+11)
952 a(3,n3)=a(3,n3)+phi3*wa(myadrn+12)
953 stifn(n3)=stifn(n3)+phi3*wa(myadrn+7)
954
955 a(1,n4)=a(1,n4)+phi4*wa(myadrn+10)
956 a(2,n4)=a(2,n4)+phi4*wa(myadrn+11)
957 a(3,n4)=a(3,n4)+phi4*wa(myadrn+12)
958 stifn(n4)=stifn(n4)+phi4*wa(myadrn+7)
959C
960 ELSE
961C---------------
962C------ Hexa --
963C---------------
964 nsol=sph2sol(n)
965C
966 n1=ixs(2,nsol)
967 n2=ixs(3,nsol)
968 n3=ixs(4,nsol)
969 n4=ixs(5,nsol)
970 n5=ixs(6,nsol)
971 n6=ixs(7,nsol)
972 n7=ixs(8,nsol)
973 n8=ixs(9,nsol)
974C
975 ir=irst(1,n-first_sphsol+1)
976 is=irst(2,n-first_sphsol+1)
977 it=irst(3,n-first_sphsol+1)
978 nsphdir=nint((sol2sph(2,nsol)-sol2sph(1,nsol))**third)
979C
980 ksi = a_gauss(ir,nsphdir)
981 eta = a_gauss(is,nsphdir)
982 zeta = a_gauss(it,nsphdir)
983C
984 phi1=one_over_8*(one-ksi)*(one-eta)*(one-zeta)
985 phi2=one_over_8*(one-ksi)*(one-eta)*(one+zeta)
986 phi3=one_over_8*(one+ksi)*(one-eta)*(one+zeta)
987 phi4=one_over_8*(one+ksi)*(one-eta)*(one-zeta)
988 phi5=one_over_8*(one-ksi)*(one+eta)*(one-zeta)
989 phi6=one_over_8*(one-ksi)*(one+eta)*(one+zeta)
990 phi7=one_over_8*(one+ksi)*(one+eta)*(one+zeta)
991 phi8=one_over_8*(one+ksi)*(one+eta)*(one-zeta)
992C
993 a(1,n1)=a(1,n1)+phi1*wa(myadrn+10)
994 a(2,n1)=a(2,n1)+phi1*wa(myadrn+11)
995 a(3,n1)=a(3,n1)+phi1*wa(myadrn+12)
996 stifn(n1)=stifn(n1)+phi1*wa(myadrn+7)
997
998 a(1,n2)=a(1,n2)+phi2*wa(myadrn+10)
999 a(2,n2)=a(2,n2)+phi2*wa(myadrn+11)
1000 a(3,n2)=a(3,n2)+phi2*wa(myadrn+12)
1001 stifn(n2)=stifn(n2)+phi2*wa(myadrn+7)
1002
1003 a(1,n3)=a(1,n3)+phi3*wa(myadrn+10)
1004 a(2,n3)=a(2,n3)+phi3*wa(myadrn+11)
1005 a(3,n3)=a(3,n3)+phi3*wa(myadrn+12)
1006 stifn(n3)=stifn(n3)+phi3*wa(myadrn+7)
1007
1008 a(1,n4)=a(1,n4)+phi4*wa(myadrn+10)
1009 a(2,n4)=a(2,n4)+phi4*wa(myadrn+11)
1010 a(3,n4)=a(3,n4)+phi4*wa(myadrn+12)
1011 stifn(n4)=stifn(n4)+phi4*wa(myadrn+7)
1012
1013 a(1,n5)=a(1,n5)+phi5*wa(myadrn+10)
1014 a(2,n5)=a(2,n5)+phi5*wa(myadrn+11)
1015 a(3,n5)=a(3,n5)+phi5*wa(myadrn+12)
1016 stifn(n5)=stifn(n5)+phi5*wa(myadrn+7)
1017
1018 a(1,n6)=a(1,n6)+phi6*wa(myadrn+10)
1019 a(2,n6)=a(2,n6)+phi6*wa(myadrn+11)
1020 a(3,n6)=a(3,n6)+phi6*wa(myadrn+12)
1021 stifn(n6)=stifn(n6)+phi6*wa(myadrn+7)
1022
1023 a(1,n7)=a(1,n7)+phi7*wa(myadrn+10)
1024 a(2,n7)=a(2,n7)+phi7*wa(myadrn+11)
1025 a(3,n7)=a(3,n7)+phi7*wa(myadrn+12)
1026 stifn(n7)=stifn(n7)+phi7*wa(myadrn+7)
1027
1028 a(1,n8)=a(1,n8)+phi8*wa(myadrn+10)
1029 a(2,n8)=a(2,n8)+phi8*wa(myadrn+11)
1030 a(3,n8)=a(3,n8)+phi8*wa(myadrn+12)
1031 stifn(n8)=stifn(n8)+phi8*wa(myadrn+7)
1032C
1033 END IF
1034 ENDDO
1035 ELSE
1036 IF(itask==0)THEN
1037 nski=0
1038 DO ns=1,nsphact
1039 n=waspact(ns)
1040 myadrn =kwasph*(n-1)
1041 IF(sph2sol(n)==0)THEN
1042 inod=kxsp(3,n)
1043 nski=nski+1
1044 fskyi(nisky_l+nski,1)=wa(myadrn+10)
1045 fskyi(nisky_l+nski,2)=wa(myadrn+11)
1046 fskyi(nisky_l+nski,3)=wa(myadrn+12)
1047 fskyi(nisky_l+nski,4)=wa(myadrn+7)
1048 isky(nisky_l+nski) =inod
1049 ELSEIF (sol2sph_typ(sph2sol(n))==4) THEN
1050C---------------
1051C------ Tetra --
1052C---------------
1053 nsol=sph2sol(n)
1054C
1055 k1=iads(1,nsol)
1056 k2=iads(3,nsol)
1057 k3=iads(6,nsol)
1058 k4=iads(5,nsol)
1059C
1060 ir=irst(1,n-first_sphsol+1)
1061 is=irst(2,n-first_sphsol+1)
1062 it=irst(3,n-first_sphsol+1)
1063 nsphdir=igeo(37,ixs(10,nsol))
1064C
1065 ksi = a_gauss_tetra(ir,nsphdir)
1066 eta = a_gauss_tetra(is,nsphdir)
1067 zeta = a_gauss_tetra(it,nsphdir)
1068C
1069 phi1=ksi
1070 phi2=eta
1071 phi3=zeta
1072 phi4=1-ksi-eta-zeta
1073C
1074 fsky(1,k1)=fsky(1,k1)+phi1*wa(myadrn+10)
1075 fsky(2,k1)=fsky(2,k1)+phi1*wa(myadrn+11)
1076 fsky(3,k1)=fsky(3,k1)+phi1*wa(myadrn+12)
1077 fsky(4,k1)=fsky(4,k1)+phi1*wa(myadrn+7)
1078
1079 fsky(1,k2)=fsky(1,k2)+phi2*wa(myadrn+10)
1080 fsky(2,k2)=fsky(2,k2)+phi2*wa(myadrn+11)
1081 fsky(3,k2)=fsky(3,k2)+phi2*wa(myadrn+12)
1082 fsky(4,k2)=fsky(4,k2)+phi2*wa(myadrn+7)
1083
1084 fsky(1,k3)=fsky(1,k3)+phi3*wa(myadrn+10)
1085 fsky(2,k3)=fsky(2,k3)+phi3*wa(myadrn+11)
1086 fsky(3,k3)=fsky(3,k3)+phi3*wa(myadrn+12)
1087 fsky(4,k3)=fsky(4,k3)+phi3*wa(myadrn+7)
1088
1089 fsky(1,k4)=fsky(1,k4)+phi4*wa(myadrn+10)
1090 fsky(2,k4)=fsky(2,k4)+phi4*wa(myadrn+11)
1091 fsky(3,k4)=fsky(3,k4)+phi4*wa(myadrn+12)
1092 fsky(4,k4)=fsky(4,k4)+phi4*wa(myadrn+7)
1093C
1094 ELSE
1095C---------------
1096C------ Hexa --
1097C---------------
1098 nsol=sph2sol(n)
1099C
1100 k1=iads(1,nsol)
1101 k2=iads(2,nsol)
1102 k3=iads(3,nsol)
1103 k4=iads(4,nsol)
1104 k5=iads(5,nsol)
1105 k6=iads(6,nsol)
1106 k7=iads(7,nsol)
1107 k8=iads(8,nsol)
1108C
1109 ir=irst(1,n-first_sphsol+1)
1110 is=irst(2,n-first_sphsol+1)
1111 it=irst(3,n-first_sphsol+1)
1112C
1113 nsphdir=nint((sol2sph(2,nsol)-sol2sph(1,nsol))**third)
1114 ksi = a_gauss(ir,nsphdir)
1115 eta = a_gauss(is,nsphdir)
1116 zeta = a_gauss(it,nsphdir)
1117C
1118 phi1=one_over_8*(one-ksi)*(one-eta)*(one-zeta)
1119 phi2=one_over_8*(one-ksi)*(one-eta)*(one+zeta)
1120 phi3=one_over_8*(one+ksi)*(one-eta)*(one+zeta)
1121 phi4=one_over_8*(one+ksi)*(one-eta)*(one-zeta)
1122 phi5=one_over_8*(one-ksi)*(one+eta)*(one-zeta)
1123 phi6=one_over_8*(one-ksi)*(one+eta)*(one+zeta)
1124 phi7=one_over_8*(one+ksi)*(one+eta)*(one+zeta)
1125 phi8=one_over_8*(one+ksi)*(one+eta)*(one-zeta)
1126C
1127 fsky(1,k1)=fsky(1,k1)+phi1*wa(myadrn+10)
1128 fsky(2,k1)=fsky(2,k1)+phi1*wa(myadrn+11)
1129 fsky(3,k1)=fsky(3,k1)+phi1*wa(myadrn+12)
1130 fsky(4,k1)=fsky(4,k1)+phi1*wa(myadrn+7)
1131
1132 fsky(1,k2)=fsky(1,k2)+phi2*wa(myadrn+10)
1133 fsky(2,k2)=fsky(2,k2)+phi2*wa(myadrn+11)
1134 fsky(3,k2)=fsky(3,k2)+phi2*wa(myadrn+12)
1135 fsky(4,k2)=fsky(4,k2)+phi2*wa(myadrn+7)
1136
1137 fsky(1,k3)=fsky(1,k3)+phi3*wa(myadrn+10)
1138 fsky(2,k3)=fsky(2,k3)+phi3*wa(myadrn+11)
1139 fsky(3,k3)=fsky(3,k3)+phi3*wa(myadrn+12)
1140 fsky(4,k3)=fsky(4,k3)+phi3*wa(myadrn+7)
1141
1142 fsky(1,k4)=fsky(1,k4)+phi4*wa(myadrn+10)
1143 fsky(2,k4)=fsky(2,k4)+phi4*wa(myadrn+11)
1144 fsky(3,k4)=fsky(3,k4)+phi4*wa(myadrn+12)
1145 fsky(4,k4)=fsky(4,k4)+phi4*wa(myadrn+7)
1146
1147 fsky(1,k5)=fsky(1,k5)+phi5*wa(myadrn+10)
1148 fsky(2,k5)=fsky(2,k5)+phi5*wa(myadrn+11)
1149 fsky(3,k5)=fsky(3,k5)+phi5*wa(myadrn+12)
1150 fsky(4,k5)=fsky(4,k5)+phi5*wa(myadrn+7)
1151
1152 fsky(1,k6)=fsky(1,k6)+phi6*wa(myadrn+10)
1153 fsky(2,k6)=fsky(2,k6)+phi6*wa(myadrn+11)
1154 fsky(3,k6)=fsky(3,k6)+phi6*wa(myadrn+12)
1155 fsky(4,k6)=fsky(4,k6)+phi6*wa(myadrn+7)
1156
1157 fsky(1,k7)=fsky(1,k7)+phi7*wa(myadrn+10)
1158 fsky(2,k7)=fsky(2,k7)+phi7*wa(myadrn+11)
1159 fsky(3,k7)=fsky(3,k7)+phi7*wa(myadrn+12)
1160 fsky(4,k7)=fsky(4,k7)+phi7*wa(myadrn+7)
1161
1162 fsky(1,k8)=fsky(1,k8)+phi8*wa(myadrn+10)
1163 fsky(2,k8)=fsky(2,k8)+phi8*wa(myadrn+11)
1164 fsky(3,k8)=fsky(3,k8)+phi8*wa(myadrn+12)
1165 fsky(4,k8)=fsky(4,k8)+phi8*wa(myadrn+7)
1166C
1167 END IF
1168 ENDDO
1169C
1170 nisky = nisky + nski
1171C
1172 END IF ! IF(ITASK==0)THEN
1173 END IF
1174 END IF
1175C--------
1176C for artificial viscosity force work (per cell).
1177 dt05=half*dt1
1178 DO ns=itask+1,nsphact,nthread
1179 n=waspact(ns)
1180 spbuf(10,n)=spbuf(10,n)+dt05*spbuf(11,n)
1181 ENDDO
1182C-----------------------------------------------
1183C SPH gauges
1184C-----------------------------------------------
1185 CALL spgauge(lgauge ,gauge ,kxsp ,ixsp ,
1186 1 spbuf ,iparg ,elbuf_tab,ispsym ,xspsym,
1187 2 nod2sp ,x ,itask ,wa ,sph_work%WASIGSM,
1188 3 sph_work%WAR ,sphg_f6)
1189C----------------------------------
1190C
1191 ngdone = 1
1192
1193C /---------------/
1194 CALL my_barrier
1195C /---------------/
1196C
1197C SPMD DEALLOGATION TO MAKE IN MEMORY SHARED
1198C
1199 IF(itask==0) DEALLOCATE(sph_work%STAB, sph_work%WASIGSM)
1200 IF(itask==0 .AND. nspmd > 1)THEN
1201 DEALLOCATE(sph_work%WAR, sph_work%WTR, sph_work%WGR, sph_work%LAMBDR, sph_work%WAR2)
1202 END IF
1203C
1204C--------
1205 IF(nodadt==1.AND.
1206 . (idtmin(51)==1
1207 . .OR.idtmin(51)==2
1208 . .OR.idtmin(51)==5))THEN
1209400 CONTINUE
1210#include "lockon.inc"
1211 IF(ngdone>ngroup) THEN
1212#include "lockoff.inc"
1213 GOTO 401
1214 ENDIF
1215 ng=ngdone
1216 ngdone = ng + 1
1217#include "lockoff.inc"
1218C--------
1219 IF(iparg(8,ng)==1)GOTO 400
1220 IF (iddw>0) CALL startimeg(ng)
1221 DO nelem = 1,iparg(2,ng),nvsiz
1222 offset = nelem - 1
1223 CALL initbuf(iparg ,ng ,
1224 2 mtn ,nel ,nft ,kad ,ity ,
1225 3 npt ,jale ,ismstr ,jeul ,jtur ,
1226 4 jthe ,jlag ,jmult ,jhbe ,jivf ,
1227 5 nvaux ,jpor ,jcvt ,jclose ,jplasol ,
1228 6 irep ,iint ,igtyp ,israt ,isrot ,
1229 7 icsen ,isorth ,isorthg ,ifailure,jsms )
1230 lft=1
1231 llt=min(nvsiz,nel-nelem+1)
1232 IF(ity==51) THEN
1233 gbuf => elbuf_tab(ng)%GBUF
1234 DO 500 k=lft,llt
1235 n=nft+k
1236 IF(kxsp(2,n)<=0)GOTO 500
1237 inod=kxsp(3,n)
1238 adrn=kwasph*(n-1)+7
1239 dtx =dtfac1(51)*sqrt(two*ms(inod)/max(em20,wa(adrn)))
1240 IF(dtx>dtmin1(51)) GO TO 500
1241 IF(idtmin(51)==1)THEN
1242 tstop = tt
1243#include "lockon.inc"
1244 WRITE(iout,*)
1245 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SPH PARTICLE'
1246 WRITE(istdo,*)
1247 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SPH PARTICLE'
1248#include "lockoff.inc"
1249 ELSEIF(idtmin(51)==2)THEN
1250 IF (gbuf%OFF(k)/=zero)THEN
1251 gbuf%OFF(k) = zero
1252 kxsp(2,n) = 0
1253#include "lockon.inc"
1254 isphbuc =1
1255 idel7nok=1
1256 WRITE(iout,*)
1257 . ' -- DELETE SPH PARTICLE',kxsp(nisp,n)
1258 WRITE(istdo,*)
1259 . ' -- DELETE SPH PARTICLE',kxsp(nisp,n)
1260#include "lockoff.inc"
1261 END IF
1262 ELSEIF(idtmin(51)==5)THEN
1263 mstop=2
1264#include "lockon.inc"
1265 WRITE(iout,*)
1266 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SPH PARTICLE'
1267 WRITE(istdo,*)
1268 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SPH PARTICLE'
1269#include "lockoff.inc"
1270 ENDIF
1271 500 CONTINUE
1272 ENDIF
1273 END DO
1274 IF (iddw>0) CALL stoptimeg(ng)
1275 GOTO 400
1276C--------
1277 401 CONTINUE
1278C
1279C /---------------/
1280 CALL my_barrier
1281C /---------------/
1282
1283 ngdone = 1
1284C
1285 ENDIF
1286C-----------------------------------------------
1287C variable search distances
1288C After calculation of the forces (CPU optimization if CSPH).
1289C-----------------------------------------------
1290 CALL spadah(
1291 1 x ,v ,ms ,spbuf ,itab ,
1292 2 kxsp ,ixsp ,nod2sp ,wa ,waspact ,
1293 3 itask ,ipartsp ,ipart)
1294C-----------------------------------------------
1295C
1296C /---------------/
1297 CALL my_barrier
1298C /---------------/
1299C-----------
1300 RETURN
1301 END
my_real
#define my_real
Definition cppsort.cpp:32
soltosphp
subroutine soltosphp(x, spbuf, ixs, kxsp, ipartsp, irst, elbuf_tab, iparg, ngrounc, igrounc, sol2sph, wa, pm)
Definition soltosph.F:527
startimeg
subroutine startimeg(ng)
Definition timer.F:1371
stoptimeg
subroutine stoptimeg(ng)
Definition timer.F:1419
forintp
subroutine forintp(timers, pm, geo, x, a, v, ms, w, elbuf_tab, wa, fv, stifn, pld, bufmat, partsav, nloc_dmg, fsav, dt2t, iads, iparg, npc, neltst, ityptst, ipart, itab, isky, bufgeo, fskyi, xframe, kxsp, ixsp, nod2sp, ipartsp, spbuf, ispcond, ispsym, xspsym, vspsym, wasph, lprtsph, lonfsph, waspact, isphio, vsphio, sphveln, itask, ipm, gresav, grth, igrth, table, lgauge, gauge, ngrounc, igrounc, ixs, irst, sol2sph, sph2sol, fskyv, fsky, igeo, temp, fthe, ftheskyi, sphg_f6, wsmcomp, sol2sph_typ, mat_elem, output, sph_iord1, snpc, stf, sbufmat, idtmins, nsvois, iresp, maxfunc, imon_mat, userl_avail, impl_s, idyna, dt, glob_therm, sph_work, wfext, sensors)
Definition forintp.F:88
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
dt_mod
Definition dt_mod.F:80
initbuf_mod
Definition initbuf.F:160
initbuf_mod::initbuf
subroutine initbuf(iparg, ng, mtn, llt, nft, iad, ity, npt, jale, ismstr, jeul, jtur, jthe, jlag, jmult, jhbe, jivf, mid, jpor, jcvt, jclose, jpla, irep, iint, igtyp, israt, isrot, icsen, isorth, isorthg, ifailure, jsms)
Definition initbuf.F:261
nlocal_reg_mod
Definition nlocal_reg_mod.F:119
sphbox::nsphr
integer nsphr
Definition sphbox.F:83
table_mod
Definition table_mod.F:113
spadah
subroutine spadah(x, v, ms, spbuf, itab, kxsp, ixsp, nod2sp, wa, waspact, itask, ipartsp, ipart)
Definition spadah.F:37
spscomp
subroutine spscomp(ispsym, wacomp, ispcond, xframe, wsmcomp, geo, ipart, ipartsp, waspact, itask)
Definition spcompl.F:1172
spdens
subroutine spdens(x, v, ms, spbuf, itab, kxsp, ixsp, nod2sp, ispsym, xspsym, vspsym, iparg, wa, wacomp)
Definition spdens.F:36
spforcp
subroutine spforcp(pm, geo, x, v, ms, spbuf, itab, pld, bufmat, bufgeo, partsav, fsav, dt2t, iparg, npc, kxsp, ixsp, nod2sp, neltst, ityptst, ipart, ipartsp, ispcond, xframe, ispsym, xspsym, vspsym, wa, wasigsm, wacomp, wsmcomp, waspact, war, stab, wfext)
Definition spforcp.F:40
spgauge
subroutine spgauge(lgauge, gauge, kxsp, ixsp, spbuf, iparg, elbuf_tab, ispsym, xspsym, nod2sp, x, itask, wa, wasigsm, war, sphg_f6)
Definition spgauge.F:40
spmd_sphgetwa
subroutine spmd_sphgetwa(wa, war2, kxsp)
Definition spmd_sph.F:1499
spmd_sphgett
subroutine spmd_sphgett(wt, wtr, lambda, lambdr)
Definition spmd_sph.F:1198
spmd_sphgetg
subroutine spmd_sphgetg(wgradt, wacomp, wgr, sph_iord1)
Definition spmd_sph.F:1301
spmd_sphgetw
subroutine spmd_sphgetw(spbuf, wacomp, wa, war, sph_iord1)
Definition spmd_sph.F:487
spmd_sphgetstb
subroutine spmd_sphgetstb(stab, stabr)
Definition spmd_sph.F:719
sponfprs
subroutine sponfprs(x, v, a, ms, spbuf, itab, kxsp, ixsp, nod2sp, isphio, vsphio, npc, pld, pm, iparg, elbuf_tab, ipart, ipartsp, waspact, vnormal, wa, sphveln, war, wfext)
Definition sponfprs.F:39
sponfro
subroutine sponfro(x, v, a, ms, spbuf, itab, kxsp, ixsp, nod2sp, isphio, ipart, ipartsp, waspact, wa_epsd, vnormal, war2)
Definition sponfro.F:34
spsgsym
subroutine spsgsym(ispcond, xframe, ispsym, xspsym, vspsym, wa, wasigsm, waspact, war)
Definition spsgsym.F:33
spstabw
subroutine spstabw(itask, iparg, ngrounc, igrounc, kxsp, ispcond, ispsym, waspact, sph2sol, wa, wasigsm, war, stab, ixsp, nod2sp, spbuf, x, ipart, ipartsp, xspsym)
Definition spstab.F:37
spstabs
subroutine spstabs(itask, iparg, ngrounc, igrounc, kxsp, ispcond, ispsym, waspact, sph2sol, wa, wasigsm, war, stab, ixsp, nod2sp, spbuf, x)
Definition spstab.F:150
spstres
subroutine spstres(timers, elbuf_tab, ng, pm, geo, x, v, ms, w, spbuf, wa, nloc_dmg, itab, pld, bufmat, bufgeo, partsav, fsav, dt2t, iparg, npc, kxsp, ixsp, nod2sp, neltst, ityptst, ipart, ipartsp, fv, nel, ipm, gresav, grth, igrth, table, istrain, voln, igeo, iexpan, temp, itask, sph2sol, mat_elem, h3d_strain, output, snpc, stf, sbufmat, svis, nsvois, idtmins, iresp, idel7ng, idel7nok, idtmin, maxfunc, lipart1, imon_mat, userl_avail, impl_s, idyna, dt, glob_therm, sensors)
Definition spstres.F:67
spgtsym
subroutine spgtsym(ispcond, xframe, ispsym, xspsym, wgradt, wgradtsm, waspact, wgr, lft, llt, nft)
Definition sptemp.F:661
spgradt
subroutine spgradt(x, ms, spbuf, kxsp, ixsp, nod2sp, ispsym, xspsym, wa, wacomp, wtemp, wtr, wgradt, lft, llt, nft)
Definition sptemp.F:37
splaplt
subroutine splaplt(x, ms, spbuf, kxsp, ixsp, nod2sp, ispsym, xspsym, wa, wacomp, wgradt, wgr, wgradtsm, wlaplt, wsmcomp, lambda, lambdr, lft, llt, nft)
Definition sptemp.F:247
my_barrier
subroutine my_barrier
Definition machine.F:31
startime
subroutine startime(event, itask)
Definition timer.F:93
stoptime
subroutine stoptime(event, itask)
Definition timer.F:135
table_mod::ttable
Definition table_mod.F:126