OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
horipoly.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!|| horiedge ../starter/source/airbag/horipoly.F
25!||--- called by ------------------------------------------------------
26!|| fvmesh1 ../starter/source/airbag/fvmesh.F
27!||====================================================================
28 SUBROUTINE horiedge(IPOLY , RPOLY , NX , NY , NZ ,
29 . NBNEDGE, INEDGE, RNEDGE, X0 , Y0 ,
30 . Z0 , INZ , NNS3 , NREF, AREF,
31 . NNSP )
32C-----------------------------------------------
33C I m p l i c i t T y p e s
34C-----------------------------------------------
35#include "implicit_f.inc"
36C-----------------------------------------------
37C D u m m y A r g u m e n t s
38C-----------------------------------------------
39 INTEGER IPOLY(*), NBNEDGE, INEDGE(6,*), INZ, NNS3, NREF(2,*), NNSP
40 my_real
41 . RPOLY(*), NX, NY, NZ, RNEDGE(6,*), X0, Y0, Z0, AREF(4,*)
42C-----------------------------------------------
43C L o c a l V a r i a b l e s
44C-----------------------------------------------
45 INTEGER NN, I, II
46 my_real
47 . X1, Y1, Z1, X2, Y2, Z2, ZL1, ZL2, TOLE, DD
48C
49#ifdef MYREAL8
50 tole=em10
51#else
52 tole=em5
53#endif
54 nn=ipoly(2)
55 nnsp=0
56 DO i=1,nn
57 ii=i+1
58 IF (i==nn) ii=1
59 x1=rpoly(4+3*(i-1)+1)
60 y1=rpoly(4+3*(i-1)+2)
61 z1=rpoly(4+3*(i-1)+3)
62 x2=rpoly(4+3*(ii-1)+1)
63 y2=rpoly(4+3*(ii-1)+2)
64 z2=rpoly(4+3*(ii-1)+3)
65 dd=(x1-x2)**2+(y1-y2)**2+(z1-z2)**2
66 zl1=(x1-x0)*nx+(y1-y0)*ny+(z1-z0)*nz
67 zl2=(x2-x0)*nx+(y2-y0)*ny+(z2-z0)*nz
68 IF (zl1==zero.AND.zl2==zero.AND.dd>=tole) THEN
69 nbnedge=nbnedge+1
70C
71 nnsp=nnsp+1
72 nref(1,nnsp)=ipoly(6+i)
73 nref(2,nnsp)=ipoly(6+ii)
74 aref(1,nnsp)=one
75 aref(2,nnsp)=x1
76 aref(3,nnsp)=y1
77 aref(4,nnsp)=z1
78 nnsp=nnsp+1
79 nref(1,nnsp)=ipoly(6+i)
80 nref(2,nnsp)=ipoly(6+ii)
81 aref(1,nnsp)=zero
82 aref(2,nnsp)=x2
83 aref(3,nnsp)=y2
84 aref(4,nnsp)=z2
85C
86 inedge(1,nbnedge)=ipoly(1)
87 inedge(2,nbnedge)=nns3+nnsp-1
88 inedge(3,nbnedge)=nns3+nnsp
89 inedge(4,nbnedge)=ipoly(3)
90 inedge(5,nbnedge)=ipoly(4)
91 inedge(6,nbnedge)=inz
92C
93 rnedge(1,nbnedge)=rpoly(4+3*(i-1)+1)
94 rnedge(2,nbnedge)=rpoly(4+3*(i-1)+2)
95 rnedge(3,nbnedge)=rpoly(4+3*(i-1)+3)
96 rnedge(4,nbnedge)=rpoly(4+3*(ii-1)+1)
97 rnedge(5,nbnedge)=rpoly(4+3*(ii-1)+2)
98 rnedge(6,nbnedge)=rpoly(4+3*(ii-1)+3)
99 ENDIF
100 ENDDO
101C
102 RETURN
103 END
104!||====================================================================
105!|| horipoly ../starter/source/airbag/horipoly.F
106!||--- called by ------------------------------------------------------
107!|| fvmesh1 ../starter/source/airbag/fvmesh.F
108!||====================================================================
109 SUBROUTINE horipoly(INEDGE, RNEDGE, LEDGE , NEDGE, IPOLY,
110 . RPOLY , IZ , IELNOD, NPOLY, NX ,
111 . NY , NZ , INZ , IBRIC, NEL ,
112 . TAGELA )
113C-----------------------------------------------
114C I m p l i c i t T y p e s
115C-----------------------------------------------
116#include "implicit_f.inc"
117C-----------------------------------------------
118C D u m m y A r g u m e n t s
119C-----------------------------------------------
120 INTEGER INEDGE(6,*), LEDGE(*), NEDGE, IPOLY(6+2*NEDGE+1+NEDGE,*),
121 . IZ(3,*), IELNOD(NEDGE,*), NPOLY, INZ, IBRIC, NEL,
122 . TAGELA(*)
123 my_real
124 . RNEDGE(6,*), RPOLY(4+6*NEDGE+3*NEDGE,*), NX, NY, NZ
125C-----------------------------------------------
126C L o c a l V a r i a b l e s
127C-----------------------------------------------
128 INTEGER NN, I, II, TNOD(2*NEDGE), TSEG(3,NEDGE), NN_OLD, JFOUND,
129 . j, jj, redir1(2*nedge), redir2(2*nedge), itag(2*nedge),
130 . itagseg(nedge+1), istop, iclose, n1, n2, in, nnp,
131 . poly(2*nedge,nedge), iseg, lenpoly(nedge), ifound,
132 . iadhol(nedge), nhol, k, kk, ipout, m, mm, nsec, ksmin,
133 . nedge_old, tseg_old(3,nedge), redir(nedge), jseg,
134 . jtagseg(nedge), jtag(2*nedge)
135 my_real
136 . tole, xnod(3,2*nedge), xx1, yy1, zz1, xx2, yy2, zz2,
137 . dd, xloc(2,nedge), xsec(nedge), phol(3,nedge), alpha,
138 . x1, y1, z1, vx1, vy1, vz1, vx2, vy2, vz2, nr1, nr2,
139 . ss, vvx, vvy, vvz, ss1, x2, y2, z2, ylmin, ylmax, xsmin1,
140 . xsmin2, xx, yy, zz, ylsec, xs, ys, lmax, ll
141C
142C
143#ifdef MYREAL8
144 tole=em10
145#else
146 tole=em5
147#endif
148C Creation de la liste des noeuds
149 nn=0
150 lmax=zero
151 ksmin = 0
152 DO i=1,nedge
153 ii=ledge(i)
154 nn=nn+1
155 tnod(nn)=inedge(2,ii)
156 xnod(1,nn)=rnedge(1,ii)
157 xnod(2,nn)=rnedge(2,ii)
158 xnod(3,nn)=rnedge(3,ii)
159 tseg(1,i)=nn
160 nn=nn+1
161 tnod(nn)=inedge(3,ii)
162 xnod(1,nn)=rnedge(4,ii)
163 xnod(2,nn)=rnedge(5,ii)
164 xnod(3,nn)=rnedge(6,ii)
165 tseg(2,i)=nn
166 tseg(3,i)=inedge(1,ii)
167 jj=inedge(4,ii)
168 IF(inedge(1,ii)==1 .AND. tagela(jj) > nel) THEN
169 tseg(3,i)=3
170 ENDIF
171 ll=(rnedge(1,ii)-rnedge(4,ii))**2+
172 . (rnedge(2,ii)-rnedge(5,ii))**2+
173 . (rnedge(3,ii)-rnedge(6,ii))**2
174 lmax=max(lmax,ll)
175 ENDDO
176* TOLE=TOLE*LMAX
177C Elimination des noeuds doubles
178 DO i=1,2*nedge
179 redir1(i)=0
180 redir2(i)=0
181 ENDDO
182 nn_old=nn
183 nn=0
184 DO i=1,nn_old
185 xx1=xnod(1,i)
186 yy1=xnod(2,i)
187 zz1=xnod(3,i)
188 jfound=0
189 DO j=1,nn
190 jj=redir1(j)
191 xx2=xnod(1,jj)
192 yy2=xnod(2,jj)
193 zz2=xnod(3,jj)
194 dd=sqrt((xx1-xx2)**2+(yy1-yy2)**2+(zz1-zz2)**2)
195 IF (dd<=tole) jfound=j
196 ENDDO
197 IF (jfound==0) THEN
198 nn=nn+1
199 redir1(nn)=i
200 redir2(i)=nn
201 ELSE
202 redir2(i)=jfound
203 ENDIF
204 ENDDO
205C
206 DO i=1,nedge
207 n1=tseg(1,i)
208 n2=tseg(2,i)
209 tseg(1,i)=redir2(n1)
210 tseg(2,i)=redir2(n2)
211 ENDDO
212C Elimination des segments de longueur nulle
213 nedge_old=nedge
214 nedge=0
215 DO i=1,nedge_old
216 tseg_old(1,i)=tseg(1,i)
217 tseg_old(2,i)=tseg(2,i)
218 tseg_old(3,i)=tseg(3,i)
219 ENDDO
220 DO i=1,nedge_old
221 IF (tseg_old(1,i)/=tseg_old(2,i)) THEN
222 nedge=nedge+1
223 tseg(1,nedge)=tseg_old(1,i)
224 tseg(2,nedge)=tseg_old(2,i)
225 tseg(3,nedge)=tseg_old(3,i)
226 ENDIF
227 ENDDO
228C Mise des segments internes en tete
229 j=0
230 DO i=1,nedge
231 IF(tseg(3,i) /= 3) cycle
232 j=j+1
233 redir(j)=i
234 ENDDO
235 DO i=1,nedge
236 IF(tseg(3,i) == 3) cycle
237 j=j+1
238 redir(j)=i
239 ENDDO
240C Construction des polygones
241 DO i=1,nn
242 itag(i)=0
243 jtag(i)=0
244 ENDDO
245 DO i=1,nedge
246 itagseg(i)=0
247 jtagseg(i)=1
248 j=redir(i)
249 IF(tseg(3,j)==3) jtagseg(i)=2
250 n1=tseg(1,j)
251 n2=tseg(2,j)
252 jtag(n1)=jtag(n1)+jtagseg(i)
253 jtag(n2)=jtag(n2)+jtagseg(i)
254 ENDDO
255 itagseg(nedge+1)=0
256 npoly=1
257 istop=0
258 DO WHILE (istop==0)
259 i=1
260 DO WHILE (itagseg(i)==1.AND.i<=nedge)
261 i=i+1
262 ENDDO
263 IF (i==nedge+1) THEN
264 istop=1
265 cycle
266 ENDIF
267C
268 iclose=0
269 iseg=i
270 jseg=redir(i)
271 itagseg(iseg)=1
272 n1=tseg(1,jseg)
273 n2=tseg(2,jseg)
274 itag(n1)=1
275 itag(n2)=1
276 in=n2
277 nnp=1
278 poly(1,npoly)=redir1(n1)
279 DO WHILE (iclose==0)
280 ifound=0
281 i=0
282 DO WHILE (ifound==0)
283 i=i+1
284 IF (itagseg(i) == 1) cycle
285 j=redir(i)
286 n1=tseg(1,j)
287 n2=tseg(2,j)
288 IF (n1==in) THEN
289 ifound=1
290 IF (itag(n2) == 1) iclose=1
291 iseg=i
292 in=n2
293 nnp=nnp+1
294 poly(nnp,npoly)=redir1(n1)
295 itag(n2)=1
296 ELSEIF (n2==in) THEN
297 ifound=1
298 IF (itag(n1) == 1) iclose=1
299 iseg=i
300 in=n1
301 nnp=nnp+1
302 poly(nnp,npoly)=redir1(n2)
303 itag(n1)=1
304 ENDIF
305 ENDDO
306 itagseg(iseg)=1
307 ENDDO
308C
309 IF (iclose==1) THEN
310 lenpoly(npoly)=nnp
311 npoly=npoly+1
312 DO i=1,nedge
313 jtagseg(i)=jtagseg(i)-itagseg(i)
314 ENDDO
315 DO i=1,nedge
316 itagseg(i)=0
317 IF(jtagseg(i) <= 0) itagseg(i)=1
318 ENDDO
319 DO i=1,nn
320 jtag(i)=jtag(i)-2*itag(i)
321 ENDDO
322 DO i=1,nn
323 itag(i)=0
324 IF(jtag(i) <= 0) itag(i)=1
325 ENDDO
326 ENDIF
327 ENDDO
328 npoly=npoly-1
329C Creation des tableaux de sortie
330 DO i=1,npoly
331 ipoly(1,i)=2
332 ipoly(2,i)=lenpoly(i)
333 ipoly(3,i)=ibric
334 ipoly(4,i)=ibric
335 ipoly(5,i)=0
336 ipoly(6,i)=0
337 rpoly(1,i)=zero
338 rpoly(2,i)=nx
339 rpoly(3,i)=ny
340 rpoly(4,i)=nz
341 DO j=1,lenpoly(i)
342 jj=poly(j,i)
343 ipoly(6+j,i)=-tnod(jj)
344 rpoly(4+3*(j-1)+1,i)=xnod(1,jj)
345 rpoly(4+3*(j-1)+2,i)=xnod(2,jj)
346 rpoly(4+3*(j-1)+3,i)=xnod(3,jj)
347 ielnod(j,i)=-1
348 ENDDO
349 ipoly(6+lenpoly(i)+1,i)=0
350C
351 iz(1,i)=2
352 iz(2,i)=inz
353 iz(3,i)=inz+1
354 ENDDO
355C---------------------------------------
356C Recherche des trous dans les polygones
357C---------------------------------------
358 DO i=1,npoly
359C
360 nhol=0
361 DO j=1,npoly
362 iadhol(j)=0
363 ENDDO
364C
365 DO j=1,npoly
366 IF (i==j) cycle
367 alpha=zero
368 x1=rpoly(5,j)
369 y1=rpoly(6,j)
370 z1=rpoly(7,j)
371 DO k=1,lenpoly(i)
372 kk=k+1
373 IF (k==lenpoly(i)) kk=1
374 xx1=rpoly(4+3*(k-1)+1,i)
375 yy1=rpoly(4+3*(k-1)+2,i)
376 zz1=rpoly(4+3*(k-1)+3,i)
377 xx2=rpoly(4+3*(kk-1)+1,i)
378 yy2=rpoly(4+3*(kk-1)+2,i)
379 zz2=rpoly(4+3*(kk-1)+3,i)
380 vx1=xx1-x1
381 vy1=yy1-y1
382 vz1=zz1-z1
383 vx2=xx2-x1
384 vy2=yy2-y1
385 vz2=zz2-z1
386 nr1=sqrt(vx1**2+vy1**2+vz1**2)
387 nr2=sqrt(vx2**2+vy2**2+vz2**2)
388 IF(nr1 > zero) THEN
389 vx1=vx1/nr1
390 vy1=vy1/nr1
391 vz1=vz1/nr1
392 ELSE
393 cycle
394 ENDIF
395 IF(nr2 > zero) THEN
396 vx2=vx2/nr2
397 vy2=vy2/nr2
398 vz2=vz2/nr2
399 ELSE
400 cycle
401 ENDIF
402 ss=vx1*vx2+vy1*vy2+vz1*vz2
403 vvx=vy1*vz2-vz1*vy2
404 vvy=vz1*vx2-vx1*vz2
405 vvz=vx1*vy2-vy1*vx2
406 ss1=nx*vvx+ny*vvy+nz*vvz
407 IF(ss < -one) ss=-one
408 IF(ss > one) ss= one
409 IF (ss1>=zero) THEN
410 alpha=alpha+acos(ss)
411 ELSE
412 alpha=alpha-acos(ss)
413 ENDIF
414 ENDDO
415C
416 IF (abs(alpha)>=two*pi) THEN
417C---------------------------------------------------------------
418C Le premier point du polygone j est a l'interieur du polygone i
419C On teste tous les autres
420C---------------------------------------------------------------
421 ipout=0
422 DO k=2,lenpoly(j)
423 x2=rpoly(4+3*(k-1)+1,j)
424 y2=rpoly(4+3*(k-1)+2,j)
425 z2=rpoly(4+3*(k-1)+3,j)
426 alpha=zero
427 DO m=1,lenpoly(i)
428 mm=m+1
429 IF (m==lenpoly(i)) mm=1
430 xx1=rpoly(4+3*(m-1)+1,i)
431 yy1=rpoly(4+3*(m-1)+2,i)
432 zz1=rpoly(4+3*(m-1)+3,i)
433 xx2=rpoly(4+3*(mm-1)+1,i)
434 yy2=rpoly(4+3*(mm-1)+2,i)
435 zz2=rpoly(4+3*(mm-1)+3,i)
436 vx1=xx1-x1
437 vy1=yy1-y1
438 vz1=zz1-z1
439 vx2=xx2-x1
440 vy2=yy2-y1
441 vz2=zz2-z1
442 nr1=sqrt(vx1**2+vy1**2+vz1**2)
443 nr2=sqrt(vx2**2+vy2**2+vz2**2)
444 IF(nr1 > zero) THEN
445 vx1=vx1/nr1
446 vy1=vy1/nr1
447 vz1=vz1/nr1
448 ELSE
449 cycle
450 ENDIF
451 IF(nr2 > zero) THEN
452 vx2=vx2/nr2
453 vy2=vy2/nr2
454 vz2=vz2/nr2
455 ELSE
456 cycle
457 ENDIF
458 ss=vx1*vx2+vy1*vy2+vz1*vz2
459 vvx=vy1*vz2-vz1*vy2
460 vvy=vz1*vx2-vx1*vz2
461 vvz=vx1*vy2-vy1*vx2
462 ss1=nx*vvx+ny*vvy+nz*vvz
463 IF(ss < -one) ss=-one
464 IF(ss > one) ss= one
465 IF (ss1>=zero) THEN
466 alpha=alpha+acos(ss)
467 ELSE
468 alpha=alpha-acos(ss)
469 ENDIF
470 ENDDO
471 IF (abs(alpha)<two*pi) ipout=1
472 ENDDO
473C
474 IF (ipout==1) cycle
475C---------------------------------------------
476C Le polygone j est un trou dans le polygone i
477C---------------------------------------------
478 ipoly(1,j)=-1
479 nhol=nhol+1
480 iadhol(nhol)=lenpoly(i)
481 DO k=1,lenpoly(j)
482 ipoly(6+iadhol(nhol)+k,i)=ipoly(6+k,j)
483 ielnod(iadhol(nhol)+k,i)=ielnod(k,j)
484 rpoly(4+3*iadhol(nhol)+3*(k-1)+1,i)=
485 . rpoly(4+3*(k-1)+1,j)
486 rpoly(4+3*iadhol(nhol)+3*(k-1)+2,i)=
487 . rpoly(4+3*(k-1)+2,j)
488 rpoly(4+3*iadhol(nhol)+3*(k-1)+3,i)=
489 . rpoly(4+3*(k-1)+3,j)
490 ENDDO
491 lenpoly(i)=lenpoly(i)+lenpoly(j)
492C Point interieur polygone j
493 vx1=rpoly(5,j)-rpoly(8,j)
494 vy1=rpoly(6,j)-rpoly(9,j)
495 vz1=rpoly(7,j)-rpoly(10,j)
496 ss=sqrt(vx1**2+vy1**2+vz1**2)
497 vx1=vx1/ss
498 vy1=vy1/ss
499 vz1=vz1/ss
500 vx2=ny*vz1-nz*vy1
501 vy2=nz*vx1-nx*vz1
502 vz2=nx*vy1-ny*vx1
503 x1=rpoly(5,j)
504 y1=rpoly(6,j)
505 z1=rpoly(7,j)
506 xloc(1,1)=zero
507 xloc(2,1)=zero
508 ylmin=ep30
509 ylmax=-ep30
510 DO k=2,lenpoly(j)
511 xx=rpoly(4+3*(k-1)+1,j)
512 yy=rpoly(4+3*(k-1)+2,j)
513 zz=rpoly(4+3*(k-1)+3,j)
514 vvx=xx-x1
515 vvy=yy-y1
516 vvz=zz-z1
517 xloc(1,k)=vvx*vx1+vvy*vy1+vvz*vz1
518 xloc(2,k)=vvx*vx2+vvy*vy2+vvz*vz2
519 IF (xloc(2,k)<ylmin) ylmin=xloc(2,k)
520 IF (xloc(2,k)>ylmax) ylmax=xloc(2,k)
521 ENDDO
522 ylsec=half*(ylmin+ylmax)
523C
524 nsec=0
525 DO k=1,lenpoly(j)
526 kk=k+1
527 IF (k==lenpoly(j)) kk=1
528 x1=xloc(1,k)
529 y1=xloc(2,k)
530 x2=xloc(1,kk)
531 y2=xloc(2,kk)
532 IF (y1-y2/=zero) THEN
533 alpha=(ylsec-y2)/(y1-y2)
534 IF (alpha>=zero.AND.alpha<=one) THEN
535 nsec=nsec+1
536 xsec(nsec)=alpha*x1+(one-alpha)*x2
537 ENDIF
538 ELSE
539 IF (y1==ylsec) THEN
540 nsec=nsec+1
541 xsec(nsec)=x1
542 nsec=nsec+1
543 xsec(nsec)=x2
544 ENDIF
545 ENDIF
546 ENDDO
547C
548 xsmin1=ep30
549 DO k=1,nsec
550 IF (xsec(k)<xsmin1) THEN
551 xsmin1=xsec(k)
552 ksmin=k
553 ENDIF
554 ENDDO
555 xsmin2=ep30
556 DO k=1,nsec
557 IF (k==ksmin) cycle
558 IF (xsec(k)<xsmin2) xsmin2=xsec(k)
559 ENDDO
560C
561 xs=half*(xsmin1+xsmin2)
562 ys=ylsec
563 phol(1,nhol)=rpoly(5,j)+xs*vx1+ys*vx2
564 phol(2,nhol)=rpoly(6,j)+xs*vy1+ys*vy2
565 phol(3,nhol)=rpoly(7,j)+xs*vz1+ys*vz2
566 ENDIF
567 ENDDO
568C
569 ipoly(2,i)=lenpoly(i)
570 ipoly(6+lenpoly(i)+1,i)=nhol
571 DO j=1,nhol
572 ipoly(6+lenpoly(i)+1+j,i)=iadhol(j)
573 rpoly(4+3*lenpoly(i)+3*(j-1)+1,i)=phol(1,j)
574 rpoly(4+3*lenpoly(i)+3*(j-1)+2,i)=phol(2,j)
575 rpoly(4+3*lenpoly(i)+3*(j-1)+3,i)=phol(3,j)
576 ENDDO
577 ENDDO
578C
579 RETURN
580 END
alpha
#define alpha
Definition eval.h:35
horiedge
subroutine horiedge(ipoly, rpoly, nx, ny, nz, nbnedge, inedge, rnedge, x0, y0, z0, inz, nns3, nref, aref, nnsp)
Definition horipoly.F:32
horipoly
subroutine horipoly(inedge, rnedge, ledge, nedge, ipoly, rpoly, iz, ielnod, npoly, nx, ny, nz, inz, ibric, nel, tagela)
Definition horipoly.F:113
max
#define max(a, b)
Definition macros.h:21