correct merge nodes algorithm

libs/MeshKernel/include/MeshKernel/Mesh.hpp

@@ -228,6 +228,12 @@ namespace meshkernel
         /// @brief Compute the lengths of all edges in one go
         void ComputeEdgesLengths();
 
+        /// @brief Compute the minimum edge length of the edges included in the polygon.
+        /// An edge is considered included if one of the two nodes is inside the polygon.
+        /// @param[in] polygon The polygon for considering an edge included
+        /// @return The minimum edge length
+        [[nodiscard]] double ComputeMinEdgeLength(const Polygons& polygon) const;
+
         /// @brief Computes the edges centers  in one go
         void ComputeEdgesCenters();
 
@@ -299,6 +305,12 @@ namespace meshkernel
         /// @return The vector with the mesh locations.
         [[nodiscard]] std::vector<Point> ComputeLocations(Location location) const;
 
+        /// @brief Computes if a location is in polygon.
+        /// @param[in] polygon The input polygon.
+        /// @param[in] location The mesh location (e.g. nodes, edge centers or face circumcenters).
+        /// @return A vector of booleans indicating if a location is in a polygon or not.
+        [[nodiscard]] std::vector<bool> IsLocationInPolygon(const Polygons& polygon, Location location) const;
+
         /// @brief Add meshes: result is a mesh composed of the additions
         /// firstMesh += secondmesh results in the second mesh being added to firstMesh
         /// @param[in] rhs The mesh to add

libs/MeshKernel/src/Mesh.cpp

@@ -274,20 +274,29 @@ void Mesh::MergeTwoNodes(UInt firstNodeIndex, UInt secondNodeIndex)
 void Mesh::MergeNodesInPolygon(const Polygons& polygon, double mergingDistance)
 {
     // first filter the nodes in polygon
-    std::vector<Point> filteredNodes(GetNumNodes());
-    std::vector<UInt> originalNodeIndices(GetNumNodes(), constants::missing::uintValue);
-    UInt index = 0;
-    for (UInt i = 0; i < GetNumNodes(); ++i)
+    const auto numNodes = GetNumNodes();
+    std::vector<Point> filteredNodes(numNodes);
+    std::vector<UInt> originalNodeIndices(numNodes, constants::missing::uintValue);
+    const auto isNodeInPolygon = IsLocationInPolygon(polygon, Location::Nodes);
+
+    UInt filteredNodeCount = 0;
+    for (UInt i = 0; i < numNodes; ++i)
     {
-        const bool inPolygon = polygon.IsPointInPolygon(m_nodes[i], 0);
-        if (inPolygon)
+        if (isNodeInPolygon[i])
         {
-            filteredNodes[index] = m_nodes[i];
-            originalNodeIndices[index] = i;
-            index++;
+            filteredNodes[filteredNodeCount] = m_nodes[i];
+            originalNodeIndices[filteredNodeCount] = i;
+            filteredNodeCount++;
         }
     }
-    filteredNodes.resize(index);
+
+    // no node to merge
+    if (filteredNodeCount == 0)
+    {
+        return;
+    }
+
+    filteredNodes.resize(filteredNodeCount);
 
     // Update the R-Tree of the mesh nodes
     RTree nodesRtree;
@@ -395,6 +404,23 @@ void Mesh::ComputeEdgesLengths()
     }
 }
 
+double Mesh::ComputeMinEdgeLength(const Polygons& polygon) const
+{
+    auto const numEdges = GetNumEdges();
+    auto result = std::numeric_limits<double>::max();
+
+    const auto isNodeInPolygon = IsLocationInPolygon(polygon, Location::Nodes);
+    for (UInt e = 0; e < numEdges; e++)
+    {
+        const auto& [firstNode, secondNode] = m_edges[e];
+        if (isNodeInPolygon[firstNode] || isNodeInPolygon[secondNode])
+        {
+            result = std::min(result, m_edgeLengths[e]);
+        }
+    }
+    return result;
+}
+
 void Mesh::ComputeEdgesCenters()
 {
     m_edgesCenters = ComputeEdgeCenters(m_nodes, m_edges);
@@ -459,7 +485,7 @@ meshkernel::UInt Mesh::FindNodeCloseToAPoint(Point const& point, double searchRa
         return GetLocationsIndices(0, Location::Nodes);
     }
 
-    throw AlgorithmError("Could not find the node index close to a point.");
+    return constants::missing::uintValue;
 }
 
 meshkernel::UInt Mesh::FindNodeCloseToAPoint(Point point, const std::vector<bool>& oneDNodeMask)
@@ -878,6 +904,18 @@ std::vector<meshkernel::Point> Mesh::ComputeLocations(Location location) const
     return result;
 }
 
+std::vector<bool> Mesh::IsLocationInPolygon(const Polygons& polygon, Location location) const
+{
+    const auto locations = ComputeLocations(location);
+    std::vector<bool> result(locations.size(), false);
+    for (UInt i = 0; i < result.size(); ++i)
+    {
+        result[i] = polygon.IsPointInPolygon(locations[i], 0);
+    }
+
+    return result;
+}
+
 Mesh& Mesh::operator+=(Mesh const& rhs)
 {
     if (m_projection != rhs.m_projection)

libs/MeshKernel/src/Utilities/RTree.cpp

@@ -113,7 +113,7 @@ void RTree::DeleteNode(UInt position)
     const auto numberRemoved = m_rtree2D.remove(m_points[position]);
     if (numberRemoved != 1)
     {
-        throw std::invalid_argument("DeleteNode: Could not remove node at given position.");
+        return;
     }
     m_points[position] = {Point2D{constants::missing::doubleValue, constants::missing::doubleValue}, std::numeric_limits<UInt>::max()};
 }

libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp

@@ -1061,12 +1061,18 @@ namespace meshkernelapi
         MKERNEL_API int mkernel_mesh2d_make_rectangular_mesh_on_extension(int meshKernelId,
                                                                           const meshkernel::MakeGridParameters& makeGridParameters);
 
+        /// @brief Merges the mesh2d nodes within a distance of 0.001 m, effectively removing all small edges
+        /// @param[in] meshKernelId   The id of the mesh state
+        /// @param[in] geometryListIn The polygon defining the area where the operation will be performed
+        /// @returns Error code
+        MKERNEL_API int mkernel_mesh2d_merge_nodes(int meshKernelId, const GeometryList& geometryListIn);
+
         /// @brief Merges the mesh2d nodes within a distance of 0.001 m, effectively removing all small edges
         /// @param[in] meshKernelId   The id of the mesh state
         /// @param[in] geometryListIn The polygon defining the area where the operation will be performed
         /// @param[in] mergingDistance The distance below which two nodes will be merged
         /// @returns Error code
-        MKERNEL_API int mkernel_mesh2d_merge_nodes(int meshKernelId, const GeometryList& geometryListIn, double mergingDistance);
+        MKERNEL_API int mkernel_mesh2d_merge_nodes_with_merging_distance(int meshKernelId, const GeometryList& geometryListIn, double mergingDistance);
 
         /// @brief Merges two mesh2d nodes into one.
         /// @param[in] meshKernelId The id of the mesh state

libs/MeshKernelApi/src/MeshKernel.cpp

@@ -1175,7 +1175,33 @@ namespace meshkernelapi
         return lastExitCode;
     }
 
-    MKERNEL_API int mkernel_mesh2d_merge_nodes(int meshKernelId, const GeometryList& geometryListIn, double mergingDistance)
+    MKERNEL_API int mkernel_mesh2d_merge_nodes(int meshKernelId, const GeometryList& geometryListIn)
+    {
+        lastExitCode = meshkernel::ExitCode::Success;
+        try
+        {
+            if (!meshKernelState.contains(meshKernelId))
+            {
+                throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist.");
+            }
+
+            auto const polygonVector = ConvertGeometryListToPointVector(geometryListIn);
+
+            const meshkernel::Polygons polygon(polygonVector, meshKernelState[meshKernelId].m_mesh2d->m_projection);
+            meshKernelState[meshKernelId].m_mesh2d->ComputeEdgesLengths();
+
+            const auto minEdgeLength = meshKernelState[meshKernelId].m_mesh2d->ComputeMinEdgeLength(polygon);
+            const auto searchRadius = std::max(1e-6, minEdgeLength * 0.1);
+            meshKernelState[meshKernelId].m_mesh2d->MergeNodesInPolygon(polygon, searchRadius);
+        }
+        catch (...)
+        {
+            lastExitCode = HandleException();
+        }
+        return lastExitCode;
+    }
+
+    MKERNEL_API int mkernel_mesh2d_merge_nodes_with_merging_distance(int meshKernelId, const GeometryList& geometryListIn, double mergingDistance)
     {
         lastExitCode = meshkernel::ExitCode::Success;
         try

libs/MeshKernelApi/tests/src/ApiTest.cpp

@@ -219,9 +219,38 @@ TEST_F(CartesianApiTestFixture, MergeNodesThroughApi)
     MakeMesh();
     auto const meshKernelId = GetMeshKernelId();
     meshkernelapi::GeometryList geometry_list{};
+    std::vector xCoordinates{-0.5, 2.5, 2.5, -0.5, -0.5};
+    std::vector yCoordinates{-0.5, -0.5, 2.5, 2.5, -0.5};
+    std::vector zCoordinates{0.0, 0.0, 0.0, 0.5, 0.5};
+
+    geometry_list.geometry_separator = meshkernel::constants::missing::doubleValue;
+    geometry_list.coordinates_x = xCoordinates.data();
+    geometry_list.coordinates_y = yCoordinates.data();
+    geometry_list.values = zCoordinates.data();
+    geometry_list.num_coordinates = static_cast<int>(xCoordinates.size());
+
+    // Execute
+    auto errorCode = mkernel_mesh2d_merge_nodes(meshKernelId, geometry_list);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    meshkernelapi::Mesh2D mesh2d{};
+    errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, mesh2d);
+
+    // Assert (nothing is done, just check that the api communication works)
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+    ASSERT_EQ(12, mesh2d.num_nodes);
+    ASSERT_EQ(17, mesh2d.num_edges);
+}
+
+TEST_F(CartesianApiTestFixture, MergeNodesWithMergingDistanceThroughApi)
+{
+    // Prepare
+    MakeMesh(2, 3, 1.0);
+    auto const meshKernelId = GetMeshKernelId();
+    meshkernelapi::GeometryList geometry_list{};
 
     // Execute
-    auto errorCode = mkernel_mesh2d_merge_nodes(meshKernelId, geometry_list, 0.001);
+    auto errorCode = mkernel_mesh2d_merge_nodes_with_merging_distance(meshKernelId, geometry_list, 0.001);
     ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
 
     meshkernelapi::Mesh2D mesh2d{};