1 subroutine decr05 ( tyconf, homolo,
3 > hetare, filare, posifa, facare,
4 > hettri, aretri, voltri,
5 > hetqua, arequa, volqua,
6 > tritet, quahex, coquhe,
9 > ulsort, langue, codret )
10 c ______________________________________________________________________
14 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
16 c Version originale enregistree le 18 juin 1996 sous le numero 96036
17 c aupres des huissiers de justice Simart et Lavoir a Clamart
18 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
19 c aupres des huissiers de justice
20 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
22 c HOMARD est une marque deposee d'Electricite de France
28 c ______________________________________________________________________
30 c traitement des DEcisions - Contraintes de Raffinement - 05
32 c Pas de segments decoupes sans sa face voisine, ni de face decoupee
33 c sans son volume voisin
34 c Il faut faire ce controle a la fin de l'algorithme sur la
35 c propagation du raffinement, car on ne peut pas prevoir au depart
36 c tout ce qui va se passer. En particulier dans des cas bizarres pour
37 c lesquels on aurait plusieurs boites.
38 c ______________________________________________________________________
40 c . nom . e/s . taille . description .
41 c .____________________________________________________________________.
42 c . tyconf . e . 1 . 0 : conforme (defaut) .
43 c . . . . 1 : non-conforme avec au minimum 2 aretes .
44 c . . . . non decoupees en 2 .
45 c . . . . 2 : non-conforme avec 1 seul noeud .
46 c . . . . pendant par arete .
47 c . . . . 3 : non-conforme fidele a l'indicateur .
48 c . . . . -1 : conforme, avec des boites pour les .
49 c . . . . quadrangles, hexaedres et pentaedres .
50 c . . . . -2 : non-conforme avec au maximum 1 arete .
51 c . . . . decoupee en 2 (boite pour les .
52 c . . . . quadrangles, hexaedres et pentaedres) .
53 c . . . . -2 : non-conforme avec au maximum 1 arete .
54 c . . . . decoupee en 2 (boite pour les .
55 c . . . . quadrangles, hexaedres et pentaedres) .
56 c . homolo . e . 1 . type de relations par homologues .
57 c . . . . 0 : pas d'homologues .
58 c . . . . 1 : relations sur les noeuds .
59 c . . . . 2 : relations sur les noeuds et les aretes .
60 c . . . . 3 : relations sur les noeuds, les aretes .
61 c . . . . et les triangles .
62 c . decfac . e . -nbquto. decision sur les faces (quad. + tri.) .
64 c . decare . es . nbarto . decisions des aretes .
65 c . hetare . e . nbarto . historique de l'etat des aretes .
66 c . filare . e . nbarto . premiere fille des aretes .
67 c . posifa . e . nbarto . pointeur sur tableau facare .
68 c . facare . e . nbfaar . liste des faces contenant une arete .
69 c . hettri . e . nbtrto . historique de l'etat des triangles .
70 c . aretri . e .nbtrto*3. numeros des 3 aretes des triangles .
71 c . voltri . e .2*nbtrto. numeros des 2 volumes par triangle .
72 c . . . . voltri(i,k) definit le i-eme voisin de k .
73 c . . . . 0 : pas de voisin .
74 c . . . . j>0 : tetraedre j .
75 c . . . . j<0 : pyramide/pentaedre dans pypetr(1/2,j).
76 c . hetqua . e . nbquto . historique de l'etat des quadrangles .
77 c . arequa . e .nbquto*4. numeros des 4 aretes des quadrangles .
78 c . volqua . e .2*nbquto. numeros des 2 volumes par quadrangle .
79 c . . . . volqua(i,k) definit le i-eme voisin de k .
80 c . . . . 0 : pas de voisin .
81 c . . . . j>0 : hexaedre j .
82 c . . . . j<0 : pyramide/pentaedre dans pypequ(1/2,j).
83 c . tritet . e .nbtecf*4. numeros des 4 triangles des tetraedres .
84 c . quahex . e .nbhecf*6. numeros des 6 quadrangles des hexaedres .
85 c . coquhe . e .nbhecf*6. codes des 6 quadrangles des hexaedres .
86 c . arehom . e . nbarto . ensemble des aretes homologues .
87 c . afaire . es . 1 . que faire a la sortie .
88 c . . . . 0 : aucune action .
89 c . . . . 1 : refaire une iteration de l'algorithme .
90 c . ulsort . e . 1 . numero d'unite logique de la liste standard.
91 c . langue . e . 1 . langue des messages .
92 c . . . . 1 : francais, 2 : anglais .
93 c . codret . es . 1 . code de retour des modules .
94 c . . . . 0 : pas de probleme .
95 c . . . . sinon : probleme .
96 c ______________________________________________________________________
99 c 0. declarations et dimensionnement
102 c 0.1. ==> generalites
108 parameter ( nompro = 'DECR05' )
125 integer tyconf, homolo
126 integer decfac(-nbquto:nbtrto), decare(0:nbarto)
127 integer hetare(nbarto), filare(nbarto)
128 integer posifa(0:nbarto), facare(nbfaar)
129 integer hettri(nbtrto), aretri(nbtrto,3), voltri(2,nbtrto)
130 integer hetqua(nbquto), arequa(nbquto,4), volqua(2,nbquto)
131 integer tritet(nbtecf,4)
132 integer quahex(nbhecf,6), coquhe(nbhecf,6)
133 integer arehom(nbarto)
137 integer ulsort, langue, codret
139 c 0.4. ==> variables locales
142 integer iaux, jaux, kaux
144 integer etatar, etatfa
145 integer larete, laret1, larelo, laface, iface, letetr, lehexa
146 integer nbarpb, nbfapb
147 integer nbaret, listar(12)
148 #ifdef _DEBUG_HOMARD_
153 parameter ( nbmess = 30 )
154 character*80 texte(nblang,nbmess)
156 c 0.5. ==> initialisations
157 c ______________________________________________________________________
165 #ifdef _DEBUG_HOMARD_
166 write (ulsort,texte(langue,1)) 'Entree', nompro
171 >'(5x,''Pas de maille de bord decoupe sans son voisin.'',/)'
172 texte(1,5) = '(7x,''Nombre de '',a,''a reconsiderer :'',i6,/)'
173 texte(1,6) = '(7x,''Aucun changement.'')'
174 texte(1,7) = '(7x,''Apres l''''analyse '',a)'
175 texte(1,8) = '(a,''numero '',i8,'' : decision ='',i2)'
177 > '(a,''numero '',i8,'' : decision ='',i2,'', etat ='',i5)'
178 texte(1,10) = '(/,i1,''. Examen des'',i10,1x,a,)'
181 > '(5x,''No border mesh cut without its neighbour.'',/)'
182 texte(2,5) = '(7x,''Number of'',a,''to deal with :'',i6,/)'
183 texte(2,6) = '(7x,''No modification.'')'
184 texte(2,7) = '(7x,''After analysis '',a)'
185 texte(2,8) = '(a,''#'',i8,'' : decision ='',i25)'
187 > '(a,''#'',i8,'' : decision ='',i2,'', status ='',i5)'
188 texte(2,10) = '(/,''Examination of the'',i10,1x,a,)'
196 write (ulsort,texte(langue,4))
202 c 2. on interdit les situations ou on aurait un segment decoupe alors
203 c qu'aucune des faces auxquelles il appartient ne le serait.
204 c Cela peut arriver si on a fait du decoupage selon une zone
205 c geometrique et que cette zone incluait une serie d'aretes.
206 c Ou avec un indicateur sur aretes ou noeuds et que ce seul
207 c segment a ete retenu.
209 #ifdef _DEBUG_HOMARD_
210 write (ulsort,90002) '2. face/arete ; codret', codret
213 #ifdef _DEBUG_HOMARD_
214 write (ulsort,texte(langue,10)) 1, nbarto, mess14(langue,3,1)
217 do 20 , larete = 1 , nbarto
219 #ifdef _DEBUG_HOMARD_
220 write (ulsort,texte(langue,8)) mess14(langue,1,1),
221 > larete,decare(larete)
224 if ( decare(larete).eq.2 ) then
226 c 2.1. ==> on parcourt chacune des faces voisines de l'arete
227 c on compte le nombre de faces a couper ou a reactualiser
228 c s'il y a des equivalences, il faut traiter ensemble une
229 c arete et son homologue
235 if ( homolo.ge.2 ) then
236 laret1 = arehom(larete)
237 if ( laret1.ne.0 ) then
238 listar(2) = abs(laret1)
243 do 211 , iaux = 1 , nbaret
245 laret1 = listar(iaux)
247 ideb = posifa(laret1-1)+1
248 ifin = posifa(laret1)
250 do 2111 , ipos = ideb , ifin
253 if ( iface.gt.0 ) then
254 etatfa = mod( hettri(iface) , 10 )
256 etatfa = mod( hetqua(-iface) , 100 )
258 #ifdef _DEBUG_HOMARD_
259 write (ulsort,texte(langue,9))'.. '//mess14(langue,1,8),
260 > abs(iface),decfac(iface),etatfa
262 if ( etatfa.ne.0 .and. decfac(iface).ne.-1 ) then
264 else if ( decfac(iface).eq.4 ) then
274 c 2.2. ==> aucune face n'est a couper ou a reactualiser, on ne doit pas
277 if ( kaux.gt.0 ) then
279 do 22 , iaux = 1 , nbaret
280 laret1 = listar(iaux)
282 #ifdef _DEBUG_HOMARD_
283 write (ulsort,texte(langue,30))'decare', laret1,decare(laret1),' '
293 #ifdef _DEBUG_HOMARD_
294 write (ulsort,texte(langue,7)) 'face/arete'
295 if ( nbarpb.gt.0 .or. nbfapb.gt.0 ) then
296 write (ulsort,texte(langue,5)) mess14(langue,3,1), nbarpb
297 write (ulsort,texte(langue,5)) mess14(langue,3,8), nbfapb
299 write (ulsort,texte(langue,6))
304 c 3. on interdit les situations ou on aurait un triangle decoupe alors
305 c qu'aucun de ses tetraedres voisins ne le serait.
306 c Cela peut arriver si on a fait du decoupage selon une zone
307 c geometrique et que cette zone incluait une zone purement 2D.
308 c Ou avec un indicateur sur faces ou noeuds et que cette seule
309 c face a ete retenue.
311 #ifdef _DEBUG_HOMARD_
312 write (ulsort,90002) '3. tetr/tria ; codret', codret
315 if ( nbteto.ne.0 ) then
317 #ifdef _DEBUG_HOMARD_
318 write (ulsort,texte(langue,10)) 2, nbtrto, mess14(langue,3,2)
321 do 30 , laface = 1 , nbtrto
323 #ifdef _DEBUG_HOMARD_
324 write (ulsort,texte(langue,8)) mess14(langue,1,2),
325 > laface,decfac(laface)
328 if ( decfac(laface).eq.4 ) then
332 c 3.1. ==> on parcourt chacun des tetraedres voisins du triangle
333 c un tetraedre sera coupe si au moins une autre de ses faces
335 c ATTENTION A FAIRE COMME LES HEXAS
339 letetr = voltri(iaux,laface)
341 if ( letetr.gt.0 ) then
343 do 311 , jaux = 1 , 4
345 iface = tritet(letetr,jaux)
347 if ( iface.ne.laface ) then
349 etatfa = mod( hettri(iface) , 10 )
350 #ifdef _DEBUG_HOMARD_
351 write (ulsort,texte(langue,9))'.. '//mess14(langue,1,2),
352 > iface,decfac(iface),etatfa
354 if ( etatfa.ne.0 .and. decfac(iface).ne.-1 ) then
356 else if ( decfac(iface).eq.4 ) then
370 c 3.2. ==> aucun tetraedre n'est a couper ou a reactualiser, on ne doit
371 c pas couper le triangle, ni ses aretes
373 if ( kaux.gt.0 ) then
376 #ifdef _DEBUG_HOMARD_
377 write (ulsort,texte(langue,30))'decfac', laface,decfac(laface),' '
379 do 32 , larelo = 1 , 3
380 larete = aretri(laface,larelo)
381 if ( decare(larete).eq.2 ) then
384 #ifdef _DEBUG_HOMARD_
385 write (ulsort,texte(langue,30))'decare', larete,decare(larete),' '
390 #ifdef _DEBUG_HOMARD_
400 #ifdef _DEBUG_HOMARD_
401 write (ulsort,texte(langue,7)) 'tetr/tria'
402 if ( nbarpb.gt.0 .or. nbfapb.gt.0 ) then
403 write (ulsort,texte(langue,5)) mess14(langue,3,1), nbarpb
404 write (ulsort,texte(langue,5)) mess14(langue,3,8), nbfapb
406 write (ulsort,texte(langue,6))
411 c 4. on interdit les situations ou on aurait un quadrangle decoupe alors
412 c qu'aucun de ses hexaedres voisins ne le serait.
413 c Cela peut arriver si on a fait du decoupage selon une zone
414 c geometrique et que cette zone incluait une zone purement 2D.
415 c Ou avec un indicateur sur faces ou noeuds et que cette seule
416 c face a ete retenue.
417 c Cela peut aussi arriver par contamination entre faces.
419 #ifdef _DEBUG_HOMARD_
420 write (ulsort,90002) '4. hexa/quad ; codret', codret
423 if ( nbheto.ne.0 ) then
425 #ifdef _DEBUG_HOMARD_
426 write (ulsort,texte(langue,10)) 3, nbquto, mess14(langue,3,4)
429 do 40 , laface = 1 , nbquto
431 #ifdef _DEBUG_HOMARD_
432 if ( laface.eq.215996 .or.
433 > laface.eq.66980 ) then
438 if ( glop.gt.0 ) then
439 write (ulsort,texte(langue,8)) mess14(langue,1,4),
440 > laface,decfac(-laface)
441 write (ulsort,*) ' volqua(*,laface) : ',
442 > volqua(1,laface),volqua(2,laface)
443 do 401 , iaux = 1 , 2
444 lehexa = volqua(iaux,laface)
445 if ( lehexa.gt.0 ) then
446 write (ulsort,*)'.. hexaedre ', lehexa
447 do 4011 , jaux = 1 , 6
448 iface = quahex(lehexa,jaux)
449 etatfa = mod( hetqua(iface) , 100 )
450 write (ulsort,texte(langue,9))'.. '//mess14(langue,1,4),
451 > iface,decfac(-iface),etatfa
458 if ( decfac(-laface).eq.4 ) then
460 c 4.1. ==> on parcourt chacun des hexaedres voisins du quadrangle
461 c un hexaedre sera coupe si toutes ses faces le sont
462 c kaux = nombre de faces coupees ou a couper pour le
467 lehexa = volqua(iaux,laface)
471 if ( lehexa.gt.0 ) then
472 #ifdef _DEBUG_HOMARD_
473 if ( glop.gt.0 ) then
474 write (ulsort,*)'.. hexaedre ', lehexa
478 do 411 , jaux = 1 , 6
480 iface = quahex(lehexa,jaux)
482 if ( iface.ne.laface ) then
484 etatfa = mod( hetqua(iface) , 100 )
485 #ifdef _DEBUG_HOMARD_
486 if ( glop.gt.0 ) then
487 write (ulsort,texte(langue,9))'.. '//mess14(langue,1,4),
488 > iface,decfac(-iface),etatfa
491 if ( etatfa.ne.0 .and. decfac(-iface).ne.-1 ) then
493 else if ( decfac(-iface).eq.4 ) then
501 c les 6 faces de l'hexaedre seront coupees, donc RAS
503 if ( kaux.eq.6 ) then
511 c 4.2. ==> si on arrive ici, c'est qu'aucun des hexaedres voisins
513 c 2 cas se presentent :
514 c A. . si on est en mode non-conforme fidele a l'indicateur
515 c . ou si les aretes de chacun des hexaedres voisins ne
516 c sont pas decoupees plus d'une fois
517 c ==> ne pas couper le quadrangle courant, ni ses aretes
518 c B. . si on n'est pas en non-conforme fidele a l'indicateur
519 c . et si au moins une des aretes des hexaedres voisins
520 c a une de ses filles a couper
521 c ==> couper toutes les faces et toutes les aretes du ou des
522 c hexaedres voisins dont une face aura une fille coupee
524 c Remarque : le cas A apparait dans le cas d'une contamination
525 c par la regle des 2 voisins
526 c le cas B apparait dans le cas d'une contamination
527 c par la regle des ecarts de niveau
529 #ifdef _DEBUG_HOMARD_
530 if ( glop.gt.0 ) then
532 write (ulsort,texte(langue,9))'.. '//mess14(langue,1,4),
533 > laface,decfac(-laface),mod( hetqua(laface) , 100 )
538 if ( tyconf.lt.3 ) then
542 lehexa = volqua(iaux,laface)
544 if ( lehexa.gt.0 ) then
545 #ifdef _DEBUG_HOMARD_
546 if ( glop.gt.0 ) then
547 write (ulsort,*)'.. hexaedre voisin : ', lehexa
551 call utarhe ( lehexa,
553 > arequa, quahex, coquhe,
556 do 421 , jaux = 1 , 12
558 larete = listar(jaux)
560 etatar = mod( hetare(larete) , 10 )
561 #ifdef _DEBUG_HOMARD_
562 if ( glop.gt.0 ) then
563 write (ulsort,texte(langue,9))'.... '//mess14(langue,1,1),
564 > larete,decare(larete),etatar
567 if ( etatar.eq.2 ) then
568 if ( decare(filare(larete)) .eq.2 .or.
569 > decare(filare(larete)+1).eq.2 ) then
583 c 4.3. ==> modification des decisions
587 c 4.3.1. ==> Cas A : ne pas couper le quadrangle courant, ni ses aretes
589 if ( kaux.eq.0 ) then
593 #ifdef _DEBUG_HOMARD_
594 if ( glop.gt.0 ) then
595 write (ulsort,texte(langue,30))'decfac',laface,decfac(-laface),' '
598 do 431 , larelo = 1 , 4
599 larete = arequa(laface,larelo)
600 #ifdef _DEBUG_HOMARD_
601 if ( glop.gt.0 ) then
602 write (ulsort,texte(langue,9))'. . '//mess14(langue,1,1),
603 > larete,decare(larete),mod( hetare(larete) , 10 )
606 if ( decare(larete).eq.2 ) then
609 #ifdef _DEBUG_HOMARD_
610 if ( glop.gt.0 ) then
611 write (ulsort,texte(langue,30))'decare', larete,decare(larete),' '
618 c 4.3.2. ==> Cas B : couper toutes les faces et toutes les aretes du ou
619 c des hexaedres voisins
623 do 432 , iaux = 1 , 2
625 lehexa = volqua(iaux,laface)
627 if ( lehexa.gt.0 ) then
628 #ifdef _DEBUG_HOMARD_
629 if ( glop.gt.0 ) then
630 write (ulsort,*)'.. hexaedre voisin : ', lehexa
633 do 4321 , jaux = 1 , 6
635 iface = quahex(lehexa,jaux)
636 etatfa = mod( hetqua(iface) , 100 )
637 if ( etatfa.ne.0 .and. decfac(-iface).eq.-1 ) then
640 #ifdef _DEBUG_HOMARD_
641 if ( glop.gt.0 ) then
642 write (ulsort,texte(langue,30))'decfac',iface,decfac(-iface),' '
645 else if ( etatfa.eq.0 .and. decfac(-iface).eq.0 ) then
648 #ifdef _DEBUG_HOMARD_
649 if ( glop.gt.0 ) then
650 write (ulsort,texte(langue,30))'decfac',iface,decfac(-iface),' '
657 call utarhe ( lehexa,
659 > arequa, quahex, coquhe,
662 do 4322 , jaux = 1 , 12
664 larete = listar(jaux)
665 etatar = mod( hetare(larete) , 10 )
666 if ( etatar.eq.2 .and. decare(larete).eq.-1 ) then
669 #ifdef _DEBUG_HOMARD_
670 if ( glop.gt.0 ) then
671 write (ulsort,texte(langue,30)) 'decare',larete,decare(larete)
674 elseif ( etatar.eq.0 .and. decare(larete).eq.0 ) then
677 #ifdef _DEBUG_HOMARD_
678 if ( glop.gt.0 ) then
679 write (ulsort,texte(langue,30)) 'decare',larete,decare(larete)
697 #ifdef _DEBUG_HOMARD_
698 write (ulsort,texte(langue,7)) 'hexa/quad'
699 if ( nbarpb.gt.0 .or. nbfapb.gt.0 ) then
700 write (ulsort,texte(langue,5)) mess14(langue,3,1), nbarpb
701 write (ulsort,texte(langue,5)) mess14(langue,3,8), nbfapb
703 write (ulsort,texte(langue,6))
708 c 5. Les suppressions de decoupage d'aretes peuvent etre nefastes
709 c pour les faces voisines. Il faut parcourir toutes les faces a
710 c couper et controler que toutes leurs aretes sont soit coupees, soit
713 #ifdef _DEBUG_HOMARD_
714 write (ulsort,90002) '5. menage ; codret', codret
717 if ( codret.eq.0 ) then
719 do 50 , laface = -nbquto, nbtrto
720 #ifdef _DEBUG_HOMARD_
721 if ( laface.eq.-215996 .or.
722 > laface.eq.20633 ) then
727 if ( glop.gt.0 ) then
728 write (ulsort,texte(langue,8)) mess14(langue,1,4),
729 > -laface,decfac(laface)
730 write (ulsort,*) ' volqua(*,laface) : ',
731 > volqua(1,-laface),volqua(2,-laface)
732 do 501 , larelo = 1 , 4
733 larete = arequa(-laface,larelo)
734 write (ulsort,texte(langue,9))'.. '//mess14(langue,1,1),
735 > larete,decare(larete),hetare(larete)
740 if ( decfac(laface).eq.4 ) then
742 if ( laface.lt.0 ) then
744 do 51 , larelo = 1 , 4
745 larete = arequa(-laface,larelo)
746 if ( decare(larete).eq.0 ) then
747 if ( mod(hetare(larete),10).eq.0 ) then
756 do 52 , larelo = 1 , 3
757 larete = aretri(laface,larelo)
758 if ( decare(larete).eq.0 ) then
759 if ( mod(hetare(larete),10).eq.0 ) then
774 #ifdef _DEBUG_HOMARD_
775 write (ulsort,texte(langue,7)) 'de coherence'
776 if ( nbarpb.gt.0 .or. nbfapb.gt.0 ) then
777 write (ulsort,texte(langue,5)) mess14(langue,3,1), nbarpb
778 write (ulsort,texte(langue,5)) mess14(langue,3,8), nbfapb
780 write (ulsort,texte(langue,6))
787 #ifdef _DEBUG_HOMARD_
788 write (ulsort,90002) '6. bilan ; codret', codret
791 if ( nbarpb.gt.0 .or. nbfapb.gt.0 ) then
794 write (ulsort,texte(langue,5)) mess14(langue,3,1), nbarpb
795 write (ulsort,texte(langue,5)) mess14(langue,3,8), nbfapb
797 #ifdef _DEBUG_HOMARD_
799 write (ulsort,texte(langue,6))
807 if ( codret.ne.0 ) then
811 write (ulsort,texte(langue,1)) 'Sortie', nompro
812 write (ulsort,texte(langue,2)) codret
816 #ifdef _DEBUG_HOMARD_
817 write (ulsort,texte(langue,1)) 'Sortie', nompro