Merge pull request #3557 from mfem/mpi_h1lor

Parallel H1 HO <-> LOR Transfer Operator

.gitignore

@@ -296,6 +296,7 @@ miniapps/tools/display-basis
 miniapps/tools/load-dc
 miniapps/tools/convert-dc
 miniapps/tools/lor-transfer
+miniapps/tools/plor-transfer
 miniapps/tools/get-values
 miniapps/tools/check-tmop-metric
 miniapps/tools/tmop-metric-magnitude

CHANGELOG

@@ -45,6 +45,10 @@ New and updated examples and miniapps
 
 - Moved the distance solver methods from miniapps/shifted to miniapps/common.
 
+- Added a new parallel LOR transfer miniapp, miniapps/tools/plor-transfer, which
+  mirrors the functionality of the serial LOR transfer miniapp,
+  miniapps/tools/lor-transfer
+
 Meshing improvements
 --------------------
 - Added support for free connectivity of NURBS patches allowing for more complex
@@ -81,6 +85,9 @@ Discretization improvements
 
 - Added support for p-refined meshes in FindPointsGSLIB.
 
+- Support for parallel transfer of H1 fields using the low-order refined (LOR)
+  transfer operators in L2ProjectionGridTransfer
+
 Linear and nonlinear solvers
 ----------------------------
 - Updated interface to MUMPS direct solver to support multiple right-hand

fem/transfer.cpp

@@ -546,6 +546,234 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
    const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_)
    : L2Projection(fes_ho_, fes_lor_)
 {
+   std::unique_ptr<SparseMatrix> R_mat, M_LH_mat;
+   std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH();
+
+   FiniteElementSpace fes_ho_scalar(fes_ho.GetMesh(), fes_ho.FEColl(), 1);
+   FiniteElementSpace fes_lor_scalar(fes_lor.GetMesh(), fes_lor.FEColl(), 1);
+
+   const SparseMatrix *P_ho = fes_ho_scalar.GetConformingProlongation();
+   const SparseMatrix *P_lor = fes_lor_scalar.GetConformingProlongation();
+
+   if (P_ho || P_lor)
+   {
+      if (P_ho && P_lor)
+      {
+         R_mat.reset(RAP(*P_lor, *R_mat, *P_ho));
+         M_LH_mat.reset(RAP(*P_lor, *M_LH_mat, *P_ho));
+      }
+      else if (P_ho)
+      {
+         R_mat.reset(mfem::Mult(*R_mat, *P_ho));
+         M_LH_mat.reset(mfem::Mult(*M_LH_mat, *P_ho));
+      }
+      else // P_lor != nullptr
+      {
+         R_mat.reset(mfem::Mult(*P_lor, *R_mat));
+         M_LH_mat.reset(mfem::Mult(*P_lor, *M_LH_mat));
+      }
+   }
+
+   SparseMatrix *RTxM_LH_mat = TransposeMult(*R_mat, *M_LH_mat);
+   precon.reset(new DSmoother(*RTxM_LH_mat));
+
+   // Set ownership
+   RTxM_LH.reset(RTxM_LH_mat);
+   R = std::move(R_mat);
+   M_LH = std::move(M_LH_mat);
+
+   SetupPCG();
+}
+
+#ifdef MFEM_USE_MPI
+
+L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
+   const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor)
+   : L2Projection(pfes_ho, pfes_lor),
+     pcg(pfes_ho.GetComm())
+{
+   std::tie(R, M_LH) = ComputeSparseRAndM_LH();
+
+   ParFiniteElementSpace pfes_ho_scalar(pfes_ho.GetParMesh(),
+                                        pfes_ho.FEColl(), 1);
+   ParFiniteElementSpace pfes_lor_scalar(pfes_lor.GetParMesh(),
+                                         pfes_lor.FEColl(), 1);
+
+   HypreParMatrix R_local = HypreParMatrix(pfes_ho.GetComm(),
+                                           pfes_lor_scalar.GlobalVSize(),
+                                           pfes_ho_scalar.GlobalVSize(),
+                                           pfes_lor_scalar.GetDofOffsets(),
+                                           pfes_ho_scalar.GetDofOffsets(),
+                                           static_cast<SparseMatrix*>(R.get()));
+   HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho.GetComm(),
+                                              pfes_lor_scalar.GlobalVSize(),
+                                              pfes_ho_scalar.GlobalVSize(),
+                                              pfes_lor_scalar.GetDofOffsets(),
+                                              pfes_ho_scalar.GetDofOffsets(),
+                                              static_cast<SparseMatrix*>(M_LH.get()));
+
+   HypreParMatrix *R_mat = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(),
+                               &R_local, pfes_ho_scalar.Dof_TrueDof_Matrix());
+   HypreParMatrix *M_LH_mat = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(),
+                                  &M_LH_local, pfes_ho_scalar.Dof_TrueDof_Matrix());
+
+   std::unique_ptr<HypreParMatrix> R_T(R_mat->Transpose());
+   HypreParMatrix *RTxM_LH_mat = ParMult(R_T.get(), M_LH_mat, true);
+
+   HypreBoomerAMG *amg = new HypreBoomerAMG(*RTxM_LH_mat);
+   amg->SetPrintLevel(0);
+
+   R.reset(R_mat);
+   M_LH.reset(M_LH_mat);
+   RTxM_LH.reset(RTxM_LH_mat);
+   precon.reset(amg);
+
+   SetupPCG();
+   pcg.SetPreconditioner(*precon);
+   pcg.SetOperator(*RTxM_LH);
+}
+
+#endif
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::SetupPCG()
+{
+   // Basic PCG solver setup
+   pcg.SetPrintLevel(0);
+   // pcg.SetPrintLevel(IterativeSolver::PrintLevel().Summary());
+   pcg.SetMaxIter(1000);
+   // initial values for relative and absolute tolerance
+   pcg.SetRelTol(1e-13);
+   pcg.SetAbsTol(1e-13);
+   pcg.SetPreconditioner(*precon);
+   pcg.SetOperator(*RTxM_LH);
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult(
+   const Vector& x, Vector& y) const
+{
+   Vector X(fes_ho.GetTrueVSize());
+   Vector X_dim(R->Width());
+
+   Vector Y_dim(R->Height());
+   Vector Y(fes_lor.GetTrueVSize());
+
+   Array<int> vdofs_list;
+
+   GetTDofs(fes_ho, x, X);
+
+   for (int d = 0; d < fes_ho.GetVDim(); ++d)
+   {
+      TDofsListByVDim(fes_ho, d, vdofs_list);
+      X.GetSubVector(vdofs_list, X_dim);
+      R->Mult(X_dim, Y_dim);
+      TDofsListByVDim(fes_lor, d, vdofs_list);
+      Y.SetSubVector(vdofs_list, Y_dim);
+   }
+
+   SetFromTDofs(fes_lor, Y, y);
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose(
+   const Vector& x, Vector& y) const
+{
+   Vector X(fes_lor.GetTrueVSize());
+   Vector X_dim(R->Height());
+
+   Vector Y_dim(R->Width());
+   Vector Y(fes_ho.GetTrueVSize());
+
+   Array<int> vdofs_list;
+
+   GetTDofsTranspose(fes_lor, x, X);
+
+   for (int d = 0; d < fes_ho.GetVDim(); ++d)
+   {
+      TDofsListByVDim(fes_lor, d, vdofs_list);
+      X.GetSubVector(vdofs_list, X_dim);
+      R->MultTranspose(X_dim, Y_dim);
+      TDofsListByVDim(fes_ho, d, vdofs_list);
+      Y.SetSubVector(vdofs_list, Y_dim);
+   }
+
+   SetFromTDofsTranspose(fes_ho, Y, y);
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate(
+   const Vector& x, Vector& y) const
+{
+   Vector X(fes_lor.GetTrueVSize());
+   Vector X_dim(M_LH->Height());
+   Vector Xbar(pcg.Width());
+
+   Vector Y_dim(pcg.Height());
+   Vector Y(fes_ho.GetTrueVSize());
+
+   Array<int> vdofs_list;
+
+   GetTDofs(fes_lor, x, X);
+
+   for (int d = 0; d < fes_ho.GetVDim(); ++d)
+   {
+      TDofsListByVDim(fes_lor, d, vdofs_list);
+      X.GetSubVector(vdofs_list, X_dim);
+      // Compute y = P x = (R^T M_LH)^(-1) M_LH^T X = (R^T M_LH)^(-1) Xbar
+      M_LH->MultTranspose(X_dim, Xbar);
+      Y_dim = 0.0;
+      pcg.Mult(Xbar, Y_dim);
+      TDofsListByVDim(fes_ho, d, vdofs_list);
+      Y.SetSubVector(vdofs_list, Y_dim);
+   }
+
+   SetFromTDofs(fes_ho, Y, y);
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose(
+   const Vector& x, Vector& y) const
+{
+   Vector X(fes_ho.GetTrueVSize());
+   Vector X_dim(pcg.Width());
+   Vector Xbar(pcg.Height());
+
+   Vector Y_dim(M_LH->Height());
+   Vector Y(fes_lor.GetTrueVSize());
+
+   Array<int> vdofs_list;
+
+   GetTDofsTranspose(fes_ho, x, X);
+
+   for (int d = 0; d < fes_ho.GetVDim(); ++d)
+   {
+      TDofsListByVDim(fes_ho, d, vdofs_list);
+      X.GetSubVector(vdofs_list, X_dim);
+      // Compute y = P^T x = M_LH (R^T M_LH)^(-1) X = M_LH Xbar
+      Xbar = 0.0;
+      pcg.Mult(X_dim, Xbar);
+      M_LH->Mult(Xbar, Y_dim);
+      TDofsListByVDim(fes_lor, d, vdofs_list);
+      Y.SetSubVector(vdofs_list, Y_dim);
+   }
+
+   SetFromTDofsTranspose(fes_lor, Y, y);
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_)
+{
+   pcg.SetRelTol(p_rtol_);
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_)
+{
+   pcg.SetAbsTol(p_atol_);
+}
+
+std::pair<
+std::unique_ptr<SparseMatrix>,
+std::unique_ptr<SparseMatrix>>
+                            L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH()
+{
+   std::pair<std::unique_ptr<SparseMatrix>,
+       std::unique_ptr<SparseMatrix>> r_and_mlh;
+
    Mesh* mesh_ho = fes_ho.GetMesh();
    Mesh* mesh_lor = fes_lor.GetMesh();
    int nel_ho = mesh_ho->GetNE();
@@ -553,7 +781,7 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
    int ndof_lor = fes_lor.GetNDofs();
 
    // If the local mesh is empty, skip all computations
-   if (nel_ho == 0) { return; }
+   if (nel_ho == 0) { return {nullptr, nullptr}; }
 
    const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
 
@@ -611,18 +839,26 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
       }
    }
    // DOF by DOF inverse of non-zero entries
-   for (int i = 0; i < ndof_lor; ++i)
-   {
-      ML_inv[i] = 1.0 / ML_inv[i];
-   }
+   LumpedMassInverse(ML_inv);
 
    // Compute sparsity pattern for R = M_L^(-1) M_LH and allocate
-   AllocR();
+   r_and_mlh.first = AllocR();
    // Allocate M_LH (same sparsity pattern as R)
    // L refers to the low-order refined mesh (DOFs correspond to rows)
    // H refers to the higher-order mesh (DOFs correspond to columns)
-   M_LH = SparseMatrix(R.GetI(), R.GetJ(), NULL,
-                       R.Height(), R.Width(), false, true, true);
+   Memory<int> I(r_and_mlh.first->Height() + 1);
+   for (int icol = 0; icol < r_and_mlh.first->Height() + 1; ++icol)
+   {
+      I[icol] = r_and_mlh.first->GetI()[icol];
+   }
+   Memory<int> J(r_and_mlh.first->NumNonZeroElems());
+   for (int jcol = 0; jcol < r_and_mlh.first->NumNonZeroElems(); ++jcol)
+   {
+      J[jcol] = r_and_mlh.first->GetJ()[jcol];
+   }
+   r_and_mlh.second = std::unique_ptr<SparseMatrix>(new SparseMatrix(
+                                                       I, J, NULL,
+                                                       r_and_mlh.first->Height(), r_and_mlh.first->Width(), true, true, true));
 
    IntegrationPointTransformation ip_tr;
    IsoparametricTransformation& emb_tr = ip_tr.Transf;
@@ -667,131 +903,118 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
          }
          Array<int> dofs_ho(nedof_ho);
          fes_ho.GetElementDofs(iho, dofs_ho);
-         M_LH.AddSubMatrix(dofs_lor, dofs_ho, M_LH_el);
-         R.AddSubMatrix(dofs_lor, dofs_ho, R_el);
+         r_and_mlh.second->AddSubMatrix(dofs_lor, dofs_ho, M_LH_el);
+         r_and_mlh.first->AddSubMatrix(dofs_lor, dofs_ho, R_el);
       }
    }
 
-   // Create PCG solver
-   RTxM_LH = TransposeMult(R, M_LH);
-   pcg.SetPrintLevel(0);
-   pcg.SetMaxIter(1000);
-   // initial values for relative and absolute tolerance
-   SetRelTol(1e-13);
-   SetAbsTol(1e-13);
-   Ds = DSmoother(*RTxM_LH);
-   pcg.SetPreconditioner(Ds);
-   pcg.SetOperator(*RTxM_LH);
+   return r_and_mlh;
 }
 
-L2ProjectionGridTransfer::L2ProjectionH1Space::~L2ProjectionH1Space()
+void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofs(
+   const FiniteElementSpace& fes, const Vector& x, Vector& X) const
 {
-   delete RTxM_LH;
-}
-
-void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult(
-   const Vector& x, Vector& y) const
-{
-   int vdim = fes_ho.GetVDim();
-   const int ndof_ho = fes_ho.GetNDofs();
-   const int ndof_lor = fes_lor.GetNDofs();
-   Array<int> dofs_ho(ndof_ho);
-   Array<int> dofs_lor(ndof_lor);
-   Vector x_dim(ndof_ho);
-   Vector y_dim(ndof_lor);
-
-   for (int d = 0; d < vdim; ++d)
+   const Operator* res = fes.GetRestrictionOperator();
+   if (res)
    {
-      fes_ho.GetVDofs(d, dofs_ho);
-      fes_lor.GetVDofs(d, dofs_lor);
-      x.GetSubVector(dofs_ho, x_dim);
-      R.Mult(x_dim, y_dim);
-      y.SetSubVector(dofs_lor, y_dim);
+      res->Mult(x, X);
+   }
+   else
+   {
+      X = x;
    }
 }
 
-void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose(
-   const Vector& x, Vector& y) const
+void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofs(
+   const FiniteElementSpace& fes, const Vector &X, Vector& x) const
 {
-   int vdim = fes_ho.GetVDim();
-   const int ndof_ho = fes_ho.GetNDofs();
-   const int ndof_lor = fes_lor.GetNDofs();
-   Array<int> dofs_ho(ndof_ho);
-   Array<int> dofs_lor(ndof_lor);
-   Vector x_dim(ndof_lor);
-   Vector y_dim(ndof_ho);
-
-   for (int d = 0; d < vdim; ++d)
+   const Operator* P = fes.GetProlongationMatrix();
+   if (P)
    {
-      fes_ho.GetVDofs(d, dofs_ho);
-      fes_lor.GetVDofs(d, dofs_lor);
-      x.GetSubVector(dofs_lor, x_dim);
-      R.MultTranspose(x_dim, y_dim);
-      y.SetSubVector(dofs_ho, y_dim);
+      P->Mult(X, x);
+   }
+   else
+   {
+      x = X;
    }
 }
 
-void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate(
-   const Vector& x, Vector& y) const
+void L2ProjectionGridTransfer::L2ProjectionH1Space::GetTDofsTranspose(
+   const FiniteElementSpace& fes, const Vector& x, Vector& X) const
 {
-   int vdim = fes_ho.GetVDim();
-   const int ndof_ho = fes_ho.GetNDofs();
-   const int ndof_lor = fes_lor.GetNDofs();
-   Array<int> dofs_ho(ndof_ho);
-   Array<int> dofs_lor(ndof_lor);
-   Vector x_dim(ndof_lor);
-   Vector y_dim(ndof_ho);
-   Vector xbar(ndof_ho);
-
-   for (int d = 0; d < vdim; ++d)
+   const Operator* P = fes.GetProlongationMatrix();
+   if (P)
    {
-      fes_lor.GetVDofs(d, dofs_lor);
-      x.GetSubVector(dofs_lor, x_dim);
-      // Compute y = P x = (R^T M_LH)^(-1) M_LH^T x = (R^T M_LH)^(-1) xbar
-      M_LH.MultTranspose(x_dim, xbar);
-      y_dim = 0.0;
-      pcg.Mult(xbar, y_dim);
-      fes_ho.GetVDofs(d, dofs_ho);
-      y.SetSubVector(dofs_ho, y_dim);
+      P->MultTranspose(x, X);
+   }
+   else
+   {
+      X = x;
    }
 }
 
-void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose(
-   const Vector& x, Vector& y) const
+void L2ProjectionGridTransfer::L2ProjectionH1Space::SetFromTDofsTranspose(
+   const FiniteElementSpace& fes, const Vector &X, Vector& x) const
 {
-   int vdim = fes_ho.GetVDim();
-   const int ndof_ho = fes_ho.GetNDofs();
-   const int ndof_lor = fes_lor.GetNDofs();
-   Array<int> dofs_ho(ndof_ho);
-   Array<int> dofs_lor(ndof_lor);
-   Vector x_dim(ndof_ho);
-   Vector y_dim(ndof_lor);
-   Vector xbar(ndof_ho);
-
-   for (int d = 0; d < vdim; ++d)
+   const Operator *R_op = fes.GetRestrictionOperator();
+   if (R_op)
    {
-      fes_ho.GetVDofs(d, dofs_ho);
-      x.GetSubVector(dofs_ho, x_dim);
-      // Compute y = P^T x = M_LH (R^T M_LH)^(-1) x = M_LH xbar
-      xbar = 0.0;
-      pcg.Mult(x_dim, xbar);
-      M_LH.Mult(xbar, y_dim);
-      fes_lor.GetVDofs(d, dofs_lor);
-      y.SetSubVector(dofs_lor, y_dim);
+      R_op->MultTranspose(X, x);
+   }
+   else
+   {
+      x = X;
    }
 }
 
-void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_)
+void L2ProjectionGridTransfer::L2ProjectionH1Space::TDofsListByVDim(
+   const FiniteElementSpace& fes, int vdim, Array<int>& vdofs_list) const
 {
-   pcg.SetRelTol(p_rtol_);
+   const SparseMatrix *R_mat = fes.GetRestrictionMatrix();
+   if (R_mat)
+   {
+      Array<int> x_vdofs_list(fes.GetNDofs());
+      Array<int> x_vdofs_marker(fes.GetVSize());
+      Array<int> X_vdofs_marker(fes.GetTrueVSize());
+      fes.GetVDofs(vdim, x_vdofs_list);
+      FiniteElementSpace::ListToMarker(x_vdofs_list, fes.GetVSize(), x_vdofs_marker);
+      R_mat->BooleanMult(x_vdofs_marker, X_vdofs_marker);
+      FiniteElementSpace::MarkerToList(X_vdofs_marker, vdofs_list);
+   }
+   else
+   {
+      vdofs_list.SetSize(fes.GetNDofs());
+      fes.GetVDofs(vdim, vdofs_list);
+   }
 }
 
-void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_)
+void L2ProjectionGridTransfer::L2ProjectionH1Space::LumpedMassInverse(
+   Vector& ML_inv) const
 {
-   pcg.SetAbsTol(p_atol_);
+   Vector ML_inv_full(fes_lor.GetVSize());
+   // set ML_inv on dofs for vdim = 0
+   Array<int> vdofs_list(fes_lor.GetNDofs());
+   fes_lor.GetVDofs(0, vdofs_list);
+   ML_inv_full.SetSubVector(vdofs_list, ML_inv);
+
+   Vector ML_inv_true(fes_lor.GetTrueVSize());
+   const Operator *P = fes_lor.GetProlongationMatrix();
+   if (P) { P->MultTranspose(ML_inv_full, ML_inv_true); }
+   else { ML_inv_true = ML_inv_full; }
+
+   for (int i = 0; i < ML_inv_true.Size(); ++i)
+   {
+      ML_inv_true[i] = 1.0 / ML_inv_true[i];
+   }
+
+   if (P) { P->Mult(ML_inv_true, ML_inv_full); }
+   else { ML_inv_full = ML_inv_true; }
+
+   ML_inv_full.GetSubVector(vdofs_list, ML_inv);
 }
 
-void L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR()
+std::unique_ptr<SparseMatrix>
+L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR()
 {
    const Table& elem_dof_ho = fes_ho.GetElementToDofTable();
    const Table& elem_dof_lor = fes_lor.GetElementToDofTable();
@@ -871,11 +1094,13 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR()
    dof_lor_dof_ho.SortRows();
    double* data = Memory<double>(dof_dofI[ndof_lor]);
 
-   R = SparseMatrix(dof_dofI, dof_dofJ, data, ndof_lor, ndof_ho,
-                    true, true, true);
-   R = 0.0;
+   std::unique_ptr<SparseMatrix> R_local(new SparseMatrix(
+                                            dof_dofI, dof_dofJ, data, ndof_lor, ndof_ho, true, true, true));
+   (*R_local) = 0.0;
 
    dof_lor_dof_ho.LoseData();
+
+   return R_local;
 }
 
 L2ProjectionGridTransfer::~L2ProjectionGridTransfer()
@@ -905,7 +1130,20 @@ void L2ProjectionGridTransfer::BuildF()
    if (!force_l2_space &&
        dom_fes.FEColl()->GetContType() == FiniteElementCollection::CONTINUOUS)
    {
-      F = new L2ProjectionH1Space(dom_fes, ran_fes);
+      if (!Parallel())
+      {
+         F = new L2ProjectionH1Space(dom_fes, ran_fes);
+      }
+      else
+      {
+#ifdef MFEM_USE_MPI
+         const mfem::ParFiniteElementSpace& dom_pfes =
+            static_cast<mfem::ParFiniteElementSpace&>(dom_fes);
+         const mfem::ParFiniteElementSpace& ran_pfes =
+            static_cast<mfem::ParFiniteElementSpace&>(ran_fes);
+         F = new L2ProjectionH1Space(dom_pfes, ran_pfes);
+#endif
+      }
    }
    else
    {

fem/transfer.hpp

@@ -180,9 +180,15 @@ protected:
    public:
       virtual void Prolongate(const Vector& x, Vector& y) const = 0;
       virtual void ProlongateTranspose(const Vector& x, Vector& y) const = 0;
-      /// Sets relative tolerance and absolute tolerance in preconditioned
-      /// conjugate gradient solver. Only used for H1 spaces.
+      /// @brief Sets relative tolerance in preconditioned conjugate gradient
+      /// solver.
+      ///
+      /// Only used for H1 spaces.
       virtual void SetRelTol(double p_rtol_) = 0;
+      /// @brief Sets absolute tolerance in preconditioned conjugate gradient
+      /// solver.
+      ///
+      /// Only used for H1 spaces.
       virtual void SetAbsTol(double p_atol_) = 0;
    protected:
       const FiniteElementSpace& fes_ho;
@@ -249,29 +255,22 @@ protected:
       /// conservative left-inverse prolongation operation. This functionality
       /// is also provided as an Operator by L2Prolongation.
       virtual void ProlongateTranspose(const Vector& x, Vector& y) const;
-      virtual void SetRelTol(double p_rtol_) {}
-      virtual void SetAbsTol(double p_atol_) {}
+      virtual void SetRelTol(double p_rtol_) { } ///< No-op.
+      virtual void SetAbsTol(double p_atol_) { } ///< No-op.
    };
 
-   /** Class for projection operator between a H1 high-order finite element
-       space on a coarse mesh, and a H1 low-order finite element space on a
-       refined mesh (LOR). */
+   /** Projection operator between a H1 high-order finite element space on a
+       coarse mesh, and a H1 low-order finite element space on a refined mesh
+       (LOR). */
    class L2ProjectionH1Space : public L2Projection
    {
-      // The restriction operator is represented as a SparseMatrix R. The
-      // prolongation operator is a dense matrix computed as the inverse of (R^T
-      // M_L R), and hence, is not stored.
-      SparseMatrix R;
-      // Used to compute P = (RTxM_LH)^(-1) M_LH^T
-      SparseMatrix M_LH;
-      SparseMatrix* RTxM_LH;
-      CGSolver pcg;
-      DSmoother Ds;
-
    public:
-      L2ProjectionH1Space(const FiniteElementSpace& fes_ho_,
-                          const FiniteElementSpace& fes_lor_);
-      virtual ~L2ProjectionH1Space();
+      L2ProjectionH1Space(const FiniteElementSpace &fes_ho_,
+                          const FiniteElementSpace &fes_lor_);
+#ifdef MFEM_USE_MPI
+      L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_,
+                          const ParFiniteElementSpace &pfes_lor_);
+#endif
       /// Maps <tt>x</tt>, primal field coefficients defined on a coarse mesh
       /// with a higher order H1 finite element space, to <tt>y</tt>, primal
       /// field coefficients defined on a refined mesh with a low order H1
@@ -305,11 +304,51 @@ protected:
       virtual void ProlongateTranspose(const Vector& x, Vector& y) const;
       virtual void SetRelTol(double p_rtol_);
       virtual void SetAbsTol(double p_atol_);
-   private:
-      /// Computes sparsity pattern and initializes R matrix. Based on
-      /// BilinearForm::AllocMat() except maps between HO elements and LOR
-      /// elements.
-      void AllocR();
+   protected:
+      /// Sets up the PCG solver (sets parameters, operator, and preconditioner)
+      void SetupPCG();
+      /// Computes on-rank R and M_LH matrices.
+      std::pair<std::unique_ptr<SparseMatrix>,
+          std::unique_ptr<SparseMatrix>> ComputeSparseRAndM_LH();
+      /// @brief Recovers vector of tdofs given a vector of dofs and a finite
+      /// element space
+      void GetTDofs(const FiniteElementSpace& fes, const Vector& x, Vector& X) const;
+      /// Sets dof values given a vector of tdofs and a finite element space
+      void SetFromTDofs(const FiniteElementSpace& fes,
+                        const Vector& X,
+                        Vector& x) const;
+      /// @brief Recovers a vector of dual field coefficients on the tdofs given
+      /// a vector of dual coefficients and a finite element space
+      void GetTDofsTranspose(const FiniteElementSpace& fes,
+                             const Vector& x,
+                             Vector& X) const;
+      /// @brief Sets dual field coefficients given a vector of dual field
+      /// coefficients on the tdofs and a finite element space
+      void SetFromTDofsTranspose(const FiniteElementSpace& fes,
+                                 const Vector& X,
+                                 Vector& x) const;
+      /// @brief Fills the vdofs_list array with a list of vdofs for a given
+      /// vdim and a given finite element space
+      void TDofsListByVDim(const FiniteElementSpace& fes,
+                           int vdim,
+                           Array<int>& vdofs_list) const;
+      /// Returns the inverse of an on-rank lumped mass matrix
+      void LumpedMassInverse(Vector& ML_inv) const;
+      /// @brief Computes sparsity pattern and initializes R matrix.
+      ///
+      /// Based on BilinearForm::AllocMat(), except maps between coarse HO
+      /// elements and refined LOR elements.
+      std::unique_ptr<SparseMatrix> AllocR();
+
+      CGSolver pcg;
+      std::unique_ptr<Solver> precon;
+      // The restriction operator is represented as an Operator R. The
+      // prolongation operator is a dense matrix computed as the inverse of (R^T
+      // M_L R), and hence, is not stored.
+      std::unique_ptr<Operator> R;
+      // Used to compute P = (RT*M_LH)^(-1) M_LH^T
+      std::unique_ptr<Operator> M_LH;
+      std::unique_ptr<Operator> RTxM_LH;
    };
 
    /** Mass-conservative prolongation operator going in the opposite direction

miniapps/tools/CMakeLists.txt

@@ -28,5 +28,10 @@ add_mfem_miniapp(convert-dc
 add_mfem_miniapp(lor-transfer
   MAIN lor-transfer.cpp LIBRARIES mfem)
 
+if (MFEM_USE_MPI)
+  add_mfem_miniapp(plor-transfer
+    MAIN plor-transfer.cpp LIBRARIES mfem)
+endif()
+
 add_mfem_miniapp(check-tmop-metric
   MAIN check-tmop-metric.cpp LIBRARIES mfem)

miniapps/tools/lor-transfer.cpp

@@ -95,13 +95,7 @@ int main(int argc, char *argv[])
    args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
                   "-no-t", "--dont-use-pointwise-transfer",
                   "Use pointwise transfer operators instead of L2 projection.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-   args.PrintOptions(cout);
+   args.ParseCheck();
 
    // Read the mesh from the given mesh file.
    Mesh mesh(mesh_file, 1, 1);
@@ -157,6 +151,10 @@ int main(int argc, char *argv[])
    direction = "HO -> LOR @ HO";
    FunctionCoefficient RHO(RHO_exact);
    rho.ProjectCoefficient(RHO);
+   // Make sure AMR constraints are satisfied
+   rho.SetTrueVector();
+   rho.SetFromTrueVector();
+
    double ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO       ");
    if (vis) { visualize(HO_dc, "HO", Wx, Wy); Wx += offx; }
 
@@ -193,17 +191,13 @@ int main(int argc, char *argv[])
    }
 
    // HO* to LOR* dual fields
-   GridFunction ones(&fespace), ones_lor(&fespace_lor);
-   ones = 1.0;
-   ones_lor = 1.0;
    LinearForm M_rho(&fespace), M_rho_lor(&fespace_lor);
    if (!use_pointwise_transfer && gt->SupportsBackwardsOperator())
    {
       const Operator &P = gt->BackwardOperator();
       M_ho.Mult(rho, M_rho);
       P.MultTranspose(M_rho, M_rho_lor);
-      cout << "HO -> LOR dual field: " << fabs(M_rho(ones)-M_rho_lor(ones_lor))
-           << endl << endl;
+      cout << "HO -> LOR dual field: " << abs(M_rho.Sum()-M_rho_lor.Sum()) << "\n\n";
    }
 
    // LOR projections
@@ -239,8 +233,7 @@ int main(int argc, char *argv[])
    {
       M_lor.Mult(rho_lor, M_rho_lor);
       R.MultTranspose(M_rho_lor, M_rho);
-      cout << "LOR -> HO dual field: " << fabs(M_rho(ones)-M_rho_lor(ones_lor))
-           << '\n';
+      cout << "LOR -> HO dual field: " << abs(M_rho.Sum() - M_rho_lor.Sum()) << '\n';
    }
 
    delete fec;
@@ -288,14 +281,11 @@ double compute_mass(FiniteElementSpace *L2, double massL2,
                     VisItDataCollection &dc, string prefix)
 {
    ConstantCoefficient one(1.0);
-   BilinearForm ML2(L2);
-   ML2.AddDomainIntegrator(new MassIntegrator(one));
-   ML2.Assemble();
+   LinearForm lf(L2);
+   lf.AddDomainIntegrator(new DomainLFIntegrator(one));
+   lf.Assemble();
 
-   GridFunction rhoone(L2);
-   rhoone = 1.0;
-
-   double newmass = ML2.InnerProduct(*dc.GetField("density"),rhoone);
+   double newmass = lf(*dc.GetField("density"));
    cout.precision(18);
    cout << space << " " << prefix << " mass   = " << newmass;
    if (massL2 >= 0)

miniapps/tools/makefile

@@ -28,7 +28,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer \
 	check-tmop-metric tmop-metric-magnitude
-PAR_MINIAPPS =
+PAR_MINIAPPS = plor-transfer
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -81,9 +81,9 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 	@$(call mfem-test,$<,, Tools miniapp)
 
 # Testing: Specific execution options
-# Do not test: display-basis, load-dc, convert-dc, get-values, lor-transfer
+# Do not test: display-basis, load-dc, convert-dc, get-values, lor-transfer, plor-transfer
 NO_TEST_APPS = display-basis load-dc convert-dc get-values lor-transfer \
-	check-tmop-metric tmop-metric-magnitude
+	plor-transfer check-tmop-metric tmop-metric-magnitude
 $(foreach app,$(NO_TEST_APPS),$(app)-test-seq $(app)-test-par):
 	@true
 

miniapps/tools/plor-transfer.cpp

@@ -0,0 +1,360 @@
+// Copyright (c) 2010-2023, 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.
+//
+//   -----------------------------------------------------------------------
+//   Parallel LOR Transfer Miniapp:  Map functions between HO and LOR spaces
+//   -----------------------------------------------------------------------
+//
+// This miniapp visualizes the maps between a high-order (HO) finite element
+// space, typically using high-order functions on a high-order mesh, and a
+// low-order refined (LOR) finite element space, typically defined by 0th or 1st
+// order functions on a low-order refinement of the HO mesh.
+//
+// The grid transfer operators are represented using either
+// InterpolationGridTransfer or L2ProjectionGridTransfer (depending on the
+// options requested by the user). The two transfer operators are then:
+//
+//  1. R: HO -> LOR, defined by GridTransfer::ForwardOperator
+//  2. P: LOR -> HO, defined by GridTransfer::BackwardOperator
+//
+// While defined generally, these operators have some nice properties for
+// particular finite element spaces. For example they satisfy PR=I, plus mass
+// conservation in both directions for L2 fields.
+//
+// Compile with: make plor-transfer
+//
+// Sample runs:  plor-transfer
+//               plor-transfer -h1
+//               plor-transfer -t
+//               plor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4
+//               plor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2
+//               plor-transfer -lref 4 -o 4 -lo 0 -p 1
+//               plor-transfer -lref 5 -o 4 -lo 0 -p 1
+//               plor-transfer -lref 5 -o 4 -lo 3 -p 2
+//               plor-transfer -lref 5 -o 4 -lo 0 -p 3
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int problem = 1; // problem type
+
+int Wx = 0, Wy = 0; // window position
+int Ww = 350, Wh = 350; // window size
+int offx = Ww+5, offy = Wh+25; // window offsets
+
+string space;
+string direction;
+
+// Exact functions to project
+double RHO_exact(const Vector &x);
+
+// Helper functions
+void visualize(VisItDataCollection &, string, int, int);
+double compute_mass(ParFiniteElementSpace *, double, VisItDataCollection &,
+                    string);
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 3;
+   int lref = order+1;
+   int lorder = 0;
+   bool vis = true;
+   bool useH1 = false;
+   bool use_pointwise_transfer = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem type (see the RHO_exact function).");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&lref, "-lref", "--lor-ref-level", "LOR refinement level.");
+   args.AddOption(&lorder, "-lo", "--lor-order",
+                  "LOR space order (polynomial degree, zero by default).");
+   args.AddOption(&vis, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&useH1, "-h1", "--use-h1", "-l2", "--use-l2",
+                  "Use H1 spaces instead of L2.");
+   args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
+                  "-no-t", "--dont-use-pointwise-transfer",
+                  "Use pointwise transfer operators instead of L2 projection.");
+   args.ParseCheck();
+
+   // Read the mesh from the given mesh file.
+   Mesh serial_mesh(mesh_file, 1, 1);
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+   int dim = mesh.Dimension();
+
+   // Create the low-order refined mesh
+   int basis_lor = BasisType::GaussLobatto; // BasisType::ClosedUniform;
+   ParMesh mesh_lor = ParMesh::MakeRefined(mesh, lref, basis_lor);
+
+   // Create spaces
+   FiniteElementCollection *fec, *fec_lor;
+   if (useH1)
+   {
+      space = "H1";
+      if (lorder == 0)
+      {
+         lorder = 1;
+         if (Mpi::Root())
+         {
+            cerr << "Switching the H1 LOR space order from 0 to 1\n";
+         }
+      }
+      fec = new H1_FECollection(order, dim);
+      fec_lor = new H1_FECollection(lorder, dim);
+   }
+   else
+   {
+      space = "L2";
+      fec = new L2_FECollection(order, dim);
+      fec_lor = new L2_FECollection(lorder, dim);
+   }
+
+   ParFiniteElementSpace fespace(&mesh, fec);
+   ParFiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
+
+   ParGridFunction rho(&fespace);
+   ParGridFunction rho_lor(&fespace_lor);
+
+   // Data collections for vis/analysis
+   VisItDataCollection HO_dc(MPI_COMM_WORLD, "HO", &mesh);
+   HO_dc.RegisterField("density", &rho);
+   VisItDataCollection LOR_dc(MPI_COMM_WORLD, "LOR", &mesh_lor);
+   LOR_dc.RegisterField("density", &rho_lor);
+
+   ParBilinearForm M_ho(&fespace);
+   M_ho.AddDomainIntegrator(new MassIntegrator);
+   M_ho.Assemble();
+   M_ho.Finalize();
+   HypreParMatrix* M_ho_tdof = M_ho.ParallelAssemble();
+
+   ParBilinearForm M_lor(&fespace_lor);
+   M_lor.AddDomainIntegrator(new MassIntegrator);
+   M_lor.Assemble();
+   M_lor.Finalize();
+   HypreParMatrix* M_lor_tdof = M_lor.ParallelAssemble();
+
+   // HO projections
+   direction = "HO -> LOR @ HO";
+   FunctionCoefficient RHO(RHO_exact);
+   rho.ProjectCoefficient(RHO);
+   // Make sure AMR constraints are satisfied
+   rho.SetTrueVector();
+   rho.SetFromTrueVector();
+
+   double ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO       ");
+   if (vis) { visualize(HO_dc, "HO", Wx, Wy); Wx += offx; }
+
+   GridTransfer *gt;
+   if (use_pointwise_transfer)
+   {
+      gt = new InterpolationGridTransfer(fespace, fespace_lor);
+   }
+   else
+   {
+      gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
+   }
+   const Operator &R = gt->ForwardOperator();
+
+   // HO->LOR restriction
+   direction = "HO -> LOR @ LOR";
+   R.Mult(rho, rho_lor);
+   compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO)    ");
+   if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy); Wx += offx; }
+   auto global_max = [](const Vector& v)
+   {
+      double max = v.Normlinf();
+      MPI_Allreduce(MPI_IN_PLACE, &max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+      return max;
+   };
+
+   if (gt->SupportsBackwardsOperator())
+   {
+      const Operator &P = gt->BackwardOperator();
+      // LOR->HO prolongation
+      direction = "HO -> LOR @ HO";
+      ParGridFunction rho_prev = rho;
+      P.Mult(rho_lor, rho);
+      compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
+      if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy); Wx = 0; Wy += offy; }
+
+      rho_prev -= rho;
+      Vector rho_prev_true(fespace.GetTrueVSize());
+      rho_prev.GetTrueDofs(rho_prev_true);
+      double l_inf = global_max(rho_prev_true);
+      if (Mpi::Root())
+      {
+         cout.precision(12);
+         cout << "|HO - P(R(HO))|_∞   = " << l_inf << endl;
+      }
+   }
+
+   // HO* to LOR* dual fields
+   ParLinearForm M_rho(&fespace), M_rho_lor(&fespace_lor);
+   auto global_sum = [](const Vector& v)
+   {
+      double sum = v.Sum();
+      MPI_Allreduce(MPI_IN_PLACE, &sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+      return sum;
+   };
+   if (!use_pointwise_transfer && gt->SupportsBackwardsOperator())
+   {
+      Vector M_rho_true(fespace.GetTrueVSize());
+      M_ho_tdof->Mult(rho.GetTrueVector(), M_rho_true);
+      fespace.GetRestrictionOperator()->MultTranspose(M_rho_true, M_rho);
+      const Operator &P = gt->BackwardOperator();
+      P.MultTranspose(M_rho, M_rho_lor);
+      double ho_dual_mass = global_sum(M_rho);
+      double lor_dual_mass = global_sum(M_rho_lor);
+      if (Mpi::Root())
+      {
+         cout << "HO -> LOR dual field: " << abs(ho_dual_mass - lor_dual_mass) << "\n\n";
+      }
+   }
+
+   // LOR projections
+   direction = "LOR -> HO @ LOR";
+   rho_lor.ProjectCoefficient(RHO);
+   ParGridFunction rho_lor_prev = rho_lor;
+   double lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR      ");
+   if (vis) { visualize(LOR_dc, "LOR", Wx, Wy); Wx += offx; }
+
+   if (gt->SupportsBackwardsOperator())
+   {
+      const Operator &P = gt->BackwardOperator();
+      // Prolongate to HO space
+      direction = "LOR -> HO @ HO";
+      P.Mult(rho_lor, rho);
+      compute_mass(&fespace, lor_mass, HO_dc, "P(LOR)   ");
+      if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy); Wx += offx; }
+
+      // Restrict back to LOR space. This won't give the original function because
+      // the rho_lor doesn't necessarily live in the range of R.
+      direction = "LOR -> HO @ LOR";
+      R.Mult(rho, rho_lor);
+      compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
+      if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy); }
+
+      rho_lor_prev -= rho_lor;
+      Vector rho_lor_prev_true(fespace_lor.GetTrueVSize());
+      rho_lor_prev.GetTrueDofs(rho_lor_prev_true);
+      double l_inf = global_max(rho_lor_prev_true);
+      if (Mpi::Root())
+      {
+         cout.precision(12);
+         cout << "|LOR - R(P(LOR))|_∞ = " << l_inf << endl;
+      }
+   }
+
+   // LOR* to HO* dual fields
+   if (!use_pointwise_transfer)
+   {
+      Vector M_rho_lor_true(fespace_lor.GetTrueVSize());
+      M_lor_tdof->Mult(rho_lor.GetTrueVector(), M_rho_lor_true);
+      fespace_lor.GetRestrictionOperator()->MultTranspose(M_rho_lor_true,
+                                                          M_rho_lor);
+      R.MultTranspose(M_rho_lor, M_rho);
+      double ho_dual_mass = global_sum(M_rho);
+      double lor_dual_mass = global_sum(M_rho_lor);
+
+      cout << lor_dual_mass << '\n';
+      cout << ho_dual_mass << '\n';
+
+      if (Mpi::Root())
+      {
+         cout << "LOR -> HO dual field: " << abs(ho_dual_mass - lor_dual_mass) << '\n';
+      }
+   }
+
+   delete fec;
+   delete fec_lor;
+   delete M_ho_tdof;
+   delete M_lor_tdof;
+   delete gt;
+
+   return 0;
+}
+
+
+double RHO_exact(const Vector &x)
+{
+   switch (problem)
+   {
+      case 1: // smooth field
+         return x(1)+0.25*cos(2*M_PI*x.Norml2());
+      case 2: // cubic function
+         return x(1)*x(1)*x(1) + 2*x(0)*x(1) + x(0);
+      case 3: // sharp gradient
+         return M_PI/2-atan(5*(2*x.Norml2()-1));
+      case 4: // basis function
+         return (x.Norml2() < 0.1) ? 1 : 0;
+      default:
+         return 1.0;
+   }
+}
+
+
+void visualize(VisItDataCollection &dc, string prefix, int x, int y)
+{
+   int w = Ww, h = Wh;
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   socketstream sol_sockL2(vishost, visport);
+   sol_sockL2 << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() <<
+              "\n";
+   sol_sockL2.precision(8);
+   sol_sockL2 << "solution\n" << *dc.GetMesh() << *dc.GetField("density")
+              << "window_geometry " << x << " " << y << " " << w << " " << h
+              << "plot_caption '" << space << " " << prefix << " Density'"
+              << "window_title '" << direction << "'" << flush;
+}
+
+
+double compute_mass(ParFiniteElementSpace *L2, double massL2,
+                    VisItDataCollection &dc, string prefix)
+{
+   ConstantCoefficient one(1.0);
+   ParLinearForm lf(L2);
+   lf.AddDomainIntegrator(new DomainLFIntegrator(one));
+   lf.Assemble();
+
+   double newmass = lf(*dc.GetParField("density"));
+   if (Mpi::Root())
+   {
+      cout.precision(18);
+      cout << space << " " << prefix << " mass   = " << newmass;
+      if (massL2 >= 0)
+      {
+         cout.precision(4);
+         cout << " ("  << fabs(newmass-massL2)*100/massL2 << "%)";
+      }
+      cout << endl;
+   }
+   return newmass;
+}