From 2a44019ca525bca03dff8861096086fd7d1a55dd Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 2 Jul 2026 17:49:51 +0200 Subject: [PATCH 01/12] 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/12] 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 a1e171ff260c9df3e69eb09cb1e08758b9b0b2d9 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 2 Jul 2026 17:59:26 +0200 Subject: [PATCH 03/12] 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/12] 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 6bdf18625779f40262a0cb0e9eceafb5dae082ef Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 15:01:12 +0200 Subject: [PATCH 05/12] 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 3f22acc8777fca912f8f2ec916defc3f0e3abed0 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 15:08:26 +0200 Subject: [PATCH 06/12] 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 d50466dc880f7996deaabfa574159cfea423ce66 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 15:09:13 +0200 Subject: [PATCH 07/12] 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 58c7e9d12608f473a292305f520c548739169151 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 15:30:04 +0200 Subject: [PATCH 08/12] 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 4e7b8f6fcbff2a696e502fb395aa24fcc9ceb770 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 6 Jul 2026 18:21:45 +0200 Subject: [PATCH 09/12] 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 869865f2b0d79b158f649bf9a356bb3f335a7ad3 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Thu, 9 Jul 2026 18:28:58 +0200 Subject: [PATCH 10/12] 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 393c34965dc9942a11e82003b57fe59f7ada5e4a Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 13 Jul 2026 11:15:59 +0200 Subject: [PATCH 11/12] 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 3d41b0bef1f859d8f472cbd117719fff6484e2d3 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 13 Jul 2026 12:56:20 +0200 Subject: [PATCH 12/12] GRIDEDIT-2102 Fix compilation under macos --- cmake/compiler_config.cmake | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/cmake/compiler_config.cmake b/cmake/compiler_config.cmake index ad5a12ff5..0e8a1a8e0 100644 --- a/cmake/compiler_config.cmake +++ b/cmake/compiler_config.cmake @@ -15,6 +15,13 @@ if(APPLE) message(STATUS "Configuring build for macOS with ${CMAKE_CXX_COMPILER_ID} (${CMAKE_CXX_COMPILER_VERSION}).") # Common warning and visibility flags add_compile_options("-fvisibility=hidden;-Wall;-Wextra;-pedantic;-Werror;-Wno-unused-function") + + cmake_host_system_information(RESULT os_release QUERY OS_RELEASE) + + if(os_release VERSION_GREATER_EQUAL "26.0") + add_compile_options("-Wno-character-conversion") + endif() + # Conditionally suppress unused parameter warnings (used for platform-specific compiler issues) if(SUPPRESS_UNUSED_PARAMETER_WARNING) add_compile_options("-Wno-unused-parameter")