1 subroutine deinse ( typenh,
4 > typseh, typseb, seuilh, seuilb, nbsoci,
6 > ulsort, langue, codret)
7 c ______________________________________________________________________
11 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
13 c Version originale enregistree le 18 juin 1996 sous le numero 96036
14 c aupres des huissiers de justice Simart et Lavoir a Clamart
15 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
16 c aupres des huissiers de justice
17 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
19 c HOMARD est une marque deposee d'Electricite de France
25 c ______________________________________________________________________
27 c traitement des DEcisions - INitialisation des SEuils
29 c ______________________________________________________________________
31 c . nom . e/s . taille . description .
32 c .____________________________________________________________________.
33 c . typenh . e . 1 . type d'entites concernees .
34 c . . . . -1 : noeuds .
35 c . . . . 1 : aretes .
36 c . . . . 2 : triangles .
37 c . . . . 3 : tetraedres .
38 c . . . . 4 : quadrangles .
39 c . . . . 5 : pyramides .
40 c . . . . 6 : hexaedres .
41 c . . . . 7 : pentaedres .
42 c . seuihe . s . 1 . borne superieure absolue de l'erreur entite.
43 c . seuibe . s . 1 . borne inferieure absolue de l'erreur entite.
44 c . pilraf . e . 1 . pilotage du raffinement .
45 c . . . . -1 : raffinement uniforme .
46 c . . . . 0 : pas de raffinement .
47 c . . . . 1 : raffinement libre .
48 c . . . . 2 : raff. libre homogene en type d'element.
49 c . pilder . e . 1 . pilotage du deraffinement .
50 c . . . . 0 : pas de deraffinement .
51 c . . . . 1 : deraffinement libre .
52 c . . . . -1 : deraffinement uniforme .
53 c . typseh . e . 1 . type de seuil haut .
54 c . . . . 1 : absolu .
55 c . . . . 2 : relatif .
56 c . . . . 3 : pourcentage d'entites .
57 c . . . . 4 : moyenne + nh*ecart-type .
58 c . . . . 5 : cible en nombre de noeuds .
59 c . typseb . e . 1 . type de seuil bas .
60 c . . . . 1 : absolu .
61 c . . . . 2 : relatif .
62 c . . . . 3 : pourcentage d'entites .
63 c . . . . 4 : moyenne - nb*ecart-type .
64 c . seuilh . e . 1 . borne superieure de l'erreur (absolue, .
65 c . . . . relatif, pourcentage d'entites ou nh) .
66 c . seuilb . e . 1 . borne inferieure de l'erreur (absolue, .
67 c . . . . relatif, pourcentage d'entites ou nb) .
68 c . nbsoci . e . 1 . cible en nombre de sommets (-1 si non) .
69 c . indtab . e . 1 . dernier indice affecte dans tabind .
70 c . tabind . e . indtab . tableau de l'indicateur .
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 . . . . 4 : nombres d'entites incoherents .
77 c . . . . 2 : probleme dans le traitement .
78 c . . . . 3 : les seuils sont mal definis .
79 c . . . . 5 : mauvaise cible .
80 c ______________________________________________________________________
83 c 0. declarations et dimensionnement
86 c 0.1. ==> generalites
92 parameter ( nompro = 'DEINSE' )
112 integer pilraf, pilder
113 integer typseh, typseb
117 integer ulsort, langue, codret
119 double precision seuibe, seuihe
120 double precision seuilb, seuilh
121 double precision tabind(indtab)
123 c 0.4. ==> variables locales
130 double precision vmin, vmax
131 double precision vmoy, sigma
132 double precision daux
135 cgn character*8 saux08
136 cgn character*80 repere
138 logical lgaux1, lgaux2, lgaux3
141 parameter (nbmess = 16 )
142 character*80 texte(nblang,nbmess)
143 c ______________________________________________________________________
151 #ifdef _DEBUG_HOMARD_
152 write (ulsort,texte(langue,1)) 'Entree', nompro
158 cgn write (ulsort,90002) 'typenh', typenh
159 cgn write (ulsort,90002) 'pilraf', pilraf
160 cgn write (ulsort,90002) 'pilder', pilder
161 cgn write (ulsort,90002) 'typseh', typseh
162 cgn write (ulsort,90004) 'seuilh', seuilh
163 cgn write (ulsort,90002) 'typseb', typseb
164 cgn write (ulsort,90004) 'seuilb', seuilb
165 cgn 1400 format(5(i5,' : ',i11,' |'))
166 cgn 1401 format(5(i5,' : ',g11.4,' |'))
168 texte(1,4) = '(''Le seuil haut n''''est pas defini.'')'
169 texte(1,5) = '(''Le seuil bas n''''est pas defini.'')'
170 texte(1,6) = '(''Entite '',i10)'
171 texte(1,7) = '(''. Nombre d''''entites actives :'',i10)'
173 >'(''. Nombre d''''entites designees par le support :'',i10)'
174 texte(1,9) = '(5x,a14,'' : seuil haut ='',g13.5,/)'
175 texte(1,10) = '(5x,a14,'' : seuil bas ='',g13.5,/)'
176 texte(1,11) = '(''Recherche des seuils pour les '',a))'
177 texte(1,12) = '(''On prend la valeur brute de l''''indicateur.'')'
179 > '(''On prend la valeur absolue de l''''indicateur.'')'
180 texte(1,14) = '(''Nombre de sommets actuel :'',i10)'
181 texte(1,15) = '(''Nombre de sommets voulu :'',i10)'
182 texte(1,16) = '(''Impossible'')'
184 texte(2,4) = '(''Upper threshold is not defined.'')'
185 texte(2,5) = '(''Lower threshold is not defined.'')'
186 texte(2,6) = '(''Entity '',i10)'
187 texte(2,7) = '(''. Number of active entities :'',i10)'
189 >'(''. Number of entities declared by support of error :'',i10)'
190 texte(2,9) = '(5x,a14,'': Upper threshold ='',g13.5,/)'
191 texte(2,10) = '(5x,a14,'': Lower threshold ='',g13.5,/)'
192 texte(2,11) = '(''Thresholds for the '',a))'
193 texte(2,12) = '(''Inlet value for indicator is taken.'')'
194 texte(2,13) = '(''Absolute value for indicator is taken.'')'
195 texte(2,14) = '(''Number of vertices :'',i10)'
196 texte(2,15) = '(''Targetted number of vertices:'',i10)'
197 texte(2,16) = '(''Impossible'')'
204 #ifdef _DEBUG_HOMARD_
205 write (ulsort,90002) 'typseh', typseh
206 write (ulsort,90002) 'typseb', typseb
209 if ( pilraf.gt.0 .and. typseh.eq.0 .and. nbsoci.le.0 ) then
210 write (ulsort,texte(langue,4))
214 if ( nbiter.gt.0 ) then
216 if ( pilder.gt.0 .and. typseb.eq.0 ) then
217 write (ulsort,texte(langue,5))
222 c 2.2. ==> Par defaut, on prend des valeurs extremes inhibant toute
228 c 2.3. ==> Pour une cible, on va estimer un pourcentage de mailles
229 c 2.3.1. ==> Au premier passage, on va estimer un pourcentage de mailles
231 if ( typseh.eq.5 ) then
233 #ifdef _DEBUG_HOMARD_
234 write (ulsort,texte(langue,14)) nbnop1
235 write (ulsort,texte(langue,15)) nbsoci
238 if ( nbsoci.lt.nbnop1 ) then
239 write (ulsort,texte(langue,14)) nbnop1
240 write (ulsort,texte(langue,15)) nbsoci
241 write (ulsort,texte(langue,16))
245 if ( codret.eq.0 ) then
247 daux = dble(nbsoci)/dble(nbnop1)
248 cgn write (ulsort,90004) 'nbsoci/nbnop1', daux
251 if ( mdim.eq.1 ) then
253 elseif ( mdim.eq.2 ) then
258 cgn write (ulsort,90004) 'daux', daux
259 seuilh = 100.d0 * daux
260 seuilh = min(seuilh, 100.d0)
261 #ifdef _DEBUG_HOMARD_
262 write (ulsort,90004) 'seuilh', seuilh
269 c 2.3.2. ==> Ensuite, on transfere
271 if ( codret.eq.0 ) then
273 if ( ( typseh.eq.0 ) .and. (nbsoci.gt.0 ) ) then
281 #ifdef _DEBUG_HOMARD_
282 write (ulsort,texte(langue,11)) mess14(langue,3,typenh)
286 c 3. si les seuils sont definis par des valeurs absolues
289 #ifdef _DEBUG_HOMARD_
290 write (ulsort,90002) '3. seuils en valeur absolue ; codret',codret
293 if ( codret.eq.0 ) then
295 if ( pilraf.gt.0 .and. typseh.eq.1 ) then
299 if ( pilder.gt.0 .and. typseb.eq.1 ) then
306 c 4. determination des seuils si :
307 c . un des seuils est fourni en relatif
308 c . un des seuils est fourni en pourcentage d'entites
309 c . un des seuils est fourni en mu+n.sigma
310 c . un nombre de noeuds cibles est recherche
313 #ifdef _DEBUG_HOMARD_
314 write (ulsort,90002) '4. determination des seuils ; codret',codret
317 if ( ( pilraf.gt.0 .and. typseh.ge.2 .and. typseh.le.5 ) .or.
318 > ( pilder.gt.0 .and. typseb.ge.2 .and. typseb.le.4 ) ) then
320 c 4.1. ==> allocation du tableau de travail pour uttris
322 if ( codret.eq.0 ) then
324 #ifdef _DEBUG_HOMARD_
325 write (ulsort,90015) 'Allocation pour',
326 > indtab, ' '//mess14(langue,3,typenh)
329 call gmalot ( ntrav1, 'entier ', indtab, ptrav1, codre0 )
331 codret = max ( abs(codre0), codret )
336 c On a besoin de la valeur max dans les cas suivants :
337 c - raffinement ou deraffinement libre, seuil exprime
338 c en relatif et valant plus de 0%
339 c - raffinement libre, seuil exprime en pourcentage
340 c d'elements et valant moins de 0%
341 c - deraffinement libre, seuil exprime en pourcentage
342 c d'elements et valant plus de 100%
343 c On a besoin de la valeur min dans les cas suivants :
344 c - raffinement ou deraffinement libre, seuil exprime
345 c en relatif et valant moins de 100%
346 c - raffinement libre, seuil exprime en pourcentage
347 c d'elements et valant plus de 100%
348 c - deraffinement libre, seuil exprime en pourcentage
349 c d'elements et valant moins de 0%
350 c On a besoin de la valeur moy et de l'ecart-type dans les
352 c - raffinement ou deraffinement libre, seuil exprime
353 c en moyenne + coeff*(ecart-type)
355 c lgaux1 = calcul de la valeur minimale
356 c lgaux2 = calcul de la valeur maximale
357 c lgaux3 = calcul de la valeur moyenne et de l'ecart-type
359 if ( codret.eq.0 ) then
365 c 4.2.1. ==> examen du raffinement
367 if ( pilraf.gt.0 ) then
369 #ifdef _DEBUG_HOMARD_
370 write (ulsort,93030) '4.2.1. Examen du raffinement'
373 if ( typseh.eq.2 ) then
374 if ( abs(seuilh).le.epsima ) then
376 elseif ( abs(seuilh-100.d0).le.epsima ) then
382 c pourcentage d'entites
383 elseif ( typseh.eq.3 .or. typseh.eq.5 ) then
384 if ( abs(seuilh).le.epsima ) then
386 elseif ( abs(seuilh-100.d0).le.epsima ) then
392 elseif ( typseh.eq.4 ) then
398 c 4.2.2. ==> examen du deraffinement
400 if ( pilder.gt.0 ) then
402 #ifdef _DEBUG_HOMARD_
403 write (ulsort,93030) '4.2.2. Examen du deraffinement'
406 if ( typseb.eq.2 ) then
407 if ( abs(seuilb).le.epsima ) then
409 elseif ( abs(seuilb-100.d0).le.epsima ) then
415 c pourcentage d'entites
416 elseif ( typseb.eq.3 ) then
417 if ( abs(seuilb).le.epsima ) then
419 elseif ( abs(seuilb-100.d0).le.epsima ) then
425 elseif ( typseb.eq.4 ) then
430 #ifdef _DEBUG_HOMARD_
431 write (ulsort,99001) 'lgaux1', lgaux1
432 write (ulsort,99001) 'lgaux2', lgaux2
433 write (ulsort,99001) 'lgaux3', lgaux3
440 if ( codret.eq.0 ) then
442 c 4.3.1. ==> Mini/maxi
444 if ( lgaux1 .or. lgaux2 ) then
446 #ifdef _DEBUG_HOMARD_
447 write (ulsort,93030) '4.3.1. Mini/maxi'
452 if ( lgaux1 .and. lgaux2 ) then
453 do 4311 , iaux = 1, indtab
454 vmin = min(vmin,tabind(iaux))
455 vmax = max(vmax,tabind(iaux))
457 elseif ( lgaux2 ) then
458 do 4312 , iaux = 1, indtab
459 vmax = max(vmax,tabind(iaux))
462 do 4313 , iaux = 1, indtab
463 vmin = min(vmin,tabind(iaux))
466 #ifdef _DEBUG_HOMARD_
467 write (ulsort,90004) 'vmin', vmin
468 write (ulsort,90004) 'vmax', vmax
471 c 4.3.2. ==> Moyenne et ecart-type
473 elseif ( lgaux3 ) then
475 #ifdef _DEBUG_HOMARD_
476 write (ulsort,93030) '4.3.2. Moyenne et ecart-type'
481 do 432 , iaux = 1, indtab
482 vmoy = vmoy + tabind(iaux)
483 daux = daux + tabind(iaux)**2
485 vmoy = vmoy/dble(indtab)
486 daux = daux/dble(indtab)
487 sigma = sqrt(daux - vmoy**2)
488 #ifdef _DEBUG_HOMARD_
489 write(ulsort,90004) 'vmoy ', vmoy
490 write(ulsort,90004) 'sigma', sigma
497 c 4.4. ==> Deduction des seuils si exprime en pourcentage d'entites
498 c 4.4.1. ==> si le seuil haut est exprime en pourcentage d'entites,
499 c strictement compris entre 0 et 100, on repere la valeur
502 if ( codret.eq.0 ) then
504 if ( pilraf.gt.0 .and.
505 > ( typseh.eq.3 .or. typseh.eq.5 ) ) then
507 if ( abs(seuilh).gt.epsima .and.
508 > abs(seuilh-100.d0).gt.epsima ) then
510 cgn write (ulsort,1401)(iaux,tabind(iaux),iaux=1,indtab)
512 #ifdef _DEBUG_HOMARD_
513 write (ulsort,texte(langue,3)) 'UTTRIS', nompro
515 call uttris ( seuihe,
516 > iaux, imem(ptrav1),
517 > seuilh, indtab, tabind,
518 > ulsort, langue, codret )
520 cgn codre0 = nint(seuilh*dble(indtab)/100.d0)
521 cgn if ( indtab.ge.2 ) codre0 = min(codre0+1,indtab)
522 cgn write (ulsort,*) '================== ptrav1 ========='
523 cgn write (ulsort,1400) (iaux,imem(ptrav1+iaux-1),iaux=1,codre0)
524 cgn write (ulsort,*) '=================='
525 cgn write (ulsort,*) '=========== ptrav1 trie pour haut ========'
526 cgn write (ulsort,1401)
527 cgn >(iaux,tabind(imem(ptrav1+iaux-1)),iaux=1,codre0)
528 cgn write (ulsort,90004) '==> seuihe',seuihe
535 c 4.4.2. ==> si le seuil bas est exprime en pourcentage d'entites,
536 c strictement compris entre 0 et 100, on repere la valeur
539 if ( codret.eq.0 ) then
541 if ( pilder.gt.0 .and. typseh.eq.3 ) then
543 if ( abs(seuilh).gt.epsima .and.
544 > abs(seuilh-100.d0).gt.epsima ) then
547 #ifdef _DEBUG_HOMARD_
548 write (ulsort,texte(langue,3)) 'UTTRIS', nompro
550 call uttris ( seuibe,
551 > iaux, imem(ptrav1),
552 > seuilb, indtab , tabind,
553 > ulsort, langue, codret )
555 cgn codre0 = nint(seuilb*dble(indtab)/100.d0)
556 cgn if ( indtab.ge.2 ) codre0 = min(codre0+1,indtab)
557 cgn write (ulsort,*) '================== ptrav1 ========='
558 cgn write (ulsort,1400) (iaux,imem(ptrav1+iaux-1),iaux=1,codre0)
559 cgn write (ulsort,*) '=================='
560 cgn write (ulsort,*) '=========== ptrav1 trie pour bas ========'
561 cgn write (ulsort,1401)
562 cgn >(iaux,tabind(imem(ptrav1+iaux-1)),iaux=1,codre0)
563 cgn write (ulsort,90004) '==> seuibe',seuibe
570 c 4.5. ==> Les seuils definitifs
572 if ( codret.eq.0 ) then
574 c 4.5.1. ==> en raffinement
576 if ( pilraf.gt.0 ) then
578 #ifdef _DEBUG_HOMARD_
579 write (ulsort,93030) '4.5.1. Seuil en raffinement'
583 if ( typseh.eq.2 ) then
584 if ( abs(vmax-vmin).le.epsima ) then
585 seuihe = vmax + epsima
586 elseif ( abs(seuilh).le.epsima ) then
587 seuihe = 0.999d0*vmin
588 elseif ( abs(seuilh-100.d0).le.epsima ) then
591 seuihe = vmin + seuilh*(vmax-vmin)/100.d0
593 c pourcentage d'entites
594 elseif ( typseh.eq.3 .or. typseh.eq.5 ) then
595 if ( abs(vmax-vmin).le.epsima ) then
596 seuihe = vmax + epsima
597 elseif ( abs(seuilh).le.epsima ) then
599 elseif ( abs(seuilh-100.d0).le.epsima ) then
600 seuihe = 0.999d0*vmin
602 c moyenne et ecart-type
603 elseif ( typseh.eq.4 ) then
604 if ( abs(sigma).le.epsima ) then
605 seuihe = vmoy + epsima
607 seuihe = vmoy + seuilh*sigma
613 c 4.5.2. ==> en deraffinement
615 if ( pilder.gt.0 ) then
617 #ifdef _DEBUG_HOMARD_
618 write (ulsort,93030) '4.5.2. Seuil en deraffinement'
622 if ( typseb.eq.2 ) then
623 if ( abs(vmax-vmin).le.epsima ) then
624 seuibe = vmin - epsima
625 elseif ( abs(seuilb).le.epsima ) then
626 seuibe = 0.999d0*vmin
627 elseif ( abs(seuilb-100.d0).le.epsima ) then
630 seuibe = vmin + seuilb*(vmax-vmin)/100.d0
632 c pourcentage d'entites
633 elseif ( typseb.eq.3 ) then
634 if ( abs(vmax-vmin).le.epsima ) then
635 seuibe = vmin - epsima
636 elseif ( abs(seuilb).le.epsima ) then
637 seuibe = 0.999d0*vmin
638 elseif ( abs(seuilb-100.d0).le.epsima ) then
641 c moyenne et ecart-type
642 elseif ( typseh.eq.4 ) then
643 if ( abs(sigma).le.epsima ) then
644 seuibe = vmoy - epsima
646 seuibe = vmoy - seuilb*sigma
654 c 4.6. ==> liberation des tableaux temporaires
656 if ( codret.eq.0 ) then
658 call gmlboj ( ntrav1 , codre0 )
660 codret = max ( abs(codre0), codret )
667 c 5. Ecriture sur la sortie standard et sur le fichier recapitulatif
670 #ifdef _DEBUG_HOMARD_
671 write (ulsort,90002) '5. Ecriture standard ; codret',codret
674 if ( pilraf.gt.0 .and.
675 > ( ( typseh.ge.1 .and. typseh.le.5 ) .or.
676 > ( typseh.eq.0 .and. nbsoci.gt.0 ) ) ) then
678 if ( codret.eq.0 ) then
680 write (ulsort,texte(langue,9)) mess14(langue,4,typenh), seuihe
685 cgn repere(1:jaux) = 'Seuil haut '//mess14(langue,4,typenh)
686 cgn#ifdef _DEBUG_HOMARD_
687 cgn write (ulsort,texte(langue,3)) 'UTSYNT', nompro
689 cgn call utsynt ( repere, jaux,
690 cgn > iaux, jaux, seuihe, saux08, jaux,
691 cgn > ulsort, langue, codret )
697 if ( nbiter.gt.0 ) then
699 if ( pilder.gt.0 .and. typseb.ge.1 .and. typseb.le.4 ) then
701 if ( codret.eq.0 ) then
703 write (ulsort,texte(langue,10)) mess14(langue,4,typenh), seuibe
708 cgn repere(1:jaux) = 'Seuil bas '//mess14(langue,4,typenh)
709 cgn#ifdef _DEBUG_HOMARD_
710 cgn write (ulsort,texte(langue,3)) 'UTSYNT', nompro
712 cgn call utsynt ( repere, jaux,
713 cgn > iaux, jaux, seuibe, saux08, jaux,
714 cgn > ulsort, langue, codret )
721 if ( typseh.eq.5 ) then
725 #ifdef _DEBUG_HOMARD_
726 write (ulsort,90004) '==> seuihe', seuihe
727 write (ulsort,90004) '==> seuibe', seuibe
734 if ( codret.ne.0 ) then
738 write (ulsort,texte(langue,1)) 'Sortie', nompro
739 write (ulsort,texte(langue,2)) codret
743 #ifdef _DEBUG_HOMARD_
744 write (ulsort,texte(langue,1)) 'Sortie', nompro