src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx

   1 // Copyright (C) 2009-2016  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 //  File   : SMESH_HexaBlocks.cxx
  20 //  Author :
  21 //  Module : SMESH
  22 //
  23
  24 #include <sstream>
  25 #include <algorithm>
  26
  27 // CasCade includes
  28
  29 #include <Precision.hxx>
  30 #include <BRep_Tool.hxx>
  31 #include <BRepTools.hxx>
  32 #include <BRep_Builder.hxx>
  33 #include <BRepAdaptor_Curve.hxx>
  34 #include <BRepBuilderAPI_MakeFace.hxx>
  35
  36 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
  37 #include <GCPnts_AbscissaPoint.hxx>
  38 #include <TopoDS_Wire.hxx>
  39
  40 #include <TopoDS.hxx>
  41 #include <TopoDS_Shape.hxx>
  42 #include <TopoDS_Compound.hxx>
  43 #include <gp_Pln.hxx>
  44 #include <gp_Pnt.hxx>
  45 #include <gp_Dir.hxx>
  46 #include <gp_Lin.hxx>
  47 #include <IntCurvesFace_ShapeIntersector.hxx>
  48
  49 // SMESH includes
  50 #include "SMDS_MeshNode.hxx"
  51 #include "SMDS_MeshVolume.hxx"
  52 #include "SMESH_Gen.hxx"
  53 #include "SMESH_MesherHelper.hxx"
  54 #include "SMESHDS_Group.hxx"
  55
  56 // HEXABLOCK includes
  57 #include "HexDocument.hxx"
  58 #include "HexVertex.hxx"
  59 #include "HexEdge.hxx"
  60 #include "HexQuad.hxx"
  61 #include "HexHexa.hxx"
  62 #include "HexPropagation.hxx"
  63 #include "HexGroup.hxx"
  64 #include "hexa_base.hxx"
  65 #include "HexAssoEdge.hxx"
  66
  67 // HEXABLOCKPLUGIN includes
  68 #include "HEXABLOCKPlugin_mesh.hxx"
  69 #include "HEXABLOCKPlugin_FromSkin_3D.hxx"
  70
  71 // other includes
  72 #include "Basics_Utils.hxx"
  73 #include "utilities.h"
  74
  75 #ifdef WNT
  76 #include <process.h>
  77 #else
  78 #include <unistd.h>
  79 #endif
  80
  81 #include <stdexcept>
  82
  83 #ifndef EXCEPTION
  84 #define EXCEPTION(TYPE, MSG) {\
  85   std::ostringstream aStream;\
  86   aStream<<__FILE__<<"["<<__LINE__<<"]::"<<MSG;\
  87   throw TYPE(aStream.str());\
  88 }
  89 #endif
  90
  91 #ifdef _DEBUG_
  92 static int  MYDEBUG = HEXA_NS::on_debug ();
  93 #else
  94 static int  MYDEBUG = 0;
  95 #endif
  96
  97
  98 static double HEXA_EPS      = 1.0e-3; //1E-3;
  99 static double HEXA_QUAD_WAY = M_PI/4.; //3.*PI/8.;
 100
 101 // ============================================================ string2shape
 102 TopoDS_Shape string2shape( const string& brep )
 103 {
 104   TopoDS_Shape shape;
 105   istringstream streamBrep(brep);
 106   BRep_Builder aBuilder;
 107   BRepTools::Read(shape, streamBrep, aBuilder);
 108   return shape;
 109 }
 110
 111 // ============================================================ shape2coord
 112 bool shape2coord(const TopoDS_Shape& aShape, double& x, double& y, double& z)
 113 {
 114    if ( aShape.ShapeType() == TopAbs_VERTEX ){
 115       TopoDS_Vertex aPoint;
 116        aPoint = TopoDS::Vertex( aShape );
 117        gp_Pnt aPnt = BRep_Tool::Pnt( aPoint );
 118        x = aPnt.X();
 119        y = aPnt.Y();
 120        z = aPnt.Z();
 121        return(1);
 122    } else {
 123        return(0);
 124    };
 125 }
 126 // ============================================================= edge_length
 127 double edge_length(const TopoDS_Edge & E)
 128 {
 129   double UMin = 0, UMax = 0;
 130   if (BRep_Tool::Degenerated(E))
 131     return 0;
 132   TopLoc_Location L;
 133   Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, UMin, UMax);
 134   GeomAdaptor_Curve AdaptCurve(C);
 135   double length = GCPnts_AbscissaPoint::Length(AdaptCurve, UMin, UMax);
 136   return length;
 137 }
 138 // =============================================================== make_curve
 139 BRepAdaptor_Curve* make_curve (gp_Pnt& pstart, gp_Pnt& pend,
 140                                double& length, double& u_start,
 141                                HEXA_NS::Edge& edge, int nro)
 142 {
 143    Hex::AssoEdge*     asso       = edge.getAssociation (nro);
 144    BRepAdaptor_Curve* the_curve  = asso->getCurve();
 145    length    = asso->length ();
 146    if (length <= 1e-6)
 147        return NULL;
 148
 149    const double* point1 = asso->getOrigin ();
 150    const double* point2 = asso->getExtrem ();
 151    pstart  = gp_Pnt (point1[0],  point1[1],  point1[2]);
 152    pend    = gp_Pnt (point2[0],  point2[1],  point2[2]);
 153    u_start = asso->getUstart();
 154    return the_curve;
 155 }
 156 // ============================================================== Constructeur
 157 // SMESH_HexaBlocks::SMESH_HexaBlocks( SMESH_Mesh* theMesh ):
 158 SMESH_HexaBlocks::SMESH_HexaBlocks(SMESH_Mesh& theMesh):
 159   _total(0),
 160   _found(0),
 161   _notFound(0),
 162   _computeVertexOK(false),
 163   _computeEdgeOK(false),
 164   _computeQuadOK(false),
 165   _theMesh(&theMesh),  //groups creation
 166   _theMeshDS(theMesh.GetMeshDS()) //meshing
 167 {
 168 }
 169
 170
 171 SMESH_HexaBlocks::~SMESH_HexaBlocks()
 172 {
 173 }
 174
 175
 176 // --------------------------------------------------------------
 177 //                      PUBLIC METHODS
 178 // --------------------------------------------------------------
 179
 180 // --------------------------------------------------------------
 181 //                     Vertex computing
 182 // --------------------------------------------------------------
 183 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
 184 // ============================================================== computeVertex
 185 bool SMESH_HexaBlocks::computeVertex(HEXA_NS::Vertex& vx)
 186 {
 187   double           px, py, pz;
 188   vx.getAssoCoord (px, py, pz);
 189
 190   SMDS_MeshNode* new_node = _theMeshDS->AddNode (px, py, pz);
 191   _vertex [new_node] = &vx;
 192   _node   [&vx]      = new_node;    //needed in computeEdge()
 193   _computeVertexOK   = true;
 194   return true;
 195 }
 196 // --------------------------------------------------------------
 197 //                      Edge computing
 198 // --------------------------------------------------------------
 199 bool SMESH_HexaBlocks::computeEdge(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
 200 {
 201   bool ok = false;
 202
 203   ok = computeEdgeByAssoc( edge, law);
 204   if ( ok == false ){
 205     ok = computeEdgeByShortestWire( edge, law);
 206   }
 207   if ( ok == false ){
 208     ok = computeEdgeByPlanWire( edge, law);
 209   }
 210   if ( ok == false ){
 211     ok = computeEdgeByIsoWire( edge, law);
 212   }
 213   if ( ok == false ){
 214     ok = computeEdgeBySegment( edge, law);
 215   }
 216   if (ok == true){
 217     _computeEdgeOK = true;
 218   }
 219   return ok;
 220 }
 221
 222
 223
 224 bool SMESH_HexaBlocks::computeEdgeByAssoc( HEXA_NS::Edge& edge, HEXA_NS::Law& law )
 225 {
 226   if(MYDEBUG) MESSAGE("computeEdgeByAssoc(edgeID = "<<edge.getId()<<"): begin   <<<<<<");
 227   ASSERT( _computeVertexOK );
 228   bool ok = true;
 229
 230   if (NOT edge.isAssociated())
 231      return false;
 232
 233   //vertex from edge
 234   HEXA_NS::Vertex* vx0 = NULL;
 235   HEXA_NS::Vertex* vx1 = NULL;
 236
 237   // way of discretization
 238   if (edge.getWay() == true){
 239     vx0 = edge.getVertex(0);
 240     vx1 = edge.getVertex(1);
 241   } else {
 242     vx0 = edge.getVertex(1);
 243     vx1 = edge.getVertex(0);
 244   }
 245   // nodes on mesh
 246   SMDS_MeshNode* FIRST_NODE = _node[vx0];
 247   SMDS_MeshNode* LAST_NODE  = _node[vx1];
 248
 249
 250   // A) Build myCurve_list
 251   std::list< BRepAdaptor_Curve* >        myCurve_list;
 252   double                                 myCurve_tot_len;
 253   std::map< BRepAdaptor_Curve*, double>  myCurve_lengths;
 254   std::map< BRepAdaptor_Curve*, bool>    myCurve_ways;
 255   std::map< BRepAdaptor_Curve*, double>  myCurve_starts;
 256
 257   gp_Pnt myCurve_pt_start(FIRST_NODE->X(), FIRST_NODE->Y(), FIRST_NODE->Z());
 258   gp_Pnt myCurve_pt_end  (LAST_NODE->X(),  LAST_NODE->Y(),  LAST_NODE->Z());
 259
 260   _buildMyCurve(
 261       myCurve_pt_start,
 262       myCurve_pt_end,
 263       myCurve_list,
 264       myCurve_tot_len,
 265       myCurve_lengths,
 266       myCurve_ways,
 267       myCurve_starts,
 268       edge
 269   );
 270
 271
 272   // B) Build nodes and edges on mesh from myCurve_list
 273   SMDS_MeshNode* node_a  = NULL;
 274   SMDS_MeshNode* node_b  = NULL;
 275   SMDS_MeshEdge* edge_ab = NULL;
 276   SMESHNodes     nodesOnEdge;
 277   SMESHEdges     edgesOnEdge; //backup for group creation
 278 //   Xx             nodesXxOnEdge;
 279
 280   node_a = FIRST_NODE;
 281   nodesOnEdge.push_back(FIRST_NODE);
 282 //   nodesXxOnEdge.push_back(0.);
 283   // _nodeXx[FIRST_NODE] = 0.;
 284
 285   gp_Pnt ptOnMyCurve;
 286   double u, myCurve_u;
 287   double myCurve_start_u = 0.;
 288   int nbNodes = law.getNodes(); //law of discretization
 289   if (myCurve_list.size()==0)
 290      {
 291      PutData (edge.getName());
 292      nbNodes = 0;
 293      }
 294   if (MYDEBUG) MESSAGE("nbNodes -> "<<nbNodes);
 295   for (int i = 0; i < nbNodes; ++i){
 296       u = _Xx(i, law, nbNodes); //u between [0,1]
 297       myCurve_u = u*myCurve_tot_len;
 298       if (MYDEBUG) {
 299         MESSAGE("u -> "<<u);
 300         MESSAGE("myCurve_u  -> "<<myCurve_u);
 301         MESSAGE("myCurve_tot_len -> "<<myCurve_tot_len);
 302       }
 303       ptOnMyCurve = _getPtOnMyCurve( myCurve_u,
 304                                      myCurve_ways,
 305                                      myCurve_lengths,
 306                                      myCurve_starts,
 307                                      myCurve_list,
 308                                      myCurve_start_u
 309                                      );
 310
 311       node_b = _theMeshDS->AddNode( ptOnMyCurve.X(), ptOnMyCurve.Y(), ptOnMyCurve.Z() );
 312       edge_ab     = _theMeshDS->AddEdge( node_a, node_b );
 313       nodesOnEdge.push_back( node_b );
 314       edgesOnEdge.push_back( edge_ab );
 315       if  (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
 316       _nodeXx[node_b] = u;
 317       node_a = node_b;
 318   }
 319   edge_ab      = _theMeshDS->AddEdge( node_a, LAST_NODE );
 320   nodesOnEdge.push_back( LAST_NODE );
 321   edgesOnEdge.push_back( edge_ab );
 322   _nodesOnEdge[&edge] = nodesOnEdge;
 323   _edgesOnEdge[&edge] = edgesOnEdge;
 324
 325
 326   if(MYDEBUG) MESSAGE("computeEdgeByAssoc() : end  >>>>>>>>");
 327   return ok;
 328 }
 329
 330
 331
 332
 333 bool SMESH_HexaBlocks::computeEdgeByShortestWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
 334 {
 335   if(MYDEBUG) MESSAGE("computeEdgeByShortestWire() not implemented");
 336   return false;
 337 }
 338
 339 bool SMESH_HexaBlocks::computeEdgeByPlanWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
 340 {
 341   if(MYDEBUG) MESSAGE("computeEdgeByPlanWire() not implemented");
 342   return false;
 343 }
 344
 345 bool SMESH_HexaBlocks::computeEdgeByIsoWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
 346 {
 347   if(MYDEBUG) MESSAGE("computeEdgeByIsoWire() not implemented");
 348   return false;
 349 }
 350
 351
 352 bool SMESH_HexaBlocks::computeEdgeBySegment(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
 353 {
 354   if(MYDEBUG) MESSAGE("computeEdgeBySegment() : : begin   <<<<<<");
 355   ASSERT( _computeVertexOK );
 356   bool ok = true;
 357
 358   //vertex from edge
 359   HEXA_NS::Vertex* vx0 = NULL;
 360   HEXA_NS::Vertex* vx1 = NULL;
 361
 362   // way of discretization
 363   if (edge.getWay() == true){
 364     vx0 = edge.getVertex(0);
 365     vx1 = edge.getVertex(1);
 366   } else {
 367     vx0 = edge.getVertex(1);
 368     vx1 = edge.getVertex(0);
 369   }
 370
 371   // nodes on mesh
 372   SMDS_MeshNode* FIRST_NODE = _node[vx0];
 373   SMDS_MeshNode* LAST_NODE  = _node[vx1];
 374   SMDS_MeshNode* node_a = NULL; //new node (to be added)
 375   SMDS_MeshNode* node_b = NULL; //new node (to be added)
 376
 377   // node and edge creation
 378   SMESHNodes nodesOnEdge;
 379   SMESHEdges edgesOnEdge;
 380
 381   double u; //
 382   double newNodeX, newNodeY, newNodeZ;
 383   SMDS_MeshEdge* newEdge = NULL;
 384
 385   node_a = FIRST_NODE;
 386   nodesOnEdge.push_back(FIRST_NODE);
 387
 388   //law of discretization
 389   int nbNodes = law.getNodes();
 390   if (MYDEBUG) MESSAGE("nbNodes -> "<<nbNodes);
 391   for (int i = 0; i < nbNodes; ++i){
 392     u = _Xx(i, law, nbNodes);
 393     newNodeX = FIRST_NODE->X() + u * ( LAST_NODE->X() - FIRST_NODE->X() );
 394     newNodeY = FIRST_NODE->Y() + u * ( LAST_NODE->Y() - FIRST_NODE->Y() );
 395     newNodeZ = FIRST_NODE->Z() + u * ( LAST_NODE->Z() - FIRST_NODE->Z() );
 396     node_b = _theMeshDS->AddNode(newNodeX, newNodeY, newNodeZ);
 397     newEdge = _theMeshDS->AddEdge(node_a, node_b);
 398     edgesOnEdge.push_back(newEdge);
 399     nodesOnEdge.push_back(node_b);
 400     if  (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
 401     _nodeXx[ node_b ] = u;
 402     if(MYDEBUG) MESSAGE("_nodeXx <-"<<u);
 403     node_a = node_b;
 404   }
 405   newEdge = _theMeshDS->AddEdge(node_a, LAST_NODE);
 406   nodesOnEdge.push_back(LAST_NODE);
 407   edgesOnEdge.push_back(newEdge);
 408
 409   _nodesOnEdge[&edge] = nodesOnEdge;
 410   _edgesOnEdge[&edge] = edgesOnEdge;
 411
 412   if(MYDEBUG) MESSAGE("computeEdgeBySegment() : end  >>>>>>>>");
 413   return ok;
 414 }
 415
 416
 417 // --------------------------------------------------------------
 418 //                        Quad computing
 419 // --------------------------------------------------------------
 420 std::map<HEXA_NS::Quad*, bool>  SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document* doc )
 421 {
 422   std::map<HEXA_NS::Quad*, bool>  quadWays;
 423   std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
 424   std::list<HEXA_NS::Quad*>       skinQuad;
 425   std::list<HEXA_NS::Quad*>       workingQuad;
 426   HEXA_NS::Quad* first_q = NULL;
 427   HEXA_NS::Quad* q = NULL;
 428   HEXA_NS::Edge* e = NULL;
 429   HEXA_NS::Vertex *e_0, *e_1 = NULL;
 430
 431   // FIRST STEP: eliminate free quad + internal quad
 432   int nTotalQuad = doc->countUsedQuad();
 433   for (int i=0; i < nTotalQuad; ++i ){
 434     q = doc->getUsedQuad(i);
 435     switch ( q->getNbrParents() ){ // parent == hexaedron
 436       case 0: case 2: quadWays[q] = true; break;
 437       case 1: skinQuad.push_back(q); break;
 438       default: if ( q->getNbrParents() > 2 ) ASSERT(false);
 439     }
 440   }
 441
 442   // SECOND STEP: setting edges ways
 443   while ( skinQuad.size()>0 ){
 444     if(MYDEBUG) MESSAGE("SEARCHING INITIAL QUAD ..." );
 445     for ( std::list<HEXA_NS::Quad*>::iterator it = skinQuad.begin(); it != skinQuad.end(); it++ ){
 446         _searchInitialQuadWay( *it, e_0, e_1 );
 447         if ( e_0 != NULL && e_1 != NULL ){
 448           q = first_q = *it;
 449           break;
 450         }
 451     }
 452     if ( e_0 == NULL && e_1 == NULL ) ASSERT(false);// should never happened,
 453     if(MYDEBUG) MESSAGE("INITIAL QUAD FOUND!" );
 454     for ( int j=0 ; j < 4 ; ++j ){
 455       e = q->getEdge(j);
 456       if  (    ((e_0 == e->getVertex(0)) && (e_1 == e->getVertex(1)))
 457             || ((e_0 == e->getVertex(1)) && (e_1 == e->getVertex(0))) ){
 458         break;
 459       }
 460     }
 461     if(MYDEBUG) MESSAGE("INITIAL EDGE WAY FOUND!" );
 462
 463     edgeWays[e] = std::make_pair( e_0, e_1 );
 464     workingQuad.push_back(q);
 465
 466     while ( workingQuad.size() > 0 ){
 467         if(MYDEBUG) MESSAGE("COMPUTE QUAD WAY ... ID ="<< q->getId());
 468         HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
 469         int i = 0;
 470         std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
 471         while ( localEdgeWays.size() != 4 ){
 472             HEXA_NS::Edge* e = q->getEdge(i%4);
 473             if ( lastVertex == NULL ){
 474                 if ( edgeWays.count(e) == 1 ){
 475                   if ( q == first_q ){
 476                     localEdgeWays[e] = std::make_pair( edgeWays[e].first, edgeWays[e].second );
 477                   } else {
 478                     localEdgeWays[e] = std::make_pair( edgeWays[e].second, edgeWays[e].first);
 479                   }
 480                   firstVertex = localEdgeWays[e].first;
 481                   lastVertex  = localEdgeWays[e].second;
 482                 }
 483             } else {
 484               HEXA_NS::Vertex* e_0 = e->getVertex(0);
 485               HEXA_NS::Vertex* e_1 = e->getVertex(1);
 486
 487               if ( lastVertex == e_0 ){
 488                 firstVertex = e_0; lastVertex = e_1;
 489               } else if ( lastVertex == e_1 ){
 490                 firstVertex = e_1; lastVertex = e_0;
 491               } else if ( firstVertex == e_0 ) {
 492                 firstVertex = e_1; lastVertex = e_0;
 493               } else if ( firstVertex == e_1 ) {
 494                 firstVertex = e_0; lastVertex = e_1;
 495               } else {
 496                 ASSERT(false);
 497               }
 498               localEdgeWays[e] = std::make_pair( firstVertex, lastVertex );
 499               if ( edgeWays.count(e) == 0 ){ // keep current value if present otherwise add it
 500                 edgeWays[e] = localEdgeWays[e];
 501               }
 502             }
 503             ++i;
 504         }
 505
 506
 507         //add new working quad
 508         for ( int i=0 ; i < 4; ++i ){
 509             HEXA_NS::Quad* next_q = NULL;
 510             HEXA_NS::Edge* e = q->getEdge(i);
 511             for ( int j=0 ; j < e->getNbrParents(); ++j ){
 512                 next_q = e->getParent(j);
 513                 if ( next_q == q ){
 514                   next_q = NULL;
 515                 }
 516                 int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
 517                 if (fromSkin != 0){
 518                   int fromWorkingQuad = std::count( workingQuad.begin(), workingQuad.end(), next_q );
 519                     if ( fromWorkingQuad == 0 ){
 520                         workingQuad.push_front( next_q );
 521                     }
 522                 }
 523             }
 524         }
 525
 526         // setting quad way
 527         HEXA_NS::Edge* e0 = q->getEdge(0);
 528         HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
 529
 530         if (  e0_0 == localEdgeWays[ e0 ].first ){
 531             quadWays[q] = true;
 532         } else if ( e0_0 == localEdgeWays[ e0 ].second ){
 533             quadWays[q] = false;
 534         } else {
 535           ASSERT(false);
 536         }
 537         workingQuad.remove( q );
 538         skinQuad.remove( q );
 539         q = workingQuad.back();
 540     }
 541   }
 542
 543   return quadWays;
 544 }
 545
 546
 547
 548
 549
 550 // std::map<HEXA_NS::Quad*, bool>  SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document& doc, std::map<HEXA_NS::Quad*, bool>  initQuads )
 551 // {
 552 //   std::map<HEXA_NS::Quad*, bool>  quadWays;
 553 // //   std::map<HEXA_NS::Edge*, bool>  edgeWays;
 554 // //   std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
 555 //   std::map<HEXA_NS::Quad*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > workingQuads;
 556 //
 557 //   std::list<HEXA_NS::Quad*>       skinQuad;
 558 //   std::list<HEXA_NS::Quad*>       notSkinQuad;
 559 // //   std::list<HEXA_NS::Quad*>       workingQuad;
 560 //   HEXA_NS::Quad* first_q = NULL;
 561 //   HEXA_NS::Quad* q = NULL;
 562 //   HEXA_NS::Edge* e = NULL;
 563 //   HEXA_NS::Vertex *e_0, *e_1 = NULL;
 564 //
 565 //   // FIRST STEP: eliminate free quad + internal quad
 566 //   int nTotalQuad = doc.countQuad();
 567 //   for (int i=0; i < nTotalQuad; ++i ){
 568 //     q = doc.getQuad(i);
 569 //     switch ( q->getNbrParents() ){ // parent == hexaedron
 570 //       case 0: case 2: quadWays[q] = true; break;
 571 // //       case 0: case 2: notSkinQuad.push_back(q); break; //CS_TEST
 572 //       case 1: skinQuad.push_back(q); break;
 573 //       default: if ( q->getNbrParents() > 2 ) ASSERT(false);
 574 //     }
 575 //   }
 576 //
 577 //
 578 //   // SECOND STEP
 579 //   q = first_q = skinQuad.front();
 580 //   e = q->getUsedEdge(0);
 581 //   e_0 = e->getVertex(0);
 582 //   e_1 = e->getVertex(1);
 583 //
 584 //   workingQuads[q] = std::make_pair( e_0, e_1 );
 585 //
 586 //   while ( workingQuads.size() > 0 ){
 587 //       MESSAGE("CURRENT QUAD ------>"<< q->getId());
 588 //       HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
 589 //       int i = 0;
 590 //
 591 //       std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
 592 //       while ( localEdgeWays.size() != 4 ){
 593 //           HEXA_NS::Edge* e = q->getUsedEdge(i%4);
 594 //           if ( lastVertex == NULL ){
 595 //               HEXA_NS::Vertex* e_0 = e->getVertex(0);
 596 //               HEXA_NS::Vertex* e_1 = e->getVertex(1);
 597 //
 598 //               if ( (workingQuads[q].first == e_0 and workingQuads[q].second == e_1)
 599 //                     or (workingQuads[q].first == e_1 and workingQuads[q].second == e_0) ){
 600 //                 if ( q == first_q ){
 601 //                   localEdgeWays[e] = std::make_pair( workingQuads[q].first, workingQuads[q].second );
 602 //                   MESSAGE("FIRST QUAD ");
 603 //                 } else {
 604 //                   localEdgeWays[e] = std::make_pair( workingQuads[q].second, workingQuads[q].first);
 605 //                   MESSAGE("NOT FIRST QUAD ");
 606 //                 }
 607 //                 firstVertex = localEdgeWays[e].first;
 608 //                 lastVertex  = localEdgeWays[e].second;
 609 //               }
 610 //           } else {
 611 //             HEXA_NS::Vertex* e_0 = e->getVertex(0);
 612 //             HEXA_NS::Vertex* e_1 = e->getVertex(1);
 613 //             if ( lastVertex == e_0 ){
 614 //               localEdgeWays[e] = std::make_pair( e_0, e_1 );
 615 //               firstVertex = e_0;
 616 //               lastVertex = e_1;
 617 //             } else if ( lastVertex == e_1 ){
 618 //               localEdgeWays[e] = std::make_pair( e_1, e_0 );
 619 //               firstVertex = e_1;
 620 //               lastVertex = e_0;
 621 //             } else if ( firstVertex == e_0 ) {
 622 //               localEdgeWays[e] = std::make_pair( e_1, e_0 );
 623 //               firstVertex = e_1;
 624 //               lastVertex = e_0;
 625 //             } else if ( firstVertex == e_1 ) {
 626 //               localEdgeWays[e] = std::make_pair( e_0, e_1 );
 627 //               firstVertex = e_0;
 628 //               lastVertex = e_1;
 629 //             } else {
 630 //               ASSERT(false);
 631 //             }
 632 //           }
 633 //           ++i;
 634 //       }
 635 //
 636 //
 637 //       //add new working quad
 638 //       for ( int i=0 ; i < 4; ++i ){
 639 //           HEXA_NS::Quad* next_q = NULL;
 640 //           HEXA_NS::Edge* e = q->getUsedEdge(i);
 641 //           MESSAGE("NB PARENTS ="<< e->getNbrParents() );
 642 //           for ( int j=0 ; j < e->getNbrParents(); ++j ){
 643 //               next_q = e->getParent(j);
 644 //               if ( next_q == q ){
 645 //                 next_q = NULL;
 646 //               }
 647 //               int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
 648 //               if (fromSkin != 0){
 649 // //                 int fromWorkingQuad = std::count( workingQuads.begin(), workingQuads.end(), next_q );
 650 //                 int fromWorkingQuad = workingQuads.count( next_q );
 651 // //             MESSAGE("CHECK QUAD:"<< newWorkingQuad->getId());
 652 //                   if ( fromWorkingQuad == 0 ){
 653 // //                       workingQuads.push_front( next_q );
 654 //                       workingQuads[ next_q ] = localEdgeWays[e];
 655 // //                   MESSAGE("EDGE :"<<e->getId()<<"ADD QUAD :"<< newWorkingQuad->getId());
 656 //                   }
 657 //               }
 658 //           }
 659 //       }
 660 //
 661 //       //setting quad way
 662 //       HEXA_NS::Edge* e0 = q->getUsedEdge(0);
 663 //       HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
 664 //
 665 //       if (  e0_0 == localEdgeWays[ e0 ].first ){
 666 //           quadWays[q] = true;
 667 //       } else if ( e0_0 == localEdgeWays[ e0 ].second ){
 668 //           quadWays[q] = false;
 669 //       } else {
 670 //         ASSERT(false);
 671 //       }
 672 //       MESSAGE("quadWays ID ="<< q->getId() << ", WAY = " << quadWays[q] );
 673 //
 674 // //       workingQuad.remove( q );
 675 //       workingQuads.erase( q );
 676 //       skinQuad.remove( q );
 677 //       *workingQuads.begin();
 678 //       q = (*workingQuads.begin()).first;
 679 //   }
 680 //   return quadWays;
 681 // }
 682
 683
 684 bool SMESH_HexaBlocks::computeQuad( HEXA_NS::Quad& quad, bool way )
 685 {
 686   bool ok = false;
 687
 688   ok = computeQuadByAssoc(quad, way);
 689   if ( ok == false ){
 690     ok = computeQuadByFindingGeom(quad, way);
 691   }
 692   if ( ok == false ){
 693     ok = computeQuadByLinearApproximation(quad, way);
 694   }
 695   if (ok == true){
 696     _computeQuadOK = true;
 697   }
 698   return ok;
 699 }
 700
 701
 702 bool SMESH_HexaBlocks::computeQuadByAssoc( HEXA_NS::Quad& quad, bool way  )
 703 {
 704 //   int id = quad.getId();
 705 //   if ( id != 11 )  return false; //7
 706   if (MYDEBUG){
 707     MESSAGE("computeQuadByLinearApproximation() : : begin   <<<<<<");
 708     MESSAGE("quadID = "<<quad.getId());
 709   }
 710   ASSERT( _computeEdgeOK );
 711   bool ok = true;
 712
 713   ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
 714   SMESHFaces      facesOnQuad;
 715   SMDS_MeshFace*  newFace = NULL;
 716   std::vector<double> xx, yy;
 717
 718   // Elements for quad computation
 719   SMDS_MeshNode *S1, *S2, *S4, *S3;
 720
 721 //   bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
 722   bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
 723   if ( initOk == false ){
 724     return false;
 725   }
 726
 727   int nbass  = quad.countAssociation ();
 728   if (nbass==0)
 729      {
 730      if(MYDEBUG) MESSAGE("computeQuadByAssoc() : end  >>>>>>>>");
 731      return false;
 732      }
 733
 734   TopoDS_Shape shapeOrCompound = getFaceShapes ( quad );
 735
 736
 737   std::map<SMDS_MeshNode*, gp_Pnt> interpolatedPoints;
 738   int iSize = nodesOnQuad.size();
 739   int jSize = nodesOnQuad[0].size();
 740
 741   S1 = nodesOnQuad[0][0];
 742   S2 = nodesOnQuad[iSize-1][0];
 743   S4 = nodesOnQuad[0][jSize-1];
 744   S3 = nodesOnQuad[iSize-1][jSize-1];
 745
 746
 747   for (int j = 1; j < jSize; ++j){
 748     for (int i = 1; i < iSize; ++i){
 749         SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
 750         SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
 751         SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
 752         SMDS_MeshNode* n4 = nodesOnQuad[i][j];
 753
 754         if ( n4 == NULL ){
 755             double newNodeX, newNodeY, newNodeZ;
 756             SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1];   //dNodes[h_i];
 757             SMDS_MeshNode* Pb = nodesOnQuad[i][0];   //bNodes[b_i];
 758             SMDS_MeshNode* Pg = nodesOnQuad[0][j];   //gNodes[g_j];
 759             SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j];  //dNodes[d_j];
 760             double u = xx[i];
 761             double v = yy[j];
 762
 763             _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
 764               gp_Pnt newPt = gp_Pnt( newNodeX, newNodeY, newNodeZ );//interpolated point
 765               gp_Pnt pt1;
 766               gp_Pnt pt3;
 767               if ( interpolatedPoints.count(n1) > 0 ){
 768                       pt1 = interpolatedPoints[n1];
 769               } else {
 770                       pt1 = gp_Pnt( n1->X(), n1->Y(), n1->Z() );
 771               }
 772               if ( interpolatedPoints.count(n3) > 0 ){
 773                       pt3 = interpolatedPoints[n3];
 774               } else {
 775                       pt3 = gp_Pnt( n3->X(), n3->Y(), n3->Z() );
 776               }
 777               gp_Vec vec1( newPt, pt1 );
 778               gp_Vec vec2( newPt, pt3 );
 779
 780               gp_Pnt ptOnShape = _intersect(newPt, vec1, vec2, shapeOrCompound);
 781               newNodeX = ptOnShape.X();
 782               newNodeY = ptOnShape.Y();
 783               newNodeZ = ptOnShape.Z();
 784               n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
 785               nodesOnQuad[i][j] = n4;
 786               interpolatedPoints[ n4 ] = newPt;
 787
 788               if (MYDEBUG) {
 789                   MESSAGE("u parameter is "<<u);
 790                   MESSAGE("v parameter is "<<v);
 791                   MESSAGE("point interpolated ("<<newPt.X()<<","<<newPt.Y()<<","<<newPt.Z()<<" )");
 792                   MESSAGE("point on shape     ("<<newNodeX<<","<<newNodeY<<","<<newNodeZ<<" )");
 793               }
 794         }
 795
 796         if (MYDEBUG){
 797           MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
 798           MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
 799           MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
 800           MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
 801         }
 802
 803         if ( way == true ){
 804             if (MYDEBUG) MESSAGE("AddFace( n1, n2, n3, n4 )");
 805             newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
 806         } else {
 807             if (MYDEBUG) MESSAGE("AddFace( n4, n3, n2, n1 )");
 808             newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
 809         }
 810         facesOnQuad.push_back(newFace);
 811       }
 812   }
 813   _quadNodes[ &quad ] = nodesOnQuad;
 814   _facesOnQuad[&quad] = facesOnQuad;
 815
 816   if(MYDEBUG) MESSAGE("computeQuadByLinearApproximation() : end  >>>>>>>>");
 817   return ok;
 818 }
 819
 820
 821 bool SMESH_HexaBlocks::computeQuadByFindingGeom( HEXA_NS::Quad& quad, bool way )
 822 {
 823   if(MYDEBUG) MESSAGE("computeQuadByFindingGeom() not implemented");
 824   return false;
 825 }
 826
 827 bool SMESH_HexaBlocks::_computeQuadInit(
 828   HEXA_NS::Quad& quad,
 829   ArrayOfSMESHNodes& nodesOnQuad,
 830   std::vector<double>& xx, std::vector<double>& yy)
 831 {
 832   if(MYDEBUG) MESSAGE("_computeQuadInit() : begin ---------------");
 833   bool ok = true;
 834
 835   SMDS_MeshNode *S1, *S2, *S4, *S3;
 836   HEXA_NS::Edge *eh, *eb, *eg, *ed;
 837   HEXA_NS::Edge *e1, *e2, *e3, *e4;
 838   HEXA_NS::Vertex *e1_0, *e1_1, *e2_0, *e2_1, *e3_0, *e3_1, *e4_0, *e4_1;
 839
 840   e1 = quad.getEdge(0);
 841   e2 = quad.getEdge(1);
 842   e3 = quad.getEdge(2);
 843   e4 = quad.getEdge(3);
 844
 845   e1_0 = e1->getVertex(0); e1_1 = e1->getVertex(1);
 846   e2_0 = e2->getVertex(0); e2_1 = e2->getVertex(1);
 847   e3_0 = e3->getVertex(0); e3_1 = e3->getVertex(1);
 848   e4_0 = e4->getVertex(0); e4_1 = e4->getVertex(1);
 849
 850   //S1, S2
 851   S1 = _node[e1_0]; S2 = _node[e1_1];
 852   eb = e1; eh = e3;
 853   //S4
 854   if ( e1_0 == e2_0 ){
 855     S4 = _node[e2_1];
 856     eg = e2; ed = e4;
 857   } else if ( e1_0 == e2_1 ){
 858     S4 = _node[e2_0];
 859     eg = e2; ed = e4;
 860   } else if ( e1_0 == e4_0 ){
 861     S4 = _node[e4_1];
 862     eg = e4; ed = e2;
 863   } else if ( e1_0 == e4_1 ){
 864     S4 = _node[e4_0];
 865     eg = e4; ed = e2;
 866   } else {
 867     ASSERT(false);
 868   }
 869   //S3
 870   if ( S4 == _node[e3_0] ){
 871     S3 = _node[e3_1];
 872   } else if ( S4 == _node[e3_1] ){
 873     S3 = _node[e3_0];
 874   } else {
 875     ASSERT(false);
 876   }
 877
 878   SMESHNodes hNodes = _nodesOnEdge[eh];
 879   SMESHNodes bNodes = _nodesOnEdge[eb];
 880   SMESHNodes gNodes = _nodesOnEdge[eg];
 881   SMESHNodes dNodes = _nodesOnEdge[ed];
 882   nodesOnQuad.resize( bNodes.size(), SMESHNodes(gNodes.size(), static_cast<SMDS_MeshNode*>(NULL)) );
 883
 884
 885   int i, j, _i, _j;
 886 //   int &b_i = i, &h_i = i, &g_j = j, &d_j = j;
 887   int *b_i = &i, *h_i = &i, *g_j = &j, *d_j = &j;
 888   bool uWay = true, vWay = true;
 889
 890   if ( bNodes[0] != S1 ){
 891     b_i = &_i;
 892     uWay = false;
 893     ASSERT( bNodes[bNodes.size()-1] == S1 );
 894   } else {
 895     ASSERT( bNodes[0] == S1);
 896   }
 897   if ( hNodes[0] != S4 ){
 898     h_i = &_i;
 899   } else {
 900     ASSERT( hNodes[0] == S4 );
 901   }
 902   if ( gNodes[0] != S1 ){
 903     g_j = &_j;
 904     vWay = false;
 905   } else {
 906     ASSERT( gNodes[0] == S1 );
 907   }
 908   if ( dNodes[0] != S2 ){
 909     d_j = &_j;
 910   } else {
 911     ASSERT( dNodes[0] == S2 );
 912   }
 913
 914   //bNodes, hNodes
 915   double u;
 916   for (i = 0, _i = bNodes.size()-1; i < bNodes.size(); ++i, --_i){
 917     nodesOnQuad[i][0]                = bNodes[*b_i];
 918     nodesOnQuad[i][gNodes.size()-1 ] = hNodes[*h_i];
 919
 920     u = _nodeXx[ bNodes[*b_i] ];
 921     if ( uWay == true ){
 922       xx.push_back(u);
 923     } else {
 924       xx.push_back(1.-u);
 925     }
 926   }
 927   if ( S1 != nodesOnQuad[0][0] ){
 928     if(MYDEBUG) MESSAGE("ZZZZZZZZZZZZZZZZ quadID = "<<quad.getId());
 929   }
 930 //   ASSERT( S1 == nodesOnQuad[0][0] );
 931
 932   //gNodes, dNodes
 933   double v;
 934   for (j = 0, _j = gNodes.size()-1; j < gNodes.size(); ++j, --_j){
 935     nodesOnQuad[0][j] = gNodes[*g_j];
 936     if ( S1 != nodesOnQuad[0][0] ){
 937       if(MYDEBUG) MESSAGE("XXXXXXXXXXXXXXXX quadID = "<<quad.getId());
 938     }
 939 //     ASSERT( S1 == nodesOnQuad[0][0] );
 940     nodesOnQuad[bNodes.size()-1][j] = dNodes[*d_j];
 941     v = _nodeXx[ gNodes[*g_j] ];
 942     if ( vWay == true ){
 943       yy.push_back(v);
 944     } else {
 945       yy.push_back(1.-v);
 946     }
 947   }
 948
 949
 950   int iSize = nodesOnQuad.size();
 951   int jSize = nodesOnQuad[0].size();
 952   ASSERT( iSize = bNodes.size() );
 953   ASSERT( jSize = gNodes.size() );
 954
 955 //   ASSERT( S1 == nodesOnQuad[0][0] );
 956 //   ASSERT( S2 == nodesOnQuad[iSize-1][0]);
 957 //   ASSERT( S4 == nodesOnQuad[0][jSize-1]);
 958 //   ASSERT( S3 == nodesOnQuad[iSize-1][jSize-1]);
 959
 960   return ok;
 961 }
 962
 963
 964 bool SMESH_HexaBlocks::computeQuadByLinearApproximation( HEXA_NS::Quad& quad, bool way )
 965 {
 966 //   int id = quad.getId();
 967 //   if ( quad.getId() != 66 )  return false; //7, 41
 968   if (MYDEBUG){
 969     MESSAGE("computeQuadByLinearApproximation() : : begin   <<<<<<");
 970     MESSAGE("quadID = "<<quad.getId());
 971   }
 972   ASSERT( _computeEdgeOK );
 973   bool ok = true;
 974
 975   ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
 976   SMESHFaces      facesOnQuad;
 977   SMDS_MeshFace*  newFace = NULL;
 978   std::vector<double> xx, yy;
 979
 980   // Elements for quad computation
 981   SMDS_MeshNode *S1, *S2, *S4, *S3;
 982
 983 //   bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
 984   bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
 985   if ( initOk == false ){
 986     return false;
 987   }
 988
 989   int iSize = nodesOnQuad.size();
 990   int jSize = nodesOnQuad[0].size();
 991
 992   S1 = nodesOnQuad[0][0];
 993 //   S2 = nodesOnQuad[bNodes.size()-1][0];
 994   S2 = nodesOnQuad[iSize-1][0];
 995   S4 = nodesOnQuad[0][jSize-1];
 996   S3 = nodesOnQuad[iSize-1][jSize-1];
 997
 998   for (int j = 1; j < jSize; ++j){
 999     for (int i = 1; i < iSize; ++i){
1000         SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
1001         SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
1002         SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
1003         SMDS_MeshNode* n4 = nodesOnQuad[i][j];
1004
1005         if ( n4 == NULL ){
1006             double newNodeX, newNodeY, newNodeZ;
1007             SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1];   //dNodes[h_i];
1008             SMDS_MeshNode* Pb = nodesOnQuad[i][0];   //bNodes[b_i];
1009             SMDS_MeshNode* Pg = nodesOnQuad[0][j];   //gNodes[g_j];
1010             SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j];  //dNodes[d_j];
1011             double u = xx[i];
1012             double v = yy[j];
1013
1014             _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
1015             n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
1016             nodesOnQuad[i][j] = n4;
1017         }
1018
1019         if (MYDEBUG){
1020           MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
1021           MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
1022           MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
1023           MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
1024         }
1025
1026         if ( way == true ){
1027             if (MYDEBUG) MESSAGE("AddFace( n1, n2, n3, n4 )");
1028             newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
1029         } else {
1030             if (MYDEBUG) MESSAGE("AddFace( n4, n3, n2, n1 )");
1031             newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
1032         }
1033         facesOnQuad.push_back(newFace);
1034       }
1035   }
1036   _quadNodes[ &quad ] = nodesOnQuad;
1037   _facesOnQuad[&quad] = facesOnQuad;
1038
1039   if(MYDEBUG) MESSAGE("computeQuadByLinearApproximation() : end  >>>>>>>>");
1040   return ok;
1041 }
1042
1043
1044 // --------------------------------------------------------------
1045 //                      Hexa computing
1046 // --------------------------------------------------------------
1047 bool SMESH_HexaBlocks::computeHexa( HEXA_NS::Document* doc )
1048 {
1049   if(MYDEBUG) MESSAGE("computeHexa() : : begin   <<<<<<");
1050   bool ok=false;
1051
1052   SMESH_MesherHelper aHelper(*_theMesh);
1053   TopoDS_Shape shape = _theMesh->GetShapeToMesh();
1054   aHelper.SetSubShape( shape );
1055   aHelper.SetElementsOnShape( true );
1056
1057   SMESH_Gen* gen = _theMesh->GetGen();
1058   SMESH_HexaFromSkin_3D algo( gen->GetANewId(), 0, gen, doc );
1059   algo.InitComputeError();
1060   try {
1061       ok = algo.Compute( *_theMesh, &aHelper, _volumesOnHexa, _node );
1062   } catch(...) {
1063     if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D error!!! ");
1064   }
1065   if (MYDEBUG){
1066     MESSAGE("SMESH_HexaFromSkin_3D.comment = "<<algo.GetComputeError()->myComment);
1067     MESSAGE("computeHexa() : end  >>>>>>>>");
1068   }
1069   return ok;
1070 }
1071
1072
1073
1074 // --------------------------------------------------------------
1075 //                Document computing
1076 // --------------------------------------------------------------
1077 bool SMESH_HexaBlocks::computeDoc(  HEXA_NS::Document* doc )
1078 {
1079   if(MYDEBUG) MESSAGE("computeDoc() : : begin   <<<<<<");
1080   bool ok = true;
1081
1082   // A) Vertex computation
1083
1084   doc->lockDump ();
1085   int nVertex = doc->countUsedVertex();
1086   HEXA_NS::Vertex* vertex = NULL;
1087
1088   for (int j=0; j <nVertex; ++j ){ //Computing each vertex of the document
1089     vertex = doc->getUsedVertex(j);
1090     ok = computeVertex(*vertex);
1091   }
1092
1093   // B) Edges computation
1094   int nbPropa = 0;
1095   HEXA_NS::Propagation* propa = NULL;
1096   HEXA_NS::Law*         law   = NULL;
1097   HEXA_NS::Edges edges;
1098
1099   nbPropa = doc->countPropagation();
1100   for (int j=0; j < nbPropa; ++j ){//Computing each edge's propagations of the document
1101     propa = doc->getPropagation(j);
1102     edges = propa->getEdges();
1103     law   = propa->getLaw();
1104 //     ASSERT( law );
1105     if (law == NULL){
1106       law = doc->getLaw(0); // default law
1107     }
1108     for( HEXA_NS::Edges::const_iterator iter = edges.begin();
1109         iter != edges.end();
1110         ++iter ){
1111         ok = computeEdge(**iter, *law);
1112     }
1113   }
1114   // C) Quad computation
1115   std::map<HEXA_NS::Quad*, bool>  quadWays = computeQuadWays(doc);
1116   int nQuad = doc->countUsedQuad();
1117   HEXA_NS::Quad* q = NULL;
1118   for (int j=0; j <nQuad; ++j ){ //Computing each quad of the document
1119     q = doc->getUsedQuad(j);
1120     int id = q->getId();
1121     if ( quadWays.count(q) > 0 )
1122       ok = computeQuad( *q, quadWays[q] );
1123     else
1124       if(MYDEBUG) MESSAGE("NO QUAD WAY ID = "<<id);
1125
1126   }
1127
1128   // D) Hexa computation: Calling HexaFromSkin algo
1129   ok = computeHexa(doc);
1130
1131   if(MYDEBUG) MESSAGE("computeDoc() : end  >>>>>>>>");
1132   doc->lockDump ();
1133   return ok;
1134 }
1135
1136
1137 void SMESH_HexaBlocks::buildGroups(HEXA_NS::Document* doc)
1138 {
1139   if (MYDEBUG){
1140     MESSAGE("_addGroups() : : begin   <<<<<<");
1141     MESSAGE("_addGroups() : : nb. hexas= " << doc->countUsedHexa());
1142     MESSAGE("_addGroups() : : nb. quads= " << doc->countUsedQuad());
1143     MESSAGE("_addGroups() : : nb. edges= " << doc->countUsedEdge());
1144     MESSAGE("_addGroups() : : nb. nodes= " << doc->countUsedVertex());
1145   }
1146   // Looping on each groups of the document
1147   for ( int i=0; i < doc->countGroup(); i++ ){
1148       _fillGroup( doc->getGroup(i) );
1149   };
1150
1151   if(MYDEBUG) MESSAGE("_addGroups() : end  >>>>>>>>");
1152 }
1153
1154 // --------------------------------------------------------------
1155 //                      PRIVATE METHODS
1156 // --------------------------------------------------------------
1157 double SMESH_HexaBlocks::_Xx( double i, HEXA_NS::Law law, double nbNodes) //, double pos0 )
1158 {
1159   double result;
1160   double u0;
1161
1162   HEXA_NS::KindLaw k = law.getKind();
1163   double coeff       = law.getCoefficient();
1164   switch (k){
1165     case HEXA_NS::Uniform:
1166         result = (i+1)/(nbNodes+1);
1167         if(MYDEBUG) MESSAGE( "_Xx():" << " Uniform u("<<i<< ")"<< " = " << result);
1168         break;
1169     case HEXA_NS::Arithmetic:
1170         u0 = 1./(nbNodes + 1.) - (coeff*nbNodes)/2.;
1171 //         ASSERT(u0>0);
1172         if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Arithmetic discretization : check coefficient")));
1173         if (i==0){
1174           result = u0;
1175         } else {
1176           result = (i + 1.)*u0 + coeff*i*(i+1.)/2.;
1177         };
1178         if(MYDEBUG) MESSAGE( "_Xx():" << " Arithmetic u("<<i<< ")"<< " = " << result);
1179         break;
1180     case HEXA_NS::Geometric:
1181         u0 = (1.-coeff)/(1.-pow(coeff, nbNodes + 1) )  ;
1182 //         ASSERT(u0>0);
1183         if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Geometric discretization : check coefficient")));
1184         if (i==0){
1185           result = u0;
1186         } else {
1187           result = u0 * (1.- pow(coeff, i + 1) )/(1.-coeff) ;
1188         };
1189         if(MYDEBUG) MESSAGE( "_Xx():" << " Geometric u("<<i<< ")"<< " = " << result);
1190         break;
1191   }
1192   return result;
1193 }
1194
1195
1196 // ============================================================= _edgeLength
1197 double SMESH_HexaBlocks::_edgeLength(const TopoDS_Edge & E)
1198 {
1199   if(MYDEBUG) MESSAGE("_edgeLength() : : begin   <<<<<<");
1200   double UMin = 0, UMax = 0;
1201   if (BRep_Tool::Degenerated(E))
1202     return 0;
1203   TopLoc_Location L;
1204   Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, UMin, UMax);
1205   GeomAdaptor_Curve AdaptCurve(C);
1206   double length = GCPnts_AbscissaPoint::Length(AdaptCurve, UMin, UMax);
1207   if(MYDEBUG) MESSAGE("_edgeLength() : end  >>>>>>>>");
1208   return length;
1209 }
1210 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1211 // ============================================================== _buildMyCurve
1212 // === construire ma courbe a moi
1213 void SMESH_HexaBlocks::_buildMyCurve(
1214     const gp_Pnt&                               myCurve_pt_start,  //IN
1215     const gp_Pnt&                               myCurve_pt_end,    //IN
1216     std::list< BRepAdaptor_Curve* >&            myCurve_list,        //INOUT
1217     double&                                     myCurve_tot_len, //INOUT
1218     std::map< BRepAdaptor_Curve*, double>&      myCurve_lengths,//INOUT
1219     std::map< BRepAdaptor_Curve*, bool>&        myCurve_ways,   //INOUT
1220     std::map< BRepAdaptor_Curve*, double>&      myCurve_starts,   //INOUT
1221     HEXA_NS::Edge&                              edge) // For error diagnostic
1222 {
1223     if(MYDEBUG) MESSAGE("_buildMyCurve() : : begin   <<<<<<");
1224     bool current_way  = true;
1225     myCurve_tot_len    = 0.;
1226     BRepAdaptor_Curve* thePreviousCurve = NULL;
1227     BRepAdaptor_Curve* theCurve         = NULL;
1228
1229     gp_Pnt  curv_start, curv_end;
1230     gp_Pnt  p_curv_start , p_curv_end;
1231     double u_start     = 0;
1232     double curv_length = 0;
1233
1234     int nbass = edge.countAssociation();
1235     for (int nro = 0 ; nro < nbass ; ++nro)
1236         {
1237         theCurve = make_curve (curv_start, curv_end, curv_length, u_start,
1238                                edge, nro);
1239         if (theCurve != NULL)
1240            {
1241             bool    sens   = true;
1242             if ( thePreviousCurve == NULL )
1243                {
1244                // setting current_way and first curve way
1245                if ( myCurve_pt_start.IsEqual(curv_start, HEXA_EPS) )
1246                   {
1247                   current_way = true;
1248                   sens        = true;
1249                   }
1250                else if ( myCurve_pt_start.IsEqual(curv_end, HEXA_EPS) )
1251                   {
1252                   current_way = true;
1253                   sens        = false;
1254                   }
1255                else if ( myCurve_pt_end.IsEqual(curv_end, HEXA_EPS) )
1256                   {
1257                   current_way = false;
1258                   sens        = true;
1259                   }
1260                else if ( myCurve_pt_end.IsEqual(curv_start, HEXA_EPS) )
1261                   {
1262                   current_way = false;
1263                   sens        = false;
1264                   }
1265                else
1266                   {
1267                  if (MYDEBUG)
1268                      MESSAGE("SOMETHING WRONG on edge association... Bad script?");
1269                  edge.dumpAsso();
1270                  throw (SALOME_Exception (LOCALIZED("Edge association : check association parameters ( first, last ) between HEXA model and CAO")));
1271                    }
1272
1273                }
1274                 // it is not the first or last curve.
1275                 // ways are calculated between previous and new one.
1276             else
1277                {
1278                if (    p_curv_end.IsEqual( curv_end, HEXA_EPS  )
1279                     || p_curv_start.IsEqual( curv_start, HEXA_EPS ) )
1280                   {
1281                   sens  = NOT myCurve_ways[thePreviousCurve];// opposite WAY
1282                   }
1283                else if (   p_curv_end.IsEqual( curv_start, HEXA_EPS )
1284                         || p_curv_start.IsEqual( curv_end, HEXA_EPS ) )
1285                   {
1286                   sens  = myCurve_ways[thePreviousCurve];// same WAY
1287                   }
1288                else
1289                   {
1290                   if (MYDEBUG)
1291                      MESSAGE("SOMETHING WRONG on edge association... bad script?");
1292 //                     ASSERT(false);
1293                   edge.dumpAsso();
1294                   throw (SALOME_Exception(LOCALIZED("Edge association : Check association parameters ( first, last ) between HEXA model and CAO")));
1295                 }
1296             }
1297
1298             myCurve_list.push_back( theCurve );
1299             myCurve_ways   [theCurve] = sens;
1300             myCurve_starts [theCurve] = u_start;
1301             myCurve_lengths[theCurve] = curv_length;
1302             myCurve_tot_len          += curv_length;
1303
1304             thePreviousCurve  = theCurve;
1305             p_curv_start = curv_start;
1306             p_curv_end   = curv_end;
1307            }//if ( theCurveLength > 0 )
1308     }// for
1309
1310
1311     if ( NOT current_way ){
1312         std::list< BRepAdaptor_Curve* > tmp( myCurve_list.size() );
1313         std::copy( myCurve_list.rbegin(), myCurve_list.rend(), tmp.begin() );
1314         myCurve_list = tmp;
1315     }
1316
1317     if (MYDEBUG) {
1318       MESSAGE("current_way  was :" << current_way);
1319       MESSAGE("_buildMyCurve() : end  >>>>>>>>");
1320     }
1321 }
1322
1323
1324
1325
1326 gp_Pnt SMESH_HexaBlocks::_getPtOnMyCurve(
1327     const double&                             myCurve_u,      //IN
1328     std::map< BRepAdaptor_Curve*, bool>&      myCurve_ways,   //IN
1329     std::map< BRepAdaptor_Curve*, double>&    myCurve_lengths,//IN
1330     std::map< BRepAdaptor_Curve*, double>&    myCurve_starts, //IN
1331     std::list< BRepAdaptor_Curve* >&          myCurve_list,        //INOUT
1332     double&                                   myCurve_start ) //INOUT
1333 //     std::map< BRepAdaptor_Curve*, double>&  myCurve_firsts,
1334 //     std::map< BRepAdaptor_Curve*, double>&  myCurve_lasts,
1335 {
1336   if(MYDEBUG) MESSAGE("_getPtOnMyCurve() : : begin   <<<<<<");
1337   gp_Pnt ptOnMyCurve;
1338
1339   // looking for curve which contains parameter myCurve_u
1340   BRepAdaptor_Curve* curve      = myCurve_list.front();
1341   double            curve_start = myCurve_start;
1342   double            curve_end   = curve_start + myCurve_lengths[curve];
1343   double            curve_u;
1344   GCPnts_AbscissaPoint discret;
1345
1346   if (MYDEBUG){
1347     MESSAGE("looking for curve               = "<<(long) curve);
1348     MESSAGE("looking for curve: curve_u      = "<<myCurve_u);
1349     MESSAGE("looking for curve: curve_start  = "<<curve_start);
1350     MESSAGE("looking for curve: curve_end    = "<<curve_end);
1351     MESSAGE("looking for curve: curve_lenght = "<<myCurve_lengths[curve]);
1352     MESSAGE("looking for curve: curve.size _lenght= "<<myCurve_list.size());
1353   }
1354   while ( !( (myCurve_u >= curve_start) &&  (myCurve_u <= curve_end) ) ) {
1355
1356     ASSERT( myCurve_list.size() != 0 );
1357     myCurve_list.pop_front();
1358     curve       = myCurve_list.front();
1359     curve_start = curve_end;
1360     curve_end   = curve_start + myCurve_lengths[curve];
1361     if (MYDEBUG){
1362       MESSAGE("go next curve: curve_lenght = "<<myCurve_lengths[curve]);
1363       MESSAGE("go next curve: curve_start  = "<<curve_start);
1364       MESSAGE("go next curve: curve_end    = "<<curve_end);
1365       MESSAGE("go next curve: myCurve_u    = "<<myCurve_u);
1366     }
1367   }
1368   myCurve_start = curve_start;
1369
1370   // compute point
1371   if ( myCurve_ways[curve] ){
1372     discret = GCPnts_AbscissaPoint( *curve, (myCurve_u - curve_start), myCurve_starts[curve] );
1373   } else {
1374     discret = GCPnts_AbscissaPoint( *curve, myCurve_lengths[curve]- (myCurve_u - curve_start), myCurve_starts[curve] );
1375   }
1376   // PutData (discret);
1377   ASSERT(discret.IsDone());
1378   curve_u = discret.Parameter();
1379   ptOnMyCurve = curve->Value( curve_u );
1380
1381   if (MYDEBUG){
1382     MESSAGE("curve found!");
1383     MESSAGE("curve_u = "<< curve_u);
1384     MESSAGE("curve way = "<< myCurve_ways[curve]);
1385     MESSAGE("_getPtOnMyCurve() : end  >>>>>>>>");
1386   }
1387   return ptOnMyCurve;
1388 }
1389
1390
1391
1392
1393 void SMESH_HexaBlocks::_nodeInterpolationUV(double u, double v,
1394     SMDS_MeshNode* Pg, SMDS_MeshNode* Pd, SMDS_MeshNode* Ph, SMDS_MeshNode* Pb,
1395     SMDS_MeshNode* S1, SMDS_MeshNode* S2, SMDS_MeshNode* S3, SMDS_MeshNode* S4,
1396     double& xOut, double& yOut, double& zOut )
1397 {
1398   if (MYDEBUG){
1399     MESSAGE("_nodeInterpolationUV() IN:");
1400     MESSAGE("u ( "<< u <<" )");
1401     MESSAGE("v ( "<< v <<" )");
1402
1403     MESSAGE("S1 (" << S1->X() << "," << S1->Y() << "," << S1->Z() << ")");
1404     MESSAGE("S2 (" << S2->X() << "," << S2->Y() << "," << S2->Z() << ")");
1405     MESSAGE("S4 (" << S4->X() << "," << S4->Y() << "," << S4->Z() << ")");
1406     MESSAGE("S3 (" << S3->X() << "," << S3->Y() << "," << S3->Z() << ")");
1407
1408     MESSAGE("Pg (" << Pg->X() << "," << Pg->Y() << "," << Pg->Z() << ")");
1409     MESSAGE("Pd (" << Pd->X() << "," << Pd->Y() << "," << Pd->Z() << ")");
1410     MESSAGE("Ph (" << Ph->X() << "," << Ph->Y() << "," << Ph->Z() << ")");
1411     MESSAGE("Pb (" << Pb->X() << "," << Pb->Y() << "," << Pb->Z() << ")");
1412   }
1413
1414   xOut = ((1.-u)*Pg->X() + v*Ph->X() + u*Pd->X() + (1.-v)*Pb->X()) - (1.-u)*(1.-v)*S1->X() - u*(1.-v)*S2->X() - u*v*S3->X() - (1.-u)*v*S4->X();
1415   yOut = ((1.-u)*Pg->Y() + v*Ph->Y() + u*Pd->Y() + (1.-v)*Pb->Y()) - (1.-u)*(1.-v)*S1->Y() - u*(1.-v)*S2->Y() - u*v*S3->Y() - (1.-u)*v*S4->Y();
1416   zOut = ((1.-u)*Pg->Z() + v*Ph->Z() + u*Pd->Z() + (1.-v)*Pb->Z()) - (1.-u)*(1.-v)*S1->Z() - u*(1.-v)*S2->Z() - u*v*S3->Z() - (1.-u)*v*S4->Z();
1417
1418   if (MYDEBUG){
1419     MESSAGE("_nodeInterpolationUV() OUT("<<xOut<<","<<yOut<<","<<zOut<<" )");
1420   }
1421 }
1422
1423 // =========================================================== getFaceShapes
1424 TopoDS_Shape SMESH_HexaBlocks::getFaceShapes (Hex::Quad& quad)
1425 {
1426    int             nbass = quad.countAssociation ();
1427    Hex::FaceShape* face  = quad.getAssociation (0);
1428    if (nbass==1)
1429        return face->getShape ();
1430
1431    TopoDS_Compound compound;
1432    BRep_Builder    builder;
1433    builder.MakeCompound (compound);
1434
1435    for (int nro=0 ; nro < nbass ; nro++)
1436        {
1437        face = quad.getAssociation (nro);
1438        TopoDS_Shape shape = face->getShape ();
1439        builder.Add (compound, shape);
1440        }
1441
1442    return compound;
1443 }
1444
1445
1446 // ================================================== carre
1447 inline double carre (double val)
1448 {
1449   return val*val;
1450 }
1451 // ================================================== dist2
1452 inline double dist2 (const gp_Pnt& pt1, const gp_Pnt& pt2)
1453 {
1454    double dist = carre (pt2.X()-pt1.X()) + carre (pt2.Y()-pt1.Y())
1455                                          + carre (pt2.Z()-pt1.Z());
1456    return dist;
1457 }
1458 // ================================================== _intersect
1459 gp_Pnt SMESH_HexaBlocks::_intersect( const gp_Pnt& Pt,
1460                                      const gp_Vec& u, const gp_Vec& v,
1461                                      const TopoDS_Shape& shape,
1462                                      Standard_Real tol )
1463 {
1464   gp_Pnt result;
1465
1466   gp_Vec normale = u^v;
1467   gp_Dir dir(normale);
1468   gp_Lin li( Pt, dir );
1469
1470   Standard_Real s = -Precision::Infinite();
1471   Standard_Real e = +Precision::Infinite();
1472
1473   IntCurvesFace_ShapeIntersector inter;
1474   inter.Load(shape, tol);
1475 //   inter.Load(S, tol);
1476   inter.Perform(li, s, e);//inter.PerformNearest(li, s, e);
1477
1478 /***********************************************  Abu 2011-11-04 */
1479   if ( inter.IsDone() )
1480      {
1481      result = inter.Pnt(1);//first
1482      int nbrpts = inter.NbPnt();
1483      if (nbrpts>1)
1484         {
1485         double d0 = dist2 (result, Pt);
1486         for (int i=2; i <= inter.NbPnt(); ++i )
1487             {
1488             double d1 = dist2 (Pt, inter.Pnt(i));
1489             if (d1<d0)
1490                {
1491                d0 = d1;
1492                result = inter.Pnt (i);
1493                }
1494             }
1495         }
1496 /**********************************************  Abu 2011-11-04 (getEnd()) */
1497     if (MYDEBUG){
1498       MESSAGE("_intersect() : OK");
1499       for ( int i=1; i <= inter.NbPnt(); ++i ){
1500         gp_Pnt tmp = inter.Pnt(i);
1501         MESSAGE("_intersect() : pnt("<<i<<") = ("<<tmp.X()<<","<<tmp.Y()<<","<<tmp.Z()<<" )");
1502       }
1503     }
1504     _found +=1;
1505   } else {
1506     if(MYDEBUG) MESSAGE("_intersect() : KO");
1507     result = Pt;
1508     _notFound +=1;
1509   }
1510   _total+=1;
1511
1512   return result;
1513 }
1514
1515 // parameters q : IN,  v0: INOUT, v1: INOUT
1516 void SMESH_HexaBlocks::_searchInitialQuadWay( HEXA_NS::Quad* q, HEXA_NS::Vertex*& v0, HEXA_NS::Vertex*& v1 )
1517 {
1518   if(MYDEBUG) MESSAGE("_searchInitialQuadWay() : begin");
1519   v0 = NULL; v1 = NULL;
1520   if ( q->getNbrParents() != 1 ) return; // q must be a skin quad
1521
1522   HEXA_NS::Vertex* qA = q->getVertex(0);
1523   HEXA_NS::Vertex* qB = q->getVertex(1);
1524   HEXA_NS::Vertex* qC = q->getVertex(2);
1525   HEXA_NS::Vertex* qD = q->getVertex(3);
1526
1527   // searching for vertex on opposed quad
1528   HEXA_NS::Vertex *qAA = NULL, *qBB = NULL, *qCC = NULL, *qDD = NULL;
1529   HEXA_NS::Hexa* h = q->getParent(0);
1530   for( int i=0; i < h->countEdge(); ++i  ){
1531     HEXA_NS::Edge* e = h->getEdge(i);
1532     HEXA_NS::Vertex* e0 = e->getVertex(0);
1533     HEXA_NS::Vertex* e1 = e->getVertex(1);
1534
1535     if ( e0 == qA && e1 != qB && e1 != qC && e1 != qD ){
1536       qAA = e1;
1537     } else if ( e1 == qA && e0 != qB && e0 != qC && e0 != qD ){
1538       qAA = e0;
1539     } else if ( e0 == qB && e1 != qA && e1 != qC && e1 != qD ){
1540       qBB = e1;
1541     } else if ( e1 == qB && e0 != qA && e0 != qC && e0 != qD ){
1542       qBB = e0;
1543     } else if ( e0 == qC && e1 != qA && e1 != qB && e1 != qD ){
1544       qCC = e1;
1545     } else if ( e1 == qC && e0 != qA && e0 != qB && e0 != qD ){
1546       qCC = e0;
1547     } else if ( e0 == qD && e1 != qA && e1 != qB && e1 != qC ){
1548       qDD = e1;
1549     } else if ( e1 == qD && e0 != qA && e0 != qB && e0 != qC ){
1550       qDD = e0;
1551     }
1552   }
1553
1554   // working on final value ( point on CAO ), not on model
1555   SMDS_MeshNode *nA = _node[qA], *nAA = _node[qAA];
1556   SMDS_MeshNode *nB = _node[qB], *nBB = _node[qBB];
1557   SMDS_MeshNode *nC = _node[qC], *nCC = _node[qCC];
1558   SMDS_MeshNode *nD = _node[qD], *nDD = _node[qDD];
1559
1560   gp_Pnt pA( nA->X(), nA->Y(), nA->Z() );
1561   gp_Pnt pB( nB->X(), nB->Y(), nB->Z() );
1562   gp_Pnt pC( nC->X(), nC->Y(), nC->Z() );
1563   gp_Pnt pD( nD->X(), nD->Y(), nD->Z() );
1564
1565   gp_Pnt pAA( nAA->X(), nAA->Y(), nAA->Z() );
1566   gp_Pnt pBB( nBB->X(), nBB->Y(), nBB->Z() );
1567   gp_Pnt pCC( nCC->X(), nCC->Y(), nCC->Z() );
1568   gp_Pnt pDD( nDD->X(), nDD->Y(), nDD->Z() );
1569
1570   gp_Vec AB( pA, pB );
1571   gp_Vec AC( pA, pC );
1572   gp_Vec normP = AB^AC;
1573   gp_Dir dirP( normP );
1574
1575   // building plane for point projection
1576   gp_Pln plnP( gp_Pnt(nA->X(), nA->Y(), nA->Z()), dirP);
1577   TopoDS_Shape sPlnP = BRepBuilderAPI_MakeFace(plnP).Face();
1578
1579   // PAAA is the result of PAA projection
1580   gp_Pnt pAAA = _intersect( pAA, AB, AC, sPlnP );
1581   gp_Pnt pBBB = _intersect( pBB, AB, AC, sPlnP );
1582   gp_Pnt pCCC = _intersect( pCC, AB, AC, sPlnP );
1583   gp_Pnt pDDD = _intersect( pDD, AB, AC, sPlnP );
1584
1585   gp_Dir AA( gp_Vec(pAA, pAAA) );
1586   gp_Dir BB( gp_Vec(pBB, pBBB) );
1587   gp_Dir CC( gp_Vec(pCC, pCCC) );
1588   gp_Dir DD( gp_Vec(pDD, pDDD) );
1589
1590   // eventually, we are able to know if the input quad is a good client!
1591   // exit the fonction otherwise
1592   if ( AA.IsOpposite(BB, HEXA_QUAD_WAY) ) return;
1593   if ( BB.IsOpposite(CC, HEXA_QUAD_WAY) ) return;
1594   if ( CC.IsOpposite(DD, HEXA_QUAD_WAY) ) return;
1595
1596   // ok, give the input quad the good orientation by
1597   // setting 2 vertex
1598   if ( !dirP.IsOpposite(AA, HEXA_QUAD_WAY) ) { //OK
1599       v0 = qA; v1 = qB;
1600   } else {
1601       v0 = qB; v1 = qA;
1602   }
1603
1604   if(MYDEBUG) MESSAGE("_searchInitialQuadWay() : end");
1605 }
1606
1607 SMESH_Group* SMESH_HexaBlocks::_createGroup(HEXA_NS::Group& grHex)
1608 {
1609   if(MYDEBUG) MESSAGE("_createGroup() : : begin   <<<<<<");
1610
1611   std::string aGrName           = grHex.getName();
1612   HEXA_NS::EnumGroup grHexKind  = grHex.getKind();
1613
1614   if(MYDEBUG) MESSAGE("aGrName"<<aGrName);
1615
1616   SMDSAbs_ElementType aGrType;
1617   switch ( grHexKind ){
1618     case HEXA_NS::HexaCell   : aGrType = SMDSAbs_Volume; break;
1619     case HEXA_NS::QuadCell   : aGrType = SMDSAbs_Face  ; break;
1620     case HEXA_NS::EdgeCell   : aGrType = SMDSAbs_Edge  ; break;
1621     case HEXA_NS::HexaNode   : aGrType = SMDSAbs_Node  ; break;
1622     case HEXA_NS::QuadNode   : aGrType = SMDSAbs_Node  ; break;
1623     case HEXA_NS::EdgeNode   : aGrType = SMDSAbs_Node  ; break;
1624     case HEXA_NS::VertexNode : aGrType = SMDSAbs_Node  ; break;
1625     default : ASSERT(false);
1626   }
1627
1628   int aId;
1629   SMESH_Group* aGr = _theMesh->AddGroup(aGrType, aGrName.c_str(), aId);
1630
1631   if(MYDEBUG) MESSAGE("_createGroup() : end  >>>>>>>>");
1632   return aGr;
1633 }
1634
1635 void SMESH_HexaBlocks::_fillGroup(HEXA_NS::Group* grHex )
1636 {
1637   if(MYDEBUG) MESSAGE("_fillGroup() : : begin   <<<<<<");
1638
1639   SMESH_Group* aGr = _createGroup( *grHex );
1640   HEXA_NS::EltBase*  grHexElt   = NULL;
1641   HEXA_NS::EnumGroup grHexKind  = grHex->getKind();
1642   int                grHexNbElt = grHex->countElement();
1643
1644   if(MYDEBUG) MESSAGE("_fillGroup() : kind = " << grHexKind);
1645   if(MYDEBUG) MESSAGE("_fillGroup() : count= " << grHexNbElt);
1646
1647   // A)Looking for elements ID
1648   std::vector<const SMDS_MeshElement*> aGrEltIDs;
1649
1650   for ( int n=0; n<grHexNbElt; ++n ){
1651       grHexElt = grHex->getElement(n);
1652
1653       switch ( grHexKind ){
1654         case HEXA_NS::HexaCell:
1655         {
1656             HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1657             if ( _volumesOnHexa.count(h)>0 ){
1658               SMESHVolumes volumes = _volumesOnHexa[h];
1659               for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1660                   aGrEltIDs.push_back(*aVolume);
1661               }
1662             } else {
1663               if(MYDEBUG) MESSAGE("GROUP OF VOLUME: volume for hexa (id = "<<h->getId()<<") not found");
1664             }
1665         }
1666         break;
1667         case HEXA_NS::QuadCell:
1668         {
1669             HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1670             if ( _facesOnQuad.count(q)>0 ){
1671               SMESHFaces faces = _facesOnQuad[q];
1672               for ( SMESHFaces::iterator aFace = faces.begin(); aFace != faces.end(); ++aFace ){
1673                   aGrEltIDs.push_back(*aFace);
1674               }
1675             } else {
1676               if(MYDEBUG) MESSAGE("GROUP OF FACE: face for quad (id = "<<q->getId()<<") not found");
1677             }
1678         }
1679         break;
1680         case HEXA_NS::EdgeCell:
1681         {
1682             HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1683             if ( _edgesOnEdge.count(e)>0 ){
1684               SMESHEdges edges = _edgesOnEdge[e];
1685               for ( SMESHEdges::iterator anEdge = edges.begin(); anEdge != edges.end(); ++anEdge ){
1686                   aGrEltIDs.push_back(*anEdge);
1687               }
1688             } else {
1689               if(MYDEBUG) MESSAGE("GROUP OF Edge: edge for edge (id = "<<e->getId()<<") not found");
1690             }
1691         }
1692         break;
1693         case HEXA_NS::HexaNode:
1694         {
1695             HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1696             if ( _volumesOnHexa.count(h)>0 ){
1697               SMESHVolumes volumes = _volumesOnHexa[h];
1698               for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1699                 SMDS_ElemIteratorPtr aNodeIter = (*aVolume)->nodesIterator();
1700                 while( aNodeIter->more() ){
1701                   const SMDS_MeshNode* aNode =
1702                     dynamic_cast<const SMDS_MeshNode*>( aNodeIter->next() );
1703                   if ( aNode ){
1704                       aGrEltIDs.push_back(aNode);
1705                   }
1706                 }
1707               }
1708             } else {
1709               if(MYDEBUG) MESSAGE("GROUP OF HEXA NODES: nodes on hexa  (id = "<<h->getId()<<") not found");
1710             }
1711         }
1712         break;
1713         case HEXA_NS::QuadNode:
1714         {
1715             HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1716             if ( _quadNodes.count(q)>0 ){
1717               ArrayOfSMESHNodes nodesOnQuad = _quadNodes[q];
1718               for ( ArrayOfSMESHNodes::iterator nodes = nodesOnQuad.begin(); nodes != nodesOnQuad.end(); ++nodes){
1719                 for ( SMESHNodes::iterator aNode = nodes->begin(); aNode != nodes->end(); ++aNode){
1720                   aGrEltIDs.push_back(*aNode);
1721                 }
1722               }
1723             } else {
1724               if(MYDEBUG) MESSAGE("GROUP OF QUAD NODES: nodes on quad (id = "<<q->getId()<<") not found");
1725             }
1726         }
1727         break;
1728         case HEXA_NS::EdgeNode:
1729         {
1730             HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1731             if ( _nodesOnEdge.count(e)>0 ){
1732               SMESHNodes nodes = _nodesOnEdge[e];
1733               for ( SMESHNodes::iterator aNode = nodes.begin(); aNode != nodes.end(); ++aNode){
1734                 aGrEltIDs.push_back(*aNode);
1735               }
1736             } else {
1737               if(MYDEBUG) MESSAGE("GROUP OF EDGE NODES: nodes on edge (id = "<<e->getId()<<") not found");
1738             }
1739         }
1740         break;
1741         case HEXA_NS::VertexNode:
1742         {
1743           HEXA_NS::Vertex* v = reinterpret_cast<HEXA_NS::Vertex*>(grHexElt);
1744             if ( _node.count(v)>0 ){
1745               aGrEltIDs.push_back(_node[v]);
1746             } else {
1747               if(MYDEBUG) MESSAGE("GROUP OF VERTEX NODES: nodes for vertex (id = "<<v->getId()<<") not found");
1748             }
1749         }
1750         break;
1751         default : ASSERT(false);
1752       }
1753   }
1754
1755   // B)Filling the group on SMESH
1756   SMESHDS_Group* aGroupDS = dynamic_cast<SMESHDS_Group*>( aGr->GetGroupDS() );
1757
1758   for ( int i=0; i < aGrEltIDs.size(); i++ ) {
1759     aGroupDS->SMDSGroup().Add( aGrEltIDs[i] );
1760   };
1761
1762   if(MYDEBUG) MESSAGE("_fillGroup() : end  >>>>>>>>");
1763 }