Merge pull request #12680 from KratosMultiphysics/core/obj--improvement

[Core] Extend `ObjIO` with new functionality

kratos/input_output/obj_io.h

@@ -130,6 +130,29 @@ class KRATOS_API(KRATOS_CORE) ObjIO
      */
     void WriteModelPart(const ModelPart& rThisModelPart) override;
 
+    /**
+     * @brief This method clean up the problematic geometries (null area geometries) in the mesh.
+     * @details For the moment it will be assumed that the problematic geometries are triangles with two nodes in the same position. This means several this:
+     * - The triangle is degenerated, one of the nodes is redundant and should be cleaned up.
+     * - The triangle must be removed as well. The sides must be reconnected, so the mesh is water-tight.
+     * - The nodes and entities must be renumbered
+     * @param rThisModelPart Reference to the model part to clean up into.
+     * @param rEntityType The entity type to create in the model part. Can be "element", "condition" or "geometry".
+     * @param FirstNodeId The first node ID to create in the model part.
+     * @param FirstELementId The first element ID to create in the model part.
+     * @param FirstConditionId The first condition ID to create in the model part.
+     * @param AreaTolerance The tolerance to consider a geometry as problematic.
+     * @todo This could be moved to an independent utility. It is a temporary solution to clean up problematic geometries in the mesh that can be found in the OBJ files
+     */
+    static void CleanUpProblematicGeometriesInMesh(
+        ModelPart& rThisModelPart,
+        const std::string& rEntityType = "element",
+        const IndexType FirstNodeId = 1,
+        const IndexType FirstElementId = 1,
+        const IndexType FirstConditionId = 1,
+        const double AreaTolerance = 1.0e-6
+        );
+
     ///@}
     ///@name Input and output
     ///@{
@@ -148,13 +171,18 @@ class KRATOS_API(KRATOS_CORE) ObjIO
     ///@name Protected member Variables
     ///@{
 
-    Parameters mParameters;         /// The configuration parameters
+    Parameters mParameters;          /// The configuration parameters
 
-    IndexType mFirstNodeId = 0;     /// The first node ID
-    IndexType mNextNodeId = 0;      /// The next node ID
-    IndexType mNextElementId = 0;   /// The next element ID
-    IndexType mNextConditionId = 0; /// The next condition ID
-    IndexType mNormalCounter = 0;   /// The normal counter
+    IndexType mFirstNodeId = 0;      /// The first node ID
+    IndexType mNextNodeId = 0;       /// The next node ID
+
+    IndexType mFirstElementId = 0;   /// The first element ID
+    IndexType mNextElementId = 0;    /// The next element ID
+
+    IndexType mFirstConditionId = 0; /// The first condition ID
+    IndexType mNextConditionId = 0;  /// The next condition ID
+
+    IndexType mNormalCounter = 0;    /// The normal counter
 
     ///@}
 private:
@@ -174,11 +202,13 @@ class KRATOS_API(KRATOS_CORE) ObjIO
      * @param rThisModelPart Reference to the model part to read into.
      * @param rEntityType The entity type to create in the model part. Can be "element" or "condition".
      * @param NormalAsHistoricalVariable If true, the normals are stored as historical variables in the nodes.
+     * @param DecomposeQuadrilateral If true, the quadrilateral faces are decomposed into triangles.
      */
     void ReadVerticesAndFaces(
         ModelPart& rThisModelPart,
         const std::string& rEntityType = "element",
-        const bool NormalAsHistoricalVariable = false
+        const bool NormalAsHistoricalVariable = false,
+        const bool DecomposeQuadrilateral = false
         );
 
     /**
@@ -211,11 +241,13 @@ class KRATOS_API(KRATOS_CORE) ObjIO
      * @param rThisModelPart Reference to the model part.
      * @param rLine The line containing the face data.
      * @param rEntityType The entity type to create in the model part. Can be "element" or "condition".
+     * @param DecomposeQuadrilateral If true, the quadrilateral faces are decomposed into triangles.
      */
     void ParseFaceLine(
         ModelPart& rThisModelPart,
         const std::string& rLine,
-        const std::string& rEntityType = "element"
+        const std::string& rEntityType = "element",
+        const bool DecomposeQuadrilateral = false
         );
 
     /**
@@ -225,12 +257,50 @@ class KRATOS_API(KRATOS_CORE) ObjIO
      */
     std::vector<std::string> Tokenize(const std::string& rLine);
 
+    /**
+     * @brief Clean up the problematic geometries (null area geometries) in the mesh.
+     * @details For the moment it will be assumed that the problematic geometries are triangles with two nodes in the same position. This means several this:
+     * - The triangle is degenerated, one of the nodes is redundant and should be cleaned up.
+     * - The triangle must be removed as well. The sides must be reconnected, so the mesh is water-tight.
+     * - The nodes and entities must be renumbered
+     * @param rThisModelPart Reference to the model part to clean up into.
+     * @param rEntityType The entity type to create in the model part. Can be "element", "condition" or "geometry".
+     * @param FirstNodeId The first node ID to create in the model part.
+     * @param FirstEntityId The first entity ID to create in the model part.
+     * @param AreaTolerance The tolerance to consider a geometry as problematic.
+     * @tparam TEntityType The entity type to create in the model part. Can be Element, Condition or Geometry.
+     */
+    template <typename TEntityType>
+    static void CleanUpProblematicGeometries(
+        ModelPart& rThisModelPart,
+        const IndexType FirstNodeId = 1,
+        const IndexType FirstEntityId = 1,
+        const double AreaTolerance = 1.0e-6
+        );
+
+    /// Helper function declarations
+
+    /**
+     * @brief Get the nodes of the model part.
+     * @param rThisModelPart Reference to the model part.
+     * @return The nodes of the model part.
+     * @tparam TEntityType The entity type to get from the model part. Can be Element,  Condition or Geometry.
+     */
+    template <typename TEntityType>
+    static void RemoveEntitiesAndNodes(ModelPart& rThisModelPart);
+
     ///@}
 
 }; // Class ObjIO
 
-///@}
+// Helper function definitions specialization
+template <>
+void ObjIO::RemoveEntitiesAndNodes<Element>(ModelPart& rThisModelPart);
+
+template <>
+void ObjIO::RemoveEntitiesAndNodes<Condition>(ModelPart& rThisModelPart);
 
+///@}
 ///@name Input and output
 ///@{
 
@@ -248,4 +318,4 @@ inline std::ostream& operator << (std::ostream& rOStream,
 
 ///@} addtogroup block
 
-}  // namespace Kratos.
+}  // namespace Kratos.

kratos/python/add_io_to_python.cpp

@@ -186,15 +186,19 @@ void  AddIOToPython(pybind11::module& m)
     py::class_<StlIO, StlIO::Pointer, IO>(m, "StlIO")
         .def(py::init<std::filesystem::path const& >())
         .def(py::init<std::filesystem::path const&, Parameters>())
-        .def("ReadModelPart", &StlIO::ReadModelPart)
-        .def("WriteModelPart", &StlIO::WriteModelPart)
         ;
 
     py::class_<ObjIO, ObjIO::Pointer, IO>(m, "ObjIO")
         .def(py::init<std::filesystem::path const& >())
         .def(py::init<std::filesystem::path const&, Parameters>())
-        .def("ReadModelPart", &ObjIO::ReadModelPart)
-        .def("WriteModelPart", &ObjIO::WriteModelPart)
+        .def_static("CleanUpProblematicGeometriesInMesh", &ObjIO::CleanUpProblematicGeometriesInMesh,
+          py::arg("rThisModelPart"),
+          py::arg("rEntityType") = "element",
+          py::arg("FirstNodeId") = 1,
+          py::arg("FirstElementId") = 1,
+          py::arg("FirstConditionId") = 1,
+          py::arg("AreaTolerance") = 1.0e-6,
+          "Clean up the problematic geometries (null area geometries) in the mesh.")
         ;
 
     // Import of CAD models to the model part

kratos/python_scripts/modelers/import_obj_modeler.py

@@ -38,7 +38,7 @@ def __init__(self, model, settings):
         super().__init__(model, settings)
 
         # Cannot validate as settings may differ among input types
-        settings.AddMissingParameters(self.__GetDefaultSettings())
+        settings.RecursivelyAddMissingParameters(self.__GetDefaultSettings())
 
         # Declare required member variables
         self.model_part = None
@@ -102,9 +102,12 @@ def __GetDefaultSettings(cls):
             "input_filename"  : "",
             "model_part_name" : "",
             "obj_io_settings"  : {
-                "open_mode"            : "read",
-                "entity_type"          : "element",
-                "normal_as_historical" : false
+                "open_mode"                       : "read",
+                "entity_type"                     : "element",
+                "decompose_quads_into_triangles"  : false,
+                "normal_as_historical"            : false,
+                "clean_up_problematic_geometries" : true,
+                "area_tolerance"                  : 1.0e-6
             }
         }''')
         return default_settings

kratos/tests/auxiliar_files_for_python_unittest/obj_files/cube.obj

@@ -12,10 +12,10 @@ v -0.5 0.5 -0.5  # 8
 
 # Faces
 f 1 2 3 4    # front
-f 5 6 7 8    # back
-f 1 5 8 4    # left
+f 8 7 6 5    # back
+f 4 8 5 1    # left
 f 2 6 7 3    # right
-f 1 2 6 5    # bottom
+f 5 6 2 1    # bottom
 f 4 3 7 8    # top
 
 # Normal

kratos/tests/auxiliar_files_for_python_unittest/obj_files/cube_degenerated.obj

@@ -0,0 +1,24 @@
+# Simple cube
+
+# Vertices
+v -0.5 -0.5 0.5
+v 0.5 -0.5 0.5
+v 0.5 0.5 0.5
+v -0.5 0.5 0.5
+v -0.5 -0.5 -0.5
+v 0.5 -0.5 -0.5
+v 0.5 0.5 -0.5
+v -0.5 0.5 -0.5
+v -0.5 -0.5 0.5
+v -0.5 0.5 -0.5
+
+# Faces
+f 1 2 3 4
+f 8 7 6 5
+f 4 8 5 1
+f 2 6 7 3
+f 5 6 2 1
+f 4 3 7 8
+f 1 2 9
+f 8 7 10
+

kratos/tests/test_obj_io.py

@@ -24,17 +24,28 @@ def debug_vtk(model_part):
     It sets up the necessary parameters, including the model part name and the nodal data value variables.
     The VTK output process is then executed through its initialization, pre-solution loop, and solution step phases.
     """
+    # Compute normals
+    for cond in model_part.Conditions:
+        geom = cond.GetGeometry()
+        normal = geom.UnitNormal()
+        cond.SetValue(KratosMultiphysics.NORMAL, normal)
+
+    # Generate VTK output
     from KratosMultiphysics.vtk_output_process import VtkOutputProcess
     parameters = KratosMultiphysics.Parameters("""{
-        "model_part_name"            : "",
-        "nodal_data_value_variables" : ["NORMAL"]
+        "model_part_name"                : "",
+        "nodal_data_value_variables"     : ["NORMAL"],
+        "condition_data_value_variables" : ["NORMAL"]
     }
     """)
     parameters["model_part_name"].SetString(model_part.Name)
     output = VtkOutputProcess(model_part.GetModel(), parameters)
     output.ExecuteInitialize()
     output.ExecuteBeforeSolutionLoop()
     output.ExecuteInitializeSolutionStep()
+    output.PrintOutput()
+    output.ExecuteFinalizeSolutionStep()
+    output.ExecuteFinalize()
 
 def calculate_analytical_normal(node):
     """
@@ -111,7 +122,7 @@ def WriteModelPartToOBJ(model_part, obj_file):
     return obj_file
 
 # Define a function to read a Kratos model part from an OBJ file
-def ReadModelPartFromOBJ(model_part, obj_file):
+def ReadModelPartFromOBJ(model_part, obj_file, decompose_quad=False):
     """
     Read a Kratos model part from an OBJ file.
 
@@ -120,7 +131,12 @@ def ReadModelPartFromOBJ(model_part, obj_file):
         data_comm (KratosMultiphysics.DataCommunicator): The data communicator.
         obj_file (str): The name of the OBJ file to read from.
     """
-    read_settings = KratosMultiphysics.Parameters("""{"open_mode" : "read", "entity_type" : "element"}""")
+    read_settings = KratosMultiphysics.Parameters("""{
+        "open_mode"                      : "read",
+        "entity_type"                    : "condition",
+        "decompose_quads_into_triangles" : false
+    }""")
+    read_settings["decompose_quads_into_triangles"].SetBool(decompose_quad)
     obj_io = KratosMultiphysics.ObjIO(obj_file, read_settings)
     obj_io.ReadModelPart(model_part)
 
@@ -161,7 +177,7 @@ def test_ReadObjIO(self):
 
         # Assert that the model part has the correct number of nodes and elements
         self.assertEqual(self.model_part.NumberOfNodes(), 8)
-        self.assertEqual(self.model_part.NumberOfElements(), 6)
+        self.assertEqual(self.model_part.NumberOfConditions(), 6)
 
         # Assert that the normals are the same as the nodes coordinates
         for node in self.model_part.Nodes:
@@ -170,6 +186,53 @@ def test_ReadObjIO(self):
             self.assertAlmostEqual(normal[1], node.Y)
             self.assertAlmostEqual(normal[2], node.Z)
 
+    def test_ReadObjIOTriangles(self):
+        """
+        Test the ReadModelPart function from ObjIO
+        """
+        # Create a model part and set the domain size
+        self.current_model = KratosMultiphysics.Model()
+        self.model_part = self.current_model.CreateModelPart("Main")
+        self.model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] = 3
+
+        # Read a model part from an OBJ file
+        obj_name = GetFilePath("auxiliar_files_for_python_unittest/obj_files/cube.obj")
+        ReadModelPartFromOBJ(self.model_part, obj_name, True)
+
+        # # Debug
+        # debug_vtk(self.model_part)
+
+        # Assert that the model part has the correct number of nodes and elements
+        self.assertEqual(self.model_part.NumberOfNodes(), 8)
+        self.assertEqual(self.model_part.NumberOfConditions(), 12)
+
+        # Assert that the normals are the same as the nodes coordinates
+        for node in self.model_part.Nodes:
+            normal = node.GetValue(KratosMultiphysics.NORMAL)
+            self.assertAlmostEqual(normal[0], node.X)
+            self.assertAlmostEqual(normal[1], node.Y)
+            self.assertAlmostEqual(normal[2], node.Z)
+
+    def test_ReadObjIOTrianglesDegenerated(self):
+        """
+        Test the ReadModelPart function from ObjIO
+        """
+        # Create a model part and set the domain size
+        self.current_model = KratosMultiphysics.Model()
+        self.model_part = self.current_model.CreateModelPart("Main")
+        self.model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] = 3
+
+        # Read a model part from an OBJ file
+        obj_name = GetFilePath("auxiliar_files_for_python_unittest/obj_files/cube_degenerated.obj")
+        ReadModelPartFromOBJ(self.model_part, obj_name, True)
+
+        # # Debug
+        # debug_vtk(self.model_part)
+
+        # Assert that the model part has the correct number of nodes and elements
+        self.assertEqual(self.model_part.NumberOfNodes(), 8)
+        self.assertEqual(self.model_part.NumberOfConditions(), 12)
+
     def test_WriteObjIO(self):
         """
         Test the WriteModelPart function form ObjIO
@@ -206,13 +269,13 @@ def test_WriteObjIO(self):
         # Compute resulting area
         number_of_conditions = self.model_part.NumberOfConditions()
         area_2 = 0.0
-        for elem in obj_model_part.Elements:
+        for elem in obj_model_part.Conditions:
             geom = elem.GetGeometry()
             area_2 += geom.Area()
 
         # Assert number of nodes and elements
         self.assertEqual(obj_model_part.NumberOfNodes(), self.model_part.NumberOfNodes())
-        self.assertEqual(number_of_conditions, obj_model_part.NumberOfElements())
+        self.assertEqual(number_of_conditions, obj_model_part.NumberOfConditions())
 
         # Assert that the areas match approximately
         self.assertAlmostEqual(area_1, area_2)