diff --git a/test/field_test_utils.h b/test/field_test_utils.h
index 5aaac78e..2ddd2f45 100644
--- a/test/field_test_utils.h
+++ b/test/field_test_utils.h
@@ -3,6 +3,7 @@
 
 #include <catch2/catch_approx.hpp>
 #include <Kokkos_Core.hpp>
+#include "pcms/configuration.h"
 #include "pcms/field/field.h"
 #include "pcms/field/field_data.h"
 #include "pcms/field/field_evaluator_factory.h"
@@ -11,10 +12,22 @@
 #include "pcms/field/out_of_bounds_policy.h"
 #include "pcms/coupler/field_serializer.h"
 #include "pcms/field/coordinate_system.h"
-#include "pcms/localization/adj_search.hpp"
 #include "pcms/utility/arrays.h"
 #include "pcms/utility/memory_spaces.h"
+#ifdef PCMS_ENABLE_OMEGA_H
+#include <Omega_h_build.hpp>
+#include <Omega_h_mesh.hpp>
+#include <Omega_h_shape.hpp>
+#include "pcms/field/function_space/lagrange.h"
+#include "pcms/localization/adj_search.hpp"
+#endif
+#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS)
+#include "pcms/transfer/linear_form_integrator.hpp"
+#endif
 #include <cmath>
+#include <memory>
+#include <string>
+#include <type_traits>
 #include <unordered_map>
 #include <vector>
 
@@ -36,6 +49,62 @@ inline std::vector<Real> StandardEvalCoords2D()
   return {0.1, 0.2, 0.5, 0.5, 0.7, 0.3, 0.9, 0.1, 0.2, 0.8};
 }
 
+// Points strictly outside a unit [0,1]^2 box mesh.
+inline std::vector<Real> StandardOutsideCoords2D()
+{
+  return {-0.5, 0.5, 1.5, 0.5, 0.5, -0.5, 0.5, 1.5};
+}
+
+#ifdef PCMS_ENABLE_OMEGA_H
+// Builds a unit-square 2D simplex mesh from the given element connectivity and
+// adds the geometric classification tags required to build an Omega_h-backed
+// Lagrange function space.
+inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib,
+                                     const Omega_h::LOs& ev2v)
+{
+  const Omega_h::Reals coords({
+    0.0, 0.0, // v0
+    1.0, 0.0, // v1
+    1.0, 1.0, // v2
+    0.0, 1.0  // v3
+  });
+  Omega_h::Mesh mesh(&lib);
+  Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords);
+  for (Omega_h::Int dim = 0; dim <= 2; ++dim) {
+    mesh.add_tag<Omega_h::I8>(
+      dim, "class_dim", 1,
+      Omega_h::Read<Omega_h::I8>(mesh.nents(dim), Omega_h::I8(dim)));
+    mesh.add_tag<Omega_h::ClassId>(
+      dim, "class_id", 1,
+      Omega_h::Read<Omega_h::ClassId>(mesh.nents(dim), Omega_h::ClassId(0)));
+  }
+  return mesh;
+}
+
+// Convenience overload: diagonal=0 splits the square along vertices (1,3),
+// diagonal=1 along (0,2).
+inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal)
+{
+  return BuildUnitSquare(lib, (diagonal == 0)
+                                ? Omega_h::LOs({0, 1, 3, 1, 2, 3})
+                                : Omega_h::LOs({0, 1, 2, 0, 2, 3}));
+}
+
+inline std::shared_ptr<LagrangeFunctionSpace> MakeP1Space(
+  Omega_h::Mesh& mesh, const std::string& global_id_name = "global")
+{
+  return LagrangeFunctionSpace::FromMesh(
+    mesh, 1, 1, CoordinateSystem::Cartesian, global_id_name,
+    LagrangeFunctionSpace::Backend::OmegaH);
+}
+
+inline std::shared_ptr<LagrangeFunctionSpace> MakeP0Space(Omega_h::Mesh& mesh)
+{
+  return LagrangeFunctionSpace::FromMesh(
+    mesh, 0, 1, CoordinateSystem::Cartesian, "global",
+    LagrangeFunctionSpace::Backend::OmegaH);
+}
+
 inline bool AreArraysEqualUnordered(
   const Omega_h::HostRead<Omega_h::LO>& array1,
   const Omega_h::HostRead<Omega_h::LO>& array2, int start, int end)
@@ -76,6 +145,40 @@ inline std::vector<Real> CopyOmegaHRealsToVector(const Omega_h::Reals& coords)
                            coords_read.data() + coords_read.size());
 }
 
+inline double IntegrateP0Field(Omega_h::Mesh& mesh, const Field<Real>& field)
+{
+  const auto values = FlattenToRank1View(field.GetDOFHolderDataHost());
+  const auto measures = Omega_h::measure_elements_real(&mesh);
+  const auto measures_h = Omega_h::HostRead<Omega_h::Real>(measures);
+
+  double integral = 0.0;
+  for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) {
+    integral += measures_h[e] * values[e];
+  }
+  return integral;
+}
+
+inline double IntegrateP1Field(Omega_h::Mesh& mesh, const Field<Real>& field)
+{
+  const auto values = FlattenToRank1View(field.GetDOFHolderDataHost());
+  const auto measures = Omega_h::measure_elements_real(&mesh);
+  const auto measures_h = Omega_h::HostRead<Omega_h::Real>(measures);
+  const auto elem_verts_h =
+    Omega_h::HostRead<Omega_h::LO>(mesh.ask_elem_verts());
+  const int verts_per_elem = mesh.dim() + 1;
+
+  double integral = 0.0;
+  for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) {
+    double avg = 0.0;
+    for (int k = 0; k < verts_per_elem; ++k) {
+      avg += values[elem_verts_h[verts_per_elem * e + k]];
+    }
+    integral += measures_h[e] * (avg / verts_per_elem);
+  }
+  return integral;
+}
+#endif // PCMS_ENABLE_OMEGA_H
+
 // Copy coordinates from device memory to a host view.
 // This handles potential layout mismatches between host and device memory
 // spaces.
@@ -86,14 +189,8 @@ inline Kokkos::View<const Real**, HostMemorySpace> CopyCoordinatesToHost(
   auto coords_view =
     Kokkos::View<Real**, HostMemorySpace>("coords_view", nents, dim);
   auto coords_view_device =
-    Kokkos::create_mirror(DeviceMemorySpace(), coords_view);
-  Kokkos::parallel_for(
-    "copy_coords_to_host_view", Kokkos::RangePolicy<>(0, nents),
-    KOKKOS_LAMBDA(int i) {
-      for (int d = 0; d < dim; ++d) {
-        coords_view_device(i, d) = coords_device(i, d);
-      }
-    });
+    Kokkos::create_mirror_view(DeviceMemorySpace(), coords_view);
+  ConvertMismatchLayoutView2D(coords_view_device, coords_device);
   Kokkos::deep_copy(coords_view, coords_view_device);
   return coords_view;
 }
@@ -122,44 +219,49 @@ inline std::vector<Real> EvaluateReferenceFunction(const std::vector<Real>& pts,
                            expected_host.data() + expected_host.extent(0));
 }
 
-template <typename MemorySpace, typename Func>
+template <typename DataView, typename CoordsView, typename Func>
 struct SetFieldFunctor
 {
-  Kokkos::View<Real*, MemorySpace> data;
-  Kokkos::View<Real* [2], MemorySpace> coords;
+  DataView data;
+  CoordsView coords;
   Func f;
 
-  SetFieldFunctor(Kokkos::View<Real*, MemorySpace> data_,
-                  Kokkos::View<Real* [2], MemorySpace> coords_, Func f_)
+  SetFieldFunctor(DataView data_, CoordsView coords_, Func f_)
     : data(data_), coords(coords_), f(f_)
   {
   }
 
   KOKKOS_INLINE_FUNCTION
-  void operator()(int i) const { data(i) = f(coords(i, 0), coords(i, 1)); }
+  void operator()(int i) const
+  {
+    if constexpr (std::is_invocable_v<Func, Real, Real, Real>) {
+      data(i) = f(coords(i, 0), coords(i, 1), coords(i, 2));
+    } else {
+      data(i) = f(coords(i, 0), coords(i, 1));
+    }
+  }
 };
 
-// Set scalar DOF data by sampling func at each DOF-holder coordinate.
+// Set scalar DOF data by sampling func at each DOF-holder coordinate. The
+// arity of func selects the spatial dimension: func(x, y) for 2D layouts,
+// func(x, y, z) for 3D.
 template <typename ExecutionSpace = DefaultExecutionSpace, typename Func>
 inline void SetField(const FieldLayout& layout, FieldData<Real>& field,
                      Func func)
 {
   using MemorySpace = typename ExecutionSpace::memory_space;
 
+  static_assert(std::is_invocable_v<Func, Real, Real> ||
+                  std::is_invocable_v<Func, Real, Real, Real>,
+                "SetField requires func(x, y) or func(x, y, z)");
+
   auto dof_coords = layout.GetDOFHolderCoordinates().GetValues();
   int n = static_cast<int>(dof_coords.extent(0));
-  Kokkos::View<Real* [2], MemorySpace> coords_device("coords_device", n);
-  Kokkos::parallel_for(
-    "field_test_utils_copy_coords", Kokkos::RangePolicy<ExecutionSpace>(0, n),
-    KOKKOS_LAMBDA(int i) {
-      coords_device(i, 0) = dof_coords(i, 0);
-      coords_device(i, 1) = dof_coords(i, 1);
-    });
 
   Kokkos::View<Real*, MemorySpace> data_device("data_device", n);
   Kokkos::parallel_for("field_test_utils_set_field",
                        Kokkos::RangePolicy<ExecutionSpace>(0, n),
-                       SetFieldFunctor{data_device, coords_device, func});
+                       SetFieldFunctor{data_device, dof_coords, func});
 
   auto data_host =
     Kokkos::create_mirror_view_and_copy(HostMemorySpace(), data_device);
@@ -253,12 +355,8 @@ void CheckEvaluation(const PointEvaluator<Real>& evaluator,
 {
   int n = static_cast<int>(pts.size()) / 2;
 
-  Kokkos::View<Real*, DeviceMemorySpace> out_device("out_device", n);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<DeviceMemorySpace>;
-  Rank2View<Real, DeviceMemorySpace, LayoutPolicy> out(out_device.data(), n, 1);
-  // Rank2View<Real, HostMemorySpace> out(eval.data(), n, 1);
-  evaluator.Evaluate(field, out);
+  Kokkos::View<Real**, DeviceMemorySpace> out_device("out_device", n, 1);
+  evaluator.Evaluate(field, MakeRank2View(out_device));
   auto out_host =
     Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device);
 
@@ -267,8 +365,8 @@ void CheckEvaluation(const PointEvaluator<Real>& evaluator,
   for (int i = 0; i < n; ++i) {
     INFO("Point " << i << " (" << pts[2 * static_cast<size_t>(i)] << ", "
                   << pts[2 * static_cast<size_t>(i) + 1] << ")"
-                  << " got=" << out_host(i) << " expected=" << expected[i]);
-    REQUIRE(out_host(i) == Catch::Approx(expected[i]).margin(abs_tol));
+                  << " got=" << out_host(i, 0) << " expected=" << expected[i]);
+    REQUIRE(out_host(i, 0) == Catch::Approx(expected[i]).margin(abs_tol));
   }
 }
 
@@ -295,16 +393,13 @@ inline void CheckFillMode(const PointEvaluator<Real>& evaluator,
 {
   int n = static_cast<int>(outside_pts.size()) / 2;
 
-  Kokkos::View<Real*, DeviceMemorySpace> out_device("out_device", n);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<DeviceMemorySpace>;
-  Rank2View<Real, DeviceMemorySpace, LayoutPolicy> out(out_device.data(), n, 1);
-  evaluator.Evaluate(field, out);
+  Kokkos::View<Real**, DeviceMemorySpace> out_device("out_device", n, 1);
+  evaluator.Evaluate(field, MakeRank2View(out_device));
   auto out_host =
     Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device);
 
   for (int i = 0; i < n; ++i) {
-    REQUIRE(out_host(i) == fill_value);
+    REQUIRE(out_host(i, 0) == fill_value);
   }
 }
 
@@ -339,11 +434,8 @@ void CheckEvaluationWithFill(const Factory& factory, const Field<Real>& field,
   auto evaluator = factory->template CreatePointEvaluator<Real>(
     EvaluationRequest::FromCoordinates(device_coords.coordinate_view, policy));
 
-  Kokkos::View<Real*, DeviceMemorySpace> out_device("out_device", n);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<DeviceMemorySpace>;
-  Rank2View<Real, DeviceMemorySpace, LayoutPolicy> out(out_device.data(), n, 1);
-  evaluator->Evaluate(field, out);
+  Kokkos::View<Real**, DeviceMemorySpace> out_device("out_device", n, 1);
+  evaluator->Evaluate(field, MakeRank2View(out_device));
   auto out_host =
     Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device);
 
@@ -352,16 +444,34 @@ void CheckEvaluationWithFill(const Factory& factory, const Field<Real>& field,
   for (int i = 0; i < n; ++i) {
     INFO("Point " << i << " (" << pts[2 * static_cast<size_t>(i)] << ", "
                   << pts[2 * static_cast<size_t>(i) + 1] << ")"
-                  << " got=" << out_host(i));
+                  << " got=" << out_host(i, 0));
     if (is_inside[i]) {
       INFO(" expected=" << expected[i]);
-      REQUIRE(out_host(i) == Catch::Approx(expected[i]).margin(abs_tol));
+      REQUIRE(out_host(i, 0) == Catch::Approx(expected[i]).margin(abs_tol));
     } else {
-      REQUIRE(out_host(i) == fill_value);
+      REQUIRE(out_host(i, 0) == fill_value);
     }
   }
 }
 
+#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS)
+// Evaluates source_field at the integrator's sample points and assembles the
+// load vector.
+inline void EvaluateAndAssemble(
+  LinearFormIntegrator& integrator,
+  const std::shared_ptr<const FunctionSpace>& source_space,
+  const Field<Real>& source_field)
+{
+  const auto& pts = integrator.GetIntegrationPoints();
+  const std::size_t npts = pts.GetValues().extent(0);
+  auto evaluator = source_space->CreatePointEvaluator<Real>(
+    EvaluationRequest::FromCoordinates(pts));
+  Kokkos::View<Real**, DeviceMemorySpace> sampled("sampled", npts, 1);
+  evaluator->Evaluate(source_field, MakeRank2View(sampled));
+  integrator.Assemble(MakeConstRank2View(sampled));
+}
+#endif // PCMS_ENABLE_PETSC && PCMS_ENABLE_MESHFIELDS
+
 } // namespace pcms::test
 
 #endif // PCMS_TEST_FIELD_TEST_UTILS_H
diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp
index d546fbbe..c51f46b1 100644
--- a/test/test_eqdsk.cpp
+++ b/test/test_eqdsk.cpp
@@ -7,6 +7,7 @@
 #include <Kokkos_Core.hpp>
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/catch_approx.hpp>
+#include "field_test_utils.h"
 #include <iostream>
 #include <cmath>
 
@@ -212,42 +213,24 @@ TEST_CASE("EQDSKData with SplineFunctionSpace")
 
     const int num_eval_points = 3;
 
-    auto eval_coords_host = Kokkos::View<pcms::Real**, pcms::HostMemorySpace>(
-      "eval_coords_host", num_eval_points, 2);
-    for (size_t i = 0; i < num_eval_points; ++i) {
-      eval_coords_host(i, 0) = eval_coords[2 * i];     // R
-      eval_coords_host(i, 1) = eval_coords[2 * i + 1]; // Z
-    }
-
-    auto eval_coords_device =
-      Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace>("eval_coords_device",
-                                                          num_eval_points, 2);
-    pcms::DeepCopyMismatchLayouts(eval_coords_device, eval_coords_host);
-
-    auto coords_view = pcms::MakeRank2View(eval_coords_device);
-    auto coord_view = pcms::CoordinateView<pcms::DeviceMemorySpace>{
-      CoordinateSystem::Cartesian, coords_view};
-    auto eval_request = pcms::EvaluationRequest::FromCoordinates(coord_view);
+    auto device_coords = pcms::test::CreateDeviceCoordinateView(
+      eval_coords, CoordinateSystem::Cartesian);
+    auto eval_request =
+      pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view);
 
     auto evaluator =
       spline_space->CreatePointEvaluator<pcms::Real>(eval_request);
 
-    auto eval_results_1d = Kokkos::View<pcms::Real*, pcms::DeviceMemorySpace>(
-      "eval_results", num_eval_points);
-
-    using LayoutPolicy =
-      pcms::detail::default_layout_for_memory_space_t<pcms::DeviceMemorySpace>;
-    pcms::Rank2View<pcms::Real, pcms::DeviceMemorySpace, LayoutPolicy>
-      eval_results(eval_results_1d.data(), num_eval_points, 1);
-
-    evaluator->Evaluate(psi_field, eval_results);
+    Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace> eval_results(
+      "eval_results", num_eval_points, 1);
+    evaluator->Evaluate(psi_field, pcms::MakeRank2View(eval_results));
 
     auto results_host = Kokkos::create_mirror_view_and_copy(
-      pcms::HostMemorySpace(), eval_results_1d);
+      pcms::HostMemorySpace(), eval_results);
 
     // Verify that all results are finite
     for (int i = 0; i < num_eval_points; ++i) {
-      REQUIRE(std::isfinite(results_host(i)));
+      REQUIRE(std::isfinite(results_host(i, 0)));
     }
   }
 }
diff --git a/test/test_field_evaluation.cpp b/test/test_field_evaluation.cpp
index ff238e11..85af38b4 100644
--- a/test/test_field_evaluation.cpp
+++ b/test/test_field_evaluation.cpp
@@ -34,10 +34,10 @@ TEST_CASE("evaluate linear 2d omega_h_field")
 
   pcms::test::SetField(
     field.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   pcms::test::CheckEvaluation(
     factory, field, pcms::test::StandardEvalCoords2D(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 }
 
 #ifdef PCMS_ENABLE_MESHFIELDS
@@ -50,35 +50,11 @@ TEST_CASE("evaluate quadratic 2d meshfields_field")
     mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::MeshFields);
 
-  // Quadratic DOF holders span vertices and edge midpoints; set them inline.
-  const auto nverts = mesh.nents(0);
-  const auto nedges = mesh.nents(1);
-  auto mesh_coords = mesh.coords();
-  auto edge_verts = mesh.ask_verts_of(1);
-
-  Omega_h::Write<Real> test_f(nverts + nedges);
-  Omega_h::parallel_for(
-    nverts, OMEGA_H_LAMBDA(int i) {
-      test_f[i] = sin_f(mesh_coords[2 * static_cast<size_t>(i)],
-                        mesh_coords[2 * static_cast<size_t>(i) + 1]);
-    });
-  Omega_h::parallel_for(
-    nedges, OMEGA_H_LAMBDA(int i) {
-      auto ep = Omega_h::gather_verts<2>(edge_verts, i);
-      Real cx = (mesh_coords[2 * static_cast<size_t>(ep[0])] +
-                 mesh_coords[2 * static_cast<size_t>(ep[1])]) /
-                2;
-      Real cy = (mesh_coords[2 * static_cast<size_t>(ep[0]) + 1] +
-                 mesh_coords[2 * static_cast<size_t>(ep[1]) + 1]) /
-                2;
-      test_f[nverts + i] = sin_f(cx, cy);
-    });
-
-  Omega_h::HostWrite<Real> test_f_host(test_f);
+  // Quadratic DOF holders span vertices and edge midpoints; the layout's DOF
+  // coordinates cover both, so SetField samples sin_f at every holder.
   auto field = factory->CreateFunction<Real>();
-  field.GetData().SetDOFHolderDataHost(
-    pcms::Rank2View<const Real, pcms::HostMemorySpace>(test_f_host.data(),
-                                                       test_f_host.size(), 1));
+  pcms::test::SetField(
+    field, OMEGA_H_LAMBDA(Real x, Real y) { return sin_f(x, y); });
 
   pcms::test::CheckEvaluation(
     factory, field, kEvalCoords,
diff --git a/test/test_field_exchange_planner.cpp b/test/test_field_exchange_planner.cpp
index 61b5bc54..99e6181b 100644
--- a/test/test_field_exchange_planner.cpp
+++ b/test/test_field_exchange_planner.cpp
@@ -51,9 +51,7 @@ TEST_CASE("GID messages insert headers around compact field payloads",
     "gid_message",
     plan.msg_size + 2 * static_cast<std::size_t>(pcms::ent_offsets_len));
   pcms::GenericFieldExchangePlanner planner;
-  planner.FillGidMessage(layout, plan,
-                         pcms::Rank1View<pcms::GO, pcms::HostMemorySpace>(
-                           message.data(), message.size()));
+  planner.FillGidMessage(layout, plan, pcms::MakeRank1View(message));
 
   const pcms::GO expected[] = {
     0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 1, 3,
diff --git a/test/test_field_interpolation.cpp b/test/test_field_interpolation.cpp
index bc88ab97..f8e2432d 100644
--- a/test/test_field_interpolation.cpp
+++ b/test/test_field_interpolation.cpp
@@ -55,35 +55,11 @@ TEST_CASE("interpolate quadratic 2d meshfields_field")
   auto factory2 = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::MeshFields);
-  auto layout = factory2->GetLayout();
-  const auto nverts = mesh.nents(0);
-  const auto nedges = mesh.nents(1);
-  auto mesh_coords = mesh.coords();
-  auto edge_verts = mesh.ask_verts_of(1);
-  Omega_h::Write<Real> test_f(nverts + nedges);
-  Omega_h::parallel_for(
-    nverts, OMEGA_H_LAMBDA(int i) {
-      Real x = mesh_coords[2 * i + 0];
-      Real y = mesh_coords[2 * i + 1];
-      test_f[i] = interpolation_linear_f(x, y);
-    });
-  Omega_h::parallel_for(
-    nedges, OMEGA_H_LAMBDA(int i) {
-      auto endpoints = Omega_h::gather_verts<2>(edge_verts, i);
-      Real x0 = mesh_coords[2 * endpoints[0] + 0];
-      Real y0 = mesh_coords[2 * endpoints[0] + 1];
-      Real x1 = mesh_coords[2 * endpoints[1] + 0];
-      Real y1 = mesh_coords[2 * endpoints[1] + 1];
-      Real cx = (x0 + x1) / 2;
-      Real cy = (y0 + y1) / 2;
-      test_f[nverts + i] = interpolation_linear_f(cx, cy);
-    });
-
-  Omega_h::HostWrite<Real> test_f_host(test_f);
   auto field = factory2->CreateFunction<Real>();
   auto interpolated = factory2->CreateFunction<Real>();
-  field.SetDOFHolderDataHost(pcms::Rank2View<const Real, pcms::HostMemorySpace>(
-    test_f_host.data(), test_f_host.size(), 1));
+  pcms::test::SetField(
+    field,
+    OMEGA_H_LAMBDA(Real x, Real y) { return interpolation_linear_f(x, y); });
 
   pcms::Interpolator<Real> interp(*factory2, *factory2);
   interp.Apply(field, interpolated);
diff --git a/test/test_interpolation_class.cpp b/test/test_interpolation_class.cpp
index 7a34af4d..878f6288 100644
--- a/test/test_interpolation_class.cpp
+++ b/test/test_interpolation_class.cpp
@@ -10,36 +10,10 @@
 #include <pcms/transfer/interpolation_base.h>
 #include <pcms/utility/mesh_geometry.h>
 #include <pcms/utility/print.h>
+#include "field_test_utils.h"
 
 #include <vector>
 #include <iostream>
-#include <unordered_map>
-
-bool areArraysEqualUnordered(const Omega_h::HostRead<Omega_h::LO>& array1,
-                             const Omega_h::HostRead<Omega_h::LO>& array2,
-                             int start, int end)
-{
-  // Ensure the indices are valid
-  assert(start >= 0 && end <= array1.size() && start <= end);
-  assert(start >= 0 && end <= array2.size() && start <= end);
-
-  // Use frequency maps to count occurrences of each value
-  std::unordered_map<Omega_h::LO, int> freq1, freq2;
-
-  for (int i = start; i < end; ++i) {
-    freq1[array1[i]]++;
-    freq2[array2[i]]++;
-  }
-
-  // Compare the frequency maps
-  if (freq1 != freq2) {
-    pcms::printError("[ERROR] Arrays differ in the range [%d, %d)\n", start,
-                     end);
-    return false;
-  }
-
-  return true;
-}
 
 void translate_mesh(Omega_h::Mesh* mesh, Omega_h::Vector<2> translation_vector)
 {
@@ -135,28 +109,13 @@ TEST_CASE("Test MLSMeshInterpolation")
     auto mls_single =
       pcms::MLSMeshInterpolation(source_mesh, 0.12, 15, 3, true, 0.0, 5.0);
 
-    auto source_points_reals =
-      pcms::get_entity_centroids(source_mesh, Omega_h::FACE);
-    auto source_points_host =
-      Omega_h::HostRead<Omega_h::Real>(source_points_reals);
-    auto source_points_host_write =
-      Omega_h::HostWrite<Omega_h::Real>(source_points_host.size());
-    for (int i = 0; i < source_points_host.size(); i++) {
-      source_points_host_write[i] = source_points_host[i];
-    }
-    auto source_points_view = pcms::Rank1View<double, pcms::HostMemorySpace>(
-      source_points_host_write.data(), source_points_host_write.size());
-
-    auto target_points_reals = source_mesh.coords();
-    auto target_points_host =
-      Omega_h::HostRead<Omega_h::Real>(target_points_reals);
-    auto target_points_host_write =
-      Omega_h::HostWrite<Omega_h::Real>(target_points_host.size());
-    for (int i = 0; i < target_points_host.size(); i++) {
-      target_points_host_write[i] = target_points_host[i];
-    }
-    auto target_points_view = pcms::Rank1View<double, pcms::HostMemorySpace>(
-      target_points_host_write.data(), target_points_host_write.size());
+    auto source_points_vec = pcms::test::CopyOmegaHRealsToVector(
+      pcms::get_entity_centroids(source_mesh, Omega_h::FACE));
+    auto source_points_view = pcms::make_array_view(source_points_vec);
+
+    auto target_points_vec =
+      pcms::test::CopyOmegaHRealsToVector(source_mesh.coords());
+    auto target_points_view = pcms::make_array_view(target_points_vec);
     REQUIRE(source_mesh.dim() == 2);
     pcms::printInfo("Point cloud based search...\n");
     auto point_mls = pcms::MLSPointCloudInterpolation(
@@ -173,14 +132,11 @@ TEST_CASE("Test MLSMeshInterpolation")
       source_mesh.nverts());
     Omega_h::HostWrite<double> exact_values_at_nodes(source_sinxcosy_node);
 
-    pcms::Rank1View<double, pcms::HostMemorySpace> sourceArrayView(
-      source_data_host_write.data(), source_data_host_write.size());
-    pcms::Rank1View<double, pcms::HostMemorySpace> interpolatedArrayView(
-      interpolated_data_hwrite.data(), interpolated_data_hwrite.size());
-    pcms::Rank1View<double, pcms::HostMemorySpace>
-      point_cloud_interpolatedArrayView(
-        point_cloud_interpolated_data_hwrite.data(),
-        point_cloud_interpolated_data_hwrite.size());
+    auto sourceArrayView = pcms::make_array_view(source_data_host_write);
+    auto interpolatedArrayView =
+      pcms::make_array_view(interpolated_data_hwrite);
+    auto point_cloud_interpolatedArrayView =
+      pcms::make_array_view(point_cloud_interpolated_data_hwrite);
 
     OMEGA_H_CHECK_PRINTF(sourceArrayView.size() == mls_single.getSourceSize(),
                          "Source size mismatch: %zu vs %zu\n",
@@ -222,33 +178,8 @@ TEST_CASE("Test MLSMeshInterpolation")
     ///*****************************//
     auto mesh_based_supports = mls_single.getSupports();
     auto point_cloud_based_supports = point_mls.getSupports();
-    auto mesh_based_support_ptr_host =
-      Omega_h::HostRead<Omega_h::LO>(mesh_based_supports.supports_ptr);
-    auto mesh_based_support_idx_host =
-      Omega_h::HostRead<Omega_h::LO>(mesh_based_supports.supports_idx);
-    auto point_cloud_based_support_ptr_host =
-      Omega_h::HostRead<Omega_h::LO>(point_cloud_based_supports.supports_ptr);
-    auto point_cloud_based_support_idx_host =
-      Omega_h::HostRead<Omega_h::LO>(point_cloud_based_supports.supports_idx);
-
-    REQUIRE(point_cloud_based_support_idx_host.size() ==
-            mesh_based_support_idx_host.size());
-    REQUIRE(point_cloud_based_support_ptr_host.size() ==
-            mesh_based_support_ptr_host.size());
-    for (int i = 0; i < mesh_based_support_ptr_host.size(); i++) {
-      REQUIRE(point_cloud_based_support_ptr_host[i] ==
-              mesh_based_support_ptr_host[i]);
-    }
-
-    for (int i = 0; i < mesh_based_support_ptr_host.size() - 1; i++) {
-      auto start = mesh_based_support_ptr_host[i];
-      auto end = mesh_based_support_ptr_host[i + 1];
-
-      bool isEqual =
-        areArraysEqualUnordered(mesh_based_support_idx_host,
-                                point_cloud_based_support_idx_host, start, end);
-      REQUIRE(isEqual);
-    }
+    pcms::test::CheckSupportResultsEquivalent(point_cloud_based_supports,
+                                              mesh_based_supports);
 
     // Check if the point cloud interpolation is same as the MLS interpolation
     pcms::printDebugInfo("Interpolated data size: %d\n",
@@ -289,10 +220,9 @@ TEST_CASE("Test MLSMeshInterpolation")
     Omega_h::HostWrite<double> interpolated_data_hwrite(
       mls_double.getTargetSize());
 
-    pcms::Rank1View<double, pcms::HostMemorySpace> sourceArrayView(
-      source_data_host_write.data(), source_data_host_write.size());
-    pcms::Rank1View<double, pcms::HostMemorySpace> interpolatedArrayView(
-      interpolated_data_hwrite.data(), interpolated_data_hwrite.size());
+    auto sourceArrayView = pcms::make_array_view(source_data_host_write);
+    auto interpolatedArrayView =
+      pcms::make_array_view(interpolated_data_hwrite);
 
     mls_double.eval(sourceArrayView, interpolatedArrayView);
 
@@ -326,10 +256,8 @@ TEST_CASE("MLSPointCloudInterpolation honors provided dimension in eval")
   }
   auto target_points = source_points;
 
-  auto source_points_view = pcms::Rank1View<double, pcms::HostMemorySpace>(
-    source_points.data(), source_points.size());
-  auto target_points_view = pcms::Rank1View<double, pcms::HostMemorySpace>(
-    target_points.data(), target_points.size());
+  auto source_points_view = pcms::make_array_view(source_points);
+  auto target_points_view = pcms::make_array_view(target_points);
 
   // Degree-1 polynomial that depends on z to catch accidental 2D behavior.
   auto source_values = Omega_h::HostWrite<Omega_h::Real>(27);
@@ -341,10 +269,8 @@ TEST_CASE("MLSPointCloudInterpolation honors provided dimension in eval")
   }
 
   auto output_values = Omega_h::HostWrite<Omega_h::Real>(27, "output_values");
-  auto source_values_view = pcms::Rank1View<double, pcms::HostMemorySpace>(
-    source_values.data(), source_values.size());
-  auto output_values_view = pcms::Rank1View<double, pcms::HostMemorySpace>(
-    output_values.data(), output_values.size());
+  auto source_values_view = pcms::make_array_view(source_values);
+  auto output_values_view = pcms::make_array_view(output_values);
 
   auto mls = pcms::MLSPointCloudInterpolation(
     source_points_view, target_points_view, 3, 2.5, 10, 1, true, 0.0, 5.0);
diff --git a/test/test_interpolation_on_ltx_mesh.cpp b/test/test_interpolation_on_ltx_mesh.cpp
index a348acae..141ebfa8 100644
--- a/test/test_interpolation_on_ltx_mesh.cpp
+++ b/test/test_interpolation_on_ltx_mesh.cpp
@@ -83,16 +83,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]")
   printf("[INFO] Degas2 Mesh loaded from %s with %d elements\n",
          degas2_mesh_filename.c_str(), degas2_num_elems);
   const auto degas2_mesh_centroids_view =
-    pcms::Rank1View<const double, pcms::HostMemorySpace>(
-      degas2_mesh_centroids_host.data(), degas2_mesh_centroids_host.size());
+    pcms::make_array_view(degas2_mesh_centroids_host);
 
   auto xgc_mesh_points = read_xgc_mesh_nodes(ltx_mesh_base_filename + ".node");
   const int xgc_num_nodes = xgc_mesh_points.size() / 2;
   printf("[INFO] XGC Mesh loaded from %s with %d points\n",
          ltx_mesh_base_filename.c_str(), xgc_num_nodes);
-  const auto xgc_mesh_points_view =
-    pcms::Rank1View<const double, pcms::HostMemorySpace>(
-      xgc_mesh_points.data(), xgc_mesh_points.size());
+  const auto xgc_mesh_points_view = pcms::make_array_view(xgc_mesh_points);
 
   auto xgc_to_degas2_interpolator = pcms::MLSPointCloudInterpolation(
     xgc_mesh_points_view, degas2_mesh_centroids_view, 2, 0.000001, 10, 1, true,
@@ -120,18 +117,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]")
 
   // ------------------ First Interpolation ------------------ //
   const auto density_at_xgc_nodes_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(density_at_xgc_nodes.data(),
-                                                   density_at_xgc_nodes.size());
-  const auto temp_at_xgc_nodes_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(temp_at_xgc_nodes.data(),
-                                                   temp_at_xgc_nodes.size());
+    pcms::make_array_view(density_at_xgc_nodes);
+  const auto temp_at_xgc_nodes_view = pcms::make_array_view(temp_at_xgc_nodes);
 
   const auto density_at_degas2_centroids_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      density_at_degas2_centroids.data(), density_at_degas2_centroids.size());
+    pcms::make_array_view(density_at_degas2_centroids);
   const auto temp_at_degas2_centroids_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      temp_at_degas2_centroids.data(), temp_at_degas2_centroids.size());
+    pcms::make_array_view(temp_at_degas2_centroids);
 
   Omega_h::HostWrite<Omega_h::Real> interpolated_xgc_density(degas2_num_elems);
   Omega_h::HostWrite<Omega_h::Real> interpolated_xgc_temp(degas2_num_elems);
@@ -139,17 +131,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]")
   Omega_h::HostWrite<Omega_h::Real> interpolated_degas2_temp(xgc_num_nodes);
 
   const auto interpolated_xgc_density_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      interpolated_xgc_density.data(), interpolated_xgc_density.size());
+    pcms::make_array_view(interpolated_xgc_density);
   const auto interpolated_xgc_temp_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      interpolated_xgc_temp.data(), interpolated_xgc_temp.size());
+    pcms::make_array_view(interpolated_xgc_temp);
   const auto interpolated_degas2_density_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      interpolated_degas2_density.data(), interpolated_degas2_density.size());
+    pcms::make_array_view(interpolated_degas2_density);
   const auto interpolated_degas2_temp_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      interpolated_degas2_temp.data(), interpolated_degas2_temp.size());
+    pcms::make_array_view(interpolated_degas2_temp);
 
   xgc_to_degas2_interpolator.eval(density_at_xgc_nodes_view,
                                   interpolated_xgc_density_view);
@@ -174,21 +162,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]")
     degas2_num_elems);
 
   const auto interpolated_back_density_at_xgc_nodes_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      interpolated_back_density_at_xgc_nodes.data(),
-      interpolated_back_density_at_xgc_nodes.size());
+    pcms::make_array_view(interpolated_back_density_at_xgc_nodes);
   const auto interpolated_back_temp_at_xgc_nodes_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      interpolated_back_temp_at_xgc_nodes.data(),
-      interpolated_back_temp_at_xgc_nodes.size());
+    pcms::make_array_view(interpolated_back_temp_at_xgc_nodes);
   const auto interpolated_back_density_at_degas2_centroids_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      interpolated_back_density_at_degas2_centroids.data(),
-      interpolated_back_density_at_degas2_centroids.size());
+    pcms::make_array_view(interpolated_back_density_at_degas2_centroids);
   const auto interpolated_back_temp_at_degas2_centroids_view =
-    pcms::Rank1View<double, pcms::HostMemorySpace>(
-      interpolated_back_temp_at_degas2_centroids.data(),
-      interpolated_back_temp_at_degas2_centroids.size());
+    pcms::make_array_view(interpolated_back_temp_at_degas2_centroids);
 
   degas2_to_xgc_interpolator.eval(interpolated_xgc_density_view,
                                   interpolated_back_density_at_xgc_nodes_view);
diff --git a/test/test_localization_factory.cpp b/test/test_localization_factory.cpp
index 571ac175..d65e8fc0 100644
--- a/test/test_localization_factory.cpp
+++ b/test/test_localization_factory.cpp
@@ -14,6 +14,7 @@
 #include "pcms/localization/point_cloud_localization.h"
 #include "pcms/discretization/discretization/omega_h.hpp"
 #include "pcms/utility/arrays.h"
+#include "pcms/utility/omega_h_array_utils.h"
 #include "field_test_utils.h"
 
 namespace
@@ -46,17 +47,7 @@ TEST_CASE(
   auto target_device = pcms::test::CreateDeviceCoordinateView(
     target_coords, pcms::CoordinateSystem::Cartesian, dim);
 
-  auto coords_read = Omega_h::HostRead<Omega_h::Real>(source_coords);
-
-  // auto coords_dev = Kokkos::create_mirror_view_and_copy(
-  //   Kokkos::DefaultExecutionSpace{}, coords_host);
-  auto coords_dev = Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace>(
-    "point_cloud_coords", mesh.nverts(), dim);
-  auto coords_host = Kokkos::create_mirror(pcms::HostMemorySpace(), coords_dev);
-  for (int i = 0; i < mesh.nverts(); ++i)
-    for (int d = 0; d < dim; ++d)
-      coords_host(i, d) = coords_read[i * dim + d];
-  Kokkos::deep_copy(coords_dev, coords_host);
+  auto coords_dev = pcms::ConvertCoordsTo2D(source_coords, mesh.nverts(), dim);
 
   auto layout = std::make_shared<pcms::PointCloudLayout>(
     dim, coords_dev, pcms::CoordinateSystem::Cartesian);
diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp
index fbb34ddb..a627ab81 100644
--- a/test/test_mesh_intersection_field_transfer.cpp
+++ b/test/test_mesh_intersection_field_transfer.cpp
@@ -7,53 +7,14 @@
 
 #include <pcms/transfer/omega_h_conservative_projection.hpp>
 #include <pcms/field/function_space/lagrange.h>
+#include <pcms/utility/mesh_geometry.h>
 #include "field_test_utils.h"
 
-#include <array>
 #include <vector>
 
 namespace
 {
 
-// Builds a unit-square 2D simplex mesh from the given element connectivity and
-// adds the geometric classification tags required to build an Omega_h-backed
-// Lagrange function space.
-Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, const Omega_h::LOs& ev2v)
-{
-  const Omega_h::Reals coords({
-    0.0, 0.0, // v0
-    1.0, 0.0, // v1
-    1.0, 1.0, // v2
-    0.0, 1.0  // v3
-  });
-  Omega_h::Mesh mesh(&lib);
-  Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords);
-  for (Omega_h::Int dim = 0; dim <= 2; ++dim) {
-    mesh.add_tag<Omega_h::I8>(
-      dim, "class_dim", 1,
-      Omega_h::Read<Omega_h::I8>(mesh.nents(dim), Omega_h::I8(dim)));
-    mesh.add_tag<Omega_h::ClassId>(
-      dim, "class_id", 1,
-      Omega_h::Read<Omega_h::ClassId>(mesh.nents(dim), Omega_h::ClassId(0)));
-  }
-  return mesh;
-}
-
-std::shared_ptr<pcms::LagrangeFunctionSpace> MakeP1Space(
-  Omega_h::Mesh& mesh, const std::string& global_id_name = "global")
-{
-  return pcms::LagrangeFunctionSpace::FromMesh(
-    mesh, 1, 1, pcms::CoordinateSystem::Cartesian, global_id_name,
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
-}
-
-std::shared_ptr<pcms::LagrangeFunctionSpace> MakeP0Space(Omega_h::Mesh& mesh)
-{
-  return pcms::LagrangeFunctionSpace::FromMesh(
-    mesh, 0, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
-}
-
 void AddReorderedVertexGlobalIds(Omega_h::Mesh& mesh)
 {
   mesh.add_tag<Omega_h::GO>(
@@ -68,64 +29,6 @@ void AddSparseVertexGlobalIds(Omega_h::Mesh& mesh)
     Omega_h::Read<Omega_h::GO>({102, 7, 41, 19}, "sparse_global"));
 }
 
-// Integral over the mesh of an order-0 (piecewise-constant) field whose values
-// are stored in element order. Exact by construction.
-double IntegrateP0Field(Omega_h::Mesh& mesh,
-                        const pcms::Field<pcms::Real>& field)
-{
-  const auto values = pcms::FlattenToRank1View(field.GetDOFHolderDataHost());
-  const auto elem_areas = Omega_h::measure_elements_real(&mesh);
-  const auto elem_areas_h = Omega_h::HostRead<Omega_h::Real>(elem_areas);
-
-  double integral = 0.0;
-  for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) {
-    integral += elem_areas_h[e] * values[e];
-  }
-  return integral;
-}
-
-// Per-element centroid coordinates in element order.
-std::vector<std::array<double, 2>> ElementCentroids(Omega_h::Mesh& mesh)
-{
-  const auto coords_h = Omega_h::HostRead<Omega_h::Real>(mesh.coords());
-  const auto elem_verts_h =
-    Omega_h::HostRead<Omega_h::LO>(mesh.ask_elem_verts());
-  std::vector<std::array<double, 2>> centroids(mesh.nelems());
-  for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) {
-    double cx = 0.0, cy = 0.0;
-    for (int k = 0; k < 3; ++k) {
-      const Omega_h::LO v = elem_verts_h[3 * e + k];
-      cx += coords_h[2 * v + 0];
-      cy += coords_h[2 * v + 1];
-    }
-    centroids[e] = {cx / 3.0, cy / 3.0};
-  }
-  return centroids;
-}
-
-// Integral over the mesh of an order-1 nodal field whose values are stored in
-// vertex order (as produced by Field<Real>::GetDOFHolderDataHost for the
-// Omega_h backend), computed exactly via per-element vertex averaging.
-double IntegrateP1Field(Omega_h::Mesh& mesh,
-                        const pcms::Field<pcms::Real>& field)
-{
-  const auto values = pcms::FlattenToRank1View(field.GetDOFHolderDataHost());
-  const auto elem_areas = Omega_h::measure_elements_real(&mesh);
-  const auto elem_verts = mesh.ask_elem_verts();
-  const auto elem_areas_h = Omega_h::HostRead<Omega_h::Real>(elem_areas);
-  const auto elem_verts_h = Omega_h::HostRead<Omega_h::LO>(elem_verts);
-
-  double integral = 0.0;
-  for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) {
-    const Omega_h::LO v0 = elem_verts_h[3 * e + 0];
-    const Omega_h::LO v1 = elem_verts_h[3 * e + 1];
-    const Omega_h::LO v2 = elem_verts_h[3 * e + 2];
-    const double avg = (values[v0] + values[v1] + values[v2]) / 3.0;
-    integral += elem_areas_h[e] * avg;
-  }
-  return integral;
-}
-
 } // namespace
 
 // Fields that already live in the target order-1 space (constants and affine
@@ -139,12 +42,12 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields",
 
   // Source and target triangulate the same square along opposite diagonals.
   Omega_h::Mesh source_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
   Omega_h::Mesh target_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
 
-  auto source_space = MakeP1Space(source_mesh);
-  auto target_space = MakeP1Space(target_mesh);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   auto source = source_space->CreateFunction<pcms::Real>();
   auto target = target_space->CreateFunction<pcms::Real>();
@@ -166,8 +69,9 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields",
     for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) {
       REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-10));
     }
-    REQUIRE(IntegrateP1Field(target_mesh, target) ==
-            Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-10));
+    REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) ==
+            Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source))
+              .margin(1e-10));
   }
 
   SECTION("linear field is reproduced on target vertices and conserved")
@@ -179,34 +83,31 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields",
 
     const auto target_values =
       pcms::FlattenToRank1View(target.GetDOFHolderDataHost());
-    const auto tgt_coords_h =
-      Omega_h::HostRead<Omega_h::Real>(target_mesh.coords());
+    const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost(
+      pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(),
+      2);
     for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) {
-      const double expected = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1];
+      const double expected = tgt_coords_h(i, 0) + tgt_coords_h(i, 1);
       REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-9));
     }
-    REQUIRE(IntegrateP1Field(target_mesh, target) ==
-            Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9));
+    REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) ==
+            Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source))
+              .margin(1e-9));
   }
 }
 
-// For a field that does not live in the target space the projection is not an
-// interpolation, but conservation of the integral is the defining property of
-// the conservative (Galerkin) projection and must still hold exactly. This
-// also exercises a second Apply with a different source field to confirm the
-// cached integrators/evaluator/factorization are reused without stale state.
 TEST_CASE("OmegaHConservativeProjection conserves the integral",
           "[transfer][mesh_intersection]")
 {
   Omega_h::Library lib;
 
   Omega_h::Mesh source_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
   Omega_h::Mesh target_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
 
-  auto source_space = MakeP1Space(source_mesh);
-  auto target_space = MakeP1Space(target_mesh);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   auto source = source_space->CreateFunction<pcms::Real>();
   auto target = target_space->CreateFunction<pcms::Real>();
@@ -218,8 +119,9 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral",
       return x * x + x * y + 0.5 * y * y;
     });
   projection.Apply(source, target);
-  REQUIRE(IntegrateP1Field(target_mesh, target) ==
-          Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12));
+  REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) ==
+          Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source))
+            .margin(1e-12));
 
   // Second Apply with a different source field — verifies cached state is
   // reused correctly and is not stale.
@@ -227,8 +129,9 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral",
     source,
     OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return 2.0 * x - y + 0.5; });
   projection.Apply(source, target);
-  REQUIRE(IntegrateP1Field(target_mesh, target) ==
-          Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12));
+  REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) ==
+          Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source))
+            .margin(1e-12));
 }
 
 TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local "
@@ -238,13 +141,13 @@ TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local "
   Omega_h::Library lib;
 
   Omega_h::Mesh source_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
   Omega_h::Mesh target_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
   AddReorderedVertexGlobalIds(target_mesh);
 
-  auto source_space = MakeP1Space(source_mesh);
-  auto target_space = MakeP1Space(target_mesh, "reordered_global");
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh, "reordered_global");
 
   auto source = source_space->CreateFunction<pcms::Real>();
   auto target = target_space->CreateFunction<pcms::Real>();
@@ -277,13 +180,13 @@ TEST_CASE("OmegaHConservativeProjection maps sparse target GIDs to active "
   Omega_h::Library lib;
 
   Omega_h::Mesh source_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
   Omega_h::Mesh target_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
   AddSparseVertexGlobalIds(target_mesh);
 
-  auto source_space = MakeP1Space(source_mesh);
-  auto target_space = MakeP1Space(target_mesh, "sparse_global");
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh, "sparse_global");
 
   auto source = source_space->CreateFunction<pcms::Real>();
   auto target = target_space->CreateFunction<pcms::Real>();
@@ -318,12 +221,12 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target",
   Omega_h::Library lib;
 
   Omega_h::Mesh source_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
   Omega_h::Mesh target_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
 
-  auto source_space = MakeP0Space(source_mesh);
-  auto target_space = MakeP1Space(target_mesh);
+  auto source_space = pcms::test::MakeP0Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   auto source = source_space->CreateFunction<pcms::Real>();
   auto target = target_space->CreateFunction<pcms::Real>();
@@ -345,8 +248,9 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target",
     for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) {
       REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-10));
     }
-    REQUIRE(IntegrateP1Field(target_mesh, target) ==
-            Catch::Approx(IntegrateP0Field(source_mesh, source)).margin(1e-10));
+    REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) ==
+            Catch::Approx(pcms::test::IntegrateP0Field(source_mesh, source))
+              .margin(1e-10));
   }
 
   SECTION("non-constant P0 source conserves the integral")
@@ -357,8 +261,9 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target",
 
     projection.Apply(source, target);
 
-    REQUIRE(IntegrateP1Field(target_mesh, target) ==
-            Catch::Approx(IntegrateP0Field(source_mesh, source)).margin(1e-10));
+    REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) ==
+            Catch::Approx(pcms::test::IntegrateP0Field(source_mesh, source))
+              .margin(1e-10));
   }
 }
 
@@ -372,12 +277,12 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target",
   Omega_h::Library lib;
 
   Omega_h::Mesh source_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3}));
   Omega_h::Mesh target_mesh =
-    BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
+    pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3}));
 
-  auto source_space = MakeP1Space(source_mesh);
-  auto target_space = MakeP0Space(target_mesh);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP0Space(target_mesh);
 
   auto source = source_space->CreateFunction<pcms::Real>();
   auto target = target_space->CreateFunction<pcms::Real>();
@@ -399,8 +304,9 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target",
     for (Omega_h::LO e = 0; e < target_mesh.nelems(); ++e) {
       REQUIRE(target_values[e] == Catch::Approx(c).margin(1e-10));
     }
-    REQUIRE(IntegrateP0Field(target_mesh, target) ==
-            Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-10));
+    REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) ==
+            Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source))
+              .margin(1e-10));
   }
 
   SECTION("linear field projects to cell average and conserves the integral")
@@ -414,13 +320,16 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target",
 
     const auto target_values =
       pcms::FlattenToRank1View(target.GetDOFHolderDataHost());
-    const auto centroids = ElementCentroids(target_mesh);
-    REQUIRE(static_cast<std::size_t>(target_values.size()) == centroids.size());
-    for (std::size_t e = 0; e < centroids.size(); ++e) {
-      const double expected = f(centroids[e][0], centroids[e][1]);
+    const auto centroids_h = Omega_h::HostRead<Omega_h::Real>(
+      pcms::get_entity_centroids(target_mesh, target_mesh.dim()));
+    REQUIRE(static_cast<Omega_h::LO>(target_values.size()) ==
+            target_mesh.nelems());
+    for (Omega_h::LO e = 0; e < target_mesh.nelems(); ++e) {
+      const double expected = f(centroids_h[2 * e + 0], centroids_h[2 * e + 1]);
       REQUIRE(target_values[e] == Catch::Approx(expected).margin(1e-9));
     }
-    REQUIRE(IntegrateP0Field(target_mesh, target) ==
-            Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9));
+    REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) ==
+            Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source))
+              .margin(1e-9));
   }
 }
diff --git a/test/test_omega_h_intersection_rhs_integrator.cpp b/test/test_omega_h_intersection_rhs_integrator.cpp
index 18ea11bb..f5d667d9 100644
--- a/test/test_omega_h_intersection_rhs_integrator.cpp
+++ b/test/test_omega_h_intersection_rhs_integrator.cpp
@@ -21,37 +21,6 @@
 namespace
 {
 
-// Build a unit-square 2D simplex mesh.
-//   diagonal=0: T0=(0,1,3), T1=(1,2,3)
-//   diagonal=1: T0=(0,1,2), T1=(0,2,3)
-Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal)
-{
-  const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0});
-  Omega_h::LOs ev2v = (diagonal == 0) ? Omega_h::LOs({0, 1, 3, 1, 2, 3})
-                                      : Omega_h::LOs({0, 1, 2, 0, 2, 3});
-  Omega_h::Mesh mesh(&lib);
-  Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords);
-  mesh.add_tag<Omega_h::I8>(
-    0, "class_dim", 1,
-    Omega_h::Read<Omega_h::I8>(mesh.nverts(), Omega_h::I8(0)));
-  mesh.add_tag<Omega_h::ClassId>(
-    0, "class_id", 1,
-    Omega_h::Read<Omega_h::ClassId>(mesh.nverts(), Omega_h::ClassId(0)));
-  mesh.add_tag<Omega_h::I8>(
-    1, "class_dim", 1,
-    Omega_h::Read<Omega_h::I8>(mesh.nedges(), Omega_h::I8(1)));
-  mesh.add_tag<Omega_h::ClassId>(
-    1, "class_id", 1,
-    Omega_h::Read<Omega_h::ClassId>(mesh.nedges(), Omega_h::ClassId(0)));
-  mesh.add_tag<Omega_h::I8>(
-    2, "class_dim", 1,
-    Omega_h::Read<Omega_h::I8>(mesh.nelems(), Omega_h::I8(2)));
-  mesh.add_tag<Omega_h::ClassId>(
-    2, "class_id", 1,
-    Omega_h::Read<Omega_h::ClassId>(mesh.nelems(), Omega_h::ClassId(0)));
-  return mesh;
-}
-
 // Independent reference for the assembled conservative load vector
 //   b_j = \int phi_j^target f dx
 // when f is exactly representable in the target order-1 space (constant or
@@ -100,17 +69,7 @@ void CheckAssembledLoadMatches(
   const pcms::Field<pcms::Real>& source_field,
   const std::unordered_map<pcms::GO, pcms::Real>& expected)
 {
-  const auto pts = integrator.GetIntegrationPoints();
-  const std::size_t npts = pts.GetValues().extent(0);
-
-  auto evaluator = source_space->CreatePointEvaluator<pcms::Real>(
-    pcms::EvaluationRequest::FromCoordinates(pts));
-
-  Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace> sampled("sampled", npts,
-                                                              1);
-  evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled));
-
-  integrator.Assemble(pcms::MakeConstRank2View(sampled));
+  pcms::test::EvaluateAndAssemble(integrator, source_space, source_field);
   Vec vec = integrator.GetVector();
 
   const PetscScalar* vec_array = nullptr;
@@ -135,26 +94,19 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: integration points lie inside "
           "[rhs_integrator]")
 {
   Omega_h::Library lib;
-  auto source_mesh = BuildUnitSquare(lib, 1);
-  auto target_mesh = BuildUnitSquare(lib, 0);
+  auto source_mesh = pcms::test::BuildUnitSquare(lib, 1);
+  auto target_mesh = pcms::test::BuildUnitSquare(lib, 0);
 
-  auto source_space = pcms::LagrangeFunctionSpace::FromMesh(
-    source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto target_space = pcms::LagrangeFunctionSpace::FromMesh(
-    target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   auto integrator =
     pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space);
 
   const auto raw_coords = integrator->GetIntegrationPoints().GetValues();
-  auto raw_coords_view = Kokkos::View<const pcms::Real**, Kokkos::LayoutRight,
-                                      pcms::DeviceMemorySpace,
-                                      Kokkos::MemoryTraits<Kokkos::Unmanaged>>(
-    raw_coords.data_handle(), raw_coords.extent(0), raw_coords.extent(1));
-  auto raw_coords_host =
-    Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, raw_coords_view);
+  auto raw_coords_host = pcms::test::CopyCoordinatesToHost(
+    raw_coords, static_cast<int>(raw_coords.extent(0)),
+    static_cast<int>(raw_coords.extent(1)));
   const std::size_t n = raw_coords_host.extent(0);
 
   REQUIRE(n > 0);
@@ -172,32 +124,18 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: zero source field gives zero "
           "[rhs_integrator]")
 {
   Omega_h::Library lib;
-  auto source_mesh = BuildUnitSquare(lib, 1);
-  auto target_mesh = BuildUnitSquare(lib, 0);
+  auto source_mesh = pcms::test::BuildUnitSquare(lib, 1);
+  auto target_mesh = pcms::test::BuildUnitSquare(lib, 0);
 
-  auto source_space = pcms::LagrangeFunctionSpace::FromMesh(
-    source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto target_space = pcms::LagrangeFunctionSpace::FromMesh(
-    target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   auto source_field = source_space->CreateFunction<pcms::Real>();
   // Default-constructed field has zero data.
 
   auto integrator =
     pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space);
-  const auto& pts = integrator->GetIntegrationPoints();
-  const std::size_t npts = pts.GetValues().extent(0);
-
-  auto evaluator = source_space->CreatePointEvaluator<pcms::Real>(
-    pcms::EvaluationRequest::FromCoordinates(pts));
-
-  Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace> sampled("sampled", npts,
-                                                              1);
-  evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled));
-
-  integrator->Assemble(pcms::MakeConstRank2View(sampled));
+  pcms::test::EvaluateAndAssemble(*integrator, source_space, source_field);
   Vec vec = integrator->GetVector();
 
   PetscReal norm = 0.0;
@@ -214,15 +152,11 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: constant field matches "
   // For a constant source field the assembled load vector must equal the
   // target consistent mass matrix applied to the constant nodal vector.
   Omega_h::Library lib;
-  auto source_mesh = BuildUnitSquare(lib, 1);
-  auto target_mesh = BuildUnitSquare(lib, 0);
+  auto source_mesh = pcms::test::BuildUnitSquare(lib, 1);
+  auto target_mesh = pcms::test::BuildUnitSquare(lib, 0);
 
-  auto source_space = pcms::LagrangeFunctionSpace::FromMesh(
-    source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto target_space = pcms::LagrangeFunctionSpace::FromMesh(
-    target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   const double c = 2.0;
   auto source_field = source_space->CreateFunction<pcms::Real>();
@@ -251,15 +185,11 @@ TEST_CASE(
   // vector must equal the target consistent mass matrix applied to the nodal
   // values of x + y.
   Omega_h::Library lib;
-  auto source_mesh = BuildUnitSquare(lib, 1);
-  auto target_mesh = BuildUnitSquare(lib, 0);
+  auto source_mesh = pcms::test::BuildUnitSquare(lib, 1);
+  auto target_mesh = pcms::test::BuildUnitSquare(lib, 0);
 
-  auto source_space = pcms::LagrangeFunctionSpace::FromMesh(
-    source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto target_space = pcms::LagrangeFunctionSpace::FromMesh(
-    target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   auto source_field = source_space->CreateFunction<pcms::Real>();
   pcms::test::SetField(
@@ -268,11 +198,11 @@ TEST_CASE(
   const auto target_layout =
     std::dynamic_pointer_cast<const pcms::OmegaHLagrangeLayout>(
       target_space->GetLayout());
-  const auto tgt_coords_h =
-    Omega_h::HostRead<Omega_h::Real>(target_mesh.coords());
+  const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost(
+    pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), 2);
   std::vector<pcms::Real> g(target_mesh.nverts());
   for (int i = 0; i < target_mesh.nverts(); ++i) {
-    g[i] = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1];
+    g[i] = tgt_coords_h(i, 0) + tgt_coords_h(i, 1);
   }
   const auto expected = ExpectedLoadByGid(target_mesh, *target_layout, g);
 
@@ -285,26 +215,22 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts",
           "[rhs_integrator]")
 {
   Omega_h::Library lib;
-  auto source_mesh = BuildUnitSquare(lib, 1);
-  auto target_mesh = BuildUnitSquare(lib, 0);
+  auto source_mesh = pcms::test::BuildUnitSquare(lib, 1);
+  auto target_mesh = pcms::test::BuildUnitSquare(lib, 0);
 
   SECTION("multi-component source space throws")
   {
     auto source_space = pcms::LagrangeFunctionSpace::FromMesh(
       source_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
-    auto target_space = pcms::LagrangeFunctionSpace::FromMesh(
-      target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-      pcms::LagrangeFunctionSpace::Backend::OmegaH);
+    auto target_space = pcms::test::MakeP1Space(target_mesh);
     REQUIRE_THROWS(
       pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space));
   }
 
   SECTION("multi-component target space throws")
   {
-    auto source_space = pcms::LagrangeFunctionSpace::FromMesh(
-      source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-      pcms::LagrangeFunctionSpace::Backend::OmegaH);
+    auto source_space = pcms::test::MakeP1Space(source_mesh);
     auto target_space = pcms::LagrangeFunctionSpace::FromMesh(
       target_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
@@ -317,18 +243,14 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts",
     auto source_space = pcms::LagrangeFunctionSpace::FromMesh(
       source_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
-    auto target_space = pcms::LagrangeFunctionSpace::FromMesh(
-      target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-      pcms::LagrangeFunctionSpace::Backend::OmegaH);
+    auto target_space = pcms::test::MakeP1Space(target_mesh);
     REQUIRE_THROWS(
       pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space));
   }
 
   SECTION("non-Cartesian target coordinate system throws")
   {
-    auto source_space = pcms::LagrangeFunctionSpace::FromMesh(
-      source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-      pcms::LagrangeFunctionSpace::Backend::OmegaH);
+    auto source_space = pcms::test::MakeP1Space(source_mesh);
     auto target_space = pcms::LagrangeFunctionSpace::FromMesh(
       target_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
diff --git a/test/test_omega_h_lagrange_field.cpp b/test/test_omega_h_lagrange_field.cpp
index c17d45b1..0dd5ffc7 100644
--- a/test/test_omega_h_lagrange_field.cpp
+++ b/test/test_omega_h_lagrange_field.cpp
@@ -155,10 +155,10 @@ TEST_CASE("OmegaHLagrangeField order-1: linear function evaluation")
 
   pcms::test::SetField(
     field.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   pcms::test::CheckEvaluation(
     factory, field, pcms::test::StandardEvalCoords2D(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 }
 
 // The same linear evaluation test run on the MeshFields-backed order-1 field
@@ -173,10 +173,10 @@ TEST_CASE("MeshFieldsAdapter order-1: linear function evaluation (shared util)")
 
   pcms::test::SetField(
     field.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   pcms::test::CheckEvaluation(
     factory, field, pcms::test::StandardEvalCoords2D(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 }
 
 TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode")
@@ -189,10 +189,10 @@ TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode")
 
   pcms::test::SetField(
     field.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 
   Real fill_value = -999.0;
-  std::vector<Real> outside{-0.5, 0.5, 1.5, 0.5, 0.5, -0.5, 0.5, 1.5};
+  auto outside = pcms::test::StandardOutsideCoords2D();
   pcms::test::CheckFillMode(factory, field, fill_value, outside);
 }
 
@@ -206,7 +206,7 @@ TEST_CASE("OmegaHLagrangeField order-1: serialize / deserialize round-trip")
 
   pcms::test::SetField(
     field.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   pcms::test::CheckSerializeDeserialize(*factory->GetLayout(), field.GetData());
 }
 
@@ -270,24 +270,9 @@ TEST_CASE("OmegaHLagrangeField order-0: constant field evaluation")
                                                           1);
   field.GetData().SetDOFHolderDataHost(view);
 
-  auto pts = pcms::test::StandardEvalCoords2D();
-  int n = static_cast<int>(pts.size()) / 2;
-  auto device_coords =
-    pcms::test::CreateDeviceCoordinateView(pts, factory->GetCoordinateSystem());
-  auto evaluator = factory->CreatePointEvaluator<Real>(
-    pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
-
-  Kokkos::View<Real*, pcms::DeviceMemorySpace> eval_device("eval", n);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<pcms::DeviceMemorySpace>;
-  auto out = pcms::Rank2View<Real, pcms::DeviceMemorySpace, LayoutPolicy>(
-    eval_device.data(), n, 1);
-  evaluator->Evaluate(field, out);
-  auto eval_host =
-    Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), eval_device);
-
-  for (int i = 0; i < n; ++i)
-    REQUIRE(eval_host(i) == Catch::Approx(kValue));
+  pcms::test::CheckEvaluation(
+    factory, field, pcms::test::StandardEvalCoords2D(),
+    OMEGA_H_LAMBDA(Real, Real) { return kValue; });
 }
 
 TEST_CASE("OmegaHLagrangeField order-0: out-of-bounds FILL mode")
@@ -357,7 +342,8 @@ TEST_CASE("OmegaHLagrangeField: field valid after layout destruction")
   } // factory goes out of scope; field keeps layout alive
 
   pcms::test::SetField(
-    *field, OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    *field,
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   // Just verify data was set correctly (no evaluator needed for this lifetime
   // test)
   auto data = field->GetDOFHolderDataHost();
diff --git a/test/test_omega_h_mass_integrator.cpp b/test/test_omega_h_mass_integrator.cpp
index c9013f9d..c49a2e7d 100644
--- a/test/test_omega_h_mass_integrator.cpp
+++ b/test/test_omega_h_mass_integrator.cpp
@@ -8,6 +8,7 @@
 #include <pcms/field/layout/omega_h_lagrange.h>
 #include <pcms/transfer/mass_matrix_integrator.hpp>
 #include <pcms/transfer/omega_h_mass_integrator.hpp>
+#include "field_test_utils.h"
 #include <KokkosController.hpp>
 #include <MeshField.hpp>
 #include <MeshField_Element.hpp>
@@ -23,33 +24,6 @@
 namespace
 {
 
-Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib)
-{
-  const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0});
-  Omega_h::LOs ev2v({0, 1, 3, 1, 2, 3});
-  Omega_h::Mesh mesh(&lib);
-  Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords);
-  mesh.add_tag<Omega_h::I8>(
-    0, "class_dim", 1,
-    Omega_h::Read<Omega_h::I8>(mesh.nverts(), Omega_h::I8(0)));
-  mesh.add_tag<Omega_h::ClassId>(
-    0, "class_id", 1,
-    Omega_h::Read<Omega_h::ClassId>(mesh.nverts(), Omega_h::ClassId(0)));
-  mesh.add_tag<Omega_h::I8>(
-    1, "class_dim", 1,
-    Omega_h::Read<Omega_h::I8>(mesh.nedges(), Omega_h::I8(1)));
-  mesh.add_tag<Omega_h::ClassId>(
-    1, "class_id", 1,
-    Omega_h::Read<Omega_h::ClassId>(mesh.nedges(), Omega_h::ClassId(0)));
-  mesh.add_tag<Omega_h::I8>(
-    2, "class_dim", 1,
-    Omega_h::Read<Omega_h::I8>(mesh.nelems(), Omega_h::I8(2)));
-  mesh.add_tag<Omega_h::ClassId>(
-    2, "class_id", 1,
-    Omega_h::Read<Omega_h::ClassId>(mesh.nelems(), Omega_h::ClassId(0)));
-  return mesh;
-}
-
 std::map<std::pair<pcms::GO, pcms::GO>, pcms::Real> BuildReferenceMassMap(
   Omega_h::Mesh& mesh, const pcms::FunctionSpace& space)
 {
@@ -98,11 +72,9 @@ TEST_CASE(
   "[mass_integrator]")
 {
   Omega_h::Library lib;
-  auto mesh = BuildUnitSquare(lib);
+  auto mesh = pcms::test::BuildUnitSquare(lib, 0);
 
-  auto space = pcms::LagrangeFunctionSpace::FromMesh(
-    mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
+  auto space = pcms::test::MakeP1Space(mesh);
 
   const auto ref = BuildReferenceMassMap(mesh, *space);
 
@@ -128,11 +100,9 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass",
   // each node: row_sum(i) = integral(N_i dx). For a regular mesh of unit
   // area with uniform nodal distribution the sum over all nodes equals 1.
   Omega_h::Library lib;
-  auto mesh = BuildUnitSquare(lib);
+  auto mesh = pcms::test::BuildUnitSquare(lib, 0);
 
-  auto space = pcms::LagrangeFunctionSpace::FromMesh(
-    mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
+  auto space = pcms::test::MakeP1Space(mesh);
 
   const auto layout =
     std::dynamic_pointer_cast<const pcms::OmegaHLagrangeLayout>(
@@ -169,7 +139,7 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass",
 TEST_CASE("OmegaHMassIntegrator: rejects invalid layouts", "[mass_integrator]")
 {
   Omega_h::Library lib;
-  auto mesh = BuildUnitSquare(lib);
+  auto mesh = pcms::test::BuildUnitSquare(lib, 0);
 
   SECTION("multi-component space throws")
   {
diff --git a/test/test_omega_h_mc_rhs_integrator.cpp b/test/test_omega_h_mc_rhs_integrator.cpp
index 43172fce..4cf8a651 100644
--- a/test/test_omega_h_mc_rhs_integrator.cpp
+++ b/test/test_omega_h_mc_rhs_integrator.cpp
@@ -15,69 +15,6 @@
 #include <petscksp.h>
 #include <vector>
 
-// ---------------------------------------------------------------------------
-// Helpers
-// ---------------------------------------------------------------------------
-
-namespace
-{
-
-// Build a unit-square 2D simplex mesh.
-//   diagonal=0: T0=(0,1,3), T1=(1,2,3)
-//   diagonal=1: T0=(0,1,2), T1=(0,2,3)
-Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal)
-{
-  const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0});
-  Omega_h::LOs ev2v = (diagonal == 0) ? Omega_h::LOs({0, 1, 3, 1, 2, 3})
-                                      : Omega_h::LOs({0, 1, 2, 0, 2, 3});
-  Omega_h::Mesh mesh(&lib);
-  Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords);
-  mesh.add_tag<Omega_h::I8>(
-    0, "class_dim", 1,
-    Omega_h::Read<Omega_h::I8>(mesh.nverts(), Omega_h::I8(0)));
-  mesh.add_tag<Omega_h::ClassId>(
-    0, "class_id", 1,
-    Omega_h::Read<Omega_h::ClassId>(mesh.nverts(), Omega_h::ClassId(0)));
-  mesh.add_tag<Omega_h::I8>(
-    1, "class_dim", 1,
-    Omega_h::Read<Omega_h::I8>(mesh.nedges(), Omega_h::I8(1)));
-  mesh.add_tag<Omega_h::ClassId>(
-    1, "class_id", 1,
-    Omega_h::Read<Omega_h::ClassId>(mesh.nedges(), Omega_h::ClassId(0)));
-  mesh.add_tag<Omega_h::I8>(
-    2, "class_dim", 1,
-    Omega_h::Read<Omega_h::I8>(mesh.nelems(), Omega_h::I8(2)));
-  mesh.add_tag<Omega_h::ClassId>(
-    2, "class_id", 1,
-    Omega_h::Read<Omega_h::ClassId>(mesh.nelems(), Omega_h::ClassId(0)));
-  return mesh;
-}
-
-std::shared_ptr<pcms::LagrangeFunctionSpace> MakeP1Space(Omega_h::Mesh& mesh)
-{
-  return pcms::LagrangeFunctionSpace::FromMesh(
-    mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
-}
-
-// Evaluates source_field at the integrator's sample points and assembles.
-void EvaluateAndAssemble(
-  pcms::LinearFormIntegrator& integrator,
-  const std::shared_ptr<const pcms::FunctionSpace>& source_space,
-  const pcms::Field<pcms::Real>& source_field)
-{
-  const auto& pts = integrator.GetIntegrationPoints();
-  const std::size_t npts = pts.GetValues().extent(0);
-  auto evaluator = source_space->CreatePointEvaluator<pcms::Real>(
-    pcms::EvaluationRequest::FromCoordinates(pts));
-  Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace> sampled("sampled", npts,
-                                                              1);
-  evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled));
-  integrator.Assemble(pcms::MakeConstRank2View(sampled));
-}
-
-} // namespace
-
 // ---------------------------------------------------------------------------
 // Monte Carlo RHS integrator
 // ---------------------------------------------------------------------------
@@ -86,20 +23,17 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: sample points lie inside the domain",
           "[mc_rhs_integrator]")
 {
   Omega_h::Library lib;
-  auto target_mesh = BuildUnitSquare(lib, 0);
-  auto target_space = MakeP1Space(target_mesh);
+  auto target_mesh = pcms::test::BuildUnitSquare(lib, 0);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   const int samples_per_element = 16;
   for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) {
     pcms::OmegaHMonteCarloRHSIntegrator integrator(
       *target_space, samples_per_element, sampling);
     const auto raw_coords = integrator.GetIntegrationPoints().GetValues();
-    auto coords_view = Kokkos::View<const pcms::Real**, Kokkos::LayoutRight,
-                                    pcms::DeviceMemorySpace,
-                                    Kokkos::MemoryTraits<Kokkos::Unmanaged>>(
-      raw_coords.data_handle(), raw_coords.extent(0), raw_coords.extent(1));
-    auto coords_h =
-      Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coords_view);
+    auto coords_h = pcms::test::CopyCoordinatesToHost(
+      raw_coords, static_cast<int>(raw_coords.extent(0)),
+      static_cast<int>(raw_coords.extent(1)));
 
     REQUIRE(
       coords_h.extent(0) ==
@@ -120,10 +54,10 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly",
   // For f = c the estimator sums to c * |domain| for any sample placement,
   // because the P1 basis functions partition unity at every sample point.
   Omega_h::Library lib;
-  auto source_mesh = BuildUnitSquare(lib, 1);
-  auto target_mesh = BuildUnitSquare(lib, 0);
-  auto source_space = MakeP1Space(source_mesh);
-  auto target_space = MakeP1Space(target_mesh);
+  auto source_mesh = pcms::test::BuildUnitSquare(lib, 1);
+  auto target_mesh = pcms::test::BuildUnitSquare(lib, 0);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   auto source_field = source_space->CreateFunction<pcms::Real>();
   pcms::test::SetField(
@@ -131,7 +65,7 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly",
 
   for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) {
     pcms::OmegaHMonteCarloRHSIntegrator integrator(*target_space, 8, sampling);
-    EvaluateAndAssemble(integrator, source_space, source_field);
+    pcms::test::EvaluateAndAssemble(integrator, source_space, source_field);
 
     PetscScalar sum = 0.0;
     VecSum(integrator.GetVector(), &sum);
@@ -152,10 +86,10 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target "
   // of the source field) equals the source field, so the sampled residual is
   // identically zero and the projection is exact regardless of sample count.
   Omega_h::Library lib;
-  auto source_mesh = BuildUnitSquare(lib, 1);
-  auto target_mesh = BuildUnitSquare(lib, 0);
-  auto source_space = MakeP1Space(source_mesh);
-  auto target_space = MakeP1Space(target_mesh);
+  auto source_mesh = pcms::test::BuildUnitSquare(lib, 1);
+  auto target_mesh = pcms::test::BuildUnitSquare(lib, 0);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   auto source_field = source_space->CreateFunction<pcms::Real>();
   auto target_field = target_space->CreateFunction<pcms::Real>();
@@ -170,11 +104,11 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target "
 
   const auto values =
     pcms::FlattenToRank1View(target_field.GetDOFHolderDataHost());
-  const auto coords_h = Omega_h::HostRead<Omega_h::Real>(target_mesh.coords());
+  const auto coords_h = pcms::test::CopyCoordinatesToHost(
+    pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), 2);
   REQUIRE(static_cast<Omega_h::LO>(values.size()) == target_mesh.nverts());
   for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) {
-    const double expected =
-      3.0 * coords_h[2 * i + 0] - coords_h[2 * i + 1] + 0.25;
+    const double expected = 3.0 * coords_h(i, 0) - coords_h(i, 1) + 0.25;
     CAPTURE(i, expected, values[i]);
     CHECK(values[i] == Catch::Approx(expected).margin(1e-9));
   }
@@ -188,10 +122,10 @@ TEST_CASE("OmegaHControlVariateProjection: reduces error vs plain Monte Carlo",
   // seeds the control-variate projection must be closer to the exact
   // (intersection-quadrature) projection than the plain Monte Carlo one.
   Omega_h::Library lib;
-  auto source_mesh = BuildUnitSquare(lib, 1);
-  auto target_mesh = BuildUnitSquare(lib, 0);
-  auto source_space = MakeP1Space(source_mesh);
-  auto target_space = MakeP1Space(target_mesh);
+  auto source_mesh = pcms::test::BuildUnitSquare(lib, 1);
+  auto target_mesh = pcms::test::BuildUnitSquare(lib, 0);
+  auto source_space = pcms::test::MakeP1Space(source_mesh);
+  auto target_space = pcms::test::MakeP1Space(target_mesh);
 
   auto source_field = source_space->CreateFunction<pcms::Real>();
   pcms::test::SetField(
diff --git a/test/test_point_evaluator.cpp b/test/test_point_evaluator.cpp
index dcc57087..e9e235a9 100644
--- a/test/test_point_evaluator.cpp
+++ b/test/test_point_evaluator.cpp
@@ -34,7 +34,7 @@ TEST_CASE("PointEvaluator: OmegaH order-1 linear evaluation")
   auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
     field_data.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 
   auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
@@ -44,7 +44,7 @@ TEST_CASE("PointEvaluator: OmegaH order-1 linear evaluation")
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   pcms::test::CheckEvaluation(
     *evaluator, field_data, pts,
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 }
 
 // ============================================================================
@@ -67,45 +67,24 @@ TEST_CASE("PointEvaluator: same evaluator reused for two FieldData objects")
   // field_a: linear_f;  field_b: constant 42
   pcms::test::SetField(
     field_a.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   pcms::test::SetField(
     field_b.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real, Real) { return Real(42); });
 
   auto pts = pcms::test::StandardEvalCoords2D();
-  int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
 
-  // Create the PointEvaluator once
+  // Create the PointEvaluator once and reuse it for both fields.
   auto evaluator = factory->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
-  Kokkos::View<Real*, pcms::DeviceMemorySpace> out_a_device("out_a", n);
-  Kokkos::View<Real*, pcms::DeviceMemorySpace> out_b_device("out_b", n);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<pcms::DeviceMemorySpace>;
-  auto view_a = pcms::Rank2View<Real, pcms::DeviceMemorySpace, LayoutPolicy>(
-    out_a_device.data(), n, 1);
-  auto view_b = pcms::Rank2View<Real, pcms::DeviceMemorySpace, LayoutPolicy>(
-    out_b_device.data(), n, 1);
-
-  // Evaluate field_a then field_b with the same evaluator
-  evaluator->Evaluate(field_a, view_a);
-  evaluator->Evaluate(field_b, view_b);
-
-  auto out_a_host =
-    Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_a_device);
-  auto out_b_host =
-    Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_b_device);
-
-  for (int i = 0; i < n; ++i) {
-    Real x = pts[2 * static_cast<size_t>(i)],
-         y = pts[2 * static_cast<size_t>(i) + 1];
-    REQUIRE(out_a_host(i) ==
-            Catch::Approx(pcms::test::linear_f(x, y)).margin(1e-10));
-    REQUIRE(out_b_host(i) == Catch::Approx(42.0).margin(1e-10));
-  }
+  pcms::test::CheckEvaluation(
+    *evaluator, field_a, pts,
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
+  pcms::test::CheckEvaluation(
+    *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return Real(42); });
 }
 
 // ============================================================================
@@ -125,11 +104,10 @@ TEST_CASE("PointEvaluator: OmegaH order-1 out-of-bounds fill")
   auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
     field_data.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 
   // Points clearly outside [0,1]^2
-  const std::vector<Real> outside_pts = {-0.5, 0.5,  1.5, 0.5,
-                                         0.5,  -0.5, 0.5, 1.5};
+  const auto outside_pts = pcms::test::StandardOutsideCoords2D();
   auto device_coords = pcms::test::CreateDeviceCoordinateView(
     outside_pts, CoordinateSystem::Cartesian);
   pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0};
@@ -158,7 +136,7 @@ TEST_CASE("PointEvaluator: UniformGrid order-1 linear evaluation")
   auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
     field_data.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
@@ -167,7 +145,8 @@ TEST_CASE("PointEvaluator: UniformGrid order-1 linear evaluation")
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   pcms::test::CheckEvaluation(
     *evaluator, field_data, pts,
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 1e-8);
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); },
+    1e-8);
 }
 
 TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation")
@@ -184,7 +163,7 @@ TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation")
   auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
     field_data.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
@@ -193,7 +172,8 @@ TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation")
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   pcms::test::CheckEvaluation(
     *evaluator, field_data, pts,
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 1e-8);
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); },
+    1e-8);
 }
 
 // ============================================================================
@@ -308,7 +288,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 linear evaluation")
   auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
     field_data.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 
   auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
@@ -318,7 +298,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 linear evaluation")
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   pcms::test::CheckEvaluation(
     *evaluator, field_data, pts,
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 }
 
 TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill")
@@ -334,10 +314,9 @@ TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill")
   auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
     field_data.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
 
-  const std::vector<Real> outside_pts = {-0.5, 0.5,  1.5, 0.5,
-                                         0.5,  -0.5, 0.5, 1.5};
+  const auto outside_pts = pcms::test::StandardOutsideCoords2D();
   auto device_coords = pcms::test::CreateDeviceCoordinateView(
     outside_pts, CoordinateSystem::Cartesian);
   pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0};
@@ -362,43 +341,24 @@ TEST_CASE(
   auto field_b = factory->CreateFunction<Real>();
   pcms::test::SetField(
     field_a.GetData(), *factory->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   pcms::test::SetField(
     field_b.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real, Real) { return Real(42); });
 
   auto pts = pcms::test::StandardEvalCoords2D();
-  int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
 
+  // Create the PointEvaluator once and reuse it for both fields.
   auto evaluator = factory->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
-  Kokkos::View<Real*, pcms::DeviceMemorySpace> out_a_device("out_a", n);
-  Kokkos::View<Real*, pcms::DeviceMemorySpace> out_b_device("out_b", n);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<pcms::DeviceMemorySpace>;
-  auto view_a = pcms::Rank2View<Real, pcms::DeviceMemorySpace, LayoutPolicy>(
-    out_a_device.data(), n, 1);
-  auto view_b = pcms::Rank2View<Real, pcms::DeviceMemorySpace, LayoutPolicy>(
-    out_b_device.data(), n, 1);
-
-  evaluator->Evaluate(field_a, view_a);
-  evaluator->Evaluate(field_b, view_b);
-
-  auto out_a_host =
-    Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_a_device);
-  auto out_b_host =
-    Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_b_device);
-
-  for (int i = 0; i < n; ++i) {
-    Real x = pts[2 * static_cast<size_t>(i)],
-         y = pts[2 * static_cast<size_t>(i) + 1];
-    REQUIRE(out_a_host(i) ==
-            Catch::Approx(pcms::test::linear_f(x, y)).margin(1e-10));
-    REQUIRE(out_b_host(i) == Catch::Approx(42.0).margin(1e-10));
-  }
+  pcms::test::CheckEvaluation(
+    *evaluator, field_a, pts,
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
+  pcms::test::CheckEvaluation(
+    *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return Real(42); });
 }
 
 TEST_CASE("LagrangeFunctionSpace: MeshFields rejects multi-component fields")
diff --git a/test/test_polynomial_reconstruction_mls_evaluation.cpp b/test/test_polynomial_reconstruction_mls_evaluation.cpp
index 888a83b1..415d4968 100644
--- a/test/test_polynomial_reconstruction_mls_evaluation.cpp
+++ b/test/test_polynomial_reconstruction_mls_evaluation.cpp
@@ -77,12 +77,6 @@ pcms::MLSOptions DefaultTestOptions3D()
   return opts;
 }
 
-// Interior query points for a unit box — same as StandardEvalCoords2D.
-std::vector<Real> QueryPoints()
-{
-  return pcms::test::StandardEvalCoords2D();
-}
-
 pcms::MLSOptions SweepTestOptions(unsigned degree,
                                   pcms::RadialBasisFunction basis)
 {
@@ -122,7 +116,7 @@ void CheckPolynomialReproduction(unsigned degree,
   auto field = fs->CreateFunction<Real>();
   pcms::test::SetField(field.GetData(), *fs->GetLayout(), func);
 
-  auto pts = QueryPoints();
+  auto pts = pcms::test::StandardEvalCoords2D();
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
   auto evaluator = fs->CreatePointEvaluator<Real>(
@@ -149,7 +143,9 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: reproduces "
     CheckPolynomialReproduction(
       0, basis, OMEGA_H_LAMBDA(Real, Real) { return 3.14; }, 5e-3);
     CheckPolynomialReproduction(
-      1, basis, OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 5e-3);
+      1, basis,
+      OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); },
+      5e-3);
     CheckPolynomialReproduction(
       2, basis,
       OMEGA_H_LAMBDA(Real x, Real y) { return x * x + x * y + 2.0 * y * y; },
@@ -175,43 +171,26 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: same PointEvaluator "
 
   pcms::test::SetField(
     field_a.GetData(), *fs->GetLayout(),
-    OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); });
   const Real cval = 7.0;
   pcms::test::SetField(
     field_b.GetData(), *fs->GetLayout(),
     OMEGA_H_LAMBDA(Real, Real) { return cval; });
 
-  auto pts = QueryPoints();
-  int n = static_cast<int>(pts.size()) / 2;
+  auto pts = pcms::test::StandardEvalCoords2D();
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
+  // Create the PointEvaluator once and reuse it for both fields.
   auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
-  Kokkos::View<Real*, DeviceMemorySpace> out_a_device("out_a", n);
-  Kokkos::View<Real*, DeviceMemorySpace> out_b_device("out_b", n);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<DeviceMemorySpace>;
-  Rank2View<Real, DeviceMemorySpace, LayoutPolicy> view_a(out_a_device.data(),
-                                                          n, 1);
-  Rank2View<Real, DeviceMemorySpace, LayoutPolicy> view_b(out_b_device.data(),
-                                                          n, 1);
-
-  evaluator->Evaluate(field_a, view_a);
-  evaluator->Evaluate(field_b, view_b);
-
-  auto out_a_host =
-    Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_a_device);
-  auto out_b_host =
-    Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_b_device);
-
-  for (int i = 0; i < n; ++i) {
-    Real x = pts[2 * static_cast<size_t>(i)],
-         y = pts[2 * static_cast<size_t>(i) + 1];
-    REQUIRE(out_a_host(i) ==
-            Catch::Approx(pcms::test::linear_f(x, y)).margin(5e-3));
-    REQUIRE(out_b_host(i) == Catch::Approx(cval).margin(5e-3));
-  }
+  pcms::test::CheckEvaluation(
+    *evaluator, field_a, pts,
+    OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); },
+    5e-3);
+  pcms::test::CheckEvaluation(
+    *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return cval; },
+    5e-3);
 }
 
 // ============================================================================
@@ -228,7 +207,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: Evaluate throws for "
     coords_view, CoordinateSystem::Cartesian);
   auto field = fs->CreateFunction<Real>();
 
-  auto pts = QueryPoints();
+  auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
@@ -236,13 +215,8 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: Evaluate throws for "
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
   // Two-component output — must throw
-  Kokkos::View<Real*, DeviceMemorySpace> out_device("out",
-                                                    static_cast<size_t>(n) * 2);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<DeviceMemorySpace>;
-  auto out_view =
-    Rank2View<Real, DeviceMemorySpace, LayoutPolicy>(out_device.data(), n, 2);
-  REQUIRE_THROWS(evaluator->Evaluate(field, out_view));
+  Kokkos::View<Real**, DeviceMemorySpace> out_device("out", n, 2);
+  REQUIRE_THROWS(evaluator->Evaluate(field, pcms::MakeRank2View(out_device)));
 }
 
 // ============================================================================
@@ -263,24 +237,20 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: default MLSOptions — "
     field.GetData(), *fs->GetLayout(),
     OMEGA_H_LAMBDA(Real, Real) { return Real(1.0); });
 
-  auto pts = QueryPoints();
+  auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
   auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
-  Kokkos::View<Real*, DeviceMemorySpace> out_device("out", n);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<DeviceMemorySpace>;
-  auto out_view =
-    Rank2View<Real, DeviceMemorySpace, LayoutPolicy>(out_device.data(), n, 1);
+  Kokkos::View<Real**, DeviceMemorySpace> out_device("out", n, 1);
   // Just verify it runs without error and returns finite values
-  REQUIRE_NOTHROW(evaluator->Evaluate(field, out_view));
+  REQUIRE_NOTHROW(evaluator->Evaluate(field, pcms::MakeRank2View(out_device)));
   auto out_host =
     Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device);
   for (int i = 0; i < n; ++i)
-    REQUIRE(std::isfinite(out_host(i)));
+    REQUIRE(std::isfinite(out_host(i, 0)));
 }
 
 TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator "
@@ -293,7 +263,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator "
   auto fs = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cylindrical, DefaultTestOptions());
 
-  auto pts = QueryPoints();
+  auto pts = pcms::test::StandardEvalCoords2D();
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
   REQUIRE_THROWS(fs->CreatePointEvaluator<Real>(
@@ -310,7 +280,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator "
   auto fs = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cylindrical, DefaultTestOptions());
 
-  auto pts = QueryPoints();
+  auto pts = pcms::test::StandardEvalCoords2D();
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cylindrical);
   REQUIRE_THROWS(fs->CreatePointEvaluator<Real>(
@@ -342,16 +312,12 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: radius option is "
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
   auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
-  Kokkos::View<Real*, DeviceMemorySpace> out_device("out", 1);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<DeviceMemorySpace>;
-  Rank2View<Real, DeviceMemorySpace, LayoutPolicy> out_view(out_device.data(),
-                                                            1, 1);
-  evaluator->Evaluate(field, out_view);
+  Kokkos::View<Real**, DeviceMemorySpace> out_device("out", 1, 1);
+  evaluator->Evaluate(field, pcms::MakeRank2View(out_device));
   auto out_host =
     Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device);
 
-  REQUIRE(out_host(0) == Catch::Approx(1.0).margin(1e-8));
+  REQUIRE(out_host(0, 0) == Catch::Approx(1.0).margin(1e-8));
 }
 
 TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve "
@@ -375,34 +341,24 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve "
   }
 
   auto field = fs->CreateFunction<Real>();
-  std::vector<Real> dof_values(27);
-  for (int i = 0; i < 27; ++i) {
-    const Real x = src[3 * i + 0];
-    const Real y = src[3 * i + 1];
-    const Real z = src[3 * i + 2];
-    dof_values[i] = x + 2.0 * y + 3.0 * z;
-  }
-  field.GetData().SetDOFHolderDataHost(Rank2View<const Real, HostMemorySpace>(
-    dof_values.data(), static_cast<pcms::LO>(dof_values.size()), 1));
+  pcms::test::SetField(
+    field,
+    KOKKOS_LAMBDA(Real x, Real y, Real z) { return x + 2.0 * y + 3.0 * z; });
 
   std::vector<Real> pts{0.5, 0.5, 0.25, 0.5, 0.5, 0.75};
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian, 3);
   auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
-  Kokkos::View<Real*, DeviceMemorySpace> out_device("out", 2);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<DeviceMemorySpace>;
-  Rank2View<Real, DeviceMemorySpace, LayoutPolicy> out_view(out_device.data(),
-                                                            2, 1);
-  evaluator->Evaluate(field, out_view);
+  Kokkos::View<Real**, DeviceMemorySpace> out_device("out", 2, 1);
+  evaluator->Evaluate(field, pcms::MakeRank2View(out_device));
   auto out_host =
     Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device);
 
   // This 3D case is a small, low-resolution support cloud with MLS weights
   // built from a finite-radius neighborhood rather than exact nodal lookup, so
   // it is checked with a slightly looser tolerance than the denser 2D tests.
-  REQUIRE(out_host(0) == Catch::Approx(2.25).margin(1e-2));
-  REQUIRE(out_host(1) == Catch::Approx(3.75).margin(1e-2));
-  REQUIRE(out_host(1) - out_host(0) == Catch::Approx(1.5).margin(1e-2));
+  REQUIRE(out_host(0, 0) == Catch::Approx(2.25).margin(1e-2));
+  REQUIRE(out_host(1, 0) == Catch::Approx(3.75).margin(1e-2));
+  REQUIRE(out_host(1, 0) - out_host(0, 0) == Catch::Approx(1.5).margin(1e-2));
 }
diff --git a/test/test_spr_meshfields.cpp b/test/test_spr_meshfields.cpp
index 42a6acf0..7ad4470a 100644
--- a/test/test_spr_meshfields.cpp
+++ b/test/test_spr_meshfields.cpp
@@ -3,6 +3,7 @@
 #include <pcms/localization/adj_search.hpp>
 #include <pcms/field/evaluator/mls_interpolation.hpp>
 #include <pcms/field/evaluator/pcms_interpolator_aliases.hpp>
+#include <pcms/utility/mesh_geometry.h>
 #include <Omega_h_mesh.hpp>
 #include <Omega_h_build.hpp>
 #include <Omega_h_file.hpp>
@@ -127,22 +128,7 @@ TEST_CASE("meshfields_spr_test")
 
   const auto& ntargets = mesh.nverts();
 
-  Write<Omega_h::Real> source_coordinates(
-    dim * nfaces, 0, "stores coordinates of cell centroid of each tri element");
-
-  const auto& faces2nodes = mesh.ask_down(FACE, VERT).ab2b;
-
-  Kokkos::parallel_for(
-    "calculate the centroid in each tri element", nfaces,
-    OMEGA_H_LAMBDA(const LO id) {
-      const auto current_el_verts = gather_verts<3>(faces2nodes, id);
-      const Omega_h::Few<Omega_h::Vector<2>, 3> current_el_vert_coords =
-        gather_vectors<3, 2>(target_coordinates, current_el_verts);
-      auto centroid = average(current_el_vert_coords);
-      int index = 2 * id;
-      source_coordinates[index] = centroid[0];
-      source_coordinates[index + 1] = centroid[1];
-    });
+  const auto source_coordinates = pcms::get_entity_centroids(mesh, FACE);
 
   pcms::Points source_points;
 
diff --git a/test/test_svd_serial.cpp b/test/test_svd_serial.cpp
index bdc23290..1774e29b 100644
--- a/test/test_svd_serial.cpp
+++ b/test/test_svd_serial.cpp
@@ -38,7 +38,7 @@ TEST_CASE("test_serial_svd")
 
   Kokkos::deep_copy(A_data, host_A_data);
 
-  Kokkos::View<double*> rhs_data("Device rhs data", row, column);
+  Kokkos::View<double*> rhs_data("Device rhs data", row);
   auto host_rhs_data = Kokkos::create_mirror_view(rhs_data);
 
   host_rhs_data(0) = 1.28571;
@@ -54,7 +54,7 @@ TEST_CASE("test_serial_svd")
   {
 
     Kokkos::View<double**> result("result", row, row);
-    Kokkos::View<double**> transpose_expected("result", row, row);
+    Kokkos::View<double**> transpose_expected("transpose_expected", row, row);
     Kokkos::deep_copy(result, 0.0);
     Kokkos::deep_copy(transpose_expected, 0.0);
     team_policy tp(1, Kokkos::AUTO);
diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp
index 632bd4db..699ec26a 100644
--- a/test/test_uniform_grid_field.cpp
+++ b/test/test_uniform_grid_field.cpp
@@ -72,7 +72,7 @@ void VerifyUniformGridFieldValues(
       int vertex_id = j * (grid.divisions[0] + 1) + i;
       pcms::Real x = ug_coords.GetValues()(vertex_id, 0);
       pcms::Real y = ug_coords.GetValues()(vertex_id, 1);
-      pcms::Real expected = x + 2.0 * y;
+      pcms::Real expected = pcms::test::linear_f(x, y);
       pcms::Real actual = ug_field_data[vertex_id];
       REQUIRE(std::abs(expected - actual) <= 1e-10);
     }
@@ -151,21 +151,16 @@ TEST_CASE("UniformGrid order-0 field creation and evaluation")
   auto evaluator = eval_factory.CreatePointEvaluator(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
-  Kokkos::View<pcms::Real*, pcms::HostMemorySpace> results_host("results_host",
-                                                                4);
-  Kokkos::View<pcms::Real*, pcms::DeviceMemorySpace> results_device(
-    "results_device", 4);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<pcms::DeviceMemorySpace>;
-  pcms::Rank2View<pcms::Real, pcms::DeviceMemorySpace, LayoutPolicy> out(
-    results_device.data(), 4, 1);
-  evaluator->Evaluate(field, out);
-  Kokkos::deep_copy(results_host, results_device);
-
-  REQUIRE(results_host(0) == Catch::Approx(1.0));
-  REQUIRE(results_host(1) == Catch::Approx(2.0));
-  REQUIRE(results_host(2) == Catch::Approx(3.0));
-  REQUIRE(results_host(3) == Catch::Approx(4.0));
+  Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace> results_device(
+    "results_device", 4, 1);
+  evaluator->Evaluate(field, pcms::MakeRank2View(results_device));
+  auto results_host = Kokkos::create_mirror_view_and_copy(
+    pcms::HostMemorySpace(), results_device);
+
+  REQUIRE(results_host(0, 0) == Catch::Approx(1.0));
+  REQUIRE(results_host(1, 0) == Catch::Approx(2.0));
+  REQUIRE(results_host(2, 0) == Catch::Approx(3.0));
+  REQUIRE(results_host(3, 0) == Catch::Approx(4.0));
 }
 
 TEST_CASE("UniformGrid field data operations", "[uniform_grid_field]")
@@ -236,26 +231,21 @@ TEST_CASE("UniformGrid field evaluation - piecewise constant")
   auto evaluator = eval_factory.CreatePointEvaluator(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
-  Kokkos::View<pcms::Real*, pcms::HostMemorySpace> results_host("results_host",
-                                                                4);
-  Kokkos::View<pcms::Real*, pcms::DeviceMemorySpace> results_device(
-    "results_device", 4);
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<pcms::DeviceMemorySpace>;
-  pcms::Rank2View<pcms::Real, pcms::DeviceMemorySpace, LayoutPolicy> out(
-    results_device.data(), 4, 1);
-  evaluator->Evaluate(field, out);
-  Kokkos::deep_copy(results_host, results_device);
+  Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace> results_device(
+    "results_device", 4, 1);
+  evaluator->Evaluate(field, pcms::MakeRank2View(results_device));
+  auto results_host = Kokkos::create_mirror_view_and_copy(
+    pcms::HostMemorySpace(), results_device);
 
   // Check results - interpolated from vertices
   // Cell 0 center (2.5, 2.5): avg of v0,v1,v3,v4 = (1.0+1.5+2.0+2.5)/4 = 1.75
   // Cell 1 center (7.5, 2.5): avg of v1,v2,v4,v5 = (1.5+2.0+2.5+3.0)/4 = 2.25
   // Cell 2 center (2.5, 7.5): avg of v3,v4,v6,v7 = (2.0+2.5+3.0+3.5)/4 = 2.75
   // Cell 3 center (7.5, 7.5): avg of v4,v5,v7,v8 = (2.5+3.0+3.5+4.0)/4 = 3.25
-  REQUIRE(std::abs(results_host(0) - 1.75) < 1e-10);
-  REQUIRE(std::abs(results_host(1) - 2.25) < 1e-10);
-  REQUIRE(std::abs(results_host(2) - 2.75) < 1e-10);
-  REQUIRE(std::abs(results_host(3) - 3.25) < 1e-10);
+  REQUIRE(std::abs(results_host(0, 0) - 1.75) < 1e-10);
+  REQUIRE(std::abs(results_host(1, 0) - 2.25) < 1e-10);
+  REQUIRE(std::abs(results_host(2, 0) - 2.75) < 1e-10);
+  REQUIRE(std::abs(results_host(3, 0) - 3.25) < 1e-10);
 }
 
 TEST_CASE("UniformGrid field serialization")
@@ -328,8 +318,9 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field")
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
   auto omega_h_field = omega_h_factory->CreateFunction<pcms::Real>();
   pcms::test::SetField(
-    omega_h_field,
-    OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; });
+    omega_h_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) {
+      return pcms::test::linear_f(x, y);
+    });
 
   pcms::UniformGrid<2> grid;
   grid.edge_length = {1.0, 1.0};
@@ -354,7 +345,7 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field")
   for (int i = 0; i < num_ug_nodes; ++i) {
     pcms::Real x = ug_coords_host(i, 0);
     pcms::Real y = ug_coords_host(i, 1);
-    pcms::Real expected = x + 2.0 * y;
+    pcms::Real expected = pcms::test::linear_f(x, y);
     REQUIRE(std::abs(transferred_data[i] - expected) < 1e-6);
   }
 }
@@ -641,8 +632,9 @@ TEST_CASE("UniformGrid workflow")
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
   auto omega_h_field = omega_h_factory->CreateFunction<pcms::Real>();
   pcms::test::SetField(
-    omega_h_field,
-    OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; });
+    omega_h_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) {
+      return pcms::test::linear_f(x, y);
+    });
 
   auto ug_factory = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, pcms::CoordinateSystem::Cartesian);
@@ -658,23 +650,11 @@ TEST_CASE("UniformGrid workflow")
   auto ug_coords_host_view =
     pcms::test::CopyCoordinatesToHost(ug_coords_device_view, 25, 2);
 
-  auto ug_field_data_device = ug_field.GetDOFHolderData();
-  Kokkos::View<pcms::Real*, pcms::DeviceMemorySpace> ug_field_data_device_view(
-    "", 25);
-  Kokkos::parallel_for(
-    "CopyFieldDataToView", 25, KOKKOS_LAMBDA(int i) {
-      ug_field_data_device_view(i) = ug_field_data_device(i, 0);
-    });
-  auto ug_field_data_host_view =
-    Kokkos::View<pcms::Real*, pcms::HostMemorySpace>("", 25);
-  Kokkos::deep_copy(ug_field_data_host_view, ug_field_data_device_view);
-
-  pcms::Rank2View<const pcms::Real, pcms::HostMemorySpace> ug_coords(
-    ug_coords_host_view.data(), 25, 2);
   pcms::CoordinateView<pcms::HostMemorySpace> ug_coords_view(
-    pcms::CoordinateSystem::Cartesian, ug_coords);
-  pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace> ug_field_data(
-    ug_field_data_host_view.data(), 25);
+    pcms::CoordinateSystem::Cartesian,
+    pcms::MakeConstRank2View(ug_coords_host_view));
+  const auto ug_field_data =
+    pcms::FlattenToRank1View(ug_field.GetDOFHolderDataHost());
   VerifyUniformGridFieldValues(grid, ug_coords_view, ug_field_data);
 
   VerifyMaskFieldValues(grid, mask_field);
diff --git a/test/test_xgc_reverse_classification.cpp b/test/test_xgc_reverse_classification.cpp
index 4b87fcb9..5603b52b 100644
--- a/test/test_xgc_reverse_classification.cpp
+++ b/test/test_xgc_reverse_classification.cpp
@@ -39,7 +39,7 @@ TEST_CASE("reverse classification")
   }
   auto vec = rc.Serialize();
   pcms::ReverseClassificationVertex rc_deserialized;
-  pcms::Rank1View<pcms::LO, pcms::HostMemorySpace> av{vec.data(), vec.size()};
+  auto av = pcms::make_array_view(vec);
   rc_deserialized.Deserialize(av);
   REQUIRE(rc_deserialized == rc);
 }