OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sph_nodseg.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!|| sph_nodseg ../engine/source/elements/sph/sph_nodseg.F
25!||--- called by ------------------------------------------------------
26!|| sphreqs ../engine/source/elements/sph/sphreq.F
27!|| sponof2 ../engine/source/elements/sph/sponof2.F
28!||====================================================================
29 SUBROUTINE sph_nodseg(XI,YI,ZI,XX,TFLAG,NP,LONFSPH,IXP,DPS,WNORMAL,FLAG)
30C-----------------------------------------------
31C I m p l i c i t T y p e s
32C-----------------------------------------------
33#include "implicit_f.inc"
34C-----------------------------------------------
35C D u m m y A r g u m e n t s
36C-----------------------------------------------
37 INTEGER TFLAG,LONFSPH(*),IXP(*),NP,FLAG
38 my_real
39 . xi,yi,zi,xx(12),dps(*),wnormal(3,*)
40C-----------------------------------------------
41C L o c a l V a r i a b l e s
42C-----------------------------------------------
43 INTEGER K,NORMSEG
44 my_real
45 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,xc,yc,zc,
46 . x12,y12,z12,x23,y23,z23,x31,y31,z31,
47 . xh,yh,zh,x1h,y1h,z1h,nxh,nyh,nzh,
48 . nn,nx1,ny1,nz1,nx2,ny2,nz2,nx3,ny3,nz3,nx4,ny4,nz4,
49 . d1,d2,d3,d4,d12,d23,d34,d41,ps,n12,
50 . nx,ny,nz
51C-----------------------------------------------
52 d2 = -huge(d2)
53 d3 = -huge(d3)
54 d4 = -huge(d4)
55 IF (flag==1) dps(np) = 1.e+20
56 normseg=0
57 x1 = xx(1)
58 y1 = xx(2)
59 z1 = xx(3)
60 x2 = xx(4)
61 y2 = xx(5)
62 z2 = xx(6)
63 x3 = xx(7)
64 y3 = xx(8)
65 z3 = xx(9)
66 x4 = xx(10)
67 y4 = xx(11)
68 z4 = xx(12)
69C----------
70 IF (tflag==0) THEN
71 xc=fourth*(x1+x2+x3+x4)
72 yc=fourth*(y1+y2+y3+y4)
73 zc=fourth*(z1+z2+z3+z4)
74 ELSE
75 xc=x3
76 yc=y3
77 zc=z3
78 ENDIF
79C----------
80C triangle 1,2,C
81 k =0
82 x12=x2-x1
83 y12=y2-y1
84 z12=z2-z1
85 x23=xc-x2
86 y23=yc-y2
87 z23=zc-z2
88 x31=x1-xc
89 y31=y1-yc
90 z31=z1-zc
91 nx1=-y12*z31+z12*y31
92 ny1=-z12*x31+x12*z31
93 nz1=-x12*y31+y12*x31
94 nn =one/max(em20,sqrt(nx1*nx1+ny1*ny1+nz1*nz1))
95 nx1=nn*nx1
96 ny1=nn*ny1
97 nz1=nn*nz1
98 d1=nx1*(xi-x1)+ny1*(yi-y1)+nz1*(zi-z1)
99 xh=xi-d1*nx1
100 yh=yi-d1*ny1
101 zh=zi-d1*nz1
102 x1h=xh-x1
103 y1h=yh-y1
104 z1h=zh-z1
105 nxh=y12*z1h-z12*y1h
106 nyh=z12*x1h-x12*z1h
107 nzh=x12*y1h-y12*x1h
108 IF(nxh*nx1+nyh*ny1+nzh*nz1>=zero)k=k+1
109 x1h=xh-x2
110 y1h=yh-y2
111 z1h=zh-z2
112 nxh=y23*z1h-z23*y1h
113 nyh=z23*x1h-x23*z1h
114 nzh=x23*y1h-y23*x1h
115 IF(nxh*nx1+nyh*ny1+nzh*nz1>=zero)k=k+1
116 x1h=xh-xc
117 y1h=yh-yc
118 z1h=zh-zc
119 nxh=y31*z1h-z31*y1h
120 nyh=z31*x1h-x31*z1h
121 nzh=x31*y1h-y31*x1h
122 IF(nxh*nx1+nyh*ny1+nzh*nz1>=zero)k=k+1
123 IF(k==3)THEN
124C XH interieur a 1,2,C
125 IF(abs(d1)<abs(dps(np)))THEN
126 dps(np)=d1
127 normseg=1
128C WPROJ(1,LONFSPH(IXP(NP)))=XH
129C WPROJ(2,LONFSPH(IXP(NP)))=YH
130C WPROJ(3,LONFSPH(IXP(NP)))=ZH
131 ENDIF
132 GOTO 101
133 ENDIF
134C
135 IF (tflag==0) THEN
136C----------
137C triangle 2,3,C
138 k =0
139 x12=x3-x2
140 y12=y3-y2
141 z12=z3-z2
142 x23=xc-x3
143 y23=yc-y3
144 z23=zc-z3
145 x31=x2-xc
146 y31=y2-yc
147 z31=z2-zc
148 nx2=-y12*z31+z12*y31
149 ny2=-z12*x31+x12*z31
150 nz2=-x12*y31+y12*x31
151 nn =one/max(em20,sqrt(nx2*nx2+ny2*ny2+nz2*nz2))
152 nx2=nn*nx2
153 ny2=nn*ny2
154 nz2=nn*nz2
155 d2=nx2*(xi-x2)+ny2*(yi-y2)+nz2*(zi-z2)
156 xh=xi-d2*nx2
157 yh=yi-d2*ny2
158 zh=zi-d2*nz2
159 x1h=xh-x2
160 y1h=yh-y2
161 z1h=zh-z2
162 nxh=y12*z1h-z12*y1h
163 nyh=z12*x1h-x12*z1h
164 nzh=x12*y1h-y12*x1h
165 IF(nxh*nx2+nyh*ny2+nzh*nz2>=zero)k=k+1
166 x1h=xh-x3
167 y1h=yh-y3
168 z1h=zh-z3
169 nxh=y23*z1h-z23*y1h
170 nyh=z23*x1h-x23*z1h
171 nzh=x23*y1h-y23*x1h
172 IF(nxh*nx2+nyh*ny2+nzh*nz2>=zero)k=k+1
173 x1h=xh-xc
174 y1h=yh-yc
175 z1h=zh-zc
176 nxh=y31*z1h-z31*y1h
177 nyh=z31*x1h-x31*z1h
178 nzh=x31*y1h-y31*x1h
179 IF(nxh*nx2+nyh*ny2+nzh*nz2>=zero)k=k+1
180 IF(k==3)THEN
181C XH interieur a 2,3,C
182 IF(abs(d2)<abs(dps(np)))THEN
183 dps(np)=d2
184 normseg=1
185C WPROJ(1,LONFSPH(IXP(NP)))=XH
186C WPROJ(2,LONFSPH(IXP(NP)))=YH
187C WPROJ(3,LONFSPH(IXP(NP)))=ZH
188 ENDIF
189 GOTO 101
190 ENDIF
191C----------
192C triangle 3,4,C
193 k =0
194 x12=x4-x3
195 y12=y4-y3
196 z12=z4-z3
197 x23=xc-x4
198 y23=yc-y4
199 z23=zc-z4
200 x31=x3-xc
201 y31=y3-yc
202 z31=z3-zc
203 nx3=-y12*z31+z12*y31
204 ny3=-z12*x31+x12*z31
205 nz3=-x12*y31+y12*x31
206 nn =one/max(em20,sqrt(nx3*nx3+ny3*ny3+nz3*nz3))
207 nx3=nn*nx3
208 ny3=nn*ny3
209 nz3=nn*nz3
210 d3=nx3*(xi-x3)+ny3*(yi-y3)+nz3*(zi-z3)
211 xh=xi-d3*nx3
212 yh=yi-d3*ny3
213 zh=zi-d3*nz3
214 x1h=xh-x3
215 y1h=yh-y3
216 z1h=zh-z3
217 nxh=y12*z1h-z12*y1h
218 nyh=z12*x1h-x12*z1h
219 nzh=x12*y1h-y12*x1h
220 IF(nxh*nx3+nyh*ny3+nzh*nz3>=zero)k=k+1
221 x1h=xh-x4
222 y1h=yh-y4
223 z1h=zh-z4
224 nxh=y23*z1h-z23*y1h
225 nyh=z23*x1h-x23*z1h
226 nzh=x23*y1h-y23*x1h
227 IF(nxh*nx3+nyh*ny3+nzh*nz3>=zero)k=k+1
228 x1h=xh-xc
229 y1h=yh-yc
230 z1h=zh-zc
231 nxh=y31*z1h-z31*y1h
232 nyh=z31*x1h-x31*z1h
233 nzh=x31*y1h-y31*x1h
234 IF(nxh*nx3+nyh*ny3+nzh*nz3>=zero)k=k+1
235 IF(k==3)THEN
236C XH interieur a 3,4,C
237 IF(abs(d3)<abs(dps(np)))THEN
238 dps(np)=d3
239 normseg=1
240C WPROJ(1,LONFSPH(IXP(NP)))=XH
241C WPROJ(2,LONFSPH(IXP(NP)))=YH
242C WPROJ(3,LONFSPH(IXP(NP)))=ZH
243 ENDIF
244 GOTO 101
245 ENDIF
246C----------
247C triangle 4,1,C
248 k =0
249 x12=x1-x4
250 y12=y1-y4
251 z12=z1-z4
252 x23=xc-x1
253 y23=yc-y1
254 z23=zc-z1
255 x31=x4-xc
256 y31=y4-yc
257 z31=z4-zc
258 nx4=-y12*z31+z12*y31
259 ny4=-z12*x31+x12*z31
260 nz4=-x12*y31+y12*x31
261 nn =one/max(em20,sqrt(nx4*nx4+ny4*ny4+nz4*nz4))
262 nx4=nn*nx4
263 ny4=nn*ny4
264 nz4=nn*nz4
265 d4=nx4*(xi-x4)+ny4*(yi-y4)+nz4*(zi-z4)
266 xh=xi-d4*nx4
267 yh=yi-d4*ny4
268 zh=zi-d4*nz4
269 x1h=xh-x4
270 y1h=yh-y4
271 z1h=zh-z4
272 nxh=y12*z1h-z12*y1h
273 nyh=z12*x1h-x12*z1h
274 nzh=x12*y1h-y12*x1h
275 IF(nxh*nx4+nyh*ny4+nzh*nz4>=zero)k=k+1
276 x1h=xh-x1
277 y1h=yh-y1
278 z1h=zh-z1
279 nxh=y23*z1h-z23*y1h
280 nyh=z23*x1h-x23*z1h
281 nzh=x23*y1h-y23*x1h
282 IF(nxh*nx4+nyh*ny4+nzh*nz4>=zero)k=k+1
283 x1h=xh-xc
284 y1h=yh-yc
285 z1h=zh-zc
286 nxh=y31*z1h-z31*y1h
287 nyh=z31*x1h-x31*z1h
288 nzh=x31*y1h-y31*x1h
289 IF(nxh*nx4+nyh*ny4+nzh*nz4>=zero)k=k+1
290 IF(k==3)THEN
291C XH interieur a 4,1,C
292 IF(abs(d4)<abs(dps(np)))THEN
293 dps(np)=d4
294 normseg=1
295C WPROJ(1,LONFSPH(IXP(NP)))=XH
296C WPROJ(2,LONFSPH(IXP(NP)))=YH
297C WPROJ(3,LONFSPH(IXP(NP)))=ZH
298 ENDIF
299 GOTO 101
300 ENDIF
301C
302 ENDIF
303C----------
304C distances aux cotes.
305 x12=x2-x1
306 y12=y2-y1
307 z12=z2-z1
308 n12=sqrt(x12*x12+y12*y12+z12*z12)
309 nn =one/max(em20,n12)
310 x12=x12*nn
311 y12=y12*nn
312 z12=z12*nn
313 ps=(xi-x1)*x12+(yi-y1)*y12+(zi-z1)*z12
314 IF(ps>=zero.AND.ps<=n12)THEN
315 xh =x1+ps*x12
316 yh =y1+ps*y12
317 zh =z1+ps*z12
318 x1h=xi-xh
319 y1h=yi-yh
320 z1h=zi-zh
321 ELSEIF(ps<zero)THEN
322 x1h=xi-x1
323 y1h=yi-y1
324 z1h=zi-z1
325 ELSE
326 x1h=xi-x2
327 y1h=yi-y2
328 z1h=zi-z2
329 ENDIF
330 d12=sqrt(x1h*x1h+y1h*y1h+z1h*z1h)
331 IF(d12<=onep001*abs(dps(np)).AND.d1/=zero)THEN
332 IF(d12<zep999*abs(dps(np)))THEN
333 normseg=1
334 ELSEIF (flag==0) THEN
335 normseg=2
336 ENDIF
337 IF(d12<abs(dps(np)))dps(np)=sign(d12,d1)
338C WPROJ(1,LONFSPH(IXP(NP)))=XI-X1H
339C WPROJ(2,LONFSPH(IXP(NP)))=YI-Y1H
340C WPROJ(3,LONFSPH(IXP(NP)))=ZI-Z1H
341 ENDIF
342C
343 x12=x3-x2
344 y12=y3-y2
345 z12=z3-z2
346 n12=sqrt(x12*x12+y12*y12+z12*z12)
347 nn =one/max(em20,n12)
348 x12=x12*nn
349 y12=y12*nn
350 z12=z12*nn
351 ps=(xi-x2)*x12+(yi-y2)*y12+(zi-z2)*z12
352 IF(ps>=zero.AND.ps<=n12)THEN
353 xh =x2+ps*x12
354 yh =y2+ps*y12
355 zh =z2+ps*z12
356 x1h=xi-xh
357 y1h=yi-yh
358 z1h=zi-zh
359 ELSEIF(ps<zero)THEN
360 x1h=xi-x2
361 y1h=yi-y2
362 z1h=zi-z2
363 ELSE
364 x1h=xi-x3
365 y1h=yi-y3
366 z1h=zi-z3
367 ENDIF
368 d23=sqrt(x1h*x1h+y1h*y1h+z1h*z1h)
369 IF(d23<=onep001*abs(dps(np)).AND.d2/=zero)THEN
370 IF(d23<zep999*abs(dps(np)))THEN
371 normseg=1
372 ELSEIF (flag==0) THEN
373 normseg=2
374 ENDIF
375 IF(d23<abs(dps(np)))dps(np)=sign(d23,d2)
376C WPROJ(1,LONFSPH(IXP(NP)))=XI-X1H
377C WPROJ(2,LONFSPH(IXP(NP)))=YI-Y1H
378C WPROJ(3,LONFSPH(IXP(NP)))=ZI-Z1H
379 ENDIF
380C
381 IF (tflag==0) THEN
382 x12=x4-x3
383 y12=y4-y3
384 z12=z4-z3
385 n12=sqrt(x12*x12+y12*y12+z12*z12)
386 nn =one/max(em20,n12)
387 x12=x12*nn
388 y12=y12*nn
389 z12=z12*nn
390 ps=(xi-x3)*x12+(yi-y3)*y12+(zi-z3)*z12
391 IF(ps>=zero.AND.ps<=n12)THEN
392 xh =x3+ps*x12
393 yh =y3+ps*y12
394 zh =z3+ps*z12
395 x1h=xi-xh
396 y1h=yi-yh
397 z1h=zi-zh
398 ELSEIF(ps<zero)THEN
399 x1h=xi-x3
400 y1h=yi-y3
401 z1h=zi-z3
402 ELSE
403 x1h=xi-x4
404 y1h=yi-y4
405 z1h=zi-z4
406 ENDIF
407 d34=sqrt(x1h*x1h+y1h*y1h+z1h*z1h)
408 IF(d34<=onep001*abs(dps(np)).AND.d3/=zero)THEN
409 IF(d34<zep999*abs(dps(np)))THEN
410 normseg=1
411 ELSEIF (flag==0) THEN
412 normseg=2
413 ENDIF
414 IF(d34<abs(dps(np)))dps(np)=sign(d34,d3)
415C WPROJ(1,LONFSPH(IXP(NP)))=XI-X1H
416C WPROJ(2,LONFSPH(IXP(NP)))=YI-Y1H
417C WPROJ(3,LONFSPH(IXP(NP)))=ZI-Z1H
418 ENDIF
419 ENDIF
420C
421 x12=x1-x4
422 y12=y1-y4
423 z12=z1-z4
424 n12=sqrt(x12*x12+y12*y12+z12*z12)
425 nn =one/max(em20,n12)
426 x12=x12*nn
427 y12=y12*nn
428 z12=z12*nn
429 ps=(xi-x4)*x12+(yi-y4)*y12+(zi-z4)*z12
430 IF(ps>=zero.AND.ps<=n12)THEN
431 xh =x4+ps*x12
432 yh =y4+ps*y12
433 zh =z4+ps*z12
434 x1h=xi-xh
435 y1h=yi-yh
436 z1h=zi-zh
437 ELSEIF(ps<zero)THEN
438 x1h=xi-x4
439 y1h=yi-y4
440 z1h=zi-z4
441 ELSE
442 x1h=xi-x1
443 y1h=yi-y1
444 z1h=zi-z1
445 ENDIF
446 d41=sqrt(x1h*x1h+y1h*y1h+z1h*z1h)
447 IF(d41<=onep001*abs(dps(np)).AND.d4/=zero)THEN
448 IF(d41<zep999*abs(dps(np)))THEN
449 normseg=1
450 ELSEIF (flag==0) THEN
451 normseg=2
452 ENDIF
453 IF(d41<abs(dps(np)))dps(np)=sign(d41,d4)
454C WPROJ(1,LONFSPH(IXP(NP)))=XI-X1H
455C WPROJ(2,LONFSPH(IXP(NP)))=YI-Y1H
456C WPROJ(3,LONFSPH(IXP(NP)))=ZI-Z1H
457 ENDIF
458 101 CONTINUE
459 IF(normseg==1)THEN
460C normale sortante au segment.
461 x12=half*(x2+x3-x1-x4)
462 y12=half*(y2+y3-y1-y4)
463 z12=half*(z2+z3-z1-z4)
464 x23=half*(x3+x4-x1-x2)
465 y23=half*(y3+y4-y1-y2)
466 z23=half*(z3+z4-z1-z2)
467 nx = y12*z23-z12*y23
468 ny =-x12*z23+z12*x23
469 nz = x12*y23-y12*x23
470 nn =one/max(em20,sqrt(nx*nx+ny*ny+nz*nz))
471 wnormal(1,lonfsph(ixp(np)))=-nn*nx
472 wnormal(2,lonfsph(ixp(np)))=-nn*ny
473 wnormal(3,lonfsph(ixp(np)))=-nn*nz
474 ELSEIF(normseg==2)THEN
475C cumule normale sortante au segment.
476 x12=half*(x2+x3-x1-x4)
477 y12=half*(y2+y3-y1-y4)
478 z12=half*(z2+z3-z1-z4)
479 x23=half*(x3+x4-x1-x2)
480 y23=half*(y3+y4-y1-y2)
481 z23=half*(z3+z4-z1-z2)
482 nx = y12*z23-z12*y23
483 ny =-x12*z23+z12*x23
484 nz = x12*y23-y12*x23
485 nn =one/max(em20,sqrt(nx*nx+ny*ny+nz*nz))
486 wnormal(1,lonfsph(ixp(np)))=
487 . wnormal(1,lonfsph(ixp(np)))-nn*nx
488 wnormal(2,lonfsph(ixp(np)))=
489 . wnormal(2,lonfsph(ixp(np)))-nn*ny
490 wnormal(3,lonfsph(ixp(np)))=
491 . wnormal(3,lonfsph(ixp(np)))-nn*nz
492 ENDIF
493
494C-----------------------------------------------
495 RETURN
496 END
my_real
#define my_real
Definition cppsort.cpp:32
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