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 diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp index a3bae2d93..ce3014599 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp @@ -27,6 +27,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 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 /// @param[in] numRows The number of columns in y direction @@ -120,13 +133,44 @@ 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; + + // 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, + const double originY, + const double angle, + 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 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 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 9ae97c052..ab265e4a2 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridRectangular.cpp @@ -62,6 +62,9 @@ namespace meshkernel if (m_projection == Projection::spherical) { + // 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, @@ -124,8 +127,9 @@ 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, @@ -332,8 +336,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,4 +373,256 @@ namespace meshkernel throw NotImplementedError("Projection value {} not supported", static_cast(m_projection)); } + Point CurvilinearGridRectangular::RotateByAngle(const double originX, const double originY, + const double upperRightX, const double upperRightY, + const double cosAngle, + const double sinAngle) const + { + + if (m_projection == Projection::cartesian) + { + 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}); + + // 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}; + + // k \cdot v + double dot = k.x * point3d.x + k.y * point3d.y + k.z * point3d.z; + + // k \cross v + Cartesian3DPoint crossProd = VectorProduct(k, point3d); + + // 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)}; + + Point pointOnSphere = Cartesian3DToSpherical(rotatedPoint3d, upperRightX); + + return pointOnSphere; + } + } + + 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, + originX, + originY, + 0.0 * 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::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::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; + + // double half_width = (numRows * d_lon_rad) / 2.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 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; + + // 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); + + // 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); + 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, + 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"); + + 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"); + } + + if (numRows <= 0) + { + throw AlgorithmError("Number of rows cannot be <= 0"); + } + + 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, + 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/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 b501457ca..a498fbf04 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridRectangularTests.cpp @@ -43,6 +43,8 @@ #include #include +#include + using namespace meshkernel; TEST(CurvilinearGridUniform, CurvilinearGridRectangular_WithPolygon_ShouldComputeCurvilinearGrid) @@ -694,3 +696,96 @@ TEST_P(CurvilinearGridUniformTests, parameters) ASSERT_EQ(numValidNodes, expectedNumNodes); } INSTANTIATE_TEST_SUITE_P(curvilinearGridDeletionTests, CurvilinearGridUniformTests, ::testing::ValuesIn(CurvilinearGridUniformTests::GetData())); + +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 = 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; + + // [[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; + +#else + + double originX = 10.0; + double originY = 10.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(); + + 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}; + + ASSERT_EQ(nodes.size(), expectedX.size()); + + 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); + } + +#endif +} 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