Merge branch 'master' into lor-consistent-mass-dev

fem/transfer.cpp

@@ -1030,12 +1030,42 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::EAProlongateTranspose(
    BatchedLinAlg::MultTranspose(P_dt, x, y);
 }
 
+L2ProjectionGridTransfer::L2ProjectionH1Space::H1ConsistentMassOperator::
+H1ConsistentMassOperator(const Operator &M_LH_, const Solver &M_L_solver_)
+   : Operator(M_LH_.Height(), M_LH_.Width()),
+     M_LH(M_LH_),
+     M_L_solver(M_L_solver_)
+{
+   MFEM_VERIFY(M_LH.Height() == M_L_solver.Height() &&
+               M_LH.Height() == M_L_solver.Width(),
+               "incompatible consistent mass operator dimensions");
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::H1ConsistentMassOperator::
+Mult(const Vector &x, Vector &y) const
+{
+   Vector tmp(M_LH.Height());
+   M_LH.Mult(x, tmp);
+   M_L_solver.Mult(tmp, y);
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::H1ConsistentMassOperator::
+MultTranspose(const Vector &x, Vector &y) const
+{
+   Vector tmp(M_LH.Height());
+   M_L_solver.Mult(x, tmp);
+   M_LH.MultTranspose(tmp, y);
+}
+
 L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
    const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_,
-   const bool use_ea_, MemoryType d_mt_)
+   const bool use_ea_, const bool use_consistent_mass_, MemoryType d_mt_)
    : L2Projection(fes_ho_, fes_lor_, d_mt_),
-     use_ea(use_ea_)
+     use_ea(use_ea_),
+     use_consistent_mass(use_consistent_mass_)
 {
+   MFEM_VERIFY(!(use_ea && use_consistent_mass),
+               "consistent mass is not supported with element assembly");
 
    // need scalar to keep dimensions matching (operators are built to apply
    // individually on each vdim)
@@ -1053,7 +1083,7 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
 
    std::unique_ptr<SparseMatrix> R_mat, M_LH_mat;
 
-   std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH();
+   std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH(!use_consistent_mass);
 
    const SparseMatrix *P_ho = fes_ho_scalar->GetConformingProlongation();
    const SparseMatrix *P_lor = fes_lor_scalar->GetConformingProlongation();
@@ -1062,40 +1092,71 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
    {
       if (P_ho && P_lor)
       {
-         R_mat.reset(RAP(*P_lor, *R_mat, *P_ho));
+         if (R_mat) { 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));
+         if (R_mat) { 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));
+         if (R_mat) { 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));
+   if (use_consistent_mass)
+   {
+      BilinearForm M_lor(fes_lor_scalar.get());
+      M_lor.AddDomainIntegrator(new MassIntegrator);
+      M_lor.Assemble();
+      M_lor.Finalize();
+      SparseMatrix *M_L_mat = M_lor.LoseMat();
 
-   // Set ownership
-   RTxM_LH.reset(RTxM_LH_mat);
-   R = std::move(R_mat);
-   M_LH = std::move(M_LH_mat);
+      ML_precon.reset(new DSmoother(*M_L_mat));
+      ML_pcg.SetPrintLevel(0);
+      ML_pcg.SetMaxIter(1000);
+      ML_pcg.SetRelTol(1e-13);
+      ML_pcg.SetAbsTol(1e-13);
+      ML_pcg.SetPreconditioner(*ML_precon);
+      ML_pcg.SetOperator(*M_L_mat);
+      // Start each solve from zero so repeated Operator::Mult() calls do not
+      // depend on the output vector contents supplied by the caller.
+      ML_pcg.iterative_mode = false;
 
-   SetupPCG();
+      M_L.reset(M_L_mat);
+      M_LH = std::move(M_LH_mat);
+      R.reset(new H1ConsistentMassOperator(*M_LH, ML_pcg));
+   }
+   else
+   {
+      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,
-   const bool use_ea_, MemoryType d_mt_)
+   const bool use_ea_, const bool use_consistent_mass_, MemoryType d_mt_)
    : L2Projection(pfes_ho, pfes_lor, d_mt_),
-     use_ea(use_ea_), pcg(pfes_ho.GetComm())
+     use_ea(use_ea_),
+     use_consistent_mass(use_consistent_mass_),
+     ML_pcg(pfes_ho.GetComm()),
+     pcg(pfes_ho.GetComm())
 {
+   MFEM_VERIFY(!(use_ea && use_consistent_mass),
+               "consistent mass is not supported with element assembly");
 
    // need scalar to keep dimensions matching (operators are built to apply
    // individually on each vdim)
@@ -1111,8 +1172,42 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
       return;
    }
 
-   std::tie(R, M_LH) = ComputeSparseRAndM_LH();
+   std::tie(R, M_LH) = ComputeSparseRAndM_LH(!use_consistent_mass);
 
+   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 *M_LH_mat = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(),
+                                  &M_LH_local, pfes_ho_scalar->Dof_TrueDof_Matrix());
+
+   if (use_consistent_mass)
+   {
+      ParBilinearForm M_lor(pfes_lor_scalar.get());
+      M_lor.AddDomainIntegrator(new MassIntegrator);
+      M_lor.Assemble();
+      M_lor.Finalize();
+      HypreParMatrix *M_L_mat = M_lor.ParallelAssemble();
+
+      M_L.reset(M_L_mat);
+      M_LH.reset(M_LH_mat);
+      HypreDiagScale *ML_hypre_precon = new HypreDiagScale(*M_L_mat);
+      HyprePCG *ML_hypre_pcg = new HyprePCG(*M_L_mat);
+      ML_hypre_pcg->SetPrintLevel(0);
+      ML_hypre_pcg->SetMaxIter(1000);
+      ML_hypre_pcg->SetTol(1e-13);
+      ML_hypre_pcg->SetAbsTol(1e-13);
+      ML_hypre_pcg->SetPreconditioner(*ML_hypre_precon);
+      // Start each solve from zero so repeated Operator::Mult() calls do not
+      // depend on the output vector contents supplied by the caller.
+      ML_hypre_pcg->SetZeroInitialIterate();
+      ML_precon.reset(ML_hypre_precon);
+      ML_solver.reset(ML_hypre_pcg);
+      R.reset(new H1ConsistentMassOperator(*M_LH, *ML_solver));
+      return;
+   }
 
    HypreParMatrix R_local = HypreParMatrix(pfes_ho.GetComm(),
                                            pfes_lor_scalar->GlobalVSize(),
@@ -1120,17 +1215,9 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
                                            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);
@@ -1438,6 +1525,8 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose(
 void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate(
    const Vector& x, Vector& y) const
 {
+   MFEM_VERIFY(!use_consistent_mass,
+               "BackwardOperator is not supported with consistent mass");
 
    Vector X(fes_lor.GetTrueVSize());
    Vector X_dim(M_LH->Height());
@@ -1469,6 +1558,9 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate(
 void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose(
    const Vector& x, Vector& y) const
 {
+   MFEM_VERIFY(!use_consistent_mass,
+               "BackwardOperator is not supported with consistent mass");
+
    Vector X(fes_ho.GetTrueVSize());
    Vector X_dim(pcg.Width());
    Vector Xbar(pcg.Height());
@@ -1499,17 +1591,34 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose(
 void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(real_t p_rtol_)
 {
    pcg.SetRelTol(p_rtol_);
+   ML_pcg.SetRelTol(p_rtol_);
+#ifdef MFEM_USE_MPI
+   if (ML_solver)
+   {
+      HyprePCG *hypre_pcg = dynamic_cast<HyprePCG*>(ML_solver.get());
+      if (hypre_pcg) { hypre_pcg->SetTol(p_rtol_); }
+   }
+#endif
 }
 
 void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(real_t p_atol_)
 {
    pcg.SetAbsTol(p_atol_);
+   ML_pcg.SetAbsTol(p_atol_);
+#ifdef MFEM_USE_MPI
+   if (ML_solver)
+   {
+      HyprePCG *hypre_pcg = dynamic_cast<HyprePCG*>(ML_solver.get());
+      if (hypre_pcg) { hypre_pcg->SetAbsTol(p_atol_); }
+   }
+#endif
 }
 
 std::pair<
 std::unique_ptr<SparseMatrix>,
 std::unique_ptr<SparseMatrix>>
-                            L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH()
+                            L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH(
+                               bool build_R)
 {
    std::pair<std::unique_ptr<SparseMatrix>,
        std::unique_ptr<SparseMatrix>> r_and_mlh;
@@ -1523,10 +1632,10 @@ std::unique_ptr<SparseMatrix>>
    // If the local mesh is empty, skip all computations
    if (nel_ho == 0)
    {
-      return std::make_pair(
-                std::unique_ptr<SparseMatrix>(new SparseMatrix),
-                std::unique_ptr<SparseMatrix>(new SparseMatrix)
-             );
+      std::unique_ptr<SparseMatrix> R_empty;
+      if (build_R) { R_empty.reset(new SparseMatrix); }
+      std::unique_ptr<SparseMatrix> M_LH_empty(new SparseMatrix);
+      return std::make_pair(std::move(R_empty), std::move(M_LH_empty));
    }
 
    const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
@@ -1542,69 +1651,76 @@ std::unique_ptr<SparseMatrix>>
 
    BuildHo2Lor(nel_ho, nel_lor, cf_tr);
 
-   // ML_inv contains the inverse lumped (row sum) mass matrix. Note that the
-   // method will also work with a full (consistent) mass matrix, though this is
-   // not implemented here. L refers to the low-order refined mesh
    Vector ML_inv(ndof_lor);
-   ML_inv = 0.0;
-
-   // Compute ML_inv
-   for (int iho = 0; iho < nel_ho; ++iho)
+   if (build_R)
    {
-      Array<int> lor_els;
-      ho2lor.GetRow(iho, lor_els);
-      int nref = ho2lor.RowSize(iho);
+      // ML_inv contains the inverse lumped (row sum) mass matrix. L refers to
+      // the low-order refined mesh.
+      ML_inv = 0.0;
 
-      Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
-      const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
-      int nedof_lor = fe_lor.GetDof();
-
-      // Instead of using a MassIntegrator, manually loop over integration
-      // points so we can row sum and store the diagonal as a Vector.
-      Vector ML_el(nedof_lor);
-      Vector shape_lor(nedof_lor);
-      Array<int> dofs_lor(nedof_lor);
-
-      for (int iref = 0; iref < nref; ++iref)
+      // Compute ML_inv
+      for (int iho = 0; iho < nel_ho; ++iho)
       {
-         int ilor = lor_els[iref];
-         ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
+         Array<int> lor_els;
+         ho2lor.GetRow(iho, lor_els);
+         int nref = ho2lor.RowSize(iho);
 
-         int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
-         const IntegrationRule* ir = &IntRules.Get(geom, order);
-         ML_el = 0.0;
-         for (int i = 0; i < ir->GetNPoints(); ++i)
+         Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
+         const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
+         int nedof_lor = fe_lor.GetDof();
+
+         // Instead of using a MassIntegrator, manually loop over integration
+         // points so we can row sum and store the diagonal as a Vector.
+         Vector ML_el(nedof_lor);
+         Vector shape_lor(nedof_lor);
+         Array<int> dofs_lor(nedof_lor);
+
+         for (int iref = 0; iref < nref; ++iref)
          {
-            const IntegrationPoint& ip_lor = ir->IntPoint(i);
-            fe_lor.CalcShape(ip_lor, shape_lor);
-            el_tr->SetIntPoint(&ip_lor);
-            ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
+            int ilor = lor_els[iref];
+            ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
+
+            int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
+            const IntegrationRule* ir = &IntRules.Get(geom, order);
+            ML_el = 0.0;
+            for (int i = 0; i < ir->GetNPoints(); ++i)
+            {
+               const IntegrationPoint& ip_lor = ir->IntPoint(i);
+               fe_lor.CalcShape(ip_lor, shape_lor);
+               el_tr->SetIntPoint(&ip_lor);
+               ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
+            }
+            fes_lor.GetElementDofs(ilor, dofs_lor);
+            ML_inv.AddElementVector(dofs_lor, ML_el);
          }
-         fes_lor.GetElementDofs(ilor, dofs_lor);
-         ML_inv.AddElementVector(dofs_lor, ML_el);
       }
+      // DOF by DOF inverse of non-zero entries
+      LumpedMassInverse(ML_inv);
    }
-   // DOF by DOF inverse of non-zero entries
-   LumpedMassInverse(ML_inv);
 
    // Compute sparsity pattern for R = M_L^(-1) M_LH and allocate
-   r_and_mlh.first = AllocR();
+   std::unique_ptr<SparseMatrix> pattern = AllocR();
+   if (build_R)
+   {
+      r_and_mlh.first = std::move(pattern);
+   }
    // 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)
-   Memory<int> I(r_and_mlh.first->Height() + 1);
-   for (int icol = 0; icol < r_and_mlh.first->Height() + 1; ++icol)
+   SparseMatrix &pattern_mat = build_R ? *r_and_mlh.first : *pattern;
+   Memory<int> I(pattern_mat.Height() + 1);
+   for (int icol = 0; icol < pattern_mat.Height() + 1; ++icol)
    {
-      I[icol] = r_and_mlh.first->GetI()[icol];
+      I[icol] = pattern_mat.GetI()[icol];
    }
-   Memory<int> J(r_and_mlh.first->NumNonZeroElems());
-   for (int jcol = 0; jcol < r_and_mlh.first->NumNonZeroElems(); ++jcol)
+   Memory<int> J(pattern_mat.NumNonZeroElems());
+   for (int jcol = 0; jcol < pattern_mat.NumNonZeroElems(); ++jcol)
    {
-      J[jcol] = r_and_mlh.first->GetJ()[jcol];
+      J[jcol] = pattern_mat.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));
+                         new SparseMatrix(I, J, NULL, pattern_mat.Height(),
+                                          pattern_mat.Width(), true, true, true));
 
    IntegrationPointTransformation ip_tr;
    IsoparametricTransformation& emb_tr = ip_tr.Transf;
@@ -1647,15 +1763,21 @@ std::unique_ptr<SparseMatrix>>
          Array<int> dofs_lor(nedof_lor);
          fes_lor.GetElementDofs(ilor, dofs_lor);
          Vector R_row;
-         for (int i = 0; i < nedof_lor; ++i)
+         if (build_R)
          {
-            M_LH_el.GetRow(i, R_row);
-            R_el.SetRow(i, R_row.Set(ML_inv[dofs_lor[i]], R_row));
+            for (int i = 0; i < nedof_lor; ++i)
+            {
+               M_LH_el.GetRow(i, R_row);
+               R_el.SetRow(i, R_row.Set(ML_inv[dofs_lor[i]], R_row));
+            }
          }
          Array<int> dofs_ho(nedof_ho);
          fes_ho.GetElementDofs(iho, dofs_ho);
          r_and_mlh.second->AddSubMatrix(dofs_lor, dofs_ho, M_LH_el);
-         r_and_mlh.first->AddSubMatrix(dofs_lor, dofs_ho, R_el);
+         if (build_R)
+         {
+            r_and_mlh.first->AddSubMatrix(dofs_lor, dofs_ho, R_el);
+         }
 
       }
    }
@@ -2009,6 +2131,8 @@ const Operator &L2ProjectionGridTransfer::ForwardOperator()
 
 const Operator &L2ProjectionGridTransfer::BackwardOperator()
 {
+   MFEM_VERIFY(!UsesH1ConsistentMass(),
+               "BackwardOperator is not supported with consistent mass");
    if (!B)
    {
       if (!F) { BuildF(); }
@@ -2017,15 +2141,30 @@ const Operator &L2ProjectionGridTransfer::BackwardOperator()
    return *B;
 }
 
+void L2ProjectionGridTransfer::UseConsistentMass(bool use_consistent_mass_)
+{
+   MFEM_VERIFY(!F && !B,
+               "UseConsistentMass must be called before constructing operators");
+   use_consistent_mass = use_consistent_mass_;
+}
+
+bool L2ProjectionGridTransfer::UsesH1ConsistentMass() const
+{
+   return use_consistent_mass && !force_l2_space &&
+          dom_fes.FEColl()->GetContType() == FiniteElementCollection::CONTINUOUS;
+}
+
 void L2ProjectionGridTransfer::BuildF()
 {
    if (!force_l2_space &&
        dom_fes.FEColl()->GetContType() == FiniteElementCollection::CONTINUOUS)
    {
+      MFEM_VERIFY(!(use_ea && use_consistent_mass),
+                  "consistent mass is not supported with element assembly");
       if (!Parallel())
       {
          F = new L2ProjectionH1Space(dom_fes, ran_fes,
-                                     use_ea, d_mt);
+                                     use_ea, use_consistent_mass, d_mt);
       }
       else
       {
@@ -2035,7 +2174,7 @@ void L2ProjectionGridTransfer::BuildF()
          const mfem::ParFiniteElementSpace& ran_pfes =
             static_cast<mfem::ParFiniteElementSpace&>(ran_fes);
          F = new L2ProjectionH1Space(dom_pfes, ran_pfes,
-                                     use_ea, d_mt);
+                                     use_ea, use_consistent_mass, d_mt);
 #endif
       }
    }
@@ -2048,6 +2187,7 @@ void L2ProjectionGridTransfer::BuildF()
 
 bool L2ProjectionGridTransfer::SupportsBackwardsOperator() const
 {
+   if (UsesH1ConsistentMass()) { return false; }
    return ran_fes.GetTrueVSize() >= dom_fes.GetTrueVSize();
 }
 

fem/transfer.hpp

@@ -169,10 +169,12 @@ public:
     is the forward transfer matrix, and M_f is the mass matrix on the coarse
     element. For L2 spaces, M_f is the mass matrix on the union of all fine
     elements comprising the coarse element. For H1 spaces, M_f is a diagonal
-    (lumped) mass matrix computed through row-summation. Note that the backward
-    transfer operator, B, is a left inverse of the forward transfer operator, F,
-    i.e. B F = I. Both F and B are defined in physical space and, generally for
-    L2 spaces, vary between different mesh elements.
+    (lumped) mass matrix computed through row-summation, unless
+    UseConsistentMass() is enabled for the forward H1 operator. When the
+    backward transfer operator, B, is supported, it is a left inverse of the
+    forward transfer operator, F, i.e. B F = I. Both F and B are defined in
+    physical space and, generally for L2 spaces, vary between different mesh
+    elements.
 
     This class supports H1 and L2 finite element spaces. Fine meshes are a
     uniform refinement of the coarse mesh, usually created through
@@ -352,16 +354,21 @@ public:
    class L2ProjectionH1Space : public L2Projection
    {
       const bool use_ea;
+      /// Use the consistent low-order mass matrix in non-EA H1 Mult() and
+      /// MultTranspose().
+      const bool use_consistent_mass;
 
    public:
       L2ProjectionH1Space(const FiniteElementSpace &fes_ho_,
                           const FiniteElementSpace &fes_lor_,
                           const bool use_ea_,
+                          const bool use_consistent_mass_,
                           MemoryType d_mt_ = Device::GetHostMemoryType());
 #ifdef MFEM_USE_MPI
       L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_,
                           const ParFiniteElementSpace &pfes_lor_,
                           const bool use_ea_,
+                          const bool use_consistent_mass_,
                           MemoryType d_mt_ = Device::GetHostMemoryType());
 #endif
       /// Same as above but assembles action of R through 4 parts:
@@ -417,13 +424,33 @@ public:
       void SetAbsTol(real_t p_atol_) override;
 
    protected:
+      /// Applies the H1 transfer R = M_L^{-1} M_LH and its transpose, where
+      /// M_L is the consistent low-order mass matrix.
+      class H1ConsistentMassOperator : public Operator
+      {
+      private:
+         const Operator &M_LH;
+         const Solver &M_L_solver;
+
+      public:
+         H1ConsistentMassOperator(const Operator &M_LH_,
+                                  const Solver &M_L_solver_);
+
+         void Mult(const Vector &x, Vector &y) const override;
+
+         void MultTranspose(const Vector &x, Vector &y) const override;
+      };
+
       /// Sets up the PCG solver (sets parameters, operator, and preconditioner)
       void SetupPCG();
 
-      /// @brief Computes on-rank R and M_LH matrices. If true, computes mixed mass and/or
-      /// inverse lumped mass matrix error when compared to device implementation.
+      /** @brief Computes on-rank R and M_LH matrices.
+
+          If build_R is true, the returned pair contains both R and M_LH. If
+          build_R is false, the first pointer is null and only M_LH is built. */
       std::pair<std::unique_ptr<SparseMatrix>,
-          std::unique_ptr<SparseMatrix>> ComputeSparseRAndM_LH();
+          std::unique_ptr<SparseMatrix>> ComputeSparseRAndM_LH(
+             bool build_R = true);
 
       /// @brief Recovers vector of tdofs given a vector of dofs and a finite
       /// element space
@@ -453,20 +480,30 @@ public:
       /// elements and refined LOR elements.
       std::unique_ptr<SparseMatrix> AllocR();
 
-      CGSolver pcg;
-      std::unique_ptr<Solver> precon;
+      /// Consistent low-order mass matrix used when use_consistent_mass is true.
+      std::unique_ptr<Operator> M_L;
+      // Used to compute P = (RT*M_LH)^(-1) M_LH^T
+      std::unique_ptr<Operator> M_LH;
+      // Lumped M_L inverse operator built via EA. Wrapped with restriction maps
+      // to multiply with scalar TDof LOR vectors.
+      std::unique_ptr<Operator> ML_inv_vea;
+      /// Preconditioner for applying the inverse consistent low-order mass
+      /// matrix.
+      std::unique_ptr<Solver> ML_precon;
+      /// Serial PCG solver for applying the inverse consistent low-order mass
+      /// matrix in H1 Mult() and MultTranspose().
+      CGSolver ML_pcg;
+      /// Solver used by H1ConsistentMassOperator to apply M_L^{-1}.
+      std::unique_ptr<Solver> ML_solver;
       // 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.
       // If element assembly is enabled
       std::unique_ptr<Operator> R;
-      // Used to compute P = (RT*M_LH)^(-1) M_LH^T
-      std::unique_ptr<Operator> M_LH;
       // Inverted operator in P = (RT*M_LH)^(-1) M_LH^T. Used to compute P via PCG.
       std::unique_ptr<Operator> RTxM_LH;
-      // Lumped M_L inverse operator built via EA. Wrapped with restriction maps
-      // to multiply with scalar TDof LOR vectors.
-      std::unique_ptr<Operator> ML_inv_vea;
+      std::unique_ptr<Solver> precon;
+      CGSolver pcg;
       // LDof Mixed mass operator built via EA. Wrapped with restriction maps to send
       // scalar LDof HO vectors to LDof LOR vectors.
       Operator *M_LH_local_op;
@@ -478,7 +515,6 @@ public:
       Vector M_LH_ea;
       // Element Assembled lumped M_L inverse built via EA. Stores diagonal as a Ldof vector.
       Vector ML_inv_ea;
-
 #ifdef MFEM_USE_MPI
       std::unique_ptr<ParFiniteElementSpace> pfes_ho_scalar;
       std::unique_ptr<ParFiniteElementSpace> pfes_lor_scalar;
@@ -511,6 +547,9 @@ public:
    L2Projection   *F; ///< Forward, coarse-to-fine, operator
    L2Prolongation *B; ///< Backward, fine-to-coarse, operator
    bool force_l2_space;
+   /// Use the consistent low-order mass matrix for non-EA H1 Mult() and
+   /// MultTranspose().
+   bool use_consistent_mass;
 
 public:
    L2ProjectionGridTransfer(FiniteElementSpace &coarse_fes_,
@@ -518,16 +557,26 @@ public:
                             bool force_l2_space_ = false,
                             MemoryType d_mt_ = Device::GetHostMemoryType()) // move to method
       : GridTransfer(coarse_fes_, fine_fes_),
-        F(NULL), B(NULL), force_l2_space(force_l2_space_)
+        F(NULL), B(NULL), force_l2_space(force_l2_space_),
+        use_consistent_mass(false)
    { }
    virtual ~L2ProjectionGridTransfer();
 
+   /** @brief Use the consistent low-order mass matrix in H1 non-EA Mult() and
+       MultTranspose().
+
+       This option must be set before constructing the transfer operators. It
+       only affects H1 transfer, is not supported with element assembly, and
+       disables BackwardOperator(). */
+   void UseConsistentMass(bool use_consistent_mass_ = true);
+
    const Operator &ForwardOperator() override;
 
    const Operator &BackwardOperator() override;
 
    bool SupportsBackwardsOperator() const override;
 private:
+   bool UsesH1ConsistentMass() const;
    void BuildF();
 };
 

tests/unit/fem/test_transfer.cpp

@@ -440,6 +440,99 @@ TEST_CASE("Variable Order True Transfer", "[Transfer][VariableOrder]")
    delete c_fec;
 }
 
+TEST_CASE("H1 L2 transfer with consistent mass", "[Transfer]")
+{
+   auto vectorspace = GENERATE(VecSpace::H1, VecSpace::VectorH1nodes,
+                               VecSpace::VectorH1vdim);
+   dimension = GENERATE(2, 3);
+
+   const int order = 2;
+   const int ne = 2;
+   const int vdim = (vectorspace == VecSpace::VectorH1nodes
+                     || vectorspace == VecSpace::VectorH1vdim) ? dimension : 1;
+   Ordering::Type ordering = (vectorspace == VecSpace::VectorH1vdim)
+                             ? Ordering::byVDIM : Ordering::byNODES;
+
+   CAPTURE(VecSpaceName(vectorspace), dimension, order);
+
+   Mesh mesh;
+   if (dimension == 2)
+   {
+      mesh = Mesh::MakeCartesian2D(ne, ne, Element::QUADRILATERAL,
+                                   1, 1.0, 1.0);
+   }
+   else
+   {
+      mesh = Mesh::MakeCartesian3D(ne, ne, ne, Element::HEXAHEDRON,
+                                   1.0, 1.0, 1.0);
+   }
+
+   Mesh fineMesh(mesh);
+   fineMesh.UniformRefinement();
+
+   H1_FECollection fec(order, dimension);
+   FiniteElementSpace c_fespace(&mesh, &fec, vdim, ordering);
+   FiniteElementSpace f_fespace(&fineMesh, &fec, vdim, ordering);
+
+   L2ProjectionGridTransfer transfer(c_fespace, f_fespace);
+   transfer.UseConsistentMass();
+   const Operator &R = transfer.ForwardOperator();
+
+   GridFunction X(&c_fespace);
+   GridFunction Y(&f_fespace);
+   GridFunction Y_ref(&f_fespace);
+   coeff_order = 1;
+
+   LinearForm rhs(&f_fespace);
+   BilinearForm mass(&f_fespace);
+   FunctionCoefficient funcCoeff(&coeff);
+   VectorFunctionCoefficient vecCoeff(dimension, &vectorcoeff);
+   if (vectorspace == VecSpace::H1)
+   {
+      X.ProjectCoefficient(funcCoeff);
+      rhs.AddDomainIntegrator(new DomainLFIntegrator(funcCoeff));
+      mass.AddDomainIntegrator(new MassIntegrator);
+   }
+   else
+   {
+      X.ProjectCoefficient(vecCoeff);
+      rhs.AddDomainIntegrator(new VectorDomainLFIntegrator(vecCoeff));
+      mass.AddDomainIntegrator(new VectorMassIntegrator);
+   }
+
+   rhs.Assemble();
+   mass.Assemble();
+   SparseMatrix M;
+   Array<int> empty;
+   mass.FormSystemMatrix(empty, M);
+
+   GSSmoother M_prec(M);
+   Y_ref = 0.0;
+   PCG(M, M_prec, rhs, Y_ref, 0, 500, 1e-24, 0.0);
+
+   Y = 0.0;
+   R.Mult(X, Y);
+   Y -= Y_ref;
+   REQUIRE(Y.Norml2() < 1e-11 * Y_ref.Norml2());
+
+   Vector x(c_fespace.GetVSize());
+   Vector y(f_fespace.GetVSize());
+   Vector Ry(f_fespace.GetVSize());
+   Vector Rtx(c_fespace.GetVSize());
+   x.Randomize(1);
+   y.Randomize(2);
+
+   R.Mult(x, Ry);
+   R.MultTranspose(y, Rtx);
+
+   const real_t ip1 = InnerProduct(Ry, y);
+   const real_t ip2 = InnerProduct(x, Rtx);
+   REQUIRE(std::abs(ip1 - ip2) <
+           1e-10 * std::max(std::abs(ip1), std::abs(ip2)));
+
+   REQUIRE_FALSE(transfer.SupportsBackwardsOperator());
+}
+
 TEST_CASE("Restriction Transpose Operator")
 {
    int order = GENERATE(1, 2);
@@ -797,6 +890,59 @@ TEST_CASE("Parallel Transfer", "[Transfer][Parallel]")
    delete pmesh;
 }
 
+TEST_CASE("Parallel H1 L2 transfer with consistent mass",
+          "[Transfer][Parallel]")
+{
+   dimension = GENERATE(2, 3);
+
+   const int order = 2;
+   const int ne = 2;
+   const int vdim = 1;
+
+   CAPTURE(dimension, order);
+
+   Mesh mesh;
+   if (dimension == 2)
+   {
+      mesh = Mesh::MakeCartesian2D(ne, ne, Element::QUADRILATERAL,
+                                   1, 1.0, 1.0);
+   }
+   else
+   {
+      mesh = Mesh::MakeCartesian3D(ne, ne, ne, Element::HEXAHEDRON,
+                                   1.0, 1.0, 1.0);
+   }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   ParMesh pfineMesh(MPI_COMM_WORLD, mesh);
+   pfineMesh.UniformRefinement();
+
+   H1_FECollection fec(order, dimension);
+   ParFiniteElementSpace c_fespace(&pmesh, &fec, vdim);
+   ParFiniteElementSpace f_fespace(&pfineMesh, &fec, vdim);
+
+   L2ProjectionGridTransfer transfer(c_fespace, f_fespace);
+   transfer.UseConsistentMass();
+   const Operator &R = transfer.TrueForwardOperator();
+
+   Vector x(c_fespace.GetTrueVSize());
+   Vector y(f_fespace.GetTrueVSize());
+   Vector Rx(f_fespace.GetTrueVSize());
+   Vector Rty(c_fespace.GetTrueVSize());
+   x.Randomize(1);
+   y.Randomize(2);
+
+   R.Mult(x, Rx);
+   R.MultTranspose(y, Rty);
+
+   const real_t ip1 = InnerProduct(MPI_COMM_WORLD, Rx, y);
+   const real_t ip2 = InnerProduct(MPI_COMM_WORLD, x, Rty);
+   REQUIRE(std::abs(ip1 - ip2) <
+           1e-10 * std::max(std::abs(ip1), std::abs(ip2)));
+
+   REQUIRE_FALSE(transfer.SupportsBackwardsOperator());
+}
+
 TEST_CASE("Trace PRefinement Parallel TrueTransfer", "[Transfer][Parallel]")
 {
    auto simplex = GENERATE(true, false);