Feature/gridedit 2293 separate polygon loops

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")

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -448,7 +448,8 @@ namespace meshkernel
         void WalkBoundaryFromNode(const Polygon& polygon,
                                   std::vector<bool>& isVisited,
                                   UInt& currentNode,
-                                  std::vector<Point>& meshBoundaryPolygon) const;
+                                  std::vector<Point>& meshBoundaryPolygon,
+                                  std::vector<UInt>& boundaryPolygonFaceId) const;
 
         /// @brief Constructs a polygon or polygons from the meshboundary, by walking through the mesh
         ///
@@ -461,6 +462,18 @@ namespace meshkernel
                                        std::vector<Point>& subSequence,
                                        std::vector<Point>& illegalCells) const;
 
+        /// @brief Reconstruct the invalid cell 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
+        /// polygonal sub-sequence forms a external boundary
+        [[nodiscard]] std::tuple<std::vector<meshkernel::Point>, std::vector<bool>> GetAllBoundaryPolygons(const std::vector<Point>& polygon);
+
         /// @brief Ensure that all polynomials are orientated in the ACW direction.
         void OrientatePolygonsAntiClockwise(std::vector<Point>& polygonNodes) const;
 

libs/MeshKernel/include/MeshKernel/Operations.hpp

@@ -621,4 +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<const Point> 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<std::vector<Point>>, std::vector<UInt>> splitMultiplePolygons(std::span<const Point> boundary, std::span<UInt> elementIds);
+
 } // namespace meshkernel

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.

libs/MeshKernel/include/MeshKernel/Polygon.hpp

@@ -27,6 +27,7 @@
 
 #pragma once
 
+#include <span>
 #include <vector>
 
 #include "MeshKernel/BoundingBox.hpp"
@@ -54,6 +55,10 @@ namespace meshkernel
         /// @brief Default move constructor.
         Polygon(Polygon&& copy) = default;
 
+        /// @brief Constructor
+        Polygon(std::span<const Point> points,
+                Projection projection);
+
         /// @brief Constructor
         Polygon(const std::vector<Point>& points,
                 Projection projection);

libs/MeshKernel/src/Mesh2D.cpp

@@ -1511,14 +1511,25 @@ std::vector<meshkernel::UInt> Mesh2D::SortedFacesAroundNode(UInt node) const
 }
 
 std::vector<meshkernel::Point> Mesh2D::ComputeBoundaryPolygons(const std::vector<Point>& polygonNodes)
+{
+    Administrate();
+    auto [boundaryPoints, isEnclosingBoundary] = GetAllBoundaryPolygons(polygonNodes);
+    return boundaryPoints;
+}
+
+std::tuple<std::vector<meshkernel::Point>, std::vector<bool>> Mesh2D::GetAllBoundaryPolygons(const std::vector<Point>& polygonNodes)
 {
     const Polygon polygon(polygonNodes, m_projection);
 
-    // Find faces
-    Administrate();
     std::vector<bool> isVisited(GetNumEdges(), false);
     std::vector<Point> meshBoundaryPolygon;
+    std::vector<Point> boundaryPolygon;
+    // Elements connected to the boundary
+    std::vector<UInt> boundaryPolygonFaceId;
+    std::vector<bool> isEnclosingBoundary;
+
     meshBoundaryPolygon.reserve(GetNumNodes());
+    boundaryPolygon.reserve(GetNumNodes());
 
     for (UInt e = 0; e < GetNumEdges(); e++)
     {
@@ -1527,10 +1538,10 @@ std::vector<meshkernel::Point> 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]);
@@ -1541,38 +1552,77 @@ std::vector<meshkernel::Point> Mesh2D::ComputeBoundaryPolygons(const std::vector
         }
 
         // Start a new polyline
-        if (!meshBoundaryPolygon.empty())
-        {
-            meshBoundaryPolygon.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue);
-        }
+        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<UInt>(meshBoundaryPolygon.size());
+        const auto numNodesFirstTail = static_cast<UInt>(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(boundaryPolygon.begin() + numNodesFirstTail, boundaryPolygon.end());
+            boundaryPolygon.push_back(boundaryPolygon.front());
+        }
+
+        std::span<const Point> 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<const Point> currentPolygonSpan(multiBoundaryPolygonNodes[i]);
+                Polygon currentPolygon(currentPolygonSpan, m_projection);
+
+                isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[multiBoundaryElementIds[i]]));
+            }
+        }
+        else
         {
-            std::reverse(meshBoundaryPolygon.begin() + numNodesFirstTail, meshBoundaryPolygon.end());
-            meshBoundaryPolygon.push_back(meshBoundaryPolygon.front());
+            meshBoundaryPolygon.insert(meshBoundaryPolygon.end(), boundaryPolygon.begin(), boundaryPolygon.end());
+
+            std::span<const Point> currentPolygonSpan(boundaryPolygon);
+            Polygon currentPolygon(currentPolygonSpan, m_projection);
+
+            isEnclosingBoundary.push_back(currentPolygon.Contains(m_facesMassCenters[boundaryPolygonFaceId[0]]));
         }
     }
-    return meshBoundaryPolygon;
+
+    return {meshBoundaryPolygon, isEnclosingBoundary};
 }
 
 std::vector<meshkernel::Point> Mesh2D::ComputeInnerBoundaryPolygons() const
@@ -1778,7 +1828,8 @@ std::vector<meshkernel::Point> Mesh2D::RemoveOuterDomainBoundaryPolygon(const st
 void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon,
                                   std::vector<bool>& isVisited,
                                   UInt& currentNode,
-                                  std::vector<Point>& meshBoundaryPolygon) const
+                                  std::vector<Point>& meshBoundaryPolygon,
+                                  std::vector<UInt>& boundaryPolygonFaceId) const
 {
     UInt e = 0;
     bool currentNodeInPolygon = false;
@@ -1802,6 +1853,7 @@ void Mesh2D::WalkBoundaryFromNode(const Polygon& polygon,
         }
 
         currentNode = OtherNodeOfEdge(m_edges[currentEdge], currentNode);
+        boundaryPolygonFaceId.push_back(m_edgesFaces[currentEdge][0]);
         e = 0;
         currentNodeInPolygon = false;
 
@@ -2183,6 +2235,45 @@ std::unique_ptr<meshkernel::UndoAction> Mesh2D::DeleteMeshFacesInPolygon(const P
     return UpdateFaceInformation(faceIndices, appendDeletedFaces);
 }
 
+void Mesh2D::ReconstructInvalidCellsPolygon()
+{
+    auto [boundaryPoints, isEnclosingBoundary] = GetAllBoundaryPolygons(std::vector<meshkernel::Point>());
+
+    Polygons polygons(boundaryPoints, m_projection);
+
+    if (polygons.GetNumPolygons() <= 1)
+    {
+        // There are no interior boundary polygons
+        return;
+    }
+
+    std::vector<Point> innerBoundaryPolygons;
+    innerBoundaryPolygons.reserve(m_invalidCellPolygons.size());
+    bool firstElement = true;
+
+    for (UInt p = 0; p < polygons.GetNumPolygons(); ++p)
+    {
+        if (isEnclosingBoundary[p])
+        {
+            continue;
+        }
+
+        if (!firstElement)
+        {
+            innerBoundaryPolygons.push_back(Point(constants::missing::doubleValue, constants::missing::doubleValue));
+        }
+        else
+        {
+            firstElement = false;
+        }
+
+        const std::vector<Point>& polygonPoints(polygons.Enclosure(p).Outer().Nodes());
+        innerBoundaryPolygons.insert(innerBoundaryPolygons.end(), polygonPoints.begin(), polygonPoints.end());
+    }
+
+    m_invalidCellPolygons = std::move(innerBoundaryPolygons);
+}
+
 std::unique_ptr<meshkernel::UndoAction> Mesh2D::UpdateFaceInformation(const std::vector<UInt>& faceIndices, const bool appendDeletedFaces)
 {
     std::vector<UInt> facesToDelete;
@@ -2273,6 +2364,8 @@ std::unique_ptr<meshkernel::UndoAction> Mesh2D::UpdateFaceInformation(const std:
         m_faceArea.erase(m_faceArea.begin() + faceId);
     }
 
+    ReconstructInvalidCellsPolygon();
+
     return deleteMeshAction;
 }
 

libs/MeshKernel/src/Operations.cpp

@@ -25,6 +25,7 @@
 //
 //------------------------------------------------------------------------------
 
+#include <set>
 #include <span>
 
 #include "MeshKernel/Cartesian3DPoint.hpp"
@@ -1723,4 +1724,75 @@ 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<const Point> boundary)
+    {
+
+        if (boundary.size() < 4)
+        {
+            return false;
+        }
+
+        std::set<Point> uniquePoints;
+
+        for (const Point& p : boundary)
+        {
+            if (!uniquePoints.insert(p).second)
+            {
+                return true;
+            }
+        }
+
+        return false;
+    }
+
+    std::tuple<std::vector<std::vector<Point>>, std::vector<UInt>> splitMultiplePolygons(std::span<const Point> boundaryPoints, std::span<UInt> elementIds)
+    {
+        std::vector<std::vector<Point>> completedPolygons;
+        std::vector<UInt> firstElementIds;
+
+        std::vector<std::pair<Point, int>> pointFaceStack;
+        // Maps a point to its current index in the pointFaceStack
+        std::map<Point, size_t> 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<Point> subPolygon;
+                subPolygon.reserve((pointFaceStack.size() - loopStartIndex) + 1);
+
+                int first_edge_id = pointFaceStack[loopStartIndex].second;
+
+                for (size_t j = loopStartIndex; j < pointFaceStack.size(); ++j)
+                {
+                    subPolygon.push_back(pointFaceStack[j].first);
+                    activePoints.erase(pointFaceStack[j].first);
+                }
+
+                // close the polygon
+                if (!subPolygon.empty())
+                {
+                    subPolygon.push_back(subPolygon.front());
+                }
+
+                completedPolygons.push_back(subPolygon);
+                firstElementIds.push_back(first_edge_id);
+
+                pointFaceStack.resize(loopStartIndex);
+            }
+
+            int currentEdge = (i < elementIds.size()) ? elementIds[i] : -1;
+
+            activePoints[current_point] = pointFaceStack.size();
+            pointFaceStack.emplace_back(current_point, currentEdge);
+        }
+
+        return {completedPolygons, firstElementIds};
+    }
+
 } // namespace meshkernel

libs/MeshKernel/src/Polygon.cpp

@@ -42,6 +42,12 @@ meshkernel::Polygon::Polygon(const std::vector<Point>& points,
     Initialise();
 }
 
+meshkernel::Polygon::Polygon(std::span<const Point> points,
+                             Projection projection) : m_nodes(points.begin(), points.end()), m_projection(projection)
+{
+    Initialise();
+}
+
 meshkernel::Polygon::Polygon(std::vector<Point>&& points,
                              Projection projection) : m_nodes(points), m_projection(projection)
 {