From ec0ca918267f97ab4c0ce8b5a326af4b2be5d5fd Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Thu, 6 Aug 2026 09:36:25 -0400 Subject: [PATCH 1/4] test: share common setup helpers via field_test_utils.h Factor the repeated mesh construction, function-space creation, field sampling and evaluation boilerplate out of the individual test translation units and into field_test_utils.h: - BuildUnitSquare / MakeP1Space / MakeP0Space for the unit-square simplex meshes used across the transfer and evaluation tests - IntegrateP0Field / IntegrateP1Field for conservation checks - EvaluateAndAssemble for the sample-then-assemble pattern shared by the linear form integrator tests - StandardOutsideCoords2D for out-of-bounds policy tests SetField now dispatches on the arity of the supplied callable so the same helper covers 2D and 3D layouts, and it samples the layout's DOF coordinates directly instead of copying them into a fixed-rank view. The evaluation helpers use rank-2 output views through MakeRank2View rather than hand-built Rank2Views over rank-1 storage. No behavioral change to the tests themselves. --- test/field_test_utils.h | 202 +++++++++++++---- test/test_eqdsk.cpp | 37 +--- test/test_field_evaluation.cpp | 36 +-- test/test_field_exchange_planner.cpp | 4 +- test/test_field_interpolation.cpp | 30 +-- test/test_interpolation_class.cpp | 118 ++-------- test/test_interpolation_on_ltx_mesh.cpp | 48 ++-- test/test_localization_factory.cpp | 13 +- .../test_mesh_intersection_field_transfer.cpp | 209 +++++------------- ...st_omega_h_intersection_rhs_integrator.cpp | 138 +++--------- test/test_omega_h_lagrange_field.cpp | 38 +--- test/test_omega_h_mass_integrator.cpp | 42 +--- test/test_omega_h_mc_rhs_integrator.cpp | 108 ++------- test/test_point_evaluator.cpp | 96 +++----- ...lynomial_reconstruction_mls_evaluation.cpp | 112 +++------- test/test_spr_meshfields.cpp | 18 +- test/test_svd_serial.cpp | 4 +- test/test_uniform_grid_field.cpp | 82 +++---- test/test_xgc_reverse_classification.cpp | 2 +- 19 files changed, 440 insertions(+), 897 deletions(-) diff --git a/test/field_test_utils.h b/test/field_test_utils.h index 5aaac78e..2ddd2f45 100644 --- a/test/field_test_utils.h +++ b/test/field_test_utils.h @@ -3,6 +3,7 @@ #include #include +#include "pcms/configuration.h" #include "pcms/field/field.h" #include "pcms/field/field_data.h" #include "pcms/field/field_evaluator_factory.h" @@ -11,10 +12,22 @@ #include "pcms/field/out_of_bounds_policy.h" #include "pcms/coupler/field_serializer.h" #include "pcms/field/coordinate_system.h" -#include "pcms/localization/adj_search.hpp" #include "pcms/utility/arrays.h" #include "pcms/utility/memory_spaces.h" +#ifdef PCMS_ENABLE_OMEGA_H +#include +#include +#include +#include "pcms/field/function_space/lagrange.h" +#include "pcms/localization/adj_search.hpp" +#endif +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) +#include "pcms/transfer/linear_form_integrator.hpp" +#endif #include +#include +#include +#include #include #include @@ -36,6 +49,62 @@ inline std::vector StandardEvalCoords2D() return {0.1, 0.2, 0.5, 0.5, 0.7, 0.3, 0.9, 0.1, 0.2, 0.8}; } +// Points strictly outside a unit [0,1]^2 box mesh. +inline std::vector StandardOutsideCoords2D() +{ + return {-0.5, 0.5, 1.5, 0.5, 0.5, -0.5, 0.5, 1.5}; +} + +#ifdef PCMS_ENABLE_OMEGA_H +// Builds a unit-square 2D simplex mesh from the given element connectivity and +// adds the geometric classification tags required to build an Omega_h-backed +// Lagrange function space. +inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, + const Omega_h::LOs& ev2v) +{ + const Omega_h::Reals coords({ + 0.0, 0.0, // v0 + 1.0, 0.0, // v1 + 1.0, 1.0, // v2 + 0.0, 1.0 // v3 + }); + Omega_h::Mesh mesh(&lib); + Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); + for (Omega_h::Int dim = 0; dim <= 2; ++dim) { + mesh.add_tag( + dim, "class_dim", 1, + Omega_h::Read(mesh.nents(dim), Omega_h::I8(dim))); + mesh.add_tag( + dim, "class_id", 1, + Omega_h::Read(mesh.nents(dim), Omega_h::ClassId(0))); + } + return mesh; +} + +// Convenience overload: diagonal=0 splits the square along vertices (1,3), +// diagonal=1 along (0,2). +inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) +{ + return BuildUnitSquare(lib, (diagonal == 0) + ? Omega_h::LOs({0, 1, 3, 1, 2, 3}) + : Omega_h::LOs({0, 1, 2, 0, 2, 3})); +} + +inline std::shared_ptr MakeP1Space( + Omega_h::Mesh& mesh, const std::string& global_id_name = "global") +{ + return LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, CoordinateSystem::Cartesian, global_id_name, + LagrangeFunctionSpace::Backend::OmegaH); +} + +inline std::shared_ptr MakeP0Space(Omega_h::Mesh& mesh) +{ + return LagrangeFunctionSpace::FromMesh( + mesh, 0, 1, CoordinateSystem::Cartesian, "global", + LagrangeFunctionSpace::Backend::OmegaH); +} + inline bool AreArraysEqualUnordered( const Omega_h::HostRead& array1, const Omega_h::HostRead& array2, int start, int end) @@ -76,6 +145,40 @@ inline std::vector CopyOmegaHRealsToVector(const Omega_h::Reals& coords) coords_read.data() + coords_read.size()); } +inline double IntegrateP0Field(Omega_h::Mesh& mesh, const Field& field) +{ + const auto values = FlattenToRank1View(field.GetDOFHolderDataHost()); + const auto measures = Omega_h::measure_elements_real(&mesh); + const auto measures_h = Omega_h::HostRead(measures); + + double integral = 0.0; + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + integral += measures_h[e] * values[e]; + } + return integral; +} + +inline double IntegrateP1Field(Omega_h::Mesh& mesh, const Field& field) +{ + const auto values = FlattenToRank1View(field.GetDOFHolderDataHost()); + const auto measures = Omega_h::measure_elements_real(&mesh); + const auto measures_h = Omega_h::HostRead(measures); + const auto elem_verts_h = + Omega_h::HostRead(mesh.ask_elem_verts()); + const int verts_per_elem = mesh.dim() + 1; + + double integral = 0.0; + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + double avg = 0.0; + for (int k = 0; k < verts_per_elem; ++k) { + avg += values[elem_verts_h[verts_per_elem * e + k]]; + } + integral += measures_h[e] * (avg / verts_per_elem); + } + return integral; +} +#endif // PCMS_ENABLE_OMEGA_H + // Copy coordinates from device memory to a host view. // This handles potential layout mismatches between host and device memory // spaces. @@ -86,14 +189,8 @@ inline Kokkos::View CopyCoordinatesToHost( auto coords_view = Kokkos::View("coords_view", nents, dim); auto coords_view_device = - Kokkos::create_mirror(DeviceMemorySpace(), coords_view); - Kokkos::parallel_for( - "copy_coords_to_host_view", Kokkos::RangePolicy<>(0, nents), - KOKKOS_LAMBDA(int i) { - for (int d = 0; d < dim; ++d) { - coords_view_device(i, d) = coords_device(i, d); - } - }); + Kokkos::create_mirror_view(DeviceMemorySpace(), coords_view); + ConvertMismatchLayoutView2D(coords_view_device, coords_device); Kokkos::deep_copy(coords_view, coords_view_device); return coords_view; } @@ -122,44 +219,49 @@ inline std::vector EvaluateReferenceFunction(const std::vector& pts, expected_host.data() + expected_host.extent(0)); } -template +template struct SetFieldFunctor { - Kokkos::View data; - Kokkos::View coords; + DataView data; + CoordsView coords; Func f; - SetFieldFunctor(Kokkos::View data_, - Kokkos::View coords_, Func f_) + SetFieldFunctor(DataView data_, CoordsView coords_, Func f_) : data(data_), coords(coords_), f(f_) { } KOKKOS_INLINE_FUNCTION - void operator()(int i) const { data(i) = f(coords(i, 0), coords(i, 1)); } + void operator()(int i) const + { + if constexpr (std::is_invocable_v) { + data(i) = f(coords(i, 0), coords(i, 1), coords(i, 2)); + } else { + data(i) = f(coords(i, 0), coords(i, 1)); + } + } }; -// Set scalar DOF data by sampling func at each DOF-holder coordinate. +// Set scalar DOF data by sampling func at each DOF-holder coordinate. The +// arity of func selects the spatial dimension: func(x, y) for 2D layouts, +// func(x, y, z) for 3D. template inline void SetField(const FieldLayout& layout, FieldData& field, Func func) { using MemorySpace = typename ExecutionSpace::memory_space; + static_assert(std::is_invocable_v || + std::is_invocable_v, + "SetField requires func(x, y) or func(x, y, z)"); + auto dof_coords = layout.GetDOFHolderCoordinates().GetValues(); int n = static_cast(dof_coords.extent(0)); - Kokkos::View coords_device("coords_device", n); - Kokkos::parallel_for( - "field_test_utils_copy_coords", Kokkos::RangePolicy(0, n), - KOKKOS_LAMBDA(int i) { - coords_device(i, 0) = dof_coords(i, 0); - coords_device(i, 1) = dof_coords(i, 1); - }); Kokkos::View data_device("data_device", n); Kokkos::parallel_for("field_test_utils_set_field", Kokkos::RangePolicy(0, n), - SetFieldFunctor{data_device, coords_device, func}); + SetFieldFunctor{data_device, dof_coords, func}); auto data_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), data_device); @@ -253,12 +355,8 @@ void CheckEvaluation(const PointEvaluator& evaluator, { int n = static_cast(pts.size()) / 2; - Kokkos::View out_device("out_device", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out(out_device.data(), n, 1); - // Rank2View out(eval.data(), n, 1); - evaluator.Evaluate(field, out); + Kokkos::View out_device("out_device", n, 1); + evaluator.Evaluate(field, MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); @@ -267,8 +365,8 @@ void CheckEvaluation(const PointEvaluator& evaluator, for (int i = 0; i < n; ++i) { INFO("Point " << i << " (" << pts[2 * static_cast(i)] << ", " << pts[2 * static_cast(i) + 1] << ")" - << " got=" << out_host(i) << " expected=" << expected[i]); - REQUIRE(out_host(i) == Catch::Approx(expected[i]).margin(abs_tol)); + << " got=" << out_host(i, 0) << " expected=" << expected[i]); + REQUIRE(out_host(i, 0) == Catch::Approx(expected[i]).margin(abs_tol)); } } @@ -295,16 +393,13 @@ inline void CheckFillMode(const PointEvaluator& evaluator, { int n = static_cast(outside_pts.size()) / 2; - Kokkos::View out_device("out_device", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out(out_device.data(), n, 1); - evaluator.Evaluate(field, out); + Kokkos::View out_device("out_device", n, 1); + evaluator.Evaluate(field, MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); for (int i = 0; i < n; ++i) { - REQUIRE(out_host(i) == fill_value); + REQUIRE(out_host(i, 0) == fill_value); } } @@ -339,11 +434,8 @@ void CheckEvaluationWithFill(const Factory& factory, const Field& field, auto evaluator = factory->template CreatePointEvaluator( EvaluationRequest::FromCoordinates(device_coords.coordinate_view, policy)); - Kokkos::View out_device("out_device", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out(out_device.data(), n, 1); - evaluator->Evaluate(field, out); + Kokkos::View out_device("out_device", n, 1); + evaluator->Evaluate(field, MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); @@ -352,16 +444,34 @@ void CheckEvaluationWithFill(const Factory& factory, const Field& field, for (int i = 0; i < n; ++i) { INFO("Point " << i << " (" << pts[2 * static_cast(i)] << ", " << pts[2 * static_cast(i) + 1] << ")" - << " got=" << out_host(i)); + << " got=" << out_host(i, 0)); if (is_inside[i]) { INFO(" expected=" << expected[i]); - REQUIRE(out_host(i) == Catch::Approx(expected[i]).margin(abs_tol)); + REQUIRE(out_host(i, 0) == Catch::Approx(expected[i]).margin(abs_tol)); } else { - REQUIRE(out_host(i) == fill_value); + REQUIRE(out_host(i, 0) == fill_value); } } } +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) +// Evaluates source_field at the integrator's sample points and assembles the +// load vector. +inline void EvaluateAndAssemble( + LinearFormIntegrator& integrator, + const std::shared_ptr& source_space, + const Field& source_field) +{ + const auto& pts = integrator.GetIntegrationPoints(); + const std::size_t npts = pts.GetValues().extent(0); + auto evaluator = source_space->CreatePointEvaluator( + EvaluationRequest::FromCoordinates(pts)); + Kokkos::View sampled("sampled", npts, 1); + evaluator->Evaluate(source_field, MakeRank2View(sampled)); + integrator.Assemble(MakeConstRank2View(sampled)); +} +#endif // PCMS_ENABLE_PETSC && PCMS_ENABLE_MESHFIELDS + } // namespace pcms::test #endif // PCMS_TEST_FIELD_TEST_UTILS_H diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp index d546fbbe..c51f46b1 100644 --- a/test/test_eqdsk.cpp +++ b/test/test_eqdsk.cpp @@ -7,6 +7,7 @@ #include #include #include +#include "field_test_utils.h" #include #include @@ -212,42 +213,24 @@ TEST_CASE("EQDSKData with SplineFunctionSpace") const int num_eval_points = 3; - auto eval_coords_host = Kokkos::View( - "eval_coords_host", num_eval_points, 2); - for (size_t i = 0; i < num_eval_points; ++i) { - eval_coords_host(i, 0) = eval_coords[2 * i]; // R - eval_coords_host(i, 1) = eval_coords[2 * i + 1]; // Z - } - - auto eval_coords_device = - Kokkos::View("eval_coords_device", - num_eval_points, 2); - pcms::DeepCopyMismatchLayouts(eval_coords_device, eval_coords_host); - - auto coords_view = pcms::MakeRank2View(eval_coords_device); - auto coord_view = pcms::CoordinateView{ - CoordinateSystem::Cartesian, coords_view}; - auto eval_request = pcms::EvaluationRequest::FromCoordinates(coord_view); + auto device_coords = pcms::test::CreateDeviceCoordinateView( + eval_coords, CoordinateSystem::Cartesian); + auto eval_request = + pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view); auto evaluator = spline_space->CreatePointEvaluator(eval_request); - auto eval_results_1d = Kokkos::View( - "eval_results", num_eval_points); - - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View - eval_results(eval_results_1d.data(), num_eval_points, 1); - - evaluator->Evaluate(psi_field, eval_results); + Kokkos::View eval_results( + "eval_results", num_eval_points, 1); + evaluator->Evaluate(psi_field, pcms::MakeRank2View(eval_results)); auto results_host = Kokkos::create_mirror_view_and_copy( - pcms::HostMemorySpace(), eval_results_1d); + pcms::HostMemorySpace(), eval_results); // Verify that all results are finite for (int i = 0; i < num_eval_points; ++i) { - REQUIRE(std::isfinite(results_host(i))); + REQUIRE(std::isfinite(results_host(i, 0))); } } } diff --git a/test/test_field_evaluation.cpp b/test/test_field_evaluation.cpp index ff238e11..85af38b4 100644 --- a/test/test_field_evaluation.cpp +++ b/test/test_field_evaluation.cpp @@ -34,10 +34,10 @@ TEST_CASE("evaluate linear 2d omega_h_field") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckEvaluation( factory, field, pcms::test::StandardEvalCoords2D(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } #ifdef PCMS_ENABLE_MESHFIELDS @@ -50,35 +50,11 @@ TEST_CASE("evaluate quadratic 2d meshfields_field") mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::MeshFields); - // Quadratic DOF holders span vertices and edge midpoints; set them inline. - const auto nverts = mesh.nents(0); - const auto nedges = mesh.nents(1); - auto mesh_coords = mesh.coords(); - auto edge_verts = mesh.ask_verts_of(1); - - Omega_h::Write test_f(nverts + nedges); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - test_f[i] = sin_f(mesh_coords[2 * static_cast(i)], - mesh_coords[2 * static_cast(i) + 1]); - }); - Omega_h::parallel_for( - nedges, OMEGA_H_LAMBDA(int i) { - auto ep = Omega_h::gather_verts<2>(edge_verts, i); - Real cx = (mesh_coords[2 * static_cast(ep[0])] + - mesh_coords[2 * static_cast(ep[1])]) / - 2; - Real cy = (mesh_coords[2 * static_cast(ep[0]) + 1] + - mesh_coords[2 * static_cast(ep[1]) + 1]) / - 2; - test_f[nverts + i] = sin_f(cx, cy); - }); - - Omega_h::HostWrite test_f_host(test_f); + // Quadratic DOF holders span vertices and edge midpoints; the layout's DOF + // coordinates cover both, so SetField samples sin_f at every holder. auto field = factory->CreateFunction(); - field.GetData().SetDOFHolderDataHost( - pcms::Rank2View(test_f_host.data(), - test_f_host.size(), 1)); + pcms::test::SetField( + field, OMEGA_H_LAMBDA(Real x, Real y) { return sin_f(x, y); }); pcms::test::CheckEvaluation( factory, field, kEvalCoords, diff --git a/test/test_field_exchange_planner.cpp b/test/test_field_exchange_planner.cpp index 61b5bc54..99e6181b 100644 --- a/test/test_field_exchange_planner.cpp +++ b/test/test_field_exchange_planner.cpp @@ -51,9 +51,7 @@ TEST_CASE("GID messages insert headers around compact field payloads", "gid_message", plan.msg_size + 2 * static_cast(pcms::ent_offsets_len)); pcms::GenericFieldExchangePlanner planner; - planner.FillGidMessage(layout, plan, - pcms::Rank1View( - message.data(), message.size())); + planner.FillGidMessage(layout, plan, pcms::MakeRank1View(message)); const pcms::GO expected[] = { 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 1, 3, diff --git a/test/test_field_interpolation.cpp b/test/test_field_interpolation.cpp index bc88ab97..f8e2432d 100644 --- a/test/test_field_interpolation.cpp +++ b/test/test_field_interpolation.cpp @@ -55,35 +55,11 @@ TEST_CASE("interpolate quadratic 2d meshfields_field") auto factory2 = pcms::LagrangeFunctionSpace::FromMesh( mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::MeshFields); - auto layout = factory2->GetLayout(); - const auto nverts = mesh.nents(0); - const auto nedges = mesh.nents(1); - auto mesh_coords = mesh.coords(); - auto edge_verts = mesh.ask_verts_of(1); - Omega_h::Write test_f(nverts + nedges); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - Real x = mesh_coords[2 * i + 0]; - Real y = mesh_coords[2 * i + 1]; - test_f[i] = interpolation_linear_f(x, y); - }); - Omega_h::parallel_for( - nedges, OMEGA_H_LAMBDA(int i) { - auto endpoints = Omega_h::gather_verts<2>(edge_verts, i); - Real x0 = mesh_coords[2 * endpoints[0] + 0]; - Real y0 = mesh_coords[2 * endpoints[0] + 1]; - Real x1 = mesh_coords[2 * endpoints[1] + 0]; - Real y1 = mesh_coords[2 * endpoints[1] + 1]; - Real cx = (x0 + x1) / 2; - Real cy = (y0 + y1) / 2; - test_f[nverts + i] = interpolation_linear_f(cx, cy); - }); - - Omega_h::HostWrite test_f_host(test_f); auto field = factory2->CreateFunction(); auto interpolated = factory2->CreateFunction(); - field.SetDOFHolderDataHost(pcms::Rank2View( - test_f_host.data(), test_f_host.size(), 1)); + pcms::test::SetField( + field, + OMEGA_H_LAMBDA(Real x, Real y) { return interpolation_linear_f(x, y); }); pcms::Interpolator interp(*factory2, *factory2); interp.Apply(field, interpolated); diff --git a/test/test_interpolation_class.cpp b/test/test_interpolation_class.cpp index 7a34af4d..878f6288 100644 --- a/test/test_interpolation_class.cpp +++ b/test/test_interpolation_class.cpp @@ -10,36 +10,10 @@ #include #include #include +#include "field_test_utils.h" #include #include -#include - -bool areArraysEqualUnordered(const Omega_h::HostRead& array1, - const Omega_h::HostRead& array2, - int start, int end) -{ - // Ensure the indices are valid - assert(start >= 0 && end <= array1.size() && start <= end); - assert(start >= 0 && end <= array2.size() && start <= end); - - // Use frequency maps to count occurrences of each value - std::unordered_map freq1, freq2; - - for (int i = start; i < end; ++i) { - freq1[array1[i]]++; - freq2[array2[i]]++; - } - - // Compare the frequency maps - if (freq1 != freq2) { - pcms::printError("[ERROR] Arrays differ in the range [%d, %d)\n", start, - end); - return false; - } - - return true; -} void translate_mesh(Omega_h::Mesh* mesh, Omega_h::Vector<2> translation_vector) { @@ -135,28 +109,13 @@ TEST_CASE("Test MLSMeshInterpolation") auto mls_single = pcms::MLSMeshInterpolation(source_mesh, 0.12, 15, 3, true, 0.0, 5.0); - auto source_points_reals = - pcms::get_entity_centroids(source_mesh, Omega_h::FACE); - auto source_points_host = - Omega_h::HostRead(source_points_reals); - auto source_points_host_write = - Omega_h::HostWrite(source_points_host.size()); - for (int i = 0; i < source_points_host.size(); i++) { - source_points_host_write[i] = source_points_host[i]; - } - auto source_points_view = pcms::Rank1View( - source_points_host_write.data(), source_points_host_write.size()); - - auto target_points_reals = source_mesh.coords(); - auto target_points_host = - Omega_h::HostRead(target_points_reals); - auto target_points_host_write = - Omega_h::HostWrite(target_points_host.size()); - for (int i = 0; i < target_points_host.size(); i++) { - target_points_host_write[i] = target_points_host[i]; - } - auto target_points_view = pcms::Rank1View( - target_points_host_write.data(), target_points_host_write.size()); + auto source_points_vec = pcms::test::CopyOmegaHRealsToVector( + pcms::get_entity_centroids(source_mesh, Omega_h::FACE)); + auto source_points_view = pcms::make_array_view(source_points_vec); + + auto target_points_vec = + pcms::test::CopyOmegaHRealsToVector(source_mesh.coords()); + auto target_points_view = pcms::make_array_view(target_points_vec); REQUIRE(source_mesh.dim() == 2); pcms::printInfo("Point cloud based search...\n"); auto point_mls = pcms::MLSPointCloudInterpolation( @@ -173,14 +132,11 @@ TEST_CASE("Test MLSMeshInterpolation") source_mesh.nverts()); Omega_h::HostWrite exact_values_at_nodes(source_sinxcosy_node); - pcms::Rank1View sourceArrayView( - source_data_host_write.data(), source_data_host_write.size()); - pcms::Rank1View interpolatedArrayView( - interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); - pcms::Rank1View - point_cloud_interpolatedArrayView( - point_cloud_interpolated_data_hwrite.data(), - point_cloud_interpolated_data_hwrite.size()); + auto sourceArrayView = pcms::make_array_view(source_data_host_write); + auto interpolatedArrayView = + pcms::make_array_view(interpolated_data_hwrite); + auto point_cloud_interpolatedArrayView = + pcms::make_array_view(point_cloud_interpolated_data_hwrite); OMEGA_H_CHECK_PRINTF(sourceArrayView.size() == mls_single.getSourceSize(), "Source size mismatch: %zu vs %zu\n", @@ -222,33 +178,8 @@ TEST_CASE("Test MLSMeshInterpolation") ///*****************************// auto mesh_based_supports = mls_single.getSupports(); auto point_cloud_based_supports = point_mls.getSupports(); - auto mesh_based_support_ptr_host = - Omega_h::HostRead(mesh_based_supports.supports_ptr); - auto mesh_based_support_idx_host = - Omega_h::HostRead(mesh_based_supports.supports_idx); - auto point_cloud_based_support_ptr_host = - Omega_h::HostRead(point_cloud_based_supports.supports_ptr); - auto point_cloud_based_support_idx_host = - Omega_h::HostRead(point_cloud_based_supports.supports_idx); - - REQUIRE(point_cloud_based_support_idx_host.size() == - mesh_based_support_idx_host.size()); - REQUIRE(point_cloud_based_support_ptr_host.size() == - mesh_based_support_ptr_host.size()); - for (int i = 0; i < mesh_based_support_ptr_host.size(); i++) { - REQUIRE(point_cloud_based_support_ptr_host[i] == - mesh_based_support_ptr_host[i]); - } - - for (int i = 0; i < mesh_based_support_ptr_host.size() - 1; i++) { - auto start = mesh_based_support_ptr_host[i]; - auto end = mesh_based_support_ptr_host[i + 1]; - - bool isEqual = - areArraysEqualUnordered(mesh_based_support_idx_host, - point_cloud_based_support_idx_host, start, end); - REQUIRE(isEqual); - } + pcms::test::CheckSupportResultsEquivalent(point_cloud_based_supports, + mesh_based_supports); // Check if the point cloud interpolation is same as the MLS interpolation pcms::printDebugInfo("Interpolated data size: %d\n", @@ -289,10 +220,9 @@ TEST_CASE("Test MLSMeshInterpolation") Omega_h::HostWrite interpolated_data_hwrite( mls_double.getTargetSize()); - pcms::Rank1View sourceArrayView( - source_data_host_write.data(), source_data_host_write.size()); - pcms::Rank1View interpolatedArrayView( - interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); + auto sourceArrayView = pcms::make_array_view(source_data_host_write); + auto interpolatedArrayView = + pcms::make_array_view(interpolated_data_hwrite); mls_double.eval(sourceArrayView, interpolatedArrayView); @@ -326,10 +256,8 @@ TEST_CASE("MLSPointCloudInterpolation honors provided dimension in eval") } auto target_points = source_points; - auto source_points_view = pcms::Rank1View( - source_points.data(), source_points.size()); - auto target_points_view = pcms::Rank1View( - target_points.data(), target_points.size()); + auto source_points_view = pcms::make_array_view(source_points); + auto target_points_view = pcms::make_array_view(target_points); // Degree-1 polynomial that depends on z to catch accidental 2D behavior. auto source_values = Omega_h::HostWrite(27); @@ -341,10 +269,8 @@ TEST_CASE("MLSPointCloudInterpolation honors provided dimension in eval") } auto output_values = Omega_h::HostWrite(27, "output_values"); - auto source_values_view = pcms::Rank1View( - source_values.data(), source_values.size()); - auto output_values_view = pcms::Rank1View( - output_values.data(), output_values.size()); + auto source_values_view = pcms::make_array_view(source_values); + auto output_values_view = pcms::make_array_view(output_values); auto mls = pcms::MLSPointCloudInterpolation( source_points_view, target_points_view, 3, 2.5, 10, 1, true, 0.0, 5.0); diff --git a/test/test_interpolation_on_ltx_mesh.cpp b/test/test_interpolation_on_ltx_mesh.cpp index a348acae..141ebfa8 100644 --- a/test/test_interpolation_on_ltx_mesh.cpp +++ b/test/test_interpolation_on_ltx_mesh.cpp @@ -83,16 +83,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") printf("[INFO] Degas2 Mesh loaded from %s with %d elements\n", degas2_mesh_filename.c_str(), degas2_num_elems); const auto degas2_mesh_centroids_view = - pcms::Rank1View( - degas2_mesh_centroids_host.data(), degas2_mesh_centroids_host.size()); + pcms::make_array_view(degas2_mesh_centroids_host); auto xgc_mesh_points = read_xgc_mesh_nodes(ltx_mesh_base_filename + ".node"); const int xgc_num_nodes = xgc_mesh_points.size() / 2; printf("[INFO] XGC Mesh loaded from %s with %d points\n", ltx_mesh_base_filename.c_str(), xgc_num_nodes); - const auto xgc_mesh_points_view = - pcms::Rank1View( - xgc_mesh_points.data(), xgc_mesh_points.size()); + const auto xgc_mesh_points_view = pcms::make_array_view(xgc_mesh_points); auto xgc_to_degas2_interpolator = pcms::MLSPointCloudInterpolation( xgc_mesh_points_view, degas2_mesh_centroids_view, 2, 0.000001, 10, 1, true, @@ -120,18 +117,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") // ------------------ First Interpolation ------------------ // const auto density_at_xgc_nodes_view = - pcms::Rank1View(density_at_xgc_nodes.data(), - density_at_xgc_nodes.size()); - const auto temp_at_xgc_nodes_view = - pcms::Rank1View(temp_at_xgc_nodes.data(), - temp_at_xgc_nodes.size()); + pcms::make_array_view(density_at_xgc_nodes); + const auto temp_at_xgc_nodes_view = pcms::make_array_view(temp_at_xgc_nodes); const auto density_at_degas2_centroids_view = - pcms::Rank1View( - density_at_degas2_centroids.data(), density_at_degas2_centroids.size()); + pcms::make_array_view(density_at_degas2_centroids); const auto temp_at_degas2_centroids_view = - pcms::Rank1View( - temp_at_degas2_centroids.data(), temp_at_degas2_centroids.size()); + pcms::make_array_view(temp_at_degas2_centroids); Omega_h::HostWrite interpolated_xgc_density(degas2_num_elems); Omega_h::HostWrite interpolated_xgc_temp(degas2_num_elems); @@ -139,17 +131,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") Omega_h::HostWrite interpolated_degas2_temp(xgc_num_nodes); const auto interpolated_xgc_density_view = - pcms::Rank1View( - interpolated_xgc_density.data(), interpolated_xgc_density.size()); + pcms::make_array_view(interpolated_xgc_density); const auto interpolated_xgc_temp_view = - pcms::Rank1View( - interpolated_xgc_temp.data(), interpolated_xgc_temp.size()); + pcms::make_array_view(interpolated_xgc_temp); const auto interpolated_degas2_density_view = - pcms::Rank1View( - interpolated_degas2_density.data(), interpolated_degas2_density.size()); + pcms::make_array_view(interpolated_degas2_density); const auto interpolated_degas2_temp_view = - pcms::Rank1View( - interpolated_degas2_temp.data(), interpolated_degas2_temp.size()); + pcms::make_array_view(interpolated_degas2_temp); xgc_to_degas2_interpolator.eval(density_at_xgc_nodes_view, interpolated_xgc_density_view); @@ -174,21 +162,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") degas2_num_elems); const auto interpolated_back_density_at_xgc_nodes_view = - pcms::Rank1View( - interpolated_back_density_at_xgc_nodes.data(), - interpolated_back_density_at_xgc_nodes.size()); + pcms::make_array_view(interpolated_back_density_at_xgc_nodes); const auto interpolated_back_temp_at_xgc_nodes_view = - pcms::Rank1View( - interpolated_back_temp_at_xgc_nodes.data(), - interpolated_back_temp_at_xgc_nodes.size()); + pcms::make_array_view(interpolated_back_temp_at_xgc_nodes); const auto interpolated_back_density_at_degas2_centroids_view = - pcms::Rank1View( - interpolated_back_density_at_degas2_centroids.data(), - interpolated_back_density_at_degas2_centroids.size()); + pcms::make_array_view(interpolated_back_density_at_degas2_centroids); const auto interpolated_back_temp_at_degas2_centroids_view = - pcms::Rank1View( - interpolated_back_temp_at_degas2_centroids.data(), - interpolated_back_temp_at_degas2_centroids.size()); + pcms::make_array_view(interpolated_back_temp_at_degas2_centroids); degas2_to_xgc_interpolator.eval(interpolated_xgc_density_view, interpolated_back_density_at_xgc_nodes_view); diff --git a/test/test_localization_factory.cpp b/test/test_localization_factory.cpp index 571ac175..d65e8fc0 100644 --- a/test/test_localization_factory.cpp +++ b/test/test_localization_factory.cpp @@ -14,6 +14,7 @@ #include "pcms/localization/point_cloud_localization.h" #include "pcms/discretization/discretization/omega_h.hpp" #include "pcms/utility/arrays.h" +#include "pcms/utility/omega_h_array_utils.h" #include "field_test_utils.h" namespace @@ -46,17 +47,7 @@ TEST_CASE( auto target_device = pcms::test::CreateDeviceCoordinateView( target_coords, pcms::CoordinateSystem::Cartesian, dim); - auto coords_read = Omega_h::HostRead(source_coords); - - // auto coords_dev = Kokkos::create_mirror_view_and_copy( - // Kokkos::DefaultExecutionSpace{}, coords_host); - auto coords_dev = Kokkos::View( - "point_cloud_coords", mesh.nverts(), dim); - auto coords_host = Kokkos::create_mirror(pcms::HostMemorySpace(), coords_dev); - for (int i = 0; i < mesh.nverts(); ++i) - for (int d = 0; d < dim; ++d) - coords_host(i, d) = coords_read[i * dim + d]; - Kokkos::deep_copy(coords_dev, coords_host); + auto coords_dev = pcms::ConvertCoordsTo2D(source_coords, mesh.nverts(), dim); auto layout = std::make_shared( dim, coords_dev, pcms::CoordinateSystem::Cartesian); diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp index fbb34ddb..a627ab81 100644 --- a/test/test_mesh_intersection_field_transfer.cpp +++ b/test/test_mesh_intersection_field_transfer.cpp @@ -7,53 +7,14 @@ #include #include +#include #include "field_test_utils.h" -#include #include namespace { -// Builds a unit-square 2D simplex mesh from the given element connectivity and -// adds the geometric classification tags required to build an Omega_h-backed -// Lagrange function space. -Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, const Omega_h::LOs& ev2v) -{ - const Omega_h::Reals coords({ - 0.0, 0.0, // v0 - 1.0, 0.0, // v1 - 1.0, 1.0, // v2 - 0.0, 1.0 // v3 - }); - Omega_h::Mesh mesh(&lib); - Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); - for (Omega_h::Int dim = 0; dim <= 2; ++dim) { - mesh.add_tag( - dim, "class_dim", 1, - Omega_h::Read(mesh.nents(dim), Omega_h::I8(dim))); - mesh.add_tag( - dim, "class_id", 1, - Omega_h::Read(mesh.nents(dim), Omega_h::ClassId(0))); - } - return mesh; -} - -std::shared_ptr MakeP1Space( - Omega_h::Mesh& mesh, const std::string& global_id_name = "global") -{ - return pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, global_id_name, - pcms::LagrangeFunctionSpace::Backend::OmegaH); -} - -std::shared_ptr MakeP0Space(Omega_h::Mesh& mesh) -{ - return pcms::LagrangeFunctionSpace::FromMesh( - mesh, 0, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); -} - void AddReorderedVertexGlobalIds(Omega_h::Mesh& mesh) { mesh.add_tag( @@ -68,64 +29,6 @@ void AddSparseVertexGlobalIds(Omega_h::Mesh& mesh) Omega_h::Read({102, 7, 41, 19}, "sparse_global")); } -// Integral over the mesh of an order-0 (piecewise-constant) field whose values -// are stored in element order. Exact by construction. -double IntegrateP0Field(Omega_h::Mesh& mesh, - const pcms::Field& field) -{ - const auto values = pcms::FlattenToRank1View(field.GetDOFHolderDataHost()); - const auto elem_areas = Omega_h::measure_elements_real(&mesh); - const auto elem_areas_h = Omega_h::HostRead(elem_areas); - - double integral = 0.0; - for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { - integral += elem_areas_h[e] * values[e]; - } - return integral; -} - -// Per-element centroid coordinates in element order. -std::vector> ElementCentroids(Omega_h::Mesh& mesh) -{ - const auto coords_h = Omega_h::HostRead(mesh.coords()); - const auto elem_verts_h = - Omega_h::HostRead(mesh.ask_elem_verts()); - std::vector> centroids(mesh.nelems()); - for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { - double cx = 0.0, cy = 0.0; - for (int k = 0; k < 3; ++k) { - const Omega_h::LO v = elem_verts_h[3 * e + k]; - cx += coords_h[2 * v + 0]; - cy += coords_h[2 * v + 1]; - } - centroids[e] = {cx / 3.0, cy / 3.0}; - } - return centroids; -} - -// Integral over the mesh of an order-1 nodal field whose values are stored in -// vertex order (as produced by Field::GetDOFHolderDataHost for the -// Omega_h backend), computed exactly via per-element vertex averaging. -double IntegrateP1Field(Omega_h::Mesh& mesh, - const pcms::Field& field) -{ - const auto values = pcms::FlattenToRank1View(field.GetDOFHolderDataHost()); - const auto elem_areas = Omega_h::measure_elements_real(&mesh); - const auto elem_verts = mesh.ask_elem_verts(); - const auto elem_areas_h = Omega_h::HostRead(elem_areas); - const auto elem_verts_h = Omega_h::HostRead(elem_verts); - - double integral = 0.0; - for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { - const Omega_h::LO v0 = elem_verts_h[3 * e + 0]; - const Omega_h::LO v1 = elem_verts_h[3 * e + 1]; - const Omega_h::LO v2 = elem_verts_h[3 * e + 2]; - const double avg = (values[v0] + values[v1] + values[v2]) / 3.0; - integral += elem_areas_h[e] * avg; - } - return integral; -} - } // namespace // Fields that already live in the target order-1 space (constants and affine @@ -139,12 +42,12 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", // Source and target triangulate the same square along opposite diagonals. Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -166,8 +69,9 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-10)); } - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-10)); } SECTION("linear field is reproduced on target vertices and conserved") @@ -179,34 +83,31 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", const auto target_values = pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); - const auto tgt_coords_h = - Omega_h::HostRead(target_mesh.coords()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), + 2); for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { - const double expected = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1]; + const double expected = tgt_coords_h(i, 0) + tgt_coords_h(i, 1); REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-9)); } - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-9)); } } -// For a field that does not live in the target space the projection is not an -// interpolation, but conservation of the integral is the defining property of -// the conservative (Galerkin) projection and must still hold exactly. This -// also exercises a second Apply with a different source field to confirm the -// cached integrators/evaluator/factorization are reused without stale state. TEST_CASE("OmegaHConservativeProjection conserves the integral", "[transfer][mesh_intersection]") { Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -218,8 +119,9 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral", return x * x + x * y + 0.5 * y * y; }); projection.Apply(source, target); - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-12)); // Second Apply with a different source field — verifies cached state is // reused correctly and is not stale. @@ -227,8 +129,9 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral", source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return 2.0 * x - y + 0.5; }); projection.Apply(source, target); - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-12)); } TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local " @@ -238,13 +141,13 @@ TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local " Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); AddReorderedVertexGlobalIds(target_mesh); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh, "reordered_global"); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh, "reordered_global"); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -277,13 +180,13 @@ TEST_CASE("OmegaHConservativeProjection maps sparse target GIDs to active " Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); AddSparseVertexGlobalIds(target_mesh); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh, "sparse_global"); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh, "sparse_global"); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -318,12 +221,12 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - auto source_space = MakeP0Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_space = pcms::test::MakeP0Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -345,8 +248,9 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-10)); } - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP0Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP0Field(source_mesh, source)) + .margin(1e-10)); } SECTION("non-constant P0 source conserves the integral") @@ -357,8 +261,9 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", projection.Apply(source, target); - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP0Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP0Field(source_mesh, source)) + .margin(1e-10)); } } @@ -372,12 +277,12 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target", Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP0Space(target_mesh); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP0Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -399,8 +304,9 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target", for (Omega_h::LO e = 0; e < target_mesh.nelems(); ++e) { REQUIRE(target_values[e] == Catch::Approx(c).margin(1e-10)); } - REQUIRE(IntegrateP0Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-10)); } SECTION("linear field projects to cell average and conserves the integral") @@ -414,13 +320,16 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target", const auto target_values = pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); - const auto centroids = ElementCentroids(target_mesh); - REQUIRE(static_cast(target_values.size()) == centroids.size()); - for (std::size_t e = 0; e < centroids.size(); ++e) { - const double expected = f(centroids[e][0], centroids[e][1]); + const auto centroids_h = Omega_h::HostRead( + pcms::get_entity_centroids(target_mesh, target_mesh.dim())); + REQUIRE(static_cast(target_values.size()) == + target_mesh.nelems()); + for (Omega_h::LO e = 0; e < target_mesh.nelems(); ++e) { + const double expected = f(centroids_h[2 * e + 0], centroids_h[2 * e + 1]); REQUIRE(target_values[e] == Catch::Approx(expected).margin(1e-9)); } - REQUIRE(IntegrateP0Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9)); + REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-9)); } } diff --git a/test/test_omega_h_intersection_rhs_integrator.cpp b/test/test_omega_h_intersection_rhs_integrator.cpp index 18ea11bb..f5d667d9 100644 --- a/test/test_omega_h_intersection_rhs_integrator.cpp +++ b/test/test_omega_h_intersection_rhs_integrator.cpp @@ -21,37 +21,6 @@ namespace { -// Build a unit-square 2D simplex mesh. -// diagonal=0: T0=(0,1,3), T1=(1,2,3) -// diagonal=1: T0=(0,1,2), T1=(0,2,3) -Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) -{ - const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); - Omega_h::LOs ev2v = (diagonal == 0) ? Omega_h::LOs({0, 1, 3, 1, 2, 3}) - : Omega_h::LOs({0, 1, 2, 0, 2, 3}); - Omega_h::Mesh mesh(&lib); - Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); - mesh.add_tag( - 0, "class_dim", 1, - Omega_h::Read(mesh.nverts(), Omega_h::I8(0))); - mesh.add_tag( - 0, "class_id", 1, - Omega_h::Read(mesh.nverts(), Omega_h::ClassId(0))); - mesh.add_tag( - 1, "class_dim", 1, - Omega_h::Read(mesh.nedges(), Omega_h::I8(1))); - mesh.add_tag( - 1, "class_id", 1, - Omega_h::Read(mesh.nedges(), Omega_h::ClassId(0))); - mesh.add_tag( - 2, "class_dim", 1, - Omega_h::Read(mesh.nelems(), Omega_h::I8(2))); - mesh.add_tag( - 2, "class_id", 1, - Omega_h::Read(mesh.nelems(), Omega_h::ClassId(0))); - return mesh; -} - // Independent reference for the assembled conservative load vector // b_j = \int phi_j^target f dx // when f is exactly representable in the target order-1 space (constant or @@ -100,17 +69,7 @@ void CheckAssembledLoadMatches( const pcms::Field& source_field, const std::unordered_map& expected) { - const auto pts = integrator.GetIntegrationPoints(); - const std::size_t npts = pts.GetValues().extent(0); - - auto evaluator = source_space->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(pts)); - - Kokkos::View sampled("sampled", npts, - 1); - evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); - - integrator.Assemble(pcms::MakeConstRank2View(sampled)); + pcms::test::EvaluateAndAssemble(integrator, source_space, source_field); Vec vec = integrator.GetVector(); const PetscScalar* vec_array = nullptr; @@ -135,26 +94,19 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: integration points lie inside " "[rhs_integrator]") { Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto integrator = pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space); const auto raw_coords = integrator->GetIntegrationPoints().GetValues(); - auto raw_coords_view = Kokkos::View>( - raw_coords.data_handle(), raw_coords.extent(0), raw_coords.extent(1)); - auto raw_coords_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, raw_coords_view); + auto raw_coords_host = pcms::test::CopyCoordinatesToHost( + raw_coords, static_cast(raw_coords.extent(0)), + static_cast(raw_coords.extent(1))); const std::size_t n = raw_coords_host.extent(0); REQUIRE(n > 0); @@ -172,32 +124,18 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: zero source field gives zero " "[rhs_integrator]") { Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); // Default-constructed field has zero data. auto integrator = pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space); - const auto& pts = integrator->GetIntegrationPoints(); - const std::size_t npts = pts.GetValues().extent(0); - - auto evaluator = source_space->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(pts)); - - Kokkos::View sampled("sampled", npts, - 1); - evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); - - integrator->Assemble(pcms::MakeConstRank2View(sampled)); + pcms::test::EvaluateAndAssemble(*integrator, source_space, source_field); Vec vec = integrator->GetVector(); PetscReal norm = 0.0; @@ -214,15 +152,11 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: constant field matches " // For a constant source field the assembled load vector must equal the // target consistent mass matrix applied to the constant nodal vector. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); const double c = 2.0; auto source_field = source_space->CreateFunction(); @@ -251,15 +185,11 @@ TEST_CASE( // vector must equal the target consistent mass matrix applied to the nodal // values of x + y. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); pcms::test::SetField( @@ -268,11 +198,11 @@ TEST_CASE( const auto target_layout = std::dynamic_pointer_cast( target_space->GetLayout()); - const auto tgt_coords_h = - Omega_h::HostRead(target_mesh.coords()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), 2); std::vector g(target_mesh.nverts()); for (int i = 0; i < target_mesh.nverts(); ++i) { - g[i] = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1]; + g[i] = tgt_coords_h(i, 0) + tgt_coords_h(i, 1); } const auto expected = ExpectedLoadByGid(target_mesh, *target_layout, g); @@ -285,26 +215,22 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts", "[rhs_integrator]") { Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); SECTION("multi-component source space throws") { auto source_space = pcms::LagrangeFunctionSpace::FromMesh( source_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::test::MakeP1Space(target_mesh); REQUIRE_THROWS( pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space)); } SECTION("multi-component target space throws") { - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); auto target_space = pcms::LagrangeFunctionSpace::FromMesh( target_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); @@ -317,18 +243,14 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts", auto source_space = pcms::LagrangeFunctionSpace::FromMesh( source_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::test::MakeP1Space(target_mesh); REQUIRE_THROWS( pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space)); } SECTION("non-Cartesian target coordinate system throws") { - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); auto target_space = pcms::LagrangeFunctionSpace::FromMesh( target_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); diff --git a/test/test_omega_h_lagrange_field.cpp b/test/test_omega_h_lagrange_field.cpp index c17d45b1..0dd5ffc7 100644 --- a/test/test_omega_h_lagrange_field.cpp +++ b/test/test_omega_h_lagrange_field.cpp @@ -155,10 +155,10 @@ TEST_CASE("OmegaHLagrangeField order-1: linear function evaluation") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckEvaluation( factory, field, pcms::test::StandardEvalCoords2D(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } // The same linear evaluation test run on the MeshFields-backed order-1 field @@ -173,10 +173,10 @@ TEST_CASE("MeshFieldsAdapter order-1: linear function evaluation (shared util)") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckEvaluation( factory, field, pcms::test::StandardEvalCoords2D(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode") @@ -189,10 +189,10 @@ TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); Real fill_value = -999.0; - std::vector outside{-0.5, 0.5, 1.5, 0.5, 0.5, -0.5, 0.5, 1.5}; + auto outside = pcms::test::StandardOutsideCoords2D(); pcms::test::CheckFillMode(factory, field, fill_value, outside); } @@ -206,7 +206,7 @@ TEST_CASE("OmegaHLagrangeField order-1: serialize / deserialize round-trip") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckSerializeDeserialize(*factory->GetLayout(), field.GetData()); } @@ -270,24 +270,9 @@ TEST_CASE("OmegaHLagrangeField order-0: constant field evaluation") 1); field.GetData().SetDOFHolderDataHost(view); - auto pts = pcms::test::StandardEvalCoords2D(); - int n = static_cast(pts.size()) / 2; - auto device_coords = - pcms::test::CreateDeviceCoordinateView(pts, factory->GetCoordinateSystem()); - auto evaluator = factory->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - - Kokkos::View eval_device("eval", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto out = pcms::Rank2View( - eval_device.data(), n, 1); - evaluator->Evaluate(field, out); - auto eval_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), eval_device); - - for (int i = 0; i < n; ++i) - REQUIRE(eval_host(i) == Catch::Approx(kValue)); + pcms::test::CheckEvaluation( + factory, field, pcms::test::StandardEvalCoords2D(), + OMEGA_H_LAMBDA(Real, Real) { return kValue; }); } TEST_CASE("OmegaHLagrangeField order-0: out-of-bounds FILL mode") @@ -357,7 +342,8 @@ TEST_CASE("OmegaHLagrangeField: field valid after layout destruction") } // factory goes out of scope; field keeps layout alive pcms::test::SetField( - *field, OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + *field, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); // Just verify data was set correctly (no evaluator needed for this lifetime // test) auto data = field->GetDOFHolderDataHost(); diff --git a/test/test_omega_h_mass_integrator.cpp b/test/test_omega_h_mass_integrator.cpp index c9013f9d..c49a2e7d 100644 --- a/test/test_omega_h_mass_integrator.cpp +++ b/test/test_omega_h_mass_integrator.cpp @@ -8,6 +8,7 @@ #include #include #include +#include "field_test_utils.h" #include #include #include @@ -23,33 +24,6 @@ namespace { -Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib) -{ - const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); - Omega_h::LOs ev2v({0, 1, 3, 1, 2, 3}); - Omega_h::Mesh mesh(&lib); - Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); - mesh.add_tag( - 0, "class_dim", 1, - Omega_h::Read(mesh.nverts(), Omega_h::I8(0))); - mesh.add_tag( - 0, "class_id", 1, - Omega_h::Read(mesh.nverts(), Omega_h::ClassId(0))); - mesh.add_tag( - 1, "class_dim", 1, - Omega_h::Read(mesh.nedges(), Omega_h::I8(1))); - mesh.add_tag( - 1, "class_id", 1, - Omega_h::Read(mesh.nedges(), Omega_h::ClassId(0))); - mesh.add_tag( - 2, "class_dim", 1, - Omega_h::Read(mesh.nelems(), Omega_h::I8(2))); - mesh.add_tag( - 2, "class_id", 1, - Omega_h::Read(mesh.nelems(), Omega_h::ClassId(0))); - return mesh; -} - std::map, pcms::Real> BuildReferenceMassMap( Omega_h::Mesh& mesh, const pcms::FunctionSpace& space) { @@ -98,11 +72,9 @@ TEST_CASE( "[mass_integrator]") { Omega_h::Library lib; - auto mesh = BuildUnitSquare(lib); + auto mesh = pcms::test::BuildUnitSquare(lib, 0); - auto space = pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto space = pcms::test::MakeP1Space(mesh); const auto ref = BuildReferenceMassMap(mesh, *space); @@ -128,11 +100,9 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass", // each node: row_sum(i) = integral(N_i dx). For a regular mesh of unit // area with uniform nodal distribution the sum over all nodes equals 1. Omega_h::Library lib; - auto mesh = BuildUnitSquare(lib); + auto mesh = pcms::test::BuildUnitSquare(lib, 0); - auto space = pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto space = pcms::test::MakeP1Space(mesh); const auto layout = std::dynamic_pointer_cast( @@ -169,7 +139,7 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass", TEST_CASE("OmegaHMassIntegrator: rejects invalid layouts", "[mass_integrator]") { Omega_h::Library lib; - auto mesh = BuildUnitSquare(lib); + auto mesh = pcms::test::BuildUnitSquare(lib, 0); SECTION("multi-component space throws") { diff --git a/test/test_omega_h_mc_rhs_integrator.cpp b/test/test_omega_h_mc_rhs_integrator.cpp index 43172fce..4cf8a651 100644 --- a/test/test_omega_h_mc_rhs_integrator.cpp +++ b/test/test_omega_h_mc_rhs_integrator.cpp @@ -15,69 +15,6 @@ #include #include -// --------------------------------------------------------------------------- -// Helpers -// --------------------------------------------------------------------------- - -namespace -{ - -// Build a unit-square 2D simplex mesh. -// diagonal=0: T0=(0,1,3), T1=(1,2,3) -// diagonal=1: T0=(0,1,2), T1=(0,2,3) -Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) -{ - const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); - Omega_h::LOs ev2v = (diagonal == 0) ? Omega_h::LOs({0, 1, 3, 1, 2, 3}) - : Omega_h::LOs({0, 1, 2, 0, 2, 3}); - Omega_h::Mesh mesh(&lib); - Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); - mesh.add_tag( - 0, "class_dim", 1, - Omega_h::Read(mesh.nverts(), Omega_h::I8(0))); - mesh.add_tag( - 0, "class_id", 1, - Omega_h::Read(mesh.nverts(), Omega_h::ClassId(0))); - mesh.add_tag( - 1, "class_dim", 1, - Omega_h::Read(mesh.nedges(), Omega_h::I8(1))); - mesh.add_tag( - 1, "class_id", 1, - Omega_h::Read(mesh.nedges(), Omega_h::ClassId(0))); - mesh.add_tag( - 2, "class_dim", 1, - Omega_h::Read(mesh.nelems(), Omega_h::I8(2))); - mesh.add_tag( - 2, "class_id", 1, - Omega_h::Read(mesh.nelems(), Omega_h::ClassId(0))); - return mesh; -} - -std::shared_ptr MakeP1Space(Omega_h::Mesh& mesh) -{ - return pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); -} - -// Evaluates source_field at the integrator's sample points and assembles. -void EvaluateAndAssemble( - pcms::LinearFormIntegrator& integrator, - const std::shared_ptr& source_space, - const pcms::Field& source_field) -{ - const auto& pts = integrator.GetIntegrationPoints(); - const std::size_t npts = pts.GetValues().extent(0); - auto evaluator = source_space->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(pts)); - Kokkos::View sampled("sampled", npts, - 1); - evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); - integrator.Assemble(pcms::MakeConstRank2View(sampled)); -} - -} // namespace - // --------------------------------------------------------------------------- // Monte Carlo RHS integrator // --------------------------------------------------------------------------- @@ -86,20 +23,17 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: sample points lie inside the domain", "[mc_rhs_integrator]") { Omega_h::Library lib; - auto target_mesh = BuildUnitSquare(lib, 0); - auto target_space = MakeP1Space(target_mesh); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto target_space = pcms::test::MakeP1Space(target_mesh); const int samples_per_element = 16; for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) { pcms::OmegaHMonteCarloRHSIntegrator integrator( *target_space, samples_per_element, sampling); const auto raw_coords = integrator.GetIntegrationPoints().GetValues(); - auto coords_view = Kokkos::View>( - raw_coords.data_handle(), raw_coords.extent(0), raw_coords.extent(1)); - auto coords_h = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coords_view); + auto coords_h = pcms::test::CopyCoordinatesToHost( + raw_coords, static_cast(raw_coords.extent(0)), + static_cast(raw_coords.extent(1))); REQUIRE( coords_h.extent(0) == @@ -120,10 +54,10 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly", // For f = c the estimator sums to c * |domain| for any sample placement, // because the P1 basis functions partition unity at every sample point. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); pcms::test::SetField( @@ -131,7 +65,7 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly", for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) { pcms::OmegaHMonteCarloRHSIntegrator integrator(*target_space, 8, sampling); - EvaluateAndAssemble(integrator, source_space, source_field); + pcms::test::EvaluateAndAssemble(integrator, source_space, source_field); PetscScalar sum = 0.0; VecSum(integrator.GetVector(), &sum); @@ -152,10 +86,10 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target " // of the source field) equals the source field, so the sampled residual is // identically zero and the projection is exact regardless of sample count. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); auto target_field = target_space->CreateFunction(); @@ -170,11 +104,11 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target " const auto values = pcms::FlattenToRank1View(target_field.GetDOFHolderDataHost()); - const auto coords_h = Omega_h::HostRead(target_mesh.coords()); + const auto coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), 2); REQUIRE(static_cast(values.size()) == target_mesh.nverts()); for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { - const double expected = - 3.0 * coords_h[2 * i + 0] - coords_h[2 * i + 1] + 0.25; + const double expected = 3.0 * coords_h(i, 0) - coords_h(i, 1) + 0.25; CAPTURE(i, expected, values[i]); CHECK(values[i] == Catch::Approx(expected).margin(1e-9)); } @@ -188,10 +122,10 @@ TEST_CASE("OmegaHControlVariateProjection: reduces error vs plain Monte Carlo", // seeds the control-variate projection must be closer to the exact // (intersection-quadrature) projection than the plain Monte Carlo one. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); pcms::test::SetField( diff --git a/test/test_point_evaluator.cpp b/test/test_point_evaluator.cpp index dcc57087..e9e235a9 100644 --- a/test/test_point_evaluator.cpp +++ b/test/test_point_evaluator.cpp @@ -34,7 +34,7 @@ TEST_CASE("PointEvaluator: OmegaH order-1 linear evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; @@ -44,7 +44,7 @@ TEST_CASE("PointEvaluator: OmegaH order-1 linear evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } // ============================================================================ @@ -67,45 +67,24 @@ TEST_CASE("PointEvaluator: same evaluator reused for two FieldData objects") // field_a: linear_f; field_b: constant 42 pcms::test::SetField( field_a.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::SetField( field_b.GetData(), *factory->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); auto pts = pcms::test::StandardEvalCoords2D(); - int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); - // Create the PointEvaluator once + // Create the PointEvaluator once and reuse it for both fields. auto evaluator = factory->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_a_device("out_a", n); - Kokkos::View out_b_device("out_b", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto view_a = pcms::Rank2View( - out_a_device.data(), n, 1); - auto view_b = pcms::Rank2View( - out_b_device.data(), n, 1); - - // Evaluate field_a then field_b with the same evaluator - evaluator->Evaluate(field_a, view_a); - evaluator->Evaluate(field_b, view_b); - - auto out_a_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_a_device); - auto out_b_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_b_device); - - for (int i = 0; i < n; ++i) { - Real x = pts[2 * static_cast(i)], - y = pts[2 * static_cast(i) + 1]; - REQUIRE(out_a_host(i) == - Catch::Approx(pcms::test::linear_f(x, y)).margin(1e-10)); - REQUIRE(out_b_host(i) == Catch::Approx(42.0).margin(1e-10)); - } + pcms::test::CheckEvaluation( + *evaluator, field_a, pts, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); + pcms::test::CheckEvaluation( + *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); } // ============================================================================ @@ -125,11 +104,10 @@ TEST_CASE("PointEvaluator: OmegaH order-1 out-of-bounds fill") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); // Points clearly outside [0,1]^2 - const std::vector outside_pts = {-0.5, 0.5, 1.5, 0.5, - 0.5, -0.5, 0.5, 1.5}; + const auto outside_pts = pcms::test::StandardOutsideCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView( outside_pts, CoordinateSystem::Cartesian); pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0}; @@ -158,7 +136,7 @@ TEST_CASE("PointEvaluator: UniformGrid order-1 linear evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = @@ -167,7 +145,8 @@ TEST_CASE("PointEvaluator: UniformGrid order-1 linear evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 1e-8); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 1e-8); } TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation") @@ -184,7 +163,7 @@ TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = @@ -193,7 +172,8 @@ TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 1e-8); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 1e-8); } // ============================================================================ @@ -308,7 +288,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 linear evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; @@ -318,7 +298,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 linear evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill") @@ -334,10 +314,9 @@ TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); - const std::vector outside_pts = {-0.5, 0.5, 1.5, 0.5, - 0.5, -0.5, 0.5, 1.5}; + const auto outside_pts = pcms::test::StandardOutsideCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView( outside_pts, CoordinateSystem::Cartesian); pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0}; @@ -362,43 +341,24 @@ TEST_CASE( auto field_b = factory->CreateFunction(); pcms::test::SetField( field_a.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::SetField( field_b.GetData(), *factory->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); auto pts = pcms::test::StandardEvalCoords2D(); - int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); + // Create the PointEvaluator once and reuse it for both fields. auto evaluator = factory->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_a_device("out_a", n); - Kokkos::View out_b_device("out_b", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto view_a = pcms::Rank2View( - out_a_device.data(), n, 1); - auto view_b = pcms::Rank2View( - out_b_device.data(), n, 1); - - evaluator->Evaluate(field_a, view_a); - evaluator->Evaluate(field_b, view_b); - - auto out_a_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_a_device); - auto out_b_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_b_device); - - for (int i = 0; i < n; ++i) { - Real x = pts[2 * static_cast(i)], - y = pts[2 * static_cast(i) + 1]; - REQUIRE(out_a_host(i) == - Catch::Approx(pcms::test::linear_f(x, y)).margin(1e-10)); - REQUIRE(out_b_host(i) == Catch::Approx(42.0).margin(1e-10)); - } + pcms::test::CheckEvaluation( + *evaluator, field_a, pts, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); + pcms::test::CheckEvaluation( + *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); } TEST_CASE("LagrangeFunctionSpace: MeshFields rejects multi-component fields") diff --git a/test/test_polynomial_reconstruction_mls_evaluation.cpp b/test/test_polynomial_reconstruction_mls_evaluation.cpp index 888a83b1..415d4968 100644 --- a/test/test_polynomial_reconstruction_mls_evaluation.cpp +++ b/test/test_polynomial_reconstruction_mls_evaluation.cpp @@ -77,12 +77,6 @@ pcms::MLSOptions DefaultTestOptions3D() return opts; } -// Interior query points for a unit box — same as StandardEvalCoords2D. -std::vector QueryPoints() -{ - return pcms::test::StandardEvalCoords2D(); -} - pcms::MLSOptions SweepTestOptions(unsigned degree, pcms::RadialBasisFunction basis) { @@ -122,7 +116,7 @@ void CheckPolynomialReproduction(unsigned degree, auto field = fs->CreateFunction(); pcms::test::SetField(field.GetData(), *fs->GetLayout(), func); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); auto evaluator = fs->CreatePointEvaluator( @@ -149,7 +143,9 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: reproduces " CheckPolynomialReproduction( 0, basis, OMEGA_H_LAMBDA(Real, Real) { return 3.14; }, 5e-3); CheckPolynomialReproduction( - 1, basis, OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 5e-3); + 1, basis, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 5e-3); CheckPolynomialReproduction( 2, basis, OMEGA_H_LAMBDA(Real x, Real y) { return x * x + x * y + 2.0 * y * y; }, @@ -175,43 +171,26 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: same PointEvaluator " pcms::test::SetField( field_a.GetData(), *fs->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); const Real cval = 7.0; pcms::test::SetField( field_b.GetData(), *fs->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return cval; }); - auto pts = QueryPoints(); - int n = static_cast(pts.size()) / 2; + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); + // Create the PointEvaluator once and reuse it for both fields. auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_a_device("out_a", n); - Kokkos::View out_b_device("out_b", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View view_a(out_a_device.data(), - n, 1); - Rank2View view_b(out_b_device.data(), - n, 1); - - evaluator->Evaluate(field_a, view_a); - evaluator->Evaluate(field_b, view_b); - - auto out_a_host = - Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_a_device); - auto out_b_host = - Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_b_device); - - for (int i = 0; i < n; ++i) { - Real x = pts[2 * static_cast(i)], - y = pts[2 * static_cast(i) + 1]; - REQUIRE(out_a_host(i) == - Catch::Approx(pcms::test::linear_f(x, y)).margin(5e-3)); - REQUIRE(out_b_host(i) == Catch::Approx(cval).margin(5e-3)); - } + pcms::test::CheckEvaluation( + *evaluator, field_a, pts, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 5e-3); + pcms::test::CheckEvaluation( + *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return cval; }, + 5e-3); } // ============================================================================ @@ -228,7 +207,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: Evaluate throws for " coords_view, CoordinateSystem::Cartesian); auto field = fs->CreateFunction(); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); @@ -236,13 +215,8 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: Evaluate throws for " pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); // Two-component output — must throw - Kokkos::View out_device("out", - static_cast(n) * 2); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto out_view = - Rank2View(out_device.data(), n, 2); - REQUIRE_THROWS(evaluator->Evaluate(field, out_view)); + Kokkos::View out_device("out", n, 2); + REQUIRE_THROWS(evaluator->Evaluate(field, pcms::MakeRank2View(out_device))); } // ============================================================================ @@ -263,24 +237,20 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: default MLSOptions — " field.GetData(), *fs->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return Real(1.0); }); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_device("out", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto out_view = - Rank2View(out_device.data(), n, 1); + Kokkos::View out_device("out", n, 1); // Just verify it runs without error and returns finite values - REQUIRE_NOTHROW(evaluator->Evaluate(field, out_view)); + REQUIRE_NOTHROW(evaluator->Evaluate(field, pcms::MakeRank2View(out_device))); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); for (int i = 0; i < n; ++i) - REQUIRE(std::isfinite(out_host(i))); + REQUIRE(std::isfinite(out_host(i, 0))); } TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator " @@ -293,7 +263,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator " auto fs = pcms::PolynomialReconstructionFunctionSpace::Create( coords_view, CoordinateSystem::Cylindrical, DefaultTestOptions()); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); REQUIRE_THROWS(fs->CreatePointEvaluator( @@ -310,7 +280,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator " auto fs = pcms::PolynomialReconstructionFunctionSpace::Create( coords_view, CoordinateSystem::Cylindrical, DefaultTestOptions()); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cylindrical); REQUIRE_THROWS(fs->CreatePointEvaluator( @@ -342,16 +312,12 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: radius option is " pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_device("out", 1); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out_view(out_device.data(), - 1, 1); - evaluator->Evaluate(field, out_view); + Kokkos::View out_device("out", 1, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); - REQUIRE(out_host(0) == Catch::Approx(1.0).margin(1e-8)); + REQUIRE(out_host(0, 0) == Catch::Approx(1.0).margin(1e-8)); } TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve " @@ -375,34 +341,24 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve " } auto field = fs->CreateFunction(); - std::vector dof_values(27); - for (int i = 0; i < 27; ++i) { - const Real x = src[3 * i + 0]; - const Real y = src[3 * i + 1]; - const Real z = src[3 * i + 2]; - dof_values[i] = x + 2.0 * y + 3.0 * z; - } - field.GetData().SetDOFHolderDataHost(Rank2View( - dof_values.data(), static_cast(dof_values.size()), 1)); + pcms::test::SetField( + field, + KOKKOS_LAMBDA(Real x, Real y, Real z) { return x + 2.0 * y + 3.0 * z; }); std::vector pts{0.5, 0.5, 0.25, 0.5, 0.5, 0.75}; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian, 3); auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_device("out", 2); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out_view(out_device.data(), - 2, 1); - evaluator->Evaluate(field, out_view); + Kokkos::View out_device("out", 2, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); // This 3D case is a small, low-resolution support cloud with MLS weights // built from a finite-radius neighborhood rather than exact nodal lookup, so // it is checked with a slightly looser tolerance than the denser 2D tests. - REQUIRE(out_host(0) == Catch::Approx(2.25).margin(1e-2)); - REQUIRE(out_host(1) == Catch::Approx(3.75).margin(1e-2)); - REQUIRE(out_host(1) - out_host(0) == Catch::Approx(1.5).margin(1e-2)); + REQUIRE(out_host(0, 0) == Catch::Approx(2.25).margin(1e-2)); + REQUIRE(out_host(1, 0) == Catch::Approx(3.75).margin(1e-2)); + REQUIRE(out_host(1, 0) - out_host(0, 0) == Catch::Approx(1.5).margin(1e-2)); } diff --git a/test/test_spr_meshfields.cpp b/test/test_spr_meshfields.cpp index 42a6acf0..7ad4470a 100644 --- a/test/test_spr_meshfields.cpp +++ b/test/test_spr_meshfields.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -127,22 +128,7 @@ TEST_CASE("meshfields_spr_test") const auto& ntargets = mesh.nverts(); - Write source_coordinates( - dim * nfaces, 0, "stores coordinates of cell centroid of each tri element"); - - const auto& faces2nodes = mesh.ask_down(FACE, VERT).ab2b; - - Kokkos::parallel_for( - "calculate the centroid in each tri element", nfaces, - OMEGA_H_LAMBDA(const LO id) { - const auto current_el_verts = gather_verts<3>(faces2nodes, id); - const Omega_h::Few, 3> current_el_vert_coords = - gather_vectors<3, 2>(target_coordinates, current_el_verts); - auto centroid = average(current_el_vert_coords); - int index = 2 * id; - source_coordinates[index] = centroid[0]; - source_coordinates[index + 1] = centroid[1]; - }); + const auto source_coordinates = pcms::get_entity_centroids(mesh, FACE); pcms::Points source_points; diff --git a/test/test_svd_serial.cpp b/test/test_svd_serial.cpp index bdc23290..1774e29b 100644 --- a/test/test_svd_serial.cpp +++ b/test/test_svd_serial.cpp @@ -38,7 +38,7 @@ TEST_CASE("test_serial_svd") Kokkos::deep_copy(A_data, host_A_data); - Kokkos::View rhs_data("Device rhs data", row, column); + Kokkos::View rhs_data("Device rhs data", row); auto host_rhs_data = Kokkos::create_mirror_view(rhs_data); host_rhs_data(0) = 1.28571; @@ -54,7 +54,7 @@ TEST_CASE("test_serial_svd") { Kokkos::View result("result", row, row); - Kokkos::View transpose_expected("result", row, row); + Kokkos::View transpose_expected("transpose_expected", row, row); Kokkos::deep_copy(result, 0.0); Kokkos::deep_copy(transpose_expected, 0.0); team_policy tp(1, Kokkos::AUTO); diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp index 632bd4db..699ec26a 100644 --- a/test/test_uniform_grid_field.cpp +++ b/test/test_uniform_grid_field.cpp @@ -72,7 +72,7 @@ void VerifyUniformGridFieldValues( int vertex_id = j * (grid.divisions[0] + 1) + i; pcms::Real x = ug_coords.GetValues()(vertex_id, 0); pcms::Real y = ug_coords.GetValues()(vertex_id, 1); - pcms::Real expected = x + 2.0 * y; + pcms::Real expected = pcms::test::linear_f(x, y); pcms::Real actual = ug_field_data[vertex_id]; REQUIRE(std::abs(expected - actual) <= 1e-10); } @@ -151,21 +151,16 @@ TEST_CASE("UniformGrid order-0 field creation and evaluation") auto evaluator = eval_factory.CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View results_host("results_host", - 4); - Kokkos::View results_device( - "results_device", 4); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View out( - results_device.data(), 4, 1); - evaluator->Evaluate(field, out); - Kokkos::deep_copy(results_host, results_device); - - REQUIRE(results_host(0) == Catch::Approx(1.0)); - REQUIRE(results_host(1) == Catch::Approx(2.0)); - REQUIRE(results_host(2) == Catch::Approx(3.0)); - REQUIRE(results_host(3) == Catch::Approx(4.0)); + Kokkos::View results_device( + "results_device", 4, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(results_device)); + auto results_host = Kokkos::create_mirror_view_and_copy( + pcms::HostMemorySpace(), results_device); + + REQUIRE(results_host(0, 0) == Catch::Approx(1.0)); + REQUIRE(results_host(1, 0) == Catch::Approx(2.0)); + REQUIRE(results_host(2, 0) == Catch::Approx(3.0)); + REQUIRE(results_host(3, 0) == Catch::Approx(4.0)); } TEST_CASE("UniformGrid field data operations", "[uniform_grid_field]") @@ -236,26 +231,21 @@ TEST_CASE("UniformGrid field evaluation - piecewise constant") auto evaluator = eval_factory.CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View results_host("results_host", - 4); - Kokkos::View results_device( - "results_device", 4); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View out( - results_device.data(), 4, 1); - evaluator->Evaluate(field, out); - Kokkos::deep_copy(results_host, results_device); + Kokkos::View results_device( + "results_device", 4, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(results_device)); + auto results_host = Kokkos::create_mirror_view_and_copy( + pcms::HostMemorySpace(), results_device); // Check results - interpolated from vertices // Cell 0 center (2.5, 2.5): avg of v0,v1,v3,v4 = (1.0+1.5+2.0+2.5)/4 = 1.75 // Cell 1 center (7.5, 2.5): avg of v1,v2,v4,v5 = (1.5+2.0+2.5+3.0)/4 = 2.25 // Cell 2 center (2.5, 7.5): avg of v3,v4,v6,v7 = (2.0+2.5+3.0+3.5)/4 = 2.75 // Cell 3 center (7.5, 7.5): avg of v4,v5,v7,v8 = (2.5+3.0+3.5+4.0)/4 = 3.25 - REQUIRE(std::abs(results_host(0) - 1.75) < 1e-10); - REQUIRE(std::abs(results_host(1) - 2.25) < 1e-10); - REQUIRE(std::abs(results_host(2) - 2.75) < 1e-10); - REQUIRE(std::abs(results_host(3) - 3.25) < 1e-10); + REQUIRE(std::abs(results_host(0, 0) - 1.75) < 1e-10); + REQUIRE(std::abs(results_host(1, 0) - 2.25) < 1e-10); + REQUIRE(std::abs(results_host(2, 0) - 2.75) < 1e-10); + REQUIRE(std::abs(results_host(3, 0) - 3.25) < 1e-10); } TEST_CASE("UniformGrid field serialization") @@ -328,8 +318,9 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") mesh, 1, 1, pcms::CoordinateSystem::Cartesian); auto omega_h_field = omega_h_factory->CreateFunction(); pcms::test::SetField( - omega_h_field, - OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + omega_h_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { + return pcms::test::linear_f(x, y); + }); pcms::UniformGrid<2> grid; grid.edge_length = {1.0, 1.0}; @@ -354,7 +345,7 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") for (int i = 0; i < num_ug_nodes; ++i) { pcms::Real x = ug_coords_host(i, 0); pcms::Real y = ug_coords_host(i, 1); - pcms::Real expected = x + 2.0 * y; + pcms::Real expected = pcms::test::linear_f(x, y); REQUIRE(std::abs(transferred_data[i] - expected) < 1e-6); } } @@ -641,8 +632,9 @@ TEST_CASE("UniformGrid workflow") mesh, 1, 1, pcms::CoordinateSystem::Cartesian); auto omega_h_field = omega_h_factory->CreateFunction(); pcms::test::SetField( - omega_h_field, - OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + omega_h_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { + return pcms::test::linear_f(x, y); + }); auto ug_factory = pcms::LagrangeFunctionSpace::FromUniformGrid( grid, 1, pcms::CoordinateSystem::Cartesian); @@ -658,23 +650,11 @@ TEST_CASE("UniformGrid workflow") auto ug_coords_host_view = pcms::test::CopyCoordinatesToHost(ug_coords_device_view, 25, 2); - auto ug_field_data_device = ug_field.GetDOFHolderData(); - Kokkos::View ug_field_data_device_view( - "", 25); - Kokkos::parallel_for( - "CopyFieldDataToView", 25, KOKKOS_LAMBDA(int i) { - ug_field_data_device_view(i) = ug_field_data_device(i, 0); - }); - auto ug_field_data_host_view = - Kokkos::View("", 25); - Kokkos::deep_copy(ug_field_data_host_view, ug_field_data_device_view); - - pcms::Rank2View ug_coords( - ug_coords_host_view.data(), 25, 2); pcms::CoordinateView ug_coords_view( - pcms::CoordinateSystem::Cartesian, ug_coords); - pcms::Rank1View ug_field_data( - ug_field_data_host_view.data(), 25); + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(ug_coords_host_view)); + const auto ug_field_data = + pcms::FlattenToRank1View(ug_field.GetDOFHolderDataHost()); VerifyUniformGridFieldValues(grid, ug_coords_view, ug_field_data); VerifyMaskFieldValues(grid, mask_field); diff --git a/test/test_xgc_reverse_classification.cpp b/test/test_xgc_reverse_classification.cpp index 4b87fcb9..5603b52b 100644 --- a/test/test_xgc_reverse_classification.cpp +++ b/test/test_xgc_reverse_classification.cpp @@ -39,7 +39,7 @@ TEST_CASE("reverse classification") } auto vec = rc.Serialize(); pcms::ReverseClassificationVertex rc_deserialized; - pcms::Rank1View av{vec.data(), vec.size()}; + auto av = pcms::make_array_view(vec); rc_deserialized.Deserialize(av); REQUIRE(rc_deserialized == rc); } From ae79ccbf1647a6ecde2a71429088d3487ffbf3a0 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Wed, 15 Jul 2026 21:35:08 -0400 Subject: [PATCH 2/4] Add missing deep_copy in omega_h_lagrange. --- src/pcms/field/layout/omega_h_lagrange.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/pcms/field/layout/omega_h_lagrange.cpp b/src/pcms/field/layout/omega_h_lagrange.cpp index a325140f..c9d7bdb9 100644 --- a/src/pcms/field/layout/omega_h_lagrange.cpp +++ b/src/pcms/field/layout/omega_h_lagrange.cpp @@ -171,6 +171,7 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout( owned_ = BuildOwned(mesh_, entity_dim, owned_mask); owned_host_ = Kokkos::View("owned_host", owned_.size()); + Kokkos::deep_copy(owned_host_, owned_); class_ids_ = Omega_h::Read( mesh_.get_array(entity_dim, "class_id")); From cf1018726a10e4d65a1eec7de33151e62be28276 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Thu, 6 Aug 2026 09:37:32 -0400 Subject: [PATCH 3/4] transfer: support 3D tetrahedral meshes in conservative transfer Generalize the intersection-based and Monte Carlo transfer machinery from 2D triangles to linear simplices in both 2D and 3D: - get_vert_coords_of_elem and adjBasedIntersectSearch are templated on the spatial dimension and gather Dim+1 vertices per element - the form integrator utilities compute reference-element mappings, barycentric coordinates and quadrature on triangles or tetrahedra - MassMatrixIntegrator sizes its element blocks from FieldElement::MeshEntDim + 1 rather than a hard coded 3, and integrates against the unsigned volume element so tets with negative orientation contribute positively - the MC RHS integrator carries the number of barycentric DOFs per sample instead of assuming 3 The BFS ring buffers in queue_visited.hpp are enlarged for 3D, where a target tet can overlap many more source elements than a triangle, and their capacities are now exposed as the PCMS_INTERSECTION_QUEUE_SIZE and PCMS_INTERSECTION_TRACK_SIZE cache variables so device builds can lower the per-thread stack footprint. --- CMakeLists.txt | 9 + src/pcms/configuration.h.in | 3 + src/pcms/localization/queue_visited.hpp | 19 +- src/pcms/transfer/mass_matrix_integrator.hpp | 19 +- src/pcms/transfer/mesh_intersection.cpp | 126 ++++++-- src/pcms/transfer/mesh_intersection.hpp | 34 +- .../omega_h_form_integrator_utils.hpp | 292 ++++++++++++++---- .../omega_h_intersection_rhs_integrator.cpp | 83 +++-- src/pcms/transfer/omega_h_mass_integrator.cpp | 159 ++++++---- .../transfer/omega_h_mc_rhs_integrator.cpp | 143 ++++++--- .../transfer/omega_h_mc_rhs_integrator.hpp | 4 +- 11 files changed, 636 insertions(+), 255 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d407a53..40be9ed6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,6 +49,15 @@ option(PCMS_ENABLE_Python "Enable pcms Python api" OFF) option(PCMS_ENABLE_PRINT "PCMS print statements enabled" ON) +# Per-thread on-stack ring-buffer capacities for the adjacency-based +# intersection/adjacency BFS (queue_visited.hpp). Defaults are sized for 3D +# tetrahedral meshes on the serial/OpenMP backends. Lower them for device +# (GPU) builds, where each search thread's stack must hold these arrays. +set(PCMS_INTERSECTION_QUEUE_SIZE 1024 CACHE STRING + "Max BFS queue size per target element in the intersection search") +set(PCMS_INTERSECTION_TRACK_SIZE 2048 CACHE STRING + "Max visited-set size per target element in the intersection search") + option(PETSC_LINK_STATIC "Use pkg-config --static results for PETSc" ${_pcms_link_petsc_static_default}) diff --git a/src/pcms/configuration.h.in b/src/pcms/configuration.h.in index 51859ec0..24585146 100644 --- a/src/pcms/configuration.h.in +++ b/src/pcms/configuration.h.in @@ -7,3 +7,6 @@ #cmakedefine PCMS_ENABLE_Fortran #cmakedefine PCMS_ENABLE_MESHFIELDS #cmakedefine PCMS_ENABLE_PETSC + +#cmakedefine PCMS_INTERSECTION_QUEUE_SIZE @PCMS_INTERSECTION_QUEUE_SIZE@ +#cmakedefine PCMS_INTERSECTION_TRACK_SIZE @PCMS_INTERSECTION_TRACK_SIZE@ diff --git a/src/pcms/localization/queue_visited.hpp b/src/pcms/localization/queue_visited.hpp index 8d768eec..3db919e8 100644 --- a/src/pcms/localization/queue_visited.hpp +++ b/src/pcms/localization/queue_visited.hpp @@ -7,8 +7,23 @@ #include #include #include -#define MAX_SIZE_QUEUE 500 -#define MAX_SIZE_TRACK 800 +#include "pcms/configuration.h" +// Per-thread BFS ring-buffer capacities for the adjacency-based intersection +// search. Sized for 3D: a target tetrahedron can overlap many more source +// elements than a 2D triangle, and the visited set additionally holds the +// non-intersecting neighbors probed along the way. These are on-stack arrays, +// so raising them raises per-thread stack usage; lower them for device (GPU) +// builds via the PCMS_INTERSECTION_{QUEUE,TRACK}_SIZE CMake cache variables, +// which flow in through pcms/configuration.h. The fallbacks below apply only if +// the configured header is unavailable (e.g. a header-only consumer). +#ifndef PCMS_INTERSECTION_QUEUE_SIZE +#define PCMS_INTERSECTION_QUEUE_SIZE 1024 +#endif +#ifndef PCMS_INTERSECTION_TRACK_SIZE +#define PCMS_INTERSECTION_TRACK_SIZE 2048 +#endif +#define MAX_SIZE_QUEUE PCMS_INTERSECTION_QUEUE_SIZE +#define MAX_SIZE_TRACK PCMS_INTERSECTION_TRACK_SIZE namespace pcms { diff --git a/src/pcms/transfer/mass_matrix_integrator.hpp b/src/pcms/transfer/mass_matrix_integrator.hpp index 7d36c975..8eabf475 100644 --- a/src/pcms/transfer/mass_matrix_integrator.hpp +++ b/src/pcms/transfer/mass_matrix_integrator.hpp @@ -17,16 +17,19 @@ template class MassMatrixIntegrator : public MeshField::Integrator { public: + // Linear simplex: numNodes = spatial dim + 1 (3 for triangles, 4 for tets). + static constexpr int numNodes = FieldElement::MeshEntDim + 1; + MassMatrixIntegrator(Omega_h::Mesh& mesh_in, FieldElement& fe_in, int order = 2) : mesh(mesh_in), fe(fe_in), - subMatrixSize(3 * 3), // FIXME remove hard coded size - elmMassMatrix("elmMassMatrix", mesh_in.nelems() * 3 * 3), + subMatrixSize(numNodes * numNodes), + elmMassMatrix("elmMassMatrix", mesh_in.nelems() * numNodes * numNodes), Integrator(order) { Kokkos::deep_copy(elmMassMatrix, 0); - assert(mesh.dim() == 2); // TODO support 1d,2d,3d + assert(mesh.dim() == 2 || mesh.dim() == 3); assert(mesh.family() == OMEGA_H_SIMPLEX); } void atPoints(Kokkos::View p, @@ -58,12 +61,18 @@ class MassMatrixIntegrator : public MeshField::Integrator } const auto N = shapeFn.getValues(localCoord); const auto wPt = w(pt); - const auto dVPt = dV(pt); + // Use the unsigned volume element: MeshField returns a signed + // Jacobian determinant, which is negative for tetrahedra whose vertex + // ordering has negative orientation. A mass matrix integrates against + // the positive volume measure, so take the magnitude (a no-op in 2D + // where the differential area is already positive). + const auto dVPt = Kokkos::fabs(dV(pt)); // printf("Shape Functions: %f, %f, %f \n", N[0], N[1], N[2]); // printf("wPt, dVPt: %f, %f \n", wPt, dVPt); for (auto i = 0; i < N.size(); i++) { for (auto j = 0; j < N.size(); j++) { - massMatrix(elm * subMat + i * 3 + j) += N[i] * N[j] * wPt * dVPt; + massMatrix(elm * subMat + i * numNodes + j) += + N[i] * N[j] * wPt * dVPt; } } } diff --git a/src/pcms/transfer/mesh_intersection.cpp b/src/pcms/transfer/mesh_intersection.cpp index 564ae350..06f26d20 100644 --- a/src/pcms/transfer/mesh_intersection.cpp +++ b/src/pcms/transfer/mesh_intersection.cpp @@ -4,51 +4,66 @@ namespace pcms { +namespace +{ +// Construct the source-mesh containing-element search appropriate to the +// spatial dimension: a uniform 20^Dim background grid over the source mesh. +template +auto MakeGridPointSearch(Omega_h::Mesh& source_mesh) +{ + if constexpr (Dim == 3) { + return pcms::GridPointSearch3D(source_mesh, 20, 20, 20); + } else { + return pcms::GridPointSearch2D(source_mesh, 20, 20); + } +} +} // namespace + +template void FindIntersections::adjBasedIntersectSearch( const Omega_h::LOs& tgt2src_offsets, Omega_h::Write& nIntersections, Omega_h::Write& tgt2src_indices, bool is_count_only) { - + // Element entity dimension equals the spatial dimension (FACE for 2D, REGION + // for 3D); measures are triangle areas (2D) or tet volumes (3D). const auto& tgt_coords = target_mesh_.coords(); const auto& src_coords = source_mesh_.coords(); - const auto& tgt_faces2nodes = - target_mesh_.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; - const auto& src_faces2nodes = - source_mesh_.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; - const auto& src_elem_areas = measure_elements_real(&source_mesh_); - const auto& tgt_elem_areas = measure_elements_real(&target_mesh_); + const auto& tgt_elems2nodes = target_mesh_.ask_down(Dim, Omega_h::VERT).ab2b; + const auto& src_elems2nodes = source_mesh_.ask_down(Dim, Omega_h::VERT).ab2b; + const auto& src_elem_measures = measure_elements_real(&source_mesh_); + const auto& tgt_elem_measures = measure_elements_real(&target_mesh_); const auto& t2t = source_mesh_.ask_dual(); // gives connected element neighbors const auto& t2tt = t2t.a2ab; const auto& tt2t = t2t.ab2b; - const auto flat_centroids = - pcms::get_entity_centroids(target_mesh_, Omega_h::FACE); + const auto flat_centroids = pcms::get_entity_centroids(target_mesh_, Dim); // Convert layout_right 1D Omega_h array to 2D Kokkos view with correct layout - auto centroids = ConvertCoordsTo2D(flat_centroids, target_mesh_.nfaces(), 2); + auto centroids = + ConvertCoordsTo2D(flat_centroids, target_mesh_.nelems(), Dim); - pcms::GridPointSearch2D search_cell(source_mesh_, 20, 20); + auto search_cell = MakeGridPointSearch(source_mesh_); auto results = search_cell(centroids); auto owning_cell_ids = search_cell.GetOwningElementIds(results); - auto nfaces_target = target_mesh_.nfaces(); + auto nelems_target = target_mesh_.nelems(); Omega_h::parallel_for( - nfaces_target, + nelems_target, OMEGA_H_LAMBDA(const Omega_h::LO id) { Queue queue; Track visited; auto current_cell_id = owning_cell_ids(id); - auto current_tgt_elm_area = tgt_elem_areas[id]; + auto current_tgt_elm_measure = tgt_elem_measures[id]; OMEGA_H_CHECK_PRINTF(current_cell_id >= 0, "ERROR: source cell id not found for given target " - "centroid %d (%f, %f)\n", - id, centroids(id, 0), centroids(id, 1)); + "centroid %d\n", + id); auto tgt_elm_vert_coords = - get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, id); + get_vert_coords_of_elem(tgt_coords, tgt_elems2nodes, id); Omega_h::LO start_counter; if (!is_count_only) { @@ -76,22 +91,42 @@ void FindIntersections::adjBasedIntersectSearch( auto neighborElmId = tt2t[i]; if (visited.notVisited(neighborElmId)) { - visited.push_back(neighborElmId); - auto elm_vert_coords = get_vert_coords_of_elem( - src_coords, src_faces2nodes, neighborElmId); - r3d::Polytope<2> intersection; + // If the visited buffer is full, skip this neighbor so the BFS + // terminates. Without this, an unrecorded neighbor stays "not + // visited" and is re-queued forever (infinite loop). Mirrors the + // guard in adj_search.cpp. Raise MAX_SIZE_TRACK if this fires. + if (!visited.push_back(neighborElmId)) { + printf("ERROR: visited buffer full (MAX_SIZE_TRACK=%d) for " + "target %d; some intersections may be missed\n", + MAX_SIZE_TRACK, id); + continue; + } + auto elm_vert_coords = get_vert_coords_of_elem( + src_coords, src_elems2nodes, neighborElmId); + r3d::Polytope intersection; r3d::intersect_simplices(intersection, tgt_elm_vert_coords, elm_vert_coords); - auto intersected_area = r3d::measure(intersection); - auto current_src_elm_area = src_elem_areas[neighborElmId]; + // Take the magnitude: r3d::measure is signed by the orientation of + // the target simplex used to initialize the polytope, which can be + // negative for tetrahedra. This mirrors the fabs applied in the + // sub-simplex decomposition and mass assembly; without it a + // negatively-oriented target element would reject all of its real + // overlaps and break conservation. + auto intersected_measure = Kokkos::fabs(r3d::measure(intersection)); + auto current_src_elm_measure = src_elem_measures[neighborElmId]; auto scale = - Kokkos::fmax(current_tgt_elm_area, current_src_elm_area); + Kokkos::fmax(current_tgt_elm_measure, current_src_elm_measure); auto eps = Kokkos::fmax(abs_tol, rel_tol * scale); - if (intersection.nverts >= 3 && intersected_area >= eps) { + // A valid intersection is a non-degenerate simplex-simplex overlap: + // at least Dim+1 vertices (a polygon in 2D, a polyhedron in 3D). + if (intersection.nverts >= Dim + 1 && intersected_measure >= eps) { count++; OMEGA_H_CHECK_PRINTF( - count < 500, "WARNING: count exceeds 500 for target %d", id); + count < kMaxIntersectionsPerTarget, + "intersection count for target %d reached the cap %d; raise " + "kMaxIntersectionsPerTarget/MAX_SIZE_QUEUE", + id, kMaxIntersectionsPerTarget); queue.push_back(neighborElmId); @@ -113,20 +148,32 @@ void FindIntersections::adjBasedIntersectSearch( }, // end of lambda "count the number of intersections for each target element"); } -IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh, - Omega_h::Mesh& target_mesh) + +// Explicit instantiations for the supported spatial dimensions. +template void FindIntersections::adjBasedIntersectSearch<2>( + const Omega_h::LOs&, Omega_h::Write&, + Omega_h::Write&, bool); +template void FindIntersections::adjBasedIntersectSearch<3>( + const Omega_h::LOs&, Omega_h::Write&, + Omega_h::Write&, bool); + +namespace +{ +template +IntersectionResults intersectTargetsImpl(Omega_h::Mesh& source_mesh, + Omega_h::Mesh& target_mesh) { FindIntersections intersect(source_mesh, target_mesh); - auto nfaces_target = target_mesh.nfaces(); + auto nelems_target = target_mesh.nelems(); Omega_h::Write nIntersections( - nfaces_target, 0, "number of intersections in each target vertex"); + nelems_target, 0, "number of intersections in each target element"); Omega_h::Write tgt2src_indices; - intersect.adjBasedIntersectSearch(Omega_h::LOs(), nIntersections, - tgt2src_indices, true); + intersect.adjBasedIntersectSearch(Omega_h::LOs(), nIntersections, + tgt2src_indices, true); Kokkos::fence(); auto tgt2src_offsets = Omega_h::offset_scan(Omega_h::Read(nIntersections), @@ -139,9 +186,20 @@ IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh, ntotal_intersections, 0, "indices of the source elements that intersect the given target element"); - intersect.adjBasedIntersectSearch(tgt2src_offsets, nIntersections, - tgt2src_indices, false); + intersect.adjBasedIntersectSearch(tgt2src_offsets, nIntersections, + tgt2src_indices, false); return {.tgt2src_offsets = tgt2src_offsets, .tgt2src_indices = Omega_h::read(tgt2src_indices)}; } +} // namespace + +IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh, + Omega_h::Mesh& target_mesh) +{ + OMEGA_H_CHECK(source_mesh.dim() == target_mesh.dim()); + if (source_mesh.dim() == 3) { + return intersectTargetsImpl<3>(source_mesh, target_mesh); + } + return intersectTargetsImpl<2>(source_mesh, target_mesh); +} } // namespace pcms diff --git a/src/pcms/transfer/mesh_intersection.hpp b/src/pcms/transfer/mesh_intersection.hpp index 84a52114..f72b5c86 100644 --- a/src/pcms/transfer/mesh_intersection.hpp +++ b/src/pcms/transfer/mesh_intersection.hpp @@ -15,19 +15,29 @@ namespace pcms constexpr static double abs_tol = 1e-18; /// abs tolerance constexpr static double rel_tol = 1e-12; /// rel tolerance -[[nodiscard]] OMEGA_H_INLINE r3d::Few, 3> +// Upper bound on intersecting source elements recorded per target element. +// Must not exceed the BFS queue capacity (MAX_SIZE_QUEUE); raised for 3D where +// a single target tet can overlap many source tets. +constexpr static int kMaxIntersectionsPerTarget = MAX_SIZE_QUEUE; + +// Gather the vertex coordinates of a simplex element (triangle for Dim==2, +// tetrahedron for Dim==3) into an r3d simplex, ready for +// r3d::intersect_simplices. +template +[[nodiscard]] OMEGA_H_INLINE r3d::Few, Dim + 1> get_vert_coords_of_elem(const Omega_h::Reals& coords, - const Omega_h::LOs& faces2nodes, const int id) + const Omega_h::LOs& elems2nodes, const int id) { - const auto elm_verts = Omega_h::gather_verts<3>(faces2nodes, id); + const auto elm_verts = Omega_h::gather_verts(elems2nodes, id); - const Omega_h::Matrix<2, 3> elm_vert_coords = - Omega_h::gather_vectors<3, 2>(coords, elm_verts); + const Omega_h::Matrix elm_vert_coords = + Omega_h::gather_vectors(coords, elm_verts); - r3d::Few, 3> r3d_vector; - for (int i = 0; i < 3; ++i) { - r3d_vector[i][0] = elm_vert_coords[i][0]; - r3d_vector[i][1] = elm_vert_coords[i][1]; + r3d::Few, Dim + 1> r3d_vector; + for (int i = 0; i < Dim + 1; ++i) { + for (int d = 0; d < Dim; ++d) { + r3d_vector[i][d] = elm_vert_coords[i][d]; + } } return r3d_vector; @@ -78,11 +88,13 @@ class FindIntersections * @param is_count_only If true, only counts intersections; if false, also * fills tgt2src_indices. * - * @note This method assumes 2D linear triangles and uses - * `r3d::intersect_simplices` for geometric intersection. + * @note Templated on spatial dimension `Dim`: linear triangles (Dim==2) or + * linear tetrahedra (Dim==3), using `r3d::intersect_simplices` for geometric + * intersection. * * @see r3d::intersect_simplices, intersectTargets */ + template void adjBasedIntersectSearch(const Omega_h::LOs& tgt2src_offsets, Omega_h::Write& nIntersections, Omega_h::Write& tgt2src_indices, diff --git a/src/pcms/transfer/omega_h_form_integrator_utils.hpp b/src/pcms/transfer/omega_h_form_integrator_utils.hpp index 23f5e35d..ee2b5194 100644 --- a/src/pcms/transfer/omega_h_form_integrator_utils.hpp +++ b/src/pcms/transfer/omega_h_form_integrator_utils.hpp @@ -13,9 +13,10 @@ namespace pcms::detail { -// Shared checks for a scalar Cartesian Lagrange space on a 2D simplex mesh, -// independent of order. Order is validated separately by the callers below. -inline void CheckOmegaHScalarSimplex2DLayout( +// Shared checks for a scalar Cartesian Lagrange space on a simplex mesh +// (triangles in 2D, tetrahedra in 3D), independent of order. Order is validated +// separately by the callers below. +inline void CheckOmegaHScalarSimplexLayout( CoordinateSystem coordinate_system, const std::shared_ptr& layout, const char* context, const char* role) @@ -33,12 +34,13 @@ inline void CheckOmegaHScalarSimplex2DLayout( " space must use Cartesian coordinates"); } const Omega_h::Mesh& mesh = layout->GetMesh(); - if (mesh.dim() != 2) { - throw pcms_error(std::string(context) + ": " + role + " mesh must be 2D"); + if (mesh.dim() != 2 && mesh.dim() != 3) { + throw pcms_error(std::string(context) + ": " + role + + " mesh must be 2D or 3D"); } if (mesh.family() != OMEGA_H_SIMPLEX) { throw pcms_error(std::string(context) + ": " + role + - " mesh must be a simplex (triangle) mesh"); + " mesh must be a simplex (triangle/tetrahedron) mesh"); } } @@ -48,7 +50,7 @@ inline void CheckOmegaHScalarP1Layout( const std::shared_ptr& layout, const char* context, const char* role) { - CheckOmegaHScalarSimplex2DLayout(coordinate_system, layout, context, role); + CheckOmegaHScalarSimplexLayout(coordinate_system, layout, context, role); if (layout->GetOrder() != 1) { throw pcms_error(std::string(context) + ": " + role + " space must be order-1"); @@ -63,7 +65,7 @@ inline void CheckOmegaHScalarLagrangeLayout( const std::shared_ptr& layout, const char* context, const char* role) { - CheckOmegaHScalarSimplex2DLayout(coordinate_system, layout, context, role); + CheckOmegaHScalarSimplexLayout(coordinate_system, layout, context, role); const int order = layout->GetOrder(); if (order != 0 && order != 1) { throw pcms_error(std::string(context) + ": " + role + @@ -71,28 +73,23 @@ inline void CheckOmegaHScalarLagrangeLayout( } } -[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<2> GlobalFromBarycentric( - const MeshField::Vector3& barycentric_coord, - const Omega_h::Few, 3>& verts_coord) +// Map barycentric coordinates on a simplex (Dim+1 barycentric components) to +// the global Cartesian point, given the simplex's Dim+1 vertex coordinates. +template +[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector GlobalFromBarycentric( + const Omega_h::Vector& barycentric_coord, + const Omega_h::Few, Dim + 1>& verts_coord) { - Omega_h::Vector<2> real_coords = {0.0, 0.0}; - for (int i = 0; i < 3; ++i) { - real_coords[0] += barycentric_coord[i] * verts_coord[i][0]; - real_coords[1] += barycentric_coord[i] * verts_coord[i][1]; + Omega_h::Vector real_coords; + for (int d = 0; d < Dim; ++d) { + real_coords[d] = 0.0; } - return real_coords; -} - -[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<3> EvaluateBarycentric( - const Omega_h::Vector<2>& point, - const r3d::Few, 3>& verts_coord) -{ - Omega_h::Few, 3> omegah_vector; - for (int i = 0; i < 3; ++i) { - omegah_vector[i][0] = verts_coord[i][0]; - omegah_vector[i][1] = verts_coord[i][1]; + for (int i = 0; i < Dim + 1; ++i) { + for (int d = 0; d < Dim; ++d) { + real_coords[d] += barycentric_coord[i] * verts_coord[i][d]; + } } - return Omega_h::barycentric_from_global<2, 2>(point, omegah_vector); + return real_coords; } [[nodiscard]] OMEGA_H_INLINE int RemoveDuplicateVerticesAndFixLinks( @@ -188,22 +185,24 @@ inline void CheckOmegaHScalarLagrangeLayout( return new_n; } -template -OMEGA_H_INLINE void ForEachIntersectionSubtriangle( +// 2D: fan the clipped intersection polygon into triangles anchored at vertex 0, +// invoking op(sub_triangle, src_elm, area) for each non-degenerate piece. +template +OMEGA_H_INLINE void ForEachIntersectionSubtriangleImpl( const int elm, const IntersectionResults& intersection, const Omega_h::Reals& tgt_coords, const Omega_h::Reals& src_coords, - const Omega_h::LOs& tgt_faces2nodes, const Omega_h::LOs& src_faces2nodes, - TriangleOp&& op) + const Omega_h::LOs& tgt_elems2nodes, const Omega_h::LOs& src_elems2nodes, + SimplexOp&& op) { auto tgt_elm_vert_coords = - get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, elm); + get_vert_coords_of_elem<2>(tgt_coords, tgt_elems2nodes, elm); const int start = intersection.tgt2src_offsets[elm]; const int end = intersection.tgt2src_offsets[elm + 1]; for (int i = start; i < end; ++i) { const int current_src_elm = intersection.tgt2src_indices[i]; auto src_elm_vert_coords = - get_vert_coords_of_elem(src_coords, src_faces2nodes, current_src_elm); + get_vert_coords_of_elem<2>(src_coords, src_elems2nodes, current_src_elm); r3d::Polytope<2> poly; r3d::intersect_simplices(poly, tgt_elm_vert_coords, src_elm_vert_coords); auto nverts = RemoveDuplicateVerticesAndFixLinks(poly, 1e-12); @@ -231,40 +230,201 @@ OMEGA_H_INLINE void ForEachIntersectionSubtriangle( continue; } - op(tri_coords, tgt_elm_vert_coords, src_elm_vert_coords, current_src_elm, - area); + op(tri_coords, current_src_elm, area); } } } -// Barycentric integration points and weights for a reference triangle, taken -// from MeshField's predefined triangle quadrature rules and staged on device -// for use in element integration kernels. +// Walk each face of a clipped r3d polyhedron exactly once and fan it into +// triangles, invoking op(v0, v1, v2) for each triangle (v0 is the face's anchor +// vertex, so a face with k vertices yields k-2 triangles). Every vertex of an +// r3d clipped Polytope<3> has exactly three face-neighbors (pnbrs); marking +// each directed edge as it is consumed guarantees every face is emitted once. +// This mirrors the edge-marking traversal buried inside r3d::reduce, which r3d +// does not expose for reuse, so the traversal is reproduced here once and +// shared. Vertices arrive as r3d::Vector<3> (indexable [0..2]). +template +OMEGA_H_INLINE void ForEachPolytopeFaceTriangle(const r3d::Polytope<3>& poly, + TriangleOp&& op) +{ + // emarks[v][p] == 1 once the directed edge (v, pnbr p) has been consumed. + int emarks[r3d::Polytope<3>::max_verts][3] = {{}}; + for (int vstart = 0; vstart < poly.nverts; ++vstart) { + for (int pstart = 0; pstart < 3; ++pstart) { + if (emarks[vstart][pstart]) { + continue; + } + int pnext = pstart; + int vcur = vstart; + emarks[vcur][pnext] = 1; + int vnext = poly.verts[vcur].pnbrs[pnext]; + const auto face_v0 = poly.verts[vcur].pos; + + // Move to the second edge of this face. + int np = 0; + for (np = 0; np < 3; ++np) { + if (poly.verts[vnext].pnbrs[np] == vcur) { + break; + } + } + vcur = vnext; + pnext = (np + 1) % 3; + emarks[vcur][pnext] = 1; + vnext = poly.verts[vcur].pnbrs[pnext]; + + // Fan the face into triangles anchored at face_v0. + while (vnext != vstart) { + op(face_v0, poly.verts[vnext].pos, poly.verts[vcur].pos); + + // Advance around the face. + for (np = 0; np < 3; ++np) { + if (poly.verts[vnext].pnbrs[np] == vcur) { + break; + } + } + vcur = vnext; + pnext = (np + 1) % 3; + emarks[vcur][pnext] = 1; + vnext = poly.verts[vcur].pnbrs[pnext]; + } + } + } +} + +// 3D: star-decompose the clipped intersection polyhedron into tetrahedra from +// its centroid. The centroid lies strictly inside the convex intersection, so +// lifting each boundary-face triangle (enumerated by +// ForEachPolytopeFaceTriangle) to the centroid tiles the polyhedron without +// overlap; op(sub_tet, src_elm, volume) fires for each non-degenerate piece. +template +OMEGA_H_INLINE void ForEachIntersectionSubtetImpl( + const int elm, const IntersectionResults& intersection, + const Omega_h::Reals& tgt_coords, const Omega_h::Reals& src_coords, + const Omega_h::LOs& tgt_elems2nodes, const Omega_h::LOs& src_elems2nodes, + SimplexOp&& op) +{ + auto tgt_elm_vert_coords = + get_vert_coords_of_elem<3>(tgt_coords, tgt_elems2nodes, elm); + const int start = intersection.tgt2src_offsets[elm]; + const int end = intersection.tgt2src_offsets[elm + 1]; + + for (int i = start; i < end; ++i) { + const int current_src_elm = intersection.tgt2src_indices[i]; + auto src_elm_vert_coords = + get_vert_coords_of_elem<3>(src_coords, src_elems2nodes, current_src_elm); + r3d::Polytope<3> poly; + r3d::intersect_simplices(poly, tgt_elm_vert_coords, src_elm_vert_coords); + if (poly.nverts < 4) { + continue; + } + const double poly_vol = Kokkos::fabs(r3d::measure(poly)); + const double eps_vol = abs_tol + rel_tol * poly_vol; + + // Centroid of the (convex) intersection polyhedron: interior apex. + Omega_h::Vector<3> apex = {0.0, 0.0, 0.0}; + for (int v = 0; v < poly.nverts; ++v) { + apex[0] += poly.verts[v].pos[0]; + apex[1] += poly.verts[v].pos[1]; + apex[2] += poly.verts[v].pos[2]; + } + apex[0] /= poly.nverts; + apex[1] /= poly.nverts; + apex[2] /= poly.nverts; + + // Lift each boundary-face triangle to the interior centroid to form a tet. + ForEachPolytopeFaceTriangle(poly, [&](const r3d::Vector<3>& a, + const r3d::Vector<3>& b, + const r3d::Vector<3>& c) { + Omega_h::Few, 4> tet_coords; + tet_coords[0] = apex; + tet_coords[1] = {a[0], a[1], a[2]}; + tet_coords[2] = {b[0], b[1], b[2]}; + tet_coords[3] = {c[0], c[1], c[2]}; + + Omega_h::Few, 3> basis; + basis[0] = tet_coords[1] - tet_coords[0]; + basis[1] = tet_coords[2] - tet_coords[0]; + basis[2] = tet_coords[3] - tet_coords[0]; + + const Omega_h::Real vol = + Kokkos::fabs(Omega_h::tet_volume_from_basis(basis)); + if (vol > eps_vol) { + op(tet_coords, current_src_elm, vol); + } + }); + } +} + +// Dimension-generic driver over the sub-simplices (triangles in 2D, tets in 3D) +// that tile each target element's intersection with the source mesh. Invokes +// op(sub_simplex_coords, src_elm, measure) for every non-degenerate piece. +template +OMEGA_H_INLINE void ForEachIntersectionSubsimplex( + const int elm, const IntersectionResults& intersection, + const Omega_h::Reals& tgt_coords, const Omega_h::Reals& src_coords, + const Omega_h::LOs& tgt_elems2nodes, const Omega_h::LOs& src_elems2nodes, + SimplexOp&& op) +{ + if constexpr (Dim == 3) { + ForEachIntersectionSubtetImpl(elm, intersection, tgt_coords, src_coords, + tgt_elems2nodes, src_elems2nodes, op); + } else { + ForEachIntersectionSubtriangleImpl(elm, intersection, tgt_coords, + src_coords, tgt_elems2nodes, + src_elems2nodes, op); + } +} + +// Maps spatial dimension to the MeshField simplex topology whose reference +// quadrature rules we use (triangle in 2D, tetrahedron in 3D). +template +struct SimplexTopology; +template <> +struct SimplexTopology<2> +{ + static constexpr MeshField::Mesh_Topology value = MeshField::Triangle; +}; +template <> +struct SimplexTopology<3> +{ + static constexpr MeshField::Mesh_Topology value = MeshField::Tetrahedron; +}; + +// Barycentric integration points and weights for a reference simplex (triangle +// in 2D, tetrahedron in 3D), taken from MeshField's predefined quadrature rules +// and staged on device for use in element integration kernels. // // MeshField::getIntegrationPoints returns a host std::vector, which cannot be // dereferenced inside a device kernel, so the (tiny) rule is copied into device -// Kokkos views once at construction. +// Kokkos views once at construction. Each barycentric point has Dim+1 +// components. // // The quadrature order is a runtime argument because the required polynomial // accuracy depends on the source and target element orders (degree = // source_order + target_order), which are only known at construction. +template struct IntegrationData { - Kokkos::View bary_coords; // barycentric coordinates - Kokkos::View weights; // quadrature weights + Kokkos::View + bary_coords; // barycentric coordinates + Kokkos::View weights; // quadrature weights explicit IntegrationData(int order) { - auto ip_vec = MeshField::getIntegrationPoints(order); + auto ip_vec = + MeshField::getIntegrationPoints::value>(order); const std::size_t num_ip = ip_vec.size(); - bary_coords = Kokkos::View("bary_coords", num_ip); + bary_coords = + Kokkos::View("bary_coords", num_ip); weights = Kokkos::View("weights", num_ip); auto bary_coords_host = Kokkos::create_mirror_view(bary_coords); auto weights_host = Kokkos::create_mirror_view(weights); for (std::size_t i = 0; i < num_ip; ++i) { - bary_coords_host(i) = ip_vec[i].param; + for (int d = 0; d < Dim + 1; ++d) { + bary_coords_host(i, d) = ip_vec[i].param[d]; + } weights_host(i) = ip_vec[i].weight; } Kokkos::deep_copy(bary_coords, bary_coords_host); @@ -274,9 +434,10 @@ struct IntegrationData int size() const { return bary_coords.extent(0); } }; -// Target Lagrange basis on a triangle, parameterized by element order, for the +// Target Lagrange basis on a simplex (triangle in 2D, tetrahedron in 3D), +// parameterized by spatial dimension and element order, for the // conservative-projection RHS assembly. Order 0 is a single element-constant -// DOF; order 1 is the three vertex (barycentric) DOFs. Higher orders slot in as +// DOF; order 1 is the Dim+1 vertex (barycentric) DOFs. Higher orders slot in as // additional specializations, kept in lock-step with element_dispatch.h. // // Each specialization provides, for a target element `elm` with local vertex @@ -284,40 +445,39 @@ struct IntegrationData // ndof number of local target DOFs // Index(...) active PETSc row for local dof k // Values(pt, ...) basis values at the (global) integration point pt -template -struct TargetTriBasis; +template +struct TargetSimplexBasis; -template <> -struct TargetTriBasis<0> +template +struct TargetSimplexBasis { static constexpr int ndof = 1; template - KOKKOS_INLINE_FUNCTION static LO Index(const Permutation& permutation, - int elm, - const Omega_h::Few&, - int /*k*/) + KOKKOS_INLINE_FUNCTION static LO Index( + const Permutation& permutation, int elm, + const Omega_h::Few&, int /*k*/) { return permutation(elm); } KOKKOS_INLINE_FUNCTION static void Values( - const Omega_h::Vector<2>&, const Omega_h::Few, 3>&, - Omega_h::Real out[ndof]) + const Omega_h::Vector&, + const Omega_h::Few, Dim + 1>&, Omega_h::Real out[ndof]) { out[0] = 1.0; } }; -template <> -struct TargetTriBasis<1> +template +struct TargetSimplexBasis { - static constexpr int ndof = 3; + static constexpr int ndof = Dim + 1; template KOKKOS_INLINE_FUNCTION static LO Index( const Permutation& permutation, int /*elm*/, - const Omega_h::Few& verts, int k) + const Omega_h::Few& verts, int k) { return permutation(verts[k]); } @@ -325,14 +485,14 @@ struct TargetTriBasis<1> // P1 basis functions are the barycentric coordinates of the target element // evaluated at the (global) integration point. KOKKOS_INLINE_FUNCTION static void Values( - const Omega_h::Vector<2>& pt, - const Omega_h::Few, 3>& tgt_verts, + const Omega_h::Vector& pt, + const Omega_h::Few, Dim + 1>& tgt_verts, Omega_h::Real out[ndof]) { - const auto bary = Omega_h::barycentric_from_global<2, 2>(pt, tgt_verts); - out[0] = bary[0]; - out[1] = bary[1]; - out[2] = bary[2]; + const auto bary = Omega_h::barycentric_from_global(pt, tgt_verts); + for (int i = 0; i < ndof; ++i) { + out[i] = bary[i]; + } } }; diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp index ae9b293f..f3c16df1 100644 --- a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp @@ -25,7 +25,7 @@ struct Data PetscInt num_target_dofs = 0; }; -template +template Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, const OmegaHLagrangeLayout& target_layout, int quad_order); @@ -41,24 +41,39 @@ Data BuildData(const std::shared_ptr& source_layout, target_coordinate_system, target_layout, "OmegaHIntersectionRHSIntegrator", "target"); + const int dim = target_layout->GetMesh().dim(); + if (source_layout->GetMesh().dim() != dim) { + throw pcms_error("OmegaHIntersectionRHSIntegrator: source and target mesh " + "dimensions differ"); + } + // The integrand f_src * phi_target has polynomial degree source_order + - // target_order on each intersection subtriangle; integrate it exactly (with a + // target_order on each intersection sub-simplex; integrate it exactly (with a // 1-point floor so a P0->P0 pair still gets a valid rule). const int quad_order = std::max(1, source_layout->GetOrder() + target_layout->GetOrder()); return detail::DispatchByOrder(target_layout->GetOrder(), [&](auto order_c) { constexpr int TgtOrder = decltype(order_c)::value; - return BuildDataImpl(*source_layout, *target_layout, quad_order); + if (dim == 3) { + return BuildDataImpl<3, TgtOrder>(*source_layout, *target_layout, + quad_order); + } + return BuildDataImpl<2, TgtOrder>(*source_layout, *target_layout, + quad_order); }); } -template +template Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, const OmegaHLagrangeLayout& target_layout, int quad_order) { - using Basis = detail::TargetTriBasis; + using Basis = detail::TargetSimplexBasis; constexpr int ndof = Basis::ndof; + // Reference-to-physical Jacobian factor for a simplex: the reference simplex + // measure is 1/Dim! (1/2 in 2D, 1/6 in 3D), so a physical sub-simplex of + // measure `m` scales the reference quadrature weights by Dim! * m. + constexpr Omega_h::Real ref_factor = (Dim == 3) ? 6.0 : 2.0; Omega_h::Mesh& source_mesh = source_layout.GetMesh(); Omega_h::Mesh& target_mesh = target_layout.GetMesh(); @@ -66,13 +81,11 @@ Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, const auto intersections = intersectTargets(source_mesh, target_mesh); const auto& tgt_coords = target_mesh.coords(); - const auto& tgt_faces2nodes = - target_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& tgt_elems2nodes = target_mesh.ask_down(Dim, Omega_h::VERT).ab2b; const auto& src_coords = source_mesh.coords(); - const auto& src_faces2nodes = - source_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& src_elems2nodes = source_mesh.ask_down(Dim, Omega_h::VERT).ab2b; - detail::IntegrationData ip_data(quad_order); + detail::IntegrationData ip_data(quad_order); const int npts = ip_data.size(); auto bary_coords = ip_data.bary_coords; // device view auto weights = ip_data.weights; // device view @@ -88,12 +101,10 @@ Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, Kokkos::parallel_for( "rhs_count", nelems, KOKKOS_LAMBDA(int elm) { int count = 0; - detail::ForEachIntersectionSubtriangle( + detail::ForEachIntersectionSubsimplex( elm, {tgt2src_offsets, tgt2src_indices}, tgt_coords, src_coords, - tgt_faces2nodes, src_faces2nodes, - [&](const Omega_h::Few, 3>&, - const r3d::Few, 3>&, - const r3d::Few, 3>&, int, + tgt_elems2nodes, src_elems2nodes, + [&](const Omega_h::Few, Dim + 1>&, int, Omega_h::Real) { count += npts; }); ip_counts[elm] = count; }); @@ -105,7 +116,7 @@ Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, // Pass 2: fill coords, node_gids, and coeffs on device. node_gids/coeffs hold // ndof (target DOFs per element) entries per integration point. - Kokkos::View coords("rhs_coords", num_pts, 2); + Kokkos::View coords("rhs_coords", num_pts, Dim); Kokkos::View node_gids( "rhs_node_gids", static_cast(num_pts) * ndof); Kokkos::View coeffs( @@ -113,37 +124,45 @@ Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, Kokkos::parallel_for( "rhs_fill", nelems, KOKKOS_LAMBDA(int elm) { - const auto tgt_verts = Omega_h::gather_verts<3>(tgt_faces2nodes, elm); - const Omega_h::Matrix<2, 3> tgt_vert_mat = - Omega_h::gather_vectors<3, 2>(tgt_coords, tgt_verts); - Omega_h::Few, 3> tgt_omh; - for (int i = 0; i < 3; ++i) - tgt_omh[i] = {tgt_vert_mat[i][0], tgt_vert_mat[i][1]}; + const auto tgt_verts = + Omega_h::gather_verts(tgt_elems2nodes, elm); + const Omega_h::Matrix tgt_vert_mat = + Omega_h::gather_vectors(tgt_coords, tgt_verts); + Omega_h::Few, Dim + 1> tgt_omh; + for (int i = 0; i < Dim + 1; ++i) { + for (int d = 0; d < Dim; ++d) { + tgt_omh[i][d] = tgt_vert_mat[i][d]; + } + } int ip_local = 0; const int offset = ip_offsets[elm]; - detail::ForEachIntersectionSubtriangle( + detail::ForEachIntersectionSubsimplex( elm, {tgt2src_offsets, tgt2src_indices}, tgt_coords, src_coords, - tgt_faces2nodes, src_faces2nodes, - [&](const Omega_h::Few, 3>& tri, - const r3d::Few, 3>&, - const r3d::Few, 3>&, int, Omega_h::Real area) { + tgt_elems2nodes, src_elems2nodes, + [&](const Omega_h::Few, Dim + 1>& sub, int, + Omega_h::Real measure) { for (int ip_idx = 0; ip_idx < npts; ++ip_idx) { - const auto bary = bary_coords(ip_idx); + Omega_h::Vector bary; + for (int d = 0; d < Dim + 1; ++d) { + bary[d] = bary_coords(ip_idx, d); + } const double w = weights(ip_idx); - const auto pt = detail::GlobalFromBarycentric(bary, tri); + const auto pt = detail::GlobalFromBarycentric(bary, sub); Omega_h::Real basis[ndof]; Basis::Values(pt, tgt_omh, basis); const int global_ip = offset + ip_local; - coords(global_ip, 0) = pt[0]; - coords(global_ip, 1) = pt[1]; + for (int d = 0; d < Dim; ++d) { + coords(global_ip, d) = pt[d]; + } for (int k = 0; k < ndof; ++k) { node_gids(global_ip * ndof + k) = static_cast( Basis::Index(global_to_local, elm, tgt_verts, k)); - coeffs(global_ip * ndof + k) = basis[k] * w * 2.0 * area; + coeffs(global_ip * ndof + k) = + basis[k] * w * ref_factor * measure; } ++ip_local; } diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp index 29cd90e2..99edb7eb 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.cpp +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -17,10 +17,24 @@ namespace pcms namespace { +// Measure (area in 2D, volume in 3D) of a simplex from its vertex-difference +// basis. +template +KOKKOS_INLINE_FUNCTION Omega_h::Real SimplexMeasure( + const Omega_h::Few, Dim>& basis) +{ + if constexpr (Dim == 3) { + return Kokkos::fabs(Omega_h::tet_volume_from_basis(basis)); + } else { + return Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + } +} + // Fills the diagonal COO entries of a P0 (piecewise-constant) mass matrix: one -// entry per element on its own DOF, valued at the element area. +// entry per element on its own DOF, valued at the element measure. +template void FillP0MassCoo( - int nelems, const Omega_h::Reals& coords, const Omega_h::LOs& faces2nodes, + int nelems, const Omega_h::Reals& coords, const Omega_h::LOs& elems2nodes, const Kokkos::View& global_to_local, const Kokkos::View& coo_rows, const Kokkos::View& coo_cols, @@ -28,35 +42,37 @@ void FillP0MassCoo( { Kokkos::parallel_for( "mass_p0_diag", nelems, KOKKOS_LAMBDA(int e) { - const auto verts = Omega_h::gather_verts<3>(faces2nodes, e); - const Omega_h::Matrix<2, 3> vm = - Omega_h::gather_vectors<3, 2>(coords, verts); - Omega_h::Few, 2> basis; - basis[0] = vm[1] - vm[0]; - basis[1] = vm[2] - vm[0]; - const Omega_h::Real area = - Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + const auto verts = Omega_h::gather_verts(elems2nodes, e); + const Omega_h::Matrix vm = + Omega_h::gather_vectors(coords, verts); + Omega_h::Few, Dim> basis; + for (int d = 0; d < Dim; ++d) { + basis[d] = vm[d + 1] - vm[0]; + } + const Omega_h::Real measure = SimplexMeasure(basis); const PetscInt g = static_cast(global_to_local(e)); coo_rows(e) = g; coo_cols(e) = g; - vals(e) = static_cast(area); + vals(e) = static_cast(measure); }); } -// Fills the 3x3-block COO sparsity pattern of a P1 (linear) mass matrix: each -// element contributes a dense block coupling its three vertex DOFs. +// Fills the (Dim+1)x(Dim+1)-block COO sparsity pattern of a P1 (linear) mass +// matrix: each element contributes a dense block coupling its vertex DOFs. +template void FillP1MassCooPattern( - int nelems, const Omega_h::LOs& faces2nodes, + int nelems, const Omega_h::LOs& elems2nodes, const Kokkos::View& global_to_local, const Kokkos::View& coo_rows, const Kokkos::View& coo_cols) { + constexpr int nv = Dim + 1; Kokkos::parallel_for( "mass_coo_pattern", nelems, KOKKOS_LAMBDA(int e) { - const auto verts = Omega_h::gather_verts<3>(faces2nodes, e); - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { - const int idx = e * 9 + i * 3 + j; + const auto verts = Omega_h::gather_verts(elems2nodes, e); + for (int i = 0; i < nv; ++i) { + for (int j = 0; j < nv; ++j) { + const int idx = e * (nv * nv) + i * nv + j; coo_rows(idx) = static_cast(global_to_local(verts[i])); coo_cols(idx) = static_cast(global_to_local(verts[j])); } @@ -64,92 +80,115 @@ void FillP1MassCooPattern( }); } -} // namespace - -OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space) - : OmegaHMassIntegrator(std::dynamic_pointer_cast( - target_space.GetLayout()), - target_space.GetCoordinateSystem()) +// Assembles the target-space mass matrix for spatial dimension Dim (triangles +// for Dim==2, tetrahedra for Dim==3) and returns the owned PETSc matrix. +template +Mat BuildOmegaHMassMatrixImpl(Omega_h::Mesh& mesh, + const OmegaHLagrangeLayout& target_layout) { -} - -OmegaHMassIntegrator::OmegaHMassIntegrator( - std::shared_ptr target_layout, - CoordinateSystem coordinate_system) -{ - detail::CheckOmegaHScalarLagrangeLayout(coordinate_system, target_layout, - "OmegaHMassIntegrator", "target"); - - Omega_h::Mesh& mesh = target_layout->GetMesh(); - const auto global_to_local = target_layout->GetGlobalToLocalPermutation(); + const auto global_to_local = target_layout.GetGlobalToLocalPermutation(); const PetscInt num_dofs = - static_cast(target_layout->GetNumOwnedDofHolder()); + static_cast(target_layout.GetNumOwnedDofHolder()); const int nelems = mesh.nelems(); + Mat mat = nullptr; - if (target_layout->GetOrder() == 0) { + if (target_layout.GetOrder() == 0) { // P0 target: piecewise-constant basis functions have disjoint support, so - // the mass matrix is diagonal with M_ee = area(e). One COO entry per + // the mass matrix is diagonal with M_ee = measure(e). One COO entry per // element on its own (element-id) diagonal. const auto& coords = mesh.coords(); - const auto& faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& elems2nodes = mesh.ask_down(Dim, Omega_h::VERT).ab2b; const PetscInt nnz = static_cast(nelems); Kokkos::View coo_rows("mass_coo_rows", nnz); Kokkos::View coo_cols("mass_coo_cols", nnz); Kokkos::View vals("mass_vals", nnz); - FillP0MassCoo(nelems, coords, faces2nodes, global_to_local, coo_rows, - coo_cols, vals); + FillP0MassCoo(nelems, coords, elems2nodes, global_to_local, coo_rows, + coo_cols, vals); PetscErrorCode ierr = - createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); + createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat); CHKERRABORT(PETSC_COMM_WORLD, ierr); auto coo_rows_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_rows); auto coo_cols_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_cols); - ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows_host.data(), + ierr = MatSetPreallocationCOO(mat, nnz, coo_rows_host.data(), coo_cols_host.data()); CHKERRABORT(PETSC_COMM_WORLD, ierr); - ierr = MatSetValuesCOO(mat_, vals.data(), INSERT_VALUES); + ierr = MatSetValuesCOO(mat, vals.data(), INSERT_VALUES); CHKERRABORT(PETSC_COMM_WORLD, ierr); - return; + return mat; } - // P1 target: consistent mass matrix assembled from MeshField per-element 3x3 - // blocks. (Higher MeshField orders extend this branch via - // getTriangleElement.) - MeshField::OmegahMeshField omf(mesh); auto coordField = omf.getCoordField(); - const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); - MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); - auto elm_mass_dev = buildElementMassMatrix(mesh, coordFe); + Kokkos::View elm_mass_dev; + if constexpr (Dim == 3) { + const auto [shp, map] = MeshField::Omegah::getTetrahedronElement<1>(mesh); + MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); + elm_mass_dev = buildElementMassMatrix(mesh, coordFe); + } else { + const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); + MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); + elm_mass_dev = buildElementMassMatrix(mesh, coordFe); + } - // Build COO sparsity pattern on device: each element contributes a 3x3 block. - const auto& faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + // Build COO sparsity pattern on device: each element contributes a + // (Dim+1)x(Dim+1) block. + const auto& elems2nodes = mesh.ask_down(Dim, Omega_h::VERT).ab2b; - const PetscInt nnz = static_cast(nelems) * 9; + constexpr int nv = Dim + 1; + const PetscInt nnz = static_cast(nelems) * (nv * nv); Kokkos::View coo_rows("mass_coo_rows", nnz); Kokkos::View coo_cols("mass_coo_cols", nnz); - FillP1MassCooPattern(nelems, faces2nodes, global_to_local, coo_rows, - coo_cols); + FillP1MassCooPattern(nelems, elems2nodes, global_to_local, coo_rows, + coo_cols); // Create sparse matrix, preallocate with COO pattern, then bulk-set values. // elm_mass_dev is in the same element-major order as coo_rows/coo_cols, so // it can be passed directly to MatSetValuesCOO — no host copy needed. PetscErrorCode ierr = - createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); + createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat); CHKERRABORT(PETSC_COMM_WORLD, ierr); auto coo_rows_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_rows); auto coo_cols_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_cols); - ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows_host.data(), + ierr = MatSetPreallocationCOO(mat, nnz, coo_rows_host.data(), coo_cols_host.data()); CHKERRABORT(PETSC_COMM_WORLD, ierr); - ierr = MatSetValuesCOO(mat_, elm_mass_dev.data(), INSERT_VALUES); + ierr = MatSetValuesCOO(mat, elm_mass_dev.data(), INSERT_VALUES); CHKERRABORT(PETSC_COMM_WORLD, ierr); + return mat; +} + +} // namespace + +OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space) + : OmegaHMassIntegrator(std::dynamic_pointer_cast( + target_space.GetLayout()), + target_space.GetCoordinateSystem()) +{ +} + +OmegaHMassIntegrator::OmegaHMassIntegrator( + std::shared_ptr target_layout, + CoordinateSystem coordinate_system) +{ + detail::CheckOmegaHScalarLagrangeLayout(coordinate_system, target_layout, + "OmegaHMassIntegrator", "target"); + + Omega_h::Mesh& mesh = target_layout->GetMesh(); + if (mesh.dim() == 3) { + mat_ = BuildOmegaHMassMatrixImpl<3>(mesh, *target_layout); + } else { + mat_ = BuildOmegaHMassMatrixImpl<2>(mesh, *target_layout); + } } OmegaHMassIntegrator::~OmegaHMassIntegrator() diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp index c57ce9f1..05016379 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp @@ -12,9 +12,9 @@ namespace pcms namespace { -// Maps a uniform point on the unit square to uniform barycentric coordinates -// Shape Distributions (ACM Transactions on Graphics, Vol. 21, No. 4, October -// 2002.) page 814 Eq 1 +// Maps a uniform point on the unit square to uniform barycentric coordinates on +// a triangle. Shape Distributions (ACM Transactions on Graphics, Vol. 21, +// No. 4, October 2002.) page 814 Eq 1. KOKKOS_INLINE_FUNCTION Omega_h::Vector<3> UniformTriangleBarycentric(Real u, Real v) { @@ -22,51 +22,96 @@ KOKKOS_INLINE_FUNCTION Omega_h::Vector<3> UniformTriangleBarycentric(Real u, return {1.0 - s, s * (1.0 - v), s * v}; } +// Maps three uniform draws on the unit cube to uniform barycentric coordinates +// on a tetrahedron via the cut-and-fold method (Rocchini & Cignoni, +// "Generating Random Points in a Tetrahedron", J. Graphics Tools 2000). +KOKKOS_INLINE_FUNCTION Omega_h::Vector<4> UniformTetBarycentric(Real s, Real t, + Real u) +{ + if (s + t > 1.0) { // fold the cube into a prism + s = 1.0 - s; + t = 1.0 - t; + } + if (t + u > 1.0) { // fold the prism into a tetrahedron + const Real tmp = u; + u = 1.0 - s - t; + t = 1.0 - tmp; + } else if (s + t + u > 1.0) { + const Real tmp = u; + u = s + t + u - 1.0; + s = 1.0 - t - tmp; + } + const Real a = 1.0 - s - t - u; + return {a, s, t, u}; +} + +template +KOKKOS_INLINE_FUNCTION Omega_h::Vector UniformSimplexBarycentric( + const Real r[Dim]) +{ + if constexpr (Dim == 3) { + return UniformTetBarycentric(r[0], r[1], r[2]); + } else { + return UniformTriangleBarycentric(r[0], r[1]); + } +} + // Fills coords, node_gids, and coeffs for all samples of one target element. -// unit_sample_at(s, u, v) provides the s-th unit-square sample. -template +// unit_sample_at(s, r) fills the s-th sample's Dim unit-hypercube draws. +template KOKKOS_INLINE_FUNCTION void FillElementSamples( const int elm, const int samples_per_element, - const Omega_h::Reals& mesh_coords, const Omega_h::LOs& faces2nodes, + const Omega_h::Reals& mesh_coords, const Omega_h::LOs& elems2nodes, const Kokkos::View& global_to_local, const Kokkos::View& coords, const Kokkos::View& node_gids, const Kokkos::View& coeffs, const UnitSampleAt& unit_sample_at) { - const auto verts = Omega_h::gather_verts<3>(faces2nodes, elm); - const auto vert_coords = Omega_h::gather_vectors<3, 2>(mesh_coords, verts); - Omega_h::Few, 2> basis; - basis[0] = vert_coords[1] - vert_coords[0]; - basis[1] = vert_coords[2] - vert_coords[0]; - const Real area = Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); - const Real weight = area / samples_per_element; + constexpr int nv = Dim + 1; + const auto verts = Omega_h::gather_verts(elems2nodes, elm); + const auto vert_coords = Omega_h::gather_vectors(mesh_coords, verts); + Omega_h::Few, Dim> basis; + for (int d = 0; d < Dim; ++d) { + basis[d] = vert_coords[d + 1] - vert_coords[0]; + } + Real measure; + if constexpr (Dim == 3) { + measure = Kokkos::fabs(Omega_h::tet_volume_from_basis(basis)); + } else { + measure = Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + } + const Real weight = measure / samples_per_element; for (int s = 0; s < samples_per_element; ++s) { - Real u = 0.0; - Real v = 0.0; - unit_sample_at(s, u, v); - const auto bary = UniformTriangleBarycentric(u, v); + Real r[Dim]; + unit_sample_at(s, r); + const auto bary = UniformSimplexBarycentric(r); const int i = elm * samples_per_element + s; - Real x = 0.0; - Real y = 0.0; - for (int k = 0; k < 3; ++k) { - x += bary[k] * vert_coords[k][0]; - y += bary[k] * vert_coords[k][1]; - node_gids(i * 3 + k) = static_cast(global_to_local(verts[k])); - coeffs(i * 3 + k) = bary[k] * weight; + Omega_h::Vector x; + for (int d = 0; d < Dim; ++d) { + x[d] = 0.0; + } + for (int k = 0; k < nv; ++k) { + for (int d = 0; d < Dim; ++d) { + x[d] += bary[k] * vert_coords[k][d]; + } + node_gids(i * nv + k) = static_cast(global_to_local(verts[k])); + coeffs(i * nv + k) = bary[k] * weight; + } + for (int d = 0; d < Dim; ++d) { + coords(i, d) = x[d]; } - coords(i, 0) = x; - coords(i, 1) = y; } } -// Samples samples_per_element from a uniform random distribution over the +// Samples samples_per_element from a uniform random distribution over each // target element, writing their coordinates, node GIDs, and coefficients. +template void FillElementSamples( int nelems, int samples_per_element, const Omega_h::Reals& mesh_coords, - const Omega_h::LOs& faces2nodes, + const Omega_h::LOs& elems2nodes, const Kokkos::View& global_to_local, const Kokkos::View& coords, const Kokkos::View& node_gids, @@ -77,12 +122,13 @@ void FillElementSamples( "mc_rhs_fill_random", Kokkos::RangePolicy(0, nelems), KOKKOS_LAMBDA(int elm) { auto gen = pool.get_state(); - FillElementSamples(elm, samples_per_element, mesh_coords, faces2nodes, - global_to_local, coords, node_gids, coeffs, - [&](int /*s*/, Real& u, Real& v) { - u = gen.drand(); - v = gen.drand(); - }); + FillElementSamples(elm, samples_per_element, mesh_coords, + elems2nodes, global_to_local, coords, node_gids, + coeffs, [&](int /*s*/, Real r[Dim]) { + for (int d = 0; d < Dim; ++d) { + r[d] = gen.drand(); + } + }); pool.free_state(gen); }); Kokkos::fence(); @@ -113,22 +159,29 @@ OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( } Omega_h::Mesh& mesh = target_layout->GetMesh(); + const int dim = mesh.dim(); + nbary_ = dim + 1; const int nelems = mesh.nelems(); const int num_samples = nelems * samples_per_element; const auto mesh_coords = mesh.coords(); - const auto faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto elems2nodes = mesh.ask_down(dim, Omega_h::VERT).ab2b; const auto global_to_local = target_layout->GetGlobalToLocalPermutation(); Kokkos::View coords("mc_rhs_coords", num_samples, - 2); + dim); Kokkos::View node_gids( - "mc_rhs_node_gids", static_cast(num_samples) * 3); + "mc_rhs_node_gids", static_cast(num_samples) * nbary_); Kokkos::View coeffs( - "mc_rhs_coeffs", static_cast(num_samples) * 3); + "mc_rhs_coeffs", static_cast(num_samples) * nbary_); - FillElementSamples(nelems, samples_per_element, mesh_coords, faces2nodes, - global_to_local, coords, node_gids, coeffs, seed); + if (dim == 3) { + FillElementSamples<3>(nelems, samples_per_element, mesh_coords, elems2nodes, + global_to_local, coords, node_gids, coeffs, seed); + } else { + FillElementSamples<2>(nelems, samples_per_element, mesh_coords, elems2nodes, + global_to_local, coords, node_gids, coeffs, seed); + } coords_ = std::move(coords); node_gids_ = std::move(node_gids); @@ -173,8 +226,9 @@ Vec OmegaHMonteCarloRHSIntegrator::GetVector() const noexcept void OmegaHMonteCarloRHSIntegrator::Assemble( Rank2View sampled_values) { + const int nbary = nbary_; const std::size_t num_samples = - static_cast(node_gids_.extent(0) / 3); + static_cast(node_gids_.extent(0) / nbary); PCMS_ALWAYS_ASSERT(static_cast(sampled_values.extent(0)) == num_samples); PCMS_ALWAYS_ASSERT(sampled_values.extent(1) >= 1); @@ -187,13 +241,14 @@ void OmegaHMonteCarloRHSIntegrator::Assemble( sampled_values.data_handle(), sampled_values.extent(0), sampled_values.extent(1)); Kokkos::View coo_vals("mc_rhs_coo_vals", - num_samples * 3); + num_samples * nbary); auto coeffs = coeffs_; Kokkos::parallel_for( "mc_rhs_coo_vals", static_cast(num_samples), KOKKOS_LAMBDA(int i) { const PetscScalar f = static_cast(sv(i, 0)); - for (int j = 0; j < 3; ++j) { - coo_vals(i * 3 + j) = static_cast(coeffs(i * 3 + j)) * f; + for (int j = 0; j < nbary; ++j) { + coo_vals(i * nbary + j) = + static_cast(coeffs(i * nbary + j)) * f; } }); diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp index c9b533b4..5d8946eb 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.hpp @@ -14,7 +14,8 @@ namespace pcms { // Monte Carlo RHS integrator for conservative L2 projection onto order-1 -// Lagrange spaces on Omega_h 2D simplex meshes. +// Lagrange spaces on Omega_h simplex meshes (triangles in 2D, tetrahedra in +// 3D). // // Instead of intersection-based quadrature, each load-vector entry // b_j = \int phi_j f dx @@ -57,6 +58,7 @@ class OmegaHMonteCarloRHSIntegrator : public LinearFormIntegrator coords_; // [num_pts][dim] sample coordinates Kokkos::View node_gids_; // COO indices Kokkos::View coeffs_; // basis * |T| / N + int nbary_ = 3; // barycentric DOFs per sample (Dim+1): 3 in 2D, 4 in 3D }; } // namespace pcms From c1dbc2a4ea1ecb4cb6071f756ff6654793e4593b Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Thu, 6 Aug 2026 09:37:48 -0400 Subject: [PATCH 4/4] test: add 3D intersection and conservative projection tests Cover the tetrahedral paths added to the transfer machinery: - test_omega_h_3d_conservative_projection.cpp exercises conservative L2 projection between 3D tet meshes and checks that the integral of the field is preserved - test_omega_h_form_integrator_utils.cpp unit tests the reference mapping and barycentric helpers in 2D and 3D - test_intersections.cpp is updated for the dimension-templated intersection search Both new tests use the shared helpers in field_test_utils.h. --- test/CMakeLists.txt | 3 +- test/test_intersections.cpp | 4 +- ...est_omega_h_3d_conservative_projection.cpp | 227 ++++++++++++++++++ test/test_omega_h_form_integrator_utils.cpp | 118 +++++++++ 4 files changed, 349 insertions(+), 3 deletions(-) create mode 100644 test/test_omega_h_3d_conservative_projection.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5efdc81e..06a24e7a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -413,7 +413,8 @@ if(Catch2_FOUND) test_omega_h_intersection_rhs_integrator.cpp test_omega_h_mass_integrator.cpp test_omega_h_mc_rhs_integrator.cpp - test_mesh_intersection_field_transfer.cpp) + test_mesh_intersection_field_transfer.cpp + test_omega_h_3d_conservative_projection.cpp) endif() add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES}) diff --git a/test/test_intersections.cpp b/test/test_intersections.cpp index 02598f0a..36ccb766 100644 --- a/test/test_intersections.cpp +++ b/test/test_intersections.cpp @@ -111,7 +111,7 @@ TEST_CASE("Mesh intersection test with source and target", "[intersection]") Omega_h::parallel_for( ntgt, OMEGA_H_LAMBDA(int t) { auto tgt_vert_coords = - get_vert_coords_of_elem(tgt_coords, tgt_faces2verts, t); + get_vert_coords_of_elem<2>(tgt_coords, tgt_faces2verts, t); int start = intersection.tgt2src_offsets[t]; int end = intersection.tgt2src_offsets[t + 1]; @@ -120,7 +120,7 @@ TEST_CASE("Mesh intersection test with source and target", "[intersection]") for (int i = start; i < end; ++i) { int sid = intersection.tgt2src_indices[i]; auto src_vert_coords = - get_vert_coords_of_elem(src_coords, src_faces2verts, sid); + get_vert_coords_of_elem<2>(src_coords, src_faces2verts, sid); r3d::Polytope<2> poly; r3d::intersect_simplices(poly, tgt_vert_coords, src_vert_coords); diff --git a/test/test_omega_h_3d_conservative_projection.cpp b/test/test_omega_h_3d_conservative_projection.cpp new file mode 100644 index 00000000..d2783ae3 --- /dev/null +++ b/test/test_omega_h_3d_conservative_projection.cpp @@ -0,0 +1,227 @@ +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include "field_test_utils.h" + +#include + +namespace +{ + +// Build a tetrahedral box mesh of [0,1]^3 at the given resolution. build_box +// produces a classified simplex mesh, which is what the Omega_h Lagrange +// function space requires. +Omega_h::Mesh BuildUnitCube(Omega_h::Library& lib, int n) +{ + return Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, 1.0, n, n, + n); +} + +} // namespace + +TEST_CASE("OmegaHConservativeProjection (3D tets) reproduces constant and " + "linear fields", + "[transfer][mesh_intersection][3d]") +{ + Omega_h::Library lib; + + // Two independent tessellations of the same unit cube. + Omega_h::Mesh source_mesh = BuildUnitCube(lib, 1); + Omega_h::Mesh target_mesh = BuildUnitCube(lib, 2); + + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); + + auto source = source_space->CreateFunction(); + auto target = target_space->CreateFunction(); + + pcms::OmegaHConservativeProjection projection(*source_space, *target_space); + + SECTION("constant field is preserved and conserved") + { + const double c = 2.0; + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real, pcms::Real, pcms::Real) { return c; }); + + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + REQUIRE(static_cast(target_values.size()) == + target_mesh.nverts()); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-9)); + } + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-9)); + } + + SECTION("linear field is reproduced on target vertices and conserved") + { + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 3), target_mesh.nverts(), + 3); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = 1.0 + tgt_coords_h(i, 0) + + 2.0 * tgt_coords_h(i, 1) + + 3.0 * tgt_coords_h(i, 2); + REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-8)); + } + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-8)); + } +} + +TEST_CASE("OmegaHConservativeProjection (3D tets) conserves the integral for a " + "P0 target", + "[transfer][mesh_intersection][3d]") +{ + Omega_h::Library lib; + + Omega_h::Mesh source_mesh = BuildUnitCube(lib, 2); + Omega_h::Mesh target_mesh = BuildUnitCube(lib, 1); + + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP0Space(target_mesh); + + auto source = source_space->CreateFunction(); + auto target = target_space->CreateFunction(); + + pcms::OmegaHConservativeProjection projection(*source_space, *target_space); + + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + projection.Apply(source, target); + + REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-8)); +} + +TEST_CASE("Copy transfer (3D tets) reproduces the source field", + "[transfer][copy][3d]") +{ + Omega_h::Library lib; + Omega_h::Mesh mesh = BuildUnitCube(lib, 2); + auto space = pcms::test::MakeP1Space(mesh); + + auto source = space->CreateFunction(); + auto target = space->CreateFunction(); + + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + pcms::Copy copy(*space, *space); + copy.Apply(source, target); + + const auto sv = pcms::FlattenToRank1View(source.GetDOFHolderDataHost()); + const auto tv = pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + REQUIRE(tv.size() == sv.size()); + for (std::size_t i = 0; i < tv.size(); ++i) { + REQUIRE(tv[i] == Catch::Approx(sv[i])); + } +} + +// Point interpolation evaluates the source field at each target DOF site. A +// linear field is reproduced exactly at the target vertices in 3D. +TEST_CASE("Interpolation transfer (3D tets) reproduces a linear field", + "[transfer][interpolation][3d]") +{ + Omega_h::Library lib; + Omega_h::Mesh source_mesh = BuildUnitCube(lib, 2); + Omega_h::Mesh target_mesh = BuildUnitCube(lib, 3); + + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); + + auto source = source_space->CreateFunction(); + auto target = target_space->CreateFunction(); + + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + pcms::Interpolator interp(*source_space, *target_space); + interp.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 3), target_mesh.nverts(), 3); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = 1.0 + tgt_coords_h(i, 0) + + 2.0 * tgt_coords_h(i, 1) + 3.0 * tgt_coords_h(i, 2); + REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-8)); + } +} + +// The Monte-Carlo/control-variate projection uses the source field interpolated +// onto the target space as a control variate, so a field already representable +// in the target P1 space (an affine function) is reproduced exactly and its +// integral conserved, even with very few stochastic samples. +TEST_CASE("OmegaHControlVariateProjection (3D tets) is exact for target-space " + "fields", + "[transfer][monte_carlo][3d]") +{ + Omega_h::Library lib; + + Omega_h::Mesh source_mesh = BuildUnitCube(lib, 1); + Omega_h::Mesh target_mesh = BuildUnitCube(lib, 2); + + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); + + auto source = source_space->CreateFunction(); + auto target = target_space->CreateFunction(); + + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return 1.0 + x + 2.0 * y + 3.0 * z; + }); + + pcms::OmegaHControlVariateProjection projection( + *source_space, *target_space, /*samples_per_element=*/8, + pcms::MonteCarloSampling::UniformRandom, /*seed=*/12345); + projection.Apply(source, target); + + const auto target_values = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 3), target_mesh.nverts(), 3); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = 1.0 + tgt_coords_h(i, 0) + + 2.0 * tgt_coords_h(i, 1) + 3.0 * tgt_coords_h(i, 2); + REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-8)); + } + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-8)); +} diff --git a/test/test_omega_h_form_integrator_utils.cpp b/test/test_omega_h_form_integrator_utils.cpp index 91cedac7..785fd19a 100644 --- a/test/test_omega_h_form_integrator_utils.cpp +++ b/test/test_omega_h_form_integrator_utils.cpp @@ -2,6 +2,10 @@ #include #include +#include +#include +#include + namespace { @@ -78,3 +82,117 @@ TEST_CASE( REQUIRE(n == 3); REQUIRE(poly.nverts == 3); } + +// --------------------------------------------------------------------------- +// ForEachPolytopeFaceTriangle: the 3D face-walk that ForEachIntersectionSubtet +// relies on to star-decompose a clipped intersection polyhedron. +// --------------------------------------------------------------------------- + +namespace +{ + +double TetVolume(const r3d::Vector<3>& a, const r3d::Vector<3>& b, + const r3d::Vector<3>& c, const r3d::Vector<3>& d) +{ + const double bx = b[0] - a[0], by = b[1] - a[1], bz = b[2] - a[2]; + const double cx = c[0] - a[0], cy = c[1] - a[1], cz = c[2] - a[2]; + const double dx = d[0] - a[0], dy = d[1] - a[1], dz = d[2] - a[2]; + const double triple = bx * (cy * dz - cz * dy) - by * (cx * dz - cz * dx) + + bz * (cx * dy - cy * dx); + return std::abs(triple) / 6.0; +} + +// Star-decompose `poly` from its centroid exactly as ForEachIntersectionSubtet +// does, summing the sub-tet volumes and counting the emitted face triangles. +std::pair DecomposeFromCentroid(const r3d::Polytope<3>& poly) +{ + r3d::Vector<3> apex; + apex[0] = apex[1] = apex[2] = 0.0; + for (int v = 0; v < poly.nverts; ++v) { + apex[0] += poly.verts[v].pos[0]; + apex[1] += poly.verts[v].pos[1]; + apex[2] += poly.verts[v].pos[2]; + } + apex[0] /= poly.nverts; + apex[1] /= poly.nverts; + apex[2] /= poly.nverts; + + double total_vol = 0.0; + int ntri = 0; + pcms::detail::ForEachPolytopeFaceTriangle(poly, [&](const r3d::Vector<3>& a, + const r3d::Vector<3>& b, + const r3d::Vector<3>& c) { + total_vol += TetVolume(apex, a, b, c); + ++ntri; + }); + return {total_vol, ntri}; +} + +r3d::Few, 4> MakeTetVerts( + std::initializer_list> pts) +{ + r3d::Few, 4> verts; + int i = 0; + for (const auto& p : pts) { + verts[i][0] = p[0]; + verts[i][1] = p[1]; + verts[i][2] = p[2]; + ++i; + } + return verts; +} + +} // namespace + +TEST_CASE("ForEachPolytopeFaceTriangle: a tetrahedron yields 4 face triangles " + "that tile its volume", + "[form_integrator_utils]") +{ + // A tetrahedron has 4 triangular faces, so the walk must emit exactly 4 + // triangles (= 2*(nverts-2) for nverts==4), and the centroid star-tiling must + // reproduce the tet's volume (1/6 for the reference tet). + r3d::Polytope<3> poly; + r3d::init(poly, MakeTetVerts({{{0.0, 0.0, 0.0}}, + {{1.0, 0.0, 0.0}}, + {{0.0, 1.0, 0.0}}, + {{0.0, 0.0, 1.0}}})); + + const auto [vol, ntri] = DecomposeFromCentroid(poly); + REQUIRE(ntri == 4); + REQUIRE(ntri == 2 * (poly.nverts - 2)); + REQUIRE(vol == Catch::Approx(1.0 / 6.0).epsilon(1e-12)); + REQUIRE(vol == Catch::Approx(std::abs(r3d::measure(poly))).epsilon(1e-12)); +} + +TEST_CASE( + "ForEachPolytopeFaceTriangle: a clipped polyhedron with quad faces is " + "tiled exactly", + "[form_integrator_utils]") +{ + // Clip the reference tet with the plane x <= 0.5. This truncates the corner + // at (1,0,0), producing a polyhedron with a quadrilateral face, so at least + // one face must fan into more than one triangle. The removed piece is a tet + // similar to the original at scale 0.5 (volume (1/6)*0.5^3 = 1/48), leaving + // 1/6 - 1/48 = 7/48. + r3d::Polytope<3> poly; + r3d::init(poly, MakeTetVerts({{{0.0, 0.0, 0.0}}, + {{1.0, 0.0, 0.0}}, + {{0.0, 1.0, 0.0}}, + {{0.0, 0.0, 1.0}}})); + r3d::Few, 1> planes; + planes[0].n[0] = -1.0; // keep -x + 0.5 >= 0, i.e. x <= 0.5 + planes[0].n[1] = 0.0; + planes[0].n[2] = 0.0; + planes[0].d = 0.5; + r3d::clip(poly, planes); + + REQUIRE(poly.nverts > 4); // truncation added vertices / a quad face + const double measure = std::abs(r3d::measure(poly)); + REQUIRE(measure == Catch::Approx(7.0 / 48.0).epsilon(1e-12)); + + const auto [vol, ntri] = DecomposeFromCentroid(poly); + // Every face walked exactly once <=> the genus-0 fan invariant holds. + REQUIRE(ntri == 2 * (poly.nverts - 2)); + REQUIRE(ntri > 4); // a quad face fans into 2+ triangles + REQUIRE(vol == Catch::Approx(measure).epsilon(1e-12)); +}