src/tool/AP_Conversion/pcseq3.F

   1       subroutine pcseq3 ( etanp1, quhn, quhnp1, typint,
   2      >                    prfcan, prfcap,
   3      >                    coonoe,
   4      >                    somare,
   5      >                    arequa, filqua,
   6      >                    nbanqu, anfiqu,
   7      >                    nqueca, nqusca,
   8      >                    aretri,
   9      >                            ntrsca,
  10      >                    nbfonc, vafoen, vafott,
  11      >                            vatrtt,
  12      >                            prftrp,
  13      >                    ulsort, langue, codret )
  14 c ______________________________________________________________________
  15 c
  16 c                             H O M A R D
  17 c
  18 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
  19 c
  20 c Version originale enregistree le 18 juin 1996 sous le numero 96036
  21 c aupres des huissiers de justice Simart et Lavoir a Clamart
  22 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
  23 c aupres des huissiers de justice
  24 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
  25 c
  26 c    HOMARD est une marque deposee d'Electricite de France
  27 c
  28 c Copyright EDF 1996
  29 c Copyright EDF 1998
  30 c Copyright EDF 2002
  31 c Copyright EDF 2020
  32 c ______________________________________________________________________
  33 c
  34 c    aPres adaptation - Conversion de Solution Elements de volume -
  35 c     -                 -             -        -
  36 c                       Quadrangles d'etat anterieur 4
  37 c                       -
  38 c ______________________________________________________________________
  39 c .        .     .        .                                            .
  40 c .  nom   . e/s . taille .           description                      .
  41 c .____________________________________________________________________.
  42 c . etanp1 . e   .    1   . ETAt du quadrangle a l'iteration N+1       .
  43 c . quhn   . e   .    1   . Quadrangle courant en numerotation Homard  .
  44 c .        .     .        . a l'iteration N                            .
  45 c . quhnp1 . e   .    1   . Quadrangle courant en numerotation Homard  .
  46 c .        .     .        . a l'iteration N+1                          .
  47 c . typint . e   .   1    . type d'interpolation                       .
  48 c .        .     .        .  0, si automatique                         .
  49 c .        .     .        .  elements : 0 si intensif, sans orientation.
  50 c .        .     .        .             1 si extensif, sans orientation.
  51 c .        .     .        .             2 si intensif, avec orientation.
  52 c .        .     .        .             3 si extensif, avec orientation.
  53 c .        .     .        .  noeuds : 1 si degre 1                     .
  54 c .        .     .        .           2 si degre 2                     .
  55 c .        .     .        .           3 si iso-P2                      .
  56 c . typint . e   .   1    . type d'interpolation                       .
  57 c .        .     .        .  0, si automatique                         .
  58 c .        .     .        .  elements : 0 si intensif, sans orientation.
  59 c .        .     .        .             1 si extensif, sans orientation.
  60 c .        .     .        .             2 si intensif, avec orientation.
  61 c .        .     .        .             3 si extensif, avec orientation.
  62 c .        .     .        .  noeuds : 1 si degre 1                     .
  63 c .        .     .        .           2 si degre 2                     .
  64 c .        .     .        .           3 si iso-P2                      .
  65 c . deraff . e   .    1   . vrai, s'il y a eu du deraffinement en      .
  66 c .        .     .        . passant de l'iteration n a n+1 ; faux sinon.
  67 c . prfcan . e   .   *    . En numero du calcul a l'iteration n :      .
  68 c .        .     .        . 0 : l'entite est absente du profil         .
  69 c .        .     .        . i : l'entite est au rang i dans le profil  .
  70 c . prfcap . es  .   *    . En numero du calcul a l'iteration n+1 :    .
  71 c .        .     .        . 0 : l'entite est absente du profil         .
  72 c .        .     .        . 1 : l'entite est presente dans le profil   .
  73 c . coonoe . e   . nbnoto . coordonnees des noeuds                     .
  74 c .        .     . * sdim .                                            .
  75 c . somare . e   .2*nbarto. numeros des extremites d'arete             .
  76 c . arequa . e   .nbquto*4. numeros des 4 aretes des quadrangles       .
  77 c . filqua . e   . nbquto . premier fils des quadrangles               .
  78 c . nbanqu . e   .   1    . nombre de quadrangles decoupes par         .
  79 c .        .     .        . conformite sur le maillage avant adaptation.
  80 c . anfiqu . e   . nbanqu . tableau filqua du maillage de l'iteration n.
  81 c . nqueca . e   .   *    . nro des quadrangles dans le calcul en ent. .
  82 c . nqusca . e   . rsquto . numero des quadrangles du calcul           .
  83 c . aretri . e   .nbtrto*3. numeros des 3 aretes des triangles         .
  84 c . ntrsca . e   . rstrto . numero des triangles du calcul             .
  85 c . nbfonc . e   .    1   . nombre de fonctions elements de volume     .
  86 c . vafoen . e   . nbfonc*. variables en entree de l'adaptation        .
  87 c .        .     .    *   .                                            .
  88 c . vafott .   a . nbfonc*. tableau temporaire de la solution          .
  89 c .        .     .    *   .                                            .
  90 c . vatrtt .   a . nbfonc*. tableau temporaire de la solution pour     .
  91 c .        .     .    *   . les triangles de conformite                .
  92 c . prftrp . es  .   *    . En numero du calcul a l'iteration n+1 :    .
  93 c .        .     .        . 0 : le triangle est absent du profil       .
  94 c .        .     .        . 1 : le triangle est present dans le profil .
  95 c . ulsort . e   .   1    . numero d'unite logique de la liste standard.
  96 c . langue . e   .    1   . langue des messages                        .
  97 c .        .     .        . 1 : francais, 2 : anglais                  .
  98 c . codret . es  .    1   . code de retour des modules                 .
  99 c .        .     .        . 0 : pas de probleme                        .
 100 c .        .     .        . 1 : probleme                               .
 101 c ______________________________________________________________________
 102 c
 103 c====
 104 c 0. declarations et dimensionnement
 105 c====
 106 c
 107 c 0.1. ==> generalites
 108 c
 109       implicit none
 110       save
 111 c
 112       character*6 nompro
 113       parameter ( nompro = 'PCSEQ3' )
 114 c
 115 #include "nblang.h"
 116 #include "fractc.h"
 117 c
 118 c 0.2. ==> communs
 119 c
 120 #include "envca1.h"
 121 #include "nombno.h"
 122 #include "nombar.h"
 123 #include "nombtr.h"
 124 #include "nombqu.h"
 125 #include "nombsr.h"
 126 #include "nomber.h"
 127 c
 128 c 0.3. ==> arguments
 129 c
 130       integer etanp1, quhn, quhnp1
 131       integer typint
 132       integer nbfonc
 133       integer prfcan(*), prfcap(*)
 134       integer somare(2,nbarto)
 135       integer arequa(nbquto,4), filqua(nbquto)
 136       integer nbanqu, anfiqu(nbanqu)
 137       integer nqueca(requto), nqusca(rsquto)
 138       integer aretri(nbtrto,3)
 139       integer ntrsca(rstrto)
 140       integer prftrp(*)
 141 c
 142       double precision coonoe(nbnoto,sdim)
 143       double precision vafoen(nbfonc,*)
 144       double precision vafott(nbfonc,*)
 145       double precision vatrtt(nbfonc,*)
 146 c
 147       integer ulsort, langue, codret
 148 c
 149 c 0.4. ==> variables locales
 150 c
 151 c     qucnp1 = QUadrangle courant en numerotation Calcul a l'it. N+1
 152 c
 153       integer qucnp1
 154 c
 155 c     f1hp = Fils 1er du quadrangle en numerotation Homard a l'it. N+1
 156 c     f1cp = Fils 1er du quadrangle en numerotation Calcul a l'it. N+1
 157 c     f2cp = Fils 2eme du quadrangle en numerotation Calcul a l'it. N+1
 158 c     f3cp = Fils 3eme du quadrangle en numerotation Calcul a l'it. N+1
 159 c
 160       integer f1hp
 161       integer f1cp, f2cp, f3cp
 162 c
 163 c     f1hn = Fils 1er du quadrangle en numerotation Homard a l'it. N
 164 c     f1cn = Fils 1er du quadrangle en numerotation Calcul a l'it. N
 165 c     f2cn = Fils 2eme du quadrangle en numerotation Calcul a l'it. N
 166 c     f3cn = Fils 3eme du quadrangle en numerotation Calcul a l'it. N
 167 c     f4cn = Fils 4eme du quadrangle en numerotation Calcul a l'it. N
 168 c
 169       integer f1hn
 170       integer f1cn, f2cn, f3cn, f4cn
 171 c
 172       integer coderr
 173       integer nrofon
 174       integer iaux
 175 c
 176       double precision daux
 177       double precision daux0, daux1, daux2, daux3
 178 c
 179       integer nbmess
 180       parameter ( nbmess = 10 )
 181       character*80 texte(nblang,nbmess)
 182 c
 183 c 0.5. ==> initialisations
 184 c ______________________________________________________________________
 185 c
 186 c====
 187 c 1. initialisations
 188 c====
 189 c
 190 #include "pcimp0.h"
 191 #include "impr01.h"
 192 #include "impr03.h"
 193 c
 194 #ifdef _DEBUG_HOMARD_
 195       write (ulsort,texte(langue,1)) 'Entree', nompro
 196       call dmflsh (iaux)
 197 #endif
 198 c
 199       texte(1,4) =
 200      >'(''Quad. en cours : nro a l''''iteration '',a3,'' :'',i10)'
 201       texte(1,5) = '( 16x,''etat a l''''iteration '',a3,'' : '',i4)'
 202 c
 203       texte(2,4) =
 204      >'(''Current quadrangle : # at iteration '',a3,''     : '',i10)'
 205       texte(2,5) = '( 17x,''status at iteration '',a3,'' : '',i4)'
 206 c
 207       coderr = 0
 208 c
 209 c 1.2. ==> on repere les numeros dans le calcul pour les fils
 210 c          a l'iteration n
 211 c
 212       f1hn = anfiqu(quhn)
 213       f1cn = nqueca(f1hn)
 214       f2cn = nqueca(f1hn+1)
 215       f3cn = nqueca(f1hn+2)
 216       f4cn = nqueca(f1hn+3)
 217 c
 218       if ( prfcan(f1cn).gt.0 .and. prfcan(f2cn).gt.0 .and.
 219      >     prfcan(f3cn).gt.0 .and. prfcan(f3cn).gt.0 ) then
 220 c
 221 c====
 222 c 2. le quadrangle etait coupe en 4 quadrangles
 223 c====
 224 c 2.1. ==> etanp1 = 0 : le quadrangle est actif ; il est reactive.
 225 c          on lui attribue la valeur moyenne des quatre anciens fils.
 226 c          remarque : cela arrive seulement avec du deraffinement.
 227 c           .................         .................
 228 c           .       .       .         .               .
 229 c           .       .       .         .               .
 230 c           .       .       .         .               .
 231 c           .................  ===>   .               .
 232 c           .       .       .         .               .
 233 c           .       .       .         .               .
 234 c           .       .       .         .               .
 235 c           .................         .................
 236 c
 237       if ( etanp1.eq.0 ) then
 238 cgn      write(ulsort,*)'... quadrangle reactive'
 239 c
 240         qucnp1 = nqusca(quhnp1)
 241         prfcap(qucnp1) = 1
 242 c
 243         if ( typint.eq.0 ) then
 244           daux1 = unsqu
 245         else
 246           daux1 = 1.d0
 247         endif
 248         do 21 , nrofon = 1 , nbfonc
 249           daux = daux1 * ( vafoen(nrofon,prfcan(f1cn))
 250      >                   + vafoen(nrofon,prfcan(f2cn))
 251      >                   + vafoen(nrofon,prfcan(f3cn))
 252      >                   + vafoen(nrofon,prfcan(f4cn)) )
 253 cgn            write(ulsort,90004) 'daux', daux
 254           vafott(nrofon,qucnp1) = daux
 255    21   continue
 256 c
 257 c 2.2. ==> etanp1 = 21/22 : le quadrangle est decoupe en
 258 c          deux quadrangles
 259 c             On donne la valeur moyenne de la fonction sur les deux
 260 c             anciens fils a chaque nouveau fils.
 261 c             remarque : on pourrait certainement faire mieux, avec des
 262 c                        moyennes ponderees en fonction du recouvrement
 263 c                        des anciens et nouveaux fils. c'est trop
 264 c                        complique pour que cela vaille le coup.
 265 c           .................         .................
 266 c           .       .       .         .       .       .
 267 c           .       .       .         .       .       .
 268 c           .       .       .         .       .       .
 269 c           .................  ===>   .       .       .
 270 c           .       .       .         .       .       .
 271 c           .       .       .         .       .       .
 272 c           .       .       .         .       .       .
 273 c           .................         .................
 274       elseif ( etanp1.ge.21 .and. etanp1.le.22 ) then
 275 cgn      write(ulsort,*)'... quadrangle coupe en 2'
 276 c
 277         f1hp = filqua(quhnp1)
 278         f1cp = nqusca(f1hp)
 279         f2cp = nqusca(f1hp+1)
 280         prfcap(f1cp) = 1
 281         prfcap(f2cp) = 1
 282 c
 283         if ( typint.eq.0 ) then
 284           do 221 , nrofon = 1 , nbfonc
 285             daux = unsqu * ( vafoen(nrofon,prfcan(f1cn))
 286      >                     + vafoen(nrofon,prfcan(f2cn))
 287      >                     + vafoen(nrofon,prfcan(f3cn))
 288      >                     + vafoen(nrofon,prfcan(f4cn)) )
 289             vafott(nrofon,f1cp) = daux
 290             vafott(nrofon,f2cp) = daux
 291   221     continue
 292         else
 293           call utqqua ( f1hp  , daux, daux1, coonoe, somare, arequa )
 294           call utqqua ( f1hp+1, daux, daux2, coonoe, somare, arequa )
 295           daux0 = daux1 + daux2
 296           daux1 = daux1 / daux0
 297           daux2 = daux2 / daux0
 298           do 222 , nrofon = 1 , nbfonc
 299             daux = vafoen(nrofon,prfcan(f1cn))
 300      >           + vafoen(nrofon,prfcan(f2cn))
 301      >           + vafoen(nrofon,prfcan(f3cn))
 302      >           + vafoen(nrofon,prfcan(f4cn))
 303             vafott(nrofon,f1cp) = daux * daux1
 304             vafott(nrofon,f2cp) = daux * daux2
 305   222     continue
 306         endif
 307 c
 308 c 2.3. ==> etanp1 = 31, 32, 33 ou 34 : le quadrangle est
 309 c             decoupe en trois triangles.
 310 c             on attribue la valeur moyenne sur les quatre anciens
 311 c             fils a chacune des trois nouveaux fils.
 312 c             remarque : on pourrait certainement faire mieux, avec des
 313 c                        moyennes ponderees en fonction du recouvrement
 314 c                        des anciens et nouveaux fils. c'est trop
 315 c                        complique pour que cela vaille le coup.
 316 c             remarque : cela arrive seulement avec du deraffinement.
 317 c           .................         .................
 318 c           .       .       .         .      . .      .
 319 c           .       .       .         .     .   .     .
 320 c           .       .       .         .    .     .    .
 321 c           .................  ===>   .   .       .   .
 322 c           .       .       .         .  .         .  .
 323 c           .       .       .         . .           . .
 324 c           .       .       .         ..             ..
 325 c           .................         .................
 326 c
 327       elseif ( etanp1.ge.31 .and. etanp1.le.34 ) then
 328 cgn      write(ulsort,*)'... quadrangle coupe en 3 triangles'
 329 c
 330         f1hp = -filqua(quhnp1)
 331         f1cp = ntrsca(f1hp)
 332         f2cp = ntrsca(f1hp+1)
 333         f3cp = ntrsca(f1hp+2)
 334         prftrp(f1cp) = 1
 335         prftrp(f2cp) = 1
 336         prftrp(f3cp) = 1
 337         if ( typint.eq.0 ) then
 338           do 231 , nrofon = 1 , nbfonc
 339             daux = unsqu * ( vafoen(nrofon,prfcan(f1cn))
 340      >                     + vafoen(nrofon,prfcan(f2cn))
 341      >                     + vafoen(nrofon,prfcan(f3cn))
 342      >                     + vafoen(nrofon,prfcan(f4cn)) )
 343             vatrtt(nrofon,f1cp) = daux
 344             vatrtt(nrofon,f2cp) = daux
 345             vatrtt(nrofon,f3cp) = daux
 346   231     continue
 347         else
 348           call utqtri ( f1hp  , daux, daux1, coonoe, somare, aretri )
 349           call utqtri ( f1hp+1, daux, daux2, coonoe, somare, aretri )
 350           call utqtri ( f1hp+2, daux, daux3, coonoe, somare, aretri )
 351           daux0 = daux1 + daux2 + daux3
 352           daux1 = daux1 / daux0
 353           daux2 = daux2 / daux0
 354           daux3 = daux3 / daux0
 355           do 232 , nrofon = 1 , nbfonc
 356             daux = vafoen(nrofon,prfcan(f1cn))
 357      >           + vafoen(nrofon,prfcan(f2cn))
 358      >           + vafoen(nrofon,prfcan(f3cn))
 359      >           + vafoen(nrofon,prfcan(f4cn))
 360             vatrtt(nrofon,f1cp) = daux * daux1
 361             vatrtt(nrofon,f2cp) = daux * daux2
 362             vatrtt(nrofon,f3cp) = daux * daux3
 363   232     continue
 364         endif
 365 c
 366 c 2.4. ==> etanp1 = 41, 42, 43 ou 43 : le quadrangle est decoupe en
 367 c          trois quadrangles
 368 c             c'est ce qui se passe quand un decoupage de conformite
 369 c             est supprime au debut des algorithmes d'adaptation. il y
 370 c             a du deraffinement dans la zone qui induisait le decoupage
 371 c             de conformite et raffinement sur une autre zone.
 372 c             On donne la valeur moyenne de la fonction sur les trois
 373 c             anciens fils a chaque nouveau fils.
 374 c             remarque : on pourrait certainement faire mieux, avec des
 375 c                        moyennes ponderees en fonction du recouvrement
 376 c                        des anciens et nouveaux fils. c'est trop
 377 c                        complique pour que cela vaille le coup.
 378 c             remarque : cela arrive seulement avec du deraffinement.
 379 c
 380 c           .................         .................
 381 c           .       .       .         .       .       .
 382 c           .       .       .         .       .       .
 383 c           .       .       .         .       .       .
 384 c           .................  ===>   .........       .
 385 c           .       .       .         .         .     .
 386 c           .       .       .         .           .   .
 387 c           .       .       .         .             . .
 388 c           .................         .................
 389 c
 390       elseif ( etanp1.ge.41 .and. etanp1.le.44 ) then
 391 cgn         print *,'... le quadrangle est coupe en 3 quadrangles'
 392 c
 393         f1hp = filqua(quhnp1)
 394         f1cp = nqusca(f1hp)
 395         f2cp = nqusca(f1hp+1)
 396         f3cp = nqusca(f1hp+2)
 397         prfcap(f1cp) = 1
 398         prfcap(f2cp) = 1
 399         prfcap(f3cp) = 1
 400 c
 401         if ( typint.eq.0 ) then
 402 c
 403           do 241 , nrofon = 1 , nbfonc
 404             daux = unsqu * ( vafoen(nrofon,prfcan(f1cn))
 405      >                     + vafoen(nrofon,prfcan(f2cn))
 406      >                     + vafoen(nrofon,prfcan(f3cn))
 407      >                     + vafoen(nrofon,prfcan(f4cn)) )
 408 cgn            write(ulsort,90004) 'daux', daux
 409             vafott(nrofon,f1cp) = daux
 410             vafott(nrofon,f2cp) = daux
 411             vafott(nrofon,f3cp) = daux
 412   241     continue
 413 c
 414         else
 415 c
 416           call utqqua ( f1hp  , daux, daux1, coonoe, somare, arequa )
 417           call utqqua ( f1hp+1, daux, daux2, coonoe, somare, arequa )
 418           call utqqua ( f1hp+2, daux, daux3, coonoe, somare, arequa )
 419           daux0 = daux1 + daux2 + daux3
 420           daux1 = daux1 / daux0
 421           daux2 = daux2 / daux0
 422           daux3 = daux3 / daux0
 423           do 242 , nrofon = 1 , nbfonc
 424             daux = vafoen(nrofon,prfcan(f1cn))
 425      >           + vafoen(nrofon,prfcan(f2cn))
 426      >           + vafoen(nrofon,prfcan(f3cn))
 427      >           + vafoen(nrofon,prfcan(f4cn))
 428 cgn            write(ulsort,90004) 'unsqu*daux', unsqu*daux
 429 cgn            write(ulsort,90004) 'trshu*daux', trshu*daux
 430             vafott(nrofon,f1cp) = daux * daux1
 431             vafott(nrofon,f2cp) = daux * daux2
 432             vafott(nrofon,f3cp) = daux * daux3
 433   242     continue
 434 c
 435         endif
 436 c
 437       endif
 438 c
 439       endif
 440 c
 441 c====
 442 c 3. la fin
 443 c====
 444 c
 445       if ( coderr.ne.0 ) then
 446 c
 447       write (ulsort,texte(langue,1)) 'Sortie', nompro
 448       write (ulsort,texte(langue,2)) coderr
 449       codret = codret + coderr
 450 c
 451       endif
 452 c
 453 #ifdef _DEBUG_HOMARD_
 454       write (ulsort,texte(langue,1)) 'Sortie', nompro
 455       call dmflsh (iaux)
 456 #endif
 457 c
 458       end