1 subroutine derco2 ( tyconf, niveau,
4 > hettri, aretri, filtri, nivtri,
6 > hetqua, arequa, filqua, nivqua,
12 > ulsort, langue, codret )
13 c ______________________________________________________________________
17 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
19 c Version originale enregistree le 18 juin 1996 sous le numero 96036
20 c aupres des huissiers de justice Simart et Lavoir a Clamart
21 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
22 c aupres des huissiers de justice
23 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
25 c HOMARD est une marque deposee d'Electricite de France
31 c ______________________________________________________________________
33 c traitement des DEcisions - Raffinement : COntamination - option 2
35 c Application de la regle des ecarts de niveau, tout type de raffinement
36 c ______________________________________________________________________
38 c . nom . e/s . taille . description .
39 c .____________________________________________________________________.
40 c . tyconf . e . 1 . 0 : conforme (defaut) .
41 c . . . . 1 : non-conforme avec au minimum 2 aretes .
42 c . . . . non decoupees en 2 .
43 c . . . . 2 : non-conforme avec 1 seul noeud .
44 c . . . . pendant par arete .
45 c . . . . 3 : non-conforme fidele a l'indicateur .
46 c . . . . -1 : conforme, avec des boites pour les .
47 c . . . . quadrangles, hexaedres et pentaedres .
48 c . . . . -2 : non-conforme avec au maximum 1 arete .
49 c . . . . decoupee en 2 (boite pour les .
50 c . . . . quadrangles, hexaedres et pentaedres) .
51 c . . . . -2 : non-conforme avec au maximum 1 arete .
52 c . . . . decoupee en 2 (boite pour les .
53 c . . . . quadrangles, hexaedres et pentaedres) .
54 c . niveau . e . 1 . niveau en cours d'examen .
55 c . decare . es . nbarto . decisions des aretes .
56 c . decfac . es . -nbquto. decision sur les faces (quad. + tri.) .
58 c . hetare . e . nbarto . historique de l'etat des aretes .
59 c . filare . e . nbarto . premiere fille des aretes .
60 c . hettri . e . nbtrto . historique de l'etat des triangles .
61 c . aretri . e . nbtrto . numeros des 3 aretes des triangles .
62 c . filtri . e . nbtrto . premier fils des triangles .
63 c . nivtri . e . nbtrto . niveau des triangles .
64 c . voltri . e .2*nbtrto. numeros des 2 volumes par triangle .
65 c . . . . voltri(i,k) definit le i-eme voisin de k .
66 c . . . . 0 : pas de voisin .
67 c . . . . j>0 : tetraedre j .
68 c . . . . j<0 : pyramide/pentaedre dans pypetr(1/2,j).
69 c . pypetr . e .2*lgpype. pypetr(1,j) = numero de la pyramide voisine.
70 c . . . . du triangle k tel que voltri(1/2,k) = -j .
71 c . . . . pypetr(2,j) = numero du pentaedre voisin .
72 c . . . . du triangle k tel que voltri(1/2,k) = -j .
73 c . hetqua . e . nbquto . historique de l'etat des quadrangles .
74 c . arequa . e . nbquto . numeros des 4 aretes des quadrangles .
75 c . filqua . e . nbquto . premier fils des quadrangles .
76 c . nivqua . e . nbquto . niveau des quadrangles .
77 c . volqua . e .2*nbquto. numeros des 2 volumes par quadrangle .
78 c . . . . volqua(i,k) definit le i-eme voisin de k .
79 c . . . . 0 : pas de voisin .
80 c . . . . j>0 : hexaedre j .
81 c . . . . j<0 : pyramide/pentaedre dans pypequ(1/2,j).
82 c . pypequ . e .2*lgpype. pypequ(1,j) = numero de la pyramide voisine.
83 c . . . . du quadrangle k tel que volqua(1/2,k) = -j .
84 c . . . . pypequ(2,j) = numero du pentaedre voisin .
85 c . . . . du quadrangle k tel que volqua(1/2,k) = -j .
86 c . tritet . e .nbtecf*4. numeros des 4 triangles des tetraedres .
87 c . quahex . e .nbhecf*6. numeros des 6 quadrangles des hexaedres .
88 c . coquhe . e .nbhecf*6. codes des 6 quadrangles des hexaedres .
89 c . facpyr . e .nbpycf*5. numeros des 5 faces des pyramides .
90 c . cofapy . e .nbpycf*5. codes des faces des pyramides .
91 c . facpen . e .nbpecf*5. numeros des faces des pentaedres .
92 c . cofape . e .nbpecf*5. code des 5 faces des pentaedres .
93 c . ulsort . e . 1 . numero d'unite logique de la liste standard.
94 c . langue . e . 1 . langue des messages .
95 c . . . . 1 : francais, 2 : anglais .
96 c . codret . es . 1 . code de retour des modules .
97 c . . . . 0 : pas de probleme .
98 c ______________________________________________________________________
101 c 0. declarations et dimensionnement
104 c 0.1. ==> generalites
110 parameter ( nompro = 'DERCO2' )
129 integer decare(0:nbarto)
130 integer decfac(-nbquto:nbtrto)
131 integer hetare(nbarto), filare(nbarto)
132 integer hettri(nbtrto), aretri(nbtrto,3), filtri(nbtrto)
133 integer nivtri(nbtrto)
134 integer voltri(2,nbtrto), pypetr(2,*)
135 integer hetqua(nbquto), arequa(nbquto,4), filqua(nbquto)
136 integer nivqua(nbquto)
137 integer volqua(2,nbquto), pypequ(2,*)
138 integer tritet(nbtecf,4)
139 integer quahex(nbhecf,6), coquhe(nbhecf,6)
140 integer facpyr(nbpycf,5), cofapy(nbpycf,5)
141 integer facpen(nbpecf,5), cofape(nbpecf,5)
143 integer ulsort, langue, codret
145 c 0.4. ==> variables locales
147 integer iaux, jaux, kaux
149 integer laface, lafac2
152 integer nbaret, listar(12), nbface, listfa(12)
153 integer nbvolu, listvo(2), typevo(2)
154 integer nbvotr, nbvoqu, nbvoto
161 parameter ( nbmess = 30 )
162 character*80 texte(nblang,nbmess)
164 c 0.5. ==> initialisations
165 c ______________________________________________________________________
173 #ifdef _DEBUG_HOMARD_
174 write (ulsort,texte(langue,1)) 'Entree', nompro
182 #ifdef _DEBUG_HOMARD_
183 write (ulsort,texte(langue,12)) niveau
185 #ifdef _DEBUG_HOMARD_
186 write (ulsort,90002) 'tyconf', tyconf
189 #ifdef _DEBUG_HOMARD_
190 write (ulsort,*) 'entree de ',nompro
191 do 1105 , iaux = 1 , -nbquto
192 write (ulsort,90001) 'decision quadrangle', iaux,decfac(-iaux)
193 cgn write (ulsort,90001) 'quadrangle', iaux,
194 cgn > arequa(iaux,1), arequa(iaux,2),
195 cgn > arequa(iaux,3), arequa(iaux,4)
198 #ifdef _DEBUG_HOMARD_
199 if ( nbquto.gt.0 ) then
200 iaux = min(nbquto,38)
201 write (ulsort,90001) 'quadrangle', iaux,
202 > arequa(iaux,1), arequa(iaux,2),
203 > arequa(iaux,3), arequa(iaux,4)
204 write (ulsort,90001) 'quadrangle', iaux,
205 > decare(arequa(iaux,1)), decare(arequa(iaux,2)),
206 > decare(arequa(iaux,3)), decare(arequa(iaux,4)), decfac(-iaux)
209 #ifdef _DEBUG_HOMARD_
210 if ( nbquto.gt.0 ) then
211 iaux = min(nbquto,10)
212 write (ulsort,90001) 'quadrangle', iaux,
213 > arequa(iaux,1), arequa(iaux,2),
214 > arequa(iaux,3), arequa(iaux,4)
215 write (ulsort,90001) 'quadrangle', iaux,
216 > decare(arequa(iaux,1)), decare(arequa(iaux,2)),
217 > decare(arequa(iaux,3)), decare(arequa(iaux,4)), decfac(-iaux)
218 iaux = min(nbquto,19)
219 write (ulsort,90001) 'quadrangle', iaux,
220 > arequa(iaux,1), arequa(iaux,2),
221 > arequa(iaux,3), arequa(iaux,4)
222 write (ulsort,90001) 'quadrangle', iaux,
223 > decare(arequa(iaux,1)), decare(arequa(iaux,2)),
224 > decare(arequa(iaux,3)), decare(arequa(iaux,4)), decfac(-iaux)
230 nbvoto = nbteto + nbpyto + nbheto + nbpeto
232 c nombre maximum de volumes par triangle ou quadrangle
234 if ( nbteto.eq.0 .and. nbpyto.eq.0 .and. nbpeto.eq.0 ) then
240 if ( nbheto.eq.0 .and. nbpyto.eq.0 .and. nbpeto.eq.0 ) then
247 c 2. Application de la regle des ecarts de niveau aux faces
250 do 2 , laface = -nbquto , nbtrto
251 cgn print *,'debut boucle 2 : decfac(',laface,') :',decfac(laface)
253 c 2.1. ==> On s'interesse aux faces du niveau courant :
255 c ou . inactives a garder et bord de volume
256 c ou . inactives a reactiver
258 if ( laface.gt.0 ) then
260 if ( nivtri(laface).eq.niveau ) then
261 etat = mod( hettri(laface) , 10 )
262 cgn write (ulsort,texte(langue,29))'Triangle', laface,
263 cgn > nivtri(laface), hettri(laface), decfac(laface)
268 elseif ( laface.lt.0 ) then
271 if ( nivqua(iaux).eq.niveau ) then
272 etat = mod( hetqua(iaux) , 100 )
273 cgn write (ulsort,texte(langue,29))'Quadrangle', -laface,
274 cgn > nivqua(-laface), hetqua(-laface), decfac(laface)
286 if ( etat.eq.0 ) then
287 if ( decfac(laface).eq.0 ) then
290 elseif ( etat.eq.4 .or.
291 > etat.eq.6 .or. etat.eq.7 .or. etat.eq.8 ) then
292 if ( decfac(laface).eq.0 .and. nbvoto.gt.0 ) then
294 elseif ( decfac(laface).eq.-1 ) then
295 if ( nbvoto.eq.0 ) then
302 cgn write (ulsort,*) 'Face', laface, ', choix = ', choix
304 c 2.2. ==> Liste des aretes de la face
306 if ( choix.gt.0 ) then
308 if ( laface.gt.0 ) then
310 do 221 , iarelo = 1 , nbaret
311 listar(iarelo) = aretri(laface,iarelo)
316 do 222 , iarelo = 1 , nbaret
317 listar(iarelo) = arequa(iaux,iarelo)
321 #ifdef _DEBUG_HOMARD_
323 write (ulsort,texte(langue,15))
327 c 2.3. ==> Cas du raffinement a propager par voisinage
329 if ( choix.eq.1 ) then
331 c 2.3.1. ==> Decompte des aretes coupees en 2 avec une fille a couper :
332 c . celles d'etat > 0
333 c . et avec une fille de decision > 0
334 c S'il n'y en a pas, rien n'est a faire
337 do 231, iaux = 1 , nbaret
338 larete = listar(iaux)
339 cgn if ( larete.eq.1661 ) then
340 cgn write(ulsort,90002) '.... arete possible', larete
342 if ( mod(hetare(larete),10).gt.0 ) then
343 jaux = filare(larete)
344 if ( decare(jaux).gt.0 .or. decare(jaux+1).gt.0 ) then
352 c 2.3.2. ==> Propagation du raffinement sur la face et ses
358 #ifdef _DEBUG_HOMARD_
359 write (ulsort,texte(langue,30))'decfac', laface,decfac(laface),' '
361 do 232 , iaux = 1 , nbaret
362 larete = listar(iaux)
363 if ( decare(larete).eq.0 ) then
364 if ( mod(hetare(larete),10).eq.0 ) then
366 #ifdef _DEBUG_HOMARD_
367 if ( larete.eq.1661 ) then
368 write (ulsort,texte(langue,30))' decare',larete,decare(larete),' '
379 c 2.4. ==> Cas du deraffinement a inhiber par voisinage
381 if ( choix.eq.2 .or. choix.eq.4 ) then
383 c 2.4.1. ==> Existe-t-il des aretes coupees en 2 avec une fille coupee
384 c qui doit etre coupee ?
387 do 241, iaux = 1 , nbaret
388 larete = listar(iaux)
389 cgn write (ulsort,*) larete, decare(larete)
390 jaux = filare(larete)
391 if ( decare(jaux).gt.0 .or. decare(jaux+1).gt.0 ) then
398 c 2.4.2. ==> Inhibition du raffinement sur la face et ses aretes
403 #ifdef _DEBUG_HOMARD_
404 write (ulsort,texte(langue,30))'decfac', laface,decfac(laface),' '
406 do 242 , iaux = 1 , nbaret
407 larete = listar(iaux)
408 if ( decare(larete).eq.-1 ) then
410 #ifdef _DEBUG_HOMARD_
411 write (ulsort,texte(langue,30))' decare',larete,decare(larete),' '
420 c 2.5. ==> Cas du raffinement a propager ou du deraffinement a inhiber
421 c par l'interieur de volumes
423 if ( choix.ge.3 ) then
425 c 2.5.1. ==> Pour chaque face, on regarde si une arete tracee sur
426 c la face est coupee ou va etre coupee.
427 c . Pour un triangle, ces aretes sont celles qui definissent
428 c la fille face centrale (cf. cmrdtr)
429 c . Pour un quadrangle, ces aretes sont la 2eme et le 3eme
430 c du premier et du troisieme fils (cf. cmrdqu)
431 c S'il n'y en a pas, rien n'est a faire
433 if ( laface.gt.0 ) then
434 jaux = filtri(laface)
436 do 2511 , iarelo = 1 , nbaret
437 listar(iarelo) = aretri(jaux,iarelo)
440 jaux = filqua(-laface)
442 listar(1) = arequa(jaux ,2)
443 listar(2) = arequa(jaux ,3)
444 listar(3) = arequa(jaux+2,2)
445 listar(4) = arequa(jaux+2,3)
449 do 2513 , iarelo = 1 , nbaret
450 cgn write (ulsort,*) 'hetare, decare(',listar(iarelo),') =',
451 cgn >hetare(listar(iarelo)), decare(listar(iarelo))
452 if ( decare(listar(iarelo)).gt.0 .or.
453 > mod(hetare(listar(iarelo)),10).eq.2 ) then
459 cgn write (ulsort,99001) 'afaire', afaire
461 c 2.5.2. ==> La face retenue borde-t-elle un volume ?
467 if ( laface.gt.0 ) then
469 do 2521, iaux = 1 , nbvotr
470 jaux = voltri(iaux,laface)
471 if ( jaux.gt.0 ) then
473 listvo(nbvolu) = jaux
475 elseif ( jaux.lt.0 ) then
476 if ( pypetr(1,-jaux).ne.0 ) then
478 listvo(nbvolu) = pypetr(1,-jaux)
481 if ( pypetr(2,-jaux).ne.0 ) then
483 listvo(nbvolu) = pypetr(2,-jaux)
491 do 2522, iaux = 1 , nbvoqu
492 jaux = volqua(iaux,-laface)
493 if ( jaux.gt.0 ) then
495 listvo(nbvolu) = jaux
497 elseif ( jaux.lt.0 ) then
498 if ( pypequ(1,-jaux).ne.0 ) then
500 listvo(nbvolu) = pypequ(1,-jaux)
503 if ( pypequ(2,-jaux).ne.0 ) then
505 listvo(nbvolu) = pypequ(2,-jaux)
512 cgn write (ulsort,*)nbvolu,'volumes', (listvo(iaux),iaux=1,nbvolu)
513 cgn write (ulsort,*)nbvolu,'types ', (typevo(iaux),iaux=1,nbvolu)
517 c 2.5.3. ==> Une des aretes tracees sur laface sera coupee. Il faut que
518 c le ou les volumes s'appuyant sur laface soient coupes
520 if ( nbvolu.gt.0 ) then
522 c 2.5.3. ==> Recherche des faces concernees
525 do 2531 , iaux = 1 , nbvolu
527 cgn write (ulsort,*)'Volume', jaux,' de type',typevo(iaux)
528 if ( typevo(iaux).eq.3 ) then
529 do 25311 , kaux = 1 , 4
531 listfa(nbface) = tritet(jaux,kaux)
533 elseif ( typevo(iaux).eq.5 ) then
534 listfa(1) = facpyr(jaux,1)
535 listfa(2) = facpyr(jaux,2)
536 listfa(3) = facpyr(jaux,3)
537 listfa(4) = facpyr(jaux,4)
538 listfa(5) = -facpyr(jaux,5)
540 elseif ( typevo(iaux).eq.6 ) then
541 do 25313 , kaux = 1 , 6
543 listfa(nbface) = -quahex(jaux,kaux)
545 elseif ( typevo(iaux).eq.7 ) then
546 listfa(1) = facpen(jaux,1)
547 listfa(2) = facpen(jaux,2)
548 listfa(3) = -facpen(jaux,3)
549 listfa(4) = -facpen(jaux,4)
550 listfa(5) = -facpen(jaux,5)
554 cgn write (ulsort,1000)nbface,' faces :',
555 cgn > (listfa(iaux),iaux=1,nbface)
556 cgn 1000 format(i2,a,12i5)
558 do 2532 , iaux = 1 , nbface
560 lafac2 = listfa(iaux)
561 cgn if ( lafac2.gt.0 ) then
562 cgn write (ulsort,texte(langue,29))'Triangle', lafac2,
563 cgn > nivtri(lafac2), hettri(lafac2), decfac(lafac2)
565 cgn write (ulsort,texte(langue,29))'Quadrangle', -lafac2,
566 cgn > nivqua(-lafac2), hetqua(-lafac2), decfac(lafac2)
568 if ( decfac(lafac2).eq.-1 ) then
570 #ifdef _DEBUG_HOMARD_
571 write (ulsort,texte(langue,30))'decfac', lafac2,decfac(lafac2),' '
573 elseif ( decfac(lafac2).eq.0 ) then
574 if ( lafac2.gt.0 ) then
575 if ( mod(hettri(lafac2),10).eq.0 ) then
577 #ifdef _DEBUG_HOMARD_
578 write (ulsort,texte(langue,30))'decfac', lafac2,decfac(lafac2),' '
582 if ( mod(hetqua(-lafac2),100).eq.0 ) then
584 #ifdef _DEBUG_HOMARD_
585 write (ulsort,texte(langue,30))'decfac', lafac2,decfac(lafac2),' '
591 if ( lafac2.gt.0 ) then
593 do 2533 , iarelo = 1 , nbaret
594 listar(iarelo) = aretri(lafac2,iarelo)
598 do 2534 , iarelo = 1 , nbaret
599 listar(iarelo) = arequa(-lafac2,iarelo)
603 do 2535 , jaux = 1 , nbaret
604 larete = listar(jaux)
605 if ( decare(larete).eq.0 ) then
606 if ( mod(hetare(larete),10).eq.0 ) then
608 #ifdef _DEBUG_HOMARD_
609 if ( larete.eq.1661)then
610 write (ulsort,texte(langue,30))' decare',larete,decare(larete),' '
614 elseif ( decare(larete).eq.-1 ) then
616 #ifdef _DEBUG_HOMARD_
617 write (ulsort,texte(langue,30))' decare',larete,decare(larete),' '
631 c 3. Transfert via les volumes ayant des quadrangles comme faces
632 c Si une fille de l'une de ses aretes est a couper, le volume
633 c doit l'etre entierement : on le declare par ses aretes.
634 c 1/12/16 : c'est trop. Le decoupage est assure par l'etape 2.3.
636 #ifdef _DEBUG_HOMARD_
637 write(ulsort,90002) '3. Transfert ; codret', codret
639 #ifdef _DEBUG_HOMARD_
640 if ( nbquto.gt.0 ) then
641 iaux = min(nbquto,38)
642 write (ulsort,90001) 'quadrangle', iaux,
643 > arequa(iaux,1), arequa(iaux,2),
644 > arequa(iaux,3), arequa(iaux,4)
645 write (ulsort,90001) 'quadrangle', iaux,
646 > decare(arequa(iaux,1)), decare(arequa(iaux,2)),
647 > decare(arequa(iaux,3)), decare(arequa(iaux,4)), decfac(-iaux)
651 if ( tyconf.eq.1789 ) then
653 if ( codret.eq.0 ) then
655 write (ulsort,*) 'ATTENTION'
656 write (ulsort,texte(langue,3)) 'DERCO9', nompro
657 call derco9 ( niveau,
665 > ulsort, langue, codret )
671 #ifdef _DEBUG_HOMARD_
672 write (ulsort,*) 'sortie de ',nompro
673 do 11060 , iaux = 1 , nbarto
674 if ( iaux.eq.-17735 .or. iaux.le.-877 ) then
675 write (ulsort,90001) '.. arete e/d', iaux,
676 > hetare(iaux), decare(iaux)
680 #ifdef _DEBUG_HOMARD_
681 do 1106 , iaux = 1 , nbquto
682 write (ulsort,90001) 'decision quadrangle', iaux,decfac(-iaux)
683 cgn write (ulsort,90001) 'quadrangle', iaux,
684 cgn > arequa(iaux,1), arequa(iaux,2),
685 cgn > arequa(iaux,3), arequa(iaux,4)
688 #ifdef _DEBUG_HOMARD_
689 if ( nbquto.gt.0 ) then
690 iaux = min(nbquto,12)
691 write (ulsort,90001) 'quadrangle', iaux,
692 > arequa(iaux,1), arequa(iaux,2),
693 > arequa(iaux,3), arequa(iaux,4)
694 write (ulsort,90001) 'quadrangle', iaux,
695 > decare(arequa(iaux,1)), decare(arequa(iaux,2)),
696 > decare(arequa(iaux,3)), decare(arequa(iaux,4)), decfac(-iaux)
697 iaux = min(nbquto,10)
698 write (ulsort,90001) 'quadrangle', iaux,
699 > arequa(iaux,1), arequa(iaux,2),
700 > arequa(iaux,3), arequa(iaux,4)
701 write (ulsort,90001) 'quadrangle', iaux,
702 > decare(arequa(iaux,1)), decare(arequa(iaux,2)),
703 > decare(arequa(iaux,3)), decare(arequa(iaux,4)), decfac(-iaux)
704 iaux = min(nbquto,19)
705 write (ulsort,90001) 'quadrangle', iaux,
706 > arequa(iaux,1), arequa(iaux,2),
707 > arequa(iaux,3), arequa(iaux,4)
708 write (ulsort,90001) 'quadrangle', iaux,
709 > decare(arequa(iaux,1)), decare(arequa(iaux,2)),
710 > decare(arequa(iaux,3)), decare(arequa(iaux,4)), decfac(-iaux)
713 #ifdef _DEBUG_HOMARD_
714 if ( nbquto.gt.0 ) then
715 iaux = min(nbquto,38)
716 write (ulsort,90001) 'quadrangle', iaux,
717 > arequa(iaux,1), arequa(iaux,2),
718 > arequa(iaux,3), arequa(iaux,4)
719 write (ulsort,90001) 'quadrangle', iaux,
720 > decare(arequa(iaux,1)), decare(arequa(iaux,2)),
721 > decare(arequa(iaux,3)), decare(arequa(iaux,4)), decfac(-iaux)
729 if ( codret.ne.0 ) then
733 write (ulsort,texte(langue,1)) 'Sortie', nompro
734 write (ulsort,texte(langue,2)) codret
738 #ifdef _DEBUG_HOMARD_
739 write (ulsort,texte(langue,1)) 'Sortie', nompro