From 636116137c976d6e053e2804a32e96d3a0789455 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 2 Jul 2026 17:49:51 +0200 Subject: [PATCH 01/34] GRIDEDIT-2292 Removed debugging output when saving vtk files --- libs/MeshKernel/src/Utilities/Utilities.cpp | 24 --------------------- 1 file changed, 24 deletions(-) diff --git a/libs/MeshKernel/src/Utilities/Utilities.cpp b/libs/MeshKernel/src/Utilities/Utilities.cpp index 28cdec407..ada3a13b1 100644 --- a/libs/MeshKernel/src/Utilities/Utilities.cpp +++ b/libs/MeshKernel/src/Utilities/Utilities.cpp @@ -103,9 +103,6 @@ void meshkernel::SaveVtk(const std::vector& nodes, const std::vector unsavedElements; - - unsavedElements.fill(0); UInt numberOfElements = 0; @@ -116,16 +113,8 @@ void meshkernel::SaveVtk(const std::vector& nodes, const std::vector& nodesX, const std::vector unsavedElements; - - unsavedElements.fill(0); UInt numberOfElements = 0; @@ -270,16 +256,6 @@ void meshkernel::SaveVtk(const std::vector& nodesX, const std::vector Date: Thu, 2 Jul 2026 17:50:48 +0200 Subject: [PATCH 02/34] GRIDEDIT-2292 Reconstructed the invalid cell polygons to unify polygons that are covering the same patch --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 8 ++- libs/MeshKernel/src/Mesh2D.cpp | 53 ++++++++++++++++++- 2 files changed, 58 insertions(+), 3 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index f1f163975..b24cdc645 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -252,7 +252,7 @@ namespace meshkernel /// @brief Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol) /// @param[in] polygon The polygon where the operation is performed; only boundary segments intersecting the polygon are included /// @return The resulting polygon mesh boundary - [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon); + [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon, const bool doAdministrate = true); /// @brief Convert all mesh boundaries to a vector of polygon nodes /// @return The resulting set of polygons, describing interior mesh boundaries @@ -461,6 +461,12 @@ namespace meshkernel std::vector& subSequence, std::vector& illegalCells) const; + /// @brief Reconstruct the invalid cell polygons + /// + /// When constructing the invalid cell polygons, they can be compujted with many smaller polygons. + /// If these smaller polygons form a single patch on the domain, then they need to be combined + void ReconstructInvalidCellsPolygon(); + /// @brief Ensure that all polynomials are orientated in the ACW direction. void OrientatePolygonsAntiClockwise(std::vector& polygonNodes) const; diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index e5fa3f6c7..51314ade4 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1510,12 +1510,16 @@ std::vector Mesh2D::SortedFacesAroundNode(UInt node) const return result; } -std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector& polygonNodes) +std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector& polygonNodes, const bool doAdministrate) { const Polygon polygon(polygonNodes, m_projection); // Find faces - Administrate(); + if (doAdministrate) + { + Administrate(); + } + std::vector isVisited(GetNumEdges(), false); std::vector meshBoundaryPolygon; meshBoundaryPolygon.reserve(GetNumNodes()); @@ -1572,6 +1576,7 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector meshBoundaryPolygon.push_back(meshBoundaryPolygon.front()); } } + return meshBoundaryPolygon; } @@ -2183,6 +2188,48 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygon(const P return UpdateFaceInformation(faceIndices, appendDeletedFaces); } +void Mesh2D::ReconstructInvalidCellsPolygon() +{ + std::vector allBoundaries(ComputeBoundaryPolygons(std::vector(), false)); + + if (allBoundaries.GetNumPolygons() <= 1) + { + // There are no interior boundary polygons + return; + } + + Polygons polygons(allBoundaries, m_projection); + std::vector polygonAreas(polygons.GetNumPolygons(), 0.0); + + for (UInt p = 0; p < polygons.GetNumPolygons(); ++p) + { + auto [area, centre, direction] = polygons.Enclosure(p).Outer().FaceAreaAndCenterOfMass(); + polygonAreas[p] = area; + } + + auto maxAreaIter = std::ranges::max_element(polygonAreas); + size_t maxAreaIndex = std::distance(polygonAreas.begin(), maxAreaIter); + + std::vector innerBoundaryPolygons; + innerBoundaryPolygons.reserve(m_invalidCellPolygons.size()); + + for (size_t p = 0; p < polygons.GetNumPolygons(); ++p) + { + if (p != maxAreaIndex) + { + const std::vector& polygonPoints(polygons.Enclosure(p).Outer().Nodes()); + innerBoundaryPolygons.insert(innerBoundaryPolygons.end(), polygonPoints.begin(), polygonPoints.end()); + + if ((p < maxAreaIndex && ((maxAreaIndex + 1) != polygons.GetNumPolygons())) || (p > maxAreaIndex && (p + 1) != polygons.GetNumPolygons())) + { + innerBoundaryPolygons.push_back(Point(constants::missing::doubleValue, constants::missing::doubleValue)); + } + } + } + + m_invalidCellPolygons = innerBoundaryPolygons; +} + std::unique_ptr Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, const bool appendDeletedFaces) { std::vector facesToDelete; @@ -2262,6 +2309,8 @@ std::unique_ptr Mesh2D::UpdateFaceInformation(const std: } } + ReconstructInvalidCellsPolygon(); + // Shift face connectivity in arrays where deleted faces have been removed from arrays // Loop must be iterated in reverse order, from highest face id value to lowest. for (UInt faceId : facesToDelete | std::views::reverse) From c006610688fd61b59890d0bd90b3469dab285cef Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 2 Jul 2026 17:59:26 +0200 Subject: [PATCH 03/34] GRIDEDIT-2292 Fixed doxygen and clang formatting warnings --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 1 + libs/MeshKernel/src/Utilities/Utilities.cpp | 2 -- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index b24cdc645..61206cac7 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -251,6 +251,7 @@ namespace meshkernel /// @brief Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol) /// @param[in] polygon The polygon where the operation is performed; only boundary segments intersecting the polygon are included + /// @param[in] doAdministrate Indicate that an adminstrate is required. /// @return The resulting polygon mesh boundary [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon, const bool doAdministrate = true); diff --git a/libs/MeshKernel/src/Utilities/Utilities.cpp b/libs/MeshKernel/src/Utilities/Utilities.cpp index ada3a13b1..32151435e 100644 --- a/libs/MeshKernel/src/Utilities/Utilities.cpp +++ b/libs/MeshKernel/src/Utilities/Utilities.cpp @@ -112,10 +112,8 @@ void meshkernel::SaveVtk(const std::vector& nodes, const std::vector Date: Mon, 6 Jul 2026 13:27:21 +0200 Subject: [PATCH 04/34] GRIDEDIT-2292 Refactored method for determining enclosing and non-enclosing polygons --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 8 +- .../MeshKernel/include/MeshKernel/Polygon.hpp | 5 ++ libs/MeshKernel/src/Mesh2D.cpp | 83 +++++++++++-------- libs/MeshKernel/src/Polygon.cpp | 6 ++ .../tests/src/MeshRefinementTests.cpp | 68 +++------------ .../tests/src/InvalidCellsPolygonsTests.cpp | 16 +++- 6 files changed, 90 insertions(+), 96 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 61206cac7..04a156888 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -253,7 +253,7 @@ namespace meshkernel /// @param[in] polygon The polygon where the operation is performed; only boundary segments intersecting the polygon are included /// @param[in] doAdministrate Indicate that an adminstrate is required. /// @return The resulting polygon mesh boundary - [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon, const bool doAdministrate = true); + [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon); /// @brief Convert all mesh boundaries to a vector of polygon nodes /// @return The resulting set of polygons, describing interior mesh boundaries @@ -468,6 +468,12 @@ namespace meshkernel /// If these smaller polygons form a single patch on the domain, then they need to be combined void ReconstructInvalidCellsPolygon(); + /// @brief Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol) + /// + /// @return a sequence of boundary points, which may be separated by the invalid point, and a matching sequence of Boolean values indicating which + /// polygnal sub-sequence forms a external boudary + [[nodiscard]] std::tuple, std::vector> GetAllBoundaryPolygons(const std::vector& polygon); + /// @brief Ensure that all polynomials are orientated in the ACW direction. void OrientatePolygonsAntiClockwise(std::vector& polygonNodes) const; diff --git a/libs/MeshKernel/include/MeshKernel/Polygon.hpp b/libs/MeshKernel/include/MeshKernel/Polygon.hpp index 2c390c295..8b7edb995 100644 --- a/libs/MeshKernel/include/MeshKernel/Polygon.hpp +++ b/libs/MeshKernel/include/MeshKernel/Polygon.hpp @@ -27,6 +27,7 @@ #pragma once +#include #include #include "MeshKernel/BoundingBox.hpp" @@ -54,6 +55,10 @@ namespace meshkernel /// @brief Default move constructor. Polygon(Polygon&& copy) = default; + /// @brief Constructor + Polygon(std::span points, + Projection projection); + /// @brief Constructor Polygon(const std::vector& points, Projection projection); diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 51314ade4..631e48f1b 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1510,18 +1510,21 @@ std::vector Mesh2D::SortedFacesAroundNode(UInt node) const return result; } -std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector& polygonNodes, const bool doAdministrate) +std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector& polygonNodes) { - const Polygon polygon(polygonNodes, m_projection); + Administrate(); + auto [boundaryPoints, isEnclosingBoundary] = GetAllBoundaryPolygons(polygonNodes); + return boundaryPoints; +} - // Find faces - if (doAdministrate) - { - Administrate(); - } +std::tuple, std::vector> Mesh2D::GetAllBoundaryPolygons(const std::vector& polygonNodes) +{ + const Polygon polygon(polygonNodes, m_projection); std::vector isVisited(GetNumEdges(), false); std::vector meshBoundaryPolygon; + std::vector isEnclosingBoundary; + meshBoundaryPolygon.reserve(GetNumNodes()); for (UInt e = 0; e < GetNumEdges(); e++) @@ -1531,10 +1534,10 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector continue; } - const auto firstNodeIndex = m_edges[e].first; - const auto secondNodeIndex = m_edges[e].second; - const auto firstNode = m_nodes[firstNodeIndex]; - const auto secondNode = m_nodes[secondNodeIndex]; + const UInt firstNodeIndex = m_edges[e].first; + const UInt secondNodeIndex = m_edges[e].second; + const Point firstNode = m_nodes[firstNodeIndex]; + const Point secondNode = m_nodes[secondNodeIndex]; bool firstNodeInPolygon = polygon.Contains(m_nodes[firstNodeIndex]); bool secondNodeInPolygon = polygon.Contains(m_nodes[secondNodeIndex]); @@ -1550,9 +1553,12 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); } + const size_t boundaryPolygonStartIndex = meshBoundaryPolygon.size(); + // Put the current edge on the mesh boundary, mark it as visited meshBoundaryPolygon.emplace_back(firstNode); meshBoundaryPolygon.emplace_back(secondNode); + isVisited[e] = true; // walk the current mesh boundary @@ -1575,9 +1581,17 @@ std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector std::reverse(meshBoundaryPolygon.begin() + numNodesFirstTail, meshBoundaryPolygon.end()); meshBoundaryPolygon.push_back(meshBoundaryPolygon.front()); } + + const size_t boundaryPolygonEndIndex = meshBoundaryPolygon.size(); + + std::span currentPolygonSpan(std::span(meshBoundaryPolygon.data() + boundaryPolygonStartIndex, + meshBoundaryPolygon.data() + boundaryPolygonEndIndex)); + Polygon currentPolygon(currentPolygonSpan, m_projection); + + isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[m_edgesFaces[e][0]])); } - return meshBoundaryPolygon; + return {meshBoundaryPolygon, isEnclosingBoundary}; } std::vector Mesh2D::ComputeInnerBoundaryPolygons() const @@ -2190,44 +2204,41 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygon(const P void Mesh2D::ReconstructInvalidCellsPolygon() { - std::vector allBoundaries(ComputeBoundaryPolygons(std::vector(), false)); + auto [boundaryPoints, isEnclosingBoundary] = GetAllBoundaryPolygons(std::vector()); + + Polygons polygons(boundaryPoints, m_projection); - if (allBoundaries.GetNumPolygons() <= 1) + if (polygons.GetNumPolygons() <= 1) { // There are no interior boundary polygons return; } - Polygons polygons(allBoundaries, m_projection); - std::vector polygonAreas(polygons.GetNumPolygons(), 0.0); - - for (UInt p = 0; p < polygons.GetNumPolygons(); ++p) - { - auto [area, centre, direction] = polygons.Enclosure(p).Outer().FaceAreaAndCenterOfMass(); - polygonAreas[p] = area; - } - - auto maxAreaIter = std::ranges::max_element(polygonAreas); - size_t maxAreaIndex = std::distance(polygonAreas.begin(), maxAreaIter); - std::vector innerBoundaryPolygons; innerBoundaryPolygons.reserve(m_invalidCellPolygons.size()); + bool firstElement = true; for (size_t p = 0; p < polygons.GetNumPolygons(); ++p) { - if (p != maxAreaIndex) + if (isEnclosingBoundary[p]) { - const std::vector& polygonPoints(polygons.Enclosure(p).Outer().Nodes()); - innerBoundaryPolygons.insert(innerBoundaryPolygons.end(), polygonPoints.begin(), polygonPoints.end()); + continue; + } - if ((p < maxAreaIndex && ((maxAreaIndex + 1) != polygons.GetNumPolygons())) || (p > maxAreaIndex && (p + 1) != polygons.GetNumPolygons())) - { - innerBoundaryPolygons.push_back(Point(constants::missing::doubleValue, constants::missing::doubleValue)); - } + if (!firstElement) + { + innerBoundaryPolygons.push_back(Point(constants::missing::doubleValue, constants::missing::doubleValue)); } + else + { + firstElement = false; + } + + const std::vector& polygonPoints(polygons.Enclosure(p).Outer().Nodes()); + innerBoundaryPolygons.insert(innerBoundaryPolygons.end(), polygonPoints.begin(), polygonPoints.end()); } - m_invalidCellPolygons = innerBoundaryPolygons; + m_invalidCellPolygons = std::move(innerBoundaryPolygons); } std::unique_ptr Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, const bool appendDeletedFaces) @@ -2309,8 +2320,6 @@ std::unique_ptr Mesh2D::UpdateFaceInformation(const std: } } - ReconstructInvalidCellsPolygon(); - // Shift face connectivity in arrays where deleted faces have been removed from arrays // Loop must be iterated in reverse order, from highest face id value to lowest. for (UInt faceId : facesToDelete | std::views::reverse) @@ -2322,6 +2331,8 @@ std::unique_ptr Mesh2D::UpdateFaceInformation(const std: m_faceArea.erase(m_faceArea.begin() + faceId); } + ReconstructInvalidCellsPolygon(); + return deleteMeshAction; } diff --git a/libs/MeshKernel/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp index f72fd156a..119a9186d 100644 --- a/libs/MeshKernel/src/Polygon.cpp +++ b/libs/MeshKernel/src/Polygon.cpp @@ -42,6 +42,12 @@ meshkernel::Polygon::Polygon(const std::vector& points, Initialise(); } +meshkernel::Polygon::Polygon(std::span points, + Projection projection) : m_nodes(points.begin(), points.end()), m_projection(projection) +{ + Initialise(); +} + meshkernel::Polygon::Polygon(std::vector&& points, Projection projection) : m_nodes(points), m_projection(projection) { diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index e22d177eb..948780dd8 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -3082,67 +3082,23 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon std::vector innerBoundaryPoints = mesh2.GetInnerBoundaryPolygons(); // The expected number of points, should not include any land boundary points - UInt expectedNumberOfNodes = 26; + constexpr UInt expectedNumberOfNodes = 22; ASSERT_EQ(expectedNumberOfNodes, innerBoundaryPoints.size()); // The edge of one of the deleted elements lies on the boundary, so will be not be part of the // interior set of polygons - std::vector expectedXPoints{ - 95.0, - 105.0, - 100.0, - 95.0, - constants::missing::doubleValue, - 85.0, - 85.0, - 95.0, - 95.0, - 85.0, - constants::missing::doubleValue, - 80.0, - 75.0, - 85.0, - 80.0, - constants::missing::doubleValue, - 120.0, - 115.0, - 125.0, - 120.0, - constants::missing::doubleValue, - 125.0, - 125.0, - 135.0, - 135.0, - 125.0}; - - std::vector expectedYPoints{ - 15.0, - 15.0, - 0.0, - 15.0, - constants::missing::doubleValue, - 15.0, - 25.0, - 25.0, - 15.0, - 15.0, - constants::missing::doubleValue, - 0.0, - 15.0, - 15.0, - 0.0, - constants::missing::doubleValue, - 0.0, - 15.0, - 15.0, - 0.0, - constants::missing::doubleValue, - 125.0, - 135.0, - 135.0, - 125.0, - 125.0}; + std::vector expectedXPoints{80.0, 85.0, 95.0, 100.0, 105.0, 95.0, 95.0, 85.0, 85.0, 75.0, 80.0, + constants::missing::doubleValue, + 120.0, 125.0, 115.0, 120.0, + constants::missing::doubleValue, + 125.0, 125.0, 135.0, 135.0, 125.0}; + + std::vector expectedYPoints{0.0, 15.0, 15.0, 0.0, 15.0, 15.0, 25.0, 25.0, 15.0, 15.0, 0.0, + constants::missing::doubleValue, + 0.0, 15.0, 15.0, 0.0, + constants::missing::doubleValue, + 125.0, 135.0, 135.0, 125.0, 125.0}; for (size_t i = 0; i < expectedXPoints.size(); ++i) { diff --git a/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp index 33c5987c4..97d77b1fd 100644 --- a/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp +++ b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp @@ -191,15 +191,25 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) //----------------------- // Check the inner boundary polygon are correct - std::vector expectedInnerX{80.0, 85.0, 75.0, 80.0, -999.0, 100.0, 105.0, 95.0, 100.0, -999.0, 120.0, 125.0, 115.0, 120.0, -999.0, 180.0, 200.0, 200.0, 180.0, 180.0, -999.0, 95.0, 95.0, 85.0, 85.0, 95.0, -999.0, 125.0, 135.0, 135.0, 125.0, 125.0}; - std::vector expectedInnerY{0.0, 15.0, 15.0, 0.0, -999.0, 0.0, 15.0, 15.0, 0.0, -999.0, 0.0, 15.0, 15.0, 0.0, -999.0, 140.0, 140.0, 160.0, 160.0, 140.0, -999.0, 15.0, 25.0, 25.0, 15.0, 15.0, -999.0, 125.0, 125.0, 135.0, 135.0, 125.0}; + std::vector expectedInnerX{80.0, 85.0, 95.0, 100.0, 105.0, 95.0, 95.0, 85.0, 85.0, 75.0, 80.0, + meshkernel::constants::missing::doubleValue, + 120.0, 125.0, 115.0, 120.0, + meshkernel::constants::missing::doubleValue, + 125.0, 125.0, 135.0, 135.0, 125.0}; + + std::vector expectedInnerY{0.0, 15.0, 15.0, 0.0, 15.0, 15.0, 25.0, 25.0, 15.0, 15.0, 0.0, + meshkernel::constants::missing::doubleValue, + 0.0, 15.0, 15.0, 0.0, + meshkernel::constants::missing::doubleValue, + 125.0, 135.0, 135.0, 125.0, 125.0}; int innerPolygonSize = 0; meshkernelapi::GeometryList innerPolygon; errorCode = meshkernelapi::mkernel_mesh2d_get_mesh_inner_boundaries_as_polygons_dimension(meshKernelId, innerPolygonSize); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); - ASSERT_EQ(innerPolygonSize, 32); + + ASSERT_EQ(innerPolygonSize, 22); innerPolygon.num_coordinates = innerPolygonSize; std::vector xInner(innerPolygon.num_coordinates); From c0d9aabcf883152f124769c485c58737504471b8 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 15:01:12 +0200 Subject: [PATCH 05/34] GRIDEDIT-2292 Fix windows build --- libs/MeshKernel/src/Mesh2D.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 631e48f1b..f763b0ba2 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -2218,7 +2218,7 @@ void Mesh2D::ReconstructInvalidCellsPolygon() innerBoundaryPolygons.reserve(m_invalidCellPolygons.size()); bool firstElement = true; - for (size_t p = 0; p < polygons.GetNumPolygons(); ++p) + for (UInt p = 0; p < polygons.GetNumPolygons(); ++p) { if (isEnclosingBoundary[p]) { From a455ef6e5db65244ae21671a836ed23b66ee4f73 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 15:08:26 +0200 Subject: [PATCH 06/34] GRIDEDIT-2292 Fixed doxygen spelling warning --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 04a156888..795cf0ac6 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -251,7 +251,7 @@ namespace meshkernel /// @brief Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol) /// @param[in] polygon The polygon where the operation is performed; only boundary segments intersecting the polygon are included - /// @param[in] doAdministrate Indicate that an adminstrate is required. + /// @param[in] doAdministrate Indicate that an administrate is required. /// @return The resulting polygon mesh boundary [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon); @@ -471,7 +471,7 @@ namespace meshkernel /// @brief Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol) /// /// @return a sequence of boundary points, which may be separated by the invalid point, and a matching sequence of Boolean values indicating which - /// polygnal sub-sequence forms a external boudary + /// polygnal sub-sequence forms a external boundary [[nodiscard]] std::tuple, std::vector> GetAllBoundaryPolygons(const std::vector& polygon); /// @brief Ensure that all polynomials are orientated in the ACW direction. From f68995d70a146a60f5fae10bab09e3e6a52356d3 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 15:09:13 +0200 Subject: [PATCH 07/34] GRIDEDIT-2292 Fixed doxygen warning --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 795cf0ac6..233770424 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -251,7 +251,6 @@ namespace meshkernel /// @brief Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol) /// @param[in] polygon The polygon where the operation is performed; only boundary segments intersecting the polygon are included - /// @param[in] doAdministrate Indicate that an administrate is required. /// @return The resulting polygon mesh boundary [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon); From 817075b927dff9386fde2ff223c1ebbaa84effd3 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 15:30:04 +0200 Subject: [PATCH 08/34] GRIDEDIT-2292 Fixed spelling errors --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 233770424..899188f22 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -463,14 +463,14 @@ namespace meshkernel /// @brief Reconstruct the invalid cell polygons /// - /// When constructing the invalid cell polygons, they can be compujted with many smaller polygons. + /// When constructing the invalid cell polygons, they can be computed with many smaller polygons. /// If these smaller polygons form a single patch on the domain, then they need to be combined void ReconstructInvalidCellsPolygon(); /// @brief Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol) /// /// @return a sequence of boundary points, which may be separated by the invalid point, and a matching sequence of Boolean values indicating which - /// polygnal sub-sequence forms a external boundary + /// polygonal sub-sequence forms a external boundary [[nodiscard]] std::tuple, std::vector> GetAllBoundaryPolygons(const std::vector& polygon); /// @brief Ensure that all polynomials are orientated in the ACW direction. From 7b2b9a0e26c2bc23119f6f2655ce26153d463693 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 18:21:45 +0200 Subject: [PATCH 09/34] GRIDEDIT-2293 First attempt at separating the multi polygons --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 3 +- .../include/MeshKernel/Operations.hpp | 6 + libs/MeshKernel/include/MeshKernel/Point.hpp | 8 ++ libs/MeshKernel/src/Mesh2D.cpp | 120 ++++++++++++++--- libs/MeshKernel/src/Operations.cpp | 124 ++++++++++++++++++ .../tests/src/MeshRefinementTests.cpp | 15 +++ 6 files changed, 255 insertions(+), 21 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 899188f22..c8ce399c4 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -448,7 +448,8 @@ namespace meshkernel void WalkBoundaryFromNode(const Polygon& polygon, std::vector& isVisited, UInt& currentNode, - std::vector& meshBoundaryPolygon) const; + std::vector& meshBoundaryPolygon, + std::vector& boundaryPolygonFaceId) const; /// @brief Constructs a polygon or polygons from the meshboundary, by walking through the mesh /// diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 8b9f5d07e..63da7b6f9 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -621,4 +621,10 @@ namespace meshkernel } } + bool isMultiPolygon(std::span boundary); + + std::tuple>, std::vector> splitMultiplePolygons(std::span boundary, std::span elementIds); + + std::vector> splitMultiplePolygons(std::span boundary); + } // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Point.hpp b/libs/MeshKernel/include/MeshKernel/Point.hpp index e40775df0..824c38401 100644 --- a/libs/MeshKernel/include/MeshKernel/Point.hpp +++ b/libs/MeshKernel/include/MeshKernel/Point.hpp @@ -125,6 +125,14 @@ namespace meshkernel return !isInvalid; } + + /// @brief Required by std::set + bool operator<(const Point& other) const + { + if (x != other.x) + return x < other.x; + return y < other.y; + } }; /// @brief Compute the dot product of a point with itself. diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index f763b0ba2..522f38de3 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1523,9 +1523,13 @@ std::tuple, std::vector> Mesh2D::GetAllBoun std::vector isVisited(GetNumEdges(), false); std::vector meshBoundaryPolygon; + std::vector boundaryPolygon; + // Elements connected to the boundary + std::vector boundaryPolygonFaceId; std::vector isEnclosingBoundary; meshBoundaryPolygon.reserve(GetNumNodes()); + boundaryPolygon.reserve(GetNumNodes()); for (UInt e = 0; e < GetNumEdges(); e++) { @@ -1548,49 +1552,123 @@ std::tuple, std::vector> Mesh2D::GetAllBoun } // Start a new polyline - if (!meshBoundaryPolygon.empty()) - { - meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); - } - - const size_t boundaryPolygonStartIndex = meshBoundaryPolygon.size(); + boundaryPolygon.clear(); + boundaryPolygonFaceId.clear(); // Put the current edge on the mesh boundary, mark it as visited - meshBoundaryPolygon.emplace_back(firstNode); - meshBoundaryPolygon.emplace_back(secondNode); + boundaryPolygon.emplace_back(firstNode); + boundaryPolygon.emplace_back(secondNode); + boundaryPolygonFaceId.push_back(m_edgesFaces[e][0]); isVisited[e] = true; // walk the current mesh boundary auto currentNode = secondNodeIndex; - WalkBoundaryFromNode(polygon, isVisited, currentNode, meshBoundaryPolygon); + WalkBoundaryFromNode(polygon, isVisited, currentNode, boundaryPolygon, boundaryPolygonFaceId); - const auto numNodesFirstTail = static_cast(meshBoundaryPolygon.size()); + const auto numNodesFirstTail = static_cast(boundaryPolygon.size()); // if the boundary polygon is not closed if (currentNode != firstNodeIndex) { // Now grow a polyline starting at the other side of the original link L, i.e., the second tail currentNode = firstNodeIndex; - WalkBoundaryFromNode(polygon, isVisited, currentNode, meshBoundaryPolygon); + WalkBoundaryFromNode(polygon, isVisited, currentNode, boundaryPolygon, boundaryPolygonFaceId); } // There is a nonempty second tail: reverse the second tail so that the tails connect and close the polygon. - if (meshBoundaryPolygon.size() > numNodesFirstTail) + if (boundaryPolygon.size() > numNodesFirstTail) { - std::reverse(meshBoundaryPolygon.begin() + numNodesFirstTail, meshBoundaryPolygon.end()); - meshBoundaryPolygon.push_back(meshBoundaryPolygon.front()); + std::reverse(boundaryPolygon.begin() + numNodesFirstTail, boundaryPolygon.end()); + boundaryPolygon.push_back(boundaryPolygon.front()); } - const size_t boundaryPolygonEndIndex = meshBoundaryPolygon.size(); - - std::span currentPolygonSpan(std::span(meshBoundaryPolygon.data() + boundaryPolygonStartIndex, - meshBoundaryPolygon.data() + boundaryPolygonEndIndex)); + std::span currentPolygonSpan(boundaryPolygon); Polygon currentPolygon(currentPolygonSpan, m_projection); - isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[m_edgesFaces[e][0]])); + if (!meshBoundaryPolygon.empty()) + { + meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); + } + + if (isMultiPolygon(currentPolygonSpan)) + { + std::cout << " ++++++++++++++++++++++++++++++++ " << std::endl; + std::cout << " currentPolygonSpan "; + + for (size_t j = 0; j < currentPolygonSpan.size(); ++j) + { + std::cout << "{" << currentPolygonSpan[j].x << ", " << currentPolygonSpan[j].y << "}, "; + } + + std::cout << std::endl; + std::cout << " boundaryPolygonFaceId "; + + for (size_t j = 0; j < boundaryPolygonFaceId.size(); ++j) + { + std::cout << boundaryPolygonFaceId[j] << ", "; + } + + std::cout << std::endl; + + auto [multiBoundaryPolygonNodes, multiBoundaryElementIds] = (splitMultiplePolygons(currentPolygonSpan, boundaryPolygonFaceId)); + // std::vector> multiBoundaryPolygonNodes(splitMultiplePolygons(currentPolygonSpan)); + // size_t faceIndex = 0; + + std::cout << " sub-boundaryPolygonFaceId " << multiBoundaryElementIds.size() << " "; + + for (size_t i = 0; i < multiBoundaryElementIds.size(); ++i) + { + std::cout << multiBoundaryElementIds[i] << ", "; + + for (size_t j = 0; j < multiBoundaryPolygonNodes[i].size(); ++j) + { + std::cout << "{" << multiBoundaryPolygonNodes[i][j].x << ", " << multiBoundaryPolygonNodes[i][j].y << "}, "; + } + + std::cout << std::endl; + } + + std::cout << std::endl; + + for (size_t i = 0; i < multiBoundaryPolygonNodes.size(); ++i) + { + if (i > 0) + { + meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); + } + + meshBoundaryPolygon.insert(meshBoundaryPolygon.end(), multiBoundaryPolygonNodes[i].begin(), multiBoundaryPolygonNodes[i].end()); + + std::span currentPolygonSpan(multiBoundaryPolygonNodes[i]); + Polygon currentPolygon(currentPolygonSpan, m_projection); + + isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[multiBoundaryElementIds[i]])); + // faceIndex += multiBoundaryPolygonNodes[i].size(); + } + } + else + { + meshBoundaryPolygon.insert(meshBoundaryPolygon.end(), boundaryPolygon.begin(), boundaryPolygon.end()); + + std::span currentPolygonSpan(boundaryPolygon); + Polygon currentPolygon(currentPolygonSpan, m_projection); + + isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[boundaryPolygonFaceId[0]])); + } + + // std::cout << " boundaryPolygonFaceId "; + + // for (size_t i = 0; i < boundaryPolygonFaceId.size(); ++i) + // { + // std::cout << boundaryPolygonFaceId[i] << ", "; + // } + + // std::cout << std::endl; } + std::cout << "--------------------------------" << std::endl; + return {meshBoundaryPolygon, isEnclosingBoundary}; } @@ -1797,7 +1875,8 @@ std::vector Mesh2D::RemoveOuterDomainBoundaryPolygon(const st void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, std::vector& isVisited, UInt& currentNode, - std::vector& meshBoundaryPolygon) const + std::vector& meshBoundaryPolygon, + std::vector& boundaryPolygonFaceId) const { UInt e = 0; bool currentNodeInPolygon = false; @@ -1821,6 +1900,7 @@ void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, } currentNode = OtherNodeOfEdge(m_edges[currentEdge], currentNode); + boundaryPolygonFaceId.push_back(m_edgesFaces[currentEdge][0]); e = 0; currentNodeInPolygon = false; diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 1ede27d7c..63331fcf9 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -25,6 +25,7 @@ // //------------------------------------------------------------------------------ +#include #include #include "MeshKernel/Cartesian3DPoint.hpp" @@ -1723,4 +1724,127 @@ namespace meshkernel return (matCoefficients[0] * x[0] + matCoefficients[1] * x[1]) * y[0] + (matCoefficients[2] * x[0] + matCoefficients[3] * x[1]) * y[1]; } + bool isMultiPolygon(std::span boundary) + { + + if (boundary.size() < 4) + { + return false; + } + + std::set uniquePoints; + + for (const Point& p : boundary) + { + if (!uniquePoints.insert(p).second) + { + return true; + } + } + + return false; + } + + std::tuple>, std::vector> splitMultiplePolygons(std::span boundary, std::span elementIds) + { + std::vector> completedPolygons; + std::vector pointStack; + std::vector elementStack; + std::map stackRegistry; // Maps point to its current index in pointStack + + for (const Point& currentPoint : boundary) + { + auto it = stackRegistry.find(currentPoint); + + if (it != stackRegistry.end()) + { + // A duplicate point is found, which means a loop is closed! + size_t loopStartIndex = it->second; + std::vector newPolygon; + + // Extract everything from the start of the loop to the end of the stack + for (size_t i = loopStartIndex; i < pointStack.size(); ++i) + { + newPolygon.push_back(pointStack[i]); + stackRegistry.erase(pointStack[i]); // Clean registry for reused points + } + + // Add the closing point to complete the loop topology + newPolygon.push_back(currentPoint); + completedPolygons.push_back(newPolygon); + + // Shrink the stack back down, discarding the extracted loop + pointStack.resize(loopStartIndex); + elementStack.push_back(elementIds[loopStartIndex]); + } + + // Push the current point onto the active path stack + stackRegistry[currentPoint] = pointStack.size(); + pointStack.push_back(currentPoint); + } + + // Wrap up any remaining points left on the main outer path + if (pointStack.size() > 2) + { + // Ensure it self-closes if the original boundary loop was implicit + if (!(pointStack.front() == pointStack.back())) + { + pointStack.push_back(pointStack.front()); + } + completedPolygons.push_back(pointStack); + } + + return {completedPolygons, elementStack}; + } + + std::vector> splitMultiplePolygons(std::span boundary) + { + std::vector> completedPolygons; + std::vector pointStack; + std::map stackRegistry; // Maps point to its current index in pointStack + + for (const Point& currentPoint : boundary) + { + auto it = stackRegistry.find(currentPoint); + + if (it != stackRegistry.end()) + { + // A duplicate point is found, which means a loop is closed! + size_t loopStartIndex = it->second; + std::vector newPolygon; + + // Extract everything from the start of the loop to the end of the stack + for (size_t i = loopStartIndex; i < pointStack.size(); ++i) + { + newPolygon.push_back(pointStack[i]); + stackRegistry.erase(pointStack[i]); // Clean registry for reused points + } + + // Add the closing point to complete the loop topology + newPolygon.push_back(currentPoint); + completedPolygons.push_back(newPolygon); + + // Shrink the stack back down, discarding the extracted loop + pointStack.resize(loopStartIndex); + } + + // Push the current point onto the active path stack + stackRegistry[currentPoint] = pointStack.size(); + pointStack.push_back(currentPoint); + } + + // Wrap up any remaining points left on the main outer path + if (pointStack.size() > 2) + { + // Ensure it self-closes if the original boundary loop was implicit + if (!(pointStack.front() == pointStack.back())) + { + pointStack.push_back(pointStack.front()); + } + completedPolygons.push_back(pointStack); + } + + return completedPolygons; + } + } // namespace meshkernel diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 948780dd8..6b23ce73d 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -3018,6 +3018,11 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); + auto node71 = mesh.FindNodeCloseToAPoint({75.0, 25.0}, 1.0e-5); + auto node72 = mesh.FindNodeCloseToAPoint({85.0, 25.0}, 1.0e-5); + auto node73 = mesh.FindNodeCloseToAPoint({85.0, 35.0}, 1.0e-5); + auto node74 = mesh.FindNodeCloseToAPoint({75.0, 35.0}, 1.0e-5); + std::vector boundaryNodes; std::vector elementNodes1{{constants::missing::doubleValue, constants::missing::doubleValue}, @@ -3059,6 +3064,13 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon mesh.Node(node63), mesh.Node(node61)}; + std::vector elementNodes7{{constants::missing::doubleValue, constants::missing::doubleValue}, + mesh.Node(node71), + mesh.Node(node72), + mesh.Node(node73), + mesh.Node(node74), + mesh.Node(node71)}; + // Combine all nodes to form a sequence of polygons boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); @@ -3066,6 +3078,7 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); + boundaryNodes.insert(boundaryNodes.end(), elementNodes7.begin(), elementNodes7.end()); Polygons boundaryWithMissingElements(boundaryNodes, Projection::cartesian); @@ -3077,9 +3090,11 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon // Call administrate, to re-compute the face-node and face-edge connectivity. // This should not fill in the holes in the mesh mesh2.Administrate(); + meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "mesh6.vtu"); // Get interior boundary polygon points std::vector innerBoundaryPoints = mesh2.GetInnerBoundaryPolygons(); + return; // The expected number of points, should not include any land boundary points constexpr UInt expectedNumberOfNodes = 22; From 25c8e82d299a29bcbe01f26c664101a1634eeb75 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Jul 2026 18:28:58 +0200 Subject: [PATCH 10/34] GRIDEDIT-2293 Separated multiple overlapping polygons and included Boolean indicating if the polygon was exterior of not --- libs/MeshKernel/src/Mesh2D.cpp | 47 ------------- libs/MeshKernel/src/Operations.cpp | 67 +++++++++---------- .../tests/src/MeshRefinementTests.cpp | 29 +++----- 3 files changed, 44 insertions(+), 99 deletions(-) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 522f38de3..966eee6a6 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1593,43 +1593,8 @@ std::tuple, std::vector> Mesh2D::GetAllBoun if (isMultiPolygon(currentPolygonSpan)) { - std::cout << " ++++++++++++++++++++++++++++++++ " << std::endl; - std::cout << " currentPolygonSpan "; - - for (size_t j = 0; j < currentPolygonSpan.size(); ++j) - { - std::cout << "{" << currentPolygonSpan[j].x << ", " << currentPolygonSpan[j].y << "}, "; - } - - std::cout << std::endl; - std::cout << " boundaryPolygonFaceId "; - - for (size_t j = 0; j < boundaryPolygonFaceId.size(); ++j) - { - std::cout << boundaryPolygonFaceId[j] << ", "; - } - - std::cout << std::endl; auto [multiBoundaryPolygonNodes, multiBoundaryElementIds] = (splitMultiplePolygons(currentPolygonSpan, boundaryPolygonFaceId)); - // std::vector> multiBoundaryPolygonNodes(splitMultiplePolygons(currentPolygonSpan)); - // size_t faceIndex = 0; - - std::cout << " sub-boundaryPolygonFaceId " << multiBoundaryElementIds.size() << " "; - - for (size_t i = 0; i < multiBoundaryElementIds.size(); ++i) - { - std::cout << multiBoundaryElementIds[i] << ", "; - - for (size_t j = 0; j < multiBoundaryPolygonNodes[i].size(); ++j) - { - std::cout << "{" << multiBoundaryPolygonNodes[i][j].x << ", " << multiBoundaryPolygonNodes[i][j].y << "}, "; - } - - std::cout << std::endl; - } - - std::cout << std::endl; for (size_t i = 0; i < multiBoundaryPolygonNodes.size(); ++i) { @@ -1644,7 +1609,6 @@ std::tuple, std::vector> Mesh2D::GetAllBoun Polygon currentPolygon(currentPolygonSpan, m_projection); isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[multiBoundaryElementIds[i]])); - // faceIndex += multiBoundaryPolygonNodes[i].size(); } } else @@ -1656,19 +1620,8 @@ std::tuple, std::vector> Mesh2D::GetAllBoun isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[boundaryPolygonFaceId[0]])); } - - // std::cout << " boundaryPolygonFaceId "; - - // for (size_t i = 0; i < boundaryPolygonFaceId.size(); ++i) - // { - // std::cout << boundaryPolygonFaceId[i] << ", "; - // } - - // std::cout << std::endl; } - std::cout << "--------------------------------" << std::endl; - return {meshBoundaryPolygon, isEnclosingBoundary}; } diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 63331fcf9..2bc5ace76 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -1745,56 +1745,55 @@ namespace meshkernel return false; } - std::tuple>, std::vector> splitMultiplePolygons(std::span boundary, std::span elementIds) + std::tuple>, std::vector> splitMultiplePolygons(std::span boundaryPoints, std::span elementIds) { std::vector> completedPolygons; - std::vector pointStack; - std::vector elementStack; - std::map stackRegistry; // Maps point to its current index in pointStack + std::vector firstElementIds; - for (const Point& currentPoint : boundary) + std::vector> pointFaceStack; + // Maps a point to its current index in the pointFaceStack + std::map activePoints; + + for (size_t i = 0; i < boundaryPoints.size(); ++i) { - auto it = stackRegistry.find(currentPoint); + const Point& current_point = boundaryPoints[i]; - if (it != stackRegistry.end()) + // Check if this node closes a loop with a previously visited point + if (auto it = activePoints.find(current_point); it != activePoints.end()) { - // A duplicate point is found, which means a loop is closed! - size_t loopStartIndex = it->second; - std::vector newPolygon; + size_t loop_start_idx = it->second; - // Extract everything from the start of the loop to the end of the stack - for (size_t i = loopStartIndex; i < pointStack.size(); ++i) + std::vector sub_polygon; + sub_polygon.reserve((pointFaceStack.size() - loop_start_idx) + 1); + + int first_edge_id = pointFaceStack[loop_start_idx].second; + + for (size_t j = loop_start_idx; j < pointFaceStack.size(); ++j) { - newPolygon.push_back(pointStack[i]); - stackRegistry.erase(pointStack[i]); // Clean registry for reused points + sub_polygon.push_back(pointFaceStack[j].first); + activePoints.erase(pointFaceStack[j].first); } - // Add the closing point to complete the loop topology - newPolygon.push_back(currentPoint); - completedPolygons.push_back(newPolygon); + // Explicitly close the polygon + if (!sub_polygon.empty()) + { + sub_polygon.push_back(sub_polygon.front()); + } - // Shrink the stack back down, discarding the extracted loop - pointStack.resize(loopStartIndex); - elementStack.push_back(elementIds[loopStartIndex]); + completedPolygons.push_back(sub_polygon); + firstElementIds.push_back(first_edge_id); + + // Pop the loop off the stack + pointFaceStack.resize(loop_start_idx); } - // Push the current point onto the active path stack - stackRegistry[currentPoint] = pointStack.size(); - pointStack.push_back(currentPoint); - } + int current_edge = (i < elementIds.size()) ? elementIds[i] : -1; - // Wrap up any remaining points left on the main outer path - if (pointStack.size() > 2) - { - // Ensure it self-closes if the original boundary loop was implicit - if (!(pointStack.front() == pointStack.back())) - { - pointStack.push_back(pointStack.front()); - } - completedPolygons.push_back(pointStack); + activePoints[current_point] = pointFaceStack.size(); + pointFaceStack.emplace_back(current_point, current_edge); } - return {completedPolygons, elementStack}; + return {completedPolygons, firstElementIds}; } std::vector> splitMultiplePolygons(std::span boundary) diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 6b23ce73d..c71641600 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -3018,11 +3018,6 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); - auto node71 = mesh.FindNodeCloseToAPoint({75.0, 25.0}, 1.0e-5); - auto node72 = mesh.FindNodeCloseToAPoint({85.0, 25.0}, 1.0e-5); - auto node73 = mesh.FindNodeCloseToAPoint({85.0, 35.0}, 1.0e-5); - auto node74 = mesh.FindNodeCloseToAPoint({75.0, 35.0}, 1.0e-5); - std::vector boundaryNodes; std::vector elementNodes1{{constants::missing::doubleValue, constants::missing::doubleValue}, @@ -3064,13 +3059,6 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon mesh.Node(node63), mesh.Node(node61)}; - std::vector elementNodes7{{constants::missing::doubleValue, constants::missing::doubleValue}, - mesh.Node(node71), - mesh.Node(node72), - mesh.Node(node73), - mesh.Node(node74), - mesh.Node(node71)}; - // Combine all nodes to form a sequence of polygons boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); @@ -3078,7 +3066,6 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); - boundaryNodes.insert(boundaryNodes.end(), elementNodes7.begin(), elementNodes7.end()); Polygons boundaryWithMissingElements(boundaryNodes, Projection::cartesian); @@ -3090,26 +3077,32 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon // Call administrate, to re-compute the face-node and face-edge connectivity. // This should not fill in the holes in the mesh mesh2.Administrate(); - meshkernel::SaveVtk(mesh.Nodes(), mesh.m_facesNodes, "mesh6.vtu"); // Get interior boundary polygon points std::vector innerBoundaryPoints = mesh2.GetInnerBoundaryPolygons(); - return; // The expected number of points, should not include any land boundary points - constexpr UInt expectedNumberOfNodes = 22; + constexpr UInt expectedNumberOfNodes = 26; ASSERT_EQ(expectedNumberOfNodes, innerBoundaryPoints.size()); // The edge of one of the deleted elements lies on the boundary, so will be not be part of the // interior set of polygons - std::vector expectedXPoints{80.0, 85.0, 95.0, 100.0, 105.0, 95.0, 95.0, 85.0, 85.0, 75.0, 80.0, + std::vector expectedXPoints{95.0, 100.0, 105.0, 95.0, + constants::missing::doubleValue, + 85.0, 95.0, 95.0, 85.0, 85.0, + constants::missing::doubleValue, + 80.0, 85.0, 75.0, 80.0, constants::missing::doubleValue, 120.0, 125.0, 115.0, 120.0, constants::missing::doubleValue, 125.0, 125.0, 135.0, 135.0, 125.0}; - std::vector expectedYPoints{0.0, 15.0, 15.0, 0.0, 15.0, 15.0, 25.0, 25.0, 15.0, 15.0, 0.0, + std::vector expectedYPoints{15.0, 0.0, 15.0, 15.0, + constants::missing::doubleValue, + 15.0, 15.0, 25.0, 25.0, 15.0, + constants::missing::doubleValue, + 0.0, 15.0, 15.0, 0.0, constants::missing::doubleValue, 0.0, 15.0, 15.0, 0.0, constants::missing::doubleValue, From 8968d1b4cd145b0c521aa68e116a159e13b0bdbe Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 13 Jul 2026 11:15:59 +0200 Subject: [PATCH 11/34] GRIDEDIT-2293 Added comment describing the limitation of the alogorithm and fixed the unit test --- .../include/MeshKernel/Operations.hpp | 8 +- libs/MeshKernel/src/Operations.cpp | 77 ++++--------------- .../tests/src/InvalidCellsPolygonsTests.cpp | 46 +++++++---- 3 files changed, 51 insertions(+), 80 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 63da7b6f9..9ec570886 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -621,10 +621,14 @@ namespace meshkernel } } + /// @brief Determine is a polygon is comprised of multiple polygons intersecting at the control points only + /// + /// @note This currently works only for polygons that intersect at the control points bool isMultiPolygon(std::span boundary); + /// @brief Split a single polygon line comprised of multiple polygons into separate polygons. + /// + /// @note This currently works only for polygons that intersect at the control points std::tuple>, std::vector> splitMultiplePolygons(std::span boundary, std::span elementIds); - std::vector> splitMultiplePolygons(std::span boundary); - } // namespace meshkernel diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 2bc5ace76..0fb145962 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -1761,89 +1761,38 @@ namespace meshkernel // Check if this node closes a loop with a previously visited point if (auto it = activePoints.find(current_point); it != activePoints.end()) { - size_t loop_start_idx = it->second; + size_t loopStartIndex = it->second; - std::vector sub_polygon; - sub_polygon.reserve((pointFaceStack.size() - loop_start_idx) + 1); + std::vector subPolygon; + subPolygon.reserve((pointFaceStack.size() - loopStartIndex) + 1); - int first_edge_id = pointFaceStack[loop_start_idx].second; + int first_edge_id = pointFaceStack[loopStartIndex].second; - for (size_t j = loop_start_idx; j < pointFaceStack.size(); ++j) + for (size_t j = loopStartIndex; j < pointFaceStack.size(); ++j) { - sub_polygon.push_back(pointFaceStack[j].first); + subPolygon.push_back(pointFaceStack[j].first); activePoints.erase(pointFaceStack[j].first); } - // Explicitly close the polygon - if (!sub_polygon.empty()) + // close the polygon + if (!subPolygon.empty()) { - sub_polygon.push_back(sub_polygon.front()); + subPolygon.push_back(subPolygon.front()); } - completedPolygons.push_back(sub_polygon); + completedPolygons.push_back(subPolygon); firstElementIds.push_back(first_edge_id); - // Pop the loop off the stack - pointFaceStack.resize(loop_start_idx); + pointFaceStack.resize(loopStartIndex); } - int current_edge = (i < elementIds.size()) ? elementIds[i] : -1; + int currentEdge = (i < elementIds.size()) ? elementIds[i] : -1; activePoints[current_point] = pointFaceStack.size(); - pointFaceStack.emplace_back(current_point, current_edge); + pointFaceStack.emplace_back(current_point, currentEdge); } return {completedPolygons, firstElementIds}; } - std::vector> splitMultiplePolygons(std::span boundary) - { - std::vector> completedPolygons; - std::vector pointStack; - std::map stackRegistry; // Maps point to its current index in pointStack - - for (const Point& currentPoint : boundary) - { - auto it = stackRegistry.find(currentPoint); - - if (it != stackRegistry.end()) - { - // A duplicate point is found, which means a loop is closed! - size_t loopStartIndex = it->second; - std::vector newPolygon; - - // Extract everything from the start of the loop to the end of the stack - for (size_t i = loopStartIndex; i < pointStack.size(); ++i) - { - newPolygon.push_back(pointStack[i]); - stackRegistry.erase(pointStack[i]); // Clean registry for reused points - } - - // Add the closing point to complete the loop topology - newPolygon.push_back(currentPoint); - completedPolygons.push_back(newPolygon); - - // Shrink the stack back down, discarding the extracted loop - pointStack.resize(loopStartIndex); - } - - // Push the current point onto the active path stack - stackRegistry[currentPoint] = pointStack.size(); - pointStack.push_back(currentPoint); - } - - // Wrap up any remaining points left on the main outer path - if (pointStack.size() > 2) - { - // Ensure it self-closes if the original boundary loop was implicit - if (!(pointStack.front() == pointStack.back())) - { - pointStack.push_back(pointStack.front()); - } - completedPolygons.push_back(pointStack); - } - - return completedPolygons; - } - } // namespace meshkernel diff --git a/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp index 97d77b1fd..b00ad89fa 100644 --- a/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp +++ b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp @@ -82,18 +82,28 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) //-------------------------------- - std::vector invalidCellsX{80.0, 85.0, 75.0, 80.0, meshkernel::constants::missing::doubleValue, - 180.0, 200.0, 200.0, 180.0, 180.0, meshkernel::constants::missing::doubleValue, - 125.0, 135.0, 135.0, 125.0, 125.0, meshkernel::constants::missing::doubleValue, - 85.0, 95.0, 95.0, 85.0, 85.0, meshkernel::constants::missing::doubleValue, - 100.0, 105.0, 95.0, 100.0, meshkernel::constants::missing::doubleValue, + std::vector invalidCellsX{80.0, 85.0, 75.0, 80.0, + meshkernel::constants::missing::doubleValue, + 180.0, 200.0, 200.0, 180.0, 180.0, + meshkernel::constants::missing::doubleValue, + 125.0, 135.0, 135.0, 125.0, 125.0, + meshkernel::constants::missing::doubleValue, + 85.0, 95.0, 95.0, 85.0, 85.0, + meshkernel::constants::missing::doubleValue, + 100.0, 105.0, 95.0, 100.0, + meshkernel::constants::missing::doubleValue, 120.0, 125.0, 115.0, 120.0}; - std::vector invalidCellsY{0.0, 15.0, 15.0, 0.0, meshkernel::constants::missing::doubleValue, - 140.0, 140.0, 160.0, 160.0, 140.0, meshkernel::constants::missing::doubleValue, - 125.0, 125.0, 135.0, 135.0, 125.0, meshkernel::constants::missing::doubleValue, - 15.0, 15.0, 25.0, 25.0, 15.0, meshkernel::constants::missing::doubleValue, - 0.0, 15.0, 15.0, 0.0, meshkernel::constants::missing::doubleValue, + std::vector invalidCellsY{0.0, 15.0, 15.0, 0.0, + meshkernel::constants::missing::doubleValue, + 140.0, 140.0, 160.0, 160.0, 140.0, + meshkernel::constants::missing::doubleValue, + 125.0, 125.0, 135.0, 135.0, 125.0, + meshkernel::constants::missing::doubleValue, + 15.0, 15.0, 25.0, 25.0, 15.0, + meshkernel::constants::missing::doubleValue, + 0.0, 15.0, 15.0, 0.0, + meshkernel::constants::missing::doubleValue, 0.0, 15.0, 15.0, 0.0}; meshkernelapi::GeometryList invalidCells; @@ -126,7 +136,7 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) ASSERT_EQ(mesh2d.num_nodes, 3130); ASSERT_EQ(mesh2d.num_edges, 6383); - ASSERT_EQ(mesh2d.num_faces, 3248); + ASSERT_EQ(mesh2d.num_faces, 3232); int whichMeshkernelId = -1; bool isUndone = false; @@ -191,13 +201,21 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) //----------------------- // Check the inner boundary polygon are correct - std::vector expectedInnerX{80.0, 85.0, 95.0, 100.0, 105.0, 95.0, 95.0, 85.0, 85.0, 75.0, 80.0, + std::vector expectedInnerX{95.0, 100.0, 105.0, 95.0, + meshkernel::constants::missing::doubleValue, + 85.0, 95.0, 95.0, 85.0, 85.0, + meshkernel::constants::missing::doubleValue, + 80.0, 85.0, 75.0, 80.0, meshkernel::constants::missing::doubleValue, 120.0, 125.0, 115.0, 120.0, meshkernel::constants::missing::doubleValue, 125.0, 125.0, 135.0, 135.0, 125.0}; - std::vector expectedInnerY{0.0, 15.0, 15.0, 0.0, 15.0, 15.0, 25.0, 25.0, 15.0, 15.0, 0.0, + std::vector expectedInnerY{15.0, 0.0, 15.0, 15.0, + meshkernel::constants::missing::doubleValue, + 15.0, 15.0, 25.0, 25.0, 15.0, + meshkernel::constants::missing::doubleValue, + 0.0, 15.0, 15.0, 0.0, meshkernel::constants::missing::doubleValue, 0.0, 15.0, 15.0, 0.0, meshkernel::constants::missing::doubleValue, @@ -209,7 +227,7 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) errorCode = meshkernelapi::mkernel_mesh2d_get_mesh_inner_boundaries_as_polygons_dimension(meshKernelId, innerPolygonSize); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); - ASSERT_EQ(innerPolygonSize, 22); + ASSERT_EQ(innerPolygonSize, 26); innerPolygon.num_coordinates = innerPolygonSize; std::vector xInner(innerPolygon.num_coordinates); From 12b9b7df090fda33c6a186e0e2ef909dd130af69 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Wed, 29 Jul 2026 19:15:15 +0200 Subject: [PATCH 12/34] GRIDEDIT-2293 Passing to other computer, another attempt at extracting the domain boundaries --- libs/MeshKernel/CMakeLists.txt | 2 + .../MeshKernel/MeshBoundaryExtractor.hpp | 73 ++ .../include/MeshKernel/Operations.hpp | 5 + libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 199 ++++ libs/MeshKernel/src/Operations.cpp | 52 +- .../tests/src/MeshRefinementTests.cpp | 910 ++++++++++++++++++ 6 files changed, 1239 insertions(+), 2 deletions(-) create mode 100644 libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp create mode 100644 libs/MeshKernel/src/MeshBoundaryExtractor.cpp diff --git a/libs/MeshKernel/CMakeLists.txt b/libs/MeshKernel/CMakeLists.txt index 2fb3902ad..e2dd97211 100644 --- a/libs/MeshKernel/CMakeLists.txt +++ b/libs/MeshKernel/CMakeLists.txt @@ -45,6 +45,7 @@ set( ${SRC_DIR}/LandBoundaries.cpp ${SRC_DIR}/LandBoundary.cpp ${SRC_DIR}/Mesh.cpp + ${SRC_DIR}/MeshBoundaryExtractor.cpp ${SRC_DIR}/MeshEdgeLength.cpp ${SRC_DIR}/Mesh1D.cpp ${SRC_DIR}/Mesh2D.cpp @@ -181,6 +182,7 @@ set( ${DOMAIN_INC_DIR}/LandBoundaries.hpp ${DOMAIN_INC_DIR}/LandBoundary.hpp ${DOMAIN_INC_DIR}/Mesh.hpp + ${DOMAIN_INC_DIR}/MeshBoundaryExtractor.hpp ${DOMAIN_INC_DIR}/MeshEdgeLength.hpp ${DOMAIN_INC_DIR}/Mesh1D.hpp ${DOMAIN_INC_DIR}/Mesh2D.hpp diff --git a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp new file mode 100644 index 000000000..f3b40ce1d --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp @@ -0,0 +1,73 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2026. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program 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 General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include +#include + +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/Point.hpp" + +namespace meshkernel +{ + + /// @brief Extract the boundary polygon from the mesh + class MeshBoundaryExtractor + { + public: + /// @brief Extract all boundaries, concatinated as a single sequence of points, separated by an invalid point + static std::vector ExtractAll(const Mesh2D& mesh); + + /// @brief Extract all boundaries keeping them separated and + static std::tuple>, std::vector> Extract(const Mesh2D& mesh); + + private: + /// @brief Temporary struct, used when computing the boundaries + struct BoundaryEdge + { + UInt edgeId; + UInt neighbourNode; + UInt leftFace; // Store face mapping on the edge for easy retrieval during loop trace + double angle; // Angle of the edge pointing away from the pivot node + }; + + /// @brief Ensure the angle lies between 0 .. 2 pi. + static double NormalizeAngle(double angle); + + /// @brief Find boundary loops + /// + /// Any boundary loops found may need to be processed further as they may themselves contain sub-loops + static void FindBoundaryLoops(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::vector>& allLoops, + std::vector>& allTouchedFaces); + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 9ec570886..77a025d87 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -631,4 +631,9 @@ namespace meshkernel /// @note This currently works only for polygons that intersect at the control points std::tuple>, std::vector> splitMultiplePolygons(std::span boundary, std::span elementIds); + /// @brief Split a single polygon line comprised of multiple polygons into separate polygons. + /// + /// @note This currently works only for polygons that intersect at the control points + std::vector> splitMultiplePolygons(std::span boundary); + } // namespace meshkernel diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp new file mode 100644 index 000000000..69c15cfd0 --- /dev/null +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -0,0 +1,199 @@ +#include "MeshKernel/MeshBoundaryExtractor.hpp" + +#include + +#include "MeshKernel/Operations.hpp" + +std::vector meshkernel::MeshBoundaryExtractor::ExtractAll(const Mesh2D& mesh) +{ + + auto [boundarySequences, isExterior] = Extract(mesh); + std::vector allPoints; + + bool isFirst = true; + + for (const std::vector& loop : boundarySequences) + { + if (!isFirst) + { + allPoints.push_back({constants::missing::doubleValue, constants::missing::doubleValue}); + } + + allPoints.insert(allPoints.end(), loop.begin(), loop.end()); + isFirst = false; + } + + return allPoints; +} + +std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh) +{ + + std::vector> allBoundaryLoops; + std::vector> allTouchedFaces; + + std::vector> separatedBoundaryLoops; + std::vector isExterior; + + const std::vector& meshNodes(mesh.Nodes()); + + FindBoundaryLoops(meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryLoops, allTouchedFaces); + + for (size_t i = 0; i < allBoundaryLoops.size(); ++i) + { + if (isMultiPolygon(allBoundaryLoops[i])) + { + auto [individualBoundaryPolygons, firstElement] = splitMultiplePolygons(allBoundaryLoops[i], allTouchedFaces[i]); + + for (size_t i = 0; i < individualBoundaryPolygons.size(); ++i) + { + Point centre = mesh.m_facesMassCenters[firstElement[i]]; + isExterior.push_back(IsPointInPolygonNodes(centre, individualBoundaryPolygons[i], mesh.m_projection)); + separatedBoundaryLoops.push_back(std::move(individualBoundaryPolygons[i])); + } + } + else + { + separatedBoundaryLoops.push_back(std::move(allBoundaryLoops[i])); + } + } + + return {separatedBoundaryLoops, isExterior}; +} + +double meshkernel::MeshBoundaryExtractor::NormalizeAngle(double angle) +{ + while (angle < 0) + { + angle += 2.0 * std::numbers::pi; + } + + while (angle >= 2.0 * std::numbers::pi) + { + angle -= 2.0 * std::numbers::pi; + } + + return angle; +} + +void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::vector>& allLoops, + std::vector>& allTouchedFaces) +{ + allLoops.clear(); + allTouchedFaces.clear(); + + // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle + std::unordered_map> boundaryAdjacency; + + // May be better to use meshkernel::Boolean + std::vector edgeVisited(edges.size(), false); + + // Collect all boundary edges and compute the angle + for (UInt count = 0; count < edges.size(); ++count) + { + + if (edgesFaces[count][1] == constants::missing::uintValue) + { + double dx = nodes[edges[count].second].x - nodes[edges[count].first].x; + double dy = nodes[edges[count].second].y - nodes[edges[count].first].y; + + // Compute angle for "pinch point" sorting + // Pinch points are when the boundary polygon intersects itself, the points will have identical values. + double angleStartToEnd = NormalizeAngle(std::atan2(dy, dx)); + double angleEndToStart = NormalizeAngle(std::atan2(-dy, -dx)); + + // Note: edgesFaces[count][0] is the internal mesh face valid for both directions of boundary traversal + boundaryAdjacency[edges[count].first].push_back({count, edges[count].second, edgesFaces[count][0], angleStartToEnd}); + boundaryAdjacency[edges[count].second].push_back({count, edges[count].first, edgesFaces[count][0], angleEndToStart}); + } + } + + auto boundaryAngleLessThan = [](const BoundaryEdge& a, const BoundaryEdge& b) + { return a.angle < b.angle; }; + + // Sort outgoing edges anti-clockwise + for (auto& [node_id, connected_edges] : boundaryAdjacency) + { + std::sort(connected_edges.begin(), connected_edges.end(), boundaryAngleLessThan); + } + + // Trace boundary polygons + for (UInt count = 0; count < edges.size(); ++count) + { + if (edgesFaces[count][1] != constants::missing::uintValue || edgeVisited[count]) + { + continue; + } + + std::vector currentNodes; + std::vector currentFaces; + + UInt prevNodeIndex = edges[count].first; + UInt currentNodeIndex = edges[count].second; + + currentNodes.push_back(nodes[prevNodeIndex]); + currentFaces.push_back(edgesFaces[count][0]); // First face touched by the initial edge + edgeVisited[count] = true; // Mark edge as having been visited + + double dx = nodes[currentNodeIndex].x - nodes[prevNodeIndex].x; + double dy = nodes[currentNodeIndex].y - nodes[prevNodeIndex].y; + double incomingAngle = NormalizeAngle(std::atan2(dy, dx)); + + // Loop until we find the start node, making a closed boundary polygon + while (currentNodeIndex != edges[count].first) + { + currentNodes.push_back(nodes[currentNodeIndex]); + + const std::vector& boundaryEdges = boundaryAdjacency[currentNodeIndex]; + + // The smallest visited angle. + // All edge angles must be in interval [-2pi, 2pi], + double deltaAngle = 3.0 * std::numbers::pi; + UInt chosenEdgeIndex = constants::missing::uintValue; + + // Find the boundary edge that tracks closest clockwise to keep empty space on the right + for (size_t e = 0; e < boundaryEdges.size(); ++e) + { + if (edgeVisited[boundaryEdges[e].edgeId]) + { + continue; + } + + double delta = NormalizeAngle(incomingAngle - boundaryEdges[e].angle); + + if (delta < deltaAngle) + { + deltaAngle = delta; + chosenEdgeIndex = e; + } + } + + // No unvisited edges were found + // So eigher the boundary loop was completed or a dead-end reached. + if (chosenEdgeIndex == constants::missing::uintValue) + { + break; + } + + const BoundaryEdge& chosenEdge = boundaryEdges[chosenEdgeIndex]; + edgeVisited[chosenEdge.edgeId] = true; + + // Record the face touched by this next step in the loop + currentFaces.push_back(chosenEdge.leftFace); + + prevNodeIndex = currentNodeIndex; + currentNodeIndex = chosenEdge.neighbourNode; + incomingAngle = chosenEdge.angle; + } + + if (currentNodes.size() >= 3) + { + currentNodes.push_back(currentNodes.front()); + allLoops.push_back(std::move(currentNodes)); + allTouchedFaces.push_back(std::move(currentFaces)); + } + } +} diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 0fb145962..6bee46021 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -1766,7 +1766,7 @@ namespace meshkernel std::vector subPolygon; subPolygon.reserve((pointFaceStack.size() - loopStartIndex) + 1); - int first_edge_id = pointFaceStack[loopStartIndex].second; + int firstEdgeId = pointFaceStack[loopStartIndex].second; for (size_t j = loopStartIndex; j < pointFaceStack.size(); ++j) { @@ -1781,7 +1781,7 @@ namespace meshkernel } completedPolygons.push_back(subPolygon); - firstElementIds.push_back(first_edge_id); + firstElementIds.push_back(firstEdgeId); pointFaceStack.resize(loopStartIndex); } @@ -1795,4 +1795,52 @@ namespace meshkernel return {completedPolygons, firstElementIds}; } + std::vector> splitMultiplePolygons(std::span boundaryPoints) + { + // std::vector dummyConnectedFaces (boundaryPoints.size () - 1); + // std::ranges::iota (dummyConnectedFaces, 0); + // auto [completedPolygons, dummyFirstTouchedFaces] = splitMultiplePolygons (boundaryPoints, dummyConnectedFaces); + + std::vector> completedPolygons; + + std::vector pointFaceStack; + // Maps a point to its current index in the pointFaceStack + std::map activePoints; + + for (size_t i = 0; i < boundaryPoints.size(); ++i) + { + const Point& current_point = boundaryPoints[i]; + + // Check if this node closes a loop with a previously visited point + if (auto it = activePoints.find(current_point); it != activePoints.end()) + { + size_t loopStartIndex = it->second; + + std::vector subPolygon; + subPolygon.reserve((pointFaceStack.size() - loopStartIndex) + 1); + + for (size_t j = loopStartIndex; j < pointFaceStack.size(); ++j) + { + subPolygon.push_back(pointFaceStack[j]); + activePoints.erase(pointFaceStack[j]); + } + + // close the polygon + if (!subPolygon.empty()) + { + subPolygon.push_back(subPolygon.front()); + } + + completedPolygons.push_back(subPolygon); + + pointFaceStack.resize(loopStartIndex); + } + + activePoints[current_point] = pointFaceStack.size(); + pointFaceStack.emplace_back(current_point); + } + + return completedPolygons; + } + } // namespace meshkernel diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index c71641600..df618741d 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -36,6 +36,7 @@ #include "MeshKernel/CasulliDeRefinement.hpp" #include "MeshKernel/CasulliRefinement.hpp" #include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/MeshBoundaryExtractor.hpp" #include "MeshKernel/MeshEdgeLength.hpp" #include "MeshKernel/MeshFaceCenters.hpp" #include "MeshKernel/MeshRefinement.hpp" @@ -2962,6 +2963,754 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa } } +#include +#include +#include +#include + +struct Edge2 +{ + meshkernel::UInt id; + meshkernel::UInt startNode; + meshkernel::UInt endNode; + meshkernel::UInt leftFace; + meshkernel::UInt rightFace; +}; + +struct BoundaryEdge +{ + meshkernel::UInt edgeId; + meshkernel::UInt neighbourNode; + meshkernel::UInt leftFace; // Store face mapping on the edge for easy retrieval during loop trace + double angle; // Angle of the edge pointing away from the pivot node +}; + +double normalizeAngle(double angle) +{ + while (angle < 0) + { + angle += 2.0 * std::numbers::pi; + } + + while (angle >= 2.0 * std::numbers::pi) + { + angle -= 2.0 * std::numbers::pi; + } + + return angle; +} + +namespace meshkernel +{ + static const UInt nullValue = constants::missing::uintValue; + +} + +// Main function to extract loops +std::vector> findBoundaryLoops(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces [[maybe_unused]]) +{ + + // // Step 1 & 2: Build a radially sorted boundary adjacency list + // // Node ID -> list of connected boundary edges + // std::unordered_map> boundaryAdjacency; + // std::unordered_map edgeVisited; + + // for (meshkernel::UInt edgeId = 0; const auto& edge : edges) + // { + // if (edgesFaces[edgeId][1] == nullValue) + // { + // edgeVisited[edgeId] = false; + + // // Calculate geometric angles for pinch-point sorting + // double dx = nodes[edge.second].x - nodes[edge.first].x; + // double dy = nodes[edge.second].y - nodes[edge.first].y; + + // double angle_start_to_end = std::atan2(dy, dx); + // double angle_end_to_start = std::atan2(-dy, -dx); + + // boundaryAdjacency[edge.first].push_back({edgeId, edge.second, angle_start_to_end}); + // boundaryAdjacency[edge.second].push_back({edgeId, edge.first, angle_end_to_start}); + // } + + // ++edgeId; + // } + + // // Sort outgoing boundary edges at each node counter-clockwise + // for (auto& [node_id, connected_edges] : boundaryAdjacency) + // { + // std::sort(connected_edges.begin(), connected_edges.end(), [](const BoundaryEdge& a, const BoundaryEdge& b) + // { return a.angle < b.angle; }); + // } + + // std::vector> allLoops; + + // // Step 3: Trace loops + // for (meshkernel::UInt edgeId = 0; const auto& edge : edges) + // { + // if (edgesFaces[edgeId][1] != nullValue || edgeVisited[edgeId]) + // continue; + + // std::vector currentLoop; + // meshkernel::UInt currentNode = edge.first; + // meshkernel::UInt startNode = currentNode; + + // while (true) + // { + // currentLoop.push_back(currentNode); + + // // Find the next unvisited boundary edge radiating from currentNode + // auto& options = boundaryAdjacency[currentNode]; + // meshkernel::UInt nextNode = nullValue; + // meshkernel::UInt chosenEdgeId = nullValue; + + // for (const auto& option : options) + // { + // if (!edgeVisited[option.edgeId]) + // { + // chosenEdgeId = option.edgeId; + // nextNode = option.neighbourNode; + // break; + // } + // } + + // if (chosenEdgeId == nullValue) + // break; // Loop fully completed or dead end reached + + // edgeVisited[chosenEdgeId] = true; + // currentNode = nextNode; + + // if (currentNode == startNode) + // { + // break; // Loop closed successfully + // } + // } + + // if (!currentLoop.empty()) + // { + // currentLoop.push_back(currentLoop.front()); + + // allLoops.push_back(currentLoop); + // } + + // ++edgeId; + // } + + //-------------------------------- + + // Node ID -> list of connected boundary edges + std::unordered_map> boundaryAdjacency; + std::unordered_map edgeVisited; + + // 1. Collect all boundary edges (where rightFace is nullValue) + for (const auto& edge : edges) + { + if (edge.rightFace == constants::missing::uintValue) + { + edgeVisited[edge.id] = false; + + double dx = nodes[edge.endNode].x - nodes[edge.startNode].x; + double dy = nodes[edge.endNode].y - nodes[edge.startNode].y; + + // Angular direction of the edge + double angle_start_to_end = normalizeAngle(std::atan2(dy, dx)); + double angle_end_to_start = normalizeAngle(std::atan2(-dy, -dx)); + + boundaryAdjacency[edge.startNode].push_back({edge.id, edge.endNode, 0, angle_start_to_end}); + boundaryAdjacency[edge.endNode].push_back({edge.id, edge.startNode, 0, angle_end_to_start}); + } + } + + // 2. Sort outgoing edges counter-clockwise around every node + for (auto& [node_id, connected_edges] : boundaryAdjacency) + { + std::sort(connected_edges.begin(), connected_edges.end(), [](const BoundaryEdge& a, const BoundaryEdge& b) + { return a.angle < b.angle; }); + } + + std::vector> allLoops; + // 3. Trace loops tracking the incoming direction + for (const auto& edge : edges) + { + if (edge.rightFace != constants::missing::uintValue || edgeVisited[edge.id]) + continue; + + std::vector currentLoop; + + // Start tracing this edge from startNode to endNode + meshkernel::UInt prevNode = edge.startNode; + meshkernel::UInt currentNode = edge.endNode; + + currentLoop.push_back(nodes[prevNode]); + edgeVisited[edge.id] = true; + + // Calculate initial incoming angle into the second node + double dx = nodes[currentNode].x - nodes[prevNode].x; + double dy = nodes[currentNode].y - nodes[prevNode].y; + double incomingAngle = normalizeAngle(std::atan2(dy, dx)); + + while (currentNode != edge.startNode) + { + currentLoop.push_back(nodes[currentNode]); + + const auto& options = boundaryAdjacency[currentNode]; + double bestDelta = 3.0 * std::numbers::pi; // Higher than max possible delta (2*PI) + size_t chosenOptionIdx = static_cast(-1); + + // Find the boundary edge that tracks closest clockwise to keep empty space on the right + for (size_t i = 0; i < options.size(); ++i) + { + if (edgeVisited[options[i].edgeId]) + continue; + + double delta = incomingAngle - options[i].angle; + if (delta < 0) + delta += 2.0 * std::numbers::pi; + + if (delta < bestDelta) + { + bestDelta = delta; + chosenOptionIdx = i; + } + } + + // If no unvisited edges are found, break out safely + if (chosenOptionIdx == static_cast(-1)) + break; + + const auto& chosen = options[chosenOptionIdx]; + edgeVisited[chosen.edgeId] = true; + + // Step forward + prevNode = currentNode; + currentNode = chosen.neighbourNode; + incomingAngle = chosen.angle; + } + if (currentLoop.size() >= 3) + { + currentLoop.push_back(currentLoop.front()); + allLoops.push_back(currentLoop); + } + } + + return allLoops; +} + +// void findBoundaryLoops(const std::vector& nodes, +// const std::vector& edges, +// std::vector>& allLoops, // Out: The boundary node sequences +// std::vector& firstTouchedFaces) // Out: The associated face IDs +// { +// // Clear outputs to ensure clean state +// allLoops.clear(); +// firstTouchedFaces.clear(); + +// std::unordered_map> boundaryAdjacency; +// std::unordered_map edgeVisited; + +// // 1. Collect all boundary edges +// for (const auto& edge : edges) +// { +// if (edge.rightFace == nullValue) +// { +// edgeVisited[edge.id] = false; + +// double dx = nodes[edge.endNode].x - nodes[edge.startNode].x; +// double dy = nodes[edge.endNode].y - nodes[edge.startNode].y; + +// double angle_start_to_end = normalizeAngle(std::atan2(dy, dx)); +// double angle_end_to_start = normalizeAngle(std::atan2(-dy, -dx)); + +// boundaryAdjacency[edge.startNode].push_back({edge.id, edge.endNode, 0, angle_start_to_end}); +// boundaryAdjacency[edge.endNode].push_back({edge.id, edge.startNode, 0, angle_end_to_start}); +// } +// } + +// // 2. Sort outgoing edges counter-clockwise +// for (auto& [node_id, connected_edges] : boundaryAdjacency) +// { +// std::sort(connected_edges.begin(), connected_edges.end(), [](const BoundaryEdge& a, const BoundaryEdge& b) +// { return a.angle < b.angle; }); +// } + +// // 3. Trace loops +// for (const auto& edge : edges) +// { +// if (edge.rightFace != nullValue || edgeVisited[edge.id]) +// continue; + +// // Capture the face ID of the first edge that starts this loop +// meshkernel::UInt associatedFaceId = edge.leftFace; + +// std::vector currentNodes; +// meshkernel::UInt prevNode = edge.startNode; +// meshkernel::UInt currentNode = edge.endNode; + +// currentNodes.push_back(prevNode); +// edgeVisited[edge.id] = true; + +// double dx = nodes[currentNode].x - nodes[prevNode].x; +// double dy = nodes[currentNode].y - nodes[prevNode].y; +// double incomingAngle = normalizeAngle(std::atan2(dy, dx)); + +// while (currentNode != edge.startNode) +// { +// currentNodes.push_back(currentNode); + +// const auto& options = boundaryAdjacency[currentNode]; +// double bestDelta = 3.0 * std::numbers::pi; +// size_t chosenOptionIdx = static_cast(-1); + +// for (size_t i = 0; i < options.size(); ++i) +// { +// if (edgeVisited[options[i].edgeId]) +// continue; + +// double delta = incomingAngle - options[i].angle; +// if (delta < 0) +// delta += 2.0 * std::numbers::pi; + +// if (delta < bestDelta) +// { +// bestDelta = delta; +// chosenOptionIdx = i; +// } +// } + +// if (chosenOptionIdx == static_cast(-1)) +// break; + +// const auto& chosen = options[chosenOptionIdx]; +// edgeVisited[chosen.edgeId] = true; + +// prevNode = currentNode; +// currentNode = chosen.neighbourNode; +// incomingAngle = chosen.angle; +// } + +// if (currentNodes.size() >= 3) +// { +// allLoops.push_back(std::move(currentNodes)); +// firstTouchedFaces.push_back(associatedFaceId); +// } +// } +// } + +void findBoundaryLoops(const std::vector& nodes, + const std::vector& edges, + std::vector>& allLoops, // Out: The boundary node sequences + std::vector>& allTouchedFaces) // Out: All faces touched per loop in step order +{ + // Clear outputs to ensure clean state + allLoops.clear(); + allTouchedFaces.clear(); + + // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle + std::unordered_map> boundaryAdjacency; + + // May be better to use meshkernel::Boolean + std::vector edgeVisited(edges.size(), false); + + // Collect all boundary edges and compute the angle + for (const auto& edge : edges) + { + if (edge.rightFace == constants::missing::uintValue) + { + double dx = nodes[edge.endNode].x - nodes[edge.startNode].x; + double dy = nodes[edge.endNode].y - nodes[edge.startNode].y; + + // Compute angle for "pinch point" sorting + // Pinch points are when the boundary polygon intersects itself, the points will have identical values. + double angleStartToEnd = normalizeAngle(std::atan2(dy, dx)); + double angleEndToStart = normalizeAngle(std::atan2(-dy, -dx)); + + // Note: edge.leftFace is the internal mesh face valid for both directions of boundary traversal + boundaryAdjacency[edge.startNode].push_back({edge.id, edge.endNode, edge.leftFace, angleStartToEnd}); + boundaryAdjacency[edge.endNode].push_back({edge.id, edge.startNode, edge.leftFace, angleEndToStart}); + } + } + + // Sort outgoing edges anti-clockwise + for (auto& [node_id, connected_edges] : boundaryAdjacency) + { + std::sort(connected_edges.begin(), connected_edges.end(), [](const BoundaryEdge& a, const BoundaryEdge& b) + { return a.angle < b.angle; }); + } + + // Trace boundary polygons + for (const auto& edge : edges) + { + if (edge.rightFace != constants::missing::uintValue || edgeVisited[edge.id]) + { + continue; + } + + std::vector currentNodes; + std::vector currentFaces; + + meshkernel::UInt prevNode = edge.startNode; + meshkernel::UInt currentNode = edge.endNode; + + currentNodes.push_back(prevNode); + currentFaces.push_back(edge.leftFace); // First face touched by the initial edge + edgeVisited[edge.id] = true; // Mark edge as having been visited + + double dx = nodes[currentNode].x - nodes[prevNode].x; + double dy = nodes[currentNode].y - nodes[prevNode].y; + double incomingAngle = normalizeAngle(std::atan2(dy, dx)); + + // Loop until we find the start node, making a closed boundary polygon + while (currentNode != edge.startNode) + { + currentNodes.push_back(currentNode); + + const std::vector& boundaryEdges = boundaryAdjacency[currentNode]; + + // The smallest visited angle. + // All edge angles must be in interval [-2pi, 2pi], + // So 3pi is used to indicate a suitable edge angle has been found + double bestDelta = 3.0 * std::numbers::pi; + meshkernel::UInt chosenEdgeIndex = constants::missing::uintValue; + + // Find the boundary edge that tracks closest clockwise to keep empty space on the right + for (size_t i = 0; i < boundaryEdges.size(); ++i) + { + if (edgeVisited[boundaryEdges[i].edgeId]) + { + continue; + } + + double delta = normalizeAngle(incomingAngle - boundaryEdges[i].angle); + + if (delta < bestDelta) + { + bestDelta = delta; + chosenEdgeIndex = i; + } + } + + // No unvisited edges were found + // So eigher the boundary loop was fully completed or a dead-end reached. + if (chosenEdgeIndex == constants::missing::uintValue) + { + break; + } + + const BoundaryEdge& chosenEdge = boundaryEdges[chosenEdgeIndex]; + edgeVisited[chosenEdge.edgeId] = true; + + // Record the face touched by this next step in the loop + currentFaces.push_back(chosenEdge.leftFace); + + prevNode = currentNode; + currentNode = chosenEdge.neighbourNode; + incomingAngle = chosenEdge.angle; + } + + if (currentNodes.size() >= 3) + { + allLoops.push_back(std::move(currentNodes)); + allTouchedFaces.push_back(std::move(currentFaces)); + } + } +} + +// in: nodes; +// in edges +// in edgesFaces +// Out: The boundary node sequences +// Out: All faces touched per loop in step order +void findBoundaryLoops(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::vector>& allLoops, + std::vector>& allTouchedFaces) +{ + allLoops.clear(); + allTouchedFaces.clear(); + + // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle + std::unordered_map> boundaryAdjacency; + + // May be better to use meshkernel::Boolean + std::vector edgeVisited(edges.size(), false); + + // Collect all boundary edges and compute the angle + for (meshkernel::UInt count = 0; count < edges.size(); ++count) + // for (meshkernel::UInt count = 0; const auto& edge : edges) + { + // const meshkernel::Edge& edge = edges[count]; + + if (edgesFaces[count][1] == constants::missing::uintValue) + { + double dx = nodes[edges[count].second].x - nodes[edges[count].first].x; + double dy = nodes[edges[count].second].y - nodes[edges[count].first].y; + + // Compute angle for "pinch point" sorting + // Pinch points are when the boundary polygon intersects itself, the points will have identical values. + double angleStartToEnd = normalizeAngle(std::atan2(dy, dx)); + double angleEndToStart = normalizeAngle(std::atan2(-dy, -dx)); + + // Note: edgesFaces[count][0] is the internal mesh face valid for both directions of boundary traversal + boundaryAdjacency[edges[count].first].push_back({count, edges[count].second, edgesFaces[count][0], angleStartToEnd}); + boundaryAdjacency[edges[count].second].push_back({count, edges[count].first, edgesFaces[count][0], angleEndToStart}); + } + } + + auto boundaryAngleLessThan = [](const BoundaryEdge& a, const BoundaryEdge& b) + { return a.angle < b.angle; }; + + // Sort outgoing edges anti-clockwise + for (auto& [node_id, connected_edges] : boundaryAdjacency) + { + std::sort(connected_edges.begin(), connected_edges.end(), boundaryAngleLessThan); + } + + // Trace boundary polygons + for (meshkernel::UInt count = 0; count < edges.size(); ++count) + { + if (edgesFaces[count][1] != constants::missing::uintValue || edgeVisited[count]) + { + continue; + } + + std::vector currentNodes; + std::vector currentFaces; + + meshkernel::UInt prevNode = edges[count].first; + meshkernel::UInt currentNode = edges[count].second; + + currentNodes.push_back(prevNode); + currentFaces.push_back(edgesFaces[count][0]); // First face touched by the initial edge + edgeVisited[count] = true; // Mark edge as having been visited + + double dx = nodes[currentNode].x - nodes[prevNode].x; + double dy = nodes[currentNode].y - nodes[prevNode].y; + double incomingAngle = normalizeAngle(std::atan2(dy, dx)); + + // Loop until we find the start node, making a closed boundary polygon + while (currentNode != edges[count].first) + { + currentNodes.push_back(currentNode); + + const std::vector& boundaryEdges = boundaryAdjacency[currentNode]; + + // The smallest visited angle. + // All edge angles must be in interval [-2pi, 2pi], + double deltaAngle = 3.0 * std::numbers::pi; + meshkernel::UInt chosenEdgeIndex = constants::missing::uintValue; + + // Find the boundary edge that tracks closest clockwise to keep empty space on the right + for (size_t e = 0; e < boundaryEdges.size(); ++e) + { + if (edgeVisited[boundaryEdges[e].edgeId]) + { + continue; + } + + double delta = normalizeAngle(incomingAngle - boundaryEdges[e].angle); + + if (delta < deltaAngle) + { + deltaAngle = delta; + chosenEdgeIndex = e; + } + } + + // No unvisited edges were found + // So eigher the boundary loop was completed or a dead-end reached. + if (chosenEdgeIndex == constants::missing::uintValue) + { + break; + } + + const BoundaryEdge& chosenEdge = boundaryEdges[chosenEdgeIndex]; + edgeVisited[chosenEdge.edgeId] = true; + + // Record the face touched by this next step in the loop + currentFaces.push_back(chosenEdge.leftFace); + + prevNode = currentNode; + currentNode = chosenEdge.neighbourNode; + incomingAngle = chosenEdge.angle; + } + + if (currentNodes.size() >= 3) + { + currentNodes.push_back(currentNodes.front()); + allLoops.push_back(std::move(currentNodes)); + allTouchedFaces.push_back(std::move(currentFaces)); + } + } +} + +void findBoundaryLoops(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + const meshkernel::Polygons& clippingPolygon, + std::vector>& allLoops, + std::vector>& allTouchedFaces) +{ + allLoops.clear(); + allTouchedFaces.clear(); + + // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle + std::unordered_map> boundaryAdjacency; + + // May be better to use meshkernel::Boolean + std::vector edgeVisited(edges.size(), false); + + // Collect all boundary edges and compute the angle + for (meshkernel::UInt count = 0; count < edges.size(); ++count) + // for (meshkernel::UInt count = 0; const auto& edge : edges) + { + // const meshkernel::Edge& edge = edges[count]; + + if (edgesFaces[count][1] == constants::missing::uintValue) + { + double dx = nodes[edges[count].second].x - nodes[edges[count].first].x; + double dy = nodes[edges[count].second].y - nodes[edges[count].first].y; + + // Compute angle for "pinch point" sorting + // Pinch points are when the boundary polygon intersects itself, the points will have identical values. + double angleStartToEnd = normalizeAngle(std::atan2(dy, dx)); + double angleEndToStart = normalizeAngle(std::atan2(-dy, -dx)); + + // Note: edgesFaces[count][0] is the internal mesh face valid for both directions of boundary traversal + boundaryAdjacency[edges[count].first].push_back({count, edges[count].second, edgesFaces[count][0], angleStartToEnd}); + boundaryAdjacency[edges[count].second].push_back({count, edges[count].first, edgesFaces[count][0], angleEndToStart}); + } + } + + auto boundaryAngleLessThan = [](const BoundaryEdge& a, const BoundaryEdge& b) + { return a.angle < b.angle; }; + + // Sort outgoing edges anti-clockwise + for (auto& [node_id, connected_edges] : boundaryAdjacency) + { + std::sort(connected_edges.begin(), connected_edges.end(), boundaryAngleLessThan); + } + + // Trace boundary polygons + for (meshkernel::UInt count = 0; count < edges.size(); ++count) + { + if (edgesFaces[count][1] != constants::missing::uintValue || edgeVisited[count]) + { + continue; + } + + std::vector currentNodes; + std::vector currentFaces; + + meshkernel::UInt prevNode = edges[count].first; + meshkernel::UInt currentNode = edges[count].second; + + currentNodes.push_back(prevNode); + currentFaces.push_back(edgesFaces[count][0]); // First face touched by the initial edge + edgeVisited[count] = true; // Mark edge as having been visited + + double dx = nodes[currentNode].x - nodes[prevNode].x; + double dy = nodes[currentNode].y - nodes[prevNode].y; + double incomingAngle = normalizeAngle(std::atan2(dy, dx)); + + // Loop until we find the start node, making a closed boundary polygon + while (currentNode != edges[count].first) + { + currentNodes.push_back(currentNode); + + const std::vector& boundaryEdges = boundaryAdjacency[currentNode]; + + // The smallest visited angle. + // All edge angles must be in interval [-2pi, 2pi], + double deltaAngle = 3.0 * std::numbers::pi; + meshkernel::UInt chosenEdgeIndex = constants::missing::uintValue; + + // Find the boundary edge that tracks closest clockwise to keep empty space on the right + for (size_t e = 0; e < boundaryEdges.size(); ++e) + { + if (edgeVisited[boundaryEdges[e].edgeId]) + { + continue; + } + + double delta = normalizeAngle(incomingAngle - boundaryEdges[e].angle); + + if (delta < deltaAngle) + { + deltaAngle = delta; + chosenEdgeIndex = e; + } + } + + // No unvisited edges were found + // So eigher the boundary loop was completed or a dead-end reached. + if (chosenEdgeIndex == constants::missing::uintValue) + { + break; + } + + const BoundaryEdge& chosenEdge = boundaryEdges[chosenEdgeIndex]; + edgeVisited[chosenEdge.edgeId] = true; + + // Record the face touched by this next step in the loop + currentFaces.push_back(chosenEdge.leftFace); + + prevNode = currentNode; + currentNode = chosenEdge.neighbourNode; + incomingAngle = chosenEdge.angle; + } + + if (currentNodes.size() >= 3) + { + std::vector activeSegment; + std::vector activeFaces; + + for (size_t i = 0; i < currentNodes.size(); ++i) + { + + // Check if this specific node falls inside the user's bounding box/polygon + meshkernel::Point aNode = nodes[currentNodes[i]]; + auto [isIn, whichPoly] = clippingPolygon.IsPointInPolygons(aNode); + + if (isIn) + { + activeSegment.push_back(currentNodes[i]); + + if (i < currentFaces.size()) + { + activeFaces.push_back(currentFaces[i]); + } + } + else + { + // Node fell outside, push the collected segment if it holds meaningful data + if (activeSegment.size() >= 2) + { + std::cout << " activeSegment size: " << activeSegment.size() << std::endl; + activeSegment.push_back(activeSegment.front()); + allLoops.push_back(activeSegment); + allTouchedFaces.push_back(activeFaces); + } + + activeSegment.clear(); + activeFaces.clear(); + } + } + + // Push any remaining segment if the loop ended while still inside the polygon + if (activeSegment.size() >= 2) + { + std::cout << " activeSegment size: " << activeSegment.size() << std::endl; + activeSegment.push_back(activeSegment.front()); + // activeFaces.push_back(activeFaces.front()); + allLoops.push_back(std::move(activeSegment)); + allTouchedFaces.push_back(std::move(activeFaces)); + } + } + } +} + TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygons) { @@ -3078,8 +3827,169 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon // This should not fill in the holes in the mesh mesh2.Administrate(); + meshkernel::SaveVtk(mesh2.Nodes(), mesh2.m_facesNodes, "mesh10.vtu"); + + std::vector edges2(mesh2.Edges().size()); + + for (meshkernel::UInt i = 0; i < mesh2.Edges().size(); ++i) + { + edges2[i] = Edge2{i, mesh2.Edges()[i].first, mesh2.Edges()[i].second, mesh2.m_edgesFaces[i][0], mesh2.m_edgesFaces[i][1]}; + } + + auto bLoops = findBoundaryLoops(mesh2.Nodes(), edges2, mesh2.m_edgesFaces); + + [[maybe_unused]] const auto& meshNodes = mesh2.Nodes(); + + // for (size_t i = 0; i < bLoops.size(); ++i) + for (const auto& loop : bLoops) + { + std::cout << "boudnary loop: "; + + for (size_t j = 0; j < loop.size(); ++j) + { + std::cout << "{" << loop[j].x << ", " << loop[j].y << "}, "; + } + + std::cout << std::endl; + } + + std::vector> allLoops; + std::vector> allTouchedFaces; + + std::vector clippingpolygonNodes{{-5.0, -5.0}, {140.0, 0.0}, {140.0, 35.0}, {-5.0, 35.0}, {-5.0, -5.0}}; + meshkernel::Polygons clippingPolygon(clippingpolygonNodes, mesh2.m_projection); + + findBoundaryLoops(mesh2.Nodes(), mesh2.Edges(), mesh2.m_edgesFaces, allLoops, allTouchedFaces); + // findBoundaryLoops(mesh2.Nodes(), mesh2.Edges(), mesh2.m_edgesFaces, clippingPolygon, allLoops, allTouchedFaces); + + std::cout << std::endl; + std::cout << "--------------------------------" << std::endl; + + for (meshkernel::UInt count = 0; const auto& loop : bLoops) + { + std::cout << "boudnary loop: " << " " << count << " " << loop.size() << ": " << ": "; + // std::cout << "boudnary loop: " << allTouchedFaces[count].size() << " " << loop.size() << ": " << allTouchedFaces[count][0] << ": "; + + auto boundaryPolygons = splitMultiplePolygons(loop); + + std::cout << "boudnary loop: " << loop.size() << " " << boundaryPolygons.size() << " "; + + for (size_t j = 0; j < loop.size(); ++j) + { + std::cout << "{" << loop[j].x << ", " << loop[j].y << "}, "; + } + + std::cout << std::endl; + + for (size_t j = 0; j < boundaryPolygons.size(); ++j) + { + std::cout << " sub-loop " << boundaryPolygons[j].size() << " "; + + for (size_t k = 0; k < boundaryPolygons[j].size(); ++k) + { + std::cout << "{" << boundaryPolygons[j][k].x << ", " << boundaryPolygons[j][k].y << "}, "; + } + + std::cout << std::endl; + } + + std::cout << std::endl; + ++count; + } + + std::cout << std::endl; + std::cout << "----------- clipping ---------------------" << std::endl; + + for (meshkernel::UInt count = 0; const auto& loop : allLoops) + { + + std::vector loopPnts(loop.size()); + + for (size_t i = 0; i < loopPnts.size(); ++i) + { + loopPnts[i] = meshNodes[loop[i]]; + } + + // std::tuple>, std::vector> splitMultiplePolygons(std::span boundary, std::span elementIds); + auto [boundaryPolygons, firstTouchedFaces] = splitMultiplePolygons(loopPnts, allTouchedFaces[count]); + + std::cout << "boudnary loop: " << firstTouchedFaces.size() << " " << loop.size() << " " << boundaryPolygons.size() << " "; + + for (size_t j = 0; j < loop.size(); ++j) + { + std::cout << "{" << loopPnts[j].x << ", " << loopPnts[j].y << "}, "; + } + + std::cout << std::endl; + std::cout << "face ids: "; + + for (size_t j = 0; j < allTouchedFaces[count].size(); ++j) + { + std::cout << allTouchedFaces[count][j] << ", "; + } + + std::cout << std::endl; + + for (size_t j = 0; j < boundaryPolygons.size(); ++j) + { + Polygons pgs(boundaryPolygons[j], mesh2.m_projection); + std::cout << " sub-loop " << boundaryPolygons[j].size() << " " << firstTouchedFaces[j] << " " << std::boolalpha << pgs.IsPointInAnyPolygon(mesh2.m_facesMassCenters[firstTouchedFaces[j]]) << ": "; + + for (size_t k = 0; k < boundaryPolygons[j].size(); ++k) + { + std::cout << "{" << boundaryPolygons[j][k].x << ", " << boundaryPolygons[j][k].y << "}, "; + } + + std::cout << std::endl; + } + + std::cout << std::endl; + ++count; + } + // Get interior boundary polygon points std::vector innerBoundaryPoints = mesh2.GetInnerBoundaryPolygons(); + std::vector outerBoundaryPoints = mesh2.ComputeBoundaryPolygons(std::vector()); + + std::cout << std::endl; + std::cout << "--------------------------------" << std::endl; + std::cout << "boudnary loop: "; + + for (size_t j = 0; j < outerBoundaryPoints.size(); ++j) + { + if (outerBoundaryPoints[j].IsValid()) + { + std::cout << "{" << outerBoundaryPoints[j].x << ", " << outerBoundaryPoints[j].y << "}, "; + } + else + { + std::cout << std::endl; + std::cout << "boundary loop"; + } + } + + std::cout << std::endl; + + MeshBoundaryExtractor extractor; + + auto allBoundaryPoints = extractor.ExtractAll(mesh2); + + Polygons allBoundaryPolygons(allBoundaryPoints, mesh2.m_projection); + + for (UInt i = 0; i < allBoundaryPolygons.GetNumPolygons(); ++i) + { + const auto& enclosure = allBoundaryPolygons.Enclosure(i); + const auto& enclosurePoints = enclosure.Outer().Nodes(); + + std::cout << "enclosure points: "; + + for (size_t j = 0; j < enclosurePoints.size(); ++j) + { + std::cout << "{" << enclosurePoints[j].x << ", " << enclosurePoints[j].y << "}, "; + } + + std::cout << std::endl; + } // The expected number of points, should not include any land boundary points constexpr UInt expectedNumberOfNodes = 26; From 92e3d0da24c06d2c9a55e99856fb575ec2f377ca Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 3 Aug 2026 18:14:46 +0200 Subject: [PATCH 13/34] GRIDEDIT-2293 Passing to other computer. --- .../include/MeshKernel/Definitions.hpp | 8 + .../MeshKernel/MeshBoundaryExtractor.hpp | 14 +- .../include/MeshKernel/Operations.hpp | 41 +- libs/MeshKernel/src/Mesh2D.cpp | 76 +- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 120 ++- libs/MeshKernel/src/Operations.cpp | 156 ++- libs/MeshKernel/src/Polygon.cpp | 92 +- libs/MeshKernel/tests/src/Mesh2DTest.cpp | 46 +- .../tests/src/MeshRefinementTests.cpp | 942 +----------------- .../tests/src/InvalidCellsPolygonsTests.cpp | 6 +- 10 files changed, 379 insertions(+), 1122 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Definitions.hpp b/libs/MeshKernel/include/MeshKernel/Definitions.hpp index c02d6c250..8cc1ba3c3 100644 --- a/libs/MeshKernel/include/MeshKernel/Definitions.hpp +++ b/libs/MeshKernel/include/MeshKernel/Definitions.hpp @@ -174,4 +174,12 @@ namespace meshkernel AllNetlinksLoop }; + /// \brief Boundary selection indicator + enum class BoundarySelection + { + ExteriorOnly, ///< Exterior boundaries only + InteriorOnly, ///< Interior boundaries only + All ///< All boundaries, exterior and interior + }; + } // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp index f3b40ce1d..416714b06 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp @@ -41,10 +41,14 @@ namespace meshkernel class MeshBoundaryExtractor { public: - /// @brief Extract all boundaries, concatinated as a single sequence of points, separated by an invalid point - static std::vector ExtractAll(const Mesh2D& mesh); + /// @brief Extract all boundaries as a single sequence of points, separated by an invalid point + static std::vector ExtractConcatenated(const Mesh2D& mesh, BoundarySelection boundaryType = BoundarySelection::All); /// @brief Extract all boundaries keeping them separated and + /// + /// The result consists of an array of each of the boundary polygons + /// Additionally, an array indicating if the boundary polygon is a exterior boudnary or not. + /// True => is-exterior, False => otherwise static std::tuple>, std::vector> Extract(const Mesh2D& mesh); private: @@ -60,6 +64,12 @@ namespace meshkernel /// @brief Ensure the angle lies between 0 .. 2 pi. static double NormalizeAngle(double angle); + /// @brief Construct mapping from node-id to all impinging edges for boundary all edges + static void FindAllBoundarEdges(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::unordered_map>& boundaryAdjacency); + /// @brief Find boundary loops /// /// Any boundary loops found may need to be processed further as they may themselves contain sub-loops diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 77a025d87..7d6ec8d01 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -298,6 +298,12 @@ namespace meshkernel const std::span triangleNodes, const Projection& projection); + /// \brief Compute the area and the centre of mass of a polygonal area + /// + /// The area of the polygon will be positive is the boundary points traverse the boundary in a anti clockwise direction, + /// If the are is negative then the points traverse the boundary in a clockwise direction. + std::tuple ComputePolygonAreaAndCentre(const std::vector& boundaryPolygon, const Projection& projection); + /// @brief Computes three base components void ComputeThreeBaseComponents(const Point& point, std::array& exxp, std::array& eyyp, std::array& ezzp); @@ -624,16 +630,39 @@ namespace meshkernel /// @brief Determine is a polygon is comprised of multiple polygons intersecting at the control points only /// /// @note This currently works only for polygons that intersect at the control points - bool isMultiPolygon(std::span boundary); + bool IsMultiPolygon(std::span boundary); /// @brief Split a single polygon line comprised of multiple polygons into separate polygons. /// /// @note This currently works only for polygons that intersect at the control points - std::tuple>, std::vector> splitMultiplePolygons(std::span boundary, std::span elementIds); + std::tuple>, std::vector> SplitMultiplePolygons(std::span boundary, std::span elementIds); - /// @brief Split a single polygon line comprised of multiple polygons into separate polygons. - /// - /// @note This currently works only for polygons that intersect at the control points - std::vector> splitMultiplePolygons(std::span boundary); + /// @brief Concaenate vectors of vectors of points to a vector of points separated by the invalid point + template + std::vector ConcatenatePointVectors(const std::vector>& pointVectors, Predicate predicate) + { + + std::vector combinedPoints; + bool isFirst = true; + + for (size_t count = 0; const auto& points : pointVectors) + { + + if (predicate(count)) + { + if (!isFirst) + { + combinedPoints.push_back({constants::missing::doubleValue, constants::missing::doubleValue}); + } + + combinedPoints.insert(combinedPoints.end(), points.begin(), points.end()); + isFirst = false; + } + + ++count; + } + + return combinedPoints; + } } // namespace meshkernel diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 966eee6a6..077a224f8 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -33,6 +33,7 @@ #include "MeshKernel/Exceptions.hpp" #include "MeshKernel/Mesh2D.hpp" #include "MeshKernel/Mesh2DIntersections.hpp" +#include "MeshKernel/MeshBoundaryExtractor.hpp" #include "MeshKernel/MeshFaceCenters.hpp" #include "MeshKernel/MeshOrthogonality.hpp" #include "MeshKernel/Operations.hpp" @@ -1521,6 +1522,64 @@ std::tuple, std::vector> Mesh2D::GetAllBoun { const Polygon polygon(polygonNodes, m_projection); +#if 1 + MeshBoundaryExtractor meshBoundaryExtractor; + + auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this); + + if (polygonNodes.size() == 0) + { + + auto isTrue = [](size_t idx [[maybe_unused]]) + { + return true; + }; + + return {ConcatenatePointVectors(boundaryPoints, isTrue), isEnclosingBoundary}; + } + + std::vector containedIsEnclosingBoundary; + + auto removeExteriorPoints = [&polygon](const Point& p) + { + return !polygon.Contains(p); + }; + + for (size_t i = 0; i < boundaryPoints.size(); ++i) + { + std::erase_if(boundaryPoints[i], removeExteriorPoints); + } + + auto addNonEmpty = [&boundaryPoints, &isEnclosingBoundary, &containedIsEnclosingBoundary](size_t idx) mutable + { + if (boundaryPoints[idx].size() > 0) + { + containedIsEnclosingBoundary.push_back(isEnclosingBoundary[idx]); + return true; + } + + return false; + }; + + std::vector containedBoundaryPoints = ConcatenatePointVectors(boundaryPoints, addNonEmpty); + + return {containedBoundaryPoints, containedIsEnclosingBoundary}; + +#if 0 + for (size_t i = 0; i < boundaryPoints.size(); ++i) + { + + if (boundaryPoints[i].size() > 0) + { + containedBoundaryPoints.insert(containedBoundaryPoints.end(), boundaryPoints[i].begin(), boundaryPoints[i].end()); + containedIsEnclosingBoundary.push_back(isEnclosingBoundary[i]); + } + } + + return {containedBoundaryPoints, containedIsEnclosingBoundary}; +#endif + +#else std::vector isVisited(GetNumEdges(), false); std::vector meshBoundaryPolygon; std::vector boundaryPolygon; @@ -1591,10 +1650,10 @@ std::tuple, std::vector> Mesh2D::GetAllBoun meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); } - if (isMultiPolygon(currentPolygonSpan)) + if (IsMultiPolygon(currentPolygonSpan)) { - auto [multiBoundaryPolygonNodes, multiBoundaryElementIds] = (splitMultiplePolygons(currentPolygonSpan, boundaryPolygonFaceId)); + auto [multiBoundaryPolygonNodes, multiBoundaryElementIds] = (SplitMultiplePolygons(currentPolygonSpan, boundaryPolygonFaceId)); for (size_t i = 0; i < multiBoundaryPolygonNodes.size(); ++i) { @@ -1623,6 +1682,7 @@ std::tuple, std::vector> Mesh2D::GetAllBoun } return {meshBoundaryPolygon, isEnclosingBoundary}; +#endif } std::vector Mesh2D::ComputeInnerBoundaryPolygons() const @@ -2237,11 +2297,11 @@ std::unique_ptr Mesh2D::DeleteMeshFacesInPolygon(const P void Mesh2D::ReconstructInvalidCellsPolygon() { - auto [boundaryPoints, isEnclosingBoundary] = GetAllBoundaryPolygons(std::vector()); + MeshBoundaryExtractor meshBoundaryExtractor; - Polygons polygons(boundaryPoints, m_projection); + auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this); - if (polygons.GetNumPolygons() <= 1) + if (boundaryPoints.size() <= 1) { // There are no interior boundary polygons return; @@ -2251,8 +2311,9 @@ void Mesh2D::ReconstructInvalidCellsPolygon() innerBoundaryPolygons.reserve(m_invalidCellPolygons.size()); bool firstElement = true; - for (UInt p = 0; p < polygons.GetNumPolygons(); ++p) + for (UInt p = 0; p < boundaryPoints.size(); ++p) { + if (isEnclosingBoundary[p]) { continue; @@ -2267,8 +2328,7 @@ void Mesh2D::ReconstructInvalidCellsPolygon() firstElement = false; } - const std::vector& polygonPoints(polygons.Enclosure(p).Outer().Nodes()); - innerBoundaryPolygons.insert(innerBoundaryPolygons.end(), polygonPoints.begin(), polygonPoints.end()); + innerBoundaryPolygons.insert(innerBoundaryPolygons.end(), boundaryPoints[p].begin(), boundaryPoints[p].end()); } m_invalidCellPolygons = std::move(innerBoundaryPolygons); diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index 69c15cfd0..6af50b626 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -1,19 +1,48 @@ #include "MeshKernel/MeshBoundaryExtractor.hpp" +#include #include #include "MeshKernel/Operations.hpp" -std::vector meshkernel::MeshBoundaryExtractor::ExtractAll(const Mesh2D& mesh) +std::vector meshkernel::MeshBoundaryExtractor::ExtractConcatenated(const Mesh2D& mesh, BoundarySelection boundaryType) { auto [boundarySequences, isExterior] = Extract(mesh); std::vector allPoints; + auto isBoundarySelection = [boundaryType, &isExterior](size_t idx) + { + using enum BoundarySelection; + + return (boundaryType == All) || + (boundaryType == ExteriorOnly && isExterior[idx]) || + (boundaryType == InteriorOnly && !isExterior[idx]); + }; + + return ConcatenatePointVectors(boundarySequences, isBoundarySelection); + +#if 0 bool isFirst = true; - for (const std::vector& loop : boundarySequences) + auto selectedBoundary = [boundaryType](bool isExterior) + { + using enum BoundarySelection; + + return (boundaryType == All) || + (boundaryType == ExteriorOnly && isExterior) || + (boundaryType == InteriorOnly && !isExterior); + }; + + for (size_t i = 0; i < boundarySequences.size(); ++i) { + const std::vector& loop = boundarySequences[i]; + + if (!selectedBoundary(isExterior[i])) + { + continue; + } + if (!isFirst) { allPoints.push_back({constants::missing::doubleValue, constants::missing::doubleValue}); @@ -24,6 +53,7 @@ std::vector meshkernel::MeshBoundaryExtractor::ExtractAll(con } return allPoints; +#endif } std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh) @@ -41,13 +71,24 @@ std::tuple>, std::vector> meshk for (size_t i = 0; i < allBoundaryLoops.size(); ++i) { - if (isMultiPolygon(allBoundaryLoops[i])) + + if (IsMultiPolygon(allBoundaryLoops[i])) { - auto [individualBoundaryPolygons, firstElement] = splitMultiplePolygons(allBoundaryLoops[i], allTouchedFaces[i]); + // It can be that some of the boundaries that are found are composed of multiple sub-boundaries. + // Some of the sub-boundaries may be combined in a non conforming way. + // In either case, the boundaries are separated into distinct boundary polygons. + auto [individualBoundaryPolygons, firstElement] = SplitMultiplePolygons(allBoundaryLoops[i], allTouchedFaces[i]); for (size_t i = 0; i < individualBoundaryPolygons.size(); ++i) { Point centre = mesh.m_facesMassCenters[firstElement[i]]; + + // If the area is calculated ot be less than zero, i.e. the boundary is traversed in clockwise direction and needs to be reversed + if (auto [area, centreOfMass] = ComputePolygonAreaAndCentre(individualBoundaryPolygons[i], mesh.m_projection); area < 0.0) + { + std::ranges::reverse(individualBoundaryPolygons[i]); + } + isExterior.push_back(IsPointInPolygonNodes(centre, individualBoundaryPolygons[i], mesh.m_projection)); separatedBoundaryLoops.push_back(std::move(individualBoundaryPolygons[i])); } @@ -76,39 +117,33 @@ double meshkernel::MeshBoundaryExtractor::NormalizeAngle(double angle) return angle; } -void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector& nodes, - const std::vector& edges, - const std::vector>& edgesFaces, - std::vector>& allLoops, - std::vector>& allTouchedFaces) +void meshkernel::MeshBoundaryExtractor::FindAllBoundarEdges(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::unordered_map>& boundaryAdjacency) { - allLoops.clear(); - allTouchedFaces.clear(); - - // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle - std::unordered_map> boundaryAdjacency; - - // May be better to use meshkernel::Boolean - std::vector edgeVisited(edges.size(), false); // Collect all boundary edges and compute the angle for (UInt count = 0; count < edges.size(); ++count) { - if (edgesFaces[count][1] == constants::missing::uintValue) + if (!IsValidEdge(edges[count]) || edgesFaces[count][1] != constants::missing::uintValue) { - double dx = nodes[edges[count].second].x - nodes[edges[count].first].x; - double dy = nodes[edges[count].second].y - nodes[edges[count].first].y; + // Edge is either invalid or not on boundaryy + continue; + } - // Compute angle for "pinch point" sorting - // Pinch points are when the boundary polygon intersects itself, the points will have identical values. - double angleStartToEnd = NormalizeAngle(std::atan2(dy, dx)); - double angleEndToStart = NormalizeAngle(std::atan2(-dy, -dx)); + double dx = nodes[edges[count].second].x - nodes[edges[count].first].x; + double dy = nodes[edges[count].second].y - nodes[edges[count].first].y; - // Note: edgesFaces[count][0] is the internal mesh face valid for both directions of boundary traversal - boundaryAdjacency[edges[count].first].push_back({count, edges[count].second, edgesFaces[count][0], angleStartToEnd}); - boundaryAdjacency[edges[count].second].push_back({count, edges[count].first, edgesFaces[count][0], angleEndToStart}); - } + // Compute angle for "crossong/pinch point" sorting + // Pinch points are when the boundary polygon intersects itself, the points will have identical values. + double angleStartToEnd = NormalizeAngle(std::atan2(dy, dx)); + double angleEndToStart = NormalizeAngle(std::atan2(-dy, -dx)); + + // Note: edgesFaces[count][0] is the internal mesh face valid for both directions of boundary traversal + boundaryAdjacency[edges[count].first].push_back({count, edges[count].second, edgesFaces[count][0], angleStartToEnd}); + boundaryAdjacency[edges[count].second].push_back({count, edges[count].first, edgesFaces[count][0], angleEndToStart}); } auto boundaryAngleLessThan = [](const BoundaryEdge& a, const BoundaryEdge& b) @@ -119,10 +154,33 @@ void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::vector>& allLoops, + std::vector>& allTouchedFaces) +{ + allLoops.clear(); + allTouchedFaces.clear(); + + // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle + std::unordered_map> boundaryAdjacency; + + FindAllBoundarEdges(nodes, edges, edgesFaces, boundaryAdjacency); + + std::vector edgeVisited(edges.size(), false); // Trace boundary polygons for (UInt count = 0; count < edges.size(); ++count) { + + if (!IsValidEdge(edges[count])) + { + continue; + } + if (edgesFaces[count][1] != constants::missing::uintValue || edgeVisited[count]) { continue; @@ -149,8 +207,8 @@ void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector& boundaryEdges = boundaryAdjacency[currentNodeIndex]; - // The smallest visited angle. - // All edge angles must be in interval [-2pi, 2pi], + // Find the smallest visited angle. + // All edge angles must be in interval [0, 2pi], initialise with number greater than 2pi double deltaAngle = 3.0 * std::numbers::pi; UInt chosenEdgeIndex = constants::missing::uintValue; @@ -181,7 +239,7 @@ void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector ComputePolygonAreaAndCentre(const std::vector& polygon, const Projection& projection) + { + + if (polygon.size() < constants::geometric::numNodesInTriangle) + { + throw std::invalid_argument("FaceAreaAndCenterOfMass: The polygon has less than 3 unique nodes."); + } + + Point centreOfMass(0.0, 0.0); + double area = 0.0; + + const double minArea = 1e-8; + const Point reference = ReferencePoint(polygon, projection); + const auto numberOfPointsOpenedPolygon = static_cast(polygon.size()) - 1; + + if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInTriangle) + { + Vector delta1 = GetDelta(reference, polygon[0], projection); + Vector delta2 = GetDelta(reference, polygon[1], projection); + Vector delta3 = GetDelta(reference, polygon[2], projection); + + Vector middle1 = 0.5 * (delta1 + delta2); + Vector middle2 = 0.5 * (delta2 + delta3); + Vector middle3 = 0.5 * (delta3 + delta1); + + delta1 = GetDelta(polygon[0], polygon[1], projection); + delta2 = GetDelta(polygon[1], polygon[2], projection); + delta3 = GetDelta(polygon[2], polygon[0], projection); + + double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); + double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); + double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); + + area = 0.5 * (xds1 + xds2 + xds3); + + centreOfMass += xds1 * middle1; + centreOfMass += xds2 * middle2; + centreOfMass += xds3 * middle3; + } + else if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInQuadrilateral) + { + Vector delta1 = GetDelta(reference, polygon[0], projection); + Vector delta2 = GetDelta(reference, polygon[1], projection); + Vector delta3 = GetDelta(reference, polygon[2], projection); + Vector delta4 = GetDelta(reference, polygon[3], projection); + + Vector middle1 = 0.5 * (delta1 + delta2); + Vector middle2 = 0.5 * (delta2 + delta3); + Vector middle3 = 0.5 * (delta3 + delta4); + Vector middle4 = 0.5 * (delta4 + delta1); + + delta1 = GetDelta(polygon[0], polygon[1], projection); + delta2 = GetDelta(polygon[1], polygon[2], projection); + delta3 = GetDelta(polygon[2], polygon[3], projection); + delta4 = GetDelta(polygon[3], polygon[0], projection); + + double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); + double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); + double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); + double xds4 = delta4.y() * middle4.x() - delta4.x() * middle4.y(); + + area = 0.5 * (xds1 + xds2 + xds3 + xds4); + + centreOfMass += xds1 * middle1; + centreOfMass += xds2 * middle2; + centreOfMass += xds3 * middle3; + centreOfMass += xds4 * middle4; + } + else + { + + for (UInt n = 0; n < numberOfPointsOpenedPolygon; ++n) + { + const auto nextNode = NextCircularForwardIndex(n, numberOfPointsOpenedPolygon); + + Vector delta = GetDelta(reference, polygon[n], projection); + Vector deltaNext = GetDelta(reference, polygon[nextNode], projection); + Vector middle = 0.5 * (delta + deltaNext); + delta = GetDelta(polygon[n], polygon[nextNode], projection); + + // Rotate by 3pi/2 + Vector normal(delta.y(), -delta.x()); + double xds = dot(normal, middle); + area += 0.5 * xds; + + centreOfMass += xds * middle; + } + } + + area = std::abs(area) < minArea ? minArea : area; + centreOfMass *= 1.0 / (3.0 * area); + + // TODO SHould this also apply to spheciral accurate? + if (projection == Projection::spherical) + { + centreOfMass.y /= (constants::geometric::earth_radius * constants::conversion::degToRad); + centreOfMass.x /= (constants::geometric::earth_radius * constants::conversion::degToRad * std::cos((centreOfMass.y + reference.y) * constants::conversion::degToRad)); + } + + centreOfMass += reference; + + return {area, centreOfMass}; + } + void ComputeThreeBaseComponents(const Point& point, std::array& exxp, std::array& eyyp, std::array& ezzp) { const double phi0 = point.y * constants::conversion::degToRad; @@ -1724,7 +1828,7 @@ namespace meshkernel return (matCoefficients[0] * x[0] + matCoefficients[1] * x[1]) * y[0] + (matCoefficients[2] * x[0] + matCoefficients[3] * x[1]) * y[1]; } - bool isMultiPolygon(std::span boundary) + bool IsMultiPolygon(std::span boundary) { if (boundary.size() < 4) @@ -1745,7 +1849,7 @@ namespace meshkernel return false; } - std::tuple>, std::vector> splitMultiplePolygons(std::span boundaryPoints, std::span elementIds) + std::tuple>, std::vector> SplitMultiplePolygons(std::span boundaryPoints, std::span elementIds) { std::vector> completedPolygons; std::vector firstElementIds; @@ -1795,52 +1899,4 @@ namespace meshkernel return {completedPolygons, firstElementIds}; } - std::vector> splitMultiplePolygons(std::span boundaryPoints) - { - // std::vector dummyConnectedFaces (boundaryPoints.size () - 1); - // std::ranges::iota (dummyConnectedFaces, 0); - // auto [completedPolygons, dummyFirstTouchedFaces] = splitMultiplePolygons (boundaryPoints, dummyConnectedFaces); - - std::vector> completedPolygons; - - std::vector pointFaceStack; - // Maps a point to its current index in the pointFaceStack - std::map activePoints; - - for (size_t i = 0; i < boundaryPoints.size(); ++i) - { - const Point& current_point = boundaryPoints[i]; - - // Check if this node closes a loop with a previously visited point - if (auto it = activePoints.find(current_point); it != activePoints.end()) - { - size_t loopStartIndex = it->second; - - std::vector subPolygon; - subPolygon.reserve((pointFaceStack.size() - loopStartIndex) + 1); - - for (size_t j = loopStartIndex; j < pointFaceStack.size(); ++j) - { - subPolygon.push_back(pointFaceStack[j]); - activePoints.erase(pointFaceStack[j]); - } - - // close the polygon - if (!subPolygon.empty()) - { - subPolygon.push_back(subPolygon.front()); - } - - completedPolygons.push_back(subPolygon); - - pointFaceStack.resize(loopStartIndex); - } - - activePoints[current_point] = pointFaceStack.size(); - pointFaceStack.emplace_back(current_point); - } - - return completedPolygons; - } - } // namespace meshkernel diff --git a/libs/MeshKernel/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp index 119a9186d..c9f045dbd 100644 --- a/libs/MeshKernel/src/Polygon.cpp +++ b/libs/MeshKernel/src/Polygon.cpp @@ -754,100 +754,10 @@ std::tuple meshkernel throw std::invalid_argument("FaceAreaAndCenterOfMass: The polygon has less than 3 unique nodes."); } - Point centreOfMass(0.0, 0.0); - double area = 0.0; - - const double minArea = 1e-8; - const Point reference = ReferencePoint(polygon, projection); - const auto numberOfPointsOpenedPolygon = static_cast(polygon.size()) - 1; - - if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInTriangle) - { - Vector delta1 = GetDelta(reference, polygon[0], projection); - Vector delta2 = GetDelta(reference, polygon[1], projection); - Vector delta3 = GetDelta(reference, polygon[2], projection); - - Vector middle1 = 0.5 * (delta1 + delta2); - Vector middle2 = 0.5 * (delta2 + delta3); - Vector middle3 = 0.5 * (delta3 + delta1); - - delta1 = GetDelta(polygon[0], polygon[1], projection); - delta2 = GetDelta(polygon[1], polygon[2], projection); - delta3 = GetDelta(polygon[2], polygon[0], projection); - - double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); - double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); - double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); - - area = 0.5 * (xds1 + xds2 + xds3); - - centreOfMass += xds1 * middle1; - centreOfMass += xds2 * middle2; - centreOfMass += xds3 * middle3; - } - else if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInQuadrilateral) - { - Vector delta1 = GetDelta(reference, polygon[0], projection); - Vector delta2 = GetDelta(reference, polygon[1], projection); - Vector delta3 = GetDelta(reference, polygon[2], projection); - Vector delta4 = GetDelta(reference, polygon[3], projection); - - Vector middle1 = 0.5 * (delta1 + delta2); - Vector middle2 = 0.5 * (delta2 + delta3); - Vector middle3 = 0.5 * (delta3 + delta4); - Vector middle4 = 0.5 * (delta4 + delta1); - - delta1 = GetDelta(polygon[0], polygon[1], projection); - delta2 = GetDelta(polygon[1], polygon[2], projection); - delta3 = GetDelta(polygon[2], polygon[3], projection); - delta4 = GetDelta(polygon[3], polygon[0], projection); - - double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); - double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); - double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); - double xds4 = delta4.y() * middle4.x() - delta4.x() * middle4.y(); - - area = 0.5 * (xds1 + xds2 + xds3 + xds4); - - centreOfMass += xds1 * middle1; - centreOfMass += xds2 * middle2; - centreOfMass += xds3 * middle3; - centreOfMass += xds4 * middle4; - } - else - { - - for (UInt n = 0; n < numberOfPointsOpenedPolygon; ++n) - { - const auto nextNode = NextCircularForwardIndex(n, numberOfPointsOpenedPolygon); - - Vector delta = GetDelta(reference, polygon[n], projection); - Vector deltaNext = GetDelta(reference, polygon[nextNode], projection); - Vector middle = 0.5 * (delta + deltaNext); - delta = GetDelta(polygon[n], polygon[nextNode], projection); - - // Rotate by 3pi/2 - Vector normal(delta.y(), -delta.x()); - double xds = dot(normal, middle); - area += 0.5 * xds; - - centreOfMass += xds * middle; - } - } + auto [area, centreOfMass] = ComputePolygonAreaAndCentre(polygon, projection); TraversalDirection direction = area > 0.0 ? TraversalDirection::AntiClockwise : TraversalDirection::Clockwise; - area = std::abs(area) < minArea ? minArea : area; - centreOfMass *= 1.0 / (3.0 * area); - - // TODO SHould this also apply to spheciral accurate? - if (projection == Projection::spherical) - { - centreOfMass.y /= (constants::geometric::earth_radius * constants::conversion::degToRad); - centreOfMass.x /= (constants::geometric::earth_radius * constants::conversion::degToRad * std::cos((centreOfMass.y + reference.y) * constants::conversion::degToRad)); - } - - centreOfMass += reference; return {std::abs(area), centreOfMass, direction}; } diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index f4c8e5a8d..6c5bf145e 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -242,9 +242,9 @@ TEST(Mesh2D, MeshBoundaryToPolygon) ASSERT_NEAR(0.0, meshBoundaryPolygon[4].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].y, tolerance); - ASSERT_NEAR(5.0, meshBoundaryPolygon[1].y, tolerance); + ASSERT_NEAR(-5.0, meshBoundaryPolygon[1].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[2].y, tolerance); - ASSERT_NEAR(-5.0, meshBoundaryPolygon[3].y, tolerance); + ASSERT_NEAR(5.0, meshBoundaryPolygon[3].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[4].y, tolerance); } @@ -307,25 +307,51 @@ TEST(Mesh2D, MeshBoundaryToPolygonWithSelection) // 2 Execution const auto meshBoundaryPolygon = mesh.ComputeBoundaryPolygons(polygonNodes); + meshkernel::Print(mesh.Nodes(), mesh.Edges()); + + std::cout << std::endl; + + for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) + { + std::cout << meshBoundaryPolygon[i].x << ", "; + } + + std::cout << std::endl; + + for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) + { + std::cout << meshBoundaryPolygon[i].y << ", "; + } + + std::cout << std::endl; + // 3 Validation const double tolerance = 1e-5; ASSERT_EQ(9, meshBoundaryPolygon.size()); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].y, tolerance); - ASSERT_NEAR(0.0, meshBoundaryPolygon[1].x, tolerance); - ASSERT_NEAR(10.0, meshBoundaryPolygon[1].y, tolerance); - ASSERT_NEAR(0.0, meshBoundaryPolygon[2].x, tolerance); - ASSERT_NEAR(20.0, meshBoundaryPolygon[2].y, tolerance); - ASSERT_NEAR(0.0, meshBoundaryPolygon[3].x, tolerance); - ASSERT_NEAR(30.0, meshBoundaryPolygon[3].y, tolerance); + + ASSERT_NEAR(10.0, meshBoundaryPolygon[7].x, tolerance); + ASSERT_NEAR(00.0, meshBoundaryPolygon[7].y, tolerance); + ASSERT_NEAR(10.0, meshBoundaryPolygon[4].x, tolerance); ASSERT_NEAR(30.0, meshBoundaryPolygon[4].y, tolerance); + + ASSERT_NEAR(0.0, meshBoundaryPolygon[3].x, tolerance); + ASSERT_NEAR(30.0, meshBoundaryPolygon[3].y, tolerance); + + ASSERT_NEAR(0.0, meshBoundaryPolygon[2].x, tolerance); + ASSERT_NEAR(20.0, meshBoundaryPolygon[2].y, tolerance); + + ASSERT_NEAR(0.0, meshBoundaryPolygon[1].x, tolerance); + ASSERT_NEAR(10.0, meshBoundaryPolygon[1].y, tolerance); + ASSERT_NEAR(20.0, meshBoundaryPolygon[5].x, tolerance); ASSERT_NEAR(30.0, meshBoundaryPolygon[5].y, tolerance); + ASSERT_NEAR(20.0, meshBoundaryPolygon[6].x, tolerance); ASSERT_NEAR(00.0, meshBoundaryPolygon[6].y, tolerance); - ASSERT_NEAR(10.0, meshBoundaryPolygon[7].x, tolerance); - ASSERT_NEAR(00.0, meshBoundaryPolygon[7].y, tolerance); + ASSERT_NEAR(0.0, meshBoundaryPolygon[8].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[8].y, tolerance); } diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index df618741d..20fbad288 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -2963,754 +2963,6 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa } } -#include -#include -#include -#include - -struct Edge2 -{ - meshkernel::UInt id; - meshkernel::UInt startNode; - meshkernel::UInt endNode; - meshkernel::UInt leftFace; - meshkernel::UInt rightFace; -}; - -struct BoundaryEdge -{ - meshkernel::UInt edgeId; - meshkernel::UInt neighbourNode; - meshkernel::UInt leftFace; // Store face mapping on the edge for easy retrieval during loop trace - double angle; // Angle of the edge pointing away from the pivot node -}; - -double normalizeAngle(double angle) -{ - while (angle < 0) - { - angle += 2.0 * std::numbers::pi; - } - - while (angle >= 2.0 * std::numbers::pi) - { - angle -= 2.0 * std::numbers::pi; - } - - return angle; -} - -namespace meshkernel -{ - static const UInt nullValue = constants::missing::uintValue; - -} - -// Main function to extract loops -std::vector> findBoundaryLoops(const std::vector& nodes, - const std::vector& edges, - const std::vector>& edgesFaces [[maybe_unused]]) -{ - - // // Step 1 & 2: Build a radially sorted boundary adjacency list - // // Node ID -> list of connected boundary edges - // std::unordered_map> boundaryAdjacency; - // std::unordered_map edgeVisited; - - // for (meshkernel::UInt edgeId = 0; const auto& edge : edges) - // { - // if (edgesFaces[edgeId][1] == nullValue) - // { - // edgeVisited[edgeId] = false; - - // // Calculate geometric angles for pinch-point sorting - // double dx = nodes[edge.second].x - nodes[edge.first].x; - // double dy = nodes[edge.second].y - nodes[edge.first].y; - - // double angle_start_to_end = std::atan2(dy, dx); - // double angle_end_to_start = std::atan2(-dy, -dx); - - // boundaryAdjacency[edge.first].push_back({edgeId, edge.second, angle_start_to_end}); - // boundaryAdjacency[edge.second].push_back({edgeId, edge.first, angle_end_to_start}); - // } - - // ++edgeId; - // } - - // // Sort outgoing boundary edges at each node counter-clockwise - // for (auto& [node_id, connected_edges] : boundaryAdjacency) - // { - // std::sort(connected_edges.begin(), connected_edges.end(), [](const BoundaryEdge& a, const BoundaryEdge& b) - // { return a.angle < b.angle; }); - // } - - // std::vector> allLoops; - - // // Step 3: Trace loops - // for (meshkernel::UInt edgeId = 0; const auto& edge : edges) - // { - // if (edgesFaces[edgeId][1] != nullValue || edgeVisited[edgeId]) - // continue; - - // std::vector currentLoop; - // meshkernel::UInt currentNode = edge.first; - // meshkernel::UInt startNode = currentNode; - - // while (true) - // { - // currentLoop.push_back(currentNode); - - // // Find the next unvisited boundary edge radiating from currentNode - // auto& options = boundaryAdjacency[currentNode]; - // meshkernel::UInt nextNode = nullValue; - // meshkernel::UInt chosenEdgeId = nullValue; - - // for (const auto& option : options) - // { - // if (!edgeVisited[option.edgeId]) - // { - // chosenEdgeId = option.edgeId; - // nextNode = option.neighbourNode; - // break; - // } - // } - - // if (chosenEdgeId == nullValue) - // break; // Loop fully completed or dead end reached - - // edgeVisited[chosenEdgeId] = true; - // currentNode = nextNode; - - // if (currentNode == startNode) - // { - // break; // Loop closed successfully - // } - // } - - // if (!currentLoop.empty()) - // { - // currentLoop.push_back(currentLoop.front()); - - // allLoops.push_back(currentLoop); - // } - - // ++edgeId; - // } - - //-------------------------------- - - // Node ID -> list of connected boundary edges - std::unordered_map> boundaryAdjacency; - std::unordered_map edgeVisited; - - // 1. Collect all boundary edges (where rightFace is nullValue) - for (const auto& edge : edges) - { - if (edge.rightFace == constants::missing::uintValue) - { - edgeVisited[edge.id] = false; - - double dx = nodes[edge.endNode].x - nodes[edge.startNode].x; - double dy = nodes[edge.endNode].y - nodes[edge.startNode].y; - - // Angular direction of the edge - double angle_start_to_end = normalizeAngle(std::atan2(dy, dx)); - double angle_end_to_start = normalizeAngle(std::atan2(-dy, -dx)); - - boundaryAdjacency[edge.startNode].push_back({edge.id, edge.endNode, 0, angle_start_to_end}); - boundaryAdjacency[edge.endNode].push_back({edge.id, edge.startNode, 0, angle_end_to_start}); - } - } - - // 2. Sort outgoing edges counter-clockwise around every node - for (auto& [node_id, connected_edges] : boundaryAdjacency) - { - std::sort(connected_edges.begin(), connected_edges.end(), [](const BoundaryEdge& a, const BoundaryEdge& b) - { return a.angle < b.angle; }); - } - - std::vector> allLoops; - // 3. Trace loops tracking the incoming direction - for (const auto& edge : edges) - { - if (edge.rightFace != constants::missing::uintValue || edgeVisited[edge.id]) - continue; - - std::vector currentLoop; - - // Start tracing this edge from startNode to endNode - meshkernel::UInt prevNode = edge.startNode; - meshkernel::UInt currentNode = edge.endNode; - - currentLoop.push_back(nodes[prevNode]); - edgeVisited[edge.id] = true; - - // Calculate initial incoming angle into the second node - double dx = nodes[currentNode].x - nodes[prevNode].x; - double dy = nodes[currentNode].y - nodes[prevNode].y; - double incomingAngle = normalizeAngle(std::atan2(dy, dx)); - - while (currentNode != edge.startNode) - { - currentLoop.push_back(nodes[currentNode]); - - const auto& options = boundaryAdjacency[currentNode]; - double bestDelta = 3.0 * std::numbers::pi; // Higher than max possible delta (2*PI) - size_t chosenOptionIdx = static_cast(-1); - - // Find the boundary edge that tracks closest clockwise to keep empty space on the right - for (size_t i = 0; i < options.size(); ++i) - { - if (edgeVisited[options[i].edgeId]) - continue; - - double delta = incomingAngle - options[i].angle; - if (delta < 0) - delta += 2.0 * std::numbers::pi; - - if (delta < bestDelta) - { - bestDelta = delta; - chosenOptionIdx = i; - } - } - - // If no unvisited edges are found, break out safely - if (chosenOptionIdx == static_cast(-1)) - break; - - const auto& chosen = options[chosenOptionIdx]; - edgeVisited[chosen.edgeId] = true; - - // Step forward - prevNode = currentNode; - currentNode = chosen.neighbourNode; - incomingAngle = chosen.angle; - } - if (currentLoop.size() >= 3) - { - currentLoop.push_back(currentLoop.front()); - allLoops.push_back(currentLoop); - } - } - - return allLoops; -} - -// void findBoundaryLoops(const std::vector& nodes, -// const std::vector& edges, -// std::vector>& allLoops, // Out: The boundary node sequences -// std::vector& firstTouchedFaces) // Out: The associated face IDs -// { -// // Clear outputs to ensure clean state -// allLoops.clear(); -// firstTouchedFaces.clear(); - -// std::unordered_map> boundaryAdjacency; -// std::unordered_map edgeVisited; - -// // 1. Collect all boundary edges -// for (const auto& edge : edges) -// { -// if (edge.rightFace == nullValue) -// { -// edgeVisited[edge.id] = false; - -// double dx = nodes[edge.endNode].x - nodes[edge.startNode].x; -// double dy = nodes[edge.endNode].y - nodes[edge.startNode].y; - -// double angle_start_to_end = normalizeAngle(std::atan2(dy, dx)); -// double angle_end_to_start = normalizeAngle(std::atan2(-dy, -dx)); - -// boundaryAdjacency[edge.startNode].push_back({edge.id, edge.endNode, 0, angle_start_to_end}); -// boundaryAdjacency[edge.endNode].push_back({edge.id, edge.startNode, 0, angle_end_to_start}); -// } -// } - -// // 2. Sort outgoing edges counter-clockwise -// for (auto& [node_id, connected_edges] : boundaryAdjacency) -// { -// std::sort(connected_edges.begin(), connected_edges.end(), [](const BoundaryEdge& a, const BoundaryEdge& b) -// { return a.angle < b.angle; }); -// } - -// // 3. Trace loops -// for (const auto& edge : edges) -// { -// if (edge.rightFace != nullValue || edgeVisited[edge.id]) -// continue; - -// // Capture the face ID of the first edge that starts this loop -// meshkernel::UInt associatedFaceId = edge.leftFace; - -// std::vector currentNodes; -// meshkernel::UInt prevNode = edge.startNode; -// meshkernel::UInt currentNode = edge.endNode; - -// currentNodes.push_back(prevNode); -// edgeVisited[edge.id] = true; - -// double dx = nodes[currentNode].x - nodes[prevNode].x; -// double dy = nodes[currentNode].y - nodes[prevNode].y; -// double incomingAngle = normalizeAngle(std::atan2(dy, dx)); - -// while (currentNode != edge.startNode) -// { -// currentNodes.push_back(currentNode); - -// const auto& options = boundaryAdjacency[currentNode]; -// double bestDelta = 3.0 * std::numbers::pi; -// size_t chosenOptionIdx = static_cast(-1); - -// for (size_t i = 0; i < options.size(); ++i) -// { -// if (edgeVisited[options[i].edgeId]) -// continue; - -// double delta = incomingAngle - options[i].angle; -// if (delta < 0) -// delta += 2.0 * std::numbers::pi; - -// if (delta < bestDelta) -// { -// bestDelta = delta; -// chosenOptionIdx = i; -// } -// } - -// if (chosenOptionIdx == static_cast(-1)) -// break; - -// const auto& chosen = options[chosenOptionIdx]; -// edgeVisited[chosen.edgeId] = true; - -// prevNode = currentNode; -// currentNode = chosen.neighbourNode; -// incomingAngle = chosen.angle; -// } - -// if (currentNodes.size() >= 3) -// { -// allLoops.push_back(std::move(currentNodes)); -// firstTouchedFaces.push_back(associatedFaceId); -// } -// } -// } - -void findBoundaryLoops(const std::vector& nodes, - const std::vector& edges, - std::vector>& allLoops, // Out: The boundary node sequences - std::vector>& allTouchedFaces) // Out: All faces touched per loop in step order -{ - // Clear outputs to ensure clean state - allLoops.clear(); - allTouchedFaces.clear(); - - // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle - std::unordered_map> boundaryAdjacency; - - // May be better to use meshkernel::Boolean - std::vector edgeVisited(edges.size(), false); - - // Collect all boundary edges and compute the angle - for (const auto& edge : edges) - { - if (edge.rightFace == constants::missing::uintValue) - { - double dx = nodes[edge.endNode].x - nodes[edge.startNode].x; - double dy = nodes[edge.endNode].y - nodes[edge.startNode].y; - - // Compute angle for "pinch point" sorting - // Pinch points are when the boundary polygon intersects itself, the points will have identical values. - double angleStartToEnd = normalizeAngle(std::atan2(dy, dx)); - double angleEndToStart = normalizeAngle(std::atan2(-dy, -dx)); - - // Note: edge.leftFace is the internal mesh face valid for both directions of boundary traversal - boundaryAdjacency[edge.startNode].push_back({edge.id, edge.endNode, edge.leftFace, angleStartToEnd}); - boundaryAdjacency[edge.endNode].push_back({edge.id, edge.startNode, edge.leftFace, angleEndToStart}); - } - } - - // Sort outgoing edges anti-clockwise - for (auto& [node_id, connected_edges] : boundaryAdjacency) - { - std::sort(connected_edges.begin(), connected_edges.end(), [](const BoundaryEdge& a, const BoundaryEdge& b) - { return a.angle < b.angle; }); - } - - // Trace boundary polygons - for (const auto& edge : edges) - { - if (edge.rightFace != constants::missing::uintValue || edgeVisited[edge.id]) - { - continue; - } - - std::vector currentNodes; - std::vector currentFaces; - - meshkernel::UInt prevNode = edge.startNode; - meshkernel::UInt currentNode = edge.endNode; - - currentNodes.push_back(prevNode); - currentFaces.push_back(edge.leftFace); // First face touched by the initial edge - edgeVisited[edge.id] = true; // Mark edge as having been visited - - double dx = nodes[currentNode].x - nodes[prevNode].x; - double dy = nodes[currentNode].y - nodes[prevNode].y; - double incomingAngle = normalizeAngle(std::atan2(dy, dx)); - - // Loop until we find the start node, making a closed boundary polygon - while (currentNode != edge.startNode) - { - currentNodes.push_back(currentNode); - - const std::vector& boundaryEdges = boundaryAdjacency[currentNode]; - - // The smallest visited angle. - // All edge angles must be in interval [-2pi, 2pi], - // So 3pi is used to indicate a suitable edge angle has been found - double bestDelta = 3.0 * std::numbers::pi; - meshkernel::UInt chosenEdgeIndex = constants::missing::uintValue; - - // Find the boundary edge that tracks closest clockwise to keep empty space on the right - for (size_t i = 0; i < boundaryEdges.size(); ++i) - { - if (edgeVisited[boundaryEdges[i].edgeId]) - { - continue; - } - - double delta = normalizeAngle(incomingAngle - boundaryEdges[i].angle); - - if (delta < bestDelta) - { - bestDelta = delta; - chosenEdgeIndex = i; - } - } - - // No unvisited edges were found - // So eigher the boundary loop was fully completed or a dead-end reached. - if (chosenEdgeIndex == constants::missing::uintValue) - { - break; - } - - const BoundaryEdge& chosenEdge = boundaryEdges[chosenEdgeIndex]; - edgeVisited[chosenEdge.edgeId] = true; - - // Record the face touched by this next step in the loop - currentFaces.push_back(chosenEdge.leftFace); - - prevNode = currentNode; - currentNode = chosenEdge.neighbourNode; - incomingAngle = chosenEdge.angle; - } - - if (currentNodes.size() >= 3) - { - allLoops.push_back(std::move(currentNodes)); - allTouchedFaces.push_back(std::move(currentFaces)); - } - } -} - -// in: nodes; -// in edges -// in edgesFaces -// Out: The boundary node sequences -// Out: All faces touched per loop in step order -void findBoundaryLoops(const std::vector& nodes, - const std::vector& edges, - const std::vector>& edgesFaces, - std::vector>& allLoops, - std::vector>& allTouchedFaces) -{ - allLoops.clear(); - allTouchedFaces.clear(); - - // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle - std::unordered_map> boundaryAdjacency; - - // May be better to use meshkernel::Boolean - std::vector edgeVisited(edges.size(), false); - - // Collect all boundary edges and compute the angle - for (meshkernel::UInt count = 0; count < edges.size(); ++count) - // for (meshkernel::UInt count = 0; const auto& edge : edges) - { - // const meshkernel::Edge& edge = edges[count]; - - if (edgesFaces[count][1] == constants::missing::uintValue) - { - double dx = nodes[edges[count].second].x - nodes[edges[count].first].x; - double dy = nodes[edges[count].second].y - nodes[edges[count].first].y; - - // Compute angle for "pinch point" sorting - // Pinch points are when the boundary polygon intersects itself, the points will have identical values. - double angleStartToEnd = normalizeAngle(std::atan2(dy, dx)); - double angleEndToStart = normalizeAngle(std::atan2(-dy, -dx)); - - // Note: edgesFaces[count][0] is the internal mesh face valid for both directions of boundary traversal - boundaryAdjacency[edges[count].first].push_back({count, edges[count].second, edgesFaces[count][0], angleStartToEnd}); - boundaryAdjacency[edges[count].second].push_back({count, edges[count].first, edgesFaces[count][0], angleEndToStart}); - } - } - - auto boundaryAngleLessThan = [](const BoundaryEdge& a, const BoundaryEdge& b) - { return a.angle < b.angle; }; - - // Sort outgoing edges anti-clockwise - for (auto& [node_id, connected_edges] : boundaryAdjacency) - { - std::sort(connected_edges.begin(), connected_edges.end(), boundaryAngleLessThan); - } - - // Trace boundary polygons - for (meshkernel::UInt count = 0; count < edges.size(); ++count) - { - if (edgesFaces[count][1] != constants::missing::uintValue || edgeVisited[count]) - { - continue; - } - - std::vector currentNodes; - std::vector currentFaces; - - meshkernel::UInt prevNode = edges[count].first; - meshkernel::UInt currentNode = edges[count].second; - - currentNodes.push_back(prevNode); - currentFaces.push_back(edgesFaces[count][0]); // First face touched by the initial edge - edgeVisited[count] = true; // Mark edge as having been visited - - double dx = nodes[currentNode].x - nodes[prevNode].x; - double dy = nodes[currentNode].y - nodes[prevNode].y; - double incomingAngle = normalizeAngle(std::atan2(dy, dx)); - - // Loop until we find the start node, making a closed boundary polygon - while (currentNode != edges[count].first) - { - currentNodes.push_back(currentNode); - - const std::vector& boundaryEdges = boundaryAdjacency[currentNode]; - - // The smallest visited angle. - // All edge angles must be in interval [-2pi, 2pi], - double deltaAngle = 3.0 * std::numbers::pi; - meshkernel::UInt chosenEdgeIndex = constants::missing::uintValue; - - // Find the boundary edge that tracks closest clockwise to keep empty space on the right - for (size_t e = 0; e < boundaryEdges.size(); ++e) - { - if (edgeVisited[boundaryEdges[e].edgeId]) - { - continue; - } - - double delta = normalizeAngle(incomingAngle - boundaryEdges[e].angle); - - if (delta < deltaAngle) - { - deltaAngle = delta; - chosenEdgeIndex = e; - } - } - - // No unvisited edges were found - // So eigher the boundary loop was completed or a dead-end reached. - if (chosenEdgeIndex == constants::missing::uintValue) - { - break; - } - - const BoundaryEdge& chosenEdge = boundaryEdges[chosenEdgeIndex]; - edgeVisited[chosenEdge.edgeId] = true; - - // Record the face touched by this next step in the loop - currentFaces.push_back(chosenEdge.leftFace); - - prevNode = currentNode; - currentNode = chosenEdge.neighbourNode; - incomingAngle = chosenEdge.angle; - } - - if (currentNodes.size() >= 3) - { - currentNodes.push_back(currentNodes.front()); - allLoops.push_back(std::move(currentNodes)); - allTouchedFaces.push_back(std::move(currentFaces)); - } - } -} - -void findBoundaryLoops(const std::vector& nodes, - const std::vector& edges, - const std::vector>& edgesFaces, - const meshkernel::Polygons& clippingPolygon, - std::vector>& allLoops, - std::vector>& allTouchedFaces) -{ - allLoops.clear(); - allTouchedFaces.clear(); - - // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle - std::unordered_map> boundaryAdjacency; - - // May be better to use meshkernel::Boolean - std::vector edgeVisited(edges.size(), false); - - // Collect all boundary edges and compute the angle - for (meshkernel::UInt count = 0; count < edges.size(); ++count) - // for (meshkernel::UInt count = 0; const auto& edge : edges) - { - // const meshkernel::Edge& edge = edges[count]; - - if (edgesFaces[count][1] == constants::missing::uintValue) - { - double dx = nodes[edges[count].second].x - nodes[edges[count].first].x; - double dy = nodes[edges[count].second].y - nodes[edges[count].first].y; - - // Compute angle for "pinch point" sorting - // Pinch points are when the boundary polygon intersects itself, the points will have identical values. - double angleStartToEnd = normalizeAngle(std::atan2(dy, dx)); - double angleEndToStart = normalizeAngle(std::atan2(-dy, -dx)); - - // Note: edgesFaces[count][0] is the internal mesh face valid for both directions of boundary traversal - boundaryAdjacency[edges[count].first].push_back({count, edges[count].second, edgesFaces[count][0], angleStartToEnd}); - boundaryAdjacency[edges[count].second].push_back({count, edges[count].first, edgesFaces[count][0], angleEndToStart}); - } - } - - auto boundaryAngleLessThan = [](const BoundaryEdge& a, const BoundaryEdge& b) - { return a.angle < b.angle; }; - - // Sort outgoing edges anti-clockwise - for (auto& [node_id, connected_edges] : boundaryAdjacency) - { - std::sort(connected_edges.begin(), connected_edges.end(), boundaryAngleLessThan); - } - - // Trace boundary polygons - for (meshkernel::UInt count = 0; count < edges.size(); ++count) - { - if (edgesFaces[count][1] != constants::missing::uintValue || edgeVisited[count]) - { - continue; - } - - std::vector currentNodes; - std::vector currentFaces; - - meshkernel::UInt prevNode = edges[count].first; - meshkernel::UInt currentNode = edges[count].second; - - currentNodes.push_back(prevNode); - currentFaces.push_back(edgesFaces[count][0]); // First face touched by the initial edge - edgeVisited[count] = true; // Mark edge as having been visited - - double dx = nodes[currentNode].x - nodes[prevNode].x; - double dy = nodes[currentNode].y - nodes[prevNode].y; - double incomingAngle = normalizeAngle(std::atan2(dy, dx)); - - // Loop until we find the start node, making a closed boundary polygon - while (currentNode != edges[count].first) - { - currentNodes.push_back(currentNode); - - const std::vector& boundaryEdges = boundaryAdjacency[currentNode]; - - // The smallest visited angle. - // All edge angles must be in interval [-2pi, 2pi], - double deltaAngle = 3.0 * std::numbers::pi; - meshkernel::UInt chosenEdgeIndex = constants::missing::uintValue; - - // Find the boundary edge that tracks closest clockwise to keep empty space on the right - for (size_t e = 0; e < boundaryEdges.size(); ++e) - { - if (edgeVisited[boundaryEdges[e].edgeId]) - { - continue; - } - - double delta = normalizeAngle(incomingAngle - boundaryEdges[e].angle); - - if (delta < deltaAngle) - { - deltaAngle = delta; - chosenEdgeIndex = e; - } - } - - // No unvisited edges were found - // So eigher the boundary loop was completed or a dead-end reached. - if (chosenEdgeIndex == constants::missing::uintValue) - { - break; - } - - const BoundaryEdge& chosenEdge = boundaryEdges[chosenEdgeIndex]; - edgeVisited[chosenEdge.edgeId] = true; - - // Record the face touched by this next step in the loop - currentFaces.push_back(chosenEdge.leftFace); - - prevNode = currentNode; - currentNode = chosenEdge.neighbourNode; - incomingAngle = chosenEdge.angle; - } - - if (currentNodes.size() >= 3) - { - std::vector activeSegment; - std::vector activeFaces; - - for (size_t i = 0; i < currentNodes.size(); ++i) - { - - // Check if this specific node falls inside the user's bounding box/polygon - meshkernel::Point aNode = nodes[currentNodes[i]]; - auto [isIn, whichPoly] = clippingPolygon.IsPointInPolygons(aNode); - - if (isIn) - { - activeSegment.push_back(currentNodes[i]); - - if (i < currentFaces.size()) - { - activeFaces.push_back(currentFaces[i]); - } - } - else - { - // Node fell outside, push the collected segment if it holds meaningful data - if (activeSegment.size() >= 2) - { - std::cout << " activeSegment size: " << activeSegment.size() << std::endl; - activeSegment.push_back(activeSegment.front()); - allLoops.push_back(activeSegment); - allTouchedFaces.push_back(activeFaces); - } - - activeSegment.clear(); - activeFaces.clear(); - } - } - - // Push any remaining segment if the loop ended while still inside the polygon - if (activeSegment.size() >= 2) - { - std::cout << " activeSegment size: " << activeSegment.size() << std::endl; - activeSegment.push_back(activeSegment.front()); - // activeFaces.push_back(activeFaces.front()); - allLoops.push_back(std::move(activeSegment)); - allTouchedFaces.push_back(std::move(activeFaces)); - } - } - } -} - TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygons) { @@ -3763,10 +3015,6 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon auto node52 = mesh.FindNodeCloseToAPoint({105.0, 15.0}, 1.0e-5); auto node53 = mesh.FindNodeCloseToAPoint({95.0, 15.0}, 1.0e-5); - auto node61 = mesh.FindNodeCloseToAPoint({120.0, 0.0}, 1.0e-5); - auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); - auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); - std::vector boundaryNodes; std::vector elementNodes1{{constants::missing::doubleValue, constants::missing::doubleValue}, @@ -3802,19 +3050,12 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon mesh.Node(node53), mesh.Node(node51)}; - std::vector elementNodes6{{constants::missing::doubleValue, constants::missing::doubleValue}, - mesh.Node(node61), - mesh.Node(node62), - mesh.Node(node63), - mesh.Node(node61)}; - // Combine all nodes to form a sequence of polygons boundaryNodes.insert(boundaryNodes.end(), elementNodes1.begin(), elementNodes1.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes2.begin(), elementNodes2.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes3.begin(), elementNodes3.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes4.begin(), elementNodes4.end()); boundaryNodes.insert(boundaryNodes.end(), elementNodes5.begin(), elementNodes5.end()); - boundaryNodes.insert(boundaryNodes.end(), elementNodes6.begin(), elementNodes6.end()); Polygons boundaryWithMissingElements(boundaryNodes, Projection::cartesian); @@ -3827,174 +3068,33 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon // This should not fill in the holes in the mesh mesh2.Administrate(); - meshkernel::SaveVtk(mesh2.Nodes(), mesh2.m_facesNodes, "mesh10.vtu"); - - std::vector edges2(mesh2.Edges().size()); - - for (meshkernel::UInt i = 0; i < mesh2.Edges().size(); ++i) - { - edges2[i] = Edge2{i, mesh2.Edges()[i].first, mesh2.Edges()[i].second, mesh2.m_edgesFaces[i][0], mesh2.m_edgesFaces[i][1]}; - } - - auto bLoops = findBoundaryLoops(mesh2.Nodes(), edges2, mesh2.m_edgesFaces); - - [[maybe_unused]] const auto& meshNodes = mesh2.Nodes(); - - // for (size_t i = 0; i < bLoops.size(); ++i) - for (const auto& loop : bLoops) - { - std::cout << "boudnary loop: "; - - for (size_t j = 0; j < loop.size(); ++j) - { - std::cout << "{" << loop[j].x << ", " << loop[j].y << "}, "; - } - - std::cout << std::endl; - } - - std::vector> allLoops; - std::vector> allTouchedFaces; - - std::vector clippingpolygonNodes{{-5.0, -5.0}, {140.0, 0.0}, {140.0, 35.0}, {-5.0, 35.0}, {-5.0, -5.0}}; - meshkernel::Polygons clippingPolygon(clippingpolygonNodes, mesh2.m_projection); - - findBoundaryLoops(mesh2.Nodes(), mesh2.Edges(), mesh2.m_edgesFaces, allLoops, allTouchedFaces); - // findBoundaryLoops(mesh2.Nodes(), mesh2.Edges(), mesh2.m_edgesFaces, clippingPolygon, allLoops, allTouchedFaces); - - std::cout << std::endl; - std::cout << "--------------------------------" << std::endl; - - for (meshkernel::UInt count = 0; const auto& loop : bLoops) - { - std::cout << "boudnary loop: " << " " << count << " " << loop.size() << ": " << ": "; - // std::cout << "boudnary loop: " << allTouchedFaces[count].size() << " " << loop.size() << ": " << allTouchedFaces[count][0] << ": "; - - auto boundaryPolygons = splitMultiplePolygons(loop); - - std::cout << "boudnary loop: " << loop.size() << " " << boundaryPolygons.size() << " "; - - for (size_t j = 0; j < loop.size(); ++j) - { - std::cout << "{" << loop[j].x << ", " << loop[j].y << "}, "; - } - - std::cout << std::endl; + // After having computed an administrate, the deleted elememnts will have ebeen found and removed again + // by the illegal cells polygons. - for (size_t j = 0; j < boundaryPolygons.size(); ++j) - { - std::cout << " sub-loop " << boundaryPolygons[j].size() << " "; - - for (size_t k = 0; k < boundaryPolygons[j].size(); ++k) - { - std::cout << "{" << boundaryPolygons[j][k].x << ", " << boundaryPolygons[j][k].y << "}, "; - } - - std::cout << std::endl; - } - - std::cout << std::endl; - ++count; - } - - std::cout << std::endl; - std::cout << "----------- clipping ---------------------" << std::endl; - - for (meshkernel::UInt count = 0; const auto& loop : allLoops) - { - - std::vector loopPnts(loop.size()); - - for (size_t i = 0; i < loopPnts.size(); ++i) - { - loopPnts[i] = meshNodes[loop[i]]; - } - - // std::tuple>, std::vector> splitMultiplePolygons(std::span boundary, std::span elementIds); - auto [boundaryPolygons, firstTouchedFaces] = splitMultiplePolygons(loopPnts, allTouchedFaces[count]); - - std::cout << "boudnary loop: " << firstTouchedFaces.size() << " " << loop.size() << " " << boundaryPolygons.size() << " "; - - for (size_t j = 0; j < loop.size(); ++j) - { - std::cout << "{" << loopPnts[j].x << ", " << loopPnts[j].y << "}, "; - } - - std::cout << std::endl; - std::cout << "face ids: "; - - for (size_t j = 0; j < allTouchedFaces[count].size(); ++j) - { - std::cout << allTouchedFaces[count][j] << ", "; - } - - std::cout << std::endl; - - for (size_t j = 0; j < boundaryPolygons.size(); ++j) - { - Polygons pgs(boundaryPolygons[j], mesh2.m_projection); - std::cout << " sub-loop " << boundaryPolygons[j].size() << " " << firstTouchedFaces[j] << " " << std::boolalpha << pgs.IsPointInAnyPolygon(mesh2.m_facesMassCenters[firstTouchedFaces[j]]) << ": "; - - for (size_t k = 0; k < boundaryPolygons[j].size(); ++k) - { - std::cout << "{" << boundaryPolygons[j][k].x << ", " << boundaryPolygons[j][k].y << "}, "; - } - - std::cout << std::endl; - } - - std::cout << std::endl; - ++count; - } - - // Get interior boundary polygon points - std::vector innerBoundaryPoints = mesh2.GetInnerBoundaryPolygons(); - std::vector outerBoundaryPoints = mesh2.ComputeBoundaryPolygons(std::vector()); - - std::cout << std::endl; - std::cout << "--------------------------------" << std::endl; - std::cout << "boudnary loop: "; + auto node61 = mesh.FindNodeCloseToAPoint({120.0, 0.0}, 1.0e-5); + auto node62 = mesh.FindNodeCloseToAPoint({125.0, 15.0}, 1.0e-5); + auto node63 = mesh.FindNodeCloseToAPoint({115.0, 15.0}, 1.0e-5); - for (size_t j = 0; j < outerBoundaryPoints.size(); ++j) - { - if (outerBoundaryPoints[j].IsValid()) - { - std::cout << "{" << outerBoundaryPoints[j].x << ", " << outerBoundaryPoints[j].y << "}, "; - } - else - { - std::cout << std::endl; - std::cout << "boundary loop"; - } - } + std::vector elementNodes6{{constants::missing::doubleValue, constants::missing::doubleValue}, + mesh.Node(node61), + mesh.Node(node62), + mesh.Node(node63), + mesh.Node(node61)}; - std::cout << std::endl; + // Now delete another cell and compute another administrate. + // The original deletd cells should remain deleted, and with the addiitonal deleted cell. + Polygons boundaryWithMissingElements2(elementNodes6, Projection::cartesian); + auto deleteMeshFacesUndoAction2 = mesh2.DeleteMeshFacesInPolygon(boundaryWithMissingElements2); + mesh2.Administrate(); MeshBoundaryExtractor extractor; - auto allBoundaryPoints = extractor.ExtractAll(mesh2); - - Polygons allBoundaryPolygons(allBoundaryPoints, mesh2.m_projection); - - for (UInt i = 0; i < allBoundaryPolygons.GetNumPolygons(); ++i) - { - const auto& enclosure = allBoundaryPolygons.Enclosure(i); - const auto& enclosurePoints = enclosure.Outer().Nodes(); - - std::cout << "enclosure points: "; - - for (size_t j = 0; j < enclosurePoints.size(); ++j) - { - std::cout << "{" << enclosurePoints[j].x << ", " << enclosurePoints[j].y << "}, "; - } - - std::cout << std::endl; - } + auto interiorBoundaryPoints = extractor.ExtractConcatenated(mesh2, meshkernel::BoundarySelection::InteriorOnly); // The expected number of points, should not include any land boundary points constexpr UInt expectedNumberOfNodes = 26; - ASSERT_EQ(expectedNumberOfNodes, innerBoundaryPoints.size()); + ASSERT_EQ(expectedNumberOfNodes, interiorBoundaryPoints.size()); // The edge of one of the deleted elements lies on the boundary, so will be not be part of the // interior set of polygons @@ -4006,7 +3106,7 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon constants::missing::doubleValue, 120.0, 125.0, 115.0, 120.0, constants::missing::doubleValue, - 125.0, 125.0, 135.0, 135.0, 125.0}; + 125.0, 135.0, 135.0, 125.0, 125.0}; std::vector expectedYPoints{15.0, 0.0, 15.0, 15.0, constants::missing::doubleValue, @@ -4016,11 +3116,11 @@ TEST(MeshRefinement, MeshWithHole_ShouldConstructMeshWithInteriorBoundaryPolygon constants::missing::doubleValue, 0.0, 15.0, 15.0, 0.0, constants::missing::doubleValue, - 125.0, 135.0, 135.0, 125.0, 125.0}; + 125.0, 125.0, 135.0, 135.0, 125.0}; for (size_t i = 0; i < expectedXPoints.size(); ++i) { - EXPECT_EQ(expectedXPoints[i], innerBoundaryPoints[i].x); - EXPECT_EQ(expectedYPoints[i], innerBoundaryPoints[i].y); + EXPECT_EQ(expectedXPoints[i], interiorBoundaryPoints[i].x); + EXPECT_EQ(expectedYPoints[i], interiorBoundaryPoints[i].y); } } diff --git a/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp index b00ad89fa..5e799a062 100644 --- a/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp +++ b/libs/MeshKernelApi/tests/src/InvalidCellsPolygonsTests.cpp @@ -136,7 +136,7 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) ASSERT_EQ(mesh2d.num_nodes, 3130); ASSERT_EQ(mesh2d.num_edges, 6383); - ASSERT_EQ(mesh2d.num_faces, 3232); + ASSERT_EQ(mesh2d.num_faces, 3248); int whichMeshkernelId = -1; bool isUndone = false; @@ -209,7 +209,7 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) meshkernel::constants::missing::doubleValue, 120.0, 125.0, 115.0, 120.0, meshkernel::constants::missing::doubleValue, - 125.0, 125.0, 135.0, 135.0, 125.0}; + 125.0, 135.0, 135.0, 125.0, 125.0}; std::vector expectedInnerY{15.0, 0.0, 15.0, 15.0, meshkernel::constants::missing::doubleValue, @@ -219,7 +219,7 @@ TEST(InvalidCellsPolygonsTests, MeshHolesAreMainainedAfterRefinement) meshkernel::constants::missing::doubleValue, 0.0, 15.0, 15.0, 0.0, meshkernel::constants::missing::doubleValue, - 125.0, 135.0, 135.0, 125.0, 125.0}; + 125.0, 125.0, 135.0, 135.0, 125.0}; int innerPolygonSize = 0; meshkernelapi::GeometryList innerPolygon; From 3134e43e614d75800831f4e44b271f4ae4b7158d Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 14:10:11 +0200 Subject: [PATCH 14/34] GRIDEDIT-2293 Refactored boundary finding --- .../MeshKernel/MeshBoundaryExtractor.hpp | 20 +++ .../include/MeshKernel/Operations.hpp | 2 +- libs/MeshKernel/src/Mesh2D.cpp | 59 +------- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 136 +++++++++--------- libs/MeshKernel/src/Operations.cpp | 7 +- libs/MeshKernel/tests/src/Mesh2DTest.cpp | 9 +- 6 files changed, 102 insertions(+), 131 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp index 416714b06..7fb4d6b67 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp @@ -64,12 +64,23 @@ namespace meshkernel /// @brief Ensure the angle lies between 0 .. 2 pi. static double NormalizeAngle(double angle); + /// @brief Find the edge that has the smallest angle [0 .. 2pi), to the incident edge + static UInt FindEdgeWithMinumumAngle(const std::vector& boundaryEdges, + const std::vector& edgeVisited, + const double incomingAngle); + /// @brief Construct mapping from node-id to all impinging edges for boundary all edges static void FindAllBoundarEdges(const std::vector& nodes, const std::vector& edges, const std::vector>& edgesFaces, std::unordered_map>& boundaryAdjacency); + static void Append(const Point& centre, + const Projection projection, + std::vector& boundaryLoop, + std::vector& isExterior, + std::vector>& separatedBoundaryLoops); + /// @brief Find boundary loops /// /// Any boundary loops found may need to be processed further as they may themselves contain sub-loops @@ -78,6 +89,15 @@ namespace meshkernel const std::vector>& edgesFaces, std::vector>& allLoops, std::vector>& allTouchedFaces); + + /// @brief Separate polygons that contains multiple sub-polygons and determine externality + /// + /// allBoundaryLoops is not const because it may be updated. + /// @note allBoundaryLoops should not be accessed after calling this function + static std::tuple>, std::vector> + SeparateAndDetermineExternality(const Mesh2D& mesh, + std::vector>& allBoundaryLoops, + const std::vector>& allTouchedFaces); }; } // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 7d6ec8d01..57ffdda79 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -635,7 +635,7 @@ namespace meshkernel /// @brief Split a single polygon line comprised of multiple polygons into separate polygons. /// /// @note This currently works only for polygons that intersect at the control points - std::tuple>, std::vector> SplitMultiplePolygons(std::span boundary, std::span elementIds); + std::tuple>, std::vector> SplitMultiplePolygons(std::span boundary, std::span elementIds); /// @brief Concaenate vectors of vectors of points to a vector of points separated by the invalid point template diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 077a224f8..30cf6759c 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1522,7 +1522,7 @@ std::tuple, std::vector> Mesh2D::GetAllBoun { const Polygon polygon(polygonNodes, m_projection); -#if 1 +#if 0 MeshBoundaryExtractor meshBoundaryExtractor; auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this); @@ -1530,26 +1530,16 @@ std::tuple, std::vector> Mesh2D::GetAllBoun if (polygonNodes.size() == 0) { - auto isTrue = [](size_t idx [[maybe_unused]]) + auto alwaysTrue = [](size_t idx [[maybe_unused]]) { return true; }; - return {ConcatenatePointVectors(boundaryPoints, isTrue), isEnclosingBoundary}; + return {ConcatenatePointVectors(boundaryPoints, alwaysTrue), isEnclosingBoundary}; } std::vector containedIsEnclosingBoundary; - auto removeExteriorPoints = [&polygon](const Point& p) - { - return !polygon.Contains(p); - }; - - for (size_t i = 0; i < boundaryPoints.size(); ++i) - { - std::erase_if(boundaryPoints[i], removeExteriorPoints); - } - auto addNonEmpty = [&boundaryPoints, &isEnclosingBoundary, &containedIsEnclosingBoundary](size_t idx) mutable { if (boundaryPoints[idx].size() > 0) @@ -1565,20 +1555,6 @@ std::tuple, std::vector> Mesh2D::GetAllBoun return {containedBoundaryPoints, containedIsEnclosingBoundary}; -#if 0 - for (size_t i = 0; i < boundaryPoints.size(); ++i) - { - - if (boundaryPoints[i].size() > 0) - { - containedBoundaryPoints.insert(containedBoundaryPoints.end(), boundaryPoints[i].begin(), boundaryPoints[i].end()); - containedIsEnclosingBoundary.push_back(isEnclosingBoundary[i]); - } - } - - return {containedBoundaryPoints, containedIsEnclosingBoundary}; -#endif - #else std::vector isVisited(GetNumEdges(), false); std::vector meshBoundaryPolygon; @@ -1652,7 +1628,6 @@ std::tuple, std::vector> Mesh2D::GetAllBoun if (IsMultiPolygon(currentPolygonSpan)) { - auto [multiBoundaryPolygonNodes, multiBoundaryElementIds] = (SplitMultiplePolygons(currentPolygonSpan, boundaryPolygonFaceId)); for (size_t i = 0; i < multiBoundaryPolygonNodes.size(); ++i) @@ -2307,31 +2282,9 @@ void Mesh2D::ReconstructInvalidCellsPolygon() return; } - std::vector innerBoundaryPolygons; - innerBoundaryPolygons.reserve(m_invalidCellPolygons.size()); - bool firstElement = true; - - for (UInt p = 0; p < boundaryPoints.size(); ++p) - { - - if (isEnclosingBoundary[p]) - { - continue; - } - - if (!firstElement) - { - innerBoundaryPolygons.push_back(Point(constants::missing::doubleValue, constants::missing::doubleValue)); - } - else - { - firstElement = false; - } - - innerBoundaryPolygons.insert(innerBoundaryPolygons.end(), boundaryPoints[p].begin(), boundaryPoints[p].end()); - } - - m_invalidCellPolygons = std::move(innerBoundaryPolygons); + auto isInteriorBoundary = [&isEnclosingBoundary](size_t i) + { return !isEnclosingBoundary[i]; }; + m_invalidCellPolygons = ConcatenatePointVectors(boundaryPoints, isInteriorBoundary); } std::unique_ptr Mesh2D::UpdateFaceInformation(const std::vector& faceIndices, const bool appendDeletedFaces) diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index 6af50b626..bab5c14a9 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -2,6 +2,7 @@ #include #include +#include #include "MeshKernel/Operations.hpp" @@ -21,54 +22,32 @@ std::vector meshkernel::MeshBoundaryExtractor::ExtractConcate }; return ConcatenatePointVectors(boundarySequences, isBoundarySelection); +} -#if 0 - bool isFirst = true; - - auto selectedBoundary = [boundaryType](bool isExterior) - { - using enum BoundarySelection; - - return (boundaryType == All) || - (boundaryType == ExteriorOnly && isExterior) || - (boundaryType == InteriorOnly && !isExterior); - }; +void meshkernel::MeshBoundaryExtractor::Append(const Point& centre, + const Projection projection, + std::vector& boundaryLoop, + std::vector& isExterior, + std::vector>& separatedBoundaryLoops) +{ - for (size_t i = 0; i < boundarySequences.size(); ++i) + if (auto [area, centreOfMass] = ComputePolygonAreaAndCentre(boundaryLoop, projection); area < 0.0) { - const std::vector& loop = boundarySequences[i]; - - if (!selectedBoundary(isExterior[i])) - { - continue; - } - - if (!isFirst) - { - allPoints.push_back({constants::missing::doubleValue, constants::missing::doubleValue}); - } - - allPoints.insert(allPoints.end(), loop.begin(), loop.end()); - isFirst = false; + std::ranges::reverse(boundaryLoop); } - return allPoints; -#endif + isExterior.push_back(IsPointInPolygonNodes(centre, boundaryLoop, projection)); + separatedBoundaryLoops.push_back(std::move(boundaryLoop)); } -std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh) +std::tuple>, std::vector> +meshkernel::MeshBoundaryExtractor::SeparateAndDetermineExternality(const Mesh2D& mesh, + std::vector>& allBoundaryLoops, + const std::vector>& allTouchedFaces) { - - std::vector> allBoundaryLoops; - std::vector> allTouchedFaces; - std::vector> separatedBoundaryLoops; std::vector isExterior; - const std::vector& meshNodes(mesh.Nodes()); - - FindBoundaryLoops(meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryLoops, allTouchedFaces); - for (size_t i = 0; i < allBoundaryLoops.size(); ++i) { @@ -82,26 +61,31 @@ std::tuple>, std::vector> meshk for (size_t i = 0; i < individualBoundaryPolygons.size(); ++i) { Point centre = mesh.m_facesMassCenters[firstElement[i]]; - - // If the area is calculated ot be less than zero, i.e. the boundary is traversed in clockwise direction and needs to be reversed - if (auto [area, centreOfMass] = ComputePolygonAreaAndCentre(individualBoundaryPolygons[i], mesh.m_projection); area < 0.0) - { - std::ranges::reverse(individualBoundaryPolygons[i]); - } - - isExterior.push_back(IsPointInPolygonNodes(centre, individualBoundaryPolygons[i], mesh.m_projection)); - separatedBoundaryLoops.push_back(std::move(individualBoundaryPolygons[i])); + Append(centre, mesh.m_projection, individualBoundaryPolygons[i], isExterior, separatedBoundaryLoops); } } else { - separatedBoundaryLoops.push_back(std::move(allBoundaryLoops[i])); + Point centre = mesh.m_facesMassCenters[allTouchedFaces[i][0]]; + Append(centre, mesh.m_projection, allBoundaryLoops[i], isExterior, separatedBoundaryLoops); } } return {separatedBoundaryLoops, isExterior}; } +std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh) +{ + + std::vector> allBoundaryLoops; + std::vector> allTouchedFaces; + + const std::vector& meshNodes(mesh.Nodes()); + + FindBoundaryLoops(meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryLoops, allTouchedFaces); + return SeparateAndDetermineExternality(mesh, allBoundaryLoops, allTouchedFaces); +} + double meshkernel::MeshBoundaryExtractor::NormalizeAngle(double angle) { while (angle < 0) @@ -156,6 +140,34 @@ void meshkernel::MeshBoundaryExtractor::FindAllBoundarEdges(const std::vector& boundaryEdges, + const std::vector& edgeVisited, + const double incomingAngle) +{ + + // Find the smallest visited angle. + // All edge angles must be in interval [0, 2pi], initialise with number greater than 2pi + double deltaAngle = 3.0 * std::numbers::pi; + UInt edgeIndex = constants::missing::uintValue; + + // Find the boundary edge that tracks closest clockwise to keep empty space on the right + for (size_t e = 0; e < boundaryEdges.size(); ++e) + { + if (edgeVisited[boundaryEdges[e].edgeId]) + { + continue; + } + + if (double delta = NormalizeAngle(incomingAngle - boundaryEdges[e].angle); delta < deltaAngle) + { + deltaAngle = delta; + edgeIndex = e; + } + } + + return edgeIndex; +} + void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector& nodes, const std::vector& edges, const std::vector>& edgesFaces, @@ -207,39 +219,18 @@ void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector& boundaryEdges = boundaryAdjacency[currentNodeIndex]; - // Find the smallest visited angle. - // All edge angles must be in interval [0, 2pi], initialise with number greater than 2pi - double deltaAngle = 3.0 * std::numbers::pi; - UInt chosenEdgeIndex = constants::missing::uintValue; - - // Find the boundary edge that tracks closest clockwise to keep empty space on the right - for (size_t e = 0; e < boundaryEdges.size(); ++e) - { - if (edgeVisited[boundaryEdges[e].edgeId]) - { - continue; - } - - double delta = NormalizeAngle(incomingAngle - boundaryEdges[e].angle); - - if (delta < deltaAngle) - { - deltaAngle = delta; - chosenEdgeIndex = e; - } - } + UInt edgeIndex = FindEdgeWithMinumumAngle(boundaryEdges, edgeVisited, incomingAngle); // No unvisited edges were found // So eigher the boundary loop was completed or a dead-end reached. - if (chosenEdgeIndex == constants::missing::uintValue) + if (edgeIndex == constants::missing::uintValue) { break; } - const BoundaryEdge& chosenEdge = boundaryEdges[chosenEdgeIndex]; + const BoundaryEdge& chosenEdge = boundaryEdges[edgeIndex]; edgeVisited[chosenEdge.edgeId] = true; - // Save the face touched by this next step in the loop currentFaces.push_back(chosenEdge.leftFace); prevNodeIndex = currentNodeIndex; @@ -249,8 +240,9 @@ void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector= 3) { + // Close the polygon currentNodes.push_back(currentNodes.front()); - allLoops.push_back(std::move(currentNodes)); + allLoops.emplace_back(std::move(currentNodes)); allTouchedFaces.push_back(std::move(currentFaces)); } } diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index fd82d83d5..b1bc75d27 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -1838,9 +1838,10 @@ namespace meshkernel std::set uniquePoints; - for (const Point& p : boundary) + // Since polygons are always closed (first = last) then start from 1 after the first + for (size_t i = 1; i < boundary.size(); ++i) { - if (!uniquePoints.insert(p).second) + if (!uniquePoints.insert(boundary[i]).second) { return true; } @@ -1849,7 +1850,7 @@ namespace meshkernel return false; } - std::tuple>, std::vector> SplitMultiplePolygons(std::span boundaryPoints, std::span elementIds) + std::tuple>, std::vector> SplitMultiplePolygons(std::span boundaryPoints, std::span elementIds) { std::vector> completedPolygons; std::vector firstElementIds; diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index 6c5bf145e..2c118afc2 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -233,6 +233,11 @@ TEST(Mesh2D, MeshBoundaryToPolygon) // 2 Execution const auto meshBoundaryPolygon = mesh.ComputeBoundaryPolygons(polygonNodes); + for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) + { + std::cout << "{" << meshBoundaryPolygon[i].x << ", " << meshBoundaryPolygon[i].y << "}" << std::endl; + } + // 3 Validation const double tolerance = 1e-5; ASSERT_NEAR(0.0, meshBoundaryPolygon[0].x, tolerance); @@ -242,9 +247,9 @@ TEST(Mesh2D, MeshBoundaryToPolygon) ASSERT_NEAR(0.0, meshBoundaryPolygon[4].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].y, tolerance); - ASSERT_NEAR(-5.0, meshBoundaryPolygon[1].y, tolerance); + ASSERT_NEAR(5.0, meshBoundaryPolygon[1].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[2].y, tolerance); - ASSERT_NEAR(5.0, meshBoundaryPolygon[3].y, tolerance); + ASSERT_NEAR(-5.0, meshBoundaryPolygon[3].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[4].y, tolerance); } From c56259ad99943a0894b3b3ee6dae4db373301f8c Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 14:18:45 +0200 Subject: [PATCH 15/34] GRIDEDIT-2293 Fixed docygen warnings --- .../MeshKernel/MeshBoundaryExtractor.hpp | 31 ++++++----- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 52 +++++++++---------- 2 files changed, 43 insertions(+), 40 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp index 7fb4d6b67..725798adc 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp @@ -55,10 +55,10 @@ namespace meshkernel /// @brief Temporary struct, used when computing the boundaries struct BoundaryEdge { - UInt edgeId; - UInt neighbourNode; - UInt leftFace; // Store face mapping on the edge for easy retrieval during loop trace - double angle; // Angle of the edge pointing away from the pivot node + UInt edgeId; ///< Id of edge + UInt neighbourNode; ///< Id of node at opposite end of edge + UInt leftFace; ///< Store face mapping on the edge for easy retrieval during polygon trace + double angle; ///< Angle of the edge pointing away from the pivot node }; /// @brief Ensure the angle lies between 0 .. 2 pi. @@ -75,28 +75,31 @@ namespace meshkernel const std::vector>& edgesFaces, std::unordered_map>& boundaryAdjacency); + /// @brief Append the boundary polygon to the set of all boundary polygon + /// + /// The boundary polygons may be reversed if found to be in clockwise direction static void Append(const Point& centre, const Projection projection, - std::vector& boundaryLoop, + std::vector& boundaryPolygon, std::vector& isExterior, - std::vector>& separatedBoundaryLoops); + std::vector>& separatedBoundaryPolygons); /// @brief Find boundary loops /// /// Any boundary loops found may need to be processed further as they may themselves contain sub-loops - static void FindBoundaryLoops(const std::vector& nodes, - const std::vector& edges, - const std::vector>& edgesFaces, - std::vector>& allLoops, - std::vector>& allTouchedFaces); + static void FindBoundaryPolygons(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::vector>& allPolygons, + std::vector>& allTouchedFaces); /// @brief Separate polygons that contains multiple sub-polygons and determine externality /// - /// allBoundaryLoops is not const because it may be updated. - /// @note allBoundaryLoops should not be accessed after calling this function + /// allBoundaryPolygons is not const because it may be updated. + /// @note allBoundaryPolygons should not be accessed after calling this function static std::tuple>, std::vector> SeparateAndDetermineExternality(const Mesh2D& mesh, - std::vector>& allBoundaryLoops, + std::vector>& allBoundaryPolygons, const std::vector>& allTouchedFaces); }; diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index bab5c14a9..b74edf731 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -26,64 +26,64 @@ std::vector meshkernel::MeshBoundaryExtractor::ExtractConcate void meshkernel::MeshBoundaryExtractor::Append(const Point& centre, const Projection projection, - std::vector& boundaryLoop, + std::vector& boundaryPolygon, std::vector& isExterior, - std::vector>& separatedBoundaryLoops) + std::vector>& separatedBoundaryPolygons) { - if (auto [area, centreOfMass] = ComputePolygonAreaAndCentre(boundaryLoop, projection); area < 0.0) + if (auto [area, centreOfMass] = ComputePolygonAreaAndCentre(boundaryPolygon, projection); area < 0.0) { - std::ranges::reverse(boundaryLoop); + std::ranges::reverse(boundaryPolygon); } - isExterior.push_back(IsPointInPolygonNodes(centre, boundaryLoop, projection)); - separatedBoundaryLoops.push_back(std::move(boundaryLoop)); + isExterior.push_back(IsPointInPolygonNodes(centre, boundaryPolygon, projection)); + separatedBoundaryPolygons.push_back(std::move(boundaryPolygon)); } std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::SeparateAndDetermineExternality(const Mesh2D& mesh, - std::vector>& allBoundaryLoops, + std::vector>& allBoundaryPolygons, const std::vector>& allTouchedFaces) { - std::vector> separatedBoundaryLoops; + std::vector> separatedBoundaryPolygons; std::vector isExterior; - for (size_t i = 0; i < allBoundaryLoops.size(); ++i) + for (size_t i = 0; i < allBoundaryPolygons.size(); ++i) { - if (IsMultiPolygon(allBoundaryLoops[i])) + if (IsMultiPolygon(allBoundaryPolygons[i])) { // It can be that some of the boundaries that are found are composed of multiple sub-boundaries. // Some of the sub-boundaries may be combined in a non conforming way. // In either case, the boundaries are separated into distinct boundary polygons. - auto [individualBoundaryPolygons, firstElement] = SplitMultiplePolygons(allBoundaryLoops[i], allTouchedFaces[i]); + auto [individualBoundaryPolygons, firstElement] = SplitMultiplePolygons(allBoundaryPolygons[i], allTouchedFaces[i]); for (size_t i = 0; i < individualBoundaryPolygons.size(); ++i) { Point centre = mesh.m_facesMassCenters[firstElement[i]]; - Append(centre, mesh.m_projection, individualBoundaryPolygons[i], isExterior, separatedBoundaryLoops); + Append(centre, mesh.m_projection, individualBoundaryPolygons[i], isExterior, separatedBoundaryPolygons); } } else { Point centre = mesh.m_facesMassCenters[allTouchedFaces[i][0]]; - Append(centre, mesh.m_projection, allBoundaryLoops[i], isExterior, separatedBoundaryLoops); + Append(centre, mesh.m_projection, allBoundaryPolygons[i], isExterior, separatedBoundaryPolygons); } } - return {separatedBoundaryLoops, isExterior}; + return {separatedBoundaryPolygons, isExterior}; } std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh) { - std::vector> allBoundaryLoops; + std::vector> allBoundaryPolygons; std::vector> allTouchedFaces; const std::vector& meshNodes(mesh.Nodes()); - FindBoundaryLoops(meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryLoops, allTouchedFaces); - return SeparateAndDetermineExternality(mesh, allBoundaryLoops, allTouchedFaces); + FindBoundaryPolygons(meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryPolygons, allTouchedFaces); + return SeparateAndDetermineExternality(mesh, allBoundaryPolygons, allTouchedFaces); } double meshkernel::MeshBoundaryExtractor::NormalizeAngle(double angle) @@ -168,13 +168,13 @@ meshkernel::UInt meshkernel::MeshBoundaryExtractor::FindEdgeWithMinumumAngle(con return edgeIndex; } -void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector& nodes, - const std::vector& edges, - const std::vector>& edgesFaces, - std::vector>& allLoops, - std::vector>& allTouchedFaces) +void meshkernel::MeshBoundaryExtractor::FindBoundaryPolygons(const std::vector& nodes, + const std::vector& edges, + const std::vector>& edgesFaces, + std::vector>& allPolygons, + std::vector>& allTouchedFaces) { - allLoops.clear(); + allPolygons.clear(); allTouchedFaces.clear(); // Mapping from mesh node-id to a sequence of boundary edges, the boundary edges will be sorted by angle @@ -212,7 +212,7 @@ void meshkernel::MeshBoundaryExtractor::FindBoundaryLoops(const std::vector Date: Tue, 4 Aug 2026 14:26:35 +0200 Subject: [PATCH 16/34] GRIDEDIT-2293 Fixed spelling errors in comments and fixed macos build --- .../include/MeshKernel/MeshBoundaryExtractor.hpp | 2 +- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 9 ++++++++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp index 725798adc..e77261457 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp @@ -47,7 +47,7 @@ namespace meshkernel /// @brief Extract all boundaries keeping them separated and /// /// The result consists of an array of each of the boundary polygons - /// Additionally, an array indicating if the boundary polygon is a exterior boudnary or not. + /// Additionally, an array indicating if the boundary polygon is a exterior boundary or not. /// True => is-exterior, False => otherwise static std::tuple>, std::vector> Extract(const Mesh2D& mesh); diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index b74edf731..4d61a7fc6 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -9,7 +9,14 @@ std::vector meshkernel::MeshBoundaryExtractor::ExtractConcatenated(const Mesh2D& mesh, BoundarySelection boundaryType) { - auto [boundarySequences, isExterior] = Extract(mesh); + // The use of std::tie instead of a structured inding is due to limitations in the macos compiler + // + // auto [boundarySequences, isExterior] = Extract(mesh); + // Replace the 3 lines below with the line above + std::vector> boundarySequences; + std::vector isExterior; + std::tie (boundarySequences, isExterior) = Extract(mesh); + std::vector allPoints; auto isBoundarySelection = [boundaryType, &isExterior](size_t idx) From fd5100e92cbe550dcfb450b497cb60a6e93cc13d Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 14:27:13 +0200 Subject: [PATCH 17/34] GRIDEDIT-2293 Fixed comment --- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index 4d61a7fc6..97a15f9a0 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -9,7 +9,7 @@ std::vector meshkernel::MeshBoundaryExtractor::ExtractConcatenated(const Mesh2D& mesh, BoundarySelection boundaryType) { - // The use of std::tie instead of a structured inding is due to limitations in the macos compiler + // The use of std::tie instead of a structured binding is due to limitations in the macos compiler // // auto [boundarySequences, isExterior] = Extract(mesh); // Replace the 3 lines below with the line above From 822ab059baab9f1048d1a1b9e55b79846d6ce2bc Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 14:27:44 +0200 Subject: [PATCH 18/34] GRIDEDIT-2293 Fixed formattng warning --- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index 97a15f9a0..886b17806 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -15,7 +15,7 @@ std::vector meshkernel::MeshBoundaryExtractor::ExtractConcate // Replace the 3 lines below with the line above std::vector> boundarySequences; std::vector isExterior; - std::tie (boundarySequences, isExterior) = Extract(mesh); + std::tie(boundarySequences, isExterior) = Extract(mesh); std::vector allPoints; From d36cefc79d4826a6af36deb7b0c07ad76ac7df5a Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 14:34:16 +0200 Subject: [PATCH 19/34] GRIDEDIT-2293 Fixed another macos error --- libs/MeshKernel/src/Mesh2D.cpp | 9 ++++++++- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 1 + 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 30cf6759c..a0738dc6a 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -2274,7 +2274,14 @@ void Mesh2D::ReconstructInvalidCellsPolygon() { MeshBoundaryExtractor meshBoundaryExtractor; - auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this); + // The use of std::tie instead of a structured binding is due to limitations in the macos compiler + // The compiler error "error: capturing a structured binding is not yet supported in OpenMP" + // + // auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this); + // Replace the 3 lines below with the line above + std::vector> boundaryPoints; + std::vector isEnclosingBoundary; + std::tie(boundaryPoints, isEnclosingBoundary) = meshBoundaryExtractor.Extract(*this); if (boundaryPoints.size() <= 1) { diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index 886b17806..9b0c96ef8 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -10,6 +10,7 @@ std::vector meshkernel::MeshBoundaryExtractor::ExtractConcate { // The use of std::tie instead of a structured binding is due to limitations in the macos compiler + // The compiler error "error: capturing a structured binding is not yet supported in OpenMP" // // auto [boundarySequences, isExterior] = Extract(mesh); // Replace the 3 lines below with the line above From 58d3d1066a3a45df590464e41ee298d56f55e6d4 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 13 Jul 2026 11:28:44 +0200 Subject: [PATCH 20/34] GRIDEDIT-2292 Fix compilation under macos --- libs/MeshKernel/tests/CMakeLists.txt | 4 ++++ libs/MeshKernelApi/tests/CMakeLists.txt | 4 ++++ 2 files changed, 8 insertions(+) diff --git a/libs/MeshKernel/tests/CMakeLists.txt b/libs/MeshKernel/tests/CMakeLists.txt index 81cca6b35..4d85973e8 100644 --- a/libs/MeshKernel/tests/CMakeLists.txt +++ b/libs/MeshKernel/tests/CMakeLists.txt @@ -67,6 +67,10 @@ set( # add sources to target target_sources(${TARGET_NAME} PRIVATE ${SRC_LIST}) +if(APPLE) + add_compile_options("-Wnocharacter-conversion") +endif () + # Should be linked to the main library, as well as the google test library target_link_libraries( ${TARGET_NAME} diff --git a/libs/MeshKernelApi/tests/CMakeLists.txt b/libs/MeshKernelApi/tests/CMakeLists.txt index f642be1ba..2d10cb031 100644 --- a/libs/MeshKernelApi/tests/CMakeLists.txt +++ b/libs/MeshKernelApi/tests/CMakeLists.txt @@ -72,6 +72,10 @@ add_test( COMMAND ${TARGET_NAME} ) +if(APPLE) + add_compile_options("-Wnocharacter-conversion") +endif () + # Copy the MeshKernel shared library to the target directory add_custom_command( TARGET ${TARGET_NAME} From f5834a938a398074fb0d75eca928140f6e59872a Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 13 Jul 2026 12:24:47 +0200 Subject: [PATCH 21/34] GRIDEDIT-2292 Fix compilation under macos --- libs/MeshKernel/tests/CMakeLists.txt | 2 +- libs/MeshKernelApi/tests/CMakeLists.txt | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/libs/MeshKernel/tests/CMakeLists.txt b/libs/MeshKernel/tests/CMakeLists.txt index 4d85973e8..6b6eecb96 100644 --- a/libs/MeshKernel/tests/CMakeLists.txt +++ b/libs/MeshKernel/tests/CMakeLists.txt @@ -68,7 +68,7 @@ set( target_sources(${TARGET_NAME} PRIVATE ${SRC_LIST}) if(APPLE) - add_compile_options("-Wnocharacter-conversion") + add_compile_options("-Wno-character-conversion") endif () # Should be linked to the main library, as well as the google test library diff --git a/libs/MeshKernelApi/tests/CMakeLists.txt b/libs/MeshKernelApi/tests/CMakeLists.txt index 2d10cb031..cb6f97a02 100644 --- a/libs/MeshKernelApi/tests/CMakeLists.txt +++ b/libs/MeshKernelApi/tests/CMakeLists.txt @@ -73,7 +73,7 @@ add_test( ) if(APPLE) - add_compile_options("-Wnocharacter-conversion") + add_compile_options("-Wno-character-conversion") endif () # Copy the MeshKernel shared library to the target directory From 5596606753e6a6d3618ab40783a0aefbe0fcc89e Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 13 Jul 2026 12:27:59 +0200 Subject: [PATCH 22/34] GRIDEDIT-2292 Fix compilation under macos --- libs/MeshKernel/tests/CMakeLists.txt | 4 ---- libs/MeshKernelApi/tests/CMakeLists.txt | 4 ---- 2 files changed, 8 deletions(-) diff --git a/libs/MeshKernel/tests/CMakeLists.txt b/libs/MeshKernel/tests/CMakeLists.txt index 6b6eecb96..81cca6b35 100644 --- a/libs/MeshKernel/tests/CMakeLists.txt +++ b/libs/MeshKernel/tests/CMakeLists.txt @@ -67,10 +67,6 @@ set( # add sources to target target_sources(${TARGET_NAME} PRIVATE ${SRC_LIST}) -if(APPLE) - add_compile_options("-Wno-character-conversion") -endif () - # Should be linked to the main library, as well as the google test library target_link_libraries( ${TARGET_NAME} diff --git a/libs/MeshKernelApi/tests/CMakeLists.txt b/libs/MeshKernelApi/tests/CMakeLists.txt index cb6f97a02..f642be1ba 100644 --- a/libs/MeshKernelApi/tests/CMakeLists.txt +++ b/libs/MeshKernelApi/tests/CMakeLists.txt @@ -72,10 +72,6 @@ add_test( COMMAND ${TARGET_NAME} ) -if(APPLE) - add_compile_options("-Wno-character-conversion") -endif () - # Copy the MeshKernel shared library to the target directory add_custom_command( TARGET ${TARGET_NAME} From ec39844f502f0bfd10975015dbc2e7754297aba5 Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 14:57:25 +0200 Subject: [PATCH 23/34] GRIDEDIT-2293 Clear interior boundary polygons array if no interior polygons were found --- libs/MeshKernel/src/Mesh2D.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index a0738dc6a..f334725a0 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -2286,6 +2286,7 @@ void Mesh2D::ReconstructInvalidCellsPolygon() if (boundaryPoints.size() <= 1) { // There are no interior boundary polygons + m_invalidCellPolygons.clear (); return; } From 888bd26765b963b1ee597592949555ca0366ce33 Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 14:59:11 +0200 Subject: [PATCH 24/34] GRIDEDIT-2293 Fixed formatting warning --- libs/MeshKernel/src/Mesh2D.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index f334725a0..befaf8a86 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -2286,7 +2286,7 @@ void Mesh2D::ReconstructInvalidCellsPolygon() if (boundaryPoints.size() <= 1) { // There are no interior boundary polygons - m_invalidCellPolygons.clear (); + m_invalidCellPolygons.clear(); return; } From afae0521139f4d9f5fb7b96b0dd8ec769ddea4e6 Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 15:10:59 +0200 Subject: [PATCH 25/34] GRIDEDIT-2293 Added check for invalid point when splitting multi-polygon sequences --- libs/MeshKernel/src/Operations.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index b1bc75d27..e9c0f5044 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -1861,10 +1861,15 @@ namespace meshkernel for (size_t i = 0; i < boundaryPoints.size(); ++i) { - const Point& current_point = boundaryPoints[i]; + const Point& currentPoint = boundaryPoints[i]; + + if (!currentPoint.IsValid()) + { + continue; + } // Check if this node closes a loop with a previously visited point - if (auto it = activePoints.find(current_point); it != activePoints.end()) + if (auto it = activePoints.find(currentPoint); it != activePoints.end()) { size_t loopStartIndex = it->second; @@ -1893,8 +1898,8 @@ namespace meshkernel int currentEdge = (i < elementIds.size()) ? elementIds[i] : -1; - activePoints[current_point] = pointFaceStack.size(); - pointFaceStack.emplace_back(current_point, currentEdge); + activePoints[currentPoint] = pointFaceStack.size(); + pointFaceStack.emplace_back(currentPoint, currentEdge); } return {completedPolygons, firstElementIds}; From 6309b923c7dab60682cbfef4aea2429751c76da5 Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 15:13:16 +0200 Subject: [PATCH 26/34] GRIDEDIT-2293 Fixed windows build --- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index 9b0c96ef8..b912ce7ca 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -66,10 +66,10 @@ meshkernel::MeshBoundaryExtractor::SeparateAndDetermineExternality(const Mesh2D& // In either case, the boundaries are separated into distinct boundary polygons. auto [individualBoundaryPolygons, firstElement] = SplitMultiplePolygons(allBoundaryPolygons[i], allTouchedFaces[i]); - for (size_t i = 0; i < individualBoundaryPolygons.size(); ++i) + for (size_t j = 0; j < individualBoundaryPolygons.size(); ++j) { - Point centre = mesh.m_facesMassCenters[firstElement[i]]; - Append(centre, mesh.m_projection, individualBoundaryPolygons[i], isExterior, separatedBoundaryPolygons); + Point centre = mesh.m_facesMassCenters[firstElement[j]]; + Append(centre, mesh.m_projection, individualBoundaryPolygons[j], isExterior, separatedBoundaryPolygons); } } else @@ -159,7 +159,7 @@ meshkernel::UInt meshkernel::MeshBoundaryExtractor::FindEdgeWithMinumumAngle(con UInt edgeIndex = constants::missing::uintValue; // Find the boundary edge that tracks closest clockwise to keep empty space on the right - for (size_t e = 0; e < boundaryEdges.size(); ++e) + for (UInt e = 0; e < boundaryEdges.size(); ++e) { if (edgeVisited[boundaryEdges[e].edgeId]) { From f56976216b2fa268d37c1c1508e41a5bc95d7321 Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 15:27:19 +0200 Subject: [PATCH 27/34] GRIDEDIT-2293 Removed cout in unittest, added copyright header to implementation file --- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 27 +++++++++++++++++++ libs/MeshKernel/tests/src/Mesh2DTest.cpp | 20 -------------- 2 files changed, 27 insertions(+), 20 deletions(-) diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index b912ce7ca..976ffcda9 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -1,3 +1,30 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2026. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program 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 General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + #include "MeshKernel/MeshBoundaryExtractor.hpp" #include diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index 2c118afc2..91a9a04fc 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -233,11 +233,6 @@ TEST(Mesh2D, MeshBoundaryToPolygon) // 2 Execution const auto meshBoundaryPolygon = mesh.ComputeBoundaryPolygons(polygonNodes); - for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) - { - std::cout << "{" << meshBoundaryPolygon[i].x << ", " << meshBoundaryPolygon[i].y << "}" << std::endl; - } - // 3 Validation const double tolerance = 1e-5; ASSERT_NEAR(0.0, meshBoundaryPolygon[0].x, tolerance); @@ -314,21 +309,6 @@ TEST(Mesh2D, MeshBoundaryToPolygonWithSelection) meshkernel::Print(mesh.Nodes(), mesh.Edges()); - std::cout << std::endl; - - for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) - { - std::cout << meshBoundaryPolygon[i].x << ", "; - } - - std::cout << std::endl; - - for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) - { - std::cout << meshBoundaryPolygon[i].y << ", "; - } - - std::cout << std::endl; // 3 Validation const double tolerance = 1e-5; From 0e0af25bf676e20fe438baff20d8282bab1d30fb Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 16:19:52 +0200 Subject: [PATCH 28/34] GRIDEDIT-2293 Replaced function with new code for extracting boundaries --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 17 -- libs/MeshKernel/src/Mesh2D.cpp | 234 +----------------- .../tests/src/MeshRefinementTests.cpp | 19 +- 3 files changed, 12 insertions(+), 258 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index c8ce399c4..b4e283ee4 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -254,10 +254,6 @@ namespace meshkernel /// @return The resulting polygon mesh boundary [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon); - /// @brief Convert all mesh boundaries to a vector of polygon nodes - /// @return The resulting set of polygons, describing interior mesh boundaries - std::vector ComputeInnerBoundaryPolygons() const; - /// @brief Gets the hanging edges /// @return A vector with the indices of the hanging edges [[nodiscard]] std::vector GetHangingEdges() const; @@ -451,16 +447,6 @@ namespace meshkernel std::vector& meshBoundaryPolygon, std::vector& boundaryPolygonFaceId) const; - /// @brief Constructs a polygon or polygons from the meshboundary, by walking through the mesh - /// - /// If there are multiple polygons connected by a single node, then these will be separated into individual polygons - void WalkMultiBoundaryFromNode(std::vector& edgeIsVisited, - std::vector& nodeIsVisited, - UInt& currentNode, - std::vector& meshBoundaryPolygon, - std::vector& nodeIds, - std::vector& subSequence, - std::vector& illegalCells) const; /// @brief Reconstruct the invalid cell polygons /// @@ -474,9 +460,6 @@ namespace meshkernel /// polygonal sub-sequence forms a external boundary [[nodiscard]] std::tuple, std::vector> GetAllBoundaryPolygons(const std::vector& polygon); - /// @brief Ensure that all polynomials are orientated in the ACW direction. - void OrientatePolygonsAntiClockwise(std::vector& polygonNodes) const; - /// @brief Removes the outer domain boundary polygon from the set of polygons /// /// It is assumed that the outer domain polygon contains the most nodes diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index befaf8a86..24e368482 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -218,7 +218,8 @@ void Mesh2D::DoAdministrationGivenFaceNodesMapping(const std::vector, std::vector> Mesh2D::GetAllBoun #endif } -std::vector Mesh2D::ComputeInnerBoundaryPolygons() const -{ - if (GetNumFaces() == 0) - { - return std::vector(); - } - - std::vector illegalCells; - illegalCells.reserve(GetNumNodes()); - std::vector meshBoundaryPolygon; - meshBoundaryPolygon.reserve(GetNumNodes()); - std::vector subSequence; - subSequence.reserve(GetNumNodes()); - - std::vector edgeIsVisited(GetNumEdges(), false); - std::vector nodeIsVisited(GetNumNodes(), false); - - std::vector nodeIds; - nodeIds.reserve(GetNumNodes()); - - for (UInt e = 0; e < GetNumEdges(); e++) - { - if (edgeIsVisited[e] || !IsEdgeOnBoundary(e)) - { - continue; - } - - const auto firstNodeIndex = m_edges[e].first; - const auto secondNodeIndex = m_edges[e].second; - const auto firstNode = m_nodes[firstNodeIndex]; - const auto secondNode = m_nodes[secondNodeIndex]; - - // Start a new polyline - if (!subSequence.empty()) - { - subSequence.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); - nodeIds.emplace_back(constants::missing::uintValue); - } - - // Put the current edge on the mesh boundary, mark it as visited - const auto startPolygonEdges = static_cast(subSequence.size()); - subSequence.emplace_back(firstNode); - subSequence.emplace_back(secondNode); - nodeIds.emplace_back(firstNodeIndex); - nodeIds.emplace_back(secondNodeIndex); - edgeIsVisited[e] = true; - nodeIsVisited[firstNodeIndex] = true; - nodeIsVisited[secondNodeIndex] = true; - - // walk the current mesh boundary - auto currentNode = secondNodeIndex; - WalkMultiBoundaryFromNode(edgeIsVisited, nodeIsVisited, currentNode, subSequence, nodeIds, meshBoundaryPolygon, illegalCells); - - const auto numNodesFirstTail = static_cast(subSequence.size()); - - // if the boundary polygon is not closed - if (currentNode != firstNodeIndex) - { - // Now grow a polyline starting at the other side of the original link L, i.e., the second tail - currentNode = firstNodeIndex; - WalkMultiBoundaryFromNode(edgeIsVisited, nodeIsVisited, currentNode, subSequence, nodeIds, meshBoundaryPolygon, illegalCells); - } - - // There is a nonempty second tail, so reverse the first tail, so that they connect. - if (subSequence.size() > numNodesFirstTail) - { - const auto start = startPolygonEdges + static_cast(std::ceil((numNodesFirstTail - startPolygonEdges + static_cast(1)) * 0.5)); - - for (auto n = start; n < numNodesFirstTail; n++) - { - const auto backupPoint = subSequence[n]; - const auto replaceIndex = numNodesFirstTail - n + firstNodeIndex; - subSequence[n] = subSequence[replaceIndex]; - subSequence[replaceIndex] = backupPoint; - - const UInt backupPointIndex = nodeIds[n]; - nodeIds[n] = nodeIds[replaceIndex]; - nodeIds[replaceIndex] = backupPointIndex; - } - } - } - - OrientatePolygonsAntiClockwise(illegalCells); - - return illegalCells; -} - -void Mesh2D::OrientatePolygonsAntiClockwise(std::vector& polygonNodes) const -{ - UInt polygonStart = 0; - UInt polygonLength = 0; - UInt index = 0; - - while (index < polygonNodes.size()) - { - polygonStart = index; - polygonLength = 0; - - for (UInt i = polygonStart; i < polygonNodes.size(); ++i) - { - ++index; - - if (!polygonNodes[i].IsValid()) - { - polygonLength = i - polygonStart; - break; - } - - if (index == polygonNodes.size()) - { - ++polygonLength; - } - } - - if (polygonLength > 0) - { - const Point inValidPoint = {constants::missing::doubleValue, constants::missing::doubleValue}; - Point zeroPoint{0.0, 0.0}; - - Point midPoint = std::accumulate(polygonNodes.begin() + polygonStart, polygonNodes.begin() + polygonStart + polygonLength - 1, zeroPoint) / static_cast(polygonLength - 1); - - if (!IsPointInPolygonNodes(midPoint, polygonNodes, m_projection, inValidPoint, polygonStart, polygonStart + polygonLength)) - { - // reverse order of polygon nodes - if (polygonLength - 1 == 3) - { - // Only the second and third points need be swapped to reverse the points in a triangle polygon - std::swap(polygonNodes[polygonStart + 1], polygonNodes[polygonStart + 2]); - } - else if (polygonLength - 1 == 4) - { - // Only the second and fourth points need be swapped to reverse the points in a quadrilateral polygon - std::swap(polygonNodes[polygonStart + 1], polygonNodes[polygonStart + 3]); - } - else - { - std::reverse(polygonNodes.begin() + polygonStart, polygonNodes.begin() + polygonStart + polygonLength); - } - } - } - } -} - std::vector Mesh2D::RemoveOuterDomainBoundaryPolygon(const std::vector& polygonNodes) const { // Remove outer boundary. @@ -1897,94 +1755,6 @@ void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, } } -void Mesh2D::WalkMultiBoundaryFromNode(std::vector& edgeIsVisited, - std::vector& nodeIsVisited, - UInt& currentNode, - std::vector& subSequence, - std::vector& nodeIds, - std::vector& meshBoundaryPolygon, - std::vector& illegalCells) const -{ - UInt e = 0; - - while (e < m_nodesNumEdges[currentNode]) - { - const auto currentEdge = m_nodesEdges[currentNode][e]; - - if (edgeIsVisited[currentEdge] || !IsEdgeOnBoundary(currentEdge)) - { - e++; - continue; - } - - UInt nextNode = OtherNodeOfEdge(m_edges[currentEdge], currentNode); - e = 0; - - if (nodeIsVisited[nextNode]) - { - UInt lastIndex = constants::missing::uintValue; - - // Find index of last time node was added - for (size_t ii = nodeIds.size(); ii >= 1; --ii) - { - UInt i = static_cast(ii) - 1; - - if (nodeIds[i] == constants::missing::uintValue) - { - break; - } - - if (nodeIds[i] == nextNode) - { - lastIndex = i; - break; - } - } - - if (lastIndex != constants::missing::uintValue) - { - size_t start = meshBoundaryPolygon.size(); - - if (!meshBoundaryPolygon.empty()) - { - meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); - ++start; - } - - meshBoundaryPolygon.insert(meshBoundaryPolygon.end(), subSequence.begin() + lastIndex, subSequence.end()); - meshBoundaryPolygon.emplace_back(subSequence[lastIndex]); - - // the points making up the last found polygon - std::span currentPolygon(meshBoundaryPolygon.data() + start, meshBoundaryPolygon.data() + meshBoundaryPolygon.size()); - // Since the edge lies on a boundary, there will be only 1 attached element. - // This element will be in the 0th position - UInt connectedFace = m_edgesFaces[currentEdge][0]; - - if (!IsPointInPolygonNodes(m_facesMassCenters[connectedFace], currentPolygon, m_projection)) - { - // If the centre of this element does not lie within the polygon, then the polygon defines a hole in the mesh. - - if (!illegalCells.empty()) - { - // If illegal cells array is not empty then add the polygon separator - illegalCells.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); - } - - illegalCells.insert(illegalCells.end(), currentPolygon.begin(), currentPolygon.end()); - } - - subSequence.resize(lastIndex); - nodeIds.resize(lastIndex); - } - } - - currentNode = nextNode; - subSequence.emplace_back(m_nodes[currentNode]); - edgeIsVisited[currentEdge] = true; - nodeIsVisited[currentNode] = true; - nodeIds.emplace_back(currentNode); - } -} std::vector Mesh2D::GetHangingEdges() const { diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 20fbad288..b4cfdd19f 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -2891,7 +2891,8 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa auto deleteMeshFacesUndoAction = mesh.DeleteMeshFacesInPolygon(boundaryWithMissingElements); // Compute interior boundary polygon points - std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); + // std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); + std::vector boundaryNodes2 = MeshBoundaryExtractor::ExtractConcatenated (mesh, BoundarySelection::InteriorOnly); // The expected number of points include the land boundary points UInt expectedNumberOfNodes = 26; @@ -2902,43 +2903,43 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa // interior set of polygons std::vector expectedXPoints{ 95.0, - 105.0, 100.0, + 105.0, 95.0, constants::missing::doubleValue, 85.0, - 85.0, 95.0, 95.0, 85.0, + 85.0, constants::missing::doubleValue, 80.0, - 75.0, 85.0, + 75.0, 80.0, constants::missing::doubleValue, 120.0, - 115.0, 125.0, + 115.0, 120.0, constants::missing::doubleValue, 125.0, - 125.0, 135.0, 135.0, + 125.0, 125.0}; std::vector expectedYPoints{ - 15.0, 15.0, 0.0, 15.0, + 15.0, constants::missing::doubleValue, 15.0, + 15.0, 25.0, 25.0, 15.0, - 15.0, constants::missing::doubleValue, 0.0, 15.0, @@ -2951,9 +2952,9 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa 0.0, constants::missing::doubleValue, 125.0, + 125.0, 135.0, 135.0, - 125.0, 125.0}; for (size_t i = 0; i < expectedXPoints.size(); ++i) From 1b68e924a3a17bec913fc339fe320969d75b71f9 Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 16:20:43 +0200 Subject: [PATCH 29/34] GRIDEDIT-2293 Removed commented out code --- libs/MeshKernel/tests/src/MeshRefinementTests.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index b4cfdd19f..b16edbc28 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -2891,7 +2891,6 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa auto deleteMeshFacesUndoAction = mesh.DeleteMeshFacesInPolygon(boundaryWithMissingElements); // Compute interior boundary polygon points - // std::vector boundaryNodes2 = mesh.ComputeInnerBoundaryPolygons(); std::vector boundaryNodes2 = MeshBoundaryExtractor::ExtractConcatenated (mesh, BoundarySelection::InteriorOnly); // The expected number of points include the land boundary points From 56e2d171d4a55936d3bf32cd474c445a54160b6f Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 4 Aug 2026 16:30:39 +0200 Subject: [PATCH 30/34] GRIDEDIT-2293 Fixed formatting warning --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 1 - libs/MeshKernel/src/Mesh2D.cpp | 3 +-- libs/MeshKernel/tests/src/Mesh2DTest.cpp | 1 - libs/MeshKernel/tests/src/MeshRefinementTests.cpp | 2 +- 4 files changed, 2 insertions(+), 5 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index b4e283ee4..f7855a9ec 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -447,7 +447,6 @@ namespace meshkernel std::vector& meshBoundaryPolygon, std::vector& boundaryPolygonFaceId) const; - /// @brief Reconstruct the invalid cell polygons /// /// When constructing the invalid cell polygons, they can be computed with many smaller polygons. diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 24e368482..238ba9378 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -219,7 +219,7 @@ void Mesh2D::DoAdministrationGivenFaceNodesMapping(const std::vector Mesh2D::GetHangingEdges() const { std::vector result; diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index 91a9a04fc..994eb6198 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -309,7 +309,6 @@ TEST(Mesh2D, MeshBoundaryToPolygonWithSelection) meshkernel::Print(mesh.Nodes(), mesh.Edges()); - // 3 Validation const double tolerance = 1e-5; ASSERT_EQ(9, meshBoundaryPolygon.size()); diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index b16edbc28..a92a046a4 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -2891,7 +2891,7 @@ TEST(MeshRefinement, MeshWithHole_ShouldGenerateInteriorBoundaryPolygonsForSixFa auto deleteMeshFacesUndoAction = mesh.DeleteMeshFacesInPolygon(boundaryWithMissingElements); // Compute interior boundary polygon points - std::vector boundaryNodes2 = MeshBoundaryExtractor::ExtractConcatenated (mesh, BoundarySelection::InteriorOnly); + std::vector boundaryNodes2 = MeshBoundaryExtractor::ExtractConcatenated(mesh, BoundarySelection::InteriorOnly); // The expected number of points include the land boundary points UInt expectedNumberOfNodes = 26; From e36b89cb2ff224c4b760a5afa044c546c7faa556 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Wed, 5 Aug 2026 15:59:41 +0200 Subject: [PATCH 31/34] GRIDEDIT-2293 Added finding of boundary within a bounding polygon --- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 7 - .../MeshKernel/MeshBoundaryExtractor.hpp | 19 ++- .../MeshKernel/include/MeshKernel/Polygon.hpp | 8 + libs/MeshKernel/src/Mesh2D.cpp | 147 +----------------- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 78 ++++++++-- libs/MeshKernel/tests/src/Mesh2DTest.cpp | 70 ++++++--- libs/MeshKernel/tests/src/MeshTests.cpp | 4 +- 7 files changed, 146 insertions(+), 187 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index f7855a9ec..ae7250e58 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -440,13 +440,6 @@ namespace meshkernel /// @brief Find the mesh faces that lie entirely within the polygon. std::vector FindFacesEntirelyInsidePolygon(const std::vector& isNodeInsidePolygon) const; - /// @brief Constructs a polygon from the meshboundary, by walking through the mesh - void WalkBoundaryFromNode(const Polygon& polygon, - std::vector& isVisited, - UInt& currentNode, - std::vector& meshBoundaryPolygon, - std::vector& boundaryPolygonFaceId) const; - /// @brief Reconstruct the invalid cell polygons /// /// When constructing the invalid cell polygons, they can be computed with many smaller polygons. diff --git a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp index e77261457..ca5ff6f17 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp @@ -33,6 +33,7 @@ #include "MeshKernel/Definitions.hpp" #include "MeshKernel/Mesh2D.hpp" #include "MeshKernel/Point.hpp" +#include "MeshKernel/Polygon.hpp" namespace meshkernel { @@ -44,13 +45,20 @@ namespace meshkernel /// @brief Extract all boundaries as a single sequence of points, separated by an invalid point static std::vector ExtractConcatenated(const Mesh2D& mesh, BoundarySelection boundaryType = BoundarySelection::All); - /// @brief Extract all boundaries keeping them separated and + /// @brief Extract all boundaries keeping them separated /// /// The result consists of an array of each of the boundary polygons /// Additionally, an array indicating if the boundary polygon is a exterior boundary or not. /// True => is-exterior, False => otherwise static std::tuple>, std::vector> Extract(const Mesh2D& mesh); + /// @brief Extract all boundaries contained within a constraining polygon + /// + /// The result consists of an array of each of the boundary polygons + /// Additionally, an array indicating if the boundary polygon is a exterior boundary or not. + /// True => is-exterior, False => otherwise + static std::tuple>, std::vector> Extract(const Mesh2D& mesh, const Polygon& polygon); + private: /// @brief Temporary struct, used when computing the boundaries struct BoundaryEdge @@ -84,10 +92,17 @@ namespace meshkernel std::vector& isExterior, std::vector>& separatedBoundaryPolygons); + /// @brief Clip a boundary polygon node sequence to be contains within a constraining polygon + /// + /// Edges that cross the constraining polygons are also included, i.e. edges that have 1 one node + /// contained with the polygons and another not. + static void ClipToConstrainingPolygon(const Polygon& polygon, std::vector& nodes); + /// @brief Find boundary loops /// /// Any boundary loops found may need to be processed further as they may themselves contain sub-loops - static void FindBoundaryPolygons(const std::vector& nodes, + static void FindBoundaryPolygons(const Polygon& polygon, + const std::vector& nodes, const std::vector& edges, const std::vector>& edgesFaces, std::vector>& allPolygons, diff --git a/libs/MeshKernel/include/MeshKernel/Polygon.hpp b/libs/MeshKernel/include/MeshKernel/Polygon.hpp index 8b7edb995..487d95ab2 100644 --- a/libs/MeshKernel/include/MeshKernel/Polygon.hpp +++ b/libs/MeshKernel/include/MeshKernel/Polygon.hpp @@ -78,6 +78,9 @@ namespace meshkernel /// @brief Move assignment operator Polygon& operator=(Polygon&& copy); + /// @brief Determine if the polygon is empty, i.e. has no control points. + bool IsEmpty() const; + /// @brief Return the number of points in the polygon UInt Size() const; @@ -223,6 +226,11 @@ inline meshkernel::UInt meshkernel::Polygon::Size() const return static_cast(m_nodes.size()); } +inline bool meshkernel::Polygon::IsEmpty() const +{ + return m_nodes.empty(); +} + inline const std::vector& meshkernel::Polygon::Nodes() const { return m_nodes; diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 238ba9378..4efb7b33b 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1523,12 +1523,11 @@ std::tuple, std::vector> Mesh2D::GetAllBoun { const Polygon polygon(polygonNodes, m_projection); -#if 0 MeshBoundaryExtractor meshBoundaryExtractor; - auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this); + auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this, polygon); - if (polygonNodes.size() == 0) + if (polygon.IsEmpty()) { auto alwaysTrue = [](size_t idx [[maybe_unused]]) @@ -1541,6 +1540,7 @@ std::tuple, std::vector> Mesh2D::GetAllBoun std::vector containedIsEnclosingBoundary; + // NOTE: addNonEmpty lambda also updates containedIsEnclosingBoundary auto addNonEmpty = [&boundaryPoints, &isEnclosingBoundary, &containedIsEnclosingBoundary](size_t idx) mutable { if (boundaryPoints[idx].size() > 0) @@ -1555,110 +1555,6 @@ std::tuple, std::vector> Mesh2D::GetAllBoun std::vector containedBoundaryPoints = ConcatenatePointVectors(boundaryPoints, addNonEmpty); return {containedBoundaryPoints, containedIsEnclosingBoundary}; - -#else - std::vector isVisited(GetNumEdges(), false); - std::vector meshBoundaryPolygon; - std::vector boundaryPolygon; - // Elements connected to the boundary - std::vector boundaryPolygonFaceId; - std::vector isEnclosingBoundary; - - meshBoundaryPolygon.reserve(GetNumNodes()); - boundaryPolygon.reserve(GetNumNodes()); - - for (UInt e = 0; e < GetNumEdges(); e++) - { - if (isVisited[e] || !IsEdgeOnBoundary(e)) - { - continue; - } - - const UInt firstNodeIndex = m_edges[e].first; - const UInt secondNodeIndex = m_edges[e].second; - const Point firstNode = m_nodes[firstNodeIndex]; - const Point secondNode = m_nodes[secondNodeIndex]; - - bool firstNodeInPolygon = polygon.Contains(m_nodes[firstNodeIndex]); - bool secondNodeInPolygon = polygon.Contains(m_nodes[secondNodeIndex]); - - if (!firstNodeInPolygon && !secondNodeInPolygon) - { - continue; - } - - // Start a new polyline - boundaryPolygon.clear(); - boundaryPolygonFaceId.clear(); - - // Put the current edge on the mesh boundary, mark it as visited - boundaryPolygon.emplace_back(firstNode); - boundaryPolygon.emplace_back(secondNode); - boundaryPolygonFaceId.push_back(m_edgesFaces[e][0]); - - isVisited[e] = true; - - // walk the current mesh boundary - auto currentNode = secondNodeIndex; - WalkBoundaryFromNode(polygon, isVisited, currentNode, boundaryPolygon, boundaryPolygonFaceId); - - const auto numNodesFirstTail = static_cast(boundaryPolygon.size()); - - // if the boundary polygon is not closed - if (currentNode != firstNodeIndex) - { - // Now grow a polyline starting at the other side of the original link L, i.e., the second tail - currentNode = firstNodeIndex; - WalkBoundaryFromNode(polygon, isVisited, currentNode, boundaryPolygon, boundaryPolygonFaceId); - } - - // There is a nonempty second tail: reverse the second tail so that the tails connect and close the polygon. - if (boundaryPolygon.size() > numNodesFirstTail) - { - std::reverse(boundaryPolygon.begin() + numNodesFirstTail, boundaryPolygon.end()); - boundaryPolygon.push_back(boundaryPolygon.front()); - } - - std::span currentPolygonSpan(boundaryPolygon); - Polygon currentPolygon(currentPolygonSpan, m_projection); - - if (!meshBoundaryPolygon.empty()) - { - meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); - } - - if (IsMultiPolygon(currentPolygonSpan)) - { - auto [multiBoundaryPolygonNodes, multiBoundaryElementIds] = (SplitMultiplePolygons(currentPolygonSpan, boundaryPolygonFaceId)); - - for (size_t i = 0; i < multiBoundaryPolygonNodes.size(); ++i) - { - if (i > 0) - { - meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue); - } - - meshBoundaryPolygon.insert(meshBoundaryPolygon.end(), multiBoundaryPolygonNodes[i].begin(), multiBoundaryPolygonNodes[i].end()); - - std::span currentPolygonSpan(multiBoundaryPolygonNodes[i]); - Polygon currentPolygon(currentPolygonSpan, m_projection); - - isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[multiBoundaryElementIds[i]])); - } - } - else - { - meshBoundaryPolygon.insert(meshBoundaryPolygon.end(), boundaryPolygon.begin(), boundaryPolygon.end()); - - std::span currentPolygonSpan(boundaryPolygon); - Polygon currentPolygon(currentPolygonSpan, m_projection); - - isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[boundaryPolygonFaceId[0]])); - } - } - - return {meshBoundaryPolygon, isEnclosingBoundary}; -#endif } std::vector Mesh2D::RemoveOuterDomainBoundaryPolygon(const std::vector& polygonNodes) const @@ -1718,43 +1614,6 @@ std::vector Mesh2D::RemoveOuterDomainBoundaryPolygon(const st return innerBoundaryNodes; } -void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon, - std::vector& isVisited, - UInt& currentNode, - std::vector& meshBoundaryPolygon, - std::vector& boundaryPolygonFaceId) const -{ - UInt e = 0; - bool currentNodeInPolygon = false; - while (e < m_nodesNumEdges[currentNode]) - { - if (!currentNodeInPolygon) - { - currentNodeInPolygon = polygon.Contains(m_nodes[currentNode]); - } - - if (!currentNodeInPolygon) - { - break; - } - - const auto currentEdge = m_nodesEdges[currentNode][e]; - if (isVisited[currentEdge] || !IsEdgeOnBoundary(currentEdge)) - { - e++; - continue; - } - - currentNode = OtherNodeOfEdge(m_edges[currentEdge], currentNode); - boundaryPolygonFaceId.push_back(m_edgesFaces[currentEdge][0]); - e = 0; - currentNodeInPolygon = false; - - meshBoundaryPolygon.emplace_back(m_nodes[currentNode]); - isVisited[currentEdge] = true; - } -} - std::vector Mesh2D::GetHangingEdges() const { std::vector result; diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index 976ffcda9..cb0e4998d 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -59,6 +59,24 @@ std::vector meshkernel::MeshBoundaryExtractor::ExtractConcate return ConcatenatePointVectors(boundarySequences, isBoundarySelection); } +std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh) +{ + const Polygon emptyPolygon; + return Extract(mesh, emptyPolygon); +} + +std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh, const Polygon& polygon) +{ + + std::vector> allBoundaryPolygons; + std::vector> allTouchedFaces; + + const std::vector& meshNodes(mesh.Nodes()); + + FindBoundaryPolygons(polygon, meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryPolygons, allTouchedFaces); + return SeparateAndDetermineExternality(mesh, allBoundaryPolygons, allTouchedFaces); +} + void meshkernel::MeshBoundaryExtractor::Append(const Point& centre, const Projection projection, std::vector& boundaryPolygon, @@ -109,18 +127,6 @@ meshkernel::MeshBoundaryExtractor::SeparateAndDetermineExternality(const Mesh2D& return {separatedBoundaryPolygons, isExterior}; } -std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::Extract(const Mesh2D& mesh) -{ - - std::vector> allBoundaryPolygons; - std::vector> allTouchedFaces; - - const std::vector& meshNodes(mesh.Nodes()); - - FindBoundaryPolygons(meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryPolygons, allTouchedFaces); - return SeparateAndDetermineExternality(mesh, allBoundaryPolygons, allTouchedFaces); -} - double meshkernel::MeshBoundaryExtractor::NormalizeAngle(double angle) { while (angle < 0) @@ -146,7 +152,7 @@ void meshkernel::MeshBoundaryExtractor::FindAllBoundarEdges(const std::vector& nodes, +void meshkernel::MeshBoundaryExtractor::ClipToConstrainingPolygon(const Polygon& polygon, std::vector& nodes) +{ + if (polygon.IsEmpty() || nodes.size() <= 2) + { + return; + } + + std::vector clippedPolygon; + std::vector nodeIsContained(nodes.size()); + UInt numberOfNodesContained = 0; + + for (size_t i = 0; i < nodes.size(); ++i) + { + nodeIsContained[i] = polygon.Contains(nodes[i]); + numberOfNodesContained += nodeIsContained[i] ? 1 : 0; + } + + if (numberOfNodesContained <= 1) + { + // Cannot make a reasonable polygon with only 1 node (and the closing node htat will be added later, making 2 nodes) + nodes.clear(); + return; + } + + clippedPolygon.reserve(nodes.size()); + + // At this point the boundary polygon is open, i.e. the last point in the sequence is not the same as the first. + for (size_t i = 0; i < nodes.size(); ++i) + { + UInt nextNodeId = (i + 1) % nodes.size(); + UInt previousNodeId = (i + nodes.size() - 1) % nodes.size(); + + if (nodeIsContained[i] || nodeIsContained[nextNodeId] || nodeIsContained[previousNodeId]) + { + clippedPolygon.push_back(nodes[i]); + } + } + + nodes = std::move(clippedPolygon); +} + +void meshkernel::MeshBoundaryExtractor::FindBoundaryPolygons(const Polygon& polygon, + const std::vector& nodes, const std::vector& edges, const std::vector>& edgesFaces, std::vector>& allPolygons, @@ -273,6 +321,8 @@ void meshkernel::MeshBoundaryExtractor::FindBoundaryPolygons(const std::vector

= 3) { // Close the polygon diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index 994eb6198..69cc53a30 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -242,9 +242,9 @@ TEST(Mesh2D, MeshBoundaryToPolygon) ASSERT_NEAR(0.0, meshBoundaryPolygon[4].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].y, tolerance); - ASSERT_NEAR(5.0, meshBoundaryPolygon[1].y, tolerance); + ASSERT_NEAR(-5.0, meshBoundaryPolygon[1].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[2].y, tolerance); - ASSERT_NEAR(-5.0, meshBoundaryPolygon[3].y, tolerance); + ASSERT_NEAR(5.0, meshBoundaryPolygon[3].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[4].y, tolerance); } @@ -298,6 +298,7 @@ TEST(Mesh2D, MeshBoundaryToPolygonWithSelection) auto mesh = meshkernel::Mesh2D(edges, nodes, meshkernel::Projection::cartesian); std::vector polygonNodes; + polygonNodes.push_back({-5.0, 35.0}); polygonNodes.push_back({15.0, 35.0}); polygonNodes.push_back({15.0, -5.0}); @@ -305,9 +306,7 @@ TEST(Mesh2D, MeshBoundaryToPolygonWithSelection) polygonNodes.push_back({-5.0, 35.0}); // 2 Execution - const auto meshBoundaryPolygon = mesh.ComputeBoundaryPolygons(polygonNodes); - - meshkernel::Print(mesh.Nodes(), mesh.Edges()); + auto meshBoundaryPolygon = mesh.ComputeBoundaryPolygons(polygonNodes); // 3 Validation const double tolerance = 1e-5; @@ -315,31 +314,66 @@ TEST(Mesh2D, MeshBoundaryToPolygonWithSelection) ASSERT_NEAR(0.0, meshBoundaryPolygon[0].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].y, tolerance); - ASSERT_NEAR(10.0, meshBoundaryPolygon[7].x, tolerance); - ASSERT_NEAR(00.0, meshBoundaryPolygon[7].y, tolerance); + ASSERT_NEAR(10.0, meshBoundaryPolygon[1].x, tolerance); + ASSERT_NEAR(0.0, meshBoundaryPolygon[1].y, tolerance); - ASSERT_NEAR(10.0, meshBoundaryPolygon[4].x, tolerance); - ASSERT_NEAR(30.0, meshBoundaryPolygon[4].y, tolerance); + ASSERT_NEAR(20.0, meshBoundaryPolygon[2].x, tolerance); + ASSERT_NEAR(0.0, meshBoundaryPolygon[2].y, tolerance); - ASSERT_NEAR(0.0, meshBoundaryPolygon[3].x, tolerance); + ASSERT_NEAR(20.0, meshBoundaryPolygon[3].x, tolerance); ASSERT_NEAR(30.0, meshBoundaryPolygon[3].y, tolerance); - ASSERT_NEAR(0.0, meshBoundaryPolygon[2].x, tolerance); - ASSERT_NEAR(20.0, meshBoundaryPolygon[2].y, tolerance); - - ASSERT_NEAR(0.0, meshBoundaryPolygon[1].x, tolerance); - ASSERT_NEAR(10.0, meshBoundaryPolygon[1].y, tolerance); + ASSERT_NEAR(10.0, meshBoundaryPolygon[4].x, tolerance); + ASSERT_NEAR(30.0, meshBoundaryPolygon[4].y, tolerance); - ASSERT_NEAR(20.0, meshBoundaryPolygon[5].x, tolerance); + ASSERT_NEAR(0.0, meshBoundaryPolygon[5].x, tolerance); ASSERT_NEAR(30.0, meshBoundaryPolygon[5].y, tolerance); - ASSERT_NEAR(20.0, meshBoundaryPolygon[6].x, tolerance); - ASSERT_NEAR(00.0, meshBoundaryPolygon[6].y, tolerance); + ASSERT_NEAR(0.0, meshBoundaryPolygon[6].x, tolerance); + ASSERT_NEAR(20.0, meshBoundaryPolygon[6].y, tolerance); + + ASSERT_NEAR(0.0, meshBoundaryPolygon[7].x, tolerance); + ASSERT_NEAR(10.0, meshBoundaryPolygon[7].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[8].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[8].y, tolerance); } +TEST(Mesh2D, MeshBoundaryToPolygonWithAnotherSelection) +{ + // 1 Setup + auto mesh = MakeRectangularMeshForTesting(7, 7, 5.0, meshkernel::Projection::cartesian); + + std::vector polygonNodes; + + polygonNodes.push_back({-5.0, 35.0}); + polygonNodes.push_back({15.0, 35.0}); + polygonNodes.push_back({15.0, -5.0}); + polygonNodes.push_back({-5.0, -5.0}); + // Add small kink in the polygon, it should miss out a small section of the boundary + polygonNodes.push_back({-1.0, 7.5}); + polygonNodes.push_back({10.0, 12.0}); + polygonNodes.push_back({-1.0, 22.0}); + polygonNodes.push_back({-5.0, 35.0}); + + // 2 Execution + auto meshBoundaryPolygon = mesh->ComputeBoundaryPolygons(polygonNodes); + + // Note missing point at (0.0, 15.0) + std::vector expectedXs{0.0, 5.0, 10.0, 15.0, 20.0, 20.0, 15.0, 10.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + std::vector expectedYs{0.0, 0.0, 0.0, 0.0, 0.0, 30.0, 30.0, 30.0, 30.0, 30.0, 25.0, 20.0, 10.0, 5.0, 0.0}; + + // 3 Validation + const double tolerance = 1e-5; + ASSERT_EQ(15, meshBoundaryPolygon.size()); + + for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) + { + EXPECT_NEAR(expectedXs[i], meshBoundaryPolygon[i].x, tolerance); + EXPECT_NEAR(expectedYs[i], meshBoundaryPolygon[i].y, tolerance); + } +} + TEST(Mesh2D, HangingEdge) { // 1 Setup diff --git a/libs/MeshKernel/tests/src/MeshTests.cpp b/libs/MeshKernel/tests/src/MeshTests.cpp index d89e2e041..aef9995a0 100644 --- a/libs/MeshKernel/tests/src/MeshTests.cpp +++ b/libs/MeshKernel/tests/src/MeshTests.cpp @@ -271,9 +271,9 @@ TEST(Mesh, MeshBoundaryToPolygon) ASSERT_NEAR(0.0, meshBoundaryPolygon[4].x, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[0].y, tolerance); - ASSERT_NEAR(5.0, meshBoundaryPolygon[1].y, tolerance); + ASSERT_NEAR(-5.0, meshBoundaryPolygon[1].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[2].y, tolerance); - ASSERT_NEAR(-5.0, meshBoundaryPolygon[3].y, tolerance); + ASSERT_NEAR(5.0, meshBoundaryPolygon[3].y, tolerance); ASSERT_NEAR(0.0, meshBoundaryPolygon[4].y, tolerance); } From 829d087e39d7e262ba49e0a593cce34e5c02d579 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Wed, 5 Aug 2026 16:58:15 +0200 Subject: [PATCH 32/34] GRIDEDIT-2293 Fixed windows and macos builds --- libs/MeshKernel/src/Mesh2D.cpp | 10 +++++++++- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 4 ++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 4efb7b33b..1621e73f7 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1525,7 +1525,14 @@ std::tuple, std::vector> Mesh2D::GetAllBoun MeshBoundaryExtractor meshBoundaryExtractor; - auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this, polygon); + // The use of std::tie instead of a structured binding is due to limitations in the macos compiler + // The compiler error "error: capturing a structured binding is not yet supported in OpenMP" + // + // auto [boundaryPoints, isEnclosingBoundary] = meshBoundaryExtractor.Extract(*this, polygon); + // Replace the 3 lines below with the line above + std::vector> boundaryPoints; + std::vector isEnclosingBoundary; + std::tie(boundaryPoints, isEnclosingBoundary) = meshBoundaryExtractor.Extract(*this, polygon); if (polygon.IsEmpty()) { @@ -1920,6 +1927,7 @@ void Mesh2D::ReconstructInvalidCellsPolygon() auto isInteriorBoundary = [&isEnclosingBoundary](size_t i) { return !isEnclosingBoundary[i]; }; + m_invalidCellPolygons = ConcatenatePointVectors(boundaryPoints, isInteriorBoundary); } diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index cb0e4998d..7f186d959 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -238,8 +238,8 @@ void meshkernel::MeshBoundaryExtractor::ClipToConstrainingPolygon(const Polygon& // At this point the boundary polygon is open, i.e. the last point in the sequence is not the same as the first. for (size_t i = 0; i < nodes.size(); ++i) { - UInt nextNodeId = (i + 1) % nodes.size(); - UInt previousNodeId = (i + nodes.size() - 1) % nodes.size(); + size_t nextNodeId = (i + 1) % nodes.size(); + size_t previousNodeId = (i + nodes.size() - 1) % nodes.size(); if (nodeIsContained[i] || nodeIsContained[nextNodeId] || nodeIsContained[previousNodeId]) { From ba642b297f5642da863add05fd75b50407e99f33 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 6 Aug 2026 11:07:58 +0200 Subject: [PATCH 33/34] GRIDEDIT-2293 Clipping of boundary polygons is now done after separation of multiple sub-polygons --- .../MeshKernel/MeshBoundaryExtractor.hpp | 13 +- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 125 ++++++++++-------- libs/MeshKernel/tests/src/Mesh2DTest.cpp | 29 ++++ 3 files changed, 109 insertions(+), 58 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp index ca5ff6f17..41fc0c4a4 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp @@ -60,6 +60,9 @@ namespace meshkernel static std::tuple>, std::vector> Extract(const Mesh2D& mesh, const Polygon& polygon); private: + /// @brief The minimum number of points in a polygon, excluding the closing point + static constexpr size_t MinimumNumberOfPoints = 3; + /// @brief Temporary struct, used when computing the boundaries struct BoundaryEdge { @@ -85,8 +88,10 @@ namespace meshkernel /// @brief Append the boundary polygon to the set of all boundary polygon /// - /// The boundary polygons may be reversed if found to be in clockwise direction - static void Append(const Point& centre, + /// The boundary polygons may be reversed if found to be in clockwise direction and clipped to be + /// inside a constraining polygon + static void Append(const Polygon& polygon, + const Point& centre, const Projection projection, std::vector& boundaryPolygon, std::vector& isExterior, @@ -101,8 +106,7 @@ namespace meshkernel /// @brief Find boundary loops /// /// Any boundary loops found may need to be processed further as they may themselves contain sub-loops - static void FindBoundaryPolygons(const Polygon& polygon, - const std::vector& nodes, + static void FindBoundaryPolygons(const std::vector& nodes, const std::vector& edges, const std::vector>& edgesFaces, std::vector>& allPolygons, @@ -114,6 +118,7 @@ namespace meshkernel /// @note allBoundaryPolygons should not be accessed after calling this function static std::tuple>, std::vector> SeparateAndDetermineExternality(const Mesh2D& mesh, + const Polygon& polygon, std::vector>& allBoundaryPolygons, const std::vector>& allTouchedFaces); }; diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index 7f186d959..4f6d95c0e 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -73,11 +73,57 @@ std::tuple>, std::vector> meshk const std::vector& meshNodes(mesh.Nodes()); - FindBoundaryPolygons(polygon, meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryPolygons, allTouchedFaces); - return SeparateAndDetermineExternality(mesh, allBoundaryPolygons, allTouchedFaces); + FindBoundaryPolygons(meshNodes, mesh.Edges(), mesh.m_edgesFaces, allBoundaryPolygons, allTouchedFaces); + return SeparateAndDetermineExternality(mesh, polygon, allBoundaryPolygons, allTouchedFaces); } -void meshkernel::MeshBoundaryExtractor::Append(const Point& centre, +void meshkernel::MeshBoundaryExtractor::ClipToConstrainingPolygon(const Polygon& polygon, std::vector& nodes) +{ + if (polygon.IsEmpty() || nodes.size() <= MinimumNumberOfPoints) + { + return; + } + + const size_t uniquePointCount = nodes.size() - 1; + std::vector clippedPolygon; + std::vector nodeIsContained(uniquePointCount); + UInt numberOfNodesContained = 0; + + for (size_t i = 0; i < uniquePointCount; ++i) + { + nodeIsContained[i] = polygon.Contains(nodes[i]); + numberOfNodesContained += nodeIsContained[i] ? 1 : 0; + } + + if (numberOfNodesContained < MinimumNumberOfPoints) + { + // Cannot make a reasonable polygon with only 2 nodes (and the closing node htat will be added later, making 3 nodes) + nodes.clear(); + return; + } + + clippedPolygon.reserve(nodes.size()); + + // At this point the boundary polygon is closed, so we have to ignore the last point on the sequence. + for (size_t i = 0; i < uniquePointCount; ++i) + { + size_t nextNodeId = (i + 1) % uniquePointCount; + size_t previousNodeId = (i + uniquePointCount - 1) % uniquePointCount; + + if (nodeIsContained[i] || nodeIsContained[nextNodeId] || nodeIsContained[previousNodeId]) + { + clippedPolygon.push_back(nodes[i]); + } + } + + // Now re-close the polygon + clippedPolygon.push_back(clippedPolygon[0]); + + nodes = std::move(clippedPolygon); +} + +void meshkernel::MeshBoundaryExtractor::Append(const Polygon& polygon, + const Point& centre, const Projection projection, std::vector& boundaryPolygon, std::vector& isExterior, @@ -90,11 +136,22 @@ void meshkernel::MeshBoundaryExtractor::Append(const Point& centre, } isExterior.push_back(IsPointInPolygonNodes(centre, boundaryPolygon, projection)); - separatedBoundaryPolygons.push_back(std::move(boundaryPolygon)); + + // Clipping to boundary polygon must be done after determining if the boundary-points form an exterior or interior boundary. + // Because the element centre (taken from the first element the boundary-polygon touches) may no longer be inside the clipped + // boundary-polygon + ClipToConstrainingPolygon(polygon, boundaryPolygon); + + if (boundaryPolygon.size() > MinimumNumberOfPoints) + { + // Only add the polygon if it has a sufficient number of points + separatedBoundaryPolygons.push_back(std::move(boundaryPolygon)); + } } std::tuple>, std::vector> meshkernel::MeshBoundaryExtractor::SeparateAndDetermineExternality(const Mesh2D& mesh, + const Polygon& polygon, std::vector>& allBoundaryPolygons, const std::vector>& allTouchedFaces) { @@ -113,14 +170,16 @@ meshkernel::MeshBoundaryExtractor::SeparateAndDetermineExternality(const Mesh2D& for (size_t j = 0; j < individualBoundaryPolygons.size(); ++j) { - Point centre = mesh.m_facesMassCenters[firstElement[j]]; - Append(centre, mesh.m_projection, individualBoundaryPolygons[j], isExterior, separatedBoundaryPolygons); + const Point centre = mesh.m_facesMassCenters[firstElement[j]]; + + Append(polygon, centre, mesh.m_projection, individualBoundaryPolygons[j], isExterior, separatedBoundaryPolygons); } } else { - Point centre = mesh.m_facesMassCenters[allTouchedFaces[i][0]]; - Append(centre, mesh.m_projection, allBoundaryPolygons[i], isExterior, separatedBoundaryPolygons); + const Point centre = mesh.m_facesMassCenters[allTouchedFaces[i][0]]; + + Append(polygon, centre, mesh.m_projection, allBoundaryPolygons[i], isExterior, separatedBoundaryPolygons); } } @@ -209,49 +268,7 @@ meshkernel::UInt meshkernel::MeshBoundaryExtractor::FindEdgeWithMinumumAngle(con return edgeIndex; } -void meshkernel::MeshBoundaryExtractor::ClipToConstrainingPolygon(const Polygon& polygon, std::vector& nodes) -{ - if (polygon.IsEmpty() || nodes.size() <= 2) - { - return; - } - - std::vector clippedPolygon; - std::vector nodeIsContained(nodes.size()); - UInt numberOfNodesContained = 0; - - for (size_t i = 0; i < nodes.size(); ++i) - { - nodeIsContained[i] = polygon.Contains(nodes[i]); - numberOfNodesContained += nodeIsContained[i] ? 1 : 0; - } - - if (numberOfNodesContained <= 1) - { - // Cannot make a reasonable polygon with only 1 node (and the closing node htat will be added later, making 2 nodes) - nodes.clear(); - return; - } - - clippedPolygon.reserve(nodes.size()); - - // At this point the boundary polygon is open, i.e. the last point in the sequence is not the same as the first. - for (size_t i = 0; i < nodes.size(); ++i) - { - size_t nextNodeId = (i + 1) % nodes.size(); - size_t previousNodeId = (i + nodes.size() - 1) % nodes.size(); - - if (nodeIsContained[i] || nodeIsContained[nextNodeId] || nodeIsContained[previousNodeId]) - { - clippedPolygon.push_back(nodes[i]); - } - } - - nodes = std::move(clippedPolygon); -} - -void meshkernel::MeshBoundaryExtractor::FindBoundaryPolygons(const Polygon& polygon, - const std::vector& nodes, +void meshkernel::MeshBoundaryExtractor::FindBoundaryPolygons(const std::vector& nodes, const std::vector& edges, const std::vector>& edgesFaces, std::vector>& allPolygons, @@ -321,11 +338,11 @@ void meshkernel::MeshBoundaryExtractor::FindBoundaryPolygons(const Polygon& poly incomingAngle = chosenEdge.angle; } - ClipToConstrainingPolygon(polygon, currentNodes); - - if (currentNodes.size() >= 3) + if (currentNodes.size() >= MinimumNumberOfPoints) { // Close the polygon + // The polygon needs to be closed here. If it is composed of multiple sub-polygons then closing here + // makes separating them later a much easier task currentNodes.push_back(currentNodes.front()); allPolygons.emplace_back(std::move(currentNodes)); allTouchedFaces.push_back(std::move(currentFaces)); diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index 69cc53a30..05b6cc4e6 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -374,6 +374,35 @@ TEST(Mesh2D, MeshBoundaryToPolygonWithAnotherSelection) } } +TEST(Mesh2D, MeshBoundaryToPolygonWithSelectionExcludingFirstFoundPoint) +{ + // 1 Setup + auto mesh = MakeRectangularMeshForTesting(7, 7, 5.0, meshkernel::Projection::cartesian); + + // The way that the mesh is constructed means that when the full boundary polygon is computed + // the origin (0.0) and connected edges will be part of this polygon. This point is also used + // to close the polygon, since it is the first in the seuqnece. + // However this clipping polygons defined below, this point and its neighbouring points will outside of this clipping polygon + // The test is to check that this is handled correctly and that the boudnary polygon is then closed correctly. + std::vector polygonNodes{{7.5, -1.0}, {18.0, -1.0}, {18.0, 31.0}, {-1.0, 31.0}, {-1.0, 14.0}, {7.5, -1.0}}; + + // 2 Execution + auto meshBoundaryPolygon = mesh->ComputeBoundaryPolygons(polygonNodes); + + std::vector expectedXs{5.0, 10.0, 15.0, 20.0, 20.0, 15.0, 10.0, 5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0}; + std::vector expectedYs{0.0, 0.0, 0.0, 0.0, 30.0, 30.0, 30.0, 30.0, 30.0, 25.0, 20.0, 15.0, 10.0, 0.0}; + + // 3 Validation + const double tolerance = 1e-5; + ASSERT_EQ(14, meshBoundaryPolygon.size()); + + for (size_t i = 0; i < meshBoundaryPolygon.size(); ++i) + { + EXPECT_NEAR(expectedXs[i], meshBoundaryPolygon[i].x, tolerance); + EXPECT_NEAR(expectedYs[i], meshBoundaryPolygon[i].y, tolerance); + } +} + TEST(Mesh2D, HangingEdge) { // 1 Setup From d175dee23b68954d49f4505227c6e636212ed488 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 6 Aug 2026 11:40:22 +0200 Subject: [PATCH 34/34] GRIDEDIT-2293 Add the exterior/interior indicator only if the polygon is added --- libs/MeshKernel/src/MeshBoundaryExtractor.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp index 4f6d95c0e..ccd28fbe8 100644 --- a/libs/MeshKernel/src/MeshBoundaryExtractor.cpp +++ b/libs/MeshKernel/src/MeshBoundaryExtractor.cpp @@ -135,7 +135,7 @@ void meshkernel::MeshBoundaryExtractor::Append(const Polygon& polygon, std::ranges::reverse(boundaryPolygon); } - isExterior.push_back(IsPointInPolygonNodes(centre, boundaryPolygon, projection)); + const bool isExteriorBoundary = IsPointInPolygonNodes(centre, boundaryPolygon, projection); // Clipping to boundary polygon must be done after determining if the boundary-points form an exterior or interior boundary. // Because the element centre (taken from the first element the boundary-polygon touches) may no longer be inside the clipped @@ -144,7 +144,8 @@ void meshkernel::MeshBoundaryExtractor::Append(const Polygon& polygon, if (boundaryPolygon.size() > MinimumNumberOfPoints) { - // Only add the polygon if it has a sufficient number of points + // Only add the polygon and the exterior/interior indicator if the polygon has a sufficient number of points + isExterior.push_back(isExteriorBoundary); separatedBoundaryPolygons.push_back(std::move(boundaryPolygon)); } }