From fcb88ba9687bf55e78027a9085f5cd83f4d55845 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Tue, 28 Jul 2026 23:34:40 -0400 Subject: [PATCH] apply meshFields 1.0.0 API --- .../field/evaluator/mesh_fields_backend.h | 60 +++++++++++++++++-- src/pcms/transfer/mass_matrix_integrator.hpp | 11 +++- .../omega_h_form_integrator_utils.hpp | 24 +++++++- src/pcms/transfer/omega_h_mass_integrator.cpp | 4 ++ test/test_omega_h_mass_integrator.cpp | 5 ++ 5 files changed, 98 insertions(+), 6 deletions(-) diff --git a/src/pcms/field/evaluator/mesh_fields_backend.h b/src/pcms/field/evaluator/mesh_fields_backend.h index 84eb206d..19e44d1f 100644 --- a/src/pcms/field/evaluator/mesh_fields_backend.h +++ b/src/pcms/field/evaluator/mesh_fields_backend.h @@ -6,6 +6,7 @@ #include #include +#include #include #include "pcms/field/layout/mesh_fields.h" @@ -52,8 +53,13 @@ class MeshFieldBackendImpl : public MeshFieldBackend Kokkos::View offsets) const override { auto self = const_cast*>(this); +#if MeshFields_VERSION < 10000 return self->mesh_field_.triangleLocalPointEval(localCoords, offsets, shape_field_); +#else + return self->mesh_field_.triangleLocalPointEval(localCoords, offsets, + shape_field_.field); +#endif } void SetData(Rank1View data, size_t num_nodes, @@ -65,7 +71,11 @@ class MeshFieldBackendImpl : public MeshFieldBackend mesh_.nents(dim), KOKKOS_CLASS_LAMBDA(size_t ent) { for (size_t n = 0; n < num_nodes; ++n) { for (size_t c = 0; c < num_components; ++c) { +#if MeshFields_VERSION < 10000 shape_field_(ent, n, c, topo) = +#else + shape_field_.field(ent, n, c, topo) = +#endif data[ent * stride + n * num_components + c]; } } @@ -82,7 +92,11 @@ class MeshFieldBackendImpl : public MeshFieldBackend for (size_t n = 0; n < num_nodes; ++n) { for (size_t c = 0; c < num_components; ++c) { data[ent * stride + n * num_components + c] = +#if MeshFields_VERSION < 10000 shape_field_(ent, n, c, topo); +#else + shape_field_.field(ent, n, c, topo); +#endif } } }); @@ -91,9 +105,8 @@ class MeshFieldBackendImpl : public MeshFieldBackend private: Omega_h::Mesh& mesh_; MeshField::OmegahMeshField mesh_field_; - using ShapeField = - decltype(mesh_field_.template CreateLagrangeField()); - ShapeField shape_field_; + using FWC = decltype(mesh_field_.template CreateLagrangeField()); + FWC shape_field_; // FWC = FieldWithController; keeps ctrlr alive }; // --------------------------------------------------------------------------- @@ -240,7 +253,11 @@ struct FillCoordinatesAndIndicesFunctor const auto owner_idx = owning_elem_ids_(i); LO count = Kokkos::atomic_sub_fetch(&elem_counts_(owner_idx), 1); LO index = offsets_(owner_idx) + count - 1; +#if MeshFields_VERSION < 10000 for (int j = 0; j < (dim_ + 1); ++j) { +#else + for (int j = 0; j < dim_; ++j) { +#endif coordinates_(index, j) = coord[j]; } indices_(index) = i; @@ -404,8 +421,18 @@ struct FillCoordinatesDeviceFunctor global_coords_(orig_idx, 1)}; const auto local = Omega_h::barycentric_from_global<2, 2>(point, vertex_coords); +#if MeshFields_VERSION < 10000 for (int j = 0; j < (dim_ + 1); ++j) coordinates_(index, j) = local[j]; +#else + // Computation in barycentric_from_global sometimes produce values + // slightly outside [0, 1] for points on or near element boundaries. + for (int j = 0; j < dim_; ++j) { + coordinates_(index, j) = local[j] < 0.0 ? 0.0 + : local[j] > 1.0 ? 1.0 + : local[j]; + } +#endif indices_(index) = orig_idx; } }; @@ -494,8 +521,13 @@ struct MeshFieldsAdapter2LocalizationHint } offsets_ = Kokkos::View("offsets", mesh.nelems() + 1); +#if MeshFields_VERSION < 10000 coordinates_ = Kokkos::View( "coordinates", num_valid_, mesh.dim() + 1); +#else + coordinates_ = Kokkos::View( + "coordinates", num_valid_, mesh.dim()); +#endif indices_ = Kokkos::View("indices", num_valid_); if (num_missing_ > 0) { @@ -537,8 +569,18 @@ struct MeshFieldsAdapter2LocalizationHint global_coords(orig_idx, 1)}; const auto local = Omega_h::barycentric_from_global<2, 2>(point, vertex_coords); +#if MeshFields_VERSION < 10000 for (int j = 0; j < (mesh.dim() + 1); ++j) coordinates_(index, j) = local[j]; +#else + // Computation in barycentric_from_global sometimes produce values + // slightly outside [0, 1] for points on or near element boundaries. + for (int j = 0; j < mesh.dim(); ++j) { + coordinates_(index, j) = local[j] < 0.0 ? 0.0 + : local[j] > 1.0 ? 1.0 + : local[j]; + } +#endif indices_(index) = static_cast(orig_idx); } @@ -631,8 +673,13 @@ struct MeshFieldsAdapter2LocalizationHint set_final_functor); // Step 6: Fill coordinates and indices on device +#if MeshFields_VERSION < 10000 coordinates_d_ = Kokkos::View("coordinates_d", num_valid_, mesh.dim() + 1); +#else + coordinates_d_ = + Kokkos::View("coordinates_d", num_valid_, mesh.dim()); +#endif indices_d_ = Kokkos::View("indices_d", num_valid_); const auto tris2verts = mesh.ask_elem_verts(); const auto mesh_coords = mesh.coords(); @@ -653,8 +700,13 @@ struct MeshFieldsAdapter2LocalizationHint // Create host mirrors for compatibility (lazy copy - only if needed) offsets_ = Kokkos::create_mirror_view(offsets_d_); +#if MeshFields_VERSION < 10000 coordinates_ = Kokkos::View( "coordinates_", num_valid_, mesh.dim() + 1); +#else + coordinates_ = Kokkos::View( + "coordinates_", num_valid_, mesh.dim()); +#endif DeepCopyMismatchLayouts(coordinates_, coordinates_d_); indices_ = Kokkos::create_mirror_view(indices_d_); if (num_missing_ > 0) { @@ -681,4 +733,4 @@ struct MeshFieldsAdapter2LocalizationHint } // namespace pcms -#endif // PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_BACKEND_H +#endif // PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_BACKEND_H \ No newline at end of file diff --git a/src/pcms/transfer/mass_matrix_integrator.hpp b/src/pcms/transfer/mass_matrix_integrator.hpp index 7d36c975..d44e513c 100644 --- a/src/pcms/transfer/mass_matrix_integrator.hpp +++ b/src/pcms/transfer/mass_matrix_integrator.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -38,7 +39,11 @@ class MassMatrixIntegrator : public MeshField::Integrator // std::cerr << " Number points per Elem : " << numPtsPerElem << "\n"; assert(numPtsPerElem >= 1); const size_t ptDim = p.extent(1); +#if MeshFields_VERSION < 10000 assert(ptDim == fe.MeshEntDim + 1); +#else + assert(ptDim == fe.MeshEntDim); +#endif // Copy values needed in the kernel to avoid capturing host references // (mesh and fe are host objects and cannot be dereferenced on the device) const auto numElems = mesh.nelems(); @@ -50,8 +55,12 @@ class MassMatrixIntegrator : public MeshField::Integrator const auto first = elm * numPtsPerElem; const auto last = first + numPtsPerElem; for (auto pt = first; pt < last; pt++) { - // FIXME better way to fill? pass kokkos::subview to getValues? + // FIXME better way to fill? pass kokkos::subview to getValues? +#if MeshFields_VERSION < 10000 Kokkos::Array +#else + Kokkos::Array +#endif localCoord; for (auto i = 0; i < localCoord.size(); i++) { localCoord[i] = p(pt, i); diff --git a/src/pcms/transfer/omega_h_form_integrator_utils.hpp b/src/pcms/transfer/omega_h_form_integrator_utils.hpp index 23f5e35d..b9f43047 100644 --- a/src/pcms/transfer/omega_h_form_integrator_utils.hpp +++ b/src/pcms/transfer/omega_h_form_integrator_utils.hpp @@ -7,6 +7,7 @@ #include "pcms/utility/assert.h" #include #include +#include #include #include @@ -72,14 +73,27 @@ inline void CheckOmegaHScalarLagrangeLayout( } [[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<2> GlobalFromBarycentric( +#if MeshFields_VERSION < 10000 const MeshField::Vector3& barycentric_coord, +#else + const MeshField::Vector2& barycentric_coord, +#endif const Omega_h::Few, 3>& verts_coord) { Omega_h::Vector<2> real_coords = {0.0, 0.0}; +#if MeshFields_VERSION < 10000 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]; } +#else + 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]; + } +#endif return real_coords; } @@ -250,15 +264,23 @@ OMEGA_H_INLINE void ForEachIntersectionSubtriangle( // source_order + target_order), which are only known at construction. struct IntegrationData { +#if MeshFields_VERSION < 10000 Kokkos::View bary_coords; // barycentric coordinates - Kokkos::View weights; // quadrature weights +#else + Kokkos::View bary_coords; // barycentric coordinates +#endif + Kokkos::View weights; // quadrature weights explicit IntegrationData(int order) { auto ip_vec = MeshField::getIntegrationPoints(order); const std::size_t num_ip = ip_vec.size(); +#if MeshFields_VERSION < 10000 bary_coords = Kokkos::View("bary_coords", num_ip); +#else + bary_coords = Kokkos::View("bary_coords", num_ip); +#endif weights = Kokkos::View("weights", num_ip); auto bary_coords_host = Kokkos::create_mirror_view(bary_coords); diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp index 29cd90e2..d0b10ee1 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.cpp +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -123,7 +123,11 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( omf(mesh); auto coordField = omf.getCoordField(); const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); +#if MeshFields_VERSION < 10000 MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); +#else + MeshField::FieldElement coordFe(mesh.nelems(), coordField.field, shp, map); +#endif auto elm_mass_dev = buildElementMassMatrix(mesh, coordFe); // Build COO sparsity pattern on device: each element contributes a 3x3 block. diff --git a/test/test_omega_h_mass_integrator.cpp b/test/test_omega_h_mass_integrator.cpp index 939283d1..a1de7e76 100644 --- a/test/test_omega_h_mass_integrator.cpp +++ b/test/test_omega_h_mass_integrator.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -60,7 +61,11 @@ std::map, pcms::Real> BuildReferenceMassMap( omf(mesh); auto coordField = omf.getCoordField(); const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); +#if MeshFields_VERSION < 10000 MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); +#else + MeshField::FieldElement coordFe(mesh.nelems(), coordField.field, shp, map); +#endif auto elm_mass_dev = buildElementMassMatrix(mesh, coordFe); auto elm_mass_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, elm_mass_dev);