From 1e6996fecb8d158b31aceb3881a5b95afcaa735b Mon Sep 17 00:00:00 2001 From: Luca Carniato Date: Wed, 24 Jul 2024 11:35:46 +0200 Subject: [PATCH] Create a MeshKernelAPI for generating curvilinear grid boundary polygons (#351 | GRIDEDIT-1317) --- .../CurvilinearGrid/CurvilinearGrid.hpp | 45 +++-- .../CurvilinearGridNodeIndices.hpp | 18 +- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 2 +- .../src/CurvilinearGrid/CurvilinearGrid.cpp | 115 ++++++++++-- .../CurvilinearGridLineShift.cpp | 2 +- .../CurvilinearGridSnapping.cpp | 2 +- libs/MeshKernel/src/LandBoundaries.cpp | 2 +- libs/MeshKernel/src/Mesh2D.cpp | 2 +- libs/MeshKernel/tests/CMakeLists.txt | 1 + .../tests/src/CurvilinearGridBasicTests.cpp | 1 - .../CurvilinearGridBoundaryPolygonTests.cpp | 173 ++++++++++++++++++ .../src/CurvilinearGridSnappingTests.cpp | 1 - libs/MeshKernel/tests/src/Mesh2DTest.cpp | 2 +- libs/MeshKernel/tests/src/MeshTests.cpp | 2 +- libs/MeshKernel/tests/src/PolygonsTests.cpp | 2 +- .../include/MeshKernelApi/MeshKernel.hpp | 22 +++ libs/MeshKernelApi/src/MeshKernel.cpp | 81 +++++++- .../tests/src/CurvilinearGridTests.cpp | 99 ++++++++-- 18 files changed, 506 insertions(+), 66 deletions(-) create mode 100644 libs/MeshKernel/tests/src/CurvilinearGridBoundaryPolygonTests.cpp diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp index 2d99fcf4c..490b24222 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp @@ -27,6 +27,7 @@ #pragma once +#include #include #include @@ -37,7 +38,6 @@ #include #include #include -#include #include #include #include @@ -56,6 +56,9 @@ namespace meshkernel /// @brief Typedef for face indices using CurvilinearFaceNodeIndices = std::array; + /// @brief Typedef for defining a curvilinear edge + using CurvilinearEdge = std::pair; + /// @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 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 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(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 ComputeEdgeIndices() const; - /// @brief Compute the face indices. - /// @returns The mapping (m and n indices for each node of the face) - [[nodiscard]] std::vector 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 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); diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridNodeIndices.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridNodeIndices.hpp index 8e68ed079..a10060add 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridNodeIndices.hpp +++ b/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()) diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 2f0751ee9..e3b30fc80 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/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 MeshBoundaryToPolygon(const std::vector& polygon); + [[nodiscard]] std::vector ComputeBoundaryPolygons(const std::vector& polygon); /// @brief Constructs a polygon from the meshboundary, by walking through the mesh /// @param[in] polygon The input polygon diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp index 2ce222477..55cb4c76c 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp @@ -1283,41 +1283,118 @@ std::vector CurvilinearGrid::Comput return result; } -std::vector CurvilinearGrid::ComputeFaceIndices() const +std::set CurvilinearGrid::ComputeBoundaryEdges(const CurvilinearGridNodeIndices& lowerLeft, const CurvilinearGridNodeIndices& upperRight) const { - const auto numFaces = (NumM() - 1) * (NumN() - 1); + std::vector result; + std::set 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 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 CurvilinearGrid::ComputeBoundaryPolygons(const CurvilinearGridNodeIndices& lowerLeft, const CurvilinearGridNodeIndices& upperRight) const +{ + std::vector 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 CurvilinearGrid::ComputeFaceCenters() const { - const auto faceIndices = ComputeFaceIndices(); + std::vector result; + result.reserve((NumN() - 1) * (NumM() - 1)); - std::vector 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; } diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineShift.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineShift.cpp index a1f796159..e744f87c5 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridLineShift.cpp +++ b/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)) diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSnapping.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSnapping.cpp index 456d7c462..896e72aa5 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSnapping.cpp +++ b/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})); diff --git a/libs/MeshKernel/src/LandBoundaries.cpp b/libs/MeshKernel/src/LandBoundaries.cpp index d415532c2..d5d198be5 100644 --- a/libs/MeshKernel/src/LandBoundaries.cpp +++ b/libs/MeshKernel/src/LandBoundaries.cpp @@ -64,7 +64,7 @@ void LandBoundaries::Administrate() // Mesh boundary to polygon const std::vector 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(); diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index fe234e802..6eda77622 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -1505,7 +1505,7 @@ std::vector Mesh2D::SortedFacesAroundNode(UInt node) const return result; } -std::vector Mesh2D::MeshBoundaryToPolygon(const std::vector& polygonNodes) +std::vector Mesh2D::ComputeBoundaryPolygons(const std::vector& polygonNodes) { Polygon polygon(polygonNodes, m_projection); diff --git a/libs/MeshKernel/tests/CMakeLists.txt b/libs/MeshKernel/tests/CMakeLists.txt index 6875e9fcb..f9942a055 100644 --- a/libs/MeshKernel/tests/CMakeLists.txt +++ b/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 diff --git a/libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp index 97a0a979b..def12db84 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp +++ b/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" diff --git a/libs/MeshKernel/tests/src/CurvilinearGridBoundaryPolygonTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridBoundaryPolygonTests.cpp new file mode 100644 index 000000000..cdd0891b9 --- /dev/null +++ b/libs/MeshKernel/tests/src/CurvilinearGridBoundaryPolygonTests.cpp @@ -0,0 +1,173 @@ +#include "MakeCurvilinearGrids.hpp" +#include "MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp" + +#include +#include + +using namespace meshkernel; + +TEST(CurvilinearGridBoundaryPolygon, ComputeBoundaryToPolygon_OnValidCurvilinearGrid_ShouldComputeBoundaryPolygon) + +{ + // Prepare + constexpr double originX = 0.0; + constexpr double originY = 0.0; + + constexpr double deltaX = 1.0; + constexpr double deltaY = 1.0; + + constexpr UInt nx = 3; + constexpr UInt ny = 3; + + std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); + CurvilinearGridNodeIndices lowerLeft(0, 0); + CurvilinearGridNodeIndices upperRight(ny - 1, nx - 1); + + // Execute + + const auto boundaryPolygon = grid->ComputeBoundaryPolygons(lowerLeft, upperRight); + + // Assert + std::vector expectedBoundaryPolygon{{0.0, 0.0}, + {1.0, 0.0}, + {2.0, 0.0}, + {2.0, 1.0}, + {2.0, 2.0}, + {1.0, 2.0}, + {0.0, 2.0}, + {0.0, 1.0}, + {0.0, 0.0}}; + ASSERT_THAT(boundaryPolygon, ::testing::ContainerEq(expectedBoundaryPolygon)); +} + +TEST(CurvilinearGridBoundaryPolygon, ComputeBoundaryToPolygon_OnValidCurvilinearGridWithInvalidNodes_ShouldComputeBoundaryPolygon) + +{ + // Prepare + constexpr double originX = 0.0; + constexpr double originY = 0.0; + + constexpr double deltaX = 1.0; + constexpr double deltaY = 1.0; + + constexpr UInt nx = 4; + constexpr UInt ny = 4; + + std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); + auto undoAction = grid->DeleteNode({3.0, 0.0}); + undoAction = grid->DeleteNode({3.0, 3.0}); + CurvilinearGridNodeIndices lowerLeft(0, 0); + CurvilinearGridNodeIndices upperRight(ny - 1, nx - 1); + + // Execute + const auto boundaryPolygon = grid->ComputeBoundaryPolygons(lowerLeft, upperRight); + + // Assert + std::vector expectedBoundaryPolygon{{0.0, 0.0}, + {1.0, 0.0}, + {2.0, 0.0}, + {2.0, 1.0}, + {3.0, 1.0}, + {3.0, 2.0}, + {2.0, 2.0}, + {2.0, 3.0}, + {1.0, 3.0}, + {0.0, 3.0}, + {0.0, 2.0}, + {0.0, 1.0}, + {0.0, 0.0}}; + + ASSERT_THAT(boundaryPolygon, ::testing::ContainerEq(expectedBoundaryPolygon)); +} + +TEST(CurvilinearGridBoundaryPolygon, ComputeBoundaryToPolygon_OnValidCurvilinearGridWithDisconnectedCurvilinearGrid_ShouldComputeBoundaryPolygons) + +{ + // Prepare + constexpr double originX = 0.0; + constexpr double originY = 0.0; + + constexpr double deltaX = 1.0; + constexpr double deltaY = 1.0; + + constexpr UInt nx = 9; + constexpr UInt ny = 3; + + std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); + auto undoAction = grid->DeleteNode({3.0, 1.0}); + undoAction = grid->DeleteNode({4.0, 1.0}); + CurvilinearGridNodeIndices lowerLeft(0, 0); + CurvilinearGridNodeIndices upperRight(ny - 1, nx - 1); + + // Execute + const auto boundaryPolygon = grid->ComputeBoundaryPolygons(lowerLeft, upperRight); + + // Assert + std::vector expectedBoundaryPolygon{{0.0, 0.0}, + {1.0, 0.0}, + {2.0, 0.0}, + {2.0, 1.0}, + {2.0, 2.0}, + {1.0, 2.0}, + {0.0, 2.0}, + {0.0, 1.0}, + {0.0, 0.0}, + {-999.0, -999.0}, + {5.0, 0.0}, + {6.0, 0.0}, + {7.0, 0.0}, + {8.0, 0.0}, + {8.0, 1.0}, + {8.0, 2.0}, + {7.0, 2.0}, + {6.0, 2.0}, + {5.0, 2.0}, + {5.0, 1.0}, + {5.0, 0.0}}; + + ASSERT_THAT(boundaryPolygon, ::testing::ContainerEq(expectedBoundaryPolygon)); +} + +TEST(CurvilinearGridBoundaryPolygon, ComputeBoundaryToPolygon_OnValidCurvilinearGridWithDisconnectedCurvilinearGridAndSelection_ShouldComputeBoundaryPolygons) + +{ + // Prepare + constexpr double originX = 0.0; + constexpr double originY = 0.0; + + constexpr double deltaX = 1.0; + constexpr double deltaY = 1.0; + + constexpr UInt nx = 9; + constexpr UInt ny = 3; + + std::unique_ptr grid = MakeCurvilinearGrid(originX, originY, deltaX, deltaY, nx, ny); + auto undoAction = grid->DeleteNode({3.0, 1.0}); + undoAction = grid->DeleteNode({4.0, 1.0}); + CurvilinearGridNodeIndices lowerLeft(0, 0); + CurvilinearGridNodeIndices upperRight(1, nx - 1); // not the entire clg + + // Execute + const auto boundaryPolygon = grid->ComputeBoundaryPolygons(lowerLeft, upperRight); + + // Assert + std::vector expectedBoundaryPolygon{{0.0, 0.0}, + {1.0, 0.0}, + {2.0, 0.0}, + {2.0, 1.0}, + {1.0, 1.0}, + {0.0, 1.0}, + {0.0, 0.0}, + {-999.0, -999.0}, + {5.0, 0.0}, + {6.0, 0.0}, + {7.0, 0.0}, + {8.0, 0.0}, + {8.0, 1.0}, + {7.0, 1.0}, + {6.0, 1.0}, + {5.0, 1.0}, + {5.0, 0.0}}; + + ASSERT_THAT(boundaryPolygon, ::testing::ContainerEq(expectedBoundaryPolygon)); +} diff --git a/libs/MeshKernel/tests/src/CurvilinearGridSnappingTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridSnappingTests.cpp index d963b11fe..b46848ed4 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridSnappingTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridSnappingTests.cpp @@ -4,7 +4,6 @@ #include #include #include -#include #include #include diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index 95452459c..47d6f2ad3 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -194,7 +194,7 @@ TEST(Mesh2D, MeshBoundaryToPolygon) auto mesh = meshkernel::Mesh2D(edges, nodes, meshkernel::Projection::cartesian); std::vector polygonNodes; - const auto meshBoundaryPolygon = mesh.MeshBoundaryToPolygon(polygonNodes); + const auto meshBoundaryPolygon = mesh.ComputeBoundaryPolygons(polygonNodes); const double tolerance = 1e-5; ASSERT_NEAR(0.0, meshBoundaryPolygon[0].x, tolerance); diff --git a/libs/MeshKernel/tests/src/MeshTests.cpp b/libs/MeshKernel/tests/src/MeshTests.cpp index 18ca86e05..6efd778f0 100644 --- a/libs/MeshKernel/tests/src/MeshTests.cpp +++ b/libs/MeshKernel/tests/src/MeshTests.cpp @@ -177,7 +177,7 @@ TEST(Mesh, MeshBoundaryToPolygon) auto mesh = meshkernel::Mesh2D(edges, nodes, meshkernel::Projection::cartesian); std::vector polygonNodes; - const auto meshBoundaryPolygon = mesh.MeshBoundaryToPolygon(polygonNodes); + const auto meshBoundaryPolygon = mesh.ComputeBoundaryPolygons(polygonNodes); const double tolerance = 1e-5; ASSERT_NEAR(0.0, meshBoundaryPolygon[0].x, tolerance); diff --git a/libs/MeshKernel/tests/src/PolygonsTests.cpp b/libs/MeshKernel/tests/src/PolygonsTests.cpp index 70449915b..c50e1bab1 100644 --- a/libs/MeshKernel/tests/src/PolygonsTests.cpp +++ b/libs/MeshKernel/tests/src/PolygonsTests.cpp @@ -17,7 +17,7 @@ TEST(Polygons, MeshBoundaryToPolygon) auto mesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/SmallTriangularGrid_net.nc"); std::vector polygonNodes; - const auto meshBoundaryPolygon = mesh->MeshBoundaryToPolygon(polygonNodes); + const auto meshBoundaryPolygon = mesh->ComputeBoundaryPolygons(polygonNodes); ASSERT_EQ(8, meshBoundaryPolygon.size()); diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index a21fff069..f23dfcecc 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -305,6 +305,28 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_curvilinear_get_data(int meshKernelId, CurvilinearGrid& curvilinearGrid); + /// @brief Gets the boundary polygon of a curvilinear grid, nodes with invalid coordinates are excluded + /// + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] lowerLeftN The n index of the lower left corner + /// @param[in] lowerLeftM The m index of the lower left corner + /// @param[in] upperRightN The n index of the upper right corner + /// @param[in] upperRightM The m index of the upper right corner + /// @param[out] boundaryPolygons The geometry list containing the boundary polygons + /// @returns Error code + MKERNEL_API int mkernel_curvilinear_get_boundaries_as_polygons(int meshKernelId, int lowerLeftN, int lowerLeftM, int upperRightN, int upperRightM, GeometryList& boundaryPolygons); + + /// @brief Count the number of nodes in curvilinear grid boundary polygons. + /// + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] lowerLeftN The n index of the lower left corner + /// @param[in] lowerLeftM The m index of the lower left corner + /// @param[in] upperRightN The n index of the upper right corner + /// @param[in] upperRightM The m index of the upper right corner + /// @param[out] numberOfPolygonNodes The number of polygon nodes + /// @returns Error code + MKERNEL_API int mkernel_curvilinear_count_boundaries_as_polygons(int meshKernelId, int lowerLeftN, int lowerLeftM, int upperRightN, int upperRightM, int& numberOfPolygonNodes); + /// @brief Gets the curvilinear grid dimensions as a CurvilinearGrid struct (converted as set of edges and nodes). /// /// The integer parameters of the CurvilinearGrid struct are set to the corresponding dimensions diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index df6802566..157b13dfd 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -687,6 +687,83 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_curvilinear_get_boundaries_as_polygons(int meshKernelId, int lowerLeftN, int lowerLeftM, int upperRightN, int upperRightM, GeometryList& boundaryPolygons) + { + lastExitCode = meshkernel::ExitCode::Success; + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + if (!meshKernelState[meshKernelId].m_curvilinearGrid->IsValid()) + { + throw meshkernel::MeshKernelError("Invalid curvilinear grid"); + } + + auto lowerLeftNUnsigned = static_cast(lowerLeftN); + auto lowerLeftMUnsigned = static_cast(lowerLeftM); + auto upperRightNUnsigned = static_cast(upperRightN); + auto upperRightMUnsigned = static_cast(upperRightM); + + const auto minN = std::min(lowerLeftNUnsigned, upperRightNUnsigned); + const auto maxN = std::max(lowerLeftNUnsigned, upperRightNUnsigned); + const auto minM = std::min(lowerLeftMUnsigned, upperRightMUnsigned); + const auto maxM = std::max(lowerLeftMUnsigned, upperRightMUnsigned); + + const auto boundaryPolygon = meshKernelState[meshKernelId].m_curvilinearGrid->ComputeBoundaryPolygons({minN, minM}, + {maxN, maxM}); + ConvertPointVectorToGeometryList(boundaryPolygon, boundaryPolygons); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_curvilinear_count_boundaries_as_polygons(int meshKernelId, int lowerLeftN, int lowerLeftM, int upperRightN, int upperRightM, int& numberOfPolygonNodes) + { + lastExitCode = meshkernel::ExitCode::Success; + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + if (!meshKernelState[meshKernelId].m_curvilinearGrid->IsValid()) + { + throw meshkernel::MeshKernelError("Invalid curvilinear grid"); + } + + const auto lowerLeftNUnsigned = static_cast(lowerLeftN); + const auto lowerLeftMUnsigned = static_cast(lowerLeftM); + const auto upperRightNUnsigned = static_cast(upperRightN); + const auto upperRightMUnsigned = static_cast(upperRightM); + + const auto minN = std::min(lowerLeftNUnsigned, upperRightNUnsigned); + const auto maxN = std::max(lowerLeftNUnsigned, upperRightNUnsigned); + const auto minM = std::min(lowerLeftMUnsigned, upperRightMUnsigned); + const auto maxM = std::max(lowerLeftMUnsigned, upperRightMUnsigned); + + const auto boundaryPolygon = meshKernelState[meshKernelId].m_curvilinearGrid->ComputeBoundaryPolygons({minN, minM}, + {maxN, maxM}); + numberOfPolygonNodes = static_cast(boundaryPolygon.size()); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + MKERNEL_API int mkernel_contacts_get_dimensions(int meshKernelId, Contacts& contacts) { lastExitCode = meshkernel::ExitCode::Success; @@ -1425,7 +1502,7 @@ namespace meshkernelapi } const std::vector polygonNodes; - const auto meshBoundaryPolygon = meshKernelState[meshKernelId].m_mesh2d->MeshBoundaryToPolygon(polygonNodes); + const auto meshBoundaryPolygon = meshKernelState[meshKernelId].m_mesh2d->ComputeBoundaryPolygons(polygonNodes); ConvertPointVectorToGeometryList(meshBoundaryPolygon, boundaryPolygons); } @@ -1447,7 +1524,7 @@ namespace meshkernelapi } const std::vector polygonNodes; - const auto meshBoundaryPolygon = meshKernelState[meshKernelId].m_mesh2d->MeshBoundaryToPolygon(polygonNodes); + const auto meshBoundaryPolygon = meshKernelState[meshKernelId].m_mesh2d->ComputeBoundaryPolygons(polygonNodes); numberOfPolygonNodes = static_cast(meshBoundaryPolygon.size()); // last value is a separator } catch (...) diff --git a/libs/MeshKernelApi/tests/src/CurvilinearGridTests.cpp b/libs/MeshKernelApi/tests/src/CurvilinearGridTests.cpp index 0203a2f7a..3c7a6ac91 100644 --- a/libs/MeshKernelApi/tests/src/CurvilinearGridTests.cpp +++ b/libs/MeshKernelApi/tests/src/CurvilinearGridTests.cpp @@ -1203,7 +1203,7 @@ TEST(CurvilinearGrid, SnapToLandBoundary) auto errorCode = meshkernelapi::mkernel_allocate_state(0, meshKernelId); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); - meshkernel::MakeGridParameters makeGridParameters; + MakeGridParameters makeGridParameters; makeGridParameters.num_columns = 10; makeGridParameters.num_rows = 10; @@ -1217,14 +1217,14 @@ TEST(CurvilinearGrid, SnapToLandBoundary) ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); meshkernelapi::CurvilinearGrid curvilinearGrid{}; - errorCode = meshkernelapi::mkernel_curvilinear_get_dimensions(meshKernelId, curvilinearGrid); + errorCode = mkernel_curvilinear_get_dimensions(meshKernelId, curvilinearGrid); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); std::vector node_x(curvilinearGrid.num_m * curvilinearGrid.num_n); std::vector node_y(curvilinearGrid.num_m * curvilinearGrid.num_n); curvilinearGrid.node_x = node_x.data(); curvilinearGrid.node_y = node_y.data(); - errorCode = meshkernelapi::mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); + errorCode = mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); // Make copy of node values. @@ -1280,14 +1280,14 @@ TEST(CurvilinearGrid, SnapToLandBoundary) double sectionControlPoint2x = 90.0; double sectionControlPoint2y = 100.0; - errorCode = meshkernelapi::mkernel_curvilinear_snap_to_landboundary(meshKernelId, land, - sectionControlPoint1x, - sectionControlPoint1y, - sectionControlPoint2x, - sectionControlPoint2y); + errorCode = mkernel_curvilinear_snap_to_landboundary(meshKernelId, land, + sectionControlPoint1x, + sectionControlPoint1y, + sectionControlPoint2x, + sectionControlPoint2y); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); - errorCode = meshkernelapi::mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); + errorCode = mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); for (int i = 0; i < curvilinearGrid.num_m * curvilinearGrid.num_n; ++i) @@ -1305,7 +1305,7 @@ TEST(CurvilinearGrid, SnapToLandBoundary) ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); EXPECT_TRUE(didUndoOfDeleteNode); - errorCode = meshkernelapi::mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); + errorCode = mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); for (int i = 0; i < curvilinearGrid.num_m * curvilinearGrid.num_n; ++i) @@ -1324,7 +1324,7 @@ TEST(CurvilinearGrid, SnapToSpline) auto errorCode = meshkernelapi::mkernel_allocate_state(0, meshKernelId); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); - meshkernel::MakeGridParameters makeGridParameters; + MakeGridParameters makeGridParameters; makeGridParameters.num_columns = 10; makeGridParameters.num_rows = 10; @@ -1338,14 +1338,14 @@ TEST(CurvilinearGrid, SnapToSpline) ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); meshkernelapi::CurvilinearGrid curvilinearGrid{}; - errorCode = meshkernelapi::mkernel_curvilinear_get_dimensions(meshKernelId, curvilinearGrid); + errorCode = mkernel_curvilinear_get_dimensions(meshKernelId, curvilinearGrid); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); std::vector node_x(curvilinearGrid.num_m * curvilinearGrid.num_n); std::vector node_y(curvilinearGrid.num_m * curvilinearGrid.num_n); curvilinearGrid.node_x = node_x.data(); curvilinearGrid.node_y = node_y.data(); - errorCode = meshkernelapi::mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); + errorCode = mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); // Make copy of node values. @@ -1399,14 +1399,14 @@ TEST(CurvilinearGrid, SnapToSpline) double sectionControlPoint2x = 10.0; double sectionControlPoint2y = 1.0; - errorCode = meshkernelapi::mkernel_curvilinear_snap_to_spline(meshKernelId, spline, - sectionControlPoint1x, - sectionControlPoint1y, - sectionControlPoint2x, - sectionControlPoint2y); + errorCode = mkernel_curvilinear_snap_to_spline(meshKernelId, spline, + sectionControlPoint1x, + sectionControlPoint1y, + sectionControlPoint2x, + sectionControlPoint2y); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); - errorCode = meshkernelapi::mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); + errorCode = mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); for (int i = 0; i < curvilinearGrid.num_m * curvilinearGrid.num_n; ++i) @@ -1424,7 +1424,7 @@ TEST(CurvilinearGrid, SnapToSpline) ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); EXPECT_TRUE(didUndoOfDeleteNode); - errorCode = meshkernelapi::mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); + errorCode = mkernel_curvilinear_get_data(meshKernelId, curvilinearGrid); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); for (int i = 0; i < curvilinearGrid.num_m * curvilinearGrid.num_n; ++i) @@ -1433,3 +1433,62 @@ TEST(CurvilinearGrid, SnapToSpline) EXPECT_NEAR(curvilinearGrid.node_y[i], originalNodeY[i], tolerance); } } + +class CurvilineartBoundariesAsPolygonsTests : public testing::TestWithParam, std::vector>> +{ +public: + [[nodiscard]] static std::vector, std::vector>> GetData() + { + return { + std::make_tuple, std::vector>(0, 0, 9, 9, std::vector{0, 10, 20, 30}, std::vector{0, 0, 0, 0}), + std::make_tuple, std::vector>(1, 1, 8, 8, std::vector{10, 20, 30, 40}, std::vector{10, 10, 10, 10})}; + } +}; + +TEST_P(CurvilineartBoundariesAsPolygonsTests, GetLocationIndex_OnACurvilinearGrid_ShouldGetTheLocationIndex) +{ + // Prepare + auto const& [lowerLeftN, lowerLeftM, upperRightN, upperRightM, ValidCoordinatesX, ValidCoordinatesY] = GetParam(); + + // Prepare + int meshKernelId; + auto errorCode = meshkernelapi::mkernel_allocate_state(0, meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + MakeGridParameters makeGridParameters; + + makeGridParameters.num_columns = 10; + makeGridParameters.num_rows = 10; + makeGridParameters.angle = 0.0; + makeGridParameters.origin_x = 0.0; + makeGridParameters.origin_y = 0.0; + makeGridParameters.block_size_x = 10.0; + makeGridParameters.block_size_y = 10.0; + + errorCode = meshkernelapi::mkernel_curvilinear_compute_rectangular_grid(meshKernelId, makeGridParameters); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + int numberOfPolygonNodes; + errorCode = meshkernelapi::mkernel_curvilinear_count_boundaries_as_polygons(meshKernelId, lowerLeftN, lowerLeftM, upperRightN, upperRightM, numberOfPolygonNodes); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + meshkernelapi::GeometryList boundaryPolygon; + std::vector coordinates_x(numberOfPolygonNodes); + std::vector coordinates_y(numberOfPolygonNodes); + + boundaryPolygon.coordinates_x = coordinates_x.data(); + boundaryPolygon.coordinates_y = coordinates_y.data(); + boundaryPolygon.num_coordinates = numberOfPolygonNodes; + boundaryPolygon.geometry_separator = constants::missing::doubleValue; + + errorCode = mkernel_curvilinear_get_boundaries_as_polygons(meshKernelId, lowerLeftN, lowerLeftM, upperRightN, upperRightM, boundaryPolygon); + + std::vector firstFourCoordinatesX(coordinates_x.begin(), coordinates_x.begin() + 4); + ASSERT_THAT(firstFourCoordinatesX, ::testing::ContainerEq(ValidCoordinatesX)); + + std::vector firstFourCoordinatesY(coordinates_y.begin(), coordinates_y.begin() + 4); + ASSERT_THAT(firstFourCoordinatesY, ::testing::ContainerEq(ValidCoordinatesY)); + + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); +} +INSTANTIATE_TEST_SUITE_P(CurvilineartBoundariesAsPolygonsTests, CurvilineartBoundariesAsPolygonsTests, ::testing::ValuesIn(CurvilineartBoundariesAsPolygonsTests::GetData()));