Add mfem adapter

CMakeLists.txt

@@ -44,6 +44,7 @@ option(PCMS_ENABLE_CLIENT "enable the coupling client implementation" ON)
 
 option(PCMS_ENABLE_XGC "enable xgc field adapter" ON)
 option(PCMS_ENABLE_OMEGA_H "enable Omega_h field adapter" OFF)
+option(PCMS_ENABLE_MFEM "enable MFEM field adapter" OFF)
 option(PCMS_ENABLE_C "Enable pcms C api" ON)
 option(PCMS_ENABLE_Python "Enable pcms Python api" OFF)
 
@@ -119,6 +120,14 @@ if(PCMS_ENABLE_OMEGA_H)
     message(FATAL_ERROR "Omega_h must be built with MPI enabled.")
   endif()
 endif()
+
+if(PCMS_ENABLE_MFEM)
+  find_package(MFEM REQUIRED)
+  message(STATUS "Found MFEM: ${MFEM_DIR} (found version ${MFEM_VERSION})")
+  if(NOT MFEM_USE_MPI)
+    message(FATAL_ERROR "MFEM must be built with MPI enabled.")
+  endif()
+endif()
 # adios2 adds C and Fortran depending on how it was built
 find_package(ADIOS2 CONFIG 2.10.2 REQUIRED)
 find_package(Kokkos CONFIG 4.5 REQUIRED)

src/pcms/capi/client.cpp

@@ -1,5 +1,6 @@
 #include "client.h"
 #include "pcms.h"
+#include "pcms/field/function_space.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)
   {
@@ -251,7 +252,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>>(

src/pcms/coupler/field_exchange_planner.cpp

@@ -12,6 +12,10 @@ namespace pcms
 namespace
 {
 
+// Sentinel in a receive permutation: this local DOF's GID was not present in
+// the received message. The deserializer skips it (negative => not received).
+constexpr redev::LO kUnreceivedDof = -1;
+
 struct PartitionMapping
 {
   std::vector<LO> indices;
@@ -134,8 +138,16 @@ static redev::LOs ConstructPermutation(
     const auto start = ent_offsets[e];
     const auto end = ent_offsets[e + 1];
 
-    for (size_t i = start; i < end; ++i)
-      permutation.push_back(gid_to_buffer_index[e][local_gids[i]]);
+    for (size_t i = start; i < end; ++i) {
+      // A local DOF whose GID was not in the received message (e.g. it lies
+      // outside the sender's overlap mask) gets the sentinel kUnreceivedDof.
+      // The deserializer must skip these and preserve the field's existing
+      // value, rather than reading buffer[0]. Use find() rather than
+      // operator[] so missing keys are not silently inserted as 0.
+      const auto it = gid_to_buffer_index[e].find(local_gids[i]);
+      permutation.push_back(
+        it != gid_to_buffer_index[e].end() ? it->second : kUnreceivedDof);
+    }
   }
 
   REDEV_ALWAYS_ASSERT(permutation.size() == local_gids.size());
@@ -275,7 +287,8 @@ ExchangePlan GenericFieldExchangePlanner::BuildReceivePlan(
 
 void GenericFieldExchangePlanner::FillGidMessage(
   const FieldLayout& layout, const ExchangePlan& plan,
-  Rank1View<GO, HostMemorySpace> gid_message) const
+  Rank1View<GO, HostMemorySpace> gid_message,
+  const OverlapMask* overlap_mask) const
 {
   PCMS_FUNCTION_TIMER;
   PCMS_ALWAYS_ASSERT(static_cast<size_t>(gid_message.size()) == plan.msg_size);
@@ -286,11 +299,19 @@ void GenericFieldExchangePlanner::FillGidMessage(
   auto offsets = Rank1View<const redev::LO, HostMemorySpace>(
     plan.offsets.data(), plan.offsets.size());
 
+  // Participation must match BuildReversePartitionMap (owned AND in overlap);
+  // otherwise owned-but-non-overlap DOFs inflate the per-rank counts written
+  // here and corrupt the message. A default all-true mask is used when none is
+  // provided so behavior is unchanged for callers without an overlap mask.
+  OverlapMask default_mask(static_cast<size_t>(gids.size()));
+  auto overlap =
+    (overlap_mask ? *overlap_mask : default_mask).GetMask(layout);
+
   std::vector<EntOffsetsArray> per_rank_offsets(plan.dest_ranks.size());
 
   for (LO local_index = 0; local_index < static_cast<LO>(gids.size());
        ++local_index) {
-    if (!owned[local_index])
+    if (!owned[local_index] || !overlap[local_index])
       continue;
 
     LO perm_index = plan.permutation[local_index];

src/pcms/coupler/field_exchange_planner.h

@@ -32,7 +32,8 @@ class FieldExchangePlanner
 
   virtual void FillGidMessage(
     const FieldLayout& layout, const ExchangePlan& plan,
-    Rank1View<GO, HostMemorySpace> gid_message) const = 0;
+    Rank1View<GO, HostMemorySpace> gid_message,
+    const OverlapMask* overlap_mask = nullptr) const = 0;
 
   virtual ~FieldExchangePlanner() noexcept = default;
 };
@@ -51,7 +52,8 @@ class GenericFieldExchangePlanner : public FieldExchangePlanner
 
   void FillGidMessage(
     const FieldLayout& layout, const ExchangePlan& plan,
-    Rank1View<GO, HostMemorySpace> gid_message) const override;
+    Rank1View<GO, HostMemorySpace> gid_message,
+    const OverlapMask* overlap_mask = nullptr) const override;
 };
 
 } // namespace pcms

src/pcms/coupler/field_layout_communicator.cpp

@@ -78,7 +78,8 @@ void FieldLayoutCommunicator::UpdateLayout()
     std::vector<GO> gid_message(plan_.msg_size);
     planner_->FillGidMessage(
       layout_, plan_,
-      Rank1View<GO, HostMemorySpace>(gid_message.data(), gid_message.size()));
+      Rank1View<GO, HostMemorySpace>(gid_message.data(), gid_message.size()),
+      overlap_mask_.get());
 
     channel_.BeginSendCommunicationPhase();
     gid_comm_.Send(gid_message.data());

src/pcms/coupler/field_serializer.h

@@ -34,10 +34,19 @@ class FieldSerializer
     Rank1View<const T, HostMemorySpace> buffer,
     Rank1View<const LO, HostMemorySpace> permutation) const
   {
+    // Seed from the field's current values so DOFs that were not received are
+    // preserved rather than overwritten. This matters for masked coupling: a
+    // permutation entry < 0 (kUnreceivedDof) means the sender did not include
+    // that DOF's GID, so its existing value must be left untouched instead of
+    // reading buffer[0].
+    auto current = field.GetDOFHolderDataHost();
     Kokkos::View<T*, HostMemorySpace> sorted("sorted", permutation.size());
+    for (LO i = 0; i < static_cast<LO>(sorted.size()); ++i) {
+      sorted[i] = current[i];
+    }
     auto owned = layout.GetOwnedHost();
     for (LO i = 0; i < static_cast<LO>(sorted.size()); ++i) {
-      if (owned[i])
+      if (owned[i] && permutation[i] >= 0)
         sorted[i] = buffer[permutation[i]];
     }
     field.SetDOFHolderDataHost(make_const_array_view(sorted));

src/pcms/coupler/overlap_mask.h

@@ -3,6 +3,7 @@
 
 #include "pcms/field/field_layout.h"
 #include "pcms/utility/arrays.h"
+#include "pcms/utility/assert.h"
 #include <Omega_h_array.hpp>
 #include <functional>
 
@@ -27,6 +28,15 @@ struct OverlapMask
     Kokkos::deep_copy(is_overlap_, true);
   }
 
+  // Construct from a precomputed per-DOF-holder host mask (e.g. an MFEM
+  // subdomain selected by element attribute). Indexed by local DOF holder.
+  OverlapMask(size_t size, Kokkos::View<bool*, HostMemorySpace> is_overlap_host)
+    : is_overlap_("overlap_info", size)
+  {
+    PCMS_ALWAYS_ASSERT(is_overlap_host.extent(0) == size);
+    Kokkos::deep_copy(is_overlap_, is_overlap_host);
+  }
+
   // Construct from Omega_h host array
   OverlapMask(size_t size, Omega_h::HostRead<Omega_h::I8> is_overlap_host)
     : is_overlap_("overlap_info", size)

src/pcms/coupler/serializer/xgc.h

@@ -64,7 +64,9 @@ class XGCFieldSerializer : public FieldSerializer<T>
     }
     if (rank_participates_) {
       for (LO i = 0; i < static_cast<LO>(current.size()); ++i) {
-        if (owned[i]) {
+        // permutation[i] < 0 (kUnreceivedDof) => this DOF's GID was not in the
+        // received message; preserve its current value (and avoid buffer[-1]).
+        if (owned[i] && permutation[i] >= 0) {
           full_data[i] = buffer[permutation[i]];
         }
       }

src/pcms/field/CMakeLists.txt

@@ -9,6 +9,7 @@ set(PCMS_FIELD_HEADERS
         layout/empty.h
         out_of_bounds_policy.h
         field_data.h
+        field_factory.h
         function_space.h
         point_evaluator.h
         data/simple.h
@@ -54,6 +55,18 @@ if (PCMS_ENABLE_OMEGA_H)
     )
 endif ()
 
+if (PCMS_ENABLE_MFEM)
+    list(APPEND PCMS_FIELD_SOURCES
+            layout/mfem.cpp
+            data/mfem.cpp
+    )
+    list(APPEND PCMS_FIELD_HEADERS
+            layout/mfem.h
+            data/mfem.h
+            function_space/mfem.h
+    )
+endif ()
+
 if (PCMS_ENABLE_MESHFIELDS)
     list(APPEND PCMS_FIELD_SOURCES
             layout/mesh_fields.cpp
@@ -97,6 +110,10 @@ if(PCMS_ENABLE_OMEGA_H)
     target_link_libraries(pcms_field PRIVATE Kokkos::kokkoskernels)
 endif()
 
+if(PCMS_ENABLE_MFEM)
+    target_link_libraries(pcms_field PUBLIC mfem)
+endif()
+
 if (PCMS_ENABLE_MESHFIELDS)
     target_link_libraries(pcms_field PUBLIC meshfields::meshfields)
 endif ()

src/pcms/field/data/mfem.cpp

@@ -0,0 +1,89 @@
+#include "pcms/field/data/mfem.h"
+
+#include "pcms/utility/arrays.h"
+#include "pcms/utility/assert.h"
+#include "pcms/utility/profile.h"
+
+namespace pcms
+{
+
+MFEMVertexFieldData::MFEMVertexFieldData(mfem::ParFiniteElementSpace& pfes,
+                                         mfem::ParGridFunction& gf,
+                                         FieldMetadata metadata)
+  : pfes_(pfes),
+    gf_(gf),
+    metadata_(metadata),
+    host_data_("mfem_field_data_host", pfes.GetNDofs()),
+    device_data_("mfem_field_data_device", pfes.GetNDofs())
+{
+  PCMS_ALWAYS_ASSERT(pfes_.GetVDim() == 1);
+}
+
+const FieldMetadata& MFEMVertexFieldData::GetMetadata() const
+{
+  return metadata_;
+}
+
+Rank1View<const Real, HostMemorySpace>
+MFEMVertexFieldData::GetDOFHolderDataHost() const
+{
+  PCMS_FUNCTION_TIMER;
+
+  // Pull owner values into a true-DOF vector via the parallel restriction so
+  // every shared vertex reports its owner's value.
+  mfem::Vector true_values(pfes_.GetTrueVSize());
+  gf_.GetTrueDofs(true_values);
+
+  mfem::Array<int> vdofs;
+  const int nv = pfes_.GetNDofs();
+  for (int v = 0; v < nv; ++v) {
+    pfes_.GetVertexDofs(v, vdofs);
+    PCMS_ALWAYS_ASSERT(vdofs.Size() == 1);
+    const int lt = pfes_.GetLocalTDofNumber(vdofs[0]);
+    host_data_(v) = (lt >= 0) ? static_cast<Real>(true_values[lt]) : Real{0};
+  }
+
+  return make_const_array_view(host_data_);
+}
+
+void MFEMVertexFieldData::SetDOFHolderDataHost(
+  Rank1View<const Real, HostMemorySpace> data)
+{
+  PCMS_FUNCTION_TIMER;
+  PCMS_ALWAYS_ASSERT(static_cast<int>(data.size()) == pfes_.GetNDofs());
+
+  // Scatter received owner values into a true-DOF vector, then distribute to
+  // all local DOFs (including shared, non-owned vertices) via the parallel
+  // prolongation.
+  mfem::Vector true_values(pfes_.GetTrueVSize());
+  mfem::Array<int> vdofs;
+  const int nv = pfes_.GetNDofs();
+  for (int v = 0; v < nv; ++v) {
+    pfes_.GetVertexDofs(v, vdofs);
+    PCMS_ALWAYS_ASSERT(vdofs.Size() == 1);
+    const int lt = pfes_.GetLocalTDofNumber(vdofs[0]);
+    if (lt >= 0) {
+      true_values[lt] = static_cast<double>(data[v]);
+    }
+  }
+
+  gf_.SetFromTrueDofs(true_values);
+}
+
+Rank1View<const Real, DeviceMemorySpace>
+MFEMVertexFieldData::GetDOFHolderData() const
+{
+  GetDOFHolderDataHost();
+  Kokkos::deep_copy(device_data_, host_data_);
+  return make_const_array_view(device_data_);
+}
+
+void MFEMVertexFieldData::SetDOFHolderData(
+  Rank1View<const Real, DeviceMemorySpace> data)
+{
+  PCMS_ALWAYS_ASSERT(static_cast<int>(data.size()) == pfes_.GetNDofs());
+  CopyDeviceRank1ViewToHostView(host_data_, data);
+  SetDOFHolderDataHost(make_const_array_view(host_data_));
+}
+
+} // namespace pcms

src/pcms/field/data/mfem.h

@@ -0,0 +1,47 @@
+#ifndef PCMS_FIELD_DATA_MFEM_H
+#define PCMS_FIELD_DATA_MFEM_H
+
+#include "pcms/field/field_data.h"
+#include "pcms/field/field_metadata.h"
+#include "pcms/utility/arrays.h"
+
+#include <mfem.hpp>
+
+namespace pcms
+{
+
+// FieldData backend for an MFEM order-1 H1 (vertex) scalar field. The
+// coefficient store is the live mfem::ParGridFunction: Get/Set operate on it
+// directly so the coupler and the MFEM solver share state.
+//
+// DOF-holder ordering is the local vertex ordering, matching MFEMLayout. To
+// stay consistent across process boundaries, Get/Set round-trip through the
+// FE space true DOFs: Get reads owner values via the parallel restriction and
+// Set distributes received owner values via the parallel prolongation, filling
+// shared (non-owned) vertices on this rank.
+class MFEMVertexFieldData : public FieldData<Real>
+{
+public:
+  MFEMVertexFieldData(mfem::ParFiniteElementSpace& pfes,
+                      mfem::ParGridFunction& gf, FieldMetadata metadata = {});
+
+  const FieldMetadata& GetMetadata() const override;
+
+  Rank1View<const Real, HostMemorySpace> GetDOFHolderDataHost() const override;
+  void SetDOFHolderDataHost(
+    Rank1View<const Real, HostMemorySpace> data) override;
+
+  Rank1View<const Real, DeviceMemorySpace> GetDOFHolderData() const override;
+  void SetDOFHolderData(Rank1View<const Real, DeviceMemorySpace> data) override;
+
+private:
+  mfem::ParFiniteElementSpace& pfes_;
+  mfem::ParGridFunction& gf_;
+  FieldMetadata metadata_;
+  mutable Kokkos::View<Real*, HostMemorySpace> host_data_;
+  mutable Kokkos::View<Real*, DeviceMemorySpace> device_data_;
+};
+
+} // namespace pcms
+
+#endif // PCMS_FIELD_DATA_MFEM_H

src/pcms/field/data/xgc.h

@@ -25,7 +25,7 @@ class XGCFieldData : public FieldData<T>
                        layout_->GetFullDataSize());
   }
 
-  // Self-allocating constructor: XGCFunctionSpace::CreateFieldImpl uses this
+  // Self-allocating constructor: XGCFieldFactory::CreateFieldImpl uses this
   // to produce a field with internally-managed storage.
   XGCFieldData(std::shared_ptr<const XGCFieldLayout> layout,
                FieldMetadata metadata)