src/tool/AP_Conversion/pcsptd.F

   1       subroutine pcsptd ( etan, etanp1, trhn, trhnp1,
   2      >                    prfcan, prfcap,
   3      >                    hettri, filtri, nbantr, anfitr,
   4      >                    ntreca, ntrsca,
   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                       Triangles d'etat anterieur Deux
  30 c                       -                          -
  31 c ______________________________________________________________________
  32 c .        .     .        .                                            .
  33 c .  nom   . e/s . taille .           description                      .
  34 c .____________________________________________________________________.
  35 c . etan   . e   .    1   . ETAt du triangle a l'iteration N           .
  36 c . etanp1 . e   .    1   . ETAt du triangle a l'iteration N+1         .
  37 c . trhn   . e   .    1   . TRiangle courant en numerotation Homard    .
  38 c .        .     .        . a l'iteration N                            .
  39 c . trhnp1 . e   .    1   . TRiangle 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 . hettri . e   . nbtrto . historique de l'etat des triangles         .
  48 c . filtri . e   . nbtrto . premier fils des triangles                 .
  49 c . nbantr . e   .   1    . nombre de triangles decoupes par           .
  50 c .        .     .        . conformite sur le maillage avant adaptation.
  51 c . anfitr . e   . nbantr . tableau filtri du maillage de l'iteration n.
  52 c . ntreca . e   . retrto . numero des triangles dans le calcul entree .
  53 c . ntrsca . e   . rstrto . numero des triangles du calcul en 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*. tableau temporaire de la solution          .
  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 = 'PCSPTD' )
  81 c
  82 #include "nblang.h"
  83 c
  84 c 0.2. ==> communs
  85 c
  86 #include "nombtr.h"
  87 #include "nombsr.h"
  88 #include "nomber.h"
  89 c
  90 c 0.3. ==> arguments
  91 c
  92       integer etan, etanp1, trhn, trhnp1
  93       integer nbfonc, ngauss
  94       integer prfcan(*), prfcap(*)
  95       integer hettri(nbtrto), filtri(nbtrto)
  96       integer nbantr
  97       integer anfitr(nbantr)
  98       integer ntreca(retrto), ntrsca(rstrto)
  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     trcnp1 = TRiangle courant en numerotation Calcul a l'iteration N+1
 108 c
 109       integer trcnp1
 110 c
 111 c     f1hp = Fils 1er du triangle en numerotation Homard a l'it. N+1
 112 c     f1cp = Fils 1er du triangle en numerota. du Calcul a l'it. N+1
 113 c     f2cp = Fils 2eme du triangle en numerota. du Calcul a l'it. N+1
 114 c
 115       integer f1hp, fihp
 116       integer f1cp, f2cp
 117 c
 118 c     f1hn = Fils 1er du triangle en numerotation Homard a l'it. N
 119 c     f1cn = Fils 1er du triangle en numerotation du Calcul a l'it. N
 120 c     f2cn = Fils 2eme du triangle en numerotation du Calcul a l'it. N
 121 c
 122       integer f1hn
 123       integer f1cn, f2cn
 124 c
 125       integer iaux
 126       integer lglist, nrlist
 127       integer list(30)
 128 c
 129       integer coderr
 130       integer nrofon, nugaus
 131 c
 132       double precision daux
 133       double precision daux1
 134 c
 135       integer nbmess
 136       parameter ( nbmess = 10 )
 137       character*80 texte(nblang,nbmess)
 138 c
 139 c 0.5. ==> initialisations
 140 c ______________________________________________________________________
 141 c
 142 c====
 143 c 1. initialisations
 144 c====
 145 c
 146 #include "impr01.h"
 147 c
 148 #ifdef _DEBUG_HOMARD_
 149       write (ulsort,texte(langue,1)) 'Entree', nompro
 150       call dmflsh (iaux)
 151 #endif
 152 c
 153       texte(1,4) =
 154      >'(/,''Triangle en cours : nro a l''''iteration '',a3,'' : '',i10)'
 155       texte(1,5) =
 156      > '(  ''                etat a l''''iteration '',a3,''   : '',i4)'
 157 c
 158       texte(2,4) =
 159      >'(/,''Current triangle : # at iteration '',a3,''     : '',i10)'
 160       texte(2,5) =
 161      > '(  ''              status at iteration '',a3,'' : '',i4)'
 162 c
 163 #include "impr03.h"
 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 = anfitr(trhn)
 171       f1cn = ntreca(f1hn)
 172       f2cn = ntreca(f1hn+1)
 173 c
 174 c====
 175 c 2. etan = 1, 2, 3 : le triangle etait actif
 176 c    On explore tous les etats du triangle a l'iteration n+1
 177 c====
 178 cgn      write (ulsort,90002) 'etanp1', etanp1
 179 c
 180       if ( prfcan(f1cn).gt.0 .and. prfcan(f2cn).gt.0 ) then
 181 c
 182 c ===> etanp1 = 0 : le triangle est actif
 183 c             Cela veut dire qu'il est reactive.
 184 c             on lui attribue la valeur moyenne sur les deux anciens
 185 c             fils.
 186 c             remarque : cela arrive seulement avec du deraffinement.
 187 c
 188       if ( etanp1.eq.0 ) then
 189 c
 190         trcnp1 = ntrsca(trhnp1)
 191         prfcap(trcnp1) = 1
 192 c
 193         daux1 = 1.d0/dble(2*ngauss)
 194         do 21 , nrofon = 1 , nbfonc
 195           daux = 0.d0
 196           do 211 , nugaus = 1, ngauss
 197             daux = daux + vafoen(nrofon,nugaus,prfcan(f1cn))
 198      >                  + vafoen(nrofon,nugaus,prfcan(f2cn))
 199   211     continue
 200           daux = daux*daux1
 201           do 212 , nugaus = 1 , ngauss
 202             vafott(nrofon,nugaus,trcnp1) = daux
 203   212     continue
 204    21   continue
 205 c
 206 c       etanp1 = etan : le triangle est decoupe en deux
 207 c                       selon le meme decoupage.
 208 c             c'est ce qui se passe quand un decoupage de conformite
 209 c             est supprime au debut des algorithmes d'adaptation,
 210 c             puis reproduit a la creation du maillage car les faces
 211 c             autour n'ont pas change entre les deux iterations.
 212 c             le fils prend la valeur de la fonction sur l'ancien
 213 c             fils qui etait au meme endroit. comme la procedure de
 214 c             numerotation est la meme (voir cmcdtr), le premier fils
 215 c             est toujours le meme, le second egalement. on prendra
 216 c             alors la valeur sur le fils de rang identique a
 217 c             l'iteration n.
 218 c
 219       elseif ( etanp1.eq.etan ) then
 220 c
 221         f1hp = filtri(trhnp1)
 222         f1cp = ntrsca(f1hp)
 223         f2cp = ntrsca(f1hp+1)
 224         prfcap(f1cp) = 1
 225         prfcap(f2cp) = 1
 226         do 22 , nrofon = 1, nbfonc
 227           do 221 , nugaus = 1 , ngauss
 228             vafott(nrofon,nugaus,f1cp) =
 229      >                               vafoen(nrofon,nugaus,prfcan(f1cn))
 230             vafott(nrofon,nugaus,f2cp) =
 231      >                               vafoen(nrofon,nugaus,prfcan(f2cn))
 232   221     continue
 233    22   continue
 234 c
 235 c    etanp1 = 1, 2, 3 et different de etan : le triangle est decoupe
 236 c             en deux mais par un autre decoupage.
 237 c             c'est ce qui se passe quand un decoupage de conformite
 238 c             est supprime au debut des algorithmes d'adaptation. il y
 239 c             a du deraffinement dans la zone qui induisait le decoupage
 240 c             de conformite et raffinement sur une autre zone.
 241 c             on donne la valeur moyenne de la fonction sur les deux
 242 c             anciens fils a chaque nouveau fils.
 243 c             remarque : cela arrive seulement avec du deraffinement.
 244 c
 245       elseif ( etanp1.ge.1 .and. etanp1.le.3 ) then
 246 c
 247         f1hp = filtri(trhnp1)
 248         f1cp = ntrsca(f1hp)
 249         f2cp = ntrsca(f1hp+1)
 250         prfcap(f1cp) = 1
 251         prfcap(f2cp) = 1
 252         daux1 = 1.d0/dble(2*ngauss)
 253         do 23 , nrofon = 1 , nbfonc
 254           daux = 0.d0
 255           do 231 , nugaus = 1, ngauss
 256             daux = daux + vafoen(nrofon,nugaus,prfcan(f1cn))
 257      >                  + vafoen(nrofon,nugaus,prfcan(f2cn))
 258   231     continue
 259           daux = daux*daux1
 260           do 232 , nugaus = 1 , ngauss
 261             vafott(nrofon,nugaus,trcnp1) = daux
 262   232     continue
 263    23   continue
 264 c
 265 c    etanp1 = 4, 6, 7, 8 : le triangle est decoupe en 4
 266 c             on donne la valeur moyenne de la fonction sur les deux
 267 c             anciens fils a chaque nouveau fils.
 268 c
 269       elseif ( etanp1.eq.4 .or.
 270      >         ( etanp1.ge.6 .and. etanp1.le.8 ) ) then
 271 c
 272         f1hp = filtri(trhnp1)
 273         lglist = 0
 274         do 241 , nrlist = 1 , 4
 275           fihp = f1hp+nrlist-1
 276           iaux = mod(hettri(fihp),10)
 277           if ( iaux.eq.0 ) then
 278             lglist = lglist + 1
 279             list(lglist) = ntrsca(fihp)
 280           elseif ( iaux.ge.1 .and. iaux.le.3 ) then
 281             lglist = lglist + 1
 282             list(lglist) = ntrsca(filtri(fihp))
 283             lglist = lglist + 1
 284             list(lglist) = ntrsca(filtri(fihp)+1)
 285           else
 286             coderr = 1
 287           endif
 288   241   continue
 289 c
 290         do 242 , nrlist = 1 , lglist
 291           prfcap(list(nrlist)) = 1
 292   242   continue
 293 c
 294         daux1 = 1.d0/dble(2*ngauss)
 295         do 24 , nrofon = 1 , nbfonc
 296           daux = 0.d0
 297           do 243 , nugaus = 1, ngauss
 298             daux = daux + vafoen(nrofon,nugaus,prfcan(f1cn))
 299      >                  + vafoen(nrofon,nugaus,prfcan(f2cn))
 300   243     continue
 301           daux = daux*daux1
 302           do 244 , nugaus = 1, ngauss
 303             do 245 , nrlist = 1 , lglist
 304               vafott(nrofon,nugaus,list(nrlist)) = daux
 305   245       continue
 306   244     continue
 307    24   continue
 308 c
 309 c
 310 c doc.0.erreur. ==> aucun autre etat sur le  courant elgnairt
 311 c                   n'est possible
 312 c
 313       else
 314 c
 315         coderr = 1
 316         write (ulsort,texte(langue,4)) 'n  ', trhn
 317         write (ulsort,texte(langue,5)) 'n  ', etan
 318         write (ulsort,texte(langue,4)) 'n+1', trhnp1
 319         write (ulsort,texte(langue,5)) 'n+1', etanp1
 320 c
 321       endif
 322 c
 323       endif
 324 c
 325 c====
 326 c 3. la fin
 327 c====
 328 c
 329       if ( coderr.ne.0 ) then
 330 c
 331       write (ulsort,texte(langue,1)) 'Sortie', nompro
 332       write (ulsort,texte(langue,2)) coderr
 333       codret = codret + 1
 334 c
 335       endif
 336 c
 337       end