src/HEXABLOCK/HexHexa.cxx

   1
   2 // C++ : Gestion des Hexaedres
   3
   4 // Copyright (C) 2009-2013  CEA/DEN, EDF R&D
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 #include "HexHexa.hxx"
  23 #include "HexQuad.hxx"
  24
  25 #include "HexVertex.hxx"
  26 #include "HexDocument.hxx"
  27 #include "HexEdge.hxx"
  28
  29 #include "HexGlobale.hxx"
  30 #include "HexXmlWriter.hxx"
  31 #include "HexDiagnostics.hxx"
  32 #include "HexGlobale.hxx"
  33 #include "HexMatrix.hxx"
  34 #include "HexElements.hxx"
  35
  36 BEGIN_NAMESPACE_HEXA
  37
  38 // ======================================================== Constructeur
  39 Hexa::Hexa (Vertex* va, Vertex* vb, Vertex* vc, Vertex* vd,
  40             Vertex* ve, Vertex* vf, Vertex* vg, Vertex* vh)
  41     : EltBase (va->dad(), EL_HEXA)
  42 {
  43    h_vertex [V_ACE] = va;
  44    h_vertex [V_ACF] = vb;   // = vc ; Modif Abu 30/08/2011
  45    h_vertex [V_ADE] = vc;   // = vb ; Modif Abu 30/08/2011
  46    h_vertex [V_ADF] = vd;
  47
  48    h_vertex [V_BCE] = ve;
  49    h_vertex [V_BCF] = vf;   // = vg ; Modif Abu 30/08/2011
  50    h_vertex [V_BDE] = vg;   // = vf ; Modif Abu 30/08/2011
  51    h_vertex [V_BDF] = vh;
  52
  53    h_clone          = NULL;
  54
  55    controlerSommets ();
  56
  57    Globale* glob = Globale::getInstance ();
  58
  59    for (int nro=0 ; nro<HE_MAXI ; nro++)
  60        {
  61        h_edge[nro] = new Edge (h_vertex[glob->EdgeVertex (nro, V_AMONT)],
  62                                h_vertex[glob->EdgeVertex (nro, V_AVAL)]);
  63        }
  64
  65    for (int nro=0 ; nro<HQ_MAXI ; nro++)
  66        {
  67        h_quad[nro] = new Quad (h_edge [glob->QuadEdge (nro, E_A)],
  68                                h_edge [glob->QuadEdge (nro, E_B)],
  69                                h_edge [glob->QuadEdge (nro, E_C)],
  70                                h_edge [glob->QuadEdge (nro, E_D)]);
  71
  72        h_quad[nro] -> addParent (this);
  73        }
  74    majReferences ();
  75    if (el_root != NULL && el_status==HOK)
  76        el_root->addHexa (this);
  77 }
  78 // ======================================================== Constructeur 2
  79 Hexa::Hexa (Quad* qa, Quad* qb, Quad* qc, Quad* qd, Quad* qe, Quad* qf)
  80     : EltBase (qa->dad(), EL_HEXA)
  81 {
  82    h_quad [Q_A] = qa;
  83    h_quad [Q_B] = qb;
  84    h_quad [Q_C] = qc;
  85    h_quad [Q_D] = qd;
  86    h_quad [Q_E] = qe;
  87    h_quad [Q_F] = qf;
  88    h_clone      = NULL;
  89
  90    for (int nb=0 ; nb<HE_MAXI ; nb++) h_edge   [nb] = NULL;
  91    for (int nb=0 ; nb<HV_MAXI ; nb++) h_vertex [nb] = NULL;
  92
  93    controlerFaces  ();
  94
  95    verifierAretes  ();
  96    verifierSommets ();
  97    majReferences   ();
  98
  99    if (el_status != HOK)
 100       {
 101       printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
 102       printf (" +++ Heaedre impossible \n");
 103       printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
 104       // el_root->dump ();
 105       // printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
 106       dumpFull ();
 107       for (int nro=0; nro<HQ_MAXI; nro++) HexDump (h_quad[nro]);
 108       for (int nro=0; nro<HE_MAXI; nro++) HexDump (h_edge[nro]);
 109       for (int nro=0; nro<HV_MAXI; nro++) HexDump (h_vertex[nro]);
 110       printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
 111       // exit (1);
 112       }
 113    else
 114        el_root->addHexa (this);
 115 }
 116 // =========================================================  Constructeur 3
 117 // === a utiliser uniquement pour le clonage
 118 Hexa::Hexa (Hexa* other)
 119     : EltBase (other->el_root, EL_HEXA)
 120 {
 121    for (int nro=0 ; nro<HQ_MAXI ; nro++) h_quad   [nro] = NULL;
 122    for (int nro=0 ; nro<HE_MAXI ; nro++) h_edge   [nro] = NULL;
 123    for (int nro=0 ; nro<HV_MAXI ; nro++) h_vertex [nro] = NULL;
 124    h_clone = NULL;
 125
 126    if (el_root != NULL && el_status==HOK)
 127        el_root->addHexa (this);
 128 }
 129 // ======================================================== majReferences
 130 void Hexa::majReferences  ()
 131 {
 132    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 133        h_quad[nro]->addParent (this);
 134 }
 135 // ======================================================== controlerArete
 136 void Hexa::controlerArete  (int arete, int face1, int face2)
 137 {
 138    h_edge [arete] = h_quad[face1]->commonEdge (h_quad[face2]);
 139    if (h_edge [arete]==NULL)
 140        {
 141        el_root->putError (W_H_BAD_EDGE,
 142                           el_root->glob->namofHexaEdge (arete));
 143        el_status = 888;
 144        }
 145 }
 146 // ======================================================== controlerSomet
 147 void Hexa::controlerSommet  (int node, int ne1, int ne2, int ne3)
 148 {
 149     if (h_edge[ne1] == NULL || h_edge[ne2] == NULL || h_edge[ne3] == NULL)
 150        return;
 151
 152     Vertex* hv = h_edge[ne1]->commonVertex (h_edge[ne2]);
 153     h_vertex [node] = hv;
 154     if (hv == NULL)
 155        {
 156        el_root->putError (W_H_BAD_VERTEX,
 157                           el_root->glob->namofHexaVertex(node));
 158        el_status = 888;
 159        }
 160     else if (hv == NULL || h_edge[ne3]->index (hv)<0)
 161        {
 162        char b[10];
 163        el_root->putError (W_H_BAD_VERTEX,
 164                           el_root->glob->namofHexaVertex(node));
 165        el_status = 888;
 166        printf ("%s = %s\n", el_root->glob->namofHexaVertex(node),
 167                             hv->getName(b));
 168        printf ("%s = %s\n", el_root->glob->namofHexaEdge(ne1),
 169                             h_edge[ne1]->getName(b));
 170        printf ("%s = %s\n", el_root->glob->namofHexaEdge(ne2),
 171                             h_edge[ne2]->getName(b));
 172        printf ("%s = %s\n", el_root->glob->namofHexaEdge(ne3),
 173                             h_edge[ne3]->getName(b));
 174        }
 175 }
 176 // ======================================================== controlerFaces
 177 void Hexa::controlerFaces  ()
 178 {
 179    for (int n1=0 ; n1<HQ_MAXI ; n1++)
 180        {
 181        if (BadElement (h_quad [n1]))
 182           {
 183           el_root->putError (W_H_NULL_QUAD,
 184                              el_root->glob->namofHexaQuad (n1));
 185           setError (886);
 186           return;
 187           }
 188        for (int n2=n1+1 ; n2<HQ_MAXI ; n2++)
 189            if (h_quad [n1] == h_quad[n2])
 190               {
 191               el_root->putError (W_H_EQ_QUAD,
 192                          el_root->glob->namofHexaQuad (n1),
 193                          el_root->glob->namofHexaQuad (n2));
 194               setError (888);
 195               }
 196        }
 197 }
 198 // ======================================================== controlerSommets
 199 void Hexa::controlerSommets  ()
 200 {
 201    for (int n1=0 ; n1<HV_MAXI ; n1++)
 202        {
 203        if (BadElement (h_vertex [n1]))
 204           {
 205           el_root->putError (W_H_NULL_QUAD,
 206                              el_root->glob->namofHexaVertex (n1));
 207           setError (886);
 208           return;
 209           }
 210        for (int n2=n1+1 ; n2<HQ_MAXI ; n2++)
 211            if (h_vertex [n1] == h_vertex[n2])
 212               {
 213               el_root->putError (W_H_EQ_QUAD,
 214                          el_root->glob->namofHexaVertex (n1),
 215                          el_root->glob->namofHexaVertex (n2));
 216               setError (888);
 217               }
 218        }
 219 }
 220 // ======================================================== verifierAretes
 221 void Hexa::verifierAretes  ()
 222 {
 223    for (int nro=0 ; nro<HE_MAXI ; nro++) h_edge [nro] = NULL;
 224
 225    controlerArete ( E_AC, Q_A, Q_C);
 226    controlerArete ( E_AF, Q_A, Q_F);
 227    controlerArete ( E_AD, Q_A, Q_D);
 228    controlerArete ( E_AE, Q_A, Q_E);
 229
 230    controlerArete ( E_BC, Q_B, Q_C);
 231    controlerArete ( E_BF, Q_B, Q_F);
 232    controlerArete ( E_BD, Q_B, Q_D);
 233    controlerArete ( E_BE, Q_B, Q_E);
 234
 235    controlerArete ( E_CE, Q_C, Q_E);
 236    controlerArete ( E_CF, Q_C, Q_F);
 237    controlerArete ( E_DF, Q_D, Q_F);
 238    controlerArete ( E_DE, Q_D, Q_E);
 239 }
 240 // ======================================================== verifierSommets
 241 void Hexa::verifierSommets  ()
 242 {
 243    for (int nro=0 ; nro<HV_MAXI ; nro++) h_vertex [nro] = NULL;
 244
 245    controlerSommet (V_ACE, E_AC, E_AE, E_CE);
 246    controlerSommet (V_ACF, E_AC, E_AF, E_CF);
 247    controlerSommet (V_ADF, E_AD, E_AF, E_DF);
 248    controlerSommet (V_ADE, E_AD, E_AE, E_DE);
 249
 250    controlerSommet (V_BCE, E_BC, E_BE, E_CE);
 251    controlerSommet (V_BCF, E_BC, E_BF, E_CF);
 252    controlerSommet (V_BDF, E_BD, E_BF, E_DF);
 253    controlerSommet (V_BDE, E_BD, E_BE, E_DE);
 254 }
 255 // ======================================================== Inter
 256 Quad* Hexa::Inter (Hexa* other)
 257 {
 258    for (int nf1=0 ; nf1 < HQ_MAXI ; nf1++)
 259        for (int nf2=0 ; nf2 < HQ_MAXI ; nf2++)
 260            if (h_quad[nf1] == other->h_quad[nf2])
 261               return h_quad[nf1];
 262
 263    return NULL;
 264 }
 265 // -------------------------------------------------------------------
 266 //                         Debug
 267 // -------------------------------------------------------------------
 268 // ======================================================= razNodes
 269 void Hexa::razNodes ()
 270 {
 271    for (int nb=0 ; nb<HV_MAXI ; nb++)
 272        h_vertex[nb]->setMark (NO_COUNTED);
 273 }
 274 // ======================================================= countNodes
 275 int Hexa::countNodes ()
 276 {
 277    int nombre = 0;
 278    for (int nb=0 ; nb<HV_MAXI ; nb++)
 279        if (h_vertex[nb]->getMark () == NO_COUNTED)
 280           {
 281           h_vertex[nb]->setMark (NO_USED);
 282           nombre ++;
 283           }
 284
 285    return nombre;
 286 }
 287 // ======================================================= printNodes
 288 void Hexa::printNodes (pfile vtk, int& compteur)
 289 {
 290    const double minvtk = 1e-30;
 291    Real3 koord;
 292
 293    for (int nb=0 ; nb<HV_MAXI ; nb++)
 294        if (h_vertex[nb]->getMark()==NO_USED)
 295           {
 296           h_vertex[nb]->getPoint (koord);
 297           for (int nc=dir_x ; nc<=dir_z ; nc++)
 298               if (koord [nc] < minvtk &&  koord [nc] > -minvtk)
 299                   koord[nc] = 0.0;
 300
 301           fprintf (vtk, "%g %g %g\n", koord[dir_x], koord[dir_y], koord[dir_z]);
 302           h_vertex[nb]->setMark (compteur);
 303           compteur ++;
 304           }
 305 }
 306 // ======================================================= colorNodes
 307 void Hexa::colorNodes (pfile vtk)
 308 {
 309    for (int nb=0 ; nb<HV_MAXI ; nb++)
 310        if (h_vertex[nb]->getMark()>=0)
 311           {
 312           double color = 100*(h_vertex[nb]->getScalar()+1);
 313           fprintf (vtk, "%g\n", color);
 314           h_vertex[nb]->setMark (NO_COUNTED);
 315           }
 316 }
 317 // ======================================================= moveNodes
 318 void Hexa::moveNodes (Matrix* matrice)
 319 {
 320    for (int nb=0 ; nb<HV_MAXI ; nb++)
 321        if (h_vertex[nb]->getMark()!=IS_USED)
 322           {
 323           matrice->perform (h_vertex[nb]);
 324           h_vertex[nb]->setMark (IS_USED);
 325           }
 326 }
 327 // ======================================================= transform
 328 void Hexa::transform (Matrix* matrice)
 329 {
 330    for (int nb=0 ; nb<HV_MAXI ; nb++)
 331        matrice->perform (h_vertex[nb]);
 332 }
 333 // ======================================================= printHexa
 334 void Hexa::printHexa  (pfile vtk)
 335 {
 336    fprintf (vtk, "%d", HV_MAXI);
 337
 338    fprintf (vtk, " %d", h_vertex[V_ACE]->getMark ());
 339    fprintf (vtk, " %d", h_vertex[V_ACF]->getMark ());
 340    fprintf (vtk, " %d", h_vertex[V_ADF]->getMark ());
 341    fprintf (vtk, " %d", h_vertex[V_ADE]->getMark ());
 342
 343    fprintf (vtk, " %d", h_vertex[V_BCE]->getMark ());
 344    fprintf (vtk, " %d", h_vertex[V_BCF]->getMark ());
 345    fprintf (vtk, " %d", h_vertex[V_BDF]->getMark ());
 346    fprintf (vtk, " %d", h_vertex[V_BDE]->getMark ());
 347
 348    fprintf (vtk, "\n");
 349 }
 350 // ======================================================= printHexaVtk
 351 void Hexa::printHexaVtk (pfile vtk)
 352 {
 353    fprintf (vtk, "%d", HV_MAXI);
 354
 355    fprintf (vtk, " %d", h_vertex[V_ACE]->getId ());
 356    fprintf (vtk, " %d", h_vertex[V_ACF]->getId ());
 357    fprintf (vtk, " %d", h_vertex[V_ADF]->getId ());
 358    fprintf (vtk, " %d", h_vertex[V_ADE]->getId ());
 359
 360    fprintf (vtk, " %d", h_vertex[V_BCE]->getId ());
 361    fprintf (vtk, " %d", h_vertex[V_BCF]->getId ());
 362    fprintf (vtk, " %d", h_vertex[V_BDF]->getId ());
 363    fprintf (vtk, " %d", h_vertex[V_BDE]->getId ());
 364
 365    fprintf (vtk, "\n");
 366 }
 367 // ======================================================== hasFreEdges
 368 bool Hexa::hasFreeEdges ()
 369 {
 370    if (isDeleted())
 371        return false;
 372
 373    for (int nro=0; nro<HE_MAXI ; nro++)
 374        if (h_edge[nro]->getMark()<0)
 375           return true;
 376
 377    return false;
 378 }
 379 // ======================================================== propager
 380 void Hexa::propager (Propagation* prop, int nro)
 381 {
 382    if (isDeleted())
 383        return;
 384
 385    for (int nro=0; nro<HE_MAXI ; nro++)
 386        if (h_edge[nro]->getMark()<0)
 387            h_edge[nro]->propager (prop, nro);
 388 }
 389 // ========================================================= saveXml
 390 void Hexa::saveXml (XmlWriter* xml)
 391 {
 392    char ident[12];
 393    string quads;
 394
 395    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 396        {
 397        if (nro>0) quads += " ";
 398        quads += h_quad[nro]->getName(ident);
 399        }
 400
 401    getName (ident);
 402    xml->openMark     ("Hexa");
 403    xml->addAttribute ("id",    ident);
 404    xml->addAttribute ("quads", quads);
 405    if (el_name!=ident)
 406        xml->addAttribute ("name", el_name);
 407    xml->closeMark ();
 408 }
 409 // ========================================================= findQuad
 410 int Hexa::findQuad (Quad* element)
 411 {
 412    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 413        if (h_quad[nro]==element)
 414           return nro;
 415
 416    return NOTHING;
 417 }
 418 // ========================================================= findEdge
 419 int Hexa::findEdge (Edge* element)
 420 {
 421    for (int nro=0 ; nro<HE_MAXI ; nro++)
 422        if (h_edge[nro]==element)
 423           return nro;
 424
 425    return NOTHING;
 426 }
 427 // ========================================================= findVertex
 428 int Hexa::findVertex (Vertex* element)
 429 {
 430    for (int nro=0 ; nro<HV_MAXI ; nro++)
 431        if (h_vertex[nro]==element)
 432           return nro;
 433
 434    return NOTHING;
 435 }
 436 // ========================================================= disconnectQuad
 437 Elements* Hexa::disconnectQuad (Quad* quad)
 438 {
 439    if (quad==NULL)
 440       return NULL;
 441
 442    if (debug())
 443       {
 444       printf (" ... Avant disconnectQuad, quad=");
 445       quad->printName ("\n");
 446       dumpFull ();
 447       }
 448
 449    int nface = findQuad (quad);
 450    if (nface==NOTHING)
 451       return NULL;
 452                                        // Face opposee : replace
 453    // int nfopp = (nface MODULO 2==0) ? nface+1 : nface-1;
 454
 455    int  ind_edge  [QUAD4], ind_opp_quad [QUAD4];
 456    bool make_quad [QUAD4], make_edge [QUAD4];
 457
 458    for (int nro=0 ; nro<QUAD4 ; nro++)
 459        make_quad [nro] = make_edge[nro] = false;
 460
 461    for (int nro=0 ; nro<QUAD4 ; nro++)
 462        {
 463        int nro1  = (nro+1) MODULO QUAD4;
 464        int pedge = findEdge   (quad->getEdge   (nro));
 465        int pnode = findVertex (quad->getVertex (nro));
 466        int oppq  = findOpposedQuad (quad, quad->getEdge (nro));
 467
 468        ind_edge [nro]     = pedge;
 469        ind_opp_quad [nro] = oppq;
 470
 471        if (pedge==NOTHING || pnode==NOTHING || oppq==NOTHING)
 472           return NULL;
 473
 474        make_quad [nro]  = h_quad[oppq]->getNbrParents() == 2;
 475        make_edge [nro ] = make_edge [nro ] || make_quad [nro];
 476        make_edge [nro1] = make_edge [nro1] || make_quad [nro];
 477
 478        if (debug())
 479           {
 480           printf (" Sommet nro %d : ", nro);
 481           quad->getVertex(nro)->printName (", ");
 482           printf (" edge = ");
 483           quad->getEdge(nro)->printName (", ");
 484           printf (" quad oppose = ");
 485           h_quad[oppq]->printName("");
 486           if (make_quad [nro])
 487              printf (" a dissocier\n");
 488           else
 489              printf ("\n");
 490
 491           }
 492        }
 493
 494    Vertex* new_node     [QUAD4];
 495    Edge*   new_opp_edge [QUAD4];
 496    Edge*   old_opp_edge [QUAD4];
 497
 498    for (int nro=0 ; nro<QUAD4 ; nro++)
 499        {
 500        old_opp_edge [nro] = NULL;
 501        new_opp_edge [nro] = NULL;
 502        Vertex*  o_v0  = quad->getVertex (nro);
 503        new_node [nro] = new Vertex (o_v0);
 504        if (debug())
 505           {
 506           printf (" quad.vertex [%d] = ", nro);
 507           quad->getVertex (nro)->printName (" --> ");
 508           new_node [nro]->printName ("\n");
 509           }
 510
 511        if (make_edge[nro])
 512           {
 513           Quad*   pface  = h_quad [ind_opp_quad [nro]];
 514           int     bid;
 515           int     ncut = pface->inter (quad, bid);
 516           Edge*   ecut = pface->getEdge ((ncut+1) MODULO QUAD4);
 517           Vertex* vopp = ecut->getVertex(V_AMONT);
 518           if (vopp==o_v0)
 519               vopp = ecut->getVertex (V_AVAL);
 520           else if (o_v0 != ecut->getVertex (V_AVAL));
 521               {
 522               ecut = pface->getEdge ((ncut+3) MODULO QUAD4);
 523               vopp = ecut->getVertex(V_AMONT);
 524               if (vopp==o_v0)
 525                   vopp = ecut->getVertex (V_AVAL);
 526               else if (o_v0 != ecut->getVertex (V_AVAL))
 527                   return NULL;
 528               }
 529
 530           old_opp_edge [nro] = ecut;
 531           new_opp_edge [nro] = new Edge (new_node[nro], vopp);
 532           if (debug())
 533              {
 534              printf (" quad.opp_edge [%d] = ", nro);
 535              old_opp_edge [nro]->printName (" --> ");
 536              new_opp_edge [nro]->printName ("\n");
 537              }
 538           }
 539        }
 540
 541    Quad* new_quad = new Quad (new_node[0], new_node[1], new_node[2],
 542                                                         new_node[3]);
 543
 544    Quad* new_opp_quad [QUAD4];
 545    Quad* old_opp_quad [QUAD4];
 546    for (int nro=0 ; nro<QUAD4 ; nro++)
 547        {
 548        old_opp_quad [nro] = NULL;
 549        new_opp_quad [nro] = NULL;
 550        if (make_quad[nro])
 551           {
 552           int nro1 = (nro+1) MODULO QUAD4;
 553
 554           Edge* n_edge0 = new_quad->getEdge (nro);
 555           Edge* n_edge1 = new_opp_edge [nro];
 556           Edge* n_edge3 = new_opp_edge [nro1];
 557
 558           int iv1 = n_edge1->inter (n_edge0);
 559           int iv3 = n_edge3->inter (n_edge0);
 560           if (iv1 <0 || iv3 <0)
 561              return NULL;
 562
 563           Quad* o_face = h_quad [ind_opp_quad [nro]];
 564           Edge* edge2  = o_face->findEdge (n_edge1->getVertex (1-iv1),
 565                                            n_edge3->getVertex (1-iv3));
 566           if (edge2==NULL)
 567              return NULL;
 568           // Edge* o_edge0 = h_edge [ind_edge     [nro]];
 569           // int sens = 1;
 570           // Edge* edge2   = o_face->getOpposEdge (o_edge0, sens);
 571
 572           old_opp_quad [nro] = o_face;
 573           if (debug())
 574              printf (" -------- Quad oppose nro %d\n", nro);
 575           new_opp_quad [nro] = new Quad (n_edge0, n_edge1, edge2, n_edge3);
 576           }
 577        }
 578
 579    replaceQuad (quad, new_quad);
 580    for (int nro=0 ; nro<QUAD4 ; nro++)
 581        if (make_quad[nro])
 582           replaceQuad (old_opp_quad [nro], new_opp_quad [nro]);
 583
 584    for (int nro=0 ; nro<QUAD4 ; nro++)
 585        {
 586        replaceEdge (h_edge[ind_edge[nro]], new_quad->getEdge (nro));
 587        if (make_edge[nro])
 588           replaceEdge (old_opp_edge [nro], new_opp_edge [nro]);
 589        }
 590
 591    for (int nro=0 ; nro<QUAD4 ; nro++)
 592        {
 593        replaceVertex (quad->getVertex(nro), new_node[nro]);
 594        }
 595
 596
 597    h_quad [nface] = new_quad;
 598    if (debug())
 599       {
 600       printf (" ... Apres disconnectQuad, new_quad=");
 601       new_quad->printName ("\n");
 602       dumpFull ();
 603       }
 604
 605    Elements* nouveaux  = new Elements (el_root);
 606    nouveaux->addQuad   (new_quad);
 607    for (int nro=0 ; nro<QUAD4 ; nro++)
 608        {
 609        nouveaux->addEdge   (new_quad->getEdge   (nro));
 610        nouveaux->addVertex (new_quad->getVertex (nro));
 611        if (make_edge[nro])
 612           nouveaux->addEdge (new_opp_edge [nro]);
 613        if (make_quad[nro])
 614           nouveaux->addQuad (new_opp_quad [nro]);
 615        }
 616    nouveaux->moveDisco (this);
 617    return nouveaux;
 618 }
 619 // ========================================================= disconnectEdge
 620 Elements* Hexa::disconnectEdge (Edge* arete)
 621 {
 622    int nedge  = findEdge   (arete);
 623    int namont = findVertex (arete->getVertex(V_AMONT));
 624    int naval  = findVertex (arete->getVertex(V_AVAL ));
 625
 626    if (nedge==NOTHING || namont==NOTHING || naval==NOTHING)
 627       return NULL;
 628
 629    if (debug())
 630       {
 631       printf (" ... Avant disconnectEdge, arete=");
 632       arete->printName ("\n");
 633       dumpFull ();
 634       }
 635
 636    Edge*   n_edge   [HE_MAXI];
 637    Quad*   n_quad   [HQ_MAXI];
 638    Vertex* n_vertex [HV_MAXI];
 639
 640    for (int nro=0 ; nro<HQ_MAXI ; nro++) n_quad [nro]   = NULL;
 641    for (int nro=0 ; nro<HE_MAXI ; nro++) n_edge [nro]   = NULL;
 642    for (int nro=0 ; nro<HV_MAXI ; nro++) n_vertex [nro] = NULL;
 643
 644    Vertex* old_amont = arete->getVertex (V_AMONT);
 645    Vertex* old_aval  = arete->getVertex (V_AVAL );
 646    Vertex* new_amont = n_vertex [namont] = new Vertex (old_amont);
 647    Vertex* new_aval  = n_vertex [naval]  = new Vertex (old_aval);
 648    n_edge [nedge]    = new Edge   (new_amont, new_aval);
 649
 650    // Un edge non remplace, qui contient un vertex remplace
 651    //         commun a plus de 2 faces (donc appartenant a un autre hexa)
 652    //         doit etre duplique
 653
 654    for (int nro=0 ; nro<HE_MAXI ; nro++)
 655        {
 656        if (   n_edge[nro]==NULL && h_edge[nro] != NULL
 657            && h_edge[nro]->getNbrParents()>2)
 658           {
 659           Vertex* va =  h_edge[nro]->getVertex (V_AMONT);
 660           Vertex* vb =  h_edge[nro]->getVertex (V_AVAL);
 661
 662           if (va==old_amont)
 663              n_edge [nro] = new Edge (new_amont, vb);
 664           else if (va==old_aval)
 665              n_edge [nro] = new Edge (new_aval,  vb);
 666           else if (vb==old_amont)
 667              n_edge [nro] = new Edge (va, new_amont);
 668           else if (vb==old_aval)
 669              n_edge [nro] = new Edge (va, new_aval);
 670           }
 671        }
 672
 673    // Un quad non remplace, qui contient un edge remplace
 674    //         commun a plus de 2 Hexas
 675    //         doit etre duplique
 676
 677    Globale* glob = Globale::getInstance();
 678    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 679        if (   n_quad[nro]==NULL && h_quad[nro] != NULL
 680            && h_quad[nro]->getNbrParents()>1)
 681           {
 682           Edge* qedge[QUAD4];
 683           bool  duplic = false;
 684           for (int ned=0 ; ned<QUAD4 ; ned++)
 685               {
 686               int ndup = glob->QuadEdge (nro, (EnumQuad)ned);
 687               if (n_edge [ndup] ==NULL)
 688                  qedge [ned] = h_edge[ndup];
 689               else
 690                  {
 691                  qedge [ned] = n_edge[ndup];
 692                  duplic = true;
 693                  }
 694               }
 695           if (duplic)
 696              n_quad[nro] = new Quad (qedge[Q_A], qedge[Q_B],
 697                                      qedge[Q_C], qedge[Q_D]);
 698           }
 699
 700    Elements* nouveaux  = new Elements (el_root);
 701    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 702        if (n_quad[nro]!=NULL)
 703           {
 704           replaceQuad (h_quad[nro], n_quad[nro]);
 705           nouveaux->addQuad  (n_quad[nro]);
 706           }
 707
 708    for (int nro=0 ; nro<HE_MAXI ; nro++)
 709        if (n_edge[nro]!=NULL)
 710           {
 711           replaceEdge (h_edge[nro], n_edge[nro]);
 712           nouveaux->addEdge  (n_edge[nro]);
 713           }
 714
 715    for (int nro=0 ; nro<HV_MAXI ; nro++)
 716        if (n_vertex[nro]!=NULL)
 717           {
 718           replaceVertex (h_vertex[nro], n_vertex[nro]);
 719           nouveaux->addVertex  (n_vertex[nro]);
 720           }
 721
 722    if (debug())
 723       {
 724       printf (" ... Apres disconnectEdge\n");
 725       dumpFull ();
 726       }
 727
 728    nouveaux->moveDisco (this);
 729    return nouveaux;
 730 }
 731 // ========================================================= disconnectVertex
 732 Elements* Hexa::disconnectVertex (Vertex* noeud)
 733 {
 734    if (debug())
 735       {
 736       printf (" ... Avant disconnectVertex, vertex=");
 737       noeud->printName ("\n");
 738       dumpFull ();
 739       }
 740
 741    int node = findVertex (noeud);
 742    if (node==NOTHING)
 743       return NULL;
 744
 745    Vertex* new_node = new Vertex (noeud);
 746    Quad*   new_quad [HQ_MAXI];
 747    Edge*   new_edge [HE_MAXI];
 748
 749    for (int nro=0 ; nro<HE_MAXI ; nro++) new_edge [nro] = NULL;
 750    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 751        {
 752        new_quad [nro] = NULL;
 753             // Cete face contient le sommet et est commune a 2 hexas
 754        if (   h_quad[nro]->indexVertex(noeud) >= 0
 755            && h_quad[nro]->getNbrParents  ()  >= 2)
 756            {
 757            int nbmod = 0;
 758            Edge* tedge [QUAD4];
 759            for (int qed=0 ; qed<QUAD4 ; qed++)
 760                {
 761                Edge* arete = tedge[qed] = h_quad[nro]->getEdge (qed);
 762                int   indv  = arete->index (noeud);
 763                if (indv>=0)
 764                   {
 765                   nbmod++;
 766                   int hed = findEdge (arete);
 767                   if (hed<0)
 768                      return NULL;
 769                   if (new_edge [hed]==NULL)
 770                       new_edge [hed] = new Edge (new_node,
 771                                                  arete->getVertex(1-indv));
 772                   tedge [qed] = new_edge [hed];
 773                   }
 774                }
 775            if (nbmod!=2)
 776               return NULL;
 777            new_quad [nro] = new Quad (tedge[0], tedge[1], tedge[2], tedge[3]);
 778            }
 779        }
 780
 781    Elements* nouveaux  = new Elements (el_root);
 782
 783    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 784        if (new_quad [nro] != NULL)
 785           {
 786           replaceQuad (h_quad [nro], new_quad [nro]);
 787           nouveaux->addQuad (new_quad[nro]);
 788           }
 789
 790    for (int nro=0 ; nro<HE_MAXI ; nro++)
 791        if (new_edge [nro] != NULL)
 792           {
 793           replaceEdge (h_edge [nro], new_edge [nro]);
 794           nouveaux->addEdge (new_edge[nro]);
 795           }
 796
 797    replaceVertex (noeud, new_node);
 798    nouveaux->addVertex (new_node);
 799
 800
 801    if (debug())
 802       {
 803       printf (" ... Apres disconnectVertex\n");
 804       dumpFull ();
 805       }
 806
 807    nouveaux->moveDisco (this);
 808    return nouveaux;
 809 }
 810 // ========================================================= getBase
 811 int Hexa::getBase (Vertex* orig, Edge* normale)
 812 {
 813    for (int nq=0 ; nq<HQ_MAXI ; nq++)
 814        {
 815        if (   h_quad[nq]->indexVertex(orig)    >= 0
 816            && h_quad[nq]->indexEdge  (normale) < 0)
 817            return nq;
 818        }
 819    return NOTHING;
 820 }
 821 // ======================================================== replaceQuad
 822 void Hexa::replaceQuad (Quad* old, Quad* par)
 823 {
 824    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 825        {
 826        if (h_quad[nro]==old)
 827            {
 828            h_quad[nro] = par;
 829            if (debug())
 830               {
 831               printf (" Dans ");
 832               printName ();
 833               printf (" [%d], ", nro);
 834               old->printName (" est remplace par ");
 835               par->printName ("\n");
 836               }
 837            }
 838        }
 839
 840 }
 841 // ======================================================== replaceEdge
 842 void Hexa::replaceEdge (Edge* old, Edge* par)
 843 {
 844    for (int nro=0 ; nro<HE_MAXI ; nro++)
 845        {
 846        if (h_edge[nro]==old)
 847            {
 848            h_edge[nro] = par;
 849            if (debug())
 850               {
 851               printf (" Dans ");
 852               printName ();
 853               printf (" [%d], ", nro);
 854               old->printName (" est remplace par ");
 855               par->printName ("\n");
 856               }
 857            }
 858        }
 859
 860    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 861        {
 862        h_quad[nro]->replaceEdge (old, par);
 863        }
 864 }
 865 // ======================================================== replaceVertex
 866 void Hexa::replaceVertex (Vertex* old, Vertex* par)
 867 {
 868    for (int nro=0 ; nro<HV_MAXI ; nro++)
 869        {
 870        if (h_vertex [nro]==old)
 871            {
 872            h_vertex [nro] = par;
 873            if (debug())
 874               {
 875               printf (" Dans ");
 876               printName ();
 877               printf (" [%d], ", nro);
 878               old->printName (" est remplace par ");
 879               par->printName ("\n");
 880               }
 881            }
 882        }
 883
 884    for (int nro=0 ; nro<HE_MAXI ; nro++)
 885        {
 886        h_edge[nro]->replaceVertex (old, par);
 887        }
 888
 889    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 890        {
 891        h_quad[nro]->replaceVertex (old, par);
 892        }
 893 }
 894 // ======================================================== removeConnected
 895 void Hexa::removeConnected ()
 896 {
 897
 898    if (el_type == EL_REMOVED)
 899       return;
 900
 901    remove();
 902
 903    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 904        {
 905        Quad*  face = h_quad [nro];
 906        int nbhexas = face->getNbrParents ();
 907
 908        for (int nc=0 ; nc<nbhexas ; nc++)
 909            {
 910            Hexa* cell = face->getParent(nc);
 911            if (cell!=NULL && cell->isValid ())
 912               cell->removeConnected ();
 913            }
 914        }
 915
 916    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 917        h_quad [nro]->remove();
 918    for (int nro=0 ; nro<HE_MAXI ; nro++)
 919        h_edge [nro]->remove();
 920    for (int nro=0 ; nro<HV_MAXI ; nro++)
 921        h_vertex [nro]->remove();
 922 }
 923 // ======================================================== findOpposedQuad
 924 int Hexa::findOpposedQuad (Quad* face, Edge* arete)
 925 {
 926    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 927        {
 928        Quad*  quad = h_quad [nro];
 929        if (quad!=face && quad->indexEdge (arete) >=0)
 930           return nro;
 931        }
 932
 933    return NOTHING;
 934 }
 935 // ========================================================= dump
 936 void Hexa::dump ()
 937 {
 938    printName(" = (");
 939    if (NOT isHere ())
 940       {
 941       printf ("*** deleted ***)\n");
 942       return;
 943       }
 944
 945    for (int nro=0; nro<HQ_MAXI ; nro++)
 946         PrintName (h_quad[nro]);
 947    printf (")\n");
 948
 949    printf ("      = (");
 950
 951    for (int nro=0; nro<HE_MAXI ; nro++)
 952         {
 953         PrintName (h_edge[nro]);
 954         if (nro==3 || nro ==7)
 955            printf ("\n         ");
 956         }
 957    printf (")\n");
 958
 959    printf ("      = (");
 960    for (int nro=0; nro<HV_MAXI ; nro++)
 961         PrintName (h_vertex[nro]);
 962    printf (")\n");
 963    Real3 cg;
 964    getCenter (cg);
 965    printf ("cg    = (%g, %g, %g)\n", cg[0], cg[1], cg[2]);
 966
 967 }
 968 // ======================================================== dumpPlus
 969 void Hexa::dumpPlus ()
 970 {
 971    dump ();
 972    for (int nro=0 ; nro < HV_MAXI ; nro++)
 973        {
 974        Vertex* pv = h_vertex[nro];
 975        printf ( "    ");
 976        if (pv!=NULL)
 977           {
 978           pv->printName ("");
 979           printf ( " (%g, %g, %g)\n", pv->getX(),  pv->getY(),  pv->getZ());
 980           }
 981        else
 982           {
 983           printf ( "NULL\n");
 984           }
 985        }
 986 }
 987 // ======================================================== dumpFull
 988 void Hexa::dumpFull ()
 989 {
 990    dump ();
 991    Globale* glob = Globale::getInstance ();
 992
 993    printf ("\n");
 994    for (int nro=0; nro<HQ_MAXI ; nro++)
 995        {
 996        printf (" quad(%s) = ", glob->namofHexaQuad(nro));
 997        if (h_quad[nro] ==NULL)
 998            printf (" NULL\n");
 999        else
1000            {
1001            h_quad[nro]->printName (" = (");
1002            for (int nc=0; nc<QUAD4 ; nc++)
1003                 h_quad[nro]->getEdge(nc)->printName ();
1004            printf (")\n");
1005            printf ("                   = (");
1006            for (int nc=0; nc<QUAD4 ; nc++)
1007                 h_quad[nro]->getVertex(nc)->printName ();
1008            printf (")\n");
1009            }
1010        }
1011
1012    printf ("\n");
1013    for (int nro=0; nro<HE_MAXI ; nro++)
1014        {
1015        printf (" edge(%s) = ", glob->namofHexaEdge(nro));
1016        if (h_edge[nro] ==NULL)
1017            printf (" NULL\n");
1018        else
1019            {
1020            h_edge[nro]->printName (" = (");
1021            for (int nc=0; nc<V_TWO ; nc++)
1022                 h_edge[nro]->getVertex(nc)->printName ();
1023            printf (")\n");
1024            }
1025        }
1026    printf ("\n");
1027
1028    for (int nro=0; nro<HV_MAXI ; nro++)
1029        {
1030        Vertex* pv = h_vertex[nro];
1031        printf (" vertex(%s) = ", glob->namofHexaVertex(nro));
1032        if (pv ==NULL)
1033            printf (" NULL");
1034        else
1035            {
1036            pv->printName (" = (");
1037            printf ("%g, %g, %g)\n", pv->getX(), pv->getY(), pv->getZ());
1038            }
1039        }
1040    printf ("\n");
1041 }
1042 // ======================================================== getOpposedQuad
1043 Quad* Hexa::getOpposedQuad (Quad* face)
1044 {
1045    if      (face == h_quad [Q_A]) return h_quad [Q_B];
1046    else if (face == h_quad [Q_B]) return h_quad [Q_A];
1047    else if (face == h_quad [Q_C]) return h_quad [Q_D];
1048    else if (face == h_quad [Q_D]) return h_quad [Q_C];
1049    else if (face == h_quad [Q_E]) return h_quad [Q_F];
1050    else if (face == h_quad [Q_F]) return h_quad [Q_F];
1051    else                           return NULL;
1052 }
1053 // ========================================================= findQuad
1054 Quad* Hexa::findQuad (Edge* ed1, Edge* ed2)
1055 {
1056    for (int nro=0 ; nro<HQ_MAXI ; nro++)
1057        {
1058        if (   h_quad[nro]->indexEdge (ed1) >= 0
1059            && h_quad[nro]->indexEdge (ed2) >= 0)
1060           return h_quad [nro];
1061        }
1062
1063    return NULL;
1064 }
1065 // ========================================================= findEdge
1066 Edge* Hexa::findEdge (Vertex* v1, Vertex* v2)
1067 {
1068    for (int nro=0 ; nro<HE_MAXI ; nro++)
1069        {
1070        if (   h_edge[nro]->index (v1) >= 0
1071            && h_edge[nro]->index (v2) >= 0)
1072           return h_edge [nro];
1073        }
1074
1075    return NULL;
1076 }
1077 // ====================================================== getPerpendicularEdge
1078 Edge* Hexa::getPerpendicularEdge (Quad* quad, Vertex* vertex)
1079 {
1080    if (quad==NULL || vertex==NULL)
1081       return NULL;
1082
1083    int nv = quad->indexVertex (vertex);
1084    int nq = findQuad (quad);
1085    if (nv<0 || nq<0)
1086       return NULL;
1087
1088    for (int nro=0 ; nro<HE_MAXI ; nro++)
1089        {
1090        if (quad->indexEdge (h_edge[nro])<0 && h_edge[nro]->index(vertex)>=0)
1091           return h_edge [nro];
1092        }
1093
1094    return NULL;
1095 }
1096 // ============================================================  getQuad
1097 Quad* Hexa::getQuad (int nro)
1098 {
1099    Quad* elt = NULL;
1100    if (nro >=0 && nro < HQ_MAXI && el_status == HOK && h_quad [nro]->isValid())
1101       elt = h_quad [nro];
1102
1103    return elt;
1104 }
1105 // ============================================================  getEdge
1106 Edge* Hexa::getEdge (int nro)
1107 {
1108    Edge* elt = NULL;
1109    if (nro >=0 && nro < HE_MAXI && el_status == HOK && h_edge [nro]->isValid())
1110       elt = h_edge [nro];
1111
1112    return elt;
1113 }
1114 // ============================================================  getVertex
1115 Vertex* Hexa::getVertex (int nro)
1116 {
1117    Vertex* elt = NULL;
1118    if (nro >=0 && nro <  HV_MAXI && el_status == HOK && h_vertex [nro]->isValid())
1119       elt = h_vertex [nro];
1120
1121    return elt;
1122 }
1123 // ============================================================  getCenter
1124 double* Hexa::getCenter (double centre[])
1125 {
1126    centre [dir_x] = centre [dir_y] = centre [dir_z] = 0;
1127
1128    for (int nv=0 ; nv<HV_MAXI ; nv++)
1129        {
1130        centre [dir_x] += h_vertex[nv]->getX ();
1131        centre [dir_y] += h_vertex[nv]->getY ();
1132        centre [dir_z] += h_vertex[nv]->getZ ();
1133        }
1134
1135    centre [dir_x] /= HV_MAXI;
1136    centre [dir_y] /= HV_MAXI;
1137    centre [dir_z] /= HV_MAXI;
1138    return centre;
1139 }
1140 // =============================================================== definedBy
1141 bool Hexa::definedBy  (Vertex* v1, Vertex* v2)
1142 {
1143    for (int n1=0 ; n1< HV_MAXI ; n1++)
1144        {
1145 //              (   Diagonale        )  Dessus
1146        int n2 = (n1 + 2) MODULO HV_MAXI + HV_MAXI;
1147        if (   (v1 == h_vertex[n1] && v2 == h_vertex[n2])
1148            || (v1 == h_vertex[n2] && v2 == h_vertex[n1])) return true;
1149        }
1150    return false;
1151 }
1152 // =============================================================== definedBy
1153 bool Hexa::definedBy  (Quad* qa, Quad* qb)
1154 {
1155    if (qa==qb || BadElement (qa) || BadElement (qb))
1156        return false;
1157
1158    bool p1 = false, p2 = false;
1159    for (int nq=0 ; nq< HQ_MAXI ; nq++)
1160        {
1161        if (qa == h_quad[nq])
1162           p1 = true;
1163        else if (qb == h_quad[nq])
1164           p2 = true;
1165        }
1166    return p1 && p2;
1167 }
1168 // =============================================================== setColor
1169 void Hexa::setColor  (double val)
1170 {
1171    for (int nc=0 ; nc< HV_MAXI ; nc++)
1172        h_vertex[nc] -> setColor (val);
1173 }
1174 // ============================================================== markElements
1175 void Hexa::markElements  (int marque)
1176 {
1177    for (int nc=0 ; nc< HQ_MAXI ; nc++) h_quad  [nc] -> setMark (marque);
1178    for (int nc=0 ; nc< HE_MAXI ; nc++) h_edge  [nc] -> setMark (marque);
1179    for (int nc=0 ; nc< HV_MAXI ; nc++) h_vertex[nc] -> setMark (marque);
1180 }
1181 // =============================================================== duplicate
1182 void Hexa::duplicate  ()
1183 {
1184    h_clone = new Hexa (GetClone (h_quad [Q_A]),
1185                        GetClone (h_quad [Q_B]),
1186                        GetClone (h_quad [Q_C]),
1187                        GetClone (h_quad [Q_D]),
1188                        GetClone (h_quad [Q_E]),
1189                        GetClone (h_quad [Q_F]));
1190 }
1191 END_NAMESPACE_HEXA