diff --git a/src/pcms/capi/client.cpp b/src/pcms/capi/client.cpp
index 3d320094..485e3c12 100644
--- a/src/pcms/capi/client.cpp
+++ b/src/pcms/capi/client.cpp
@@ -1,5 +1,6 @@
 #include "client.h"
 #include "pcms.h"
+#include "pcms/field/field_factory.h"
 #include "pcms/field/function_space/xgc.h"
 #include "pcms/field/data/xgc.h"
 #include "pcms/field/layout/xgc.h"
@@ -23,11 +24,11 @@ namespace detail
 template <typename T>
 struct XGCFieldRegistration
 {
-  XGCFunctionSpace function_space;
+  XGCFieldFactory function_space;
   MPI_Comm plane_comm;
   Rank1View<T, HostMemorySpace> data;
 
-  XGCFieldRegistration(XGCFunctionSpace fs, MPI_Comm comm,
+  XGCFieldRegistration(XGCFieldFactory fs, MPI_Comm comm,
                        Rank1View<T, HostMemorySpace> d)
     : function_space(std::move(fs)), plane_comm(comm), data(d)
   {
@@ -37,10 +38,15 @@ struct XGCFieldRegistration
 struct DummyFieldRegistration
 {};
 
-class EmptyFunctionSpace : public FunctionSpace
+class EmptyFieldFactory : public FieldFactory
 {
 public:
-  EmptyFunctionSpace() : layout_(std::make_shared<EmptyFieldLayout>()) {}
+  explicit EmptyFieldFactory(std::string layout_name = "")
+  {
+    auto layout = std::make_shared<EmptyFieldLayout>();
+    layout->SetName(std::move(layout_name));
+    layout_ = layout;
+  }
 
   [[nodiscard]] std::shared_ptr<const FieldLayout> GetLayout()
     const noexcept override
@@ -48,11 +54,6 @@ class EmptyFunctionSpace : public FunctionSpace
     return layout_;
   }
 
-  [[nodiscard]] CoordinateSystem GetCoordinateSystem() const noexcept override
-  {
-    return CoordinateSystem::Cartesian;
-  }
-
 protected:
   [[nodiscard]] FieldVariant CreateFieldImpl(
     Type value_type, FieldMetadata metadata) const override
@@ -75,12 +76,12 @@ class EmptyFunctionSpace : public FunctionSpace
         PCMS_ALWAYS_ASSERT(fd != nullptr);
         if (dynamic_cast<const SimpleFieldData<T>*>(fd.get()) == nullptr) {
           throw pcms_error(
-            "EmptyFunctionSpace::CreateField: requires SimpleFieldData");
+            "EmptyFieldFactory::CreateField: requires SimpleFieldData");
         }
         if (fd->GetDOFHolderDataHost().size() !=
             detail::ExpectedFlatFieldDataSize(*layout_)) {
           throw pcms_error(
-            "EmptyFunctionSpace::CreateField: field data size does not match "
+            "EmptyFieldFactory::CreateField: field data size does not match "
             "layout");
         }
         return WrapField<T>(layout_, std::forward<decltype(fd)>(fd));
@@ -88,12 +89,6 @@ class EmptyFunctionSpace : public FunctionSpace
       std::move(data));
   }
 
-  [[nodiscard]] PointEvaluatorVariant CreatePointEvaluatorImpl(
-    Type /*value_type*/, const EvaluationRequest& /*request*/) const override
-  {
-    throw pcms_error("EmptyFunctionSpace does not support point evaluation");
-  }
-
 private:
   std::shared_ptr<const EmptyFieldLayout> layout_;
 };
@@ -130,29 +125,31 @@ ClientState::HandleVariant RegisterField(
   Application& app, std::string name,
   const detail::XGCFieldRegistration<T>& registration, bool participates)
 {
+  // The layout's coupling identity is the field name
+  // not the arbitrary adapter name given at creation.
+  registration.function_space.SetLayoutName(name);
   auto field = registration.function_space.template CreateField<T>(
-    std::make_unique<XGCFieldData<T>>(
-      registration.function_space.GetXGCLayout(), FieldMetadata{},
-      registration.data));
-  app.AddLayout(name, registration.function_space.GetLayout(), participates);
+    std::move(name), std::make_unique<XGCFieldData<T>>(
+                       registration.function_space.GetXGCLayout(),
+                       FieldMetadata{}, registration.data));
   std::unique_ptr<FieldSerializer<T>> serializer =
     std::make_unique<XGCFieldSerializer<T>>(registration.plane_comm,
                                             participates);
-  return ClientState::HandleVariant{app.AddField(
-    std::move(name), std::move(field), std::move(serializer), participates)};
+  return ClientState::HandleVariant{
+    app.AddField(std::move(field), std::move(serializer), participates)};
 }
 
 inline ClientState::HandleVariant RegisterField(
   Application& app, std::string name, const detail::DummyFieldRegistration&,
   bool participates)
 {
-  auto function_space = detail::EmptyFunctionSpace{};
-  app.AddLayout(name, function_space.GetLayout(), participates);
-  auto field = function_space.CreateField<int>(FieldMetadata{});
+  auto function_space = detail::EmptyFieldFactory{name};
+  auto field =
+    function_space.CreateField<int>(std::move(name), FieldMetadata{});
   std::unique_ptr<FieldSerializer<int>> serializer =
     std::make_unique<FieldSerializer<int>>();
-  return ClientState::HandleVariant{app.AddField(
-    std::move(name), std::move(field), std::move(serializer), participates)};
+  return ClientState::HandleVariant{
+    app.AddField(std::move(field), std::move(serializer), participates)};
 }
 
 } // namespace pcms
@@ -251,7 +248,7 @@ void pcms_create_xgc_field_adapter_t(
 {
   PCMS_ALWAYS_ASSERT((size > 0) ? (data != nullptr) : true);
   auto function_space =
-    pcms::XGCFunctionSpace(reverse_classification, in_overlap, size);
+    pcms::XGCFieldFactory(reverse_classification, in_overlap, size);
   pcms::Rank1View<T, pcms::HostMemorySpace> data_view(
     reinterpret_cast<T*>(data), size);
   field_adapter.emplace<pcms::detail::XGCFieldRegistration<T>>(
diff --git a/src/pcms/coupler/coupler.cpp b/src/pcms/coupler/coupler.cpp
index a0c83e7d..91a1bf80 100644
--- a/src/pcms/coupler/coupler.cpp
+++ b/src/pcms/coupler/coupler.cpp
@@ -3,32 +3,23 @@
 namespace pcms
 {
 
-FieldLayoutCommunicator& Application::GetLayoutCommunicator(
-  const FieldLayout& layout)
+FieldLayoutCommunicator& Application::GetOrCreateLayoutCommunicator(
+  const FieldLayout& layout, bool participates)
 {
   PCMS_FUNCTION_TIMER;
-  auto it = field_layout_communicators_.find(&layout);
+  const std::string& name = layout.GetName();
+  if (name.empty()) {
+    throw pcms_error(
+      "Application: the field's layout has no name; name it at construction "
+      "(From*(..., layout_name)) before registering it with a coupler");
+  }
+
+  // Fields whose layouts share a name share one GID-exchange communicator.
+  auto it = field_layout_communicators_.find(name);
   if (it != field_layout_communicators_.end()) {
     return *it->second;
-  } else {
-    throw pcms_error("Field added with unregistered layout. Call AddLayout() "
-                     "before AddField().");
   }
-}
 
-const FieldLayout& Application::AddLayout(
-  std::string name, std::shared_ptr<const FieldLayout> layout,
-  bool participates)
-{
-  return AddLayout(std::move(name), std::move(layout),
-                   std::make_unique<GenericFieldExchangePlanner>(),
-                   participates);
-}
-
-const FieldLayout& Application::AddLayout(
-  std::string name, std::shared_ptr<const FieldLayout> layout,
-  std::unique_ptr<FieldExchangePlanner> planner, bool participates)
-{
   MPI_Comm mpi_comm_subset = MPI_COMM_NULL;
   bool own_mpi_comm = false;
   PCMS_ALWAYS_ASSERT((mpi_comm_ == MPI_COMM_NULL) ? (!participates) : true);
@@ -39,21 +30,19 @@ const FieldLayout& Application::AddLayout(
                    &mpi_comm_subset);
     own_mpi_comm = true;
   }
-  layouts_.push_back(std::move(layout));
-  const FieldLayout& layout_ref = *layouts_.back();
 
-  // Check if there's an overlap mask for this layout
   const OverlapMask* overlap_mask = nullptr;
   auto mask_it = layout_overlap_masks_.find(name);
   if (mask_it != layout_overlap_masks_.end()) {
     overlap_mask = mask_it->second.get();
   }
 
-  field_layout_communicators_.emplace(
-    &layout_ref, std::make_unique<FieldLayoutCommunicator>(
-                   name, mpi_comm_subset, redev_, channel_, layout_ref,
-                   std::move(planner), own_mpi_comm, overlap_mask));
-  return layout_ref;
+  auto [it2, inserted] = field_layout_communicators_.emplace(
+    name, std::make_unique<FieldLayoutCommunicator>(
+            name, mpi_comm_subset, redev_, channel_, layout,
+            std::make_unique<GenericFieldExchangePlanner>(), own_mpi_comm,
+            overlap_mask));
+  return *it2->second;
 }
 
 void Application::SetLayoutOverlapMask(
diff --git a/src/pcms/coupler/coupler.hpp b/src/pcms/coupler/coupler.hpp
index d2d12cb1..eef0f564 100644
--- a/src/pcms/coupler/coupler.hpp
+++ b/src/pcms/coupler/coupler.hpp
@@ -31,11 +31,41 @@ class FieldHandle
   void Receive(redev::Mode mode = redev::Mode::Synchronous) const;
   [[nodiscard]] Field<T>& GetField() const;
 
+  [[nodiscard]] const std::string& GetName() const noexcept { return name_; }
+
 private:
   Application* app_;
   std::string name_;
 };
 
+// FunctionHandle is a FieldHandle that additionally carries the function space,
+// so it can be evaluated / used to build transfer operators. Only AddFunction
+// produces one; CreateTransfer accepts FunctionHandle (not FieldHandle), which
+// makes passing a comm-only field to a transfer a compile error.
+template <typename T>
+class FunctionHandle : public FieldHandle<T>
+{
+public:
+  FunctionHandle(Application* app, std::string name,
+                 std::shared_ptr<const FunctionSpace> space)
+    : FieldHandle<T>(app, std::move(name)), space_(std::move(space))
+  {
+  }
+
+  [[nodiscard]] const FunctionSpace& GetSpace() const noexcept
+  {
+    return *space_;
+  }
+  [[nodiscard]] std::shared_ptr<const FunctionSpace> GetSpacePtr()
+    const noexcept
+  {
+    return space_;
+  }
+
+private:
+  std::shared_ptr<const FunctionSpace> space_;
+};
+
 class Application
 {
 public:
@@ -50,27 +80,33 @@ class Application
     PCMS_FUNCTION_TIMER;
   }
 
-  const FieldLayout& AddLayout(std::string name,
-                               std::shared_ptr<const FieldLayout> layout,
-                               bool participates = true);
-  const FieldLayout& AddLayout(std::string name,
-                               std::shared_ptr<const FieldLayout> layout,
-                               std::unique_ptr<FieldExchangePlanner> planner,
-                               bool participates = true);
-
-  // Set the overlap mask for a specific layout by name
+  // Set the overlap mask for a field's layout by name. Must be set before the
+  // AddField/AddFunction that first creates that layout's communicator.
   void SetLayoutOverlapMask(const std::string& layout_name,
                             std::unique_ptr<OverlapMask> overlap_mask);
 
+  // Register a comm-only field, identified by its own name. The layout's
+  // communicator is created on first use and shared by any later
+  // field/function whose layout has the same name.
   template <typename T>
-  FieldHandle<T> AddField(std::string name, Field<T>&& field,
-                          bool participates = true);
+  FieldHandle<T> AddField(Field<T>&& field, bool participates = true);
 
   template <typename T>
-  FieldHandle<T> AddField(std::string name, Field<T>&& field,
+  FieldHandle<T> AddField(Field<T>&& field,
                           std::unique_ptr<FieldSerializer<T>> serializer,
                           bool participates = true);
 
+  // Register a transferable field: like AddField but retains the function space
+  // (via the stored Function) and returns a FunctionHandle usable in transfers.
+  template <typename T>
+  FunctionHandle<T> AddFunction(Function<T>&& function,
+                                bool participates = true);
+
+  template <typename T>
+  FunctionHandle<T> AddFunction(Function<T>&& function,
+                                std::unique_ptr<FieldSerializer<T>> serializer,
+                                bool participates = true);
+
   void SendField(const std::string& name,
                  redev::Mode mode = redev::Mode::Synchronous)
   {
@@ -144,20 +180,29 @@ class Application
   [[nodiscard]] Field<T>& GetField(const std::string& name);
 
 private:
-  FieldLayoutCommunicator& GetLayoutCommunicator(const FieldLayout& layout);
+  // Returns the communicator for `layout`, creating it (MPI split + overlap
+  // mask lookup + planner) on first use and reusing it for any later
+  // field/function whose layout has the same name.
+  FieldLayoutCommunicator& GetOrCreateLayoutCommunicator(
+    const FieldLayout& layout, bool participates);
+
+  template <typename T>
+  void RegisterFieldCommunicator(Field<T>& field_obj,
+                                 std::unique_ptr<FieldSerializer<T>> serializer,
+                                 bool participates);
 
   MPI_Comm mpi_comm_;
   redev::Redev& redev_;
   redev::Channel channel_;
-  std::vector<std::shared_ptr<const FieldLayout>> layouts_;
   std::map<std::string, FieldVariant> fields_;
+  std::map<std::string, FunctionVariant> functions_;
   std::vector<std::unique_ptr<FieldLayoutCommunicator>>
     owned_field_layout_communicators_;
   // map is used rather than unordered_map because we give pointers to the
   // internal data and rehash of unordered_map can cause pointer invalidation.
   // map is less cache friendly, but pointers are not invalidated.
   std::map<std::string, FieldCommunicator2Ptr> field_communicators_;
-  std::map<const FieldLayout*, std::unique_ptr<FieldLayoutCommunicator>>
+  std::map<std::string, std::unique_ptr<FieldLayoutCommunicator>>
     field_layout_communicators_;
   std::map<std::string, std::unique_ptr<OverlapMask>> layout_overlap_masks_;
 };
@@ -238,47 +283,106 @@ pcms::Field<T>& pcms::FieldHandle<T>::GetField() const
 template <typename T>
 pcms::Field<T>& pcms::Application::GetField(const std::string& name)
 {
-  auto* field = std::get_if<Field<T>>(&detail::find_or_error(name, fields_));
-  if (field == nullptr) {
-    throw pcms_error("Field stored with different type than requested");
+  auto field_it = fields_.find(name);
+  if (field_it != fields_.end()) {
+    auto* field = std::get_if<Field<T>>(&field_it->second);
+    if (field == nullptr) {
+      throw pcms_error("Field stored with different type than requested");
+    }
+    return *field;
+  }
+  auto fn_it = functions_.find(name);
+  if (fn_it != functions_.end()) {
+    auto* function = std::get_if<Function<T>>(&fn_it->second);
+    if (function == nullptr) {
+      throw pcms_error("Field stored with different type than requested");
+    }
+    return *function;
+  }
+  throw pcms_error("Field '" + name + "' not found");
+}
+
+template <typename T>
+void pcms::Application::RegisterFieldCommunicator(
+  Field<T>& field_obj, std::unique_ptr<FieldSerializer<T>> serializer,
+  bool participates)
+{
+  const std::string& name = field_obj.GetName();
+  FieldLayoutCommunicator& layout_communicator =
+    GetOrCreateLayoutCommunicator(field_obj.GetLayout(), participates);
+  FieldCommunicator2Ptr field_communicator =
+    std::make_unique<FieldCommunicator<T>>(name, layout_communicator, field_obj,
+                                           std::move(serializer));
+  auto [it, inserted] =
+    field_communicators_.emplace(name, std::move(field_communicator));
+  if (!inserted) {
+    throw pcms_error("Field with this name already exists");
   }
-  return *field;
 }
 
 template <typename T>
-pcms::FieldHandle<T> pcms::Application::AddField(std::string name,
-                                                 Field<T>&& field,
+pcms::FieldHandle<T> pcms::Application::AddField(Field<T>&& field,
                                                  bool participates)
 {
-  return AddField(std::move(name), std::move(field),
-                  std::make_unique<FieldSerializer<T>>(), participates);
+  return AddField(std::move(field), std::make_unique<FieldSerializer<T>>(),
+                  participates);
 }
 
 template <typename T>
 pcms::FieldHandle<T> pcms::Application::AddField(
-  std::string name, Field<T>&& field,
-  std::unique_ptr<FieldSerializer<T>> serializer, bool participates)
+  Field<T>&& field, std::unique_ptr<FieldSerializer<T>> serializer,
+  bool participates)
 {
   PCMS_FUNCTION_TIMER;
-  (void)participates;
+  std::string name = field.GetName();
+  if (name.empty()) {
+    throw pcms_error("Application::AddField: field has no name");
+  }
   auto [field_it, field_inserted] = fields_.emplace(name, std::move(field));
   if (!field_inserted) {
     throw pcms_error("Field with this name already exists");
   }
-  auto& field_obj = std::get<Field<T>>(field_it->second);
-  const FieldLayout& layout = field_obj.GetLayout();
-  FieldLayoutCommunicator& layout_communicator = GetLayoutCommunicator(layout);
-  FieldCommunicator2Ptr field_communicator =
-    std::make_unique<FieldCommunicator<T>>(name, layout_communicator, field_obj,
-                                           std::move(serializer));
-
-  auto [it, inserted] =
-    field_communicators_.emplace(name, std::move(field_communicator));
-  if (!inserted) {
+  try {
+    RegisterFieldCommunicator<T>(std::get<Field<T>>(field_it->second),
+                                 std::move(serializer), participates);
+  } catch (...) {
     fields_.erase(field_it);
-    throw pcms_error("Field with this name already exists");
+    throw;
   }
   return FieldHandle<T>{this, std::move(name)};
 }
 
+template <typename T>
+pcms::FunctionHandle<T> pcms::Application::AddFunction(Function<T>&& function,
+                                                       bool participates)
+{
+  return AddFunction(std::move(function),
+                     std::make_unique<FieldSerializer<T>>(), participates);
+}
+
+template <typename T>
+pcms::FunctionHandle<T> pcms::Application::AddFunction(
+  Function<T>&& function, std::unique_ptr<FieldSerializer<T>> serializer,
+  bool participates)
+{
+  PCMS_FUNCTION_TIMER;
+  std::string name = function.GetName();
+  if (name.empty()) {
+    throw pcms_error("Application::AddFunction: function has no name");
+  }
+  auto space = function.GetSpacePtr();
+  auto [fn_it, fn_inserted] = functions_.emplace(name, std::move(function));
+  if (!fn_inserted) {
+    throw pcms_error("Field with this name already exists");
+  }
+  try {
+    RegisterFieldCommunicator<T>(std::get<Function<T>>(fn_it->second),
+                                 std::move(serializer), participates);
+  } catch (...) {
+    functions_.erase(fn_it);
+    throw;
+  }
+  return FunctionHandle<T>{this, std::move(name), std::move(space)};
+}
+
 #endif // COUPLER2_H_
diff --git a/src/pcms/coupler/field_communicator.hpp b/src/pcms/coupler/field_communicator.hpp
index b1874507..fa731487 100644
--- a/src/pcms/coupler/field_communicator.hpp
+++ b/src/pcms/coupler/field_communicator.hpp
@@ -34,7 +34,8 @@ class FieldCommunicator
       field_(field),
       serializer_(std::move(serializer))
   {
-    PCMS_ALWAYS_ASSERT(&layout_comm.GetLayout() == &field.GetLayout());
+    PCMS_ALWAYS_ASSERT(layout_comm.GetLayout().GetName() ==
+                       field.GetLayout().GetName());
     // The exchange plan is per DOF holder; the field payload carries
     // num_components values per holder, so the buffer is scaled accordingly.
     comm_buffer_.resize(
diff --git a/src/pcms/field/eqdsk_field.h b/src/pcms/field/eqdsk_field.h
index bb4188c8..05e8308d 100644
--- a/src/pcms/field/eqdsk_field.h
+++ b/src/pcms/field/eqdsk_field.h
@@ -12,7 +12,7 @@ namespace pcms
 
 struct EQDSKField
 {
-  SplineFunctionSpace space;
+  std::shared_ptr<SplineFunctionSpace> space;
   Field<Real> field;
 };
 
@@ -20,7 +20,7 @@ struct EQDSKField
 struct EQDSKFieldWithData
 {
   EQDSKData data;
-  SplineFunctionSpace space;
+  std::shared_ptr<SplineFunctionSpace> space;
   Field<Real> field;
 };
 
@@ -43,7 +43,7 @@ inline EQDSKField MakeEQDSKField(const EQDSKData& eqdsk_data)
 {
   auto space = SplineFunctionSpace::FromUniformGrid(
     eqdsk_data.grid, CoordinateSystem::Cartesian);
-  auto field = space.CreateField<Real>();
+  Field<Real> field = space->CreateFunction<Real>();
   field.GetData().SetDOFHolderData(Rank2View<const Real, DeviceMemorySpace>(
     eqdsk_data.PSIZR.data(), eqdsk_data.PSIZR.extent(0), 1));
   return {std::move(space), std::move(field)};
diff --git a/src/pcms/field/field.h b/src/pcms/field/field.h
index 433717d0..6da4261d 100644
--- a/src/pcms/field/field.h
+++ b/src/pcms/field/field.h
@@ -2,12 +2,13 @@
 #define PCMS_COUPLING_FIELD_H
 
 #include "field_data.h"
-#include "field_evaluator_factory.h"
 #include "field_layout.h"
 #include "pcms/utility/arrays.h"
 #include "pcms/utility/memory_spaces.h"
 #include "pcms/utility/types.h"
 #include <memory>
+#include <string>
+#include <utility>
 #include <variant>
 
 namespace pcms
@@ -15,12 +16,13 @@ namespace pcms
 
 class FunctionSpace;
 
-// Field<T> is a composed per-field object: it owns coefficient data and holds
-// a shared reference to the evaluator factory so the function space stays alive
-// as long as the field does.
+// Field<T> is the (layout + data) bundle: it pairs coefficient data with the
+// FieldLayout that gives the DOF-holder ordering its meaning. It is the
+// currency of the communication path and of TransferOperator::Apply; it is NOT
+// evaluatable on its own (evaluation semantics live on FunctionSpace).
 //
-// Fields are typically created via LagrangeFunctionSpace::CreateField().
-// They are move-only (unique_ptr member).
+// Fields are created via FunctionSpace::CreateField(). They are move-only
+// (unique_ptr member).
 template <typename T>
 class Field
 {
@@ -29,12 +31,16 @@ class Field
   Field& operator=(Field&&) = default;
   Field(const Field&) = delete;
   Field& operator=(const Field&) = delete;
+  ~Field() = default;
 
   FieldData<T>& GetData() noexcept { return *data_; }
   const FieldData<T>& GetData() const noexcept { return *data_; }
 
   const FieldLayout& GetLayout() const { return *layout_; }
 
+  // Optional identifier for this field. Empty by default; set at creation.
+  [[nodiscard]] const std::string& GetName() const noexcept { return name_; }
+
   Rank2View<const T, HostMemorySpace> GetDOFHolderDataHost() const
   {
     return data_->GetDOFHolderDataHost();
@@ -55,32 +61,69 @@ class Field
     data_->SetDOFHolderData(v);
   }
 
-private:
-  class CtorKey
-  {
-    CtorKey() = default;
-    friend class FunctionSpace;
-  };
-
-  Field(CtorKey, std::shared_ptr<const FieldLayout> layout,
-        std::shared_ptr<const FieldEvaluatorFactory<Real>> evaluator_factory,
+protected:
+  Field(std::string name, std::shared_ptr<const FieldLayout> layout,
         std::unique_ptr<FieldData<T>> data)
-    : layout_(std::move(layout)),
-      evaluator_factory_(std::move(evaluator_factory)),
-      data_(std::move(data))
+    : name_(std::move(name)), layout_(std::move(layout)), data_(std::move(data))
   {
   }
 
+  // FieldFactory constructs Fields (via WrapField); FunctionSpace constructs
+  // Fields the same way and moves a Field's name, layout and data out to build
+  // a Function (via CreateFunction). Both stamp the field name.
+  friend class FieldFactory;
   friend class FunctionSpace;
 
+private:
+  std::string name_;
   std::shared_ptr<const FieldLayout> layout_;
-  std::shared_ptr<const FieldEvaluatorFactory<Real>> evaluator_factory_;
   std::unique_ptr<FieldData<T>> data_;
 };
 
+// Function<T> is a Field that additionally knows its FunctionSpace, so it can
+// be evaluated and used to build transfer operators. It is-a Field so the
+// communication path (FieldCommunicator holds Field&) and
+// TransferOperator::Apply(const Field&) bind a Function slice-free. It adds no
+// virtuals: moving a Function into a Field slot (e.g. AddField(Field&&))
+// intentionally slices off the space reference, leaving a valid comm-only
+// Field.
+template <typename T>
+class Function : public Field<T>
+{
+public:
+  Function(Function&&) = default;
+  Function& operator=(Function&&) = default;
+  Function(const Function&) = delete;
+  Function& operator=(const Function&) = delete;
+
+  const FunctionSpace& GetSpace() const noexcept { return *space_; }
+  std::shared_ptr<const FunctionSpace> GetSpacePtr() const noexcept
+  {
+    return space_;
+  }
+
+protected:
+  Function(std::string name, std::shared_ptr<const FieldLayout> layout,
+           std::unique_ptr<FieldData<T>> data,
+           std::shared_ptr<const FunctionSpace> space)
+    : Field<T>(std::move(name), std::move(layout), std::move(data)),
+      space_(std::move(space))
+  {
+  }
+
+  friend class FunctionSpace;
+
+private:
+  std::shared_ptr<const FunctionSpace> space_;
+};
+
 using FieldVariant = std::variant<Field<int8_t>, Field<int32_t>, Field<int64_t>,
                                   Field<float>, Field<double>>;
 
+using FunctionVariant =
+  std::variant<Function<int8_t>, Function<int32_t>, Function<int64_t>,
+               Function<float>, Function<double>>;
+
 } // namespace pcms
 
 #endif // PCMS_COUPLING_FIELD_H
diff --git a/src/pcms/field/field_factory.h b/src/pcms/field/field_factory.h
new file mode 100644
index 00000000..7903c72e
--- /dev/null
+++ b/src/pcms/field/field_factory.h
@@ -0,0 +1,102 @@
+#ifndef PCMS_FIELD_FACTORY_H
+#define PCMS_FIELD_FACTORY_H
+
+#include "field.h"
+#include "field_data.h"
+#include "field_layout.h"
+#include "field_metadata.h"
+#include "pcms/utility/types.h"
+#include <memory>
+#include <string>
+
+namespace pcms
+{
+
+namespace detail
+{
+
+inline size_t ExpectedFlatFieldDataSize(const FieldLayout& layout)
+{
+  return static_cast<size_t>(layout.GetNumOwnedDofHolder()) *
+         static_cast<size_t>(layout.GetNumComponents());
+}
+
+} // namespace detail
+
+// Compile-time gate: true only for the five supported field value types.
+template <typename T>
+inline constexpr bool is_supported_field_type_v =
+  std::is_same_v<T, int8_t> || std::is_same_v<T, int32_t> ||
+  std::is_same_v<T, int64_t> || std::is_same_v<T, float> ||
+  std::is_same_v<T, double>;
+
+// FieldFactory constructs Field<T> bundles ({layout, data}) with no evaluation
+// capability and no shared-ownership requirement. Comm-only backends (e.g.
+// XGCFieldFactory) derive it directly. FunctionSpace is a separate abstraction
+// (not a FieldFactory) that additionally supports evaluation and produces
+// Functions; a FieldFactory can never be passed where a FunctionSpace is
+// required, which is the compile-time form of "comm-only cannot be evaluated".
+class FieldFactory
+{
+public:
+  virtual std::shared_ptr<const FieldLayout> GetLayout() const noexcept = 0;
+
+  virtual ~FieldFactory() noexcept = default;
+
+  // Create a new named field with freshly allocated data. The name is optional
+  // (empty by default) and identifies the field to consumers such as a coupler.
+  template <typename T>
+  [[nodiscard]] Field<T> CreateField(std::string name = "",
+                                     FieldMetadata metadata = {}) const;
+
+  // Expert API: wrap externally constructed field data into a Field. The
+  // concrete factory validates backend-specific field-data type and storage
+  // size compatibility.
+  template <typename T>
+  [[nodiscard]] Field<T> CreateField(std::string name,
+                                     std::unique_ptr<FieldData<T>> data) const;
+
+protected:
+  template <typename T>
+  static Field<T> WrapField(std::shared_ptr<const FieldLayout> layout,
+                            std::unique_ptr<FieldData<T>> data)
+  {
+    return Field<T>(std::string{}, std::move(layout), std::move(data));
+  }
+
+  virtual FieldVariant CreateFieldImpl(Type value_type,
+                                       FieldMetadata metadata) const = 0;
+
+  virtual FieldVariant CreateFieldImpl(FieldDataVariant data) const = 0;
+};
+
+template <typename T>
+Field<T> FieldFactory::CreateField(std::string name,
+                                   FieldMetadata metadata) const
+{
+  static_assert(is_supported_field_type_v<T>,
+                "T is not a supported field type");
+  Field<T> field =
+    std::get<Field<T>>(CreateFieldImpl(TypeEnumFromType<T>(), metadata));
+  field.name_ = std::move(name);
+  return field;
+}
+
+template <typename T>
+Field<T> FieldFactory::CreateField(std::string name,
+                                   std::unique_ptr<FieldData<T>> data) const
+{
+  static_assert(is_supported_field_type_v<T>,
+                "T is not a supported field type");
+  if (!data) {
+    throw pcms_error("FieldFactory::CreateField: data must not be null");
+  }
+  Field<T> field =
+    std::get<Field<T>>(CreateFieldImpl(FieldDataVariant{std::move(data)}));
+  field.name_ = std::move(name);
+  return field;
+}
+
+} // namespace pcms
+
+#endif // PCMS_FIELD_FACTORY_H
diff --git a/src/pcms/field/field_layout.h b/src/pcms/field/field_layout.h
index d01a13a0..b4551e4c 100644
--- a/src/pcms/field/field_layout.h
+++ b/src/pcms/field/field_layout.h
@@ -2,6 +2,7 @@
 #define PCMS_FIELD_LAYOUT_H
 #include <map>
 #include <memory>
+#include <string>
 #include <vector>
 #include "pcms/discretization/discretization.h"
 #include "pcms/utility/types.h"
@@ -19,6 +20,12 @@ using ReversePartitionMap = std::map<pcms::LO, std::vector<pcms::LO>>;
 class FieldLayout
 {
 public:
+  // Optional identifier for this layout. Set once at construction (via a From*
+  // factory) and empty by default. Consumers that need to identify a layout by
+  // a stable name use it; it carries no meaning within the field layer itself.
+  [[nodiscard]] const std::string& GetName() const noexcept { return name_; }
+  void SetName(std::string name) { name_ = std::move(name); }
+
   virtual std::shared_ptr<const Discretization> GetDiscretization()
     const noexcept = 0;
 
@@ -77,6 +84,7 @@ class FieldLayout
   void BuildGlobalToLocalPermutation();
 
 private:
+  std::string name_;
   Kokkos::View<LO*, HostMemorySpace> global_to_local_host_;
   Kokkos::View<LO*, DeviceMemorySpace> global_to_local_;
 };
diff --git a/src/pcms/field/function_space.h b/src/pcms/field/function_space.h
index a5e6f867..26ef4fad 100644
--- a/src/pcms/field/function_space.h
+++ b/src/pcms/field/function_space.h
@@ -5,7 +5,7 @@
 #include "evaluation_request.h"
 #include "field.h"
 #include "field_data.h"
-#include "field_evaluator_factory.h"
+#include "field_factory.h"
 #include "field_layout.h"
 #include "field_metadata.h"
 #include "out_of_bounds_policy.h"
@@ -20,33 +20,18 @@
 namespace pcms
 {
 
-namespace detail
-{
-
-inline size_t ExpectedFlatFieldDataSize(const FieldLayout& layout)
-{
-  return static_cast<size_t>(layout.GetNumOwnedDofHolder()) *
-         static_cast<size_t>(layout.GetNumComponents());
-}
-
-} // namespace detail
-
-// Compile-time gate: true only for the five supported field value types.
-template <typename T>
-inline constexpr bool is_supported_field_type_v =
-  std::is_same_v<T, int8_t> || std::is_same_v<T, int32_t> ||
-  std::is_same_v<T, int64_t> || std::is_same_v<T, float> ||
-  std::is_same_v<T, double>;
-
 // FunctionSpace is an abstract interface representing an evaluatable field
-// space: layout, evaluation rules, and coordinate interpretation.
+// space: layout, evaluation rules, and coordinate interpretation. It produces
+// Functions (a Field bound to this space).
 //
 // Concrete implementations (e.g. LagrangeFunctionSpace) provide backends for
-// specific discretizations or mesh types.
-//
-// FunctionSpace is used as the parameter type for operation objects such as
-// Interpolator<T>, so that operations are not coupled to a specific backend.
-class FunctionSpace
+// specific discretizations or mesh types, and their From* factories return a
+// shared_ptr. A Function co-owns its space via shared_ptr (obtained through
+// enable_shared_from_this), which is why concrete spaces are always
+// shared-owned; a stack-allocated FunctionSpace is not constructible (the
+// backend constructors are private and reachable only via the shared_ptr
+// From* factories).
+class FunctionSpace : public std::enable_shared_from_this<FunctionSpace>
 {
 public:
   virtual std::shared_ptr<const Discretization> GetDiscretization()
@@ -61,17 +46,16 @@ class FunctionSpace
 
   virtual ~FunctionSpace() noexcept = default;
 
-  // Create a new field with freshly allocated data for this function space.
-  // Compile-time error for unsupported T; runtime error for T unsupported by
-  // the concrete backend.
   template <typename T>
-  [[nodiscard]] Field<T> CreateField(FieldMetadata metadata = {}) const;
+  [[nodiscard]] Function<T> CreateFunction(std::string name = "",
+                                           FieldMetadata metadata = {}) const;
 
-  // Expert API: wrap externally constructed field data into a Field for this
-  // function space. The concrete function space validates backend-specific
-  // field-data type and storage size compatibility.
+  // Expert API: wrap externally constructed field data into a Function for this
+  // space. The concrete space validates backend-specific field-data type and
+  // storage size compatibility.
   template <typename T>
-  [[nodiscard]] Field<T> CreateField(std::unique_ptr<FieldData<T>> data) const;
+  [[nodiscard]] Function<T> CreateFunction(
+    std::string name, std::unique_ptr<FieldData<T>> data) const;
 
   // Create a point evaluator for the given evaluation request.
   // Compile-time error for unsupported T; runtime error for T or capability
@@ -89,12 +73,19 @@ class FunctionSpace
 protected:
   template <typename T>
   static Field<T> WrapField(std::shared_ptr<const FieldLayout> layout,
-                            std::unique_ptr<FieldData<T>> data,
-                            std::shared_ptr<const FieldEvaluatorFactory<Real>>
-                              evaluator_factory = nullptr)
+                            std::unique_ptr<FieldData<T>> data)
+  {
+    return Field<T>(std::string{}, std::move(layout), std::move(data));
+  }
+
+  template <typename T>
+  static Function<T> WrapFunction(std::string name,
+                                  std::shared_ptr<const FieldLayout> layout,
+                                  std::unique_ptr<FieldData<T>> data,
+                                  std::shared_ptr<const FunctionSpace> space)
   {
-    return Field<T>(typename Field<T>::CtorKey{}, std::move(layout),
-                    std::move(evaluator_factory), std::move(data));
+    return Function<T>(std::move(name), std::move(layout), std::move(data),
+                       std::move(space));
   }
 
   virtual FieldVariant CreateFieldImpl(Type value_type,
@@ -107,22 +98,30 @@ class FunctionSpace
 };
 
 template <typename T>
-Field<T> FunctionSpace::CreateField(FieldMetadata metadata) const
+Function<T> FunctionSpace::CreateFunction(std::string name,
+                                          FieldMetadata metadata) const
 {
   static_assert(is_supported_field_type_v<T>,
                 "T is not a supported field type");
-  return std::get<Field<T>>(CreateFieldImpl(TypeEnumFromType<T>(), metadata));
+  Field<T> field =
+    std::get<Field<T>>(CreateFieldImpl(TypeEnumFromType<T>(), metadata));
+  return WrapFunction<T>(std::move(name), std::move(field.layout_),
+                         std::move(field.data_), shared_from_this());
 }
 
 template <typename T>
-Field<T> FunctionSpace::CreateField(std::unique_ptr<FieldData<T>> data) const
+Function<T> FunctionSpace::CreateFunction(
+  std::string name, std::unique_ptr<FieldData<T>> data) const
 {
   static_assert(is_supported_field_type_v<T>,
                 "T is not a supported field type");
   if (!data) {
-    throw pcms_error("FunctionSpace::CreateField: data must not be null");
+    throw pcms_error("FunctionSpace::CreateFunction: data must not be null");
   }
-  return std::get<Field<T>>(CreateFieldImpl(FieldDataVariant{std::move(data)}));
+  Field<T> field =
+    std::get<Field<T>>(CreateFieldImpl(FieldDataVariant{std::move(data)}));
+  return WrapFunction<T>(std::move(name), std::move(field.layout_),
+                         std::move(field.data_), shared_from_this());
 }
 
 template <typename T>
diff --git a/src/pcms/field/function_space/lagrange.cpp b/src/pcms/field/function_space/lagrange.cpp
index 8328e266..8d15a6ec 100644
--- a/src/pcms/field/function_space/lagrange.cpp
+++ b/src/pcms/field/function_space/lagrange.cpp
@@ -54,7 +54,7 @@ void ValidateLagrangeWrappedFieldData(const FieldLayout& layout,
 } // namespace
 
 LagrangeFunctionSpace::LagrangeFunctionSpace(
-  std::shared_ptr<const FieldLayout> layout,
+  Key, std::shared_ptr<const FieldLayout> layout,
   std::function<FieldDataVariant(Type, FieldMetadata)> create_field_data_fn,
   std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory) noexcept
   : layout_(std::move(layout)),
@@ -63,10 +63,10 @@ LagrangeFunctionSpace::LagrangeFunctionSpace(
 {
 }
 
-LagrangeFunctionSpace LagrangeFunctionSpace::FromMesh(
+std::shared_ptr<LagrangeFunctionSpace> LagrangeFunctionSpace::FromMesh(
   Omega_h::Mesh& mesh, int order, int num_components,
   CoordinateSystem coordinate_system, std::string global_id_name,
-  Backend backend)
+  Backend backend, std::string layout_name)
 {
   // https://github.com/SCOREC/meshFields/issues/88
   if (backend == Backend::MeshFields && order == 0) {
@@ -89,10 +89,11 @@ LagrangeFunctionSpace LagrangeFunctionSpace::FromMesh(
     auto mesh_layout = std::make_shared<MeshFieldsAdapterLayout>(
       mesh, nodes_per_dim, num_components, coordinate_system,
       std::move(global_id_name));
+    mesh_layout->SetName(layout_name);
     auto eval_factory =
       std::make_shared<MeshFieldsEvaluatorFactory<Real>>(mesh_layout);
-    return LagrangeFunctionSpace(
-      mesh_layout,
+    return std::make_shared<LagrangeFunctionSpace>(
+      Key{}, mesh_layout,
       [mesh_layout](Type t, FieldMetadata metadata) -> FieldDataVariant {
         if (t == Type::Float) {
           if constexpr (std::is_same_v<MeshField::Real4, float> ||
@@ -126,10 +127,11 @@ LagrangeFunctionSpace LagrangeFunctionSpace::FromMesh(
   }
   auto layout = std::make_shared<OmegaHLagrangeLayout>(
     mesh, order, num_components, coordinate_system, std::move(global_id_name));
+  layout->SetName(std::move(layout_name));
   auto eval_factory =
     std::make_shared<OmegaHLagrangeEvaluatorFactory<Real>>(layout);
-  return LagrangeFunctionSpace(
-    layout,
+  return std::make_shared<LagrangeFunctionSpace>(
+    Key{}, layout,
     [layout](Type t, FieldMetadata metadata) -> FieldDataVariant {
       return apply_to_type(t, [&](auto tag) -> FieldDataVariant {
         using T = typename decltype(tag)::type;
@@ -143,10 +145,10 @@ LagrangeFunctionSpace LagrangeFunctionSpace::FromMesh(
     std::move(eval_factory));
 }
 
-LagrangeFunctionSpace LagrangeFunctionSpace::FromMesh(
+std::shared_ptr<LagrangeFunctionSpace> LagrangeFunctionSpace::FromMesh(
   Omega_h::Mesh& mesh, int order, int num_components,
   CoordinateSystem coordinate_system, Omega_h::Read<Omega_h::I8> owned_mask,
-  std::string global_id_name, Backend backend)
+  std::string global_id_name, Backend backend, std::string layout_name)
 {
   if (backend == Backend::MeshFields) {
     throw pcms_error(
@@ -161,10 +163,11 @@ LagrangeFunctionSpace LagrangeFunctionSpace::FromMesh(
   auto layout = std::make_shared<OmegaHLagrangeLayout>(
     mesh, order, num_components, coordinate_system, std::move(owned_mask),
     std::move(global_id_name));
+  layout->SetName(std::move(layout_name));
   auto eval_factory =
     std::make_shared<OmegaHLagrangeEvaluatorFactory<Real>>(layout);
-  return LagrangeFunctionSpace(
-    layout,
+  return std::make_shared<LagrangeFunctionSpace>(
+    Key{}, layout,
     [layout](Type t, FieldMetadata metadata) -> FieldDataVariant {
       return apply_to_type(t, [&](auto tag) -> FieldDataVariant {
         using T = typename decltype(tag)::type;
@@ -178,9 +181,9 @@ LagrangeFunctionSpace LagrangeFunctionSpace::FromMesh(
     std::move(eval_factory));
 }
 
-LagrangeFunctionSpace LagrangeFunctionSpace::FromUniformGrid(
+std::shared_ptr<LagrangeFunctionSpace> LagrangeFunctionSpace::FromUniformGrid(
   const UniformGrid<2>& grid, int num_components,
-  CoordinateSystem coordinate_system, int order)
+  CoordinateSystem coordinate_system, int order, std::string layout_name)
 {
   if (order != 0 && order != 1) {
     throw std::invalid_argument("LagrangeFunctionSpace::FromUniformGrid: only "
@@ -188,10 +191,11 @@ LagrangeFunctionSpace LagrangeFunctionSpace::FromUniformGrid(
   }
   auto ug_layout = std::make_shared<UniformGridFieldLayout<2>>(
     grid, num_components, coordinate_system, order);
+  ug_layout->SetName(std::move(layout_name));
   auto eval_factory =
     std::make_shared<UniformGridEvaluatorFactory<2>>(ug_layout);
-  return LagrangeFunctionSpace(
-    ug_layout,
+  return std::make_shared<LagrangeFunctionSpace>(
+    Key{}, ug_layout,
     [ug_layout](Type t, FieldMetadata metadata) -> FieldDataVariant {
       return apply_to_type(t, [&](auto tag) -> FieldDataVariant {
         using T = typename decltype(tag)::type;
@@ -205,9 +209,9 @@ LagrangeFunctionSpace LagrangeFunctionSpace::FromUniformGrid(
     std::move(eval_factory));
 }
 
-LagrangeFunctionSpace LagrangeFunctionSpace::FromUniformGrid(
+std::shared_ptr<LagrangeFunctionSpace> LagrangeFunctionSpace::FromUniformGrid(
   const UniformGrid<3>& grid, int num_components,
-  CoordinateSystem coordinate_system, int order)
+  CoordinateSystem coordinate_system, int order, std::string layout_name)
 {
   if (order != 0 && order != 1) {
     throw std::invalid_argument("LagrangeFunctionSpace::FromUniformGrid: only "
@@ -215,10 +219,11 @@ LagrangeFunctionSpace LagrangeFunctionSpace::FromUniformGrid(
   }
   auto ug_layout = std::make_shared<UniformGridFieldLayout<3>>(
     grid, num_components, coordinate_system, order);
+  ug_layout->SetName(std::move(layout_name));
   auto eval_factory =
     std::make_shared<UniformGridEvaluatorFactory<3>>(ug_layout);
-  return LagrangeFunctionSpace(
-    ug_layout,
+  return std::make_shared<LagrangeFunctionSpace>(
+    Key{}, ug_layout,
     [ug_layout](Type t, FieldMetadata metadata) -> FieldDataVariant {
       return apply_to_type(t, [&](auto tag) -> FieldDataVariant {
         using T = typename decltype(tag)::type;
@@ -251,8 +256,7 @@ FieldVariant LagrangeFunctionSpace::CreateFieldImpl(
     [this](auto&& fd) -> FieldVariant {
       using FD = std::decay_t<decltype(fd)>;
       using T = typename FD::element_type::value_type;
-      return WrapField<T>(layout_, std::forward<decltype(fd)>(fd),
-                          evaluator_factory_);
+      return WrapField<T>(layout_, std::forward<decltype(fd)>(fd));
     },
     std::move(field_data));
 }
@@ -269,8 +273,7 @@ FieldVariant LagrangeFunctionSpace::CreateFieldImpl(FieldDataVariant data) const
           "LagrangeFunctionSpace: int8_t is not a supported field type");
       } else {
         ValidateLagrangeWrappedFieldData(*layout_, *fd);
-        return WrapField<T>(layout_, std::forward<decltype(fd)>(fd),
-                            evaluator_factory_);
+        return WrapField<T>(layout_, std::forward<decltype(fd)>(fd));
       }
     },
     std::move(data));
diff --git a/src/pcms/field/function_space/lagrange.h b/src/pcms/field/function_space/lagrange.h
index 54b2f695..d729d628 100644
--- a/src/pcms/field/function_space/lagrange.h
+++ b/src/pcms/field/function_space/lagrange.h
@@ -27,7 +27,19 @@ namespace pcms
 
 class LagrangeFunctionSpace : public FunctionSpace
 {
+  // Passkey: only members (the From* factories) can create a Key, so only they
+  // can construct a space, while std::make_shared can still call the ctor.
+  struct Key
+  {
+    explicit Key() = default;
+  };
+
 public:
+  LagrangeFunctionSpace(
+    Key, std::shared_ptr<const FieldLayout> layout,
+    std::function<FieldDataVariant(Type, FieldMetadata)> create_field_data_fn,
+    std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory) noexcept;
+
   enum class Backend
   {
     MeshFields,
@@ -41,25 +53,28 @@ class LagrangeFunctionSpace : public FunctionSpace
 #endif
 
   // Unstructured mesh — dispatches to MeshFields or native Omega_h backend
-  [[nodiscard]] static LagrangeFunctionSpace FromMesh(
+  [[nodiscard]] static std::shared_ptr<LagrangeFunctionSpace> FromMesh(
     Omega_h::Mesh& mesh, int order, int num_components,
     CoordinateSystem coordinate_system, std::string global_id_name = "global",
-    Backend backend = DefaultBackend);
+    Backend backend = DefaultBackend, std::string layout_name = "");
 
-  [[nodiscard]] static LagrangeFunctionSpace FromMesh(
+  [[nodiscard]] static std::shared_ptr<LagrangeFunctionSpace> FromMesh(
     Omega_h::Mesh& mesh, int order, int num_components,
     CoordinateSystem coordinate_system, Omega_h::Read<Omega_h::I8> owned_mask,
-    std::string global_id_name = "global", Backend backend = DefaultBackend);
+    std::string global_id_name = "global", Backend backend = DefaultBackend,
+    std::string layout_name = "");
 
   // Structured uniform grid — order-1 H1-conforming nodal field on a regular
   // grid
-  [[nodiscard]] static LagrangeFunctionSpace FromUniformGrid(
+  [[nodiscard]] static std::shared_ptr<LagrangeFunctionSpace> FromUniformGrid(
     const UniformGrid<2>& grid, int num_components,
-    CoordinateSystem coordinate_system, int order = 1);
+    CoordinateSystem coordinate_system, int order = 1,
+    std::string layout_name = "");
 
-  [[nodiscard]] static LagrangeFunctionSpace FromUniformGrid(
+  [[nodiscard]] static std::shared_ptr<LagrangeFunctionSpace> FromUniformGrid(
     const UniformGrid<3>& grid, int num_components,
-    CoordinateSystem coordinate_system, int order = 1);
+    CoordinateSystem coordinate_system, int order = 1,
+    std::string layout_name = "");
 
   [[nodiscard]] std::shared_ptr<const FieldLayout> GetLayout()
     const noexcept override;
@@ -77,11 +92,6 @@ class LagrangeFunctionSpace : public FunctionSpace
     Type value_type, const EvaluationRequest& request) const override;
 
 private:
-  explicit LagrangeFunctionSpace(
-    std::shared_ptr<const FieldLayout> layout,
-    std::function<FieldDataVariant(Type, FieldMetadata)> create_field_data_fn,
-    std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory) noexcept;
-
   std::shared_ptr<const FieldLayout> layout_;
   std::function<FieldDataVariant(Type, FieldMetadata)> create_field_data_fn_;
   std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory_;
diff --git a/src/pcms/field/function_space/polynomial_reconstruction.cpp b/src/pcms/field/function_space/polynomial_reconstruction.cpp
index 060b9c24..397e5f7c 100644
--- a/src/pcms/field/function_space/polynomial_reconstruction.cpp
+++ b/src/pcms/field/function_space/polynomial_reconstruction.cpp
@@ -17,13 +17,13 @@ namespace pcms
 {
 
 PolynomialReconstructionFunctionSpace::PolynomialReconstructionFunctionSpace(
-  std::shared_ptr<const FieldLayout> layout,
+  Key, std::shared_ptr<const FieldLayout> layout,
   std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory) noexcept
   : layout_(std::move(layout)), evaluator_factory_(std::move(evaluator_factory))
 {
 }
 
-PolynomialReconstructionFunctionSpace
+std::shared_ptr<PolynomialReconstructionFunctionSpace>
 PolynomialReconstructionFunctionSpace::Create(
   Rank2View<Real, HostMemorySpace> coords, CoordinateSystem coordinate_system,
   MLSOptions options)
@@ -40,11 +40,11 @@ PolynomialReconstructionFunctionSpace::Create(
     std::make_shared<PointCloudLocalizationFactory>(pc_layout, options);
   auto eval_factory = std::make_shared<PointCloudEvaluatorFactory>(
     pc_layout, localization, options);
-  return PolynomialReconstructionFunctionSpace(pc_layout,
-                                               std::move(eval_factory));
+  return std::make_shared<PolynomialReconstructionFunctionSpace>(
+    Key{}, pc_layout, std::move(eval_factory));
 }
 
-PolynomialReconstructionFunctionSpace
+std::shared_ptr<PolynomialReconstructionFunctionSpace>
 PolynomialReconstructionFunctionSpace::FromMesh(
   Omega_h::Mesh& mesh, int source_entity_dim,
   CoordinateSystem coordinate_system, MLSOptions options)
@@ -66,8 +66,8 @@ PolynomialReconstructionFunctionSpace::FromMesh(
     mesh, source_entity_dim, options);
   auto eval_factory = std::make_shared<PointCloudEvaluatorFactory>(
     mesh_layout, localization, options);
-  return PolynomialReconstructionFunctionSpace(mesh_layout,
-                                               std::move(eval_factory));
+  return std::make_shared<PolynomialReconstructionFunctionSpace>(
+    Key{}, mesh_layout, std::move(eval_factory));
 }
 
 std::shared_ptr<const FieldLayout>
@@ -93,8 +93,7 @@ FieldVariant PolynomialReconstructionFunctionSpace::CreateFieldImpl(
         "supported");
     } else {
       return WrapField<double>(
-        layout_, std::make_unique<SimpleFieldData<double>>(layout_, metadata),
-        evaluator_factory_);
+        layout_, std::make_unique<SimpleFieldData<double>>(layout_, metadata));
     }
   });
 }
@@ -120,7 +119,7 @@ FieldVariant PolynomialReconstructionFunctionSpace::CreateFieldImpl(
       "PolynomialReconstructionFunctionSpace::CreateField: field data size "
       "does not match layout");
   }
-  return WrapField<double>(layout_, std::move(fd), evaluator_factory_);
+  return WrapField<double>(layout_, std::move(fd));
 }
 
 PointEvaluatorVariant
diff --git a/src/pcms/field/function_space/polynomial_reconstruction.hpp b/src/pcms/field/function_space/polynomial_reconstruction.hpp
index 0222c282..890f61a7 100644
--- a/src/pcms/field/function_space/polynomial_reconstruction.hpp
+++ b/src/pcms/field/function_space/polynomial_reconstruction.hpp
@@ -23,14 +23,25 @@ class FieldEvaluatorFactory;
 
 class PolynomialReconstructionFunctionSpace : public FunctionSpace
 {
+  // Passkey: only members (the Create/FromMesh factories) can create a Key, so
+  // only they can construct a space, while std::make_shared can call the ctor.
+  struct Key
+  {
+    explicit Key() = default;
+  };
+
 public:
-  [[nodiscard]] static PolynomialReconstructionFunctionSpace Create(
-    Rank2View<Real, HostMemorySpace> coords, CoordinateSystem coordinate_system,
-    MLSOptions options = {});
+  PolynomialReconstructionFunctionSpace(
+    Key, std::shared_ptr<const FieldLayout> layout,
+    std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory) noexcept;
+
+  [[nodiscard]] static std::shared_ptr<PolynomialReconstructionFunctionSpace>
+  Create(Rank2View<Real, HostMemorySpace> coords,
+         CoordinateSystem coordinate_system, MLSOptions options = {});
 
-  [[nodiscard]] static PolynomialReconstructionFunctionSpace FromMesh(
-    Omega_h::Mesh& mesh, int source_entity_dim,
-    CoordinateSystem coordinate_system, MLSOptions options = {});
+  [[nodiscard]] static std::shared_ptr<PolynomialReconstructionFunctionSpace>
+  FromMesh(Omega_h::Mesh& mesh, int source_entity_dim,
+           CoordinateSystem coordinate_system, MLSOptions options = {});
 
   [[nodiscard]] std::shared_ptr<const FieldLayout> GetLayout()
     const noexcept override;
@@ -48,10 +59,6 @@ class PolynomialReconstructionFunctionSpace : public FunctionSpace
     Type value_type, const EvaluationRequest& request) const override;
 
 private:
-  explicit PolynomialReconstructionFunctionSpace(
-    std::shared_ptr<const FieldLayout> layout,
-    std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory) noexcept;
-
   std::shared_ptr<const FieldLayout> layout_;
   std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory_;
 };
diff --git a/src/pcms/field/function_space/spline.cpp b/src/pcms/field/function_space/spline.cpp
index f4530fc8..c746a8d2 100644
--- a/src/pcms/field/function_space/spline.cpp
+++ b/src/pcms/field/function_space/spline.cpp
@@ -11,20 +11,21 @@ namespace pcms
 {
 
 SplineFunctionSpace::SplineFunctionSpace(
-  std::shared_ptr<const FieldLayout> layout,
+  Key, std::shared_ptr<const FieldLayout> layout,
   std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory) noexcept
   : layout_(std::move(layout)), evaluator_factory_(std::move(evaluator_factory))
 {
 }
 
-SplineFunctionSpace SplineFunctionSpace::FromUniformGrid(
+std::shared_ptr<SplineFunctionSpace> SplineFunctionSpace::FromUniformGrid(
   const UniformGrid<2>& grid, CoordinateSystem coordinate_system)
 {
   auto layout =
     std::make_shared<UniformGridFieldLayout<2>>(grid, 1, coordinate_system, 1);
   auto evaluator_factory =
     std::make_shared<UniformGridSplineEvaluatorFactory2D>(layout);
-  return SplineFunctionSpace(layout, std::move(evaluator_factory));
+  return std::make_shared<SplineFunctionSpace>(Key{}, layout,
+                                               std::move(evaluator_factory));
 }
 
 std::shared_ptr<const FieldLayout> SplineFunctionSpace::GetLayout()
@@ -47,8 +48,7 @@ FieldVariant SplineFunctionSpace::CreateFieldImpl(Type value_type,
       throw pcms_error("SplineFunctionSpace: only double (Real) is supported");
     } else {
       return WrapField<double>(
-        layout_, std::make_unique<SimpleFieldData<double>>(layout_, metadata),
-        evaluator_factory_);
+        layout_, std::make_unique<SimpleFieldData<double>>(layout_, metadata));
     }
   });
 }
@@ -70,7 +70,7 @@ FieldVariant SplineFunctionSpace::CreateFieldImpl(FieldDataVariant data) const
       "SplineFunctionSpace::CreateField: field data size does not match "
       "layout");
   }
-  return WrapField<double>(layout_, std::move(fd), evaluator_factory_);
+  return WrapField<double>(layout_, std::move(fd));
 }
 
 PointEvaluatorVariant SplineFunctionSpace::CreatePointEvaluatorImpl(
diff --git a/src/pcms/field/function_space/spline.h b/src/pcms/field/function_space/spline.h
index 15c7d712..f8e467d7 100644
--- a/src/pcms/field/function_space/spline.h
+++ b/src/pcms/field/function_space/spline.h
@@ -21,8 +21,19 @@ namespace pcms
 
 class SplineFunctionSpace : public FunctionSpace
 {
+  // Passkey: only members (the From* factories) can create a Key, so only they
+  // can construct a space, while std::make_shared can still call the ctor.
+  struct Key
+  {
+    explicit Key() = default;
+  };
+
 public:
-  [[nodiscard]] static SplineFunctionSpace FromUniformGrid(
+  SplineFunctionSpace(
+    Key, std::shared_ptr<const FieldLayout> layout,
+    std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory) noexcept;
+
+  [[nodiscard]] static std::shared_ptr<SplineFunctionSpace> FromUniformGrid(
     const UniformGrid<2>& grid, CoordinateSystem coordinate_system);
 
   [[nodiscard]] std::shared_ptr<const FieldLayout> GetLayout()
@@ -41,10 +52,6 @@ class SplineFunctionSpace : public FunctionSpace
     Type value_type, const EvaluationRequest& request) const override;
 
 private:
-  explicit SplineFunctionSpace(
-    std::shared_ptr<const FieldLayout> layout,
-    std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory) noexcept;
-
   std::shared_ptr<const FieldLayout> layout_;
   std::shared_ptr<FieldEvaluatorFactory<Real>> evaluator_factory_;
 };
diff --git a/src/pcms/field/function_space/xgc.h b/src/pcms/field/function_space/xgc.h
index 164abc7a..cbea70b4 100644
--- a/src/pcms/field/function_space/xgc.h
+++ b/src/pcms/field/function_space/xgc.h
@@ -1,9 +1,9 @@
-#ifndef PCMS_XGC_FUNCTION_SPACE_H
-#define PCMS_XGC_FUNCTION_SPACE_H
+#ifndef PCMS_XGC_FIELD_FACTORY_H
+#define PCMS_XGC_FIELD_FACTORY_H
 
 #include "pcms/field/data/xgc.h"
 #include "pcms/field/field.h"
-#include "pcms/field/function_space.h"
+#include "pcms/field/field_factory.h"
 #include "pcms/field/field_metadata.h"
 #include "pcms/utility/common.h"
 #include <functional>
@@ -12,15 +12,22 @@
 namespace pcms
 {
 
-class XGCFunctionSpace : public FunctionSpace
+// XGC is a comm-only backend
+class XGCFieldFactory : public FieldFactory
 {
 public:
-  XGCFunctionSpace(const ReverseClassificationVertex& reverse_classification,
-                   std::function<int8_t(int, int)> in_overlap,
-                   LO num_plane_nodes)
+  XGCFieldFactory(const ReverseClassificationVertex& reverse_classification,
+                  std::function<int8_t(int, int)> in_overlap,
+                  LO num_plane_nodes, std::string layout_name = "")
     : layout_(std::make_shared<XGCFieldLayout>(
         reverse_classification, std::move(in_overlap), num_plane_nodes))
   {
+    layout_->SetName(std::move(layout_name));
+  }
+
+  void SetLayoutName(std::string name) const
+  {
+    layout_->SetName(std::move(name));
   }
 
   [[nodiscard]] std::shared_ptr<const FieldLayout> GetLayout()
@@ -29,9 +36,10 @@ class XGCFunctionSpace : public FunctionSpace
     return layout_;
   }
 
-  [[nodiscard]] CoordinateSystem GetCoordinateSystem() const noexcept override
+  [[nodiscard]] std::shared_ptr<const XGCFieldLayout> GetXGCLayout()
+    const noexcept
   {
-    return CoordinateSystem::XGC;
+    return layout_;
   }
 
 protected:
@@ -55,12 +63,12 @@ class XGCFunctionSpace : public FunctionSpace
         PCMS_ALWAYS_ASSERT(fd != nullptr);
         if (dynamic_cast<const XGCFieldData<T>*>(fd.get()) == nullptr) {
           throw pcms_error(
-            "XGCFunctionSpace::CreateField: requires XGCFieldData");
+            "XGCFieldFactory::CreateField: requires XGCFieldData");
         }
         if (fd->GetDOFHolderDataHost().size() !=
             static_cast<size_t>(layout_->GetFullDataSize())) {
           throw pcms_error(
-            "XGCFunctionSpace::CreateField: field data size does not match "
+            "XGCFieldFactory::CreateField: field data size does not match "
             "layout");
         }
         return WrapField<T>(layout_, std::move(fd));
@@ -68,23 +76,10 @@ class XGCFunctionSpace : public FunctionSpace
       std::move(data));
   }
 
-  [[nodiscard]] PointEvaluatorVariant CreatePointEvaluatorImpl(
-    Type /*value_type*/, const EvaluationRequest& /*request*/) const override
-  {
-    throw pcms_error("XGCFunctionSpace does not support point evaluation");
-  }
-
-public:
-  [[nodiscard]] std::shared_ptr<const XGCFieldLayout> GetXGCLayout()
-    const noexcept
-  {
-    return layout_;
-  }
-
 private:
-  std::shared_ptr<const XGCFieldLayout> layout_;
+  std::shared_ptr<XGCFieldLayout> layout_;
 };
 
 } // namespace pcms
 
-#endif // PCMS_XGC_FUNCTION_SPACE_H
+#endif // PCMS_XGC_FIELD_FACTORY_H
diff --git a/src/pcms/field/uniform_grid_binary_field.h b/src/pcms/field/uniform_grid_binary_field.h
index 6a71de7a..bce42439 100644
--- a/src/pcms/field/uniform_grid_binary_field.h
+++ b/src/pcms/field/uniform_grid_binary_field.h
@@ -39,9 +39,9 @@ CreateUniformGridBinaryField(Omega_h::Mesh& mesh, const UniformGrid<Dim>& grid)
   auto function_space = LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, CoordinateSystem::Cartesian);
   auto layout = std::dynamic_pointer_cast<const UniformGridFieldLayout<Dim>>(
-    function_space.GetLayout());
+    function_space->GetLayout());
   PCMS_ALWAYS_ASSERT(layout != nullptr);
-  auto field = function_space.template CreateField<Real>(FieldMetadata{});
+  Field<Real> field = function_space->template CreateFunction<Real>("");
 
   auto coord_view = layout->GetDOFHolderCoordinates();
   auto coords = coord_view.GetValues();
diff --git a/src/pcms/pythonapi/bind_field_base.cpp b/src/pcms/pythonapi/bind_field_base.cpp
index 77be9871..eca1c98a 100644
--- a/src/pcms/pythonapi/bind_field_base.cpp
+++ b/src/pcms/pythonapi/bind_field_base.cpp
@@ -286,10 +286,12 @@ void bind_create_field_module(py::module& m)
       },
       "DOF holder coordinates as a 2D numpy array (num_dof_holders × dim)");
 
-  py::class_<FunctionSpace>(m, "FunctionSpace")
+  py::class_<FunctionSpace, std::shared_ptr<FunctionSpace>>(m, "FunctionSpace")
     .def(
       "create_field",
-      [](const FunctionSpace& self) { return self.CreateField<Real>(); },
+      [](const FunctionSpace& self) -> Field<Real> {
+        return self.CreateFunction<Real>();
+      },
       "Create a Field<Real> for this function space.")
     .def(
       "create_point_evaluator",
@@ -317,8 +319,9 @@ void bind_create_field_module(py::module& m)
       "only).");
 
   // Bind LagrangeFunctionSpace as a concrete FunctionSpace subtype.
-  py::class_<LagrangeFunctionSpace, FunctionSpace> lagrange_space(
-    m, "LagrangeFunctionSpace");
+  py::class_<LagrangeFunctionSpace, FunctionSpace,
+             std::shared_ptr<LagrangeFunctionSpace>>
+    lagrange_space(m, "LagrangeFunctionSpace");
 
   // Backend selecting the underlying field layout. The conservative-projection
   // transfer operators require the native Omega_h backend.
@@ -386,7 +389,8 @@ void bind_create_field_module(py::module& m)
 
   // Bind PolynomialReconstructionFunctionSpace as a concrete FunctionSpace
   // subtype.
-  py::class_<PolynomialReconstructionFunctionSpace, FunctionSpace>(
+  py::class_<PolynomialReconstructionFunctionSpace, FunctionSpace,
+             std::shared_ptr<PolynomialReconstructionFunctionSpace>>(
     m, "PolynomialReconstructionFunctionSpace")
     .def_static(
       "from_coords",
diff --git a/src/pcms/transfer/omega_h_control_variate_projection.cpp b/src/pcms/transfer/omega_h_control_variate_projection.cpp
index bce1a89a..aea49cb2 100644
--- a/src/pcms/transfer/omega_h_control_variate_projection.cpp
+++ b/src/pcms/transfer/omega_h_control_variate_projection.cpp
@@ -18,7 +18,7 @@ OmegaHControlVariateProjection::OmegaHControlVariateProjection(
       target_layout_, target_space.GetCoordinateSystem(), samples_per_element,
       sampling, seed)),
     interpolator_(source_space, target_space),
-    control_variate_(target_space.CreateField<Real>())
+    control_variate_(target_space.CreateFunction<Real>())
 {
   const auto sample_coords =
     rhs_integrator_->GetIntegrationPoints().GetCoordinates();
diff --git a/test/field_test_utils.h b/test/field_test_utils.h
index 3fbb624f..5aaac78e 100644
--- a/test/field_test_utils.h
+++ b/test/field_test_utils.h
@@ -282,8 +282,8 @@ void CheckEvaluation(const Factory& factory, const Field<Real>& field,
                      double abs_tol = 1e-10)
 {
   auto device_coords =
-    CreateDeviceCoordinateView(pts, factory.GetCoordinateSystem());
-  auto evaluator = factory.template CreatePointEvaluator<Real>(
+    CreateDeviceCoordinateView(pts, factory->GetCoordinateSystem());
+  auto evaluator = factory->template CreatePointEvaluator<Real>(
     EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   CheckEvaluation<ExecutionSpace>(*evaluator, field, pts, func, abs_tol);
 }
@@ -314,9 +314,9 @@ void CheckFillMode(const Factory& factory, const Field<Real>& field,
                    Real fill_value, const std::vector<Real>& outside_pts)
 {
   auto device_coords =
-    CreateDeviceCoordinateView(outside_pts, factory.GetCoordinateSystem());
+    CreateDeviceCoordinateView(outside_pts, factory->GetCoordinateSystem());
   OutOfBoundsPolicy policy{OutOfBoundsMode::FILL, fill_value};
-  auto evaluator = factory.template CreatePointEvaluator<Real>(
+  auto evaluator = factory->template CreatePointEvaluator<Real>(
     EvaluationRequest::FromCoordinates(device_coords.coordinate_view, policy));
   CheckFillMode(*evaluator, field, fill_value, outside_pts);
 }
@@ -334,9 +334,9 @@ void CheckEvaluationWithFill(const Factory& factory, const Field<Real>& field,
   REQUIRE(static_cast<int>(is_inside.size()) == n);
 
   auto device_coords =
-    CreateDeviceCoordinateView(pts, factory.GetCoordinateSystem());
+    CreateDeviceCoordinateView(pts, factory->GetCoordinateSystem());
   OutOfBoundsPolicy policy{OutOfBoundsMode::FILL, fill_value};
-  auto evaluator = factory.template CreatePointEvaluator<Real>(
+  auto evaluator = factory->template CreatePointEvaluator<Real>(
     EvaluationRequest::FromCoordinates(device_coords.coordinate_view, policy));
 
   Kokkos::View<Real*, DeviceMemorySpace> out_device("out_device", n);
diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp
index 40150574..d546fbbe 100644
--- a/test/test_eqdsk.cpp
+++ b/test/test_eqdsk.cpp
@@ -182,10 +182,10 @@ TEST_CASE("EQDSKData with SplineFunctionSpace")
   {
     // auto spline_space = pcms::SplineFunctionSpace::FromEQDSK(eqdsk_data);
 
-    // auto psi_field = spline_space.CreateField<pcms::Real>();
+    // auto psi_field = spline_space->CreateFunction<pcms::Real>();
 
     auto [spline_space, psi_field] = pcms::MakeEQDSKField(eqdsk_data);
-    auto layout = spline_space.GetLayout();
+    auto layout = spline_space->GetLayout();
     const int layout_size = layout->OwnedSize();
     const int psizr_size = static_cast<int>(eqdsk_data.PSIZR.extent(0));
 
@@ -230,7 +230,7 @@ TEST_CASE("EQDSKData with SplineFunctionSpace")
     auto eval_request = pcms::EvaluationRequest::FromCoordinates(coord_view);
 
     auto evaluator =
-      spline_space.CreatePointEvaluator<pcms::Real>(eval_request);
+      spline_space->CreatePointEvaluator<pcms::Real>(eval_request);
 
     auto eval_results_1d = Kokkos::View<pcms::Real*, pcms::DeviceMemorySpace>(
       "eval_results", num_eval_points);
diff --git a/test/test_field_communication.cpp b/test/test_field_communication.cpp
index 83b8bc5b..dcee3252 100644
--- a/test/test_field_communication.cpp
+++ b/test/test_field_communication.cpp
@@ -78,15 +78,12 @@ static void test_shared_layout(Omega_h::Library& lib,
                       redev::Partition{partition});
     auto* app = cpl.AddApplication("shared_layout");
     auto factory = pcms::LagrangeFunctionSpace::FromMesh(
-      mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-    auto layout = factory.GetLayout();
-    app->AddLayout("shared", layout);
+      mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
+      pcms::LagrangeFunctionSpace::DefaultBackend, "shared");
+    auto f1 = app->AddField(factory->CreateFunction<Real>("field_a"));
     PCMS_ALWAYS_ASSERT(app->GetLayoutCommunicatorCount() == 1);
-    auto f1 = app->AddField("field_a", factory.CreateField<Real>());
-    PCMS_ALWAYS_ASSERT(app->GetLayoutCommunicatorCount() ==
-                       1); // still 1 after adding field
-    auto f2 = app->AddField("field_b", factory.CreateField<Real>());
-    PCMS_ALWAYS_ASSERT(app->GetLayoutCommunicatorCount() == 1); // still 1
+    auto f2 = app->AddField(factory->CreateFunction<Real>("field_b"));
+    PCMS_ALWAYS_ASSERT(app->GetLayoutCommunicatorCount() == 1); // shared layout
     app->ReceivePhase([&]() {
       f1.Receive();
       f2.Receive();
@@ -100,11 +97,10 @@ static void test_shared_layout(Omega_h::Library& lib,
                       redev::Partition{});
     auto* app = cpl.AddApplication("shared_layout");
     auto factory = pcms::LagrangeFunctionSpace::FromMesh(
-      mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-    auto layout = factory.GetLayout();
-    app->AddLayout("shared", layout);
-    auto f1 = app->AddField("field_a", factory.CreateField<Real>());
-    auto f2 = app->AddField("field_b", factory.CreateField<Real>());
+      mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
+      pcms::LagrangeFunctionSpace::DefaultBackend, "shared");
+    auto f1 = app->AddField(factory->CreateFunction<Real>("field_a"));
+    auto f2 = app->AddField(factory->CreateFunction<Real>("field_b"));
     app->SendPhase([&]() {
       f1.Send();
       f2.Send();
@@ -125,7 +121,7 @@ void client1(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name,
 
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, order, 1, pcms::CoordinateSystem::Cartesian);
-  auto layout = factory.GetLayout();
+  auto layout = factory->GetLayout();
   auto gids = layout->GetGidsHost();
   const auto n = layout->GetNumOwnedDofHolder();
   Omega_h::HostWrite<Real> ids(n);
@@ -134,7 +130,7 @@ void client1(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name,
     "id gid", Kokkos::RangePolicy<pcms::HostMemorySpace::execution_space>(0, n),
     [=](int i) { ids[i] = gids[i]; });
 
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
   field.SetDOFHolderDataHost(
     pcms::Rank2View<const Real, pcms::HostMemorySpace>(ids.data(), n, 1));
 
@@ -157,11 +153,11 @@ void client2(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name,
 
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, order, 1, pcms::CoordinateSystem::Cartesian);
-  auto layout = factory.GetLayout();
+  auto layout = factory->GetLayout();
   auto gids = layout->GetGidsHost();
   const auto n = layout->GetNumOwnedDofHolder();
 
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
   pcms::FieldLayoutCommunicator layout_comm(comm_name + "2", comm, rdv, channel,
                                             *layout);
   pcms::FieldCommunicator<pcms::Real> field_comm(layout_comm.GetName(),
@@ -217,14 +213,14 @@ void server(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name,
 
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, order, 1, pcms::CoordinateSystem::Cartesian);
-  auto layout = factory.GetLayout();
+  auto layout = factory->GetLayout();
   const auto n = layout->GetNumOwnedDofHolder();
   Omega_h::HostWrite<Real> ids(n);
   Kokkos::parallel_for(
     "id 0", Kokkos::RangePolicy<pcms::HostMemorySpace::execution_space>(0, n),
     [=](int i) { ids[i] = 0; });
 
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
   pcms::FieldLayoutCommunicator layout_comm1(comm_name + "1", comm, rdv,
                                              channel1, *layout);
   pcms::FieldLayoutCommunicator layout_comm2(comm_name + "2", comm, rdv,
diff --git a/test/test_field_copy.cpp b/test/test_field_copy.cpp
index 9906bbbc..46a1b94e 100644
--- a/test/test_field_copy.cpp
+++ b/test/test_field_copy.cpp
@@ -22,20 +22,20 @@ void test_copy(Omega_h::CommPtr world, int dim, int order, int num_components)
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, order, num_components, pcms::CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto layout = factory.GetLayout();
+  auto layout = factory->GetLayout();
   int ndata = layout->GetNumOwnedDofHolder() * num_components;
   Omega_h::HostWrite<Real> ids(ndata);
   Kokkos::parallel_for(
     Kokkos::RangePolicy<pcms::HostMemorySpace::execution_space>(0, ndata),
     [=](int i) { ids[i] = i; });
 
-  auto original = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto original = factory->CreateFunction<Real>();
   original.SetDOFHolderDataHost(
     pcms::Rank2View<const Real, pcms::HostMemorySpace>(
       ids.data(), layout->GetNumOwnedDofHolder(), num_components));
 
-  auto copied = factory.CreateField<Real>(pcms::FieldMetadata{});
-  pcms::Copy<Real> copy(factory, factory);
+  auto copied = factory->CreateFunction<Real>();
+  pcms::Copy<Real> copy(*factory, *factory);
   copy.Apply(original, copied);
   auto copied_array = pcms::FlattenToRank1View(copied.GetDOFHolderDataHost());
 
diff --git a/test/test_field_evaluation.cpp b/test/test_field_evaluation.cpp
index c4a29378..ff238e11 100644
--- a/test/test_field_evaluation.cpp
+++ b/test/test_field_evaluation.cpp
@@ -30,10 +30,10 @@ TEST_CASE("evaluate linear 2d omega_h_field")
                                  100, 0, false);
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
   pcms::test::SetField(
-    field.GetData(), *factory.GetLayout(),
+    field.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
   pcms::test::CheckEvaluation(
     factory, field, pcms::test::StandardEvalCoords2D(),
@@ -75,7 +75,7 @@ TEST_CASE("evaluate quadratic 2d meshfields_field")
     });
 
   Omega_h::HostWrite<Real> test_f_host(test_f);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
   field.GetData().SetDOFHolderDataHost(
     pcms::Rank2View<const Real, pcms::HostMemorySpace>(test_f_host.data(),
                                                        test_f_host.size(), 1));
diff --git a/test/test_field_interpolation.cpp b/test/test_field_interpolation.cpp
index fdada727..bc88ab97 100644
--- a/test/test_field_interpolation.cpp
+++ b/test/test_field_interpolation.cpp
@@ -27,12 +27,12 @@ TEST_CASE("interpolate linear 2d omega_h_field")
     Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 100, 100, 0, false);
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
-  auto interpolated = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
+  auto interpolated = factory->CreateFunction<Real>();
   pcms::test::SetField(
     field, OMEGA_H_LAMBDA(Real x, Real y) { return -0.3 * x + 0.5 * y; });
 
-  pcms::Interpolator<Real> interp(factory, factory);
+  pcms::Interpolator<Real> interp(*factory, *factory);
   interp.Apply(field, interpolated);
   auto interpolated_dof =
     pcms::FlattenToRank1View(interpolated.GetDOFHolderDataHost());
@@ -55,7 +55,7 @@ TEST_CASE("interpolate quadratic 2d meshfields_field")
   auto factory2 = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::MeshFields);
-  auto layout = factory2.GetLayout();
+  auto layout = factory2->GetLayout();
   const auto nverts = mesh.nents(0);
   const auto nedges = mesh.nents(1);
   auto mesh_coords = mesh.coords();
@@ -80,12 +80,12 @@ TEST_CASE("interpolate quadratic 2d meshfields_field")
     });
 
   Omega_h::HostWrite<Real> test_f_host(test_f);
-  auto field = factory2.CreateField<Real>(pcms::FieldMetadata{});
-  auto interpolated = factory2.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory2->CreateFunction<Real>();
+  auto interpolated = factory2->CreateFunction<Real>();
   field.SetDOFHolderDataHost(pcms::Rank2View<const Real, pcms::HostMemorySpace>(
     test_f_host.data(), test_f_host.size(), 1));
 
-  pcms::Interpolator<Real> interp(factory2, factory2);
+  pcms::Interpolator<Real> interp(*factory2, *factory2);
   interp.Apply(field, interpolated);
 
   auto interpolated_dof =
@@ -126,13 +126,13 @@ TEST_CASE("Interpolator: construct once, apply twice with different data")
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
 
-  auto source = factory.CreateField<Real>(pcms::FieldMetadata{});
-  auto target = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto source = factory->CreateFunction<Real>();
+  auto target = factory->CreateFunction<Real>();
 
   // Construct Interpolator once — localization happens here.
   // Because source and target share the same layout (same factory),
   // the target DOF holder coordinates are the source mesh node coordinates.
-  pcms::Interpolator<pcms::Real> interp(factory, factory);
+  pcms::Interpolator<pcms::Real> interp(*factory, *factory);
 
   // First application: linear function f(x,y) = -0.3*x + 0.5*y
   pcms::test::SetField(
diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp
index 54b16038..fbb34ddb 100644
--- a/test/test_mesh_intersection_field_transfer.cpp
+++ b/test/test_mesh_intersection_field_transfer.cpp
@@ -39,7 +39,7 @@ Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, const Omega_h::LOs& ev2v)
   return mesh;
 }
 
-pcms::LagrangeFunctionSpace MakeP1Space(
+std::shared_ptr<pcms::LagrangeFunctionSpace> MakeP1Space(
   Omega_h::Mesh& mesh, const std::string& global_id_name = "global")
 {
   return pcms::LagrangeFunctionSpace::FromMesh(
@@ -47,7 +47,7 @@ pcms::LagrangeFunctionSpace MakeP1Space(
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
 }
 
-pcms::LagrangeFunctionSpace MakeP0Space(Omega_h::Mesh& mesh)
+std::shared_ptr<pcms::LagrangeFunctionSpace> MakeP0Space(Omega_h::Mesh& mesh)
 {
   return pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 0, 1, pcms::CoordinateSystem::Cartesian, "global",
@@ -146,10 +146,10 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields",
   auto source_space = MakeP1Space(source_mesh);
   auto target_space = MakeP1Space(target_mesh);
 
-  auto source = source_space.CreateField<pcms::Real>();
-  auto target = target_space.CreateField<pcms::Real>();
+  auto source = source_space->CreateFunction<pcms::Real>();
+  auto target = target_space->CreateFunction<pcms::Real>();
 
-  pcms::OmegaHConservativeProjection projection(source_space, target_space);
+  pcms::OmegaHConservativeProjection projection(*source_space, *target_space);
 
   SECTION("constant field is preserved and conserved")
   {
@@ -208,10 +208,10 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral",
   auto source_space = MakeP1Space(source_mesh);
   auto target_space = MakeP1Space(target_mesh);
 
-  auto source = source_space.CreateField<pcms::Real>();
-  auto target = target_space.CreateField<pcms::Real>();
+  auto source = source_space->CreateFunction<pcms::Real>();
+  auto target = target_space->CreateFunction<pcms::Real>();
 
-  pcms::OmegaHConservativeProjection projection(source_space, target_space);
+  pcms::OmegaHConservativeProjection projection(*source_space, *target_space);
 
   pcms::test::SetField(
     source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) {
@@ -246,18 +246,18 @@ TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local "
   auto source_space = MakeP1Space(source_mesh);
   auto target_space = MakeP1Space(target_mesh, "reordered_global");
 
-  auto source = source_space.CreateField<pcms::Real>();
-  auto target = target_space.CreateField<pcms::Real>();
+  auto source = source_space->CreateFunction<pcms::Real>();
+  auto target = target_space->CreateFunction<pcms::Real>();
 
   pcms::test::SetField(
     source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; });
 
-  pcms::OmegaHConservativeProjection projection(source_space, target_space);
+  pcms::OmegaHConservativeProjection projection(*source_space, *target_space);
   projection.Apply(source, target);
 
   const auto target_values = target.GetDOFHolderDataHost();
   const auto target_coords =
-    target_space.GetLayout()->GetDOFHolderCoordinates().GetValues();
+    target_space->GetLayout()->GetDOFHolderCoordinates().GetValues();
   const auto target_coords_h = pcms::test::CopyCoordinatesToHost(
     target_coords, target_mesh.nverts(), target_mesh.dim());
   REQUIRE(static_cast<Omega_h::LO>(target_values.extent(0)) ==
@@ -285,18 +285,18 @@ TEST_CASE("OmegaHConservativeProjection maps sparse target GIDs to active "
   auto source_space = MakeP1Space(source_mesh);
   auto target_space = MakeP1Space(target_mesh, "sparse_global");
 
-  auto source = source_space.CreateField<pcms::Real>();
-  auto target = target_space.CreateField<pcms::Real>();
+  auto source = source_space->CreateFunction<pcms::Real>();
+  auto target = target_space->CreateFunction<pcms::Real>();
 
   pcms::test::SetField(
     source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; });
 
-  pcms::OmegaHConservativeProjection projection(source_space, target_space);
+  pcms::OmegaHConservativeProjection projection(*source_space, *target_space);
   projection.Apply(source, target);
 
   const auto target_values = target.GetDOFHolderDataHost();
   const auto target_coords =
-    target_space.GetLayout()->GetDOFHolderCoordinates().GetValues();
+    target_space->GetLayout()->GetDOFHolderCoordinates().GetValues();
   const auto target_coords_h = pcms::test::CopyCoordinatesToHost(
     target_coords, target_mesh.nverts(), target_mesh.dim());
   REQUIRE(static_cast<Omega_h::LO>(target_values.extent(0)) ==
@@ -325,10 +325,10 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target",
   auto source_space = MakeP0Space(source_mesh);
   auto target_space = MakeP1Space(target_mesh);
 
-  auto source = source_space.CreateField<pcms::Real>();
-  auto target = target_space.CreateField<pcms::Real>();
+  auto source = source_space->CreateFunction<pcms::Real>();
+  auto target = target_space->CreateFunction<pcms::Real>();
 
-  pcms::OmegaHConservativeProjection projection(source_space, target_space);
+  pcms::OmegaHConservativeProjection projection(*source_space, *target_space);
 
   SECTION("constant field is preserved and conserved")
   {
@@ -379,10 +379,10 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target",
   auto source_space = MakeP1Space(source_mesh);
   auto target_space = MakeP0Space(target_mesh);
 
-  auto source = source_space.CreateField<pcms::Real>();
-  auto target = target_space.CreateField<pcms::Real>();
+  auto source = source_space->CreateFunction<pcms::Real>();
+  auto target = target_space->CreateFunction<pcms::Real>();
 
-  pcms::OmegaHConservativeProjection projection(source_space, target_space);
+  pcms::OmegaHConservativeProjection projection(*source_space, *target_space);
 
   SECTION("constant field is preserved and conserved")
   {
diff --git a/test/test_omega_h_field2_outofbounds.cpp b/test/test_omega_h_field2_outofbounds.cpp
index da0cd338..8581dcb5 100644
--- a/test/test_omega_h_field2_outofbounds.cpp
+++ b/test/test_omega_h_field2_outofbounds.cpp
@@ -19,9 +19,9 @@ TEST_CASE("omega_h_field2 out of bounds FILL mode")
     Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 10, 10, 0, false);
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
   pcms::test::SetField(
-    field.GetData(), *factory.GetLayout(),
+    field.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + y; });
 
   Real fill_value = -999.0;
diff --git a/test/test_omega_h_intersection_rhs_integrator.cpp b/test/test_omega_h_intersection_rhs_integrator.cpp
index b778f2ae..fcbb3aed 100644
--- a/test/test_omega_h_intersection_rhs_integrator.cpp
+++ b/test/test_omega_h_intersection_rhs_integrator.cpp
@@ -95,14 +95,15 @@ std::unordered_map<pcms::GO, pcms::Real> ExpectedLoadByGid(
 // Assembles the integrator's load vector for source_field and compares it,
 // per global DOF id, against the supplied expected values.
 void CheckAssembledLoadMatches(
-  pcms::FunctionSpace& source_space, pcms::LinearFormIntegrator& integrator,
+  const std::shared_ptr<const pcms::FunctionSpace>& source_space,
+  pcms::LinearFormIntegrator& integrator,
   const pcms::Field<pcms::Real>& source_field,
   const std::unordered_map<pcms::GO, pcms::Real>& expected)
 {
   const auto& pts = integrator.GetIntegrationPoints();
   const std::size_t npts = pts.GetCoordinates().GetValues().extent(0);
 
-  auto evaluator = source_space.CreatePointEvaluator<pcms::Real>(
+  auto evaluator = source_space->CreatePointEvaluator<pcms::Real>(
     pcms::EvaluationRequest::FromCoordinates(pts.GetCoordinates()));
   Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace> sampled("sampled", npts,
                                                               1);
@@ -144,7 +145,7 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: integration points lie inside "
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
 
   auto integrator =
-    pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space);
+    pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space);
 
   const auto raw_coords =
     integrator->GetIntegrationPoints().GetCoordinates().GetValues();
@@ -181,15 +182,15 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: zero source field gives zero "
     target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
 
-  auto source_field = source_space.CreateField<pcms::Real>();
+  auto source_field = source_space->CreateFunction<pcms::Real>();
   // Default-constructed field has zero data.
 
   auto integrator =
-    pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space);
+    pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space);
   const auto& pts = integrator->GetIntegrationPoints();
   const std::size_t npts = pts.GetCoordinates().GetValues().extent(0);
 
-  auto evaluator = source_space.CreatePointEvaluator<pcms::Real>(
+  auto evaluator = source_space->CreatePointEvaluator<pcms::Real>(
     pcms::EvaluationRequest::FromCoordinates(pts.GetCoordinates()));
 
   Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace> sampled("sampled", npts,
@@ -224,18 +225,18 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: constant field matches "
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
 
   const double c = 2.0;
-  auto source_field = source_space.CreateField<pcms::Real>();
+  auto source_field = source_space->CreateFunction<pcms::Real>();
   pcms::test::SetField(
     source_field, OMEGA_H_LAMBDA(pcms::Real, pcms::Real) { return c; });
 
   const auto target_layout =
     std::dynamic_pointer_cast<const pcms::OmegaHLagrangeLayout>(
-      target_space.GetLayout());
+      target_space->GetLayout());
   const std::vector<pcms::Real> g(target_mesh.nverts(), c);
   const auto expected = ExpectedLoadByGid(target_mesh, *target_layout, g);
 
   auto integrator =
-    pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space);
+    pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space);
   CheckAssembledLoadMatches(source_space, *integrator, source_field, expected);
 }
 
@@ -260,13 +261,13 @@ TEST_CASE(
     target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
 
-  auto source_field = source_space.CreateField<pcms::Real>();
+  auto source_field = source_space->CreateFunction<pcms::Real>();
   pcms::test::SetField(
     source_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + y; });
 
   const auto target_layout =
     std::dynamic_pointer_cast<const pcms::OmegaHLagrangeLayout>(
-      target_space.GetLayout());
+      target_space->GetLayout());
   const auto tgt_coords_h =
     Omega_h::HostRead<Omega_h::Real>(target_mesh.coords());
   std::vector<pcms::Real> g(target_mesh.nverts());
@@ -276,7 +277,7 @@ TEST_CASE(
   const auto expected = ExpectedLoadByGid(target_mesh, *target_layout, g);
 
   auto integrator =
-    pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space);
+    pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space);
   CheckAssembledLoadMatches(source_space, *integrator, source_field, expected);
 }
 
@@ -296,7 +297,7 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts",
       target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
     REQUIRE_THROWS(
-      pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space));
+      pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space));
   }
 
   SECTION("multi-component target space throws")
@@ -308,7 +309,7 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts",
       target_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
     REQUIRE_THROWS(
-      pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space));
+      pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space));
   }
 
   SECTION("non-Cartesian source coordinate system throws")
@@ -320,7 +321,7 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts",
       target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
     REQUIRE_THROWS(
-      pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space));
+      pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space));
   }
 
   SECTION("non-Cartesian target coordinate system throws")
@@ -332,6 +333,6 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts",
       target_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
     REQUIRE_THROWS(
-      pcms::BuildOmegaHConservativeRHSIntegrator(source_space, target_space));
+      pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space));
   }
 }
diff --git a/test/test_omega_h_lagrange_field.cpp b/test/test_omega_h_lagrange_field.cpp
index 259a0165..c17d45b1 100644
--- a/test/test_omega_h_lagrange_field.cpp
+++ b/test/test_omega_h_lagrange_field.cpp
@@ -115,8 +115,8 @@ TEST_CASE("OmegaHLagrangeLayout layout sharing")
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
 
-  auto f1 = factory.CreateField<Real>(pcms::FieldMetadata{});
-  auto f2 = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto f1 = factory->CreateFunction<Real>();
+  auto f2 = factory->CreateFunction<Real>();
 
   REQUIRE(&f1.GetLayout() == &f2.GetLayout());
 }
@@ -129,9 +129,9 @@ TEST_CASE("OmegaHLagrangeField order-1: set/get DOF data round-trip")
   auto mesh = MakeBox2D(lib.world());
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
-  int n = factory.GetLayout()->GetNumOwnedDofHolder();
+  int n = factory->GetLayout()->GetNumOwnedDofHolder();
   std::vector<Real> data(n);
   for (int i = 0; i < n; ++i)
     data[i] = static_cast<Real>(i);
@@ -151,10 +151,10 @@ TEST_CASE("OmegaHLagrangeField order-1: linear function evaluation")
   auto mesh = MakeBox2D(lib.world(), 20, 20);
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
   pcms::test::SetField(
-    field.GetData(), *factory.GetLayout(),
+    field.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
   pcms::test::CheckEvaluation(
     factory, field, pcms::test::StandardEvalCoords2D(),
@@ -169,10 +169,10 @@ TEST_CASE("MeshFieldsAdapter order-1: linear function evaluation (shared util)")
   auto mesh = MakeBox2D(lib.world(), 20, 20);
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
   pcms::test::SetField(
-    field.GetData(), *factory.GetLayout(),
+    field.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
   pcms::test::CheckEvaluation(
     factory, field, pcms::test::StandardEvalCoords2D(),
@@ -185,10 +185,10 @@ TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode")
   auto mesh = MakeBox2D(lib.world());
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
   pcms::test::SetField(
-    field.GetData(), *factory.GetLayout(),
+    field.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
 
   Real fill_value = -999.0;
@@ -202,12 +202,12 @@ TEST_CASE("OmegaHLagrangeField order-1: serialize / deserialize round-trip")
   auto mesh = MakeBox2D(lib.world());
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
   pcms::test::SetField(
-    field.GetData(), *factory.GetLayout(),
+    field.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
-  pcms::test::CheckSerializeDeserialize(*factory.GetLayout(), field.GetData());
+  pcms::test::CheckSerializeDeserialize(*factory->GetLayout(), field.GetData());
 }
 
 TEST_CASE(
@@ -221,9 +221,9 @@ TEST_CASE(
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, nc, pcms::CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
-  const int n = factory.GetLayout()->GetNumOwnedDofHolder();
+  const int n = factory->GetLayout()->GetNumOwnedDofHolder();
   std::vector<Real> data(static_cast<size_t>(n) * nc);
   for (int i = 0; i < n; ++i)
     for (int c = 0; c < nc; ++c)
@@ -231,7 +231,7 @@ TEST_CASE(
   field.GetData().SetDOFHolderDataHost(
     pcms::Rank2View<const Real, pcms::HostMemorySpace>(data.data(), n, nc));
 
-  pcms::test::CheckSerializeDeserialize(*factory.GetLayout(), field.GetData());
+  pcms::test::CheckSerializeDeserialize(*factory->GetLayout(), field.GetData());
 }
 
 // ---- Order-0 field tests ----------------------------------------------------
@@ -242,9 +242,9 @@ TEST_CASE("OmegaHLagrangeField order-0: set/get DOF data round-trip")
   auto mesh = MakeBox2D(lib.world());
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 0, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
-  int n = factory.GetLayout()->GetNumOwnedDofHolder();
+  int n = factory->GetLayout()->GetNumOwnedDofHolder();
   std::vector<Real> data(n, 3.14);
   pcms::Rank2View<const Real, pcms::HostMemorySpace> view(data.data(), n, 1);
   field.GetData().SetDOFHolderDataHost(view);
@@ -261,7 +261,7 @@ TEST_CASE("OmegaHLagrangeField order-0: constant field evaluation")
   auto mesh = MakeBox2D(lib.world(), 10, 10);
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 0, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
   const Real kValue = 42.0;
   int nelems = mesh.nelems();
@@ -273,8 +273,8 @@ TEST_CASE("OmegaHLagrangeField order-0: constant field evaluation")
   auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
-    pcms::test::CreateDeviceCoordinateView(pts, factory.GetCoordinateSystem());
-  auto evaluator = factory.CreatePointEvaluator<Real>(
+    pcms::test::CreateDeviceCoordinateView(pts, factory->GetCoordinateSystem());
+  auto evaluator = factory->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
   Kokkos::View<Real*, pcms::DeviceMemorySpace> eval_device("eval", n);
@@ -296,7 +296,7 @@ TEST_CASE("OmegaHLagrangeField order-0: out-of-bounds FILL mode")
   auto mesh = MakeBox2D(lib.world());
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 0, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
   int nelems = mesh.nelems();
   std::vector<Real> data(nelems, 1.0);
@@ -315,7 +315,7 @@ TEST_CASE("OmegaHLagrangeField order-0: serialize / deserialize round-trip")
   auto mesh = MakeBox2D(lib.world());
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 0, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
   int nelems = mesh.nelems();
   std::vector<Real> data(nelems);
@@ -325,7 +325,7 @@ TEST_CASE("OmegaHLagrangeField order-0: serialize / deserialize round-trip")
                                                           1);
   field.GetData().SetDOFHolderDataHost(view);
 
-  pcms::test::CheckSerializeDeserialize(*factory.GetLayout(), field.GetData());
+  pcms::test::CheckSerializeDeserialize(*factory->GetLayout(), field.GetData());
 }
 
 // ---- Layout sharing communicator contract -----------------------------------
@@ -353,7 +353,7 @@ TEST_CASE("OmegaHLagrangeField: field valid after layout destruction")
   {
     auto factory = pcms::LagrangeFunctionSpace::FromMesh(
       mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-    field.emplace(factory.CreateField<Real>(pcms::FieldMetadata{}));
+    field.emplace(factory->CreateFunction<Real>());
   } // factory goes out of scope; field keeps layout alive
 
   pcms::test::SetField(
diff --git a/test/test_omega_h_mass_integrator.cpp b/test/test_omega_h_mass_integrator.cpp
index b916651c..f77fd85f 100644
--- a/test/test_omega_h_mass_integrator.cpp
+++ b/test/test_omega_h_mass_integrator.cpp
@@ -102,9 +102,9 @@ TEST_CASE(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
 
-  const auto ref = BuildReferenceMassMap(mesh, space);
+  const auto ref = BuildReferenceMassMap(mesh, *space);
 
-  auto integrator = pcms::BuildOmegaHMassIntegrator(space);
+  auto integrator = pcms::BuildOmegaHMassIntegrator(*space);
   Mat mat = integrator->GetMatrix();
 
   for (const auto& [key, ref_val] : ref) {
@@ -134,11 +134,11 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass",
 
   const auto layout =
     std::dynamic_pointer_cast<const pcms::OmegaHLagrangeLayout>(
-      space.GetLayout());
+      space->GetLayout());
   const auto& gids = layout->GetGidsHost();
   const int nverts = mesh.nverts();
 
-  auto integrator = pcms::BuildOmegaHMassIntegrator(space);
+  auto integrator = pcms::BuildOmegaHMassIntegrator(*space);
   Mat mat = integrator->GetMatrix();
 
   // Sum all row sums.
@@ -174,7 +174,7 @@ TEST_CASE("OmegaHMassIntegrator: rejects invalid layouts", "[mass_integrator]")
     auto space = pcms::LagrangeFunctionSpace::FromMesh(
       mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
-    REQUIRE_THROWS(pcms::BuildOmegaHMassIntegrator(space));
+    REQUIRE_THROWS(pcms::BuildOmegaHMassIntegrator(*space));
   }
 
   SECTION("non-Cartesian coordinate system throws")
@@ -182,6 +182,6 @@ TEST_CASE("OmegaHMassIntegrator: rejects invalid layouts", "[mass_integrator]")
     auto space = pcms::LagrangeFunctionSpace::FromMesh(
       mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global",
       pcms::LagrangeFunctionSpace::Backend::OmegaH);
-    REQUIRE_THROWS(pcms::BuildOmegaHMassIntegrator(space));
+    REQUIRE_THROWS(pcms::BuildOmegaHMassIntegrator(*space));
   }
 }
diff --git a/test/test_omega_h_mc_rhs_integrator.cpp b/test/test_omega_h_mc_rhs_integrator.cpp
index a26095ce..eabe5404 100644
--- a/test/test_omega_h_mc_rhs_integrator.cpp
+++ b/test/test_omega_h_mc_rhs_integrator.cpp
@@ -53,7 +53,7 @@ Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal)
   return mesh;
 }
 
-pcms::LagrangeFunctionSpace MakeP1Space(Omega_h::Mesh& mesh)
+std::shared_ptr<pcms::LagrangeFunctionSpace> MakeP1Space(Omega_h::Mesh& mesh)
 {
   return pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
@@ -61,13 +61,14 @@ pcms::LagrangeFunctionSpace MakeP1Space(Omega_h::Mesh& mesh)
 }
 
 // Evaluates source_field at the integrator's sample points and assembles.
-void EvaluateAndAssemble(pcms::LinearFormIntegrator& integrator,
-                         const pcms::FunctionSpace& source_space,
-                         const pcms::Field<pcms::Real>& source_field)
+void EvaluateAndAssemble(
+  pcms::LinearFormIntegrator& integrator,
+  const std::shared_ptr<const pcms::FunctionSpace>& source_space,
+  const pcms::Field<pcms::Real>& source_field)
 {
   const auto& pts = integrator.GetIntegrationPoints();
   const std::size_t npts = pts.GetCoordinates().GetValues().extent(0);
-  auto evaluator = source_space.CreatePointEvaluator<pcms::Real>(
+  auto evaluator = source_space->CreatePointEvaluator<pcms::Real>(
     pcms::EvaluationRequest::FromCoordinates(pts.GetCoordinates()));
   Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace> sampled("sampled", npts,
                                                               1);
@@ -91,7 +92,7 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: sample points lie inside the domain",
   const int samples_per_element = 16;
   for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) {
     pcms::OmegaHMonteCarloRHSIntegrator integrator(
-      target_space, samples_per_element, sampling);
+      *target_space, samples_per_element, sampling);
     const auto raw_coords =
       integrator.GetIntegrationPoints().GetCoordinates().GetValues();
     auto coords_view = Kokkos::View<const pcms::Real**, Kokkos::LayoutRight,
@@ -125,12 +126,12 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly",
   auto source_space = MakeP1Space(source_mesh);
   auto target_space = MakeP1Space(target_mesh);
 
-  auto source_field = source_space.CreateField<pcms::Real>();
+  auto source_field = source_space->CreateFunction<pcms::Real>();
   pcms::test::SetField(
     source_field, OMEGA_H_LAMBDA(pcms::Real, pcms::Real) { return 2.0; });
 
   for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) {
-    pcms::OmegaHMonteCarloRHSIntegrator integrator(target_space, 8, sampling);
+    pcms::OmegaHMonteCarloRHSIntegrator integrator(*target_space, 8, sampling);
     EvaluateAndAssemble(integrator, source_space, source_field);
 
     PetscScalar sum = 0.0;
@@ -157,14 +158,14 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target "
   auto source_space = MakeP1Space(source_mesh);
   auto target_space = MakeP1Space(target_mesh);
 
-  auto source_field = source_space.CreateField<pcms::Real>();
-  auto target_field = target_space.CreateField<pcms::Real>();
+  auto source_field = source_space->CreateFunction<pcms::Real>();
+  auto target_field = target_space->CreateFunction<pcms::Real>();
   pcms::test::SetField(
     source_field,
     OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return 3.0 * x - y + 0.25; });
 
   pcms::OmegaHControlVariateProjection projection(
-    source_space, target_space, /*samples_per_element=*/4,
+    *source_space, *target_space, /*samples_per_element=*/4,
     pcms::MonteCarloSampling::UniformRandom);
   projection.Apply(source_field, target_field);
 
@@ -193,16 +194,16 @@ TEST_CASE("OmegaHControlVariateProjection: reduces error vs plain Monte Carlo",
   auto source_space = MakeP1Space(source_mesh);
   auto target_space = MakeP1Space(target_mesh);
 
-  auto source_field = source_space.CreateField<pcms::Real>();
+  auto source_field = source_space->CreateFunction<pcms::Real>();
   pcms::test::SetField(
     source_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) {
       return x * x + x * y + 0.5 * y * y;
     });
 
   // Reference: exact mesh-intersection quadrature projection.
-  auto reference_field = target_space.CreateField<pcms::Real>();
-  pcms::OmegaHConservativeProjection reference_projection(source_space,
-                                                          target_space);
+  auto reference_field = target_space->CreateFunction<pcms::Real>();
+  pcms::OmegaHConservativeProjection reference_projection(*source_space,
+                                                          *target_space);
   reference_projection.Apply(source_field, reference_field);
   const auto reference =
     pcms::FlattenToRank1View(reference_field.GetDOFHolderDataHost());
@@ -212,20 +213,20 @@ TEST_CASE("OmegaHControlVariateProjection: reduces error vs plain Monte Carlo",
 
   // Plain Monte Carlo projection.
   pcms::OmegaHMonteCarloRHSIntegrator mc_integrator(
-    target_space, samples_per_element, pcms::MonteCarloSampling::UniformRandom,
+    *target_space, samples_per_element, pcms::MonteCarloSampling::UniformRandom,
     seed);
-  auto evaluator = source_space.CreatePointEvaluator<pcms::Real>(
+  auto evaluator = source_space->CreatePointEvaluator<pcms::Real>(
     pcms::EvaluationRequest::FromCoordinates(
       mc_integrator.GetIntegrationPoints().GetCoordinates()));
-  pcms::OmegaHMassIntegrator mass_integrator(target_space);
+  pcms::OmegaHMassIntegrator mass_integrator(*target_space);
   pcms::GalerkinProjectionSolver solver(mass_integrator, mc_integrator);
   const auto mc_projected = solver.Solve(*evaluator, source_field);
   const auto mc_projected_h = Omega_h::HostRead<Omega_h::Real>(mc_projected);
 
   // Control-variate projection with the same sampling configuration.
-  auto cv_field = target_space.CreateField<pcms::Real>();
+  auto cv_field = target_space->CreateFunction<pcms::Real>();
   pcms::OmegaHControlVariateProjection cv_projection(
-    source_space, target_space, samples_per_element,
+    *source_space, *target_space, samples_per_element,
     pcms::MonteCarloSampling::UniformRandom, seed);
   cv_projection.Apply(source_field, cv_field);
   const auto cv_values =
diff --git a/test/test_point_evaluator.cpp b/test/test_point_evaluator.cpp
index fa0fc2f5..dcc57087 100644
--- a/test/test_point_evaluator.cpp
+++ b/test/test_point_evaluator.cpp
@@ -31,16 +31,16 @@ TEST_CASE("PointEvaluator: OmegaH order-1 linear evaluation")
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
 
-  auto field_data = factory.CreateField<Real>();
+  auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
-    field_data.GetData(), *factory.GetLayout(),
+    field_data.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
 
   auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  auto evaluator = factory.CreatePointEvaluator<Real>(
+  auto evaluator = factory->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   pcms::test::CheckEvaluation(
     *evaluator, field_data, pts,
@@ -61,15 +61,15 @@ TEST_CASE("PointEvaluator: same evaluator reused for two FieldData objects")
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
 
-  auto field_a = factory.CreateField<Real>();
-  auto field_b = factory.CreateField<Real>();
+  auto field_a = factory->CreateFunction<Real>();
+  auto field_b = factory->CreateFunction<Real>();
 
   // field_a: linear_f;  field_b: constant 42
   pcms::test::SetField(
-    field_a.GetData(), *factory.GetLayout(),
+    field_a.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
   pcms::test::SetField(
-    field_b.GetData(), *factory.GetLayout(),
+    field_b.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real, Real) { return Real(42); });
 
   auto pts = pcms::test::StandardEvalCoords2D();
@@ -78,7 +78,7 @@ TEST_CASE("PointEvaluator: same evaluator reused for two FieldData objects")
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
 
   // Create the PointEvaluator once
-  auto evaluator = factory.CreatePointEvaluator<Real>(
+  auto evaluator = factory->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
   Kokkos::View<Real*, pcms::DeviceMemorySpace> out_a_device("out_a", n);
@@ -122,9 +122,9 @@ TEST_CASE("PointEvaluator: OmegaH order-1 out-of-bounds fill")
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
 
-  auto field_data = factory.CreateField<Real>();
+  auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
-    field_data.GetData(), *factory.GetLayout(),
+    field_data.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
 
   // Points clearly outside [0,1]^2
@@ -133,9 +133,9 @@ TEST_CASE("PointEvaluator: OmegaH order-1 out-of-bounds fill")
   auto device_coords = pcms::test::CreateDeviceCoordinateView(
     outside_pts, CoordinateSystem::Cartesian);
   pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0};
-  auto evaluator =
-    factory.CreatePointEvaluator<Real>(pcms::EvaluationRequest::FromCoordinates(
-      device_coords.coordinate_view, policy));
+  auto evaluator = factory->CreatePointEvaluator<Real>(
+    pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view,
+                                             policy));
   pcms::test::CheckFillMode(*evaluator, field_data, -999.0, outside_pts);
 }
 
@@ -155,15 +155,15 @@ TEST_CASE("PointEvaluator: UniformGrid order-1 linear evaluation")
   auto factory = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, CoordinateSystem::Cartesian, 1);
 
-  auto field_data = factory.CreateField<Real>();
+  auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
-    field_data.GetData(), *factory.GetLayout(),
+    field_data.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
   auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  auto evaluator = factory.CreatePointEvaluator<Real>(
+  auto evaluator = factory->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   pcms::test::CheckEvaluation(
     *evaluator, field_data, pts,
@@ -181,15 +181,15 @@ TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation")
   auto factory = pcms::SplineFunctionSpace::FromUniformGrid(
     grid, CoordinateSystem::Cartesian);
 
-  auto field_data = factory.CreateField<Real>();
+  auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
-    field_data.GetData(), *factory.GetLayout(),
+    field_data.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
   auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  auto evaluator = factory.CreatePointEvaluator<Real>(
+  auto evaluator = factory->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   pcms::test::CheckEvaluation(
     *evaluator, field_data, pts,
@@ -209,7 +209,7 @@ TEST_CASE("FieldLayout: metadata queries")
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto layout = factory.GetLayout();
+  auto layout = factory->GetLayout();
 
   auto coords = layout->GetDOFHolderCoordinates();
   REQUIRE(coords.GetCoordinateSystem() == CoordinateSystem::Cartesian);
@@ -230,9 +230,9 @@ TEST_CASE("FieldData: layout metadata queries")
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto field_data = factory.CreateField<Real>();
+  auto field_data = factory->CreateFunction<Real>();
 
-  auto coords = factory.GetLayout()->GetDOFHolderCoordinates();
+  auto coords = factory->GetLayout()->GetDOFHolderCoordinates();
   REQUIRE(coords.GetCoordinateSystem() == CoordinateSystem::Cartesian);
   REQUIRE(coords.GetValues().extent(0) > 0);
   REQUIRE(coords.GetValues().extent(1) == 2);
@@ -251,9 +251,9 @@ TEST_CASE("SimpleFieldData: set and get DOF holder data round-trip")
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto field_data = factory.CreateField<Real>();
+  auto field_data = factory->CreateFunction<Real>();
 
-  auto& layout = *factory.GetLayout();
+  auto& layout = *factory->GetLayout();
   int n = layout.GetNumOwnedDofHolder();
   REQUIRE(n > 0);
 
@@ -288,7 +288,7 @@ TEST_CASE("FieldLayout: MeshFields metadata queries")
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::MeshFields);
 
-  auto layout = factory.GetLayout();
+  auto layout = factory->GetLayout();
   auto coords = layout->GetDOFHolderCoordinates();
   REQUIRE(coords.GetCoordinateSystem() == CoordinateSystem::Cartesian);
   REQUIRE(coords.GetValues().extent(0) > 0);
@@ -305,16 +305,16 @@ TEST_CASE("PointEvaluator: MeshFields order-1 linear evaluation")
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::MeshFields);
 
-  auto field_data = factory.CreateField<Real>();
+  auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
-    field_data.GetData(), *factory.GetLayout(),
+    field_data.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
 
   auto pts = pcms::test::StandardEvalCoords2D();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  auto evaluator = factory.CreatePointEvaluator<Real>(
+  auto evaluator = factory->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   pcms::test::CheckEvaluation(
     *evaluator, field_data, pts,
@@ -331,9 +331,9 @@ TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill")
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::MeshFields);
 
-  auto field_data = factory.CreateField<Real>();
+  auto field_data = factory->CreateFunction<Real>();
   pcms::test::SetField(
-    field_data.GetData(), *factory.GetLayout(),
+    field_data.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
 
   const std::vector<Real> outside_pts = {-0.5, 0.5,  1.5, 0.5,
@@ -341,9 +341,9 @@ TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill")
   auto device_coords = pcms::test::CreateDeviceCoordinateView(
     outside_pts, CoordinateSystem::Cartesian);
   pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0};
-  auto evaluator =
-    factory.CreatePointEvaluator<Real>(pcms::EvaluationRequest::FromCoordinates(
-      device_coords.coordinate_view, policy));
+  auto evaluator = factory->CreatePointEvaluator<Real>(
+    pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view,
+                                             policy));
   pcms::test::CheckFillMode(*evaluator, field_data, -999.0, outside_pts);
 }
 
@@ -358,13 +358,13 @@ TEST_CASE(
     mesh, 1, 1, CoordinateSystem::Cartesian, "global",
     pcms::LagrangeFunctionSpace::Backend::MeshFields);
 
-  auto field_a = factory.CreateField<Real>();
-  auto field_b = factory.CreateField<Real>();
+  auto field_a = factory->CreateFunction<Real>();
+  auto field_b = factory->CreateFunction<Real>();
   pcms::test::SetField(
-    field_a.GetData(), *factory.GetLayout(),
+    field_a.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
   pcms::test::SetField(
-    field_b.GetData(), *factory.GetLayout(),
+    field_b.GetData(), *factory->GetLayout(),
     OMEGA_H_LAMBDA(Real, Real) { return Real(42); });
 
   auto pts = pcms::test::StandardEvalCoords2D();
@@ -372,7 +372,7 @@ TEST_CASE(
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
 
-  auto evaluator = factory.CreatePointEvaluator<Real>(
+  auto evaluator = factory->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
   Kokkos::View<Real*, pcms::DeviceMemorySpace> out_a_device("out_a", n);
diff --git a/test/test_polynomial_reconstruction_function_space.cpp b/test/test_polynomial_reconstruction_function_space.cpp
index 1f3b675a..0281c13f 100644
--- a/test/test_polynomial_reconstruction_function_space.cpp
+++ b/test/test_polynomial_reconstruction_function_space.cpp
@@ -37,7 +37,7 @@ TEST_CASE(
 
   auto factory = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian);
-  auto layout = factory.GetLayout();
+  auto layout = factory->GetLayout();
 
   REQUIRE(layout->GetNumComponents() == 1);
   REQUIRE(layout->GetNumOwnedDofHolder() == 4);
@@ -64,8 +64,8 @@ TEST_CASE("PolynomialReconstructionFunctionSpace fields share layout")
 
   auto factory = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian);
-  auto source = factory.CreateField<Real>(pcms::FieldMetadata{});
-  auto target = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto source = factory->CreateFunction<Real>();
+  auto target = factory->CreateFunction<Real>();
 
   REQUIRE(&source.GetLayout() == &target.GetLayout());
 }
@@ -78,7 +78,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace point-cloud field set/get DOF "
 
   auto field = pcms::PolynomialReconstructionFunctionSpace::Create(
                  coords_view, CoordinateSystem::Cartesian)
-                 .CreateField<Real>(pcms::FieldMetadata{});
+                 ->CreateFunction<Real>();
 
   std::vector<Real> data{1.0, 2.0, 3.0, 4.0};
   Rank2View<const Real, HostMemorySpace> data_view(
@@ -100,14 +100,14 @@ TEST_CASE("PolynomialReconstructionFunctionSpace point-cloud field serialize / "
 
   auto factory = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<Real>();
 
   std::vector<Real> data{5.0, 6.0, 7.0, 8.0};
   Rank2View<const Real, HostMemorySpace> data_view(
     data.data(), static_cast<LO>(data.size()), 1);
   field.GetData().SetDOFHolderDataHost(data_view);
 
-  pcms::test::CheckSerializeDeserialize(*factory.GetLayout(), field.GetData());
+  pcms::test::CheckSerializeDeserialize(*factory->GetLayout(), field.GetData());
 }
 
 TEST_CASE("PolynomialReconstructionFunctionSpace field keeps layout alive "
@@ -119,7 +119,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace field keeps layout alive "
   auto field = [&]() {
     auto factory = pcms::PolynomialReconstructionFunctionSpace::Create(
       coords_view, CoordinateSystem::Cartesian);
-    return factory.CreateField<Real>(pcms::FieldMetadata{});
+    return factory->CreateFunction<Real>();
   }();
 
   auto point_cloud_layout =
@@ -148,15 +148,15 @@ TEST_CASE("Different layouts on the same mesh report SameEntities")
   auto lagrange = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, CoordinateSystem::Cartesian);
 
-  auto nodal_disc = nodal.GetLayout()->GetDiscretization();
-  auto lagrange_disc = lagrange.GetLayout()->GetDiscretization();
+  auto nodal_disc = nodal->GetLayout()->GetDiscretization();
+  auto lagrange_disc = lagrange->GetLayout()->GetDiscretization();
 
   REQUIRE(nodal_disc != nullptr);
   REQUIRE(lagrange_disc != nullptr);
   REQUIRE(nodal_disc->SameEntities(*lagrange_disc));
 
-  REQUIRE(nodal.GetLayout()->GetNumOwnedDofHolder() == mesh.nfaces());
-  REQUIRE(lagrange.GetLayout()->GetNumOwnedDofHolder() == mesh.nverts());
+  REQUIRE(nodal->GetLayout()->GetNumOwnedDofHolder() == mesh.nfaces());
+  REQUIRE(lagrange->GetLayout()->GetNumOwnedDofHolder() == mesh.nverts());
 
   REQUIRE(nodal_disc->GetNumEntities(pcms::Face) == mesh.nfaces());
   REQUIRE(lagrange_disc->GetNumEntities(pcms::Vertex) == mesh.nverts());
@@ -175,8 +175,8 @@ TEST_CASE("Layouts on different meshes do not report SameEntities")
   auto nodal_b = pcms::PolynomialReconstructionFunctionSpace::FromMesh(
     mesh_b, pcms::Vertex, CoordinateSystem::Cartesian);
 
-  auto disc_a = nodal_a.GetLayout()->GetDiscretization();
-  auto disc_b = nodal_b.GetLayout()->GetDiscretization();
+  auto disc_a = nodal_a->GetLayout()->GetDiscretization();
+  auto disc_b = nodal_b->GetLayout()->GetDiscretization();
 
   REQUIRE_FALSE(disc_a->SameEntities(*disc_b));
 }
@@ -196,8 +196,8 @@ TEST_CASE(
   auto standalone = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian);
 
-  auto mesh_disc = lagrange.GetLayout()->GetDiscretization();
-  auto point_cloud_disc = standalone.GetLayout()->GetDiscretization();
+  auto mesh_disc = lagrange->GetLayout()->GetDiscretization();
+  auto point_cloud_disc = standalone->GetLayout()->GetDiscretization();
 
   REQUIRE_FALSE(mesh_disc->SameEntities(*point_cloud_disc));
   REQUIRE_FALSE(point_cloud_disc->SameEntities(*mesh_disc));
diff --git a/test/test_polynomial_reconstruction_mls_evaluation.cpp b/test/test_polynomial_reconstruction_mls_evaluation.cpp
index 8b979e68..888a83b1 100644
--- a/test/test_polynomial_reconstruction_mls_evaluation.cpp
+++ b/test/test_polynomial_reconstruction_mls_evaluation.cpp
@@ -119,13 +119,13 @@ void CheckPolynomialReproduction(unsigned degree,
 
   auto fs = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian, SweepTestOptions(degree, basis));
-  auto field = fs.CreateField<Real>(pcms::FieldMetadata{});
-  pcms::test::SetField(field.GetData(), *fs.GetLayout(), func);
+  auto field = fs->CreateFunction<Real>();
+  pcms::test::SetField(field.GetData(), *fs->GetLayout(), func);
 
   auto pts = QueryPoints();
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  auto evaluator = fs.CreatePointEvaluator<Real>(
+  auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   pcms::test::CheckEvaluation(*evaluator, field, pts, func, abs_tol);
 }
@@ -170,22 +170,22 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: same PointEvaluator "
 
   auto fs = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian, DefaultTestOptions());
-  auto field_a = fs.CreateField<Real>(pcms::FieldMetadata{});
-  auto field_b = fs.CreateField<Real>(pcms::FieldMetadata{});
+  auto field_a = fs->CreateFunction<Real>();
+  auto field_b = fs->CreateFunction<Real>();
 
   pcms::test::SetField(
-    field_a.GetData(), *fs.GetLayout(),
+    field_a.GetData(), *fs->GetLayout(),
     OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; });
   const Real cval = 7.0;
   pcms::test::SetField(
-    field_b.GetData(), *fs.GetLayout(),
+    field_b.GetData(), *fs->GetLayout(),
     OMEGA_H_LAMBDA(Real, Real) { return cval; });
 
   auto pts = QueryPoints();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  auto evaluator = fs.CreatePointEvaluator<Real>(
+  auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
   Kokkos::View<Real*, DeviceMemorySpace> out_a_device("out_a", n);
@@ -226,13 +226,13 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: Evaluate throws for "
 
   auto fs = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian);
-  auto field = fs.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = fs->CreateFunction<Real>();
 
   auto pts = QueryPoints();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  auto evaluator = fs.CreatePointEvaluator<Real>(
+  auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
   // Two-component output — must throw
@@ -258,16 +258,16 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: default MLSOptions — "
   // Use default options — no third argument
   auto fs = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian);
-  auto field = fs.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = fs->CreateFunction<Real>();
   pcms::test::SetField(
-    field.GetData(), *fs.GetLayout(),
+    field.GetData(), *fs->GetLayout(),
     OMEGA_H_LAMBDA(Real, Real) { return Real(1.0); });
 
   auto pts = QueryPoints();
   int n = static_cast<int>(pts.size()) / 2;
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  auto evaluator = fs.CreatePointEvaluator<Real>(
+  auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
 
   Kokkos::View<Real*, DeviceMemorySpace> out_device("out", n);
@@ -296,7 +296,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator "
   auto pts = QueryPoints();
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  REQUIRE_THROWS(fs.CreatePointEvaluator<Real>(
+  REQUIRE_THROWS(fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)));
 }
 
@@ -313,7 +313,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator "
   auto pts = QueryPoints();
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cylindrical);
-  REQUIRE_THROWS(fs.CreatePointEvaluator<Real>(
+  REQUIRE_THROWS(fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)));
 }
 
@@ -332,7 +332,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: radius option is "
 
   auto fs = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian, opts);
-  auto field = fs.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = fs->CreateFunction<Real>();
   std::vector<Real> dof_values{1.0, 5.0};
   field.GetData().SetDOFHolderDataHost(Rank2View<const Real, HostMemorySpace>(
     dof_values.data(), static_cast<pcms::LO>(dof_values.size()), 1));
@@ -340,7 +340,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: radius option is "
   std::vector<Real> pts{0.0, 0.0};
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian);
-  auto evaluator = fs.CreatePointEvaluator<Real>(
+  auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   Kokkos::View<Real*, DeviceMemorySpace> out_device("out", 1);
   using LayoutPolicy =
@@ -364,7 +364,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve "
   auto fs = pcms::PolynomialReconstructionFunctionSpace::Create(
     coords_view, CoordinateSystem::Cartesian, DefaultTestOptions3D());
   auto layout_coords_device =
-    fs.GetLayout()->GetDOFHolderCoordinates().GetValues();
+    fs->GetLayout()->GetDOFHolderCoordinates().GetValues();
   auto layout_coords =
     pcms::test::CopyCoordinatesToHost(layout_coords_device, 27, 3);
   REQUIRE(static_cast<int>(layout_coords.extent(1)) == 3);
@@ -374,7 +374,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve "
     REQUIRE(layout_coords(i, 2) == Catch::Approx(src[3 * i + 2]));
   }
 
-  auto field = fs.CreateField<Real>(pcms::FieldMetadata{});
+  auto field = fs->CreateFunction<Real>();
   std::vector<Real> dof_values(27);
   for (int i = 0; i < 27; ++i) {
     const Real x = src[3 * i + 0];
@@ -388,7 +388,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve "
   std::vector<Real> pts{0.5, 0.5, 0.25, 0.5, 0.5, 0.75};
   auto device_coords =
     pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian, 3);
-  auto evaluator = fs.CreatePointEvaluator<Real>(
+  auto evaluator = fs->CreatePointEvaluator<Real>(
     pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view));
   Kokkos::View<Real*, DeviceMemorySpace> out_device("out", 2);
   using LayoutPolicy =
diff --git a/test/test_proxy_coupling.cpp b/test/test_proxy_coupling.cpp
index 455fc183..6e166ad6 100644
--- a/test/test_proxy_coupling.cpp
+++ b/test/test_proxy_coupling.cpp
@@ -90,20 +90,20 @@ void xgc_delta_f(MPI_Comm comm, Omega_h::Mesh& mesh)
   pcms::Application* app = coupler.AddApplication("proxy_couple_xgc_delta_f");
 
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
-    mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  app->AddLayout("gids", factory.GetLayout());
+    mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
+    pcms::LagrangeFunctionSpace::DefaultBackend, "gids");
 
-  auto gids_field = factory.CreateField<pcms::Real>(pcms::FieldMetadata{});
-  auto gids2_field = factory.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto gids_field = factory->CreateFunction<pcms::Real>("gids");
+  auto gids2_field = factory->CreateFunction<pcms::Real>("gids2");
 
   auto* gids_ptr = &gids_field.GetData();
   auto* gids2_ptr = &gids2_field.GetData();
 
-  initializeFieldWithGids(*factory.GetLayout(), gids_ptr, 1.0);
-  initializeFieldWithGids(*factory.GetLayout(), gids2_ptr, 2.0);
+  initializeFieldWithGids(*factory->GetLayout(), gids_ptr, 1.0);
+  initializeFieldWithGids(*factory->GetLayout(), gids2_ptr, 2.0);
 
-  app->AddField("gids", std::move(gids_field));
-  app->AddField("gids2", std::move(gids2_field));
+  app->AddField(std::move(gids_field));
+  app->AddField(std::move(gids2_field));
 
   do {
     for (int i = 0; i < COMM_ROUNDS; ++i) {
@@ -114,7 +114,7 @@ void xgc_delta_f(MPI_Comm comm, Omega_h::Mesh& mesh)
       app->ReceiveField("gids"); //(Alt) df_gid_field->Receive();
       app->EndReceivePhase();
 
-      if (!validateField(*factory.GetLayout(), gids_ptr, "gids", rank, 1.0)) {
+      if (!validateField(*factory->GetLayout(), gids_ptr, "gids", rank, 1.0)) {
         std::cerr << "xgc_delta_f: Field validation failed at round " << i
                   << std::endl;
         exit(EXIT_FAILURE);
@@ -131,16 +131,16 @@ void xgc_total_f(MPI_Comm comm, Omega_h::Mesh& mesh)
   pcms::Application* app = coupler.AddApplication("proxy_couple_xgc_total_f");
 
   auto factory = pcms::LagrangeFunctionSpace::FromMesh(
-    mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  app->AddLayout("gids", factory.GetLayout());
+    mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
+    pcms::LagrangeFunctionSpace::DefaultBackend, "gids");
 
-  auto gids_field = factory.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto gids_field = factory->CreateFunction<pcms::Real>("gids");
 
   auto* gids_ptr = &gids_field.GetData();
 
-  initializeFieldWithGids(*factory.GetLayout(), gids_ptr, 10.0);
+  initializeFieldWithGids(*factory->GetLayout(), gids_ptr, 10.0);
 
-  app->AddField("gids", std::move(gids_field));
+  app->AddField(std::move(gids_field));
 
   do {
     for (int i = 0; i < COMM_ROUNDS; ++i) {
@@ -151,7 +151,7 @@ void xgc_total_f(MPI_Comm comm, Omega_h::Mesh& mesh)
       app->ReceiveField("gids"); //(Alt) tf_gid_field->Receive();
       app->EndReceivePhase();
 
-      if (!validateField(*factory.GetLayout(), gids_ptr, "gids", rank, 10.0)) {
+      if (!validateField(*factory->GetLayout(), gids_ptr, "gids", rank, 10.0)) {
         std::cerr << "xgc_total_f: Field validation failed at round " << i
                   << std::endl;
         exit(EXIT_FAILURE);
@@ -173,27 +173,24 @@ void xgc_coupler(MPI_Comm comm, Omega_h::Mesh& mesh, std::string_view cpn_file)
   auto* total_f = cpl.AddApplication("proxy_couple_xgc_total_f");
   auto* delta_f = cpl.AddApplication("proxy_couple_xgc_delta_f");
   auto factory_total = pcms::LagrangeFunctionSpace::FromMesh(
-    mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  total_f->AddLayout("gids", factory_total.GetLayout());
+    mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
+    pcms::LagrangeFunctionSpace::DefaultBackend, "gids");
 
   auto factory_delta = pcms::LagrangeFunctionSpace::FromMesh(
-    mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  delta_f->AddLayout("gids", factory_delta.GetLayout());
+    mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global",
+    pcms::LagrangeFunctionSpace::DefaultBackend, "gids");
   // TODO, fields should have a transfer policy rather than parameters
-  auto total_gids_field =
-    factory_total.CreateField<pcms::Real>(pcms::FieldMetadata{});
-  auto delta_gids_field =
-    factory_delta.CreateField<pcms::Real>(pcms::FieldMetadata{});
-  auto delta_gids2_field =
-    factory_delta.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto total_gids_field = factory_total->CreateFunction<pcms::Real>("gids");
+  auto delta_gids_field = factory_delta->CreateFunction<pcms::Real>("gids");
+  auto delta_gids2_field = factory_delta->CreateFunction<pcms::Real>("gids2");
 
   auto* total_gids_ptr = &total_gids_field.GetData();
   auto* delta_gids_ptr = &delta_gids_field.GetData();
   auto* delta_gids2_ptr = &delta_gids2_field.GetData();
 
-  total_f->AddField("gids", std::move(total_gids_field));
-  delta_f->AddField("gids", std::move(delta_gids_field));
-  delta_f->AddField("gids2", std::move(delta_gids2_field));
+  total_f->AddField(std::move(total_gids_field));
+  delta_f->AddField(std::move(delta_gids_field));
+  delta_f->AddField(std::move(delta_gids2_field));
 
   do {
     for (int i = 0; i < COMM_ROUNDS; ++i) {
@@ -201,7 +198,7 @@ void xgc_coupler(MPI_Comm comm, Omega_h::Mesh& mesh, std::string_view cpn_file)
       total_f->ReceiveField("gids");
       total_f->EndReceivePhase();
 
-      if (!validateField(*factory_total.GetLayout(), total_gids_ptr, "gids",
+      if (!validateField(*factory_total->GetLayout(), total_gids_ptr, "gids",
                          rank, 10.0)) {
         std::cerr << "xgc_coupler: total_f field validation failed at round "
                   << i << std::endl;
@@ -212,7 +209,7 @@ void xgc_coupler(MPI_Comm comm, Omega_h::Mesh& mesh, std::string_view cpn_file)
       delta_f->ReceiveField("gids");
       delta_f->EndReceivePhase();
 
-      if (!validateField(*factory_delta.GetLayout(), delta_gids_ptr, "gids",
+      if (!validateField(*factory_delta->GetLayout(), delta_gids_ptr, "gids",
                          rank, 1.0)) {
         std::cerr << "xgc_coupler: delta_f field validation failed at round "
                   << i << std::endl;
diff --git a/test/test_proxy_coupling_xgc_server.cpp b/test/test_proxy_coupling_xgc_server.cpp
index 4a76b748..ada0b759 100644
--- a/test/test_proxy_coupling_xgc_server.cpp
+++ b/test/test_proxy_coupling_xgc_server.cpp
@@ -51,17 +51,17 @@ void xgc_coupler(MPI_Comm comm, Omega_h::Mesh& mesh, std::string_view cpn_file)
     // FIXME: The current C/Fortran proxy API couples layout registration to
     // field registration, so each XGC plane is registered as a separate layout
     // communicator even though the layouts are geometrically identical.
-    auto function_space = pcms::XGCFunctionSpace(
-      rc, ts::IsModelEntInOverlap{}, static_cast<pcms::LO>(mesh.nverts()));
+    auto function_space =
+      pcms::XGCFieldFactory(rc, ts::IsModelEntInOverlap{},
+                            static_cast<pcms::LO>(mesh.nverts()), ss.str());
     auto field = function_space.CreateField<pcms::GO>(
-      std::make_unique<pcms::XGCFieldData<pcms::GO>>(
-        function_space.GetXGCLayout(), pcms::FieldMetadata{},
-        make_array_view(data[i])));
-    application->AddLayout(ss.str(), function_space.GetLayout());
+      ss.str(), std::make_unique<pcms::XGCFieldData<pcms::GO>>(
+                  function_space.GetXGCLayout(), pcms::FieldMetadata{},
+                  make_array_view(data[i])));
     std::unique_ptr<pcms::FieldSerializer<GO>> serializer =
       std::make_unique<pcms::XGCFieldSerializer<GO>>(comm);
     fields.push_back(
-      application->AddField(ss.str(), std::move(field), std::move(serializer)));
+      application->AddField(std::move(field), std::move(serializer)));
   }
 
   do {
@@ -172,15 +172,14 @@ void omegah_coupler(MPI_Comm comm, Omega_h::Mesh& mesh,
     // communicator even though the layouts are geometrically identical.
     auto factory = pcms::LagrangeFunctionSpace::FromMesh(
       mesh, 1, 1, pcms::CoordinateSystem::Cartesian, numbering,
-      pcms::LagrangeFunctionSpace::Backend::OmegaH);
-    application->AddLayout(ss.str(), factory.GetLayout());
-    auto field =
-      factory.CreateField<GO>(std::make_unique<pcms::SimpleFieldData<GO>>(
-        factory.GetLayout(), pcms::FieldMetadata{}));
+      pcms::LagrangeFunctionSpace::Backend::OmegaH, ss.str());
+    auto field = factory->CreateFunction<GO>(
+      ss.str(), std::make_unique<pcms::SimpleFieldData<GO>>(
+                  factory->GetLayout(), pcms::FieldMetadata{}));
     std::unique_ptr<pcms::FieldSerializer<GO>> serializer =
       std::make_unique<pcms::FieldSerializer<GO>>();
     fields.push_back(
-      application->AddField(ss.str(), std::move(field), std::move(serializer)));
+      application->AddField(std::move(field), std::move(serializer)));
   }
   do {
     application->ReceivePhase([&]() {
diff --git a/test/test_proxy_coupling_xgc_server_overlap.cpp b/test/test_proxy_coupling_xgc_server_overlap.cpp
index 0161ee4e..4bfafc20 100644
--- a/test/test_proxy_coupling_xgc_server_overlap.cpp
+++ b/test/test_proxy_coupling_xgc_server_overlap.cpp
@@ -54,21 +54,20 @@ void xgc_coupler_with_overlap(MPI_Comm comm, Omega_h::Mesh& mesh,
 
     application->SetLayoutOverlapMask(ss.str(), std::move(overlap_mask));
 
-    auto function_space = pcms::XGCFunctionSpace(
-      rc, ts::IsModelEntInOverlap{}, static_cast<pcms::LO>(mesh.nverts()));
-
-    application->AddLayout(ss.str(), function_space.GetLayout());
+    auto function_space =
+      pcms::XGCFieldFactory(rc, ts::IsModelEntInOverlap{},
+                            static_cast<pcms::LO>(mesh.nverts()), ss.str());
 
     auto field = function_space.CreateField<pcms::GO>(
-      std::make_unique<pcms::XGCFieldData<pcms::GO>>(
-        function_space.GetXGCLayout(), pcms::FieldMetadata{},
-        make_array_view(data[i])));
+      ss.str(), std::make_unique<pcms::XGCFieldData<pcms::GO>>(
+                  function_space.GetXGCLayout(), pcms::FieldMetadata{},
+                  make_array_view(data[i])));
 
     std::unique_ptr<pcms::FieldSerializer<GO>> serializer =
       std::make_unique<pcms::XGCFieldSerializer<GO>>(comm);
 
     fields.push_back(
-      application->AddField(ss.str(), std::move(field), std::move(serializer)));
+      application->AddField(std::move(field), std::move(serializer)));
   }
 
   do {
@@ -175,19 +174,17 @@ void omegah_coupler_with_overlap(MPI_Comm comm, Omega_h::Mesh& mesh,
 
     auto factory = pcms::LagrangeFunctionSpace::FromMesh(
       mesh, 1, 1, pcms::CoordinateSystem::Cartesian, numbering,
-      pcms::LagrangeFunctionSpace::Backend::OmegaH);
-
-    application->AddLayout(ss.str(), factory.GetLayout());
+      pcms::LagrangeFunctionSpace::Backend::OmegaH, ss.str());
 
-    auto field =
-      factory.CreateField<GO>(std::make_unique<pcms::SimpleFieldData<GO>>(
-        factory.GetLayout(), pcms::FieldMetadata{}));
+    auto field = factory->CreateFunction<GO>(
+      ss.str(), std::make_unique<pcms::SimpleFieldData<GO>>(
+                  factory->GetLayout(), pcms::FieldMetadata{}));
 
     std::unique_ptr<pcms::FieldSerializer<GO>> serializer =
       std::make_unique<pcms::FieldSerializer<GO>>();
 
     fields.push_back(
-      application->AddField(ss.str(), std::move(field), std::move(serializer)));
+      application->AddField(std::move(field), std::move(serializer)));
   }
 
   do {
diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp
index b179167d..632bd4db 100644
--- a/test/test_uniform_grid_field.cpp
+++ b/test/test_uniform_grid_field.cpp
@@ -109,7 +109,7 @@ TEST_CASE("UniformGrid field creation")
 
   auto field_space = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = field_space.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto field = field_space->CreateFunction<pcms::Real>();
   REQUIRE(field.GetDOFHolderDataHost().size() ==
           static_cast<size_t>(layout->OwnedSize()));
 }
@@ -125,7 +125,7 @@ TEST_CASE("UniformGrid order-0 field creation and evaluation")
     grid, 1, pcms::CoordinateSystem::Cartesian, 0);
   auto field_space = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, pcms::CoordinateSystem::Cartesian, 0);
-  auto field = field_space.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto field = field_space->CreateFunction<pcms::Real>();
   pcms::UniformGridEvaluatorFactory<2> eval_factory(layout);
 
   REQUIRE(layout->GetOrder() == 0);
@@ -179,7 +179,7 @@ TEST_CASE("UniformGrid field data operations", "[uniform_grid_field]")
     grid, 1, pcms::CoordinateSystem::Cartesian);
   auto field_space = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = field_space.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto field = field_space->CreateFunction<pcms::Real>();
 
   std::vector<pcms::Real> data(25);
   for (size_t i = 0; i < 25; ++i)
@@ -206,7 +206,7 @@ TEST_CASE("UniformGrid field evaluation - piecewise constant")
     grid, 1, pcms::CoordinateSystem::Cartesian);
   auto field_space = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = field_space.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto field = field_space->CreateFunction<pcms::Real>();
   pcms::UniformGridEvaluatorFactory<2> eval_factory(layout);
 
   // Set vertex values for a 2x2 cell grid (3x3 = 9 vertices)
@@ -269,7 +269,7 @@ TEST_CASE("UniformGrid field serialization")
     grid, 1, pcms::CoordinateSystem::Cartesian);
   auto field_space = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = field_space.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto field = field_space->CreateFunction<pcms::Real>();
 
   std::vector<pcms::Real> data(16);
   for (size_t i = 0; i < 16; ++i)
@@ -303,13 +303,13 @@ TEST_CASE("UniformGrid field copy")
 
   auto factory = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, pcms::CoordinateSystem::Cartesian);
-  auto field = factory.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto field = factory->CreateFunction<pcms::Real>();
   field.SetDOFHolderDataHost(
     pcms::Rank2View<const pcms::Real, pcms::HostMemorySpace>(
       data.data(), static_cast<pcms::LO>(data.size()), 1));
 
-  auto field2 = factory.CreateField<pcms::Real>(pcms::FieldMetadata{});
-  pcms::Copy<pcms::Real> copy(factory, factory);
+  auto field2 = factory->CreateFunction<pcms::Real>();
+  pcms::Copy<pcms::Real> copy(*factory, *factory);
   copy.Apply(field, field2);
 
   auto copied_data = pcms::FlattenToRank1View(field2.GetDOFHolderDataHost());
@@ -326,8 +326,7 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field")
 
   auto omega_h_factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto omega_h_field =
-    omega_h_factory.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto omega_h_field = omega_h_factory->CreateFunction<pcms::Real>();
   pcms::test::SetField(
     omega_h_field,
     OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; });
@@ -338,15 +337,15 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field")
   grid.divisions = {2, 2};
   auto ug_factory = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, pcms::CoordinateSystem::Cartesian);
-  auto ug_field = ug_factory.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto ug_field = ug_factory->CreateFunction<pcms::Real>();
 
-  pcms::Interpolator<pcms::Real> interp(omega_h_factory, ug_factory);
+  pcms::Interpolator<pcms::Real> interp(*omega_h_factory, *ug_factory);
   interp.Apply(omega_h_field, ug_field);
 
   auto transferred_data =
     pcms::FlattenToRank1View(ug_field.GetDOFHolderDataHost());
-  auto ug_coords = ug_factory.GetLayout()->GetDOFHolderCoordinates();
-  int num_ug_nodes = ug_factory.GetLayout()->GetNumOwnedDofHolder();
+  auto ug_coords = ug_factory->GetLayout()->GetDOFHolderCoordinates();
+  int num_ug_nodes = ug_factory->GetLayout()->GetNumOwnedDofHolder();
 
   // set up_coords to host
   auto ug_coords_host =
@@ -640,23 +639,22 @@ TEST_CASE("UniformGrid workflow")
 
   auto omega_h_factory = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian);
-  auto omega_h_field =
-    omega_h_factory.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto omega_h_field = omega_h_factory->CreateFunction<pcms::Real>();
   pcms::test::SetField(
     omega_h_field,
     OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; });
 
   auto ug_factory = pcms::LagrangeFunctionSpace::FromUniformGrid(
     grid, 1, pcms::CoordinateSystem::Cartesian);
-  auto ug_field = ug_factory.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto ug_field = ug_factory->CreateFunction<pcms::Real>();
 
   auto [mask_layout, mask_field] =
     pcms::CreateUniformGridBinaryField<2>(mesh, grid);
 
-  pcms::Interpolator<pcms::Real> interp(omega_h_factory, ug_factory);
+  pcms::Interpolator<pcms::Real> interp(*omega_h_factory, *ug_factory);
   interp.Apply(omega_h_field, ug_field);
   auto ug_coords_device_view =
-    ug_factory.GetLayout()->GetDOFHolderCoordinates().GetValues();
+    ug_factory->GetLayout()->GetDOFHolderCoordinates().GetValues();
   auto ug_coords_host_view =
     pcms::test::CopyCoordinatesToHost(ug_coords_device_view, 25, 2);
 
diff --git a/test/test_xgc_field_data.cpp b/test/test_xgc_field_data.cpp
index 234bd2d1..30a95163 100644
--- a/test/test_xgc_field_data.cpp
+++ b/test/test_xgc_field_data.cpp
@@ -104,29 +104,18 @@ TEST_CASE("XGC FieldData serializer preserves inactive entries")
   }
 }
 
-TEST_CASE("XGCFunctionSpace creates fields and rejects evaluator access")
+TEST_CASE("XGCFieldFactory creates fields and rejects evaluator access")
 {
   static constexpr int data_size = 16;
   auto rc = create_dummy_rc(data_size);
-  pcms::XGCFunctionSpace function_space(rc, in_overlap, data_size);
+  pcms::XGCFieldFactory function_space(rc, in_overlap, data_size);
 
   std::vector<pcms::Real> data(data_size);
   std::iota(data.begin(), data.end(), 0.0);
   auto field = function_space.CreateField<pcms::Real>(
-    std::make_unique<pcms::XGCFieldData<pcms::Real>>(
-      function_space.GetXGCLayout(), pcms::FieldMetadata{},
-      pcms::make_array_view(data)));
+    "", std::make_unique<pcms::XGCFieldData<pcms::Real>>(
+          function_space.GetXGCLayout(), pcms::FieldMetadata{},
+          pcms::make_array_view(data)));
 
   REQUIRE(&field.GetLayout() == function_space.GetLayout().get());
-  REQUIRE(function_space.GetCoordinateSystem() == pcms::CoordinateSystem::XGC);
-  // XGC does not support point evaluation; CreatePointEvaluator throws.
-  using LayoutPolicy =
-    pcms::detail::default_layout_for_memory_space_t<pcms::DeviceMemorySpace>;
-  pcms::Rank2View<const pcms::Real, pcms::DeviceMemorySpace, LayoutPolicy>
-    empty_coords{nullptr, 0, 2};
-  REQUIRE_THROWS_AS(function_space.CreatePointEvaluator<pcms::Real>(
-                      pcms::EvaluationRequest::FromCoordinates(
-                        pcms::CoordinateView<pcms::DeviceMemorySpace>{
-                          pcms::CoordinateSystem::XGC, empty_coords})),
-                    pcms::pcms_error);
 }
diff --git a/test/xgc_n0_coupling_server.cpp b/test/xgc_n0_coupling_server.cpp
index ddbac140..0484c810 100644
--- a/test/xgc_n0_coupling_server.cpp
+++ b/test/xgc_n0_coupling_server.cpp
@@ -43,16 +43,15 @@ struct RegisteredField
 [[nodiscard]]
 static RegisteredField AddField(
   pcms::Application* application,
-  const pcms::LagrangeFunctionSpace& function_space, const std::string& name,
-  const std::string& path, int plane)
+  const std::shared_ptr<pcms::LagrangeFunctionSpace>& function_space,
+  const std::string& name, const std::string& path, int plane)
 {
   PCMS_ALWAYS_ASSERT(application != nullptr);
   auto field_name = MakeFieldName(name, plane);
-  auto field = function_space.CreateField<pcms::Real>(pcms::FieldMetadata{});
+  auto field = function_space->CreateFunction<pcms::Real>(path + field_name);
   std::unique_ptr<pcms::FieldSerializer<pcms::Real>> serializer =
     std::make_unique<pcms::FieldSerializer<pcms::Real>>();
-  auto handle = application->AddField(path + field_name, std::move(field),
-                                      std::move(serializer));
+  auto handle = application->AddField(std::move(field), std::move(serializer));
   return {std::move(handle)};
 }
 
@@ -220,10 +219,7 @@ void omegah_coupler(MPI_Comm comm, Omega_h::Mesh& mesh,
     });
   auto function_space = pcms::LagrangeFunctionSpace::FromMesh(
     mesh, 1, 1, pcms::CoordinateSystem::Cartesian, is_overlap, numbering,
-    pcms::LagrangeFunctionSpace::Backend::OmegaH);
-  auto layout = function_space.GetLayout();
-  core->AddLayout("core_layout", layout);
-  edge->AddLayout("edge_layout", layout);
+    pcms::LagrangeFunctionSpace::Backend::OmegaH, "n0_layout");
   auto time2 = std::chrono::steady_clock::now();
   elapsed_seconds = time2 - time1;
   ts::timeMinMaxAvg(elapsed_seconds.count(), min, max, avg);