[WIP] Improve integrals of scalar functions for postprocessing

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()];
@@ -214,6 +215,7 @@ std::unique_ptr<ceed::Operator> DiscreteLinearOperator::PartialAssemble() const
   // This should work fine if some threads create an empty operator (no elements or bounday
   // 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()];

palace/fem/coefficient.hpp

@@ -731,6 +731,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

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,

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.
@@ -302,7 +316,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;
@@ -317,7 +331,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;

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()];