src/HEXABLOCK/HexQuad.cxx

   1
   2 // C++ : Gestion des Quadrangles
   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 "HexQuad.hxx"
  23
  24 #include "HexDocument.hxx"
  25 #include "HexHexa.hxx"
  26 #include "HexElements.hxx"
  27 #include "HexGlobale.hxx"
  28
  29 #include "HexXmlWriter.hxx"
  30 #include "HexNewShape.hxx"
  31 #include "HexFaceShape.hxx"
  32
  33 BEGIN_NAMESPACE_HEXA
  34
  35 bool db = false;
  36
  37 // ======================================================== Constructeur
  38 Quad::Quad (Vertex* va, Vertex* vb, Vertex* vc, Vertex* vd)
  39     : EltBase (va->dad(), EL_QUAD)
  40 {
  41    q_vertex [E_A] = va;
  42    q_vertex [E_B] = vb;
  43    q_vertex [E_C] = vc;
  44    q_vertex [E_D] = vd;
  45    q_clone        = NULL;
  46    q_orientation  = Q_UNDEFINED;
  47
  48    for (int nro=0 ; nro<QUAD4 ; nro++)
  49        {
  50        q_edge [nro] = new Edge (q_vertex[nro],
  51                                 q_vertex[(nro+1) MODULO QUAD4]);
  52
  53        if (BadElement (q_vertex [nro]) || BadElement (q_edge [nro]))
  54           setError ();
  55        else
  56           for (int nv=nro+1 ; nv<QUAD4 ; nv++)
  57               if (q_vertex[nv] == q_vertex[nro])
  58                  setError ();
  59        }
  60
  61    majReferences ();
  62 }
  63 // ======================================================== Constructeur bis
  64 Quad::Quad (Edge* ea, Edge* eb, Edge* ec, Edge* ed)
  65     : EltBase (ea->dad(), EL_QUAD)
  66 {
  67    q_edge [E_A] = ea;
  68    q_edge [E_B] = eb;
  69    q_edge [E_C] = ec;
  70    q_edge [E_D] = ed;
  71    q_clone       = NULL;
  72    q_orientation = Q_UNDEFINED;
  73
  74    for (int nro=0 ; nro<QUAD4 ; nro++)
  75        {
  76        int prec = (nro+1) MODULO QUAD4;
  77        Vertex* node = NULL;
  78
  79        if (BadElement (q_edge[nro]))
  80           setError ();
  81        else
  82           {
  83           for (int nv=nro+1 ; nv<QUAD4 ; nv++)
  84               if (q_edge[nv] == q_edge[nro])
  85                   setError ();
  86           int nc  = q_edge[nro] -> inter (q_edge[prec]);
  87           if (nc>=0)
  88              node = q_edge[nro]->getVertex (nc);
  89           else
  90              setError (888);
  91           }
  92        q_vertex [prec] = node;
  93        }
  94
  95    if (isBad())
  96       {
  97       printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
  98       printf (" +++ Quadrangle impossible \n");
  99       printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
 100       dump ();
 101       printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
 102       // el_root->dump ();
 103       for (int ned=0; ned<QUAD4; ned++)
 104           {
 105           q_edge[ned]->dumpPlus ();
 106           }
 107       HexDump (q_vertex[0]);
 108       HexDump (q_vertex[1]);
 109       HexDump (q_vertex[2]);
 110       HexDump (q_vertex[3]);
 111
 112       printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
 113       fatal_error ("Quadrangle impossible");
 114       }
 115
 116    majReferences ();
 117 }
 118 // ======================================================== Constructeur ter
 119 Quad::Quad (Quad* other)
 120     : EltBase (other->dad(), EL_QUAD)
 121 {
 122    for (int nro=0 ; nro<QUAD4 ; nro++)
 123        {
 124        q_edge   [nro] = NULL;
 125        q_vertex [nro] = NULL;
 126        }
 127    q_orientation = Q_UNDEFINED;
 128    q_clone       = NULL;
 129 }
 130 // ============================================================  getEdge
 131 Edge* Quad::getEdge (int nro)
 132 {
 133    Edge* elt = NULL;
 134    if (nro >=0 && nro < QUAD4 && el_status == HOK && q_edge [nro]->isValid())
 135       elt = q_edge [nro];
 136
 137    DumpStart  ("getEdge", nro);
 138    DumpReturn (elt);
 139    return elt;
 140 }
 141 // ============================================================  getVertex
 142 Vertex* Quad::getVertex (int nro)
 143 {
 144    Vertex* elt = NULL;
 145    if (nro >=0 && nro < QUAD4 && el_status == HOK && q_vertex [nro]->isValid())
 146       elt = q_vertex [nro];
 147
 148    DumpStart  ("getVertex", nro);
 149    DumpReturn (elt);
 150    return elt;
 151 }
 152 // ========================================================= majReferences
 153 void Quad::majReferences ()
 154 {
 155    for (int nro=0 ; nro<QUAD4 ; nro++)
 156        q_edge [nro] -> addParent (this);
 157 }
 158 // ========================================================= getParent
 159 Hexa* Quad::getParent  (int nro)
 160 {
 161    return static_cast <Hexa*> (getFather (nro));
 162 }
 163 // ======================================================== anaMerge
 164 int Quad::anaMerge (Vertex* v1, Vertex* v2, Vertex* tv1[], Edge* te1[])
 165 {
 166    int orig = NOTHING;
 167    for (int nro=0 ; orig == NOTHING && nro < QUAD4 ; nro++)
 168        if (q_vertex [nro] == v1)
 169            orig = nro;
 170
 171    if (orig==NOTHING)
 172       return HERR;
 173
 174    int nsp1 = (orig+1)       MODULO QUAD4;
 175    int nsm1 = (orig+QUAD4-1) MODULO QUAD4;
 176
 177    if (q_vertex [nsp1] == v2)
 178       {
 179       for (int nro=0 ; nro < QUAD4 ; nro++)
 180           {
 181           tv1 [nro] = q_vertex [(orig+nro) MODULO QUAD4];
 182           te1 [nro] = q_edge   [(orig+nro) MODULO QUAD4];
 183           }
 184       }
 185    else if (q_vertex [nsm1] == v2)
 186       {
 187       for (int nro=0 ; nro < QUAD4 ; nro++)
 188           {
 189           tv1 [nro] = q_vertex [(orig+QUAD4-nro) MODULO QUAD4];
 190           te1 [nro] = q_edge   [(orig+QUAD4-nro) MODULO QUAD4];
 191           }
 192       }
 193    else
 194       return 588;
 195
 196    return HOK;
 197 }
 198 // ======================================================== ordoVertex
 199 int Quad::ordoVertex (Vertex* v1, Vertex* v2, Vertex* tv1[])
 200 {
 201    int orig = NOTHING;
 202    for (int nro=0 ; orig == NOTHING && nro < QUAD4 ; nro++)
 203        if (q_vertex [nro] == v1)
 204            orig = nro;
 205
 206    if (orig==NOTHING)
 207       return HERR;
 208
 209    int nsp1 = (orig+1)       MODULO QUAD4;
 210    int nsm1 = (orig+QUAD4-1) MODULO QUAD4;
 211
 212    if (q_vertex [nsp1] == v2)
 213       {
 214       for (int nro=0 ; nro < QUAD4 ; nro++)
 215           tv1 [nro] = q_vertex [(orig+nro) MODULO QUAD4];
 216       }
 217    else if (q_vertex [nsm1] == v2)
 218       {
 219       for (int nro=0 ; nro < QUAD4 ; nro++)
 220           tv1 [nro] = q_vertex [(orig+QUAD4-nro) MODULO QUAD4];
 221       }
 222    else
 223       return 588;
 224
 225    return HOK;
 226 }
 227 // ======================================================== getBrother
 228 Quad* Quad::getBrother (StrOrient* orient)
 229 {
 230 /* *****************************
 231    printf (" getBrother ");
 232    dump ();
 233    printf (" .. Base  : ");
 234    orient->v21->printName();
 235    orient->v22->printName();
 236    printf ("dir=%d, arete=", orient->dir);
 237   ***************************** */
 238
 239    int n21 = indexVertex (orient->v21);
 240    int n22 = indexVertex (orient->v22);
 241
 242    int sens  = n22 - n21;
 243    if (sens >  1) sens -= QUAD4;
 244    if (sens < -1) sens += QUAD4;
 245    if (sens*sens !=1) return NULL;
 246
 247    switch (orient->dir)
 248       {
 249       case OR_LEFT  : n22 = n21 - sens;
 250            break;
 251       case OR_RIGHT : n21 = n22 + sens;
 252            break;
 253       case OR_FRONT : n21 += 2;
 254                       n22 += 2;
 255            break;
 256       default : ;
 257       }
 258
 259    n21 = (n21 + QUAD4) MODULO QUAD4;
 260    n22 = (n22 + QUAD4) MODULO QUAD4;
 261
 262    orient->v21 = q_vertex [n21];
 263    orient->v22 = q_vertex [n22];
 264
 265    Edge* arete  = findEdge (orient->v21, orient->v22);
 266    // arete->printName("\n");
 267
 268    int nbfreres = arete->getNbrParents ();
 269
 270    for (int nq = 0 ; nq < nbfreres ; nq++)
 271        {
 272        Quad* next = arete->getParent (nq);
 273        if (next!=NULL && next != this )
 274           {
 275           int   nbp   = next->getNbrParents();
 276           Hexa* dad   = next->getParent(0);
 277           int   mark  = next->getMark();
 278           int   mark2 = dad ? dad->getMark() : IS_NONE;
 279
 280           if (nbp  <= 1  && mark2 != IS_MARRIED && mark == IS_NONE)
 281              return next;
 282           // if (nbp  <= 1  && mark == IS_NONE)
 283              // return next;
 284           }
 285        }
 286    return NULL;
 287 }
 288 // ======================================================== coupler
 289 int Quad::coupler (Quad* other, StrOrient* orient, Elements* table)
 290 {
 291    if (other==NULL)
 292       return HERR;
 293
 294    Hexa* hexa = other->getParent(0);
 295
 296    setMark (IS_MARRIED);
 297    other->setMark (IS_MARRIED);
 298    if (hexa != NULL)
 299        hexa->setMark (IS_MARRIED);
 300
 301    for (int ned = 0 ; ned < QUAD4 ; ned++)
 302        {
 303        Edge* arete  = q_edge[ned];
 304        int nbfreres = arete ->getNbrParents ();
 305        for (int nq = 0 ; nq < nbfreres ; nq++)
 306            {
 307            Quad* next = arete->getParent (nq);
 308            if (next!=NULL && next != this && next->getMark() > 0)
 309               {
 310               StrOrient new_ori (orient);
 311               new_ori.dir = OR_FRONT;
 312               Vertex* va  = arete->getVertex (V_AMONT);
 313               Vertex* vb  = arete->getVertex (V_AVAL);
 314
 315 //    On voit si un point de repere est conserve
 316               if (va == orient->v11)
 317                  {
 318                  new_ori.v12 = vb;
 319                  new_ori.dir += OR_LEFT;
 320                  }
 321               else if (vb == orient->v11)
 322                  {
 323                  new_ori.v12 = va;
 324                  new_ori.dir += OR_LEFT;
 325                  }
 326
 327               if (va == orient->v12)
 328                  {
 329                  new_ori.v11 = vb;
 330                  new_ori.dir += OR_RIGHT;
 331                  }
 332               else if (vb == orient->v12)
 333                  {
 334                  new_ori.v11 = va;
 335                  new_ori.dir += OR_RIGHT;
 336                  }
 337
 338               if (new_ori.dir == OR_FRONT)
 339                  {
 340                  if (definedBy (va, orient->v11))
 341                     {
 342                     new_ori.v11 = va;
 343                     new_ori.v12 = vb;
 344                     }
 345                  else
 346                     {
 347                     new_ori.v11 = vb;
 348                     new_ori.v12 = va;
 349                     }
 350                  }
 351
 352               int nro = next->getMark ();
 353               Quad* beauf = other->getBrother (&new_ori);
 354               int ier = table->coupler (nro, beauf, &new_ori);
 355               if (ier != HOK)
 356                  return ier;
 357               ier = next->coupler (beauf, &new_ori, table);
 358               if (ier != HOK)
 359                  return ier;
 360               }
 361            }
 362        }
 363    return HOK;
 364 }
 365 // ======================================================== getOpposVertex
 366 Vertex* Quad::getOpposVertex (Vertex* start)
 367 {
 368    int  na = indexVertex (start);
 369    if (na==NOTHING)
 370       return NULL;
 371    return  q_vertex [(na+2) MODULO QUAD4];
 372 }
 373 // ======================================================== getOpposEdge
 374 Edge* Quad::getOpposEdge (Edge* start, int& sens)
 375 {
 376    sens = 1;
 377    int  na = indexVertex (start->getVertex (V_AMONT));
 378    int  nb = indexVertex (start->getVertex (V_AVAL));
 379
 380    Vertex* vaprim = q_vertex [(nb+2) MODULO QUAD4];
 381    Vertex* vbprim = q_vertex [(na+2) MODULO QUAD4];
 382
 383    for (int ned = 0 ; ned < QUAD4 ; ned++)
 384        {
 385        if (   q_edge[ned]->getVertex(V_AMONT) == vaprim
 386            && q_edge[ned]->getVertex(V_AVAL ) == vbprim)
 387            {
 388            sens = 1;
 389            return q_edge[ned];
 390            }
 391        else if (   q_edge[ned]->getVertex(V_AMONT) == vbprim
 392                 && q_edge[ned]->getVertex(V_AVAL ) == vaprim)
 393            {
 394            sens = -1;
 395            return q_edge[ned];
 396            }
 397        }
 398    //             TODO : traiter l'erreur
 399    cout << " ... Probleme dans Quad::getOpposedEdge :" << endl;
 400    HexDisplay (el_name);
 401    PutName (start);
 402    PutName (vaprim);
 403    PutName (vbprim);
 404    HexDisplay (na);
 405    HexDisplay (nb);
 406    dumpPlus ();
 407
 408    for (int ned = 0 ; ned < QUAD4 ; ned++)
 409        q_edge[ned]->dump();
 410
 411    return NULL;
 412 }
 413 // ========================================================= saveXml
 414 void Quad::saveXml (XmlWriter* xml)
 415 {
 416    char buffer[12];
 417    string edges;
 418
 419    for (int nro=0 ; nro<QUAD4 ; nro++)
 420        {
 421        if (nro>0) edges += " ";
 422        edges += q_edge[nro]->getName(buffer);
 423        }
 424
 425    xml->openMark     ("Quad");
 426    xml->addAttribute ("id",    getName (buffer));
 427    xml->addAttribute ("edges", edges);
 428    if (el_name!=buffer)
 429        xml->addAttribute ("name", el_name);
 430    xml->closeMark ();
 431
 432    int nbass = tab_assoc.size();
 433    for (int nro=0 ; nro<nbass ; nro++)
 434        if (tab_assoc[nro] != NULL)
 435            tab_assoc[nro]->saveXml (xml);
 436 }
 437 // ======================================================== replaceEdge
 438 void Quad::replaceEdge (Edge* old, Edge* par)
 439 {
 440    for (int nro=0 ; nro<QUAD4 ; nro++)
 441        {
 442        if (q_edge[nro]==old)
 443            {
 444            q_edge[nro] = par;
 445            if (debug())
 446               {
 447               printf (" Dans ");
 448               printName ();
 449               printf (" [%d], ", nro);
 450               old->printName (" est remplace par ");
 451               par->printName ("\n");
 452               }
 453            }
 454        }
 455 }
 456 // ======================================================== replaceVertex
 457 void Quad::replaceVertex (Vertex* old, Vertex* par)
 458 {
 459    for (int nro=0 ; nro<QUAD4 ; nro++)
 460        {
 461        if (q_vertex [nro]==old)
 462           {
 463           q_vertex [nro] = par;
 464            if (debug())
 465               {
 466               printf (" Dans ");
 467               printName ();
 468               printf (" [%d], ", nro);
 469               old->printName (" est remplace par ");
 470               par->printName ("\n");
 471               }
 472           }
 473        }
 474 }
 475 // ======================================================== dump
 476 void Quad::dump ()
 477 {
 478    printName(" = (");
 479    if (NOT isHere ())
 480       {
 481       printf ("*** deleted ***)\n");
 482       return;
 483       }
 484
 485    for (int nro=0 ; nro<QUAD4 ; nro++)
 486         PrintName (q_edge[nro]);
 487    printf (")\n");
 488
 489    printf ("        (");
 490    for (int nro=0 ; nro<QUAD4 ; nro++)
 491         PrintName (q_vertex[nro]);
 492    printf (")");
 493
 494    dumpRef ();
 495 }
 496 // ======================================================== dumpPlus
 497 void Quad::dumpPlus ()
 498 {
 499    dump ();
 500    if (NOT isHere ())
 501       return;
 502
 503    for (int nro=0 ; nro < QUAD4 ; nro++)
 504        {
 505        Vertex* pv = q_vertex[nro];
 506        printf ( "    ");
 507        if (pv!=NULL)
 508           {
 509           pv->printName ("");
 510           printf ( " (%g, %g, %g)\n", pv->getX(),  pv->getY(),  pv->getZ());
 511           }
 512        else
 513           {
 514           printf ( "NULL\n");
 515           }
 516        }
 517 }
 518 // ======================================================== getOpposEdge (2)
 519 Edge* Quad::getOpposEdge (Edge* start)
 520 {
 521    int  na = indexEdge (start);
 522    if (na<0)
 523       return NULL;
 524    return q_edge [(na+2) MODULO QUAD4];
 525 }
 526 // ======================================================== getPerpendicular
 527 Edge* Quad::getPerpendicular (Edge* arete, Vertex* node)
 528 {
 529    int na = indexEdge (arete);
 530    if (na<0)
 531       return NULL;
 532
 533    int nv = arete->index (node);
 534    if (nv<0)
 535       return NULL;
 536    Edge* perp = q_edge [(na+1) MODULO QUAD4];
 537
 538    nv = perp->index (node);
 539    if (nv>=0)
 540       return perp;
 541
 542    perp = q_edge [(na+3) MODULO QUAD4];
 543    nv   = perp->index (node);
 544    if (nv>=0)
 545       return perp;
 546    else
 547       return NULL;
 548 }
 549 static cpchar t_ori[] = {"Q_INSIDE", "Q_DIRECT", "Q_INVERSE", "Q_UNDEF"};
 550 // ======================================================== setOrientation
 551 void Quad::setOrientation (int ori)
 552 {
 553     q_orientation = ori;
 554     if (db && (ori==Q_DIRECT || ori==Q_INVERSE))
 555        printf (" %s = %s\n", el_name.c_str(), t_ori [ q_orientation ]);
 556 }
 557 // ======================================================== setOrientation
 558 int Quad::setOrientation ()
 559 {
 560     q_orientation = Q_INSIDE;
 561     if (getNbrParents() != 1)
 562        return q_orientation;
 563
 564     Real3 cg, orig, pi, pj, vi, vj, vk;
 565
 566     Hexa* hexa = getParent(0);
 567     hexa->getCenter (cg);
 568
 569 /********************************************************************
 570     printName (" = ");
 571
 572     for (int np=0 ; np < QUAD4 ; np++)
 573         {
 574         q_vertex [np        ] -> getPoint (orig);
 575         q_vertex [(np+1) % 4] -> getPoint (pi);
 576         q_vertex [(np+3) % 4] -> getPoint (pj);
 577
 578         calc_vecteur (orig, pi, vi);
 579         calc_vecteur (orig, pj, vj);
 580         calc_vecteur (orig, cg, vk);
 581         double pmixte = prod_mixte (vi, vj, vk);
 582         q_orientation = pmixte > ZEROR ? Q_DIRECT : Q_INVERSE;
 583         if (pmixte>0) printf (">");
 584            else       printf ("<");
 585         }
 586
 587     printf ("\n");
 588     return;
 589   ******************************************************************* */
 590     q_vertex [0] -> getPoint (orig);
 591     q_vertex [1] -> getPoint (pi);
 592     q_vertex [3] -> getPoint (pj);
 593
 594     calc_vecteur (orig, pi, vi);
 595     calc_vecteur (orig, pj, vj);
 596     calc_vecteur (cg, orig, vk);
 597
 598     double pmixte = prod_mixte (vi, vj, vk);
 599     q_orientation = pmixte > ZEROR ? Q_DIRECT : Q_INVERSE;
 600     if (db)
 601        printf (" %s = %s\n", el_name.c_str(), t_ori [ q_orientation ]);
 602     return q_orientation;
 603 }
 604 // ========================================================== clearAssociation
 605 void Quad::clearAssociation ()
 606 {
 607    tab_assoc.clear ();
 608    is_associated = false;
 609 }
 610 // ========================================================== addAssociation
 611 int Quad::addAssociation (NewShape* geom, int subid)
 612 {
 613    if (geom == NULL)
 614       return HERR;
 615
 616    FaceShape* face = geom->findFace (subid);
 617    int ier = addAssociation (face);
 618
 619    return ier;
 620 }
 621 // ========================================================== addAssociation
 622 int Quad::addAssociation (FaceShape* face)
 623 {
 624    if (face == NULL)
 625       return HERR;
 626
 627    face->addAssociation (this);
 628    tab_assoc.push_back (face);
 629    is_associated = true;
 630    return HOK;
 631 }
 632 // ========================================================== getAssociation
 633 FaceShape* Quad::getAssociation (int nro)
 634 {
 635    if (nro < 0 || nro >= tab_assoc.size())
 636       return NULL;
 637
 638    return tab_assoc [nro];
 639 }
 640 END_NAMESPACE_HEXA