src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx

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