OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
horipoly.F File Reference
#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine horiedge (ipoly, rpoly, nx, ny, nz, nbnedge, inedge, rnedge, x0, y0, z0, inz, nns3, nref, aref, nnsp)
subroutine horipoly (inedge, rnedge, ledge, nedge, ipoly, rpoly, iz, ielnod, npoly, nx, ny, nz, inz, ibric, nel, tagela)

Function/Subroutine Documentation

◆ horiedge()

subroutine horiedge ( integer, dimension(*) ipoly,
rpoly,
nx,
ny,
nz,
integer nbnedge,
integer, dimension(6,*) inedge,
rnedge,
x0,
y0,
z0,
integer inz,
integer nns3,
integer, dimension(2,*) nref,
aref,
integer nnsp )

Definition at line 28 of file horipoly.F.

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
my_real
#define my_real
Definition cppsort.cpp:32

◆ horipoly()

subroutine horipoly ( integer, dimension(6,*) inedge,
rnedge,
integer, dimension(*) ledge,
integer nedge,
integer, dimension(6+2*nedge+1+nedge,*) ipoly,
rpoly,
integer, dimension(3,*) iz,
integer, dimension(nedge,*) ielnod,
integer npoly,
nx,
ny,
nz,
integer inz,
integer ibric,
integer nel,
integer, dimension(*) tagela )

Definition at line 109 of file horipoly.F.

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
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21