src/tool/Creation_Maillage/cmch64.F

   1       subroutine cmch64 ( lehexa, listar, listso,
   2      >                    indare, indtri, indpyr,
   3      >                    indptp,
   4      >                    hetare, somare,
   5      >                    filare, merare, famare,
   6      >                    hettri, aretri,
   7      >                    filtri, pertri, famtri,
   8      >                    nivtri,
   9      >                    filqua,
  10      >                    hetpyr, facpyr, cofapy,
  11      >                    filpyr, perpyr, fampyr,
  12      >                    quahex, coquhe,
  13      >                    famhex, cfahex,
  14      >                    ulsort, langue, codret )
  15 c ______________________________________________________________________
  16 c
  17 c                             H O M A R D
  18 c
  19 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
  20 c
  21 c Version originale enregistree le 18 juin 1996 sous le numero 96036
  22 c aupres des huissiers de justice Simart et Lavoir a Clamart
  23 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
  24 c aupres des huissiers de justice
  25 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
  26 c
  27 c    HOMARD est une marque deposee d'Electricite de France
  28 c
  29 c Copyright EDF 1996
  30 c Copyright EDF 1998
  31 c Copyright EDF 2002
  32 c Copyright EDF 2020
  33 c ______________________________________________________________________
  34 c
  35 c    Creation du Maillage - Conformite - decoupage des Hexaedres
  36 c    -           -          -                          -
  37 c                         - par 1 Arete - etat 64
  38 c                                              --
  39 c ______________________________________________________________________
  40 c .        .     .        .                                            .
  41 c .  nom   . e/s . taille .           description                      .
  42 c .____________________________________________________________________.
  43 c . lehexa . e   .   1    . hexaedre a decouper                        .
  44 c . listar . e   .   12   . liste des aretes de l'hexaedre a decouper  .
  45 c . listso . e   .    8   . liste des sommets de l'hexaedre a decouper .
  46 c . indare . es  .   1    . indice de la derniere arete creee          .
  47 c . indtri . es  .   1    . indice du dernier triangle cree            .
  48 c . indpyr . es  .   1    . indice de la derniere pyramide creee       .
  49 c . indptp . e   .   1    . indice du dernier pere enregistre          .
  50 c . hetare . es  . nouvar . historique de l'etat des aretes            .
  51 c . somare . es  .2*nouvar. numeros des extremites d'arete             .
  52 c . filare . es  . nouvar . premiere fille des aretes                  .
  53 c . merare . es  . nouvar . mere des aretes                            .
  54 c . famare .     . nouvar . famille des aretes                         .
  55 c . hettri . es  . nouvtr . historique de l'etat des triangles         .
  56 c . aretri . es  .nouvtr*3. numeros des 3 aretes des triangles         .
  57 c . filtri . es  . nouvtr . premier fils des triangles                 .
  58 c . pertri . es  . nouvtr . pere des triangles                         .
  59 c . famtri . es  . nouvtr . famille des triangles                      .
  60 c . nivtri . es  . nouvtr . niveau des triangles                       .
  61 c . filqua . e   . nouvqu . premier fils des quadrangles               .
  62 c . hetpyr . e   . nouvpy . historique de l'etat des pyramides         .
  63 c . facpyr . e   .nouvyf*5. numeros des 5 faces des pyramides          .
  64 c . cofapy . e   .nouvyf*5. codes des faces des pyramides              .
  65 c . filpyr . e   . nouvpy . premier fils des pyramides                 .
  66 c . perpyr . e   . nouvpy . pere des pyramides                         .
  67 c .        .     .        . si perpyr(i) > 0 : numero de la pyramide   .
  68 c .        .     .        . si perpyr(i) < 0 : -numero dans pphepe     .
  69 c . fampyr . e   . nouvpy . famille des pyramides                      .
  70 c . quahex . e   .nouvhf*6. numeros des 6 quadrangles des hexaedres    .
  71 c . coquhe . e   .nouvhf*6. codes des 6 quadrangles des hexaedres      .
  72 c . famhex . e   . nouvhe . famille des hexaedres                      .
  73 c . cfahex .     . nctfhe. codes des familles des hexaedres            .
  74 c .        .     . nbfhex .   1 : famille MED                          .
  75 c .        .     .        .   2 : type d'hexaedres                     .
  76 c .        .     .        .   3 : famille des tetraedres de conformite .
  77 c .        .     .        .   4 : famille des pyramides de conformite  .
  78 c . ulsort . e   .   1    . unite logique de la sortie generale        .
  79 c . langue . e   .    1   . langue des messages                        .
  80 c .        .     .        . 1 : francais, 2 : anglais                  .
  81 c . codret . es  .    1   . code de retour des modules                 .
  82 c .        .     .        . 0 : pas de probleme                        .
  83 c .        .     .        . 1 : aucune arete ne correspond             .
  84 c ______________________________________________________________________
  85 c
  86 c====
  87 c 0. declarations et dimensionnement
  88 c====
  89 c
  90 c 0.1. ==> generalites
  91 c
  92       implicit none
  93       save
  94 c
  95       character*6 nompro
  96       parameter ( nompro = 'CMCH64' )
  97 c
  98 #include "nblang.h"
  99 c
 100 c 0.2. ==> communs
 101 c
 102 #include "envex1.h"
 103 c
 104 #include "dicfen.h"
 105 #include "nbfami.h"
 106 #include "nouvnb.h"
 107 #include "cofpfh.h"
 108 c
 109 c 0.3. ==> arguments
 110 c
 111       integer lehexa
 112       integer listar(12), listso(8)
 113       integer indare, indtri, indpyr
 114       integer indptp
 115       integer hetare(nouvar), somare(2,nouvar)
 116       integer filare(nouvar), merare(nouvar), famare(nouvar)
 117       integer hettri(nouvtr), aretri(nouvtr,3)
 118       integer filtri(nouvtr), pertri(nouvtr), famtri(nouvtr)
 119       integer nivtri(nouvtr)
 120       integer filqua(nouvqu)
 121       integer hetpyr(nouvpy), facpyr(nouvyf,5), cofapy(nouvyf,5)
 122       integer filpyr(nouvpy), perpyr(nouvpy), fampyr(nouvpy)
 123       integer quahex(nouvhf,6), coquhe(nouvhf,6)
 124       integer famhex(nouvhe), cfahex(nctfhe,nbfhex)
 125 c
 126       integer ulsort, langue, codret
 127 c
 128 c 0.4. ==> variables locales
 129 c
 130       integer iaux, jaux
 131       integer nlarco, nuarco
 132       integer noemil, somm(2)
 133       integer areint(2)
 134       integer areqtr(2,2)
 135       integer triint(5)
 136       integer f1, f2, f3, f4, f5, f6
 137       integer cf1, cf2, cf3, cf4, cf5, cf6
 138       integer trifad(2,0:2), cotrvo(2,0:2)
 139       integer niveau
 140       integer laface, coface
 141 c
 142       integer nbmess
 143       parameter ( nbmess = 10 )
 144       character*80 texte(nblang,nbmess)
 145 c
 146 c 0.5. ==> initialisations
 147 c ______________________________________________________________________
 148 c
 149 c====
 150 c 1. initialisations
 151 c====
 152 c
 153 c 1.1. ==> messages
 154 c
 155 #include "impr01.h"
 156 c
 157 #ifdef _DEBUG_HOMARD_
 158       write (ulsort,texte(langue,1)) 'Entree', nompro
 159       call dmflsh (iaux)
 160 #endif
 161 c
 162 #ifdef _DEBUG_HOMARD_
 163  1789 format(5(a,i5,', '))
 164 #endif
 165 c
 166       codret = 0
 167 c
 168 c 1.2. ==> grandeurs independantes du cas traite (phase 1)
 169 c          les faces de l'hexaedre et leurs codes
 170 c
 171       f1 = quahex(lehexa,1)
 172       f2 = quahex(lehexa,2)
 173       f3 = quahex(lehexa,3)
 174       f4 = quahex(lehexa,4)
 175       f5 = quahex(lehexa,5)
 176       f6 = quahex(lehexa,6)
 177       cf1 = coquhe(lehexa,1)
 178       cf2 = coquhe(lehexa,2)
 179       cf3 = coquhe(lehexa,3)
 180       cf4 = coquhe(lehexa,4)
 181       cf5 = coquhe(lehexa,5)
 182       cf6 = coquhe(lehexa,6)
 183 c
 184 c 1.3. ==> grandeurs dependant du cas traite
 185 c     nlarco = numero local de l'arete coupee
 186       nlarco = 4
 187 c
 188 c     nuarco = numero global de l'arete coupee
 189       nuarco = listar(nlarco)
 190 c
 191 c     noemil = noeud milieu de l'arete coupee
 192       noemil = somare(2,filare(nuarco))
 193 c
 194 c     somm(1) = sommet a joindre au milieu de l'arete coupee pour
 195 c               definir la 1ere arete interne
 196       somm(1) = listso(5)
 197 c     somm(2) = sommet a joindre au milieu de l'arete coupee pour
 198 c               definir la 2nde arete interne
 199       somm(2) = listso(6)
 200 #ifdef _DEBUG_HOMARD_
 201       write(ulsort,2000) 'listso', listso
 202       write(ulsort,2000) 'nuarco', nuarco
 203       write(ulsort,2000) 'noemil', noemil
 204       write(ulsort,2001) 'somm(1)', somm(1),'somm(2)', somm(2)
 205  2000 format(a,10i10)
 206  2001 format(a,i10,', ',a,i10)
 207 #endif
 208 c
 209 c 2.2.6. ==> Triangles et aretes tracees sur les faces coupees en 3
 210 c            L'arete coupee s'appuie sur deux faces de l'hexaedre.
 211 c            trifad(1,*)  se rapporte a celle de plus petit numero local
 212 c            trifad(2,*)  se rapporte a celle de plus grand numero local
 213 c     trifad(p,0) : triangle central de ce decoupage
 214 c     trifad(p,1) : triangle contenant le sommet de l'arete coupee qui a
 215 c                   le plus petit numero local
 216 c     trifad(p,2) : triangle contenant le sommet de l'arete coupee qui a
 217 c                   le plus petit numero local
 218 c     cotrvo(p,0/1/2) : futur code du triangle trifad(p,0/1/2) dans la
 219 c                       description de la pyramide
 220 c     areqtr(p,1) : arete interne au quadrangle de bord et bordant le
 221 c                   triangle trifad(p,1)
 222 c     areqtr(p,2) : arete interne au quadrangle de bord et bordant le
 223 c                   triangle trifad(p,2)
 224 c
 225 c     trifad(1,0) = triangle central de la face 1 : FF1
 226 c     trifad(1,1) = triangle de la face 1 du cote de S2 : FF1 + 1/2
 227 c     trifad(1,2) = triangle de la face 1 du cote de S3 : FF1 + 2/1
 228 c     areqtr(1,1) : AS2N4
 229 c     areqtr(1,2) : AS1N4
 230       laface = f1
 231       coface = cf1
 232       trifad(1,0) = -filqua(laface)
 233       if ( coface.lt.5 ) then
 234         cotrvo(1,0) = 2
 235         trifad(1,1) = trifad(1,0) + 1
 236         cotrvo(1,1) = 2
 237         trifad(1,2) = trifad(1,0) + 2
 238         cotrvo(1,2) = 1
 239         areqtr(1,1) = aretri(trifad(1,0),1)
 240         areqtr(1,2) = aretri(trifad(1,0),3)
 241       else
 242         cotrvo(1,0) = 4
 243         trifad(1,1) = trifad(1,0) + 2
 244         cotrvo(1,1) = 6
 245         trifad(1,2) = trifad(1,0) + 1
 246         cotrvo(1,2) = 4
 247         areqtr(1,1) = aretri(trifad(1,0),3)
 248         areqtr(1,2) = aretri(trifad(1,0),1)
 249       endif
 250 #ifdef _DEBUG_HOMARD_
 251       write(ulsort,1789) 'laface = ', laface,', coface = ', coface
 252       write(ulsort,1789) 'trifad(1,0) = ', trifad(1,0),
 253      >                   'trifad(1,1) = ', trifad(1,1),
 254      >                   'trifad(1,2) = ', trifad(1,2)
 255       write(ulsort,1789) 'cotrvo(1,0) = ', cotrvo(1,0),
 256      >                   'cotrvo(1,1) = ', cotrvo(1,1),
 257      >                   'cotrvo(1,2) = ', cotrvo(1,2)
 258       write(ulsort,1789) 'areqtr(1,1) = ', areqtr(1,1),
 259      >                   ' de ',somare(1,areqtr(1,1)),
 260      >                   ' a ',somare(2,areqtr(1,1))
 261       write(ulsort,1789) 'areqtr(1,2) = ', areqtr(1,2),
 262      >                   ' de ',somare(1,areqtr(1,2)),
 263      >                   ' a ',somare(2,areqtr(1,2))
 264 #endif
 265 c
 266 c     trifad(2,0) = triangle central de la face 2 : FF5
 267 c     trifad(2,1) = triangle de la face 2 du cote de S3 : FF5 + 2/1
 268 c     trifad(2,2) = triangle de la face 2 du cote de S4 : FF5 + 1/2
 269 c     areqtr(2,1) : AS8N4
 270 c     areqtr(2,2) : AS7N4
 271       laface = f5
 272       coface = cf5
 273       trifad(2,0) = -filqua(laface)
 274       if ( coface.lt.5 ) then
 275         cotrvo(2,0) = 2
 276         trifad(2,1) = trifad(2,0) + 2
 277         cotrvo(2,1) = 1
 278         trifad(2,2) = trifad(2,0) + 1
 279         cotrvo(2,2) = 1
 280         areqtr(2,1) = aretri(trifad(2,0),3)
 281         areqtr(2,2) = aretri(trifad(2,0),1)
 282       else
 283         cotrvo(2,0) = 4
 284         trifad(2,1) = trifad(2,0) + 1
 285         cotrvo(2,1) = 4
 286         trifad(2,2) = trifad(2,0) + 2
 287         cotrvo(2,2) = 4
 288         areqtr(2,1) = aretri(trifad(2,0),1)
 289         areqtr(2,2) = aretri(trifad(2,0),3)
 290       endif
 291 #ifdef _DEBUG_HOMARD_
 292       write(ulsort,1789) 'laface = ', laface,', coface = ', coface
 293       write(ulsort,1789) 'trifad(2,0) = ', trifad(2,0),
 294      >                   'trifad(2,1) = ', trifad(2,1),
 295      >                   'trifad(2,2) = ', trifad(2,2)
 296       write(ulsort,1789) 'cotrvo(2,0) = ', cotrvo(2,0),
 297      >                   'cotrvo(2,1) = ', cotrvo(2,1),
 298      >                   'cotrvo(2,2) = ', cotrvo(2,2)
 299       write(ulsort,1789) 'areqtr(2,1) = ', areqtr(2,1),
 300      >                   ' de ',somare(1,areqtr(2,1)),
 301      >                   ' a ',somare(2,areqtr(2,1))
 302       write(ulsort,1789) 'areqtr(2,2) = ', areqtr(2,2),
 303      >                   ' de ',somare(1,areqtr(2,2)),
 304      >                   ' a ',somare(2,areqtr(2,2))
 305 #endif
 306 c
 307 c 1.4. ==> grandeurs independantes du cas traite (phase 2)
 308 c
 309 c     niveau = niveau des triangles des conformites des faces
 310       niveau = nivtri(trifad(1,0))
 311 #ifdef _DEBUG_HOMARD_
 312       write(ulsort,1400) niveau
 313  1400 format('niveau =',i3)
 314 #endif
 315 c
 316 c====
 317 c 2. Creation des deux aretes internes
 318 c    noemil : N4
 319 c    somm(1) : S5
 320 c    somm(2) : S6
 321 c    areint(1) : AS5N4
 322 c    areint(2) : AS6N4
 323 c====
 324 c
 325       if ( codret.eq.0 ) then
 326 c
 327       do 21 , iaux = 1 , 2
 328 c
 329         indare = indare + 1
 330         areint(iaux) = indare
 331 c
 332         somare(1,areint(iaux)) = min ( noemil , somm(iaux) )
 333         somare(2,areint(iaux)) = max ( noemil , somm(iaux) )
 334 c
 335         famare(areint(iaux)) = 1
 336         hetare(areint(iaux)) = 50
 337         merare(areint(iaux)) = 0
 338         filare(areint(iaux)) = 0
 339 #ifdef _DEBUG_HOMARD_
 340       write(ulsort,1789) 'areint(iaux) = ', areint(iaux),
 341      >                   ' de ',somare(1,areint(iaux)),
 342      >                   ' a ',somare(2,areint(iaux))
 343 #endif
 344 c
 345    21 continue
 346 c
 347       endif
 348 c
 349 c====
 350 c 3. Creation des cinq triangles internes
 351 c    areqtr(1,1) : AS2N4
 352 c    areqtr(1,2) : AS1N4
 353 c    areqtr(2,1) : AS8N4
 354 c    areqtr(2,2) : AS7N4
 355 c    areint(1) : AS5N4
 356 c    areint(2) : AS6N4
 357 c    triint(1) : le triangle contenant l'arete areqtr(1,1)
 358 c    triint(3) : le triangle contenant l'arete areqtr(1,2)
 359 c    triint(2) : le triangle contenant l'arete areqtr(2,1)
 360 c    triint(4) : le triangle contenant l'arete areqtr(2,2)
 361 c    triint(5) : le triangle qui s'appuie sur l'arete opposee a l'arete
 362 c                coupee ; il ne touche donc pas les faces coupees
 363 c    triint(1) : FN4A6
 364 c    triint(2) : FN4A5
 365 c    triint(3) : FN4A11
 366 c    triint(4) : FN4A10
 367 c    triint(5) : FN4A9
 368 c     par convention, le niveau est le meme que les triangles fils
 369 c     sur l'exterieur
 370 c====
 371 c
 372       iaux = 1
 373 c
 374 #ifdef _DEBUG_HOMARD_
 375       write (ulsort,texte(langue,3)) 'CMCTRI_64', nompro
 376       write (ulsort,3000) indtri+1, indtri+5
 377  3000 format('.. triangles de',i10,' a',i10)
 378 #endif
 379       triint(1) = indtri + 1
 380       call cmctri ( aretri, famtri, hettri,
 381      >              filtri, pertri, nivtri,
 382      >              triint(1), listar( 6), areqtr(1,1), areint(1),
 383      >              iaux, niveau )
 384 c
 385       triint(2) = indtri + 2
 386       call cmctri ( aretri, famtri, hettri,
 387      >              filtri, pertri, nivtri,
 388      >              triint(2), listar( 5), areqtr(1,2), areint(2),
 389      >              iaux, niveau )
 390 c
 391       triint(3) = indtri + 3
 392       call cmctri ( aretri, famtri, hettri,
 393      >              filtri, pertri, nivtri,
 394      >              triint(3), listar(11), areint(1), areqtr(2,1),
 395      >              iaux, niveau )
 396 c
 397       triint(4) = indtri + 4
 398       call cmctri ( aretri, famtri, hettri,
 399      >              filtri, pertri, nivtri,
 400      >              triint(4), listar(10), areint(2), areqtr(2,2),
 401      >              iaux, niveau )
 402 c
 403       triint(5) = indtri + 5
 404       call cmctri ( aretri, famtri, hettri,
 405      >              filtri, pertri, nivtri,
 406      >              triint(5), listar( 9), areint(1), areint(2),
 407      >              iaux, niveau )
 408 c
 409       indtri = triint(5)
 410 c
 411 c====
 412 c 4. Creation des quatre pyramides
 413 c    Elles arrivent dans l'ordre de numerotation locale de leur
 414 c    quadrangle dans l'hexaedre
 415 c     trifad(1,0) : FF1
 416 c     trifad(1,1) : FF1 + 1/2
 417 c     trifad(1,2) : FF1 + 2/1
 418 c     trifad(2,0) : FF5
 419 c     trifad(2,1) : FF5 + 2/1
 420 c     trifad(2,2) : FF5 + 1/2
 421 c    triint(1) : FN4A6
 422 c    triint(2) : FN4A5
 423 c    triint(3) : FN4A11
 424 c    triint(4) : FN4A10
 425 c    triint(5) : FN4A9
 426 c====
 427 c
 428       jaux = cfahex(cofpfh,famhex(lehexa))
 429       iaux = -indptp
 430 c
 431 #ifdef _DEBUG_HOMARD_
 432       write (ulsort,texte(langue,3)) 'CMCPYR_64', nompro
 433       write (ulsort,4000) indpyr+1, indpyr+4
 434  4000 format('.. pyramides de',i10,' a',i10)
 435 #endif
 436       indpyr = indpyr + 1
 437       call cmcpyr ( facpyr, cofapy, fampyr, hetpyr, filpyr, perpyr,
 438      >              trifad(1,0), cotrvo(1,0),
 439      >                triint(1),           3,
 440      >                triint(5),           3,
 441      >                triint(2),           6,
 442      >                       f2,         cf2,
 443      >              iaux,  jaux,   indpyr )
 444 c
 445       indpyr = indpyr + 1
 446       call cmcpyr ( facpyr, cofapy, fampyr, hetpyr, filpyr, perpyr,
 447      >              trifad(1,2), cotrvo(1,2),
 448      >                triint(2),           3,
 449      >                triint(4),           3,
 450      >              trifad(2,2), cotrvo(2,2),
 451      >                       f3,         cf3,
 452      >              iaux,  jaux,   indpyr )
 453 c
 454       indpyr = indpyr + 1
 455       call cmcpyr ( facpyr, cofapy, fampyr, hetpyr, filpyr, perpyr,
 456      >              trifad(1,1), cotrvo(1,1),
 457      >              trifad(2,1), cotrvo(2,1),
 458      >                triint(3),           5,
 459      >                triint(1),           6,
 460      >                       f4,         cf4,
 461      >              iaux,  jaux,   indpyr )
 462 c
 463       indpyr = indpyr + 1
 464       call cmcpyr ( facpyr, cofapy, fampyr, hetpyr, filpyr, perpyr,
 465      >                triint(5),           5,
 466      >                triint(3),           3,
 467      >              trifad(2,0), cotrvo(2,0),
 468      >                triint(4),           6,
 469      >                       f6,         cf6,
 470      >              iaux,  jaux,   indpyr )
 471 c
 472 c====
 473 c 5. la fin
 474 c====
 475 c
 476       if ( codret.ne.0 ) then
 477 c
 478 #include "envex2.h"
 479 c
 480       write (ulsort,texte(langue,1)) 'Sortie', nompro
 481       write (ulsort,texte(langue,2)) codret
 482       write (ulsort,texte(langue,4))
 483 c
 484       endif
 485 c
 486 #ifdef _DEBUG_HOMARD_
 487       write (ulsort,texte(langue,1)) 'Sortie', nompro
 488       call dmflsh (iaux)
 489 #endif
 490 c
 491       end