src/HEXABLOCK/HexHexa.cxx

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