- added mma test

- eleminated redundancies

examples/elasticity.cpp

@@ -0,0 +1,362 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+// Problem 1 - L-shaped domain with homogeneous dirichlet and inhomogeneous traction
+// make elasticity -j && ./elasticity -rs 5 -o 2 -m lshaped.mesh -prob 1
+// Problem 2 - Plate with hole domain with homogeneous dirichlet and inhomogeneous traction
+// make elasticity -j && ./elasticity -rs 3 -o 2 -m platewhole.mesh -prob 2
+// Problem 3 - Plate with hole but homogeneous and inhomogeneous dirichlet
+// make elasticity -j && ./elasticity -rs 3 -o 2 -m platewhole.mesh -prob 3
+using namespace mfem;
+using namespace std;
+
+//----------------------------------------------------------------
+// Forcing function: Zero body force.
+void ZeroForce(const Vector &x, Vector &f)
+{
+   f.SetSize(x.Size());
+   f = 0.0;
+}
+
+//----------------------------------------------------------------
+// Traction function: Apply a constant (nonzero) traction.
+// For example, here we apply a constant upward traction.
+// (You can adjust the constant vector as needed.)
+void ConstantTractionY(const Vector &x, Vector &g)
+{
+   g.SetSize(x.Size());
+   g = 0.0;
+   // Let the vertical component be nonzero.
+   // For a 2D problem, we set g = [0, 1] (traction in the y-direction).
+   g[1] = 1.0;
+}
+
+void ConstantTractionX(const Vector &x, Vector &g)
+{
+   g.SetSize(x.Size());
+   g = 0.0;
+   // Let the vertical component be nonzero.
+   // For a 2D problem, we set g = [0, 1] (traction in the y-direction).
+   g[0] = 1.0;
+}
+
+void ConstantTractionNegX(const Vector &x, Vector &g)
+{
+   g.SetSize(x.Size());
+   g = 0.0;
+   // Let the vertical component be nonzero.
+   // For a 2D problem, we set g = [0, 1] (traction in the y-direction).
+   g[0] = -1.0;
+}
+
+// Unit displacement on x=1 smoothly going to zero at y=0 and y=1
+void ConstantDisplacementPlateHoleX1(const Vector &x, Vector &g)
+{
+   g.SetSize(x.Size());
+   g = 0.0;
+   double alpha = 10;
+   double dx1 = std::pow(x[0]-1.0,2.0);
+   double dx2 = std::pow(x[0]-0.0,2.0);
+   g[0] = 5.0*std::exp(-alpha*(x[0]-0.5)*(x[0]-0.5))*x[1]*(1.0-x[1]);
+   g[0] = std::tanh(alpha*(x[0]-0.5))*x[1]*(1.0-x[1]);
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   OptionsParser args(argc, argv);
+   const char *mesh_file = "lshaped.mesh"; // Your L-shaped mesh file.
+   int order = 2;            // Finite element order
+   int ref_levels = 3;       // Number of uniform refinements
+   bool visualization = true;
+   int problem = 3;
+
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Polynomial degree of mesh finite element space.");
+   args.AddOption(&ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&problem, "-prob", "--problem type"," ");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Uniform mesh refinements.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define a vector finite element space (each displacement component is in H1).
+   FiniteElementCollection *fec = new H1_FECollection(order, dim);
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec, dim);
+
+   // 5. Mark essential (Dirichlet) boundaries.
+   // Here we assume boundary attribute 1 is for Dirichlet conditions (clamped boundaries).
+   Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   ess_bdr = 0;
+   // ess_bdr[3] = 1;
+   if (problem == 1)
+   {
+      ess_bdr[3] = 1;
+   }
+   if (problem == 3)
+   {
+      ess_bdr = 1;
+      ess_bdr[4] = 0.0;
+      // ess_bdr[4] = 1;
+   }
+   Array<int> ess_tdof_list;
+   fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+   // 6. Define the GridFunction (solution) and initialize it to zero.
+   GridFunction x(fespace);
+   x = 0.0;
+   VectorFunctionCoefficient dbc_coeff(dim, ConstantDisplacementPlateHoleX1);
+   if (problem == 3)
+   {
+      x.ProjectBdrCoefficient(dbc_coeff, ess_bdr);
+   }
+
+   // 7. Set up the right-hand side linear form.
+   LinearForm *b = new LinearForm(fespace);
+
+   // Domain integrator for the body force: here, zero.
+   VectorFunctionCoefficient f_coeff(dim, ZeroForce);
+   b->AddDomainIntegrator(new VectorDomainLFIntegrator(f_coeff));
+
+   // Neumann boundary: apply a nonzero (constant) traction on attribute 2.
+   Array<int> neumann_bdr(mesh->bdr_attributes.Max());
+   neumann_bdr = 0;
+   if (problem == 1 || problem == 2)
+   {
+      neumann_bdr[1] = 1;
+   }
+   VectorFunctionCoefficient g_coeff(dim, ConstantTractionX);
+   b->AddBdrFaceIntegrator(new VectorBoundaryLFIntegrator(g_coeff), neumann_bdr);
+
+   Array<int> neumann_bdr_negx(mesh->bdr_attributes.Max());
+   neumann_bdr_negx = 0;
+   neumann_bdr_negx[3] = 1;
+   VectorFunctionCoefficient g_coeff_negx(dim, ConstantTractionNegX);
+   if (problem == 2)
+   {
+      b->AddBdrFaceIntegrator(new VectorBoundaryLFIntegrator(g_coeff_negx), neumann_bdr_negx);
+   }
+
+   b->Assemble();
+
+   // 8. Set up the bilinear form corresponding to the linear elasticity operator.
+   BilinearForm *a = new BilinearForm(fespace);
+
+   // Define the material properties.
+   // For a given Young's modulus E and Poisson's ratio nu:
+   double E = 1.0;    // Young's modulus (you can adjust as needed)
+   double nu = 0.3;   // Poisson's ratio
+   double mu = E / (2.0 * (1.0 + nu));
+   double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0*nu));
+   ConstantCoefficient lc(lambda);
+   ConstantCoefficient muc(mu);
+
+   // MFEM's ElasticityIntegrator implements:
+   //   a(u,v) = ∫ [λ (div u)(div v) + 2μ ε(u):ε(v)] dx.
+   a->AddDomainIntegrator(new ElasticityIntegrator(lc, muc));
+   a->Assemble();
+
+   // 9. Form the linear system A X = B, eliminating the Dirichlet DOFs.
+   SparseMatrix A;
+   Vector X, Bvec;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, Bvec);
+
+   // 10. Solve the linear system.
+   // Here we use the Conjugate Gradient (CG) method with a Gauss-Seidel smoother.
+   GSSmoother M(A);
+   PCG(A, M, Bvec, X, 1, 1000, 1e-12, 0.0);
+
+   // 11. Recover the finite element solution.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 12. Save the solution for visualization.
+   {
+      ofstream sol_ofs("solution.gf");
+      x.Save(sol_ofs);
+   }
+
+   // Optionally, visualize the solution in GLVis.
+   if (visualization)
+   {
+      osockstream sock(19916, "localhost");
+      sock << "solution\n";
+      mesh->Print(sock);
+      x.Save(sock);
+      sock.send();
+      sock << "window_title 'Elasticity: Displacement'\n"
+           << "window_geometry "
+           << 0 << " " << 0 << " " << 400 << " " << 400 << "\n"
+           << "keys jRmclA" << endl;
+   }
+
+   FiniteElementCollection *fec_dc = new L2_FECollection(order-1, dim);
+   FiniteElementSpace *fespace_dc = new FiniteElementSpace(mesh, fec_dc);
+   GridFunction stress_gf(fespace_dc);
+   GridFunction trace_strain(fespace_dc);
+
+   FiniteElementCollection *fec_dc2 = new L2_FECollection(0, dim);
+   FiniteElementSpace *fespace_dc2 = new FiniteElementSpace(mesh, fec_dc2);
+   GridFunction stress_gf2(fespace_dc2);
+   stress_gf = 0.0;
+
+   // Loop over elements.
+   for (int i = 0; i < mesh->GetNE(); i++)
+   {
+      const FiniteElement *fe = fespace->GetFE(i);
+      const FiniteElement *fe2 = fespace_dc->GetFE(i);
+      ElementTransformation *T = mesh->GetElementTransformation(i);
+      int nd = fe->GetDof();
+      int nd2 = fe2->GetDof();
+      int order = fe->GetOrder();
+      int int_order = 2*order + 3;
+      const IntegrationRule *ir = &IntRules.Get(fe->GetGeomType(), int_order);
+      Vector shape(nd2);
+      Array<int> dofs;
+      Vector vals;
+      fespace->GetElementVDofs(i, dofs);
+      x.GetSubVector(dofs, vals);
+      DenseMatrix dshape(nd, dim);
+      Vector elvec(nd2);
+      elvec = 0.0;
+      Vector elvec2(nd2);
+      elvec2 = 0.0;
+      for (int q = 0; q < ir->GetNPoints(); q++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(q);
+         T->SetIntPoint(&ip);
+         double w = ip.weight;
+         fe->CalcDShape(ip, dshape);
+         DenseMatrix grad(dim, dim);
+         grad = 0.0;
+         for (int a = 0; a < nd; a++)
+         {
+            for (int j = 0; j < dim; j++)
+            {
+               for (int k = 0; k < dim; k++)
+               {
+                  // In the vector FE, the a-th node has "dim" components stored consecutively.
+                  // du_j/dr_k
+                  grad(j,k) += vals[a + j*nd] * dshape(a,k);
+               }
+            }
+         }
+
+         DenseMatrix Jac(dim);
+         Jac = T->Jacobian();
+         DenseMatrix invJac(dim);
+         invJac = Jac;
+         invJac.Invert();
+         DenseMatrix grad_phys(dim, dim);
+         Mult(grad, invJac, grad_phys);
+
+         DenseMatrix strain(dim, dim);
+         strain = 0.0;
+         for (int j = 0; j < dim; j++)
+         {
+            for (int k = 0; k < dim; k++)
+            {
+               strain(j,k) = 0.5*(grad_phys(j,k) + grad_phys(k,j));
+            }
+         }
+
+         // Compute the stress tensor: σ = λ (tr ε) I + 2 μ ε.
+         DenseMatrix stress(dim, dim);
+         double trace = strain(0,0) + strain(1,1);
+         stress(0,0) = lambda * trace + 2 * mu * strain(0,0);
+         stress(1,1) = lambda * trace + 2 * mu * strain(1,1);
+         stress(0,1) = 2 * mu * strain(0,1);
+         stress(1,0) = stress(0,1);
+
+         // Compute the von Mises stress (for plane stress/strain in 2D)
+         double sigma_vm = sqrt( 0.5*stress(0,0)*stress(0,0)
+                               + 0.5*stress(1,1)*stress(1,1)
+                     + 0.5*(stress(0,0)-stress(1,1))*(stress(0,0)-stress(1,1))
+                           + 3.0*stress(0,1)*stress(0,1) );
+         // Store the computed stress in the DG GridFunction.
+         fe2->CalcShape(ip, shape);
+         elvec.Add(ip.weight*sigma_vm, shape);
+         elvec2.Add(ip.weight*trace, shape);
+      }
+      fespace_dc->GetElementDofs(i, dofs);
+      stress_gf.AddElementVector(dofs, elvec);
+      trace_strain.AddElementVector(dofs, elvec2);
+      stress_gf2(i) = elvec.Sum();
+   }
+
+   ConstantCoefficient zero(0.0);
+   std::cout << stress_gf2.ComputeL2Error(zero) << " k10l2\n";
+
+   // 14. Save the stress field.
+   {
+      ofstream stress_ofs("stress.gf");
+      stress_gf.Save(stress_ofs);
+   }
+
+   // 15. Visualize the displacement solution and the stress field using GLVis.
+   if (visualization)
+   {
+      osockstream sock(19916, "localhost");
+      sock << "solution\n";
+      mesh->Print(sock);
+      stress_gf.Save(sock);
+      sock.send();
+      sock << "window_title 'Elasticity: Von Mises Stress'\n"
+           << "window_geometry "
+           << 470 << " " << 0 << " " << 400 << " " << 400 << "\n"
+           << "keys jRmclA" << endl;
+   }
+
+   if (visualization)
+   {
+      osockstream sock(19916, "localhost");
+      sock << "solution\n";
+      mesh->Print(sock);
+      stress_gf2.Save(sock);
+      sock.send();
+      sock << "window_title 'Elasticity: Von Mises Stress (element-total)'\n"
+           << "window_geometry "
+           << 940 << " " << 0 << " " << 400 << " " << 400 << "\n"
+           << "keys jRmclA" << endl;
+   }
+
+
+   if (visualization)
+   {
+      osockstream sock(19916, "localhost");
+      sock << "solution\n";
+      mesh->Print(sock);
+      trace_strain.Save(sock);
+      sock.send();
+      sock << "window_title 'Elasticity: Trace strain'\n"
+           << "window_geometry "
+           << 470 << " " << 500 << " " << 400 << " " << 400 << "\n"
+           << "keys jRmclA" << endl;
+   }
+
+   // 13. Free the allocated memory.
+   delete a;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete mesh;
+
+   return 0;
+}
+

fem/CMakeLists.txt

@@ -147,6 +147,7 @@ set(SRCS
   tmop/tmop_pa_w3_c0.cpp
   tmop_tools.cpp
   tmop_amr.cpp
+  tmop_ad_err.cpp
   gslib.cpp
   transfer.cpp
   hyperbolic.cpp
@@ -245,6 +246,7 @@ set(HDRS
   tmop/tmop_pa.hpp
   tmop_tools.hpp
   tmop_amr.hpp
+  tmop_ad_err.hpp
   gslib.hpp
   transfer.hpp
   hyperbolic.hpp

fem/fem.hpp

@@ -38,6 +38,7 @@
 #include "tmop.hpp"
 #include "tmop_tools.hpp"
 #include "tmop_amr.hpp"
+#include "tmop_ad_err.hpp"
 #include "gslib.hpp"
 #include "restriction.hpp"
 #include "quadinterpolator.hpp"

fem/gridfunc.cpp

@@ -2892,6 +2892,37 @@ real_t GridFunction::ComputeL2Error(
    return (error < 0.0) ? -sqrt(-error) : sqrt(error);
 }
 
+real_t GridFunction::ComputeIntegral(const IntegrationRule *irs[],
+   const Array<int> *elems) const
+{
+   real_t sumv = 0.0;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   DenseMatrix vals, exact_vals;
+   Vector loc_errs;
+
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      if (elems != NULL && (*elems)[i] == 0) { continue; }
+      fe = fes->GetFE(i);
+      int intorder = 2*fe->GetOrder() + 3; // <----------
+      const IntegrationRule *ir;
+      IntegrationRules IntRulesGLL(0, Quadrature1D::GaussLobatto);
+      ir = &(IntRulesGLL.Get(fe->GetGeomType(), intorder));
+      T = fes->GetElementTransformation(i);
+      GetVectorValues(*T, *ir, vals);
+      Vector valsv(vals.GetData(), vals.TotalSize());
+      for (int j = 0; j < ir->GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(j);
+         T->SetIntPoint(&ip);
+         sumv += ip.weight * T->Weight() * valsv(j);
+      }
+   }
+
+   return sumv;
+}
+
 real_t GridFunction::ComputeElementGradError(int ielem,
                                              VectorCoefficient *exgrad,
                                              const IntegrationRule *irs[]) const
@@ -2930,6 +2961,17 @@ real_t GridFunction::ComputeElementGradError(int ielem,
    return (error < 0.0) ? -sqrt(-error) : sqrt(error);
 }
 
+void GridFunction::ComputeElementGradErrors(VectorCoefficient *exgrad,
+                                              Vector &error,
+                                             const IntegrationRule *irs[]) const
+{
+   error = 0.0;
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      error(i) = ComputeElementGradError(i, exgrad, irs);
+   }
+}
+
 real_t GridFunction::ComputeGradError(VectorCoefficient *exgrad,
                                       const IntegrationRule *irs[]) const
 {

fem/gridfunc.hpp

@@ -508,6 +508,9 @@ public:
    virtual real_t ComputeElementGradError(int ielem, VectorCoefficient *exgrad,
                                           const IntegrationRule *irs[] = NULL) const;
 
+   virtual real_t ComputeIntegral(const IntegrationRule *irs[] = NULL,
+                                  const Array<int> *elems = NULL) const;
+
    /// Returns ||u_ex - u_h||_L2 for H1 or L2 elements
    /* The @a elems input variable expects a list of markers:
       an elem marker equal to 1 will compute the L2 error on that element
@@ -640,6 +643,11 @@ public:
                                        ) const
    { ComputeElementLpErrors(infinity(), exsol, error, NULL, irs); }
 
+   /// Returns ||grad u_ex - grad u_h||_L2 in element ielem for H1 or L2 elements
+   virtual void ComputeElementGradErrors(VectorCoefficient *exgrad,
+                                           Vector &error,
+                                           const IntegrationRule *irs[] = NULL) const;
+
    /** When given a vector weight, compute the pointwise (scalar) error as the
        dot product of the vector error with the vector weight. Otherwise, the
        scalar error is the l_2 norm of the vector error. */

fem/gslib.cpp

@@ -1180,7 +1180,6 @@ void FindPointsGSLIB::DistributePointInfoToOwningMPIRanks(
 
    // Pack data to send via crystal router
    struct gslib::array *outpt = new gslib::array;
-
    struct out_pt { double rst[3]; uint index, elem, proc, code; };
    struct out_pt *pt;
    array_init(struct out_pt, outpt, points_cnt);

fem/pgridfunc.hpp

@@ -288,6 +288,14 @@ public:
                           pfes->GetComm());
    }
 
+   real_t ComputeIntegral(const IntegrationRule *irs[] = NULL,
+    const Array<int> *elems = NULL) const override
+    {
+        double local_sumv = GridFunction::ComputeIntegral(irs, elems);
+        MPI_Allreduce(MPI_IN_PLACE, &local_sumv, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+        return local_sumv;
+    }
+
    real_t ComputeL2Error(Coefficient &exsol,
                          const IntegrationRule *irs[] = NULL,
                          const Array<int> *elems = NULL) const override

fem/tmop.cpp

@@ -14,12 +14,497 @@
 #include "pgridfunc.hpp"
 #include "tmop_tools.hpp"
 #include "../general/forall.hpp"
+#include "../linalg/dual.hpp"
 
 namespace mfem
 {
 
+/* AD related definitions below ========================================*/
+
+/// MFEM native AD-type for first derivatives
+using AD1Type = internal::dual<real_t, real_t>;
+/// MFEM native AD-type for second derivatives
+using AD2Type = internal::dual<AD1Type, AD1Type>;
+
+/*
+Functions for 2x2 DenseMatrix cast as std::vector<type>, assuming column-major storage
+*/
+template <typename type>
+type fnorm2_2D(const std::vector<type> &u)
+{
+   return u[0]*u[0] + u[1]*u[1] + u[2]*u[2] + u[3]*u[3];
+}
+
+template <typename type>
+type fnorm2_3D(const std::vector<type> &u)
+{
+   return u[0]*u[0] + u[1]*u[1] + u[2]*u[2] + u[3]*u[3] + u[4]*u[4] +
+          u[5]*u[5] + u[6]*u[6] + u[7]*u[7] + u[8]*u[8];
+}
+
+template <typename type>
+type det_2D(const std::vector<type> &u)
+{
+   return u[0]*u[3] - u[1]*u[2];
+}
+
+template <typename type>
+type det_3D(const std::vector<type> &u)
+{
+   return u[0]*(u[4]*u[8] - u[5]*u[7]) -
+          u[1]*(u[3]*u[8] - u[5]*u[6]) +
+          u[2]*(u[3]*u[7] - u[4]*u[6]);
+}
+
+template <typename type>
+void mult_2D(const std::vector<type> &u, const std::vector<type> &M,
+             std::vector<type> &mat)
+{
+   mat.resize(u.size());
+
+   mat[0] = u[0]*M[0] + u[2]*M[1];
+   mat[1] = u[1]*M[0] + u[3]*M[1];
+   mat[2] = u[0]*M[2] + u[2]*M[3];
+   mat[3] = u[1]*M[2] + u[3]*M[3];
+}
+
+template <typename type>
+void mult_aTa_2D(const std::vector<type> &in, std::vector<type> &outm)
+{
+   outm.resize(in.size());
+   outm[0] = in[0]*in[0];
+   outm[1] = in[0]*in[2] + in[1]*in[3];
+   outm[2] = in[0]*in[2] + in[1]*in[3];
+   outm[3] = in[3]*in[3];
+}
+
+template <typename scalartype, typename type>
+void add_2D(const scalartype &scalar, const std::vector<type> &u,
+            const DenseMatrix *M, std::vector<type> &mat)
+{
+   mat.resize(u.size());
+   mat[0] = u[0] + scalar * M->Elem(0,0);
+   mat[1] = u[1] + scalar * M->Elem(1,0);
+   mat[2] = u[2] + scalar * M->Elem(0,1);
+   mat[3] = u[3] + scalar * M->Elem(1,1);
+}
+
+template <typename scalartype, typename type>
+void add_2D(const scalartype &scalar, const std::vector<type> &u,
+            const std::vector<type> &M, std::vector<type> &mat)
+{
+   mat.resize(M.size());
+   mat[0] = u[0] + scalar * M[0];
+   mat[1] = u[1] + scalar * M[1];
+   mat[2] = u[2] + scalar * M[2];
+   mat[3] = u[3] + scalar * M[3];
+}
+
+template <typename type>
+void adjoint_2D(const std::vector<type> &in, std::vector<type> &outm)
+{
+   outm.resize(in.size());
+   outm[0] = in[3];
+   outm[1] = -in[1];
+   outm[2] = -in[2];
+   outm[3] = in[0];
+}
+
+template <typename type>
+void transpose_2D(const std::vector<type> &in, std::vector<type> &outm)
+{
+   outm.resize(in.size());
+   outm[0] = in[0];
+   outm[1] = in[2];
+   outm[2] = in[1];
+   outm[3] = in[3];
+}
+
+template <typename scalartype, typename type>
+void add_3D(const scalartype &scalar, const std::vector<type> &u,
+            const DenseMatrix *M, std::vector<type> &mat)
+{
+   mat.resize(u.size());
+   mat[0] = u[0] + scalar * M->Elem(0,0);
+   mat[1] = u[1] + scalar * M->Elem(1,0);
+   mat[2] = u[2] + scalar * M->Elem(2,0);
+
+   mat[3] = u[3] + scalar * M->Elem(0,1);
+   mat[4] = u[4] + scalar * M->Elem(1,1);
+   mat[5] = u[5] + scalar * M->Elem(2,1);
+
+   mat[6] = u[6] + scalar * M->Elem(0,2);
+   mat[7] = u[7] + scalar * M->Elem(1,2);
+   mat[8] = u[8] + scalar * M->Elem(2,2);
+}
+
+/* Metric definitions */
+
+// W = |T-T'|^2, where T'= |T|*I/sqrt(2).
+template <typename type>
+type mu85_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   auto fnorm = sqrt(fnorm2_2D(T));
+   return T[1]*T[1] + T[2]*T[2] +
+          (T[0] - fnorm/sqrt(2))*(T[0] - fnorm/sqrt(2)) +
+          (T[3] - fnorm/sqrt(2))*(T[3] - fnorm/sqrt(2));
+};
+
+// W = 1/tau |T-I|^2.
+template <typename type>
+type mu98_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   DenseMatrix Id(2,2); Id = 0.0;
+   Id(0,0) = 1; Id(1,1) = 1;
+
+   std::vector<type> Mat;
+   add_2D(real_t{-1.0}, T, &Id, Mat);
+
+   return fnorm2_2D(Mat)/det_2D(T);
+};
+
+// W = 1/(tau^0.5) |T-I|^2.
+template <typename type>
+type mu342_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   DenseMatrix Id(3,3); Id = 0.0;
+   Id(0,0) = 1; Id(1,1) = 1; Id(2,2) = 1;
+
+   std::vector<type> Mat;
+   add_3D(real_t{-1.0}, T, &Id, Mat);
+
+   return fnorm2_3D(Mat)/sqrt(det_3D(T));
+};
+
+// (1/4 alpha) | A - (adj A)^t W^t W / omega |^2
+template <typename type>
+type nu11_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   std::vector<type> A;   // T*W = A
+   std::vector<type> AdjA,AdjAt, WtW, WRK, WRK2;
+
+   mult_2D(T,W,A); // We assume that both A and W are nonsingular.
+
+   auto alpha = det_2D(A);
+   auto omega = det_2D(W);
+   adjoint_2D(A, AdjA);
+   transpose_2D(AdjA, AdjAt);
+
+   mult_aTa_2D(W, WtW);
+   mult_2D(AdjAt, WtW, WRK);
+
+   add_2D(-1.0/omega, A, WRK, WRK2);
+   auto fnorm =  fnorm2_2D(WRK2);
+
+   return 0.25 / (alpha) * fnorm;
+};
+
+// 0.5 * ( sqrt(alpha/omega) - sqrt(omega/alpha) )^2
+template <typename type>
+type nu14_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   std::vector<type> A;   // T*W = A
+   mult_2D(T,W,A);
+
+   auto sqalpha = sqrt(det_2D(A));
+   auto sqomega = sqrt(det_2D(W));
+
+   return 0.5*pow(sqalpha/sqomega - sqomega/sqalpha, 2.0);
+};
+
+// (1/alpha) | A - W |^2
+template <typename type>
+type nu36_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   std::vector<type> A;   // T*W = A
+   std::vector<type> AminusW;  // A-W
+
+   mult_2D(T,W,A);
+   add_2D(-1.0,A,W,AminusW);
+   auto fnorm =  fnorm2_2D(AminusW);
+
+   return 1.0 / (det_2D(A)) * fnorm;
+};
+
+// [ 1.0 - cos( phi_A - phi_W ) ] / (sin phi_A * sin phi_W)
+template <typename type>
+type nu50_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   // We assume that both A and W are nonsingular.
+   std::vector<type> A;
+   mult_2D(T,W,A);
+   auto l1_A = sqrt(A[0]*A[0] + A[1]*A[1]);
+   auto l2_A = sqrt(A[2]*A[2] + A[3]*A[3]);
+   auto prod_A = l1_A*l2_A;
+   auto det_A = A[0]*A[3] - A[1]*A[2];
+   auto sin_A = det_A/prod_A;
+   auto cos_A = (A[0]*A[2] + A[1]*A[3])/prod_A;
+
+   auto l1_W = sqrt(W[0]*W[0] + W[1]*W[1]);
+   auto l2_W = sqrt(W[2]*W[2] + W[3]*W[3]);
+   auto prod_W = l1_W*l2_W;
+   auto det_W = W[0]*W[3] - W[1]*W[2];
+   auto sin_W = det_W/prod_W;
+   auto cos_W = (W[0]*W[2] + W[1]*W[3])/prod_W;
+
+   return (1.0 - cos_A*cos_W - sin_A*sin_W)/(sin_A*sin_W);
+};
+
+// [ 0.5 * (ups_A / ups_W + ups_W / ups_A) - cos(phi_A - phi_W) ] /
+// (sin phi_A * sin phi_W), where ups = l_1 l_2 sin(phi)
+template <typename type>
+type nu51_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   std::vector<type> A;
+   mult_2D(T,W,A);
+   // We assume that both A and W are nonsingular.
+   auto l1_A = sqrt(A[0]*A[0] + A[1]*A[1]);
+   auto l2_A = sqrt(A[2]*A[2] + A[3]*A[3]);
+   auto prod_A = l1_A*l2_A;
+   auto det_A = A[0]*A[3] - A[1]*A[2];
+   auto sin_A = det_A/prod_A;
+   auto cos_A = (A[0]*A[2] + A[1]*A[3])/prod_A;
+   auto ups_A = l1_A*l2_A*sin_A;
+
+   auto l1_W = sqrt(W[0]*W[0] + W[1]*W[1]);
+   auto l2_W = sqrt(W[2]*W[2] + W[3]*W[3]);
+   auto prod_W = l1_W*l2_W;
+   auto det_W = W[0]*W[3] - W[1]*W[2];
+   auto sin_W = det_W/prod_W;
+   auto cos_W = (W[0]*W[2] + W[1]*W[3])/prod_W;
+   auto ups_W = l1_W*l2_W*sin_W;
+
+   return (0.5 * (ups_A / ups_W + ups_W / ups_A) - cos_A*cos_W - sin_A*sin_W) /
+          (sin_A*sin_W);
+};
+
+// (1/2 alpha) | A - (|A|/|W|) W |^2
+template <typename type>
+type nu107_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   std::vector<type> A;   // T*W = A
+   std::vector<type> Mat;  // A-W
+   mult_2D(T,W,A);
+
+   auto alpha = det_2D(A);
+   auto aw = sqrt(fnorm2_2D(A))/sqrt(fnorm2_2D(W));
+
+   add_2D(-aw, A, W, Mat);
+   return (0.5/alpha)*fnorm2_2D(Mat);
+};
+
+// 0.5[ 1.0 - cos( phi_A - phi_W ) ]
+template <typename type>
+type skew2D_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   // We assume that both A and W are nonsingular.
+   std::vector<type> A;
+   mult_2D(T,W,A);
+   auto l1_A = sqrt(A[0]*A[0] + A[1]*A[1]);
+   auto l2_A = sqrt(A[2]*A[2] + A[3]*A[3]);
+   auto prod_A = l1_A*l2_A;
+   auto det_A = A[0]*A[3] - A[1]*A[2];
+   auto sin_A = det_A/prod_A;
+   auto cos_A = (A[0]*A[2] + A[1]*A[3])/prod_A;
+
+   auto l1_W = sqrt(W[0]*W[0] + W[1]*W[1]);
+   auto l2_W = sqrt(W[2]*W[2] + W[3]*W[3]);
+   auto prod_W = l1_W*l2_W;
+   auto det_W = W[0]*W[3] - W[1]*W[2];
+   auto sin_W = det_W/prod_W;
+   auto cos_W = (W[0]*W[2] + W[1]*W[3])/prod_W;
+
+   return 0.5*(1.0 - cos_A*cos_W - sin_A*sin_W);
+};
+
+// 2 - cos( theta_A - theta_W) - cos( theta_A - theta_W + phi_A - phi_W )
+// [2 - cos(theta_A-theta_W) - cos(phi_A - phi_W)]/sin_Asin_W
+template <typename type>
+type nuOQ_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   std::vector<type> A;
+   mult_2D(T,W,A);
+   // We assume that both A and W are nonsingular.
+   auto l1_A = sqrt(A[0]*A[0] + A[1]*A[1]);
+   auto l2_A = sqrt(A[2]*A[2] + A[3]*A[3]);
+   auto prod_A = l1_A*l2_A;
+   auto det_A = A[0]*A[3] - A[1]*A[2];
+   auto sin_A = det_A/prod_A;
+   auto cos_A = (A[0]*A[2] + A[1]*A[3])/prod_A;
+   auto ups_A = l1_A*l2_A*sin_A;
+   auto lenA = sqrt(A[0]*A[0] + A[1]*A[1]);
+   auto cos_theta_A = A[0]/lenA;
+   auto sin_theta_A = A[1]/lenA;
+
+   auto l1_W = sqrt(W[0]*W[0] + W[1]*W[1]);
+   auto l2_W = sqrt(W[2]*W[2] + W[3]*W[3]);
+   auto prod_W = l1_W*l2_W;
+   auto det_W = W[0]*W[3] - W[1]*W[2];
+   auto sin_W = det_W/prod_W;
+   auto cos_W = (W[0]*W[2] + W[1]*W[3])/prod_W;
+   auto ups_W = l1_W*l2_W*sin_W;
+   auto lenW = sqrt(W[0]*W[0] + W[1]*W[1]);
+   auto cos_theta_W = W[0]/lenW;
+   auto sin_theta_W = W[1]/lenW;
+
+   return (2.0 - cos_theta_A*cos_theta_W - sin_theta_A*sin_theta_W -
+          cos_A*cos_W - sin_A*sin_W)/(sin_A*sin_W);
+};
+
+// 2 - cos( theta_A - theta_W) - cos( theta_A - theta_W + phi_A - phi_W )
+// 2 - cos(theta_A)*cos(theta_W) - sin_theta_A*sin_theta_W - [ cos(theta_A)*cos(theta_W)*cos_A*cos_W + cos(theta_A)*cos(theta_W)*sin_A*sin_W + sin_theta_A*sin_theta_W*cos_A*cos_W + sin_theta_A*sin_theta_W*sin_A*sin_W - sin_theta_A*cos(theta_W)*sin_A*cos_W + sin_theta_A*cos(theta_W)*cos_A*sin_W + cos(theta_A)*sin_theta_W*sin_A*cos_W - cos(theta_A)*sin_theta_W*cos_A*sin_W ]
+template <typename type>
+type nuOQ2_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   std::vector<type> A;
+   mult_2D(T,W,A);
+   // We assume that both A and W are nonsingular.
+   auto l1_A = sqrt(A[0]*A[0] + A[1]*A[1]);
+   auto l2_A = sqrt(A[2]*A[2] + A[3]*A[3]);
+   auto prod_A = l1_A*l2_A;
+   auto det_A = A[0]*A[3] - A[1]*A[2];
+   auto sin_A = det_A/prod_A;
+   auto cos_A = (A[0]*A[2] + A[1]*A[3])/prod_A;
+   auto ups_A = l1_A*l2_A*sin_A;
+   auto lenA = sqrt(A[0]*A[0] + A[1]*A[1]);
+   auto cos_theta_A = A[0]/lenA;
+   auto sin_theta_A = A[1]/lenA;
+
+   auto l1_W = sqrt(W[0]*W[0] + W[1]*W[1]);
+   auto l2_W = sqrt(W[2]*W[2] + W[3]*W[3]);
+   auto prod_W = l1_W*l2_W;
+   auto det_W = W[0]*W[3] - W[1]*W[2];
+   auto sin_W = det_W/prod_W;
+   auto cos_W = (W[0]*W[2] + W[1]*W[3])/prod_W;
+   auto ups_W = l1_W*l2_W*sin_W;
+   auto lenW = sqrt(W[0]*W[0] + W[1]*W[1]);
+   auto cos_theta_W = W[0]/lenW;
+   auto sin_theta_W = W[1]/lenW;
+
+   return (2 - cos_theta_A*cos_theta_W - sin_theta_A*sin_theta_W -
+            (cos_theta_A*cos_theta_W*cos_A*cos_W + cos_theta_A*cos_theta_W*sin_A*sin_W + sin_theta_A*sin_theta_W*cos_A*cos_W + sin_theta_A*sin_theta_W*sin_A*sin_W - sin_theta_A*cos_theta_W*sin_A*cos_W + sin_theta_A*cos_theta_W*cos_A*sin_W + cos_theta_A*sin_theta_W*sin_A*cos_W - cos_theta_A*sin_theta_W*cos_A*sin_W))/(sin_A*sin_W);
+};
+
+// Given mu(X,Y), compute dmu/dX or dmu/dY. Y is an optional parameter when
+// computing dmu/dX.
+void ADGrad(std::function<AD1Type(std::vector<AD1Type>&,
+                                  std::vector<AD1Type>&)>mu_ad,
+            DenseMatrix &dmu, //output
+            const DenseMatrix &X, // parameter 1
+            const DenseMatrix *Y = nullptr, //parameter 2
+            const bool dX = true /*derivative with respect to X*/)
+{
+   int matsize = X.TotalSize();
+   std::vector<AD1Type> adX(matsize), adY(matsize);
+
+   for (int i=0; i<matsize; i++) { adX[i] = AD1Type{X.GetData()[i], 0.0}; }
+   if (Y)
+   {
+      for (int i=0; i<matsize; i++) { adY[i] = AD1Type{Y->GetData()[i], 0.0}; }
+   }
+
+   if (dX)
+   {
+      for (int i=0; i<matsize; i++)
+      {
+         adX[i] = AD1Type{X.GetData()[i], 1.0};
+         AD1Type rez = mu_ad(adX, adY);
+         dmu.GetData()[i] = rez.gradient;
+         adX[i] = AD1Type{X.GetData()[i], 0.0};
+      }
+   }
+   else
+   {
+      MFEM_VERIFY(Y, "Y cannot be nullptr when dX = false.");
+      for (int i=0; i<matsize; i++)
+      {
+         adY[i] = AD1Type{Y->GetData()[i], 1.0};
+         AD1Type rez = mu_ad(adX,adY);
+         dmu.GetData()[i] = rez.gradient;
+         adY[i] = AD1Type{Y->GetData()[i], 0.0};
+      }
+   }
+}
+
+// Given mu(X,Y), compute d2mu/dX2, where Y is an optional parameter.
+void ADHessian(std::function<AD2Type(std::vector<AD2Type>&,
+                                     std::vector<AD2Type>&)> mu_ad,
+               DenseTensor &d2mu_dX2,
+               const DenseMatrix &X,
+               const DenseMatrix *Y = nullptr)
+{
+   const int matsize = X.TotalSize();
+
+   //use forward-forward mode
+   std::vector<AD2Type> aduu(matsize), adY(matsize);
+   for (int ii = 0; ii < matsize; ii++)
+   {
+      aduu[ii].value = AD1Type{X.GetData()[ii], 0.0};
+      aduu[ii].gradient = AD1Type{0.0, 0.0};
+   }
+   if (Y)
+   {
+      for (int ii=0; ii<matsize; ii++)
+      {
+         adY[ii].value = AD1Type{Y->GetData()[ii], 0.0};
+         adY[ii].gradient = AD1Type{0.0, 0.0};
+      }
+   }
+
+   for (int ii = 0; ii < matsize; ii++)
+   {
+      aduu[ii].value = AD1Type{X.GetData()[ii], 1.0};
+      for (int jj = 0; jj < (ii + 1); jj++)
+      {
+         aduu[jj].gradient = AD1Type{1.0, 0.0};
+         AD2Type rez = mu_ad(aduu, adY);
+         d2mu_dX2(ii).GetData()[jj] = rez.gradient.gradient;
+         d2mu_dX2(jj).GetData()[ii] = rez.gradient.gradient;
+         aduu[jj].gradient = AD1Type{0.0, 0.0};
+      }
+      aduu[ii].value = AD1Type{X.GetData()[ii], 0.0};
+   }
+   return;
+}
+/* end AD related definitions ========================================*/
+
 // Target-matrix optimization paradigm (TMOP) mesh quality metrics.
 
+void TMOP_QualityMetric::DefaultAssembleH(const DenseTensor &H,
+                                          const DenseMatrix &DS,
+                                          const real_t weight,
+                                          DenseMatrix &A) const
+{
+   const int dof = DS.Height(), dim = DS.Width();
+
+   // The first two go over the rows and cols of dP_dJ where P = dW_dJ.
+   for (int r = 0; r < dim; r++)
+   {
+      for (int c = 0; c < dim; c++)
+      {
+         DenseMatrix Hrc = H(r+c*dim);
+
+         // Compute each entry of d(Prc)_dJ.
+         for (int rr = 0; rr < dim; rr++)
+         {
+            for (int cc = 0; cc < dim; cc++)
+            {
+               const double entry_rr_cc = Hrc(rr, cc);
+
+               for (int i = 0; i < dof; i++)
+               {
+                  for (int j = 0; j < dof; j++)
+                  {
+                     A(i+r*dof, j+rr*dof) +=
+                        weight * DS(i, c) * DS(j, cc) * entry_rr_cc;
+                  }
+               }
+            }
+         }
+      }
+   }
+}
+
 real_t TMOP_Combo_QualityMetric::EvalWMatrixForm(const DenseMatrix &Jpt) const
 {
    real_t metric = 0.;
@@ -752,38 +1237,55 @@ void TMOP_Metric_077::AssembleH(const DenseMatrix &Jpt,
 }
 
 // mu_85 = |T-T'|^2, where T'= |T|*I/sqrt(2)
-real_t TMOP_Metric_085::EvalW(const DenseMatrix &Jpt) const
+real_t TMOP_Metric_085::EvalWMatrixForm(const DenseMatrix &Jpt) const
 {
-   MFEM_VERIFY(Jtr != NULL,
-               "Requires a target Jacobian, use SetTargetJacobian().");
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++) { T[i] = AD1Type{Jpt.GetData()[i], 0.0}; }
+   return mu85_ad(T, W).value;
+}
 
-   DenseMatrix Id(2,2);
-   DenseMatrix Mat(2,2);
-   Mat = Jpt;
+void TMOP_Metric_085::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(mu85_ad<AD1Type>, P, Jpt);
+   return;
+}
 
-   Id(0,0) = 1; Id(0,1) = 0;
-   Id(1,0) = 0; Id(1,1) = 1;
-   Id *= Mat.FNorm()/pow(2,0.5);
-
-   Mat.Add(-1.,Id);
-   return Mat.FNorm2();
+void TMOP_Metric_085::AssembleH(const DenseMatrix &Jpt,
+                                const DenseMatrix &DS,
+                                const real_t weight,
+                                DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(mu85_ad<AD2Type>, H, Jpt);
+   this->DefaultAssembleH(H,DS,weight,A);
 }
 
 // mu_98 = 1/(tau)|T-I|^2
-real_t TMOP_Metric_098::EvalW(const DenseMatrix &Jpt) const
+real_t TMOP_Metric_098::EvalWMatrixForm(const DenseMatrix &Jpt) const
 {
-   MFEM_VERIFY(Jtr != NULL,
-               "Requires a target Jacobian, use SetTargetJacobian().");
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++) { T[i] = AD1Type{Jpt.GetData()[i], 0.0}; }
+   return mu98_ad(T, W).value;
+}
 
-   DenseMatrix Id(2,2);
+void TMOP_Metric_098::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(mu98_ad<AD1Type>, P, Jpt);
+   return;
+}
 
-   Id(0,0) = 1; Id(0,1) = 0;
-   Id(1,0) = 0; Id(1,1) = 1;
-
-   DenseMatrix Mat(2,2);
-   Mat = Jpt;
-   Mat.Add(-1,Id);
-   return Mat.FNorm2()/Jtr->Det();
+void TMOP_Metric_098::AssembleH(const DenseMatrix &Jpt,
+                                const DenseMatrix &DS,
+                                const real_t weight,
+                                DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(mu98_ad<AD2Type>, H, Jpt);
+   this->DefaultAssembleH(H,DS,weight,A);
 }
 
 real_t TMOP_Metric_211::EvalW(const DenseMatrix &Jpt) const
@@ -1437,84 +1939,301 @@ void TMOP_Metric_360::AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
    ie.Assemble_ddI3b(-weight, A.GetData());
 }
 
-real_t TMOP_AMetric_011::EvalW(const DenseMatrix &Jpt) const
+real_t TMOP_AMetric_011::EvalWMatrixForm(const DenseMatrix &Jpt) const
 {
    MFEM_VERIFY(Jtr != NULL,
                "Requires a target Jacobian, use SetTargetJacobian().");
-
-   int dim = Jpt.Size();
-
-   DenseMatrix Jpr(dim, dim);
-   Mult(Jpt, *Jtr, Jpr);
-
-   real_t alpha = Jpr.Det(),
-          omega = Jtr->Det();
-
-   DenseMatrix AdjAt(dim), WtW(dim), WRK(dim), Jtrt(dim);
-   CalcAdjugateTranspose(Jpr, AdjAt);
-   Jtrt.Transpose(*Jtr);
-   MultAAt(Jtrt, WtW);
-   WtW *= 1./omega;
-   Mult(AdjAt, WtW, WRK);
-
-   WRK -= Jpr;
-   WRK *= -1.;
-
-   return (0.25/alpha)*WRK.FNorm2();
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++)
+   {
+      T[i] = AD1Type{Jpt.GetData()[i], 0.0};
+      W[i] = AD1Type{Jtr->GetData()[i], 0.0};
+   }
+   return nu11_ad(T, W).value;
 }
 
-real_t TMOP_AMetric_014a::EvalW(const DenseMatrix &Jpt) const
+void TMOP_AMetric_011::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(nu11_ad<AD1Type>, P, Jpt, Jtr);
+   return;
+}
+
+void TMOP_AMetric_011::EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
+{
+   ADGrad(nu11_ad<AD1Type>, PW, Jpt, Jtr, false);
+   return;
+}
+
+void TMOP_AMetric_011::AssembleH(const DenseMatrix &Jpt,
+                                 const DenseMatrix &DS,
+                                 const real_t weight,
+                                 DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(nu11_ad<AD2Type>, H, Jpt, Jtr);
+   this->DefaultAssembleH(H,DS,weight,A);
+}
+
+real_t TMOP_AMetric_014::EvalWMatrixForm(const DenseMatrix &Jpt) const
 {
    MFEM_VERIFY(Jtr != NULL,
                "Requires a target Jacobian, use SetTargetJacobian().");
-
-   int dim = Jpt.Size();
-
-   DenseMatrix Jpr(dim, dim);
-   Mult(Jpt, *Jtr, Jpr);
-
-   real_t sqalpha = pow(Jpr.Det(), 0.5),
-          sqomega = pow(Jtr->Det(), 0.5);
-
-   return 0.5*pow(sqalpha/sqomega - sqomega/sqalpha, 2.);
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++)
+   {
+      T[i] = AD1Type{Jpt.GetData()[i], 0.0};
+      W[i] = AD1Type{Jtr->GetData()[i], 0.0};
+   }
+   return nu14_ad(T, W).value;
 }
 
-real_t TMOP_AMetric_036::EvalW(const DenseMatrix &Jpt) const
+void TMOP_AMetric_014::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(nu14_ad<AD1Type>, P, Jpt, Jtr);
+   return;
+}
+
+void TMOP_AMetric_014::EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
+{
+   ADGrad(nu14_ad<AD1Type>, PW, Jpt, Jtr, false);
+   return;
+}
+
+void TMOP_AMetric_014::AssembleH(const DenseMatrix &Jpt,
+                                 const DenseMatrix &DS,
+                                 const real_t weight,
+                                 DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(nu14_ad<AD2Type>, H, Jpt, Jtr);
+   this->DefaultAssembleH(H,DS,weight,A);
+}
+
+real_t TMOP_AMetric_036::EvalWMatrixForm(const DenseMatrix &Jpt) const
 {
    MFEM_VERIFY(Jtr != NULL,
                "Requires a target Jacobian, use SetTargetJacobian().");
-
-   int dim = Jpt.Size();
-
-   DenseMatrix Jpr(dim, dim);
-   Mult(Jpt, *Jtr, Jpr); // T*W = A
-
-   real_t alpha = Jpr.Det(); // det(A)
-   Jpr -= *Jtr; // A-W
-
-   return (1./alpha)*(Jpr.FNorm2()); //(1/alpha)*(|A-W|^2)
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++)
+   {
+      T[i] = AD1Type{Jpt.GetData()[i], 0.0};
+      W[i] = AD1Type{Jtr->GetData()[i], 0.0};
+   }
+   return nu36_ad(T, W).value;
 }
 
-real_t TMOP_AMetric_107a::EvalW(const DenseMatrix &Jpt) const
+void TMOP_AMetric_036::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(nu36_ad<AD1Type>, P, Jpt, Jtr);
+   return;
+}
+
+void TMOP_AMetric_036::EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
+{
+   ADGrad(nu36_ad<AD1Type>, PW, Jpt, Jtr, false);
+   return;
+}
+
+void TMOP_AMetric_036::AssembleH(const DenseMatrix &Jpt,
+                                 const DenseMatrix &DS,
+                                 const real_t weight,
+                                 DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(nu36_ad<AD2Type>, H, Jpt, Jtr);
+   this->DefaultAssembleH(H,DS,weight,A);
+}
+
+real_t TMOP_AMetric_050::EvalWMatrixForm(const DenseMatrix &Jpt) const
 {
    MFEM_VERIFY(Jtr != NULL,
                "Requires a target Jacobian, use SetTargetJacobian().");
-
-   int dim = Jpt.Size();
-
-   DenseMatrix Jpr(dim, dim);
-   Mult(Jpt, *Jtr, Jpr);
-
-   real_t alpha = Jpr.Det(),
-          aw = Jpr.FNorm()/Jtr->FNorm();
-
-   DenseMatrix W = *Jtr;
-   W *= aw;
-   Jpr -= W;
-
-   return (0.5/alpha)*Jpr.FNorm2();
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++)
+   {
+      T[i] = AD1Type{Jpt.GetData()[i], 0.0};
+      W[i] = AD1Type{Jtr->GetData()[i], 0.0};
+   }
+   return nu50_ad(T, W).value;
 }
 
+void TMOP_AMetric_050::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(nu50_ad<AD1Type>, P, Jpt, Jtr);
+   return;
+}
+
+void TMOP_AMetric_050::EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
+{
+   ADGrad(nu50_ad<AD1Type>, PW, Jpt, Jtr, false);
+   return;
+}
+
+void TMOP_AMetric_050::AssembleH(const DenseMatrix &Jpt,
+                                 const DenseMatrix &DS,
+                                 const real_t weight,
+                                 DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(nu50_ad<AD2Type>, H, Jpt, Jtr);
+   this->DefaultAssembleH(H,DS,weight,A);
+}
+
+real_t TMOP_AMetric_051::EvalWMatrixForm(const DenseMatrix &Jpt) const
+{
+   MFEM_VERIFY(Jtr != NULL,
+               "Requires a target Jacobian, use SetTargetJacobian().");
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++)
+   {
+      T[i] = AD1Type{Jpt.GetData()[i], 0.0};
+      W[i] = AD1Type{Jtr->GetData()[i], 0.0};
+   }
+   return nu51_ad(T, W).value;
+}
+
+void TMOP_AMetric_051::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(nu51_ad<AD1Type>, P, Jpt, Jtr);
+   return;
+}
+
+void TMOP_AMetric_051::EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
+{
+   ADGrad(nu51_ad<AD1Type>, PW, Jpt, Jtr, false);
+   return;
+}
+
+void TMOP_AMetric_051::AssembleH(const DenseMatrix &Jpt,
+                                 const DenseMatrix &DS,
+                                 const real_t weight,
+                                 DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(nu51_ad<AD2Type>, H, Jpt, Jtr);
+   this->DefaultAssembleH(H,DS,weight,A);
+}
+
+real_t TMOP_AMetric_107::EvalWMatrixForm(const DenseMatrix &Jpt) const
+{
+   MFEM_VERIFY(Jtr != NULL,
+               "Requires a target Jacobian, use SetTargetJacobian().");
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++)
+   {
+      T[i] = AD1Type{Jpt.GetData()[i], 0.0};
+      W[i] = AD1Type{Jtr->GetData()[i], 0.0};
+   }
+   return nu107_ad(T, W).value;
+}
+
+void TMOP_AMetric_107::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(nu107_ad<AD1Type>, P, Jpt, Jtr);
+   return;
+}
+
+void TMOP_AMetric_107::EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
+{
+   ADGrad(nu107_ad<AD1Type>, PW, Jpt, Jtr, false);
+   return;
+}
+
+void TMOP_AMetric_107::AssembleH(const DenseMatrix &Jpt,
+                                 const DenseMatrix &DS,
+                                 const real_t weight,
+                                 DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(nu107_ad<AD2Type>, H, Jpt, Jtr);
+   this->DefaultAssembleH(H,DS,weight,A);
+}
+
+real_t TMOP_AMetric_OQ::EvalWMatrixForm(const DenseMatrix &Jpt) const
+{
+   MFEM_VERIFY(Jtr != NULL,
+               "Requires a target Jacobian, use SetTargetJacobian().");
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++)
+   {
+      T[i] = AD1Type{Jpt.GetData()[i], 0.0};
+      W[i] = AD1Type{Jtr->GetData()[i], 0.0};
+   }
+   return nuOQ_ad(T, W).value;
+}
+
+void TMOP_AMetric_OQ::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(nuOQ_ad<AD1Type>, P, Jpt, Jtr);
+   return;
+}
+
+void TMOP_AMetric_OQ::EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
+{
+   ADGrad(nuOQ_ad<AD1Type>, PW, Jpt, Jtr, false);
+   return;
+}
+
+void TMOP_AMetric_OQ::AssembleH(const DenseMatrix &Jpt,
+                                 const DenseMatrix &DS,
+                                 const real_t weight,
+                                 DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(nuOQ_ad<AD2Type>, H, Jpt, Jtr);
+   this->DefaultAssembleH(H,DS,weight,A);
+}
+
+real_t TMOP_AMetric_OQ2::EvalWMatrixForm(const DenseMatrix &Jpt) const
+{
+   MFEM_VERIFY(Jtr != NULL,
+               "Requires a target Jacobian, use SetTargetJacobian().");
+   int matsize = Jpt.TotalSize();
+   std::vector<AD1Type> T(matsize), W(matsize);
+   for (int i=0; i<matsize; i++)
+   {
+      T[i] = AD1Type{Jpt.GetData()[i], 0.0};
+      W[i] = AD1Type{Jtr->GetData()[i], 0.0};
+   }
+   return nuOQ2_ad(T, W).value;
+}
+
+void TMOP_AMetric_OQ2::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   ADGrad(nuOQ2_ad<AD1Type>, P, Jpt, Jtr);
+   return;
+}
+
+void TMOP_AMetric_OQ2::EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
+{
+   ADGrad(nuOQ2_ad<AD1Type>, PW, Jpt, Jtr, false);
+   return;
+}
+
+void TMOP_AMetric_OQ2::AssembleH(const DenseMatrix &Jpt,
+                                 const DenseMatrix &DS,
+                                 const real_t weight,
+                                 DenseMatrix &A) const
+{
+   const int dim = Jpt.Height();
+   DenseTensor H(dim, dim, dim*dim); H = 0.0;
+   ADHessian(nuOQ2_ad<AD2Type>, H, Jpt, Jtr);
+   this->DefaultAssembleH(H,DS,weight,A);
+}
 
 void TargetConstructor::ComputeAvgVolume() const
 {
@@ -3704,6 +4423,22 @@ void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
 
          d_detW_dx += d_Winv_dx;
          AddMultVWt(shape, d_detW_dx, PMatO);
+
+         // dmu/dW:dW/dx_i
+         DenseMatrix PW(dim);
+         Vector dmudxw(dim);
+         metric->EvalPW(Jpt, PW);
+         DenseMatrix Prod(dim);
+
+         for (int d = 0; d < dim; d++)
+         {
+            const DenseMatrix &dJtr_q = dJtr(q + d*nqp);
+            Prod = 0.0;
+            MultAtB(PW, dJtr_q, Prod); // dmu/dW:dW/dx_i
+            dmudxw(d) = Prod.Trace();
+         }
+         dmudxw *= weight_m;
+         AddMultVWt(shape, dmudxw, PMatO);
       }
 
       if (lim_coeff)

fem/tmop.hpp

@@ -30,6 +30,19 @@ protected:
        for TMOP_QualityMetric%s, because it is not used. */
    void SetTransformation(ElementTransformation &) { }
 
+   /** @brief Default function for assembling the AD computed derivatives into the local gradient matrix 'A'.
+    *
+       @param[in] H       Dense tensor holding the AD computed local gradients.
+       @param[in] DS      Gradient of the basis matrix (dof x dim).
+       @param[in] weight  Quadrature weight coefficient for the point.
+       @param[in,out]  A  Local gradient matrix where the contribution from this
+                          point will be added.
+
+       Computes weight * d(dW_dxi)_d(xj) computed by AD at the current point,
+       for all i and j, where x1 ... xn are the FE dofs. . */
+   void DefaultAssembleH(const DenseTensor &H, const DenseMatrix &DS,
+                         const real_t weight, DenseMatrix &A) const;
+
 public:
    TMOP_QualityMetric() : Jtr(NULL) { }
    virtual ~TMOP_QualityMetric() { }
@@ -59,6 +72,9 @@ public:
        @param[out]  P  The evaluated 1st Piola-Kirchhoff stress tensor. */
    virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const = 0;
 
+   virtual void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
+   { PW = 0.0;}
+
    /** @brief Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor
        and assemble its contribution to the local gradient matrix 'A'.
        @param[in] Jpt     Represents the target->physical transformation
@@ -572,14 +588,14 @@ class TMOP_Metric_085 : public TMOP_QualityMetric
 {
 public:
    // W = |T-T'|^2, where T'= |T|*I/sqrt(2).
-   real_t EvalW(const DenseMatrix &Jpt) const override;
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
 
-   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override
-   { MFEM_ABORT("Not implemented"); }
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
 
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                  const real_t weight, DenseMatrix &A) const override
-   { MFEM_ABORT("Not implemented"); }
+                  const real_t weight, DenseMatrix &A) const override;
 };
 
 /// 2D compound barrier Shape+Size (VS) metric (balanced).
@@ -629,14 +645,14 @@ class TMOP_Metric_098 : public TMOP_QualityMetric
 {
 public:
    // W = 1/tau |T-I|^2.
-   real_t EvalW(const DenseMatrix &Jpt) const override;
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
 
-   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override
-   { MFEM_ABORT("Not implemented"); }
+   virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
 
-   void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                  const real_t weight, DenseMatrix &A) const override
-   { MFEM_ABORT("Not implemented"); }
+   virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
+                          const real_t weight, DenseMatrix &A) const override;
 };
 
 /// 2D untangling metric.
@@ -1096,32 +1112,36 @@ protected:
 
 public:
    // (1/4 alpha) | A - (adj A)^t W^t W / omega |^2
-   real_t EvalW(const DenseMatrix &Jpt) const override;
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
 
-   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override
-   { MFEM_ABORT("Not implemented"); }
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
+
+   virtual void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const override;
 
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                  const real_t weight, DenseMatrix &A) const override
-   { MFEM_ABORT("Not implemented"); }
+                  const real_t weight, DenseMatrix &A) const override;
 };
 
 /// 2D barrier Size (V) metric (polyconvex).
-class TMOP_AMetric_014a : public TMOP_QualityMetric
+class TMOP_AMetric_014 : public TMOP_QualityMetric
 {
 protected:
    mutable InvariantsEvaluator3D<real_t> ie;
 
 public:
    // 0.5 * ( sqrt(alpha/omega) - sqrt(omega/alpha) )^2
-   real_t EvalW(const DenseMatrix &Jpt) const override;
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
 
-   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override
-   { MFEM_ABORT("Not implemented"); }
+   virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
 
-   void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                  const real_t weight, DenseMatrix &A) const override
-   { MFEM_ABORT("Not implemented"); }
+   virtual void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const override;
+
+   virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
+                          const real_t weight, DenseMatrix &A) const override;
 };
 
 /// 2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
@@ -1132,32 +1152,78 @@ protected:
 
 public:
    // (1/alpha) | A - W |^2
-   real_t EvalW(const DenseMatrix &Jpt) const override;
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
 
-   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override
-   { MFEM_ABORT("Not implemented"); }
+   virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
+
+   virtual void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const override;
+
+   virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
+                          const real_t weight, DenseMatrix &A) const override;
+};
+
+/// 2D barrier Skew (Q) metric.
+class TMOP_AMetric_050 : public TMOP_QualityMetric
+{
+protected:
+   mutable InvariantsEvaluator2D<real_t> ie;
+
+public:
+   // [ 1.0 - cos( phi_A - phi_W ) ] / (sin phi_A * sin phi_W)
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
+
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
+
+   void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const override;
 
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                  const real_t weight, DenseMatrix &A) const override
-   { MFEM_ABORT("Not implemented"); }
+                  const real_t weight, DenseMatrix &A) const override;
+};
+
+/// 2D barrier Size+Skew (VQ) metric.
+class TMOP_AMetric_051 : public TMOP_QualityMetric
+{
+protected:
+   mutable InvariantsEvaluator2D<real_t> ie;
+
+public:
+   // [ 0.5 * (ups_A / ups_W + ups_W / ups_A) - cos(phi_A - phi_W) ] /
+   // (sin phi_A * sin phi_W)
+   // where ups = l_1 l_2 sin(phi)
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
+
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
+
+   void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const override;
+
+   void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
+                  const real_t weight, DenseMatrix &A) const override;
 };
 
 /// 2D barrier Shape+Orientation (OS) metric (polyconvex).
-class TMOP_AMetric_107a : public TMOP_QualityMetric
+class TMOP_AMetric_107 : public TMOP_QualityMetric
 {
 protected:
-   mutable InvariantsEvaluator3D<real_t> ie;
+   mutable InvariantsEvaluator2D<real_t> ie;
 
 public:
    // (1/2 alpha) | A - (|A|/|W|) W |^2
-   real_t EvalW(const DenseMatrix &Jpt) const override;
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
 
-   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override
-   { MFEM_ABORT("Not implemented"); }
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
+
+   void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const override;
 
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                  const real_t weight, DenseMatrix &A) const override
-   { MFEM_ABORT("Not implemented"); }
+                  const real_t weight, DenseMatrix &A) const override;
 };
 
 /// 2D barrier Shape+Size (VS) metric (polyconvex).
@@ -1169,7 +1235,7 @@ protected:
 
 public:
    TMOP_AMetric_126(real_t gamma)
-      : sh_metric(new TMOP_AMetric_011), sz_metric(new TMOP_AMetric_014a)
+      : sh_metric(new TMOP_AMetric_011), sz_metric(new TMOP_AMetric_014)
    {
       // (1-gamma) nu_11 + gamma nu_14
       AddQualityMetric(sh_metric, 1.-gamma);
@@ -1179,6 +1245,72 @@ public:
    virtual ~TMOP_AMetric_126() { delete sh_metric; delete sz_metric; }
 };
 
+/// 2D barrier Shape+Skew (SQ) metric.
+/// gamma is recommended to be in (0, 0.9) as a pure skew metric has poor
+/// convergence properties.
+class TMOP_AMetric_049 : public TMOP_Combo_QualityMetric
+{
+protected:
+   mutable InvariantsEvaluator2D<real_t> ie;
+   TMOP_QualityMetric *sh_metric, *sk_metric;
+
+public:
+   TMOP_AMetric_049(real_t gamma)
+      : sh_metric(new TMOP_Metric_002), sk_metric(new TMOP_AMetric_050)
+   {
+      // (1-gamma) mu_2 + gamma nu_50
+      AddQualityMetric(sh_metric, 1.0 - gamma);
+      AddQualityMetric(sk_metric, gamma);
+   }
+
+   int Id() const override { return 49; }
+   real_t GetGamma() const { return wt_arr[1]; }
+
+   virtual ~TMOP_AMetric_049() { delete sh_metric; delete sk_metric; }
+};
+
+/// 2D barrier Skew+Orientation (OQ) metric (polyconvex).
+class TMOP_AMetric_OQ : public TMOP_QualityMetric
+{
+protected:
+   mutable InvariantsEvaluator2D<real_t> ie;
+
+public:
+// 2 - cos( theta_A - theta_W) - cos( theta_A - theta_W + phi_A - phi_W )
+// [2 - cos(theta_A-theta_W) - cos(phi_A - phi_W)]/sin(phi_A)sin(phi_W)
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
+
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
+
+   void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const override;
+
+   void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
+                  const real_t weight, DenseMatrix &A) const override;
+};
+
+/// 2D barrier Skew+Orientation (OQ) metric (polyconvex).
+class TMOP_AMetric_OQ2 : public TMOP_QualityMetric
+{
+protected:
+   mutable InvariantsEvaluator2D<real_t> ie;
+
+public:
+// 2 - cos( theta_A - theta_W) - cos( theta_A - theta_W + phi_A - phi_W )
+// [2 - cos(theta_A-theta_W) - cos(phi_A - phi_W)]/sin(phi_A)sin(phi_W)
+   real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
+   real_t EvalW(const DenseMatrix &Jpt) const override
+   { return EvalWMatrixForm(Jpt); }
+
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
+
+   void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const override;
+
+   void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
+                  const real_t weight, DenseMatrix &A) const override;
+};
+
 /// Base class for limiting functions to be used in class TMOP_Integrator.
 /** This class represents a scalar function f(x, x0, d), where x and x0 are
     positions in physical space, and d is a reference physical distance

fem/tmop_tools.hpp

@@ -14,8 +14,13 @@
 
 #include "bilinearform.hpp"
 #include "pbilinearform.hpp"
+#include "linearform.hpp"
+#include "plinearform.hpp"
 #include "tmop.hpp"
 #include "gslib.hpp"
+#include "../linalg/mma.hpp"
+#include "tmop_ad_err.hpp"
+#include "datacollection.hpp"
 
 namespace mfem
 {
@@ -228,6 +233,9 @@ public:
    /// Update (i) discrete functions at new nodal positions, and
    /// (ii) surface fitting weight.
    void ProcessNewState(const Vector &x) const override;
+   Vector GetProlongedVector(const Vector &x) const;
+   Vector GetRestrictionTransposeVector(const Vector &x) const;
+   Vector GetProlongedTransposeVector(const Vector &x) const;
 
    /** @name Methods for adaptive surface fitting.
        \brief These methods control the behavior of the weight and the
@@ -356,6 +364,99 @@ void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm,
                        char *title, int position);
 #endif
 
+class TMOP_MMA : public MMAOpt, public TMOPNewtonSolver
+{
+protected:
+   // const Operator *oper;
+   // int height; ///< Dimension of the output / number of rows in the matrix.
+   // int width;  ///< Dimension of the input / number of columns in the matrix.
+   // mutable Vector r, c;
+   // /// Limit for the number of iterations the solver is allowed to do
+   // int max_iter;
+   mutable double dlower = 1e-8,
+                  dupper = 1e-8;
+   Vector true_dofs;
+   QuantityOfInterest *qoi = nullptr;
+   Diffusion_Solver *ds = nullptr;
+   Diffusion_Solver *dss = nullptr; //strong diffusion solver
+   VectorHelmholtz *filter = nullptr;
+   double weight = 1.0;
+   real_t ls_norm_fac = 1.2;
+   real_t ls_energy_fac = 1.1;
+   VisItDataCollection *dc = nullptr;
+   VisItDataCollection *dcminl2 = nullptr;
+   VisItDataCollection *dcmingrad = nullptr;
+   ParMesh *pmesh = nullptr;
+   ParGridFunction *l2error_gf = nullptr;
+   ParGridFunction *graderror_gf = nullptr;
+   int ofq = 1;
+   real_t min_l2_err = 100.0;
+   real_t min_grad_err = 100.0;
+   int min_err_iter = -1;
+   int cycle_count = 1;
+
+public:
+   TMOP_MMA(int nVar, int nCon, Vector xval, const IntegrationRule &irule) :
+      MMAOpt(nVar, nCon, xval), TMOPNewtonSolver(irule) {}
+
+#ifdef MFEM_USE_MPI
+   TMOP_MMA(MPI_Comm comm_, int nVar, int nCon, Vector xval,
+            const IntegrationRule &irule) :
+      MMAOpt(comm_, nVar, nCon, xval), TMOPNewtonSolver(comm_, irule) {}
+#endif
+
+   using TMOPNewtonSolver::SetOperator;
+   // void SetOperator(const Operator &op);
+
+   using TMOPNewtonSolver::Mult;
+   void Mult(Vector &x);
+   void MultFilter(Vector &x);
+
+
+   real_t ComputeScalingFactor2(const Vector &x, const Vector &b) const;
+   real_t ComputeScalingFactor2Filter(const Vector &x_orig,
+                                      const Vector &dx,
+                                      const Vector &x_old,
+                                            Vector &fdx) const;
+
+   // void SetMaxIter(int max_it) { max_iter = max_it; }
+
+   void SetLowerBound(const double dlower_) { dlower = dlower_; }
+   void SetUpperBound(const double dupper_) { dupper = dupper_; }
+
+   void SetTrueDofs(Vector &tvec) { true_dofs = tvec; }
+   void SetQuantityOfInterest(QuantityOfInterest *qoi_) { qoi = qoi_; }
+   void SetDiffusionSolver(Diffusion_Solver *ds_) { ds = ds_; }
+   void SetDiffusionSolverStrong(Diffusion_Solver *dss_) { dss = dss_; }
+   void SetVectorHelmholtzFilter(VectorHelmholtz *filter_) { filter = filter_; }
+   void SetQoIWeight(double w) { weight = w; }
+   void SetLineSearchEnergyFactor(double ls_energy_fac_) { ls_energy_fac = ls_energy_fac_;}
+   void SetLineSearchNormFactor(double ls_norm_fac_) { ls_norm_fac = ls_norm_fac_; }
+   real_t GetEnergy(const Vector &x, bool include_qoi = true);
+   void SetDataCollectionObjectandMesh(VisItDataCollection *vdc, ParMesh *pm, int output_freq=1)
+   {
+      dc = vdc;
+      pmesh = pm;
+      ofq = output_freq;
+   }
+   void SetVIsitDataCollectionForMinErrors(VisItDataCollection *vdcminl2,
+                                           VisItDataCollection *vdcmingrad,
+                                           ParGridFunction *l2error_gf_,ParGridFunction *graderror_gf_)
+   {
+      dcminl2 = vdcminl2;
+      dcmingrad = vdcmingrad;
+      l2error_gf = l2error_gf_;
+      graderror_gf = graderror_gf_;
+   }
+   void GetMinErrInfo(double &min_l2, double &min_grad, int &min_iter)
+   {
+      min_l2 = min_l2_err;
+      min_grad = min_grad_err;
+      min_iter = min_err_iter;
+   }
+   int GetVisCount() { return cycle_count; }
+};
+
 }
 
 #endif

linalg/CMakeLists.txt

@@ -26,6 +26,7 @@ list(APPEND SRCS
   symmat.cpp
   handle.cpp
   matrix.cpp
+  mma.cpp
   ode.cpp
   operator.cpp
   solvers.cpp
@@ -57,6 +58,7 @@ list(APPEND HDRS
   lapack.hpp
   linalg.hpp
   matrix.hpp
+  mma.hpp
   ode.hpp
   operator.hpp
   solvers.hpp

linalg/densemat.hpp

@@ -100,6 +100,7 @@ public:
 
    /// For backward compatibility define Size to be synonym of Width()
    int Size() const { return Width(); }
+   int TotalSize() const { return Width()*Height(); }
 
    /// Change the size of the DenseMatrix to s x s.
    void SetSize(int s) { SetSize(s, s); }

linalg/dual.hpp

@@ -345,6 +345,15 @@ dual<value_type, gradient_type> tan(dual<value_type, gradient_type> a)
    return {f, a.gradient * (value_type{1} + f * f)};
 }
 
+/** @brief implementation of tanh for dual numbers */
+template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
+dual<value_type, gradient_type> tanh(dual<value_type, gradient_type> a)
+{
+   using std::tanh;
+   value_type f = tanh(a.value);
+   return {f, a.gradient * (value_type{1} - f * f)};
+}
+
 /** @brief implementation of atan for dual numbers */
 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
 dual<value_type, gradient_type> atan(dual<value_type, gradient_type> a)

linalg/linalg.hpp

@@ -32,6 +32,7 @@
 #include "invariants.hpp"
 #include "constraints.hpp"
 #include "auxiliary.hpp"
+#include "mma.hpp"
 #include "batched/batched.hpp"
 #include "batched/gpu_blas.hpp"
 #include "batched/solver.hpp"

linalg/mma.hpp

@@ -0,0 +1,272 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_MMA
+#define MFEM_MMA
+
+#include <iostream>
+#include "../config/config.hpp"
+#include "vector.hpp"
+#include "densemat.hpp"
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#endif
+
+namespace mfem
+{
+
+class MMA
+{
+public:
+
+   /// Serial constructor:
+   /// nVar - number of design parameters;
+   /// nCon - number of constraints;
+   /// xval[nVar] - initial parameter values
+   MMA(int nVar, int nCon, real_t *xval);
+
+#ifdef MFEM_USE_MPI
+   /// Parallel constructor:
+   /// nVar - number of design parameters;
+   /// nCon - number of constraints;
+   /// xval[nVar] - initial parameter values
+   MMA(MPI_Comm comm_, int nVar, int nCon, real_t *xval);
+#endif
+
+   /// Destructor
+   ~MMA();
+
+   /// Update the optimization parameters
+   /// iter - current iteration number
+   /// dfdx[nVar] - gradients of the objective
+   /// gx[nCon] - values of the constraints
+   /// dgdx[nCon*nVar] - gradients of the constraints
+   /// xxmin[nVar] - lower bounds
+   /// xxmax[nVar] - upper bounds
+   /// xval[nVar] - (input: current optimization parameters)
+   /// xval[nVar] - (output: updated optimization parameters)
+   void Update(int iter, const real_t* dfdx, const real_t* gx,
+               const real_t* dgdx,
+               const real_t* xxmin,const real_t* xxmax,
+               real_t* xval);
+
+
+   /// Dump the internal state into a file
+   /// xval[nVar] - current optimization parameters
+   /// iter - current interation
+   /// fname, fpath - file name and file path
+   void WriteState(real_t* xval, int iter,
+                   std::string fname,
+                   std::string fpath = "./");
+
+   /// Load the internal state
+   void LoadState(real_t* xval, int iter,
+                  std::string fname,
+                  std::string fpath = "./");
+
+protected:
+   // Local vectors
+   real_t *a, *b, *c, *d;
+   real_t a0, machineEpsilon, epsimin;
+   real_t z, zet;
+   int nCon, nVar;
+
+   // Global: Asymptotes, bounds, objective approx., constraint approx.
+   real_t *low, *upp;
+   real_t *x, *y, *xsi, *eta, *lam, *mu, *s;
+
+
+private:
+
+   // MMA-specific
+   real_t asyinit, asyincr, asydecr;
+   real_t xmamieps, lowmin, lowmax, uppmin, uppmax, zz;
+   real_t *factor;
+
+   /// values from the previous two iterations
+   real_t *xo1, *xo2;
+
+   /// KKT norm
+   real_t kktnorm;
+
+   /// intialization state
+   bool isInitialized = false;
+
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm;
+#endif
+
+   /// Allocate the memory for MMA
+   void AllocData(int nVar, int nCon);
+
+   /// Free the memory for MMA
+   void FreeData();
+
+   /// Initialize data
+   void  InitData(real_t *xval);
+
+   /// Subproblem base class
+   class MMASubBase
+   {
+   public:
+      /// Constructor
+      MMASubBase(MMA* mma) {mma_ptr=mma;}
+
+      /// Destructor
+      virtual ~MMASubBase() {}
+
+      /// Update the optimization parameters
+      virtual
+      void Update(const real_t* dfdx,
+                  const real_t* gx,
+                  const real_t* dgdx,
+                  const real_t* xmin,
+                  const real_t* xmax,
+                  const real_t* xval)=0;
+
+   protected:
+      MMA* mma_ptr;
+
+   };
+
+   MMASubBase* mSubProblem;
+
+   friend class MMASubParallel;
+
+   class MMASubParallel:public MMASubBase
+   {
+   public:
+
+      /// Constructor
+      MMASubParallel(MMA* mma, int nVar, int nCon):MMASubBase(mma)
+      {
+         AllocSubData(nVar,nCon);
+
+         nVar_global = nVar;
+#ifdef MFEM_USE_MPI
+         MPI_Allreduce(&nVar, &nVar_global, 1, MPI_INT, MPI_SUM, mma->comm);
+#endif
+      }
+
+      /// Destructor
+      virtual
+      ~MMASubParallel()
+      {
+         FreeSubData();
+      }
+
+      /// Update the optimization parameters
+      virtual
+      void Update(const real_t* dfdx,
+                  const real_t* gx,
+                  const real_t* dgdx,
+                  const real_t* xmin,
+                  const real_t* xmax,
+                  const real_t* xval);
+
+   private:
+      int ittt, itto, itera, nVar_global;
+
+      real_t epsi, delz, dz, dzet, stmxx, stmalfa, stmbeta,
+             sum, stminv, steg, zold, zetold,
+             residunorm, residumax, resinew, raa0, albefa, move, xmamieps;
+
+      real_t *sum1, *ux1, *xl1, *plam, *qlam, *gvec, *residu, *GG, *delx, *dely,
+             *dellam,
+             *dellamyi, *diagx, *diagy, *diaglamyi, *bb, *bb1, *Alam, *AA, *AA1,
+             *dlam, *dx, *dy, *dxsi, *deta, *dmu, *Axx, *axz, *ds, *xx, *dxx, *stepxx,
+             *stepalfa, *stepbeta, *xold, *yold,
+             *lamold, *xsiold, *etaold, *muold, *sold, *p0, *q0, *P, *Q, *alfa, *beta,
+             *xmami, *b;
+
+      // parallel helper variables
+      real_t global_max = 0.0;
+      real_t global_norm = 0.0;
+      real_t stmxx_global = 0.0;
+      real_t stmalfa_global = 0.0;
+      real_t stmbeta_global = 0.0;
+
+      real_t *b_local, *gvec_local, *Alam_local, *sum_local, *sum_global;
+
+      /// Allocate the memory for the subproblem
+      void AllocSubData(int nVar, int nCon);
+
+      /// Free the memeory for the subproblem
+      void FreeSubData();
+   };
+};
+
+/// Native MFEM MMA interface
+class MMAOpt
+{
+public:
+   /// Default constructor
+   MMAOpt(int nVar, int nCon, mfem::Vector& xval)
+   {
+      opt=new mfem::MMA(nVar,nCon, xval.GetData());
+   }
+
+#ifdef MFEM_USE_MPI
+   MMAOpt(MPI_Comm comm_, int nVar, int nCon, mfem::Vector& xval)
+   {
+      int rank = 0;
+      MPI_Comm_rank(comm_, &rank);
+
+      // create new communicator
+      int colour;
+
+      if ( 0 != nVar)
+      {
+         colour = 0;
+      }
+      else
+      {
+         colour = MPI_UNDEFINED;
+      }
+
+      // Split the global communicator
+      MPI_Comm_split(comm_, colour, rank, &new_comm);
+
+      opt=new mfem::MMA(new_comm, nVar,nCon, xval.GetData());
+   }
+#endif
+
+   /// Destructor
+   ~MMAOpt()
+   {
+      delete opt;
+   }
+
+   /// Design update
+   void Update(int iter, const mfem::Vector& dfdx,
+               const mfem::Vector& gx, const mfem::Vector& dgdx,
+               const mfem::Vector& xmin, const mfem::Vector& xmax,
+               mfem::Vector& xval)
+   {
+      opt->Update(iter, dfdx.GetData(),
+                  gx.GetData(),dgdx.GetData(),
+                  xmin.GetData(), xmax.GetData(),
+                  xval.GetData());
+   }
+
+private:
+   mfem::MMA* opt;// the actual mma optimizer
+
+#ifdef MFEM_USE_MPI
+   MPI_Comm new_comm;
+#endif
+
+};
+
+} // mfem namespace
+
+#endif // MFEM_MMA

miniapps/meshing/CMakeLists.txt

@@ -111,12 +111,33 @@ endif()
 
 # Parallel apps.
 if (MFEM_USE_MPI)
+
+  list(APPEND PARMESHING_COMMON_SOURCES
+    mesh-optimizer_using_NLP.cpp
+    MMA.cpp)
+  list(APPEND PARMESHING_COMMON_HEADERS
+    mesh-optimizer_using_NLP.hpp
+    MMA.hpp)
+
+  convert_filenames_to_full_paths(PARMESHING_COMMON_SOURCES)
+  convert_filenames_to_full_paths(PARMESHING_COMMON_HEADERS)
+
+  set(PARMESHING_COMMON_FILES
+	  EXTRA_SOURCES ${PARMESHING_COMMON_SOURCES} 
+	  EXTRA_HEADERS ${PARMESHING_COMMON_HEADERS})
+
   add_mfem_miniapp(pmesh-optimizer
     MAIN pmesh-optimizer.cpp
     ${MFEM_MINIAPPS_COMMON_HEADERS}
     LIBRARIES mfem-common)
   add_dependencies(pmesh-optimizer copy_miniapps_meshing_data)
 
+  add_mfem_miniapp(pmesh-optimizer_NLP
+    MAIN pmesh-optimizer_NLP.cpp
+    ${PARMESHING_COMMON_FILES}
+    LIBRARIES mfem mfem-common)
+  add_dependencies(pmesh-optimizer_NLP copy_miniapps_meshing_data)
+
   add_mfem_miniapp(pmesh-fitting
     MAIN pmesh-fitting.cpp
     ${MFEM_MINIAPPS_COMMON_HEADERS}

miniapps/meshing/MMA.cpp

@@ -0,0 +1,995 @@
+
+#include "MMA.hpp"
+#include <iostream>
+#include <math.h>
+
+#ifdef MFEM_USE_PETSC
+
+/* -----------------------------------------------------------------------------
+Authors: Niels Aage
+ Copyright (C) 2013-2019,
+This MMA implementation is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+This Module is distributed in the hope that it will be useful,implementation
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Lesser General Public License for more details.
+You should have received a copy of the GNU Lesser General Public
+License along with this Module; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+-------------------------------------------------------------------------- */
+
+MMA::MMA(PetscInt nn, PetscInt mm, PetscInt kk, Vec xo1t, Vec xo2t, Vec Ut, Vec Lt, PetscScalar* at, PetscScalar* ct,
+         PetscScalar* dt) {
+    n = nn;
+    m = mm;
+    k = kk;
+    if (k < 3) {
+        PetscPrintf(PETSC_COMM_WORLD, "NOT A LEGAL RESTART POINT (k<3): EXPECT BREAKDOWN\n");
+    }
+
+    asyminit = 0.5;
+    asymdec  = 0.7;
+    asyminc  = 1.2;
+
+    NonLinConstraints      = PETSC_TRUE;
+    constraintModification = PETSC_FALSE;
+    RobustAsymptotesType   = 0;
+
+    a = new PetscScalar[m];
+    c = new PetscScalar[m];
+    d = new PetscScalar[m];
+
+    memcpy(a, at, mm * sizeof(PetscScalar));
+    memcpy(c, ct, mm * sizeof(PetscScalar));
+    memcpy(d, dt, mm * sizeof(PetscScalar));
+
+    y   = new PetscScalar[m];
+    lam = new PetscScalar[m];
+
+    VecDuplicate(xo1t, &L);
+    VecDuplicate(xo1t, &U);
+
+    VecDuplicate(xo1t, &alpha);
+    VecDuplicate(xo1t, &beta);
+
+    VecDuplicate(xo1t, &p0);
+    VecDuplicate(xo1t, &q0);
+    VecDuplicateVecs(xo1t, m, &pij);
+    VecDuplicateVecs(xo1t, m, &qij);
+
+    b = new PetscScalar[m];
+
+    VecDuplicate(xo1t, &xo1);
+    VecDuplicate(xo1t, &xo2);
+
+    grad = new PetscScalar[m];
+    mu   = new PetscScalar[m];
+    s    = new PetscScalar[2 * m];
+
+    Hess = new PetscScalar[m * m];
+
+    // Now insert the values into xo1,xo2,U,L
+    PetscInt nloc;
+    VecGetLocalSize(xo1t, &nloc);
+
+    // input
+    PetscScalar *pxo1t, *pxo2t, *pUt, *pLt;
+    VecGetArray(xo1t, &pxo1t);
+    VecGetArray(xo2t, &pxo2t);
+    VecGetArray(Ut, &pUt);
+    VecGetArray(Lt, &pLt);
+
+    // internal
+    PetscScalar *pxo1, *pxo2, *pU, *pL;
+    VecGetArray(xo1, &pxo1);
+    VecGetArray(xo2, &pxo2);
+    VecGetArray(U, &pU);
+    VecGetArray(L, &pL);
+
+    // Copy data
+    memcpy(pxo1, pxo1t, nloc * sizeof(PetscScalar));
+    memcpy(pxo2, pxo2t, nloc * sizeof(PetscScalar));
+    memcpy(pU, pUt, nloc * sizeof(PetscScalar));
+    memcpy(pL, pLt, nloc * sizeof(PetscScalar));
+
+    // Restore arrays
+    VecRestoreArray(xo1t, &pxo1t);
+    VecRestoreArray(xo2t, &pxo2t);
+    VecRestoreArray(Ut, &pUt);
+    VecRestoreArray(Lt, &pLt);
+
+    VecRestoreArray(xo1, &pxo1);
+    VecRestoreArray(xo2, &pxo2);
+    VecRestoreArray(U, &pU);
+    VecRestoreArray(L, &pL);
+}
+
+MMA::MMA(PetscInt nn, PetscInt mm, PetscInt kk, Vec xo1t, Vec xo2t, Vec Ut, Vec Lt) {
+    n = nn;
+    m = mm;
+    k = kk;
+    if (k < 3) {
+        PetscPrintf(PETSC_COMM_WORLD, "NOT A LEGAL RESTART POINT (k<3): EXPECT BREAKDOWN\n");
+    }
+
+    asyminit = 0.5;
+    asymdec  = 0.7;
+    asyminc  = 1.2;
+
+    NonLinConstraints      = PETSC_TRUE;
+    constraintModification = PETSC_FALSE;
+    RobustAsymptotesType   = 0;
+
+    a = new PetscScalar[m];
+    c = new PetscScalar[m];
+    d = new PetscScalar[m];
+
+    for (PetscInt i = 0; i < m; i++) {
+        a[i] = 0.0;
+        c[i] = 1000.0;
+        d[i] = 0.0;
+    }
+
+    y   = new PetscScalar[m];
+    lam = new PetscScalar[m];
+
+    VecDuplicate(xo1t, &L);
+    VecDuplicate(xo1t, &U);
+
+    VecDuplicate(xo1t, &alpha);
+    VecDuplicate(xo1t, &beta);
+
+    VecDuplicate(xo1t, &p0);
+    VecDuplicate(xo1t, &q0);
+    VecDuplicateVecs(xo1t, m, &pij);
+    VecDuplicateVecs(xo1t, m, &qij);
+
+    b = new PetscScalar[m];
+
+    VecDuplicate(xo1t, &xo1);
+    VecDuplicate(xo1t, &xo2);
+
+    grad = new PetscScalar[m];
+    mu   = new PetscScalar[m];
+    s    = new PetscScalar[2 * m];
+
+    Hess = new PetscScalar[m * m];
+
+    // Now insert the values into xo1,xo2,U,L
+    PetscInt nloc;
+    VecGetLocalSize(xo1t, &nloc);
+
+    // input
+    PetscScalar *pxo1t, *pxo2t, *pUt, *pLt;
+    VecGetArray(xo1t, &pxo1t);
+    VecGetArray(xo2t, &pxo2t);
+    VecGetArray(Ut, &pUt);
+    VecGetArray(Lt, &pLt);
+
+    // internal
+    PetscScalar *pxo1, *pxo2, *pU, *pL;
+    VecGetArray(xo1, &pxo1);
+    VecGetArray(xo2, &pxo2);
+    VecGetArray(U, &pU);
+    VecGetArray(L, &pL);
+
+    // Copy data
+    memcpy(pxo1, pxo1t, nloc * sizeof(PetscScalar));
+    memcpy(pxo2, pxo2t, nloc * sizeof(PetscScalar));
+    memcpy(pU, pUt, nloc * sizeof(PetscScalar));
+    memcpy(pL, pLt, nloc * sizeof(PetscScalar));
+
+    // Restore arrays
+    VecRestoreArray(xo1t, &pxo1t);
+    VecRestoreArray(xo2t, &pxo2t);
+    VecRestoreArray(Ut, &pUt);
+    VecRestoreArray(Lt, &pLt);
+
+    VecRestoreArray(xo1, &pxo1);
+    VecRestoreArray(xo2, &pxo2);
+    VecRestoreArray(U, &pU);
+    VecRestoreArray(L, &pL);
+}
+
+MMA::MMA(PetscInt nn, PetscInt mm, Vec x, PetscScalar* at, PetscScalar* ct, PetscScalar* dt) {
+    n = nn;
+    m = mm;
+
+    asyminit = 0.5;
+    asymdec  = 0.7;
+    asyminc  = 1.2;
+
+    NonLinConstraints      = PETSC_TRUE;
+    constraintModification = PETSC_FALSE;
+    RobustAsymptotesType   = 0;
+
+    k = 0;
+
+    a = new PetscScalar[m];
+    c = new PetscScalar[m];
+    d = new PetscScalar[m];
+
+    memcpy(a, at, mm * sizeof(PetscScalar));
+    memcpy(c, ct, mm * sizeof(PetscScalar));
+    memcpy(d, dt, mm * sizeof(PetscScalar));
+
+    y   = new PetscScalar[m];
+    lam = new PetscScalar[m];
+
+    VecDuplicate(x, &L);
+    VecDuplicate(x, &U);
+
+    VecDuplicate(x, &alpha);
+    VecDuplicate(x, &beta);
+
+    VecDuplicate(x, &p0);
+    VecDuplicate(x, &q0);
+    VecDuplicateVecs(x, m, &pij);
+    VecDuplicateVecs(x, m, &qij);
+
+    b = new PetscScalar[m];
+
+    VecDuplicate(x, &xo1);
+    VecDuplicate(x, &xo2);
+
+    grad = new PetscScalar[m];
+    mu   = new PetscScalar[m];
+    s    = new PetscScalar[2 * m];
+
+    Hess = new PetscScalar[m * m];
+}
+
+MMA::MMA(PetscInt nn, PetscInt mm, Vec x) {
+    n = nn;
+    m = mm;
+
+    asyminit = 0.5;
+    asymdec  = 0.7;
+    asyminc  = 1.2;
+
+    NonLinConstraints      = PETSC_TRUE;
+    constraintModification = PETSC_FALSE;
+
+    RobustAsymptotesType = 0;
+
+    k = 0;
+
+    a = new PetscScalar[m];
+    c = new PetscScalar[m];
+    d = new PetscScalar[m];
+
+    for (PetscInt i = 0; i < m; i++) {
+        a[i] = 0.0;
+        c[i] = 1000.0;
+        d[i] = 0.0;
+    }
+
+    y   = new PetscScalar[m];
+    lam = new PetscScalar[m];
+
+    VecDuplicate(x, &L);
+    VecDuplicate(x, &U);
+
+    VecDuplicate(x, &alpha);
+    VecDuplicate(x, &beta);
+
+    VecDuplicate(x, &p0);
+    VecDuplicate(x, &q0);
+    VecDuplicateVecs(x, m, &pij);
+    VecDuplicateVecs(x, m, &qij);
+    b = new PetscScalar[m];
+
+    VecDuplicate(x, &xo1);
+    VecDuplicate(x, &xo2);
+
+    grad = new PetscScalar[m];
+    mu   = new PetscScalar[m];
+    s    = new PetscScalar[2 * m];
+
+    Hess = new PetscScalar[m * m];
+}
+
+MMA::~MMA() {
+
+    delete[] a;
+    delete[] b;
+    delete[] c;
+    delete[] d;
+    delete[] y;
+    delete[] lam;
+    VecDestroy(&L);
+    VecDestroy(&U);
+    VecDestroy(&alpha);
+    VecDestroy(&beta);
+    VecDestroy(&p0);
+    VecDestroy(&q0);
+    VecDestroyVecs(m, &pij);
+    VecDestroyVecs(m, &qij);
+    VecDestroy(&xo1);
+    VecDestroy(&xo2);
+    delete[] grad;
+    delete[] mu;
+    delete[] s;
+    delete[] Hess;
+}
+
+// restart method
+
+PetscErrorCode MMA::Restart(Vec xo1t, Vec xo2t, Vec Ut, Vec Lt) {
+
+    PetscErrorCode ierr = 0;
+
+    // Insert values into xo1t,xo2t,Ut,Lt
+    PetscInt nloc;
+    VecGetLocalSize(xo1t, &nloc);
+
+    // input
+    PetscScalar *pxo1t, *pxo2t, *pUt, *pLt;
+    VecGetArray(xo1t, &pxo1t);
+    VecGetArray(xo2t, &pxo2t);
+    VecGetArray(Ut, &pUt);
+    VecGetArray(Lt, &pLt);
+
+    // internal
+    PetscScalar *pxo1, *pxo2, *pU, *pL;
+    VecGetArray(xo1, &pxo1);
+    VecGetArray(xo2, &pxo2);
+    VecGetArray(U, &pU);
+    VecGetArray(L, &pL);
+
+    // Copy data
+    memcpy(pxo1t, pxo1, nloc * sizeof(PetscScalar));
+    memcpy(pxo2t, pxo2, nloc * sizeof(PetscScalar));
+    memcpy(pUt, pU, nloc * sizeof(PetscScalar));
+    memcpy(pLt, pL, nloc * sizeof(PetscScalar));
+
+    // Restore arrays
+    VecRestoreArray(xo1t, &pxo1t);
+    VecRestoreArray(xo2t, &pxo2t);
+    VecRestoreArray(Ut, &pUt);
+    VecRestoreArray(Lt, &pLt);
+
+    VecRestoreArray(xo1, &pxo1);
+    VecRestoreArray(xo2, &pxo2);
+    VecRestoreArray(U, &pU);
+    VecRestoreArray(L, &pL);
+
+    return (ierr);
+}
+
+// Set the aggresivity of the moving asymptotes
+PetscErrorCode MMA::SetAsymptotes(PetscScalar init, PetscScalar decrease, PetscScalar increase) {
+    PetscErrorCode ierr = 0;
+
+    // asymptotes initialization and increase/decrease
+    asyminit = init;
+    asymdec  = decrease;
+    asyminc  = increase;
+    return ierr;
+}
+
+PetscErrorCode MMA::SetRobustAsymptotesType(PetscInt val) {
+
+    PetscErrorCode ierr = 0;
+
+    RobustAsymptotesType = val;
+
+    if (RobustAsymptotesType == 0 || RobustAsymptotesType == 1) {
+    } else {
+        RobustAsymptotesType = 0;
+        PetscPrintf(PETSC_COMM_WORLD, "ERROR in MMA.cc/h: RobustAsymptotesType cannot be set to: %d \n", val);
+    }
+    return ierr;
+}
+
+PetscErrorCode MMA::SetOuterMovelimit(PetscScalar Xmin, PetscScalar Xmax, PetscScalar movlim, Vec x, Vec xmin,
+                                      Vec xmax) {
+
+    PetscErrorCode ierr = 0;
+
+    PetscScalar *xv, *xmiv, *xmav;
+    PetscInt     nloc;
+    VecGetLocalSize(x, &nloc);
+    VecGetArray(x, &xv);
+    VecGetArray(xmin, &xmiv);
+    VecGetArray(xmax, &xmav);
+    for (PetscInt i = 0; i < nloc; i++) {
+        xmav[i] = Min(Xmax, xv[i] + movlim);
+        xmiv[i] = Max(Xmin, xv[i] - movlim);
+    }
+    VecRestoreArray(x, &xv);
+    VecRestoreArray(xmin, &xmiv);
+    VecRestoreArray(xmax, &xmav);
+    return ierr;
+}
+
+PetscScalar MMA::DesignChange(Vec x, Vec xold) {
+
+    PetscScalar *xv, *xo;
+    PetscInt     nloc;
+    VecGetLocalSize(x, &nloc);
+    VecGetArray(x, &xv);
+    VecGetArray(xold, &xo);
+    PetscScalar ch = 0.0;
+    for (PetscInt i = 0; i < nloc; i++) {
+        ch    = PetscMax(ch, PetscAbsReal(xv[i] - xo[i]));
+        xo[i] = xv[i];
+    }
+    PetscScalar tmp;
+    MPI_Allreduce(&ch, &tmp, 1, MPIU_SCALAR, MPI_MAX, PETSC_COMM_WORLD);
+    ch = tmp;
+    VecRestoreArray(x, &xv);
+    VecRestoreArray(xold, &xo);
+
+    return (ch);
+}
+
+PetscErrorCode MMA::KKTresidual(Vec x, Vec dfdx, PetscScalar* fx, Vec* dgdx, Vec xmin, Vec xmax, PetscScalar* norm2,
+                                PetscScalar* normInf) {
+    PetscErrorCode ierr = 0;
+
+    if (!NonLinConstraints) {
+        PetscErrorPrintf("MMA->KKTresidual called WITH constraints but object was "
+                         "allocated WITHOUT !\n");
+        return -1;
+    }
+
+    PetscScalar *xp, *xminp, *xmaxp, *df0dxp, **dfdxp;
+    PetscInt     locsiz;
+    VecGetLocalSize(x, &locsiz);
+
+    VecGetArray(x, &xp);
+    VecGetArray(xmin, &xminp);
+    VecGetArray(xmax, &xmaxp);
+    VecGetArray(dfdx, &df0dxp);
+    VecGetArrays(dgdx, m, &dfdxp);
+
+    PetscScalar resi, ri, mu_min, mu_max;
+
+    norm2[0]   = 0;
+    normInf[0] = 0;
+    for (PetscInt i = 0; i < locsiz; i++) {
+        ri = df0dxp[i];
+        for (PetscInt j = 0; j < m; j++) {
+            ri += lam[j] * dfdxp[j][i];
+        }
+        mu_min = 0.0;
+        if (xp[i] < xminp[i] + 1.0e-5 && ri > 0.0) {
+            mu_min = ri;
+        }
+        mu_max = 0.0;
+        if (xp[i] > xmaxp[i] - 1.0e-5 && ri < 0.0) {
+            mu_max = -ri;
+        }
+        ri += -mu_min + mu_max;
+        norm2[0] += pow(ri, 2.0);
+        normInf[0] = Max(Abs(ri), normInf[0]);
+        resi       = mu_min * (xp[i] - xminp[i]);
+        norm2[0] += pow(resi, 2.0);
+        normInf[0] = Max(Abs(resi), normInf[0]);
+        resi       = mu_max * (xmaxp[i] - xp[i]);
+        norm2[0] += pow(resi, 2.0);
+        normInf[0] = Max(Abs(resi), normInf[0]);
+    }
+
+    VecRestoreArray(x, &xp);
+    VecRestoreArray(xmin, &xminp);
+    VecRestoreArray(xmax, &xmaxp);
+    VecRestoreArray(dfdx, &df0dxp);
+    VecRestoreArrays(dgdx, m, &dfdxp);
+    PetscScalar n2tmp = norm2[0];
+    PetscScalar nItmp = normInf[0];
+    norm2[0]          = 0.0;
+    normInf[0]        = 0.0;
+    MPI_Allreduce(&n2tmp, norm2, 1, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+    MPI_Allreduce(&nItmp, normInf, 1, MPIU_SCALAR, MPI_MAX, PETSC_COMM_WORLD);
+    ri = 0.0;
+    for (PetscInt j = 0; j < m; j++) {
+        ri += lam[j] * (a[j] * z + y[j] - fx[j]);
+    }
+    norm2[0] += pow(ri, 2.0);
+    normInf[0] = Max(Abs(ri), normInf[0]);
+    norm2[0]   = sqrt(norm2[0]);
+
+    return ierr;
+}
+
+// Set and solve a subproblem: return new xval
+PetscErrorCode MMA::Update(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax) {
+    PetscErrorCode ierr = 0;
+
+    if (!NonLinConstraints) {
+        PetscErrorPrintf("MMA->Update called WITH constraints but object was "
+                         "allocated WITHOUT !\n");
+        return -1;
+    }
+
+    // Generate the subproblem
+    GenSub(xval, dfdx, gx, dgdx, xmin, xmax);
+
+    // Update xolds
+    VecCopy(xo1, xo2);
+    VecCopy(xval, xo1);
+
+    // Solve the dual with an interior point method
+    SolveDIP(xval);
+    return ierr;
+}
+
+// PRIVATE METHODS
+
+PetscErrorCode MMA::GenSub(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax) {
+    PetscErrorCode ierr = 0;
+
+    PetscScalar gamma, helpvar;
+
+    k++;
+    PetscInt nloc;
+    VecGetLocalSize(xval, &nloc);
+    PetscScalar *xv, *Lv, *Uv, *x1v, *x2v, *xminv, *xmaxv;
+    PetscScalar *alf, *bet, *dfdxv, *p0v, *q0v, **dgdxv, **pijv, **qijv;
+    if (k < 3) {
+        VecAXPBYPCZ(L, (PetscScalar)1.0, -asyminit, (PetscScalar)0.0, xval, xmax);
+        VecAXPY(L, asyminit, xmin);
+        VecAXPBYPCZ(U, (PetscScalar)1.0, +asyminit, (PetscScalar)0.0, xval, xmax);
+        VecAXPY(U, -asyminit, xmin);
+    }
+    VecGetArray(xval, &xv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    VecGetArray(xo1, &x1v);
+    VecGetArray(xo2, &x2v);
+    VecGetArray(xmin, &xminv);
+    VecGetArray(xmax, &xmaxv);
+
+    VecGetArray(alpha, &alf);
+    VecGetArray(beta, &bet);
+    VecGetArray(dfdx, &dfdxv);
+    VecGetArray(p0, &p0v);
+    VecGetArray(q0, &q0v);
+
+    VecGetArrays(dgdx, m, &dgdxv);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    if (k > 2) {
+        for (PetscInt i = 0; i < nloc; i++) {
+            helpvar = (xv[i] - x1v[i]) * (x1v[i] - x2v[i]);
+            if (helpvar < 0.0) {
+                gamma = asymdec;
+            } else if (helpvar > 0.0) {
+                gamma = asyminc;
+            } else {
+                gamma = 1.0;
+            }
+            Lv[i] = xv[i] - gamma * (x1v[i] - Lv[i]);
+            Uv[i] = xv[i] + gamma * (Uv[i] - x1v[i]);
+            PetscScalar xmi, xma;
+            xmi = Max(1.0e-5, xmaxv[i] - xminv[i]);
+            if (RobustAsymptotesType == 0) {
+                Lv[i] = Max(Lv[i], xv[i] - 10.0 * xmi);
+                Lv[i] = Min(Lv[i], xv[i] - 0.01 * xmi);
+                Uv[i] = Max(Uv[i], xv[i] + 0.01 * xmi);
+                Uv[i] = Min(Uv[i], xv[i] + 10.0 * xmi);
+            } else if (RobustAsymptotesType == 1) {
+                Lv[i] = Max(Lv[i], xv[i] - 100.0 * xmi);
+                Lv[i] = Min(Lv[i], xv[i] - 1.0e-4 * xmi);
+                Uv[i] = Max(Uv[i], xv[i] + 1.0e-4 * xmi);
+                Uv[i] = Min(Uv[i], xv[i] + 100.0 * xmi);
+                xmi   = xminv[i] - 1.0e-5;
+                xma   = xmaxv[i] + 1.0e-5;
+                if (xv[i] < xmi) {
+                    Lv[i] = xv[i] - (xma - xv[i]) / 0.9;
+                    Uv[i] = xv[i] + (xma - xv[i]) / 0.9;
+                }
+                if (xv[i] > xma) {
+                    Lv[i] = xv[i] - (xv[i] - xmi) / 0.9;
+                    Uv[i] = xv[i] + (xv[i] - xmi) / 0.9;
+                }
+            }
+        }
+    }
+    PetscScalar dfdxp, dfdxm;
+    PetscScalar feps = 1.0e-6;
+    for (PetscInt i = 0; i < nloc; i++) {
+        alf[i] = Max(xminv[i], 0.9 * Lv[i] + 0.1 * xv[i]);
+        bet[i] = Min(xmaxv[i], 0.9 * Uv[i] + 0.1 * xv[i]);
+        dfdxp  = Max(0.0, dfdxv[i]);
+        dfdxm  = Max(0.0, -1.0 * dfdxv[i]);
+        p0v[i] = pow(Uv[i] - xv[i], 2.0) * (dfdxp + 0.001 * Abs(dfdxv[i]) + 0.5 * feps / (Uv[i] - Lv[i]));
+        q0v[i] = pow(xv[i] - Lv[i], 2.0) * (dfdxm + 0.001 * Abs(dfdxv[i]) + 0.5 * feps / (Uv[i] - Lv[i]));
+        for (PetscInt j = 0; j < m; j++) {
+            dfdxp = Max(0.0, dgdxv[j][i]);
+            dfdxm = Max(0.0, -1.0 * dgdxv[j][i]);
+            if (constraintModification) {
+                pijv[j][i] =
+                    pow(Uv[i] - xv[i], 2.0) * (dfdxp + 0.001 * Abs(dgdxv[j][i]) + 0.5 * feps / (Uv[i] - Lv[i]));
+                qijv[j][i] =
+                    pow(xv[i] - Lv[i], 2.0) * (dfdxm + 0.001 * Abs(dgdxv[j][i]) + 0.5 * feps / (Uv[i] - Lv[i]));
+            } else {
+                pijv[j][i] = pow(Uv[i] - xv[i], 2.0) * (dfdxp);
+                qijv[j][i] = pow(xv[i] - Lv[i], 2.0) * (dfdxm);
+            }
+        }
+    }
+    for (PetscInt j = 0; j < m; j++) {
+        b[j] = 0.0;
+        for (PetscInt i = 0; i < nloc; i++) {
+            b[j] += pijv[j][i] / (Uv[i] - xv[i]) + qijv[j][i] / (xv[i] - Lv[i]);
+        }
+    }
+    {
+        PetscScalar* tmp = new PetscScalar[m];
+        for (PetscInt i = 0; i < m; i++) {
+            tmp[i] = 0.0;
+        }
+        MPI_Allreduce(b, tmp, m, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+        memcpy(b, tmp, sizeof(PetscScalar) * m);
+        delete[] tmp;
+    }
+    for (PetscInt j = 0; j < m; j++) {
+        b[j] += -gx[j];
+    }
+    VecRestoreArray(xval, &xv);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    VecRestoreArray(xo1, &x1v);
+    VecRestoreArray(xo2, &x2v);
+    VecRestoreArray(xmin, &xminv);
+    VecRestoreArray(xmax, &xmaxv);
+
+    VecRestoreArray(alpha, &alf);
+    VecRestoreArray(beta, &bet);
+    VecRestoreArray(dfdx, &dfdxv);
+
+    VecRestoreArrays(dgdx, m, &dgdxv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    return ierr;
+}
+
+PetscErrorCode MMA::SolveDIP(Vec x) {
+    PetscErrorCode ierr = 0;
+
+    for (PetscInt j = 0; j < m; j++) {
+        lam[j] = c[j] / 2.0;
+        mu[j]  = 1.0;
+    }
+    PetscScalar tol  = 1.0e-9 * sqrt(m + n);
+    PetscScalar epsi = 1.0;
+    PetscScalar err  = 1.0;
+    PetscInt    loop;
+    while (epsi > tol) {
+
+        loop = 0;
+        while (err > 0.9 * epsi && loop < 100) {
+            loop++;
+            XYZofLAMBDA(x);
+            DualGrad(x);
+            for (PetscInt j = 0; j < m; j++) {
+                grad[j] = -1.0 * grad[j] - epsi / lam[j];
+            }
+            DualHess(x);
+            Factorize(Hess, m);
+            Solve(Hess, grad, m);
+            for (PetscInt j = 0; j < m; j++) {
+                s[j] = grad[j];
+            }
+            for (PetscInt i = 0; i < m; i++) {
+                s[m + i] = -mu[i] + epsi / lam[i] - s[i] * mu[i] / lam[i];
+            }
+            DualLineSearch();
+            XYZofLAMBDA(x);
+            err = DualResidual(x, epsi);
+        }
+        epsi = epsi * 0.1;
+    }
+    return ierr;
+}
+
+PetscErrorCode MMA::XYZofLAMBDA(Vec x) {
+    PetscErrorCode ierr = 0;
+
+    PetscInt nloc;
+    VecGetLocalSize(x, &nloc);
+    PetscScalar *xv, **pijv, **qijv, *p0v, *q0v, *alf, *bet, *Lv, *Uv;
+    VecGetArray(x, &xv);
+    VecGetArray(p0, &p0v);
+    VecGetArray(q0, &q0v);
+    VecGetArray(alpha, &alf);
+    VecGetArray(beta, &bet);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    PetscScalar lamai = 0.0;
+    for (PetscInt i = 0; i < m; i++) {
+        if (lam[i] < 0.0) {
+            lam[i] = 0;
+        }
+        y[i] = Max(0.0, lam[i] - c[i]);
+        lamai += lam[i] * a[i];
+    }
+    z = Max(0.0, 10.0 * (lamai - 1.0)); // SINCE a0 = 1.0
+    PetscScalar pjlam, qjlam;
+    for (PetscInt i = 0; i < nloc; i++) {
+        pjlam = p0v[i];
+        qjlam = q0v[i];
+        for (PetscInt j = 0; j < m; j++) {
+            pjlam += pijv[j][i] * lam[j];
+            qjlam += qijv[j][i] * lam[j];
+        }
+        xv[i] = (sqrt(pjlam) * Lv[i] + sqrt(qjlam) * Uv[i]) / (sqrt(pjlam) + sqrt(qjlam));
+        if (xv[i] < alf[i]) {
+            xv[i] = alf[i];
+        }
+        if (xv[i] > bet[i]) {
+            xv[i] = bet[i];
+        }
+    }
+    VecRestoreArray(x, &xv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    VecRestoreArray(p0, &p0v);
+    VecRestoreArray(q0, &q0v);
+    VecRestoreArray(alpha, &alf);
+    VecRestoreArray(beta, &bet);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    return ierr;
+}
+
+PetscErrorCode MMA::DualGrad(Vec x) {
+    PetscErrorCode ierr = 0;
+
+    PetscInt nloc;
+    VecGetLocalSize(x, &nloc);
+    PetscScalar *xv, *Lv, *Uv, **pijv, **qijv;
+    VecGetArray(x, &xv);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    for (PetscInt j = 0; j < m; j++) {
+        grad[j] = 0.0;
+        for (PetscInt i = 0; i < nloc; i++) {
+            grad[j] += pijv[j][i] / (Uv[i] - xv[i]) + qijv[j][i] / (xv[i] - Lv[i]);
+        }
+    }
+    {
+        PetscScalar* tmp = new PetscScalar[m];
+        for (PetscInt i = 0; i < m; i++) {
+            tmp[i] = 0.0;
+        }
+        MPI_Allreduce(grad, tmp, m, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+        memcpy(grad, tmp, sizeof(PetscScalar) * m);
+        delete[] tmp;
+    }
+    for (PetscInt j = 0; j < m; j++) {
+        grad[j] += -b[j] - a[j] * z - y[j];
+    }
+    VecRestoreArray(x, &xv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    return ierr;
+}
+
+PetscErrorCode MMA::DualHess(Vec x) {
+    PetscErrorCode ierr = 0;
+
+    PetscInt nloc;
+    VecGetLocalSize(x, &nloc);
+    PetscScalar *xv, *Lv, *Uv, **pijv, **qijv, *alf, *bet, *p0v, *q0v;
+    VecGetArray(x, &xv);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    VecGetArray(alpha, &alf);
+    VecGetArray(beta, &bet);
+    VecGetArray(p0, &p0v);
+    VecGetArray(q0, &q0v);
+    PetscScalar* df2 = new PetscScalar[nloc];
+    PetscScalar* PQ  = new PetscScalar[nloc * m];
+    PetscScalar  pjlam, qjlam;
+    for (PetscInt i = 0; i < nloc; i++) {
+        pjlam = p0v[i];
+        qjlam = q0v[i];
+        for (PetscInt j = 0; j < m; j++) {
+            pjlam += pijv[j][i] * lam[j];
+            qjlam += qijv[j][i] * lam[j];
+            PQ[i * m + j] = pijv[j][i] / pow(Uv[i] - xv[i], 2.0) - qijv[j][i] / pow(xv[i] - Lv[i], 2.0);
+        }
+        df2[i]         = -1.0 / (2.0 * pjlam / pow(Uv[i] - xv[i], 3.0) + 2.0 * qjlam / pow(xv[i] - Lv[i], 3.0));
+        PetscScalar xp = (sqrt(pjlam) * Lv[i] + sqrt(qjlam) * Uv[i]) / (sqrt(pjlam) + sqrt(qjlam));
+        if (xp < alf[i]) {
+            df2[i] = 0.0;
+        }
+        if (xp > bet[i]) {
+            df2[i] = 0.0;
+        }
+    }
+    PetscScalar* tmp = new PetscScalar[n * m];
+    for (PetscInt j = 0; j < m; j++) {
+        for (PetscInt i = 0; i < nloc; i++) {
+            tmp[j * nloc + i] = 0.0;
+            tmp[j * nloc + i] += PQ[i * m + j] * df2[i];
+        }
+    }
+    for (PetscInt i = 0; i < m; i++) {
+        for (PetscInt j = 0; j < m; j++) {
+            Hess[i * m + j] = 0.0;
+            for (PetscInt k = 0; k < nloc; k++) {
+                Hess[i * m + j] += tmp[i * nloc + k] * PQ[k * m + j];
+            }
+        }
+    }
+    {
+        PetscScalar* tmpp = new PetscScalar[m * m];
+        for (PetscInt i = 0; i < m * m; i++) {
+            tmpp[i] = Hess[i];
+        }
+        MPI_Allreduce(Hess, tmpp, m * m, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+        memcpy(Hess, tmpp, sizeof(PetscScalar) * m * m);
+        delete[] tmpp;
+    }
+    PetscScalar lamai = 0.0;
+    for (PetscInt j = 0; j < m; j++) {
+        if (lam[j] < 0.0) {
+            lam[j] = 0.0;
+        }
+        lamai += lam[j] * a[j];
+        if (lam[j] > c[j]) {
+            Hess[j * m + j] += -1.0;
+        }
+        Hess[j * m + j] += -mu[j] / lam[j];
+    }
+    if (lamai > 0.0) {
+        for (PetscInt j = 0; j < m; j++) {
+            for (PetscInt k = 0; k < m; k++) {
+                Hess[j * m + k] += -10.0 * a[j] * a[k];
+            }
+        }
+    }
+    PetscScalar HessTrace = 0.0;
+    for (PetscInt i = 0; i < m; i++) {
+        HessTrace += Hess[i * m + i];
+    }
+    PetscScalar HessCorr = 1e-4 * HessTrace / m;
+    if (-1.0 * HessCorr < 1.0e-7) {
+        HessCorr = -1.0e-7;
+    }
+    for (PetscInt i = 0; i < m; i++) {
+        Hess[i * m + i] += HessCorr;
+    }
+    VecRestoreArray(x, &xv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    VecRestoreArray(q0, &q0v);
+    VecRestoreArray(p0, &p0v);
+    VecRestoreArray(alpha, &alf);
+    VecRestoreArray(beta, &bet);
+    delete[] df2;
+    delete[] PQ;
+    delete[] tmp;
+    return ierr;
+}
+
+PetscErrorCode MMA::DualLineSearch() {
+    PetscErrorCode ierr = 0;
+
+    PetscScalar theta = 1.005;
+    for (PetscInt i = 0; i < m; i++) {
+        if (theta < -1.01 * s[i] / lam[i]) {
+            theta = -1.01 * s[i] / lam[i];
+        }
+        if (theta < -1.01 * s[i + m] / mu[i]) {
+            theta = -1.01 * s[i + m] / mu[i];
+        }
+    }
+    theta = 1.0 / theta;
+    for (PetscInt i = 0; i < m; i++) {
+        lam[i] = lam[i] + theta * s[i];
+        mu[i]  = mu[i] + theta * s[i + m];
+    }
+    return ierr;
+}
+
+PetscScalar MMA::DualResidual(Vec x, PetscScalar epsi) {
+
+    PetscInt nloc;
+    VecGetLocalSize(x, &nloc);
+    PetscScalar* res = new PetscScalar[2 * m];
+    PetscScalar *xv, *Lv, *Uv, **pijv, **qijv;
+    VecGetArray(x, &xv);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    for (PetscInt j = 0; j < m; j++) {
+        res[j]     = 0.0;
+        res[j + m] = 0.0;
+        for (PetscInt i = 0; i < nloc; i++) {
+            res[j] += pijv[j][i] / (Uv[i] - xv[i]) + qijv[j][i] / (xv[i] - Lv[i]);
+        }
+    }
+    {
+        PetscScalar* tmp = new PetscScalar[2 * m];
+        for (PetscInt i = 0; i < 2 * m; i++) {
+            tmp[i] = 0.0;
+        }
+        MPI_Allreduce(res, tmp, 2 * m, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+        memcpy(res, tmp, sizeof(PetscScalar) * 2 * m);
+        delete[] tmp;
+    }
+    for (PetscInt j = 0; j < m; j++) {
+        res[j] += -b[j] - a[j] * z - y[j] + mu[j];
+        res[j + m] += mu[j] * lam[j] - epsi;
+    }
+    PetscScalar nrI = 0.0;
+    for (PetscInt i = 0; i < 2 * m; i++) {
+        if (nrI < Abs(res[i])) {
+            nrI = Abs(res[i]);
+        }
+    }
+    delete[] res;
+    VecRestoreArray(x, &xv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    return nrI;
+}
+
+PetscErrorCode MMA::Factorize(PetscScalar* K, PetscInt nn) {
+    PetscErrorCode ierr = 0;
+
+    for (PetscInt ss = 0; ss < nn - 1; ss++) {
+        for (PetscInt i = ss + 1; i < nn; i++) {
+            K[i * nn + ss] = K[i * nn + ss] / K[ss * nn + ss];
+            for (PetscInt j = ss + 1; j < nn; j++) {
+                K[i * nn + j] = K[i * nn + j] - K[i * nn + ss] * K[ss * nn + j];
+            }
+        }
+    }
+    return ierr;
+}
+
+PetscErrorCode MMA::Solve(PetscScalar* K, PetscScalar* x, PetscInt nn) {
+    PetscErrorCode ierr = 0;
+
+    for (PetscInt i = 1; i < nn; i++) {
+        PetscScalar a = 0.0;
+        for (PetscInt j = 0; j < i; j++) {
+            a = a - K[i * nn + j] * x[j];
+        }
+        x[i] = x[i] + a;
+    }
+    x[nn - 1] = x[nn - 1] / K[(nn - 1) * nn + (nn - 1)];
+    for (PetscInt i = nn - 2; i >= 0; i--) {
+        PetscScalar a = x[i];
+        for (PetscInt j = i + 1; j < nn; j++) {
+            a = a - K[i * nn + j] * x[j];
+        }
+        x[i] = a / K[i * nn + i];
+    }
+    return ierr;
+}
+
+PetscScalar MMA::Min(PetscScalar d1, PetscScalar d2) { return d1 < d2 ? d1 : d2; }
+
+PetscScalar MMA::Max(PetscScalar d1, PetscScalar d2) { return d1 > d2 ? d1 : d2; }
+
+PetscInt MMA::Min(PetscInt d1, PetscInt d2) { return d1 < d2 ? d1 : d2; }
+
+PetscInt MMA::Max(PetscInt d1, PetscInt d2) { return d1 > d2 ? d1 : d2; }
+
+PetscScalar MMA::Abs(PetscScalar d1) { return d1 > 0 ? d1 : -1.0 * d1; }
+
+#endif

miniapps/meshing/MMA.hpp

@@ -0,0 +1,242 @@
+#ifndef __MMA__HPP
+#define __MMA__HPP
+
+#include "mfem.hpp"
+
+#ifdef MFEM_USE_PETSC
+#include <petsc.h>
+
+/* -----------------------------------------------------------------------------
+Authors: Niels Aage
+ Copyright (C) 2013-2020,
+This MMA implementation is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+This Module is distributed in the hope that it will be useful,implementation
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Lesser General Public License for more details.
+You should have received a copy of the GNU Lesser General Public
+License along with this Module; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+-------------------------------------------------------------------------- */
+class MMA {
+  public:
+    // Construct using defaults subproblem penalization
+    MMA(PetscInt n, PetscInt m, Vec x);
+    // User defined subproblem penalization
+    MMA(PetscInt n, PetscInt m, Vec x, PetscScalar* a, PetscScalar* c, PetscScalar* d);
+    // Initialize with restart from itr
+    MMA(PetscInt n, PetscInt m, PetscInt itr, Vec xo1, Vec xo2, Vec U, Vec L);
+    // Initialize with restart and specify subproblem parameters
+    MMA(PetscInt n, PetscInt m, PetscInt itr, Vec xo1, Vec xo2, Vec U, Vec L, PetscScalar* a, PetscScalar* c,
+        PetscScalar* d);
+    // Destructor
+    ~MMA();
+
+    // Set and solve a subproblem: return new xval
+    PetscErrorCode Update(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax);
+
+    // Return necessary data for possible restart
+    PetscErrorCode Restart(Vec xo1, Vec xo2, Vec U, Vec L);
+
+    // Set the aggresivity of the moving asymptotes
+    PetscErrorCode SetAsymptotes(PetscScalar init, PetscScalar decrease, PetscScalar increase);
+
+    // do/don't add convexity approx to constraints: default=false
+    PetscErrorCode ConstraintModification(PetscBool conMod) {
+        constraintModification = conMod;
+        return 0;
+    };
+
+    // val=0: default, val=1: increase robustness, i.e
+    // control the spacing between L < alp < x < beta < U,
+    PetscErrorCode SetRobustAsymptotesType(PetscInt val);
+
+    // Sets outer movelimits on all primal design variables
+    // This is often requires to prevent the solver from oscilating
+    PetscErrorCode SetOuterMovelimit(PetscScalar Xmin, PetscScalar Xmax, PetscScalar movelim, Vec x, Vec xmin,
+                                     Vec xmax);
+
+    // Return KKT residual norms (norm2 and normInf)
+    PetscErrorCode KKTresidual(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax, PetscScalar* norm2,
+                               PetscScalar* normInf);
+
+    // Inf norm on diff between two vectors: SHOULD NOT BE HERE - USE BASIC
+    // PETSc!!!!!
+    PetscScalar DesignChange(Vec x, Vec xold);
+
+  private:
+    // Set up the MMA subproblem based on old x's and xval
+    PetscErrorCode GenSub(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax);
+
+    // Interior point solver for the subproblem
+    PetscErrorCode SolveDIP(Vec xval);
+
+    // Compute primal vars based on dual solution
+    PetscErrorCode XYZofLAMBDA(Vec x);
+
+    // Dual gradient
+    PetscErrorCode DualGrad(Vec x);
+
+    // Dual Hessian
+    PetscErrorCode DualHess(Vec x);
+
+    // Dual line search
+    PetscErrorCode DualLineSearch();
+
+    // Dual residual
+    PetscScalar DualResidual(Vec x, PetscScalar epsi);
+
+    // Problem size and iteration counter
+    PetscInt n, m, k;
+
+    // "speed-control" for the asymptotes
+    PetscScalar asyminit, asymdec, asyminc;
+
+    // do/don't add convexity constraint approximation in subproblem
+    PetscBool constraintModification; // default = FALSE
+
+    // Bool specifying if non lin constraints are included or not
+    PetscBool NonLinConstraints;
+
+    // 0: (default) span between alp L x U beta,
+    // 1: increase the span for further robustness
+    PetscInt RobustAsymptotesType;
+
+    // Local vectors: penalty numbers for subproblem
+    PetscScalar *a, *c, *d;
+
+    // Local vectors: elastic variables
+    PetscScalar* y;
+    PetscScalar  z;
+
+    // Local vectors: Lagrange multipliers:
+    PetscScalar *lam, *mu, *s;
+
+    // Global: Asymptotes, bounds, objective approx., constraint approx.
+    Vec L, U, alpha, beta, p0, q0, *pij, *qij;
+
+    // Local: subproblem constant terms, dual gradient, dual hessian
+    PetscScalar *b, *grad, *Hess;
+
+    // Global: Old design variables
+    Vec xo1, xo2;
+
+    // Math helpers
+    PetscErrorCode Factorize(PetscScalar* K, PetscInt nn);
+    PetscErrorCode Solve(PetscScalar* K, PetscScalar* x, PetscInt nn);
+    PetscScalar    Min(PetscScalar d1, PetscScalar d2);
+    PetscScalar    Max(PetscScalar d1, PetscScalar d2);
+    PetscInt       Min(PetscInt d1, PetscInt d2);
+    PetscInt       Max(PetscInt d1, PetscInt d2);
+    PetscScalar    Abs(PetscScalar d1);
+};
+
+
+namespace mfem {
+
+class NativeMMA
+{
+public:
+    // User defined subproblem penalization
+    NativeMMA(MPI_Comm comm_, int m, mfem::Vector& x, double* a, double* c, double* d){
+
+        comm=comm_;
+        num_con=m;
+        Vec pv;
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&pv);
+
+        PetscScalar* ap=new PetscScalar[m];
+        PetscScalar* cp=new PetscScalar[m];
+        PetscScalar* dp=new PetscScalar[m];
+        for(int i=0;i<m;i++)
+        {
+            ap[i]=a[i];
+            cp[i]=c[i];
+            dp[i]=c[i];
+        }
+
+        PetscInt nn;
+        VecGetSize(pv,&nn);
+
+        mma=new ::MMA(nn, m, pv, ap, cp, dp);
+
+        delete [] ap;
+        delete [] cp;
+        delete [] dp;
+
+        VecDestroy(&pv);
+
+        //allocate the native PETSc objects necessary for the subproblems
+
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&xval);
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&dfdx);
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&xmin);
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&xmax);
+
+        VecDuplicateVecs(xval,num_con, &dgdx);
+    }
+
+    ~NativeMMA()
+    {
+        delete mma;
+
+        VecDestroy(&xval);
+        VecDestroy(&dfdx);
+        VecDestroy(&xmin);
+        VecDestroy(&xmax);
+        VecDestroyVecs(num_con,&dgdx);
+    }
+
+    // Set and solve a subproblem: return new xval
+    void Update(Vector& xval_, Vector& dfdx_, double* gx_, Vector* dgdx_, Vector& xmin_, Vector& xmax_){
+        //copy data
+        double* data;
+        VecGetArray(xval,&data);
+        for(int i=0;i<xval_.Size();i++){data[i]=xval_[i];}
+        VecRestoreArray(xval,&data);
+        //dfdx
+        VecGetArray(dfdx,&data);
+        for(int i=0;i<xval_.Size();i++){data[i]=dfdx_[i];}
+        VecRestoreArray(dfdx,&data);
+        //dgdx
+        for(int j=0;j<num_con;j++){
+            VecGetArray(dgdx[j],&data);
+            for(int i=0;i<xval_.Size();i++){data[i]=(dgdx_[j])[i];}
+            VecRestoreArray(dgdx[j],&data);
+        }
+        //xmin
+        VecGetArray(xmin,&data);
+        for(int i=0;i<xval_.Size();i++){data[i]=xmin_[i];}
+        VecRestoreArray(xmin,&data);
+        //xmax
+        VecGetArray(xmax,&data);
+        for(int i=0;i<xval_.Size();i++){data[i]=xmax_[i];}
+        VecRestoreArray(xmax,&data);
+
+        mma->Update(xval,dfdx,gx_,dgdx,xmin,xmax);
+
+        VecGetArray(xval,&data);
+        for(int i=0;i<xval_.Size();i++){xval_[i]=data[i];}
+        VecRestoreArray(xval,&data);
+    }
+private:
+    ::MMA* mma;
+    int num_con;
+    MPI_Comm comm;
+
+    Vec xval;
+    Vec dfdx;
+    Vec* dgdx;
+    Vec xmin;
+    Vec xmax;
+
+
+};
+
+}
+
+#endif
+#endif

miniapps/meshing/makefile

@@ -27,7 +27,8 @@ MFEM_LIB_FILE = mfem_is_not_built
 SEQ_MINIAPPS = mobius-strip klein-bottle toroid trimmer twist mesh-explorer\
 	shaper extruder mesh-optimizer minimal-surface polar-nc reflector\
 	mesh-quality
-PAR_MINIAPPS = pmesh-optimizer pminimal-surface pmesh-fitting fit-node-position
+PAR_MINIAPPS = pmesh-optimizer pminimal-surface pmesh-fitting fit-node-position\
+               pmesh-optimizer_NLP
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -58,11 +59,24 @@ all: $(MINIAPPS)
 lib-common:
 	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
 
+# Rule for building pmesh-optimizer_NLP
+pmesh-optimizer_NLP: pmesh-optimizer_NLP.o mesh-optimizer_using_NLP.o \
+        $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_FLAGS) -o pmesh-optimizer_NLP pmesh-optimizer_NLP.o\
+        mesh-optimizer_using_NLP.o $(COMMON_LIB) $(MFEM_LIBS)
+
+mesh-optimizer_using_NLP.o: mesh-optimizer_using_NLP.cpp \
+    mesh-optimizer_using_NLP.hpp
+	$(MFEM_CXX) $(MFEM_FLAGS) -c mesh-optimizer_using_NLP.cpp
+
+pmesh-optimizer_NLP.o: pmesh-optimizer_NLP.cpp mesh-optimizer_using_NLP.cpp
+	$(MFEM_CXX) $(MFEM_FLAGS) -c pmesh-optimizer_NLP.cpp
+
 # Rules to copy the *.mesh files - needed for running the sample runs when
 # building out-of-source:
 ifneq ($(SRC),)
 MESH_FILES = amr-quad-q2.mesh blade.mesh cube.mesh cube-tet.mesh icf.mesh\
-jagged.mesh square01.mesh square01-tri.mesh stretched2D.mesh 
+jagged.mesh square01.mesh square01-tri.mesh stretched2D.mesh
 $(MESH_FILES): %: $(SRC)%
 	ln -sf $(<) .
 mesh-optimizer pmesh-optimizer pmesh-fitting fit-node-position: | $(MESH_FILES)

miniapps/meshing/mesh-optimizer.cpp

@@ -474,7 +474,6 @@ int main(int argc, char *argv[])
       // A-metrics
       case 11: metric = new TMOP_AMetric_011; break;
       case 36: metric = new TMOP_AMetric_036; break;
-      case 107: metric = new TMOP_AMetric_107a; break;
       case 126: metric = new TMOP_AMetric_126(0.9); break;
       default:
          cout << "Unknown metric_id: " << metric_id << endl;

miniapps/meshing/mesh-optimizer_using_NLP.cpp

@@ -0,0 +1,932 @@
+#include "mesh-optimizer_using_NLP.hpp"
+#include "mfem.hpp"
+
+#ifdef MFEM_USE_PETSC
+#include "petsc.h"
+#endif
+
+namespace mfem {
+
+// void IdentityMatrix(int dim, mfem::DenseMatrix &I)
+// {
+//   I.SetSize(dim, dim);
+//   I = 0.0;
+//   for (int i = 0; i < dim; i++) {
+//     I(i, i) = 1.0;
+//   }
+// }
+
+// void Vectorize(const mfem::DenseMatrix &A, mfem::Vector &a)
+// {
+//   int m = A.NumRows();
+//   int n = A.NumCols();
+//   a.SetSize(m * n);
+//   int k = 0;
+//   for (int j = 0; j < n; j++) {
+//     for (int i = 0; i < m; i++) {
+//       a(k++) = A(i, j);
+//     }
+//   }
+// }
+
+// double MatrixInnerProduct(const mfem::DenseMatrix &A, const mfem::DenseMatrix &B)
+// {
+//   double inner_product = 0.0;
+//   for (int i = 0; i < A.NumRows(); i++) {
+//     for (int j = 0; j < A.NumCols(); j++) {
+//       inner_product += A(i, j) * B(i, j);
+//     }
+//   }
+//   return inner_product;
+// }
+
+// void ConjugationProduct(const mfem::DenseMatrix &A, const mfem::DenseMatrix &B, const mfem::DenseMatrix &C, mfem::DenseMatrix &D)
+// {
+//   mfem::DenseMatrix CBt(C.NumRows(), B.NumRows());
+//   mfem::MultABt(C, B, CBt);
+//   mfem::Mult(A, CBt, D);
+// }
+
+//  LFNodeCoordinateSensitivityIntegrator::LFNodeCoordinateSensitivityIntegrator( int Index1, int Index2,
+//     int IntegrationOrder)
+//   : Index1_(Index1), Index2_(Index2), IntegrationOrder_(IntegrationOrder)
+// {}
+
+// void LFNodeCoordinateSensitivityIntegrator::AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T,
+//     mfem::Vector &elvect)
+// {
+//   // grab sizes
+//   int dof = el.GetDof();
+//   int dim = el.GetDim();
+
+//   // initialize storage
+//   mfem::DenseMatrix dN(dof, dim);
+//   mfem::DenseMatrix NxPhix(dof, dim);
+//   mfem::DenseMatrix B(dof, dim);
+//   mfem::DenseMatrix IxB(dof, dim);
+//   mfem::Vector IxBTvec(dof * dim);
+//   mfem::Vector IxN_vec(dof * dim);
+//   mfem::Vector N(dof);
+//   mfem::DenseMatrix graduDerivxB(dof, dim);
+//   mfem::Vector graduDerivxBvec(dof * dim);
+
+//   // output vector
+//   elvect.SetSize(dim * dof);
+//   elvect = 0.0;
+
+//   int integrationOrder = 2 * el.GetOrder() + 2;
+//   if (IntegrationOrder_ != INT_MAX) {
+//     integrationOrder = IntegrationOrder_;
+//   }
+
+//   // set integration rule
+//   const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), integrationOrder);
+
+//   // identity tensor
+//   mfem::DenseMatrix I;
+//   IdentityMatrix(dim, I);
+
+//   // loop over integration points
+//   for (int i = 0; i < ir->GetNPoints(); i++) {
+//     // set current integration point
+//     const ::mfem::IntegrationPoint &ip = ir->IntPoint(i);
+//     T.SetIntPoint(&ip);
+
+//     // evaluate gaussian integration weight
+//     double w = ip.weight * T.Weight();
+
+//     // evaluate shape function derivative
+//     el.CalcDShape(ip, dN);
+//     el.CalcShape(ip, N);
+//     mfem::DenseMatrix matN(N.GetData(), dof, 1);
+
+//     // get inverse jacobian
+//     mfem::DenseMatrix Jinv = T.InverseJacobian();
+//     mfem::Mult(dN, Jinv, B);
+
+//     // term 1
+//     mfem::Mult(B, QoI_->gradTimesexplicitSolutionGradientDerivative(T, ip), graduDerivxB);
+//     Vectorize(graduDerivxB, graduDerivxBvec);
+//     elvect.Add( -1.0 * w , graduDerivxBvec);
+
+//     // term 2
+//     mfem::Mult(B, I, IxB);
+//     Vectorize(IxB, IxBTvec);
+//     elvect.Add( w * QoI_->Eval(T, ip), IxBTvec);
+
+//     // term 3
+//     Mult(matN, QoI_->explicitShapeDerivative(T, ip), NxPhix);
+//     Vectorize(NxPhix, IxN_vec);
+//     elvect.Add(w , IxN_vec);
+//   }
+// }
+
+// LFErrorIntegrator::LFErrorIntegrator( int IntegrationOrder)
+//   : IntegrationOrder_(IntegrationOrder)
+// {}
+
+// void LFErrorIntegrator::AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T,
+//     mfem::Vector &elvect)
+// {
+//   // grab sizes
+//   int dof = el.GetDof();
+//   //int dim = el.GetDim();
+
+//   // initialize storage
+//   mfem::Vector N(dof);
+
+//   // output vector
+//   elvect.SetSize( dof);
+//   elvect = 0.0;
+
+//   int integrationOrder = 2 * el.GetOrder() + 2;
+//   if (IntegrationOrder_ != INT_MAX) {
+//     integrationOrder = IntegrationOrder_;
+//   }
+
+//   // set integration rule
+//   const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), integrationOrder);
+
+//   // loop over integration points
+//   for (int i = 0; i < ir->GetNPoints(); i++) {
+//     // set current integration point
+//     const ::mfem::IntegrationPoint &ip = ir->IntPoint(i);
+//     T.SetIntPoint(&ip);
+
+//     // evaluate gaussian integration weight
+//     double w = ip.weight * T.Weight();
+
+//     el.CalcShape(ip, N);
+//     elvect.Add(w  * QoI_->Eval(T, ip), N);
+//   }
+// }
+
+// LFErrorDerivativeIntegrator::LFErrorDerivativeIntegrator( int IntegrationOrder)
+//   : IntegrationOrder_(IntegrationOrder)
+// {}
+
+// void LFErrorDerivativeIntegrator::AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T,
+//     mfem::Vector &elvect)
+// {
+//   // grab sizes
+//   int dof = el.GetDof();
+//   int dim = el.GetDim();
+
+//   // initialize storage
+//   mfem::Vector N(dof);
+//   mfem::DenseMatrix dN(dof, dim);
+//   mfem::DenseMatrix B(dof, dim);
+//   mfem::DenseMatrix BT(dim, dof);
+//   mfem::Vector DQdgradxdNdxvec(dof);
+//   mfem::DenseMatrix DQdgradxdNdx(1, dof);
+
+//   // output vector
+//   elvect.SetSize( dof);
+//   elvect = 0.0;
+
+//   int integrationOrder = 2 * el.GetOrder() + 2;
+//   if (IntegrationOrder_ != INT_MAX) {
+//     integrationOrder = IntegrationOrder_;
+//   }
+
+//   // set integration rule
+//   const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), integrationOrder);
+
+//   // loop over integration points
+//   for (int i = 0; i < ir->GetNPoints(); i++) {
+//     // set current integration point
+//     const ::mfem::IntegrationPoint &ip = ir->IntPoint(i);
+//     T.SetIntPoint(&ip);
+
+//     // evaluate gaussian integration weight
+//     double w = ip.weight * T.Weight();
+
+//     el.CalcDShape(ip, dN);
+//     el.CalcShape(ip, N);
+//     // get inverse jacobian
+//     mfem::DenseMatrix Jinv = T.InverseJacobian();
+//     mfem::Mult(dN, Jinv, B);
+
+//     // term 1
+//     const mfem::DenseMatrix & SolGradDeriv = QoI_->explicitSolutionGradientDerivative(T, ip);
+
+//     BT.Transpose(B);
+
+//     mfem::Mult(SolGradDeriv, BT, DQdgradxdNdx);
+//     Vectorize(DQdgradxdNdx, DQdgradxdNdxvec);
+//     elvect.Add(w , DQdgradxdNdxvec);
+
+//     //term 2
+//     const mfem::DenseMatrix & derivVal = QoI_->explicitSolutionDerivative(T, ip);
+//     double val = derivVal.Elem(0,0);
+
+//     elvect.Add(w  * val, N);
+//   }
+// }
+
+// LFErrorDerivativeIntegrator_2::LFErrorDerivativeIntegrator_2( ParFiniteElementSpace * fespace, Array<int> count,int IntegrationOrder)
+//   : fespace_(fespace), count_(count), IntegrationOrder_(IntegrationOrder)
+// {}
+
+// void LFErrorDerivativeIntegrator_2::AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T,
+//     mfem::Vector &elvect)
+// {
+//   // grab sizes
+//   int dof = el.GetDof();
+//   int dim = el.GetDim();
+
+//   int integrationOrder = 2 * el.GetOrder() + 2;
+//   if (IntegrationOrder_ != INT_MAX) {
+//     integrationOrder = IntegrationOrder_;
+//   }
+
+//   // set integration rule
+//   const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), integrationOrder);
+
+//   // initialize storage
+//   mfem::Vector N(dof);
+//   mfem::Vector shape(dof);
+//   mfem::DenseMatrix dN(dof, dim);
+//   mfem::DenseMatrix B(dof, dim);
+//   mfem::DenseMatrix BT(dim, dof);
+//   mfem::Vector DQdgradxdNdxvec(dof);
+//   mfem::DenseMatrix DQdgradxdNdx(1, dof);
+
+//   mfem::DenseMatrix SumdNdxatNodes(dof, dof*dim); SumdNdxatNodes = 0.0;
+
+//   //-----------------------------------------------------------------------------------
+
+//   const mfem::IntegrationRule *ir_p = &el.GetNodes();
+//   int fnd = ir_p->GetNPoints();
+//   //flux.SetSize( fnd * spaceDim );
+
+//   DenseMatrix dshape(dof,dim), invdfdx(dim, dim);
+
+//   int EleNo = T.ElementNo;
+
+//   Array<int> fdofs;
+//   fespace_->GetElementVDofs(EleNo, fdofs);
+
+//   for (int i = 0; i < fnd; i++)
+//   {
+//       const IntegrationPoint &ip = ir_p->IntPoint(i);
+//       el.CalcDShape(ip, dshape);
+//       el.CalcShape(ip, shape);
+//       // dshape.MultTranspose(u, vec);
+
+//       T.SetIntPoint (&ip);
+//       CalcInverse(T.Jacobian(), invdfdx);
+//       mfem::Mult(dshape, invdfdx, B);
+//       //invdfdx.MultTranspose(vec, vecdxt);
+
+//       std::cout<<"dshape H | W: "<<dshape.Height()<<" | "<<dshape.Width()<<std::endl;
+//       std::cout<<"B H | W: "<<B.Height()<<" | "<<B.Width()<<std::endl;
+//       std::cout<<"shape S: "<<shape.Size()<< std::endl;
+//       std::cout<<"count: "<<count_[fdofs[i]]<< std::endl;
+
+//     //  B.Print();
+//     //  SumdNdxatNodes.Print();
+
+//       for( int jj = 0; jj < dof; jj++)
+//       {
+
+//         // deveide here by node weights  // FIXME
+//         SumdNdxatNodes(jj,i) = B(jj,0) / count_[fdofs[i]];
+//         SumdNdxatNodes(jj,dof+i) = B(jj,1) / count_[fdofs[i]];
+
+//         // SumdNdxatNodes(jj,i) = B(jj,0) / 4.0;
+//         // SumdNdxatNodes(jj,dof+i) = B(jj,1) / 4.0;
+//       }
+
+//       // shape.Print();
+
+//   }
+
+//   // SumdNdxatNodes.Print();
+
+//    //mfem_error("stop in LFErrorDerivativeIntegrator_2");
+
+
+
+
+
+//   //-----------------------------------------------------------------------------------
+
+
+
+//   // output vector
+//   elvect.SetSize( dof);
+//   elvect = 0.0;
+
+//   // loop over integration points
+//   for (int i = 0; i < ir->GetNPoints(); i++) {
+//     // set current integration point
+//     const ::mfem::IntegrationPoint &ip = ir->IntPoint(i);
+//     T.SetIntPoint(&ip);
+
+//     // evaluate gaussian integration weight
+//     double w = ip.weight * T.Weight();
+
+//     el.CalcDShape(ip, dN);
+//     el.CalcShape(ip, N);
+//     // get inverse jacobian
+//     mfem::DenseMatrix Jinv = T.InverseJacobian();
+//     mfem::Mult(dN, Jinv, B);
+
+//     // term 1
+//     const mfem::DenseMatrix & SolGradDeriv = QoI_->explicitSolutionGradientDerivative(T, ip);
+
+//     BT.Transpose(B);
+
+//     mfem::Mult(SolGradDeriv, BT, DQdgradxdNdx);
+//     Vectorize(DQdgradxdNdx, DQdgradxdNdxvec);
+//     elvect.Add(w , DQdgradxdNdxvec);
+
+//     //term 2
+//     mfem::DenseMatrix Nhelper(dof*dim, dim); Nhelper = 0.0;
+//     mfem::DenseMatrix NtimesSumdNdx(dof, dim); NtimesSumdNdx = 0.0;
+//     mfem::DenseMatrix NtimesSumdNdxT(dim, dof); NtimesSumdNdxT = 0.0;
+
+//       for( int jj = 0; jj < dof; jj++)
+//       {
+//         Nhelper(jj,0) =     N(jj);
+//         Nhelper(dof+jj,1) = N(jj);
+//       }
+
+//     mfem::Mult(SumdNdxatNodes, Nhelper, NtimesSumdNdx);
+
+
+//     NtimesSumdNdxT.Transpose(NtimesSumdNdx);
+//     mfem::Mult(SolGradDeriv, NtimesSumdNdxT, DQdgradxdNdx);
+//     Vectorize(DQdgradxdNdx, DQdgradxdNdxvec);
+//     elvect.Add(-1.0*w , DQdgradxdNdxvec);
+
+
+
+//     //term 2
+//     // const mfem::DenseMatrix & derivVal = QoI_->explicitSolutionDerivative(T, ip);
+//     // double val = derivVal.Elem(0,0);
+
+//     // elvect.Add(w  * val, N);
+//   }
+// }
+
+// ThermalConductivityShapeSensitivityIntegrator::ThermalConductivityShapeSensitivityIntegrator(
+//   mfem::Coefficient &conductivity, const mfem::ParGridFunction &t_primal, const mfem::ParGridFunction &t_adjoint)
+//   : k_(&conductivity), t_primal_(&t_primal), t_adjoint_(&t_adjoint)
+// {}
+
+// void ThermalConductivityShapeSensitivityIntegrator::AssembleRHSElementVect(const mfem::FiniteElement &el,
+//     mfem::ElementTransformation &T, mfem::Vector &elvect)
+// {
+//   // grab sizes
+//   int dof = el.GetDof();
+//   int dim = el.GetDim();
+
+//   // initialize storage
+//   mfem::DenseMatrix dN(dof, dim);
+//   mfem::DenseMatrix B(dof, dim);
+
+//   mfem::DenseMatrix BBT(dof, dof);
+//   mfem::DenseMatrix dBBT(dof, dof);
+//   mfem::DenseMatrix BdBT(dof, dof);
+
+//   mfem::DenseMatrix dX_dXk(dim, dof);
+//   mfem::DenseMatrix dJ_dXk(dim, dim);
+//   mfem::DenseMatrix dJinv_dXk(dim, dim);
+//   mfem::DenseMatrix dB_dXk(dof, dim);
+
+//   mfem::DenseMatrix dK_dXk(dof, dof);
+
+//   mfem::Vector te_primal(dof);
+//   mfem::Vector te_adjoint(dof);
+
+//   mfem::Array<int> vdofs;
+//   t_primal_->ParFESpace()->GetElementVDofs(T.ElementNo, vdofs);
+//   t_primal_->GetSubVector(vdofs, te_primal);
+
+//   t_adjoint_->ParFESpace()->GetElementVDofs(T.ElementNo, vdofs);
+//   t_adjoint_->GetSubVector(vdofs, te_adjoint);
+
+//   // output matrix
+//   elvect.SetSize(dim * dof);
+//   elvect = 0.0;
+
+//   // set integration rule
+//   const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * T.OrderGrad(&el));
+
+//   // loop over nodal coordinates (X_k)
+//   for (int m = 0; m < dim; m++) {
+//     for (int n = 0; n < dof; n++) {
+//       dX_dXk = 0.0;
+//       dX_dXk(m, n) = 1.0;
+
+//       dK_dXk = 0.0;
+//       // loop over integration points
+//       for (int i = 0; i < ir->GetNPoints(); i++) {
+//         const ::mfem::IntegrationPoint &ip = ir->IntPoint(i);
+//         T.SetIntPoint(&ip);
+
+//         double w = ip.weight * T.Weight();
+//         el.CalcDShape(ip, dN);
+//         mfem::Mult(dN, T.InverseJacobian(), B);
+
+//         // compute derivative of Jacobian w.r.t. nodal coordinate
+//         mfem::Mult(dX_dXk, dN, dJ_dXk);
+
+//         // compute derivative of J^(-1)
+//         mfem::DenseMatrix JinvT = T.InverseJacobian();
+//         JinvT.Transpose();
+//         ConjugationProduct(T.InverseJacobian(), JinvT, dJ_dXk, dJinv_dXk);
+//         dJinv_dXk *= -1.0;
+
+//         // compute derivative of B w.r.t. nodal coordinate
+//         mfem::Mult(dN, dJinv_dXk, dB_dXk);
+
+//         // compute derivative of stiffness matrix w.r.t. X_k
+//         mfem::MultAAt(B, BBT);
+//         mfem::MultABt(dB_dXk, B, dBBT);
+//         mfem::MultABt(B, dB_dXk, BdBT);
+
+//         // compute derivative of integration weight w.r.t. X_k
+//         double dw_dXk = w * MatrixInnerProduct(JinvT, dJ_dXk);
+
+//         // put together all terms of product rule
+//         double k = k_->Eval(T, ip);
+//         dK_dXk.Add(w * k, dBBT);
+//         dK_dXk.Add(w * k, BdBT);
+//         dK_dXk.Add(dw_dXk * k, BBT);
+//       }
+//       elvect(n + m * dof) += dK_dXk.InnerProduct(te_primal, te_adjoint);
+//     }
+//   }
+// }
+
+// ThermalHeatSourceShapeSensitivityIntegrator::ThermalHeatSourceShapeSensitivityIntegrator(mfem::Coefficient &heatSource,
+//     const mfem::ParGridFunction &t_adjoint,
+//     int oa, int ob)
+//   : Q_(&heatSource), t_adjoint_(&t_adjoint), oa_(oa), ob_(ob)
+// {}
+
+// void ThermalHeatSourceShapeSensitivityIntegrator::AssembleRHSElementVect(const mfem::FiniteElement &el,
+//     mfem::ElementTransformation &T, mfem::Vector &elvect)
+// {
+//   // grab sizes
+//   int dof = el.GetDof();
+//   int dim = el.GetDim();
+
+//   // initialize storage
+//   mfem::DenseMatrix dX_dXk(dim, dof);
+//   mfem::DenseMatrix dJ_dXk(dim, dim);
+//   mfem::Vector N(dof);
+//   mfem::DenseMatrix dN(dof, dim);
+//   mfem::Vector dp_dXk(dof);
+
+//   mfem::Vector te_adjoint(dof);
+
+//   mfem::Array<int> vdofs;
+//   t_adjoint_->ParFESpace()->GetElementVDofs(T.ElementNo, vdofs);
+//   t_adjoint_->GetSubVector(vdofs, te_adjoint);
+
+//   // output vector
+//   elvect.SetSize(dim * dof);
+//   elvect = 0.0;
+
+//   // set integration rule
+//   const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), oa_ * el.GetOrder() + ob_);
+
+//   // loop over nodal coordinates (X_k)
+//   for (int m = 0; m < dim; m++) {
+//     for (int n = 0; n < dof; n++) {
+//       dX_dXk = 0.0;
+//       dX_dXk(m, n) = 1.0;
+
+//       dp_dXk = 0.0;
+//       // loop over integration points
+//       for (int i = 0; i < ir->GetNPoints(); i++) {
+//         // set integration point
+//         const ::mfem::IntegrationPoint &ip = ir->IntPoint(i);
+//         T.SetIntPoint(&ip);
+//         double w = ip.weight * T.Weight();
+
+//         // evaluate shape functions and their derivatives
+//         el.CalcShape(ip, N);
+//         el.CalcDShape(ip, dN);
+
+//         // compute derivative of Jacobian w.r.t. nodal coordinate
+//         mfem::Mult(dX_dXk, dN, dJ_dXk);
+
+//         // compute derivative of J^(-1)
+//         mfem::DenseMatrix JinvT = T.InverseJacobian();
+//         JinvT.Transpose();
+
+//         // compute derivative of integration weight w.r.t. X_k
+//         double dw_dXk = w * MatrixInnerProduct(JinvT, dJ_dXk);
+
+//         // add integration point's derivative contribution
+//         dp_dXk.Add(dw_dXk * Q_->Eval(T, ip), N);
+//       }
+//       elvect(n + m * dof) += InnerProduct(te_adjoint, dp_dXk);
+//     }
+//   }
+// }
+
+// void QuantityOfInterest::UpdateMesh(mfem::Vector const &U)
+// {
+//   ::mfem::Vector Xi = X0_;
+//   Xi += U;
+//   coord_fes_->GetParMesh()->SetNodes(Xi);
+
+//   coord_fes_->GetParMesh()->DeleteGeometricFactors();
+// }
+
+// void Diffusion_Solver::UpdateMesh(mfem::Vector const &U)
+// {
+//   ::mfem::Vector Xi = X0_;
+//   Xi += U;
+//   coord_fes_->GetParMesh()->SetNodes(Xi);
+
+//   coord_fes_->GetParMesh()->DeleteGeometricFactors();
+// }
+
+void NodeAwareTMOPQuality::UpdateMesh(mfem::Vector const &U)
+{
+  ::mfem::Vector Xi = X0_;
+  Xi += U;
+  coord_fes_->GetParMesh()->SetNodes(Xi);
+
+  coord_fes_->GetParMesh()->DeleteGeometricFactors();
+}
+
+// double QuantityOfInterest::EvalQoI()
+// {
+//   this->UpdateMesh(designVar);
+
+//   // make \nabla T vector coefficient
+
+//   ::mfem::ParGridFunction oneGridFunction = ::mfem::ParGridFunction(temp_fes_);
+//   oneGridFunction = 1.0;
+//   ConstantCoefficient one(1.0);
+//   BilinearFormIntegrator *integ = nullptr;
+//   ParGridFunction flux(coord_fes_);
+
+//   switch (qoiType_) {
+//   case 0:
+//      if( trueSolution_ == nullptr ){ mfem_error("true solution not set.");}
+//      ErrorCoefficient_ = std::make_shared<Error_QoI>(&solgf_, trueSolution_);
+//      break;
+//   case 1:
+//     if( trueSolutionGrad_ == nullptr ){ mfem_error("true solution not set.");}
+//     ErrorCoefficient_ = std::make_shared<H1Error_QoI>(&solgf_, trueSolutionGrad_);
+//     break;
+//   case 2:
+//     integ = new DiffusionIntegrator(one);
+//     solgf_.ComputeFlux(*integ, flux, false);
+//     trueSolutionGrad_ = new mfem::VectorGridFunctionCoefficient( &flux);
+//     ErrorCoefficient_ = std::make_shared<ZZError_QoI>(&solgf_, trueSolutionGrad_);
+//     break;
+//   default:
+//     std::cout << "Unknown Error Coeff: " << qoiType_ << std::endl;
+//   }
+
+//   ::mfem::ParGridFunction ErrorGF = ::mfem::ParGridFunction(temp_fes_);
+
+//   ::mfem::ParLinearForm scalarErrorForm(temp_fes_);
+//   LFErrorIntegrator *lfi = new LFErrorIntegrator;
+//   lfi->SetQoI(ErrorCoefficient_);
+//   lfi->SetIntRule(&mfem::IntRules.Get(temp_fes_->GetFE(0)->GetGeomType(), 8));
+//   scalarErrorForm.AddDomainIntegrator(lfi);
+//   scalarErrorForm.Assemble();
+
+//   return scalarErrorForm(oneGridFunction);
+// }
+
+// void QuantityOfInterest::EvalQoIGrad()
+// {
+//   this->UpdateMesh(designVar);
+
+//   ConstantCoefficient one(1.0);
+//   BilinearFormIntegrator *integ = nullptr;
+//   ParGridFunction flux(coord_fes_);
+
+//   switch (qoiType_) {
+//     case 0:
+//       if( trueSolution_ == nullptr ){ mfem_error("true solution not set.");}
+//       ErrorCoefficient_ = std::make_shared<Error_QoI>(&solgf_, trueSolution_);
+//       break;
+//     case 1:
+//       if( trueSolutionGrad_ == nullptr ){ mfem_error("true solution not set.");}
+//       ErrorCoefficient_ = std::make_shared<H1Error_QoI>(&solgf_, trueSolutionGrad_);
+//       break;
+//     case 2:
+//       integ = new DiffusionIntegrator(one);
+//       solgf_.ComputeFlux(*integ, flux, false);
+//       trueSolutionGrad_ = new mfem::VectorGridFunctionCoefficient( &flux);
+//       ErrorCoefficient_ = std::make_shared<ZZError_QoI>(&solgf_, trueSolutionGrad_);
+//       break;
+//     default:
+//       std::cout << "Unknown Error Coeff: " << qoiType_ << std::endl;
+//    }
+
+//   if(qoiType_ == QoIType::ZZ_ERROR)
+//   {
+//     // evaluate grad wrt temp
+//     {
+//         int nfe = temp_fes_->GetNE();
+
+
+//         Array<int> fdofs;
+//         DofTransformation *fdoftrans;
+//         Array<int> count(solgf_.Size()); count = 0;
+
+//         for (int i = 0; i < nfe; i++)
+//         {
+//           fdoftrans = temp_fes_->GetElementVDofs(i, fdofs);
+//           FiniteElementSpace::AdjustVDofs(fdofs);
+
+//           for (int j = 0; j < fdofs.Size(); j++)
+//           {
+//             count[fdofs[j]]++;
+//           }
+//         }
+
+//       ::mfem::ParLinearForm T_gradForm(temp_fes_);
+//       LFErrorDerivativeIntegrator_2 *lfi = new LFErrorDerivativeIntegrator_2(temp_fes_, count);
+//       lfi->SetQoI(ErrorCoefficient_);
+//       lfi->SetIntRule(&mfem::IntRules.Get(temp_fes_->GetFE(0)->GetGeomType(), 8));
+//       T_gradForm.AddDomainIntegrator(lfi);
+//       T_gradForm.Assemble();
+//       *dQdu_ = 0.0;
+//       dQdu_->Add( 1.0, T_gradForm);
+
+//       //solgf_.ComputeFlux(*integ, flux, false);
+//     }
+//       // evaluate grad wrt coord
+//     {
+//       LFNodeCoordinateSensitivityIntegrator *lfi = new LFNodeCoordinateSensitivityIntegrator;
+//       lfi->SetQoI(ErrorCoefficient_);
+//       lfi->SetIntRule(&mfem::IntRules.Get(coord_fes_->GetFE(0)->GetGeomType(), 8));
+
+//       ::mfem::ParLinearForm ud_gradForm(coord_fes_);
+//       ud_gradForm.AddDomainIntegrator(lfi);
+//       ud_gradForm.Assemble();
+//       *dQdx_ = 0.0;
+//       dQdx_->Add(1.0, ud_gradForm);
+//     }
+
+//   }
+//   else
+//   {
+//     // evaluate grad wrt temp
+//     {
+//       ::mfem::ParLinearForm T_gradForm(temp_fes_);
+//       LFErrorDerivativeIntegrator *lfi = new LFErrorDerivativeIntegrator;
+//       lfi->SetQoI(ErrorCoefficient_);
+//       lfi->SetIntRule(&mfem::IntRules.Get(temp_fes_->GetFE(0)->GetGeomType(), 8));
+//       T_gradForm.AddDomainIntegrator(lfi);
+//       T_gradForm.Assemble();
+//       *dQdu_ = 0.0;
+//       dQdu_->Add( 1.0, T_gradForm);
+//     }
+
+//     // evaluate grad wrt coord
+//     {
+//       LFNodeCoordinateSensitivityIntegrator *lfi = new LFNodeCoordinateSensitivityIntegrator;
+//       lfi->SetQoI(ErrorCoefficient_);
+//       lfi->SetIntRule(&mfem::IntRules.Get(coord_fes_->GetFE(0)->GetGeomType(), 8));
+
+//       ::mfem::ParLinearForm ud_gradForm(coord_fes_);
+//       ud_gradForm.AddDomainIntegrator(lfi);
+//       ud_gradForm.Assemble();
+//       *dQdx_ = 0.0;
+//       dQdx_->Add(1.0, ud_gradForm);
+//     }
+//   }
+// }
+
+double NodeAwareTMOPQuality::EvalQoI()
+{
+  this->UpdateMesh(designVar);
+
+  ::mfem::Vector Xi(X0_.Size()); Xi = X0_;
+  Xi += designVar;
+
+  // int targetId = 1;
+  // int metricId = 2;
+  int quadOrder = 8;
+
+  // Setup the mesh quality metric
+  // mfem::TMOP_QualityMetric *metric = nullptr;
+  // switch (metricId) {
+  //   // T-metrics
+  //   case 1:
+  //     metric = new mfem::TMOP_Metric_001;
+  //     break;
+  //   case 2:
+  //     metric = new mfem::TMOP_Metric_002;
+  //     break;
+  //   case 7:
+  //     metric = new mfem::TMOP_Metric_007;
+  //     break;
+  //   case 9:
+  //     metric = new mfem::TMOP_Metric_009;
+  //     break;
+  //   default:
+  //     std::cout << "Unknown metricId_: " << metricId << std::endl;
+  // }
+
+  // mfem::TargetConstructor::TargetType targetT;
+  // mfem::TargetConstructor *targetC = nullptr;
+  // switch (targetId) {
+  //   case 1:
+  //     targetT = mfem::TargetConstructor::IDEAL_SHAPE_UNIT_SIZE;
+  //     break;
+  //   case 2:
+  //     targetT = mfem::TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE;
+  //     break;
+  //   case 3:
+  //     targetT = mfem::TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE;
+  //     break;
+  //   default:
+  //     std::cout << "Unknown targetId: " << targetId << "\n";
+  // }
+  // if (nullptr == targetC) {
+  //   targetC = new mfem::TargetConstructor(targetT, MPI_COMM_WORLD);
+  // }
+  //targetC->SetNodes(X0_);
+
+  mfem::TMOP_Integrator *TMOPInteg = new mfem::TMOP_Integrator(metric_in, target_in);
+  // std::cout << metric_in << " " << target_in << " k10info\n";
+  // mfem::TMOP_Integrator *TMOPInteg = new mfem::TMOP_Integrator(metric, targetC);
+  mfem::IntegrationRules *irules = nullptr;
+  mfem::IntegrationRules IntRulesLo(0, mfem::Quadrature1D::GaussLobatto);
+  irules = &IntRulesLo;
+  TMOPInteg->SetIntegrationRules(*irules, quadOrder);
+
+  mfem::ParNonlinearForm a(coord_fes_);
+  a.AddDomainIntegrator(TMOPInteg);
+
+  double finalTMOPEnergy = a.GetParGridFunctionEnergy(Xi);
+
+  return finalTMOPEnergy;
+}
+
+void NodeAwareTMOPQuality::EvalQoIGrad()
+{
+  this->UpdateMesh(designVar);
+
+  ::mfem::ParGridFunction Xi(coord_fes_); Xi = X0_;
+  Xi += designVar;
+
+  int targetId = 1;
+  int metricId = 2;
+  int quadOrder = 8;
+
+  // // Setup the mesh quality metric
+  // mfem::TMOP_QualityMetric *metric = nullptr;
+  // switch (metricId) {
+  //   // T-metrics
+  //   case 1:
+  //     metric = new mfem::TMOP_Metric_001;
+  //     break;
+  //   case 2:
+  //     metric = new mfem::TMOP_Metric_002;
+  //     break;
+  //   case 7:
+  //     metric = new mfem::TMOP_Metric_007;
+  //     break;
+  //   case 9:
+  //     metric = new mfem::TMOP_Metric_009;
+  //     break;
+  //   default:
+  //     std::cout << "Unknown metricId_: " << metricId << std::endl;
+  // }
+
+  // mfem::TargetConstructor::TargetType targetT;
+  // mfem::TargetConstructor *targetC = nullptr;
+  // switch (targetId) {
+  //   case 1:
+  //     targetT = mfem::TargetConstructor::IDEAL_SHAPE_UNIT_SIZE;
+  //     break;
+  //   case 2:
+  //     targetT = mfem::TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE;
+  //     break;
+  //   case 3:
+  //     targetT = mfem::TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE;
+  //     break;
+  //   default:
+  //     std::cout << "Unknown targetId: " << targetId << "\n";
+  // }
+  // if (nullptr == targetC) {
+  //   targetC = new mfem::TargetConstructor(targetT, MPI_COMM_WORLD);
+  // }
+  //targetC->SetNodes(X0_);
+
+  mfem::TMOP_Integrator *TMOPInteg = new mfem::TMOP_Integrator(metric_in, target_in);
+  // mfem::TMOP_Integrator *TMOPInteg = new mfem::TMOP_Integrator(metric, targetC);
+  mfem::IntegrationRules *irules = nullptr;
+  mfem::IntegrationRules IntRulesLo(0, mfem::Quadrature1D::GaussLobatto);
+  irules = &IntRulesLo;
+  TMOPInteg->SetIntegrationRules(*irules, quadOrder);
+
+  mfem::ParNonlinearForm a(coord_fes_);
+  a.AddDomainIntegrator(TMOPInteg);
+
+  *dQdx_ = 0.0;
+
+  Xi.SetTrueVector();
+  mfem::Vector& trueXi = Xi.GetTrueVector();
+  ::mfem::ParGridFunction dQdx_GF(coord_fes_);
+  dQdx_GF.SetTrueVector();
+  mfem::Vector& truedQdx_GF = dQdx_GF.GetTrueVector();
+
+  a.Mult(trueXi, truedQdx_GF);
+
+  coord_fes_->GetRestrictionTransposeOperator()->Mult(truedQdx_GF, *dQdx_ );
+}
+
+// void Diffusion_Solver::FSolve()
+// {
+//   this->UpdateMesh(designVar);
+
+//   ::mfem::Array<int> ess_tdof_list(ess_tdof_list_);
+
+//   // assemble LHS matrix
+//   ::mfem::ConstantCoefficient kCoef(1.0);
+
+//   ::mfem::ParBilinearForm kForm(temp_fes_);
+//   ::mfem::ParLinearForm QForm(temp_fes_);
+//   kForm.AddDomainIntegrator(new ::mfem::DiffusionIntegrator(kCoef));
+
+//   QForm.AddDomainIntegrator(new ::mfem::DomainLFIntegrator(*QCoef_));
+
+//   kForm.Assemble();
+//   QForm.Assemble();
+
+//   // solve for temperature
+//   ::mfem::ParGridFunction &T = solgf;
+
+//   ::mfem::HypreParMatrix A;
+//   ::mfem::Vector X, B;
+//   kForm.FormLinearSystem(ess_tdof_list, T, QForm, A, X, B);
+
+//   ::mfem::HypreBoomerAMG amg(A);
+//   amg.SetPrintLevel(0);
+
+//   ::mfem::CGSolver cg(temp_fes_->GetParMesh()->GetComm());
+//   cg.SetRelTol(1e-10);
+//   cg.SetMaxIter(500);
+//   cg.SetPreconditioner(amg);
+//   cg.SetOperator(A);
+//   cg.Mult(B, X);
+
+//   kForm.RecoverFEMSolution(X, QForm, T);
+// }
+
+// void Diffusion_Solver::ASolve( mfem::Vector & rhs )
+// {
+//     // the nodal coordinates will default to the initial mesh
+//     this->UpdateMesh(designVar);
+
+//     ::mfem::Array<int> ess_tdof_list(ess_tdof_list_);
+
+//     // assemble LHS matrix
+//     ::mfem::ConstantCoefficient kCoef(1.0);
+
+//     ::mfem::ParBilinearForm kForm(temp_fes_);
+//     kForm.AddDomainIntegrator(new ::mfem::DiffusionIntegrator(kCoef));
+//     kForm.Assemble();
+
+//     // solve adjoint problem
+//     ::mfem::ParGridFunction adj_sol(temp_fes_);
+//     adj_sol = 0.0;
+
+//     ::mfem::HypreParMatrix A;
+//     ::mfem::Vector X, B;
+//     kForm.FormLinearSystem(ess_tdof_list, adj_sol, rhs, A, X, B);
+
+//     ::mfem::HypreBoomerAMG amg(A);
+//     amg.SetPrintLevel(0);
+
+//     ::mfem::CGSolver cg(temp_fes_->GetParMesh()->GetComm());
+//     cg.SetRelTol(1e-10);
+//     cg.SetMaxIter(500);
+//     cg.SetPreconditioner(amg);
+//     cg.SetOperator(A);
+//     cg.Mult(B, X);
+
+//     kForm.RecoverFEMSolution(X, rhs, adj_sol);
+
+//     ::mfem::ParLinearForm LHS_sensitivity(coord_fes_);
+//     LHS_sensitivity.AddDomainIntegrator(new ThermalConductivityShapeSensitivityIntegrator(kCoef, solgf, adj_sol));
+//     LHS_sensitivity.Assemble();
+
+//     ::mfem::ParLinearForm RHS_sensitivity(coord_fes_);
+//     RHS_sensitivity.AddDomainIntegrator(new ThermalHeatSourceShapeSensitivityIntegrator(*QCoef_, adj_sol));
+//     RHS_sensitivity.Assemble();
+
+//     *dQdx_ = 0.0;
+//     dQdx_->Add(-1.0, LHS_sensitivity);
+//     dQdx_->Add( 1.0, RHS_sensitivity);
+// }
+}

miniapps/meshing/mesh-optimizer_using_NLP.hpp

@@ -0,0 +1,709 @@
+#ifndef DIFFUSION_SOLVER_HPP
+#define DIFFUSION_SOLVER_HPP
+
+#include "mfem.hpp"
+
+
+namespace mfem{
+
+// void IdentityMatrix(int dim, mfem::DenseMatrix &I);
+
+// void Vectorize(const mfem::DenseMatrix &A, mfem::Vector &a);
+
+// double MatrixInnerProduct(const mfem::DenseMatrix &A, const mfem::DenseMatrix &B);
+
+// void ConjugationProduct(const mfem::DenseMatrix &A, const mfem::DenseMatrix &B, const mfem::DenseMatrix &C, mfem::DenseMatrix &D);
+
+// enum QoIType
+// {
+//   L2_ERROR,
+//   H1_ERROR,
+//   ZZ_ERROR
+// };
+
+// class QoIBaseCoefficient : public mfem::Coefficient {
+// public:
+//   QoIBaseCoefficient() {};
+
+//   virtual ~QoIBaseCoefficient() {};
+
+//   virtual const mfem::DenseMatrix &explicitSolutionDerivative(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
+//     = 0;
+
+//   virtual const mfem::DenseMatrix &explicitSolutionGradientDerivative(mfem::ElementTransformation &T,
+//       const mfem::IntegrationPoint &ip)
+//     = 0;
+
+//   virtual const mfem::DenseMatrix &gradTimesexplicitSolutionGradientDerivative(mfem::ElementTransformation &T,
+//       const mfem::IntegrationPoint &ip)
+//     = 0;
+
+//   virtual const mfem::DenseMatrix &explicitShapeDerivative(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) = 0;
+// private:
+// };
+
+// class Error_QoI : public QoIBaseCoefficient
+// {
+// public:
+//   Error_QoI(mfem::ParGridFunction * solutionField, mfem::Coefficient * trueSolution)
+//     : solutionField_(solutionField), trueSolution_(trueSolution)
+//   {};
+
+//   ~Error_QoI() {};
+
+//   double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
+//   {
+//     double fieldVal = solutionField_->GetValue( T, ip );
+//     double trueVal = trueSolution_->Eval( T, ip );
+
+//     double squaredError = std::pow( fieldVal-trueVal, 2.0);
+
+//     return squaredError;
+//   };
+
+//   const mfem::DenseMatrix &explicitSolutionDerivative(mfem::ElementTransformation & T, const mfem::IntegrationPoint & ip)
+//   {
+//     dtheta_dU.SetSize(1);
+
+//     double val = 2.0* (solutionField_->GetValue( T, ip ) - trueSolution_->Eval( T, ip ));
+
+//     double & matVal = dtheta_dU.Elem(0,0);
+//     matVal = val;
+//     return dtheta_dU;
+//   };
+
+//   const mfem::DenseMatrix &explicitSolutionGradientDerivative(mfem::ElementTransformation & /*T*/,
+//       const mfem::IntegrationPoint & /*ip*/)
+//   {
+//     dtheta_dGradU.SetSize(1, Dim_);
+//     dtheta_dGradU = 0.0;
+
+//     return dtheta_dGradU;
+//   };
+
+//   const mfem::DenseMatrix &explicitShapeDerivative(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
+//   {
+//     dtheta_dX.SetSize(1, Dim_);
+//     dtheta_dX = 0.0;
+
+//     return dtheta_dX;
+//   };
+
+//   virtual const mfem::DenseMatrix &gradTimesexplicitSolutionGradientDerivative(mfem::ElementTransformation &T,
+//       const mfem::IntegrationPoint &ip)
+//   {
+//     dtheta_dX.SetSize(Dim_, Dim_);
+//     dtheta_dX = 0.0;
+
+//     return dtheta_dX;
+//   };
+
+// private:
+
+//   mfem::ParGridFunction * solutionField_;
+//   mfem::Coefficient * trueSolution_;
+
+//   int Dim_ = 2;
+
+//   double theta = 0.0;
+//   mfem::DenseMatrix dtheta_dX;
+//   mfem::DenseMatrix dtheta_dU;
+//   mfem::DenseMatrix dtheta_dGradU;
+// };
+
+// class H1Error_QoI : public QoIBaseCoefficient {
+// public:
+//   H1Error_QoI(mfem::ParGridFunction * solutionField, mfem::VectorCoefficient * trueSolution)
+//     : solutionField_(solutionField), trueSolution_(trueSolution)
+//   {};
+
+//   ~H1Error_QoI() {};
+
+//   double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
+//   {
+//     mfem::Vector grad;
+//     mfem::Vector trueGrad;
+//     trueSolution_->Eval (trueGrad, T, ip);
+//     solutionField_->GetGradient (T, grad);
+
+//     grad -= trueGrad;
+
+//     double val = grad.Norml2();
+//     val = val * val;
+
+//     return val;
+//   };
+
+//   const mfem::DenseMatrix &explicitSolutionDerivative(mfem::ElementTransformation & /*T*/, const mfem::IntegrationPoint & /*ip*/)
+//   {
+//     dtheta_dU.SetSize(1);
+//     dtheta_dU = 0.0;
+
+//     return dtheta_dU;
+//   };
+
+//   const mfem::DenseMatrix &explicitSolutionGradientDerivative(mfem::ElementTransformation & T,
+//       const mfem::IntegrationPoint & ip)
+//   {
+//     mfem::Vector grad(Dim_);
+//     mfem::Vector trueGrad(Dim_);
+//     mfem::Vector gradMinusTrueGrad(Dim_);
+//     trueSolution_->Eval (trueGrad, T, ip);
+//     solutionField_->GetGradient (T, grad);
+//     gradMinusTrueGrad = grad;
+//     gradMinusTrueGrad -= trueGrad;
+//     gradMinusTrueGrad *= 2.0;
+
+//     dtheta_dGradU.SetSize(1, Dim_);
+//     dtheta_dGradU = 0.0;
+
+//     dtheta_dGradU(0,0) = gradMinusTrueGrad[0];
+//     dtheta_dGradU(0,1) = gradMinusTrueGrad[1];
+
+
+//     return dtheta_dGradU;
+//   };
+
+//   const mfem::DenseMatrix &explicitShapeDerivative(mfem::ElementTransformation & /*T*/, const mfem::IntegrationPoint &/*ip*/)
+//   {
+//     dtheta_dX.SetSize(1, Dim_);
+//     dtheta_dX = 0.0;
+
+//     return dtheta_dX;
+//   };
+
+//   virtual const mfem::DenseMatrix &gradTimesexplicitSolutionGradientDerivative(mfem::ElementTransformation &T,
+//       const mfem::IntegrationPoint &ip)
+//   {
+//     mfem::Vector grad(Dim_);
+//     mfem::Vector trueGrad(Dim_);
+//     mfem::Vector gradMinusTrueGrad(Dim_);
+//     trueSolution_->Eval (trueGrad, T, ip);
+//     solutionField_->GetGradient (T, grad);
+//     gradMinusTrueGrad = grad;
+//     gradMinusTrueGrad -= trueGrad;
+//     gradMinusTrueGrad *= 2.0;
+
+//     dUXdtheta_dGradU.SetSize(Dim_, Dim_);
+//     dUXdtheta_dGradU = 0.0;
+
+//     dUXdtheta_dGradU(0,0) =  grad[0] * gradMinusTrueGrad[0];
+//     dUXdtheta_dGradU(1,0) =  grad[1] * gradMinusTrueGrad[0];
+//     dUXdtheta_dGradU(0,1) =  grad[0] * gradMinusTrueGrad[1];
+//     dUXdtheta_dGradU(1,1) =  grad[1] * gradMinusTrueGrad[1];
+
+//     dUXdtheta_dGradU.Transpose();
+
+
+//     return dUXdtheta_dGradU;
+//   };
+// private:
+
+
+//   mfem::ParGridFunction * solutionField_;
+//   mfem::VectorCoefficient * trueSolution_;
+
+//   int Dim_ = 2;
+
+//   double theta = 0.0;
+//   mfem::DenseMatrix dtheta_dX;
+//   mfem::DenseMatrix dtheta_dU;
+//   mfem::DenseMatrix dtheta_dGradU;
+//   mfem::DenseMatrix dUXdtheta_dGradU;
+// };
+
+// class ZZError_QoI : public QoIBaseCoefficient {
+// public:
+//   ZZError_QoI(mfem::ParGridFunction * solutionField, mfem::VectorCoefficient * trueSolution)
+//     : solutionField_(solutionField), trueSolution_(trueSolution)
+//   {};
+
+//   ~ZZError_QoI() {};
+
+//   double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip)
+//   {
+//     mfem::Vector grad;
+//     mfem::Vector trueGrad;
+//     trueSolution_->Eval (trueGrad, T, ip);
+//     solutionField_->GetGradient (T, grad);
+
+//     grad -= trueGrad;
+
+//     double val = grad.Norml2();
+//     val = val * val;
+
+//     return val;
+//   };
+
+//   const mfem::DenseMatrix &explicitSolutionDerivative(mfem::ElementTransformation & /*T*/, const mfem::IntegrationPoint & /*ip*/)
+//   {
+//     dtheta_dU.SetSize(1);
+//     dtheta_dU = 0.0;
+
+//     return dtheta_dU;
+//   };
+
+//   const mfem::DenseMatrix &explicitSolutionGradientDerivative(mfem::ElementTransformation & T,
+//       const mfem::IntegrationPoint & ip)
+//   {
+//     mfem::Vector grad(Dim_);
+//     mfem::Vector trueGrad(Dim_);
+//     mfem::Vector gradMinusTrueGrad(Dim_);
+//     trueSolution_->Eval (trueGrad, T, ip);
+//     solutionField_->GetGradient (T, grad);
+//     gradMinusTrueGrad = grad;
+//     gradMinusTrueGrad -= trueGrad;
+//     gradMinusTrueGrad *= 2.0;
+
+//     dtheta_dGradU.SetSize(1, Dim_);
+//     dtheta_dGradU = 0.0;
+
+//     dtheta_dGradU(0,0) = gradMinusTrueGrad[0];
+//     dtheta_dGradU(0,1) = gradMinusTrueGrad[1];
+
+
+//     return dtheta_dGradU;
+//   };
+
+//   const mfem::DenseMatrix &explicitShapeDerivative(mfem::ElementTransformation & /*T*/, const mfem::IntegrationPoint &/*ip*/)
+//   {
+//     dtheta_dX.SetSize(1, Dim_);
+//     dtheta_dX = 0.0;
+
+//     return dtheta_dX;
+//   };
+
+//   virtual const mfem::DenseMatrix &gradTimesexplicitSolutionGradientDerivative(mfem::ElementTransformation &T,
+//       const mfem::IntegrationPoint &ip)
+//   {
+//     mfem::Vector grad(Dim_);
+//     mfem::Vector trueGrad(Dim_);
+//     mfem::Vector gradMinusTrueGrad(Dim_);
+//     trueSolution_->Eval (trueGrad, T, ip);
+//     solutionField_->GetGradient (T, grad);
+//     gradMinusTrueGrad = grad;
+//     gradMinusTrueGrad -= trueGrad;
+//     gradMinusTrueGrad *= 2.0;
+
+//     dUXdtheta_dGradU.SetSize(Dim_, Dim_);
+//     dUXdtheta_dGradU = 0.0;
+
+//     dUXdtheta_dGradU(0,0) =  grad[0] * gradMinusTrueGrad[0];
+//     dUXdtheta_dGradU(1,0) =  grad[1] * gradMinusTrueGrad[0];
+//     dUXdtheta_dGradU(0,1) =  grad[0] * gradMinusTrueGrad[1];
+//     dUXdtheta_dGradU(1,1) =  grad[1] * gradMinusTrueGrad[1];
+
+//     dUXdtheta_dGradU.Transpose();
+
+
+//     return dUXdtheta_dGradU;
+//   };
+// private:
+
+
+//   mfem::ParGridFunction * solutionField_;
+//   mfem::VectorCoefficient * trueSolution_;
+
+//   int Dim_ = 2;
+
+//   double theta = 0.0;
+//   mfem::DenseMatrix dtheta_dX;
+//   mfem::DenseMatrix dtheta_dU;
+//   mfem::DenseMatrix dtheta_dGradU;
+//   mfem::DenseMatrix dUXdtheta_dGradU;
+// };
+
+// //   const mfem::DenseMatrix &explicitShapeDerivative(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip);
+// // private:
+
+// //   mfem::ParGridFunction * solutionField_;
+// //   mfem::Coefficient * trueSolution_;
+
+// //   int Dim_ = 2;
+
+// //   double theta = 0.0;
+// //   mfem::DenseMatrix dtheta_dX;
+// //   mfem::DenseMatrix dtheta_dU;
+// //   mfem::DenseMatrix dtheta_dGradU;
+// // };
+
+// // class ExplicitPhysicsAware_QoI : public QoIBaseCoefficient {
+// // public:
+// //   ExplicitPhysicsAware_QoI(mfem::ParGridFunction * solutionField, mfem::Coefficient * trueSolution)
+// //     : solutionField_(solutionField), trueSolution_(trueSolution)
+// //   {};
+
+// //   ~ExplicitPhysicsAware_QoI() {};
+
+// //   double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip);
+
+// //   const mfem::DenseMatrix &explicitSolutionDerivative(mfem::ElementTransformation & T, const mfem::IntegrationPoint & ip)
+// //   {
+// //     dtheta_dU.SetSize(1);
+
+// //     double val = 2.0* (solutionField_->GetValue( T, ip ) - trueSolution_->Eval( T, ip ));
+
+// //     double & matVal = dtheta_dU.Elem(0,0);
+// //     matVal = val;
+// //     return dtheta_dU;
+// //   };
+
+// //   const mfem::DenseMatrix &explicitSolutionGradientDerivative(mfem::ElementTransformation & /*T*/,
+// //       const mfem::IntegrationPoint & /*ip*/)
+// //   {
+// //     return dtheta_dGradU;
+// //   };
+
+// //   const mfem::DenseMatrix &explicitShapeDerivative(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip);
+// // private:
+
+// //   mfem::ParGridFunction * solutionField_;
+// //   mfem::Coefficient * trueSolution_;
+
+// //   int Dim_ = 2;
+
+// //   double theta = 0.0;
+// //   mfem::DenseMatrix dtheta_dX;
+// //   mfem::DenseMatrix dtheta_dU;
+// //   mfem::DenseMatrix dtheta_dGradU;
+// // };
+
+// class LFNodeCoordinateSensitivityIntegrator : public mfem::LinearFormIntegrator {
+// public:
+//   LFNodeCoordinateSensitivityIntegrator(int Index1 = INT_MAX, int Index2 = INT_MAX,
+//                                         int IntegrationOrder = INT_MAX);
+//   ~LFNodeCoordinateSensitivityIntegrator() {};
+//   void AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T, mfem::Vector &elvect);
+//   void SetQoI(std::shared_ptr<QoIBaseCoefficient> QoI) { QoI_ = QoI; };
+// private:
+//   std::shared_ptr<QoIBaseCoefficient> QoIFactoryFunction(const int dim);
+
+//   const int Index1_;
+//   const int Index2_;
+//   const int IntegrationOrder_;
+
+//   std::shared_ptr<QoIBaseCoefficient> QoI_ = nullptr;
+// };
+
+// class LFErrorIntegrator : public mfem::LinearFormIntegrator {
+// public:
+//   LFErrorIntegrator( int IntegrationOrder = INT_MAX);
+//   ~LFErrorIntegrator() {};
+//   void AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T, mfem::Vector &elvect);
+//   void SetQoI(std::shared_ptr<QoIBaseCoefficient> QoI) { QoI_ = QoI; };
+// private:
+//   std::shared_ptr<QoIBaseCoefficient> QoIFactoryFunction(const int dim);
+
+//   const int IntegrationOrder_;
+
+//   std::shared_ptr<QoIBaseCoefficient> QoI_ = nullptr;
+// };
+
+// class LFErrorDerivativeIntegrator : public mfem::LinearFormIntegrator {
+// public:
+//   LFErrorDerivativeIntegrator( int IntegrationOrder = INT_MAX);
+//   ~LFErrorDerivativeIntegrator() {};
+//   void AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T, mfem::Vector &elvect);
+//   void SetQoI(std::shared_ptr<QoIBaseCoefficient> QoI) { QoI_ = QoI; };
+// private:
+//   std::shared_ptr<QoIBaseCoefficient> QoIFactoryFunction(const int dim);
+
+//   const int IntegrationOrder_;
+
+//   std::shared_ptr<QoIBaseCoefficient> QoI_ = nullptr;
+// };
+
+// class LFErrorDerivativeIntegrator_2 : public mfem::LinearFormIntegrator {
+// public:
+//   LFErrorDerivativeIntegrator_2( ParFiniteElementSpace * fespace, Array<int> count, int IntegrationOrder = INT_MAX);
+//   ~LFErrorDerivativeIntegrator_2() {};
+//   void AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T, mfem::Vector &elvect);
+//   void SetQoI(std::shared_ptr<QoIBaseCoefficient> QoI) { QoI_ = QoI; };
+// private:
+//   std::shared_ptr<QoIBaseCoefficient> QoIFactoryFunction(const int dim);
+
+//   ParFiniteElementSpace * fespace_ = nullptr;
+//   Array<int> count_;
+//   const int IntegrationOrder_;
+
+//   std::shared_ptr<QoIBaseCoefficient> QoI_ = nullptr;
+// };
+
+// class ThermalConductivityShapeSensitivityIntegrator : public mfem::LinearFormIntegrator {
+// public:
+//   ThermalConductivityShapeSensitivityIntegrator(mfem::Coefficient &conductivity, const mfem::ParGridFunction &t_primal,
+//       const mfem::ParGridFunction &t_adjoint);
+//   void AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T, mfem::Vector &elvect);
+// private:
+//   mfem::Coefficient *k_;
+//   const mfem::ParGridFunction *t_primal_;
+//   const mfem::ParGridFunction *t_adjoint_;
+// };
+
+// class ThermalHeatSourceShapeSensitivityIntegrator : public mfem::LinearFormIntegrator {
+// public:
+//   ThermalHeatSourceShapeSensitivityIntegrator(mfem::Coefficient &heatSource, const mfem::ParGridFunction &t_adjoint, int oa = 2,
+//       int ob = 2);
+//   void AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &T, mfem::Vector &elvect);
+// private:
+//   mfem::Coefficient *Q_;
+//   const mfem::ParGridFunction *t_adjoint_;
+//   int oa_, ob_;
+// };
+
+// class QuantityOfInterest
+// {
+// public:
+//     QuantityOfInterest(mfem::ParMesh* mesh_, enum QoIType qoiType, int order_=1)
+//     : pmesh(mesh_), qoiType_(qoiType)
+//     {
+//         int dim=pmesh->Dimension();
+
+//         pmesh->GetNodes(X0_);
+
+//         fec = new H1_FECollection(order_,dim);
+//         temp_fes_ = new ParFiniteElementSpace(pmesh,fec);
+//         coord_fes_ = new ParFiniteElementSpace(pmesh,fec,dim);
+
+//         solgf_.SetSpace(temp_fes_);
+
+//         dQdu_ = new mfem::ParLinearForm(temp_fes_);
+//         dQdx_ = new mfem::ParLinearForm(coord_fes_);
+//     }
+
+//     ~QuantityOfInterest()
+//     {
+//         delete temp_fes_;
+//         delete coord_fes_;
+//         delete fec;
+
+//         delete dQdu_;
+//         delete dQdx_;
+//     }
+
+//     void setTrueSolCoeff( mfem::Coefficient * trueSolution ){ trueSolution_ = trueSolution; };
+//     void setTrueSolGradCoeff( mfem::VectorCoefficient * trueSolutionGrad ){ trueSolutionGrad_ = trueSolutionGrad; };
+//     void SetDesign( mfem::Vector & design){ designVar = design; };
+//     void SetDiscreteSol( mfem::ParGridFunction & sol){ solgf_ = sol; };
+//     void UpdateMesh(mfem::Vector const &U);
+//     double EvalQoI();
+//     void EvalQoIGrad();
+//     mfem::ParLinearForm * GetDQDu(){ return dQdu_; };
+//     mfem::ParLinearForm * GetDQDx(){ return dQdx_; };
+
+// private:
+//     mfem::Coefficient * trueSolution_ = nullptr;
+//     mfem::VectorCoefficient * trueSolutionGrad_ = nullptr;
+
+//     mfem::ParMesh* pmesh;
+//     enum QoIType qoiType_;
+
+//     mfem::Vector X0_;
+//     mfem::Vector designVar;
+
+//     mfem::FiniteElementCollection *fec;
+//     mfem::ParFiniteElementSpace	  *temp_fes_;
+//     mfem::ParFiniteElementSpace	  *coord_fes_;
+
+//     mfem::ParLinearForm * dQdu_;
+//     mfem::ParLinearForm * dQdx_;
+
+//     mfem::ParGridFunction solgf_;
+
+//     std::shared_ptr<QoIBaseCoefficient> ErrorCoefficient_ = nullptr;
+// };
+
+
+
+class NodeAwareTMOPQuality
+{
+public:
+    NodeAwareTMOPQuality(mfem::ParMesh* mesh_, int order_, TMOP_QualityMetric *metric, TargetConstructor *target_c)
+    {
+        pmesh=mesh_;
+        int dim=pmesh->Dimension();
+
+        fec = new H1_FECollection(order_,dim);
+        coord_fes_ = new ParFiniteElementSpace(pmesh,fec,dim);
+
+        X0_.SetSpace(coord_fes_);
+        designVar.SetSpace(coord_fes_);
+        mfem::Vector tempX0_;
+        pmesh->GetNodes(tempX0_);
+
+        X0_ = tempX0_;
+
+        dQdx_ = new mfem::ParLinearForm(coord_fes_);
+        metric_in = metric;
+        target_in = target_c;
+    }
+
+    ~NodeAwareTMOPQuality()
+    {
+    }
+
+    void UpdateMesh(mfem::Vector const &U);
+    double EvalQoI();
+    void EvalQoIGrad();
+
+    mfem::ParLinearForm * GetDQDx(){ return dQdx_; };
+
+    void SetDesign( mfem::ParGridFunction & design){ designVar = design; };
+
+    private:
+
+    mfem::ParMesh* pmesh;
+    mfem::ParGridFunction X0_;
+    mfem::ParGridFunction designVar;
+
+    mfem::FiniteElementCollection *fec;
+    mfem::ParFiniteElementSpace	  *coord_fes_;
+
+    mfem::ParLinearForm * dQdx_;
+    TMOP_QualityMetric *metric_in = nullptr;
+    TargetConstructor *target_in = nullptr;
+};
+
+// class Diffusion_Solver
+// {
+// public:
+//     Diffusion_Solver(mfem::ParMesh* mesh_, std::vector<std::pair<int, double>> ess_bdr, int order_=2)
+//     {
+//         pmesh=mesh_;
+//         int dim=pmesh->Dimension();
+
+//         pmesh->GetNodes(X0_);
+
+//         fec = new H1_FECollection(order_,dim);
+//         temp_fes_ = new ParFiniteElementSpace(pmesh,fec);
+//         coord_fes_ = new ParFiniteElementSpace(pmesh,fec,dim);
+
+//         sol.SetSize(temp_fes_->GetTrueVSize()); sol=0.0;
+//         rhs.SetSize(temp_fes_->GetTrueVSize()); rhs=0.0;
+//         adj.SetSize(temp_fes_->GetTrueVSize()); adj=0.0;
+
+//         solgf.SetSpace(temp_fes_);
+//         adjgf.SetSpace(temp_fes_);
+
+//         dQdu_ = new mfem::ParLinearForm(temp_fes_);
+//         dQdx_ = new mfem::ParLinearForm(coord_fes_);
+
+//         SetLinearSolver();
+
+//         // store list of essential dofs
+//         int maxAttribute = pmesh->bdr_attributes.Max();
+//         ::mfem::Array<int> bdr_attr_is_ess(maxAttribute);
+//         ess_tdof_list_.DeleteAll();
+//         ::mfem::Vector ess_bc(temp_fes_->GetTrueVSize());
+//         ess_bc = 0.0;
+
+//         // loop over input attribute, value pairs
+//         for (const auto &bc: ess_bdr)
+//         {
+//             int attribute = bc.first;
+
+//             // get dofs associated with this attribute, component pair
+//             bdr_attr_is_ess = 0;
+//             bdr_attr_is_ess[attribute - 1] = 1; // mfem attributes 1-indexed, arrays 0-indexed
+//             ::mfem::Array<int> temp_tdofs;
+//             temp_fes_->GetEssentialTrueDofs(bdr_attr_is_ess, temp_tdofs);
+
+//             // append to global dof list
+//             ess_tdof_list_.Append(temp_tdofs);
+
+//             // set value in grid function
+//             double value = bc.second;
+//             ess_bc.SetSubVector(temp_tdofs, value);
+//         }
+//         bcGridFunc_.SetSpace(temp_fes_);
+//         bcGridFunc_.SetFromTrueDofs(ess_bc);
+//     }
+
+//     ~Diffusion_Solver(){
+//         delete temp_fes_;
+//         delete coord_fes_;
+//         delete fec;
+
+//         delete dQdu_;
+//         delete dQdx_;
+//     }
+
+//     void UpdateMesh(mfem::Vector const &U);
+
+//     double Eval_QoI();
+
+//     void Eval_QoI_Grad();
+
+//     /// Set the Linear Solver
+//     void SetLinearSolver(double rtol=1e-8, double atol=1e-12, int miter=2000)
+//     {
+//         linear_rtol=rtol;
+//         linear_atol=atol;
+//         linear_iter=miter;
+//     }
+
+//     /// Solves the forward problem.
+//     void FSolve();
+
+//     void ASolve( mfem::Vector & rhs );
+
+//     void SetDesign( mfem::Vector & design)
+//     {
+//         designVar = design;
+//     };
+
+//     void SetManufacturedSolution( mfem::Coefficient * QCoef )
+//     {
+//       QCoef_ = QCoef;
+//     }
+
+//     /// Returns the solution
+//     mfem::ParGridFunction& GetSolution(){return solgf;}
+
+//     /// Returns the solution vector.
+//     mfem::Vector& GetSol(){return sol;}
+
+//     /// Returns the adjoint solution vector.
+//     mfem::Vector& GetAdj(){return adj;}
+
+//     mfem::ParLinearForm * GetImplicitDqDx(){ return dQdx_; };
+
+// private:
+//     mfem::ParMesh* pmesh;
+
+//     mfem::Vector X0_;
+//     mfem::Vector designVar;
+
+//     //solution true vector
+//     mfem::Vector sol;
+//     mfem::Vector adj;
+//     mfem::Vector rhs;
+//     mfem::ParGridFunction solgf;
+//     mfem::ParGridFunction adjgf;
+//     mfem::ParGridFunction bcGridFunc_;
+
+//     mfem::ParLinearForm * dQdu_;
+//     mfem::ParLinearForm * dQdx_;
+
+//     mfem::FiniteElementCollection *fec;
+//     mfem::ParFiniteElementSpace	  *temp_fes_;
+//     mfem::ParFiniteElementSpace	  *coord_fes_;
+
+//     //Linear solver parameters
+//     double linear_rtol;
+//     double linear_atol;
+//     int linear_iter;
+
+//     int print_level = 1;
+
+//     // holds NBC in coefficient form
+//     std::map<int, mfem::Coefficient*> ncc;
+
+//     mfem::Array<int> ess_tdof_list_;
+
+//     mfem::Coefficient * QCoef_ = nullptr;
+// };
+
+}
+
+#endif

miniapps/meshing/pmesh-optimizer.cpp

@@ -36,8 +36,8 @@
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -bnd -qt 1 -qo 8
 //   Adapted analytic size+orientation:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 200 -bnd -qt 1 -qo 8
-//   Adapted analytic shape+orientation:
-//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
+//   Adapted analytic shape+orientation (AD):
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8
 //
 //   Adapted analytic shape and/or size with hr-adaptivity:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 9  -ni 50 -li 20 -hmid 55 -mid 7 -hr
@@ -59,8 +59,8 @@
 //   Adapted discrete size+aspect_ratio:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
-//   Adapted discrete size+orientation:
-//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -fd -nor
+//   Adapted discrete size+orientation (AD):
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -nor -ex
 //   Adapted discrete aspect ratio (3D):
 //     mpirun -np 4 pmesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20 -bnd -qt 1 -qo 8
 //
@@ -73,6 +73,8 @@
 //
 //   Blade shape:
 //     mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8
+//   Blade shape (AD):
+//     mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 11 -ex -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8
 //     (requires CUDA):
 //   * mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8 -d cuda
 //   Blade shape with FD-based solver:
@@ -107,6 +109,19 @@
 //   Kershaw transformation can be imposed using the transformation ('t') feature in the mesh-explorer miniapp.
 //   * mpirun - np 6 pmesh-optimizer -m kershaw-24x24x24.mesh -mid 303 -tid 1 -bnd -ni 100 -art 1 -ls 3 -qo 8 -li 40 -o 2 -qo 8 -ker -pa
 
+//    AD-based runs
+//    mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8
+
+// COMPARISON WITH MMA
+/////// BLADE
+// Blade with Newton
+// mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 1000 -ls 3 -bnd -qt 1 -qo 8 -vl 2 -st 0
+// Blade with LBFGS
+// mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 1000 -ls 3 -bnd -qt 1 -qo 8 -vl 2 -st 1
+// Blade with MMA
+// mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 1000 -ls 3 -bnd -qt 1 -qo 8 -vl 2 -mma -ch 1e-3
+// grep -i "Energy decrease:" comp_newton.out | awk '{print $2}'
+
 #include "mfem.hpp"
 #include "../common/mfem-common.hpp"
 #include <iostream>
@@ -164,6 +179,8 @@ int main (int argc, char *argv[])
    int mesh_node_ordering = 0;
    int barrier_type       = 0;
    int worst_case_type    = 0;
+   bool mma               = false;
+  double max_ch=0.0001; //max design change
 
    // 2. Parse command-line options.
    OptionsParser args(argc, argv);
@@ -320,6 +337,11 @@ int main (int argc, char *argv[])
                   "0 - None,"
                   "1 - Beta,"
                   "2 - PMean.");
+   args.AddOption(&mma, "-mma", "--mma", "-no-mma",
+                  "--no-mma",
+                  "Enable or disable mma solver.");
+   args.AddOption(&max_ch, "-ch", "--max-ch",
+                  "max node movement");
 
    args.Parse();
    if (!args.Good())
@@ -378,6 +400,10 @@ int main (int argc, char *argv[])
    //    changing x automatically changes the shapes of the mesh elements.
    ParGridFunction x(pfespace);
    pmesh->SetNodalGridFunction(&x);
+  ParGridFunction gridfuncLSBoundIndicator(pfespace);
+  gridfuncLSBoundIndicator = 0.0;
+  ParGridFunction gridfuncOptVar(pfespace);
+  gridfuncOptVar = 0.0;
 
    // 8. Define a vector representing the minimal local mesh size in the mesh
    //    nodes. We index the nodes using the scalar version of the degrees of
@@ -499,7 +525,7 @@ int main (int argc, char *argv[])
       // A-metrics
       case 11: metric = new TMOP_AMetric_011; break;
       case 36: metric = new TMOP_AMetric_036; break;
-      case 107: metric = new TMOP_AMetric_107a; break;
+      case 107: metric = new TMOP_AMetric_107; break;
       case 126: metric = new TMOP_AMetric_126(0.9); break;
       default:
          if (myid == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
@@ -1056,12 +1082,18 @@ int main (int argc, char *argv[])
          if (attr == 1) // Fix x components.
          {
             for (int j = 0; j < nd; j++)
-            { ess_vdofs[n++] = vdofs[j]; }
+            {
+               ess_vdofs[n++] = vdofs[j];
+               if (mma) { gridfuncLSBoundIndicator[ vdofs[j] ] = 1.0; }
+            }
          }
          else if (attr == 2) // Fix y components.
          {
             for (int j = 0; j < nd; j++)
-            { ess_vdofs[n++] = vdofs[j+nd]; }
+            {
+               ess_vdofs[n++] = vdofs[j+nd];
+               if (mma) { gridfuncLSBoundIndicator[ vdofs[j+nd] ] = 1.0; }
+            }
          }
          else if (attr == 3) // Fix z components.
          {
@@ -1076,6 +1108,9 @@ int main (int argc, char *argv[])
       }
       a.SetEssentialVDofs(ess_vdofs);
    }
+  gridfuncLSBoundIndicator.SetTrueVector();
+  gridfuncOptVar.SetTrueVector();
+  Vector & trueOptvar = gridfuncOptVar.GetTrueVector();
 
    // As we use the inexact Newton method to solve the resulting nonlinear
    // system, here we setup the linear solver for the system's Jacobian.
@@ -1138,43 +1173,87 @@ int main (int argc, char *argv[])
    //
    const IntegrationRule &ir =
       irules->Get(pfespace->GetFE(0)->GetGeomType(), quad_order);
-   TMOPNewtonSolver solver(pfespace->GetComm(), ir, solver_type);
-   // Provide all integration rules in case of a mixed mesh.
-   solver.SetIntegrationRules(*irules, quad_order);
-   // Specify linear solver when we use a Newton-based solver.
-   if (solver_type == 0) { solver.SetPreconditioner(*S); }
-   // For untangling, the solver will update the min det(T) values.
-   solver.SetMinDetPtr(&min_detJ);
-   solver.SetMaxIter(solver_iter);
-   solver.SetRelTol(solver_rtol);
-   solver.SetAbsTol(0.0);
-   if (solver_art_type > 0)
+   if (mma)
    {
-      solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
+      TMOP_MMA *tmma = new TMOP_MMA(MPI_COMM_WORLD, trueOptvar.Size(), 0,
+                                    trueOptvar, ir);
+      x.SetTrueVector();
+      double init_energy = a.GetParGridFunctionEnergy(x);
+      IterativeSolver::PrintLevel newton_print;
+      if (verbosity_level > 0)
+      { newton_print.Errors().Warnings().Iterations(); }
+      tmma->SetPrintLevel(newton_print);
+      // newton_print.Errors().Warnings().Iterations();
+      // set the TMOP Integrator
+      tmma->SetOperator(a);
+      // Set change limits on dx
+      tmma->SetUpperBound(max_ch);
+      tmma->SetLowerBound(max_ch);
+      // Set true vector so that it can be zeroed out
+      {
+         Vector & trueBounds = gridfuncLSBoundIndicator.GetTrueVector();
+         tmma->SetTrueDofs(trueBounds);
+      }
+      // Set QoI and Solver and weight
+      // tmma->SetQuantityOfInterest(&QoIEvaluator);
+      // tmma->SetDiffusionSolver(&solver);
+      // tmma->SetQoIWeight(weight_1);
+
+      // Set max # iterations
+      tmma->SetMaxIter(solver_iter);
+
+      tmma->Mult(x.GetTrueVector());
+      x.SetFromTrueVector();
+
+      // Visualize the mesh displacement.
+      if (visualization)
+      {
+         x0 -= x;
+         socketstream vis;
+         common::VisualizeField(vis, "localhost", 19916, x0,
+                                 "Displacements", 600, 000, 500, 500, "jRmclAppppppppppppp");
+      }
    }
-   // Level of output.
-   IterativeSolver::PrintLevel newton_print;
-   if (verbosity_level > 0)
-   { newton_print.Errors().Warnings().Iterations(); }
-   solver.SetPrintLevel(newton_print);
-   // hr-adaptivity solver.
-   // If hr-adaptivity is disabled, r-adaptivity is done once using the
-   // TMOPNewtonSolver.
-   // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
-   // "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity.
-   // The solver terminates if an h-adaptivity iteration does not modify
-   // any element in the mesh.
-   TMOPHRSolver hr_solver(*pmesh, a, solver,
-                          x, move_bnd, hradaptivity,
-                          mesh_poly_deg, h_metric_id,
-                          n_hr_iter, n_h_iter);
-   hr_solver.AddGridFunctionForUpdate(&x0);
-   if (adapt_lim_const > 0.)
+   else
    {
-      hr_solver.AddGridFunctionForUpdate(&adapt_lim_gf0);
-      hr_solver.AddFESpaceForUpdate(&ind_fes);
+      TMOPNewtonSolver solver(pfespace->GetComm(), ir, solver_type);
+      // Provide all integration rules in case of a mixed mesh.
+      solver.SetIntegrationRules(*irules, quad_order);
+      // Specify linear solver when we use a Newton-based solver.
+      if (solver_type == 0) { solver.SetPreconditioner(*S); }
+      // For untangling, the solver will update the min det(T) values.
+      solver.SetMinDetPtr(&min_detJ);
+      solver.SetMaxIter(solver_iter);
+      solver.SetRelTol(solver_rtol);
+      solver.SetAbsTol(0.0);
+      if (solver_art_type > 0)
+      {
+         solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
+      }
+      // Level of output.
+      IterativeSolver::PrintLevel newton_print;
+      if (verbosity_level > 0)
+      { newton_print.Errors().Warnings().Iterations(); }
+      solver.SetPrintLevel(newton_print);
+      // hr-adaptivity solver.
+      // If hr-adaptivity is disabled, r-adaptivity is done once using the
+      // TMOPNewtonSolver.
+      // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
+      // "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity.
+      // The solver terminates if an h-adaptivity iteration does not modify
+      // any element in the mesh.
+      TMOPHRSolver hr_solver(*pmesh, a, solver,
+                           x, move_bnd, hradaptivity,
+                           mesh_poly_deg, h_metric_id,
+                           n_hr_iter, n_h_iter);
+      hr_solver.AddGridFunctionForUpdate(&x0);
+      if (adapt_lim_const > 0.)
+      {
+         hr_solver.AddGridFunctionForUpdate(&adapt_lim_gf0);
+         hr_solver.AddFESpaceForUpdate(&ind_fes);
+      }
+      hr_solver.Mult();
    }
-   hr_solver.Mult();
 
    // 16. Save the optimized mesh to a file. This output can be viewed later
    //     using GLVis: "glvis -m optimized -np num_mpi_tasks".

miniapps/tools/tmop-check-metric.cpp

@@ -95,12 +95,23 @@ int main(int argc, char *argv[])
       // A-metrics
       case 11: metric = new TMOP_AMetric_011; break;
       case 36: metric = new TMOP_AMetric_036; break;
-      case 107: metric = new TMOP_AMetric_107a; break;
       case 126: metric = new TMOP_AMetric_126(0.9); break;
       default: cout << "Unknown metric_id: " << metric_id << endl; return 3;
    }
 
    const int dim = (metric_id < 300) ? 2 : 3;
+   Mesh *mesh;
+   if (dim == 2)
+   {
+      mesh = new Mesh(Mesh::MakeCartesian2D(1, 1, Element::QUADRILATERAL));
+   }
+   else
+   {
+      mesh = new Mesh(Mesh::MakeCartesian3D(1, 1, 1, Element::HEXAHEDRON));
+   }
+   H1_FECollection fec(2, dim);
+   FiniteElementSpace fespace(mesh, &fec, dim);
+
    DenseMatrix T(dim);
    Vector T_vec(T.GetData(), dim * dim);
 
@@ -113,6 +124,10 @@ int main(int argc, char *argv[])
       T(0, 0) += T_vec.Max();
       if (T.Det() <= 0.0) { continue; }
 
+      const DenseMatrix &Wideal = Geometries.GetGeomToPerfGeomJac(fespace.GetFE(
+                                                                     0)->GetGeomType());
+      metric->SetTargetJacobian(Wideal);
+
       const real_t i_form = metric->EvalW(T),
                    m_form = metric->EvalWMatrixForm(T);
       const real_t diff = std::abs(i_form - m_form) / std::abs(m_form);
@@ -131,17 +146,6 @@ int main(int argc, char *argv[])
    cout << "--- EvalW:     " << bad_cnt << " errors out of "
         << valid_cnt << " comparisons with det(T) > 0.\n";
 
-   Mesh *mesh;
-   if (dim == 2)
-   {
-      mesh = new Mesh(Mesh::MakeCartesian2D(1, 1, Element::QUADRILATERAL));
-   }
-   else
-   {
-      mesh = new Mesh(Mesh::MakeCartesian3D(1, 1, 1, Element::HEXAHEDRON));
-   }
-   H1_FECollection fec(2, dim);
-   FiniteElementSpace fespace(mesh, &fec, dim);
    NonlinearForm a(&fespace);
    mesh->SetNodalFESpace(&fespace);
    GridFunction x(&fespace);

miniapps/tools/tmop-metric-magnitude.cpp

@@ -89,7 +89,6 @@ int main(int argc, char *argv[])
       // A-metrics
       case 11: metric = new TMOP_AMetric_011; break;
       case 36: metric = new TMOP_AMetric_036; break;
-      case 107: metric = new TMOP_AMetric_107a; break;
       default: cout << "Unknown metric_id: " << metric_id << endl; return 3;
    }
 

tests/unit/CMakeLists.txt

@@ -40,6 +40,7 @@ set(UNIT_TESTS_SRCS
   linalg/test_matrix_rectangular.cpp
   linalg/test_matrix_sparse.cpp
   linalg/test_matrix_square.cpp
+  linalg/test_mma.cpp
   linalg/test_ode.cpp
   linalg/test_ode2.cpp
   linalg/test_operator.cpp

tests/unit/linalg/test_mma.cpp

@@ -0,0 +1,241 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "unit_tests.hpp"
+#include "mfem.hpp"
+
+using namespace mfem;
+
+real_t obj0(mfem::Vector& x)
+{
+   const int n=x.Size();
+   real_t rez=0.0;
+   for (int i=0; i<n; i++)
+   {
+      rez=rez+x[i]*x[i];
+   }
+
+#ifdef MFEM_USE_MPI
+   real_t grez;
+   MPI_Allreduce(&rez, &grez, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+   rez = grez;
+#endif
+
+   return rez;
+}
+
+real_t dobj0(mfem::Vector& x, mfem::Vector& dx)
+{
+   const int n=x.Size();
+   real_t rez=0.0;
+   for (int i=0; i<n; i++)
+   {
+      rez=rez+x[i]*x[i];
+      dx[i]=2.0*x[i];
+   }
+#ifdef MFEM_USE_MPI
+   real_t grez;
+   MPI_Allreduce(&rez, &grez, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+   rez = grez;
+#endif
+
+   return rez;
+}
+
+real_t g0(mfem::Vector& x)
+{
+   int n=x.Size();
+   real_t rez=0.0;
+   for (int i=0; i<n; i++)
+   {
+      rez=rez+x[i];
+   }
+
+   int gn = n;
+#ifdef MFEM_USE_MPI
+   real_t grez;
+   MPI_Allreduce(&n, &gn, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+   MPI_Allreduce(&rez, &grez, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+   rez = grez;
+#endif
+
+   rez=rez/gn;
+   return rez-2.0;
+}
+
+real_t dg0(mfem::Vector& x, mfem::Vector& dx)
+{
+   const int n=x.Size();
+
+   int gn = n;
+#ifdef MFEM_USE_MPI
+   MPI_Allreduce(&n, &gn, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+#endif
+
+   real_t rez=0.0;
+   for (int i=0; i<n; i++)
+   {
+      rez=rez+x[i];
+      dx[i]=1.0/gn;
+   }
+
+#ifdef MFEM_USE_MPI
+   real_t grez;
+   MPI_Allreduce(&rez, &grez, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+   rez = grez;
+#endif
+
+   rez=rez/gn;
+   return rez-1.0;
+}
+
+#ifdef MFEM_USE_MPI
+TEST_CASE("MMA Test 1", "[Parallel], [MMA_1]")
+{
+#else
+TEST_CASE("MMA Test 1", "[MMA_1]")
+{
+#endif
+   int world_size = 1;
+#ifdef MFEM_USE_MPI
+   MPI_Comm_size(MPI_COMM_WORLD, &world_size);
+#endif
+
+   int num_var=12 / world_size;
+
+   mfem::Vector x(num_var);
+   mfem::Vector dx(num_var);
+   mfem::Vector xmin(num_var); xmin=-1.0;
+   mfem::Vector xmax(num_var); xmax=2.0;
+   x=xmin; x+=0.5;
+
+   mfem::MMAOpt* mma = nullptr;
+
+#ifdef MFEM_USE_MPI
+   mma = new mfem::MMAOpt(MPI_COMM_WORLD,num_var,2,x);
+#else
+   mma = new mfem::MMAOpt(num_var,2,x);
+#endif
+
+   mfem::Vector g(2); g=-1.0;
+   mfem::Vector dg(num_var); dg=0.0;
+      mfem::Vector dg2(2*num_var); dg2=0.0;
+
+   real_t o;
+   for (int it=0; it<30; it++)
+   {
+      o=dobj0(x,dx);
+      g[0]=dg0(x,dg);
+
+         for (int i=0; i<num_var; i++)
+         {
+            dg2[i]=dg[i];
+         }
+
+      std::cout<<"it="<<it<<" o="<<o<<" g="<<g[0]<<std::endl;
+
+      for (int i=0; i<num_var; i++)
+      {
+         std::cout<<" "<<x[i];
+      }
+      std::cout<<std::endl;
+      for (int i=0; i<num_var; i++)
+      {
+         std::cout<<" "<<dx[i];
+      }
+      std::cout<<std::endl;
+
+      mma->Update(it,dx,g,dg2,xmin,xmax,x);
+      std::cout<<std::endl;
+   }
+
+   for (int i=0; i<num_var; i++)
+   {
+      std::cout<<" "<<x[i];
+   }
+   std::cout<<std::endl;
+
+   o=obj0(x);
+   std::cout<<"Final o="<<o<<std::endl;
+
+   delete mma;
+
+   REQUIRE( std::fabs(o - 0.0005790847638021212) < 1e-12 );
+}
+
+
+#ifdef MFEM_USE_MPI
+TEST_CASE("MMA Test", "[Parallel], [MMA]")
+{
+#else
+TEST_CASE("MMA Test", "[MMA]")
+{
+#endif
+   int world_size = 1;
+#ifdef MFEM_USE_MPI
+   MPI_Comm_size(MPI_COMM_WORLD, &world_size);
+#endif
+
+   int num_var=12 / world_size;
+
+   mfem::Vector x(num_var);
+   mfem::Vector dx(num_var);
+   mfem::Vector xmin(num_var); xmin=-1.0;
+   mfem::Vector xmax(num_var); xmax=2.0;
+   x=xmin; x+=0.5;
+
+   mfem::MMAOpt* mma = nullptr;
+
+#ifdef MFEM_USE_MPI
+   mma = new mfem::MMAOpt(MPI_COMM_WORLD,num_var,1,x);
+#else
+   mma = new mfem::MMAOpt(num_var,1,x);
+#endif
+
+   mfem::Vector g(1); g=-1.0;
+   mfem::Vector dg(num_var); dg=0.0;
+
+   real_t o;
+   for (int it=0; it<30; it++)
+   {
+      o=dobj0(x,dx);
+      g[0]=dg0(x,dg);
+
+      std::cout<<"it="<<it<<" o="<<o<<" g="<<g[0]<<std::endl;
+
+      for (int i=0; i<num_var; i++)
+      {
+         std::cout<<" "<<x[i];
+      }
+      std::cout<<std::endl;
+      for (int i=0; i<num_var; i++)
+      {
+         std::cout<<" "<<dx[i];
+      }
+      std::cout<<std::endl;
+
+      mma->Update(it,dx,g,dg,xmin,xmax,x);
+      std::cout<<std::endl;
+   }
+
+   for (int i=0; i<num_var; i++)
+   {
+      std::cout<<" "<<x[i];
+   }
+   std::cout<<std::endl;
+
+   o=obj0(x);
+   std::cout<<"Final o="<<o<<std::endl;
+
+   delete mma;
+
+   REQUIRE( std::fabs(o - 0.0005790847638021212) < 1e-12 );
+}