From 998dfd0112ad88ea36c7c7f0c03d11f02f14e874 Mon Sep 17 00:00:00 2001 From: Ahmad El Sayed Date: Tue, 18 Jul 2023 08:46:53 +0200 Subject: [PATCH] More replacemenmts in CurvilinearGrid --- .../CurvilinearGrid/CurvilinearGrid.hpp | 2 +- .../MeshKernel/Utilities/LinearAlgebra.hpp | 4 +- .../src/CurvilinearGrid/CurvilinearGrid.cpp | 65 ++++++++++--------- .../CurvilinearGridOrthogonalization.cpp | 32 ++++----- .../CurvilinearGridSmoothing.cpp | 36 +++++----- 5 files changed, 69 insertions(+), 70 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp index 74edfc5ab..4f9ab52e2 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp @@ -183,7 +183,7 @@ namespace meshkernel UInt m_numN = 0; ///< The number of n coordinates (horizontal lines) std::vector> m_gridNodes; ///< Member variable storing the grid lin_alg::MatrixRowMajor m_gridFacesMask; ///< The mask of the grid faces (true/false) - std::vector> m_gridNodesTypes; ///< The grid node types + lin_alg::MatrixRowMajor m_gridNodesTypes; ///< The grid node types std::vector m_gridIndices; ///< The original mapping of the flatten nodes in the curvilinear grid private: diff --git a/libs/MeshKernel/include/MeshKernel/Utilities/LinearAlgebra.hpp b/libs/MeshKernel/include/MeshKernel/Utilities/LinearAlgebra.hpp index ae263ceec..cb4bbc11b 100644 --- a/libs/MeshKernel/include/MeshKernel/Utilities/LinearAlgebra.hpp +++ b/libs/MeshKernel/include/MeshKernel/Utilities/LinearAlgebra.hpp @@ -2,8 +2,6 @@ #include "MeshKernel/Definitions.hpp" -// Not allowed, must be signed! -// #define EIGEN_DEFAULT_DENSE_INDEX_TYPE UInt #include namespace meshkernel @@ -29,7 +27,7 @@ namespace meshkernel /// @param[in] preserve Optional parameter that determines if the existing elements in /// the matrix should be presrved or overwritten by the provided /// or deflault (if not provided) fill value - /// @param[in] fill_value The otional the fill value (defaults to the template class constructor) + /// @param[in] fill_value The optional the fill value (defaults to the template class constructor) template inline static void ResizeAndFillMatrix(MatrixRowMajor& matrix, UInt rows, diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp index 91dd57ddf..3e529689a 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp @@ -27,6 +27,7 @@ #include #include +#include #include #include @@ -358,7 +359,7 @@ namespace meshkernel void CurvilinearGrid::ComputeGridNodeTypes() { RemoveInvalidNodes(true); - ResizeAndFill2DVector(m_gridNodesTypes, m_numM, m_numN, true, NodeType::Invalid); + lin_alg::ResizeAndFillMatrix(m_gridNodesTypes, m_numM, m_numN, true, NodeType::Invalid); // Flag faces based on boundaries for (UInt m = 0; m < m_numM; ++m) @@ -374,85 +375,85 @@ namespace meshkernel // Left side if (m == 0 && n == 0) { - m_gridNodesTypes[m][n] = NodeType::BottomLeft; + m_gridNodesTypes(m, n) = NodeType::BottomLeft; continue; } if (m == 0 && n == m_numN - 1) { - m_gridNodesTypes[m][n] = NodeType::UpperLeft; + m_gridNodesTypes(m, n) = NodeType::UpperLeft; continue; } if (m == 0 && !m_gridNodes[m][n - 1].IsValid()) { - m_gridNodesTypes[m][n] = NodeType::BottomLeft; + m_gridNodesTypes(m, n) = NodeType::BottomLeft; continue; } if (m == 0 && !m_gridNodes[m][n + 1].IsValid()) { - m_gridNodesTypes[m][n] = NodeType::UpperLeft; + m_gridNodesTypes(m, n) = NodeType::UpperLeft; continue; } if (m == 0) { - m_gridNodesTypes[m][n] = NodeType::Left; + m_gridNodesTypes(m, n) = NodeType::Left; continue; } // Right side if (m == m_numM - 1 && n == 0) { - m_gridNodesTypes[m][n] = NodeType::BottomRight; + m_gridNodesTypes(m, n) = NodeType::BottomRight; continue; } if (m == m_numM - 1 && n == m_numN - 1) { - m_gridNodesTypes[m][n] = NodeType::UpperRight; + m_gridNodesTypes(m, n) = NodeType::UpperRight; continue; } if (m == m_numM - 1 && !m_gridNodes[m][n - 1].IsValid()) { - m_gridNodesTypes[m][n] = NodeType::BottomRight; + m_gridNodesTypes(m, n) = NodeType::BottomRight; continue; } if (m == m_numM - 1 && !m_gridNodes[m][n + 1].IsValid()) { - m_gridNodesTypes[m][n] = NodeType::UpperRight; + m_gridNodesTypes(m, n) = NodeType::UpperRight; continue; } if (m == m_numM - 1) { - m_gridNodesTypes[m][n] = NodeType::Right; + m_gridNodesTypes(m, n) = NodeType::Right; continue; } // Bottom side if (n == 0 && !m_gridNodes[m - 1][n].IsValid()) { - m_gridNodesTypes[m][n] = NodeType::BottomLeft; + m_gridNodesTypes(m, n) = NodeType::BottomLeft; continue; } if (n == 0 && !m_gridNodes[m + 1][n].IsValid()) { - m_gridNodesTypes[m][n] = NodeType::BottomRight; + m_gridNodesTypes(m, n) = NodeType::BottomRight; continue; } if (n == 0) { - m_gridNodesTypes[m][n] = NodeType::Bottom; + m_gridNodesTypes(m, n) = NodeType::Bottom; continue; } // Upper side if (n == m_numN - 1 && !m_gridNodes[m - 1][n].IsValid()) { - m_gridNodesTypes[m][n] = NodeType::UpperLeft; + m_gridNodesTypes(m, n) = NodeType::UpperLeft; continue; } if (n == m_numN - 1 && !m_gridNodes[m + 1][n].IsValid()) { - m_gridNodesTypes[m][n] = NodeType::UpperRight; + m_gridNodesTypes(m, n) = NodeType::UpperRight; continue; } if (n == m_numN - 1) { - m_gridNodesTypes[m][n] = NodeType::Up; + m_gridNodesTypes(m, n) = NodeType::Up; continue; } @@ -466,7 +467,7 @@ namespace meshkernel isBottomLeftFaceValid && isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::InternalValid; + m_gridNodesTypes(m, n) = NodeType::InternalValid; continue; } if (!isTopRightFaceValid && @@ -474,7 +475,7 @@ namespace meshkernel isBottomLeftFaceValid && isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::BottomLeft; + m_gridNodesTypes(m, n) = NodeType::BottomLeft; continue; } if (isTopRightFaceValid && @@ -482,7 +483,7 @@ namespace meshkernel isBottomLeftFaceValid && isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::BottomRight; + m_gridNodesTypes(m, n) = NodeType::BottomRight; continue; } if (isTopRightFaceValid && @@ -490,7 +491,7 @@ namespace meshkernel !isBottomLeftFaceValid && isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::UpperRight; + m_gridNodesTypes(m, n) = NodeType::UpperRight; continue; } if (isTopRightFaceValid && @@ -498,7 +499,7 @@ namespace meshkernel isBottomLeftFaceValid && !isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::UpperLeft; + m_gridNodesTypes(m, n) = NodeType::UpperLeft; continue; } @@ -507,7 +508,7 @@ namespace meshkernel !isBottomLeftFaceValid && !isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::Bottom; + m_gridNodesTypes(m, n) = NodeType::Bottom; continue; } if (isTopRightFaceValid && @@ -515,7 +516,7 @@ namespace meshkernel !isBottomLeftFaceValid && isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::Left; + m_gridNodesTypes(m, n) = NodeType::Left; continue; } @@ -524,7 +525,7 @@ namespace meshkernel isBottomLeftFaceValid && isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::Up; + m_gridNodesTypes(m, n) = NodeType::Up; continue; } @@ -533,7 +534,7 @@ namespace meshkernel isBottomLeftFaceValid && !isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::Right; + m_gridNodesTypes(m, n) = NodeType::Right; continue; } @@ -542,7 +543,7 @@ namespace meshkernel !isBottomLeftFaceValid && !isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::BottomLeft; + m_gridNodesTypes(m, n) = NodeType::BottomLeft; continue; } @@ -551,7 +552,7 @@ namespace meshkernel !isBottomLeftFaceValid && !isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::BottomRight; + m_gridNodesTypes(m, n) = NodeType::BottomRight; continue; } @@ -560,7 +561,7 @@ namespace meshkernel isBottomLeftFaceValid && !isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::UpperRight; + m_gridNodesTypes(m, n) = NodeType::UpperRight; continue; } @@ -569,7 +570,7 @@ namespace meshkernel !isBottomLeftFaceValid && isBottomRightFaceValid) { - m_gridNodesTypes[m][n] = NodeType::UpperLeft; + m_gridNodesTypes(m, n) = NodeType::UpperLeft; } } } @@ -641,8 +642,8 @@ namespace meshkernel CurvilinearGrid::BoundaryGridLineType CurvilinearGrid::GetBoundaryGridLineType(CurvilinearGridNodeIndices const& firstNode, CurvilinearGridNodeIndices const& secondNode) const { - auto const firstNodeType = m_gridNodesTypes[firstNode.m_m][firstNode.m_n]; - auto const secondNodeType = m_gridNodesTypes[secondNode.m_m][secondNode.m_n]; + auto const firstNodeType = m_gridNodesTypes(firstNode.m_m, firstNode.m_n); + auto const secondNodeType = m_gridNodesTypes(secondNode.m_m, secondNode.m_n); if (firstNodeType == NodeType::InternalValid || firstNodeType == NodeType::Invalid || secondNodeType == NodeType::InternalValid || secondNodeType == NodeType::Invalid) diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp index 3fbbb6261..4bb5a5452 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp @@ -106,7 +106,7 @@ void CurvilinearGridOrthogonalization::ProjectHorizontalBoundaryGridNodes() int nextVertical = 0; for (UInt m = 0; m < m_grid.m_numM; ++m) { - const auto nodeType = m_grid.m_gridNodesTypes[m][n]; + const auto nodeType = m_grid.m_gridNodesTypes(m, n); if (nodeType == CurvilinearGrid::NodeType::BottomLeft || nodeType == CurvilinearGrid::NodeType::UpperLeft) { startM = m; @@ -136,7 +136,7 @@ void CurvilinearGridOrthogonalization::ProjectHorizontalBoundaryGridNodes() { continue; } - if (m_grid.m_gridNodesTypes[mm][n] == CurvilinearGrid::NodeType::Invalid) + if (m_grid.m_gridNodesTypes(mm, n) == CurvilinearGrid::NodeType::Invalid) { continue; } @@ -185,7 +185,7 @@ void CurvilinearGridOrthogonalization::ProjectVerticalBoundariesGridNodes() int nextHorizontal = 0; for (UInt n = 0; n < m_grid.m_numN; ++n) { - const auto nodeType = m_grid.m_gridNodesTypes[m][n]; + const auto nodeType = m_grid.m_gridNodesTypes(m, n); if (nodeType == CurvilinearGrid::NodeType::BottomLeft || nodeType == CurvilinearGrid::NodeType::BottomRight) { startN = n; @@ -215,7 +215,7 @@ void CurvilinearGridOrthogonalization::ProjectVerticalBoundariesGridNodes() { continue; } - if (m_grid.m_gridNodesTypes[m][nn] == CurvilinearGrid::NodeType::Invalid) + if (m_grid.m_gridNodesTypes(m, nn) == CurvilinearGrid::NodeType::Invalid) { continue; } @@ -275,7 +275,7 @@ void CurvilinearGridOrthogonalization::Solve() { for (auto n = minNInternal; n < maxNInternal; ++n) { - if (m_grid.m_gridNodesTypes[m][n] != CurvilinearGrid::NodeType::InternalValid) + if (m_grid.m_gridNodesTypes(m, n) != CurvilinearGrid::NodeType::InternalValid) { continue; } @@ -390,7 +390,7 @@ void CurvilinearGridOrthogonalization::ComputeCoefficients() { for (auto n = m_lowerLeft.m_n + 1; n < m_upperRight.m_n; ++n) { - if (m_grid.m_gridNodesTypes[m][n] != CurvilinearGrid::NodeType::InternalValid) + if (m_grid.m_gridNodesTypes(m, n) != CurvilinearGrid::NodeType::InternalValid) { continue; } @@ -520,19 +520,19 @@ CurvilinearGridOrthogonalization::ComputeInvalidHorizontalBoundaryNodes() const for (auto n = m_lowerLeft.m_n + 1; n < m_upperRight.m_n; ++n) { int step = 0; - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::BottomLeft) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::BottomLeft) { step = -1; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::BottomRight) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::BottomRight) { step = 1; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::UpperRight) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::UpperRight) { step = 1; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::UpperLeft) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::UpperLeft) { step = -1; } @@ -544,7 +544,7 @@ CurvilinearGridOrthogonalization::ComputeInvalidHorizontalBoundaryNodes() const auto lastValidM = m + step; while (lastValidM > 0 && lastValidM < m_grid.m_numM && - m_grid.m_gridNodesTypes[lastValidM][n] == CurvilinearGrid::NodeType::InternalValid) + m_grid.m_gridNodesTypes(lastValidM, n) == CurvilinearGrid::NodeType::InternalValid) { lastValidM += step; } @@ -569,19 +569,19 @@ CurvilinearGridOrthogonalization::ComputeInvalidVerticalBoundaryNodes() const for (auto n = m_lowerLeft.m_n + 1; n < m_upperRight.m_n; ++n) { int step = 0; - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::BottomLeft) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::BottomLeft) { step = -1; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::BottomRight) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::BottomRight) { step = -1; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::UpperRight) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::UpperRight) { step = 1; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::UpperLeft) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::UpperLeft) { step = 1; } @@ -592,7 +592,7 @@ CurvilinearGridOrthogonalization::ComputeInvalidVerticalBoundaryNodes() const auto lastValidN = n + step; while (lastValidN > 0 && lastValidN < m_grid.m_numN && - m_grid.m_gridNodesTypes[m][lastValidN] == CurvilinearGrid::NodeType::InternalValid) + m_grid.m_gridNodesTypes(m, lastValidN) == CurvilinearGrid::NodeType::InternalValid) { lastValidN += step; } diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSmoothing.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSmoothing.cpp index 55a2fe007..a7992a729 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSmoothing.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSmoothing.cpp @@ -109,14 +109,14 @@ void CurvilinearGridSmoothing::SolveDirectional() { if (m_lines[0].IsMGridLine()) { - return m_grid.m_gridNodesTypes[m][n] != CurvilinearGrid::NodeType::InternalValid && - m_grid.m_gridNodesTypes[m][n] != CurvilinearGrid::NodeType::Bottom && - m_grid.m_gridNodesTypes[m][n] != CurvilinearGrid::NodeType::Up; + return m_grid.m_gridNodesTypes(m, n) != CurvilinearGrid::NodeType::InternalValid && + m_grid.m_gridNodesTypes(m, n) != CurvilinearGrid::NodeType::Bottom && + m_grid.m_gridNodesTypes(m, n) != CurvilinearGrid::NodeType::Up; } - return m_grid.m_gridNodesTypes[m][n] != CurvilinearGrid::NodeType::InternalValid && - m_grid.m_gridNodesTypes[m][n] != CurvilinearGrid::NodeType::Left && - m_grid.m_gridNodesTypes[m][n] != CurvilinearGrid::NodeType::Right; + return m_grid.m_gridNodesTypes(m, n) != CurvilinearGrid::NodeType::InternalValid && + m_grid.m_gridNodesTypes(m, n) != CurvilinearGrid::NodeType::Left && + m_grid.m_gridNodesTypes(m, n) != CurvilinearGrid::NodeType::Right; }; // Apply smoothing @@ -190,17 +190,17 @@ void CurvilinearGridSmoothing::Solve() { // It is invalid or a corner point, skip smoothing - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::Invalid || - m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::BottomLeft || - m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::UpperLeft || - m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::BottomRight || - m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::UpperRight) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::Invalid || + m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::BottomLeft || + m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::UpperLeft || + m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::BottomRight || + m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::UpperRight) { continue; } // Compute new position based on a smoothing operator - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::InternalValid) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::InternalValid) { m_grid.m_gridNodes[m][n] = m_gridNodesCache[m][n] * a + (m_gridNodesCache[m - 1][n] + m_gridNodesCache[m + 1][n]) * 0.25 * b + (m_gridNodesCache[m][n - 1] + m_gridNodesCache[m][n + 1]) * 0.25 * b; @@ -209,19 +209,19 @@ void CurvilinearGridSmoothing::Solve() // For the point on the boundaries first computed the new position Point newNodePosition; - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::Bottom) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::Bottom) { newNodePosition = m_gridNodesCache[m][n] * a + (m_gridNodesCache[m - 1][n] + m_gridNodesCache[m + 1][n] + m_gridNodesCache[m][n + 1]) * constants::numeric::oneThird * b; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::Up) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::Up) { newNodePosition = m_gridNodesCache[m][n] * a + (m_gridNodesCache[m - 1][n] + m_gridNodesCache[m + 1][n] + m_gridNodesCache[m][n - 1]) * constants::numeric::oneThird * b; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::Right) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::Right) { newNodePosition = m_gridNodesCache[m][n] * a + (m_gridNodesCache[m][n - 1] + m_gridNodesCache[m][n + 1] + m_gridNodesCache[m - 1][n]) * constants::numeric::oneThird * b; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::Left) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::Left) { newNodePosition = m_gridNodesCache[m][n] * a + (m_gridNodesCache[m][n - 1] + m_gridNodesCache[m][n + 1] + m_gridNodesCache[m + 1][n]) * constants::numeric::oneThird * b; } @@ -236,12 +236,12 @@ void CurvilinearGridSmoothing::ProjectPointOnClosestGridBoundary(Point const& po // Project the new position on the original boundary segment Point previousNode; Point nextNode; - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::Bottom || m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::Up) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::Bottom || m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::Up) { previousNode = m_gridNodesCache[m - 1][n]; nextNode = m_gridNodesCache[m + 1][n]; } - if (m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::Right || m_grid.m_gridNodesTypes[m][n] == CurvilinearGrid::NodeType::Left) + if (m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::Right || m_grid.m_gridNodesTypes(m, n) == CurvilinearGrid::NodeType::Left) { previousNode = m_gridNodesCache[m][n - 1]; nextNode = m_gridNodesCache[m][n + 1];