src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx

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