From 4f87cbea1d3d5a44a58bc8a794858cb6762930cf Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Wed, 19 Aug 2026 09:18:20 +0200 Subject: [PATCH 1/7] GRIDEDIT-2303 Passing to other computer --- .../CurvilinearGridRectangular.hpp | 21 +++ .../CurvilinearGridRectangular.cpp | 142 ++++++++++++++++++ .../src/CurvilinearGridRectangularTests.cpp | 42 ++++++ 3 files changed, 205 insertions(+) diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp index a3bae2d93..b9aba8054 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp @@ -28,6 +28,7 @@ #pragma once #include +#include #include #include @@ -90,7 +91,19 @@ namespace meshkernel const double upperRightX, const double upperRightY) const; + std::unique_ptr Compute(const double originX, + const double originY, + const double blockSizeX, + const double blockSizeY, + const double upperRightX, + const double upperRightY, + const double angle) const; + private: + + Point RotateUpperRightByAngle (const double originX, const double originY, const double upperRightX, const double upperRightY, const double angle) const; + + /// @brief Compute a rectangular curvilinear grid on cartesian coordinates. /// @param[in] numColumns The number of columns in x direction /// @param[in] numRows The number of columns in y direction @@ -128,6 +141,14 @@ namespace meshkernel const double blockSizeX, const double blockSizeY); + lin_alg::Matrix ComputeSpherical2(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const; + /// @brief Compute the adjusted latitude for keeping an aspect ratio of 1, considering the spherical coordinates /// @param[in] blockSize The grid block size in y dimension /// @param[in] aspectRatio The expected element aspect ratio, defined as being: blockSizeY / blockSizeX diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp index 9ae97c052..fd0f303b2 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp @@ -371,4 +371,146 @@ namespace meshkernel throw NotImplementedError("Projection value {} not supported", static_cast(m_projection)); } + Point CurvilinearGridRectangular::RotateUpperRightByAngle (const double originX, const double originY, const double upperRightX, const double upperRightY, const double angle) const + { + + const double cosA = std::cos(angle * constants::conversion::degToRad); + const double sinA = std::sin(angle * constants::conversion::degToRad); + + if (m_projection == Projection::cartesian) + { + Point translated (upperRightX - originX, upperRightY - originY); + return {cosA * translated.x - sinA * translated.y, sinA * translated.x + cosA * translated.y}; + } + else + { + Cartesian3DPoint rotationPoint = SphericalToCartesian3D ({originX, originY}); + Cartesian3DPoint point3d = SphericalToCartesian3D ({upperRightX, upperRightY}); + + // Normalize the rotation axis (Rodrigues' formula requires a unit vector) + // Points are on Earth's surface, so the magnitude is earth_radius + Cartesian3DPoint k = {rotationPoint.x / constants::geometric::earth_radius, + rotationPoint.y / constants::geometric::earth_radius, + rotationPoint.z / constants::geometric::earth_radius}; + + // std::cout << "rotationPoint " << rotationPoint.x << ", " << rotationPoint.y << ", " << rotationPoint.z << std::endl; + // std::cout << "point3d " << point3d.x << ", " << point3d.y << ", " << point3d.z << std::endl; + // std::cout << "k (unit axis) " << k.x << ", " << k.y << ", " << k.z << std::endl; + + // k · v + double dot = k.x * point3d.x + k.y * point3d.y + k.z * point3d.z; + + // k × v + Cartesian3DPoint crossProd = VectorProduct(k, point3d); + // std::cout << "crossProd " << crossProd.x << ", " << crossProd.y << ", " << crossProd.z << std::endl; + + // Rodrigues' formula: v_rot = v·cos(θ) + (k × v)·sin(θ) + k·(k·v)·(1 - cos(θ)) + Cartesian3DPoint rotatedPoint3d = {point3d.x * cosA + crossProd.x * sinA + k.x * dot * (1.0 - cosA), + point3d.y * cosA + crossProd.y * sinA + k.y * dot * (1.0 - cosA), + point3d.z * cosA + crossProd.z * sinA + k.z * dot * (1.0 - cosA)}; + + // std::cout << "rotatedPoint3d " << rotatedPoint3d.x << ", " << rotatedPoint3d.y << ", " << rotatedPoint3d.z << std::endl; + + Point pointOnSphere = Cartesian3DToSpherical (rotatedPoint3d, upperRightX); + + std::cout << "pointOnSphere "<< pointOnSphere.x << ", " << pointOnSphere.y << std::endl; + + return pointOnSphere; + } + + } + + lin_alg::Matrix CurvilinearGridRectangular::ComputeSpherical2(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const + { + lin_alg::Matrix result = ComputeCartesian(numColumns, + numRows, + originX, + originY, + 0.0 * angle, + blockSizeX, + blockSizeY); + + const auto numM = result.cols(); + const auto numN = result.rows(); + + for (Eigen::Index n = 0; n < numN; ++n) + { + + for (Eigen::Index m = 0; m < numM; ++m) + { + result (n,m) = RotateUpperRightByAngle (originX, originY, result (n,m).x, result (n,m).y, angle); + } + } + + + return result; + } + + + std::unique_ptr CurvilinearGridRectangular::Compute(const double originX, + const double originY, + const double blockSizeX, + const double blockSizeY, + const double upperRightX, + const double upperRightY, + const double angle) const + { + range_check::CheckGreater(blockSizeX, 0.0, "X block size"); + range_check::CheckGreater(blockSizeY, 0.0, "Y block size"); + + Point rotatedUpperRight = RotateUpperRightByAngle (originX, originY, upperRightX, upperRightY, -angle); + + std::cout << "rotated distance: " << rotatedUpperRight.x - originX << ", " << rotatedUpperRight.y - originY << std::endl; + std::cout << "rotated deltas : " << (rotatedUpperRight.x - originX) / blockSizeX << ", " << (rotatedUpperRight.y - originY) / blockSizeY << std::endl; + + + const int numColumns = static_cast(std::ceil((rotatedUpperRight.x - originX) / blockSizeX)); + + if (numColumns <= 0) + { + throw AlgorithmError("Number of columns cannot be <= 0"); + } + + std::cout << "rotatedUpperRight "<< rotatedUpperRight.x << ", " << rotatedUpperRight.y << std::endl; + + const int numRows = static_cast(std::ceil((rotatedUpperRight.y - originY) / blockSizeY)); + // const int numRows = ComputeNumRows(originY, rotatedUpperRight.y, blockSizeX, blockSizeY, m_projection); + std::cout << "num rows, cols "<< numRows << " " << numColumns << " "<< originX << ", " << originY << " "<< rotatedUpperRight.x << ", " << rotatedUpperRight.y << std::endl; + + if (m_projection == Projection::spherical) + { + // double deltax = (rotatedUpperRight.x - originX) / blockSizeX; + // double bsx = + auto grid = std::make_unique(ComputeSpherical(numColumns, + numRows, + originX, + originY, + angle, + blockSizeX, + blockSizeY), + m_projection); + + return grid; + } + if (m_projection == Projection::cartesian) + { + auto grid = std::make_unique(ComputeCartesian(numColumns, + numRows, + originX, + originY, + angle, + blockSizeX, + blockSizeY), + m_projection); + return grid; + } + throw NotImplementedError("Projection value {} not supported", static_cast(m_projection)); + } + } // namespace meshkernel diff --git a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp index b501457ca..d463168ca 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp @@ -694,3 +694,45 @@ TEST_P(CurvilinearGridUniformTests, parameters) ASSERT_EQ(numValidNodes, expectedNumNodes); } INSTANTIATE_TEST_SUITE_P(curvilinearGridDeletionTests, CurvilinearGridUniformTests, ::testing::ValuesIn(CurvilinearGridUniformTests::GetData())); + + +TEST(CurvilinearGridUniform, WithAngle) +{ + + double originX = 10.0; + double originY = 10.0; + double blockSizeX= 1.0; + double blockSizeY = 0.5; + // double upperRightX = 1.283012701892219e+01; + // double upperRightY = 1.509807621135332e+01; + + + double upperRightX = -0.7503; + double upperRightY = 19.2242; + // double upperRightX = 6.8579; + // double upperRightY = 14.8570; + double angle = 90.0; + + // double upperRightX = 12.3963; + // double upperRightY = 14.1524; + // double angle = 30.0; + + + + + // double originX = 0.0; + // double originY = 0.0; + // double blockSizeX= 1.0; + // double blockSizeY = 0.5; + // double upperRightX = 0.0; + // double upperRightY = 5.0 * 1.414213562373095; + // double angle = 45.0; + + + // Execution + CurvilinearGridRectangular curvilinearGridRectangular(Projection::spherical); + auto mesh = curvilinearGridRectangular.Compute (originX, originY, blockSizeX, blockSizeY, upperRightX, upperRightY, angle); + + meshkernel::Print (mesh->ComputeNodes(), mesh->ComputeEdges ()); + +} From 6e84df10b543cf4b19a7c9102f60f1cbd51eaf09 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Wed, 19 Aug 2026 13:06:33 +0200 Subject: [PATCH 2/7] GRIDEDIT-2303 First working version of grid on extension with non zero angle --- .../CurvilinearGridRectangular.hpp | 22 ++--- .../MeshKernel/Utilities/Utilities.hpp | 6 ++ .../CurvilinearGridRectangular.cpp | 98 ++++++++----------- libs/MeshKernel/src/Utilities/Utilities.cpp | 33 +++++++ .../src/CurvilinearGridRectangularTests.cpp | 56 +++++------ .../include/MeshKernelApi/Utils.hpp | 26 +++-- 6 files changed, 137 insertions(+), 104 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp index b9aba8054..54a984136 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp @@ -27,8 +27,8 @@ #pragma once -#include #include +#include #include #include @@ -100,9 +100,9 @@ namespace meshkernel const double angle) const; private: - - Point RotateUpperRightByAngle (const double originX, const double originY, const double upperRightX, const double upperRightY, const double angle) const; - + Point RotateByAngle(const double originX, const double originY, + const double upperRightX, const double upperRightY, + const double cosAngle, const double sinAngle) const; /// @brief Compute a rectangular curvilinear grid on cartesian coordinates. /// @param[in] numColumns The number of columns in x direction @@ -141,13 +141,13 @@ namespace meshkernel const double blockSizeX, const double blockSizeY); - lin_alg::Matrix ComputeSpherical2(const int numColumns, - const int numRows, - const double originX, - const double originY, - const double angle, - const double blockSizeX, - const double blockSizeY) const; + lin_alg::Matrix ComputeSphericalOnExtension(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const; /// @brief Compute the adjusted latitude for keeping an aspect ratio of 1, considering the spherical coordinates /// @param[in] blockSize The grid block size in y dimension diff --git a/libs/MeshKernel/include/MeshKernel/Utilities/Utilities.hpp b/libs/MeshKernel/include/MeshKernel/Utilities/Utilities.hpp index b50fe5525..edffda009 100644 --- a/libs/MeshKernel/include/MeshKernel/Utilities/Utilities.hpp +++ b/libs/MeshKernel/include/MeshKernel/Utilities/Utilities.hpp @@ -31,6 +31,7 @@ #include #include +#include "MeshKernel/Cartesian3DPoint.hpp" #include "MeshKernel/Entities.hpp" #include "MeshKernel/Mesh2D.hpp" #include "MeshKernel/Point.hpp" @@ -43,6 +44,11 @@ namespace meshkernel /// Only nodes and node connectivity need be printed to visualise the graph. void Print(const std::vector& nodes, const std::vector& edges, std::ostream& out = std::cout); + /// @brief Print the (simplified) graph in a form that can be loaded into matlab/octave. + /// + /// Only nodes and node connectivity need be printed to visualise the graph. + void Print(const std::vector& nodes, const std::vector& edges, std::ostream& out = std::cout); + /// @brief Print the (simplified) graph in a form that can be loaded into matlab/octave. /// /// Only nodes and node connectivity need be printed to visualise the graph. diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp index fd0f303b2..f42100813 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp @@ -332,8 +332,6 @@ namespace meshkernel const double upperRightX, const double upperRightY) const { - range_check::CheckGreater(blockSizeX, 0.0, "X block size"); - range_check::CheckGreater(blockSizeY, 0.0, "Y block size"); const int numColumns = static_cast(std::ceil((upperRightX - originX) / blockSizeX)); if (numColumns <= 0) @@ -371,21 +369,21 @@ namespace meshkernel throw NotImplementedError("Projection value {} not supported", static_cast(m_projection)); } - Point CurvilinearGridRectangular::RotateUpperRightByAngle (const double originX, const double originY, const double upperRightX, const double upperRightY, const double angle) const + Point CurvilinearGridRectangular::RotateByAngle(const double originX, const double originY, + const double upperRightX, const double upperRightY, + const double cosAngle, + const double sinAngle) const { - const double cosA = std::cos(angle * constants::conversion::degToRad); - const double sinA = std::sin(angle * constants::conversion::degToRad); - if (m_projection == Projection::cartesian) { - Point translated (upperRightX - originX, upperRightY - originY); - return {cosA * translated.x - sinA * translated.y, sinA * translated.x + cosA * translated.y}; + Point translated(upperRightX - originX, upperRightY - originY); + return {cosAngle * translated.x - sinAngle * translated.y, sinAngle * translated.x + cosAngle * translated.y}; } else { - Cartesian3DPoint rotationPoint = SphericalToCartesian3D ({originX, originY}); - Cartesian3DPoint point3d = SphericalToCartesian3D ({upperRightX, upperRightY}); + Cartesian3DPoint rotationPoint = SphericalToCartesian3D({originX, originY}); + Cartesian3DPoint point3d = SphericalToCartesian3D({upperRightX, upperRightY}); // Normalize the rotation axis (Rodrigues' formula requires a unit vector) // Points are on Earth's surface, so the magnitude is earth_radius @@ -393,40 +391,30 @@ namespace meshkernel rotationPoint.y / constants::geometric::earth_radius, rotationPoint.z / constants::geometric::earth_radius}; - // std::cout << "rotationPoint " << rotationPoint.x << ", " << rotationPoint.y << ", " << rotationPoint.z << std::endl; - // std::cout << "point3d " << point3d.x << ", " << point3d.y << ", " << point3d.z << std::endl; - // std::cout << "k (unit axis) " << k.x << ", " << k.y << ", " << k.z << std::endl; - - // k · v + // k \cdot v double dot = k.x * point3d.x + k.y * point3d.y + k.z * point3d.z; - // k × v + // k \cross v Cartesian3DPoint crossProd = VectorProduct(k, point3d); - // std::cout << "crossProd " << crossProd.x << ", " << crossProd.y << ", " << crossProd.z << std::endl; - // Rodrigues' formula: v_rot = v·cos(θ) + (k × v)·sin(θ) + k·(k·v)·(1 - cos(θ)) - Cartesian3DPoint rotatedPoint3d = {point3d.x * cosA + crossProd.x * sinA + k.x * dot * (1.0 - cosA), - point3d.y * cosA + crossProd.y * sinA + k.y * dot * (1.0 - cosA), - point3d.z * cosA + crossProd.z * sinA + k.z * dot * (1.0 - cosA)}; + // Rodrigues formula: v_rot = v·cos(θ) + (k × v)·sin(θ) + k·(k·v)·(1 - cos(θ)) + Cartesian3DPoint rotatedPoint3d = {point3d.x * cosAngle + crossProd.x * sinAngle + k.x * dot * (1.0 - cosAngle), + point3d.y * cosAngle + crossProd.y * sinAngle + k.y * dot * (1.0 - cosAngle), + point3d.z * cosAngle + crossProd.z * sinAngle + k.z * dot * (1.0 - cosAngle)}; - // std::cout << "rotatedPoint3d " << rotatedPoint3d.x << ", " << rotatedPoint3d.y << ", " << rotatedPoint3d.z << std::endl; - - Point pointOnSphere = Cartesian3DToSpherical (rotatedPoint3d, upperRightX); - - std::cout << "pointOnSphere "<< pointOnSphere.x << ", " << pointOnSphere.y << std::endl; + Point pointOnSphere = Cartesian3DToSpherical(rotatedPoint3d, upperRightX); return pointOnSphere; } - } - lin_alg::Matrix CurvilinearGridRectangular::ComputeSpherical2(const int numColumns, - const int numRows, - const double originX, - const double originY, - const double angle, - const double blockSizeX, - const double blockSizeY) const + lin_alg::Matrix CurvilinearGridRectangular::ComputeSphericalOnExtension(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const { lin_alg::Matrix result = ComputeCartesian(numColumns, numRows, @@ -439,20 +427,21 @@ namespace meshkernel const auto numM = result.cols(); const auto numN = result.rows(); + const double cosAngle = std::cos(angle * constants::conversion::degToRad); + const double sinAngle = std::sin(angle * constants::conversion::degToRad); + for (Eigen::Index n = 0; n < numN; ++n) { for (Eigen::Index m = 0; m < numM; ++m) { - result (n,m) = RotateUpperRightByAngle (originX, originY, result (n,m).x, result (n,m).y, angle); + result(n, m) = RotateByAngle(originX, originY, result(n, m).x, result(n, m).y, cosAngle, sinAngle); } } - return result; } - std::unique_ptr CurvilinearGridRectangular::Compute(const double originX, const double originY, const double blockSizeX, @@ -464,36 +453,35 @@ namespace meshkernel range_check::CheckGreater(blockSizeX, 0.0, "X block size"); range_check::CheckGreater(blockSizeY, 0.0, "Y block size"); - Point rotatedUpperRight = RotateUpperRightByAngle (originX, originY, upperRightX, upperRightY, -angle); - - std::cout << "rotated distance: " << rotatedUpperRight.x - originX << ", " << rotatedUpperRight.y - originY << std::endl; - std::cout << "rotated deltas : " << (rotatedUpperRight.x - originX) / blockSizeX << ", " << (rotatedUpperRight.y - originY) / blockSizeY << std::endl; + const double cosAngle = std::cos(-angle * constants::conversion::degToRad); + const double sinAngle = std::sin(-angle * constants::conversion::degToRad); + // rotate the upper right, by -angle so that the grid is aligned with the axis + Point rotatedUpperRight = RotateByAngle(originX, originY, upperRightX, upperRightY, cosAngle, sinAngle); + // Now the number of cells in each direction can be computed. const int numColumns = static_cast(std::ceil((rotatedUpperRight.x - originX) / blockSizeX)); + const int numRows = ComputeNumRows(originY, rotatedUpperRight.y, blockSizeX, blockSizeY, m_projection); if (numColumns <= 0) { throw AlgorithmError("Number of columns cannot be <= 0"); } - std::cout << "rotatedUpperRight "<< rotatedUpperRight.x << ", " << rotatedUpperRight.y << std::endl; - - const int numRows = static_cast(std::ceil((rotatedUpperRight.y - originY) / blockSizeY)); - // const int numRows = ComputeNumRows(originY, rotatedUpperRight.y, blockSizeX, blockSizeY, m_projection); - std::cout << "num rows, cols "<< numRows << " " << numColumns << " "<< originX << ", " << originY << " "<< rotatedUpperRight.x << ", " << rotatedUpperRight.y << std::endl; + if (numRows <= 0) + { + throw AlgorithmError("Number of rows cannot be <= 0"); + } if (m_projection == Projection::spherical) { - // double deltax = (rotatedUpperRight.x - originX) / blockSizeX; - // double bsx = - auto grid = std::make_unique(ComputeSpherical(numColumns, - numRows, - originX, - originY, - angle, - blockSizeX, - blockSizeY), + auto grid = std::make_unique(ComputeSphericalOnExtension(numColumns, + numRows, + originX, + originY, + angle, + blockSizeX, + blockSizeY), m_projection); return grid; diff --git a/libs/MeshKernel/src/Utilities/Utilities.cpp b/libs/MeshKernel/src/Utilities/Utilities.cpp index 32151435e..4387e1851 100644 --- a/libs/MeshKernel/src/Utilities/Utilities.cpp +++ b/libs/MeshKernel/src/Utilities/Utilities.cpp @@ -58,6 +58,39 @@ void meshkernel::Print(const std::vector& nodes, const std::vector& } } +void meshkernel::Print(const std::vector& nodes, const std::vector& edges, std::ostream& out) +{ + out << "nullId = " << constants::missing::uintValue << ";" << std::endl; + out << "nullValue = " << constants::missing::doubleValue << ";" << std::endl; + out << "nodex = zeros ( " << nodes.size() << ", 1);" << std::endl; + out << "nodey = zeros ( " << nodes.size() << ", 1);" << std::endl; + out << "nodez = zeros ( " << nodes.size() << ", 1);" << std::endl; + out << "edges = zeros ( " << edges.size() << ", 2);" << std::endl; + + for (UInt i = 0; i < nodes.size(); ++i) + { + out << "nodex (" << i + 1 << " ) = " << nodes[i].x << ";" << std::endl; + } + + for (UInt i = 0; i < nodes.size(); ++i) + { + out << "nodey (" << i + 1 << " ) = " << nodes[i].y << ";" << std::endl; + } + + for (UInt i = 0; i < nodes.size(); ++i) + { + out << "nodez (" << i + 1 << " ) = " << nodes[i].z << ";" << std::endl; + } + + out << "edges = zeros ( " << edges.size() << ", 2 );" << std::endl; + + for (UInt i = 0; i < edges.size(); ++i) + { + out << "edges ( " << i + 1 << ", 1 ) = " << edges[i].first + 1 << ";" << std::endl; + out << "edges ( " << i + 1 << ", 2 ) = " << edges[i].second + 1 << ";" << std::endl; + } +} + void meshkernel::Print(const std::vector& xNodes, const std::vector& yNodes, const std::vector& edges, diff --git a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp index d463168ca..8fe6ad4c7 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp @@ -695,44 +695,42 @@ TEST_P(CurvilinearGridUniformTests, parameters) } INSTANTIATE_TEST_SUITE_P(curvilinearGridDeletionTests, CurvilinearGridUniformTests, ::testing::ValuesIn(CurvilinearGridUniformTests::GetData())); - -TEST(CurvilinearGridUniform, WithAngle) +TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_ShouldComputeCurvilinearGrid) { double originX = 10.0; double originY = 10.0; - double blockSizeX= 1.0; - double blockSizeY = 0.5; - // double upperRightX = 1.283012701892219e+01; - // double upperRightY = 1.509807621135332e+01; - - - double upperRightX = -0.7503; - double upperRightY = 19.2242; - // double upperRightX = 6.8579; - // double upperRightY = 14.8570; - double angle = 90.0; - - // double upperRightX = 12.3963; - // double upperRightY = 14.1524; - // double angle = 30.0; - + double blockSizeX = 2.5; + double blockSizeY = 2.5; + double upperRightX = 13.3295; + double upperRightY = 23.5238; + double angle = 30.0; + // Execution + CurvilinearGridRectangular curvilinearGridRectangular(Projection::spherical); + auto mesh = curvilinearGridRectangular.Compute(originX, originY, blockSizeX, blockSizeY, upperRightX, upperRightY, angle); + auto nodes = mesh->ComputeNodes(); - // double originX = 0.0; - // double originY = 0.0; - // double blockSizeX= 1.0; - // double blockSizeY = 0.5; - // double upperRightX = 0.0; - // double upperRightY = 5.0 * 1.414213562373095; - // double angle = 45.0; + std::vector expectedX{10.0, 12.16888013786, 14.3468065085138, 16.5359263214781, 18.7384207202835, 20.9565003943212, 8.72159804065856, 10.8828919460438, 13.0544076217419, + 15.2383407185244, 17.4369344284732, 19.6524758047091, 7.42224872698887, 9.57361905652058, 11.7363446582283, 13.9126732433567, 16.1049144936475, + 18.3154373418827, 6.09772407721882, 8.23669775793749, 10.388104315838, 12.5542463310128, 14.737503640127, 16.940331857308, 4.74354270565264, + 6.8674833359555, 9.0048626713965, 11.1580398538707, 13.329467428263, 15.5216914097679, 3.35492174705887, 5.46099947138391, 7.58143163777757, + 9.71863415455067, 11.8751332977645, 14.0535676553407}; + std::vector expectedY{10.0, 11.2317992926809, 12.4640006025882, 13.6948207289256, 14.9224495909188, 16.1450456084977, 12.1624690481718, 13.3931999286175, + 14.6253766863979, 15.8572318061806, 17.0869674455067, 18.3127505245545, 14.3189673215657, 15.5479790641342, 16.7794891846917, 18.0117518427796, + 19.2429876591642, 20.4713785004972, 16.4682677192212, 17.6948692356156, 18.9250292278688, 20.1570297171979, 21.3891162846345, 22.6194925439959, + 18.609057890161, 19.8325082588868, 21.0605833410122, 22.2915994013818, 23.523833753982, 24.7555189266162, 20.7399226564049, 21.9594204751702, + 23.1846133882867, 24.4138584680854, 25.6454718020147, 26.8777223650489}; - // Execution - CurvilinearGridRectangular curvilinearGridRectangular(Projection::spherical); - auto mesh = curvilinearGridRectangular.Compute (originX, originY, blockSizeX, blockSizeY, upperRightX, upperRightY, angle); + ASSERT_EQ(nodes.size(), expectedX.size()); - meshkernel::Print (mesh->ComputeNodes(), mesh->ComputeEdges ()); + constexpr double tolerance = 1.0e-10; + for (size_t i = 0; i < nodes.size(); ++i) + { + EXPECT_NEAR(nodes[i].x, expectedX[i], tolerance); + EXPECT_NEAR(nodes[i].y, expectedY[i], tolerance); + } } diff --git a/libs/MeshKernelApi/include/MeshKernelApi/Utils.hpp b/libs/MeshKernelApi/include/MeshKernelApi/Utils.hpp index 47d6c8ba9..59d027316 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/Utils.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/Utils.hpp @@ -531,17 +531,25 @@ namespace meshkernelapi { const meshkernel::CurvilinearGridRectangular grid(projection); - if (!meshkernel::IsEqual(makeGridParameters.angle, 0.0)) + if (makeGridParameters.angle == 0.0) { - throw meshkernel::AlgorithmError("When generating an uniform grid on an defined extension, the grid angle must be equal to 0"); + return grid.Compute(makeGridParameters.origin_x, + makeGridParameters.origin_y, + makeGridParameters.block_size_x, + makeGridParameters.block_size_y, + makeGridParameters.upper_right_x, + makeGridParameters.upper_right_y); + } + else + { + return grid.Compute(makeGridParameters.origin_x, + makeGridParameters.origin_y, + makeGridParameters.block_size_x, + makeGridParameters.block_size_y, + makeGridParameters.upper_right_x, + makeGridParameters.upper_right_y, + makeGridParameters.angle); } - - return grid.Compute(makeGridParameters.origin_x, - makeGridParameters.origin_y, - makeGridParameters.block_size_x, - makeGridParameters.block_size_y, - makeGridParameters.upper_right_x, - makeGridParameters.upper_right_y); } template From c4ebbea7cc4587988f0fd71d6da04276af69f8e3 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 19 Aug 2026 11:24:01 +0000 Subject: [PATCH 3/7] ci: retry macOS artifact upload on transient network failures Co-authored-by: BillSenior <89970704+BillSenior@users.noreply.github.com> --- .github/workflows/build-and-test-workflow.yml | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/.github/workflows/build-and-test-workflow.yml b/.github/workflows/build-and-test-workflow.yml index 8859aa844..2a01b66ec 100644 --- a/.github/workflows/build-and-test-workflow.yml +++ b/.github/workflows/build-and-test-workflow.yml @@ -156,9 +156,20 @@ jobs: # Step: Upload artifacts - name: Upload artifacts + id: upload_artifacts uses: actions/upload-artifact@v7 if: always() + continue-on-error: true with: name: meshkernel-${{ inputs.platform }}-${{ inputs.build_type }} path: ${{ steps.paths.outputs.install_dir }} if-no-files-found: error + + - name: Retry upload artifacts on transient failure + uses: actions/upload-artifact@v7 + if: always() && steps.upload_artifacts.outcome == 'failure' + with: + name: meshkernel-${{ inputs.platform }}-${{ inputs.build_type }} + path: ${{ steps.paths.outputs.install_dir }} + if-no-files-found: error + overwrite: true From 96c86e5c3c1f738804efd5af9bc971fe75634a5f Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Thu, 27 Aug 2026 11:39:28 +0200 Subject: [PATCH 4/7] Passing to other computer --- .../CurvilinearGridRectangular.hpp | 37 ++- .../CurvilinearGridRectangular.cpp | 236 +++++++++++++++++- .../src/CurvilinearGridRectangularTests.cpp | 52 ++++ 3 files changed, 314 insertions(+), 11 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp index 54a984136..f3c49860c 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp @@ -133,14 +133,19 @@ namespace meshkernel /// @param[in] blockSizeX The grid block size in x dimension /// @param[in] blockSizeY The grid block size in y dimension /// @returns[in] The coordinates of the grid point - static lin_alg::Matrix ComputeSpherical(const int numColumns, + lin_alg::Matrix ComputeSpherical(const int numColumns, const int numRows, const double originX, const double originY, const double angle, const double blockSizeX, - const double blockSizeY); + const double blockSizeY) const; + // Preserves Orthogonality: Standard flat 2D rotations introduce a slight + // shearing skew over a curved globe because degrees of longitude shrink + // away from the equator. Rodrigues' 3D rotation preserves the exact + // geometry of your initial unrotated grid matrix relative to the surface + // plane. lin_alg::Matrix ComputeSphericalOnExtension(const int numColumns, const int numRows, const double originX, @@ -149,6 +154,34 @@ namespace meshkernel const double blockSizeX, const double blockSizeY) const; + // Not really RGF grid algorithm, there seems to be a step missing + lin_alg::Matrix ComputeSphericalRgfGrid(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const; + + // Generate the rotated grid using fixed delta-x and delta-y (metres) when generating + lin_alg::Matrix ComputeSphericalBoostGrid(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const; + + static Cartesian3DPoint RotateVectorRodrigues(const Cartesian3DPoint& v, const Cartesian3DPoint& k, double theta_rad); + + lin_alg::Matrix ComputeSphericalMercator(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const; + /// @brief Compute the adjusted latitude for keeping an aspect ratio of 1, considering the spherical coordinates /// @param[in] blockSize The grid block size in y dimension /// @param[in] aspectRatio The expected element aspect ratio, defined as being: blockSizeY / blockSizeX diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp index f42100813..7708d9a54 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp @@ -25,6 +25,14 @@ // //------------------------------------------------------------------------------ +// include boost +#define BOOST_ALLOW_DEPRECATED_HEADERS +#include +#include +#include +#undef BOOST_ALLOW_DEPRECATED_HEADERS + + #include "MeshKernel/Exceptions.hpp" #include @@ -62,7 +70,11 @@ namespace meshkernel if (m_projection == Projection::spherical) { - return std::make_unique(ComputeSpherical(numColumns, + return std::make_unique(ComputeSphericalMercator(numColumns, + // return std::make_unique(ComputeSphericalBoostGrid(numColumns, + // return std::make_unique(ComputeSphericalRgfGrid(numColumns, + // return std::make_unique(ComputeSphericalOnExtension(numColumns, + // return std::make_unique(ComputeSpherical(numColumns, numRows, originX, originY, @@ -124,15 +136,24 @@ namespace meshkernel const double originY, const double angle, const double blockSizeX, - const double blockSizeY) + const double blockSizeY) const { - lin_alg::Matrix result = ComputeCartesian(numColumns, - numRows, - originX, - originY, - angle, - blockSizeX, - blockSizeY); + + lin_alg::Matrix result = ComputeSphericalOnExtension (numColumns, + numRows, + originX, + originY, + angle, + blockSizeX, + blockSizeY); + + // lin_alg::Matrix result = ComputeCartesian(numColumns, + // numRows, + // originX, + // originY, + // angle, + // blockSizeX, + // blockSizeY); const auto numM = result.cols(); const auto numN = result.rows(); @@ -442,6 +463,203 @@ namespace meshkernel return result; } + lin_alg::Matrix CurvilinearGridRectangular::ComputeSphericalRgfGrid(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const + { + lin_alg::Matrix result = ComputeCartesian(numColumns, + numRows, + originX, + originY, + angle, + blockSizeX, + blockSizeY); + + // const auto numM = result.cols(); + // const auto numN = result.rows(); + + // const double cosAngle = std::cos(angle * constants::conversion::degToRad); + // const double sinAngle = std::sin(angle * constants::conversion::degToRad); + + // for (Eigen::Index n = 0; n < numN; ++n) + // { + + // for (Eigen::Index m = 0; m < numM; ++m) + // { + // result(n, m) = RotateByAngle(originX, originY, result(n, m).x, result(n, m).y, cosAngle, sinAngle); + // } + // } + + return result; + } + + + lin_alg::Matrix CurvilinearGridRectangular::ComputeSphericalBoostGrid(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeXDeg, + const double blockSizeYDeg) const + { + namespace bg = boost::geometry; + + const int numM = numColumns + 1; + const int numN = numRows + 1; + + using Point2D = bg::model::point>; + + lin_alg::Matrix result(numN, numM); + + double angle_y_rad = angle * constants::conversion::degToRad; + double angle_x_rad = (angle + 90.0) * constants::conversion::degToRad; + + bg::srs::spheroid earth_spheroid(constants::geometric::earth_radius, constants::geometric::earth_radius); + using thomas_type = bg::formula::thomas_direct; + + double blockSizeX = blockSizeXDeg * 111000.0 * std::cos (originY * constants::conversion::degToRad); + double blockSizeY = blockSizeYDeg * 111000.0 * std::cos (originY * constants::conversion::degToRad); + + // using FormulaStrategy = bg::strategy::formula::thomas; + + for (Eigen::Index n = 0; n < numN; ++n) + { + + Point2D current_row_origin; + + if (n == 0) { + current_row_origin = Point2D(originX, originY); + } else { + // Step the next row along the rotated Y-axis vector (heading = angle_y_rad) + double prev_row_lon_rad = result(n-1, 0).x * bg::math::d2r(); + double prev_row_lat_rad = result(n-1, 0).y * bg::math::d2r(); + + auto dir_y = thomas_type::apply (prev_row_lon_rad, prev_row_lat_rad, blockSizeY, angle_y_rad, earth_spheroid); + + current_row_origin = Point2D(dir_y.lon2 * bg::math::r2d(), + dir_y.lat2 * bg::math::r2d()); + } + + for (Eigen::Index m = 0; m < numM; ++m) + { + if (m == 0) { + result(n, m).x = bg::get<0>(current_row_origin); + result(n, m).y = bg::get<1>(current_row_origin); + } else { + + [[maybe_unused]] Point cur = result(n, m-1); + + // Step columns outward along the perpendicular rotated X-axis vector (heading = angle_x_rad) + double current_lon_rad = result(n, m-1).x * bg::math::d2r(); + double current_lat_rad = result(n, m-1).y * bg::math::d2r(); + + + auto dir_x = thomas_type::apply (current_lon_rad, current_lat_rad, blockSizeX, angle_x_rad, earth_spheroid); + + result(n, m).x = dir_x.lon2 * bg::math::r2d(); + result(n, m).y = dir_x.lat2 * bg::math::r2d(); + } + } + } + + return result; + } + + Cartesian3DPoint CurvilinearGridRectangular::RotateVectorRodrigues(const Cartesian3DPoint& v, const Cartesian3DPoint& k, double theta_rad) { + + double cos_t = std::cos(theta_rad); + double sin_t = std::sin(theta_rad); + + Cartesian3DPoint cross = { k.y * v.z - k.z * v.y, k.z * v.x - k.x * v.z, k.x * v.y - k.y * v.x }; + + double dot = k.x * v.x + k.y * v.y + k.z * v.z; + + return { + v.x * cos_t + cross.x * sin_t + k.x * dot * (1.0 - cos_t), + v.y * cos_t + cross.y * sin_t + k.y * dot * (1.0 - cos_t), + v.z * cos_t + cross.z * sin_t + k.z * dot * (1.0 - cos_t) + }; + } + + lin_alg::Matrix CurvilinearGridRectangular::ComputeSphericalMercator(const int numColumns, + const int numRows, + const double origin_lon, + const double origin_lat, + const double rotation_deg, + const double d_lon, + const double d_lat) const + { + const int ny = numColumns + 1; + const int nx = numRows + 1; + + lin_alg::Matrix result(ny, nx); + + const double theta_rad = rotation_deg * (M_PI / 180.0); + Cartesian3DPoint k = SphericalToCartesian3D(Point(origin_lon, origin_lat)); + + double kLength = std::sqrt (k.x * k.x + k.y * k.y + k.z * k.z); + + k.x /= kLength; + k.y /= kLength; + k.z /= kLength; + + // 1. Transform starting origin latitude into Conformal Mercator space + double origin_lat_rad = origin_lat * (M_PI / 180.0); + double origin_y_mercator = std::log(std::tan(M_PI / 4.0 + origin_lat_rad / 2.0)); + + // 2. Compute the precise scaling factor for your rectangular dimensions. + // Instead of forcing square steps, we scale the vertical Mercator increment + // by the exact ratio of your desired d_lat to d_lon. + double d_step_lon_rad = d_lon * (M_PI / 180.0); + double d_step_lat_conformal = d_step_lon_rad * (d_lat / d_lon); + + for (int j = 0; j < ny; ++j) { + // Step along the scaled conformal Y axis + double current_y_mercator = origin_y_mercator + (j * d_step_lat_conformal); + + // Inverse Mercator equation back to physical latitude + double unrotated_lat_rad = 2.0 * std::atan(std::exp(current_y_mercator)) - M_PI / 2.0; + double unrotated_lat = unrotated_lat_rad * (180.0 / M_PI); + + for (int i = 0; i < nx; ++i) { + // Step uniformly along the X axis + double unrotated_lon = origin_lon + (i * d_lon); + + // Convert conformal unrotated point into 3D space + Cartesian3DPoint p_initial = SphericalToCartesian3D(Point(unrotated_lon, unrotated_lat)); + + // Apply 3D Rodrigues rotation + Cartesian3DPoint p_rotated = RotateVectorRodrigues(p_initial, k, theta_rad); + + // Project back to degrees and update the Eigen Matrix + double final_lon, final_lat; + + final_lat = std::asin(std::max(-1.0, std::min (1.0, p_rotated.z))) * (180.0 / M_PI); + final_lon = std::atan2(p_rotated.y, p_rotated.x) * (180.0 / M_PI); + + // std::cout << current_y_mercator << " "<< k.x << ", " << k.y << ", " << k.z << " " << final_lat << " " << final_lon << " " << rotation_deg << " " << d_lon << " " << d_lat << std::endl; + // std::cout << current_y_mercator << " "<< p_initial.x << ", " << p_initial.y << ", " << p_initial.z << " " << final_lat << " " << final_lon << " " << rotation_deg << " " << d_lon << " " << d_lat << std::endl; + // std::cout << current_y_mercator << " "<< p_rotated.x << ", " << p_rotated.y << ", " << p_rotated.z << " " << final_lat << " " << final_lon << " " << rotation_deg << " " << d_lon << " " << d_lat << " " << theta_rad << std::endl; + + //CartesianToLatLon(p_rotated, final_lon, final_lat); + + result(j, i).x = final_lon; + result(j, i).y = final_lat; + } + } + + + return result; + } + + + + + std::unique_ptr CurvilinearGridRectangular::Compute(const double originX, const double originY, const double blockSizeX, diff --git a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp index 8fe6ad4c7..fb32031e8 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp @@ -43,6 +43,9 @@ #include #include +#include + + using namespace meshkernel; TEST(CurvilinearGridUniform, CurvilinearGridRectangular_WithPolygon_ShouldComputeCurvilinearGrid) @@ -698,6 +701,53 @@ INSTANTIATE_TEST_SUITE_P(curvilinearGridDeletionTests, CurvilinearGridUniformTes TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_ShouldComputeCurvilinearGrid) { + [[maybe_unused]] double originX = 105.7; + [[maybe_unused]] double originY = 18.2; + [[maybe_unused]] double blockSizeX = 0.1; + [[maybe_unused]] double blockSizeY = 0.1; + [[maybe_unused]] double upperRightX = 13.3295; + [[maybe_unused]] double upperRightY = 23.5238; + [[maybe_unused]] double angle = -43.0; + + // [[maybe_unused]] double originX = 10.0; + // [[maybe_unused]] double originY = 10.0; + // [[maybe_unused]] double blockSizeX = 2.5; + // [[maybe_unused]] double blockSizeY = 2.5; + // [[maybe_unused]] double upperRightX = 13.3295; + // [[maybe_unused]] double upperRightY = 23.5238; + // [[maybe_unused]] double angle = 60.0; + + // Execution + CurvilinearGridRectangular curvilinearGridRectangular(Projection::spherical); + auto mesh = curvilinearGridRectangular.Compute(45, 45, originX, originY, angle, blockSizeX, blockSizeY); + // auto mesh = curvilinearGridRectangular.Compute(originX, originY, blockSizeX, blockSizeY, upperRightX, upperRightY, angle); + + auto nodes = mesh->ComputeNodes(); + + meshkernel::Print (mesh->ComputeNodes (), mesh->ComputeEdges ()); + + meshkernel::Mesh2D mesh2d (mesh->ComputeEdges (), mesh->ComputeNodes (), mesh->projection ()); + auto ortho = meshkernel::MeshOrthogonality::Compute (mesh2d); + std::ranges::sort (ortho); + + for (size_t i = 0; i < ortho.size (); ++i) + { + + if (i % 10 == 0) std::cout << std::endl; + + std::cout << std::setw (15) << " " << ortho[i] << " "; + + + } + + + double max = *std::max_element (ortho.begin (), ortho.end ()); + + std::cout << "Maximum orthoginality = " << max << std::endl; + + return; + /* + double originX = 10.0; double originY = 10.0; double blockSizeX = 2.5; @@ -712,6 +762,7 @@ TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_Sho auto nodes = mesh->ComputeNodes(); + std::vector expectedX{10.0, 12.16888013786, 14.3468065085138, 16.5359263214781, 18.7384207202835, 20.9565003943212, 8.72159804065856, 10.8828919460438, 13.0544076217419, 15.2383407185244, 17.4369344284732, 19.6524758047091, 7.42224872698887, 9.57361905652058, 11.7363446582283, 13.9126732433567, 16.1049144936475, 18.3154373418827, 6.09772407721882, 8.23669775793749, 10.388104315838, 12.5542463310128, 14.737503640127, 16.940331857308, 4.74354270565264, @@ -733,4 +784,5 @@ TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_Sho EXPECT_NEAR(nodes[i].x, expectedX[i], tolerance); EXPECT_NEAR(nodes[i].y, expectedY[i], tolerance); } + */ } From 0f83495309c9266581007275c8e8b67842cc41ec Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 31 Aug 2026 12:55:07 +0200 Subject: [PATCH 5/7] GRIDEDIT-2303 Removed unused algorithms --- .../CurvilinearGridRectangular.hpp | 31 +-- .../CurvilinearGridRectangular.cpp | 210 ++---------------- .../src/CurvilinearGridRectangularTests.cpp | 25 ++- 3 files changed, 34 insertions(+), 232 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp index f3c49860c..c2ce25466 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp @@ -134,12 +134,12 @@ namespace meshkernel /// @param[in] blockSizeY The grid block size in y dimension /// @returns[in] The coordinates of the grid point lin_alg::Matrix ComputeSpherical(const int numColumns, - const int numRows, - const double originX, - const double originY, - const double angle, - const double blockSizeX, - const double blockSizeY) const; + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const; // Preserves Orthogonality: Standard flat 2D rotations introduce a slight // shearing skew over a curved globe because degrees of longitude shrink @@ -163,25 +163,6 @@ namespace meshkernel const double blockSizeX, const double blockSizeY) const; - // Generate the rotated grid using fixed delta-x and delta-y (metres) when generating - lin_alg::Matrix ComputeSphericalBoostGrid(const int numColumns, - const int numRows, - const double originX, - const double originY, - const double angle, - const double blockSizeX, - const double blockSizeY) const; - - static Cartesian3DPoint RotateVectorRodrigues(const Cartesian3DPoint& v, const Cartesian3DPoint& k, double theta_rad); - - lin_alg::Matrix ComputeSphericalMercator(const int numColumns, - const int numRows, - const double originX, - const double originY, - const double angle, - const double blockSizeX, - const double blockSizeY) const; - /// @brief Compute the adjusted latitude for keeping an aspect ratio of 1, considering the spherical coordinates /// @param[in] blockSize The grid block size in y dimension /// @param[in] aspectRatio The expected element aspect ratio, defined as being: blockSizeY / blockSizeX diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp index 7708d9a54..6f22b40d2 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp @@ -25,14 +25,6 @@ // //------------------------------------------------------------------------------ -// include boost -#define BOOST_ALLOW_DEPRECATED_HEADERS -#include -#include -#include -#undef BOOST_ALLOW_DEPRECATED_HEADERS - - #include "MeshKernel/Exceptions.hpp" #include @@ -70,17 +62,14 @@ namespace meshkernel if (m_projection == Projection::spherical) { - return std::make_unique(ComputeSphericalMercator(numColumns, - // return std::make_unique(ComputeSphericalBoostGrid(numColumns, - // return std::make_unique(ComputeSphericalRgfGrid(numColumns, - // return std::make_unique(ComputeSphericalOnExtension(numColumns, // return std::make_unique(ComputeSpherical(numColumns, - numRows, - originX, - originY, - angle, - blockSizeX, - blockSizeY), + return std::make_unique(ComputeSphericalOnExtension(numColumns, + numRows, + originX, + originY, + angle, + blockSizeX, + blockSizeY), m_projection); } if (m_projection == Projection::cartesian) @@ -139,21 +128,13 @@ namespace meshkernel const double blockSizeY) const { - lin_alg::Matrix result = ComputeSphericalOnExtension (numColumns, - numRows, - originX, - originY, - angle, - blockSizeX, - blockSizeY); - - // lin_alg::Matrix result = ComputeCartesian(numColumns, - // numRows, - // originX, - // originY, - // angle, - // blockSizeX, - // blockSizeY); + lin_alg::Matrix result = ComputeCartesian(numColumns, + numRows, + originX, + originY, + angle, + blockSizeX, + blockSizeY); const auto numM = result.cols(); const auto numN = result.rows(); @@ -497,169 +478,6 @@ namespace meshkernel return result; } - - lin_alg::Matrix CurvilinearGridRectangular::ComputeSphericalBoostGrid(const int numColumns, - const int numRows, - const double originX, - const double originY, - const double angle, - const double blockSizeXDeg, - const double blockSizeYDeg) const - { - namespace bg = boost::geometry; - - const int numM = numColumns + 1; - const int numN = numRows + 1; - - using Point2D = bg::model::point>; - - lin_alg::Matrix result(numN, numM); - - double angle_y_rad = angle * constants::conversion::degToRad; - double angle_x_rad = (angle + 90.0) * constants::conversion::degToRad; - - bg::srs::spheroid earth_spheroid(constants::geometric::earth_radius, constants::geometric::earth_radius); - using thomas_type = bg::formula::thomas_direct; - - double blockSizeX = blockSizeXDeg * 111000.0 * std::cos (originY * constants::conversion::degToRad); - double blockSizeY = blockSizeYDeg * 111000.0 * std::cos (originY * constants::conversion::degToRad); - - // using FormulaStrategy = bg::strategy::formula::thomas; - - for (Eigen::Index n = 0; n < numN; ++n) - { - - Point2D current_row_origin; - - if (n == 0) { - current_row_origin = Point2D(originX, originY); - } else { - // Step the next row along the rotated Y-axis vector (heading = angle_y_rad) - double prev_row_lon_rad = result(n-1, 0).x * bg::math::d2r(); - double prev_row_lat_rad = result(n-1, 0).y * bg::math::d2r(); - - auto dir_y = thomas_type::apply (prev_row_lon_rad, prev_row_lat_rad, blockSizeY, angle_y_rad, earth_spheroid); - - current_row_origin = Point2D(dir_y.lon2 * bg::math::r2d(), - dir_y.lat2 * bg::math::r2d()); - } - - for (Eigen::Index m = 0; m < numM; ++m) - { - if (m == 0) { - result(n, m).x = bg::get<0>(current_row_origin); - result(n, m).y = bg::get<1>(current_row_origin); - } else { - - [[maybe_unused]] Point cur = result(n, m-1); - - // Step columns outward along the perpendicular rotated X-axis vector (heading = angle_x_rad) - double current_lon_rad = result(n, m-1).x * bg::math::d2r(); - double current_lat_rad = result(n, m-1).y * bg::math::d2r(); - - - auto dir_x = thomas_type::apply (current_lon_rad, current_lat_rad, blockSizeX, angle_x_rad, earth_spheroid); - - result(n, m).x = dir_x.lon2 * bg::math::r2d(); - result(n, m).y = dir_x.lat2 * bg::math::r2d(); - } - } - } - - return result; - } - - Cartesian3DPoint CurvilinearGridRectangular::RotateVectorRodrigues(const Cartesian3DPoint& v, const Cartesian3DPoint& k, double theta_rad) { - - double cos_t = std::cos(theta_rad); - double sin_t = std::sin(theta_rad); - - Cartesian3DPoint cross = { k.y * v.z - k.z * v.y, k.z * v.x - k.x * v.z, k.x * v.y - k.y * v.x }; - - double dot = k.x * v.x + k.y * v.y + k.z * v.z; - - return { - v.x * cos_t + cross.x * sin_t + k.x * dot * (1.0 - cos_t), - v.y * cos_t + cross.y * sin_t + k.y * dot * (1.0 - cos_t), - v.z * cos_t + cross.z * sin_t + k.z * dot * (1.0 - cos_t) - }; - } - - lin_alg::Matrix CurvilinearGridRectangular::ComputeSphericalMercator(const int numColumns, - const int numRows, - const double origin_lon, - const double origin_lat, - const double rotation_deg, - const double d_lon, - const double d_lat) const - { - const int ny = numColumns + 1; - const int nx = numRows + 1; - - lin_alg::Matrix result(ny, nx); - - const double theta_rad = rotation_deg * (M_PI / 180.0); - Cartesian3DPoint k = SphericalToCartesian3D(Point(origin_lon, origin_lat)); - - double kLength = std::sqrt (k.x * k.x + k.y * k.y + k.z * k.z); - - k.x /= kLength; - k.y /= kLength; - k.z /= kLength; - - // 1. Transform starting origin latitude into Conformal Mercator space - double origin_lat_rad = origin_lat * (M_PI / 180.0); - double origin_y_mercator = std::log(std::tan(M_PI / 4.0 + origin_lat_rad / 2.0)); - - // 2. Compute the precise scaling factor for your rectangular dimensions. - // Instead of forcing square steps, we scale the vertical Mercator increment - // by the exact ratio of your desired d_lat to d_lon. - double d_step_lon_rad = d_lon * (M_PI / 180.0); - double d_step_lat_conformal = d_step_lon_rad * (d_lat / d_lon); - - for (int j = 0; j < ny; ++j) { - // Step along the scaled conformal Y axis - double current_y_mercator = origin_y_mercator + (j * d_step_lat_conformal); - - // Inverse Mercator equation back to physical latitude - double unrotated_lat_rad = 2.0 * std::atan(std::exp(current_y_mercator)) - M_PI / 2.0; - double unrotated_lat = unrotated_lat_rad * (180.0 / M_PI); - - for (int i = 0; i < nx; ++i) { - // Step uniformly along the X axis - double unrotated_lon = origin_lon + (i * d_lon); - - // Convert conformal unrotated point into 3D space - Cartesian3DPoint p_initial = SphericalToCartesian3D(Point(unrotated_lon, unrotated_lat)); - - // Apply 3D Rodrigues rotation - Cartesian3DPoint p_rotated = RotateVectorRodrigues(p_initial, k, theta_rad); - - // Project back to degrees and update the Eigen Matrix - double final_lon, final_lat; - - final_lat = std::asin(std::max(-1.0, std::min (1.0, p_rotated.z))) * (180.0 / M_PI); - final_lon = std::atan2(p_rotated.y, p_rotated.x) * (180.0 / M_PI); - - // std::cout << current_y_mercator << " "<< k.x << ", " << k.y << ", " << k.z << " " << final_lat << " " << final_lon << " " << rotation_deg << " " << d_lon << " " << d_lat << std::endl; - // std::cout << current_y_mercator << " "<< p_initial.x << ", " << p_initial.y << ", " << p_initial.z << " " << final_lat << " " << final_lon << " " << rotation_deg << " " << d_lon << " " << d_lat << std::endl; - // std::cout << current_y_mercator << " "<< p_rotated.x << ", " << p_rotated.y << ", " << p_rotated.z << " " << final_lat << " " << final_lon << " " << rotation_deg << " " << d_lon << " " << d_lat << " " << theta_rad << std::endl; - - //CartesianToLatLon(p_rotated, final_lon, final_lat); - - result(j, i).x = final_lon; - result(j, i).y = final_lat; - } - } - - - return result; - } - - - - - std::unique_ptr CurvilinearGridRectangular::Compute(const double originX, const double originY, const double blockSizeX, diff --git a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp index fb32031e8..f218582ad 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp @@ -45,7 +45,6 @@ #include - using namespace meshkernel; TEST(CurvilinearGridUniform, CurvilinearGridRectangular_WithPolygon_ShouldComputeCurvilinearGrid) @@ -701,12 +700,17 @@ INSTANTIATE_TEST_SUITE_P(curvilinearGridDeletionTests, CurvilinearGridUniformTes TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_ShouldComputeCurvilinearGrid) { +#if 0 + + [[maybe_unused]] double blockSizeX = 0.02; + [[maybe_unused]] double blockSizeY = 0.02; + [[maybe_unused]] double originX = 105.7; - [[maybe_unused]] double originY = 18.2; - [[maybe_unused]] double blockSizeX = 0.1; - [[maybe_unused]] double blockSizeY = 0.1; - [[maybe_unused]] double upperRightX = 13.3295; - [[maybe_unused]] double upperRightY = 23.5238; + [[maybe_unused]] double originY = 18.2; + [[maybe_unused]] double upperRightX = 106.79; + [[maybe_unused]] double upperRightY = 18.46; + // [[maybe_unused]] double upperRightX = 13.3295; + // [[maybe_unused]] double upperRightY = 23.5238; [[maybe_unused]] double angle = -43.0; // [[maybe_unused]] double originX = 10.0; @@ -719,8 +723,8 @@ TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_Sho // Execution CurvilinearGridRectangular curvilinearGridRectangular(Projection::spherical); - auto mesh = curvilinearGridRectangular.Compute(45, 45, originX, originY, angle, blockSizeX, blockSizeY); - // auto mesh = curvilinearGridRectangular.Compute(originX, originY, blockSizeX, blockSizeY, upperRightX, upperRightY, angle); + // auto mesh = curvilinearGridRectangular.Compute(45, 45, originX, originY, angle, blockSizeX, blockSizeY); + auto mesh = curvilinearGridRectangular.Compute(originX, originY, blockSizeX, blockSizeY, upperRightX, upperRightY, angle); auto nodes = mesh->ComputeNodes(); @@ -746,7 +750,8 @@ TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_Sho std::cout << "Maximum orthoginality = " << max << std::endl; return; - /* + +#endif double originX = 10.0; double originY = 10.0; @@ -762,7 +767,6 @@ TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_Sho auto nodes = mesh->ComputeNodes(); - std::vector expectedX{10.0, 12.16888013786, 14.3468065085138, 16.5359263214781, 18.7384207202835, 20.9565003943212, 8.72159804065856, 10.8828919460438, 13.0544076217419, 15.2383407185244, 17.4369344284732, 19.6524758047091, 7.42224872698887, 9.57361905652058, 11.7363446582283, 13.9126732433567, 16.1049144936475, 18.3154373418827, 6.09772407721882, 8.23669775793749, 10.388104315838, 12.5542463310128, 14.737503640127, 16.940331857308, 4.74354270565264, @@ -784,5 +788,4 @@ TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_Sho EXPECT_NEAR(nodes[i].x, expectedX[i], tolerance); EXPECT_NEAR(nodes[i].y, expectedY[i], tolerance); } - */ } From 6a8a1c600b90062d96d263ddea4f53e0a7da084e Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 31 Aug 2026 15:36:06 +0200 Subject: [PATCH 6/7] GRIDEDIT-2303 Implemented a fixed delta grid generation with rotation --- .../CurvilinearGridRectangular.hpp | 9 ++ .../CurvilinearGridRectangular.cpp | 105 ++++++++++++++++-- 2 files changed, 106 insertions(+), 8 deletions(-) diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp index c2ce25466..ce3014599 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp @@ -163,6 +163,15 @@ namespace meshkernel const double blockSizeX, const double blockSizeY) const; + // Generate the rotated grid using fixed delta-x and delta-y when generating + lin_alg::Matrix ComputeSphericalFixedDelta(const int numColumns, + const int numRows, + const double originX, + const double originY, + const double angle, + const double blockSizeX, + const double blockSizeY) const; + /// @brief Compute the adjusted latitude for keeping an aspect ratio of 1, considering the spherical coordinates /// @param[in] blockSize The grid block size in y dimension /// @param[in] aspectRatio The expected element aspect ratio, defined as being: blockSizeY / blockSizeX diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp index 6f22b40d2..7ba86b491 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp @@ -62,14 +62,16 @@ namespace meshkernel if (m_projection == Projection::spherical) { - // return std::make_unique(ComputeSpherical(numColumns, - return std::make_unique(ComputeSphericalOnExtension(numColumns, - numRows, - originX, - originY, - angle, - blockSizeX, - blockSizeY), + // return std::make_unique(ComputeSphericalFixedDelta(numColumns, + // return std::make_unique(ComputeSphericalRgfGrid(numColumns, + // return std::make_unique(ComputeSphericalOnExtension(numColumns, + return std::make_unique(ComputeSpherical(numColumns, + numRows, + originX, + originY, + angle, + blockSizeX, + blockSizeY), m_projection); } if (m_projection == Projection::cartesian) @@ -478,6 +480,91 @@ namespace meshkernel return result; } + lin_alg::Matrix CurvilinearGridRectangular::ComputeSphericalFixedDelta(const int numColumns, + const int numRows, + const double origin_lon, + const double origin_lat, + const double rotation_deg, + const double d_lon, + const double d_lat) const + { + const int numM = numColumns + 1; + const int numN = numRows + 1; + + lin_alg::Matrix result(numN, numM); + + // Convert all angles to radians up front + const double lon_orig_rad = origin_lon * (M_PI / 180.0); + const double lat_orig_rad = origin_lat * (M_PI / 180.0); + const double rot_rad = rotation_deg * (M_PI / 180.0); + const double d_lon_rad = d_lon * (M_PI / 180.0); + const double d_lat_rad = d_lat * (M_PI / 180.0); + + // Conformal FixedDelta step size calculation relative to an unrotated equator (lat = 0) + // At the equator, cos(0) = 1, so the conformal step matches the physical radian step. + double d_step_lat_conformal = d_lat_rad; + + // We center our unrotated grid around (0,0) space + // so that it spins perfectly around its true local center. + // double half_width = (numRows * d_lon_rad) / 2.0; + + // Conformal latitude center tracker + // double origin_y_mercator = std::log(std::tan(M_PI / 4.0 + 0.0 / 2.0)); // Equator = 0 + + for (int j = 0; j < numN; ++j) + { + // Local relative latitude centered around 0 (Equator) + double current_y_mercator = 0.0 + ((j - numColumns / 2.0) * d_step_lat_conformal); + double local_lat_rad = 2.0 * std::atan(std::exp(current_y_mercator)) - M_PI / 2.0; + + for (int i = 0; i < numM; ++i) + { + // Local relative longitude centered around 0 + double local_lon_rad = (i - numRows / 2.0) * d_lon_rad; + + // 1. Convert local relative point to 3D Cartesian + double cos_local_lat = std::cos(local_lat_rad); + Cartesian3DPoint p; + p.x = cos_local_lat * std::cos(local_lon_rad); + p.y = cos_local_lat * std::sin(local_lon_rad); + p.z = std::sin(local_lat_rad); + + // 2. Twist the grid locally by rotation_deg around its own center axis (1, 0, 0) + // Since the local center is at (0,0), its Cartesian vector is pointing down the X axis! + if (rotation_deg != 0.0) + { + double ty = p.y * std::cos(rot_rad) - p.z * std::sin(rot_rad); + double tz = p.y * std::sin(rot_rad) + p.z * std::cos(rot_rad); + p.y = ty; + p.z = tz; + } + + // 3. Pitch: Move up/down to the target latitude. + // We rotate around the Y axis to change the latitude (X and Z change). + double px1 = p.x * std::cos(lat_orig_rad) - p.z * std::sin(lat_orig_rad); + double pz1 = p.x * std::sin(lat_orig_rad) + p.z * std::cos(lat_orig_rad); + p.x = px1; + p.z = pz1; + + // 4. Yaw: Spin over to the target longitude. + // We rotate around the Z axis to change longitude (X and Y change). + double px2 = p.x * std::cos(lon_orig_rad) - p.y * std::sin(lon_orig_rad); + double py2 = p.x * std::sin(lon_orig_rad) + p.y * std::cos(lon_orig_rad); + p.x = px2; + p.y = py2; + + // 5. Project back to geodetic degrees + double final_lat = std::asin(std::max(-1.0, std::min(1.0, p.z))) * (180.0 / M_PI); + double final_lon = std::atan2(p.y, p.x) * (180.0 / M_PI); + + result(j, i).x = final_lon; + result(j, i).y = final_lat; + } + } + + return result; + } + std::unique_ptr CurvilinearGridRectangular::Compute(const double originX, const double originY, const double blockSizeX, @@ -511,6 +598,8 @@ namespace meshkernel if (m_projection == Projection::spherical) { + // auto grid = std::make_unique(ComputeSphericalFixedDelta(numColumns, + // auto grid = std::make_unique(ComputeSpherical(numColumns, auto grid = std::make_unique(ComputeSphericalOnExtension(numColumns, numRows, originX, From 98c91f1f22cbb79b8c2ec1f6f96763125222c3f2 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Tue, 1 Sep 2026 08:45:24 +0200 Subject: [PATCH 7/7] GRIDEDIT-2303 Passing to other computer --- .../CurvilinearGridRectangular.cpp | 17 +++++----- .../src/CurvilinearGridRectangularTests.cpp | 34 +++++++++---------- 2 files changed, 25 insertions(+), 26 deletions(-) diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp index 7ba86b491..ab265e4a2 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp @@ -504,33 +504,32 @@ namespace meshkernel // At the equator, cos(0) = 1, so the conformal step matches the physical radian step. double d_step_lat_conformal = d_lat_rad; - // We center our unrotated grid around (0,0) space - // so that it spins perfectly around its true local center. // double half_width = (numRows * d_lon_rad) / 2.0; - - // Conformal latitude center tracker - // double origin_y_mercator = std::log(std::tan(M_PI / 4.0 + 0.0 / 2.0)); // Equator = 0 + double origin_y_mercator = std::log(std::tan(M_PI / 4.0)); for (int j = 0; j < numN; ++j) { // Local relative latitude centered around 0 (Equator) - double current_y_mercator = 0.0 + ((j - numColumns / 2.0) * d_step_lat_conformal); + // double current_y_mercator = 0.0 + ((j - numColumns / 2.0) * d_step_lat_conformal); + // double local_lat_rad = 2.0 * std::atan(std::exp(current_y_mercator)) - M_PI / 2.0; + + double current_y_mercator = origin_y_mercator + ((j - numColumns / 2.0) * d_step_lat_conformal); double local_lat_rad = 2.0 * std::atan(std::exp(current_y_mercator)) - M_PI / 2.0; + // double local_lon_rad = (i * d_lon_rad) - half_width; for (int i = 0; i < numM; ++i) { // Local relative longitude centered around 0 double local_lon_rad = (i - numRows / 2.0) * d_lon_rad; - // 1. Convert local relative point to 3D Cartesian + // Convert local relative point to 3D Cartesian double cos_local_lat = std::cos(local_lat_rad); Cartesian3DPoint p; p.x = cos_local_lat * std::cos(local_lon_rad); p.y = cos_local_lat * std::sin(local_lon_rad); p.z = std::sin(local_lat_rad); - // 2. Twist the grid locally by rotation_deg around its own center axis (1, 0, 0) - // Since the local center is at (0,0), its Cartesian vector is pointing down the X axis! + // Rotate the grid locally by rotation_deg around its origin if (rotation_deg != 0.0) { double ty = p.y * std::cos(rot_rad) - p.z * std::sin(rot_rad); diff --git a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp index f218582ad..a498fbf04 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp @@ -706,9 +706,9 @@ TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_Sho [[maybe_unused]] double blockSizeY = 0.02; [[maybe_unused]] double originX = 105.7; - [[maybe_unused]] double originY = 18.2; - [[maybe_unused]] double upperRightX = 106.79; - [[maybe_unused]] double upperRightY = 18.46; + [[maybe_unused]] double originY = 88.2; + [[maybe_unused]] double upperRightX = 106.79; + [[maybe_unused]] double upperRightY = 18.46; // [[maybe_unused]] double upperRightX = 13.3295; // [[maybe_unused]] double upperRightY = 23.5238; [[maybe_unused]] double angle = -43.0; @@ -723,35 +723,33 @@ TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_Sho // Execution CurvilinearGridRectangular curvilinearGridRectangular(Projection::spherical); - // auto mesh = curvilinearGridRectangular.Compute(45, 45, originX, originY, angle, blockSizeX, blockSizeY); - auto mesh = curvilinearGridRectangular.Compute(originX, originY, blockSizeX, blockSizeY, upperRightX, upperRightY, angle); + auto mesh = curvilinearGridRectangular.Compute(45, 45, originX, originY, angle, blockSizeX, blockSizeY); + // auto mesh = curvilinearGridRectangular.Compute(originX, originY, blockSizeX, blockSizeY, upperRightX, upperRightY, angle); auto nodes = mesh->ComputeNodes(); - meshkernel::Print (mesh->ComputeNodes (), mesh->ComputeEdges ()); + meshkernel::Print(mesh->ComputeNodes(), mesh->ComputeEdges()); - meshkernel::Mesh2D mesh2d (mesh->ComputeEdges (), mesh->ComputeNodes (), mesh->projection ()); - auto ortho = meshkernel::MeshOrthogonality::Compute (mesh2d); - std::ranges::sort (ortho); + meshkernel::Mesh2D mesh2d(mesh->ComputeEdges(), mesh->ComputeNodes(), mesh->projection()); + auto ortho = meshkernel::MeshOrthogonality::Compute(mesh2d); + std::ranges::sort(ortho); - for (size_t i = 0; i < ortho.size (); ++i) + for (size_t i = 0; i < ortho.size(); ++i) { - if (i % 10 == 0) std::cout << std::endl; - - std::cout << std::setw (15) << " " << ortho[i] << " "; - + if (i % 10 == 0) + std::cout << std::endl; + std::cout << std::setw(15) << " " << ortho[i] << " "; } - - double max = *std::max_element (ortho.begin (), ortho.end ()); + double max = *std::max_element(ortho.begin(), ortho.end()); std::cout << "Maximum orthoginality = " << max << std::endl; return; -#endif +#else double originX = 10.0; double originY = 10.0; @@ -788,4 +786,6 @@ TEST(CurvilinearGridUniform, CurvilinearGridRectangularOnExtension_WithAngle_Sho EXPECT_NEAR(nodes[i].x, expectedX[i], tolerance); EXPECT_NEAR(nodes[i].y, expectedY[i], tolerance); } + +#endif }