diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp index fd5537e98..af8633e93 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh.hpp +++ b/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 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 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 diff --git a/libs/MeshKernel/src/Mesh.cpp b/libs/MeshKernel/src/Mesh.cpp index 8022f023b..75c6a172d 100644 --- a/libs/MeshKernel/src/Mesh.cpp +++ b/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 filteredNodes(GetNumNodes()); - std::vector originalNodeIndices(GetNumNodes(), constants::missing::uintValue); - UInt index = 0; - for (UInt i = 0; i < GetNumNodes(); ++i) + const auto numNodes = GetNumNodes(); + std::vector filteredNodes(numNodes); + std::vector 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::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& oneDNodeMask) @@ -878,6 +904,18 @@ std::vector Mesh::ComputeLocations(Location location) const return result; } +std::vector Mesh::IsLocationInPolygon(const Polygons& polygon, Location location) const +{ + const auto locations = ComputeLocations(location); + std::vector 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) diff --git a/libs/MeshKernel/src/Utilities/RTree.cpp b/libs/MeshKernel/src/Utilities/RTree.cpp index a22cf8261..7b04a8c9d 100644 --- a/libs/MeshKernel/src/Utilities/RTree.cpp +++ b/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::max()}; } diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 0150f3447..9717a633b 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/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 diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index d0285c451..81071d60c 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/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 diff --git a/libs/MeshKernelApi/tests/src/ApiTest.cpp b/libs/MeshKernelApi/tests/src/ApiTest.cpp index 50ad9d816..ce84f64d5 100644 --- a/libs/MeshKernelApi/tests/src/ApiTest.cpp +++ b/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(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{};