anderson accelation

miniapps/mtop/CMakeLists.txt

@@ -27,4 +27,20 @@ if(MFEM_USE_MPI AND MFEM_USE_DOUBLE)
 					 ${PAR_MTOP_COMMON_FILES}
 					 LIBRARIES mfem)
 
+	add_mfem_miniapp(mtop_test_dfem_elast
+				MAIN	mtop_test_dfem_elast.cpp
+				${PAR_MTOP_COMMON_FILES}
+                                 LIBRARIES mfem)
+
+	add_mfem_miniapp(mtop_test_filt
+				MAIN mtop_test_filt.cpp
+				${PAR_MTOP_COMMON_FILES}
+				LIBRARIES mfem)
+
+	add_mfem_miniapp(anderson_fp_test
+                                MAIN anderson_fp_test.cpp
+                                ${PAR_MTOP_COMMON_FILES}
+                                LIBRARIES mfem)
+
+
 endif (MFEM_USE_MPI AND MFEM_USE_DOUBLE)

miniapps/mtop/anderson_fp_solver.hpp

@@ -0,0 +1,552 @@
+
+#ifndef MFEM_ANDERSON_FP_SOLVER_HPP
+#define MFEM_ANDERSON_FP_SOLVER_HPP
+#include "mfem.hpp"
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#endif
+
+#include <vector>
+#include <cmath>
+#include <algorithm>
+#include <numeric>
+#include <iomanip>
+
+// Parallel Anderson acceleration for fixed point x = G(x)
+// LS via Gram matrix: (A^T A + lambda I) gamma = A^T f
+// All dot products use MFEM's device-capable inner product (x * y).
+//
+// Update form (difference form):
+//   f_k = G(x_k) - x_k
+//   A = [Δf_{k-p}, ..., Δf_{k-1}]   (n x p)
+//   solve gamma = argmin || f_k - A gamma ||_2
+//   x_{k+1} = x_k + beta * ( f_k - (Δg) gamma )
+//
+// Notes:
+// - LS communication: ONE MPI_Allreduce on packed (G,b).
+// - Local dots use x*y (Vector::operator*(Vector)), device-capable when vectors are UseDevice(true).
+
+class AndersonFixedPointSolverParGramDeviceIP : public mfem::Solver
+{
+public:
+#ifdef MFEM_USE_MPI
+   AndersonFixedPointSolverParGramDeviceIP(MPI_Comm comm, int m = 5)
+      : mfem::Solver(0, /*iter_mode=*/false), comm_(comm), m_(m)
+   {
+      MFEM_VERIFY(m_ >= 0, "Anderson: depth m must be >= 0");
+      MPI_Comm_rank(comm_, &myid_);
+      MPI_Comm_size(comm_, &nprocs_);
+   }
+#else
+   AndersonFixedPointSolverParGramDeviceIP(int m = 5)
+      : mfem::Solver(0, /*iter_mode=*/false), m_(m)
+   {
+      MFEM_VERIFY(m_ >= 0, "Anderson: depth m must be >= 0");
+      myid_ = 0; nprocs_ = 1;
+   }
+#endif
+
+   void SetOperator(const mfem::Operator &op) override
+   {
+      Gmap_ = &op;
+      height = op.Height();
+      width  = op.Width();
+      MFEM_VERIFY(height == width, "Anderson: Operator must be square.");
+
+      const int n = height;
+
+      xk_.SetSize(n);
+      xkp1_.SetSize(n);
+      gx_.SetSize(n);
+      fk_.SetSize(n);
+
+      g_prev_.SetSize(n);
+      f_prev_.SetSize(n);
+
+      df_.SetSize(n);
+      dg_.SetSize(n);
+      corr_.SetSize(n);
+
+      // Enable device semantics for all vectors that participate in dot-products / axpys.
+      // (This does not force GPU; it enables MFEM Device execution where available.) :contentReference[oaicite:3]{index=3}
+      EnableDeviceVectors_();
+
+      // Ring buffers
+      const int cap = std::max(1, m_);
+      dF_.resize(cap);
+      dG_.resize(cap);
+      for (int i = 0; i < cap; ++i)
+      {
+         dF_[i].SetSize(n);
+         dG_[i].SetSize(n);
+         dF_[i].UseDevice(true);
+         dG_[i].UseDevice(true);
+      }
+
+      // Small dense work (max m x m)
+      Gsmall_.SetSize(cap, cap);
+      Asmall_.SetSize(cap, cap);
+      Vsmall_.SetSize(cap, cap);
+
+      eval_.assign(cap, mfem::real_t(0));
+
+      bsmall_.SetSize(cap);
+      ysmall_.SetSize(cap);
+      gamma_.SetSize(cap);
+
+      ResetHistory_();
+   }
+
+   void Mult(const mfem::Vector &x0, mfem::Vector &x) const override
+   {
+      MFEM_VERIFY(Gmap_ != nullptr, "Anderson: call SetOperator() first.");
+      MFEM_VERIFY(x0.Size() == height, "Anderson: bad input size.");
+
+      ResetHistory_();
+
+      if (iterative_mode) { xk_ = x; }
+      else                { xk_ = x0; }
+
+      double res0 = -1.0;
+      double res  = 0.0;
+
+      if (print_level_ > 0 && myid_ == 0)
+      {
+         mfem::out << "AndersonFixedPointSolverParGramDeviceIP: m=" << m_
+                   << " beta=" << beta_
+                   << " max_it=" << max_iter_
+                   << " rtol=" << rel_tol_
+                   << " atol=" << abs_tol_
+                   << " reg_rel=" << reg_rel_
+                   << " rcond=" << rcond_
+                   << "\n";
+      }
+
+      for (int it = 0; it < max_iter_; ++it)
+      {
+         // gx = G(xk)
+         Gmap_->Mult(xk_, gx_);
+
+         // fk = gx - xk
+         fk_ = gx_;
+         fk_ -= xk_;
+
+         // global ||fk||
+         res = std::sqrt((double)DotGlobal_(fk_, fk_));
+         if (it == 0) { res0 = res; }
+
+         if (print_level_ > 0 && myid_ == 0)
+         {
+            mfem::out << "  it " << std::setw(4) << it
+                      << "  ||G(x)-x|| = " << std::scientific << res
+                      << "  depth=" << p_
+                      << "\n";
+         }
+
+         const double tol = std::max(abs_tol_, rel_tol_ * res0);
+         if (res <= tol)
+         {
+            final_iter_ = it;
+            final_norm_ = res;
+            x.SetSize(height, xk_); // keep memory type consistent with xk_
+            x = xk_;
+            return;
+         }
+
+         // No prev or m=0 -> damped fixed point
+         if (!has_prev_ || m_ == 0)
+         {
+            xkp1_ = xk_;
+            xkp1_.Add(beta_, fk_);
+
+            g_prev_ = gx_;
+            f_prev_ = fk_;
+            has_prev_ = true;
+
+            xk_ = xkp1_;
+            continue;
+         }
+
+         // df = fk - f_prev, dg = gx - g_prev
+         df_ = fk_; df_ -= f_prev_;
+         dg_ = gx_; dg_ -= g_prev_;
+
+         PushHistory_(df_, dg_); // updates p_ (<=m_) and ring start_
+
+         // Solve LS via Gram system
+         const bool ok = SolveLeastSquares_GramEigen_();
+         if (!ok)
+         {
+            if (print_level_ > 0 && myid_ == 0)
+            {
+               mfem::out << "  LS solve failed -> restarting history, plain step.\n";
+            }
+            ResetHistory_();
+
+            xkp1_ = xk_;
+            xkp1_.Add(beta_, fk_);
+
+            g_prev_ = gx_;
+            f_prev_ = fk_;
+            has_prev_ = true;
+
+            xk_ = xkp1_;
+            continue;
+         }
+
+         // corr = dG * gamma
+         corr_ = 0.0;
+         for (int j = 0; j < p_; ++j) { corr_.Add(gamma_(j), DGcol_(j)); }
+
+         // x_{k+1} = x_k + beta * ( f_k - corr )
+         xkp1_ = xk_;
+         xkp1_.Add(beta_, fk_);
+         xkp1_.Add(-beta_, corr_);
+
+         // update previous
+         g_prev_ = gx_;
+         f_prev_ = fk_;
+         xk_ = xkp1_;
+      }
+
+      final_iter_ = max_iter_;
+      final_norm_ = res;
+      x.SetSize(height, xk_);
+      x = xk_;
+   }
+
+   // ---------------- parameters ----------------
+   void SetMaxIter(int max_it) { max_iter_ = max_it; }
+   void SetRelTol(double rtol) { rel_tol_ = rtol; }
+   void SetAbsTol(double atol) { abs_tol_ = atol; }
+   void SetBeta(double beta)   { beta_ = beta; }
+   void SetPrintLevel(int pl)  { print_level_ = pl; }
+
+   // lambda = reg_rel * trace(G)/p
+   void SetRegularizationRel(double reg_rel) { reg_rel_ = reg_rel; }
+
+   // keep eigenvalues >= rcond * lambda_max (set rcond<0 to disable)
+   void SetRcond(double rcond) { rcond_ = rcond; }
+
+   void SetDepth(int m)
+   {
+      MFEM_VERIFY(m >= 0, "Anderson: m must be >= 0");
+      m_ = m;
+      if (Gmap_ != nullptr) { SetOperator(*Gmap_); }
+   }
+
+   int    GetNumIterations() const { return final_iter_; }
+   double GetFinalNorm()     const { return final_norm_; }
+
+private:
+   const mfem::Operator *Gmap_ = nullptr;
+
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm_ = MPI_COMM_WORLD;
+#endif
+   int myid_ = 0;
+   int nprocs_ = 1;
+
+   int m_ = 5;
+   int max_iter_ = 50;
+   double rel_tol_ = 1e-8;
+   double abs_tol_ = 0.0;
+   double beta_ = 1.0;
+
+   double reg_rel_ = 1e-12;
+   double rcond_   = 1e-12;
+
+   int print_level_ = 0;
+
+   // History ring buffer state
+   mutable bool has_prev_ = false;
+   mutable int p_ = 0;
+   mutable int start_ = 0;
+
+   // Iteration vectors
+   mutable mfem::Vector xk_, xkp1_, gx_, fk_;
+   mutable mfem::Vector g_prev_, f_prev_;
+   mutable mfem::Vector df_, dg_, corr_;
+
+   // Stored columns Δf, Δg
+   mutable std::vector<mfem::Vector> dF_, dG_;
+
+   // Small dense work buffers (allocated as max m x m, used as p x p)
+   mutable mfem::DenseMatrix Gsmall_, Asmall_, Vsmall_;
+   mutable std::vector<mfem::real_t> eval_;
+   mutable mfem::Vector bsmall_, ysmall_, gamma_;
+
+   // Stats
+   mutable int final_iter_ = 0;
+   mutable double final_norm_ = 0.0;
+
+private:
+   void EnableDeviceVectors_() const
+   {
+      xk_.UseDevice(true);
+      xkp1_.UseDevice(true);
+      gx_.UseDevice(true);
+      fk_.UseDevice(true);
+
+      g_prev_.UseDevice(true);
+      f_prev_.UseDevice(true);
+
+      df_.UseDevice(true);
+      dg_.UseDevice(true);
+      corr_.UseDevice(true);
+   }
+
+   void ResetHistory_() const
+   {
+      has_prev_ = false;
+      p_ = 0;
+      start_ = 0;
+   }
+
+   int RingIndex_(int pos) const
+   {
+      return (start_ + pos) % std::max(1, m_);
+   }
+
+   const mfem::Vector& DFcol_(int pos) const { return dF_[RingIndex_(pos)]; }
+   const mfem::Vector& DGcol_(int pos) const { return dG_[RingIndex_(pos)]; }
+
+   void PushHistory_(const mfem::Vector &df, const mfem::Vector &dg) const
+   {
+      if (m_ == 0) { return; }
+      const int M = m_;
+
+      if (p_ < M)
+      {
+         const int idx = (start_ + p_) % M;
+         dF_[idx] = df;
+         dG_[idx] = dg;
+         ++p_;
+      }
+      else
+      {
+         const int idx = start_;
+         dF_[idx] = df;
+         dG_[idx] = dg;
+         start_ = (start_ + 1) % M;
+         p_ = M;
+      }
+   }
+
+   // --------- device-capable dot products ----------
+   static mfem::real_t DotLocal_(const mfem::Vector &a,
+                                                   const mfem::Vector &b)
+   {
+      // This is x*y (Vector::operator*(Vector)), which is the device-capable inner product. :contentReference[oaicite:4]{index=4}
+      MFEM_ASSERT(a.Size() == b.Size(), "Dot: size mismatch.");
+      return a * b;
+   }
+
+   mfem::real_t DotGlobal_(const mfem::Vector &a, const mfem::Vector &b) const
+   {
+#ifdef MFEM_USE_MPI
+      // MFEM provides an inline MPI_Comm overload: loc = a*b; MPI_Allreduce(...). :contentReference[oaicite:5]{index=5}
+      return mfem::InnerProduct(comm_, a, b);
+#else
+      return mfem::InnerProduct(a, b);
+#endif
+   }
+
+   // Small symmetric Jacobi eigensolver for p x p matrix (p<=~20 typical).
+   static void JacobiEigenSymmetric_(mfem::DenseMatrix &A, mfem::DenseMatrix &V,
+                                    mfem::real_t *eval, int p,
+                                    int max_sweeps = 50, mfem::real_t tol = 1e-14)
+   {
+      for (int i = 0; i < p; ++i)
+      {
+         for (int j = 0; j < p; ++j) { V(i,j) = (i == j) ? 1.0 : 0.0; }
+      }
+
+      auto offdiag_max = [&]() {
+         mfem::real_t mx = 0.0;
+         for (int i = 0; i < p; ++i)
+            for (int j = i+1; j < p; ++j)
+               mx = std::max(mx, (mfem::real_t)std::abs(A(i,j)));
+         return mx;
+      };
+
+      for (int sweep = 0; sweep < max_sweeps; ++sweep)
+      {
+         const mfem::real_t mx = offdiag_max();
+         if (mx < tol) { break; }
+
+         for (int q = 1; q < p; ++q)
+         {
+            for (int r = 0; r < q; ++r)
+            {
+               const mfem::real_t a_rq = A(r,q);
+               if (std::abs(a_rq) < tol) { continue; }
+
+               const mfem::real_t a_rr = A(r,r);
+               const mfem::real_t a_qq = A(q,q);
+
+               const mfem::real_t tau = (a_qq - a_rr) / (2.0 * a_rq);
+               mfem::real_t t;
+               if (tau >= 0.0) { t =  1.0 / (tau + (mfem::real_t)std::sqrt(1.0 + tau*tau)); }
+               else            { t = -1.0 / (-tau + (mfem::real_t)std::sqrt(1.0 + tau*tau)); }
+
+               const mfem::real_t c = 1.0 / (mfem::real_t)std::sqrt(1.0 + t*t);
+               const mfem::real_t s = t * c;
+
+               for (int k = 0; k < p; ++k)
+               {
+                  if (k == r || k == q) { continue; }
+                  const mfem::real_t a_kr = A(k,r);
+                  const mfem::real_t a_kq = A(k,q);
+
+                  const mfem::real_t nr = c*a_kr - s*a_kq;
+                  const mfem::real_t nq = s*a_kr + c*a_kq;
+
+                  A(k,r) = A(r,k) = nr;
+                  A(k,q) = A(q,k) = nq;
+               }
+
+               const mfem::real_t a_rr_new = c*c*a_rr - 2.0*s*c*a_rq + s*s*a_qq;
+               const mfem::real_t a_qq_new = s*s*a_rr + 2.0*s*c*a_rq + c*c*a_qq;
+
+               A(r,r) = a_rr_new;
+               A(q,q) = a_qq_new;
+               A(r,q) = A(q,r) = 0.0;
+
+               for (int k = 0; k < p; ++k)
+               {
+                  const mfem::real_t v_kr = V(k,r);
+                  const mfem::real_t v_kq = V(k,q);
+                  V(k,r) = c*v_kr - s*v_kq;
+                  V(k,q) = s*v_kr + c*v_kq;
+               }
+            }
+         }
+      }
+
+      for (int i = 0; i < p; ++i) { eval[i] = A(i,i); }
+   }
+
+   bool SolveLeastSquares_GramEigen_() const
+   {
+      const int p = p_;
+      if (p <= 0) { return false; }
+
+      // pack upper triangle of G plus b: size = p*(p+1)/2 + p
+      const int tri = p*(p+1)/2;
+      std::vector<mfem::real_t> pack(tri + p, mfem::real_t(0));
+
+      // Local assembly: all dots are local, device-capable x*y (no MPI yet)
+      int idx = 0;
+      for (int j = 0; j < p; ++j)
+      {
+         const mfem::Vector &aj = DFcol_(j);
+         for (int i = 0; i <= j; ++i)
+         {
+            pack[idx++] = DotLocal_(DFcol_(i), aj);
+         }
+      }
+      for (int j = 0; j < p; ++j)
+      {
+         pack[idx++] = DotLocal_(DFcol_(j), fk_);
+      }
+
+#ifdef MFEM_USE_MPI
+      // One batched allreduce
+      MPI_Allreduce(MPI_IN_PLACE, pack.data(), (int)pack.size(),
+                    MFEM_MPI_REAL_T, MPI_SUM, comm_);
+#endif
+
+      // Unpack into Gsmall and b
+      idx = 0;
+      for (int j = 0; j < p; ++j)
+      {
+         for (int i = 0; i <= j; ++i)
+         {
+            const mfem::real_t val = pack[idx++];
+            Gsmall_(i,j) = val;
+            Gsmall_(j,i) = val;
+         }
+      }
+      for (int j = 0; j < p; ++j) { bsmall_(j) = pack[idx++]; }
+
+      // lambda = reg_rel * trace(G)/p
+      double trace = 0.0;
+      for (int i = 0; i < p; ++i) { trace += (double)Gsmall_(i,i); }
+      const double lambda = (reg_rel_ > 0.0) ? (reg_rel_ * (trace / (double)p)) : 0.0;
+
+      // Asmall = G + lambda I
+      for (int i = 0; i < p; ++i)
+      {
+         for (int j = 0; j < p; ++j) { Asmall_(i,j) = Gsmall_(i,j); }
+      }
+      if (lambda > 0.0)
+      {
+         for (int i = 0; i < p; ++i) { Asmall_(i,i) += (mfem::real_t)lambda; }
+      }
+
+      // Eigendecomposition (Jacobi)
+      JacobiEigenSymmetric_(Asmall_, Vsmall_, eval_.data(), p);
+
+      // Sort eigenpairs descending
+      std::vector<int> perm(p);
+      std::iota(perm.begin(), perm.end(), 0);
+      std::stable_sort(perm.begin(), perm.end(),
+                       [&](int a, int b) { return eval_[a] > eval_[b]; });
+
+      mfem::DenseMatrix Vsorted(p, p);
+      std::vector<mfem::real_t> eval_sorted(p, mfem::real_t(0));
+      for (int j = 0; j < p; ++j)
+      {
+         const int oj = perm[j];
+         eval_sorted[j] = eval_[oj];
+         for (int i = 0; i < p; ++i) { Vsorted(i,j) = Vsmall_(i,oj); }
+      }
+      for (int j = 0; j < p; ++j)
+      {
+         eval_[j] = eval_sorted[j];
+         for (int i = 0; i < p; ++i) { Vsmall_(i,j) = Vsorted(i,j); }
+      }
+
+      mfem::real_t lambda_max = 0.0;
+      for (int i = 0; i < p; ++i) { lambda_max = std::max(lambda_max, eval_[i]); }
+
+      if (lambda_max <= 0.0)
+      {
+         for (int i = 0; i < p; ++i) { gamma_(i) = 0.0; }
+         return true;
+      }
+
+      const mfem::real_t cutoff = (rcond_ >= 0.0) ? (mfem::real_t)(rcond_ * (double)lambda_max) : mfem::real_t(0);
+
+      // y = V^T b
+      for (int i = 0; i < p; ++i)
+      {
+         mfem::real_t s = 0.0;
+         for (int k = 0; k < p; ++k) { s += Vsmall_(k,i) * bsmall_(k); }
+         ysmall_(i) = s;
+      }
+
+      // z_i = y_i / eval_i with truncation
+      for (int i = 0; i < p; ++i)
+      {
+         const mfem::real_t ei = eval_[i];
+         if (ei < cutoff || ei <= 0.0) { ysmall_(i) = 0.0; }
+         else                          { ysmall_(i) = ysmall_(i) / ei; }
+      }
+
+      // gamma = V z
+      for (int k = 0; k < p; ++k)
+      {
+         mfem::real_t s = 0.0;
+         for (int i = 0; i < p; ++i) { s += Vsmall_(k,i) * ysmall_(i); }
+         gamma_(k) = s;
+      }
+      for (int k = p; k < gamma_.Size(); ++k) { gamma_(k) = 0.0; }
+
+      return true;
+   }
+};
+
+#endif // MFEM_ANDERSON_FP_SOLVER_HPP

miniapps/mtop/anderson_fp_test.cpp

@@ -0,0 +1,236 @@
+#include "anderson_fp_solver.hpp"
+
+#include "mtop_solvers.hpp"
+
+using namespace std;
+using namespace mfem;
+
+constexpr auto MESH_TRI = MFEM_SOURCE_DIR "/miniapps/mtop/sq_2D_9_tri.mesh";
+constexpr auto MESH_QUAD = MFEM_SOURCE_DIR "/miniapps/mtop/sq_2D_9_quad.mesh";
+
+
+/// Fixed-point operator for testing Anderson acceleration
+class FPOperator : public mfem::Operator
+{
+public:
+   FPOperator(const HypreParMatrix& A_, const Vector& rhs_) :
+                        mfem::Operator(A_.Height()), A(&A_), rhs(rhs_)  
+   {
+       Vector diag(A_.Width());
+       tmp.SetSize(A_.Width()); 
+       tmp=1.0;
+       //A_.AbsMult(tmp, diag);
+       A->GetDiag(diag);
+       smoother = std::make_unique<OperatorJacobiSmoother>(diag, Array<int>(), 1.0);
+
+   }
+
+   virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const override
+   {
+      A->Mult(x, tmp);
+      tmp.Neg();
+      tmp.Add(1.0, rhs);
+      smoother->Mult(tmp, y); y.Add(1.0, x);
+   }
+
+private:
+   const HypreParMatrix *A;
+   const Vector& rhs;
+   mutable Vector tmp;
+   std::unique_ptr<OperatorJacobiSmoother> smoother;
+
+};
+
+
+class DensCoeff : public mfem::Coefficient
+{
+private:
+    real_t l;
+public:
+    DensCoeff(real_t d=1.0) : l(d) {}
+
+    virtual real_t Eval(mfem::ElementTransformation &T,
+                        const mfem::IntegrationPoint &ip)
+    {
+        Vector x;
+        T.Transform(ip, x);
+        real_t r = x.Norml2();
+        r=sin(M_PI*r/l);
+        if(r>0.5)
+            r=1.0;
+        else
+            r=0.0;
+        return r;
+    }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init();
+   Hypre::Init();
+
+   // Parse command-line options.
+    const char *mesh_file = MESH_QUAD;
+   const char *device_config = "cpu";
+   int order = 2;
+   bool pa = false;
+   bool dfem = false;
+   bool mesh_tri = false;
+   bool mesh_quad = false;
+   int par_ref_levels = 1;
+   bool paraview = false;
+   bool visualization = true;
+   int m=1;
+   real_t beta=1.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&dfem, "-dfem", "--dFEM", "-no-dfem", "--no-dFEM",
+                  "Enable or not dFEM.");
+    args.AddOption(&mesh_tri, "-tri", "--triangular", "-no-tri",
+                  "--no-triangular", "Enable or not triangular mesh.");
+    args.AddOption(&mesh_quad, "-quad", "--quadrilateral", "-no-quad",
+                  "--no-quadrilateral", "Enable or not quadrilateral mesh.");
+    args.AddOption(&par_ref_levels, "-prl", "--par-ref-levels",
+                  "Number of times to refine the mesh uniformly in parallel.");
+    args.AddOption(&paraview, "-pv", "--paraview", "-no-pv", "--no-paraview",
+                  "Enable or not Paraview visualization");
+    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+    args.AddOption(&m, "-ma", "--m-accel",
+                  "Anderson acceleration parameter m.");
+    args.AddOption(&beta, "-beta", "--beta",
+                  "Anderson acceleration relaxation parameter beta.");
+    args.ParseCheck();
+
+    // Enable hardware devices such as GPUs, and programming models such as
+    // CUDA, OCCA, RAJA and OpenMP based on command line options.
+    Device device(device_config);
+    if (Mpi::Root()) { device.Print(); }
+
+    // Read the (serial) mesh from the given mesh file on all processors.  We
+    // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+    // and volume meshes with the same code.
+    Mesh mesh(mesh_tri ? MESH_TRI : mesh_quad ? MESH_QUAD : mesh_file, 1, 1);
+    const int dim = mesh.Dimension();
+
+    // Refine the serial mesh on all processors 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 1000 elements.
+    {
+       const int ref_levels =
+          (int)floor(log(1000. / mesh.GetNE()) / log(2.) / dim);
+       for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
+    }
+    if (Mpi::Root())
+    {
+       std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
+    }
+
+    // Define a parallel mesh by a partitioning of the serial mesh. Refine
+    // this mesh further in parallel to increase the resolution. Once the
+    // parallel mesh is defined, the serial mesh can be deleted.
+    ParMesh pmesh(MPI_COMM_WORLD, mesh);
+    mesh.Clear();
+    for (int l = 0; l < par_ref_levels; l++) { pmesh.UniformRefinement(); }   
+
+    // Define a finite element space on the mesh. Here we use continuous
+    // Lagrange finite elements of the specified order.
+    H1_FECollection fec(order,dim);
+    ParFiniteElementSpace fespace(&pmesh, &fec, 1);
+
+    // Define the solution vector x as a finite element grid function
+    ParGridFunction x(&fespace); x=0.0;
+
+    std::unique_ptr<HypreParMatrix> A;
+    Vector rhs; rhs.SetSize(fespace.GetTrueVSize());
+    {
+        // Set up the linear system Ax=b
+        ParBilinearForm a(&fespace);
+        ConstantCoefficient diff_coeff(0.01);
+        a.AddDomainIntegrator(new DiffusionIntegrator(diff_coeff));
+        a.AddDomainIntegrator(new MassIntegrator());
+        a.Assemble();
+        a.Finalize();   
+
+        A.reset(a.ParallelAssemble());
+
+        ParLinearForm b(&fespace);
+        DensCoeff dens_coeff(1.0);
+        b.AddDomainIntegrator(new DomainLFIntegrator(dens_coeff));
+        b.Assemble();
+        b.ParallelAssemble(rhs);
+    }
+
+
+    FPOperator G(*A, rhs);
+
+    Vector res;
+    Vector x_new;
+    x_new.SetSize(fespace.GetTrueVSize());
+    res.SetSize(fespace.GetTrueVSize());
+
+    for(int i=0; i<50; i++)
+    {
+        G.Mult(x.GetTrueVector(), x_new);
+        add(-1.0, x_new, 1.0, x.GetTrueVector(), res);
+        real_t res_norm = InnerProduct(pmesh.GetComm(), res, res);
+        real_t x_norm = InnerProduct(pmesh.GetComm(), x_new, x_new);
+        if(0==pmesh.GetMyRank()){    
+            std::cout << "Iter " << i << " : Residual norm = " << res_norm << "  Solution norm = " << x_norm << std::endl;    
+        }
+        //x.SetFromTrueDofs(x_new); x.SetTrueVector();
+        x.GetTrueVector()=x_new;
+    }
+
+    AndersonFixedPointSolverParGramDeviceIP aa(pmesh.GetComm(),/*m=*/m);
+    aa.SetOperator(G);
+
+    // Typical knobs:
+    aa.SetMaxIter(200);
+    aa.SetRelTol(1e-10);
+    aa.SetAbsTol(0.0);
+    aa.SetBeta(beta);
+
+    // Default recommendation knobs:
+    aa.SetRegularizationRel(1e-12); // try 1e-10 .. 1e-6 if coefficients blow up
+    aa.SetRcond(1e-12);             // try 1e-10 if history gets near-dependent
+
+    aa.SetPrintLevel(1);
+
+    x.GetTrueVector()=0.0;
+    aa.Mult(x.GetTrueVector(), x_new);
+    real_t x_norm = InnerProduct(pmesh.GetComm(), x_new, x_new);
+    if(0==pmesh.GetMyRank()){    
+            std::cout  <<  "  Solution norm = " << x_norm << std::endl;    
+    }
+
+    x.SetFromTrueDofs(x_new); x.SetTrueVector();
+
+    if (paraview)
+    {
+        ParaViewDataCollection paraview_dc("anderson", &pmesh);
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetLevelsOfDetail(order);
+        paraview_dc.SetDataFormat(VTKFormat::BINARY);
+        paraview_dc.SetHighOrderOutput(true);
+        paraview_dc.SetCycle(0);
+        paraview_dc.SetTime(0.0);
+        paraview_dc.RegisterField("filt", &x);
+        paraview_dc.Save();
+    }
+
+   return EXIT_SUCCESS;
+
+}

miniapps/mtop/mtop_solvers.cpp

@@ -338,6 +338,7 @@ void IsoLinElasticSolver::SetEssTDofs(Vector &bsol, Array<int> &ess_dofs)
 
 void IsoLinElasticSolver::Mult(const Vector &x, Vector &y) const
 {
+
    // the rhs x is assumed to have the contribution of the BC set in advance
    // the BC values are not modified here
    ls->Mult(x, y);
@@ -547,3 +548,607 @@ void IsoLinElasticSolver::FSolve()
    delete lf;
    lf = nullptr;
 }
+
+template <int DIM, typename scalar_t=real_t> struct FilterQFunction
+{
+   using matd_t = tensor<real_t, DIM, DIM>;
+   using vecd_t = tensor<scalar_t, DIM>;
+
+
+   struct Diffusion
+   {
+      real_t diff_coeff=1.0;
+
+      void SetDiffusion(real_t val_)
+      {
+         diff_coeff=val_;
+      }
+
+      MFEM_HOST_DEVICE inline auto operator()(const vecd_t &dfdxi,
+                                              //const real_t &diff,
+                                              const matd_t &J,
+                                              const real_t &w) const
+      {
+         const auto invJ = mfem::future::inv(J);
+         const auto TinJ = mfem::future::transpose(invJ);
+         const auto detJ = mfem::future::det(J);
+         return tuple{ (dfdxi * invJ) * TinJ * detJ * w * diff_coeff};
+      }
+   };
+
+   struct Mass
+   {
+      real_t density=1.0;
+
+      void SetDensity(real_t val_)
+      {
+         density=val_;
+      }
+
+      MFEM_HOST_DEVICE inline auto operator()(const scalar_t &frho,
+                                              //const real_t &diff,
+                                              const matd_t &J,
+                                              const real_t &w) const
+      {
+         const auto detJ = mfem::future::det(J);
+         return tuple{ density * frho * detJ * w };
+      }
+   };
+
+
+   struct RHS
+   {
+      MFEM_HOST_DEVICE inline auto operator()(const scalar_t &frho,
+                                              const scalar_t &urho,
+                                              //const real_t &diff,
+                                              const matd_t &J,
+                                              const real_t &w) const
+      {
+         const auto detJ = mfem::future::det(J);
+         return tuple{ urho* detJ * w };
+      }
+   };
+};
+
+PDEFilter::PDEFilter(ParMesh *mesh, real_t r, int order):
+   pmesh(mesh),
+   dim(mesh->Dimension()),
+   spaceDim(mesh->SpaceDimension()),
+   filter_radius(r),
+   ffec(new H1_FECollection(order, dim)),
+   ffes(new ParFiniteElementSpace(pmesh, ffec, 1, Ordering::byNODES)),
+   ifec(new L2_FECollection(order-1, dim)),
+   ifes(new ParFiniteElementSpace(pmesh, ifec, 1, Ordering::byNODES)),
+   filtered_field(ffes),
+   input_field(ifes),
+   h1_gradient(ffes),
+   prec(nullptr),
+   ls(nullptr),
+   ess_tdofv(),
+   fe(ffes->GetFE(0)),
+   nodes((pmesh->EnsureNodes(),
+          static_cast<ParGridFunction *>(pmesh->GetNodes()))),
+   mfes(nodes->ParFESpace()),
+   ir(IntRules.Get(fe->GetGeomType(),
+                   fe->GetOrder() + fe->GetOrder() + fe->GetDim() - 1)),
+   qs(*pmesh, ir),
+   diff_ps(*pmesh, ir, 1),
+   K(nullptr)
+{
+   filtered_field = 0.0;
+   input_field = 0.0;   
+   h1_gradient = 0.0;
+
+   SetLinearSolver();
+
+   Operator::width = ifes->GetTrueVSize();
+   Operator::height = ffes->GetTrueVSize();
+
+
+   if (pmesh->attributes.Size() > 0)
+   {
+      domain_attributes.SetSize(pmesh->attributes.Max());
+      domain_attributes = 1;
+   }
+}
+
+PDEFilter::PDEFilter(ParFiniteElementSpace *fespace, real_t r, int order):
+   pmesh(fespace->GetParMesh()),
+   dim(fespace->GetParMesh()->Dimension()),
+   spaceDim(fespace->GetParMesh()->SpaceDimension()),
+   filter_radius(r),
+   ffec(new H1_FECollection(order, dim)),
+   ffes(new ParFiniteElementSpace(pmesh, ffec, 1, Ordering::byNODES)),
+   ifec(nullptr),
+   ifes(new ParFiniteElementSpace(*fespace)),
+   filtered_field(ffes),
+   input_field(ifes),
+   h1_gradient(ffes),
+   prec(nullptr),
+   ls(nullptr),
+   ess_tdofv(),
+   fe(ffes->GetFE(0)),
+   nodes((pmesh->EnsureNodes(),
+          static_cast<ParGridFunction *>(pmesh->GetNodes()))),
+   mfes(nodes->ParFESpace()),
+   ir(IntRules.Get(fe->GetGeomType(),
+                   fe->GetOrder() + fe->GetOrder() + fe->GetDim() - 1)),
+   qs(*pmesh, ir),
+   diff_ps(*pmesh, ir, 1),
+   K(nullptr)
+{
+
+   filtered_field = 0.0;
+   input_field = 0.0;   
+   h1_gradient = 0.0;
+
+   SetLinearSolver();
+
+   Operator::width = ifes->GetTrueVSize();
+   Operator::height = ffes->GetTrueVSize();
+
+
+   if (pmesh->attributes.Size() > 0)
+   {
+      domain_attributes.SetSize(pmesh->attributes.Max());
+      domain_attributes = 1;
+   }
+}
+
+
+PDEFilter::~PDEFilter()
+{
+   delete prec;
+   delete ls;
+
+   delete ffes;
+   delete ffec;
+
+   delete ifes;
+   delete ifec;
+
+   delete K;
+}  
+
+void PDEFilter::SetFilterRadius(real_t r)
+{
+   filter_radius = r;
+}
+
+void PDEFilter::SetLinearSolver(real_t rtol,
+                                 real_t atol,
+                                 int miter)
+{
+   linear_rtol = rtol;
+   linear_atol = atol;
+   linear_iter = miter;
+}
+
+void PDEFilter::Assemble()
+{
+   delete prec; prec=nullptr;
+   delete ls; ls=nullptr;
+
+   // define the differentiable operator
+   // defined the matrix vector product
+   dop= std::make_unique<mfem::future::DifferentiableOperator>(
+      std::vector<mfem::future::FieldDescriptor> {{ FSol, ffes }},
+      std::vector<mfem::future::FieldDescriptor>
+      {
+         //{ DiffCoeff, &diff_ps },
+         { Coords, mfes }
+      },
+      *pmesh);
+
+   ConstantCoefficient diffusion_coeff(filter_radius*filter_radius); 
+   // sample filter coefficient on the integration points
+   diff_cv = std::make_unique<CoefficientVector>(diffusion_coeff, qs);
+
+   // set the parameters of the differentiable operator
+   //dop->SetParameters({ diff_cv.get(), nodes });
+   dop->SetParameters({nodes});
+
+   // define the q-function for dimensions 2 and 3
+   // diffusion term
+   const auto dinputs =
+      mfem::future::tuple{ mfem::future::Gradient<FSol>{},
+                           //mfem::future::Identity<DiffCoeff>{},
+                           mfem::future::Gradient<Coords>{},
+                           mfem::future::Weight{} };
+   // mass term
+   const auto minputs =
+      mfem::future::tuple{ mfem::future::Value<FSol>{},
+                           //mfem::future::Identity<DiffCoeff>{},
+                           mfem::future::Gradient<Coords>{},
+                           mfem::future::Weight{} };
+
+   // output of the diffusion term
+   const auto doutput = mfem::future::tuple{ mfem::future::Gradient<FSol>{} };
+
+   // output of the mass term
+   const auto moutput = mfem::future::tuple{ mfem::future::Value<FSol>{} };  
+
+   
+   if (2 == spaceDim)
+   {
+     typename FilterQFunction<2>::Diffusion diff_qf;
+     diff_qf.SetDiffusion(filter_radius*filter_radius);
+     dop->AddDomainIntegrator(diff_qf, dinputs, doutput, ir, 
+                               domain_attributes);
+
+     typename FilterQFunction<2>::Mass mass_qf;
+     mass_qf.SetDensity(1.0);
+     dop->AddDomainIntegrator(mass_qf, minputs, moutput, ir, 
+                                 domain_attributes);
+   }
+   else if (3 == spaceDim)
+   {
+      typename FilterQFunction<3>::Diffusion diff_qf;
+      diff_qf.SetDiffusion(filter_radius*filter_radius);
+      dop->AddDomainIntegrator(diff_qf, dinputs, doutput, ir,
+                                 domain_attributes);
+
+      typename FilterQFunction<3>::Mass mass_qf;
+      mass_qf.SetDensity(1.0);
+      dop->AddDomainIntegrator(mass_qf, minputs, moutput, ir, 
+                                 domain_attributes);
+   }
+   else { MFEM_ABORT("Space dimension not supported"); }
+   
+
+   // set BC
+   // TODO:: Do not forget to set the BCs here 
+
+   Operator *Kop;
+   dop->FormSystemOperator(ess_tdofv, Kop);
+   Kh = std::make_unique<OperatorHandle>(Kop);
+   Kc = dynamic_cast<mfem::ConstrainedOperator*>(Kop);
+
+   // delete old assembled Jacobian if it exists
+   delete K; K=nullptr;
+   //LOR Preconditioner
+   {
+      std::unique_ptr<mfem::ParLORDiscretization> lor_disc;
+      lor_disc = std::make_unique<ParLORDiscretization>(*ffes);
+      ParFiniteElementSpace &lor_space = lor_disc->GetParFESpace();
+      ParMesh &lor_mesh = *lor_space.GetParMesh();
+      lor_mesh.EnsureNodes();
+      ParGridFunction* lor_nodes=static_cast<ParGridFunction *>(lor_mesh.GetNodes());
+      ParFiniteElementSpace* lor_nodes_fes = lor_nodes->ParFESpace();
+
+      /*
+      // Get the LOR Jacobian of the differentiable operator
+      std::unique_ptr<mfem::future::DifferentiableOperator> dopd;
+      // define the differentiable operator
+      dopd= std::make_unique<mfem::future::DifferentiableOperator>(
+         std::vector<mfem::future::FieldDescriptor> {{FSol, &lor_space}},
+         //std::vector<mfem::future::FieldDescriptor> {{FSol, ffes}},
+         std::vector<mfem::future::FieldDescriptor>
+         {
+            //{ DiffCoeff, &diff_ps },
+            { Coords, lor_nodes_fes }
+         },
+         lor_mesh);
+
+      // set the parameters of the differentiable operator
+      //dopd->SetParameters({ diff_cv.get(), nodes });
+      dopd->SetParameters({lor_nodes});
+      
+      auto derivatives = std::integer_sequence<size_t, FSol> {};
+      // define the q-function for dimensions 2 and 3
+      using mfem::future::dual;
+      using dual_t = dual<real_t, real_t>;
+      if (2 == spaceDim)
+      {
+         typename FilterQFunction<2,dual_t>::Diffusion diff_qfd;
+         diff_qfd.SetDiffusion(filter_radius*filter_radius);
+         dopd->AddDomainIntegrator(diff_qfd, dinputs, doutput, ir,
+                                    domain_attributes, derivatives);
+
+         typename FilterQFunction<2,dual_t>::Mass mass_qfd;
+         mass_qfd.SetDensity(1.0);
+         dopd->AddDomainIntegrator(mass_qfd, minputs, moutput, ir, 
+                                    domain_attributes, derivatives);
+
+         //typename FilterQFunction<2,dual_t>::Diffusion diff_qfd;
+         //diff_qfd.SetDiffusion(filter_radius*filter_radius);
+         //dopd->AddDomainIntegrator(diff_qfd, dinputs, doutput, ir,
+         //                           domain_attributes, derivatives);
+      }
+      else if (3 == spaceDim)
+      {
+         typename FilterQFunction<3,dual_t>::Diffusion diff_qfd;
+         diff_qfd.SetDiffusion(filter_radius*filter_radius);
+         dopd->AddDomainIntegrator(diff_qfd, dinputs, doutput, ir,
+                                    domain_attributes, derivatives);
+
+         typename FilterQFunction<3,dual_t>::Mass mass_qfd;
+         mass_qfd.SetDensity(1.0);
+         dopd->AddDomainIntegrator(mass_qfd, minputs, moutput, ir, 
+                                    domain_attributes, derivatives);
+      }
+
+      std::shared_ptr<mfem::future::DerivativeOperator> dres_du;
+      // set parameters using grid functions
+      ParGridFunction lor_gf(&lor_space); lor_gf=0.0;
+
+      dres_du = dopd->GetDerivative(FSol, {&lor_gf}, 
+                              { lor_nodes });
+
+      // get the Jacobian
+      dres_du->Assemble(K);
+      */      
+
+      
+      ParBilinearForm bf_lor(&lor_space);
+      ConstantCoefficient diff_coeff_lor(filter_radius*filter_radius);
+      bf_lor.AddDomainIntegrator(new DiffusionIntegrator(diff_coeff_lor));
+      bf_lor.AddDomainIntegrator(new MassIntegrator());
+      bf_lor.Assemble();
+      bf_lor.Finalize();
+      K=bf_lor.ParallelAssemble();
+      
+
+      K->EliminateBC(ess_tdofv,Operator::DiagonalPolicy::DIAG_ONE);
+   }
+
+    // set the linear solver
+   ls = new CGSolver(pmesh->GetComm());
+   prec = new HypreBoomerAMG();
+   ls->SetOperator(*Kh->Ptr());
+   ls->SetPrintLevel(1); 
+   ls->SetPreconditioner(*prec);
+
+   ls->SetRelTol(linear_rtol);
+   ls->SetAbsTol(linear_atol);
+   ls->SetMaxIter(linear_iter);
+
+   // set the preconditioner for the linear solver
+   prec->SetOperator(*K);
+
+   // defined the RHS operator
+   {
+      drh = std::make_unique<mfem::future::DifferentiableOperator>(
+         std::vector<mfem::future::FieldDescriptor> {{ FSol, ffes }},
+         std::vector<mfem::future::FieldDescriptor>
+         {
+            { USol, ifes },
+            { Coords, mfes }
+         },
+         *pmesh);
+
+      // set the parameters of the differentiable operator
+      drh->SetParameters({&input_field,nodes});
+      // define the q-function for dimensions 2 and 3
+      // mass term
+      const auto rhsinputs =
+         mfem::future::tuple{ mfem::future::Value<FSol>{},
+                              mfem::future::Value<USol>{},
+                              mfem::future::Gradient<Coords>{},
+                              mfem::future::Weight{} };
+      // output of the mass term
+      const auto rhsoutput = mfem::future::tuple{ mfem::future::Value<FSol>{} };
+
+      if (2 == spaceDim)
+      {
+         typename FilterQFunction<2>::RHS rhs_qf;
+         drh->AddDomainIntegrator(rhs_qf, rhsinputs, rhsoutput, ir,
+                                    domain_attributes);
+      }
+      else if (3 == spaceDim)
+      {
+         typename FilterQFunction<3>::RHS rhs_qf;
+         drh->AddDomainIntegrator(rhs_qf, rhsinputs, rhsoutput, ir,
+                                    domain_attributes);
+      }
+      else { MFEM_ABORT("Space dimension not supported"); }
+
+      Operator *Kop;
+      drh->FormSystemOperator(ess_tdofv, Kop);
+      Rh = std::make_unique<OperatorHandle>(Kop);
+      Rc = dynamic_cast<mfem::ConstrainedOperator*>(Kop);
+   }
+
+}
+
+void PDEFilter::Mult(const Vector &x, Vector &y) const
+{
+   input_field.SetFromTrueDofs(x);
+   dop->SetParameters({ nodes });
+   drh->SetParameters({ &input_field, nodes });
+   rhs.SetSize(ffes->GetTrueVSize());
+   drh->Mult(y, rhs);
+
+   {
+      ConstantCoefficient diffusion_coeff(filter_radius*filter_radius);
+      ParBilinearForm bf(ffes);
+      bf.AddDomainIntegrator(
+         new mfem::DiffusionIntegrator(diffusion_coeff));
+      bf.AddDomainIntegrator(new mfem::MassIntegrator());
+      bf.Assemble();
+      bf.Finalize();
+      std::unique_ptr<mfem::HypreParMatrix> Kb(bf.ParallelAssemble());
+      Kb->EliminateBC(ess_tdofv,Operator::DiagonalPolicy::DIAG_ONE);
+
+
+      Vector ost(ffes->GetTrueVSize()); ost.Randomize();
+      Vector tst1(ffes->GetTrueVSize());
+      Vector tst2(ffes->GetTrueVSize());
+      Kb->Mult(ost, tst1);
+      dop->Mult(ost,tst2);
+
+      real_t norm1 = 0.0;
+      real_t norm2 = 0.0;
+      norm1= InnerProduct(pmesh->GetComm(), tst1, tst1);
+      norm2= InnerProduct(pmesh->GetComm(), tst2, tst2);
+
+      if(0==pmesh->GetMyRank())
+      {   
+         std::cout<< "PDE Filter check: ||K*o|| = " << sqrt(norm1)
+               << ", ||Dop*o|| = " << sqrt(norm2) << std::endl;
+      }
+
+   }
+
+   // set BC to the solution vector y
+   // set BC in the rhs
+   // Kc->EliminateRHS(y, rhs);
+
+   ls->SetRelTol(linear_rtol);
+   ls->SetAbsTol(linear_atol);
+   ls->SetMaxIter(linear_iter);
+
+   ls->Mult(rhs, y);
+
+   filtered_field.SetFromTrueDofs(y);
+}
+
+void PDEFilter::MultTranspose(const Vector &x,
+                                        Vector &y) const
+{
+}
+
+
+DFEMLinElasticSolver::DFEMLinElasticSolver(ParMesh *mesh, int vorder)
+   : Operator(0, 0),
+     pmesh(mesh),
+     dim(mesh->Dimension()),
+     spaceDim(mesh->SpaceDimension()),
+     vfec(std::make_unique<H1_FECollection>(vorder, dim)),
+     // dfem path in upstream forces Ordering::byNODES
+     vfes(std::make_unique<ParFiniteElementSpace>(pmesh, vfec.get(), dim,
+                                                  Ordering::byNODES)),
+     fdisp(vfes.get())
+{
+   height = width = vfes->GetTrueVSize();
+   sol.SetSize(width);
+   rhs.SetSize(width);
+   sol = 0.0;
+   rhs = 0.0;
+
+   // default body force = 0
+   vol_force_vec.SetSize(spaceDim);
+   vol_force_vec = 0.0;
+
+}
+
+DFEMLinElasticSolver::~DFEMLinElasticSolver()
+{
+   // clean any owned BC coefficients
+   for (auto &kv : disp_bcs)
+   {
+      if (kv.second.owned && kv.second.coeff) { delete kv.second.coeff; }
+   }
+   disp_bcs.clear();
+
+   delete ls;
+   ls = nullptr;
+}
+
+void DFEMLinElasticSolver::SetLinearSolver(real_t rtol, real_t atol, int miter)
+{
+   linear_rtol = rtol;
+   linear_atol = atol;
+   linear_iter = miter;
+}
+
+
+void DFEMLinElasticSolver::SetVolForce(real_t fx, real_t fy, real_t fz)
+{
+   vol_force_vec.SetSize(spaceDim);
+   vol_force_vec = 0.0;
+   vol_force_vec[0] = fx;
+   vol_force_vec[1] = fy;
+   if (spaceDim == 3) { vol_force_vec[2] = fz; }
+
+   volforce_owned = std::make_unique<VectorConstantCoefficient>(vol_force_vec);
+   volforce = volforce_owned.get();
+}
+
+void DFEMLinElasticSolver::SetVolForce(VectorCoefficient &ff)
+{
+   volforce_owned.reset();
+   volforce = &ff;
+}
+
+void DFEMLinElasticSolver::AddSurfLoad(int id, real_t fx, real_t fy, real_t fz)
+{
+ 
+}
+
+void DFEMLinElasticSolver::AddSurfLoad(int id, VectorCoefficient &ff)
+{
+
+}
+
+void DFEMLinElasticSolver::AddDispBC(int id, int dir, real_t val)
+{
+   auto *c = new ConstantCoefficient(val);
+   disp_bcs.emplace(id, DispBC{dir, c, true});
+}
+
+void DFEMLinElasticSolver::AddDispBC(int id, int dir, Coefficient &val)
+{
+   disp_bcs.emplace(id, DispBC{dir, &val, false});
+}
+
+void DFEMLinElasticSolver::DelDispBC()
+{
+   for (auto &kv : disp_bcs)
+   {
+      if (kv.second.owned && kv.second.coeff) { delete kv.second.coeff; }
+   }
+   disp_bcs.clear();
+}
+
+void DFEMLinElasticSolver::SetEssTDofs(Vector &bsol, Array<int> &ess_dofs) const
+{
+  
+}
+
+void DFEMLinElasticSolver::Assemble()
+{
+
+
+}
+
+void DFEMLinElasticSolver::FSolve()
+{
+   MFEM_VERIFY(ls != nullptr && Kc != nullptr, "Call Assemble() before FSolve().");
+
+   ls->SetAbsTol(linear_atol);
+   ls->SetRelTol(linear_rtol);
+   ls->SetMaxIter(linear_iter);
+
+   // Eliminate RHS using constrained operator (dfem path)
+   Kc->EliminateRHS(sol, rhs);
+
+   ls->Mult(rhs, sol);
+}
+
+void DFEMLinElasticSolver::Mult(const Vector &x, Vector &y) const
+{
+   MFEM_VERIFY(ls != nullptr, "Call Assemble() before Mult().");
+   ls->Mult(x, y);
+
+   int N = ess_tdofv.Size();
+   real_t *yp = y.ReadWrite();
+   const real_t *sp = sol.Read();
+   const int *ep = ess_tdofv.Read();
+   mfem::forall(N, [=] MFEM_HOST_DEVICE(int i) { yp[ep[i]] = sp[ep[i]]; });
+}
+
+void DFEMLinElasticSolver::MultTranspose(const Vector &x, Vector &y) const
+{
+   ls->Mult(x, y);
+
+   int N = ess_tdofv.Size();
+   ess_tdofv.Read();
+
+   auto yp = y.Write();
+   const auto ep = ess_tdofv.Read();
+
+   mfem::forall(N, [=] MFEM_HOST_DEVICE(int i) { yp[ep[i]] = 0.0; });
+}
+
+

miniapps/mtop/mtop_solvers.hpp

@@ -326,3 +326,241 @@ private:
       }
    };
 };
+
+class PDEFilter : public mfem::Operator
+{
+public:
+   /// Construct the PDE filter for a given mesh and discretization
+   /// order.
+   PDEFilter(mfem::ParMesh *mesh, real_t r = 1.0, int order = 1);
+
+   /// Construct the PDE filter for a given input finite element space,
+   /// discretization order and filter radius.
+   PDEFilter(mfem::ParFiniteElementSpace *fespace, real_t r = 1.0, int order = 1); 
+   
+   /// Destructor of the filter.
+   virtual ~PDEFilter();
+
+   /// Set the filter radius
+   void SetFilterRadius(real_t r);
+
+   /// Set the linear solver relative tolerance (rtol),
+   /// absolute tolerance (atol) and maximum number of
+   /// iterations miter.
+   void SetLinearSolver(real_t rtol = 1e-8,
+                        real_t atol = 1e-12,
+                        int miter = 200);
+
+   /// Assemble the filter operator
+   void Assemble();
+
+   /// Apply the filter to input x, output y.
+   virtual
+   void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+   /// Adjoint apply the filter to input x, output y.
+   virtual
+   void MultTranspose(const mfem::Vector &x, mfem::Vector &y) const override;
+
+   /// Returns the filtered field.
+   mfem::ParGridFunction &GetFilteredField()
+   {
+      return filtered_field;
+   }
+
+   /// Returns the filtered finite element space
+   mfem::ParFiniteElementSpace &GetFilteredFESpace()
+   {
+      return *ffes;
+   }
+
+   /// Returns the input finite element space
+   mfem::ParFiniteElementSpace &GetInputFESpace()
+   {
+      return *ifes;
+   }
+
+   class NqptUniformParameterSpace : public
+      mfem::future::UniformParameterSpace
+   {
+   public:
+      NqptUniformParameterSpace(mfem::ParMesh &mesh,
+                                const mfem::IntegrationRule &ir,
+                                int vdim) :
+         mfem::future::UniformParameterSpace(mesh, ir, vdim, false)
+      {
+         dtq.nqpt = ir.GetNPoints();
+      }
+   };
+
+private:
+   mfem::ParMesh *pmesh;
+   const int dim; 
+   const int spaceDim;
+
+   real_t filter_radius;
+
+   mfem::FiniteElementCollection *ffec;
+   mfem::ParFiniteElementSpace *ffes;
+
+   mfem::FiniteElementCollection *ifec;
+   mfem::ParFiniteElementSpace *ifes;
+
+   mutable mfem::ParGridFunction filtered_field;
+   mutable mfem::ParGridFunction input_field;
+   mutable mfem::Vector rhs;
+
+   // H1 gradient of the input field: for more info see the paper
+   // "A Simple Introduciton to the SiMPL-method for density-based topology optimization"
+   // by D. Kim, B. Lazarov, T. Surowiec, B. Keith, 
+   // Structural and Multidisciplinary Optimization, 2025, 68.
+   mutable mfem::ParGridFunction h1_gradient;
+
+   // Linear solver parameters
+   real_t linear_rtol;
+   real_t linear_atol;
+   int linear_iter;
+
+   mfem::HypreBoomerAMG *prec; // preconditioner
+   mfem::CGSolver *ls;         // linear solver
+
+   // holds the displacement contrained DOFs
+   mfem::Array<int> ess_tdofv;
+
+   // boundary conditions
+   std::map<int, mfem::ConstantCoefficient> bcr;
+
+   // DFEM related definitions
+   std::unique_ptr<mfem::future::DifferentiableOperator> dop;
+   // RHS of the PDE filter
+   std::unique_ptr<mfem::future::DifferentiableOperator> drh;
+   const mfem::FiniteElement *fe;
+   mfem::ParGridFunction *nodes;
+   mfem::ParFiniteElementSpace *mfes;
+   const mfem::IntegrationRule &ir;
+   mfem::QuadratureSpace qs;
+   mfem::Array<int> domain_attributes;
+   NqptUniformParameterSpace diff_ps;
+   std::unique_ptr<mfem::CoefficientVector> diff_cv;
+   // Fsol - filtered solution, USol - unfiltered solution
+   // Coords - nodal coordinates
+   static constexpr int FSol = 0, USol=1, Coords = 2, DiffCoeff = 3;
+
+   mfem::HypreParMatrix *K;
+   mfem::ConstrainedOperator *Kc;
+   std::unique_ptr<mfem::OperatorHandle> Kh;
+   // RHS operators
+   mfem::ConstrainedOperator *Rc;
+   std::unique_ptr<mfem::OperatorHandle> Rh;
+};
+
+
+// A dFEM-only version of the mtop IsoLinElasticSolver:
+// - removes partial assembly option
+// - removes full/classical assembly option
+// - always builds the operator using mfem::future::DifferentiableOperator
+
+class DFEMLinElasticSolver : public mfem::Operator
+{
+public:
+   DFEMLinElasticSolver(mfem::ParMesh *mesh, int vorder = 1);
+   ~DFEMLinElasticSolver();
+
+   void SetLinearSolver(mfem::real_t rtol = 1e-8,
+                        mfem::real_t atol = 1e-12,
+                        int miter = 200);
+
+
+   // Volumetric force
+   void SetVolForce(mfem::real_t fx, mfem::real_t fy, mfem::real_t fz = 0.0);
+   void SetVolForce(mfem::VectorCoefficient &ff);
+
+   // Displacement BCs
+   void AddDispBC(int bdr_attr, int dir, mfem::real_t val);
+   void AddDispBC(int bdr_attr, int dir, mfem::Coefficient &val);
+   void DelDispBC();
+
+   // Surface loads (same interface pattern as mtop)
+   void AddSurfLoad(int bdr_attr, mfem::real_t fx, mfem::real_t fy, mfem::real_t fz = 0.0);
+   void AddSurfLoad(int bdr_attr, mfem::VectorCoefficient &ff);
+
+   // Build operator/preconditioner (call after setting material and BCs).
+   void Assemble();
+
+   // Solve Ku = f (forward)
+   void FSolve();
+
+   // Operator interface: y = K^{-1} x (solve with x as RHS)
+   void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+   void MultTranspose(const mfem::Vector &x, mfem::Vector &y) const override;
+
+   mfem::ParGridFunction &GetDisplacements(){
+         fdisp.SetFromTrueDofs(sol); return fdisp;
+   }
+
+   mfem::Vector          &GetSolutionVector() { return sol; }
+
+private:
+   // Essential dofs helper
+   void SetEssTDofs(mfem::Vector &bsol, mfem::Array<int> &ess_dofs) const;
+  
+
+   mfem::ParMesh *pmesh = nullptr;
+   int dim = 0;
+   int spaceDim = 0;
+
+   std::unique_ptr<mfem::FiniteElementCollection> vfec;
+   std::unique_ptr<mfem::ParFiniteElementSpace>   vfes;
+
+   // State solutions
+   mfem::ParGridFunction fdisp;
+   mfem::Vector sol, rhs;
+
+   // Body force
+   mfem::Vector vol_force_vec;
+   std::unique_ptr<mfem::VectorCoefficient> volforce_owned;
+   mfem::VectorCoefficient *volforce = nullptr;
+
+
+   // Dirichlet BC storage (attr -> per-direction coeff)
+   struct DispBC
+   {
+      int dir = -1; // 0,1,2 or -1(all)
+      mfem::Coefficient *coeff = nullptr; // may be owned or external
+      bool owned = false;
+   };
+
+   std::multimap<int, DispBC> disp_bcs;
+
+   mfem::Array<int> ess_tdofv;
+
+   // dFEM operator objects (future::DifferentiableOperator path)
+   static constexpr int U = 0, Coords = 1, LCoeff = 2, MuCoeff = 3;
+
+   // quadrature / coefficient sampling
+   mfem::QuadratureSpace *qs = nullptr;
+   std::unique_ptr<mfem::CoefficientVector> Lambda_cv;
+   std::unique_ptr<mfem::CoefficientVector> Mu_cv;
+
+   // The differentiable operator itself
+   std::unique_ptr<mfem::future::DifferentiableOperator> dop;
+
+   // Constrained operator handle (for eliminate RHS)
+   std::unique_ptr<mfem::OperatorHandle> Kh;
+   mfem::ConstrainedOperator *Kc = nullptr;
+
+   // Linear solver + preconditioner
+   mfem::CGSolver *ls = nullptr;
+
+   // LOR preconditioner 
+   std::unique_ptr<mfem::ParLORDiscretization> lor_disc;
+   std::unique_ptr<mfem::ParFiniteElementSpace> lor_scalar_fespace;
+
+   std::unique_ptr<mfem::HypreParMatrix>  lor_mat;
+   std::unique_ptr<mfem::HypreBoomerAMG>  lor_amg;
+
+   // Solver params
+   mfem::real_t linear_rtol = 1e-8;
+   mfem::real_t linear_atol = 1e-12;
+   int linear_iter = 200;
+};

miniapps/mtop/mtop_test_dfem_elast.cpp

@@ -0,0 +1,155 @@
+// 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.
+//
+// Sample runs:
+//    mpirun -np 4 mtop_test_iso_elasticity -tri -o 2
+//
+//    mpirun -np 4 mtop_test_iso_elasticity -quad -o 2
+//
+// Device sample runs:
+//    mpirun -np 4 mtop_test_iso_elasticity -d gpu -quad -o 2
+
+#include "mtop_solvers.hpp"
+
+using namespace std;
+using namespace mfem;
+
+constexpr auto MESH_TRI = MFEM_SOURCE_DIR "/miniapps/mtop/sq_2D_9_tri.mesh";
+constexpr auto MESH_QUAD = MFEM_SOURCE_DIR "/miniapps/mtop/sq_2D_9_quad.mesh";
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init();
+   Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = MESH_QUAD;
+   const char *device_config = "cpu";
+   int order = 2;
+   bool pa = false;
+   bool dfem = false;
+   bool mesh_tri = false;
+   bool mesh_quad = false;
+   int par_ref_levels = 1;
+   bool paraview = false;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&dfem, "-dfem", "--dFEM", "-no-dfem", "--no-dFEM",
+                  "Enable or not dFEM.");
+   args.AddOption(&mesh_tri, "-tri", "--triangular", "-no-tri",
+                  "--no-triangular", "Enable or not triangular mesh.");
+   args.AddOption(&mesh_quad, "-quad", "--quadrilateral", "-no-quad",
+                  "--no-quadrilateral", "Enable or not quadrilateral mesh.");
+   args.AddOption(&par_ref_levels, "-prl", "--par-ref-levels",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&paraview, "-pv", "--paraview", "-no-pv", "--no-paraview",
+                  "Enable or not Paraview visualization");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+   MFEM_VERIFY(!(pa && dfem), "pa and dfem cannot be both set");
+
+   // Enable hardware devices such as GPUs, and programming models such as
+   // CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   // Read the (serial) mesh from the given mesh file on all processors.  We
+   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   // and volume meshes with the same code.
+   Mesh mesh(mesh_tri ? MESH_TRI : mesh_quad ? MESH_QUAD : mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+
+   // Refine the serial mesh on all processors 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 1000 elements.
+   {
+      const int ref_levels =
+         (int)floor(log(1000. / mesh.GetNE()) / log(2.) / dim);
+      for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
+   }
+   if (Mpi::Root())
+   {
+      std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
+   }
+
+   // Define a parallel mesh by a partitioning of the serial mesh. Refine
+   // this mesh further in parallel to increase the resolution. Once the
+   // parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   for (int l = 0; l < par_ref_levels; l++) { pmesh.UniformRefinement(); }
+
+   // Create the solver
+   DFEMLinElasticSolver elsolver(&pmesh, order);
+   if (Mpi::Root())
+   {
+      std::cout << "Number of unknowns: "
+                << elsolver.GetSolutionVector().Size() << std::endl;
+   }
+
+   // set boundary conditions
+   elsolver.AddDispBC(2, -1, 0.0);
+   elsolver.AddDispBC(5, -1, 0.0);
+
+   // set material properties
+   ConstantCoefficient E(1.0), nu(0.2);
+   //elsolver.SetMaterial(E, nu);
+
+   // set surface load
+   //elsolver.AddSurfLoad(1, 0.0, 1.0);
+
+   // set convergence tolerances and max iterations
+   elsolver.SetLinearSolver(1e-6,1e-8,100);
+
+   // assemble the discrete system
+   elsolver.Assemble();
+
+   // solve the system
+   elsolver.FSolve();
+
+   // extract the solution
+   ParGridFunction &sol = elsolver.GetDisplacements();
+
+   if (paraview)
+   {
+      ParaViewDataCollection paraview_dc("isoel", &pmesh);
+      paraview_dc.SetPrefixPath("ParaView");
+      paraview_dc.SetLevelsOfDetail(order);
+      paraview_dc.SetDataFormat(VTKFormat::BINARY);
+      paraview_dc.SetHighOrderOutput(true);
+      paraview_dc.SetCycle(0);
+      paraview_dc.SetTime(0.0);
+      paraview_dc.RegisterField("disp", &sol);
+      paraview_dc.Save();
+   }
+
+   if (socketstream glvis; visualization &&
+       (glvis.open("localhost", 19916), glvis.is_open()))
+   {
+      glvis << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n";
+      glvis << "solution\n" << pmesh << sol << std::flush;
+   }
+
+   return EXIT_SUCCESS;
+}

miniapps/mtop/mtop_test_filt.cpp

@@ -0,0 +1,159 @@
+// 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.
+//
+// Sample runs:
+//    mpirun -np 4 mtop_test_iso_elasticity -tri -o 2
+//
+//    mpirun -np 4 mtop_test_iso_elasticity -quad -o 2
+//
+// Device sample runs:
+//    mpirun -np 4 mtop_test_iso_elasticity -d gpu -quad -o 2
+
+#include "mtop_solvers.hpp"
+
+using namespace std;
+using namespace mfem;
+
+constexpr auto MESH_TRI = MFEM_SOURCE_DIR "/miniapps/mtop/sq_2D_9_tri.mesh";
+constexpr auto MESH_QUAD = MFEM_SOURCE_DIR "/miniapps/mtop/sq_2D_9_quad.mesh";
+
+class DensCoeff : public mfem::Coefficient
+{
+private:
+    real_t l;
+public:
+    DensCoeff(real_t d=1.0) : l(d) {}
+
+    virtual real_t Eval(mfem::ElementTransformation &T,
+                        const mfem::IntegrationPoint &ip)
+    {
+        Vector x;
+        T.Transform(ip, x);
+        real_t r = x.Norml2();
+        r=sin(M_PI*r/l);
+        if(r>0.5)
+            r=1.0;
+        else
+            r=0.0;
+        return r;
+    }
+};
+
+
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init();
+   Hypre::Init();
+
+    // Parse command-line options.
+    const char *mesh_file = MESH_QUAD;
+    const char *device_config = "cpu";
+    int order = 2;
+    bool mesh_tri = false;
+    bool mesh_quad = false;
+    int par_ref_levels = 1;
+    bool paraview = false;
+    bool visualization = true;
+    OptionsParser args(argc, argv);
+    args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+    args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+    args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+    args.AddOption(&mesh_tri, "-tri", "--triangular", "-no-tri",
+                  "--no-triangular", "Enable or not triangular mesh.");
+    args.AddOption(&mesh_quad, "-quad", "--quadrilateral", "-no-quad",
+                  "--no-quadrilateral", "Enable or not quadrilateral mesh.");
+    args.AddOption(&par_ref_levels, "-prl", "--par-ref-levels",
+                  "Number of times to refine the mesh uniformly in parallel.");
+    args.AddOption(&paraview, "-pv", "--paraview", "-no-pv", "--no-paraview",
+                  "Enable or not Paraview visualization");
+    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+    args.ParseCheck();
+
+    // Enable hardware devices such as GPUs, and programming models such as
+    // CUDA, OCCA, RAJA and OpenMP based on command line options.
+    Device device(device_config);
+    if (Mpi::Root()) { device.Print(); }
+
+    // Read the (serial) mesh from the given mesh file on all processors.  We
+    // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+    // and volume meshes with the same code.
+    Mesh mesh(mesh_tri ? MESH_TRI : mesh_quad ? MESH_QUAD : mesh_file, 1, 1);
+    const int dim = mesh.Dimension();
+
+    // Refine the serial mesh on all processors 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 1000 elements.
+    {
+       const int ref_levels =
+          (int)floor(log(1000. / mesh.GetNE()) / log(2.) / dim);
+       for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
+    }
+    if (Mpi::Root())
+    {
+       std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
+    }
+
+    // Define a parallel mesh by a partitioning of the serial mesh. Refine
+    // this mesh further in parallel to increase the resolution. Once the
+    // parallel mesh is defined, the serial mesh can be deleted.
+    ParMesh pmesh(MPI_COMM_WORLD, mesh);
+    mesh.Clear();
+    for (int l = 0; l < par_ref_levels; l++) { pmesh.UniformRefinement(); }
+
+
+    PDEFilter* filt=new PDEFilter(&pmesh, 0.1, order);
+
+    filt->Assemble();
+
+    // define two grid functions: one for the input and one for the output
+    ParGridFunction filt_gf(&filt->GetFilteredFESpace());
+    ParGridFunction orig_gf(&filt->GetInputFESpace());
+    QuadratureFunction orig_qf;
+                    
+
+    filt_gf.GetTrueVector()=0.0;
+
+    DensCoeff dens_coeff(1.0);
+    orig_gf.ProjectCoefficient(dens_coeff);
+
+    filt->Mult(orig_gf.GetTrueVector(), filt_gf.GetTrueVector());
+    filt_gf.SetFromTrueVector();
+
+
+
+    if (paraview)
+    {
+        ParaViewDataCollection paraview_dc("isoel", &pmesh);
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetLevelsOfDetail(order);
+        paraview_dc.SetDataFormat(VTKFormat::BINARY);
+        paraview_dc.SetHighOrderOutput(true);
+        paraview_dc.SetCycle(0);
+        paraview_dc.SetTime(0.0);
+        paraview_dc.RegisterField("filt", &filt_gf);
+        paraview_dc.RegisterField("orig", &orig_gf);
+        paraview_dc.Save();
+    }
+
+
+    delete filt;
+    return EXIT_SUCCESS;
+}
+

miniapps/mtop/test_dfem.cpp

@@ -0,0 +1,641 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#include "mtop_solvers.hpp"
+
+using namespace mfem;
+
+using mfem::future::dual;
+using mfem::future::tuple;
+using mfem::future::tensor;
+
+using mfem::future::Weight;
+using mfem::future::Gradient;
+using mfem::future::Identity;
+
+
+
+///////////////////////////////////////////////////////////////////////////////
+/// \brief The SQFunction struct defining the Stokes operator at
+/// integration points which is valid in 2D and 3D
+template <int DIM, typename scalar_t=real_t> struct SQFunction
+{
+   using mati_t = tensor<scalar_t, DIM, DIM>;
+
+   struct Stokes
+   {
+      MFEM_HOST_DEVICE inline auto operator()(const mati_t &dudxi,
+                                              const real_t &M, // viscosity
+                                              const mati_t &J,
+                                              const real_t &w) const
+      {
+         /*
+         mati_t invJ = mfem::future::inv<scalar_t>(J);
+         const mati_t JxW = transpose(invJ) * det(J) * w;
+         const auto eps = sym(dudxi * invJ);
+         return tuple{(2.0 * M * eps) * JxW};
+         */
+
+         const mati_t invJ=mfem::future::inv<scalar_t>(J);
+         const mati_t JxW = transpose(invJ) * det(J) * w;
+         const auto eps = mfem::future::sym(dudxi * invJ);
+         return tuple{(2.0 * M * eps) * JxW};
+
+      }
+   };
+};
+
+///////////////////////////////////////////////////////////////////////////////
+/// \brief The Mass MQFunction struct defining Vector Mass operator
+/// at integration points which is valid in 2D and 3D
+template <int DIM,typename scalar_t=real_t> struct MQFunction
+{
+   using mati_t = tensor<scalar_t, DIM, DIM>;
+   using veci_t = tensor<scalar_t, DIM>;
+   struct Mass
+   {
+      MFEM_HOST_DEVICE inline auto operator()(const veci_t &u,
+                                              const scalar_t &M, // mass coefficient
+                                              const mati_t &J,
+                                              const real_t &w) const
+      {
+         return tuple{(M * u) * det(J) * w};
+      }
+   };
+};
+
+
+///////////////////////////////////////////////////////////////////////////////
+/// \brief The Anisotropic AEQFunction struct defining the Elasticity operator
+/// at integration points which is valid in 3D
+template <int NMAT,typename scalar_t=real_t> struct AEQFunction3D
+{
+   // Dimension
+   static constexpr int DIM = 3;
+
+   // Number of independent components in a symmetric DIM×DIM matrix
+   static constexpr int NVOIGT = DIM * (DIM + 1) / 2;
+
+   // Number of entries in the material matrix
+   static constexpr int NMAT_ENTRIES = NVOIGT * (NVOIGT+1) / 2;
+
+   // Total size of the flat array
+   static constexpr int SIZE = NMAT * NMAT_ENTRIES;
+
+   // The actual packed data: [NMAT][NVOIGT] in row-major order
+   real_t data[SIZE];
+
+   // ------------------------------------------------------------------
+   // Mapping (i,j) → Voigt index for order: 11,22,33, 23,13,12
+   // ------------------------------------------------------------------
+   //  (0,0)→0   (1,1)→1   (2,2)→2
+   //  (1,2)→3   (0,2)→4   (0,1)→5   (and symmetric)
+   MFEM_HOST_DEVICE inline
+   static constexpr int voigt_index(int i, int j)
+   {
+      if (i > j) { int tmp = i; i = j; j = tmp; }  // ensure i <= j
+
+      if (i == 0 && j == 0) { return 0; }
+      if (i == 1 && j == 1) { return 1; }
+      if (i == 2 && j == 2) { return 2; }
+      if (i == 1 && j == 2) { return 3; }
+      if (i == 0 && j == 2) { return 4; }
+      if (i == 0 && j == 1) { return 5; }
+
+      return -1; // unreachable
+   }
+
+   // ------------------------------------------------------------------
+   // Accessors
+   // ------------------------------------------------------------------
+   MFEM_HOST_DEVICE inline
+   static constexpr int mat_index(int i, int j, int k, int l)
+   {
+      int ij = voigt_index(i,j);
+      int kl = voigt_index(k,l);
+      if (ij > kl) { int tmp = ij; ij = kl; kl = tmp; }  // ensure ij <= kl
+      return (kl * (kl + 1)) / 2 + ij;
+   }
+
+   MFEM_HOST_DEVICE inline
+   static constexpr int mat_index(int i, int j)
+   {
+      if (i>j) { int tmp = i; i = j; j = tmp; }
+      {
+         return (j * (j + 1)) / 2 + i;
+      }
+   }
+
+   using mati_t = tensor<scalar_t, DIM, DIM>;
+   using veci_t = tensor<scalar_t, NMAT>;
+
+   struct Elasticity
+   {
+      MFEM_HOST_DEVICE inline auto operator()(const mati_t &dudxi,
+                                              const veci_t &rhoi,
+                                              const mati_t &J,
+                                              const real_t &w) const
+      {
+         const mati_t invJ = mfem::future::inv<scalar_t>(J);
+         const mati_t JxW = transpose(invJ) * det(J) * w;
+
+         const auto eps = mfem::future::sym(dudxi * invJ);
+         auto str= 0.0 * eps; // initialize to zero the stress tensor
+         for (int im=0; im<NMAT; im++)
+         {
+            for (int i=0; i<DIM; i++)
+            {
+               for (int j=0; j<DIM; j++)
+               {
+                  for (int k=0; k<DIM; k++)
+                  {
+                     for (int l=0; l<DIM; l++)
+                     {
+                        // Voigt notation mapping
+                        str[i][j] += rhoi[im] * data[im*NMAT_ENTRIES + mat_index(i,j,k,l)] * eps[k][l];
+                     }
+                  }
+               }
+            }
+         }
+
+         return tuple{str*JxW};
+      }
+   };
+};
+
+
+
+
+
+class NqptUniformParameterSpace : public
+   mfem::future::UniformParameterSpace
+{
+public:
+   NqptUniformParameterSpace(mfem::ParMesh &mesh,
+                             const mfem::IntegrationRule &ir,
+                             int vdim) :
+      mfem::future::UniformParameterSpace(mesh, ir, vdim, false)
+   {
+      dtq.nqpt = ir.GetNPoints();
+   }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   mfem::Mpi::Init(argc, argv);
+   int myrank = mfem::Mpi::WorldRank();
+   mfem::Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = "./dfg_bench_flow_tri.msh";
+   int order = 2;
+   bool static_cond = false;
+   int ser_ref_levels = 1;
+   int par_ref_levels = 1;
+   real_t newton_rel_tol = 1e-7;
+   real_t newton_abs_tol = 1e-12;
+   int newton_iter = 10;
+   int print_level = 1;
+   bool visualization = false;
+
+   mfem::OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels,
+                  "-rs",
+                  "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels,
+                  "-rp",
+                  "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&visualization,
+                  "-vis",
+                  "--visualization",
+                  "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&newton_rel_tol,
+                  "-rel",
+                  "--relative-tolerance",
+                  "Relative tolerance for the Newton solve.");
+   args.AddOption(&newton_abs_tol,
+                  "-abs",
+                  "--absolute-tolerance",
+                  "Absolute tolerance for the Newton solve.");
+   args.AddOption(&newton_iter,
+                  "-it",
+                  "--newton-iterations",
+                  "Maximum iterations for the Newton solve.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(std::cout);
+      }
+      return 1;
+   }
+
+   if (myrank == 0)
+   {
+      args.PrintOptions(std::cout);
+   }
+
+   // Read the (serial) mesh from the given mesh file on all processors.  We
+   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   // and volume meshes with the same code.
+   mfem::Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int spaceDim = mesh.SpaceDimension();
+
+
+   // Refine the serial mesh on all processors 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 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(1000./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // Define a parallel mesh by a partitioning of the serial mesh. Refine
+   // this mesh further in parallel to increase the resolution. Once the
+   // parallel mesh is defined, the serial mesh can be deleted.
+   mfem::ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+
+   std::cout<<"My rank="<<pmesh.GetMyRank()<<std::endl;
+   pmesh.PrintInfo(std::cout);
+
+   H1_FECollection* vfec=new H1_FECollection(order, dim);
+   H1_FECollection* pfec=new H1_FECollection(order-1, dim);
+
+
+   //construct the FEM spaces
+   ParFiniteElementSpace* vfes= new ParFiniteElementSpace(&pmesh, vfec, dim,
+                                                          Ordering::byNODES);
+   ParFiniteElementSpace* pfes=new ParFiniteElementSpace(&pmesh, pfec);
+
+   Vector y; y.SetSize(vfes->TrueVSize()); y.Randomize();
+   Vector x; x.SetSize(vfes->TrueVSize()); x.Randomize();
+   ParGridFunction xgf(vfes);
+   vfes->GetProlongationMatrix()->Mult(x, xgf);
+
+
+   ConstantCoefficient zerocoef(0.0);
+   ConstantCoefficient visc(1.0);
+
+   Array<int> ess_tdofv;
+
+   ParBilinearForm* stokes_bf=new ParBilinearForm(vfes);
+   stokes_bf->AddDomainIntegrator(new ElasticityIntegrator(zerocoef,visc));
+   stokes_bf->Assemble();
+   stokes_bf->Finalize();
+   std::unique_ptr<HypreParMatrix> A(stokes_bf->ParallelAssemble());
+   delete stokes_bf;
+
+   //partial assembly
+   stokes_bf=new ParBilinearForm(vfes);
+   stokes_bf->AddDomainIntegrator(new ElasticityIntegrator(zerocoef,visc));
+   stokes_bf->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+   stokes_bf->Assemble();
+   mfem::ConstrainedOperator *Kc;
+   std::unique_ptr<mfem::OperatorHandle> Kh;
+   {
+      Operator *Kop;
+      stokes_bf->FormSystemOperator(ess_tdofv, Kop);
+      Kh = std::make_unique<OperatorHandle>(Kop);
+      Kc = dynamic_cast<mfem::ConstrainedOperator*>(Kop);
+   }
+
+   //dfem operator
+   static constexpr int U = 0, Coords = 1, LCoeff = 2, MuCoeff = 3;
+   std::unique_ptr<mfem::CoefficientVector> Mu_cv;
+   std::unique_ptr<mfem::future::DifferentiableOperator> dop;
+
+   pmesh.EnsureNodes();
+   ParGridFunction* nodes(static_cast<ParGridFunction *>(pmesh.GetNodes()));
+   ParFiniteElementSpace* mfes(static_cast<ParFiniteElementSpace*>
+                               (nodes->ParFESpace()));
+
+   const mfem::FiniteElement *fe(vfes->GetFE(0));
+   //const mfem::IntegrationRule &ir(IntRules.Get(fe->GetGeomType(),
+   //                fe->GetOrder() + fe->GetOrder() + fe->GetDim() - 1 + 5));
+   const mfem::IntegrationRule &ir(IntRules.Get(fe->GetGeomType(),
+                                                fe->GetOrder() + fe->GetOrder()));
+   mfem::QuadratureSpace qs(pmesh, ir);
+   NqptUniformParameterSpace Mu_ps(pmesh, ir, 1) ;
+
+   std::cout<<"rank="<<myrank<<" qspace.size="<<qs.GetSize()<<std::endl;
+
+   std::cout<<"rank="<<myrank<<" mu_ps.size="<<Mu_ps.GetTrueVSize()<<std::endl;
+
+   Array<int> domain_attributes;
+   if (pmesh.attributes.Size() > 0)
+   {
+      domain_attributes.SetSize(pmesh.attributes.Max());
+      domain_attributes = 1;
+   }
+
+   // sample mu on the integration points
+   Mu_cv = std::make_unique<CoefficientVector>(visc, qs);
+
+   std::cout<<"rank="<<myrank<<" m_cv.size="<<Mu_cv->Size()<<std::endl;
+
+
+   // define the differentiable operator
+   dop = std::make_unique<mfem::future::DifferentiableOperator>(
+   std::vector<mfem::future::FieldDescriptor> {{ U, vfes }},
+   std::vector<mfem::future::FieldDescriptor>
+   {
+      { MuCoeff, &Mu_ps},
+      { Coords, mfes }
+   },
+   pmesh);
+
+   dop->SetParameters({ Mu_cv.get(), nodes });
+
+   const auto inputs =
+      mfem::future::tuple{ Gradient<U>{},
+                           Identity<MuCoeff>{},
+                           Gradient<Coords>{},
+                           Weight{} };
+
+   const auto output = mfem::future::tuple{ Gradient<U>{} };
+
+   //const auto output = mfem::future::tuple{ Identity<MuCoeff>{} };
+
+   //define the q-function
+   if (2 == spaceDim)
+   {
+      typename SQFunction<2>::Stokes s2qf;
+      dop->AddDomainIntegrator(s2qf, inputs, output, ir, domain_attributes);
+   }
+   else if (3 == spaceDim)
+   {
+      typename SQFunction<3>::Stokes s3qf;
+      dop->AddDomainIntegrator(s3qf, inputs, output, ir, domain_attributes);
+   }
+   else { MFEM_ABORT("Space dimension not supported"); }
+
+   mfem::ConstrainedOperator *Kdc;
+   std::unique_ptr<mfem::OperatorHandle> Kdh;
+   {
+      Operator *Kdop;
+      dop->FormSystemOperator(ess_tdofv, Kdop);
+      Kdh = std::make_unique<OperatorHandle>(Kdop);
+      Kdc = dynamic_cast<mfem::ConstrainedOperator*>(Kdop);
+   }
+
+
+   //differentiable operator
+   std::unique_ptr<mfem::future::DifferentiableOperator> dopd;
+   //define the differentiable operator
+   dopd = std::make_unique<mfem::future::DifferentiableOperator>(
+   std::vector<mfem::future::FieldDescriptor> {{ U, vfes }},
+   std::vector<mfem::future::FieldDescriptor>
+   {
+      { MuCoeff, &Mu_ps},
+      { Coords, mfes }
+   },
+   pmesh);
+
+   //dopd->SetParameters({ Mu_cv.get(), nodes });
+   //define the q-function
+   if (2 == spaceDim)
+   {
+      using mfem::future::dual;
+      using dual_t = dual<real_t, real_t>;
+      typename SQFunction<2,dual_t>::Stokes s2qf;
+      auto derivatives = std::integer_sequence<size_t, U, Coords> {};
+      dopd->AddDomainIntegrator(s2qf, inputs, output, ir, domain_attributes,
+                                derivatives);
+
+   }
+   else if (3 == spaceDim)
+   {
+      using mfem::future::dual;
+      using dual_t = dual<real_t, real_t>;
+      typename SQFunction<3,dual_t>::Stokes s3qf;
+      auto derivatives = std::integer_sequence<size_t, U, Coords> {};
+      dopd->AddDomainIntegrator(s3qf, inputs, output, ir, domain_attributes,
+                                derivatives);
+   }
+   else { MFEM_ABORT("Space dimension not supported"); }
+
+
+
+   std::shared_ptr<mfem::future::DerivativeOperator> dres_du;
+   //the parameters should be set from grid functions
+   dres_du=dopd->GetDerivative(U, {&xgf}, {Mu_cv.get(), nodes});
+
+
+   std::shared_ptr<mfem::future::DerivativeOperator> dres_dcoor;
+   dres_dcoor=dopd->GetDerivative(Coords, {&xgf}, {Mu_cv.get(), nodes});
+
+
+   double dt;
+   int maxnit=2;
+   //test the full assembly matrix-vector product
+   tic();
+   for (int i=0; i<maxnit; i++) {  A->Mult(x, y);   }
+   dt=toc();
+   real_t normy1=InnerProduct(pmesh.GetComm(), y, y);
+   if (0==myrank)
+   {
+      std::cout << "y = Af*x computed: " << normy1 << " dt=" << dt <<std::endl;
+   }
+
+   //test partial matrix-vector product
+   tic();
+   for (int i=0; i<maxnit; i++) {   Kc->Mult(x, y);}
+   dt=toc();
+   real_t normy2=InnerProduct(pmesh.GetComm(), y, y);
+   if (0==myrank)
+   {
+      std::cout << "y = Ap*x computed: " << normy2 << " dt=" << dt <<std::endl;
+      //  std::cout << "Difference norm = " << fabs(normy1 - normy2) <<std::endl;
+   }
+
+   //test dfem matrix-vector product
+   tic();
+   for (int i=0; i<maxnit; i++) {   Kdc->Mult(x, y); }
+   dt=toc();
+   real_t normy3=InnerProduct(pmesh.GetComm(), y, y);
+   if (0==myrank)
+   {
+      std::cout << "y = Ad*x computed: " << normy3 << " dt=" << dt <<std::endl;
+      // std::cout << "Difference norm = " << fabs(normy1 - normy3) <<std::endl;
+   }
+
+   //test dfem derivative operator matrix-vector product
+   tic();
+   for (int i=0; i<maxnit; i++) {   dres_du->Mult(x, y);}
+   dt=toc();
+   real_t normy4=InnerProduct(pmesh.GetComm(), y, y);
+   if (0==myrank)
+   {
+      std::cout << "y = Af*x computed: " << normy4 << " dt=" << dt <<std::endl;
+      //  std::cout << "Difference norm = " << fabs(normy1 - normy4) <<std::endl;
+   }
+
+   Vector resc; resc.SetSize(mfes->TrueVSize());
+   Vector inpv; inpv.SetSize(mfes->TrueVSize()); inpv.Randomize();
+   dres_dcoor->Mult(inpv, resc);
+   real_t normc=InnerProduct(pmesh.GetComm(), resc, resc);
+   real_t normi=InnerProduct(pmesh.GetComm(), inpv, inpv);
+   if (0==myrank)
+   {
+      std::cout << "y = Adcoor*x computed: " << normc <<" inp=" <<normi <<std::endl;
+   }
+
+
+   delete stokes_bf;
+
+   //vector mass tests
+   ConstantCoefficient mass_coef(1.0);
+   const auto mass_inputs =
+      mfem::future::tuple{ mfem::future::Value<U>{},
+                           mfem::future::Identity<MuCoeff>{},
+                           mfem::future::Gradient<Coords>{},
+                           mfem::future::Weight{} };
+
+   const auto mass_output = mfem::future::tuple{ mfem::future::Value<U>{} };
+   // sample \rho on the integration points
+   Mu_cv = std::make_unique<CoefficientVector>(mass_coef, qs);
+   //differentiable operator
+   std::unique_ptr<mfem::future::DifferentiableOperator> mopd;
+   //define the differentiable operator
+   mopd = std::make_unique<mfem::future::DifferentiableOperator>(
+   std::vector<mfem::future::FieldDescriptor> {{ U, vfes }},
+   std::vector<mfem::future::FieldDescriptor>
+   {
+      { MuCoeff, &Mu_ps}, //same dimmension as viscoity for Stokes
+      { Coords, mfes }
+   },
+   pmesh);
+
+   //define the q-function
+   mopd->SetParameters({ Mu_cv.get(), nodes });
+   if (2 == spaceDim)
+   {
+      using mfem::future::dual;
+      using dual_t = dual<real_t, real_t>;
+      typename MQFunction<2,dual_t>::Mass m2qf;
+      auto derivatives = std::integer_sequence<size_t, U, Coords> {};
+      mopd->AddDomainIntegrator(m2qf, mass_inputs, mass_output, ir, domain_attributes,
+                                derivatives);
+
+   }
+   else if (3 == spaceDim)
+   {
+      using mfem::future::dual;
+      using dual_t = dual<real_t, real_t>;
+      typename MQFunction<3,dual_t>::Mass m3qf;
+      auto derivatives = std::integer_sequence<size_t, U, Coords> {};
+      mopd->AddDomainIntegrator(m3qf, mass_inputs, mass_output, ir, domain_attributes,
+                                derivatives);
+   }
+   else { MFEM_ABORT("Space dimension not supported"); }
+
+   std::shared_ptr<mfem::future::DerivativeOperator> dmass_du;
+   //the parameters should be set from grid functions
+   dmass_du=mopd->GetDerivative(U, {&xgf}, {Mu_cv.get(), nodes});
+   std::shared_ptr<mfem::future::DerivativeOperator> dmass_dcoor;
+   dmass_dcoor=mopd->GetDerivative(Coords, {&xgf}, {Mu_cv.get(), nodes});
+
+   //test dfem derivative operator matrix-vector product
+   tic();
+   for (int i=0; i<maxnit; i++) {   dmass_du->Mult(x, y); }
+   dt=toc();
+   real_t normm1=InnerProduct(pmesh.GetComm(), y, y);
+   if (0==myrank)
+   {
+      std::cout << "y = Mf*x computed: " << normm1 << " dt=" << dt <<std::endl;
+   }
+
+   mfem::ConstrainedOperator *Mdc;
+   std::unique_ptr<mfem::OperatorHandle> Mdh;
+   {
+      Operator *Mdop;
+      mopd->FormSystemOperator(ess_tdofv, Mdop);
+      Mdh = std::make_unique<OperatorHandle>(Mdop);
+      Mdc = dynamic_cast<mfem::ConstrainedOperator*>(Mdop);
+   }
+
+   //test dfem derivative operator matrix-vector product
+   tic();
+   for (int i=0; i<maxnit; i++) {   Mdc->Mult(x, y); }
+   dt=toc();
+   normm1=InnerProduct(pmesh.GetComm(), y, y);
+   if (0==myrank)
+   {
+      std::cout << "y = Mf*x computed: " << normm1 << " dt=" << dt <<std::endl;
+   }
+
+   stokes_bf=new ParBilinearForm(vfes);
+   stokes_bf->AddDomainIntegrator(new VectorMassIntegrator(mass_coef));
+   stokes_bf->Assemble();
+   {
+      Operator *Kop;
+      stokes_bf->FormSystemOperator(ess_tdofv, Kop);
+      Kh = std::make_unique<OperatorHandle>(Kop);
+      Kc = dynamic_cast<mfem::ConstrainedOperator*>(Kop);
+   }
+
+   //test partial matrix-vector product
+   tic();
+   for (int i=0; i<maxnit; i++) {   Kc->Mult(x, y);}
+   dt=toc();
+   real_t normm2=InnerProduct(pmesh.GetComm(), y, y);
+   if (0==myrank)
+   {
+      std::cout << "y = Mp*x computed: " << normm2 << " dt=" << dt <<std::endl;
+      std::cout << "Difference norm = " << fabs(normm1 - normm2) <<std::endl;
+   }
+
+   mfem::PowerMethod pm(pmesh.GetComm());
+   Vector v0(x.Size()); v0.Randomize();
+   double eigenvalue=pm.EstimateLargestEigenvalue(*Kdc, v0, 20, 1e-12);
+   if (0==myrank)
+   {
+      std::cout << "Largest eigenvalue of dfem mass matrix: " << eigenvalue
+                <<std::endl;
+   }
+
+
+
+
+
+
+   delete stokes_bf;
+
+
+   delete pfes;
+   delete vfes;
+   delete pfec;
+   delete vfec;
+
+   MPI::Finalize();
+   return 0;
+}