OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
polyhedr1.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!|| polyhedr1 ../starter/source/airbag/polyhedr1.F
25!||--- called by ------------------------------------------------------
26!|| fvmesh1 ../starter/source/airbag/fvmesh.F
27!||--- calls -----------------------------------------------------
28!|| arret ../starter/source/system/arret.F
29!|| c_tricall ../starter/stub/fvmbags_stub.F
30!||====================================================================
31 SUBROUTINE polyhedr1(IPOLY, RPOLY , POLB , NPOLB, POLH,
32 . NPOLH, NRPMAX , NPHMAX, IBRIC, LMIN,
33 . INFO , NPOLHMAX, NPPMAX, NEL, INZ ,
34 . TAGELA )
35C-----------------------------------------------
36C I m p l i c i t T y p e s
37C-----------------------------------------------
38#include "implicit_f.inc"
39C-----------------------------------------------
40C C o m m o n B l o c k s
41C-----------------------------------------------
42#include "units_c.inc"
43C-----------------------------------------------
44C D u m m y A r g u m e n t s
45C-----------------------------------------------
46 INTEGER NPPMAX, IPOLY(6+NPPMAX,*), POLB(*), NPOLB, NPHMAX,
47 . POLH(NPHMAX+2,*),NPOLH, NRPMAX, IBRIC, INFO, NPOLHMAX
48 INTEGER TAGELA(*)
49 INTEGER NEL, INZ
50 my_real
51 . rpoly(nrpmax,*), lmin
52C-----------------------------------------------
53C L o c a l V a r i a b l e s
54C-----------------------------------------------
55 INTEGER I, ICMAX, II, J, JJ, K, KK, IEL
56 INTEGER L, LL, ITY, JMIN, PMIN, POLD, NEDGE, NVERTEX
57 INTEGER IMIN, IMAX, I1, I2, I3, I4, N1, N2, K1, K2
58 INTEGER M1, M2, L1, L2
59 INTEGER ITAG(2,NPOLB), POLEDG(NPOLB,NPPMAX+1), ITYP(NPOLB)
60 INTEGER NTYP(3), REDIR(4), TEMP1(4), TEMP2(4)
61 INTEGER NNP, NHOL, NSEG, NELP, JTAG(3), NTRI3(NPOLB,NPPMAX)
62 INTEGER, ALLOCATABLE :: EDGPOL(:,:)
63 INTEGER, ALLOCATABLE :: EDGVER(:,:)
64 INTEGER, ALLOCATABLE :: EDGLOC(:,:)
65 INTEGER, ALLOCATABLE :: PSEG(:,:),PTRI(:,:)
66 my_real
67 . x1, y1, z1, x2, y2, z2, xx1, yy1, zz1, xx2, yy2, zz2,
68 . dd11, dd12, dd21, dd22, tole
69 my_real
70 . tst12, tst13, tst21, tst32, tst14,
71 . nx, ny, nz, nx1, ny1, nz1
72 my_real
73 . x0, y0, z0, nrm1, vx1, vy1, vz1, vx2, vy2, vz2,
74 . vx, vy, vz
75 my_real
76 . alpha13, alpha14
77 my_real, ALLOCATABLE :: xx(:), yy(:), zz(:)
78 my_real, ALLOCATABLE :: pnodes(:,:), pholes(:,:)
79C
80 tole=epsilon(zero)*0.5*lmin*lmin
81C--------------------------------
82C Triangularisation des polygones
83C--------------------------------
84 DO ii=1,npolb
85 i=polb(ii)
86 nnp=ipoly(2,i)
87 nhol=0
88 ALLOCATE(pnodes(2,nnp), pseg(2,nnp),ptri(3,2*nnp),pholes(2,1))
89 ptri(1:3,1:2*nnp) = 0
90C
91C Coordonnees des sommets dans le plan du polygone
92 nx=rpoly(2,i)
93 ny=rpoly(3,i)
94 nz=rpoly(4,i)
95C
96 x0=rpoly(5,i)
97 y0=rpoly(6,i)
98 z0=rpoly(7,i)
99 x1=rpoly(8,i)
100 y1=rpoly(9,i)
101 z1=rpoly(10,i)
102 nrm1=sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
103 vx1=(x1-x0)/nrm1
104 vy1=(y1-y0)/nrm1
105 vz1=(z1-z0)/nrm1
106 vx2=ny*vz1-nz*vy1
107 vy2=nz*vx1-nx*vz1
108 vz2=nx*vy1-ny*vx1
109C
110 pnodes(1,1)=zero
111 pnodes(2,1)=zero
112 DO j=2,nnp
113 x1=rpoly(4+3*(j-1)+1,i)
114 y1=rpoly(4+3*(j-1)+2,i)
115 z1=rpoly(4+3*(j-1)+3,i)
116 vx=x1-x0
117 vy=y1-y0
118 vz=z1-z0
119 pnodes(1,j)=vx*vx1+vy*vy1+vz*vz1
120 pnodes(2,j)=vx*vx2+vy*vy2+vz*vz2
121 pseg(1,j-1)=j-1
122 pseg(2,j-1)=j
123 ENDDO
124 pseg(1,nnp)=nnp
125 pseg(2,nnp)=1
126 nseg=nnp
127C
128 nelp = 0
129 CALL c_tricall(pnodes, pseg, pholes, ptri, nnp,
130 . nseg, nhol, nelp )
131CFA IF(NELP > NNP) THEN
132CFA WRITE(IOUT,'(A,I5,A,I2,2(I5,A))') 'BRICK',IBRIC,
133CFA . ' POLYGONE', I, NNP,' SOMMETS',NELP,' TRIANGLES'
134CFA DO J=1,NELP
135CFA WRITE(IOUT,'(4I5)') J,PTRI(1,J),PTRI(2,J),PTRI(3,J)
136CFA ENDDO
137CFA ENDIF
138C--------------------------------------------------------------------------------
139C Calcul pour chaque arete du 3eme noeud du triangle associe, stockage dans NTRI3
140C--------------------------------------------------------------------------------
141 DO j=1,nnp
142 n1=j
143 n2=j+1
144 IF(j==nnp) n2=1
145 DO k=1,nelp
146 jtag(k)=0
147 DO l=1,3
148 IF(ptri(l,k) == n1) jtag(k)=jtag(k)+1
149 IF(ptri(l,k) == n2) jtag(k)=jtag(k)+1
150 ENDDO
151 ENDDO
152 DO k=1,nelp
153 IF(jtag(k) /= 2) cycle
154 DO l=1,3
155 IF(ptri(l,k) == n1 .OR. ptri(l,k) == n2) cycle
156 ntri3(ii,j) = ptri(l,k)
157 ENDDO
158 ENDDO
159 ENDDO ! II=1,NPOLB
160C
161 DEALLOCATE(pnodes, pseg, ptri, pholes)
162 ENDDO
163C--------------------------------------------------------------
164C Construction polygone-arete
165C POLEDG(i,1) donne le nombre d'aretes connectees au polygone i
166C POLEDG(i,j+1) donne le numero de la jeme arete
167C EDGVER(i,1) numero du premier sommet de l'arete i
168C EDGVER(i,2) numero du deuxieme sommet de l'arete i
169C--------------------------------------------------------------
170 ALLOCATE (xx(npolb*nppmax*2))
171 ALLOCATE (yy(npolb*nppmax*2))
172 ALLOCATE (zz(npolb*nppmax*2))
173 ALLOCATE (edgver(npolb*nppmax,2))
174 DO i=1,npolb
175 DO j=1,nppmax+1
176 poledg(i,j)=0
177 ENDDO
178 ENDDO
179 nedge=0
180 nvertex=0
181 DO i=1,npolb
182 ii=polb(i)
183 poledg(i,1)=ipoly(2,ii)
184 DO j=1,ipoly(2,ii)
185 jj=j+1
186 IF(poledg(i,jj) > 0) cycle
187 nedge=nedge+1
188 poledg(i,jj)=nedge
189 x1=rpoly(4+3*(j-1)+1,ii)
190 y1=rpoly(4+3*(j-1)+2,ii)
191 z1=rpoly(4+3*(j-1)+3,ii)
192 nvertex=nvertex+1
193 edgver(nedge,1)=nvertex
194 xx(nvertex)=x1
195 yy(nvertex)=y1
196 zz(nvertex)=z1
197 IF (j==ipoly(2,ii)) jj=1
198 x2=rpoly(4+3*(jj-1)+1,ii)
199 y2=rpoly(4+3*(jj-1)+2,ii)
200 z2=rpoly(4+3*(jj-1)+3,ii)
201 nvertex=nvertex+1
202 edgver(nedge,2)=nvertex
203 xx(nvertex)=x2
204 yy(nvertex)=y2
205 zz(nvertex)=z2
206 DO k=1,npolb
207 IF (k==i) cycle
208 kk=polb(k)
209 DO l=1,ipoly(2,kk)
210 ll=l+1
211 IF (l==ipoly(2,kk)) ll=1
212 xx1=rpoly(4+3*(l-1)+1,kk)
213 yy1=rpoly(4+3*(l-1)+2,kk)
214 zz1=rpoly(4+3*(l-1)+3,kk)
215 xx2=rpoly(4+3*(ll-1)+1,kk)
216 yy2=rpoly(4+3*(ll-1)+2,kk)
217 zz2=rpoly(4+3*(ll-1)+3,kk)
218 dd11=(xx1-x1)**2+(yy1-y1)**2+(zz1-z1)**2
219 dd21=(xx2-x1)**2+(yy2-y1)**2+(zz2-z1)**2
220 dd12=(xx1-x2)**2+(yy1-y2)**2+(zz1-z2)**2
221 dd22=(xx2-x2)**2+(yy2-y2)**2+(zz2-z2)**2
222 IF ((dd11<=tole.AND.dd22<=tole).OR.
223 . (dd21<=tole.AND.dd12<=tole)) THEN
224 poledg(k,l+1)=nedge
225 ENDIF
226 ENDDO
227 ENDDO
228 ENDDO
229 ENDDO
230C---------------------------------------------------------------
231C Construction arete-polygone
232C EDGPOL(i,1) donne le nombre de polygones connectes a l'arete i
233C EDGPOL(i,j+1) donne le numero du jeme polygone
234C EDGLOC(i,k) donne la position de l'arete i dans POLEDG
235C---------------------------------------------------------------
236 ALLOCATE(edgpol(nedge,5))
237 ALLOCATE(edgloc(nedge,4))
238 DO i=1,nedge
239 DO j=1,4
240 edgpol(i,j)=0
241 edgloc(i,j)=0
242 ENDDO
243 edgpol(i,5)=0
244 ENDDO
245 DO i=1,npolb
246 DO j=1,poledg(i,1)
247 k=poledg(i,j+1)
248 IF(k==0) cycle
249 kk=edgpol(k,1)
250 kk=kk+1
251 IF(kk > 4) THEN
252 n1=edgver(k,1)
253 n2=edgver(k,2)
254 WRITE(iout,'(A,I5,A /A,1P3G20.13,A/A,1P3G20.13,A)')
255 . ' fvmbag mesh error : edge n1n2 is connected to ',KK,
256 . ' polygons but the limit is 4',
257 . ' n1 = (',XX(N1),YY(N1),ZZ(N1),')',
258 . ' n2 = (',XX(N2),YY(N2),ZZ(N2),')'
259 CALL ARRET(2)
260 ENDIF
261 EDGPOL(K,1)=KK
262 EDGPOL(K,KK+1)=I
263 EDGLOC(K,KK)=J
264 ENDDO
265 ENDDO
266C
267 DO I=1,NPOLB
268 ITAG(1,I)=0
269 ITAG(2,I)=0
270 II=POLB(I)
271 ITY=IPOLY(1,II)
272 IEL=IPOLY(3,II)
273.AND. IF(ITY == 1 TAGELA(IEL) > NEL) ITY=3
274 ITYP(I)=ITY
275 ENDDO
276C-------
277C Checks
278C-------
279 DO J=1,NEDGE
280 N1=EDGVER(J,1)
281 N2=EDGVER(J,2)
282 DO K=1,3
283 NTYP(K)=0
284 ENDDO
285 L=EDGPOL(J,1)
286 DO K=1,L
287 II=EDGPOL(J,K+1)
288 IF(ITYP(II) == 1) NTYP(1)=NTYP(1)+1
289 IF(ITYP(II) == 2) NTYP(2)=NTYP(2)+1
290 IF(ITYP(II) == 3) NTYP(3)=NTYP(3)+1
291 ENDDO
292 IF(L == 1) THEN
293 WRITE(IOUT,'(/2(a,i5),a,2(/a,1p3g20.13,a))')
294 . ' fvmbag mesh error : brick ',IBRIC,' layer ',INZ,
295 . ' edge n1n2 is connected to only 1 polygon',
296 . ' n1 = (',XX(N1),YY(N1),ZZ(N1),')',
297 . ' n2 = (',XX(N2),YY(N2),ZZ(N2),')'
298 CALL ARRET(2)
299.AND. ELSEIF(L == 2 NTYP(3) == 1) THEN
300 WRITE(IOUT,'(/2(a,i5),2a,2(/a,1p3g20.13,a))')
301 . ' fvmbag mesh error : brick ',IBRIC,' layer ',INZ,
302 . ' edge n1n2 is connected to 2 polygons but only 1',
303 . ' internal polygon',
304 . ' n1 = (',XX(N1),YY(N1),ZZ(N1),')',
305 . ' n2 = (',XX(N2),YY(N2),ZZ(N2),')'
306 CALL ARRET(2)
307.AND. ELSEIF(L == 3 NTYP(3) == 2) THEN
308 WRITE(IOUT,'(/2(a,i5),2a,2(/a,1p3g20.13,a))')
309 . ' fvmbag mesh error : brick ',IBRIC,' layer ',INZ,
310 . ' edge n1n2 is connected to 3 polygons 2 being',
311 . ' internal polygons',
312 . ' n1 = (',XX(N1),YY(N1),ZZ(N1),')',
313 . ' n2 = (',XX(N2),YY(N2),ZZ(N2),')'
314 CALL arret(2)
315 ELSEIF(l == 4) THEN
316 IF(ntyp(3) == 1) THEN
317 WRITE(iout,'(/2(A,I5),2A,2(/A,1P3G20.13,A))')
318 . ' FVMBAG MESH ERROR : BRICK ',ibric,' LAYER ',inz,
319 . ' EDGE N1N2 IS CONNECTED TO 4 POLYGONS ONLY 1 BEING',
320 . ' INTERNAL POLYGON',
321 . ' N1 = (',xx(n1),yy(n1),zz(n1),')',
322 . ' N2 = (',xx(n2),yy(n2),zz(n2),')'
323 CALL arret(2)
324 ELSEIF( ntyp(3) >= 3) THEN
325 WRITE(iout,'(/2(A,I5),2A,2(/A,1P3G20.13,A))')
326 . ' FVMBAG MESH ERROR : BRICK ',ibric,' LAYER ',inz,
327 . ' EDGE N1N2 IS CONNECTED TO 4 POLYGONS 3 BEING',
328 . ' INTERNAL POLYGONS',
329 . ' N1 = (',xx(n1),yy(n1),zz(n1),')',
330 . ' N2 = (',xx(n2),yy(n2),zz(n2),')'
331 CALL arret(2)
332 ENDIF
333 ELSEIF(l >= 5) THEN
334 WRITE(iout,'(/2(A,I5),A,2(/A,1P3G20.13,A))')
335 . ' fvmbag mesh error : brick ',IBRIC,' layer ',INZ,
336 . ' edge n1n2 is connected to more than 4 polygons',
337 . ' n1 = (',XX(N1),YY(N1),ZZ(N1),')',
338 . ' n2 = (',XX(N2),YY(N2),ZZ(N2),')'
339 CALL ARRET(2)
340 ENDIF
341 ENDDO ! J=1,NEDGE
342C---------------------------------------------------
343C Mettre les polygones internes en premiere position
344C---------------------------------------------------
345 DO J=1,NEDGE
346 IF(EDGPOL(J,1) == 2 ) CYCLE
347 DO K=1,EDGPOL(J,1)
348 TEMP1(K)=EDGPOL(J,K+1)
349 TEMP2(K)=EDGLOC(J,K)
350 ENDDO
351 I=1
352 DO K=1,EDGPOL(J,1)
353 L=EDGPOL(J,K+1)
354 IF(ITYP(L) == 3) THEN
355 REDIR(I)=K
356 I=I+1
357 ENDIF
358 ENDDO
359 DO K=1,EDGPOL(J,1)
360 L=EDGPOL(J,K+1)
361 IF(ITYP(L) /= 3) THEN
362 REDIR(I)=K
363 I=I+1
364 ENDIF
365 ENDDO
366 DO K=1,EDGPOL(J,1)
367 EDGPOL(J,K+1)=TEMP1(REDIR(K))
368 EDGLOC(J,K) =TEMP2(REDIR(K))
369 ENDDO
370 ENDDO ! J=1,NEDGE
371C--------------------------------------------------------------------
372C Construction des polyedres
373C ITAG(1,i) le polyedre est du cote oppose a la normale au polygone
374C ITAG(2,i) le polyedre est du cote de la normale au polygone
375C--------------------------------------------------------------------
376 ICMAX=0
377C------------------------------
378C Arete connectee a 4 polygones
379C------------------------------
380 DO J=1,NEDGE
381 IF(EDGPOL(J,1) < 4) CYCLE
382 I1=EDGPOL(J,2)
383 I2=EDGPOL(J,3)
384 I3=EDGPOL(J,4)
385 I4=EDGPOL(J,5)
386 II=POLB(I1)
387 NX=RPOLY(2,II)
388 NY=RPOLY(3,II)
389 NZ=RPOLY(4,II)
390 N1=EDGVER(J,1)
391 X1=XX(N1)
392 Y1=YY(N1)
393 Z1=ZZ(N1)
394 JJ=EDGLOC(J,2)
395 JJ=NTRI3(I2,JJ)
396 II=POLB(I2)
397 X2=RPOLY(4+3*(JJ-1)+1,II)
398 Y2=RPOLY(4+3*(JJ-1)+2,II)
399 Z2=RPOLY(4+3*(JJ-1)+3,II)
400 TST12=NX*(X2-X1)+NY*(Y2-Y1)+NZ*(Z2-Z1)
401C
402 JJ=EDGLOC(J,3)
403 JJ=NTRI3(I3,JJ)
404 II=POLB(I3)
405 NX1=RPOLY(2,II)
406 NY1=RPOLY(3,II)
407 NZ1=RPOLY(4,II)
408 X2=RPOLY(4+3*(JJ-1)+1,II)
409 Y2=RPOLY(4+3*(JJ-1)+2,II)
410 Z2=RPOLY(4+3*(JJ-1)+3,II)
411 VX=X2-X1
412 VY=Y2-Y1
413 VZ=Z2-Z1
414 TST13=NX*VX+NY*VY+NZ*VZ
415 ALPHA13=ACOS(NX*NX1+NY*NY1+NZ*NZ1)
416C
417 JJ=EDGLOC(J,4)
418 JJ=NTRI3(I4,JJ)
419 II=POLB(I4)
420 NX1=RPOLY(2,II)
421 NY1=RPOLY(3,II)
422 NZ1=RPOLY(4,II)
423 X2=RPOLY(4+3*(JJ-1)+1,II)
424 Y2=RPOLY(4+3*(JJ-1)+2,II)
425 Z2=RPOLY(4+3*(JJ-1)+3,II)
426 VX=X2-X1
427 VY=Y2-Y1
428 VZ=Z2-Z1
429 TST14=NX*VX+NY*VY+NZ*VZ
430 ALPHA14=ACOS(NX*NX1+NY*NY1+NZ*NZ1)
431C
432 II=POLB(I2)
433 NX=RPOLY(2,II)
434 NY=RPOLY(3,II)
435 NZ=RPOLY(4,II)
436 JJ=EDGLOC(J,1)
437 JJ=NTRI3(I1,JJ)
438 II=POLB(I1)
439 X2=RPOLY(4+3*(JJ-1)+1,II)
440 Y2=RPOLY(4+3*(JJ-1)+2,II)
441 Z2=RPOLY(4+3*(JJ-1)+3,II)
442 TST21=NX*(X2-X1)+NY*(Y2-Y1)+NZ*(Z2-Z1)
443C
444 IF(TST12*TST21 < ZERO) THEN
445 N2=EDGVER(J,2)
446 WRITE(IOUT,'(/2(a,i5),2a,2(/a,1p3g20.13,a))')
447 . ' fvmbag mesh error : brick ',IBRIC,' layer ',INZ,
448 . ' wrong normal orientation of internal surface elements',
449 . ' connected to edge n1n2',
450 . ' n1 = (',XX(N1),YY(N1),ZZ(N1),')',
451 . ' n2 = (',XX(N2),YY(N2),ZZ(N2),')'
452 CALL ARRET(2)
453 ELSE
454 ICMAX=ICMAX+1
455 ITAG(1,I1)=ICMAX
456 ITAG(1,I2)=ICMAX
457 ENDIF
458C
459.AND. IF(TST13 > ZERO TST14 > ZERO) THEN
460 IF(ALPHA13 < ALPHA14) THEN
461 ICMAX=ICMAX+1
462 ITAG(2,I1)=ICMAX
463 ITAG(1,I3)=ICMAX
464 ICMAX=ICMAX+1
465 ITAG(2,I2)=ICMAX
466 ITAG(1,I4)=ICMAX
467 ELSE
468 ICMAX=ICMAX+1
469 ITAG(2,I1)=ICMAX
470 ITAG(1,I4)=ICMAX
471 ICMAX=ICMAX+1
472 ITAG(2,I2)=ICMAX
473 ITAG(1,I3)=ICMAX
474 ENDIF
475.AND. ELSEIF(TST13 < ZERO TST14 < ZERO) THEN
476 IF(ALPHA13 < ALPHA14) THEN
477 ICMAX=ICMAX+1
478 ITAG(2,I1)=ICMAX
479 ITAG(1,I4)=ICMAX
480 ICMAX=ICMAX+1
481 ITAG(2,I2)=ICMAX
482 ITAG(1,I3)=ICMAX
483 ELSE
484 ICMAX=ICMAX+1
485 ITAG(2,I1)=ICMAX
486 ITAG(1,I3)=ICMAX
487 ICMAX=ICMAX+1
488 ITAG(2,I2)=ICMAX
489 ITAG(1,I4)=ICMAX
490 ENDIF
491.AND. ELSEIF(TST13 > ZERO TST14 < ZERO) THEN
492 ICMAX=ICMAX+1
493 ITAG(2,I1)=ICMAX
494 ITAG(1,I3)=ICMAX
495 ICMAX=ICMAX+1
496 ITAG(2,I2)=ICMAX
497 ITAG(1,I4)=ICMAX
498.AND. ELSEIF(TST13 < ZERO TST14 > ZERO) THEN
499 ICMAX=ICMAX+1
500 ITAG(2,I1)=ICMAX
501 ITAG(1,I4)=ICMAX
502 ICMAX=ICMAX+1
503 ITAG(2,I2)=ICMAX
504 ITAG(1,I3)=ICMAX
505 ENDIF
506 ENDDO ! J=1,NEDGE
507C---------------------------------------
508C Arete connectee a 3 polygones internes
509C---------------------------------------
510 DO J=1,NEDGE
511 IF(EDGPOL(J,1) /= 3) CYCLE
512 I1=EDGPOL(J,2)
513 I2=EDGPOL(J,3)
514 IF(ITYP(I2) /= 3) CYCLE
515 I3=EDGPOL(J,4)
516 N1=EDGVER(J,1)
517 X1=XX(N1)
518 Y1=YY(N1)
519 Z1=ZZ(N1)
520C Position du polygone I2 par rapport au polygone I1
521 II=POLB(I1)
522 NX=RPOLY(2,II)
523 NY=RPOLY(3,II)
524 NZ=RPOLY(4,II)
525 JJ=EDGLOC(J,2)
526 JJ=NTRI3(I2,JJ)
527 II=POLB(I2)
528 X2=RPOLY(4+3*(JJ-1)+1,II)
529 Y2=RPOLY(4+3*(JJ-1)+2,II)
530 Z2=RPOLY(4+3*(JJ-1)+3,II)
531 TST12=NX*(X2-X1)+NY*(Y2-Y1)+NZ*(Z2-Z1)
532 IF(TST12 < ZERO) THEN
533 K1=1
534 K2=2
535 ELSE
536 K1=2
537 K2=1
538 ENDIF
539C Position du polygone I1 par rapport au polygone I2
540 II=POLB(I2)
541 NX=RPOLY(2,II)
542 NY=RPOLY(3,II)
543 NZ=RPOLY(4,II)
544 JJ=EDGLOC(J,2)
545 JJ=NTRI3(I1,JJ)
546 II=POLB(I1)
547 X2=RPOLY(4+3*(JJ-1)+1,II)
548 Y2=RPOLY(4+3*(JJ-1)+2,II)
549 Z2=RPOLY(4+3*(JJ-1)+3,II)
550 TST21=NX*(X2-X1)+NY*(Y2-Y1)+NZ*(Z2-Z1)
551 IF(TST21 < ZERO) THEN
552 L1=1
553 L2=2
554 ELSE
555 L1=2
556 L2=1
557 ENDIF
558C Position du polygone I2 par rapport au polygone I3
559 II=POLB(I3)
560 NX=RPOLY(2,II)
561 NY=RPOLY(3,II)
562 NZ=RPOLY(4,II)
563 JJ=EDGLOC(J,2)
564 JJ=NTRI3(I2,JJ)
565 II=POLB(I2)
566 X2=RPOLY(4+3*(JJ-1)+1,II)
567 Y2=RPOLY(4+3*(JJ-1)+2,II)
568 Z2=RPOLY(4+3*(JJ-1)+3,II)
569 TST32=NX*(X2-X1)+NY*(Y2-Y1)+NZ*(Z2-Z1)
570 IF(TST32 < ZERO) THEN
571 M1=1
572 M2=2
573 ELSE
574 M1=2
575 M2=1
576 ENDIF
577C Premier polyedre
578.AND. IF (ITAG(K1,I1) == 0 ITAG(L1,I2) == 0) THEN
579 ICMAX=ICMAX+1
580 ITAG(K1,I1)=ICMAX
581 ITAG(L1,I2)=ICMAX
582.AND. ELSEIF(ITAG(K1,I1) == 0 ITAG(L1,I2) /= 0) THEN
583 ITAG(K1,I1)=ITAG(L1,I2)
584 ELSEIF(ITAG(K1,I1) /= 0. AND. ITAG(L1,I2) == 0) THEN
585 ITAG(L1,I2)=ITAG(K1,I1)
586.AND. ELSEIF(ITAG(K1,I1) /= 0 ITAG(L1,I2) /= 0) THEN
587 IF(ITAG(K1,I1) /= ITAG(L1,I2)) THEN
588 IMIN=MIN(ITAG(L1,I2),ITAG(K1,I1))
589 IMAX=MAX(ITAG(L1,I2),ITAG(K1,I1))
590 DO I=1,NPOLB
591 IF(ITAG(1,I) == IMAX) ITAG(1,I)=IMIN
592 IF(ITYP(I) /= 3) CYCLE
593 IF(ITAG(2,I) == IMAX) ITAG(2,I)=IMIN
594 ENDDO
595 ENDIF
596 ENDIF
597C Deuxieme polyedre
598.AND. IF(ITAG(M1,I3) == 0 ITAG(L2,I2) == 0) THEN
599 ICMAX=ICMAX+1
600 ITAG(M1,I3)=ICMAX
601 ITAG(L2,I2)=ICMAX
602.AND. ELSEIF(ITAG(M1,I3) /= 0 ITAG(L2,I2) == 0) THEN
603 ITAG(L2,I2)=ITAG(M1,I3)
604.AND. ELSEIF(ITAG(M1,I3) == 0 ITAG(L2,I2) /= 0) THEN
605 ITAG(M1,I3)=ITAG(L2,I2)
606.AND. ELSEIF(ITAG(M1,I3) /= 0 ITAG(L2,I2) /= 0) THEN
607 IF(ITAG(M1,I3) /= ITAG(L2,I2)) THEN
608 IMIN=MIN(ITAG(M1,I3),ITAG(L2,I2))
609 IMAX=MAX(ITAG(M1,I3),ITAG(L2,I2))
610 DO I=1,NPOLB
611 IF(ITAG(1,I) == IMAX) ITAG(1,I)=IMIN
612 IF(ITYP(I) /= 3) CYCLE
613 IF(ITAG(2,I) == IMAX) ITAG(2,I)=IMIN
614 ENDDO
615 ENDIF
616 ENDIF
617C Troisieme polyedre
618.AND. IF(ITAG(M2,I3) == 0 ITAG(K2,I1) == 0) THEN
619 ICMAX=ICMAX+1
620 ITAG(M2,I3)=ICMAX
621 ITAG(K2,I1)=ICMAX
622.AND. ELSEIF(ITAG(M2,I3) /= 0 ITAG(K2,I1) == 0) THEN
623 ITAG(K2,I1)=ITAG(M2,I3)
624.AND. ELSEIF(ITAG(M2,I3) == 0 ITAG(K2,I1) /= 0) THEN
625 ITAG(M2,I3)=ITAG(K2,I1)
626.AND. ELSEIF(ITAG(M2,I3) /= 0 ITAG(K2,I1) /= 0) THEN
627 IF(ITAG(M2,I3) /= ITAG(K2,I1)) THEN
628 IMIN=MIN(ITAG(M2,I3),ITAG(K2,I1))
629 IMAX=MAX(ITAG(M2,I3),ITAG(K2,I1))
630 DO I=1,NPOLB
631 IF(ITAG(1,I) == IMAX) ITAG(1,I)=IMIN
632 IF(ITYP(I) /= 3) CYCLE
633 IF(ITAG(2,I) == IMAX) ITAG(2,I)=IMIN
634 ENDDO
635 ENDIF
636 ENDIF
637 ENDDO ! J=1,NEDGE
638C-------------------------------------------------------
639C Arete connectee a 3 polygones (1 interne + 2 externes)
640C-------------------------------------------------------
641 DO J=1,NEDGE
642 IF(EDGPOL(J,1) /= 3) CYCLE
643 I1=EDGPOL(J,2)
644 I2=EDGPOL(J,3)
645 IF(ITYP(I2) == 3) CYCLE
646 I3=EDGPOL(J,4)
647 II=POLB(I1)
648 NX=RPOLY(2,II)
649 NY=RPOLY(3,II)
650 NZ=RPOLY(4,II)
651 N1=EDGVER(J,1)
652 X1=XX(N1)
653 Y1=YY(N1)
654 Z1=ZZ(N1)
655 JJ=EDGLOC(J,2)
656 JJ=NTRI3(I2,JJ)
657 II=POLB(I2)
658 X2=RPOLY(4+3*(JJ-1)+1,II)
659 Y2=RPOLY(4+3*(JJ-1)+2,II)
660 Z2=RPOLY(4+3*(JJ-1)+3,II)
661 TST12=NX*(X2-X1)+NY*(Y2-Y1)+NZ*(Z2-Z1)
662 IF(TST12 < ZERO) THEN
663 K1=1
664 K2=2
665 ELSE
666 K1=2
667 K2=1
668 ENDIF
669.AND. IF (ITAG(K1,I1) == 0 ITAG(1,I2) == 0) THEN
670 ICMAX=ICMAX+1
671 ITAG(K1,I1)=ICMAX
672 ITAG(1,I2)=ICMAX
673 IF(ITAG(1,I3) == 0) THEN
674 ICMAX=ICMAX+1
675 ITAG(K2,I1)=ICMAX
676 ITAG(1,I3)=ICMAX
677 ELSE
678 ITAG(K2,I1)=ITAG(1,I3)
679 ENDIF
680.AND. ELSEIF(ITAG(K1,I1) == 0 ITAG(1,I2) /= 0) THEN
681 ITAG(K1,I1)=ITAG(1,I2)
682 IF(ITAG(1,I3) == 0) THEN
683 ICMAX=ICMAX+1
684 ITAG(K2,I1)=ICMAX
685 ITAG(1,I3)=ICMAX
686 ELSE
687 ITAG(K2,I1)=ITAG(1,I3)
688 ENDIF
689 ELSEIF(ITAG(K1,I1) /= 0. AND. ITAG(1,I2) == 0) THEN
690 ITAG(1,I2)=ITAG(K1,I1)
691 IF(ITAG(1,I3) == 0) THEN
692 ITAG(1,I3)=ITAG(K2,I1)
693 ELSEIF(ITAG(1,I3) /= ITAG(K2,I1)) THEN
694 IMIN=MIN(ITAG(1,I3),ITAG(K2,I1))
695 IMAX=MAX(ITAG(1,I3),ITAG(K2,I1))
696 DO I=1,NPOLB
697 IF(ITAG(1,I) == IMAX) ITAG(1,I)=IMIN
698 IF(ITYP(I) /= 3) CYCLE
699 IF(ITAG(2,I) == IMAX) ITAG(2,I)=IMIN
700 ENDDO
701 ENDIF
702.AND. ELSEIF(ITAG(K1,I1) /= 0 ITAG(1,I2) /= 0) THEN
703 IF(ITAG(K1,I1) /= ITAG(1,I2)) THEN
704 IMIN=MIN(ITAG(1,I2),ITAG(K1,I1))
705 IMAX=MAX(ITAG(1,I2),ITAG(K1,I1))
706 DO I=1,NPOLB
707 IF(ITAG(1,I) == IMAX) ITAG(1,I)=IMIN
708 IF(ITYP(I) /= 3) CYCLE
709 IF(ITAG(2,I) == IMAX) ITAG(2,I)=IMIN
710 ENDDO
711 ENDIF
712 IF(ITAG(1,I3) == 0) THEN
713 ITAG(1,I3)=ITAG(K2,I1)
714 ELSEIF(ITAG(1,I3) /= ITAG(K2,I1)) THEN
715 IMIN=MIN(ITAG(1,I3),ITAG(K2,I1))
716 IMAX=MAX(ITAG(1,I3),ITAG(K2,I1))
717 DO I=1,NPOLB
718 IF(ITAG(1,I) == IMAX) ITAG(1,I)=IMIN
719 IF(ITYP(I) /= 3) CYCLE
720 IF(ITAG(2,I) == IMAX) ITAG(2,I)=IMIN
721 ENDDO
722 ENDIF
723 ENDIF
724 ENDDO ! J=1,NEDGE
725C------------------------------------------------
726C Arete connectee a 2 polygones
727C Les 2 polygones appartiennent au meme polyedre
728C------------------------------------------------
729 DO J=1,NEDGE
730 IF(EDGPOL(J,1) > 2) CYCLE
731 I1=EDGPOL(J,2)
732 I2=EDGPOL(J,3)
733.AND..OR. IF((ITYP(I1) /= 3 ITYP(I2) /= 3)
734.AND. . (ITYP(I1) == 3 ITYP(I2) == 3) ) THEN
735.AND. IF (ITAG(1,I1) == 0 ITAG(1,I2) == 0) THEN
736 ICMAX=ICMAX+1
737 ITAG(1,I1)=ICMAX
738 ITAG(1,I2)=ICMAX
739.AND. ELSEIF(ITAG(1,I1) == 0 ITAG(1,I2) /= 0) THEN
740 ITAG(1,I1)=ITAG(1,I2)
741 ELSEIF(ITAG(1,I1) /= 0. AND. ITAG(1,I2) == 0) THEN
742 ITAG(1,I2)=ITAG(1,I1)
743 ELSEIF(ITAG(1,I1) /= ITAG(1,I2)) THEN
744 IMIN=MIN(ITAG(1,I1),ITAG(1,I2))
745 IMAX=MAX(ITAG(1,I1),ITAG(1,I2))
746 DO I=1,NPOLB
747 IF(ITAG(1,I) == IMAX) ITAG(1,I)=IMIN
748 IF(ITYP(I) /= 3) CYCLE
749 IF(ITAG(2,I) == IMAX) ITAG(2,I)=IMIN
750 ENDDO
751 ENDIF
752 ENDIF
753C
754.AND. IF(ITYP(I1) == 3 ITYP(I2) == 3) THEN
755.AND. IF (ITAG(2,I1) == 0 ITAG(2,I2) == 0) THEN
756 ICMAX=ICMAX+1
757 ITAG(2,I1)=ICMAX
758 ITAG(2,I2)=ICMAX
759.AND. ELSEIF(ITAG(2,I1) == 0 ITAG(2,I2) /= 0) THEN
760 ITAG(2,I1)=ITAG(2,I2)
761 ELSEIF(ITAG(2,I1) /= 0. AND. ITAG(2,I2) == 0) THEN
762 ITAG(2,I2)=ITAG(2,I1)
763 ELSEIF(ITAG(2,I1) /= ITAG(2,I2)) THEN
764 IMIN=MIN(ITAG(2,I1),ITAG(2,I2))
765 IMAX=MAX(ITAG(2,I1),ITAG(2,I2))
766 DO I=1,NPOLB
767 IF(ITAG(1,I) == IMAX) ITAG(1,I)=IMIN
768 IF(ITYP(I) /= 3) CYCLE
769 IF(ITAG(2,I) == IMAX) ITAG(2,I)=IMIN
770 ENDDO
771 ENDIF
772 ENDIF
773 ENDDO ! J=1,NEDGE
774C
775 DEALLOCATE (EDGPOL)
776 DEALLOCATE (EDGLOC)
777 DEALLOCATE (EDGVER)
778 DEALLOCATE (XX,YY,ZZ)
779C
780 NPOLH=0
781 DO I=1,ICMAX
782 II=0
783 DO J=1,NPOLB
784.OR. IF (ITAG(1,J) == I ITAG(2,J) == I) II=II+1
785 ENDDO
786 IF (II/=0) NPOLH=NPOLH+1
787 ENDDO
788 IF (NPOLH>NPOLHMAX) THEN
789 INFO=1
790 NPOLHMAX=NPOLH
791 RETURN
792 ENDIF
793C
794 NPOLH=0
795 DO I=1,ICMAX
796 II=0
797 DO J=1,NPOLB
798.OR. IF (ITAG(1,J) == I ITAG(2,J) == I) II=II+1
799 ENDDO
800 IF (II/=0) THEN
801 NPOLH=NPOLH+1
802 POLH(1,NPOLH)=II
803 POLH(2,NPOLH)=IBRIC
804 II=0
805 DO J=1,NPOLB
806.OR. IF (ITAG(1,J) == I ITAG(2,J) == I) THEN
807 II=II+1
808 JJ=POLB(J)
809 POLH(2+II,NPOLH)=JJ
810 ITY=ITYP(J)
811 IF (ITY == 1) THEN
812 IPOLY(5,JJ)=NPOLH
813 ELSEIF(ITY == 2) THEN
814 IF (IPOLY(5,JJ)==0) THEN
815 IPOLY(5,JJ)=NPOLH
816 ELSE
817 IPOLY(6,JJ)=NPOLH
818 ENDIF
819 ELSEIF(ITY == 3) THEN
820 IF(ITAG(1,J) == I) THEN
821 IPOLY(5,JJ)=NPOLH
822 ELSEIF(ITAG(2,J) == I) THEN
823 IPOLY(6,JJ)=NPOLH
824 ENDIF
825 ENDIF
826 ENDIF
827 ENDDO
828C-------------------------------------------------------
829C Tri de POLH par ordre croissant de numero de polygones
830C-------------------------------------------------------
831 DO J=1,POLH(1,NPOLH)
832 JMIN=J
833 PMIN=POLH(2+J,NPOLH)
834 DO K=J+1,POLH(1,NPOLH)
835 IF (POLH(2+K,NPOLH)<PMIN) THEN
836 JMIN=K
837 PMIN=POLH(2+K,NPOLH)
838 ENDIF
839 ENDDO
840 POLD=POLH(2+J,NPOLH)
841 POLH(2+J,NPOLH)=PMIN
842 POLH(2+JMIN,NPOLH)=POLD
843 ENDDO !J=1,NPOLB
844 ENDIF
845 ENDDO !I=1,ICMAX
846C
847 INFO=0
848 RETURN
849 END
the
end diagonal values have been computed in the(sparse) matrix id.SOL
polyhedr1
subroutine polyhedr1(ipoly, rpoly, polb, npolb, polh, npolh, nrpmax, nphmax, ibric, lmin, info, npolhmax, nppmax, nel, inz, tagela)
Definition polyhedr1.F:35
arret
subroutine arret(nn)
Definition arret.F:87
c_tricall
subroutine c_tricall(pnodes, pseg, pholes, ptri, nnp, nseg, nhol, nelp)
Definition fvmbags_stub.F:34