src/tool/Decision/deisv1.F

   1       subroutine deisv1 ( ncmpin, usacmp, nbvoto, typenh,
   2      >                    vosupp, voindi,
   3      >                    hettri, filtri, pertri,
   4      >                    hetqua, filqua, perqua,
   5      >                    hetvol, facvol,
   6      >                    voltri, volqua,
   7      >                    tabent, voinin,
   8      >                    ulsort, langue, codret )
   9 c ______________________________________________________________________
  10 c
  11 c                             H O M A R D
  12 c
  13 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
  14 c
  15 c Version originale enregistree le 18 juin 1996 sous le numero 96036
  16 c aupres des huissiers de justice Simart et Lavoir a Clamart
  17 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
  18 c aupres des huissiers de justice
  19 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
  20 c
  21 c    HOMARD est une marque deposee d'Electricite de France
  22 c
  23 c Copyright EDF 1996
  24 c Copyright EDF 1998
  25 c Copyright EDF 2002
  26 c Copyright EDF 2020
  27 c ______________________________________________________________________
  28 c
  29 c    traitement des DEcisions - Initialisations - par Saut - Volumes - 1
  30 c                   --          -                     -      -         -
  31 c
  32 c    On traite ici uniquement les sauts entre volumes pour les cas ou
  33 c    il n'y a que des tetraedres ou que des hexaedres concernes par
  34 c    un indicateur d'erreur.
  35 c    Quand l'indicateur est reparti sur des entites de nature
  36 c    differente, on gerera le saut dans deisv2.
  37 c ______________________________________________________________________
  38 c .        .     .        .                                            .
  39 c .  nom   . e/s . taille .           description                      .
  40 c .____________________________________________________________________.
  41 c . ncmpin .  e  .   1    . nombre de composantes de l'indicateur      .
  42 c . usacmp . e   .   1    . usage des composantes de l'indicateur      .
  43 c .        .     .        . 0 : norme L2                               .
  44 c .        .     .        . 1 : norme infinie -max des valeurs absolues.
  45 c .        .     .        . 2 : valeur relative si une seule composante.
  46 c . nbvoto . e   .   1    . nombre de volumes total                    .
  47 c . typenh . e   .   1    . 3 : tetraedres                             .
  48 c .        .     .        . 6 : hexaedres                              .
  49 c . vosupp . e   . nbvoto . support pour les volumes                   .
  50 c . voindi .  s  . nbvoto . valeurs reelles pour les volumes           .
  51 c . hettri . e   . nbtrto . historique de l'etat des triangles         .
  52 c . filtri . e   . nbtrto . premier fils des triangles                 .
  53 c . pertri . e   . nbtrto . pere des triangles                         .
  54 c . hetqua . e   . nbquto . historique de l'etat des quadrangles       .
  55 c . filqua . e   . nbquto . fils des quadrangles                       .
  56 c . perqua . e   . nbquto . pere des quadrangles                       .
  57 c . hetvol . e   . nbvoto . historique de l'etat des volumes           .
  58 c . facvol . e   .nbvoto*n. numeros des n faces des volumes            .
  59 c . voltri . es  .2*nbtrto. numeros des 2 volumes par triangle         .
  60 c .        .     .        . voltri(i,k) definit le i-eme voisin de k   .
  61 c .        .     .        .   0 : pas de voisin                        .
  62 c .        .     .        . j>0 : tetraedre j                          .
  63 c .        .     .        . j<0 : pyramide/pentaedre dans pypetr(1/2,j).
  64 c . volqua . e   .2*nbquto. numeros des 2 volumes par quadrangle       .
  65 c .        .     .        . volqua(i,k) definit le i-eme voisin de k   .
  66 c .        .     .        .   0 : pas de voisin                        .
  67 c .        .     .        . j>0 : hexaedre j                           .
  68 c .        .     .        . j<0 : pyramide/pentaedre dans pypequ(1/2,j).
  69 c . tabent . aux .   *    . tableau auxiliaire entier                  .
  70 c . voinin . aux .   *    . sauvegarde des valeurs des indicateurs     .
  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 : mauvais typenh                         .
  77 c ______________________________________________________________________
  78 c
  79 c====
  80 c 0. declarations et dimensionnement
  81 c====
  82 c
  83 c 0.1. ==> generalites
  84 c
  85       implicit none
  86       save
  87 c
  88       character*6 nompro
  89       parameter ( nompro = 'DEISV1' )
  90 c
  91 #include "nblang.h"
  92 c
  93 c 0.2. ==> communs
  94 c
  95 #include "envex1.h"
  96 c
  97 #include "infini.h"
  98 #include "impr02.h"
  99 #include "nombtr.h"
 100 #include "nombqu.h"
 101 c
 102 c 0.3. ==> arguments
 103 c
 104       integer ncmpin
 105       integer usacmp
 106       integer nbvoto, typenh
 107 c
 108       integer vosupp(nbvoto)
 109       integer hettri(nbtrto), filtri(nbtrto), pertri(nbtrto)
 110       integer hetqua(nbquto), filqua(nbquto), perqua(nbquto)
 111       integer voltri(2,nbtrto)
 112       integer volqua(2,nbquto)
 113       integer hetvol(nbvoto), facvol(nbvoto,*)
 114       integer tabent(2,*)
 115 c
 116       integer ulsort, langue, codret
 117 c
 118       double precision voindi(nbvoto,ncmpin)
 119       double precision voinin(nbvoto,ncmpin)
 120 c
 121 c 0.4. ==> variables locales
 122 c
 123       integer iaux, jaux, kaux
 124       integer lgpile, nupile
 125       integer levolu
 126       integer lafavo, laface, tyface, typfac
 127       integer etat
 128       integer merefa
 129       integer nbfavo
 130       integer nrcomp
 131 cgn      integer glop
 132 c
 133       double precision daux1, daux2
 134       integer lgdaux
 135       parameter( lgdaux = 100 )
 136       double precision daux(lgdaux), vect(lgdaux)
 137 c
 138       logical calcul
 139 c
 140       integer nbmess
 141       parameter (nbmess = 11 )
 142       character*80 texte(nblang,nbmess)
 143 c ______________________________________________________________________
 144 c
 145 c====
 146 c 1. initialisation
 147 c====
 148 c
 149 #include "impr01.h"
 150 c
 151 #ifdef _DEBUG_HOMARD_
 152       write (ulsort,texte(langue,1)) 'Entree', nompro
 153       call dmflsh (iaux)
 154 #endif
 155 c
 156       texte(1,4) = '(''. Saut a la traversee des '',a)'
 157       texte(1,5) =
 158      > '(''On veut'',i6,'' composantes, mais taille de daux ='',i6)'
 159       texte(1,6) = ' (''Type d''''entite incorrect :'',i10)'
 160       texte(1,9) = '(''. Norme L2 des composantes.'')'
 161       texte(1,10) = '(''. Norme infinie des composantes.'')'
 162       texte(1,11) = '(''. Valeur relative de la composante.'')'
 163 c
 164       texte(2,4) = '(''. Jump through the '',a)'
 165       texte(2,5) =
 166      > '(i6,''components are requested, but size of daux equals'',i6)'
 167       texte(2,6) = ' (''Uncorrect type of entity: :'',i10)'
 168       texte(2,9) = '(''. L2 norm of components.'')'
 169       texte(2,10) = '(''. Infinite norm of components.'')'
 170       texte(2,11) = '(''. Relative value for the component.'')'
 171 c
 172 #include "impr03.h"
 173 c
 174       codret = 0
 175 c
 176 c 1.2. ==> controle
 177 c
 178       if ( ncmpin.gt.lgdaux ) then
 179         write (ulsort,texte(langue,5)) ncmpin, lgdaux
 180         codret = 1
 181       endif
 182 c
 183 #ifdef _DEBUG_HOMARD_
 184       write (ulsort,texte(langue,9+usacmp))
 185 #endif
 186 c
 187 c 1.3. ==> Les types de faces : triangle (2) ou quadrangle (4)
 188 c
 189       if ( codret.eq.0 ) then
 190 c
 191       if ( typenh.eq.3 ) then
 192         nbfavo = 4
 193         typfac = 2
 194       elseif ( typenh.eq.6 ) then
 195         nbfavo = 6
 196         typfac = 4
 197       else
 198         codret = 1
 199         write (ulsort,texte(langue,6)) typenh
 200       endif
 201 c
 202 #ifdef _DEBUG_HOMARD_
 203       write (ulsort,texte(langue,4)) mess14(langue,3,typfac)
 204 #endif
 205 c
 206       endif
 207 c
 208 c====
 209 c 2. On parcourt tous les volumes.
 210 c    On calcule l'ecart entre la valeur de l'indicateur sur le volume
 211 c    courant et sur les voisins.
 212 c    On garde le max au sens de la norme voulue
 213 c    levolu = numero local dans la categorie
 214 c====
 215 #ifdef _DEBUG_HOMARD_
 216       write (ulsort,90002) '2. parcours des volumes ; codret', codret
 217 #endif
 218 c
 219       do 2 , levolu = 1 , nbvoto
 220 c
 221         if ( codret.eq.0 ) then
 222 c
 223         if ( vosupp(levolu).ne.0 ) then
 224 c
 225 c 2.1. ==> Exploration de toutes les faces du volume
 226 c
 227           daux1 = vinfne
 228 c
 229           do 21 , iaux = 1 , nbfavo
 230 c
 231 c 2.1.1. ==> pour chaque face du volume, on stocke les numeros
 232 c            des faces voisines, en descendant les parentes.
 233 c            Au final, on stocke la premiere face mere active
 234 c
 235             lafavo = facvol(levolu,iaux)
 236             tyface = typfac
 237 cgn            if ( glop.eq.1) then
 238 cgn              write(ulsort,90002) '. face de rang',iaux
 239 cgn            endif
 240 c
 241 #ifdef _DEBUG_HOMARD_
 242         write (ulsort,texte(langue,3)) 'DEISV3', nompro
 243 #endif
 244             call deisv3 ( lafavo, tyface,
 245      >                    hettri, filtri,
 246      >                    hetqua, filqua,
 247      >                    lgpile, tabent,
 248      >                    ulsort, langue, codret )
 249 c
 250 c 2.1.2. ==> pour chaque face de la pile : si elle est active, on
 251 c            cherche le max de l'ecart entre la valeur de l'indicateur
 252 c            sur le volume voisin et celle sur le volume courant
 253 c            on sait que les faces sont toutes de meme type
 254 c
 255             do 212 , nupile = 1 , lgpile
 256 c
 257               laface = tabent(1,nupile)
 258 c
 259 c 2.1.2.1. ==> reperage selon la face
 260 c
 261               if ( typfac.eq.2 ) then
 262                 etat = mod(hettri(laface),10)
 263               else
 264                 etat = mod(hetqua(laface),100)
 265               endif
 266 c
 267 c 2.1.2.2. ==> traitement
 268 c
 269               if ( etat.eq.0 ) then
 270 cgn        if ( glop.eq.1) then
 271 cgn       write(ulsort,90002)'.. Voisinage par le '//
 272 cgn     >            mess14(langue,1,typfac),laface
 273 cgn        endif
 274                 do 2122 , jaux = 1 , 2
 275 c
 276                   calcul = .false.
 277 c
 278                   if ( typfac.eq.2 ) then
 279                     kaux = voltri(jaux,laface)
 280                   else
 281                     kaux = volqua(jaux,laface)
 282                   endif
 283 c
 284                   if ( kaux.gt.0 ) then
 285 c
 286                     if ( vosupp(kaux).ne.0 ) then
 287 c
 288                       if ( kaux.ne.levolu ) then
 289                         calcul = .true.
 290                         do 21221 , nrcomp = 1 , ncmpin
 291                           daux(nrcomp) = voinin(kaux,nrcomp)
 292      >                                 - voinin(levolu,nrcomp)
 293 21221                   continue
 294 cgn        if ( glop.eq.1) then
 295 cgn        write(ulsort,90054)'...... ==> ecart avec  ', kaux,' : ',
 296 cgn     >          (daux(nrcomp),nrcomp=1,ncmpin)
 297 cgn        endif
 298 c
 299                       endif
 300 c
 301                     endif
 302 c
 303                   elseif ( kaux.lt.0 ) then
 304 cgn        if ( glop.eq.1) then
 305 cgn       write(ulsort,*)'.... Autre type de volume voisin.'
 306 cgn        endif
 307                     codret = codret + 1
 308                   endif
 309 c
 310 c                 calcul de la norme de l'ecart
 311 c                 si on a passe le max, on stocke
 312 c
 313                   if ( calcul ) then
 314 c
 315                     if ( usacmp.eq.0 ) then
 316                       daux2 = daux(1)**2
 317                       do 21222 , nrcomp = 2 , ncmpin
 318                         daux2 = daux2 + daux(nrcomp)**2
 319 21222                 continue
 320                     elseif ( usacmp.eq.1 ) then
 321                       daux2 = abs(daux(1))
 322                       do 21223 , nrcomp = 2 , ncmpin
 323                         daux2 = max(daux2,abs(daux(nrcomp)))
 324 21223                 continue
 325                     else
 326                       daux2 = daux(1)
 327                     endif
 328                     if ( daux2.gt.daux1 ) then
 329                       daux1 = daux2
 330                       do 21224 , nrcomp = 1 , ncmpin
 331                         vect(nrcomp) = daux(nrcomp)
 332 21224                 continue
 333                     endif
 334 c
 335                   endif
 336 c
 337  2122           continue
 338 c
 339               endif
 340 c
 341   212       continue
 342 c
 343 c 2.3. ==> on remonte la parente pour pieger les non-conformites
 344 c
 345             laface = lafavo
 346 c
 347   231       continue
 348 c
 349 cgn        if ( glop.eq.1) then
 350 cgn          write(ulsort,90002)'.... Parente du '//
 351 cgn     >            mess14(langue,1,typfac),laface
 352 cgn        endif
 353             if ( typfac.eq.2 ) then
 354               merefa = pertri(laface)
 355             else
 356               merefa = perqua(laface)
 357             endif
 358 c
 359             if ( merefa.gt.0 ) then
 360 cgn        if ( glop.eq.1) then
 361 cgn          write(ulsort,90006)'.... le '//
 362 cgn     >       mess14(langue,1,typfac), laface,' a une mere : ',merefa
 363 cgn        endif
 364 c
 365               if ( typfac.eq.2 ) then
 366                 kaux = voltri(2,merefa)
 367               else
 368                 kaux = volqua(2,merefa)
 369               endif
 370 c
 371               calcul = .false.
 372               if ( kaux.eq.0 ) then
 373                 laface = merefa
 374                 goto 231
 375               elseif ( kaux.gt.0 ) then
 376                 do 232 , jaux = 1 , 2
 377                   if ( typfac.eq.2 ) then
 378                     kaux = voltri(jaux,merefa)
 379                   else
 380                     kaux = volqua(jaux,merefa)
 381                   endif
 382                   if ( kaux.gt.0 ) then
 383                     if ( mod(hetvol(kaux),100).eq.0 ) then
 384                       calcul = .true.
 385                       do 2311 , nrcomp = 1 , ncmpin
 386                         daux(nrcomp) = voinin(kaux,nrcomp)
 387      >                               - voinin(levolu,nrcomp)
 388  2311                 continue
 389 cgn        if ( glop.eq.1) then
 390 cgn        write(ulsort,90054)'...... ==> ecart avec ', kaux,' : ',
 391 cgn     >          (daux(nrcomp),nrcomp=1,ncmpin)
 392 cgn        endif
 393                     endif
 394                   elseif ( kaux.lt.0 ) then
 395                     codret = codret + 1
 396                   endif
 397   232           continue
 398               elseif ( kaux.lt.0 ) then
 399                 codret = codret + 1
 400               endif
 401               if ( calcul ) then
 402                 if ( usacmp.eq.0 ) then
 403                   daux2 = daux(1)**2
 404                   do 23111 , nrcomp = 2 , ncmpin
 405                     daux2 = daux2 + daux(nrcomp)**2
 406 23111             continue
 407                 elseif ( usacmp.eq.1 ) then
 408                   daux2 = abs(daux(1))
 409                   do 23112 , nrcomp = 2 , ncmpin
 410                     daux2 = max(daux2,abs(daux(nrcomp)))
 411 23112             continue
 412                 else
 413                   daux2 = daux(1)
 414                 endif
 415                 if ( daux2.gt.daux1 ) then
 416                   daux1 = daux2
 417                   do 2312 , nrcomp = 1 , ncmpin
 418                     vect(nrcomp) = daux(nrcomp)
 419  2312             continue
 420                 endif
 421               endif
 422 c
 423             endif
 424 c
 425    21     continue
 426 c
 427 c 2.4. ==> stockage
 428 c
 429 cgn        if ( glop.eq.1 ) then
 430 cgn       write(ulsort,90054) 'Final    '//
 431 cgn     > mess14(langue,1,typenh), levolu, ' : ',
 432 cgn     > (vect(nrcomp),nrcomp=1,ncmpin)
 433 cgn        endif
 434 c
 435           do 241 , nrcomp = 1 , ncmpin
 436             voindi(levolu,nrcomp) = vect(nrcomp)
 437   241     continue
 438 c
 439         endif
 440 c
 441         endif
 442 c
 443     2 continue
 444 c
 445 c====
 446 c 3. la fin
 447 c====
 448 c
 449       if ( codret.ne.0 ) then
 450 c
 451 #include "envex2.h"
 452 c
 453       write (ulsort,texte(langue,1)) 'Sortie', nompro
 454       write (ulsort,texte(langue,2)) codret
 455       if ( codret.ne.1 ) then
 456       write (ulsort,texte(langue,4)) mess14(langue,3,typenh)
 457       endif
 458 c
 459       endif
 460 c
 461 #ifdef _DEBUG_HOMARD_
 462       write (ulsort,texte(langue,1)) 'Sortie', nompro
 463       call dmflsh (iaux)
 464 #endif
 465 c
 466       end