Print the parts of the dof_truedof matrix

Change default command line argument values
Use different manufactured solution
2026-01-12 18:50:34 -08:00 · 2025-09-08 18:13:21 -07:00 · 2025-09-04 20:29:30 -07:00 · 2025-06-15 19:42:36 -07:00 · 2024-07-28 18:41:55 -07:00 · 2024-05-27 21:22:11 -07:00
20 changed files with 1757 additions and 253 deletions
@@ -103,6 +103,7 @@ New and updated examples and miniapps
  saddle-point solver methodology, suitable for high-order discretizations and
  for GPU acceleration. Examples illustrating the solution of Darcy and grad-div
  problems are included. See miniapps/hdiv-linear-solver.
+>>>>>>> master

 - Added new Discontinuous Petrov-Galerkin (DPG) miniapp which includes serial
  and parallel examples for diffusion, convection-diffusion, acoustics and
@@ -394,8 +394,8 @@ void BilinearForm::Assemble(int skip_zeros)
   }

   ElementTransformation *eltrans;
-   DofTransformation * doftrans;
-   Mesh *mesh = fes -> GetMesh();
+   DofTransformation *doftrans;
+   Mesh *mesh = fes->GetMesh();
   DenseMatrix elmat, *elmat_p;

   if (mat == NULL)
@@ -432,8 +432,9 @@ void BilinearForm::Assemble(int skip_zeros)
      }

      // Element-wise integration
-      for (int i = 0; i < fes -> GetNE(); i++)
+      for (int i = 0; i < fes->GetNE(); i++)
      {
+         doftrans = fes->GetElementVDofs(i, vdofs);
         if (element_matrices)
         {
            elmat_p = &(*element_matrices)(i);
@@ -441,9 +442,7 @@ void BilinearForm::Assemble(int skip_zeros)
         else
         {
            const int elem_attr = fes->GetMesh()->GetAttribute(i);
-            doftrans = fes->GetElementVDofs(i, vdofs);
            eltrans = fes->GetElementTransformation(i);
-
            elmat.SetSize(0);
            for (int k = 0; k < domain_integs.Size(); k++)
            {
@@ -494,7 +493,7 @@ void BilinearForm::Assemble(int skip_zeros)
      // Patch-wise integration
      if (fes->GetNURBSext())
      {
-         for (int p=0; p<mesh->NURBSext->GetNP(); ++p)
+         for (int p = 0; p < mesh->NURBSext->GetNP(); ++p)
         {
            bool vdofsSet = false;
            for (int k = 0; k < domain_integs.Size(); k++)
@@ -507,15 +506,15 @@ void BilinearForm::Assemble(int skip_zeros)
                     vdofsSet = true;
                  }

-                  SparseMatrix* spmat = nullptr;
+                  SparseMatrix *spmat = nullptr;
                  domain_integs[k]->AssemblePatchMatrix(p, *fes, spmat);
                  Array<int> cols;
                  Vector srow;

-                  for (int r=0; r<spmat->Height(); ++r)
+                  for (int r = 0; r < spmat->Height(); ++r)
                  {
                     spmat->GetRow(r, cols, srow);
-                     for (int i=0; i<cols.Size(); ++i)
+                     for (int i = 0; i < cols.Size(); ++i)
                     {
                        cols[i] = vdofs[cols[i]];
                     }
@@ -552,14 +551,14 @@ void BilinearForm::Assemble(int skip_zeros)
         }
      }

-      for (int i = 0; i < fes -> GetNBE(); i++)
+      for (int i = 0; i < fes->GetNBE(); i++)
      {
         const int bdr_attr = mesh->GetBdrAttribute(i);
         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }

         const FiniteElement &be = *fes->GetBE(i);
-         doftrans = fes -> GetBdrElementVDofs (i, vdofs);
-         eltrans = fes -> GetBdrElementTransformation (i);
+         doftrans = fes->GetBdrElementVDofs(i, vdofs);
+         eltrans = fes->GetBdrElementTransformation(i);
         int k = 0;
         for (; k < boundary_integs.Size(); k++)
         {
@@ -606,11 +605,11 @@ void BilinearForm::Assemble(int skip_zeros)
      int nfaces = mesh->GetNumFaces();
      for (int i = 0; i < nfaces; i++)
      {
-         tr = mesh -> GetInteriorFaceTransformations (i);
+         tr = mesh->GetInteriorFaceTransformations(i);
         if (tr != NULL)
         {
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
-            fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
+            fes->GetElementVDofs(tr->Elem1No, vdofs);
+            fes->GetElementVDofs(tr->Elem2No, vdofs2);
            vdofs.Append (vdofs2);
            for (int k = 0; k < interior_face_integs.Size(); k++)
            {
@@ -618,7 +617,7 @@ void BilinearForm::Assemble(int skip_zeros)
               AssembleFaceMatrix(*fes->GetFE(tr->Elem1No),
                                  *fes->GetFE(tr->Elem2No),
                                  *tr, elemmat);
-               mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
+               mat->AddSubMatrix(vdofs, vdofs, elemmat, skip_zeros);
            }
         }
      }
@@ -650,16 +649,16 @@ void BilinearForm::Assemble(int skip_zeros)
         }
      }

-      for (int i = 0; i < fes -> GetNBE(); i++)
+      for (int i = 0; i < fes->GetNBE(); i++)
      {
         const int bdr_attr = mesh->GetBdrAttribute(i);
         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }

-         tr = mesh -> GetBdrFaceTransformations (i);
+         tr = mesh->GetBdrFaceTransformations(i);
         if (tr != NULL)
         {
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
-            fe1 = fes -> GetFE (tr -> Elem1No);
+            fes->GetElementVDofs(tr->Elem1No, vdofs);
+            fe1 = fes->GetFE(tr->Elem1No);
            // The fe2 object is really a dummy and not used on the boundaries,
            // but we can't dereference a NULL pointer, and we don't want to
            // actually make a fake element.
@@ -670,9 +669,9 @@ void BilinearForm::Assemble(int skip_zeros)
                   (*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
               { continue; }

-               boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
-                                                              elemmat);
-               mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
+               boundary_face_integs[k]->AssembleFaceMatrix(*fe1, *fe2, *tr,
+                                                           elemmat);
+               mat->AddSubMatrix(vdofs, vdofs, elemmat, skip_zeros);
            }
         }
      }
@@ -1206,6 +1205,7 @@ MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
   mat = NULL;
   mat_e = NULL;
   extern_bfs = 0;
+   element_matrices = NULL;
   assembly = AssemblyLevel::LEGACY;
   ext = NULL;
 }
@@ -1220,6 +1220,7 @@ MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
   mat = NULL;
   mat_e = NULL;
   extern_bfs = 1;
+   element_matrices = NULL;
   ext = NULL;

   // Copy the pointers to the integrators
@@ -1405,17 +1406,26 @@ void MixedBilinearForm::Assemble(int skip_zeros)
   }

   ElementTransformation *eltrans;
-   DofTransformation * dom_dof_trans;
-   DofTransformation * ran_dof_trans;
-   DenseMatrix elmat;
+   DofTransformation *dom_dof_trans;
+   DofTransformation *ran_dof_trans;
+   DenseMatrix elmat, *elmat_p;

-   Mesh *mesh = test_fes -> GetMesh();
+   Mesh *mesh = test_fes->GetMesh();

   if (mat == NULL)
   {
      mat = new SparseMatrix(height, width);
   }

+#ifdef MFEM_USE_LEGACY_OPENMP
+   int free_element_matrices = 0;
+   if (!element_matrices)
+   {
+      ComputeElementMatrices();
+      free_element_matrices = 1;
+   }
+#endif
+
   if (domain_integs.Size())
   {
      for (int k = 0; k < domain_integs.Size(); k++)
@@ -1429,31 +1439,38 @@ void MixedBilinearForm::Assemble(int skip_zeros)
         }
      }

-      for (int i = 0; i < test_fes -> GetNE(); i++)
+      for (int i = 0; i < test_fes->GetNE(); i++)
      {
         const int elem_attr = mesh->GetAttribute(i);
-         dom_dof_trans = trial_fes -> GetElementVDofs (i, trial_vdofs);
-         ran_dof_trans = test_fes  -> GetElementVDofs (i, test_vdofs);
-         eltrans = test_fes -> GetElementTransformation (i);
-
-         elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
-         elmat = 0.0;
-         for (int k = 0; k < domain_integs.Size(); k++)
+         dom_dof_trans = trial_fes->GetElementVDofs(i, trial_vdofs);
+         ran_dof_trans = test_fes->GetElementVDofs(i, test_vdofs);
+         if (element_matrices)
         {
-            if (domain_integs_marker[k] == NULL ||
-                (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+            elmat_p = &(*element_matrices)(i);
+         }
+         else
+         {
+            eltrans = test_fes->GetElementTransformation(i);
+            elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
+            elmat = 0.0;
+            for (int k = 0; k < domain_integs.Size(); k++)
            {
-               domain_integs[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
-                                                           *test_fes  -> GetFE(i),
+               if (domain_integs_marker[k] == NULL ||
+                   (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+               {
+                  domain_integs[k]->AssembleElementMatrix2(*trial_fes->GetFE(i),
+                                                           *test_fes->GetFE(i),
                                                           *eltrans, elemmat);
-               elmat += elemmat;
+                  elmat += elemmat;
+               }
            }
+            if (ran_dof_trans || dom_dof_trans)
+            {
+               TransformDual(ran_dof_trans, dom_dof_trans, elmat);
+            }
+            elmat_p = &elmat;
         }
-         if (ran_dof_trans || dom_dof_trans)
-         {
-            TransformDual(ran_dof_trans, dom_dof_trans, elmat);
-         }
-         mat -> AddSubMatrix (test_vdofs, trial_vdofs, elmat, skip_zeros);
+         mat->AddSubMatrix(test_vdofs, trial_vdofs, *elmat_p, skip_zeros);
      }
   }

@@ -1480,14 +1497,14 @@ void MixedBilinearForm::Assemble(int skip_zeros)
         }
      }

-      for (int i = 0; i < test_fes -> GetNBE(); i++)
+      for (int i = 0; i < test_fes->GetNBE(); i++)
      {
         const int bdr_attr = mesh->GetBdrAttribute(i);
         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }

-         dom_dof_trans = trial_fes -> GetBdrElementVDofs (i, trial_vdofs);
-         ran_dof_trans = test_fes  -> GetBdrElementVDofs (i, test_vdofs);
-         eltrans = test_fes -> GetBdrElementTransformation (i);
+         dom_dof_trans = trial_fes->GetBdrElementVDofs(i, trial_vdofs);
+         ran_dof_trans = test_fes->GetBdrElementVDofs(i, test_vdofs);
+         eltrans = test_fes->GetBdrElementTransformation(i);

         elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
         elmat = 0.0;
@@ -1496,16 +1513,17 @@ void MixedBilinearForm::Assemble(int skip_zeros)
            if (boundary_integs_marker[k] &&
                (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }

-            boundary_integs[k]->AssembleElementMatrix2 (*trial_fes -> GetBE(i),
-                                                        *test_fes  -> GetBE(i),
-                                                        *eltrans, elemmat);
+            boundary_integs[k]->AssembleElementMatrix2(*trial_fes->GetBE(i),
+                                                       *test_fes->GetBE(i),
+                                                       *eltrans, elemmat);
            elmat += elemmat;
         }
         if (ran_dof_trans || dom_dof_trans)
         {
            TransformDual(ran_dof_trans, dom_dof_trans, elmat);
         }
-         mat -> AddSubMatrix (test_vdofs, trial_vdofs, elmat, skip_zeros);
+         elmat_p = &elmat;
+         mat->AddSubMatrix(test_vdofs, trial_vdofs, *elmat_p, skip_zeros);
      }
   }

@@ -1603,6 +1621,13 @@ void MixedBilinearForm::Assemble(int skip_zeros)
         }
      }
   }
+
+#ifdef MFEM_USE_LEGACY_OPENMP
+   if (free_element_matrices)
+   {
+      FreeElementMatrices();
+   }
+#endif
 }

 void MixedBilinearForm::AssembleDiagonal_ADAt(const Vector &D,
@@ -1688,6 +1713,13 @@ void MixedBilinearForm::ConformingAssemble()

 void MixedBilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
 {
+   if (element_matrices)
+   {
+      elmat.SetSize(element_matrices->SizeI(), element_matrices->SizeJ());
+      elmat = element_matrices->GetData(i);
+      return;
+   }
+
   if (domain_integs.Size())
   {
      const FiniteElement &trial_fe = *trial_fes->GetFE(i);
@@ -1774,6 +1806,50 @@ void MixedBilinearForm::AssembleBdrElementMatrix(
   mat->AddSubMatrix(test_vdofs_, trial_vdofs_, elmat, skip_zeros);
 }

+void MixedBilinearForm::ComputeElementMatrices()
+{
+   if (element_matrices || domain_integs.Size() == 0 || trial_fes->GetNE() == 0)
+   {
+      return;
+   }
+
+   int num_elements = trial_fes->GetNE();
+   int trial_dofs_per_el = trial_fes->GetFE(0)->GetDof() * trial_fes->GetVDim();
+   int test_dofs_per_el = test_fes->GetFE(0)->GetDof() * test_fes->GetVDim();
+
+   element_matrices = new DenseTensor(test_dofs_per_el, trial_dofs_per_el,
+                                      num_elements);
+   DenseMatrix tmp;
+   IsoparametricTransformation eltrans;
+#ifdef MFEM_USE_LEGACY_OPENMP
+   #pragma omp parallel for private(tmp,eltrans)
+#endif
+   for (int i = 0; i < num_elements; i++)
+   {
+      DenseMatrix elmat(element_matrices->GetData(i),
+                        test_dofs_per_el, trial_dofs_per_el);
+      const FiniteElement &trial_fe = *trial_fes->GetFE(i);
+      const FiniteElement &test_fe = *test_fes->GetFE(i);
+#ifdef MFEM_DEBUG
+      if (trial_dofs_per_el != trial_fe.GetDof() * trial_fes->GetVDim())
+         mfem_error("MixedBilinearForm::ComputeElementMatrices:"
+                    " all elements must have same number of dofs");
+#endif
+      test_fes->GetElementTransformation(i, &eltrans);
+
+      domain_integs[0]->AssembleElementMatrix2(trial_fe, test_fe, eltrans,
+                                               elmat);
+      for (int k = 1; k < domain_integs.Size(); k++)
+      {
+         // note: some integrators may not be thread-safe
+         domain_integs[k]->AssembleElementMatrix2(trial_fe, test_fe, eltrans,
+                                                  tmp);
+         elmat += tmp;
+      }
+      elmat.ClearExternalData();
+   }
+}
+
 void MixedBilinearForm::EliminateTrialDofs (
   const Array<int> &bdr_attr_is_ess, const Vector &sol, Vector &rhs )
 {
@@ -739,6 +739,7 @@ protected:
   Array<Array<int>*> boundary_trace_face_integs_marker;

   DenseMatrix elemmat;
+   DenseTensor *element_matrices; ///< Owned.
   Array<int>  trial_vdofs, test_vdofs;

 private:
@@ -895,7 +896,18 @@ public:
       MixedBilinearForm becomes an operator on the conforming FE spaces. */
   void ConformingAssemble();

+   /// Compute and store internally all element matrices.
+   void ComputeElementMatrices();
+
+   /// Free the memory used by the element matrices.
+   void FreeElementMatrices()
+   { delete element_matrices; element_matrices = NULL; }
+
   /// Compute the element matrix of the given element
+   /** The element matrix is computed by calling the domain integrators
+       or the one stored internally by a prior call of ComputeElementMatrices()
+       is returned when available.
+   */
   void ComputeElementMatrix(int i, DenseMatrix &elmat);

   /// Compute the boundary element matrix of the given boundary element
@@ -120,6 +120,10 @@ public:
   /// Return the serial hybridized matrix.
   SparseMatrix &GetMatrix() { return *H; }

+   /// Return the transpose of the serial constraint matrix.
+   SparseMatrix &GetConstraintMatrixTranspose()
+   { if (!Ct) mfem_error("Ct has not been constructed!"); return *Ct; }
+
 #ifdef MFEM_USE_MPI
   /// Return the parallel hybridized matrix.
   HypreParMatrix &GetParallelMatrix() { return *pH.Is<HypreParMatrix>(); }
@@ -266,6 +266,9 @@ public:
      test_pfes  = test_fes;
   }

+   /// Return the test parallel FE space associated with the ParMixedBilinearForm.
+   ParFiniteElementSpace *TestParFESpace() const { return test_pfes; }
+
   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
   HypreParMatrix *ParallelAssemble();

@@ -744,6 +744,8 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
      // Safe to assume 1-1 correspondence between shared dofs
      int ldof  = GetVSize();
      int ltdof = TrueVSize();
+      std::cout << MyRank << " (ldof, ltdof) = ("
+                << ldof << "," << ltdof << ")" << std::endl;

      HYPRE_Int *i_diag = Memory<HYPRE_Int>(ldof+1);
      HYPRE_Int *j_diag = Memory<HYPRE_Int>(ltdof);
@@ -756,7 +758,11 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
      HYPRE_BigInt *cmap   = Memory<HYPRE_BigInt>(ldof-ltdof);

      HYPRE_BigInt *col_starts = GetTrueDofOffsets();
+      std::cout << MyRank << " col_starts: [" << col_starts[0]
+                << "," << col_starts[1] << "]" << std::endl;
      HYPRE_BigInt *row_starts = GetDofOffsets();
+      std::cout << MyRank << " row_starts: [" << row_starts[0]
+                << "," << row_starts[1] << "]" << std::endl;

      Array<Pair<HYPRE_BigInt, int> > cmap_j_offd(ldof-ltdof);

@@ -778,6 +784,18 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
         i_diag[i+1] = diag_counter;
         i_offd[i+1] = offd_counter;
      }
+      std::cout << MyRank << " i_diag: ";
+      for (int i = 0; i < ldof+1; ++i)
+      {
+         std::cout << i_diag[i] << " ";
+      }
+      std::cout << std::endl;
+      std::cout << MyRank << " j_diag: ";
+      for (int i = 0; i < ltdof; ++i)
+      {
+         std::cout << j_diag[i] << " ";
+      }
+      std::cout << std::endl;

      SortPairs<HYPRE_BigInt, int>(cmap_j_offd, offd_counter);

@@ -786,6 +804,24 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
         cmap[i] = cmap_j_offd[i].one;
         j_offd[cmap_j_offd[i].two] = i;
      }
+      std::cout << MyRank << " i_offd: ";
+      for (int i = 0; i < ldof+1; ++i)
+      {
+         std::cout << i_offd[i] << " ";
+      }
+      std::cout << std::endl;
+      std::cout << MyRank << " j_offd: ";
+      for (int i = 0; i < ldof-ltdof; ++i)
+      {
+         std::cout << j_offd[i] << " ";
+      }
+      std::cout << std::endl;
+      std::cout << MyRank << " cmap: ";
+      for (int i = 0; i < ldof-ltdof; ++i)
+      {
+         std::cout << cmap[i] << " ";
+      }
+      std::cout << std::endl;

      P = new HypreParMatrix(MyComm, MyRank, NRanks, row_starts, col_starts,
                             i_diag, j_diag, i_offd, j_offd,
@@ -365,6 +365,44 @@ void SecondOrderTimeDependentOperator::ImplicitSolve(const double dt0,
   mfem_error("SecondOrderTimeDependentOperator::ImplicitSolve() is not overridden!");
 }

+SumOperator::SumOperator(const Operator *A, const double alpha,
+                         const Operator *B, const double beta,
+                         bool ownA, bool ownB)
+   : Operator(A->Height(), A->Width()),
+     A(A), B(B), alpha(alpha), beta(beta), ownA(ownA), ownB(ownB),
+     a(A->Width()), b(B->Width())
+{
+   MFEM_VERIFY(A->Width() == B->Width(),
+               "incompatible Operators: different widths\n"
+               << "A->Width() = " << A->Width()
+               << ", B->Width() = " << B->Width() );
+   MFEM_VERIFY(A->Height() == B->Height(),
+               "incompatible Operators: different heights\n"
+               << "A->Height() = " << A->Height()
+               << ", B->Height() = " << B->Height() );
+
+   {
+      const Solver* SolverA = dynamic_cast<const Solver*>(A);
+      const Solver* SolverB = dynamic_cast<const Solver*>(B);
+      if (SolverA)
+      {
+         MFEM_VERIFY(!(SolverA->iterative_mode),
+                     "Operator A of a SumOperator should not be in iterative mode");
+      }
+      if (SolverB)
+      {
+         MFEM_VERIFY(!(SolverB->iterative_mode),
+                     "Operator B of a SumOperator should not be in iterative mode");
+      }
+   }
+
+}
+
+SumOperator::~SumOperator()
+{
+   if (ownA) { delete A; }
+   if (ownB) { delete B; }
+}

 ProductOperator::ProductOperator(const Operator *A, const Operator *B,
                                 bool ownA, bool ownB)
@@ -769,6 +769,28 @@ public:
   { A.Mult(x, y); }
 };

+/// General linear combination operator: x -> a A(x) + b B(x).
+class SumOperator : public Operator
+{
+   const Operator *A, *B;
+   const double alpha, beta;
+   bool ownA, ownB;
+   mutable Vector a, b;
+
+public:
+   SumOperator(
+      const Operator *A, const double alpha,
+      const Operator *B, const double beta,
+      bool ownA, bool ownB);
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   { A->Mult(x, a); B->Mult(x, b); add(alpha, a, beta, b, y); }
+
+   virtual void MultTranspose(const Vector &x, Vector &y) const
+   { A->MultTranspose(x, a); B->MultTranspose(x, b); add(alpha, a, beta, b, y); }
+
+   virtual ~SumOperator();
+};

 /// General product operator: x -> (A*B)(x) = A(B(x)).
 class ProductOperator : public Operator
@@ -12,8 +12,10 @@
 if (MFEM_USE_MPI)
  add_mfem_miniapp(block-solvers
    MAIN block-solvers.cpp
-    EXTRA_SOURCES div_free_solver.cpp
-    EXTRA_HEADERS div_free_solver.hpp
+    EXTRA_SOURCES darcy_solver.cpp div_free_solver.cpp bramble_pasciak.cpp
+     block_hybridization.cpp
+    EXTRA_HEADERS darcy_solver.hpp div_free_solver.hpp bramble_pasciak.hpp
+     block_hybridization.hpp
    LIBRARIES mfem)

  add_mfem_miniapp(plor_solvers
@@ -26,16 +28,16 @@ if (MFEM_USE_MPI)
    add_test(NAME block-solvers-constant_np${MFEM_MPI_NP}
      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:block-solvers> -r 2 -no-vis
+      $<TARGET_FILE:block-solvers> -sr 1 -pr 1 -no-vis
      ${MPIEXEC_POSTFLAGS})

    add_test(NAME block-solvers-anisotropic_np${MFEM_MPI_NP}
      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:block-solvers> -r 2
+      $<TARGET_FILE:block-solvers> -sr 1 -pr 1
      -m ${CMAKE_CURRENT_SOURCE_DIR}/anisotropic.mesh
      -c ${CMAKE_CURRENT_SOURCE_DIR}/anisotropic.coeff
-      -eb ${CMAKE_CURRENT_SOURCE_DIR}/anisotropic.brd
+      -eb ${CMAKE_CURRENT_SOURCE_DIR}/anisotropic.bdr
      -no-vis
      ${MPIEXEC_POSTFLAGS})

@@ -12,8 +12,8 @@ generated by MFEM.

 # Block Solvers

-The miniapp block-solvers compares the relative performance of some solvers for
-the discrete saddle point problems arising from the mixed finite element
+The miniapp block-solvers compares the relative performance of some solvers
+for the discrete saddle point problems arising from the mixed finite element
 discretization of the second order scalar elliptic problem
                   -div(K grad p) = r.h.s.
 The user can customize the problem settings by providing the mesh, the
@@ -21,7 +21,7 @@ coefficient K, and the essential/natural boundary assignment.

 The discrete saddle point problem has a block structure of the form

-                   [ M  B^T ] [u] = [g]
+              Ax = [ M  B^T ] [u] = [g] = b
                   [ B   0  ] [p] = [f]

 The solvers for the above block system include:
@@ -51,7 +51,28 @@ The solvers for the above block system include:
       using constrained first-order system least squares (CFOSLS).
       J. Comput. Phys. 373: 863-876, 2018.

-2. MINRES preconditioned by a block diagonal preconditioner.
+2. A block hybridization solver.
+
+   The flux u belongs to the Raviart-Thomas finite element space, where the
+   normal derivative of u is continuous across interior element faces. In
+   hybridization, the continuity requirement is broken and enforced through
+   a Lagrange multiplier instead.
+
+   A reduced system, for the multiplier only, is constructed using element
+   matrices. Its coefficient matrix H is symmetric positive definite, so
+   the conjugate gradient method is used to solve the system. In MFEM the
+   default constraint matrix C has a proper scaling that ensures the near
+   nullspace of H consists of constant vectors. Therefore, BoomerAMG is a
+   suitable preconditioner.
+
+   [3] Dobrev, Kolev, Lee, Tomov, Vassilevski, Algebraic Hybridization and Static
+       Condensation with Application to Scalable H(div) Preconditioning.
+       SIAM J. Sci. Comput. 41: B425-B447, 2019.
+   [4] Lee, Vassilevski, Parallel Solver for H(div) Problems Using Hybridization
+       and AMG, in Domain Decomposition Methods in Science and Engineering 23:
+       69-80, 2017.
+
+3. MINRES preconditioned by a block diagonal preconditioner.

   This solver is the one used in examples/ex5p.cpp.

@@ -62,6 +83,26 @@ The solvers for the above block system include:
             P^{-1} = [ diag(M)^{-1}     0    ]
                      [  0             AMG(S) ]

+3. CG with Bramble-Pasciak transformation
+
+   The solver explores two approaches. Firstly, we consider a preconditioner Q
+   such that M - Q is still s.p.d.. The transformed system XA x = X b, where
+                  X = [ M*Q^{-1} - I    0 ]
+                      [   B*Q^{-1}     -I ],
+   can be solver with a standard PGC solver, where we use the above block
+   diagonal preconditioner.
+   Secondly, we consider the particular preconditioner
+             H^{-1} = [ (M - Q)^{-1}     0     ]
+                      [     0         P_2^{-1} ].
+   This preconditioner enables a more efficient implementation of PCG what we
+   refer to as Bramble-Pasciak CG (BPCG).
+
+   For more details see
+
+   [1] Bramble, James H., and Joseph E. Pasciak. A preconditioning technique
+       for indefinite systems resulting from mixed approximations of elliptic
+       problems. Mathematics of Computation 50.181 (1988): 1-17. 
+
 # Low-Order Refined Solvers

 The miniapp `lor` (and its parallel counterpart `plor`) demonstrate the use of
@@ -26,33 +26,37 @@
 // polynomials (pressure p).
 //
 // The solvers being compared include:
-//    1. The divergence free solver (couple and decoupled modes)
-//    2. MINRES preconditioned by a block diagonal preconditioner
+//    1. MINRES preconditioned by a block diagonal preconditioner
+//    2. The divergence free solver (couple and decoupled modes)
+//    3. A block hybridization solver
+//    4. The Bramble-Pasciak solver (using BPCG or regular PCG)
 //
 // We recommend viewing example 5 before viewing this miniapp.
 //
 // Sample runs:
 //
-//    mpirun -np 8 block-solvers -r 2 -o 0
+//    mpirun -np 8 block-solvers -pr 2 -o 0
 //    mpirun -np 8 block-solvers -m anisotropic.mesh -c anisotropic.coeff -eb anisotropic.bdr
 //
 //
 // NOTE:  The coefficient file (provided through -c) defines a piecewise constant
-//        scalar coefficient k. The number of entries in this file should equal
-//        to the number of "element attributes" in the mesh file. The value of
-//        the coefficient in elements with the i-th attribute is given by the
-//        i-th entry of the coefficient file.
+//        scalar coefficient k. The number of entries in this file must equal the
+//        number of "element attributes" in the mesh file. The value of the
+//        coefficient in elements with the i-th attribute is given by the i-th
+//        entry of the coefficient file.
 //
 //
 // NOTE:  The essential boundary attribute file (provided through -eb) defines
-//        which attributes to impose essential boundary condition (on u). The
-//        number of entries in this file should equal to the number of "boundary
+//        which attributes to impose an essential boundary condition (on u).
+//        The number of entries in this file must equal the number of "boundary
 //        attributes" in the mesh file. If the i-th entry of the file is nonzero
-//        (respectively 0), essential (respectively natural) boundary condition
-//        will be imposed on boundary with the i-th attribute.
+//        (respectively zero), then the essential (respectively natural) boundary
+//        condition will be imposed on boundary elements with the i-th attribute.

 #include "mfem.hpp"
+#include "bramble_pasciak.hpp"
 #include "div_free_solver.hpp"
+#include "block_hybridization.hpp"
 #include <fstream>
 #include <iostream>
 #include <memory>
@@ -62,11 +66,12 @@ using namespace mfem;
 using namespace blocksolvers;

 // Exact solution, u and p, and r.h.s., f and g.
-void u_exact(const Vector & x, Vector & u);
-double p_exact(const Vector & x);
-void f_exact(const Vector & x, Vector & f);
-double g_exact(const Vector & x);
-double natural_bc(const Vector & x);
+void u_exact(const Vector &x, Vector &u);
+double p_exact(const Vector &x);
+void f_exact(const Vector &x, Vector &f);
+double g_exact(const Vector &x);
+double natural_bc(const Vector &x);
+double pi = M_PI;

 /** Wrapper for assembling the discrete Darcy problem (ex5p)
                     [ M  B^T ] [u] = [f]
@@ -81,35 +86,47 @@ double natural_bc(const Vector & x);
       D: subset of the boundary where natural boundary condition is imposed. */
 class DarcyProblem
 {
-   OperatorPtr M_;
-   OperatorPtr B_;
+   Array<int> offsets_;
   Vector rhs_;
   Vector ess_data_;
   ParGridFunction u_;
   ParGridFunction p_;
   ParMesh mesh_;
+   shared_ptr<ParBilinearForm> Mform_;
+   shared_ptr<ParMixedBilinearForm> Bform_;
+   unique_ptr<PWConstCoefficient> mass_coeff_;
   VectorFunctionCoefficient ucoeff_;
   FunctionCoefficient pcoeff_;
+   const Array<int> &ess_bdr_;
   DFSSpaces dfs_spaces_;
   const IntegrationRule *irs_[Geometry::NumGeom];
-public:
-   DarcyProblem(Mesh &mesh, int num_refines, int order, const char *coef_file,
-                Array<int> &ess_bdr, DFSParameters param);

-   HypreParMatrix& GetM() { return *M_.As<HypreParMatrix>(); }
-   HypreParMatrix& GetB() { return *B_.As<HypreParMatrix>(); }
-   const Vector& GetRHS() { return rhs_; }
-   const Vector& GetEssentialBC() { return ess_data_; }
-   const DFSData& GetDFSData() const { return dfs_spaces_.GetDFSData(); }
+   void Distribute(const Vector &x);
+public:
+   DarcyProblem(MPI_Comm comm, Mesh &mesh, int num_refines, int order,
+                const char *coef_file, Array<int> &ess_bdr, DFSParameters param);
+
+   void GetParallelSystems(shared_ptr<HypreParMatrix> &M,
+                           shared_ptr<HypreParMatrix> &B,
+                           shared_ptr<HypreParMatrix> &M_e,
+                           shared_ptr<HypreParMatrix> &B_e,
+                           mfem::Array<int> &ess_tdof_list) const;
+
   void ShowError(const Vector &sol, bool verbose);
   void VisualizeSolution(const Vector &sol, string tag);
+
+   const Vector &GetRHS() const { return rhs_; }
+   const Vector &GetEssentialBC() const { return ess_data_; }
+   shared_ptr<ParBilinearForm> GetMform() const { return Mform_; }
+   shared_ptr<ParMixedBilinearForm> GetBform() const { return Bform_; }
+   const DFSData &GetDFSData() const { return dfs_spaces_.GetDFSData(); }
 };

-DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
+DarcyProblem::DarcyProblem(MPI_Comm comm, Mesh &mesh, int num_refs, int order,
                           const char *coef_file, Array<int> &ess_bdr,
                           DFSParameters dfs_param)
-   : mesh_(MPI_COMM_WORLD, mesh), ucoeff_(mesh.Dimension(), u_exact),
-     pcoeff_(p_exact), dfs_spaces_(order, num_refs, &mesh_, ess_bdr, dfs_param)
+   : mesh_(comm, mesh), ucoeff_(mesh.Dimension(), u_exact), pcoeff_(p_exact),
+     ess_bdr_(ess_bdr), dfs_spaces_(order, num_refs, &mesh_, ess_bdr, dfs_param)
 {
   for (int l = 0; l < num_refs; l++)
   {
@@ -124,51 +141,54 @@ DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
      ifstream coef_str(coef_file);
      coef_vector.Load(coef_str, mesh.GetNE());
   }
-   PWConstCoefficient mass_coeff(coef_vector);
+   mass_coeff_.reset(new PWConstCoefficient(coef_vector));
+
   VectorFunctionCoefficient fcoeff(mesh_.Dimension(), f_exact);
   FunctionCoefficient natcoeff(natural_bc);
   FunctionCoefficient gcoeff(g_exact);

-   u_.SetSpace(dfs_spaces_.GetHdivFES());
-   p_.SetSpace(dfs_spaces_.GetL2FES());
+   ParFiniteElementSpace* u_fes = dfs_spaces_.GetHdivFES();
+   ParFiniteElementSpace* p_fes = dfs_spaces_.GetL2FES();
+   offsets_.SetSize(3, 0);
+   offsets_[1] = u_fes->GetTrueVSize();
+   offsets_[2] = offsets_[1] + p_fes->GetTrueVSize();
+
+   u_.SetSpace(u_fes);
+   p_.SetSpace(p_fes);
   p_ = 0.0;
   u_ = 0.0;
   u_.ProjectBdrCoefficientNormal(ucoeff_, ess_bdr);

-   ParLinearForm fform(dfs_spaces_.GetHdivFES());
-   fform.AddDomainIntegrator(new VectorFEDomainLFIntegrator(fcoeff));
+   ParLinearForm fform(u_fes);
+   // fform.AddDomainIntegrator(new VectorFEDomainLFIntegrator(fcoeff));
   fform.AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(natcoeff));
   fform.Assemble();

-   ParLinearForm gform(dfs_spaces_.GetL2FES());
+   ParLinearForm gform(p_fes);
   gform.AddDomainIntegrator(new DomainLFIntegrator(gcoeff));
   gform.Assemble();

-   ParBilinearForm mVarf(dfs_spaces_.GetHdivFES());
-   ParMixedBilinearForm bVarf(dfs_spaces_.GetHdivFES(), dfs_spaces_.GetL2FES());
+   rhs_.SetSize(offsets_[2]);
+   BlockVector blk_rhs(rhs_, offsets_);
+   fform.ParallelAssemble(blk_rhs.GetBlock(0));
+   gform.ParallelAssemble(blk_rhs.GetBlock(1));

-   mVarf.AddDomainIntegrator(new VectorFEMassIntegrator(mass_coeff));
-   mVarf.Assemble();
-   mVarf.EliminateEssentialBC(ess_bdr, u_, fform);
-   mVarf.Finalize();
-   M_.Reset(mVarf.ParallelAssemble());
+   Mform_ = make_shared<ParBilinearForm>(u_fes);
+   Mform_->AddDomainIntegrator(new VectorFEMassIntegrator(*mass_coeff_));
+   Mform_->ComputeElementMatrices();
+   Mform_->Assemble();
+   Mform_->Finalize();

-   bVarf.AddDomainIntegrator(new VectorFEDivergenceIntegrator);
-   bVarf.Assemble();
-   bVarf.SpMat() *= -1.0;
-   bVarf.EliminateTrialDofs(ess_bdr, u_, gform);
-   bVarf.Finalize();
-   B_.Reset(bVarf.ParallelAssemble());
+   Bform_ = make_shared<ParMixedBilinearForm>(u_fes, p_fes);
+   Bform_->AddDomainIntegrator(new VectorFEDivergenceIntegrator());
+   Bform_->ComputeElementMatrices();
+   Bform_->Assemble();
+   Bform_->SpMat() *= -1.0;
+   Bform_->Finalize();

-   rhs_.SetSize(M_->NumRows() + B_->NumRows());
-   Vector rhs_block0(rhs_.GetData(), M_->NumRows());
-   Vector rhs_block1(rhs_.GetData()+M_->NumRows(), B_->NumRows());
-   fform.ParallelAssemble(rhs_block0);
-   gform.ParallelAssemble(rhs_block1);
-
-   ess_data_.SetSize(M_->NumRows() + B_->NumRows());
+   ess_data_.SetSize(rhs_.Size());
   ess_data_ = 0.0;
-   Vector ess_data_block0(ess_data_.GetData(), M_->NumRows());
+   Vector ess_data_block0(ess_data_.GetData(), offsets_[1]);
   u_.ParallelProject(ess_data_block0);

   int order_quad = max(2, 2*order+1);
@@ -178,11 +198,29 @@ DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
   }
 }

+void DarcyProblem::Distribute(const Vector &x)
+{
+   BlockVector blk_x(x.GetData(), offsets_);
+   u_.Distribute(blk_x.GetBlock(0));
+   p_.Distribute(blk_x.GetBlock(1));
+}
+
+void DarcyProblem::GetParallelSystems(shared_ptr<HypreParMatrix> &M,
+                                      shared_ptr<HypreParMatrix> &B,
+                                      shared_ptr<HypreParMatrix> &M_e,
+                                      shared_ptr<HypreParMatrix> &B_e,
+                                      mfem::Array<int> &ess_tdof_list) const
+{
+   dfs_spaces_.GetHdivFES()->GetEssentialTrueDofs(ess_bdr_, ess_tdof_list);
+   M.reset(Mform_->ParallelAssemble());
+   M_e.reset(M->EliminateRowsCols(ess_tdof_list));
+   B.reset(Bform_->ParallelAssemble());
+   B_e.reset(B->EliminateCols(ess_tdof_list));
+}
+
 void DarcyProblem::ShowError(const Vector& sol, bool verbose)
 {
-   u_.Distribute(Vector(sol.GetData(), M_->NumRows()));
-   p_.Distribute(Vector(sol.GetData()+M_->NumRows(), B_->NumRows()));
-
+   Distribute(sol);
   double err_u  = u_.ComputeL2Error(ucoeff_, irs_);
   double norm_u = ComputeGlobalLpNorm(2, ucoeff_, mesh_, irs_);
   double err_p  = p_.ComputeL2Error(pcoeff_, irs_);
@@ -199,9 +237,7 @@ void DarcyProblem::VisualizeSolution(const Vector& sol, string tag)
   MPI_Comm_size(mesh_.GetComm(), &num_procs);
   MPI_Comm_rank(mesh_.GetComm(), &myid);

-   u_.Distribute(Vector(sol.GetData(), M_->NumRows()));
-   p_.Distribute(Vector(sol.GetData()+M_->NumRows(), B_->NumRows()));
-
+   Distribute(sol);
   const char vishost[] = "localhost";
   const int  visport   = 19916;
   socketstream u_sock(vishost, visport);
@@ -243,16 +279,24 @@ int main(int argc, char *argv[])
   const char *coef_file = "";
   const char *ess_bdr_attr_file = "";
   int order = 0;
-   int par_ref_levels = 2;
-   bool show_error = false;
+   int ser_ref_levels = 0;
+   int par_ref_levels = 0;
+   bool show_error = true;
   bool visualization = false;
+
   DFSParameters param;
+#ifdef MFEM_USE_LAPACK
+   BPSParameters bps_param;
+#endif
+
   OptionsParser args(argc, argv);
   args.AddOption(&mesh_file, "-m", "--mesh",
                  "Mesh file to use.");
   args.AddOption(&order, "-o", "--order",
                  "Finite element order (polynomial degree).");
-   args.AddOption(&par_ref_levels, "-r", "--ref",
+   args.AddOption(&ser_ref_levels, "-sr", "--serial-ref",
+                  "Number of serial refinement steps.");
+   args.AddOption(&par_ref_levels, "-pr", "--parallel-ref",
                  "Number of parallel refinement steps.");
   args.AddOption(&coef_file, "-c", "--coef",
                  "Coefficient file to use.");
@@ -280,16 +324,25 @@ int main(int argc, char *argv[])

   // Initialize the mesh, boundary attributes, and solver parameters
   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   // Mesh *mesh = new Mesh(Mesh::MakeCartesian2D(1, 1, Element::TRIANGLE, 1));
   int dim = mesh->Dimension();
-   int ser_ref_lvls =
-      (int)ceil(log(Mpi::WorldSize()/mesh->GetNE())/log(2.)/dim);
-   for (int i = 0; i < ser_ref_lvls; ++i)
+
+   for (int i = 0; i < ser_ref_levels; ++i)
   {
      mesh->UniformRefinement();
   }

+   if (Mpi::Root())
+   {
+      MFEM_ASSERT(Mpi::WorldSize() < mesh->GetNE(),
+                  "Not enough elements in the mesh to be distributed:\n"
+                  << "Number of processors: " << Mpi::WorldSize() << "\n"
+                  << "Number of elements:   " << mesh->GetNE());
+   }
+
   Array<int> ess_bdr(mesh->bdr_attributes.Max());
   ess_bdr = 0;
+   // ess_bdr[0] = 0;
   if (std::strcmp(ess_bdr_attr_file, ""))
   {
      ifstream ess_bdr_attr_str(ess_bdr_attr_file);
@@ -312,9 +365,11 @@ int main(int argc, char *argv[])
   ResetTimer();

   // Generate components of the saddle point problem
-   DarcyProblem darcy(*mesh, par_ref_levels, order, coef_file, ess_bdr, param);
-   HypreParMatrix& M = darcy.GetM();
-   HypreParMatrix& B = darcy.GetB();
+   mfem::Array<int> ess_tdof_list;
+   shared_ptr<HypreParMatrix> M, B, M_e, B_e;
+   DarcyProblem darcy(MPI_COMM_WORLD, *mesh, par_ref_levels, order,
+                      coef_file, ess_bdr, param);
+   darcy.GetParallelSystems(M, B, M_e, B_e, ess_tdof_list);
   const DFSData& DFS_data = darcy.GetDFSData();
   delete mesh;

@@ -322,7 +377,7 @@ int main(int argc, char *argv[])
   {
      cout << line << "System assembled in " << chrono.RealTime() << "s.\n";
      cout << "Dimension of the physical space: " << dim << "\n";
-      cout << "Size of the discrete Darcy system: " << M.M() + B.M() << "\n";
+      cout << "Size of the discrete Darcy system: " << M->M() + B->M() << "\n";
      if (par_ref_levels > 0)
      {
         cout << "Dimension of the divergence free subspace: "
@@ -332,23 +387,57 @@ int main(int argc, char *argv[])

   // Setup various solvers for the discrete problem
   std::map<const DarcySolver*, double> setup_time;
+
   ResetTimer();
-   BDPMinresSolver bdp(M, B, param);
+   BlockHybridizationSolver bh(darcy.GetMform(), darcy.GetBform(), param, ess_bdr);
+   bh.SetEliminatedSystems(M_e, B_e, ess_tdof_list);
+   setup_time[&bh] = chrono.RealTime();
+   /*
+   ResetTimer();
+   BDPMinresSolver bdp(*M, *B, param);
+   bdp.SetEliminatedSystems(M_e, B_e, ess_tdof_list);
   setup_time[&bdp] = chrono.RealTime();

   ResetTimer();
-   DivFreeSolver dfs_dm(M, B, DFS_data);
+   DivFreeSolver dfs_dm(*M, *B, DFS_data);
+   dfs_dm.SetEliminatedSystems(M_e, B_e, ess_tdof_list);
   setup_time[&dfs_dm] = chrono.RealTime();

   ResetTimer();
   const_cast<bool&>(DFS_data.param.coupled_solve) = true;
-   DivFreeSolver dfs_cm(M, B, DFS_data);
+   DivFreeSolver dfs_cm(*M, *B, DFS_data);
+   dfs_cm.SetEliminatedSystems(M_e, B_e, ess_tdof_list);
   setup_time[&dfs_cm] = chrono.RealTime();

+#ifdef MFEM_USE_LAPACK
+   ResetTimer();
+   BramblePasciakSolver bp_bpcg(darcy.GetMform(), darcy.GetBform(), ess_tdof_list,
+                                bps_param);
+   bp_bpcg.SetEliminatedSystems(M_e, B_e, ess_tdof_list);
+   setup_time[&bp_bpcg] = chrono.RealTime();
+
+   ResetTimer();
+   bps_param.use_bpcg = false;
+   BramblePasciakSolver bp_pcg(darcy.GetMform(), darcy.GetBform(), ess_tdof_list,
+                               bps_param);
+   bp_pcg.SetEliminatedSystems(M_e, B_e, ess_tdof_list);
+   setup_time[&bp_pcg] = chrono.RealTime();
+#else
+   MFEM_WARNING("BramblePasciakSolver class unavailable: Compiled without LAPACK");
+#endif
+   */
+
   std::map<const DarcySolver*, std::string> solver_to_name;
+   solver_to_name[&bh] = "Block hybridization";
+   /*
   solver_to_name[&bdp] = "Block-diagonal-preconditioned MINRES";
   solver_to_name[&dfs_dm] = "Divergence free (decoupled mode)";
   solver_to_name[&dfs_cm] = "Divergence free (coupled mode)";
+#ifdef MFEM_USE_LAPACK
+   solver_to_name[&bp_bpcg] = "Bramble Pasciak CG (using BPCG)";
+   solver_to_name[&bp_pcg] = "Bramble Pasciak CG (using regular PCG)";
+#endif
+   */

   // Solve the problem using all solvers
   for (const auto& solver_pair : solver_to_name)
@@ -377,7 +466,7 @@ int main(int argc, char *argv[])
      else if (show_error && Mpi::Root())
      {
         cout << "Exact solution is unknown for coefficient '" << coef_file
-              << "'.\nApproximation error is computed in this case!\n\n";
+              << "'.\nApproximation error is not computed in this case!\n\n";
      }

      if (visualization) { darcy.VisualizeSolution(sol, name); }
@@ -387,40 +476,33 @@ int main(int argc, char *argv[])
   return 0;
 }

-void u_exact(const Vector & x, Vector & u)
+void u_exact(const Vector &x, Vector &u)
 {
   double xi(x(0));
-   double yi(x(1));
-   double zi(x.Size() == 3 ? x(2) : 0.0);
-
-   u(0) = - exp(xi)*sin(yi)*cos(zi);
-   u(1) = - exp(xi)*cos(yi)*cos(zi);
-   if (x.Size() == 3)
-   {
-      u(2) = exp(xi)*sin(yi)*sin(zi);
-   }
+   double xj(x(1));
+   u(0) = -exp(xi) * sin(xj);
+   u(1) = -exp(xi) * cos(xj);
 }

-double p_exact(const Vector & x)
+double p_exact(const Vector &x)
 {
   double xi(x(0));
-   double yi(x(1));
-   double zi(x.Size() == 3 ? x(2) : 0.0);
-   return exp(xi)*sin(yi)*cos(zi);
+   double xj(x(1));
+   return exp(xi) * sin(xj);
 }

-void f_exact(const Vector & x, Vector & f)
+void f_exact(const Vector &x, Vector &f)
 {
-   f = 0.0;
+   u_exact(x, f);
+   f(0) -= pi * cos(pi*x(0));
 }

-double g_exact(const Vector & x)
+double g_exact(const Vector &x)
 {
-   if (x.Size() == 3) { return -p_exact(x); }
-   return 0;
+   return 0.0;
 }

-double natural_bc(const Vector & x)
+double natural_bc(const Vector &x)
 {
-   return (-p_exact(x));
+   return -p_exact(x);
 }
@@ -0,0 +1,443 @@
+// 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.
+
+#include "block_hybridization.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace blocksolvers;
+
+namespace mfem
+{
+namespace blocksolvers
+{
+void BlockHybridizationSolver::Init(const int ne)
+{
+   hat_offsets.SetSize(ne + 1);
+   hat_offsets[0] = 0;
+   for (int i = 0; i < ne; ++i)
+   {
+      hat_offsets[i + 1] = trial_space.GetFE(i)->GetDof();
+   }
+   hat_offsets.PartialSum();
+
+   data_offsets.SetSize(ne + 1);
+   data_offsets[0] = 0;
+
+   ipiv_offsets.SetSize(ne + 1);
+   ipiv_offsets[0] = 0;
+
+   test_offsets.SetSize(ne + 1);
+   test_offsets[0] = 0;
+
+   for (int i = 0; i < ne; ++i)
+   {
+      const int trial_size = trial_space.GetFE(i)->GetDof();
+      test_offsets[i + 1] = test_space.GetFE(i)->GetDof();
+      const int matrix_size = trial_size + test_offsets[i + 1];
+
+      data_offsets[i + 1] = data_offsets[i] + matrix_size*matrix_size;
+      ipiv_offsets[i + 1] = ipiv_offsets[i] + matrix_size;
+   }
+   test_offsets.PartialSum();
+
+   data = new double[data_offsets.Last()]();
+   ipiv = new int[ipiv_offsets.Last()];
+
+   mixed_dofs.Reserve(ipiv_offsets.Last());
+}
+
+void BlockHybridizationSolver::ConstructCt(const ParFiniteElementSpace &c_space)
+{
+   ParMesh &pmesh(*trial_space.GetParMesh());
+   Ct = new SparseMatrix(hat_offsets.Last(), c_space.GetNDofs());
+   Array<int> dofs, c_dofs;
+   const double eps = 1e-12;
+   DenseMatrix elmat;
+   FaceElementTransformations *FTr;
+   NormalTraceJumpIntegrator c_int;
+   const int num_faces = pmesh.GetNumFaces();
+
+   for (int i = 0; i < num_faces; ++i)
+   {
+      FTr = pmesh.GetInteriorFaceTransformations(i);
+      if (!FTr)
+      {
+         continue;
+      }
+
+      int o1 = hat_offsets[FTr->Elem1No];
+      int s1 = hat_offsets[FTr->Elem1No + 1] - o1;
+      int o2 = hat_offsets[FTr->Elem2No];
+      int s2 = hat_offsets[FTr->Elem2No + 1] - o2;
+
+      dofs.SetSize(s1 + s2);
+      for (int j = 0; j < s1; ++j)
+      {
+         dofs[j] = o1 + j;
+      }
+      for (int j = 0; j < s2; ++j)
+      {
+         dofs[s1 + j] = o2 + j;
+      }
+      c_space.GetFaceDofs(i, c_dofs);
+      c_int.AssembleFaceMatrix(*c_space.GetFaceElement(i),
+                               *trial_space.GetFE(FTr->Elem1No),
+                               *trial_space.GetFE(FTr->Elem2No),
+                               *FTr,
+                               elmat);
+      elmat.Threshold(eps * elmat.MaxMaxNorm());
+      Ct->AddSubMatrix(dofs, c_dofs, elmat);
+   }
+
+   const int num_shared_faces = pmesh.GetNSharedFaces();
+   for (int i = 0; i < num_shared_faces; ++i)
+   {
+      const int face_no = pmesh.GetSharedFace(i);
+      FTr = pmesh.GetFaceElementTransformations(face_no);
+      c_space.GetFaceDofs(face_no, c_dofs);
+      const FiniteElement *face_fe(c_space.GetFaceElement(face_no));
+      const FiniteElement *fe(trial_space.GetFE(FTr->Elem1No));
+
+      int o1 = hat_offsets[FTr->Elem1No];
+      int s1 = hat_offsets[FTr->Elem1No + 1] - o1;
+
+      dofs.SetSize(s1);
+      for (int j = 0; j < s1; ++j)
+      {
+         dofs[j] = o1 + j;
+      }
+      c_int.AssembleFaceMatrix(*face_fe, *fe, *fe, *FTr, elmat);
+      elmat.Threshold(eps * elmat.MaxMaxNorm());
+      Ct->AddSubMatrix(dofs, c_dofs, elmat);
+   }
+   Ct->Finalize();
+}
+
+void BlockHybridizationSolver::ConstructH(const shared_ptr<ParBilinearForm> &a,
+                                          const shared_ptr<ParMixedBilinearForm> &b,
+                                          const Array<int> &marker,
+                                          const ParFiniteElementSpace &c_space)
+{
+   ParMesh &pmesh(*trial_space.GetParMesh());
+   const int ne = pmesh.GetNE();
+   const double eps = 1e-12;
+   SparseMatrix H(Ct->Width());
+   DenseMatrix Ct_local, Minv_Ct_local, H_local;
+
+   Array<int> dofs, c_dofs;
+   Array<int> c_dof_marker(Ct->Width());
+   c_dof_marker = -1;
+   int c_mark_start = 0;
+
+   for (int i = 0; i < ne; ++i)
+   {
+      trial_space.GetElementDofs(i, dofs);
+      const int trial_size = dofs.Size();
+
+      DenseMatrix A(trial_size);
+      a->ComputeElementMatrix(i, A);
+      A.Threshold(eps * A.MaxMaxNorm());
+
+      const int test_offset = test_offsets[i];
+      const int test_size = test_offsets[i + 1] - test_offset;
+
+      DenseMatrix B(test_size, trial_size);
+      b->ComputeElementMatrix(i, B);
+      B.Neg();
+      B.Threshold(eps * B.MaxMaxNorm());
+
+      const int matrix_size = trial_size + test_size;
+      DenseMatrix M(data + data_offsets[i], matrix_size, matrix_size);
+
+      M.CopyMN(A, 0, 0);
+      M.CopyMN(B, trial_size, 0);
+      M.CopyMNt(B, 0, trial_size);
+
+      if (elimination_)
+      {
+         FiniteElementSpace::AdjustVDofs(dofs);
+         for (int j = 0; j < trial_size; ++j)
+         {
+            if (marker[dofs[j]])
+            {
+               for (int k = 0; k < matrix_size; ++k)
+               {
+                  if (k == j)
+                  {
+                     M(k, k) = 1.0;
+                  }
+                  else
+                  {
+                     M(k, j) = 0.0;
+                     M(j, k) = 0.0;
+                  }
+               }
+            }
+         }
+      }
+
+      c_dofs.SetSize(0);
+      dofs.SetSize(trial_size);
+      const int hat_offset = hat_offsets[i];
+
+      for (int j = 0; j < trial_size; ++j)
+      {
+         const int row = hat_offset + j;
+         const int ncols = Ct->RowSize(row);
+         const int *cols = Ct->GetRowColumns(row);
+         for (int l = 0; l < ncols; ++l)
+         {
+            const int c_dof = cols[l];
+            if (c_dof_marker[c_dof] < c_mark_start)
+            {
+               c_dof_marker[c_dof] = c_mark_start + c_dofs.Size();
+               c_dofs.Append(c_dof);
+            }
+         }
+         dofs[j] = row;
+      }
+
+      mixed_dofs.Append(dofs);
+      for (int j = 0; j < test_size; ++j)
+      {
+         mixed_dofs.Append(hat_offsets.Last() + test_offset + j);
+      }
+
+      Ct_local.SetSize(M.Height(), c_dofs.Size()); // Ct_local = [C 0]^T
+      Ct_local = 0.0;
+      for (int j = 0; j < trial_size; ++j)
+      {
+         const int row = dofs[j];
+         const int ncols = Ct->RowSize(row);
+         const int *cols = Ct->GetRowColumns(row);
+         const double *vals = Ct->GetRowEntries(row);
+         for (int l = 0; l < ncols; ++l)
+         {
+            const int loc = c_dof_marker[cols[l]] - c_mark_start;
+            Ct_local(j, loc) = vals[l];
+         }
+      }
+
+      LUFactors Minv(data + data_offsets[i], ipiv + ipiv_offsets[i]);
+      Minv.Factor(matrix_size);
+      Minv_Ct_local = Ct_local;
+      Minv.Solve(Ct_local.Height(), Ct_local.Width(), Minv_Ct_local.Data());
+
+      H_local.SetSize(Ct_local.Width());
+
+      MultAtB(Ct_local, Minv_Ct_local, H_local);
+      H.AddSubMatrix(c_dofs, c_dofs, H_local);
+
+      c_mark_start += c_dofs.Size();
+      MFEM_VERIFY(c_mark_start >= 0, "overflow");
+   }
+   H.Finalize(1, true);  // skip_zeros = 1 (default), fix_empty_rows = true
+
+   OperatorPtr pP(Operator::Hypre_ParCSR);
+   pP.ConvertFrom(c_space.Dof_TrueDof_Matrix());
+   OperatorPtr dH(pP.Type());
+   dH.MakeSquareBlockDiag(c_space.GetComm(), c_space.GlobalVSize(),
+                          c_space.GetDofOffsets(), &H);
+   OperatorPtr AP(ParMult(dH.As<HypreParMatrix>(), pP.As<HypreParMatrix>()));
+   OperatorPtr R(pP.As<HypreParMatrix>()->Transpose());
+   pH = ParMult(R.As<HypreParMatrix>(), AP.As<HypreParMatrix>(), true);
+}
+
+BlockHybridizationSolver::BlockHybridizationSolver(const
+                                                   shared_ptr<ParBilinearForm> &a,
+                                                   const shared_ptr<ParMixedBilinearForm> &b,
+                                                   const IterSolveParameters &param,
+                                                   const Array<int> &ess_bdr_attr)
+   : DarcySolver(a->ParFESpace()->GetTrueVSize(),
+                 b->TestParFESpace()->GetTrueVSize()),
+     trial_space(*a->ParFESpace()), test_space(*b->TestParFESpace()),
+     elimination_(false),
+     solver_(a->ParFESpace()->GetComm())
+{
+   ParMesh &pmesh(*trial_space.GetParMesh());
+   const int ne = pmesh.GetNE();
+
+   StopWatch chrono;
+
+   chrono.Start();
+   Init(ne);
+   chrono.Stop();
+   cout << "init time: " << chrono.RealTime() << endl;
+
+   Array<int> ess_dof_marker;
+   for (int attr : ess_bdr_attr)
+   {
+      if (attr)
+      {
+         elimination_ = true;
+         trial_space.GetEssentialVDofs(ess_bdr_attr, ess_dof_marker);
+         break;
+      }
+   }
+
+   chrono.Restart();
+   const int order = trial_space.FEColl()->GetOrder()-1;
+   DG_Interface_FECollection fec(order, pmesh.Dimension());
+   c_fes = new ParFiniteElementSpace(&pmesh, &fec);
+   ParFiniteElementSpace &c_space(*c_fes);
+
+   ConstructCt(c_space);
+   chrono.Stop();
+   cout << "ct time: " << chrono.RealTime() << endl;
+
+   chrono.Restart();
+   ConstructH(a, b, ess_dof_marker, c_space);
+   chrono.Stop();
+   cout << "h time: " << chrono.RealTime() << endl;
+
+   M = new HypreBoomerAMG(*pH);
+   M->SetPrintLevel(0);
+
+   SetOptions(solver_, param);
+   solver_.SetPreconditioner(*M);
+   solver_.SetOperator(*pH);
+}
+
+BlockHybridizationSolver::~BlockHybridizationSolver()
+{
+   delete M;
+   delete [] ipiv;
+   delete [] data;
+   delete Ct;
+   delete c_fes;
+}
+
+void BlockHybridizationSolver::ReduceRHS(const Vector &b,
+                                         const Vector &sol,
+                                         BlockVector &rhs,
+                                         Vector &b_r) const
+{
+   const SparseMatrix &R = *trial_space.GetRestrictionMatrix();
+   Vector x_e(b);
+   if (elimination_) { EliminateEssentialBC(sol, x_e);}
+   Vector x0(R.Width());
+   BlockVector block_x(x_e.GetData(), offsets_);
+   R.MultTranspose(block_x.GetBlock(0), x0);
+
+   ParMesh &pmesh(*trial_space.GetParMesh());
+   const int ne = pmesh.GetNE();
+
+   rhs.SetVector(block_x.GetBlock(1), hat_offsets.Last());
+   Array<bool> dof_marker(x0.Size());
+   dof_marker = false;
+   Array<int> dofs;
+   Vector Minv_sub_vec, g_i;
+
+   for (int i = 0; i < ne; ++i)
+   {
+      trial_space.GetElementDofs(i, dofs);
+      const int trial_size = dofs.Size();
+      g_i.MakeRef(rhs, hat_offsets[i], trial_size);
+      x0.GetSubVector(dofs, g_i);  // reverses the sign if dof < 0
+
+      trial_space.AdjustVDofs(dofs);
+      for (int j = 0; j < trial_size; ++j)
+      {
+         int dof = dofs[j];
+         if (dof_marker[dof])
+         {
+            g_i(j) = 0.0;
+         }
+         else
+         {
+            dof_marker[dof] = true;
+         }
+      }
+
+      const int ipiv_offset = ipiv_offsets[i];
+      const int matrix_size = ipiv_offsets[i + 1] - ipiv_offset;
+
+      mixed_dofs.GetSubArray(ipiv_offset, matrix_size, dofs);
+
+      LUFactors Minv(data + data_offsets[i], ipiv + ipiv_offset);
+      rhs.GetSubVector(dofs, Minv_sub_vec);
+      Minv.Solve(Minv_sub_vec.Size(), 1, Minv_sub_vec.GetData());
+      rhs.SetSubVector(dofs, Minv_sub_vec);
+   }
+   Ct->MultTranspose(rhs.GetBlock(0), b_r);
+}
+
+void BlockHybridizationSolver::ComputeSolution(Vector &y,
+                                               BlockVector &rhs,
+                                               const Vector &rhs_r,
+                                               Array<int> &block_offsets) const
+{
+   BlockVector Ct_lambda(block_offsets);
+   Ct->Mult(rhs_r, Ct_lambda.GetBlock(0));
+   Ct_lambda.GetBlock(1) = 0.0;
+
+   ParMesh &pmesh(*trial_space.GetParMesh());
+   const int ne = pmesh.GetNE();
+   Array<int> dofs;
+   Vector Minv_sub_vec;
+
+   for (int i = 0; i < ne; ++i)
+   {
+      const int ipiv_offset = ipiv_offsets[i];
+      const int matrix_size = ipiv_offsets[i + 1] - ipiv_offset;
+
+      mixed_dofs.GetSubArray(ipiv_offset, matrix_size, dofs);
+
+      LUFactors Minv(data + data_offsets[i], ipiv + ipiv_offset);
+      Ct_lambda.GetSubVector(dofs, Minv_sub_vec);
+      Minv.Solve(Minv_sub_vec.Size(), 1, Minv_sub_vec.GetData());
+      Minv_sub_vec.Neg();
+      rhs.AddElementVector(dofs, Minv_sub_vec);
+   }
+
+   const SparseMatrix &R = *trial_space.GetRestrictionMatrix();
+   Vector x0(R.Width());
+   Vector sub_vec;
+   for (int i = 0; i < ne; ++i)
+   {
+      trial_space.GetElementDofs(i, dofs);
+      sub_vec.MakeRef(rhs.GetBlock(0), hat_offsets[i], dofs.Size());
+      x0.SetSubVector(dofs, sub_vec);
+   }
+   BlockVector block_y(y, offsets_);
+   R.Mult(x0, block_y.GetBlock(0));
+   y.SetVector(rhs.GetBlock(1), offsets_[1]);
+}
+
+void BlockHybridizationSolver::Mult(const Vector &x, Vector &y) const
+{
+   Array<int> block_offsets(3);
+   block_offsets[0] = 0;
+   block_offsets[1] = hat_offsets.Last();
+   block_offsets[2] = offsets_[2] - offsets_[1];
+   block_offsets.PartialSum();
+
+   BlockVector rhs(block_offsets);
+
+   Vector rhs_r(Ct->Width());
+   rhs_r = 0.0;
+   ReduceRHS(x, y, rhs, rhs_r);
+
+   Vector rhs_true(pH->Height());
+   const Operator &P(*c_fes->GetProlongationMatrix());
+   P.MultTranspose(rhs_r, rhs_true);
+
+   Vector lambda_true(rhs_true.Size());
+   lambda_true = 0.0;
+
+   solver_.Mult(rhs_true, lambda_true);
+
+   P.Mult(lambda_true, rhs_r);
+   ComputeSolution(y, rhs, rhs_r, block_offsets);
+}
+} // namespace blocksolvers
+} // namespace mfem
@@ -0,0 +1,62 @@
+// 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.
+
+#ifndef MFEM_BLOCK_HYBRIDIZATION_SOLVER_HPP
+#define MFEM_BLOCK_HYBRIDIZATION_SOLVER_HPP
+
+#include "darcy_solver.hpp"
+
+namespace mfem
+{
+namespace blocksolvers
+{
+
+
+class BlockHybridizationSolver : public DarcySolver
+{
+   ParFiniteElementSpace trial_space, test_space;
+   ParFiniteElementSpace *c_fes;
+   Array<int> hat_offsets, test_offsets, data_offsets, ipiv_offsets, mixed_dofs;
+   double *data;
+   int *ipiv;
+   bool elimination_;
+   SparseMatrix *Ct;
+   HypreBoomerAMG *M;
+   HypreParMatrix *pH;
+   CGSolver solver_;
+
+   void Init(const int ne);
+   void ConstructCt(const ParFiniteElementSpace &c_space);
+   void ConstructH(const std::shared_ptr<ParBilinearForm> &a,
+                   const std::shared_ptr<ParMixedBilinearForm> &b,
+                   const Array<int> &marker,
+                   const ParFiniteElementSpace &c_space);
+   void ReduceRHS(const Vector &b, const Vector &sol, BlockVector &rhs,
+                  Vector &b_r) const;
+   void ComputeSolution(Vector &y,
+                        BlockVector &rhs,
+                        const Vector &rhs_r,
+                        Array<int> &block_offsets) const;
+public:
+   BlockHybridizationSolver(const std::shared_ptr<ParBilinearForm> &a,
+                            const std::shared_ptr<ParMixedBilinearForm> &b,
+                            const IterSolveParameters &param,
+                            const Array<int> &ess_bdr_attr);
+   ~BlockHybridizationSolver();
+   virtual void Mult(const Vector &x, Vector &y) const;
+   virtual void SetOperator(const Operator &op) { }
+   virtual int GetNumIterations() const { return solver_.GetNumIterations(); }
+};
+
+} // namespace blocksolvers
+} // namespace mfem
+
+#endif // MFEM_BLOCK_HYBRIDIZATION_SOLVER_HPP
@@ -0,0 +1,400 @@
+// Copyright (c) 2023-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.
+
+
+#include "bramble_pasciak.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace blocksolvers;
+
+namespace mfem
+{
+namespace blocksolvers
+{
+/// Bramble-Pasciak Solver
+BramblePasciakSolver::BramblePasciakSolver(
+   const std::shared_ptr<ParBilinearForm> &mVarf,
+   const std::shared_ptr<ParMixedBilinearForm> &bVarf,
+   const Array<int>& ess_tdof_list,
+   const BPSParameters &param)
+   : DarcySolver(mVarf->ParFESpace()->GetTrueVSize(),
+                 bVarf->TestFESpace()->GetTrueVSize())
+{
+   // Cf. DarcyProblem::GetParallelSystems
+   M_.reset(mVarf->ParallelAssemble());
+   M_e_.reset(M_->EliminateRowsCols(ess_tdof_list));
+   B_.reset(bVarf->ParallelAssemble());
+   B_e_.reset(B_->EliminateCols(ess_tdof_list));
+   Q_.reset(ConstructMassPreconditioner(*mVarf, param.q_scaling));
+   Q_e_.reset(Q_->EliminateRowsCols(ess_tdof_list));
+
+   Vector diagM;
+   M_->GetDiag(diagM);
+   auto BT = B_->Transpose();
+   auto invDBt = new HypreParMatrix(*BT);
+   invDBt->InvScaleRows(diagM);
+   auto S = ParMult(B_.get(), invDBt);
+   M0_.Reset(new HypreDiagScale(*M_));
+   M1_.Reset(new HypreBoomerAMG(*S));
+   M1_.As<HypreBoomerAMG>()->SetPrintLevel(0);
+
+   Init(*M_, *B_, *Q_, *M0_.As<Solver>(), *M1_.As<Solver>(), param);
+}
+
+BramblePasciakSolver::BramblePasciakSolver(
+   const std::shared_ptr<ParBilinearForm> &mVarf,
+   const std::shared_ptr<ParMixedBilinearForm> &bVarf,
+   const BPSParameters &param)
+   : DarcySolver(mVarf->ParFESpace()->GetTrueVSize(),
+                 bVarf->TestFESpace()->GetTrueVSize())
+{
+   M_.reset(mVarf->ParallelAssemble());
+   B_.reset(bVarf->ParallelAssemble());
+   Q_.reset(ConstructMassPreconditioner(*mVarf, param.q_scaling));
+
+   Vector diagM;
+   M_->GetDiag(diagM);
+   auto BT = B_->Transpose();
+   auto invDBt = new HypreParMatrix(*BT);
+   invDBt->InvScaleRows(diagM);
+   auto S = ParMult(B_.get(), invDBt);
+   M0_.Reset(new HypreDiagScale(*M_));
+   M1_.Reset(new HypreBoomerAMG(*S));
+   M1_.As<HypreBoomerAMG>()->SetPrintLevel(0);
+
+   Init(*M_, *B_, *Q_, *M0_.As<Solver>(), *M1_.As<Solver>(), param);
+}
+
+BramblePasciakSolver::BramblePasciakSolver(
+   HypreParMatrix &M, HypreParMatrix &B, HypreParMatrix &Q,
+   Solver &M0, Solver &M1,
+   const BPSParameters &param)
+   : DarcySolver(M.NumRows(), B.NumRows())
+{
+   Init(M, B, Q, M0, M1, param);
+}
+
+void BramblePasciakSolver::Init(
+   HypreParMatrix &M, HypreParMatrix &B, HypreParMatrix &Q,
+   Solver &M0, Solver &M1,
+   const BPSParameters &param)
+{
+   auto Bt = new TransposeOperator(&B);
+   auto invQ = new HypreDiagScale(Q);
+
+   use_bpcg = param.use_bpcg;
+
+   if (use_bpcg)
+   {
+      oop_ = new BlockOperator(offsets_);
+      oop_->owns_blocks = false;
+      oop_->SetBlock(0, 0, &M);
+      oop_->SetBlock(0, 1, Bt);
+      oop_->SetBlock(1, 0, &B);
+
+      // cpc_ unused in bpcg
+      auto temp_cpc = new BlockDiagonalPreconditioner(offsets_);
+      temp_cpc->owns_blocks = true;
+      temp_cpc->SetDiagonalBlock(0, invQ);
+      temp_cpc->SetDiagonalBlock(1, &M1);
+      // tri(1,0) = B M0 = B invQ
+      auto id_m = new IdentityOperator(M.NumRows());
+      auto id_b = new IdentityOperator(B.NumRows());
+      auto BinvQ = new ProductOperator(&B, invQ, false, false);
+      // tri
+      auto temp_tri = new BlockOperator(offsets_);
+      temp_tri->owns_blocks = true;
+      temp_tri->SetBlock(0, 0, id_m);
+      temp_tri->SetBlock(1, 1, id_b, -1.0);
+      temp_tri->SetBlock(1, 0, BinvQ);
+
+      ppc_ = new ProductOperator(temp_cpc, temp_tri, true, true);
+
+      ipc_ = new BlockOperator(offsets_);
+      ipc_->owns_blocks = false;
+      ipc_->SetDiagonalBlock(0, invQ);
+
+      // bpcg
+      solver_.reset(new BPCGSolver(M.GetComm(), *ipc_, *ppc_));
+      solver_->SetOperator(*oop_);
+   }
+   else
+   {
+      // oop_ unused in cg
+      auto temp_oop = new BlockOperator(offsets_);
+      temp_oop->owns_blocks = false;
+      temp_oop->SetBlock(0, 0, &M);
+      temp_oop->SetBlock(0, 1, Bt);
+      temp_oop->SetBlock(1, 0, &B);
+
+      // ipc_ unused in cg
+      auto temp_ipc = new BlockOperator(offsets_);
+      temp_ipc->owns_blocks = false;
+      temp_ipc->SetDiagonalBlock(0, invQ);
+
+      // temp_AN = temp_oop * temp_ipc
+      auto temp_AN = new ProductOperator(temp_oop, temp_ipc, true, true);
+
+      // Required for updating the RHS
+      auto id = new IdentityOperator(M.NumRows()+B.NumRows());
+      map_ = new SumOperator(temp_AN, 1.0, id, -1.0, true, true);
+
+      mop_ = new ProductOperator(map_, temp_oop, false, true);
+
+      cpc_ = new BlockDiagonalPreconditioner(offsets_);
+      cpc_->owns_blocks = true;
+      cpc_->SetDiagonalBlock(0, &M0);
+      cpc_->SetDiagonalBlock(1, &M1);
+
+      // (P)CG
+      solver_.reset(new CGSolver(M.GetComm()));
+      solver_->SetOperator(*mop_);
+      solver_->SetPreconditioner(*cpc_);
+   }
+   SetOptions(*solver_, param);
+}
+
+HypreParMatrix *BramblePasciakSolver::ConstructMassPreconditioner(
+   ParBilinearForm &mVarf, double q_scaling)
+{
+#ifdef MFEM_USE_LAPACK
+   MFEM_ASSERT((q_scaling > 0.0) && (q_scaling < 1.0),
+               "Invalid Q-scaling factor: q_scaling = " << q_scaling );
+   ParBilinearForm qVarf(mVarf.ParFESpace());
+   for (int i = 0; i < mVarf.ParFESpace()->GetNE(); ++i)
+   {
+      DenseMatrix M_i, Q_i, evec;
+      Vector eval, diag_i;
+      double scaling = 0.0;
+
+      mVarf.ComputeElementMatrix(i, M_i);
+      M_i.GetDiag(diag_i);
+      // M_i <- D^{-1/2} M_i D^{-1/2}, where D = diag(M_i)
+      M_i.InvSymmetricScaling(diag_i);
+      // M_i x = ev diag(M_i) x
+      M_i.Eigenvalues(eval, evec);
+
+      scaling = q_scaling*eval.Min();
+      diag_i.Set(scaling, diag_i);
+      Q_i.Diag(diag_i.GetData(), diag_i.Size());
+      qVarf.AssembleElementMatrix(i, Q_i, 1);
+   }
+   qVarf.Finalize();
+   return qVarf.ParallelAssemble();
+#else
+   MFEM_CONTRACT_VAR(mVarf);
+   MFEM_CONTRACT_VAR(q_scaling);
+   mfem_error("BramblePasciakSolver::ConstructMassPreconditioner: Compiled without LAPACK");
+   return nullptr;
+#endif
+}
+
+void BramblePasciakSolver::Mult(const Vector & x, Vector & y) const
+{
+   Vector x_e(x);
+   if (rhs_needs_elimination_) { EliminateEssentialBC(y, x_e);}
+   if (!use_bpcg)
+   {
+      Vector transformed_rhs(x_e.Size());
+      map_->Mult(x_e, transformed_rhs);
+      solver_->Mult(transformed_rhs, y);
+   }
+   else
+   {
+      solver_->Mult(x_e, y);
+   }
+   for (int dof : ess_zero_dofs_) { y[dof] = 0.0; }
+}
+
+/// Bramble-Pasciak CG
+void BPCGSolver::UpdateVectors()
+{
+   MemoryType mt = GetMemoryType(oper->GetMemoryClass());
+
+   r.SetSize(width, mt); r.UseDevice(true);
+   p.SetSize(width, mt); p.UseDevice(true);
+   g.SetSize(width, mt); g.UseDevice(true);
+   t.SetSize(width, mt); t.UseDevice(true);
+   r_bar.SetSize(width, mt); r_bar.UseDevice(true);
+   r_red.SetSize(width, mt); r_red.UseDevice(true);
+   g_red.SetSize(width, mt); g_red.UseDevice(true);
+}
+
+void BPCGSolver::Mult(const Vector &b, Vector &x) const
+{
+   int i;
+   double delta, delta0, del0;
+   double alpha, beta, gamma;
+
+   // Initialization
+   x.UseDevice(true);
+   if (iterative_mode)
+   {
+      oper->Mult(x, r);
+      subtract(b, r, r); // r = b - A x
+   }
+   else
+   {
+      r = b;
+      x = 0.0;
+   }
+
+   pprec->Mult(r,r_bar);  // r_bar = P r
+   p = r_bar;
+   oper->Mult(p, g);      // g = A p
+   oper->Mult(r_bar, t);  // t = A r_bar
+   iprec->Mult(r, r_red); // r_red = N r
+
+   delta = delta0 = Dot(t, r_red) - Dot(r_bar, r); // Dot(Pr, r)
+   if (delta0 >= 0.0) { initial_norm = sqrt(delta0); }
+   MFEM_ASSERT(IsFinite(delta), "norm = " << delta);
+   if (print_options.iterations || print_options.first_and_last)
+   {
+      mfem::out << "   Iteration : " << setw(3) << 0 << "  (P r, r) = "
+                << delta << (print_options.first_and_last ? " ...\n" : "\n");
+   }
+   Monitor(0, delta, r, x);
+
+   if (delta < 0.0)
+   {
+      if (print_options.warnings)
+      {
+         mfem::out << "BPCG: The preconditioner is not positive definite. (Pr, r) = "
+                   << delta << '\n';
+      }
+      converged = false;
+      final_iter = 0;
+      initial_norm = delta;
+      final_norm = delta;
+      return;
+   }
+   del0 = std::max(delta*rel_tol*rel_tol, abs_tol*abs_tol);
+   if (delta <= del0)
+   {
+      converged = true;
+      final_iter = 0;
+      final_norm = sqrt(delta);
+      return;
+   }
+
+   iprec->Mult(g, g_red);
+   gamma = Dot(g, g_red) - Dot(g,p); // Dot(Ap, p)
+   MFEM_ASSERT(IsFinite(gamma), "den (gamma) = " << gamma);
+   if (gamma <= 0.0)
+   {
+      if (Dot(r_bar, r_bar) > 0.0 && print_options.warnings)
+      {
+         mfem::out << "BPCG: The operator is not positive definite. (Ar, r) = "
+                   << gamma << '\n';
+      }
+      if (gamma == 0.0)
+      {
+         converged = false;
+         final_iter = 0;
+         final_norm = sqrt(delta);
+         return;
+      }
+   }
+
+   // Start iteration
+   converged = false;
+   final_iter = max_iter;
+   for (i = 1; true; )
+   {
+      alpha = delta0/gamma;
+      add(x,  alpha, p, x);  // x = x + alpha p
+      add(r, -alpha, g, r);  // r = r - alpha g
+
+      pprec->Mult(r, r_bar); // r_bar = P r
+      iprec->Mult(r, r_red); // r_red = N r
+      oper->Mult(r_bar, t);  // t = A r_bar
+      delta = Dot(t, r_red) - Dot(r_bar,r);
+
+      // Check
+      MFEM_ASSERT(IsFinite(delta), "norm = " << delta);
+      if (delta < 0.0)
+      {
+         if (print_options.warnings)
+         {
+            mfem::out << "BPCG: The preconditioner is not positive definite. (Pr, r) = "
+                      << delta << '\n';
+         }
+         converged = false;
+         final_iter = i;
+         break;
+      }
+      if (print_options.iterations)
+      {
+         mfem::out << "   Iteration : " << setw(3) << i << "  (Pr, r) = "
+                   << delta << std::endl;
+      }
+      Monitor(i, delta, r, x);
+      if (delta <= del0)
+      {
+         converged = true;
+         final_iter = i;
+         break;
+      }
+      if (++i > max_iter)
+      {
+         break;
+      }
+      // End check
+
+      beta = delta/delta0;
+      add(r_bar, beta, p, p); // p = r_bar + beta p
+      add(t, beta, g, g);     // g = t + beta g
+
+      delta0 = delta;
+      iprec->Mult(g, g_red);
+      gamma = Dot(g, g_red) - Dot(g,p); // Dot(Ap, p)
+      MFEM_ASSERT(IsFinite(gamma), "den (gamma) = " << gamma);
+      if (gamma <= 0.0)
+      {
+         if (Dot(r_bar, r_bar) > 0.0 && print_options.warnings)
+         {
+            mfem::out << "BPCG: The operator is not positive definite. (Ar, r) = "
+                      << gamma << '\n';
+         }
+         if (gamma == 0.0)
+         {
+            final_iter = i;
+            break;
+         }
+      }
+   }
+
+   if (print_options.first_and_last && !print_options.iterations)
+   {
+      mfem::out << "   Iteration : " << setw(3) << final_iter << "  (Pr, r) = "
+                << delta << '\n';
+   }
+   if (print_options.summary || (print_options.warnings && !converged))
+   {
+      mfem::out << "BPCG: Number of iterations: " << final_iter << '\n';
+   }
+   if (print_options.summary || print_options.iterations ||
+       print_options.first_and_last)
+   {
+      const auto arf = pow (gamma/delta0, 0.5/final_iter);
+      mfem::out << "Average reduction factor = " << arf << '\n';
+   }
+   if (print_options.warnings && !converged)
+   {
+      mfem::out << "BPCG: No convergence!" << '\n';
+   }
+
+   final_norm = sqrt(delta);
+   Monitor(final_iter, final_norm, r, x, true);
+}
+} // namespace blocksolvers
+} // namespace mfem
@@ -0,0 +1,189 @@
+// Copyright (c) 2023-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.
+//
+//          ----------------------------------------------------------
+//               Bramble-Pasciak preconditioning for Darcy problem
+//          ----------------------------------------------------------
+//
+// Main idea is to transform the block system
+//                 Ax = [ M  B^T ] [u] = [f] = b
+//                      [ B   0  ] [p] = [g]
+//     where:
+//        M = \int_\Omega (k u_h) \cdot v_h dx,
+//        B = -\int_\Omega (div_h u_h) q_h dx,
+//        u_h, v_h \in R_h (Raviart-Thomas finite element space),
+//        q_h \in W_h (piecewise discontinuous polynomials),
+// with a block transformation of the form X = A*N - Id
+//                  X = [ M*invQ - Id    0   ]
+//                      [     B*invQ    -Id  ]
+// where N is defined by
+//                  N = [ invQ    0 ]
+//                      [   0     0 ]
+// and Q is constructed such that Q and M-Q are both s.p.d..
+//
+// The solution x is then obtained by solving XAx = Xb with PCG as XA is s.p.d.
+//
+// The codes allows the user to provide such Q, or to construct it from the
+// element matrices M_T. Moreover, the user can provide a block preconditioner
+//                  P = [ M_0    0  ]
+//                      [  0    M_1 ]
+// for the transformed system XA.
+//
+// The code also allows the user to use BPCG, which is a special implementation
+// of the PCG iteration with the particular preconditioner H, defined as
+//                  H = [ M - Q    0  ]
+//                      [  0      M_1 ]
+// BPCG is efficient as it avoids the direct application of invH and X.
+
+#ifndef MFEM_BP_SOLVER_HPP
+#define MFEM_BP_SOLVER_HPP
+
+#include "darcy_solver.hpp"
+#include <memory>
+
+namespace mfem
+{
+namespace blocksolvers
+{
+
+/// Parameters for the BramblePasciakSolver method
+struct BPSParameters : IterSolveParameters
+{
+   bool use_bpcg = true;   // whether to use BPCG
+   double q_scaling = 0.5; // scaling (> 0 and < 1) of the Q preconditioner
+};
+
+/// Bramble-Pasciak Conjugate Gradient
+class BPCGSolver : public IterativeSolver
+{
+protected:
+   mutable Vector r, p, g, t, r_bar, r_red, g_red;
+   /// Remaining required operators
+   /*  Operator list
+    *  From IterativeSolver:
+    *  *oper  -> A  = [M, Bt; B, 0]
+    *  From this class:
+    *  *iprec -> N  = diag(M0, 0)
+    *  *pprec -> P' = diag(M0, M1) * [Id, 0; B*M0, -Id]
+    */
+   const Operator *iprec, *pprec;
+   void UpdateVectors();
+
+public:
+   BPCGSolver() { }
+   BPCGSolver(const Operator &ipc, const Operator &ppc) { pprec = &ppc; iprec = &ipc; }
+
+#ifdef MFEM_USE_MPI
+   BPCGSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
+   BPCGSolver(MPI_Comm comm_, const Operator &ipc, const Operator &ppc)
+      : IterativeSolver(comm_) { pprec = &ppc; iprec = &ipc; }
+#endif
+
+   virtual void SetOperator(const Operator &op)
+   { IterativeSolver::SetOperator(op); UpdateVectors(); }
+
+   virtual void SetPreconditioner(Solver &pc)
+   { if (Mpi::Root()) { MFEM_WARNING("SetPreconditioner does NO effect to BPCGSolver.\n"); } }
+
+   virtual void SetIncompletePreconditioner(const Operator &ipc)
+   { iprec = &ipc; }
+
+   virtual void SetParticularPreconditioner(const Operator &ppc)
+   { pprec = &ppc; }
+
+   virtual void Mult(const Vector &b, Vector &x) const;
+};
+
+/// Bramble-Pasciak Solver for Darcy equation.
+/** Bramble-Pasciak Solver for Darcy equation.
+ *  The basic idea is to precondition the mass matrix M with a s.p.d. matrix Q
+ *  such that M - Q remains s.p.d. Then we can transform the block operator into a
+ *  s.p.d. operator under a modified inner product.
+ *  In particular, this enable us to implement modified versions of CG iterations,
+ *  that rely on efficient applications of the required transformations.
+ *
+ *  We offer a mass preconditioner based on a rescalling of the diagonal of the
+ *  element mass matrices M_T.
+ *  We consider Q_T := alpha * lambda_min * D_T, where D_T := diag(M_T), and
+ *  lambda_min is the smallest eigenvalue of the following problem
+ *                M_T x = lambda * D_T x.
+ *  alpha is a parameter that is stricly between 0 and 1.
+ *
+ *  For more details, see:
+ *  1. Vassilevski, Multilevel Block Factorization Preconditioners (Appendix F.3),
+ *     Springer, 2008.
+ *  2. James H. Bramble and Joseph E. Pasciak.
+ *     A Preconditioning Technique for Indefinite Systems Resulting From Mixed
+ *     Approximations of Elliptic Problems. Mathematics of Computation, 50:1–17, 1988.
+ */
+class BramblePasciakSolver : public DarcySolver
+{
+   mutable bool use_bpcg;
+   std::unique_ptr<IterativeSolver> solver_;
+   BlockOperator *oop_, *ipc_;
+   ProductOperator *mop_;
+   SumOperator *map_;
+   ProductOperator *ppc_;
+   BlockDiagonalPreconditioner *cpc_;
+   std::unique_ptr<HypreParMatrix> M_;
+   std::unique_ptr<HypreParMatrix> B_;
+   std::unique_ptr<HypreParMatrix> Q_;
+   std::unique_ptr<HypreParMatrix> Q_e_;
+   OperatorPtr M0_;
+   OperatorPtr M1_;
+   Array<int> ess_zero_dofs_;
+
+   void Init(HypreParMatrix &M, HypreParMatrix &B,
+             HypreParMatrix &Q,
+             Solver &M0, Solver &M1,
+             const BPSParameters &param);
+public:
+   /// System and mass preconditioner are constructed from bilinear forms.
+   /// Elimination done internally.
+   BramblePasciakSolver(
+      const std::shared_ptr<ParBilinearForm> &mVarf,
+      const std::shared_ptr<ParMixedBilinearForm> &bVarf,
+      const Array<int>& ess_tdof_list,
+      const BPSParameters &param);
+
+   /// System and mass preconditioner are constructed from bilinear forms
+   BramblePasciakSolver(
+      const std::shared_ptr<ParBilinearForm> &mVarf,
+      const std::shared_ptr<ParMixedBilinearForm> &bVarf,
+      const BPSParameters &param);
+
+   /// System and mass preconditioner are user-provided
+   BramblePasciakSolver(
+      HypreParMatrix &M, HypreParMatrix &B, HypreParMatrix &Q,
+      Solver &M0, Solver &M1,
+      const BPSParameters &param);
+
+   /// Assemble a preconditioner for the mass matrix
+   /** Mass preconditioner corresponds to a local re-scaling
+    * based on the smallest eigenvalue of the generalized
+    * eigenvalue problem locally on each element T:
+    *         M_T x_T = lambda_T diag(M_T) x_T
+    * and we set Q_T = alpha * min(lambda_T) * diag(M_T),
+    * 0 < alpha < 1.
+   */
+   static HypreParMatrix *ConstructMassPreconditioner(ParBilinearForm &mVarf,
+                                                      double alpha = 0.5);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+   virtual void SetOperator(const Operator &op) { }
+   void SetEssZeroDofs(const Array<int>& dofs) { dofs.Copy(ess_zero_dofs_); }
+   virtual int GetNumIterations() const { return solver_->GetNumIterations(); }
+};
+
+} // namespace blocksolvers
+} // namespace mfem
+
+#endif // MFEM_BP_SOLVER_HPP
@@ -0,0 +1,88 @@
+// 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.
+
+#include "darcy_solver.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace blocksolvers;
+
+namespace mfem
+{
+namespace blocksolvers
+{
+void SetOptions(IterativeSolver& solver, const IterSolveParameters& param)
+{
+   solver.SetPrintLevel(param.print_level);
+   solver.SetMaxIter(param.max_iter);
+   solver.SetAbsTol(param.abs_tol);
+   solver.SetRelTol(param.rel_tol);
+}
+
+void DarcySolver::SetEliminatedSystems(std::shared_ptr<HypreParMatrix> M_e,
+                                       std::shared_ptr<HypreParMatrix> B_e,
+                                       const Array<int>& ess_tdof_list)
+{
+   M_e_ = M_e;
+   B_e_ = B_e;
+   rhs_needs_elimination_ = true;
+   ess_tdof_list.Copy(ess_tdof_list_);
+}
+
+void DarcySolver::EliminateEssentialBC(const Vector &ess_data,
+                                       Vector &rhs) const
+{
+   BlockVector blk_ess_data(ess_data.GetData(), offsets_);
+   BlockVector blk_rhs(rhs, offsets_);
+   M_e_->Mult(-1.0, blk_ess_data.GetBlock(0), 1.0, blk_rhs.GetBlock(0));
+   B_e_->Mult(-1.0, blk_ess_data.GetBlock(0), 1.0, blk_rhs.GetBlock(1));
+   for (int dof : ess_tdof_list_) { rhs[dof] = ess_data[dof]; }
+}
+
+/// Wrapper Block Diagonal Preconditioned MINRES (ex5p)
+/** Wrapper for assembling the discrete Darcy problem (ex5p)
+                     [ M  B^T ] [u] = [f]
+                     [ B   0  ] [p] = [g]
+**/
+BDPMinresSolver::BDPMinresSolver(HypreParMatrix& M, HypreParMatrix& B,
+                                 IterSolveParameters param)
+   : DarcySolver(M.NumRows(), B.NumRows()), op_(offsets_), prec_(offsets_),
+     BT_(B.Transpose()), solver_(M.GetComm())
+{
+   op_.SetBlock(0,0, &M);
+   op_.SetBlock(0,1, BT_.As<HypreParMatrix>());
+   op_.SetBlock(1,0, &B);
+
+   Vector Md;
+   M.GetDiag(Md);
+   BT_.As<HypreParMatrix>()->InvScaleRows(Md);
+   S_.Reset(ParMult(&B, BT_.As<HypreParMatrix>()));
+   BT_.As<HypreParMatrix>()->ScaleRows(Md);
+
+   prec_.SetDiagonalBlock(0, new HypreDiagScale(M));
+   prec_.SetDiagonalBlock(1, new HypreBoomerAMG(*S_.As<HypreParMatrix>()));
+   static_cast<HypreBoomerAMG&>(prec_.GetDiagonalBlock(1)).SetPrintLevel(0);
+   prec_.owns_blocks = true;
+
+   SetOptions(solver_, param);
+   solver_.SetOperator(op_);
+   solver_.SetPreconditioner(prec_);
+}
+
+void BDPMinresSolver::Mult(const Vector & x, Vector & y) const
+{
+   Vector x_e(x);
+   if (rhs_needs_elimination_) { EliminateEssentialBC(y, x_e);}
+   solver_.Mult(x_e, y);
+   for (int dof : ess_zero_dofs_) { y[dof] = 0.0; }
+}
+} // namespace blocksolvers
+} // namespace mfem
@@ -0,0 +1,73 @@
+// 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.
+
+#ifndef MFEM_DARCY_SOLVER_HPP
+#define MFEM_DARCY_SOLVER_HPP
+
+#include "mfem.hpp"
+#include <memory>
+#include <vector>
+
+namespace mfem
+{
+namespace blocksolvers
+{
+struct IterSolveParameters
+{
+   int print_level = 0;
+   int max_iter = 500;
+   double abs_tol = 1e-12;
+   double rel_tol = 1e-9;
+};
+
+void SetOptions(IterativeSolver& solver, const IterSolveParameters& param);
+
+/// Abstract solver class for Darcy's flow
+class DarcySolver : public Solver
+{
+protected:
+   Array<int> offsets_;
+   bool rhs_needs_elimination_;
+   Array<int> ess_tdof_list_;
+   std::shared_ptr<HypreParMatrix> M_e_;
+   std::shared_ptr<HypreParMatrix> B_e_;
+public:
+   DarcySolver(int size0, int size1)
+      : Solver(size0 + size1), offsets_(3), rhs_needs_elimination_(false)
+   { offsets_[0] = 0; offsets_[1] = size0; offsets_[2] = height; }
+   virtual int GetNumIterations() const = 0;
+   void SetEliminatedSystems(std::shared_ptr<HypreParMatrix> M_e,
+                             std::shared_ptr<HypreParMatrix> B_e,
+                             const Array<int>& ess_tdof_list);
+   void EliminateEssentialBC(const Vector &ess_data, Vector &rhs) const;
+};
+
+/// Wrapper for the block-diagonal-preconditioned MINRES defined in ex5p.cpp
+class BDPMinresSolver : public DarcySolver
+{
+   BlockOperator op_;
+   BlockDiagonalPreconditioner prec_;
+   OperatorPtr BT_;
+   OperatorPtr S_;   // S_ = B diag(M)^{-1} B^T
+   MINRESSolver solver_;
+   Array<int> ess_zero_dofs_;
+public:
+   BDPMinresSolver(HypreParMatrix& M, HypreParMatrix& B,
+                   IterSolveParameters param);
+   virtual void Mult(const Vector & x, Vector & y) const;
+   virtual void SetOperator(const Operator &op) { }
+   void SetEssZeroDofs(const Array<int>& dofs) { dofs.Copy(ess_zero_dofs_); }
+   virtual int GetNumIterations() const { return solver_.GetNumIterations(); }
+};
+} // namespace blocksolvers
+} // namespace mfem
+
+#endif // MFEM_DARCY_SOLVER_HPP
@@ -15,14 +15,10 @@ using namespace std;
 using namespace mfem;
 using namespace blocksolvers;

-void SetOptions(IterativeSolver& solver, const IterSolveParameters& param)
+namespace mfem
+{
+namespace blocksolvers
 {
-   solver.SetPrintLevel(param.print_level);
-   solver.SetMaxIter(param.max_iter);
-   solver.SetAbsTol(param.abs_tol);
-   solver.SetRelTol(param.rel_tol);
-}
-
 HypreParMatrix* TwoStepsRAP(const HypreParMatrix& Rt, const HypreParMatrix& A,
                            const HypreParMatrix& P)
 {
@@ -167,13 +163,14 @@ void DFSSpaces::CollectDFSData()

   GetP(data_.P_hcurl[level_], coarse_hcurl_fes_, *hcurl_fes_, true);

-   Vector trash1(hcurl_fes_->GetVSize()), trash2(hdiv_fes_->GetVSize());
   ParDiscreteLinearOperator curl(hcurl_fes_.get(), hdiv_fes_.get());
   curl.AddDomainInterpolator(new CurlInterpolator);
   curl.Assemble();
-   curl.EliminateTrialDofs(ess_bdr_attr_, trash1, trash2);
   curl.Finalize();
   data_.C[level_+1].Reset(curl.ParallelAssemble());
+   mfem::Array<int> ess_hcurl_tdof;
+   hcurl_fes_->GetEssentialTrueDofs(ess_bdr_attr_, ess_hcurl_tdof);
+   data_.C[level_+1].As<HypreParMatrix>()->EliminateCols(ess_hcurl_tdof);

   ++level_;

@@ -314,37 +311,6 @@ void SaddleSchwarzSmoother::Mult(const Vector & x, Vector & y) const
   blk_y.GetBlock(1) -= coarse_l2_projection;
 }

-BDPMinresSolver::BDPMinresSolver(HypreParMatrix& M, HypreParMatrix& B,
-                                 IterSolveParameters param)
-   : DarcySolver(M.NumRows(), B.NumRows()), op_(offsets_), prec_(offsets_),
-     BT_(B.Transpose()), solver_(M.GetComm())
-{
-   op_.SetBlock(0,0, &M);
-   op_.SetBlock(0,1, BT_.As<HypreParMatrix>());
-   op_.SetBlock(1,0, &B);
-
-   Vector Md;
-   M.GetDiag(Md);
-   BT_.As<HypreParMatrix>()->InvScaleRows(Md);
-   S_.Reset(ParMult(&B, BT_.As<HypreParMatrix>()));
-   BT_.As<HypreParMatrix>()->ScaleRows(Md);
-
-   prec_.SetDiagonalBlock(0, new HypreDiagScale(M));
-   prec_.SetDiagonalBlock(1, new HypreBoomerAMG(*S_.As<HypreParMatrix>()));
-   static_cast<HypreBoomerAMG&>(prec_.GetDiagonalBlock(1)).SetPrintLevel(0);
-   prec_.owns_blocks = true;
-
-   SetOptions(solver_, param);
-   solver_.SetOperator(op_);
-   solver_.SetPreconditioner(prec_);
-}
-
-void BDPMinresSolver::Mult(const Vector & x, Vector & y) const
-{
-   solver_.Mult(x, y);
-   for (int dof : ess_zero_dofs_) { y[dof] = 0.0; }
-}
-
 DivFreeSolver::DivFreeSolver(const HypreParMatrix &M, const HypreParMatrix& B,
                             const DFSData& data)
   : DarcySolver(M.NumRows(), B.NumRows()), data_(data), param_(data.param),
@@ -540,13 +506,16 @@ void DivFreeSolver::Mult(const Vector & x, Vector & y) const
   MFEM_VERIFY(x.Size() == offsets_[2], "MLDivFreeSolver: x size is invalid");
   MFEM_VERIFY(y.Size() == offsets_[2], "MLDivFreeSolver: y size is invalid");

-   if (ops_.Size() == 1) { smoothers_[0]->Mult(x, y); return; }
+   Vector x_e(x);
+   if (rhs_needs_elimination_) { EliminateEssentialBC(y, x_e);}
+
+   if (ops_.Size() == 1) { smoothers_[0]->Mult(x_e, y); return; }

   BlockVector blk_y(y, offsets_);

   BlockVector resid(offsets_);
   ops_.Last()->Mult(y, resid);
-   add(1.0, x, -1.0, resid, resid);
+   add(1.0, x_e, -1.0, resid, resid);

   BlockVector correction(offsets_);
   correction = 0.0;
@@ -573,7 +542,7 @@ void DivFreeSolver::Mult(const Vector & x, Vector & y) const
      ch.Start();

      ops_.Last()->Mult(y, resid);
-      add(1.0, x, -1.0, resid, resid);
+      add(1.0, x_e, -1.0, resid, resid);

      SolveDivFree(resid.GetBlock(0), correction.GetBlock(0));
      blk_y.GetBlock(0) += correction.GetBlock(0);
@@ -606,3 +575,6 @@ int DivFreeSolver::GetNumIterations() const
   }
   return solver_.As<IterativeSolver>()->GetNumIterations();
 }
+
+} // namespace blocksolvers
+} // namespace mfem
@@ -12,24 +12,12 @@
 #ifndef MFEM_DIVFREE_SOLVER_HPP
 #define MFEM_DIVFREE_SOLVER_HPP

-#include "mfem.hpp"
-#include <memory>
-#include <vector>
+#include "darcy_solver.hpp"

 namespace mfem
 {
 namespace blocksolvers
 {
-
-/// Parameters for iterative solver
-struct IterSolveParameters
-{
-   int print_level = 0;
-   int max_iter = 500;
-   double abs_tol = 1e-12;
-   double rel_tol = 1e-9;
-};
-
 /// Parameters for the divergence free solver
 struct DFSParameters : IterSolveParameters
 {
@@ -58,7 +46,7 @@ struct DFSData
   DFSParameters param;
 };

-/// Finite element spaces concerning divergence free solver.
+/// Finite element spaces concerning divergence free solvers
 /// The main usage of this class is to collect data needed for the solver.
 class DFSSpaces
 {
@@ -99,17 +87,7 @@ public:
   ParFiniteElementSpace* GetL2FES() const { return l2_fes_.get(); }
 };

-/// Abstract solver class for Darcy's flow
-class DarcySolver : public Solver
-{
-protected:
-   Array<int> offsets_;
-public:
-   DarcySolver(int size0, int size1) : Solver(size0 + size1), offsets_(3)
-   { offsets_[0] = 0; offsets_[1] = size0; offsets_[2] = height; }
-   virtual int GetNumIterations() const = 0;
-};
-
+/// Solvers for DFS
 /// Solver for B * B^T
 /// Compute the product B * B^T and solve it with CG preconditioned by BoomerAMG
 class BBTSolver : public Solver
@@ -180,26 +158,8 @@ public:
   virtual void SetOperator(const Operator &op) { }
 };

-/// Wrapper for the block-diagonal-preconditioned MINRES defined in ex5p.cpp
-class BDPMinresSolver : public DarcySolver
-{
-   BlockOperator op_;
-   BlockDiagonalPreconditioner prec_;
-   OperatorPtr BT_;
-   OperatorPtr S_;   // S_ = B diag(M)^{-1} B^T
-   MINRESSolver solver_;
-   Array<int> ess_zero_dofs_;
-public:
-   BDPMinresSolver(HypreParMatrix& M, HypreParMatrix& B,
-                   IterSolveParameters param);
-   virtual void Mult(const Vector & x, Vector & y) const;
-   virtual void SetOperator(const Operator &op) { }
-   void SetEssZeroDofs(const Array<int>& dofs) { dofs.Copy(ess_zero_dofs_); }
-   virtual int GetNumIterations() const { return solver_.GetNumIterations(); }
-};
-
-/** Divergence free solver.
-    The basic idea of the solver is to exploit a multilevel decomposition of
+/// Divergence free solver.
+/** The basic idea of the solver is to exploit a multilevel decomposition of
    Raviart-Thomas space to find a particular solution satisfying the divergence
    constraint, and then solve the remaining (divergence-free) component in the
    kernel space of the discrete divergence operator.
@@ -234,9 +194,7 @@ public:
   virtual void SetOperator(const Operator &op) { }
   virtual int GetNumIterations() const;
 };
-
 } // namespace blocksolvers
-
 } // namespace mfem

 #endif // MFEM_DIVFREE_SOLVER_HPP
@@ -21,7 +21,9 @@ CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)

-BLOCK_SOLVERS_SRC = div_free_solver.cpp block-solvers.cpp
+BLOCK_SOLVERS_SRC = block-solvers.cpp darcy_solver.cpp \
+                    bramble_pasciak.cpp block_hybridization.cpp \
+                    div_free_solver.cpp
 BLOCK_SOLVERS_OBJ = $(BLOCK_SOLVERS_SRC:.cpp=.o)

 SEQ_MINIAPPS = lor_solvers
@@ -51,7 +53,7 @@ plor_solvers.o: $(SRC)lor_mms.hpp
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(MFEM_LIBS)

 block-solvers: $(BLOCK_SOLVERS_OBJ)
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(BLOCK_SOLVERS_OBJ) $(MFEM_LIBS)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $^ $(MFEM_LIBS)

 %.o: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
@@ -63,9 +65,9 @@ include $(MFEM_TEST_MK)
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 block-solvers-test-par: block-solvers-constant block-solvers-anisotropic
 block-solvers-constant: block-solvers
-	@$(call mfem-test,$<, $(RUN_MPI), BlockSolver miniapp,-r 2)
+	@$(call mfem-test,$<, $(RUN_MPI), BlockSolver miniapp,-sr 1 -pr 1)
 block-solvers-anisotropic: block-solvers
-	@$(call mfem-test,$<, $(RUN_MPI), BlockSolver miniapp,-r 2\
+	@$(call mfem-test,$<, $(RUN_MPI), BlockSolver miniapp,-sr 1 -pr 1\
 	  -m $(SRC)anisotropic.mesh -c $(SRC)anisotropic.coeff\
 	  -eb $(SRC)anisotropic.bdr)
 lor_solvers-test-seq: lor_solvers
Author	SHA1	Message	Date
Barry Fadness	ec3e3e2db0	Print the parts of the dof_truedof matrix	2026-01-12 18:50:34 -08:00
Barry Fadness	eae8e49444	Change default command line argument values	2025-09-08 18:13:21 -07:00
Barry Fadness	037dc77bc6	Use different manufactured solution	2025-09-04 20:29:30 -07:00
Barry Fadness	bece01008c	Add stopwatch to time block-hybridization solver	2025-06-15 19:42:36 -07:00
Barry Fadness	f3490c0cab	Use automatic variable for prerequisites	2024-07-28 18:41:55 -07:00
Barry Fadness	df04a9235d	Change manufactured solution	2024-05-27 21:22:11 -07:00
Barry Fadness	33372c6154	Merge branch 'master' into block-hybridization-dev	2023-12-14 18:31:14 -08:00
Barry Fadness	7da4561e09	Remove redundant print statements	2023-10-30 12:31:27 -07:00
Barry Fadness	320b1f4022	Improve formatting consistency	2023-10-27 17:58:39 -07:00
Barry Fadness	c317677de2	Delete unncecessary pointer assignment	2023-10-27 17:33:33 -07:00
Barry Fadness	f983922aba	Fix mistake move GetElementVDofs out of else-statement	2023-10-27 17:32:36 -07:00
Barry Fadness	095ea3ab0f	Merge branch 'master' into block-hybridization-dev	2023-10-27 17:25:11 -07:00
Barry Fadness	62b3b06d33	Make pH owner of its factors	2023-10-25 00:02:50 -07:00
Barry Fadness	3080cba39d	Fix mistake move AddSubMatrix	2023-10-24 18:18:21 -07:00
Barry Fadness	1a77427557	Add MixedBilinearForm::ComputeElementMatrices() This new class member function is based on BilinearForm::ComputeElementMatrices(). It was also necessary to make changes to both MixedBilinearForm::Assemble() and MixedBilinearForm::ComputeElementMatrix().	2023-10-24 13:09:48 -07:00
Chak Shing Lee	a2941b0589	fix random initial values in ess_data_.Block(1)	2023-10-23 13:25:07 -07:00
Chak Shing Lee	d362363b6f	Merge branch 'block-hybridization-dev' of https://github.com/mfem/mfem into block-hybridization-dev	2023-10-23 11:09:30 -07:00
Chak Shing Lee	ebec0673ac	fix the way we eliminate Hcurl essential dofs	2023-10-23 11:09:12 -07:00
Barry Fadness	55f497f3db	Replace PtAP with two-step RAP	2023-10-16 13:45:10 -07:00
Barry Fadness	4f2d638993	Fix text output typo	2023-09-29 21:31:25 -07:00
Barry Fadness	6faef3873b	Extract block hybridization code	2023-09-20 17:29:59 -07:00
Gabriel Esteban Pinochet Soto	4b7fdcd10a	Update BCs in BramblerPasciakSolver - Create new constructor that ussses the Essential BC for defining the linear system and the preconditioner. - Make style	2023-09-14 15:25:39 -07:00
Gabriel Esteban Pinochet Soto	1c303ee013	Merge branch 'block-solvers/bramble-pasciak-dev' of github.com:mfem/mfem into block-hybridization-dev	2023-09-13 14:33:28 -07:00
Gabriel Esteban Pinochet Soto	1364914e7d	Small fix. Use MFEM_LAPACK macro.	2023-09-13 12:22:42 -07:00
Gabriel Esteban Pinochet Soto	84a248e42c	Add return nullptr	2023-09-13 11:55:33 -07:00
Gabriel Esteban Pinochet Soto	9d2c4c899c	Use MFEM_USE_LAPACK macro.	2023-09-13 10:30:04 -07:00
Gabriel Esteban Pinochet Soto	08b87fd9fe	Disable iterative mode in SumOperator	2023-09-12 14:15:12 -07:00
Gabriel Esteban Pinochet Soto	41131c3e0c	Remove unused varaible in BramblePasciakSolver	2023-09-12 13:42:36 -07:00
Gabriel Esteban Pinochet Soto	bfcb60ccec	Override SetPreconditioner.	2023-09-12 12:59:44 -07:00
Gabriel Esteban Pinochet Soto	26a6c3c735	Update tests options in makefiles.	2023-09-12 12:05:08 -07:00
Gabriel Esteban Pinochet Soto	0ec2ea20de	Make style.	2023-09-12 10:52:09 -07:00
Gabriel Esteban Pinochet Soto	1a4d4d4a61	Merge branch 'master' of github.com:mfem/mfem into block-solvers/bramble-pasciak-dev	2023-09-12 10:42:56 -07:00
Gabriel Esteban Pinochet Soto	ffb58713b1	Small change to SumOperator constructor. - To decide if we can add two iterative solvers.	2023-09-12 10:32:13 -07:00
Gabriel Esteban Pinochet Soto	c90d52d32b	Modify makefile. Minor fixes. - Move assert into ConstructMassPreconditioner - Update makefile	2023-09-11 17:58:16 -07:00
Gabriel Esteban Pinochet Soto	05ace32f42	Cleaner code. - Add description to README - Remove commentaries in BPCG - Move assert into ConstructMassPreconditioner	2023-09-11 11:21:27 -07:00
Chak Shing Lee	bca1919f39	fix typo in documentation	2023-09-06 14:57:35 -07:00
Chak Shing Lee	ac495b7490	Merge branch 'block-solvers/bramble-pasciak-dev' of https://github.com/mfem/mfem into block-hybridization-dev	2023-09-06 14:55:47 -07:00
Chak Shing Lee	fcc35e6c6f	remove one target in makefile	2023-09-06 14:54:28 -07:00
Chak Shing Lee	ed98b4bf4a	precompute element matrices of Mform_ in DarcyProblem	2023-09-06 14:42:20 -07:00
Chak Shing Lee	8f11035617	Merge branch 'block-solvers/bramble-pasciak-dev' of https://github.com/mfem/mfem into block-hybridization-dev	2023-09-06 14:21:10 -07:00
Chak Shing Lee	e60b446438	remove unused lines	2023-09-06 14:18:30 -07:00
Chak Shing Lee	ece621c512	resolve conflicts from merging with block-solvers/bramble-pasciak-dev	2023-09-06 14:16:50 -07:00
Chak Shing Lee	1253fe5b16	Merge branch 'block-hybridization-dev' of https://github.com/mfem/mfem into block-hybridization-dev	2023-09-06 11:54:29 -07:00
Chak Shing Lee	78d869d6fe	Merge branch 'csl/bilinearform-vdofs-fix' of https://github.com/mfem/mfem into block-solvers/bramble-pasciak-dev	2023-09-06 11:13:01 -07:00
Chak Shing Lee	2d8d46a208	use M instead of Q to construct the (1,1)-block preconditioner in BramblePasciakSolver. adjust documentation and print	2023-09-05 11:22:24 -07:00
Chak Shing Lee	40abb09c99	resolve conflicts with master	2023-09-01 16:58:24 -07:00
Chak Shing Lee	bb681b5172	Merge branch 'master' of https://github.com/mfem/mfem into block-hybridization-dev	2023-09-01 16:43:25 -07:00
Gabriel Esteban Pinochet Soto	bf9bdddfdc	Restructure of the code. - Remove overloaded BramblePasciakSolver::Init function - Move MFEM pre-defined BPSolver setting to the constructor - Add pointers to solvers in BramblePasciakSolver class - Simplify options - Remove redundant option in Parameter structs - Rename AddOperator to SumOperator - Wrap *.cpp files into namespace scopes - Add some documentation - Make style	2023-09-01 14:15:47 -07:00
Gabriel Esteban Pinochet Soto	2b88ef81ae	Add serial and parallel refinement flag.	2023-08-28 16:50:53 -07:00
Gabriel Esteban Pinochet Soto	5fc779e128	Reorganize functions - Move general-purpose functions to block-solvers.cpp - Reorganize code - TODO Serial and parallel refinement flags	2023-08-28 14:01:19 -07:00
Gabriel Esteban Pinochet Soto	2ca61c74b0	Remove warnings. - Remove warnings concerning unset/unused operators.	2023-08-22 17:29:46 -07:00
Gabriel Esteban Pinochet Soto	98e2af65da	Modify documentation. - Remove characters.	2023-08-21 10:34:54 -07:00
Gabriel Esteban Pinochet Soto	46a8bca7f5	Add documentation	2023-08-21 10:30:56 -07:00
Gabriel Esteban Pinochet Soto	cb95ddc422	Remove draft elasticity. Add documentation. - Remove draft. Comment CMake files. - Add some explanation regarding BPCG and BP transform. - Make style.	2023-08-21 10:12:52 -07:00
Gabriel Esteban Pinochet Soto	aae0f2b7f2	Compartmentalization of block-solvers. - Add bramble_pasciack.xpp, darcy_solver.xpp - Move solvers to their respective files and headders - Update makefile and CMakeList.txt - (WIP) Add draft elasticity solver (similar to darcy_solver) - (WIP) Add specific classes to handle the elasticity FES - (WIP) Define basic block structure	2023-08-17 17:55:16 -07:00
Gabriel Esteban Pinochet Soto	31b17b2913	Condense Init functions. - Define operator pointers when required. - Use OperatorPtr for solver_ in BPSClass. - Add MFEM Warnings. - Add H preconditioner.	2023-08-14 14:44:26 -07:00
Gabriel Esteban Pinochet Soto	6f3cf11d58	Indentation fix (make style).	2023-08-10 16:25:23 -07:00
Gabriel Esteban Pinochet Soto	0631cea458	Update BramblePasciakSolver::GetNumIterations.	2023-08-10 16:15:35 -07:00
Gabriel Esteban Pinochet Soto	79f2d77612	Add BPCGSolver. Modify BramblePasciakSolver class. - Add class BPCGSolver as a derived class of CGSolver. - Add bool use_bpcg in BramblePasciakSolver. - Storage main matrices in BramblePasciakSolver. - Remaining operators stored as pointers. - Update linalg/operator.hpp (make style). - TODO Fix bug final_iter.	2023-08-10 15:51:36 -07:00
Gabriel Esteban Pinochet Soto	08deeaacd9	Add AddOperator operator.	2023-08-07 10:35:26 -07:00
Gabriel Esteban Pinochet Soto	b4a0df6f5c	Add parameter control. Add some documentation. - Add parameter control for the default preconditioner. - Add brief explanation of the current class. - Include some relevant citations.	2023-08-07 10:29:28 -07:00
Chak Shing Lee	17271c5d27	fix doc	2023-08-02 18:57:15 -07:00
Chak Shing Lee	da5222a860	some minor adjustment, and renaming of variables to match with other DarcySolver	2023-08-02 18:42:47 -07:00
Gabriel Esteban Pinochet Soto	b765c24d1c	Check style with make style	2023-07-31 15:17:47 -07:00
Gabriel Esteban Pinochet Soto	3056cddd22	Add Bramble-Pasciak solver. - Modify DarcyProblem class. - Add Bramble-Pasciak class. - Define Bramble-Pasciak constructors. - Define Bramble-Pasciak Mult. - Update header accordingly.	2023-07-31 15:17:47 -07:00
Barry Fadness	e83a2b4dfa	Split Mult into member functions	2023-07-18 19:03:59 -07:00
Barry Fadness	aa6ee1e659	Split constructor into member functions	2023-07-17 23:13:02 -07:00
Barry Fadness	637de53e10	Change default boundary conditions	2023-02-12 18:34:48 -08:00
Barry Fadness	f347adafd7	Store and use mixed_dofs array	2023-02-12 18:13:35 -08:00
Barry Fadness	5cb5bdfa47	Make small style changes	2023-02-12 00:58:36 -08:00
Barry Fadness	b6b515a5cc	Remove unnecessary zero initializations	2023-02-12 00:57:10 -08:00
Barry Fadness	89c9ac4a3e	Move test_size init before use	2023-02-11 20:46:39 -08:00
Barry Fadness	2b238b3ccc	Reuse Minv_sub_vec vector	2023-02-11 20:42:58 -08:00
Barry Fadness	115cf9cc35	Move vector initialization outside of loops	2023-02-11 20:41:55 -08:00
Barry Fadness	030f6b0adf	Utilize offsets arrays	2023-02-11 20:37:40 -08:00
Barry Fadness	6ab54fae6a	Remove unnecessary test_dofs array	2023-02-11 20:32:44 -08:00
Barry Fadness	98f8c25102	Use existing function to restrict u hat	2023-02-11 19:51:22 -08:00
Barry Fadness	cb3649b0e7	Reuse rhs_r vector	2023-02-11 19:45:23 -08:00
Barry Fadness	8affea4ff7	Delete num_hat_dofs variable	2023-02-11 19:27:06 -08:00
Barry Fadness	52ce79500b	Make Ct_local code more consistent	2023-02-11 19:07:42 -08:00
Barry Fadness	f7b685a949	Initialize dense matrices outside of loop	2023-02-11 19:04:50 -08:00
Barry Fadness	561f9495c8	Correct the order of the contraint space If the FE space is RT_p, then GetOrder() returns p + 1. See fe_coll.hpp L360.	2023-02-03 20:32:37 -08:00
Barry Fadness	640502cd19	Use Last Array class member function	2023-02-03 19:16:08 -08:00
Barry Fadness	1e1524ef0a	Store the number of elements in const int	2023-02-03 19:07:52 -08:00
Barry Fadness	60f6b4138f	Use test_offsets array	2023-01-30 21:23:16 -08:00
Barry Fadness	8799a5a6c1	Move offsets code closer to Init	2023-01-30 20:37:12 -08:00
Barry Fadness	b6a7aa9e7e	Make lines more similar to hybridization.cpp	2023-01-29 21:24:46 -08:00
Barry Fadness	873fd364bd	Edit whitespace	2023-01-29 20:54:26 -08:00
Barry Fadness	8bdf7a5bde	Increase the value of epsilon	2023-01-29 20:52:55 -08:00
Barry Fadness	7586d62bee	Move creation of dofs array	2023-01-29 20:51:59 -08:00
Barry Fadness	ebafb82533	Remove get calls from for loops	2023-01-29 20:49:56 -08:00
Barry Fadness	de77303ca3	Move elimination step so that dofs are correct The element matrix elimination previously used the hat_dof_marker array. Now we need adjusted dofs instead, which are obtained from the signed dofs.	2023-01-29 17:29:18 -08:00
Barry Fadness	ffafd46e06	Remove the unnecessary hat_dof_marker array	2023-01-28 16:51:22 -08:00
Barry Fadness	5ffb1421b3	Move ess_dof_marker creation	2023-01-28 16:48:03 -08:00
Barry Fadness	bf0ebd3b4a	Remove skip_zeros variable and use default value	2023-01-28 16:44:58 -08:00
Barry Fadness	4483127d46	Remove num_hat_dofs variable in constructor	2023-01-28 16:44:10 -08:00
Barry Fadness	1c02f91d61	Fix bug in serial refinement number Each number in the fraction was an integer. If neither the numerator nor the denominator is a float, then only the quotient is returned. As a result the logarithm will be less than it should be, and the ceiling will not round up.	2022-12-23 21:14:21 -08:00
Barry Fadness	bc8d146070	Add description of block hybridization solver	2022-12-20 13:39:14 -08:00
Barry Fadness	d48a9fb683	Add brackets to free array memory in destructor	2022-12-20 13:38:34 -08:00
Barry Fadness	15666ae446	Fix small typo and edit preamble notes	2022-12-05 11:52:16 -08:00
Barry Fadness	c234398430	Prevent the mass coefficient from being destroyed An integrator does not own its coefficient. If a coefficient is created in a scoped region using stack memory, then it will expire even though a bilinear form may rely on it. Make the mass coefficient a unique pointer owned by the Darcy problem class. Then the coefficient will persist until a Darcy problem instance reaches the end of the main function.	2022-11-21 00:26:28 -08:00
Barry Fadness	842b8e81a5	Enable the hypre preconditioner again No error message is encountered for hypre versions < 22100.	2022-11-21 00:12:56 -08:00
Barry Fadness	95131ae6df	Fix dof true dof problems We need to go back and forth between dofs and true dofs. The offsets_ array should be true dofs, so some code could be deleted. All array sizes may be okay now.	2022-10-20 17:06:42 -07:00
Barry Fadness	7b7fd1d7d1	Fix essential bc elimination in right-hand side	2022-10-17 01:51:11 -07:00
Barry Fadness	51bf50bd1f	Edit comment and remove forgotten line	2022-10-17 01:49:17 -07:00
Barry Fadness	ba77f139b2	Simplify setting the right-hand side pressure block	2022-10-17 01:46:08 -07:00
Barry Fadness	c0db35f4be	Unset preconditioner for now because of seg faults	2022-09-29 13:04:54 -07:00
Barry Fadness	fad4206ccb	Clean up formatting to be consistent	2022-09-29 12:53:29 -07:00
Barry Fadness	8844cbbd87	Substitute offsets to avoid computing sizes multiple times	2022-09-29 12:31:01 -07:00
Barry Fadness	5da2f32570	Add essential bc elimination - incorrect	2022-09-29 12:23:47 -07:00
Barry Fadness	1d6b812c22	Change order of variable initialization to address warning	2022-09-26 00:34:04 -07:00
Barry Fadness	0c10422f95	Remove polynomial order from parameter list because unnecessary	2022-09-26 00:29:18 -07:00
Barry Fadness	836bc54397	Fix typo in anisotropic boundary file path	2022-09-26 00:06:56 -07:00
Barry Fadness	0fe70454b5	Change array name to know it is a list	2022-09-26 00:00:56 -07:00
Barry Fadness	3d420605fe	Allocate the preconditioner on the heap	2022-09-22 00:51:06 -07:00
Barry Fadness	f182173df4	Change tolerance value eps and delete print statement	2022-09-21 23:18:03 -07:00
Barry Fadness	feba2c4bfe	Set preconditioner and then operator (as recommended)	2022-09-21 23:15:18 -07:00
Barry Fadness	25f972fcdb	Convert between true dofs and local dofs	2022-09-21 23:09:50 -07:00
Barry Fadness	dd1d6f4e26	Fix obtaining the hat dofs for solution recovery	2022-09-21 23:07:57 -07:00
Barry Fadness	cd11ddf27a	Remove second memory allocation of private members	2022-09-21 22:59:21 -07:00
Barry Fadness	96f79c97b0	Make parallel mesh alias instead of a copy	2022-09-21 22:53:47 -07:00
Barry Fadness	d6574909ee	Make the constraint FE space a class member The constraint space owns its action-only prolongation operator. This operator is used in the Mult() function, so we must prevent the space from being destroyed after the constructor is called.	2022-09-21 22:04:54 -07:00
Barry Fadness	bed8da810d	Remove vector mass coefficient temporarily The vector mass piecewise constant coefficient created in the DarcyProblem class constructor is destroyed. It is needed in the hybridization solver to compute the element matrices. We will deal with this later.	2022-09-21 21:31:11 -07:00
Barry Fadness	4096ed429d	Make changes to match block-solvers and just to compile	2022-09-19 13:47:13 -07:00
Barry Fadness	513dd5914d	Add block hybridization solver class	2022-09-19 13:45:20 -07:00
Barry Fadness	c72c5f6570	Add block hybridization solver	2022-09-19 13:44:26 -07:00
Barry Fadness	b2146af63a	Add getter for parallel test FE space	2022-09-19 13:44:03 -07:00
Barry Fadness	e5f597a556	Cut and paste right-hand side code The code to reduce the right-hand side was intertwined with the construction of the left-hand side. Move the lines of code from the constructor to the Mult function.	2022-09-18 18:07:11 -07:00
Barry Fadness	1dd16d99f2	Copy and paste block hybridization code	2022-09-18 13:37:31 -07:00
chakshinglee	5f92c58dfe	make style	2022-09-12 22:14:29 -07:00
chakshinglee	ce8f6bb37d	shorter name for the block-diagonally-preconditioned MINRES solver	2022-09-12 21:57:11 -07:00
chakshinglee	26c86d7b8c	add the elimination step as an option for DarcySolver	2022-09-12 21:09:50 -07:00
chakshinglee	1fdc4b7dae	store both the right hand sides before and after elimination in DarcyProblem; add getter for the bilinear forms for hybridization solver	2022-09-12 18:02:03 -07:00
chakshinglee	3dfa563818	resolve conflicts with master	2022-09-12 16:16:03 -07:00
Lee	de57fe137e	Adding an access function to the serial constraint matrix (transpose) Ct in the Hybridization class. This will be useful in hybridizing block systems.	2017-06-12 09:50:58 -07:00
Lee	45167598a0	First commit. Adding functions to compute element matrices (only for domain integrators) for MixedBilinearForm. These are needed for forming local block systems.	2017-06-12 09:46:45 -07:00