src/tool/Decision/deisfa.F

   1       subroutine deisfa ( ncmpin, usacmp,
   2      >                    trsupp, trindi, qusupp, quindi,
   3      >                    hetare, filare, merare,
   4      >                    posifa, facare,
   5      >                    hettri, aretri, hetqua, arequa,
   6      >                    tabent, tabree,
   7      >                    ulsort, langue, codret)
   8 c ______________________________________________________________________
   9 c
  10 c                             H O M A R D
  11 c
  12 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
  13 c
  14 c Version originale enregistree le 18 juin 1996 sous le numero 96036
  15 c aupres des huissiers de justice Simart et Lavoir a Clamart
  16 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
  17 c aupres des huissiers de justice
  18 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
  19 c
  20 c    HOMARD est une marque deposee d'Electricite de France
  21 c
  22 c Copyright EDF 1996
  23 c Copyright EDF 1998
  24 c Copyright EDF 2002
  25 c Copyright EDF 2020
  26 c ______________________________________________________________________
  27 c
  28 c    traitement des DEcisions - Initialisations - par Saut - FAces
  29 c                   --          -                     -      --
  30 c ______________________________________________________________________
  31 c .        .     .        .                                            .
  32 c .  nom   . e/s . taille .           description                      .
  33 c .____________________________________________________________________.
  34 c . ncmpin .  e  .   1    . nombre de composantes de l'indicateur      .
  35 c . usacmp . e   .   1    . usage des composantes de l'indicateur      .
  36 c .        .     .        . 0 : norme L2                               .
  37 c .        .     .        . 1 : norme infinie -max des valeurs absolues.
  38 c .        .     .        . 2 : valeur relative si une seule composante.
  39 c . trsupp . e   . nbtrto . support pour les triangles                 .
  40 c . trindi . es  . nbtrto . valeurs reelles pour les triangles         .
  41 c . qusupp . e   . nbquto . support pour les quadrangles               .
  42 c . quindi . es  . nbquto . valeurs reelles pour les quadrangles       .
  43 c . hetare . e   . nbarto . historique de l'etat des aretes            .
  44 c . filare . e   . nbarto . fils des aretes                            .
  45 c . merare . e   . nbarto . pere des aretes                            .
  46 c . posifa . e   . nbarto . pointeur sur tableau facare                .
  47 c . facare . e   . nbfaar . liste des faces contenant une arete        .
  48 c . hettri . e   . nbtrto . historique de l'etat des triangles         .
  49 c . aretri . e   .nbtrto*3. numeros des 3 aretes des triangles         .
  50 c . hetqua . e   . nbquto . historique de l'etat des quadrangles       .
  51 c . arequa . e   .nbquto*4. numeros des 4 aretes des quadrangles       .
  52 c . tabent . aux .   *    . tableau auxiliaire entier                  .
  53 c . tabree . aux .   *    . tableau auxiliaire reel                    .
  54 c . ulsort . e   .   1    . numero d'unite logique de la liste standard.
  55 c . langue . e   .    1   . langue des messages                        .
  56 c .        .     .        . 1 : francais, 2 : anglais                  .
  57 c . codret . es  .    1   . code de retour des modules                 .
  58 c .        .     .        . 0 : pas de probleme                        .
  59 c .        .     .        . 2 : probleme dans le traitement            .
  60 c ______________________________________________________________________
  61 c
  62 c====
  63 c 0. declarations et dimensionnement
  64 c====
  65 c
  66 c 0.1. ==> generalites
  67 c
  68       implicit none
  69       save
  70 c
  71       character*6 nompro
  72       parameter ( nompro = 'DEISFA' )
  73 c
  74 #include "nblang.h"
  75 c
  76 c 0.2. ==> communs
  77 c
  78 #include "envex1.h"
  79 c
  80 #include "nombar.h"
  81 #include "nombtr.h"
  82 #include "nombqu.h"
  83 #include "infini.h"
  84 #include "impr02.h"
  85 c
  86 c 0.3. ==> arguments
  87 c
  88       integer ncmpin
  89       integer usacmp
  90       integer trsupp(nbtrto), qusupp(nbquto)
  91       integer hetare(nbarto), filare(nbarto), merare(nbarto)
  92       integer posifa(0:nbarto), facare(nbfaar)
  93       integer hettri(nbtrto), aretri(nbtrto,3)
  94       integer hetqua(nbquto), arequa(nbquto,4)
  95       integer tabent(*)
  96 c
  97       integer ulsort, langue, codret
  98 c
  99       double precision trindi(nbtrto,ncmpin), quindi(nbquto,ncmpin)
 100       double precision tabree(-nbquto:nbtrto,ncmpin)
 101 c
 102 c 0.4. ==> variables locales
 103 c
 104       integer iaux, jaux, kaux
 105       integer jdeb, jfin
 106       integer lgpile, nupile
 107       integer laface, larefa, larete
 108       integer merear
 109       integer nbarfa
 110       integer nrcomp
 111 cgn      integer glop
 112 c
 113       double precision daux1, daux2
 114       integer lgdaux
 115       parameter( lgdaux = 100 )
 116       double precision daux(lgdaux), vect(lgdaux)
 117 c
 118       logical afaire
 119       logical calcul
 120 c
 121       integer nbmess
 122       parameter (nbmess = 11 )
 123       character*80 texte(nblang,nbmess)
 124 c ______________________________________________________________________
 125 c
 126 c====
 127 c 1. initialisation
 128 c====
 129 c
 130 c 1.1. ==> Les messages
 131 c
 132 #include "impr01.h"
 133 c
 134 #ifdef _DEBUG_HOMARD_
 135       write (ulsort,texte(langue,1)) 'Entree', nompro
 136       call dmflsh (iaux)
 137 #endif
 138 c
 139       texte(1,4) = '(''. Saut a la traversee des '',a)'
 140       texte(1,5) =
 141      > '(''On veut'',i6,'' composantes, mais taille de daux ='',i6)'
 142       texte(1,9) = '(''. Norme L2 des composantes.'')'
 143       texte(1,10) = '(''. Norme infinie des composantes.'')'
 144       texte(1,11) = '(''. Valeur relative de la composante.'')'
 145 c
 146       texte(2,4) = '(''. Jump through the '',a)'
 147       texte(2,5) =
 148      > '(i6,''components are requested, but size of daux equals'',i6)'
 149       texte(2,9) = '(''. L2 norm of components.'')'
 150       texte(2,10) = '(''. Infinite norm of components.'')'
 151       texte(2,11) = '(''. Relative value for the component.'')'
 152 cgn      print *, hettri
 153 cgn      print *,  aretri
 154 cgn      print *, hetqua
 155 cgn      print *, arequa
 156 c
 157 #include "impr03.h"
 158 c
 159       codret = 0
 160 c
 161 c 1.2. ==> controle
 162 c
 163       if ( ncmpin.gt.lgdaux ) then
 164         write (ulsort,texte(langue,5)) ncmpin, lgdaux
 165         codret = 1
 166       endif
 167 c
 168 #ifdef _DEBUG_HOMARD_
 169       write (ulsort,texte(langue,9+usacmp))
 170 #endif
 171 c
 172 #ifdef _DEBUG_HOMARD_
 173       write (ulsort,texte(langue,4)) mess14(langue,3,1)
 174 #endif
 175 c
 176 c====
 177 c 2. sauvegarde de l'indicateur
 178 c====
 179 c
 180       if ( codret.eq.0 ) then
 181 c
 182       do 211 , iaux = 1 , nbtrto
 183         if ( trsupp(iaux).ne.0 ) then
 184           do 2110 , nrcomp = 1 , ncmpin
 185             tabree(iaux,nrcomp) = trindi(iaux,nrcomp)
 186  2110     continue
 187         endif
 188   211 continue
 189 c
 190       do 212 , iaux = 1 , nbquto
 191         if ( qusupp(iaux).ne.0 ) then
 192           do 2120 , nrcomp = 1 , ncmpin
 193             tabree(-iaux,nrcomp) = quindi(iaux,nrcomp)
 194  2120     continue
 195         endif
 196   212 continue
 197 c
 198 c====
 199 c 3. traitement des indicateurs portant sur les faces
 200 c====
 201 c
 202       do 3 , laface = -nbquto , nbtrto
 203 c
 204 cgn        glop=0
 205 c
 206 c 3.0. ==> On y va ?
 207 c
 208         afaire = .false.
 209         if ( laface.lt.0 ) then
 210           if ( qusupp(-laface).ne.0 ) then
 211             afaire = .true.
 212             nbarfa = 4
 213           endif
 214         elseif ( laface.gt.0 ) then
 215           if ( trsupp(laface).ne.0 ) then
 216             afaire = .true.
 217             nbarfa = 3
 218           endif
 219         endif
 220 c
 221         if ( afaire ) then
 222 c
 223 cgn        if ( laface.ge.-42 .and. laface.le.49 ) then
 224 cgn        glop = 1
 225 cgn        endif
 226 cgn        if ( glop.eq.1) then
 227 cgn        print *,'==========================='
 228 cgn        print *,'FACE = ',laface,', d''indic ',
 229 cgn     > (tabree(laface,nrcomp),nrcomp=1,ncmpin)
 230 cgn        endif
 231 c
 232           daux1 = vinfne
 233 c
 234           do 31 , iaux = 1 , nbarfa
 235 c
 236 c 3.1. ==> pour une des aretes de la face, on stocke les numeros
 237 c            des faces voisines, en descendant les parentes.
 238 c            ensuite, on stocke la premiere arete mere dont une des
 239 c            faces voisines est active
 240 c
 241 c
 242             if ( laface.gt.0 ) then
 243               larefa = aretri(laface,iaux)
 244             else
 245               larefa = arequa(-laface,iaux)
 246             endif
 247 cgn        if ( glop.eq.1) then
 248 cgn          print *,'.', iaux,'-ieme arete de la face : ',larefa
 249 cgn        endif
 250             lgpile = 1
 251             tabent(lgpile) = larefa
 252             nupile = 1
 253 c
 254   310       continue
 255 c
 256             larete = tabent(nupile)
 257             if ( mod(hetare(larete),10).ne.0 ) then
 258 cgn        if ( glop.eq.1) then
 259 cgn          print *,'..  des filles'
 260 cgn        endif
 261               lgpile = lgpile + 1
 262               tabent(lgpile) = filare(larete)
 263               lgpile = lgpile + 1
 264               tabent(lgpile) = filare(larete) + 1
 265 cgn        if ( glop.eq.1) then
 266 cgn          print *,'.... ajout de ',tabent(lgpile), ' a la pile'
 267 cgn        endif
 268             endif
 269 c
 270             nupile = nupile + 1
 271             if ( nupile.le.lgpile ) then
 272               goto 310
 273             endif
 274 c
 275 c 3.2. ==> pour chaque arete de la pile : si elle est active, on
 276 c            cherche le max de l'ecart entre la valeur de l'indicateur
 277 c            sur la face voisine et celle sur la face courante
 278 c
 279             do 32 , nupile = 1 , lgpile
 280               larete = tabent(nupile)
 281               if ( mod(hetare(larete),10).eq.0 ) then
 282 cgn        if ( glop.eq.1) then
 283 cgn       print *,'...... Examen de la pile, arete : ',larete
 284 cgn        endif
 285                 jdeb = posifa(larete-1)+1
 286                 jfin = posifa(larete)
 287                 do 321 , jaux = jdeb, jfin
 288                   kaux = facare(jaux)
 289                   if ( kaux.ne.laface ) then
 290 cgn        if ( glop.eq.1) then
 291 cgn          print *,'........ ==> face ', kaux,' : ',
 292 cgn     > (tabree(kaux,nrcomp),nrcomp=1,ncmpin)
 293 cgn        endif
 294                     calcul = .false.
 295                     if ( kaux.gt.0 ) then
 296                       if ( trsupp(kaux).ne.0 ) then
 297                         calcul = .true.
 298                         do 3211 , nrcomp = 1 , ncmpin
 299                           daux(nrcomp) = tabree(kaux,nrcomp)
 300      >                                 - tabree(laface,nrcomp)
 301  3211                   continue
 302                       endif
 303                     elseif ( kaux.lt.0 ) then
 304                       if ( qusupp(-kaux).ne.0 ) then
 305                         calcul = .true.
 306                         do 3212 , nrcomp = 1 , ncmpin
 307                           daux(nrcomp) = tabree(kaux,nrcomp)
 308      >                                 - tabree(laface,nrcomp)
 309  3212                   continue
 310                       endif
 311                     endif
 312                     if ( calcul ) then
 313 c            calcul de la norme de l'ecart
 314 c            si on a passe le max, on stocke
 315                       if ( usacmp.eq.0 ) then
 316                         daux2 = daux(1)**2
 317                         do 32111 , nrcomp = 2 , ncmpin
 318                           daux2 = daux2 + daux(nrcomp)**2
 319 32111                   continue
 320                       elseif ( usacmp.eq.1 ) then
 321                         daux2 = abs(daux(1))
 322                         do 32112 , nrcomp = 2 , ncmpin
 323                           daux2 = max(daux2,abs(daux(nrcomp)))
 324 32112                   continue
 325                       else
 326                         daux2 = daux(1)
 327                       endif
 328                       if ( daux2.gt.daux1 ) then
 329                         daux1 = daux2
 330                         do 3213 , nrcomp = 1 , ncmpin
 331                           vect(nrcomp) = daux(nrcomp)
 332  3213                   continue
 333                       endif
 334                     endif
 335 cgn        if ( glop.eq.1) then
 336 cgn          print *,'........ daux1 ', daux1
 337 cgn        endif
 338                   endif
 339   321           continue
 340 cgn         if ( glop.eq.1) then
 341 cgn          print *,'...... ==> valeur finale =', daux1
 342 cgn        endif
 343              endif
 344    32      continue
 345 c
 346 c 3.3. ==> on remonte la parente pour pieger les non-conformites
 347 c
 348             larete = larefa
 349 c
 350    33       continue
 351 c
 352             merear = merare(larete)
 353             if ( merear.gt.0 ) then
 354 cgn        if ( glop.eq.1) then
 355 cgn          print *,'......', larete,' a une mere : ',merear
 356 cgn        endif
 357               jdeb = posifa(merear-1)+1
 358               jfin = posifa(merear)
 359               if ( jdeb.gt.jfin ) then
 360                 larete = merear
 361                 goto 33
 362               else
 363                 do 331 , jaux = jdeb, jfin
 364                   kaux = facare(jaux)
 365                   if ( kaux.ne.laface ) then
 366 cgn        if ( glop.eq.1) then
 367 cgn          print *,'.......... ==> face ', kaux,' : ',
 368 cgn     > (tabree(kaux,nrcomp),nrcomp=1,ncmpin)
 369 cgn        endif
 370                     calcul = .false.
 371                     if ( kaux.gt.0 ) then
 372                       if ( mod(hettri(kaux),10).eq.0 ) then
 373                         if ( trsupp(kaux).ne.0 ) then
 374                           calcul = .true.
 375                           do 3311 , nrcomp = 1 , ncmpin
 376                             daux(nrcomp) = tabree(kaux,nrcomp)
 377      >                                   - tabree(laface,nrcomp)
 378  3311                     continue
 379                         endif
 380                       endif
 381                     elseif ( kaux.lt.0 ) then
 382                       if ( mod(hetqua(-kaux),100).eq.0 ) then
 383                         if ( qusupp(-kaux).ne.0 ) then
 384                           calcul = .true.
 385                           do 3312 , nrcomp = 1 , ncmpin
 386                             daux(nrcomp) = tabree(kaux,nrcomp)
 387      >                                   - tabree(laface,nrcomp)
 388  3312                     continue
 389                         endif
 390                       endif
 391                     endif
 392                     if ( calcul ) then
 393                       if ( usacmp.eq.0 ) then
 394                         daux2 = daux(1)**2
 395                         do 33111 , nrcomp = 2 , ncmpin
 396                           daux2 = daux2 + daux(nrcomp)**2
 397 33111                   continue
 398                       elseif ( usacmp.eq.1 ) then
 399                         daux2 = abs(daux(1))
 400                         do 33112 , nrcomp = 2 , ncmpin
 401                           daux2 = max(daux2,abs(daux(nrcomp)))
 402 33112                   continue
 403                       else
 404                         daux2 = daux(1)
 405                       endif
 406                       if ( daux2.gt.daux1 ) then
 407                         daux1 = daux2
 408                         do 3313 , nrcomp = 1 , ncmpin
 409                           vect(nrcomp) = daux(nrcomp)
 410  3313                   continue
 411                       endif
 412                     endif
 413 cgn        if ( glop.eq.1) then
 414 cgn          print *,'.......... daux1 ', daux1
 415 cgn        endif
 416                   endif
 417   331           continue
 418 cgn         if ( glop.eq.1) then
 419 cgn          print *,'........ ==> valeur finale =', daux1
 420 cgn        endif
 421               endif
 422             endif
 423 c
 424    31     continue
 425 c
 426 c 3.4. ==> stockage
 427 c
 428 cgn      if ( glop.eq.1) then
 429 cgn       write(ulsort,20000) 'Final    '//
 430 cgn     > 'face', laface,' : ',
 431 cgn     > mess14(langue,1,typenh), laface,' : ',
 432 cgn     > (vect(nrcomp),nrcomp=1,ncmpin)
 433 cgn      endif
 434           if ( laface.gt.0 ) then
 435             do 341 , nrcomp = 1 , ncmpin
 436               trindi(laface,nrcomp) = vect(nrcomp)
 437   341       continue
 438           else
 439             do 342 , nrcomp = 1 , ncmpin
 440               quindi(-laface,nrcomp) = vect(nrcomp)
 441   342       continue
 442           endif
 443 c
 444         endif
 445 c
 446     3 continue
 447 c
 448       endif
 449 c
 450 c====
 451 c 4. la fin
 452 c====
 453 c
 454       if ( codret.ne.0 ) then
 455 c
 456 #include "envex2.h"
 457 c
 458       write (ulsort,texte(langue,1)) 'Sortie', nompro
 459       write (ulsort,texte(langue,2)) codret
 460       write (ulsort,texte(langue,4)) mess14(langue,3,1)
 461 c
 462       endif
 463 c
 464 #ifdef _DEBUG_HOMARD_
 465       write (ulsort,texte(langue,1)) 'Sortie', nompro
 466       call dmflsh (iaux)
 467 #endif
 468 c
 469       end