src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx

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