OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
spbuc3.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!|| mod_spbuc3 ../engine/source/elements/sph/spbuc3.F
25!||--- called by ------------------------------------------------------
26!|| spbuc3 ../engine/source/elements/sph/spbuc3.F
27!||====================================================================
28 MODULE mod_spbuc3
29 implicit none
30 INTEGER*8, DIMENSION(:), ALLOCATABLE :: iv
31 INTEGER, DIMENSION(:), ALLOCATABLE :: kv,nv,ivs,
32 . iaux,kxspr
33 INTEGER, DIMENSION(:,:), ALLOCATABLE :: ixspr
34
35 END MODULE mod_spbuc3
36!||====================================================================
37!|| spbuc3 ../engine/source/elements/sph/spbuc3.F
38!||--- called by ------------------------------------------------------
39!|| sphtri ../engine/source/elements/sph/sphtri.F
40!||--- calls -----------------------------------------------------
41!|| ancmsg ../engine/source/output/message/message.F
42!|| arret ../engine/source/system/arret.F
43!|| my_barrier ../engine/source/system/machine.F
44!|| myqsort ../common_source/tools/sort/myqsort.F
45!|| sptrivox ../engine/source/elements/sph/sptrivox.F
46!||--- uses -----------------------------------------------------
47!|| message_mod ../engine/share/message_module/message_mod.F
48!|| mod_spbuc3 ../engine/source/elements/sph/spbuc3.F
49!|| sph_struct_mod ../engine/share/modules/sph_struct_mod.F
50!|| sphbox ../engine/share/modules/sphbox.F
51!|| tri7box ../engine/share/modules/tri7box.F
52!||====================================================================
53 SUBROUTINE spbuc3(X ,KXSP ,IXSP ,NOD2SP,NSN ,
54 2 SPBUF ,MA ,JVOIS ,JSTOR ,JPERM ,
55 3 DVOIS ,IREDUCE,BMINMA,NSNR ,NSP2SORTF,
56 4 NSP2SORTL,ITASK,KREDUCE,LGAUGE ,GAUGE )
57C============================================================================
58C M o d u l e s
59C-----------------------------------------------
60 USE tri7box
61 USE sphbox
62 USE mod_spbuc3
63 USE message_mod
64 USE sph_struct_mod
65C-----------------------------------------------
66C I m p l i c i t T y p e s
67C-----------------------------------------------
68#include "implicit_f.inc"
69C-----------------------------------------------
70C C o m m o n B l o c k s
71C-----------------------------------------------
72#include "param_c.inc"
73#include "sphcom.inc"
74#include "task_c.inc"
75C-----------------------------------------------
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 INTEGER NSN, NSNR,NSP2SORTF,NSP2SORTL
79C REAL
80 INTEGER KXSP(NISP,*), IXSP(KVOISPH,*), NOD2SP(*),
81 . ma(*), jvois(*), jstor(*), jperm(*), ireduce,itask,
82 . kreduce(*),lgauge(3,*)
83C REAL
84 my_real
85 . x(3,*),spbuf(nspbuf,*),dvois(*), bminma(12), gauge(llgauge,*)
86C-----------------------------------------------
87C L o c a l V a r i a b l e s
88C-----------------------------------------------
89 INTEGER :: I, N, INOD,JNOD,NVOIS,NVOIS1,NVOIS2,K,M,NN,NS,MS,JK,IERROR,JTASK,MY_ADRV,NVOISS
90 my_real :: MARGE, AAA, AAA2, XI,YI,ZI,XJ,YJ,ZJ, D2,D1X,D1Y,D1Z,DK
91 INTEGER NBX,NBY,NBZ
92 INTEGER (KIND=8) :: NBX8,NBY8,NBZ8,RES8,LVOXEL8
93C-----------------------------------------------
94c!!!! MARGE = TZINF-MAX(GAP,PMAX_GAP) ! il s agit bien de la marge
95
96 aaa = sqrt(nsn /
97 . ((bminma(7)-bminma(10))*(bminma(8)-bminma(11))
98 . +(bminma(8)-bminma(11))*(bminma(9)-bminma(12))
99 . +(bminma(9)-bminma(12))*(bminma(7)-bminma(10))))
100
101 aaa = 0.75*aaa
102
103 nbx = nint(aaa*(bminma(7)-bminma(10)))
104 nby = nint(aaa*(bminma(8)-bminma(11)))
105 nbz = nint(aaa*(bminma(9)-bminma(12)))
106 nbx = max(nbx,1)
107 nby = max(nby,1)
108 nbz = max(nbz,1)
109
110 nbx8=nbx
111 nby8=nby
112 nbz8=nbz
113 res8=(nbx8+2)*(nby8+2)*(nbz8+2)
114 lvoxel8 = lvoxel
115
116 IF(res8 > lvoxel8) THEN
117 aaa = lvoxel
118 aaa = aaa/((nbx8+2)*(nby8+2)*(nbz8+2))
119 aaa = aaa**(third)
120 nbx = int((nbx+2)*aaa)-2
121 nby = int((nby+2)*aaa)-2
122 nbz = int((nbz+2)*aaa)-2
123 nbx = max(nbx,1)
124 nby = max(nby,1)
125 nbz = max(nbz,1)
126 ENDIF
127
128 nbx8=nbx
129 nby8=nby
130 nbz8=nbz
131 res8=(nbx8+2)*(nby8+2)*(nbz8+2)
132
133 IF(res8 > lvoxel8) THEN
134 nbx = min(100,max(nbx8,1))
135 nby = min(100,max(nby8,1))
136 nbz = min(100,max(nbz8,1))
137 ENDIF
138
139C initialisation complete de VOXEL
140C (en // SMP il y a possibilite de redondance de traitement mais no pb)
141 DO i=inivoxel,(nbx+2)*(nby+2)*(nbz+2)
142 voxel1(i)=0
143 ENDDO
144 inivoxel = max(inivoxel,(nbx+2)*(nby+2)*(nbz+2)+1)
145
146 IF(itask==0)THEN
147 ALLOCATE(kxspr(nsnr),ixspr(kvoisph,nsnr),stat=ierror)
148 IF(ierror/=0)THEN
149 CALL ancmsg(msgid=248,anmode=aninfo_blind)
150 CALL arret(2)
151 END IF
152 kxspr(1:nsnr)=0
153 END IF
154C
155 CALL my_barrier
156C--------------------------------------------------
157C VOXEL SORT
158C--------------------------------------------------
159 CALL sptrivox(
160 1 nsn ,nsnr ,x ,bminma ,nod2sp ,
161 2 nbx ,nby ,nbz ,marge ,itask ,
162 3 ma ,spbuf ,jvois ,jstor ,jperm ,
163 4 dvois ,ireduce,nsp2sortf,nsp2sortl,voxel1 ,
164 5 kxsp ,ixsp ,kreduce ,lgauge ,gauge ,
165 6 kxspr ,ixspr )
166
167C--------------------------------------------------
168C SYMETRISE VOISINS
169C--------------------------------------------------
170c CALL MY_BARRIER fait dans SPTRIVOX
171 IF(itask==0)THEN
172
173 ALLOCATE(nv(numsph*nthread),iv(numsph*nthread),ivs(numsph),
174 . kv(numsph),iaux(kvoisph*(numsph+nsnr)),
175 . stat=ierror)
176
177 IF(ierror/=0)THEN
178 CALL ancmsg(msgid=248,anmode=aninfo_blind)
179 CALL arret(2)
180 END IF
181 END IF
182C
183 CALL my_barrier
184C
185 nv(itask*numsph+1:(itask+1)*numsph)=0
186
187 DO ns=nsp2sortf,nsp2sortl
188 n =ma(ns)
189 kv(n) =ns
190 END DO
191C-------
192C
193 CALL my_barrier
194C
195C-------
196 my_adrv=itask*numsph
197 DO ns=nsp2sortf,nsp2sortl
198 n =ma(ns)
199 inod=kxsp(3,n)
200 nvois=kxsp(5,n)
201 DO k=1,nvois
202 jnod=ixsp(k,n)
203 IF(jnod > 0) THEN
204 m=nod2sp(jnod)
205 IF(spbuf(1,m) < spbuf(1,n) .OR.
206 . (spbuf(1,m) == spbuf(1,n) .AND. kxsp(8,m) < kxsp(8,n)) )THEN
207
208 ms=kv(m)
209 nv(my_adrv+ms)=nv(my_adrv+ms)+1
210
211 END IF
212 ELSE
213C on ne fait rien pour le remote
214 ENDIF
215 END DO
216 END DO
217
218 DO n=itask+1,nsnr,nthread
219 nvois=kxspr(n)
220 DO k=1,nvois
221 jnod=ixspr(k,n)
222 IF(jnod > 0) THEN
223 m=nod2sp(jnod)
224 ELSE
225 print *,'internal error'
226 END IF
227
228 ms=kv(m)
229 nv(my_adrv+ms)=nv(my_adrv+ms)+1
230
231 END DO
232 END DO
233C-------
234C
235 CALL my_barrier
236C
237C-------
238 IF(itask==0)THEN
239
240 iv(1)=1
241 DO ns=1,nsp2sort-1
242 iv(ns+1)=iv(ns)
243 DO jtask=1,nthread
244 iv(ns+1)=iv(ns+1)+nv((jtask-1)*numsph+ns)
245 END DO
246 END DO
247
248 END IF
249C-------
250C
251 CALL my_barrier
252C
253C-------
254 DO ns=nsp2sortf,nsp2sortl
255 ivs(ns)=iv(ns)
256 END DO
257
258 DO jtask=1,nthread-1
259 DO ns=nsp2sortf,nsp2sortl
260 iv(jtask*numsph+ns)=iv((jtask-1)*numsph+ns)+nv((jtask-1)*numsph+ns)
261 END DO
262 END DO
263C-------
264C
265 CALL my_barrier
266C
267C-------
268 DO ns=nsp2sortf,nsp2sortl
269 n =ma(ns)
270 inod=kxsp(3,n)
271 nvois=kxsp(5,n)
272 DO k=1,nvois
273 jnod=ixsp(k,n)
274 IF(jnod > 0) THEN
275 m=nod2sp(jnod)
276 IF(spbuf(1,m) < spbuf(1,n) .OR.
277 . (spbuf(1,m) == spbuf(1,n) .AND. kxsp(8,m) < kxsp(8,n)) )THEN
278
279 ms=kv(m)
280 iaux(iv(itask*numsph+ms))=inod
281 iv(itask*numsph+ms)=iv(itask*numsph+ms)+1
282
283 END IF
284 ELSE
285C on ne fait rien pour le remote
286 ENDIF
287 END DO
288 END DO
289
290 DO n=itask+1,nsnr,nthread
291 nvois=kxspr(n)
292 DO k=1,nvois
293 jnod=ixspr(k,n)
294 IF(jnod > 0) THEN
295 m=nod2sp(jnod)
296 ELSE
297 print *,'internal error'
298 END IF
299
300 ms=kv(m)
301 iaux(iv(itask*numsph+ms))=-n
302 iv(itask*numsph+ms)=iv(itask*numsph+ms)+1
303
304 END DO
305 END DO
306C-------
307C
308 CALL my_barrier
309C
310C-------
311 DO ns=nsp2sortf,nsp2sortl
312 n =ma(ns)
313
314 inod=kxsp(3,n)
315 xi = x(1,inod)
316 yi = x(2,inod)
317 zi = x(3,inod)
318
319 nvois=kxsp(5,n)
320 DO k=1,nvois
321 jnod=ixsp(k,n)
322 IF(jnod > 0) THEN
323 m=nod2sp(jnod)
324 xj = x(1,jnod)
325 yj = x(2,jnod)
326 zj = x(3,jnod)
327 aaa = spbuf(1,n)+spbuf(1,m)
328 ELSE
329 nn=-jnod
330 xj = xsphr(3,nn)
331 yj = xsphr(4,nn)
332 zj = xsphr(5,nn)
333 aaa = spbuf(1,n)+xsphr(2,nn)
334 ENDIF
335
336 aaa2 = aaa*aaa
337
338C a faire : symetrie parfaite <=> X(hmin)-X(hmax)
339 d1x = xi - xj
340 d1y = yi - yj
341 d1z = zi - zj
342 d2 = d1x*d1x+d1y*d1y+d1z*d1z
343
344 jvois(k)=jnod
345 dvois(k)=d2/aaa2
346
347 END DO
348
349 nvoiss=0
350 DO jtask=1,nthread
351 nvoiss=nvoiss+nv((jtask-1)*numsph+ns)
352 END DO
353 nvois=nvois+nvoiss
354 IF(nvois>kvoisph)THEN
355 ireduce=1
356 kreduce(n)=1
357 END IF
358
359 DO k=1,nvoiss
360 jnod=iaux(ivs(ns)+k-1)
361 IF(jnod > 0) THEN
362 m=nod2sp(jnod)
363 xj = x(1,jnod)
364 yj = x(2,jnod)
365 zj = x(3,jnod)
366 aaa = spbuf(1,n)+spbuf(1,m)
367 ELSE
368 nn=-jnod
369 xj = xsphr(3,nn)
370 yj = xsphr(4,nn)
371 zj = xsphr(5,nn)
372 aaa = spbuf(1,n)+xsphr(2,nn)
373 ENDIF
374
375 aaa2 = aaa*aaa
376
377C a faire : symetrie parfaite <=> X(hmin)-X(hmax)
378 d1x = xi - xj
379 d1y = yi - yj
380 d1z = zi - zj
381 d2 = d1x*d1x+d1y*d1y+d1z*d1z
382
383 jvois(kxsp(5,n)+k)=jnod
384 dvois(kxsp(5,n)+k)=d2/aaa2
385
386 END DO
387
388 bool_sph_sort(n) = .false.
389 IF(kreduce(n)/=0 .AND. nvois > kvoisph) THEN
390 bool_sph_sort(n) = .true.
391
392 CALL myqsort(nvois,dvois,jperm,ierror)
393 DO k=1,nvois
394 jstor(k)=jvois(k)
395 ENDDO
396 DO k=1,kvoisph
397 jvois(k)=jstor(jperm(k))
398 ENDDO
399
400 dk=dvois(kvoisph)
401C-----------------
402C Choix des cellules a conserver tq distance < DK pour eviter pb de parith/on
403 nvois=0
404 DO k=1,kvoisph
405 IF(dvois(k)<dk)THEN
406 nvois=nvois+1
407 END IF
408 END DO
409 END IF
410
411 nvois=min(nvois,kvoisph)
412 kxsp(5,n)=nvois
413
414 nvois1=0
415 nvois2=nvois+1
416 DO k=1,nvois
417 jk =jvois(k)
418 dk =dvois(k)
419
420 IF(dk < one)THEN
421 nvois1=nvois1+1
422 ixsp(nvois1,n)=jk
423 ELSE
424 nvois2=nvois2-1
425 ixsp(nvois2,n)=jk
426 END IF
427 ENDDO
428 kxsp(4,n)=nvois1
429 IF(nvois1>lvoisph)ireduce=1
430
431 END DO
432C-------
433C
434 CALL my_barrier
435C
436 IF(itask==0) DEALLOCATE(kxspr,ixspr,iv,ivs,kv,nv,iaux)
437C--------------------------------------------------
438 RETURN
439 END
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
myqsort
subroutine myqsort(n, a, perm, error)
Definition myqsort.F:51
message_mod
Definition message_mod.F:1249
mod_spbuc3
Definition spbuc3.F:28
mod_spbuc3::iv
integer *8, dimension(:), allocatable iv
Definition spbuc3.F:30
mod_spbuc3::iaux
integer, dimension(:), allocatable iaux
Definition spbuc3.F:31
mod_spbuc3::kv
integer, dimension(:), allocatable kv
Definition spbuc3.F:31
mod_spbuc3::ixspr
integer, dimension(:,:), allocatable ixspr
Definition spbuc3.F:33
mod_spbuc3::kxspr
integer, dimension(:), allocatable kxspr
Definition spbuc3.F:31
mod_spbuc3::ivs
integer, dimension(:), allocatable ivs
Definition spbuc3.F:31
mod_spbuc3::nv
integer, dimension(:), allocatable nv
Definition spbuc3.F:31
sph_struct_mod
Definition sph_struct_mod.F:30
sph_struct_mod::bool_sph_sort
logical, dimension(:), allocatable bool_sph_sort
Definition sph_struct_mod.F:34
tri7box::voxel1
integer, dimension(lvoxel) voxel1
Definition tri7box.F:53
tri7box::inivoxel
integer inivoxel
Definition tri7box.F:53
tri7box::lvoxel
integer lvoxel
Definition tri7box.F:51
spbuc3
subroutine spbuc3(x, kxsp, ixsp, nod2sp, nsn, spbuf, ma, jvois, jstor, jperm, dvois, ireduce, bminma, nsnr, nsp2sortf, nsp2sortl, itask, kreduce, lgauge, gauge)
Definition spbuc3.F:57
sptrivox
subroutine sptrivox(nsn, x, bminma, nod2sp, nbx, nby, nbz, nlist, spbuf, jvois, jstor, jperm, dvois, ireduce, nsphactf, nsphactl, voxel, kxsp, ixsp, kreduce, ipartsp, sz_intp_dist, max_intp_dist_part, pre_search)
Definition sptrivox.F:40
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87
my_barrier
subroutine my_barrier
Definition machine.F:31