OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
facepoly.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!|| facepoly ../starter/source/airbag/facepoly.F
25!||--- called by ------------------------------------------------------
26!|| fvmesh1 ../starter/source/airbag/fvmesh.F
27!||--- calls -----------------------------------------------------
28!|| arret ../starter/source/system/arret.F
29!||====================================================================
30 SUBROUTINE facepoly(QUAD , TRI , NTRI , IPOLY, RPOLY ,
31 . N , NORMF, NORMT , NSMAX, NNP ,
32 . NRPMAX, NB , NV , LMIN , NPOLMAX,
33 . NPPMAX, INFO , IELNOD, X , IFAC ,
34 . ILVOUT, MVID )
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 Recensement des aretes de triangles
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 Calcul des intersections des aretes avec le plan de la facette
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 Calcul des intersections des triangles avec le plan de la facette (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 Intersection (eventuelle) de chaque segment avec les aretes de la facette
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 Intersection eventuelle des segments internes avec les segments externes
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 SI SS1=SS2=0 les segments sont colineaires (cas non traite)
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 Segments additionels sur aretes
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 Tri des points sur les bords de l'arete
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 des doubles
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 des segments de longueur nulle
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 du tableau REDIR
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 des segments sur les aretes
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 Recensement des noeuds
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 sur triangle interne
625 iseg(4,i)=2
626 iseg(5,i)=1
627 algo=1
628 ELSE
629C Segment sur triangle airbag ou brique additionelle
630 iseg(4,i)=1
631 iseg(5,i)=0
632 ENDIF
633 ELSE
634C Segment sur arete de boite
635 iseg(4,i)=1
636 iseg(5,i)=0
637 ENDIF
638 ENDDO
639C
640C Elimination des segments doubles
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 des polygones
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 Tableau noeuds-segments
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 Tri des noeuds
819C 4 segments : 2 internes + 2 externes
820C 3 segments : 3 internes
821C 3 segments : 1 interne + 2 externes
822C Test et mise des segments internes en tete
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 des polygones (nouvel algorithme v14)
974C JTAG(1,i) le polygone interne est du cote oppose a la normale au segment
975C JTAG(2,i) le polygone interne est du cote de la normale au segment
976C---------------------------------------------------------------------------
977 DO i=1,nseg
978 jtag(1,i)=0
979 jtag(2,i)=0
980 ENDDO
981 icmax=0
982C-----------------------------
983C Noeud connecte a 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 Traitement des 2 segments internes
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 Traitement des 2 segments externes
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 Noeud connecte a 3 segments : 1 interne et 2 externes
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 Normale au segment interne
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 Le cas peu probable alpha1=alpha2 est a completer
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 Noeud connecte a 3 segments internes
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 Normale au 1er segment interne
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 Noeud connecte a 2 segments internes
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 Noeud connecte a 2 segments externes
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 polygones
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 IF(temp(k) == isseg .OR. 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 Recherche des trous
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 Le premier point du polygone i est a l'interieur du polygone j
1577C Un point est interne si ABS(ALPHA)=2PI!
1578C On teste tous les autres
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 Le polygone j est un trou dans le polygone 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 Point interieur polygone 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 IF (alpha>=zero.AND.alpha<=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
1751 END
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
facepoly
subroutine facepoly(quad, tri, ntri, ipoly, rpoly, n, normf, normt, nsmax, nnp, nrpmax, nb, nv, lmin, npolmax, nppmax, info, ielnod, x, ifac, ilvout, mvid)
Definition facepoly.F:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
arret
subroutine arret(nn)
Definition arret.F:87