diff --git a/CMakeLists.txt b/CMakeLists.txt
index d9c519ca..72cfaf09 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -24,6 +24,7 @@ option(PCMS_ENABLE_CLIENT "enable the coupling client implementation" ON)
 option(PCMS_ENABLE_XGC "enable xgc field adapter" ON)
 option(PCMS_ENABLE_OMEGA_H "enable Omega_h field adapter" OFF)
 option(PCMS_ENABLE_C "Enable pcms C api" ON)
+option(PCMS_ENABLE_Python "Enable pcms Python api" OFF)
 
 option(PCMS_ENABLE_PRINT "PCMS print statements enabled" ON)
 option(PCMS_ENABLE_SPDLOG "use spdlog for logging" ON)
@@ -37,6 +38,9 @@ endif()
 # broken
 find_package(redev 4.3.0 REQUIRED)
 
+if (PCMS_ENABLE_Python)
+  set(CMAKE_POSITION_INDEPENDENT_CODE TRUE)
+endif ()
 if(PCMS_ENABLE_C)
   enable_language(C)
   option(PCMS_ENABLE_Fortran "Enable pcms fortran api" ON)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 0d174940..18b6bc61 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -54,7 +54,11 @@ if(PCMS_ENABLE_XGC)
   list(APPEND PCMS_HEADERS pcms/adapter/xgc/xgc_reverse_classification.h)
 endif()
 if(PCMS_ENABLE_OMEGA_H)
-  list(APPEND PCMS_SOURCES pcms/point_search.cpp)
+  list(APPEND PCMS_SOURCES 
+    pcms/point_search.cpp
+    pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp
+    pcms/adapter/uniform_grid/uniform_grid_field.cpp
+  )
   list(
     APPEND
     PCMS_HEADERS
@@ -63,6 +67,8 @@ if(PCMS_ENABLE_OMEGA_H)
     pcms/transfer_field2.h
     pcms/uniform_grid.h
     pcms/point_search.h
+    pcms/adapter/uniform_grid/uniform_grid_field_layout.h
+    pcms/adapter/uniform_grid/uniform_grid_field.h
   )
 endif()
 
@@ -138,6 +144,10 @@ install(
 add_library(pcms_pcms INTERFACE)
 target_link_libraries(pcms_pcms INTERFACE pcms::core)
 set_target_properties(pcms_pcms PROPERTIES EXPORT_NAME pcms)
+if (PCMS_ENABLE_Python)
+  add_subdirectory(pcms/pythonapi)
+  #target_link_libraries(pcms_pcms INTERFACE pcms::pythonapi)
+endif()
 if(PCMS_ENABLE_C)
   add_subdirectory(pcms/capi)
   target_link_libraries(pcms_pcms INTERFACE pcms::capi)
diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp b/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp
new file mode 100644
index 00000000..dc584b74
--- /dev/null
+++ b/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp
@@ -0,0 +1,270 @@
+#include "uniform_grid_field.h"
+#include "pcms/profile.h"
+#include "pcms/interpolator/linear_interpolant.hpp"
+#include "pcms/interpolator/multidimarray.hpp"
+#include "pcms/assert.h"
+#include <memory>
+
+namespace pcms
+{
+
+// Localization hint structure for UniformGrid
+template <unsigned Dim>
+struct UniformGridFieldLocalizationHint
+{
+  UniformGridFieldLocalizationHint(
+    Kokkos::View<LO*, HostMemorySpace> cell_indices,
+    Kokkos::View<Real**, HostMemorySpace> coordinates)
+    : cell_indices_(cell_indices), coordinates_(coordinates)
+  {
+  }
+
+  // Cell indices for each point as 1d array
+  Kokkos::View<LO*, HostMemorySpace> cell_indices_;
+  // Coordinates of each point
+  Kokkos::View<Real**, HostMemorySpace> coordinates_;
+};
+
+template <unsigned Dim>
+UniformGridField<Dim>::UniformGridField(
+  const UniformGridFieldLayout<Dim>& layout)
+  : layout_(layout),
+    grid_(layout.GetGrid()),
+    dof_holder_data_("dof_holder_data", static_cast<size_t>(layout.OwnedSize()))
+{
+  PCMS_FUNCTION_TIMER;
+}
+
+template <unsigned Dim>
+Rank1View<const Real, HostMemorySpace> UniformGridField<Dim>::GetDOFHolderData()
+  const
+{
+  PCMS_FUNCTION_TIMER;
+  return make_const_array_view(dof_holder_data_);
+}
+
+template <unsigned Dim>
+void UniformGridField<Dim>::SetDOFHolderData(
+  Rank1View<const Real, HostMemorySpace> data)
+{
+  PCMS_FUNCTION_TIMER;
+  PCMS_ALWAYS_ASSERT(data.size() == dof_holder_data_.size());
+
+  for (size_t i = 0; i < data.size(); ++i) {
+    dof_holder_data_[i] = data[i];
+  }
+}
+
+template <unsigned Dim>
+View<Dim, Real, HostMemorySpace> UniformGridField<Dim>::to_mdspan()
+{
+  PCMS_FUNCTION_TIMER;
+
+  if constexpr (Dim == 2) {
+    return View<Dim, Real, HostMemorySpace>(
+      dof_holder_data_.data(), grid_.divisions[0] + 1,
+      grid_.divisions[1] + 1);
+  } else if constexpr (Dim == 3) {
+    return View<Dim, Real, HostMemorySpace>(
+      dof_holder_data_.data(), grid_.divisions[0] + 1,
+      grid_.divisions[1] + 1, grid_.divisions[2] + 1);
+  } else {
+    static_assert(Dim == 2 || Dim == 3,
+                  "to_mdspan only supports 2D or 3D uniform grids");
+  }
+}
+
+template <unsigned Dim>
+View<Dim, const Real, HostMemorySpace> UniformGridField<Dim>::to_mdspan() const
+{
+  PCMS_FUNCTION_TIMER;
+
+  if constexpr (Dim == 2) {
+    return View<Dim, const Real, HostMemorySpace>(
+      dof_holder_data_.data(), grid_.divisions[0] + 1,
+      grid_.divisions[1] + 1);
+  } else if constexpr (Dim == 3) {
+    return View<Dim, const Real, HostMemorySpace>(
+      dof_holder_data_.data(), grid_.divisions[0] + 1,
+      grid_.divisions[1] + 1, grid_.divisions[2] + 1);
+  } else {
+    static_assert(Dim == 2 || Dim == 3,
+                  "to_mdspan only supports 2D or 3D uniform grids");
+  }
+}
+
+template <unsigned Dim>
+LocalizationHint UniformGridField<Dim>::GetLocalizationHint(
+  CoordinateView<HostMemorySpace> coordinate_view) const
+{
+  PCMS_FUNCTION_TIMER;
+
+  if (coordinate_view.GetCoordinateSystem() !=
+      layout_.GetDOFHolderCoordinates().GetCoordinateSystem()) {
+    throw std::runtime_error(
+      "Coordinate system mismatch in GetLocalizationHint");
+  }
+
+  auto coordinates = coordinate_view.GetCoordinates();
+  LO num_points = coordinates.extent(0);
+
+  Kokkos::View<LO*, HostMemorySpace> cell_indices("cell_indices", num_points);
+  Kokkos::View<Real**, HostMemorySpace> coords_copy("coords_copy", num_points,
+                                                    Dim);
+
+  // Find which cell each point belongs to
+  for (LO i = 0; i < num_points; ++i) {
+    Omega_h::Vector<Dim> point;
+    for (unsigned d = 0; d < Dim; ++d) {
+      point[d] = coordinates(i, d);
+      coords_copy(i, d) = coordinates(i, d);
+    }
+    cell_indices[i] = grid_.ClosestCellID(point);
+  }
+
+  auto hint = std::make_shared<UniformGridFieldLocalizationHint<Dim>>(
+    cell_indices, coords_copy);
+  return LocalizationHint{hint};
+}
+
+template <unsigned Dim>
+void UniformGridField<Dim>::Evaluate(
+  LocalizationHint location, FieldDataView<Real, HostMemorySpace> results) const
+{
+  PCMS_FUNCTION_TIMER;
+
+  if (results.GetCoordinateSystem() !=
+      layout_.GetDOFHolderCoordinates().GetCoordinateSystem()) {
+    throw std::runtime_error("Coordinate system mismatch in Evaluate");
+  }
+
+  UniformGridFieldLocalizationHint<Dim> hint =
+    *reinterpret_cast<UniformGridFieldLocalizationHint<Dim>*>(
+      location.data.get());
+
+  auto values = dof_holder_data_;
+  auto coordinates = hint.coordinates_;
+  auto cell_indices = hint.cell_indices_;
+  LO num_points = coordinates.extent(0);
+
+  Kokkos::View<LO**, HostMemorySpace> cell_dimensioned_indices(
+    "cell_dimensioned_indices", num_points, Dim);
+  // Convert cell indices to multi-dimensional indices
+  for (LO i = 0; i < num_points; ++i) {
+    auto dim_indices = grid_.GetDimensionedIndex(cell_indices[i]);
+    for (unsigned d = 0; d < Dim; ++d) {
+      cell_dimensioned_indices(i, d) = dim_indices[d];
+    }
+  }
+
+  // Dimensions for vertex grid (m+1 vertices per dimension for m cells)
+  auto cell_divisions = layout_.GetGrid().divisions;
+  IntVecView dimensions_view("dimensions", Dim);
+  auto dimensions_view_host = Kokkos::create_mirror_view(dimensions_view);
+  for (unsigned d = 0; d < Dim; ++d) {
+    dimensions_view_host(d) = cell_divisions[d] + 1;
+  }
+  Kokkos::deep_copy(dimensions_view, dimensions_view_host);
+
+  // Compute parametric coordinates for each point
+  RealMatView parametric_coords("parametric_coords", num_points, Dim);
+  auto parametric_coords_host = Kokkos::create_mirror_view(parametric_coords);
+
+  for (LO i = 0; i < num_points; ++i) {
+    auto cell_bbox = grid_.GetCellBBOX(cell_indices[i]);
+    for (unsigned d = 0; d < Dim; ++d) {
+      Real coord = coordinates(i, d);
+      Real cell_min = cell_bbox.center[d] - cell_bbox.half_width[d];
+      Real cell_max = cell_bbox.center[d] + cell_bbox.half_width[d];
+      parametric_coords_host(i, d) = (coord - cell_min) / (cell_max - cell_min);
+    }
+  }
+
+  Kokkos::deep_copy(parametric_coords, parametric_coords_host);
+
+  // Convert cell indices and values to the required view types
+  IntMatView cell_indices_interp("cell_indices_interp", num_points, Dim);
+  auto cell_indices_interp_host =
+    Kokkos::create_mirror_view(cell_indices_interp);
+  for (LO i = 0; i < num_points; ++i) {
+    for (unsigned d = 0; d < Dim; ++d) {
+      cell_indices_interp_host(i, d) = cell_dimensioned_indices(i, d);
+    }
+  }
+  Kokkos::deep_copy(cell_indices_interp, cell_indices_interp_host);
+
+  RealVecView values_interp("values_interp", values.extent(0));
+  Kokkos::deep_copy(values_interp, values);
+
+  auto interpolator = RegularGridInterpolator(
+    parametric_coords, values_interp, cell_indices_interp, dimensions_view);
+  auto evaluated_values = results.GetValues();
+  auto interpolated_values = interpolator.linear_interpolation();
+  auto interpolated_values_host =
+    Kokkos::create_mirror_view(interpolated_values);
+  Kokkos::deep_copy(interpolated_values_host, interpolated_values);
+  for (LO i = 0; i < evaluated_values.size(); ++i) {
+    evaluated_values[i] = interpolated_values_host[i];
+  }
+}
+
+template <unsigned Dim>
+void UniformGridField<Dim>::EvaluateGradient(
+  FieldDataView<Real, HostMemorySpace>)
+{
+  throw std::runtime_error("Not implemented");
+}
+
+template <unsigned Dim>
+const FieldLayout& UniformGridField<Dim>::GetLayout() const
+{
+  return layout_;
+}
+
+template <unsigned Dim>
+bool UniformGridField<Dim>::CanEvaluateGradient()
+{
+  return false;
+}
+
+template <unsigned Dim>
+int UniformGridField<Dim>::Serialize(
+  Rank1View<Real, pcms::HostMemorySpace> buffer,
+  Rank1View<const pcms::LO, pcms::HostMemorySpace> permutation) const
+{
+  PCMS_FUNCTION_TIMER;
+
+  const auto array_h = GetDOFHolderData();
+  if (buffer.size() > 0) {
+    PCMS_ALWAYS_ASSERT(buffer.size() == array_h.size());
+    for (LO i = 0; i < array_h.size(); ++i) {
+      buffer[permutation[i]] = array_h[i];
+    }
+  }
+  return array_h.size();
+}
+
+template <unsigned Dim>
+void UniformGridField<Dim>::Deserialize(
+  Rank1View<const Real, pcms::HostMemorySpace> buffer,
+  Rank1View<const pcms::LO, pcms::HostMemorySpace> permutation)
+{
+  PCMS_FUNCTION_TIMER;
+
+  Kokkos::View<Real*, HostMemorySpace> sorted_buffer("sorted_buffer",
+                                                     permutation.size());
+  auto owned = layout_.GetOwned();
+
+  for (LO i = 0; i < sorted_buffer.size(); ++i) {
+    PCMS_ALWAYS_ASSERT(owned[i]);
+    sorted_buffer[i] = buffer[permutation[i]];
+  }
+
+  SetDOFHolderData(pcms::make_const_array_view(sorted_buffer));
+}
+
+// Explicit template instantiations
+template class UniformGridField<2>;
+template class UniformGridField<3>;
+
+} // namespace pcms
diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field.h b/src/pcms/adapter/uniform_grid/uniform_grid_field.h
new file mode 100644
index 00000000..eba81048
--- /dev/null
+++ b/src/pcms/adapter/uniform_grid/uniform_grid_field.h
@@ -0,0 +1,57 @@
+#ifndef PCMS_UNIFORM_GRID_FIELD_H
+#define PCMS_UNIFORM_GRID_FIELD_H
+
+#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h"
+#include "pcms/types.h"
+#include "pcms/field.h"
+#include "pcms/coordinate_system.h"
+#include "pcms/uniform_grid.h"
+#include <memory>
+
+namespace pcms
+{
+template <unsigned Dim = 2>
+class UniformGridField : public FieldT<Real>
+{
+public:
+  UniformGridField(const UniformGridFieldLayout<Dim>& layout);
+
+  LocalizationHint GetLocalizationHint(
+    CoordinateView<HostMemorySpace> coordinate_view) const override;
+
+  void Evaluate(LocalizationHint location,
+                FieldDataView<Real, HostMemorySpace> results) const override;
+
+  void EvaluateGradient(FieldDataView<Real, HostMemorySpace> results) override;
+
+  const FieldLayout& GetLayout() const override;
+
+  bool CanEvaluateGradient() override;
+
+  int Serialize(Rank1View<Real, pcms::HostMemorySpace> buffer,
+                Rank1View<const pcms::LO, pcms::HostMemorySpace> permutation)
+    const override;
+
+  void Deserialize(
+    Rank1View<const Real, pcms::HostMemorySpace> buffer,
+    Rank1View<const pcms::LO, pcms::HostMemorySpace> permutation) override;
+
+  Rank1View<const Real, HostMemorySpace> GetDOFHolderData() const override;
+  void SetDOFHolderData(Rank1View<const Real, HostMemorySpace> data) override;
+
+  View<Dim, Real, HostMemorySpace> to_mdspan();
+  View<Dim, const Real, HostMemorySpace> to_mdspan() const;
+
+  ~UniformGridField() noexcept = default;
+
+private:
+  const UniformGridFieldLayout<Dim>& layout_;
+  UniformGrid<Dim>& grid_;
+  Kokkos::View<Real*, HostMemorySpace> dof_holder_data_;
+};
+
+using UniformGridField2D = UniformGridField<2>;
+
+} // namespace pcms
+
+#endif // PCMS_UNIFORM_GRID_FIELD_H
diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp
new file mode 100644
index 00000000..e8a52315
--- /dev/null
+++ b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp
@@ -0,0 +1,162 @@
+#include "uniform_grid_field_layout.h"
+#include "uniform_grid_field.h"
+#include "pcms/profile.h"
+#include <memory>
+
+namespace pcms
+{
+
+template <unsigned Dim>
+UniformGridFieldLayout<Dim>::UniformGridFieldLayout(
+  UniformGrid<Dim>& grid, int num_components,
+  CoordinateSystem coordinate_system)
+  : grid_(grid),
+    num_components_(num_components),
+    coordinate_system_(coordinate_system),
+    gids_("gids", GetNumVertices()),
+    dof_holder_coords_("dof_holder_coords", GetNumVertices(), Dim),
+    owned_("owned", GetNumVertices())
+{
+  PCMS_FUNCTION_TIMER;
+
+  LO num_vertices = GetNumVertices();
+
+  // Initialize global IDs and ownership
+  for (LO i = 0; i < num_vertices; ++i) {
+    gids_[i] = static_cast<GO>(i);
+    owned_[i] = true;
+  }
+
+  // Initialize DOF holder coordinates at grid vertices
+  Real vertex_spacing[Dim];
+  for (unsigned d = 0; d < Dim; ++d) {
+    vertex_spacing[d] = grid_.edge_length[d] / grid_.divisions[d];
+  }
+
+  if constexpr (Dim == 2) {
+    LO vertex_idx = 0;
+    for (LO j = 0; j <= grid_.divisions[1]; ++j) {
+      for (LO i = 0; i <= grid_.divisions[0]; ++i) {
+        dof_holder_coords_(vertex_idx, 0) =
+          grid_.bot_left[0] + i * vertex_spacing[0];
+        dof_holder_coords_(vertex_idx, 1) =
+          grid_.bot_left[1] + j * vertex_spacing[1];
+        ++vertex_idx;
+      }
+    }
+  } else if constexpr (Dim == 3) {
+    LO vertex_idx = 0;
+    for (LO k = 0; k <= grid_.divisions[2]; ++k) {
+      for (LO j = 0; j <= grid_.divisions[1]; ++j) {
+        for (LO i = 0; i <= grid_.divisions[0]; ++i) {
+          dof_holder_coords_(vertex_idx, 0) =
+            grid_.bot_left[0] + i * vertex_spacing[0];
+          dof_holder_coords_(vertex_idx, 1) =
+            grid_.bot_left[1] + j * vertex_spacing[1];
+          dof_holder_coords_(vertex_idx, 2) =
+            grid_.bot_left[2] + k * vertex_spacing[2];
+          ++vertex_idx;
+        }
+      }
+    }
+  }
+}
+
+template <unsigned Dim>
+std::unique_ptr<FieldT<Real>> UniformGridFieldLayout<Dim>::CreateField() const
+{
+  return std::make_unique<UniformGridField<Dim>>(*this);
+}
+
+template <unsigned Dim>
+int UniformGridFieldLayout<Dim>::GetNumComponents() const
+{
+  return num_components_;
+}
+
+template <unsigned Dim>
+LO UniformGridFieldLayout<Dim>::GetNumOwnedDofHolder() const
+{
+  return GetNumVertices();
+}
+
+template <unsigned Dim>
+GO UniformGridFieldLayout<Dim>::GetNumGlobalDofHolder() const
+{
+  return GetNumVertices();
+}
+
+template <unsigned Dim>
+Rank1View<const bool, HostMemorySpace> UniformGridFieldLayout<Dim>::GetOwned()
+  const
+{
+  return make_const_array_view(owned_);
+}
+
+template <unsigned Dim>
+GlobalIDView<HostMemorySpace> UniformGridFieldLayout<Dim>::GetGids() const
+{
+  return GlobalIDView<HostMemorySpace>(gids_.data(), gids_.size());
+}
+
+template <unsigned Dim>
+CoordinateView<HostMemorySpace>
+UniformGridFieldLayout<Dim>::GetDOFHolderCoordinates() const
+{
+  Rank2View<const Real, HostMemorySpace> coords_view(
+    dof_holder_coords_.data(), dof_holder_coords_.extent(0), Dim);
+  return CoordinateView<HostMemorySpace>{coordinate_system_, coords_view};
+}
+
+template <unsigned Dim>
+bool UniformGridFieldLayout<Dim>::IsDistributed()
+{
+  return false;
+}
+
+template <unsigned Dim>
+UniformGrid<Dim>& UniformGridFieldLayout<Dim>::GetGrid() const
+{
+  return grid_;
+}
+
+template <unsigned Dim>
+LO UniformGridFieldLayout<Dim>::GetNumCells() const
+{
+  return grid_.GetNumCells();
+}
+
+template <unsigned Dim>
+LO UniformGridFieldLayout<Dim>::GetNumVertices() const
+{
+  LO num_vertices = 1;
+  for (unsigned d = 0; d < Dim; ++d) {
+    num_vertices *= (grid_.divisions[d] + 1);
+  }
+  return num_vertices;
+}
+
+template <unsigned Dim>
+EntOffsetsArray UniformGridFieldLayout<Dim>::GetEntOffsets() const
+{
+  EntOffsetsArray offsets{};
+  offsets[0] = 0;
+  offsets[1] = grid_.GetNumCells();
+  offsets[2] = grid_.GetNumCells();
+  offsets[3] = grid_.GetNumCells();
+  offsets[4] = grid_.GetNumCells();
+  return offsets;
+}
+
+template <unsigned Dim>
+ReversePartitionMap2 UniformGridFieldLayout<Dim>::GetReversePartitionMap(
+  const redev::Partition& partition) const
+{
+  throw std::runtime_error("Unimplemented");
+}
+
+// Explicit template instantiations
+template class UniformGridFieldLayout<2>;
+template class UniformGridFieldLayout<3>;
+
+} // namespace pcms
diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h
new file mode 100644
index 00000000..c87c7733
--- /dev/null
+++ b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h
@@ -0,0 +1,54 @@
+#ifndef PCMS_UNIFORM_GRID_FIELD_LAYOUT_H
+#define PCMS_UNIFORM_GRID_FIELD_LAYOUT_H
+
+#include "pcms/arrays.h"
+#include "pcms/field_layout.h"
+#include "pcms/coordinate_system.h"
+#include "pcms/field.h"
+#include "pcms/uniform_grid.h"
+
+#include <array>
+
+namespace pcms
+{
+template <unsigned Dim = 2>
+class UniformGridFieldLayout : public FieldLayout
+{
+public:
+  UniformGridFieldLayout(UniformGrid<Dim>& grid, int num_components,
+                         CoordinateSystem coordinate_system);
+
+  std::unique_ptr<FieldT<Real>> CreateField() const override;
+
+  int GetNumComponents() const override;
+  LO GetNumOwnedDofHolder() const override;
+  GO GetNumGlobalDofHolder() const override;
+
+  Rank1View<const bool, HostMemorySpace> GetOwned() const override;
+  GlobalIDView<HostMemorySpace> GetGids() const override;
+  CoordinateView<HostMemorySpace> GetDOFHolderCoordinates() const override;
+
+  bool IsDistributed() override;
+
+  EntOffsetsArray GetEntOffsets() const override;
+
+  ReversePartitionMap2 GetReversePartitionMap(
+    const redev::Partition& partition) const override;
+
+  UniformGrid<Dim>& GetGrid() const;
+  LO GetNumCells() const;
+  LO GetNumVertices() const;
+
+private:
+  UniformGrid<Dim>& grid_;
+  int num_components_;
+  CoordinateSystem coordinate_system_;
+  Kokkos::View<GO*, HostMemorySpace> gids_;
+  Kokkos::View<Real**, HostMemorySpace> dof_holder_coords_;
+  Kokkos::View<bool*, HostMemorySpace> owned_;
+};
+
+using UniformGridFieldLayout2D = UniformGridFieldLayout<2>;
+
+} // namespace pcms
+#endif // PCMS_UNIFORM_GRID_FIELD_LAYOUT_H
diff --git a/src/pcms/create_field.cpp b/src/pcms/create_field.cpp
index 95db5d8f..c864ddcd 100644
--- a/src/pcms/create_field.cpp
+++ b/src/pcms/create_field.cpp
@@ -1,7 +1,12 @@
 #include "create_field.h"
 #include "adapter/omega_h/omega_h_field2.h"
 #include "adapter/omega_h/omega_h_field_layout.h"
+#include "adapter/uniform_grid/uniform_grid_field.h"
+#include "adapter/uniform_grid/uniform_grid_field_layout.h"
+#include "uniform_grid.h"
+#include "point_search.h"
 
+#include <Kokkos_Core.hpp>
 #include <utility>
 
 namespace pcms
@@ -24,4 +29,93 @@ std::unique_ptr<FieldLayout> CreateLagrangeLayout(
                                              num_components, coordinate_system);
 }
 
+template <>
+std::pair<std::unique_ptr<UniformGridFieldLayout<2>>,
+          std::unique_ptr<UniformGridField<2>>>
+CreateUniformGridBinaryFieldFromGrid<2>(Omega_h::Mesh& mesh,
+                                         UniformGrid<2>& grid)
+{
+  constexpr unsigned dim = 2;
+
+  // Get total number of vertices
+  const LO num_vertices = (grid.divisions[0] + 1) * (grid.divisions[1] + 1);
+
+  // Create GridPointSearch for point-in-mesh queries
+  GridPointSearch point_search(mesh, grid.divisions[0], grid.divisions[1]);
+
+  // Create array of grid vertex points
+  Kokkos::View<Real* [dim]> vertices("vertices", num_vertices);
+  auto vertices_h = Kokkos::create_mirror_view(vertices);
+
+  // Fill vertex coordinates
+  const Real dx = grid.edge_length[0] / grid.divisions[0];
+  const Real dy = grid.edge_length[1] / grid.divisions[1];
+  
+  for (LO j = 0; j <= grid.divisions[1]; ++j) {
+    for (LO i = 0; i <= grid.divisions[0]; ++i) {
+      const LO vertex_id = j * (grid.divisions[0] + 1) + i;
+      vertices_h(vertex_id, 0) = grid.bot_left[0] + i * dx;
+      vertices_h(vertex_id, 1) = grid.bot_left[1] + j * dy;
+    }
+  }
+
+  // Copy to device
+  Kokkos::deep_copy(vertices, vertices_h);
+
+  // Perform point search
+  auto results = point_search(vertices);
+  fprintf(stderr, "Completed point-in-mesh search for %d vertices.\n",
+          num_vertices);
+
+  // Copy results back to host
+  auto results_h = Kokkos::create_mirror_view(results);
+  Kokkos::deep_copy(results_h, results);
+  fprintf(stderr, "Copied point search results back to host.\n");
+
+  // Create UniformGridFieldLayout and field
+  auto layout = std::make_unique<UniformGridFieldLayout<2>>(
+    grid, 1, CoordinateSystem::Cartesian);
+  auto field = std::make_unique<UniformGridField<2>>(*layout);
+
+  // Create binary data as Real values (0.0 or 1.0)
+  Kokkos::View<Real*, HostMemorySpace> binary_data("binary_data", num_vertices);
+  for (LO i = 0; i < num_vertices; ++i) {
+    binary_data(i) = (results_h(i).tri_id >= 0) ? 1.0 : 0.0;
+  }
+  fprintf(stderr, "Generated binary inside/outside data for grid vertices.\n");
+
+  // Set the DOF holder data
+  Rank1View<const Real, HostMemorySpace> data_view(
+    binary_data.data(), binary_data.extent(0));
+  field->SetDOFHolderData(data_view);
+  fprintf(stderr, "Set binary data on UniformGridField.\n");
+
+  return {std::move(layout), std::move(field)};
+}
+
+template <>
+std::pair<std::unique_ptr<UniformGridFieldLayout<2>>,
+          std::unique_ptr<UniformGridField<2>>>
+CreateUniformGridBinaryField<2>(Omega_h::Mesh& mesh,
+                                 const std::array<LO, 2>& divisions)
+{
+  constexpr unsigned dim = 2;
+
+  // Create the uniform grid from the mesh
+  auto grid = CreateUniformGridFromMesh<dim>(mesh, divisions);
+
+  // Delegate to the grid-based implementation
+  return CreateUniformGridBinaryFieldFromGrid<2>(mesh, grid);
+}
+
+template <>
+std::pair<std::unique_ptr<UniformGridFieldLayout<2>>,
+          std::unique_ptr<UniformGridField<2>>>
+CreateUniformGridBinaryField<2>(Omega_h::Mesh& mesh, LO cells_per_dim)
+{
+  std::array<LO, 2> divisions;
+  divisions.fill(cells_per_dim);
+  return CreateUniformGridBinaryField<2>(mesh, divisions);
+}
+
 } // namespace pcms
diff --git a/src/pcms/create_field.h b/src/pcms/create_field.h
index 8cf7ba35..0cc8b993 100644
--- a/src/pcms/create_field.h
+++ b/src/pcms/create_field.h
@@ -2,18 +2,89 @@
 #define CREATE_FIELD_H_
 
 #include "adapter/omega_h/omega_h_field_layout.h"
+#include "adapter/uniform_grid/uniform_grid_field.h"
+#include "adapter/uniform_grid/uniform_grid_field_layout.h"
 #include "field_layout.h"
 #include "field.h"
 #include "coordinate_system.h"
+#include "types.h"
+#include "uniform_grid.h"
 
 #include <Omega_h_mesh.hpp>
+#include <array>
 #include <memory>
+#include <vector>
 
 namespace pcms
 {
 std::unique_ptr<FieldLayout> CreateLagrangeLayout(
   Omega_h::Mesh& mesh, int order, int num_components,
   CoordinateSystem coordinate_system);
+
+/**
+ * \brief Create a binary field on a uniform grid indicating inside/outside mesh
+ *
+ * This function creates a uniform grid from the mesh and generates a binary
+ * field where each grid vertex is assigned:
+ *   - 1 if the vertex is inside the original mesh
+ *   - 0 if the vertex is outside the original mesh
+ *
+ * Uses GridPointSearch to determine if a point lies within the mesh domain.
+ * Currently only supports 2D meshes.
+ *
+ * \tparam dim Spatial dimension of the mesh (currently only dim=2 is supported)
+ * \param mesh The Omega_h mesh to create the grid from
+ * \param divisions Array specifying the number of cells in each dimension
+ * \return std::pair containing the layout and field (layout must outlive field)
+ *
+ * \note The returned field has vertex-centered data with values 0.0 or 1.0.
+ *       Vertex values are ordered according to the grid's internal indexing.
+ *       The layout must be kept alive as long as the field is used.
+ */
+template <unsigned dim = 2>
+std::pair<std::unique_ptr<UniformGridFieldLayout<dim>>,
+          std::unique_ptr<UniformGridField<dim>>>
+CreateUniformGridBinaryField(Omega_h::Mesh& mesh,
+                              const std::array<LO, dim>& divisions);
+
+/**
+ * \brief Create a binary field with equal divisions in all dimensions
+ *
+ * Convenience function that creates a uniform grid with the same number of
+ * cells in each dimension and generates the binary inside/outside field.
+ *
+ * \tparam dim Spatial dimension of the mesh (currently only dim=2 is supported)
+ * \param mesh The Omega_h mesh to create the grid from
+ * \param cells_per_dim Number of cells per dimension (same for all dimensions)
+ * \return std::pair containing the layout and field (layout must outlive field)
+ */
+template <unsigned dim = 2>
+std::pair<std::unique_ptr<UniformGridFieldLayout<dim>>,
+          std::unique_ptr<UniformGridField<dim>>>
+CreateUniformGridBinaryField(Omega_h::Mesh& mesh, LO cells_per_dim);
+
+/**
+ * \brief Create a binary field on a given uniform grid
+ *
+ * This function takes a pre-defined uniform grid and generates a binary field
+ * where each grid vertex is assigned:
+ *   - 1 if the vertex is inside the mesh
+ *   - 0 if the vertex is outside the mesh
+ *
+ * This allows testing with custom grids that may extend beyond the mesh
+ * boundaries.
+ *
+ * \tparam dim Spatial dimension (currently only dim=2 is supported)
+ * \param mesh The Omega_h mesh to test against
+ * \param grid The uniform grid to evaluate
+ * \return std::pair containing the layout and field (layout must outlive field)
+ */
+template <unsigned dim = 2>
+std::pair<std::unique_ptr<UniformGridFieldLayout<dim>>,
+          std::unique_ptr<UniformGridField<dim>>>
+CreateUniformGridBinaryFieldFromGrid(Omega_h::Mesh& mesh,
+                                      UniformGrid<dim>& grid);
+
 } // namespace pcms
 
 #endif // CREATE_FIELD_H_
diff --git a/src/pcms/pythonapi/CMakeLists.txt b/src/pcms/pythonapi/CMakeLists.txt
new file mode 100644
index 00000000..c33b2c70
--- /dev/null
+++ b/src/pcms/pythonapi/CMakeLists.txt
@@ -0,0 +1,71 @@
+find_package(Python COMPONENTS Interpreter Development REQUIRED)
+find_package(pybind11 REQUIRED)
+include(FindPythonModule.cmake)
+message(STATUS "Python_EXECUTABLE ${Python_EXECUTABLE}")
+find_python_module(mpi4py REQUIRED) 
+
+if(mpi4py_FOUND)
+  execute_process(
+    COMMAND
+      "${Python_EXECUTABLE}" "-c"
+      "import mpi4py as m; print(m.__version__); print(m.get_include());"
+    RESULT_VARIABLE
+      _mpi4py_SEARCH_SUCCESS
+    OUTPUT_VARIABLE
+      _mpi4py_VALUES
+    ERROR_VARIABLE
+      _mpi4py_ERROR_VALUE
+    OUTPUT_STRIP_TRAILING_WHITESPACE
+    )
+
+  # Convert the process output into a list
+  string(REGEX REPLACE ";" "\\\\;" _mpi4py_VALUES ${_mpi4py_VALUES})
+  string(REGEX REPLACE "\n" ";" _mpi4py_VALUES ${_mpi4py_VALUES})
+  list(GET _mpi4py_VALUES 0 mpi4py_VERSION)
+  list(GET _mpi4py_VALUES 1 mpi4py_INCLUDE_DIRS)
+
+  # Make sure all directory separators are '/'
+  string(REGEX REPLACE "\\\\" "/" mpi4py_INCLUDE_DIRS ${mpi4py_INCLUDE_DIRS})
+
+  # Get the major and minor version numbers
+  string(REGEX REPLACE "\\." ";" _mpi4py_VERSION_LIST ${mpi4py_VERSION})
+  list(GET _mpi4py_VERSION_LIST 0 mpi4py_VERSION_MAJOR)
+  list(GET _mpi4py_VERSION_LIST 1 mpi4py_VERSION_MINOR)
+  list(GET _mpi4py_VERSION_LIST 2 mpi4py_VERSION_PATCH)
+  string(REGEX MATCH "[0-9]*" mpi4py_VERSION_PATCH ${mpi4py_VERSION_PATCH})
+  math(EXPR mpi4py_VERSION_DECIMAL
+      "(${mpi4py_VERSION_MAJOR} * 10000) + (${mpi4py_VERSION_MINOR} * 100) + ${mpi4py_VERSION_PATCH}")
+endif()
+
+pybind11_add_module(py_pcms 
+  pythonapi.cpp
+  bind_field_base.cpp
+  bind_field_layout.cpp
+  bind_omega_h_field2.cpp
+  bind_omega_h_field_layout.cpp
+  bind_uniform_grid_field.cpp
+  bind_uniform_grid_field_layout.cpp
+  bind_omega_h.cpp
+  bind_transfer_field2.cpp
+)
+target_include_directories(py_pcms
+  SYSTEM PRIVATE
+    ${mpi4py_INCLUDE_DIRS}
+)
+set_target_properties(py_pcms PROPERTIES CXX_STANDARD 17)
+target_link_libraries(py_pcms PRIVATE Kokkos::kokkos pybind11::module MPI::MPI_C pcms::core pcms::interpolator)
+
+
+
+install(
+        TARGETS py_pcms
+        EXPORT pcms_python-targets
+        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
+        INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+        PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/pcms/pythonapi/)
+install(
+        EXPORT pcms_python-targets
+        NAMESPACE pcms::
+        DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/pcms)
\ No newline at end of file
diff --git a/src/pcms/pythonapi/FindPythonModule.cmake b/src/pcms/pythonapi/FindPythonModule.cmake
new file mode 100644
index 00000000..5583442b
--- /dev/null
+++ b/src/pcms/pythonapi/FindPythonModule.cmake
@@ -0,0 +1,82 @@
+# Find if a Python module is installed.
+# Usage: find_python_module(<package> [[ATLEAST | EXACT] version] [QUIET] [REQUIRED])
+
+macro(find_python_module module)  
+  cmake_parse_arguments(ARG
+    "QUIET;REQUIRED"  # options
+    "ATLEAST;EXACT"  # one-value arguments
+    ""  # multi-value arguments
+    ${ARGN}  # everything else
+    )
+
+  if(ARG_QUIET)
+    set(${module}_FIND_QUIETLY TRUE)
+  endif()
+
+  if(ARG_REQUIRED)
+    set(${module}_FIND_REQUIRED TRUE)
+  endif()
+
+  if(ARG_ATLEAST AND ARG_EXACT)
+    message(FATAL_ERROR "Can't be both ATLEAST and EXACT")
+  endif()
+  if(ARG_ATLEAST)
+    set(_op ">=")
+    set(${module}_tgtver ${ARG_ATLEAST})
+  elseif(ARG_EXACT)
+    set(_op "==")
+    set(${module}_tgtver ${ARG_EXACT})
+  else()
+    # deceive handle_standard_arguments into not caring about version
+    set(_${module}_requested_version_found "${Python_EXECUTABLE}")
+  endif()
+
+  unset(PY_${module} CACHE)
+  unset(${module}_VERSION CACHE)
+
+  execute_process(
+    COMMAND "${Python_EXECUTABLE}" "-c"
+            "import re; \
+             import ${module}; \
+             print(re.compile('/__init__.py.*').sub('', ${module}.__file__))"
+    RESULT_VARIABLE _${module}_status
+    OUTPUT_VARIABLE _${module}_location
+    ERROR_QUIET
+    OUTPUT_STRIP_TRAILING_WHITESPACE
+    )
+  if(NOT ${_${module}_status})
+    set(PY_${module} ${_${module}_location} CACHE STRING "Location of Python module ${module}" FORCE)
+    execute_process(
+      COMMAND "${Python_EXECUTABLE}" "-c"
+              "import sys; \
+               import ${module}; \
+               print(${module}.__version__)"
+      RESULT_VARIABLE _${module}_ver_status
+      OUTPUT_VARIABLE _${module}_version
+      ERROR_QUIET
+      OUTPUT_STRIP_TRAILING_WHITESPACE
+      )
+    if(NOT ${_${module}_ver_status})
+      set(${module}_VERSION ${_${module}_version} CACHE STRING
+          "Version of Python module ${module}" FORCE)
+
+      if(${module}_tgtver)
+        execute_process(
+          COMMAND "${Python_EXECUTABLE}" "-c"
+                  "from pkg_resources import parse_version; \
+                   print(parse_version('${${module}_VERSION}') ${_op} parse_version('${${module}_tgtver}'))"
+          RESULT_VARIABLE _${module}_verenuf_status
+          OUTPUT_VARIABLE _${module}_verenuf
+          ERROR_QUIET
+          OUTPUT_STRIP_TRAILING_WHITESPACE
+          )
+        if(NOT ${_${module}_verenuf_status})
+          if(${_${module}_verenuf} STREQUAL "True")
+            set(_${module}_requested_version_found "${Python_EXECUTABLE}")
+          endif()
+        endif()
+      endif()
+    endif()
+  endif()
+  find_package_handle_standard_args(${module} DEFAULT_MSG PY_${module} _${module}_requested_version_found)
+endmacro()
\ No newline at end of file
diff --git a/src/pcms/pythonapi/bind_field_base.cpp b/src/pcms/pythonapi/bind_field_base.cpp
new file mode 100644
index 00000000..d98e0fb3
--- /dev/null
+++ b/src/pcms/pythonapi/bind_field_base.cpp
@@ -0,0 +1,282 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/numpy.h>
+#include "pcms/coordinate_system.h"
+#include "pcms/coordinate.h"
+#include "pcms/create_field.h"
+#include "pcms/uniform_grid.h"
+#include "numpy_array_transform.h"
+
+namespace py = pybind11;
+
+namespace pcms {
+
+void bind_coordinate_system_module(py::module& m) {
+  // Bind CoordinateSystem enum
+  py::enum_<CoordinateSystem>(m, "CoordinateSystem")
+    .value("Cartesian", CoordinateSystem::Cartesian)
+    .value("Cylindrical", CoordinateSystem::Cylindrical)
+    .value("XGC", CoordinateSystem::XGC)
+    .value("GNET", CoordinateSystem::GNET)
+    .value("BEAMS3D", CoordinateSystem::BEAMS3D)
+    .export_values();
+
+  // Bind CoordinateView for HostMemorySpace
+  py::class_<CoordinateView<HostMemorySpace>>(m, "CoordinateView")
+    .def(py::init([](CoordinateSystem cs, py::array_t<Real> coords) {
+      // Convert 2D numpy array to Rank2View
+      auto buf = coords.request();
+      if (buf.ndim != 2) {
+        throw std::runtime_error("Coordinates must be a 2D array");
+      }
+      // Create a view from the numpy array
+      Rank2View<const Real, HostMemorySpace> coords_view(
+        static_cast<Real*>(buf.ptr), buf.shape[0], buf.shape[1]);
+      return CoordinateView<HostMemorySpace>(cs, coords_view);
+    }),
+    py::arg("coordinate_system"),
+    py::arg("coordinates"),
+    "Constructor for CoordinateView")
+
+    .def("get_coordinate_system", &CoordinateView<HostMemorySpace>::GetCoordinateSystem,
+         "Get the coordinate system")
+
+    .def("set_coordinate_system", &CoordinateView<HostMemorySpace>::SetCoordinateSystem,
+         py::arg("cs"),
+         "Set the coordinate system")
+
+    .def("get_coordinates", [](const CoordinateView<HostMemorySpace>& self) {
+      auto coords = self.GetCoordinates();
+      // Convert to numpy array
+      py::array_t<Real> result({static_cast<py::ssize_t>(coords.extent(0)),
+                                 static_cast<py::ssize_t>(coords.extent(1))});
+      auto buf = result.request();
+      Real* ptr = static_cast<Real*>(buf.ptr);
+      for (size_t i = 0; i < coords.extent(0); ++i) {
+        for (size_t j = 0; j < coords.extent(1); ++j) {
+          ptr[i * coords.extent(1) + j] = coords(i, j);
+        }
+      }
+      return result;
+    },
+    "Get the coordinates as numpy array");
+
+  // Bind CoordinateTransformation abstract base class
+  py::class_<CoordinateTransformation, std::shared_ptr<CoordinateTransformation>>(
+    m, "CoordinateTransformation")
+    .def("evaluate", &CoordinateTransformation::Evaluate,
+         py::arg("from"),
+         py::arg("to"),
+         "Evaluate the coordinate transformation");
+}
+
+void bind_coordinate_module(py::module& m) {
+  // Bind Coordinate template for common cases
+  // 3D Cartesian coordinates with Real type
+  py::class_<Coordinate<CoordinateSystem, Real, 3>>(m, "Coordinate3D")
+    .def(py::init<Real, Real, Real>(),
+         py::arg("x"),
+         py::arg("y"),
+         py::arg("z"),
+         "Constructor for 3D coordinate")
+
+    .def("values", [](const Coordinate<CoordinateSystem, Real, 3>& self) {
+      auto vals = self.Values();
+      py::array_t<Real> result(3);
+      auto buf = result.request();
+      Real* ptr = static_cast<Real*>(buf.ptr);
+      ptr[0] = vals[0];
+      ptr[1] = vals[1];
+      ptr[2] = vals[2];
+      return result;
+    },
+    "Get coordinate values as numpy array")
+
+    .def("__getitem__", &Coordinate<CoordinateSystem, Real, 3>::operator[],
+         py::arg("index"),
+         "Get coordinate value by index");
+
+  // 2D coordinates
+  py::class_<Coordinate<CoordinateSystem, Real, 2>>(m, "Coordinate2D")
+    .def(py::init<Real, Real>(),
+         py::arg("x"),
+         py::arg("y"),
+         "Constructor for 2D coordinate")
+
+    .def("values", [](const Coordinate<CoordinateSystem, Real, 2>& self) {
+      auto vals = self.Values();
+      py::array_t<Real> result(2);
+      auto buf = result.request();
+      Real* ptr = static_cast<Real*>(buf.ptr);
+      ptr[0] = vals[0];
+      ptr[1] = vals[1];
+      return result;
+    },
+    "Get coordinate values as numpy array")
+
+    .def("__getitem__", &Coordinate<CoordinateSystem, Real, 2>::operator[],
+         py::arg("index"),
+         "Get coordinate value by index");
+
+  // Bind CoordinateElement for common types
+  py::class_<CoordinateElement<CoordinateSystem, Real>>(m, "CoordinateElement")
+    .def(py::init<Real>(),
+         py::arg("data"),
+         "Constructor for CoordinateElement")
+
+    .def("underlying", 
+         py::overload_cast<>(&CoordinateElement<CoordinateSystem, Real>::underlying, py::const_),
+         "Get the underlying value");
+}
+
+void bind_create_field_module(py::module& m) {
+  // Bind CreateLagrangeLayout function with shared_ptr wrapper
+  // pybind11 handles shared_ptr better than unique_ptr for Python ownership
+  m.def("create_lagrange_layout", 
+        [](Omega_h::Mesh& mesh, int order, int num_components,
+           CoordinateSystem coordinate_system) {
+          return std::shared_ptr<FieldLayout>(
+            CreateLagrangeLayout(mesh, order, num_components, coordinate_system));
+        },
+        py::arg("mesh"),
+        py::arg("order"),
+        py::arg("num_components"),
+        py::arg("coordinate_system"));
+
+  // Bind CreateUniformGridFromMesh for 2D
+  m.def("create_uniform_grid_from_mesh",
+        [](Omega_h::Mesh& mesh, const std::array<LO, 2>& divisions) {
+          return CreateUniformGridFromMesh<2>(mesh, divisions);
+        },
+        py::arg("mesh"),
+        py::arg("divisions"),
+        "Create a 2D uniform grid from an Omega_h mesh");
+
+  // Bind CreateUniformGridFromMesh for 3D
+  m.def("create_uniform_grid_from_mesh",
+        [](Omega_h::Mesh& mesh, const std::array<LO, 3>& divisions) {
+          return CreateUniformGridFromMesh<3>(mesh, divisions);
+        },
+        py::arg("mesh"),
+        py::arg("divisions"),
+        "Create a 3D uniform grid from an Omega_h mesh");
+
+  // Bind CreateUniformGridBinaryField for 2D
+  m.def("create_uniform_grid_binary_field",
+        [](Omega_h::Mesh& mesh, const std::array<LO, 2>& divisions) {
+          auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, divisions);
+          // Wrap in shared_ptr for proper Python ownership and lifetime management
+          return py::make_tuple(std::shared_ptr<UniformGridFieldLayout<2>>(std::move(layout)),
+                                std::shared_ptr<UniformGridField<2>>(std::move(field)));
+        },
+        py::arg("mesh"),
+        py::arg("divisions"),
+        "Create a 2D binary field on a uniform grid indicating inside/outside mesh. "
+        "Returns tuple of (layout, field). Layout lifetime is properly managed via shared_ptr.");
+}
+
+template<typename T>
+void bind_field_t(py::module& m, const std::string& type_suffix) {
+  std::string class_name = "FieldT_" + type_suffix;
+
+  py::class_<FieldT<T>, std::shared_ptr<FieldT<T>>>(m, class_name.c_str())
+    .def("get_coordinate_system", &FieldT<T>::GetCoordinateSystem,
+         "Get the coordinate system of the field")
+
+    .def("get_localization_hint", &FieldT<T>::GetLocalizationHint,
+         py::arg("coordinates"),
+         "Get a localization hint for a set of coordinates")
+
+    .def("evaluate", [](const FieldT<T>& self,
+                       LocalizationHint hint,
+                       py::array_t<T> results_array,
+                       CoordinateSystem coord_sys) {
+      auto results_view = numpy_to_view<T>(results_array);
+      FieldDataView<T, HostMemorySpace> results(results_view, coord_sys);
+      self.Evaluate(hint, results);
+    },
+    py::arg("hint"),
+    py::arg("results"),
+    py::arg("coordinate_system"),
+    "Evaluate the field at given locations")
+
+    .def("evaluate_gradient", [](FieldT<T>& self,
+                                 py::array_t<T> results_array,
+                                 CoordinateSystem coord_sys) {
+      auto results_view = numpy_to_view<T>(results_array);
+      FieldDataView<T, HostMemorySpace> results(results_view, coord_sys);
+      self.EvaluateGradient(results);
+    },
+    py::arg("results"),
+    py::arg("coordinate_system"),
+    "Evaluate the gradient of the field")
+
+    .def("get_dof_holder_data", [](const FieldT<T>& self) {
+      auto data = self.GetDOFHolderData();
+      return view_to_numpy(data);
+    },
+    "Get the DOF holder data")
+
+    .def("set_dof_holder_data", [](FieldT<T>& self, py::array_t<const T> data) {
+      auto data_view = numpy_to_view<const T>(data);
+      self.SetDOFHolderData(data_view);
+    },
+    py::arg("data"),
+    "Set the DOF holder data")
+
+    .def("get_layout", &FieldT<T>::GetLayout,
+         py::return_value_policy::reference,
+         "Get the field layout")
+
+    .def("can_evaluate_gradient", &FieldT<T>::CanEvaluateGradient,
+         "Check if the field can evaluate gradients")
+
+    .def("serialize", [](const FieldT<T>& self,
+                        py::array_t<T> buffer,
+                        py::array_t<const LO> permutation) {
+      auto buffer_view = numpy_to_view<T>(buffer);
+      auto perm_view = numpy_to_view<const LO>(permutation);
+      return self.Serialize(buffer_view, perm_view);
+    },
+    py::arg("buffer"),
+    py::arg("permutation"),
+    "Serialize the field data")
+
+    .def("deserialize", [](FieldT<T>& self,
+                          py::array_t<const T> buffer,
+                          py::array_t<const LO> permutation) {
+      auto buffer_view = numpy_to_view<const T>(buffer);
+      auto perm_view = numpy_to_view<const LO>(permutation);
+      self.Deserialize(buffer_view, perm_view);
+    },
+    py::arg("buffer"),
+    py::arg("permutation"),
+    "Deserialize field data from buffer");
+}
+
+void bind_field_module(py::module& m) {
+  // Bind LocalizationHint (opaque type - data member is internal only)
+  py::class_<LocalizationHint>(m, "LocalizationHint")
+    .def(py::init<>());
+
+  // Bind FieldDataView for common types
+  py::class_<FieldDataView<Real, HostMemorySpace>>(m, "FieldDataView_Real")
+    .def(py::init<Rank1View<Real, HostMemorySpace>, CoordinateSystem>(),
+         py::arg("values"),
+         py::arg("coordinate_system"))
+    .def("size", &FieldDataView<Real, HostMemorySpace>::Size)
+    .def("get_coordinate_system", &FieldDataView<Real, HostMemorySpace>::GetCoordinateSystem)
+    // TODO: const GetValues?
+    .def("get_values", [](FieldDataView<Real, HostMemorySpace>& self) {
+      return view_to_numpy(self.GetValues());
+    });
+
+  // Bind FieldT for common types
+  // bind_field_t<Real>(m, "Real");
+  // bind_field_t<float>(m, "Float");
+
+  // Only bind double separately if Real is not already double
+  bind_field_t<double>(m, "Double");
+}
+
+} // namespace pcms
\ No newline at end of file
diff --git a/src/pcms/pythonapi/bind_field_layout.cpp b/src/pcms/pythonapi/bind_field_layout.cpp
new file mode 100644
index 00000000..0f7c925b
--- /dev/null
+++ b/src/pcms/pythonapi/bind_field_layout.cpp
@@ -0,0 +1,88 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/numpy.h>
+#include "pcms/field_layout.h"
+#include "numpy_array_transform.h"
+
+namespace py = pybind11;
+
+namespace pcms {
+
+void bind_field_layout_module(py::module& m) {
+  // Bind the base FieldLayout class as an abstract base
+  py::class_<FieldLayout, std::shared_ptr<FieldLayout>>(m, "FieldLayout")
+    .def("create_field", &FieldLayout::CreateField,
+         "Create a field with this layout")
+
+    .def("get_num_components", &FieldLayout::GetNumComponents,
+         "Get the number of components in the field")
+
+    .def("get_num_owned_dof_holder", &FieldLayout::GetNumOwnedDofHolder,
+         "Get the number of owned DOF holders")
+
+    .def("get_num_global_dof_holder", &FieldLayout::GetNumGlobalDofHolder,
+         "Get the number of global DOF holders")
+
+    .def("owned_size", &FieldLayout::OwnedSize,
+         "Get the owned size (num_components * num_owned_dof_holder)")
+
+    .def("global_size", &FieldLayout::GlobalSize,
+         "Get the global size (num_components * num_global_dof_holder)")
+
+    .def("get_owned", [](const FieldLayout& self) {
+      auto owned = self.GetOwned();
+      // Convert to numpy array
+      py::array_t<bool> result(owned.extent(0));
+      auto buf = result.request();
+      bool* ptr = static_cast<bool*>(buf.ptr);
+      for (size_t i = 0; i < owned.extent(0); ++i) {
+        ptr[i] = owned[i];
+      }
+      return result;
+    },
+    "Get the owned mask array")
+
+    .def("get_gids", [](const FieldLayout& self) {
+      auto gids = self.GetGids();
+      // Convert to numpy array
+      py::array_t<GO> result(gids.extent(0));
+      auto buf = result.request();
+      GO* ptr = static_cast<GO*>(buf.ptr);
+      for (size_t i = 0; i < gids.extent(0); ++i) {
+        ptr[i] = gids[i];
+      }
+      return result;
+    },
+    "Get the global IDs array")
+
+    .def("is_distributed", &FieldLayout::IsDistributed,
+         "Check if the field layout is distributed")
+
+    .def("get_ent_offsets", [](const FieldLayout& self) {
+      auto offsets = self.GetEntOffsets();
+      // Convert std::array to list
+      py::list result;
+      for (size_t i = 0; i < offsets.size(); ++i) {
+        result.append(offsets[i]);
+      }
+      return result;
+    },
+    "Get the entity offsets array")
+
+    .def("get_dof_holder_coordinates", [](const FieldLayout& self) {
+      auto coords = self.GetDOFHolderCoordinates().GetCoordinates();
+      // Convert to numpy array (2D)
+      py::array_t<Real> result({coords.extent(0), coords.extent(1)});
+      auto buf = result.request();
+      Real* ptr = static_cast<Real*>(buf.ptr);
+      for (size_t i = 0; i < coords.extent(0); ++i) {
+        for (size_t j = 0; j < coords.extent(1); ++j) {
+          ptr[i * coords.extent(1) + j] = coords(i, j);
+        }
+      }
+      return result;
+    },
+    "Get the DOF holder coordinates");
+}
+
+} // namespace pcms
\ No newline at end of file
diff --git a/src/pcms/pythonapi/bind_omega_h.cpp b/src/pcms/pythonapi/bind_omega_h.cpp
new file mode 100644
index 00000000..eaed81c8
--- /dev/null
+++ b/src/pcms/pythonapi/bind_omega_h.cpp
@@ -0,0 +1,526 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/numpy.h>
+#include <Omega_h_mesh.hpp>
+#include <Omega_h_file.hpp>
+#include <Omega_h_library.hpp>
+#include <Omega_h_build.hpp>
+#include <Omega_h_comm.hpp>
+#ifdef OMEGA_H_USE_ADIOS2
+#include <Omega_h_adios2.hpp>
+#endif
+#include "numpy_array_transform.h"
+
+namespace py = pybind11;
+
+namespace pcms {
+
+void bind_omega_h_mesh_module(py::module& m) {
+  // Bind Omega_h::Library
+  py::class_<Omega_h::Library, std::shared_ptr<Omega_h::Library>>(m, "OmegaHLibrary")
+    .def(py::init<>(), "Default constructor")
+    .def("world", &Omega_h::Library::world,
+         py::return_value_policy::reference,
+         "Get the world communicator");
+
+  // Bind the Comm class
+  py::class_<Omega_h::Comm, std::shared_ptr<Omega_h::Comm>>(m, "Comm")
+      // Constructors
+#ifdef OMEGA_H_USE_MPI
+      .def(py::init<Omega_h::Library*, MPI_Comm>(),
+            py::arg("library"),
+            py::arg("impl"))
+      .def(py::init([](Omega_h::Library* library, MPI_Comm impl, py::array_t<const Omega_h::I32> srcs,
+                       py::array_t<const Omega_h::I32> dsts) {
+        auto srcs_view = numpy_to_omega_h_read<Omega_h::I32>(srcs);
+        auto dsts_view = numpy_to_omega_h_read<Omega_h::I32>(dsts);
+        return new Omega_h::Comm(library, impl, srcs_view, dsts_view);
+      }))
+      .def("get_impl", &Omega_h::Comm::get_impl,
+            "Get the underlying MPI communicator")
+#else
+      .def(py::init<Omega_h::Library*, bool, bool>(),
+            py::arg("library"),
+            py::arg("is_graph"),
+            py::arg("sends_to_self"))
+#endif
+      // Methods
+      .def("library", &Omega_h::Comm::library,
+            py::return_value_policy::reference,
+            "Get the library pointer")
+      .def("rank", &Omega_h::Comm::rank,
+            "Get the rank of this process")
+      .def("size", &Omega_h::Comm::size,
+            "Get the total number of processes")
+      .def("dup", &Omega_h::Comm::dup,
+            "Duplicate the communicator")
+      .def("split", &Omega_h::Comm::split,
+            py::arg("color"),
+            py::arg("key"),
+            "Split the communicator")
+      .def("graph", &Omega_h::Comm::graph,
+            py::arg("dsts"),
+            "Create a graph communicator")
+      .def("graph_adjacent", &Omega_h::Comm::graph_adjacent,
+            py::arg("srcs"),
+            py::arg("dsts"),
+            "Create an adjacent graph communicator")
+      .def("graph_inverse", &Omega_h::Comm::graph_inverse,
+            "Get the inverse graph communicator")
+      .def("sources", &Omega_h::Comm::sources,
+            "Get source ranks")
+      .def("destinations", &Omega_h::Comm::destinations,
+            "Get destination ranks")
+      .def("reduce_or", &Omega_h::Comm::reduce_or,
+            py::arg("x"),
+            "Reduce using logical OR")
+      .def("reduce_and", &Omega_h::Comm::reduce_and,
+            py::arg("x"),
+            "Reduce using logical AND")
+      .def("add_int128", &Omega_h::Comm::add_int128,
+            py::arg("x"),
+            "Add Int128 values across processes")
+      .def("bcast_string", &Omega_h::Comm::bcast_string,
+            py::arg("s"),
+            py::arg("root_rank") = 0,
+            "Broadcast a string")
+      .def("barrier", &Omega_h::Comm::barrier,
+            "Synchronize all processes");
+  m.def("build_box", &Omega_h::build_box,
+        py::arg("comm"),
+        py::arg("family"),
+        py::arg("x"),
+        py::arg("y"),
+        py::arg("z"),
+        py::arg("nx"),
+        py::arg("ny"),
+        py::arg("nz"),
+        py::arg("symmetric") = false,
+        "Build a box mesh with specified dimensions and discretization",
+        py::return_value_policy::move);
+
+  // Bind Omega_h::Mesh
+  py::class_<Omega_h::Mesh, std::shared_ptr<Omega_h::Mesh>>(m, "OmegaHMesh")
+    .def(py::init<>(), "Default constructor")
+    .def(py::init<Omega_h::Library*>(), py::arg("library"),
+         "Constructor with library")
+
+    .def("set_library", &Omega_h::Mesh::set_library,
+         py::arg("library"),
+         "Set the library")
+
+    .def("library", &Omega_h::Mesh::library,
+         py::return_value_policy::reference,
+         "Get the library")
+
+    .def("set_comm", &Omega_h::Mesh::set_comm,
+         py::arg("comm"),
+         "Set the communicator")
+
+    .def("comm", &Omega_h::Mesh::comm,
+         "Get the communicator")
+
+    .def("set_dim", &Omega_h::Mesh::set_dim,
+         py::arg("dim"),
+         "Set mesh dimension")
+
+    .def("dim", &Omega_h::Mesh::dim,
+         "Get mesh dimension")
+
+    .def("set_family", &Omega_h::Mesh::set_family,
+         py::arg("family"),
+         "Set mesh family (simplex, hypercube, etc.)")
+
+    .def("family", &Omega_h::Mesh::family,
+         "Get mesh family")
+
+    .def("nverts", &Omega_h::Mesh::nverts,
+         "Get number of vertices")
+
+    .def("nedges", &Omega_h::Mesh::nedges,
+         "Get number of edges")
+
+    .def("nfaces", &Omega_h::Mesh::nfaces,
+         "Get number of faces")
+
+    .def("nregions", &Omega_h::Mesh::nregions,
+         "Get number of regions")
+
+    .def("nelems", &Omega_h::Mesh::nelems,
+         "Get number of elements")
+
+    .def("nents", &Omega_h::Mesh::nents,
+         py::arg("ent_dim"),
+         "Get number of entities of given dimension")
+
+    .def("nglobal_ents", &Omega_h::Mesh::nglobal_ents,
+         py::arg("ent_dim"),
+         "Get global number of entities")
+
+    .def("coords", [](const Omega_h::Mesh& mesh) {
+     auto coords = mesh.coords();
+     return omega_h_read_to_numpy(coords);
+    },
+    "Get mesh coordinates as numpy array")
+
+    .def("set_coords", [](Omega_h::Mesh& mesh, py::array_t<Omega_h::Real> coords) {
+      auto coords_write = numpy_to_omega_h_write<Omega_h::Real>(coords);
+      mesh.set_coords(Omega_h::Reals(coords_write));
+    },
+    py::arg("coords"),
+    "Set mesh coordinates from numpy array")
+
+    .def("add_coords", [](Omega_h::Mesh& mesh, py::array_t<Omega_h::Real> coords) {
+      auto coords_write = numpy_to_omega_h_write<Omega_h::Real>(coords);
+      mesh.add_coords(Omega_h::Reals(coords_write));
+    },
+    py::arg("coords"),
+    "Add mesh coordinates from numpy array")
+
+    .def("has_tag", &Omega_h::Mesh::has_tag,
+         py::arg("ent_dim"),
+         py::arg("name"),
+         "Check if mesh has a tag")
+
+    .def("ntags", &Omega_h::Mesh::ntags,
+         py::arg("ent_dim"),
+         "Get number of tags for entity dimension")
+
+    .def("remove_tag", &Omega_h::Mesh::remove_tag,
+         py::arg("ent_dim"),
+         py::arg("name"),
+         "Remove a tag")
+
+    .def("has_ents", &Omega_h::Mesh::has_ents,
+         py::arg("ent_dim"),
+         "Check if mesh has entities of given dimension")
+
+    .def("ask_lengths", [](Omega_h::Mesh& mesh) {
+      auto lengths = mesh.ask_lengths();
+      return omega_h_read_to_numpy(lengths);
+    },
+    "Get edge lengths")
+
+    .def("ask_qualities", [](Omega_h::Mesh& mesh) {
+      auto qualities = mesh.ask_qualities();
+      return omega_h_read_to_numpy(qualities);
+    },
+    "Get element qualities")
+
+    .def("min_quality", &Omega_h::Mesh::min_quality,
+         "Get minimum element quality")
+
+    .def("max_length", &Omega_h::Mesh::max_length,
+         "Get maximum edge length")
+
+    .def("balance", py::overload_cast<bool>(&Omega_h::Mesh::balance),
+         py::arg("predictive") = false,
+         "Balance the mesh across processors")
+
+    .def("set_parting", 
+         py::overload_cast<Omega_h_Parting, bool>(&Omega_h::Mesh::set_parting),
+         py::arg("parting"),
+         py::arg("verbose") = false,
+         "Set mesh partitioning")
+
+    .def("parting", &Omega_h::Mesh::parting,
+         "Get mesh partitioning type")
+
+    .def("nghost_layers", &Omega_h::Mesh::nghost_layers,
+         "Get number of ghost layers")
+
+    .def("owned", [](Omega_h::Mesh& mesh, Omega_h::Int ent_dim) {
+      auto owned = mesh.owned(ent_dim);
+      return omega_h_read_to_numpy(owned);
+    },
+    py::arg("ent_dim"),
+    "Get ownership flags for entities");
+
+  // Bind mesh I/O functions
+  
+  // General mesh file reader (detects format)
+  m.def("read_mesh_file",
+        [](const std::string& filepath, std::shared_ptr<Omega_h::Comm> comm) {
+          return Omega_h::read_mesh_file(filepath, comm);
+        },
+        py::arg("filepath"),
+        py::arg("comm"),
+        "Read mesh from file (auto-detects format)",
+        py::return_value_policy::move);
+
+  // Binary format I/O
+  m.def("read_mesh_binary", 
+        [](const std::string& filepath, Omega_h::Library* lib) {
+          Omega_h::Mesh mesh(lib);
+          Omega_h::binary::read(filepath, lib->world(), &mesh);
+          return mesh;
+        },
+        py::arg("filepath"),
+        py::arg("library"),
+        "Read mesh from binary file",
+        py::return_value_policy::move);
+
+  m.def("read_mesh_binary", 
+        [](const std::string& filepath, std::shared_ptr<Omega_h::Comm> comm) {
+          return Omega_h::binary::read(filepath, comm);
+        },
+        py::arg("filepath"),
+        py::arg("comm"),
+        "Read mesh from binary file",
+        py::return_value_policy::move);
+
+  m.def("write_mesh_binary",
+        [](const std::string& filepath, Omega_h::Mesh& mesh) {
+          Omega_h::binary::write(filepath, &mesh);
+        },
+        py::arg("filepath"),
+        py::arg("mesh"),
+        "Write mesh to binary file");
+
+  // Gmsh format I/O
+  m.def("read_mesh_gmsh",
+        [](const std::string& filepath, std::shared_ptr<Omega_h::Comm> comm) {
+          return Omega_h::gmsh::read(filepath, comm);
+        },
+        py::arg("filepath"),
+        py::arg("comm"),
+        "Read mesh from Gmsh file",
+        py::return_value_policy::move);
+
+  m.def("write_mesh_gmsh",
+        [](const std::string& filepath, Omega_h::Mesh& mesh) {
+          Omega_h::gmsh::write(filepath, &mesh);
+        },
+        py::arg("filepath"),
+        py::arg("mesh"),
+        "Write mesh to Gmsh file");
+
+#ifdef OMEGA_H_USE_GMSH
+  m.def("read_mesh_gmsh_parallel",
+        [](const std::string& filepath, std::shared_ptr<Omega_h::Comm> comm) {
+          return Omega_h::gmsh::read_parallel(filepath, comm);
+        },
+        py::arg("filepath"),
+        py::arg("comm"),
+        "Read parallel Gmsh mesh (MSH 4.1+)",
+        py::return_value_policy::move);
+
+  m.def("write_mesh_gmsh_parallel",
+        [](const std::string& filepath, Omega_h::Mesh& mesh) {
+          Omega_h::gmsh::write_parallel(filepath, mesh);
+        },
+        py::arg("filepath"),
+        py::arg("mesh"),
+        "Write parallel Gmsh mesh (MSH 4.1)");
+#endif
+
+  // VTK format I/O
+  m.def("write_mesh_vtu",
+        [](const std::string& filepath, Omega_h::Mesh& mesh, bool compress) {
+          Omega_h::vtk::write_vtu(filepath, &mesh, compress);
+        },
+        py::arg("filepath"),
+        py::arg("mesh"),
+        py::arg("compress") = OMEGA_H_DEFAULT_COMPRESS,
+        "Write mesh to VTU file");
+
+  m.def("write_mesh_vtu",
+        [](const std::string& filepath, Omega_h::Mesh& mesh, Omega_h::Int cell_dim, bool compress) {
+          Omega_h::vtk::write_vtu(filepath, &mesh, cell_dim, compress);
+        },
+        py::arg("filepath"),
+        py::arg("mesh"),
+        py::arg("cell_dim"),
+        py::arg("compress") = OMEGA_H_DEFAULT_COMPRESS,
+        "Write mesh to VTU file with specified cell dimension");
+
+  m.def("write_mesh_parallel_vtk",
+        [](const std::string& filepath, Omega_h::Mesh& mesh, bool compress) {
+          Omega_h::vtk::write_parallel(filepath, &mesh, compress);
+        },
+        py::arg("filepath"),
+        py::arg("mesh"),
+        py::arg("compress") = OMEGA_H_DEFAULT_COMPRESS,
+        "Write mesh to parallel VTK files");
+
+  m.def("read_mesh_parallel_vtk",
+        [](const std::string& pvtupath, std::shared_ptr<Omega_h::Comm> comm) {
+          Omega_h::Mesh mesh;
+          Omega_h::vtk::read_parallel(pvtupath, comm, &mesh);
+          return mesh;
+        },
+        py::arg("pvtupath"),
+        py::arg("comm"),
+        "Read parallel VTK mesh",
+        py::return_value_policy::move);
+
+#ifdef OMEGA_H_USE_SEACASEXODUS
+  // Exodus ClassifyWith enum
+  py::enum_<Omega_h::exodus::ClassifyWith>(m, "ExodusClassifyWith")
+    .value("NODE_SETS", Omega_h::exodus::NODE_SETS)
+    .value("SIDE_SETS", Omega_h::exodus::SIDE_SETS)
+    .export_values();
+
+  // Exodus file operations
+  m.def("exodus_open",
+        [](const std::string& filepath, bool verbose) {
+          return Omega_h::exodus::open(filepath, verbose);
+        },
+        py::arg("filepath"),
+        py::arg("verbose") = false,
+        "Open an Exodus file and return file handle");
+
+  m.def("exodus_close",
+        [](int exodus_file) {
+          Omega_h::exodus::close(exodus_file);
+        },
+        py::arg("exodus_file"),
+        "Close an Exodus file");
+
+  m.def("exodus_get_num_time_steps",
+        [](int exodus_file) {
+          return Omega_h::exodus::get_num_time_steps(exodus_file);
+        },
+        py::arg("exodus_file"),
+        "Get the number of time steps in an Exodus file");
+
+  m.def("read_mesh_exodus",
+        [](int exodus_file, Omega_h::Mesh& mesh, bool verbose, int classify_with) {
+          Omega_h::exodus::read_mesh(exodus_file, &mesh, verbose, classify_with);
+        },
+        py::arg("exodus_file"),
+        py::arg("mesh"),
+        py::arg("verbose") = false,
+        py::arg("classify_with") = Omega_h::exodus::NODE_SETS | Omega_h::exodus::SIDE_SETS,
+        "Read mesh from an open Exodus file");
+
+  m.def("read_exodus_nodal_fields",
+        [](int exodus_file, Omega_h::Mesh& mesh, int time_step, 
+           const std::string& prefix, const std::string& postfix, bool verbose) {
+          Omega_h::exodus::read_nodal_fields(exodus_file, &mesh, time_step, 
+                                              prefix, postfix, verbose);
+        },
+        py::arg("exodus_file"),
+        py::arg("mesh"),
+        py::arg("time_step"),
+        py::arg("prefix") = "",
+        py::arg("postfix") = "",
+        py::arg("verbose") = false,
+        "Read nodal fields from an open Exodus file at a specific time step");
+
+  m.def("read_exodus_element_fields",
+        [](int exodus_file, Omega_h::Mesh& mesh, int time_step,
+           const std::string& prefix, const std::string& postfix, bool verbose) {
+          Omega_h::exodus::read_element_fields(exodus_file, &mesh, time_step,
+                                                prefix, postfix, verbose);
+        },
+        py::arg("exodus_file"),
+        py::arg("mesh"),
+        py::arg("time_step"),
+        py::arg("prefix") = "",
+        py::arg("postfix") = "",
+        py::arg("verbose") = false,
+        "Read element fields from an open Exodus file at a specific time step");
+
+  m.def("read_mesh_exodus_sliced",
+        [](const std::string& filepath, std::shared_ptr<Omega_h::Comm> comm,
+           bool verbose, int classify_with, int time_step) {
+          return Omega_h::exodus::read_sliced(filepath, comm, verbose, 
+                                               classify_with, time_step);
+        },
+        py::arg("filepath"),
+        py::arg("comm"),
+        py::arg("verbose") = false,
+        py::arg("classify_with") = Omega_h::exodus::NODE_SETS | Omega_h::exodus::SIDE_SETS,
+        py::arg("time_step") = -1,
+        "Read sliced Exodus mesh in parallel",
+        py::return_value_policy::move);
+
+  m.def("write_mesh_exodus",
+        [](const std::string& filepath, Omega_h::Mesh& mesh, bool verbose, int classify_with) {
+          Omega_h::exodus::write(filepath, &mesh, verbose, classify_with);
+        },
+        py::arg("filepath"),
+        py::arg("mesh"),
+        py::arg("verbose") = false,
+        py::arg("classify_with") = Omega_h::exodus::NODE_SETS | Omega_h::exodus::SIDE_SETS,
+        "Write mesh to Exodus file");
+#endif
+
+#ifdef OMEGA_H_USE_LIBMESHB
+  // MESHB format I/O
+  m.def("read_mesh_meshb",
+        [](Omega_h::Mesh& mesh, const std::string& filepath) {
+          Omega_h::meshb::read(&mesh, filepath);
+        },
+        py::arg("mesh"),
+        py::arg("filepath"),
+        "Read mesh from MESHB file");
+
+  m.def("write_mesh_meshb",
+        [](Omega_h::Mesh& mesh, const std::string& filepath, int version) {
+          Omega_h::meshb::write(&mesh, filepath, version);
+        },
+        py::arg("mesh"),
+        py::arg("filepath"),
+        py::arg("version") = 2,
+        "Write mesh to MESHB file");
+
+  m.def("read_meshb_sol",
+        [](Omega_h::Mesh& mesh, const std::string& filepath, const std::string& sol_name) {
+          Omega_h::meshb::read_sol(&mesh, filepath, sol_name);
+        },
+        py::arg("mesh"),
+        py::arg("filepath"),
+        py::arg("sol_name"),
+        "Read solution/resolution data from MESHB .sol file");
+
+  m.def("write_meshb_sol",
+        [](Omega_h::Mesh& mesh, const std::string& filepath, const std::string& sol_name, int version) {
+          Omega_h::meshb::write_sol(&mesh, filepath, sol_name, version);
+        },
+        py::arg("mesh"),
+        py::arg("filepath"),
+        py::arg("sol_name"),
+        py::arg("version") = 2,
+        "Write solution/resolution data to MESHB .sol file");
+#endif
+
+#ifdef OMEGA_H_USE_ADIOS2
+  // ADIOS2 format I/O
+  m.def("read_mesh_adios2",
+        [](const std::string& filepath, Omega_h::Library* lib, const std::string& prefix) {
+          return Omega_h::adios::read(filepath, lib, prefix);
+        },
+        py::arg("filepath"),
+        py::arg("library"),
+        py::arg("prefix") = "",
+        "Read mesh from ADIOS2 file",
+        py::return_value_policy::move);
+
+  m.def("write_mesh_adios2",
+        [](const std::string& filepath, Omega_h::Mesh& mesh, const std::string& prefix) {
+          Omega_h::adios::write(filepath, &mesh, prefix);
+        },
+        py::arg("filepath"),
+        py::arg("mesh"),
+        py::arg("prefix") = "",
+        "Write mesh to ADIOS2 file");
+#endif
+
+  // Other utilities
+  py::enum_<Omega_h_Family>(m, "Family")
+        .value("SIMPLEX", OMEGA_H_SIMPLEX)
+        .value("HYPERCUBE", OMEGA_H_HYPERCUBE)
+        .value("MIXED", OMEGA_H_MIXED)
+        .export_values();  // This allows using values without the class prefix
+
+  // Bind Omega_h_Parting enum
+  py::enum_<Omega_h_Parting>(m, "OmegaHParting")
+    .value("ELEM_BASED", OMEGA_H_ELEM_BASED)
+    .value("VERT_BASED", OMEGA_H_VERT_BASED)
+    .value("GHOSTED", OMEGA_H_GHOSTED)
+    .export_values();
+}
+
+} // namespace pcms
\ No newline at end of file
diff --git a/src/pcms/pythonapi/bind_omega_h_field2.cpp b/src/pcms/pythonapi/bind_omega_h_field2.cpp
new file mode 100644
index 00000000..46d79fb1
--- /dev/null
+++ b/src/pcms/pythonapi/bind_omega_h_field2.cpp
@@ -0,0 +1,127 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/functional.h>
+#include <pybind11/numpy.h>
+#include "pcms/arrays.h"
+#include "pcms/adapter/omega_h/omega_h_field2.h"
+#include "pcms/adapter/omega_h/omega_h_field_layout.h"
+#include "numpy_array_transform.h"
+
+namespace py = pybind11;
+
+namespace pcms {
+
+void bind_omega_h_field2(py::module& m) {
+
+  // Bind OmegaHField2 class
+  py::class_<OmegaHField2, FieldT<Real>, std::shared_ptr<OmegaHField2>>(m, "OmegaHField2")
+    .def(py::init<const OmegaHFieldLayout&>(),
+         py::arg("layout"),
+         "Constructor for OmegaHField2")
+
+    .def("get_localization_hint", [](const OmegaHField2& self,
+                                     py::array_t<Real> coordinates,
+                                     const CoordinateSystem& coord_system) {
+      // Create CoordinateView from numpy array
+      auto coords_view = numpy_to_view_2d<const Real>(coordinates);
+      CoordinateView<HostMemorySpace> coord_view(coord_system, coords_view);
+      return self.GetLocalizationHint(coord_view);
+    },
+    py::arg("coordinates"),
+    py::arg("coordinate_system"),
+    "Get localization hint for given coordinates")
+
+    .def("evaluate", [](const OmegaHField2& self,
+                       const LocalizationHint& location,
+                       py::array_t<Real> values,
+                       const CoordinateSystem& coord_system) {
+      // Create FieldDataView
+      auto values_view = numpy_to_view<Real>(values);
+      FieldDataView<Real, HostMemorySpace> field_data_view(values_view, coord_system);
+      self.Evaluate(location, field_data_view);
+    },
+    py::arg("location"),
+    py::arg("values"),
+    py::arg("coordinate_system"),
+    "Evaluate field at locations specified by localization hint")
+
+    .def("evaluate_gradient", [](OmegaHField2& self,
+                                py::array_t<Real> gradients,
+                                const CoordinateSystem& coord_system) {
+      auto gradients_view = numpy_to_view<Real>(gradients);
+      FieldDataView<Real, HostMemorySpace> field_data_view(gradients_view, coord_system);
+      self.EvaluateGradient(field_data_view);
+    },
+    py::arg("gradients"),
+    py::arg("coordinate_system"),
+    "Evaluate gradient of the field")
+
+    .def("get_layout", &OmegaHField2::GetLayout,
+         py::return_value_policy::reference,
+         "Get the field layout")
+
+    .def("can_evaluate_gradient", &OmegaHField2::CanEvaluateGradient,
+         "Check if gradient evaluation is supported")
+
+    .def("serialize", [](const OmegaHField2& self,
+                        py::array_t<Real> buffer,
+                        py::array_t<const pcms::LO> permutation) {
+      auto buffer_view = numpy_to_view<Real>(buffer);
+      auto perm_view = numpy_to_view<const pcms::LO>(permutation);
+      return self.Serialize(buffer_view, perm_view);
+    },
+    py::arg("buffer"),
+    py::arg("permutation"),
+    "Serialize field data into buffer")
+
+    .def("deserialize", [](OmegaHField2& self,
+                          py::array_t<const Real> buffer,
+                          py::array_t<const pcms::LO> permutation) {
+      auto buffer_view = numpy_to_view<const Real>(buffer);
+      auto perm_view = numpy_to_view<const pcms::LO>(permutation);
+      self.Deserialize(buffer_view, perm_view);
+    },
+    py::arg("buffer"),
+    py::arg("permutation"),
+    "Deserialize field data from buffer")
+
+    .def("get_dof_holder_data", [](const OmegaHField2& self) {
+      auto const_data = self.GetDOFHolderData();
+      // Create a numpy array that owns its own data
+      return view_to_numpy<const Real>(const_data);
+    },
+    "Get the DOF holder data")
+
+    .def("set_dof_holder_data", [](OmegaHField2& self,
+                                   py::array_t<Real> data) {
+      // Ensure array is contiguous
+      auto contiguous_data = py::array_t<Real>(data);
+      auto data_view = numpy_to_view<Real>(contiguous_data);
+      // Create const view wrapper
+      Rank1View<const Real, HostMemorySpace> const_view(data_view.data_handle(), data_view.size());
+      self.SetDOFHolderData(const_view);
+    },
+    py::arg("data"),
+    "Set the DOF holder data");
+
+  // Helper functions for creating views (if needed for testing)
+  m.def("create_coordinate_view", [](py::array_t<Real> coordinates,
+                                     const CoordinateSystem& coord_system) {
+    auto coords_view = numpy_to_view_2d<const Real>(coordinates);
+    return CoordinateView<HostMemorySpace>(coord_system, coords_view);
+  },
+  py::arg("coordinates"),
+  py::arg("coordinate_system"),
+  "Create a CoordinateView from numpy array");
+
+  m.def("create_field_data_view", [](py::array_t<Real> values,
+                                     const CoordinateSystem& coord_system) {
+    auto values_view = numpy_to_view<Real>(values);
+    return FieldDataView<Real, HostMemorySpace>(values_view, coord_system);
+  },
+  py::arg("values"),
+  py::arg("coordinate_system"),
+  "Create a FieldDataView from numpy array");
+}
+
+} // namespace pcms
\ No newline at end of file
diff --git a/src/pcms/pythonapi/bind_omega_h_field_layout.cpp b/src/pcms/pythonapi/bind_omega_h_field_layout.cpp
new file mode 100644
index 00000000..e228a43e
--- /dev/null
+++ b/src/pcms/pythonapi/bind_omega_h_field_layout.cpp
@@ -0,0 +1,119 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/numpy.h>
+#include "pcms/adapter/omega_h/omega_h_field_layout.h"
+#include "pcms/field.h"
+#include "pcms/field_layout.h"
+#include "numpy_array_transform.h"
+
+namespace py = pybind11;
+
+namespace pcms {
+
+void bind_omega_h_field_layout_module(py::module& m) {
+  // Bind the OmegaHFieldLayout class
+  py::class_<OmegaHFieldLayout, FieldLayout, std::shared_ptr<OmegaHFieldLayout>>(m, "OmegaHFieldLayout")
+    .def(py::init<Omega_h::Mesh&, std::array<int, 4>, int, CoordinateSystem, std::string>(),
+         py::arg("mesh"),
+         py::arg("nodes_per_dim"),
+         py::arg("num_components"),
+         py::arg("coordinate_system"),
+         py::arg("global_id_name") = "global",
+         "Constructor for OmegaHFieldLayout")
+
+    .def("create_field", [](OmegaHFieldLayout& self) {
+           return std::shared_ptr<FieldT<Real>>(self.CreateField());
+         },
+         "Create a field with this layout")
+
+    .def("get_num_components", &OmegaHFieldLayout::GetNumComponents,
+         "Get the number of components in the field")
+
+    .def("get_num_owned_dof_holder", &OmegaHFieldLayout::GetNumOwnedDofHolder,
+         "Get the number of owned DOF holders")
+
+    .def("get_num_global_dof_holder", &OmegaHFieldLayout::GetNumGlobalDofHolder,
+         "Get the number of global DOF holders")
+
+    .def("get_owned", [](const OmegaHFieldLayout& self) {
+      auto owned = self.GetOwned();
+      // Convert to numpy array
+      py::array_t<bool> result(owned.extent(0));
+      auto buf = result.request();
+      bool* ptr = static_cast<bool*>(buf.ptr);
+      for (size_t i = 0; i < owned.extent(0); ++i) {
+        ptr[i] = owned[i];
+      }
+      return result;
+    },
+    "Get the owned mask array")
+
+    .def("get_gids", [](const OmegaHFieldLayout& self) {
+      auto gids = self.GetGids();
+      // Convert to numpy array
+      py::array_t<GO> result(gids.extent(0));
+      auto buf = result.request();
+      GO* ptr = static_cast<GO*>(buf.ptr);
+      for (size_t i = 0; i < gids.extent(0); ++i) {
+        ptr[i] = gids[i];
+      }
+      return result;
+    },
+    "Get the global IDs array")
+
+    .def("get_dof_holder_coordinates", [](const OmegaHFieldLayout& self) {
+      auto coords = self.GetDOFHolderCoordinates().GetCoordinates();
+      // Convert to numpy array (2D)
+      py::array_t<Real> result({coords.extent(0), coords.extent(1)});
+      auto buf = result.request();
+      Real* ptr = static_cast<Real*>(buf.ptr);
+      for (size_t i = 0; i < coords.extent(0); ++i) {
+        for (size_t j = 0; j < coords.extent(1); ++j) {
+          ptr[i * coords.extent(1) + j] = coords(i, j);
+        }
+      }
+      return result;
+    },
+    "Get the DOF holder coordinates")
+
+    .def("is_distributed", &OmegaHFieldLayout::IsDistributed,
+         "Check if the field layout is distributed")
+
+    .def("get_ent_offsets", [](const OmegaHFieldLayout& self) {
+      auto offsets = self.GetEntOffsets();
+      // Convert std::array to list/tuple
+      py::list result;
+      for (size_t i = 0; i < offsets.size(); ++i) {
+        result.append(offsets[i]);
+      }
+      return result;
+    },
+    "Get the entity offsets array")
+
+    .def("get_nodes_per_dim", [](const OmegaHFieldLayout& self) {
+      auto nodes = self.GetNodesPerDim();
+      py::list result;
+      for (size_t i = 0; i < nodes.size(); ++i) {
+        result.append(nodes[i]);
+      }
+      return result;
+    },
+    "Get the nodes per dimension array")
+
+    .def("get_num_ents", &OmegaHFieldLayout::GetNumEnts,
+         "Get the total number of entities")
+
+    .def("get_mesh", [](OmegaHFieldLayout& self) -> Omega_h::Mesh& {
+      return self.GetMesh();
+    },
+    py::return_value_policy::reference,
+    "Get the underlying Omega_h mesh")
+
+    .def("owned_size", &OmegaHFieldLayout::OwnedSize,
+         "Get the owned size (num_components * num_owned_dof_holder)")
+
+    .def("global_size", &OmegaHFieldLayout::GlobalSize,
+         "Get the global size (num_components * num_global_dof_holder)");
+}
+
+} // namespace pcms
\ No newline at end of file
diff --git a/src/pcms/pythonapi/bind_transfer_field2.cpp b/src/pcms/pythonapi/bind_transfer_field2.cpp
new file mode 100644
index 00000000..3f9c28b7
--- /dev/null
+++ b/src/pcms/pythonapi/bind_transfer_field2.cpp
@@ -0,0 +1,62 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include "pcms/transfer_field2.h"
+#include "pcms/field.h"
+
+namespace py = pybind11;
+
+namespace pcms {
+
+template<typename T>
+void bind_transfer_field2(py::module& m, const std::string& type_suffix) {
+  std::string copy_name = "copy_field2_" + type_suffix;
+  std::string interpolate_name = "interpolate_field2_" + type_suffix;
+
+  m.def(copy_name.c_str(),
+        [](const FieldT<T>& source, FieldT<T>& target) {
+          copy_field2(source, target);
+        },
+        py::arg("source"),
+        py::arg("target"),
+        "Copy field data from source to target. Source and target must be of the same type.");
+
+  m.def(interpolate_name.c_str(),
+        [](const FieldT<T>& source, FieldT<T>& target) {
+          interpolate_field2(source, target);
+        },
+        py::arg("source"),
+        py::arg("target"),
+        "Interpolate field from source to target. Coordinate systems must match.");
+}
+
+void bind_transfer_field2_module(py::module& m) {
+  // Bind for Real type (typically double)
+  bind_transfer_field2<Real>(m, "Real");
+
+  // If Real is not double, also bind double explicitly
+  if constexpr (!std::is_same_v<Real, double>) {
+    bind_transfer_field2<double>(m, "Double");
+  }
+
+  // Bind for float if needed
+  bind_transfer_field2<float>(m, "Float");
+
+  // Add convenience aliases for the most common case (Real)
+  m.def("interpolate_field",
+        [](const FieldT<Real>& source, FieldT<Real>& target) {
+          interpolate_field2(source, target);
+        },
+        py::arg("source"),
+        py::arg("target"),
+        "Interpolate field from source to target. Coordinate systems must match.");
+
+  m.def("copy_field",
+        [](const FieldT<Real>& source, FieldT<Real>& target) {
+          copy_field2(source, target);
+        },
+        py::arg("source"),
+        py::arg("target"),
+        "Copy field data from source to target. Source and target must be of the same type.");
+}
+
+} // namespace pcms
\ No newline at end of file
diff --git a/src/pcms/pythonapi/bind_uniform_grid_field.cpp b/src/pcms/pythonapi/bind_uniform_grid_field.cpp
new file mode 100644
index 00000000..3f9436f5
--- /dev/null
+++ b/src/pcms/pythonapi/bind_uniform_grid_field.cpp
@@ -0,0 +1,230 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/functional.h>
+#include <pybind11/numpy.h>
+#include "pcms/arrays.h"
+#include "pcms/adapter/uniform_grid/uniform_grid_field.h"
+#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h"
+#include "numpy_array_transform.h"
+
+namespace py = pybind11;
+
+namespace pcms {
+
+void bind_uniform_grid_field_module(py::module& m) {
+
+  // Bind UniformGridField class for 2D
+  py::class_<UniformGridField<2>, FieldT<Real>, std::shared_ptr<UniformGridField<2>>>(
+    m, "UniformGridField2D")
+    .def(py::init<const UniformGridFieldLayout<2>&>(),
+         py::arg("layout"),
+         "Constructor for UniformGridField2D")
+
+    .def("get_localization_hint", [](const UniformGridField<2>& self,
+                                     py::array_t<Real> coordinates,
+                                     const CoordinateSystem& coord_system) {
+      // Create CoordinateView from numpy array
+      auto coords_view = numpy_to_view_2d<const Real>(coordinates);
+      CoordinateView<HostMemorySpace> coord_view(coord_system, coords_view);
+      return self.GetLocalizationHint(coord_view);
+    },
+    py::arg("coordinates"),
+    py::arg("coordinate_system"),
+    "Get localization hint for given coordinates")
+
+    .def("evaluate", [](const UniformGridField<2>& self,
+                       const LocalizationHint& location,
+                       py::array_t<Real> values,
+                       const CoordinateSystem& coord_system) {
+      // Create FieldDataView
+      auto values_view = numpy_to_view<Real>(values);
+      FieldDataView<Real, HostMemorySpace> field_data_view(values_view, coord_system);
+      self.Evaluate(location, field_data_view);
+    },
+    py::arg("location"),
+    py::arg("values"),
+    py::arg("coordinate_system"),
+    "Evaluate field at locations specified by localization hint")
+
+    .def("evaluate_gradient", [](UniformGridField<2>& self,
+                                py::array_t<Real> gradients,
+                                const CoordinateSystem& coord_system) {
+      auto gradients_view = numpy_to_view<Real>(gradients);
+      FieldDataView<Real, HostMemorySpace> field_data_view(gradients_view, coord_system);
+      self.EvaluateGradient(field_data_view);
+    },
+    py::arg("gradients"),
+    py::arg("coordinate_system"),
+    "Evaluate gradient of the field")
+
+    .def("get_layout", &UniformGridField<2>::GetLayout,
+         py::return_value_policy::reference,
+         "Get the field layout")
+
+    .def("can_evaluate_gradient", &UniformGridField<2>::CanEvaluateGradient,
+         "Check if gradient evaluation is supported")
+
+    .def("serialize", [](const UniformGridField<2>& self,
+                        py::array_t<Real> buffer,
+                        py::array_t<const pcms::LO> permutation) {
+      auto buffer_view = numpy_to_view<Real>(buffer);
+      auto perm_view = numpy_to_view<const pcms::LO>(permutation);
+      return self.Serialize(buffer_view, perm_view);
+    },
+    py::arg("buffer"),
+    py::arg("permutation"),
+    "Serialize field data into buffer")
+
+    .def("deserialize", [](UniformGridField<2>& self,
+                          py::array_t<const Real> buffer,
+                          py::array_t<const pcms::LO> permutation) {
+      auto buffer_view = numpy_to_view<const Real>(buffer);
+      auto perm_view = numpy_to_view<const pcms::LO>(permutation);
+      self.Deserialize(buffer_view, perm_view);
+    },
+    py::arg("buffer"),
+    py::arg("permutation"),
+    "Deserialize field data from buffer")
+
+    .def("get_dof_holder_data", [](const UniformGridField<2>& self) {
+      auto const_data = self.GetDOFHolderData();
+      // Create a numpy array that owns its own data
+      return view_to_numpy<const Real>(const_data);
+    },
+    "Get the DOF holder data")
+
+    .def("set_dof_holder_data", [](UniformGridField<2>& self,
+                                   py::array_t<Real> data) {
+      // Ensure array is contiguous
+      auto contiguous_data = py::array_t<Real>(data);
+      auto data_view = numpy_to_view<Real>(contiguous_data);
+      // Create const view wrapper
+      Rank1View<const Real, HostMemorySpace> const_view(data_view.data_handle(), data_view.size());
+      self.SetDOFHolderData(const_view);
+    },
+    py::arg("data"),
+    "Set the DOF holder data")
+
+    .def("to_mdspan", [](const UniformGridField<2>& self) {
+      return mdspan_view_to_numpy(self.to_mdspan());
+    },
+    "Get the DOF holder data as a 2D numpy array (copy)");
+
+  // Bind UniformGridField class for 3D
+  py::class_<UniformGridField<3>, FieldT<Real>, std::shared_ptr<UniformGridField<3>>>(
+    m, "UniformGridField3D")
+    .def(py::init<const UniformGridFieldLayout<3>&>(),
+         py::arg("layout"),
+         "Constructor for UniformGridField3D")
+
+    .def("get_localization_hint", [](const UniformGridField<3>& self,
+                                     py::array_t<Real> coordinates,
+                                     const CoordinateSystem& coord_system) {
+      // Create CoordinateView from numpy array
+      auto coords_view = numpy_to_view_2d<const Real>(coordinates);
+      CoordinateView<HostMemorySpace> coord_view(coord_system, coords_view);
+      return self.GetLocalizationHint(coord_view);
+    },
+    py::arg("coordinates"),
+    py::arg("coordinate_system"),
+    "Get localization hint for given coordinates")
+
+    .def("evaluate", [](const UniformGridField<3>& self,
+                       const LocalizationHint& location,
+                       py::array_t<Real> values,
+                       const CoordinateSystem& coord_system) {
+      // Create FieldDataView
+      auto values_view = numpy_to_view<Real>(values);
+      FieldDataView<Real, HostMemorySpace> field_data_view(values_view, coord_system);
+      self.Evaluate(location, field_data_view);
+    },
+    py::arg("location"),
+    py::arg("values"),
+    py::arg("coordinate_system"),
+    "Evaluate field at locations specified by localization hint")
+
+    .def("evaluate_gradient", [](UniformGridField<3>& self,
+                                py::array_t<Real> gradients,
+                                const CoordinateSystem& coord_system) {
+      auto gradients_view = numpy_to_view<Real>(gradients);
+      FieldDataView<Real, HostMemorySpace> field_data_view(gradients_view, coord_system);
+      self.EvaluateGradient(field_data_view);
+    },
+    py::arg("gradients"),
+    py::arg("coordinate_system"),
+    "Evaluate gradient of the field")
+
+    .def("get_layout", &UniformGridField<3>::GetLayout,
+         py::return_value_policy::reference,
+         "Get the field layout")
+
+    .def("can_evaluate_gradient", &UniformGridField<3>::CanEvaluateGradient,
+         "Check if gradient evaluation is supported")
+
+    .def("serialize", [](const UniformGridField<3>& self,
+                        py::array_t<Real> buffer,
+                        py::array_t<const pcms::LO> permutation) {
+      auto buffer_view = numpy_to_view<Real>(buffer);
+      auto perm_view = numpy_to_view<const pcms::LO>(permutation);
+      return self.Serialize(buffer_view, perm_view);
+    },
+    py::arg("buffer"),
+    py::arg("permutation"),
+    "Serialize field data into buffer")
+
+    .def("deserialize", [](UniformGridField<3>& self,
+                          py::array_t<const Real> buffer,
+                          py::array_t<const pcms::LO> permutation) {
+      auto buffer_view = numpy_to_view<const Real>(buffer);
+      auto perm_view = numpy_to_view<const pcms::LO>(permutation);
+      self.Deserialize(buffer_view, perm_view);
+    },
+    py::arg("buffer"),
+    py::arg("permutation"),
+    "Deserialize field data from buffer")
+
+    .def("get_dof_holder_data", [](const UniformGridField<3>& self) {
+      auto const_data = self.GetDOFHolderData();
+      // Create a numpy array that owns its own data
+      return view_to_numpy<const Real>(const_data);
+    },
+    "Get the DOF holder data")
+
+    .def("set_dof_holder_data", [](UniformGridField<3>& self,
+                                   py::array_t<Real> data) {
+      // Ensure array is contiguous
+      auto contiguous_data = py::array_t<Real>(data);
+      auto data_view = numpy_to_view<Real>(contiguous_data);
+      // Create const view wrapper
+      Rank1View<const Real, HostMemorySpace> const_view(data_view.data_handle(), data_view.size());
+      self.SetDOFHolderData(const_view);
+    },
+    py::arg("data"),
+    "Set the DOF holder data")
+
+    .def("to_mdspan", [](const UniformGridField<3>& self) {
+      return mdspan_view_to_numpy(self.to_mdspan());
+    },
+    "Get the DOF holder data as a 3D numpy array (copy)");
+
+  // Helper functions for creating views (if needed for testing)
+  m.def("create_coordinate_view", [](py::array_t<Real> coordinates,
+                                     const CoordinateSystem& coord_system) {
+    auto coords_view = numpy_to_view_2d<const Real>(coordinates);
+    return CoordinateView<HostMemorySpace>(coord_system, coords_view);
+  },
+  py::arg("coordinates"),
+  py::arg("coordinate_system"),
+  "Create a CoordinateView from numpy array");
+
+  m.def("create_field_data_view", [](py::array_t<Real> values,
+                                     const CoordinateSystem& coord_system) {
+    auto values_view = numpy_to_view<Real>(values);
+    return FieldDataView<Real, HostMemorySpace>(values_view, coord_system);
+  },
+  py::arg("values"),
+  py::arg("coordinate_system"),
+  "Create a FieldDataView from numpy array");
+}
+
+} // namespace pcms
diff --git a/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp b/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp
new file mode 100644
index 00000000..2d8329ad
--- /dev/null
+++ b/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp
@@ -0,0 +1,222 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/numpy.h>
+#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h"
+#include "pcms/field.h"
+#include "pcms/field_layout.h"
+#include "pcms/uniform_grid.h"
+#include "numpy_array_transform.h"
+
+namespace py = pybind11;
+
+namespace pcms {
+
+void bind_uniform_grid_field_layout_module(py::module& m) {
+  // Bind UniformGrid structure for 2D
+  py::class_<UniformGrid<2>>(m, "UniformGrid2D")
+    .def(py::init<>())
+    .def_readwrite("edge_length", &UniformGrid<2>::edge_length,
+                   "Edge length in each dimension")
+    .def_readwrite("bot_left", &UniformGrid<2>::bot_left,
+                   "Bottom-left corner coordinates")
+    .def_readwrite("divisions", &UniformGrid<2>::divisions,
+                   "Number of divisions in each dimension")
+    .def("get_num_cells", &UniformGrid<2>::GetNumCells,
+         "Get the total number of cells in the grid")
+    .def("closest_cell_id", [](const UniformGrid<2>& self, py::array_t<Real> point) {
+           if (point.size() != 2) {
+             throw std::runtime_error("Point must have 2 coordinates for 2D grid");
+           }
+           auto buf = point.request();
+           Real* ptr = static_cast<Real*>(buf.ptr);
+           Omega_h::Vector<2> omega_point({ptr[0], ptr[1]});
+           return self.ClosestCellID(omega_point);
+         },
+         py::arg("point"),
+         "Get the cell ID that contains or is closest to the given point")
+    .def("get_cell_bbox", &UniformGrid<2>::GetCellBBOX,
+         py::arg("cell_index"),
+         "Get the bounding box of a cell");
+
+  // Bind UniformGrid structure for 3D
+  py::class_<UniformGrid<3>>(m, "UniformGrid3D")
+    .def(py::init<>())
+    .def_readwrite("edge_length", &UniformGrid<3>::edge_length,
+                   "Edge length in each dimension")
+    .def_readwrite("bot_left", &UniformGrid<3>::bot_left,
+                   "Bottom-left corner coordinates")
+    .def_readwrite("divisions", &UniformGrid<3>::divisions,
+                   "Number of divisions in each dimension")
+    .def("get_num_cells", &UniformGrid<3>::GetNumCells,
+         "Get the total number of cells in the grid")
+    .def("closest_cell_id", [](const UniformGrid<3>& self, py::array_t<Real> point) {
+           if (point.size() != 3) {
+             throw std::runtime_error("Point must have 3 coordinates for 3D grid");
+           }
+           auto buf = point.request();
+           Real* ptr = static_cast<Real*>(buf.ptr);
+           Omega_h::Vector<3> omega_point({ptr[0], ptr[1], ptr[2]});
+           return self.ClosestCellID(omega_point);
+         },
+         py::arg("point"),
+         "Get the cell ID that contains or is closest to the given point")
+    .def("get_cell_bbox", &UniformGrid<3>::GetCellBBOX,
+         py::arg("cell_index"),
+         "Get the bounding box of a cell");
+
+  // Bind the UniformGridFieldLayout class for 2D
+  py::class_<UniformGridFieldLayout<2>, FieldLayout, std::shared_ptr<UniformGridFieldLayout<2>>>(
+    m, "UniformGridFieldLayout2D")
+    .def(py::init<UniformGrid<2>&, int, CoordinateSystem>(),
+         py::arg("grid"),
+         py::arg("num_components"),
+         py::arg("coordinate_system"),
+         "Constructor for UniformGridFieldLayout2D")
+
+    .def("create_field", [](UniformGridFieldLayout<2>& self) {
+           return std::shared_ptr<FieldT<Real>>(self.CreateField());
+         },
+         "Create a field with this layout")
+
+    .def("get_num_components", &UniformGridFieldLayout<2>::GetNumComponents,
+         "Get the number of components in the field")
+
+    .def("get_num_owned_dof_holder", &UniformGridFieldLayout<2>::GetNumOwnedDofHolder,
+         "Get the number of owned DOF holders")
+
+    .def("get_num_global_dof_holder", &UniformGridFieldLayout<2>::GetNumGlobalDofHolder,
+         "Get the number of global DOF holders")
+
+    .def("get_owned", [](const UniformGridFieldLayout<2>& self) {
+      auto owned = self.GetOwned();
+      // Convert to numpy array
+      py::array_t<bool> result(owned.extent(0));
+      auto buf = result.request();
+      bool* ptr = static_cast<bool*>(buf.ptr);
+      for (size_t i = 0; i < owned.extent(0); ++i) {
+        ptr[i] = owned[i];
+      }
+      return result;
+    },
+    "Get the owned mask array")
+
+    .def("get_gids", [](const UniformGridFieldLayout<2>& self) {
+      auto gids = self.GetGids();
+      // Convert to numpy array
+      py::array_t<GO> result(gids.extent(0));
+      auto buf = result.request();
+      GO* ptr = static_cast<GO*>(buf.ptr);
+      for (size_t i = 0; i < gids.extent(0); ++i) {
+        ptr[i] = gids[i];
+      }
+      return result;
+    },
+    "Get the global IDs array")
+
+    .def("get_dof_holder_coordinates", [](const UniformGridFieldLayout<2>& self) {
+      auto coords = self.GetDOFHolderCoordinates().GetCoordinates();
+      // Convert to numpy array (2D)
+      py::array_t<Real> result({coords.extent(0), coords.extent(1)});
+      auto buf = result.request();
+      Real* ptr = static_cast<Real*>(buf.ptr);
+      for (size_t i = 0; i < coords.extent(0); ++i) {
+        for (size_t j = 0; j < coords.extent(1); ++j) {
+          ptr[i * coords.extent(1) + j] = coords(i, j);
+        }
+      }
+      return result;
+    },
+    "Get the DOF holder coordinates")
+
+    .def("is_distributed", &UniformGridFieldLayout<2>::IsDistributed,
+         "Check if the field layout is distributed")
+
+    .def("get_grid", &UniformGridFieldLayout<2>::GetGrid,
+         py::return_value_policy::reference,
+         "Get the underlying uniform grid")
+
+    .def("get_num_cells", &UniformGridFieldLayout<2>::GetNumCells,
+         "Get the number of cells in the grid")
+
+    .def("get_num_vertices", &UniformGridFieldLayout<2>::GetNumVertices,
+         "Get the number of vertices in the grid");
+
+  // Bind the UniformGridFieldLayout class for 3D
+  py::class_<UniformGridFieldLayout<3>, FieldLayout, std::shared_ptr<UniformGridFieldLayout<3>>>(
+    m, "UniformGridFieldLayout3D")
+    .def(py::init<UniformGrid<3>&, int, CoordinateSystem>(),
+         py::arg("grid"),
+         py::arg("num_components"),
+         py::arg("coordinate_system"),
+         "Constructor for UniformGridFieldLayout3D")
+
+    .def("create_field", [](UniformGridFieldLayout<3>& self) {
+           return std::shared_ptr<FieldT<Real>>(self.CreateField());
+         },
+         "Create a field with this layout")
+
+    .def("get_num_components", &UniformGridFieldLayout<3>::GetNumComponents,
+         "Get the number of components in the field")
+
+    .def("get_num_owned_dof_holder", &UniformGridFieldLayout<3>::GetNumOwnedDofHolder,
+         "Get the number of owned DOF holders")
+
+    .def("get_num_global_dof_holder", &UniformGridFieldLayout<3>::GetNumGlobalDofHolder,
+         "Get the number of global DOF holders")
+
+    .def("get_owned", [](const UniformGridFieldLayout<3>& self) {
+      auto owned = self.GetOwned();
+      // Convert to numpy array
+      py::array_t<bool> result(owned.extent(0));
+      auto buf = result.request();
+      bool* ptr = static_cast<bool*>(buf.ptr);
+      for (size_t i = 0; i < owned.extent(0); ++i) {
+        ptr[i] = owned[i];
+      }
+      return result;
+    },
+    "Get the owned mask array")
+
+    .def("get_gids", [](const UniformGridFieldLayout<3>& self) {
+      auto gids = self.GetGids();
+      // Convert to numpy array
+      py::array_t<GO> result(gids.extent(0));
+      auto buf = result.request();
+      GO* ptr = static_cast<GO*>(buf.ptr);
+      for (size_t i = 0; i < gids.extent(0); ++i) {
+        ptr[i] = gids[i];
+      }
+      return result;
+    },
+    "Get the global IDs array")
+
+    .def("get_dof_holder_coordinates", [](const UniformGridFieldLayout<3>& self) {
+      auto coords = self.GetDOFHolderCoordinates().GetCoordinates();
+      // Convert to numpy array (2D)
+      py::array_t<Real> result({coords.extent(0), coords.extent(1)});
+      auto buf = result.request();
+      Real* ptr = static_cast<Real*>(buf.ptr);
+      for (size_t i = 0; i < coords.extent(0); ++i) {
+        for (size_t j = 0; j < coords.extent(1); ++j) {
+          ptr[i * coords.extent(1) + j] = coords(i, j);
+        }
+      }
+      return result;
+    },
+    "Get the DOF holder coordinates")
+
+    .def("is_distributed", &UniformGridFieldLayout<3>::IsDistributed,
+         "Check if the field layout is distributed")
+
+    .def("get_grid", &UniformGridFieldLayout<3>::GetGrid,
+         py::return_value_policy::reference,
+         "Get the underlying uniform grid")
+
+    .def("get_num_cells", &UniformGridFieldLayout<3>::GetNumCells,
+         "Get the number of cells in the grid")
+
+    .def("get_num_vertices", &UniformGridFieldLayout<3>::GetNumVertices,
+         "Get the number of vertices in the grid");
+}
+
+} // namespace pcms
diff --git a/src/pcms/pythonapi/numpy_array_transform.h b/src/pcms/pythonapi/numpy_array_transform.h
new file mode 100644
index 00000000..b1d6203d
--- /dev/null
+++ b/src/pcms/pythonapi/numpy_array_transform.h
@@ -0,0 +1,232 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/numpy.h>
+#include <Omega_h_array.hpp>
+#include "pcms/arrays.h"
+#include "pcms/types.h"
+
+
+namespace py = pybind11;
+
+
+namespace pcms {
+
+
+// Helper function to convert numpy array to Rank1 View
+template<typename T>
+Rank1View<T, HostMemorySpace> numpy_to_view(py::array_t<T> arr) {
+  py::buffer_info buf = arr.request();
+  if (buf.ndim != 1) {
+    throw std::runtime_error("Number of dimensions must be 1");
+  }
+  T* data_ptr = reinterpret_cast<T*>(buf.ptr);
+  Rank1View<T, HostMemorySpace> view(data_ptr, buf.shape[0]);
+  return view;
+}
+
+
+// Helper function to convert numpy array to Kokkos View
+template<typename T>
+Kokkos::View<T*, HostMemorySpace> numpy_to_kokkos_view(py::array_t<T> arr) {
+  py::buffer_info buf = arr.request();
+  if (buf.ndim != 1) {
+    throw std::runtime_error("Number of dimensions must be 1");
+  }
+  Kokkos::View<T*, HostMemorySpace> view(
+    reinterpret_cast<T*>(buf.ptr), buf.shape[0]);
+  return view;
+}
+
+
+// Helper function to convert Rank1 View to array (creates a copy)
+template<typename T>
+py::array_t<typename std::remove_const<T>::type> view_to_numpy(Rank1View<T, HostMemorySpace> view) {
+  // Create a copy to avoid lifetime issues and mdspan type registration
+  using NonConstT = typename std::remove_const<T>::type;
+  py::array_t<NonConstT> result(view.size());
+  py::buffer_info buf = result.request();
+  NonConstT* ptr = static_cast<NonConstT*>(buf.ptr);
+  // Use element-wise copy to ensure proper mdspan access
+  for (size_t i = 0; i < view.size(); ++i) {
+    ptr[i] = view[i];
+  }
+  return result;
+}
+
+
+// Helper function to convert Kokkos View to array (creates a reference)
+template<typename T>
+py::array_t<T> kokkos_view_to_numpy(Kokkos::View<T*, HostMemorySpace> view) {
+  return py::array_t<T>(
+    {static_cast<py::ssize_t>(view.extent(0))},  // shape
+    {sizeof(T)},                                   // strides
+    view.data(),                                   // data pointer
+    py::cast(view)                                 // base object to manage lifetime
+  );
+}
+
+
+// Helper to convert 2D numpy array to Kokkos View
+template<typename T>
+Kokkos::View<T**, HostMemorySpace> numpy_to_kokkos_view_2d(py::array_t<T> arr) {
+  py::buffer_info buf = arr.request();
+  if (buf.ndim != 2) {
+    throw std::runtime_error("Number of dimensions must be 2");
+  }
+  Kokkos::View<T**, HostMemorySpace> view(
+    reinterpret_cast<T*>(buf.ptr), buf.shape[0], buf.shape[1]);
+  return view;
+}
+
+
+// Helper to convert 2D Kokkos View to numpy array (creates a reference)
+template<typename T>
+py::array_t<T> kokkos_view_2d_to_numpy(Kokkos::View<T**, HostMemorySpace> view) {
+  return py::array_t<T>(
+    {static_cast<py::ssize_t>(view.extent(0)), static_cast<py::ssize_t>(view.extent(1))},  // shape
+    {sizeof(T) * view.extent(1), sizeof(T)},                                                 // strides (row-major)
+    view.data(),                                                                             // data pointer
+    py::cast(view)                                                                           // base object to manage lifetime
+  );
+}
+
+
+// Helper to convert 1D numpy array to Omega_h::Read
+template<typename T>
+Omega_h::Read<T> numpy_to_omega_h_read(py::array_t<T> arr) {
+  py::buffer_info buf = arr.request();
+  if (buf.ndim != 1) {
+    throw std::runtime_error("Number of dimensions must be 1");
+  }
+  Kokkos::View<T*, Kokkos::DefaultExecutionSpace::memory_space,
+    Kokkos::MemoryTraits<Kokkos::Unmanaged>> view(
+    reinterpret_cast<T*>(buf.ptr), buf.shape[0]);
+  Omega_h::Write<T> write_view(view);
+  Omega_h::Read<T> read_view(write_view);
+  return read_view;
+}
+
+
+// Helper to convert Omega_h::Read to numpy array (creates a copy)
+template<typename T>
+py::array_t<T> omega_h_read_to_numpy(Omega_h::Read<T> read_view) {
+  auto read_view_host = Omega_h::HostRead<T>(read_view);
+  py::array_t<T> result(read_view_host.size());
+  py::buffer_info buf = result.request();
+  T* ptr = static_cast<T*>(buf.ptr);
+  for (Omega_h::LO i = 0; i < read_view_host.size(); ++i) {
+    ptr[i] = read_view_host[i];
+  }
+  return result;
+}
+
+
+// Helper to convert 1D numpy array to Omega_h::Write
+template<typename T>
+Omega_h::Write<T> numpy_to_omega_h_write(py::array_t<T> arr) {
+  py::buffer_info buf = arr.request();
+  if (buf.ndim != 1) {
+    throw std::runtime_error("Number of dimensions must be 1");
+  }
+  // Get host mirror and copy data
+  auto write_view_host = Omega_h::HostWrite<T>(buf.shape[0]);
+  T* ptr = static_cast<T*>(buf.ptr);
+  for (Omega_h::LO i = 0; i < buf.shape[0]; ++i) {
+    write_view_host[i] = ptr[i];
+  }
+  auto write_view = Omega_h::Write<T>(write_view_host);
+  return write_view;
+}
+
+
+// Helper to convert Omega_h::Write to numpy array (creates a reference)
+template<typename T>
+py::array_t<T> omega_h_write_to_numpy(Omega_h::Write<T> write_view) {
+  return py::array_t<T>(
+    {static_cast<py::ssize_t>(write_view.size())},  // shape
+    {sizeof(T)},                                     // strides
+    write_view.data(),                               // data pointer
+    py::cast(write_view)                             // base object to manage lifetime
+  );
+}
+
+
+// Helper function to convert numpy array to Rank2 View
+template<typename T>
+Rank2View<T, HostMemorySpace> numpy_to_view_2d(py::array_t<T> arr) {
+  py::buffer_info buf = arr.request();
+  if (buf.ndim != 2) {
+    throw std::runtime_error("Number of dimensions must be 2");
+  }
+  Rank2View<T, HostMemorySpace> view(
+    reinterpret_cast<T*>(buf.ptr), buf.shape[0], buf.shape[1]);
+  return view;
+}
+
+
+// Helper function to convert Rank2 View to array (creates a reference)
+template<typename T>
+py::array_t<T> view_2d_to_numpy(Rank2View<T, HostMemorySpace> view) {
+  return py::array_t<T>(
+    {static_cast<py::ssize_t>(view.extent(0)), static_cast<py::ssize_t>(view.extent(1))},  // shape
+    {sizeof(T) * view.extent(1), sizeof(T)},                                                 // strides (row-major)
+    view.data_handle(),                                                                      // data pointer
+    py::cast(view)                                                                           // base object to manage lifetime
+  );
+}
+
+
+// Helper function to convert Rank3 View to array (creates a reference)
+template<typename T>
+py::array_t<T> view_3d_to_numpy(Rank3View<T, HostMemorySpace> view) {
+  return py::array_t<T>(
+    {static_cast<py::ssize_t>(view.extent(0)),
+     static_cast<py::ssize_t>(view.extent(1)),
+     static_cast<py::ssize_t>(view.extent(2))},  // shape
+    {sizeof(T) * view.extent(1) * view.extent(2),
+     sizeof(T) * view.extent(2),
+     sizeof(T)},                                 // strides (row-major)
+    view.data_handle(),                          // data pointer
+    py::cast(view)                               // base object to manage lifetime
+  );
+}
+
+
+// Helper function to convert mdspan View to array (creates a reference)
+template<typename T>
+py::array_t<typename std::remove_const<T>::type> mdspan_view_to_numpy(
+  View<2, T, HostMemorySpace> view) {
+  using NonConstT = typename std::remove_const<T>::type;
+  py::array_t<NonConstT> result(
+    {static_cast<py::ssize_t>(view.extent(0)),
+     static_cast<py::ssize_t>(view.extent(1))});
+  auto out = result.template mutable_unchecked<2>();
+  for (py::ssize_t i = 0; i < out.shape(0); ++i) {
+    for (py::ssize_t j = 0; j < out.shape(1); ++j) {
+      out(i, j) = view(i, j);
+    }
+  }
+  return result;
+}
+
+template<typename T>
+py::array_t<typename std::remove_const<T>::type> mdspan_view_to_numpy(
+  View<3, T, HostMemorySpace> view) {
+  using NonConstT = typename std::remove_const<T>::type;
+  py::array_t<NonConstT> result(
+    {static_cast<py::ssize_t>(view.extent(0)),
+     static_cast<py::ssize_t>(view.extent(1)),
+     static_cast<py::ssize_t>(view.extent(2))});
+  auto out = result.template mutable_unchecked<3>();
+  for (py::ssize_t i = 0; i < out.shape(0); ++i) {
+    for (py::ssize_t j = 0; j < out.shape(1); ++j) {
+      for (py::ssize_t k = 0; k < out.shape(2); ++k) {
+        out(i, j, k) = view(i, j, k);
+      }
+    }
+  }
+  return result;
+}
+
+
+} // namespace pcms
\ No newline at end of file
diff --git a/src/pcms/pythonapi/pythonapi.cpp b/src/pcms/pythonapi/pythonapi.cpp
new file mode 100644
index 00000000..7fb274ed
--- /dev/null
+++ b/src/pcms/pythonapi/pythonapi.cpp
@@ -0,0 +1,53 @@
+
+#include <pybind11/pybind11.h>
+
+namespace py = pybind11;
+
+namespace pcms {
+
+void bind_omega_h_field2(py::module& m);
+
+void bind_transfer_field2_module(py::module& m);
+
+void bind_omega_h_mesh_module(py::module& m);
+
+// void bind_interpolation_base_module(py::module& m);
+
+void bind_coordinate_system_module(py::module& m);
+
+void bind_coordinate_module(py::module& m);
+
+void bind_create_field_module(py::module& m);
+
+void bind_field_layout_module(py::module& m);
+
+void bind_field_module(py::module& m);
+
+void bind_omega_h_field_layout_module(py::module& m);
+
+void bind_uniform_grid_field_layout_module(py::module& m);
+
+void bind_uniform_grid_field_module(py::module& m);
+
+}
+PYBIND11_MODULE(py_pcms, m) {
+  // Bind fundamental types first (coordinate systems, etc.)
+  pcms::bind_coordinate_system_module(m);
+  pcms::bind_coordinate_module(m);
+
+  // Bind mesh and field infrastructure
+  pcms::bind_omega_h_mesh_module(m);
+  pcms::bind_field_layout_module(m);
+  pcms::bind_field_module(m);
+  pcms::bind_omega_h_field_layout_module(m);
+  pcms::bind_uniform_grid_field_layout_module(m);
+  pcms::bind_create_field_module(m);
+
+  // Bind field types
+  pcms::bind_omega_h_field2(m);
+  pcms::bind_uniform_grid_field_module(m);
+
+  // Bind field operations
+  pcms::bind_transfer_field2_module(m);
+  // pcms::bind_interpolation_base_module(m);
+}
\ No newline at end of file
diff --git a/src/pcms/pythonapi/test_field_copy.py b/src/pcms/pythonapi/test_field_copy.py
new file mode 100644
index 00000000..36d2ee23
--- /dev/null
+++ b/src/pcms/pythonapi/test_field_copy.py
@@ -0,0 +1,84 @@
+import py_pcms
+import numpy as np
+
+def test_copy(world, dim, order, num_components):
+    """Test copying omega_h_field2 data"""
+    nx = 100
+    ny = 100 if dim > 1 else 0
+    nz = 100 if dim > 2 else 0
+    print(f"\nStarting test: dim={dim}, order={order}, num_components={num_components}")
+
+    # Build mesh
+    mesh = py_pcms.build_box(
+        world, 
+        py_pcms.Family.SIMPLEX, 
+        1.0, 1.0, 1.0, 
+        nx, ny, nz, 
+        False
+    )
+    print(f"  Built {dim}D mesh with {nx}x{ny}x{nz} elements")
+    print(f"  Mesh type: {type(mesh)}")
+    print(f"  Mesh object: {mesh}")
+
+    # Create layout
+    print("  About to create layout...")
+    layout = py_pcms.create_lagrange_layout(
+        mesh, 
+        order, 
+        num_components,
+        py_pcms.CoordinateSystem.Cartesian
+    )
+    print(f"  Layout created successfully")
+    print(f"Testing dim={dim}, order={order}, num_components={num_components}...")
+
+    # Get number of data points
+    ndata = layout.get_num_owned_dof_holder() * num_components
+    print(f"  Number of data points: {ndata}")
+
+    # Create sequential array of IDs
+    ids = np.arange(ndata, dtype=np.float64)
+    print(f"  Created array of IDs from 0 to {ndata-1}")
+
+    # Create original field and set data
+    original = layout.create_field()
+    print("  Created original field")
+    original.set_dof_holder_data(ids)
+    print("  Set data in original field")
+
+    # Create copied field and copy data
+    copied = layout.create_field()
+    py_pcms.copy_field2_Real(original, copied)
+    print("  Copied data to new field")
+
+    # Get copied data
+    copied_array = copied.get_dof_holder_data()
+
+    # Verify the copy
+    assert len(copied_array) == ndata, f"Expected {ndata} elements, got {len(copied_array)}"
+
+    # Check that all values match
+    differences = np.abs(ids - copied_array)
+    num_matches = np.sum(differences < 1e-12)
+
+    assert num_matches == ndata, f"Only {num_matches}/{ndata} elements matched"
+
+    print(f"✓ Test passed: dim={dim}, order={order}, num_components={num_components}")
+
+def main():
+    """Run all test cases"""
+    print("Testing copy omega_h_field2 data...")
+
+    # Initialize Omega_h library
+    lib = py_pcms.OmegaHLibrary()
+    world = lib.world()
+    print("Initialized Omega_h library and world")
+
+    # Run test cases
+    test_copy(world, 2, 1, 1)
+    print("first passed")
+    test_copy(world, 2, 2, 1)
+
+    print("\nAll tests passed!")
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file
diff --git a/src/pcms/pythonapi/test_file_io.py b/src/pcms/pythonapi/test_file_io.py
new file mode 100644
index 00000000..a688e7a4
--- /dev/null
+++ b/src/pcms/pythonapi/test_file_io.py
@@ -0,0 +1,452 @@
+#!/usr/bin/env python3
+"""
+Test Omega_h file I/O Python bindings
+"""
+
+import py_pcms
+import numpy as np
+import os
+import shutil
+
+def test_binary_io(lib, world):
+    """Test binary file format I/O"""
+    print("\n=== Testing Binary Format I/O ===")
+    
+    # Build a simple 3D box mesh
+    print("Creating test mesh...")
+    mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        1.0, 1.0, 1.0,  # x, y, z dimensions
+        4, 4, 4,         # nx, ny, nz divisions
+        False            # symmetric
+    )
+    
+    # Get initial mesh properties
+    nverts_orig = mesh.nverts()
+    nelems_orig = mesh.nelems()
+    dim_orig = mesh.dim()
+    print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}")
+    
+    # Create temporary directory for test files
+    test_dir = "test_data"
+    os.makedirs(test_dir, exist_ok=True)
+    try:
+        # Test binary write/read
+        binary_file = os.path.join(test_dir, "test_mesh.osh")
+        print(f"Writing to binary file: {binary_file}")
+        py_pcms.write_mesh_binary(binary_file, mesh)
+        
+        print(f"Reading from binary file: {binary_file}")
+        mesh_read = py_pcms.read_mesh_binary(binary_file, lib)
+        
+        # Verify mesh properties
+        assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}"
+        assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}"
+        assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}"
+        
+        print("✓ Binary I/O test passed")
+        
+    finally:
+        # Clean up temporary files
+        if os.path.exists(test_dir):
+            for item in os.listdir(test_dir):
+                item_path = os.path.join(test_dir, item)
+                if os.path.isfile(item_path):
+                    os.remove(item_path)
+                elif os.path.isdir(item_path):
+                    shutil.rmtree(item_path)
+            os.rmdir(test_dir)
+            print(f"Cleaned up test directory: {test_dir}")
+
+
+def test_gmsh_io(lib, world):
+    """Test Gmsh file format I/O"""
+    print("\n=== Testing Gmsh Format I/O ===")
+    
+    # Build a simple 2D box mesh
+    print("Creating test mesh...")
+    mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        2.0, 1.0, 0.0,  # x, y, z dimensions (z=0 for 2D)
+        5, 3, 0,         # nx, ny, nz divisions (nz=0 for 2D)
+        False            # symmetric
+    )
+    
+    # Get initial mesh properties
+    nverts_orig = mesh.nverts()
+    nelems_orig = mesh.nelems()
+    dim_orig = mesh.dim()
+    print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}")
+    
+    # Create temporary directory for test files
+    test_dir = "test_data"
+    os.makedirs(test_dir, exist_ok=True)
+    try:
+        # Test Gmsh write/read
+        gmsh_file = os.path.join(test_dir, "test_mesh.msh")
+        print(f"Writing to Gmsh file: {gmsh_file}")
+        py_pcms.write_mesh_gmsh(gmsh_file, mesh)
+        
+        print(f"Reading from Gmsh file: {gmsh_file}")
+        mesh_read = py_pcms.read_mesh_gmsh(gmsh_file, world)
+        
+        # Verify mesh properties
+        assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}"
+        assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}"
+        assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}"
+        
+        print("✓ Gmsh I/O test passed")
+        
+    finally:
+        # Clean up temporary files
+        if os.path.exists(test_dir):
+            for item in os.listdir(test_dir):
+                item_path = os.path.join(test_dir, item)
+                if os.path.isfile(item_path):
+                    os.remove(item_path)
+                elif os.path.isdir(item_path):
+                    shutil.rmtree(item_path)
+            os.rmdir(test_dir)
+            print(f"Cleaned up test directory: {test_dir}")
+
+
+def test_vtk_io(lib, world):
+    """Test VTK file format I/O"""
+    print("\n=== Testing VTK Format I/O ===")
+    
+    # Build a simple 3D box mesh
+    print("Creating test mesh...")
+    mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        1.5, 1.0, 0.5,  # x, y, z dimensions
+        3, 3, 2,         # nx, ny, nz divisions
+        False            # symmetric
+    )
+    
+    # Get initial mesh properties
+    nverts_orig = mesh.nverts()
+    nelems_orig = mesh.nelems()
+    dim_orig = mesh.dim()
+    coords_orig = mesh.coords()
+    print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}")
+    print(f"Coordinates shape: {coords_orig.shape}")
+    
+    # Create temporary directory for test files
+    test_dir = "test_data"
+    os.makedirs(test_dir, exist_ok=True)
+    try:
+        # Test VTU write (VTU is write-only in the API, typically for visualization)
+        vtu_file = os.path.join(test_dir, "test_mesh.vtu")
+        print(f"Writing to VTU file: {vtu_file}")
+        py_pcms.write_mesh_vtu(vtu_file, mesh, compress=False)
+        
+        # Verify file was created
+        assert os.path.exists(vtu_file), f"VTU file was not created: {vtu_file}"
+        file_size = os.path.getsize(vtu_file)
+        print(f"VTU file created successfully (size: {file_size} bytes)")
+        
+        # Test compressed VTU write
+        vtu_compressed = os.path.join(test_dir, "test_mesh_compressed.vtu")
+        print(f"Writing compressed VTU file: {vtu_compressed}")
+        py_pcms.write_mesh_vtu(vtu_compressed, mesh, compress=True)
+        
+        assert os.path.exists(vtu_compressed), f"Compressed VTU file was not created"
+        compressed_size = os.path.getsize(vtu_compressed)
+        print(f"Compressed VTU file created (size: {compressed_size} bytes)")
+        
+        print("✓ VTK I/O test passed")
+        
+    finally:
+        # Clean up temporary files
+        if os.path.exists(test_dir):
+            for item in os.listdir(test_dir):
+                item_path = os.path.join(test_dir, item)
+                if os.path.isfile(item_path):
+                    os.remove(item_path)
+                elif os.path.isdir(item_path):
+                    shutil.rmtree(item_path)
+            os.rmdir(test_dir)
+            print(f"Cleaned up test directory: {test_dir}")
+
+
+def test_meshb_io(lib, world):
+    """Test MESHB file format I/O"""
+    print("\n=== Testing MESHB Format I/O ===")
+    
+    # Check if MESHB support is available
+    if not hasattr(py_pcms, 'write_mesh_meshb'):
+        print("⊘ MESHB support not available (OMEGA_H_USE_LIBMESHB not enabled)")
+        return
+    
+    # Build a simple 3D box mesh
+    print("Creating test mesh...")
+    mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        1.0, 1.0, 1.0,  # x, y, z dimensions
+        3, 3, 3,         # nx, ny, nz divisions
+        False            # symmetric
+    )
+    
+    # Get initial mesh properties
+    nverts_orig = mesh.nverts()
+    nelems_orig = mesh.nelems()
+    dim_orig = mesh.dim()
+    print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}")
+    
+    # Create temporary directory for test files
+    test_dir = "test_data"
+    os.makedirs(test_dir, exist_ok=True)
+    try:
+        # Test MESHB write/read
+        meshb_file = os.path.join(test_dir, "test_mesh.mesh")
+        print(f"Writing to MESHB file: {meshb_file}")
+        py_pcms.write_mesh_meshb(mesh, meshb_file, version=2)
+        
+        print(f"Reading from MESHB file: {meshb_file}")
+        # MESHB read requires pre-created mesh object
+        mesh_read = py_pcms.OmegaHMesh(lib)
+        mesh_read.set_comm(world)
+        py_pcms.read_mesh_meshb(mesh_read, meshb_file)
+        
+        # Verify mesh properties
+        assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}"
+        assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}"
+        assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}"
+        
+        print("✓ MESHB I/O test passed")
+        
+    finally:
+        # Clean up temporary files
+        if os.path.exists(test_dir):
+            for item in os.listdir(test_dir):
+                item_path = os.path.join(test_dir, item)
+                if os.path.isfile(item_path):
+                    os.remove(item_path)
+                elif os.path.isdir(item_path):
+                    shutil.rmtree(item_path)
+            os.rmdir(test_dir)
+            print(f"Cleaned up test directory: {test_dir}")
+
+
+def test_exodus_io(lib, world):
+    """Test Exodus file format I/O"""
+    print("\n=== Testing Exodus Format I/O ===")
+    
+    # Check if Exodus support is available
+    if not hasattr(py_pcms, 'write_mesh_exodus'):
+        print("⊘ Exodus support not available (OMEGA_H_USE_SEACASEXODUS not enabled)")
+        return
+    
+    # Build a simple 3D box mesh
+    print("Creating test mesh...")
+    mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        1.0, 1.0, 1.0,  # x, y, z dimensions
+        3, 3, 3,         # nx, ny, nz divisions
+        False            # symmetric
+    )
+    
+    # Get initial mesh properties
+    nverts_orig = mesh.nverts()
+    nelems_orig = mesh.nelems()
+    dim_orig = mesh.dim()
+    print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}")
+    
+    # Create temporary directory for test files
+    test_dir = "test_data"
+    os.makedirs(test_dir, exist_ok=True)
+    try:
+        # Test Exodus write
+        exodus_file = os.path.join(test_dir, "test_mesh.exo")
+        print(f"Writing to Exodus file: {exodus_file}")
+        py_pcms.write_mesh_exodus(exodus_file, mesh, verbose=False)
+        
+        # Test Exodus file handle API
+        if hasattr(py_pcms, 'exodus_open'):
+            print(f"Testing Exodus file handle API with: {exodus_file}")
+            
+            # Open the file
+            exo_handle = py_pcms.exodus_open(exodus_file, verbose=False)
+            print(f"Opened Exodus file with handle: {exo_handle}")
+            
+            # Get number of time steps
+            num_steps = py_pcms.exodus_get_num_time_steps(exo_handle)
+            print(f"Number of time steps: {num_steps}")
+            
+            # Read mesh using file handle
+            mesh_from_handle = py_pcms.OmegaHMesh(lib)
+            mesh_from_handle.set_comm(world)
+            py_pcms.read_mesh_exodus(exo_handle, mesh_from_handle, verbose=False)
+            
+            # Verify mesh properties
+            assert mesh_from_handle.dim() == dim_orig, f"Dimension mismatch: {mesh_from_handle.dim()} != {dim_orig}"
+            assert mesh_from_handle.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_from_handle.nverts()} != {nverts_orig}"
+            assert mesh_from_handle.nelems() == nelems_orig, f"Element count mismatch: {mesh_from_handle.nelems()} != {nelems_orig}"
+            
+            # Close the file
+            py_pcms.exodus_close(exo_handle)
+        
+        print("✓ Exodus I/O test passed")
+        
+    finally:
+        # Clean up temporary files
+        if os.path.exists(test_dir):
+            for item in os.listdir(test_dir):
+                item_path = os.path.join(test_dir, item)
+                if os.path.isfile(item_path):
+                    os.remove(item_path)
+                elif os.path.isdir(item_path):
+                    shutil.rmtree(item_path)
+            os.rmdir(test_dir)
+            print(f"Cleaned up test directory: {test_dir}")
+
+
+def test_adios2_io(lib, world):
+    """Test ADIOS2 file format I/O"""
+    print("\n=== Testing ADIOS2 Format I/O ===")
+    
+    # Check if ADIOS2 support is available
+    if not hasattr(py_pcms, 'write_mesh_adios2'):
+        print("⊘ ADIOS2 support not available (OMEGA_H_USE_ADIOS2 not enabled)")
+        return
+    
+    # Build a simple 3D box mesh
+    print("Creating test mesh...")
+    mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        1.0, 1.0, 1.0,  # x, y, z dimensions
+        3, 3, 3,         # nx, ny, nz divisions
+        False            # symmetric
+    )
+    
+    # Get initial mesh properties
+    nverts_orig = mesh.nverts()
+    nelems_orig = mesh.nelems()
+    dim_orig = mesh.dim()
+    print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}")
+    
+    # Create temporary directory for test files
+    test_dir = "test_data"
+    os.makedirs(test_dir, exist_ok=True)
+    try:
+        # Test ADIOS2 write/read
+        adios2_file = os.path.join(test_dir, "test_mesh.bp")
+        print(f"Writing to ADIOS2 file: {adios2_file}")
+        py_pcms.write_mesh_adios2(adios2_file, mesh, prefix="")
+        
+        print(f"Reading from ADIOS2 file: {adios2_file}")
+        mesh_read = py_pcms.read_mesh_adios2(adios2_file, lib, prefix="")
+        
+        # Verify mesh properties
+        assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}"
+        assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}"
+        assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}"
+        
+        print("✓ ADIOS2 I/O test passed")
+        
+    finally:
+        # Clean up temporary files
+        if os.path.exists(test_dir):
+            for item in os.listdir(test_dir):
+                item_path = os.path.join(test_dir, item)
+                if os.path.isfile(item_path):
+                    os.remove(item_path)
+                elif os.path.isdir(item_path):
+                    shutil.rmtree(item_path)
+            os.rmdir(test_dir)
+            print(f"Cleaned up test directory: {test_dir}")
+
+
+def test_read_mesh_file_auto_detect(lib, world):
+    """Test automatic format detection with read_mesh_file"""
+    print("\n=== Testing Automatic Format Detection ===")
+    
+    # Build a test mesh
+    print("Creating test mesh...")
+    mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        1.0, 1.0, 1.0,
+        3, 3, 3,
+        False
+    )
+    
+    nverts_orig = mesh.nverts()
+    nelems_orig = mesh.nelems()
+    
+    # Create temporary directory for test files
+    test_dir = "test_data"
+    os.makedirs(test_dir, exist_ok=True)
+    try:
+        # Write in different formats
+        binary_file = os.path.join(test_dir, "mesh.osh")
+        gmsh_file = os.path.join(test_dir, "mesh.msh")
+        
+        py_pcms.write_mesh_binary(binary_file, mesh)
+        py_pcms.write_mesh_gmsh(gmsh_file, mesh)
+        
+        # Test auto-detection for binary format
+        print(f"Auto-detecting binary file: {binary_file}")
+        mesh_binary = py_pcms.read_mesh_file(binary_file, world)
+        assert mesh_binary.nverts() == nverts_orig
+        assert mesh_binary.nelems() == nelems_orig
+        print("✓ Binary auto-detection passed")
+        
+        # Test auto-detection for Gmsh format
+        print(f"Auto-detecting Gmsh file: {gmsh_file}")
+        mesh_gmsh = py_pcms.read_mesh_file(gmsh_file, world)
+        assert mesh_gmsh.nverts() == nverts_orig
+        assert mesh_gmsh.nelems() == nelems_orig
+        print("✓ Gmsh auto-detection passed")
+        
+    finally:
+        # Clean up temporary files
+        if os.path.exists(test_dir):
+            for item in os.listdir(test_dir):
+                item_path = os.path.join(test_dir, item)
+                if os.path.isfile(item_path):
+                    os.remove(item_path)
+                elif os.path.isdir(item_path):
+                    shutil.rmtree(item_path)
+            os.rmdir(test_dir)
+            print(f"Cleaned up test directory: {test_dir}")
+
+if __name__ == "__main__":
+    print("=" * 60)
+    print("Omega_h File I/O Python Binding Tests")
+    print("=" * 60)
+    
+    # Create library and world communicator once for all tests
+    lib = py_pcms.OmegaHLibrary()
+    world = lib.world()
+    
+    try:
+        test_binary_io(lib, world)
+        test_gmsh_io(lib, world)
+        test_vtk_io(lib, world)
+        test_meshb_io(lib, world)
+        test_exodus_io(lib, world)
+        test_adios2_io(lib, world)
+        test_read_mesh_file_auto_detect(lib, world)
+        
+        print("\n" + "=" * 60)
+        print("✓ All tests passed!")
+        print("=" * 60)
+        
+    except Exception as e:
+        print("\n" + "=" * 60)
+        print(f"✗ Test failed with error: {e}")
+        print("=" * 60)
+        import traceback
+        traceback.print_exc()
+        exit(1)
+    finally:
+        # Explicitly delete objects to avoid an MPI finalizing issue
+        del world
+        del lib
diff --git a/src/pcms/pythonapi/test_omega_h_field.py b/src/pcms/pythonapi/test_omega_h_field.py
new file mode 100644
index 00000000..c0c2053a
--- /dev/null
+++ b/src/pcms/pythonapi/test_omega_h_field.py
@@ -0,0 +1,276 @@
+#!/usr/bin/env python3
+"""
+Test OmegaHFieldLayout Python bindings
+"""
+
+import py_pcms
+import numpy as np
+
+def test_layout_methods(world, dim, order, num_components):
+    """Test OmegaHFieldLayout methods"""
+    nx = 10
+    ny = 10 if dim > 1 else 0
+    nz = 10 if dim > 2 else 0
+
+    print(f"\nTesting layout: dim={dim}, order={order}, num_components={num_components}")
+
+    # Build mesh
+    mesh = py_pcms.build_box(
+        world, 
+        py_pcms.Family.SIMPLEX, 
+        1.0, 1.0, 1.0, 
+        nx, ny, nz, 
+        False
+    )
+
+    # Create layout
+    layout = py_pcms.create_lagrange_layout(
+        mesh, 
+        order, 
+        num_components,
+        py_pcms.CoordinateSystem.Cartesian
+    )
+
+    # Test layout methods
+    num_comp = layout.get_num_components()
+    assert num_comp == num_components, f"Expected {num_components} components, got {num_comp}"
+
+    num_owned = layout.get_num_owned_dof_holder()
+    print(f"  Number of owned DOF holders: {num_owned}")
+    assert num_owned > 0, "Expected at least one owned DOF holder"
+
+    num_global = layout.get_num_global_dof_holder()
+    print(f"  Number of global DOF holders: {num_global}")
+    assert num_global >= num_owned, "Global DOF holders should be >= owned"
+
+    # Get ownership mask
+    owned = layout.get_owned()
+    print(f"  Ownership mask length: {len(owned)}")
+    assert len(owned) > 0, "Ownership mask should not be empty"
+
+    # Get global IDs
+    gids = layout.get_gids()
+    print(f"  Global IDs length: {len(gids)}")
+    assert len(gids) > 0, "Global IDs should not be empty"
+
+    # Get coordinates of DOF holders (2D parametric coordinates)
+    coords = layout.get_dof_holder_coordinates()
+    print(f"  Coordinates shape: {coords.shape}")
+    assert coords.shape[0] > 0, "Should have coordinates"
+    assert coords.shape[1] == 2, "Coordinates should be 2D (parametric)"
+
+    # Check if distributed
+    is_distributed = layout.is_distributed()
+    print(f"  Is distributed: {is_distributed}")
+
+    # Get entity offsets
+    ent_offsets = layout.get_ent_offsets()
+    print(f"  Entity offsets: {ent_offsets}")
+    assert len(ent_offsets) == 5, f"Expected 5 entity offsets"
+
+    # Get nodes per dimension
+    nodes = layout.get_nodes_per_dim()
+    print(f"  Nodes per dim: {nodes}")
+    assert len(nodes) == 4, "Should have 4 entries (verts, edges, faces, regions)"
+
+    # Get number of entities
+    num_ents = layout.get_num_ents()
+    print(f"  Number of entities: {num_ents}")
+    assert num_ents > 0, "Should have at least one entity"
+
+    # Get sizes
+    owned_size = layout.owned_size()
+    global_size = layout.global_size()
+    print(f"  Owned size: {owned_size}")
+    print(f"  Global size: {global_size}")
+    assert owned_size == num_owned * num_components
+    assert global_size == num_global * num_components
+
+    # Create a field from the layout
+    field = layout.create_field()
+    print(f"  Created field: {type(field)}")
+
+    # Test setting and getting data
+    ndata = num_owned * num_components
+    test_data = np.arange(ndata, dtype=np.float64)
+    field.set_dof_holder_data(test_data)
+
+    retrieved_data = field.get_dof_holder_data()
+    assert len(retrieved_data) == ndata, f"Expected {ndata} elements, got {len(retrieved_data)}"
+
+    # Verify data matches
+    # debug print
+    for i in range(min(10, ndata)):
+        print(f"  Data[{i}]: set={test_data[i]}, got={retrieved_data[i]}")
+    differences = np.abs(test_data - retrieved_data)
+    assert np.all(differences < 1e-12), "Data mismatch after set/get"
+
+    print(f"✓ Test passed: dim={dim}, order={order}, num_components={num_components}")
+
+def test_field_evaluation(world, dim, order, num_components):
+    """Test OmegaHField evaluation at points"""
+    nx = 10
+    ny = 10 if dim > 1 else 0
+    nz = 10 if dim > 2 else 0
+
+    print(f"\nTesting field evaluation: dim={dim}, order={order}, num_components={num_components}")
+
+    # Build mesh
+    mesh = py_pcms.build_box(
+        world, 
+        py_pcms.Family.SIMPLEX, 
+        1.0, 1.0, 1.0, 
+        nx, ny, nz, 
+        False
+    )
+
+    # Create layout
+    layout = py_pcms.create_lagrange_layout(
+        mesh, 
+        order, 
+        num_components,
+        py_pcms.CoordinateSystem.Cartesian
+    )
+
+    # Create field
+    field = layout.create_field()
+
+    # Set up test data - use a simple function: f(x,y,z) = sin(x*y) for component 0, etc.
+    print(f"  Setting up field data...")
+    # Get mesh vertex coordinates to set field values
+    mesh_coords = mesh.coords()
+    num_verts = mesh.nverts()
+    print(f"  Mesh has {num_verts} vertices")
+
+    # For order 1, we only need vertex values
+    # For order 2, we also need edge midpoint values
+    num_owned = layout.get_num_owned_dof_holder()
+    ndata = num_owned * num_components
+    print(f"  Setting up field data for {ndata} DOFs")
+
+    # Define test function
+    def test_func(x, y, z, component):
+        return np.sin(5.0 * x * y) + float(component)
+
+    # Create test data array
+    test_data = np.zeros(ndata, dtype=np.float64)
+
+    # Set vertex values
+    for i in range(min(num_verts, num_owned)):
+        x = mesh_coords[i * dim + 0] if dim >= 1 else 0.0
+        y = mesh_coords[i * dim + 1] if dim >= 2 else 0.0
+        z = mesh_coords[i * dim + 2] if dim >= 3 else 0.0
+        for c in range(num_components):
+            idx = i * num_components + c
+            test_data[idx] = test_func(x, y, z, c)
+
+    print(f"  Set vertex values for field data")
+
+    # For order 2, set edge midpoint values (simplified - would need proper edge coordinates)
+    if order == 2 and num_owned > num_verts:
+        for i in range(num_verts, num_owned):
+            for c in range(num_components):
+                idx = i * num_components + c
+                # Use simplified values for edges
+                test_data[idx] = float(c) + 0.5
+
+    field.set_dof_holder_data(test_data)
+    print(f"  Set field data with test function")
+
+    # Define evaluation points (physical coordinates)
+    if dim == 2:
+        eval_coords = np.array([
+            [0.5, 0.5],
+            [0.25, 0.25],
+            [0.75, 0.75],
+            [0.1, 0.9],
+            [0.9, 0.1]
+        ], dtype=np.float64)
+    else:  # dim == 3
+        eval_coords = np.array([
+            [0.5, 0.5, 0.5],
+            [0.25, 0.25, 0.25],
+            [0.75, 0.75, 0.75],
+            [0.1, 0.9, 0.1],
+            [0.9, 0.1, 0.9]
+        ], dtype=np.float64)
+
+    num_eval_points = eval_coords.shape[0]
+    print(f"  Evaluating at {num_eval_points} points")
+
+    # Get localization hint for the coordinates
+    coord_system = py_pcms.CoordinateSystem.Cartesian
+    location_hint = field.get_localization_hint(eval_coords, coord_system)
+    print(f"  Got localization hint")
+
+    # Create output buffer for evaluation results
+    eval_values = np.zeros(num_eval_points * num_components, dtype=np.float64)
+
+    # Evaluate the field
+    field.evaluate(location_hint, eval_values, coord_system)
+    print(f"  Field evaluated successfully")
+
+    # Print and verify results
+    for i in range(num_eval_points):
+        coords_str = ", ".join([f"{eval_coords[i, j]:.3f}" for j in range(dim)])
+        values_str = ", ".join([f"{eval_values[i * num_components + c]:.4f}" 
+                               for c in range(num_components)])
+        print(f"  Point {i} ({coords_str}): values = [{values_str}]")
+
+        # Verify we got valid values (not NaN or infinity)
+        for c in range(num_components):
+            val = eval_values[i * num_components + c]
+            assert not np.isnan(val), f"Got NaN value at point {i}, component {c}"
+            assert not np.isinf(val), f"Got inf value at point {i}, component {c}"
+
+    # Test gradient evaluation if available
+    if field.can_evaluate_gradient():
+        print(f"  Testing gradient evaluation...")
+        gradient_values = np.zeros(num_eval_points * num_components * dim, dtype=np.float64)
+
+        try:
+            field.evaluate_gradient(gradient_values, coord_system)
+            print(f"  Gradient evaluated successfully")
+
+            # Print gradient results
+            for i in range(min(3, num_eval_points)):
+                for c in range(num_components):
+                    grad_components = []
+                    for d in range(dim):
+                        idx = i * num_components * dim + c * dim + d
+                        grad_components.append(f"{gradient_values[idx]:.4f}")
+                    grad_str = ", ".join(grad_components)
+                    print(f"    Point {i}, component {c}: gradient = [{grad_str}]")
+        except Exception as e:
+            print(f"  Gradient evaluation failed: {e}")
+    else:
+        print(f"  Gradient evaluation not supported")
+
+    print(f"✓ Field evaluation test passed: dim={dim}, order={order}, num_components={num_components}")
+
+def main():
+    """Run all test cases"""
+    print("Testing OmegaHFieldLayout Python bindings...")
+
+    # Initialize Omega_h library
+    lib = py_pcms.OmegaHLibrary()
+    world = lib.world()
+    print("Initialized Omega_h library and world")
+
+    # Test different configurations
+    test_layout_methods(world, 2, 1, 1)
+    test_layout_methods(world, 2, 2, 1)
+    # test_layout_methods(world, 2, 1, 3)
+    # test_layout_methods(world, 3, 1, 1)
+
+    # Test field evaluation
+    print("\n" + "="*60)
+    print("Testing field evaluation...")
+    print("="*60)
+    test_field_evaluation(world, 2, 1, 1)
+    test_field_evaluation(world, 2, 2, 1)
+
+    print("\n✓ All tests passed!")
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file
diff --git a/src/pcms/pythonapi/test_uniform_grid_field.py b/src/pcms/pythonapi/test_uniform_grid_field.py
new file mode 100644
index 00000000..06109c78
--- /dev/null
+++ b/src/pcms/pythonapi/test_uniform_grid_field.py
@@ -0,0 +1,476 @@
+"""
+Example test file for uniform grid field Python bindings.
+
+This demonstrates how to use the UniformGridField and UniformGridFieldLayout
+classes in Python, similar to how OmegaHField is used.
+"""
+import numpy as np
+import py_pcms
+
+
+def test_uniform_grid_field_creation():
+    """Test creating a 2D uniform grid field."""
+    # Create a 2D uniform grid
+    grid = py_pcms.UniformGrid2D()
+    grid.bot_left = [0.0, 0.0]
+    grid.edge_length = [10.0, 10.0]
+    grid.divisions = [4, 4]  # 4x4 cells
+    
+    print(f"Grid cells: {grid.get_num_cells()}")  # Should be 16
+    
+    # Create field layout with 1 component (scalar field)
+    layout = py_pcms.UniformGridFieldLayout2D(
+        grid, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    
+    print(f"Num components: {layout.get_num_components()}")
+    print(f"Num owned DOF holders: {layout.get_num_owned_dof_holder()}")
+    print(f"Num vertices: {layout.get_num_vertices()}")  # 5x5 = 25 vertices
+    
+    # Create field
+    field = layout.create_field()
+    print(f"Field created: {field is not None}")
+    
+    return grid, layout, field
+
+
+def test_uniform_grid_field_data_operations():
+    """Test setting and getting field data."""
+    # Create a simple grid
+    grid = py_pcms.UniformGrid2D()
+    grid.bot_left = [0.0, 0.0]
+    grid.edge_length = [10.0, 10.0]
+    grid.divisions = [2, 2]  # 2x2 cells, 3x3 vertices
+    
+    layout = py_pcms.UniformGridFieldLayout2D(
+        grid, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    ug_field = layout.create_field()
+    
+    # Get the number of DOF holders (vertices)
+    num_vertices = layout.get_num_vertices()
+    print(f"Number of vertices: {num_vertices}")  # Should be 9
+    
+    # Create data for vertices (initialize with values)
+    data = np.arange(num_vertices, dtype=np.float64)
+    print(f"Setting data: {data}")
+    
+    # Set field data
+    ug_field.set_dof_holder_data(data)
+    
+    # Get field data back
+    retrieved_data = ug_field.get_dof_holder_data()
+    print(f"Retrieved data: {retrieved_data}")
+    
+    # Verify they match
+    assert np.allclose(data, retrieved_data)
+    print("Data successfully set and retrieved!")
+
+
+def test_uniform_grid_field_mdspan_2d():
+    """Test 2D mdspan output as numpy array."""
+    grid = py_pcms.UniformGrid2D()
+    grid.bot_left = [0.0, 0.0]
+    grid.edge_length = [10.0, 10.0]
+    grid.divisions = [2, 3]  # 3x4 vertices
+
+    layout = py_pcms.UniformGridFieldLayout2D(
+        grid, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    ug_field = layout.create_field()
+
+    num_vertices = layout.get_num_vertices()
+    data = np.arange(num_vertices, dtype=np.float64)
+    ug_field.set_dof_holder_data(data)
+
+    mdspan = ug_field.to_mdspan()
+    expected = data.reshape((grid.divisions[0] + 1, grid.divisions[1] + 1))
+    assert type(mdspan) == np.ndarray
+
+    assert mdspan.shape == expected.shape
+    assert np.allclose(mdspan, expected)
+    print("2D mdspan output verified!")
+
+
+def test_uniform_grid_field_evaluation():
+    """Test evaluating field at specific coordinates."""
+    # Create grid
+    grid = py_pcms.UniformGrid2D()
+    grid.bot_left = [0.0, 0.0]
+    grid.edge_length = [10.0, 10.0]
+    grid.divisions = [2, 2]
+    
+    layout = py_pcms.UniformGridFieldLayout2D(
+        grid, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    ug_field = layout.create_field()
+    
+    # Set simple linear field data: f(x,y) = x + y
+    num_vertices = layout.get_num_vertices()
+    coords = layout.get_dof_holder_coordinates()
+    data = np.zeros(num_vertices)
+    for i in range(num_vertices):
+        data[i] = coords[i, 0] + coords[i, 1]
+    
+    ug_field.set_dof_holder_data(data)
+    
+    # Evaluate at some points
+    eval_coords = np.array([[2.5, 2.5], [7.5, 7.5]])
+    
+    # Get localization hint
+    hint = ug_field.get_localization_hint(
+        eval_coords, py_pcms.CoordinateSystem.Cartesian
+    )
+    
+    # Evaluate field
+    results = np.zeros(len(eval_coords))
+    ug_field.evaluate(
+        hint, results, py_pcms.CoordinateSystem.Cartesian
+    )
+    
+    print(f"Evaluation results: {results}")
+    print(f"Expected (approximately): [5.0, 15.0]")
+
+
+def test_uniform_grid_closest_cell():
+    """Test finding closest cell to a point."""
+    grid = py_pcms.UniformGrid2D()
+    grid.bot_left = [0.0, 0.0]
+    grid.edge_length = [10.0, 10.0]
+    grid.divisions = [4, 4]
+    
+    # Test point in middle of grid
+    point = np.array([5.0, 5.0])
+    cell_id = grid.closest_cell_id(point)
+    print(f"Point {point} is in cell {cell_id}")
+    
+    # Test point outside grid (should clamp to closest cell)
+    point_outside = np.array([-1.0, -1.0])
+    cell_id_outside = grid.closest_cell_id(point_outside)
+    print(f"Point {point_outside} (outside) maps to cell {cell_id_outside}")
+
+
+def test_3d_uniform_grid():
+    """Test 3D uniform grid field."""
+    grid = py_pcms.UniformGrid3D()
+    grid.bot_left = [0.0, 0.0, 0.0]
+    grid.edge_length = [10.0, 10.0, 10.0]
+    grid.divisions = [2, 2, 2]  # 2x2x2 cells
+    
+    print(f"3D Grid cells: {grid.get_num_cells()}")  # Should be 8
+    
+    layout = py_pcms.UniformGridFieldLayout3D(
+        grid, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    
+    print(f"3D Grid vertices: {layout.get_num_vertices()}")  # 3x3x3 = 27
+    
+    ug_field = layout.create_field()
+    
+    # Set and get data
+    num_vertices = layout.get_num_vertices()
+    data = np.ones(num_vertices) * 42.0
+    ug_field.set_dof_holder_data(data)
+    
+    retrieved = ug_field.get_dof_holder_data()
+    assert np.allclose(data, retrieved)
+    print("3D field data successfully set and retrieved!")
+
+
+def test_uniform_grid_field_mdspan_3d():
+    """Test 3D mdspan output as numpy array."""
+    grid = py_pcms.UniformGrid3D()
+    grid.bot_left = [0.0, 0.0, 0.0]
+    grid.edge_length = [10.0, 10.0, 10.0]
+    grid.divisions = [2, 1, 3]  # 3x2x4 vertices
+
+    layout = py_pcms.UniformGridFieldLayout3D(
+        grid, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    ug_field = layout.create_field()
+
+    num_vertices = layout.get_num_vertices()
+    data = np.arange(num_vertices, dtype=np.float64)
+    ug_field.set_dof_holder_data(data)
+
+    mdspan = ug_field.to_mdspan()
+    expected = data.reshape(
+        (grid.divisions[0] + 1, grid.divisions[1] + 1, grid.divisions[2] + 1)
+    )
+
+    assert mdspan.shape == expected.shape
+    assert np.allclose(mdspan, expected)
+    print("3D mdspan output verified!")
+
+
+def test_uniform_grid_workflow(world):
+    """
+    Test complete UniformGrid workflow with Omega_h mesh and field interpolation.
+    
+    This test:
+    1. Creates an Omega_h mesh
+    2. Creates a uniform grid from the mesh
+    3. Creates a binary mask field
+    4. Creates and initializes an Omega_h field with f(x,y) = x + 2*y
+    5. Transfers the field to uniform grid via interpolation
+    6. Verifies the transferred values
+    """
+    # Create a simple 2D box mesh: 1.0 x 1.0 domain with 4x4 elements
+    mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        1.0, 1.0, 0.0,
+        4, 4, 0,
+        False
+    )
+    
+    # Create uniform grid from mesh with 4x4 divisions
+    grid = py_pcms.create_uniform_grid_from_mesh(mesh, [4, 4])
+    print(f"Created uniform grid with {grid.get_num_cells()} cells")
+    
+    # Create binary mask field (returns tuple of (layout, field))
+    mask_layout, mask_field = py_pcms.create_uniform_grid_binary_field(mesh, [4, 4])
+    mask_data = mask_field.get_dof_holder_data()
+    print(f"Created mask field with {len(mask_data)} vertices")
+    
+    # Create Omega_h field layout with linear elements
+    omega_h_layout = py_pcms.create_lagrange_layout(
+        mesh, 1, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    omega_h_field = omega_h_layout.create_field()
+    
+    # Initialize omega_h field with f(x,y) = x + 2*y
+    coords = omega_h_layout.get_dof_holder_coordinates()
+    num_nodes = omega_h_layout.get_num_owned_dof_holder()
+    omega_h_data = np.zeros(num_nodes)
+    
+    for i in range(num_nodes):
+        x = coords[i, 0]
+        y = coords[i, 1]
+        omega_h_data[i] = x + 2.0 * y
+    
+    omega_h_field.set_dof_holder_data(omega_h_data)
+    print(f"Initialized Omega_h field with {num_nodes} nodes")
+    
+    # Create uniform grid field layout
+    ug_layout = py_pcms.UniformGridFieldLayout2D(
+        grid, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    ug_field = ug_layout.create_field()
+    
+    # Transfer from omega_h field to uniform grid field using interpolation
+    py_pcms.interpolate_field(omega_h_field, ug_field)
+    print("Field interpolation completed")
+    
+    # Get uniform grid field data and coordinates
+    ug_field_data = ug_field.get_dof_holder_data()
+    ug_coords = ug_layout.get_dof_holder_coordinates()
+    
+    # Verify ug_field values
+    print("\nVerifying uniform grid field values:")
+    num_errors = 0
+    for j in range(grid.divisions[1] + 1):
+        for i in range(grid.divisions[0] + 1):
+            vertex_id = j * (grid.divisions[0] + 1) + i
+            x = ug_coords[vertex_id, 0]
+            y = ug_coords[vertex_id, 1]
+            expected = x + 2.0 * y
+            actual = ug_field_data[vertex_id]
+            
+            error = abs(expected - actual)
+            if error > 1e-10:
+                num_errors += 1
+                print(f"  Vertex ({i}, {j}) at ({x:.2f}, {y:.2f}): "
+                      f"expected {expected:.4f}, got {actual:.4f}, error {error:.2e}")
+            
+    if num_errors == 0:
+        print("All uniform grid field values verified successfully!")
+    else:
+        print(f"Found {num_errors} verification errors")
+        
+    # Verify mask field values (all should be 1 for vertices inside the mesh)
+    print("\nVerifying mask field values:")
+    mask_errors = 0
+    nx, ny = grid.divisions[0] + 1, grid.divisions[1] + 1
+    print(f"\nBinary mask field ({nx}x{ny} vertices) visualization:")
+    print("(Each position shows the mask value at that vertex)")
+    for j in range(ny - 1, -1, -1):  # Print from top to bottom
+        row = []
+        for i in range(nx):
+            vertex_id = j * nx + i
+            mask_value = mask_data[vertex_id]
+            row.append(str(int(mask_value)))
+            if mask_value != 1.0:
+                mask_errors += 1
+        print(f"  Row {j}: [{' '.join(row)}]")
+    
+    if mask_errors == 0:
+        print("All mask field values verified successfully!")
+    else:
+        print(f"Found {mask_errors} mask errors")
+    
+    # Serialize the uniform grid field
+    print("\nSerializing uniform grid field:")
+    num_vertices = ug_layout.get_num_vertices()
+    buffer = np.zeros(num_vertices)
+    permutation = np.arange(num_vertices, dtype=np.int32)  # Identity permutation
+    
+    bytes_written = ug_field.serialize(buffer, permutation)
+    print(f"Serialized {bytes_written} values to buffer")
+    
+    print(f"\nSerialized uniform grid field data ({grid.divisions[0]+1}x{grid.divisions[1]+1} vertices) visualization:")
+    print("(Each vertex shows its field value)")
+    for j in range(grid.divisions[1], -1, -1):  # Print from top to bottom
+        row_values = []
+        for i in range(grid.divisions[0] + 1):
+            vertex_id = j * (grid.divisions[0] + 1) + i
+            value = buffer[vertex_id]
+            row_values.append(f"{value:6.3f}")
+        print(f"  Row {j}: [{' '.join(row_values)}]")
+    
+    # Also print with coordinates for reference
+    print(f"\nDetailed vertex information:")
+    for j in range(grid.divisions[1] + 1):
+        for i in range(grid.divisions[0] + 1):
+            vertex_id = j * (grid.divisions[0] + 1) + i
+            x = ug_coords[vertex_id, 0]
+            y = ug_coords[vertex_id, 1]
+            value = buffer[vertex_id]
+            print(f"  V[{i},{j}] (id={vertex_id:2d}) at ({x:.3f}, {y:.3f}): value = {value:.6f}")
+    
+    # Test deserialization
+    print("\nTesting deserialization:")
+    ug_field_copy = ug_layout.create_field()
+    ug_field_copy.deserialize(buffer, permutation)
+    
+    ug_field_copy_data = ug_field_copy.get_dof_holder_data()
+    if np.allclose(ug_field_data, ug_field_copy_data):
+        print("✓ Deserialization successful - data matches!")
+    else:
+        print("✗ Deserialization failed - data mismatch!")
+        print(f"  Max difference: {np.max(np.abs(ug_field_data - ug_field_copy_data))}")
+        
+    assert num_errors == 0, f"Field verification failed with {num_errors} errors"
+    assert mask_errors == 0, f"Mask verification failed with {mask_errors} errors"
+
+
+def test_omega_h_to_omega_h_transfer_workflow(world):
+    """
+    Test field transfer from one Omega_h mesh to another Omega_h mesh.
+
+    This test:
+    1. Creates a source Omega_h mesh
+    2. Creates a target Omega_h mesh
+    3. Creates Lagrange layouts and fields on both
+    4. Initializes the source with f(x,y) = x + 2*y
+    5. Transfers the field to the target mesh
+    6. Verifies the transferred values at target DOF holders
+    """
+    # Source mesh (coarser)
+    src_mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        1.0, 1.0, 0.0,
+        4, 4, 0,
+        False
+    )
+
+    # Target mesh (finer)
+    tgt_mesh = py_pcms.build_box(
+        world,
+        py_pcms.Family.SIMPLEX,
+        1.0, 1.0, 0.0,
+        8, 8, 0,
+        False
+    )
+
+    # Create Lagrange layouts and fields
+    src_layout = py_pcms.create_lagrange_layout(
+        src_mesh, 1, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    tgt_layout = py_pcms.create_lagrange_layout(
+        tgt_mesh, 1, 1, py_pcms.CoordinateSystem.Cartesian
+    )
+    src_field = src_layout.create_field()
+    tgt_field = tgt_layout.create_field()
+
+    # Initialize source field with f(x,y) = x + 2*y
+    src_coords = src_layout.get_dof_holder_coordinates()
+    src_num_nodes = src_layout.get_num_owned_dof_holder()
+    src_data = np.zeros(src_num_nodes)
+    for i in range(src_num_nodes):
+        x = src_coords[i, 0]
+        y = src_coords[i, 1]
+        src_data[i] = x + 2.0 * y
+    src_field.set_dof_holder_data(src_data)
+
+    # Transfer field to target mesh
+    py_pcms.interpolate_field(src_field, tgt_field)
+
+    # Verify target field values at target DOF holders
+    tgt_coords = tgt_layout.get_dof_holder_coordinates()
+    tgt_data = tgt_field.get_dof_holder_data()
+    errors = 0
+    for i in range(tgt_layout.get_num_owned_dof_holder()):
+        x = tgt_coords[i, 0]
+        y = tgt_coords[i, 1]
+        expected = x + 2.0 * y
+        actual = tgt_data[i]
+        if abs(expected - actual) > 1e-10:
+            errors += 1
+    assert errors == 0, f"Omega_h transfer verification failed with {errors} errors"
+
+
+if __name__ == "__main__":
+    lib = py_pcms.OmegaHLibrary()
+    world = lib.world()
+    print("=" * 60)
+    print("Testing UniformGrid Field Creation")
+    print("=" * 60)
+    test_uniform_grid_field_creation()
+    
+    print("\n" + "=" * 60)
+    print("Testing UniformGrid Field Data Operations")
+    print("=" * 60)
+    test_uniform_grid_field_data_operations()
+
+    print("\n" + "=" * 60)
+    print("Testing UniformGrid Field mdspan (2D)")
+    print("=" * 60)
+    test_uniform_grid_field_mdspan_2d()
+    
+    print("\n" + "=" * 60)
+    print("Testing UniformGrid Field Evaluation")
+    print("=" * 60)
+    test_uniform_grid_field_evaluation()
+    
+    print("\n" + "=" * 60)
+    print("Testing Closest Cell ID")
+    print("=" * 60)
+    test_uniform_grid_closest_cell()
+    
+    print("\n" + "=" * 60)
+    print("Testing 3D UniformGrid")
+    print("=" * 60)
+    test_3d_uniform_grid()
+
+    print("\n" + "=" * 60)
+    print("Testing UniformGrid Field mdspan (3D)")
+    print("=" * 60)
+    test_uniform_grid_field_mdspan_3d()
+    
+    print("\n" + "=" * 60)
+    print("Testing UniformGrid Workflow")
+    print("=" * 60)
+    test_uniform_grid_workflow(world)
+
+    print("\n" + "=" * 60)
+    print("Testing Omega_h to Omega_h Transfer Workflow")
+    print("=" * 60)
+    test_omega_h_to_omega_h_transfer_workflow(world)
+    
+    print("\n" + "=" * 60)
+    print("All tests passed!")
+    print("=" * 60)
+
+    del world
diff --git a/src/pcms/pythonapi/uniform_grid_workflow_guide.md b/src/pcms/pythonapi/uniform_grid_workflow_guide.md
new file mode 100644
index 00000000..8614e4da
--- /dev/null
+++ b/src/pcms/pythonapi/uniform_grid_workflow_guide.md
@@ -0,0 +1,264 @@
+# Unstructured mesh field transfer to structured mesh field workflow guide
+
+This guide demonstrates how to transfer unstructured mesh fields to uniform grid fields in PCMS Python API, including mesh creation, field interpolation, and data serialization.
+
+## Overview
+
+The workflow allows you to:
+1. Create a structured uniform grid as the bounding box as an unstructured Omega_h mesh
+2. Create binary mask fields indicating which cells are inside the mesh
+3. Interpolate field data from Omega_h mesh fields to uniform grid fields
+4. Serialize and deserialize field data for I/O operations
+
+## Complete Workflow Example
+
+### Step 1: Initialize Omega_h Library
+
+Initialize the Omega_h library which manages parallel communication.
+
+```python
+import py_pcms
+import numpy as np
+
+# Create library and world communicator
+lib = py_pcms.OmegaHLibrary()
+world = lib.world()
+```
+
+
+---
+
+### Step 2: Create an Omega_h Mesh
+
+You can either create a mesh or read an existing mesh from a file.
+
+#### Option A: Create Mesh Programmatically
+
+```python
+# Create a 2D box mesh: 1.0 x 1.0 domain with 4x4 elements
+mesh = py_pcms.build_box(
+    world,                          # Communicator
+    py_pcms.Family.SIMPLEX,        # Element type
+    1.0, 1.0, 0.0,                 # Domain size: x, y, z
+    4, 4, 0,                       # Number of elements: nx, ny, nz
+    False                          # Symmetric flag
+)
+```
+
+#### Option B: Read Mesh from File
+
+```python
+# Auto-detect file format and read
+mesh = py_pcms.read_mesh_file("my_mesh.osh", world)
+
+# Or use format-specific readers for explicit control:
+
+# Binary format (.osh)
+mesh = py_pcms.read_mesh_binary("mesh.osh", lib)
+
+# Gmsh format (.msh)
+mesh = py_pcms.read_mesh_gmsh("mesh.msh", world)
+
+# VTK format (.pvtu for parallel files)
+mesh = py_pcms.read_mesh_parallel_vtk("mesh.pvtu", world)
+
+# Exodus format (.exo) - if SEACAS support is enabled
+# Using file handle API
+exo_handle = py_pcms.exodus_open("mesh.exo", verbose=False)
+num_steps = py_pcms.exodus_get_num_time_steps(exo_handle)
+mesh = py_pcms.OmegaHMesh(lib)
+mesh.set_comm(world)
+py_pcms.read_mesh_exodus(exo_handle, mesh, verbose=False)
+py_pcms.exodus_close(exo_handle)
+
+# ADIOS2 format (.bp) - if ADIOS2 support is enabled
+mesh = py_pcms.read_mesh_adios2("mesh.bp", lib, prefix="")
+
+# MESHB format (.mesh) - if LIBMESHB support is enabled
+mesh = py_pcms.OmegaHMesh(lib)
+mesh.set_comm(world)
+py_pcms.read_mesh_meshb(mesh, "mesh.mesh")
+# Read solution/resolution data if available
+py_pcms.read_meshb_sol(mesh, "mesh.sol", sol_name="resolution")
+```
+
+**PS**: File formats such as Exodus or libMeshB require that PCMS be built with the respective libraries enabled.
+
+Checking Mesh Properties:
+
+```python
+print(f"Mesh dimension: {mesh.dim()}")
+print(f"Number of vertices: {mesh.nverts()}")
+print(f"Number of elements: {mesh.nelems()}")
+print(f"Mesh family: {mesh.family()}")
+```
+
+---
+
+### Step 3: Create Uniform Grid from Mesh
+
+Generate a structured uniform grid that covers the bounding box of the mesh.
+
+```python
+# Create uniform grid with 4x4 cell divisions
+grid = py_pcms.create_uniform_grid_from_mesh(mesh, [4, 4])
+print(f"Created uniform grid with {grid.get_num_cells()} cells")
+```
+
+---
+
+### Step 4: Create Binary Mask Field
+
+Create a mask where each uniform grid vertex is marked as 1 (inside mesh) or 0 (outside mesh).
+
+```python
+# Create binary mask indicating which vertices are inside the mesh
+# Returns tuple of (layout, field) - layout must be kept alive while using field
+mask_layout, mask_field = py_pcms.create_uniform_grid_binary_field(mesh, [4, 4])
+print(f"Created mask field with {mask_layout.get_num_vertices()} vertices")
+```
+
+**Accessing Values**:
+```python
+# Access mask value for vertex (i, j)
+vertex_id = j * (grid.divisions[0] + 1) + i
+mask_data = mask_field.get_dof_holder_data()
+mask_value = mask_data[vertex_id]  # Returns 0.0 or 1.0
+```
+
+---
+
+### Step 5: Create and Initialize Omega_h Field
+
+Create a field on the unstructured mesh and initialize it with data.
+
+```python
+# Create field layout with linear (order=1) elements and 1 component (scalar)
+omega_h_layout = py_pcms.create_lagrange_layout(
+    mesh, 
+    1,
+    1,
+    py_pcms.CoordinateSystem.Cartesian
+)
+omega_h_field = omega_h_layout.create_field()
+
+# Initialize field with f(x,y) = x + 2*y
+coords = omega_h_layout.get_dof_holder_coordinates()
+num_nodes = omega_h_layout.get_num_owned_dof_holder()
+omega_h_data = np.zeros(num_nodes)
+
+for i in range(num_nodes):
+    x = coords[i, 0]
+    y = coords[i, 1]
+    omega_h_data[i] = x + 2.0 * y
+
+omega_h_field.set_dof_holder_data(omega_h_data)
+```
+
+---
+
+### Step 6: Create Uniform Grid Field Layout
+
+Set up the data structure for storing field values on the uniform grid vertices.
+
+```python
+# Create uniform grid field layout
+ug_layout = py_pcms.UniformGridFieldLayout2D(
+    grid, 
+    1,                                      # Number of components
+    py_pcms.CoordinateSystem.Cartesian
+)
+ug_field = ug_layout.create_field()
+```
+
+---
+
+### Step 7: Interpolate Field Data
+
+Interpolate field values from the unstructured mesh to the structured uniform grid vertices.
+
+```python
+# Transfer field from Omega_h mesh to uniform grid
+py_pcms.interpolate_field(omega_h_field, ug_field)
+```
+
+---
+
+### Step 8: Access Field Data
+
+Retrieve interpolated values and verify correctness.
+
+```python
+# Get field data and coordinates
+ug_field_data = ug_field.get_dof_holder_data()
+ug_coords = ug_layout.get_dof_holder_coordinates()
+
+# Access data at vertex (i, j)
+vertex_id = j * (grid.divisions[0] + 1) + i
+value = ug_field_data[vertex_id]
+x, y = ug_coords[vertex_id, 0], ug_coords[vertex_id, 1]
+```
+
+---
+
+### Step 9: Export Field Data with `to_mdspan`
+
+Convert field data to a structured $x \times y$ (or $x \times y \times z$) array
+for downstream analyses.
+
+```python
+# Convert field data to a structured array
+grid_values = ug_field.to_mdspan()  # 2D: (nx+1, ny+1), 3D: (nx+1, ny+1, nz+1)
+mask_values = mask_field.to_mdspan()  # 2D: (nx+1, ny+1), 3D: (nx+1, ny+1, nz+1)
+
+# Save to file (example)
+np.save('field_data.npy', grid_values)
+```
+
+---
+
+## API Reference Summary
+
+### Grid Creation
+- `create_uniform_grid_from_mesh(mesh, divisions)` - Create grid from mesh
+- `create_uniform_grid_binary_field(mesh, divisions)` - Create inside/outside mask
+
+### Field Layout
+- `create_lagrange_layout(mesh, order, num_components, coord_system)` - Omega_h field layout
+- `UniformGridFieldLayout2D(grid, num_components, coord_system)` - 2D grid layout
+- `UniformGridFieldLayout3D(grid, num_components, coord_system)` - 3D grid layout
+
+### Field Operations
+- `layout.create_field()` - Create field from layout
+- `field.set_dof_holder_data(data)` - Set field values
+- `field.get_dof_holder_data()` - Get field values
+- `field.to_mdspan()` - Get field values as a structured array
+- `interpolate_field(source_field, target_field)` - Interpolate between fields
+
+### Exodus I/O (if OMEGA_H_USE_SEACASEXODUS enabled)
+- `exodus_open(filepath, verbose=False)` - Open Exodus file, returns file handle
+- `exodus_close(exodus_file)` - Close Exodus file handle
+- `exodus_get_num_time_steps(exodus_file)` - Get number of time steps in file
+- `read_mesh_exodus(exodus_file, mesh, verbose=False, classify_with=...)` - Read mesh from file handle
+- `read_exodus_nodal_fields(exodus_file, mesh, time_step, prefix="", postfix="", verbose=False)` - Read nodal fields
+- `read_exodus_element_fields(exodus_file, mesh, time_step, prefix="", postfix="", verbose=False)` - Read element fields
+- `read_mesh_exodus_sliced(filepath, comm, verbose=False, classify_with=..., time_step=-1)` - Read mesh in parallel
+- `write_mesh_exodus(filepath, mesh, verbose=False, classify_with=...)` - Write mesh to Exodus file
+- `ExodusClassifyWith.NODE_SETS` - Classify using node sets
+- `ExodusClassifyWith.SIDE_SETS` - Classify using side sets
+
+### MESHB I/O (if OMEGA_H_USE_LIBMESHB enabled)
+- `read_mesh_meshb(mesh, filepath)` - Read mesh from MESHB file
+- `write_mesh_meshb(mesh, filepath, version)` - Write mesh to MESHB file
+- `read_meshb_sol(mesh, filepath, sol_name)` - Read solution/resolution from .sol file
+- `write_meshb_sol(mesh, filepath, sol_name, version)` - Write solution/resolution to .sol file
+
+### Grid Properties
+- `grid.get_num_cells()` - Total number of cells
+- `grid.divisions` - Cell divisions in each dimension
+- `layout.get_num_vertices()` - Total number of vertices
+- `layout.get_dof_holder_coordinates()` - Vertex coordinates
+
+---
+
+For more examples and advanced usage, see the test files in the `pcms/pythonapi/` directory.
diff --git a/src/pcms/transfer_field2.h b/src/pcms/transfer_field2.h
index 60a0266f..df5712ef 100644
--- a/src/pcms/transfer_field2.h
+++ b/src/pcms/transfer_field2.h
@@ -34,8 +34,8 @@ void interpolate_field2(const FieldT<T>& source, FieldT<T>& target)
   }
 
   auto coords = target.GetLayout().GetDOFHolderCoordinates();
-  std::vector<Real> evaluation(coords.GetCoordinates().size() / 2);
-  FieldDataView<Real, HostMemorySpace> data_view{make_array_view(evaluation),
+  std::vector<T> evaluation(coords.GetCoordinates().size() / 2);
+  FieldDataView<T, HostMemorySpace> data_view{make_array_view(evaluation),
                                                  source.GetCoordinateSystem()};
   auto locale = source.GetLocalizationHint(coords);
   source.Evaluate(locale, data_view);
diff --git a/src/pcms/uniform_grid.h b/src/pcms/uniform_grid.h
index 5ec13d19..f1320f1d 100644
--- a/src/pcms/uniform_grid.h
+++ b/src/pcms/uniform_grid.h
@@ -2,6 +2,8 @@
 #define PCMS_COUPLING_UNIFORM_GRID_H
 #include "pcms/bounding_box.h"
 #include "Omega_h_vector.hpp"
+#include "Omega_h_bbox.hpp"
+#include "Omega_h_mesh.hpp"
 #include <numeric>
 namespace pcms
 {
@@ -119,6 +121,60 @@ struct UniformGrid
 
 using Uniform2DGrid = UniformGrid<>;
 
+/**
+ * \brief Create a uniform grid layout from an Omega_h mesh as a bounding box
+ *
+ * This function computes the bounding box of the given mesh and creates a
+ * uniform Cartesian grid that covers the entire mesh domain. The grid cells
+ * are uniformly distributed across each dimension.
+ *
+ * \tparam dim Spatial dimension of the mesh (2 or 3)
+ * \param mesh The Omega_h mesh to create the grid from
+ * \param divisions Array specifying the number of cells in each dimension
+ * \return UniformGrid<dim> A uniform grid covering the mesh bounding box
+ *
+ */
+template <unsigned dim = 2>
+UniformGrid<dim> CreateUniformGridFromMesh(Omega_h::Mesh& mesh,
+                                           const std::array<LO, dim>& divisions)
+{
+  // Get the bounding box of the mesh
+  auto bbox = Omega_h::get_bounding_box<dim>(&mesh);
+
+  // Calculate edge lengths and bottom-left corner
+  std::array<Real, dim> edge_length;
+  std::array<Real, dim> bot_left;
+
+  for (unsigned i = 0; i < dim; ++i) {
+    bot_left[i] = bbox.min[i];
+    edge_length[i] = bbox.max[i] - bbox.min[i];
+  }
+
+  return UniformGrid<dim>{
+    .edge_length = edge_length, .bot_left = bot_left, .divisions = divisions};
+}
+
+/**
+ * \brief Create a uniform grid with equal divisions in all dimensions
+ *
+ * This is a convenience function that creates a uniform grid with the same
+ * number of cells in each dimension.
+ *
+ * \tparam dim Spatial dimension of the mesh (2 or 3)
+ * \param mesh The Omega_h mesh to create the grid from
+ * \param cells_per_dim Number of cells per dimension (same for all dimensions)
+ * \return UniformGrid<dim> A uniform grid covering the mesh bounding box
+ *
+ */
+template <unsigned dim = 2>
+UniformGrid<dim> CreateUniformGridFromMesh(Omega_h::Mesh& mesh,
+                                           LO cells_per_dim)
+{
+  std::array<LO, dim> divisions;
+  divisions.fill(cells_per_dim);
+  return CreateUniformGridFromMesh<dim>(mesh, divisions);
+}
+
 } // namespace pcms
 
 #endif // PCMS_COUPLING_UNIFORM_GRID_H
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 4325dab2..bfe767af 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -394,6 +394,7 @@ if(Catch2_FOUND)
       test_normalisation.cpp
       test_spline_interpolator.cpp
       test_svd_serial.cpp
+      test_uniform_grid_field.cpp
       test_interpolation_class.cpp)
   endif()
   add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES})
diff --git a/test/test_uniform_grid.cpp b/test/test_uniform_grid.cpp
index 9cdc3f39..df145357 100644
--- a/test/test_uniform_grid.cpp
+++ b/test/test_uniform_grid.cpp
@@ -1,8 +1,12 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/catch_approx.hpp>
 #include <pcms/uniform_grid.h>
+#include <Omega_h_build.hpp>
+#include <Omega_h_library.hpp>
 
+using pcms::CreateUniformGridFromMesh;
 using pcms::Uniform2DGrid;
+using pcms::UniformGrid;
 
 TEST_CASE("uniform grid")
 {
@@ -50,3 +54,119 @@ TEST_CASE("uniform grid")
     REQUIRE(119 == uniform_grid.GetCellIndex({11, 9}));
   }
 }
+
+TEST_CASE("Create uniform grid from omega_h mesh")
+{
+  // Initialize Omega_h library
+  auto lib = Omega_h::Library{};
+  auto world = lib.world();
+
+  // Create a simple 2D box mesh: 1.0 x 1.0 domain with 4x4 elements
+  auto mesh =
+    Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 4, 4, 0, false);
+
+  SECTION("Create grid with specified divisions per dimension")
+  {
+    // Create a 10x12 uniform grid over the mesh
+    auto grid = CreateUniformGridFromMesh<2>(mesh, {10, 12});
+
+    // Verify the grid properties
+    REQUIRE(grid.GetNumCells() == 120);
+    REQUIRE(grid.edge_length[0] == Catch::Approx(1.0));
+    REQUIRE(grid.edge_length[1] == Catch::Approx(1.0));
+    REQUIRE(grid.bot_left[0] == Catch::Approx(0.0));
+    REQUIRE(grid.bot_left[1] == Catch::Approx(0.0));
+    REQUIRE(grid.divisions[0] == 10);
+    REQUIRE(grid.divisions[1] == 12);
+
+    // Test cell bounding boxes
+    auto bbox_first = grid.GetCellBBOX(0);
+    REQUIRE(bbox_first.center[0] == Catch::Approx(0.05)); // 1.0/10/2
+    REQUIRE(bbox_first.center[1] == Catch::Approx(0.5 / 12.0));
+
+    auto bbox_last = grid.GetCellBBOX(119);
+    REQUIRE(bbox_last.center[0] == Catch::Approx(0.95));
+    REQUIRE(bbox_last.center[1] == Catch::Approx(11.5 / 12.0));
+  }
+
+  SECTION("Create grid with equal divisions")
+  {
+    // Create a 15x15 uniform grid
+    auto grid = CreateUniformGridFromMesh<2>(mesh, 15);
+
+    REQUIRE(grid.GetNumCells() == 225);
+    REQUIRE(grid.divisions[0] == 15);
+    REQUIRE(grid.divisions[1] == 15);
+
+    // Verify all cells have equal size
+    auto cell_width = grid.edge_length[0] / grid.divisions[0];
+    auto cell_height = grid.edge_length[1] / grid.divisions[1];
+    REQUIRE(cell_width == Catch::Approx(1.0 / 15.0));
+    REQUIRE(cell_height == Catch::Approx(1.0 / 15.0));
+  }
+
+  SECTION("Test with non-unit domain")
+  {
+    // Create a larger domain: 5.0 x 3.0
+    auto large_mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 5.0, 3.0, 0.0,
+                                         10, 6, 0, false);
+
+    auto grid = CreateUniformGridFromMesh<2>(large_mesh, {20, 30});
+
+    REQUIRE(grid.GetNumCells() == 600);
+    REQUIRE(grid.edge_length[0] == Catch::Approx(5.0));
+    REQUIRE(grid.edge_length[1] == Catch::Approx(3.0));
+    REQUIRE(grid.bot_left[0] == Catch::Approx(0.0));
+    REQUIRE(grid.bot_left[1] == Catch::Approx(0.0));
+
+    // Test cell sizes
+    auto bbox = grid.GetCellBBOX(0);
+    REQUIRE(bbox.half_width[0] == Catch::Approx(5.0 / 20.0 / 2.0));
+    REQUIRE(bbox.half_width[1] == Catch::Approx(3.0 / 30.0 / 2.0));
+  }
+
+  SECTION("Test point location in grid cells")
+  {
+    auto grid = CreateUniformGridFromMesh<2>(mesh, {10, 10});
+
+    // Test corner points
+    REQUIRE(grid.ClosestCellID(Omega_h::Vector<2>{0.0, 0.0}) == 0);
+    REQUIRE(grid.ClosestCellID(Omega_h::Vector<2>{1.0, 1.0}) == 99);
+
+    // Test middle point
+    REQUIRE(grid.ClosestCellID(Omega_h::Vector<2>{0.5, 0.5}) == 55);
+  }
+}
+
+TEST_CASE("Create uniform grid from 3D mesh")
+{
+  auto lib = Omega_h::Library{};
+  auto world = lib.world();
+
+  // Create a 3D box mesh: 2.0 x 3.0 x 1.5 domain
+  auto mesh =
+    Omega_h::build_box(world, OMEGA_H_SIMPLEX, 2.0, 3.0, 1.5, 4, 6, 3, false);
+
+  SECTION("Create 3D grid with specified divisions")
+  {
+    auto grid = CreateUniformGridFromMesh<3>(mesh, {10, 15, 8});
+
+    REQUIRE(grid.GetNumCells() == 1200); // 10 * 15 * 8
+    REQUIRE(grid.edge_length[0] == Catch::Approx(2.0));
+    REQUIRE(grid.edge_length[1] == Catch::Approx(3.0));
+    REQUIRE(grid.edge_length[2] == Catch::Approx(1.5));
+    REQUIRE(grid.divisions[0] == 10);
+    REQUIRE(grid.divisions[1] == 15);
+    REQUIRE(grid.divisions[2] == 8);
+  }
+
+  SECTION("Create 3D grid with equal divisions")
+  {
+    auto grid = CreateUniformGridFromMesh<3>(mesh, 12);
+
+    REQUIRE(grid.GetNumCells() == 1728); // 12^3
+    REQUIRE(grid.divisions[0] == 12);
+    REQUIRE(grid.divisions[1] == 12);
+    REQUIRE(grid.divisions[2] == 12);
+  }
+}
diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp
new file mode 100644
index 00000000..f04ec029
--- /dev/null
+++ b/test/test_uniform_grid_field.cpp
@@ -0,0 +1,740 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/catch_approx.hpp>
+#include <Kokkos_Core.hpp>
+#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h"
+#include "pcms/adapter/uniform_grid/uniform_grid_field.h"
+#include "pcms/uniform_grid.h"
+#include "Omega_h_library.hpp"
+#include "Omega_h_build.hpp"
+#include "pcms/transfer_field2.h"
+#include "pcms/adapter/omega_h/omega_h_field_layout.h"
+#include "pcms/adapter/omega_h/omega_h_field2.h"
+#include "pcms/create_field.h"
+#include "pcms/arrays.h"
+#include <cmath>
+
+using pcms::CreateUniformGridBinaryField;
+using pcms::CreateUniformGridFromMesh;
+
+// Helper function to initialize omega_h field data with f(x,y) = x + 2*y
+std::vector<pcms::Real> CreateOmegaHFieldData(
+  const pcms::CoordinateView<pcms::HostMemorySpace>& coords, int num_nodes) {
+  auto coords_data = coords.GetCoordinates();
+  std::vector<pcms::Real> omega_h_data(num_nodes);
+  for (int i = 0; i < num_nodes; ++i) {
+    pcms::Real x = coords_data(i, 0);
+    pcms::Real y = coords_data(i, 1);
+    omega_h_data[i] = x + 2.0 * y; // f(x,y) = x + 2y
+  }
+  return omega_h_data;
+}
+
+// Helper function to verify ug_field values
+void VerifyUniformGridFieldValues(
+  const pcms::UniformGrid<2>& grid,
+  const pcms::CoordinateView<pcms::HostMemorySpace>& ug_coords,
+  const pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace>& ug_field_data) {
+  fprintf(stderr, "\nVerifying ug_field values:\n");
+  for (int j = 0; j <= grid.divisions[1]; ++j) {
+    for (int i = 0; i <= grid.divisions[0]; ++i) {
+      int vertex_id = j * (grid.divisions[0] + 1) + i;
+      pcms::Real x = ug_coords.GetCoordinates()(vertex_id, 0);
+      pcms::Real y = ug_coords.GetCoordinates()(vertex_id, 1);
+      pcms::Real expected = x + 2.0 * y;
+      pcms::Real actual = ug_field_data(vertex_id);
+      fprintf(
+        stderr,
+        "ug_field Vertex (%d, %d) at (%.2f, %.2f): expected %.4f, got %.4f\n",
+        i, j, x, y, expected, actual);
+      REQUIRE(std::abs(expected - actual) <= 1e-10);
+    }
+  }
+}
+
+// Helper function to verify mask field values
+void VerifyMaskFieldValues(
+  const pcms::UniformGrid<2>& grid,
+  const pcms::UniformGridField<2>& mask_field) {
+  fprintf(stderr, "\nVerifying mask_field values:\n");
+  auto mask_data = mask_field.GetDOFHolderData();
+  for (int j = 0; j <= grid.divisions[1]; ++j) {
+    for (int i = 0; i <= grid.divisions[0]; ++i) {
+      int vertex_id = j * (grid.divisions[0] + 1) + i;
+      pcms::Real mask_value = mask_data(vertex_id);
+      fprintf(stderr, "Vertex (%d, %d): mask = %.0f\n", i, j, mask_value);
+      REQUIRE(mask_value == 1.0);
+    }
+  }
+}
+
+TEST_CASE("UniformGrid field creation")
+{
+  // Create a simple 2D uniform grid
+  pcms::UniformGrid<2> grid;
+  grid.bot_left = {0.0, 0.0};
+  grid.edge_length = {10.0, 10.0};
+  grid.divisions = {5, 5};
+
+  // Create field layout with 1 component (scalar field)
+  pcms::UniformGridFieldLayout<2> layout(grid, 1,
+                                         pcms::CoordinateSystem::Cartesian);
+
+  REQUIRE(layout.GetNumComponents() == 1);
+  REQUIRE(layout.GetNumOwnedDofHolder() == 36); // (5+1)x(5+1) = 36 vertices
+  REQUIRE(layout.GetNumGlobalDofHolder() == 36);
+  REQUIRE_FALSE(layout.IsDistributed());
+
+  // Create field
+  auto field = layout.CreateField();
+  REQUIRE(field != nullptr);
+}
+
+TEST_CASE("UniformGrid field data operations", "[uniform_grid_field]")
+{
+  pcms::UniformGrid<2> grid;
+  grid.bot_left = {0.0, 0.0};
+  grid.edge_length = {10.0, 10.0};
+  grid.divisions = {4, 4}; // 4x4 = 16 cells
+
+  pcms::UniformGridFieldLayout<2> layout(grid, 1,
+                                         pcms::CoordinateSystem::Cartesian);
+  auto field = layout.CreateField();
+
+  // Initialize field data with vertex indices (5x5 = 25 vertices for 4x4 cells)
+  std::vector<pcms::Real> data(25);
+  for (size_t i = 0; i < 25; ++i) {
+    data[i] = static_cast<pcms::Real>(i);
+  }
+
+  auto data_view = pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace>(
+    data.data(), data.size());
+  field->SetDOFHolderData(data_view);
+
+  // Retrieve data
+  auto retrieved = field->GetDOFHolderData();
+  REQUIRE(retrieved.size() == 25);
+
+  for (size_t i = 0; i < 25; ++i) {
+    REQUIRE(retrieved[i] == static_cast<pcms::Real>(i));
+  }
+}
+
+TEST_CASE("UniformGrid field evaluation - piecewise constant")
+{
+  pcms::UniformGrid<2> grid;
+  grid.bot_left = {0.0, 0.0};
+  grid.edge_length = {10.0, 10.0};
+  grid.divisions = {2, 2}; // 2x2 = 4 cells
+
+  pcms::UniformGridFieldLayout<2> layout(grid, 1,
+                                         pcms::CoordinateSystem::Cartesian);
+  auto field = layout.CreateField();
+
+  // Set vertex values for a 2x2 cell grid (3x3 = 9 vertices)
+  // Vertex layout:
+  //   v6---v7---v8
+  //   |  2 |  3 |
+  //   v3---v4---v5
+  //   |  0 |  1 |
+  //   v0---v1---v2
+  std::vector<pcms::Real> data = {
+    1.0, 1.5, 2.0, // v0, v1, v2 (bottom row, y=0)
+    2.0, 2.5, 3.0, // v3, v4, v5 (middle row, y=5)
+    3.0, 3.5, 4.0  // v6, v7, v8 (top row, y=10)
+  };
+  auto data_view = pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace>(
+    data.data(), data.size());
+  field->SetDOFHolderData(data_view);
+
+  // Evaluate at cell centers
+  std::vector<pcms::Real> eval_coords = {
+    2.5, 2.5, // Cell 0 center
+    7.5, 2.5, // Cell 1 center
+    2.5, 7.5, // Cell 2 center
+    7.5, 7.5  // Cell 3 center
+  };
+
+  auto coords_view = pcms::Rank2View<const pcms::Real, pcms::HostMemorySpace>(
+    eval_coords.data(), 4, 2);
+  auto coord_view = pcms::CoordinateView<pcms::HostMemorySpace>(
+    pcms::CoordinateSystem::Cartesian, coords_view);
+
+  auto hint = field->GetLocalizationHint(coord_view);
+
+  std::vector<pcms::Real> results(4);
+  auto results_view =
+    pcms::Rank1View<pcms::Real, pcms::HostMemorySpace>(results.data(), 4);
+  auto results_field_view =
+    pcms::FieldDataView<pcms::Real, pcms::HostMemorySpace>(
+      results_view, pcms::CoordinateSystem::Cartesian);
+
+  field->Evaluate(hint, results_field_view);
+
+  // 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[0] - 1.75) < 1e-10);
+  REQUIRE(std::abs(results[1] - 2.25) < 1e-10);
+  REQUIRE(std::abs(results[2] - 2.75) < 1e-10);
+  REQUIRE(std::abs(results[3] - 3.25) < 1e-10);
+}
+
+TEST_CASE("UniformGrid field serialization")
+{
+  pcms::UniformGrid<2> grid;
+  grid.bot_left = {0.0, 0.0};
+  grid.edge_length = {10.0, 10.0};
+  grid.divisions = {3, 3}; // 9 cells
+
+  pcms::UniformGridFieldLayout<2> layout(grid, 1,
+                                         pcms::CoordinateSystem::Cartesian);
+  auto field = layout.CreateField();
+
+  // Set field data (4x4 = 16 vertices for 3x3 cells)
+  std::vector<pcms::Real> data(16);
+  for (size_t i = 0; i < 16; ++i) {
+    data[i] = static_cast<pcms::Real>(i * 10);
+  }
+
+  auto data_view = pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace>(
+    data.data(), data.size());
+  field->SetDOFHolderData(data_view);
+
+  // Create identity permutation
+  std::vector<pcms::LO> permutation(16);
+  for (size_t i = 0; i < 16; ++i) {
+    permutation[i] = i;
+  }
+
+  auto perm_view = pcms::Rank1View<const pcms::LO, pcms::HostMemorySpace>(
+    permutation.data(), 16);
+
+  // Serialize
+  std::vector<pcms::Real> buffer(16);
+  auto buffer_view =
+    pcms::Rank1View<pcms::Real, pcms::HostMemorySpace>(buffer.data(), 16);
+  int size = field->Serialize(buffer_view, perm_view);
+
+  REQUIRE(size == 16);
+
+  // Verify serialized data
+  for (size_t i = 0; i < 16; ++i) {
+    REQUIRE(buffer[i] == data[i]);
+  }
+
+  // Create new field and deserialize
+  auto field2 = layout.CreateField();
+  auto buffer_const_view =
+    pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace>(buffer.data(), 16);
+  field2->Deserialize(buffer_const_view, perm_view);
+
+  // Verify deserialized data
+  auto retrieved = field2->GetDOFHolderData();
+  for (size_t i = 0; i < 16; ++i) {
+    REQUIRE(retrieved[i] == data[i]);
+  }
+}
+
+TEST_CASE("UniformGrid field copy")
+{
+  pcms::UniformGrid<2> grid;
+  grid.bot_left = {0.0, 0.0};
+  grid.edge_length = {10.0, 10.0};
+  grid.divisions = {2, 2}; // 2x2 grid
+
+  pcms::UniformGridFieldLayout<2> layout(grid, 1,
+                                         pcms::CoordinateSystem::Cartesian);
+  auto field = layout.CreateField();
+
+  // Set vertex values with f(x,y) = x + y at 3x3 vertex positions
+  // Vertices at: (0,0), (5,0), (10,0), (0,5), (5,5), (10,5), (0,10), (5,10),
+  // (10,10)
+  std::vector<pcms::Real> data = {
+    0.0,  5.0,  10.0, // y=0:  v0(0,0)=0,   v1(5,0)=5,   v2(10,0)=10
+    5.0,  10.0, 15.0, // y=5:  v3(0,5)=5,   v4(5,5)=10,  v5(10,5)=15
+    10.0, 15.0, 20.0  // y=10: v6(0,10)=10, v7(5,10)=15, v8(10,10)=20
+  };
+  auto data_view = pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace>(
+    data.data(), data.size());
+  field->SetDOFHolderData(data_view);
+
+  // Test 1: Evaluate at cell centers
+  std::vector<pcms::Real> eval_coords = {
+    2.5, 2.5, // Cell 0 center
+    7.5, 2.5, // Cell 1 center
+    2.5, 7.5, // Cell 2 center
+    7.5, 7.5  // Cell 3 center
+  };
+
+  auto coords_view = pcms::Rank2View<const pcms::Real, pcms::HostMemorySpace>(
+    eval_coords.data(), 4, 2);
+  auto coord_view = pcms::CoordinateView<pcms::HostMemorySpace>(
+    pcms::CoordinateSystem::Cartesian, coords_view);
+
+  auto hint = field->GetLocalizationHint(coord_view);
+
+  std::vector<pcms::Real> results(coord_view.GetCoordinates().size() / 2);
+  auto results_view =
+    pcms::Rank1View<pcms::Real, pcms::HostMemorySpace>(results.data(), 4);
+  auto results_field_view =
+    pcms::FieldDataView<pcms::Real, pcms::HostMemorySpace>(
+      results_view, pcms::CoordinateSystem::Cartesian);
+
+  field->Evaluate(hint, results_field_view);
+
+  REQUIRE(std::abs(results[0] - 5.0) < 1e-10);
+  REQUIRE(std::abs(results[1] - 10.0) < 1e-10);
+  REQUIRE(std::abs(results[2] - 10.0) < 1e-10);
+  REQUIRE(std::abs(results[3] - 15.0) < 1e-10);
+
+  // Test 2: test copy_field2
+  auto field2 = layout.CreateField();
+  pcms::copy_field2(*field, *field2);
+
+  auto copied_data = field2->GetDOFHolderData();
+  REQUIRE(copied_data.size() == data.size());
+  for (size_t i = 0; i < data.size(); ++i) {
+    REQUIRE(copied_data[i] == data[i]);
+  }
+}
+
+TEST_CASE("Transfer from OmegaH field to UniformGrid field")
+{
+  Omega_h::Library lib;
+
+  // Create a simple omega_h mesh (2x2 box)
+  auto mesh = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 2,
+                                 2, 0, false);
+
+  // Create OmegaH field layout with linear elements
+  auto omega_h_layout =
+    pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
+  auto omega_h_field = omega_h_layout->CreateField();
+
+  // Initialize omega_h field with a simple function f(x,y) = x + 2*y
+  auto coords = omega_h_layout->GetDOFHolderCoordinates();
+  int num_nodes = omega_h_layout->GetNumOwnedDofHolder();
+  std::vector<pcms::Real> omega_h_data = CreateOmegaHFieldData(coords, num_nodes);
+
+  auto omega_h_data_view =
+    pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace>(
+      omega_h_data.data(), omega_h_data.size());
+  omega_h_field->SetDOFHolderData(omega_h_data_view);
+
+  // Create a uniform grid field covering the same domain [0,1] x [0,1]
+  pcms::UniformGrid<2> grid;
+  grid.bot_left = {0.0, 0.0};
+  grid.edge_length = {1.0, 1.0};
+  grid.divisions = {2, 2}; // 4x4 grid for finer resolution
+
+  pcms::UniformGridFieldLayout<2> ug_layout(grid, 1,
+                                            pcms::CoordinateSystem::Cartesian);
+  auto ug_field = ug_layout.CreateField();
+
+  // Transfer from omega_h field to uniform grid field using interpolation
+  auto coords_interpolation = ug_layout.GetDOFHolderCoordinates();
+  std::vector<pcms::Real> evaluation(
+    coords_interpolation.GetCoordinates().size() / 2);
+  auto evaluation_view = pcms::Rank1View<pcms::Real, pcms::HostMemorySpace>(
+    evaluation.data(), evaluation.size());
+  pcms::FieldDataView<pcms::Real, pcms::HostMemorySpace> data_view{
+    evaluation_view, omega_h_field->GetCoordinateSystem()};
+  auto locale = omega_h_field->GetLocalizationHint(coords_interpolation);
+  omega_h_field->Evaluate(locale, data_view);
+  auto evaluation_view_const =
+    pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace>(evaluation.data(),
+                                                             evaluation.size());
+  ug_field->SetDOFHolderData(evaluation_view_const);
+
+  // Verify the transferred data at uniform grid vertices
+  auto transferred_data = ug_field->GetDOFHolderData();
+  auto ug_coords = ug_layout.GetDOFHolderCoordinates();
+  auto ug_coords_data = ug_coords.GetCoordinates();
+  int num_ug_nodes = ug_layout.GetNumOwnedDofHolder(); // 5x5 = 25 vertices
+
+  // Check a few sample points
+  for (int i = 0; i < num_ug_nodes; ++i) {
+    pcms::Real x = ug_coords_data(i, 0);
+    pcms::Real y = ug_coords_data(i, 1);
+    pcms::Real expected = x + 2.0 * y;
+    pcms::Real actual = transferred_data[i];
+
+    // Allow some tolerance for interpolation
+    REQUIRE(std::abs(actual - expected) < 1e-6);
+  }
+}
+
+TEST_CASE("Create binary field from uniform grid")
+{
+  auto lib = Omega_h::Library{};
+  auto world = lib.world();
+
+  SECTION("Simple 2D box mesh - all vertices inside")
+  {
+    // Create a mesh that fills the domain [0,1] x [0,1]
+    auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 10,
+                                   10, 0, false);
+
+    // Create a 5x5 grid (coarser than mesh)
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {5, 5});
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 36);  // (5+1) * (5+1) = 36 vertices
+
+    // All grid vertices should be inside the mesh
+    pcms::Real sum = 0.0;
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      sum += field_data(i);
+    }
+    REQUIRE(sum == 36.0);
+  }
+
+  SECTION("Binary field with custom divisions")
+  {
+    auto mesh =
+      Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 8, 8, 0, false);
+
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {10, 8});
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 99);  // (10+1) * (8+1) = 99 vertices
+
+    // Most vertices should be inside
+    pcms::Real inside_count = 0.0;
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      inside_count += field_data(i);
+    }
+    REQUIRE(inside_count > 0.0);
+    REQUIRE(inside_count <= 99.0);
+  }
+
+  SECTION("Binary field with equal divisions convenience function")
+  {
+    auto mesh =
+      Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 5, 5, 0, false);
+
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 8);
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 81);  // (8+1) * (8+1) = 81 vertices
+
+    pcms::Real inside_count = 0.0;
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      inside_count += field_data(i);
+    }
+    REQUIRE(inside_count > 0.0);
+  }
+
+  SECTION("Fine grid over coarse mesh")
+  {
+    // Create a simple mesh
+    auto mesh =
+      Omega_h::build_box(world, OMEGA_H_SIMPLEX, 2.0, 2.0, 0.0, 4, 4, 0, false);
+
+    // Create a fine grid (20x20)
+    auto grid = CreateUniformGridFromMesh<2>(mesh, {20, 20});
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {20, 20});
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 441);  // (20+1) * (20+1) = 441 vertices
+
+    // Verify consistency: check some specific vertices
+    // Center vertex should be inside (vertex at i=10, j=10)
+    int center_idx = 10 * 21 + 10;  // 21 vertices per row
+    REQUIRE(field_data(center_idx) == 1.0);
+
+    // Corner vertex should be inside
+    int corner_idx = 0;  // vertex (0, 0)
+    REQUIRE(field_data(corner_idx) == 1.0);
+  }
+
+  SECTION("Verify field values are binary")
+  {
+    auto mesh =
+      Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 5, 5, 0, false);
+
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 10);
+    auto field_data = field->GetDOFHolderData();
+
+    // All values should be 0 or 1
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      pcms::Real val = field_data(i);
+      REQUIRE((val == 0.0 || val == 1.0));
+    }
+  }
+
+  SECTION("Grid extends beyond mesh - vertices outside should be marked 0")
+  {
+    // Create a small mesh in the center of a domain
+    auto mesh =
+      Omega_h::build_box(world, OMEGA_H_SIMPLEX, 0.5, 0.5, 0.0, 5, 5, 0, false);
+
+    // The mesh occupies [0, 0.5] x [0, 0.5]
+    // Create a grid that would cover this
+    auto grid = CreateUniformGridFromMesh<2>(mesh, {10, 10});
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {10, 10});
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 121);  // (10+1) * (10+1) = 121 vertices
+
+    // All vertices should be inside since grid is exactly on mesh bbox
+    pcms::Real inside_count = 0.0;
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      inside_count += field_data(i);
+    }
+    REQUIRE(inside_count > 0.0);
+  }
+
+  SECTION("Test with different aspect ratio")
+  {
+    // Create a rectangular mesh
+    auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 3.0, 1.0, 0.0, 12, 4,
+                                   0, false);
+
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {30, 10});
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 341);  // (30+1) * (10+1) = 341 vertices
+
+    pcms::Real inside_count = 0.0;
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      inside_count += field_data(i);
+    }
+    REQUIRE(inside_count > 0.0);
+
+    // Calculate percentage inside
+    double inside_percent = 100.0 * inside_count / field_data.extent(0);
+    // Most vertices should be inside
+    REQUIRE(inside_percent > 50.0);
+  }
+}
+
+TEST_CASE("Binary field integration with grid methods")
+{
+  auto lib = Omega_h::Library{};
+  auto world = lib.world();
+
+  auto mesh =
+    Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 10, 10, 0, false);
+
+  SECTION("Query field value at specific grid vertex")
+  {
+    auto grid = CreateUniformGridFromMesh<2>(mesh, {8, 8});
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {8, 8});
+    auto field_data = field->GetDOFHolderData();
+
+    // Get field value for a specific vertex (middle vertex at i=4, j=4)
+    pcms::LO vertex_id = 4 * 9 + 4;  // 9 = (8+1) vertices per row
+    REQUIRE(field_data(vertex_id) == 1.0);  // Should be inside
+
+    // Compute vertex position
+    pcms::Real dx = grid.edge_length[0] / grid.divisions[0];
+    pcms::Real dy = grid.edge_length[1] / grid.divisions[1];
+    pcms::Real x = grid.bot_left[0] + 4 * dx;
+    pcms::Real y = grid.bot_left[1] + 4 * dy;
+    
+    // Verify it's roughly in the middle
+    REQUIRE(x == Catch::Approx(0.5).margin(0.1));
+    REQUIRE(y == Catch::Approx(0.5).margin(0.1));
+  }
+
+  SECTION("Count vertices by region")
+  {
+    auto grid = CreateUniformGridFromMesh<2>(mesh, 10);
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 10);
+    auto field_data = field->GetDOFHolderData();
+
+    // Count vertices in different quadrants
+    int q1 = 0, q2 = 0, q3 = 0, q4 = 0; // quadrants
+
+    pcms::Real dx = grid.edge_length[0] / grid.divisions[0];
+    pcms::Real dy = grid.edge_length[1] / grid.divisions[1];
+    
+    for (int j = 0; j <= grid.divisions[1]; ++j) {
+      for (int i = 0; i <= grid.divisions[0]; ++i) {
+        pcms::LO vertex_id = j * (grid.divisions[0] + 1) + i;
+        if (field_data(vertex_id) == 1.0) {
+          pcms::Real x = grid.bot_left[0] + i * dx;
+          pcms::Real y = grid.bot_left[1] + j * dy;
+          
+          if (x < 0.5 && y < 0.5)
+            q1++;
+          else if (x >= 0.5 && y < 0.5)
+            q2++;
+          else if (x < 0.5 && y >= 0.5)
+            q3++;
+          else
+            q4++;
+        }
+      }
+    }
+
+    // All quadrants should have some inside vertices
+    REQUIRE(q1 > 0);
+    REQUIRE(q2 > 0);
+    REQUIRE(q3 > 0);
+    REQUIRE(q4 > 0);
+
+    // For 11x11 vertices, boundary at x=0.5, y=0.5 splits asymmetrically:
+    // q1 (x<0.5, y<0.5): 5x5=25, q2 (x>=0.5, y<0.5): 6x5=30
+    // q3 (x<0.5, y>=0.5): 5x6=30, q4 (x>=0.5, y>=0.5): 6x6=36
+    int total = q1 + q2 + q3 + q4;
+    REQUIRE(total == 121); // All vertices inside
+    REQUIRE(q1 == 25);
+    REQUIRE(q2 == 30);
+    REQUIRE(q3 == 30);
+    REQUIRE(q4 == 36);
+  }
+}
+
+TEST_CASE("Performance and edge cases")
+{
+  auto lib = Omega_h::Library{};
+  auto world = lib.world();
+
+  SECTION("Very fine grid")
+  {
+    auto mesh =
+      Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 5, 5, 0, false);
+
+    // Create a very fine grid
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 50);
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 2601);  // (50+1) * (50+1) = 2601 vertices
+
+    pcms::Real inside_count = 0.0;
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      inside_count += field_data(i);
+    }
+    REQUIRE(inside_count > 0.0);
+  }
+
+  SECTION("Coarse grid")
+  {
+    auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 10,
+                                   10, 0, false);
+
+    // Very coarse grid
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {2, 2});
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 9);  // (2+1) * (2+1) = 9 vertices
+
+    // All vertices should be inside for this configuration
+    pcms::Real inside_count = 0.0;
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      inside_count += field_data(i);
+    }
+    REQUIRE(inside_count > 0.0);
+  }
+
+  SECTION("Non-square domain")
+  {
+    auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 5.0, 2.0, 0.0, 20, 8,
+                                   0, false);
+
+    auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {25, 10});
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 286);  // (25+1) * (10+1) = 286 vertices
+
+    pcms::Real inside_count = 0.0;
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      inside_count += field_data(i);
+    }
+    REQUIRE(inside_count > 0.0);
+  }
+
+  SECTION("Grid larger than mesh - vertices outside marked as 0")
+  {
+    // Create a mesh covering [0, 0.5] x [0, 0.5]
+    auto mesh =
+      Omega_h::build_box(world, OMEGA_H_SIMPLEX, 0.5, 0.5, 0.0, 5, 5, 0, false);
+
+    // Manually create a larger grid covering [0, 1] x [0, 1]
+    pcms::UniformGrid<2> grid;
+    grid.edge_length = {1.0, 1.0};
+    grid.bot_left = {0.0, 0.0};
+    grid.divisions = {10, 10};
+
+    // Create binary field on the larger grid
+    auto [layout, field] = pcms::CreateUniformGridBinaryFieldFromGrid<2>(mesh, grid);
+    auto field_data = field->GetDOFHolderData();
+
+    REQUIRE(field_data.extent(0) == 121);  // (10+1) * (10+1) = 121 vertices
+
+    // Count vertices inside and outside
+    pcms::Real inside_count = 0.0;
+    for (size_t i = 0; i < field_data.extent(0); ++i) {
+      inside_count += field_data(i);
+    }
+    pcms::Real outside_count = field_data.extent(0) - inside_count;
+
+    // Should have both inside (1) and outside (0) vertices
+    REQUIRE(inside_count > 0.0);
+    REQUIRE(outside_count > 0.0);
+
+    // Check specific vertices
+    // Bottom-left corner should be inside (mesh covers [0, 0.5])
+    int corner_id = 0 * 11 + 0;  // vertex (0, 0)
+    REQUIRE(field_data(corner_id) == 1.0);
+
+    // Top-right corner should be outside (mesh ends at 0.5)
+    corner_id = 10 * 11 + 10;  // vertex (10, 10)
+    REQUIRE(field_data(corner_id) == 0.0);
+
+    // Vertex at i=6, j=6 (coords ~0.6, ~0.6) should be outside
+    corner_id = 6 * 11 + 6;
+    REQUIRE(field_data(corner_id) == 0.0);
+
+    // Vertex at i=2, j=2 (coords ~0.2, ~0.2) should be inside
+    auto center_id = 2 * 11 + 2;
+    REQUIRE(field_data(center_id) == 1.0);
+  }
+}
+
+TEST_CASE("UniformGrid workflow")
+{
+  auto lib = Omega_h::Library{};
+  auto world = lib.world();
+
+  // Create a simple 2D box mesh: 1.0 x 1.0 domain with 4x4 elements
+  auto mesh =
+    Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 4, 4, 0, false);
+  auto grid = pcms::CreateUniformGridFromMesh<2>(mesh, {4, 4});
+  auto [mask_layout, mask_field] = pcms::CreateUniformGridBinaryField<2>(mesh, {4, 4});
+
+  // Create OmegaH field layout with linear elements
+  auto omega_h_layout =
+    pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
+  auto omega_h_field = omega_h_layout->CreateField();
+
+  // Initialize omega_h field with a simple function f(x,y) = x + 2*y
+  auto coords = omega_h_layout->GetDOFHolderCoordinates();
+  int num_nodes = omega_h_layout->GetNumOwnedDofHolder();
+  std::vector<pcms::Real> omega_h_data = CreateOmegaHFieldData(coords, num_nodes);
+
+  auto omega_h_data_view =
+    pcms::Rank1View<const pcms::Real, pcms::HostMemorySpace>(
+      omega_h_data.data(), omega_h_data.size());
+  omega_h_field->SetDOFHolderData(omega_h_data_view);
+
+  // Create uniform grid field layout
+  pcms::UniformGridFieldLayout<2> ug_layout(grid, 1,
+                                            pcms::CoordinateSystem::Cartesian);
+  auto ug_field = ug_layout.CreateField();
+
+  // Transfer from omega_h field to uniform grid field using interpolation
+  pcms::interpolate_field2(*omega_h_field, *ug_field);
+  auto ug_coords = ug_layout.GetDOFHolderCoordinates();
+
+  // Verify ug_field values directly from the field object
+  auto ug_field_data = ug_field->GetDOFHolderData();
+  VerifyUniformGridFieldValues(grid, ug_coords, ug_field_data);
+
+  // Verify mask field values
+  VerifyMaskFieldValues(grid, *mask_field);
+}