src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx

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