Add comment and delete constant error epsilon

Add block hybridization solver. This version constructs
local facet matrices. Only the class constructor and
destructor have been implemented. The member function
Mult currently does nothing, but the reduced matrix
seems to be correct in both serial and parallel.

fem/pbilinearform.hpp

@@ -297,6 +297,8 @@ public:
       test_pfes  = test_fes;
    }
 
+   ParFiniteElementSpace *TestParFESpace() const { return test_pfes; }
+
    /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
    HypreParMatrix *ParallelAssemble();
 

miniapps/solvers/CMakeLists.txt

@@ -12,8 +12,10 @@
 if (MFEM_USE_MPI)
   add_mfem_miniapp(block-solvers
     MAIN block-solvers.cpp
-    EXTRA_SOURCES darcy_solver.cpp div_free_solver.cpp bramble_pasciak.cpp
-    EXTRA_HEADERS darcy_solver.hpp div_free_solver.hpp bramble_pasciak.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
@@ -43,7 +45,7 @@ if (MFEM_USE_MPI)
         $<TARGET_FILE:block-solvers> -rs 1 -rp 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})
     endif()

miniapps/solvers/block-solvers.cpp

@@ -14,11 +14,11 @@
 //          ----------------------------------------------------------
 //
 // This miniapp compares various linear solvers for the saddle point system
-// obtained from mixed finite element discretization of the simple mixed Darcy
-// problem in ex5p
+// obtained from the mixed finite element discretization of the simple mixed
+// Darcy problem in ex5p
 //
 //                            k*u + grad p = f
-//                           - div u      = g
+//                            - div u      = g
 //
 // with natural boundary condition -p = <given pressure>. We use a given exact
 // solution (u,p) and compute the corresponding r.h.s. (f,g). We discretize
@@ -29,30 +29,31 @@
 //    1. MINRES preconditioned by a block diagonal preconditioner
 //    2. The divergence free solver (couple and decoupled modes)
 //    3. The Bramble-Pasciak solver (using BPCG or regular PCG)
+//    4. A block hybridization solver
 //
 // 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 -rp 2 -se
 //    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. Each attribute may
+//        correspond to many elements. 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
-//        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.
+//        an essential boundary condition (on velocity 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 zero),
+//        then an essential (respectively a natural) boundary condition will be
+//        imposed on the part of the boundary with the i-th attribute.
 
 #include "mfem.hpp"
+#include "block_hybridization.hpp"
 #include "bramble_pasciak.hpp"
 #include "div_free_solver.hpp"
 #include <fstream>
@@ -64,11 +65,11 @@ 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);
-real_t p_exact(const Vector & x);
-void f_exact(const Vector & x, Vector & f);
-real_t g_exact(const Vector & x);
-real_t natural_bc(const Vector & x);
+void u_exact(const Vector &x, Vector &u);
+real_t p_exact(const Vector &x);
+void f_exact(const Vector &x, Vector &f);
+real_t g_exact(const Vector &x);
+real_t natural_bc(const Vector &x);
 
 /** Wrapper for assembling the discrete Darcy problem (ex5p)
                      [ M  B^T ] [u] = [f]
@@ -81,6 +82,7 @@ real_t natural_bc(const Vector & x);
        u_h, v_h \in R_h (Raviart-Thomas finite element space),
        q_h \in W_h (piecewise discontinuous polynomials),
        D: subset of the boundary where natural boundary condition is imposed. */
+
 class DarcyProblem
 {
    OperatorPtr M_;
@@ -95,7 +97,7 @@ class DarcyProblem
    VectorFunctionCoefficient ucoeff_;
    FunctionCoefficient pcoeff_;
    DFSSpaces dfs_spaces_;
-   PWConstCoefficient mass_coeff;
+   PWConstCoefficient mass_coeff_;
    const IntegrationRule *irs_[Geometry::NumGeom];
 public:
    DarcyProblem(Mesh &mesh, int num_refines, int order, const char *coef_file,
@@ -117,9 +119,9 @@ DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
                            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),
-     mass_coeff()
+     mass_coeff_()
 {
-   for (int l = 0; l < num_refs; l++)
+   for (int l = 0; l < num_refs; ++l)
    {
       mesh_.UniformRefinement();
       dfs_spaces_.CollectDFSData();
@@ -133,7 +135,7 @@ DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
       coef_vector.Load(coef_str, mesh.GetNE());
    }
 
-   mass_coeff.UpdateConstants(coef_vector);
+   mass_coeff_.UpdateConstants(coef_vector);
    VectorFunctionCoefficient fcoeff(mesh_.Dimension(), f_exact);
    FunctionCoefficient natcoeff(natural_bc);
    FunctionCoefficient gcoeff(g_exact);
@@ -157,7 +159,7 @@ DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
    bVarf_ = new ParMixedBilinearForm(dfs_spaces_.GetHdivFES(),
                                      dfs_spaces_.GetL2FES());
 
-   mVarf_->AddDomainIntegrator(new VectorFEMassIntegrator(mass_coeff));
+   mVarf_->AddDomainIntegrator(new VectorFEMassIntegrator(mass_coeff_));
    mVarf_->ComputeElementMatrices();
    mVarf_->Assemble();
    mVarf_->EliminateEssentialBC(ess_bdr, u_, fform);
@@ -174,7 +176,7 @@ DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
 
    rhs_.SetSize(M_->NumRows() + B_->NumRows());
    Vector rhs_block0(rhs_.GetData(), M_->NumRows());
-   Vector rhs_block1(rhs_.GetData()+M_->NumRows(), B_->NumRows());
+   Vector rhs_block1(rhs_.GetData() + M_->NumRows(), B_->NumRows());
    fform.ParallelAssemble(rhs_block0);
    gform.ParallelAssemble(rhs_block1);
 
@@ -183,17 +185,17 @@ DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
    Vector ess_data_block0(ess_data_.GetData(), M_->NumRows());
    u_.ParallelProject(ess_data_block0);
 
-   int order_quad = max(2, 2*order+1);
-   for (int i=0; i < Geometry::NumGeom; ++i)
+   int order_quad = max(2, 2 * order + 1);
+   for (int i = 0; i < Geometry::NumGeom; ++i)
    {
       irs_[i] = &(IntRules.Get(i, order_quad));
    }
 }
 
-void DarcyProblem::ShowError(const Vector& sol, bool verbose)
+void DarcyProblem::ShowError(const Vector &sol, bool verbose)
 {
    u_.Distribute(Vector(sol.GetData(), M_->NumRows()));
-   p_.Distribute(Vector(sol.GetData()+M_->NumRows(), B_->NumRows()));
+   p_.Distribute(Vector(sol.GetData() + M_->NumRows(), B_->NumRows()));
 
    real_t err_u  = u_.ComputeL2Error(ucoeff_, irs_);
    real_t norm_u = ComputeGlobalLpNorm(2, ucoeff_, mesh_, irs_);
@@ -205,7 +207,7 @@ void DarcyProblem::ShowError(const Vector& sol, bool verbose)
    mfem::out << "|| p_h - p_ex || / || p_ex || = " << err_p / norm_p << "\n";
 }
 
-void DarcyProblem::VisualizeSolution(const Vector& sol, string tag,
+void DarcyProblem::VisualizeSolution(const Vector &sol, string tag,
                                      int visport)
 {
    int num_procs, myid;
@@ -213,7 +215,7 @@ void DarcyProblem::VisualizeSolution(const Vector& sol, string tag,
    MPI_Comm_rank(mesh_.GetComm(), &myid);
 
    u_.Distribute(Vector(sol.GetData(), M_->NumRows()));
-   p_.Distribute(Vector(sol.GetData()+M_->NumRows(), B_->NumRows()));
+   p_.Distribute(Vector(sol.GetData() + M_->NumRows(), B_->NumRows()));
 
    const char vishost[] = "localhost";
    socketstream u_sock(vishost, visport);
@@ -229,7 +231,7 @@ void DarcyProblem::VisualizeSolution(const Vector& sol, string tag,
           << tag << " solver)'" << endl;
 }
 
-bool IsAllNeumannBoundary(const Array<int>& ess_bdr_attr)
+bool IsAllNeumannBoundary(const Array<int> &ess_bdr_attr)
 {
    for (int attr : ess_bdr_attr) { if (attr == 0) { return false; } }
    return true;
@@ -341,9 +343,9 @@ int main(int argc, char *argv[])
 
    // 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();
-   const DFSData& DFS_data = darcy.GetDFSData();
+   HypreParMatrix &M = darcy.GetM();
+   HypreParMatrix &B = darcy.GetB();
+   const DFSData &DFS_data = darcy.GetDFSData();
    delete mesh;
 
    if (Mpi::Root())
@@ -360,6 +362,10 @@ int main(int argc, char *argv[])
 
    // Setup various solvers for the discrete problem
    std::map<const DarcySolver*, real_t> setup_time;
+   chrono.Restart();
+   BlockHybridizationSolver bh(darcy.GetMform(), darcy.GetBform(), param);
+   setup_time[&bh] = chrono.RealTime();
+
    chrono.Restart();
    BDPMinresSolver bdp(M, B, param);
    setup_time[&bdp] = chrono.RealTime();
@@ -384,16 +390,17 @@ int main(int argc, char *argv[])
 
    std::map<const DarcySolver*, std::string> solver_to_name;
    solver_to_name[&bdp] = "Block-diagonal-preconditioned MINRES";
+   solver_to_name[&bh] = "Block hybridization";
    solver_to_name[&dfs_dm] = "Divergence free (decoupled mode)";
    solver_to_name[&dfs_cm] = "Divergence free (coupled mode)";
    solver_to_name[&bp_bpcg] = "Bramble Pasciak CG (using BPCG)";
    solver_to_name[&bp_pcg] = "Bramble Pasciak CG (using regular PCG)";
 
    // Solve the problem using all solvers
-   for (const auto& solver_pair : solver_to_name)
+   for (const auto &solver_pair : solver_to_name)
    {
-      auto& solver = solver_pair.first;
-      auto& name = solver_pair.second;
+      auto &solver = solver_pair.first;
+      auto &name = solver_pair.second;
 
       Vector sol = darcy.GetEssentialBC();
 
@@ -409,6 +416,7 @@ int main(int argc, char *argv[])
                    << setup_time[solver] + chrono.RealTime() << "s.\n"
                    << "   Iteration count: " << solver->GetNumIterations() <<"\n\n";
       }
+
       if (show_error && std::strcmp(coef_file, "") == 0)
       {
          darcy.ShowError(sol, Mpi::Root());
@@ -419,47 +427,51 @@ int main(int argc, char *argv[])
                    << "'.\nApproximation error is computed in this case!\n\n";
       }
 
-      if (visualization) { darcy.VisualizeSolution(sol, name, visport); }
-
+      if (visualization)
+      {
+         darcy.VisualizeSolution(sol, name, visport);
+      }
    }
-
    return 0;
 }
 
-void u_exact(const Vector & x, Vector & u)
+void u_exact(const Vector &x, Vector &u)
 {
    real_t xi(x(0));
-   real_t yi(x(1));
-   real_t zi(x.Size() == 3 ? x(2) : 0.0);
+   real_t xj(x(1));
+   real_t xk(x.Size() == 3 ? x(2) : 0.0);
 
-   u(0) = - exp(xi)*sin(yi)*cos(zi);
-   u(1) = - exp(xi)*cos(yi)*cos(zi);
+   u(0) = -exp(xi) * sin(xj) * cos(xk);
+   u(1) = -exp(xi) * cos(xj) * cos(xk);
    if (x.Size() == 3)
    {
-      u(2) = exp(xi)*sin(yi)*sin(zi);
+      u(2) = exp(xi) * sin(xj) * sin(xk);
    }
 }
 
-real_t p_exact(const Vector & x)
+real_t p_exact(const Vector &x)
 {
    real_t xi(x(0));
-   real_t yi(x(1));
-   real_t zi(x.Size() == 3 ? x(2) : 0.0);
-   return exp(xi)*sin(yi)*cos(zi);
+   real_t xj(x(1));
+   real_t xk(x.Size() == 3 ? x(2) : 0.0);
+   return exp(xi) * sin(xj) * cos(xk);
 }
 
-void f_exact(const Vector & x, Vector & f)
+void f_exact(const Vector &x, Vector &f)
 {
    f = 0.0;
 }
 
-real_t g_exact(const Vector & x)
+real_t g_exact(const Vector &x)
 {
-   if (x.Size() == 3) { return -p_exact(x); }
+   if (x.Size() == 3)
+   {
+      return -p_exact(x);
+   }
    return 0;
 }
 
-real_t natural_bc(const Vector & x)
+real_t natural_bc(const Vector &x)
 {
-   return (-p_exact(x));
+   return -p_exact(x);
 }

miniapps/solvers/block_hybridization.cpp

@@ -0,0 +1,288 @@
+#include "block_hybridization.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+namespace blocksolvers
+{
+
+/// Block hybridization solver
+BlockHybridizationSolver::BlockHybridizationSolver(
+   ParBilinearForm *mVarf,
+   ParMixedBilinearForm *bVarf,
+   IterSolveParameters param)
+   : DarcySolver(mVarf->ParFESpace()->GetTrueVSize(),
+                 bVarf->TestParFESpace()->GetTrueVSize()),
+     hdiv_space(mVarf->ParFESpace()),
+     l2_space(bVarf->TestParFESpace()),
+     solver_(hdiv_space->GetComm())
+{
+   ParMesh *mesh(hdiv_space->GetParMesh());
+   mesh->ExchangeFaceNbrData();
+   const int dim = mesh->Dimension();
+   const int num_elements(mesh->GetNE());
+
+   DG_Interface_FECollection fec(hdiv_space->FEColl()->GetOrder()-1,
+                                 mesh->Dimension());
+   multiplier_space = new ParFiniteElementSpace(mesh, &fec);
+   SparseMatrix reduced_matrix(multiplier_space->GetNDofs());
+
+   FaceElementTransformations *trans(nullptr);
+   NormalTraceIntegrator integ;
+   const FiniteElement *element(nullptr);
+
+   Table element_to_hdiv_dof = hdiv_space->GetElementToDofTable();
+   int *hdiv_offsets = element_to_hdiv_dof.GetI();
+
+   Table element_to_l2_dof = l2_space->GetElementToDofTable();
+   int *l2_offsets = element_to_l2_dof.GetI();
+
+   Array<int> hdiv_sizes(num_elements+1);
+   Array<int> l2_sizes(hdiv_sizes.Size());
+   Array<int> mixed_sizes(l2_sizes.Size());
+
+   hdiv_sizes[0] = 0;
+   l2_sizes[0] = 0;
+   mixed_sizes[0] = 0;
+
+   for (int i = 0; i < num_elements; ++i)
+   {
+      const int m = hdiv_offsets[i+1] - hdiv_offsets[i];
+      const int n = l2_offsets[i+1] - l2_offsets[i];
+      hdiv_sizes[i+1] = m * m;
+      hdiv_sizes[i+1] += hdiv_sizes[i];
+      l2_sizes[i+1] = n * n;
+      l2_sizes[i+1] += l2_sizes[i];
+      mixed_sizes[i+1] = m * n;
+      mixed_sizes[i+1] += mixed_sizes[i];
+   }
+
+   hdiv_data = new real_t[hdiv_sizes[num_elements]];
+   l2_data = new real_t[l2_sizes[num_elements]];
+   mixed_data = new real_t[mixed_sizes[num_elements]];
+   interior_indices = new Array<int>[num_elements];
+
+   // change the table name or use switch case?
+   Table element_to_facet_table;
+   if (2 == dim)
+   {
+      element_to_facet_table = mesh->ElementToEdgeTable();
+   }
+   else
+   {
+      element_to_facet_table = mesh->ElementToFaceTable();
+   }
+
+   // indirectly mark the interior faces using the be_to_face table
+   Array<int> interior_face_marker(mesh->GetNumFaces());
+   interior_face_marker = 1;
+   for (int i = 0; i < mesh->GetNBE(); ++i)
+   {
+      interior_face_marker[mesh->GetBdrElementFaceIndex(i)] = 0;
+   }
+
+   int *I = element_to_facet_table.GetI();
+   int *J = element_to_facet_table.GetJ();
+
+   // create element-to-interior-face table
+   // at some point really make this into a table
+   // instead of just arrays
+   int *nI = new int[num_elements+1];
+   nI[0] = 0;
+
+   // we already know the number of local interior faces summed over all elements
+   // this should still work for mixed meshes
+   int nnz(I[num_elements]-mesh->GetNBE());
+   int *nJ = new int[nnz];
+
+   int counter(0);
+   for (int elem_idx = 0; elem_idx < num_elements; ++elem_idx)
+   {
+      for (int j = I[elem_idx]; j < I[elem_idx+1]; ++j)
+      {
+         const int face_idx(J[j]);
+         if (interior_face_marker[face_idx]) // we check every face index usually twice
+                                             // can this be made more efficient?
+         {
+            nJ[counter++] = face_idx;
+         }
+      }
+      nI[elem_idx+1] = counter;
+   }
+
+   const Table &face_to_dof_table(multiplier_space->GetFaceToDofTable());
+   const int *faceI = face_to_dof_table.GetI();
+   const int *faceJ = face_to_dof_table.GetJ();
+
+   int *interior_faceI = new int[faceI[mesh->GetNumFaces()]];
+   interior_faceI[0] = 0;
+   int total_dofs(0);
+   for (int elem_idx = 0; elem_idx < num_elements; ++elem_idx)
+   {
+      for (int k = nI[elem_idx]; k < nI[elem_idx+1]; ++k)
+      {
+         const int face_idx = nJ[k];
+         total_dofs += faceI[face_idx+1] - faceI[face_idx];
+      }
+      interior_faceI[elem_idx+1] = total_dofs;
+   }
+
+   int *interior_faceJ = new int[interior_faceI[num_elements]];
+   counter = 0;
+   for (int elem_idx = 0; elem_idx < num_elements; ++elem_idx)
+   {
+      for (int k = nI[elem_idx]; k < nI[elem_idx+1]; ++k)
+      {
+         const int face_idx = nJ[k];
+         for (int dof_idx = faceI[face_idx]; dof_idx < faceI[face_idx+1]; ++dof_idx)
+         {
+            interior_faceJ[counter++] = faceJ[dof_idx];
+         }
+      }
+   }
+
+   for (int elem_idx = 0; elem_idx < num_elements; ++elem_idx)
+   {
+      const int m = hdiv_offsets[elem_idx+1] - hdiv_offsets[elem_idx];
+      DenseMatrix hdiv_matrix(hdiv_data+hdiv_sizes[elem_idx], m, m);
+      mVarf->ComputeElementMatrix(elem_idx, hdiv_matrix);
+      hdiv_matrix.Invert(); // overwrite hdiv_matrix
+
+      const int n = l2_offsets[elem_idx+1] - l2_offsets[elem_idx];
+      DenseMatrix mixed_matrix(mixed_data+mixed_sizes[elem_idx], n, m);
+      bVarf->ComputeElementMatrix(elem_idx, mixed_matrix);
+
+      DenseMatrix product_matrix(mixed_matrix.Height(), hdiv_matrix.Width());
+      mfem::Mult(mixed_matrix, hdiv_matrix, product_matrix); // BA^{-1}
+
+      DenseMatrix l2_matrix(l2_data+l2_sizes[elem_idx], n, n);
+      MultABt(product_matrix, mixed_matrix, l2_matrix); // BA^{-1}B^t = -S
+      l2_matrix.Invert(); // (BA^{-1}B^t)^{-1} = -S^{-1}
+      l2_matrix.Neg(); // -(BA^{-1}B^t)^{-1} = S^{-1}
+
+      mfem::Mult(l2_matrix, product_matrix, mixed_matrix); // overwrite mixed matrix with S^{-1}BA^{-1}
+      DenseMatrix temp_matrix(product_matrix.Width(), mixed_matrix.Width());
+      MultAtB(product_matrix, mixed_matrix, temp_matrix); // A^{-T}B^TS^{-1}BA^{-1}
+      hdiv_matrix += temp_matrix;
+      mixed_matrix.Neg(); // -S^{-1}BA^{-1}
+
+      element = hdiv_space->GetFE(elem_idx);
+      const FiniteElement *face(nullptr);
+
+      // exploit that dense matrices are in column-major format
+      // we do not copy each face matrix into the boundary matrix
+      const int num_multiplier_dofs = interior_faceI[elem_idx+1] - interior_faceI[elem_idx];
+      Array<int> multiplier_dofs;
+      multiplier_dofs.MakeRef(interior_faceJ+interior_faceI[elem_idx], num_multiplier_dofs);
+      DenseMatrix bdr_matrix(hdiv_matrix.Height(), num_multiplier_dofs);
+      DenseMatrix face_matrix;
+      int data_offset = 0;
+      const int num_faces = nI[elem_idx+1] - nI[elem_idx];
+      const int height = bdr_matrix.Height();
+      const int width = bdr_matrix.Width() / num_faces; // number of dofs per face
+
+      for (int k = nI[elem_idx]; k < nI[elem_idx+1]; ++k)
+      {
+         const int face_idx = nJ[k];
+         trans = mesh->GetFaceElementTransformations(face_idx);
+         face = multiplier_space->GetFaceElement(face_idx);
+         face_matrix.UseExternalData(bdr_matrix.Data()+data_offset, height, width);
+         integ.AssembleTraceFaceMatrix(elem_idx, *face, *element, *trans,
+                                       face_matrix);
+         data_offset += height*width;
+      }
+
+      temp_matrix.SetSize(height, bdr_matrix.Width());
+      mfem::Mult(hdiv_matrix, bdr_matrix, temp_matrix);
+      product_matrix.SetSize(bdr_matrix.Width());
+      MultAtB(bdr_matrix, temp_matrix, product_matrix);
+      reduced_matrix.AddSubMatrix(multiplier_dofs, multiplier_dofs, product_matrix);
+
+   /*
+      Array<int> face_indices_array;
+      element_to_facet_table.GetRow(elem_idx, face_indices_array);
+      DenseMatrix *face_matrices = new DenseMatrix[face_indices_array.Size()];
+      Array<int> *face_dofs = new Array<int>[face_indices_array.Size()];
+
+      for (int local_index = 0; local_index < face_indices_array.Size();
+           ++local_index)
+      {
+         const int face_index(face_indices_array[local_index]);
+         trans = mesh->GetFaceElementTransformations(face_index);
+         if (!mesh->FaceIsTrueInterior(face_index))
+         {
+            continue;
+         }
+         interior_indices[elem_idx].Append(local_index);
+         multiplier_space->GetFaceVDofs(face_index, face_dofs[local_index]);
+         face = multiplier_space->GetFaceElement(face_index);
+         integ.AssembleTraceFaceMatrix(elem_idx, *face, *element, *trans,
+                                       face_matrices[local_index]);
+      }
+      for (int column_index : interior_indices[elem_idx])
+      {
+         DenseMatrix temp_matrix(saved_hdiv_matrices[elem_idx].Height(),
+                                 face_matrices[column_index].Width());
+         mfem::Mult(saved_hdiv_matrices[elem_idx], face_matrices[column_index],
+                    temp_matrix);
+         for (int row_index : interior_indices[elem_idx])
+         {
+            DenseMatrix product_matrix(face_matrices[row_index].Width(),
+                                       temp_matrix.Width());
+            MultAtB(face_matrices[row_index], temp_matrix, product_matrix);
+            reduced_matrix.AddSubMatrix(face_dofs[row_index], face_dofs[column_index],
+                                        product_matrix);
+         }
+      }
+      delete []face_dofs;
+      delete []face_matrices;
+   */
+   }
+   delete []interior_faceI;
+   delete []interior_faceJ;
+   delete []nI;
+   delete []nJ;
+
+   reduced_matrix.Finalize(1, true);
+
+   HypreParMatrix *P(multiplier_space->Dof_TrueDof_Matrix());
+   HypreParMatrix *dH = new HypreParMatrix(multiplier_space->GetComm(),
+                                           multiplier_space->GlobalVSize(),
+                                           multiplier_space->GetDofOffsets(), &reduced_matrix);
+   HypreParMatrix *dHP = ParMult(dH, P);
+   HypreParMatrix *Pt(P->Transpose());
+   pH = ParMult(Pt, dHP, true);
+
+   delete dH;
+   /*
+      OperatorPtr pP(Operator::Hypre_ParCSR);
+      pP.ConvertFrom(multiplier_space->Dof_TrueDof_Matrix());
+      OperatorPtr dH(pP.Type());
+      dH.MakeSquareBlockDiag(multiplier_space->GetComm(), multiplier_space->GlobalVSize(),
+                             multiplier_space->GetDofOffsets(), &reduced_matrix);
+      OperatorPtr AP(ParMult(dH.As<HypreParMatrix>(), pP.As<HypreParMatrix>()));
+      OperatorPtr R(pP.As<HypreParMatrix>()->Transpose());
+      pH = ParMult(R.As<HypreParMatrix>(), AP.As<HypreParMatrix>(), true);
+   */
+   preconditioner = new HypreBoomerAMG(*pH);
+   preconditioner->SetPrintLevel(0);
+
+   SetOptions(solver_, param);
+   solver_.SetPreconditioner(*preconditioner);
+   solver_.SetOperator(*pH);
+}
+
+BlockHybridizationSolver::~BlockHybridizationSolver()
+{
+   delete pH;
+   delete []hdiv_data;
+   delete []l2_data;
+   delete []mixed_data;
+   delete []interior_indices;
+   delete preconditioner;
+   delete multiplier_space;
+}
+
+} // namespace blocksolvers
+} // namespace mfem

miniapps/solvers/block_hybridization.hpp

@@ -0,0 +1,40 @@
+#ifndef MFEM_BH_SOLVER_HPP
+#define MFEM_BH_SOLVER_HPP
+
+#include "darcy_solver.hpp"
+
+namespace mfem
+{
+namespace blocksolvers
+{
+
+class BlockHybridizationSolver : public DarcySolver
+{
+   ParFiniteElementSpace *hdiv_space;
+   ParFiniteElementSpace *l2_space;
+   ParFiniteElementSpace *multiplier_space;
+
+   CGSolver solver_;
+
+   real_t *hdiv_data;
+   real_t *l2_data;
+   real_t *mixed_data;
+   Array<int> *interior_indices;
+
+   HypreParMatrix *pH;
+   HypreBoomerAMG *preconditioner;
+
+public:
+   BlockHybridizationSolver(ParBilinearForm *mVarf,
+                            ParMixedBilinearForm *bVarf,
+                            IterSolveParameters param);
+   ~BlockHybridizationSolver();
+   void Mult(const Vector &x, Vector&y) const override { }
+   void SetOperator(const Operator &op) override { }
+   int GetNumIterations() const override { return solver_.GetNumIterations(); }
+};
+
+} // namespace blocksolvers
+} // namespace mfem
+
+#endif // MFEM_BH_SOLVER_HPP

miniapps/solvers/darcy_solver.cpp

@@ -17,7 +17,7 @@ namespace mfem
 {
 namespace blocksolvers
 {
-void SetOptions(IterativeSolver& solver, const IterSolveParameters& param)
+void SetOptions(IterativeSolver &solver, const IterSolveParameters &param)
 {
    solver.SetPrintLevel(param.print_level);
    solver.SetMaxIter(param.max_iter);
@@ -30,8 +30,8 @@ void SetOptions(IterativeSolver& solver, const IterSolveParameters& param)
                      [ M  B^T ] [u] = [f]
                      [ B   0  ] [p] = [g]
 */
-BDPMinresSolver::BDPMinresSolver(const HypreParMatrix& M,
-                                 const HypreParMatrix& B,
+BDPMinresSolver::BDPMinresSolver(const HypreParMatrix &M,
+                                 const HypreParMatrix &B,
                                  IterSolveParameters param)
    : DarcySolver(M.NumRows(), B.NumRows()), op_(offsets_), prec_(offsets_),
      BT_(B.Transpose()), solver_(M.GetComm())
@@ -56,7 +56,7 @@ BDPMinresSolver::BDPMinresSolver(const HypreParMatrix& M,
    solver_.SetPreconditioner(prec_);
 }
 
-void BDPMinresSolver::Mult(const Vector & x, Vector & y) const
+void BDPMinresSolver::Mult(const Vector &x, Vector &y) const
 {
    solver_.Mult(x, y);
    for (int dof : ess_zero_dofs_) { y[dof] = 0.0; }

miniapps/solvers/darcy_solver.hpp

@@ -37,7 +37,7 @@ struct IterSolveParameters
 #endif
 };
 
-void SetOptions(IterativeSolver& solver, const IterSolveParameters& param);
+void SetOptions(IterativeSolver &solver, const IterSolveParameters &param);
 
 /// Abstract solver class for Darcy's flow
 class DarcySolver : public Solver
@@ -60,10 +60,10 @@ class BDPMinresSolver : public DarcySolver
    MINRESSolver solver_;
    Array<int> ess_zero_dofs_;
 public:
-   BDPMinresSolver(const HypreParMatrix& M,
-                   const HypreParMatrix& B,
+   BDPMinresSolver(const HypreParMatrix &M,
+                   const HypreParMatrix &B,
                    IterSolveParameters param);
-   void Mult(const Vector & x, Vector & y) const override;
+   void Mult(const Vector &x, Vector &y) const override;
    void SetOperator(const Operator &op) override { }
    void SetEssZeroDofs(const Array<int>& dofs) { dofs.Copy(ess_zero_dofs_); }
    int GetNumIterations() const override { return solver_.GetNumIterations(); }

miniapps/solvers/makefile

@@ -20,7 +20,8 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-BLOCK_SOLVERS_SRC = block-solvers.cpp darcy_solver.cpp bramble_pasciak.cpp div_free_solver.cpp
+BLOCK_SOLVERS_SRC = block-solvers.cpp darcy_solver.cpp bramble_pasciak.cpp \
+                    div_free_solver.cpp block_hybridization.cpp
 BLOCK_SOLVERS_OBJ = $(BLOCK_SOLVERS_SRC:.cpp=.o)
 
 LOR_ELAST_SRC = block_fespace_operator.cpp lor_elast.cpp