Merge branch 'master' into feature/GRIDEDIT-1244_split_row_column

libs/MeshKernel/include/MeshKernel/Mesh2DToCurvilinear.hpp

@@ -52,6 +52,129 @@ namespace meshkernel
         std::unique_ptr<CurvilinearGrid> Compute(const Point& point);
 
     private:
+        /// @brief Matrix class supporting negative indices, mimicking behavior similar to Fortran.
+        /// This class allows matrix access with negative indices, which is not possible in C++ by default.
+        class MatrixWithNegativeIndices
+        {
+        public:
+            /// @brief Retrieves the value at the given position (j, i) in the matrix.
+            /// @param j Row index, can be negative.
+            /// @param i Column index, can be negative.
+            /// @return Value at the specified position.
+            int getValue(int j, int i) const
+            {
+                const auto jIndex = j - m_minJ;
+                const auto iIndex = i - m_minI;
+                const auto value = m_matrix(jIndex, iIndex);
+                return value;
+            }
+            /// @brief Sets the value at the given position (j, i) in the matrix.
+            /// @param j Row index, can be negative.
+            /// @param i Column index, can be negative.
+            /// @param value Value to be set at the specified position.
+            void setValue(int j, int i, int value)
+            {
+                const auto jIndex = j - m_minJ;
+                const auto iIndex = i - m_minI;
+                m_matrix(jIndex, iIndex) = value;
+            }
+
+            /// @brief Checks if the position (j, i) in the matrix is valid (not equal to missing value).
+            /// @param j Row index, can be negative.
+            /// @param i Column index, can be negative.
+            /// @return True if the value at the position is not missing, otherwise false.
+            bool IsValid(int j, int i) const
+            {
+                const auto jIndex = j - m_minJ;
+                const auto iIndex = i - m_minI;
+                return m_matrix(jIndex, iIndex) != missing::intValue;
+            }
+
+            /// @brief Gets the number of rows in the matrix.
+            /// @return Number of rows.
+            [[nodiscard]] int rows() const
+            {
+                return static_cast<int>(m_matrix.rows());
+            }
+
+            /// @brief Gets the number of columns in the matrix.
+            /// @return Number of columns.
+            [[nodiscard]] int cols() const
+            {
+                return static_cast<int>(m_matrix.cols());
+            }
+
+            /// @brief Gets the minimum column index (i) allowed in the matrix.
+            /// @return Minimum column index.
+            [[nodiscard]] int minCol() const
+            {
+                return m_minI;
+            }
+
+            /// @brief Gets the maximum column index (i) allowed in the matrix.
+            /// @return Maximum column index.
+            [[nodiscard]] int maxCol() const
+            {
+                return m_maxI;
+            }
+
+            /// @brief Gets the minimum row index (j) allowed in the matrix.
+            /// @return Minimum row index.
+            [[nodiscard]] int minRow() const
+            {
+                return m_minJ;
+            }
+
+            /// @brief Gets the maximum row index (j) allowed in the matrix.
+            /// @return Maximum row index.
+            [[nodiscard]] int maxRow() const
+            {
+                return m_maxJ;
+            }
+
+            /// @brief Resizes the matrix to accommodate the new bounds.
+            /// @param minJ New minimum row index.
+            /// @param minI New minimum column index.
+            /// @param maxJ New maximum row index.
+            /// @param maxI New maximum column index.
+            void resize(int minJ, int minI, int maxJ, int maxI)
+            {
+                // Determine the size change needed
+                const int extraRowsTop = std::max(m_minJ - minJ, 0);
+                const int extraRowsBottom = std::max(maxJ - m_maxJ, 0);
+                const int extraColsLeft = std::max(m_minI - minI, 0);
+                const int extraColsRight = std::max(maxI - m_maxI, 0);
+
+                if (extraRowsTop == 0 && extraRowsBottom == 0 && extraColsLeft == 0 && extraColsRight == 0)
+                {
+                    return;
+                }
+
+                // Create new matrix with the new size and initialize to missing value
+                const auto newRows = static_cast<int>(m_matrix.rows()) + extraRowsTop + extraRowsBottom;
+                const auto newCols = static_cast<int>(m_matrix.cols()) + extraColsLeft + extraColsRight;
+                lin_alg::Matrix<int> newMatrix(newRows, newCols);
+                newMatrix.setConstant(missing::intValue);
+
+                // Copy the existing matrix into the new one
+                newMatrix.block(extraRowsTop, extraColsLeft, m_matrix.rows(), m_matrix.cols()) = m_matrix;
+
+                // Update matrix and bounds
+                m_matrix.swap(newMatrix);
+                m_minI = std::min(m_minI, minI);
+                m_minJ = std::min(m_minJ, minJ);
+                m_maxI = std::max(m_maxI, maxI);
+                m_maxJ = std::max(m_maxJ, maxJ);
+            }
+
+        private:
+            lin_alg::Matrix<int> m_matrix = lin_alg::Matrix<int>(1, 1); ///< Underlying matrix storage.
+            int m_minI = 0;                                             ///< Minimum column index.
+            int m_minJ = 0;                                             ///< Minimum row index.
+            int m_maxI = 0;                                             ///< Maximum column index.
+            int m_maxJ = 0;                                             ///< Maximum row index.
+        };
+
         /// @brief Computes the local mapping of the nodes composing the face
         [[nodiscard]] Eigen::Matrix<UInt, 2, 2> ComputeLocalNodeMapping(UInt face) const;
 
@@ -61,6 +184,12 @@ namespace meshkernel
                                                         const UInt d,
                                                         const std::vector<bool>& visitedFace);
 
+        /// @brief Checks the valid node and the candidate node are connected and are part of a quadrangular face
+        [[nodiscard]] bool CheckGridLine(const UInt validNode, const UInt candidateNode) const;
+
+        /// @brief Computes the final curvilinear matrix
+        bool IsConnectionValid(const UInt candidateNode, const int iCandidate, const int jCandidate);
+
         /// @brief Computes the final curvilinear matrix
         [[nodiscard]] lin_alg::Matrix<Point> ComputeCurvilinearMatrix();
 
@@ -84,6 +213,8 @@ namespace meshkernel
                                                                        {0, 1}}}; ///< increments for the new nodes depending on the node direction
 
         const int n_maxNumRowsColumns = 1000000; ///< The maximum number of allowed rows or columns
+
+        MatrixWithNegativeIndices m_mapping; ///< Unstructured node indices in the curvilinear grid
     };
 
 } // namespace meshkernel

libs/MeshKernel/src/Mesh2DToCurvilinear.cpp

@@ -83,28 +83,35 @@ std::unique_ptr<CurvilinearGrid> Mesh2DToCurvilinear::Compute(const Point& point
     const auto thirdEdge = m_mesh.m_facesEdges[initialFaceIndex][2];
     const auto fourthEdge = m_mesh.m_facesEdges[initialFaceIndex][3];
 
+    m_mapping.resize(-3, -3, 3, 3);
+
     const auto firstNodeIndex = m_mesh.FindCommonNode(firstEdge, secondEdge);
-    m_i[firstNodeIndex] = 0;
     m_j[firstNodeIndex] = 0;
+    m_i[firstNodeIndex] = 0;
+
+    m_mapping.setValue(0, 0, firstNodeIndex);
 
     const auto secondNodeIndex = m_mesh.FindCommonNode(secondEdge, thirdEdge);
-    m_i[secondNodeIndex] = 1;
     m_j[secondNodeIndex] = 0;
+    m_i[secondNodeIndex] = 1;
+    m_mapping.setValue(0, 1, secondNodeIndex);
 
     const auto thirdNodeIndex = m_mesh.FindCommonNode(thirdEdge, fourthEdge);
     m_i[thirdNodeIndex] = 1;
     m_j[thirdNodeIndex] = 1;
+    m_mapping.setValue(1, 1, thirdNodeIndex);
 
     const auto fourthNodeIndex = m_mesh.FindCommonNode(fourthEdge, firstEdge);
-    m_i[fourthNodeIndex] = 0;
     m_j[fourthNodeIndex] = 1;
+    m_i[fourthNodeIndex] = 0;
+
+    m_mapping.setValue(1, 0, fourthNodeIndex);
 
     // 4. Grow the front using the breath first search algorithm
     const auto numFaces = m_mesh.GetNumFaces();
     std::vector visitedFace(numFaces, false);
     std::queue<UInt> q;
     q.push(initialFaceIndex);
-
     while (!q.empty())
     {
         const auto face = q.front();
@@ -122,7 +129,7 @@ std::unique_ptr<CurvilinearGrid> Mesh2DToCurvilinear::Compute(const Point& point
         }
 
         const auto localNodeMapping = ComputeLocalNodeMapping(face);
-        for (UInt d = 0; d < geometric::numNodesInQuadrilateral; ++d)
+        for (auto d = 0u; d < geometric::numNodesInQuadrilateral; ++d)
         {
             const auto newFaceIndex = ComputeNeighbouringFaceNodes(face, localNodeMapping, d, visitedFace);
             if (newFaceIndex != missing::uintValue)
@@ -217,15 +224,15 @@ UInt Mesh2DToCurvilinear::ComputeNeighbouringFaceNodes(const UInt face,
         }
     }
     auto nextEdgeIndexInNewFace = edgeIndexInNewFace + 1;
-    nextEdgeIndexInNewFace = nextEdgeIndexInNewFace == 4 ? 0 : nextEdgeIndexInNewFace;
+    nextEdgeIndexInNewFace = nextEdgeIndexInNewFace == geometric::numNodesInQuadrilateral ? 0 : nextEdgeIndexInNewFace;
     const auto nextEdgeInNewFace = m_mesh.m_facesEdges[newFace][nextEdgeIndexInNewFace];
     const auto firstCommonNode = m_mesh.FindCommonNode(edgeIndex, nextEdgeInNewFace);
     const auto firstOtherNode = OtherNodeOfEdge(m_mesh.GetEdge(nextEdgeInNewFace), firstCommonNode);
     const auto iFirstOtherNode = m_i[firstCommonNode] + m_directionsDeltas[d][0];
     const auto jFirstOtherNode = m_j[firstCommonNode] + m_directionsDeltas[d][1];
 
     auto previousEdgeIndexInNewFace = edgeIndexInNewFace - 1;
-    previousEdgeIndexInNewFace = previousEdgeIndexInNewFace == -1 ? 3 : previousEdgeIndexInNewFace;
+    previousEdgeIndexInNewFace = previousEdgeIndexInNewFace == -1 ? geometric::numNodesInQuadrilateral - 1 : previousEdgeIndexInNewFace;
     const auto previousEdgeInNewFace = m_mesh.m_facesEdges[newFace][previousEdgeIndexInNewFace];
     const auto secondCommonNode = m_mesh.FindCommonNode(edgeIndex, previousEdgeInNewFace);
     const auto secondOtherNode = OtherNodeOfEdge(m_mesh.GetEdge(previousEdgeInNewFace), secondCommonNode);
@@ -237,60 +244,108 @@ UInt Mesh2DToCurvilinear::ComputeNeighbouringFaceNodes(const UInt face,
                          (m_i[secondOtherNode] != missing::intValue && m_i[secondOtherNode] != iSecondCommonNode) ||
                          (m_j[secondOtherNode] != missing::intValue && m_j[secondOtherNode] != jSecondCommonNode);
 
-    if (!invalid)
+    if (invalid)
     {
-        m_i[firstOtherNode] = iFirstOtherNode;
-        m_j[firstOtherNode] = jFirstOtherNode;
-        m_i[secondOtherNode] = iSecondCommonNode;
-        m_j[secondOtherNode] = jSecondCommonNode;
-        return newFace;
+        return missing::uintValue;
     }
-    return missing::uintValue;
+    if (!IsConnectionValid(firstOtherNode, iFirstOtherNode, jFirstOtherNode))
+    {
+        return missing::uintValue;
+    }
+    if (!IsConnectionValid(secondOtherNode, iSecondCommonNode, jSecondCommonNode))
+    {
+        return missing::uintValue;
+    }
+
+    m_i[firstOtherNode] = iFirstOtherNode;
+    m_j[firstOtherNode] = jFirstOtherNode;
+    m_i[secondOtherNode] = iSecondCommonNode;
+    m_j[secondOtherNode] = jSecondCommonNode;
+    m_mapping.setValue(jFirstOtherNode, iFirstOtherNode, firstOtherNode);
+    m_mapping.setValue(jSecondCommonNode, iSecondCommonNode, secondOtherNode);
+    return newFace;
 }
 
-lin_alg::Matrix<Point> Mesh2DToCurvilinear::ComputeCurvilinearMatrix()
+bool Mesh2DToCurvilinear::IsConnectionValid(const UInt candidateNode, const int iCandidate, const int jCandidate)
 {
+    const int iLeft = iCandidate - 1;
+    const int iRight = iCandidate + 1;
+    const int jBottom = jCandidate - 1;
+    const int jUp = jCandidate + 1;
+    m_mapping.resize(jBottom, iLeft, jUp, iRight);
 
-    int minI = n_maxNumRowsColumns;
-    int minJ = n_maxNumRowsColumns;
-    const auto numNodes = m_mesh.GetNumNodes();
-    for (UInt n = 0; n < numNodes; ++n)
+    if (m_mapping.IsValid(jCandidate, iLeft) && !CheckGridLine(m_mapping.getValue(jCandidate, iLeft), candidateNode))
     {
-        if (m_i[n] != missing::intValue)
-        {
-            minI = std::min(minI, m_i[n]);
-        }
-        if (m_j[n] != missing::intValue)
-        {
-            minJ = std::min(minJ, m_j[n]);
-        }
+        return false;
     }
 
-    int maxI = -n_maxNumRowsColumns;
-    int maxJ = -n_maxNumRowsColumns;
-    for (UInt n = 0; n < numNodes; ++n)
+    if (m_mapping.IsValid(jCandidate, iRight) && !CheckGridLine(m_mapping.getValue(jCandidate, iRight), candidateNode))
     {
-        if (m_i[n] != missing::intValue)
+        return false;
+    }
+
+    if (m_mapping.IsValid(jBottom, iCandidate) && !CheckGridLine(m_mapping.getValue(jBottom, iCandidate), candidateNode))
+    {
+        return false;
+    }
+
+    if (m_mapping.IsValid(jUp, iCandidate) && !CheckGridLine(m_mapping.getValue(jUp, iCandidate), candidateNode))
+    {
+        return false;
+    }
+
+    return true;
+}
+
+bool Mesh2DToCurvilinear::CheckGridLine(const UInt validNode, const UInt candidateNode) const
+{
+    bool valid = false;
+    for (auto e = 0u; e < m_mesh.m_nodesEdges[candidateNode].size(); ++e)
+    {
+        const auto edgeIndex = m_mesh.m_nodesEdges[candidateNode][e];
+        const auto otherNode = OtherNodeOfEdge(m_mesh.GetEdge(edgeIndex), candidateNode);
+
+        bool doCheck = m_mesh.GetNumFaceEdges(m_mesh.m_edgesFaces[edgeIndex][0]) == geometric::numNodesInQuadrilateral;
+
+        if (m_mesh.m_edgesNumFaces[edgeIndex] == 2)
         {
-            m_i[n] = m_i[n] - minI;
-            maxI = std::max(maxI, m_i[n]);
+            doCheck = doCheck || m_mesh.GetNumFaceEdges(m_mesh.m_edgesFaces[edgeIndex][1]) == geometric::numNodesInQuadrilateral;
         }
-        if (m_j[n] != missing::intValue)
+
+        if (doCheck && otherNode == validNode)
         {
-            m_j[n] = m_j[n] - minJ;
-            maxJ = std::max(maxJ, m_j[n]);
+            valid = true;
+            break;
         }
     }
 
-    lin_alg::Matrix<Point> result(maxI + 1, maxJ + 1);
-    for (UInt n = 0; n < numNodes; ++n)
+    return valid;
+}
+
+lin_alg::Matrix<Point> Mesh2DToCurvilinear::ComputeCurvilinearMatrix()
+{
+    auto validI = m_i | std::views::filter([](const auto& x)
+                                           { return x != missing::intValue; });
+    const auto [minI, maxI] = std::ranges::minmax_element(validI);
+
+    auto validJ = m_j | std::views::filter([](const auto& x)
+                                           { return x != missing::intValue; });
+    const auto [minJ, maxJ] = std::ranges::minmax_element(validJ);
+
+    const auto numM = *maxI - *minI + 1;
+    const auto numN = *maxJ - *minJ + 1;
+
+    lin_alg::Matrix<Point> result(numN, numM);
+
+    for (auto n = 0u; n < m_mesh.GetNumNodes(); ++n)
     {
-        const auto i = m_i[n];
-        const auto j = m_j[n];
-        if (i != missing::intValue && j != missing::intValue)
+        if (m_j[n] != missing::intValue && m_i[n] != missing::intValue)
         {
-            result(i, j) = m_mesh.Node(n);
+            const int j = m_j[n] - *minJ;
+            const int i = m_i[n] - *minI;
+            result(j, i) = m_mesh.Node(n);
         }
     }
+
     return result;
 }