Merge branch 'dg-agglom' of github.com:mfem/mfem into dg-agglom

.gitignore

@@ -373,6 +373,9 @@ miniapps/solvers/ParaView
 miniapps/solvers/mesh.*
 miniapps/solvers/sol.*
 
+miniapps/dg-agglomeration/ParaView
+miniapps/dg-agglomeration/dg_agglom
+
 miniapps/hdiv-linear-solver/darcy
 miniapps/hdiv-linear-solver/grad_div
 

linalg/solvers.cpp

@@ -780,7 +780,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
    if (print_options.iterations || print_options.first_and_last)
    {
       mfem::out << "   Iteration : " << setw(3) << 0 << "  (B r, r) = "
-                << nom << (print_options.first_and_last ? " ...\n" : "\n");
+                << r.Norml2() << (print_options.first_and_last ? " ...\n" : "\n");
    }
 
    if (nom < 0.0)
@@ -864,7 +864,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
       if (print_options.iterations)
       {
          mfem::out << "   Iteration : " << setw(3) << i << "  (B r, r) = "
-                   << betanom << std::endl;
+                   << r.Norml2() << std::endl;
       }
 
       if (Monitor(i, betanom, r, x) || betanom <= r0)
@@ -909,7 +909,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
    if (print_options.first_and_last && !print_options.iterations)
    {
       mfem::out << "   Iteration : " << setw(3) << final_iter << "  (B r, r) = "
-                << betanom << '\n';
+                << r.Norml2() << '\n';
    }
    if (print_options.summary || (print_options.warnings && !converged))
    {
@@ -2894,20 +2894,329 @@ void MinimumDiscardedFillOrdering(SparseMatrix &C, Array<int> &p)
    }
 }
 
-BlockILU::BlockILU(int block_size_,
+Blockl1Jacobi::Blockl1Jacobi(const Operator &op, int block_size_, real_t damping_)
+   : Solver(0),
+     block_size(block_size_),
+     damping(damping_)
+{
+   SetOperator(op);
+}
+
+void Blockl1Jacobi::SetOperator(const Operator &op)
+{
+   const SparseMatrix *A = NULL;
+   if (A == NULL)
+   {
+      A = dynamic_cast<const SparseMatrix *>(&op);
+      if (A == NULL)
+      {
+         MFEM_ABORT("Blockl1Jacobi must be created with a SparseMatrix or HypreParMatrix");
+      }
+   }
+   height = op.Height();
+   width = op.Width();
+   MFEM_VERIFY(A->Finalized(), "Matrix must be finalized.");
+   GetDiagonalBlocks(*A);
+}
+
+void Blockl1Jacobi::GetDiagonalBlocks(const SparseMatrix &A)
+{
+   if (A.Height() % block_size != 0)
+   {
+      MFEM_ABORT("Blockl1Jacobi: block size must evenly divide the matrix size");
+   }
+
+   const int nrows = A.Height();
+   const int nblockrows = nrows / block_size;
+
+   const int *I = A.HostReadI();
+   const int *J = A.HostReadJ();
+   const real_t *V = A.HostReadData();
+
+   DB.SetSize(block_size, block_size, nblockrows);
+   DB = 0.0;
+   real_t row_sum;
+   ipiv.SetSize(block_size*nblockrows);
+
+   for (int iblock = 0; iblock < nblockrows; ++iblock)
+   {
+      for (int bi = 0; bi < block_size; ++bi)
+      {
+         int i = iblock * block_size + bi;
+         row_sum = 0;
+         for (int k = I[i]; k < I[i+1]; ++k)
+         {
+            const int j = J[k];
+            real_t val = V[k];
+            if (j >= iblock*block_size && j < (iblock + 1)*block_size)
+            {
+               const int bj = j - iblock*block_size;
+               DB(bi, bj, iblock) = val;
+            }
+            row_sum = row_sum+val;
+         }
+         DB(bi, bi, iblock) += row_sum;
+      }
+      LUFactors factorization(DB.GetData(iblock), &ipiv[iblock*block_size]);
+      factorization.Factor(block_size);
+   }
+}
+
+void Blockl1Jacobi::Mult(const Vector &b, Vector &x) const
+{
+   MFEM_VERIFY(height > 0, "Blockl1Jacobi preconditioner is not constructed");
+   const int nblockrows = Height()/block_size;
+   Vector xi;
+
+   for (int i = 0; i<nblockrows; ++i)
+   {
+      xi.SetDataAndSize(&x[i*block_size], block_size);
+      for (int ib=0; ib<block_size; ++ib)
+      {
+         xi[ib] = b[ib + i*block_size];
+      }
+      // x_i = D_ii^{-1} x_i
+      LUFactors A_ii_inv(DB.GetData(i), &ipiv[i*block_size]);
+      A_ii_inv.Solve(block_size, 1, xi.GetData());
+   }
+   x *= damping;
+}
+
+void Blockl1Jacobi::MultTranspose(const Vector &b, Vector &x) const
+{
+   MFEM_VERIFY(height > 0, "Blockl1Jacobi preconditioner is not constructed");
+   const int nblockrows = Height()/block_size;
+
+   Vector xi;
+   for (int i = nblockrows-1; i>=0; --i)
+   {
+      xi.SetDataAndSize(&x[i*block_size], block_size);
+      for (int ib=0; ib<block_size; ++ib)
+      {
+         xi[ib] = b[ib + i*block_size];
+      }
+      // x_i = D_ii^{-1} x_i
+      LUFactors A_ii_inv(DB.GetData(i), &ipiv[i*block_size]);
+      A_ii_inv.Solve(block_size, 1, xi.GetData());
+   }
+   x *= damping;
+}
+
+BlockJacobi::BlockJacobi(const Operator &op, SparseMatrix &Block_Diag_Mat_)
+   : Solver(0),
+     Block_Diag_Mat(Block_Diag_Mat_)
+{
+   SetOperator(op);
+}
+
+void BlockJacobi::SetOperator(const Operator &op)
+{
+   const SparseMatrix *A = NULL;
+   if (A == NULL)
+   {
+      A = dynamic_cast<const SparseMatrix *>(&op);
+      if (A == NULL)
+      {
+         MFEM_ABORT("BlockJacobi must be created with a SparseMatrix or HypreParMatrix");
+      }
+   }
+   height = op.Height();
+   width = op.Width();
+   MFEM_VERIFY(A->Finalized(), "Matrix must be finalized.");
+   // CreateBlockPatternAndFactorize(*A);
+}
+
+
+void BlockJacobi::Mult(const Vector &b, Vector &x) const
+{
+   MFEM_VERIFY(height > 0, "BlockJacobi preconditioner is not constructed");
+   Block_Diag_Mat.Mult(b, x);
+}
+
+void BlockJacobi::MultTranspose(const Vector &b, Vector &x) const
+{
+   MFEM_VERIFY(height > 0, "BlockJacobi preconditioner is not constructed");
+   Block_Diag_Mat.MultTranspose(b, x);
+}
+
+BlockGS::BlockGS(const Operator &op, int block_size_, real_t damping_)
+   : Solver(0),
+     block_size(block_size_),
+     damping(damping_)
+{
+   SetOperator(op);
+}
+
+
+void BlockGS::SetOperator(const Operator &op)
+{
+   const SparseMatrix *A = NULL;
+   if (A == NULL)
+   {
+      A = dynamic_cast<const SparseMatrix *>(&op);
+      if (A == NULL)
+      {
+         MFEM_ABORT("BlockGS must be created with a SparseMatrix or HypreParMatrix");
+      }
+   }
+   height = op.Height();
+   width = op.Width();
+   MFEM_VERIFY(A->Finalized(), "Matrix must be finalized.");
+   CreateBlockPatternAndFactorize(*A);
+}
+
+void BlockGS::CreateBlockPatternAndFactorize(const SparseMatrix &A)
+{
+   if (A.Height() % block_size != 0)
+   {
+      MFEM_ABORT("BlockGS: block size must evenly divide the matrix size");
+   }
+
+   const int nrows = A.Height();
+   const int nblockrows = nrows / block_size;
+
+   const int *I = A.HostReadI();
+   const int *J = A.HostReadJ();
+   const real_t *V = A.HostReadData();
+
+   int nnz = 0;
+
+   std::vector<std::set<int>> unique_block_cols(nblockrows);
+
+   for (int iblock = 0; iblock < nblockrows; ++iblock)
+   {
+      for (int bi = 0; bi < block_size; ++bi)
+      {
+         int i = iblock * block_size + bi;
+         for (int k = I[i]; k < I[i + 1]; ++k)
+         {
+            unique_block_cols[iblock].insert(J[k] / block_size);
+         }
+      }
+      nnz += static_cast<int>(unique_block_cols[iblock].size());
+   }
+
+   ID.SetSize(nblockrows);
+   IB.SetSize(nblockrows + 1);
+   IB[0] = 0;
+   JB.SetSize(nnz);
+   AB.SetSize(block_size, block_size, nnz);
+   DB.SetSize(block_size, block_size, nblockrows);
+   AB = 0.0;
+   DB = 0.0;
+   ipiv.SetSize(block_size*nblockrows);
+   int counter = 0;
+   for (int iblock = 0; iblock < nblockrows; ++iblock)
+   {
+      for (int jblock : unique_block_cols[iblock])
+      {
+         JB[counter] = jblock;
+         if (iblock == jblock)
+         {
+            ID[iblock] = counter;
+         }
+         for (int bi = 0; bi < block_size; ++bi)
+         {
+            const int i = iblock*block_size + bi;
+            for (int k = I[i]; k < I[i + 1]; ++k)
+            {
+               const int j = J[k];
+               if (j >= jblock*block_size && j < (jblock + 1)*block_size)
+               {
+                  const int bj = j - jblock*block_size;
+                  real_t val = V[k];
+                  AB(bi, bj, counter) = val;
+                  // Extract the diagonal
+                  if (iblock == jblock)
+                  {
+                     DB(bi, bj, iblock) = val;
+                  }
+               }
+            }
+         }
+         ++counter;
+      }
+      IB[iblock + 1] = counter;
+
+      LUFactors factorization(DB.GetData(iblock), &ipiv[iblock*block_size]);
+      factorization.Factor(block_size);
+   }
+}
+
+void BlockGS::Mult(const Vector &b, Vector &x) const
+{
+   MFEM_VERIFY(height > 0, "BlockGS preconditioner is not constructed");
+   const int nblockrows = Height()/block_size;
+   DenseMatrix L_ij;
+   Vector xi, xj;
+
+   for (int i = 0; i<nblockrows; ++i)
+   {
+      xi.SetDataAndSize(&x[i*block_size], block_size);
+      for (int ib=0; ib<block_size; ++ib)
+      {
+         xi[ib] = b[ib + i*block_size];
+      }
+      for (int k=IB[i]; k<ID[i]; ++k)
+      {
+         const int j = JB[k];
+         xj.SetDataAndSize(&x[j*block_size], block_size);
+         const DenseMatrix &U_ij = AB(k);
+         U_ij.AddMult_a(-1.0, xj, xi);
+      }
+
+      LUFactors A_ii_inv(DB.GetData(i), &ipiv[i*block_size]);
+      A_ii_inv.Solve(block_size, 1, xi.GetData());
+   }
+   x *= damping;
+}
+
+void BlockGS::MultTranspose(const Vector &b, Vector &x) const 
+{
+   MFEM_VERIFY(height > 0, "BlockGS preconditioner is not constructed");
+   const int nblockrows = Height()/block_size;
+   DenseMatrix L_ij;
+   Vector xi, xj;
+
+   for (int i = nblockrows-1; i>=0; --i)
+   {
+      xi.SetDataAndSize(&x[i*block_size], block_size);
+      for (int ib=0; ib<block_size; ++ib)
+      {
+         xi[ib] = b[ib + i*block_size];
+      }
+      for (int k=ID[i] + 1; k<IB[i+1]; ++k)
+      {
+         const int j = JB[k];
+         xj.SetDataAndSize(&x[j*block_size], block_size);
+         const DenseMatrix &L_ij = AB(k);
+         L_ij.AddMult_a(-1.0, xj, xi);
+      }
+      // x_i = D_ii^{-1} x_i
+      LUFactors A_ii_inv(DB.GetData(i), &ipiv[i*block_size]);
+      A_ii_inv.Solve(block_size, 1, xi.GetData());
+   }
+   x *= damping;
+}
+
+
+
+
+BlockILU::BlockILU(int block_size_, real_t damping_,
                    Reordering reordering_,
                    int k_fill_)
    : Solver(0),
      block_size(block_size_),
+     damping(damping_),
      k_fill(k_fill_),
      reordering(reordering_)
 { }
 
 BlockILU::BlockILU(const Operator &op,
-                   int block_size_,
+                   int block_size_, real_t damping_,
                    Reordering reordering_,
                    int k_fill_)
-   : BlockILU(block_size_, reordering_, k_fill_)
+   : BlockILU(block_size_, damping_, reordering_, k_fill_)
 {
    SetOperator(op);
 }
@@ -3191,6 +3500,61 @@ void BlockILU::Mult(const Vector &b, Vector &x) const
       // x_i = D_ii^{-1} x_i
       A_ii_inv.Solve(block_size, 1, xi.GetData());
    }
+   x *= damping;
+}
+
+void BlockILU::MultTranspose(const Vector &b, Vector &x) const
+{
+   // MFEM_VERIFY(height > 0, "BlockILU(0) preconditioner is not constructed");
+   // int nblockrows = Height()/block_size;
+   // y.SetSize(Height());
+
+   // DenseMatrix B;
+   // Vector yi, yj, xi, xj;
+   // Vector tmp(block_size);
+   // y = 0.0;
+   // // Forward substitute to solve (U^T)y = b
+   // for (int i=0; i<nblockrows; ++i)
+   // {
+   //    yi.SetDataAndSize(&y[i*block_size], block_size);
+   //    for (int ib=0; ib<block_size; ++ib)
+   //    {
+   //       yi[ib] = b[ib + P[i]*block_size];
+   //    }
+   //    for (int k=ID[i]+1; k<IB[i+1]; ++k)
+   //    {
+   //       int j = JB[k];
+   //       const DenseMatrix &L_ij = AB(k);
+   //       yj.SetDataAndSize(&y[j*block_size], block_size);
+   //       // y_i = y_i - L_ij*y_j
+   //       L_ij.AddMult_a(-1.0, yj, yi);
+   //    }
+   // }
+   // // Backward substitution to solve (L^T)x = y
+   // for (int i=nblockrows-1; i >= 0; --i)
+   // {
+   //    xi.SetDataAndSize(&x[P[i]*block_size], block_size);
+   //    for (int ib=0; ib<block_size; ++ib)
+   //    {
+   //       xi[ib] = y[ib + i*block_size];
+   //    }
+   //    for (int k=IB[i]; k<ID[i]; ++k)
+   //    {
+   //       int j = JB[k];
+   //       const DenseMatrix &U_ij = AB(k);
+   //       xj.SetDataAndSize(&x[P[j]*block_size], block_size);
+   //       // x_i = x_i - U_ij*x_j
+   //       U_ij.AddMult_a(-1.0, xj, xi);
+   //    }
+   //    LUFactors A_ii_inv(&DB(0,0,i), &ipiv[i*block_size]);
+   //    // x_i = D_ii^{-1} x_i
+   //    A_ii_inv.Solve(block_size, 1, xi.GetData());
+   // }
+
+   // TEMPORARY HACK: assuming this operator is symmetric, Mult and
+   // MultTranspose have the same action.
+   Mult(b, x);
+   x *= damping;
 }
 
 

linalg/solvers.hpp

@@ -1005,6 +1005,94 @@ public:
    void Mult(const Vector &xt, Vector &x) const override;
 };
 
+class BlockGS : public Solver
+{
+public:
+   /** Find L^{-1} for the matrix @a op.
+    *  @a op should be a SparseMatrix.
+    */
+   BlockGS(const Operator &op, int block_size_ = 1, real_t damping_ = 1.0);
+
+   void SetOperator(const Operator &op);
+
+   /// Solve L^{-1} x  = b
+   void Mult(const Vector &b, Vector &x) const;
+
+   /// Solve U^{-1} x = b
+   virtual void MultTranspose(const Vector &b, Vector &x) const override;
+
+
+private:
+   /// @brief Set up the block CSR structure corresponding to a sparse matrix @a A and factorize the diagonal blocks.
+   void CreateBlockPatternAndFactorize(const class SparseMatrix &A);
+
+   const int block_size;
+
+   real_t damping;
+
+   /// Temporary vector used in the Mult() function.
+   mutable Vector y;
+
+   // Block CSR storage
+   Array<int> IB, ID, JB;
+   DenseTensor AB;
+   mutable DenseTensor DB;
+   mutable Array<int> ipiv;
+};
+
+
+class BlockJacobi : public Solver
+{
+public:
+   BlockJacobi(const Operator &op, SparseMatrix &Block_Diag_Mat);
+
+   void SetOperator(const Operator &op);
+
+   void Mult(const Vector &b, Vector &x) const;
+
+   void MultTranspose(const Vector &b, Vector &x) const;
+
+
+private:
+   SparseMatrix Block_Diag_Mat;
+};
+
+class Blockl1Jacobi : public Solver
+{
+public:
+   /** Find L^{-1} for the matrix @a op.
+    *  @a op should be a SparseMatrix.
+    */
+   Blockl1Jacobi(const Operator &op, int block_size_ = 1, real_t damping_ = 1.0);
+
+   void SetOperator(const Operator &op);
+
+   /// Solve L^{-1} x  = b
+   void Mult(const Vector &b, Vector &x) const;
+
+   /// Solve U^{-1} x = b
+   void MultTranspose(const Vector &b, Vector &x) const;
+
+
+private:
+   /// @brief Set up the block CSR structure corresponding to a sparse matrix @a A and factorize the diagonal blocks.
+   void GetDiagonalBlocks(const class SparseMatrix &A);
+
+   const int block_size;
+
+   real_t damping;
+
+   /// Temporary vector used in the Mult() function.
+   mutable Vector y;
+
+   // Block CSR storage
+   Array<int> IB, ID, JB;
+   DenseTensor AB;
+   mutable DenseTensor DB;
+   mutable Array<int> ipiv;
+};
+
+
 /** Block ILU solver:
  *  Performs a block ILU(k) approximate factorization with specified block
  *  size. Currently only k=0 is supported. This is useful as a preconditioner
@@ -1035,7 +1123,7 @@ public:
    /** Create an "empty" BlockILU solver. SetOperator must be called later to
     *  actually form the factorization
     */
-   BlockILU(int block_size_,
+   BlockILU(int block_size_, real_t damping_ = 1.0,
             Reordering reordering_ = Reordering::MINIMUM_DISCARDED_FILL,
             int k_fill_ = 0);
 
@@ -1044,7 +1132,7 @@ public:
     *  case that @a op is a HypreParMatrix, the ILU factorization is performed
     *  on the diagonal blocks of the parallel decomposition.
     */
-   BlockILU(const Operator &op, int block_size_ = 1,
+   BlockILU(const Operator &op, int block_size_ = 1, real_t damping_ = 1.0,
             Reordering reordering_ = Reordering::MINIMUM_DISCARDED_FILL,
             int k_fill_ = 0);
 
@@ -1057,6 +1145,10 @@ public:
    /// Solve the system `LUx = b`, where `L` and `U` are the block ILU factors.
    void Mult(const Vector &b, Vector &x) const;
 
+   /// @brief Solve the system `(LU)^T x = b`, where `L` and `U` are the block
+   /// ILU factors.
+   void MultTranspose(const Vector &b, Vector &x) const;
+
    /** Get the I array for the block CSR representation of the factorization.
     *  Similar to SparseMatrix::GetI(). Mostly used for testing.
     */
@@ -1072,6 +1164,7 @@ public:
     */
    real_t *GetBlockData() { return AB.Data(); }
 
+
 private:
    /// Set up the block CSR structure corresponding to a sparse matrix @a A
    void CreateBlockPattern(const class SparseMatrix &A);
@@ -1086,6 +1179,8 @@ private:
 
    Reordering reordering;
 
+   real_t damping;
+
    /// Temporary vector used in the Mult() function.
    mutable Vector y;
 

miniapps/dg-agglomeration/aux_cg_dg.cpp

@@ -0,0 +1,321 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "mg_agglom.hpp"
+
+using namespace std;
+using namespace mfem;
+
+struct OswaldOperator : Operator
+{
+   const FiniteElementSpace &fes_aux;
+   const FiniteElementSpace &fes;
+   const Array<int> &ess_dofs;
+   Vector multiplicity;
+   mutable Vector z;
+   OswaldOperator(const FiniteElementSpace &fes_aux_,
+                  const FiniteElementSpace &fes_,
+                  const Array<int> &ess_dofs_)
+      : Operator(fes_.GetTrueVSize(), fes_aux_.GetTrueVSize()),
+        fes_aux(fes_aux_),
+        fes(fes_),
+        ess_dofs(ess_dofs_)
+   {
+      const auto ordering = ElementDofOrdering::LEXICOGRAPHIC;
+      const Operator *restr_op = fes.GetElementRestriction(ordering);
+      const auto *restr = dynamic_cast<const ElementRestriction*>(restr_op);
+      MFEM_VERIFY(restr, "");
+
+      multiplicity.SetSize(restr->Width());
+      {
+         Vector ones(restr->Height());
+         ones = 1.0;
+         restr->MultTransposeUnsigned(ones, multiplicity);
+      }
+   }
+
+   SparseMatrix Assemble() const
+   {
+      SparseMatrix R(fes.GetTrueVSize(), fes_aux.GetTrueVSize());
+
+      for (int e = 0; e < fes_aux.GetNE(); ++e)
+      {
+         DenseMatrix I;
+         {
+            const auto T = fes.GetElementTransformation(e);
+            fes.GetFE(e)->Project(*fes_aux.GetFE(e), *T, I);
+         }
+
+         Array<int> vdofs, vdofs_aux;
+         fes_aux.GetElementVDofs(e, vdofs_aux);
+         fes.GetElementVDofs(e, vdofs);
+         R.AddSubMatrix(vdofs, vdofs_aux, I);
+      }
+
+      Vector m = multiplicity;
+      m.Reciprocal();
+      R.ScaleRows(m);
+
+      R.Finalize();
+
+      for (int i : ess_dofs)
+      {
+         R.EliminateRow(i);
+      }
+
+      return R;
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      const int ne = fes_aux.GetNE();
+      y = 0.0;
+
+      for (int e = 0; e < ne; ++e)
+      {
+         DenseMatrix I;
+         {
+            const auto T = fes.GetElementTransformation(e);
+            fes.GetFE(e)->Project(*fes_aux.GetFE(e), *T, I);
+         }
+
+         Array<int> vdofs, vdofs_aux;
+         fes_aux.GetElementVDofs(e, vdofs_aux);
+         fes.GetElementVDofs(e, vdofs);
+         Vector x_e(vdofs_aux.Size()), y_e(vdofs.Size());
+         x.GetSubVector(vdofs_aux, x_e);
+
+         I.Mult(x_e, y_e);
+
+         y.AddElementVector(vdofs, y_e);
+      }
+      for (int i = 0; i < y.Size(); ++i)
+      {
+         y[i] /= multiplicity[i];
+      }
+      for (int i : ess_dofs)
+      {
+         y[i] = 0.0;
+      }
+   }
+
+   void MultTranspose(const Vector &x, Vector &y) const override
+   {
+      const int ne = fes.GetNE();
+      y = 0.0;
+      z = x;
+      for (int i : ess_dofs)
+      {
+         z[i] = 0.0;
+      }
+      for (int i = 0; i < z.Size(); ++i)
+      {
+         z[i] /= multiplicity[i];
+      }
+
+      for (int e = 0; e < ne; ++e)
+      {
+         DenseMatrix I;
+         {
+            const auto T = fes.GetElementTransformation(e);
+            fes.GetFE(e)->Project(*fes_aux.GetFE(e), *T, I);
+         }
+
+         Array<int> vdofs, vdofs_aux;
+         fes.GetElementVDofs(e, vdofs);
+         fes_aux.GetElementVDofs(e, vdofs_aux);
+
+         Vector x_e(vdofs.Size()), y_e(vdofs_aux.Size());
+         z.GetSubVector(vdofs, x_e);
+
+         I.MultTranspose(x_e, y_e);
+
+         y.AddElementVector(vdofs_aux, y_e);
+      }
+   }
+};
+
+struct AuxiliarySolver : Solver
+{
+   const Solver &A_hat_inv;
+   const Operator &R;
+   const Array<int> ess_dofs;
+   const Solver *D;
+   mutable Vector z1, z2, z3;
+   // mutable Vector z;
+
+   AuxiliarySolver(const Solver &A_hat_inv_, const Operator &R_,
+                   const Array<int> &ess_dofs_,
+                   const Solver *D_)
+      : Solver(R_.Height()),
+        A_hat_inv(A_hat_inv_),
+        R(R_),
+        ess_dofs(ess_dofs_),
+        D(D_)
+   { }
+
+   void SetOperator(const Operator &op) { }
+
+   void Mult(const Vector &b, Vector &x) const
+   {
+      z1.SetSize(R.Width());
+      z2.SetSize(R.Width());
+
+      R.MultTranspose(b, z1);
+      A_hat_inv.Mult(z1, z2);
+      R.Mult(z2, x);
+      // A_hat_inv.Mult(b, x);
+
+      if (D)
+      {
+         z3.SetSize(x.Size());
+         D->Mult(b, z3);
+         x += z3;
+      }
+
+      for (int i : ess_dofs)
+      {
+         x[i] = b[i];
+      }
+   }
+};
+
+class CompositeAuxiliaryAgglomerationSolver : public Solver
+{
+   AgglomerationMultigrid a_mg;
+   TruncatedMultigrid t_mg;
+   OswaldOperator oswald;
+   unique_ptr<SparseMatrix> RP;
+   GSSmoother S;
+   AuxiliarySolver aux;
+public:
+   CompositeAuxiliaryAgglomerationSolver(
+      FiniteElementSpace &fes,
+      SparseMatrix &A,
+      FiniteElementSpace &fes_aux,
+      SparseMatrix &A_aux,
+      Array<int> &ess_dofs,
+      int ncoarse,
+      int num_levels,
+      int smoother_choice)
+      : a_mg(fes_aux, A_aux, ncoarse, num_levels, smoother_choice, false),
+        t_mg(a_mg),
+        oswald(fes_aux, fes, ess_dofs),
+        RP(mfem::Mult(oswald.Assemble(), a_mg.GetFinestProlongation())),
+        S(A),
+        aux(t_mg, *RP, ess_dofs, &S)
+   { }
+
+   void SetOperator(const Operator &op) { }
+
+   void Mult(const Vector &b, Vector &x) const
+   {
+      aux.Mult(b, x);
+   }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options.
+   string mesh_file = "../../data/inline-tri.mesh";
+   int order = 1;
+   int ref = 2;
+   real_t kappa_0 = 1.0;
+   int cf = 4;
+   int smoother = 0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree.");
+   args.AddOption(&ref, "-r", "--refine", "Number of refinements.");
+   args.AddOption(&kappa_0, "-k", "--kappa", "Penalty factor.");
+   args.AddOption(&cf, "-cf", "--coarse_factor", "Coarsening Factor.");
+   args.AddOption(&smoother, "-s", "--smoother", "Smoother Choice.");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file, and refine once uniformly.
+   Mesh mesh(mesh_file);
+   for (int i = 0; i < ref; ++i) { mesh.UniformRefinement(); }
+   int dim = mesh.Dimension();
+
+   // 3. Define a finite element space on the mesh. Here we use H1 continuous
+   //    high-order Lagrange finite elements of the given order.
+   H1_FECollection h1_fec(order, mesh.Dimension());
+   FiniteElementSpace h1_fes(&mesh, &h1_fec);
+
+   DG_FECollection dg_fec(order, mesh.Dimension(), BasisType::GaussLobatto);
+   FiniteElementSpace dg_fes(&mesh, &dg_fec);
+
+   cout << "Number of H1 unknowns: " << h1_fes.GetTrueVSize() << endl;
+   cout << "Number of DG unknowns: " << dg_fes.GetTrueVSize() << endl;
+
+   // 4. Extract the list of all the boundary DOFs. These will be marked as
+   //    Dirichlet in order to enforce zero boundary conditions.
+   Array<int> ess_dofs;
+   h1_fes.GetBoundaryTrueDofs(ess_dofs);
+
+   // 5. Define the solution x as a finite element grid function in h1_fes. Set
+   //    the initial guess to zero, which also sets the boundary conditions.
+   GridFunction x(&h1_fes);
+   x = 0.0;
+
+   // 6. Set up the linear form b(.) corresponding to the right-hand side.
+   ConstantCoefficient one(1.0);
+   LinearForm b(&h1_fes);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   BilinearForm a(&h1_fes);
+   a.AddDomainIntegrator(new DiffusionIntegrator);
+   a.AddDomainIntegrator(new MassIntegrator);
+   a.Assemble();
+
+   // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   const real_t sigma = -1.0;
+   const real_t kappa = kappa_0*(order+1)*(order+dim)/dim;
+   BilinearForm a_aux(&dg_fes);
+   a_aux.AddDomainIntegrator(new DiffusionIntegrator);
+   a_aux.AddDomainIntegrator(new MassIntegrator);
+   a_aux.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a_aux.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a_aux.Assemble();
+   a_aux.Finalize();
+
+   SparseMatrix &A_cg= a.SpMat();
+   SparseMatrix &A_dg= a_aux.SpMat();
+   OswaldOperator R_op(dg_fes, h1_fes, ess_dofs);
+   CompositeAuxiliaryAgglomerationSolver prec(h1_fes, A_cg, dg_fes, A_dg, ess_dofs, cf, ref-1, smoother);
+
+   // 8. Form the linear system A X = B. This includes eliminating boundary
+   //    conditions, applying AMR constraints, and other transformations.
+   SparseMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(ess_dofs, x, b, A, X, B);
+
+   {
+      ofstream f("A.txt");
+      A.PrintMatlab(f);
+   }
+   {
+      ofstream f("prec.txt");
+      prec.PrintMatlab(f);
+   }
+
+   CGSolver cg;
+   cg.SetRelTol(1e-7);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   cg.SetPreconditioner(prec);
+   cg.SetOperator(A);
+   cg.Mult(B, X);
+
+   // 10. Recover the solution x as a grid function and save to file. The output
+   //     can be viewed using GLVis as follows: "glvis -m mesh.mesh -g sol.gf"
+   a.RecoverFEMSolution(X, b, x);
+   x.Save("sol.gf");
+   mesh.Save("mesh.mesh");
+
+   return 0;
+}

miniapps/dg-agglomeration/aux_cg_dg_ad_diff.cpp

@@ -0,0 +1,5 @@
+// Understand smoothed aggregation paper and see if we can do something similar.
+
+// no intra aggregate coarsening
+
+// write down summary of paper

miniapps/dg-agglomeration/aux_cg_dg_elasticity.cpp

@@ -0,0 +1,487 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "mg_agglom.hpp"
+
+using namespace std;
+using namespace mfem;
+
+struct OswaldOperator : Operator
+{
+   const FiniteElementSpace &fes_aux;
+   const FiniteElementSpace &fes;
+   const Array<int> &ess_dofs;
+   Vector multiplicity;
+   mutable Vector z;
+   OswaldOperator(const FiniteElementSpace &fes_aux_,
+                  const FiniteElementSpace &fes_,
+                  const Array<int> &ess_dofs_)
+      : Operator(fes_.GetTrueVSize(), fes_aux_.GetTrueVSize()),
+        fes_aux(fes_aux_),
+        fes(fes_),
+        ess_dofs(ess_dofs_)
+   {
+      const auto ordering = ElementDofOrdering::LEXICOGRAPHIC;
+      const Operator *restr_op = fes.GetElementRestriction(ordering);
+      const auto *restr = dynamic_cast<const ElementRestriction*>(restr_op);
+      MFEM_VERIFY(restr, "");
+
+      multiplicity.SetSize(restr->Width());
+      {
+         Vector ones(restr->Height());
+         ones = 1.0;
+         restr->MultTransposeUnsigned(ones, multiplicity);
+      }
+   }
+
+   SparseMatrix Assemble() const
+   {
+      SparseMatrix R(fes.GetTrueVSize(), fes_aux.GetTrueVSize());
+
+      for (int e = 0; e < fes_aux.GetNE(); ++e)
+      {
+         DenseMatrix I;
+         {
+            const auto T = fes.GetElementTransformation(e);
+            fes.GetFE(e)->Project(*fes_aux.GetFE(e), *T, I);
+         }
+
+         Array<int> vdofs, vdofs_aux;
+         fes_aux.GetElementVDofs(e, vdofs_aux);
+         int nodes_per_dim = vdofs_aux.Size()/fes.GetVDim();
+         fes.GetElementVDofs(e, vdofs);
+         for (int vd=0; vd < fes.GetVDim(); vd++)
+         {
+            Array<int> sub_vdofs, sub_vdofs_aux;
+            vdofs.GetSubArray(vd*nodes_per_dim, nodes_per_dim, sub_vdofs);
+            vdofs_aux.GetSubArray(vd*nodes_per_dim, nodes_per_dim, sub_vdofs_aux);
+            R.AddSubMatrix(sub_vdofs, sub_vdofs_aux, I);
+         }
+      }
+
+      Vector m = multiplicity;
+      m.Reciprocal();
+      R.ScaleRows(m);
+
+      R.Finalize();
+
+      for (int i : ess_dofs)
+      {
+         R.EliminateRow(i);
+      }
+      std::cout << "num rows R = " << R.NumRows() << std::endl;
+      std::cout << "num cols R = " << R.NumCols() << std::endl;
+      return R;
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      const int ne = fes_aux.GetNE();
+      y = 0.0;
+
+      for (int e = 0; e < ne; ++e)
+      {
+         DenseMatrix I;
+         {
+            const auto T = fes.GetElementTransformation(e);
+            fes.GetFE(e)->Project(*fes_aux.GetFE(e), *T, I);
+         }
+
+         Array<int> vdofs, vdofs_aux;
+         fes_aux.GetElementVDofs(e, vdofs_aux);
+         fes.GetElementVDofs(e, vdofs);
+         Vector x_e(vdofs_aux.Size()), y_e(vdofs.Size());
+         x.GetSubVector(vdofs_aux, x_e);
+
+         I.Mult(x_e, y_e);
+
+         y.AddElementVector(vdofs, y_e);
+      }
+      for (int i = 0; i < y.Size(); ++i)
+      {
+         y[i] /= multiplicity[i];
+      }
+      for (int i : ess_dofs)
+      {
+         y[i] = 0.0;
+      }
+   }
+
+   void MultTranspose(const Vector &x, Vector &y) const override
+   {
+      const int ne = fes.GetNE();
+      y = 0.0;
+      z = x;
+      for (int i : ess_dofs)
+      {
+         z[i] = 0.0;
+      }
+      for (int i = 0; i < z.Size(); ++i)
+      {
+         z[i] /= multiplicity[i];
+      }
+
+      for (int e = 0; e < ne; ++e)
+      {
+         DenseMatrix I;
+         {
+            const auto T = fes.GetElementTransformation(e);
+            fes.GetFE(e)->Project(*fes_aux.GetFE(e), *T, I);
+         }
+
+         Array<int> vdofs, vdofs_aux;
+         fes.GetElementVDofs(e, vdofs);
+         fes_aux.GetElementVDofs(e, vdofs_aux);
+         Vector x_e(vdofs.Size()), y_e(vdofs_aux.Size());
+         z.GetSubVector(vdofs, x_e);
+
+         I.MultTranspose(x_e, y_e);
+
+         y.AddElementVector(vdofs_aux, y_e);
+      }
+   }
+};
+
+struct AuxiliarySolver : Solver
+{
+   const Solver &A_hat_inv;
+   const Operator &R;
+   const Array<int> ess_dofs;
+   const Solver *D;
+   mutable Vector z1, z2, z3;
+   // mutable Vector z;
+
+   AuxiliarySolver(const Solver &A_hat_inv_, const Operator &R_,
+                   const Array<int> &ess_dofs_,
+                   const Solver *D_)
+      : Solver(R_.Height()),
+        A_hat_inv(A_hat_inv_),
+        R(R_),
+        ess_dofs(ess_dofs_),
+        D(D_)
+   { }
+
+   void SetOperator(const Operator &op) { }
+
+   void Mult(const Vector &b, Vector &x) const
+   {
+      z1.SetSize(R.Width());
+      z2.SetSize(R.Width());
+
+      R.MultTranspose(b, z1);
+      A_hat_inv.Mult(z1, z2);
+      R.Mult(z2, x);
+      // A_hat_inv.Mult(b, x);
+
+      if (D)
+      {
+         z3.SetSize(x.Size());
+         D->Mult(b, z3);
+         x += z3;
+      }
+
+      for (int i : ess_dofs)
+      {
+         x[i] = b[i];
+      }
+   }
+};
+
+class CompositeAuxiliaryAgglomerationSolver : public Solver
+{
+   AgglomerationMultigrid a_mg;
+   TruncatedMultigrid t_mg;
+   OswaldOperator oswald;
+   unique_ptr<SparseMatrix> RP;
+   GSSmoother S;
+   AuxiliarySolver aux;
+public:
+   CompositeAuxiliaryAgglomerationSolver(
+      FiniteElementSpace &fes,
+      SparseMatrix &A,
+      FiniteElementSpace &fes_aux,
+      SparseMatrix &A_aux,
+      Array<int> &ess_dofs,
+      int ncoarse,
+      int num_levels,
+      int smoother_choice)
+      : a_mg(fes_aux, A_aux, ncoarse, num_levels, smoother_choice, false),
+        t_mg(a_mg),
+        oswald(fes_aux, fes, ess_dofs),
+        RP(mfem::Mult(oswald.Assemble(), a_mg.GetFinestProlongation())),
+        S(A),
+        aux(t_mg, *RP, ess_dofs, &S)
+   { }
+
+   void SetOperator(const Operator &op) { }
+
+   void Mult(const Vector &b, Vector &x) const
+   {
+      aux.Mult(b, x);
+   }
+};
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/beam-tri.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool visualization = true;
+   real_t kappa_0 = 10.0;
+   int ref_levels = 4;
+
+   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(&ref_levels, "-r", "--ref_levels",
+                  "Number of times to refine mesh");
+   args.AddOption(&kappa_0, "-k", "--kappa_0",
+                  "DG Penalty Param");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral or hexahedral elements with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
+   {
+      cerr << "\nInput mesh should have at least two materials and "
+           << "two boundary attributes! (See schematic in ex2.cpp)\n"
+           << endl;
+      return 3;
+   }
+
+   // 3. Select the order of the finite element discretization space. For NURBS
+   //    meshes, we increase the order by degree elevation.
+   if (mesh->NURBSext)
+   {
+      mesh->DegreeElevate(order, order);
+   }
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
+   //    largest number that gives a final mesh with no more than 5,000
+   //    elements.
+   // int ref_levels_old = (int)floor(log(5000./mesh->GetNE())/log(2.)/dim);
+   // std::cout << "num ref old = " << ref_levels_old << std::endl;
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a finite element space on the mesh. Here we use vector finite
+   //    elements, i.e. dim copies of a scalar finite element space. The vector
+   //    dimension is specified by the last argument of the FiniteElementSpace
+   //    constructor. For NURBS meshes, we use the (degree elevated) NURBS space
+   //    associated with the mesh nodes.
+   FiniteElementCollection *fec;
+   FiniteElementSpace *fespace;
+   FiniteElementCollection *dg_fec;
+   FiniteElementSpace *dg_fes;
+   if (mesh->NURBSext)
+   {
+      fec = NULL;
+      fespace = mesh->GetNodes()->FESpace();
+   }
+   else
+   {
+      fec = new H1_FECollection(order, dim);
+      fespace = new FiniteElementSpace(mesh, fec, dim);
+      dg_fec = new DG_FECollection(order, dim, BasisType::GaussLobatto);
+      dg_fes = new FiniteElementSpace(mesh, dg_fec, dim);
+   }
+   std::cout << "num elements \n" << mesh->GetNE() << std::endl;
+   cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
+        << endl << "Assembling: " << flush;
+
+   std::cout << "num dg unknowns \n" << dg_fes->GetTrueVSize() << std::endl;
+   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking only
+   //    boundary attribute 1 from the mesh as essential and converting it to a
+   //    list of true dofs.
+   // Array<int> ess_tdof_list;
+   // fespace->GetBoundaryTrueDofs(ess_tdof_list);
+   Array<int> ess_tdof_list, ess_bdr(mesh->bdr_attributes.Max());
+   ess_bdr = 1;
+   // ess_bdr[0] = 1;
+   fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+   // 7. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system. In this case, b_i equals the boundary integral
+   //    of f*phi_i where f represents a "pull down" force on the Neumann part
+   //    of the boundary and phi_i are the basis functions in the finite element
+   //    fespace. The force is defined by the VectorArrayCoefficient object f,
+   //    which is a vector of Coefficient objects. The fact that f is non-zero
+   //    on boundary attribute 2 is indicated by the use of piece-wise constants
+   //    coefficient for its last component.
+   VectorArrayCoefficient f(dim);
+   for (int i = 0; i < dim-1; i++)
+   {
+      f.Set(i, new ConstantCoefficient(1.0));
+   }
+   {
+      Vector pull_force(mesh->bdr_attributes.Max());
+      pull_force = 0.0;
+      pull_force(1) = -1.0e-2;
+      f.Set(dim-1, new PWConstCoefficient(pull_force));
+   }
+
+   Vector g_vec(dim);
+   g_vec = 1.0;
+   VectorConstantCoefficient g(g_vec); 
+   
+
+
+   LinearForm b(fespace);
+   b.AddDomainIntegrator(new VectorDomainLFIntegrator(g));
+   b.AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(f));
+   cout << "r.h.s. ... " << flush;
+   b.Assemble();
+   // b->Finalize();
+
+   // 8. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   GridFunction x(fespace);
+   x = 0.0;
+
+   // 9. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the linear elasticity integrator with piece-wise
+   //    constants coefficient lambda and mu.
+   Vector lambda(mesh->attributes.Max());
+   lambda = 1e5;
+   lambda(0) = lambda(1);
+   PWConstCoefficient lambda_func(lambda);
+   Vector mu(mesh->attributes.Max());
+   mu = 1.0;
+   mu(0) = mu(1);
+   PWConstCoefficient mu_func(mu);
+
+   BilinearForm *a = new BilinearForm(fespace);
+   a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func,mu_func));
+   cout << "matrix ... " << flush;
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+   // a->Finalize();
+
+   // const real_t sigma = -1.0;
+   const real_t kappa = kappa_0*(order+1)*(order+dim)/dim;
+   BilinearForm *a_aux = new BilinearForm(dg_fes);
+   a_aux->AddDomainIntegrator(new ElasticityIntegrator(lambda_func, mu_func));
+   a_aux->AddInteriorFaceIntegrator(
+      new DGElasticityIntegrator(lambda_func, mu_func, -1.0, kappa));
+   a_aux->AddBdrFaceIntegrator(
+      new DGElasticityIntegrator(lambda_func, mu_func, -1.0, kappa), ess_bdr);
+   a_aux->Assemble();
+   a_aux->Finalize();
+
+   SparseMatrix &A_cg= a->SpMat();
+   SparseMatrix &A_dg= a_aux->SpMat();
+   std::cout << "num cols Acg = " << A_cg.NumCols() << std::endl;
+   std::cout << "num cols dg = " << A_dg.NumCols() << std::endl;
+   OswaldOperator R_op(*dg_fes, *fespace, ess_tdof_list);
+   SparseMatrix R_mat = R_op.Assemble();
+   // UMFPackSolver A_dg_inv(A_dg);
+   // AgglomerationMultigrid prec(fespace, A_dg, 4, 2, 0, false);
+   //AuxiliarySolver prec(A_dg_inv, R_mat, ess_tdof_list, nullptr);
+   CompositeAuxiliaryAgglomerationSolver prec(*fespace, A_cg, *dg_fes, A_dg, ess_tdof_list, 4, ref_levels-1, 0);
+
+      // 8. Form the linear system A X = B. This includes eliminating boundary
+   //    conditions, applying AMR constraints, and other transformations.
+   SparseMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   // {
+   //    ofstream f("A.txt");
+   //    A.PrintMatlab(f);
+   // }
+   // {
+   //    ofstream f("prec.txt");
+   //    prec.PrintMatlab(f);
+   // }
+
+
+   CGSolver cg;
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   cg.SetPreconditioner(prec);
+   cg.SetOperator(A);
+   cg.Mult(B, X);
+
+   // 10. Recover the solution x as a grid function and save to file. The output
+   //     can be viewed using GLVis as follows: "glvis -m mesh.mesh -g sol.gf"
+   a->RecoverFEMSolution(X, b, x);
+   
+   // 13. For non-NURBS meshes, make the mesh curved based on the finite element
+   //     space. This means that we define the mesh elements through a fespace
+   //     based transformation of the reference element. This allows us to save
+   //     the displaced mesh as a curved mesh when using high-order finite
+   //     element displacement field. We assume that the initial mesh (read from
+   //     the file) is not higher order curved mesh compared to the chosen FE
+   //     space.
+   if (!mesh->NURBSext)
+   {
+      mesh->SetNodalFESpace(fespace);
+   }
+
+   // 14. Save the displaced mesh and the inverted solution (which gives the
+   //     backward displacements to the original grid). This output can be
+   //     viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
+   {
+      GridFunction *nodes = mesh->GetNodes();
+      *nodes += x;
+      x *= -1;
+      ofstream mesh_ofs("displaced.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+      ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 15. Send the above data by socket to a GLVis server. Use the "n" and "b"
+   //     keys in GLVis to visualize the displacements.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   // 16. Free the used memory.
+   delete a;
+   // delete b;
+   if (fec)
+   {
+      delete fespace;
+      delete fec;
+   }
+   delete mesh;
+
+   return 0;
+}
+
+
+// try with lame parameters = 1 and all dirichlet bcs - behavior should be very similar to poisson 
+// if different - there is a bug
+// agglom only same material?
+// do bigger blocks in block gs

miniapps/dg-agglomeration/dg_agglom.cpp

@@ -0,0 +1,90 @@
+// Copyright (c) 2010-2025, 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.
+//
+//                     -----------------------
+//                     DG Agglomeration Solver
+//                     -----------------------
+
+#include "mfem.hpp"
+#include <iostream>
+#include <memory>
+
+#include "mg_agglom.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   const char *mesh_file = "../../data/square-mixed.mesh";
+   int ref_levels = 2;
+   int order = 1;
+   real_t kappa_0 = 1.0;
+   int ncoarse = 4; 
+   int num_levels = 2;
+   int smoother = 0; // 0 - Block GS, 1 - Block L1 Jacobi, 2 - Block ILU
+   bool paraview_vis = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file.");
+   args.AddOption(&ref_levels, "-r", "--refine", "Refinement levels.");
+   args.AddOption(&order, "-o", "--order", "Polynomial degree.");
+   args.AddOption(&kappa_0, "-k", "--kappa", "DG penalty parameter.");
+   args.AddOption(&ncoarse, "-nc", "--ncoarse", "Number of Fine Elements per Coarse.");
+   args.AddOption(&num_levels, "-nl", "--levels", "Number of Multigrid Levels.");
+   args.AddOption(&smoother, "-s", "--smoother", "Choice of Multigrid Smoother.");
+   args.AddOption(&paraview_vis, "-pv", "--paraview", "-npv", "--no-paraview", "Enable ParaView visualization.");
+   args.ParseCheck();
+
+   Mesh mesh(mesh_file);
+   const int dim = mesh.Dimension();
+
+   for (int i = 0; i < ref_levels; ++i) { mesh.UniformRefinement(); }
+
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fespace(&mesh, &fec);
+   cout << "Number of unknowns: " << fespace.GetVSize() << endl;
+
+   const real_t sigma = -1.0;
+   const real_t kappa = kappa_0 * (order + 1) * (order + 1);
+
+   LinearForm b(&fespace);
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(0.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(zero, one, sigma, kappa));
+   b.Assemble();
+
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a.Assemble();
+   a.Finalize();
+
+   SparseMatrix &A = a.SpMat();
+
+   AgglomerationMultigrid mg(fespace, A, ncoarse, num_levels, smoother, paraview_vis);
+
+   CGSolver cg;
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(500);
+   cg.SetPrintLevel(1);
+   cg.SetOperator(A);
+   cg.SetPreconditioner(mg);
+   cg.Mult(b, x);
+
+   return 0;
+}

miniapps/dg-agglomeration/makefile

@@ -0,0 +1,60 @@
+# Copyright (c) 2010-2025, 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.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+MFEM_INSTALL_DIR ?= ../../mfem
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/dg-agglomeration%/,)
+CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
+   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+DG_HEADERS = partition.hpp mg_agglom.hpp
+DG_SRC = partition.cpp mg_agglom.cpp
+DG_OBJ = $(DG_SRC:.cpp=.o)
+
+MINIAPPS = dg_agglom aux_cg_dg aux_cg_dg_elasticity smooth_agg_gmg_ex
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+.PRECIOUS: %.o
+
+all: $(MINIAPPS)
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+$(MINIAPPS):%: %.o $(DG_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $(DG_OBJ) $< -o $@ $(MFEM_LIBS)
+
+%.o: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf mesh.* sol.* ParaView

miniapps/dg-agglomeration/mg_agglom.cpp

@@ -0,0 +1,771 @@
+
+// Copyright (c) 2010-2025, 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 "mg_agglom.hpp"
+#include "partition.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+TruncatedMultigrid::TruncatedMultigrid(const AgglomerationMultigrid &other)
+{
+   MFEM_VERIFY(other.NumLevels() >= 2, "");
+   const int nlevels = other.NumLevels() - 1;
+
+   operators = other.operators;
+   operators.DeleteLast();
+   ownedOperators.SetSize(nlevels, false);
+
+   smoothers = other.smoothers;
+   smoothers.DeleteLast();
+   ownedSmoothers.SetSize(nlevels, false);
+
+   prolongations = other.prolongations;
+   prolongations.DeleteLast();
+   ownedProlongations.SetSize(nlevels - 1, false);
+}
+
+const SparseMatrix& AgglomerationMultigrid::GetFinestProlongation() const
+{
+   return static_cast<SparseMatrix&>(*prolongations.Last());
+
+}
+
+std::vector<std::vector<std::vector<int>>> Agglomerate(Mesh &mesh, int ncoarse, int num_levels)
+{
+   const int ne = mesh.GetNE();
+
+   const int ne_coarsest = ne / pow(ncoarse, num_levels-1);
+
+   // const int num_partitions = std::ceil(std::log(ne)/std::log(ncoarse));
+
+   std::cout << "number of fine elements: " << ne << std::endl;
+   std::cout << "number of elements per macro elements: " << ncoarse << std::endl;
+   std::cout << "number of coarsest level elements: " << ne_coarsest << std::endl;
+
+   // E is a 3-dimensional data structure which describes how the mesh is partitioned at each level.
+   // The first axis of E refers to each level. So E[i] is the partition information for level i. Note that the 
+   // indexing goes from coarsest to finest. 
+   // Let j be the index of a macro-element on level i. E[i][j] lists the indices of elements on level i+1 belonging to macro-element j.
+   // 'E' data is formed using METIS partitioning. 
+
+   std::vector<std::vector<std::vector<int>>> E(num_levels);
+
+
+   DG_FECollection fec(0, mesh.Dimension());
+   FiniteElementSpace fes(&mesh, &fec);
+
+   // Partition the coarsest mesh.
+   GridFunction p(&fes);
+   p = 0;
+   Array<int> partitioning = PartitionMesh(mesh, ne_coarsest, 4);
+   for (int i = 0; i < p.Size(); ++i)
+   {
+      p[i] = partitioning[i];
+   }
+
+   // Store partitioning into E
+   std::vector<int> coarse_vec;
+   for (int i = 0; i < ne_coarsest; ++i)
+   {
+      coarse_vec.push_back(0);
+   }
+   E[0].push_back(coarse_vec);
+   int total_new_macros = ne_coarsest;
+
+   // Iterate through each level, and populate E. 
+   int j = 1;
+   int ej = ne_coarsest;
+   std::vector<std::vector<int>> macro_el_last;
+   while (j < num_levels)
+   {
+      // If j is >= num_partitions, but we still have not fully refined the mesh, resize E
+      if(j >= num_levels+1){E.resize(E.size()+1);}
+
+      // macro_elements is a data structure which, for each macro element idx i,
+      // lists the indices of all the fine mesh elements which belong to it
+      std::vector<std::vector<int>> macro_elements(total_new_macros);
+      for (int i = 0; i < p.Size(); ++i)
+      {
+         int k = p[i];
+         macro_elements[k].push_back(i);
+      }
+
+      // for each macro_element, partition it
+      total_new_macros = 0;
+      for (int e = 0; e < ej; ++e)
+      {
+         int num_fine_in_macro = macro_elements[e].size(); // number of fine elements in this macro element
+         Array<int> subset(num_fine_in_macro);
+         for (int i=0; i<num_fine_in_macro; i++) {subset[i] = macro_elements[e][i];}
+
+         Array<int> partitioning = PartitionMesh(mesh, ncoarse, 3, subset); // partition this macro element
+
+         // number elements this macro element was partitioned into
+         int new_macros_in_subset = (j+1 != num_levels) ? partitioning.Max() + 1 : macro_elements[e].size(); 
+
+         // fill in data structure E.
+         for (int ip = 0; ip < partitioning.Size(); ++ip)
+         {
+            int i = partitioning[ip];
+            p[subset[ip]] = i + total_new_macros;
+         }
+         std::vector<int> macro_vec;
+         for (int k = 0; k < new_macros_in_subset; ++k) 
+         {
+            macro_vec.push_back(k + total_new_macros);
+         }
+         E[j].push_back(macro_vec);
+
+         total_new_macros += new_macros_in_subset; // record number of elements in next level
+      }
+      j = j+1;
+      ej = total_new_macros;
+
+      // If at the finest level, form macro-elements data structure.
+      if (j == num_levels)
+      {
+         int macro_elements_size = macro_elements.size();
+         for (int ee = 0; ee < macro_elements_size; ee++)
+         {
+            macro_el_last.push_back(macro_elements[ee]);
+         }
+      }
+   }
+
+   // At finest level, fill in E using correct element indices.
+   int num_macro_last_level = E[j-1].size();
+   for(int i = 0; i < num_macro_last_level; i++)
+   {
+      int num_el_in_macro = E[j-1][i].size();
+      for(int k = 0; k < num_el_in_macro; k++)
+      {
+         int bb = macro_el_last[i][k];
+         E[j-1][i][k] = bb;
+      }
+      
+   }
+
+   return E;
+}
+
+
+// void GetMeshSubsetBoundingBox(const mfem::Mesh &mesh, const std::vector<int> &E2l, int c, mfem::Vector &min_coords, mfem::Vector &max_coords) {
+//    int dim = mesh.SpaceDimension();
+//    min_coords.SetSize(dim);
+//    max_coords.SetSize(dim);
+
+//    // Initialize min/max with extreme values
+//    for (int i = 0; i < dim; ++i) {
+//        min_coords(i) = infinity();
+//        max_coords(i) = -infinity();
+//    }
+//    int E2l_size = E2l.size();
+//    for (int i = 0; i < E2l_size; ++i) {
+//        if (E2l[i] == c) {
+//            mfem::Array<int> vert_indices;
+//            mesh.GetElementVertices(i, vert_indices);
+//            for (int j = 0; j < vert_indices.Size(); ++j) {
+//                const double* coords = mesh.GetVertex(vert_indices[j]);
+//                for (int k = 0; k < dim; ++k) {
+//                    if (coords[k] < min_coords(k)) {min_coords(k) = coords[k];};
+//                    if (coords[k] > max_coords(k)) {max_coords(k) = coords[k];};
+//                }
+//            }
+//        }
+//    }
+// }
+
+// SparseMatrix *CreateNodalProlongation(
+//    const std::vector<std::vector<std::vector<int>>> &E, FiniteElementSpace &fes)
+// {
+//    Mesh &mesh = *fes.GetMesh();
+//    int vdim = fes.GetVDim();
+//    int ne = mesh.GetNE();
+//    int dim = mesh.Dimension();
+//    FiniteElementSpace nodal_fes(&mesh, fes.FEColl(), dim);
+//    GridFunction nodes(&nodal_fes);
+//    mesh.GetNodes(nodes);
+//    int nnodes = nodes.Size()/dim*vdim;
+//    int nodes_per_dim = nodes.Size()/dim;
+//    const int n = E.size();
+//    const int p = fes.GetOrder(0);
+//    int d = (dim == 2) ? (p+1)*(p+2)/2 : (p+1)*(p+2)*(p+3)/6;
+//    int num_el_coarse = E[n-1].size();
+//    int ncc = d*E[n-1].size();
+//    int nc = ncc*vdim;
+//    int nr = nnodes;
+//    std::cout << "nodal prolongation rows: " << nr << std::endl;
+//    std::cout << "nodal prolongation columns: " << nc << std::endl;
+//    SparseMatrix *P = new SparseMatrix(nr, nc);
+//    for (int me = 0; me < num_el_coarse; ++me)
+//    {
+
+//       std::vector<int> macro_element = E[n-1][me];
+//       int macro_size = macro_element.size();
+//       for (int el_idx = 0; el_idx < macro_size; el_idx++)
+//       {
+//          int el = macro_element[el_idx];
+//          Vector bb_min_coarse; Vector bb_max_coarse;
+//          GetMeshSubsetBoundingBox(mesh, macro_element, el_idx, bb_min_coarse, bb_max_coarse);
+//          Array<int> local_element_dof_indices;
+//          fes.GetElementDofs(el, local_element_dof_indices);
+//          int num_el_dofs = local_element_dof_indices.Size();
+//          for (int i=0; i < num_el_dofs; ++i)
+//          {
+//             int dof_idx = local_element_dof_indices[i];
+//             double x_phys = nodes(dof_idx); double y_phys = nodes(nodes_per_dim + dof_idx);
+//             double x_ref = (x_phys - bb_min_coarse(0)) / (bb_max_coarse(0) - bb_min_coarse(0));
+//             double y_ref = (y_phys - bb_min_coarse(1)) / (bb_max_coarse(1) - bb_min_coarse(1));
+//             IntegrationPoint ip;
+//             Vector shape_vec(d);
+//             if (dim == 3){
+//                L2_TetrahedronElement rfe(p);
+//                double z_phys = nodes(2*nodes_per_dim + dof_idx);
+//                double z_ref = (z_phys - bb_min_coarse(2)) / (bb_max_coarse(2) - bb_min_coarse(2));
+//                ip.Set(x_ref, y_ref, z_ref, 1);
+//                rfe.CalcShape(ip, shape_vec);
+//             }
+//             else
+//             {
+//                L2_TriangleElement rfe(p);
+//                ip.Set2w(x_ref, y_ref, 1);
+//                rfe.CalcShape(ip, shape_vec);
+//             }
+//             for (int k = 0; k < d; k++)
+//             {
+//                for(int vd=0; vd < vdim; vd++)
+//                {
+//                   P -> Set(fes.DofToVDof(dof_idx, vd), (ncc*vd) + d*me + k, shape_vec(k));
+//                }
+//             }
+//          }  
+//       }
+//    }
+//    P -> Finalize();
+//    return P;
+// }
+
+// SparseMatrix *CreateInclusionProlongation(
+//    int l, const std::vector<std::vector<int>> &E, const std::vector<std::vector<int>> &E2, FiniteElementSpace &fes)
+// {
+//    Mesh &mesh = *fes.GetMesh();
+//    int vdim = fes.GetVDim();
+//    int dim = mesh.Dimension();
+//    FiniteElementSpace nodal_fes(&mesh, fes.FEColl(), dim);
+//    const int p = fes.GetOrder(0);
+//    int d = (dim == 2) ? (p+1)*(p+2)/2 : (p+1)*(p+2)*(p+3)/6;
+//    GridFunction nodes(&nodal_fes);
+//    int ncc = (l == 0) ? d: d*E[l-1].size();
+//    int nc = vdim*ncc;
+//    int nrr = d*E[l].size();
+//    int nr = vdim*nrr;
+//    std::cout << "inclusion prolongation rows: " << nr << std::endl;
+//    std::cout << "inclusion prolongation columns: " << nc << std::endl;
+//    SparseMatrix *P = new SparseMatrix(nr, nc);
+//    int el_size = E[l].size();
+//    int num_el_coarse = E[l-1].size();
+//    for (int me = 0; me < num_el_coarse; ++me)
+//     {
+//       std::vector<int> macro_element = E[l][me];
+//       int macro_size = macro_element.size();
+//       Vector bb_min_coarse; Vector bb_max_coarse;
+//       GetMeshSubsetBoundingBox(mesh, macro_element, me, bb_min_coarse, bb_max_coarse);
+//       for (int el_idx = 0; el_idx < macro_size; el_idx++)
+//       { 
+
+//       }
+//       Vector bb_min_fine; Vector bb_max_fine;
+//       GetMeshSubsetBoundingBox(mesh, E2[l+1], e, bb_min_fine, bb_max_fine);
+//       if (dim == 2)
+//       {
+//          L2_TriangleElement rfe(p);
+//          const IntegrationRule rfe_nodes = rfe.GetNodes();
+//          for (int i = 0; i < rfe_nodes.Size(); i++)
+//          {
+//             IntegrationPoint ip = rfe_nodes.IntPoint(i);
+//             Vector small_bb_map(2);
+//             small_bb_map(0) = (bb_max_fine(0) - bb_min_fine(0))*ip.x + bb_min_fine(0);
+//             small_bb_map(1) = (bb_max_fine(1) - bb_min_fine(1))*ip.y + bb_min_fine(1);
+//             Vector ref_coord_big(2); 
+//             ref_coord_big(0) = (small_bb_map(0) - bb_min_coarse(0))/(bb_max_coarse(0) - bb_min_coarse(0));
+//             ref_coord_big(1) = (small_bb_map(1) - bb_min_coarse(1))/(bb_max_coarse(1) - bb_min_coarse(1));
+//             IntegrationPoint ip2;
+//             ip2.Set2w(ref_coord_big(0), ref_coord_big(1), 1);
+//             Vector shape_vec(4);
+//             rfe.CalcShape(ip2, shape_vec);
+//             for (int k = 0; k < d; k++)
+//             {
+//                for (int vd = 0; vd < vdim; vd++)
+//                {
+//                   P -> Set((nrr*vd) + d*e + i, (ncc*vd) + d*c + k, shape_vec(k));
+//                }
+//             } 
+//          }
+//       }
+//       else
+//       {
+//          L2_TetrahedronElement rfe(p);
+//          const IntegrationRule rfe_nodes = rfe.GetNodes();
+//          for (int i = 0; i < rfe_nodes.Size(); i++)
+//          {
+//             IntegrationPoint ip = rfe_nodes.IntPoint(i);
+//             Vector small_bb_map(3);
+//             small_bb_map(0) = (bb_max_fine(0) - bb_min_fine(0))*ip.x + bb_min_fine(0);
+//             small_bb_map(1) = (bb_max_fine(1) - bb_min_fine(1))*ip.y + bb_min_fine(1);
+//             small_bb_map(2) = (bb_max_fine(2) - bb_min_fine(2))*ip.z + bb_min_fine(2);
+//             Vector ref_coord_big(3); 
+//             ref_coord_big(0) = (small_bb_map(0) - bb_min_coarse(0))/(bb_max_coarse(0) - bb_min_coarse(0));
+//             ref_coord_big(1) = (small_bb_map(1) - bb_min_coarse(1))/(bb_max_coarse(1) - bb_min_coarse(1));
+//             ref_coord_big(2) = (small_bb_map(2) - bb_min_coarse(2))/(bb_max_coarse(2) - bb_min_coarse(2));
+//             IntegrationPoint ip2;
+//             ip2.Set3(ref_coord_big(0), ref_coord_big(1), ref_coord_big(2));
+//             Vector shape_vec(8);
+//             rfe.CalcShape(ip2, shape_vec);
+//             for (int k = 0; k < d; k++)
+//             {
+//                for (int vd = 0; vd < vdim; vd++)
+//                {
+//                   P -> Set((nrr*vd) + d*e + i, (ncc*vd) + d*c + k, shape_vec(k));
+//                }
+//             } 
+//          }
+//       }
+//    }
+//    P->Finalize();
+//    return P;
+// }
+
+// AgglomerationMultigrid::AgglomerationMultigrid(
+//    FiniteElementSpace &fes, SparseMatrix &Af, int ncoarse, int num_levels, int smoother_choice, bool paraview_vis)
+// {
+//    Mesh &mesh = *fes.GetMesh();
+//    const int ne = mesh.GetNE();
+
+//    // Create the mesh hierarchy
+//    // E2 is a data structure such that E2_ij = gives the level i index for fine element j
+//    auto E = Agglomerate(*fes.GetMesh(), ncoarse, num_levels);
+//    // vector<vector<int>> E2(E.size());
+//    // E2.back() = E.back();
+//    // for (int i = E.size() - 2; i >= 0; --i)
+//    // {
+//    //    E2[i].resize(ne);
+//    //    for (int e = 0; e < ne; ++e)
+//    //    {
+//    //       const int m_e = E2[i+1][e];
+//    //       E2[i][e] = E[i][m_e];
+//    //    }
+//    // }
+//    // vector<int> E2_finest_level(ne);
+//    // for (int i = 0; i < ne; ++i) 
+//    // {
+//    //    E2_finest_level[i] = i; 
+//    // }
+//    // E2.push_back(E2_finest_level);
+
+//    // // output a paraview visualization of the partition, if desired
+//    // if (paraview_vis)
+//    // {
+//    //    L2_FECollection l2_fec(0, mesh.Dimension());
+//    //    FiniteElementSpace l2_fes(&mesh, &l2_fec);
+//    //    GridFunction p_gf(&l2_fes);
+//    //    ParaViewDataCollection pv("Agglomeration", &mesh);
+//    //    pv.SetPrefixPath("ParaView");
+//    //    pv.RegisterField("p", &p_gf);
+//    //    int E_size = E.size();
+//    //    for (int i = 0; i < E_size; ++i)
+//    //    {
+//    //       for (int e = 0; e < ne; ++e)
+//    //       {
+//    //          p_gf[e] = E2[i][e];
+//    //       }
+//    //       pv.SetCycle(i);
+//    //       pv.SetTime(i);
+//    //       pv.Save();
+//    //    }
+//    //    for (int e = 0; e < ne; ++e)
+//    //    {
+//    //       p_gf[e] = e;
+//    //       pv.SetCycle(E.size());
+//    //       pv.SetTime(E.size());
+//    //       pv.Save();
+//    //    }
+//    // }
+
+//    std::cout << "num levels: " << num_levels << std::endl;
+
+//    // Populate the arrays: operators, smoothers, ownedOperators, ownedSmoothers
+//    // from the MultigridBase class. (All smoothers are owned, all operators
+//    // except the finest are owned).
+//    operators.SetSize(num_levels);
+//    smoothers.SetSize(num_levels);
+//    ownedOperators.SetSize(num_levels);
+//    ownedSmoothers.SetSize(num_levels);
+//    prolongations.SetSize(num_levels-1);
+//    ownedProlongations.SetSize(num_levels-1);
+
+//    //Set the ownership
+//    for (int l = 0; l < num_levels-1; ++l)
+//    {
+//       ownedOperators[l] = true;
+//       ownedSmoothers[l] = true;
+//       ownedProlongations[l] = true;
+//    }
+//    ownedOperators[num_levels-1] = false;
+//    ownedSmoothers[num_levels-1] = true;
+
+//    // Populate the arrays: prolongations, ownedProlongations from the Multigrid
+//    // class. All prolongations are owned.
+//    // Create the prolongations using 'E' using the SparseMatrix class
+//    operators[num_levels - 1] = &Af;
+//    int k = num_levels;
+//    for (int l = num_levels - 2; l >= 0; --l)
+//    {
+//       SparseMatrix *P;
+//       if (l < num_levels - 2)
+//       {
+//          // P = CreateInclusionProlongation(k, E, E2, fes);
+//          SparseMatrix &A_prev = static_cast<SparseMatrix&>(*operators[l + 1]);
+//          // unique_ptr<SparseMatrix> AP(mfem::Mult(A_prev, *P));
+//          // unique_ptr<SparseMatrix> Pt(Transpose(*P));
+//          // operators[l] = mfem::Mult(*Pt, *AP);
+//          // prolongations[l] = P;
+//       }
+//       else
+//       {
+//          P = CreateNodalProlongation(E, fes);
+//          SparseMatrix &A_prev = static_cast<SparseMatrix&>(*operators[l + 1]);
+//          std::cout << "num cols Af = " << A_prev.NumCols() << std::endl;
+//          unique_ptr<SparseMatrix> AP(mfem::Mult(A_prev, *P));
+//          unique_ptr<SparseMatrix> Pt(Transpose(*P));
+//          operators[l] = mfem::Mult(*Pt, *AP);
+//          prolongations[l] = P;
+//       }
+//       k = k-1;
+//    }
+
+//    // Create the smoothers remember block size is num degrees of freedom per element
+//    SparseMatrix &Ac = static_cast<SparseMatrix&>(*operators[0]);
+//    smoothers[0] = new UMFPackSolver(Ac);
+//    int block_size = 3;
+//    real_t damping = 1.0;
+//    for (int l=1; l < num_levels; l++)
+//    {
+//       if (smoother_choice == 0)
+//       {
+//          smoothers[l] = new BlockGS(*operators[l], block_size, damping);
+//       }
+//       else if (smoother_choice == 1)
+//       {
+//          smoothers[l] = new Blockl1Jacobi(*operators[l], block_size, damping);
+//       }
+//       else if (smoother_choice == 2)
+//       {
+//          smoothers[l] = new BlockILU(*operators[l], block_size, damping);
+//       }
+//       else
+//       {
+//          MFEM_ABORT("Unknown Smoother.")
+//       }
+//    }
+// }
+
+/** Construct Block Jacobi Smoother
+ *  Extracts the block diagonal of A, given a block size. Then inverts the matrix by inverting
+ *  each of the blocks.
+ */
+SparseMatrix ExtractBlockDiagonalInverse(SparseMatrix &A, Array<int> block_sizes)
+{
+   SparseMatrix Block_Diag_Mat(A.Height(), A.Width());
+   const int nblockrows = block_sizes.Size();
+   const int *I = A.HostReadI();
+   const int *J = A.HostReadJ();
+   const real_t *V = A.HostReadData();
+   int i = 0;
+   for (int iblock = 0; iblock < nblockrows; ++iblock)
+   {
+      int blocksize = block_sizes[iblock];
+      DenseMatrix Block(blocksize, blocksize);
+      Block = 0.0;
+      int start_j = i;
+      // Extract the block from the matrix A
+      for (int bi = 0; bi < blocksize; ++bi)
+      {
+         for (int k = I[i]; k < I[i+1]; ++k)
+         {
+            int j = J[k];
+            real_t val = V[k];
+            if (j >= start_j && j < start_j + blocksize)
+            {
+               int bj = j - start_j;
+               Block(bi, bj) = val;
+            }
+         }
+         i = i+1;
+      }
+      // Invert the block
+      Block.Invert();
+      // Insert the inverted block into the block diagonal matrix
+      for (int r = 0; r < blocksize; r++)
+      {
+         for(int c = 0; c < blocksize; c++)
+         {
+            int row_idx = start_j + r;
+            int col_idx = start_j + c;
+            Block_Diag_Mat.Set(row_idx, col_idx, Block(r, c));
+         
+         }
+      }
+   }
+   Block_Diag_Mat.Finalize();
+   return Block_Diag_Mat;
+}
+
+/** Adaptive Smoother
+ *  Returns the block Jacobi Smoother with appropriate damping coefficient. In addition, also smooths
+ *  the columns of the matrix B.
+ */
+SparseMatrix AdaptiveSmoother(Operator &op, SparseMatrix &A, Array<int> block_sizes, Vector &x_random, DenseMatrix &B, int level)
+{
+   // Find optimal damping coefficient through power iteration.
+   Vector PinvAx = x_random;
+   SparseMatrix smoother = ExtractBlockDiagonalInverse(A, block_sizes);
+   for(int i = 0; i < 4; i++)
+   {
+      Vector Ax(PinvAx.Size());
+      A.Mult(PinvAx, Ax);
+      smoother.Mult(Ax, PinvAx);
+      if (i != 3){PinvAx /= PinvAx.Norml2();}
+   }
+   real_t rho = PinvAx.Norml2();
+   real_t w = 4.0/3.0/rho;
+   smoother *= w;
+
+   // Smooth the columns of B.
+   // real_t tol = 1.03;
+   int num_B_cols = B.Width();
+   int num_B_rows = B.Height();
+   for(int j = 0; j < num_B_cols; j++)
+   {
+      Vector b(num_B_rows);
+      Vector b_prev(num_B_rows);
+      Vector Ab(num_B_rows);
+      B.GetColumn(j, b);
+
+      // real_t norm_prev = b.Norml2();
+      // real_t ratio;
+      int it = 0;
+      do
+      {
+         b_prev = b;
+         A.Mult(b, Ab);
+         smoother.Mult(Ab, b);
+         b *= -1.0;
+         b += b_prev;
+         // real_t norm = b.Norml2();
+         // ratio = norm_prev / norm;
+         // norm_prev = norm;
+         it += 1;
+      }
+      while(it < 80);
+
+      for(int i = 0; i < num_B_rows; i++){B(i, j) = b(i);}
+   }
+   return smoother;
+}
+
+SmoothedAggregationGMG::SmoothedAggregationGMG(FiniteElementSpace &fes, SparseMatrix &Af, int ncoarse, int num_levels, bool paraview_vis)
+{
+   Mesh &mesh = *fes.GetMesh();
+   int vdim = fes.GetVDim();
+   int ne = mesh.GetNE();
+   int dim = mesh.Dimension();
+   int n_cut = (dim == 2) ? 2*2 - 2 + 1: 2*2*2 - 3 + 1;
+   FiniteElementSpace nodal_fes(&mesh, fes.FEColl(), dim);
+   GridFunction nodes(&nodal_fes);
+   mesh.GetNodes(nodes);
+   int nnodes = nodes.Size()/dim*vdim;
+
+   // Create the mesh hierarchy.
+   auto E = Agglomerate(*fes.GetMesh(), ncoarse, num_levels);
+   std::cout << "num agglomerated levels " << E.size() <<  std::endl;
+
+   // Initialize B
+   int num_samp = ncoarse*ncoarse; // 16 if 2-dimensional, 64 if 3-dimensional.
+   DenseMatrix B(nnodes, num_samp); 
+   std::random_device rd;
+   std::mt19937 gen(rd()); 
+   std::normal_distribution<double> dist(0.0, 1.0);
+   for (int i = 0; i < B.Height(); i++)
+   {
+      for (int j = 0; j < B.Width(); j++)
+      {
+         B(i,j) = dist(gen);
+      }
+   }
+
+   // Populate the arrays: operators, smoothers, ownedOperators, ownedSmoothers
+   // from the MultigridBase class. (All smoothers are owned, all operators
+   // except the finest are owned).
+   operators.SetSize(num_levels);
+   smoothers.SetSize(num_levels);
+   ownedOperators.SetSize(num_levels);
+   ownedSmoothers.SetSize(num_levels);
+   prolongations.SetSize(num_levels-1);
+   ownedProlongations.SetSize(num_levels-1);
+   //Set the ownership
+   for (int l = 0; l < num_levels-1; ++l)
+   {
+      ownedOperators[l] = true;
+      ownedSmoothers[l] = true;
+      ownedProlongations[l] = true;
+   }
+   ownedOperators[num_levels-1] = false;
+   ownedSmoothers[num_levels-1] = true;
+   operators[num_levels - 1] = &Af;
+
+   // prev_dof_idx is an array initialized to be of length num_elements in fine mesh + 1.
+   // prev_dof_idx[i]:prev_dof_idx[i+1]-1 is range containing indices for dofs in element i.
+   // curr_dof_idx (see below) is similar in structure, but prev_dof_idx is for the "previous" level
+   // and curr_dof_idx is for the "current" level.
+   Array<int> prev_dof_idx(ne+1);
+   prev_dof_idx = 0;
+   Array<int> block_sizes(ne);
+   for (int i = 0; i<ne; i++)
+   {
+      const FiniteElement *fe = fes.GetFE(i);
+      int num_nodes = fe->GetDof();
+      prev_dof_idx[i+1] = prev_dof_idx[i]+num_nodes;
+      block_sizes[i] = num_nodes;
+   }
+
+
+   for (int k = num_levels - 2; k >= 0; --k)
+   {
+      SparseMatrix &A_prev = static_cast<SparseMatrix&>(*operators[k + 1]);
+
+      // Generate a random vector, x, to initialize power iteration.
+      Vector x_random(A_prev.Height());
+      for (int i = 0; i < x_random.Size(); i++){x_random(i) = dist(gen);}
+      // x_random /= x_random.Norml2();
+
+
+      // Find the damping coefficient. Smooth columns of B.
+      // int block_size_in = (k == num_levels-2) ? ncoarse : n_cut;
+      SparseMatrix A_tilde_inv = AdaptiveSmoother(*operators[k + 1], A_prev, block_sizes, x_random, B, k);
+
+      int num_el_part = E[k+1].size()+1;
+
+      Array<int> curr_dof_idx(num_el_part);
+      curr_dof_idx = 0;
+
+      // Parr is a vector, members of which are size-3 mfem Vectors. These mfem vectors contain
+      // the column idx, the row idx, and value for cells in the matrix P.
+      std::vector<Vector> Parr;
+      block_sizes.SetSize(num_el_part-1);
+      std::unique_ptr<SparseMatrix> P(new SparseMatrix(prev_dof_idx.Last()));
+      // Loop over elements in current level
+      for (int j = 0; j < num_el_part - 1; j++)
+      {
+         // Get relevant dofs for current coarse element.
+         std::vector<int> indices;
+         int num_el_in_macro = E[k+1][j].size();
+         for (int i = 0; i < num_el_in_macro; i++)
+         {
+            int kk = E[k+1][j][i];
+            for (int mm=prev_dof_idx[kk]; mm<prev_dof_idx[kk+1]; mm++)
+            {
+               indices.push_back(mm);
+            }
+         }
+
+         // Get relevant submatrix of B
+         DenseMatrix B_sub(indices.size(), B.Width());
+         int indices_size = indices.size();
+         for (int i = 0; i < indices_size; i++)
+         {
+            for (int c = 0; c <  B.Width(); c++)
+            {
+               B_sub(i, c) = B(indices[i], c);
+            }
+         }
+
+         // Perform SVD on the B submatrix. Grab the left singular vectors.
+         DenseMatrixSVD Bsvd(B_sub, 'A', 'N'); 
+         Bsvd.Eval(B_sub);
+         DenseMatrix U = Bsvd.LeftSingularvectors(); 
+         Vector S = Bsvd.Singularvalues();
+
+         // Fill in Parr. 
+         int num_keep = indices.size()/n_cut + (indices.size() % n_cut != 0);
+
+         block_sizes[j] = num_keep;
+         for (int r = 0; r < U.Height(); r++)
+         {
+            for (int c = 0; c < num_keep; c++)
+            {
+               real_t val = U(r,c);
+               P -> Set(indices[r], curr_dof_idx[j] + c, val);
+            }
+         }
+
+         // update curr_dof_idx
+         curr_dof_idx[j+1] = curr_dof_idx[j] + num_keep; 
+      }
+      P->OverrideSize(prev_dof_idx.Last(), curr_dof_idx.Last());
+      P->Finalize();
+
+      // Construct the prolongation matrix, T, by computing T= (I - \tilde{A}^{-1} A)P
+      std::unique_ptr<SparseMatrix> A_til_inv_A_mat_ptr(mfem::Mult(A_tilde_inv, A_prev));
+      SparseMatrix A_til_inv_A_mat = *A_til_inv_A_mat_ptr; // \tilde{A}^{-1} A
+      Vector ones(A_prev.Height());
+      ones = 1.0;
+      SparseMatrix Id(ones);
+      A_til_inv_A_mat *= -1; // -\tilde{A}^{-1} A
+      A_til_inv_A_mat.Add(1.0, Id); // I - \tilde{A}^{-1} A
+      SparseMatrix *T = mfem::Mult(A_til_inv_A_mat, *P); // T = (I - \tilde{A}^{-1} A)P
+      T->Threshold(1e-6, false);
+      prolongations[k] = T;
+
+      // Get the coarse operator by perfoming R A T, with R = T^T
+      unique_ptr<SparseMatrix> AP(mfem::Mult(A_prev, *T));
+      unique_ptr<SparseMatrix> Pt(Transpose(*T));
+      SparseMatrix *Anew = mfem::Mult(*Pt, *AP);
+      Anew->Threshold(1e-6, false);
+      operators[k] = Anew;
+
+      // Set the smoother. Note: I would get an error when assigning the SparseMatrix A_tilde_inv
+      // to smoothers[k+1]. So, I created a "dummy" solver object BlockJacobi which takes in A_tilde_inv and essentially
+      // performs the action of multiplying with A_tilde_inv.
+      smoothers[k+1] = new BlockJacobi(*operators[k+1], A_tilde_inv);
+      DenseMatrix *B_new = mfem::Mult(*Pt, B); // Initialize B for next level 
+      B = *B_new;
+      delete B_new;
+      prev_dof_idx = curr_dof_idx; // set prev_dof_idx for next level 
+   }
+   SparseMatrix &Ac = static_cast<SparseMatrix&>(*operators[0]);
+   smoothers[0] = new UMFPackSolver(Ac);
+
+   for (int i = num_levels - 1; i >= 0; --i)
+   {
+      std::cout << "Level " << i << ": " << operators[i]->Height() << " DOFs. ";
+      std::cout << "nnz = " << static_cast<SparseMatrix*>
+                (operators[i])->NumNonZeroElems() << ".\n";
+   }
+}
+} // namespace mfem
+
+
+

miniapps/dg-agglomeration/mg_agglom.hpp

@@ -0,0 +1,44 @@
+// Copyright (c) 2010-2025, 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_MG_AGGLOM_HPP
+#define MFEM_MG_AGGLOM_HPP
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+std::vector<std::vector<int>> Agglomerate(Mesh &mesh);
+
+class AgglomerationMultigrid : public Multigrid
+{
+friend class TruncatedMultigrid;
+public:
+   AgglomerationMultigrid(FiniteElementSpace &fes, SparseMatrix &Af, int ncoarse, int num_levels, int smoother_choice, bool paraview_vis);
+
+   const SparseMatrix& GetFinestProlongation() const;
+};
+
+class TruncatedMultigrid : public Multigrid
+{
+public:
+   TruncatedMultigrid(const AgglomerationMultigrid &other);
+};
+
+class SmoothedAggregationGMG : public Multigrid
+{
+   public:
+   SmoothedAggregationGMG(FiniteElementSpace &fes, SparseMatrix &Af, int ncoarse, int num_levels, bool paraview_vis);
+};
+} // namespace mfem
+
+#endif

miniapps/dg-agglomeration/optref.hpp

@@ -0,0 +1,39 @@
+// Copyright (c) 2010-2025, 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_OPTREF_HPP
+#define MFEM_OPTREF_HPP
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+class NullOptT { };
+
+template <typename T>
+class OptRef
+{
+   T *ptr = nullptr;
+public:
+   OptRef() = default;
+   OptRef(NullOptT) { }
+   OptRef(T &t) : ptr(&t) { }
+   operator bool() const { return ptr; }
+   T &operator*() const { return *ptr; }
+   T *operator->() const { return ptr; }
+};
+
+static constexpr NullOptT NullOpt;
+
+} // namespace mfem
+
+#endif

miniapps/dg-agglomeration/partition.cpp

@@ -0,0 +1,201 @@
+#include "partition.hpp"
+
+#ifdef MFEM_USE_METIS_5
+#include "metis.h"
+#else
+#error "METIS is required"
+#endif
+
+namespace mfem
+{
+
+Array<idx_t> PartitionMesh(Mesh &mesh, const int npart, const int part_method,
+                           OptRef<Array<int>> subset)
+{
+   // const int part_method = 3;
+
+   const int ne = subset ? subset->Size() : mesh.GetNE();
+
+   // if (subset){
+   //    for (int l = 0; l < 5; l++){
+   //       std::cout << "sub[l]: " << (*subset)[l] << std::endl;
+   //    }
+   // }
+
+   idx_t mpart = npart;
+   Array<idx_t> p(ne);
+
+   // Early return special cases for one partition requested, or more partitions
+   // then elements.
+   if (npart == 1)
+   {
+      for (int i = 0; i < ne; i++)
+      {
+         p[i] = 0;
+      }
+      return p;
+   }
+   else if (ne <= npart)
+   {
+      for (int i = 0; i < ne; i++)
+      {
+         p[i] = i;
+      }
+      return p;
+   }
+
+   const Table &e2e = mesh.ElementToElementTable();
+   Array<idx_t> I, J;
+   {
+      const int *iI = e2e.HostReadI();
+      const int *iJ = e2e.HostReadJ();
+      const int m = iI[ne];
+      I.SetSize(ne+1);
+      if (subset){
+         I[0] = 0;
+         for (int r=0; r<ne; r++){
+            I[r+1] = I[r];
+            for (int c=0; c<ne; c++){
+               bool non_zero = false; 
+               for (int i=iI[(*subset)[r]]; i < iI[(*subset)[r]+1]; i++){
+                  non_zero = (iJ[i] == (*subset)[c] ? true : non_zero);
+               }
+               if (non_zero){I[r+1] += 1; J.Append(c);}
+            }
+         }
+      //for (int k = 0; k < ne+1; k++){std::cout << "k: " << k << std::endl; std::cout << "I[k]: " << I[k] << std::endl;}
+      //for (int k = 0; k < J.Size(); k++){std::cout << "k: " << k << std::endl; std::cout << "J[k]: " << J[k] << std::endl;}
+      }
+      else{
+         J.SetSize(m);
+         for (int k = 0; k < ne + 1; k++) { I[k] = iI[k];}
+         for (int k = 0; k < m; k++) { J[k] = iJ[k];}
+      }
+   }
+
+   idx_t options[40];
+   METIS_SetDefaultOptions(options);
+   options[METIS_OPTION_CONTIG] = 1; // set METIS_OPTION_CONTIG
+
+   // If the mesh is disconnected, disable METIS_OPTION_CONTIG.
+   // {
+   //    Array<int> part(partitioning, ne);
+   //    part = 0; // single part for the whole mesh
+   //    Array<int> component; // size will be set to num. elem.
+   //    Array<int> num_comp;  // size will be set to num. parts (1)
+   //    mesh.FindPartitioningComponents(*el_to_el, part, component, num_comp);
+   //    if (num_comp[0] > 1) { options[METIS_OPTION_CONTIG] = 0; }
+   // }
+
+   // Sort the neighbor lists
+   if (part_method >= 0 && part_method <= 2)
+   {
+      for (int i = 0; i < ne; i++)
+      {
+         // Sort in increasing order.
+         // std::sort(J+I[i], J+I[i+1]);
+
+         // Sort in decreasing order, as in previous versions of MFEM.
+         std::sort(J+I[i], J+I[i+1], std::greater<idx_t>());
+      }
+   }
+   // This function should be used to partition a graph into a small
+   // number of partitions (less than 8).
+   if (part_method == 0 || part_method == 3)
+   {
+      idx_t n = ne;
+      idx_t ncon = 1;
+      idx_t edgecut;
+      const idx_t err = METIS_PartGraphRecursive(
+                           &n, &ncon, I, J, NULL, NULL, NULL, &mpart, NULL,
+                           NULL, options, &edgecut, p.HostWrite());
+      MFEM_VERIFY(err == 1, "Error in METIS_PartGraphRecursive");
+   }
+
+   // This function should be used to partition a graph into a large
+   // number of partitions (greater than 8).
+   if (part_method == 1 || part_method == 4)
+   {
+      idx_t n = ne;
+      idx_t ncon = 1;
+      idx_t edgecut;
+      const idx_t err = METIS_PartGraphKway(
+                           &n, &ncon, I, J, NULL, NULL, NULL, &mpart, NULL,
+                           NULL, options, &edgecut, p.HostWrite());
+      MFEM_VERIFY(err == 1, "Error in METIS_PartGraphKway");
+   }
+
+   // Check for empty partitionings (a "feature" in METIS)
+   // if (npart > 1 && ne > npart)
+   // {
+   //    Array< Pair<int,int> > psize(npart);
+   //    int empty_parts;
+
+   //    // Count how many elements are in each partition, and store the result in
+   //    // psize, where psize[i].one is the number of elements, and psize[i].two
+   //    // is partition index. Keep track of the number of empty parts.
+   //    auto count_partition_elements = [&]()
+   //    {
+   //       for (int i = 0; i < npart; i++)
+   //       {
+   //          psize[i].one = 0;
+   //          psize[i].two = i;
+   //       }
+
+   //       for (int i = 0; i < ne; i++)
+   //       {
+   //          psize[partitioning[i]].one++;
+   //       }
+
+   //       empty_parts = 0;
+   //       for (int i = 0; i < npart; i++)
+   //       {
+   //          if (psize[i].one == 0) { empty_parts++; }
+   //       }
+   //    };
+
+   //    count_partition_elements();
+
+   //    // This code just split the largest partitionings in two.
+   //    // Do we need to replace it with something better?
+   //    while (empty_parts)
+   //    {
+   //       if (print_messages)
+   //       {
+   //          mfem::err << "Mesh::GeneratePartitioning(...): METIS returned "
+   //                    << empty_parts << " empty parts!"
+   //                    << " Applying a simple fix ..." << endl;
+   //       }
+
+   //       SortPairs<int,int>(psize, npart);
+
+   //       for (int i = npart-1; i > npart-1-empty_parts; i--)
+   //       {
+   //          psize[i].one /= 2;
+   //       }
+
+   //       for (int j = 0; j < ne; j++)
+   //       {
+   //          for (int i = npart-1; i > npart-1-empty_parts; i--)
+   //          {
+   //             if (psize[i].one == 0 || partitioning[j] != psize[i].two)
+   //             {
+   //                continue;
+   //             }
+   //             else
+   //             {
+   //                partitioning[j] = psize[npart-1-i].two;
+   //                psize[i].one--;
+   //             }
+   //          }
+   //       }
+
+   //       // Check for empty partitionings again
+   //       count_partition_elements();
+   //    }
+   // }
+
+   return p;
+}
+
+}

miniapps/dg-agglomeration/partition.hpp

@@ -0,0 +1,27 @@
+// Copyright (c) 2010-2025, 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_PARTITION_HPP
+#define MFEM_PARTITION_HPP
+
+#include "mfem.hpp"
+#include "optref.hpp"
+
+namespace mfem
+{
+
+Array<int> PartitionMesh(Mesh &mesh,
+                         const int npart, const int part_method,
+                         OptRef<Array<int>> subset = NullOpt);
+
+} // namespace mfem
+
+#endif

miniapps/dg-agglomeration/smooth_agg_gmg_ex.cpp

@@ -0,0 +1,241 @@
+// Copyright (c) 2010-2025, 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.
+//
+//                     ------------------------
+//                     DG Smooth-Agg GMG Solver
+//                     ------------------------
+
+#include "mfem.hpp"
+#include <iostream>
+#include <memory>
+
+#include "mg_agglom.hpp"
+
+using namespace std;
+using namespace mfem;
+
+// RHS
+real_t rhs_function(const Vector &x);
+
+// true solution
+real_t u_true(const Vector &x);
+
+int main(int argc, char *argv[])
+{
+   const char *mesh_file = "../../data/inline-hex.mesh";
+   int order = 1;
+   real_t kappa_0 = 1.0;
+   int num_levels = 2;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file.");
+   // args.AddOption(&ref_levels, "-r", "--refine", "Refinement levels.");
+   args.AddOption(&order, "-o", "--order", "Polynomial degree.");
+   args.AddOption(&kappa_0, "-k", "--kappa", "DG penalty parameter.");
+   // args.AddOption(&ncoarse, "-nc", "--ncoarse", "Number of Fine Elements per Coarse.");
+   args.AddOption(&num_levels, "-nl", "--levels", "Number of Multigrid Levels.");
+   args.ParseCheck();
+
+   Mesh mesh(mesh_file);
+   const int dim = mesh.Dimension();
+   int ncoarse = pow(2, dim); 
+   int ref_levels = num_levels - 2;
+
+
+   for (int i = 0; i < ref_levels; ++i) { mesh.UniformRefinement(); }
+
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fespace(&mesh, &fec);
+   cout << "Number of unknowns: " << fespace.GetVSize() << endl;
+
+   int ne = mesh.GetNE();
+
+   const real_t sigma = -1.0;
+   const real_t kappa = kappa_0 * (order + 1) * (order + 1) / 2;
+
+   // Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   // ess_bdr = 0;
+   // ess_bdr[0] = 1;
+
+   // // output a paraview visualization of the partition, if desired
+   // bool paraview_vis = false;
+   // if (paraview_vis)
+   // {
+   //    L2_FECollection l2_fec(0, mesh.Dimension());
+   //    FiniteElementSpace l2_fes(&mesh, &l2_fec);
+   //    GridFunction p_gf(&l2_fes);
+   //    ParaViewDataCollection pv("Indexing", &mesh);
+   //    pv.SetPrefixPath("ParaView");
+   //    pv.RegisterField("p", &p_gf);
+   //    for (int i = 0; i < ne; ++i)
+   //    {
+   //       p_gf[i] = i;
+   //       pv.SetCycle(0);
+   //       pv.SetTime(0);
+   //       pv.Save();
+   //    }
+   // }
+
+   // Array<int> dof_indices;
+   // for (int e = 0; e < ne; ++e)
+   // {
+   //    fespace.GetElementDofs(e, dof_indices);
+   //        std::cout << "Element " << e << " DoF indices: ";
+   //    for (int j = 0; j < dof_indices.Size(); ++j) {
+   //       std::cout << dof_indices[j] << " ";
+   //    }
+   //    std::cout << std::endl;
+   // }
+   // std::string file_name_t2t = "../../../adaptiveMG/t2t_mfem.txt";
+
+   // {
+   //    const auto &e2e = mesh.ElementToElementTable();
+   //    Array<int> fn;
+
+   //    std::ofstream f(file_name_t2t);
+   //    for (int e = 0; e < mesh.GetNE(); ++e)
+   //    {
+   //       e2e.GetRow(e, fn);
+   //       int i = 0;
+   //       for (; i < fn.Size(); ++i)
+   //       {
+   //          if (fn[i] != e)
+   //          {
+   //             f << (fn[i] + 1) << " ";
+   //          }
+   //       }
+   //       const int nf = dim == 2 ? Geometry::NumEdges[mesh.GetElementGeometry(e)]
+   //                      : Geometry::NumFaces[mesh.GetElementGeometry(e)];
+   //       for (; i < nf; ++i)
+   //       {
+   //          f << -1 << " ";
+   //       }
+   //       f << '\n';
+   //    }
+   // }
+
+   // mfem::Array<int> dbc_marker(6);
+   // dbc_marker = 1;
+   // dbc_marker[4] = 0;
+
+   // LinearForm b(&fespace);
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient mone(-1.0);
+   ConstantCoefficient zero(0.0);
+   FunctionCoefficient rhs(rhs_function);
+   // b.AddDomainIntegrator(new DomainLFIntegrator(rhs));
+   // b.AddBdrFaceIntegrator(
+   //    new DGDirichletLFIntegrator(zero, one, sigma, kappa), dbc_marker); //make this -1 on just the left bc
+   // b.Assemble();
+
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a.Assemble();
+   a.Finalize();
+
+   SparseMatrix &A = a.SpMat();
+
+   // std::string file_name = "A_mfem.mtx";
+   // std::ofstream ofs1("../../../adaptiveMG/" + file_name);
+   // A.PrintMM(ofs1); 
+   // ofs1.close();
+
+   {
+      std::ofstream f("A.txt");
+      A.PrintMatlab(f);
+   }
+   {
+      const auto &e2e = mesh.ElementToElementTable();
+      Array<int> fn;
+
+      std::ofstream f("t2t.txt");
+      for (int e = 0; e < mesh.GetNE(); ++e)
+      {
+         e2e.GetRow(e, fn);
+         int i = 0;
+         for (; i < fn.Size(); ++i)
+         {
+            if (fn[i] != e)
+            {
+               f << (fn[i] + 1) << " ";
+            }
+         }
+         const int nf = dim == 2 ? Geometry::NumEdges[mesh.GetElementGeometry(e)]
+                        : Geometry::NumFaces[mesh.GetElementGeometry(e)];
+         for (; i < nf; ++i)
+         {
+            f << -1 << " ";
+         }
+         f << '\n';
+      }
+   }
+
+   SmoothedAggregationGMG mg(fespace, A, ncoarse, num_levels, false);
+   mg.SetCycleType(mfem::MultigridBase::CycleType::VCYCLE, 3, 3);
+
+   CGSolver cg;
+   cg.SetRelTol(1e-7);
+   cg.SetMaxIter(500);
+   cg.SetPrintLevel(1);
+   cg.SetOperator(A);
+   cg.SetPreconditioner(mg);
+   Vector bb(fespace.GetVSize());
+   bb = 1.0;
+   x = 0.0;
+   cg.Mult(bb, x);
+
+   // FunctionCoefficient true_solution(u_true);
+   // GridFunction true_sol_gf(&fespace);
+   // true_sol_gf.ProjectCoefficient(true_solution);
+   // double l2_error = x.ComputeL2Error(true_solution);
+
+   // std::cout << "True  rel L2 Error: " << l2_error/true_sol_gf.Norml2() << std::endl;
+
+   return 0;
+}
+
+// Initial condition
+real_t rhs_function(const Vector &x)
+{
+   // int dim = x.Size();
+
+   real_t px = M_PI*x(0);
+   real_t py = M_PI*x(1);
+   real_t pz = M_PI*x(2);
+
+   real_t pi_s = M_PI*M_PI;
+
+   real_t sss = sin(px)*sin(py)*sin(pz);
+
+   real_t css = cos(px)*sin(py)*sin(pz)*cos(px)*sin(py)*sin(pz);
+   real_t scs = sin(px)*cos(py)*sin(pz)*sin(px)*cos(py)*sin(pz);
+   real_t ssc = sin(px)*sin(py)*cos(pz)*sin(px)*sin(py)*cos(pz);
+   return -pi_s*exp(sss)*(-3*sss + css + scs + ssc);
+}
+
+
+// Initial condition
+real_t u_true(const Vector &x)
+{
+   // int dim = x.Size();
+
+   real_t px = M_PI*x(0);
+   real_t py = M_PI*x(1);
+   real_t pz = M_PI*x(2);
+
+   real_t sss = sin(px)*sin(py)*sin(pz);
+   return exp(sss) - 1;
+}