1 subroutine pppma3 ( nbtrvi, nbquvi,
4 > posini, xyzfac, tabaux, nivsup,
5 > ulsort, langue, codret )
6 c ______________________________________________________________________
10 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
12 c Version originale enregistree le 18 juin 1996 sous le numero 96036
13 c aupres des huissiers de justice Simart et Lavoir a Clamart
14 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
15 c aupres des huissiers de justice
16 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
18 c HOMARD est une marque deposee d'Electricite de France
24 c ______________________________________________________________________
26 c Post-Processeur - Preparation du MAillage - phase 3
28 c ______________________________________________________________________
30 c but : recherche de l'ordre d'affichage des faces pour les cacher
32 c par defaut, on affiche les objets vus par l'observateur avec
33 c l'axe (oz+) dans l'oeil, donc regardant de z>0 vers z<0.
34 c on utilise l'algorithme dit du peintre, ou encore du z-buffer
35 c utilise par les affichages graphiques standards, consistant a
36 c imprimer les objets dans l'ordre inverse de leur eloignement,
37 c i.e. a imprimer en dernier, et donc par dessus le reste, les
38 c objets les plus proches. (donc avec les z les plus grands)
39 c ______________________________________________________________________
41 c . nom . e/s . taille . description .
42 c .____________________________________________________________________.
43 c . nbtrvi . e . 1 . nombre triangles visualisables .
44 c . nbquvi . e . 1 . nombre de quadrangles visualisables .
45 c . nntrvi . e .10nbtrvi. 1 : niveau du triangle a afficher .
46 c . . . . 2 : numero HOMARD du triangle .
47 c . . . . 3, 4, 5 : numeros des noeuds p1 .
48 c . . . . 6 : famille du triangle .
49 c . . . . 7, 8, 9 : numeros des noeuds p2 .
50 c . . . . 10 : numero du noeud interne .
51 c . nnquvi . e .12nbquvi. 1 : niveau du quadrangle a afficher .
52 c . . . . 2 : numero HOMARD du quadrangle .
53 c . . . . 3, 4, 5, 6 : numeros des noeuds p1 .
54 c . . . . 7 : famille du quadrangle .
55 c . . . . 8, 9, 10, 11 : numeros des noeuds p2 .
56 c . . . . 12 : numero du noeud interne .
57 c . coopro . e . 3* . coordonnees projetees de : .
58 c . . .nbnot+12. le triedre : -8:O ; -9:I ; -10:J ; -11:K .
59 c . . . . la fenetre de zoom : de -7 a 0 en 3D ou .
60 c . . . . de -3 a 0 en 2D .
61 c . . . . les noeuds de 1 a nbnoto .
62 c . posini . aux . nbquvi . tableau auxiliaire de renumerotation des .
63 c . . .+nbtrvi . faces en fonction de l'affichage .
64 c . xyzfac . / .nbtrvi+ . tableau de travail reel .
66 c . tabaux . aux . nbquvi . tableau auxiliaire .
68 c . nivsup . e . 1 . niveau superieur present dans le maillage .
69 c . ulsort . e . 1 . unite logique de la sortie generale .
70 c . langue . e . 1 . langue des messages .
71 c . . . . 1 : francais, 2 : anglais .
72 c . codret . s . 1 . code de retour des modules .
73 c . . . . 0 : pas de probleme .
74 c ._____________________________________________________________________
77 c 0. declarations et dimensionnement
80 c 0.1. ==> generalites
86 parameter ( nompro = 'PPPMA3' )
99 integer nbtrvi, nbquvi
100 integer nntrvi(10,nbtrvi)
101 integer nnquvi(12,nbquvi)
102 integer posini(nbtrvi+nbquvi)
103 integer tabaux(nbtrvi+nbquvi)
105 double precision coopro(3,-11:nbnoto)
106 double precision xyzfac(nbtrvi+nbquvi,9)
108 integer ulsort, langue, codret
110 c 0.4. ==> variables locales
112 integer nutrvi, nuquvi
114 integer rangfa, rangqu
115 integer nbfast, nbfas0
119 double precision daux
120 double precision borne
121 double precision xminfa, xmaxfa
122 double precision yminfa, ymaxfa
123 double precision zminfa, zmaxfa
124 double precision coface(2,3)
125 double precision v1(3), v2(3), v3(3)
128 parameter ( nbmess = 10 )
129 character*80 texte(nblang,nbmess)
131 c 0.5. ==> initialisations
135 c 0.5. ==> initialisations
136 c_______________________________________________________________________
144 #ifdef _DEBUG_HOMARD_
145 write (ulsort,texte(langue,1)) 'Entree', nompro
151 #ifdef _DEBUG_HOMARD_
152 write (ulsort,90002) 'nbtrvi', nbtrvi
153 write (ulsort,90002) 'nbquvi', nbquvi
154 write (ulsort,90002) 'sdim ', sdim
155 write (ulsort,90002) 'nivsup', nivsup
161 c 2. En dimension 2, les faces sont reclassses par niveau croissant
162 c A la fin de cette etape, posini(1) contient l'indice de la face a
163 c tracer en premier, parce que la plus loin, posini(2) contient
164 c l'indice de la suivante a tracer et ainsi de suite, jusqu'a
165 c posini(nbfavi) qui contient l'indice de la derniere face
166 c a tracer parce que la plus proche de l'observateur.
167 c Attention : posini ne contient pas les numeros des faces mais une
168 c indirection dans la liste des faces a traiter.
170 cgn do 34 , jaux = 1 , nbnoto
171 cgn write(*,1797)jaux,(coopro(iaux,jaux),iaux=1,sdim)
174 if ( sdim.le.2 ) then
178 do 21 , iaux = 0 , nivsup+1
180 do 211 , nutrvi = 1 , nbtrvi
181 if ( nntrvi(1,nutrvi).eq.iaux ) then
183 posini(jaux) = nutrvi
187 do 212 , nuquvi = 1 , nbquvi
188 if ( nnquvi(1,nuquvi).eq.iaux ) then
190 posini(jaux) = nbtrvi+nuquvi
197 c 2. En dimension 3, on applique l'algorithme du z-buffer
203 c 3. Calcul des dimensions extremes des faces et des vecteurs normaux
205 c 3.1. ==> Extrema pour un triangle
207 do 31 , nutrvi = 1 , nbtrvi
209 xminfa = min(coopro(1,nntrvi(3,nutrvi)),
210 > coopro(1,nntrvi(4,nutrvi)),coopro(1,nntrvi(5,nutrvi)))
211 xmaxfa = max(coopro(1,nntrvi(3,nutrvi)),
212 > coopro(1,nntrvi(4,nutrvi)),coopro(1,nntrvi(5,nutrvi)))
213 xyzfac(nutrvi,1) = xminfa + borne * (xmaxfa - xminfa)
214 xyzfac(nutrvi,2) = xmaxfa - borne * (xmaxfa - xminfa)
216 yminfa = min(coopro(2,nntrvi(3,nutrvi)),
217 > coopro(2,nntrvi(4,nutrvi)),coopro(2,nntrvi(5,nutrvi)))
218 ymaxfa = max(coopro(2,nntrvi(3,nutrvi)),
219 > coopro(2,nntrvi(4,nutrvi)),coopro(2,nntrvi(5,nutrvi)))
220 xyzfac(nutrvi,3) = yminfa + borne * (ymaxfa - yminfa)
221 xyzfac(nutrvi,4) = ymaxfa - borne * (ymaxfa - yminfa)
223 zminfa = min(coopro(3,nntrvi(3,nutrvi)),
224 > coopro(3,nntrvi(4,nutrvi)),coopro(3,nntrvi(5,nutrvi)))
225 zmaxfa = max(coopro(3,nntrvi(3,nutrvi)),
226 > coopro(3,nntrvi(4,nutrvi)),coopro(3,nntrvi(5,nutrvi)))
227 xyzfac(nutrvi,5) = zminfa + borne * (zmaxfa - zminfa)
228 xyzfac(nutrvi,6) = zmaxfa - borne * (zmaxfa - zminfa)
229 cgn write(*,1798)nutrvi,(nntrvi(iaux,nutrvi),iaux=2,5),
230 cgn > (xyzfac(nutrvi,iaux),iaux=1,6)
234 c 3.2. ==> Extrema pour un quadrangle
236 do 32 , nuquvi = 1 , nbquvi
238 xminfa = coopro(1,nnquvi(3,nuquvi))
240 yminfa = coopro(2,nnquvi(3,nuquvi))
242 zminfa = coopro(3,nnquvi(3,nuquvi))
244 do 321 , iaux = 4 , 6
245 daux = coopro(1,nnquvi(iaux,nuquvi))
246 xminfa = min ( daux , xminfa )
247 xmaxfa = max ( daux , xmaxfa )
248 daux = coopro(2,nnquvi(iaux,nuquvi))
249 yminfa = min ( daux , yminfa )
250 ymaxfa = max ( daux , ymaxfa )
251 daux = coopro(3,nnquvi(iaux,nuquvi))
252 zminfa = min ( daux , zminfa )
253 zmaxfa = max ( daux , zmaxfa )
255 rangfa = nbtrvi+nuquvi
256 xyzfac(rangfa,1) = xminfa + borne * (xmaxfa - xminfa)
257 xyzfac(rangfa,2) = xmaxfa - borne * (xmaxfa - xminfa)
258 xyzfac(rangfa,3) = yminfa + borne * (ymaxfa - yminfa)
259 xyzfac(rangfa,4) = ymaxfa - borne * (ymaxfa - yminfa)
260 xyzfac(rangfa,5) = zminfa + borne * (zmaxfa - zminfa)
261 xyzfac(rangfa,6) = zmaxfa - borne * (zmaxfa - zminfa)
262 cgn write(*,1799)nuquvi,(nnquvi(iaux,nuquvi),iaux=2,6),
263 cgn > (xyzfac(nuquvi,iaux),iaux=1,6)
267 c 3.3. ==> Vecteurs normaux
269 nbfavi = nbtrvi + nbquvi
270 cgn write (ulsort,1795) 'nbfavi', nbfavi
271 do 33 , rangfa = 1 , nbfavi
273 if ( rangfa.le.nbtrvi ) then
274 v1(1) = coopro(1,nntrvi(3,rangfa))
275 v1(2) = coopro(2,nntrvi(3,rangfa))
276 v1(3) = coopro(3,nntrvi(3,rangfa))
277 v2(1) = coopro(1,nntrvi(4,rangfa))
278 v2(2) = coopro(2,nntrvi(4,rangfa))
279 v2(3) = coopro(3,nntrvi(4,rangfa))
280 v3(1) = coopro(1,nntrvi(5,rangfa))
281 v3(2) = coopro(2,nntrvi(5,rangfa))
282 v3(3) = coopro(3,nntrvi(5,rangfa))
284 rangqu = rangfa - nbtrvi
285 v1(1) = coopro(1,nnquvi(3,rangqu))
286 v1(2) = coopro(2,nnquvi(3,rangqu))
287 v1(3) = coopro(3,nnquvi(3,rangqu))
288 v2(1) = coopro(1,nnquvi(4,rangqu))
289 v2(2) = coopro(2,nnquvi(4,rangqu))
290 v2(3) = coopro(3,nnquvi(4,rangqu))
291 v3(1) = coopro(1,nnquvi(5,rangqu))
292 v3(2) = coopro(2,nnquvi(5,rangqu))
293 v3(3) = coopro(3,nnquvi(5,rangqu))
295 cgn write(*,1796)'v1',v1
296 cgn write(*,1796)'v2',v2
297 cgn write(*,1796)'v3',v3
299 coface(1,1) = v2(1) - v1(1)
300 coface(1,2) = v2(2) - v1(2)
301 coface(1,3) = v2(3) - v1(3)
302 coface(2,1) = v3(1) - v1(1)
303 coface(2,2) = v3(2) - v1(2)
304 coface(2,3) = v3(3) - v1(3)
305 cgn write(*,1796)'v1',v1
306 cgn write(*,1796)'v2',v2
307 cgn write(*,1796)'v3',v3
309 xyzfac(rangfa,7) = coface(1,2)*coface(2,3)
310 > - coface(1,3)*coface(2,2)
311 xyzfac(rangfa,8) = coface(1,3)*coface(2,1)
312 > - coface(1,1)*coface(2,3)
313 xyzfac(rangfa,9) = coface(1,1)*coface(2,2)
314 > - coface(1,2)*coface(2,1)
315 cgn write(*,1796)'normale',(xyzfac(rangfa,iaux),iaux=7,9)
319 1795 format(3(a,' =',i5,' , '))
320 1796 format(a,6f12.5)
321 1797 format(i5,' *',6f12.5)
322 1798 format(i4,' :',i5,' *',3i4,' *',6f12.5)
323 1799 format(i4,' :',i5,' *',4i4,' *',6f12.5)
326 c 4. On parcourt toutes les faces a afficher et on les range de la
327 c plus eloignee a la plus proche du point de vue.
328 c Groso modo, les plus proches sont avec les z les plus grands.
329 c A la fin de cette etape, posini(1) contient l'indice de la face a
330 c tracer en premier, parce que la plus loin, posini(2) contient
331 c l'indice de la suivante a tracer et ainsi de suite, jusqu'a
332 c posini(nbfavi) qui contient l'indice de la derniere face
333 c a tracer parce que la plus proche de l'observateur.
334 c Attention : posini ne contient pas les numeros des faces mais une
335 c indirection dans la liste des faces a traiter.
338 c 4.1. ==> A priori, les faces ne sont en contact avec aucune autre
340 do 41 , rangfa = 1 , nbfavi
344 c 4.2. ==> A priori, la 1ere face est en premiere position
350 c 4.3. ==> On examine toutes les faces en ne classant que celles qui
351 c ont une partie commune avec celles deja classees
352 c On n'examine que les faces non encore stockees
353 c On boucle jusqu'a ce qu'il n'y en ait plus
355 c . tabaux(i) = 1 : la face i est stockee
356 c . tabaux(i) = 0 : la face i n'est pas stockee
362 cgn write (ulsort,*) '---------- NOUVELLE SERIE ---------'
363 cgn write (ulsort,1795) 'Depart avec rangfd', rangfd,'nbfast',nbfast
367 do 43 , rangfa = rangfd , nbfavi
369 if ( codret.eq.0 ) then
371 if ( tabaux(rangfa).eq.0 ) then
374 #ifdef _DEBUG_HOMARD_
375 write (ulsort,texte(langue,3)) 'PPPMA4', nompro
377 call pppma4 ( iaux, nbfast,
381 > posini, xyzfac, tabaux,
382 > ulsort, langue, codret )
384 if ( rangfa.eq.-70) stop
393 if ( codret.eq.0 ) then
395 cgn write (ulsort,*) 'nbfas0, nbfast =', nbfas0, nbfast
396 cgn if ( nbfast.eq.222 ) then
401 c . Si le nombre de faces stockees est egal au nombre de faces a
402 c visualiser, c'est fini
404 cgn write (ulsort,*) 'Bilan'
405 if ( nbfast.eq.nbfavi ) then
406 cgn write (ulsort,1795) 'Fin avec nbfast', nbfast
409 c . Si le nombre de faces stockees a change, c'est que l'on
410 c se trouve encore dans une suite de faces qui se superposent
411 c . Sinon, c'est qu'il ne reste plus de faces se superposant
412 c a l'une de la serie en cours. Il faut commencer une
414 c Dans les deux cas, on va refaire un tour en partant de la
415 c premiere face non superposee.
418 cgn do 998 , iaux = 1 , nbfavi
419 cgn if ( iaux.le.nbtrvi ) then
420 cgn jaux = nntrvi(2,iaux)
422 cgn jaux = -nnquvi(2,iaux-nbtrvi)
424 cgn if ( tabaux(iaux).ne.0 ) then
425 cgn write (*,1999) iaux,tabaux(iaux),jaux
426 cgn 1999 format('tabaux(',i5,') =',i5,', face',i5)
430 do 44 , rangfa = rangfd , nbfavi
431 if ( tabaux(rangfa).eq.0 ) then
439 cgn write (ulsort,1795) 'nbfas0', nbfas0, 'nbfast', nbfast,
440 cgn > 'rangfd', rangfd
441 if ( nbfast.eq.nbfas0 ) then
442 cgn write (ulsort,*) '---------- NOUVELLE SERIE ---------'
443 cgn write (ulsort,*) '---------- FACE', -nnquvi(2,rangfd-nbtrvi)
445 cgn write (ulsort,888)
446 cgn > (rangfa,posini(rangfa),rangfa =1,nbtrvi+nbquvi)
447 cgn 888 format(5(i3,i4,' * '))
448 posini(nbfast) = rangfd
453 cgn write (ulsort,*) 'goto 430 goto 430 goto 430'
460 c 4.5. ==> C'est fini
464 cgn nbfavi=nbtrvi+nbquvi
465 cgn do 999 , iaux = 1 , nbfast
466 cgn if ( posini(iaux).le.nbtrvi ) then
467 cgn jaux = nntrvi(2,posini(iaux))
469 cgn jaux = -nnquvi(2,posini(iaux)-nbtrvi)
471 cgn write (*,*) iaux,posini(iaux),'-eme face :',jaux
480 if ( codret.ne.0 ) then
484 write (ulsort,texte(langue,1)) 'Sortie', nompro
485 write (ulsort,texte(langue,2)) codret
489 #ifdef _DEBUG_HOMARD_
490 write (ulsort,texte(langue,1)) 'Sortie', nompro