Adding unit test of GridFunction integrals in 1D

# Conflicts:
#	miniapps/common/CMakeLists.txt
#	miniapps/common/mesh_extras.cpp
#	miniapps/common/mesh_extras.hpp
#	miniapps/tools/makefile

fem/gridfunc.cpp

@@ -3653,6 +3653,278 @@ void GridFunction::ComputeElementLpErrors(const double p,
    }
 }
 
+double GridFunction::ComputeIntegral(Coefficient &weight,
+                                     int ir_order)
+{
+   double val = 0.0;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   Vector vals;
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      fe = fes->GetFE(i);
+      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
+      T = fes->GetElementTransformation(i);
+      GetValues(*T, *ir, vals);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+
+         val += ip.weight * T->Weight() * vals[j] * weight.Eval(*T, ip);
+      }
+   }
+
+   return val;
+}
+
+double GridFunction::ComputeIntegral(const Array<int> &attr_marker,
+                                     Coefficient &weight,
+                                     int ir_order)
+{
+   double val = 0.0;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   Vector vals;
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      int attr = fes->GetAttribute(i);
+      if (attr_marker[attr-1] == 0) { continue; }
+
+      fe = fes->GetFE(i);
+      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
+      T = fes->GetElementTransformation(i);
+      GetValues(*T, *ir, vals);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+
+         val += ip.weight * T->Weight() * vals[j] * weight.Eval(*T, ip);
+      }
+   }
+
+   return val;
+}
+
+double GridFunction::ComputeIntegral(VectorCoefficient &v_weight,
+                                     int ir_order)
+{
+   double val = 0.0;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   DenseMatrix vals, v_weights;
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      fe = fes->GetFE(i);
+      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
+      T = fes->GetElementTransformation(i);
+      GetVectorValues(*T, *ir, vals);
+
+      v_weight.Eval(v_weights, *T, *ir);
+
+      // column-wise dot product of the vector field values (in vals) and the
+      // vector weights (in v_weight)
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+
+         double val_ip = 0.0;
+         for (int d = 0; d < vals.Height(); d++)
+         {
+            val_ip += vals(d,j)*v_weights(d,j);
+         }
+
+         val += ip.weight * T->Weight() * val_ip;
+      }
+   }
+
+   return val;
+}
+
+double GridFunction::ComputeIntegral(const Array<int> &attr_marker,
+                                     VectorCoefficient &v_weight,
+                                     int ir_order)
+{
+   double val = 0.0;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   DenseMatrix vals, v_weights;
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      int attr = fes->GetAttribute(i);
+      if (attr_marker[attr-1] == 0) { continue; }
+
+      fe = fes->GetFE(i);
+      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
+      T = fes->GetElementTransformation(i);
+      GetVectorValues(*T, *ir, vals);
+
+      v_weight.Eval(v_weights, *T, *ir);
+
+      // column-wise dot product of the vector field values (in vals) and the
+      // vector weights (in v_weight)
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+
+         double val_ip = 0.0;
+         for (int d = 0; d < vals.Height(); d++)
+         {
+            val_ip += vals(d,j)*v_weights(d,j);
+         }
+
+         val += ip.weight * T->Weight() * val_ip;
+      }
+   }
+
+   return val;
+}
+
+double GridFunction::ComputeBdrIntegral(Coefficient &weight,
+                                        int ir_order)
+{
+   double val = 0.0;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   Vector vals;
+
+   for (int i = 0; i < fes->GetNBE(); i++)
+   {
+      fe = fes->GetBE(i);
+      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
+      T = fes->GetBdrElementTransformation(i);
+      GetValues(*T, *ir, vals);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+
+         val += ip.weight * T->Weight() * vals[j] * weight.Eval(*T, ip);
+      }
+   }
+
+   return val;
+}
+
+double GridFunction::ComputeBdrIntegral(const Array<int> &bdr_attr_marker,
+                                        Coefficient &weight,
+                                        int ir_order)
+{
+   double val = 0.0;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   Vector vals;
+
+   for (int i = 0; i < fes->GetNBE(); i++)
+   {
+      int attr = fes->GetBdrAttribute(i);
+      if (bdr_attr_marker[attr-1] == 0) { continue; }
+
+      fe = fes->GetBE(i);
+      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
+      T = fes->GetBdrElementTransformation(i);
+      GetValues(*T, *ir, vals);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+
+         val += ip.weight * T->Weight() * vals[j] * weight.Eval(*T, ip);
+      }
+   }
+
+   return val;
+}
+
+double GridFunction::ComputeBdrIntegral(VectorCoefficient &v_weight,
+                                        int ir_order)
+{
+   double val = 0.0;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   DenseMatrix vals, v_weights;
+
+   for (int i = 0; i < fes->GetNBE(); i++)
+   {
+      fe = fes->GetBE(i);
+      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
+      T = fes->GetBdrElementTransformation(i);
+      GetVectorValues(*T, *ir, vals);
+
+      v_weight.Eval(v_weights, *T, *ir);
+
+      // column-wise dot product of the vector field values (in vals) and the
+      // vector weights (in v_weight)
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+
+         double val_ip = 0.0;
+         for (int d = 0; d < vals.Height(); d++)
+         {
+            val_ip += vals(d,j)*v_weights(d,j);
+         }
+
+         val += ip.weight * T->Weight() * val_ip;
+      }
+   }
+
+   return val;
+}
+
+double GridFunction::ComputeBdrIntegral(const Array<int> &bdr_attr_marker,
+                                        VectorCoefficient &v_weight,
+                                        int ir_order)
+{
+   double val = 0.0;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   DenseMatrix vals, v_weights;
+
+   for (int i = 0; i < fes->GetNBE(); i++)
+   {
+      int attr = fes->GetBdrAttribute(i);
+      if (bdr_attr_marker[attr-1] == 0) { continue; }
+
+      fe = fes->GetBE(i);
+      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
+      T = fes->GetBdrElementTransformation(i);
+      GetVectorValues(*T, *ir, vals);
+
+      v_weight.Eval(v_weights, *T, *ir);
+
+      // column-wise dot product of the vector field values (in vals) and the
+      // vector weights (in v_weight)
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+
+         double val_ip = 0.0;
+         for (int d = 0; d < vals.Height(); d++)
+         {
+            val_ip += vals(d,j)*v_weights(d,j);
+         }
+
+         val += ip.weight * T->Weight() * val_ip;
+      }
+   }
+
+   return val;
+}
+
 GridFunction & GridFunction::operator=(double value)
 {
    Vector::operator=(value);

fem/gridfunc.hpp

@@ -638,6 +638,34 @@ public:
                                        ) const
    { ComputeElementLpErrors(infinity(), exsol, error, NULL, NULL, irs); }
 
+   virtual double ComputeIntegral(Coefficient &weight,
+                                  int ir_order);
+
+   virtual double ComputeIntegral(const Array<int> &attr_marker,
+                                  Coefficient &weight,
+                                  int ir_order);
+
+   virtual double ComputeIntegral(VectorCoefficient &v_weight,
+                                  int ir_order);
+
+   virtual double ComputeIntegral(const Array<int> &attr_marker,
+                                  VectorCoefficient &v_weight,
+                                  int ir_order);
+
+   virtual double ComputeBdrIntegral(Coefficient &weight,
+                                     int ir_order);
+
+   virtual double ComputeBdrIntegral(const Array<int> &bdr_attr_marker,
+                                     Coefficient &weight,
+                                     int ir_order);
+
+   virtual double ComputeBdrIntegral(VectorCoefficient &v_weight,
+                                     int ir_order);
+
+   virtual double ComputeBdrIntegral(const Array<int> &bdr_attr_marker,
+                                     VectorCoefficient &v_weight,
+                                     int ir_order);
+
    virtual void ComputeFlux(BilinearFormIntegrator &blfi,
                             GridFunction &flux,
                             bool wcoef = true, int subdomain = -1);

miniapps/common/makefile

@@ -49,7 +49,7 @@ SEQ_MINIOBJS = mesh_extras.o fem_extras.o
 ifeq ($(MFEM_USE_MPI),NO)
    MINIOBJS = $(SEQ_MINIOBJS)
 else
-   MINIOBJS = $(SEQ_MINIOBJS) pfem_extras.o
+   MINIOBJS = $(SEQ_MINIOBJS) pmesh_extras.o pfem_extras.o
 endif
 
 .SUFFIXES:

miniapps/common/mesh_extras.cpp

@@ -136,6 +136,324 @@ ElementMeshStream::ElementMeshStream(Element::Type e)
 
 }
 
+double ComputeVolume(const Mesh &mesh, int ir_order)
+{
+   double vol = 0.0;
+
+   IsoparametricTransformation T;
+
+   for (int i=0; i<mesh.GetNE(); i++)
+   {
+      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
+      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
+      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T.SetIntPoint(&ip);
+
+         double w = T.Weight() * ip.weight;
+
+         vol += w;
+      }
+   }
+   return vol;
+}
+
+double ComputeVolume(const Mesh &mesh, const Array<int> &attr_marker,
+                     int ir_order)
+{
+   double vol = 0.0;
+
+   IsoparametricTransformation T;
+
+   for (int i=0; i<mesh.GetNE(); i++)
+   {
+      int attr = mesh.GetAttribute(i);
+      if (attr_marker[attr-1] == 0) { continue; }
+
+      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
+      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
+      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T.SetIntPoint(&ip);
+
+         double w = T.Weight() * ip.weight;
+
+         vol += w;
+      }
+   }
+   return vol;
+}
+
+double ComputeSurfaceArea(const Mesh &mesh, int ir_order)
+{
+   double area = 0.0;
+
+   IsoparametricTransformation T;
+
+   for (int i=0; i<mesh.GetNBE(); i++)
+   {
+      // ElementTransformation *T = mesh.GetBdrElementTransformation(i);
+      const_cast<Mesh&>(mesh).GetBdrElementTransformation(i, &T);
+      Geometry::Type geom = mesh.GetBdrElementBaseGeometry(i);
+      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T.SetIntPoint(&ip);
+
+         double w = T.Weight() * ip.weight;
+
+         area += w;
+      }
+   }
+   return area;
+}
+
+double ComputeSurfaceArea(const Mesh &mesh, const Array<int> &bdr_attr_marker,
+                          int ir_order)
+{
+   double area = 0.0;
+
+   IsoparametricTransformation T;
+
+   for (int i=0; i<mesh.GetNBE(); i++)
+   {
+      int attr = mesh.GetBdrAttribute(i);
+      if (bdr_attr_marker[attr-1] == 0) { continue; }
+
+      const_cast<Mesh&>(mesh).GetBdrElementTransformation(i, &T);
+      Geometry::Type geom = mesh.GetBdrElementBaseGeometry(i);
+      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T.SetIntPoint(&ip);
+
+         double w = T.Weight() * ip.weight;
+
+         area += w;
+      }
+   }
+   return area;
+}
+
+double ComputeZerothMoment(const Mesh &mesh, Coefficient &rho,
+                           int ir_order)
+{
+   double mom = 0.0;
+
+   IsoparametricTransformation T;
+
+   for (int i=0; i<mesh.GetNE(); i++)
+   {
+      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
+      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
+      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T.SetIntPoint(&ip);
+
+         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
+
+         mom += w;
+      }
+   }
+
+   return mom;
+}
+
+void ComputeElementZerothMoments(const Mesh &mesh, Coefficient &rho,
+                                 int ir_order, GridFunction &m)
+{
+   MFEM_ASSERT(m.Size() == mesh.GetNE(), "Invalid GridFunction.  "
+               "Must have a length equal to the number of mesh elements.");
+
+   FiniteElementSpace * fes = m.FESpace();
+
+   Array<int> vdofs;
+
+   IsoparametricTransformation T;
+
+   for (int i=0; i<mesh.GetNE(); i++)
+   {
+      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
+      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
+      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
+
+      fes->GetElementVDofs(i, vdofs);
+
+      double mom = 0.0;
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T.SetIntPoint(&ip);
+
+         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
+
+         mom += w;
+      }
+
+      m[vdofs[0]] = mom;
+   }
+}
+
+double ComputeFirstMoment(const Mesh &mesh, Coefficient &rho,
+                          int ir_order, Vector &mom)
+{
+   double mom0 = 0.0;
+
+   int sdim = mesh.SpaceDimension();
+
+   mom.SetSize(sdim);
+   mom = 0.0;
+
+   double x_data[3];
+   Vector x(x_data, sdim);
+
+   IsoparametricTransformation T;
+
+   for (int i=0; i<mesh.GetNE(); i++)
+   {
+      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
+      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
+      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T.SetIntPoint(&ip);
+         T.Transform(ip, x);
+
+         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
+
+         mom0 += w;
+         mom.Add(w, x);
+      }
+   }
+
+   return mom0;
+}
+
+double ComputeSecondMoment(const Mesh &mesh, Coefficient &rho,
+                           const Vector &center,
+                           int ir_order, DenseMatrix &mom)
+{
+   double mom0 = 0.0;
+
+   int sdim = mesh.SpaceDimension();
+
+   mom.SetSize(sdim);
+   mom = 0.0;
+
+   double x_data[3];
+   Vector x(x_data, sdim);
+
+   IsoparametricTransformation T;
+
+   for (int i=0; i<mesh.GetNE(); i++)
+   {
+      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
+      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
+      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T.SetIntPoint(&ip);
+         T.Transform(ip, x);
+         x.Add(-1.0, center);
+
+         double r2 = x * x;
+
+         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
+
+         mom0 += w;
+
+         for (int k=0; k<sdim; k++)
+         {
+            mom(k,k) += w * r2;
+            for (int l=k; l<sdim; l++)
+            {
+               mom(k,l) -= w * x[k] * x[l];
+            }
+         }
+      }
+   }
+
+   for (int k=0; k<sdim; k++)
+   {
+      for (int l=0; l<k; l++)
+      {
+         mom(k,l) = mom(l,k);
+      }
+   }
+
+   return mom0;
+}
+
+void ComputeElementCentersOfMass(const Mesh &mesh, Coefficient &rho,
+                                 int ir_order, GridFunction &c)
+{
+   int sdim = mesh.SpaceDimension();
+
+   MFEM_ASSERT(c.Size() == mesh.GetNE() * sdim,
+               "Invalid GridFunction.  Must have a length equal to the "
+               "number of mesh elements times the spatial dimension.");
+
+   FiniteElementSpace * fes = c.FESpace();
+
+   Array<int> vdofs;
+
+   double x_data[3];
+   Vector x(x_data, sdim);
+
+   c = 0.0;
+
+   IsoparametricTransformation T;
+
+   for (int i=0; i<mesh.GetNE(); i++)
+   {
+      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
+      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
+      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
+
+      fes->GetElementVDofs(i, vdofs);
+
+      double mom0 = 0.0;
+
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T.SetIntPoint(&ip);
+         T.Transform(ip, x);
+
+         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
+
+         mom0 += w;
+
+         for (int k=0; k<sdim; k++)
+         {
+            c[vdofs[k]] += w * x[k];
+         }
+      }
+      for (int k=0; k<sdim; k++)
+      {
+         c[vdofs[k]] /= mom0;
+      }
+   }
+}
+
 void
 MergeMeshNodes(Mesh * mesh, int logging)
 {

miniapps/common/mesh_extras.hpp

@@ -27,6 +27,42 @@ public:
    ElementMeshStream(Element::Type e);
 };
 
+double ComputeVolume(const Mesh &mesh, int ir_order);
+double ComputeVolume(const Mesh &mesh, const Array<int> &attr_marker,
+                     int ir_order);
+double ComputeSurfaceArea(const Mesh &mesh, int ir_order);
+double ComputeSurfaceArea(const Mesh &mesh, const Array<int> &bdr_attr_marker,
+                          int ir_order);
+
+double ComputeZerothMoment(const Mesh &mesh, Coefficient &rho,
+                           int ir_order);
+double ComputeFirstMoment(const Mesh &mesh, Coefficient &rho,
+                          int ir_order, Vector &mom);
+double ComputeSecondMoment(const Mesh &mesh, Coefficient &rho,
+                           const Vector &center,
+                           int ir_order, DenseMatrix &mom);
+
+inline void ComputeNormalizedFirstMoment(const Mesh &mesh, Coefficient &rho,
+                                         int ir_order, Vector &mom)
+{
+   double mom0 = ComputeFirstMoment(mesh, rho, ir_order, mom);
+   mom /= mom0;
+}
+
+inline void ComputeNormalizedSecondMoment(const Mesh &mesh, Coefficient &rho,
+                                          const Vector &center,
+                                          int ir_order, DenseMatrix &mom)
+{
+   double mom0 = ComputeSecondMoment(mesh, rho, center, ir_order, mom);
+   mom *= 1.0 / mom0;
+}
+
+void ComputeElementZerothMoments(const Mesh &mesh, Coefficient &rho,
+                                 int ir_order, GridFunction &m);
+
+void ComputeElementCentersOfMass(const Mesh &mesh, Coefficient &rho,
+                                 int ir_order, GridFunction &c);
+
 /// Merges vertices which lie at the same location
 void MergeMeshNodes(Mesh * mesh, int logging);
 

miniapps/common/pmesh_extras.cpp

@@ -0,0 +1,109 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "pmesh_extras.hpp"
+#include "mesh_extras.hpp"
+
+#ifdef MFEM_USE_MPI
+
+using namespace std;
+
+namespace mfem
+{
+
+namespace common
+{
+
+double ComputeVolume(const ParMesh &pmesh, const Array<int> &attr_marker,
+                     int ir_order)
+{
+   double loc_vol = ComputeVolume(dynamic_cast<const Mesh&>(pmesh),
+                                  attr_marker, ir_order);
+   double glb_vol = 0.0;
+   MPI_Allreduce(&loc_vol, &glb_vol, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
+   return glb_vol;
+}
+
+double ComputeVolume(const ParMesh &pmesh, int ir_order)
+{
+   double loc_vol = ComputeVolume(dynamic_cast<const Mesh&>(pmesh), ir_order);
+   double glb_vol = 0.0;
+   MPI_Allreduce(&loc_vol, &glb_vol, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
+   return glb_vol;
+}
+
+double ComputeSurfaceArea(const ParMesh &pmesh, int ir_order)
+{
+   double loc_area = ComputeSurfaceArea(dynamic_cast<const Mesh&>(pmesh),
+                                        ir_order);
+   double glb_area = 0.0;
+   MPI_Allreduce(&loc_area, &glb_area, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
+   return glb_area;
+}
+
+double ComputeSurfaceArea(const ParMesh &pmesh,
+                          const Array<int> &bdr_attr_marker, int ir_order)
+{
+   double loc_area = ComputeSurfaceArea(dynamic_cast<const Mesh&>(pmesh),
+                                        bdr_attr_marker, ir_order);
+   double glb_area = 0.0;
+   MPI_Allreduce(&loc_area, &glb_area, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
+   return glb_area;
+}
+
+double ComputeZerothMoment(const ParMesh &pmesh, Coefficient &rho,
+                           int ir_order)
+{
+   double loc_mom = ComputeZerothMoment(dynamic_cast<const Mesh&>(pmesh),
+                                        rho, ir_order);
+   double glb_mom = 0.0;
+   MPI_Allreduce(&loc_mom, &glb_mom, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
+   return glb_mom;
+}
+
+double ComputeFirstMoment(const ParMesh &pmesh, Coefficient &rho,
+                          int ir_order, Vector &mom)
+{
+   int sdim = pmesh.SpaceDimension();
+   mom.SetSize(sdim);
+   double loc_mom_data[3];
+   Vector loc_mom(loc_mom_data, sdim);
+   double loc_mom0 = ComputeFirstMoment(dynamic_cast<const Mesh&>(pmesh),
+                                        rho, ir_order, loc_mom);
+   double glb_mom0 = 0.0;
+   MPI_Allreduce(&loc_mom0, &glb_mom0, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
+   MPI_Allreduce(loc_mom_data, mom.GetData(), sdim, MPI_DOUBLE, MPI_SUM,
+                 pmesh.GetComm());
+   return glb_mom0;
+}
+
+double ComputeSecondMoment(const ParMesh &pmesh, Coefficient &rho,
+                           const Vector &center,
+                           int ir_order, DenseMatrix &mom)
+{
+   int sdim = pmesh.SpaceDimension();
+   mom.SetSize(sdim);
+   double loc_mom_data[9];
+   DenseMatrix loc_mom(loc_mom_data, sdim, sdim);
+   double loc_mom0 = ComputeSecondMoment(dynamic_cast<const Mesh&>(pmesh),
+                                         rho, center, ir_order, loc_mom);
+   double glb_mom0 = 0.0;
+   MPI_Allreduce(&loc_mom0, &glb_mom0, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
+   MPI_Allreduce(loc_mom_data, mom.GetData(), sdim * sdim, MPI_DOUBLE, MPI_SUM,
+                 pmesh.GetComm());
+   return glb_mom0;
+}
+
+} // namespace common
+
+} // namespace mfem
+
+#endif

miniapps/common/pmesh_extras.hpp

@@ -0,0 +1,63 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_PMESH_EXTRAS
+#define MFEM_PMESH_EXTRAS
+
+#include "mfem.hpp"
+
+#ifdef MFEM_USE_MPI
+
+namespace mfem
+{
+
+namespace common
+{
+
+double ComputeVolume(const ParMesh &pmesh, int ir_order);
+double ComputeVolume(const ParMesh &mesh, const Array<int> &attr_marker,
+                     int ir_order);
+double ComputeSurfaceArea(const ParMesh &pmesh, int ir_order);
+double ComputeSurfaceArea(const ParMesh &mesh,
+                          const Array<int> &bdr_attr_marker,
+                          int ir_order);
+
+double ComputeZerothMoment(const ParMesh &pmesh, Coefficient &rho,
+                           int ir_order);
+double ComputeFirstMoment(const ParMesh &pmesh, Coefficient &rho,
+                          int ir_order, Vector &mom);
+double ComputeSecondMoment(const ParMesh &pmesh, Coefficient &rho,
+                           const Vector &center,
+                           int ir_order, DenseMatrix &mom);
+
+inline void ComputeNormalizedFirstMoment(const ParMesh &pmesh,
+                                         Coefficient &rho,
+                                         int ir_order, Vector &mom)
+{
+   double mom0 = ComputeFirstMoment(pmesh, rho, ir_order, mom);
+   mom /= mom0;
+}
+
+inline void ComputeNormalizedSecondMoment(const ParMesh &pmesh,
+                                          Coefficient &rho,
+                                          const Vector &center,
+                                          int ir_order, DenseMatrix &mom)
+{
+   double mom0 = ComputeSecondMoment(pmesh, rho, center, ir_order, mom);
+   mom *= 1.0 / mom0;
+}
+
+} // namespace common
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+#endif

miniapps/tools/makefile

@@ -25,8 +25,10 @@ include $(DEFAULTS_MK)
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer
-PAR_MINIAPPS =
+
+SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer \
+	seq-moment
+PAR_MINIAPPS = par-moment
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -38,6 +40,7 @@ endif
 .PHONY: all lib-common clean clean-build clean-exec
 .PRECIOUS: %.o
 
+#COMMON_O=../common/fem_extras.o ../common/mesh_extras.o ../common/pmesh_extras.o
 COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
@@ -60,6 +63,15 @@ display-basis: %: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
 
+
+seq-moment: %: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
+
+par-moment: %: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
+
 # Rule for building lib-common
 lib-common:
 	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common

miniapps/tools/par-moment.cpp

@@ -0,0 +1,181 @@
+#include "mfem.hpp"
+#include "../common/mesh_extras.hpp"
+#include "../common/pmesh_extras.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::common;
+
+static int prob_ = -1;
+
+double densityFunc(const Vector &x)
+{
+   switch (prob_)
+   {
+      case 0:
+         // Linear in the radius
+         return sqrt(x * x);
+      case 1:
+         // Hydrostatic equilibrium (1D)
+         return pow(cosh(sqrt(0.5 * (x * x))), -2.0);
+      case 2:
+         // Hydrostatic equilibrium (2D)
+         return pow(1.0 + 0.125 * (x * x), -2.0);
+      case 3:
+         // Hydrostatic equilibrium (3D), approx.
+         return pow(1.0 + (x * x) / 9.0, -1.5);
+      case 4:
+         // Off-center Gaussian
+      {
+         double s_data[3];
+         Vector s(s_data, x.Size());
+         s = 0.0; s[0] = -1.0; s += x;
+         return exp(-(s * s));
+      }
+   }
+   // Default to homogeneous
+   return 1.0;
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/ball-nurbs.mesh";
+   int ser_ref_levels = 2;
+   int par_ref_levels = 0;
+   int ir_order = 2;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&ir_order, "-o", "--order",
+                  "Integration rule order.");
+   args.AddOption(&prob_, "-d", "--density",
+                  "Density profile:\n"
+                  "  0 - Linear increase with radius,\n"
+                  "  1 - Hydrostatic equilibrium in 1D,\n"
+                  "  2 - Hydrostatic equilibrium in 2D,\n"
+                  "  3 - Hydrostatic equilibrium in 3D (approximately),\n"
+                  "  4 - Gaussian centered at x = 1,\n"
+                  "  Default - homogeneous.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement where 'ref_levels' is user defined.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   FunctionCoefficient density(densityFunc);
+
+   L2_FECollection fec(0, dim);
+   ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
+   ParFiniteElementSpace *vfes = new ParFiniteElementSpace(pmesh, &fec, sdim);
+   ParGridFunction elemMass(fes);
+   ParGridFunction elemCent(vfes);
+
+   double vol = ComputeVolume(*pmesh, ir_order);
+   double area = ComputeSurfaceArea(*pmesh, ir_order);
+
+   Vector mom1(sdim);
+   Vector cent(sdim); cent = 0.0;
+   DenseMatrix mom2(sdim);
+
+   double mass  = ComputeZerothMoment(*pmesh, density, ir_order);
+   double mass1 = ComputeFirstMoment(*pmesh, density, ir_order, mom1);
+   double mass2 = ComputeSecondMoment(*pmesh, density, cent,
+                                      ir_order, mom2);
+
+   if (myid == 0)
+   {
+      cout << "Volume:       " << vol << endl;
+      cout << "Surface Area: " << area << endl;
+      cout << "Mass:         " << mass << endl;
+      cout << "Mass:         " << mass1 << endl;
+      cout << "Mass:         " << mass2 << endl;
+
+      cout << "First Moment:   "; mom1.Print(cout);
+      mom1 /= mass1;
+      cout << "Center of Mass: "; mom1.Print(cout);
+      cout << "Second Moment (moment of inertia):\n";
+      mom2.Print(cout);
+   }
+
+   ComputeElementZerothMoments(*pmesh, density, ir_order, elemMass);
+   ComputeElementCentersOfMass(*pmesh, density, ir_order, elemCent);
+
+   // 15. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream m_sock(vishost, visport);
+      m_sock << "parallel " << num_procs << " " << myid << "\n";
+      m_sock.precision(8);
+      m_sock << "solution\n" << *pmesh << elemMass
+             << "window_title 'Element Masses'" << flush;
+
+      socketstream c_sock(vishost, visport);
+      c_sock << "parallel " << num_procs << " " << myid << "\n";
+      c_sock.precision(8);
+      c_sock << "solution\n" << *pmesh << elemCent
+             << "window_title 'Element Centers'"
+             << "keys vvv" << flush;
+   }
+
+   // 16. Free the used memory.
+   delete fes;
+   delete vfes;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}

miniapps/tools/seq-moment.cpp

@@ -0,0 +1,144 @@
+#include "mfem.hpp"
+#include "../common/mesh_extras.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::common;
+
+static int prob_ = -1;
+
+double densityFunc(const Vector &x)
+{
+   switch (prob_)
+   {
+      case 0:
+         // Linear in the radius
+         return sqrt(x * x);
+      case 1:
+         // Hydrostatic equilibrium (1D)
+         return pow(cosh(sqrt(0.5 * (x * x))), -2.0);
+      case 2:
+         // Hydrostatic equilibrium (2D)
+         return pow(1.0 + 0.125 * (x * x), -2.0);
+      case 3:
+         // Hydrostatic equilibrium (3D), approx.
+         return pow(1.0 + (x * x) / 9.0, -1.5);
+      case 4:
+         // Off-center Gaussian
+      {
+         double s_data[3];
+         Vector s(s_data, x.Size());
+         s = 0.0; s[0] = -1.0; s += x;
+         return exp(-(s * s));
+      }
+   }
+   // Default to homogeneous
+   return 1.0;
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/ball-nurbs.mesh";
+   int ref_levels = 2;
+   int ir_order = 2;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&ir_order, "-o", "--order",
+                  "Integration rule order.");
+   args.AddOption(&prob_, "-d", "--density",
+                  "Density profile:\n"
+                  "  0 - Linear increase with radius,\n"
+                  "  1 - Hydrostatic equilibrium in 1D,\n"
+                  "  2 - Hydrostatic equilibrium in 2D,\n"
+                  "  3 - Hydrostatic equilibrium in 3D (approximately),\n"
+                  "  4 - Gaussian centered at x = 1,\n"
+                  "  Default - homogeneous.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement where 'ref_levels' is user defined.
+   {
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   FunctionCoefficient density(densityFunc);
+
+   L2_FECollection fec(0, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+   FiniteElementSpace vfes(&mesh, &fec, sdim);
+   GridFunction elemMass(&fes);
+   GridFunction elemCent(&vfes);
+
+   double vol = ComputeVolume(mesh, ir_order);
+   double area = ComputeSurfaceArea(mesh, ir_order);
+
+   Vector mom1(sdim);
+   Vector cent(sdim); cent = 0.0;
+   DenseMatrix mom2(sdim);
+
+   double mass  = ComputeZerothMoment(mesh, density, ir_order);
+   double mass1 = ComputeFirstMoment(mesh, density, ir_order, mom1);
+   double mass2 = ComputeSecondMoment(mesh, density, cent,
+                                      ir_order, mom2);
+
+   cout << "Volume:       " << vol << endl;
+   cout << "Surface Area: " << area << endl;
+   cout << "Mass:         " << mass << endl;
+   cout << "Mass:         " << mass1 << endl;
+   cout << "Mass:         " << mass2 << endl;
+
+   cout << "First Moment:   "; mom1.Print(cout);
+   mom1 /= mass1;
+   cout << "Center of Mass: "; mom1.Print(cout);
+   cout << "Second Moment (moment of inertia):\n";
+   mom2.Print(cout);
+
+   ComputeElementZerothMoments(mesh, density, ir_order, elemMass);
+   ComputeElementCentersOfMass(mesh, density, ir_order, elemCent);
+
+   // 15. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream m_sock(vishost, visport);
+      m_sock.precision(8);
+      m_sock << "solution\n" << mesh << elemMass
+             << "window_title 'Element Masses'" << flush;
+
+      socketstream c_sock(vishost, visport);
+      c_sock.precision(8);
+      c_sock << "solution\n" << mesh << elemCent
+             << "window_title 'Element Centers'"
+             << "keys vvv" << flush;
+   }
+
+   return 0;
+}

tests/unit/CMakeLists.txt

@@ -67,6 +67,7 @@ set(UNIT_TESTS_SRCS
   fem/test_fe.cpp
   fem/test_get_value.cpp
   fem/test_getderivative.cpp
+  fem/test_gridfunction_integrals.cpp
   fem/test_intrules.cpp
   fem/test_intruletypes.cpp
   fem/test_inversetransform.cpp

tests/unit/fem/test_gridfunction_integrals.cpp

@@ -0,0 +1,232 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+namespace gridfunction_integrals
+{
+
+double func_1D_lin(const Vector &x)
+{
+   return x[0] - 0.5;
+}
+
+double wgt_1D_lin(const Vector &x)
+{
+   return 7.0 - 3.0 * x[0];
+}
+
+double func_1D_gaussian(const Vector &x)
+{
+   return sqrt(50.0/M_PI) * exp(-50.0 * pow(x[0] - 0.5, 2.0));
+}
+
+class FourierCosine1D : public Coefficient
+{
+private:
+   int n_;
+   double l_;
+   mutable Vector x_;
+
+public:
+   FourierCosine1D(double l) : n_(0), l_(l), x_(1) {}
+
+   void SetMode(int n) { n_ = n; }
+
+   double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      T.Transform(ip, x_);
+
+      return cos(2.0 * M_PI * (double)n_ * x_[0] / l_);
+   }
+};
+
+TEST_CASE("1D GridFunction::ComputeIntegral",
+          "[GridFunction]"
+          "[GridFunction::ComputeIntegral]")
+{
+   int n = 10;
+   int dim = 1;
+   int order = 1;
+
+   FunctionCoefficient funcCoef(func_1D_lin);
+   ConstantCoefficient wgt0Coef(1.0);
+   FunctionCoefficient wgt1Coef(wgt_1D_lin);
+
+   for (int type = (int)Element::SEGMENT;
+        type <= (int)Element::SEGMENT; type++)
+   {
+      Mesh mesh = Mesh::MakeCartesian1D(n, 2.0);
+
+      H1_FECollection h1_fec(order, dim);
+      DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
+                              FiniteElement::VALUE);
+      DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
+                              FiniteElement::INTEGRAL);
+
+      FiniteElementSpace h1_fespace(&mesh, &h1_fec);
+      FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
+      FiniteElementSpace dgi_fespace(&mesh, &dgi_fec);
+
+      GridFunction h1_x(&h1_fespace);
+      GridFunction dgv_x(&dgv_fespace);
+      GridFunction dgi_x(&dgi_fespace);
+
+      h1_x.ProjectCoefficient(funcCoef);
+      dgv_x.ProjectCoefficient(funcCoef);
+      dgi_x.ProjectCoefficient(funcCoef);
+
+      // First integrate with a weight of 1
+      double w0 = 1.0;
+
+      double w0_h1 = h1_x.ComputeIntegral(wgt0Coef, 0);
+      double w0_dgv = dgv_x.ComputeIntegral(wgt0Coef, 0);
+      double w0_dgi = dgi_x.ComputeIntegral(wgt0Coef, 0);
+
+      REQUIRE(w0_h1  == MFEM_Approx(w0));
+      REQUIRE(w0_dgv == MFEM_Approx(w0));
+      REQUIRE(w0_dgi == MFEM_Approx(w0));
+
+      // Integrate with a linear weight function
+      double w1 = 2.0;
+
+      double w1_h1  = h1_x.ComputeIntegral(wgt1Coef, 2);
+      double w1_dgv = dgv_x.ComputeIntegral(wgt1Coef, 2);
+      double w1_dgi = dgi_x.ComputeIntegral(wgt1Coef, 2);
+
+      REQUIRE(w1_h1  == MFEM_Approx(w1));
+      REQUIRE(w1_dgv == MFEM_Approx(w1));
+      REQUIRE(w1_dgi == MFEM_Approx(w1));
+   }
+}
+
+TEST_CASE("1D GridFunction::ComputeIntegral (Fourier)",
+          "[GridFunction]"
+          "[GridFunction::ComputeIntegral]")
+{
+   int n = 20;
+   int dim = 1;
+   int order = 3;
+   double l = 1.0;
+
+   FunctionCoefficient funcCoef(func_1D_gaussian);
+   FourierCosine1D wgtCoef(l);
+
+   for (int type = (int)Element::SEGMENT;
+        type <= (int)Element::SEGMENT; type++)
+   {
+      Mesh mesh = Mesh::MakeCartesian1D(n, l);
+
+      H1_FECollection h1_fec(order, dim);
+      DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
+                              FiniteElement::VALUE);
+      DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
+                              FiniteElement::INTEGRAL);
+
+      FiniteElementSpace h1_fespace(&mesh, &h1_fec);
+      FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
+      FiniteElementSpace dgi_fespace(&mesh, &dgi_fec);
+
+      GridFunction h1_x(&h1_fespace);
+      GridFunction dgv_x(&dgv_fespace);
+      GridFunction dgi_x(&dgi_fespace);
+
+      h1_x.ProjectCoefficient(funcCoef);
+      dgv_x.ProjectCoefficient(funcCoef);
+      dgi_x.ProjectCoefficient(funcCoef);
+
+      LinearForm h1_lf(&h1_fespace);
+      LinearForm dgv_lf(&dgv_fespace);
+      LinearForm dgi_lf(&dgi_fespace);
+
+      h1_lf.AddDomainIntegrator(new DomainLFIntegrator(wgtCoef));
+      dgv_lf.AddDomainIntegrator(new DomainLFIntegrator(wgtCoef));
+      dgi_lf.AddDomainIntegrator(new DomainLFIntegrator(wgtCoef));
+
+      // First integrate with a weight of 1
+      double w0 = 1.0;
+
+      h1_lf.Assemble();
+      dgi_lf.Assemble();
+      dgv_lf.Assemble();
+
+      double i1_w0_h1  = h1_lf(h1_x);
+      double i1_w0_dgv = dgv_lf(dgv_x);
+      double i1_w0_dgi = dgi_lf(dgi_x);
+
+      REQUIRE(i1_w0_h1  == MFEM_Approx(w0, 1e-6));
+      REQUIRE(i1_w0_dgv == MFEM_Approx(w0, 1e-6));
+      REQUIRE(i1_w0_dgi == MFEM_Approx(w0, 1e-6));
+
+      double i2_w0_h1  = h1_x.ComputeIntegral(wgtCoef, 2 * order + 1);
+      double i2_w0_dgv = dgv_x.ComputeIntegral(wgtCoef, 2 * order + 1);
+      double i2_w0_dgi = dgi_x.ComputeIntegral(wgtCoef, 2 * order + 1);
+
+      REQUIRE(i2_w0_h1  == MFEM_Approx(i1_w0_h1));
+      REQUIRE(i2_w0_dgv == MFEM_Approx(i1_w0_dgv));
+      REQUIRE(i2_w0_dgi == MFEM_Approx(i1_w0_dgi));
+
+      // Integrate with a weight of cos(2 pi x)
+      wgtCoef.SetMode(1);
+
+      double w1 = -exp(-M_PI * M_PI / 50.0);
+
+      h1_lf.Assemble();
+      dgi_lf.Assemble();
+      dgv_lf.Assemble();
+
+      double i1_w1_h1  = h1_lf(h1_x);
+      double i1_w1_dgv = dgv_lf(dgv_x);
+      double i1_w1_dgi = dgi_lf(dgi_x);
+
+      REQUIRE(i1_w1_h1  == MFEM_Approx(w1, 1e-6));
+      REQUIRE(i1_w1_dgv == MFEM_Approx(w1, 1e-6));
+      REQUIRE(i1_w1_dgi == MFEM_Approx(w1, 1e-6));
+
+      double i2_w1_h1  = h1_x.ComputeIntegral(wgtCoef, 2 * order + 1);
+      double i2_w1_dgv = dgv_x.ComputeIntegral(wgtCoef, 2 * order + 1);
+      double i2_w1_dgi = dgi_x.ComputeIntegral(wgtCoef, 2 * order + 1);
+
+      REQUIRE(i2_w1_h1  == MFEM_Approx(i1_w1_h1));
+      REQUIRE(i2_w1_dgv == MFEM_Approx(i1_w1_dgv));
+      REQUIRE(i2_w1_dgi == MFEM_Approx(i1_w1_dgi));
+
+      // Integrate with a weight of cos(4 pi x)
+      wgtCoef.SetMode(2);
+
+      double w2 = exp(-4.0 * M_PI * M_PI / 50.0);
+
+      h1_lf.Assemble();
+      dgi_lf.Assemble();
+      dgv_lf.Assemble();
+
+      double i1_w2_h1  = h1_lf(h1_x);
+      double i1_w2_dgv = dgv_lf(dgv_x);
+      double i1_w2_dgi = dgi_lf(dgi_x);
+
+      REQUIRE(i1_w2_h1  == MFEM_Approx(w2, 1e-6));
+      REQUIRE(i1_w2_dgv == MFEM_Approx(w2, 1e-6));
+      REQUIRE(i1_w2_dgi == MFEM_Approx(w2, 1e-6));
+
+      double i2_w2_h1  = h1_x.ComputeIntegral(wgtCoef, 2 * order + 1);
+      double i2_w2_dgv = dgv_x.ComputeIntegral(wgtCoef, 2 * order + 1);
+      double i2_w2_dgi = dgi_x.ComputeIntegral(wgtCoef, 2 * order + 1);
+
+      REQUIRE(i2_w2_h1  == MFEM_Approx(i1_w2_h1));
+      REQUIRE(i2_w2_dgv == MFEM_Approx(i1_w2_dgv));
+      REQUIRE(i2_w2_dgi == MFEM_Approx(i1_w2_dgi));
+   }
+}
+
+} // gridfunction_integrals