OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sphreq.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/.
23C-----------------------------------------------
24C Build buckets from segments & particles.
25C-----------------------------------------------
26!||====================================================================
27!|| sphreqs ../engine/source/elements/sph/sphreq.F
28!||--- called by ------------------------------------------------------
29!|| sponof1 ../engine/source/elements/sph/sponof1.F
30!|| sponof2 ../engine/source/elements/sph/sponof2.F
31!||--- calls -----------------------------------------------------
32!|| sph_nodseg ../engine/source/elements/sph/sph_nodseg.F
33!||====================================================================
34 SUBROUTINE sphreqs(NSEG ,ISEG ,X ,NCEL ,IXP ,
35 2 LONFSPH ,KXSP ,NSX ,ISX ,NSY ,
36 3 ISY ,NSZ ,ISZ ,NPX ,IPX ,
37 4 NPY ,IPY ,NPZ ,IPZ ,DPS ,
38 5 WNORMAL ,DT_OLD,MASS_CUM,V,SPBUF ,
39 6 ITYPE)
40C-----------------------------------------------
41C I m p l i c i t T y p e s
42C-----------------------------------------------
43#include "implicit_f.inc"
44C-----------------------------------------------
45C C o m m o n B l o c k s
46C-----------------------------------------------
47#include "sphcom.inc"
48#include "param_c.inc"
49C-----------------------------------------------
50C D u m m y A r g u m e n t s
51C-----------------------------------------------
52 INTEGER NSEG,ISEG(NIBSPH,*),NCEL,IXP(*),
53 . LONFSPH(*),KXSP(NISP,*),
54 . NSX(0:*),ISX(*),NSY(0:*),ISY(*),NSZ(0:*),ISZ(*),
55 . NPX(0:*),IPX(*),NPY(0:*),IPY(*),NPZ(0:*),IPZ(*),ITYPE
56 my_real
57 . x(3,*),dps(*),wnormal(3,*),dt_old,mass_cum,v(3,*),spbuf(nspbuf,*)
58C-----------------------------------------------
59C L o c a l V a r i a b l e s
60C-----------------------------------------------
61 INTEGER N,IN,
62 . NBX,NBY,NBZ,NBAUX,NBAUY,NBAUZ,
63 . NBX1,NBY1,NBZ1,NBX2,NBY2,NBZ2,
64 . NBX3,NBY3,NBZ3,NBX4,NBY4,NBZ4,
65 . NP,NS,KP,KS,IB,IBX,IBY,IBZ,TFLAG,J
66 my_real
67 . dboitex,dboitey,dboitez,xi,yi,zi,nn,xx(12),dps_old
68C-----------------------------------------------
69 dps_old = -huge(dps_old)
70C if more than 1 box, each box size >= DBUCS
71C NBOX=MAX(IUN,INT((XBMAX-XBMIN)/DBUCS))
72C NBOY=MAX(IUN,INT((YBMAX-YBMIN)/DBUCS))
73C NBOZ=MAX(IUN,INT((ZBMAX-ZBMIN)/DBUCS))
74 dboitex=max(em20,(xbmax-xbmin)/nbox)
75 dboitey=max(em20,(ybmax-ybmin)/nboy)
76 dboitez=max(em20,(zbmax-zbmin)/nboz)
77C-----
78 DO ib=0,nbox+1
79 npx(ib)=0
80 ENDDO
81C all points inside (0<=NBX<=NBOX,0<=NBY<=NBOY,0<=NBZ<=NBOZ)
82 DO n=1,ncel
83 in=kxsp(3,lonfsph(ixp(n)))
84 xi =x(1,in)
85 nbx=int((xi-xbmin)/dboitex)+1
86C bande NBX uniquement.
87 IF(nbx>=1.AND.nbx<=nbox+1)THEN
88 npx(nbx)=npx(nbx)+1
89 ENDIF
90 ENDDO
91 DO ib=1,nbox+1
92 npx(ib)=npx(ib)+npx(ib-1)
93 ENDDO
94 DO ib=nbox+1,1,-1
95 npx(ib)=npx(ib-1)
96 ENDDO
97 DO n=1,ncel
98 in=kxsp(3,lonfsph(ixp(n)))
99 xi =x(1,in)
100 nbx=int((xi-xbmin)/dboitex)+1
101C bande NBX uniquement.
102 IF(nbx>=1.AND.nbx<=nbox+1)THEN
103 npx(nbx)=npx(nbx)+1
104 ipx(npx(nbx))=n
105 ENDIF
106 ENDDO
107C-----
108 DO ib=0,nbox+1
109 nsx(ib)=0
110 ENDDO
111C all points inside (0<=NBX<=NBOX,0<=NBY<=NBOY,0<=NBZ<=NBOZ)
112 DO n=1,nseg
113 in=iseg(1,n)
114 xi =x(1,in)
115 nbx1=int((xi-xbmin)/dboitex)+1
116 in=iseg(2,n)
117 xi =x(1,in)
118 nbx2=int((xi-xbmin)/dboitex)+1
119 in=iseg(3,n)
120 xi =x(1,in)
121 nbx3=int((xi-xbmin)/dboitex)+1
122 in=iseg(4,n)
123 xi =x(1,in)
124 nbx4=int((xi-xbmin)/dboitex)+1
125C bandes NBX1,NBX2,NBX3,NBX4 & bandes adjacentes.
126 DO nbaux=min(nbx1,nbx2,nbx3,nbx4)-1,
127 . max(nbx1,nbx2,nbx3,nbx4)+1
128 IF(nbaux>=1.AND.nbaux<=nbox+1)THEN
129 nsx(nbaux)=nsx(nbaux)+1
130 ENDIF
131 ENDDO
132 ENDDO
133 DO ib=1,nbox+1
134 nsx(ib)=nsx(ib)+nsx(ib-1)
135 ENDDO
136 DO ib=nbox+1,1,-1
137 nsx(ib)=nsx(ib-1)
138 ENDDO
139 DO n=1,nseg
140 in=iseg(1,n)
141 xi =x(1,in)
142 nbx1=int((xi-xbmin)/dboitex)+1
143 in=iseg(2,n)
144 xi =x(1,in)
145 nbx2=int((xi-xbmin)/dboitex)+1
146 in=iseg(3,n)
147 xi =x(1,in)
148 nbx3=int((xi-xbmin)/dboitex)+1
149 in=iseg(4,n)
150 xi =x(1,in)
151 nbx4=int((xi-xbmin)/dboitex)+1
152C bandes NBX1,NBX2,NBX3,NBX4 & bandes adjacentes.
153 DO nbaux=min(nbx1,nbx2,nbx3,nbx4)-1,
154 . max(nbx1,nbx2,nbx3,nbx4)+1
155 IF(nbaux>=1.AND.nbaux<=nbox+1)THEN
156 nsx(nbaux)=nsx(nbaux)+1
157 isx(nsx(nbaux))=n
158 ENDIF
159 ENDDO
160 ENDDO
161C-----
162 DO ibx=1,nbox+1
163 DO iby=0,nboy+1
164 npy(iby)=0
165 ENDDO
166 DO ks=npx(ibx-1)+1,npx(ibx)
167 n =ipx(ks)
168 in=kxsp(3,lonfsph(ixp(n)))
169 yi =x(2,in)
170 nby=int((yi-ybmin)/dboitey)+1
171C bande NBY uniquement.
172 IF(nby>=1.AND.nby<=nboy+1)THEN
173 npy(nby)=npy(nby)+1
174 ENDIF
175 ENDDO
176 DO iby=1,nboy+1
177 npy(iby)=npy(iby)+npy(iby-1)
178 ENDDO
179 DO iby=nboy+1,1,-1
180 npy(iby)=npy(iby-1)
181 ENDDO
182 DO ks=npx(ibx-1)+1,npx(ibx)
183 n =ipx(ks)
184 in=kxsp(3,lonfsph(ixp(n)))
185 yi =x(2,in)
186 nby=int((yi-ybmin)/dboitey)+1
187C bande NBY uniquement.
188 IF(nby>=1.AND.nby<=nboy+1)THEN
189 npy(nby)=npy(nby)+1
190 ipy(npy(nby))=n
191 ENDIF
192 ENDDO
193C------
194 DO iby=0,nboy+1
195 nsy(iby)=0
196 ENDDO
197 DO ks=nsx(ibx-1)+1,nsx(ibx)
198 n=isx(ks)
199 in=iseg(1,n)
200 yi =x(2,in)
201 nby1=int((yi-ybmin)/dboitey)+1
202 in=iseg(2,n)
203 yi =x(2,in)
204 nby2=int((yi-ybmin)/dboitey)+1
205 in=iseg(3,n)
206 yi =x(2,in)
207 nby3=int((yi-ybmin)/dboitey)+1
208 in=iseg(4,n)
209 yi =x(2,in)
210 nby4=int((yi-ybmin)/dboitey)+1
211C bande NBY & bandes adjacentes.
212 DO nbauy=min(nby1,nby2,nby3,nby4)-1,
213 . max(nby1,nby2,nby3,nby4)+1
214 IF(nbauy>=1.AND.nbauy<=nboy+1)THEN
215 nsy(nbauy)=nsy(nbauy)+1
216 ENDIF
217 ENDDO
218 ENDDO
219C------
220 DO iby=1,nboy+1
221 nsy(iby)=nsy(iby)+nsy(iby-1)
222 ENDDO
223 DO iby=nboy+1,1,-1
224 nsy(iby)=nsy(iby-1)
225 ENDDO
226C------
227 DO ks=nsx(ibx-1)+1,nsx(ibx)
228 n=isx(ks)
229 in=iseg(1,n)
230 yi =x(2,in)
231 nby1=int((yi-ybmin)/dboitey)+1
232 in=iseg(2,n)
233 yi =x(2,in)
234 nby2=int((yi-ybmin)/dboitey)+1
235 in=iseg(3,n)
236 yi =x(2,in)
237 nby3=int((yi-ybmin)/dboitey)+1
238 in=iseg(4,n)
239 yi =x(2,in)
240 nby4=int((yi-ybmin)/dboitey)+1
241C bande NBY & bandes adjacentes.
242 DO nbauy=min(nby1,nby2,nby3,nby4)-1,
243 . max(nby1,nby2,nby3,nby4)+1
244 IF(nbauy>=1.AND.nbauy<=nboy+1)THEN
245 nsy(nbauy)=nsy(nbauy)+1
246 isy(nsy(nbauy))=n
247 ENDIF
248 ENDDO
249 ENDDO
250C------
251 DO iby=1,nboy+1
252C-------
253 DO ibz=0,nboz+1
254 npz(ibz)=0
255 ENDDO
256 DO ks=npy(iby-1)+1,npy(iby)
257 n =ipy(ks)
258 in=kxsp(3,lonfsph(ixp(n)))
259 zi =x(3,in)
260 nbz=int((zi-zbmin)/dboitez)+1
261C bande NBZ uniquement.
262 IF(nbz>=1.AND.nbz<=nboz+1)THEN
263 npz(nbz)=npz(nbz)+1
264 ENDIF
265 ENDDO
266 DO ibz=1,nboz+1
267 npz(ibz)=npz(ibz)+npz(ibz-1)
268 ENDDO
269 DO ibz=nboz+1,1,-1
270 npz(ibz)=npz(ibz-1)
271 ENDDO
272 DO ks=npy(iby-1)+1,npy(iby)
273 n =ipy(ks)
274 in=kxsp(3,lonfsph(ixp(n)))
275 zi =x(3,in)
276 nbz=int((zi-zbmin)/dboitez)+1
277C bande NBZ uniquement.
278 IF(nbz>=1.AND.nbz<=nboz+1)THEN
279 npz(nbz)=npz(nbz)+1
280 ipz(npz(nbz))=n
281 ENDIF
282 ENDDO
283C-------
284 DO ibz=0,nboz+1
285 nsz(ibz)=0
286 ENDDO
287 DO ks=nsy(iby-1)+1,nsy(iby)
288 n=isy(ks)
289 in=iseg(1,n)
290 zi =x(3,in)
291 nbz1=int((zi-zbmin)/dboitez)+1
292 in=iseg(2,n)
293 zi =x(3,in)
294 nbz2=int((zi-zbmin)/dboitez)+1
295 in=iseg(3,n)
296 zi =x(3,in)
297 nbz3=int((zi-zbmin)/dboitez)+1
298 in=iseg(4,n)
299 zi =x(3,in)
300 nbz4=int((zi-zbmin)/dboitez)+1
301C bande NBZ & bandes adjacentes.
302 DO nbauz=min(nbz1,nbz2,nbz3,nbz4)-1,
303 . max(nbz1,nbz2,nbz3,nbz4)+1
304 IF(nbauz>=1.AND.nbauz<=nboz+1)THEN
305 nsz(nbauz)=nsz(nbauz)+1
306 ENDIF
307 ENDDO
308 ENDDO
309 DO ibz=1,nboz+1
310 nsz(ibz)=nsz(ibz)+nsz(ibz-1)
311 ENDDO
312 DO ibz=nboz+1,1,-1
313 nsz(ibz)=nsz(ibz-1)
314 ENDDO
315 DO ks=nsy(iby-1)+1,nsy(iby)
316 n=isy(ks)
317 in=iseg(1,n)
318 zi =x(3,in)
319 nbz1=int((zi-zbmin)/dboitez)+1
320 in=iseg(2,n)
321 zi =x(3,in)
322 nbz2=int((zi-zbmin)/dboitez)+1
323 in=iseg(3,n)
324 zi =x(3,in)
325 nbz3=int((zi-zbmin)/dboitez)+1
326 in=iseg(4,n)
327 zi =x(3,in)
328 nbz4=int((zi-zbmin)/dboitez)+1
329C bande NBZ & bandes adjacentes.
330 DO nbauz=min(nbz1,nbz2,nbz3,nbz4)-1,
331 . max(nbz1,nbz2,nbz3,nbz4)+1
332 IF(nbauz>=1.AND.nbauz<=nboz+1)THEN
333 nsz(nbauz)=nsz(nbauz)+1
334 isz(nsz(nbauz))=n
335 ENDIF
336 ENDDO
337 ENDDO
338C-------
339 DO ibz=1,nboz+1
340 DO kp=npz(ibz-1)+1,npz(ibz)
341 np =ipz(kp)
342C
343 IF (itype>1) THEN
344 dps(np)= 1.e+20
345 wnormal(1,lonfsph(ixp(np)))=zero
346 wnormal(2,lonfsph(ixp(np)))=zero
347 wnormal(3,lonfsph(ixp(np)))=zero
348 in=kxsp(3,lonfsph(ixp(np)))
349 xi =x(1,in)-dt_old*v(1,in)
350 yi =x(2,in)-dt_old*v(2,in)
351 zi =x(3,in)-dt_old*v(3,in)
352C
353 DO ks=nsz(ibz-1)+1,nsz(ibz)
354 ns=isz(ks)
355 DO j=1,4
356 in=iseg(j,ns)
357 xx(3*(j-1)+1)=x(1,in)-dt_old*v(1,in)
358 xx(3*(j-1)+2)=x(2,in)-dt_old*v(2,in)
359 xx(3*(j-1)+3)=x(3,in)-dt_old*v(3,in)
360 END DO
361 IF (iseg(3,ns)/=iseg(4,ns)) THEN
362 tflag = 0
363 ELSE
364 tflag = 1
365 ENDIF
366 CALL sph_nodseg(xi,yi,zi,xx,tflag,np,lonfsph,ixp,dps,wnormal,0)
367 END DO
368 dps_old=dps(np)
369 ENDIF
370C
371 dps(np)= 1.e+20
372 wnormal(1,lonfsph(ixp(np)))=zero
373 wnormal(2,lonfsph(ixp(np)))=zero
374 wnormal(3,lonfsph(ixp(np)))=zero
375 in=kxsp(3,lonfsph(ixp(np)))
376 xi =x(1,in)
377 yi =x(2,in)
378 zi =x(3,in)
379C
380 DO ks=nsz(ibz-1)+1,nsz(ibz)
381 ns=isz(ks)
382 DO j=1,4
383 in=iseg(j,ns)
384 xx(3*(j-1)+1)=x(1,in)
385 xx(3*(j-1)+2)=x(2,in)
386 xx(3*(j-1)+3)=x(3,in)
387 END DO
388 IF (iseg(3,ns)/=iseg(4,ns)) THEN
389 tflag = 0
390 ELSE
391 tflag = 1
392 ENDIF
393 CALL sph_nodseg(xi,yi,zi,xx,tflag,np,lonfsph,ixp,dps,wnormal,0)
394 END DO
395C
396 IF (itype>1) THEN
397C detection if the particle is crossing the surface
398 IF (dps_old*dps(np)<zero) THEN
399C pur outlet/silent boundaries on permute le sens de comptage
400 IF (itype/=4) dps_old = -dps_old
401C la particule a traverse la surface dans le sens +
402 IF (dps_old>zero) mass_cum = mass_cum-spbuf(12,n)
403C la particule a traverse la surface dans le sens -
404 IF (dps_old<zero) mass_cum = mass_cum+spbuf(12,n)
405 ENDIF
406 ENDIF
407C
408 nn=wnormal(1,lonfsph(ixp(np)))*wnormal(1,lonfsph(ixp(np)))
409 . +wnormal(2,lonfsph(ixp(np)))*wnormal(2,lonfsph(ixp(np)))
410 . +wnormal(3,lonfsph(ixp(np)))*wnormal(3,lonfsph(ixp(np)))
411 nn=one/max(em20,sqrt(nn))
412 wnormal(1,lonfsph(ixp(np)))=wnormal(1,lonfsph(ixp(np)))*nn
413 wnormal(2,lonfsph(ixp(np)))=wnormal(2,lonfsph(ixp(np)))*nn
414 wnormal(3,lonfsph(ixp(np)))=wnormal(3,lonfsph(ixp(np)))*nn
415C
416 ENDDO
417 ENDDO
418 ENDDO
419 ENDDO
420C-----------------------------------------------
421 RETURN
422 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sph_nodseg
subroutine sph_nodseg(xi, yi, zi, xx, tflag, np, lonfsph, ixp, dps, wnormal, flag)
Definition sph_nodseg.F:30
sphreqs
subroutine sphreqs(nseg, iseg, x, ncel, ixp, lonfsph, kxsp, nsx, isx, nsy, isy, nsz, isz, npx, ipx, npy, ipy, npz, ipz, dps, wnormal, dt_old, mass_cum, v, spbuf, itype)
Definition sphreq.F:40