GRIDEDIT-1663: Fix line mirror after generation of curvilinear grid f…

…rom splines (#421)
Deltares · Feb 20, 2025 · 283a36c · 283a36c
1 parent 001945e
commit 283a36c
Show file tree

Hide file tree

Showing 4 changed files with 184 additions and 7 deletions.
diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp
@@ -548,6 +548,9 @@ namespace meshkernel
         /// @brief Get the node type of the interior nodes
         NodeType GetInteriorNodeType(const UInt n, const UInt m) const;
 
+        /// @brief Computes the indices of the first valid row and columns of a grid node matrix
+        std::tuple<UInt, UInt, UInt, UInt> TrimGridNodes(const lin_alg::Matrix<Point>& gridNodes) const;
+
         Projection m_projection;                               ///< The curvilinear grid projection
         lin_alg::Matrix<Point> m_gridNodes;                    ///< Member variable storing the grid
         lin_alg::Matrix<bool> m_gridFacesMask;                 ///< The mask of the grid faces (true/false)

diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp
@@ -136,12 +136,28 @@ CurvilinearGrid& CurvilinearGrid::operator=(const CurvilinearGrid& copy)
 
 void CurvilinearGrid::SetGridNodes(const lin_alg::Matrix<Point>& gridNodes)
 {
-    if (gridNodes.rows() <= 1 || gridNodes.cols() <= 1)
+    const auto [firstValidRow, lastValidRow, firstValidCol, lastValidCol] = TrimGridNodes(gridNodes);
+
+    if (lastValidRow < firstValidRow || lastValidCol < firstValidCol)
     {
-        throw std::invalid_argument("CurvilinearGrid::CurvilinearGrid: Invalid curvilinear grid nodes");
+        throw std::invalid_argument("CurvilinearGrid::SetGridNodes: Invalid curvilinear grid nodes");
     }
 
-    m_gridNodes = gridNodes;
+    // If the entire grid is valid, move it directly
+    if (firstValidRow == 0 && lastValidRow == gridNodes.rows() - 1 &&
+        firstValidCol == 0 && lastValidCol == gridNodes.cols() - 1)
+    {
+        m_gridNodes = gridNodes;
+    }
+    else
+    {
+        UInt newRows = lastValidRow - firstValidRow + 1;
+        UInt newCols = lastValidCol - firstValidCol + 1;
+        m_gridNodes = gridNodes.block(firstValidRow,
+                                      firstValidCol,
+                                      newRows,
+                                      newCols);
+    }
 
     m_nodesRTreeRequiresUpdate = true;
     m_edgesRTreeRequiresUpdate = true;
@@ -152,20 +168,129 @@ void CurvilinearGrid::SetGridNodes(const lin_alg::Matrix<Point>& gridNodes)
 
 void CurvilinearGrid::SetGridNodes(lin_alg::Matrix<Point>&& gridNodes)
 {
-    if (gridNodes.rows() <= 1 || gridNodes.cols() <= 1)
+    const auto [firstValidRow, lastValidRow, firstValidCol, lastValidCol] = TrimGridNodes(gridNodes);
+
+    if (lastValidRow < firstValidRow || lastValidCol < firstValidCol)
     {
-        throw std::invalid_argument("CurvilinearGrid::CurvilinearGrid: Invalid curvilinear grid nodes");
+        throw std::invalid_argument("CurvilinearGrid::SetGridNodes: Invalid curvilinear grid nodes");
     }
 
-    m_gridNodes = std::move(gridNodes);
+    // If the entire grid is valid, move it directly
+    if (firstValidRow == 0 && lastValidRow == gridNodes.rows() - 1 &&
+        firstValidCol == 0 && lastValidCol == gridNodes.cols() - 1)
+    {
+        m_gridNodes = std::move(gridNodes);
+    }
+    else
+    {
+        UInt newRows = lastValidRow - firstValidRow + 1;
+        UInt newCols = lastValidCol - firstValidCol + 1;
+        m_gridNodes = gridNodes.block(firstValidRow,
+                                      firstValidCol,
+                                      newRows,
+                                      newCols);
+    }
 
+    // Mark R-Trees as requiring updates
     m_nodesRTreeRequiresUpdate = true;
     m_edgesRTreeRequiresUpdate = true;
     m_facesRTreeRequiresUpdate = true;
 
+    // Compute new indices
     m_gridIndices = ComputeNodeIndices();
 }
 
+std::tuple<meshkernel::UInt,
+           meshkernel::UInt,
+           meshkernel::UInt,
+           meshkernel::UInt>
+CurvilinearGrid::TrimGridNodes(const lin_alg::Matrix<Point>& gridNodes) const
+{
+    const auto rows = static_cast<UInt>(gridNodes.rows());
+    const auto cols = static_cast<UInt>(gridNodes.cols());
+
+    // Initialize valid row/column bounds
+    UInt firstValidRow = 0;
+    UInt lastValidRow = rows - 1;
+    UInt firstValidCol = 0;
+    UInt lastValidCol = cols - 1;
+
+    bool foundFirstRow = false;
+    bool foundFirstCol = false;
+    bool foundLastRow = false;
+    bool foundLastCol = false;
+
+    // Find first valid row
+    for (UInt r = 0; r < rows; ++r)
+    {
+        for (UInt c = 0; c < cols; ++c)
+        {
+            if (gridNodes(r, c).x != constants::missing::doubleValue &&
+                gridNodes(r, c).y != constants::missing::doubleValue)
+            {
+                firstValidRow = r;
+                foundFirstRow = true;
+                break;
+            }
+        }
+        if (foundFirstRow)
+            break;
+    }
+
+    // Find last valid row
+    for (UInt r = rows - 1; r >= firstValidRow; --r)
+    {
+        for (UInt c = 0; c < cols; ++c)
+        {
+            if (gridNodes(r, c).x != constants::missing::doubleValue &&
+                gridNodes(r, c).y != constants::missing::doubleValue)
+            {
+                lastValidRow = r;
+                foundLastRow = true;
+                break;
+            }
+        }
+        if (foundLastRow)
+            break;
+    }
+
+    // Find first valid column
+    for (UInt c = 0; c < cols; ++c)
+    {
+        for (UInt r = firstValidRow; r <= lastValidRow; ++r)
+        {
+            if (gridNodes(r, c).x != constants::missing::doubleValue &&
+                gridNodes(r, c).y != constants::missing::doubleValue)
+            {
+                firstValidCol = c;
+                foundFirstCol = true;
+                break;
+            }
+        }
+        if (foundFirstCol)
+            break;
+    }
+
+    // Find last valid column
+    for (UInt c = cols - 1; c >= firstValidCol; --c)
+    {
+        for (UInt r = firstValidRow; r <= lastValidRow; ++r)
+        {
+            if (gridNodes(r, c).x != constants::missing::doubleValue &&
+                gridNodes(r, c).y != constants::missing::doubleValue)
+            {
+                lastValidCol = c;
+                foundLastCol = true;
+                break;
+            }
+        }
+        if (foundLastCol)
+            break;
+    }
+
+    return {firstValidRow, lastValidRow, firstValidCol, lastValidCol};
+}
+
 void CurvilinearGrid::Delete(std::shared_ptr<Polygons> polygons, UInt polygonIndex)
 {
     // no polygons available

diff --git a/libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridBasicTests.cpp
@@ -971,3 +971,52 @@ TEST(CurvilinearBasicTests, InsertMultipleFacesAlongAllBoundaryEdges)
 
     //--------------------------------
 }
+
+class CurvilinearSetGridTests : public ::testing::TestWithParam<bool>
+{
+};
+
+TEST_P(CurvilinearSetGridTests, SetGridNodes_ShouldTrimMatrix)
+{
+    // SetUp
+    meshkernel::CurvilinearGrid curvilinearGrid;
+    lin_alg::Matrix<meshkernel::Point> gridNodes(3, 3);
+
+    // Fill the matrix with valid Point values
+    for (int r = 0; r < 3; ++r)
+    {
+        for (int c = 0; c < 3; ++c)
+        {
+            gridNodes(r, c) = meshkernel::Point{static_cast<double>(r), static_cast<double>(c)};
+        }
+    }
+    for (int r = 0; r < 3; ++r)
+    {
+        gridNodes(r, 0) = meshkernel::Point{meshkernel::constants::missing::doubleValue,
+                                            meshkernel::constants::missing::doubleValue};
+        gridNodes(r, 2) = meshkernel::Point{meshkernel::constants::missing::doubleValue,
+                                            meshkernel::constants::missing::doubleValue};
+    }
+
+    // Execute
+    if (GetParam())
+    {
+        // lvalue
+        curvilinearGrid.SetGridNodes(gridNodes);
+    }
+    else
+    {
+        // rvalue
+        curvilinearGrid.SetGridNodes(std::move(gridNodes));
+    }
+
+    // Assert: Ensure only one column remains
+    EXPECT_EQ(curvilinearGrid.NumN(), 3);
+    EXPECT_EQ(curvilinearGrid.NumM(), 1);
+}
+
+// Run the test for both lvalue (true) and rvalue (false) cases
+INSTANTIATE_TEST_SUITE_P(
+    CurvilinearGrid,
+    CurvilinearSetGridTests,
+    ::testing::Values(true, false));
diff --git a/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp
@@ -1285,7 +1285,7 @@ TEST(CurvilinearGridFromSplines, GridFromSeventySplines)
     }
 
     EXPECT_EQ(grid.NumN(), 36);
-    EXPECT_EQ(grid.NumM(), 306);
+    EXPECT_EQ(grid.NumM(), 301);
     EXPECT_EQ(validPointCount, 4941);
 }