Reintroduce mesh deletion option with faces circumcenters (#303 | GRI…

…DEDIT-917)

libs/MeshKernel/include/MeshKernel/Definitions.hpp

@@ -48,6 +48,9 @@ namespace meshkernel
     /// @brief Gets the valid projectionbs as vector of integers
     const std::vector<int>& GetValidProjections();
 
+    /// @brief Gets the valid deletion options as vector of integers
+    const std::vector<int>& GetValidDeletionOptions();
+
     /// @brief Convert an integer value to the Projection enumeration type
     ///
     /// If the integer projection value does not correspond to an enumeration

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -56,7 +56,8 @@ namespace meshkernel
         enum DeleteMeshOptions
         {
             InsideNotIntersected = 0,
-            InsideAndIntersected = 1
+            InsideAndIntersected = 1,
+            FacesWithIncludedCircumcenters = 2
         };
 
         /// Enumerator describing the different node types
@@ -271,7 +272,7 @@ namespace meshkernel
         ///                           If this Polygons instance contains multiple polygons, the first one will be taken.
         /// @param[in] deletionOption The deletion option
         /// @param[in] invertDeletion Inverts the selected node to delete (instead of outside the polygon, inside the polygon)
-        void DeleteMesh(const Polygons& polygon, int deletionOption, bool invertDeletion);
+        void DeleteMesh(const Polygons& polygon, DeleteMeshOptions deletionOption, bool invertDeletion);
 
         /// @brief Inquire if a segment is crossing a face
         /// @param[in] firstPoint The first point of the segment
@@ -338,6 +339,12 @@ namespace meshkernel
         /// @brief Bounded array for storing hanging node indices.
         using HangingNodeIndexArray = std::array<UInt, m_maximumNumberOfHangingNodesAlongEdge>;
 
+        /// @brief Deletes the mesh faces inside a polygon
+        /// @param[in] polygon        The polygon where to perform the operation
+        ///                           If this Polygons instance contains multiple polygons, the first one will be taken.
+        /// @param[in] invertDeletion Inverts the selected node to delete (instead of outside the polygon, inside the polygon)
+        void DeleteMeshFaces(const Polygons& polygon, bool invertDeletion);
+
         /// @brief Find cells recursive, works with an arbitrary number of edges
         /// @param[in] startNode The starting node
         /// @param[in] node The current node

libs/MeshKernel/src/Definitions.cpp

@@ -1,5 +1,6 @@
 #include "MeshKernel/Definitions.hpp"
 #include "MeshKernel/Exceptions.hpp"
+#include "MeshKernel/Mesh2D.hpp"
 
 namespace meshkernel
 {
@@ -12,6 +13,14 @@ namespace meshkernel
         return validProjections;
     }
 
+    const std::vector<int>& GetValidDeletionOptions()
+    {
+        static std::vector validMesh2DDeletionOptions{static_cast<int>(Mesh2D::DeleteMeshOptions::InsideNotIntersected),
+                                                      static_cast<int>(Mesh2D::DeleteMeshOptions::InsideAndIntersected),
+                                                      static_cast<int>(Mesh2D::DeleteMeshOptions::FacesWithIncludedCircumcenters)};
+        return validMesh2DDeletionOptions;
+    }
+
     Projection GetProjectionValue(const int projection)
     {
         static const int Cartesian = static_cast<int>(Projection::cartesian);

libs/MeshKernel/src/Mesh2D.cpp

@@ -1597,8 +1597,13 @@ std::vector<meshkernel::UInt> Mesh2D::GetHangingEdges() const
     return result;
 }
 
-void Mesh2D::DeleteMesh(const Polygons& polygon, int deletionOption, bool invertDeletion)
+void Mesh2D::DeleteMesh(const Polygons& polygon, DeleteMeshOptions deletionOption, bool invertDeletion)
 {
+    if (deletionOption == FacesWithIncludedCircumcenters)
+    {
+        DeleteMeshFaces(polygon, invertDeletion);
+        return;
+    }
     // Find crossed faces
     Mesh2DIntersections mesh2DIntersections(*this);
     mesh2DIntersections.Compute(polygon);
@@ -1675,6 +1680,65 @@ void Mesh2D::DeleteMesh(const Polygons& polygon, int deletionOption, bool invert
     Administrate();
 }
 
+void Mesh2D::DeleteMeshFaces(const Polygons& polygon, bool invertDeletion)
+{
+
+    Administrate();
+
+    for (UInt e = 0u; e < GetNumEdges(); ++e)
+    {
+        bool allFaceCircumcentersInPolygon = true;
+
+        for (UInt f = 0u; f < GetNumEdgesFaces(e); ++f)
+        {
+            const auto faceIndex = m_edgesFaces[e][f];
+            if (faceIndex == constants::missing::uintValue)
+            {
+                continue;
+            }
+
+            auto [isInPolygon, polygonIndex] = polygon.IsPointInPolygons(m_facesCircumcenters[faceIndex]);
+            if (invertDeletion)
+            {
+                isInPolygon = !isInPolygon;
+            }
+            if (!isInPolygon)
+            {
+                allFaceCircumcentersInPolygon = false;
+                break;
+            }
+        }
+
+        // 2D edge without surrounding faces.
+        if (GetNumEdgesFaces(e) == 0)
+        {
+            const auto firstNodeIndex = m_edges[e].first;
+            const auto secondNodeIndex = m_edges[e].second;
+
+            if (firstNodeIndex == constants::missing::uintValue || secondNodeIndex == constants::missing::uintValue)
+            {
+                continue;
+            }
+
+            const auto edgeCenter = (m_nodes[firstNodeIndex] + m_nodes[secondNodeIndex]) / 2.0;
+
+            auto [isInPolygon, polygonIndex] = polygon.IsPointInPolygons(edgeCenter);
+            allFaceCircumcentersInPolygon = isInPolygon;
+            if (invertDeletion)
+            {
+                allFaceCircumcentersInPolygon = !allFaceCircumcentersInPolygon;
+            }
+        }
+
+        if (allFaceCircumcentersInPolygon)
+        {
+            m_edges[e].first = constants::missing::uintValue;
+            m_edges[e].second = constants::missing::uintValue;
+        }
+    }
+    Administrate();
+}
+
 void Mesh2D::DeleteHangingEdges()
 {
     const auto hangingEdges = GetHangingEdges();

libs/MeshKernel/tests/src/Mesh2DTest.cpp

@@ -973,9 +973,10 @@ TEST(Mesh2D, DeleteMesh_WhenFacesAreIntersected_ShouldNotDeleteFaces)
         {1.87622950819672, -0.299180327868853}};
 
     auto polygon = meshkernel::Polygons(polygonNodes, meshkernel::Projection::cartesian);
+    const auto deletion_option = meshkernel::Mesh2D::DeleteMeshOptions::InsideNotIntersected;
 
     // Execute
-    mesh->DeleteMesh(polygon, 0, false);
+    mesh->DeleteMesh(polygon, deletion_option, false);
 
     // Assert
     EXPECT_EQ(mesh->GetNumFaces(), 9);
@@ -999,9 +1000,10 @@ TEST(Mesh2D, DeleteMesh_WhenFacesAreIntersectedSpherical_ShouldNotDeleteFaces)
         {1.87622950819672, -0.299180327868853}};
 
     auto polygon = meshkernel::Polygons(polygonNodes, meshkernel::Projection::spherical);
+    const auto deletion_option = meshkernel::Mesh2D::DeleteMeshOptions::InsideNotIntersected;
 
     // Execute
-    mesh->DeleteMesh(polygon, 0, false);
+    mesh->DeleteMesh(polygon, deletion_option, false);
 
     // Assert
     EXPECT_EQ(mesh->GetNumFaces(), 9);
@@ -1026,10 +1028,38 @@ TEST(Mesh2D, DeleteMesh_WithLargeSphericalPolygon_ShouldDeleteInnerMeshFaces)
         {-2.29490103397341, 50.0126381093058}};
 
     auto polygon = meshkernel::Polygons(polygonNodes, meshkernel::Projection::spherical);
+    const auto deletion_option = meshkernel::Mesh2D::DeleteMeshOptions::InsideNotIntersected;
 
     // Execute
-    mesh->DeleteMesh(polygon, 0, false);
+    mesh->DeleteMesh(polygon, deletion_option, false);
 
     // Assert
     EXPECT_EQ(mesh->GetNumFaces(), 7);
 }
+
+TEST(Mesh2D, DeleteMesh_WithPolygonAndIncludedCircumcenters_ShouldDeleteInnerFaces)
+{
+    // Prepare
+    const auto mesh = MakeRectangularMeshForTesting(5,
+                                                    5,
+                                                    8.0,
+                                                    8.0,
+                                                    meshkernel::Projection::cartesian);
+
+    // a large polygon
+    std::vector<meshkernel::Point> polygonNodes{
+        {2, 2},
+        {6, 2},
+        {6, 6},
+        {2, 6},
+        {2, 2}};
+
+    auto polygon = meshkernel::Polygons(polygonNodes, meshkernel::Projection::cartesian);
+    const auto deletion_option = meshkernel::Mesh2D::DeleteMeshOptions::FacesWithIncludedCircumcenters;
+
+    // Execute
+    mesh->DeleteMesh(polygon, deletion_option, false);
+
+    // Assert
+    EXPECT_EQ(mesh->GetNumFaces(), 12);
+}

libs/MeshKernelApi/src/MeshKernel.cpp

@@ -176,12 +176,15 @@ namespace meshkernelapi
             {
                 throw meshkernel::ConstraintError("The 2d mesh contains no nodes.");
             }
+            meshkernel::range_check::CheckOneOf<int>(deletionOption, meshkernel::GetValidDeletionOptions(), "Deletion");
 
             const auto polygonPoints = ConvertGeometryListToPointVector(polygon);
 
             const bool invertDeletionBool = invertDeletion == 1 ? true : false;
             const meshkernel::Polygons meshKernelPolygon(polygonPoints, meshKernelState[meshKernelId].m_mesh2d->m_projection);
-            meshKernelState[meshKernelId].m_mesh2d->DeleteMesh(meshKernelPolygon, deletionOption, invertDeletionBool);
+            const auto deletionOptionEnum = static_cast<meshkernel::Mesh2D::DeleteMeshOptions>(deletionOption);
+
+            meshKernelState[meshKernelId].m_mesh2d->DeleteMesh(meshKernelPolygon, deletionOptionEnum, invertDeletionBool);
         }
         catch (...)
         {

libs/MeshKernelApi/tests/src/ApiTest.cpp

@@ -1150,7 +1150,7 @@ TEST_F(CartesianApiTestFixture, MakeCurvilinearGridFromTriangleThroughApi)
     ASSERT_EQ(40, mesh2d.num_edges);
 }
 
-TEST_F(CartesianApiTestFixture, DeleteMesh2D_WithEmptyPolygon_ShouldDeleteMesh2D)
+TEST_F(CartesianApiTestFixture, Delete_WithEmptyPolygon_ShouldDeleteMesh2D)
 {
     // Prepare
     MakeMesh();
@@ -1170,6 +1170,34 @@ TEST_F(CartesianApiTestFixture, DeleteMesh2D_WithEmptyPolygon_ShouldDeleteMesh2D
     ASSERT_EQ(0, mesh2d.num_edges);
 }
 
+TEST_F(CartesianApiTestFixture, DeleteFaces_WithPolygon_ShouldDeleteMesh2D)
+{
+    // Prepare
+    MakeMesh(4, 4, 2);
+    auto const meshKernelId = GetMeshKernelId();
+
+    // By using an empty list, all nodes will be selected
+    meshkernelapi::GeometryList geometryList{};
+
+    geometryList.num_coordinates = 5;
+    geometryList.geometry_separator = meshkernel::constants::missing::doubleValue;
+    std::vector<double> xCoordinatesOut{2, 6, 6, 2, 2};
+    std::vector<double> yCoordinatesOut{2, 2, 6, 6, 2};
+    geometryList.coordinates_x = xCoordinatesOut.data();
+    geometryList.coordinates_y = yCoordinatesOut.data();
+
+    // Execute
+    int deletionOption = 2;
+    auto errorCode = mkernel_mesh2d_delete(meshKernelId, geometryList, deletionOption, false);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    meshkernelapi::Mesh2D mesh2d{};
+    errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, mesh2d);
+
+    // Assert
+    ASSERT_EQ(12, mesh2d.num_faces);
+}
+
 TEST_F(CartesianApiTestFixture, GetDimensionsMesh1D_WithMesh1D_ShouldGetDimensionsMesh1D)
 {
     // Prepare