OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fvmesh.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!|| fvmesh1 ../starter/source/airbag/fvmesh.F
25!||--- called by ------------------------------------------------------
26!|| fvmesh0 ../starter/source/airbag/fvmesh0.F
27!||--- calls -----------------------------------------------------
28!|| ancmsg ../starter/source/output/message/message.F
29!|| c_tricall ../starter/stub/fvmbags_stub.F
30!|| coorloc ../starter/source/airbag/fvmesh.F
31!|| facepoly ../starter/source/airbag/facepoly.F
32!|| fvbag_gpolcut ../starter/source/airbag/gpolcut.F
33!|| fvnormal ../starter/source/airbag/fvmbag1.F
34!|| horiedge ../starter/source/airbag/horipoly.F
35!|| horipoly ../starter/source/airbag/horipoly.F
36!|| itribox ../starter/source/airbag/itribox.F
37!|| my_exit ../starter/source/output/analyse/analyse.c
38!|| polyhedr ../starter/source/airbag/polyhedr.F
39!|| polyhedr1 ../starter/source/airbag/polyhedr1.F
40!|| tribox3 ../starter/stub/fvmbags_stub.F
41!|| tritri3 ../starter/stub/fvmbags_stub.F
42!||--- uses -----------------------------------------------------
43!|| fvbag_gpolcut_mod ../starter/source/airbag/gpolcut.F
44!|| fvbag_mod ../starter/share/modules1/fvbag_mod.F
45!|| message_mod ../starter/share/message_module/message_mod.F
46!|| reorder_mod ../starter/share/modules1/reorder_mod.F
47!||====================================================================
48 SUBROUTINE fvmesh1(IBUF , ELEM , X , IVOLU, BRIC ,
49 . XB , RVOLU , NEL , NELI, NBRIC, TBRIC,
50 . SFAC , DXM , NBA , NELA ,
51 . TBA , TFACA , TAGELS, IBUFA,
52 . ELEMA, TAGELA, IXS ,ID ,TITR, NB_NODE, ITYP)
53C-----------------------------------------------
54C M o d u l e s
55C-----------------------------------------------
56 USE message_mod
57 USE fvbag_mod
58 USE reorder_mod
59 USE names_and_titles_mod , ONLY : nchartitle
60 USE fvbag_gpolcut_mod, ONLY : fvbag_gpolcut
61C-----------------------------------------------
62C I m p l i c i t T y p e s
63C-----------------------------------------------
64#include "implicit_f.inc"
65C-----------------------------------------------
66C C o m m o n B l o c k s
67C-----------------------------------------------
68#include "units_c.inc"
69#include "sysunit.inc"
70C-----------------------------------------------
71C D u m m y A r g u m e n t s
72C-----------------------------------------------
73 INTEGER IBUF(*), ELEM(3,*), IVOLU(*), BRIC(8,*),
74 . NEL, NELI, NBRIC, TBRIC(13,*), NBA, NELA, TBA(2,*),
75 . TFACA(12,*), TAGELS(*), IBUFA(*), ELEMA(3,*), TAGELA(*),
76 . IXS(NIXS,*), NB_NODE, ITYP
77 my_real x(3,*), xb(3,*), rvolu(*), sfac(6,4,*), dxm
78 INTEGER ID
79 CHARACTER(len=nchartitle) :: TITR
80C-----------------------------------------------
81C L o c a l V a r i a b l e s
82C-----------------------------------------------
83 INTEGER ILVOUT, NLAYER, NFACMAX, NPPMAX, I, J, IEL, N1, N2, N3,
84 . NN1, NN2, NN3, IAD, NPOLY, NNS, ITAGB(NBRIC),
85 . NVMAX, INFO, ICUT, NVERTS, NF1, NF2, NF3, NF4, NS,
86 . NSMAX, NV, K, KK, NNP, NPOLMAX, NRPMAX, N, M, MM, NPOLH,
87 . npolb, ity, nn, nphmax, npolhmax, jj, jjj, ii,
88 . nns_old, nnp_old, npa, jmax, imax, jjmax, np, nntr,
89 . npoly_old, ipsurf, ic1, ic2, nhol, nelp, npolh_old,
90 . l, ll, ifv, nns2, npoly2, npl, polc, nspoly, ncpoly,
91 . npolhf, lenp, lenh, ip, lenimax,
92 . lenrmax, kkk, nseg, imin, nfac, n4, nn4,
93 . itagba(nba), ibsa(nba), itagn(nb_node), nall,
94 . nns_anim, nbz, nbnedge, nns3, npoly_new, nnsp, nedge,
95 . ityz, inz, j0, nns1, k0, i1, i2, idep, nstr, nctr, iiz,
96 . nns_old_save,ityzint,iad1, iad2, itmp1
97 INTEGER, DIMENSION(:), ALLOCATABLE :: IFLAG
98 PARAMETER (NVMAX=12)
99 my_real x1, y1, z1, x2, y2, z2, x3, y3, z3, x12, y12, z12,
100 . x13, y13, z13, nrx, nry, nrz, area2, elarea(nel+neli),
101 . norm(3,nel+neli), xbmax, ybmax, zbmax, xbmin, ybmin, zbmin,
102 . xc, yc, zc, xx, yy, zz, pp(3,3), xl7(3), bcenter(3),
103 . bhalfsize(3), xtmax, ytmax, ztmax, xtmin, ytmin, ztmin,
104 . tverts(9), tmptri(3,3), tmpbox(3,8), tmpnorm(3,6), tole,
105 . xg, yg, zg, fv0(3), fv1(3), fv2(3), fu0(3), fu1(3),
106 . fu2(3), quad(3,4), nr(3), area, nx, ny, nz, nnx,
107 . nny, nnz, x0, y0, z0, lmax2, tole2, dd, vm, volumin,
108 . areamax, volph, areap, areael, vpx, vpy, vpz,
109 . v1x, v1y, v1z, v2x, v2y, v2z, nrm1, vx, vy, vz, ss,
110 . normf(3,4), ksi, eta, vx1, vy1, vz1, vx2, vy2, vz2,
111 . vmin, x4, y4, z4, norma(3,nela),
112 . dd2, xn(3), zlmin, zlmax, zl, vnorm, vx3, vy3, vz3, lz,
113 . dz, alpha, gamai, cpai, cpbi, cpci, rmwi, pini, ti, cpi,
114 . cvi, rhoi, zl1, zl2, zl3, zlc, val, tmass,
115 . ti2, cpdi, cpei, cpfi, r1
116 INTEGER :: COUNT_ELEM_INT, COUNT_ELEM_EXT, COUNT_TRIA_INTERNE
117 INTEGER :: COUNT_ELEM_INT_OLD, COUNT_ELEM_EXT_OLD, COUNT_TRIA_INTERNE_OLD
118 INTEGER, ALLOCATABLE :: FACET(:,:), IPOLY(:,:),
119 . IELNOD(:,:), POLH(:,:), IPOLY_F(:,:),
120 . polh_f(:,:), ifvnod(:), ifvnod_old(:),
121 . redir_poly(:), pseg(:,:), redir(:),
122 . ptri(:,:), redir_polh(:), polb(:),
123 . ipoly_old(:), polh_new(:,:), itagt(:),
124 . tri(:,:), adr(:), adr_old(:), imerged(:),
125 . ipoly_f_old(:,:), inedge(:,:), nref(:,:),
126 . iz(:,:), ledge(:), ifvnod2(:,:)
127 my_real, ALLOCATABLE :: normt(:,:), poly(:,:), rpoly(:,:),
128 . rpoly_f(:,:), volu(:), pnodes(:,:),
129 . pholes(:,:), rpoly_old(:),
130 . volu_old(:), rpoly_f_old(:,:),
131 . rnedge(:,:),
132 . aref(:,:), rfvnod2(:,:), xns(:,:),
133 . xns2(:,:), areatr(:)
134
135 INTEGER FAC(4,6), FACNOR(4,6), FAC4(3,4), FAC6(4,5)
136 INTEGER FAC5(4,5), NFACE(4), NTYPE, IQUAD
137 DATA FAC /1,4,3,2,
138 . 5,6,7,8,
139 . 1,2,6,5,
140 . 2,3,7,6,
141 . 3,4,8,7,
142 . 4,1,5,8/
143 DATA facnor /3,4,5,6,
144 . 3,4,5,6,
145 . 1,4,2,6,
146 . 1,5,2,3,
147 . 1,6,2,4,
148 . 1,3,2,5/
149 DATA fac4 /1,5,3,
150 . 3,5,6,
151 . 6,5,1,
152 . 1,3,6/
153 DATA fac6 /1,3,2,0,
154 . 5,6,7,0,
155 . 1,2,6,5,
156 . 2,3,7,6,
157 . 3,4,8,7/
158 DATA fac5 /1,2,5,0,
159 . 2,3,5,0,
160 . 3,4,5,0,
161 . 4,1,5,0,
162 . 1,4,3,2/
163 DATA nface/6,4,5,5/
164C-------------------------
165C liste chainee polygones
166 TYPE polygone
167 INTEGER IZ(3), NNSP
168 INTEGER, DIMENSION(:), POINTER :: IPOLY, IELNOD
169 INTEGER, DIMENSION(:,:), POINTER :: NREF
170 my_real, DIMENSION(:), POINTER :: rpoly
171 my_real, DIMENSION(:,:), POINTER :: aref
172 TYPE(polygone), POINTER :: PTR
173 END TYPE polygone
174 TYPE(polygone), POINTER :: FIRST, PTR_PREC, PTR_CUR, PTR_TMP,
175 . ptr_old, first2, ptr_cur2, ptr_prec2,
176 . ptr_tmp2
177 TYPE(polygone), TARGET :: NOTHING
178C liste chainee polyedres
179 TYPE polyhedre
180 INTEGER RANK
181 INTEGER, DIMENSION(:), POINTER :: POLH
182 TYPE(POLYHEDRE), POINTER :: PTR
183 END TYPE polyhedre
184 TYPE(polyhedre), POINTER :: PHFIRST, PH_PREC, PH_CUR, PH_TMP
185
186 my_real, DIMENSION(:), ALLOCATABLE :: radius
187 INTEGER, DIMENSION(:), ALLOCATABLE :: IDX
188C
189C Structure de donnees maillage volumes finis
190C - NNS : Nombre de nouveaux noeuds ( = IVOLU(46) )
191C - NNTR : Nombre de nouveaux triangles ( = IVOLU(47) )
192C - NPOLY : Nombre de polygones ( = IVOLU(48) )
193C - NPOLH : Nombre de polyedres ( volumes finis) ( = IVOLU(49) )
194C
195C - IFVNOD(NNS) : Triangle airbag de reference pour le nouveau noeud
196C - RFVNOD(2,NNS) : Coordonnees locales du noeud dans le triangle
197C
198C - IFVTRI(6,NNTR) : Connectivite nouveaux triangles (1->3)
199C Element de reference si appuye sur surface (4)
200C Polyedres connectes si interne (5->6)
201C
202C - IFVPOLY(NNTR) : Attribution triangles-polygones
203C - IFVTADR(NPOLY+1) : Adresses dans le tableau precedent
204C
205C - IFVPOLH(NPOLY) : Attribution polygones-polyedres
206C - IFVPADR(NPOLH+1) : Adresses dans le tableau precedent
207C
208C
209 nothing%PTR => nothing
210 ptr_prec => nothing
211 iquad = 0
212 n4 = 0
213 nx = 0
214 ny = 0
215 nz = 0
216 ipsurf = 0
217 NULLIFY(first)
218 NULLIFY(phfirst)
219 NULLIFY(ptr_cur)
220C
221 ilvout=ivolu(44)
222 IF (ilvout >=1.AND.nba == 0) WRITE(istdo,'(A19,I10)')
223 . ' --> FVMBAG ID: ',ivolu(1)
224C
225 nlayer=ivolu(40)
226 nfacmax=ivolu(41)
227 nppmax=ivolu(42)
228C
229 DO i=1,nbric
230 tbric(7,i)=0
231 DO j=1,6
232 tbric(7+j,i)=0
233 ENDDO
234 ENDDO
235C---------------------------------------------------
236C Normales elementaires
237C---------------------------------------------------
238 DO iel=1,nel+neli
239 n1=elem(1,iel)
240 n2=elem(2,iel)
241 n3=elem(3,iel)
242 nn1=ibuf(n1)
243 nn2=ibuf(n2)
244 nn3=ibuf(n3)
245 CALL fvnormal(x,nn1,nn2,nn3,0,nrx,nry,nrz)
246 area2=sqrt(nrx**2+nry**2+nrz**2)
247 elarea(iel)=half*area2
248 norm(1,iel)=nrx/area2
249 norm(2,iel)=nry/area2
250 norm(3,iel)=nrz/area2
251 ENDDO
252C
253 DO iel=1,nela
254 n1=elema(1,iel)
255 n2=elema(2,iel)
256 n3=elema(3,iel)
257 IF (tagela(iel) > 0) THEN
258 nn1=ibuf(n1)
259 nn2=ibuf(n2)
260 nn3=ibuf(n3)
261 ELSEIF (tagela(iel) < 0) THEN
262 nn1=ibufa(n1)
263 nn2=ibufa(n2)
264 nn3=ibufa(n3)
265 ENDIF
266 CALL fvnormal(x,nn1,nn2,nn3,0,nrx,nry,nrz)
267 area2=sqrt(nrx**2+nry**2+nrz**2)
268 norma(1,iel)=nrx/area2
269 norma(2,iel)=nry/area2
270 norma(3,iel)=nrz/area2
271 ENDDO
272C---------------------------------------------------
273C Creation de tous les polygones
274C---------------------------------------------------
275 npoly=0
276 nns=0
277 npoly2=0
278 nns2=0
279 IF (nela == 0) THEN
280 npoly=0
281 npolh=0
282 GOTO 370
283 ENDIF
284C
285 DO i=1,nbric
286 itagb(i)=0
287 ENDDO
288C 1. Facettes airbag coupees
289C nullify sinon test sur associated vrai second appel fvmesh1
290 NULLIFY(first2)
291 IF (ilvout >=1) WRITE(istdo,'(8X,A)') 'BUILDING POLYGONS'
292C
293 DO i=1,nbric
294 !IF (ILVOUT >=1) CALL PROGALE_C(I, NBRIC, 1) !dynamic screen output
295C
296 ptr_old => ptr_cur
297 npoly_old=npoly
298 nns_old=nns
299 100 CONTINUE
300C
301 ALLOCATE(facet(6,1+nfacmax), tri(nfacmax,5),
302 . normt(3,nfacmax), poly(3,nvmax))
303 DO j=1,6
304 facet(j,1)=0
305 ENDDO
306C
307 xbmax=-ep30
308 ybmax=-ep30
309 zbmax=-ep30
310 xbmin=ep30
311 ybmin=ep30
312 zbmin=ep30
313 xc=zero
314 yc=zero
315 zc=zero
316 DO j=1,8
317 xx=xb(1,bric(j,i))
318 yy=xb(2,bric(j,i))
319 zz=xb(3,bric(j,i))
320 xbmax=max(xbmax,xx)
321 ybmax=max(ybmax,yy)
322 zbmax=max(zbmax,zz)
323 xbmin=min(xbmin,xx)
324 ybmin=min(ybmin,yy)
325 zbmin=min(zbmin,zz)
326 xc=xc+one_over_8*xx
327 yc=yc+one_over_8*yy
328 zc=zc+one_over_8*zz
329 ENDDO
330C Passage repere local brique
331 pp(1,1)=sfac(4,2,i)
332 pp(2,1)=sfac(4,3,i)
333 pp(3,1)=sfac(4,4,i)
334 pp(1,2)=sfac(5,2,i)
335 pp(2,2)=sfac(5,3,i)
336 pp(3,2)=sfac(5,4,i)
337 pp(1,3)=sfac(2,2,i)
338 pp(2,3)=sfac(2,3,i)
339 pp(3,3)=sfac(2,4,i)
340C
341 xx=xb(1,bric(7,i))
342 yy=xb(2,bric(7,i))
343 zz=xb(3,bric(7,i))
344 xl7(1)=pp(1,1)*(xx-xc)+pp(2,1)*(yy-yc)+pp(3,1)*(zz-zc)
345 xl7(2)=pp(1,2)*(xx-xc)+pp(2,2)*(yy-yc)+pp(3,2)*(zz-zc)
346 xl7(3)=pp(1,3)*(xx-xc)+pp(2,3)*(yy-yc)+pp(3,3)*(zz-zc)
347 bcenter(1)=zero
348 bcenter(2)=zero
349 bcenter(3)=zero
350 bhalfsize(1)=xl7(1)
351 bhalfsize(2)=xl7(2)
352 bhalfsize(3)=xl7(3)
353 info=0
354 DO iel=1,nela
355 n1=elema(1,iel)
356 n2=elema(2,iel)
357 n3=elema(3,iel)
358 IF (tagela(iel) > 0) THEN
359 nn1=ibuf(n1)
360 nn2=ibuf(n2)
361 nn3=ibuf(n3)
362 ELSEIF (tagela(iel) < 0) THEN
363 nn1=ibufa(n1)
364 nn2=ibufa(n2)
365 nn3=ibufa(n3)
366 ENDIF
367 x1=x(1,nn1)
368 y1=x(2,nn1)
369 z1=x(3,nn1)
370 x2=x(1,nn2)
371 y2=x(2,nn2)
372 z2=x(3,nn2)
373 x3=x(1,nn3)
374 y3=x(2,nn3)
375 z3=x(3,nn3)
376 xtmax=max(x1,x2,x3)
377 ytmax=max(y1,y2,y3)
378 ztmax=max(z1,z2,z3)
379 xtmin=min(x1,x2,x3)
380 ytmin=min(y1,y2,y3)
381 ztmin=min(z1,z2,z3)
382 IF (xbmax < xtmin.OR.ybmax < ytmin.OR.zbmax < ztmin.OR.
383 . xbmin > xtmax.OR.ybmin > ytmax.OR.zbmin > ztmax)
384 . cycle
385C
386 tverts(1)=pp(1,1)*(x1-xc)+pp(2,1)*(y1-yc)+pp(3,1)*(z1-zc)
387 tverts(2)=pp(1,2)*(x1-xc)+pp(2,2)*(y1-yc)+pp(3,2)*(z1-zc)
388 tverts(3)=pp(1,3)*(x1-xc)+pp(2,3)*(y1-yc)+pp(3,3)*(z1-zc)
389 tverts(4)=pp(1,1)*(x2-xc)+pp(2,1)*(y2-yc)+pp(3,1)*(z2-zc)
390 tverts(5)=pp(1,2)*(x2-xc)+pp(2,2)*(y2-yc)+pp(3,2)*(z2-zc)
391 tverts(6)=pp(1,3)*(x2-xc)+pp(2,3)*(y2-yc)+pp(3,3)*(z2-zc)
392 tverts(7)=pp(1,1)*(x3-xc)+pp(2,1)*(y3-yc)+pp(3,1)*(z3-zc)
393 tverts(8)=pp(1,2)*(x3-xc)+pp(2,2)*(y3-yc)+pp(3,2)*(z3-zc)
394 tverts(9)=pp(1,3)*(x3-xc)+pp(2,3)*(y3-yc)+pp(3,3)*(z3-zc)
395C
396 CALL tribox3(icut, bcenter, bhalfsize, tverts)
397C Liste elargie pour faces laterales
398 IF (icut == 0) THEN
399 tole=em5
400* TOLE=ZERO
401 xg=third*(x1+x2+x3)
402 yg=third*(y1+y2+y3)
403 zg=third*(z1+z2+z3)
404 x1=xg+(one+tole)*(x1-xg)
405 y1=yg+(one+tole)*(y1-yg)
406 z1=zg+(one+tole)*(z1-zg)
407 x2=xg+(one+tole)*(x2-xg)
408 y2=yg+(one+tole)*(y2-yg)
409 z2=zg+(one+tole)*(z2-zg)
410 x3=xg+(one+tole)*(x3-xg)
411 y3=yg+(one+tole)*(y3-yg)
412 z3=zg+(one+tole)*(z3-zg)
413C
414 tverts(1)=pp(1,1)*(x1-xc)+pp(2,1)*(y1-yc)+pp(3,1)*(z1-zc)
415 tverts(2)=pp(1,2)*(x1-xc)+pp(2,2)*(y1-yc)+pp(3,2)*(z1-zc)
416 tverts(3)=pp(1,3)*(x1-xc)+pp(2,3)*(y1-yc)+pp(3,3)*(z1-zc)
417 tverts(4)=pp(1,1)*(x2-xc)+pp(2,1)*(y2-yc)+pp(3,1)*(z2-zc)
418 tverts(5)=pp(1,2)*(x2-xc)+pp(2,2)*(y2-yc)+pp(3,2)*(z2-zc)
419 tverts(6)=pp(1,3)*(x2-xc)+pp(2,3)*(y2-yc)+pp(3,3)*(z2-zc)
420 tverts(7)=pp(1,1)*(x3-xc)+pp(2,1)*(y3-yc)+pp(3,1)*(z3-zc)
421 tverts(8)=pp(1,2)*(x3-xc)+pp(2,2)*(y3-yc)+pp(3,2)*(z3-zc)
422 tverts(9)=pp(1,3)*(x3-xc)+pp(2,3)*(y3-yc)+pp(3,3)*(z3-zc)
423C
424 CALL tribox3(icut, bcenter, bhalfsize, tverts)
425 IF (icut == 1) icut=2
426 ENDIF
427C
428 IF (icut >= 1) THEN
429 IF (icut == 1) THEN
430 tbric(7,i)=2
431 tmptri(1,1)=x1
432 tmptri(2,1)=y1
433 tmptri(3,1)=z1
434 tmptri(1,2)=x2
435 tmptri(2,2)=y2
436 tmptri(3,2)=z2
437 tmptri(1,3)=x3
438 tmptri(2,3)=y3
439 tmptri(3,3)=z3
440 DO j=1,8
441 tmpbox(1,j)=xb(1,bric(j,i))
442 tmpbox(2,j)=xb(2,bric(j,i))
443 tmpbox(3,j)=xb(3,bric(j,i))
444 ENDDO
445 DO j=1,6
446 tmpnorm(1,j)=-sfac(j,2,i)
447 tmpnorm(2,j)=-sfac(j,3,i)
448 tmpnorm(3,j)=-sfac(j,4,i)
449 ENDDO
450 CALL itribox(tmptri, tmpbox, tmpnorm, nverts, poly,
451 . nvmax )
452 IF (nverts > 2) THEN
453 npoly=npoly+1
454 IF (.NOT.ASSOCIATED(first)) THEN
455 ALLOCATE(first)
456 ptr_cur => first
457 ELSE
458 ALLOCATE(ptr_cur)
459 ptr_prec%PTR => ptr_cur
460 ENDIF
461 ALLOCATE(ptr_cur%IPOLY(6+nverts),
462 . ptr_cur%IELNOD(nverts),
463 . ptr_cur%RPOLY(4+3*nverts))
464 ptr_cur%IPOLY(1)=1
465 ptr_cur%IPOLY(2)=nverts
466 ptr_cur%IPOLY(3)=iel
467 ptr_cur%IPOLY(4)=i
468 ptr_cur%IPOLY(5)=0
469 ptr_cur%IPOLY(6)=0
470 ptr_cur%RPOLY(1)=zero
471 ptr_cur%RPOLY(2)=norma(1,iel)
472 ptr_cur%RPOLY(3)=norma(2,iel)
473 ptr_cur%RPOLY(4)=norma(3,iel)
474 DO j=1,nverts
475 nns=nns+1
476 ptr_cur%IPOLY(6+j)=nns
477 ptr_cur%IELNOD(j)=iel
478 ptr_cur%RPOLY(4+3*(j-1)+1)=poly(1,j)
479 ptr_cur%RPOLY(4+3*(j-1)+2)=poly(2,j)
480 ptr_cur%RPOLY(4+3*(j-1)+3)=poly(3,j)
481 ENDDO
482 ptr_prec => ptr_cur
483 ENDIF
484 ENDIF
485C Facettes coupees (1 quadrangle=2 triangles)
486 tole=em5
487* TOLE=ZERO
488 xg=third*(x1+x2+x3)
489 yg=third*(y1+y2+y3)
490 zg=third*(z1+z2+z3)
491 x1=xg+(one+tole)*(x1-xg)
492 y1=yg+(one+tole)*(y1-yg)
493 z1=zg+(one+tole)*(z1-zg)
494 x2=xg+(one+tole)*(x2-xg)
495 y2=yg+(one+tole)*(y2-yg)
496 z2=zg+(one+tole)*(z2-zg)
497 x3=xg+(one+tole)*(x3-xg)
498 y3=yg+(one+tole)*(y3-yg)
499 z3=zg+(one+tole)*(z3-zg)
500C
501 fv0(1)=x1
502 fv0(2)=y1
503 fv0(3)=z1
504 fv1(1)=x2
505 fv1(2)=y2
506 fv1(3)=z2
507 fv2(1)=x3
508 fv2(2)=y3
509 fv2(3)=z3
510 DO j=1,6
511 nf1=bric(fac(1,j),i)
512 nf2=bric(fac(2,j),i)
513 nf3=bric(fac(3,j),i)
514 nf4=bric(fac(4,j),i)
515 fu0(1)=xb(1,nf1)
516 fu0(2)=xb(2,nf1)
517 fu0(3)=xb(3,nf1)
518 fu1(1)=xb(1,nf2)
519 fu1(2)=xb(2,nf2)
520 fu1(3)=xb(3,nf2)
521 fu2(1)=xb(1,nf3)
522 fu2(2)=xb(2,nf3)
523 fu2(3)=xb(3,nf3)
524C
525 CALL tritri3(icut, fv0, fv1, fv2, fu0, fu1, fu2)
526 IF (icut == 1) THEN
527 tbric(7+j,i)=2
528 ns=facet(j,1)
529 ns=ns+1
530 facet(j,1)=ns
531 IF (ns > nfacmax) THEN
532 info=1
533 ELSE
534 facet(j,1+ns)=iel
535 ENDIF
536 cycle
537 ENDIF
538C
539 fu1(1)=xb(1,nf3)
540 fu1(2)=xb(2,nf3)
541 fu1(3)=xb(3,nf3)
542 fu2(1)=xb(1,nf4)
543 fu2(2)=xb(2,nf4)
544 fu2(3)=xb(3,nf4)
545C
546 CALL tritri3(icut, fv0, fv1, fv2, fu0, fu1, fu2)
547 IF (icut == 1) THEN
548 tbric(7+j,i)=2
549 ns=facet(j,1)
550 ns=ns+1
551 facet(j,1)=ns
552 IF (ns > nfacmax) THEN
553 info=1
554 ELSE
555 facet(j,1+ns)=iel
556 ENDIF
557 ENDIF
558 ENDDO
559 ENDIF
560 ENDDO
561C
562 IF (info == 1) THEN
563 nsmax=0
564 DO j=1,6
565 nsmax=max(nsmax,facet(j,1))
566 ENDDO
567 nfacmax=nsmax
568 IF(ilvout >=1) WRITE(iout,'(A,I10,A,I10)')
569 . ' MONITORED VOLUME ID: ',ivolu(1),
570 . ' NFACMAX IS RESET TO ',nsmax
571C Destruction des polygones crees pour rien
572 IF (npoly_old > 0) THEN
573 ptr_cur => ptr_old%PTR
574 ELSE
575 ptr_cur => first
576 ENDIF
577 DO j=1,npoly-npoly_old
578 ptr_tmp => ptr_cur%PTR
579 DEALLOCATE(ptr_cur%IPOLY, ptr_cur%RPOLY,
580 . ptr_cur%IELNOD)
581 DEALLOCATE(ptr_cur)
582 ptr_cur => ptr_tmp
583 ENDDO
584 IF (npoly_old > 0) THEN
585 ptr_prec => ptr_old
586 ELSE
587 NULLIFY(first)
588 ENDIF
589 npoly=npoly_old
590 nns=nns_old
591C
592C
593 DEALLOCATE(facet, tri, normt, poly)
594 GOTO 100
595 ENDIF
596C 2. Faces laterales
597 IF (tbric(7,i) <= 1) THEN
598 DEALLOCATE(facet, tri, normt, poly)
599 cycle
600 END IF
601C Debut boucle sur les 6 faces de la brique i
602 DO j=1,6
603 nv=tbric(j,i)
604 IF (nv == 0) cycle
605 IF (itagb(nv) == 1) cycle
606 IF (tbric(7+j,i) == 2) THEN
607 DO k=1,4
608 kk=fac(k,j)
609 quad(1,k)=xb(1,bric(kk,i))
610 quad(2,k)=xb(2,bric(kk,i))
611 quad(3,k)=xb(3,bric(kk,i))
612 ENDDO
613 ns=facet(j,1)
614 DO k=1,ns
615 iel=facet(j,1+k) ! 1<IEL< NELA
616 IF (tagela(iel) > 0) THEN
617 DO l=1,3
618 ll=ibuf(elema(l,iel))
619 tri(k,l)=ll
620 ENDDO
621 IF(tagela(iel) <= nel) THEN
622C Triangle sur airbag
623 tri(k,5)=1
624 ELSE
625C Triangle sur structure interne
626 tri(k,5)=3
627 ENDIF
628 ELSEIF (tagela(iel) < 0) THEN
629C Triangle sur brique additionelle
630 DO l=1,3
631 ll=ibufa(elema(l,iel))
632 tri(k,l)=ll
633 ENDDO
634 tri(k,5)=2
635 ENDIF
636 tri(k,4)=iel
637 normt(1,k)=norma(1,iel)
638 normt(2,k)=norma(2,iel)
639 normt(3,k)=norma(3,iel)
640 ENDDO
641 DO k=1,4
642 normf(1,k)=sfac(facnor(k,j),2,i)
643 normf(2,k)=sfac(facnor(k,j),3,i)
644 normf(3,k)=sfac(facnor(k,j),4,i)
645 ENDDO
646 nr(1)=sfac(j,2,i)
647 nr(2)=sfac(j,3,i)
648 nr(3)=sfac(j,4,i)
649C
650 nnp=0
651 npolmax=nlayer
652C
653 200 CONTINUE
654 nrpmax=nppmax*3+4
655 ALLOCATE(ipoly(6+nppmax+1+npolmax,npolmax),
656 . rpoly(nrpmax+3*npolmax,npolmax),
657 . ielnod(nppmax,npolmax))
658C
659 CALL facepoly(quad, tri, ns, ipoly, rpoly,
660 . nr, normf, normt, nfacmax, nnp,
661 . nrpmax, i, nv, dxm, npolmax,
662 . nppmax, info, ielnod, x, j,
663 . ilvout, ivolu(1) )
664 IF (info == 1) THEN
665 DEALLOCATE(ipoly, rpoly, ielnod)
666 GOTO 200
667 ENDIF
668C
669 DO n=1,nnp
670 IF (ipoly(1,n) == -1) cycle
671 npoly2=npoly2+1
672 IF (.NOT.ASSOCIATED(first2)) THEN
673 ALLOCATE(first2)
674 ptr_cur2 => first2
675 ELSE
676 ALLOCATE(ptr_cur2)
677 ptr_prec2%PTR => ptr_cur2
678 ENDIF
679 nn=ipoly(2,n)
680 nhol=ipoly(6+nn+1,n)
681 ALLOCATE(ptr_cur2%IPOLY(6+nn+1+nhol),
682 . ptr_cur2%IELNOD(nn),
683 . ptr_cur2%RPOLY(4+3*nn+3*nhol))
684 DO m=1,6
685 ptr_cur2%IPOLY(m)=ipoly(m,n)
686 ENDDO
687 DO m=1,4
688 ptr_cur2%RPOLY(m)=rpoly(m,n)
689 ENDDO
690 DO m=1,ipoly(2,n)
691 nns2=nns2+1
692 ptr_cur2%IPOLY(6+m)=nns2
693 mm=ielnod(m,n)
694 ptr_cur2%IELNOD(m)=facet(j,1+mm)
695 ptr_cur2%RPOLY(4+3*(m-1)+1)=rpoly(4+3*(m-1)+1,n)
696 ptr_cur2%RPOLY(4+3*(m-1)+2)=rpoly(4+3*(m-1)+2,n)
697 ptr_cur2%RPOLY(4+3*(m-1)+3)=rpoly(4+3*(m-1)+3,n)
698 ENDDO
699 ptr_cur2%IPOLY(6+nn+1)=nhol
700 DO m=1,nhol
701 ptr_cur2%IPOLY(6+nn+1+m)=ipoly(6+nn+1+m,n)
702 ptr_cur2%RPOLY(4+3*nn+3*(m-1)+1)=
703 . rpoly(4+3*nn+3*(m-1)+1,n)
704 ptr_cur2%RPOLY(4+3*nn+3*(m-1)+2)=
705 . rpoly(4+3*nn+3*(m-1)+2,n)
706 ptr_cur2%RPOLY(4+3*nn+3*(m-1)+3)=
707 . rpoly(4+3*nn+3*(m-1)+3,n)
708 ENDDO
709 ptr_prec2 => ptr_cur2
710 ENDDO
711C
712 DEALLOCATE(ipoly, rpoly, ielnod)
713 ENDIF
714 ENDDO
715C Fin de boucle sur les faces des briques
716 itagb(i)=1
717C
718 DEALLOCATE(facet, tri, normt, poly)
719 ENDDO
720C Fin de boucle sur les briques
721C Fin de creation des polygones
722C-----------------------------------------------------------------------------------
723C Concatenation des deux listes
724 ptr_cur2 => first2
725 ptr_prec => ptr_cur
726 DO i=1,npoly2
727 npoly=npoly+1
728 ALLOCATE(ptr_cur)
729 ptr_prec%PTR => ptr_cur
730 nn=ptr_cur2%IPOLY(2)
731 nhol=ptr_cur2%IPOLY(6+nn+1)
732 ALLOCATE(ptr_cur%IPOLY(6+nn+1+nhol),
733 . ptr_cur%IELNOD(nn),
734 . ptr_cur%RPOLY(4+3*nn+3*nhol))
735 DO j=1,6
736 ptr_cur%IPOLY(j)=ptr_cur2%IPOLY(j)
737 ENDDO
738 DO j=1,4
739 ptr_cur%RPOLY(j)=ptr_cur2%RPOLY(j)
740 ENDDO
741 DO j=1,nn
742 ptr_cur%IPOLY(6+j)=nns+ptr_cur2%IPOLY(6+j)
743 ptr_cur%IELNOD(j)=ptr_cur2%IELNOD(j)
744 ptr_cur%RPOLY(4+3*(j-1)+1)=ptr_cur2%RPOLY(4+3*(j-1)+1)
745 ptr_cur%RPOLY(4+3*(j-1)+2)=ptr_cur2%RPOLY(4+3*(j-1)+2)
746 ptr_cur%RPOLY(4+3*(j-1)+3)=ptr_cur2%RPOLY(4+3*(j-1)+3)
747 ENDDO
748 ptr_cur%IPOLY(6+nn+1)=nhol
749 DO j=1,nhol
750 ptr_cur%IPOLY(6+nn+1+j)=ptr_cur2%IPOLY(6+nn+1+j)
751 ptr_cur%RPOLY(4+3*nn+3*(j-1)+1)=
752 . ptr_cur2%RPOLY(4+3*nn+3*(j-1)+1)
753 ptr_cur%RPOLY(4+3*nn+3*(j-1)+2)=
754 . ptr_cur2%RPOLY(4+3*nn+3*(j-1)+2)
755 ptr_cur%RPOLY(4+3*nn+3*(j-1)+3)=
756 . ptr_cur2%RPOLY(4+3*nn+3*(j-1)+3)
757 ENDDO
758 ptr_tmp2 => ptr_cur2%PTR
759 DEALLOCATE(ptr_cur2%IPOLY, ptr_cur2%RPOLY, ptr_cur2%IELNOD)
760 DEALLOCATE(ptr_cur2)
761 ptr_cur2 => ptr_tmp2
762 ptr_prec => ptr_cur
763 ENDDO
764C
765 nns=nns+nns2
766C---------------------------------------------------
767C Decoupage des polygones selon la verticale locale
768C---------------------------------------------------
769 nbz=ivolu(65)
770 IF (nbz > 1) THEN
771 ptr_cur => first
772 DO i=1,npoly
773 ptr_cur%NNSP=0
774 ptr_cur => ptr_cur%PTR
775 ENDDO
776C
777 vx3=rvolu(35)
778 vy3=rvolu(36)
779 vz3=rvolu(37)
780 vx1=rvolu(38)
781 vy1=rvolu(39)
782 vz1=rvolu(40)
783C Repere local
784 vnorm=sqrt(vx3**2+vy3**2+vz3**2)
785 vx3=vx3/vnorm
786 vy3=vy3/vnorm
787 vz3=vz3/vnorm
788 ss=vx3*vx1+vy3*vy1+vz3*vz1
789 vx1=vx1-ss*vx3
790 vy1=vy1-ss*vy3
791 vz1=vz1-ss*vz3
792 vnorm=sqrt(vx1**2+vy1**2+vz1**2)
793 vx1=vx1/vnorm
794 vy1=vy1/vnorm
795 vz1=vz1/vnorm
796 vx2=vy3*vz1-vz3*vy1
797 vy2=vz3*vx1-vx3*vz1
798 vz2=vx3*vy1-vy3*vx1
799C
800 x0=rvolu(41)
801 y0=rvolu(42)
802 z0=rvolu(43)
803 lz=rvolu(53)
804 dz=two*lz/nbz
805 dxm=min(dxm,dz)
806C
807 zbmin=ep30
808 zbmax=-ep30
809 DO iel=1,nel
810 n1=elem(1,iel)
811 n2=elem(2,iel)
812 n3=elem(3,iel)
813 nn1=ibuf(n1)
814 nn2=ibuf(n2)
815 nn3=ibuf(n3)
816 x1=x(1,nn1)
817 x2=x(1,nn2)
818 x3=x(1,nn3)
819 y1=x(2,nn1)
820 y2=x(2,nn2)
821 y3=x(2,nn3)
822 z1=x(3,nn1)
823 z2=x(3,nn2)
824 z3=x(3,nn3)
825 x12=x2-x1
826 y12=y2-y1
827 z12=z2-z1
828 x13=x3-x1
829 y13=y3-y1
830 z13=z3-z1
831 zl1=(x1-x0)*vx3+(y1-y0)*vy3+(z1-z0)*vz3
832 zl2=(x2-x0)*vx3+(y2-y0)*vy3+(z2-z0)*vz3
833 zl3=(x3-x0)*vx3+(y3-y0)*vy3+(z3-z0)*vz3
834 zbmin=min(zbmin,zl1)
835 zbmin=min(zbmin,zl2)
836 zbmin=min(zbmin,zl3)
837 zbmax=max(zbmax,zl1)
838 zbmax=max(zbmax,zl2)
839 zbmax=max(zbmax,zl3)
840 ENDDO
841C Clipping des polygones
842 x0=x0-lz*vx3
843 y0=y0-lz*vy3
844 z0=z0-lz*vz3
845 zlc=-lz
846 ALLOCATE(inedge(6,npoly*nppmax), rnedge(6,npoly*nppmax))
847 nbnedge=0
848 nns3=0
849 npoly_new=npoly
850 iiz=0
851 DO i=1,nbz+1
852 IF (zlc < zbmin.OR.zlc > zbmax) THEN
853 x0=x0+dz*vx3
854 y0=y0+dz*vy3
855 z0=z0+dz*vz3
856 zlc=zlc+dz
857 cycle
858 ENDIF
859 iiz=iiz+1
860 ptr_cur => first
861 DO j=1,npoly
862 zlmin=ep30
863 zlmax=-ep30
864 DO k=1,ptr_cur%IPOLY(2)
865 xx=ptr_cur%RPOLY(4+3*(k-1)+1)
866 yy=ptr_cur%RPOLY(4+3*(k-1)+2)
867 zz=ptr_cur%RPOLY(4+3*(k-1)+3)
868 zl=(xx-x0)*vx3+(yy-y0)*vy3+(zz-z0)*vz3
869 zlmin=min(zlmin,zl)
870 zlmax=max(zlmax,zl)
871 ENDDO
872C
873 IF (zlmin*zlmax < zero) THEN
874 ity=ptr_cur%IPOLY(1)
875 nn=ptr_cur%IPOLY(2)
876 nhol=0
877 IF (ity == 2) nhol=ptr_cur%IPOLY(6+nn+1)
878 ALLOCATE(ipoly(6+2*nn+1+nhol,nn),
879 . rpoly(4+3*2*nn+3*nhol,nn),
880 . ielnod(2*nn,nn),
881 . nref(2,nn), aref(4,nn),
882 . iz(3,nn))
883 CALL fvbag_gpolcut(
884 . ipoly, rpoly, ptr_cur%IPOLY, ptr_cur%RPOLY, inedge,
885 . rnedge, nbnedge, vx3, vy3, vz3,
886 . x0, y0, z0, nref, aref,
887 . nn, nhol, iiz, iz, nns3,
888 . nnp , nnsp, ielnod, ptr_cur%IELNOD)
889C
890 IF(nnp >= 1) THEN
891 ptr_tmp => ptr_cur%PTR
892 npoly_new=npoly_new-1
893 DO n=1,nnp
894 npoly_new=npoly_new+1
895 IF (n == 1) THEN
896 DEALLOCATE(ptr_cur%IPOLY, ptr_cur%RPOLY,
897 . ptr_cur%IELNOD)
898 ELSE
899 ALLOCATE(ptr_cur)
900 ptr_prec%PTR => ptr_cur
901 ptr_cur%NNSP=0
902 ENDIF
903 nn=ipoly(2,n)
904 nhol=ipoly(6+nn+1,n)
905 ALLOCATE(ptr_cur%IPOLY(6+nn+1+nhol),
906 . ptr_cur%IELNOD(nn),
907 . ptr_cur%RPOLY(4+3*nn+3*nhol))
908 DO m=1,6
909 ptr_cur%IPOLY(m)=ipoly(m,n)
910 ENDDO
911 DO m=1,4
912 ptr_cur%RPOLY(m)=rpoly(m,n)
913 ENDDO
914 DO m=1,nn
915 mm=ipoly(6+m,n)
916 ptr_cur%IPOLY(6+m)=mm
917 ptr_cur%IELNOD(m)=ielnod(m,n)
918 ptr_cur%RPOLY(4+3*(m-1)+1)=rpoly(4+3*(m-1)+1,n)
919 ptr_cur%RPOLY(4+3*(m-1)+2)=rpoly(4+3*(m-1)+2,n)
920 ptr_cur%RPOLY(4+3*(m-1)+3)=rpoly(4+3*(m-1)+3,n)
921 ENDDO
922 nhol=ipoly(6+nn+1,n)
923 ptr_cur%IPOLY(6+nn+1)=nhol
924 DO m=1,nhol
925 ptr_cur%IPOLY(6+nn+1+m)=ipoly(6+nn+1+m,n)
926 ptr_cur%RPOLY(4+3*nn+3*(m-1)+1)=
927 . rpoly(4+3*nn+3*(m-1)+1,n)
928 ptr_cur%RPOLY(4+3*nn+3*(m-1)+2)=
929 . rpoly(4+3*nn+3*(m-1)+2,n)
930 ptr_cur%RPOLY(4+3*nn+3*(m-1)+3)=
931 . rpoly(4+3*nn+3*(m-1)+3,n)
932 ENDDO
933 ptr_cur%IZ(1)=1
934 ptr_cur%IZ(2)=iz(2,n)
935 IF (n == nnp) THEN
936 ptr_cur%NNSP=nnsp
937 ALLOCATE(ptr_cur%NREF(3,nnsp),
938 . ptr_cur%AREF(4,nnsp))
939 DO m=1,nnsp
940 nns3=nns3+1
941 ptr_cur%NREF(1,m)=nns3
942 ptr_cur%NREF(2,m)=nref(1,m)
943 ptr_cur%NREF(3,m)=nref(2,m)
944 ptr_cur%AREF(1,m)=aref(1,m)
945 ptr_cur%AREF(2,m)=aref(2,m)
946 ptr_cur%AREF(3,m)=aref(3,m)
947 ptr_cur%AREF(4,m)=aref(4,m)
948 ENDDO
949 ENDIF
950 ptr_prec => ptr_cur
951 IF (n == nnp) ptr_cur%PTR => ptr_tmp
952 ENDDO
953 ENDIF
954C
955 DEALLOCATE(ipoly, rpoly, ielnod, nref, aref, iz)
956 ELSE
957C
958 IF (zlmin == zero.AND.zlmax > zero) THEN
959 nn=ptr_cur%IPOLY(2)
960 ALLOCATE(nref(2,2*nn), aref(4,2*nn))
961 CALL horiedge(
962 . ptr_cur%IPOLY, ptr_cur%RPOLY, vx3, vy3, vz3,
963 . nbnedge, inedge, rnedge, x0, y0,
964 . z0, iiz , nns3, nref, aref,
965 . nnsp )
966C
967 IF (nnsp > 0) THEN
968 ALLOCATE(ptr_cur%NREF(3,nnsp),
969 . ptr_cur%AREF(4,nnsp))
970 ptr_cur%NNSP=nnsp
971 DO n=1,nnsp
972 nns3=nns3+1
973 ptr_cur%NREF(1,n)=nns3
974 ptr_cur%NREF(2,n)=nref(1,n)
975 ptr_cur%NREF(3,n)=nref(2,n)
976 ptr_cur%AREF(1,n)=aref(1,n)
977 ptr_cur%AREF(2,n)=aref(2,n)
978 ptr_cur%AREF(3,n)=aref(3,n)
979 ptr_cur%AREF(4,n)=aref(4,n)
980 ENDDO
981 ENDIF
982 DEALLOCATE(nref, aref)
983 ENDIF
984C
985 IF (zlmin >= zero) THEN
986 ptr_cur%IZ(1)=1
987 ptr_cur%IZ(2)=iiz+1
988 ELSEIF (iiz == 1.AND.zlmax <= zero) THEN
989 ptr_cur%IZ(1)=1
990 ptr_cur%IZ(2)=1
991 ENDIF
992 ptr_prec => ptr_cur
993 ENDIF
994 ptr_cur => ptr_cur%PTR
995 ENDDO
996 npoly=npoly_new
997 x0=x0+dz*vx3
998 y0=y0+dz*vy3
999 z0=z0+dz*vz3
1000 zlc=zlc+dz
1001 ENDDO
1002C
1003 ALLOCATE(redir(nns))
1004 ptr_cur => first
1005 nns=0
1006 DO i=1,npoly
1007 nn=ptr_cur%IPOLY(2)
1008 DO j=1,nn
1009 jj=ptr_cur%IPOLY(6+j)
1010 IF (jj > 0) THEN
1011 nns=nns+1
1012 redir(jj)=nns
1013 ENDIF
1014 ENDDO
1015 ptr_cur => ptr_cur%PTR
1016 ENDDO
1017C
1018 ptr_cur => first
1019 DO i=1,npoly
1020 nnsp=ptr_cur%NNSP
1021 IF (nnsp > 0) THEN
1022 DO j=1,nnsp
1023 n1=ptr_cur%NREF(2,j)
1024 n2=ptr_cur%NREF(3,j)
1025 IF (n1 > 0) ptr_cur%NREF(2,j)=redir(n1)
1026 IF (n2 > 0) ptr_cur%NREF(3,j)=redir(n2)
1027 ENDDO
1028 ENDIF
1029 ptr_cur => ptr_cur%PTR
1030 ENDDO
1031 DEALLOCATE(redir)
1032C Creation des nouveaux polygones horizontaux
1033 ptr_cur => first
1034 DO i=1,npoly
1035 ptr_prec => ptr_cur
1036 ptr_cur => ptr_cur%PTR
1037 ENDDO
1038C
1039 ALLOCATE(ledge(nbnedge))
1040 DO i=1,iiz
1041 DO j=1,nbric
1042 nedge=0
1043 DO k=1,nbnedge
1044 IF (inedge(6,k) /= i) cycle
1045 IF (inedge(1,k) == 1.AND.inedge(5,k) /= j) cycle
1046 IF (inedge(1,k) == 2.AND.inedge(4,k) /= j.AND.
1047 . inedge(5,k) /= j) cycle
1048 nedge=nedge+1
1049 ledge(nedge)=k
1050 ENDDO
1051 IF (nedge == 0) cycle
1052 ALLOCATE(ipoly(6+2*nedge+1+nedge,nedge),
1053 . rpoly(4+6*nedge+3*nedge,nedge),
1054 . iz(3,nedge), ielnod(nedge,nedge))
1055C
1056 CALL horipoly(inedge, rnedge, ledge, nedge, ipoly,
1057 . rpoly, iz, ielnod, nnp, vx3,
1058 . vy3, vz3, i , j , nel,
1059 . tagela )
1060C
1061 DO n=1,nnp
1062 IF (ipoly(1,n) == -1) cycle
1063 npoly=npoly+1
1064 ALLOCATE(ptr_cur)
1065 ptr_prec%PTR => ptr_cur
1066 nn=ipoly(2,n)
1067 nhol=ipoly(6+nn+1,n)
1068 ALLOCATE(ptr_cur%IPOLY(6+nn+1+nhol),
1069 . ptr_cur%RPOLY(4+3*nn+3*nhol),
1070 . ptr_cur%IELNOD(nn))
1071 ptr_cur%NNSP=0
1072 DO m=1,6
1073 ptr_cur%IPOLY(m)=ipoly(m,n)
1074 ENDDO
1075 DO m=1,4
1076 ptr_cur%RPOLY(m)=rpoly(m,n)
1077 ENDDO
1078 DO m=1,nn
1079 ptr_cur%IPOLY(6+m)=ipoly(6+m,n)
1080 ptr_cur%IELNOD(m)=ielnod(m,n)
1081 ptr_cur%RPOLY(4+3*(m-1)+1)=rpoly(4+3*(m-1)+1,n)
1082 ptr_cur%RPOLY(4+3*(m-1)+2)=rpoly(4+3*(m-1)+2,n)
1083 ptr_cur%RPOLY(4+3*(m-1)+3)=rpoly(4+3*(m-1)+3,n)
1084 ENDDO
1085 nhol=ipoly(6+nn+1,n)
1086 ptr_cur%IPOLY(6+nn+1)=nhol
1087 DO m=1,nhol
1088 ptr_cur%IPOLY(6+nn+1+m)=ipoly(6+nn+1+m,n)
1089 ptr_cur%RPOLY(4+3*nn+3*(m-1)+1)=
1090 . rpoly(4+3*nn+3*(m-1)+1,n)
1091 ptr_cur%RPOLY(4+3*nn+3*(m-1)+2)=
1092 . rpoly(4+3*nn+3*(m-1)+2,n)
1093 ptr_cur%RPOLY(4+3*nn+3*(m-1)+3)=
1094 . rpoly(4+3*nn+3*(m-1)+3,n)
1095 ENDDO
1096 DO m=1,3
1097 ptr_cur%IZ(m)=iz(m,n)
1098 ENDDO
1099 ptr_prec => ptr_cur
1100 ENDDO
1101C
1102 DEALLOCATE(ipoly, rpoly, iz, ielnod)
1103 ENDDO
1104 ENDDO
1105 DEALLOCATE(ledge)
1106 ELSE
1107 ptr_cur => first
1108 DO i=1,npoly
1109 ptr_cur%NNSP=0
1110 ptr_cur%IZ(1)=1
1111 ptr_cur%IZ(2)=1
1112 ptr_cur => ptr_cur%PTR
1113C
1114 iiz=0
1115 ENDDO
1116 ENDIF
1117C---------------------------------------------------
1118C Construction des polyedres
1119C---------------------------------------------------
1120 npolh=0
1121C
1122 IF (ilvout >=1) WRITE(istdo,'(8X,A)') 'BUILDING FINITE VOLUMES'
1123 DO i=1,nbric
1124 !IF (ILVOUT >=1) CALL PROGALE_C(I, NBRIC, 2) !dynamic screen output
1125C
1126 IF (tbric(7,i) /= 2) cycle
1127C Polygones de la brique
1128 npolb=0
1129 nppmax=0
1130 ptr_cur => first
1131 DO j=1,npoly
1132 ity=ptr_cur%IPOLY(1)
1133 IF ((ity == 1.AND.ptr_cur%IPOLY(4) == i).OR.
1134 . (ity == 2.AND.
1135 . (ptr_cur%IPOLY(3) == i.OR.ptr_cur%IPOLY(4) == i))) THEN
1136 npolb=npolb+1
1137 nppmax=max(nppmax,ptr_cur%IPOLY(2))
1138 ENDIF
1139 ptr_cur => ptr_cur%PTR
1140 ENDDO
1141 IF (npolb == 0) cycle
1142C
1143 nrpmax=4+3*nppmax
1144 ALLOCATE(polb(npolb), ipoly(6+nppmax,npolb),
1145 . rpoly(nrpmax,npolb), redir(npolb))
1146 DO inz=1,iiz+1
1147 npolb=0
1148 ityzint=0
1149 ptr_cur => first
1150 DO j=1,npoly
1151 ity=ptr_cur%IPOLY(1)
1152 ityz=ptr_cur%IZ(1)
1153 IF (((ity == 1.AND.ptr_cur%IPOLY(4) == i).OR.
1154 . (ity == 2.AND.(ptr_cur%IPOLY(3) == i.OR.
1155 . ptr_cur%IPOLY(4) == i)))
1156 . .AND.
1157 . ((ityz == 1.AND.ptr_cur%IZ(2) == inz).OR.
1158 . (ityz == 2.AND.(ptr_cur%IZ(2) == inz.OR.
1159 . ptr_cur%IZ(3) == inz)))) THEN
1160 npolb=npolb+1
1161 redir(npolb)=j
1162 nn=ptr_cur%IPOLY(2)
1163 IF(ity==1.AND.ityzint==0) THEN
1164 iel=ptr_cur%IPOLY(3) ! 1<IEL<NELA
1165 IF(tagela(iel) > nel) ityzint=1
1166 ENDIF
1167 DO k=1,6+nn
1168 ipoly(k,npolb)=ptr_cur%IPOLY(k)
1169 ENDDO
1170 DO k=1,4+3*nn
1171 rpoly(k,npolb)=ptr_cur%RPOLY(k)
1172 ENDDO
1173 polb(npolb)=npolb
1174 ENDIF
1175 ptr_cur => ptr_cur%PTR
1176 ENDDO
1177 IF (npolb == 0) cycle
1178C
1179 nphmax=npolb
1180 npolhmax=nlayer
1181C
1182 300 CONTINUE
1183 IF (ALLOCATED(polh)) DEALLOCATE(polh)
1184 ALLOCATE(polh(nphmax+2,npolhmax))
1185C
1186 IF(ityzint==0) THEN
1187 CALL polyhedr(ipoly, rpoly , polb, npolb, polh,
1188 . nnp, nrpmax , nphmax, i, dxm ,
1189 . info , npolhmax, nppmax )
1190 ELSE
1191 CALL polyhedr1(ipoly, rpoly , polb, npolb, polh,
1192 . nnp, nrpmax , nphmax, i, dxm ,
1193 . info , npolhmax, nppmax, nel , inz ,
1194 . tagela )
1195 ENDIF
1196 IF (info == 1) GOTO 300
1197C Modifications des polygones
1198 ptr_cur => first
1199 npl=1
1200 polc=redir(polb(npl))
1201 DO j=1,npoly
1202 IF (j == polc) THEN
1203 ity=ptr_cur%IPOLY(1)
1204 IF (ity == 1) THEN
1205 ptr_cur%IPOLY(5)=npolh+ipoly(5,npl)
1206 iel=ptr_cur%IPOLY(3)
1207 IF(tagela(iel) > nel) THEN
1208 ptr_cur%IPOLY(6)=npolh+ipoly(6,npl)
1209 ENDIF
1210 ELSEIF (ity == 2) THEN
1211 IF (ptr_cur%IPOLY(5) == 0) THEN
1212 ptr_cur%IPOLY(5)=npolh+ipoly(5,npl)
1213 ELSE
1214 ptr_cur%IPOLY(6)=npolh+ipoly(6,npl)
1215 ENDIF
1216 ENDIF
1217 IF (npl == npolb) GOTO 350
1218 npl=npl+1
1219 polc=redir(polb(npl))
1220 ENDIF
1221 ptr_cur => ptr_cur%PTR
1222 ENDDO
1223 350 CONTINUE
1224C
1225 DO n=1,nnp
1226 npolh=npolh+1
1227
1228 IF (.NOT.ASSOCIATED(phfirst)) THEN
1229 ALLOCATE(phfirst)
1230 ph_cur => phfirst
1231 ELSE
1232 ALLOCATE(ph_cur)
1233 ph_prec%PTR => ph_cur
1234 ENDIF
1235 nn=polh(1,n)
1236 ALLOCATE(ph_cur%POLH(2+nn))
1237 ph_cur%RANK=npolh
1238 ph_cur%POLH(1)=polh(1,n)
1239 ph_cur%POLH(2)=polh(2,n)
1240 DO m=1,nn
1241 ph_cur%POLH(2+m)=redir(polh(2+m,n))
1242 ENDDO
1243 ph_prec => ph_cur
1244 ENDDO
1245C
1246 ENDDO !INZ=1,IIZ+1
1247 DEALLOCATE(ipoly, rpoly, polb, polh, redir)
1248 ENDDO !I=1,NBRIC
1249C---------------------------------------------------
1250C Passage listes chainees -> tableaux classiques pour la suite des traitements
1251C---------------------------------------------------
1252 370 CONTINUE
1253C
1254 ptr_cur => first
1255 lenimax=0
1256 lenrmax=0
1257 nns=0
1258 nns2=0
1259 DO i=1,npoly
1260 nn=ptr_cur%IPOLY(2)
1261 nns=nns+nn
1262 nns2=nns2+ptr_cur%NNSP
1263 IF (ptr_cur%IPOLY(1) == 1) THEN
1264 lenimax=max(lenimax,6+nn+1)
1265 lenrmax=max(lenrmax,4+3*nn)
1266 ELSEIF (ptr_cur%IPOLY(1) == 2) THEN
1267 nhol=ptr_cur%IPOLY(6+nn+1)
1268 lenimax=max(lenimax,6+nn+1+nhol)
1269 lenrmax=max(lenrmax,4+3*nn+3*nhol)
1270 ENDIF
1271 ptr_cur => ptr_cur%PTR
1272 ENDDO
1273 ph_cur => phfirst
1274 nphmax=0
1275C
1276 DO i=1,npolh
1277 nn=ph_cur%POLH(1)
1278 nphmax=max(nphmax,nn)
1279 ph_cur => ph_cur%PTR
1280 ENDDO
1281C
1282 npoly_old=npoly
1283 npolh_old=npolh
1284C--------------------------------------------
1285C Briques additionelles
1286C--------------------------------------------
1287 IF (nba > 0) THEN
1288 npolh=npolh+nba
1289 DO i=1,nba
1290 itagba(i)=0
1291 ENDDO
1292C
1293 DO i=1,nb_node
1294 itagn(i)=0
1295 ENDDO
1296 DO iel=1,nel
1297 n1=ibuf(elem(1,iel))
1298 n2=ibuf(elem(2,iel))
1299 n3=ibuf(elem(3,iel))
1300 itagn(n1)=1
1301 itagn(n2)=1
1302 itagn(n3)=1
1303 ENDDO
1304C
1305 count_elem_int_old = 0
1306 count_elem_ext_old = 0
1307 count_tria_interne_old = 0
1308 DO i=1,nba
1309 ntype=tba(2,i)
1310 nfac=nface(ntype)
1311 nnp=0
1312 DO j=1,nfac
1313 ity=tfaca(2*(j-1)+1,i)
1314 IF (ity == 1 .OR. ity == -2) THEN
1315C Face interne brique/brique
1316 nv=tfaca(2*(j-1)+2,i)
1317 IF (itagba(nv) == 0) THEN
1318 IF (ity == 1) count_tria_interne_old = count_tria_interne_old + 1
1319 IF (ity == -2) count_elem_int_old = count_elem_int_old + 1
1320 lenimax=max(lenimax,6+3+1)
1321 lenrmax=max(lenrmax,6+3*3)
1322 IF (ntype==2) THEN
1323 npoly=npoly+1
1324 nns=nns+3
1325 ELSEIF (ntype==3) THEN
1326 IF(j <= 2) THEN
1327 npoly=npoly+1
1328 nns=nns+3
1329 ELSE
1330 npoly=npoly+2
1331 nns=nns+6
1332 ENDIF
1333 ELSEIF (ntype==4) THEN
1334 IF(j <= 4) THEN
1335 npoly=npoly+1
1336 nns=nns+3
1337 ELSE
1338 npoly=npoly+2
1339 nns=nns+6
1340 ENDIF
1341 ELSEIF (ntype==1) THEN
1342 npoly=npoly+2
1343 nns=nns+6
1344 ENDIF
1345 ENDIF
1346C
1347 IF (ntype==2) THEN
1348 nnp=nnp+1
1349 ELSEIF (ntype==3) THEN
1350 IF(j <= 2) THEN
1351 nnp=nnp+1
1352 ELSE
1353 nnp=nnp+2
1354 ENDIF
1355 ELSEIF (ntype==4) THEN
1356 IF(j <= 4) THEN
1357 nnp=nnp+1
1358 ELSE
1359 nnp=nnp+2
1360 ENDIF
1361 ELSEIF (ntype==1) THEN
1362 nnp=nnp+2
1363 ENDIF
1364C
1365 ELSEIF (ity == 2) THEN
1366C Face brique appuyee sur surface externe (ity == 2)
1367 count_elem_ext_old = count_elem_ext_old + 1
1368 IF (ntype==2) THEN
1369 npoly=npoly+1
1370 nns=nns+3
1371 nnp=nnp+1
1372 ELSEIF (ntype==3) THEN
1373 IF(j <= 2) THEN
1374 npoly=npoly+1
1375 nns=nns+3
1376 nnp=nnp+1
1377 ELSE
1378 npoly=npoly+2
1379 nns=nns+6
1380 nnp=nnp+2
1381 ENDIF
1382 ELSEIF (ntype==4) THEN
1383 IF(j <= 4) THEN
1384 npoly=npoly+1
1385 nns=nns+3
1386 nnp=nnp+1
1387 ELSE
1388 npoly=npoly+2
1389 nns=nns+6
1390 nnp=nnp+2
1391 ENDIF
1392 ELSEIF (ntype==1) THEN
1393 npoly=npoly+2
1394 nns=nns+6
1395 nnp=nnp+2
1396 ENDIF
1397C
1398 ELSEIF (ity == 3) THEN
1399C Face interne brique/polyedres
1400 ptr_cur => first
1401 DO k=1,npoly_old
1402 IF (ptr_cur%IPOLY(1) == 1) THEN
1403 iel=ptr_cur%IPOLY(3)
1404 IF (-tagela(iel) == i) nnp=nnp+1
1405 ENDIF
1406 ptr_cur => ptr_cur%PTR
1407 ENDDO
1408 ENDIF
1409 ENDDO !J=1,NFAC
1410 itagba(i)=1
1411 nphmax=max(nphmax,nnp)
1412C Calcul de IBSA =1 si tous les noeuds de la brique sont sur l'airbag =0 autrement
1413 ii=tba(1,i)
1414 nall=1
1415 IF (ntype==2) THEN
1416 DO j=1,4
1417 jj=ixs(1+2*(j-1)+1,ii)
1418 nall=nall*itagn(jj)
1419 ENDDO
1420 ELSEIF (ntype==3) THEN
1421 DO j=1,3
1422 jj=ixs(1+j,ii)
1423 nall=nall*itagn(jj)
1424 ENDDO
1425 DO j=4,6
1426 jj=ixs(1+j+1,ii)
1427 nall=nall*itagn(jj)
1428 ENDDO
1429 ELSEIF (ntype==4) THEN
1430 DO j=1,5
1431 jj=ixs(1+j,ii)
1432 nall=nall*itagn(jj)
1433 ENDDO
1434 ELSEIF (ntype==1) THEN
1435 DO j=1,8
1436 jj=ixs(1+j,ii)
1437 nall=nall*itagn(jj)
1438 ENDDO
1439 ENDIF
1440 ibsa(i)=nall
1441 ENDDO !I=1,NBA
1442 ENDIF
1443C
1444 ALLOCATE(ipoly_f(lenimax,npoly), rpoly_f(lenrmax,npoly),
1445 . polh_f(2+nphmax,npolh), ifvnod(nns), volu(npolh),
1446 . ifvnod2(2,nns2), rfvnod2(4,nns2), xns(3,nns),
1447 . xns2(3,nns2))
1448C
1449 npoly=npoly_old
1450 npolh=npolh_old
1451
1452 volu(1:npolh) = 0
1453 imax = 0
1454 ptr_cur => first
1455 nns=0
1456 DO i=1,npoly
1457 nn=ptr_cur%IPOLY(2)
1458 DO j=1,nn
1459 jj=ptr_cur%IPOLY(6+j)
1460 IF (jj > 0) THEN
1461 nns=nns+1
1462 xns(1,nns)=ptr_cur%RPOLY(4+3*(j-1)+1)
1463 xns(2,nns)=ptr_cur%RPOLY(4+3*(j-1)+2)
1464 xns(3,nns)=ptr_cur%RPOLY(4+3*(j-1)+3)
1465 ENDIF
1466 ENDDO
1467 ptr_cur => ptr_cur%PTR
1468 ENDDO
1469C
1470 nns=0
1471 ptr_cur => first
1472 DO i=1,npoly
1473 nn=ptr_cur%IPOLY(2)
1474 DO j=1,6
1475 ipoly_f(j,i)=ptr_cur%IPOLY(j)
1476 ENDDO
1477 DO j=1,4+3*nn
1478 rpoly_f(j,i)=ptr_cur%RPOLY(j)
1479 ENDDO
1480 DO j=1,nn
1481 jj=ptr_cur%IPOLY(6+j)
1482 IF (jj > 0) THEN
1483 nns=nns+1
1484 ipoly_f(6+j,i)=nns
1485 ifvnod(nns)=ptr_cur%IELNOD(j)
1486 ELSE
1487 ipoly_f(6+j,i)=jj
1488 ENDIF
1489 ENDDO
1490 IF (ptr_cur%IPOLY(1) == 1) THEN
1491 iel=ipoly_f(3,i)
1492 IF(tagela(iel) < 0 ) THEN
1493 ipoly_f(4,i)=ipoly_f(5,i)
1494 ipoly_f(5,i)=0
1495 ELSEIF(tagela(iel) > 0) THEN
1496 IF(tagela(iel) <= nel) THEN
1497 ipoly_f(4,i)=ipoly_f(5,i)
1498 ipoly_f(5,i)=0
1499 ENDIF
1500 ENDIF
1501 ELSEIF (ptr_cur%IPOLY(1) == 2) THEN
1502 nhol=ptr_cur%IPOLY(6+nn+1)
1503 ipoly_f(6+nn+1,i)=nhol
1504 DO j=1,nhol
1505 ipoly_f(6+nn+1+j,i)=ptr_cur%IPOLY(6+nn+1+j)
1506 rpoly_f(4+3*nn+3*(j-1)+1,i)=
1507 . ptr_cur%RPOLY(4+3*nn+3*(j-1)+1)
1508 rpoly_f(4+3*nn+3*(j-1)+2,i)=
1509 . ptr_cur%RPOLY(4+3*nn+3*(j-1)+2)
1510 rpoly_f(4+3*nn+3*(j-1)+3,i)=
1511 . ptr_cur%RPOLY(4+3*nn+3*(j-1)+3)
1512 ENDDO
1513 ENDIF
1514 nnsp=ptr_cur%NNSP
1515 IF (nnsp > 0) THEN
1516 DO j=1,nnsp
1517 jj=ptr_cur%NREF(1,j)
1518 ifvnod2(1,jj)=ptr_cur%NREF(2,j)
1519 ifvnod2(2,jj)=ptr_cur%NREF(3,j)
1520 rfvnod2(1,jj)=ptr_cur%AREF(1,j)
1521 rfvnod2(2,jj)=ptr_cur%AREF(2,j)
1522 rfvnod2(3,jj)=ptr_cur%AREF(3,j)
1523 rfvnod2(4,jj)=ptr_cur%AREF(4,j)
1524 ENDDO
1525 ENDIF
1526 ptr_tmp => ptr_cur%PTR
1527 DEALLOCATE(ptr_cur%IPOLY, ptr_cur%RPOLY, ptr_cur%IELNOD)
1528 IF (nnsp > 0) DEALLOCATE(ptr_cur%NREF, ptr_cur%AREF)
1529 DEALLOCATE(ptr_cur)
1530 ptr_cur => ptr_tmp
1531 ENDDO
1532C
1533 DO i=1,nns2
1534 n1=ifvnod2(1,i)
1535 n2=ifvnod2(2,i)
1536 idep=0
1537 IF (n1 < 0) THEN
1538 ii=-n1
1539 IF (ifvnod2(1,ii) /= n2.AND.ifvnod2(2,ii) /= n2) THEN
1540 WRITE(istdo,*) 'PROBLEM DEPENDANT NODE ',i
1541 CALL my_exit(2)
1542 ENDIF
1543 n1=ifvnod2(1,ii)
1544 n2=ifvnod2(2,ii)
1545 idep=1
1546 ELSEIF (n2 < 0) THEN
1547 ii=-n2
1548 IF (ifvnod2(1,ii) /= n1.AND.ifvnod2(2,ii) /= n1) THEN
1549 WRITE(istdo,*) 'PROBLEM DEPENDANT NODE ',i
1550 CALL my_exit(2)
1551 ENDIF
1552 n1=ifvnod2(1,ii)
1553 n2=ifvnod2(2,ii)
1554 idep=1
1555 ENDIF
1556 ifvnod2(1,i)=n1
1557 ifvnod2(2,i)=n2
1558 IF (idep == 1) THEN
1559C Recalcul des coordonnees barycentriques du point dependant
1560 ii=1
1561 val=abs(xns(1,n1)-xns(1,n2))
1562 DO j=2,3
1563 IF (abs(xns(j,n1)-xns(j,n2)) > val) THEN
1564 ii=j
1565 val=abs(xns(j,n1)-xns(j,n2))
1566 ENDIF
1567 ENDDO
1568 rfvnod2(1,i)=(rfvnod2(ii+1,i)-xns(ii,n2))/
1569 . (xns(ii,n1)-xns(ii,n2))
1570 ENDIF
1571 ENDDO
1572C
1573 ph_cur => phfirst
1574 DO i=1,npolh
1575 nn=ph_cur%POLH(1)
1576 DO j=1,2+nn
1577 polh_f(j,i)=ph_cur%POLH(j)
1578 ENDDO
1579 ph_tmp => ph_cur%PTR
1580 DEALLOCATE(ph_cur%POLH)
1581 DEALLOCATE(ph_cur)
1582 ph_cur => ph_tmp
1583 ENDDO
1584C
1585 ALLOCATE(iflag(nb_node))
1586 iflag(1:nb_node) = 0
1587
1588 IF (nba > 0) THEN
1589 DO i=1,nba
1590 itagba(i)=0
1591 ENDDO
1592C
1593 count_elem_int = 0
1594 count_elem_ext = 0
1595 count_tria_interne = 0
1596 npoly_old=npoly
1597 DO i=1,nba
1598 ii=tba(1,i)
1599 ntype=tba(2,i)
1600 nfac=nface(ntype)
1601 nnp=0
1602 DO j=1,nfac
1603 ity=tfaca(2*(j-1)+1,i)
1604 IF (ity == 1) THEN
1605C Face interne brique/brique
1606 nv=tfaca(2*(j-1)+2,i)
1607 IF (itagba(nv) == 0) THEN
1608 count_tria_interne = count_tria_interne + 1
1609 IF (ntype==2) THEN
1610 iquad=0
1611 n1=fac4(1,j)
1612 n2=fac4(2,j)
1613 n3=fac4(3,j)
1614 ELSEIF (ntype==3) THEN
1615 IF(j <= 2) THEN
1616 iquad=0
1617 n1=fac6(1,j)
1618 n2=fac6(2,j)
1619 n3=fac6(3,j)
1620 ELSE
1621 iquad=1
1622 n1=fac6(1,j)
1623 n2=fac6(2,j)
1624 n3=fac6(3,j)
1625 n4=fac6(4,j)
1626 ENDIF
1627 ELSEIF (ntype==4) THEN
1628 IF(j <= 4) THEN
1629 iquad=0
1630 n1=fac5(1,j)
1631 n2=fac5(2,j)
1632 n3=fac5(3,j)
1633 ELSE
1634 iquad=1
1635 n1=fac5(1,j)
1636 n2=fac5(2,j)
1637 n3=fac5(3,j)
1638 n4=fac5(4,j)
1639 ENDIF
1640 ELSEIF (ntype==1) THEN
1641 iquad=1
1642 n1=fac(1,j)
1643 n2=fac(2,j)
1644 n3=fac(3,j)
1645 n4=fac(4,j)
1646 ENDIF
1647C
1648 IF(iquad==0) THEN
1649C Face=triangle
1650 npoly=npoly+1
1651 ipoly_f(1,npoly)=2
1652 ipoly_f(2,npoly)=3
1653 ipoly_f(3,npoly)=0
1654 ipoly_f(4,npoly)=0
1655 ipoly_f(5,npoly)=npolh+i
1656 ipoly_f(6,npoly)=npolh+nv
1657 ipoly_f(6+1,npoly)=nns+1
1658 ipoly_f(6+2,npoly)=nns+2
1659 ipoly_f(6+3,npoly)=nns+3
1660 ipoly_f(6+3+1,npoly)=0
1661C
1662 nn1=ixs(1+n1,ii)
1663 nn2=ixs(1+n2,ii)
1664 nn3=ixs(1+n3,ii)
1665 ifvnod(nns+1)=-nn1
1666 ifvnod(nns+2)=-nn2
1667 ifvnod(nns+3)=-nn3
1668 nns=nns+3
1669 x1=x(1,nn1)
1670 y1=x(2,nn1)
1671 z1=x(3,nn1)
1672 x2=x(1,nn2)
1673 y2=x(2,nn2)
1674 z2=x(3,nn2)
1675 x3=x(1,nn3)
1676 y3=x(2,nn3)
1677 z3=x(3,nn3)
1678 CALL fvnormal(x,nn1,nn2,nn3,0,nrx,nry,nrz)
1679 area2=sqrt(nrx**2+nry**2+nrz**2)
1680 rpoly_f(2,npoly)=nrx/area2
1681 rpoly_f(3,npoly)=nry/area2
1682 rpoly_f(4,npoly)=nrz/area2
1683 rpoly_f(4+1,npoly)=x1
1684 rpoly_f(4+2,npoly)=y1
1685 rpoly_f(4+3,npoly)=z1
1686 rpoly_f(4+4,npoly)=x2
1687 rpoly_f(4+5,npoly)=y2
1688 rpoly_f(4+6,npoly)=z2
1689 rpoly_f(4+7,npoly)=x3
1690 rpoly_f(4+8,npoly)=y3
1691 rpoly_f(4+9,npoly)=z3
1692C
1693 nnp=nnp+1
1694 polh_f(2+nnp,npolh+i)=npoly
1695 ELSEIF(iquad==1) THEN
1696C Face=quadrangle
1697 nn1=ixs(1+n1,ii)
1698 nn2=ixs(1+n2,ii)
1699 nn3=ixs(1+n3,ii)
1700 nn4=ixs(1+n4,ii)
1701 x1=x(1,nn1)
1702 y1=x(2,nn1)
1703 z1=x(3,nn1)
1704 x2=x(1,nn2)
1705 y2=x(2,nn2)
1706 z2=x(3,nn2)
1707 x3=x(1,nn3)
1708 y3=x(2,nn3)
1709 z3=x(3,nn3)
1710 x4=x(1,nn4)
1711 y4=x(2,nn4)
1712 z4=x(3,nn4)
1713C
1714 npoly=npoly+1
1715 ipoly_f(1,npoly)=2
1716 ipoly_f(2,npoly)=3
1717 ipoly_f(3,npoly)=0
1718 ipoly_f(4,npoly)=0
1719 ipoly_f(5,npoly)=npolh+i
1720 ipoly_f(6,npoly)=npolh+nv
1721 ipoly_f(6+1,npoly)=nns+1
1722 ipoly_f(6+2,npoly)=nns+2
1723 ipoly_f(6+3,npoly)=nns+3
1724 ipoly_f(6+3+1,npoly)=0
1725 ifvnod(nns+1)=-nn1
1726 ifvnod(nns+2)=-nn2
1727 ifvnod(nns+3)=-nn3
1728 nns=nns+3
1729 CALL fvnormal(x,nn1,nn2,nn3,0,nrx,nry,nrz)
1730 area2=sqrt(nrx**2+nry**2+nrz**2)
1731 rpoly_f(2,npoly)=nrx/area2
1732 rpoly_f(3,npoly)=nry/area2
1733 rpoly_f(4,npoly)=nrz/area2
1734 rpoly_f(4+1,npoly)=x1
1735 rpoly_f(4+2,npoly)=y1
1736 rpoly_f(4+3,npoly)=z1
1737 rpoly_f(4+4,npoly)=x2
1738 rpoly_f(4+5,npoly)=y2
1739 rpoly_f(4+6,npoly)=z2
1740 rpoly_f(4+7,npoly)=x3
1741 rpoly_f(4+8,npoly)=y3
1742 rpoly_f(4+9,npoly)=z3
1743C
1744 nnp=nnp+1
1745 polh_f(2+nnp,npolh+i)=npoly
1746C
1747 npoly=npoly+1
1748 ipoly_f(1,npoly)=2
1749 ipoly_f(2,npoly)=3
1750 ipoly_f(3,npoly)=0
1751 ipoly_f(4,npoly)=0
1752 ipoly_f(5,npoly)=npolh+i
1753 ipoly_f(6,npoly)=npolh+nv
1754 ipoly_f(6+1,npoly)=nns+1
1755 ipoly_f(6+2,npoly)=nns+2
1756 ipoly_f(6+3,npoly)=nns+3
1757 ipoly_f(6+3+1,npoly)=0
1758 ifvnod(nns+1)=-nn1
1759 ifvnod(nns+2)=-nn3
1760 ifvnod(nns+3)=-nn4
1761 nns=nns+3
1762 CALL fvnormal(x,nn1,nn3,nn4,0,nrx,nry,nrz)
1763 area2=sqrt(nrx**2+nry**2+nrz**2)
1764 rpoly_f(2,npoly)=nrx/area2
1765 rpoly_f(3,npoly)=nry/area2
1766 rpoly_f(4,npoly)=nrz/area2
1767 rpoly_f(4+1,npoly)=x1
1768 rpoly_f(4+2,npoly)=y1
1769 rpoly_f(4+3,npoly)=z1
1770 rpoly_f(4+4,npoly)=x3
1771 rpoly_f(4+5,npoly)=y3
1772 rpoly_f(4+6,npoly)=z3
1773 rpoly_f(4+7,npoly)=x4
1774 rpoly_f(4+8,npoly)=y4
1775 rpoly_f(4+9,npoly)=z4
1776C
1777 nnp=nnp+1
1778 polh_f(2+nnp,npolh+i)=npoly
1779 ENDIF
1780 ELSE
1781 DO k=1,polh_f(1,npolh+nv)
1782 kk=polh_f(2+k,npolh+nv)
1783 IF (ipoly_f(1,kk) == 2.AND.
1784 . ipoly_f(6,kk) == npolh+i) THEN
1785 nnp=nnp+1
1786 polh_f(2+nnp,npolh+i)=kk
1787 ENDIF
1788 ENDDO
1789 ENDIF
1790 ELSEIF (ity == 2) THEN
1791C Face brique appuyee sur surface externe (traitee plus loin)
1792 count_elem_ext = count_elem_ext + 1
1793 ELSEIF (ity == -2) THEN
1794C Brick face lying on an internal surface element.
1795C Need to compute the normal in order to ensure proper orientation
1796 count_elem_int = count_elem_int + 1
1797 IF (ntype==2) THEN
1798 iquad=0
1799 n1=fac4(1,j)
1800 n2=fac4(2,j)
1801 n3=fac4(3,j)
1802 ELSEIF (ntype==3) THEN
1803 IF(j <= 2) THEN
1804 iquad=0
1805 n1=fac6(1,j)
1806 n2=fac6(2,j)
1807 n3=fac6(3,j)
1808 ELSE
1809 iquad=1
1810 n1=fac6(1,j)
1811 n2=fac6(2,j)
1812 n3=fac6(3,j)
1813 n4=fac6(4,j)
1814 ENDIF
1815 ELSEIF (ntype==4) THEN
1816 IF(j <= 4) THEN
1817 iquad=0
1818 n1=fac5(1,j)
1819 n2=fac5(2,j)
1820 n3=fac5(3,j)
1821 ELSE
1822 iquad=1
1823 n1=fac5(1,j)
1824 n2=fac5(2,j)
1825 n3=fac5(3,j)
1826 n4=fac5(4,j)
1827 ENDIF
1828 ELSEIF (ntype==1) THEN
1829 iquad=1
1830 n1=fac(1,j)
1831 n2=fac(2,j)
1832 n3=fac(3,j)
1833 n4=fac(4,j)
1834 ENDIF
1835 IF (iquad == 0) THEN
1836 nn1=ixs(1+n1,ii)
1837 nn2=ixs(1+n2,ii)
1838 nn3=ixs(1+n3,ii)
1839 CALL fvnormal(x,nn1,nn2,nn3,0,nrx,nry,nrz)
1840 iflag(nn1) = 1
1841 iflag(nn2) = 1
1842 iflag(nn3) = 1
1843 DO iel = nel + 1, nel + neli
1844 IF (iflag(ibuf(elem(1, iel))) * iflag(ibuf(elem(2, iel))) *
1845 . iflag(ibuf(elem(3, iel))) == 1) THEN
1846 IF (norm(1, iel) * nrx + norm(2, iel) * nry + norm(3, iel) * nrz < 0) THEN
1847 IF (tagels(2*iel - nel - 1) == i) THEN
1848 norm(1, iel) = -norm(1, iel)
1849 norm(2, iel) = -norm(2, iel)
1850 norm(3, iel) = -norm(3, iel)
1851 ENDIF
1852 ENDIF
1853 ENDIF
1854 ENDDO
1855 iflag(nn1) = 0
1856 iflag(nn2) = 0
1857 iflag(nn3) = 0
1858 ELSE IF (iquad == 1) THEN
1859 nn1=ixs(1+n1,ii)
1860 nn2=ixs(1+n2,ii)
1861 nn3=ixs(1+n3,ii)
1862 nn4=ixs(1+n4,ii)
1863 CALL fvnormal(x,nn1,nn2,nn3,0,nrx,nry,nrz)
1864 iflag(nn1) = 1
1865 iflag(nn2) = 1
1866 iflag(nn3) = 1
1867 iflag(nn4) = 1
1868 DO iel = nel + 1, nel + neli
1869 IF (iflag(ibuf(elem(1, iel))) * iflag(ibuf(elem(2, iel))) *
1870 . iflag(ibuf(elem(3, iel))) == 1) THEN
1871 IF (norm(1, iel) * nrx + norm(2, iel) * nry + norm(3, iel) * nrz < 0) THEN
1872 IF (tagels(2*iel - nel - 1) == i) THEN
1873 norm(1, iel) = -norm(1, iel)
1874 norm(2, iel) = -norm(2, iel)
1875 norm(3, iel) = -norm(3, iel)
1876 ENDIF
1877 ENDIF
1878 ENDIF
1879 ENDDO
1880 iflag(nn1) = 0
1881 iflag(nn2) = 0
1882 iflag(nn3) = 0
1883 iflag(nn4) = 0
1884 ENDIF
1885 ELSEIF (ity == 3) THEN
1886C Face interne brique/polyedres
1887 DO k=1,npoly_old
1888 IF (ipoly_f(1,k) == 1) THEN
1889 iel=ipoly_f(3,k)
1890 IF (-tagela(iel) == i) THEN
1891 nnp=nnp+1
1892 polh_f(2+nnp,npolh+i)=k
1893 ipoly_f(1,k)=2
1894 ipoly_f(3,k)=0
1895 ip=ipoly_f(4,k)
1896 ipoly_f(4,k)=0
1897 ipoly_f(5,k)=ip
1898 ipoly_f(6,k)=npolh+i
1899 nn=ipoly_f(2,k)
1900 ipoly_f(6+nn+1,k)=0
1901 ENDIF
1902 ENDIF
1903 ENDDO
1904 ENDIF
1905 ENDDO
1906 polh_f(1,npolh+i)=nnp
1907 polh_f(2,npolh+i)=-i
1908 itagba(i)=1
1909 ENDDO ! I=1,NBA
1910C
1911 npolh=npolh+nba
1912C
1913 DO i=1,npoly_old
1914 IF (ipoly_f(1,i) == 1) THEN
1915 iel=ipoly_f(3,i)
1916 IF (tagela(iel) > 0) ipoly_f(3,i)=tagela(iel)
1917 ENDIF
1918 ENDDO
1919C External airbag surface
1920 DO iel=1,nel
1921 IF (tagels(iel) > 0) THEN
1922 npoly=npoly+1
1923 ipoly_f(1,npoly)=1
1924 ipoly_f(2,npoly)=3
1925 ipoly_f(3,npoly)=iel
1926 ipoly_f(4,npoly)=npolh_old+tagels(iel)
1927 ipoly_f(5,npoly)=0
1928 ipoly_f(6,npoly)=0
1929 ipoly_f(7,npoly)=nns+1
1930 ipoly_f(8,npoly)=nns+2
1931 ipoly_f(9,npoly)=nns+3
1932 n1=elem(1,iel)
1933 n2=elem(2,iel)
1934 n3=elem(3,iel)
1935 nn1=ibuf(n1)
1936 nn2=ibuf(n2)
1937 nn3=ibuf(n3)
1938 ifvnod(nns+1)=-nn1
1939 ifvnod(nns+2)=-nn2
1940 ifvnod(nns+3)=-nn3
1941 nns=nns+3
1942 rpoly_f(1,npoly)=elarea(iel)
1943 rpoly_f(2,npoly)=norm(1,iel)
1944 rpoly_f(3,npoly)=norm(2,iel)
1945 rpoly_f(4,npoly)=norm(3,iel)
1946 x1=x(1,nn1)
1947 y1=x(2,nn1)
1948 z1=x(3,nn1)
1949 x2=x(1,nn2)
1950 y2=x(2,nn2)
1951 z2=x(3,nn2)
1952 x3=x(1,nn3)
1953 y3=x(2,nn3)
1954 z3=x(3,nn3)
1955 rpoly_f(5,npoly)=x1
1956 rpoly_f(6,npoly)=y1
1957 rpoly_f(7,npoly)=z1
1958 rpoly_f(8,npoly)=x2
1959 rpoly_f(9,npoly)=y2
1960 rpoly_f(10,npoly)=z2
1961 rpoly_f(11,npoly)=x3
1962 rpoly_f(12,npoly)=y3
1963 rpoly_f(13,npoly)=z3
1964C
1965 nnp=polh_f(1,npolh_old+tagels(iel))
1966 nnp=nnp+1
1967 polh_f(1,npolh_old+tagels(iel))=nnp
1968 polh_f(2+nnp,npolh_old+tagels(iel))=npoly
1969 ENDIF
1970 ENDDO
1971C Internal airbag surface
1972 DO iel=nel + 1,nel + neli
1973 IF (tagels(2*iel-nel-1) > 0) THEN
1974 npoly=npoly+1
1975 ipoly_f(1,npoly)=3
1976 ipoly_f(2,npoly)=3
1977 ipoly_f(3,npoly)=iel
1978 ipoly_f(4,npoly)=npolh_old + tagels(2 * iel - nel - 1)
1979 ipoly_f(5,npoly)=npolh_old + tagels(2*iel-nel-1)
1980 ipoly_f(6,npoly)=npolh_old + tagels(2*iel-nel)
1981 ipoly_f(7,npoly)=nns+1
1982 ipoly_f(8,npoly)=nns+2
1983 ipoly_f(9,npoly)=nns+3
1984 n1=elem(1,iel)
1985 n2=elem(2,iel)
1986 n3=elem(3,iel)
1987 nn1=ibuf(n1)
1988 nn2=ibuf(n2)
1989 nn3=ibuf(n3)
1990 ifvnod(nns+1)=-nn1
1991 ifvnod(nns+2)=-nn2
1992 ifvnod(nns+3)=-nn3
1993 nns=nns+3
1994 rpoly_f(1,npoly)=elarea(iel)
1995 rpoly_f(2,npoly)=norm(1,iel)
1996 rpoly_f(3,npoly)=norm(2,iel)
1997 rpoly_f(4,npoly)=norm(3,iel)
1998 x1=x(1,nn1)
1999 y1=x(2,nn1)
2000 z1=x(3,nn1)
2001 x2=x(1,nn2)
2002 y2=x(2,nn2)
2003 z2=x(3,nn2)
2004 x3=x(1,nn3)
2005 y3=x(2,nn3)
2006 z3=x(3,nn3)
2007 rpoly_f(5,npoly)=x1
2008 rpoly_f(6,npoly)=y1
2009 rpoly_f(7,npoly)=z1
2010 rpoly_f(8,npoly)=x2
2011 rpoly_f(9,npoly)=y2
2012 rpoly_f(10,npoly)=z2
2013 rpoly_f(11,npoly)=x3
2014 rpoly_f(12,npoly)=y3
2015 rpoly_f(13,npoly)=z3
2016C Internal surface elements have two neighbours
2017 nnp=polh_f(1,npolh_old+tagels(2*iel-nel-1))
2018 nnp=nnp+1
2019 polh_f(1,npolh_old+tagels(2*iel-nel-1))=nnp
2020 polh_f(2+nnp,npolh_old+tagels(2*iel-nel-1))=npoly
2021 nnp=polh_f(1,npolh_old+tagels(2*iel-nel))
2022 nnp=nnp+1
2023 polh_f(1,npolh_old+tagels(2*iel-nel))=nnp
2024 polh_f(2+nnp,npolh_old+tagels(2*iel-nel))=npoly
2025 ENDIF
2026 ENDDO
2027 ENDIF
2028C---------------------------------------------------
2029C Surface des polygones
2030C---------------------------------------------------
2031 DO i=1,npoly
2032 ity=ipoly_f(1,i)
2033 nn=ipoly_f(2,i)
2034 nhol=0
2035 IF (ity == 2) nhol=ipoly_f(6+nn+1,i)
2036 area=zero
2037 nx=rpoly_f(2,i)
2038 ny=rpoly_f(3,i)
2039 nz=rpoly_f(4,i)
2040 IF (nhol == 0) THEN
2041 x1=rpoly_f(5,i)
2042 y1=rpoly_f(6,i)
2043 z1=rpoly_f(7,i)
2044 DO j=1,nn-2
2045 jj=j+1
2046 jjj=jj+1
2047 x2=rpoly_f(4+3*(jj-1)+1,i)
2048 y2=rpoly_f(4+3*(jj-1)+2,i)
2049 z2=rpoly_f(4+3*(jj-1)+3,i)
2050 x3=rpoly_f(4+3*(jjj-1)+1,i)
2051 y3=rpoly_f(4+3*(jjj-1)+2,i)
2052 z3=rpoly_f(4+3*(jjj-1)+3,i)
2053 x12=x2-x1
2054 y12=y2-y1
2055 z12=z2-z1
2056 x13=x3-x1
2057 y13=y3-y1
2058 z13=z3-z1
2059 nnx=y12*z13-z12*y13
2060 nny=z12*x13-x12*z13
2061 nnz=x12*y13-y12*x13
2062 area=area+(nnx*nx+nny*ny+nnz*nz)
2063 ENDDO
2064 ELSE
2065 ALLOCATE(adr(nhol+1))
2066 DO j=1,nhol
2067 adr(j)=ipoly_f(6+nn+1+j,i)
2068 ENDDO
2069 adr(nhol+1)=nn
2070 x1=rpoly_f(5,i)
2071 y1=rpoly_f(6,i)
2072 z1=rpoly_f(7,i)
2073 DO j=1,adr(1)-2
2074 jj=j+1
2075 jjj=jj+1
2076 x2=rpoly_f(4+3*(jj-1)+1,i)
2077 y2=rpoly_f(4+3*(jj-1)+2,i)
2078 z2=rpoly_f(4+3*(jj-1)+3,i)
2079 x3=rpoly_f(4+3*(jjj-1)+1,i)
2080 y3=rpoly_f(4+3*(jjj-1)+2,i)
2081 z3=rpoly_f(4+3*(jjj-1)+3,i)
2082 x12=x2-x1
2083 y12=y2-y1
2084 z12=z2-z1
2085 x13=x3-x1
2086 y13=y3-y1
2087 z13=z3-z1
2088 nnx=y12*z13-z12*y13
2089 nny=z12*x13-x12*z13
2090 nnz=x12*y13-y12*x13
2091 area=area+(nnx*nx+nny*ny+nnz*nz)
2092 ENDDO
2093 area=abs(area)
2094C On retranche la surface des trous
2095 DO j=1,nhol
2096 x1=rpoly_f(4+3*adr(j)+1,i)
2097 y1=rpoly_f(4+3*adr(j)+2,i)
2098 z1=rpoly_f(4+3*adr(j)+3,i)
2099 area2=zero
2100 DO k=adr(j)+1,adr(j+1)-2
2101 kk=k+1
2102 kkk=kk+1
2103 x2=rpoly_f(4+3*(kk-1)+1,i)
2104 y2=rpoly_f(4+3*(kk-1)+2,i)
2105 z2=rpoly_f(4+3*(kk-1)+3,i)
2106 x3=rpoly_f(4+3*(kkk-1)+1,i)
2107 y3=rpoly_f(4+3*(kkk-1)+2,i)
2108 z3=rpoly_f(4+3*(kkk-1)+3,i)
2109 x12=x2-x1
2110 y12=y2-y1
2111 z12=z2-z1
2112 x13=x3-x1
2113 y13=y3-y1
2114 z13=z3-z1
2115 nnx=y12*z13-z12*y13
2116 nny=z12*x13-x12*z13
2117 nnz=x12*y13-y12*x13
2118 area2=area2+(nnx*nx+nny*ny+nnz*nz)
2119 ENDDO
2120 area=area-abs(area2)
2121 ENDDO
2122 DEALLOCATE(adr)
2123 ENDIF
2124 rpoly_f(1,i)=half*abs(area)
2125 ENDDO
2126C---------------------------------------------------
2127C Volume des polyedres
2128C---------------------------------------------------
2129 DO i=1,npolh
2130 volu(i)=zero
2131 DO j=1,polh_f(1,i)
2132 jj=polh_f(2+j,i)
2133 area=rpoly_f(1,jj)
2134 ity=ipoly_f(1,jj)
2135 IF (ity == 1) THEN
2136 iel=ipoly_f(3,jj)
2137 IF(tagels(iel) > 0) THEN
2138 nx=rpoly_f(2,jj)
2139 ny=rpoly_f(3,jj)
2140 nz=rpoly_f(4,jj)
2141 ELSEIF(tagels(iel) == 0) THEN
2142 IF(iel <= nel) THEN
2143 nx=rpoly_f(2,jj)
2144 ny=rpoly_f(3,jj)
2145 nz=rpoly_f(4,jj)
2146 ELSEIF(iel > nel) THEN
2147 IF (ipoly_f(5,jj) == i) THEN
2148 nx=rpoly_f(2,jj)
2149 ny=rpoly_f(3,jj)
2150 nz=rpoly_f(4,jj)
2151 ELSEIF (ipoly_f(6,jj) == i) THEN
2152 nx=-rpoly_f(2,jj)
2153 ny=-rpoly_f(3,jj)
2154 nz=-rpoly_f(4,jj)
2155 ENDIF
2156 ENDIF
2157 ENDIF
2158 ELSEIF (ity == 2.OR.ity == 3) THEN
2159 IF (ipoly_f(5,jj) == i) THEN
2160 nx=rpoly_f(2,jj)
2161 ny=rpoly_f(3,jj)
2162 nz=rpoly_f(4,jj)
2163 ELSEIF (ipoly_f(6,jj) == i) THEN
2164 nx=-rpoly_f(2,jj)
2165 ny=-rpoly_f(3,jj)
2166 nz=-rpoly_f(4,jj)
2167 ENDIF
2168 ENDIF
2169 x1=rpoly_f(5,jj)
2170 y1=rpoly_f(6,jj)
2171 z1=rpoly_f(7,jj)
2172 volu(i)=volu(i)+third*area*(x1*nx+y1*ny+z1*nz)
2173 ENDDO
2174 ENDDO
2175C---------------------------------------------------
2176C Check negative volumes
2177C---------------------------------------------------
2178 IF (ityp == 8) THEN
2179 DO i = 1, npolh
2180 IF (volu(i) <= zero) THEN
2181 CALL ancmsg(msgid = 1627, msgtype = msgerror,
2182 . anmode = aninfo,
2183 . i1 = id,
2184 . c1 = titr,
2185 . i2 = ixs(11, tba(1, i)))
2186 ENDIF
2187 ENDDO
2188 ENDIF
2189C---------------------------------------------------
2190C Elimination des segments de longueur nulle dans les polygones
2191C---------------------------------------------------
2192 nns_old=nns
2193 nns_old_save=nns_old
2194 ALLOCATE(ifvnod_old(nns_old), redir(nns_old))
2195 DO i=1,nns_old
2196 ifvnod_old(i)=ifvnod(i)
2197 ENDDO
2198 nns_old=0
2199 nns=0
2200C
2201C TOLE=RVOLU(50)
2202 tole=epsilon(zero)*0.5
2203 DO i=1,npoly
2204 nnp_old=ipoly_f(2,i)
2205C
2206 lmax2=zero
2207 x0=rpoly_f(5,i)
2208 y0=rpoly_f(6,i)
2209 z0=rpoly_f(7,i)
2210 DO j=2,nnp_old
2211 x1=rpoly_f(4+3*(j-1)+1,i)
2212 y1=rpoly_f(4+3*(j-1)+2,i)
2213 z1=rpoly_f(4+3*(j-1)+3,i)
2214 lmax2=max(lmax2,(x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
2215 x0=x1
2216 y0=y1
2217 z0=z1
2218 ENDDO
2219 x1=rpoly_f(5,i)
2220 y1=rpoly_f(6,i)
2221 z1=rpoly_f(7,i)
2222 lmax2=max(lmax2,(x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
2223 tole2=tole**2*lmax2
2224C
2225 nhol=0
2226 IF (ipoly_f(1,i) == 2) nhol=ipoly_f(6+nnp_old+1,i)
2227 ALLOCATE(ipoly_old(nnp_old), rpoly_old(3*nnp_old),
2228 . adr_old(nhol+2), adr(nhol+2))
2229 DO j=1,nnp_old
2230 ipoly_old(j)=ipoly_f(6+j,i)
2231 ENDDO
2232 DO j=1,3*nnp_old
2233 rpoly_old(j)=rpoly_f(4+j,i)
2234 ENDDO
2235 nnp=0
2236 IF (nhol == 0) THEN
2237 x0=rpoly_f(5,i)
2238 y0=rpoly_f(6,i)
2239 z0=rpoly_f(7,i)
2240 IF (ipoly_old(1) > 0) THEN
2241 nns_old=nns_old+1
2242 IF (nns_old > nns_old_save) THEN
2243 CALL ancmsg(msgid=744,
2244 . msgtype=msgerror,
2245 . anmode=aninfo,
2246 . i1=id,
2247 . c1=titr)
2248 END IF
2249 nns=nns+1
2250 redir(nns_old)=nns
2251 ipoly_f(7,i)=nns
2252 ifvnod(nns)=ifvnod_old(nns_old)
2253 ELSE
2254 ipoly_f(7,i)=ipoly_old(1)
2255 ENDIF
2256 nnp=nnp+1
2257 j0=1
2258 nns1=nns
2259 DO j=2,nnp_old
2260 IF (ipoly_old(j) > 0) THEN
2261 nns_old=nns_old+1
2262 IF (nns_old > nns_old_save) THEN
2263 CALL ancmsg(msgid=744,
2264 . msgtype=msgerror,
2265 . anmode=aninfo,
2266 . i1=id,
2267 . c1=titr)
2268 END IF
2269 END IF
2270 x1=rpoly_old(3*(j-1)+1)
2271 y1=rpoly_old(3*(j-1)+2)
2272 z1=rpoly_old(3*(j-1)+3)
2273 dd=(x1-x0)**2+(y1-y0)**2+(z1-z0)**2
2274 IF (dd > tole2) THEN
2275 nnp=nnp+1
2276 IF (ipoly_old(j) > 0) THEN
2277 nns=nns+1
2278 ifvnod(nns)=ifvnod_old(nns_old)
2279 rpoly_f(4+3*(nnp-1)+1,i)=x1
2280 rpoly_f(4+3*(nnp-1)+2,i)=y1
2281 rpoly_f(4+3*(nnp-1)+3,i)=z1
2282 ipoly_f(6+nnp,i)=nns
2283 ELSE
2284 rpoly_f(4+3*(nnp-1)+1,i)=x1
2285 rpoly_f(4+3*(nnp-1)+2,i)=y1
2286 rpoly_f(4+3*(nnp-1)+3,i)=z1
2287 ipoly_f(6+nnp,i)=ipoly_old(j)
2288 ENDIF
2289 x0=x1
2290 y0=y1
2291 z0=z1
2292 ELSEIF (ipoly_old(j) > 0.AND.
2293 . ipoly_old(j0) < 0) THEN
2294 nns=nns+1
2295 ifvnod(nns)=ifvnod_old(nns_old)
2296 rpoly_f(4+3*(nnp-1)+1,i)=x1
2297 rpoly_f(4+3*(nnp-1)+2,i)=y1
2298 rpoly_f(4+3*(nnp-1)+3,i)=z1
2299 ipoly_f(6+nnp,i)=nns
2300 ENDIF
2301 IF (ipoly_old(j) > 0) redir(nns_old)=nns
2302 j0=j
2303 ENDDO
2304 x1=rpoly_f(5,i)
2305 y1=rpoly_f(6,i)
2306 z1=rpoly_f(7,i)
2307 dd=(x1-x0)**2+(y1-y0)**2+(z1-z0)**2
2308 IF (dd <= tole2.AND.ipoly_old(1) > 0) THEN
2309 redir(nns_old)=nns1
2310 nns=nns-1
2311 nnp=nnp-1
2312 ELSEIF (dd <= tole2.AND.ipoly_old(1) < 0) THEN
2313 nnp=nnp-1
2314 rpoly_f(5,i)=x0
2315 rpoly_f(6,i)=y0
2316 rpoly_f(7,i)=z0
2317 ipoly_f(7,i)=nns
2318 ENDIF
2319 ipoly_f(6+nnp+1,i)=0
2320 ELSE
2321 adr_old(1)=0
2322 adr(1)=0
2323 DO j=1,nhol
2324 adr_old(j+1)=ipoly_f(6+nnp_old+1+j,i)
2325 ENDDO
2326 adr_old(nhol+2)=nnp_old
2327 DO j=1,nhol+1
2328 x0=rpoly_f(4+3*adr(j)+1,i)
2329 y0=rpoly_f(4+3*adr(j)+2,i)
2330 z0=rpoly_f(4+3*adr(j)+3,i)
2331 IF (ipoly_old(adr_old(j)+1) > 0) THEN
2332 nns_old=nns_old+1
2333 IF (nns_old > nns_old_save) THEN
2334 CALL ancmsg(msgid=744,
2335 . msgtype=msgerror,
2336 . anmode=aninfo,
2337 . i1=id,
2338 . c1=titr)
2339 END IF
2340 nns=nns+1
2341 redir(nns_old)=nns
2342 ipoly_f(6+adr(j)+1,i)=nns
2343 ifvnod(nns)=ifvnod_old(nns_old)
2344 ELSE
2345 ipoly_f(6+adr(j)+1,i)=ipoly_old(adr_old(j)+1)
2346 ENDIF
2347 nnp=nnp+1
2348 k0=1
2349 nns1=nns
2350 DO k=adr_old(j)+2,adr_old(j+1)
2351 nns_old=nns_old+1
2352 IF (nns_old > nns_old_save) THEN
2353 CALL ancmsg(msgid=744,
2354 . msgtype=msgerror,
2355 . anmode=aninfo,
2356 . i1=id,
2357 . c1=titr)
2358 END IF
2359 x1=rpoly_old(3*(k-1)+1)
2360 y1=rpoly_old(3*(k-1)+2)
2361 z1=rpoly_old(3*(k-1)+3)
2362 dd=(x1-x0)**2+(y1-y0)**2+(z1-z0)**2
2363 IF (dd > tole2) THEN
2364 nnp=nnp+1
2365 IF (ipoly_old(k) > 0) THEN
2366 nns=nns+1
2367 ifvnod(nns)=ifvnod_old(nns_old)
2368 rpoly_f(4+3*(nnp-1)+1,i)=x1
2369 rpoly_f(4+3*(nnp-1)+2,i)=y1
2370 rpoly_f(4+3*(nnp-1)+3,i)=z1
2371 ipoly_f(6+nnp,i)=nns
2372 ELSE
2373 rpoly_f(4+3*(nnp-1)+1,i)=x1
2374 rpoly_f(4+3*(nnp-1)+2,i)=y1
2375 rpoly_f(4+3*(nnp-1)+3,i)=z1
2376 ipoly_f(6+nnp,i)=ipoly_old(k)
2377 ENDIF
2378 x0=x1
2379 y0=y1
2380 z0=z1
2381 ELSEIF (ipoly_old(k) > 0.AND.
2382 . ipoly_old(k0) < 0) THEN
2383 nns=nns+1
2384 ifvnod(nns)=ifvnod_old(nns_old)
2385 rpoly_f(4+3*(nnp-1)+1,i)=x1
2386 rpoly_f(4+3*(nnp-1)+2,i)=y1
2387 rpoly_f(4+3*(nnp-1)+3,i)=z1
2388 ipoly_f(6+nnp,i)=nns
2389 ENDIF
2390 IF (ipoly_old(k) > 0) redir(nns_old)=nns
2391 k0=k
2392 ENDDO
2393 x1=rpoly_f(4+3*(adr(j)-1)+1,i)
2394 y1=rpoly_f(4+3*(adr(j)-1)+2,i)
2395 z1=rpoly_f(4+3*(adr(j)-1)+3,i)
2396 dd=(x1-x0)**2+(y1-y0)**2+(z1-z0)**2
2397 IF (dd <= tole2.AND.ipoly_old(adr_old(j)+1) > 0) THEN
2398 redir(nns_old)=nns1
2399 nns=nns-1
2400 nnp=nnp-1
2401 ELSEIF (dd <= tole2.AND.
2402 . ipoly_old(adr_old(j)+1) < 0) THEN
2403 nnp=nnp-1
2404 rpoly_f(4+3*(adr(j)-1)+1,i)=x0
2405 rpoly_f(4+3*(adr(j)-1)+2,i)=y0
2406 rpoly_f(4+3*(adr(j)-1)+3,i)=z0
2407 ipoly_f(6+adr(j)+1,i)=nns
2408 ENDIF
2409 adr(j+1)=nnp
2410 ENDDO
2411 ipoly_f(6+nnp+1,i)=nhol
2412 DO j=1,nhol
2413 ipoly_f(6+nnp+1+j,i)=adr(j+1)
2414 rpoly_f(4+3*nnp+3*(j-1)+1,i)=
2415 . rpoly_f(4+3*nnp_old+3*(j-1)+1,i)
2416 rpoly_f(4+3*nnp+3*(j-1)+2,i)=
2417 . rpoly_f(4+3*nnp_old+3*(j-1)+2,i)
2418 rpoly_f(4+3*nnp+3*(j-1)+3,i)=
2419 . rpoly_f(4+3*nnp_old+3*(j-1)+3,i)
2420 ENDDO
2421 ENDIF
2422 ipoly_f(2,i)=nnp
2423C
2424 IF (nnp < nnp_old) THEN
2425C Nouvelle aire du polygone
2426 area=zero
2427 nx=rpoly_f(2,i)
2428 ny=rpoly_f(3,i)
2429 nz=rpoly_f(4,i)
2430 IF (nhol == 0) THEN
2431 x1=rpoly_f(5,i)
2432 y1=rpoly_f(6,i)
2433 z1=rpoly_f(7,i)
2434 DO j=1,ipoly_f(2,i)-2
2435 jj=j+1
2436 jjj=jj+1
2437 x2=rpoly_f(4+3*(jj-1)+1,i)
2438 y2=rpoly_f(4+3*(jj-1)+2,i)
2439 z2=rpoly_f(4+3*(jj-1)+3,i)
2440 x3=rpoly_f(4+3*(jjj-1)+1,i)
2441 y3=rpoly_f(4+3*(jjj-1)+2,i)
2442 z3=rpoly_f(4+3*(jjj-1)+3,i)
2443 x12=x2-x1
2444 y12=y2-y1
2445 z12=z2-z1
2446 x13=x3-x1
2447 y13=y3-y1
2448 z13=z3-z1
2449 nnx=y12*z13-z12*y13
2450 nny=z12*x13-x12*z13
2451 nnz=x12*y13-y12*x13
2452 area=area+(nnx*nx+nny*ny+nnz*nz)
2453 ENDDO
2454 ELSE
2455 x1=rpoly_f(5,i)
2456 y1=rpoly_f(6,i)
2457 z1=rpoly_f(7,i)
2458 DO j=1,adr(1)-2
2459 jj=j+1
2460 jjj=jj+1
2461 x2=rpoly_f(4+3*(jj-1)+1,i)
2462 y2=rpoly_f(4+3*(jj-1)+2,i)
2463 z2=rpoly_f(4+3*(jj-1)+3,i)
2464 x3=rpoly_f(4+3*(jjj-1)+1,i)
2465 y3=rpoly_f(4+3*(jjj-1)+2,i)
2466 z3=rpoly_f(4+3*(jjj-1)+3,i)
2467 x12=x2-x1
2468 y12=y2-y1
2469 z12=z2-z1
2470 x13=x3-x1
2471 y13=y3-y1
2472 z13=z3-z1
2473 nnx=y12*z13-z12*y13
2474 nny=z12*x13-x12*z13
2475 nnz=x12*y13-y12*x13
2476 area=area+(nnx*nx+nny*ny+nnz*nz)
2477 ENDDO
2478 ENDIF
2479C
2480 DO j=1,nhol
2481 x1=rpoly_f(4+3*adr(j+1)+1,i)
2482 y1=rpoly_f(4+3*adr(j+1)+2,i)
2483 z1=rpoly_f(4+3*adr(j+1)+3,i)
2484 DO k=adr(j+1)+1,adr(j+2)
2485 kk=k+1
2486 kkk=kk+1
2487 x2=rpoly_f(4+3*(kk-1)+1,i)
2488 y2=rpoly_f(4+3*(kk-1)+2,i)
2489 z2=rpoly_f(4+3*(kk-1)+3,i)
2490 x3=rpoly_f(4+3*(kkk-1)+1,i)
2491 y3=rpoly_f(4+3*(kkk-1)+2,i)
2492 z3=rpoly_f(4+3*(kkk-1)+3,i)
2493 x12=x2-x1
2494 y12=y2-y1
2495 z12=z2-z1
2496 x13=x3-x1
2497 y13=y3-y1
2498 z13=z3-z1
2499 nnx=y12*z13-z12*y13
2500 nny=z12*x13-x12*z13
2501 nnz=x12*y13-y12*x13
2502 area=area-(nnx*nx+nny*ny+nnz*nz)
2503 ENDDO
2504 ENDDO
2505 rpoly_f(1,i)=half*abs(area)
2506 ENDIF
2507C
2508 DEALLOCATE(ipoly_old, rpoly_old, adr_old, adr)
2509 ENDDO ! I=1,NPOLY
2510C
2511C
2512 DO i=1,nns2
2513 i1=ifvnod2(1,i)
2514 i2=ifvnod2(2,i)
2515 ifvnod2(1,i)=redir(i1)
2516 ifvnod2(2,i)=redir(i2)
2517 ENDDO
2518C
2519 DEALLOCATE(ifvnod_old, redir)
2520C---------------------------------------------------
2521C Fusion des polyedres a volume infinitesimal a leur voisin
2522C---------------------------------------------------
2523 nphmax=2*nphmax
2524 info=0
2525 ALLOCATE(volu_old(npolh), rpoly_f_old(lenrmax,npoly),
2526 . ipoly_f_old(lenimax,npoly))
2527 DO i=1,npolh
2528 volu_old(i)=volu(i)
2529 ENDDO
2530 DO i=1,npoly
2531 DO j=1,lenimax
2532 ipoly_f_old(j,i)=ipoly_f(j,i)
2533 ENDDO
2534 DO j=1,lenrmax
2535 rpoly_f_old(j,i)=rpoly_f(j,i)
2536 ENDDO
2537 ENDDO
2538 400 IF (ALLOCATED(polh_new)) DEALLOCATE(polh_new)
2539 IF (ALLOCATED(imerged)) DEALLOCATE(imerged)
2540 ALLOCATE(polh_new(2+nphmax,npolh), imerged(npolh))
2541C
2542 IF (info == 1) THEN
2543 DO i=1,npolh
2544 volu(i)=volu_old(i)
2545 ENDDO
2546 DO i=1,npoly
2547 DO j=1,lenimax
2548 ipoly_f(j,i)=ipoly_f_old(j,i)
2549 ENDDO
2550 DO j=1,lenrmax
2551 rpoly_f(j,i)=rpoly_f_old(j,i)
2552 ENDDO
2553 ENDDO
2554 ENDIF
2555C
2556 DO i=1,npolh
2557 imerged(i)=0
2558 DO j=1,2+polh_f(1,i)
2559 polh_new(j,i)=polh_f(j,i)
2560 ENDDO
2561 ENDDO
2562C
2563 DO i=1,npolh
2564 IF (polh_f(2,i) < 0) cycle ! tetraedre, ... solide additionel
2565C IF (VOLU(I) < ZERO) THEN
2566 IF (volu(i) < -em10) THEN
2567 DO j=1,polh_f(1,i)
2568 jj=polh_f(2+j,i)
2569 rpoly_f(1,jj)=-one
2570 ENDDO
2571 volu(i)=zero
2572 ENDIF
2573 ENDDO
2574C
2575C Volume moyen
2576 vm=zero
2577 npa=0
2578 DO i=1,npolh
2579 IF (volu(i) > zero) THEN
2580 vm = min(vm ,volu(i))
2581 npa=npa+1
2582 ENDIF
2583 ENDDO
2584 vm=vm/npa
2585C
2586 rvolu(33)=vm
2587C No merging will be performed during starter.
2588C Starter merging algorithm does not take internal surfaces into account and
2589C therefore may create holes through internal surfaces ...
2590 volumin = em30 * vm
2591 info=0
2592 DO i=1,npolh
2593 IF (volu(i) <= volumin) THEN
2594 areamax=zero
2595 jmax=0
2596 imax=0
2597 DO j=1,polh_new(1,i)
2598 jj=polh_new(2+j,i)
2599 ity=ipoly_f(1,jj)
2600 IF (ity == 1) cycle
2601 area=rpoly_f(1,jj)
2602 IF (area > areamax) THEN
2603 IF (ipoly_f(5,jj) == i) THEN
2604 ii=ipoly_f(6,jj)
2605 jmax=j
2606 imax=ii
2607 ELSEIF (ipoly_f(6,jj) == i) THEN
2608 ii=ipoly_f(5,jj)
2609 jmax=j
2610 imax=ii
2611 ENDIF
2612 areamax=area
2613 ENDIF
2614 ENDDO
2615C
2616 IF (jmax /= 0) THEN
2617 jjmax=polh_new(2+jmax,i)
2618 rpoly_f(1,jjmax)=-one
2619 ELSE
2620 jmax=1
2621 jjmax=polh_new(2+jmax,i)
2622 ity=ipoly_f(1,jjmax)
2623 IF (ity == 2) rpoly_f(1,jjmax)=-one
2624C On merge certain polyedres de volume < volumin avec le polyedre 1
2625C on ne devrait plus passer la
2626 imax=1
2627 ENDIF
2628 IF (imerged(imax) == 1) imax=polh_new(2,imax)
2629 np=polh_new(1,imax)
2630
2631C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
2632 IF(ilvout >= 2) THEN
2633 WRITE(iout,'(A,I8,A6,G11.4,A1,A,I8,A6,G11.4,A1)')
2634 . '** GLOBAL MERGE: MERGING FINITE VOLUME ',i,' (VOL=',
2635 . volu(i),')',' WITH FINITE VOLUME ',imax,
2636 . ' (VOL=',volu(imax),')'
2637 ENDIF
2638 volu(imax)=volu(imax)+volu(i)
2639 volu(i)=-one
2640 imerged(i)=1
2641 DO j=1,polh_new(1,i)
2642 IF (j /= jmax) THEN
2643 np=np+1
2644 IF (np > nphmax) info=1
2645 IF (info == 0) THEN
2646 jj=polh_new(2+j,i)
2647 polh_new(2+np,imax)=jj
2648 ity=ipoly_f(1,jj)
2649 iel=ipoly_f(3,jj) ! 1<IEL<NEL
2650 IF(ity == 1 .AND.iel > nel)ity=3
2651 IF (ity >= 2) THEN
2652 IF (ipoly_f(5,jj) == i) THEN
2653 ipoly_f(5,jj)=imax
2654 ELSEIF (ipoly_f(6,jj) == i) THEN
2655 ipoly_f(6,jj)=imax
2656 ENDIF
2657 ENDIF
2658 ELSE
2659 nphmax=max(nphmax,np)
2660 ENDIF
2661 ENDIF
2662 ENDDO
2663 polh_new(2,i)=imax
2664 IF (info == 0) polh_new(1,imax)=np
2665 ENDIF
2666 ENDDO !I=1,NPOLH
2667 IF (info == 1) GOTO 400
2668C
2669 vmin=ep30
2670 imin=1
2671 DO i=1,npolh
2672 IF (volu(i) <= zero) cycle
2673 IF (volu(i) < vmin) THEN
2674 vmin=volu(i)
2675 imin=i
2676 ENDIF
2677 DO j=1,polh_new(1,i)
2678 jj=polh_new(2+j,i)
2679 IF (ipoly_f(1,jj) == 1.OR.rpoly_f(1,jj) < zero) cycle
2680 IF (ipoly_f(5,jj) == ipoly_f(6,jj)) rpoly_f(1,jj)=-one
2681 ENDDO
2682 ENDDO
2683 DEALLOCATE(volu_old, rpoly_f_old, ipoly_f_old, imerged)
2684C---------------------------------------------------
2685C Verifications
2686C---------------------------------------------------
2687 volph=zero
2688 areap=zero
2689 areael=zero
2690 npolhf=0
2691 DO i=1,npolh
2692 IF (volu(i) > zero) THEN
2693 npolhf=npolhf+1
2694 volph=volph+volu(i)
2695 ENDIF
2696 ENDDO
2697 nspoly=0
2698 ncpoly=0
2699 DO i=1,npoly
2700 ity=ipoly_f(1,i)
2701 IF (ity == 1.AND.rpoly_f(1,i) > zero) THEN
2702 nspoly=nspoly+1
2703 areap=areap+rpoly_f(1,i)
2704 ELSEIF ((ity == 2 .OR. ity == 3) .AND.rpoly_f(1,i) > zero) THEN
2705 ncpoly=ncpoly+1
2706 ENDIF
2707 ENDDO
2708 DO iel=1,nel
2709 areael=areael+elarea(iel)
2710 ENDDO
2711C---------------------------------------------------
2712C Coordonnees elementaires des nouveaux noeuds
2713C---------------------------------------------------
2714 ifv=ivolu(45)
2715 ALLOCATE(fvdata(ifv)%IFVNOD(3,nns+nns2),
2716 . fvdata(ifv)%RFVNOD(2,nns+nns2))
2717C
2718 fvdata(ifv)%RFVNOD(1:2,1:nns+nns2) = 0
2719 DO i=1,nns+nns2
2720 fvdata(ifv)%IFVNOD(1,i)=0
2721 fvdata(ifv)%IFVNOD(2,i)=0
2722 fvdata(ifv)%IFVNOD(3,i)=0
2723 ENDDO
2724C
2725 nns_old=nns
2726 nns=0
2727 DO i=1,npoly
2728 nnp=ipoly_f(2,i)
2729 DO j=1,nnp
2730 IF (ipoly_f(6+j,i) > 0) THEN
2731 nns=nns+1
2732 IF (nns > nns_old) THEN
2733 CALL ancmsg(msgid=744,
2734 . msgtype=msgerror,
2735 . anmode=anstop,
2736 . i1=id,
2737 . c1=titr)
2738 END IF
2739
2740 IF (ifvnod(nns) < 0) THEN
2741C point de solide additionnel
2742 fvdata(ifv)%IFVNOD(1,nns)=2
2743 fvdata(ifv)%IFVNOD(2,nns)=-ifvnod(nns)
2744 cycle
2745 ENDIF
2746C point appuye sur airbag
2747 fvdata(ifv)%IFVNOD(1,nns)=1
2748 iel=ifvnod(nns)
2749 fvdata(ifv)%IFVNOD(2,nns)=iel
2750 xx=rpoly_f(4+3*(j-1)+1,i)
2751 yy=rpoly_f(4+3*(j-1)+2,i)
2752 zz=rpoly_f(4+3*(j-1)+3,i)
2753C
2754 n1=elema(1,iel)
2755 n2=elema(2,iel)
2756 n3=elema(3,iel)
2757 IF (tagela(iel) > 0) THEN
2758 nn1=ibuf(n1)
2759 nn2=ibuf(n2)
2760 nn3=ibuf(n3)
2761 ELSEIF (tagela(iel) < 0) THEN
2762 nn1=ibufa(n1)
2763 nn2=ibufa(n2)
2764 nn3=ibufa(n3)
2765 ENDIF
2766 x1=x(1,nn1)
2767 y1=x(2,nn1)
2768 z1=x(3,nn1)
2769 x2=x(1,nn2)
2770 y2=x(2,nn2)
2771 z2=x(3,nn2)
2772 x3=x(1,nn3)
2773 y3=x(2,nn3)
2774 z3=x(3,nn3)
2775C
2776 vpx=xx-x1
2777 vpy=yy-y1
2778 vpz=zz-z1
2779 v1x=x2-x1
2780 v1y=y2-y1
2781 v1z=z2-z1
2782 v2x=x3-x1
2783 v2y=y3-y1
2784 v2z=z3-z1
2785 CALL coorloc(vpx, vpy, vpz, v1x, v1y,
2786 . v1z, v2x, v2y, v2z, ksi,
2787 . eta)
2788C
2789 fvdata(ifv)%RFVNOD(1,nns)=ksi
2790 fvdata(ifv)%RFVNOD(2,nns)=eta
2791 ELSE
2792C point interne
2793 jj=-ipoly_f(6+j,i)
2794 fvdata(ifv)%IFVNOD(1,nns_old+jj)=3
2795 fvdata(ifv)%IFVNOD(2,nns_old+jj)=ifvnod2(1,jj)
2796 fvdata(ifv)%IFVNOD(3,nns_old+jj)=ifvnod2(2,jj)
2797 fvdata(ifv)%RFVNOD(1,nns_old+jj)=rfvnod2(1,jj)
2798 ENDIF
2799 ENDDO
2800 ENDDO
2801C
2802 DO i=1,nns2
2803 IF (fvdata(ifv)%IFVNOD(1,nns+i) == 0) THEN
2804 fvdata(ifv)%IFVNOD(1,nns+i)=3
2805 fvdata(ifv)%IFVNOD(2,nns+i)=1
2806 fvdata(ifv)%IFVNOD(3,nns+i)=2
2807 fvdata(ifv)%RFVNOD(1,nns+i)=one
2808 ENDIF
2809 ENDDO
2810C
2811 IF (nns > nns_old) THEN
2812 CALL ancmsg(msgid=744,
2813 . msgtype=msgerror,
2814 . anmode=aninfo,
2815 . i1=id,
2816 . c1=titr)
2817 END IF
2818C---------------------------------------------------
2819C Triangulation des polygones
2820C---------------------------------------------------
2821 nntr=0
2822 DO i=1,npoly
2823 nn=ipoly_f(2,i)
2824 nntr=nntr+nn
2825 ENDDO
2826C
2827 ALLOCATE(fvdata(ifv)%IFVTRI(6,nntr),
2828 . fvdata(ifv)%IFVPOLY(nntr),
2829 . fvdata(ifv)%IFVTADR(npoly+1))
2830C
2831 nntr=0
2832 npoly_old=npoly
2833 npoly=0
2834 ALLOCATE(redir_poly(npoly_old))
2835 DO i=1,npoly_old
2836 redir_poly(i)=0
2837 ENDDO
2838C
2839 nns=0
2840 DO i=1,npoly_old
2841 IF (rpoly_f(1,i) <= zero) THEN
2842 DO j=1,ipoly_f(2,i)
2843 IF (ipoly_f(6+j,i) > 0) nns=nns+1
2844 ENDDO
2845 cycle
2846 ENDIF
2847 npoly=npoly+1
2848 redir_poly(i)=npoly
2849 nnp=ipoly_f(2,i)
2850 nhol=0
2851 IF (ipoly_f(1,i) == 2) nhol=ipoly_f(6+nnp+1,i)
2852 ALLOCATE(pnodes(2,nnp), pseg(2,nnp), pholes(2,nhol),ptri(3,nnp), redir(nnp))
2853 ptri(1:3,1:nnp) = 0
2854C
2855 DO j=1,nnp
2856 IF (ipoly_f(6+j,i) > 0) THEN
2857 nns=nns+1
2858 redir(j)=nns
2859 ELSE
2860 jj=-ipoly_f(6+j,i)
2861 redir(j)=nns_old+jj
2862 ENDIF
2863 ENDDO
2864 IF (ipoly_f(1,i) == 1) THEN
2865 ipsurf=ipoly_f(3,i)
2866 IF( ipsurf > nel) THEN
2867 ipsurf=-ipsurf
2868 ic1=ipoly_f(5,i)
2869 ic2=ipoly_f(6,i)
2870 ELSE
2871 ic1=0
2872 ic2=0
2873 ENDIF
2874 ELSEIF (ipoly_f(1,i) == 2 .OR. ipoly_f(1,i) == 3) THEN
2875 ipsurf=0
2876 ic1=ipoly_f(5,i)
2877 ic2=ipoly_f(6,i)
2878 ENDIF
2879C
2880C Coordonnees des sommets dans le plan du polygone
2881 nx=rpoly_f(2,i)
2882 ny=rpoly_f(3,i)
2883 nz=rpoly_f(4,i)
2884C
2885 x0=rpoly_f(5,i)
2886 y0=rpoly_f(6,i)
2887 z0=rpoly_f(7,i)
2888 x1=rpoly_f(8,i)
2889 y1=rpoly_f(9,i)
2890 z1=rpoly_f(10,i)
2891 nrm1=sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
2892 vx1=(x1-x0)/nrm1
2893 vy1=(y1-y0)/nrm1
2894 vz1=(z1-z0)/nrm1
2895 vx2=ny*vz1-nz*vy1
2896 vy2=nz*vx1-nx*vz1
2897 vz2=nx*vy1-ny*vx1
2898C
2899 pnodes(1,1)=zero
2900 pnodes(2,1)=zero
2901 IF (nhol == 0) THEN
2902 DO j=2,nnp
2903 xx=rpoly_f(4+3*(j-1)+1,i)
2904 yy=rpoly_f(4+3*(j-1)+2,i)
2905 zz=rpoly_f(4+3*(j-1)+3,i)
2906 vx=xx-x0
2907 vy=yy-y0
2908 vz=zz-z0
2909 pnodes(1,j)=vx*vx1+vy*vy1+vz*vz1
2910 pnodes(2,j)=vx*vx2+vy*vy2+vz*vz2
2911 pseg(1,j-1)=j-1
2912 pseg(2,j-1)=j
2913 ENDDO
2914 pseg(1,nnp)=nnp
2915 pseg(2,nnp)=1
2916 nseg=nnp
2917 ELSE
2918 ALLOCATE(adr(nhol+1))
2919 DO j=1,nhol
2920 adr(j)=ipoly_f(6+nnp+1+j,i)
2921 ENDDO
2922 adr(nhol+1)=nnp
2923 nseg=0
2924 DO j=2,adr(1)
2925 xx=rpoly_f(4+3*(j-1)+1,i)
2926 yy=rpoly_f(4+3*(j-1)+2,i)
2927 zz=rpoly_f(4+3*(j-1)+3,i)
2928 vx=xx-x0
2929 vy=yy-y0
2930 vz=zz-z0
2931 pnodes(1,j)=vx*vx1+vy*vy1+vz*vz1
2932 pnodes(2,j)=vx*vx2+vy*vy2+vz*vz2
2933 nseg=nseg+1
2934 pseg(1,nseg)=j-1
2935 pseg(2,nseg)=j
2936 ENDDO
2937 nseg=nseg+1
2938 pseg(1,nseg)=adr(1)
2939 pseg(2,nseg)=1
2940 DO j=1,nhol
2941 xx=rpoly_f(4+3*adr(j)+1,i)
2942 yy=rpoly_f(4+3*adr(j)+2,i)
2943 zz=rpoly_f(4+3*adr(j)+3,i)
2944 vx=xx-x0
2945 vy=yy-y0
2946 vz=zz-z0
2947 pnodes(1,adr(j)+1)=vx*vx1+vy*vy1+vz*vz1
2948 pnodes(2,adr(j)+1)=vx*vx2+vy*vy2+vz*vz2
2949 DO k=adr(j)+2,adr(j+1)
2950 xx=rpoly_f(4+3*(k-1)+1,i)
2951 yy=rpoly_f(4+3*(k-1)+2,i)
2952 zz=rpoly_f(4+3*(k-1)+3,i)
2953 vx=xx-x0
2954 vy=yy-y0
2955 vz=zz-z0
2956 pnodes(1,k)=vx*vx1+vy*vy1+vz*vz1
2957 pnodes(2,k)=vx*vx2+vy*vy2+vz*vz2
2958 nseg=nseg+1
2959 pseg(1,nseg)=k-1
2960 pseg(2,nseg)=k
2961 ENDDO
2962 nseg=nseg+1
2963 pseg(1,nseg)=adr(j+1)
2964 pseg(2,nseg)=adr(j)+1
2965C
2966 xx=rpoly_f(4+3*nnp+3*(j-1)+1,i)
2967 yy=rpoly_f(4+3*nnp+3*(j-1)+2,i)
2968 zz=rpoly_f(4+3*nnp+3*(j-1)+3,i)
2969 vx=xx-x0
2970 vy=yy-y0
2971 vz=zz-z0
2972 pholes(1,j)=vx*vx1+vy*vy1+vz*vz1
2973 pholes(2,j)=vx*vx2+vy*vy2+vz*vz2
2974 ENDDO
2975 DEALLOCATE(adr)
2976 ENDIF
2977C
2978 nelp = 0
2979 CALL c_tricall(pnodes, pseg, pholes, ptri, nnp, nseg, nhol, nelp)
2980C
2981 fvdata(ifv)%IFVTADR(npoly)=nntr+1
2982C avoid bug with PGI => comp O1
2983C CALL BIDON1(IFV)
2984 DO j=1,nelp
2985 nntr=nntr+1
2986 n1=ptri(1,j)
2987 n2=ptri(2,j)
2988 n3=ptri(3,j)
2989C Verification de la normale
2990 x1=rpoly_f(4+3*(n1-1)+1,i)
2991 y1=rpoly_f(4+3*(n1-1)+2,i)
2992 z1=rpoly_f(4+3*(n1-1)+3,i)
2993 x2=rpoly_f(4+3*(n2-1)+1,i)
2994 y2=rpoly_f(4+3*(n2-1)+2,i)
2995 z2=rpoly_f(4+3*(n2-1)+3,i)
2996 x3=rpoly_f(4+3*(n3-1)+1,i)
2997 y3=rpoly_f(4+3*(n3-1)+2,i)
2998 z3=rpoly_f(4+3*(n3-1)+3,i)
2999 x12=x2-x1
3000 y12=y2-y1
3001 z12=z2-z1
3002 x13=x3-x1
3003 y13=y3-y1
3004 z13=z3-z1
3005 nrx=y12*z13-z12*y13
3006 nry=z12*x13-x12*z13
3007 nrz=x12*y13-y12*x13
3008 ss=nrx*nx+nry*ny+nrz*nz
3009C
3010 IF (ss > zero) THEN
3011 fvdata(ifv)%IFVTRI(1,nntr)=redir(n1)
3012 fvdata(ifv)%IFVTRI(2,nntr)=redir(n2)
3013 fvdata(ifv)%IFVTRI(3,nntr)=redir(n3)
3014 ELSE
3015 fvdata(ifv)%IFVTRI(1,nntr)=redir(n1)
3016 fvdata(ifv)%IFVTRI(2,nntr)=redir(n3)
3017 fvdata(ifv)%IFVTRI(3,nntr)=redir(n2)
3018 ENDIF
3019 fvdata(ifv)%IFVTRI(4,nntr)=ipsurf
3020 fvdata(ifv)%IFVTRI(5,nntr)=ic1
3021 fvdata(ifv)%IFVTRI(6,nntr)=ic2
3022 fvdata(ifv)%IFVPOLY(nntr)=nntr
3023 ENDDO
3024C
3025 DEALLOCATE(pnodes, pseg, pholes, ptri, redir)
3026 ENDDO
3027 fvdata(ifv)%IFVTADR(npoly+1)=nntr+1
3028C---------------------------------------------------
3029C Structure des polyedres
3030C---------------------------------------------------
3031 ALLOCATE(fvdata(ifv)%IFVPOLH(2*npoly),
3032 . fvdata(ifv)%IFVPADR(npolh+1),
3033 . fvdata(ifv)%IDPOLH(npolh),
3034 . fvdata(ifv)%IBPOLH(npolh))
3035 nnp=0
3036 npolh_old=npolh
3037 npolh=0
3038 ALLOCATE(redir_polh(npolh_old))
3039 DO i=1,npolh_old
3040 redir_polh(i)=0
3041 ENDDO
3042C
3043 DO i=1,npolh_old
3044 IF (volu(i) <= zero) cycle
3045 npolh=npolh+1
3046 redir_polh(i)=npolh
3047 fvdata(ifv)%IFVPADR(npolh)=nnp+1
3048 DO j=1,polh_new(1,i)
3049 jj=redir_poly(polh_new(2+j,i))
3050 IF (jj == 0) cycle
3051 nnp=nnp+1
3052 fvdata(ifv)%IFVPOLH(nnp)=redir_poly(polh_new(2+j,i))
3053 ENDDO
3054C
3055 fvdata(ifv)%IDPOLH(npolh)=npolh
3056C
3057 ii=polh_new(2,i)
3058 IF (ii >= 0) ii=0
3059 IF (ii < 0)THEN
3060 IF(ibsa(-ii) == 1) ii=-ii
3061 ENDIF
3062 fvdata(ifv)%IBPOLH(npolh)=ii
3063 ENDDO
3064 fvdata(ifv)%IFVPADR(npolh+1)=nnp+1
3065C
3066 DO i=1,nntr
3067 IF (fvdata(ifv)%IFVTRI(4,i) <= 0) THEN
3068 ic1=fvdata(ifv)%IFVTRI(5,i)
3069 ic2=fvdata(ifv)%IFVTRI(6,i)
3070 fvdata(ifv)%IFVTRI(5,i)=redir_polh(ic1)
3071 fvdata(ifv)%IFVTRI(6,i)=redir_polh(ic2)
3072 ENDIF
3073 ENDDO
3074C
3075 ivolu(46)=nns+nns2
3076 ivolu(47)=nntr
3077 ivolu(48)=npoly
3078 ivolu(49)=npolh
3079C
3080 fvdata(ifv)%NNS=nns+nns2
3081 fvdata(ifv)%NNTR=nntr
3082 fvdata(ifv)%LENP=fvdata(ifv)%IFVTADR(npoly+1)-1
3083 fvdata(ifv)%NPOLY=npoly
3084 fvdata(ifv)%LENH=fvdata(ifv)%IFVPADR(npolh+1)-1
3085 fvdata(ifv)%NPOLH=npolh
3086C
3087 ALLOCATE(fvdata(ifv)%MPOLH(npolh), fvdata(ifv)%QPOLH(3,npolh),
3088 . fvdata(ifv)%EPOLH(npolh), fvdata(ifv)%PPOLH(npolh),
3089 . fvdata(ifv)%RPOLH(npolh), fvdata(ifv)%GPOLH(npolh),
3090 . fvdata(ifv)%TPOLH(npolh),
3091 . fvdata(ifv)%CPAPOLH(npolh), fvdata(ifv)%CPBPOLH(npolh),
3092 . fvdata(ifv)%CPCPOLH(npolh), fvdata(ifv)%RMWPOLH(npolh),
3093 . fvdata(ifv)%VPOLH_INI(npolh), fvdata(ifv)%DTPOLH(npolh))
3094 ALLOCATE(fvdata(ifv)%CPDPOLH(npolh), fvdata(ifv)%CPEPOLH(npolh),
3095 . fvdata(ifv)%CPFPOLH(npolh))
3096C---------------------------------------------------
3097C Initialisation des quantites dans les polyedres
3098C---------------------------------------------------
3099 gamai=rvolu(1)
3100 cpai=rvolu(7)
3101 cpbi=rvolu(8)
3102 cpci=rvolu(9)
3103 rmwi=rvolu(10)
3104 pini=rvolu(12)
3105 ti=rvolu(13)
3106 cpdi=rvolu(56)
3107 cpei=rvolu(57)
3108 cpfi=rvolu(58)
3109 ti2=ti*ti
3110 rhoi=pini/(ti*rmwi)
3111C
3112 cpi=cpai+half*cpbi*ti+third*cpci*ti2
3113 IF(ivolu(2) == 8) THEN
3114 cpi=cpi+fourth*cpdi*ti2*ti-cpei/ti2+one_fifth*cpfi*ti2*ti2
3115 ENDIF
3116 cvi=cpi-rmwi
3117C
3118 tmass=0
3119 DO i=1,npolh_old
3120 ii=redir_polh(i)
3121 IF (ii == 0) cycle
3122 fvdata(ifv)%MPOLH(ii)=rhoi*volu(i)
3123 fvdata(ifv)%QPOLH(1,ii)=zero
3124 fvdata(ifv)%QPOLH(2,ii)=zero
3125 fvdata(ifv)%QPOLH(3,ii)=zero
3126 fvdata(ifv)%EPOLH(ii)=fvdata(ifv)%MPOLH(ii)*cvi*ti
3127 fvdata(ifv)%PPOLH(ii)=pini
3128 fvdata(ifv)%RPOLH(ii)=rhoi
3129 fvdata(ifv)%GPOLH(ii)=gamai
3130 fvdata(ifv)%CPAPOLH(ii)=cpai
3131 fvdata(ifv)%CPBPOLH(ii)=cpbi
3132 fvdata(ifv)%CPCPOLH(ii)=cpci
3133 fvdata(ifv)%TPOLH(ii)=ti
3134 fvdata(ifv)%CPDPOLH(ii)=cpdi
3135 fvdata(ifv)%CPEPOLH(ii)=cpei
3136 fvdata(ifv)%CPFPOLH(ii)=cpfi
3137 fvdata(ifv)%RMWPOLH(ii)=rmwi
3138 fvdata(ifv)%DTPOLH(ii)=ep30
3139 fvdata(ifv)%VPOLH_INI(ii)=volu(i)
3140 tmass=tmass+rhoi*volu(i)
3141 ENDDO
3142C
3143 imin=redir_polh(imin)
3144 DEALLOCATE(ipoly_f, rpoly_f, polh_f, volu, redir_poly, redir_polh,
3145 . polh_new, ifvnod, ifvnod2, rfvnod2, xns, xns2)
3146C---------------------------------------------------
3147C Impressions
3148C---------------------------------------------------
3149 nstr=0
3150 nctr=0
3151 DO i=1,nntr
3152 IF (fvdata(ifv)%IFVTRI(4,i) > 0) THEN
3153 nstr=nstr+1
3154 ELSE
3155 nctr=nctr+1
3156 ENDIF
3157 ENDDO
3158C
3159C NSPOLY NUMBER OF SURFACE POLYGONS
3160C NSTR NUMBER OF SURFACE TRIANGLES
3161C NCPOLY NUMBER OF COMMUNICATION POLYGONS
3162C NCTR NUMBER OF COMMUNICATION TRIANGLES
3163C NPOLHF NUMBER OF POLYHEDRA (FINITE VOLUMES)
3164C
3165 WRITE(iout,1000) ivolu(1),nspoly,nstr,ncpoly,nctr,npolhf
3166 WRITE(iout,1010) vmin,imin,volumin
3167 WRITE(iout,1020) volph,areap,tmass
3168C
3169C---------------------------------------------------
3170C Anim des volumes finis
3171C---------------------------------------------------
3172 ivolu(43)=0
3173 fvdata(ifv)%NPOLH_ANIM=npolh
3174 lenp=fvdata(ifv)%IFVTADR(npoly+1)
3175 lenh=fvdata(ifv)%IFVPADR(npolh+1)
3176 ALLOCATE(fvdata(ifv)%IFVPOLY_ANIM(lenp),
3177 . fvdata(ifv)%IFVTADR_ANIM(npoly+1),
3178 . fvdata(ifv)%IFVPOLH_ANIM(lenh),
3179 . fvdata(ifv)%IFVPADR_ANIM(npolh+1),
3180 . fvdata(ifv)%IFVTRI_ANIM(6,nntr),
3181 . fvdata(ifv)%REDIR_ANIM(nns+nns2),
3182 . fvdata(ifv)%NOD_ANIM(3,nns+nns2),
3183 . redir(nns+nns2), itagt(nntr))
3184 DO i=1,lenp
3185 fvdata(ifv)%IFVPOLY_ANIM(i)=fvdata(ifv)%IFVPOLY(i)
3186 ENDDO
3187 DO i=1,npoly+1
3188 fvdata(ifv)%IFVTADR_ANIM(i)=fvdata(ifv)%IFVTADR(i)
3189 ENDDO
3190 DO i=1,lenh
3191 fvdata(ifv)%IFVPOLH_ANIM(i)=fvdata(ifv)%IFVPOLH(i)
3192 ENDDO
3193 DO i=1,npolh+1
3194 fvdata(ifv)%IFVPADR_ANIM(i)=fvdata(ifv)%IFVPADR(i)
3195 ENDDO
3196C-----------------------------------------
3197C Elimination des noeuds doubles pour anim
3198C-----------------------------------------
3199 tole=em05*fac_length
3200 tole2=tole**2
3201 nns_anim=0
3202 ALLOCATE(pnodes(3,nns+nns2))
3203 DO i=1,nns
3204 IF (fvdata(ifv)%IFVNOD(1,i) == 1) THEN
3205 iel=fvdata(ifv)%IFVNOD(2,i)
3206 ksi=fvdata(ifv)%RFVNOD(1,i)
3207 eta=fvdata(ifv)%RFVNOD(2,i)
3208 n1=elema(1,iel)
3209 n2=elema(2,iel)
3210 n3=elema(3,iel)
3211 IF (tagela(iel) > 0) THEN
3212 n1=ibuf(n1)
3213 n2=ibuf(n2)
3214 n3=ibuf(n3)
3215 ELSEIF (tagela(iel) < 0) THEN
3216 n1=ibufa(n1)
3217 n2=ibufa(n2)
3218 n3=ibufa(n3)
3219 ENDIF
3220 x1=x(1,n1)
3221 y1=x(2,n1)
3222 z1=x(3,n1)
3223 x2=x(1,n2)
3224 y2=x(2,n2)
3225 z2=x(3,n2)
3226 x3=x(1,n3)
3227 y3=x(2,n3)
3228 z3=x(3,n3)
3229 pnodes(1,i)=(one-ksi-eta)*x1+ksi*x2+eta*x3
3230 pnodes(2,i)=(one-ksi-eta)*y1+ksi*y2+eta*y3
3231 pnodes(3,i)=(one-ksi-eta)*z1+ksi*z2+eta*z3
3232 ELSEIF (fvdata(ifv)%IFVNOD(1,i) == 2) THEN
3233 ii=fvdata(ifv)%IFVNOD(2,i)
3234 pnodes(1,i)=x(1,ii)
3235 pnodes(2,i)=x(2,ii)
3236 pnodes(3,i)=x(3,ii)
3237 ENDIF
3238 ENDDO
3239 DO i=1,nns2
3240 ii=nns+i
3241 i1=fvdata(ifv)%IFVNOD(2,ii)
3242 i2=fvdata(ifv)%IFVNOD(3,ii)
3243 alpha=fvdata(ifv)%RFVNOD(1,ii)
3244 pnodes(1,ii)=alpha*pnodes(1,i1)+(one-alpha)*pnodes(1,i2)
3245 pnodes(2,ii)=alpha*pnodes(2,i1)+(one-alpha)*pnodes(2,i2)
3246 pnodes(3,ii)=alpha*pnodes(3,i1)+(one-alpha)*pnodes(3,i2)
3247 ENDDO
3248C--------------------------
3249C Impression
3250C--------------------------
3251 IF(ilvout >=3 ) THEN
3252 ALLOCATE(areatr(nntr))
3253 DO i=1,nntr
3254 n1=fvdata(ifv)%IFVTRI(1,i)
3255 n2=fvdata(ifv)%IFVTRI(2,i)
3256 n3=fvdata(ifv)%IFVTRI(3,i)
3257 x1=pnodes(1,n1)
3258 y1=pnodes(2,n1)
3259 z1=pnodes(3,n1)
3260 x2=pnodes(1,n2)
3261 y2=pnodes(2,n2)
3262 z2=pnodes(3,n2)
3263 x3=pnodes(1,n3)
3264 y3=pnodes(2,n3)
3265 z3=pnodes(3,n3)
3266 CALL fvnormal(pnodes,n1,n2,n3,0,nrx,nry,nrz)
3267 area2=sqrt(nrx**2+nry**2+nrz**2)
3268 areatr(i)=half*area2
3269 ENDDO
3270C
3271C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
3272 WRITE(iout,'(/A,10X,3A)')' FINITE VOLUME',' BRICK VOLUME ',
3273 . ' POLYGONE TRIANGLE AREA'
3274C
3275 DO i=1,npolh
3276 n1= fvdata(ifv)%IDPOLH(i)
3277 n2= fvdata(ifv)%IBPOLH(i)
3278 n3= 0
3279 IF(n2 /= 0) THEN
3280 n3= tba(1,iabs(n2))
3281 ENDIF
3282 volph=fvdata(ifv)%VPOLH_INI(i)
3283 ii=0
3284 jjj=0
3285 kkk=0
3286 DO j=fvdata(ifv)%IFVPADR(i),fvdata(ifv)%IFVPADR(i+1)-1
3287 jjj=jjj+1
3288 jj=fvdata(ifv)%IFVPOLH(j)
3289 DO k=fvdata(ifv)%IFVTADR(jj),fvdata(ifv)%IFVTADR(jj+1)-1
3290 kkk=kkk+1
3291 kk=fvdata(ifv)%IFVPOLY(k)
3292 area=areatr(kk)
3293 iel=fvdata(ifv)%IFVTRI(4,kk)
3294 ic1=fvdata(ifv)%IFVTRI(5,kk)
3295 ic2=fvdata(ifv)%IFVTRI(6,kk)
3296 IF(kkk == 1) THEN
3297 WRITE(iout,'(3I10,2X,F10.2,4X,I5,5X,I10,5X,G14.6,3I10)')
3298 . n1,n2,n3,volph,jjj,kk,area,iel,ic1,ic2
3299 ELSE
3300 WRITE(iout,'(46X,I5,5X,I10,5X,G14.6,3I10)')
3301 . jjj,kk,area,iel,ic1,ic2
3302 ENDIF
3303 ENDDO
3304 ENDDO
3305 ENDDO
3306 DEALLOCATE(areatr)
3307 ENDIF
3308C-------------------------------------------
3309C Merge des points coincidents
3310C-------------------------------------------
3311 IF (ilvout >=1) WRITE(istdo,'(8X,A)')
3312 . 'MERGING COINCIDENT NODES FOR ANIM'
3313
3314 ALLOCATE(radius(nns + nns2), idx(nns + nns2))
3315 DO i = 1, nns + nns2
3316 xx=pnodes(1,i)
3317 yy=pnodes(2,i)
3318 zz=pnodes(3,i)
3319 radius(i) = xx*xx + yy*yy + zz*zz
3320 idx(i) = i
3321 ENDDO
3322
3323 CALL quicksort(radius, idx, 1, nns + nns2)
3324 i = 1
3325 DO WHILE(i < nns+ nns2)
3326 !IF (ILVOUT >=1) CALL PROGALE_C(I, NNS+NNS2, 3) !dynamic screen output
3327 iad1 = idx(i)
3328 xx = pnodes(1,iad1)
3329 yy = pnodes(2,iad1)
3330 zz = pnodes(3,iad1)
3331 r1 = radius(i)
3332 nns_anim = nns_anim + 1
3333 redir(iad1) = nns_anim
3334 fvdata(ifv)%REDIR_ANIM(nns_anim)=iad1
3335 fvdata(ifv)%NOD_ANIM(1,nns_anim)=xx
3336 fvdata(ifv)%NOD_ANIM(2,nns_anim)=yy
3337 fvdata(ifv)%NOD_ANIM(3,nns_anim)=zz
3338 j = i
3339 DO WHILE(j <= nns + nns2)
3340 IF (abs(r1 - radius(j)) <= tole2) THEN
3341 j = j + 1
3342 ELSE
3343 EXIT
3344 ENDIF
3345 ENDDO
3346 j = j - 1
3347
3348 iad = i + 1
3349 DO WHILE (iad <= j)
3350 iad2 = idx(iad)
3351 xn(1) = pnodes(1,iad2)
3352 xn(2) = pnodes(2,iad2)
3353 xn(3) = pnodes(3,iad2)
3354 dd2=(xx-xn(1))**2+(yy-xn(2))**2+(zz-xn(3))**2
3355 IF (dd2 <= tole2) THEN
3356C NODE IAD2 is duplicated with node IAD1
3357 redir(iad2) = redir(iad1)
3358 iad = iad + 1
3359 ELSE
3360C put it at the end
3361 itmp1 = idx(j)
3362 idx(j) = iad2
3363 idx(iad) = itmp1
3364 j = j - 1
3365 ENDIF
3366 ENDDO
3367 i = j + 1
3368 ENDDO
3369
3370 DEALLOCATE(radius, idx)
3371
3372 fvdata(ifv)%NNS_ANIM=nns_anim
3373 fvdata(ifv)%ID=ivolu(1)
3374C
3375 DO i=1,nntr
3376 n1=fvdata(ifv)%IFVTRI(1,i)
3377 n2=fvdata(ifv)%IFVTRI(2,i)
3378 n3=fvdata(ifv)%IFVTRI(3,i)
3379 fvdata(ifv)%IFVTRI_ANIM(1,i)=redir(n1)
3380 fvdata(ifv)%IFVTRI_ANIM(2,i)=redir(n2)
3381 fvdata(ifv)%IFVTRI_ANIM(3,i)=redir(n3)
3382 fvdata(ifv)%IFVTRI_ANIM(4,i)=
3383 . fvdata(ifv)%IFVTRI(4,i)
3384 fvdata(ifv)%IFVTRI_ANIM(5,i)=
3385 . fvdata(ifv)%IFVTRI(5,i)
3386 fvdata(ifv)%IFVTRI_ANIM(6,i)=
3387 . fvdata(ifv)%IFVTRI(6,i)
3388 ENDDO
3389 DEALLOCATE(pnodes)
3390C
3391 DEALLOCATE(redir, itagt)
3392 DEALLOCATE(iflag)
3393C
3394 RETURN
3395C
33961000 FORMAT(/5x,'VOLUME NUMBER ',i10,
3397 . /5x,'NUMBER OF SURFACE POLYGONS. . . . . . .=',i10,
3398 . /5x,'NUMBER OF SURFACE TRIANGLES . . . . . .=',i10,
3399 . /5x,'NUMBER OF COMMUNICATION POLYGONS. . . .=',i10,
3400 . /5x,'NUMBER OF COMMUNICATION TRIANGLES . . .=',i10,
3401 . /5x,'NUMBER OF FINITE VOLUMES. . . . . . . .=',i10)
34021010 FORMAT( 5x,'MIN FINITE VOLUME VOLUME. . . . . . . .=',1pg20.13,
3403 . 5x,'(FINITE VOLUME ID ',i10,')',
3404 . /5x,'INITIAL MERGING VOLUME. . . . . . . . .=',1pg20.13)
34051020 FORMAT( 5x,'SUM VOLUME OF FINITE VOLUMES. . . . . .=',1pg20.13,
3406 . /5x,'SUM AREA SURFACE TRIANGLES. . . . . . .=',1pg20.13,
3407 . /5x,'SUM MASS OF FINITE VOLUMES. . . . . . .=',1pg20.13)
3408C
3409 END
3410!||====================================================================
3411!|| coorloc ../starter/source/airbag/fvmesh.F
3412!||--- called by ------------------------------------------------------
3413!|| fvmesh1 ../starter/source/airbag/fvmesh.F
3414!||====================================================================
3415 SUBROUTINE coorloc(VPX, VPY, VPZ, V1X, V1Y,
3416 . V1Z, V2X, V2Y, V2Z, KSI,
3417 . ETA)
3418C-----------------------------------------------
3419C I m p l i c i t T y p e s
3420C-----------------------------------------------
3421#include "implicit_f.inc"
3422C-----------------------------------------------
3423C D u m m y A r g u m e n t s
3424C-----------------------------------------------
3425 my_real
3426 . vpx, vpy, vpz, v1x, v1y,
3427 . v1z, v2x, v2y, v2z, ksi,
3428 . eta
3429C-----------------------------------------------
3430C L o c a l V a r i a b l e s
3431C-----------------------------------------------
3432 my_real
3433 . sp1, sp2, s11, s22, s12, d
3434C
3435 sp1=vpx*v1x+vpy*v1y+vpz*v1z
3436 sp2=vpx*v2x+vpy*v2y+vpz*v2z
3437 s11=v1x*v1x+v1y*v1y+v1z*v1z
3438 s22=v2x*v2x+v2y*v2y+v2z*v2z
3439 s12=v1x*v2x+v1y*v2y+v1z*v2z
3440C
3441 d=s11*s22-s12*s12
3442C
3443 ksi=(sp1*s22-sp2*s12)/d
3444 eta=(sp2*s11-sp1*s12)/d
3445C
3446 RETURN
3447 END
3448!||====================================================================
3449!|| fvdim ../starter/source/airbag/fvmesh.F
3450!||--- called by ------------------------------------------------------
3451!|| lectur ../starter/source/starter/lectur.F
3452!||--- uses -----------------------------------------------------
3453!|| fvbag_mod ../starter/share/modules1/fvbag_mod.f
3454!|| monvol_struct_mod ../starter/share/modules1/monvol_struct_mod.F
3455!||====================================================================
3456 SUBROUTINE fvdim(T_MONVOL)
3457C-----------------------------------------------
3458C M o d u l e s
3459C-----------------------------------------------
3460 USE fvbag_mod
3461 USE monvol_struct_mod
3462C-----------------------------------------------
3463C I m p l i c i t T y p e s
3464C-----------------------------------------------
3465#include "implicit_f.inc"
3466C-----------------------------------------------
3467C C o m m o n B l o c k s
3468C-----------------------------------------------
3469#include "com04_c.inc"
3470C-----------------------------------------------
3471C D u m m y A r g u m e n t s
3472C-----------------------------------------------
3473 TYPE(monvol_struct_), DIMENSION(NVOLU), INTENT(INOUT) :: T_MONVOL
3474C-----------------------------------------------
3475C L o c a l V a r i a b l e s
3476C-----------------------------------------------
3477 INTEGER N, ITYP
3478C
3479 ALLOCATE(fvid(nvolu))
3480 fvid(:) = -1
3481 nfvbag=0
3482 DO n = 1, nvolu
3483 ityp = t_monvol(n)%TYPE
3484 IF (ityp == 6 .OR. ityp == 8) THEN
3485 nfvbag = nfvbag+1
3486 t_monvol(n)%IVOLU(45) = nfvbag
3487 t_monvol(n)%IVOLU(57) = 0
3488 fvid(nfvbag) = n
3489 ENDIF
3490 ENDDO
3491C
3492 RETURN
3493 END
3494!||====================================================================
3495!|| fvdeal ../starter/source/airbag/fvmesh.F
3496!||--- uses -----------------------------------------------------
3497!|| fvbag_mod ../starter/share/modules1/fvbag_mod.F
3498!||====================================================================
3499 SUBROUTINE fvdeal()
3500C-----------------------------------------------
3501C M o d u l e s
3502C-----------------------------------------------
3503 USE fvbag_mod
3504C-----------------------------------------------
3505C I m p l i c i t T y p e s
3506C-----------------------------------------------
3507#include "implicit_f.inc"
3508C-----------------------------------------------
3509C C o m m o n B l o c k s
3510C-----------------------------------------------
3511C-----------------------------------------------
3512C L o c a l V a r i a b l e s
3513C-----------------------------------------------
3514 INTEGER I
3515C
3516 IF (ALLOCATED(FVDATA)) THEN
3517 DO I=1,nfvbag
3518 IF (ASSOCIATED(fvdata(i)%IFVNOD)) DEALLOCATE(fvdata(i)%IFVNOD)
3519 IF (ASSOCIATED(fvdata(i)%RFVNOD)) DEALLOCATE(fvdata(i)%RFVNOD)
3520 IF (ASSOCIATED(fvdata(i)%IFVTRI)) DEALLOCATE(fvdata(i)%IFVTRI)
3521 IF (ASSOCIATED(fvdata(i)%IFVPOLY)) DEALLOCATE(fvdata(i)%IFVPOLY)
3522 IF (ASSOCIATED(fvdata(i)%IFVTADR)) DEALLOCATE(fvdata(i)%IFVTADR)
3523 IF (ASSOCIATED(fvdata(i)%IFVPOLH)) DEALLOCATE(fvdata(i)%IFVPOLH)
3524 IF (ASSOCIATED(fvdata(i)%IFVPADR)) DEALLOCATE(fvdata(i)%IFVPADR)
3525 IF (ASSOCIATED(fvdata(i)%IDPOLH)) DEALLOCATE(fvdata(i)%IDPOLH)
3526 IF (ASSOCIATED(fvdata(i)%IBPOLH)) DEALLOCATE(fvdata(i)%IBPOLH)
3527 IF (ASSOCIATED(fvdata(i)%IFVPOLY_ANIM)) DEALLOCATE(fvdata(i)%IFVPOLY_ANIM)
3528 IF (ASSOCIATED(fvdata(i)%IFVTADR_ANIM)) DEALLOCATE(fvdata(i)%IFVTADR_ANIM)
3529 IF (ASSOCIATED(fvdata(i)%IFVPOLH_ANIM)) DEALLOCATE(fvdata(i)%IFVPOLH_ANIM)
3530 IF (ASSOCIATED(fvdata(i)%IFVPADR_ANIM)) DEALLOCATE(fvdata(i)%IFVPADR_ANIM)
3531 IF (ASSOCIATED(fvdata(i)%REDIR_ANIM)) DEALLOCATE(fvdata(i)%REDIR_ANIM)
3532 IF (ASSOCIATED(fvdata(i)%NOD_ANIM)) DEALLOCATE(fvdata(i)%NOD_ANIM)
3533 IF (ASSOCIATED(fvdata(i)%MPOLH)) DEALLOCATE(fvdata(i)%MPOLH)
3534 IF (ASSOCIATED(fvdata(i)%QPOLH)) DEALLOCATE(fvdata(i)%QPOLH)
3535 IF (ASSOCIATED(fvdata(i)%EPOLH)) DEALLOCATE(fvdata(i)%EPOLH)
3536 IF (ASSOCIATED(fvdata(i)%PPOLH)) DEALLOCATE(fvdata(i)%PPOLH)
3537 IF (ASSOCIATED(fvdata(i)%RPOLH)) DEALLOCATE(fvdata(i)%RPOLH)
3538 IF (ASSOCIATED(fvdata(i)%GPOLH)) DEALLOCATE(fvdata(i)%GPOLH)
3539 IF (ASSOCIATED(fvdata(i)%CPAPOLH)) DEALLOCATE(fvdata(i)%CPAPOLH)
3540 IF (ASSOCIATED(fvdata(i)%CPBPOLH)) DEALLOCATE(fvdata(i)%CPBPOLH)
3541 IF (ASSOCIATED(fvdata(i)%CPCPOLH)) DEALLOCATE(fvdata(i)%CPCPOLH)
3542 IF (ASSOCIATED(fvdata(i)%RMWPOLH)) DEALLOCATE(fvdata(i)%RMWPOLH)
3543 IF (ASSOCIATED(fvdata(i)%VPOLH_INI)) DEALLOCATE(fvdata(i)%VPOLH_INI)
3544 IF (ASSOCIATED(fvdata(i)%DTPOLH)) DEALLOCATE(fvdata(i)%DTPOLH)
3545C
3546 IF (ASSOCIATED(fvdata(i)%BRIC)) DEALLOCATE(fvdata(i)%BRIC)
3547 IF (ASSOCIATED(fvdata(i)%TBRIC)) DEALLOCATE(fvdata(i)%TBRIC)
3548 IF (ASSOCIATED(fvdata(i)%XB)) DEALLOCATE(fvdata(i)%XB)
3549 IF (ASSOCIATED(fvdata(i)%SFAC)) DEALLOCATE(fvdata(i)%SFAC)
3550C
3551 IF (ASSOCIATED(fvdata(i)%TPOLH)) DEALLOCATE(fvdata(i)%TPOLH)
3552 IF (ASSOCIATED(fvdata(i)%CPDPOLH)) DEALLOCATE(fvdata(i)%CPDPOLH)
3553 IF (ASSOCIATED(fvdata(i)%CPEPOLH)) DEALLOCATE(fvdata(i)%CPEPOLH)
3554 IF (ASSOCIATED(fvdata(i)%CPFPOLH)) DEALLOCATE(fvdata(i)%CPFPOLH)
3555 ENDDO
3556 DEALLOCATE(fvdata)
3557 ENDIF
3558C
3559 IF (ALLOCATED(fvspmd)) THEN
3560 DO i=1,nfvbag
3561 IF (ASSOCIATED(fvspmd(i)%IBUF_L)) DEALLOCATE(fvspmd(i)%IBUF_L)
3562 IF (ASSOCIATED(fvspmd(i)%IBUFA_L)) DEALLOCATE(fvspmd(i)%IBUFA_L)
3563 IF (ASSOCIATED(fvspmd(i)%IBUFSA_L)) DEALLOCATE(fvspmd(i)%IBUFSA_L)
3564 IF (ASSOCIATED(fvspmd(i)%IXSA)) DEALLOCATE(fvspmd(i)%IXSA)
3565 IF (ASSOCIATED(fvspmd(i)%ELEMSA)) DEALLOCATE(fvspmd(i)%ELEMSA)
3566 ENDDO
3567 DEALLOCATE(fvspmd)
3568 ENDIF
3569C
3570 RETURN
3571 END
my_exit
void my_exit(int *i)
Definition analyse.c:1038
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
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
fvnormal
subroutine fvnormal(x, n1, n2, n3, n4, nx, ny, nz)
Definition fvmbag1.F:576
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
horiedge
subroutine horiedge(ipoly, rpoly, nx, ny, nz, nbnedge, inedge, rnedge, x0, y0, z0, inz, nns3, nref, aref, nnsp)
Definition horipoly.F:32
horipoly
subroutine horipoly(inedge, rnedge, ledge, nedge, ipoly, rpoly, iz, ielnod, npoly, nx, ny, nz, inz, ibric, nel, tagela)
Definition horipoly.F:113
itribox
subroutine itribox(tri, box, norm, nverts, poly, nvmax)
Definition itribox.F:31
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
fvbag_gpolcut_mod
Definition gpolcut.F:28
fvbag_gpolcut_mod::fvbag_gpolcut
subroutine fvbag_gpolcut(ipoly, rpoly, ipoly_old, rpoly_old, inedge, rnedge, nbnedge, nx, ny, nz, x0, y0, z0, ins, rns, nn, nhol, inz, iz, nns3, npoly, ns, ielnod, ielnod_old)
Definition gpolcut.F:40
fvbag_mod
Definition fvbag_mod.F:44
fvbag_mod::fvspmd
type(fvbag_spmd), dimension(:), allocatable fvspmd
Definition fvbag_mod.F:129
fvbag_mod::fvdata
type(fvbag_data), dimension(:), allocatable fvdata
Definition fvbag_mod.F:128
fvbag_mod::fvid
integer, dimension(:), allocatable fvid
Definition fvbag_mod.F:130
fvbag_mod::nfvbag
integer nfvbag
Definition fvbag_mod.F:127
message_mod
Definition message_mod.F:1249
monvol_struct_mod
Definition monvol_struct_mod.F:66
names_and_titles_mod
Definition names_and_titles_mod.F:997
names_and_titles_mod::nchartitle
integer, parameter nchartitle
Definition names_and_titles_mod.F:1003
reorder_mod
Definition reorder_mod.F:45
polyhedr1
subroutine polyhedr1(ipoly, rpoly, polb, npolb, polh, npolh, nrpmax, nphmax, ibric, lmin, info, npolhmax, nppmax, nel, inz, tagela)
Definition polyhedr1.F:35
polyhedr
subroutine polyhedr(ipoly, rpoly, polb, npolb, polh, npolh, nrpmax, nphmax, ibric, lmin, info, npolhmax, nppmax)
Definition polyhedr.F:31
quicksort
recursive subroutine quicksort(a, idx, first, last)
Definition quicksort.F:34
fvdeal
subroutine fvdeal()
Definition fvmesh.F:3500
fvmesh1
subroutine fvmesh1(ibuf, elem, x, ivolu, bric, xb, rvolu, nel, neli, nbric, tbric, sfac, dxm, nba, nela, tba, tfaca, tagels, ibufa, elema, tagela, ixs, id, titr, nb_node, ityp)
Definition fvmesh.F:53
fvdim
subroutine fvdim(t_monvol)
Definition fvmesh.F:3457
coorloc
subroutine coorloc(vpx, vpy, vpz, v1x, v1y, v1z, v2x, v2y, v2z, ksi, eta)
Definition fvmesh.F:3418
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
c_tricall
subroutine c_tricall(pnodes, pseg, pholes, ptri, nnp, nseg, nhol, nelp)
Definition fvmbags_stub.F:34
tritri3
subroutine tritri3(icut, fv0, fv1, fv2, fu0, fu1, fu2)
Definition fvmbags_stub.F:74
tribox3
subroutine tribox3(icut, bcenter, bhalfsize, tverts)
Definition fvmbags_stub.F:55
starter
program starter
Definition starter.F:39
monvol_struct_mod::monvol_struct_
Definition monvol_struct_mod.F:84
tmass
subroutine tmass(x, nc, geo, pm, ms, stifn, partsav, v, ipart, mst, stifint, stt, area, mxt, nc1, nc2, x1, x2, y1, y2, z1, z2)
Definition tmass.F:33