diff --git a/palace/fem/bilinearform.cpp b/palace/fem/bilinearform.cpp
index 6da03d964..b6e991848 100644
--- a/palace/fem/bilinearform.cpp
+++ b/palace/fem/bilinearform.cpp
@@ -49,6 +49,7 @@ BilinearForm::PartialAssemble(const FiniteElementSpace &trial_fespace,
   // This should work fine if some threads create an empty operator (no elements or boundary
   // elements).
   const std::size_t nt = ceed::internal::GetCeedObjects().size();
+  MFEM_CONTRACT_VAR(nt);
   PalacePragmaOmp(parallel if (nt > 1))
   {
     Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
@@ -215,6 +216,7 @@ std::unique_ptr<ceed::Operator> DiscreteLinearOperator::PartialAssemble() const
   // This should work fine if some threads create an empty operator (no elements or boundary
   // elements).
   const std::size_t nt = ceed::internal::GetCeedObjects().size();
+  MFEM_CONTRACT_VAR(nt);
   PalacePragmaOmp(parallel if (nt > 1))
   {
     Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
diff --git a/palace/fem/coefficient.hpp b/palace/fem/coefficient.hpp
index e4896ab30..d85b5fc69 100644
--- a/palace/fem/coefficient.hpp
+++ b/palace/fem/coefficient.hpp
@@ -691,6 +691,41 @@ class BdrFieldCoefficient : public mfem::Coefficient, public BdrGridFunctionCoef
   }
 };
 
+// Returns the inner product of a vector grid function and a vector coefficient. Intended
+// to be evaluated on boundary elements, where MFEM's GridFunctionCoefficient would not work
+// even if the fields have the right partial continuity.
+class BdrInnerProductCoefficient : public mfem::Coefficient,
+                                   public BdrGridFunctionCoefficient
+{
+private:
+  const mfem::ParGridFunction &U;
+  mfem::VectorCoefficient &coeff;
+
+public:
+  BdrInnerProductCoefficient(const mfem::ParGridFunction &U, mfem::VectorCoefficient &coeff)
+    : mfem::Coefficient(), BdrGridFunctionCoefficient(*U.ParFESpace()->GetParMesh()), U(U),
+      coeff(coeff)
+  {
+    MFEM_VERIFY(U.VectorDim() == coeff.GetVDim(),
+                "Vector dimension mismatch for BdrInnerProductCoefficient!");
+  }
+
+  double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override
+  {
+    // Get neighboring elements.
+    MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
+                "Unexpected element type in BdrInnerProductCoefficient!");
+    GetBdrElementNeighborTransformations(T.ElementNo, ip);
+
+    // Always evaluate in element 1, assuming the inner product is continuous.
+    double V_data[3], Q_data[3];
+    mfem::Vector V(V_data, T.GetSpaceDim()), Q(Q_data, T.GetSpaceDim());
+    U.GetVectorValue(*FET.Elem1, FET.Elem1->GetIntPoint(), V);
+    coeff.Eval(Q, T, ip);
+    return V * Q;
+  }
+};
+
 //
 // More helpful coefficient types. Wrapper coefficients allow additions of scalar and vector
 // or matrix coefficients. Restricted coefficients only compute the coefficient if for the
diff --git a/palace/fem/integrator.cpp b/palace/fem/integrator.cpp
index 6799329c6..dfb192c10 100644
--- a/palace/fem/integrator.cpp
+++ b/palace/fem/integrator.cpp
@@ -4,6 +4,8 @@
 #include "integrator.hpp"
 
 #include "fem/libceed/integrator.hpp"
+#include "utils/communication.hpp"
+#include "utils/omp.hpp"
 
 namespace palace
 {
@@ -20,11 +22,15 @@ int DefaultIntegrationOrder::Get(const mfem::IsoparametricTransformation &T)
 
 int DefaultIntegrationOrder::Get(const mfem::ElementTransformation &T)
 {
+#if defined(MFEM_DEBUG)
   const auto *T_iso = dynamic_cast<const mfem::IsoparametricTransformation *>(&T);
   MFEM_VERIFY(
       T_iso,
       "Unexpected non-isoparametric element transformation to calculate quadrature order!");
   return Get(*T_iso);
+#else
+  return Get(static_cast<const mfem::IsoparametricTransformation &>(T));
+#endif
 }
 
 int DefaultIntegrationOrder::Get(const mfem::Mesh &mesh, mfem::Geometry::Type geom)
@@ -35,6 +41,68 @@ int DefaultIntegrationOrder::Get(const mfem::Mesh &mesh, mfem::Geometry::Type ge
   return Get(T);
 }
 
+double IntegrateFunction(
+    const mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
+    mfem::Coefficient &Q,
+    std::function<int(const mfem::ElementTransformation &)> GetQuadratureOrder)
+{
+  double sum = IntegrateFunctionLocal(mesh, marker, bdr, Q, GetQuadratureOrder);
+  Mpi::GlobalSum(1, &sum, mesh.GetComm());
+  return sum;
+}
+
+double IntegrateFunctionLocal(
+    const mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
+    mfem::Coefficient &Q,
+    std::function<int(const mfem::ElementTransformation &)> GetQuadratureOrder)
+{
+  auto ElementIntegral = [&Q, &GetQuadratureOrder](mfem::ElementTransformation &T)
+  {
+    double sum = 0.0;
+    const mfem::IntegrationRule &ir =
+        mfem::IntRules.Get(T.GetGeometryType(), GetQuadratureOrder(T));
+    for (int j = 0; j < ir.GetNPoints(); j++)
+    {
+      const mfem::IntegrationPoint &ip = ir.IntPoint(j);
+      T.SetIntPoint(&ip);
+      sum += Q.Eval(T, ip) * ip.weight * T.Weight();
+    }
+    return sum;
+  };
+  double sum = 0.0;
+  PalacePragmaOmp(parallel reduction(+ : sum))
+  {
+    mfem::IsoparametricTransformation T;
+    if (bdr)
+    {
+      PalacePragmaOmp(for schedule(static))
+      for (int i = 0; i < mesh.GetNBE(); i++)
+      {
+        if (!marker[mesh.GetBdrAttribute(i) - 1])
+        {
+          continue;
+        }
+        mesh.GetBdrElementTransformation(i, &T);
+        sum += ElementIntegral(T);
+      }
+    }
+    else
+    {
+      PalacePragmaOmp(for schedule(static))
+      for (int i = 0; i < mesh.GetNE(); i++)
+      {
+        if (!marker[mesh.GetAttribute(i) - 1])
+        {
+          continue;
+        }
+        mesh.GetElementTransformation(i, &T);
+        sum += ElementIntegral(T);
+      }
+    }
+  }
+  return sum;
+}
+
 }  // namespace fem
 
 void DiscreteInterpolator::Assemble(Ceed ceed, CeedElemRestriction trial_restr,
diff --git a/palace/fem/integrator.hpp b/palace/fem/integrator.hpp
index 2f62a2007..31da14fe4 100644
--- a/palace/fem/integrator.hpp
+++ b/palace/fem/integrator.hpp
@@ -4,6 +4,7 @@
 #ifndef PALACE_FEM_INTEGRATOR_HPP
 #define PALACE_FEM_INTEGRATOR_HPP
 
+#include <functional>
 #include <mfem.hpp>
 #include "fem/libceed/ceed.hpp"
 
@@ -33,6 +34,19 @@ struct DefaultIntegrationOrder
   static int Get(const mfem::Mesh &mesh, mfem::Geometry::Type geom);
 };
 
+double IntegrateFunction(
+    const mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
+    mfem::Coefficient &Q,
+    std::function<int(const mfem::ElementTransformation &)> GetQuadratureOrder =
+        [](const mfem::ElementTransformation &T)
+    { return DefaultIntegrationOrder::Get(T); });
+double IntegrateFunctionLocal(
+    const mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
+    mfem::Coefficient &Q,
+    std::function<int(const mfem::ElementTransformation &)> GetQuadratureOrder =
+        [](const mfem::ElementTransformation &T)
+    { return DefaultIntegrationOrder::Get(T); });
+
 }  // namespace fem
 
 // Base class for libCEED-based bilinear form integrators.
@@ -310,7 +324,7 @@ class VectorFEBoundaryLFIntegrator : public mfem::LinearFormIntegrator
   mfem::Vector f_loc, f_hat;
 
 public:
-  VectorFEBoundaryLFIntegrator(mfem::VectorCoefficient &QG) : Q(QG) {}
+  VectorFEBoundaryLFIntegrator(mfem::VectorCoefficient &Q) : Q(Q) {}
 
   void AssembleRHSElementVect(const mfem::FiniteElement &fe, mfem::ElementTransformation &T,
                               mfem::Vector &elvect) override;
@@ -325,7 +339,7 @@ class BoundaryLFIntegrator : public mfem::LinearFormIntegrator
   mfem::Vector shape;
 
 public:
-  BoundaryLFIntegrator(mfem::Coefficient &QG) : Q(QG) {}
+  BoundaryLFIntegrator(mfem::Coefficient &Q) : Q(Q) {}
 
   void AssembleRHSElementVect(const mfem::FiniteElement &fe, mfem::ElementTransformation &T,
                               mfem::Vector &elvect) override;
diff --git a/palace/linalg/errorestimator.cpp b/palace/linalg/errorestimator.cpp
index 479818155..d96bec2ff 100644
--- a/palace/linalg/errorestimator.cpp
+++ b/palace/linalg/errorestimator.cpp
@@ -189,6 +189,7 @@ Vector ComputeErrorEstimates(const VecType &F, VecType &F_gf, VecType &G, VecTyp
   estimates.UseDevice(true);
   estimates = 0.0;
   const std::size_t nt = ceed::internal::GetCeedObjects().size();
+  MFEM_CONTRACT_VAR(nt);
   PalacePragmaOmp(parallel if (nt > 1))
   {
     Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
@@ -269,6 +270,7 @@ GradFluxErrorEstimator<VecType>::GradFluxErrorEstimator(
   // discontinuous flux is ε E = ε ∇V.
   const auto &mesh = nd_fespace.GetMesh();
   const std::size_t nt = ceed::internal::GetCeedObjects().size();
+  MFEM_CONTRACT_VAR(nt);
   PalacePragmaOmp(parallel if (nt > 1))
   {
     Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
@@ -388,6 +390,7 @@ CurlFluxErrorEstimator<VecType>::CurlFluxErrorEstimator(
   // discontinuous flux is μ⁻¹ B ≃ μ⁻¹ ∇ × E.
   const auto &mesh = rt_fespace.GetMesh();
   const std::size_t nt = ceed::internal::GetCeedObjects().size();
+  MFEM_CONTRACT_VAR(nt);
   PalacePragmaOmp(parallel if (nt > 1))
   {
     Ceed ceed = ceed::internal::GetCeedObjects()[utils::GetThreadNum()];
diff --git a/palace/models/lumpedportoperator.cpp b/palace/models/lumpedportoperator.cpp
index fb7ff196b..3c9186ceb 100644
--- a/palace/models/lumpedportoperator.cpp
+++ b/palace/models/lumpedportoperator.cpp
@@ -159,80 +159,16 @@ double LumpedPortData::GetExcitationVoltage() const
   }
 }
 
-void LumpedPortData::InitializeLinearForms(mfem::ParFiniteElementSpace &nd_fespace) const
-{
-  const auto &mesh = *nd_fespace.GetParMesh();
-  mfem::Array<int> attr_marker;
-  if (!s || !v)
-  {
-    mfem::Array<int> attr_list;
-    for (const auto &elem : elems)
-    {
-      attr_list.Append(elem->GetAttrList());
-    }
-    int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
-    mesh::AttrToMarker(bdr_attr_max, attr_list, attr_marker);
-  }
-
-  // The port S-parameter, or the projection of the field onto the port mode, is computed
-  // as: (E x H_inc) ⋅ n = E ⋅ (E_inc / Z_s), integrated over the port surface.
-  if (!s)
-  {
-    SumVectorCoefficient fb(mesh.SpaceDimension());
-    for (const auto &elem : elems)
-    {
-      const double Rs = R * GetToSquare(*elem);
-      const double Hinc = (std::abs(Rs) > 0.0)
-                              ? 1.0 / std::sqrt(Rs * elem->GetGeometryWidth() *
-                                                elem->GetGeometryLength() * elems.size())
-                              : 0.0;
-      fb.AddCoefficient(elem->GetModeCoefficient(Hinc));
-    }
-    s = std::make_unique<mfem::LinearForm>(&nd_fespace);
-    s->AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb), attr_marker);
-    s->UseFastAssembly(false);
-    s->UseDevice(false);
-    s->Assemble();
-    s->UseDevice(true);
-  }
-
-  // The voltage across a port is computed using the electric field solution.
-  // We have:
-  //             V = ∫ E ⋅ l̂ dl = 1/w ∫ E ⋅ l̂ dS  (for rectangular ports)
-  // or,
-  //             V = 1/(2π) ∫ E ⋅ r̂ / r dS        (for coaxial ports).
-  // We compute the surface integral via an inner product between the linear form with the
-  // averaging function as a vector coefficient and the solution expansion coefficients.
-  if (!v)
-  {
-    SumVectorCoefficient fb(mesh.SpaceDimension());
-    for (const auto &elem : elems)
-    {
-      fb.AddCoefficient(
-          elem->GetModeCoefficient(1.0 / (elem->GetGeometryWidth() * elems.size())));
-    }
-    v = std::make_unique<mfem::LinearForm>(&nd_fespace);
-    v->AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fb), attr_marker);
-    v->UseFastAssembly(false);
-    v->UseDevice(false);
-    v->Assemble();
-    v->UseDevice(true);
-  }
-}
-
 std::complex<double> LumpedPortData::GetPower(GridFunction &E, GridFunction &B) const
 {
-  // Compute port power, (E x H) ⋅ n = E ⋅ (-n x H), integrated over the port surface using
-  // the computed E and H = μ⁻¹ B fields, where +n is the direction of propagation (into the
-  // domain). The BdrSurfaceCurrentVectorCoefficient computes -n x H for an outward normal,
-  // so we multiply by -1. The linear form is reconstructed from scratch each time due to
-  // changing H.
+  // Compute port power, (E x H⋆) ⋅ n = E ⋅ (-n x H⋆), integrated over the port surface
+  // using the computed E and H = μ⁻¹ B fields, where +n is the direction of propagation
+  // (into the domain). The BdrSurfaceCurrentVectorCoefficient computes -n x H for an
+  // outward normal, so we multiply by -1.
   MFEM_VERIFY((E.HasImag() && B.HasImag()) || (!E.HasImag() && !B.HasImag()),
               "Mismatch between real- and complex-valued E and B fields in port power "
               "calculation!");
-  const bool has_imag = E.HasImag();
-  auto &nd_fespace = *E.ParFESpace();
-  const auto &mesh = *nd_fespace.GetParMesh();
+  const auto &mesh = *E.ParFESpace()->GetParMesh();
   SumVectorCoefficient fbr(mesh.SpaceDimension()), fbi(mesh.SpaceDimension());
   mfem::Array<int> attr_list;
   for (const auto &elem : elems)
@@ -240,7 +176,7 @@ std::complex<double> LumpedPortData::GetPower(GridFunction &E, GridFunction &B)
     fbr.AddCoefficient(
         std::make_unique<RestrictedVectorCoefficient<BdrSurfaceCurrentVectorCoefficient>>(
             elem->GetAttrList(), B.Real(), mat_op));
-    if (has_imag)
+    if (E.HasImag())
     {
       fbi.AddCoefficient(
           std::make_unique<RestrictedVectorCoefficient<BdrSurfaceCurrentVectorCoefficient>>(
@@ -250,44 +186,50 @@ std::complex<double> LumpedPortData::GetPower(GridFunction &E, GridFunction &B)
   }
   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
-  std::complex<double> dot;
-  {
-    mfem::LinearForm pr(&nd_fespace);
-    pr.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbr), attr_marker);
-    pr.UseFastAssembly(false);
-    pr.UseDevice(false);
-    pr.Assemble();
-    pr.UseDevice(true);
-    dot = -(pr * E.Real()) + (has_imag ? -1i * (pr * E.Imag()) : 0.0);
-  }
-  if (has_imag)
-  {
-    mfem::LinearForm pi(&nd_fespace);
-    pi.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(fbi), attr_marker);
-    pi.UseFastAssembly(false);
-    pi.UseDevice(false);
-    pi.Assemble();
-    pi.UseDevice(true);
-    dot += -(pi * E.Imag()) + 1i * (pi * E.Real());
-    Mpi::GlobalSum(1, &dot, E.ParFESpace()->GetComm());
+  if (E.HasImag())
+  {
+    BdrInnerProductCoefficient prr(E.Real(), fbr), pir(E.Imag(), fbr), pri(E.Real(), fbi),
+        pii(E.Imag(), fbi);
+    mfem::SumCoefficient pr(prr, pii), pi(pir, pri, 1.0, -1.0);
+    std::complex<double> dot = {-fem::IntegrateFunctionLocal(mesh, attr_marker, true, pr),
+                                -fem::IntegrateFunctionLocal(mesh, attr_marker, true, pi)};
+    Mpi::GlobalSum(1, &dot, mesh.GetComm());
     return dot;
   }
   else
   {
-    double rdot = dot.real();
-    Mpi::GlobalSum(1, &rdot, E.ParFESpace()->GetComm());
-    return rdot;
+    BdrInnerProductCoefficient pr(E.Real(), fbr);
+    double dot = -fem::IntegrateFunctionLocal(mesh, attr_marker, true, pr);
+    Mpi::GlobalSum(1, &dot, mesh.GetComm());
+    return dot;
   }
 }
 
 std::complex<double> LumpedPortData::GetSParameter(GridFunction &E) const
 {
-  // Compute port S-parameter, or the projection of the field onto the port mode.
-  InitializeLinearForms(*E.ParFESpace());
-  std::complex<double> dot((*s) * E.Real(), 0.0);
+  // The port S-parameter, or the projection of the field onto the port mode, is computed
+  // as: (E x H_inc⋆) ⋅ n = E ⋅ (E_inc / Z_s), integrated over the port surface.
+  const auto &mesh = *E.ParFESpace()->GetParMesh();
+  SumVectorCoefficient fb(mesh.SpaceDimension());
+  mfem::Array<int> attr_list;
+  for (const auto &elem : elems)
+  {
+    const double Rs = R * GetToSquare(*elem);
+    const double Hinc = (std::abs(Rs) > 0.0)
+                            ? 1.0 / std::sqrt(Rs * elem->GetGeometryWidth() *
+                                              elem->GetGeometryLength() * elems.size())
+                            : 0.0;
+    fb.AddCoefficient(elem->GetModeCoefficient(Hinc));
+    attr_list.Append(elem->GetAttrList());
+  }
+  int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
+  mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
+  BdrInnerProductCoefficient sr(E.Real(), fb);
+  std::complex<double> dot(fem::IntegrateFunctionLocal(mesh, attr_marker, true, sr), 0.0);
   if (E.HasImag())
   {
-    dot.imag((*s) * E.Imag());
+    BdrInnerProductCoefficient si(E.Imag(), fb);
+    dot.imag(fem::IntegrateFunctionLocal(mesh, attr_marker, true, si));
   }
   Mpi::GlobalSum(1, &dot, E.GetComm());
   return dot;
@@ -295,12 +237,30 @@ std::complex<double> LumpedPortData::GetSParameter(GridFunction &E) const
 
 std::complex<double> LumpedPortData::GetVoltage(GridFunction &E) const
 {
-  // Compute the average voltage across the port.
-  InitializeLinearForms(*E.ParFESpace());
-  std::complex<double> dot((*v) * E.Real(), 0.0);
+  // The voltage across a port is computed using the electric field solution.
+  // We have:
+  //             V = ∫ E ⋅ l̂ dl = 1/w ∫ E ⋅ l̂ dS  (for rectangular ports)
+  // or,
+  //             V = 1/(2π) ∫ E ⋅ r̂ / r dS        (for coaxial ports).
+  // We compute the surface integral via an inner product between the linear form with the
+  // averaging function as a vector coefficient and the solution expansion coefficients.
+  const auto &mesh = *E.ParFESpace()->GetParMesh();
+  SumVectorCoefficient fb(mesh.SpaceDimension());
+  mfem::Array<int> attr_list;
+  for (const auto &elem : elems)
+  {
+    fb.AddCoefficient(
+        elem->GetModeCoefficient(1.0 / (elem->GetGeometryWidth() * elems.size())));
+    attr_list.Append(elem->GetAttrList());
+  }
+  int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
+  mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
+  BdrInnerProductCoefficient vr(E.Real(), fb);
+  std::complex<double> dot(fem::IntegrateFunctionLocal(mesh, attr_marker, true, vr), 0.0);
   if (E.HasImag())
   {
-    dot.imag((*v) * E.Imag());
+    BdrInnerProductCoefficient vi(E.Imag(), fb);
+    dot.imag(fem::IntegrateFunctionLocal(mesh, attr_marker, true, vi));
   }
   Mpi::GlobalSum(1, &dot, E.GetComm());
   return dot;
diff --git a/palace/models/lumpedportoperator.hpp b/palace/models/lumpedportoperator.hpp
index 43593a407..86c18abaa 100644
--- a/palace/models/lumpedportoperator.hpp
+++ b/palace/models/lumpedportoperator.hpp
@@ -45,12 +45,6 @@ class LumpedPortData
   int excitation;
   bool active;
 
-private:
-  // Linear forms for postprocessing integrated quantities on the port.
-  mutable std::unique_ptr<mfem::LinearForm> s, v;
-
-  void InitializeLinearForms(mfem::ParFiniteElementSpace &nd_fespace) const;
-
 public:
   LumpedPortData(const config::LumpedPortData &data, const MaterialOperator &mat_op,
                  const mfem::ParMesh &mesh);
diff --git a/palace/models/postoperator.cpp b/palace/models/postoperator.cpp
index 81a65675b..337981835 100644
--- a/palace/models/postoperator.cpp
+++ b/palace/models/postoperator.cpp
@@ -69,7 +69,7 @@ PostOperator<solver_t>::PostOperator(const IoData &iodata, fem_op_t<solver_t> &f
                                       fem_op_.GetRTSpace());
           }
         }())),
-    surf_post_op(iodata, fem_op->GetMaterialOp(), fem_op->GetH1Space()),
+    surf_post_op(iodata, fem_op->GetMaterialOp(), fem_op->GetMesh()),
     interp_op(iodata, fem_op->GetNDSpace())
 {
   // Define primary grid-functions.
@@ -294,7 +294,7 @@ void PostOperator<solver_t>::InitializeParaviewDataCollection(
   }
 
   // Extract energy density field for electric field energy 1/2 Dᴴ E or magnetic field
-  // energy 1/2 Hᴴ B. Also Poynting vector S = E x H⋆.
+  // energy 1/2 Hᴴ B. Also Poynting vector S = Re{E x H⋆}.
   if (U_e)
   {
     paraview->RegisterCoeffField("U_e", U_e.get());
diff --git a/palace/models/surfacepostoperator.cpp b/palace/models/surfacepostoperator.cpp
index 753c1a04f..21e04560c 100644
--- a/palace/models/surfacepostoperator.cpp
+++ b/palace/models/surfacepostoperator.cpp
@@ -110,17 +110,14 @@ SurfacePostOperator::SurfaceFluxData::GetCoefficient(const mfem::ParGridFunction
   switch (type)
   {
     case SurfaceFluxType::ELECTRIC:
-      return std::make_unique<
-          RestrictedCoefficient<BdrSurfaceFluxCoefficient<SurfaceFluxType::ELECTRIC>>>(
-          attr_list, E, nullptr, mat_op, two_sided, center);
+      return std::make_unique<BdrSurfaceFluxCoefficient<SurfaceFluxType::ELECTRIC>>(
+          E, nullptr, mat_op, two_sided, center);
     case SurfaceFluxType::MAGNETIC:
-      return std::make_unique<
-          RestrictedCoefficient<BdrSurfaceFluxCoefficient<SurfaceFluxType::MAGNETIC>>>(
-          attr_list, nullptr, B, mat_op, two_sided, center);
+      return std::make_unique<BdrSurfaceFluxCoefficient<SurfaceFluxType::MAGNETIC>>(
+          nullptr, B, mat_op, two_sided, center);
     case SurfaceFluxType::POWER:
-      return std::make_unique<
-          RestrictedCoefficient<BdrSurfaceFluxCoefficient<SurfaceFluxType::POWER>>>(
-          attr_list, E, B, mat_op, two_sided, center);
+      return std::make_unique<BdrSurfaceFluxCoefficient<SurfaceFluxType::POWER>>(
+          E, B, mat_op, two_sided, center);
   }
   return {};
 }
@@ -165,32 +162,28 @@ SurfacePostOperator::InterfaceDielectricData::GetCoefficient(
   switch (type)
   {
     case InterfaceDielectricType::DEFAULT:
-      return std::make_unique<RestrictedCoefficient<
-          InterfaceDielectricCoefficient<InterfaceDielectricType::DEFAULT>>>(
-          attr_list, E, mat_op, t, epsilon);
+      return std::make_unique<
+          InterfaceDielectricCoefficient<InterfaceDielectricType::DEFAULT>>(E, mat_op, t,
+                                                                            epsilon);
     case InterfaceDielectricType::MA:
-      return std::make_unique<RestrictedCoefficient<
-          InterfaceDielectricCoefficient<InterfaceDielectricType::MA>>>(attr_list, E,
-                                                                        mat_op, t, epsilon);
+      return std::make_unique<InterfaceDielectricCoefficient<InterfaceDielectricType::MA>>(
+          E, mat_op, t, epsilon);
     case InterfaceDielectricType::MS:
-      return std::make_unique<RestrictedCoefficient<
-          InterfaceDielectricCoefficient<InterfaceDielectricType::MS>>>(attr_list, E,
-                                                                        mat_op, t, epsilon);
+      return std::make_unique<InterfaceDielectricCoefficient<InterfaceDielectricType::MS>>(
+          E, mat_op, t, epsilon);
     case InterfaceDielectricType::SA:
-      return std::make_unique<RestrictedCoefficient<
-          InterfaceDielectricCoefficient<InterfaceDielectricType::SA>>>(attr_list, E,
-                                                                        mat_op, t, epsilon);
+      return std::make_unique<InterfaceDielectricCoefficient<InterfaceDielectricType::SA>>(
+          E, mat_op, t, epsilon);
   }
   return {};  // For compiler warning
 }
 
 SurfacePostOperator::SurfacePostOperator(const IoData &iodata,
                                          const MaterialOperator &mat_op,
-                                         mfem::ParFiniteElementSpace &h1_fespace)
-  : mat_op(mat_op), h1_fespace(h1_fespace)
+                                         const mfem::ParMesh &mesh)
+  : mat_op(mat_op)
 {
   // Check that boundary attributes have been specified correctly.
-  const auto &mesh = *h1_fespace.GetParMesh();
   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   mfem::Array<int> bdr_attr_marker;
   if (!iodata.boundaries.postpro.flux.empty() ||
@@ -215,7 +208,7 @@ SurfacePostOperator::SurfacePostOperator(const IoData &iodata,
                     data.type == config::SurfaceFluxData::Type::MAGNETIC,
                 "Electric field or power surface flux postprocessing are not available "
                 "for magnetostatic problems!");
-    flux_surfs.try_emplace(idx, data, *h1_fespace.GetParMesh(), bdr_attr_marker);
+    flux_surfs.try_emplace(idx, data, mesh, bdr_attr_marker);
   }
 
   // Interface dielectric postprocessing.
@@ -225,7 +218,7 @@ SurfacePostOperator::SurfacePostOperator(const IoData &iodata,
               "magnetostatic problems!");
   for (const auto &[idx, data] : iodata.boundaries.postpro.dielectric)
   {
-    eps_surfs.try_emplace(idx, data, *h1_fespace.GetParMesh(), bdr_attr_marker);
+    eps_surfs.try_emplace(idx, data, mesh, bdr_attr_marker);
   }
 }
 
@@ -239,17 +232,17 @@ std::complex<double> SurfacePostOperator::GetSurfaceFlux(int idx, const GridFunc
   MFEM_VERIFY(it != flux_surfs.end(),
               "Unknown surface flux postprocessing index requested!");
   const bool has_imag = (E) ? E->HasImag() : B->HasImag();
-  const auto &mesh = *h1_fespace.GetParMesh();
+  const auto &mesh = (E) ? *E->ParFESpace()->GetParMesh() : *B->ParFESpace()->GetParMesh();
   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, it->second.attr_list);
   auto f =
       it->second.GetCoefficient(E ? &E->Real() : nullptr, B ? &B->Real() : nullptr, mat_op);
-  std::complex<double> dot(GetLocalSurfaceIntegral(*f, attr_marker), 0.0);
+  std::complex<double> dot(fem::IntegrateFunctionLocal(mesh, attr_marker, true, *f), 0.0);
   if (has_imag)
   {
     f = it->second.GetCoefficient(E ? &E->Imag() : nullptr, B ? &B->Imag() : nullptr,
                                   mat_op);
-    double doti = GetLocalSurfaceIntegral(*f, attr_marker);
+    double doti = fem::IntegrateFunctionLocal(mesh, attr_marker, true, *f);
     if (it->second.type == SurfaceFluxType::POWER)
     {
       dot += doti;
@@ -259,7 +252,7 @@ std::complex<double> SurfacePostOperator::GetSurfaceFlux(int idx, const GridFunc
       dot.imag(doti);
     }
   }
-  Mpi::GlobalSum(1, &dot, (E) ? E->GetComm() : B->GetComm());
+  Mpi::GlobalSum(1, &dot, mesh.GetComm());
   return dot;
 }
 
@@ -277,28 +270,13 @@ double SurfacePostOperator::GetInterfaceElectricFieldEnergy(int idx,
   auto it = eps_surfs.find(idx);
   MFEM_VERIFY(it != eps_surfs.end(),
               "Unknown interface dielectric postprocessing index requested!");
-  const auto &mesh = *h1_fespace.GetParMesh();
+  const auto &mesh = *E.ParFESpace()->GetParMesh();
   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, it->second.attr_list);
   auto f = it->second.GetCoefficient(E, mat_op);
-  double dot = GetLocalSurfaceIntegral(*f, attr_marker);
-  Mpi::GlobalSum(1, &dot, E.GetComm());
+  double dot = fem::IntegrateFunctionLocal(mesh, attr_marker, true, *f);
+  Mpi::GlobalSum(1, &dot, mesh.GetComm());
   return dot;
 }
 
-double
-SurfacePostOperator::GetLocalSurfaceIntegral(mfem::Coefficient &f,
-                                             const mfem::Array<int> &attr_marker) const
-{
-  // Integrate the coefficient over the boundary attributes making up this surface index.
-  mfem::LinearForm s(&h1_fespace);
-  s.AddBoundaryIntegrator(new BoundaryLFIntegrator(f),
-                          const_cast<mfem::Array<int> &>(attr_marker));
-  s.UseFastAssembly(false);
-  s.UseDevice(false);
-  s.Assemble();
-  s.UseDevice(true);
-  return linalg::LocalSum(s);
-}
-
 }  // namespace palace
diff --git a/palace/models/surfacepostoperator.hpp b/palace/models/surfacepostoperator.hpp
index f514bde8f..eb3fd58dc 100644
--- a/palace/models/surfacepostoperator.hpp
+++ b/palace/models/surfacepostoperator.hpp
@@ -68,20 +68,13 @@ class SurfacePostOperator
   // Reference to material property operator (not owned).
   const MaterialOperator &mat_op;
 
-  // Reference to scalar finite element space used for computing surface integrals (not
-  // owned).
-  mfem::ParFiniteElementSpace &h1_fespace;
-
-  double GetLocalSurfaceIntegral(mfem::Coefficient &f,
-                                 const mfem::Array<int> &attr_marker) const;
-
 public:
   // Data structures for postprocessing the surface with the given type.
   std::map<int, SurfaceFluxData> flux_surfs;
   std::map<int, InterfaceDielectricData> eps_surfs;
 
   SurfacePostOperator(const IoData &iodata, const MaterialOperator &mat_op,
-                      mfem::ParFiniteElementSpace &h1_fespace);
+                      const mfem::ParMesh &mesh);
 
   // Get surface integrals computing electric or magnetic field flux through a boundary.
   std::complex<double> GetSurfaceFlux(int idx, const GridFunction *E,
diff --git a/palace/models/timeoperator.cpp b/palace/models/timeoperator.cpp
index 4871cdddd..0d3678819 100644
--- a/palace/models/timeoperator.cpp
+++ b/palace/models/timeoperator.cpp
@@ -44,7 +44,7 @@ class TimeDependentFirstOrderOperator : public mfem::TimeDependentOperator
 
   // Bindings to SpaceOperator functions to get the system matrix and preconditioner, and
   // construct the linear solver.
-  std::function<void(double dt)> ConfigureLinearSolver;
+  std::function<void(double)> ConfigureLinearSolver;
 
 public:
   TimeDependentFirstOrderOperator(const IoData &iodata, SpaceOperator &space_op,
diff --git a/palace/models/waveportoperator.cpp b/palace/models/waveportoperator.cpp
index 53d934364..4f033605c 100644
--- a/palace/models/waveportoperator.cpp
+++ b/palace/models/waveportoperator.cpp
@@ -283,24 +283,33 @@ ComplexHypreParMatrix GetSystemMatrixB(const mfem::HypreParMatrix *Bttr,
 }
 
 void Normalize(const GridFunction &S0t, GridFunction &E0t, GridFunction &E0n,
-               mfem::LinearForm &sr, mfem::LinearForm &si)
+               mfem::VectorCoefficient &port_nxH0r_func,
+               mfem::VectorCoefficient &port_nxH0i_func)
 {
   // Normalize grid functions to a chosen polarization direction and unit power, |E x H⋆| ⋅
   // n, integrated over the port surface (+n is the direction of propagation). The n x H
   // coefficients are updated implicitly as the only store references to the Et, En grid
   // functions. We choose a (rather arbitrary) phase constraint to at least make results for
   // the same port consistent between frequencies/meshes.
-
-  // |E x H⋆| ⋅ n = |E ⋅ (-n x H⋆)|. This also updates the n x H coefficients depending on
-  // Et, En. Update linear forms for postprocessing too.
+  const auto &mesh = *S0t.ParFESpace()->GetParMesh();
+  mfem::Array<int> attr_marker(mesh.attributes.Max());
+  attr_marker = 1;
+  mfem::VectorGridFunctionCoefficient S0t_func(&S0t.Real()), E0tr_func(&E0t.Real()),
+      E0ti_func(&E0t.Imag());
+  mfem::InnerProductCoefficient sr(S0t_func, port_nxH0r_func),
+      si(S0t_func, port_nxH0i_func), prr(E0tr_func, port_nxH0r_func),
+      pir(E0ti_func, port_nxH0r_func), pri(E0tr_func, port_nxH0i_func),
+      pii(E0ti_func, port_nxH0i_func);
+  mfem::SumCoefficient pr(prr, pii), pi(pir, pri, 1.0, -1.0);
   std::complex<double> dot[2] = {
-      {sr * S0t.Real(), si * S0t.Real()},
-      {-(sr * E0t.Real()) - (si * E0t.Imag()), -(sr * E0t.Imag()) + (si * E0t.Real())}};
-  Mpi::GlobalSum(2, dot, S0t.ParFESpace()->GetComm());
+      fem::IntegrateFunctionLocal(mesh, attr_marker, false, sr) +
+          1i * fem::IntegrateFunctionLocal(mesh, attr_marker, false, si),
+      -fem::IntegrateFunctionLocal(mesh, attr_marker, false, pr) -
+          1i * fem::IntegrateFunctionLocal(mesh, attr_marker, false, pi)};
+  Mpi::GlobalSum(2, dot, mesh.GetComm());
   auto scale = std::abs(dot[0]) / (dot[0] * std::sqrt(std::abs(dot[1])));
   ComplexVector::AXPBY(scale, E0t.Real(), E0t.Imag(), 0.0, E0t.Real(), E0t.Imag());
   ComplexVector::AXPBY(scale, E0n.Real(), E0n.Imag(), 0.0, E0n.Real(), E0n.Imag());
-  ComplexVector::AXPBY(scale, sr, si, 0.0, sr, si);
 
   // This parallel communication is not required since wave port boundaries are true one-
   // sided boundaries.
@@ -912,32 +921,15 @@ void WavePortData::Initialize(double omega)
     port_E0n->Imag().SetFromTrueDofs(e0ni);
   }
 
-  // Configure the linear forms for computing S-parameters (projection of the field onto the
-  // port mode). Normalize the mode for a chosen polarization direction and unit power,
-  // |E x H⋆| ⋅ n, integrated over the port surface (+n is the direction of propagation).
+  // Normalize the mode for a chosen polarization direction and unit power, |E x H⋆| ⋅ n,
+  // integrated over the port surface (+n is the direction of propagation).
   {
     const auto &port_submesh = static_cast<const mfem::ParSubMesh &>(port_mesh->Get());
     BdrSubmeshHVectorCoefficient<ValueType::REAL> port_nxH0r_func(
         *port_E0t, *port_E0n, mat_op, port_submesh, submesh_parent_elems, kn0, omega0);
     BdrSubmeshHVectorCoefficient<ValueType::IMAG> port_nxH0i_func(
         *port_E0t, *port_E0n, mat_op, port_submesh, submesh_parent_elems, kn0, omega0);
-    {
-      port_sr = std::make_unique<mfem::LinearForm>(&port_nd_fespace->Get());
-      port_sr->AddDomainIntegrator(new VectorFEDomainLFIntegrator(port_nxH0r_func));
-      port_sr->UseFastAssembly(false);
-      port_sr->UseDevice(false);
-      port_sr->Assemble();
-      port_sr->UseDevice(true);
-    }
-    {
-      port_si = std::make_unique<mfem::LinearForm>(&port_nd_fespace->Get());
-      port_si->AddDomainIntegrator(new VectorFEDomainLFIntegrator(port_nxH0i_func));
-      port_si->UseFastAssembly(false);
-      port_si->UseDevice(false);
-      port_si->Assemble();
-      port_si->UseDevice(true);
-    }
-    Normalize(*port_S0t, *port_E0t, *port_E0n, *port_sr, *port_si);
+    Normalize(*port_S0t, *port_E0t, *port_E0n, port_nxH0r_func, port_nxH0i_func);
   }
 }
 
@@ -986,40 +978,25 @@ double WavePortData::GetExcitationPower() const
 
 std::complex<double> WavePortData::GetPower(GridFunction &E, GridFunction &B) const
 {
-  // Compute port power, (E x H) ⋅ n = E ⋅ (-n x H), integrated over the port surface using
-  // the computed E and H = μ⁻¹ B fields, where +n is the direction of propagation (into the
-  // domain). The BdrSurfaceCurrentVectorCoefficient computes -n x H for an outward normal,
-  // so we multiply by -1. The linear form is reconstructed from scratch each time due to
-  // changing H.
+  // Compute port power, (E x H⋆) ⋅ n = E ⋅ (-n x H⋆), integrated over the port surface
+  // using the computed E and H = μ⁻¹ B fields, where +n is the direction of propagation
+  // (into the domain). The BdrSurfaceCurrentVectorCoefficient computes -n x H for an
+  // outward normal, so we multiply by -1. This happens on the parent mesh since we don't
+  // store a transfer operator for the B-field.
   MFEM_VERIFY(E.HasImag() && B.HasImag(),
               "Wave ports expect complex-valued E and B fields in port power "
               "calculation!");
-  auto &nd_fespace = *E.ParFESpace();
-  const auto &mesh = *nd_fespace.GetParMesh();
-  BdrSurfaceCurrentVectorCoefficient nxHr_func(B.Real(), mat_op);
-  BdrSurfaceCurrentVectorCoefficient nxHi_func(B.Imag(), mat_op);
+  const auto &mesh = *E.ParFESpace()->GetParMesh();
   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
-  std::complex<double> dot;
-  {
-    mfem::LinearForm pr(&nd_fespace);
-    pr.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(nxHr_func), attr_marker);
-    pr.UseFastAssembly(false);
-    pr.UseDevice(false);
-    pr.Assemble();
-    pr.UseDevice(true);
-    dot = -(pr * E.Real()) - 1i * (pr * E.Imag());
-  }
-  {
-    mfem::LinearForm pi(&nd_fespace);
-    pi.AddBoundaryIntegrator(new VectorFEBoundaryLFIntegrator(nxHi_func), attr_marker);
-    pi.UseFastAssembly(false);
-    pi.UseDevice(false);
-    pi.Assemble();
-    pi.UseDevice(true);
-    dot += -(pi * E.Imag()) + 1i * (pi * E.Real());
-  }
-  Mpi::GlobalSum(1, &dot, nd_fespace.GetComm());
+  BdrSurfaceCurrentVectorCoefficient nxHr_func(B.Real(), mat_op),
+      nxHi_func(B.Imag(), mat_op);
+  BdrInnerProductCoefficient prr(E.Real(), nxHr_func), pir(E.Imag(), nxHr_func),
+      pri(E.Real(), nxHi_func), pii(E.Imag(), nxHi_func);
+  mfem::SumCoefficient pr(prr, pii), pi(pir, pri, 1.0, -1.0);
+  std::complex<double> dot = {-fem::IntegrateFunctionLocal(mesh, attr_marker, true, pr),
+                              -fem::IntegrateFunctionLocal(mesh, attr_marker, true, pi)};
+  Mpi::GlobalSum(1, &dot, mesh.GetComm());
   return dot;
 }
 
@@ -1032,9 +1009,22 @@ std::complex<double> WavePortData::GetSParameter(GridFunction &E) const
               "calculation!");
   port_nd_transfer->Transfer(E.Real(), port_E->Real());
   port_nd_transfer->Transfer(E.Imag(), port_E->Imag());
-  std::complex<double> dot(-((*port_sr) * port_E->Real()) - ((*port_si) * port_E->Imag()),
-                           -((*port_sr) * port_E->Imag()) + ((*port_si) * port_E->Real()));
-  Mpi::GlobalSum(1, &dot, port_nd_fespace->GetComm());
+  const auto &port_submesh = static_cast<const mfem::ParSubMesh &>(port_mesh->Get());
+  mfem::Array<int> attr_marker(port_submesh.attributes.Max());
+  attr_marker = 1;
+  BdrSubmeshHVectorCoefficient<ValueType::REAL> port_nxH0r_func(
+      *port_E0t, *port_E0n, mat_op, port_submesh, submesh_parent_elems, kn0, omega0);
+  BdrSubmeshHVectorCoefficient<ValueType::IMAG> port_nxH0i_func(
+      *port_E0t, *port_E0n, mat_op, port_submesh, submesh_parent_elems, kn0, omega0);
+  mfem::VectorGridFunctionCoefficient Er_func(&port_E->Real()), Ei_func(&port_E->Imag());
+  mfem::InnerProductCoefficient srr(Er_func, port_nxH0r_func),
+      sir(Ei_func, port_nxH0r_func), sri(Er_func, port_nxH0i_func),
+      sii(Ei_func, port_nxH0i_func);
+  mfem::SumCoefficient sr(srr, sii), si(sir, sri, 1.0, -1.0);
+  std::complex<double> dot = {
+      -fem::IntegrateFunctionLocal(port_submesh, attr_marker, false, sr),
+      -fem::IntegrateFunctionLocal(port_submesh, attr_marker, false, si)};
+  Mpi::GlobalSum(1, &dot, port_submesh.GetComm());
   return dot;
 }
 
diff --git a/palace/models/waveportoperator.hpp b/palace/models/waveportoperator.hpp
index 4e01c9e70..06bd4baf2 100644
--- a/palace/models/waveportoperator.hpp
+++ b/palace/models/waveportoperator.hpp
@@ -80,7 +80,6 @@ class WavePortData
   // objects for computing functions of the port modes for use as an excitation or in
   // postprocessing.
   std::unique_ptr<GridFunction> port_E0t, port_E0n, port_S0t, port_E;
-  std::unique_ptr<mfem::LinearForm> port_sr, port_si;
 
 public:
   WavePortData(const config::WavePortData &data, const config::SolverData &solver,
diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp
index 66488fba5..e4645ba68 100644
--- a/palace/utils/geodata.cpp
+++ b/palace/utils/geodata.cpp
@@ -16,6 +16,7 @@
 #include <utility>
 #include <Eigen/Dense>
 #include <fmt/ranges.h>
+#include "fem/integrator.hpp"
 #include "fem/interpolator.hpp"
 #include "utils/communication.hpp"
 #include "utils/diagnostic.hpp"
@@ -1559,56 +1560,18 @@ mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array<int>
 
 double GetSurfaceArea(const mfem::ParMesh &mesh, const mfem::Array<int> &marker)
 {
-  double area = 0.0;
-  PalacePragmaOmp(parallel reduction(+ : area))
-  {
-    mfem::IsoparametricTransformation T;
-    PalacePragmaOmp(for schedule(static))
-    for (int i = 0; i < mesh.GetNBE(); i++)
-    {
-      if (!marker[mesh.GetBdrAttribute(i) - 1])
-      {
-        continue;
-      }
-      mesh.GetBdrElementTransformation(i, &T);
-      const mfem::IntegrationRule &ir = mfem::IntRules.Get(T.GetGeometryType(), T.OrderJ());
-      for (int j = 0; j < ir.GetNPoints(); j++)
-      {
-        const mfem::IntegrationPoint &ip = ir.IntPoint(j);
-        T.SetIntPoint(&ip);
-        area += ip.weight * T.Weight();
-      }
-    }
-  }
-  Mpi::GlobalSum(1, &area, mesh.GetComm());
-  return area;
+  mfem::ConstantCoefficient one(1.0);
+  return fem::IntegrateFunction(mesh, marker, true, one,
+                                [](const mfem::ElementTransformation &T)
+                                { return T.OrderJ(); });
 }
 
 double GetVolume(const mfem::ParMesh &mesh, const mfem::Array<int> &marker)
 {
-  double volume = 0.0;
-  PalacePragmaOmp(parallel reduction(+ : volume))
-  {
-    mfem::IsoparametricTransformation T;
-    PalacePragmaOmp(for schedule(static))
-    for (int i = 0; i < mesh.GetNE(); i++)
-    {
-      if (!marker[mesh.GetAttribute(i) - 1])
-      {
-        continue;
-      }
-      mesh.GetElementTransformation(i, &T);
-      const mfem::IntegrationRule &ir = mfem::IntRules.Get(T.GetGeometryType(), T.OrderJ());
-      for (int j = 0; j < ir.GetNPoints(); j++)
-      {
-        const mfem::IntegrationPoint &ip = ir.IntPoint(j);
-        T.SetIntPoint(&ip);
-        volume += ip.weight * T.Weight();
-      }
-    }
-  }
-  Mpi::GlobalSum(1, &volume, mesh.GetComm());
-  return volume;
+  mfem::ConstantCoefficient one(1.0);
+  return fem::IntegrateFunction(mesh, marker, false, one,
+                                [](const mfem::ElementTransformation &T)
+                                { return T.OrderJ(); });
 }
 
 double RebalanceMesh(const IoData &iodata, std::unique_ptr<mfem::ParMesh> &mesh)