src/HEXABLOCK/HexQuad.cxx

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