documentation

.gitignore

@@ -63,7 +63,6 @@ examples/ex3[0-9]
 examples/ex3[0-9]p
 examples/ex4[0-9]
 examples/ex4[0-9]p
-examples/vector-dg-diffusion
 
 examples/refined.mesh
 examples/displaced.mesh

examples/vector-dg-diffusion.cpp

@@ -1,164 +0,0 @@
-#include "mfem.hpp"
-#include "vector-dg-diffusion.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-class RepeatedCoefficient : public VectorCoefficient
-{
-   Coefficient &coeff;
-public:
-   RepeatedCoefficient(int dim, Coefficient &coeff_)
-      : VectorCoefficient(dim), coeff(coeff_)
-   { }
-   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
-   {
-      V.SetSize(vdim);
-      V = coeff.Eval(T, ip);
-   }
-};
-
-real_t u_fn(const Vector &xvec);
-real_t f_fn(const Vector &xvec);
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int ref_levels = 0;
-   int order = 1;
-   real_t sigma = -1.0;
-   real_t kappa = -1.0;
-   const char *device_config = "cpu";
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ref_levels, "-r", "--refine",
-                  "Number of times to refine the mesh uniformly, -1 for auto.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) >= 0.");
-   args.AddOption(&sigma, "-s", "--sigma",
-                  "One of the three DG penalty parameters, typically +1/-1."
-                  " See the documentation of class DGDiffusionIntegrator.");
-   args.AddOption(&kappa, "-k", "--kappa",
-                  "One of the three DG penalty parameters, should be positive."
-                  " Negative values are replaced with (order+1)^2.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.ParseCheck();
-
-   if (kappa < 0)
-   {
-      kappa = (order+1)*(order+1);
-   }
-
-   Device device(device_config);
-   device.Print();
-
-   Mesh mesh(mesh_file);
-   const int dim = mesh.Dimension();
-
-   {
-      if (ref_levels < 0)
-      {
-         ref_levels = (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
-      }
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   DG_FECollection fec(order, dim);
-   FiniteElementSpace fespace(&mesh, &fec, dim);
-   cout << "Number of unknowns: " << fespace.GetVSize() << endl;
-
-
-   FunctionCoefficient scalar_f_coeff(f_fn);
-   FunctionCoefficient scalar_u_coeff(u_fn);
-   RepeatedCoefficient f_coeff(dim, scalar_f_coeff);
-   RepeatedCoefficient u_coeff(dim, scalar_u_coeff);
-
-   ConstantCoefficient one(1.0);
-   ConstantCoefficient zero(5.0);
-   RepeatedCoefficient zero_vec(dim, zero);
-
-   LinearForm b(&fespace);
-   b.AddDomainIntegrator(new VectorDomainLFIntegrator(f_coeff));
-   b.AddBdrFaceIntegrator(
-      new VectorDGDirichletLFIntegrator(u_coeff, one, sigma, kappa));
-   b.Assemble();
-
-   GridFunction x(&fespace);
-   x = 0.0;
-
-   BilinearForm a(&fespace);
-   a.AddDomainIntegrator(new VectorDiffusionIntegrator(one));
-   a.AddInteriorFaceIntegrator(new VectorDGDiffusionIntegrator(
-                                  one, sigma, kappa, dim));
-   a.AddBdrFaceIntegrator(new VectorDGDiffusionIntegrator(
-                             one, sigma, kappa, dim));
-   a.Assemble();
-   a.Finalize();
-
-   const SparseMatrix &A = a.SpMat();
-#ifndef MFEM_USE_SUITESPARSE
-   GSSmoother M(A);
-   if (sigma == -1.0)
-   {
-      PCG(A, M, b, x, 1, 500, 1e-12, 0.0);
-   }
-   else
-   {
-      GMRES(A, M, b, x, 1, 500, 10, 1e-12, 0.0);
-   }
-#else
-   UMFPackSolver umf_solver;
-   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
-   umf_solver.SetOperator(A);
-   umf_solver.Mult(b, x);
-#endif
-
-   ParaViewDataCollection pv("DGDiffusion", &mesh);
-   pv.SetPrefixPath("ParaView");
-   pv.SetHighOrderOutput(true);
-   pv.SetLevelsOfDetail(order);
-   pv.RegisterField("u", &x);
-   pv.SetCycle(0);
-   pv.SetTime(0.0);
-   pv.Save();
-
-   cout << "L2 error: " << x.ComputeL2Error(u_coeff) << '\n';
-
-   return 0;
-}
-
-constexpr real_t pi = M_PI;
-constexpr real_t pi2 = pi*pi;
-
-real_t u_fn(const Vector &xvec)
-{
-   int dim = xvec.Size();
-   real_t x = pi*xvec[0], y = pi*xvec[1];
-   if (dim == 2) { return sin(x)*sin(y); }
-   else { real_t z = pi*xvec[2]; return sin(x)*sin(y)*sin(z); }
-}
-
-real_t f_fn(const Vector &xvec)
-{
-   int dim = xvec.Size();
-   real_t x = pi*xvec[0], y = pi*xvec[1];
-
-   if (dim == 2)
-   {
-      return 2*pi2*sin(x)*sin(y);
-   }
-   else // dim == 3
-   {
-      real_t z = pi*xvec[2];
-      return 3*pi2*sin(x)*sin(y)*sin(z);
-   }
-}

examples/vector-dg-diffusion.hpp

@@ -1,433 +0,0 @@
-#include "mfem.hpp"
-
-namespace mfem
-{
-
-class VectorDGDiffusionIntegrator : public BilinearFormIntegrator
-{
-protected:
-   Coefficient *Q = nullptr;
-   MatrixCoefficient *MQ = nullptr;
-   real_t sigma, kappa;
-   int vdim;
-
-   // these are not thread-safe!
-   Vector shape1, shape2, dshape1dn, dshape2dn, nor, nh, ni;
-   DenseMatrix jmat, dshape1, dshape2, mq, adjJ;
-
-public:
-   VectorDGDiffusionIntegrator(real_t s, real_t k, int vd=-1)
-      : sigma(s), kappa(k), vdim(vd) { }
-   VectorDGDiffusionIntegrator(Coefficient &q, real_t s, real_t k, int vd=-1)
-      : Q(&q), sigma(s), kappa(k), vdim(vd) { }
-   VectorDGDiffusionIntegrator(MatrixCoefficient &mq, real_t s, real_t k,
-                               int vd=-1)
-      : MQ(&mq), sigma(s), kappa(k), vdim(vd) { }
-
-   using BilinearFormIntegrator::AssembleFaceMatrix;
-   virtual void AssembleFaceMatrix(const FiniteElement &el1,
-                                   const FiniteElement &el2,
-                                   FaceElementTransformations &Trans,
-                                   DenseMatrix &full_elmat);
-};
-
-class VectorDGDirichletLFIntegrator : public LinearFormIntegrator
-{
-protected:
-   VectorCoefficient &uD;
-   Coefficient *Q = nullptr;
-   MatrixCoefficient *MQ = nullptr;
-   real_t sigma, kappa;
-   int vdim;
-
-   // these are not thread-safe!
-   Vector shape, dshape_dn, nor, nh, ni, uD_vec;
-   DenseMatrix dshape, mq, adjJ;
-
-public:
-   VectorDGDirichletLFIntegrator(VectorCoefficient &u, real_t s, real_t k,
-                                 int vd=-1)
-      : uD(u), sigma(s), kappa(k), vdim(vd) { }
-   VectorDGDirichletLFIntegrator(VectorCoefficient &u, Coefficient &q, real_t s,
-                                 real_t k, int vd=-1)
-      : uD(u), Q(&q), sigma(s), kappa(k), vdim(vd) { }
-   VectorDGDirichletLFIntegrator(VectorCoefficient &u, MatrixCoefficient &mq,
-                                 real_t s, real_t k, int vd=-1)
-      : uD(u), MQ(&mq), sigma(s), kappa(k), vdim(vd) { }
-
-   using LinearFormIntegrator::AssembleRHSElementVect;
-
-   void AssembleRHSElementVect(const FiniteElement &el,
-                               ElementTransformation &Tr,
-                               Vector &elvect) override
-   { MFEM_ABORT("Not implemented."); }
-
-   void AssembleRHSElementVect(const FiniteElement &el,
-                               FaceElementTransformations &Tr,
-                               Vector &elvect) override;
-};
-
-void VectorDGDiffusionIntegrator::AssembleFaceMatrix(
-   const FiniteElement &el1, const FiniteElement &el2,
-   FaceElementTransformations &Trans, DenseMatrix &full_elmat)
-{
-   int dim, ndof1, ndof2, ndofs;
-   bool kappa_is_nonzero = (kappa != 0.);
-   real_t w, wq = 0.0;
-
-   const int sdim = Trans.GetSpaceDim();
-   if (vdim < 0) { vdim = sdim; }
-
-   dim = el1.GetDim();
-   ndof1 = el1.GetDof();
-
-   nor.SetSize(dim);
-   nh.SetSize(dim);
-   ni.SetSize(dim);
-   adjJ.SetSize(dim);
-   if (MQ)
-   {
-      mq.SetSize(dim);
-   }
-
-   shape1.SetSize(ndof1);
-   dshape1.SetSize(ndof1, dim);
-   dshape1dn.SetSize(ndof1);
-   if (Trans.Elem2No >= 0)
-   {
-      ndof2 = el2.GetDof();
-      shape2.SetSize(ndof2);
-      dshape2.SetSize(ndof2, dim);
-      dshape2dn.SetSize(ndof2);
-   }
-   else
-   {
-      ndof2 = 0;
-   }
-
-   ndofs = ndof1 + ndof2;
-   DenseMatrix elmat;
-   elmat.SetSize(ndofs);
-   elmat = 0.0;
-   if (kappa_is_nonzero)
-   {
-      jmat.SetSize(ndofs);
-      jmat = 0.;
-   }
-
-   const IntegrationRule *ir = IntRule;
-   if (ir == nullptr)
-   {
-      // a simple choice for the integration order
-      int order;
-      if (ndof2)
-      {
-         order = 2 * std::max(el1.GetOrder(), el2.GetOrder());
-      }
-      else
-      {
-         order = 2 * el1.GetOrder();
-      }
-      ir = &IntRules.Get(Trans.GetGeometryType(), order);
-   }
-
-   // assemble: < {(Q \nabla u).n},[v] >      --> elmat
-   //           kappa < {h^{-1} Q} [u],[v] >  --> jmat
-   for (int p = 0; p < ir->GetNPoints(); p++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(p);
-
-      // Set the integration point in the face and the neighboring elements
-      Trans.SetAllIntPoints(&ip);
-
-      // Access the neighboring elements' integration points
-      // Note: eip2 will only contain valid data if Elem2 exists
-      const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
-      const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
-
-      if (dim == 1)
-      {
-         nor(0) = 2 * eip1.x - 1.0;
-      }
-      else
-      {
-         CalcOrtho(Trans.Jacobian(), nor);
-      }
-
-      el1.CalcShape(eip1, shape1);
-      el1.CalcDShape(eip1, dshape1);
-      w = ip.weight / Trans.Elem1->Weight();
-      if (ndof2)
-      {
-         w /= 2;
-      }
-      if (!MQ)
-      {
-         if (Q)
-         {
-            w *= Q->Eval(*Trans.Elem1, eip1);
-         }
-         ni.Set(w, nor);
-      }
-      else
-      {
-         nh.Set(w, nor);
-         MQ->Eval(mq, *Trans.Elem1, eip1);
-         mq.MultTranspose(nh, ni);
-      }
-      CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
-      adjJ.Mult(ni, nh);
-      if (kappa_is_nonzero)
-      {
-         wq = ni * nor;
-      }
-
-      // Note: in the jump term, we use 1/h1 = |nor|/det(J1) which is
-      // independent of Loc1 and always gives the size of element 1 in
-      // direction perpendicular to the face. Indeed, for linear transformation
-      //
-      //      |nor|=measure(face)/measure(ref. face),
-      //
-      //      det(J1)=measure(element)/measure(ref. element),
-      //
-      //      and the ratios measure(ref. element)/measure(ref. face)
-      //      are compatible for all element/face pairs.
-      //
-      // For example: meas(ref. tetrahedron)/meas(ref. triangle) = 1/3, and
-      // for any tetrahedron vol(tet)=(1/3)*height*area(base).
-      //
-      // For interior faces: q_e/h_e=(q1/h1+q2/h2)/2.
-
-      dshape1.Mult(nh, dshape1dn);
-      for (int i = 0; i < ndof1; i++)
-         for (int j = 0; j < ndof1; j++)
-         {
-            elmat(i, j) += shape1(i) * dshape1dn(j);
-         }
-
-      if (ndof2)
-      {
-         el2.CalcShape(eip2, shape2);
-         el2.CalcDShape(eip2, dshape2);
-         w = ip.weight / 2 / Trans.Elem2->Weight();
-         if (!MQ)
-         {
-            if (Q)
-            {
-               w *= Q->Eval(*Trans.Elem2, eip2);
-            }
-            ni.Set(w, nor);
-         }
-         else
-         {
-            nh.Set(w, nor);
-            MQ->Eval(mq, *Trans.Elem2, eip2);
-            mq.MultTranspose(nh, ni);
-         }
-         CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
-         adjJ.Mult(ni, nh);
-         if (kappa_is_nonzero)
-         {
-            wq += ni * nor;
-         }
-
-         dshape2.Mult(nh, dshape2dn);
-
-         for (int i = 0; i < ndof1; i++)
-            for (int j = 0; j < ndof2; j++)
-            {
-               elmat(i, ndof1 + j) += shape1(i) * dshape2dn(j);
-            }
-
-         for (int i = 0; i < ndof2; i++)
-            for (int j = 0; j < ndof1; j++)
-            {
-               elmat(ndof1 + i, j) -= shape2(i) * dshape1dn(j);
-            }
-
-         for (int i = 0; i < ndof2; i++)
-            for (int j = 0; j < ndof2; j++)
-            {
-               elmat(ndof1 + i, ndof1 + j) -= shape2(i) * dshape2dn(j);
-            }
-      }
-
-      if (kappa_is_nonzero)
-      {
-         // only assemble the lower triangular part of jmat
-         wq *= kappa;
-         for (int i = 0; i < ndof1; i++)
-         {
-            const real_t wsi = wq * shape1(i);
-            for (int j = 0; j <= i; j++)
-            {
-               jmat(i, j) += wsi * shape1(j);
-            }
-         }
-         if (ndof2)
-         {
-            for (int i = 0; i < ndof2; i++)
-            {
-               const int i2 = ndof1 + i;
-               const real_t wsi = wq * shape2(i);
-               for (int j = 0; j < ndof1; j++)
-               {
-                  jmat(i2, j) -= wsi * shape1(j);
-               }
-               for (int j = 0; j <= i; j++)
-               {
-                  jmat(i2, ndof1 + j) += wsi * shape2(j);
-               }
-            }
-         }
-      }
-   }
-
-   // elmat := -elmat + sigma*elmat^t + jmat
-   if (kappa_is_nonzero)
-   {
-      for (int i = 0; i < ndofs; i++)
-      {
-         for (int j = 0; j < i; j++)
-         {
-            real_t aij = elmat(i, j), aji = elmat(j, i), mij = jmat(i, j);
-            elmat(i, j) = sigma * aji - aij + mij;
-            elmat(j, i) = sigma * aij - aji + mij;
-         }
-         elmat(i, i) = (sigma - 1.) * elmat(i, i) + jmat(i, i);
-      }
-   }
-   else
-   {
-      for (int i = 0; i < ndofs; i++)
-      {
-         for (int j = 0; j < i; j++)
-         {
-            real_t aij = elmat(i, j), aji = elmat(j, i);
-            elmat(i, j) = sigma * aji - aij;
-            elmat(j, i) = sigma * aij - aji;
-         }
-         elmat(i, i) *= (sigma - 1.);
-      }
-   }
-
-   // populate full matrix following github issue #2909
-   full_elmat.SetSize(vdim*(ndof1 + ndof2));
-   full_elmat = 0.0;
-   for (int d=0; d<vdim; ++d)
-   {
-      for (int j=0; j<ndofs; ++j)
-      {
-         int jj = (j < ndof1) ? j + d*ndof1 : j - ndof1 + d*ndof2 + vdim*ndof1;
-         for (int i=0; i<ndofs; ++i)
-         {
-            int ii = (i < ndof1) ? i + d*ndof1 : i - ndof1 + d*ndof2 + vdim*ndof1;
-            full_elmat(ii, jj) += elmat(i, j);
-         }
-      }
-   }
-};
-
-void VectorDGDirichletLFIntegrator::AssembleRHSElementVect(
-   const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
-{
-   const int dim = el.GetDim();
-   const int sdim = Tr.GetSpaceDim();
-
-   if (vdim < 0) { vdim = sdim; }
-
-   const int ndof = el.GetDof();
-
-   bool kappa_is_nonzero = (kappa != 0.);
-   real_t w;
-
-   nor.SetSize(dim);
-   nh.SetSize(dim);
-   ni.SetSize(dim);
-   adjJ.SetSize(dim);
-   if (MQ)
-   {
-      mq.SetSize(dim);
-   }
-
-   shape.SetSize(ndof);
-   dshape.SetSize(ndof, dim);
-   dshape_dn.SetSize(ndof);
-
-   elvect.SetSize(vdim * ndof);
-   elvect = 0.0;
-
-   const IntegrationRule *ir = IntRule;
-   if (ir == NULL)
-   {
-      // a simple choice for the integration order; is this OK?
-      int order = 2*el.GetOrder();
-      ir = &IntRules.Get(Tr.GetGeometryType(), order);
-   }
-
-   for (int p = 0; p < ir->GetNPoints(); p++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(p);
-
-      // Set the integration point in the face and the neighboring element
-      Tr.SetAllIntPoints(&ip);
-
-      // Access the neighboring element's integration point
-      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
-
-      uD.Eval(uD_vec, Tr, ip);
-
-      if (dim == 1)
-      {
-         nor(0) = 2*eip.x - 1.0;
-      }
-      else
-      {
-         CalcOrtho(Tr.Jacobian(), nor);
-      }
-
-      el.CalcShape(eip, shape);
-      el.CalcDShape(eip, dshape);
-
-      // compute uD through the face transformation
-      w = ip.weight / Tr.Elem1->Weight();
-      if (!MQ)
-      {
-         if (Q)
-         {
-            w *= Q->Eval(*Tr.Elem1, eip);
-         }
-         ni.Set(w, nor);
-      }
-      else
-      {
-         nh.Set(w, nor);
-         MQ->Eval(mq, *Tr.Elem1, eip);
-         mq.MultTranspose(nh, ni);
-      }
-      CalcAdjugate(Tr.Elem1->Jacobian(), adjJ);
-      adjJ.Mult(ni, nh);
-
-      dshape.Mult(nh, dshape_dn);
-
-      for (int vd = 0; vd < vdim; ++vd)
-      {
-         for (int i = 0; i < ndof; ++i)
-         {
-            elvect[i + vd*ndof] += sigma * uD_vec[vd] * dshape_dn[i];
-         }
-      }
-      if (kappa_is_nonzero)
-      {
-         for (int vd = 0; vd < vdim; ++vd)
-         {
-            for (int i = 0; i < ndof; ++i)
-            {
-               elvect[i + vd*ndof] += kappa*(ni*nor) * uD_vec[vd] * shape[i];
-            }
-         }
-      }
-   }
-}
-
-} // namespace mfem

fem/fe_coll.cpp

@@ -401,6 +401,9 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
    {
       MFEM_ABORT("unknown FiniteElementCollection: " << name);
    }
+   MFEM_VERIFY(!strcmp(fec->Name(), name), "input name: \"" << name
+               << "\" does not match the created collection name: \""
+               << fec->Name() << '"');
 
    return fec;
 }

linalg/CMakeLists.txt

@@ -23,6 +23,7 @@ list(APPEND SRCS
   complex_operator.cpp
   constraints.cpp
   densemat.cpp
+  functional.cpp
   symmat.cpp
   handle.cpp
   matrix.cpp
@@ -49,6 +50,7 @@ list(APPEND HDRS
   complex_operator.hpp
   constraints.hpp
   densemat.hpp
+  functional.hpp
   dinvariants.hpp
   symmat.hpp
   dtensor.hpp

linalg/functional.cpp

@@ -0,0 +1,68 @@
+#include "functional.hpp"
+
+namespace mfem
+{
+QuadraticFunctional::QuadraticFunctional(const Operator *A_,
+                                         const Vector *b_, const real_t beta_, const real_t c_)
+   : Functional(A_ ? A_->Width() : 0)
+   , A(A_), beta(beta_), b(b_), c(c_)
+   , aux(A_ ? A_->Width() : 0)
+{
+   MFEM_VERIFY(A_ != nullptr,
+               "QuadraticFunctional: A must not be nullptr. "
+               << "Use QuadraticFunctional() constructor to create an empty Quadratic functional.");
+   MFEM_VERIFY(A_->Width() == A_->Height(),
+               "QuadraticFunctional: A must be a square operator.");
+   MFEM_VERIFY(b_ == nullptr || A_->Width() == b_->Size(),
+               "QuadraticFunctional: A and b must have compatible sizes");
+}
+
+void QuadraticFunctional::SetOperator(const Operator &A_)
+{
+   MFEM_VERIFY(A_.Width() == A_.Height(),
+               "QuadraticFunctional: A must be a square operator.");
+   A = &A_;
+   width = A_.Width();
+   aux.SetSize(width);
+}
+
+void QuadraticFunctional::SetVector(const Vector &b_, const real_t beta_)
+{
+   MFEM_VERIFY(A != nullptr && A->Width() == b_.Size(),
+               "QuadraticFunctional: A and b must have compatible sizes.");
+   b = &b_;
+   beta = beta_;
+}
+void QuadraticFunctional::SetConstant(real_t c_) { c = c_; }
+
+void QuadraticFunctional::Mult(const Vector &x, Vector &y) const
+{
+   MFEM_ASSERT(A != nullptr, "QuadraticFunctional: A must not be nullptr");
+   A->Mult(x, aux); // aux = A(x)
+   if (b != nullptr) { aux.Add(2.0*beta, *b); } // aux = A(x) + 2*beta*b
+   real_t result = 0.0;
+#ifdef MFEM_USE_MPI
+   if (IsParallel())
+   {
+      result = InnerProduct(GetComm(), x, aux); // result = <A(x) + 2*beta*b, x>
+   }
+   else
+   {
+      result = InnerProduct(x, aux); // result = <A(x) + 2*beta*b, x>
+   }
+#else
+   result = InnerProduct(x, aux); // result  = <A(x) + 2*beta*b, x>
+#endif
+   y.SetSize(1);
+   y[0] = result*0.5 + c; // y = 0.5 * <A(x) + 2*beta*b, x> + c
+}
+
+void QuadraticFunctional::EvalGradient(const Vector &x,
+                                       Vector &y) const
+{
+   y.SetSize(A->Width());
+   A->Mult(x, y);
+   if (b) { y.Add(beta,*b); }
+}
+
+} // namespace mfem

linalg/functional.hpp

@@ -0,0 +1,744 @@
+#ifndef MFEM_FUNCTIONAL_HPP
+#define MFEM_FUNCTIONAL_HPP
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+#include "../general/communication.hpp"
+#endif
+
+#include "operator.hpp"
+#include "blockvector.hpp"
+#include "solvers.hpp"
+#include <cxxabi.h>
+#include <vector>
+
+namespace mfem
+{
+/// @brief A base class for functionals F:R^n->R
+///
+/// This class provides an interface for evaluating
+/// $ F:R^n->R, \nabla F:R^n->R^n $, and $ \nabla^2 F:R^n x R^n ->R^n $.
+/// F.Mult(x, y) evaluates the functional at a point x, and stores the result in y[0]
+/// F.GetGradient() returns an operator that evaluates the gradient
+/// F.GetGradient().GetGradient(x) returns an Hessian action operator.
+///
+/// The usual Operator::GetGradient(const Vector &x) method for this method
+/// is deprecated as $ \nabla F $ only takes a single argument x.
+/// It is redundant to use F.GetGradient(x).Mult(x, y) to evaluate the gradient.
+/// Instead, use F.GetGradient().Mult(x, y) to evaluate the gradient at x.
+///
+/// The gradient and Hessian can be defined in two ways:
+/// 1. If the gradient is available as a sperate operator,
+///    then override Functional::GetGradient().
+///    In this case, Functional::HessianMult() will not be called.
+///
+/// 2. Otherwise, override Functional::EvalGradient() and Functional::HessianMult()
+///    to evaluate the gradient and Hessian action, respectively.
+///    The helper classes, GradientOperator and HessianActionOperator,
+///    will call these methods to evaluate the gradient and Hessian action.
+///    If Hessian is a seperate operator, then you can override
+///    The GradientOperator::GetGradient(x) will call Functional::GetHessian(x)
+///
+class Functional : public Operator
+{
+   Operator * riesz_map = nullptr; ///< Riesz map operator, if available
+public:
+   /// @brief Create a Functional with optional gradient and hessian
+   /// @param n number of variables
+   Functional(int n=0)
+      : Operator(1, n)
+      , grad_operator(*this)
+      , hessian_action_operator(*this)
+   { }
+
+#ifdef MFEM_USE_MPI
+   Functional(MPI_Comm comm, int n=0)
+      : Functional(n)
+   { SetComm(comm); }
+
+   void SetComm(MPI_Comm comm_)
+   {
+      parallel = comm_ != MPI_COMM_NULL;
+      comm = comm_;
+   }
+   MPI_Comm GetComm() const { return comm; }
+   bool IsParallel() const { return parallel; }
+#else
+   constexpr bool IsParallel() const { return false; }
+#endif
+
+   void SetRieszMap(Operator &op) { riesz_map = &op; }
+   /// @brief return the GradientOperator that evaluates the gradient
+   ///        input x is not used. Use GetGradient().Mult(x,y) to evaluate the gradient
+   ///        we recommend using GetGradient() instead of GetGradient(x)
+   /// Deprecated. See Functional::GetGradient()
+   MFEM_DEPRECATED
+   Operator &GetGradient(const Vector &dummy) const override final { return GetGradient(); }
+
+   /// @brief Return the GradientOperator that wraps Functional::EvalGradient() for Mult().
+   /// @note If the functional has a corresponding standalone gradient operator,
+   /// override this method to return the gradient operator.
+   virtual Operator &GetGradient() const { return grad_operator; }
+
+   /// @brief Evaluate the functional at a point x that will be called by GradientOperator::Mult()
+   /// @note This method is not meant to be called directly. See, GradientOperator
+   virtual void EvalGradient(const Vector &x, Vector &y) const
+   {
+      MFEM_ABORT("Functional::EvalGradient() not implemented");
+   }
+   /// @brief Evaluate the Hessian action at a point x and direction d
+   /// that will be called by Functional::GetGradient().GetHessian(x).Mult(d,y)
+   /// @note This method is not meant to be called directly. See, HessianActionOperator
+   virtual void HessianMult(const Vector &x, const Vector &d, Vector &y) const
+   {
+      MFEM_ABORT("Functional::HessianMult() not implemented.");
+   }
+   /// @brief Return the HessianActionOperator at evaluation point x
+   ///        that wraps Functional::HessianMult() for Mult().
+   ///        See, HessianActionOperator and Functional::HessianMult().
+   ///
+   /// @note If the Hessian is available as a seperate operator, override this method.
+   ///
+   /// @warning If GetGradient() is overridden, this method will not be used.
+   virtual Operator &GetHessian(const Vector &x) const
+   {
+      hessian_action_operator.SetX(x);
+      return hessian_action_operator;
+   }
+private:
+#ifdef MFEM_USE_MPI
+   bool parallel=false;
+#else
+   const static bool parallel=false;
+#endif
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm;
+#endif
+   /// @brief A helper class to return an operator that evaluates the gradient
+   ///        using Functional::EvalGradient() method.
+   class GradientOperator : public Operator
+   {
+   private: const Functional &f; mutable Vector der;
+   public:
+      GradientOperator(const Functional &f) : Operator(f.Width()), f(f) {}
+      /// @brief Evaluate the gradient of Functional at a point x
+      void Mult(const Vector &x, Vector &y) const override final
+      {
+         if (f.riesz_map)
+         {
+            der.SetSize(f.Width());
+            f.EvalGradient(x, der);
+            f.riesz_map->Mult(der, y);
+         }
+         else
+         {
+            f.EvalGradient(x, y);
+         }
+      }
+      /// @brief Evaluate the Hessian of Functional at a point x
+      Operator &GetGradient(const Vector &x) const override final { return f.GetHessian(x); }
+   };
+   friend class GradientOperator;
+
+   /// @brief A helper class to return an operator that applies the Hessian action
+   ///        using Functional::HessianMult() method.
+   class HessianActionOperator : public Operator
+   {
+   private:
+      const Functional &f;
+      const Vector *x;
+   public:
+      HessianActionOperator(const Functional &f) : Operator(f.Width()), f(f) {}
+      void SetX(const Vector &new_x) { x = &new_x; }
+      void Mult(const Vector &d, Vector &y) const override { f.HessianMult(*x, d, y); }
+   };
+   friend class HessianActionOperator;
+
+   mutable GradientOperator grad_operator;
+   mutable HessianActionOperator hessian_action_operator;
+};
+
+/// @brief Stacked functioanl operator, [f1, ..., fk] where fi:R^n->R are functionals
+/*
+   Typical usage of this class is to provide a single operator for multiple constraints.
+   For example, consider a minimization problem with k constraints,
+   min f0(u) s.t. fi(u)=0, i=1,...,k.
+   The Lagrangian functional is
+   $ L(u, lambda) = F0(u) + sum_i lambda_i * fi(u) $
+   The first-order optimality conditions are
+   $ \nabla f0(u) + \sum_i lambda_i * grad fi(u) = 0 $
+   $ fi(u) = 0 $
+   where lambda_i are the Lagrange multipliers.
+   The StackedFunctional class can be used to represent the list of constraints fi(u).
+
+   StackedFunctional::Mult(u, y) will evaluate each functional y[i]=fi(u)
+
+   StackedFunctional::GetGradient(u) represents an operator, column-stacked gradient
+   That is, [grad f0(u), ..., grad fk(u)] in R^{n x k}
+   If you want to extract the gradient as a matrix, use
+   StackedFunctional::GetGradientMatrix(const Vector &x, DenseMatrix &grad)
+   As functionals are not assumed to return a sparse vector, the gradient is dense.
+
+   StackedFunctional::GetGradient(u).Mult(lambda, y) contract the gradients with the Lagrange multipliers
+   y = sum lambda_i * grad fi(u)
+   StackedFunctional::GetGradient(u).MultTranspose(d, y) return the directional derivative for each k
+   y[i] = <grad fi(u), d>
+
+   StackedFunctional::GetHessian(u, lambda).Mult(d, y) will return the contracted Hessian action
+   $ y = \sum_i \lambda_i * H_{fi}(u, d) $
+*/
+/// @warning Functionals should be all serial or all parallel.
+///
+class StackedFunctional : public Operator
+{
+public:
+   StackedFunctional(int n=0)
+      : Operator(0, n)
+      , funcs(0)
+      , grad_helper_op(*this)
+      , hessian_helper_op(*this)
+   {}
+   StackedFunctional(Functional &f)
+      : Operator(0, f.Width())
+      , grad_helper_op(*this)
+      , hessian_helper_op(*this)
+   { AddFunctional(f); }
+
+   StackedFunctional(const std::vector<Functional*> &funcs)
+      : Operator((int)funcs.size(), funcs[0]->Width())
+      , grad_helper_op(*this)
+      , hessian_helper_op(*this)
+   { for (auto &f : funcs) { AddFunctional(*f); } }
+
+   void AddFunctional(Functional &f)
+   {
+#ifdef MFEM_USE_MPI
+      if (funcs.empty()) { if (f.IsParallel()) { SetComm(f.GetComm()); } }
+#endif
+      MFEM_VERIFY(f.Width() == Width(),
+                  "StackedFunctional::AddFunctional: Functional width does not match with the operator.");
+      MFEM_VERIFY(parallel == f.IsParallel(),
+                  "StackedFunctional::AddFunctional: Parallelism mismatch.");
+      funcs.push_back(&f);
+      height++;
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      y.SetSize(Height());
+      Vector yview;
+      for (int i=0; i<Height(); i++)
+      {
+         yview.MakeRef(y, i, 1);
+         funcs[i]->Mult(x, yview);
+      }
+   }
+
+   Operator &GetGradient(const Vector &x) const override
+   {
+      grad_helper_op.SetX(x);
+      return grad_helper_op;
+   }
+
+   void GetGradientMatrix(const Vector &x, DenseMatrix &grads) const
+   {
+      grads.SetSize(Width(), Height());
+      Vector grad;
+      for (int i=0; i<Height(); i++)
+      {
+         grads.GetColumnReference(i, grad);
+         funcs[i]->GetGradient().Mult(x, grad);
+      }
+   }
+   Functional &GetFunctional(int i) const
+   {
+      MFEM_VERIFY(i >= 0 && i < Height(),
+                  "StackedFunctional::GetFunctional: Index out of bounds.");
+      return *funcs[i];
+   }
+
+   Operator &GetHessian(const Vector &x, const Vector &lambda) const
+   {
+      hessian_helper_op.SetX(x, lambda);
+      return hessian_helper_op;
+   }
+
+   bool parallel;
+   bool IsParallel() const { return parallel; }
+#ifdef MFEM_USE_MPI
+   void SetComm(MPI_Comm comm_)
+   {
+      parallel = comm != MPI_COMM_NULL;
+      comm = comm_;
+   }
+   MPI_Comm GetComm() const { return comm; }
+#endif
+
+protected:
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm;
+#endif
+   std::vector<Functional*> funcs;
+
+   class GradientOperator : public Operator
+   {
+   public:
+      GradientOperator(const StackedFunctional &op)
+         : Operator(op.Width(), op.Height())
+         , op(op)
+         , tmp_grad(op.Width())
+      {}
+      void SetX(const Vector &x) const { x_curr = &x; }
+      Operator &GetGradient(const Vector &lambda) const override
+      {
+         op.hessian_helper_op.SetX(*x_curr, lambda);
+         return op.hessian_helper_op;
+      }
+
+      void Mult(const Vector &lambda, Vector &y) const override
+      {
+         y.SetSize(op.Width());
+         y = 0.0;
+         for (int i=0; i<op.Height(); i++)
+         {
+            op.funcs[i]->GetGradient().Mult(*x_curr, tmp_grad);
+            y.Add(lambda[i], tmp_grad);
+         }
+      }
+      void MultTranspose(const Vector &x, Vector &y) const override
+      {
+         y.SetSize(op.Height());
+         for (int i=0; i<op.Height(); i++)
+         {
+            op.funcs[i]->GetGradient().Mult(x, tmp_grad);
+            y[i] = InnerProduct(tmp_grad, *x_curr);
+         }
+#ifdef MFEM_USE_MPI
+         if (op.IsParallel())
+         {
+            MPI_Allreduce(MPI_IN_PLACE, y.GetData(), op.Height(),
+                          MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                          op.GetComm());
+         }
+#endif
+      }
+
+   private:
+      const StackedFunctional &op;
+      mutable const Vector *x_curr;
+      mutable Vector tmp_grad;
+
+   };
+   class HessianActionOperator : public Operator
+   {
+   public:
+      HessianActionOperator(const StackedFunctional &op)
+         : Operator(op.Width()), op(op)
+      {}
+      void SetX(const Vector &x, const Vector &lambda) const { x_curr = &x; lambda_curr = &lambda; }
+
+      void Mult(const Vector &d, Vector &y) const override
+      {
+         y.SetSize(op.Width());
+         y = 0.0;
+         for (int i=0; i<op.Height(); i++)
+         {
+            op.funcs[i]->GetGradient().GetGradient(*x_curr).Mult(d, tmp_hessian);
+            y.Add((*lambda_curr)[i], tmp_hessian);
+         }
+      }
+   private:
+      const StackedFunctional &op;
+      mutable Vector tmp_hessian;
+      mutable const Vector *x_curr;
+      mutable const Vector *lambda_curr;
+   };
+   friend class GradientOperator;
+   friend class HessianActionOperator;
+   mutable GradientOperator grad_helper_op;
+   mutable HessianActionOperator hessian_helper_op;
+private:
+};
+
+class ConstrainedOptimizationProblem : public Functional
+{
+public:
+   ConstrainedOptimizationProblem(Functional &objective_,
+                                  Operator *eq_constraints_=nullptr,
+                                  Operator *ineq_constraints_=nullptr)
+      : Functional(objective_.Width())
+      , objective(objective_)
+      , eq_constraints(eq_constraints_)
+      , ineq_constraints(ineq_constraints_)
+   {
+      // Check Size
+      MFEM_VERIFY((eq_constraints == nullptr ||
+                   eq_constraints->Width() == objective.Width()),
+                  "ConstrainedFunctional: Equality constraints width does not match with the objective.");
+      MFEM_VERIFY((ineq_constraints == nullptr ||
+                   ineq_constraints->Width() == objective.Width()),
+                  "ConstrainedFunctional: Inequality constraints width does not match with the objective.");
+#ifdef MFEM_USE_MPI
+      if (objective.IsParallel()) { SetComm(objective.GetComm()); }
+#endif
+   }
+
+   Functional &GetObjective() { return objective; }
+   const Functional &GetObjective() const { return objective; }
+   Operator *GetEqualityConstraints() { return eq_constraints; }
+   const Operator *GetEqualityConstraints() const { return eq_constraints; }
+   Operator *GetInequalityConstraints() { return ineq_constraints; }
+   const Operator *GetInequalityConstraints() const { return ineq_constraints; }
+protected:
+   Functional &objective;
+   Operator *eq_constraints;
+   Operator *ineq_constraints;
+};
+
+/// @brief A Lagrangian functional for
+/// min F(u)
+/// subject to C(u) = 0
+/// That is, L(u, lambda) = F(u) + <lambda, C(u)>
+///
+/// We assume that $ F:R^n -> R $ is a functional,
+///                $ C:R^n -> R^k $ is an equality constraint operator,
+/// C should return a residual. That is,
+/// C(u) = c, then C.Mult(u, y) should return y[i] = C_i(u) - c_i.
+///
+/// C.GetGradient(u):R^k -> R^n that takes lambda and returns the contracted gradient at x
+/// C.GetGradient(u).Mult(lambda, y) returns y = sum lambda_i * grad C_i(u)
+///
+/// C's gradient should support MultTranspose method
+/// That is, C.GetGradient(u).MultTranspose(d, y) returns y[i] = <grad C_i(u), d>
+///
+class LagrangianFunctional : public ConstrainedOptimizationProblem
+{
+private:
+   mutable Vector eq_residual;
+public:
+   LagrangianFunctional(Functional &objective,
+                        Operator &eq_constraints)
+      : ConstrainedOptimizationProblem(objective, &eq_constraints)
+      , eq_residual(eq_constraints.Height())
+   {
+      width = objective.Width() + eq_constraints.Height();
+
+      offsets.SetSize(3);
+      offsets[0] = 0;
+      offsets[1] = objective.Width();
+      offsets[2] = eq_constraints.Height();
+      offsets.PartialSum();
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      const BlockVector input_block(const_cast<Vector&>(x), offsets);
+      const Vector &u = input_block.GetBlock(0);
+      const Vector &lambda = input_block.GetBlock(1);
+
+      y.SetSize(1);
+      y[0] = 0.0;
+      objective.Mult(u, y);
+      eq_constraints->Mult(u, eq_residual);
+      y[0] += InnerProduct(lambda, eq_residual);
+   }
+
+   void EvalGradient(const Vector &x, Vector &y) const override
+   {
+      const BlockVector input_block(const_cast<Vector&>(x), offsets);
+      const Vector &u = input_block.GetBlock(0);
+      const Vector &lambda = input_block.GetBlock(1);
+
+      y.SetSize(Width());
+      BlockVector output_block(y, offsets);
+      Vector &opt_residual = output_block.GetBlock(0);
+      eq_residual = output_block.GetBlock(1);
+      y = 0.0;
+      // grad F(u) + \sum_i lambda_i grad C_i(u)
+      objective.GetGradient().Mult(u, opt_residual);
+      eq_constraints->GetGradient(u).AddMult(lambda, opt_residual);
+      // grad C_i(u)^T
+      eq_constraints->GetGradient(u).MultTranspose(u, eq_residual);
+   }
+
+   /// @brief Evaluate the Hessian action at a point x=[u, lambda]
+   /// and direction d=[v, mu]
+   /// $ [H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + <\nabla C_i(u), mu> $
+   void HessianMult(const Vector &x, const Vector &d, Vector &y) const override
+   {
+      const BlockVector input_block(const_cast<Vector&>(x), offsets);
+      const Vector &u = input_block.GetBlock(0);
+      const Vector &lambda = input_block.GetBlock(1);
+
+      const BlockVector direction_block(const_cast<Vector&>(x), offsets);
+      const Vector &v = direction_block.GetBlock(0);
+      const Vector &mu = direction_block.GetBlock(1);
+
+      y.SetSize(Width());
+      BlockVector output_block(y, offsets);
+      Vector &opt_H = output_block.GetBlock(0); // Optimality Hessian
+      Vector &eq_H = output_block.GetBlock(1); // Equality Hessian
+
+      // H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + mu_i grad C_i(u)
+      objective.GetGradient().GetGradient(u).Mult(v, opt_H);
+      eq_constraints->GetGradient(u).GetGradient(lambda).AddMult(v, opt_H);
+      eq_constraints->GetGradient(u).Mult(mu, eq_H);
+      // <grad C_i(u), d>
+      eq_constraints->GetGradient(u).MultTranspose(d, eq_H);
+   }
+
+protected:
+   Array<int> offsets; // offsets for [x, lambda, mu]
+};
+
+/// @brief An augmented Lagrangian functional of the form
+/// F(u) + 0.5 mu * ||C(u)||^2 + <lambda, C(u)>
+/// where F is the objective functional,
+/// C is the equality constraint operator,
+/// lambda is the Lagrange multiplier vector (initialized to zero),
+/// mu is the penalty parameter (defaults to 1.0)
+///
+/// Currently, only equality constraints are supported.
+///
+/// AugLagrangianFunctional::Update() will update the penalty and Lagrange multiplier vectors
+/// By default, lambda <- lambda + mu * C(u)
+///             mu <- mu (no update)
+class AugLagrangianFunctional : public ConstrainedOptimizationProblem
+{
+public:
+   AugLagrangianFunctional(Functional &objective_,
+                           Operator &eq_constraints_)
+      : ConstrainedOptimizationProblem(objective_, &eq_constraints_)
+      , lambda(eq_constraints_.Height())
+      , mu(1.0)
+      , eq_residual(eq_constraints_.Height())
+      , eq_dir(eq_constraints_.Height())
+   {
+      lambda = 0.0;
+   }
+
+   void SetLambda(const Vector &lambda_)
+   {
+      MFEM_VERIFY(lambda_.Size() == eq_constraints->Height(),
+                  "AugLagrangianFunctional: Lambda size does not match with the equality constraints.");
+      lambda = lambda_;
+   }
+
+   void SetPenalty(real_t mu_)
+   {
+      MFEM_VERIFY(mu_ >= 0.0,
+                  "AugLagrangianFunctional: Penalty parameter mu must be non-negative.");
+      mu = mu_;
+   }
+
+   virtual void Update(const Vector &x)
+   {
+      // Update the Lagrange multipliers
+      eq_constraints->AddMult(x, lambda, mu);
+      // Update the penalty parameter
+      // Do nothing
+   }
+   const Vector &GetLambda() const { return lambda; }
+   real_t GetPenalty() const { return mu; }
+
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      y.SetSize(1);
+      objective.Mult(x, y);
+      eq_constraints->Mult(x, eq_residual);
+      y[0] += lambda*eq_residual;
+      y[0] += 0.5 * mu * (eq_residual*eq_residual);
+   }
+
+   void EvalGradient(const Vector &x, Vector &y) const override
+   {
+      y.SetSize(Width());
+      // grad F(x) + \sum_i (lambda_i + mu * C_i(x)) grad C_i(x)
+      Vector curr_lambda = lambda; // store lambda + mu * C(x)
+      objective.GetGradient().Mult(x, y);
+      eq_constraints->Mult(x, eq_residual);
+      curr_lambda.Add(mu, eq_residual);
+      eq_constraints->GetGradient(x).AddMult(curr_lambda, y);
+   }
+
+   /// @brief Evaluate the Hessian action at a point x=[u, lambda]
+   /// and direction d=[v, mu]
+   /// $ H_F(u,d) + \sum_i \lambda_i H_{C_i}(u, d) + <\nabla C_i(u), mu> $
+   void HessianMult(const Vector &x, const Vector &d, Vector &y) const override
+   {
+      // H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + mu_i grad C_i(u) <grad C_i(u), d>
+      objective.GetGradient().GetGradient(x).Mult(d, y);
+
+      Vector curr_lambda = lambda;
+      eq_constraints->Mult(x, eq_residual);
+      curr_lambda.Add(mu, eq_residual);
+
+      eq_constraints->GetGradient(x).GetGradient(curr_lambda).AddMult(d, y);
+      // eq_dir = <grad C_i(u), d>
+      eq_constraints->GetGradient(x).MultTranspose(d, eq_dir);
+      // mu_i <grad C_i(u), eq_dir>
+      eq_constraints->GetGradient(x).AddMult(eq_dir, y, mu);
+   }
+
+protected:
+   Vector lambda;
+   real_t mu;
+   mutable Vector eq_residual; // residual of the equality constraints, R^k
+   // directional derivative of the equality constraints, R^k
+   mutable Vector eq_dir;
+};
+
+/// @brief Quadratic functional of the form
+///        f(u) = 0.5 * <A u, u> + beta<b, u> + c
+///        where A is a square (possibly nonlinear) operator,
+///        beta is a scalar (defaults to 1.0, not used when b is nullptr),
+///        b is a vector (independent of u, optional),
+///        c is a constant (independent of u, optional).
+///        GetHessian() returns the operator A.
+///
+class QuadraticFunctional : public Functional
+{
+public:
+   QuadraticFunctional()
+      : Functional(0)
+      , A(nullptr), b(nullptr), c(0.0)
+   {}
+   QuadraticFunctional(const Operator *A_, const Vector *b_=nullptr,
+                       const real_t beta_=1.0,
+                       const real_t c_=0.0);
+#ifdef MFEM_USE_MPI
+   QuadraticFunctional(MPI_Comm comm_)
+      : QuadraticFunctional()
+   { SetComm(comm_); }
+
+   QuadraticFunctional(MPI_Comm comm_, const Operator *A_,
+                       const Vector *b_=nullptr,
+                       const real_t beta_=1.0, const real_t c_=0.0)
+      : QuadraticFunctional(A_, b_, beta_, c_)
+   { SetComm(comm_); }
+#endif
+
+   void SetOperator(const Operator &A_);
+   void SetVector(const Vector &b_, const real_t beta_=1.0);
+   void SetConstant(real_t c_);
+
+   void Mult(const Vector &x, Vector &y) const override;
+   void EvalGradient(const Vector &x, Vector &y) const override;
+protected:
+   const Operator *A;
+   real_t beta;
+   const Vector *b;
+   real_t c;
+   mutable Vector aux;
+protected:
+   /// @brief return the underlying Operator A
+   /// @warning Modifying the returned operator leads to undefined behavior.
+   Operator& GetHessian(const Vector &dummy) const override
+   {
+      return const_cast<Operator&>(*A);
+   }
+};
+
+class Optimizer : public IterativeSolver
+{
+public:
+   Optimizer() : IterativeSolver(), f(nullptr) { }
+#ifdef MFEM_USE_MPI
+   Optimizer(MPI_Comm comm) : IterativeSolver(comm), f(nullptr) { }
+#endif
+   // @brief Set the subproblem functional operator
+   // @param op the functional operator
+   // @note The functional will be stored in subproblem, and oper will be set to the gradient of the functional.
+   void SetOperator(const Functional &f_)
+   {
+      f = &f_;
+      IterativeSolver::SetOperator(f_.GetGradient());
+   }
+   virtual void SetLinearSolver(Solver &prec) { IterativeSolver::SetPreconditioner(prec); }
+
+   /// @brief This will abort. Should be called only with a Functional operator.
+   void SetOperator(const Operator &op) override
+   {
+      MFEM_ABORT("OptSolver::SetOperator() should not be called directly. Use SetFunctional() instead.");
+   }
+protected:
+   const Functional * f;
+};
+
+class NewtonOptimizer : public Optimizer
+{
+private:
+   real_t step_size = 1.0; // default step size
+public:
+   NewtonOptimizer() : Optimizer() { }
+#ifdef MFEM_USE_MPI
+   NewtonOptimizer(MPI_Comm comm) : Optimizer(comm) { }
+#endif
+   void SetStepSize(real_t step_size_) { step_size = step_size_; }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      dx.SetSize(x.Size());
+      y.SetSize(x.Size());
+      y = x;
+      MFEM_ASSERT(f != nullptr,
+                  "NewtonOptimizer::Mult() called without a functional operator.");
+      MFEM_ASSERT(prec != nullptr,
+                  "NewtonOptimizer::Mult() called without a linear solver.");
+      for (int i=0; i<max_iter; i++)
+      {
+         oper->Mult(y, grad);
+         Operator &hess = oper->GetGradient(y);
+         prec->SetOperator(hess);
+         prec->Mult(grad, dx);
+         y.Add(-step_size, dx);
+         if (Dot(dx, dx) < abs_tol*abs_tol)
+         {
+            break;
+         }
+      }
+   }
+private:
+   mutable Vector grad;
+   mutable Vector dx;
+};
+
+class GradientDescentOptimizer : public Optimizer
+{
+private:
+   real_t step_size = 1.0; // default step size
+public:
+   GradientDescentOptimizer() : Optimizer() { }
+#ifdef MFEM_USE_MPI
+   GradientDescentOptimizer(MPI_Comm comm) : Optimizer(comm) { }
+#endif
+   void SetStepSize(real_t step_size_) { step_size = step_size_; }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      grad.SetSize(x.Size());
+      y.SetSize(x.Size());
+      y = x;
+      MFEM_ASSERT(f != nullptr,
+                  "NewtonOptimizer::Mult() called without a functional operator.");
+      MFEM_ASSERT(prec != nullptr,
+                  "NewtonOptimizer::Mult() called without a linear solver.");
+      for (int i=0; i<max_iter; i++)
+      {
+         oper->Mult(y, grad);
+         y.Add(-step_size, grad);
+         if (Dot(grad, grad) < abs_tol*abs_tol)
+         {
+            break;
+         }
+      }
+   }
+private:
+   mutable Vector grad;
+};
+
+
+} // namespace mfem
+#endif // MFEM_FUNCTIONAL_HPP