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 }
  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 // ======================================================= printHexaVtk
 343 void Hexa::printHexaVtk (pfile vtk)
 344 {
 345    fprintf (vtk, "%d", HV_MAXI);
 346
 347    fprintf (vtk, " %d", h_vertex[V_ACE]->getId ());
 348    fprintf (vtk, " %d", h_vertex[V_ACF]->getId ());
 349    fprintf (vtk, " %d", h_vertex[V_ADF]->getId ());
 350    fprintf (vtk, " %d", h_vertex[V_ADE]->getId ());
 351
 352    fprintf (vtk, " %d", h_vertex[V_BCE]->getId ());
 353    fprintf (vtk, " %d", h_vertex[V_BCF]->getId ());
 354    fprintf (vtk, " %d", h_vertex[V_BDF]->getId ());
 355    fprintf (vtk, " %d", h_vertex[V_BDE]->getId ());
 356
 357    fprintf (vtk, "\n");
 358 }
 359 // ======================================================== hasFreEdges
 360 bool Hexa::hasFreeEdges ()
 361 {
 362    if (isDeleted())
 363        return false;
 364
 365    for (int nro=0; nro<HE_MAXI ; nro++)
 366        if (h_edge[nro]->getMark()<0)
 367           return true;
 368
 369    return false;
 370 }
 371 // ======================================================== propager
 372 void Hexa::propager (Propagation* prop, int nro)
 373 {
 374    if (isDeleted())
 375        return;
 376
 377    for (int nro=0; nro<HE_MAXI ; nro++)
 378        if (h_edge[nro]->getMark()<0)
 379            h_edge[nro]->propager (prop, nro);
 380 }
 381 // ========================================================= saveXml
 382 void Hexa::saveXml (XmlWriter* xml)
 383 {
 384    char ident[12];
 385    string quads;
 386
 387    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 388        {
 389        if (nro>0) quads += " ";
 390        quads += h_quad[nro]->getName(ident);
 391        }
 392
 393    getName (ident);
 394    xml->openMark     ("Hexa");
 395    xml->addAttribute ("id",    ident);
 396    xml->addAttribute ("quads", quads);
 397    if (el_name!=ident)
 398        xml->addAttribute ("name", el_name);
 399    xml->closeMark ();
 400 }
 401 // ========================================================= findQuad
 402 int Hexa::findQuad (Quad* element)
 403 {
 404    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 405        if (h_quad[nro]==element)
 406           return nro;
 407
 408    return NOTHING;
 409 }
 410 // ========================================================= findEdge
 411 int Hexa::findEdge (Edge* element)
 412 {
 413    for (int nro=0 ; nro<HE_MAXI ; nro++)
 414        if (h_edge[nro]==element)
 415           return nro;
 416
 417    return NOTHING;
 418 }
 419 // ========================================================= findVertex
 420 int Hexa::findVertex (Vertex* element)
 421 {
 422    for (int nro=0 ; nro<HV_MAXI ; nro++)
 423        if (h_vertex[nro]==element)
 424           return nro;
 425
 426    return NOTHING;
 427 }
 428 // ========================================================= disconnectQuad
 429 Elements* Hexa::disconnectQuad (Quad* quad)
 430 {
 431    if (quad==NULL)
 432       return NULL;
 433
 434    if (debug())
 435       {
 436       printf (" ... Avant disconnectQuad, quad=");
 437       quad->printName ("\n");
 438       dumpFull ();
 439       }
 440
 441    int nface = findQuad (quad);
 442    if (nface==NOTHING)
 443       return NULL;
 444                                        // Face opposee : replace
 445    // int nfopp = (nface MODULO 2==0) ? nface+1 : nface-1;
 446
 447    int  ind_edge  [QUAD4], ind_opp_quad [QUAD4];
 448    bool make_quad [QUAD4], make_edge [QUAD4];
 449
 450    for (int nro=0 ; nro<QUAD4 ; nro++)
 451        make_quad [nro] = make_edge[nro] = false;
 452
 453    for (int nro=0 ; nro<QUAD4 ; nro++)
 454        {
 455        int nro1  = (nro+1) MODULO QUAD4;
 456        int pedge = findEdge   (quad->getEdge   (nro));
 457        int pnode = findVertex (quad->getVertex (nro));
 458        int oppq  = findOpposedQuad (quad, quad->getEdge (nro));
 459
 460        ind_edge [nro]     = pedge;
 461        ind_opp_quad [nro] = oppq;
 462
 463        if (pedge==NOTHING || pnode==NOTHING || oppq==NOTHING)
 464           return NULL;
 465
 466        make_quad [nro]  = h_quad[oppq]->getNbrParents() == 2;
 467        make_edge [nro ] = make_edge [nro ] || make_quad [nro];
 468        make_edge [nro1] = make_edge [nro1] || make_quad [nro];
 469
 470        if (debug())
 471           {
 472           printf (" Sommet nro %d : ", nro);
 473           quad->getVertex(nro)->printName (", ");
 474           printf (" edge = ");
 475           quad->getEdge(nro)->printName (", ");
 476           printf (" quad oppose = ");
 477           h_quad[oppq]->printName("");
 478           if (make_quad [nro])
 479              printf (" a dissocier\n");
 480           else
 481              printf ("\n");
 482
 483           }
 484        }
 485
 486    Vertex* new_node     [QUAD4];
 487    Edge*   new_opp_edge [QUAD4];
 488    Edge*   old_opp_edge [QUAD4];
 489
 490    for (int nro=0 ; nro<QUAD4 ; nro++)
 491        {
 492        old_opp_edge [nro] = NULL;
 493        new_opp_edge [nro] = NULL;
 494        Vertex*  o_v0  = quad->getVertex (nro);
 495        new_node [nro] = new Vertex (o_v0);
 496        if (debug())
 497           {
 498           printf (" quad.vertex [%d] = ", nro);
 499           quad->getVertex (nro)->printName (" --> ");
 500           new_node [nro]->printName ("\n");
 501           }
 502
 503        if (make_edge[nro])
 504           {
 505           Quad*   pface  = h_quad [ind_opp_quad [nro]];
 506           int     bid;
 507           int     ncut = pface->inter (quad, bid);
 508           Edge*   ecut = pface->getEdge ((ncut+1) MODULO QUAD4);
 509           Vertex* vopp = ecut->getVertex(V_AMONT);
 510           if (vopp==o_v0)
 511               vopp = ecut->getVertex (V_AVAL);
 512           else if (o_v0 != ecut->getVertex (V_AVAL));
 513               {
 514               ecut = pface->getEdge ((ncut+3) MODULO QUAD4);
 515               vopp = ecut->getVertex(V_AMONT);
 516               if (vopp==o_v0)
 517                   vopp = ecut->getVertex (V_AVAL);
 518               else if (o_v0 != ecut->getVertex (V_AVAL))
 519                   return NULL;
 520               }
 521
 522           old_opp_edge [nro] = ecut;
 523           new_opp_edge [nro] = new Edge (new_node[nro], vopp);
 524           if (debug())
 525              {
 526              printf (" quad.opp_edge [%d] = ", nro);
 527              old_opp_edge [nro]->printName (" --> ");
 528              new_opp_edge [nro]->printName ("\n");
 529              }
 530           }
 531        }
 532
 533    Quad* new_quad = new Quad (new_node[0], new_node[1], new_node[2],
 534                                                         new_node[3]);
 535
 536    Quad* new_opp_quad [QUAD4];
 537    Quad* old_opp_quad [QUAD4];
 538    for (int nro=0 ; nro<QUAD4 ; nro++)
 539        {
 540        old_opp_quad [nro] = NULL;
 541        new_opp_quad [nro] = NULL;
 542        if (make_quad[nro])
 543           {
 544           int nro1 = (nro+1) MODULO QUAD4;
 545
 546           Edge* n_edge0 = new_quad->getEdge (nro);
 547           Edge* n_edge1 = new_opp_edge [nro];
 548           Edge* n_edge3 = new_opp_edge [nro1];
 549
 550           int iv1 = n_edge1->inter (n_edge0);
 551           int iv3 = n_edge3->inter (n_edge0);
 552           if (iv1 <0 || iv3 <0)
 553              return NULL;
 554
 555           Quad* o_face = h_quad [ind_opp_quad [nro]];
 556           Edge* edge2  = o_face->findEdge (n_edge1->getVertex (1-iv1),
 557                                            n_edge3->getVertex (1-iv3));
 558           if (edge2==NULL)
 559              return NULL;
 560           // Edge* o_edge0 = h_edge [ind_edge     [nro]];
 561           // int sens = 1;
 562           // Edge* edge2   = o_face->getOpposEdge (o_edge0, sens);
 563
 564           old_opp_quad [nro] = o_face;
 565           if (debug())
 566              printf (" -------- Quad oppose nro %d\n", nro);
 567           new_opp_quad [nro] = new Quad (n_edge0, n_edge1, edge2, n_edge3);
 568           }
 569        }
 570
 571    replaceQuad (quad, new_quad);
 572    for (int nro=0 ; nro<QUAD4 ; nro++)
 573        if (make_quad[nro])
 574           replaceQuad (old_opp_quad [nro], new_opp_quad [nro]);
 575
 576    for (int nro=0 ; nro<QUAD4 ; nro++)
 577        {
 578        replaceEdge (h_edge[ind_edge[nro]], new_quad->getEdge (nro));
 579        if (make_edge[nro])
 580           replaceEdge (old_opp_edge [nro], new_opp_edge [nro]);
 581        }
 582
 583    for (int nro=0 ; nro<QUAD4 ; nro++)
 584        {
 585        replaceVertex (quad->getVertex(nro), new_node[nro]);
 586        }
 587
 588
 589    h_quad [nface] = new_quad;
 590    if (debug())
 591       {
 592       printf (" ... Apres disconnectQuad, new_quad=");
 593       new_quad->printName ("\n");
 594       dumpFull ();
 595       }
 596
 597    Elements* nouveaux  = new Elements (el_root);
 598    nouveaux->addQuad   (new_quad);
 599    for (int nro=0 ; nro<QUAD4 ; nro++)
 600        {
 601        nouveaux->addEdge   (new_quad->getEdge   (nro));
 602        nouveaux->addVertex (new_quad->getVertex (nro));
 603        if (make_edge[nro])
 604           nouveaux->addEdge (new_opp_edge [nro]);
 605        if (make_quad[nro])
 606           nouveaux->addQuad (new_opp_quad [nro]);
 607        }
 608    nouveaux->moveDisco (this);
 609    return nouveaux;
 610 }
 611 // ========================================================= disconnectEdge
 612 Elements* Hexa::disconnectEdge (Edge* arete)
 613 {
 614    int nedge  = findEdge   (arete);
 615    int namont = findVertex (arete->getVertex(V_AMONT));
 616    int naval  = findVertex (arete->getVertex(V_AVAL ));
 617
 618    if (nedge==NOTHING || namont==NOTHING || naval==NOTHING)
 619       return NULL;
 620
 621    if (debug())
 622       {
 623       printf (" ... Avant disconnectEdge, arete=");
 624       arete->printName ("\n");
 625       dumpFull ();
 626       }
 627
 628    Edge*   n_edge   [HE_MAXI];
 629    Quad*   n_quad   [HQ_MAXI];
 630    Vertex* n_vertex [HV_MAXI];
 631
 632    for (int nro=0 ; nro<HQ_MAXI ; nro++) n_quad [nro]   = NULL;
 633    for (int nro=0 ; nro<HE_MAXI ; nro++) n_edge [nro]   = NULL;
 634    for (int nro=0 ; nro<HV_MAXI ; nro++) n_vertex [nro] = NULL;
 635
 636    Vertex* old_amont = arete->getVertex (V_AMONT);
 637    Vertex* old_aval  = arete->getVertex (V_AVAL );
 638    Vertex* new_amont = n_vertex [namont] = new Vertex (old_amont);
 639    Vertex* new_aval  = n_vertex [naval]  = new Vertex (old_aval);
 640    n_edge [nedge]    = new Edge   (new_amont, new_aval);
 641
 642    // Un edge non remplace, qui contient un vertex remplace
 643    //         commun a plus de 2 faces (donc appartenant a un autre hexa)
 644    //         doit etre duplique
 645
 646    for (int nro=0 ; nro<HE_MAXI ; nro++)
 647        {
 648        if (   n_edge[nro]==NULL && h_edge[nro] != NULL
 649            && h_edge[nro]->getNbrParents()>2)
 650           {
 651           Vertex* va =  h_edge[nro]->getVertex (V_AMONT);
 652           Vertex* vb =  h_edge[nro]->getVertex (V_AVAL);
 653
 654           if (va==old_amont)
 655              n_edge [nro] = new Edge (new_amont, vb);
 656           else if (va==old_aval)
 657              n_edge [nro] = new Edge (new_aval,  vb);
 658           else if (vb==old_amont)
 659              n_edge [nro] = new Edge (va, new_amont);
 660           else if (vb==old_aval)
 661              n_edge [nro] = new Edge (va, new_aval);
 662           }
 663        }
 664
 665    // Un quad non remplace, qui contient un edge remplace
 666    //         commun a plus de 2 Hexas
 667    //         doit etre duplique
 668
 669    Globale* glob = Globale::getInstance();
 670    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 671        if (   n_quad[nro]==NULL && h_quad[nro] != NULL
 672            && h_quad[nro]->getNbrParents()>1)
 673           {
 674           Edge* qedge[QUAD4];
 675           bool  duplic = false;
 676           for (int ned=0 ; ned<QUAD4 ; ned++)
 677               {
 678               int ndup = glob->QuadEdge (nro, (EnumQuad)ned);
 679               if (n_edge [ndup] ==NULL)
 680                  qedge [ned] = h_edge[ndup];
 681               else
 682                  {
 683                  qedge [ned] = n_edge[ndup];
 684                  duplic = true;
 685                  }
 686               }
 687           if (duplic)
 688              n_quad[nro] = new Quad (qedge[Q_A], qedge[Q_B],
 689                                      qedge[Q_C], qedge[Q_D]);
 690           }
 691
 692    Elements* nouveaux  = new Elements (el_root);
 693    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 694        if (n_quad[nro]!=NULL)
 695           {
 696           replaceQuad (h_quad[nro], n_quad[nro]);
 697           nouveaux->addQuad  (n_quad[nro]);
 698           }
 699
 700    for (int nro=0 ; nro<HE_MAXI ; nro++)
 701        if (n_edge[nro]!=NULL)
 702           {
 703           replaceEdge (h_edge[nro], n_edge[nro]);
 704           nouveaux->addEdge  (n_edge[nro]);
 705           }
 706
 707    for (int nro=0 ; nro<HV_MAXI ; nro++)
 708        if (n_vertex[nro]!=NULL)
 709           {
 710           replaceVertex (h_vertex[nro], n_vertex[nro]);
 711           nouveaux->addVertex  (n_vertex[nro]);
 712           }
 713
 714    if (debug())
 715       {
 716       printf (" ... Apres disconnectEdge\n");
 717       dumpFull ();
 718       }
 719
 720    nouveaux->moveDisco (this);
 721    return nouveaux;
 722 }
 723 // ========================================================= disconnectVertex
 724 Elements* Hexa::disconnectVertex (Vertex* noeud)
 725 {
 726    if (debug())
 727       {
 728       printf (" ... Avant disconnectVertex, vertex=");
 729       noeud->printName ("\n");
 730       dumpFull ();
 731       }
 732
 733    int node = findVertex (noeud);
 734    if (node==NOTHING)
 735       return NULL;
 736
 737    Vertex* new_node = new Vertex (noeud);
 738    Quad*   new_quad [HQ_MAXI];
 739    Edge*   new_edge [HE_MAXI];
 740
 741    for (int nro=0 ; nro<HE_MAXI ; nro++) new_edge [nro] = NULL;
 742    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 743        {
 744        new_quad [nro] = NULL;
 745             // Cete face contient le sommet et est commune a 2 hexas
 746        if (   h_quad[nro]->indexVertex(noeud) >= 0
 747            && h_quad[nro]->getNbrParents  ()  >= 2)
 748            {
 749            int nbmod = 0;
 750            Edge* tedge [QUAD4];
 751            for (int qed=0 ; qed<QUAD4 ; qed++)
 752                {
 753                Edge* arete = tedge[qed] = h_quad[nro]->getEdge (qed);
 754                int   indv  = arete->index (noeud);
 755                if (indv>=0)
 756                   {
 757                   nbmod++;
 758                   int hed = findEdge (arete);
 759                   if (hed<0)
 760                      return NULL;
 761                   if (new_edge [hed]==NULL)
 762                       new_edge [hed] = new Edge (new_node,
 763                                                  arete->getVertex(1-indv));
 764                   tedge [qed] = new_edge [hed];
 765                   }
 766                }
 767            if (nbmod!=2)
 768               return NULL;
 769            new_quad [nro] = new Quad (tedge[0], tedge[1], tedge[2], tedge[3]);
 770            }
 771        }
 772
 773    Elements* nouveaux  = new Elements (el_root);
 774
 775    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 776        if (new_quad [nro] != NULL)
 777           {
 778           replaceQuad (h_quad [nro], new_quad [nro]);
 779           nouveaux->addQuad (new_quad[nro]);
 780           }
 781
 782    for (int nro=0 ; nro<HE_MAXI ; nro++)
 783        if (new_edge [nro] != NULL)
 784           {
 785           replaceEdge (h_edge [nro], new_edge [nro]);
 786           nouveaux->addEdge (new_edge[nro]);
 787           }
 788
 789    replaceVertex (noeud, new_node);
 790    nouveaux->addVertex (new_node);
 791
 792
 793    if (debug())
 794       {
 795       printf (" ... Apres disconnectVertex\n");
 796       dumpFull ();
 797       }
 798
 799    nouveaux->moveDisco (this);
 800    return nouveaux;
 801 }
 802 // ========================================================= getBase
 803 int Hexa::getBase (Vertex* orig, Edge* normale)
 804 {
 805    for (int nq=0 ; nq<HQ_MAXI ; nq++)
 806        {
 807        if (   h_quad[nq]->indexVertex(orig)    >= 0
 808            && h_quad[nq]->indexEdge  (normale) < 0)
 809            return nq;
 810        }
 811    return NOTHING;
 812 }
 813 // ======================================================== replaceQuad
 814 void Hexa::replaceQuad (Quad* old, Quad* par)
 815 {
 816    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 817        {
 818        if (h_quad[nro]==old)
 819            {
 820            h_quad[nro] = par;
 821            if (debug())
 822               {
 823               printf (" Dans ");
 824               printName ();
 825               printf (" [%d], ", nro);
 826               old->printName (" est remplace par ");
 827               par->printName ("\n");
 828               }
 829            }
 830        }
 831
 832 }
 833 // ======================================================== replaceEdge
 834 void Hexa::replaceEdge (Edge* old, Edge* par)
 835 {
 836    for (int nro=0 ; nro<HE_MAXI ; nro++)
 837        {
 838        if (h_edge[nro]==old)
 839            {
 840            h_edge[nro] = par;
 841            if (debug())
 842               {
 843               printf (" Dans ");
 844               printName ();
 845               printf (" [%d], ", nro);
 846               old->printName (" est remplace par ");
 847               par->printName ("\n");
 848               }
 849            }
 850        }
 851
 852    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 853        {
 854        h_quad[nro]->replaceEdge (old, par);
 855        }
 856 }
 857 // ======================================================== replaceVertex
 858 void Hexa::replaceVertex (Vertex* old, Vertex* par)
 859 {
 860    for (int nro=0 ; nro<HV_MAXI ; nro++)
 861        {
 862        if (h_vertex [nro]==old)
 863            {
 864            h_vertex [nro] = par;
 865            if (debug())
 866               {
 867               printf (" Dans ");
 868               printName ();
 869               printf (" [%d], ", nro);
 870               old->printName (" est remplace par ");
 871               par->printName ("\n");
 872               }
 873            }
 874        }
 875
 876    for (int nro=0 ; nro<HE_MAXI ; nro++)
 877        {
 878        h_edge[nro]->replaceVertex (old, par);
 879        }
 880
 881    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 882        {
 883        h_quad[nro]->replaceVertex (old, par);
 884        }
 885 }
 886 // ======================================================== removeConnected
 887 void Hexa::removeConnected ()
 888 {
 889
 890    if (el_type == EL_REMOVED)
 891       return;
 892
 893    remove();
 894
 895    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 896        {
 897        Quad*  face = h_quad [nro];
 898        int nbhexas = face->getNbrParents ();
 899
 900        for (int nc=0 ; nc<nbhexas ; nc++)
 901            {
 902            Hexa* cell = face->getParent(nc);
 903            if (cell!=NULL && cell->isValid ())
 904               cell->removeConnected ();
 905            }
 906        }
 907
 908    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 909        h_quad [nro]->remove();
 910    for (int nro=0 ; nro<HE_MAXI ; nro++)
 911        h_edge [nro]->remove();
 912    for (int nro=0 ; nro<HV_MAXI ; nro++)
 913        h_vertex [nro]->remove();
 914 }
 915 // ======================================================== findOpposedQuad
 916 int Hexa::findOpposedQuad (Quad* face, Edge* arete)
 917 {
 918    for (int nro=0 ; nro<HQ_MAXI ; nro++)
 919        {
 920        Quad*  quad = h_quad [nro];
 921        if (quad!=face && quad->indexEdge (arete) >=0)
 922           return nro;
 923        }
 924
 925    return NOTHING;
 926 }
 927 // ========================================================= dump
 928 void Hexa::dump ()
 929 {
 930    printName(" = (");
 931    if (NOT isHere ())
 932       {
 933       printf ("*** deleted ***)\n");
 934       return;
 935       }
 936
 937    for (int nro=0; nro<HQ_MAXI ; nro++)
 938         PrintName (h_quad[nro]);
 939    printf (")\n");
 940
 941    printf ("      = (");
 942
 943    for (int nro=0; nro<HE_MAXI ; nro++)
 944         {
 945         PrintName (h_edge[nro]);
 946         if (nro==3 || nro ==7)
 947            printf ("\n         ");
 948         }
 949    printf (")\n");
 950
 951    printf ("      = (");
 952    for (int nro=0; nro<HV_MAXI ; nro++)
 953         PrintName (h_vertex[nro]);
 954    printf (")\n");
 955    Real3 cg;
 956    getCenter (cg);
 957    printf ("cg    = (%g, %g, %g)\n", cg[0], cg[1], cg[2]);
 958
 959 }
 960 // ======================================================== dumpPlus
 961 void Hexa::dumpPlus ()
 962 {
 963    dump ();
 964    for (int nro=0 ; nro < HV_MAXI ; nro++)
 965        {
 966        Vertex* pv = h_vertex[nro];
 967        printf ( "    ");
 968        if (pv!=NULL)
 969           {
 970           pv->printName ("");
 971           printf ( " (%g, %g, %g)\n", pv->getX(),  pv->getY(),  pv->getZ());
 972           }
 973        else
 974           {
 975           printf ( "NULL\n");
 976           }
 977        }
 978 }
 979 // ======================================================== dumpFull
 980 void Hexa::dumpFull ()
 981 {
 982    dump ();
 983    Globale* glob = Globale::getInstance ();
 984
 985    printf ("\n");
 986    for (int nro=0; nro<HQ_MAXI ; nro++)
 987        {
 988        printf (" quad(%s) = ", glob->namofHexaQuad(nro));
 989        if (h_quad[nro] ==NULL)
 990            printf (" NULL\n");
 991        else
 992            {
 993            h_quad[nro]->printName (" = (");
 994            for (int nc=0; nc<QUAD4 ; nc++)
 995                 h_quad[nro]->getEdge(nc)->printName ();
 996            printf (")\n");
 997            printf ("                   = (");
 998            for (int nc=0; nc<QUAD4 ; nc++)
 999                 h_quad[nro]->getVertex(nc)->printName ();
1000            printf (")\n");
1001            }
1002        }
1003
1004    printf ("\n");
1005    for (int nro=0; nro<HE_MAXI ; nro++)
1006        {
1007        printf (" edge(%s) = ", glob->namofHexaEdge(nro));
1008        if (h_edge[nro] ==NULL)
1009            printf (" NULL\n");
1010        else
1011            {
1012            h_edge[nro]->printName (" = (");
1013            for (int nc=0; nc<V_TWO ; nc++)
1014                 h_edge[nro]->getVertex(nc)->printName ();
1015            printf (")\n");
1016            }
1017        }
1018    printf ("\n");
1019
1020    for (int nro=0; nro<HV_MAXI ; nro++)
1021        {
1022        Vertex* pv = h_vertex[nro];
1023        printf (" vertex(%s) = ", glob->namofHexaVertex(nro));
1024        if (pv ==NULL)
1025            printf (" NULL");
1026        else
1027            {
1028            pv->printName (" = (");
1029            printf ("%g, %g, %g)\n", pv->getX(), pv->getY(), pv->getZ());
1030            }
1031        }
1032    printf ("\n");
1033 }
1034 // ======================================================== getOpposedQuad
1035 Quad* Hexa::getOpposedQuad (Quad* face)
1036 {
1037    if      (face == h_quad [Q_A]) return h_quad [Q_B];
1038    else if (face == h_quad [Q_B]) return h_quad [Q_A];
1039    else if (face == h_quad [Q_C]) return h_quad [Q_D];
1040    else if (face == h_quad [Q_D]) return h_quad [Q_C];
1041    else if (face == h_quad [Q_E]) return h_quad [Q_F];
1042    else if (face == h_quad [Q_F]) return h_quad [Q_F];
1043    else                           return NULL;
1044 }
1045 // ========================================================= findQuad
1046 Quad* Hexa::findQuad (Edge* ed1, Edge* ed2)
1047 {
1048    for (int nro=0 ; nro<HQ_MAXI ; nro++)
1049        {
1050        if (   h_quad[nro]->indexEdge (ed1) >= 0
1051            && h_quad[nro]->indexEdge (ed2) >= 0)
1052           return h_quad [nro];
1053        }
1054
1055    return NULL;
1056 }
1057 // ========================================================= findEdge
1058 Edge* Hexa::findEdge (Vertex* v1, Vertex* v2)
1059 {
1060    for (int nro=0 ; nro<HE_MAXI ; nro++)
1061        {
1062        if (   h_edge[nro]->index (v1) >= 0
1063            && h_edge[nro]->index (v2) >= 0)
1064           return h_edge [nro];
1065        }
1066
1067    return NULL;
1068 }
1069 // ====================================================== getPerpendicularEdge
1070 Edge* Hexa::getPerpendicularEdge (Quad* quad, Vertex* vertex)
1071 {
1072    if (quad==NULL || vertex==NULL)
1073       return NULL;
1074
1075    int nv = quad->indexVertex (vertex);
1076    int nq = findQuad (quad);
1077    if (nv<0 || nq<0)
1078       return NULL;
1079
1080    for (int nro=0 ; nro<HE_MAXI ; nro++)
1081        {
1082        if (quad->indexEdge (h_edge[nro])<0 && h_edge[nro]->index(vertex)>=0)
1083           return h_edge [nro];
1084        }
1085
1086    return NULL;
1087 }
1088 END_NAMESPACE_HEXA