1 subroutine pcsqug ( nbfonc, ngauss, nbnorf, typgeo, deraff,
5 > nbanqu, anfiqu, anhequ,
9 > conorf, copgrf, wipg,
10 > ulsort, langue, codret )
11 c ______________________________________________________________________
15 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
17 c Version originale enregistree le 18 juin 1996 sous le numero 96036
18 c aupres des huissiers de justice Simart et Lavoir a Clamart
19 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
20 c aupres des huissiers de justice
21 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
23 c HOMARD est une marque deposee d'Electricite de France
29 c ______________________________________________________________________
31 c aPres adaptation - Conversion de Solution -
33 c QUadrangles a plusieurs points de Gauss
35 c ______________________________________________________________________
37 c . nom . e/s . taille . description .
38 c .____________________________________________________________________.
39 c . nbfonc . e . 1 . nombre de fonctions aux points de Gauss .
40 c . ngauss . e . 1 . nbre de points de Gauss des fonctions pg .
41 c . nbnorf . e . 1 . nbre de noeuds de l'element de reference .
42 c . typgeo . e . 1 . type geometrique au sens MED .
43 c . deraff . e . 1 . vrai, s'il y a eu du deraffinement en .
44 c . . . . passant de l'iteration n a n+1 ; faux sinon.
45 c . prfcan . e . * . En numero du calcul a l'iteration n : .
46 c . . . . 0 : l'entite est absente du profil .
47 c . . . . i : l'entite est au rang i dans le profil .
48 c . prfcap . es . * . En numero du calcul a l'iteration n+1 : .
49 c . . . . 0 : l'entite est absente du profil .
50 c . . . . 1 : l'entite est presente dans le profil .
51 c . hetqua . e . nbquto . historique de l'etat des quadrangles .
52 c . filqua . e . nbquto . premier fils des quadrangles .
53 c . nbanqu . e . 1 . nombre de quadrangles decoupes par .
54 c . . . . conformite sur le maillage avant adaptation.
55 c . anfiqu . e . nbanqu . tableau filqua du maillage de l'iteration n
56 c . anhequ . e . nbanqu . tableau hetqua du maillage de l'iteration n
57 c . nqueca . e . * . nro des quadrangles dans le calcul en ent. .
58 c . nqusca . e . rsquto . numero des quadrangles du calcul .
59 c . ntreca . e . * . nro des triangles dans le calcul en entree .
60 c . ntrsca . e . rstrto . numero des triangles du calcul .
61 c . vafoen . e . nbfonc*. variables en entree de l'adaptation .
63 c . vafott . a . nbfonc*. tableau temporaire de la solution .
65 c . conorf . e . sdim* . coordonnees des noeuds de l'element de .
66 c . . . nbnorf . reference .
67 c . copgrf . e . sdim* . coordonnees des points de Gauss .
68 c . . . ngauss . de l'element de reference .
69 c . wipg . a . nbnorf*. fonctions de forme exprimees aux points de .
70 c . . . ngauss . Gauss .
71 c . ulsort . e . 1 . numero d'unite logique de la liste standard.
72 c . langue . e . 1 . langue des messages .
73 c . . . . 1 : francais, 2 : anglais .
74 c . codret . es . 1 . code de retour des modules .
75 c . . . . 0 : pas de probleme .
76 c . . . . 1 : probleme .
77 c ______________________________________________________________________
80 c 0. declarations et dimensionnement
83 c 0.1. ==> generalites
89 parameter ( nompro = 'PCSQUG' )
106 integer ngauss, nbnorf, typgeo
107 integer prfcan(*), prfcap(*)
108 integer hetqua(nbquto), ancqua(*)
109 integer filqua(nbquto)
110 integer nbanqu, anfiqu(nbanqu), anhequ(nbanqu)
111 integer nqueca(requto), nqusca(rsquto)
112 integer ntreca(retrto), ntrsca(rstrto)
114 double precision vafoen(nbfonc,ngauss,*)
115 double precision vafott(nbfonc,ngauss,*)
116 double precision conorf(sdim,nbnorf), copgrf(sdim,ngauss)
117 double precision wipg(nbnorf,ngauss)
121 integer ulsort, langue, codret
123 c 0.4. ==> variables locales
127 c qucn = QUadrangle courant en numerotation Calcul a l'it. N
128 c qucnp1 = QUadrangle courant en numerotation Calcul a l'it. N+1
129 c quhn = QUadrangle courant en numerotation Homard a l'it. N
130 c quhnp1 = QUadrangle courant en numerotation Homard a l'it. N+1
132 integer qucn, qucnp1, quhn, quhnp1
134 c f1hp = Fils 1er du quadrangle en numerotation Homard a l'it. N+1
135 c f1cp = Fils 1er du quadrangle en numerotation Calcul a l'it. N+1
136 c f2cp = Fils 2eme du quadrangle en numerotation Calcul a l'it. N+1
137 c f3cp = Fils 3eme du quadrangle en numerotation Calcul a l'it. N+1
138 c f4cp = Fils 4eme du quadrangle en numerotation Calcul a l'it. N+1
141 integer f1cp, f2cp, f3cp, f4cp
143 c f1hn = Fils 1er du quadrangle en numerotation Homard a l'it. N
144 c f1cn = Fils 1er du quadrangle en numerotation Calcul a l'it. N
145 c f2cn = Fils 2eme du quadrangle en numerotation Calcul a l'it. N
146 c f3cn = Fils 3eme du quadrangle en numerotation Calcul a l'it. N
147 c f4cn = Fils 4eme du quadrangle en numerotation Calcul a l'it. N
150 integer f1cn, f2cn, f3cn, f4cn
152 c etan = ETAt du quadrangle a l'iteration N
153 c etanp1 = ETAt du quadrangle a l'iteration N+1
157 integer nrofon, nugaus
159 double precision daux
162 parameter ( nbmess = 10 )
163 character*80 texte(nblang,nbmess)
165 c 0.5. ==> initialisations
166 c ______________________________________________________________________
175 #ifdef _DEBUG_HOMARD_
176 write (ulsort,texte(langue,1)) 'Entree', nompro
178 write(ulsort,*) 'nbfonc, ngauss, nbquto = ',nbfonc, ngauss, nbquto
182 >'(/,''Quad. en cours : nro a l''''iteration '',a3,'' :'',i10)'
184 > '( '' etat a l''''iteration '',a3,'' : '',i4)'
186 >'(/,''Quadrangle decoupe en triangles : on ne sait pas faire.'')'
189 >'(/,''Current quadrangle : # at iteration '',a3,'' : '',i10)'
191 > '( '' status at iteration '',a3,'' : '',i4)'
193 >'(/,''Quadrangle cut into triangles : not available.'')'
198 c 2. calcul des valeurs des fonctions de forme aux points de Gauss
199 c Inutile pour le moment
203 c 3. on boucle sur tous les quadrangles du maillage HOMARD n+1
204 c on trie en fonction de l'etat du quadrangle dans le maillage n
205 c on numerote les paragraphes en fonction de la documentation, a
206 c savoir : le paragraphe doc.n.p traite de la mise a jour de solution
207 c pour un quadrangle dont l'etat est passe de n a p.
208 c les autres paragraphes sont numerotes classiquement
211 #ifdef _DEBUG_HOMARD_
212 write (ulsort,*) '3. Boucle ; codret = ', codret
215 if ( codret.eq.0 ) then
217 if ( nbfonc.ne.0 ) then
219 do 30 , quhnp1 = 1 , nbquto
221 c 3.1. ==> caracteristiques du quadrangle :
222 c 3.1.1. ==> son numero homard dans le maillage precedent
225 quhn = ancqua(quhnp1)
230 c 3.1.2. ==> l'historique de son etat
231 c On rappelle que l'etat vaut :
232 c etan = 0 : le quadrangle etait actif
233 c etan = 4 : le quadrangle etait coupe en 4 ; il y a eu
235 c etan = 55 : le quadrangle n'existait pas ; il a ete produit
237 c etan = 31, 32, 33, 34 : le quadrangle etait coupe en 3
238 c quadrangles ; il y a eu deraffinement.
240 etanp1 = mod(hetqua(quhnp1),100)
241 etan = (hetqua(quhnp1)-etanp1) / 100
243 cgn write (ulsort,1792) 'Quadrangle', quhn, etan, quhnp1, etanp1
245 c=======================================================================
246 c doc.0.p. ==> etan = 0 : le quadrangle etait actif
247 c=======================================================================
249 if ( etan.eq.0 ) then
251 c on repere son ancien numero dans le calcul
255 if ( prfcan(qucn).gt.0 ) then
257 cgn print 1789,(vafoen(nrofon,nugaus,qucn), nugaus = 1 , ngauss)
258 cgn 1789 format(' Valeurs anciennes : ',5g12.5)
260 c doc.0.0. ===> etanp1 = 0 : le quadrangle etait actif et l'est encore ;
262 c c'est le cas le plus simple : on prend la valeur de la
263 c fonction associee a l'ancien numero du quadrangle.
265 c ................. .................
273 c ................. .................
275 if ( etanp1.eq.0 ) then
277 qucnp1 = nqusca(quhnp1)
280 do 321 , nrofon = 1 , nbfonc
281 cgn write(ulsort,7778)
282 cgn > (vafoen(nrofon,nugaus,prfcan(qucn)),nugaus=1,ngauss)
283 do 3211 , nugaus = 1 , ngauss
284 vafott(nrofon,nugaus,qucnp1) =
285 > vafoen(nrofon,nugaus,prfcan(qucn))
288 cgn write(31,7777) qucnp1
289 cgn write(ulsort,7777) qucn,-1,qucnp1
291 cgn7778 format(8g14.7)
293 c doc.0.1/2/3 ==> etanp1 = 31, 32, 33 ou 34 : le quadrangle etait actif
294 c et est decoupe en 3 triangles.
295 c les trois fils prennent la valeur de la fonction sur
297 c ................. .................
305 c ................. .................
307 elseif ( etanp1.ge.31 .and. etanp1.le.34 ) then
309 f1hp = -filqua(quhnp1)
311 f2cp = nqusca(f1hp+1)
312 f3cp = nqusca(f1hp+2)
316 write (ulsort,texte(langue,6))
319 c doc.0.4/6/7/8. ==> etanp1 = 4 : le quadrangle etait actif et
321 c les quaque fils prennent la valeur de la fonction
323 c ................. .................
327 c . . ===> .................
331 c ................. .................
333 elseif ( etanp1.eq.4 ) then
335 f1hp = filqua(quhnp1)
337 f2cp = nqusca(f1hp+1)
338 f3cp = nqusca(f1hp+2)
339 f4cp = nqusca(f1hp+3)
344 do 323 , nrofon = 1 , nbfonc
345 do 3231 , nugaus = 1 , ngauss
346 daux = vafoen(nrofon,nugaus,prfcan(qucn))
347 vafott(nrofon,nugaus,f1cp) = daux
348 vafott(nrofon,nugaus,f2cp) = daux
349 vafott(nrofon,nugaus,f3cp) = daux
350 vafott(nrofon,nugaus,f4cp) = daux
354 c doc.0.erreur. ==> aucun autre etat sur le quadrangle courant n'est
360 write (ulsort,texte(langue,4)) 'n ', quhn
361 write (ulsort,texte(langue,5)) 'n ', etan
362 write (ulsort,texte(langue,4)) 'n+1', quhnp1
363 write (ulsort,texte(langue,5)) 'n+1', etanp1
369 c=======================================================================
370 c doc.1/2/3.p. ==> etan = 31, 32, 33 ou 34 : le quadrangle etait coupe
372 c=======================================================================
374 elseif ( etan.ge.31 .and. etan.le.34 ) then
376 c on repere les numeros dans le calcul pour ses trois fils a
381 f2cn = ntreca(f1hn+1)
382 f3cn = ntreca(f1hn+2)
384 if ( prfcan(f1cn).gt.0 .and. prfcan(f2cn).gt.0 .and.
385 > prfcan(f3cn).gt.0 ) then
387 write (ulsort,texte(langue,6))
390 c doc.1/2/3.0. ===> etanp1 = 0 : le quadrangle est actif. il est
392 c on lui attribue la valeur moyenne sur les trois
394 c remarque : cela arrive seulement avec du deraffinement.
395 c ................. .................
403 c ................. .................
405 if ( etanp1.eq.0 ) then
407 qucnp1 = nqusca(quhnp1)
410 c doc.1/2/3.1/2/3. ===> etanp1 = etan : le quadrangle est decoupe en
411 c trois selon le meme decoupage.
412 c c'est ce qui se passe quand un decoupage de conformite
413 c est supprime au debut des algorithmes d'adaptation,
414 c puis reproduit a la creation du maillage car les
415 c quadrangles autour n'ont pas change entre les deux
417 c le fils prend la valeur de la fonction sur l'ancien
418 c fils qui etait au meme endroit. comme la procedure de
419 c numerotation est la meme (voir cmcdqu), le premier fils
420 c est toujours le meme, le 2eme et le 3eme egalement.
421 c on prendra alors la valeur sur le fils de rang identique
423 c ................. .................
427 c . . . . ===> . . . .
431 c ................. .................
433 elseif ( etanp1.eq.etan ) then
435 f1hp = -filqua(quhnp1)
437 f2cp = ntrsca(f1hp+1)
438 f3cp = ntrsca(f1hp+2)
443 c doc.1/2/3.perm(1/2/3). ===> etanp1 = 31, 32, 33 ou 34 et different de
444 c etan : le quadrangle est encore decoupe
445 c en trois, mais par un autre decoupage.
446 c c'est ce qui se passe quand un decoupage de conformite
447 c est supprime au debut des algorithmes d'adaptation. il y
448 c a du deraffinement dans la zone qui induisait le decoupage
449 c de conformite et raffinement sur une autre zone.
450 c On donne la valeur moyenne de la fonction sur les trois
451 c anciens fils a chaque nouveau fils.
452 c remarque : on pourrait certainement faire mieux, avec des
453 c moyennes ponderees en fonction du recouvrement
454 c des anciens et nouveaux fils. c'est trop
455 c complique pour que cela vaille le coup.
456 c remarque : cela arrive seulement avec du deraffinement.
458 c ................. .................
462 c . . . . ===> . . . .
466 c ................. .................
468 elseif ( etanp1.ge.31 .and. etanp1.le.34 ) then
470 f1hp = -filqua(quhnp1)
472 f2cp = ntrsca(f1hp+1)
473 f3cp = ntrsca(f1hp+2)
478 c doc.1/2/3.4/6/7/8. ===> etanp1 = 4 : le quadrangle est
480 c c'est ce qui se passe quand un decoupage de conformite
481 c est supprime au debut des algorithmes d'adaptation. il y
482 c a ensuite raffinement du quadrangle. qui plus est, par
483 c suite de la regle des ecarts de niveau, on peut avoir
484 c induit un decoupage de conformite sur l'un des fils.
486 c On donne la valeur moyenne de la fonction sur les trois
487 c anciens fils a chaque nouveau fils.
488 c remarque : on pourrait certainement faire mieux, avec des
489 c moyennes ponderees en fonction du recouvrement
490 c des anciens et nouveaux fils. c'est trop
491 c complique pour que cela vaille le coup.
493 c ................. .................
497 c . . . . ===> .................
501 c ................. .................
504 elseif ( etanp1.eq.4 ) then
506 f1hp = filqua(quhnp1)
508 f2cp = nqusca(f1hp+1)
509 f3cp = nqusca(f1hp+2)
510 f4cp = nqusca(f1hp+3)
516 c doc.1/2/3.erreur. ==> aucun autre etat sur le quadrangle courant
522 write (ulsort,texte(langue,4)) 'n ', quhn
523 write (ulsort,texte(langue,5)) 'n ', etan
524 write (ulsort,texte(langue,4)) 'n+1', quhnp1
525 write (ulsort,texte(langue,5)) 'n+1', etanp1
531 c=======================================================================
532 c doc.4. ==> le quadrangle etait coupe en 4 :
533 c=======================================================================
535 elseif ( etan.eq.4 ) then
537 c on repere les numeros dans le calcul pour ses quatre fils
542 f2cn = nqueca(f1hn+1)
543 f3cn = nqueca(f1hn+2)
544 f4cn = nqueca(f1hn+3)
546 if ( prfcan(f1cn).gt.0 .and. prfcan(f2cn).gt.0 .and.
547 > prfcan(f3cn).gt.0 .and. prfcan(f4cn).gt.0 ) then
549 c doc.4.0. ===> etanp1 = 0 : le quadrangle est actif ; il est reactive.
550 c on lui attribue la valeur moyenne sur les quatre anciens
552 c remarque : cela arrive seulement avec du deraffinement.
553 c ................. .................
557 c ................. ===> . .
561 c ................. .................
563 if ( etanp1.eq.0 ) then
565 qucnp1 = nqusca(quhnp1)
568 do 341 , nrofon = 1 , nbfonc
569 do 3411 , nugaus = 1 , ngauss
570 vafott(nrofon,nugaus,qucnp1) =
571 > unsqu * ( vafoen(nrofon,nugaus,prfcan(f1cn))
572 > + vafoen(nrofon,nugaus,prfcan(f2cn))
573 > + vafoen(nrofon,nugaus,prfcan(f3cn))
574 > + vafoen(nrofon,nugaus,prfcan(f4cn)) )
578 c doc.4.1/2/3. ===> etanp1 = 31, 32, 33 ou 34 : le quadrangle est
580 c on attribue la valeur moyenne sur les quatre anciens
581 c fils a chacune des trois nouveaux fils.
582 c remarque : on pourrait certainement faire mieux, avec des
583 c moyennes ponderees en fonction du recouvrement
584 c des anciens et nouveaux fils. c'est trop
585 c complique pour que cela vaille le coup.
586 c remarque : cela arrive seulement avec du deraffinement.
587 c ................. .................
591 c ................. ===> . . . .
595 c ................. .................
597 elseif ( etanp1.ge.31 .and. etanp1.le.34 ) then
599 f1hp = -filqua(quhnp1)
601 f2cp = ntrsca(f1hp+1)
602 f3cp = ntrsca(f1hp+2)
606 write (ulsort,texte(langue,6))
613 c=======================================================================
627 cgn do 922 , nrofon = 1 , nbfonc
628 cgn print *,'fonction numero ', nrofon
630 cgn do 9222 , quhnp1 = 1 , nbtrto
631 cgn if ( mod(hettri(quhnp1),100).eq.0 ) then
633 cgn print 1788,quhnp1,
634 cgn > (vafott(nrofon,nugaus,iaux), nugaus = 1 , ngauss)
638 if ( codret.ne.0 ) then
642 write (ulsort,texte(langue,1)) 'Sortie', nompro
643 write (ulsort,texte(langue,2)) codret
647 #ifdef _DEBUG_HOMARD_
648 write (ulsort,texte(langue,1)) 'Sortie', nompro