From ce074dc72af6b8316137efdcfc5b245267bfd4df Mon Sep 17 00:00:00 2001 From: Luca Carniato Date: Thu, 25 Jan 2024 16:08:25 +0100 Subject: [PATCH] Add mkernel_mesh2d_refine_ridges_based_on_gridded_samples API function (#287 | GRIDEDIT-873) --- .../BilinearInterpolationOnGriddedSamples.hpp | 8 +- .../include/MeshKernel/MeshRefinement.hpp | 2 +- .../tests/src/MeshRefinementTests.cpp | 145 ++++-------------- .../include/MeshKernelApi/MeshKernel.hpp | 19 +++ .../include/MeshKernelApi/Utils.hpp | 69 ++++++++- libs/MeshKernelApi/src/MeshKernel.cpp | 61 +++++++- .../tests/src/Mesh2DRefinmentTests.cpp | 80 +++++++++- tools/test_utils/CMakeLists.txt | 4 +- .../include/TestUtils/SampleFileReader.hpp | 1 - .../include/TestUtils/SampleGenerator.hpp | 47 ++++++ tools/test_utils/src/SampleGenerator.cpp | 85 ++++++++++ 11 files changed, 386 insertions(+), 135 deletions(-) create mode 100644 tools/test_utils/include/TestUtils/SampleGenerator.hpp create mode 100644 tools/test_utils/src/SampleGenerator.cpp diff --git a/libs/MeshKernel/include/MeshKernel/BilinearInterpolationOnGriddedSamples.hpp b/libs/MeshKernel/include/MeshKernel/BilinearInterpolationOnGriddedSamples.hpp index fe7e20777..bc8ecda9a 100644 --- a/libs/MeshKernel/include/MeshKernel/BilinearInterpolationOnGriddedSamples.hpp +++ b/libs/MeshKernel/include/MeshKernel/BilinearInterpolationOnGriddedSamples.hpp @@ -38,11 +38,7 @@ namespace meshkernel /// @brief Defines the iterpolatable data types. template - concept InterpolatableType = - std::same_as || - std::same_as || - std::same_as || - std::same_as; + concept InterpolatableType = std::same_as || std::same_as; /// @brief A class for performing bilinear interpolation on gridded samples template @@ -261,7 +257,7 @@ namespace meshkernel template double BilinearInterpolationOnGriddedSamples::GetGriddedValue(UInt columnIndex, UInt rowIndex) const { - const auto index = rowIndex * m_numXCoord + columnIndex; + const auto index = (m_numYCoord - rowIndex - 1) * m_numXCoord + columnIndex; return static_cast(m_values[index]); } diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index 807c925a5..f583eeb4c 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -81,6 +81,7 @@ namespace meshkernel LandWater = 3 }; + public: /// @brief Enumerator describing the different refinement types enum class RefinementType { @@ -89,7 +90,6 @@ namespace meshkernel RidgeDetection = 3 }; - public: /// @brief The constructor for refining based on samples /// @param[in] mesh The mesh to be refined /// @param[in] interpolant The averaging interpolation to use diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 77160a22d..170948dcf 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -3,17 +3,18 @@ #include -#include -#include -#include -#include -#include -#include -#include +#include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/MeshRefinement.hpp" +#include "MeshKernel/Parameters.hpp" +#include "MeshKernel/Polygons.hpp" +#include "TestUtils/Definitions.hpp" +#include "TestUtils/MakeMeshes.hpp" +#include "TestUtils/SampleFileReader.hpp" +#include "TestUtils/SampleGenerator.hpp" using namespace meshkernel; -TEST(MeshRefinement, FourByFourWithFourSamples) +TEST(MeshRefinement, MeshRefinementRefinementLevels_OnFourByFourWithFourSamples_ShouldRefinemesh) { auto mesh = MakeRectangularMeshForTesting(5, 5, 10.0, Projection::cartesian); @@ -211,7 +212,7 @@ TEST(MeshRefinement, RefinementOnAFourByFourMeshWithSamplesShouldRefine) ASSERT_EQ(10, mesh->m_edges[20].second); } -TEST(MeshRefinement, SmallTriangualMeshTwoSamples) +TEST(MeshRefinement, MeshRefinementRefinementLevels_SmallTriangualMeshTwoSamples_ShouldRefinemesh) { // Prepare auto mesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/SmallTriangularGrid_net.nc"); @@ -478,7 +479,7 @@ TEST(MeshRefinement, WindowOfRefinementFile) ASSERT_EQ(326, mesh->m_edges[915].second); } -TEST(MeshRefinement, WindowOfRefinementFileBasedOnLevels) +TEST(MeshRefinement, MeshRefinementRefinementLevels_OnWindowOfRefinementFile_ShouldRefinemesh) { // Prepare auto mesh = MakeRectangularMeshForTesting(4, 4, 40.0, Projection::cartesian, {197253.0, 442281.0}); @@ -700,7 +701,7 @@ TEST(MeshRefinement, FourByFourWithFourSamplesSpherical) ASSERT_EQ(6, mesh->m_edges[47].second); } -TEST(MeshRefinement, Refine_SphericalMesh_ShouldRefine) +TEST(MeshRefinement, RefinementFileBasedOnLevels_OnSphericalMesh_ShouldRefine) { // Prepare auto mesh = MakeRectangularMeshForTesting(6, 6, 0.0033, Projection::spherical, {41.1, 41.1}); @@ -819,9 +820,9 @@ TEST(MeshRefinement, BilinearInterpolationWithGriddedSamplesOnLandShouldNotRefin // Setup auto mesh = MakeRectangularMeshForTesting(2, 2, 10.0, Projection::cartesian); - std::vector values{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0}; + std::vector values{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0}; Point origin{-5.0, -5.0}; - auto interpolator = std::make_unique>(*mesh, 2, 2, origin, 10.0, values); + auto interpolator = std::make_unique>(*mesh, 2, 2, origin, 10.0, values); MeshRefinementParameters meshRefinementParameters; meshRefinementParameters.max_num_refinement_iterations = 1; @@ -846,9 +847,9 @@ TEST(MeshRefinement, BilinearInterpolationWithGriddedSamplesOnLandAndSeaShouldRe // Setup auto mesh = MakeRectangularMeshForTesting(2, 2, 10.0, Projection::cartesian); - std::vector values{-1.0, -2.0, 3.0, -4.0, -5.0, 6.0, 7.0, 8.0, 9.0}; + std::vector values{-1.0, -2.0, 3.0, -4.0, -5.0, 6.0, 7.0, 8.0, 9.0}; Point origin{-5.0, -5.0}; - auto interpolator = std::make_unique>(*mesh, 3, 3, origin, 10.0, values); + auto interpolator = std::make_unique>(*mesh, 3, 3, origin, 10.0, values); MeshRefinementParameters meshRefinementParameters; meshRefinementParameters.max_num_refinement_iterations = 1; @@ -873,9 +874,9 @@ TEST(MeshRefinement, BilinearInterpolationWithAllGriddedSamplesOnSeaShouldRefine // Setup auto mesh = MakeRectangularMeshForTesting(2, 2, 10.0, Projection::cartesian); - std::vector values{-1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -9.0}; + std::vector values{-1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -9.0}; Point origin{-5.0, -5.0}; - auto interpolator = std::make_unique>(*mesh, 2, 2, origin, 10.0, values); + auto interpolator = std::make_unique>(*mesh, 2, 2, origin, 10.0, values); MeshRefinementParameters meshRefinementParameters; meshRefinementParameters.max_num_refinement_iterations = 1; @@ -897,107 +898,16 @@ TEST(MeshRefinement, BilinearInterpolationWithAllGriddedSamplesOnSeaShouldRefine ASSERT_EQ(4, mesh->GetNumEdges()); } -enum class RidgeRefinementTestCase -{ - GaussianBump = 1, - GaussianWave = 2, - RidgeXDirection = 3, - ArctanFunction = 4, -}; - -std::vector generateSampleData(RidgeRefinementTestCase testcase, - UInt nx = 10, - UInt ny = 10, - double deltaX = 10.0, - double deltaY = 10.0) -{ - UInt start = 0; - UInt size = (nx - start) * (ny - start); - std::vector sampleData(size); - - std::vector sampleDataMatrix(ny, std::vector(nx)); - - const double centreX = static_cast((nx - 1) / 2) * deltaX; - const double centreY = static_cast((ny - 1) / 2) * deltaY; - - const double scale = ny / 4.0 * deltaY; - - const double r = nx / 5 * deltaX; - const double maxx = (nx - 1) * deltaX; - - std::function generateSample; - switch (testcase) - { - case RidgeRefinementTestCase::GaussianBump: - generateSample = [&](double x, double y) - { - const double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY); - return 100.0 * std::exp(-0.025 * centre); - }; - break; - case RidgeRefinementTestCase::GaussianWave: - generateSample = [&](double x, double y) - { - const double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY); - const double factor = std::max(1e-6, std::exp(-0.00025 * centre)); - return 100.0 * factor; - }; - break; - case RidgeRefinementTestCase::RidgeXDirection: - generateSample = [&](double x, double y) - { - const double sinx = std::sin(x / maxx * M_PI * 4.0); - const double xxx = scale * sinx + centreY; - return 10 * (std::atan(20.0 * (xxx - y)) + M_PI / 2.0); - }; - break; - case RidgeRefinementTestCase::ArctanFunction: - generateSample = [&](double x, double y) - { - const double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY); - return 10 * (std::atan(20.0 * (r * r - centre)) + M_PI / 2.0); - }; - break; - default: - throw std::invalid_argument("invalid ridge refinement test case"); - } - - for (int i = ny - 1; i >= 0; --i) - { - - for (UInt j = start; j < nx; ++j) - { - const double y = deltaY * i; - const double x = deltaX * j; - sampleDataMatrix[ny - 1 - i][j] = generateSample(x, y); - } - } - - UInt count = 0; - for (UInt j = start; j < nx; ++j) - { - for (int i = ny - 1; i >= 0; --i) - { - const double y = deltaY * i; - const double x = deltaX * j; - sampleData[count] = {x, y, sampleDataMatrix[ny - 1 - i][j]}; - count++; - } - } - - return sampleData; -} - -class RidgeRefinementTestCases : public testing::TestWithParam> +class RidgeRefinementTestCases : public testing::TestWithParam> { public: - [[nodiscard]] static std::vector> GetData() + [[nodiscard]] static std::vector> GetData() { return std::vector{ - std::make_tuple(RidgeRefinementTestCase::GaussianBump, 1165, 2344), - std::make_tuple(RidgeRefinementTestCase::GaussianWave, 5297, 10784), - std::make_tuple(RidgeRefinementTestCase::RidgeXDirection, 2618, 5694), - std::make_tuple(RidgeRefinementTestCase::ArctanFunction, 2309, 5028)}; + std::make_tuple(FunctionTestCase::GaussianBump, 1165, 2344), + std::make_tuple(FunctionTestCase::GaussianWave, 5297, 10784), + std::make_tuple(FunctionTestCase::RidgeXDirection, 2618, 5694), + std::make_tuple(FunctionTestCase::ArctanFunction, 2309, 5028)}; } }; @@ -1015,10 +925,9 @@ TEST_P(RidgeRefinementTestCases, expectedResults) double dimX = (nx - 1) * deltaX; double dimY = (ny - 1) * deltaY; - std::shared_ptr mesh = MakeRectangularMeshForTesting(nx, ny, dimX, dimY, Projection::cartesian); + auto mesh = MakeRectangularMeshForTesting(nx, ny, dimX, dimY, Projection::cartesian); UInt superSample = 2; - UInt sampleNx = (nx - 1) * superSample + 1; UInt sampleNy = (ny - 1) * superSample + 1; @@ -1027,10 +936,10 @@ TEST_P(RidgeRefinementTestCases, expectedResults) const auto sampleData = generateSampleData(testCase, sampleNx, sampleNy, sampleDeltaX, sampleDeltaY); - auto samples = SamplesHessianCalculator::ComputeSamplesHessian(sampleData, mesh->m_projection, 0, sampleNx, sampleNy); + auto samplesHessian = SamplesHessianCalculator::ComputeSamplesHessian(sampleData, mesh->m_projection, 0, sampleNx, sampleNy); auto interpolator = std::make_unique(*mesh, - samples, + samplesHessian, AveragingInterpolation::Method::Max, Location::Faces, 1.0, diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 74b95fb78..8dee5270a 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -1231,6 +1231,25 @@ namespace meshkernelapi int minimumNumSamples, const meshkernel::MeshRefinementParameters& meshRefinementParameters); + /// @brief Refines a mesh2d based on samples with ridge refinement. This method automatically detects the ridges in a sample set. + /// + /// The number of successive splits is indicated on the sample value. + /// For example a value of 0 means no split and hence no refinement, a value of 1 a single split (a quadrilateral face generates 4 faces), + /// a value of 2 two splits (a quadrilateral face generates 16 faces). + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] griddedSamples The gridded samples + /// @param[in] relativeSearchRadius The relative search radius relative to the face size, used for some interpolation algorithms + /// @param[in] minimumNumSamples The minimum number of samples used for some averaging algorithms + /// @param[in] numberOfSmoothingIterations The number of smoothing iterations to apply to the input sample set + /// @param[in] meshRefinementParameters The mesh refinement parameters + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_refine_ridges_based_on_gridded_samples(int meshKernelId, + const GriddedSamples& griddedSamples, + double relativeSearchRadius, + int minimumNumSamples, + int numberOfSmoothingIterations, + const meshkernel::MeshRefinementParameters& meshRefinementParameters); + /// @brief Remove any disconnected regions from a mesh2d. /// /// The assumption is that the main region of interest has the largest number of elements. diff --git a/libs/MeshKernelApi/include/MeshKernelApi/Utils.hpp b/libs/MeshKernelApi/include/MeshKernelApi/Utils.hpp index ae862b9e0..0ab5ed660 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/Utils.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/Utils.hpp @@ -36,6 +36,7 @@ #include #include +#include #include #include @@ -217,6 +218,72 @@ namespace meshkernelapi } } + /// @brief Computes the samples represented in gridded data in a vector of samples + /// @param[in] griddedSamples The gridded data to convert + /// @returns The converted vector of samples + template + static std::vector ComputeGriddedDataSamples(const GriddedSamples& griddedSamples) + { + std::vector result; + meshkernel::Point origin{griddedSamples.x_origin, griddedSamples.y_origin}; + const auto numSamples = static_cast(griddedSamples.num_x * griddedSamples.num_y); + result.resize(numSamples); + const T* valuePtr = static_cast(griddedSamples.values); + if (griddedSamples.x_coordinates == nullptr || griddedSamples.y_coordinates == nullptr) + { + meshkernel::UInt index = 0; + + for (int j = 0; j < griddedSamples.num_x; ++j) + { + for (int i = griddedSamples.num_y - 1; i >= 0; --i) + { + const auto griddedIndex = griddedSamples.num_x * i + j; + result[index].x = origin.x + j * griddedSamples.cell_size; + result[index].y = origin.y + i * griddedSamples.cell_size; + result[index].value = static_cast(valuePtr[griddedIndex]); + index++; + } + } + return result; + } + + meshkernel::UInt index = 0; + for (int j = 0; j < griddedSamples.num_x; ++j) + { + for (int i = griddedSamples.num_y - 1; i >= 0; --i) + { + const auto griddedIndex = griddedSamples.num_x * i + j; + result[index].x = origin.x + griddedSamples.x_coordinates[griddedIndex]; + result[index].y = origin.y + griddedSamples.y_coordinates[griddedIndex]; + result[index].value = static_cast(valuePtr[griddedIndex]); + index++; + } + } + return result; + } + + /// @brief Converts the samples represented in gridded data in a vector of samples + /// @param[in] griddedSamples The gridded data to convert + /// @returns The converted vector of samples + static std::vector ConvertGriddedData(const GriddedSamples& griddedSamples) + { + std::vector result; + if (griddedSamples.num_x <= 0 || griddedSamples.num_y <= 0) + { + return result; + } + + if (griddedSamples.value_type == static_cast(meshkernel::InterpolationValues::shortType)) + { + return ComputeGriddedDataSamples(griddedSamples); + } + if (griddedSamples.value_type == static_cast(meshkernel::InterpolationValues::floatType)) + { + return ComputeGriddedDataSamples(griddedSamples); + } + throw meshkernel::MeshKernelError("The value type for the gridded data samples is invalid."); + } + /// @brief Sets splines from a geometry list /// @param[in] geometryListIn The input geometry list /// @param[out] spline The spline which will be set @@ -405,7 +472,7 @@ namespace meshkernelapi const GriddedSamples& griddedSamples) { meshkernel::Point origin{griddedSamples.x_origin, griddedSamples.y_origin}; - if (griddedSamples.x_coordinates == nullptr && griddedSamples.y_coordinates == nullptr) + if (griddedSamples.x_coordinates == nullptr || griddedSamples.y_coordinates == nullptr) { return std::make_unique>(mesh2d, griddedSamples.num_x, diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index 0b5ffd5a5..19ce5dc21 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -27,6 +27,7 @@ #include "MeshKernel/CurvilinearGrid/CurvilinearGridDeleteExterior.hpp" #include "MeshKernel/Mesh2DIntersections.hpp" +#include "MeshKernel/SamplesHessianCalculator.hpp" #include "MeshKernel/CurvilinearGrid/CurvilinearGridDeleteInterior.hpp" #include @@ -1737,14 +1738,18 @@ namespace meshkernelapi // averagingMethod may be used uninitialised; meshkernel::AveragingInterpolation::Method averagingMethod; - if (meshRefinementParameters.refinement_type == 1) + if (meshRefinementParameters.refinement_type == static_cast(meshkernel::MeshRefinement::RefinementType::WaveCourant)) { averagingMethod = meshkernel::AveragingInterpolation::Method::MinAbsValue; } - if (meshRefinementParameters.refinement_type == 2) + else if (meshRefinementParameters.refinement_type == static_cast(meshkernel::MeshRefinement::RefinementType::RefinementLevels)) { averagingMethod = meshkernel::AveragingInterpolation::Method::Max; } + else + { + throw meshkernel::MeshKernelError("Invalid mesh refinement type."); + } const bool refineOutsideFace = meshRefinementParameters.account_for_samples_outside == 1 ? true : false; const bool transformSamples = meshRefinementParameters.refinement_type == 2 ? true : false; @@ -1770,6 +1775,58 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_mesh2d_refine_ridges_based_on_gridded_samples(int meshKernelId, + const GriddedSamples& samples, + double relativeSearchRadius, + int minimumNumSamples, + int numberOfSmoothingIterations, + const meshkernel::MeshRefinementParameters& meshRefinementParameters) + { + lastExitCode = meshkernel::ExitCode::Success; + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + if (meshKernelState[meshKernelId].m_mesh2d->GetNumNodes() <= 0) + { + throw meshkernel::ConstraintError("The selected mesh has no nodes."); + } + if (meshRefinementParameters.refinement_type != static_cast(meshkernel::MeshRefinement::RefinementType::RidgeDetection)) + { + throw meshkernel::MeshKernelError("The mesh refinement type in MeshRefinementParameters must be set equal to ridge refinement."); + } + + const auto samplesVector = ConvertGriddedData(samples); + + auto samplesHessian = meshkernel::SamplesHessianCalculator::ComputeSamplesHessian(samplesVector, + meshKernelState[meshKernelId].m_projection, + numberOfSmoothingIterations, + samples.num_x, + samples.num_y); + + auto averaging = std::make_unique(*meshKernelState[meshKernelId].m_mesh2d, + samplesHessian, + meshkernel::AveragingInterpolation::Method::Max, + meshkernel::Location::Faces, + relativeSearchRadius, + false, + false, + static_cast(minimumNumSamples)); + + meshkernel::MeshRefinement meshRefinement(*meshKernelState[meshKernelId].m_mesh2d, + std::move(averaging), + meshRefinementParameters); + meshRefinement.Compute(); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + MKERNEL_API int mkernel_mesh2d_refine_based_on_gridded_samples(int meshKernelId, const GriddedSamples& griddedSamples, const meshkernel::MeshRefinementParameters& meshRefinementParameters, diff --git a/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp b/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp index 5dc7982f7..e1f65f878 100644 --- a/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp +++ b/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp @@ -13,6 +13,11 @@ #include #include "CartesianApiTestFixture.hpp" +#include "MeshKernel/AveragingInterpolation.hpp" +#include "MeshKernel/MeshRefinement.hpp" +#include "MeshKernel/SamplesHessianCalculator.hpp" +#include "SampleFileWriter.hpp" +#include "SampleGenerator.hpp" TEST_F(CartesianApiTestFixture, RefineAPolygonThroughApi) { @@ -56,7 +61,7 @@ TEST_F(CartesianApiTestFixture, RefineAPolygonThroughApi) ASSERT_NEAR(92.626556, geometryListOut.coordinates_y[0], tolerance); } -TEST_F(CartesianApiTestFixture, RefineBasedOnSamples_OnAUniformMesh_shouldRefineMesh) +TEST_F(CartesianApiTestFixture, RefineBasedOnSamplesWaveCourant_OnAUniformMesh_shouldRefineMesh) { // Prepare MakeMesh(10, 10, 25); @@ -180,7 +185,7 @@ TEST_F(CartesianApiTestFixture, RefineAGridBasedOnPolygonThroughApi) ASSERT_EQ(595, mesh2d.num_edges); } -TEST(MeshRefinement, Mesh2DRefineBasedOnGriddedSamples_WithGriddedSamples_ShouldRefineMesh) +TEST(MeshRefinement, Mesh2DRefineBasedOnGriddedSamplesWaveCourant_WithGriddedSamples_ShouldRefineMesh) { // Prepare int meshKernelId; @@ -236,7 +241,7 @@ TEST(MeshRefinement, Mesh2DRefineBasedOnGriddedSamples_WithGriddedSamples_Should ASSERT_EQ(1936, mesh2dResults.num_face_nodes); } -TEST_F(CartesianApiTestFixture, Mesh2DRefineBasedOnGriddedSamples_WithNotUniformlySpacedSamples_ShouldRefineMesh) +TEST_F(CartesianApiTestFixture, Mesh2DRefineBasedOnGriddedSamplesWaveCourant_WithNotUniformlySpacedSamples_ShouldRefineMesh) { // Prepare meshkernel::UInt nRows{5}; @@ -301,7 +306,7 @@ TEST_F(CartesianApiTestFixture, Mesh2DRefineBasedOnGriddedSamples_WithNotUniform ASSERT_EQ(76, mesh2dResults.num_faces); } -TEST(MeshRefinement, RefineBasedOnGriddedSamples_WithUniformSamplesAndSphericalCoordinates_ShouldRefineMesh2d) +TEST(MeshRefinement, RefineBasedOnGriddedSamplesWaveCourant_WithUniformSamplesAndSphericalCoordinates_ShouldRefineMesh2d) { // Prepare int meshKernelId; @@ -354,7 +359,7 @@ TEST(MeshRefinement, RefineBasedOnGriddedSamples_WithUniformSamplesAndSphericalC ASSERT_EQ(21212, mesh2dResults.num_face_nodes); } -TEST(MeshRefinement, RefineBasedOnGriddedSamples_WithUniformSamplesAndSphericalCoordinatesAndLargeMinEdgeSize_ShouldNotRefineMesh2d) +TEST(MeshRefinement, RefineBasedOnGriddedSamplesWaveCourant_WithUniformSamplesAndSphericalCoordinatesAndLargeMinEdgeSize_ShouldNotRefineMesh2d) { // Prepare int meshKernelId; @@ -613,3 +618,68 @@ TEST_P(MeshRefinementSampleValueTypes, parameters) } INSTANTIATE_TEST_SUITE_P(MeshRefinement, MeshRefinementSampleValueTypes, ::testing::ValuesIn(MeshRefinementSampleValueTypes::GetData())); + +TEST_F(CartesianApiTestFixture, RefineAMeshBasedOnRidgeRefinement_OnAUniformMesh_shouldRefineMesh) +{ + // Prepare + meshkernel::UInt nRows{21}; + meshkernel::UInt nCols{41}; + const double meshDelta = 10.0; + MakeMesh(nRows, nCols, 10.0); + + // Generate gridded samples, from a gaussian distribution + const int numSamplesXCoordinates = (nCols - 1) * 2 + 1; + const int numSamplesYCoordinates = (nRows * 1) * 2 + 1; + const double deltaX = 5.0; + const double deltaY = 5.0; + const auto sampleData = generateSampleData(FunctionTestCase::GaussianBump, numSamplesXCoordinates, numSamplesYCoordinates, deltaX, deltaY); + + // Create an instance of gridded samples, + meshkernelapi::GriddedSamples griddedSamples; + griddedSamples.num_x = numSamplesXCoordinates; + griddedSamples.num_y = numSamplesYCoordinates; + griddedSamples.x_origin = 0.0; + griddedSamples.y_origin = 0.0; + griddedSamples.cell_size = meshDelta * 0.1; + + // Flatten the values before passing them to Mesh Kernel + meshkernel::UInt index = 0; + std::vector values(numSamplesXCoordinates * numSamplesYCoordinates, 0.0); + for (int j = 0; j < numSamplesYCoordinates; ++j) + { + for (int i = 0; i < numSamplesXCoordinates; ++i) + { + const auto griddedIndex = griddedSamples.num_y * i + j; + values[index] = static_cast(sampleData[griddedIndex].value); + index++; + } + } + int interpolationType; + auto errorCode = meshkernelapi::mkernel_get_interpolation_type_float(interpolationType); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + griddedSamples.values = values.data(); + griddedSamples.value_type = interpolationType; + + // Create meshrefinement parameters + meshkernel::MeshRefinementParameters meshRefinementParameters; + meshRefinementParameters.max_num_refinement_iterations = 3; + meshRefinementParameters.refine_intersected = 0; + meshRefinementParameters.use_mass_center_when_refining = 0; + meshRefinementParameters.min_edge_size = 2.0; + meshRefinementParameters.account_for_samples_outside = 0; + meshRefinementParameters.connect_hanging_nodes = 1; + meshRefinementParameters.smoothing_iterations = 0; + meshRefinementParameters.refinement_type = 3; + + // Execute + const auto meshKernelId = GetMeshKernelId(); + mkernel_mesh2d_refine_ridges_based_on_gridded_samples(meshKernelId, griddedSamples, 1.01, 1, 0, meshRefinementParameters); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Assert on the mesh values + meshkernelapi::Mesh2D mesh2d{}; + errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + ASSERT_EQ(956, mesh2d.num_nodes); + ASSERT_EQ(1872, mesh2d.num_edges); +} diff --git a/tools/test_utils/CMakeLists.txt b/tools/test_utils/CMakeLists.txt index 44606e63b..b4e1511c0 100644 --- a/tools/test_utils/CMakeLists.txt +++ b/tools/test_utils/CMakeLists.txt @@ -32,16 +32,18 @@ set( ${SRC_DIR}/MakeMeshes.cpp ${SRC_DIR}/SampleFileReader.cpp ${SRC_DIR}/SampleFileWriter.cpp + ${SRC_DIR}/SampleGenerator.cpp ) # list of target headers set( INC_LIST ${DOMAIN_INC_DIR}/Definitions.hpp - ${DOMAIN_INC_DIR}/MakeCurvilinearGrids.hpp + ${DOMAIN_INC_DIR}/MakeCurvilinearGrids.hpp ${DOMAIN_INC_DIR}/MakeMeshes.hpp ${DOMAIN_INC_DIR}/SampleFileReader.hpp ${DOMAIN_INC_DIR}/SampleFileWriter.hpp + ${DOMAIN_INC_DIR}/SampleGenerator.hpp ) # add sources to target diff --git a/tools/test_utils/include/TestUtils/SampleFileReader.hpp b/tools/test_utils/include/TestUtils/SampleFileReader.hpp index e31d6477d..5ec23e389 100644 --- a/tools/test_utils/include/TestUtils/SampleFileReader.hpp +++ b/tools/test_utils/include/TestUtils/SampleFileReader.hpp @@ -109,7 +109,6 @@ std::tuple> R numlines++; } - std::reverse(rows.begin(), rows.end()); std::vector values; for (size_t i = 0; i < rows.size(); ++i) { diff --git a/tools/test_utils/include/TestUtils/SampleGenerator.hpp b/tools/test_utils/include/TestUtils/SampleGenerator.hpp new file mode 100644 index 000000000..eda55ac83 --- /dev/null +++ b/tools/test_utils/include/TestUtils/SampleGenerator.hpp @@ -0,0 +1,47 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2021. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation version 3. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include + +#include "MeshKernel/Entities.hpp" + +enum class FunctionTestCase +{ + GaussianBump = 1, + GaussianWave = 2, + RidgeXDirection = 3, + ArctanFunction = 4, +}; + +std::vector +generateSampleData(FunctionTestCase testcase, + meshkernel::UInt nx = 10, + meshkernel::UInt ny = 10, + double deltaX = 10.0, + double deltaY = 10.0); diff --git a/tools/test_utils/src/SampleGenerator.cpp b/tools/test_utils/src/SampleGenerator.cpp new file mode 100644 index 000000000..d6e8f8959 --- /dev/null +++ b/tools/test_utils/src/SampleGenerator.cpp @@ -0,0 +1,85 @@ +#include "SampleGenerator.hpp" + +std::vector +generateSampleData(FunctionTestCase testcase, + meshkernel::UInt nx, + meshkernel::UInt ny, + double deltaX, + double deltaY) +{ + meshkernel::UInt start = 0; + meshkernel::UInt size = (nx - start) * (ny - start); + std::vector sampleData(size); + + std::vector sampleDataMatrix(ny, std::vector(nx)); + + const double centreX = static_cast((nx - 1) / 2) * deltaX; + const double centreY = static_cast((ny - 1) / 2) * deltaY; + + const double scale = ny / 4.0 * deltaY; + + const double r = nx / 5 * deltaX; + const double maxx = (nx - 1) * deltaX; + + std::function generateSample; + switch (testcase) + { + case FunctionTestCase::GaussianBump: + generateSample = [&](double x, double y) + { + const double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY); + return 100.0 * std::exp(-0.025 * centre); + }; + break; + case FunctionTestCase::GaussianWave: + generateSample = [&](double x, double y) + { + const double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY); + const double factor = std::max(1e-6, std::exp(-0.00025 * centre)); + return 100.0 * factor; + }; + break; + case FunctionTestCase::RidgeXDirection: + generateSample = [&](double x, double y) + { + const double sinx = std::sin(x / maxx * M_PI * 4.0); + const double xxx = scale * sinx + centreY; + return 10 * (std::atan(20.0 * (xxx - y)) + M_PI / 2.0); + }; + break; + case FunctionTestCase::ArctanFunction: + generateSample = [&](double x, double y) + { + const double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY); + return 10 * (std::atan(20.0 * (r * r - centre)) + M_PI / 2.0); + }; + break; + default: + throw std::invalid_argument("invalid ridge refinement test case"); + } + + for (int i = ny - 1; i >= 0; --i) + { + + for (meshkernel::UInt j = start; j < nx; ++j) + { + const double y = deltaY * i; + const double x = deltaX * j; + sampleDataMatrix[ny - 1 - i][j] = generateSample(x, y); + } + } + + meshkernel::UInt count = 0; + for (meshkernel::UInt j = start; j < nx; ++j) + { + for (int i = ny - 1; i >= 0; --i) + { + const double y = deltaY * i; + const double x = deltaX * j; + sampleData[count] = {x, y, sampleDataMatrix[ny - 1 - i][j]}; + count++; + } + } + + return sampleData; +}