]> SALOME platform Git repositories - modules/smesh.git/commitdiff
Salome HOME
Avoid anonymous namespaces and 'using' instructions in headers. Simplify tests code...
authorjfa <jfa@opencascade.com>
Fri, 19 Jul 2024 11:53:33 +0000 (12:53 +0100)
committerjfa <jfa@opencascade.com>
Fri, 19 Jul 2024 11:53:33 +0000 (12:53 +0100)
src/StdMeshers.test/HexahedronCanonicalShapesTest.cxx
src/StdMeshers.test/HexahedronTest.cxx
src/StdMeshers/StdMeshers_Cartesian_3D.cxx
src/StdMeshers/StdMeshers_Cartesian_3D_Grid.cxx
src/StdMeshers/StdMeshers_Cartesian_3D_Grid.hxx
src/StdMeshers/StdMeshers_Cartesian_3D_Hexahedron.cxx
src/StdMeshers/StdMeshers_Cartesian_3D_Hexahedron.hxx
test/data/HexahedronTest/NRTMJ4.brep

index fe3ab4c95ae10882c07300ddb8f99ef5ec9c0aa0..3f4b16f5620911127d0a35897fc186e97c87cea9 100644 (file)
 #include <iostream>
 #include <memory>
 
+using namespace StdMeshers::Cartesian3D;
+
 // Helper functions!
 // Build Grid
 //      Require building mesh
-//      Require building shape. For test load shapes from memory in .brep files seems the simplest
-//      
+//      Require building shape.
 
 /*!
   * \brief Mock mesh
@@ -69,62 +70,27 @@ struct CartesianHypo: public StdMeshers_CartesianParameters3D
 /*!
   * \brief Initialize the grid and intesersectors of grid with the geometry
   */
-void GridInitAndIntersectWithShapeGrid& grid,
+void GridInitAndIntersectWithShape (Grid& grid,
                                     double gridSpacing,
                                     double theSizeThreshold,
                                     const TopoDS_Shape theShape, 
-                                    std::map< TGeomID, vector< TGeomID > >& edge2faceIDsMap,
-                                    const int /*numOfThreads*/ )
+                                    TEdge2faceIDsMap& edge2faceIDsMap,
+                                    const int theNumOfThreads)
 {
-  std::vector< TopoDS_Shape > faceVec;
-  TopTools_MapOfShape faceMap;
-  TopExp_Explorer fExp;
-  for ( fExp.Init( theShape, TopAbs_FACE ); fExp.More(); fExp.Next() )
-  {
-    bool isNewFace = faceMap.Add( fExp.Current() );
-    if ( !grid._toConsiderInternalFaces )
-      if ( !isNewFace || fExp.Current().Orientation() == TopAbs_INTERNAL )
-        // remove an internal face
-        faceMap.Remove( fExp.Current() );
-  }
-  faceVec.reserve( faceMap.Extent() );
-  faceVec.assign( faceMap.cbegin(), faceMap.cend() );
-  
-  vector<FaceGridIntersector> facesItersectors( faceVec.size() );
-
-  Bnd_Box shapeBox;
-  for ( size_t i = 0; i < faceVec.size(); ++i )
-  {
-    facesItersectors[i]._face   = TopoDS::Face( faceVec[i] );
-    facesItersectors[i]._faceID = grid.ShapeID( faceVec[i] );
-    facesItersectors[i]._grid   = &grid;
-    shapeBox.Add( facesItersectors[i].GetFaceBndBox() );
-  }
   // Canonical axes(i,j,k)
   double axisDirs[9] = {1.,0.,0.,0.,1.,0.,0.,0.,1.};
-
-  Tools::GetExactBndBox( faceVec, axisDirs, shapeBox );
-  vector<double> xCoords, yCoords, zCoords;  
-  std::unique_ptr<CartesianHypo> myHypo( new CartesianHypo() );  
   std::vector<std::string> grdSpace = { std::to_string(gridSpacing) };
   std::vector<double> intPnts;
-  myHypo->SetGridSpacing(grdSpace, intPnts, 0 ); // Spacing in dir 0
-  myHypo->SetGridSpacing(grdSpace, intPnts, 1 ); // Spacing in dir 1
-  myHypo->SetGridSpacing(grdSpace, intPnts, 2 ); // Spacing in dir 2
-  myHypo->SetSizeThreshold(theSizeThreshold);    // set threshold
-  myHypo->GetCoordinates(xCoords, yCoords, zCoords, shapeBox);
-  grid.SetCoordinates( xCoords, yCoords, zCoords, axisDirs, shapeBox );
-
-  for ( size_t i = 0; i < facesItersectors.size(); ++i )
-    facesItersectors[i].Intersect();
-  
-  for ( size_t i = 0; i < facesItersectors.size(); ++i )
-    facesItersectors[i].StoreIntersections();
-
-  grid.ComputeNodes( *grid._helper );
-  grid.GetEdgesToImplement( edge2faceIDsMap, theShape, faceVec );
-}
 
+  std::unique_ptr<CartesianHypo> aHypo ( new CartesianHypo() );  
+  aHypo->SetGridSpacing(grdSpace, intPnts, 0 ); // Spacing in dir 0
+  aHypo->SetGridSpacing(grdSpace, intPnts, 1 ); // Spacing in dir 1
+  aHypo->SetGridSpacing(grdSpace, intPnts, 2 ); // Spacing in dir 2
+  aHypo->SetSizeThreshold(theSizeThreshold);    // set threshold
+  aHypo->SetAxisDirs(axisDirs);
+
+  grid.GridInitAndInterserctWithShape(theShape, edge2faceIDsMap, aHypo.get(), theNumOfThreads, false);
+}
 
 /*!
   * \brief Test runner
@@ -139,6 +105,7 @@ bool testShape (const TopoDS_Shape theShape,
   std::unique_ptr<SMESH_Mesh> aMesh( new SMESH_Mesh_Test() );
   aMesh->ShapeToMesh( theShape );
   SMESH_MesherHelper helper( *aMesh );
+
   Grid grid;
   grid._helper = &helper;
   grid._toAddEdges = toAddEdges;
@@ -147,11 +114,11 @@ bool testShape (const TopoDS_Shape theShape,
   grid._toUseThresholdForInternalFaces = false;
   grid._toUseQuanta = false;
   grid._sizeThreshold = theSizeThreshold;
-  grid.InitGeometry( theShape );
 
-  std::map< TGeomID, vector< TGeomID > > edge2faceIDsMap;
+  TEdge2faceIDsMap edge2faceIDsMap;
   GridInitAndIntersectWithShape( grid, theGridSpacing, theSizeThreshold,
                                  theShape, edge2faceIDsMap, 1 );
+
   Hexahedron hex( &grid );
   int nbAdded = hex.MakeElements( helper, edge2faceIDsMap, 1 );
   if (nbAdded != theNbCreatedExpected) {
index e993a96f58d75a22fa091edaff41218292496149..dd78c796e4b6aeadd591a82ea3b65224ede0b023 100644 (file)
 #include <iostream>
 #include <memory>
 
+using namespace StdMeshers::Cartesian3D;
+
 // Helper functions!
 // Build Grid
 //      Require building mesh
-//      Require building shape. For test load shapes from memory in .brep files seems the simplest
-//      
+//      Require building shape.
 
 /*!
   * \brief Mock mesh
@@ -78,60 +79,26 @@ void loadBrepShape( std::string shapeName, TopoDS_Shape & shape )
 /*!
   * \brief Initialize the grid and intesersectors of grid with the geometry
   */
-void GridInitAndIntersectWithShapeGrid& grid,
+void GridInitAndIntersectWithShape (Grid& grid,
                                     double gridSpacing,
                                     double theSizeThreshold,
                                     const TopoDS_Shape theShape, 
-                                    std::map< TGeomID, vector< TGeomID > >& edge2faceIDsMap,
-                                    const int /*numOfThreads*/ )
+                                    TEdge2faceIDsMap& edge2faceIDsMap,
+                                    const int theNumOfThreads)
 {
-  std::vector< TopoDS_Shape > faceVec;
-  TopTools_MapOfShape faceMap;
-  TopExp_Explorer fExp;
-  for ( fExp.Init( theShape, TopAbs_FACE ); fExp.More(); fExp.Next() )
-  {
-    bool isNewFace = faceMap.Add( fExp.Current() );
-    if ( !grid._toConsiderInternalFaces )
-      if ( !isNewFace || fExp.Current().Orientation() == TopAbs_INTERNAL )
-        // remove an internal face
-        faceMap.Remove( fExp.Current() );
-  }
-  faceVec.reserve( faceMap.Extent() );
-  faceVec.assign( faceMap.cbegin(), faceMap.cend() );
-  
-  vector<FaceGridIntersector> facesItersectors( faceVec.size() );
-
-  Bnd_Box shapeBox;
-  for ( size_t i = 0; i < faceVec.size(); ++i )
-  {
-    facesItersectors[i]._face   = TopoDS::Face( faceVec[i] );
-    facesItersectors[i]._faceID = grid.ShapeID( faceVec[i] );
-    facesItersectors[i]._grid   = &grid;
-    shapeBox.Add( facesItersectors[i].GetFaceBndBox() );
-  }
   // Canonical axes(i,j,k)
-  double axisDirs[9] = {1.,0.,0.,0.,1.,0.,0.,0.,1.};
-
-  Tools::GetExactBndBox( faceVec, axisDirs, shapeBox );
-  vector<double> xCoords, yCoords, zCoords;  
-  std::unique_ptr<CartesianHypo> myHypo( new CartesianHypo() );  
+  double axisDirs[9] = {1.,0.,0., 0.,1.,0., 0.,0.,1.};
   std::vector<std::string> grdSpace = { std::to_string(gridSpacing) };
   std::vector<double> intPnts;
-  myHypo->SetGridSpacing(grdSpace, intPnts, 0 ); // Spacing in dir 0
-  myHypo->SetGridSpacing(grdSpace, intPnts, 1 ); // Spacing in dir 1
-  myHypo->SetGridSpacing(grdSpace, intPnts, 2 ); // Spacing in dir 2
-  myHypo->SetSizeThreshold(theSizeThreshold);    // set threshold
-  myHypo->GetCoordinates(xCoords, yCoords, zCoords, shapeBox);
-  grid.SetCoordinates( xCoords, yCoords, zCoords, axisDirs, shapeBox );
-
-  for ( size_t i = 0; i < facesItersectors.size(); ++i )
-    facesItersectors[i].Intersect();
-  
-  for ( size_t i = 0; i < facesItersectors.size(); ++i )
-    facesItersectors[i].StoreIntersections();
-
-  grid.ComputeNodes( *grid._helper );
-  grid.GetEdgesToImplement( edge2faceIDsMap, theShape, faceVec );
+
+  std::unique_ptr<CartesianHypo> aHypo ( new CartesianHypo() );  
+  aHypo->SetAxisDirs(axisDirs);
+  aHypo->SetGridSpacing(grdSpace, intPnts, 0 ); // Spacing in dir 0
+  aHypo->SetGridSpacing(grdSpace, intPnts, 1 ); // Spacing in dir 1
+  aHypo->SetGridSpacing(grdSpace, intPnts, 2 ); // Spacing in dir 2
+  aHypo->SetSizeThreshold(theSizeThreshold);    // set threshold
+
+  grid.GridInitAndInterserctWithShape(theShape, edge2faceIDsMap, aHypo.get(), theNumOfThreads, false);
 }
 
 /*!
@@ -139,6 +106,10 @@ void GridInitAndIntersectWithShape( Grid& grid,
   */
 bool testNRTM1()
 {
+  TopoDS_Shape aShape;
+  loadBrepShape( "data/HexahedronTest/NRTM1.brep", aShape );
+  CPPUNIT_ASSERT_MESSAGE( "Could not load the brep shape!", !aShape.IsNull() );      
+
   const auto numOfCores = std::thread::hardware_concurrency() == 0 ? 1 : std::thread::hardware_concurrency();
   std::vector<int> numberOfThreads(numOfCores);
   std::iota (std::begin(numberOfThreads), std::end(numberOfThreads), 1);
@@ -147,20 +118,22 @@ bool testNRTM1()
   {
     for (size_t i = 0; i < 10 /*trials*/; i++)
     {
-      TopoDS_Shape myShape;
-      loadBrepShape( "data/HexahedronTest/NRTM1.brep", myShape );
-      CPPUNIT_ASSERT_MESSAGE( "Could not load the brep shape!", !myShape.IsNull() );      
-      std::unique_ptr<SMESH_Mesh> myMesh( new SMESH_Mesh_Test() );
-      myMesh->ShapeToMesh( myShape );
-      SMESH_MesherHelper helper( *myMesh );
+      std::unique_ptr<SMESH_Mesh> aMesh( new SMESH_Mesh_Test() );
+      aMesh->ShapeToMesh( aShape );
+      SMESH_MesherHelper helper( *aMesh );
+
       Grid grid;
       grid._helper = &helper;
-      grid._toAddEdges = false; grid._toCreateFaces = false; grid._toConsiderInternalFaces = false; grid._toUseThresholdForInternalFaces = false; grid._toUseQuanta = false;
+      grid._toAddEdges = false;
+      grid._toCreateFaces = false;
+      grid._toConsiderInternalFaces = false;
+      grid._toUseThresholdForInternalFaces = false;
+      grid._toUseQuanta = false;
       grid._sizeThreshold = 4.0;
-      grid.InitGeometry( myShape );
       
-      std::map< TGeomID, vector< TGeomID > > edge2faceIDsMap;
-      GridInitAndIntersectWithShape( grid, 1.0, 4.0, myShape, edge2faceIDsMap, nThreads );
+      TEdge2faceIDsMap edge2faceIDsMap;
+      GridInitAndIntersectWithShape( grid, 1.0, 4.0, aShape, edge2faceIDsMap, nThreads );
+
       Hexahedron hex( &grid );
       int nbAdded = hex.MakeElements( helper, edge2faceIDsMap, nThreads );
       CPPUNIT_ASSERT_MESSAGE( "Number of computed elements does not match", nbAdded == 1024 );
@@ -174,6 +147,10 @@ bool testNRTM1()
   */
 bool testNRTJ4()
 {
+  TopoDS_Shape aShape;
+  loadBrepShape( "data/HexahedronTest/NRTMJ4.brep", aShape );
+  CPPUNIT_ASSERT_MESSAGE( "Could not load the brep shape!", !aShape.IsNull() );      
+
   const auto numOfCores = std::thread::hardware_concurrency() == 0 ? 1 : std::thread::hardware_concurrency()/2;
   std::vector<int> numberOfThreads(numOfCores);
   std::iota (std::begin(numberOfThreads), std::end(numberOfThreads), 1);
@@ -183,22 +160,22 @@ bool testNRTJ4()
   {
     for (size_t i = 0; i < 10 /*trials*/; i++)
     {
-      TopoDS_Shape myShape;
-      loadBrepShape( "data/HexahedronTest/NRTMJ4.brep", myShape );
-      CPPUNIT_ASSERT_MESSAGE( "Could not load the brep shape!", !myShape.IsNull() );      
-      std::unique_ptr<SMESH_Mesh> myMesh( new SMESH_Mesh_Test() );
-      myMesh->ShapeToMesh( myShape );
-      SMESH_MesherHelper helper( *myMesh );
+      std::unique_ptr<SMESH_Mesh> aMesh( new SMESH_Mesh_Test() );
+      aMesh->ShapeToMesh( aShape );
+      SMESH_MesherHelper helper( *aMesh );
+
       Grid grid;
       grid._helper = &helper;
-      grid._toAddEdges = false; grid._toConsiderInternalFaces = false; grid._toUseThresholdForInternalFaces = false; grid._toUseQuanta = false;
+      grid._toAddEdges = false;
+      grid._toConsiderInternalFaces = false;
+      grid._toUseThresholdForInternalFaces = false;
+      grid._toUseQuanta = false;
       double testThreshold = 1.000001;
       grid._toCreateFaces = true;
       grid._sizeThreshold = testThreshold;
-      grid.InitGeometry( myShape );
       
-      std::map< TGeomID, vector< TGeomID > > edge2faceIDsMap;
-      GridInitAndIntersectWithShape( grid, 2.0, testThreshold, myShape, edge2faceIDsMap, nThreads );
+      TEdge2faceIDsMap edge2faceIDsMap;
+      GridInitAndIntersectWithShape( grid, 2.0, testThreshold, aShape, edge2faceIDsMap, nThreads );
       Hexahedron hex( &grid );
       int nbAdded = hex.MakeElements( helper, edge2faceIDsMap, nThreads );
       CPPUNIT_ASSERT_MESSAGE( "Number of computed elements does not match", nbAdded == 35150 );
@@ -210,6 +187,9 @@ bool testNRTJ4()
 // Entry point for test
 int main()
 {
-  bool isOK = testNRTM1() && testNRTJ4();
+  bool isOK = testNRTM1();
+  if (!testNRTJ4())
+    isOK = false;
+
   return isOK ? 0 : 1;
 }
index 98733b703e7b288d518260c4fdc0d406dd66c7ea..c20c5a2a079499402f9e1f9c0be2bce8d48a457c 100644 (file)
@@ -52,8 +52,9 @@
 
 using namespace std;
 using namespace SMESH;
-using namespace gridtools;
+using namespace StdMeshers::Cartesian3D;
 // using namespace facegridintersector;
+
 namespace 
 {
   /*!
@@ -229,7 +230,7 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh &         theMesh,
     }
 
     // remove free nodes
-    if ( SMESHDS_SubMesh * smDS = meshDS->MeshElements( helper.GetSubShapeID() ))
+    if ( /*SMESHDS_SubMesh * smDS = */meshDS->MeshElements( helper.GetSubShapeID() ))
     {
       std::vector< const SMDS_MeshNode* > nodesToRemove;
       // get intersection nodes
@@ -391,4 +392,4 @@ void StdMeshers_Cartesian_3D::setSubmeshesComputed(SMESH_Mesh&         theMesh,
 {
   for ( TopExp_Explorer soExp( theShape, TopAbs_SOLID ); soExp.More(); soExp.Next() )
     _EventListener::setAlwaysComputed( true, theMesh.GetSubMesh( soExp.Current() ));
-}
\ No newline at end of file
+}
index 0d02a9f730e231e23b370a494c55f493f1b0a624..273b51c195e722c5eccc14fdd7bce7eaf84c2673 100644 (file)
 #include <tbb/parallel_for.h>
 #endif
 
+using namespace std;
 using namespace SMESH;
-using namespace gridtools;
+using namespace StdMeshers::Cartesian3D;
+
 std::mutex _bMutex;
+
 //=============================================================================
 /*
   * Remove coincident intersection points
@@ -929,8 +932,10 @@ void Grid::ComputeNodes(SMESH_MesherHelper& helper)
   return;
 }
 
-bool Grid::GridInitAndInterserctWithShape( const TopoDS_Shape& theShape, std::map< TGeomID, vector< TGeomID > >& edge2faceIDsMap, 
-                                            const StdMeshers_CartesianParameters3D* hyp, const int numOfThreads, bool computeCanceled )
+bool Grid::GridInitAndInterserctWithShape( const TopoDS_Shape& theShape,
+                                           TEdge2faceIDsMap& edge2faceIDsMap, 
+                                           const StdMeshers_CartesianParameters3D* hyp,
+                                           const int /*numOfThreads*/, bool computeCanceled )
 {
   InitGeometry( theShape );
 
@@ -1088,3 +1093,351 @@ void Tools::GetExactBndBox( const vector< TopoDS_Shape >& faceVec, const double*
   return;
 }
 
+//=============================================================================
+/*
+ * Intersects TopoDS_Face with all GridLine's
+ */
+void FaceGridIntersector::Intersect()
+{
+  FaceLineIntersector intersector;
+  intersector._surfaceInt = GetCurveFaceIntersector();
+  intersector._tol        = _grid->_tol;
+  intersector._transOut   = _face.Orientation() == TopAbs_REVERSED ? Trans_IN : Trans_OUT;
+  intersector._transIn    = _face.Orientation() == TopAbs_REVERSED ? Trans_OUT : Trans_IN;
+
+  typedef void (FaceLineIntersector::* PIntFun )(const GridLine& gridLine);
+  PIntFun interFunction;
+
+  bool isDirect = true;
+  BRepAdaptor_Surface surf( _face );
+  switch ( surf.GetType() ) {
+  case GeomAbs_Plane:
+    intersector._plane = surf.Plane();
+    interFunction = &FaceLineIntersector::IntersectWithPlane;
+    isDirect = intersector._plane.Direct();
+    break;
+  case GeomAbs_Cylinder:
+    intersector._cylinder = surf.Cylinder();
+    interFunction = &FaceLineIntersector::IntersectWithCylinder;
+    isDirect = intersector._cylinder.Direct();
+    break;
+  case GeomAbs_Cone:
+    intersector._cone = surf.Cone();
+    interFunction = &FaceLineIntersector::IntersectWithCone;
+    //isDirect = intersector._cone.Direct();
+    break;
+  case GeomAbs_Sphere:
+    intersector._sphere = surf.Sphere();
+    interFunction = &FaceLineIntersector::IntersectWithSphere;
+    isDirect = intersector._sphere.Direct();
+    break;
+  case GeomAbs_Torus:
+    intersector._torus = surf.Torus();
+    interFunction = &FaceLineIntersector::IntersectWithTorus;
+    //isDirect = intersector._torus.Direct();
+    break;
+  default:
+    interFunction = &FaceLineIntersector::IntersectWithSurface;
+  }
+  if ( !isDirect )
+    std::swap( intersector._transOut, intersector._transIn );
+
+  _intersections.clear();
+  for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
+  {
+    if ( surf.GetType() == GeomAbs_Plane )
+    {
+      // check if all lines in this direction are parallel to a plane
+      if ( intersector._plane.Axis().IsNormal( _grid->_lines[iDir][0]._line.Position(),
+                                               Precision::Angular()))
+        continue;
+      // find out a transition, that is the same for all lines of a direction
+      gp_Dir plnNorm = intersector._plane.Axis().Direction();
+      gp_Dir lineDir = _grid->_lines[iDir][0]._line.Direction();
+      intersector._transition =
+        ( plnNorm * lineDir < 0 ) ? intersector._transIn : intersector._transOut;
+    }
+    if ( surf.GetType() == GeomAbs_Cylinder )
+    {
+      // check if all lines in this direction are parallel to a cylinder
+      if ( intersector._cylinder.Axis().IsParallel( _grid->_lines[iDir][0]._line.Position(),
+                                                    Precision::Angular()))
+        continue;
+    }
+
+    // intersect the grid lines with the face
+    for ( size_t iL = 0; iL < _grid->_lines[iDir].size(); ++iL )
+    {
+      GridLine& gridLine = _grid->_lines[iDir][iL];
+      if ( _bndBox.IsOut( gridLine._line )) continue;
+
+      intersector._intPoints.clear();
+      (intersector.*interFunction)( gridLine ); // <- intersection with gridLine
+      for ( size_t i = 0; i < intersector._intPoints.size(); ++i )
+        _intersections.push_back( std::make_pair( &gridLine, intersector._intPoints[i] ));
+    }
+  }
+
+  if ( _face.Orientation() == TopAbs_INTERNAL )
+  {
+    for ( size_t i = 0; i < _intersections.size(); ++i )
+      if ( _intersections[i].second._transition == Trans_IN ||
+           _intersections[i].second._transition == Trans_OUT )
+      {
+        _intersections[i].second._transition = Trans_INTERNAL;
+      }
+  }
+  return;
+}
+
+#ifdef WITH_TBB
+//================================================================================
+/*
+ * check if its face can be safely intersected in a thread
+ */
+bool FaceGridIntersector::IsThreadSafe(std::set< const Standard_Transient* >& noSafeTShapes) const
+{
+  bool isSafe = true;
+
+  // check surface
+  TopLoc_Location loc;
+  Handle(Geom_Surface) surf = BRep_Tool::Surface( _face, loc );
+  Handle(Geom_RectangularTrimmedSurface) ts =
+    Handle(Geom_RectangularTrimmedSurface)::DownCast( surf );
+  while( !ts.IsNull() ) {
+    surf = ts->BasisSurface();
+    ts = Handle(Geom_RectangularTrimmedSurface)::DownCast(surf);
+  }
+  if ( surf->IsKind( STANDARD_TYPE(Geom_BSplineSurface )) ||
+       surf->IsKind( STANDARD_TYPE(Geom_BezierSurface )))
+    if ( !noSafeTShapes.insert( _face.TShape().get() ).second )
+      isSafe = false;
+
+  double f, l;
+  TopExp_Explorer exp( _face, TopAbs_EDGE );
+  for ( ; exp.More(); exp.Next() )
+  {
+    bool edgeIsSafe = true;
+    const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
+    // check 3d curve
+    {
+      Handle(Geom_Curve) c = BRep_Tool::Curve( e, loc, f, l);
+      if ( !c.IsNull() )
+      {
+        Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(c);
+        while( !tc.IsNull() ) {
+          c = tc->BasisCurve();
+          tc = Handle(Geom_TrimmedCurve)::DownCast(c);
+        }
+        if ( c->IsKind( STANDARD_TYPE(Geom_BSplineCurve )) ||
+             c->IsKind( STANDARD_TYPE(Geom_BezierCurve )))
+          edgeIsSafe = false;
+      }
+    }
+    // check 2d curve
+    if ( edgeIsSafe )
+    {
+      Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( e, surf, loc, f, l);
+      if ( !c2.IsNull() )
+      {
+        Handle(Geom2d_TrimmedCurve) tc = Handle(Geom2d_TrimmedCurve)::DownCast(c2);
+        while( !tc.IsNull() ) {
+          c2 = tc->BasisCurve();
+          tc = Handle(Geom2d_TrimmedCurve)::DownCast(c2);
+        }
+        if ( c2->IsKind( STANDARD_TYPE(Geom2d_BSplineCurve )) ||
+             c2->IsKind( STANDARD_TYPE(Geom2d_BezierCurve )))
+          edgeIsSafe = false;
+      }
+    }
+    if ( !edgeIsSafe && !noSafeTShapes.insert( e.TShape().get() ).second )
+      isSafe = false;
+  }
+  return isSafe;
+}
+#endif
+
+//================================================================================
+/*
+ * Store an intersection if it is IN or ON the face
+ */
+void FaceLineIntersector::addIntPoint(const bool toClassify)
+{
+  if ( !toClassify || UVIsOnFace() )
+  {
+    F_IntersectPoint p;
+    p._paramOnLine = _w;
+    p._u           = _u;
+    p._v           = _v;
+    p._transition  = _transition;
+    _intPoints.push_back( p );
+  }
+}
+
+//================================================================================
+/*
+ * Intersect a line with a plane
+ */
+void FaceLineIntersector::IntersectWithPlane(const GridLine& gridLine)
+{
+  IntAna_IntConicQuad linPlane( gridLine._line, _plane, Precision::Angular());
+  _w = linPlane.ParamOnConic(1);
+  if ( isParamOnLineOK( gridLine._length ))
+  {
+    ElSLib::Parameters(_plane, linPlane.Point(1) ,_u,_v);
+    addIntPoint();
+  }
+}
+
+//================================================================================
+/*
+ * Intersect a line with a cylinder
+ */
+void FaceLineIntersector::IntersectWithCylinder(const GridLine& gridLine)
+{
+  IntAna_IntConicQuad linCylinder( gridLine._line, _cylinder );
+  if ( linCylinder.IsDone() && linCylinder.NbPoints() > 0 )
+  {
+    _w = linCylinder.ParamOnConic(1);
+    if ( linCylinder.NbPoints() == 1 )
+      _transition = Trans_TANGENT;
+    else
+      _transition = _w < linCylinder.ParamOnConic(2) ? _transIn : _transOut;
+    if ( isParamOnLineOK( gridLine._length ))
+    {
+      ElSLib::Parameters(_cylinder, linCylinder.Point(1) ,_u,_v);
+      addIntPoint();
+    }
+    if ( linCylinder.NbPoints() > 1 )
+    {
+      _w = linCylinder.ParamOnConic(2);
+      if ( isParamOnLineOK( gridLine._length ))
+      {
+        ElSLib::Parameters(_cylinder, linCylinder.Point(2) ,_u,_v);
+        _transition = ( _transition == Trans_OUT ) ? Trans_IN : Trans_OUT;
+        addIntPoint();
+      }
+    }
+  }
+}
+
+//================================================================================
+/*
+ * Intersect a line with a cone
+ */
+void FaceLineIntersector::IntersectWithCone (const GridLine& gridLine)
+{
+  IntAna_IntConicQuad linCone(gridLine._line,_cone);
+  if ( !linCone.IsDone() ) return;
+  gp_Pnt P;
+  gp_Vec du, dv, norm;
+  for ( int i = 1; i <= linCone.NbPoints(); ++i )
+  {
+    _w = linCone.ParamOnConic( i );
+    if ( !isParamOnLineOK( gridLine._length )) continue;
+    ElSLib::Parameters(_cone, linCone.Point(i) ,_u,_v);
+    if ( UVIsOnFace() )
+    {
+      ElSLib::D1( _u, _v, _cone, P, du, dv );
+      norm = du ^ dv;
+      double normSize2 = norm.SquareMagnitude();
+      if ( normSize2 > Precision::Angular() * Precision::Angular() )
+      {
+        double cos = norm.XYZ() * gridLine._line.Direction().XYZ();
+        cos /= sqrt( normSize2 );
+        if ( cos < -Precision::Angular() )
+          _transition = _transIn;
+        else if ( cos > Precision::Angular() )
+          _transition = _transOut;
+        else
+          _transition = Trans_TANGENT;
+      }
+      else
+      {
+        _transition = Trans_APEX;
+      }
+      addIntPoint( /*toClassify=*/false);
+    }
+  }
+}
+
+//================================================================================
+/*
+ * Intersect a line with a sphere
+ */
+void FaceLineIntersector::IntersectWithSphere  (const GridLine& gridLine)
+{
+  IntAna_IntConicQuad linSphere(gridLine._line,_sphere);
+  if ( linSphere.IsDone() && linSphere.NbPoints() > 0 )
+  {
+    _w = linSphere.ParamOnConic(1);
+    if ( linSphere.NbPoints() == 1 )
+      _transition = Trans_TANGENT;
+    else
+      _transition = _w < linSphere.ParamOnConic(2) ? _transIn : _transOut;
+    if ( isParamOnLineOK( gridLine._length ))
+    {
+      ElSLib::Parameters(_sphere, linSphere.Point(1) ,_u,_v);
+      addIntPoint();
+    }
+    if ( linSphere.NbPoints() > 1 )
+    {
+      _w = linSphere.ParamOnConic(2);
+      if ( isParamOnLineOK( gridLine._length ))
+      {
+        ElSLib::Parameters(_sphere, linSphere.Point(2) ,_u,_v);
+        _transition = ( _transition == Trans_OUT ) ? Trans_IN : Trans_OUT;
+        addIntPoint();
+      }
+    }
+  }
+}
+
+//================================================================================
+/*
+ * Intersect a line with a torus
+ */
+void FaceLineIntersector::IntersectWithTorus (const GridLine& gridLine)
+{
+  IntAna_IntLinTorus linTorus(gridLine._line,_torus);
+  if ( !linTorus.IsDone()) return;
+  gp_Pnt P;
+  gp_Vec du, dv, norm;
+  for ( int i = 1; i <= linTorus.NbPoints(); ++i )
+  {
+    _w = linTorus.ParamOnLine( i );
+    if ( !isParamOnLineOK( gridLine._length )) continue;
+    linTorus.ParamOnTorus( i, _u,_v );
+    if ( UVIsOnFace() )
+    {
+      ElSLib::D1( _u, _v, _torus, P, du, dv );
+      norm = du ^ dv;
+      double normSize = norm.Magnitude();
+      double cos = norm.XYZ() * gridLine._line.Direction().XYZ();
+      cos /= normSize;
+      if ( cos < -Precision::Angular() )
+        _transition = _transIn;
+      else if ( cos > Precision::Angular() )
+        _transition = _transOut;
+      else
+        _transition = Trans_TANGENT;
+      addIntPoint( /*toClassify=*/false);
+    }
+  }
+}
+
+//================================================================================
+/*
+ * Intersect a line with a non-analytical surface
+ */
+void FaceLineIntersector::IntersectWithSurface (const GridLine& gridLine)
+{
+  _surfaceInt->Perform( gridLine._line, 0.0, gridLine._length );
+  if ( !_surfaceInt->IsDone() ) return;
+  for ( int i = 1; i <= _surfaceInt->NbPnt(); ++i )
+  {
+    _transition = Transition( _surfaceInt->Transition( i ) );
+    _w = _surfaceInt->WParameter( i );
+    addIntPoint(/*toClassify=*/false);
+  }
+}
index 3f8cc68007ebd80a98e1ff39743ba0481becc9d1..3cf439aa3434b27f18610d5f90579c015eac82e9 100644 (file)
 
 // STD
 #include <algorithm>
+#include <map>
 #include <vector>
 #include <memory>
 #include <mutex>
 #include <thread>
 
-
 // SMESH
 #include "SMESH_StdMeshers.hxx"
 #include "StdMeshers_FaceSide.hxx"
 
 // All utility structs used in Grid and hexahedron class will be included here
 // Ideally each one of this should define their own testable class
-using namespace std;
-using namespace SMESH;
-namespace gridtools
+namespace StdMeshers
+{
+namespace Cartesian3D
 {
   typedef int                     TGeomID; // IDs of sub-shapes
   typedef TopTools_ShapeMapHasher TShapeHasher; // non-oriented shape hasher
   typedef std::array< int, 3 >    TIJK;
 
+  typedef std::map< TGeomID, std::vector< TGeomID > > TEdge2faceIDsMap;
+
   const TGeomID theUndefID = 1e+9;
 
   //=============================================================================
@@ -136,6 +138,7 @@ namespace gridtools
     Trans_APEX,
     Trans_INTERNAL // for INTERNAL FACE
   };
+
   // --------------------------------------------------------------------------
   /*!
    * \brief Sub-entities of a FACE neighboring its concave VERTEX.
@@ -155,6 +158,7 @@ namespace gridtools
     void SetVertex( TGeomID v  ) { ( _v1 ? _v2 : _v1 ) = v; }
   };
   typedef NCollection_DataMap< TGeomID, ConcaveFace > TConcaveVertex2Face;
+
   // --------------------------------------------------------------------------
   /*!
    * \brief Container of IDs of SOLID sub-shapes
@@ -167,7 +171,7 @@ namespace gridtools
   public:
     virtual ~Solid() {}
     virtual bool Contains( TGeomID /*subID*/ ) const { return true; }
-    virtual bool ContainsAny( const vector< TGeomID>& /*subIDs*/ ) const { return true; }
+    virtual bool ContainsAny( const std::vector< TGeomID>& /*subIDs*/ ) const { return true; }
     virtual TopAbs_Orientation Orientation( const TopoDS_Shape& s ) const { return s.Orientation(); }
     virtual bool IsOutsideOriented( TGeomID /*faceID*/ ) const { return true; }
     void SetID( TGeomID id ) { _id = id; }
@@ -179,6 +183,7 @@ namespace gridtools
     bool HasConcaveVertex() const { return !_concaveVertex.IsEmpty(); }
     const ConcaveFace* GetConcave( TGeomID V ) const { return _concaveVertex.Seek( V ); }
   };
+
   // --------------------------------------------------------------------------
   class OneOfSolids : public Solid
   {
@@ -190,7 +195,7 @@ namespace gridtools
                TopAbs_ShapeEnum    subType,
                const SMESHDS_Mesh* mesh );
     virtual bool Contains( TGeomID i ) const { return i == ID() || _subIDs.Contains( i ); }
-    virtual bool ContainsAny( const vector< TGeomID>& subIDs ) const
+    virtual bool ContainsAny( const std::vector< TGeomID>& subIDs ) const
     {
       for ( size_t i = 0; i < subIDs.size(); ++i ) if ( Contains( subIDs[ i ])) return true;
       return false;
@@ -205,6 +210,7 @@ namespace gridtools
       return faceID == 0 || _outFaceIDs.Contains( faceID );
     }
   };
+
   // --------------------------------------------------------------------------
   /*!
    * \brief Hold a vector of TGeomID and clear it at destruction
@@ -246,6 +252,7 @@ namespace gridtools
       return common;
     }
   };
+
   // --------------------------------------------------------------------------
   /*!
    * \brief Geom data
@@ -253,20 +260,20 @@ namespace gridtools
   struct Geometry
   {
     TopoDS_Shape                _mainShape;
-    vector< vector< TGeomID > > _solidIDsByShapeID;// V/E/F ID -> SOLID IDs
-    Solid                       _soleSolid;
-    map< TGeomID, OneOfSolids > _solidByID;
+    std::vector< std::vector< TGeomID > > _solidIDsByShapeID;// V/E/F ID -> SOLID IDs
+    Solid                            _soleSolid;
+    std::map< TGeomID, OneOfSolids > _solidByID;
     TColStd_MapOfInteger        _boundaryFaces; // FACEs on boundary of mesh->ShapeToMesh()
     TColStd_MapOfInteger        _strangeEdges; // EDGEs shared by strange FACEs
     TGeomID                     _extIntFaceID; // pseudo FACE - extension of INTERNAL FACE
 
     TopTools_DataMapOfShapeInteger _shape2NbNodes; // nb of pre-existing nodes on shapes
 
-    Controls::ElementsOnShape _edgeClassifier;
-    Controls::ElementsOnShape _vertexClassifier;
+    SMESH::Controls::ElementsOnShape _edgeClassifier;
+    SMESH::Controls::ElementsOnShape _vertexClassifier;
 
     bool IsOneSolid() const { return _solidByID.size() < 2; }
-    GeomIDVecHelder GetSolidIDsByShapeID( const vector< TGeomID >& shapeIDs ) const;
+    GeomIDVecHelder GetSolidIDsByShapeID( const std::vector< TGeomID >& shapeIDs ) const;
   };
 
   // --------------------------------------------------------------------------
@@ -279,16 +286,17 @@ namespace gridtools
     // See Add method modify _node and _faceIDs class members dinamicaly during execution
     // of Hexahedron.compute() method.
     // std::mutex _mutex;
-    mutable const SMDS_MeshNode* _node;
-    mutable vector< TGeomID >    _faceIDs;
+    mutable const SMDS_MeshNode*   _node;
+    mutable std::vector< TGeomID > _faceIDs;
 
     B_IntersectPoint(): _node(NULL) {}
-    bool Add( const vector< TGeomID >& fIDs, const SMDS_MeshNode* n=NULL ) const;
+    bool Add( const std::vector< TGeomID >& fIDs, const SMDS_MeshNode* n=NULL ) const;
     TGeomID HasCommonFace( const B_IntersectPoint * other, TGeomID avoidFace=-1 ) const;
     size_t GetCommonFaces( const B_IntersectPoint * other, TGeomID * commonFaces ) const;
     bool IsOnFace( TGeomID faceID ) const;
     virtual ~B_IntersectPoint() {}
   };
+
   // --------------------------------------------------------------------------
   /*!
    * \brief Data of intersection between a GridLine and a TopoDS_Face
@@ -304,6 +312,7 @@ namespace gridtools
       return _paramOnLine < o._paramOnLine;         
     }
   };
+
   // --------------------------------------------------------------------------
   /*!
    * \brief Data of intersection between GridPlanes and a TopoDS_EDGE
@@ -323,10 +332,10 @@ namespace gridtools
   {
     gp_Lin _line;
     double _length; // line length
-    multiset< F_IntersectPoint > _intPoints;
+    std::multiset< F_IntersectPoint > _intPoints;
 
     void RemoveExcessIntPoints( const double tol );
-    TGeomID GetSolidIDBefore( multiset< F_IntersectPoint >::iterator ip,
+    TGeomID GetSolidIDBefore( std::multiset< F_IntersectPoint >::iterator ip,
                               const TGeomID                          prevID,
                               const Geometry&                        geom);
   };
@@ -336,9 +345,9 @@ namespace gridtools
    */
   struct GridPlanes
   {
-    gp_XYZ           _zNorm;
-    vector< gp_XYZ > _origins; // origin points of all planes in one direction
-    vector< double > _zProjs;  // projections of origins to _zNorm
+    gp_XYZ                _zNorm;
+    std::vector< gp_XYZ > _origins; // origin points of all planes in one direction
+    std::vector< double > _zProjs;  // projections of origins to _zNorm
   };
   // --------------------------------------------------------------------------
   /*!
@@ -349,11 +358,11 @@ namespace gridtools
     size_t _size  [3];
     size_t _curInd[3];
     size_t _iVar1, _iVar2, _iConst;
-    string _name1, _name2, _nameConst;
+    std::string _name1, _name2, _nameConst;
     LineIndexer() {}
     LineIndexer( size_t sz1, size_t sz2, size_t sz3,
                  size_t iv1, size_t iv2, size_t iConst,
-                 const string& nv1, const string& nv2, const string& nConst )
+                 const std::string& nv1, const std::string& nv2, const std::string& nConst )
     {
       _size[0] = sz1; _size[1] = sz2; _size[2] = sz3;
       _curInd[0] = _curInd[1] = _curInd[2] = 0;
@@ -398,31 +407,28 @@ namespace gridtools
       * \brief computes exact bounding box with axes parallel to given ones
       */
       //================================================================================
-      static void GetExactBndBox( const vector< TopoDS_Shape >& faceVec, const double* axesDirs, Bnd_Box& shapeBox );
+      static void GetExactBndBox( const std::vector< TopoDS_Shape >& faceVec, const double* axesDirs, Bnd_Box& shapeBox );
   };
-} // end namespace gridtools
-
-using namespace gridtools;  
-class STDMESHERS_EXPORT Grid
-{
 
-public:
-    vector< double >   _coords[3]; // coordinates of grid nodes
-    gp_XYZ             _axes  [3]; // axis directions
-    vector< GridLine > _lines [3]; //    in 3 directions
-    double             _tol, _minCellSize;
-    gp_XYZ             _origin;
-    gp_Mat             _invB; // inverted basis of _axes
+  class STDMESHERS_EXPORT Grid
+  {
+    public:
+    std::vector< double >   _coords[3]; // coordinates of grid nodes
+    gp_XYZ                  _axes  [3]; // axis directions
+    std::vector< GridLine > _lines [3]; //    in 3 directions
+    double                  _tol, _minCellSize;
+    gp_XYZ                 _origin;
+    gp_Mat                 _invB; // inverted basis of _axes
 
     // index shift within _nodes of nodes of a cell from the 1st node
-    int                _nodeShift[8];
+    int                    _nodeShift[8];
 
-    vector< const SMDS_MeshNode* >    _nodes;          // mesh nodes at grid nodes
-    vector< const SMDS_MeshNode* >    _allBorderNodes; // mesh nodes between the bounding box and the geometry boundary
+    std::vector< const SMDS_MeshNode* >    _nodes;          // mesh nodes at grid nodes
+    std::vector< const SMDS_MeshNode* >    _allBorderNodes; // mesh nodes between the bounding box and the geometry boundary
 
-    vector< const F_IntersectPoint* > _gridIntP; // grid node intersection with geometry
-    ObjectPool< E_IntersectPoint >    _edgeIntPool; // intersections with EDGEs
-    ObjectPool< F_IntersectPoint >    _extIntPool; // intersections with extended INTERNAL FACEs
+    std::vector< const F_IntersectPoint* > _gridIntP; // grid node intersection with geometry
+    ObjectPool< E_IntersectPoint >        _edgeIntPool; // intersections with EDGEs
+    ObjectPool< F_IntersectPoint >        _extIntPool; // intersections with extended INTERNAL FACEs
     //list< E_IntersectPoint >          _edgeIntP; // intersections with EDGEs
 
     Geometry                          _geometry;
@@ -469,10 +475,10 @@ public:
     void InitGeometry( const TopoDS_Shape& theShape );
     void InitClassifier( const TopoDS_Shape&        mainShape,
                         TopAbs_ShapeEnum           shapeType,
-                        Controls::ElementsOnShape& classifier );
-    void GetEdgesToImplement( map< TGeomID, vector< TGeomID > > & edge2faceMap,
+                        SMESH::Controls::ElementsOnShape& classifier );
+    void GetEdgesToImplement( std::map< TGeomID, std::vector< TGeomID > > & edge2faceMap,
                               const TopoDS_Shape&                 shape,
-                              const vector< TopoDS_Shape >&       faces );
+                              const std::vector< TopoDS_Shape >&       faces );
     void SetSolidFather( const TopoDS_Shape& s, const TopoDS_Shape& theShapeToMesh );
     bool IsShared( TGeomID faceID ) const;
     bool IsAnyShared( const std::vector< TGeomID >& faceIDs ) const;
@@ -486,7 +492,7 @@ public:
     TGeomID PseudoIntExtFaceID() const { return _geometry._extIntFaceID; }
     Solid* GetSolid( TGeomID solidID = 0 );
     Solid* GetOneOfSolids( TGeomID solidID );
-    const vector< TGeomID > & GetSolidIDs( TGeomID subShapeID ) const;
+    const std::vector< TGeomID > & GetSolidIDs( TGeomID subShapeID ) const;
     bool IsCorrectTransition( TGeomID faceID, const Solid* solid );
     bool IsBoundaryFace( TGeomID face ) const { return _geometry._boundaryFaces.Contains( face ); }
     void SetOnShape( const SMDS_MeshNode* n, const F_IntersectPoint& ip,
@@ -495,19 +501,19 @@ public:
     bool IsToCheckNodePos() const { return !_toAddEdges && _toCreateFaces; }
     bool IsToRemoveExcessEntities() const { return !_toAddEdges; }
 
-    void SetCoordinates(const vector<double>& xCoords,
-                        const vector<double>& yCoords,
-                        const vector<double>& zCoords,
+    void SetCoordinates(const std::vector<double>& xCoords,
+                        const std::vector<double>& yCoords,
+                        const std::vector<double>& zCoords,
                         const double*         axesDirs,
                         const Bnd_Box&        bndBox );
     void ComputeUVW(const gp_XYZ& p, double uvw[3]);
     void ComputeNodes(SMESH_MesherHelper& helper);
-    bool GridInitAndInterserctWithShape( const TopoDS_Shape& theShape, std::map< TGeomID, vector< TGeomID > >& edge2faceIDsMap, 
-                                          const StdMeshers_CartesianParameters3D* hyp, const int numOfThreads, bool computeCanceled );
-};
-
-namespace
-{
+    bool GridInitAndInterserctWithShape( const TopoDS_Shape& theShape,
+                                         std::map< TGeomID, std::vector< TGeomID > >& edge2faceIDsMap, 
+                                         const StdMeshers_CartesianParameters3D* hyp,
+                                         const int numOfThreads,
+                                         bool computeCanceled );
+  };
 
   // Implement parallel computation of Hexa with c++ thread implementation
   template<typename Iterator, class Function>
@@ -529,6 +535,7 @@ namespace
       std::for_each(it, last, f); // last steps while we wait for other threads
       std::for_each(threads.begin(), threads.end(), [](std::thread& x){x.join();});
   }
+
   // --------------------------------------------------------------------------
   /*!
    * \brief Intersector of TopoDS_Face with all GridLine's
@@ -540,7 +547,7 @@ namespace
     Grid*       _grid;
     Bnd_Box     _bndBox;
     IntCurvesFace_Intersector* _surfaceInt;
-    vector< std::pair< GridLine*, F_IntersectPoint > > _intersections;
+    std::vector< std::pair< GridLine*, F_IntersectPoint > > _intersections;
 
     FaceGridIntersector(): _grid(0), _surfaceInt(0) {}
     void Intersect();
@@ -549,7 +556,7 @@ namespace
     {
       for ( size_t i = 0; i < _intersections.size(); ++i )
       {
-        multiset< F_IntersectPoint >::iterator ip =
+        std::multiset< F_IntersectPoint >::iterator ip =
           _intersections[i].first->_intPoints.insert( _intersections[i].second );
         ip->_faceIDs.reserve( 1 );
         ip->_faceIDs.push_back( _faceID );
@@ -572,7 +579,7 @@ namespace
       return _surfaceInt;
     }
 #ifdef WITH_TBB    
-    bool IsThreadSafe(set< const Standard_Transient* >& noSafeTShapes) const;
+    bool IsThreadSafe(std::set< const Standard_Transient* >& noSafeTShapes) const;
 #endif
   };
   
@@ -583,8 +590,8 @@ namespace
    */
   struct ParallelIntersector
   {
-    vector< FaceGridIntersector >& _faceVec;
-    ParallelIntersector( vector< FaceGridIntersector >& faceVec): _faceVec(faceVec){}
+    std::vector< FaceGridIntersector >& _faceVec;
+    ParallelIntersector( std::vector< FaceGridIntersector >& faceVec): _faceVec(faceVec){}
     void operator() ( const tbb::blocked_range<size_t>& r ) const
     {
       for ( size_t i = r.begin(); i != r.end(); ++i )
@@ -617,7 +624,7 @@ namespace
     gp_Torus    _torus;
     IntCurvesFace_Intersector* _surfaceInt;
 
-    vector< F_IntersectPoint > _intPoints;
+    std::vector< F_IntersectPoint > _intPoints;
 
     void IntersectWithPlane   (const GridLine& gridLine);
     void IntersectWithCylinder(const GridLine& gridLine);
@@ -626,7 +633,14 @@ namespace
     void IntersectWithTorus   (const GridLine& gridLine);
     void IntersectWithSurface (const GridLine& gridLine);
 
-    bool UVIsOnFace() const;
+    /*
+     * Return true if (_u,_v) is on the face
+     */
+    bool UVIsOnFace() const
+    {
+      TopAbs_State state = _surfaceInt->ClassifyUVPoint(gp_Pnt2d( _u,_v ));
+      return ( state == TopAbs_IN || state == TopAbs_ON );
+    }
     void addIntPoint(const bool toClassify=true);
     bool isParamOnLineOK( const double linLength )
     {
@@ -636,356 +650,7 @@ namespace
     ~FaceLineIntersector() { if (_surfaceInt ) delete _surfaceInt; _surfaceInt = 0; }
   };
 
-    //=============================================================================
-  /*
-   * Intersects TopoDS_Face with all GridLine's
-   */
-  void FaceGridIntersector::Intersect()
-  {
-    FaceLineIntersector intersector;
-    intersector._surfaceInt = GetCurveFaceIntersector();
-    intersector._tol        = _grid->_tol;
-    intersector._transOut   = _face.Orientation() == TopAbs_REVERSED ? Trans_IN : Trans_OUT;
-    intersector._transIn    = _face.Orientation() == TopAbs_REVERSED ? Trans_OUT : Trans_IN;
-
-    typedef void (FaceLineIntersector::* PIntFun )(const GridLine& gridLine);
-    PIntFun interFunction;
-
-    bool isDirect = true;
-    BRepAdaptor_Surface surf( _face );
-    switch ( surf.GetType() ) {
-    case GeomAbs_Plane:
-      intersector._plane = surf.Plane();
-      interFunction = &FaceLineIntersector::IntersectWithPlane;
-      isDirect = intersector._plane.Direct();
-      break;
-    case GeomAbs_Cylinder:
-      intersector._cylinder = surf.Cylinder();
-      interFunction = &FaceLineIntersector::IntersectWithCylinder;
-      isDirect = intersector._cylinder.Direct();
-      break;
-    case GeomAbs_Cone:
-      intersector._cone = surf.Cone();
-      interFunction = &FaceLineIntersector::IntersectWithCone;
-      //isDirect = intersector._cone.Direct();
-      break;
-    case GeomAbs_Sphere:
-      intersector._sphere = surf.Sphere();
-      interFunction = &FaceLineIntersector::IntersectWithSphere;
-      isDirect = intersector._sphere.Direct();
-      break;
-    case GeomAbs_Torus:
-      intersector._torus = surf.Torus();
-      interFunction = &FaceLineIntersector::IntersectWithTorus;
-      //isDirect = intersector._torus.Direct();
-      break;
-    default:
-      interFunction = &FaceLineIntersector::IntersectWithSurface;
-    }
-    if ( !isDirect )
-      std::swap( intersector._transOut, intersector._transIn );
-
-    _intersections.clear();
-    for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions
-    {
-      if ( surf.GetType() == GeomAbs_Plane )
-      {
-        // check if all lines in this direction are parallel to a plane
-        if ( intersector._plane.Axis().IsNormal( _grid->_lines[iDir][0]._line.Position(),
-                                                 Precision::Angular()))
-          continue;
-        // find out a transition, that is the same for all lines of a direction
-        gp_Dir plnNorm = intersector._plane.Axis().Direction();
-        gp_Dir lineDir = _grid->_lines[iDir][0]._line.Direction();
-        intersector._transition =
-          ( plnNorm * lineDir < 0 ) ? intersector._transIn : intersector._transOut;
-      }
-      if ( surf.GetType() == GeomAbs_Cylinder )
-      {
-        // check if all lines in this direction are parallel to a cylinder
-        if ( intersector._cylinder.Axis().IsParallel( _grid->_lines[iDir][0]._line.Position(),
-                                                      Precision::Angular()))
-          continue;
-      }
-
-      // intersect the grid lines with the face
-      for ( size_t iL = 0; iL < _grid->_lines[iDir].size(); ++iL )
-      {
-        GridLine& gridLine = _grid->_lines[iDir][iL];
-        if ( _bndBox.IsOut( gridLine._line )) continue;
-
-        intersector._intPoints.clear();
-        (intersector.*interFunction)( gridLine ); // <- intersection with gridLine
-        for ( size_t i = 0; i < intersector._intPoints.size(); ++i )
-          _intersections.push_back( make_pair( &gridLine, intersector._intPoints[i] ));
-      }
-    }
-
-    if ( _face.Orientation() == TopAbs_INTERNAL )
-    {
-      for ( size_t i = 0; i < _intersections.size(); ++i )
-        if ( _intersections[i].second._transition == Trans_IN ||
-             _intersections[i].second._transition == Trans_OUT )
-        {
-          _intersections[i].second._transition = Trans_INTERNAL;
-        }
-    }
-    return;
-  }
-  //================================================================================
-  /*
-   * Return true if (_u,_v) is on the face
-   */
-  bool FaceLineIntersector::UVIsOnFace() const
-  {
-    TopAbs_State state = _surfaceInt->ClassifyUVPoint(gp_Pnt2d( _u,_v ));
-    return ( state == TopAbs_IN || state == TopAbs_ON );
-  }
-  //================================================================================
-  /*
-   * Store an intersection if it is IN or ON the face
-   */
-  void FaceLineIntersector::addIntPoint(const bool toClassify)
-  {
-    if ( !toClassify || UVIsOnFace() )
-    {
-      F_IntersectPoint p;
-      p._paramOnLine = _w;
-      p._u           = _u;
-      p._v           = _v;
-      p._transition  = _transition;
-      _intPoints.push_back( p );
-    }
-  }
-  //================================================================================
-  /*
-   * Intersect a line with a plane
-   */
-  void FaceLineIntersector::IntersectWithPlane(const GridLine& gridLine)
-  {
-    IntAna_IntConicQuad linPlane( gridLine._line, _plane, Precision::Angular());
-    _w = linPlane.ParamOnConic(1);
-    if ( isParamOnLineOK( gridLine._length ))
-    {
-      ElSLib::Parameters(_plane, linPlane.Point(1) ,_u,_v);
-      addIntPoint();
-    }
-  }
-  //================================================================================
-  /*
-   * Intersect a line with a cylinder
-   */
-  void FaceLineIntersector::IntersectWithCylinder(const GridLine& gridLine)
-  {
-    IntAna_IntConicQuad linCylinder( gridLine._line, _cylinder );
-    if ( linCylinder.IsDone() && linCylinder.NbPoints() > 0 )
-    {
-      _w = linCylinder.ParamOnConic(1);
-      if ( linCylinder.NbPoints() == 1 )
-        _transition = Trans_TANGENT;
-      else
-        _transition = _w < linCylinder.ParamOnConic(2) ? _transIn : _transOut;
-      if ( isParamOnLineOK( gridLine._length ))
-      {
-        ElSLib::Parameters(_cylinder, linCylinder.Point(1) ,_u,_v);
-        addIntPoint();
-      }
-      if ( linCylinder.NbPoints() > 1 )
-      {
-        _w = linCylinder.ParamOnConic(2);
-        if ( isParamOnLineOK( gridLine._length ))
-        {
-          ElSLib::Parameters(_cylinder, linCylinder.Point(2) ,_u,_v);
-          _transition = ( _transition == Trans_OUT ) ? Trans_IN : Trans_OUT;
-          addIntPoint();
-        }
-      }
-    }
-  }
-  //================================================================================
-  /*
-   * Intersect a line with a cone
-   */
-  void FaceLineIntersector::IntersectWithCone (const GridLine& gridLine)
-  {
-    IntAna_IntConicQuad linCone(gridLine._line,_cone);
-    if ( !linCone.IsDone() ) return;
-    gp_Pnt P;
-    gp_Vec du, dv, norm;
-    for ( int i = 1; i <= linCone.NbPoints(); ++i )
-    {
-      _w = linCone.ParamOnConic( i );
-      if ( !isParamOnLineOK( gridLine._length )) continue;
-      ElSLib::Parameters(_cone, linCone.Point(i) ,_u,_v);
-      if ( UVIsOnFace() )
-      {
-        ElSLib::D1( _u, _v, _cone, P, du, dv );
-        norm = du ^ dv;
-        double normSize2 = norm.SquareMagnitude();
-        if ( normSize2 > Precision::Angular() * Precision::Angular() )
-        {
-          double cos = norm.XYZ() * gridLine._line.Direction().XYZ();
-          cos /= sqrt( normSize2 );
-          if ( cos < -Precision::Angular() )
-            _transition = _transIn;
-          else if ( cos > Precision::Angular() )
-            _transition = _transOut;
-          else
-            _transition = Trans_TANGENT;
-        }
-        else
-        {
-          _transition = Trans_APEX;
-        }
-        addIntPoint( /*toClassify=*/false);
-      }
-    }
-  }
-  //================================================================================
-  /*
-   * Intersect a line with a sphere
-   */
-  void FaceLineIntersector::IntersectWithSphere  (const GridLine& gridLine)
-  {
-    IntAna_IntConicQuad linSphere(gridLine._line,_sphere);
-    if ( linSphere.IsDone() && linSphere.NbPoints() > 0 )
-    {
-      _w = linSphere.ParamOnConic(1);
-      if ( linSphere.NbPoints() == 1 )
-        _transition = Trans_TANGENT;
-      else
-        _transition = _w < linSphere.ParamOnConic(2) ? _transIn : _transOut;
-      if ( isParamOnLineOK( gridLine._length ))
-      {
-        ElSLib::Parameters(_sphere, linSphere.Point(1) ,_u,_v);
-        addIntPoint();
-      }
-      if ( linSphere.NbPoints() > 1 )
-      {
-        _w = linSphere.ParamOnConic(2);
-        if ( isParamOnLineOK( gridLine._length ))
-        {
-          ElSLib::Parameters(_sphere, linSphere.Point(2) ,_u,_v);
-          _transition = ( _transition == Trans_OUT ) ? Trans_IN : Trans_OUT;
-          addIntPoint();
-        }
-      }
-    }
-  }
-  //================================================================================
-  /*
-   * Intersect a line with a torus
-   */
-  void FaceLineIntersector::IntersectWithTorus   (const GridLine& gridLine)
-  {
-    IntAna_IntLinTorus linTorus(gridLine._line,_torus);
-    if ( !linTorus.IsDone()) return;
-    gp_Pnt P;
-    gp_Vec du, dv, norm;
-    for ( int i = 1; i <= linTorus.NbPoints(); ++i )
-    {
-      _w = linTorus.ParamOnLine( i );
-      if ( !isParamOnLineOK( gridLine._length )) continue;
-      linTorus.ParamOnTorus( i, _u,_v );
-      if ( UVIsOnFace() )
-      {
-        ElSLib::D1( _u, _v, _torus, P, du, dv );
-        norm = du ^ dv;
-        double normSize = norm.Magnitude();
-        double cos = norm.XYZ() * gridLine._line.Direction().XYZ();
-        cos /= normSize;
-        if ( cos < -Precision::Angular() )
-          _transition = _transIn;
-        else if ( cos > Precision::Angular() )
-          _transition = _transOut;
-        else
-          _transition = Trans_TANGENT;
-        addIntPoint( /*toClassify=*/false);
-      }
-    }
-  }
-  //================================================================================
-  /*
-   * Intersect a line with a non-analytical surface
-   */
-  void FaceLineIntersector::IntersectWithSurface (const GridLine& gridLine)
-  {
-    _surfaceInt->Perform( gridLine._line, 0.0, gridLine._length );
-    if ( !_surfaceInt->IsDone() ) return;
-    for ( int i = 1; i <= _surfaceInt->NbPnt(); ++i )
-    {
-      _transition = Transition( _surfaceInt->Transition( i ) );
-      _w = _surfaceInt->WParameter( i );
-      addIntPoint(/*toClassify=*/false);
-    }
-  }
-
-#ifdef WITH_TBB
-  //================================================================================
-  /*
-   * check if its face can be safely intersected in a thread
-   */
-  bool FaceGridIntersector::IsThreadSafe(set< const Standard_Transient* >& noSafeTShapes) const
-  {
-    bool isSafe = true;
-
-    // check surface
-    TopLoc_Location loc;
-    Handle(Geom_Surface) surf = BRep_Tool::Surface( _face, loc );
-    Handle(Geom_RectangularTrimmedSurface) ts =
-      Handle(Geom_RectangularTrimmedSurface)::DownCast( surf );
-    while( !ts.IsNull() ) {
-      surf = ts->BasisSurface();
-      ts = Handle(Geom_RectangularTrimmedSurface)::DownCast(surf);
-    }
-    if ( surf->IsKind( STANDARD_TYPE(Geom_BSplineSurface )) ||
-         surf->IsKind( STANDARD_TYPE(Geom_BezierSurface )))
-      if ( !noSafeTShapes.insert( _face.TShape().get() ).second )
-        isSafe = false;
-
-    double f, l;
-    TopExp_Explorer exp( _face, TopAbs_EDGE );
-    for ( ; exp.More(); exp.Next() )
-    {
-      bool edgeIsSafe = true;
-      const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
-      // check 3d curve
-      {
-        Handle(Geom_Curve) c = BRep_Tool::Curve( e, loc, f, l);
-        if ( !c.IsNull() )
-        {
-          Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(c);
-          while( !tc.IsNull() ) {
-            c = tc->BasisCurve();
-            tc = Handle(Geom_TrimmedCurve)::DownCast(c);
-          }
-          if ( c->IsKind( STANDARD_TYPE(Geom_BSplineCurve )) ||
-               c->IsKind( STANDARD_TYPE(Geom_BezierCurve )))
-            edgeIsSafe = false;
-        }
-      }
-      // check 2d curve
-      if ( edgeIsSafe )
-      {
-        Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( e, surf, loc, f, l);
-        if ( !c2.IsNull() )
-        {
-          Handle(Geom2d_TrimmedCurve) tc = Handle(Geom2d_TrimmedCurve)::DownCast(c2);
-          while( !tc.IsNull() ) {
-            c2 = tc->BasisCurve();
-            tc = Handle(Geom2d_TrimmedCurve)::DownCast(c2);
-          }
-          if ( c2->IsKind( STANDARD_TYPE(Geom2d_BSplineCurve )) ||
-               c2->IsKind( STANDARD_TYPE(Geom2d_BezierCurve )))
-            edgeIsSafe = false;
-        }
-      }
-      if ( !edgeIsSafe && !noSafeTShapes.insert( e.TShape().get() ).second )
-        isSafe = false;
-    }
-    return isSafe;
-  }
-#endif
-}
+} // end namespace Cartesian3D
+} // end namespace StdMeshers
 
 #endif
index f1f8f6965a79a6f081248a164b4e9d7a5812cfa1..8a2b103008a9744ec7784dbfe683c61d413443dc 100644 (file)
@@ -23,8 +23,9 @@
 
 #include "StdMeshers_Cartesian_3D_Hexahedron.hxx"
 
+using namespace std;
 using namespace SMESH;
-using namespace gridtools;
+using namespace StdMeshers::Cartesian3D;
 
 std::mutex _eMutex;
 
@@ -856,8 +857,31 @@ void Hexahedron::computeElements( const Solid* solid, int solidIndex )
 /*!
   * \brief Compute mesh volumes resulted from intersection of the Hexahedron
   */
-void Hexahedron::defineHexahedralFaces( std::vector< _OrientedLink >& splits, std::vector<_Node*>& chainNodes, std::set< TGeomID >& concaveFaces, bool toCheckSideDivision )
+void Hexahedron::defineHexahedralFaces( const Solid* solid, const IsInternalFlag intFlag )
 {
+  std::set< TGeomID > concaveFaces; // to avoid connecting nodes laying on them
+
+  if ( solid->HasConcaveVertex() ) {
+    for ( const E_IntersectPoint* ip : _eIntPoints ) {
+      if ( const ConcaveFace* cf = solid->GetConcave( ip->_shapeID )) {
+        if ( this->hasEdgesAround( cf ))
+          concaveFaces.insert( cf->_concaveFace );
+      }
+    }
+    if ( concaveFaces.empty() || concaveFaces.size() * 3  < _eIntPoints.size() ) {
+      for ( const _Node& hexNode: _hexNodes ) {
+        if ( hexNode._node && hexNode._intPoint && hexNode._intPoint->_faceIDs.size() >= 3 )
+          if ( const ConcaveFace* cf = solid->GetConcave( hexNode._node->GetShapeID() ))
+            if ( this->hasEdgesAround( cf ))
+              concaveFaces.insert( cf->_concaveFace );
+      }
+    }
+  }
+
+  const bool toCheckSideDivision = isImplementEdges() || intFlag & IS_CUT_BY_INTERNAL_FACE;
+
+  std::vector< _OrientedLink > splits;
+  std::vector< _Node* >        chainNodes;
   for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
   {
     _Face& quad = _hexQuads[ iF ] ;
@@ -1026,6 +1050,47 @@ void Hexahedron::defineHexahedralFaces( std::vector< _OrientedLink >& splits, st
   }  //
 }
 
+//================================================================================
+/*!
+  * \brief get a remaining link to start from, one lying on minimal nb of FACEs
+  */
+Hexahedron::TFaceOfLink Hexahedron::findStartLink(const vector<_OrientedLink*>& freeLinks,
+                                                  std::set<TGeomID>&            usedFaceIDs)
+{
+  TFaceOfLink faceOfLink( -1, -1 );
+  TFaceOfLink facesOfLink[3] = { faceOfLink, faceOfLink, faceOfLink };
+  for ( size_t iL = 0; iL < freeLinks.size(); ++iL ) {
+    if ( freeLinks[ iL ] ) {
+      std::vector< TGeomID > faces = freeLinks[ iL ]->GetNotUsedFaces( usedFaceIDs );
+      if ( faces.size() == 1 ) {
+        faceOfLink = TFaceOfLink( faces[0], iL );
+        if ( !freeLinks[ iL ]->HasEdgeNodes() )
+          break;
+        facesOfLink[0] = faceOfLink;
+      }
+      else if ( facesOfLink[0].first < 0 ) {
+        faceOfLink = TFaceOfLink(( faces.empty() ? -1 : faces[0]), iL );
+        facesOfLink[ 1 + faces.empty() ] = faceOfLink;
+      }
+    }
+  }
+  for ( int i = 0; faceOfLink.first < 0 && i < 3; ++i ) {
+    faceOfLink = facesOfLink[i];
+  }
+
+  if ( faceOfLink.first < 0 ) { // all faces used
+    for ( size_t iL = 0; iL < freeLinks.size() && faceOfLink.first < 1; ++iL ) {
+      _OrientedLink* curLinkLocal = freeLinks[ iL ];
+      if ( curLinkLocal ) {
+        faceOfLink.first = curLinkLocal->FirstNode()->IsLinked( curLinkLocal->LastNode()->_intPoint );
+        faceOfLink.second = iL;
+      }
+    }
+    usedFaceIDs.clear();
+  }
+  return faceOfLink;
+}
+
 //================================================================================
 /*!
   * \brief Compute mesh volumes resulted from intersection of the Hexahedron
@@ -1041,36 +1106,10 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
   if ( intFlag & IS_CUT_BY_INTERNAL_FACE && !_grid->_toAddEdges ) // Issue #19913
     preventVolumesOverlapping();
 
-  std::set< TGeomID > concaveFaces; // to avoid connecting nodes laying on them
-
-  if ( solid->HasConcaveVertex() )
-  {
-    for ( const E_IntersectPoint* ip : _eIntPoints )
-    {
-      if ( const ConcaveFace* cf = solid->GetConcave( ip->_shapeID ))
-        if ( this->hasEdgesAround( cf ))
-          concaveFaces.insert( cf->_concaveFace );
-    }
-    if ( concaveFaces.empty() || concaveFaces.size() * 3  < _eIntPoints.size() )
-      for ( const _Node& hexNode: _hexNodes )
-      {
-        if ( hexNode._node && hexNode._intPoint && hexNode._intPoint->_faceIDs.size() >= 3 )
-          if ( const ConcaveFace* cf = solid->GetConcave( hexNode._node->GetShapeID() ))
-            if ( this->hasEdgesAround( cf ))
-              concaveFaces.insert( cf->_concaveFace );
-      }
-  }
+  const bool hasEdgeIntersections = !_eIntPoints.empty();
 
   // Create polygons from quadrangles
-  // --------------------------------
-
-  vector< _OrientedLink > splits;
-  vector<_Node*>          chainNodes;
-  _Face*                  coplanarPolyg;
-
-  const bool hasEdgeIntersections = !_eIntPoints.empty();
-  const bool toCheckSideDivision = isImplementEdges() || intFlag & IS_CUT_BY_INTERNAL_FACE;
-  defineHexahedralFaces( splits, chainNodes, concaveFaces, toCheckSideDivision );
+  defineHexahedralFaces( solid, intFlag );
 
   // Create polygons closing holes in a polyhedron
   // ----------------------------------------------
@@ -1080,24 +1119,24 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
     _intNodes[ iN ]._usedInFace = 0;
 
   // add polygons to their links and mark used nodes
-  for ( size_t iP = 0; iP < _polygons.size(); ++iP )
-  {
+  for ( size_t iP = 0; iP < _polygons.size(); ++iP ) {
     _Face& polygon = _polygons[ iP ];
-    for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
-    {
+    for ( size_t iL = 0; iL < polygon._links.size(); ++iL ) {
       polygon._links[ iL ].AddFace( &polygon );
       polygon._links[ iL ].FirstNode()->_usedInFace = &polygon;
     }
   }
+
   // find free links
   vector< _OrientedLink* > freeLinks;
   freeLinks.reserve(20);
   for ( size_t iP = 0; iP < _polygons.size(); ++iP )
   {
     _Face& polygon = _polygons[ iP ];
-    for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
+    for ( size_t iL = 0; iL < polygon._links.size(); ++iL ) {
       if ( polygon._links[ iL ].NbFaces() < 2 )
         freeLinks.push_back( & polygon._links[ iL ]);
+    }
   }
   int nbFreeLinks = freeLinks.size();
   if ( nbFreeLinks == 1 ) return false;
@@ -1123,11 +1162,11 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
   }
 
   std::set<TGeomID> usedFaceIDs;
-  std::vector< TGeomID > faces;
   TGeomID curFace = 0;
   const size_t nbQuadPolygons = _polygons.size();
   E_IntersectPoint ipTmp;
   std::map< TGeomID, std::vector< const B_IntersectPoint* > > tmpAddedFace; // face added to _intPoint
+
   // std::cout << "2\n";
   // create polygons by making closed chains of free links
   size_t iPolygon = _polygons.size();
@@ -1147,9 +1186,11 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
         ( nbFreeLinks < 4 && nbVertexNodes == 0 ))
     {
       // get a remaining link to start from
-      for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
-        if (( curLink = freeLinks[ iL ] ))
+      for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) {
+        curLink = freeLinks[ iL ];
+        if ( curLink )
           freeLinks[ iL ] = 0;
+      }
       polygon._links.push_back( *curLink );
       --nbFreeLinks;
       do
@@ -1170,109 +1211,70 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
     else // there are intersections with EDGEs
     {
       // get a remaining link to start from, one lying on minimal nb of FACEs
-      {
-        typedef pair< TGeomID, int > TFaceOfLink;
-        TFaceOfLink faceOfLink( -1, -1 );
-        TFaceOfLink facesOfLink[3] = { faceOfLink, faceOfLink, faceOfLink };
-        for ( size_t iL = 0; iL < freeLinks.size(); ++iL )
-          if ( freeLinks[ iL ] )
-          {
-            faces = freeLinks[ iL ]->GetNotUsedFace( usedFaceIDs );
-            if ( faces.size() == 1 )
-            {
-              faceOfLink = TFaceOfLink( faces[0], iL );
-              if ( !freeLinks[ iL ]->HasEdgeNodes() )
-                break;
-              facesOfLink[0] = faceOfLink;
-            }
-            else if ( facesOfLink[0].first < 0 )
-            {
-              faceOfLink = TFaceOfLink(( faces.empty() ? -1 : faces[0]), iL );
-              facesOfLink[ 1 + faces.empty() ] = faceOfLink;
-            }
-          }
-        for ( int i = 0; faceOfLink.first < 0 && i < 3; ++i )
-          faceOfLink = facesOfLink[i];
+      TFaceOfLink faceOfLink = findStartLink(freeLinks, usedFaceIDs);
+      curFace = faceOfLink.first;
+      curLink = freeLinks[ faceOfLink.second ];
+      freeLinks[ faceOfLink.second ] = 0;
 
-        if ( faceOfLink.first < 0 ) // all faces used
-        {
-          for ( size_t iL = 0; iL < freeLinks.size() && faceOfLink.first < 1; ++iL )
-            if (( curLink = freeLinks[ iL ]))
-            {
-              faceOfLink.first = 
-                curLink->FirstNode()->IsLinked( curLink->LastNode()->_intPoint );
-              faceOfLink.second = iL;
-            }
-          usedFaceIDs.clear();
-        }
-        curFace = faceOfLink.first;
-        curLink = freeLinks[ faceOfLink.second ];
-        freeLinks[ faceOfLink.second ] = 0;
-      }
       usedFaceIDs.insert( curFace );
       polygon._links.push_back( *curLink );
       --nbFreeLinks;
 
       // find all links lying on a curFace
-      do
-      {
+      do {
         // go forward from curLink
         curNode = curLink->LastNode();
         curLink = 0;
-        for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
+        for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) {
           if ( freeLinks[ iL ] &&
-                freeLinks[ iL ]->FirstNode() == curNode &&
-                freeLinks[ iL ]->LastNode()->IsOnFace( curFace ))
-          {
+               freeLinks[ iL ]->FirstNode() == curNode &&
+               freeLinks[ iL ]->LastNode()->IsOnFace( curFace )) {
             curLink = freeLinks[ iL ];
             freeLinks[ iL ] = 0;
             polygon._links.push_back( *curLink );
             --nbFreeLinks;
           }
+        }
       } while ( curLink );
 
       std::reverse( polygon._links.begin(), polygon._links.end() );
 
       curLink = & polygon._links.back();
-      do
-      {
+      do {
         // go backward from curLink
         curNode = curLink->FirstNode();
         curLink = 0;
-        for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
+        for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) {
           if ( freeLinks[ iL ] &&
-                freeLinks[ iL ]->LastNode() == curNode &&
-                freeLinks[ iL ]->FirstNode()->IsOnFace( curFace ))
-          {
+               freeLinks[ iL ]->LastNode() == curNode &&
+               freeLinks[ iL ]->FirstNode()->IsOnFace( curFace )) {
             curLink = freeLinks[ iL ];
             freeLinks[ iL ] = 0;
             polygon._links.push_back( *curLink );
             --nbFreeLinks;
           }
+        }
       } while ( curLink );
 
       curNode = polygon._links.back().FirstNode();
 
-      if ( polygon._links[0].LastNode() != curNode )
-      {
-        if ( nbVertexNodes > 0 )
-        {
+      if ( polygon._links[0].LastNode() != curNode ) {
+        if ( nbVertexNodes > 0 ) {
           // add links with _vIntNodes if not already used
-          chainNodes.clear();
-          for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN )
+          vector< _Node* > chainNodes;
+          //chainNodes.clear();
+          for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN ) {
             if ( !_vIntNodes[ iN ]->IsUsedInFace() &&
-                  _vIntNodes[ iN ]->IsOnFace( curFace ))
-            {
+                 _vIntNodes[ iN ]->IsOnFace( curFace )) {
               _vIntNodes[ iN ]->_usedInFace = &polygon;
               chainNodes.push_back( _vIntNodes[ iN ] );
             }
+          }
           if ( chainNodes.size() > 1 &&
-                curFace != _grid->PseudoIntExtFaceID() ) /////// TODO
-          {
+               curFace != _grid->PseudoIntExtFaceID() ) { /////// TODO
             sortVertexNodes( chainNodes, curNode, curFace );
           }
-          for ( size_t i = 0; i < chainNodes.size(); ++i )
-          {
+          for ( size_t i = 0; i < chainNodes.size(); ++i ) {
             polygon.AddPolyLink( chainNodes[ i ], curNode );
             curNode = chainNodes[ i ];
             freeLinks.push_back( &polygon._links.back() );
@@ -1290,7 +1292,7 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
     } // if there are intersections with EDGEs
 
     if ( polygon._links.size() < 2 ||
-          polygon._links[0].LastNode() != polygon._links.back().FirstNode() )
+         polygon._links[0].LastNode() != polygon._links.back().FirstNode() )
     {
       _polygons.clear();
       break; // closed polygon not found -> invalid polyhedron
@@ -1311,7 +1313,7 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
       if ( iPolygon == _polygons.size()-1 )
         _polygons.pop_back();
     }
-    else // polygon._links.size() >= 2
+    else // polygon._links.size() > 2
     {
       // add polygon to its links
       for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
@@ -1322,7 +1324,7 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
       if ( /*hasEdgeIntersections &&*/ iPolygon == _polygons.size() - 1 )
       {
         // check that a polygon does not lie on a hexa side
-        coplanarPolyg = 0;
+        _Face* coplanarPolyg = 0;
         for ( size_t iL = 0; iL < polygon._links.size() && !coplanarPolyg; ++iL )
         {
           if ( polygon._links[ iL ].NbFaces() < 2 )
@@ -1332,11 +1334,13 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
           size_t iL2;
           coplanarPolyg = polygon._links[ iL ]._link->_faces[0];
           for ( iL2 = iL + 1; iL2 < polygon._links.size(); ++iL2 )
+          {
             if ( polygon._links[ iL2 ]._link->_faces[0] == coplanarPolyg &&
                   !coplanarPolyg->IsPolyLink( polygon._links[ iL  ]) &&
                   !coplanarPolyg->IsPolyLink( polygon._links[ iL2 ]) &&
                   coplanarPolyg < & _polygons[ nbQuadPolygons ])
               break;
+          }
           if ( iL2 == polygon._links.size() )
             coplanarPolyg = 0;
         }
@@ -1350,11 +1354,13 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
 
           // fill freeLinks with links not shared by coplanarPolyg and polygon
           for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
+          {
             if ( polygon._links[ iL ]._link->_faces[1] &&
                   polygon._links[ iL ]._link->_faces[0] != coplanarPolyg )
             {
               _Face* p = polygon._links[ iL ]._link->_faces[0];
               for ( size_t iL2 = 0; iL2 < p->_links.size(); ++iL2 )
+              {
                 if ( p->_links[ iL2 ]._link == polygon._links[ iL ]._link )
                 {
                   freeLinks.push_back( & p->_links[ iL2 ] );
@@ -1362,8 +1368,11 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
                   freeLinks.back()->RemoveFace( &polygon );
                   break;
                 }
+              }
             }
+          }
           for ( size_t iL = 0; iL < coplanarPolyg->_links.size(); ++iL )
+          {
             if ( coplanarPolyg->_links[ iL ]._link->_faces[1] &&
                   coplanarPolyg->_links[ iL ]._link->_faces[1] != &polygon )
             {
@@ -1371,6 +1380,7 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
               if ( p == coplanarPolyg )
                 p = coplanarPolyg->_links[ iL ]._link->_faces[1];
               for ( size_t iL2 = 0; iL2 < p->_links.size(); ++iL2 )
+              {
                 if ( p->_links[ iL2 ]._link == coplanarPolyg->_links[ iL ]._link )
                 {
                   // set links of coplanarPolyg in place of used freeLinks
@@ -1395,9 +1405,12 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
                   }
                   break;
                 }
+              }
             }
+          }
           // set coplanarPolyg to be re-created next
           for ( size_t iP = 0; iP < _polygons.size(); ++iP )
+          {
             if ( coplanarPolyg == & _polygons[ iP ] )
             {
               iPolygon = iP;
@@ -1405,6 +1418,7 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
               _polygons[ iPolygon ]._polyLinks.clear();
               break;
             }
+          }
           _polygons.pop_back();
           usedFaceIDs.erase( curFace );
           continue;
@@ -1415,6 +1429,7 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
 
     } // end of case ( polygon._links.size() > 2 )
   } // while ( nbFreeLinks > 0 )
+
   // std::cout << "3\n";
   for ( auto & face_ip : tmpAddedFace )
   {
@@ -1435,13 +1450,14 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
   _hasTooSmall = ! checkPolyhedronSize( intFlag & IS_CUT_BY_INTERNAL_FACE, volSize );
 
   for ( size_t i = 0; i < 8; ++i )
+  {
     if ( _hexNodes[ i ]._intPoint == &ipTmp )
       _hexNodes[ i ]._intPoint = 0;
+  }
 
   if ( _hasTooSmall )
     return false; // too small volume
 
-
   // Try to find out names of no-name polygons (issue # 19887)
   if ( _grid->IsToRemoveExcessEntities() && _polygons.back()._name == SMESH_Block::ID_NONE )
   {
@@ -1474,8 +1490,7 @@ bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
     }
   }
 
-
-  /* This call is irrelevant here because _volumeDefs datas are were not filled! 
+  /* This call is irrelevant here because _volumeDefs datas were not filled! 
   or .... it is potentialy filled by other thread?? */
   _volumeDefs._nodes.clear();
   _volumeDefs._quantities.clear();
@@ -4223,4 +4238,4 @@ bool Hexahedron::_SplitIterator::Next()
     }
   }
   return More();
-}
\ No newline at end of file
+}
index c87877b4d79c5eb1c1a787b01f1c0eaceaccce51..ec1e1d0fe44a722f1fe7d9d6f037c1519844eb94 100644 (file)
@@ -35,9 +35,9 @@
 #include "SMESH_StdMeshers.hxx"
 #include "StdMeshers_Cartesian_3D_Grid.hxx"
 
-using namespace std;
-using namespace SMESH;
-namespace
+namespace StdMeshers
+{
+namespace Cartesian3D
 {
   // --------------------------------------------------------------------------
   /*!
@@ -49,14 +49,14 @@ namespace
     int    _dInd[4][3];
     size_t _nbCells[3];
     int    _i,_j,_k;
-    Grid*  _grid;
+    StdMeshers::Cartesian3D::Grid*  _grid;
 
-    CellsAroundLink( Grid* grid, int iDir ):
+    CellsAroundLink( StdMeshers::Cartesian3D::Grid* grid, int iDir ):
       _iDir( iDir ),
       _dInd{ {0,0,0}, {0,0,0}, {0,0,0}, {0,0,0} },
       _nbCells{ grid->_coords[0].size() - 1,
-          grid->_coords[1].size() - 1,
-          grid->_coords[2].size() - 1 },
+      grid->_coords[1].size() - 1,
+      grid->_coords[2].size() - 1 },
       _grid( grid )
     {
       const int iDirOther[3][2] = {{ 1,2 },{ 0,2 },{ 0,1 }};
@@ -85,419 +85,433 @@ namespace
       return true;
     }
   };
-}
-  
-// --------------------------------------------------------------------------
-/*!
-  * \brief Class representing topology of the hexahedron and creating a mesh
-  *        volume basing on analysis of hexahedron intersection with geometry
-  */
-class STDMESHERS_EXPORT Hexahedron
-{
-  // --------------------------------------------------------------------------------
-  struct _Face;
-  struct _Link;
-  enum IsInternalFlag { IS_NOT_INTERNAL, IS_INTERNAL, IS_CUT_BY_INTERNAL_FACE };
-  // --------------------------------------------------------------------------------
-  struct _Node //!< node either at a hexahedron corner or at intersection
-  {
-    const SMDS_MeshNode*    _node;        // mesh node at hexahedron corner
-    const SMDS_MeshNode*    _boundaryCornerNode; // missing mesh node due to hex truncation on the boundary
-    const B_IntersectPoint* _intPoint;
-    const _Face*            _usedInFace;
-    char                    _isInternalFlags;
-
-    _Node(const SMDS_MeshNode* n=0, const B_IntersectPoint* ip=0)
-      :_node(n), _intPoint(ip), _usedInFace(0), _isInternalFlags(0) {} 
-    const SMDS_MeshNode*    Node() const
-    { return ( _intPoint && _intPoint->_node ) ? _intPoint->_node : _node; }
-    const SMDS_MeshNode*    BoundaryNode() const
-    { return _node ? _node : _boundaryCornerNode; }
-    const E_IntersectPoint* EdgeIntPnt() const
-    { return static_cast< const E_IntersectPoint* >( _intPoint ); }
-    const F_IntersectPoint* FaceIntPnt() const
-    { return static_cast< const F_IntersectPoint* >( _intPoint ); }
-    const vector< TGeomID >& faces() const { return _intPoint->_faceIDs; }
-    TGeomID face(size_t i) const { return _intPoint->_faceIDs[ i ]; }
-    void SetInternal( IsInternalFlag intFlag ) { _isInternalFlags |= intFlag; }
-    bool IsCutByInternal() const { return _isInternalFlags & IS_CUT_BY_INTERNAL_FACE; }
-    bool IsUsedInFace( const _Face* polygon = 0 )
-    {
-      return polygon ? ( _usedInFace == polygon ) : bool( _usedInFace );
-    }
-    TGeomID IsLinked( const B_IntersectPoint* other,
-                      TGeomID                 avoidFace=-1 ) const // returns id of a common face
-    {
-      return _intPoint ? _intPoint->HasCommonFace( other, avoidFace ) : 0;
-    }
-    bool IsOnFace( TGeomID faceID ) const // returns true if faceID is found
-    {
-      return _intPoint ? _intPoint->IsOnFace( faceID ) : false;
-    }
-    size_t GetCommonFaces( const B_IntersectPoint * other, TGeomID* common ) const
-    {
-      return _intPoint && other ? _intPoint->GetCommonFaces( other, common ) : 0;
-    }
-    gp_Pnt Point() const
-    {
-      if ( const SMDS_MeshNode* n = Node() )
-        return SMESH_NodeXYZ( n );
-      if ( const E_IntersectPoint* eip =
-            dynamic_cast< const E_IntersectPoint* >( _intPoint ))
-        return eip->_point;
-      return gp_Pnt( 1e100, 0, 0 );
-    }
-    TGeomID ShapeID() const
-    {
-      if ( const E_IntersectPoint* eip = dynamic_cast< const E_IntersectPoint* >( _intPoint ))
-        return eip->_shapeID;
-      return 0;
-    }
 
-    void Add( const E_IntersectPoint* ip );
-  };
-  // --------------------------------------------------------------------------------
-  struct _Link // link connecting two _Node's
-  {
-    _Node* _nodes[2];
-    _Face* _faces[2]; // polygons sharing a link
-    vector< const F_IntersectPoint* > _fIntPoints; // GridLine intersections with FACEs
-    vector< _Node* >                  _fIntNodes;   // _Node's at _fIntPoints
-    vector< _Link >                   _splits;
-    _Link(): _faces{ 0, 0 } {}
-  };
-  // --------------------------------------------------------------------------------
-  struct _OrientedLink
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Class representing topology of the hexahedron and creating a mesh
+   *        volume basing on analysis of hexahedron intersection with geometry
+   */
+  class STDMESHERS_EXPORT Hexahedron
   {
-    _Link* _link;
-    bool   _reverse;
-    _OrientedLink( _Link* link=0, bool reverse=false ): _link(link), _reverse(reverse) {}
-    void Reverse() { _reverse = !_reverse; }
-    size_t NbResultLinks() const { return _link->_splits.size(); }
-    _OrientedLink ResultLink(int i) const
+    // --------------------------------------------------------------------------------
+    struct _Face;
+    struct _Link;
+    enum IsInternalFlag { IS_NOT_INTERNAL, IS_INTERNAL, IS_CUT_BY_INTERNAL_FACE };
+    // --------------------------------------------------------------------------------
+    struct _Node //!< node either at a hexahedron corner or at intersection
     {
-      return _OrientedLink(&_link->_splits[_reverse ? NbResultLinks()-i-1 : i],_reverse);
-    }
-    _Node* FirstNode() const { return _link->_nodes[ _reverse ]; }
-    _Node* LastNode()  const { return _link->_nodes[ !_reverse ]; }
-    operator bool() const { return _link; }
-    vector< TGeomID > GetNotUsedFace(const set<TGeomID>& usedIDs ) const // returns supporting FACEs
-    {
-      vector< TGeomID > faces;
-      const B_IntersectPoint *ip0, *ip1;
-      if (( ip0 = _link->_nodes[0]->_intPoint ) &&
-          ( ip1 = _link->_nodes[1]->_intPoint ))
+      const SMDS_MeshNode*    _node;        // mesh node at hexahedron corner
+      const SMDS_MeshNode*    _boundaryCornerNode; // missing mesh node due to hex truncation on the boundary
+      const StdMeshers::Cartesian3D::B_IntersectPoint* _intPoint;
+      const _Face*            _usedInFace;
+      char                    _isInternalFlags;
+
+      _Node(const SMDS_MeshNode* n=0, const StdMeshers::Cartesian3D::B_IntersectPoint* ip=0)
+        :_node(n), _intPoint(ip), _usedInFace(0), _isInternalFlags(0) {} 
+      const SMDS_MeshNode*    Node() const
+      { return ( _intPoint && _intPoint->_node ) ? _intPoint->_node : _node; }
+      const SMDS_MeshNode*    BoundaryNode() const
+      { return _node ? _node : _boundaryCornerNode; }
+      const StdMeshers::Cartesian3D::E_IntersectPoint* EdgeIntPnt() const
+      { return static_cast< const StdMeshers::Cartesian3D::E_IntersectPoint* >( _intPoint ); }
+      const StdMeshers::Cartesian3D::F_IntersectPoint* FaceIntPnt() const
+      { return static_cast< const StdMeshers::Cartesian3D::F_IntersectPoint* >( _intPoint ); }
+      const std::vector< StdMeshers::Cartesian3D::TGeomID >& faces() const { return _intPoint->_faceIDs; }
+      StdMeshers::Cartesian3D::TGeomID face(size_t i) const { return _intPoint->_faceIDs[ i ]; }
+      void SetInternal( IsInternalFlag intFlag ) { _isInternalFlags |= intFlag; }
+      bool IsCutByInternal() const { return _isInternalFlags & IS_CUT_BY_INTERNAL_FACE; }
+      bool IsUsedInFace( const _Face* polygon = 0 )
       {
-        for ( size_t i = 0; i < ip0->_faceIDs.size(); ++i )
-          if ( ip1->IsOnFace ( ip0->_faceIDs[i] ) &&
-                !usedIDs.count( ip0->_faceIDs[i] ) )
-            faces.push_back( ip0->_faceIDs[i] );
+        return polygon ? ( _usedInFace == polygon ) : bool( _usedInFace );
       }
-      return faces;
-    }
-    bool HasEdgeNodes() const
-    {
-      return ( dynamic_cast< const E_IntersectPoint* >( _link->_nodes[0]->_intPoint ) ||
-                dynamic_cast< const E_IntersectPoint* >( _link->_nodes[1]->_intPoint ));
-    }
-    int NbFaces() const
-    {
-      return !_link->_faces[0] ? 0 : 1 + bool( _link->_faces[1] );
-    }
-    void AddFace( _Face* f )
-    {
-      if ( _link->_faces[0] )
+      StdMeshers::Cartesian3D::TGeomID IsLinked( const StdMeshers::Cartesian3D::B_IntersectPoint* other,
+                                                 StdMeshers::Cartesian3D::TGeomID avoidFace=-1 ) const // returns id of a common face
       {
-        _link->_faces[1] = f;
+        return _intPoint ? _intPoint->HasCommonFace( other, avoidFace ) : 0;
       }
-      else
+      bool IsOnFace( StdMeshers::Cartesian3D::TGeomID faceID ) const // returns true if faceID is found
       {
-        _link->_faces[0] = f;
-        _link->_faces[1] = 0;
+        return _intPoint ? _intPoint->IsOnFace( faceID ) : false;
       }
-    }
-    void RemoveFace( _Face* f )
+      size_t GetCommonFaces( const StdMeshers::Cartesian3D::B_IntersectPoint * other,
+                             StdMeshers::Cartesian3D::TGeomID* common ) const
+      {
+        return _intPoint && other ? _intPoint->GetCommonFaces( other, common ) : 0;
+      }
+      gp_Pnt Point() const
+      {
+        if ( const SMDS_MeshNode* n = Node() )
+          return SMESH_NodeXYZ( n );
+        if ( const StdMeshers::Cartesian3D::E_IntersectPoint* eip =
+             dynamic_cast< const StdMeshers::Cartesian3D::E_IntersectPoint* >( _intPoint ))
+          return eip->_point;
+        return gp_Pnt( 1e100, 0, 0 );
+      }
+      StdMeshers::Cartesian3D::TGeomID ShapeID() const
+      {
+        if ( const StdMeshers::Cartesian3D::E_IntersectPoint* eip =
+             dynamic_cast< const StdMeshers::Cartesian3D::E_IntersectPoint* >( _intPoint ))
+          return eip->_shapeID;
+        return 0;
+      }
+
+      void Add( const StdMeshers::Cartesian3D::E_IntersectPoint* ip );
+    };
+
+    // --------------------------------------------------------------------------------
+    struct _Link // link connecting two _Node's
     {
-      if ( !_link->_faces[0] ) return;
+      _Node* _nodes[2];
+      _Face* _faces[2]; // polygons sharing a link
+      std::vector< const StdMeshers::Cartesian3D::F_IntersectPoint* > _fIntPoints; // GridLine intersections with FACEs
+      std::vector< _Node* > _fIntNodes;   // _Node's at _fIntPoints
+      std::vector< _Link >  _splits;
+      _Link(): _faces{ 0, 0 } {}
+    };
 
-      if ( _link->_faces[1] == f )
+    // --------------------------------------------------------------------------------
+    struct _OrientedLink
+    {
+      _Link* _link;
+      bool   _reverse;
+      _OrientedLink( _Link* link=0, bool reverse=false ): _link(link), _reverse(reverse) {}
+      void Reverse() { _reverse = !_reverse; }
+      size_t NbResultLinks() const { return _link->_splits.size(); }
+      _OrientedLink ResultLink(int i) const
       {
-        _link->_faces[1] = 0;
+        return _OrientedLink(&_link->_splits[_reverse ? NbResultLinks()-i-1 : i],_reverse);
       }
-      else if ( _link->_faces[0] == f )
+      _Node* FirstNode() const { return _link->_nodes[ _reverse ]; }
+      _Node* LastNode()  const { return _link->_nodes[ !_reverse ]; }
+      operator bool() const { return _link; }
+
+      // returns supporting FACEs
+      std::vector< StdMeshers::Cartesian3D::TGeomID > GetNotUsedFaces
+      (const std::set<StdMeshers::Cartesian3D::TGeomID>& usedIDs ) const
       {
-        _link->_faces[0] = 0;
-        if ( _link->_faces[1] )
+        std::vector< StdMeshers::Cartesian3D::TGeomID > faces;
+        const StdMeshers::Cartesian3D::B_IntersectPoint *ip0, *ip1;
+        if (( ip0 = _link->_nodes[0]->_intPoint ) &&
+            ( ip1 = _link->_nodes[1]->_intPoint ))
         {
-          _link->_faces[0] = _link->_faces[1];
+          for ( size_t i = 0; i < ip0->_faceIDs.size(); ++i )
+            if ( ip1->IsOnFace ( ip0->_faceIDs[i] ) &&
+                 !usedIDs.count( ip0->_faceIDs[i] ) )
+              faces.push_back( ip0->_faceIDs[i] );
+        }
+        return faces;
+      }
+
+      bool HasEdgeNodes() const
+      {
+        return ( dynamic_cast< const StdMeshers::Cartesian3D::E_IntersectPoint* >( _link->_nodes[0]->_intPoint ) ||
+                 dynamic_cast< const StdMeshers::Cartesian3D::E_IntersectPoint* >( _link->_nodes[1]->_intPoint ));
+      }
+      int NbFaces() const
+      {
+        return !_link->_faces[0] ? 0 : 1 + bool( _link->_faces[1] );
+      }
+      void AddFace( _Face* f )
+      {
+        if ( _link->_faces[0] ) {
+          _link->_faces[1] = f;
+        }
+        else {
+          _link->_faces[0] = f;
           _link->_faces[1] = 0;
         }
       }
-    }
-  };
-  // --------------------------------------------------------------------------------
-  struct _SplitIterator //! set to _hexLinks splits on one side of INTERNAL FACEs
-  {
-    struct _Split // data of a link split
-    {
-      int    _linkID;          // hex link ID
-      _Node* _nodes[2];
-      int    _iCheckIteration; // iteration where split is tried as Hexahedron split
-      _Link* _checkedSplit;    // split set to hex links
-      bool   _isUsed;          // used in a volume
+      void RemoveFace( const _Face* f )
+      {
+        if ( !_link->_faces[0] ) return;
 
-      _Split( _Link & split, int iLink ):
-        _linkID( iLink ), _nodes{ split._nodes[0], split._nodes[1] },
-        _iCheckIteration( 0 ), _isUsed( false )
-      {}
-      bool IsCheckedOrUsed( bool used ) const { return used ? _isUsed : _iCheckIteration > 0; }
-    };
-    _Link*                _hexLinks;
-    std::vector< _Split > _splits;
-    int                   _iterationNb;
-    size_t                _nbChecked;
-    size_t                _nbUsed;
-    std::vector< _Node* > _freeNodes; // nodes reached while composing a split set
-
-    _SplitIterator( _Link* hexLinks ):
-      _hexLinks( hexLinks ), _iterationNb(0), _nbChecked(0), _nbUsed(0)
-    {
-      _freeNodes.reserve( 12 );
-      _splits.reserve( 24 );
-      for ( int iL = 0; iL < 12; ++iL )
-        for ( size_t iS = 0; iS < _hexLinks[ iL ]._splits.size(); ++iS )
-          _splits.emplace_back( _hexLinks[ iL ]._splits[ iS ], iL );
-      Next();
-    }
-    bool More() const { return _nbUsed < _splits.size(); }
-    bool Next();
-  };
-  // --------------------------------------------------------------------------------
-  struct _Face
-  {
-    SMESH_Block::TShapeID   _name;
-    vector< _OrientedLink > _links;       // links on GridLine's
-    vector< _Link >         _polyLinks;   // links added to close a polygonal face
-    vector< _Node* >        _eIntNodes;   // nodes at intersection with EDGEs
-
-    _Face():_name( SMESH_Block::ID_NONE )
-    {}
-    bool IsPolyLink( const _OrientedLink& ol )
-    {
-      return _polyLinks.empty() ? false :
-        ( &_polyLinks[0] <= ol._link &&  ol._link <= &_polyLinks.back() );
-    }
-    void AddPolyLink(_Node* n0, _Node* n1, _Face* faceToFindEqual=0)
-    {
-      if ( faceToFindEqual && faceToFindEqual != this ) {
-        for ( size_t iL = 0; iL < faceToFindEqual->_polyLinks.size(); ++iL )
-          if ( faceToFindEqual->_polyLinks[iL]._nodes[0] == n1 &&
-                faceToFindEqual->_polyLinks[iL]._nodes[1] == n0 )
-          {
-            _links.push_back
-              ( _OrientedLink( & faceToFindEqual->_polyLinks[iL], /*reverse=*/true ));
-            return;
+        if ( _link->_faces[1] == f ) {
+          _link->_faces[1] = 0;
+        }
+        else if ( _link->_faces[0] == f ) {
+          _link->_faces[0] = 0;
+          if ( _link->_faces[1] ) {
+            _link->_faces[0] = _link->_faces[1];
+            _link->_faces[1] = 0;
           }
+        }
       }
-      _Link l;
-      l._nodes[0] = n0;
-      l._nodes[1] = n1;
-      _polyLinks.push_back( l );
-      _links.push_back( _OrientedLink( &_polyLinks.back() ));
-    }
-  };
-  // --------------------------------------------------------------------------------
-  struct _volumeDef // holder of nodes of a volume mesh element
-  {
-    typedef void* _ptr;
+    };
 
-    struct _nodeDef
+    // --------------------------------------------------------------------------------
+    struct _SplitIterator //! set to _hexLinks splits on one side of INTERNAL FACEs
     {
-      const SMDS_MeshNode*    _node; // mesh node at hexahedron corner
-      const B_IntersectPoint* _intPoint;
-
-      _nodeDef(): _node(0), _intPoint(0) {}
-      _nodeDef( _Node* n ): _node( n->_node), _intPoint( n->_intPoint ) {}
-      const SMDS_MeshNode*    Node() const
-      { return ( _intPoint && _intPoint->_node ) ? _intPoint->_node : _node; }
-      const E_IntersectPoint* EdgeIntPnt() const
-      { return static_cast< const E_IntersectPoint* >( _intPoint ); }
-      _ptr Ptr() const { return Node() ? (_ptr) Node() : (_ptr) EdgeIntPnt(); }
-      bool operator==(const _nodeDef& other ) const { return Ptr() == other.Ptr(); }
+      struct _Split // data of a link split
+      {
+        int    _linkID;          // hex link ID
+        _Node* _nodes[2];
+        int    _iCheckIteration; // iteration where split is tried as Hexahedron split
+        _Link* _checkedSplit;    // split set to hex links
+        bool   _isUsed;          // used in a volume
+
+        _Split( _Link & split, int iLink ):
+          _linkID( iLink ), _nodes{ split._nodes[0], split._nodes[1] },
+          _iCheckIteration( 0 ), _isUsed( false )
+        {}
+        bool IsCheckedOrUsed( bool used ) const { return used ? _isUsed : _iCheckIteration > 0; }
+      };
+      _Link*                _hexLinks;
+      std::vector< _Split > _splits;
+      int                   _iterationNb;
+      size_t                _nbChecked;
+      size_t                _nbUsed;
+      std::vector< _Node* > _freeNodes; // nodes reached while composing a split set
+
+      _SplitIterator( _Link* hexLinks ):
+        _hexLinks( hexLinks ), _iterationNb(0), _nbChecked(0), _nbUsed(0)
+      {
+        _freeNodes.reserve( 12 );
+        _splits.reserve( 24 );
+        for ( int iL = 0; iL < 12; ++iL )
+          for ( size_t iS = 0; iS < _hexLinks[ iL ]._splits.size(); ++iS )
+            _splits.emplace_back( _hexLinks[ iL ]._splits[ iS ], iL );
+        Next();
+      }
+      bool More() const { return _nbUsed < _splits.size(); }
+      bool Next();
     };
 
-    vector< _nodeDef >      _nodes;
-    vector< int >           _quantities;
-    _volumeDef*             _next; // to store several _volumeDefs in a chain
-    TGeomID                 _solidID;
-    double                  _size;
-    const SMDS_MeshElement* _volume; // new volume
-    std::vector<const SMDS_MeshElement*> _brotherVolume; // produced due to poly split 
-
-    vector< SMESH_Block::TShapeID > _names; // name of side a polygon originates from
-
-    _volumeDef(): _next(0), _solidID(0), _size(0), _volume(0) {}
-    ~_volumeDef() { delete _next; }
-    _volumeDef( _volumeDef& other ):
-      _next(0), _solidID( other._solidID ), _size( other._size ), _volume( other._volume )
-    { _nodes.swap( other._nodes ); _quantities.swap( other._quantities ); other._volume = 0;
-      _names.swap( other._names ); }
-
-    size_t size() const { return 1 + ( _next ? _next->size() : 0 ); } // nb _volumeDef in a chain
-    _volumeDef* at(int index)
-    { return index == 0 ? this : ( _next ? _next->at(index-1) : _next ); }
-
-    void Set( _Node** nodes, int nb )
-    { _nodes.assign( nodes, nodes + nb ); }
-
-    void SetNext( _volumeDef* vd )
-    { if ( _next ) { _next->SetNext( vd ); } else { _next = vd; }}
-
-    bool IsEmpty() const { return (( _nodes.empty() ) &&
-                                    ( !_next || _next->IsEmpty() )); }
-    bool IsPolyhedron() const { return ( !_quantities.empty() ||
-                                          ( _next && !_next->_quantities.empty() )); }
-
-
-    struct _linkDef: public std::pair<_ptr,_ptr> // to join polygons in removeExcessSideDivision()
+    // --------------------------------------------------------------------------------
+    struct _Face
     {
-      _nodeDef _node1;//, _node2;
-      mutable /*const */_linkDef *_prev, *_next;
-      size_t _loopIndex;
-
-      _linkDef():_prev(0), _next(0) {}
+      SMESH_Block::TShapeID   _name;
+      std::vector< _OrientedLink > _links;       // links on GridLine's
+      std::vector< _Link >         _polyLinks;   // links added to close a polygonal face
+      std::vector< _Node* >        _eIntNodes;   // nodes at intersection with EDGEs
 
-      void init( const _nodeDef& n1, const _nodeDef& n2, size_t iLoop )
+      _Face():_name( SMESH_Block::ID_NONE )
+      {}
+      bool IsPolyLink( const _OrientedLink& ol )
       {
-        _node1     = n1; //_node2 = n2;
-        _loopIndex = iLoop;
-        first      = n1.Ptr();
-        second     = n2.Ptr();
-        if ( first > second ) std::swap( first, second );
+        return _polyLinks.empty() ? false :
+          ( &_polyLinks[0] <= ol._link &&  ol._link <= &_polyLinks.back() );
       }
-      void setNext( _linkDef* next )
+      void AddPolyLink(_Node* n0, _Node* n1, _Face* faceToFindEqual=0)
       {
-        _next = next;
-        next->_prev = this;
+        if ( faceToFindEqual && faceToFindEqual != this ) {
+          for ( size_t iL = 0; iL < faceToFindEqual->_polyLinks.size(); ++iL )
+            if ( faceToFindEqual->_polyLinks[iL]._nodes[0] == n1 &&
+                 faceToFindEqual->_polyLinks[iL]._nodes[1] == n0 )
+            {
+              _links.push_back
+                ( _OrientedLink( & faceToFindEqual->_polyLinks[iL], /*reverse=*/true ));
+              return;
+            }
+        }
+        _Link l;
+        l._nodes[0] = n0;
+        l._nodes[1] = n1;
+        _polyLinks.push_back( l );
+        _links.push_back( _OrientedLink( &_polyLinks.back() ));
       }
     };
-  };
 
-  // topology of a hexahedron
-  _Node _hexNodes [8];
-  _Link _hexLinks [12];
-  _Face _hexQuads [6];
-
-  // faces resulted from hexahedron intersection
-  vector< _Face > _polygons;
-
-  // intresections with EDGEs
-  vector< const E_IntersectPoint* > _eIntPoints;
-
-  // additional nodes created at intersection points
-  vector< _Node > _intNodes;
-
-  // nodes inside the hexahedron (at VERTEXes) refer to _intNodes
-  vector< _Node* > _vIntNodes;
-
-  // computed volume elements
-  _volumeDef _volumeDefs;
-
-  Grid*       _grid;
-  double      _sideLength[3];
-  int         _nbCornerNodes, _nbFaceIntNodes, _nbBndNodes;
-  int         _origNodeInd; // index of _hexNodes[0] node within the _grid
-  size_t      _i,_j,_k;
-  bool        _hasTooSmall;
-  int         _cellID;
-
-public:
-  Hexahedron(Grid* grid);
-  int MakeElements(SMESH_MesherHelper&                      helper,
-                    const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap, const int numOfThreads = 1 );
-  void computeElements( const Solid* solid = 0, int solidIndex = -1 );
-
-private:
-  Hexahedron(const Hexahedron& other, size_t i, size_t j, size_t k, int cellID );
-  void init( size_t i, size_t j, size_t k, const Solid* solid=0 );
-  void init( size_t i );
-  void setIJK( size_t i );
-  /*Auxiliary methods to extract operations from monolitic compute method*/
-  void defineHexahedralFaces( std::vector< _OrientedLink >& splits, std::vector<_Node*>& chainNodes, std::set< TGeomID >& concaveFaces, bool toCheckSideDivision );
-  bool compute( const Solid* solid, const IsInternalFlag intFlag );
-  size_t getSolids( TGeomID ids[] );
-  bool isCutByInternalFace( IsInternalFlag & maxFlag );
-  void addEdges(SMESH_MesherHelper&                      helper,
-                vector< Hexahedron* >&                   intersectedHex,
-                const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap);
-  gp_Pnt findIntPoint( double u1, double proj1, double u2, double proj2,
-                        double proj, BRepAdaptor_Curve& curve,
-                        const gp_XYZ& axis, const gp_XYZ& origin );
-  int  getEntity( const E_IntersectPoint* ip, int* facets, int& sub );
-  bool addIntersection( const E_IntersectPoint* ip,
-                        vector< Hexahedron* >&  hexes,
-                        int ijk[], int dIJK[] );
-  bool isQuadOnFace( const size_t iQuad );
-  bool findChain( _Node* n1, _Node* n2, _Face& quad, vector<_Node*>& chainNodes );
-  bool closePolygon( _Face* polygon, vector<_Node*>& chainNodes ) const;
-  bool findChainOnEdge( const vector< _OrientedLink >& splits,
-                        const _OrientedLink&           prevSplit,
-                        const _OrientedLink&           avoidSplit,
-                        const std::set< TGeomID > &    concaveFaces,
-                        size_t &                       iS,
-                        _Face&                         quad,
-                        vector<_Node*>&                chn);
-  int  addVolumes(SMESH_MesherHelper& helper );
-  void addFaces( SMESH_MesherHelper&                       helper,
-                  const vector< const SMDS_MeshElement* > & boundaryVolumes );
-  void addSegments( SMESH_MesherHelper&                      helper,
-                    const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap );
-  void getVolumes( vector< const SMDS_MeshElement* > & volumes );
-  void getBoundaryElems( vector< const SMDS_MeshElement* > & boundaryVolumes );
-  void removeExcessSideDivision(const vector< Hexahedron* >& allHexa);
-  void removeExcessNodes(vector< Hexahedron* >& allHexa);
-  void preventVolumesOverlapping();
-  TGeomID getAnyFace() const;
-  void cutByExtendedInternal( std::vector< Hexahedron* >& hexes,
-                              const TColStd_MapOfInteger& intEdgeIDs );
-  gp_Pnt mostDistantInternalPnt( int hexIndex, const gp_Pnt& p1, const gp_Pnt& p2 );
-  bool isOutPoint( _Link& link, int iP, SMESH_MesherHelper& helper, const Solid* solid ) const;
-  void sortVertexNodes(vector<_Node*>& nodes, _Node* curNode, TGeomID face);
-  bool isInHole() const;
-  bool hasStrangeEdge() const;
-  bool checkPolyhedronSize( bool isCutByInternalFace, double & volSize ) const;
-  int checkPolyhedronValidity( _volumeDef* volDef, std::vector<std::vector<int>>& splitQuantities, 
-                                std::vector<std::vector<const SMDS_MeshNode*>>& splitNodes );
-  const SMDS_MeshElement* addPolyhedronToMesh( _volumeDef* volDef,  SMESH_MesherHelper& helper, const std::vector<const SMDS_MeshNode*>& nodes, 
-                                              const std::vector<int>& quantities );
-  bool addHexa ();
-  bool addTetra();
-  bool addPenta();
-  bool addPyra ();
-  bool debugDumpLink( _Link* link );
-  _Node* findEqualNode( vector< _Node* >&       nodes,
-                        const E_IntersectPoint* ip,
-                        const double            tol2 )
-  {
-    for ( size_t i = 0; i < nodes.size(); ++i )
-      if ( nodes[i]->EdgeIntPnt() == ip ||
-            nodes[i]->Point().SquareDistance( ip->_point ) <= tol2 )
-        return nodes[i];
-    return 0;
-  }
-  bool isCorner( const _Node* node ) const { return ( node >= &_hexNodes[0] &&
-                                                      node -  &_hexNodes[0] < 8 ); }
-  bool hasEdgesAround( const ConcaveFace* cf ) const;
-  bool isImplementEdges() const { return _grid->_edgeIntPool.nbElements(); }
-  bool isOutParam(const double uvw[3]) const;
-
-  typedef boost::container::flat_map< TGeomID, size_t > TID2Nb;
-  static void insertAndIncrement( TGeomID id, TID2Nb& id2nbMap )
-  {
-    TID2Nb::value_type s0( id, 0 );
-    TID2Nb::iterator id2nb = id2nbMap.insert( s0 ).first;
-    id2nb->second++;
-  }
-}; // class Hexahedron
+    // --------------------------------------------------------------------------------
+    struct _volumeDef // holder of nodes of a volume mesh element
+    {
+      typedef void* _ptr;
+
+      struct _nodeDef
+      {
+        const SMDS_MeshNode*    _node; // mesh node at hexahedron corner
+        const StdMeshers::Cartesian3D::B_IntersectPoint* _intPoint;
+
+        _nodeDef(): _node(0), _intPoint(0) {}
+        _nodeDef( _Node* n ): _node( n->_node), _intPoint( n->_intPoint ) {}
+        const SMDS_MeshNode*    Node() const
+        { return ( _intPoint && _intPoint->_node ) ? _intPoint->_node : _node; }
+        const StdMeshers::Cartesian3D::E_IntersectPoint* EdgeIntPnt() const
+        { return static_cast< const StdMeshers::Cartesian3D::E_IntersectPoint* >( _intPoint ); }
+        _ptr Ptr() const { return Node() ? (_ptr) Node() : (_ptr) EdgeIntPnt(); }
+        bool operator==(const _nodeDef& other ) const { return Ptr() == other.Ptr(); }
+      };
+
+      std::vector< _nodeDef >              _nodes;
+      std::vector< int >                   _quantities;
+      _volumeDef*                          _next; // to store several _volumeDefs in a chain
+      StdMeshers::Cartesian3D::TGeomID     _solidID;
+      double                               _size;
+      const SMDS_MeshElement*              _volume; // new volume
+      std::vector<const SMDS_MeshElement*> _brotherVolume; // produced due to poly split 
+      std::vector< SMESH_Block::TShapeID > _names; // name of side a polygon originates from
+
+      _volumeDef(): _next(0), _solidID(0), _size(0), _volume(0) {}
+      ~_volumeDef() { delete _next; }
+      _volumeDef( _volumeDef& other ):
+        _next(0), _solidID( other._solidID ), _size( other._size ), _volume( other._volume )
+      { _nodes.swap( other._nodes ); _quantities.swap( other._quantities ); other._volume = 0;
+        _names.swap( other._names ); }
+
+      size_t size() const { return 1 + ( _next ? _next->size() : 0 ); } // nb _volumeDef in a chain
+      _volumeDef* at(int index)
+      { return index == 0 ? this : ( _next ? _next->at(index-1) : _next ); }
+
+      void Set( _Node** nodes, int nb )
+      { _nodes.assign( nodes, nodes + nb ); }
+
+      void SetNext( _volumeDef* vd )
+      { if ( _next ) { _next->SetNext( vd ); } else { _next = vd; }}
+
+      bool IsEmpty() const { return (( _nodes.empty() ) &&
+                                     ( !_next || _next->IsEmpty() )); }
+      bool IsPolyhedron() const { return ( !_quantities.empty() ||
+                                           ( _next && !_next->_quantities.empty() )); }
+
+
+      struct _linkDef: public std::pair<_ptr,_ptr> // to join polygons in removeExcessSideDivision()
+      {
+        _nodeDef _node1;//, _node2;
+        mutable /*const */_linkDef *_prev, *_next;
+        size_t _loopIndex;
+
+        _linkDef():_prev(0), _next(0) {}
+
+        void init( const _nodeDef& n1, const _nodeDef& n2, size_t iLoop )
+        {
+          _node1     = n1; //_node2 = n2;
+          _loopIndex = iLoop;
+          first      = n1.Ptr();
+          second     = n2.Ptr();
+          if ( first > second ) std::swap( first, second );
+        }
+        void setNext( _linkDef* next )
+        {
+          _next = next;
+          next->_prev = this;
+        }
+      };
+    };
+
+    // topology of a hexahedron
+    _Node _hexNodes [8];
+    _Link _hexLinks [12];
+    _Face _hexQuads [6];
+
+    // faces resulted from hexahedron intersection
+    std::vector< _Face > _polygons;
+
+    // intresections with EDGEs
+    std::vector< const StdMeshers::Cartesian3D::E_IntersectPoint* > _eIntPoints;
+
+    // additional nodes created at intersection points
+    std::vector< _Node > _intNodes;
+
+    // nodes inside the hexahedron (at VERTEXes) refer to _intNodes
+    std::vector< _Node* > _vIntNodes;
+
+    // computed volume elements
+    _volumeDef _volumeDefs;
+
+    StdMeshers::Cartesian3D::Grid*       _grid;
+    double      _sideLength[3];
+    int         _nbCornerNodes, _nbFaceIntNodes, _nbBndNodes;
+    int         _origNodeInd; // index of _hexNodes[0] node within the _grid
+    size_t      _i,_j,_k;
+    bool        _hasTooSmall;
+    int         _cellID;
+
+  public:
+    Hexahedron(StdMeshers::Cartesian3D::Grid* grid);
+    int MakeElements(SMESH_MesherHelper&                      helper,
+                     const TEdge2faceIDsMap& edge2faceIDsMap,
+                     const int numOfThreads = 1 );
+    void computeElements( const Solid* solid = 0, int solidIndex = -1 );
+
+  private:
+    Hexahedron(const Hexahedron& other, size_t i, size_t j, size_t k, int cellID );
+    void init( size_t i, size_t j, size_t k, const Solid* solid=0 );
+    void init( size_t i );
+    void setIJK( size_t i );
+    /*Auxiliary methods to extract operations from monolitic compute method*/
+    void defineHexahedralFaces( const Solid* solid, const IsInternalFlag intFlag );
+    bool compute( const Solid* solid, const IsInternalFlag intFlag );
+    size_t getSolids( StdMeshers::Cartesian3D::TGeomID ids[] );
+    bool isCutByInternalFace( IsInternalFlag & maxFlag );
+    void addEdges(SMESH_MesherHelper&         helper,
+                  std::vector< Hexahedron* >& intersectedHex,
+                  const TEdge2faceIDsMap&     edge2faceIDsMap);
+    gp_Pnt findIntPoint( double u1, double proj1, double u2, double proj2,
+                         double proj, BRepAdaptor_Curve& curve,
+                         const gp_XYZ& axis, const gp_XYZ& origin );
+    int  getEntity( const StdMeshers::Cartesian3D::E_IntersectPoint* ip, int* facets, int& sub );
+    bool addIntersection( const StdMeshers::Cartesian3D::E_IntersectPoint* ip,
+                          std::vector< Hexahedron* >&  hexes,
+                          int ijk[], int dIJK[] );
+    bool isQuadOnFace( const size_t iQuad );
+    bool findChain( _Node* n1, _Node* n2, _Face& quad, std::vector<_Node*>& chainNodes );
+    bool closePolygon( _Face* polygon, std::vector<_Node*>& chainNodes ) const;
+    bool findChainOnEdge( const std::vector< _OrientedLink >& splits,
+                          const _OrientedLink&                prevSplit,
+                          const _OrientedLink&                avoidSplit,
+                          const std::set< StdMeshers::Cartesian3D::TGeomID > & concaveFaces,
+                          size_t &                            iS,
+                          _Face&                              quad,
+                          std::vector<_Node*>&                chn);
+    typedef std::pair< StdMeshers::Cartesian3D::TGeomID, int > TFaceOfLink; // (face, link)
+    static TFaceOfLink findStartLink(const std::vector< _OrientedLink* >& freeLinks,
+                                     std::set< StdMeshers::Cartesian3D::TGeomID >& usedFaceIDs);
+    int  addVolumes( SMESH_MesherHelper& helper );
+    void addFaces( SMESH_MesherHelper& helper,
+                   const std::vector< const SMDS_MeshElement* > & boundaryVolumes );
+    void addSegments( SMESH_MesherHelper&     helper,
+                      const TEdge2faceIDsMap& edge2faceIDsMap );
+    void getVolumes( std::vector< const SMDS_MeshElement* > & volumes );
+    void getBoundaryElems( std::vector< const SMDS_MeshElement* > & boundaryVolumes );
+    void removeExcessSideDivision(const std::vector< Hexahedron* >& allHexa);
+    void removeExcessNodes(std::vector< Hexahedron* >& allHexa);
+    void preventVolumesOverlapping();
+    StdMeshers::Cartesian3D::TGeomID getAnyFace() const;
+    void cutByExtendedInternal( std::vector< Hexahedron* >& hexes,
+                                const TColStd_MapOfInteger& intEdgeIDs );
+    gp_Pnt mostDistantInternalPnt( int hexIndex, const gp_Pnt& p1, const gp_Pnt& p2 );
+    bool isOutPoint( _Link& link, int iP, SMESH_MesherHelper& helper, const Solid* solid ) const;
+    void sortVertexNodes(std::vector<_Node*>& nodes,
+                         _Node* curNode,
+                         StdMeshers::Cartesian3D::TGeomID face);
+    bool isInHole() const;
+    bool hasStrangeEdge() const;
+    bool checkPolyhedronSize( bool isCutByInternalFace, double & volSize ) const;
+    int checkPolyhedronValidity( _volumeDef* volDef, std::vector<std::vector<int>>& splitQuantities, 
+                                 std::vector<std::vector<const SMDS_MeshNode*>>& splitNodes );
+    const SMDS_MeshElement* addPolyhedronToMesh( _volumeDef* volDef,
+                                                 SMESH_MesherHelper& helper,
+                                                 const std::vector<const SMDS_MeshNode*>& nodes,
+                                                 const std::vector<int>& quantities );
+    bool addHexa ();
+    bool addTetra();
+    bool addPenta();
+    bool addPyra ();
+    bool debugDumpLink( _Link* link );
+    _Node* findEqualNode( std::vector< _Node* >&       nodes,
+                          const StdMeshers::Cartesian3D::E_IntersectPoint* ip,
+                          const double            tol2 )
+    {
+      for ( size_t i = 0; i < nodes.size(); ++i )
+        if ( nodes[i]->EdgeIntPnt() == ip ||
+             nodes[i]->Point().SquareDistance( ip->_point ) <= tol2 )
+          return nodes[i];
+      return 0;
+    }
+    bool isCorner( const _Node* node ) const { return ( node >= &_hexNodes[0] &&
+                                                        node -  &_hexNodes[0] < 8 ); }
+    bool hasEdgesAround( const ConcaveFace* cf ) const;
+    bool isImplementEdges() const { return _grid->_edgeIntPool.nbElements(); }
+    bool isOutParam(const double uvw[3]) const;
+
+    typedef boost::container::flat_map< StdMeshers::Cartesian3D::TGeomID, size_t > TID2Nb;
+    static void insertAndIncrement( StdMeshers::Cartesian3D::TGeomID id, TID2Nb& id2nbMap )
+    {
+      TID2Nb::value_type s0( id, 0 );
+      TID2Nb::iterator id2nb = id2nbMap.insert( s0 ).first;
+      id2nb->second++;
+    }
+  }; // class Hexahedron
+} // end namespace Cartesian3D
+} // end namespace StdMeshers
 
 #endif
index 3ed6ccc59fc3d182886a0bce643c164010760693..7f6412ea52dde48625bfe6e2578bd38a9b794f25 100644 (file)
@@ -1,8 +1,8 @@
 DBRep_DrawableShape
 
-CASCADE Topology V3, (c) Open Cascade
+CASCADE Topology V1, (c) Matra-Datavision
 Locations 0
-Curve2ds 59
+Curve2ds 54
 1 0 0 1 0 
 1 0 0 1 0 
 1 100 0 0 -1 
@@ -11,58 +11,53 @@ Curve2ds 59
 1 0 0 1 0 
 1 0 0 0 -1 
 1 0 0 0 1 
-1 0 0 1 0 
 1 0 60 1 0 
-1 0 0 0 1 
 1 0 0 1 0 
 1 16.899999999999999 50 0 1 
 1 0 0 1 0 
-1 0 0 1 0 
 1 16.899999999999999 50 1 0 
-1 59.899999999999999 0 0 1 
 1 0 0 1 0 
 1 76.799999999999997 50 0 1 
+1 0 0 1 0 
+1 0 60 1 0 
+1 0 0 1 0 
 1 100 0 0 1 
 1 0 0 1 0 
 1 0 0 0 1 
 1 0 0 1 0 
-1 100 0 0 1 
 1 0 50 1 0 
-1 100 0 0 -
+1 100 0 0 1 
 1 60 0 0 1 
-1 0 -50 1 0 
+1 100 0 0 -1 
 1 0 60 1 0 
-1 0 0 0 1 
-1 0 50 1 0 
-1 16.899999999999999 50 0 1 
 1 0 -50 1 0 
-1 0 0 1 0 
-1 16.899999999999999 50 1 0 
-1 59.899999999999999 0 0 1 
+1 16.899999999999999 50 0 1 
 1 0 50 1 0 
+1 16.899999999999999 50 1 0 
+1 0 -50 1 0 
 1 76.799999999999997 50 0 1 
+1 0 50 1 0 
+1 0 60 1 0 
+1 0 -50 1 0 
 1 0 0 0 1 
 1 0 50 1 0 
-1 0 0 0 -1 
 1 60 0 0 1 
-1 0 0 1 0 
-2 0 0 1 0 -0 1 10
-2 25 25 -1 0 0 1 10
 1 0 0 0 -1 
-1 10 0 0 1 
+2 25 25 -1 0 0 1 10
+1 0 0 1 0 
 1 16.899999999999999 0 0 -1 
-1 0 0 0 -
+1 10 0 0 
 1 0 0 0 1 
+1 0 0 0 -1 
 1 59.899999999999999 0 0 -1 
 1 0 0 0 1 
-1 59.899999999999999 0 0 -1 
 1 10 0 0 1 
 1 76.799999999999997 0 0 -1 
 1 0 50 1 0 
 2 0 0 1 0 -0 1 10
 1 6.2831853071795862 -0 0 1 
 1 0 -0 0 1 
-Curves 25
+Curves 27
 1 0 0 0 0 0 1 
 1 0 0 100 -0 1 0 
 1 0 50 0 0 0 1 
@@ -71,6 +66,7 @@ Curves 25
 1 50 0 16.899999999999999 1 0 -0 
 1 50 0 16.899999999999999 0 0 1 
 1 50 0 76.799999999999997 1 0 -0 
+1 60 0 0 0 0 1 
 1 0 0 100 1 0 -0 
 1 0 0 0 1 0 -0 
 1 0 50 100 1 0 -0 
@@ -79,6 +75,7 @@ Curves 25
 1 50 50 16.899999999999999 1 0 -0 
 1 50 50 16.899999999999999 0 0 1 
 1 50 50 76.799999999999997 1 0 -0 
+1 60 50 0 0 0 1 
 1 0 50 0 1 0 -0 
 1 60 0 0 -0 1 0 
 2 25 25 0 0 0 -1 -1 0 -0 0 1 0 10
@@ -89,156 +86,21 @@ Curves 25
 2 25 25 -50 0 0 -1 -1 0 -0 0 1 0 10
 1 15.000000000000002 24.999999999999996 0 0 0 -1 
 Polygon3D 0
-PolygonOnTriangulations 54
-2 1 2 
-p 0.6000000008 1 0 100 
-2 3 4 
-p 0.6000000008 1 0 100 
-2 2 4 
-p 0.6000000008 1 0 50 
-2 1 2 
-p 0.6000000008 1 0 50 
-2 3 4 
-p 0.6000000008 1 0 100 
-2 3 4 
-p 0.6000000008 1 0 100 
-2 1 3 
-p 0.6000000008 1 0 50 
-2 1 2 
-p 0.6000000008 1 0 50 
-2 1 2 
-p 0.6000000008 1 0 16.9 
-2 1 2 
-p 0.6000000008 1 0 16.9 
-2 8 2 
-p 0.6000000008 1 0 10 
-2 1 3 
-p 0.6000000008 1 0 10 
-2 8 7 
-p 0.6000000008 1 0 59.9 
-2 1 2 
-p 0.6000000008 1 0 59.9 
-2 7 6 
-p 0.6000000008 1 0 10 
-2 1 3 
-p 0.6000000008 1 0 10 
-2 6 5 
-p 0.6000000008 1 76.8 100 
-2 1 2 
-p 0.6000000008 1 76.8 100 
-2 4 5 
-p 0.6000000008 1 0 60 
-2 1 3 
-p 0.6000000008 1 0 60 
-2 3 1 
-p 0.6000000008 1 0 60 
-2 1 3 
-p 0.6000000008 1 0 60 
-2 2 4 
-p 0.6000000008 1 0 60 
-2 4 5 
-p 0.6000000008 1 0 60 
-2 3 4 
-p 0.6000000008 1 0 50 
-2 2 4 
-p 0.6000000008 1 0 50 
-2 1 2 
-p 0.6000000008 1 0 16.9 
-2 3 4 
-p 0.6000000008 1 0 16.9 
-2 8 2 
-p 0.6000000008 1 0 10 
-2 2 4 
-p 0.6000000008 1 0 10 
-2 8 7 
-p 0.6000000008 1 0 59.9 
-2 3 4 
-p 0.6000000008 1 0 59.9 
-2 7 6 
-p 0.6000000008 1 0 10 
-2 2 4 
-p 0.6000000008 1 0 10 
-2 6 5 
-p 0.6000000008 1 76.8 100 
-2 3 4 
-p 0.6000000008 1 76.8 100 
-2 3 1 
-p 0.6000000008 1 0 60 
-2 2 4 
-p 0.6000000008 1 0 60 
-2 3 4 
-p 0.6000000008 1 0 50 
-2 1 3 
-p 0.6000000008 1 0 50 
-37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 
-p 0.6000000008 1 0 0.174532925199433 0.349065850398866 0.523598775598299 0.698131700797732 0.872664625997165 1.0471975511966 1.22173047639603 1.39626340159546 1.5707963267949 1.74532925199433 1.91986217719376 2.0943951023932 2.26892802759263 2.44346095279206 2.61799387799149 2.79252680319093 2.96705972839036 3.14159265358979 3.31612557878922 3.49065850398866 3.66519142918809 3.83972435438752 4.01425727958696 4.18879020478639 4.36332312998582 4.53785605518525 4.71238898038469 4.88692190558412 5.06145483078355 5.23598775598299 5.41052068118242 5.58505360638185 5.75958653158128 5.93411945678072 6.10865238198015 6.28318530717959 
-37 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 5 
-p 0.6000000008 1 0 0.174532925199433 0.349065850398866 0.523598775598299 0.698131700797732 0.872664625997165 1.0471975511966 1.22173047639603 1.39626340159546 1.5707963267949 1.74532925199433 1.91986217719376 2.0943951023932 2.26892802759263 2.44346095279206 2.61799387799149 2.79252680319093 2.96705972839036 3.14159265358979 3.31612557878922 3.49065850398866 3.66519142918809 3.83972435438752 4.01425727958696 4.18879020478639 4.36332312998582 4.53785605518525 4.71238898038469 4.88692190558412 5.06145483078355 5.23598775598299 5.41052068118242 5.58505360638185 5.75958653158128 5.93411945678072 6.10865238198015 6.28318530717959 
-2 2 4 
-p 0.6000000008 1 0 50 
-2 3 4 
-p 0.6000000008 1 0 50 
-2 1 2 
-p 0.6000000008 1 0 50 
-2 1 3 
-p 0.6000000008 1 0 50 
-2 2 4 
-p 0.6000000008 1 0 50 
-2 1 2 
-p 0.6000000008 1 0 50 
-2 3 4 
-p 0.6000000008 1 0 50 
-2 1 3 
-p 0.6000000008 1 0 50 
-37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 
-p 0.6000000008 1 0 0.174532925199433 0.349065850398866 0.523598775598299 0.698131700797732 0.872664625997165 1.0471975511966 1.22173047639603 1.39626340159546 1.5707963267949 1.74532925199433 1.91986217719376 2.0943951023932 2.26892802759263 2.44346095279206 2.61799387799149 2.79252680319093 2.96705972839036 3.14159265358979 3.31612557878922 3.49065850398866 3.66519142918809 3.83972435438752 4.01425727958696 4.18879020478639 4.36332312998582 4.53785605518525 4.71238898038469 4.88692190558412 5.06145483078355 5.23598775598299 5.41052068118242 5.58505360638185 5.75958653158128 5.93411945678072 6.10865238198015 6.28318530717959 
-37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 1 
-p 0.6000000008 1 0 0.174532925199433 0.349065850398866 0.523598775598299 0.698131700797732 0.872664625997165 1.0471975511966 1.22173047639603 1.39626340159546 1.5707963267949 1.74532925199433 1.91986217719376 2.0943951023932 2.26892802759263 2.44346095279206 2.61799387799149 2.79252680319093 2.96705972839036 3.14159265358979 3.31612557878922 3.49065850398866 3.66519142918809 3.83972435438752 4.01425727958696 4.18879020478639 4.36332312998582 4.53785605518525 4.71238898038469 4.88692190558412 5.06145483078355 5.23598775598299 5.41052068118242 5.58505360638185 5.75958653158128 5.93411945678072 6.10865238198015 6.28318530717959 
-2 38 1 
-p 0.6000000008 1 0 50 
-2 74 37 
-p 0.6000000008 1 0 50 
-Surfaces 15
+PolygonOnTriangulations 0
+Surfaces 12
 1 0 0 0 1 0 -0 0 0 1 0 -1 0 
 1 0 0 0 -0 1 0 0 0 1 1 0 -0 
 1 0 0 100 0 0 1 1 0 -0 -0 1 0 
 1 0 50 0 -0 1 0 0 0 1 1 0 -0 
 1 0 0 0 0 0 1 1 0 -0 -0 1 0 
 1 60 0 0 1 0 -0 0 0 1 0 -1 0 
-1 50 0 16.899999999999999 -0 1 0 0 0 1 1 0 -0 
 1 50 0 16.899999999999999 0 0 1 1 0 -0 -0 1 0 
 1 50 0 16.899999999999999 1 0 -0 0 0 1 0 -1 0 
 1 50 0 76.799999999999997 0 0 1 1 0 -0 -0 1 0 
-1 50 50 16.899999999999999 -0 1 0 0 0 1 1 0 -
+1 60 0 0 1 0 -0 0 0 1 0 -1 
 2 25 25 0 0 0 -1 -1 0 -0 0 1 0 10
-1 25 25 0 0 0 -1 -1 0 -0 0 1 0 
-1 60 0 16.899999999999999 1 0 -0 0 0 1 0 -1 0 
 1 25 25 -50 0 0 -1 -1 0 -0 0 1 0 
-Triangulations 12
-4 2 1 0 0
-0 0 0 0 0 100 0 50 0 0 50 100 0 0 100 0 0 -50 100 -50 2 1 3 2 3 4 
-8 6 1 0 0
-60 0 0 60 0 16.9 0 0 0 0 0 100 60 0 100 60 0 76.8 50 0 76.8 50 0 16.9 0 60 16.9 60 0 0 100 0 100 60 76.8 60 76.8 50 16.9 50 8 1 3 2 1 8 7 3 4 7 8 3 5 6 7 5 7 4 
-4 2 1 0 0
-0 0 100 0 50 100 60 0 100 60 50 100 0 0 0 50 60 0 60 50 3 2 1 4 2 3 
-8 6 1 0 0
-60 50 0 60 50 16.9 0 50 0 0 50 100 60 50 100 60 50 76.8 50 50 76.8 50 50 16.9 0 60 16.9 60 0 0 100 0 100 60 76.8 60 76.8 50 16.9 50 8 1 3 2 1 8 7 3 4 7 8 3 5 6 7 5 7 4 
-40 40 1 0 1.77635683940025e-15
-0 0 0 0 50 0 60 0 0 60 50 0 15 25 0 15.1519224698779 26.7364817766693 0 15.6030737921409 28.4202014332567 0 16.3397459621556 30 0 17.3395555688102 31.4278760968654 0 18.5721239031346 32.6604444311898 0 20 33.6602540378444 0 21.5797985667433 34.3969262078591 0 23.2635182233307 34.8480775301221 0 25 35 0 26.7364817766693 34.8480775301221 0 28.4202014332567 34.3969262078591 0 30 33.6602540378444 0 31.4278760968654 32.6604444311898 0 32.6604444311898 31.4278760968654 0 33.6602540378444 30 0 34.3969262078591 28.4202014332567 0 34.8480775301221 26.7364817766693 0 35 25 0 34.8480775301221 23.2635182233307 0 34.3969262078591 21.5797985667433 0 33.6602540378444 20 0 32.6604444311898 18.5721239031346 0 31.4278760968654 17.3395555688102 0 30 16.3397459621556 0 28.4202014332567 15.6030737921409 0 26.7364817766693 15.1519224698779 0 25 15 0 23.2635182233307 15.1519224698779 0 21.5797985667433 15.6030737921409 0 20 16.3397459621556 0 18.5721239031346 17.3395555688102 0 17.3395555688102 18.5721239031346 0 16.3397459621556 20 0 15.6030737921409 21.5797985667433 0 15.1519224698779 23.2635182233307 0 0 0 0 50 60 0 60 50 15 25 15.1519224698779 26.7364817766693 15.6030737921409 28.4202014332567 16.3397459621556 30 17.3395555688102 31.4278760968654 18.5721239031346 32.6604444311898 20 33.6602540378444 21.5797985667433 34.3969262078591 23.2635182233307 34.8480775301221 25 35 26.7364817766693 34.8480775301221 28.4202014332567 34.3969262078591 30 33.6602540378444 31.4278760968654 32.6604444311898 32.6604444311898 31.4278760968654 33.6602540378444 30 34.3969262078591 28.4202014332567 34.8480775301221 26.7364817766693 35 25 34.8480775301221 23.2635182233307 34.3969262078591 21.5797985667433 33.6602540378444 20 32.6604444311898 18.5721239031346 31.4278760968654 17.3395555688102 30 16.3397459621556 28.4202014332567 15.6030737921409 26.7364817766693 15.1519224698779 25 15 23.2635182233307 15.1519224698779 21.5797985667433 15.6030737921409 20 16.3397459621556 18.5721239031346 17.3395555688102 17.3395555688102 18.5721239031346 16.3397459621556 20 15.6030737921409 21.5797985667433 15.1519224698779 23.2635182233307 36 37 1 38 1 37 35 36 1 39 1 38 34 35 1 40 1 39 33 34 1 5 1 40 32 33 1 31 32 1 2 5 6 2 6 7 2 7 8 2 8 9 2 1 5 10 2 9 11 2 10 12 2 11 13 2 12 3 24 25 3 25 26 3 26 27 3 27 28 3 28 29 3 29 30 3 30 31 3 31 1 14 2 13 23 24 3 15 2 14 4 2 15 4 15 16 4 16 17 4 17 18 4 18 19 4 19 20 4 20 21 4 21 22 4 22 23 4 23 3 
-4 2 1 0 0
-60 0 0 60 0 16.9 60 50 0 60 50 16.9 0 0 16.9 0 0 -50 16.9 -50 1 3 4 2 1 4 
-4 2 1 0 0
-50 0 16.9 50 50 16.9 60 0 16.9 60 50 16.9 0 0 0 50 10 0 10 50 4 2 1 4 1 3 
-4 2 1 0 7.105427357601e-15
-50 0 16.9 50 0 76.8 50 50 16.9 50 50 76.8 0 0 59.9 0 0 -50 59.9 -50 2 1 3 2 3 4 
-4 2 1 0 0
-50 0 76.8 50 50 76.8 60 0 76.8 60 50 76.8 0 0 0 50 10 0 10 50 4 2 1 4 1 3 
-4 2 1 0 0
-60 0 76.8 60 0 100 60 50 76.8 60 50 100 76.8 0 100 0 76.8 -50 100 -50 1 3 4 2 1 4 
-74 72 1 0 0.0380530190825497
-15 25 -50 15.1519224698779 26.7364817766693 -50 15.6030737921409 28.4202014332567 -50 16.3397459621556 30 -50 17.3395555688102 31.4278760968654 -50 18.5721239031346 32.6604444311898 -50 20 33.6602540378444 -50 21.5797985667433 34.3969262078591 -50 23.2635182233307 34.8480775301221 -50 25 35 -50 26.7364817766693 34.8480775301221 -50 28.4202014332567 34.3969262078591 -50 30 33.6602540378444 -50 31.4278760968654 32.6604444311898 -50 32.6604444311898 31.4278760968654 -50 33.6602540378444 30 -50 34.3969262078591 28.4202014332567 -50 34.8480775301221 26.7364817766693 -50 35 25 -50 34.8480775301221 23.2635182233307 -50 34.3969262078591 21.5797985667433 -50 33.6602540378444 20 -50 32.6604444311898 18.5721239031346 -50 31.4278760968654 17.3395555688102 -50 30 16.3397459621556 -50 28.4202014332567 15.6030737921409 -50 26.7364817766693 15.1519224698779 -50 25 15 -50 23.2635182233307 15.1519224698779 -50 21.5797985667433 15.6030737921409 -50 20 16.3397459621556 -50 18.5721239031346 17.3395555688102 -50 17.3395555688102 18.5721239031346 -50 16.3397459621556 20 -50 15.6030737921409 21.5797985667433 -50 15.1519224698779 23.2635182233307 -50 15 25 -50 15 25 0 15.1519224698779 26.7364817766693 0 15.6030737921409 28.4202014332567 0 16.3397459621556 30 0 17.3395555688102 31.4278760968654 0 18.5721239031346 32.6604444311898 0 20 33.6602540378444 0 21.5797985667433 34.3969262078591 0 23.2635182233307 34.8480775301221 0 25 35 0 26.7364817766693 34.8480775301221 0 28.4202014332567 34.3969262078591 0 30 33.6602540378444 0 31.4278760968654 32.6604444311898 0 32.6604444311898 31.4278760968654 0 33.6602540378444 30 0 34.3969262078591 28.4202014332567 0 34.8480775301221 26.7364817766693 0 35 25 0 34.8480775301221 23.2635182233307 0 34.3969262078591 21.5797985667433 0 33.6602540378444 20 0 32.6604444311898 18.5721239031346 0 31.4278760968654 17.3395555688102 0 30 16.3397459621556 0 28.4202014332567 15.6030737921409 0 26.7364817766693 15.1519224698779 0 25 15 0 23.2635182233307 15.1519224698779 0 21.5797985667433 15.6030737921409 0 20 16.3397459621556 0 18.5721239031346 17.3395555688102 0 17.3395555688102 18.5721239031346 0 16.3397459621556 20 0 15.6030737921409 21.5797985667433 0 15.1519224698779 23.2635182233307 0 15 25 0 0 50 0.174532925199433 50 0.349065850398866 50 0.523598775598299 50 0.698131700797732 50 0.872664625997165 50 1.0471975511966 50 1.22173047639603 50 1.39626340159546 50 1.5707963267949 50 1.74532925199433 50 1.91986217719376 50 2.0943951023932 50 2.26892802759263 50 2.44346095279206 50 2.61799387799149 50 2.79252680319093 50 2.96705972839036 50 3.14159265358979 50 3.31612557878922 50 3.49065850398866 50 3.66519142918809 50 3.83972435438752 50 4.01425727958696 50 4.18879020478639 50 4.36332312998582 50 4.53785605518525 50 4.71238898038469 50 4.88692190558412 50 5.06145483078355 50 5.23598775598299 50 5.41052068118242 50 5.58505360638185 50 5.75958653158128 50 5.93411945678072 50 6.10865238198015 50 6.28318530717959 50 0 0 0.174532925199433 0 0.349065850398866 0 0.523598775598299 0 0.698131700797732 0 0.872664625997165 0 1.0471975511966 0 1.22173047639603 0 1.39626340159546 0 1.5707963267949 0 1.74532925199433 0 1.91986217719376 0 2.0943951023932 0 2.26892802759263 0 2.44346095279206 0 2.61799387799149 0 2.79252680319093 0 2.96705972839036 0 3.14159265358979 0 3.31612557878922 0 3.49065850398866 0 3.66519142918809 0 3.83972435438752 0 4.01425727958696 0 4.18879020478639 0 4.36332312998582 0 4.53785605518525 0 4.71238898038469 0 4.88692190558412 0 5.06145483078355 0 5.23598775598299 0 5.41052068118242 0 5.58505360638185 0 5.75958653158128 0 5.93411945678072 0 6.10865238198015 0 6.28318530717959 0 2 1 38 2 38 39 3 39 40 3 2 39 4 40 41 4 3 40 5 41 42 5 42 43 5 4 41 6 43 44 6 5 43 7 6 44 8 7 44 8 44 45 9 8 45 9 45 46 10 9 46 10 46 47 11 10 47 11 47 48 12 11 48 12 48 49 13 12 49 13 49 50 14 13 50 14 50 51 15 14 51 15 51 52 16 15 52 16 52 53 17 16 53 17 53 54 18 17 54 18 54 55 18 55 56 19 18 56 20 56 57 20 19 56 21 57 58 21 20 57 22 58 59 22 21 58 23 59 60 23 22 59 24 60 61 24 23 60 25 61 62 25 24 61 26 62 63 26 25 62 27 63 64 27 64 65 27 26 63 28 27 65 29 65 66 29 28 65 30 66 67 30 29 66 31 67 68 31 30 67 32 31 68 32 68 69 33 32 69 33 69 70 34 33 70 34 70 71 35 34 71 35 71 72 36 35 72 36 72 73 37 36 73 37 73 74 
-36 34 1 0 7.94410929039127e-15
-15 25 -50 15.1519224698779 26.7364817766693 -50 15.6030737921409 28.4202014332567 -50 16.3397459621556 30 -50 17.3395555688102 31.4278760968654 -50 18.5721239031346 32.6604444311898 -50 20 33.6602540378444 -50 21.5797985667433 34.3969262078591 -50 23.2635182233307 34.8480775301221 -50 25 35 -50 26.7364817766693 34.8480775301221 -50 28.4202014332567 34.3969262078591 -50 30 33.6602540378444 -50 31.4278760968654 32.6604444311898 -50 32.6604444311898 31.4278760968654 -50 33.6602540378444 30 -50 34.3969262078591 28.4202014332567 -50 34.8480775301221 26.7364817766693 -50 35 25 -50 34.8480775301221 23.2635182233307 -50 34.3969262078591 21.5797985667433 -50 33.6602540378444 20 -50 32.6604444311898 18.5721239031346 -50 31.4278760968654 17.3395555688102 -50 30 16.3397459621556 -50 28.4202014332567 15.6030737921409 -50 26.7364817766693 15.1519224698779 -50 25 15 -50 23.2635182233307 15.1519224698779 -50 21.5797985667433 15.6030737921409 -50 20 16.3397459621556 -50 18.5721239031346 17.3395555688102 -50 17.3395555688102 18.5721239031346 -50 16.3397459621556 20 -50 15.6030737921409 21.5797985667433 -50 15.1519224698779 23.2635182233307 -50 10 -1.77635683940025e-15 9.84807753012208 1.7364817766693 9.39692620785909 3.42020143325669 8.66025403784439 5 7.66044443118978 6.42787609686539 6.42787609686539 7.66044443118978 5 8.66025403784439 3.42020143325669 9.39692620785908 1.7364817766693 9.84807753012208 0 10 -1.7364817766693 9.84807753012208 -3.42020143325669 9.39692620785909 -5 8.66025403784439 -6.42787609686539 7.66044443118978 -7.66044443118978 6.4278760968654 -8.66025403784438 5.00000000000001 -9.39692620785908 3.4202014332567 -9.84807753012208 1.73648177666932 -10 1.4210854715202e-14 -9.84807753012208 -1.73648177666929 -9.39692620785909 -3.42020143325667 -8.6602540378444 -4.99999999999998 -7.6604444311898 -6.42787609686538 -6.42787609686541 -7.66044443118977 -5.00000000000002 -8.66025403784437 -3.42020143325671 -9.39692620785907 -1.73648177666933 -9.84807753012208 -2.8421709430404e-14 -10 1.73648177666927 -9.84807753012209 3.42020143325666 -9.3969262078591 4.99999999999997 -8.6602540378444 6.42787609686536 -7.6604444311898 7.66044443118976 -6.42787609686542 8.66025403784437 -5.00000000000004 9.39692620785907 -3.42020143325673 9.84807753012207 -1.73648177666935 22 23 24 20 21 22 28 26 27 18 19 20 18 24 25 18 22 24 18 20 22 30 28 29 31 25 26 31 26 28 31 28 30 33 31 32 13 14 15 13 15 16 13 16 17 13 17 18 35 33 34 35 31 33 36 18 25 36 25 31 36 31 35 11 12 13 11 13 18 9 10 11 2 36 1 8 9 11 8 11 18 3 36 2 4 36 3 4 7 8 4 18 36 4 8 18 6 7 4 6 4 5 
+Triangulations 0
 
 TShapes 72
 Ve
@@ -260,8 +122,6 @@ Ed
 1  1 0 0 100
 2  1 1 0 0 100
 2  2 2 0 0 100
-6  1 1 0
-6  2 2 0
 0
 
 0101000
@@ -278,8 +138,6 @@ Ed
 1  2 0 0 50
 2  3 1 0 0 50
 2  4 3 0 0 50
-6  3 1 0
-6  4 3 0
 0
 
 0101000
@@ -296,8 +154,6 @@ Ed
 1  3 0 0 100
 2  5 1 0 0 100
 2  6 4 0 0 100
-6  5 1 0
-6  6 4 0
 0
 
 0101000
@@ -307,8 +163,6 @@ Ed
 1  4 0 0 50
 2  7 1 0 0 50
 2  8 5 0 0 50
-6  7 1 0
-6  8 5 0
 0
 
 0101000
@@ -319,7 +173,7 @@ Wi
 -70 0 -68 0 +66 0 +65 0 *
 Fa
 0  1e-07 1 0
-2  1
+
 0101000
 +64 0 *
 Ve
@@ -339,10 +193,8 @@ Ve
 Ed
  1e-07 1 1 0
 1  5 0 0 16.9
-2  9 6 0 0 16.9
-2  10 2 0 0 16.9
-6  9 2 0
-6  10 6 0
+2  9 2 0 0 16.9
+2  10 6 0 0 16.9
 0
 
 0101000
@@ -357,11 +209,8 @@ Ve
 Ed
  1e-07 1 1 0
 1  6 0 0 10
-2  11 7 0 0 10
-2  12 8 0 0 10
-2  13 2 0 0 10
-6  11 2 0
-6  12 7 0
+2  11 2 0 0 10
+2  12 7 0 0 10
 0
 
 0101000
@@ -376,11 +225,8 @@ Ve
 Ed
  1e-07 1 1 0
 1  7 0 0 59.9
-2  14 9 0 0 59.9
-2  15 7 0 0 59.9
-2  16 2 0 0 59.9
-6  13 2 0
-6  14 8 0
+2  13 2 0 0 59.9
+2  14 8 0 0 59.9
 0
 
 0101000
@@ -395,11 +241,8 @@ Ve
 Ed
  1e-07 1 1 0
 1  8 0 0 10
-2  17 7 0 0 10
-2  18 10 0 0 10
-2  19 2 0 0 10
-6  15 2 0
-6  16 9 0
+2  15 2 0 0 10
+2  16 9 0 0 10
 0
 
 0101000
@@ -413,33 +256,27 @@ Ve
 *
 Ed
  1e-07 1 1 0
-1  5 0 76.8 100
-2  9 6 0 76.8 100
-2  10 2 0 76.8 100
-6  17 2 0
-6  18 10 0
+1  9 0 76.8 100
+2  17 2 0 76.8 100
+2  18 10 0 76.8 100
 0
 
 0101000
 +55 0 -53 0 *
 Ed
  1e-07 1 1 0
-1  9 0 0 60
-2  20 2 0 0 60
-2  21 3 0 0 60
-6  19 2 0
-6  20 3 0
+1  10 0 0 60
+2  19 2 0 0 60
+2  20 3 0 0 60
 0
 
 0101000
 -53 0 +72 0 *
 Ed
  1e-07 1 1 0
-1  10 0 0 60
-2  22 2 0 0 60
-2  23 5 0 0 60
-6  21 2 0
-6  22 5 0
+1  11 0 0 60
+2  21 2 0 0 60
+2  22 5 0 0 60
 0
 
 0101000
@@ -450,7 +287,7 @@ Wi
 -60 0 +58 0 -56 0 -54 0 -52 0 +51 0 +70 0 -50 0 *
 Fa
 0  1e-07 2 0
-2  2
+
 0101000
 +49 0 *
 Ve
@@ -462,22 +299,18 @@ Ve
 *
 Ed
  1e-07 1 1 0
-1  11 0 0 60
+1  12 0 0 60
+2  23 3 0 0 60
 2  24 4 0 0 60
-2  25 3 0 0 60
-6  23 3 0
-6  24 4 0
 0
 
 0101000
 -47 0 +69 0 *
 Ed
  1e-07 1 1 0
-1  12 0 0 50
-2  26 6 0 0 50
-2  27 3 0 0 50
-6  25 3 0
-6  26 10 0
+1  13 0 0 50
+2  25 3 0 0 50
+2  26 10 0 0 50
 0
 
 0101000
@@ -488,7 +321,7 @@ Wi
 -68 0 -46 0 +45 0 +51 0 *
 Fa
 0  1e-07 3 0
-2  3
+
 0101000
 +44 0 *
 Ve
@@ -507,11 +340,9 @@ Ve
 *
 Ed
  1e-07 1 1 0
-1  13 0 0 16.9
+1  14 0 0 16.9
+2  27 4 0 0 16.9
 2  28 6 0 0 16.9
-2  29 4 0 0 16.9
-6  27 4 0
-6  28 6 0
 0
 
 0101000
@@ -525,12 +356,9 @@ Ve
 *
 Ed
  1e-07 1 1 0
-1  14 0 0 10
-2  30 11 0 0 10
-2  31 8 0 0 10
-2  32 4 0 0 10
-6  29 4 0
-6  30 7 0
+1  15 0 0 10
+2  29 4 0 0 10
+2  30 7 0 0 10
 0
 
 0101000
@@ -544,12 +372,9 @@ Ve
 *
 Ed
  1e-07 1 1 0
-1  15 0 0 59.9
-2  33 9 0 0 59.9
-2  34 11 0 0 59.9
-2  35 4 0 0 59.9
-6  31 4 0
-6  32 8 0
+1  16 0 0 59.9
+2  31 4 0 0 59.9
+2  32 8 0 0 59.9
 0
 
 0101000
@@ -563,34 +388,27 @@ Ve
 *
 Ed
  1e-07 1 1 0
-1  16 0 0 10
-2  36 11 0 0 10
-2  37 10 0 0 10
-2  38 4 0 0 10
-6  33 4 0
-6  34 9 0
+1  17 0 0 10
+2  33 4 0 0 10
+2  34 9 0 0 10
 0
 
 0101000
 -35 0 +37 0 *
 Ed
  1e-07 1 1 0
-1  13 0 76.8 100
-2  28 6 0 76.8 100
-2  29 4 0 76.8 100
-6  35 4 0
-6  36 10 0
+1  18 0 76.8 100
+2  35 4 0 76.8 100
+2  36 10 0 76.8 100
 0
 
 0101000
 +35 0 -47 0 *
 Ed
  1e-07 1 1 0
-1  17 0 0 60
-2  39 4 0 0 60
-2  40 5 0 0 60
-6  37 4 0
-6  38 5 0
+1  19 0 0 60
+2  37 4 0 0 60
+2  38 5 0 0 60
 0
 
 0101000
@@ -601,16 +419,14 @@ Wi
 -40 0 +38 0 -36 0 -34 0 -33 0 +46 0 +66 0 -32 0 *
 Fa
 0  1e-07 4 0
-2  4
+
 0101000
 +31 0 *
 Ed
  1e-07 1 1 0
-1  18 0 0 50
-2  41 6 0 0 50
-2  42 5 0 0 50
-6  39 5 0
-6  40 6 0
+1  20 0 0 50
+2  39 5 0 0 50
+2  40 6 0 0 50
 0
 
 0101000
@@ -628,12 +444,9 @@ Ve
 *
 Ed
  1e-07 1 1 0
-1  19 0 0 6.28318530717959
-2  43 12 0 0 6.28318530717959
-2  44 13 0 0 6.28318530717959
-2  45 5 0 0 6.28318530717959
-6  41 11 0
-6  42 5 0
+1  21 0 0 6.28318530717959
+2  41 5 0 0 6.28318530717959
+2  42 11 0 0 6.28318530717959
 0
 
 0101000
@@ -644,17 +457,14 @@ Wi
 +26 0 *
 Fa
 0  1e-07 5 0
-2  5
+
 0101000
 +28 0 +25 0 *
 Ed
  1e-07 1 1 0
-1  20 0 0 50
-2  46 14 0 0 50
-2  47 8 0 0 50
-2  48 6 0 0 50
-6  43 6 0
-6  44 7 0
+1  22 0 0 50
+2  43 6 0 0 50
+2  44 7 0 0 50
 0
 
 0101000
@@ -665,16 +475,14 @@ Wi
 -60 0 -23 0 +40 0 +29 0 *
 Fa
 0  1e-07 6 0
-2  6
+
 0101000
 +22 0 *
 Ed
  1e-07 1 1 0
-1  21 0 0 50
-2  49 9 0 0 50
-2  50 8 0 0 50
-6  45 7 0
-6  46 8 0
+1  23 0 0 50
+2  45 7 0 0 50
+2  46 8 0 0 50
 0
 
 0101000
@@ -684,17 +492,15 @@ Wi
 0101100
 -20 0 -38 0 +23 0 +58 0 *
 Fa
-0  1e-07 8 0
-2  7
+0  1e-07 7 0
+
 0101000
 +19 0 *
 Ed
  1e-07 1 1 0
-1  22 0 0 50
-2  51 9 0 0 50
-2  52 10 0 0 50
-6  47 8 0
-6  48 9 0
+1  24 0 0 50
+2  47 8 0 0 50
+2  48 9 0 0 50
 0
 
 0101000
@@ -704,18 +510,15 @@ Wi
 0101100
 -56 0 -17 0 +36 0 +20 0 *
 Fa
-0  1e-07 9 0
-2  8
+0  1e-07 8 0
+
 0101000
 +16 0 *
 Ed
  1e-07 1 1 0
-1  23 0 0 50
-2  53 14 0 0 50
-2  54 10 0 0 50
-2  55 6 0 0 50
-6  49 9 0
-6  50 10 0
+1  25 0 0 50
+2  49 9 0 0 50
+2  50 10 0 0 50
 0
 
 0101000
@@ -725,8 +528,8 @@ Wi
 0101100
 -17 0 -34 0 +14 0 +54 0 *
 Fa
-0  1e-07 10 0
-2  9
+0  1e-07 9 0
+
 0101000
 +13 0 *
 Wi
@@ -734,8 +537,8 @@ Wi
 0101100
 -52 0 -45 0 +33 0 +14 0 *
 Fa
-0  1e-07 6 0
-2  10
+0  1e-07 10 0
+
 0101000
 +11 0 *
 Ve
@@ -747,20 +550,17 @@ Ve
 *
 Ed
  1e-07 1 1 0
-1  24 0 0 6.28318530717959
-2  56 12 0 0 6.28318530717959
-2  57 15 0 0 6.28318530717959
-6  51 11 0
-6  52 12 0
+1  26 0 0 6.28318530717959
+2  51 11 0 0 6.28318530717959
+2  52 12 0 0 6.28318530717959
 0
 
 0101000
 +9 0 -9 0 *
 Ed
  1e-07 1 1 0
-1  25 0 0 50
-3  58 59CN 12 0 0 50
-7  53 54 11 0
+1  27 0 0 50
+3  53 54CN 11 0 0 50
 0
 
 0101000
@@ -770,8 +570,8 @@ Wi
 0101100
 -8 0 +7 0 +26 0 -7 0 *
 Fa
-0  1e-07 12 0
-2  11
+0  1e-07 11 0
+
 0101000
 +6 0 *
 Wi
@@ -779,8 +579,8 @@ Wi
 0101100
 +8 0 *
 Fa
-0  1e-07 15 0
-2  12
+0  1e-07 12 0
+
 0101000
 +4 0 *
 Sh
@@ -790,7 +590,9 @@ Sh
 +5 0 +3 0 *
 So
 
-0100000
+1100000
 +2 0 *
 
-+1 0 
\ No newline at end of file
++1 0 
+0
+