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

Go to the source code of this file.

Functions/Subroutines

subroutine facepoly (quad, tri, ntri, ipoly, rpoly, n, normf, normt, nsmax, nnp, nrpmax, nb, nv, lmin, npolmax, nppmax, info, ielnod, x, ifac, ilvout, mvid)

Function/Subroutine Documentation

◆ facepoly()

subroutine facepoly ( quad,
integer, dimension(nsmax,*) tri,
integer ntri,
integer, dimension(6+nppmax+1+npolmax,*) ipoly,
rpoly,
n,
normf,
normt,
integer nsmax,
integer nnp,
integer nrpmax,
integer nb,
integer nv,
lmin,
integer npolmax,
integer nppmax,
integer info,
integer, dimension(nppmax,*) ielnod,
x,
integer ifac,
integer ilvout,
integer mvid )

Definition at line 30 of file facepoly.F.

35C
36C-----------------------------------------------
37C I m p l i c i t T y p e s
38C-----------------------------------------------
39#include "implicit_f.inc"
40C-----------------------------------------------
41C C o m m o n B l o c k s
42C-----------------------------------------------
43#include "units_c.inc"
44C-----------------------------------------------
45C D u m m y A r g u m e n t s
46C-----------------------------------------------
47 INTEGER TRI(NSMAX,*), NTRI, NPPMAX, IPOLY(6+NPPMAX+1+NPOLMAX,*),
48 . NSMAX, NRPMAX, NNP, NB, NV, NPOLMAX, INFO,
49 . IELNOD(NPPMAX,*), IFAC, ILVOUT, MVID
50 my_real
51 . quad(3,*), rpoly(nrpmax+3*npolmax,*), n(*), normf(3,*),
52 . normt(3,*), lmin, x(3,*)
53C-----------------------------------------------
54C L o c a l V a r i a b l e s
55C-----------------------------------------------
56 INTEGER I, NINTER(4), NP, J, JJ, IDBL, NPI, NSEG, REDIR(5*NTRI),
57 . K, KK, JJ1, NNO, ID1, ID2, ISEG(5,5*NTRI), NSEGF, ISTOP,
58 . ISSEG, IPSEG, ITAGSEG(5*NTRI), ITAG(5*NTRI*2), ICLOSE,
59 . INEXT, NN, NPOLY, IADPOLY(NTRI), LENPOLY(NTRI), IAD, II,
60 . POLY(2,5*NTRI+1), I0, JJJ, REDIR_TMP(5*NTRI), NPI_TMP,
61 . ELSEG(5*NTRI,2), ELINTER(4,NTRI), ELNODE(5*NTRI*2),
62 . NNS, LENMAX, NEDGE, J1, J2, IEDGE, K1, K2,
63 . TEDGE(3*NTRI,3), EDGE(3*NTRI,3), N1, N2, IPOUT, M, MM,
64 . NHOL, IADHOL(NPOLMAX), ISEG3_OLD(5*NTRI), IDIFF, NSEC,
65 . KSMIN,NOK,ISSEGOLD
66 INTEGER JSEG(4), IP0SEG
67 INTEGER IADR(4), I1, I2, I3, I4, L, ICMAX, IMIN, IMAX, ALGO
68 INTEGER TEMP(5*NTRI),JTAG(2,5*NTRI)
69 my_real
70 . tole, tole2, a(3), x1, y1, z1, x2, y2, z2, ss1, ss2,
71 . x12, y12, z12, xa1, ya1, za1, alpha, xm, ym, zm,
72 . stmp(4,3), seg(5*ntri,3,2), nx, ny, nz, xa2, ya2, za2,
73 . xa12, ya12, za12, xa1p1, ya1p1, za1p1, beta,
74 . pinter(4,ntri,4), xl(ntri), xlref, xp1, yp1, zp1, xp2,
75 . yp2, zp2, node(3,5*ntri*2), xx1, yy1, zz1, dd1, dd2,
76 . xlc, ll, tolei, t1, t2, dd, xa1p2, ya1p2, za1p2,
77 . xedge(3*ntri,4), xx2, yy2, zz2, nr1, nr2, vx1, vy1, vz1,
78 . vx2, vy2, vz2, ss, vvx, vvy, vvz, xx, yy, zz,
79 . phol(3,npolmax), xloc(2,nppmax), xsec(nppmax), ylmin,
80 . ylmax, ylsec, xsmin1, xsmin2, xs, ys
81 my_real
82 . xp12, yp12, zp12, xa2p1, ya2p1, za2p1
83 my_real
84 . xa1a2, ya1a2, za1a2, alpha1, alpha2, norm
85
86C
87 INTEGER, ALLOCATABLE :: NODSEG(:,:)
88C
89 tole2=epsilon(zero)*0.5
90 tolei=tole2*ten
91 tole=tole2*lmin
92C
93C Initialization to avoid error in debug mode
94
95 IF(ntri>0) iseg(1:5,1:5*ntri) = 0
96
97 a(1)=quad(1,1)
98 a(2)=quad(2,1)
99 a(3)=quad(3,1)
100C Census of triangle edges
101 nedge=0
102 DO i=1,ntri
103 DO j=1,3
104 jj=j+1
105 IF (j==3) jj=1
106 j1=tri(i,j)
107 j2=tri(i,jj)
108 iedge=0
109 DO k=1,nedge
110 k1=edge(k,1)
111 k2=edge(k,2)
112 IF ((k1==j1.AND.k2==j2).OR.
113 . (k2==j1.AND.k1==j2)) iedge=k
114 ENDDO
115 IF (iedge==0) THEN
116 nedge=nedge+1
117 tedge(i,j)=nedge
118 edge(nedge,1)=j1
119 edge(nedge,2)=j2
120 ELSE
121 tedge(i,j)=iedge
122 ENDIF
123 ENDDO
124 ENDDO
125C Calculation of the intersections of the edges with the facet plane
126 DO i=1,nedge
127 n1=edge(i,1)
128 n2=edge(i,2)
129 x1=x(1,n1)
130 y1=x(2,n1)
131 z1=x(3,n1)
132 x2=x(1,n2)
133 y2=x(2,n2)
134 z2=x(3,n2)
135 ss1=(x1-a(1))*n(1)+(y1-a(2))*n(2)+(z1-a(3))*n(3)
136 ss2=(x2-a(1))*n(1)+(y2-a(2))*n(2)+(z2-a(3))*n(3)
137 edge(i,3)=0
138 IF (abs(ss1)<=tole.OR.abs(ss2)<=tole) THEN
139 IF (abs(ss1)<=tole.AND.abs(ss2)<=tole) cycle
140 ELSEIF (ss1<zero.AND.ss2<zero) THEN
141 cycle
142 ELSEIF (ss1>=zero.AND.ss2>=zero) THEN
143 cycle
144 ENDIF
145 x12=x2-x1
146 y12=y2-y1
147 z12=z2-z1
148 xa1=x1-a(1)
149 ya1=y1-a(2)
150 za1=z1-a(3)
151 alpha=-(xa1*n(1)+ya1*n(2)+za1*n(3))
152 . /(x12*n(1)+y12*n(2)+z12*n(3))
153 xm=x1+alpha*x12
154 ym=y1+alpha*y12
155 zm=z1+alpha*z12
156C
157 edge(i,3)=1
158 xedge(i,1)=xm
159 xedge(i,2)=ym
160 xedge(i,3)=zm
161 xedge(i,4)=ss1*ss2
162 ENDDO
163C Calculation of the intersections of triangles with the facet plane (segments)
164 DO i=1,4
165 ninter(i)=0
166 ENDDO
167 DO i=1,ntri
168 np=0
169 elseg(i,1)=i
170 elseg(i,2)=i
171 DO j=1,3
172 iedge=tedge(i,j)
173 IF (edge(iedge,3)==0) cycle
174 np=np+1
175 stmp(1,np)=xedge(iedge,1)
176 stmp(2,np)=xedge(iedge,2)
177 stmp(3,np)=xedge(iedge,3)
178 stmp(4,np)=xedge(iedge,4)
179 ENDDO
180 IF (np==0) THEN
181 seg(i,1,1)=zero
182 seg(i,2,1)=zero
183 seg(i,3,1)=zero
184 seg(i,1,2)=zero
185 seg(i,2,2)=zero
186 seg(i,3,2)=zero
187 cycle
188 ELSEIF (np==3) THEN
189 np=0
190 idbl=0
191 DO j=1,3
192 IF (stmp(4,j)<zero) THEN
193 np=np+1
194 seg(i,1,np)=stmp(1,j)
195 seg(i,2,np)=stmp(2,j)
196 seg(i,3,np)=stmp(3,j)
197 ELSEIF (stmp(4,j)==zero) THEN
198 IF (idbl==0) THEN
199 np=np+1
200 seg(i,1,np)=stmp(1,j)
201 seg(i,2,np)=stmp(2,j)
202 seg(i,3,np)=stmp(3,j)
203 idbl=1
204 ENDIF
205 ENDIF
206 ENDDO
207 ELSEIF (np==2) THEN
208 seg(i,1,1)=stmp(1,1)
209 seg(i,2,1)=stmp(2,1)
210 seg(i,3,1)=stmp(3,1)
211 seg(i,1,2)=stmp(1,2)
212 seg(i,2,2)=stmp(2,2)
213 seg(i,3,2)=stmp(3,2)
214 ENDIF
215C (Possible) intersection of each segment with the edges of the facet
216 DO j=1,4
217 npi=ninter(j)
218 IF (j/=4) THEN
219 jj=j+1
220 ELSE
221 jj=1
222 ENDIF
223 nx=normf(1,j)
224 ny=normf(2,j)
225 nz=normf(3,j)
226 x1=seg(i,1,1)
227 y1=seg(i,2,1)
228 z1=seg(i,3,1)
229 x2=seg(i,1,2)
230 y2=seg(i,2,2)
231 z2=seg(i,3,2)
232 x12=x2-x1
233 y12=y2-y1
234 z12=z2-z1
235 xa1=quad(1,j)
236 ya1=quad(2,j)
237 za1=quad(3,j)
238 xa2=quad(1,jj)
239 ya2=quad(2,jj)
240 za2=quad(3,jj)
241 xa12=xa2-xa1
242 ya12=ya2-ya1
243 za12=za2-za1
244 xa1p1=x1-xa1
245 ya1p1=y1-ya1
246 za1p1=z1-za1
247 xa1p2=x2-xa1
248 ya1p2=y2-ya1
249 za1p2=z2-za1
250C
251 ss1=xa1p1*nx+ya1p1*ny+za1p1*nz
252 ss2=xa1p2*nx+ya1p2*ny+za1p2*nz
253 IF (ss1>zero.AND.ss2>zero) THEN
254 seg(i,1,1)=zero
255 seg(i,2,1)=zero
256 seg(i,3,1)=zero
257 seg(i,1,2)=zero
258 seg(i,2,2)=zero
259 seg(i,3,2)=zero
260 cycle
261 ENDIF
262C
263 ss2=x12*nx+y12*ny+z12*nz
264 IF (ss2==zero) cycle
265 alpha=-ss1/ss2
266 IF (alpha<-tolei.OR.alpha>one+tolei) cycle
267 xm=x1+alpha*x12
268 ym=y1+alpha*y12
269 zm=z1+alpha*z12
270 beta=((xm-xa1)*xa12+(ym-ya1)*ya12+(zm-za1)*za12)
271 . /(xa12**2+ya12**2+za12**2)
272 IF (beta<-tolei.OR.beta>one+tolei) cycle
273 ss1=(x1-xa1)*nx+(y1-ya1)*ny+(z1-za1)*nz
274 ss2=(x2-xa1)*nx+(y2-ya1)*ny+(z2-za1)*nz
275C
276 IF (ss1*ss2>zero) THEN
277 IF (abs(ss1)<=abs(ss2)) THEN
278 seg(i,1,1)=xm
279 seg(i,2,1)=ym
280 seg(i,3,1)=zm
281 ELSE
282 seg(i,1,2)=xm
283 seg(i,2,2)=ym
284 seg(i,3,2)=zm
285 ENDIF
286 ELSEIF (ss1*ss2<zero) THEN
287 IF (ss1<zero) THEN
288 seg(i,1,2)=xm
289 seg(i,2,2)=ym
290 seg(i,3,2)=zm
291 ELSEIF (ss2<zero) THEN
292 seg(i,1,1)=xm
293 seg(i,2,1)=ym
294 seg(i,3,1)=zm
295 ENDIF
296 ELSE
297 IF (ss1==zero) THEN
298 IF (ss2<zero) THEN
299 seg(i,1,1)=xm
300 seg(i,2,1)=ym
301 seg(i,3,1)=zm
302 ELSEIF (ss2>=zero) THEN
303 seg(i,1,2)=xm
304 seg(i,2,2)=ym
305 seg(i,3,2)=zm
306 ENDIF
307 ELSEIF (ss2==zero) THEN
308 IF (ss1<zero) THEN
309 seg(i,1,2)=xm
310 seg(i,2,2)=ym
311 seg(i,3,2)=zm
312 ELSEIF (ss1>=zero) THEN
313 seg(i,1,1)=xm
314 seg(i,2,1)=ym
315 seg(i,3,1)=zm
316 ENDIF
317 ENDIF
318 ENDIF
319C
320 npi=npi+1
321 ninter(j)=npi
322 pinter(j,npi,1)=xm
323 pinter(j,npi,2)=ym
324 pinter(j,npi,3)=zm
325 elinter(j,npi)=i
326 IF(tri(i,5) == 3) THEN
327 pinter(j,npi,4)=zero
328 ELSE
329 nx=normt(1,i)
330 ny=normt(2,i)
331 nz=normt(3,i)
332 ss1=xa12*nx+ya12*ny+za12*nz
333 IF (ss1<=zero) THEN
334 pinter(j,npi,4)=one
335 ELSEIF (ss1>=zero) THEN
336 pinter(j,npi,4)=-one
337 ENDIF
338 ENDIF
339 ENDDO !J=1,4
340 ENDDO !I=1,NTRI
341C
342 nseg=ntri
343C
344C Possible intersection of internal segments with external segments
345 DO i=1,ntri
346 IF(tri(i,5)/=3) cycle
347 xa1=seg(i,1,1)
348 ya1=seg(i,2,1)
349 za1=seg(i,3,1)
350 xa2=seg(i,1,2)
351 ya2=seg(i,2,2)
352 za2=seg(i,3,2)
353 xa12=xa2-xa1
354 ya12=ya2-ya1
355 za12=za2-za1
356 DO j=1,ntri
357 IF(tri(j,5) == 3) cycle
358 xp1=seg(j,1,1)
359 yp1=seg(j,2,1)
360 zp1=seg(j,3,1)
361 xp2=seg(j,1,2)
362 yp2=seg(j,2,2)
363 zp2=seg(j,3,2)
364 xp12=xp2-xp1
365 yp12=yp2-yp1
366 zp12=zp2-zp1
367C
368 xa1p1=xp1-xa1
369 ya1p1=yp1-ya1
370 za1p1=zp1-za1
371 x12=ya12*za1p1-za12*ya1p1
372 y12=za12*xa1p1-xa12*za1p1
373 z12=xa12*ya1p1-ya12*xa1p1
374 ss1=x12*n(1)+y12*n(2)+z12*n(3)
375C
376 xa1p2=xp2-xa1
377 ya1p2=yp2-ya1
378 za1p2=zp2-za1
379 x12=ya12*za1p2-za12*ya1p2
380 y12=za12*xa1p2-xa12*za1p2
381 z12=xa12*ya1p2-ya12*xa1p2
382 ss2=x12*n(1)+y12*n(2)+z12*n(3)
383C
384 IF(ss1*ss2 > zero) cycle
385 x12=yp12*za1p1-zp12*ya1p1
386 y12=zp12*xa1p1-xp12*za1p1
387 z12=xp12*ya1p1-yp12*xa1p1
388 ss1=x12*n(1)+y12*n(2)+z12*n(3)
389C
390 xa2p1=xp1-xa2
391 ya2p1=yp1-ya2
392 za2p1=zp1-za2
393 x12=yp12*za2p1-zp12*ya2p1
394 y12=zp12*xa2p1-xp12*za2p1
395 z12=xp12*ya2p1-yp12*xa2p1
396 ss2=x12*n(1)+y12*n(2)+z12*n(3)
397 IF(ss1*ss2 > zero) cycle
398 IF(ss1 == zero .AND. ss2 /= zero) THEN
399 xm=xa1
400 ym=ya1
401 zm=za1
402 ELSEIF(ss2 == zero .AND. ss1 /= zero) THEN
403 xm=xa2
404 ym=ya2
405 zm=za2
406 ELSE
407C If SS1 = SS2 = 0 The segments are colinaries (unfortunate case)
408 cycle
409 ENDIF
410 seg(j,1,2)=xm
411 seg(j,2,2)=ym
412 seg(j,3,2)=zm
413 nseg=nseg+1
414 seg(nseg,1,1)=xm
415 seg(nseg,2,1)=ym
416 seg(nseg,3,1)=zm
417 elseg(nseg,1)=j
418 seg(nseg,1,2)=xp2
419 seg(nseg,2,2)=yp2
420 seg(nseg,3,2)=zp2
421 elseg(nseg,2)=j
422 ENDDO
423 ENDDO
424C
425C Additional segments on edges
426 DO i=1,4
427 IF (i/=4) THEN
428 ii=i+1
429 ELSE
430 ii=1
431 ENDIF
432 xa1=quad(1,i)
433 ya1=quad(2,i)
434 za1=quad(3,i)
435 xa2=quad(1,ii)
436 ya2=quad(2,ii)
437 za2=quad(3,ii)
438 xa12=xa2-xa1
439 ya12=ya2-ya1
440 za12=za2-za1
441 IF (ninter(i)==0) cycle
442C Sorting of points on the edges
443 npi=ninter(i)
444 DO j=1,npi
445 xm=pinter(i,j,1)
446 ym=pinter(i,j,2)
447 zm=pinter(i,j,3)
448 xl(j)=((xm-xa1)*xa12+(ym-ya1)*ya12+(zm-za1)*za12)
449 . /(xa12**2+ya12**2+za12**2)
450 ENDDO
451 DO j=1,npi
452 redir(j)=j
453 ENDDO
454 DO j=1,npi
455 xlref=xl(redir(j))
456 DO k=j+1,npi
457 xlc=xl(redir(k))
458 IF (xlc<=xlref) THEN
459 kk=redir(k)
460 redir(k)=redir(j)
461 redir(j)=kk
462 xlref=xlc
463 ENDIF
464 ENDDO
465 ENDDO
466C Identification of duplicates
467 DO j=1,npi
468 jj=redir(j)
469 IF (jj>0) THEN
470 x1=pinter(i,jj,1)
471 y1=pinter(i,jj,2)
472 z1=pinter(i,jj,3)
473 t1=pinter(i,jj,4)
474 DO k=j+1,npi
475 kk=redir(k)
476 IF (kk>0) THEN
477 x2=pinter(i,kk,1)
478 y2=pinter(i,kk,2)
479 z2=pinter(i,kk,3)
480 t2=pinter(i,kk,4)
481 dd=(x2-x1)**2+(y2-y1)**2+(z2-z1)**2+(t2-t1)**2
482 IF (dd<=tole) redir(k)=0
483 ENDIF
484 ENDDO
485 ENDIF
486 ENDDO
487C Elimination of zero-length segments
488 DO j=1,npi
489 jj=redir(j)
490 IF (jj>0) THEN
491 x1=pinter(i,jj,1)
492 y1=pinter(i,jj,2)
493 z1=pinter(i,jj,3)
494 DO k=j+1,npi
495 kk=redir(k)
496 IF (kk>0) THEN
497 x2=pinter(i,kk,1)
498 y2=pinter(i,kk,2)
499 z2=pinter(i,kk,3)
500 dd=(x2-x1)**2+(y2-y1)**2+(z2-z1)**2
501 IF (dd<=tole) THEN
502 redir(j)=0
503 redir(k)=0
504c print *,'on ne doit pas passer la non plus...',
505c + dd,tolem,i,j,k,JJ,KK
506c call flush(6)
507c stop
508 ENDIF
509 ENDIF
510 ENDDO
511 ENDIF
512 ENDDO
513C Condensation of the REDIR table
514 DO j=1,npi
515 redir_tmp(j)=redir(j)
516 ENDDO
517 npi_tmp=npi
518 npi=0
519 DO j=1,npi_tmp
520 IF (redir_tmp(j)>0) THEN
521 npi=npi+1
522 redir(npi)=redir_tmp(j)
523 ENDIF
524 ENDDO
525 IF (npi==0) cycle
526C Creation of segments on the edges
527 IF (pinter(i,redir(1),4)==-one) THEN
528 nseg=nseg+1
529 seg(nseg,1,1)=xa1
530 seg(nseg,2,1)=ya1
531 seg(nseg,3,1)=za1
532 elseg(nseg,1)=0
533 seg(nseg,1,2)=pinter(i,redir(1),1)
534 seg(nseg,2,2)=pinter(i,redir(1),2)
535 seg(nseg,3,2)=pinter(i,redir(1),3)
536 elseg(nseg,2)=elinter(i,redir(1))
537 ENDIF
538 DO j=1,npi-1
539 jj=redir(j)
540 jj1=redir(j+1)
541 IF (pinter(i,jj,4)==-one) cycle
542 xp1=pinter(i,jj,1)
543 yp1=pinter(i,jj,2)
544 zp1=pinter(i,jj,3)
545 xp2=pinter(i,jj1,1)
546 yp2=pinter(i,jj1,2)
547 zp2=pinter(i,jj1,3)
548 nseg=nseg+1
549 seg(nseg,1,1)=xp1
550 seg(nseg,2,1)=yp1
551 seg(nseg,3,1)=zp1
552 elseg(nseg,1)=elinter(i,jj)
553 seg(nseg,1,2)=xp2
554 seg(nseg,2,2)=yp2
555 seg(nseg,3,2)=zp2
556 elseg(nseg,2)=elinter(i,jj1)
557 ENDDO
558 IF (pinter(i,redir(npi),4)==one) THEN
559 nseg=nseg+1
560 seg(nseg,1,1)=pinter(i,redir(npi),1)
561 seg(nseg,2,1)=pinter(i,redir(npi),2)
562 seg(nseg,3,1)=pinter(i,redir(npi),3)
563 elseg(nseg,1)=elinter(i,redir(npi))
564 seg(nseg,1,2)=xa2
565 seg(nseg,2,2)=ya2
566 seg(nseg,3,2)=za2
567 elseg(nseg,2)=0
568 ENDIF
569 ENDDO
570C Census of nodes
571 nno=0
572 algo=0
573 DO i=1,nseg
574 x1=seg(i,1,1)
575 y1=seg(i,2,1)
576 z1=seg(i,3,1)
577 x2=seg(i,1,2)
578 y2=seg(i,2,2)
579 z2=seg(i,3,2)
580 IF (x1==zero.AND.y1==zero.AND.z1==zero.AND.
581 . x2==zero.AND.y2==zero.AND.z2==zero) THEN
582 iseg(1,i)=0
583 iseg(2,i)=0
584 iseg(3,i)=1
585 iseg(4,i)=0
586 iseg(5,i)=0
587 cycle
588 ENDIF
589 id1=0
590 id2=0
591 DO j=1,nno
592 xx1=node(1,j)
593 yy1=node(2,j)
594 zz1=node(3,j)
595 dd1=(xx1-x1)**2+(yy1-y1)**2+(zz1-z1)**2
596 dd2=(xx1-x2)**2+(yy1-y2)**2+(zz1-z2)**2
597 IF (dd1<=tole) id1=j
598 IF (dd2<=tole) id2=j
599 ENDDO
600 ll=(x2-x1)**2+(y2-y1)**2+(z2-z1)**2
601 IF (id1==0) THEN
602 nno=nno+1
603 node(1,nno)=x1
604 node(2,nno)=y1
605 node(3,nno)=z1
606 elnode(nno)=elseg(i,1)
607 id1=nno
608 ENDIF
609 IF (ll<=tole) id2=id1
610 IF (id2==0) THEN
611 nno=nno+1
612 node(1,nno)=x2
613 node(2,nno)=y2
614 node(3,nno)=z2
615 elnode(nno)=elseg(i,2)
616 id2=nno
617 ENDIF
618 iseg(1,i)=id1
619 iseg(2,i)=id2
620 iseg(3,i)=0
621 IF (id1==id2) iseg(3,i)=1
622 IF(elseg(i,1) == elseg(i,2)) THEN
623 IF(tri(elseg(i,1),5) == 3) THEN
624C Segment on internal triangle
625 iseg(4,i)=2
626 iseg(5,i)=1
627 algo=1
628 ELSE
629C Segment on airbag triangle or additional brick
630 iseg(4,i)=1
631 iseg(5,i)=0
632 ENDIF
633 ELSE
634C Segment on box edge
635 iseg(4,i)=1
636 iseg(5,i)=0
637 ENDIF
638 ENDDO
639C
640C Elimination of double segments
641 DO i=1,nseg
642 IF (iseg(3,i)==1) cycle
643 id1=iseg(1,i)
644 id2=iseg(2,i)
645 DO j=i+1,nseg
646 IF ((iseg(1,j)==id1.AND.iseg(2,j)==id2).OR.
647 . (iseg(1,j)==id2.AND.iseg(2,j)==id1)) iseg(3,j)=1
648 ENDDO
649 ENDDO
650C
651 nsegf=0
652 DO i=1,nseg
653 IF (iseg(3,i)==0) nsegf=nsegf+1
654 ENDDO
655 IF (nsegf<=1) RETURN
656C
657 IF(algo == 0) THEN
658C-----------------------------
659C ALGORITHME V12
660C Reconstruction of polygons
661C-----------------------------
662 istop=0
663 npoly=0
664 np=0
665 issegold=0
666C
667 DO i=1,nseg
668 iseg3_old(i)=iseg(3,i)
669 ENDDO
670C
671 DO WHILE (istop==0)
672 i=1
673 IF(i <= nseg) THEN
674 DO WHILE (iseg(3,i)==1.AND.i<=nseg)
675 i=i+1
676 ENDDO
677 ENDIF
678 IF (i==nseg+1) THEN
679 istop=1
680 cycle
681 ENDIF
682 isseg=i
683 ipseg=2
684 DO j=1,nno
685 itag(j)=0
686 ENDDO
687 DO j=1,nseg
688 itagseg(j)=0
689 ENDDO
690 i0=1
691 itag(iseg(1,isseg))=i0
692 itagseg(isseg)=i0
693 iclose=0
694 nok = 0
695 DO WHILE (iclose==0.AND.i<nseg)
696 i=i+1
697 inext=0
698 DO j=1,nseg
699 IF (j==isseg.OR.iseg(3,j)==1.OR.inext/=0) cycle
700 id1=iseg(1,j)
701 id2=iseg(2,j)
702 IF (id1==iseg(ipseg,isseg)) THEN
703 inext=1
704 isseg=j
705 ipseg=2
706 i0=i0+1
707 itag(iseg(1,isseg))=i0
708 itagseg(j)=i0
709 ELSEIF (id2==iseg(ipseg,isseg)) THEN
710 inext=1
711 isseg=j
712 ipseg=1
713 i0=i0+1
714 itag(iseg(2,isseg))=i0
715 itagseg(j)=-i0
716 ENDIF
717 ENDDO
718 IF (inext==0) THEN
719 iseg(3,isseg)=1
720 i=nseg
721 nok = 1
722 ELSE
723 nn=itag(iseg(ipseg,isseg))
724 IF (nn/=0) iclose=nn
725 ENDIF
726 ENDDO
727 IF (iclose/=0) THEN
728 npoly=npoly+1
729 iadpoly(npoly)=np
730 iad=np
731 DO j=1,nseg
732 jj=itagseg(j)
733 IF (abs(jj)>=iclose) THEN
734 np=np+1
735 poly(1,iad+abs(jj)-iclose+1)=j
736 IF (jj>0) THEN
737 poly(2,iad+abs(jj)-iclose+1)=2
738 ELSEIF (jj<0) THEN
739 poly(2,iad+abs(jj)-iclose+1)=1
740 ENDIF
741 iseg(3,j)=1
742 ENDIF
743 ENDDO
744 lenpoly(npoly)=np-iadpoly(npoly)
745 ELSEIF(nok==0.AND.issegold==isseg)THEN
746 iseg(3,isseg)=1
747c print *,'cycle corrige !!!',ISSEG,i,NSEG
748c stop 204
749 ENDIF
750 issegold=isseg
751C
752 idiff=0
753 DO j=1,nseg
754 idiff=max(idiff,abs(iseg3_old(j)-iseg(3,j)))
755 ENDDO
756 IF (idiff==0) THEN
757 WRITE(istdo,*)
758 WRITE(istdo,'(A25,I8,A25)')
759 . ' ** MONITORED VOLUME ID: ',mvid,' - INFINITE LOOP DETECTED'
760 WRITE(iout,'(A25,I8,A25)')
761 . ' ** MONITORED VOLUME ID: ',mvid,' - INFINITE LOOP DETECTED'
762 IF (ilvout >=1) WRITE(iout,'(A26,I8,A16,I8)')
763 . ' CUTTING BRICK NUMBER: ',nb,' - FACE NUMBER: ',ifac
764 CALL arret(2)
765 ENDIF
766 ENDDO ! DO WHILE (ISTOP==0)
767C
768 info=0
769 IF (npoly>npolmax) THEN
770 npolmax=npoly
771C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
772 IF(ilvout >=1) WRITE(iout,'(A,I10,A,I10)')
773 . ' MONITORED VOLUME ID: ',mvid,' NLAYER IS RESET TO ',npoly
774 info=1
775 ENDIF
776 lenmax=0
777 DO i=1,npoly
778 lenmax=max(lenmax,lenpoly(i))
779 ENDDO
780 IF (lenmax>nppmax) THEN
781 nppmax=lenmax
782 IF(ilvout >=1) WRITE(iout,'(A,I10,A,I10)')
783 . ' MONITORED VOLUME ID: ',mvid,' NPPMAX IS RESET TO ',lenmax
784 info=1
785 ENDIF
786 IF (info==1) RETURN
787 nnp=npoly
788C
789 ELSE
790C-------------------------
791C ALGORITME V14
792C-------------------------
793C Node-segment table
794C-------------------------
795 ALLOCATE (nodseg(nno,5))
796 DO i=1,nno
797 DO j=1,5
798 nodseg(i,j)=0
799 ENDDO
800 ENDDO
801 DO i=1,nseg
802 IF(iseg(3,i) == 1) cycle
803 DO j=1,2
804 k=iseg(j,i)
805 m=nodseg(k,1)
806 IF( m <= 3 ) THEN
807 nodseg(k,1)=m+1
808 nodseg(k,m+2)=i
809 ELSE
810 WRITE(iout,'(A,I10,2A,3E12.4,A1)')
811 . ' ** MONITORED VOLUME ID: ',mvid,' MORE THAN 4 SEGMENTS',
812 . ' CONNECTED TO NODE [',node(1,k),node(2,k),node(3,k),']'
813 CALL arret(2)
814 ENDIF
815 ENDDO
816 ENDDO
817C-------------------------------------------
818C Sorting of nodes
819C 4 segments: 2 internal + 2 external
820C 3 segments: 3 internal
821C 3 segments: 1 internal + 2 external
822C Testing and moving internal segments to the front
823C-------------------------------------------
824 DO i=1,nno
825 IF(nodseg(i,1) /= 4) cycle
826 DO k=1,4
827 jseg(k)=nodseg(i,k+1)
828 ENDDO
829 l=0
830 DO k=1,4
831 IF(iseg(5,jseg(k)) == 0) cycle
832 l=l+1
833 iadr(l)=k
834 ENDDO
835 IF(l /= 2) THEN
836 WRITE(iout,'(A,I10,A,3E12.4,2A,I2,A)')
837 . ' ** FVMBAG ID: ',mvid,' ERROR IN POLYGON BUILDING : NODE [',
838 . node(1,i),node(2,i),node(3,i),'] IS CONNECTED TO 4',
839 . ' EDGES',l,' BEING INTERNAL EDGES'
840 ENDIF
841 DO k=1,4
842 IF(iseg(5,jseg(k)) == 1) cycle
843 l=l+1
844 iadr(l)=k
845 ENDDO
846 DO k=1,4
847 nodseg(i,k+1)=jseg(iadr(k))
848 ENDDO
849 ENDDO
850C
851 DO i=1,nno
852 IF(nodseg(i,1) <= 1) THEN
853 WRITE(iout,'(A,I10,A,3E12.4,A,I2,A)')
854 . ' ** FVMBAG ID ',mvid,' ERROR IN POLYGON BUILDING : NODE [',
855 . node(1,i),node(2,i),node(3,i),'] IS CONNECTED TO ',
856 . nodseg(i,1),' EDGE'
857CFA CALL ARRET(2)
858 ENDIF
859 IF(nodseg(i,1) /= 2) cycle
860 l=0
861 IF(iseg(5,nodseg(i,2)) == 1) l=l+1
862 IF(iseg(5,nodseg(i,3)) == 1) l=l+1
863 IF(l == 1) THEN
864 WRITE(iout,'(A,I10,A,3E12.4,2A)')
865 . ' ** FVMBAG ID ',mvid,' ERROR IN POLYGON BUILDING : NODE [',
866 . node(1,i),node(2,i),node(3,i),'] IS CONNECTED TO 2',
867 . ' EDGES : 1 INTERNAL AND 1 EXTERNAL'
868CFA CALL ARRET(2)
869 ENDIF
870 ENDDO
871C
872 DO i=1,nno
873 IF(nodseg(i,1) /= 3) cycle
874 DO k=1,3
875 jseg(k)=nodseg(i,k+1)
876 ENDDO
877 l=0
878 DO k=1,3
879 IF(iseg(5,jseg(k)) == 0) cycle
880 l=l+1
881 iadr(l)=k
882 ENDDO
883 IF(l == 2) THEN
884 WRITE(iout,'(A,I10,A,3E12.4,2A)')
885 . ' ** FVMBAG ID: ',mvid,' ERROR IN POLYGON BUILDING : NODE [',
886 . node(1,i),node(2,i),node(3,i),'] IS CONNECTED TO 3',
887 . ' EDGES 2 BEING INTERNAL EDGES'
888CFA CALL ARRET(2)
889 ELSEIF(l == 1) THEN
890 DO k=1,3
891 IF(iseg(5,jseg(k)) == 1) cycle
892 l=l+1
893 iadr(l)=k
894 ENDDO
895 ELSEIF(l == 3) THEN
896 id1=iseg(1,jseg(1))
897 id2=iseg(2,jseg(1))
898 IF(id1==i) l=id2
899 IF(id2==i) l=id1
900 xa1=node(1,i)
901 ya1=node(2,i)
902 za1=node(3,i)
903 xa2=node(1,l)
904 ya2=node(2,l)
905 za2=node(3,l)
906 xa1a2=xa2-xa1
907 ya1a2=ya2-ya1
908 za1a2=za2-za1
909 norm=sqrt(xa1a2*xa1a2+ya1a2*ya1a2+za1a2*za1a2)
910 xa1a2=xa1a2/norm
911 ya1a2=ya1a2/norm
912 za1a2=za1a2/norm
913C
914 id1=iseg(1,jseg(2))
915 id2=iseg(2,jseg(2))
916 IF(id1==i) l=id2
917 IF(id2==i) l=id1
918 xp1=node(1,l)
919 yp1=node(2,l)
920 zp1=node(3,l)
921 xa1p1=xp1-xa1
922 ya1p1=yp1-ya1
923 za1p1=zp1-za1
924 norm=sqrt(xa1p1*xa1p1+ya1p1*ya1p1+za1p1*za1p1)
925 xa1p1=xa1p1/norm
926 ya1p1=ya1p1/norm
927 za1p1=za1p1/norm
928 ss=xa1a2*xa1p1+ya1a2*ya1p1+za1a2*za1p1
929 IF(ss < -one) ss=-one
930 IF(ss > one) ss= one
931 alpha1=acos(ss)
932C
933 id1=iseg(1,jseg(3))
934 id2=iseg(2,jseg(3))
935 IF(id1==i) l=id2
936 IF(id2==i) l=id1
937 xp1=node(1,l)
938 yp1=node(2,l)
939 zp1=node(3,l)
940 xa1p1=xp1-xa1
941 ya1p1=yp1-ya1
942 za1p1=zp1-za1
943 norm=sqrt(xa1p1*xa1p1+ya1p1*ya1p1+za1p1*za1p1)
944 xa1p1=xa1p1/norm
945 ya1p1=ya1p1/norm
946 za1p1=za1p1/norm
947 ss=xa1a2*xa1p1+ya1a2*ya1p1+za1a2*za1p1
948 IF(ss < -one) ss=-one
949 IF(ss > one) ss= one
950 alpha2=acos(ss)
951C
952 IF((alpha1 < zero .AND. alpha2 > zero).OR.
953 . (alpha2 < zero .AND. alpha1 > zero)) THEN
954 iadr(1)=1
955 iadr(2)=2
956 iadr(3)=3
957 ELSEIF(abs(alpha1) < abs(alpha2)) THEN
958 iadr(1)=2
959 iadr(2)=1
960 iadr(3)=3
961 ELSEIF(abs(alpha1) > abs(alpha2)) THEN
962 iadr(1)=3
963 iadr(2)=1
964 iadr(3)=2
965 ENDIF
966 ENDIF
967C
968 DO k=1,3
969 nodseg(i,k+1)=jseg(iadr(k))
970 ENDDO
971 ENDDO ! I=1,NNO
972C---------------------------------------------------------------------------
973C Construction of polygons (new algorithm V14)
974C JTAG(1,i) the internal polygon is on the side opposite to the segment's normal
975C JTAG(2,i) the internal polygon is on the side of the segment's normal
976C---------------------------------------------------------------------------
977 DO i=1,nseg
978 jtag(1,i)=0
979 jtag(2,i)=0
980 ENDDO
981 icmax=0
982C-----------------------------
983C Node connected to 4 segments
984C-----------------------------
985 DO i=1,nno
986 IF(nodseg(i,1) /= 4) cycle
987 i1=nodseg(i,2)
988 i2=nodseg(i,3)
989 i3=nodseg(i,4)
990 i4=nodseg(i,5)
991C Treatment of the 2 internal segments
992 IF(jtag(1,i1)==0.AND.jtag(1,i2)==0) THEN
993 icmax=icmax+1
994 jtag(1,i1)=icmax
995 jtag(1,i2)=icmax
996 ELSEIF(jtag(1,i1)/=0.AND.jtag(1,i2)==0) THEN
997 jtag(1,i2)=jtag(1,i1)
998 ELSEIF(jtag(1,i1)==0.AND.jtag(1,i2)/=0) THEN
999 jtag(1,i1)=jtag(1,i2)
1000 ELSEIF(jtag(1,i1) /= jtag(1,i2)) THEN
1001 imin=min(jtag(1,i1),jtag(1,i2))
1002 imax=max(jtag(1,i1),jtag(1,i2))
1003 DO j=1,nseg
1004 IF(jtag(1,j) == imax) jtag(1,j)=imin
1005 IF(jtag(2,j) == imax) jtag(2,j)=imin
1006 ENDDO
1007 ENDIF
1008C Treatment of 2 external segments
1009 id1=iseg(1,i3)
1010 id2=iseg(2,i3)
1011 IF(id1==i) l=id2
1012 IF(id2==i) l=id1
1013 xa1=node(1,i)
1014 ya1=node(2,i)
1015 za1=node(3,i)
1016 xa2=node(1,l)
1017 ya2=node(2,l)
1018 za2=node(3,l)
1019 xa1a2=xa2-xa1
1020 ya1a2=ya2-ya1
1021 za1a2=za2-za1
1022 norm=sqrt(xa1a2*xa1a2+ya1a2*ya1a2+za1a2*za1a2)
1023 xa1a2=xa1a2/norm
1024 ya1a2=ya1a2/norm
1025 za1a2=za1a2/norm
1026C
1027 id1=iseg(1,i1)
1028 id2=iseg(2,i1)
1029 IF(id1==i) l=id2
1030 IF(id2==i) l=id1
1031 xp1=node(1,l)
1032 yp1=node(2,l)
1033 zp1=node(3,l)
1034 xa1p1=xp1-xa1
1035 ya1p1=yp1-ya1
1036 za1p1=zp1-za1
1037 norm=sqrt(xa1p1*xa1p1+ya1p1*ya1p1+za1p1*za1p1)
1038 xa1p1=xa1p1/norm
1039 ya1p1=ya1p1/norm
1040 za1p1=za1p1/norm
1041 ss=xa1a2*xa1p1+ya1a2*ya1p1+za1a2*za1p1
1042 IF(ss < -one) ss=-one
1043 IF(ss > one) ss= one
1044 alpha1=acos(ss)
1045C
1046 id1=iseg(1,i2)
1047 id2=iseg(2,i2)
1048 IF(id1==i) l=id2
1049 IF(id2==i) l=id1
1050 xp1=node(1,l)
1051 yp1=node(2,l)
1052 zp1=node(3,l)
1053 xa1p1=xp1-xa1
1054 ya1p1=yp1-ya1
1055 za1p1=zp1-za1
1056 norm=sqrt(xa1p1*xa1p1+ya1p1*ya1p1+za1p1*za1p1)
1057 xa1p1=xa1p1/norm
1058 ya1p1=ya1p1/norm
1059 za1p1=za1p1/norm
1060 ss=xa1a2*xa1p1+ya1a2*ya1p1+za1a2*za1p1
1061 IF(ss < -one) ss=-one
1062 IF(ss > one) ss= one
1063 alpha2=acos(ss)
1064C
1065 jseg(1)=i3
1066 jseg(2)=i4
1067 IF(alpha1 < alpha2) THEN
1068 jseg(3)=i1
1069 jseg(4)=i2
1070 ELSE
1071 jseg(3)=i2
1072 jseg(4)=i1
1073 ENDIF
1074C
1075 DO m=1,2
1076 i1=jseg(m)
1077 i2=jseg(m+2)
1078 IF(jtag(1,i1)==0.AND.jtag(2,i2)==0) THEN
1079 icmax=icmax+1
1080 jtag(1,i1)=icmax
1081 jtag(2,i2)=icmax
1082 ELSEIF(jtag(1,i1)/=0.AND.jtag(2,i2)==0) THEN
1083 jtag(2,i2)=jtag(1,i1)
1084 ELSEIF(jtag(1,i1)==0.AND.jtag(2,i2)/=0) THEN
1085 jtag(1,i1)=jtag(2,i2)
1086 ELSEIF(jtag(1,i1) /= jtag(2,i2)) THEN
1087 imin=min(jtag(1,i1),jtag(2,i2))
1088 imax=max(jtag(1,i1),jtag(2,i2))
1089 DO j=1,nseg
1090 IF(jtag(1,j) == imax) jtag(1,j)=imin
1091 IF(jtag(2,j) == imax) jtag(2,j)=imin
1092 ENDDO
1093 ENDIF
1094 ENDDO
1095 ENDDO
1096C------------------------------------------------------
1097C Knot connects to 3 segments: 1 internal and 2 external
1098C------------------------------------------------------
1099 DO i=1,nno
1100 IF(nodseg(i,1) /= 3) cycle
1101 jseg(1)=nodseg(i,2)
1102 jseg(2)=nodseg(i,3)
1103 IF(iseg(5,jseg(2)) == 1) cycle
1104 jseg(3)=nodseg(i,4)
1105 IF(iseg(5,jseg(3)) == 1) cycle
1106C
1107C Normal in the internal segment
1108 i1=jseg(1)
1109 nn=elseg(i1,1)
1110 nx=normt(1,nn)
1111 ny=normt(2,nn)
1112 nz=normt(3,nn)
1113 xx=ny*n(3)-nz*n(2)
1114 yy=nz*n(1)-nx*n(3)
1115 zz=nx*n(2)-ny*n(1)
1116 nx=n(2)*zz-n(3)*yy
1117 ny=n(3)*xx-n(1)*zz
1118 nz=n(1)*yy-n(2)*xx
1119C
1120 id1=iseg(1,i1)
1121 id2=iseg(2,i1)
1122 IF(id1==i) l=id2
1123 IF(id2==i) l=id1
1124 xa1=node(1,l)
1125 ya1=node(2,l)
1126 za1=node(3,l)
1127C
1128 i2=jseg(2)
1129 id1=iseg(1,i2)
1130 id2=iseg(2,i2)
1131 IF(id1==i) l=id2
1132 IF(id2==i) l=id1
1133 xp1=node(1,l)
1134 yp1=node(2,l)
1135 zp1=node(3,l)
1136 xa1p1=xp1-xa1
1137 ya1p1=yp1-ya1
1138 za1p1=zp1-za1
1139 norm=sqrt(xa1p1*xa1p1+ya1p1*ya1p1+za1p1*za1p1)
1140 xa1p1=xa1p1/norm
1141 ya1p1=ya1p1/norm
1142 za1p1=za1p1/norm
1143 ss=nx*xa1p1+ny*ya1p1+nz*za1p1
1144 IF(ss < -one) ss=-one
1145 IF(ss > one) ss= one
1146 alpha1=acos(ss)
1147C
1148 i3=jseg(3)
1149 id1=iseg(1,i3)
1150 id2=iseg(2,i3)
1151 IF(id1==i) l=id2
1152 IF(id2==i) l=id1
1153 xp1=node(1,l)
1154 yp1=node(2,l)
1155 zp1=node(3,l)
1156 xa1p1=xp1-xa1
1157 ya1p1=yp1-ya1
1158 za1p1=zp1-za1
1159 norm=sqrt(xa1p1*xa1p1+ya1p1*ya1p1+za1p1*za1p1)
1160 xa1p1=xa1p1/norm
1161 ya1p1=ya1p1/norm
1162 za1p1=za1p1/norm
1163 ss=nx*xa1p1+ny*ya1p1+nz*za1p1
1164 IF(ss < -one) ss=-one
1165 IF(ss > one) ss= one
1166 alpha2=acos(ss)
1167C
1168C The unlikely case alpha1=alpha2 needs to be completed
1169C
1170 IF(alpha1 < alpha2) THEN
1171 iadr(1)=2
1172 iadr(2)=1
1173 ELSE
1174 iadr(1)=1
1175 iadr(2)=2
1176 ENDIF
1177C
1178 DO m=1,2
1179 i1=jseg(1)
1180 i2=jseg(m+1)
1181 l =iadr(m)
1182 IF(jtag(l,i1)==0.AND.jtag(1,i2)==0) THEN
1183 icmax=icmax+1
1184 jtag(l,i1)=icmax
1185 jtag(1,i2)=icmax
1186 ELSEIF(jtag(l,i1)/=0.AND.jtag(1,i2)==0) THEN
1187 jtag(1,i2)=jtag(l,i1)
1188 ELSEIF(jtag(l,i1)==0.AND.jtag(1,i2)/=0) THEN
1189 jtag(l,i1)=jtag(1,i2)
1190 ELSEIF(jtag(l,i1) /= jtag(1,i2)) THEN
1191 imin=min(jtag(l,i1),jtag(1,i2))
1192 imax=max(jtag(l,i1),jtag(1,i2))
1193 DO j=1,nseg
1194 IF(jtag(1,j) == imax) jtag(1,j)=imin
1195 IF(jtag(2,j) == imax) jtag(2,j)=imin
1196 ENDDO
1197 ENDIF
1198 ENDDO
1199 ENDDO
1200C-------------------------------------
1201C Knot connects to 3 internal segments
1202C-------------------------------------
1203 DO i=1,nno
1204 IF(nodseg(i,1) /= 3) cycle
1205 jseg(1)=nodseg(i,2)
1206 IF(iseg(5,jseg(1)) == 0) cycle
1207 jseg(2)=nodseg(i,3)
1208 IF(iseg(5,jseg(2)) == 0) cycle
1209 jseg(3)=nodseg(i,4)
1210 IF(iseg(5,jseg(3)) == 0) cycle
1211C
1212C Normal to the 1st internal segment
1213 i1=jseg(1)
1214 nn=elseg(i1,1)
1215 nx=normt(1,nn)
1216 ny=normt(2,nn)
1217 nz=normt(3,nn)
1218 xx=ny*n(3)-nz*n(2)
1219 yy=nz*n(1)-nx*n(3)
1220 zz=nx*n(2)-ny*n(1)
1221 nx=n(2)*zz-n(3)*yy
1222 ny=n(3)*xx-n(1)*zz
1223 nz=n(1)*yy-n(2)*xx
1224C
1225 id1=iseg(1,i1)
1226 id2=iseg(2,i1)
1227 IF(id1==i) l=id2
1228 IF(id2==i) l=id1
1229 xa1=node(1,l)
1230 ya1=node(2,l)
1231 za1=node(3,l)
1232C
1233 i2=jseg(2)
1234 id1=iseg(1,i2)
1235 id2=iseg(2,i2)
1236 IF(id1==i) l=id2
1237 IF(id2==i) l=id1
1238 xp1=node(1,l)
1239 yp1=node(2,l)
1240 zp1=node(3,l)
1241 xa1p1=xp1-xa1
1242 ya1p1=yp1-ya1
1243 za1p1=zp1-za1
1244 norm=sqrt(xa1p1*xa1p1+ya1p1*ya1p1+za1p1*za1p1)
1245 xa1p1=xa1p1/norm
1246 ya1p1=ya1p1/norm
1247 za1p1=za1p1/norm
1248 ss=nx*xa1p1+ny*ya1p1+nz*za1p1
1249 IF(ss < -one) ss=-one
1250 IF(ss > one) ss= one
1251 alpha=acos(ss)
1252 IF(abs(alpha) < half*pi) THEN
1253 iadr(1)=1
1254 iadr(2)=2
1255C IADR(1)=2
1256C IADR(2)=1
1257 ELSE
1258 iadr(1)=2
1259 iadr(2)=1
1260C IADR(1)=1
1261C IADR(2)=2
1262 ENDIF
1263C
1264 DO m=1,2
1265 i1=jseg(1)
1266 i2=jseg(m+1)
1267 l =iadr(m)
1268 IF(jtag(l,i1)==0.AND.jtag(2,i2)==0) THEN
1269 icmax=icmax+1
1270 jtag(l,i1)=icmax
1271 jtag(2,i2)=icmax
1272 ELSEIF(jtag(l,i1)/=0.AND.jtag(2,i2)==0) THEN
1273 jtag(2,i2)=jtag(l,i1)
1274 ELSEIF(jtag(l,i1)==0.AND.jtag(2,i2)/=0) THEN
1275 jtag(l,i1)=jtag(2,i2)
1276 ELSEIF(jtag(l,i1) /= jtag(2,i2)) THEN
1277 imin=min(jtag(l,i1),jtag(2,i2))
1278 imax=max(jtag(l,i1),jtag(2,i2))
1279 DO j=1,nseg
1280 IF(jtag(1,j) == imax) jtag(1,j)=imin
1281 IF(jtag(2,j) == imax) jtag(2,j)=imin
1282 ENDDO
1283 ENDIF
1284 ENDDO
1285C
1286 i1=jseg(2)
1287 i2=jseg(3)
1288 IF(jtag(1,i1)==0.AND.jtag(1,i2)==0) THEN
1289 icmax=icmax+1
1290 jtag(1,i1)=icmax
1291 jtag(1,i2)=icmax
1292 ELSEIF(jtag(1,i1)/=0.AND.jtag(1,i2)==0) THEN
1293 jtag(1,i2)=jtag(1,i1)
1294 ELSEIF(jtag(1,i1)==0.AND.jtag(1,i2)/=0) THEN
1295 jtag(1,i1)=jtag(1,i2)
1296 ELSEIF(jtag(1,i1) /= jtag(1,i2)) THEN
1297 imin=min(jtag(1,i1),jtag(1,i2))
1298 imax=max(jtag(1,i1),jtag(1,i2))
1299 DO j=1,nseg
1300 IF(jtag(1,j) == imax) jtag(1,j)=imin
1301 IF(jtag(2,j) == imax) jtag(2,j)=imin
1302 ENDDO
1303 ENDIF
1304C
1305 ENDDO
1306C-------------------------------------
1307C Knot connects to 2 internal segments
1308C-------------------------------------
1309 DO i=1,nno
1310 IF(nodseg(i,1) /= 2) cycle
1311 i1=nodseg(i,2)
1312 i2=nodseg(i,3)
1313 IF(iseg(5,i1) == 0) cycle
1314C
1315 IF(jtag(2,i1)==0.AND.jtag(2,i2)==0) THEN
1316 icmax=icmax+1
1317 jtag(2,i1)=icmax
1318 jtag(2,i2)=icmax
1319 ELSEIF(jtag(2,i1)/=0.AND.jtag(2,i2)==0) THEN
1320 jtag(2,i2)=jtag(2,i1)
1321 ELSEIF(jtag(2,i1)==0.AND.jtag(2,i2)/=0) THEN
1322 jtag(2,i1)=jtag(2,i2)
1323 ELSEIF(jtag(2,i1) /= jtag(2,i2)) THEN
1324 imin=min(jtag(2,i1),jtag(2,i2))
1325 imax=max(jtag(2,i1),jtag(2,i2))
1326 DO j=1,nseg
1327 IF(jtag(1,j) == imax) jtag(1,j)=imin
1328 IF(jtag(2,j) == imax) jtag(2,j)=imin
1329 ENDDO
1330 ENDIF
1331C
1332 IF(jtag(1,i1)==0.AND.jtag(1,i2)==0) THEN
1333 icmax=icmax+1
1334 jtag(1,i1)=icmax
1335 jtag(1,i2)=icmax
1336 ELSEIF(jtag(1,i1)/=0.AND.jtag(1,i2)==0) THEN
1337 jtag(1,i2)=jtag(1,i1)
1338 ELSEIF(jtag(1,i1)==0.AND.jtag(1,i2)/=0) THEN
1339 jtag(1,i1)=jtag(1,i2)
1340 ELSEIF(jtag(1,i1) /= jtag(1,i2)) THEN
1341 imin=min(jtag(1,i1),jtag(1,i2))
1342 imax=max(jtag(1,i1),jtag(1,i2))
1343 DO j=1,nseg
1344 IF(jtag(1,j) == imax) jtag(1,j)=imin
1345 IF(jtag(2,j) == imax) jtag(2,j)=imin
1346 ENDDO
1347 ENDIF
1348 ENDDO
1349C-------------------------------------
1350C Knot connects to 2 external segments
1351C-------------------------------------
1352 DO i=1,nno
1353 IF(nodseg(i,1) /= 2) cycle
1354 i1=nodseg(i,2)
1355 i2=nodseg(i,3)
1356 IF(iseg(5,i1) == 1) cycle
1357C
1358 IF(jtag(1,i1)==0.AND.jtag(1,i2)==0) THEN
1359 icmax=icmax+1
1360 jtag(1,i1)=icmax
1361 jtag(1,i2)=icmax
1362 ELSEIF(jtag(1,i1)/=0.AND.jtag(1,i2)==0) THEN
1363 jtag(1,i2)=jtag(1,i1)
1364 ELSEIF(jtag(1,i1)==0.AND.jtag(1,i2)/=0) THEN
1365 jtag(1,i1)=jtag(1,i2)
1366 ELSEIF(jtag(1,i1) /= jtag(1,i2)) THEN
1367 imin=min(jtag(1,i1),jtag(1,i2))
1368 imax=max(jtag(1,i1),jtag(1,i2))
1369 DO j=1,nseg
1370 IF(jtag(1,j) == imax) jtag(1,j)=imin
1371 ENDDO
1372 ENDIF
1373 ENDDO
1374C
1375 DEALLOCATE (nodseg)
1376C
1377 info=0
1378 npoly=0
1379 DO i=1,icmax
1380 ii=0
1381 DO j=1,nseg
1382 IF (jtag(1,j) == i .OR. jtag(2,j) == i) ii=ii+1
1383 ENDDO
1384 IF (ii/=0) THEN
1385 npoly=npoly+1
1386 lenpoly(npoly)=ii
1387 ENDIF
1388 ENDDO
1389 IF (npoly>npolmax) THEN
1390 npolmax=npoly
1391 IF(ilvout >=1) WRITE(iout,'(A,I10,A,I10)')
1392 . ' MONITORED VOLUME ID: ',mvid,' NLAYER IS RESET TO ',npoly
1393 info=1
1394 ENDIF
1395 lenmax=0
1396 DO i=1,npoly
1397 lenmax=max(lenmax,lenpoly(i))
1398 ENDDO
1399 IF (lenmax>nppmax) THEN
1400 nppmax=lenmax
1401 IF(ilvout >=1) WRITE(iout,'(A,I10,A,I10)')
1402 . ' MONITORED VOLUME ID: ',mvid,' NPPMAX IS RESET TO ',lenmax
1403 info=1
1404 ENDIF
1405 IF (info==1) RETURN
1406C
1407 nnp=npoly
1408 npoly=0
1409 jj=0
1410 DO i=1,icmax
1411 ii=0
1412 DO j=1,nseg
1413 IF (jtag(1,j) == i .OR. jtag(2,j) == i) ii=ii+1
1414 ENDDO
1415 IF (ii/=0) THEN
1416 npoly=npoly+1
1417 lenpoly(npoly)=ii
1418 DO j=1,nseg
1419 IF (jtag(1,j) == i .OR. jtag(2,j) == i) THEN
1420 jj=jj+1
1421 poly(1,jj)=j
1422 ENDIF
1423 ENDDO
1424 ENDIF
1425 ENDDO
1426C----------------
1427C Order polygons
1428C----------------
1429 iad=0
1430 DO i=1,npoly
1431 IF(lenpoly(i) < 3) THEN
1432 WRITE(iout,'(A,I8,3(A,I5),A)')
1433 . ' CUTTING BRICK NUMBER: ',nb,' face number: ',IFAC,
1434 . ' polygone',I,' has only',LENPOLY(I),' edges'
1435 ENDIF
1436 DO J=1,LENPOLY(I)
1437 TEMP(J)=POLY(1,IAD+J)
1438 ITAGSEG(J)=0
1439 ENDDO
1440 ISSEG=TEMP(1)
1441 POLY(2,IAD+1)=1
1442 ITAGSEG(1)=1
1443 IP0SEG=2
1444 I0=ISEG(1,ISSEG)
1445C J=0
1446C ICLOSE=0
1447 DO J=1,LENPOLY(I)-1
1448C J=J+1
1449 I1=ISEG(IP0SEG,ISSEG)
1450 K=0
1451 ISTOP=0
1452 DO WHILE (ISTOP == 0)
1453 K=K+1
1454 IF( K > NSEG ) THEN
1455 WRITE(IOUT,'(a,i8,2(a,i5),a)')
1456 . ' cutting brick number: ',NB,' face number: ',IFAC,
1457 . ' polygone',I,' is not closed'
1458 CALL ARRET(2)
1459 ENDIF
1460C
1461.OR. IF(TEMP(K) == ISSEG ITAGSEG(K) == 1) CYCLE
1462 K1=ISEG(1,TEMP(K))
1463 K2=ISEG(2,TEMP(K))
1464 IF(K1 == I1) THEN
1465 IP0SEG=2
1466 ISSEG=TEMP(K)
1467 POLY(1,IAD+J+1)=ISSEG
1468 POLY(2,IAD+J+1)=1
1469 ISTOP=1
1470 ITAGSEG(K)=1
1471 ELSEIF(K2 == I1) THEN
1472 IP0SEG=1
1473 ISSEG=TEMP(K)
1474 POLY(1,IAD+J+1)=ISSEG
1475 POLY(2,IAD+J+1)=2
1476 ISTOP=1
1477 ITAGSEG(K)=1
1478 ENDIF
1479 ENDDO
1480 ENDDO
1481 IF(ISEG(IP0SEG,ISSEG) /= I0) THEN
1482 WRITE(IOUT,'(a,i8,2(a,i5),a)')
1483 . ' cutting brick number: ',NB,' face number: ',IFAC,
1484 . ' polygone',I,' is not closed'
1485 ENDIF
1486 IAD=IAD+LENPOLY(I)
1487 ENDDO
1488C
1489 ENDIF ! ALGO
1490C
1491 NNS=0
1492 IAD=0
1493 DO I=1,NPOLY
1494 IPOLY(1,I)=2
1495 IPOLY(2,I)=LENPOLY(I)
1496 IPOLY(3,I)=NB
1497 IPOLY(4,I)=NV
1498 IPOLY(5,I)=0
1499 IPOLY(6,I)=0
1500 RPOLY(1,I)=ZERO
1501 RPOLY(2,I)=N(1)
1502 RPOLY(3,I)=N(2)
1503 RPOLY(4,I)=N(3)
1504 DO J=1,LENPOLY(I)
1505 JJ =POLY(1,IAD+J)
1506 JJJ=POLY(2,IAD+J)
1507 ID1=ISEG(JJJ,JJ)
1508 NNS=NNS+1
1509 IPOLY(6+J,I)=NNS
1510 IELNOD(J,I)=ELNODE(ID1)
1511 RPOLY(4+3*(J-1)+1,I)=NODE(1,ID1)
1512 RPOLY(4+3*(J-1)+2,I)=NODE(2,ID1)
1513 RPOLY(4+3*(J-1)+3,I)=NODE(3,ID1)
1514 ENDDO
1515 IAD=IAD+LENPOLY(I)
1516 ENDDO
1517C
1518C Searching for holes
1519 INFO=0
1520 DO I=1,NPOLY
1521C
1522 NHOL=0
1523 DO J=1,NPOLY
1524 IADHOL(J)=0
1525 ENDDO
1526C
1527 DO J=1,NPOLY
1528 IF (I==J) CYCLE
1529 ALPHA=ZERO
1530 X1=RPOLY(5,J)
1531 Y1=RPOLY(6,J)
1532 Z1=RPOLY(7,J)
1533 DO K=1,LENPOLY(I)
1534 KK=K+1
1535 IF (K==LENPOLY(I)) KK=1
1536 XX1=RPOLY(4+3*(K-1)+1,I)
1537 YY1=RPOLY(4+3*(K-1)+2,I)
1538 ZZ1=RPOLY(4+3*(K-1)+3,I)
1539 XX2=RPOLY(4+3*(KK-1)+1,I)
1540 YY2=RPOLY(4+3*(KK-1)+2,I)
1541 ZZ2=RPOLY(4+3*(KK-1)+3,I)
1542 VX1=XX1-X1
1543 VY1=YY1-Y1
1544 VZ1=ZZ1-Z1
1545 VX2=XX2-X1
1546 VY2=YY2-Y1
1547 VZ2=ZZ2-Z1
1548 NR1=SQRT(VX1**2+VY1**2+VZ1**2)
1549 NR2=SQRT(VX2**2+VY2**2+VZ2**2)
1550 IF(NR1 > ZERO) THEN
1551 VX1=VX1/NR1
1552 VY1=VY1/NR1
1553 VZ1=VZ1/NR1
1554 ENDIF
1555 IF(NR2 > ZERO) THEN
1556 VX2=VX2/NR2
1557 VY2=VY2/NR2
1558 VZ2=VZ2/NR2
1559 ENDIF
1560 SS=VX1*VX2+VY1*VY2+VZ1*VZ2
1561 IF(SS < -ONE) SS=-ONE
1562 IF(SS > ONE) SS= ONE
1563 VVX=VY1*VZ2-VZ1*VY2
1564 VVY=VZ1*VX2-VX1*VZ2
1565 VVZ=VX1*VY2-VY1*VX2
1566 SS1=N(1)*VVX+N(2)*VVY+N(3)*VVZ
1567 IF (SS1>ZERO) THEN
1568 ALPHA=ALPHA+ACOS(SS)
1569 ELSEIF(SS1<ZERO) THEN
1570 ALPHA=ALPHA-ACOS(SS)
1571 ENDIF
1572 ENDDO
1573C
1574 IF (ABS(ALPHA)<TWO*PI) CYCLE
1575C IF (ABS(ALPHA)>=ONE) THEN
1576C The first point of polygon i is inside polygon j
1577C A point is internal if ABS(ALPHA)=2PI!
1578C We test all the others
1579 IPOUT=0
1580 DO K=2,LENPOLY(J)
1581 X1=RPOLY(4+3*(K-1)+1,J)
1582 Y1=RPOLY(4+3*(K-1)+2,J)
1583 Z1=RPOLY(4+3*(K-1)+3,J)
1584 ALPHA=ZERO
1585 DO M=1,LENPOLY(I)
1586 MM=M+1
1587 IF (M==LENPOLY(I)) MM=1
1588 XX1=RPOLY(4+3*(M-1)+1,I)
1589 YY1=RPOLY(4+3*(M-1)+2,I)
1590 ZZ1=RPOLY(4+3*(M-1)+3,I)
1591 XX2=RPOLY(4+3*(MM-1)+1,I)
1592 YY2=RPOLY(4+3*(MM-1)+2,I)
1593 ZZ2=RPOLY(4+3*(MM-1)+3,I)
1594 VX1=XX1-X1
1595 VY1=YY1-Y1
1596 VZ1=ZZ1-Z1
1597 VX2=XX2-X1
1598 VY2=YY2-Y1
1599 VZ2=ZZ2-Z1
1600 NR1=SQRT(VX1**2+VY1**2+VZ1**2)
1601 NR2=SQRT(VX2**2+VY2**2+VZ2**2)
1602 IF(NR1 > ZERO) THEN
1603 VX1=VX1/NR1
1604 VY1=VY1/NR1
1605 VZ1=VZ1/NR1
1606 ENDIF
1607 IF(NR2 > ZERO) THEN
1608 VX2=VX2/NR2
1609 VY2=VY2/NR2
1610 VZ2=VZ2/NR2
1611 ENDIF
1612 SS=VX1*VX2+VY1*VY2+VZ1*VZ2
1613 IF(SS < -ONE) SS=-ONE
1614 IF(SS > ONE) SS= ONE
1615 VVX=VY1*VZ2-VZ1*VY2
1616 VVY=VZ1*VX2-VX1*VZ2
1617 VVZ=VX1*VY2-VY1*VX2
1618 SS1=N(1)*VVX+N(2)*VVY+N(3)*VVZ
1619 IF (SS1>ZERO) THEN
1620 ALPHA=ALPHA+ACOS(SS)
1621 ELSEIF(SS1<ZERO) THEN
1622 ALPHA=ALPHA-ACOS(SS)
1623 ENDIF
1624 ENDDO
1625C IF (ABS(ALPHA)<ONE) IPOUT=1
1626 IF (ABS(ALPHA)<TWO*PI) IPOUT=1
1627 ENDDO
1628C
1629 IF (IPOUT==1) CYCLE
1630C Polygon j is a hole in polygon i
1631 IPOLY(1,J)=-1
1632 NHOL=NHOL+1
1633 IADHOL(NHOL)=LENPOLY(I)
1634 IF (LENPOLY(I)+LENPOLY(J)>NPPMAX) THEN
1635 NPPMAX=LENPOLY(I)+LENPOLY(J)
1636 IF(ILVOUT >=1) WRITE(IOUT,'(a,i10,a,i10)')
1637 . ' monitored volume id: ',MVID,
1638 . ' nppmax is reset to ',NPPMAX
1639 INFO=1
1640 RETURN
1641 ELSE
1642 IF (ILVOUT >=1) THEN
1643 WRITE(IOUT,'(/a25,i10,a25)')
1644 . ' ** monitored volume id: ',MVID,' - hole detected'
1645 WRITE(IOUT,'(a26,i8,a16,i8)')
1646 . ' cutting brick number: ',NB,' face number: ',IFAC
1647 ENDIF
1648 DO K=1,LENPOLY(J)
1649 IPOLY(6+IADHOL(NHOL)+K,I)=IPOLY(6+K,J)
1650 IELNOD(IADHOL(NHOL)+K,I)=IELNOD(K,J)
1651 RPOLY(4+3*IADHOL(NHOL)+3*(K-1)+1,I)=
1652 . RPOLY(4+3*(K-1)+1,J)
1653 RPOLY(4+3*IADHOL(NHOL)+3*(K-1)+2,I)=
1654 . RPOLY(4+3*(K-1)+2,J)
1655 RPOLY(4+3*IADHOL(NHOL)+3*(K-1)+3,I)=
1656 . RPOLY(4+3*(K-1)+3,J)
1657 ENDDO
1658 ENDIF
1659 LENPOLY(I)=LENPOLY(I)+LENPOLY(J)
1660C Interior polygon j.
1661 VX1=QUAD(1,2)-QUAD(1,1)
1662 VY1=QUAD(2,2)-QUAD(2,1)
1663 VZ1=QUAD(3,2)-QUAD(3,1)
1664 SS=SQRT(VX1**2+VY1**2+VZ1**2)
1665 VX1=VX1/SS
1666 VY1=VY1/SS
1667 VZ1=VZ1/SS
1668 VX2=N(2)*VZ1-N(3)*VY1
1669 VY2=N(3)*VX1-N(1)*VZ1
1670 VZ2=N(1)*VY1-N(2)*VX1
1671 X1=RPOLY(5,J)
1672 Y1=RPOLY(6,J)
1673 Z1=RPOLY(7,J)
1674 XLOC(1,1)=ZERO
1675 XLOC(2,1)=ZERO
1676 YLMIN=EP30
1677 YLMAX=-EP30
1678 DO K=2,LENPOLY(J)
1679 XX=RPOLY(4+3*(K-1)+1,J)
1680 YY=RPOLY(4+3*(K-1)+2,J)
1681 ZZ=RPOLY(4+3*(K-1)+3,J)
1682 VVX=XX-X1
1683 VVY=YY-Y1
1684 VVZ=ZZ-Z1
1685 XLOC(1,K)=VVX*VX1+VVY*VY1+VVZ*VZ1
1686 XLOC(2,K)=VVX*VX2+VVY*VY2+VVZ*VZ2
1687 IF (XLOC(2,K)<YLMIN) YLMIN=XLOC(2,K)
1688 IF (XLOC(2,K)>YLMAX) YLMAX=XLOC(2,K)
1689 ENDDO
1690 YLSEC=HALF*(YLMIN+YLMAX)
1691C
1692 NSEC=0
1693 DO K=1,LENPOLY(J)
1694 KK=K+1
1695 IF (K==LENPOLY(J)) KK=1
1696 X1=XLOC(1,K)
1697 Y1=XLOC(2,K)
1698 X2=XLOC(1,KK)
1699 Y2=XLOC(2,KK)
1700 IF (Y1-Y2/=ZERO) THEN
1701 ALPHA=(YLSEC-Y2)/(Y1-Y2)
1702.AND. IF (ALPHA>=ZEROALPHA<=ONE) THEN
1703 NSEC=NSEC+1
1704 XSEC(NSEC)=ALPHA*X1+(ONE-ALPHA)*X2
1705 ENDIF
1706 ELSE
1707 IF (Y1==YLSEC) THEN
1708 NSEC=NSEC+1
1709 XSEC(NSEC)=X1
1710 NSEC=NSEC+1
1711 XSEC(NSEC)=X2
1712 ENDIF
1713 ENDIF
1714 ENDDO
1715C
1716 XSMIN1=EP30
1717 KSMIN=-HUGE(KSMIN)
1718 DO K=1,NSEC
1719 IF (XSEC(K)<XSMIN1) THEN
1720 XSMIN1=XSEC(K)
1721 KSMIN=K
1722 ENDIF
1723 ENDDO
1724 XSMIN2=EP30
1725 DO K=1,NSEC
1726 IF (K==KSMIN) CYCLE
1727 IF (XSEC(K)<XSMIN2) XSMIN2=XSEC(K)
1728 ENDDO
1729C
1730 XS=HALF*(XSMIN1+XSMIN2)
1731 YS=YLSEC
1732 PHOL(1,NHOL)=RPOLY(5,J)+XS*VX1+YS*VX2
1733 PHOL(2,NHOL)=RPOLY(6,J)+XS*VY1+YS*VY2
1734 PHOL(3,NHOL)=RPOLY(7,J)+XS*VZ1+YS*VZ2
1735C ENDIF
1736 ENDDO
1737C
1738 IF (INFO==0) THEN
1739 IPOLY(2,I)=LENPOLY(I)
1740 IPOLY(6+LENPOLY(I)+1,I)=NHOL
1741 DO J=1,NHOL
1742 IPOLY(6+LENPOLY(I)+1+J,I)=IADHOL(J)
1743 RPOLY(4+3*LENPOLY(I)+3*(J-1)+1,I)=PHOL(1,J)
1744 RPOLY(4+3*LENPOLY(I)+3*(J-1)+2,I)=PHOL(2,J)
1745 RPOLY(4+3*LENPOLY(I)+3*(J-1)+3,I)=PHOL(3,J)
1746 ENDDO
1747 ENDIF
1748 ENDDO ! I=1,NPOLY
1749C
1750 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alpha2
#define alpha2
Definition eval.h:48
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
id
initmumps id
Definition multiplerhs_example.m:4
arret
subroutine arret(nn)
Definition arret.F:86