Merge branch 'kleontev/42217_horseshoe_with_bodyfitting_after_42002'

src/CMakeLists.txt

@@ -37,6 +37,7 @@ SET(SUBDIRS_COMMON
   SMESHClient
   SMESH_SWIG
   StdMeshers
+  StdMeshers.test
   StdMeshers_I
   SMESH_PY
   Tools

src/OBJECT/SMESH_DeviceActor.cxx

@@ -43,6 +43,7 @@
 #include <vtkObjectFactory.h>
 #include <vtkShrinkFilter.h>
 #include <vtkShrinkPolyData.h>
+#include <vtkTriangleFilter.h>
 
 #include <vtkProperty.h>
 #include <vtkPolyData.h>
@@ -109,6 +110,7 @@ SMESH_DeviceActor
   myMergeFilter = vtkMergeFilter::New();
 
   myGeomFilter = VTKViewer_GeometryFilter::New();
+  myTriangleFilter = vtkTriangleFilter::New();
 
   myTransformFilter = VTKViewer_TransformFilter::New();
 
@@ -153,6 +155,7 @@ SMESH_DeviceActor
   myFaceOrientation->Delete();
 
   myGeomFilter->Delete();
+  myTriangleFilter->Delete();
 
   myTransformFilter->Delete();
 
@@ -251,8 +254,10 @@ SMESH_DeviceActor
 
     anId++; // 3
     myGeomFilter->SetInputConnection( myPassFilter[ anId ]->GetOutputPort() );
+    myTriangleFilter->SetInputConnection(myGeomFilter->GetOutputPort());
 
     anId++; // 4
+    // myPassFilter[ anId ]->SetInputConnection( myTriangleFilter->GetOutputPort() );
     myPassFilter[ anId ]->SetInputConnection( myGeomFilter->GetOutputPort() );
     myPassFilter[ anId + 1 ]->SetInputConnection( myPassFilter[ anId ]->GetOutputPort() );
 

src/OBJECT/SMESH_DeviceActor.h

@@ -46,6 +46,7 @@ class vtkLookupTable;
 class vtkImplicitBoolean;
 class vtkPassThrough;
 class vtkPlaneCollection;
+class vtkTriangleFilter;
 
 class VTKViewer_Transform;
 class VTKViewer_TransformFilter;
@@ -182,6 +183,7 @@ class SMESHOBJECT_EXPORT SMESH_DeviceActor: public vtkLODActor{
 
   bool myStoreClippingMapping;
   VTKViewer_GeometryFilter *myGeomFilter;
+  vtkTriangleFilter* myTriangleFilter = nullptr;
   VTKViewer_TransformFilter *myTransformFilter;
   std::vector<vtkPassThrough*> myPassFilter;
 

src/SMDS/SMDS_FaceOfNodes.hxx

@@ -32,7 +32,7 @@
 class SMDS_EXPORT SMDS_FaceOfNodes: public SMDS_CellOfNodes
 {
  public:
-  void Print(std::ostream & OS) const;
+  virtual void Print(std::ostream & OS) const override;
   SMDS_FaceOfNodes(const SMDS_MeshNode* node1,
                    const SMDS_MeshNode* node2,
                    const SMDS_MeshNode* node3);

src/SMDS/SMDS_MeshElement.hxx

@@ -142,7 +142,7 @@ public:
 
   SMDS_Mesh* GetMesh() const;
 
-  void Print(std::ostream & OS) const;
+  virtual void Print(std::ostream & OS) const;
 
   friend SMDS_EXPORT std::ostream & operator <<(std::ostream & OS, const SMDS_MeshElement *);
   friend class SMDS_ElementFactory;

src/SMDS/SMDS_MeshNode.hxx

@@ -65,7 +65,7 @@ class SMDS_EXPORT SMDS_MeshNode: public SMDS_MeshElement
   virtual bool IsMediumNode(const SMDS_MeshNode* /*node*/) const  { return false; }
   virtual int  NbCornerNodes() const { return 1; }
 
-  void Print(std::ostream & OS) const;
+  virtual void Print(std::ostream & OS) const override;
 
  private:
 

src/SMDS/SMDS_PolygonalFaceOfNodes.hxx

@@ -47,7 +47,7 @@ class SMDS_EXPORT SMDS_PolygonalFaceOfNodes : public SMDS_CellOfNodes
   virtual int NbEdges() const;
   virtual int NbFaces() const;
 
-  virtual void Print (std::ostream & OS) const;
+  virtual void Print (std::ostream & OS) const override;
 
   virtual const SMDS_MeshNode* GetNode(const int ind) const;
 

src/SMDS/SMDS_VolumeOfNodes.hxx

@@ -62,7 +62,7 @@ class SMDS_EXPORT SMDS_VolumeOfNodes: public SMDS_CellOfNodes
                    const int            nbNodes);
   ~SMDS_VolumeOfNodes();
 
-  void Print(std::ostream & OS) const;
+  virtual void Print(std::ostream & OS) const override;
   int NbFaces() const;
   int NbNodes() const;
   int NbEdges() const;

src/SMESH/SMESH_Hypothesis.cxx

@@ -47,8 +47,12 @@ SMESH_Hypothesis::SMESH_Hypothesis(int hypId,
   _type           = PARAM_ALGO;
   _shapeType      = 0;  // to be set by algo with TopAbs_Enum
   _param_algo_dim = -1; // to be set by algo parameter
-  StudyContextStruct* myStudyContext = gen->GetStudyContext();
+  
-  myStudyContext->mapHypothesis[hypId] = this;
+  if ( _gen )
+  {
+    StudyContextStruct* myStudyContext = gen->GetStudyContext();
+    myStudyContext->mapHypothesis[hypId] = this;
+  }    
 }
 
 //=============================================================================
@@ -107,15 +111,17 @@ int SMESH_Hypothesis::GetShapeType() const
 void SMESH_Hypothesis::NotifySubMeshesHypothesisModification()
 {
   // for all meshes in study
-
+  if ( _gen )
-  StudyContextStruct* myStudyContext = _gen->GetStudyContext();
-  map<int, SMESH_Mesh*>::iterator itm;
-  for (itm = myStudyContext->mapMesh.begin();
-       itm != myStudyContext->mapMesh.end();
-       itm++)
   {
-    SMESH_Mesh* mesh = (*itm).second;
+    StudyContextStruct* myStudyContext = _gen->GetStudyContext();
-    mesh->NotifySubMeshesHypothesisModification( this );
+    map<int, SMESH_Mesh*>::iterator itm;
+    for (itm = myStudyContext->mapMesh.begin();
+        itm != myStudyContext->mapMesh.end();
+        itm++)
+    {
+      SMESH_Mesh* mesh = (*itm).second;
+      mesh->NotifySubMeshesHypothesisModification( this );
+    }
   }    
 }
 
@@ -148,13 +154,16 @@ void SMESH_Hypothesis::SetLibName(const char* theLibName)
 
 SMESH_Mesh* SMESH_Hypothesis::GetMeshByPersistentID(int id) const
 {
-  StudyContextStruct* myStudyContext = _gen->GetStudyContext();
+  if ( _gen )
-  map<int, SMESH_Mesh*>::iterator itm = myStudyContext->mapMesh.begin();
-  for ( ; itm != myStudyContext->mapMesh.end(); itm++)
   {
-    SMESH_Mesh* mesh = (*itm).second;
+    StudyContextStruct* myStudyContext = _gen->GetStudyContext();
-    if ( mesh->GetMeshDS()->GetPersistentId() == id )
+    map<int, SMESH_Mesh*>::iterator itm = myStudyContext->mapMesh.begin();
-      return mesh;
+    for ( ; itm != myStudyContext->mapMesh.end(); itm++)
+    {
+      SMESH_Mesh* mesh = (*itm).second;
+      if ( mesh->GetMeshDS()->GetPersistentId() == id )
+        return mesh;
+    }
   }
   return 0;
 }

src/SMESHGUI/SMESHGUI_AddMeshElementDlg.cxx

@@ -70,6 +70,8 @@
 #include <vtkPolyDataMapper.h>
 #include <vtkProperty.h>
 #include <vtkCellData.h>
+#include <vtkTriangleFilter.h>
+#include <vtkGeometryFilter.h>
 
 // Qt includes
 #include <QComboBox>
@@ -108,6 +110,8 @@ namespace SMESH
     SALOME_Actor* myFaceOrientation;
     vtkPolyDataMapper* myFaceOrientationDataMapper;
     SMESH_FaceOrientationFilter* myFaceOrientationFilter;
+    vtkSmartPointer<vtkGeometryFilter> myGeometryFilter = nullptr;
+    vtkSmartPointer<vtkTriangleFilter> myTriangleFilter = nullptr;
 
   public:
     TElementSimulation (SalomeApp_Application* theApplication)
@@ -120,9 +124,15 @@ namespace SMESH
 
       myGrid = vtkUnstructuredGrid::New();
 
+      myGeometryFilter = vtkSmartPointer<vtkGeometryFilter>::New();
+      myGeometryFilter->SetInputData(myGrid);
+
+      myTriangleFilter = vtkSmartPointer<vtkTriangleFilter>::New();
+      myTriangleFilter->SetInputConnection(myGeometryFilter->GetOutputPort());
+
       // Create and display actor
       myMapper = vtkDataSetMapper::New();
-      myMapper->SetInputData(myGrid);
+      myMapper->SetInputConnection(myTriangleFilter->GetOutputPort());
 
       myPreviewActor = SALOME_Actor::New();
       myPreviewActor->PickableOff();
@@ -148,7 +158,7 @@ namespace SMESH
 
       // Orientation of faces
       myFaceOrientationFilter = SMESH_FaceOrientationFilter::New();
-      myFaceOrientationFilter->SetInputData(myGrid);
+      myFaceOrientationFilter->SetInputConnection(myTriangleFilter->GetOutputPort());
 
       myFaceOrientationDataMapper = vtkPolyDataMapper::New();
       myFaceOrientationDataMapper->SetInputConnection(myFaceOrientationFilter->GetOutputPort());
@@ -209,6 +219,9 @@ namespace SMESH
       myGrid->InsertNextCell(theType,anIds);
       anIds->Delete();
 
+      myGeometryFilter->Update();
+      myTriangleFilter->Update();
+
       myGrid->Modified();
 
       SetVisibility(true, theActor->GetFacesOriented(), false);

src/StdMeshers.test/CMakeLists.txt

@@ -0,0 +1,70 @@
+# Copyright (C) 2012-2024  CEA, EDF, OPEN CASCADE
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+INCLUDE(tests.set)
+
+SALOME_GENERATE_TESTS_ENVIRONMENT(_test_env)
+
+SET(TEST_INSTALL_DIRECTORY ${SMESH_TEST_DIR}/other)
+
+INCLUDE_DIRECTORIES(
+  ${OpenCASCADE_INCLUDE_DIR}
+  ${Boost_INCLUDE_DIRS}
+  ${MEDCOUPLING_INCLUDE_DIRS}
+  ${CPPUNIT_INCLUDE_DIRS}
+  ${PROJECT_SOURCE_DIR}/src/StdMeshers
+  ${PROJECT_SOURCE_DIR}/src/SMESHUtils
+  ${PROJECT_SOURCE_DIR}/src/SMESH
+  ${PROJECT_SOURCE_DIR}/src/SMESHDS
+  ${PROJECT_SOURCE_DIR}/src/SMDS
+  ${PROJECT_SOURCE_DIR}/src/Controls
+)
+# additional preprocessor / compiler flags
+ADD_DEFINITIONS(
+  ${OpenCASCADE_DEFINITIONS}
+  ${BOOST_DEFINITIONS}
+  ${CPPUNIT_DEFINITIONS}
+)
+
+IF(SALOME_SMESH_USE_TBB)
+  SET(TBB_LIBS ${TBB_LIBRARIES})
+ENDIF(SALOME_SMESH_USE_TBB)
+
+FOREACH(_test ${UNIT_TESTS})
+  GET_FILENAME_COMPONENT(testname ${_test} NAME_WE)
+  SET(testname "TESTS_${testname}")
+
+  add_executable(${_test} ${_test}.cxx)
+  target_link_libraries(${_test} StdMeshers ${CPPUNIT_LIBRARIES} )
+
+  ADD_TEST(NAME ${testname}
+           COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/${_test} )
+  SET_TESTS_PROPERTIES(${testname} PROPERTIES ENVIRONMENT "${tests_env}" LABELS "tests")
+ENDFOREACH()
+
+IF(WIN32)
+  FOREACH(_test ${UNIT_TESTS})
+    INSTALL(FILES ${CMAKE_CURRENT_BINARY_DIR}/${CMAKE_BUILD_TYPE}/${_test}${CMAKE_EXECUTABLE_SUFFIX} PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ DESTINATION ${TEST_INSTALL_DIRECTORY})
+  ENDFOREACH()
+ELSE()
+  FOREACH(_test ${UNIT_TESTS})
+    INSTALL(FILES ${CMAKE_CURRENT_BINARY_DIR}/${_test} PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ DESTINATION ${TEST_INSTALL_DIRECTORY})
+  ENDFOREACH()
+ENDIF(WIN32)
+

src/StdMeshers.test/HexahedronCanonicalShapesTest.cxx

@@ -0,0 +1,245 @@
+// Copyright (C) 2016-2024  CEA, EDF
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+//  File   : HexahedronCanonicalShapesTest.cxx
+//  Module : SMESH
+//  Purpose: Implement unit tests for StdMeshers_Cartesian_3D_Hexahedron class to reproduce bugs that manifest in integration tests.
+//            The main difference between this unit test and integration tests is the fine grained control we have over the class methods and the hability to diagnose/solve bugs before the code goes into production enviroment.
+//            This test class can be used as reference for the development of future tests in other stdMesh algorithms
+
+#include "StdMeshers_Cartesian_3D_Hexahedron.hxx"
+#include "StdMeshers_CartesianParameters3D.hxx"
+
+// CPP TEST
+#include <cppunit/TestAssert.h>
+
+// OCC
+#include <BRep_Builder.hxx>
+#include <BRepTools.hxx>
+#include <BRepPrimAPI_MakeBox.hxx>
+#include <BRepPrimAPI_MakeCylinder.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+#include <BRepPrimAPI_MakeCone.hxx>
+
+#include <iostream>
+#include <memory>
+
+using namespace StdMeshers::Cartesian3D;
+
+// Helper functions!
+// Build Grid
+//      Require building mesh
+//      Require building shape.
+
+/*!
+  * \brief Mock mesh
+  */
+struct SMESH_Mesh_Test: public SMESH_Mesh
+{
+  SMESH_Mesh_Test() {
+    _isShapeToMesh = (_id = 0);
+    _meshDS  = new SMESHDS_Mesh( _id, true );
+  }
+};
+
+/*!
+  * \brief Mock Hypothesis
+  */
+struct CartesianHypo: public StdMeshers_CartesianParameters3D
+{
+  CartesianHypo() : StdMeshers_CartesianParameters3D(0/*zero hypoId*/, nullptr/*NULL generator*/)
+  {
+  }
+};
+
+/*!
+  * \brief Initialize the grid and intesersectors of grid with the geometry
+  */
+void GridInitAndIntersectWithShape (Grid& grid,
+                                    double gridSpacing,
+                                    double theSizeThreshold,
+                                    const TopoDS_Shape theShape,
+                                    TEdge2faceIDsMap& edge2faceIDsMap,
+                                    const int theNumOfThreads)
+{
+  // Canonical axes(i,j,k)
+  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;
+
+  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
+  */
+bool testShape (const TopoDS_Shape theShape,
+                const bool toAddEdges,
+                const bool toCreateFaces,
+                const double theGridSpacing,
+                const double theSizeThreshold,
+                const int theNbCreatedExpected)
+{
+  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;
+  grid._toCreateFaces = toCreateFaces;
+  grid._toConsiderInternalFaces = false;
+  grid._toUseThresholdForInternalFaces = false;
+  grid._toUseQuanta = false;
+  grid._sizeThreshold = theSizeThreshold;
+
+  TEdge2faceIDsMap edge2faceIDsMap;
+  GridInitAndIntersectWithShape( grid, theGridSpacing, theSizeThreshold,
+                                 theShape, edge2faceIDsMap, 1 );
+
+  SMESH_subMesh * aSubMesh = aMesh->GetSubMesh(theShape);
+  aSubMesh->DependsOn(); // init sub-meshes
+
+  Hexahedron hex( &grid );
+  int nbAdded = hex.MakeElements( helper, edge2faceIDsMap, 1 );
+  if (nbAdded != theNbCreatedExpected) {
+    std::stringstream buffer;
+    buffer << "Number of computed elements does not match: obtained " << nbAdded << " != expected " << theNbCreatedExpected;
+    //CPPUNIT_ASSERT_MESSAGE(buffer.str().c_str(), nbAdded == theNbCreatedExpected );
+    //MESSAGE(buffer.str().c_str());
+    //CppUnitTestFramework::Logger::WriteMessage(buffer.str().c_str());
+    std::cerr << buffer.str() << std::endl;
+    return false;
+  }
+
+  return true;
+}
+
+/*!
+  * \brief Test some primitive shapes
+  */
+bool testPrimitives()
+{
+  bool isOK = true;
+
+  // Test fitting of a box
+  BRepPrimAPI_MakeBox aMakeBox (10, 20, 30);
+  aMakeBox.Build();
+  CPPUNIT_ASSERT_MESSAGE( "Could not create the box!", aMakeBox.IsDone() );
+  TopoDS_Shape aShape = aMakeBox.Shape();
+
+  // Test exact fitting of a box
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/4, /*theNbCreatedExpected*/6))
+    isOK = false;
+  if (!testShape (aShape, /*toAddEdges*/true, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/4, /*theNbCreatedExpected*/6))
+    isOK = false;
+
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/true,
+                  /*gridSpacing*/10, /*theSizeThreshold*/4, /*theNbCreatedExpected*/6))
+    isOK = false;
+
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/5, /*theSizeThreshold*/4, /*theNbCreatedExpected*/48))
+    isOK = false;
+
+  // Test not exact fitting of a box
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/7, /*theSizeThreshold*/4, /*theNbCreatedExpected*/12))
+    isOK = false;
+
+  // Test fitting of a cylinder
+  gp_Ax2 anAxes (gp::Origin(), gp::DZ());
+  BRepPrimAPI_MakeCylinder aMakeCyl (anAxes, 20., 30.);
+  aMakeCyl.Build();
+  CPPUNIT_ASSERT_MESSAGE( "Could not create the cylinder!", aMakeCyl.IsDone() );
+  aShape = aMakeCyl.Shape();
+
+  // test for different threshold values
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/4, /*theNbCreatedExpected*/48))
+    isOK = false;
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/2, /*theNbCreatedExpected*/36))
+    isOK = false;
+
+  // Test fitting of a sphere
+  BRepPrimAPI_MakeSphere aMakeSph (anAxes, 30.);
+  aMakeSph.Build();
+  CPPUNIT_ASSERT_MESSAGE( "Could not create the sphere!", aMakeSph.IsDone() );
+  aShape = aMakeSph.Shape();
+
+  // test for different threshold values
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/4, /*theNbCreatedExpected*/136))
+    isOK = false;
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/2, /*theNbCreatedExpected*/88))
+    isOK = false;
+
+  // Test fitting of a cone
+  BRepPrimAPI_MakeCone aMakeCon (anAxes, 30., 0., 40.);
+  aMakeCon.Build();
+  CPPUNIT_ASSERT_MESSAGE( "Could not create the cone!", aMakeCon.IsDone() );
+  aShape = aMakeCon.Shape();
+
+  // test for different threshold values
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/100, /*theNbCreatedExpected*/72))
+    isOK = false;
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/4, /*theNbCreatedExpected*/40))
+    isOK = false;
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/1.5, /*theNbCreatedExpected*/32))
+    isOK = false;
+
+  // truncated cone
+  aMakeCon = BRepPrimAPI_MakeCone(anAxes, 30., 15., 20.);
+  aMakeCon.Build();
+  CPPUNIT_ASSERT_MESSAGE( "Could not create the cone!", aMakeCon.IsDone() );
+  aShape = aMakeCon.Shape();
+
+  // test for different threshold values
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/100, /*theNbCreatedExpected*/56))
+    isOK = false;
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/4, /*theNbCreatedExpected*/36))
+    isOK = false;
+  if (!testShape (aShape, /*toAddEdges*/false, /*toCreateFaces*/false,
+                  /*gridSpacing*/10, /*theSizeThreshold*/1.5, /*theNbCreatedExpected*/28))
+    isOK = false;
+
+  return isOK;
+}
+
+// Entry point for test
+int main()
+{
+  bool isOK = testPrimitives();
+  return isOK ? 0 : 1;
+}

src/StdMeshers.test/HexahedronTest.cxx

@@ -0,0 +1,203 @@
+// Copyright (C) 2016-2024  CEA, EDF
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+//  File   : HexahedronTest.cxx
+//  Module : SMESH
+//  Purpose: Implement unit tests for StdMeshers_Cartesian_3D_Hexahedron class to reproduce bugs that manifest in integration tests.
+//            The main difference between this unit test and integration tests is the fine grained control we have over the class methods and the hability to diagnose/solve bugs before the code goes into production enviroment.
+//            This test class can be used as reference for the development of future tests in other stdMesh algorithms
+
+#include "StdMeshers_Cartesian_3D_Hexahedron.hxx"
+#include "StdMeshers_CartesianParameters3D.hxx"
+
+// CPP TEST
+#include <cppunit/TestAssert.h>
+
+// OCC
+#include <BRep_Builder.hxx>
+#include <BRepTools.hxx>
+#include <BRepPrimAPI_MakeBox.hxx>
+#include <BRepPrimAPI_MakeCylinder.hxx>
+#include <BRepPrimAPI_MakeSphere.hxx>
+#include <BRepPrimAPI_MakeCone.hxx>
+
+#include <iostream>
+#include <memory>
+
+using namespace StdMeshers::Cartesian3D;
+
+// Helper functions!
+// Build Grid
+//      Require building mesh
+//      Require building shape.
+
+/*!
+  * \brief Mock mesh
+  */
+struct SMESH_Mesh_Test: public SMESH_Mesh
+{
+  SMESH_Mesh_Test() {
+    _isShapeToMesh = (_id = 0);
+    _meshDS  = new SMESHDS_Mesh( _id, true );
+  }
+};
+
+/*!
+  * \brief Mock Hypothesis
+  */
+struct CartesianHypo: public StdMeshers_CartesianParameters3D
+{
+  CartesianHypo() : StdMeshers_CartesianParameters3D(0/*zero hypoId*/, nullptr/*NULL generator*/)
+  {
+  }
+};
+
+/*!
+  * \brief Shape loader
+  */
+void loadBrepShape( std::string shapeName, TopoDS_Shape & shape )
+{
+  BRep_Builder b;
+  BRepTools::Read(shape, shapeName.c_str(), b);
+}
+
+/*!
+  * \brief Initialize the grid and intesersectors of grid with the geometry
+  */
+void GridInitAndIntersectWithShape (Grid& grid,
+                                    double gridSpacing,
+                                    double theSizeThreshold,
+                                    const TopoDS_Shape theShape,
+                                    TEdge2faceIDsMap& edge2faceIDsMap,
+                                    const int theNumOfThreads)
+{
+  // Canonical axes(i,j,k)
+  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;
+
+  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);
+}
+
+/*!
+  * \brief Reproduce conditions of TBPERF_GRIDS_PERF_SMESH_M1 test to detect and solve segfault in unit test.
+  */
+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);
+
+  for (auto nThreads : numberOfThreads )
+  {
+    for (size_t i = 0; i < 10 /*trials*/; i++)
+    {
+      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._sizeThreshold = 4.0;
+
+      TEdge2faceIDsMap edge2faceIDsMap;
+      GridInitAndIntersectWithShape( grid, 1.0, 4.0, aShape, edge2faceIDsMap, nThreads );
+
+      SMESH_subMesh * aSubMesh = aMesh->GetSubMesh(aShape);
+      aSubMesh->DependsOn(); // init sub-meshes
+
+      Hexahedron hex( &grid );
+      int nbAdded = hex.MakeElements( helper, edge2faceIDsMap, nThreads );
+      CPPUNIT_ASSERT_MESSAGE( "Number of computed elements does not match", nbAdded == 1024 );
+    }
+  }
+  return true;
+}
+
+/*!
+  * \brief Reproduce conditions of TBPERF_GRIDS_PERF_SMESH_J4 test to detect and solve segfault in unit test.
+  */
+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);
+
+  // Test with face creation
+  for (auto nThreads : numberOfThreads )
+  {
+    for (size_t i = 0; i < 10 /*trials*/; i++)
+    {
+      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;
+      double testThreshold = 1.000001;
+      grid._toCreateFaces = true;
+      grid._sizeThreshold = testThreshold;
+
+      TEdge2faceIDsMap edge2faceIDsMap;
+      GridInitAndIntersectWithShape( grid, 2.0, testThreshold,
+                                     aShape, edge2faceIDsMap, nThreads );
+
+      SMESH_subMesh * aSubMesh = aMesh->GetSubMesh(aShape);
+      aSubMesh->DependsOn(); // init sub-meshes
+
+      Hexahedron hex( &grid );
+      int nbAdded = hex.MakeElements( helper, edge2faceIDsMap, nThreads );
+      CPPUNIT_ASSERT_MESSAGE( "Number of computed elements does not match", nbAdded == 35150 );
+    }
+  }
+  return true;
+}
+
+// Entry point for test
+int main()
+{
+  bool isOK = testNRTM1();
+  if (!testNRTJ4())
+    isOK = false;
+
+  return isOK ? 0 : 1;
+}

src/StdMeshers.test/tests.set

@@ -0,0 +1,28 @@
+# Copyright (C) 2015-2024  CEA, EDF, OPEN CASCADE
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+# The following tests cannot be executed with 'make test' because they use
+# external meshing plug-ins.
+# On the other hand these tests can be executed with 'salome test'.
+# ---------------------------------------------------------------------------
+
+SET(UNIT_TESTS
+  HexahedronTest
+  HexahedronCanonicalShapesTest
+  )

src/StdMeshers/CMakeLists.txt

@@ -119,6 +119,8 @@ SET(StdMeshers_HEADERS
   StdMeshers_ViscousLayers2D.hxx
   StdMeshers_Projection_1D2D.hxx
   StdMeshers_CartesianParameters3D.hxx
+  StdMeshers_Cartesian_3D_Grid.hxx
+  StdMeshers_Cartesian_3D_Hexahedron.hxx
   StdMeshers_Cartesian_3D.hxx
   StdMeshers_Cartesian_VL.hxx
   StdMeshers_QuadFromMedialAxis_1D2D.hxx
@@ -184,6 +186,8 @@ SET(StdMeshers_SOURCES
   StdMeshers_ViscousLayers2D.cxx
   StdMeshers_Projection_1D2D.cxx
   StdMeshers_CartesianParameters3D.cxx
+  StdMeshers_Cartesian_3D_Grid.cxx
+  StdMeshers_Cartesian_3D_Hexahedron.cxx
   StdMeshers_Cartesian_3D.cxx
   StdMeshers_Cartesian_VL.cxx
   StdMeshers_Adaptive1D.cxx

src/StdMeshers/StdMeshers_CartesianParameters3D.cxx

@@ -361,9 +361,13 @@ void StdMeshers_CartesianParameters3D::ComputeCoordinates(const double    x0,
         ++iCell;
       }
     }
-    const double lastCellLen = coords.back() - coords[ coords.size() - 2 ];
+    if (coords.size() < 2)
-    if ( fabs( coords.back() - p1 ) > 0.5 * lastCellLen )
       coords.push_back ( p1 );
+    else {
+      const double lastCellLen = coords.back() - coords[ coords.size() - 2 ];
+      if ( fabs( coords.back() - p1 ) > 0.5 * lastCellLen )
+        coords.push_back ( p1 );
+    }
   }
 
   // correct coords if a forced point is too close to a neighbor node

src/StdMeshers/StdMeshers_Cartesian_3D_Grid.hxx

@@ -0,0 +1,702 @@
+// Copyright (C) 2016-2024  CEA, EDF
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+//  File   : StdMeshers_Cartesian_3D_Grid.hxx
+//  Module : SMESH
+//  Purpose: Make BodyFitting mesh algorithm more modular and testable
+//
+
+#ifndef _SMESH_Cartesian_3D_GRID_HXX_
+#define _SMESH_Cartesian_3D_GRID_HXX_
+
+#ifdef WITH_TBB
+#include <tbb/parallel_for.h>
+#endif
+
+#include <utilities.h>
+
+// STD
+#include <algorithm>
+#include <map>
+#include <vector>
+#include <memory>
+#include <mutex>
+#include <thread>
+
+// SMESH
+#include "SMESH_StdMeshers.hxx"
+#include "StdMeshers_FaceSide.hxx"
+#include "StdMeshers_CartesianParameters3D.hxx"
+
+#include <ObjectPool.hxx>
+#include <SMDS_LinearEdge.hxx>
+#include <SMDS_MeshNode.hxx>
+#include <SMDS_VolumeOfNodes.hxx>
+#include <SMDS_VolumeTool.hxx>
+#include <SMESHDS_Mesh.hxx>
+#include <SMESH_Block.hxx>
+#include <SMESH_Comment.hxx>
+#include <SMESH_ControlsDef.hxx>
+#include <SMESH_Mesh.hxx>
+#include <SMESH_MeshAlgos.hxx>
+#include <SMESH_MeshEditor.hxx>
+#include <SMESH_MesherHelper.hxx>
+#include <SMESH_subMesh.hxx>
+#include <SMESH_subMeshEventListener.hxx>
+
+#include <utilities.h>
+#include <Utils_ExceptHandlers.hxx>
+
+#include <GEOMUtils.hxx>
+
+//OCC
+#include <BRepAdaptor_Curve.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepBndLib.hxx>
+#include <BRepBuilderAPI_Copy.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepTools.hxx>
+#include <BRepTopAdaptor_FClass2d.hxx>
+#include <BRep_Builder.hxx>
+#include <BRep_Tool.hxx>
+#include <Bnd_B3d.hxx>
+#include <Bnd_Box.hxx>
+#include <ElSLib.hxx>
+#include <GCPnts_UniformDeflection.hxx>
+#include <Geom2d_BSplineCurve.hxx>
+#include <Geom2d_BezierCurve.hxx>
+#include <Geom2d_TrimmedCurve.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+#include <GeomLib.hxx>
+#include <Geom_BSplineCurve.hxx>
+#include <Geom_BSplineSurface.hxx>
+#include <Geom_BezierCurve.hxx>
+#include <Geom_BezierSurface.hxx>
+#include <Geom_RectangularTrimmedSurface.hxx>
+#include <Geom_TrimmedCurve.hxx>
+#include <IntAna_IntConicQuad.hxx>
+#include <IntAna_IntLinTorus.hxx>
+#include <IntAna_Quadric.hxx>
+#include <IntCurveSurface_TransitionOnCurve.hxx>
+#include <IntCurvesFace_Intersector.hxx>
+#include <Poly_Triangulation.hxx>
+#include <Precision.hxx>
+#include <TopExp.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopLoc_Location.hxx>
+#include <TopTools_DataMapOfShapeInteger.hxx>
+#include <TopTools_IndexedMapOfShape.hxx>
+#include <TopTools_MapOfShape.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Compound.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_TShape.hxx>
+#include <gp_Cone.hxx>
+#include <gp_Cylinder.hxx>
+#include <gp_Lin.hxx>
+#include <gp_Pln.hxx>
+#include <gp_Pnt2d.hxx>
+#include <gp_Sphere.hxx>
+#include <gp_Torus.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
+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;
+
+  //=============================================================================
+  // Definitions of internal utils
+  // --------------------------------------------------------------------------
+  enum Transition {
+    Trans_TANGENT = IntCurveSurface_Tangent,
+    Trans_IN      = IntCurveSurface_In,
+    Trans_OUT     = IntCurveSurface_Out,
+    Trans_APEX,
+    Trans_INTERNAL // for INTERNAL FACE
+  };
+
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Sub-entities of a FACE neighboring its concave VERTEX.
+   *        Help to avoid linking nodes on EDGEs that seem connected
+   *        by the concave FACE but the link actually lies outside the FACE
+   */
+  struct ConcaveFace
+  {
+    TGeomID _concaveFace;
+    TGeomID _edge1, _edge2;
+    TGeomID _v1,    _v2;
+    ConcaveFace( int f=0, int e1=0, int e2=0, int v1=0, int v2=0 )
+      : _concaveFace(f), _edge1(e1), _edge2(e2), _v1(v1), _v2(v2) {}
+    bool HasEdge( TGeomID edge ) const { return edge == _edge1 || edge == _edge2; }
+    bool HasVertex( TGeomID v  ) const { return v == _v1 || v == _v2; }
+    void SetEdge( TGeomID edge ) { ( _edge1 ? _edge2 : _edge1 ) = edge; }
+    void SetVertex( TGeomID v  ) { ( _v1 ? _v2 : _v1 ) = v; }
+  };
+  typedef NCollection_DataMap< TGeomID, ConcaveFace > TConcaveVertex2Face;
+
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Container of IDs of SOLID sub-shapes
+   */
+  class Solid // sole SOLID contains all sub-shapes
+  {
+    TGeomID             _id; // SOLID id
+    bool                _hasInternalFaces;
+    TConcaveVertex2Face _concaveVertex; // concave VERTEX -> ConcaveFace
+  public:
+    virtual ~Solid() {}
+    virtual bool Contains( TGeomID /*subID*/ ) 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; }
+    TGeomID ID() const { return _id; }
+    void SetHasInternalFaces( bool has ) { _hasInternalFaces = has; }
+    bool HasInternalFaces() const { return _hasInternalFaces; }
+    void SetConcave( TGeomID V, TGeomID F, TGeomID E1, TGeomID E2, TGeomID V1, TGeomID V2  )
+    { _concaveVertex.Bind( V, ConcaveFace{ F, E1, E2, V1, V2 }); }
+    bool HasConcaveVertex() const { return !_concaveVertex.IsEmpty(); }
+    const ConcaveFace* GetConcave( TGeomID V ) const { return _concaveVertex.Seek( V ); }
+  };
+
+  // --------------------------------------------------------------------------
+  class OneOfSolids : public Solid
+  {
+    TColStd_MapOfInteger _subIDs;
+    TopTools_MapOfShape  _faces; // keep FACE orientation
+    TColStd_MapOfInteger _outFaceIDs; // FACEs of shape_to_mesh oriented outside the SOLID
+  public:
+    void Init( const TopoDS_Shape& solid,
+               TopAbs_ShapeEnum    subType,
+               const SMESHDS_Mesh* mesh );
+    virtual bool Contains( TGeomID i ) const { return i == ID() || _subIDs.Contains( i ); }
+    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;
+    }
+    virtual TopAbs_Orientation Orientation( const TopoDS_Shape& face ) const
+    {
+      const TopoDS_Shape& sInMap = const_cast< OneOfSolids* >(this)->_faces.Added( face );
+      return sInMap.Orientation();
+    }
+    virtual bool IsOutsideOriented( TGeomID faceID ) const
+    {
+      return faceID == 0 || _outFaceIDs.Contains( faceID );
+    }
+  };
+
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Hold a vector of TGeomID and clear it at destruction
+   */
+  class GeomIDVecHelder
+  {
+    typedef std::vector< TGeomID > TVector;
+    const TVector& myVec;
+    bool           myOwn;
+
+  public:
+    GeomIDVecHelder( const TVector& idVec, bool isOwner ): myVec( idVec ), myOwn( isOwner ) {}
+    GeomIDVecHelder( const GeomIDVecHelder& holder ): myVec( holder.myVec ), myOwn( holder.myOwn )
+    {
+      const_cast< bool& >( holder.myOwn ) = false;
+    }
+    ~GeomIDVecHelder() { if ( myOwn ) const_cast<TVector&>( myVec ).clear(); }
+    size_t size() const { return myVec.size(); }
+    TGeomID operator[]( size_t i ) const { return i < size() ? myVec[i] : theUndefID; }
+    bool operator==( const GeomIDVecHelder& other ) const { return myVec == other.myVec; }
+    bool contain( const TGeomID& id ) const {
+      return std::find( myVec.begin(), myVec.end(), id ) != myVec.end();
+    }
+    TGeomID otherThan( const TGeomID& id ) const {
+      for ( const TGeomID& id2 : myVec )
+        if ( id != id2 )
+          return id2;
+      return theUndefID;
+    }
+    TGeomID oneCommon( const GeomIDVecHelder& other ) const {
+      TGeomID common = theUndefID;
+      for ( const TGeomID& id : myVec )
+        if ( other.contain( id ))
+        {
+          if ( common != theUndefID )
+            return theUndefID;
+          common = id;
+        }
+      return common;
+    }
+  };
+
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Geom data
+   */
+  struct Geometry
+  {
+    TopoDS_Shape                _mainShape;
+    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
+
+    SMESH::Controls::ElementsOnShape _edgeClassifier;
+    SMESH::Controls::ElementsOnShape _vertexClassifier;
+
+    bool IsOneSolid() const { return _solidByID.size() < 2; }
+    GeomIDVecHelder GetSolidIDsByShapeID( const std::vector< TGeomID >& shapeIDs ) const;
+  };
+
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Common data of any intersection between a Grid and a shape
+   */
+  struct B_IntersectPoint
+  {
+    // This two class members are being updated in a non thread safe way.
+    // 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 std::vector< TGeomID > _faceIDs;
+
+    B_IntersectPoint(): _node(NULL) {}
+    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
+   */
+  struct F_IntersectPoint : public B_IntersectPoint
+  {
+    double             _paramOnLine;
+    double             _u, _v;
+    mutable Transition _transition;
+    mutable size_t     _indexOnLine;
+
+    bool operator< ( const F_IntersectPoint& o ) const { 
+      return _paramOnLine < o._paramOnLine;         
+    }
+  };
+
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Data of intersection between GridPlanes and a TopoDS_EDGE
+   */
+  struct E_IntersectPoint : public B_IntersectPoint
+  {
+    gp_Pnt  _point;
+    double  _uvw[3];
+    TGeomID _shapeID; // ID of EDGE or VERTEX
+  };
+
+  // --------------------------------------------------------------------------
+  /*!
+    * \brief A line of the grid and its intersections with 2D geometry
+    */
+  struct GridLine
+  {
+    gp_Lin _line;
+    double _length; // line length
+    std::multiset< F_IntersectPoint > _intPoints;
+
+    void RemoveExcessIntPoints( const double tol );
+    TGeomID GetSolidIDBefore( std::multiset< F_IntersectPoint >::iterator ip,
+                              const TGeomID                          prevID,
+                              const Geometry&                        geom);
+  };
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Planes of the grid used to find intersections of an EDGE with a hexahedron
+   */
+  struct GridPlanes
+  {
+    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
+  };
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Iterator on the parallel grid lines of one direction
+   */
+  struct LineIndexer
+  {
+    size_t _size  [3];
+    size_t _curInd[3];
+    size_t _iVar1, _iVar2, _iConst;
+    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 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;
+      _iVar1 = iv1; _iVar2 = iv2; _iConst = iConst;
+      _name1 = nv1; _name2 = nv2; _nameConst = nConst;
+    }
+
+    size_t I() const { return _curInd[0]; }
+    size_t J() const { return _curInd[1]; }
+    size_t K() const { return _curInd[2]; }
+    void SetIJK( size_t i, size_t j, size_t k )
+    {
+      _curInd[0] = i; _curInd[1] = j; _curInd[2] = k;
+    }
+    void SetLineIndex(size_t i)
+    {
+      _curInd[_iVar2] = i / _size[_iVar1];
+      _curInd[_iVar1] = i % _size[_iVar1];
+    }
+    void operator++()
+    {
+      if ( ++_curInd[_iVar1] == _size[_iVar1] )
+        _curInd[_iVar1] = 0, ++_curInd[_iVar2];
+    }
+    bool More() const { return _curInd[_iVar2] < _size[_iVar2]; }
+    size_t LineIndex   () const { return _curInd[_iVar1] + _curInd[_iVar2]* _size[_iVar1]; }
+    size_t LineIndex10 () const { return (_curInd[_iVar1] + 1 ) + _curInd[_iVar2]* _size[_iVar1]; }
+    size_t LineIndex01 () const { return _curInd[_iVar1] + (_curInd[_iVar2] + 1 )* _size[_iVar1]; }
+    size_t LineIndex11 () const { return (_curInd[_iVar1] + 1 ) + (_curInd[_iVar2] + 1 )* _size[_iVar1]; }
+    void SetIndexOnLine (size_t i)  { _curInd[ _iConst ] = i; }
+    bool IsValidIndexOnLine (size_t i) const { return  i < _size[ _iConst ]; }
+    size_t NbLines() const { return _size[_iVar1] * _size[_iVar2]; }
+  };
+
+  class Tools 
+  { 
+    public:
+      Tools() = delete;
+      
+      //================================================================================
+      /*!
+      * \brief computes exact bounding box with axes parallel to given ones
+      */
+      //================================================================================
+      static void GetExactBndBox( const std::vector< TopoDS_Shape >& faceVec, const double* axesDirs, Bnd_Box& shapeBox );
+  };
+
+  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];
+
+    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
+
+    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;
+    bool                              _toAddEdges;
+    bool                              _toCreateFaces;
+    bool                              _toConsiderInternalFaces;
+    bool                              _toUseThresholdForInternalFaces;
+    double                            _sizeThreshold;
+    bool                              _toUseQuanta;
+    double                            _quanta;
+
+    SMESH_MesherHelper*               _helper;
+
+    size_t CellIndex( size_t i, size_t j, size_t k ) const
+    {
+      return i + j*(_coords[0].size()-1) + k*(_coords[0].size()-1)*(_coords[1].size()-1);
+    }
+    size_t NodeIndex( size_t i, size_t j, size_t k ) const
+    {
+      return i + j*_coords[0].size() + k*_coords[0].size()*_coords[1].size();
+    }
+    size_t NodeIndex( const TIJK& ijk ) const
+    {
+      return NodeIndex( ijk[0], ijk[1], ijk[2] );
+    }
+    size_t NodeIndexDX() const { return 1; }
+    size_t NodeIndexDY() const { return _coords[0].size(); }
+    size_t NodeIndexDZ() const { return _coords[0].size() * _coords[1].size(); }
+
+    LineIndexer GetLineIndexer(size_t iDir) const;
+    size_t GetLineDir( const GridLine* line, size_t & index ) const;
+
+    E_IntersectPoint* Add( const E_IntersectPoint& ip )
+    {
+      E_IntersectPoint* eip = _edgeIntPool.getNew();
+      *eip = ip;
+      return eip;
+    }
+    void Remove( E_IntersectPoint* eip ) { _edgeIntPool.destroy( eip ); }
+
+    TGeomID ShapeID( const TopoDS_Shape& s ) const;
+    const TopoDS_Shape& Shape( TGeomID id ) const;
+    TopAbs_ShapeEnum ShapeType( TGeomID id ) const { return Shape(id).ShapeType(); }
+    void InitGeometry( const TopoDS_Shape& theShape );
+    void InitClassifier( const TopoDS_Shape&        mainShape,
+                        TopAbs_ShapeEnum           shapeType,
+                        SMESH::Controls::ElementsOnShape& classifier );
+    void GetEdgesToImplement( std::map< TGeomID, std::vector< TGeomID > > & edge2faceMap,
+                              const TopoDS_Shape&                 shape,
+                              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;
+    bool IsInternal( TGeomID faceID ) const {
+      return ( faceID == PseudoIntExtFaceID() ||
+              Shape( faceID ).Orientation() == TopAbs_INTERNAL ); }
+    bool IsSolid( TGeomID shapeID ) const {
+      if ( _geometry.IsOneSolid() ) return _geometry._soleSolid.ID() == shapeID;
+      else                          return _geometry._solidByID.count( shapeID ); }
+    bool IsStrangeEdge( TGeomID id ) const { return _geometry._strangeEdges.Contains( id ); }
+    TGeomID PseudoIntExtFaceID() const { return _geometry._extIntFaceID; }
+    Solid* GetSolid( TGeomID solidID = 0 );
+    Solid* GetOneOfSolids( TGeomID solidID );
+    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,
+                    TopoDS_Vertex* vertex = nullptr, bool unset = false );
+    void UpdateFacesOfVertex( const B_IntersectPoint& ip, const TopoDS_Vertex& vertex );
+    bool IsToCheckNodePos() const { return !_toAddEdges && _toCreateFaces; }
+    bool IsToRemoveExcessEntities() const { return !_toAddEdges; }
+
+    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, 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>
+  void parallel_for(const Iterator& first, const Iterator& last, Function&& f, const unsigned int nthreads = 1)
+  {
+    MESSAGE("Start parallel computation of Hexa with c++ threads...");
+
+    assert(nthreads > 0);
+
+    const unsigned int numTasksTotal = last - first;
+    std::vector<std::thread> threads;
+    Iterator it = first;
+
+    MESSAGE("Number of elements to compute: " << numTasksTotal << "; num of threads: " << nthreads);
+
+    // Distribute tasks among threads
+    if (numTasksTotal <= nthreads)
+    {
+      // A simple case - just one task executed in one thread.
+      // TODO: check if it's faster to do it sequentially
+      threads.reserve(numTasksTotal);
+      for (; it < last; ++it)
+      {
+        threads.emplace_back(f, std::ref(*it));
+      }
+
+      // This line for debug in sequential mode
+      // std::for_each(it, last, f);
+    }
+    else
+    {
+      // Calculate how to distribute elements among threads evenly
+      const unsigned int numTasksInThread = numTasksTotal / nthreads;
+      MESSAGE("Number of tasks in thread: " << numTasksInThread);
+
+      // Store the numbers of tasks per thread
+      std::vector<unsigned int> distTasksInThreads(nthreads, numTasksInThread);
+
+      // Distribute a remainder among saved numbers
+      const unsigned int remainder = numTasksTotal % nthreads;
+      MESSAGE("Remainder of tasks " << remainder << " will be evenly distributed among threads");
+      for (unsigned int i = 0; i < remainder; ++i)
+      {
+        ++distTasksInThreads[i];
+      }
+
+      // Create threads for each number of tasks
+      threads.reserve(nthreads);
+      for (const auto i : distTasksInThreads)
+      {
+        Iterator curLast = it + i;
+
+        // Pass iterators by value and the function by reference!
+        auto lambda = [=,&f](){ std::for_each(it, curLast, f); };
+
+        // Create a thread
+        threads.emplace_back(lambda);
+
+        // Advance iterator to the next step
+        it = curLast;
+      }
+    }
+
+    std::for_each(threads.begin(), threads.end(), [](std::thread& x){ x.join(); });
+
+    MESSAGE("Parallel computation was finished successfully");
+  }
+
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Intersector of TopoDS_Face with all GridLine's
+   */
+  struct FaceGridIntersector
+  {
+    TopoDS_Face _face;
+    TGeomID     _faceID;
+    Grid*       _grid;
+    Bnd_Box     _bndBox;
+    IntCurvesFace_Intersector* _surfaceInt;
+    std::vector< std::pair< GridLine*, F_IntersectPoint > > _intersections;
+
+    FaceGridIntersector(): _grid(0), _surfaceInt(0) {}
+    void Intersect();
+
+    void StoreIntersections()
+    {
+      for ( size_t i = 0; i < _intersections.size(); ++i )
+      {
+        std::multiset< F_IntersectPoint >::iterator ip =
+          _intersections[i].first->_intPoints.insert( _intersections[i].second );
+        ip->_faceIDs.reserve( 1 );
+        ip->_faceIDs.push_back( _faceID );
+      }
+    }
+    const Bnd_Box& GetFaceBndBox()
+    {
+      GetCurveFaceIntersector();
+      return _bndBox;
+    }
+    IntCurvesFace_Intersector* GetCurveFaceIntersector()
+    {
+      if ( !_surfaceInt )
+      {
+        _surfaceInt = new IntCurvesFace_Intersector( _face, Precision::PConfusion() );
+        _bndBox     = _surfaceInt->Bounding();
+        if ( _bndBox.IsVoid() )
+          BRepBndLib::Add (_face, _bndBox);
+      }
+      return _surfaceInt;
+    }
+#ifdef WITH_TBB    
+    bool IsThreadSafe(std::set< const Standard_Transient* >& noSafeTShapes) const;
+#endif
+  };
+  
+#ifdef WITH_TBB
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Structure intersecting certain nb of faces with GridLine's in one thread
+   */
+  struct ParallelIntersector
+  {
+    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 )
+        _faceVec[i].Intersect();
+    }
+  };
+#endif
+
+  template<typename Type>
+  void computeGridIntersection( Type faceGridIntersector )
+  {
+    faceGridIntersector.Intersect();
+  }
+
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Intersector of a surface with a GridLine
+   */
+  struct FaceLineIntersector
+  {
+    double      _tol;
+    double      _u, _v, _w; // params on the face and the line
+    Transition  _transition; // transition at intersection (see IntCurveSurface.cdl)
+    Transition  _transIn, _transOut; // IN and OUT transitions depending of face orientation
+
+    gp_Pln      _plane;
+    gp_Cylinder _cylinder;
+    gp_Cone     _cone;
+    gp_Sphere   _sphere;
+    gp_Torus    _torus;
+    IntCurvesFace_Intersector* _surfaceInt;
+
+    std::vector< F_IntersectPoint > _intPoints;
+
+    void IntersectWithPlane   (const GridLine& gridLine);
+    void IntersectWithCylinder(const GridLine& gridLine);
+    void IntersectWithCone    (const GridLine& gridLine);
+    void IntersectWithSphere  (const GridLine& gridLine);
+    void IntersectWithTorus   (const GridLine& gridLine);
+    void IntersectWithSurface (const GridLine& gridLine);
+
+    /*
+     * 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 )
+    {
+      return -_tol < _w && _w < linLength + _tol;
+    }
+    FaceLineIntersector():_surfaceInt(0) {}
+    ~FaceLineIntersector() { if (_surfaceInt ) delete _surfaceInt; _surfaceInt = 0; }
+  };
+
+} // end namespace Cartesian3D
+} // end namespace StdMeshers
+
+#endif

src/StdMeshers/StdMeshers_Cartesian_3D_Hexahedron.hxx

@@ -0,0 +1,653 @@
+// Copyright (C) 2016-2024  CEA, EDF
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+//  File   : StdMeshers_Cartesian_3D_Hexahedron.hxx
+//  Module : SMESH
+//  Purpose: Make BodyFitting mesh algorithm more modular and testable
+//
+
+#ifndef _SMESH_Cartesian_3D_HEXAHEDRON_HXX_
+#define _SMESH_Cartesian_3D_HEXAHEDRON_HXX_
+
+// BOOST
+#include <boost/container/flat_map.hpp>
+
+// STD
+#include <utilities.h>
+#include <vector>
+
+// SMESH
+#include "SMESH_StdMeshers.hxx"
+#include "StdMeshers_Cartesian_3D_Grid.hxx"
+
+namespace StdMeshers
+{
+namespace Cartesian3D
+{
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief Return cells sharing a link
+   */
+  struct CellsAroundLink
+  {
+    int    _iDir;
+    int    _dInd[4][3];
+    size_t _nbCells[3];
+    int    _i,_j,_k;
+    StdMeshers::Cartesian3D::Grid*  _grid;
+
+    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( grid )
+    {
+      const int iDirOther[3][2] = {{ 1,2 },{ 0,2 },{ 0,1 }};
+      _dInd[1][ iDirOther[iDir][0] ] = -1;
+      _dInd[2][ iDirOther[iDir][1] ] = -1;
+      _dInd[3][ iDirOther[iDir][0] ] = -1; _dInd[3][ iDirOther[iDir][1] ] = -1;
+    }
+    void Init( int i, int j, int k, int link12 = 0 )
+    {
+      int iL = link12 % 4;
+      _i = i - _dInd[iL][0];
+      _j = j - _dInd[iL][1];
+      _k = k - _dInd[iL][2];
+    }
+    bool GetCell( int iL, int& i, int& j, int& k, int& cellIndex, int& linkIndex )
+    {
+      i =  _i + _dInd[iL][0];
+      j =  _j + _dInd[iL][1];
+      k =  _k + _dInd[iL][2];
+      if ( i < 0 || i >= (int)_nbCells[0] ||
+           j < 0 || j >= (int)_nbCells[1] ||
+           k < 0 || k >= (int)_nbCells[2] )
+        return false;
+      cellIndex = _grid->CellIndex( i,j,k );
+      linkIndex = iL + _iDir * 4;
+      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 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), _boundaryCornerNode(0), _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 )
+      {
+        return polygon ? ( _usedInFace == polygon ) : bool( _usedInFace );
+      }
+      StdMeshers::Cartesian3D::TGeomID IsLinked( const StdMeshers::Cartesian3D::B_IntersectPoint* other,
+                                                 StdMeshers::Cartesian3D::TGeomID avoidFace=-1 ) const // returns id of a common face
+      {
+        return _intPoint ? _intPoint->HasCommonFace( other, avoidFace ) : 0;
+      }
+      bool IsOnFace( StdMeshers::Cartesian3D::TGeomID faceID ) const // returns true if faceID is found
+      {
+        return _intPoint ? _intPoint->IsOnFace( faceID ) : false;
+      }
+      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 );
+      void clear();
+
+      friend std::ostream& operator<<(std::ostream& os, const _Node& node)
+      {
+        if (node._node)
+        {
+          os << "Node at hexahedron corner: ";
+          node._node->Print(os); 
+        }
+        else if (node._intPoint && node._intPoint->_node)
+        {
+          os << "Node at intersection point: ";
+          node._intPoint->_node->Print(os); // intersection point
+        }
+        else
+          os << "mesh node is null\n";
+
+        return os;
+      }
+    };
+
+    // --------------------------------------------------------------------------------
+    struct _Link // link connecting two _Node's
+    {
+      static const std::size_t nodesNum = 2;
+
+      _Node* _nodes[nodesNum];
+      _Face* _faces[nodesNum]; // 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(): _nodes{ 0, 0 }, _faces{ 0, 0 } {}
+
+      void clear();
+
+      friend std::ostream& operator<<(std::ostream& os, const _Link& link)
+      {
+        os << "Link:\n";
+
+        for (std::size_t i = 0; i < nodesNum; ++i)
+        {
+          if (link._nodes[i])
+            os << *link._nodes[i];
+          else
+            os << "link node with index " << i << " is null\n";
+        }
+
+        os << "_fIntPoints: " << link._fIntPoints.size() << '\n';
+        os << "_fIntNodes: " << link._fIntNodes.size() << '\n';
+        os << "_splits: " << link._splits.size() << '\n';
+
+        return os;
+      }
+    };
+
+    // --------------------------------------------------------------------------------
+    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
+      {
+        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; }
+
+      // returns supporting FACEs
+      std::vector< StdMeshers::Cartesian3D::TGeomID > GetNotUsedFaces
+      (const std::set<StdMeshers::Cartesian3D::TGeomID>& usedIDs ) const
+      {
+        std::vector< StdMeshers::Cartesian3D::TGeomID > faces;
+        const StdMeshers::Cartesian3D::B_IntersectPoint *ip0, *ip1;
+        if (( ip0 = _link->_nodes[0]->_intPoint ) &&
+            ( ip1 = _link->_nodes[1]->_intPoint ))
+        {
+          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;
+        }
+      }
+      void RemoveFace( const _Face* f )
+      {
+        if ( !_link->_faces[0] ) 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;
+          }
+        }
+      }
+
+      friend std::ostream& operator<<(std::ostream& os, const _OrientedLink& link)
+      {
+        if (link._link)
+          os << "Oriented " << *link._link;
+        else
+          os << "Oriented link is null\n";
+
+        return os;
+      }
+    };
+
+    // --------------------------------------------------------------------------------
+    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
+
+        _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;
+      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
+
+      _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;
+            }
+        }
+        _Link l;
+        l._nodes[0] = n0;
+        l._nodes[1] = n1;
+        _polyLinks.push_back( l );
+        _links.push_back( _OrientedLink( &_polyLinks.back() ));
+      }
+
+      friend std::ostream& operator<<(std::ostream& os, const _Face& face)
+      {
+        os << "Face " << face._name << '\n';
+
+        os << "Links on GridLines: \n";
+        for (const auto& link : face._links)
+        {
+          os << link;
+        }
+
+        os << "Links added to close a polygonal face: \n";
+        for (const auto& link : face._polyLinks)
+        {
+          os << link;
+        }
+
+        os << "Nodes at intersection with EDGEs: \n";
+        for (const auto node : face._eIntNodes)
+        {
+          if (node)
+          {
+            os << *node;
+          }
+        }
+
+        return os;
+      }
+    };
+
+    // --------------------------------------------------------------------------------
+    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(); }
+
+        friend std::ostream& operator<<(std::ostream& os, const _nodeDef& node)
+        {
+          if (node._node)
+          {
+            os << "Node at hexahedron corner: ";
+            node._node->Print(os); 
+          }
+          else if (node._intPoint && node._intPoint->_node)
+          {
+            os << "Node at intersection point: ";
+            node._intPoint->_node->Print(os); // intersection point
+          }
+          else
+            os << "mesh node is null\n";
+
+          return os;
+        }
+      };
+
+      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() )); }
+
+      std::vector<std::set<std::pair<int, int>>> getPolygonsEdges() const;
+      std::vector<std::set<std::pair<int, int>>> findOpenEdges() const;
+      int getStartNodeIndex(const int polygon) const;
+      std::map<int, std::vector<int>> findOverlappingPolygons() const;
+      bool divideOverlappingPolygons();
+      bool fixOpenEdgesPolygons();
+      bool capOpenEdgesPolygons(const std::vector<std::set<std::pair<int, int>>>& edgesByPolygon);
+      bool removeOpenEdgesPolygons(const std::vector<std::set<std::pair<int, int>>>& edgesByPolygon);
+
+      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;
+        }
+
+        friend std::ostream& operator<<(std::ostream& os, const _linkDef& link)
+        {
+          os << "Link def:\n";
+
+          os << link._node1;
+          if (link.first)
+            os << "first: " << link.first;
+
+          if (link.second)
+            os << "second: " << link.second;
+
+          os << "_loopIndex: " << link._loopIndex << '\n';
+
+          return os;
+        }
+      };
+    };
+
+    // topology of a hexahedron
+    static const std::size_t HEX_NODES_NUM = 8;
+    static const std::size_t HEX_LINKS_NUM = 12;
+    static const std::size_t HEX_QUADS_NUM = 6;
+    _Node _hexNodes [HEX_NODES_NUM];
+    _Link _hexLinks [HEX_LINKS_NUM];
+    _Face _hexQuads [HEX_QUADS_NUM];
+
+    // 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 clearNodesLinkedToNull(const Solid* solid, SMESH_MesherHelper& helper);
+    bool isSplittedLink(const Solid* solid, SMESH_MesherHelper& helper, const Hexahedron::_Link& linkIn) const;
+    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);
+    size_t findCoplanarPolygon
+           (const _Face& thePolygon,
+            const size_t nbQuadPolygons,
+            std::vector< _OrientedLink* >& freeLinks,
+            int& nbFreeLinks,
+            const E_IntersectPoint& ipTmp,
+            std::set< StdMeshers::Cartesian3D::TGeomID >& usedFaceIDs,
+            std::map< StdMeshers::Cartesian3D::TGeomID, std::vector< const B_IntersectPoint* > >& tmpAddedFace,
+            const StdMeshers::Cartesian3D::TGeomID& curFace);
+    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

test/CTestTestfileInstall.cmake

@@ -41,3 +41,10 @@ FOREACH(tfile ${CPP_TESTS})
   ADD_TEST(${TEST_NAME} ${BASE_NAME} )
   SET_TESTS_PROPERTIES(${TEST_NAME} PROPERTIES LABELS "${COMPONENT_NAME};${COMPONENT_NAME}_tests")
 ENDFOREACH()
+
+FOREACH(tfile ${UNIT_TESTS})
+  GET_FILENAME_COMPONENT(BASE_NAME ${tfile} NAME_WE)
+  SET(TEST_NAME SMESH_${BASE_NAME})
+  ADD_TEST(${TEST_NAME} ${BASE_NAME} )
+  SET_TESTS_PROPERTIES(${TEST_NAME} PROPERTIES LABELS "${COMPONENT_NAME};${COMPONENT_NAME}_tests")
+ENDFOREACH()

test/data/HexahedronTest/NRTM1.brep

@@ -0,0 +1,441 @@
+DBRep_DrawableShape
+
+CASCADE Topology V1, (c) Matra-Datavision
+Locations 0
+Curve2ds 36
+1 0 0 1 0 
+1 0 0 1 0 
+1 10 0 0 -1 
+1 0 0 0 1 
+1 0 -10 1 0 
+1 0 0 1 0 
+1 0 0 0 -1 
+1 0 0 0 1 
+1 0 0 0 1 
+1 0 0 1 0 
+1 0 10 1 0 
+1 0 0 1 0 
+1 10 0 0 1 
+1 0 0 1 0 
+1 0 10 1 0 
+1 10 0 0 1 
+1 10 0 0 1 
+1 10 0 0 -1 
+1 0 0 0 1 
+1 0 10 1 0 
+1 0 10 1 0 
+1 0 -10 1 0 
+1 10 0 0 1 
+1 0 0 0 -1 
+2 5 5 1 0 -0 1 1
+1 0 6 1 0 
+2 5 5 1 0 -0 1 1
+1 0 6 1 0 
+1 6.2831853071795862 -0 0 1 
+1 0 -0 0 1 
+1 0 9 1 0 
+2 0 0 1 0 -0 1 1
+1 6.2831853071795862 -0 0 1 
+1 0 -0 0 1 
+1 0 0 1 0 
+2 0 0 1 0 -0 1 1
+Curves 17
+1 -5 -5 0 0 0 1 
+1 -5 -5 10 -0 1 0 
+1 -5 5 0 0 0 1 
+1 -5 -5 0 -0 1 0 
+1 -5 -5 0 1 0 -0 
+1 5 -5 0 0 0 1 
+1 -5 -5 10 1 0 -0 
+1 -5 5 10 1 0 -0 
+1 5 -5 10 -0 1 0 
+1 -5 5 0 1 0 -0 
+1 5 5 0 0 0 1 
+1 5 -5 0 -0 1 0 
+2 0 0 0 0 0 1 1 0 -0 -0 1 0 1
+1 1 -2.4492935982947064e-16 -6 0 0 1 
+2 0 0 3 0 0 1 1 0 -0 -0 1 0 1
+1 1 -2.4492935982947064e-16 -6 0 0 1 
+2 0 0 -6 0 0 1 1 0 -0 -0 1 0 1
+Polygon3D 0
+PolygonOnTriangulations 0
+Surfaces 11
+1 -5 -5 0 1 0 -0 0 0 1 0 -1 0 
+1 -5 -5 0 -0 1 0 0 0 1 1 0 -0 
+1 -5 -5 10 0 0 1 1 0 -0 -0 1 0 
+1 -5 5 0 -0 1 0 0 0 1 1 0 -0 
+1 -5 -5 0 0 0 1 1 0 -0 -0 1 0 
+1 5 -5 0 1 0 -0 0 0 1 0 -1 0 
+2 0 0 -6 0 0 1 1 0 -0 -0 1 0 1
+1 -5 -5 0 0 0 1 1 0 -0 -0 1 0 
+2 0 0 -6 0 0 1 1 0 -0 -0 1 0 1
+1 0 0 3 0 0 1 1 0 -0 -0 1 0 
+1 0 0 -6 0 0 1 1 0 -0 -0 1 0 
+Triangulations 0
+
+TShapes 58
+Ve
+1e-07
+-5 -5 10
+0 0
+
+0101101
+*
+Ve
+1e-07
+-5 -5 0
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  1 0 0 10
+2  1 1 0 0 10
+2  2 2 0 0 10
+0
+
+0101000
+-58 0 +57 0 *
+Ve
+1e-07
+-5 5 10
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  2 0 0 10
+2  3 1 0 0 10
+2  4 3 0 0 10
+0
+
+0101000
+-55 0 +58 0 *
+Ve
+1e-07
+-5 5 0
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  3 0 0 10
+2  5 1 0 0 10
+2  6 4 0 0 10
+0
+
+0101000
+-55 0 +53 0 *
+Ed
+ 1e-07 1 1 0
+1  4 0 0 10
+2  7 1 0 0 10
+2  8 5 0 0 10
+0
+
+0101000
+-53 0 +57 0 *
+Wi
+
+0101100
+-56 0 -54 0 +52 0 +51 0 *
+Fa
+0  1e-07 1 0
+
+0101000
++50 0 *
+Ve
+1e-07
+5 -5 0
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  5 0 0 10
+2  9 2 0 0 10
+2  10 5 0 0 10
+0
+
+0101000
+-48 0 +57 0 *
+Ve
+1e-07
+5 -5 10
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  6 0 0 10
+2  11 2 0 0 10
+2  12 6 0 0 10
+0
+
+0101000
+-46 0 +48 0 *
+Ed
+ 1e-07 1 1 0
+1  7 0 0 10
+2  13 2 0 0 10
+2  14 3 0 0 10
+0
+
+0101000
+-46 0 +58 0 *
+Wi
+
+0101100
+-47 0 -45 0 +44 0 +56 0 *
+Fa
+0  1e-07 2 0
+
+0101000
++43 0 *
+Ve
+1e-07
+5 5 10
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  8 0 0 10
+2  15 3 0 0 10
+2  16 4 0 0 10
+0
+
+0101000
+-41 0 +55 0 *
+Ed
+ 1e-07 1 1 0
+1  9 0 0 10
+2  17 3 0 0 10
+2  18 6 0 0 10
+0
+
+0101000
+-41 0 +46 0 *
+Wi
+
+0101100
+-54 0 -40 0 +39 0 +44 0 *
+Fa
+0  1e-07 3 0
+
+0101000
++38 0 *
+Ve
+1e-07
+5 5 0
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  10 0 0 10
+2  19 4 0 0 10
+2  20 5 0 0 10
+0
+
+0101000
+-36 0 +53 0 *
+Ed
+ 1e-07 1 1 0
+1  11 0 0 10
+2  21 4 0 0 10
+2  22 6 0 0 10
+0
+
+0101000
+-41 0 +36 0 *
+Wi
+
+0101100
+-35 0 -34 0 +40 0 +52 0 *
+Fa
+0  1e-07 4 0
+
+0101000
++33 0 *
+Ed
+ 1e-07 1 1 0
+1  12 0 0 10
+2  23 5 0 0 10
+2  24 6 0 0 10
+0
+
+0101000
+-36 0 +48 0 *
+Wi
+
+0101100
+-51 0 -35 0 +47 0 +31 0 *
+Ve
+1.00000000244929e-07
+1 -2.44929359829471e-16 0
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  13 0 0 6.28318530717959
+2  25 5 0 0 6.28318530717959
+2  26 7 0 0 6.28318530717959
+2  27 8 0 0 6.28318530717959
+2  28 9 0 0 6.28318530717959
+0
+
+0101000
++29 0 -29 0 *
+Wi
+
+0101100
+-28 0 *
+Fa
+0  1e-07 5 0
+
+0101000
++30 0 +27 0 *
+Wi
+
+0101100
+-45 0 -39 0 +34 0 +31 0 *
+Fa
+0  1e-07 6 0
+
+0101000
++25 0 *
+Ve
+1e-07
+1 -2.44929359829471e-16 3
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  14 0 6 9
+3  29 30CN 7 0 6 9
+0
+
+0101000
++29 0 -23 0 *
+Ed
+ 1e-07 1 1 0
+1  15 0 0 6.28318530717959
+2  31 7 0 0 6.28318530717959
+2  32 10 0 0 6.28318530717959
+0
+
+0101000
++23 0 -23 0 *
+Wi
+
+0101100
++22 0 +28 0 -22 0 -21 0 *
+Fa
+0  1e-07 7 0
+
+0101000
++20 0 *
+Wi
+
+0101100
++21 0 *
+Fa
+0  1e-07 10 0
+
+0101000
++18 0 *
+Sh
+
+0101100
+-49 0 -42 0 +37 0 +32 0 -26 0 +24 0 -19 0 -17 0 *
+So
+
+0100000
++16 0 *
+Wi
+
+0101100
++28 0 *
+Fa
+0  1e-07 8 0
+
+0101000
++14 0 *
+Sh
+
+0101100
+-13 0 +19 0 +17 0 *
+So
+
+0100000
++12 0 *
+Ve
+1e-07
+1 -2.44929359829471e-16 -6
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  16 0 0 6
+3  33 34CN 9 0 0 6
+0
+
+0101000
++10 0 -29 0 *
+Ed
+ 1e-07 1 1 0
+1  17 0 0 6.28318530717959
+2  35 9 0 0 6.28318530717959
+2  36 11 0 0 6.28318530717959
+0
+
+0101000
++10 0 -10 0 *
+Wi
+
+0101100
+-28 0 +9 0 +8 0 -9 0 *
+Fa
+0  1e-07 9 0
+
+0101000
++7 0 *
+Wi
+
+0101100
+-8 0 *
+Fa
+0  1e-07 11 0
+
+0101000
+-5 0 *
+Sh
+
+0101100
++6 0 +13 0 -4 0 *
+So
+
+0100000
++3 0 *
+Co
+
+1100000
++15 0 +11 0 +2 0 *
+
++1 0 
+0
+

test/data/HexahedronTest/NRTMJ4.brep

@@ -0,0 +1,598 @@
+DBRep_DrawableShape
+
+CASCADE Topology V1, (c) Matra-Datavision
+Locations 0
+Curve2ds 54
+1 0 0 1 0 
+1 0 0 1 0 
+1 100 0 0 -1 
+1 0 0 0 1 
+1 0 -50 1 0 
+1 0 0 1 0 
+1 0 0 0 -1 
+1 0 0 0 1 
+1 0 60 1 0 
+1 0 0 1 0 
+1 16.899999999999999 50 0 1 
+1 0 0 1 0 
+1 16.899999999999999 50 1 0 
+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 0 50 1 0 
+1 100 0 0 1 
+1 60 0 0 1 
+1 100 0 0 -1 
+1 0 60 1 0 
+1 0 -50 1 0 
+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 60 0 0 1 
+1 0 0 0 -1 
+2 25 25 -1 0 0 1 10
+1 0 0 1 0 
+1 16.899999999999999 0 0 -1 
+1 10 0 0 1 
+1 0 0 0 1 
+1 0 0 0 -1 
+1 59.899999999999999 0 0 -1 
+1 0 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 27
+1 0 0 0 0 0 1 
+1 0 0 100 -0 1 0 
+1 0 50 0 0 0 1 
+1 0 0 0 -0 1 0 
+1 60 0 0 0 0 1 
+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 
+1 60 0 100 -0 1 0 
+1 60 50 0 0 0 1 
+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
+1 60 0 16.899999999999999 -0 1 0 
+1 50 0 16.899999999999999 -0 1 0 
+1 50 0 76.799999999999997 -0 1 0 
+1 60 0 76.799999999999997 -0 1 0 
+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 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 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 60 0 0 1 0 -0 0 0 1 0 -1 0 
+2 25 25 0 0 0 -1 -1 0 -0 0 1 0 10
+1 25 25 -50 0 0 -1 -1 0 -0 0 1 0 
+Triangulations 0
+
+TShapes 72
+Ve
+1e-07
+0 0 100
+0 0
+
+0101101
+*
+Ve
+1e-07
+0 0 0
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  1 0 0 100
+2  1 1 0 0 100
+2  2 2 0 0 100
+0
+
+0101000
+-72 0 +71 0 *
+Ve
+1e-07
+0 50 100
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  2 0 0 50
+2  3 1 0 0 50
+2  4 3 0 0 50
+0
+
+0101000
+-69 0 +72 0 *
+Ve
+1e-07
+0 50 0
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  3 0 0 100
+2  5 1 0 0 100
+2  6 4 0 0 100
+0
+
+0101000
+-69 0 +67 0 *
+Ed
+ 1e-07 1 1 0
+1  4 0 0 50
+2  7 1 0 0 50
+2  8 5 0 0 50
+0
+
+0101000
+-67 0 +71 0 *
+Wi
+
+0101100
+-70 0 -68 0 +66 0 +65 0 *
+Fa
+0  1e-07 1 0
+
+0101000
++64 0 *
+Ve
+1e-07
+60 0 0
+0 0
+
+0101101
+*
+Ve
+1e-07
+60 0 16.9
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  5 0 0 16.9
+2  9 2 0 0 16.9
+2  10 6 0 0 16.9
+0
+
+0101000
++62 0 -61 0 *
+Ve
+1e-07
+50 0 16.9
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  6 0 0 10
+2  11 2 0 0 10
+2  12 7 0 0 10
+0
+
+0101000
+-61 0 +59 0 *
+Ve
+1e-07
+50 0 76.8
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  7 0 0 59.9
+2  13 2 0 0 59.9
+2  14 8 0 0 59.9
+0
+
+0101000
+-57 0 +59 0 *
+Ve
+1e-07
+60 0 76.8
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  8 0 0 10
+2  15 2 0 0 10
+2  16 9 0 0 10
+0
+
+0101000
+-55 0 +57 0 *
+Ve
+1e-07
+60 0 100
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 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  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  11 0 0 60
+2  21 2 0 0 60
+2  22 5 0 0 60
+0
+
+0101000
+-62 0 +71 0 *
+Wi
+
+0101100
+-60 0 +58 0 -56 0 -54 0 -52 0 +51 0 +70 0 -50 0 *
+Fa
+0  1e-07 2 0
+
+0101000
++49 0 *
+Ve
+1e-07
+60 50 100
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  12 0 0 60
+2  23 3 0 0 60
+2  24 4 0 0 60
+0
+
+0101000
+-47 0 +69 0 *
+Ed
+ 1e-07 1 1 0
+1  13 0 0 50
+2  25 3 0 0 50
+2  26 10 0 0 50
+0
+
+0101000
+-47 0 +53 0 *
+Wi
+
+0101100
+-68 0 -46 0 +45 0 +51 0 *
+Fa
+0  1e-07 3 0
+
+0101000
++44 0 *
+Ve
+1e-07
+60 50 0
+0 0
+
+0101101
+*
+Ve
+1e-07
+60 50 16.9
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  14 0 0 16.9
+2  27 4 0 0 16.9
+2  28 6 0 0 16.9
+0
+
+0101000
++42 0 -41 0 *
+Ve
+1e-07
+50 50 16.9
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  15 0 0 10
+2  29 4 0 0 10
+2  30 7 0 0 10
+0
+
+0101000
+-41 0 +39 0 *
+Ve
+1e-07
+50 50 76.8
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  16 0 0 59.9
+2  31 4 0 0 59.9
+2  32 8 0 0 59.9
+0
+
+0101000
+-37 0 +39 0 *
+Ve
+1e-07
+60 50 76.8
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 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  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  19 0 0 60
+2  37 4 0 0 60
+2  38 5 0 0 60
+0
+
+0101000
+-42 0 +67 0 *
+Wi
+
+0101100
+-40 0 +38 0 -36 0 -34 0 -33 0 +46 0 +66 0 -32 0 *
+Fa
+0  1e-07 4 0
+
+0101000
++31 0 *
+Ed
+ 1e-07 1 1 0
+1  20 0 0 50
+2  39 5 0 0 50
+2  40 6 0 0 50
+0
+
+0101000
+-42 0 +62 0 *
+Wi
+
+0101100
+-65 0 -32 0 +50 0 +29 0 *
+Ve
+1e-07
+15 25 0
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1  21 0 0 6.28318530717959
+2  41 5 0 0 6.28318530717959
+2  42 11 0 0 6.28318530717959
+0
+
+0101000
++27 0 -27 0 *
+Wi
+
+0101100
++26 0 *
+Fa
+0  1e-07 5 0
+
+0101000
++28 0 +25 0 *
+Ed
+ 1e-07 1 1 0
+1  22 0 0 50
+2  43 6 0 0 50
+2  44 7 0 0 50
+0
+
+0101000
+-41 0 +61 0 *
+Wi
+
+0101100
+-60 0 -23 0 +40 0 +29 0 *
+Fa
+0  1e-07 6 0
+
+0101000
++22 0 *
+Ed
+ 1e-07 1 1 0
+1  23 0 0 50
+2  45 7 0 0 50
+2  46 8 0 0 50
+0
+
+0101000
+-39 0 +59 0 *
+Wi
+
+0101100
+-20 0 -38 0 +23 0 +58 0 *
+Fa
+0  1e-07 7 0
+
+0101000
++19 0 *
+Ed
+ 1e-07 1 1 0
+1  24 0 0 50
+2  47 8 0 0 50
+2  48 9 0 0 50
+0
+
+0101000
+-37 0 +57 0 *
+Wi
+
+0101100
+-56 0 -17 0 +36 0 +20 0 *
+Fa
+0  1e-07 8 0
+
+0101000
++16 0 *
+Ed
+ 1e-07 1 1 0
+1  25 0 0 50
+2  49 9 0 0 50
+2  50 10 0 0 50
+0
+
+0101000
+-35 0 +55 0 *
+Wi
+
+0101100
+-17 0 -34 0 +14 0 +54 0 *
+Fa
+0  1e-07 9 0
+
+0101000
++13 0 *
+Wi
+
+0101100
+-52 0 -45 0 +33 0 +14 0 *
+Fa
+0  1e-07 10 0
+
+0101000
++11 0 *
+Ve
+1e-07
+15 25 -50
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 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  27 0 0 50
+3  53 54CN 11 0 0 50
+0
+
+0101000
+-9 0 +27 0 *
+Wi
+
+0101100
+-8 0 +7 0 +26 0 -7 0 *
+Fa
+0  1e-07 11 0
+
+0101000
++6 0 *
+Wi
+
+0101100
++8 0 *
+Fa
+0  1e-07 12 0
+
+0101000
++4 0 *
+Sh
+
+0101100
+-63 0 -48 0 +43 0 +30 0 -24 0 +21 0 +18 0 +15 0 -12 0 +10 0 
++5 0 +3 0 *
+So
+
+1100000
++2 0 *
+
++1 0 
+0
+

test/tests.set

@@ -127,6 +127,11 @@ SET(CPP_TESTS
   SMESH_RegularGridTest
 )
 
+SET(UNIT_TESTS # Any unit test add in src names space should be added here 
+  HexahedronTest
+  HexahedronCanonicalShapesTest
+)
+
 # The following tests can be executed without driver, just by python.
 # ----------------------------------------------------------------------------