From 72d33746a5a54e094c535a849b3e07a75391147e Mon Sep 17 00:00:00 2001 From: Fuad Hasibul Hasan Date: Fri, 17 Jan 2025 19:37:03 -0500 Subject: [PATCH 1/4] interpolator base class, tests [breaking] the test is still not passing. pushed it to save it at current state --- src/pcms/interpolator/CMakeLists.txt | 2 + src/pcms/interpolator/interpolation_base.cpp | 207 +++++++++++++++++++ src/pcms/interpolator/interpolation_base.h | 73 +++++++ test/CMakeLists.txt | 1 + test/test_interpolation_class.cpp | 73 +++++++ 5 files changed, 356 insertions(+) create mode 100644 src/pcms/interpolator/interpolation_base.cpp create mode 100644 src/pcms/interpolator/interpolation_base.h create mode 100644 test/test_interpolation_class.cpp diff --git a/src/pcms/interpolator/CMakeLists.txt b/src/pcms/interpolator/CMakeLists.txt index 1e17e6c5..b60df52c 100644 --- a/src/pcms/interpolator/CMakeLists.txt +++ b/src/pcms/interpolator/CMakeLists.txt @@ -12,10 +12,12 @@ set(PCMS_FIELD_TRANSFER_HEADERS linear_interpolant.hpp multidimarray.hpp MLS_rbf_options.hpp + interpolation_base.h ) set(PCMS_FIELD_TRANSFER_SOURCES MLS_rbf_options.cpp + interpolation_base.cpp ) install(FILES ${PCMS_FIELD_TRANSFER_HEADERS} DESTINATION include/pcms/interpolator) diff --git a/src/pcms/interpolator/interpolation_base.cpp b/src/pcms/interpolator/interpolation_base.cpp new file mode 100644 index 00000000..f1e20be1 --- /dev/null +++ b/src/pcms/interpolator/interpolation_base.cpp @@ -0,0 +1,207 @@ +// +// Created by hasanm4 on 1/17/25. +// + +#include "interpolation_base.h" +#include "MLS_rbf_options.hpp" + +#include + +Omega_h::Reals getCentroids(Omega_h::Mesh& mesh) +{ + OMEGA_H_CHECK_PRINTF( + mesh.dim() == 2, "Only 2D meshes are supported but found %d\n", mesh.dim()); + + const auto& coords = mesh.coords(); + Omega_h::Write centroids(mesh.nfaces() * 2, 0.0); + + auto face2node = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + Omega_h::parallel_for( + mesh.nfaces(), OMEGA_H_LAMBDA(Omega_h::LO face) { + auto nodes = Omega_h::gather_verts<3>(face2node, face); + Omega_h::Few, 3> face_coords = + Omega_h::gather_vectors<3, 2>(coords, nodes); + Omega_h::Vector<2> centroid = Omega_h::average(face_coords); + centroids[2 * face + 0] = centroid[0]; + centroids[2 * face + 1] = centroid[1]; + }); + + return {centroids}; +} + +MLSInterpolationHandler::MLSInterpolationHandler( + const std::string& source_mesh_fname, double radius, bool adapt_radius) + : radius_(radius), adapt_radius_(adapt_radius) +{ + single_mesh_ = true; + library_ = Omega_h::Library(nullptr, nullptr); + + source_mesh_ = Omega_h::binary::read(source_mesh_fname, &library_); + source_coords_ = source_mesh_.coords(); + target_coords_ = getCentroids(source_mesh_); + + OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2, + "Only 2D meshes are supported but found %d\n", + source_mesh_.dim()); + + source_field_ = + Omega_h::HostWrite(source_mesh_.nverts(), "source field"); + + // used nfaces to interpolate to the centroid of the faces + target_field_ = + Omega_h::HostWrite(source_mesh_.nfaces(), "target field"); + + find_supports(); +} + +MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh source_mesh, + double radius, + bool adapt_radius) + : radius_(radius), adapt_radius_(adapt_radius) +{ + single_mesh_ = true; + source_mesh_ = std::move(source_mesh); + source_coords_ = source_mesh_.coords(); + target_coords_ = getCentroids(source_mesh_); + + OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2, + "Only 2D meshes are supported but found %d\n", + source_mesh_.dim()); + + source_field_ = + Omega_h::HostWrite(source_mesh_.nverts(), "source field"); + + target_field_ = + Omega_h::HostWrite(source_mesh_.nfaces(), "target field"); +} + +// TODO: delegate this constructor to a single one +MLSInterpolationHandler::MLSInterpolationHandler( + const std::string& source_mesh_fname, const std::string& target_mesh_fname, + const double radius, const bool adapt_radius) + : radius_(radius), adapt_radius_(adapt_radius) +{ + // TODO: take argc and argv + library_ = Omega_h::Library(nullptr, nullptr); + + source_mesh_ = Omega_h::binary::read(source_mesh_fname, &library_); + target_mesh_ = Omega_h::binary::read(target_mesh_fname, &library_); + + OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2 && target_mesh_.dim() == 2, + "Only 2D meshes are supported but found %d, %d\n", + source_mesh_.dim(), target_mesh_.dim()); + + source_coords_ = source_mesh_.coords(); + target_coords_ = target_mesh_.coords(); + + source_field_ = + Omega_h::HostWrite(source_mesh_.nverts(), "source field"); + target_field_ = + Omega_h::HostWrite(target_mesh_.nverts(), "target field"); + + find_supports(); +} + +MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh source_mesh, + Omega_h::Mesh target_mesh, + const double radius, + const bool adapt_radius) + : radius_(radius), adapt_radius_(adapt_radius) +{ + // TODO: check if move works for mesh + source_mesh_ = std::move(source_mesh); + target_mesh_ = std::move(target_mesh); + + OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2 && target_mesh_.dim() == 2, + "Only 2D meshes are supported but found %d, %d\n", + source_mesh_.dim(), target_mesh_.dim()); + + source_coords_ = source_mesh_.coords(); + target_coords_ = target_mesh_.coords(); + + source_field_ = + Omega_h::HostWrite(source_mesh_.nverts(), "source field"); + target_field_ = + Omega_h::HostWrite(target_mesh_.nverts(), "target field"); + + find_supports(); +} + +// TODO : find way to avoid this copy +void copyHostScalarArrayView2HostWrite( + const pcms::ScalarArrayView& source, + Omega_h::HostWrite& target) +{ + OMEGA_H_CHECK_PRINTF( + source.size() == target.size(), + "Size mismatch in copy_data_from_ScalarArray_to_HostWrite: %zu %d\n", + source.size(), target.size()); + + for (int i = 0; i < source.size(); ++i) { + target[i] = source[i]; + } +} +void copyHostWrite2ScalarArrayView( + const Omega_h::HostWrite& source, + pcms::ScalarArrayView& target) +{ + OMEGA_H_CHECK_PRINTF( + source.size() == target.size(), + "Size mismatch in copy_data_from_HostWrite_to_ScalarArray: %d %zu\n", + source.size(), target.size()); + + for (int i = 0; i < source.size(); ++i) { + target[i] = source[i]; + } +} + +void MLSInterpolationHandler::eval( + const pcms::ScalarArrayView& source_field, + pcms::ScalarArrayView& target_field) +{ + OMEGA_H_CHECK_PRINTF(source_field.size() == + source_coords_.size() / source_mesh_.dim(), + "Source Data and Source Points size mismatch: %zu %d\n", + source_field.size(), source_coords_.size()); + if (single_mesh_) { + OMEGA_H_CHECK_PRINTF( + target_field.size() == source_mesh_.nfaces(), + "Target Data and Target Points size mismatch: %zu %d\n", + target_field.size(), target_coords_.size()); + } else { + OMEGA_H_CHECK_PRINTF( + target_field.size() == target_coords_.size() / target_mesh_.dim(), + "Target Data and Target Points size mismatch: %zu %d\n", + target_field.size(), target_coords_.size()); + } + + copyHostScalarArrayView2HostWrite(source_field, source_field_); + + auto target_field_write = mls_interpolation( + Omega_h::Reals(source_field_), source_coords_, target_coords_, supports_, 2, + 4, supports_.radii2, RadialBasisFunction::RBF_GAUSSIAN); + + target_field_ = Omega_h::HostWrite(target_field_write); + copyHostWrite2ScalarArrayView(target_field_, target_field); +} + +void MLSInterpolationHandler::find_supports( + const int min_req_support) // ? Needed? +{ + if (single_mesh_) { + supports_ = + searchNeighbors(source_mesh_, radius_, min_req_support, adapt_radius_); + } else { // not single mesh + supports_ = searchNeighbors(source_mesh_, target_mesh_, radius_, + min_req_support, adapt_radius_); + } + + OMEGA_H_CHECK_PRINTF( + supports_.radii2[0] < 1e-10, + "Radius squared has to be more than zero found found [0] = %f\n", + supports_.radii2[0]); + + for (int i = 0; i < supports_.radii2.size(); ++i) { + printf("%f %d\n", supports_.radii2[i], i); + } +} \ No newline at end of file diff --git a/src/pcms/interpolator/interpolation_base.h b/src/pcms/interpolator/interpolation_base.h new file mode 100644 index 00000000..15c3e5d7 --- /dev/null +++ b/src/pcms/interpolator/interpolation_base.h @@ -0,0 +1,73 @@ +// +// Created by Fuad Hasan on 1/12/25. +// +#include "MLSInterpolation.hpp" +#include "adj_search.hpp" +#include +#include +#include +#include + +#ifndef PCMS_INTERPOLATION_BASE_H +#define PCMS_INTERPOLATION_BASE_H + +class InterpolationBase +{ + /** + * @brief Evaluate the interpolation + * @param source_field The field to interpolate from + * @param target_field The field to interpolate to + */ + virtual void eval( + const pcms::ScalarArrayView& source_field, + pcms::ScalarArrayView& target_field) = 0; +}; + +/** + * @brief Moving Least Square Radial Basis Function Interpolation + */ +class MLSInterpolationHandler : public InterpolationBase +{ + +public: + void eval( + const pcms::ScalarArrayView& source_field, + pcms::ScalarArrayView& target_field) + override; + + MLSInterpolationHandler(const std::string& source_mesh_fname, + const std::string& target_mesh_fname, double radius, + bool adapt_radius = true); + + MLSInterpolationHandler(Omega_h::Mesh source_mesh, Omega_h::Mesh target_mesh, + double radius, bool adapt_radius = true); + + MLSInterpolationHandler(const std::string& source_mesh_fname, double radius, + bool adapt_radius = true); + + MLSInterpolationHandler(Omega_h::Mesh source_mesh, double radius, + bool adapt_radius = true); + +private: + double radius_; + bool adapt_radius_; + bool single_mesh_ = false; + + std::string interpolation_type_; + + Omega_h::Library library_; + Omega_h::Mesh source_mesh_; + // TODO: handle what to do with this when only 1 mesh is provided + Omega_h::Mesh target_mesh_; + Omega_h::Reals source_coords_; + Omega_h::Reals target_coords_; + + SupportResults supports_; + + Omega_h::HostWrite target_field_; + Omega_h::HostWrite source_field_; + + void find_supports(int min_req_supports = 12); +}; + +#endif // PCMS_INTERPOLATION_BASE_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e4d8e4e7..b73a9f30 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -121,6 +121,7 @@ if(Catch2_FOUND) test_mls_basis.cpp test_rbf_interp.cpp test_linear_solver.cpp + test_interpolation_class.cpp ) endif () add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES}) diff --git a/test/test_interpolation_class.cpp b/test/test_interpolation_class.cpp new file mode 100644 index 00000000..b2191f4e --- /dev/null +++ b/test/test_interpolation_class.cpp @@ -0,0 +1,73 @@ +// +// Created by hasanm4 on 1/17/25. +// + +#include + +#include +#include +#include + +#include +#include + +void create_sinx_cosy_data( + Omega_h::Mesh& mesh, + pcms::ScalarArrayView& sinxcosy) +{ + // get the bounding box of the mesh + Omega_h::BBox<2> bb = Omega_h::get_bounding_box<2>(&mesh); + double dx = bb.max[0] - bb.min[0]; + double dy = bb.max[1] - bb.min[1]; + + auto coords = mesh.coords(); + auto nnodes = mesh.nverts(); + Omega_h::Write sinxcosytag(nnodes); + + auto assignSinCos = OMEGA_H_LAMBDA(int node) + { + auto x = (coords[2 * node + 0] - bb.min[0]) / dx * 2.0 * M_PI; + auto y = (coords[2 * node + 1] - bb.min[1]) / dy * 2.0 * M_PI; + sinxcosytag[node] = 2 + sin(x) * cos(y); + }; + Omega_h::parallel_for(nnodes, assignSinCos, "assignSinCos"); + Omega_h::HostWrite sinxcosyhost(sinxcosytag); + + OMEGA_H_CHECK_PRINTF( + nnodes = sinxcosy.size(), + "Given ScalarArrayView is not properly sized %d vs %zu\n", nnodes, + sinxcosy.size()); + + for (int i = 0; i < nnodes; ++i) { + sinxcosy[i] = sinxcosyhost[i]; + } +} + +TEST_CASE("Test MLSInterpolationHandler: Single Mesh") +{ + auto lib = Library{}; + auto world = lib.world(); + auto source_mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false); + + printf("[INFO] Mesh created with %d vertices and %d faces\n", source_mesh.nverts(), + source_mesh.nfaces()); + + auto mls_single = MLSInterpolationHandler(source_mesh, 0.09, true); + + Omega_h::HostWrite source_data_host_write(source_mesh.nverts()); + Omega_h::HostWrite interpolated_data_hwrite(source_mesh.nfaces()); + Omega_h::HostWrite target_data_expected(source_mesh.nfaces()); + + pcms::ScalarArrayView sourceArrayView( + source_data_host_write.data(), source_data_host_write.size()); + pcms::ScalarArrayView interpolatedArrayView( + interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); + pcms::ScalarArrayView + sinxcosyarrayview_target_expected(target_data_expected.data(), + target_data_expected.size()); + + create_sinx_cosy_data(source_mesh, sourceArrayView); + + mls_single.eval(sourceArrayView, interpolatedArrayView); +} \ No newline at end of file From a1a4196149724902594bbc1323f330d653c7ad8a Mon Sep 17 00:00:00 2001 From: Fuad Hasibul Hasan Date: Wed, 22 Jan 2025 18:17:49 -0500 Subject: [PATCH 2/4] working single mesh interpolation class the values are not being compared yet --- src/pcms/interpolator/interpolation_base.cpp | 142 +++++++------------ src/pcms/interpolator/interpolation_base.h | 48 ++++--- test/test_interpolation_class.cpp | 78 +++++++--- 3 files changed, 133 insertions(+), 135 deletions(-) diff --git a/src/pcms/interpolator/interpolation_base.cpp b/src/pcms/interpolator/interpolation_base.cpp index f1e20be1..15b6b17e 100644 --- a/src/pcms/interpolator/interpolation_base.cpp +++ b/src/pcms/interpolator/interpolation_base.cpp @@ -13,7 +13,7 @@ Omega_h::Reals getCentroids(Omega_h::Mesh& mesh) mesh.dim() == 2, "Only 2D meshes are supported but found %d\n", mesh.dim()); const auto& coords = mesh.coords(); - Omega_h::Write centroids(mesh.nfaces() * 2, 0.0); + Omega_h::Write centroids(mesh.nfaces() * mesh.dim(), 0.0); auto face2node = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; Omega_h::parallel_for( @@ -29,89 +29,40 @@ Omega_h::Reals getCentroids(Omega_h::Mesh& mesh) return {centroids}; } -MLSInterpolationHandler::MLSInterpolationHandler( - const std::string& source_mesh_fname, double radius, bool adapt_radius) - : radius_(radius), adapt_radius_(adapt_radius) -{ - single_mesh_ = true; - library_ = Omega_h::Library(nullptr, nullptr); - - source_mesh_ = Omega_h::binary::read(source_mesh_fname, &library_); - source_coords_ = source_mesh_.coords(); - target_coords_ = getCentroids(source_mesh_); - - OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2, - "Only 2D meshes are supported but found %d\n", - source_mesh_.dim()); - - source_field_ = - Omega_h::HostWrite(source_mesh_.nverts(), "source field"); - - // used nfaces to interpolate to the centroid of the faces - target_field_ = - Omega_h::HostWrite(source_mesh_.nfaces(), "target field"); - - find_supports(); -} - -MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh source_mesh, +MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh& source_mesh, double radius, bool adapt_radius) - : radius_(radius), adapt_radius_(adapt_radius) + : source_mesh_(source_mesh), + target_mesh_(source_mesh), + radius_(radius), + adapt_radius_(adapt_radius) { single_mesh_ = true; - source_mesh_ = std::move(source_mesh); - source_coords_ = source_mesh_.coords(); - target_coords_ = getCentroids(source_mesh_); + target_coords_ = source_mesh_.coords(); + source_coords_ = getCentroids(source_mesh_); OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2, "Only 2D meshes are supported but found %d\n", source_mesh_.dim()); source_field_ = - Omega_h::HostWrite(source_mesh_.nverts(), "source field"); + Omega_h::HostWrite(source_mesh_.nfaces(), "source field"); target_field_ = - Omega_h::HostWrite(source_mesh_.nfaces(), "target field"); -} - -// TODO: delegate this constructor to a single one -MLSInterpolationHandler::MLSInterpolationHandler( - const std::string& source_mesh_fname, const std::string& target_mesh_fname, - const double radius, const bool adapt_radius) - : radius_(radius), adapt_radius_(adapt_radius) -{ - // TODO: take argc and argv - library_ = Omega_h::Library(nullptr, nullptr); - - source_mesh_ = Omega_h::binary::read(source_mesh_fname, &library_); - target_mesh_ = Omega_h::binary::read(target_mesh_fname, &library_); - - OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2 && target_mesh_.dim() == 2, - "Only 2D meshes are supported but found %d, %d\n", - source_mesh_.dim(), target_mesh_.dim()); + Omega_h::HostWrite(source_mesh_.nverts(), "target field"); - source_coords_ = source_mesh_.coords(); - target_coords_ = target_mesh_.coords(); - - source_field_ = - Omega_h::HostWrite(source_mesh_.nverts(), "source field"); - target_field_ = - Omega_h::HostWrite(target_mesh_.nverts(), "target field"); - - find_supports(); + find_supports(10); } -MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh source_mesh, - Omega_h::Mesh target_mesh, +MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh& source_mesh, + Omega_h::Mesh& target_mesh, const double radius, const bool adapt_radius) - : radius_(radius), adapt_radius_(adapt_radius) + : source_mesh_(source_mesh), + target_mesh_(target_mesh), + radius_(radius), + adapt_radius_(adapt_radius) { - // TODO: check if move works for mesh - source_mesh_ = std::move(source_mesh); - target_mesh_ = std::move(target_mesh); - OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2 && target_mesh_.dim() == 2, "Only 2D meshes are supported but found %d, %d\n", source_mesh_.dim(), target_mesh_.dim()); @@ -124,12 +75,12 @@ MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh source_mesh, target_field_ = Omega_h::HostWrite(target_mesh_.nverts(), "target field"); - find_supports(); + find_supports(10); } // TODO : find way to avoid this copy void copyHostScalarArrayView2HostWrite( - const pcms::ScalarArrayView& source, + pcms::ScalarArrayView source, Omega_h::HostWrite& target) { OMEGA_H_CHECK_PRINTF( @@ -143,7 +94,7 @@ void copyHostScalarArrayView2HostWrite( } void copyHostWrite2ScalarArrayView( const Omega_h::HostWrite& source, - pcms::ScalarArrayView& target) + pcms::ScalarArrayView target) { OMEGA_H_CHECK_PRINTF( source.size() == target.size(), @@ -156,24 +107,18 @@ void copyHostWrite2ScalarArrayView( } void MLSInterpolationHandler::eval( - const pcms::ScalarArrayView& source_field, - pcms::ScalarArrayView& target_field) + pcms::ScalarArrayView source_field, + pcms::ScalarArrayView target_field) { - OMEGA_H_CHECK_PRINTF(source_field.size() == - source_coords_.size() / source_mesh_.dim(), - "Source Data and Source Points size mismatch: %zu %d\n", - source_field.size(), source_coords_.size()); - if (single_mesh_) { - OMEGA_H_CHECK_PRINTF( - target_field.size() == source_mesh_.nfaces(), - "Target Data and Target Points size mismatch: %zu %d\n", - target_field.size(), target_coords_.size()); - } else { - OMEGA_H_CHECK_PRINTF( - target_field.size() == target_coords_.size() / target_mesh_.dim(), - "Target Data and Target Points size mismatch: %zu %d\n", - target_field.size(), target_coords_.size()); - } + OMEGA_H_CHECK_PRINTF( + target_field.size() == target_coords_.size() / target_mesh_.dim(), + "Source Data and Source Points size mismatch: %zu %d\n", + target_field.size(), target_coords_.size() / target_mesh_.dim()); + + OMEGA_H_CHECK_PRINTF( + source_field.size() == source_coords_.size() / source_mesh_.dim(), + "Target Data and Target Points size mismatch: %zu %d\n", + source_field.size(), source_coords_.size() / source_mesh_.dim()); copyHostScalarArrayView2HostWrite(source_field, source_field_); @@ -185,23 +130,36 @@ void MLSInterpolationHandler::eval( copyHostWrite2ScalarArrayView(target_field_, target_field); } -void MLSInterpolationHandler::find_supports( - const int min_req_support) // ? Needed? +void MLSInterpolationHandler::find_supports(const int min_req_support) { if (single_mesh_) { supports_ = searchNeighbors(source_mesh_, radius_, min_req_support, adapt_radius_); - } else { // not single mesh + } else { // two mesh : vert to vert supports_ = searchNeighbors(source_mesh_, target_mesh_, radius_, min_req_support, adapt_radius_); } OMEGA_H_CHECK_PRINTF( - supports_.radii2[0] < 1e-10, + supports_.radii2[0] > 1e-10, "Radius squared has to be more than zero found found [0] = %f\n", supports_.radii2[0]); +} + +size_t MLSInterpolationHandler::getSourceSize() +{ + if (single_mesh_) { + return source_mesh_.nfaces(); + } else { + return source_mesh_.nverts(); + } +} - for (int i = 0; i < supports_.radii2.size(); ++i) { - printf("%f %d\n", supports_.radii2[i], i); +size_t MLSInterpolationHandler::getTargetSize() +{ + if (single_mesh_) { + return source_mesh_.nverts(); + } else { + return target_mesh_.nverts(); } } \ No newline at end of file diff --git a/src/pcms/interpolator/interpolation_base.h b/src/pcms/interpolator/interpolation_base.h index 15c3e5d7..d4890166 100644 --- a/src/pcms/interpolator/interpolation_base.h +++ b/src/pcms/interpolator/interpolation_base.h @@ -19,8 +19,8 @@ class InterpolationBase * @param target_field The field to interpolate to */ virtual void eval( - const pcms::ScalarArrayView& source_field, - pcms::ScalarArrayView& target_field) = 0; + pcms::ScalarArrayView source_field, + pcms::ScalarArrayView target_field) = 0; }; /** @@ -31,34 +31,42 @@ class MLSInterpolationHandler : public InterpolationBase public: void eval( - const pcms::ScalarArrayView& source_field, - pcms::ScalarArrayView& target_field) - override; + pcms::ScalarArrayView source_field, + pcms::ScalarArrayView target_field) override; - MLSInterpolationHandler(const std::string& source_mesh_fname, - const std::string& target_mesh_fname, double radius, - bool adapt_radius = true); - - MLSInterpolationHandler(Omega_h::Mesh source_mesh, Omega_h::Mesh target_mesh, - double radius, bool adapt_radius = true); + /** + * @brief Vertex to Vertex interpolation for two given meshes + * @param source_mesh The source mesh + * @param target_mesh The target mesh + * @param radius The cutoff radius for the MLS interpolation + * @param adapt_radius Whether to adapt the radius based on the local density + */ + MLSInterpolationHandler(Omega_h::Mesh& source_mesh, + Omega_h::Mesh& target_mesh, double radius, + bool adapt_radius = true); - MLSInterpolationHandler(const std::string& source_mesh_fname, double radius, - bool adapt_radius = true); + /** + * @brief Centroids to Vertices interpolation for a single mesh + * @param source_mesh The source mesh + * @param radius The cutoff radius for the MLS interpolation + * @param adapt_radius Whether to adapt the radius based on the local density + */ + MLSInterpolationHandler(Omega_h::Mesh& source_mesh, double radius, + bool adapt_radius = true); - MLSInterpolationHandler(Omega_h::Mesh source_mesh, double radius, - bool adapt_radius = true); + size_t getSourceSize(); + size_t getTargetSize(); private: double radius_; bool adapt_radius_; bool single_mesh_ = false; - std::string interpolation_type_; + //std::string interpolation_type_; - Omega_h::Library library_; - Omega_h::Mesh source_mesh_; + Omega_h::Mesh& source_mesh_; // TODO: handle what to do with this when only 1 mesh is provided - Omega_h::Mesh target_mesh_; + Omega_h::Mesh& target_mesh_; Omega_h::Reals source_coords_; Omega_h::Reals target_coords_; @@ -67,7 +75,7 @@ class MLSInterpolationHandler : public InterpolationBase Omega_h::HostWrite target_field_; Omega_h::HostWrite source_field_; - void find_supports(int min_req_supports = 12); + void find_supports(int min_req_supports = 10); }; #endif // PCMS_INTERPOLATION_BASE_H diff --git a/test/test_interpolation_class.cpp b/test/test_interpolation_class.cpp index b2191f4e..776f2a91 100644 --- a/test/test_interpolation_class.cpp +++ b/test/test_interpolation_class.cpp @@ -11,36 +11,55 @@ #include #include -void create_sinx_cosy_data( - Omega_h::Mesh& mesh, - pcms::ScalarArrayView& sinxcosy) +/** + * @brief Create sin(x)cos(y) + 2 at each vertex of the mesh + */ +void createSinxCosyAtVertex(Omega_h::Mesh& mesh, + Omega_h::Write& sinxcosy) { + OMEGA_H_CHECK_PRINTF(sinxcosy.size() == mesh.nverts(), + "Given ScalarArrayView is not properly sized %d vs %d\n", + mesh.nverts(), sinxcosy.size()); + // get the bounding box of the mesh Omega_h::BBox<2> bb = Omega_h::get_bounding_box<2>(&mesh); double dx = bb.max[0] - bb.min[0]; double dy = bb.max[1] - bb.min[1]; - auto coords = mesh.coords(); - auto nnodes = mesh.nverts(); - Omega_h::Write sinxcosytag(nnodes); + const auto& coords = mesh.coords(); auto assignSinCos = OMEGA_H_LAMBDA(int node) { auto x = (coords[2 * node + 0] - bb.min[0]) / dx * 2.0 * M_PI; auto y = (coords[2 * node + 1] - bb.min[1]) / dy * 2.0 * M_PI; - sinxcosytag[node] = 2 + sin(x) * cos(y); + sinxcosy[node] = 2 + Kokkos::sin(x) * Kokkos::cos(y); }; - Omega_h::parallel_for(nnodes, assignSinCos, "assignSinCos"); - Omega_h::HostWrite sinxcosyhost(sinxcosytag); + Omega_h::parallel_for(mesh.nverts(), assignSinCos, "assignSinCos"); +} - OMEGA_H_CHECK_PRINTF( - nnodes = sinxcosy.size(), - "Given ScalarArrayView is not properly sized %d vs %zu\n", nnodes, - sinxcosy.size()); +void node2CentroidInterpolation(Omega_h::Mesh& mesh, + Omega_h::Write& node_vals, + Omega_h::Write& centroid_vals) +{ + OMEGA_H_CHECK_PRINTF(node_vals.size() == mesh.nverts(), + "Given ScalarArrayView is not properly sized %d vs %d\n", + mesh.nverts(), node_vals.size()); + OMEGA_H_CHECK_PRINTF(centroid_vals.size() == mesh.nfaces(), + "Given ScalarArrayView is not properly sized %d vs %d\n", + mesh.nfaces(), centroid_vals.size()); - for (int i = 0; i < nnodes; ++i) { - sinxcosy[i] = sinxcosyhost[i]; - } + const auto& face2node = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + auto coords = mesh.coords(); + Omega_h::LO nfaces = mesh.nfaces(); + + auto averageSinCos = OMEGA_H_LAMBDA(Omega_h::LO face) + { + auto faceNodes = Omega_h::gather_verts<3>(face2node, face); + Omega_h::Real sum = node_vals[faceNodes[0]] + node_vals[faceNodes[1]] + + node_vals[faceNodes[2]]; + centroid_vals[face] = sum / 3.0; + }; + Omega_h::parallel_for(nfaces, averageSinCos, "averageSinCos"); } TEST_CASE("Test MLSInterpolationHandler: Single Mesh") @@ -50,14 +69,21 @@ TEST_CASE("Test MLSInterpolationHandler: Single Mesh") auto source_mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false); - printf("[INFO] Mesh created with %d vertices and %d faces\n", source_mesh.nverts(), - source_mesh.nfaces()); + printf("[INFO] Mesh created with %d vertices and %d faces\n", + source_mesh.nverts(), source_mesh.nfaces()); - auto mls_single = MLSInterpolationHandler(source_mesh, 0.09, true); + auto mls_single = MLSInterpolationHandler(source_mesh, 0.12, true); - Omega_h::HostWrite source_data_host_write(source_mesh.nverts()); - Omega_h::HostWrite interpolated_data_hwrite(source_mesh.nfaces()); - Omega_h::HostWrite target_data_expected(source_mesh.nfaces()); + Omega_h::Write sinxcosy_node(source_mesh.nverts(), + "sinxcosy_node"); + createSinxCosyAtVertex(source_mesh, sinxcosy_node); + Omega_h::Write sinxcosy_centroid(source_mesh.nfaces(), + "sinxcosy_centroid"); + node2CentroidInterpolation(source_mesh, sinxcosy_node, sinxcosy_centroid); + + Omega_h::HostWrite source_data_host_write(sinxcosy_centroid); + Omega_h::HostWrite interpolated_data_hwrite(source_mesh.nverts()); + Omega_h::HostWrite target_data_expected(source_mesh.nverts()); pcms::ScalarArrayView sourceArrayView( source_data_host_write.data(), source_data_host_write.size()); @@ -67,7 +93,13 @@ TEST_CASE("Test MLSInterpolationHandler: Single Mesh") sinxcosyarrayview_target_expected(target_data_expected.data(), target_data_expected.size()); - create_sinx_cosy_data(source_mesh, sourceArrayView); + OMEGA_H_CHECK_PRINTF(sourceArrayView.size() == mls_single.getSourceSize(), + "Source size mismatch: %zu vs %zu\n", + sourceArrayView.size(), mls_single.getSourceSize()); + OMEGA_H_CHECK_PRINTF( + interpolatedArrayView.size() == mls_single.getTargetSize(), + "Target size mismatch: %zu vs %zu\n", interpolatedArrayView.size(), + mls_single.getTargetSize()); mls_single.eval(sourceArrayView, interpolatedArrayView); } \ No newline at end of file From 3ea79c2fafea7ef7626190475bc5d1a60452a379 Mon Sep 17 00:00:00 2001 From: Fuad Hasibul Hasan Date: Wed, 22 Jan 2025 18:51:32 -0500 Subject: [PATCH 3/4] interpolator test compares each entry with expected passing with 10% error margin. [future] find best degree and min required support to minimize the epsilon of the test --- src/pcms/interpolator/interpolation_base.cpp | 22 +++++++----- src/pcms/interpolator/interpolation_base.h | 11 ++++-- test/test_interpolation_class.cpp | 35 +++++++++++++++++--- 3 files changed, 52 insertions(+), 16 deletions(-) diff --git a/src/pcms/interpolator/interpolation_base.cpp b/src/pcms/interpolator/interpolation_base.cpp index 15b6b17e..160323df 100644 --- a/src/pcms/interpolator/interpolation_base.cpp +++ b/src/pcms/interpolator/interpolation_base.cpp @@ -31,10 +31,13 @@ Omega_h::Reals getCentroids(Omega_h::Mesh& mesh) MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh& source_mesh, double radius, - bool adapt_radius) + uint min_req_support, + uint degree, bool adapt_radius) : source_mesh_(source_mesh), target_mesh_(source_mesh), radius_(radius), + min_req_supports_(min_req_support), + degree_(degree), adapt_radius_(adapt_radius) { single_mesh_ = true; @@ -51,16 +54,17 @@ MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh& source_mesh, target_field_ = Omega_h::HostWrite(source_mesh_.nverts(), "target field"); - find_supports(10); + find_supports(min_req_supports_); } -MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh& source_mesh, - Omega_h::Mesh& target_mesh, - const double radius, - const bool adapt_radius) +MLSInterpolationHandler::MLSInterpolationHandler( + Omega_h::Mesh& source_mesh, Omega_h::Mesh& target_mesh, const double radius, + uint min_req_support, uint degree, const bool adapt_radius) : source_mesh_(source_mesh), target_mesh_(target_mesh), radius_(radius), + min_req_supports_(min_req_support), + degree_(degree), adapt_radius_(adapt_radius) { OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2 && target_mesh_.dim() == 2, @@ -75,7 +79,7 @@ MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh& source_mesh, target_field_ = Omega_h::HostWrite(target_mesh_.nverts(), "target field"); - find_supports(10); + find_supports(min_req_supports_); } // TODO : find way to avoid this copy @@ -124,13 +128,13 @@ void MLSInterpolationHandler::eval( auto target_field_write = mls_interpolation( Omega_h::Reals(source_field_), source_coords_, target_coords_, supports_, 2, - 4, supports_.radii2, RadialBasisFunction::RBF_GAUSSIAN); + degree_, supports_.radii2, RadialBasisFunction::RBF_GAUSSIAN); target_field_ = Omega_h::HostWrite(target_field_write); copyHostWrite2ScalarArrayView(target_field_, target_field); } -void MLSInterpolationHandler::find_supports(const int min_req_support) +void MLSInterpolationHandler::find_supports(const uint min_req_support) { if (single_mesh_) { supports_ = diff --git a/src/pcms/interpolator/interpolation_base.h b/src/pcms/interpolator/interpolation_base.h index d4890166..deddecf2 100644 --- a/src/pcms/interpolator/interpolation_base.h +++ b/src/pcms/interpolator/interpolation_base.h @@ -39,10 +39,14 @@ class MLSInterpolationHandler : public InterpolationBase * @param source_mesh The source mesh * @param target_mesh The target mesh * @param radius The cutoff radius for the MLS interpolation + * @param min_req_supports The minimum number of source locations required for + * interpolation + * @param degree The degree of the polynomial used in the MLS interpolation * @param adapt_radius Whether to adapt the radius based on the local density */ MLSInterpolationHandler(Omega_h::Mesh& source_mesh, Omega_h::Mesh& target_mesh, double radius, + uint min_req_supports = 10, uint degree = 3, bool adapt_radius = true); /** @@ -52,6 +56,7 @@ class MLSInterpolationHandler : public InterpolationBase * @param adapt_radius Whether to adapt the radius based on the local density */ MLSInterpolationHandler(Omega_h::Mesh& source_mesh, double radius, + uint min_req_supports = 10, uint degree = 3, bool adapt_radius = true); size_t getSourceSize(); @@ -61,8 +66,10 @@ class MLSInterpolationHandler : public InterpolationBase double radius_; bool adapt_radius_; bool single_mesh_ = false; + uint degree_; + uint min_req_supports_; - //std::string interpolation_type_; + // InterpolationType interpolation_type_; Omega_h::Mesh& source_mesh_; // TODO: handle what to do with this when only 1 mesh is provided @@ -75,7 +82,7 @@ class MLSInterpolationHandler : public InterpolationBase Omega_h::HostWrite target_field_; Omega_h::HostWrite source_field_; - void find_supports(int min_req_supports = 10); + void find_supports(uint min_req_supports = 10); }; #endif // PCMS_INTERPOLATION_BASE_H diff --git a/test/test_interpolation_class.cpp b/test/test_interpolation_class.cpp index 776f2a91..c844e7a0 100644 --- a/test/test_interpolation_class.cpp +++ b/test/test_interpolation_class.cpp @@ -11,6 +11,10 @@ #include #include +bool isClose(Omega_h::HostWrite& array1, + Omega_h::HostWrite& array2, + double percent_diff = 0.1); + /** * @brief Create sin(x)cos(y) + 2 at each vertex of the mesh */ @@ -72,7 +76,7 @@ TEST_CASE("Test MLSInterpolationHandler: Single Mesh") printf("[INFO] Mesh created with %d vertices and %d faces\n", source_mesh.nverts(), source_mesh.nfaces()); - auto mls_single = MLSInterpolationHandler(source_mesh, 0.12, true); + auto mls_single = MLSInterpolationHandler(source_mesh, 0.12, 12, 3, true); Omega_h::Write sinxcosy_node(source_mesh.nverts(), "sinxcosy_node"); @@ -83,15 +87,12 @@ TEST_CASE("Test MLSInterpolationHandler: Single Mesh") Omega_h::HostWrite source_data_host_write(sinxcosy_centroid); Omega_h::HostWrite interpolated_data_hwrite(source_mesh.nverts()); - Omega_h::HostWrite target_data_expected(source_mesh.nverts()); + Omega_h::HostWrite exact_values_at_nodes(sinxcosy_node); pcms::ScalarArrayView sourceArrayView( source_data_host_write.data(), source_data_host_write.size()); pcms::ScalarArrayView interpolatedArrayView( interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); - pcms::ScalarArrayView - sinxcosyarrayview_target_expected(target_data_expected.data(), - target_data_expected.size()); OMEGA_H_CHECK_PRINTF(sourceArrayView.size() == mls_single.getSourceSize(), "Source size mismatch: %zu vs %zu\n", @@ -102,4 +103,28 @@ TEST_CASE("Test MLSInterpolationHandler: Single Mesh") mls_single.getTargetSize()); mls_single.eval(sourceArrayView, interpolatedArrayView); + + REQUIRE(isClose(exact_values_at_nodes, interpolated_data_hwrite, 10.0) == + true); +} + +bool isClose(Omega_h::HostWrite& array1, + Omega_h::HostWrite& array2, double percent_diff) +{ + if (array1.size() != array2.size()) { + fprintf(stderr, "[ERROR] Arrays are not of the same size: %d vs %d\n", + array1.size(), array2.size()); + return false; + } + + double eps = percent_diff / 100.0; + for (int i = 0; i < array1.size(); i++) { + if (std::abs(array1[i] - array2[i]) > eps * std::abs(array2[i])) { + fprintf(stderr, "[ERROR] Arrays differ at index %d: %.16f vs %.16f\n", i, + array1[i], array2[i]); + return false; + } + } + + return true; } \ No newline at end of file From c8a268e1f35023113a63c87240f2720570586b78 Mon Sep 17 00:00:00 2001 From: Fuad Hasibul Hasan Date: Wed, 22 Jan 2025 21:18:40 -0500 Subject: [PATCH 4/4] interpolator test passes for mesh to mesh with 10% error margin and [temporary] it doesn't work when a source and a target becomes coincident: discussed here: https://github.com/SCOREC/pcms/pull/148#discussion_r1926204199 --- src/pcms/interpolator/interpolation_base.cpp | 14 ++- test/test_interpolation_class.cpp | 112 ++++++++++++++----- 2 files changed, 96 insertions(+), 30 deletions(-) diff --git a/src/pcms/interpolator/interpolation_base.cpp b/src/pcms/interpolator/interpolation_base.cpp index 160323df..8da33ed8 100644 --- a/src/pcms/interpolator/interpolation_base.cpp +++ b/src/pcms/interpolator/interpolation_base.cpp @@ -126,6 +126,7 @@ void MLSInterpolationHandler::eval( copyHostScalarArrayView2HostWrite(source_field, source_field_); + // TODO: make the basis function a template or pass it as a parameter auto target_field_write = mls_interpolation( Omega_h::Reals(source_field_), source_coords_, target_coords_, supports_, 2, degree_, supports_.radii2, RadialBasisFunction::RBF_GAUSSIAN); @@ -144,10 +145,15 @@ void MLSInterpolationHandler::find_supports(const uint min_req_support) min_req_support, adapt_radius_); } - OMEGA_H_CHECK_PRINTF( - supports_.radii2[0] > 1e-10, - "Radius squared has to be more than zero found found [0] = %f\n", - supports_.radii2[0]); +#ifndef NDEBUG + Omega_h::HostRead hostRadii2(supports_.radii2); + for (size_t i = 0; i < hostRadii2.size(); ++i) { + OMEGA_H_CHECK_PRINTF( + hostRadii2[i] > 1e-10, + "Radius squared has to be more than zero found found [%zu] = %f\n", i, + hostRadii2[i]); + } +#endif } size_t MLSInterpolationHandler::getSourceSize() diff --git a/test/test_interpolation_class.cpp b/test/test_interpolation_class.cpp index c844e7a0..3d0a3e7c 100644 --- a/test/test_interpolation_class.cpp +++ b/test/test_interpolation_class.cpp @@ -11,6 +11,22 @@ #include #include +void translate_mesh(Omega_h::Mesh* mesh, Omega_h::Vector<2> translation_vector) { + auto coords = mesh->coords(); + auto nverts = mesh->nverts(); + auto out = Write(coords.size()); + + auto f = OMEGA_H_LAMBDA(LO i) { + auto coord = get_vector<2>(coords, i); + coord = coord + translation_vector; + set_vector<2>(out, i, coord); + }; + + parallel_for(nverts, f); + mesh->set_coords(Reals(out)); +} + + bool isClose(Omega_h::HostWrite& array1, Omega_h::HostWrite& array2, double percent_diff = 0.1); @@ -66,46 +82,90 @@ void node2CentroidInterpolation(Omega_h::Mesh& mesh, Omega_h::parallel_for(nfaces, averageSinCos, "averageSinCos"); } -TEST_CASE("Test MLSInterpolationHandler: Single Mesh") +TEST_CASE("Test MLSInterpolationHandler") { + fprintf(stdout, "[INFO] Starting MLS Interpolation Test...\n"); auto lib = Library{}; auto world = lib.world(); auto source_mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false); - printf("[INFO] Mesh created with %d vertices and %d faces\n", source_mesh.nverts(), source_mesh.nfaces()); - auto mls_single = MLSInterpolationHandler(source_mesh, 0.12, 12, 3, true); + Omega_h::Write source_sinxcosy_node(source_mesh.nverts(), + "source_sinxcosy_node"); + createSinxCosyAtVertex(source_mesh, source_sinxcosy_node); + + SECTION("Single Mesh") + { + fprintf(stdout, "\n-------------------- Single Mesh Interpolation Test Started --------------------\n"); + auto mls_single = MLSInterpolationHandler(source_mesh, 0.12, 12, 3, true); + + Omega_h::Write sinxcosy_centroid(source_mesh.nfaces(), + "sinxcosy_centroid"); + node2CentroidInterpolation(source_mesh, source_sinxcosy_node, + sinxcosy_centroid); + + Omega_h::HostWrite source_data_host_write(sinxcosy_centroid); + Omega_h::HostWrite interpolated_data_hwrite(source_mesh.nverts()); + Omega_h::HostWrite exact_values_at_nodes(source_sinxcosy_node); + + pcms::ScalarArrayView sourceArrayView( + source_data_host_write.data(), source_data_host_write.size()); + pcms::ScalarArrayView interpolatedArrayView( + interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); + + OMEGA_H_CHECK_PRINTF(sourceArrayView.size() == mls_single.getSourceSize(), + "Source size mismatch: %zu vs %zu\n", + sourceArrayView.size(), mls_single.getSourceSize()); + OMEGA_H_CHECK_PRINTF( + interpolatedArrayView.size() == mls_single.getTargetSize(), + "Target size mismatch: %zu vs %zu\n", interpolatedArrayView.size(), + mls_single.getTargetSize()); + + mls_single.eval(sourceArrayView, interpolatedArrayView); + + REQUIRE(isClose(exact_values_at_nodes, interpolated_data_hwrite, 10.0) == + true); + fprintf(stdout, "[****] Single Mesh Interpolation Test Passed with %.2f%% tolerance!\n", 10.0); + } + + SECTION("Double Mesh") + { + fprintf(stdout, "\n-------------------- Double Mesh Interpolation Test Started --------------------\n"); + auto target_mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 0.999, 0.999, 1, 17, 17, 0, false); + printf("[INFO] Target Mesh created with %d vertices and %d faces\n", + target_mesh.nverts(), target_mesh.nfaces()); - Omega_h::Write sinxcosy_node(source_mesh.nverts(), - "sinxcosy_node"); - createSinxCosyAtVertex(source_mesh, sinxcosy_node); - Omega_h::Write sinxcosy_centroid(source_mesh.nfaces(), - "sinxcosy_centroid"); - node2CentroidInterpolation(source_mesh, sinxcosy_node, sinxcosy_centroid); + // TODO: This is a way around. https://github.com/SCOREC/pcms/pull/148#discussion_r1926204199 + translate_mesh(&target_mesh, Omega_h::Vector<2>{(1.0-0.999)/2.0, (1.0-0.999)/2.0}); - Omega_h::HostWrite source_data_host_write(sinxcosy_centroid); - Omega_h::HostWrite interpolated_data_hwrite(source_mesh.nverts()); - Omega_h::HostWrite exact_values_at_nodes(sinxcosy_node); + auto mls_double = + MLSInterpolationHandler(source_mesh, target_mesh, 0.12, 12, 3, true); - pcms::ScalarArrayView sourceArrayView( - source_data_host_write.data(), source_data_host_write.size()); - pcms::ScalarArrayView interpolatedArrayView( - interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); + Omega_h::HostWrite source_data_host_write(source_sinxcosy_node); + Omega_h::HostWrite interpolated_data_hwrite( + mls_double.getTargetSize()); - OMEGA_H_CHECK_PRINTF(sourceArrayView.size() == mls_single.getSourceSize(), - "Source size mismatch: %zu vs %zu\n", - sourceArrayView.size(), mls_single.getSourceSize()); - OMEGA_H_CHECK_PRINTF( - interpolatedArrayView.size() == mls_single.getTargetSize(), - "Target size mismatch: %zu vs %zu\n", interpolatedArrayView.size(), - mls_single.getTargetSize()); + pcms::ScalarArrayView sourceArrayView( + source_data_host_write.data(), source_data_host_write.size()); + pcms::ScalarArrayView interpolatedArrayView( + interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); - mls_single.eval(sourceArrayView, interpolatedArrayView); + mls_double.eval(sourceArrayView, interpolatedArrayView); - REQUIRE(isClose(exact_values_at_nodes, interpolated_data_hwrite, 10.0) == - true); + Omega_h::Write exact_target_sinxcosy_node( + target_mesh.nverts(), "target_sinxcosy_node"); + createSinxCosyAtVertex(target_mesh, exact_target_sinxcosy_node); + Omega_h::HostWrite exact_target_sinxcosy_node_hwrite( + exact_target_sinxcosy_node); + + REQUIRE(isClose(interpolated_data_hwrite, exact_target_sinxcosy_node_hwrite, + 10.0) == true); + + fprintf(stdout, "[INFO] Double Mesh Interpolation Test Passed with %.2f%% tolerance!\n", 10.0); + } } bool isClose(Omega_h::HostWrite& array1,