Replace inner products with ones vector with sums

palace/fem/lumpedelement.cpp

@@ -23,10 +23,7 @@ double GetArea(mfem::ParFiniteElementSpace &fespace, mfem::Array<int> &attr_mark
   s.UseDevice(false);
   s.Assemble();
   s.UseDevice(true);
-
-  mfem::GridFunction ones(&fespace);
-  ones = 1.0;
-  return linalg::Dot<Vector>(fespace.GetComm(), s, ones);
+  return linalg::Sum<Vector>(fespace.GetComm(), s);
 }
 
 }  // namespace

palace/linalg/vector.cpp

@@ -589,6 +589,18 @@ std::complex<double> LocalDot(const ComplexVector &x, const ComplexVector &y)
   }
 }
 
+double LocalSum(const Vector &x)
+{
+  static hypre::HypreVector X;
+  X.Update(x);
+  return hypre_SeqVectorSumElts(X);
+}
+
+std::complex<double> LocalSum(const ComplexVector &x)
+{
+  return {LocalSum(x.Real()), LocalSum(x.Imag())};
+}
+
 template <>
 void AXPY(double alpha, const Vector &x, Vector &y)
 {

palace/linalg/vector.hpp

@@ -194,6 +194,19 @@ inline auto Normalize(MPI_Comm comm, VecType &x)
   return norm;
 }
 
+// Calculate the local sum of all elements in x.
+double LocalSum(const Vector &x);
+std::complex<double> LocalSum(const ComplexVector &x);
+
+// Calculate the element-wise sum.
+template <typename VecType>
+inline auto Sum(MPI_Comm comm, const VecType &x)
+{
+  auto sum = LocalSum(x);
+  Mpi::GlobalSum(1, &sum, comm);
+  return sum;
+}
+
 // Addition y += alpha * x.
 template <typename VecType, typename ScalarType>
 void AXPY(ScalarType alpha, const VecType &x, VecType &y);

palace/models/surfacepostoperator.cpp

@@ -142,12 +142,8 @@ std::unique_ptr<mfem::Coefficient> SurfacePostOperator::SurfaceFluxData::GetCoef
 SurfacePostOperator::SurfacePostOperator(const IoData &iodata,
                                          const MaterialOperator &mat_op,
                                          mfem::ParFiniteElementSpace &h1_fespace)
-  : mat_op(mat_op), ones(&h1_fespace)
+  : mat_op(mat_op), h1_fespace(h1_fespace)
 {
-  // Define a constant 1 function on the scalar finite element space for computing surface
-  // integrals.
-  ones = 1.0;
-
   // Surface dielectric loss postprocessing.
   for (const auto &[idx, data] : iodata.boundaries.postpro.dielectric)
   {
@@ -261,13 +257,13 @@ double SurfacePostOperator::GetLocalSurfaceIntegral(const SurfaceData &data,
   }
   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);
-  mfem::LinearForm s(ones.FESpace());
+  mfem::LinearForm s(&h1_fespace);
   s.AddBoundaryIntegrator(new BoundaryLFIntegrator(fb), attr_marker);
   s.UseFastAssembly(false);
   s.UseDevice(false);
   s.Assemble();
   s.UseDevice(true);
-  return linalg::LocalDot(s, ones);
+  return linalg::LocalSum(s);
 }
 
 }  // namespace palace

palace/models/surfacepostoperator.hpp

@@ -81,8 +81,9 @@ class SurfacePostOperator
   // Reference to material property operator (not owned).
   const MaterialOperator &mat_op;
 
-  // Unit function used for computing surface integrals.
-  mutable mfem::GridFunction ones;
+  // Reference to scalar finite element space used for computing surface integrals (not
+  // owned).
+  mfem::ParFiniteElementSpace &h1_fespace;
 
   double GetLocalSurfaceIntegral(const SurfaceData &data,
                                  const mfem::ParGridFunction &U) const;