Create a MeshKernelAPI for generating curvilinear grid boundary polyg…

…ons (#351 | GRIDEDIT-1317)

libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp

@@ -27,6 +27,7 @@
 
 #pragma once
 
+#include <set>
 #include <vector>
 
 #include <MeshKernel/BoundingBox.hpp>
@@ -37,7 +38,6 @@
 #include <MeshKernel/CurvilinearGrid/UndoActions/CurvilinearGridBlockUndoAction.hpp>
 #include <MeshKernel/CurvilinearGrid/UndoActions/CurvilinearGridRefinementUndoAction.hpp>
 #include <MeshKernel/CurvilinearGrid/UndoActions/ResetCurvilinearNodeAction.hpp>
-#include <MeshKernel/Entities.hpp>
 #include <MeshKernel/Exceptions.hpp>
 #include <MeshKernel/Mesh.hpp>
 #include <MeshKernel/UndoActions/UndoAction.hpp>
@@ -56,6 +56,9 @@ namespace meshkernel
         /// @brief Typedef for face indices
         using CurvilinearFaceNodeIndices = std::array<CurvilinearGridNodeIndices, 4>;
 
+        /// @brief Typedef for defining a curvilinear edge
+        using CurvilinearEdge = std::pair<CurvilinearGridNodeIndices, CurvilinearGridNodeIndices>;
+
         /// @brief An enum for boundary grid line types
         enum class BoundaryGridLineType
         {
@@ -308,13 +311,13 @@ namespace meshkernel
         /// @brief Computes the node from an index
         [[nodiscard]] const Point& Node(const UInt index) const
         {
-            const auto nodeIndex = GetNodeIndex(index);
+            const auto nodeIndex = GetCurvilinearGridNodeIndices(index);
             return GetNode(nodeIndex);
         }
 
         /// @brief Get the node CurvilinearGridNodeIndices from a position
-        /// @return The CurvilinearGridNodeIndices
-        CurvilinearGridNodeIndices GetNodeIndex(const UInt index) const
+        /// @return The CurvilinearGridNodeIndices for the position
+        CurvilinearGridNodeIndices GetCurvilinearGridNodeIndices(const UInt index) const
         {
             if (index >= NumM() * NumN())
             {
@@ -326,6 +329,18 @@ namespace meshkernel
             return {n, m};
         }
 
+        /// @brief Get the node index from a CurvilinearGridNodeIndices
+        /// @return The CurvilinearGridNodeIndices
+        UInt GetNodeIndex(const CurvilinearGridNodeIndices& nodeIndex) const
+        {
+            if (!nodeIndex.IsValid())
+            {
+                throw AlgorithmError("CurvilinearGridNodeIndices ");
+            }
+
+            return nodeIndex.m_m + nodeIndex.m_n * NumM();
+        }
+
         /// @brief Finds the index of the closest location
         /// @param[in] point the input point
         /// @param[in] location the location type
@@ -389,6 +404,12 @@ namespace meshkernel
         /// @return a vector of M nodes
         std::vector<Point> GetNodeVectorAtN(UInt n) const { return lin_alg::MatrixColToSTLVector(m_gridNodes, n + m_startOffset.m_n); }
 
+        /// @brief Compute a sequence of one or more outer boundary polygons separated by geometry separators
+        /// @param[in]  lowerLeft   The node index of the lower left corner
+        /// @param[in]  upperRight  The node index of the upper right corner
+        /// @returns The vector containing the boundary polygon points
+        [[nodiscard]] std::vector<Point> ComputeBoundaryPolygons(const CurvilinearGridNodeIndices& lowerLeft, const CurvilinearGridNodeIndices& upperRight) const;
+
         /// @brief The number of nodes M in the m dimension
         /// @return A number >= 2 for a valid curvilinear grid
         UInt FullNumM() const { return static_cast<UInt>(m_gridNodes.cols()); }
@@ -454,9 +475,11 @@ namespace meshkernel
         /// @returns The  mapping (m and n indices for each node of the edge)
         [[nodiscard]] std::vector<CurvilinearEdgeNodeIndices> ComputeEdgeIndices() const;
 
-        /// @brief Compute the face indices.
-        /// @returns The  mapping (m and n indices for each node of the face)
-        [[nodiscard]] std::vector<CurvilinearFaceNodeIndices> ComputeFaceIndices() const;
+        /// @brief Compute the boundary edges. While iterating over edges of valid faces, boundary edges are seen once, all internal edges are seen by two faces.
+        /// @param[in]  lowerLeft   The node index of the lower left corner
+        /// @param[in]  upperRight  The node index of the upper right corner
+        /// @returns The set of boundary edges
+        [[nodiscard]] std::set<CurvilinearEdge> ComputeBoundaryEdges(const CurvilinearGridNodeIndices& lowerLeft, const CurvilinearGridNodeIndices& upperRight) const;
 
         /// @brief Set the m_nodesRTreeRequiresUpdate flag
         /// @param[in] value The value of the flag
@@ -536,11 +559,11 @@ meshkernel::Point const& meshkernel::CurvilinearGrid::GetNode(const UInt n, cons
     return m_gridNodes(n + m_startOffset.m_n, m + m_startOffset.m_m);
 }
 
-meshkernel::Point& meshkernel::CurvilinearGrid::GetNode(const meshkernel::CurvilinearGridNodeIndices& index)
+meshkernel::Point& meshkernel::CurvilinearGrid::GetNode(const CurvilinearGridNodeIndices& index)
 {
     if (!index.IsValid()) [[unlikely]]
     {
-        throw meshkernel::ConstraintError("Invalid node index");
+        throw ConstraintError("Invalid node index");
     }
     m_nodesRTreeRequiresUpdate = true;
     m_edgesRTreeRequiresUpdate = true;
@@ -549,11 +572,11 @@ meshkernel::Point& meshkernel::CurvilinearGrid::GetNode(const meshkernel::Curvil
     return GetNode(index.m_n, index.m_m);
 }
 
-meshkernel::Point const& meshkernel::CurvilinearGrid::GetNode(const meshkernel::CurvilinearGridNodeIndices& index) const
+meshkernel::Point const& meshkernel::CurvilinearGrid::GetNode(const CurvilinearGridNodeIndices& index) const
 {
     if (!index.IsValid()) [[unlikely]]
     {
-        throw meshkernel::ConstraintError("Invalid node index");
+        throw ConstraintError("Invalid node index");
     }
 
     return GetNode(index.m_n, index.m_m);

libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridNodeIndices.hpp

@@ -40,7 +40,7 @@ namespace meshkernel
         /// @brief Default constructor sets the indices to invalid
         CurvilinearGridNodeIndices() : m_n(constants::missing::uintValue), m_m(constants::missing::uintValue) {}
 
-        /// @brief Constructor sets indices from values
+        /// @brief Constructor setting n and m indices
         /// @param[in] m The m index
         /// @param[in] n The n index
         CurvilinearGridNodeIndices(UInt n, UInt m) : m_n(n), m_m(m) {}
@@ -51,9 +51,6 @@ namespace meshkernel
             return m_m != missingValue && m_n != missingValue;
         }
 
-        /// @brief Overloads equality with another CurvilinearGridNodeIndices
-        bool operator==(const CurvilinearGridNodeIndices& rhs) const = default;
-
         CurvilinearGridNodeIndices& operator+=(const CurvilinearGridNodeIndices& val)
         {
             if (!IsValid())
@@ -71,6 +68,19 @@ namespace meshkernel
             return *this;
         }
 
+        /// @brief Overloads <=> operator
+        auto operator<=>(const CurvilinearGridNodeIndices& rhs) const
+        {
+            if (m_n == rhs.m_n)
+            {
+                return m_m <=> rhs.m_m; // Return result of m_m comparison if not equal
+            }
+            return m_n <=> rhs.m_n; // Compare m_n if m_m values are equal
+        }
+
+        /// @brief Overloads equality with another CurvilinearGridNodeIndices
+        bool operator==(const CurvilinearGridNodeIndices& rhs) const = default;
+
         CurvilinearGridNodeIndices& operator-=(const CurvilinearGridNodeIndices& val)
         {
             if (!IsValid())

libs/MeshKernel/include/MeshKernel/Mesh2D.hpp

@@ -257,7 +257,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
         /// @return The resulting polygon mesh boundary
-        [[nodiscard]] std::vector<Point> MeshBoundaryToPolygon(const std::vector<Point>& polygon);
+        [[nodiscard]] std::vector<Point> ComputeBoundaryPolygons(const std::vector<Point>& polygon);
 
         /// @brief Constructs a polygon from the meshboundary, by walking through the mesh
         /// @param[in] polygon The input polygon

libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp

@@ -1283,41 +1283,118 @@ std::vector<CurvilinearGrid::CurvilinearEdgeNodeIndices> CurvilinearGrid::Comput
     return result;
 }
 
-std::vector<CurvilinearGrid::CurvilinearFaceNodeIndices> CurvilinearGrid::ComputeFaceIndices() const
+std::set<CurvilinearGrid::CurvilinearEdge> CurvilinearGrid::ComputeBoundaryEdges(const CurvilinearGridNodeIndices& lowerLeft, const CurvilinearGridNodeIndices& upperRight) const
 {
-    const auto numFaces = (NumM() - 1) * (NumN() - 1);
+    std::vector<Point> result;
+    std::set<CurvilinearEdge> boundaryEdges;
+    for (UInt n = lowerLeft.m_n; n < upperRight.m_n; n++)
+    {
+        for (UInt m = lowerLeft.m_m; m < upperRight.m_m; m++)
+        {
+            CurvilinearFaceNodeIndices faceIndices;
 
-    std::vector<CurvilinearFaceNodeIndices> result(numFaces);
+            faceIndices[0] = {n, m};
+            faceIndices[1] = {n, m + 1};
+            faceIndices[2] = {n + 1, m + 1};
+            faceIndices[3] = {n + 1, m};
 
-    UInt index = 0;
-    for (UInt n = 0; n < NumN() - 1; n++)
+            if (!std::ranges::all_of(faceIndices, [this](const CurvilinearGridNodeIndices& curviNode)
+                                     { return GetNode(curviNode.m_n, curviNode.m_m).IsValid(); }))
+            {
+                continue;
+            }
+
+            for (UInt i = 0u; i < constants::geometric::numNodesInQuadrilateral; ++i)
+            {
+                const auto firstCurvilinearNodeIndex = faceIndices[i];
+                const auto nextIndex = NextCircularForwardIndex(i, constants::geometric::numNodesInQuadrilateral);
+                const auto secondCurvilinearNodeIndex = faceIndices[nextIndex];
+
+                const auto edge = std::make_pair(std::min(firstCurvilinearNodeIndex, secondCurvilinearNodeIndex),
+                                                 std::max(firstCurvilinearNodeIndex, secondCurvilinearNodeIndex));
+
+                const auto it = boundaryEdges.find(edge);
+                if (it != boundaryEdges.end())
+                {
+                    boundaryEdges.erase(it);
+                }
+                else
+                {
+                    boundaryEdges.insert(edge);
+                }
+            }
+        }
+    }
+
+    return boundaryEdges;
+}
+
+std::vector<meshkernel::Point> CurvilinearGrid::ComputeBoundaryPolygons(const CurvilinearGridNodeIndices& lowerLeft, const CurvilinearGridNodeIndices& upperRight) const
+{
+    std::vector<Point> result;
+
+    auto boundaryEdges = ComputeBoundaryEdges(lowerLeft, upperRight);
+    if (boundaryEdges.empty())
     {
-        for (UInt m = 0; m < NumM() - 1; m++)
+        return result;
+    }
+
+    auto currentEdge = boundaryEdges.begin();
+    auto startNode = currentEdge->first;
+    auto sharedNode = currentEdge->second;
+
+    // iterate over all boundary edges, until all boundary polygons are found
+    while (!boundaryEdges.empty())
+    {
+        result.push_back(GetNode(startNode.m_n, startNode.m_m));
+        // iterate over connected edges, until the start node is found and the current boundary polygon is closed
+        while (sharedNode != startNode)
         {
-            result[index][0] = {n, m};
-            result[index][1] = {n, m + 1};
-            result[index][2] = {n + 1, m + 1};
-            result[index][3] = {n + 1, m};
-            index++;
+            result.push_back(GetNode(sharedNode.m_n, sharedNode.m_m));
+
+            boundaryEdges.erase(currentEdge);
+            currentEdge = std::ranges::find_if(boundaryEdges, [&](const CurvilinearEdge& edge)
+                                               { return edge.first == sharedNode ||
+                                                        edge.second == sharedNode; });
+
+            sharedNode = currentEdge->first == sharedNode ? currentEdge->second : currentEdge->first;
+        }
+        result.push_back(GetNode(startNode.m_n, startNode.m_m));
+        boundaryEdges.erase(currentEdge);
+        if (boundaryEdges.empty())
+        {
+            return result;
         }
+
+        result.emplace_back(constants::missing::doubleValue, constants::missing::doubleValue);
+
+        // startNode node for the next boundary polygon
+        currentEdge = boundaryEdges.begin();
+        startNode = currentEdge->first;
+        sharedNode = currentEdge->second;
     }
+
     return result;
 }
 
 std::vector<meshkernel::Point> CurvilinearGrid::ComputeFaceCenters() const
 {
-    const auto faceIndices = ComputeFaceIndices();
+    std::vector<Point> result;
+    result.reserve((NumN() - 1) * (NumM() - 1));
 
-    std::vector<Point> result(faceIndices.size());
-
-    for (UInt i = 0; i < faceIndices.size(); ++i)
+    for (UInt n = 0; n < NumN() - 1; n++)
     {
-        Point massCenter{0.0, 0.0};
-        for (const auto& index : faceIndices[i])
+        for (UInt m = 0; m < NumM() - 1; m++)
         {
-            massCenter += GetNode(index.m_n, index.m_m);
+            Point massCenter{0.0, 0.0};
+
+            massCenter += GetNode(n, m);
+            massCenter += GetNode(n, m + 1);
+            massCenter += GetNode(n + 1, m + 1);
+            massCenter += GetNode(n + 1, m);
+
+            result.push_back(massCenter * 0.25);
         }
-        result[i] = massCenter * 0.25;
     }
     return result;
 }

libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineShift.cpp

@@ -161,7 +161,7 @@ meshkernel::UndoActionPtr CurvilinearGridLineShift::MoveNode(Point const& fromPo
     // Get the index in the grid of the line to be shifted
 
     auto const nodePosition = m_grid.FindLocationIndex(fromPoint, Location::Nodes);
-    auto const nodeIndex = m_grid.GetNodeIndex(nodePosition);
+    auto const nodeIndex = m_grid.GetCurvilinearGridNodeIndices(nodePosition);
 
     // Check the nodes are on the line to shift
     if (!m_lines[0].IsNodeOnLine(nodeIndex))

libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSnapping.cpp

@@ -42,7 +42,7 @@ void meshkernel::CurvilinearGridSnapping::Initialise()
         {
 
             auto const extentIndexPosition = m_grid.FindLocationIndex(m_controlPoints[2], Location::Nodes);
-            CurvilinearGridNodeIndices extentIndex = m_grid.GetNodeIndex(extentIndexPosition);
+            CurvilinearGridNodeIndices extentIndex = m_grid.GetCurvilinearGridNodeIndices(extentIndexPosition);
 
             m_indexBoxLowerLeft = CurvilinearGridNodeIndices(std::min({m_lineStartIndex.m_n, m_lineEndIndex.m_n, extentIndex.m_n}),
                                                              std::min({m_lineStartIndex.m_m, m_lineEndIndex.m_m, extentIndex.m_m}));

libs/MeshKernel/src/LandBoundaries.cpp

@@ -64,7 +64,7 @@ void LandBoundaries::Administrate()
 
     // Mesh boundary to polygon
     const std::vector<Point> polygonNodes;
-    const auto meshBoundaryPolygon = m_mesh.MeshBoundaryToPolygon(polygonNodes);
+    const auto meshBoundaryPolygon = m_mesh.ComputeBoundaryPolygons(polygonNodes);
 
     // Find the start/end node of the land boundaries
     const auto landBoundaryIndices = m_landBoundary.FindPolylineIndices();

libs/MeshKernel/src/Mesh2D.cpp

@@ -1505,7 +1505,7 @@ std::vector<meshkernel::UInt> Mesh2D::SortedFacesAroundNode(UInt node) const
     return result;
 }
 
-std::vector<meshkernel::Point> Mesh2D::MeshBoundaryToPolygon(const std::vector<Point>& polygonNodes)
+std::vector<meshkernel::Point> Mesh2D::ComputeBoundaryPolygons(const std::vector<Point>& polygonNodes)
 {
     Polygon polygon(polygonNodes, m_projection);
 

libs/MeshKernel/tests/CMakeLists.txt

@@ -16,6 +16,7 @@ set(
   ${SRC_DIR}/CompoundUndoTests.cpp
   ${SRC_DIR}/ContactsTests.cpp
   ${SRC_DIR}/CurvilinearGridBasicTests.cpp
+  ${SRC_DIR}/CurvilinearGridBoundaryPolygonTests.cpp
   ${SRC_DIR}/CurvilinearGridCurvatureTests.cpp
   ${SRC_DIR}/CurvilinearGridDeRefinementTests.cpp
   ${SRC_DIR}/CurvilinearGridFromPolygonTests.cpp

libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp

@@ -42,7 +42,6 @@
 #include "MeshKernel/CurvilinearGrid/CurvilinearGridOrthogonalization.hpp"
 #include "MeshKernel/CurvilinearGrid/CurvilinearGridRefinement.hpp"
 #include "MeshKernel/CurvilinearGrid/CurvilinearGridSmoothing.hpp"
-#include "MeshKernel/Entities.hpp"
 #include "MeshKernel/Operations.hpp"
 #include "MeshKernel/UndoActions/UndoActionStack.hpp"
 #include "MeshKernel/Utilities/LinearAlgebra.hpp"