src/tool/AP_Conversion/pcspt2.F

   1       subroutine pcspt2 ( etan, etanp1, tehn, tehnp1,
   2      >                    prfcan, prfcap,
   3      >                    hettet, filtet, nbante, anfite,
   4      >                    nteeca, ntesca,
   5      >                    nbfonc, ngauss, vafoen, vafott,
   6      >                    ulsort, langue, codret )
   7 c ______________________________________________________________________
   8 c
   9 c                             H O M A R D
  10 c
  11 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
  12 c
  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
  18 c
  19 c    HOMARD est une marque deposee d'Electricite de France
  20 c
  21 c Copyright EDF 1996
  22 c Copyright EDF 1998
  23 c Copyright EDF 2002
  24 c Copyright EDF 2020
  25 c ______________________________________________________________________
  26 c
  27 c    aPres adaptation - Conversion de Solution Points de Gauss -
  28 c     -                 -             -        -
  29 c                       Tetraedres d'etat anterieur 2
  30 c                       -                           -
  31 c ______________________________________________________________________
  32 c .        .     .        .                                            .
  33 c .  nom   . e/s . taille .           description                      .
  34 c .____________________________________________________________________.
  35 c . etan   . e   .    1   . ETAt du tetraedre a l'iteration N          .
  36 c . etanp1 . e   .    1   . ETAt du tetraedre a l'iteration N+1        .
  37 c . tehn   . e   .    1   . TEtraedre courant en numerotation Homard   .
  38 c .        .     .        . a l'iteration N                            .
  39 c . tehnp1 . e   .    1   . TEtraedre courant en numerotation Homard   .
  40 c .        .     .        . a l'iteration N+1                          .
  41 c . prfcan . e   .   *    . En numero du calcul a l'iteration n :      .
  42 c .        .     .        . 0 : l'entite est absente du profil         .
  43 c .        .     .        . i : l'entite est au rang i dans le profil  .
  44 c . prfcap . es  .   *    . En numero du calcul a l'iteration n+1 :    .
  45 c .        .     .        . 0 : l'entite est absente du profil         .
  46 c .        .     .        . 1 : l'entite est presente dans le profil   .
  47 c . hettet . e   . nbteto . historique de l'etat des tetraedres        .
  48 c . filtet . e   . nbteto . premier fils des tetraedres                .
  49 c . nbante . e   .   1    . nombre de tetraedres decoupes par          .
  50 c .        .     .        . conformite sur le maillage avant adaptation.
  51 c . anfite . e   . nbante . tableau filtet du maillage de l'iteration n.
  52 c . nteeca . e   . reteto . numero des tetraedres dans le calcul entree.
  53 c . ntesca . e   . rsteto . numero des tetraedres dans le calcul sortie.
  54 c . nbfonc . e   .    1   . nombre de fonctions elements de volume     .
  55 c . ngauss . e   .   1    . nbre de points de Gauss des fonctions pg   .
  56 c . vafoen . e   . nbfonc*. variables en entree de l'adaptation        .
  57 c .        .     . ngauss*.                                            .
  58 c .        .     . nbeven .                                            .
  59 c . vafott . es  . nbfonc*. variables en sortie de l'adaptation        .
  60 c .        .     . ngauss*.                                            .
  61 c .        .     . nbevso .                                            .
  62 c . ulsort . e   .   1    . numero d'unite logique de la liste standard.
  63 c . langue . e   .    1   . langue des messages                        .
  64 c .        .     .        . 1 : francais, 2 : anglais                  .
  65 c . codret . es  .    1   . code de retour des modules                 .
  66 c .        .     .        . 0 : pas de probleme                        .
  67 c .        .     .        . 1 : probleme                               .
  68 c ______________________________________________________________________
  69 c
  70 c====
  71 c 0. declarations et dimensionnement
  72 c====
  73 c
  74 c 0.1. ==> generalites
  75 c
  76       implicit none
  77       save
  78 c
  79       character*6 nompro
  80       parameter ( nompro = 'PCSPT2' )
  81 c
  82 #include "nblang.h"
  83 c
  84 c 0.2. ==> communs
  85 c
  86 #include "nombte.h"
  87 #include "nombsr.h"
  88 #include "nomber.h"
  89 c
  90 c 0.3. ==> arguments
  91 c
  92       integer etan, etanp1, tehn, tehnp1
  93       integer nbfonc, ngauss
  94       integer prfcan(*), prfcap(*)
  95       integer hettet(nbteto), filtet(nbteto)
  96       integer nbante
  97       integer anfite(nbante)
  98       integer nteeca(reteto), ntesca(rsteto)
  99 c
 100       double precision vafoen(nbfonc,ngauss,*)
 101       double precision vafott(nbfonc,ngauss,*)
 102 c
 103       integer ulsort, langue, codret
 104 c
 105 c 0.4. ==> variables locales
 106 c
 107 c     tecnp1 = TEtraedre courant en numerotation du Calcul a l'it. N+1
 108 c
 109       integer tecnp1
 110 c
 111 c     f1hp = Fils 1er du tetraedre en numerotation Homard a l'it. N+1
 112 c     f1cp = Fils 1er du tetraedre en numerota. du Calcul a l'it. N+1
 113 c     f2cp = Fils 2eme du tetraedre en numerota. du Calcul a l'it. N+1
 114 c     f3cp = Fils 3eme du tetraedre en numerota. du Calcul a l'it. N+1
 115 c     f4cp = Fils 4eme du tetraedre en numerota. du Calcul a l'it. N+1
 116 c
 117       integer f1hp, fihp
 118       integer f1cp, f2cp, f3cp, f4cp
 119 c
 120 c     f1hn = Fils 1er du tetraedre en numerotation Homard a l'it. N
 121 c     f1cn = Fils 1er du tetraedre en numerotation du Calcul a l'it. N
 122 c     f2cn = Fils 2eme du tetraedre en numerotation du Calcul a l'it. N
 123 c
 124       integer f1hn
 125       integer f1cn, f2cn
 126 c
 127       integer nrofon, nugaus
 128       integer coderr
 129 c
 130       integer iaux
 131       integer lglist, nrlist
 132       integer list(30)
 133 c
 134       double precision daux
 135       double precision daux1
 136 c
 137       integer nbmess
 138       parameter ( nbmess = 10 )
 139       character*80 texte(nblang,nbmess)
 140 c
 141 c 0.5. ==> initialisations
 142 c ______________________________________________________________________
 143 c
 144 c====
 145 c 1. initialisations
 146 c====
 147 c
 148 #include "impr01.h"
 149 c
 150 #ifdef _DEBUG_HOMARD_
 151       write (ulsort,texte(langue,1)) 'Entree', nompro
 152       call dmflsh (iaux)
 153 #endif
 154 c
 155       texte(1,4) =
 156      >'(/,''Tetr. en cours : numero a l''''iteration '',a3,'' : '',i10)'
 157       texte(1,5) =
 158      >'(  ''                  etat a l''''iteration '',a3,''   : '',i4)'
 159 c
 160       texte(2,4) =
 161      >'(/,''Current tetrahedron : # at iteration '',a3,''     : '',i10)'
 162       texte(2,5) =
 163      > '(  ''                     status at iteration '',a3,'' : '',i4)'
 164 c
 165       coderr = 0
 166 c
 167 c 1.2. ==> on repere les numeros dans le calcul pour ses deux fils
 168 c          a l'iteration n
 169 c
 170       f1hn = anfite(tehn)
 171       f1cn = nteeca(f1hn)
 172       f2cn = nteeca(f1hn+1)
 173 c
 174 c====
 175 c 2. etan = 21, ..., 26 : le tetraedre etait coupe en 2
 176 c    On explore tous les etats du tetraedre a l'iteration n+1
 177 c====
 178 c
 179       if ( prfcan(f1cn).gt.0 .and. prfcan(f2cn).gt.0 ) then
 180 c
 181 c ===> etanp1 = 0 : le tetraedre est actif
 182 c             Cela veut dire qu'il est reactive.
 183 c             on lui attribue la valeur moyenne sur les deux anciens
 184 c             fils.
 185 c             remarque : cela arrive seulement avec du deraffinement.
 186 c
 187       if ( etanp1.eq.0 ) then
 188 c
 189         tecnp1 = ntesca(tehnp1)
 190         prfcap(tecnp1) = 1
 191 c
 192         daux1 = 1.d0/dble(2*ngauss)
 193         do 21 , nrofon = 1, nbfonc
 194           daux = 0.d0
 195           do 211 , nugaus = 1, ngauss
 196             daux = daux + vafoen(nrofon,nugaus,prfcan(f1cn))
 197      >                  + vafoen(nrofon,nugaus,prfcan(f2cn))
 198   211     continue
 199           daux = daux*daux1
 200           do 212 , nugaus = 1 , ngauss
 201             vafott(nrofon,nugaus,tecnp1) = daux
 202   212     continue
 203    21   continue
 204 cgn        write(21,91010) tecnp1
 205 cgn        write(ulsort,91010) f1cn,f2cn,-1,tecnp1
 206 c
 207 c ===> etanp1 = etan : le tetraedre est decoupe en deux
 208 c                       selon le meme decoupage.
 209 c             c'est ce qui se passe quand un decoupage de conformite
 210 c             est supprime au debut des algorithmes d'adaptation,
 211 c             puis reproduit a la creation du maillage car les faces
 212 c             autour n'ont pas change entre les deux iterations.
 213 c             le fils prend la valeur de la fonction sur l'ancien
 214 c             fils qui etait au meme endroit. comme la procedure de
 215 c             numerotation est la meme (voir cmcdte), le premier fils
 216 c             est toujours le meme, le second egalement. on prendra
 217 c             alors la valeur sur le fils de rang identique a
 218 c             l'iteration n.
 219 c
 220       elseif ( etanp1.eq.etan ) then
 221 c
 222         f1hp = filtet(tehnp1)
 223         f1cp = ntesca(f1hp)
 224         f2cp = ntesca(f1hp+1)
 225         prfcap(f1cp) = 1
 226         prfcap(f2cp) = 1
 227         do 22 , nrofon = 1, nbfonc
 228           do 221 , nugaus = 1 , ngauss
 229             vafott(nrofon,nugaus,f1cp) =
 230      >                               vafoen(nrofon,nugaus,prfcan(f1cn))
 231             vafott(nrofon,nugaus,f2cp) =
 232      >                               vafoen(nrofon,nugaus,prfcan(f2cn))
 233   221     continue
 234    22   continue
 235 cgn        write(22,91010) f1cp,f2cp
 236 cgn        write(ulsort,91010) f1cn,f2cn,-1,f1cp,f2cp
 237 c
 238 c ===> etanp1 = 21, ..., 26 et different de
 239 c                             etan : le tetraedre est decoupe en deux
 240 c                             mais par un autre decoupage.
 241 c             c'est ce qui se passe quand un decoupage de conformite
 242 c             est supprime au debut des algorithmes d'adaptation. il y
 243 c             a du deraffinement dans la zone qui induisait le decoupage
 244 c             de conformite et raffinement sur une autre zone.
 245 c             on donne la valeur moyenne de la fonction sur les deux
 246 c             anciens fils a chaque nouveau fils.
 247 c             remarque : cela arrive seulement avec du deraffinement.
 248 c
 249       elseif ( etanp1.ge.21 .and. etanp1.le.26 ) then
 250 c
 251         f1hp = filtet(tehnp1)
 252         f1cp = ntesca(f1hp)
 253         f2cp = ntesca(f1hp+1)
 254         prfcap(f1cp) = 1
 255         prfcap(f2cp) = 1
 256         daux1 = 1.d0/dble(2*ngauss)
 257         do 23 , nrofon = 1, nbfonc
 258           daux = 0.d0
 259           do 231 , nugaus = 1, ngauss
 260             daux = daux + vafoen(nrofon,nugaus,prfcan(f1cn))
 261      >                  + vafoen(nrofon,nugaus,prfcan(f2cn))
 262   231     continue
 263           daux = daux*daux1
 264           do 232 , nugaus = 1 , ngauss
 265             vafott(nrofon,nugaus,f1cp) = daux
 266             vafott(nrofon,nugaus,f2cp) = daux
 267   232     continue
 268    23   continue
 269 cgn        write(23,91010) f1cp,f2cp
 270 cgn        write(ulsort,91010) f1cn,f2cn,-1,f1cp,f2cp
 271 c
 272 c ===> etanp1 = 41, 42, 43 ou 44 : le tetraedre est
 273 c                      decoupe en quatre par une face.
 274 c ===> etanp1 = 45, 46 ou 47 : le tetraedre est
 275 c                      decoupe en quatre par une diagonale.
 276 c             c'est ce qui se passe quand un decoupage de conformite
 277 c             est supprime au debut des algorithmes d'adaptation, puis
 278 c             remis differement car l'environnement a change.
 279 c             on donne la valeur moyenne de la fonction sur les deux
 280 c             anciens fils a chaque nouveau fils.
 281 c             remarque : on pourrait certainement faire mieux, avec des
 282 c                        moyennes ponderees en fonction du recouvrement
 283 c                        des anciens et nouveaux fils. c'est trop
 284 c                        complique pour que cela vaille le coup.
 285 c
 286       elseif ( etanp1.ge.41 .and. etanp1.le.47 ) then
 287 c
 288         f1hp = filtet(tehnp1)
 289         f1cp = ntesca(f1hp)
 290         f2cp = ntesca(f1hp+1)
 291         f3cp = ntesca(f1hp+2)
 292         f4cp = ntesca(f1hp+3)
 293         prfcap(f1cp) = 1
 294         prfcap(f2cp) = 1
 295         prfcap(f3cp) = 1
 296         prfcap(f4cp) = 1
 297         daux1 = 1.d0/dble(2*ngauss)
 298         do 24 , nrofon = 1, nbfonc
 299           daux = 0.d0
 300           do 241 , nugaus = 1, ngauss
 301             daux = daux + vafoen(nrofon,nugaus,prfcan(f1cn))
 302      >                  + vafoen(nrofon,nugaus,prfcan(f2cn))
 303   241     continue
 304           daux = daux*daux1
 305           do 242 , nugaus = 1 , ngauss
 306             vafott(nrofon,nugaus,f1cp) = daux
 307             vafott(nrofon,nugaus,f2cp) = daux
 308             vafott(nrofon,nugaus,f3cp) = daux
 309             vafott(nrofon,nugaus,f4cp) = daux
 310   242     continue
 311    24   continue
 312 cgn        write(24,91010) f1cp,f2cp,f3cp,f4cp
 313 cgn        write(ulsort,91010) f1cn,f2cn,-1,f1cp,f2cp,f3cp,f4cp
 314 c
 315 c ===> etanp1 = 85, 86 ou 87 : le tetraedre est
 316 c                      decoupe en huit par une diagonale.
 317 c             c'est ce qui se passe quand un decoupage de conformite
 318 c             est supprime au debut des algorithmes d'adaptation, puis
 319 c             remis differement car l'environnement a change.
 320 c             on donne la valeur moyenne de la fonction sur les deux
 321 c             anciens fils a chaque nouveau fils.
 322 c             attention : il est possible que les fils sur les bords
 323 c                         soient decoupes par de la conformite. Il faut
 324 c                         alors transmettre la valeur a leurs 2 ou 4
 325 c                         fils.
 326 c             attention : ce n'est pas comme en 2D ; il faut examiner
 327 c                         tous les fils, car par contamination de faces
 328 c                         coupees en 2, les fils centraux peuvent etre
 329 c                         decoupes.
 330 c             remarque : on pourrait certainement faire mieux, avec des
 331 c                        moyennes ponderees en fonction du recouvrement
 332 c                        des anciens et nouveaux fils. c'est trop
 333 c                        complique pour que cela vaille le coup.
 334 c
 335       elseif ( etanp1.ge.85 .and. etanp1.le.87 ) then
 336 c
 337         f1hp = filtet(tehnp1)
 338         lglist = 0
 339         do 251 , nrlist = 1 , 8
 340           fihp = f1hp+nrlist-1
 341           iaux = mod(hettet(fihp),100)
 342           if ( iaux.eq.0 ) then
 343             lglist = lglist + 1
 344             list(lglist) = ntesca(fihp)
 345           elseif ( iaux.ge.21 .and. iaux.le.26 ) then
 346             lglist = lglist + 1
 347             list(lglist) = ntesca(filtet(fihp))
 348             lglist = lglist + 1
 349             list(lglist) = ntesca(filtet(fihp)+1)
 350           elseif ( iaux.ge.41 .and. iaux.le.47 ) then
 351             lglist = lglist + 1
 352             list(lglist) = ntesca(filtet(fihp))
 353             lglist = lglist + 1
 354             list(lglist) = ntesca(filtet(fihp)+1)
 355             lglist = lglist + 1
 356             list(lglist) = ntesca(filtet(fihp)+2)
 357             lglist = lglist + 1
 358             list(lglist) = ntesca(filtet(fihp)+3)
 359           else
 360             coderr = 1
 361           endif
 362   251   continue
 363 c
 364         do 252 , nrlist = 1 , lglist
 365           prfcap(list(nrlist)) = 1
 366   252   continue
 367 c
 368         daux1 = 1.d0/dble(2*ngauss)
 369         do 25 , nrofon = 1, nbfonc
 370           daux = 0.d0
 371           do 253 , nugaus = 1, ngauss
 372             daux = daux + vafoen(nrofon,nugaus,prfcan(f1cn))
 373      >                  + vafoen(nrofon,nugaus,prfcan(f2cn))
 374   253     continue
 375           daux = daux*daux1
 376           do 254 , nugaus = 1 , ngauss
 377             do 255 , nrlist = 1 , lglist
 378               vafott(nrofon,nugaus,list(nrlist)) = daux
 379   255       continue
 380   254     continue
 381    25   continue
 382 cgn        write(26,91010) (list(nrlist),nrlist = 1 , lglist)
 383 cgn        write(ulsort,91010) f1cn,f2cn,-1,
 384 cgn     >                     (list(nrlist),nrlist = 1 , lglist)
 385 c
 386 c ==> aucun autre etat sur le tetraedre courant n'est possible
 387 c
 388       else
 389 c
 390         coderr = 1
 391         write (ulsort,texte(langue,4)) 'n  ', tehn
 392         write (ulsort,texte(langue,5)) 'n  ', etan
 393         write (ulsort,texte(langue,4)) 'n+1', tehnp1
 394         write (ulsort,texte(langue,5)) 'n+1', etanp1
 395 c
 396       endif
 397 c
 398       endif
 399 c
 400 c====
 401 c 3. la fin
 402 c====
 403 c
 404       if ( coderr.ne.0 ) then
 405 c
 406         write (ulsort,texte(langue,1)) 'Sortie', nompro
 407       write (ulsort,texte(langue,2)) coderr
 408       codret = codret + 1
 409 c
 410       endif
 411 c
 412       end