diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d407a53..331e4e21 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,6 +49,23 @@ 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") +set(PCMS_INTERSECTION_ABS_TOL 1e-18 CACHE STRING + "Absolute measure tolerance for accepting a clipped intersection") +set(PCMS_INTERSECTION_REL_TOL 1e-12 CACHE STRING + "Relative measure tolerance for accepting a clipped intersection") +set(PCMS_PROJECTION_KSP_RTOL 1e-14 CACHE STRING + "Default relative tolerance for the projection mass-matrix solve") +set(PCMS_PROJECTION_KSP_ATOL 1e-14 CACHE STRING + "Default absolute tolerance for the projection mass-matrix solve") + 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..51c40b44 100644 --- a/src/pcms/configuration.h.in +++ b/src/pcms/configuration.h.in @@ -7,3 +7,11 @@ #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@ + +#cmakedefine PCMS_INTERSECTION_ABS_TOL @PCMS_INTERSECTION_ABS_TOL@ +#cmakedefine PCMS_INTERSECTION_REL_TOL @PCMS_INTERSECTION_REL_TOL@ +#cmakedefine PCMS_PROJECTION_KSP_RTOL @PCMS_PROJECTION_KSP_RTOL@ +#cmakedefine PCMS_PROJECTION_KSP_ATOL @PCMS_PROJECTION_KSP_ATOL@ 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")); diff --git a/src/pcms/localization/queue_visited.hpp b/src/pcms/localization/queue_visited.hpp index 8d768eec..bab16e69 100644 --- a/src/pcms/localization/queue_visited.hpp +++ b/src/pcms/localization/queue_visited.hpp @@ -7,8 +7,7 @@ #include #include #include -#define MAX_SIZE_QUEUE 500 -#define MAX_SIZE_TRACK 800 +#include "pcms/configuration.h" namespace pcms { @@ -16,7 +15,7 @@ namespace pcms class Queue { private: - Omega_h::LO queue_array[MAX_SIZE_QUEUE]; + Omega_h::LO queue_array[PCMS_INTERSECTION_QUEUE_SIZE]; int first = 0, last = -1, count = 0; public: @@ -45,7 +44,7 @@ class Queue class Track { private: - Omega_h::LO tracking_array[MAX_SIZE_TRACK]; + Omega_h::LO tracking_array[PCMS_INTERSECTION_TRACK_SIZE]; int first = 0, last = -1, count = 0; public: @@ -68,11 +67,11 @@ class Track OMEGA_H_INLINE void Queue::push_back(const int& item) { - if (count == MAX_SIZE_QUEUE) { + if (count == PCMS_INTERSECTION_QUEUE_SIZE) { printf("queue is full %d\n", count); return; } - last = (last + 1) % MAX_SIZE_QUEUE; + last = (last + 1) % PCMS_INTERSECTION_QUEUE_SIZE; queue_array[last] = item; count++; } @@ -84,7 +83,7 @@ void Queue::pop_front() printf("queue is empty\n"); return; } - first = (first + 1) % MAX_SIZE_QUEUE; + first = (first + 1) % PCMS_INTERSECTION_QUEUE_SIZE; count--; } @@ -103,19 +102,19 @@ bool Queue::isEmpty() const OMEGA_H_INLINE bool Queue::isFull() const { - return count == MAX_SIZE_QUEUE; + return count == PCMS_INTERSECTION_QUEUE_SIZE; } OMEGA_H_INLINE bool Track::push_back(const int& item) { - if (count == MAX_SIZE_TRACK) { + if (count == PCMS_INTERSECTION_TRACK_SIZE) { // Visited buffer is full. Report failure so callers can stop expanding the // search; otherwise new vertices can never be marked visited and the BFS // re-queues them forever (infinite loop). return false; } - last = (last + 1) % MAX_SIZE_TRACK; + last = (last + 1) % PCMS_INTERSECTION_TRACK_SIZE; tracking_array[last] = item; count++; return true; @@ -126,7 +125,7 @@ bool Track::notVisited(const int& item) { int id; for (int i = 0; i < count; ++i) { - id = (first + i) % MAX_SIZE_TRACK; + id = (first + i) % PCMS_INTERSECTION_TRACK_SIZE; if (tracking_array[id] == item) { return false; } diff --git a/src/pcms/transfer/conservative_projection_solver.cpp b/src/pcms/transfer/conservative_projection_solver.cpp index bfd63f03..2452a550 100644 --- a/src/pcms/transfer/conservative_projection_solver.cpp +++ b/src/pcms/transfer/conservative_projection_solver.cpp @@ -1,5 +1,6 @@ #include +#include "pcms/configuration.h" #include "pcms/transfer/conservative_projection_solver.hpp" #include "pcms/utility/arrays.h" #include "pcms/transfer/petsc_utils.hpp" @@ -65,6 +66,10 @@ GalerkinProjectionSolver::GalerkinProjectionSolver( ierr = PCSetType(pc, PCJACOBI); CHKERRABORT(PETSC_COMM_SELF, ierr); } + ierr = + KSPSetTolerances(ksp_, PCMS_PROJECTION_KSP_RTOL, PCMS_PROJECTION_KSP_ATOL, + PETSC_DEFAULT, PETSC_DEFAULT); + CHKERRABORT(PETSC_COMM_SELF, ierr); ierr = KSPSetFromOptions(ksp_); CHKERRABORT(PETSC_COMM_SELF, ierr); ierr = KSPSetUp(ksp_); diff --git a/src/pcms/transfer/mass_matrix_integrator.hpp b/src/pcms/transfer/mass_matrix_integrator.hpp index b3031ff0..472a540f 100644 --- a/src/pcms/transfer/mass_matrix_integrator.hpp +++ b/src/pcms/transfer/mass_matrix_integrator.hpp @@ -18,16 +18,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..66477305 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,45 @@ 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 PCMS_INTERSECTION_TRACK_SIZE if + // this fires. + if (!visited.push_back(neighborElmId)) { + printf("ERROR: visited buffer full " + "(PCMS_INTERSECTION_TRACK_SIZE=%d) for " + "target %d; some intersections may be missed\n", + PCMS_INTERSECTION_TRACK_SIZE, 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); - auto eps = Kokkos::fmax(abs_tol, rel_tol * scale); - if (intersection.nverts >= 3 && intersected_area >= eps) { + Kokkos::fmax(current_tgt_elm_measure, current_src_elm_measure); + auto eps = Kokkos::fmax(PCMS_INTERSECTION_ABS_TOL, + PCMS_INTERSECTION_REL_TOL * scale); + // 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 < PCMS_INTERSECTION_QUEUE_SIZE, + "intersection count for target %d reached the cap %d; raise " + "PCMS_INTERSECTION_QUEUE_SIZE", + id, PCMS_INTERSECTION_QUEUE_SIZE); queue.push_back(neighborElmId); @@ -113,20 +151,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 +189,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..e7cc17ff 100644 --- a/src/pcms/transfer/mesh_intersection.hpp +++ b/src/pcms/transfer/mesh_intersection.hpp @@ -1,6 +1,7 @@ #ifndef PCMS_TRANSFER_MESH_INTERSECTION_HPP #define PCMS_TRANSFER_MESH_INTERSECTION_HPP +#include #include #include #include @@ -12,22 +13,24 @@ 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> +// 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 +81,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 e40768f5..b3ac305a 100644 --- a/src/pcms/transfer/omega_h_form_integrator_utils.hpp +++ b/src/pcms/transfer/omega_h_form_integrator_utils.hpp @@ -14,9 +14,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) @@ -34,12 +35,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"); } } @@ -49,7 +51,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"); @@ -64,7 +66,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 + @@ -72,30 +74,23 @@ inline void CheckOmegaHScalarLagrangeLayout( } } -[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<2> GlobalFromBarycentric( - const MeshField::Vector2& 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}; - const Omega_h::Real xi3 = 1.0 - barycentric_coord[0] - barycentric_coord[1]; - const Omega_h::Real xi[3] = {barycentric_coord[0], barycentric_coord[1], xi3}; - for (int i = 0; i < 3; ++i) { - real_coords[0] += xi[i] * verts_coord[i][0]; - real_coords[1] += xi[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( @@ -191,22 +186,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); @@ -229,45 +226,216 @@ OMEGA_H_INLINE void ForEachIntersectionSubtriangle( Omega_h::Real area = Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); - const double eps_area = abs_tol + rel_tol * poly_area; + const double eps_area = + PCMS_INTERSECTION_ABS_TOL + PCMS_INTERSECTION_REL_TOL * poly_area; if (area <= eps_area) { 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 = + PCMS_INTERSECTION_ABS_TOL + PCMS_INTERSECTION_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; + // MeshField returns points in reduced parametric coordinates: only the + // first Dim barycentric components are stored, with the last implied by + // the partition of unity. Expand to the full Dim+1 barycentric form + // consumed by GlobalFromBarycentric. + Omega_h::Real last = 1.0; + for (int d = 0; d < Dim; ++d) { + const Omega_h::Real xi = ip_vec[i].param[d]; + bary_coords_host(i, d) = xi; + last -= xi; + } + bary_coords_host(i, Dim) = last; weights_host(i) = ip_vec[i].weight; } Kokkos::deep_copy(bary_coords, bary_coords_host); @@ -277,9 +445,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 @@ -287,40 +456,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]); } @@ -328,14 +496,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 b5baad53..742b8684 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 6b37554e..4538474a 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,38 +42,40 @@ 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. When +// 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. When // lumped, every block entry is mapped onto the row's diagonal instead; PETSc's // COO assembly sums repeated indices, so the unmodified element values // accumulate into the row-sum lumped diagonal. +template void FillP1MassCooPattern( - int nelems, const Omega_h::LOs& faces2nodes, + int nelems, const Omega_h::LOs& elems2nodes, const Kokkos::View& global_to_local, bool lumped, 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; const auto row = static_cast(global_to_local(verts[i])); coo_rows(idx) = row; coo_cols(idx) = @@ -69,95 +85,120 @@ void FillP1MassCooPattern( }); } -} // namespace - -OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space, - MassMatrixType mass_type) - : OmegaHMassIntegrator(std::dynamic_pointer_cast( - target_space.GetLayout()), - target_space.GetCoordinateSystem(), mass_type) -{ -} - -OmegaHMassIntegrator::OmegaHMassIntegrator( - std::shared_ptr target_layout, - CoordinateSystem coordinate_system, MassMatrixType mass_type) +// 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, + MassMatrixType mass_type) { - 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(); - diagonal_ = - target_layout->GetOrder() == 0 || mass_type == MassMatrixType::Lumped; + 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_SELF, num_dofs, num_dofs, 0, nullptr, &mat_); + createSeqAIJMat(PETSC_COMM_SELF, num_dofs, num_dofs, 0, nullptr, &mat); CHKERRABORT(PETSC_COMM_SELF, 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_SELF, ierr); - ierr = MatSetValuesCOO(mat_, vals.data(), INSERT_VALUES); + ierr = MatSetValuesCOO(mat, vals.data(), INSERT_VALUES); CHKERRABORT(PETSC_COMM_SELF, 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.field, 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.field, shp, map); + elm_mass_dev = buildElementMassMatrix(mesh, coordFe); + } else { + const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); + MeshField::FieldElement coordFe(mesh.nelems(), coordField.field, 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, - mass_type == MassMatrixType::Lumped, coo_rows, coo_cols); + FillP1MassCooPattern(nelems, elems2nodes, global_to_local, + mass_type == MassMatrixType::Lumped, 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_SELF, num_dofs, num_dofs, 0, nullptr, &mat_); + createSeqAIJMat(PETSC_COMM_SELF, num_dofs, num_dofs, 0, nullptr, &mat); CHKERRABORT(PETSC_COMM_SELF, 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_SELF, ierr); - ierr = MatSetValuesCOO(mat_, elm_mass_dev.data(), INSERT_VALUES); + ierr = MatSetValuesCOO(mat, elm_mass_dev.data(), INSERT_VALUES); CHKERRABORT(PETSC_COMM_SELF, ierr); + return mat; +} + +} // namespace + +OmegaHMassIntegrator::OmegaHMassIntegrator(const FunctionSpace& target_space, + MassMatrixType mass_type) + : OmegaHMassIntegrator(std::dynamic_pointer_cast( + target_space.GetLayout()), + target_space.GetCoordinateSystem(), mass_type) +{ +} + +OmegaHMassIntegrator::OmegaHMassIntegrator( + std::shared_ptr target_layout, + CoordinateSystem coordinate_system, MassMatrixType mass_type) +{ + detail::CheckOmegaHScalarLagrangeLayout(coordinate_system, target_layout, + "OmegaHMassIntegrator", "target"); + + Omega_h::Mesh& mesh = target_layout->GetMesh(); + diagonal_ = + target_layout->GetOrder() == 0 || mass_type == MassMatrixType::Lumped; + if (mesh.dim() == 3) { + mat_ = BuildOmegaHMassMatrixImpl<3>(mesh, *target_layout, mass_type); + } else { + mat_ = BuildOmegaHMassMatrixImpl<2>(mesh, *target_layout, mass_type); + } } 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 1684ebbb..8bfcaeba 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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 9a7eb845..b51395eb 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -443,7 +443,9 @@ 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 + test_omega_h_3d_analytic_mass.cpp) endif() add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES}) diff --git a/test/field_test_utils.h b/test/field_test_utils.h index b579866b..68becc51 100644 --- a/test/field_test_utils.h +++ b/test/field_test_utils.h @@ -44,6 +44,11 @@ KOKKOS_INLINE_FUNCTION Real linear_f(Real x, Real y) return x + 2.0 * y; } +KOKKOS_INLINE_FUNCTION Real linear_f_3d(Real x, Real y) +{ + return x + 2.0 * y; +} + // Interior test points for a unit [0,1]^2 box mesh. inline std::vector StandardEvalCoords2D() { @@ -60,6 +65,21 @@ inline std::vector StandardOutsideCoords2D() // 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. +// + +// Geometric classification tags required by OmegaHLagrangeLayout. +inline void AddDefaultClassification(Omega_h::Mesh& mesh) +{ + for (Omega_h::Int dim = 0; dim <= mesh.dim(); ++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))); + } +} + inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, const Omega_h::LOs& ev2v) { @@ -71,14 +91,7 @@ inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, }); 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))); - } + AddDefaultClassification(mesh); return mesh; } @@ -91,6 +104,44 @@ inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) : Omega_h::LOs({0, 1, 2, 0, 2, 3})); } +// Tetrahedral mesh of [0,1]^3. build_box already adds classification. +inline 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); +} + +inline Omega_h::Mesh BuildTet(Omega_h::Library& lib, + const Omega_h::Reals& coords) +{ + + Omega_h::Mesh mesh(&lib); + const Omega_h::LOs ev2v = Omega_h::LOs({0, 1, 2, 3}); + Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 3, ev2v, coords); + AddDefaultClassification(mesh); + return mesh; +} + +// Reference tet: (0,0,0), (1,0,0), (0,1,0), (0,0,1). Volume 1/6. +// reversed=true swaps the last two verts so the signed volume is negative. +inline Omega_h::Mesh BuildReferenceTet(Omega_h::Library& lib) +{ + return BuildTet(lib, Omega_h::Reals({ + 0.0, + 0.0, + 0.0, + 1.0, + 0.0, + 0.0, + 0.0, + 1.0, + 0.0, + 0.0, + 0.0, + 1.0, + })); +} + inline std::shared_ptr MakeP1Space( Omega_h::Mesh& mesh, const std::string& global_id_name = "global") { 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_analytic_mass.cpp b/test/test_omega_h_3d_analytic_mass.cpp new file mode 100644 index 00000000..a27ef85f --- /dev/null +++ b/test/test_omega_h_3d_analytic_mass.cpp @@ -0,0 +1,134 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "field_test_utils.h" +#include + +// Reference corner tet: +// v0=(0,0,0), v1=(1,0,0), v2=(0,1,0), v3=(0,0,1) +// Volume V = 1/6. +// +// P1 consistent mass on a tet: +// M_ij = V/20 * (1 + delta_ij) +// = 1/60 on diagonal, 1/120 off-diagonal. + +TEST_CASE("OmegaHMassIntegrator (3D): P1 mass on reference tet matches " + "analytic V/20*(1+delta)", + "[mass_integrator][3d][analytic]") +{ + Omega_h::Library lib; + auto mesh = pcms::test::BuildReferenceTet(lib); + REQUIRE(mesh.dim() == 3); + REQUIRE(mesh.nverts() == 4); + REQUIRE(mesh.nelems() == 1); + + constexpr pcms::Real V = 1.0 / 6.0; + constexpr pcms::Real M_diag = V / 10.0; // 1/60 + constexpr pcms::Real M_off = V / 20.0; // 1/120 + + auto space = pcms::test::MakeP1Space(mesh); + auto integrator = pcms::BuildOmegaHMassIntegrator(*space); + Mat mat = integrator->GetMatrix(); + + const auto layout = + std::dynamic_pointer_cast( + space->GetLayout()); + REQUIRE(layout != nullptr); + + // PETSc rows are active indices = global_to_local(local_vertex). + const auto perm = layout->GetGlobalToLocalPermutationHost(); + REQUIRE(static_cast(perm.extent(0)) == 4); + + pcms::Real grand_total = 0.0; + for (int i = 0; i < 4; ++i) { + pcms::Real row_sum = 0.0; + for (int j = 0; j < 4; ++j) { + PetscInt r = static_cast(perm(i)); + PetscInt c = static_cast(perm(j)); + PetscScalar got = 0.0; + MatGetValues(mat, 1, &r, 1, &c, &got); + const pcms::Real expected = (i == j) ? M_diag : M_off; + CAPTURE(i, j, r, c, expected, got); + CHECK(static_cast(got) == + Catch::Approx(expected).epsilon(1e-12)); + row_sum += static_cast(got); + } + // Row sum =\int \lamda_i = V/4 + CHECK(row_sum == Catch::Approx(V / 4.0).epsilon(1e-12)); + grand_total += row_sum; + } + // Sum of row sums = volume + CHECK(grand_total == Catch::Approx(V).epsilon(1e-12)); +} + +TEST_CASE("OmegaHMassIntegrator (3D): P0 mass on reference tet equals volume", + "[mass_integrator][3d][analytic]") +{ + Omega_h::Library lib; + auto mesh = pcms::test::BuildReferenceTet(lib); + + constexpr pcms::Real V = 1.0 / 6.0; + + auto space = pcms::test::MakeP0Space(mesh); + auto integrator = pcms::BuildOmegaHMassIntegrator(*space); + Mat mat = integrator->GetMatrix(); + + PetscInt r = 0; + PetscInt c = 0; + PetscScalar got = 0.0; + MatGetValues(mat, 1, &r, 1, &c, &got); + CHECK(static_cast(got) == Catch::Approx(V).epsilon(1e-12)); +} + +TEST_CASE("OmegaHConservativeProjection (3D): same-mesh reference tet is " + "exact for constant and linear fields", + "[transfer][mesh_intersection][3d][analytic]") +{ + // Purpose: source == target == one tet. Intersection is trivial; P1 must + // reproduce constants/linears at vertices. If this fails, bug is in + // RHS/mass/Apply/KSP — not multi-tet clipping. + Omega_h::Library lib; + Omega_h::Mesh mesh = pcms::test::BuildReferenceTet(lib); + REQUIRE(mesh.nverts() == 4); + REQUIRE(mesh.nelems() == 1); + auto space = pcms::test::MakeP1Space(mesh); + auto source = space->CreateFunction(); + auto target = space->CreateFunction(); + pcms::OmegaHConservativeProjection projection(*space, *space); + SECTION("constant field") + { + const double c = 2.5; + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real, pcms::Real, pcms::Real) { return c; }); + projection.Apply(source, target); + const auto values = pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + REQUIRE(static_cast(values.size()) == 4); + for (Omega_h::LO i = 0; i < 4; ++i) { + CAPTURE(i, values[i]); + REQUIRE(values[i] == Catch::Approx(c).margin(1e-12)); + } + } + SECTION("linear field") + { + 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 values = pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); + const auto coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(mesh.coords(), 3), mesh.nverts(), 3); + for (Omega_h::LO i = 0; i < 4; ++i) { + const double expected = + 1.0 + coords_h(i, 0) + 2.0 * coords_h(i, 1) + 3.0 * coords_h(i, 2); + CAPTURE(i, expected, values[i]); + REQUIRE(values[i] == Catch::Approx(expected).margin(1e-12)); + } + } +} 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..9b037845 --- /dev/null +++ b/test/test_omega_h_3d_conservative_projection.cpp @@ -0,0 +1,213 @@ +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include "field_test_utils.h" + +#include + +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 = pcms::test::BuildUnitCube(lib, 1); + Omega_h::Mesh target_mesh = pcms::test::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 = pcms::test::BuildUnitCube(lib, 2); + Omega_h::Mesh target_mesh = pcms::test::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 = pcms::test::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 = pcms::test::BuildUnitCube(lib, 2); + Omega_h::Mesh target_mesh = pcms::test::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 = pcms::test::BuildUnitCube(lib, 1); + Omega_h::Mesh target_mesh = pcms::test::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)); +} diff --git a/test/test_omega_h_intersection_rhs_integrator.cpp b/test/test_omega_h_intersection_rhs_integrator.cpp index f5d667d9..f645d6d3 100644 --- a/test/test_omega_h_intersection_rhs_integrator.cpp +++ b/test/test_omega_h_intersection_rhs_integrator.cpp @@ -258,3 +258,151 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts", pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space)); } } + +TEST_CASE("OmegaHIntersectionRHSIntegrator (3D): constant load sums to c times " + "known tet-tet overlap", + "[rhs_integrator][3d]") +{ + // Source is the reference tet (volume 1/6). Target is the similar tet at the + // (1,0,0) corner, cut by the plane x=0.5: + // (1,0,0), (0.5,0,0), (0.5,0.5,0), (0.5,0,0.5) + // which sits entirely inside the source (volume (1/6)*(1/2)^3 = 1/48). + // checks PCMS search + quadrature + P1 partition of unity: + // sum_j b_j = \int_overlap f dx = c * 1/48. + Omega_h::Library lib; + auto source_mesh = pcms::test::BuildReferenceTet(lib); + auto target_mesh = pcms::test::BuildTet(lib, Omega_h::Reals({ + 1.0, + 0.0, + 0.0, + 0.5, + 0.0, + 0.0, + 0.5, + 0.5, + 0.0, + 0.5, + 0.0, + 0.5, + })); + + REQUIRE(source_mesh.nelems() == 1); + REQUIRE(target_mesh.nelems() == 1); + + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); + + constexpr double c = 3.0; + constexpr double overlap = 1.0 / 48.0; + auto source_field = source_space->CreateFunction(); + pcms::test::SetField( + source_field, + KOKKOS_LAMBDA(pcms::Real, pcms::Real, pcms::Real) { return c; }); + + auto integrator = + pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space); + + const auto raw_coords = integrator->GetIntegrationPoints().GetValues(); + REQUIRE(raw_coords.extent(0) > 0); + REQUIRE(raw_coords.extent(1) == 3); + auto coords_h = pcms::test::CopyCoordinatesToHost( + raw_coords, static_cast(raw_coords.extent(0)), + static_cast(raw_coords.extent(1))); + for (std::size_t i = 0; i < coords_h.extent(0); ++i) { + const double x = coords_h(i, 0); + const double y = coords_h(i, 1); + const double z = coords_h(i, 2); + CAPTURE(i, x, y, z); + CHECK(x >= 0.5 - 1e-12); + CHECK(y >= -1e-12); + CHECK(z >= -1e-12); + CHECK(x + y + z <= 1.0 + 1e-12); + } + + pcms::test::EvaluateAndAssemble(*integrator, source_space, source_field); + PetscScalar sum = 0.0; + VecSum(integrator->GetVector(), &sum); + CHECK(static_cast(sum) == + Catch::Approx(c * overlap).margin(1e-12)); +} + +TEST_CASE("OmegaHIntersectionRHSIntegrator (3D): dual tets overlap is an " + "octahedron; load sums to integral over the intersection", + "[rhs_integrator][3d]") +{ + // Dual tets in the cube [-1,1]^3. Intersection = {|x|+|y|+|z| <= 1}, + // volume 4/3. + Omega_h::Library lib; + auto source_mesh = pcms::test::BuildTet(lib, Omega_h::Reals({ + 1.0, + 1.0, + 1.0, + 1.0, + -1.0, + -1.0, + -1.0, + -1.0, + 1.0, + -1.0, + 1.0, + -1.0, + })); + auto target_mesh = pcms::test::BuildTet(lib, Omega_h::Reals({ + -1.0, + -1.0, + -1.0, + -1.0, + 1.0, + 1.0, + 1.0, + -1.0, + 1.0, + 1.0, + 1.0, + -1.0, + })); + + 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(); + REQUIRE(raw_coords.extent(0) > 0); + auto coords_h = pcms::test::CopyCoordinatesToHost( + raw_coords, static_cast(raw_coords.extent(0)), + static_cast(raw_coords.extent(1))); + for (std::size_t i = 0; i < coords_h.extent(0); ++i) { + const double s = std::abs(coords_h(i, 0)) + std::abs(coords_h(i, 1)) + + std::abs(coords_h(i, 2)); + CAPTURE(i, coords_h(i, 0), coords_h(i, 1), coords_h(i, 2), s); + CHECK(s <= 1.0 + 1e-10); + } + + SECTION("constant f=c => sum(b) = c * 4/3") + { + constexpr double c = 3.0; + auto source_field = source_space->CreateFunction(); + pcms::test::SetField( + source_field, + KOKKOS_LAMBDA(pcms::Real, pcms::Real, pcms::Real) { return c; }); + pcms::test::EvaluateAndAssemble(*integrator, source_space, source_field); + PetscScalar sum = 0.0; + VecSum(integrator->GetVector(), &sum); + CHECK(static_cast(sum) == + Catch::Approx(c * 4.0 / 3.0).margin(1e-10)); + } + + SECTION("f=x+y+z => sum(b) = 0 by symmetry") + { + auto source_field = source_space->CreateFunction(); + pcms::test::SetField( + source_field, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return x + y + z; + }); + pcms::test::EvaluateAndAssemble(*integrator, source_space, source_field); + PetscScalar sum = 0.0; + VecSum(integrator->GetVector(), &sum); + CHECK(static_cast(sum) == Catch::Approx(0.0).margin(1e-10)); + } +} diff --git a/test/test_omega_h_mc_rhs_integrator.cpp b/test/test_omega_h_mc_rhs_integrator.cpp index c29f4ec0..187a4d96 100644 --- a/test/test_omega_h_mc_rhs_integrator.cpp +++ b/test/test_omega_h_mc_rhs_integrator.cpp @@ -181,3 +181,120 @@ TEST_CASE("OmegaHControlVariateProjection: reduces error vs plain Monte Carlo", CAPTURE(mc_error, cv_error); CHECK(cv_error < mc_error); } + +// --------------------------------------------------------------------------- +// 3D Monte Carlo RHS integrator +// --------------------------------------------------------------------------- + +TEST_CASE("OmegaHMonteCarloRHSIntegrator (3D): sample points lie inside the " + "unit cube", + "[mc_rhs_integrator][3d]") +{ + // for Dim=3 sample coordinates inside [0,1]^3 with correct count. + Omega_h::Library lib; + auto target_mesh = pcms::test::BuildUnitCube(lib, 1); + auto target_space = pcms::test::MakeP1Space(target_mesh); + + const int samples_per_element = 16; + pcms::OmegaHMonteCarloRHSIntegrator integrator( + *target_space, samples_per_element, + pcms::MonteCarloSampling::UniformRandom); + + const auto raw_coords = integrator.GetIntegrationPoints().GetValues(); + REQUIRE(raw_coords.extent(1) == 3); + + 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) == + static_cast(target_mesh.nelems() * samples_per_element)); + + for (std::size_t i = 0; i < coords_h.extent(0); ++i) { + CAPTURE(i, coords_h(i, 0), coords_h(i, 1), coords_h(i, 2)); + CHECK(coords_h(i, 0) >= -1e-12); + CHECK(coords_h(i, 0) <= 1.0 + 1e-12); + CHECK(coords_h(i, 1) >= -1e-12); + CHECK(coords_h(i, 1) <= 1.0 + 1e-12); + CHECK(coords_h(i, 2) >= -1e-12); + CHECK(coords_h(i, 2) <= 1.0 + 1e-12); + } +} + +TEST_CASE("OmegaHMonteCarloRHSIntegrator (3D): constant field integrates " + "exactly", + "[mc_rhs_integrator][3d]") +{ + // if f=c => VecSum(b) = c * volume = c * 1 (P1 partition of unity). + Omega_h::Library lib; + auto source_mesh = pcms::test::BuildUnitCube(lib, 1); + auto target_mesh = pcms::test::BuildUnitCube(lib, 1); + 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( + source_field, + KOKKOS_LAMBDA(pcms::Real, pcms::Real, pcms::Real) { return 2.0; }); + + pcms::OmegaHMonteCarloRHSIntegrator integrator( + *target_space, /*samples_per_element=*/8, + pcms::MonteCarloSampling::UniformRandom); + + pcms::test::EvaluateAndAssemble(integrator, source_space, source_field); + + PetscScalar sum = 0.0; + VecSum(integrator.GetVector(), &sum); + CHECK(static_cast(sum) == Catch::Approx(2.0).margin(1e-12)); +} + +TEST_CASE( + "OmegaHControlVariateProjection (3D): reduces error vs plain Monte Carlo", + "[mc_rhs_integrator][control_variate][3d]") +{ + Omega_h::Library lib; + auto source_mesh = pcms::test::BuildUnitCube(lib, 1); + auto target_mesh = pcms::test::BuildUnitCube(lib, 2); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); + + auto source = source_space->CreateFunction(); + pcms::test::SetField( + source, KOKKOS_LAMBDA(pcms::Real x, pcms::Real y, pcms::Real z) { + return x * x + y * y + z * z; + }); + + auto reference = target_space->CreateFunction(); + pcms::OmegaHConservativeProjection ref_proj(*source_space, *target_space); + ref_proj.Apply(source, reference); + const auto ref_vals = + pcms::FlattenToRank1View(reference.GetDOFHolderDataHost()); + + const int nsample = 64; + const uint64_t seed = 20240611; + + pcms::OmegaHMonteCarloRHSIntegrator mc( + *target_space, nsample, pcms::MonteCarloSampling::UniformRandom, seed); + auto ev = source_space->CreatePointEvaluator( + pcms::EvaluationRequest::FromCoordinates(mc.GetIntegrationPoints())); + pcms::OmegaHMassIntegrator mass(*target_space); + pcms::GalerkinProjectionSolver solver(mass, mc); + const auto mc_x = solver.Solve(*ev, source); + const auto mc_h = Omega_h::HostRead(mc_x); + + auto cv_field = target_space->CreateFunction(); + pcms::OmegaHControlVariateProjection cv( + *source_space, *target_space, nsample, + pcms::MonteCarloSampling::UniformRandom, seed); + cv.Apply(source, cv_field); + const auto cv_vals = + pcms::FlattenToRank1View(cv_field.GetDOFHolderDataHost()); + + double mc_err = 0.0, cv_err = 0.0; + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + mc_err = std::max(mc_err, std::abs(mc_h[i] - ref_vals[i])); + cv_err = std::max(cv_err, std::abs(cv_vals[i] - ref_vals[i])); + } + CAPTURE(mc_err, cv_err); + CHECK(cv_err < mc_err); +}