implement an interface to derefine by element number

2019-11-21 16:22:14 -08:00
55 changed files with 692 additions and 2741 deletions
@@ -18,7 +18,13 @@ Improved GPU support
  list of backends is: "occa-cuda", "raja-cuda", "cuda", "hip", "occa-omp",
  "raja-omp", "omp", "occa-cpu", "raja-cpu", and "cpu".

- Improved RAJA backend and multi-GPU MPI communications.
+- Improved RAJA backend.
+
+- Improved multi-GPU MPI communication.
+
+Miscellaneous
+-------------
+- Improved element numbering after uniform mesh refinement.

 Discretization improvements
 ---------------------------
@@ -26,29 +32,6 @@ Discretization improvements
  load balancing. Anisotropic prism refinement is only available in the serial
  version at the moment.

-Meshing improvements
--------------------
- The TMOP mesh optimization algorithms were extended to support r-adaptivity.
-  Target matrices can now be constructed either via a given analytical function
-  (e.g. spatial dependence of size, aspect ratio, etc., for each element) or via
-  a (Par)GridFunction specified on the original mesh.
-
- The TMOP mesh optimization algorithms have been improved to support AMR meshes.
-
- Added support for creating refined versions of periodic meshes, making use of
-  the new L2ElementRestriction class. This class also allows for computing
-  geometric factors on periodic meshes using partial assembly.
-
- Improved element numbering after uniform mesh refinement.
-
-New and updated examples and miniapps
-------------------------------------
- The mesh-optimizer and pmesh-optimizer miniapps have been updated to
-  demonstrate the new r-adaptivity capabilities of TMOP.
-
- The (p)mesh-optimizer miniapp has been updated to demonstrate mesh
-  optimization for an AMR mesh.
-
 Miscellaneous
 -------------
 - Upgraded the SUNDIALS interface to utilize SUNDIALS version 5.0.  This
@@ -56,6 +39,10 @@ Miscellaneous
  the application level.  Example usage of this new interface can be found
  in the examples/sundials directory.

+- Added support for creating refined versions of periodic meshes, making use of
+  the new L2ElementRestriction class. This class also allows for computing
+  geometric factors on periodic meshes using partial assembly.
+

 Version 4.0, released on May 24, 2019
 =====================================
@@ -254,7 +254,7 @@ endif()

 # Axom/Sidre
 if (MFEM_USE_SIDRE)
-  find_package(Axom REQUIRED Axom)
+  find_package(Axom REQUIRED Sidre SLIC axom_utils)
 endif()

 # PUMI
@@ -383,11 +383,11 @@ MFEM_USE_MPFR = YES/NO
   see below.

 MFEM_USE_SIDRE = YES/NO
-   Sidre is a component of LLNL's axom project, https://github.com/LLNL/axom,
-   that provides an HDF5-based file format for visualization or restart
-   capability following the Conduit (https://github.com/LLNL/conduit) mesh
-   blueprint specification. When enabled, this option requires installation of
-   HDF5 (see also MFEM_USE_NETCDF), Conduit and LLNL's axom project.
+   Sidre is a component of LLNL's axom project, http://goo.gl/cZyJdn, that
+   provides an HDF5-based file format for visualization or restart capability
+   following the Conduit (https://github.com/LLNL/conduit) mesh blueprint
+   specification. When enabled, this option requires installation of HDF5 (see
+   also MFEM_USE_NETCDF), Conduit and LLNL's axom project.

 MFEM_USE_CONDUIT = YES/NO
   Enables support for converting MFEM Mesh and Grid Function objects to and
@@ -543,8 +543,7 @@ The specific libraries and their options are:
  Options: PETSC_OPT, PETSC_LIB.

 - Sidre (optional), part of LLNL's axom project, used when MFEM_USE_SIDRE = YES.
-  Starting with MFEM v4.1, Axom version 0.3.1 or later is required.
-  URL: https://github.com/LLNL/axom
+  URL: http://goo.gl/cZyJdn (axom, to be released)
       https://github.com/LLNL/conduit (Conduit)
       https://support.hdfgroup.org/HDF5 (HDF5)
  Options: SIDRE_OPT, SIDRE_LIB.
@@ -715,7 +714,6 @@ MFEM_USE_PUMI
 MFEM_USE_CUDA
 MFEM_USE_OCCA
 MFEM_USE_RAJA
-MFEM_USE_SIDRE

 The following options are CMake specific:

@@ -764,7 +762,6 @@ The CMake build system adds auto-detection for the following packages/libraries:
 - PUMI
 - OCCA
 - RAJA
- - AXOM - Used when MFEM_USE_SIDRE is enabled

 The following built-in CMake packages are also used:

@@ -18,4 +18,6 @@ include(MfemCmakeUtilities)
 # Note: components are enabled based on the find_package() parameters.
 mfem_find_package(Axom AXOM AXOM_DIR "include" "" "lib" ""
  "Paths to headers required by Axom." "Libraries required by Axom."
-  ADD_COMPONENT Axom "include" axom/config.hpp "lib" axom)
+  ADD_COMPONENT Sidre "include" sidre/sidre.hpp "lib" sidre
+  ADD_COMPONENT SLIC "include" slic/slic.hpp "lib" slic
+  ADD_COMPONENT axom_utils "include" axom_utils/Utilities.hpp "lib" axom_utils)
@@ -154,7 +154,7 @@ set(CONDUIT_DIR "${MFEM_DIR}/../conduit" CACHE PATH

 set(AXOM_DIR "${MFEM_DIR}/../axom" CACHE PATH "Path to the Axom library.")
 # May need to add "Boost" as requirement.
-set(Axom_REQUIRED_PACKAGES "Conduit/relay/blueprint" CACHE STRING
+set(Axom_REQUIRED_PACKAGES "Conduit/relay" CACHE STRING
    "Additional packages required by Axom.")

 set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
@@ -299,7 +299,7 @@ SIDRE_LIB = \
   -Wl,-rpath,$(SIDRE_DIR)/lib -L$(SIDRE_DIR)/lib \
   -Wl,-rpath,$(CONDUIT_DIR)/lib -L$(CONDUIT_DIR)/lib \
   -Wl,-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib \
-   -laxom -lconduit -lconduit_relay -lconduit_blueprint -lhdf5 $(ZLIB_LIB) -ldl
+   -lsidre -lslic -laxom_utils -lconduit -lconduit_relay -lhdf5 $(ZLIB_LIB) -ldl

 # PUMI
 # Note that PUMI_DIR is needed -- it is used to check for gmi_sim.h
@@ -32,7 +32,6 @@ set(SRCS
  nonlininteg.cpp
  staticcond.cpp
  tmop.cpp
-  tmop_tools.cpp
  )

 set(HDRS
@@ -65,7 +64,6 @@ set(HDRS
  tfespace.hpp
  tintrules.hpp
  tmop.hpp
-  tmop_tools.hpp
  )

 if (MFEM_USE_SIDRE)
@@ -27,7 +27,7 @@ static void OccaPADiffusionSetup2D(const int D1D,
                                   const int NE,
                                   const Array<double> &W,
                                   const Vector &J,
-                                   const Vector &C,
+                                   const double COEFF,
                                   Vector &op)
 {
   occa::properties props;
@@ -35,9 +35,7 @@ static void OccaPADiffusionSetup2D(const int D1D,
   props["defines/Q1D"] = Q1D;
   const occa::memory o_W = OccaMemoryRead(W.GetMemory(), W.Size());
   const occa::memory o_J = OccaMemoryRead(J.GetMemory(), J.Size());
-   const occa::memory o_C = OccaMemoryRead(C.GetMemory(), C.Size());
   occa::memory o_op = OccaMemoryWrite(op.GetMemory(), op.Size());
-   const bool const_c = C.Size() == 1;
   const occa_id_t id = std::make_pair(D1D,Q1D);
   static occa_kernel_t OccaDiffSetup2D_ker;
   if (OccaDiffSetup2D_ker.find(id) == OccaDiffSetup2D_ker.end())
@@ -47,7 +45,7 @@ static void OccaPADiffusionSetup2D(const int D1D,
                                     "DiffusionSetup2D", props);
      OccaDiffSetup2D_ker.emplace(id, DiffusionSetup2D);
   }
-   OccaDiffSetup2D_ker.at(id)(NE, o_W, o_J, o_C, o_op, const_c);
+   OccaDiffSetup2D_ker.at(id)(NE, o_W, o_J, COEFF, o_op);
 }

 static void OccaPADiffusionSetup3D(const int D1D,
@@ -55,7 +53,7 @@ static void OccaPADiffusionSetup3D(const int D1D,
                                   const int NE,
                                   const Array<double> &W,
                                   const Vector &J,
-                                   const Vector &C,
+                                   const double COEFF,
                                   Vector &op)
 {
   occa::properties props;
@@ -63,9 +61,7 @@ static void OccaPADiffusionSetup3D(const int D1D,
   props["defines/Q1D"] = Q1D;
   const occa::memory o_W = OccaMemoryRead(W.GetMemory(), W.Size());
   const occa::memory o_J = OccaMemoryRead(J.GetMemory(), J.Size());
-   const occa::memory o_C = OccaMemoryRead(C.GetMemory(), C.Size());
   occa::memory o_op = OccaMemoryWrite(op.GetMemory(), op.Size());
-   const bool const_c = C.Size() == 1;
   const occa_id_t id = std::make_pair(D1D,Q1D);
   static occa_kernel_t OccaDiffSetup3D_ker;
   if (OccaDiffSetup3D_ker.find(id) == OccaDiffSetup3D_ker.end())
@@ -75,7 +71,7 @@ static void OccaPADiffusionSetup3D(const int D1D,
                                     "DiffusionSetup3D", props);
      OccaDiffSetup3D_ker.emplace(id, DiffusionSetup3D);
   }
-   OccaDiffSetup3D_ker.at(id)(NE, o_W, o_J, o_C, o_op, const_c);
+   OccaDiffSetup3D_ker.at(id)(NE, o_W, o_J, COEFF, o_op);
 }
 #endif // MFEM_USE_OCCA

@@ -84,16 +80,14 @@ static void PADiffusionSetup2D(const int Q1D,
                               const int NE,
                               const Array<double> &w,
                               const Vector &j,
-                               const Vector &c,
-                               Vector &d)
+                               const double COEFF,
+                               Vector &op)
 {
   const int NQ = Q1D*Q1D;
-   const bool const_c = c.Size() == 1;
   auto W = w.Read();

   auto J = Reshape(j.Read(), NQ, 2, 2, NE);
-   auto C = const_c ? Reshape(c.Read(), 1, 1) : Reshape(c.Read(), NQ, NE);
-   auto D = Reshape(d.Write(), NQ, 3, NE);
+   auto y = Reshape(op.Write(), NQ, 3, NE);

   MFEM_FORALL(e, NE,
   {
@@ -103,11 +97,10 @@ static void PADiffusionSetup2D(const int Q1D,
         const double J21 = J(q,1,0,e);
         const double J12 = J(q,0,1,e);
         const double J22 = J(q,1,1,e);
-         const double coeff = const_c ? C(0,0) : C(q,e);
-         const double c_detJ = W[q] * coeff / ((J11*J22)-(J21*J12));
-         D(q,0,e) =  c_detJ * (J12*J12 + J22*J22); // 1,1
-         D(q,1,e) = -c_detJ * (J12*J11 + J22*J21); // 1,2
-         D(q,2,e) =  c_detJ * (J11*J11 + J21*J21); // 2,2
+         const double c_detJ = W[q] * COEFF / ((J11*J22)-(J21*J12));
+         y(q,0,e) =  c_detJ * (J12*J12 + J22*J22); // 1,1
+         y(q,1,e) = -c_detJ * (J12*J11 + J22*J21); // 1,2
+         y(q,2,e) =  c_detJ * (J11*J11 + J21*J21); // 2,2
      }
   });
 }
@@ -117,15 +110,13 @@ static void PADiffusionSetup3D(const int Q1D,
                               const int NE,
                               const Array<double> &w,
                               const Vector &j,
-                               const Vector &c,
-                               Vector &d)
+                               const double COEFF,
+                               Vector &op)
 {
   const int NQ = Q1D*Q1D*Q1D;
-   const bool const_c = c.Size() == 1;
   auto W = w.Read();
   auto J = Reshape(j.Read(), NQ, 3, 3, NE);
-   auto C = const_c ? Reshape(c.Read(), 1, 1) : Reshape(c.Read(), NQ, NE);
-   auto D = Reshape(d.Write(), NQ, 6, NE);
+   auto y = Reshape(op.Write(), NQ, 6, NE);
   MFEM_FORALL(e, NE,
   {
      for (int q = 0; q < NQ; ++q)
@@ -142,8 +133,7 @@ static void PADiffusionSetup3D(const int Q1D,
         const double detJ = J11 * (J22 * J33 - J32 * J23) -
         /* */               J21 * (J12 * J33 - J32 * J13) +
         /* */               J31 * (J12 * J23 - J22 * J13);
-         const double coeff = const_c ? C(0,0) : C(q,e);
-         const double c_detJ = W[q] * coeff / detJ;
+         const double c_detJ = W[q] * COEFF / detJ;
         // adj(J)
         const double A11 = (J22 * J33) - (J23 * J32);
         const double A12 = (J32 * J13) - (J12 * J33);
@@ -155,12 +145,12 @@ static void PADiffusionSetup3D(const int Q1D,
         const double A32 = (J31 * J12) - (J11 * J32);
         const double A33 = (J11 * J22) - (J12 * J21);
         // detJ J^{-1} J^{-T} = (1/detJ) adj(J) adj(J)^T
-         D(q,0,e) = c_detJ * (A11*A11 + A12*A12 + A13*A13); // 1,1
-         D(q,1,e) = c_detJ * (A11*A21 + A12*A22 + A13*A23); // 2,1
-         D(q,2,e) = c_detJ * (A11*A31 + A12*A32 + A13*A33); // 3,1
-         D(q,3,e) = c_detJ * (A21*A21 + A22*A22 + A23*A23); // 2,2
-         D(q,4,e) = c_detJ * (A21*A31 + A22*A32 + A23*A33); // 3,2
-         D(q,5,e) = c_detJ * (A31*A31 + A32*A32 + A33*A33); // 3,3
+         y(q,0,e) = c_detJ * (A11*A11 + A12*A12 + A13*A13); // 1,1
+         y(q,1,e) = c_detJ * (A11*A21 + A12*A22 + A13*A23); // 2,1
+         y(q,2,e) = c_detJ * (A11*A31 + A12*A32 + A13*A33); // 3,1
+         y(q,3,e) = c_detJ * (A21*A21 + A22*A22 + A23*A23); // 2,2
+         y(q,4,e) = c_detJ * (A21*A31 + A22*A32 + A23*A33); // 3,2
+         y(q,5,e) = c_detJ * (A31*A31 + A32*A32 + A33*A33); // 3,3
      }
   });
 }
@@ -171,8 +161,8 @@ static void PADiffusionSetup(const int dim,
                             const int NE,
                             const Array<double> &W,
                             const Vector &J,
-                             const Vector &C,
-                             Vector &D)
+                             const double COEFF,
+                             Vector &op)
 {
   if (dim == 1) { MFEM_ABORT("dim==1 not supported in PADiffusionSetup"); }
   if (dim == 2)
@@ -180,22 +170,22 @@ static void PADiffusionSetup(const int dim,
 #ifdef MFEM_USE_OCCA
      if (DeviceCanUseOcca())
      {
-         OccaPADiffusionSetup2D(D1D, Q1D, NE, W, J, C, D);
+         OccaPADiffusionSetup2D(D1D, Q1D, NE, W, J, COEFF, op);
         return;
      }
 #endif // MFEM_USE_OCCA
-      PADiffusionSetup2D(Q1D, NE, W, J, C, D);
+      PADiffusionSetup2D(Q1D, NE, W, J, COEFF, op);
   }
   if (dim == 3)
   {
 #ifdef MFEM_USE_OCCA
      if (DeviceCanUseOcca())
      {
-         OccaPADiffusionSetup3D(D1D, Q1D, NE, W, J, C, D);
+         OccaPADiffusionSetup3D(D1D, Q1D, NE, W, J, COEFF, op);
         return;
      }
 #endif // MFEM_USE_OCCA
-      PADiffusionSetup3D(Q1D, NE, W, J, C, D);
+      PADiffusionSetup3D(Q1D, NE, W, J, COEFF, op);
   }
 }

@@ -215,32 +205,11 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
   dofs1D = maps->ndof;
   quad1D = maps->nqpt;
   pa_data.SetSize(symmDims * nq * ne, Device::GetMemoryType());
-   Vector coeff;
-   if (Q == nullptr)
-   {
-      coeff.SetSize(1);
-      coeff(0) = 1.0;
-   }
-   else if (ConstantCoefficient* cQ = dynamic_cast<ConstantCoefficient*>(Q))
-   {
-      coeff.SetSize(1);
-      coeff(0) = cQ->constant;
-   }
-   else
-   {
-      coeff.SetSize(nq * ne);
-      auto C = Reshape(coeff.Write(), nq, ne);
-      for (int e = 0; e < ne; ++e)
-      {
-         ElementTransformation& T = *fes.GetElementTransformation(e);
-         for (int q = 0; q < nq; ++q)
-         {
-            C(q,e) = Q->Eval(T, ir->IntPoint(q));
-         }
-      }
-   }
-   PADiffusionSetup(dim, dofs1D, quad1D, ne, ir->GetWeights(), geom->J, coeff,
-                    pa_data);
+   ConstantCoefficient *cQ = dynamic_cast<ConstantCoefficient*>(Q);
+   MFEM_VERIFY(cQ != NULL, "only ConstantCoefficient is supported!");
+   const double coeff = cQ->constant;
+   PADiffusionSetup(dim, dofs1D, quad1D, ne, ir->GetWeights(), geom->J,
+                    coeff, pa_data);
 }

 #ifdef MFEM_USE_OCCA
@@ -38,40 +38,24 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
   dofs1D = maps->ndof;
   quad1D = maps->nqpt;
   pa_data.SetSize(ne*nq, Device::GetMemoryType());
-   Vector coeff;
-   if (Q == nullptr)
-   {
-      coeff.SetSize(1);
-      coeff(0) = 1.0;
-   }
-   else if (ConstantCoefficient* cQ = dynamic_cast<ConstantCoefficient*>(Q))
-   {
-      coeff.SetSize(1);
-      coeff(0) = cQ->constant;
-   }
-   else
-   {
-      coeff.SetSize(nq * ne);
-      auto C = Reshape(coeff.Write(), nq, ne);
-      for (int e = 0; e < ne; ++e)
-      {
-         ElementTransformation& T = *fes.GetElementTransformation(e);
-         for (int q = 0; q < nq; ++q)
-         {
-            C(q,e) = Q->Eval(T, ir->IntPoint(q));
-         }
-      }
-   }
+   ConstantCoefficient *const_coeff = dynamic_cast<ConstantCoefficient*>(Q);
+   // TODO: other types of coefficients ...
   if (dim==1) { MFEM_ABORT("Not supported yet... stay tuned!"); }
   if (dim==2)
   {
+      double constant = 0.0;
+      if (const_coeff)
+      {
+         constant = const_coeff->constant;
+      }
+      else
+      {
+         MFEM_ABORT("Coefficient type not supported");
+      }
      const int NE = ne;
      const int NQ = nq;
-      const bool const_c = coeff.Size() == 1;
      auto w = ir->GetWeights().Read();
      auto J = Reshape(geom->J.Read(), NQ,2,2,NE);
-      auto C =
-         const_c ? Reshape(coeff.Read(), 1,1) : Reshape(coeff.Read(), NQ,NE);
      auto v = Reshape(pa_data.Write(), NQ, NE);
      MFEM_FORALL(e, NE,
      {
@@ -82,20 +66,25 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
            const double J21 = J(q,0,1,e);
            const double J22 = J(q,1,1,e);
            const double detJ = (J11*J22)-(J21*J12);
-            const double coeff = const_c ? C(0,0) : C(q,e);
-            v(q,e) =  w[q] * coeff * detJ;
+            v(q,e) =  w[q] * constant * detJ;
         }
      });
   }
   if (dim==3)
   {
+      double constant = 0.0;
+      if (const_coeff)
+      {
+         constant = const_coeff->constant;
+      }
+      else
+      {
+         MFEM_ABORT("Coefficient type not supported");
+      }
      const int NE = ne;
      const int NQ = nq;
-      const bool const_c = coeff.Size() == 1;
      auto W = ir->GetWeights().Read();
      auto J = Reshape(geom->J.Read(), NQ,3,3,NE);
-      auto C =
-         const_c ? Reshape(coeff.Read(), 1,1) : Reshape(coeff.Read(), NQ,NE);
      auto v = Reshape(pa_data.Write(), NQ,NE);
      MFEM_FORALL(e, NE,
      {
@@ -107,8 +96,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
            const double detJ = J11 * (J22 * J33 - J32 * J23) -
            /* */               J21 * (J12 * J33 - J32 * J13) +
            /* */               J31 * (J12 * J23 - J22 * J13);
-            const double coeff = const_c ? C(0,0) : C(q,e);
-            v(q,e) = W[q] * coeff * detJ;
+            v(q,e) = W[q] * constant * detJ;
         }
      });
   }
@@ -9635,7 +9635,6 @@ void L2_TetrahedronElement::ProjectDelta(int vertex, Vector &dofs) const
            const IntegrationPoint &ip = Nodes.IntPoint(i);
            dofs[i] = pow(ip.y, Order);
         }
-         break;
      case 3:
         for (int i = 0; i < Dof; i++)
         {
@@ -31,7 +31,6 @@
 #include "estimators.hpp"
 #include "staticcond.hpp"
 #include "tmop.hpp"
-#include "tmop_tools.hpp"

 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
@@ -1716,7 +1716,6 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
   Array<int> vdofs;
   Vector vals;

-   HostWrite();
   // maximal element attribute for each dof
   dof_attr.SetSize(fes->GetVSize());
   dof_attr = -1;
@@ -434,8 +434,6 @@ public:
   /** The GridFunction is resized using the SetSize() method. */
   virtual void SetSpace(FiniteElementSpace *f);

-   using Vector::MakeRef;
-
   /** @brief Make the GridFunction reference external data on a new
       FiniteElementSpace. */
   /** This method changes the FiniteElementSpace associated with the
@@ -350,6 +350,7 @@ void VectorFEDomainLFIntegrator::AssembleDeltaElementVect(
   vshape.Mult(vec, elvect);
 }

+
 void VectorBoundaryFluxLFIntegrator::AssembleRHSElementVect(
   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {
@@ -396,26 +397,19 @@ void VectorFEBoundaryFluxLFIntegrator::AssembleRHSElementVect(
   if (ir == NULL)
   {
      int intorder = 2*el.GetOrder();  // <----------
-      if (F == NULL)
-      {
-         intorder -= el.GetOrder() + 1;
-      }
      ir = &IntRules.Get(el.GetGeomType(), intorder);
   }

   for (int i = 0; i < ir->GetNPoints(); i++)
   {
      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      Tr.SetIntPoint (&ip);
+      double val = ip.weight*F.Eval(Tr, ip);
+
      el.CalcShape(ip, shape);

-      double val = ip.weight;
-      if (F)
-      {
-         Tr.SetIntPoint (&ip);
-         val *= F->Eval(Tr, ip);
-      }
-
-      elvect.Add(val, shape);
+      add(elvect, val, shape, elvect);
   }
 }

@@ -279,12 +279,11 @@ public:
 class VectorFEBoundaryFluxLFIntegrator : public LinearFormIntegrator
 {
 private:
-   Coefficient *F;
+   Coefficient &F;
   Vector shape;

 public:
-   VectorFEBoundaryFluxLFIntegrator() : F(NULL) { }
-   VectorFEBoundaryFluxLFIntegrator(Coefficient &f) : F(&f) { }
+   VectorFEBoundaryFluxLFIntegrator(Coefficient &f) : F(f) { }

   virtual void AssembleRHSElementVect(const FiniteElement &el,
                                       ElementTransformation &Tr,
@@ -65,8 +65,6 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
   Vector el_x;
   const FiniteElement *fe;
   ElementTransformation *T;
-   Mesh *mesh = fes->GetMesh();
-
   double energy = 0.0;

   if (dnfi.Size())
@@ -86,81 +84,14 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const

   if (fnfi.Size())
   {
-      FaceElementTransformations *tr;
-      const FiniteElement *fe1, *fe2;
-      Array<int> vdofs2;
-
-      for (int i = 0; i < mesh->GetNumFaces(); i++)
-      {
-         tr = mesh->GetInteriorFaceTransformations(i);
-         if (tr != NULL)
-         {
-            fes->GetElementVDofs(tr->Elem1No, vdofs);
-            fes->GetElementVDofs(tr->Elem2No, vdofs2);
-            vdofs.Append (vdofs2);
-            x.GetSubVector(vdofs, el_x);
-            fe1 = fes->GetFE(tr->Elem1No);
-            fe2 = fes->GetFE(tr->Elem2No);
-            for (int k = 0; k < fnfi.Size(); k++)
-            {
-               energy += fnfi[k]->GetFaceEnergy(*fe1, *fe2, *tr, el_x);
-            }
-         }
-      }
+      MFEM_ABORT("TODO: add energy contribution from interior face terms");
   }

   if (bfnfi.Size())
   {
-      FaceElementTransformations *tr;
-      const FiniteElement *fe1, *fe2;
-
-      // Which boundary attributes need to be processed?
-      Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
-                                 mesh->bdr_attributes.Max() : 0);
-      bdr_attr_marker = 0;
-      for (int k = 0; k < bfnfi.Size(); k++)
-      {
-         if (bfnfi_marker[k] == NULL)
-         {
-            bdr_attr_marker = 1;
-            break;
-         }
-         Array<int> &bdr_marker = *bfnfi_marker[k];
-         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
-                     "invalid boundary marker for boundary face integrator #"
-                     << k << ", counting from zero");
-         for (int i = 0; i < bdr_attr_marker.Size(); i++)
-         {
-            bdr_attr_marker[i] |= bdr_marker[i];
-         }
-      }
-
-      for (int i = 0; i < fes -> GetNBE(); i++)
-      {
-         const int bdr_attr = mesh->GetBdrAttribute(i);
-         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
-
-         tr = mesh->GetBdrFaceTransformations (i);
-         if (tr != NULL)
-         {
-            fes->GetElementVDofs(tr->Elem1No, vdofs);
-            x.GetSubVector(vdofs, el_x);
-
-            fe1 = fes->GetFE(tr->Elem1No);
-            // The fe2 object is really a dummy and not used on the boundaries,
-            // but we can't dereference a NULL pointer, and we don't want to
-            // actually make a fake element.
-            fe2 = fe1;
-            for (int k = 0; k < bfnfi.Size(); k++)
-            {
-               if (bfnfi_marker[k] &&
-                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
-
-               energy += bfnfi[k]->GetFaceEnergy(*fe1, *fe2, *tr, el_x);
-            }
-         }
-      }
+      MFEM_ABORT("TODO: add energy contribution from boundary face terms");
   }
+
   return energy;
 }

@@ -111,7 +111,7 @@ public:
       be fes->GetVSize(). */
   double GetGridFunctionEnergy(const Vector &x) const;

-   /// Compute the energy corresponding to the state @a x.
+   /// Compute the enery corresponding to the state @a x.
   /** In general, @a x may have non-homogeneous essential boundary values.

       The state @a x must be a true-dof vector. */
@@ -55,14 +55,6 @@ double NonlinearFormIntegrator::GetElementEnergy(
   return 0.0;
 }

-double NonlinearFormIntegrator::GetFaceEnergy(
-   const FiniteElement &el1, const FiniteElement &el2,
-   FaceElementTransformations &Tr, const Vector &elfun)
-{
-   mfem_error("NonlinearFormIntegrator::GetFaceEnergy"
-              " is not overloaded!");
-   return 0.0;
-}

 void BlockNonlinearFormIntegrator::AssembleElementVector(
   const Array<const FiniteElement *> &el,
@@ -63,17 +63,11 @@ public:
                                 FaceElementTransformations &Tr,
                                 const Vector &elfun, DenseMatrix &elmat);

-   /// Compute the local energy/functional
+   /// Compute the local energy
   virtual double GetElementEnergy(const FiniteElement &el,
                                   ElementTransformation &Tr,
                                   const Vector &elfun);

-   /// Compute the face(s) contribution to the energy/functional
-   virtual double GetFaceEnergy(const FiniteElement &el1,
-                                const FiniteElement &el2,
-                                FaceElementTransformations &Tr,
-                                const Vector &elfun);
-
   virtual ~NonlinearFormIntegrator() { }
 };

@@ -38,24 +38,19 @@ typedef double* QLocal3D_t @dim(Q1D, Q1D, Q1D, NE);
 typedef double* Jacobian2D_t @dim(Q2D, 2, 2, NE);
 typedef double* Jacobian3D_t @dim(Q3D, 3, 3, NE);

-typedef double* Coeff2D_t @dim(Q2D, NE);
-typedef double* Coeff3D_t @dim(Q3D, NE);
-
 typedef double* SymmOperator2D_t @dim(Q2D, 3, NE);
 typedef double* SymmOperator3D_t @dim(Q3D, 6, NE);

@kernel void DiffusionSetup2D(const int NE,
                              @restrict const double *W,
                              @restrict const Jacobian2D_t J,
-                              @restrict const Coeff2D_t C,
-                              @restrict SymmOperator2D_t op,
-                              const bool const_c) {
+                              const double COEFF,
+                              @restrict SymmOperator2D_t op) {
  for (int e = 0; e < NE; ++e; @outer) {
    for (int q = 0; q < Q2D; ++q; @inner) {
      const double J11 = J(q, 0, 0, e), J12 = J(q, 1, 0, e);
      const double J21 = J(q, 0, 1, e), J22 = J(q, 1, 1, e);
-      const double coeff = const_c ? C(0,0) : C(q,e);
-      const double c_detJ = W[q] * coeff / ((J11 * J22) - (J21 * J12));
+      const double c_detJ = W[q] * COEFF / ((J11 * J22) - (J21 * J12));
      op(q, 0, e) =  c_detJ * (J21*J21 + J22*J22); // (1,1)
      op(q, 1, e) = -c_detJ * (J21*J11 + J22*J12); // (1,2), (2,1)
      op(q, 2, e) =  c_detJ * (J11*J11 + J12*J12); // (2,2)
@@ -66,9 +61,8 @@ typedef double* SymmOperator3D_t @dim(Q3D, 6, NE);
@kernel void DiffusionSetup3D(const int NE,
                              @restrict const double *W,
                              @restrict const Jacobian3D_t J,
-                              @restrict const Coeff3D_t C,
-                              @restrict SymmOperator3D_t op,
-                              const bool const_c) {
+                              const double COEFF,
+                              @restrict SymmOperator3D_t op) {
  for (int e = 0; e < NE; ++e; @outer) {
    for (int q = 0; q < Q3D; ++q; @inner) {
      const double J11 = J(q, 0, 0, e), J12 = J(q, 1, 0, e), J13 = J(q, 2, 0, e);
@@ -78,8 +72,7 @@ typedef double* SymmOperator3D_t @dim(Q3D, 6, NE);
      const double detJ = ((J11 * J22 * J33) + (J12 * J23 * J31) + (J13 * J21 * J32) -
                           (J13 * J22 * J31) - (J12 * J21 * J33) - (J11 * J23 * J32));

-      const double coeff = const_c ? C(0,0) : C(q,e);
-      const double c_detJ = W[q] * coeff / detJ;
+      const double c_detJ = W[q] * COEFF / detJ;

      // adj(J)
      const double A11 = (J22 * J33) - (J23 * J32);
@@ -225,13 +225,11 @@ void ParGridFunction::ExchangeFaceNbrData()
   MPI_Request *recv_requests = requests + num_face_nbrs;
   MPI_Status  *statuses = new MPI_Status[num_face_nbrs];

-   const double *h_data = this->HostRead();
   for (int i = 0; i < send_data.Size(); i++)
   {
-      send_data[i] = h_data[send_ldof[i]];
+      send_data[i] = data[send_ldof[i]];
   }

-   double *h_face_nbr_data = face_nbr_data.HostWrite();
   for (int fn = 0; fn < num_face_nbrs; fn++)
   {
      int nbr_rank = pmesh->GetFaceNbrRank(fn);
@@ -241,7 +239,7 @@ void ParGridFunction::ExchangeFaceNbrData()
                send_offset[fn+1] - send_offset[fn],
                MPI_DOUBLE, nbr_rank, tag, MyComm, &send_requests[fn]);

-      MPI_Irecv(&h_face_nbr_data[recv_offset[fn]],
+      MPI_Irecv(&face_nbr_data(recv_offset[fn]),
                recv_offset[fn+1] - recv_offset[fn],
                MPI_DOUBLE, nbr_rank, tag, MyComm, &recv_requests[fn]);
   }
@@ -112,8 +112,6 @@ public:
   /// Associate a new parallel space with the ParGridFunction.
   void SetSpace(ParFiniteElementSpace *f);

-   using GridFunction::MakeRef;
-
   /** @brief Make the ParGridFunction reference external data on a new
       FiniteElementSpace. */
   /** This method changes the FiniteElementSpace associated with the
@@ -46,7 +46,7 @@ double ParNonlinearForm::GetParGridFunctionEnergy(const Vector &x) const
 void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
 {
   NonlinearForm::Mult(x, y); // x --(P)--> aux1 --(A_local)--> aux2
-   Y.MakeRef(aux2, 0); // aux2 contains A_local.P.x
+   Y.SetData(aux2.GetData()); // aux2 contains A_local.P.x

   if (fnfi.Size())
   {
@@ -58,7 +58,7 @@ void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
      Array<int> vdofs1, vdofs2;
      Vector el_x, el_y;

-      X.MakeRef(aux1, 0); // aux1 contains P.x
+      X.SetData(aux1.GetData()); // aux1 contains P.x
      X.ExchangeFaceNbrData();
      const int n_shared_faces = pmesh->GetNSharedFaces();
      for (int i = 0; i < n_shared_faces; i++)
@@ -16,7 +16,9 @@

 #include "fem.hpp"

-#include <axom/sidre.hpp>
+#ifdef MFEM_USE_MPI
+#include <sidre/IOManager.hpp>
+#endif

 #include <string>
 #include <iomanip>      // for setw, setfill
@@ -202,10 +204,10 @@ SidreDataCollection::get_file_path(const std::string &filename) const

 axom::sidre::View *
 SidreDataCollection::AllocNamedBuffer(const std::string& buffer_name,
-                                      axom::sidre::IndexType sz,
+                                      axom::sidre::SidreLength sz,
                                      axom::sidre::TypeID type)
 {
-   sz = std::max(sz, sidre::IndexType(0));
+   sz = std::max(sz, sidre::SidreLength(0));
   sidre::Group *f = named_buffers_grp();
   sidre::View  *v = NULL;

@@ -823,7 +825,7 @@ void SidreDataCollection::Save(const std::string& filename,
 void SidreDataCollection::
 addScalarBasedGridFunction(const std::string &field_name, GridFunction *gf,
                           const std::string &buffer_name,
-                           axom::sidre::IndexType offset)
+                           axom::sidre::SidreLength offset)
 {
   sidre::Group* grp = m_bp_grp->getGroup("fields/" + field_name);
   MFEM_ASSERT(grp != NULL, "field " << field_name << " does not exist");
@@ -886,7 +888,7 @@ addScalarBasedGridFunction(const std::string &field_name, GridFunction *gf,
 void SidreDataCollection::
 addVectorBasedGridFunction(const std::string& field_name, GridFunction *gf,
                           const std::string &buffer_name,
-                           axom::sidre::IndexType offset)
+                           axom::sidre::SidreLength offset)
 {
   sidre::Group* grp = m_bp_grp->getGroup("fields/" + field_name);
   MFEM_ASSERT(grp != NULL, "field " << field_name << " does not exist");
@@ -1011,7 +1013,7 @@ DeregisterFieldInBPIndex(const std::string& field_name)
 void SidreDataCollection::RegisterField(const std::string &field_name,
                                        GridFunction *gf,
                                        const std::string &buffer_name,
-                                        axom::sidre::IndexType offset)
+                                        axom::sidre::SidreLength offset)
 {
   if ( field_name.empty() || buffer_name.empty() ||
        gf == NULL || gf->FESpace() == NULL )
@@ -25,7 +25,7 @@
 #  pragma GCC diagnostic ignored "-Wpedantic"
 # endif
 #endif
-#include <axom/sidre.hpp>
+#include <sidre/sidre.hpp>
 #ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
 # pragma GCC diagnostic pop
 #endif
@@ -246,7 +246,7 @@ public:
    */
   void RegisterField(const std::string &field_name, GridFunction *gf,
                      const std::string &buffer_name,
-                      axom::sidre::IndexType offset);
+                      axom::sidre::SidreLength offset);

   /// Registers an attribute field in the Sidre DataStore
   /** The registration process is similar to that of RegisterField()
@@ -385,7 +385,7 @@ public:
    */
   axom::sidre::View *
   AllocNamedBuffer(const std::string& buffer_name,
-                    axom::sidre::IndexType sz,
+                    axom::sidre::SidreLength sz,
                    axom::sidre::TypeID type =
                       axom::sidre::DOUBLE_ID);

@@ -469,7 +469,7 @@ private:
   void addScalarBasedGridFunction(const std::string& field_name,
                                   GridFunction* gf,
                                   const std::string &buffer_name,
-                                   axom::sidre::IndexType offset);
+                                   axom::sidre::SidreLength offset);

   /**
    * \brief A private helper function to set up the views associated with the
@@ -483,7 +483,7 @@ private:
   void addVectorBasedGridFunction(const std::string& field_name,
                                   GridFunction* gf,
                                   const std::string &buffer_name,
-                                   axom::sidre::IndexType offset);
+                                   axom::sidre::SidreLength offset);

   /** @brief A private helper function to set up the Views associated with
       attribute field named @a field_name */
@@ -12,7 +12,6 @@
 #include "tmop.hpp"
 #include "linearform.hpp"
 #include "pgridfunc.hpp"
-#include "tmop_tools.hpp"

 namespace mfem
 {
@@ -769,7 +768,7 @@ void TMOP_Metric_352::AssembleH(const DenseMatrix &Jpt,
 void TargetConstructor::ComputeAvgVolume() const
 {
   MFEM_VERIFY(nodes, "Nodes are not given!");
-   MFEM_ASSERT(avg_volume == 0.0, "The average volume is already computed!");
+   MFEM_ASSERT(avg_volume == 0.0, "the average volume is already computed!");

   Mesh *mesh = nodes->FESpace()->GetMesh();
   const int NE = mesh->GetNE();
@@ -788,13 +787,9 @@ void TargetConstructor::ComputeAvgVolume() const
         volume += ip.weight * Tr.Weight();
      }
   }
-
-   NCMesh *ncmesh = mesh->ncmesh;
-   if (Parallel() == false)
+   if (!Parallel())
   {
-      avg_volume = (ncmesh == NULL) ?
-                   volume / NE : volume / ncmesh->GetNumRootElements();
-
+      avg_volume = volume / NE;
   }
 #ifdef MFEM_USE_MPI
   else
@@ -802,8 +797,7 @@ void TargetConstructor::ComputeAvgVolume() const
      double area_NE[4];
      area_NE[0] = volume; area_NE[1] = NE;
      MPI_Allreduce(area_NE, area_NE + 2, 2, MPI_DOUBLE, MPI_SUM, comm);
-      avg_volume = (ncmesh == NULL) ?
-                   area_NE[2] / area_NE[3] : area_NE[2] / ncmesh->GetNumRootElements();
+      avg_volume = area_NE[2] / area_NE[3];
   }
 #endif
 }
@@ -811,7 +805,6 @@ void TargetConstructor::ComputeAvgVolume() const
 // virtual method
 void TargetConstructor::ComputeElementTargets(int e_id, const FiniteElement &fe,
                                              const IntegrationRule &ir,
-                                              const Vector &elfun,
                                              DenseTensor &Jtr) const
 {
   MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
@@ -834,15 +827,7 @@ void TargetConstructor::ComputeElementTargets(int e_id, const FiniteElement &fe,
      {
         if (avg_volume == 0.0) { ComputeAvgVolume(); }
         DenseMatrix W(Wideal.Height());
-
-         NCMesh *ncmesh = nodes->FESpace()->GetMesh()->ncmesh;
-         double el_volume = avg_volume;
-         if (ncmesh)
-         {
-            el_volume = avg_volume / ncmesh->GetElementSizeReduction(e_id);
-         }
-
-         W.Set(std::pow(volume_scale * el_volume / Wideal.Det(),
+         W.Set(std::pow(volume_scale * avg_volume / Wideal.Det(),
                        1./W.Height()), Wideal);
         for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = W; }
         break;
@@ -868,7 +853,7 @@ void TargetConstructor::ComputeElementTargets(int e_id, const FiniteElement &fe,
            if (target_type == IDEAL_SHAPE_GIVEN_SIZE)
            {
               const double det = Jtr(i).Det();
-               MFEM_VERIFY(det > 0.0, "The given mesh is inverted!");
+               MFEM_VERIFY(det > 0.0, "Initial mesh is inverted!");
               Jtr(i).Set(std::pow(det / detW, 1./dim), Wideal);
            }
         }
@@ -879,162 +864,6 @@ void TargetConstructor::ComputeElementTargets(int e_id, const FiniteElement &fe,
   }
 }

-void AnalyticAdaptTC::SetAnalyticTargetSpec(Coefficient *sspec,
-                                            VectorCoefficient *vspec,
-                                            MatrixCoefficient *mspec)
-{
-   scalar_tspec = sspec;
-   vector_tspec = vspec;
-   matrix_tspec = mspec;
-}
-
-void AnalyticAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
-                                            const IntegrationRule &ir,
-                                            const Vector &elfun,
-                                            DenseTensor &Jtr) const
-{
-   DenseMatrix point_mat;
-   point_mat.UseExternalData(elfun.GetData(), fe.GetDof(), fe.GetDim());
-
-   switch (target_type)
-   {
-      case GIVEN_FULL:
-      {
-         MFEM_VERIFY(matrix_tspec != NULL,
-                     "Target type GIVEN_FULL requires a MatrixCoefficient.");
-
-         IsoparametricTransformation Tpr;
-         Tpr.SetFE(&fe);
-         Tpr.ElementNo = e_id;
-         Tpr.GetPointMat().Transpose(point_mat);
-
-         for (int i = 0; i < ir.GetNPoints(); i++)
-         {
-            const IntegrationPoint &ip = ir.IntPoint(i);
-            Tpr.SetIntPoint(&ip);
-            matrix_tspec->Eval(Jtr(i), Tpr, ip);
-         }
-         break;
-      }
-      default:
-         MFEM_ABORT("Incompatible target type for analytic adaptation!");
-   }
-}
-
-#ifdef MFEM_USE_MPI
-void DiscreteAdaptTC::SetParDiscreteTargetSpec(ParGridFunction &tspec)
-{
-   target_spec.SetSize(tspec.Size());
-   target_spec = tspec;
-   tspec_fes   = tspec.FESpace();
-
-   // Default evaluator is based on CG advection.
-   if (adapt_eval == NULL) { adapt_eval = new AdvectorCG; }
-
-   adapt_eval->SetParMetaInfo(*tspec.ParFESpace()->GetParMesh(),
-                              *tspec.FESpace()->FEColl(),
-                              tspec.FESpace()->GetVDim());
-
-   adapt_eval->SetInitialField
-   (*tspec.FESpace()->GetMesh()->GetNodes(), target_spec);
-}
-#endif
-
-void DiscreteAdaptTC::SetSerialDiscreteTargetSpec(GridFunction &tspec)
-{
-   target_spec.SetSize(tspec.Size());
-   target_spec = tspec;
-   tspec_fes   = tspec.FESpace();
-
-   // Default evaluator is based on CG advection.
-   if (adapt_eval == NULL) { adapt_eval = new AdvectorCG; }
-
-   adapt_eval->SetSerialMetaInfo(*tspec.FESpace()->GetMesh(),
-                                 *tspec.FESpace()->FEColl(),
-                                 tspec.FESpace()->GetVDim());
-
-   adapt_eval->SetInitialField
-   (*tspec.FESpace()->GetMesh()->GetNodes(), target_spec);
-}
-
-void DiscreteAdaptTC::UpdateTargetSpecification(const Vector &new_x)
-{
-   MFEM_VERIFY(target_spec.Size() > 0, "Target specification is not set!");
-
-   adapt_eval->ComputeAtNewPosition(new_x, target_spec);
-}
-
-void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
-                                            const IntegrationRule &ir,
-                                            const Vector &elfun,
-                                            DenseTensor &Jtr) const
-{
-   MFEM_VERIFY(tspec_fes, "A call to SetDiscreteTargerSpec() is needed.");
-
-   switch (target_type)
-   {
-      case IDEAL_SHAPE_GIVEN_SIZE:
-      {
-         const DenseMatrix &Wideal =
-            Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
-         const int dim = Wideal.Height(),
-                   ntspec_dofs = tspec_fes->GetFE(0)->GetDof();
-
-         Vector shape(ntspec_dofs), tspec_vals(ntspec_dofs);
-         Array<int> dofs;
-         tspec_fes->GetElementDofs(e_id, dofs);
-         target_spec.GetSubVector(dofs, tspec_vals);
-
-         const double min_size = tspec_vals.Min();
-         MFEM_ASSERT(min_size > 0.0,
-                     "Non-positive size propagated in the target definition.");
-
-         for (int i = 0; i < ir.GetNPoints(); i++)
-         {
-            const IntegrationPoint &ip = ir.IntPoint(i);
-            tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
-            const double size = std::max(shape * tspec_vals, min_size);
-            Jtr(i).Set(std::pow(size / Wideal.Det(), 1.0/dim), Wideal);
-         }
-         break;
-      }
-      default:
-         MFEM_ABORT("Incompatible target type for analytic adaptation!");
-   }
-}
-
-void AdaptivityEvaluator::SetSerialMetaInfo(const Mesh &m,
-                                            const FiniteElementCollection &fec,
-                                            int num_comp)
-{
-   delete fes;
-   delete mesh;
-   mesh = new Mesh(m, true);
-   fes = new FiniteElementSpace(mesh, &fec, num_comp);
-}
-
-#ifdef MFEM_USE_MPI
-void AdaptivityEvaluator::SetParMetaInfo(const ParMesh &m,
-                                         const FiniteElementCollection &fec,
-                                         int num_comp)
-{
-   delete pfes;
-   delete pmesh;
-   pmesh = new ParMesh(m, true);
-   pfes  = new ParFiniteElementSpace(pmesh, &fec, num_comp);
-}
-#endif
-
-AdaptivityEvaluator::~AdaptivityEvaluator()
-{
-   delete fes;
-   delete mesh;
-#ifdef MFEM_USE_MPI
-   delete pfes;
-   delete pmesh;
-#endif
-}
-
 void TMOP_Integrator::EnableLimiting(const GridFunction &n0,
                                     const GridFunction &dist, Coefficient &w0,
                                     TMOP_LimiterFunction *lfunc)
@@ -1092,7 +921,7 @@ double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,

   energy = 0.0;
   DenseTensor Jtr(dim, dim, ir->GetNPoints());
-   targetC->ComputeElementTargets(T.ElementNo, el, *ir, elfun, Jtr);
+   targetC->ComputeElementTargets(T.ElementNo, el, *ir, Jtr);

   // Limited case.
   Vector shape, p, p0, d_vals;
@@ -1161,7 +990,6 @@ double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
      energy += weight * val;
   }
   delete Tpr;
-
   return energy;
 }

@@ -1188,7 +1016,7 @@ void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,

   elvect = 0.0;
   DenseTensor Jtr(dim, dim, ir->GetNPoints());
-   targetC->ComputeElementTargets(T.ElementNo, el, *ir, elfun, Jtr);
+   targetC->ComputeElementTargets(T.ElementNo, el, *ir, Jtr);

   // Limited case.
   DenseMatrix pos0;
@@ -1244,8 +1072,6 @@ void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
      P *= weight_m;
      AddMultABt(DS, P, PMatO);

-      // TODO: derivatives of adaptivity-based targets.
-
      if (coeff0)
      {
         el.CalcShape(ip, shape);
@@ -1281,7 +1107,7 @@ void TMOP_Integrator::AssembleElementGrad(const FiniteElement &el,

   elmat = 0.0;
   DenseTensor Jtr(dim, dim, ir->GetNPoints());
-   targetC->ComputeElementTargets(T.ElementNo, el, *ir, elfun, Jtr);
+   targetC->ComputeElementTargets(T.ElementNo, el, *ir, Jtr);

   // Limited case.
   DenseMatrix pos0, grad_grad;
@@ -1334,8 +1160,6 @@ void TMOP_Integrator::AssembleElementGrad(const FiniteElement &el,

      metric->AssembleH(Jpt, DS, weight_m, elmat);

-      // TODO: derivatives of adaptivity-based targets.
-
      if (coeff0)
      {
         el.CalcShape(ip, shape);
@@ -1410,12 +1234,11 @@ void TMOP_Integrator::ComputeNormalizationEnergies(const GridFunction &x,
   for (int i = 0; i < fes->GetNE(); i++)
   {
      fe = fes->GetFE(i);
+      targetC->ComputeElementTargets(i, *fe, *ir, Jtr);
      fes->GetElementVDofs(i, vdofs);
      x.GetSubVector(vdofs, x_vals);
      PMatI.UseExternalData(x_vals.GetData(), dof, dim);

-      targetC->ComputeElementTargets(i, *fe, *ir, x_vals, Jtr);
-
      for (int i = 0; i < ir->GetNPoints(); i++)
      {
         const IntegrationPoint &ip = ir->IntPoint(i);
@@ -1451,6 +1274,9 @@ void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric,
      const IntegrationRule &ir = metric_gf.FESpace()->GetFE(i)->GetNodes();
      const int nsp = ir.GetNPoints(), dof = fe_pos.GetDof();

+      W.SetSize(dim, dim, nsp);
+      tc.ComputeElementTargets(i, fe_pos, ir, W);
+
      dshape.SetSize(dof, dim);
      pos.SetSize(dof, dim);
      posV.SetDataAndSize(pos.Data(), dof * dim);
@@ -1459,9 +1285,6 @@ void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric,
      nodes.FESpace()->GetElementVDofs(i, pos_dofs);
      nodes.GetSubVector(pos_dofs, posV);

-      W.SetSize(dim, dim, nsp);
-      tc.ComputeElementTargets(i, fe_pos, ir, posV, W);
-
      for (int j = 0; j < nsp; j++)
      {
         const DenseMatrix &Wj = W(j);
@@ -12,6 +12,7 @@
 #ifndef MFEM_TMOP_HPP
 #define MFEM_TMOP_HPP

+#include "../config/config.hpp"
 #include "../linalg/invariants.hpp"
 #include "nonlininteg.hpp"

@@ -513,51 +514,6 @@ public:
   virtual ~TMOP_QuadraticLimiter() { }
 };

-class FiniteElementCollection;
-class FiniteElementSpace;
-class ParFiniteElementSpace;
-
-class AdaptivityEvaluator
-{
-protected:
-   // Owned.
-   Mesh *mesh;
-   FiniteElementSpace *fes;
-
-#ifdef MFEM_USE_MPI
-   // Owned.
-   ParMesh *pmesh;
-   ParFiniteElementSpace *pfes;
-#endif
-
-public:
-   AdaptivityEvaluator() : mesh(NULL), fes(NULL)
-   {
-#ifdef MFEM_USE_MPI
-      pmesh = NULL;
-      pfes = NULL;
-#endif
-   }
-   virtual ~AdaptivityEvaluator();
-
-   /** Specifies the Mesh and FiniteElementCollection of the solution that will
-       be evaluated. The given mesh will be copied into the internal object. */
-   void SetSerialMetaInfo(const Mesh &m,
-                          const FiniteElementCollection &fec, int num_comp);
-
-#ifdef MFEM_USE_MPI
-   /// Parallel version of SetSerialMetaInfo.
-   void SetParMetaInfo(const ParMesh &m,
-                       const FiniteElementCollection &fec, int num_comp);
-#endif
-
-   // TODO use GridFunctions to make clear it's on the ldofs?
-   virtual void SetInitialField(const Vector &init_nodes,
-                                const Vector &init_field) = 0;
-
-   virtual void ComputeAtNewPosition(const Vector &new_nodes,
-                                     Vector &new_field) = 0;
-};

 /** @brief Base class representing target-matrix construction algorithms for
    mesh optimization via the target-matrix optimization paradigm (TMOP). */
@@ -582,11 +538,9 @@ public:
      IDEAL_SHAPE_GIVEN_SIZE, /**<
         Ideal shape, given size/volume; the given nodes define the target
         volume at all quadrature points. */
-      GIVEN_SHAPE_AND_SIZE, /**<
+      GIVEN_SHAPE_AND_SIZE /**<
         Given shape, given size/volume; the given nodes define the exact target
         Jacobian matrix at all quadrature points. */
-      GIVEN_FULL /**<
-         Full target tensor is specified at every quadrature point. */
   };

 protected:
@@ -635,89 +589,14 @@ public:
   void SetVolumeScale(double vol_scale) { volume_scale = vol_scale; }

   /** @brief Given an element and quadrature rule, computes ref->target
-       transformation Jacobians for each quadrature point in the element.
-       The physical positions of the element's nodes are given by @a elfun. */
+       transformation Jacobians for each quadrature point in the element. */
   virtual void ComputeElementTargets(int e_id, const FiniteElement &fe,
                                      const IntegrationRule &ir,
-                                      const Vector &elfun,
-                                      DenseTensor &Jtr) const;
-};
-
-class AnalyticAdaptTC : public TargetConstructor
-{
-protected:
-   // Analytic target specification.
-   Coefficient *scalar_tspec;
-   VectorCoefficient *vector_tspec;
-   MatrixCoefficient *matrix_tspec;
-
-public:
-   AnalyticAdaptTC(TargetType ttype)
-      : TargetConstructor(ttype),
-        scalar_tspec(NULL), vector_tspec(NULL), matrix_tspec(NULL) { }
-
-   virtual void SetAnalyticTargetSpec(Coefficient *sspec,
-                                      VectorCoefficient *vspec,
-                                      MatrixCoefficient *mspec);
-
-   /** @brief Given an element and quadrature rule, computes ref->target
-       transformation Jacobians for each quadrature point in the element.
-       The physical positions of the element's nodes are given by @a elfun. */
-   virtual void ComputeElementTargets(int e_id, const FiniteElement &fe,
-                                      const IntegrationRule &ir,
-                                      const Vector &elfun,
                                      DenseTensor &Jtr) const;
 };

 class ParGridFunction;

-class DiscreteAdaptTC : public TargetConstructor
-{
-protected:
-   // Discrete target specification.
-   // Data is owned, updated by UpdateTargetSpecification.
-   Vector target_spec;
-   // Note: do not use the Nodes of this space as they may not be on the
-   // positions corresponding to the values of tspec.
-   const FiniteElementSpace *tspec_fes;
-
-   // Evaluation of the discrete target specification on different meshes.
-   // Owned.
-   AdaptivityEvaluator *adapt_eval;
-
-public:
-   DiscreteAdaptTC(TargetType ttype)
-      : TargetConstructor(ttype),
-        target_spec(), tspec_fes(NULL), adapt_eval(NULL) { }
-
-   virtual ~DiscreteAdaptTC() { delete adapt_eval; }
-
-   virtual void SetSerialDiscreteTargetSpec(GridFunction &tspec);
-#ifdef MFEM_USE_MPI
-   virtual void SetParDiscreteTargetSpec(ParGridFunction &tspec);
-#endif
-
-   /** Used to update the target specification after the mesh has changed. The
-       new mesh positions are given by new_x. */
-   void UpdateTargetSpecification(const Vector &new_x);
-
-   void SetAdaptivityEvaluator(AdaptivityEvaluator *ae)
-   {
-      if (adapt_eval) { delete adapt_eval; }
-      adapt_eval = ae;
-   }
-
-   /** @brief Given an element and quadrature rule, computes ref->target
-       transformation Jacobians for each quadrature point in the element.
-       The physical positions of the element's nodes are given by @a elfun.
-       Note that this function assumes that UpdateTargetSpecification() has
-       been called with the position vector corresponding to @a elfun. */
-   virtual void ComputeElementTargets(int e_id, const FiniteElement &fe,
-                                      const IntegrationRule &ir,
-                                      const Vector &elfun,
-                                      DenseTensor &Jtr) const;
-};
-
 /** @brief A TMOP integrator class based on any given TMOP_QualityMetric and
    TargetConstructor.

@@ -1,518 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM 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) version 2.1 dated February 1999.
-
-#include "tmop_tools.hpp"
-#include "nonlinearform.hpp"
-#include "pnonlinearform.hpp"
-#include "../general/osockstream.hpp"
-
-namespace mfem
-{
-
-using namespace mfem;
-
-void AdvectorCG::SetInitialField(const Vector &init_nodes,
-                                 const Vector &init_field)
-{
-   nodes0 = init_nodes;
-   field0 = init_field;
-}
-
-void AdvectorCG::ComputeAtNewPosition(const Vector &new_nodes,
-                                      Vector &new_field)
-{
-   int myid = 0;
-   Mesh *m = mesh;
-
-#ifdef MFEM_USE_MPI
-   if (pfes) { MPI_Comm_rank(pfes->GetComm(), &myid); }
-   if (pmesh) { m = pmesh; }
-#endif
-
-   MFEM_VERIFY(m != NULL, "No mesh has been given to the AdaptivityEvaluator.");
-
-   // This will be used to move the positions.
-   GridFunction *mesh_nodes = m->GetNodes();
-   *mesh_nodes = nodes0;
-   new_field = field0;
-
-   // Velocity of the positions.
-   GridFunction u(mesh_nodes->FESpace());
-   subtract(new_nodes, nodes0, u);
-
-   TimeDependentOperator *oper = NULL;
-   // This must be the fes of the ind, associated with the object's mesh.
-   if (fes)  { oper = new SerialAdvectorCGOper(nodes0, u, *fes); }
-#ifdef MFEM_USE_MPI
-   else if (pfes) { oper = new ParAdvectorCGOper(nodes0, u, *pfes); }
-#endif
-   MFEM_VERIFY(oper != NULL,
-               "No FE space has been given to the AdaptivityEvaluator.");
-   ode_solver.Init(*oper);
-
-   // Compute some time step [mesh_size / speed].
-   double min_h = std::numeric_limits<double>::infinity();
-   for (int i = 0; i < m->GetNE(); i++)
-   {
-      min_h = std::min(min_h, m->GetElementSize(i));
-   }
-   double v_max = 0.0;
-   const int s = u.FESpace()->GetVSize() / 2;
-   for (int i = 0; i < s; i++)
-   {
-      const double vel = u(i) * u(i) + u(i+s) * u(i+s);
-      v_max = std::max(v_max, vel);
-   }
-   if (v_max == 0.0)
-   {
-      // No need to change the field.
-      return;
-   }
-   v_max = std::sqrt(v_max);
-   double dt = 0.5 * min_h / v_max;
-   double glob_dt = dt;
-#ifdef MFEM_USE_MPI
-   if (pfes)
-   {
-      MPI_Allreduce(&dt, &glob_dt, 1, MPI_DOUBLE, MPI_MIN, pfes->GetComm());
-   }
-#endif
-
-   double t = 0.0;
-   bool last_step = false;
-   for (int ti = 1; !last_step; ti++)
-   {
-      if (t + glob_dt >= 1.0)
-      {
-#ifdef MFEM_DEBUG
-         if (myid == 0)
-         {
-            mfem::out << "Remap took " << ti << " steps." << std::endl;
-         }
-#endif
-         glob_dt = 1.0 - t;
-         last_step = true;
-      }
-      ode_solver.Step(new_field, t, glob_dt);
-   }
-
-   // Trim the overshoots and undershoots.
-   const double minv = field0.Min(), maxv = field0.Max();
-   for (int i = 0; i < new_field.Size(); i++)
-   {
-      if (new_field(i) < minv) { new_field(i) = minv; }
-      if (new_field(i) > maxv) { new_field(i) = maxv; }
-   }
-
-   nodes0 = new_nodes;
-   field0 = new_field;
-
-   delete oper;
-}
-
-SerialAdvectorCGOper::SerialAdvectorCGOper(const Vector &x_start,
-                                           GridFunction &vel,
-                                           FiniteElementSpace &fes)
-   : TimeDependentOperator(fes.GetVSize()),
-     x0(x_start), x_now(*fes.GetMesh()->GetNodes()),
-     u(vel), u_coeff(&u), M(&fes), K(&fes)
-{
-   ConvectionIntegrator *Kinteg = new ConvectionIntegrator(u_coeff);
-   K.AddDomainIntegrator(Kinteg);
-   K.Assemble(0);
-   K.Finalize(0);
-
-   MassIntegrator *Minteg = new MassIntegrator;
-   M.AddDomainIntegrator(Minteg);
-   M.Assemble();
-   M.Finalize();
-}
-
-void SerialAdvectorCGOper::Mult(const Vector &ind, Vector &di_dt) const
-{
-   // Move the mesh.
-   const double t = GetTime();
-   add(x0, t, u, x_now);
-
-   // Assemble on the new mesh.
-   K.BilinearForm::operator=(0.0);
-   K.Assemble();
-   Vector rhs(K.Size());
-   K.Mult(ind, rhs);
-   M.BilinearForm::operator=(0.0);
-   M.Assemble();
-
-   di_dt = 0.0;
-   CGSolver lin_solver;
-   DSmoother prec;
-   lin_solver.SetPreconditioner(prec);
-   lin_solver.SetOperator(M.SpMat());
-   lin_solver.SetRelTol(1e-12); lin_solver.SetAbsTol(0.0);
-   lin_solver.SetMaxIter(100);
-   lin_solver.SetPrintLevel(0);
-   lin_solver.Mult(rhs, di_dt);
-}
-
-#ifdef MFEM_USE_MPI
-ParAdvectorCGOper::ParAdvectorCGOper(const Vector &x_start,
-                                     GridFunction &vel,
-                                     ParFiniteElementSpace &pfes)
-   : TimeDependentOperator(pfes.GetVSize()),
-     x0(x_start), x_now(*pfes.GetMesh()->GetNodes()),
-     u(vel), u_coeff(&u), M(&pfes), K(&pfes)
-{
-   ConvectionIntegrator *Kinteg = new ConvectionIntegrator(u_coeff);
-   K.AddDomainIntegrator(Kinteg);
-   K.Assemble(0);
-   K.Finalize(0);
-
-   MassIntegrator *Minteg = new MassIntegrator;
-   M.AddDomainIntegrator(Minteg);
-   M.Assemble();
-   M.Finalize();
-}
-
-void ParAdvectorCGOper::Mult(const Vector &ind, Vector &di_dt) const
-{
-   // Move the mesh.
-   const double t = GetTime();
-   add(x0, t, u, x_now);
-
-   // Assemble on the new mesh.
-   K.BilinearForm::operator=(0.0);
-   K.Assemble();
-   ParGridFunction rhs(K.ParFESpace());
-   K.Mult(ind, rhs);
-   M.BilinearForm::operator=(0.0);
-   M.Assemble();
-
-   HypreParVector *RHS = rhs.ParallelAssemble();
-   HypreParVector X(K.ParFESpace());
-   X = 0.0;
-   HypreParMatrix *Mh  = M.ParallelAssemble();
-
-   CGSolver lin_solver(M.ParFESpace()->GetParMesh()->GetComm());
-   HypreSmoother prec;
-   prec.SetType(HypreSmoother::Jacobi, 1);
-   lin_solver.SetPreconditioner(prec);
-   lin_solver.SetOperator(*Mh);
-   lin_solver.SetRelTol(1e-8);
-   lin_solver.SetAbsTol(0.0);
-   lin_solver.SetMaxIter(100);
-   lin_solver.SetPrintLevel(0);
-   lin_solver.Mult(*RHS, X);
-   K.ParFESpace()->GetProlongationMatrix()->Mult(X, di_dt);
-
-   delete Mh;
-   delete RHS;
-}
-#endif
-
-double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
-                                              const Vector &b) const
-{
-   const FiniteElementSpace *fes = NULL;
-   double energy_in = 0.0;
-#ifdef MFEM_USE_MPI
-   const ParNonlinearForm *p_nlf = dynamic_cast<const ParNonlinearForm *>(oper);
-   MFEM_VERIFY(!(parallel && p_nlf == NULL), "Invalid Operator subclass.");
-   if (parallel)
-   {
-      fes = p_nlf->FESpace();
-      energy_in = p_nlf->GetEnergy(x);
-   }
-#endif
-   const bool serial = !parallel;
-   const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
-   MFEM_VERIFY(!(serial && nlf == NULL), "Invalid Operator subclass.");
-   if (serial)
-   {
-      fes = nlf->FESpace();
-      energy_in = nlf->GetEnergy(x);
-   }
-
-   const bool have_b = (b.Size() == Height());
-
-   const int NE = fes->GetMesh()->GetNE(), dim = fes->GetFE(0)->GetDim(),
-             dof = fes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
-   Array<int> xdofs(dof * dim);
-   DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
-   Vector posV(pos.Data(), dof * dim);
-
-   Vector x_out(x.Size()), x_out_loc(fes->GetVSize());
-   bool x_out_ok = false;
-   double scale = 1.0, energy_out;
-   double norm0 = Norm(r);
-
-   // Decreases the scaling of the update until the new mesh is valid.
-   for (int i = 0; i < 12; i++)
-   {
-      add(x, -scale, c, x_out);
-
-      if (serial)
-      {
-         const SparseMatrix *cP = fes->GetConformingProlongation();
-         if (!cP) {x_out_loc.SetData(x_out.GetData());}
-         else {cP->Mult(x_out,x_out_loc);}
-         energy_out = nlf->GetGridFunctionEnergy(x_out_loc);
-      }
-#ifdef MFEM_USE_MPI
-      else
-      {
-         fes->GetProlongationMatrix()->Mult(x_out, x_out_loc);
-         energy_out = p_nlf->GetParGridFunctionEnergy(x_out_loc);
-      }
-#endif
-
-      if (energy_out > 1.2*energy_in || std::isnan(energy_out) != 0)
-      {
-         if (print_level >= 0)
-         { mfem::out << "Scale = " << scale << " Increasing energy.\n"; }
-         scale *= 0.5; continue;
-      }
-
-      int jac_ok = 1;
-      for (int i = 0; i < NE; i++)
-      {
-         fes->GetElementVDofs(i, xdofs);
-         x_out_loc.GetSubVector(xdofs, posV);
-         for (int j = 0; j < nsp; j++)
-         {
-            fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
-            MultAtB(pos, dshape, Jpr);
-            if (Jpr.Det() <= 0.0) { jac_ok = 0; goto break2; }
-         }
-      }
-   break2:
-      int jac_ok_all = jac_ok;
-#ifdef MFEM_USE_MPI
-      if (parallel)
-      {
-         MPI_Allreduce(&jac_ok, &jac_ok_all, 1, MPI_INT, MPI_LAND,
-                       p_nlf->ParFESpace()->GetComm());
-      }
-#endif
-
-      if (jac_ok_all == 0)
-      {
-         if (print_level >= 0)
-         { mfem::out << "Scale = " << scale << " Neg det(J) found.\n"; }
-         scale *= 0.5; continue;
-      }
-
-      oper->Mult(x_out, r);
-      if (have_b) { r -= b; }
-      double norm = Norm(r);
-
-      if (norm > 1.2*norm0)
-      {
-         if (print_level >= 0)
-         { mfem::out << "Scale = " << scale << " Norm increased.\n"; }
-         scale *= 0.5; continue;
-      }
-      else { x_out_ok = true; break; }
-   }
-
-   if (print_level >= 0)
-   {
-      mfem::out << "Energy decrease: "
-                << (energy_in - energy_out) / energy_in * 100.0
-                << "% with " << scale << " scaling.\n";
-   }
-
-   if (x_out_ok == false) { scale = 0.0; }
-   return scale;
-}
-
-void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
-{
-   if (discr_tc)
-   {
-      if (parallel)
-      {
-#ifdef MFEM_USE_MPI
-         const ParNonlinearForm *nlf =
-            dynamic_cast<const ParNonlinearForm *>(oper);
-         Vector x_loc(nlf->ParFESpace()->GetVSize());
-         nlf->ParFESpace()->GetProlongationMatrix()->Mult(x, x_loc);
-         discr_tc->UpdateTargetSpecification(x_loc);
-#endif
-      }
-      else { discr_tc->UpdateTargetSpecification(x); }
-   }
-}
-
-double TMOPDescentNewtonSolver::ComputeScalingFactor(const Vector &x,
-                                                     const Vector &b) const
-{
-   const FiniteElementSpace *fes = NULL;
-   double energy_in = 0.0;
-#ifdef MFEM_USE_MPI
-   const ParNonlinearForm *p_nlf = dynamic_cast<const ParNonlinearForm *>(oper);
-   MFEM_VERIFY(!(parallel && p_nlf == NULL), "Invalid Operator subclass.");
-   if (parallel)
-   {
-      fes = p_nlf->FESpace();
-      energy_in = p_nlf->GetEnergy(x);
-   }
-#endif
-   const bool serial = !parallel;
-   const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
-   MFEM_VERIFY(!(serial && nlf == NULL), "Invalid Operator subclass.");
-   if (serial)
-   {
-      fes = nlf->FESpace();
-      energy_in = nlf->GetEnergy(x);
-   }
-
-   const int NE = fes->GetMesh()->GetNE(), dim = fes->GetFE(0)->GetDim(),
-             dof = fes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
-   Array<int> xdofs(dof * dim);
-   DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
-   Vector posV(pos.Data(), dof * dim);
-   Vector x_loc(fes->GetVSize());
-
-   double min_detJ = infinity();
-   for (int i = 0; i < NE; i++)
-   {
-      fes->GetElementVDofs(i, xdofs);
-      x_loc.GetSubVector(xdofs, posV);
-
-      for (int j = 0; j < nsp; j++)
-      {
-         fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
-         MultAtB(pos, dshape, Jpr);
-         min_detJ = std::min(min_detJ, Jpr.Det());
-      }
-   }
-   double min_detJ_all = min_detJ;
-#ifdef MFEM_USE_MPI
-   if (parallel)
-   {
-      MPI_Allreduce(&min_detJ, &min_detJ_all, 1, MPI_DOUBLE, MPI_MIN,
-                    p_nlf->ParFESpace()->GetComm());
-   }
-#endif
-   if (print_level >= 0)
-   {
-      mfem::out << "Minimum det(J) = " << min_detJ_all << '\n';
-   }
-
-   Vector x_out(x.Size());
-   bool x_out_ok = false;
-   double scale = 1.0, energy_out;
-
-   for (int i = 0; i < 7; i++)
-   {
-      add(x, -scale, c, x_out);
-      if (serial)
-      {
-         const SparseMatrix *cP = fes->GetConformingProlongation();
-         if (!cP) {x_loc.SetData(x_out.GetData());}
-         else {cP->Mult(x_out,x_loc);}
-         energy_out = nlf->GetGridFunctionEnergy(x_loc);
-      }
-#ifdef MFEM_USE_MPI
-      else
-      {
-         fes->GetProlongationMatrix()->Mult(x_out, x_loc);
-         energy_out = p_nlf->GetParGridFunctionEnergy(x_loc);
-      }
-#endif
-
-      if (energy_out > energy_in || std::isnan(energy_out) != 0)
-      {
-         scale *= 0.5;
-      }
-      else { x_out_ok = true; break; }
-   }
-
-   if (print_level >= 0)
-   {
-      mfem::out << "Energy decrease: "
-                << (energy_in - energy_out) / energy_in * 100.0
-                << "% with " << scale << " scaling.\n";
-   }
-
-   if (x_out_ok == false) { return 0.0; }
-
-   return scale;
-}
-
-void TMOPDescentNewtonSolver::ProcessNewState(const Vector &x) const
-{
-   if (discr_tc)
-   {
-      if (parallel)
-      {
-#ifdef MFEM_USE_MPI
-         const ParNonlinearForm *nlf =
-            dynamic_cast<const ParNonlinearForm *>(oper);
-         Vector x_loc(nlf->ParFESpace()->GetVSize());
-         nlf->ParFESpace()->GetProlongationMatrix()->Mult(x, x_loc);
-         discr_tc->UpdateTargetSpecification(x_loc);
-#endif
-      }
-      else { discr_tc->UpdateTargetSpecification(x); }
-   }
-}
-
-#ifdef MFEM_USE_MPI
-// Metric values are visualized by creating an L2 finite element functions and
-// computing the metric values at the nodes.
-void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm,
-                       const TargetConstructor &tc, ParMesh &pmesh,
-                       char *title, int position)
-{
-   L2_FECollection fec(order, pmesh.Dimension(), BasisType::GaussLobatto);
-   ParFiniteElementSpace fes(&pmesh, &fec, 1);
-   ParGridFunction metric(&fes);
-   InterpolateTMOP_QualityMetric(qm, tc, pmesh, metric);
-   socketstream sock;
-   if (pmesh.GetMyRank() == 0)
-   {
-      sock.open("localhost", 19916);
-      sock << "solution\n";
-   }
-   pmesh.PrintAsOne(sock);
-   metric.SaveAsOne(sock);
-   if (pmesh.GetMyRank() == 0)
-   {
-      sock << "window_title '"<< title << "'\n"
-           << "window_geometry "
-           << position << " " << 0 << " " << 600 << " " << 600 << "\n"
-           << "keys jRmclA\n";
-   }
-}
-#endif
-
-// Metric values are visualized by creating an L2 finite element functions and
-// computing the metric values at the nodes.
-void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm,
-                       const TargetConstructor &tc, Mesh &mesh,
-                       char *title, int position)
-{
-   L2_FECollection fec(order, mesh.Dimension(), BasisType::GaussLobatto);
-   FiniteElementSpace fes(&mesh, &fec, 1);
-   GridFunction metric(&fes);
-   InterpolateTMOP_QualityMetric(qm, tc, mesh, metric);
-   osockstream sock(19916, "localhost");
-   sock << "solution\n";
-   mesh.Print(sock);
-   metric.Save(sock);
-   sock.send();
-   sock << "window_title '"<< title << "'\n"
-        << "window_geometry "
-        << position << " " << 0 << " " << 600 << " " << 600 << "\n"
-        << "keys jRmclA\n";
-}
-
-}
@@ -1,140 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM 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) version 2.1 dated February 1999.
-
-#ifndef MFEM_TMOP_TOOLS_HPP
-#define MFEM_TMOP_TOOLS_HPP
-
-#include "bilinearform.hpp"
-#include "pbilinearform.hpp"
-#include "tmop.hpp"
-
-namespace mfem
-{
-
-// Performs the full remap advection loop.
-class AdvectorCG : public AdaptivityEvaluator
-{
-private:
-   RK4Solver ode_solver;
-   Vector nodes0;
-   Vector field0;
-
-public:
-   AdvectorCG() : AdaptivityEvaluator(), ode_solver(), nodes0(), field0() { }
-
-   virtual void SetInitialField(const Vector &init_nodes,
-                                const Vector &init_field);
-
-   virtual void ComputeAtNewPosition(const Vector &new_nodes,
-                                     Vector &new_field);
-};
-
-/// Performs a single remap advection step in serial.
-class SerialAdvectorCGOper : public TimeDependentOperator
-{
-protected:
-   const Vector &x0;
-   Vector &x_now;
-   GridFunction &u;
-   VectorGridFunctionCoefficient u_coeff;
-   mutable BilinearForm M, K;
-
-public:
-   /** Here @a fes is the FESpace of the function that will be moved. Note
-       that Mult() moves the nodes of the mesh corresponding to @a fes. */
-   SerialAdvectorCGOper(const Vector &x_start, GridFunction &vel,
-                        FiniteElementSpace &fes);
-
-   virtual void Mult(const Vector &ind, Vector &di_dt) const;
-};
-
-#ifdef MFEM_USE_MPI
-/// Performs a single remap advection step in parallel.
-class ParAdvectorCGOper : public TimeDependentOperator
-{
-protected:
-   const Vector &x0;
-   Vector &x_now;
-   GridFunction &u;
-   VectorGridFunctionCoefficient u_coeff;
-   mutable ParBilinearForm M, K;
-
-public:
-   /** Here @a pfes is the ParFESpace of the function that will be moved. Note
-       that Mult() moves the nodes of the mesh corresponding to @a pfes. */
-   ParAdvectorCGOper(const Vector &x_start, GridFunction &vel,
-                     ParFiniteElementSpace &pfes);
-
-   virtual void Mult(const Vector &ind, Vector &di_dt) const;
-};
-#endif
-
-class TMOPNewtonSolver : public NewtonSolver
-{
-private:
-   bool parallel;
-
-   // Quadrature points that are checked for negative Jacobians etc.
-   const IntegrationRule &ir;
-
-   mutable DiscreteAdaptTC *discr_tc;
-
-public:
-#ifdef MFEM_USE_MPI
-   TMOPNewtonSolver(MPI_Comm comm, const IntegrationRule &irule)
-      : NewtonSolver(comm), parallel(true), ir(irule), discr_tc(NULL) { }
-#endif
-   TMOPNewtonSolver(const IntegrationRule &irule)
-      : NewtonSolver(), parallel(false), ir(irule), discr_tc(NULL) { }
-
-   void SetDiscreteAdaptTC(DiscreteAdaptTC *tc) { discr_tc = tc; }
-
-   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
-
-   virtual void ProcessNewState(const Vector &x) const;
-};
-
-/// Allows negative Jacobians. Used for untangling.
-class TMOPDescentNewtonSolver : public NewtonSolver
-{
-private:
-   bool parallel;
-
-   // Quadrature points that are checked for negative Jacobians etc.
-   const IntegrationRule &ir;
-
-   mutable DiscreteAdaptTC *discr_tc;
-
-public:
-#ifdef MFEM_USE_MPI
-   TMOPDescentNewtonSolver(MPI_Comm comm, const IntegrationRule &irule)
-      : NewtonSolver(comm), parallel(true), ir(irule), discr_tc(NULL) { }
-#endif
-   TMOPDescentNewtonSolver(const IntegrationRule &irule)
-      : NewtonSolver(), parallel(false), ir(irule), discr_tc(NULL) { }
-
-   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
-
-   virtual void ProcessNewState(const Vector &x) const;
-};
-
-void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm,
-                       const TargetConstructor &tc, Mesh &pmesh,
-                       char *title, int position);
-#ifdef MFEM_USE_MPI
-void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm,
-                       const TargetConstructor &tc, ParMesh &pmesh,
-                       char *title, int position);
-#endif
-
-}
-
-#endif
@@ -276,7 +276,7 @@ inline T *Write(Memory<T> &mem, int size, bool on_dev = true)

 /** @brief Shortcut to Write(const Memory<T> &mem, int size, false) */
 template <typename T>
-inline T *HostWrite(Memory<T> &mem, int size)
+inline const T *HostWrite(const Memory<T> &mem, int size)
 {
   return mfem::Write(mem, size, false);
 }
@@ -11,7 +11,7 @@

 #include "../general/forall.hpp"

-#include <cstring> // std::memcpy, std::memcmp
+#include <cstring> // std::memcpy

 #include <list>
 #include <unordered_map>
@@ -57,15 +57,6 @@ MemoryClass operator*(MemoryClass mc1, MemoryClass mc2)
 }


-// Instantiate Memory<T>::PrintFlags for T = int and T = double.
-template void Memory<int>::PrintFlags() const;
-template void Memory<double>::PrintFlags() const;
-
-// Instantiate Memory<T>::CompareHostAndDevice for T = int and T = double.
-template int Memory<int>::CompareHostAndDevice(int size) const;
-template int Memory<double>::CompareHostAndDevice(int size) const;
-
-
 namespace internal
 {

@@ -704,19 +695,6 @@ void MemoryManager::CopyFromHost_(void *dest_h_ptr, const void *src_h_ptr,
                ~(dest_on_host ? Mem::VALID_DEVICE : Mem::VALID_HOST);
 }

-int MemoryManager::CompareHostAndDevice_(void *h_ptr, size_t size,
-                                         unsigned flags)
-{
-   void *d_ptr = (flags & Mem::ALIAS) ?
-                 mm.GetAliasDevicePtr(h_ptr, size, false) :
-                 mm.GetDevicePtr(h_ptr, size, false);
-   char *h_buf = new char[size];
-   CuMemcpyDtoH(h_buf, d_ptr, size);
-   int res = std::memcmp(h_ptr, h_buf, size);
-   delete [] h_buf;
-   return res;
-}
-

 void MemoryPrintFlags(unsigned flags)
 {
@@ -377,16 +377,6 @@ public:
   /// Copy @a size entries from @a *this to the host pointer @a dest.
   /** The given @a size should not exceed the Capacity() of @a *this. */
   inline void CopyToHost(T *dest, int size) const;
-
-   /// Print the internal flags.
-   /** This method can be useful for debugging. It is explicitly instantiated
-       for Memory<T> with T = int and T = double. */
-   inline void PrintFlags() const;
-
-   /// If both the host and the device data are valid, compare their contents.
-   /** This method can be useful for debugging. It is explicitly instantiated
-       for Memory<T> with T = int and T = double. */
-   inline int CompareHostAndDevice(int size) const;
 };


@@ -457,11 +447,6 @@ private:
   static void CopyFromHost_(void *dest_h_ptr, const void *src_h_ptr,
                             std::size_t size, unsigned &dest_flags);

-   // Compare the contents of the host and the device memory - useful for
-   // debugging.
-   static int CompareHostAndDevice_(void *h_ptr, size_t size, unsigned flags);
-
-
   /// Adds an address in the map
   void *Insert(void *ptr, const std::size_t bytes);

@@ -742,24 +727,10 @@ inline void Memory<T>::CopyToHost(T *dest, int size) const


 /** @brief Print the state of a Memory object based on its internal flags.
-    Useful in a debugger. See also Memory<T>::PrintFlags(). */
+    Useful in a debugger. */
 extern void MemoryPrintFlags(unsigned flags);


-template <typename T>
-inline void Memory<T>::PrintFlags() const
-{
-   MemoryPrintFlags(flags);
-}
-
-template <typename T>
-inline int Memory<T>::CompareHostAndDevice(int size) const
-{
-   if (!(flags & VALID_HOST) || !(flags & VALID_DEVICE)) { return 0; }
-   return MemoryManager::CompareHostAndDevice_(h_ptr, size*sizeof(T), flags);
-}
-
-
 /// The (single) global memory manager object
 extern MemoryManager mm;

@@ -4456,8 +4456,7 @@ const
   int n = SizeI(), ne = SizeK();
   const int *I = elem_dof.GetI(), *J = elem_dof.GetJ(), *dofs;
   const double *d_col = tdata;
-   double *yp = y.HostReadWrite();
-   double x_col;
+   double *yp = y, x_col;
   const double *xp = x;
   // the '4' here can be tuned for given platform and compiler
   if (n <= 4)
@@ -185,7 +185,7 @@ Vector * HypreParVector::GlobalVector() const

 HypreParVector& HypreParVector::operator=(double d)
 {
-   Vector::operator=(d);
+   hypre_ParVectorSetConstantValues(x,d);
   return *this;
 }

@@ -198,7 +198,10 @@ HypreParVector& HypreParVector::operator=(const HypreParVector &y)
   }
 #endif

-   Vector::operator=(y);
+   for (int i = 0; i < size; i++)
+   {
+      data[i] = y.data[i];
+   }
   return *this;
 }

@@ -1021,8 +1024,6 @@ HypreParMatrix * HypreParMatrix::Transpose() const
 HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
                               double a, double b)
 {
-   x.HostRead();
-   (b == 0.0) ? y.HostWrite() : y.HostReadWrite();
   return hypre_ParCSRMatrixMatvec(a, A, x, b, y);
 }

@@ -1034,7 +1035,7 @@ void HypreParMatrix::Mult(double a, const Vector &x, double b, Vector &y) const
               << ", expected size = " << Height());

   auto x_data = x.HostRead();
-   auto y_data = (b == 0.0) ? y.HostWrite() : y.HostReadWrite();
+   auto y_data = y.HostWrite();
   if (X == NULL)
   {
      X = new HypreParVector(A->comm,
@@ -1066,7 +1067,7 @@ void HypreParMatrix::MultTranspose(double a, const Vector &x,
   // Note: x has the dimensions of Y (height), and
   //       y has the dimensions of X (width)
   auto x_data = x.HostRead();
-   auto y_data = (b == 0.0) ? y.HostWrite() : y.HostReadWrite();
+   auto y_data = y.HostWrite();
   if (X == NULL)
   {
      X = new HypreParVector(A->comm,
@@ -2008,12 +2009,10 @@ void HypreSmoother::Mult(const HypreParVector &b, HypreParVector &x) const
      return;
   }

-   b.HostRead();
   if (!iterative_mode)
   {
      if (type == 0 && relax_times == 1)
      {
-         x.HostWrite();
         HYPRE_ParCSRDiagScale(NULL, *A, b, x);
         if (relax_weight != 1.0)
         {
@@ -2023,7 +2022,6 @@ void HypreSmoother::Mult(const HypreParVector &b, HypreParVector &x) const
      }
      x = 0.0;
   }
-   x.HostReadWrite();

   if (V == NULL)
   {
@@ -2073,25 +2071,21 @@ void HypreSmoother::Mult(const Vector &b, Vector &x) const
      mfem_error("HypreSmoother::Mult (...) : HypreParMatrix A is missing");
      return;
   }
-
-   auto b_data = b.HostRead();
-   auto x_data = iterative_mode ? x.HostReadWrite() : x.HostWrite();
-
   if (B == NULL)
   {
      B = new HypreParVector(A->GetComm(),
                             A -> GetGlobalNumRows(),
-                             const_cast<double*>(b_data),
+                             b.GetData(),
                             A -> GetRowStarts());
      X = new HypreParVector(A->GetComm(),
                             A -> GetGlobalNumCols(),
-                             x_data,
+                             x.GetData(),
                             A -> GetColStarts());
   }
   else
   {
-      B -> SetData(const_cast<double*>(b_data));
-      X -> SetData(x_data);
+      B -> SetData(b.GetData());
+      X -> SetData(x.GetData());
   }

   Mult(*B, *X);
@@ -2693,7 +2687,7 @@ void HypreEuclid::SetOperator(const Operator &op)
   if (A)
   {
      MPI_Comm comm;
-      HYPRE_ParCSRMatrixGetComm(*new_A, &comm);
+      HYPRE_ParCSRMatrixGetComm(*A, &comm);
      ResetEuclidPrecond(comm);
   }

@@ -2309,6 +2309,17 @@ void PetscLinearSolver::SetOperator(const Operator &op)
                              (dynamic_cast<const PetscParMatrix *>(&op));
   const Operator       *oA = dynamic_cast<const Operator *>(&op);

+   // Preserve Pmat if already set
+   KSP ksp = (KSP)obj;
+   Mat P = NULL;
+   PetscBool pmat;
+   ierr = KSPGetOperatorsSet(ksp,NULL,&pmat); PCHKERRQ(ksp,ierr);
+   if (pmat)
+   {
+      ierr = KSPGetOperators(ksp,NULL,&P); PCHKERRQ(ksp,ierr);
+      ierr = PetscObjectReference((PetscObject)P); PCHKERRQ(ksp,ierr);
+   }
+
   // update base classes: Operator, Solver, PetscLinearSolver
   bool delete_pA = false;
   if (!pA)
@@ -2331,7 +2342,6 @@ void PetscLinearSolver::SetOperator(const Operator &op)
   MFEM_VERIFY(pA, "Unsupported operation!");

   // Set operators into PETSc KSP
-   KSP ksp = (KSP)obj;
   Mat A = pA->A;
   if (operatorset)
   {
@@ -2352,7 +2362,15 @@ void PetscLinearSolver::SetOperator(const Operator &op)
         wrap = false;
      }
   }
-   ierr = KSPSetOperators(ksp,A,A); PCHKERRQ(ksp,ierr);
+   if (P)
+   {
+      ierr = KSPSetOperators(ksp,A,P); PCHKERRQ(ksp,ierr);
+      ierr = MatDestroy(&P); PCHKERRQ(ksp,ierr);
+   }
+   else
+   {
+      ierr = KSPSetOperators(ksp,A,A); PCHKERRQ(ksp,ierr);
+   }

   // Update PetscSolver
   operatorset = true;
@@ -3211,27 +3229,26 @@ PetscFieldSplitSolver::PetscFieldSplitSolver(MPI_Comm comm, Operator &op,
   : PetscPreconditioner(comm,op,prefix)
 {
   PC pc = (PC)obj;
-   ierr = PCSetType(pc,PCFIELDSPLIT); PCHKERRQ(pc,ierr);

   Mat pA;
   ierr = PCGetOperators(pc,&pA,NULL); PCHKERRQ(pc,ierr);

   // Check if pA is of type MATNEST
+   // (this requirement can be removed when we can pass fields).
   PetscBool isnest;
   ierr = PetscObjectTypeCompare((PetscObject)pA,MATNEST,&isnest);
+   PCHKERRQ(pA,ierr);
+   MFEM_VERIFY(isnest,
+               "PetscFieldSplitSolver needs the matrix in nested format.");

-   PetscInt nr = 0;
-   IS  *isrow = NULL;
-   if (isnest) // we now the fields
-   {
-      ierr = MatNestGetSize(pA,&nr,NULL); PCHKERRQ(pc,ierr);
-      ierr = PetscCalloc1(nr,&isrow); CCHKERRQ(PETSC_COMM_SELF,ierr);
-      ierr = MatNestGetISs(pA,isrow,NULL); PCHKERRQ(pc,ierr);
-   }
+   PetscInt nr;
+   IS  *isrow;
+   ierr = PCSetType(pc,PCFIELDSPLIT); PCHKERRQ(pc,ierr);
+   ierr = MatNestGetSize(pA,&nr,NULL); PCHKERRQ(pc,ierr);
+   ierr = PetscCalloc1(nr,&isrow); CCHKERRQ(PETSC_COMM_SELF,ierr);
+   ierr = MatNestGetISs(pA,isrow,NULL); PCHKERRQ(pc,ierr);

   // We need to customize here, before setting the index sets.
-   // This is because PCFieldSplitSetType customizes the function
-   // pointers. SubSolver options will be processed during PCApply
   Customize();

   for (PetscInt i=0; i<nr; i++)
@@ -3599,7 +3616,6 @@ void PetscODESolver::Run(Vector &x, double &t, double &dt, double t_final)
 }  // namespace mfem

 #include "petsc/private/petscimpl.h"
-#include "petsc/private/matimpl.h"

 // auxiliary functions
 static PetscErrorCode __mfem_ts_monitor(TS ts, PetscInt it, PetscReal t, Vec x,
@@ -3755,37 +3771,17 @@ static PetscErrorCode __mfem_ts_ijacobian(TS ts, PetscReal t, Vec x,
      pA->EliminateRowsCols(bchandler->GetTDofs(),dummy,dummy);
   }

-   // Get nonzerostate
-   PetscObjectState nonzerostate;
-   ierr = MatGetNonzeroState(P,&nonzerostate); CHKERRQ(ierr);
-
   // Avoid unneeded copy of the matrix by hacking
   Mat B;
   B = pA->ReleaseMat(false);
   ierr = MatHeaderReplace(P,&B); CHKERRQ(ierr);
   if (delete_pA) { delete pA; }
-
-   // Matrix-free case
   if (A && A != P)
   {
      ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
      ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
   }

-   // When using MATNEST and PCFIELDSPLIT, the second setup of the
-   // preconditioner fails because MatCreateSubMatrix_Nest does not
-   // actually return a matrix. Instead, for efficiency reasons,
-   // it returns a reference to the submatrix. The second time it
-   // is called, MAT_REUSE_MATRIX is used and MatCreateSubMatrix_Nest
-   // aborts since the two submatrices are actually different.
-   // We circumvent this issue by incrementing the nonzero state
-   // (i.e. PETSc thinks the operator sparsity pattern has changed)
-   // This does not impact performances in the case of MATNEST
-   PetscBool isnest;
-   ierr = PetscObjectTypeCompare((PetscObject)P,MATNEST,&isnest);
-   CHKERRQ(ierr);
-   if (isnest) { P->nonzerostate = nonzerostate + 1; }
-
   // Jacobian reusage
   ierr = PetscObjectStateGet((PetscObject)P,&ts_ctx->cached_ijacstate);
   CHKERRQ(ierr);
@@ -3942,7 +3938,6 @@ static PetscErrorCode __mfem_ts_computesplits(TS ts,PetscReal t,Vec x,Vec xp,
      ierr = MatAXPY(*pJxp,-1.0,*pJx,SAME_NONZERO_PATTERN); PCHKERRQ(ts,ierr);
   }

-   // Matrix-free cases
   if (Ax && Ax != Jx)
   {
      ierr = MatAssemblyBegin(Ax,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
@@ -4029,31 +4024,12 @@ static PetscErrorCode __mfem_ts_rhsjacobian(TS ts, PetscReal t, Vec x,
      pA->EliminateRowsCols(bchandler->GetTDofs(),dummy,dummy);
   }

-   // Get nonzerostate
-   PetscObjectState nonzerostate;
-   ierr = MatGetNonzeroState(P,&nonzerostate); CHKERRQ(ierr);
-
   // Avoid unneeded copy of the matrix by hacking
   Mat B;
   B = pA->ReleaseMat(false);
   ierr = MatHeaderReplace(P,&B); CHKERRQ(ierr);
   if (delete_pA) { delete pA; }

-   // When using MATNEST and PCFIELDSPLIT, the second setup of the
-   // preconditioner fails because MatCreateSubMatrix_Nest does not
-   // actually return a matrix. Instead, for efficiency reasons,
-   // it returns a reference to the submatrix. The second time it
-   // is called, MAT_REUSE_MATRIX is used and MatCreateSubMatrix_Nest
-   // aborts since the two submatrices are actually different.
-   // We circumvent this issue by incrementing the nonzero state
-   // (i.e. PETSc thinks the operator sparsity pattern has changed)
-   // This does not impact performances in the case of MATNEST
-   PetscBool isnest;
-   ierr = PetscObjectTypeCompare((PetscObject)P,MATNEST,&isnest);
-   CHKERRQ(ierr);
-   if (isnest) { P->nonzerostate = nonzerostate + 1; }
-
-   // Matrix-free case
   if (A && A != P)
   {
      ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
@@ -4157,30 +4133,10 @@ static PetscErrorCode __mfem_snes_jacobian(SNES snes, Vec x, Mat A, Mat P,
      pA->EliminateRowsCols(bchandler->GetTDofs(),dummy,dummy);
   }

-   // Get nonzerostate
-   PetscObjectState nonzerostate;
-   ierr = MatGetNonzeroState(P,&nonzerostate); CHKERRQ(ierr);
-
   // Avoid unneeded copy of the matrix by hacking
   Mat B = pA->ReleaseMat(false);
   ierr = MatHeaderReplace(P,&B); CHKERRQ(ierr);
   if (delete_pA) { delete pA; }
-
-   // When using MATNEST and PCFIELDSPLIT, the second setup of the
-   // preconditioner fails because MatCreateSubMatrix_Nest does not
-   // actually return a matrix. Instead, for efficiency reasons,
-   // it returns a reference to the submatrix. The second time it
-   // is called, MAT_REUSE_MATRIX is used and MatCreateSubMatrix_Nest
-   // aborts since the two submatrices are actually different.
-   // We circumvent this issue by incrementing the nonzero state
-   // (i.e. PETSc thinks the operator sparsity pattern has changed)
-   // This does not impact performances in the case of MATNEST
-   PetscBool isnest;
-   ierr = PetscObjectTypeCompare((PetscObject)P,MATNEST,&isnest);
-   CHKERRQ(ierr);
-   if (isnest) { P->nonzerostate = nonzerostate + 1; }
-
-   // Matrix-free case
   if (A && A != P)
   {
      ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
@@ -618,8 +618,6 @@ public:
                     const std::string &prefix = std::string());
   virtual ~PetscLinearSolver();

-   /// Sets the operator to be used for mat-vec operations and
-   /// for the construction of the preconditioner
   virtual void SetOperator(const Operator &op);

   /// Allows to prescribe a different operator (@a pop) to construct
@@ -627,7 +625,6 @@ public:
   void SetOperator(const Operator &op, const Operator &pop);

   /// Sets the solver to perform preconditioning
-   /// preserves the linear operator for the mat-vec
   void SetPreconditioner(Solver &precond);

   /// Application of the solver.
@@ -1311,8 +1311,6 @@ void NewtonSolver::Mult(const Vector &b, Vector &x) const
      }
      add(x, -c_scale, c, x);

-      ProcessNewState(x);
-
      oper->Mult(x, r);
      if (have_b)
      {
@@ -283,10 +283,6 @@ public:
       value of 0 indicates a failure, interrupting the Newton iteration. */
   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const
   { return 1.0; }
-
-   /** @brief This method can be overloaded in derived classes to perform
-       computations that need knowledge of the newest Newton state. */
-   virtual void ProcessNewState(const Vector &x) const { }
 };

 /** Adaptive restarted GMRES.
@@ -2640,11 +2640,9 @@ SparseMatrix &SparseMatrix::operator=(double a)
 {
   if (Rows == NULL)
   {
-      const int nnz = J.Capacity();
-      double *h_A = HostWrite(A, nnz);
-      for (int i = 0; i < nnz; i++)
+      for (int i = 0, nnz = I[height]; i < nnz; i++)
      {
-         h_A[i] = a;
+         A[i] = a;
      }
   }
   else
@@ -152,54 +152,6 @@ public:
   /// Return the element data, i.e. the array #A, const version.
   inline const double *GetData() const { return A; }

-   // Memory access methods for the #I array.
-   Memory<int> &GetMemoryI() { return I; }
-   const Memory<int> &GetMemoryI() const { return I; }
-   const int *ReadI(bool on_dev = true) const
-   { return mfem::Read(I, Height()+1, on_dev); }
-   int *WriteI(bool on_dev = true)
-   { return mfem::Write(I, Height()+1, on_dev); }
-   int *ReadWriteI(bool on_dev = true)
-   { return mfem::ReadWrite(I, Height()+1, on_dev); }
-   const int *HostReadI() const
-   { return mfem::Read(I, Height()+1, false); }
-   int *HostWriteI()
-   { return mfem::Write(I, Height()+1, false); }
-   int *HostReadWriteI()
-   { return mfem::ReadWrite(I, Height()+1, false); }
-
-   // Memory access methods for the #J array.
-   Memory<int> &GetMemoryJ() { return J; }
-   const Memory<int> &GetMemoryJ() const { return J; }
-   const int *ReadJ(bool on_dev = true) const
-   { return mfem::Read(J, J.Capacity(), on_dev); }
-   int *WriteJ(bool on_dev = true)
-   { return mfem::Write(J, J.Capacity(), on_dev); }
-   int *ReadWriteJ(bool on_dev = true)
-   { return mfem::ReadWrite(J, J.Capacity(), on_dev); }
-   const int *HostReadJ() const
-   { return mfem::Read(J, J.Capacity(), false); }
-   int *HostWriteJ()
-   { return mfem::Write(J, J.Capacity(), false); }
-   int *HostReadWriteJ()
-   { return mfem::ReadWrite(J, J.Capacity(), false); }
-
-   // Memory access methods for the #A array.
-   Memory<double> &GetMemoryData() { return A; }
-   const Memory<double> &GetMemoryData() const { return A; }
-   const double *ReadData(bool on_dev = true) const
-   { return mfem::Read(A, A.Capacity(), on_dev); }
-   double *WriteData(bool on_dev = true)
-   { return mfem::Write(A, A.Capacity(), on_dev); }
-   double *ReadWriteData(bool on_dev = true)
-   { return mfem::ReadWrite(A, A.Capacity(), on_dev); }
-   const double *HostReadData() const
-   { return mfem::Read(A, A.Capacity(), false); }
-   double *HostWriteData()
-   { return mfem::Write(A, A.Capacity(), false); }
-   double *HostReadWriteData()
-   { return mfem::ReadWrite(A, A.Capacity(), false); }
-
   /// Returns the number of elements in row @a i.
   int RowSize(const int i) const;

@@ -134,11 +134,6 @@ SuperLURowLocMatrix::SuperLURowLocMatrix( const HypreParMatrix & hypParMat )
   // hypre_CSRMatrix.
   hypre_CSRMatrix * csr_op = hypre_MergeDiagAndOffd(parcsr_op);
   hypre_CSRMatrixSetDataOwner(csr_op,0);
-#if MFEM_HYPRE_VERSION >= 21600
-   MFEM_VERIFY(csr_op->num_rows < INT_MAX,"SuperLU: number of local rows "
-               "is too large to store as an integer.");
-   hypre_CSRMatrixBigJtoJ(csr_op);
-#endif

   int m         = parcsr_op->global_num_rows;
   int n         = parcsr_op->global_num_cols;
@@ -834,10 +834,9 @@ double Vector::Sum() const
 {
   double sum = 0.0;

-   const double *h_data = this->HostRead();
   for (int i = 0; i < size; i++)
   {
-      sum += h_data[i];
+      sum += data[i];
   }

   return sum;
@@ -140,13 +140,6 @@ public:
       @sa NewDataAndSize(). */
   inline void NewMemoryAndSize(const Memory<double> &mem, int s, bool own_mem);

-   /// Reset the Vector to be a reference to a sub-vector of @a base.
-   inline void MakeRef(Vector &base, int offset, int size);
-
-   /** @brief Reset the Vector to be a reference to a sub-vector of @a base
-       without changing its current size. */
-   inline void MakeRef(Vector &base, int offset);
-
   /// Set the Vector data (host pointer) ownership flag.
   void MakeDataOwner() const { data.SetHostPtrOwner(true); }

@@ -462,19 +455,6 @@ inline void Vector::NewMemoryAndSize(const Memory<double> &mem, int s,
   if (!own_mem) { data.ClearOwnerFlags(); }
 }

-inline void Vector::MakeRef(Vector &base, int offset, int s)
-{
-   data.Delete();
-   size = s;
-   data.MakeAlias(base.GetMemory(), offset, s);
-}
-
-inline void Vector::MakeRef(Vector &base, int offset)
-{
-   data.Delete();
-   data.MakeAlias(base.GetMemory(), offset, size);
-}
-
 inline void Vector::Destroy()
 {
   const bool use_dev = data.UseDevice();
@@ -25,6 +25,7 @@
 #include <cstring>
 #include <ctime>
 #include <functional>
+#include <set>

 // Include the METIS header, if using version 5. If using METIS 4, the needed
 // declarations are inlined below, i.e. no header is needed.
@@ -430,8 +431,6 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
   else
   {
      const FiniteElement *bdr_el = Nodes->FESpace()->GetBE(i);
-      Nodes->HostRead();
-      const GridFunction &nodes = *Nodes;
      if (bdr_el)
      {
         Array<int> vdofs;
@@ -442,7 +441,7 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
         {
            for (int j = 0; j < n; j++)
            {
-               pm(k,j) = nodes(vdofs[n*k+j]);
+               pm(k,j) = (*Nodes)(vdofs[n*k+j]);
            }
         }
         ElTr->SetFE(bdr_el);
@@ -494,8 +493,6 @@ void Mesh::GetFaceTransformation(int FaceNo, IsoparametricTransformation *FTr)
   else // curved mesh
   {
      const FiniteElement *face_el = Nodes->FESpace()->GetFaceElement(FaceNo);
-      Nodes->HostRead();
-      const GridFunction &nodes = *Nodes;
      if (face_el)
      {
         Array<int> vdofs;
@@ -506,7 +503,7 @@ void Mesh::GetFaceTransformation(int FaceNo, IsoparametricTransformation *FTr)
         {
            for (int j = 0; j < n; j++)
            {
-               pm(i, j) = nodes(vdofs[n*i+j]);
+               pm(i, j) = (*Nodes)(vdofs[n*i+j]);
            }
         }
         FTr->SetFE(face_el);
@@ -7098,6 +7095,59 @@ bool Mesh::NonconformingDerefinement(Array<double> &elem_error,
   return true;
 }

+bool Mesh::GeneralDerefinement(Array<int>& derefs)
+{
+   MFEM_VERIFY(ncmesh, "Only supported for non-conforming meshes.");
+   MFEM_VERIFY(!NURBSext, "Derefinement of NURBS meshes is not supported. "
+               "Project the NURBS to Nodes first.");
+
+   DeleteLazyTables();
+
+   if (!derefs.Size()) { return false; }
+
+   const Table &dt = ncmesh->GetDerefinementTable();
+
+   // map elements to rows in table
+   Table dtT;
+   Transpose(dt, dtT);
+
+   // Find the set of rows in the refinement table which correspond to
+   // the specified elements to derefine
+   std::set<int> set_rows;
+   for (int i = 0; i < derefs.Size(); i++) {
+      int el = derefs[i];
+      const int *row = dtT.GetRow(el);
+      set_rows.insert(row[0]);
+   }
+
+   Array<int> deref_rows;
+   std::set<int>::iterator it;
+   for (it = set_rows.begin(); it != set_rows.end(); ++it) {
+      deref_rows.Append(*it);
+   }
+   
+   ncmesh->Derefine(deref_rows);
+
+   Mesh* mesh2 = new Mesh(*ncmesh);
+   ncmesh->OnMeshUpdated(mesh2);
+
+   Swap(*mesh2, false);
+   delete mesh2;
+
+   GenerateNCFaceInfo();
+
+   last_operation = Mesh::DEREFINE;
+   sequence++;
+
+   if (Nodes) // update/interpolate mesh curvature
+   {
+      Nodes->FESpace()->Update();
+      Nodes->Update();
+   }
+
+   return true;
+}
+
 bool Mesh::DerefineByError(Array<double> &elem_error, double threshold,
                           int nc_limit, int op)
 {
@@ -1103,6 +1103,9 @@ public:
   bool RefineByError(const Vector &elem_error, double threshold,
                      int nonconforming = -1, int nc_limit = 0);

+   /** Derefine specified elements (and associated children as required) */
+   bool GeneralDerefinement(Array<int> &el);
+
   /** Derefine the mesh based on an error measure associated with each
       element. A derefinement is performed if the sum of errors of its fine
       elements is smaller than 'threshold'. If 'nc_limit' > 0, derefinements
@@ -4307,20 +4307,6 @@ int NCMesh::GetElementDepth(int i) const
   return depth;
 }

-int NCMesh::GetElementSizeReduction(int i) const
-{
-   int elem = leaf_elements[i];
-   int parent, reduction = 1;
-   while ((parent = elements[elem].parent) != -1)
-   {
-      if (elements[parent].ref_type & 1) { reduction *= 2; }
-      if (elements[parent].ref_type & 2) { reduction *= 2; }
-      if (elements[parent].ref_type & 4) { reduction *= 2; }
-      elem = parent;
-   }
-   return reduction;
-}
-
 void NCMesh::GetElementFacesAttributes(int i,
                                       Array<int> &faces,
                                       Array<int> &fattr) const
@@ -319,21 +319,13 @@ public:
   Geometry::Type GetFaceGeometry(int index) const
   { return Geometry::Type(face_geom[index]); }

-   /// Return the number of root elements.
-   int GetNumRootElements() { return root_state.Size(); }
-
   /// Return the distance of leaf 'i' from the root.
   int GetElementDepth(int i) const;

-   /** Return the size reduction compared to the root element (ignoring local
-       stretching and curvature). */
-   int GetElementSizeReduction(int i) const;
-
   /// Return the faces and face attributes of leaf element 'i'.
   void GetElementFacesAttributes(int i, Array<int> &faces,
                                  Array<int> &fattr) const;

-
   /// I/O: Print the "vertex_parents" section of the mesh file (ver. >= 1.1).
   void PrintVertexParents(std::ostream &out) const;

@@ -4171,7 +4171,6 @@ void ParMesh::PrintAsOne(std::ostream &out)
          "# SQUARE      = 3\n"
          "# TETRAHEDRON = 4\n"
          "# CUBE        = 5\n"
-          "# PRISM       = 6\n"
          "#\n";

      out << "\ndimension\n" << Dim;
@@ -4270,15 +4269,6 @@ void ParMesh::PrintAsOne(std::ostream &out)
   {
      switch (Dim)
      {
-         case 1:
-            for (i = 0; i < svert_lvert.Size(); i++)
-            {
-               ints.Append(Geometry::POINT);
-               ints.Append(svert_lvert[i]);
-               ne++;
-            }
-            break;
-
         case 2:
            for (i = 0; i < shared_edges.Size(); i++)
            {
@@ -41,7 +41,6 @@ VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
     H1FESpace_(NULL),
     HCurlFESpace_(NULL),
     HDivFESpace_(NULL),
-     L2FESpace_(NULL),
     divEpsGrad_(NULL),
     h1Mass_(NULL),
     h1SurfMass_(NULL),
@@ -50,24 +49,20 @@ VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
     hCurlHDiv_(NULL),
     weakDiv_(NULL),
     rhod_(NULL),
-     l2_vol_int_(NULL),
-     rt_surf_int_(NULL),
     grad_(NULL),
     phi_(NULL),
-     rho_src_(NULL),
     rho_(NULL),
-     sigma_src_(NULL),
+     sigma_(NULL),
     e_(NULL),
     d_(NULL),
-     p_src_(NULL),
-     oneCoef_(1.0),
+     p_(NULL),
     epsCoef_(&epsCoef),
     phiBCCoef_(NULL),
     rhoCoef_(NULL),
     pCoef_(NULL),
-     phi_bc_func_(phi_bc),
-     rho_src_func_(rho_src),
-     p_src_func_(p_src),
+     phi_bc_(phi_bc),
+     rho_src_(rho_src),
+     p_src_(p_src),
     point_charge_params_(point_charges),
     point_charges_(0)
 {
@@ -81,38 +76,34 @@ VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
   H1FESpace_    = new H1_ParFESpace(pmesh_,order,pmesh_->Dimension());
   HCurlFESpace_ = new ND_ParFESpace(pmesh_,order,pmesh_->Dimension());
   HDivFESpace_  = new RT_ParFESpace(pmesh_,order,pmesh_->Dimension());
-   L2FESpace_    = new L2_ParFESpace(pmesh_,order-1,pmesh_->Dimension());

   // Select surface attributes for Dirichlet BCs
   ess_bdr_.SetSize(pmesh.bdr_attributes.Max());
   ess_bdr_ = 0;   // Deselect all outer surfaces
   for (int i=0; i<dbcs_->Size(); i++)
   {
-      if ((*dbcs_)[i] <= ess_bdr_.Size())
-      {
-         ess_bdr_[(*dbcs_)[i]-1] = 1;
-      }
+      ess_bdr_[(*dbcs_)[i]-1] = 1;
   }

   // Setup various coefficients

   // Potential on outer surface
-   if ( phi_bc_func_ != NULL )
+   if ( phi_bc_ != NULL )
   {
-      phiBCCoef_ = new FunctionCoefficient(*phi_bc_func_);
+      phiBCCoef_ = new FunctionCoefficient(*phi_bc_);
   }

   // Volume Charge Density
-   if ( rho_src_func_ != NULL )
+   if ( rho_src_ != NULL )
   {
-      rhoCoef_ = new FunctionCoefficient(rho_src_func_);
+      rhoCoef_ = new FunctionCoefficient(rho_src_);
   }

   // Polarization
-   if ( p_src_func_ != NULL )
+   if ( p_src_ != NULL )
   {
      pCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
-                                             p_src_func_);
+                                             p_src_);
   }

   // Bilinear Forms
@@ -125,23 +116,15 @@ VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
   hCurlHDivEps_ = new ParMixedBilinearForm(HCurlFESpace_,HDivFESpace_);
   hCurlHDivEps_->AddDomainIntegrator(new VectorFEMassIntegrator(*epsCoef_));

-   rhod_   = new ParLinearForm(H1FESpace_);
+   rhod_ = new ParLinearForm(H1FESpace_);

-   l2_vol_int_ = new ParLinearForm(L2FESpace_);
-   l2_vol_int_->AddDomainIntegrator(new DomainLFIntegrator(oneCoef_));
-
-   rt_surf_int_ = new ParLinearForm(HDivFESpace_);
-   rt_surf_int_->AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator);
-
-   // Discrete derivative operator
+   // Discrete Grad operator
   grad_ = new ParDiscreteGradOperator(H1FESpace_, HCurlFESpace_);
-   div_  = new ParDiscreteDivOperator(HDivFESpace_, L2FESpace_);

   // Build grid functions
   phi_  = new ParGridFunction(H1FESpace_);
   d_    = new ParGridFunction(HDivFESpace_);
   e_    = new ParGridFunction(HCurlFESpace_);
-   rho_  = new ParGridFunction(L2FESpace_);

   if ( point_charge_params_.Size() > 0 )
   {
@@ -166,28 +149,29 @@ VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
      }
   }

-   if ( rho_src_func_ )
+   if ( rho_src_ )
   {
-      rho_src_ = new ParGridFunction(H1FESpace_);
+      rho_ = new ParGridFunction(H1FESpace_);

      h1Mass_ = new ParBilinearForm(H1FESpace_);
      h1Mass_->AddDomainIntegrator(new MassIntegrator);
   }

-   if ( p_src_func_ )
+   if ( p_src_ )
   {
-      p_src_ = new ParGridFunction(HCurlFESpace_);
+      p_ = new ParGridFunction(HCurlFESpace_);

      hCurlHDiv_ = new ParMixedBilinearForm(HCurlFESpace_, HDivFESpace_);
      hCurlHDiv_->AddDomainIntegrator(new VectorFEMassIntegrator);

      weakDiv_ = new ParMixedBilinearForm(HCurlFESpace_, H1FESpace_);
      weakDiv_->AddDomainIntegrator(new VectorFEWeakDivergenceIntegrator);
+
   }

   if ( nbcs_->Size() > 0 )
   {
-      sigma_src_ = new ParGridFunction(H1FESpace_);
+      sigma_ = new ParGridFunction(H1FESpace_);

      h1SurfMass_ = new ParBilinearForm(H1FESpace_);
      h1SurfMass_->AddBoundaryIntegrator(new MassIntegrator);
@@ -201,18 +185,14 @@ VoltaSolver::~VoltaSolver()
   delete pCoef_;

   delete phi_;
-   delete rho_src_;
   delete rho_;
   delete rhod_;
-   delete l2_vol_int_;
-   delete rt_surf_int_;
-   delete sigma_src_;
+   delete sigma_;
   delete d_;
   delete e_;
-   delete p_src_;
+   delete p_;

   delete grad_;
-   delete div_;

   delete divEpsGrad_;
   delete h1Mass_;
@@ -225,7 +205,6 @@ VoltaSolver::~VoltaSolver()
   delete H1FESpace_;
   delete HCurlFESpace_;
   delete HDivFESpace_;
-   delete L2FESpace_;

   for (unsigned int i=0; i<point_charges_.size(); i++)
   {
@@ -251,13 +230,11 @@ VoltaSolver::PrintSizes()
   HYPRE_Int size_h1 = H1FESpace_->GlobalTrueVSize();
   HYPRE_Int size_nd = HCurlFESpace_->GlobalTrueVSize();
   HYPRE_Int size_rt = HDivFESpace_->GlobalTrueVSize();
-   HYPRE_Int size_l2 = L2FESpace_->GlobalTrueVSize();
   if (myid_ == 0)
   {
      cout << "Number of H1      unknowns: " << size_h1 << endl;
      cout << "Number of H(Curl) unknowns: " << size_nd << endl;
      cout << "Number of H(Div)  unknowns: " << size_rt << endl;
-      cout << "Number of L2      unknowns: " << size_l2 << endl;
   }
 }

@@ -277,15 +254,9 @@ void VoltaSolver::Assemble()
   *rhod_ = 0.0;
   rhod_->Assemble();

-   l2_vol_int_->Assemble();
-   rt_surf_int_->Assemble();
-
   grad_->Assemble();
   grad_->Finalize();

-   div_->Assemble();
-   div_->Finalize();
-
   if ( h1Mass_ )
   {
      h1Mass_->Assemble();
@@ -321,19 +292,15 @@ VoltaSolver::Update()
   H1FESpace_->Update(false);
   HCurlFESpace_->Update(false);
   HDivFESpace_->Update(false);
-   L2FESpace_->Update(false);

   // Inform the grid functions that the space has changed.
   phi_->Update();
   rhod_->Update();
-   l2_vol_int_->Update();
-   rt_surf_int_->Update();
   d_->Update();
   e_->Update();
-   rho_->Update();
-   if ( rho_src_   ) { rho_src_->Update(); }
-   if ( sigma_src_ ) { sigma_src_->Update(); }
-   if ( p_src_     ) { p_src_->Update(); }
+   if ( rho_   ) { rho_->Update(); }
+   if ( sigma_ ) { sigma_->Update(); }
+   if ( p_     ) { p_->Update(); }

   // Inform the bilinear forms that the space has changed.
   divEpsGrad_->Update();
@@ -347,7 +314,6 @@ VoltaSolver::Update()

   // Inform the other objects that the space has changed.
   grad_->Update();
-   div_->Update();
 }

 void
@@ -373,46 +339,40 @@ VoltaSolver::Solve()
         {
            ConstantCoefficient voltage((*dbcv_)[i]);
            dbc_bdr_attr = 0;
-            if ((*dbcs_)[i] <= dbc_bdr_attr.Size())
-            {
-               dbc_bdr_attr[(*dbcs_)[i]-1] = 1;
-            }
+            dbc_bdr_attr[(*dbcs_)[i]-1] = 1;
            phi_->ProjectBdrCoefficient(voltage, dbc_bdr_attr);
         }
      }
   }

   // Initialize the volumetric charge density
-   if ( rho_src_ )
+   if ( rho_ )
   {
-      rho_src_->ProjectCoefficient(*rhoCoef_);
-      h1Mass_->AddMult(*rho_src_, *rhod_);
+      rho_->ProjectCoefficient(*rhoCoef_);
+      h1Mass_->AddMult(*rho_, *rhod_);
   }

   // Initialize the Polarization
-   if ( p_src_ )
+   if ( p_ )
   {
-      p_src_->ProjectCoefficient(*pCoef_);
-      weakDiv_->AddMult(*p_src_, *rhod_);
+      p_->ProjectCoefficient(*pCoef_);
+      weakDiv_->AddMult(*p_, *rhod_);
   }

   // Initialize the surface charge density
-   if ( sigma_src_ )
+   if ( sigma_ )
   {
-      *sigma_src_ = 0.0;
+      *sigma_ = 0.0;

      Array<int> nbc_bdr_attr(pmesh_->bdr_attributes.Max());
      for (int i=0; i<nbcs_->Size(); i++)
      {
         ConstantCoefficient sigma_coef((*nbcv_)[i]);
         nbc_bdr_attr = 0;
-         if ((*nbcs_)[i] <= nbc_bdr_attr.Size())
-         {
-            nbc_bdr_attr[(*nbcs_)[i]-1] = 1;
-         }
-         sigma_src_->ProjectBdrCoefficient(sigma_coef, nbc_bdr_attr);
+         nbc_bdr_attr[(*nbcs_)[i]-1] = 1;
+         sigma_->ProjectBdrCoefficient(sigma_coef, nbc_bdr_attr);
      }
-      h1SurfMass_->AddMult(*sigma_src_, *rhod_);
+      h1SurfMass_->AddMult(*sigma_, *rhod_);
   }

   // Determine the essential BC degrees of freedom
@@ -464,9 +424,9 @@ VoltaSolver::Solve()

   ParGridFunction ed(HDivFESpace_);
   hCurlHDivEps_->Mult(*e_, ed);
-   if ( p_src_ )
+   if ( p_ )
   {
-      hCurlHDiv_->AddMult(*p_src_, ed, -1.0);
+      hCurlHDiv_->AddMult(*p_, ed, -1.0);
   }

   HypreParMatrix MassHDiv;
@@ -485,27 +445,8 @@ VoltaSolver::Solve()

   hDivMass_->RecoverFEMSolution(D, ed, *d_);

-   // Compute charge density from rho = Div(D)
-   div_->Mult(*d_, *rho_);
-
   if (myid_ == 0) { cout << "done." << flush; }

-   {
-      // Compute total charge as volume integral of rho
-      double charge_rho = (*l2_vol_int_)(*rho_);
-
-      // Compute total charge as surface integral of D
-      double charge_D = (*rt_surf_int_)(*d_);
-
-      if (myid_ == 0)
-      {
-         cout << endl << "Total charge: \n"
-              << "   Volume integral of charge density:   " << charge_rho
-              << "\n   Surface integral of dielectric flux: " << charge_D
-              << endl << flush;
-      }
-   }
-
   if (myid_ == 0) { cout << "Solver done. " << endl; }
 }

@@ -539,10 +480,9 @@ VoltaSolver::RegisterVisItFields(VisItDataCollection & visit_dc)
   visit_dc.RegisterField("Phi", phi_);
   visit_dc.RegisterField("D",     d_);
   visit_dc.RegisterField("E",     e_);
-   visit_dc.RegisterField("Rho", rho_);
-   if ( rho_src_   ) { visit_dc.RegisterField("Rho Source",     rho_src_); }
-   if ( p_src_     ) { visit_dc.RegisterField("P Source",         p_src_); }
-   if ( sigma_src_ ) { visit_dc.RegisterField("Sigma Source", sigma_src_); }
+   if ( rho_   ) { visit_dc.RegisterField("Rho",     rho_); }
+   if ( p_     ) { visit_dc.RegisterField("P",         p_); }
+   if ( sigma_ ) { visit_dc.RegisterField("Sigma", sigma_); }
 }

 void
@@ -575,23 +515,20 @@ VoltaSolver::InitializeGLVis()
   socks_["E"] = new socketstream;
   socks_["E"]->precision(8);

-   socks_["Rho"] = new socketstream;
-   socks_["Rho"]->precision(8);
-
-   if ( rho_src_ )
+   if ( rho_)
   {
-      socks_["RhoSrc"] = new socketstream;
-      socks_["RhoSrc"]->precision(8);
+      socks_["Rho"] = new socketstream;
+      socks_["Rho"]->precision(8);
   }
-   if ( p_src_ )
+   if ( p_)
   {
-      socks_["PSrc"] = new socketstream;
-      socks_["PSrc"]->precision(8);
+      socks_["P"] = new socketstream;
+      socks_["P"]->precision(8);
   }
-   if ( sigma_src_ )
+   if ( sigma_)
   {
-      socks_["SigmaSrc"] = new socketstream;
-      socks_["SigmaSrc"]->precision(8);
+      socks_["Sigma"] = new socketstream;
+      socks_["Sigma"]->precision(8);
   }
 }

@@ -611,36 +548,31 @@ VoltaSolver::DisplayToGLVis()
                  *phi_, "Electric Potential (Phi)", Wx, Wy, Ww, Wh);
   Wx += offx;

-   VisualizeField(*socks_["E"], vishost, visport,
-                  *e_, "Electric Field (E)", Wx, Wy, Ww, Wh);
-   Wx += offx;
-
   VisualizeField(*socks_["D"], vishost, visport,
                  *d_, "Electric Displacement (D)", Wx, Wy, Ww, Wh);
   Wx += offx;

-   VisualizeField(*socks_["Rho"], vishost, visport,
-                  *rho_, "Charge Density", Wx, Wy, Ww, Wh);
+   VisualizeField(*socks_["E"], vishost, visport,
+                  *e_, "Electric Field (E)", Wx, Wy, Ww, Wh);
+
   Wx = 0; Wy += offy; // next line

-   if ( rho_src_ )
+   if ( rho_ )
   {
-      VisualizeField(*socks_["RhoSrc"], vishost, visport,
-                     *rho_src_, "Charge Density Source (Rho)", Wx, Wy, Ww, Wh);
+      VisualizeField(*socks_["Rho"], vishost, visport,
+                     *rho_, "Charge Density (Rho)", Wx, Wy, Ww, Wh);
      Wx += offx;
   }
-   if ( p_src_ )
+   if ( p_ )
   {
-      VisualizeField(*socks_["PSrc"], vishost, visport,
-                     *p_src_, "Electric Polarization Source (P)",
-                     Wx, Wy, Ww, Wh);
+      VisualizeField(*socks_["P"], vishost, visport,
+                     *p_, "Electric Polarization (P)", Wx, Wy, Ww, Wh);
      Wx += offx;
   }
-   if ( sigma_src_ )
+   if ( sigma_ )
   {
-      VisualizeField(*socks_["SigmaSrc"], vishost, visport,
-                     *sigma_src_, "Surface Charge Density Source (Sigma)",
-                     Wx, Wy, Ww, Wh);
+      VisualizeField(*socks_["Sigma"], vishost, visport,
+                     *sigma_, "Surface Charge Density (Sigma)", Wx, Wy, Ww, Wh);
      // Wx += offx; // not used
   }
   if (myid_ == 0) { cout << " done." << endl; }
@@ -26,9 +26,7 @@ namespace mfem
 using miniapps::H1_ParFESpace;
 using miniapps::ND_ParFESpace;
 using miniapps::RT_ParFESpace;
-using miniapps::L2_ParFESpace;
 using miniapps::ParDiscreteGradOperator;
-using miniapps::ParDiscreteDivOperator;

 namespace electromagnetics
 {
@@ -86,7 +84,6 @@ private:
   H1_ParFESpace * H1FESpace_;    // Continuous space for phi
   ND_ParFESpace * HCurlFESpace_; // Tangentially continuous space for E
   RT_ParFESpace * HDivFESpace_;  // Normally continuous space for D
-   L2_ParFESpace * L2FESpace_;    // Discontinuous space for rho

   ParBilinearForm * divEpsGrad_; // Laplacian operator
   ParBilinearForm * h1Mass_;     // For Volumetric Charge Density Source
@@ -97,32 +94,26 @@ private:
   ParMixedBilinearForm * hCurlHDiv_;    // For computing D from E and P
   ParMixedBilinearForm * weakDiv_;      // For computing the source term from P

-   ParLinearForm * rhod_; // Dual of Volumetric Charge Density Source
-
-   ParLinearForm * l2_vol_int_;  // Integral of L2 field
-   ParLinearForm * rt_surf_int_; // Integral of H(Div) field over boundary
+   ParLinearForm * rhod_; // Dual of Volumetric Charge Density

   ParDiscreteGradOperator * grad_; // For Computing E from phi
-   ParDiscreteDivOperator  * div_;  // For Computing rho from D

-   ParGridFunction * phi_;       // Electric Scalar Potential
-   ParGridFunction * rho_src_;   // Volumetric Charge Density Source
-   ParGridFunction * rho_;       // Volumetric Charge Density (Div(D))
-   ParGridFunction * sigma_src_; // Surface Charge Density Source
-   ParGridFunction * e_;         // Electric Field
-   ParGridFunction * d_;         // Electric Flux Density (aka Dielectric Flux)
-   ParGridFunction * p_src_;     // Polarization Field Source
+   ParGridFunction * phi_; // Electric Scalar Potential
+   ParGridFunction * rho_; // Volumetric Charge Density
+   ParGridFunction * sigma_; // Surface Charge Density
+   ParGridFunction * e_; // Electric Field
+   ParGridFunction * d_; // Electric Flux Density (aka Dielectric Flux)
+   ParGridFunction * p_; // Polarization Field

-   ConstantCoefficient oneCoef_;   // Coefficient equal to 1
   Coefficient       * epsCoef_;   // Dielectric Permittivity Coefficient
   Coefficient       * phiBCCoef_; // Scalar Potential Boundary Condition
   Coefficient       * rhoCoef_;   // Charge Density Coefficient
   VectorCoefficient * pCoef_;     // Polarization Vector Field Coefficient

   // Source functions
-   double (*phi_bc_func_ )(const Vector&);          // Scalar Potential BC
-   double (*rho_src_func_)(const Vector&);          // Volumetric Charge Density
-   void   (*p_src_func_  )(const Vector&, Vector&); // Polarization Field
+   double (*phi_bc_ )(const Vector&);          // Scalar Potential BC
+   double (*rho_src_)(const Vector&);          // Volumetric Charge Density
+   void   (*p_src_  )(const Vector&, Vector&); // Polarization Field

   const Vector & point_charge_params_;

@@ -1,418 +0,0 @@
-MFEM mesh v1.1
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-# PRISM       = 6
-#
-
-dimension
-2
-
-elements
-28
-1 3 0 9 21 12
-1 3 12 21 11 7
-1 3 21 22 37 25
-1 3 25 37 24 11
-1 3 37 23 8 24
-1 3 22 10 23 37
-1 3 9 4 10 21
-1 3 4 13 26 10
-1 3 13 1 14 26
-1 3 26 14 5 15
-1 3 38 28 15 29
-1 3 27 26 28 38
-1 3 10 27 38 23
-1 3 23 38 29 8
-1 3 8 29 39 33
-1 3 33 39 32 18
-1 3 39 31 30 32
-1 3 29 15 31 39
-1 3 15 5 16 30
-1 3 30 16 3 17
-1 3 18 30 17 6
-1 3 34 18 6 19
-1 3 40 33 18 36
-1 3 24 8 33 40
-1 3 11 24 40 35
-1 3 35 40 36 34
-1 3 7 11 34 20
-1 3 20 34 19 2
-
-boundary
-16
-2 1 0 9
-1 1 12 0
-1 1 7 12
-2 1 9 4
-2 1 4 13
-2 1 13 1
-1 1 1 14
-1 1 14 5
-1 1 5 16
-1 1 16 3
-2 1 3 17
-2 1 17 6
-2 1 6 19
-1 1 20 7
-2 1 19 2
-1 1 2 20
-
-vertex_parents
-37
-4 0 1
-5 1 3
-6 2 3
-7 0 2
-8 4 6
-9 0 4
-10 4 8
-11 7 8
-12 0 7
-13 1 4
-14 1 5
-15 5 8
-16 3 5
-17 3 6
-18 6 8
-19 2 6
-20 2 7
-21 9 11
-22 10 21
-23 8 10
-24 8 11
-25 11 21
-26 13 15
-27 10 26
-28 15 26
-29 8 15
-30 15 17
-31 15 30
-32 18 30
-33 8 18
-34 11 19
-35 11 34
-36 18 34
-37 22 24
-38 27 29
-39 29 32
-40 24 36
-
-coarse_elements
-9
-3 2 5 4 3
-3 0 6 28 1
-3 12 11 10 13
-3 7 8 9 30
-3 14 17 16 15
-3 32 18 19 20
-3 24 23 22 25
-3 26 34 21 27
-3 29 31 33 35
-
-vertices
-41
-
-nodes
-FiniteElementSpace
-FiniteElementCollection: H1_2D_P2
-VDim: 2
-Ordering: 0
-
-0
-1
-0
-1
-0.5
-1
-0.5
-0
-0.490917
-0.25
-0.481659
-0.229025
-0
-0.75
-1
-0.75002
-1
-0.75
-0.497731
-0.25
-0
-0.230224
-0.361824
-0.498798
-0.358947
-0.256374
-0.75627
-0.609156
-0.760931
-0.614878
-0.754468
-0.729327
-0.597212
-0.479497
-0.255854
-0.252575
-0.374701
-0.365809
-0.616315
-0.620005
-0.354453
-0.125
-0.241186
-0.113998
-0
-0.256374
-0.0996625
-0
-0.297495
-0.381624
-0.298118
-0.249986
-0.354661
-0.29368
-0.249387
-0.422671
-0.50313
-0.416618
-0.423212
-0.50141
-0.375
-0.498432
-0.361824
-0.625
-0.750551
-0.609156
-0.875
-1
-0.883899
-1
-0.859365
-0.760931
-0.66728
-0.757422
-0.680277
-0.604085
-0.680261
-0.760547
-0.611485
-0.542955
-0.563608
-0.556166
-0.630915
-0.564831
-0.496409
-0.603763
-0.54025
-0.492351
-0.666944
-0.736168
-0.668618
-0.733944
-1
-0.878096
-0.729327
-1
-0.875
-0.731558
-0.597212
-0.625
-0.509544
-0.374701
-0.375
-0.262228
-0.439501
-0.435693
-0.381762
-0.375151
-0.303302
-0.243334
-0.314755
-0.256748
-0.252575
-0.100238
-0
-0.125
-0
-0.0884843
-0.133844
-0.301451
-0.319617
-0.421308
-0.439579
-0.366542
-0.638837
-0.84411
-0.887525
-0.673669
-0.684397
-0.547681
-0.553932
-0.556981
-0.569094
-0.682474
-0.680908
-0.877814
-0.843209
-0.610964
-0.363741
-0.433969
-0.438104
-0.317325
-0.315253
-0.114408
-0.0955121
-0
-0
-1
-1
-0
-0.5
-1
-0.5
-0.500169
-0
-0.230361
-0.506286
-0.25
-0
-0.25
-0.501817
-0.75
-1
-0.751223
-1
-0.75
-0.245442
-0.249774
-0.359053
-0.482625
-0.339768
-0.25622
-0.246257
-0.342759
-0.493736
-0.737081
-0.598044
-0.732223
-0.628407
-0.755429
-0.618147
-0.736346
-0.380657
-0.367339
-0.619478
-0.619469
-0
-0.112152
-0.214709
-0.125
-0.339768
-0.511468
-0.375
-0.250576
-0.312269
-0.354505
-0.283581
-0.425524
-0.495609
-0.414003
-0.367212
-0.438493
-0.482955
-0.243036
-0.314443
-0
-0.130361
-0.249774
-0
-0.1228
-0.246257
-0
-0.125
-0.243197
-0.375
-0.460998
-0.342759
-0.372076
-0.413223
-0.487797
-0.431307
-0.25198
-0.290425
-0.296938
-0.239051
-0.353925
-0.485506
-0.542001
-0.613511
-0.565178
-0.672486
-0.738741
-0.683983
-0.614503
-0.662211
-0.73167
-0.544579
-0.625
-0.717877
-0.598044
-0.875
-1
-0.883366
-0.732223
-1
-0.888992
-0.736346
-1
-0.868768
-0.606623
-0.739539
-0.671004
-0.56527
-0.611998
-0.559039
-0.741642
-0.68361
-0.618147
-0.760055
-0.625
-1
-0.875
-0.110313
-0.337719
-0.306391
-0.440854
-0.431661
-0.308523
-0.0865872
-0.124325
-0.133215
-0.332392
-0.433604
-0.313603
-0.30387
-0.419458
-0.563411
-0.693606
-0.689503
-0.553108
-0.586083
-0.880316
-0.866348
-0.884747
-0.669969
-0.569241
-0.56415
-0.679016
-0.617255
-0.87398
@@ -32,11 +32,6 @@
 // Compile with: make mesh-optimizer
 //
 // Sample runs:
-//   Adapted analytic Hessian:
-//     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
-//   Adapted discrete size:
-//     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 5 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
-//
 //   Blade shape:
 //     mesh-optimizer -m blade.mesh -o 4 -rs 0 -mid 2 -tid 1 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
 //   Blade limited shape:
@@ -53,9 +48,6 @@
 //     mesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 1000 -ls 2 -li 100 -bnd -qt 1 -qo 8 -cmb
 //   3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
 //   * mesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -rs 0 -mid 303 -tid 1 -ni 20 -ls 2 -li 500 -fix-bnd
-//   2D non-conforming shape and equal size:
-//     mesh-optimizer -m ./amr-quad-q2.mesh -o 2 -rs 1 -mid 9 -tid 2 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
-

 #include "mfem.hpp"
 #include <fstream>
@@ -66,159 +58,189 @@ using namespace std;

 double weight_fun(const Vector &x);

-double ind_values(const Vector &x)
+// Metric values are visualized by creating an L2 finite element functions and
+// computing the metric values at the nodes.
+void vis_metric(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc,
+                Mesh &mesh, char *title, int position)
 {
-   const int opt = 6;
-   const double small = 0.001, big = 0.01;
-
-   // Sine wave.
-   if (opt==1)
-   {
-      const double X = x(0), Y = x(1);
-      const double ind = std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) + 1) -
-                         std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) - 1);
-
-      return ind * small + (1.0 - ind) * big;
-   }
-
-   if (opt==2)
-   {
-      // Circle in the middle.
-      double val = 0.;
-      const double xc = x(0) - 0.5, yc = x(1) - 0.5;
-      const double r = sqrt(xc*xc + yc*yc);
-      double r1 = 0.15; double r2 = 0.35; double sf=30.0;
-      val = 0.5*(std::tanh(sf*(r-r1)) - std::tanh(sf*(r-r2)));
-      if (val > 1.) {val = 1;}
-
-      return val * small + (1.0 - val) * big;
-   }
-
-   if (opt == 3)
-   {
-      // cross
-      const double X = x(0), Y = x(1);
-      const double r1 = 0.45, r2 = 0.55;
-      const double sf = 40.0;
-
-      double val = 0.5 * ( std::tanh(sf*(X-r1)) - std::tanh(sf*(X-r2)) +
-                           std::tanh(sf*(Y-r1)) - std::tanh(sf*(Y-r2)) );
-      if (val > 1.) { val = 1.0; }
-
-      return val * small + (1.0 - val) * big;
-   }
-
-   if (opt==4)
-   {
-      // Multiple circles
-      double r1,r2,val,rval;
-      double sf = 10;
-      val = 0.;
-      // circle 1
-      r1= 0.25; r2 = 0.25; rval = 0.1;
-      double xc = x(0) - r1, yc = x(1) - r2;
-      double r = sqrt(xc*xc+yc*yc);
-      val =  0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
-                                                               (r-rval)));// std::exp(val1);
-      // circle 2
-      r1= 0.75; r2 = 0.75;
-      xc = x(0) - r1, yc = x(1) - r2;
-      r = sqrt(xc*xc+yc*yc);
-      val +=  (0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
-                                                                 (r-rval))));// std::exp(val1);
-      // circle 3
-      r1= 0.75; r2 = 0.25;
-      xc = x(0) - r1, yc = x(1) - r2;
-      r = sqrt(xc*xc+yc*yc);
-      val +=  0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
-                                                                (r-rval)));// std::exp(val1);
-      // circle 4
-      r1= 0.25; r2 = 0.75;
-      xc = x(0) - r1, yc = x(1) - r2;
-      r = sqrt(xc*xc+yc*yc);
-      val +=  0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*(r-rval)));
-      if (val > 1.0) {val = 1.;}
-      if (val < 0.0) {val = 0.;}
-
-      return val * small + (1.0 - val) * big;
-   }
-
-   if (opt==5)
-   {
-      // cross
-      double val = 0.;
-      double X = x(0)-0.5, Y = x(1)-0.5;
-      double rval = std::sqrt(X*X + Y*Y);
-      double thval = 60.*M_PI/180.;
-      double Xmod,Ymod;
-      Xmod = X*std::cos(thval) + Y*std::sin(thval);
-      Ymod= -X*std::sin(thval) + Y*std::cos(thval);
-      X = Xmod+0.5; Y = Ymod+0.5;
-      double r1 = 0.45; double r2 = 0.55; double sf=30.0;
-      val = ( 0.5*(1+std::tanh(sf*(X-r1))) - 0.5*(1+std::tanh(sf*(X-r2)))
-              + 0.5*(1+std::tanh(sf*(Y-r1))) - 0.5*(1+std::tanh(sf*(Y-r2))) );
-      if (rval > 0.4) {val = 0.;}
-      if (val > 1.0) {val = 1.;}
-      if (val < 0.0) {val = 0.;}
-
-      return val * small + (1.0 - val) * big;
-   }
-
-   if (opt==6)
-   {
-      double val = 0.;
-      const double xc = x(0) - 0.0, yc = x(1) - 0.5;
-      const double r = sqrt(xc*xc + yc*yc);
-      double r1 = 0.45; double r2 = 0.55; double sf=30.0;
-      val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
-      if (val > 1.) {val = 1;}
-      if (val < 0.) {val = 0;}
-
-      return val * small + (1.0 - val) * big;
-   }
-
-   return 0.0;
+   L2_FECollection fec(order, mesh.Dimension(), BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec, 1);
+   GridFunction metric(&fes);
+   InterpolateTMOP_QualityMetric(qm, tc, mesh, metric);
+   osockstream sock(19916, "localhost");
+   sock << "solution\n";
+   mesh.Print(sock);
+   metric.Save(sock);
+   sock.send();
+   sock << "window_title '"<< title << "'\n"
+        << "window_geometry "
+        << position << " " << 0 << " " << 600 << " " << 600 << "\n"
+        << "keys jRmclA" << endl;
 }

-class HessianCoefficient : public MatrixCoefficient
+class RelaxedNewtonSolver : public NewtonSolver
 {
 private:
-   int type;
+   // Quadrature points that are checked for negative Jacobians etc.
+   const IntegrationRule &ir;
+   FiniteElementSpace *fes;
+   mutable GridFunction x_gf;

 public:
-   HessianCoefficient(int dim, int type_)
-      : MatrixCoefficient(dim), type(type_) { }
+   RelaxedNewtonSolver(const IntegrationRule &irule, FiniteElementSpace *f)
+      : ir(irule), fes(f) { }

-   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
+};
+
+double RelaxedNewtonSolver::ComputeScalingFactor(const Vector &x,
+                                                 const Vector &b) const
+{
+   const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
+   MFEM_VERIFY(nlf != NULL, "invalid Operator subclass");
+   const bool have_b = (b.Size() == Height());
+
+   const int NE = fes->GetMesh()->GetNE(), dim = fes->GetFE(0)->GetDim(),
+             dof = fes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
+   Array<int> xdofs(dof * dim);
+   DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
+   Vector posV(pos.Data(), dof * dim);
+
+   Vector x_out(x.Size());
+   bool x_out_ok = false;
+   const double energy_in = nlf->GetEnergy(x);
+   double scale = 1.0, energy_out;
+   double norm0 = Norm(r);
+   x_gf.MakeTRef(fes, x_out, 0);
+
+   // Decreases the scaling of the update until the new mesh is valid.
+   for (int i = 0; i < 12; i++)
   {
-      Vector pos(3);
-      T.Transform(ip, pos);
+      add(x, -scale, c, x_out);
+      x_gf.SetFromTrueVector();

-      if (type == 0)
+      energy_out = nlf->GetGridFunctionEnergy(x_gf);
+      if (energy_out > 1.2*energy_in || std::isnan(energy_out) != 0)
      {
-         K(0, 0) = 1.0 + 3.0 * std::sin(M_PI*pos(0));
-         K(0, 1) = 0.0;
-         K(1, 0) = 0.0;
-         K(1, 1) = 1.0;
+         if (print_level >= 0)
+         { cout << "Scale = " << scale << " Increasing energy." << endl; }
+         scale *= 0.5; continue;
      }
-      else
+
+      int jac_ok = 1;
+      for (int i = 0; i < NE; i++)
      {
-         const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
-         const double r = sqrt(xc*xc + yc*yc);
-         double r1 = 0.15; double r2 = 0.35; double sf=30.0;
-         const double eps = 0.5;
+         fes->GetElementVDofs(i, xdofs);
+         x_gf.GetSubVector(xdofs, posV);
+         for (int j = 0; j < nsp; j++)
+         {
+            fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
+            MultAtB(pos, dshape, Jpr);
+            if (Jpr.Det() <= 0.0) { jac_ok = 0; goto break2; }
+         }
+      }
+   break2:
+      if (jac_ok == 0)
+      {
+         if (print_level >= 0)
+         { cout << "Scale = " << scale << " Neg det(J) found." << endl; }
+         scale *= 0.5; continue;
+      }

-         const double tan1 = std::tanh(sf*(r-r1)),
-                      tan2 = std::tanh(sf*(r-r2));
+      oper->Mult(x_out, r);
+      if (have_b) { r -= b; }
+      double norm = Norm(r);

-         K(0, 0) = eps + 1.0 * (tan1 - tan2);
-         K(0, 1) = 0.0;
-         K(1, 0) = 0.0;
-         K(1, 1) = 1.0;
+      if (norm > 1.2*norm0)
+      {
+         if (print_level >= 0)
+         { cout << "Scale = " << scale << " Norm increased." << endl; }
+         scale *= 0.5; continue;
+      }
+      else { x_out_ok = true; break; }
+   }
+
+   if (print_level >= 0)
+   {
+      cout << "Energy decrease: "
+           << (energy_in - energy_out) / energy_in * 100.0
+           << "% with " << scale << " scaling." << endl;
+   }
+
+   if (x_out_ok == false) { scale = 0.0; }
+
+   return scale;
+}
+
+// Allows negative Jacobians. Used in untangling metrics.
+class DescentNewtonSolver : public NewtonSolver
+{
+private:
+   // Quadrature points that are checked for negative Jacobians etc.
+   const IntegrationRule &ir;
+   FiniteElementSpace *fes;
+   mutable GridFunction x_gf;
+
+public:
+   DescentNewtonSolver(const IntegrationRule &irule, FiniteElementSpace *f)
+      : ir(irule), fes(f) { }
+
+   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
+};
+
+double DescentNewtonSolver::ComputeScalingFactor(const Vector &x,
+                                                 const Vector &b) const
+{
+   const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
+   MFEM_VERIFY(nlf != NULL, "invalid Operator subclass");
+
+   const int NE = fes->GetMesh()->GetNE(), dim = fes->GetFE(0)->GetDim(),
+             dof = fes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
+   Array<int> xdofs(dof * dim);
+   DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
+   Vector posV(pos.Data(), dof * dim);
+
+   x_gf.MakeTRef(fes, x.GetData());
+   x_gf.SetFromTrueVector();
+
+   double min_detJ = infinity();
+   for (int i = 0; i < NE; i++)
+   {
+      fes->GetElementVDofs(i, xdofs);
+      x_gf.GetSubVector(xdofs, posV);
+      for (int j = 0; j < nsp; j++)
+      {
+         fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
+         MultAtB(pos, dshape, Jpr);
+         min_detJ = min(min_detJ, Jpr.Det());
      }
   }
-};
+   cout << "Minimum det(J) = " << min_detJ << endl;
+
+   Vector x_out(x.Size());
+   bool x_out_ok = false;
+   const double energy_in = nlf->GetGridFunctionEnergy(x_gf);
+   double scale = 1.0, energy_out;
+
+   for (int i = 0; i < 7; i++)
+   {
+      add(x, -scale, c, x_out);
+
+      energy_out = nlf->GetEnergy(x_out);
+      if (energy_out > energy_in || std::isnan(energy_out) != 0)
+      {
+         scale *= 0.5;
+      }
+      else { x_out_ok = true; break; }
+   }
+
+   cout << "Energy decrease: " << (energy_in - energy_out) / energy_in * 100.0
+        << "% with " << scale << " scaling." << endl;
+
+   if (x_out_ok == false) { return 0.0; }
+
+   return scale;
+}

 // Additional IntegrationRules that can be used with the --quad-type option.
 IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto);
@@ -282,9 +304,7 @@ int main (int argc, char *argv[])
                  "Target (ideal element) type:\n\t"
                  "1: Ideal shape, unit size\n\t"
                  "2: Ideal shape, equal size\n\t"
-                  "3: Ideal shape, initial size\n\t"
-                  "4: Given full analytic Jacobian (in physical space)\n\t"
-                  "5: Ideal shape, given size (in physical space)");
+                  "3: Ideal shape, initial size");
   args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
   args.AddOption(&quad_type, "-qt", "--quad-type",
                  "Quadrature rule type:\n\t"
@@ -355,8 +375,7 @@ int main (int argc, char *argv[])
   // 6. Get the mesh nodes (vertices and other degrees of freedom in the finite
   //    element space) as a finite element grid function in fespace. Note that
   //    changing x automatically changes the shapes of the mesh elements.
-   GridFunction x(fespace);
-   mesh->SetNodalGridFunction(&x);
+   GridFunction *x = mesh->GetNodes();

   // 7. 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
@@ -406,9 +425,10 @@ int main (int argc, char *argv[])
      // Set the boundary values to zero.
      for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
   }
-   x -= rdm;
-   x.SetTrueVector();
-   x.SetFromTrueVector();
+   *x -= rdm;
+   // Set the perturbation of all nodes from the true nodes.
+   x->SetTrueVector();
+   x->SetFromTrueVector();

   // 9. Save the starting (prior to the optimization) mesh to a file. This
   //    output can be viewed later using GLVis: "glvis -m perturbed.mesh".
@@ -419,7 +439,7 @@ int main (int argc, char *argv[])

   // 10. Store the starting (prior to the optimization) positions.
   GridFunction x0(fespace);
-   x0 = x;
+   x0 = *x;

   // 11. Form the integrator that uses the chosen metric and target.
   double tauval = -0.1;
@@ -448,43 +468,15 @@ int main (int argc, char *argv[])
      default: cout << "Unknown metric_id: " << metric_id << endl; return 3;
   }
   TargetConstructor::TargetType target_t;
-   TargetConstructor *target_c = NULL;
-   HessianCoefficient *adapt_coeff = NULL;
-   H1_FECollection ind_fec(3, dim);
-   FiniteElementSpace ind_fes(mesh, &ind_fec);
-   GridFunction size;
   switch (target_id)
   {
      case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
      case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
      case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
-      case 4:
-      {
-         target_t = TargetConstructor::GIVEN_FULL;
-         AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
-         adapt_coeff = new HessianCoefficient(dim, 1);
-         tc->SetAnalyticTargetSpec(NULL, NULL, adapt_coeff);
-         target_c = tc;
-         break;
-      }
-      case 5:
-      {
-         target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE;
-         DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
-         size.SetSpace(&ind_fes);
-         FunctionCoefficient ind_coeff(ind_values);
-         size.ProjectCoefficient(ind_coeff);
-         tc->SetSerialDiscreteTargetSpec(size);
-         target_c = tc;
-         break;
-      }
-      default: cout << "Unknown target_id: " << target_id << endl; return 3;
-   }
-
-   if (target_c == NULL)
-   {
-      target_c = new TargetConstructor(target_t);
+      default: cout << "Unknown target_id: " << target_id << endl;
+         delete metric; return 3;
   }
+   TargetConstructor *target_c = new TargetConstructor(target_t);
   target_c->SetNodes(x0);
   TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c);

@@ -550,13 +542,13 @@ int main (int argc, char *argv[])
   }
   else { a.AddDomainIntegrator(he_nlf_integ); }

-   const double init_energy = a.GetGridFunctionEnergy(x);
+   const double init_energy = a.GetGridFunctionEnergy(*x);

   // 15. Visualize the starting mesh and metric values.
   if (visualization)
   {
      char title[] = "Initial metric values";
-      vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
+      vis_metric(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
   }

   // 16. Fix all boundary nodes, or fix only a given component depending on the
@@ -660,13 +652,8 @@ int main (int argc, char *argv[])
   if (tauval > 0.0)
   {
      tauval = 0.0;
-      TMOPNewtonSolver *tns = new TMOPNewtonSolver(*ir);
-      if (target_id == 5)
-      {
-         tns->SetDiscreteAdaptTC(dynamic_cast<DiscreteAdaptTC *>(target_c));
-      }
-      newton = tns;
-      cout << "TMOPNewtonSolver is used (as all det(J) > 0).\n";
+      newton = new RelaxedNewtonSolver(*ir, fespace);
+      cout << "The RelaxedNewtonSolver is used (as all det(J)>0)." << endl;
   }
   else
   {
@@ -677,8 +664,8 @@ int main (int argc, char *argv[])
         return 3;
      }
      tauval -= 0.01 * h0.Min(); // Slightly below minJ0 to avoid div by 0.
-      newton = new TMOPDescentNewtonSolver(*ir);
-      cout << "The TMOPDescentNewtonSolver is used (as some det(J) < 0).\n";
+      newton = new DescentNewtonSolver(*ir, fespace);
+      cout << "The DescentNewtonSolver is used (as some det(J)<0)." << endl;
   }
   newton->SetPreconditioner(*S);
   newton->SetMaxIter(newton_iter);
@@ -686,9 +673,8 @@ int main (int argc, char *argv[])
   newton->SetAbsTol(0.0);
   newton->SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
   newton->SetOperator(a);
-   newton->Mult(b, x.GetTrueVector());
-   x.SetFromTrueVector();
-
+   newton->Mult(b, x->GetTrueVector());
+   x->SetFromTrueVector();
   if (newton->GetConverged() == false)
   {
      cout << "NewtonIteration: rtol = " << newton_rtol << " not achieved."
@@ -705,12 +691,12 @@ int main (int argc, char *argv[])
   }

   // 21. Compute the amount of energy decrease.
-   const double fin_energy = a.GetGridFunctionEnergy(x);
+   const double fin_energy = a.GetGridFunctionEnergy(*x);
   double metric_part = fin_energy;
   if (lim_const != 0.0)
   {
      lim_coeff.constant = 0.0;
-      metric_part = a.GetGridFunctionEnergy(x);
+      metric_part = a.GetGridFunctionEnergy(*x);
      lim_coeff.constant = lim_const;
   }
   cout << "Initial strain energy: " << init_energy
@@ -726,13 +712,13 @@ int main (int argc, char *argv[])
   if (visualization)
   {
      char title[] = "Final metric values";
-      vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 600);
+      vis_metric(mesh_poly_deg, *metric, *target_c, *mesh, title, 600);
   }

   // 23. Visualize the mesh displacement.
   if (visualization)
   {
-      x0 -= x;
+      x0 -= *x;
      osockstream sock(19916, "localhost");
      sock << "solution\n";
      mesh->Print(sock);
@@ -32,11 +32,6 @@
 // Compile with: make pmesh-optimizer
 //
 // Sample runs:
-//   Adapted analytic Hessian:
-//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
-//   Adapted discrete size:
-//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 5 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
-//
 //   Blade shape:
 //     mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -rs 0 -mid 2 -tid 1 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
 //   Blade limited shape:
@@ -53,176 +48,223 @@
 //     mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 1000 -ls 2 -li 100 -bnd -qt 1 -qo 8 -cmb
 //   3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
 //   * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -rs 0 -mid 303 -tid 1 -ni 20 -ls 2 -li 500 -fix-bnd
-//   2D non-conforming shape and equal size:
-//     mpirun -np 4 pmesh-optimizer -m ./amr-quad-q2.mesh -o 2 -rs 1 -mid 9 -tid 2 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8

 #include "mfem.hpp"
-#include <iostream>
 #include <fstream>
+#include <iostream>

 using namespace mfem;
 using namespace std;

 double weight_fun(const Vector &x);

-double ind_values(const Vector &x)
+// Metric values are visualized by creating an L2 finite element functions and
+// computing the metric values at the nodes.
+void vis_metric(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc,
+                ParMesh &pmesh, char *title, int position)
 {
-   const int opt = 6;
-   const double small = 0.001, big = 0.01;
-
-   // Sine wave.
-   if (opt==1)
+   L2_FECollection fec(order, pmesh.Dimension(), BasisType::GaussLobatto);
+   ParFiniteElementSpace fes(&pmesh, &fec, 1);
+   ParGridFunction metric(&fes);
+   InterpolateTMOP_QualityMetric(qm, tc, pmesh, metric);
+   socketstream sock;
+   if (pmesh.GetMyRank() == 0)
   {
-      const double X = x(0), Y = x(1);
-      const double ind = std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) + 1) -
-                         std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) - 1);
-
-      return ind * small + (1.0 - ind) * big;
+      sock.open("localhost", 19916);
+      sock << "solution\n";
   }
-
-   if (opt==2)
+   pmesh.PrintAsOne(sock);
+   metric.SaveAsOne(sock);
+   if (pmesh.GetMyRank() == 0)
   {
-      // Circle in the middle.
-      double val = 0.;
-      const double xc = x(0) - 0.5, yc = x(1) - 0.5;
-      const double r = sqrt(xc*xc + yc*yc);
-      double r1 = 0.15; double r2 = 0.35; double sf=30.0;
-      val = 0.5*(std::tanh(sf*(r-r1)) - std::tanh(sf*(r-r2)));
-      if (val > 1.) {val = 1;}
-
-      return val * small + (1.0 - val) * big;
+      sock << "window_title '"<< title << "'\n"
+           << "window_geometry "
+           << position << " " << 0 << " " << 600 << " " << 600 << "\n"
+           << "keys jRmclA" << endl;
   }
-
-   if (opt == 3)
-   {
-      // cross
-      const double X = x(0), Y = x(1);
-      const double r1 = 0.45, r2 = 0.55;
-      const double sf = 40.0;
-
-      double val = 0.5 * ( std::tanh(sf*(X-r1)) - std::tanh(sf*(X-r2)) +
-                           std::tanh(sf*(Y-r1)) - std::tanh(sf*(Y-r2)) );
-      if (val > 1.) { val = 1.0; }
-
-      return val * small + (1.0 - val) * big;
-   }
-
-   if (opt==4)
-   {
-      // Multiple circles
-      double r1,r2,val,rval;
-      double sf = 10;
-      val = 0.;
-      // circle 1
-      r1= 0.25; r2 = 0.25; rval = 0.1;
-      double xc = x(0) - r1, yc = x(1) - r2;
-      double r = sqrt(xc*xc+yc*yc);
-      val =  0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
-                                                               (r-rval)));// std::exp(val1);
-      // circle 2
-      r1= 0.75; r2 = 0.75;
-      xc = x(0) - r1, yc = x(1) - r2;
-      r = sqrt(xc*xc+yc*yc);
-      val +=  (0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
-                                                                 (r-rval))));// std::exp(val1);
-      // circle 3
-      r1= 0.75; r2 = 0.25;
-      xc = x(0) - r1, yc = x(1) - r2;
-      r = sqrt(xc*xc+yc*yc);
-      val +=  0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
-                                                                (r-rval)));// std::exp(val1);
-      // circle 4
-      r1= 0.25; r2 = 0.75;
-      xc = x(0) - r1, yc = x(1) - r2;
-      r = sqrt(xc*xc+yc*yc);
-      val +=  0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*(r-rval)));
-      if (val > 1.0) {val = 1.;}
-      if (val < 0.0) {val = 0.;}
-
-      return val * small + (1.0 - val) * big;
-   }
-
-   if (opt==5)
-   {
-      // cross
-      double val = 0.;
-      double X = x(0)-0.5, Y = x(1)-0.5;
-      double rval = std::sqrt(X*X + Y*Y);
-      double thval = 60.*M_PI/180.;
-      double Xmod,Ymod;
-      Xmod = X*std::cos(thval) + Y*std::sin(thval);
-      Ymod= -X*std::sin(thval) + Y*std::cos(thval);
-      X = Xmod+0.5; Y = Ymod+0.5;
-      double r1 = 0.45; double r2 = 0.55; double sf=30.0;
-      val = ( 0.5*(1+std::tanh(sf*(X-r1))) - 0.5*(1+std::tanh(sf*(X-r2)))
-              + 0.5*(1+std::tanh(sf*(Y-r1))) - 0.5*(1+std::tanh(sf*(Y-r2))) );
-      if (rval > 0.4) {val = 0.;}
-      if (val > 1.0) {val = 1.;}
-      if (val < 0.0) {val = 0.;}
-
-      return val * small + (1.0 - val) * big;
-   }
-
-   if (opt==6)
-   {
-      double val = 0.;
-      const double xc = x(0) - 0.0, yc = x(1) - 0.5;
-      const double r = sqrt(xc*xc + yc*yc);
-      double r1 = 0.45; double r2 = 0.55; double sf=30.0;
-      val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
-      if (val > 1.) {val = 1;}
-      if (val < 0.) {val = 0;}
-
-      return val * small + (1.0 - val) * big;
-   }
-
-   return 0.0;
 }

-class HessianCoefficient : public MatrixCoefficient
+class RelaxedNewtonSolver : public NewtonSolver
 {
 private:
-   int type;
+   // Quadrature points that are checked for negative Jacobians etc.
+   const IntegrationRule &ir;
+   ParFiniteElementSpace *pfes;
+   mutable ParGridFunction x_gf;

 public:
-   HessianCoefficient(int dim, int type_)
-      : MatrixCoefficient(dim), type(type_) { }
+   RelaxedNewtonSolver(const IntegrationRule &irule, ParFiniteElementSpace *pf)
+      : NewtonSolver(pf->GetComm()), ir(irule), pfes(pf) { }

-   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
-                     const IntegrationPoint &ip)
+   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
+};
+
+double RelaxedNewtonSolver::ComputeScalingFactor(const Vector &x,
+                                                 const Vector &b) const
+{
+   const ParNonlinearForm *nlf = dynamic_cast<const ParNonlinearForm *>(oper);
+   MFEM_VERIFY(nlf != NULL, "invalid Operator subclass");
+   const bool have_b = (b.Size() == Height());
+
+   const int NE = pfes->GetParMesh()->GetNE(), dim = pfes->GetFE(0)->GetDim(),
+             dof = pfes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
+   Array<int> xdofs(dof * dim);
+   DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
+   Vector posV(pos.Data(), dof * dim);
+
+   Vector x_out(x.Size());
+   bool x_out_ok = false;
+   const double energy_in = nlf->GetEnergy(x);
+   double scale = 1.0, energy_out;
+   double norm0 = Norm(r);
+   x_gf.MakeTRef(pfes, x_out, 0);
+
+   // Decreases the scaling of the update until the new mesh is valid.
+   for (int i = 0; i < 12; i++)
   {
-      Vector pos(3);
-      T.Transform(ip, pos);
+      add(x, -scale, c, x_out);
+      x_gf.SetFromTrueVector();

-      if (type == 0)
+      energy_out = nlf->GetParGridFunctionEnergy(x_gf);
+      if (energy_out > 1.2*energy_in || std::isnan(energy_out) != 0)
      {
-         K(0, 0) = 1.0 + 3.0 * std::sin(M_PI*pos(0));
-         K(0, 1) = 0.0;
-         K(1, 0) = 0.0;
-         K(1, 1) = 1.0;
+         if (print_level >= 0)
+         { cout << "Scale = " << scale << " Increasing energy." << endl; }
+         scale *= 0.5; continue;
      }
-      else
+
+      int jac_ok = 1;
+      for (int i = 0; i < NE; i++)
      {
-         const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
-         const double r = sqrt(xc*xc + yc*yc);
-         double r1 = 0.15; double r2 = 0.35; double sf=30.0;
-         const double eps = 0.5;
+         pfes->GetElementVDofs(i, xdofs);
+         x_gf.GetSubVector(xdofs, posV);
+         for (int j = 0; j < nsp; j++)
+         {
+            pfes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
+            MultAtB(pos, dshape, Jpr);
+            if (Jpr.Det() <= 0.0) { jac_ok = 0; goto break2; }
+         }
+      }
+   break2:
+      int jac_ok_all;
+      MPI_Allreduce(&jac_ok, &jac_ok_all, 1, MPI_INT, MPI_LAND,
+                    pfes->GetComm());

-         const double tan1 = std::tanh(sf*(r-r1)),
-                      tan2 = std::tanh(sf*(r-r2));
+      if (jac_ok_all == 0)
+      {
+         if (print_level >= 0)
+         { cout << "Scale = " << scale << " Neg det(J) found." << endl; }
+         scale *= 0.5; continue;
+      }

-         K(0, 0) = eps + 1.0 * (tan1 - tan2);
-         K(0, 1) = 0.0;
-         K(1, 0) = 0.0;
-         K(1, 1) = 1.0;
+      oper->Mult(x_out, r);
+      if (have_b) { r -= b; }
+      double norm = Norm(r);
+
+      if (norm > 1.2*norm0)
+      {
+         if (print_level >= 0)
+         { cout << "Scale = " << scale << " Norm increased." << endl; }
+         scale *= 0.5; continue;
+      }
+      else { x_out_ok = true; break; }
+   }
+
+   if (print_level >= 0)
+   {
+      cout << "Energy decrease: "
+           << (energy_in - energy_out) / energy_in * 100.0
+           << "% with " << scale << " scaling." << endl;
+   }
+
+   if (x_out_ok == false) { scale = 0.0; }
+
+   return scale;
+}
+
+// Allows negative Jacobians. Used in untangling metrics.
+class DescentNewtonSolver : public NewtonSolver
+{
+private:
+   // Quadrature points that are checked for negative Jacobians etc.
+   const IntegrationRule &ir;
+   ParFiniteElementSpace *pfes;
+   mutable ParGridFunction x_gf;
+
+public:
+   DescentNewtonSolver(const IntegrationRule &irule, ParFiniteElementSpace *pf)
+      : NewtonSolver(pf->GetComm()), ir(irule), pfes(pf) { }
+
+   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
+};
+
+double DescentNewtonSolver::ComputeScalingFactor(const Vector &x,
+                                                 const Vector &b) const
+{
+   const ParNonlinearForm *nlf = dynamic_cast<const ParNonlinearForm *>(oper);
+   MFEM_VERIFY(nlf != NULL, "invalid Operator subclass");
+
+   const int NE = pfes->GetParMesh()->GetNE(), dim = pfes->GetFE(0)->GetDim(),
+             dof = pfes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
+   Array<int> xdofs(dof * dim);
+   DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
+   Vector posV(pos.Data(), dof * dim);
+
+   x_gf.MakeTRef(pfes, x.GetData());
+   x_gf.SetFromTrueVector();
+
+   double min_detJ = infinity();
+   for (int i = 0; i < NE; i++)
+   {
+      pfes->GetElementVDofs(i, xdofs);
+      x_gf.GetSubVector(xdofs, posV);
+      for (int j = 0; j < nsp; j++)
+      {
+         pfes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
+         MultAtB(pos, dshape, Jpr);
+         min_detJ = min(min_detJ, Jpr.Det());
      }
   }
-};
+   double min_detJ_all;
+   MPI_Allreduce(&min_detJ, &min_detJ_all, 1, MPI_DOUBLE, MPI_MIN,
+                 pfes->GetComm());
+   if (print_level >= 0)
+   { cout << "Minimum det(J) = " << min_detJ_all << endl; }
+
+   Vector x_out(x.Size());
+   bool x_out_ok = false;
+   const double energy_in = nlf->GetParGridFunctionEnergy(x_gf);
+   double scale = 1.0, energy_out;
+
+   for (int i = 0; i < 7; i++)
+   {
+      add(x, -scale, c, x_out);
+
+      energy_out = nlf->GetEnergy(x_out);
+      if (energy_out > energy_in || std::isnan(energy_out) != 0)
+      {
+         scale *= 0.5;
+      }
+      else { x_out_ok = true; break; }
+   }
+
+   if (print_level >= 0)
+   {
+      cout << "Energy decrease: "
+           << (energy_in - energy_out) / energy_in * 100.0
+           << "% with " << scale << " scaling." << endl;
+   }
+
+   if (x_out_ok == false) { return 0.0; }
+
+   return scale;
+}

 // Additional IntegrationRules that can be used with the --quad-type option.
 IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto);
 IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform);

+
 int main (int argc, char *argv[])
 {
   // 0. Initialize MPI.
@@ -289,9 +331,7 @@ int main (int argc, char *argv[])
                  "Target (ideal element) type:\n\t"
                  "1: Ideal shape, unit size\n\t"
                  "2: Ideal shape, equal size\n\t"
-                  "3: Ideal shape, initial size\n\t"
-                  "4: Given full analytic Jacobian (in physical space)\n\t"
-                  "5: Ideal shape, given size (in physical space)");
+                  "3: Ideal shape, initial size");
   args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
   args.AddOption(&quad_type, "-qt", "--quad-type",
                  "Quadrature rule type:\n\t"
@@ -341,7 +381,6 @@ int main (int argc, char *argv[])
      cout << endl;
   }
   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
-
   delete mesh;
   for (int lev = 0; lev < rp_levels; lev++) { pmesh->UniformRefinement(); }

@@ -471,47 +510,20 @@ int main (int argc, char *argv[])
         return 3;
   }
   TargetConstructor::TargetType target_t;
-   TargetConstructor *target_c = NULL;
-   HessianCoefficient *adapt_coeff = NULL;
-   H1_FECollection ind_fec(3, dim);
-   ParFiniteElementSpace ind_fes(pmesh, &ind_fec);
-   ParGridFunction size;
   switch (target_id)
   {
      case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
      case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
      case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
-      case 4:
-      {
-         target_t = TargetConstructor::GIVEN_FULL;
-         AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
-         adapt_coeff = new HessianCoefficient(dim, 1);
-         tc->SetAnalyticTargetSpec(NULL, NULL, adapt_coeff);
-         target_c = tc;
-         break;
-      }
-      case 5:
-      {
-         target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE;
-         DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
-         size.SetSpace(&ind_fes);
-         FunctionCoefficient ind_coeff(ind_values);
-         size.ProjectCoefficient(ind_coeff);
-         tc->SetParDiscreteTargetSpec(size);
-         target_c = tc;
-         break;
-      }
      default:
         if (myid == 0) { cout << "Unknown target_id: " << target_id << endl; }
         return 3;
   }
-
-   if (target_c == NULL)
-   {
-      target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
-   }
+   TargetConstructor *target_c;
+   target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
   target_c->SetNodes(x0);
-   TMOP_Integrator *he_nlf_integ= new TMOP_Integrator(metric, target_c);
+   TMOP_Integrator *he_nlf_integ;
+   he_nlf_integ = new TMOP_Integrator(metric, target_c);

   // 13. Setup the quadrature rule for the non-linear form integrator.
   const IntegrationRule *ir = NULL;
@@ -585,7 +597,7 @@ int main (int argc, char *argv[])
   if (visualization)
   {
      char title[] = "Initial metric values";
-      vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
+      vis_metric(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
   }

   // 17. Fix all boundary nodes, or fix only a given component depending on the
@@ -693,14 +705,9 @@ int main (int argc, char *argv[])
   if (tauval > 0.0)
   {
      tauval = 0.0;
-      TMOPNewtonSolver *tns = new TMOPNewtonSolver(pfespace->GetComm(), *ir);
-      if (target_id == 5)
-      {
-         tns->SetDiscreteAdaptTC(dynamic_cast<DiscreteAdaptTC *>(target_c));
-      }
-      newton = tns;
+      newton = new RelaxedNewtonSolver(*ir, pfespace);
      if (myid == 0)
-      { cout << "TMOPNewtonSolver is used (as all det(J) > 0)." << endl; }
+      { cout << "RelaxedNewtonSolver is used (as all det(J) > 0)." << endl; }
   }
   else
   {
@@ -708,15 +715,15 @@ int main (int argc, char *argv[])
           (dim == 3 && metric_id != 352) )
      {
         if (myid == 0)
-         { cout << "The mesh is inverted. Use an untangling metric.\n"; }
+         { cout << "The mesh is inverted. Use an untangling metric." << endl; }
         return 3;
      }
      double h0min = h0.Min(), h0min_all;
      MPI_Allreduce(&h0min, &h0min_all, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
      tauval -= 0.01 * h0min_all; // Slightly below minJ0 to avoid div by 0.
-      newton = new TMOPDescentNewtonSolver(pfespace->GetComm(), *ir);
+      newton = new DescentNewtonSolver(*ir, pfespace);
      if (myid == 0)
-      { cout << "TMOPDescentNewtonSolver is used (as some det(J) < 0).\n"; }
+      { cout << "DescentNewtonSolver is used (as some det(J) < 0)." << endl; }
   }
   newton->SetPreconditioner(*S);
   newton->SetMaxIter(newton_iter);
@@ -768,7 +775,7 @@ int main (int argc, char *argv[])
   if (visualization)
   {
      char title[] = "Final metric values";
-      vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 600);
+      vis_metric(mesh_poly_deg, *metric, *target_c, *pmesh, title, 600);
   }

   // 23. Visualize the mesh displacement.
@@ -798,7 +805,6 @@ int main (int argc, char *argv[])
   delete metric2;
   delete coeff1;
   delete target_c;
-   delete adapt_coeff;
   delete metric;
   delete pfespace;
   delete fec;
@@ -1,114 +0,0 @@
-MFEM mesh v1.0
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-# PRISM       = 6
-#
-
-dimension
-2
-
-elements
-16
-1 3 0 1 6 5
-1 3 1 2 7 6
-1 3 2 3 8 7
-1 3 3 4 9 8
-1 3 5 6 11 10
-1 3 6 7 12 11
-1 3 7 8 13 12
-1 3 8 9 14 13
-1 3 10 11 16 15
-1 3 11 12 17 16
-1 3 12 13 18 17
-1 3 13 14 19 18
-1 3 15 16 21 20
-1 3 16 17 22 21
-1 3 17 18 23 22
-1 3 18 19 24 23
-
-boundary
-16
-2 1 0 1
-2 1 1 2
-2 1 2 3
-2 1 3 4
-2 1 21 20
-2 1 22 21
-2 1 23 22
-2 1 24 23
-1 1 5 0
-1 1 10 5
-1 1 15 10
-1 1 20 15
-1 1 4 9
-1 1 9 14
-1 1 14 19
-1 1 19 24
-
-vertices
-25
-
-nodes
-FiniteElementSpace
-FiniteElementCollection: H1_2D_P1
-VDim: 2
-Ordering: 0
-
-0
-0.25
-0.5
-0.75
-1
-0
-0.25
-0.5
-0.75
-1
-0
-0.25
-0.5
-0.75
-1
-0
-0.25
-0.5
-0.75
-1
-0
-0.25
-0.5
-0.75
-1
-0
-0
-0
-0
-0
-0.25
-0.25
-0.25
-0.25
-0.25
-0.5
-0.5
-0.5
-0.5
-0.5
-0.75
-0.75
-0.75
-0.75
-0.75
-1
-1
-1
-1
-1
@@ -1,147 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM 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) version 2.1 dated February 1999.
-
-#include "catch.hpp"
-#include "mfem.hpp"
-
-using namespace mfem;
-
-namespace pa_coeff
-{
-
-int dimension;
-
-double coeffFunction(const Vector& x)
-{
-   if (dimension == 2)
-   {
-      return sin(8.0 * M_PI * x[0]) * cos(6.0 * M_PI * x[1]) + 2.0;
-   }
-   else
-   {
-      return sin(8.0 * M_PI * x[0]) * cos(6.0 * M_PI * x[1]) *
-             sin(4.0 * M_PI * x[2]) +
-             2.0;
-   }
-}
-
-TEST_CASE("pa_coeff")
-{
-   for (dimension = 2; dimension < 4; ++dimension)
-   {
-      for (int coeffType = 0; coeffType < 3; ++coeffType)
-      {
-         for (int integrator = 0; integrator < 2; ++integrator)
-         {
-            const int ne = 2;
-            std::cout << "Testing " << dimension << "D partial assembly with "
-                      << "coeffType " << coeffType << " and "
-                      << "integrator " << integrator << std::endl;
-            for (int order = 1; order < 5; ++order)
-            {
-               Mesh* mesh;
-               if (dimension == 2)
-               {
-                  mesh = new Mesh(ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
-               }
-               else
-               {
-                  mesh = new Mesh(ne, ne, ne, Element::HEXAHEDRON, 1, 1.0, 1.0,
-                                  1.0);
-               }
-               FiniteElementCollection* h1_fec =
-                  new H1_FECollection(order, dimension);
-               FiniteElementSpace h1_fespace(mesh, h1_fec);
-               Array<int> ess_tdof_list;
-
-               BilinearForm paform(&h1_fespace);
-               GridFunction* coeffGridFunction = nullptr;
-               Coefficient* coeff = nullptr;
-               if (coeffType == 0)
-               {
-                  coeff = new ConstantCoefficient(1.0);
-               }
-               else if (coeffType == 1)
-               {
-                  coeff = new FunctionCoefficient(&coeffFunction);
-               }
-               else if (coeffType == 2)
-               {
-                  FunctionCoefficient tmpCoeff(&coeffFunction);
-                  coeffGridFunction = new GridFunction(&h1_fespace);
-                  coeffGridFunction->ProjectCoefficient(tmpCoeff);
-                  coeff = new GridFunctionCoefficient(coeffGridFunction);
-               }
-               paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-               if (integrator < 2)
-               {
-                  paform.AddDomainIntegrator(new DiffusionIntegrator(*coeff));
-               }
-               if (integrator > 0)
-               {
-                  paform.AddDomainIntegrator(new MassIntegrator(*coeff));
-               }
-               paform.Assemble();
-               OperatorHandle paopr;
-               paform.FormSystemMatrix(ess_tdof_list, paopr);
-
-               BilinearForm assemblyform(&h1_fespace);
-               if (integrator < 2)
-               {
-                  assemblyform.AddDomainIntegrator(
-                     new DiffusionIntegrator(*coeff));
-               }
-               if (integrator > 0)
-               {
-                  assemblyform.AddDomainIntegrator(new MassIntegrator(*coeff));
-               }
-               assemblyform.SetDiagonalPolicy(Matrix::DIAG_ONE);
-               assemblyform.Assemble();
-               assemblyform.Finalize();
-               const SparseMatrix& A_explicit = assemblyform.SpMat();
-
-               Vector xin(h1_fespace.GetTrueVSize());
-               xin.Randomize();
-               Vector y_mat(xin);
-               y_mat = 0.0;
-               Vector y_assembly(xin);
-               y_assembly = 0.0;
-               Vector y_pa(xin);
-               y_pa = 0.0;
-
-               paopr->Mult(xin, y_pa);
-               assemblyform.Mult(xin, y_assembly);
-               A_explicit.Mult(xin, y_mat);
-
-               y_pa -= y_mat;
-               double pa_error = y_pa.Norml2();
-               std::cout << "  order: " << order
-                         << ", pa error norm: " << pa_error << std::endl;
-               REQUIRE(pa_error < 1.e-12);
-
-               y_assembly -= y_mat;
-               double assembly_error = y_assembly.Norml2();
-               std::cout << "  order: " << order
-                         << ", assembly error norm: " << assembly_error
-                         << std::endl;
-               REQUIRE(assembly_error < 1.e-12);
-
-               delete coeff;
-               delete coeffGridFunction;
-               delete mesh;
-               delete h1_fec;
-            }
-         }
-      }
-   }
-}
-
-} // namespace pa_coeff