src/tool/AP_Conversion/pcs2qu.F

   1       subroutine pcs2qu ( nbfop2, profho, vap2ho,
   2      >                    hetqua, arequa, filqua,
   3      >                    somare, np2are,
   4      >                    aretri )
   5 c ______________________________________________________________________
   6 c
   7 c                             H O M A R D
   8 c
   9 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
  10 c
  11 c Version originale enregistree le 18 juin 1996 sous le numero 96036
  12 c aupres des huissiers de justice Simart et Lavoir a Clamart
  13 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
  14 c aupres des huissiers de justice
  15 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
  16 c
  17 c    HOMARD est une marque deposee d'Electricite de France
  18 c
  19 c Copyright EDF 1996
  20 c Copyright EDF 1998
  21 c Copyright EDF 2002
  22 c Copyright EDF 2020
  23 c ______________________________________________________________________
  24 c
  25 c    aPres adaptation - Conversion de Solution -
  26 c     -                 -             -
  27 c    interpolation p2 sur les noeuds - decoupage des QUadrangles
  28 c                   -                                --
  29 c remarque : on devrait optimiser cela car si le quadrangle etait dans
  30 c            un etat de decoupage similaire, on recalcule une valeur
  31 c            qui est deja presente
  32 c remarque : pcs2qu et pcsiqu sont des clones
  33 c ______________________________________________________________________
  34 c .        .     .        .                                            .
  35 c .  nom   . e/s . taille .           description                      .
  36 c .____________________________________________________________________.
  37 c . nbfop2 . e   .    1   . nombre de fonctions P2                     .
  38 c . profho . es  .   *    . pour chaque entite en numerotation homard :.
  39 c .        .     .        . 0 : l'entite est absente du profil         .
  40 c .        .     .        . 1 : l'entite est presente dans le profil   .
  41 c . vap2ho . es  . nbfop2*. variables p2 numerotation homard           .
  42 c .        .     . nbnoto .                                            .
  43 c . hetqua . e   . nbquto . historique de l'etat des quadrangles       .
  44 c . arequa . e   .nbquto*4. numeros des 4 aretes des quadrangles       .
  45 c . filqua . e   . nbquto . premier fils des quadrangles               .
  46 c . somare . e   .2*nbarto. numeros des extremites d'arete             .
  47 c . np2are . e   . nbarto . numero des noeuds p2 milieux d'aretes      .
  48 c . aretri . e   .nbtrto*3. numeros des 3 aretes des triangles         .
  49 c ______________________________________________________________________
  50 c
  51 c====
  52 c 0. declarations et dimensionnement
  53 c====
  54 c
  55 c 0.1. ==> generalites
  56 c
  57       implicit none
  58       save
  59 c
  60 #include "fracta.h"
  61 #include "fractc.h"
  62 #include "fractf.h"
  63 #include "fractg.h"
  64 c
  65 c 0.2. ==> communs
  66 c
  67 #include "nombno.h"
  68 #include "nombar.h"
  69 #include "nombtr.h"
  70 #include "nombqu.h"
  71 c
  72 c 0.3. ==> arguments
  73 c
  74       integer nbfop2
  75       integer profho(nbnoto)
  76       integer hetqua(nbquto), arequa(nbquto,4), filqua(nbquto)
  77       integer somare(2,nbarto), np2are(nbarto)
  78       integer aretri(nbtrto,3)
  79 c
  80       double precision vap2ho(nbfop2,*)
  81 c
  82 c 0.4. ==> variables locales
  83 c
  84       integer iaux
  85       integer lequad
  86       integer typdec, etanp1
  87       integer iaux1, iaux2, iaux3, iaux4
  88       integer jaux1,        jaux3, jaux4
  89       integer m1, m2, m3, m4
  90       integer sm, nuv
  91 c
  92       integer listar(4), listno(8)
  93 c
  94 c     f1hp = Fils 1er du quadrangle en numerotation Homard a l'it. N+1
  95       integer f1hp
  96 c
  97       logical afaire
  98 c
  99       double precision daux
 100 c
 101 c 0.5. ==> initialisations
 102 c ______________________________________________________________________
 103 cgn       write (*,*) 'PCS2QU'
 104 c
 105       do 10 , lequad = 1, nbquto
 106 c
 107 c====
 108 c 1. interpolation p2 pour un quadrangle qui vient d'etre decoupe :
 109 c    on a une valeur a mettre sur l'eventuel noeud central et les
 110 c    milieux des aretes internes
 111 c====
 112 c
 113         iaux = lequad
 114         call pcs0qu ( iaux, profho,
 115      >                hetqua, arequa,
 116      >                somare, np2are,
 117      >                afaire, listar, listno, typdec, etanp1 )
 118 c
 119         if ( afaire ) then
 120 cgn        write(1,*) 'typdec, etanp1 =', typdec, etanp1
 121 c
 122 c====
 123 c 2. La valeur sur le noeud au centre du quadrangle
 124 c    . Soit le quadrangle vient d'etre decoupe en 4 quadrangles alors
 125 c      qu'il ne l'etait pas a l'iteration precedente
 126 c    . Soit le quadrangle vient d'etre decoupe en 3 quadrangles
 127 c    . Soit le quadrangle vient d'etre decoupe en 2 quadrangles
 128 c    Remarque : regarder cmcdqu pour les conventions
 129 c====
 130 c
 131           if ( typdec.eq.4 .or.
 132      >         typdec.eq.21 .or. typdec.eq.22 .or.
 133      >         ( typdec.ge.41 .and. typdec.le.44 ) ) then
 134 c
 135             f1hp = filqua(lequad)
 136 c
 137 c 2.1. ==> la valeur sur le noeud au centre du quadrangle
 138 c
 139             if ( typdec.eq.4 ) then
 140               sm = somare(2,arequa(f1hp,2))
 141             elseif ( typdec.eq.21 .or. typdec.eq.22 ) then
 142               sm = np2are(arequa(f1hp,4))
 143             else
 144               sm = somare(2,arequa(f1hp,3))
 145             endif
 146             profho(sm) = 1
 147 cgn        write(1,*) 'f1hp =', f1hp
 148             if ( typdec.ne.4 ) then
 149 cgn        write(1,*) 'sm =', sm
 150             endif
 151 c
 152 c         interpolation = -1/4 (u1+u2+u3+u4) + 1/2 (u5+u6+u7+u8)
 153 c
 154  1789 format( 4g13.5)
 155             do 21, nuv = 1, nbfop2
 156 cgn        write(1,1789) vap2ho(nuv,listno(1))
 157 cgn     >                                   , vap2ho(nuv,listno(2))
 158 cgn     >                                   , vap2ho(nuv,listno(3))
 159 cgn     >                                   , vap2ho(nuv,listno(4))
 160 cgn        write(1,1789) vap2ho(nuv,listno(5))
 161 cgn     >                                   , vap2ho(nuv,listno(6))
 162 cgn     >                                   , vap2ho(nuv,listno(7))
 163 cgn     >                                   , vap2ho(nuv,listno(8))
 164 c
 165               vap2ho(nuv,sm) = - unsqu * ( vap2ho(nuv,listno(1))
 166      >                                   + vap2ho(nuv,listno(2))
 167      >                                   + vap2ho(nuv,listno(3))
 168      >                                   + vap2ho(nuv,listno(4)) )
 169      >                         + unsde * ( vap2ho(nuv,listno(5))
 170      >                                   + vap2ho(nuv,listno(6))
 171      >                                   + vap2ho(nuv,listno(7))
 172      >                                   + vap2ho(nuv,listno(8)) )
 173             if ( typdec.ne.4 ) then
 174 cgn        write(1,1789) vap2ho(nuv,sm)
 175             endif
 176 c
 177    21       continue
 178 c
 179           endif
 180 c
 181 c====
 182 c 3. Les valeurs sur les noeuds au milieu des aretes tracees a
 183 c    l'interieur du quadrangle
 184 c====
 185 c 3.1. ==> Le quadrangle vient d'etre decoupe en 4 quadrangles alors
 186 c          qu'il ne l'etait pas a l'iteration precedente
 187 c          Remarque : regarder cmrdqu pour les conventions
 188 c
 189           if ( typdec.eq.4 ) then
 190 c
 191             f1hp = filqua(lequad)
 192 c
 193 c            Par convention, la deuxieme arete du i-eme fils va du
 194 c            noeud ni, milieu de l'arete ai du quadrangle pere, vers
 195 c            le noeud central.
 196 c
 197             m1 = np2are(arequa(f1hp,2))
 198             m2 = np2are(arequa(f1hp+1,2))
 199             m3 = np2are(arequa(f1hp+2,2))
 200             m4 = np2are(arequa(f1hp+3,2))
 201             profho(m1) = 1
 202             profho(m2) = 1
 203             profho(m3) = 1
 204             profho(m4) = 1
 205 c
 206 c         interpolation = -3/16 (u1+u2+u3+u4)
 207 c                        + 3/4 u5 +3/8 (u6+u8) + 1/4 u7
 208 c         avec u5 pour le noeud le plus proche, u6 et u8 pour ceux qui
 209 c         'encadrent' et u7 pour le noeud 'oppose'
 210 c
 211             do 31, nuv = 1, nbfop2
 212 c
 213               daux = - trssz * ( vap2ho(nuv,listno(1))
 214      >                         + vap2ho(nuv,listno(2))
 215      >                         + vap2ho(nuv,listno(3))
 216      >                         + vap2ho(nuv,listno(4)) )
 217 c
 218               vap2ho(nuv,m1) = daux
 219      >        + trsqu * vap2ho(nuv,listno(5))
 220      >        + trshu * ( vap2ho(nuv,listno(6)) + vap2ho(nuv,listno(8)))
 221      >        + unsqu *  vap2ho(nuv,listno(7))
 222 c
 223               vap2ho(nuv,m2) = daux
 224      >        + trsqu * vap2ho(nuv,listno(6))
 225      >        + trshu * ( vap2ho(nuv,listno(7)) + vap2ho(nuv,listno(5)))
 226      >        + unsqu *  vap2ho(nuv,listno(8))
 227 c
 228               vap2ho(nuv,m3) = daux
 229      >        + trsqu * vap2ho(nuv,listno(7))
 230      >        + trshu * ( vap2ho(nuv,listno(8)) + vap2ho(nuv,listno(6)))
 231      >        + unsqu *  vap2ho(nuv,listno(5))
 232 c
 233               vap2ho(nuv,m4) = daux
 234      >        + trsqu * vap2ho(nuv,listno(8))
 235      >        + trshu * ( vap2ho(nuv,listno(5)) + vap2ho(nuv,listno(7)))
 236      >        + unsqu *  vap2ho(nuv,listno(6))
 237 c
 238    31       continue
 239 c
 240 c 3.2. ==> Le quadrangle vient d'etre decoupe en 3 triangles
 241 c          on doit creer les valeurs sur les noeuds au milieu des
 242 c          aretes tracees
 243 c          Remarque : regarder cmcdqu pour les conventions
 244 c
 245           elseif ( typdec.ge.31 .and. typdec.le.34 ) then
 246 c
 247             f1hp = - filqua(lequad)
 248 c
 249 c            Par convention :
 250 c            . la premiere arete du 1-er fils va du noeud ni, milieu de
 251 c              l'arete ai du quadrangle pere, vers le sommet commun
 252 c              aux aretes i+1 et i+2.
 253 c            . la troisieme arete du 1-er fils va du noeud ni, milieu de
 254 c              l'arete ai du quadrangle pere, vers le sommet commun
 255 c              aux aretes i+2 et i+3.
 256 c
 257 c         interpolation :
 258 c         interpolee (ui,i=1,8) = -3/16 (u1+u2+u3+u4)
 259 c                                + 3/4 u5 +3/8 (u6+u8) + 1/4 u7
 260 c         avec u5 pour le noeud le plus proche, u6 et u8 pour ceux qui
 261 c         'encadrent' et u7 pour le noeud 'oppose'
 262 c
 263             if ( typdec.eq.31 ) then
 264               iaux1 = listno(5)
 265               iaux2 = listno(6)
 266               iaux3 = listno(7)
 267               iaux4 = listno(8)
 268             elseif ( typdec.eq.32 ) then
 269               iaux1 = listno(6)
 270               iaux2 = listno(7)
 271               iaux3 = listno(8)
 272               iaux4 = listno(5)
 273             elseif ( typdec.eq.33 ) then
 274               iaux1 = listno(7)
 275               iaux2 = listno(8)
 276               iaux3 = listno(5)
 277               iaux4 = listno(6)
 278             else
 279               iaux1 = listno(8)
 280               iaux2 = listno(5)
 281               iaux3 = listno(6)
 282               iaux4 = listno(7)
 283             endif
 284 c
 285             m1 = np2are(aretri(f1hp,1))
 286             m2 = np2are(aretri(f1hp,3))
 287             profho(m1) = 1
 288             profho(m2) = 1
 289 cgn        write(1,*) 'm1 =', m1, ', m2 =', m2
 290 c
 291             do 32 , nuv = 1, nbfop2
 292 c
 293               daux = - trssz * ( vap2ho(nuv,listno(1))
 294      >                         + vap2ho(nuv,listno(2))
 295      >                         + vap2ho(nuv,listno(3))
 296      >                         + vap2ho(nuv,listno(4)) )
 297      >               + trshu * ( vap2ho(nuv,iaux1) + vap2ho(nuv,iaux3) )
 298 c
 299               vap2ho(nuv,m1) = daux + trsqu * vap2ho(nuv,iaux2)
 300      >                       + unsqu * vap2ho(nuv,iaux4)
 301 c
 302               vap2ho(nuv,m2) = daux + trsqu * vap2ho(nuv,iaux4)
 303      >                       + unsqu * vap2ho(nuv,iaux2)
 304 cgn        write(1,1789) vap2ho(nuv,m1), vap2ho(nuv,m2)
 305 c
 306    32       continue
 307 c
 308 c 3.3. ==> Le quadrangle vient d'etre decoupe en 3 quadrangles
 309 c          on doit creer les valeurs sur les noeuds au milieu des
 310 c          aretes tracees
 311 c          Remarque : regarder cmcdqu pour les conventions
 312 c
 313           elseif ( typdec.ge.41 .and. typdec.le.44 ) then
 314 c
 315             f1hp = filqua(lequad)
 316 c
 317 c         pour les noeuds milieux des aretes entre le noeud central
 318 c         et des milieux d'aretes du quadrangle pere :
 319 c         interpolation = -3/16 (u1+u2+u3+u4)
 320 c                        + 3/4 u5 + 3/8 (u6+u8) + 1/4 u7
 321 c         avec u5 pour le noeud le plus proche, u6 et u8 pour ceux qui
 322 c         'encadrent' et u7 pour le noeud 'oppose'
 323 c
 324 c         pour le noeud milieux de l'arete entre le noeud central
 325 c         et un sommet du quadrangle pere :
 326 c         interpolation = -3/16 (u1+u3) - 1/8 u4
 327 c                        + 9/16 (u5+u6) + 3/16 (u7+u8)
 328 c         avec u1 et u3 pour les sommets qui 'encadrent', u2 pour le
 329 c         sommet le plus proche, u4 pour le sommet le plus loin,
 330 c         u5 et u6 pour les noeuds qui 'encadrent'
 331 c         et u7 et u8 pour les noeuds 'opposes'
 332 c
 333             if ( typdec.eq.41 ) then
 334               jaux1 = listno(2)
 335               jaux3 = listno(4)
 336               jaux4 = listno(1)
 337               iaux1 = listno(5)
 338               iaux2 = listno(6)
 339               iaux3 = listno(7)
 340               iaux4 = listno(8)
 341             elseif ( typdec.eq.42 ) then
 342               jaux1 = listno(3)
 343               jaux3 = listno(1)
 344               jaux4 = listno(2)
 345               iaux1 = listno(6)
 346               iaux2 = listno(7)
 347               iaux3 = listno(8)
 348               iaux4 = listno(5)
 349             elseif ( typdec.eq.43 ) then
 350               jaux1 = listno(4)
 351               jaux3 = listno(2)
 352               jaux4 = listno(3)
 353               iaux1 = listno(7)
 354               iaux2 = listno(8)
 355               iaux3 = listno(5)
 356               iaux4 = listno(6)
 357             else
 358               jaux1 = listno(1)
 359               jaux3 = listno(3)
 360               jaux4 = listno(4)
 361               iaux1 = listno(8)
 362               iaux2 = listno(5)
 363               iaux3 = listno(6)
 364               iaux4 = listno(7)
 365             endif
 366 c
 367             m1 = np2are(arequa(f1hp,4))
 368             m2 = np2are(arequa(f1hp,3))
 369             m3 = np2are(arequa(f1hp+1,3))
 370             profho(m1) = 1
 371             profho(m2) = 1
 372             profho(m3) = 1
 373 cgn        write(1,*) 'm1 =', m1, ', m2 =', m2, ', m3 =', m3
 374 cgn        write(1,*) 'listno =', listno
 375 c
 376             do 33 , nuv = 1, nbfop2
 377 c
 378               daux = - trssz * ( vap2ho(nuv,listno(1))
 379      >                         + vap2ho(nuv,listno(2))
 380      >                         + vap2ho(nuv,listno(3))
 381      >                         + vap2ho(nuv,listno(4)) )
 382 c
 383               vap2ho(nuv,m1) = daux
 384      >               + trsqu * vap2ho(nuv,iaux1)
 385      >               + trshu * ( vap2ho(nuv,iaux2) + vap2ho(nuv,iaux4) )
 386      >               + unsqu * vap2ho(nuv,iaux3)
 387 c
 388               vap2ho(nuv,m2) = daux
 389      >               + trsqu * vap2ho(nuv,iaux2)
 390      >               + trshu * ( vap2ho(nuv,iaux3) + vap2ho(nuv,iaux1) )
 391      >               + unsqu * vap2ho(nuv,iaux4)
 392 c
 393               vap2ho(nuv,m3) =
 394      >               - trssz * ( vap2ho(nuv,jaux1)
 395      >               +           vap2ho(nuv,jaux3) )
 396      >               - unshu * vap2ho(nuv,jaux4)
 397      >               + nessz * ( vap2ho(nuv,iaux3) + vap2ho(nuv,iaux4) )
 398      >               + trssz * ( vap2ho(nuv,iaux1) + vap2ho(nuv,iaux2) )
 399 cgn        write(1,1789) vap2ho(nuv,m1), vap2ho(nuv,m2),  vap2ho(nuv,m3)
 400 c
 401    33       continue
 402 c
 403           endif
 404 c
 405         endif
 406 c
 407    10 continue
 408 c
 409       end