Feature/gridedit 2293 separate polygon loops

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

libs/MeshKernel/include/MeshKernel/Definitions.hpp

@@ -171,4 +171,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

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -254,10 +254,6 @@ namespace meshkernel
         /// @return The resulting polygon mesh boundary
         [[nodiscard]] std::vector<Point> ComputeBoundaryPolygons(const std::vector<Point>& polygon);
 
-        /// @brief Convert all mesh boundaries to a vector of polygon nodes
-        /// @return The resulting set of polygons, describing interior mesh boundaries
-        std::vector<Point> ComputeInnerBoundaryPolygons() const;
-
         /// @brief Gets the hanging edges
         /// @return A vector with the indices of the hanging edges
         [[nodiscard]] std::vector<UInt> GetHangingEdges() const;
@@ -444,25 +440,17 @@ namespace meshkernel
         /// @brief Find the mesh faces that lie entirely within the polygon.
         std::vector<bool> FindFacesEntirelyInsidePolygon(const std::vector<bool>& isNodeInsidePolygon) const;
 
-        /// @brief Constructs a polygon from the meshboundary, by walking through the mesh
-        void WalkBoundaryFromNode(const Polygon& polygon,
-                                  std::vector<bool>& isVisited,
-                                  UInt& currentNode,
-                                  std::vector<Point>& meshBoundaryPolygon) 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 Constructs a polygon or polygons from the meshboundary, by walking through the mesh
+        /// @brief Convert all mesh boundaries to a vector of polygon nodes, including holes (copynetboundstopol)
         ///
-        /// If there are multiple polygons connected by a single node, then these will be separated into individual polygons
-        void WalkMultiBoundaryFromNode(std::vector<bool>& edgeIsVisited,
-                                       std::vector<bool>& nodeIsVisited,
-                                       UInt& currentNode,
-                                       std::vector<Point>& meshBoundaryPolygon,
-                                       std::vector<UInt>& nodeIds,
-                                       std::vector<Point>& subSequence,
-                                       std::vector<Point>& illegalCells) const;
-
-        /// @brief Ensure that all polynomials are orientated in the ACW direction.
-        void OrientatePolygonsAntiClockwise(std::vector<Point>& polygonNodes) const;
+        /// @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 Removes the outer domain boundary polygon from the set of polygons
         ///

libs/MeshKernel/include/MeshKernel/MeshBoundaryExtractor.hpp

@@ -0,0 +1,126 @@
+//---- 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 <http://www.gnu.org/licenses/>.
+//
+// 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 <tuple>
+#include <vector>
+
+#include "MeshKernel/Definitions.hpp"
+#include "MeshKernel/Mesh2D.hpp"
+#include "MeshKernel/Point.hpp"
+#include "MeshKernel/Polygon.hpp"
+
+namespace meshkernel
+{
+
+    /// @brief Extract the boundary polygon from the mesh
+    class MeshBoundaryExtractor
+    {
+    public:
+        /// @brief Extract all boundaries as a single sequence of points, separated by an invalid point
+        static std::vector<Point> ExtractConcatenated(const Mesh2D& mesh, BoundarySelection boundaryType = BoundarySelection::All);
+
+        /// @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<std::vector<Point>>, std::vector<bool>> 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<std::vector<Point>>, std::vector<bool>> 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
+        {
+            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.
+        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<BoundaryEdge>& boundaryEdges,
+                                             const std::vector<bool>& edgeVisited,
+                                             const double incomingAngle);
+
+        /// @brief Construct mapping from node-id to all impinging edges for boundary all edges
+        static void FindAllBoundarEdges(const std::vector<Point>& nodes,
+                                        const std::vector<Edge>& edges,
+                                        const std::vector<std::array<UInt, 2>>& edgesFaces,
+                                        std::unordered_map<UInt, std::vector<BoundaryEdge>>& 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 and clipped to be
+        /// inside a constraining polygon
+        static void Append(const Polygon& polygon,
+                           const Point& centre,
+                           const Projection projection,
+                           std::vector<Point>& boundaryPolygon,
+                           std::vector<bool>& isExterior,
+                           std::vector<std::vector<Point>>& 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<Point>& 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<Point>& nodes,
+                                         const std::vector<Edge>& edges,
+                                         const std::vector<std::array<UInt, 2>>& edgesFaces,
+                                         std::vector<std::vector<Point>>& allPolygons,
+                                         std::vector<std::vector<UInt>>& allTouchedFaces);
+
+        /// @brief Separate polygons that contains multiple sub-polygons and determine externality
+        ///
+        /// allBoundaryPolygons is not const because it may be updated.
+        /// @note allBoundaryPolygons should not be accessed after calling this function
+        static std::tuple<std::vector<std::vector<meshkernel::Point>>, std::vector<bool>>
+        SeparateAndDetermineExternality(const Mesh2D& mesh,
+                                        const Polygon& polygon,
+                                        std::vector<std::vector<Point>>& allBoundaryPolygons,
+                                        const std::vector<std::vector<UInt>>& allTouchedFaces);
+    };
+
+} // namespace meshkernel

libs/MeshKernel/include/MeshKernel/Operations.hpp

@@ -298,6 +298,12 @@ namespace meshkernel
                                          const std::span<const Point> 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<double, Point> ComputePolygonAreaAndCentre(const std::vector<Point>& boundaryPolygon, const Projection& projection);
+
     /// @brief Computes three base components
     void ComputeThreeBaseComponents(const Point& point, std::array<double, 3>& exxp, std::array<double, 3>& eyyp, std::array<double, 3>& ezzp);
 
@@ -621,4 +627,42 @@ 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<const UInt> elementIds);
+
+    /// @brief Concaenate vectors of vectors of points to a vector of points separated by the invalid point
+    template <typename Predicate>
+    std::vector<Point> ConcatenatePointVectors(const std::vector<std::vector<Point>>& pointVectors, Predicate predicate)
+    {
+
+        std::vector<Point> 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

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);
@@ -73,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;
 
@@ -218,6 +226,11 @@ inline meshkernel::UInt meshkernel::Polygon::Size() const
     return static_cast<UInt>(m_nodes.size());
 }
 
+inline bool meshkernel::Polygon::IsEmpty() const
+{
+    return m_nodes.empty();
+}
+
 inline const std::vector<meshkernel::Point>& meshkernel::Polygon::Nodes() const
 {
     return m_nodes;