diff --git a/CHANGELOG b/CHANGELOG
index f2ce986591..73ff0de018 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -18,13 +18,7 @@ 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.
-
-- Improved multi-GPU MPI communication.
-
-Miscellaneous
--------------
-- Improved element numbering after uniform mesh refinement.
+- Improved RAJA backend and multi-GPU MPI communications.
 
 Discretization improvements
 ---------------------------
@@ -32,6 +26,29 @@ 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
@@ -39,10 +56,6 @@ 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
 =====================================
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 48fe109999..5aded27c84 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -254,7 +254,7 @@ endif()
 
 # Axom/Sidre
 if (MFEM_USE_SIDRE)
-  find_package(Axom REQUIRED Sidre SLIC axom_utils)
+  find_package(Axom REQUIRED Axom)
 endif()
 
 # PUMI
diff --git a/INSTALL b/INSTALL
index 07872d4a47..0e2cfffcbf 100644
--- a/INSTALL
+++ b/INSTALL
@@ -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, 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.
+   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.
 
 MFEM_USE_CONDUIT = YES/NO
    Enables support for converting MFEM Mesh and Grid Function objects to and
@@ -543,7 +543,8 @@ 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.
-  URL: http://goo.gl/cZyJdn (axom, to be released)
+  Starting with MFEM v4.1, Axom version 0.3.1 or later is required.
+  URL: https://github.com/LLNL/axom
        https://github.com/LLNL/conduit (Conduit)
        https://support.hdfgroup.org/HDF5 (HDF5)
   Options: SIDRE_OPT, SIDRE_LIB.
@@ -714,6 +715,7 @@ MFEM_USE_PUMI
 MFEM_USE_CUDA
 MFEM_USE_OCCA
 MFEM_USE_RAJA
+MFEM_USE_SIDRE
 
 The following options are CMake specific:
 
@@ -762,6 +764,7 @@ 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:
 
diff --git a/config/cmake/modules/FindAxom.cmake b/config/cmake/modules/FindAxom.cmake
index f2cad28608..7199f67a53 100644
--- a/config/cmake/modules/FindAxom.cmake
+++ b/config/cmake/modules/FindAxom.cmake
@@ -18,6 +18,4 @@ 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 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)
+  ADD_COMPONENT Axom "include" axom/config.hpp "lib" axom)
diff --git a/config/defaults.cmake b/config/defaults.cmake
index a8ba42d635..d4525b6486 100644
--- a/config/defaults.cmake
+++ b/config/defaults.cmake
@@ -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" CACHE STRING
+set(Axom_REQUIRED_PACKAGES "Conduit/relay/blueprint" CACHE STRING
     "Additional packages required by Axom.")
 
 set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
diff --git a/config/defaults.mk b/config/defaults.mk
index e8e9fd03ce..564b1307ed 100644
--- a/config/defaults.mk
+++ b/config/defaults.mk
@@ -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 \
-   -lsidre -lslic -laxom_utils -lconduit -lconduit_relay -lhdf5 $(ZLIB_LIB) -ldl
+   -laxom -lconduit -lconduit_relay -lconduit_blueprint -lhdf5 $(ZLIB_LIB) -ldl
 
 # PUMI
 # Note that PUMI_DIR is needed -- it is used to check for gmi_sim.h
diff --git a/fem/CMakeLists.txt b/fem/CMakeLists.txt
index 59e7befcc3..2c3d1c4317 100644
--- a/fem/CMakeLists.txt
+++ b/fem/CMakeLists.txt
@@ -32,6 +32,7 @@ set(SRCS
   nonlininteg.cpp
   staticcond.cpp
   tmop.cpp
+  tmop_tools.cpp
   )
 
 set(HDRS
@@ -64,6 +65,7 @@ set(HDRS
   tfespace.hpp
   tintrules.hpp
   tmop.hpp
+  tmop_tools.hpp
   )
 
 if (MFEM_USE_SIDRE)
diff --git a/fem/bilininteg_diffusion.cpp b/fem/bilininteg_diffusion.cpp
index 8a0ae06dbe..657aac052e 100644
--- a/fem/bilininteg_diffusion.cpp
+++ b/fem/bilininteg_diffusion.cpp
@@ -27,7 +27,7 @@ static void OccaPADiffusionSetup2D(const int D1D,
                                    const int NE,
                                    const Array<double> &W,
                                    const Vector &J,
-                                   const double COEFF,
+                                   const Vector &C,
                                    Vector &op)
 {
    occa::properties props;
@@ -35,7 +35,9 @@ 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())
@@ -45,7 +47,7 @@ static void OccaPADiffusionSetup2D(const int D1D,
                                      "DiffusionSetup2D", props);
       OccaDiffSetup2D_ker.emplace(id, DiffusionSetup2D);
    }
-   OccaDiffSetup2D_ker.at(id)(NE, o_W, o_J, COEFF, o_op);
+   OccaDiffSetup2D_ker.at(id)(NE, o_W, o_J, o_C, o_op, const_c);
 }
 
 static void OccaPADiffusionSetup3D(const int D1D,
@@ -53,7 +55,7 @@ static void OccaPADiffusionSetup3D(const int D1D,
                                    const int NE,
                                    const Array<double> &W,
                                    const Vector &J,
-                                   const double COEFF,
+                                   const Vector &C,
                                    Vector &op)
 {
    occa::properties props;
@@ -61,7 +63,9 @@ 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())
@@ -71,7 +75,7 @@ static void OccaPADiffusionSetup3D(const int D1D,
                                      "DiffusionSetup3D", props);
       OccaDiffSetup3D_ker.emplace(id, DiffusionSetup3D);
    }
-   OccaDiffSetup3D_ker.at(id)(NE, o_W, o_J, COEFF, o_op);
+   OccaDiffSetup3D_ker.at(id)(NE, o_W, o_J, o_C, o_op, const_c);
 }
 #endif // MFEM_USE_OCCA
 
@@ -80,14 +84,16 @@ static void PADiffusionSetup2D(const int Q1D,
                                const int NE,
                                const Array<double> &w,
                                const Vector &j,
-                               const double COEFF,
-                               Vector &op)
+                               const Vector &c,
+                               Vector &d)
 {
    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 y = Reshape(op.Write(), NQ, 3, NE);
+   auto C = const_c ? Reshape(c.Read(), 1, 1) : Reshape(c.Read(), NQ, NE);
+   auto D = Reshape(d.Write(), NQ, 3, NE);
 
    MFEM_FORALL(e, NE,
    {
@@ -97,10 +103,11 @@ 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 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
+         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
       }
    });
 }
@@ -110,13 +117,15 @@ static void PADiffusionSetup3D(const int Q1D,
                                const int NE,
                                const Array<double> &w,
                                const Vector &j,
-                               const double COEFF,
-                               Vector &op)
+                               const Vector &c,
+                               Vector &d)
 {
    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 y = Reshape(op.Write(), NQ, 6, NE);
+   auto C = const_c ? Reshape(c.Read(), 1, 1) : Reshape(c.Read(), NQ, NE);
+   auto D = Reshape(d.Write(), NQ, 6, NE);
    MFEM_FORALL(e, NE,
    {
       for (int q = 0; q < NQ; ++q)
@@ -133,7 +142,8 @@ 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 c_detJ = W[q] * COEFF / detJ;
+         const double coeff = const_c ? C(0,0) : C(q,e);
+         const double c_detJ = W[q] * coeff / detJ;
          // adj(J)
          const double A11 = (J22 * J33) - (J23 * J32);
          const double A12 = (J32 * J13) - (J12 * J33);
@@ -145,12 +155,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
-         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
+         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
       }
    });
 }
@@ -161,8 +171,8 @@ static void PADiffusionSetup(const int dim,
                              const int NE,
                              const Array<double> &W,
                              const Vector &J,
-                             const double COEFF,
-                             Vector &op)
+                             const Vector &C,
+                             Vector &D)
 {
    if (dim == 1) { MFEM_ABORT("dim==1 not supported in PADiffusionSetup"); }
    if (dim == 2)
@@ -170,22 +180,22 @@ static void PADiffusionSetup(const int dim,
 #ifdef MFEM_USE_OCCA
       if (DeviceCanUseOcca())
       {
-         OccaPADiffusionSetup2D(D1D, Q1D, NE, W, J, COEFF, op);
+         OccaPADiffusionSetup2D(D1D, Q1D, NE, W, J, C, D);
          return;
       }
 #endif // MFEM_USE_OCCA
-      PADiffusionSetup2D(Q1D, NE, W, J, COEFF, op);
+      PADiffusionSetup2D(Q1D, NE, W, J, C, D);
    }
    if (dim == 3)
    {
 #ifdef MFEM_USE_OCCA
       if (DeviceCanUseOcca())
       {
-         OccaPADiffusionSetup3D(D1D, Q1D, NE, W, J, COEFF, op);
+         OccaPADiffusionSetup3D(D1D, Q1D, NE, W, J, C, D);
          return;
       }
 #endif // MFEM_USE_OCCA
-      PADiffusionSetup3D(Q1D, NE, W, J, COEFF, op);
+      PADiffusionSetup3D(Q1D, NE, W, J, C, D);
    }
 }
 
@@ -205,11 +215,32 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
    pa_data.SetSize(symmDims * nq * ne, Device::GetMemoryType());
-   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);
+   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);
 }
 
 #ifdef MFEM_USE_OCCA
diff --git a/fem/bilininteg_mass.cpp b/fem/bilininteg_mass.cpp
index b065f8294e..10fb8eb9f0 100644
--- a/fem/bilininteg_mass.cpp
+++ b/fem/bilininteg_mass.cpp
@@ -38,24 +38,40 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
    pa_data.SetSize(ne*nq, Device::GetMemoryType());
-   ConstantCoefficient *const_coeff = dynamic_cast<ConstantCoefficient*>(Q);
-   // TODO: other types of coefficients ...
+   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));
+         }
+      }
+   }
    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,
       {
@@ -66,25 +82,20 @@ 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);
-            v(q,e) =  w[q] * constant * detJ;
+            const double coeff = const_c ? C(0,0) : C(q,e);
+            v(q,e) =  w[q] * coeff * 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,
       {
@@ -96,7 +107,8 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
             const double detJ = J11 * (J22 * J33 - J32 * J23) -
             /* */               J21 * (J12 * J33 - J32 * J13) +
             /* */               J31 * (J12 * J23 - J22 * J13);
-            v(q,e) = W[q] * constant * detJ;
+            const double coeff = const_c ? C(0,0) : C(q,e);
+            v(q,e) = W[q] * coeff * detJ;
          }
       });
    }
diff --git a/fem/fe.cpp b/fem/fe.cpp
index 4139d8357f..776e0fccb1 100644
--- a/fem/fe.cpp
+++ b/fem/fe.cpp
@@ -9648,6 +9648,7 @@ 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++)
          {
diff --git a/fem/fem.hpp b/fem/fem.hpp
index 849dbb981a..0c9482dffb 100644
--- a/fem/fem.hpp
+++ b/fem/fem.hpp
@@ -31,6 +31,7 @@
 #include "estimators.hpp"
 #include "staticcond.hpp"
 #include "tmop.hpp"
+#include "tmop_tools.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
diff --git a/fem/gridfunc.cpp b/fem/gridfunc.cpp
index 4c061cd39d..b243b071f3 100644
--- a/fem/gridfunc.cpp
+++ b/fem/gridfunc.cpp
@@ -1716,6 +1716,7 @@ 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;
diff --git a/fem/gridfunc.hpp b/fem/gridfunc.hpp
index 84796a1f84..39f520460b 100644
--- a/fem/gridfunc.hpp
+++ b/fem/gridfunc.hpp
@@ -434,6 +434,8 @@ 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
diff --git a/fem/lininteg.cpp b/fem/lininteg.cpp
index a4d0ea72f3..ce6062a3b6 100644
--- a/fem/lininteg.cpp
+++ b/fem/lininteg.cpp
@@ -350,7 +350,6 @@ void VectorFEDomainLFIntegrator::AssembleDeltaElementVect(
    vshape.Mult(vec, elvect);
 }
 
-
 void VectorBoundaryFluxLFIntegrator::AssembleRHSElementVect(
    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {
@@ -397,19 +396,26 @@ 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);
 
-      add(elvect, val, shape, elvect);
+      double val = ip.weight;
+      if (F)
+      {
+         Tr.SetIntPoint (&ip);
+         val *= F->Eval(Tr, ip);
+      }
+
+      elvect.Add(val, shape);
    }
 }
 
diff --git a/fem/lininteg.hpp b/fem/lininteg.hpp
index 7a085f64ce..371052625a 100644
--- a/fem/lininteg.hpp
+++ b/fem/lininteg.hpp
@@ -279,11 +279,12 @@ public:
 class VectorFEBoundaryFluxLFIntegrator : public LinearFormIntegrator
 {
 private:
-   Coefficient &F;
+   Coefficient *F;
    Vector shape;
 
 public:
-   VectorFEBoundaryFluxLFIntegrator(Coefficient &f) : F(f) { }
+   VectorFEBoundaryFluxLFIntegrator() : F(NULL) { }
+   VectorFEBoundaryFluxLFIntegrator(Coefficient &f) : F(&f) { }
 
    virtual void AssembleRHSElementVect(const FiniteElement &el,
                                        ElementTransformation &Tr,
diff --git a/fem/occa.okl b/fem/occa.okl
index c9ffd296b8..7c887ff0dc 100644
--- a/fem/occa.okl
+++ b/fem/occa.okl
@@ -38,19 +38,24 @@ 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,
-                              const double COEFF,
-                              @restrict SymmOperator2D_t op) {
+                              @restrict const Coeff2D_t C,
+                              @restrict SymmOperator2D_t op,
+                              const bool const_c) {
   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 c_detJ = W[q] * COEFF / ((J11 * J22) - (J21 * J12));
+      const double coeff = const_c ? C(0,0) : C(q,e);
+      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)
@@ -61,8 +66,9 @@ typedef double* SymmOperator3D_t @dim(Q3D, 6, NE);
 @kernel void DiffusionSetup3D(const int NE,
                               @restrict const double *W,
                               @restrict const Jacobian3D_t J,
-                              const double COEFF,
-                              @restrict SymmOperator3D_t op) {
+                              @restrict const Coeff3D_t C,
+                              @restrict SymmOperator3D_t op,
+                              const bool const_c) {
   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);
@@ -72,7 +78,8 @@ 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 c_detJ = W[q] * COEFF / detJ;
+      const double coeff = const_c ? C(0,0) : C(q,e);
+      const double c_detJ = W[q] * coeff / detJ;
 
       // adj(J)
       const double A11 = (J22 * J33) - (J23 * J32);
diff --git a/fem/pgridfunc.cpp b/fem/pgridfunc.cpp
index 329967f5ef..72c2587efa 100644
--- a/fem/pgridfunc.cpp
+++ b/fem/pgridfunc.cpp
@@ -225,11 +225,13 @@ 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] = data[send_ldof[i]];
+      send_data[i] = h_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);
@@ -239,7 +241,7 @@ void ParGridFunction::ExchangeFaceNbrData()
                 send_offset[fn+1] - send_offset[fn],
                 MPI_DOUBLE, nbr_rank, tag, MyComm, &send_requests[fn]);
 
-      MPI_Irecv(&face_nbr_data(recv_offset[fn]),
+      MPI_Irecv(&h_face_nbr_data[recv_offset[fn]],
                 recv_offset[fn+1] - recv_offset[fn],
                 MPI_DOUBLE, nbr_rank, tag, MyComm, &recv_requests[fn]);
    }
diff --git a/fem/pgridfunc.hpp b/fem/pgridfunc.hpp
index 62a4f34531..9c2c12cb38 100644
--- a/fem/pgridfunc.hpp
+++ b/fem/pgridfunc.hpp
@@ -112,6 +112,8 @@ 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
diff --git a/fem/pnonlinearform.cpp b/fem/pnonlinearform.cpp
index 326606f901..123335c59c 100644
--- a/fem/pnonlinearform.cpp
+++ b/fem/pnonlinearform.cpp
@@ -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.SetData(aux2.GetData()); // aux2 contains A_local.P.x
+   Y.MakeRef(aux2, 0); // 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.SetData(aux1.GetData()); // aux1 contains P.x
+      X.MakeRef(aux1, 0); // aux1 contains P.x
       X.ExchangeFaceNbrData();
       const int n_shared_faces = pmesh->GetNSharedFaces();
       for (int i = 0; i < n_shared_faces; i++)
diff --git a/fem/sidredatacollection.cpp b/fem/sidredatacollection.cpp
index b169f82873..a280ef9a50 100644
--- a/fem/sidredatacollection.cpp
+++ b/fem/sidredatacollection.cpp
@@ -16,9 +16,7 @@
 
 #include "fem.hpp"
 
-#ifdef MFEM_USE_MPI
-#include <sidre/IOManager.hpp>
-#endif
+#include <axom/sidre.hpp>
 
 #include <string>
 #include <iomanip>      // for setw, setfill
@@ -204,10 +202,10 @@ SidreDataCollection::get_file_path(const std::string &filename) const
 
 axom::sidre::View *
 SidreDataCollection::AllocNamedBuffer(const std::string& buffer_name,
-                                      axom::sidre::SidreLength sz,
+                                      axom::sidre::IndexType sz,
                                       axom::sidre::TypeID type)
 {
-   sz = std::max(sz, sidre::SidreLength(0));
+   sz = std::max(sz, sidre::IndexType(0));
    sidre::Group *f = named_buffers_grp();
    sidre::View  *v = NULL;
 
@@ -825,7 +823,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::SidreLength offset)
+                           axom::sidre::IndexType offset)
 {
    sidre::Group* grp = m_bp_grp->getGroup("fields/" + field_name);
    MFEM_ASSERT(grp != NULL, "field " << field_name << " does not exist");
@@ -888,7 +886,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::SidreLength offset)
+                           axom::sidre::IndexType offset)
 {
    sidre::Group* grp = m_bp_grp->getGroup("fields/" + field_name);
    MFEM_ASSERT(grp != NULL, "field " << field_name << " does not exist");
@@ -1013,7 +1011,7 @@ DeregisterFieldInBPIndex(const std::string& field_name)
 void SidreDataCollection::RegisterField(const std::string &field_name,
                                         GridFunction *gf,
                                         const std::string &buffer_name,
-                                        axom::sidre::SidreLength offset)
+                                        axom::sidre::IndexType offset)
 {
    if ( field_name.empty() || buffer_name.empty() ||
         gf == NULL || gf->FESpace() == NULL )
diff --git a/fem/sidredatacollection.hpp b/fem/sidredatacollection.hpp
index 617257dda2..d2e2242748 100644
--- a/fem/sidredatacollection.hpp
+++ b/fem/sidredatacollection.hpp
@@ -25,7 +25,7 @@
 #  pragma GCC diagnostic ignored "-Wpedantic"
 # endif
 #endif
-#include <sidre/sidre.hpp>
+#include <axom/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::SidreLength offset);
+                      axom::sidre::IndexType 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::SidreLength sz,
+                    axom::sidre::IndexType 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::SidreLength offset);
+                                   axom::sidre::IndexType 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::SidreLength offset);
+                                   axom::sidre::IndexType offset);
 
    /** @brief A private helper function to set up the Views associated with
        attribute field named @a field_name */
diff --git a/fem/tmop.cpp b/fem/tmop.cpp
index 0a08b915a4..10a093f421 100644
--- a/fem/tmop.cpp
+++ b/fem/tmop.cpp
@@ -12,6 +12,7 @@
 #include "tmop.hpp"
 #include "linearform.hpp"
 #include "pgridfunc.hpp"
+#include "tmop_tools.hpp"
 
 namespace mfem
 {
@@ -768,7 +769,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();
@@ -787,9 +788,13 @@ void TargetConstructor::ComputeAvgVolume() const
          volume += ip.weight * Tr.Weight();
       }
    }
-   if (!Parallel())
+
+   NCMesh *ncmesh = mesh->ncmesh;
+   if (Parallel() == false)
    {
-      avg_volume = volume / NE;
+      avg_volume = (ncmesh == NULL) ?
+                   volume / NE : volume / ncmesh->GetNumRootElements();
+
    }
 #ifdef MFEM_USE_MPI
    else
@@ -797,7 +802,8 @@ 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 = area_NE[2] / area_NE[3];
+      avg_volume = (ncmesh == NULL) ?
+                   area_NE[2] / area_NE[3] : area_NE[2] / ncmesh->GetNumRootElements();
    }
 #endif
 }
@@ -805,6 +811,7 @@ 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, "");
@@ -827,7 +834,15 @@ void TargetConstructor::ComputeElementTargets(int e_id, const FiniteElement &fe,
       {
          if (avg_volume == 0.0) { ComputeAvgVolume(); }
          DenseMatrix W(Wideal.Height());
-         W.Set(std::pow(volume_scale * avg_volume / Wideal.Det(),
+
+         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(),
                         1./W.Height()), Wideal);
          for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = W; }
          break;
@@ -853,7 +868,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, "Initial mesh is inverted!");
+               MFEM_VERIFY(det > 0.0, "The given mesh is inverted!");
                Jtr(i).Set(std::pow(det / detW, 1./dim), Wideal);
             }
          }
@@ -864,6 +879,162 @@ 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)
@@ -921,7 +1092,7 @@ double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
 
    energy = 0.0;
    DenseTensor Jtr(dim, dim, ir->GetNPoints());
-   targetC->ComputeElementTargets(T.ElementNo, el, *ir, Jtr);
+   targetC->ComputeElementTargets(T.ElementNo, el, *ir, elfun, Jtr);
 
    // Limited case.
    Vector shape, p, p0, d_vals;
@@ -990,6 +1161,7 @@ double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
       energy += weight * val;
    }
    delete Tpr;
+
    return energy;
 }
 
@@ -1016,7 +1188,7 @@ void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
 
    elvect = 0.0;
    DenseTensor Jtr(dim, dim, ir->GetNPoints());
-   targetC->ComputeElementTargets(T.ElementNo, el, *ir, Jtr);
+   targetC->ComputeElementTargets(T.ElementNo, el, *ir, elfun, Jtr);
 
    // Limited case.
    DenseMatrix pos0;
@@ -1072,6 +1244,8 @@ 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);
@@ -1107,7 +1281,7 @@ void TMOP_Integrator::AssembleElementGrad(const FiniteElement &el,
 
    elmat = 0.0;
    DenseTensor Jtr(dim, dim, ir->GetNPoints());
-   targetC->ComputeElementTargets(T.ElementNo, el, *ir, Jtr);
+   targetC->ComputeElementTargets(T.ElementNo, el, *ir, elfun, Jtr);
 
    // Limited case.
    DenseMatrix pos0, grad_grad;
@@ -1160,6 +1334,8 @@ 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);
@@ -1234,11 +1410,12 @@ 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);
@@ -1274,9 +1451,6 @@ 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);
@@ -1285,6 +1459,9 @@ 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);
diff --git a/fem/tmop.hpp b/fem/tmop.hpp
index 5bc9470390..6190aaacc3 100644
--- a/fem/tmop.hpp
+++ b/fem/tmop.hpp
@@ -12,7 +12,6 @@
 #ifndef MFEM_TMOP_HPP
 #define MFEM_TMOP_HPP
 
-#include "../config/config.hpp"
 #include "../linalg/invariants.hpp"
 #include "nonlininteg.hpp"
 
@@ -514,6 +513,51 @@ 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). */
@@ -538,9 +582,11 @@ 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:
@@ -589,14 +635,89 @@ 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. */
+       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 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.
 
diff --git a/fem/tmop_tools.cpp b/fem/tmop_tools.cpp
new file mode 100644
index 0000000000..f67310f334
--- /dev/null
+++ b/fem/tmop_tools.cpp
@@ -0,0 +1,518 @@
+// 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";
+}
+
+}
diff --git a/fem/tmop_tools.hpp b/fem/tmop_tools.hpp
new file mode 100644
index 0000000000..bded562657
--- /dev/null
+++ b/fem/tmop_tools.hpp
@@ -0,0 +1,140 @@
+// 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
diff --git a/general/device.hpp b/general/device.hpp
index ce77bbee77..cd6cc1b706 100644
--- a/general/device.hpp
+++ b/general/device.hpp
@@ -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 const T *HostWrite(const Memory<T> &mem, int size)
+inline T *HostWrite(Memory<T> &mem, int size)
 {
    return mfem::Write(mem, size, false);
 }
diff --git a/general/mem_manager.cpp b/general/mem_manager.cpp
index 67114194fc..b45eea17f4 100644
--- a/general/mem_manager.cpp
+++ b/general/mem_manager.cpp
@@ -11,7 +11,7 @@
 
 #include "../general/forall.hpp"
 
-#include <cstring> // std::memcpy
+#include <cstring> // std::memcpy, std::memcmp
 
 #include <list>
 #include <unordered_map>
@@ -57,6 +57,15 @@ 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
 {
 
@@ -695,6 +704,19 @@ 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)
 {
diff --git a/general/mem_manager.hpp b/general/mem_manager.hpp
index 26ab2d46d9..5c4744b78b 100644
--- a/general/mem_manager.hpp
+++ b/general/mem_manager.hpp
@@ -377,6 +377,16 @@ 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;
 };
 
 
@@ -447,6 +457,11 @@ 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);
 
@@ -727,10 +742,24 @@ 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. */
+    Useful in a debugger. See also Memory<T>::PrintFlags(). */
 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;
 
diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp
index 27c367f450..bd5a8adde9 100644
--- a/linalg/densemat.cpp
+++ b/linalg/densemat.cpp
@@ -4456,7 +4456,8 @@ 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, x_col;
+   double *yp = y.HostReadWrite();
+   double x_col;
    const double *xp = x;
    // the '4' here can be tuned for given platform and compiler
    if (n <= 4)
diff --git a/linalg/hypre.cpp b/linalg/hypre.cpp
index 4153381e5a..c9d3d50e4f 100644
--- a/linalg/hypre.cpp
+++ b/linalg/hypre.cpp
@@ -185,7 +185,7 @@ Vector * HypreParVector::GlobalVector() const
 
 HypreParVector& HypreParVector::operator=(double d)
 {
-   hypre_ParVectorSetConstantValues(x,d);
+   Vector::operator=(d);
    return *this;
 }
 
@@ -198,10 +198,7 @@ HypreParVector& HypreParVector::operator=(const HypreParVector &y)
    }
 #endif
 
-   for (int i = 0; i < size; i++)
-   {
-      data[i] = y.data[i];
-   }
+   Vector::operator=(y);
    return *this;
 }
 
@@ -1024,6 +1021,8 @@ 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);
 }
 
@@ -1035,7 +1034,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 = y.HostWrite();
+   auto y_data = (b == 0.0) ? y.HostWrite() : y.HostReadWrite();
    if (X == NULL)
    {
       X = new HypreParVector(A->comm,
@@ -1067,7 +1066,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 = y.HostWrite();
+   auto y_data = (b == 0.0) ? y.HostWrite() : y.HostReadWrite();
    if (X == NULL)
    {
       X = new HypreParVector(A->comm,
@@ -2009,10 +2008,12 @@ 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)
          {
@@ -2022,6 +2023,7 @@ void HypreSmoother::Mult(const HypreParVector &b, HypreParVector &x) const
       }
       x = 0.0;
    }
+   x.HostReadWrite();
 
    if (V == NULL)
    {
@@ -2071,21 +2073,25 @@ 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(),
-                             b.GetData(),
+                             const_cast<double*>(b_data),
                              A -> GetRowStarts());
       X = new HypreParVector(A->GetComm(),
                              A -> GetGlobalNumCols(),
-                             x.GetData(),
+                             x_data,
                              A -> GetColStarts());
    }
    else
    {
-      B -> SetData(b.GetData());
-      X -> SetData(x.GetData());
+      B -> SetData(const_cast<double*>(b_data));
+      X -> SetData(x_data);
    }
 
    Mult(*B, *X);
diff --git a/linalg/petsc.cpp b/linalg/petsc.cpp
index f76e2fd37b..429d81938a 100644
--- a/linalg/petsc.cpp
+++ b/linalg/petsc.cpp
@@ -2309,17 +2309,6 @@ 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)
@@ -2342,6 +2331,7 @@ 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)
    {
@@ -2362,15 +2352,7 @@ void PetscLinearSolver::SetOperator(const Operator &op)
          wrap = false;
       }
    }
-   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);
-   }
+   ierr = KSPSetOperators(ksp,A,A); PCHKERRQ(ksp,ierr);
 
    // Update PetscSolver
    operatorset = true;
@@ -3229,26 +3211,27 @@ 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;
-   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);
+   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);
+   }
 
    // 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++)
@@ -3616,6 +3599,7 @@ 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,
@@ -3771,17 +3755,37 @@ 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);
@@ -3938,6 +3942,7 @@ 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);
@@ -4024,12 +4029,31 @@ 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);
@@ -4133,10 +4157,30 @@ 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);
diff --git a/linalg/petsc.hpp b/linalg/petsc.hpp
index 090348b033..d97a152ad4 100644
--- a/linalg/petsc.hpp
+++ b/linalg/petsc.hpp
@@ -618,6 +618,8 @@ 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
@@ -625,6 +627,7 @@ 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.
diff --git a/linalg/solvers.cpp b/linalg/solvers.cpp
index 6e240ea9dc..58e0dbaf49 100644
--- a/linalg/solvers.cpp
+++ b/linalg/solvers.cpp
@@ -1353,6 +1353,8 @@ void NewtonSolver::Mult(const Vector &b, Vector &x) const
       }
       add(x, -c_scale, c, x);
 
+      ProcessNewState(x);
+
       oper->Mult(x, r);
       if (have_b)
       {
diff --git a/linalg/solvers.hpp b/linalg/solvers.hpp
index 5de637e27d..f960590731 100644
--- a/linalg/solvers.hpp
+++ b/linalg/solvers.hpp
@@ -319,6 +319,10 @@ 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.
diff --git a/linalg/sparsemat.cpp b/linalg/sparsemat.cpp
index 9ea95f2c26..de24c9cba9 100644
--- a/linalg/sparsemat.cpp
+++ b/linalg/sparsemat.cpp
@@ -2640,9 +2640,11 @@ SparseMatrix &SparseMatrix::operator=(double a)
 {
    if (Rows == NULL)
    {
-      for (int i = 0, nnz = I[height]; i < nnz; i++)
+      const int nnz = J.Capacity();
+      double *h_A = HostWrite(A, nnz);
+      for (int i = 0; i < nnz; i++)
       {
-         A[i] = a;
+         h_A[i] = a;
       }
    }
    else
diff --git a/linalg/sparsemat.hpp b/linalg/sparsemat.hpp
index 75da18a332..c6bf5d5da0 100644
--- a/linalg/sparsemat.hpp
+++ b/linalg/sparsemat.hpp
@@ -152,6 +152,54 @@ 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;
 
diff --git a/linalg/superlu.cpp b/linalg/superlu.cpp
index fdb8d1e648..2089f88bb9 100644
--- a/linalg/superlu.cpp
+++ b/linalg/superlu.cpp
@@ -134,6 +134,11 @@ 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;
diff --git a/linalg/vector.cpp b/linalg/vector.cpp
index 30886b4761..492a2f44a6 100644
--- a/linalg/vector.cpp
+++ b/linalg/vector.cpp
@@ -834,9 +834,10 @@ double Vector::Sum() const
 {
    double sum = 0.0;
 
+   const double *h_data = this->HostRead();
    for (int i = 0; i < size; i++)
    {
-      sum += data[i];
+      sum += h_data[i];
    }
 
    return sum;
diff --git a/linalg/vector.hpp b/linalg/vector.hpp
index 64abb86c06..334bc85218 100644
--- a/linalg/vector.hpp
+++ b/linalg/vector.hpp
@@ -140,6 +140,13 @@ 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); }
 
@@ -455,6 +462,19 @@ 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();
diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
index 6e7bce2e51..b8ed54a376 100644
--- a/mesh/mesh.cpp
+++ b/mesh/mesh.cpp
@@ -430,6 +430,8 @@ 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;
@@ -440,7 +442,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);
@@ -492,6 +494,8 @@ 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;
@@ -502,7 +506,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);
diff --git a/mesh/ncmesh.cpp b/mesh/ncmesh.cpp
index b48fc0c45d..19c4e8059b 100644
--- a/mesh/ncmesh.cpp
+++ b/mesh/ncmesh.cpp
@@ -4307,6 +4307,20 @@ 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
diff --git a/mesh/ncmesh.hpp b/mesh/ncmesh.hpp
index 1fa1b238de..9a15692d9b 100644
--- a/mesh/ncmesh.hpp
+++ b/mesh/ncmesh.hpp
@@ -319,13 +319,21 @@ 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;
 
diff --git a/mesh/pmesh.cpp b/mesh/pmesh.cpp
index 06ba4e0ed5..bc97df2aac 100644
--- a/mesh/pmesh.cpp
+++ b/mesh/pmesh.cpp
@@ -4171,6 +4171,7 @@ void ParMesh::PrintAsOne(std::ostream &out)
           "# SQUARE      = 3\n"
           "# TETRAHEDRON = 4\n"
           "# CUBE        = 5\n"
+          "# PRISM       = 6\n"
           "#\n";
 
       out << "\ndimension\n" << Dim;
@@ -4269,6 +4270,15 @@ 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++)
             {
diff --git a/miniapps/electromagnetics/volta_solver.cpp b/miniapps/electromagnetics/volta_solver.cpp
index f5f477cdd8..d6eabe57d7 100644
--- a/miniapps/electromagnetics/volta_solver.cpp
+++ b/miniapps/electromagnetics/volta_solver.cpp
@@ -41,6 +41,7 @@ VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
      H1FESpace_(NULL),
      HCurlFESpace_(NULL),
      HDivFESpace_(NULL),
+     L2FESpace_(NULL),
      divEpsGrad_(NULL),
      h1Mass_(NULL),
      h1SurfMass_(NULL),
@@ -49,20 +50,24 @@ 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_(NULL),
+     sigma_src_(NULL),
      e_(NULL),
      d_(NULL),
-     p_(NULL),
+     p_src_(NULL),
+     oneCoef_(1.0),
      epsCoef_(&epsCoef),
      phiBCCoef_(NULL),
      rhoCoef_(NULL),
      pCoef_(NULL),
-     phi_bc_(phi_bc),
-     rho_src_(rho_src),
-     p_src_(p_src),
+     phi_bc_func_(phi_bc),
+     rho_src_func_(rho_src),
+     p_src_func_(p_src),
      point_charge_params_(point_charges),
      point_charges_(0)
 {
@@ -76,34 +81,38 @@ 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++)
    {
-      ess_bdr_[(*dbcs_)[i]-1] = 1;
+      if ((*dbcs_)[i] <= ess_bdr_.Size())
+      {
+         ess_bdr_[(*dbcs_)[i]-1] = 1;
+      }
    }
 
    // Setup various coefficients
 
    // Potential on outer surface
-   if ( phi_bc_ != NULL )
+   if ( phi_bc_func_ != NULL )
    {
-      phiBCCoef_ = new FunctionCoefficient(*phi_bc_);
+      phiBCCoef_ = new FunctionCoefficient(*phi_bc_func_);
    }
 
    // Volume Charge Density
-   if ( rho_src_ != NULL )
+   if ( rho_src_func_ != NULL )
    {
-      rhoCoef_ = new FunctionCoefficient(rho_src_);
+      rhoCoef_ = new FunctionCoefficient(rho_src_func_);
    }
 
    // Polarization
-   if ( p_src_ != NULL )
+   if ( p_src_func_ != NULL )
    {
       pCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
-                                             p_src_);
+                                             p_src_func_);
    }
 
    // Bilinear Forms
@@ -116,15 +125,23 @@ VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
    hCurlHDivEps_ = new ParMixedBilinearForm(HCurlFESpace_,HDivFESpace_);
    hCurlHDivEps_->AddDomainIntegrator(new VectorFEMassIntegrator(*epsCoef_));
 
-   rhod_ = new ParLinearForm(H1FESpace_);
+   rhod_   = new ParLinearForm(H1FESpace_);
 
-   // Discrete Grad operator
+   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
    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 )
    {
@@ -149,29 +166,28 @@ VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
       }
    }
 
-   if ( rho_src_ )
+   if ( rho_src_func_ )
    {
-      rho_ = new ParGridFunction(H1FESpace_);
+      rho_src_ = new ParGridFunction(H1FESpace_);
 
       h1Mass_ = new ParBilinearForm(H1FESpace_);
       h1Mass_->AddDomainIntegrator(new MassIntegrator);
    }
 
-   if ( p_src_ )
+   if ( p_src_func_ )
    {
-      p_ = new ParGridFunction(HCurlFESpace_);
+      p_src_ = 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_ = new ParGridFunction(H1FESpace_);
+      sigma_src_ = new ParGridFunction(H1FESpace_);
 
       h1SurfMass_ = new ParBilinearForm(H1FESpace_);
       h1SurfMass_->AddBoundaryIntegrator(new MassIntegrator);
@@ -185,14 +201,18 @@ VoltaSolver::~VoltaSolver()
    delete pCoef_;
 
    delete phi_;
+   delete rho_src_;
    delete rho_;
    delete rhod_;
-   delete sigma_;
+   delete l2_vol_int_;
+   delete rt_surf_int_;
+   delete sigma_src_;
    delete d_;
    delete e_;
-   delete p_;
+   delete p_src_;
 
    delete grad_;
+   delete div_;
 
    delete divEpsGrad_;
    delete h1Mass_;
@@ -205,6 +225,7 @@ VoltaSolver::~VoltaSolver()
    delete H1FESpace_;
    delete HCurlFESpace_;
    delete HDivFESpace_;
+   delete L2FESpace_;
 
    for (unsigned int i=0; i<point_charges_.size(); i++)
    {
@@ -230,11 +251,13 @@ 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;
    }
 }
 
@@ -254,9 +277,15 @@ 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();
@@ -292,15 +321,19 @@ 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();
-   if ( rho_   ) { rho_->Update(); }
-   if ( sigma_ ) { sigma_->Update(); }
-   if ( p_     ) { p_->Update(); }
+   rho_->Update();
+   if ( rho_src_   ) { rho_src_->Update(); }
+   if ( sigma_src_ ) { sigma_src_->Update(); }
+   if ( p_src_     ) { p_src_->Update(); }
 
    // Inform the bilinear forms that the space has changed.
    divEpsGrad_->Update();
@@ -314,6 +347,7 @@ VoltaSolver::Update()
 
    // Inform the other objects that the space has changed.
    grad_->Update();
+   div_->Update();
 }
 
 void
@@ -339,40 +373,46 @@ VoltaSolver::Solve()
          {
             ConstantCoefficient voltage((*dbcv_)[i]);
             dbc_bdr_attr = 0;
-            dbc_bdr_attr[(*dbcs_)[i]-1] = 1;
+            if ((*dbcs_)[i] <= dbc_bdr_attr.Size())
+            {
+               dbc_bdr_attr[(*dbcs_)[i]-1] = 1;
+            }
             phi_->ProjectBdrCoefficient(voltage, dbc_bdr_attr);
          }
       }
    }
 
    // Initialize the volumetric charge density
-   if ( rho_ )
+   if ( rho_src_ )
    {
-      rho_->ProjectCoefficient(*rhoCoef_);
-      h1Mass_->AddMult(*rho_, *rhod_);
+      rho_src_->ProjectCoefficient(*rhoCoef_);
+      h1Mass_->AddMult(*rho_src_, *rhod_);
    }
 
    // Initialize the Polarization
-   if ( p_ )
+   if ( p_src_ )
    {
-      p_->ProjectCoefficient(*pCoef_);
-      weakDiv_->AddMult(*p_, *rhod_);
+      p_src_->ProjectCoefficient(*pCoef_);
+      weakDiv_->AddMult(*p_src_, *rhod_);
    }
 
    // Initialize the surface charge density
-   if ( sigma_ )
+   if ( sigma_src_ )
    {
-      *sigma_ = 0.0;
+      *sigma_src_ = 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;
-         nbc_bdr_attr[(*nbcs_)[i]-1] = 1;
-         sigma_->ProjectBdrCoefficient(sigma_coef, nbc_bdr_attr);
+         if ((*nbcs_)[i] <= nbc_bdr_attr.Size())
+         {
+            nbc_bdr_attr[(*nbcs_)[i]-1] = 1;
+         }
+         sigma_src_->ProjectBdrCoefficient(sigma_coef, nbc_bdr_attr);
       }
-      h1SurfMass_->AddMult(*sigma_, *rhod_);
+      h1SurfMass_->AddMult(*sigma_src_, *rhod_);
    }
 
    // Determine the essential BC degrees of freedom
@@ -424,9 +464,9 @@ VoltaSolver::Solve()
 
    ParGridFunction ed(HDivFESpace_);
    hCurlHDivEps_->Mult(*e_, ed);
-   if ( p_ )
+   if ( p_src_ )
    {
-      hCurlHDiv_->AddMult(*p_, ed, -1.0);
+      hCurlHDiv_->AddMult(*p_src_, ed, -1.0);
    }
 
    HypreParMatrix MassHDiv;
@@ -445,8 +485,27 @@ 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; }
 }
 
@@ -480,9 +539,10 @@ VoltaSolver::RegisterVisItFields(VisItDataCollection & visit_dc)
    visit_dc.RegisterField("Phi", phi_);
    visit_dc.RegisterField("D",     d_);
    visit_dc.RegisterField("E",     e_);
-   if ( rho_   ) { visit_dc.RegisterField("Rho",     rho_); }
-   if ( p_     ) { visit_dc.RegisterField("P",         p_); }
-   if ( sigma_ ) { visit_dc.RegisterField("Sigma", sigma_); }
+   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_); }
 }
 
 void
@@ -515,20 +575,23 @@ VoltaSolver::InitializeGLVis()
    socks_["E"] = new socketstream;
    socks_["E"]->precision(8);
 
-   if ( rho_)
+   socks_["Rho"] = new socketstream;
+   socks_["Rho"]->precision(8);
+
+   if ( rho_src_ )
    {
-      socks_["Rho"] = new socketstream;
-      socks_["Rho"]->precision(8);
+      socks_["RhoSrc"] = new socketstream;
+      socks_["RhoSrc"]->precision(8);
    }
-   if ( p_)
+   if ( p_src_ )
    {
-      socks_["P"] = new socketstream;
-      socks_["P"]->precision(8);
+      socks_["PSrc"] = new socketstream;
+      socks_["PSrc"]->precision(8);
    }
-   if ( sigma_)
+   if ( sigma_src_ )
    {
-      socks_["Sigma"] = new socketstream;
-      socks_["Sigma"]->precision(8);
+      socks_["SigmaSrc"] = new socketstream;
+      socks_["SigmaSrc"]->precision(8);
    }
 }
 
@@ -548,31 +611,36 @@ 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_["E"], vishost, visport,
-                  *e_, "Electric Field (E)", Wx, Wy, Ww, Wh);
-
+   VisualizeField(*socks_["Rho"], vishost, visport,
+                  *rho_, "Charge Density", Wx, Wy, Ww, Wh);
    Wx = 0; Wy += offy; // next line
 
-   if ( rho_ )
+   if ( rho_src_ )
    {
-      VisualizeField(*socks_["Rho"], vishost, visport,
-                     *rho_, "Charge Density (Rho)", Wx, Wy, Ww, Wh);
+      VisualizeField(*socks_["RhoSrc"], vishost, visport,
+                     *rho_src_, "Charge Density Source (Rho)", Wx, Wy, Ww, Wh);
       Wx += offx;
    }
-   if ( p_ )
+   if ( p_src_ )
    {
-      VisualizeField(*socks_["P"], vishost, visport,
-                     *p_, "Electric Polarization (P)", Wx, Wy, Ww, Wh);
+      VisualizeField(*socks_["PSrc"], vishost, visport,
+                     *p_src_, "Electric Polarization Source (P)",
+                     Wx, Wy, Ww, Wh);
       Wx += offx;
    }
-   if ( sigma_ )
+   if ( sigma_src_ )
    {
-      VisualizeField(*socks_["Sigma"], vishost, visport,
-                     *sigma_, "Surface Charge Density (Sigma)", Wx, Wy, Ww, Wh);
+      VisualizeField(*socks_["SigmaSrc"], vishost, visport,
+                     *sigma_src_, "Surface Charge Density Source (Sigma)",
+                     Wx, Wy, Ww, Wh);
       // Wx += offx; // not used
    }
    if (myid_ == 0) { cout << " done." << endl; }
diff --git a/miniapps/electromagnetics/volta_solver.hpp b/miniapps/electromagnetics/volta_solver.hpp
index 348b8785dc..1b5a575504 100644
--- a/miniapps/electromagnetics/volta_solver.hpp
+++ b/miniapps/electromagnetics/volta_solver.hpp
@@ -26,7 +26,9 @@ 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
 {
@@ -84,6 +86,7 @@ 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
@@ -94,26 +97,32 @@ 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
+   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
 
    ParDiscreteGradOperator * grad_; // For Computing E from phi
+   ParDiscreteDivOperator  * div_;  // For Computing rho from D
 
-   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
+   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
 
+   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_ )(const Vector&);          // Scalar Potential BC
-   double (*rho_src_)(const Vector&);          // Volumetric Charge Density
-   void   (*p_src_  )(const Vector&, Vector&); // Polarization Field
+   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
 
    const Vector & point_charge_params_;
 
diff --git a/miniapps/meshing/amr-quad-q2.mesh b/miniapps/meshing/amr-quad-q2.mesh
new file mode 100644
index 0000000000..75101310f5
--- /dev/null
+++ b/miniapps/meshing/amr-quad-q2.mesh
@@ -0,0 +1,418 @@
+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
diff --git a/miniapps/meshing/mesh-optimizer.cpp b/miniapps/meshing/mesh-optimizer.cpp
index 56df9ea9ed..cfb71b5977 100644
--- a/miniapps/meshing/mesh-optimizer.cpp
+++ b/miniapps/meshing/mesh-optimizer.cpp
@@ -32,6 +32,11 @@
 // 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:
@@ -48,6 +53,9 @@
 //     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>
@@ -58,190 +66,160 @@ using namespace std;
 
 double weight_fun(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)
+double ind_values(const Vector &x)
 {
-   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;
+   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;
 }
 
-class RelaxedNewtonSolver : public NewtonSolver
+class HessianCoefficient : public MatrixCoefficient
 {
 private:
-   // Quadrature points that are checked for negative Jacobians etc.
-   const IntegrationRule &ir;
-   FiniteElementSpace *fes;
-   mutable GridFunction x_gf;
+   int type;
 
 public:
-   RelaxedNewtonSolver(const IntegrationRule &irule, FiniteElementSpace *f)
-      : ir(irule), fes(f) { }
+   HessianCoefficient(int dim, int type_)
+      : MatrixCoefficient(dim), type(type_) { }
 
-   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      Vector pos(3);
+      T.Transform(ip, pos);
+
+      if (type == 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;
+      }
+      else
+      {
+         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;
+
+         const double tan1 = std::tanh(sf*(r-r1)),
+                      tan2 = std::tanh(sf*(r-r2));
+
+         K(0, 0) = eps + 1.0 * (tan1 - tan2);
+         K(0, 1) = 0.0;
+         K(1, 0) = 0.0;
+         K(1, 1) = 1.0;
+      }
+   }
 };
 
-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++)
-   {
-      add(x, -scale, c, x_out);
-      x_gf.SetFromTrueVector();
-
-      energy_out = nlf->GetGridFunctionEnergy(x_gf);
-      if (energy_out > 1.2*energy_in || std::isnan(energy_out) != 0)
-      {
-         if (print_level >= 0)
-         { cout << "Scale = " << scale << " Increasing energy." << endl; }
-         scale *= 0.5; continue;
-      }
-
-      int jac_ok = 1;
-      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);
-            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;
-      }
-
-      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;
-   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);
 IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform);
@@ -304,7 +282,9 @@ 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");
+                  "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)");
    args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
    args.AddOption(&quad_type, "-qt", "--quad-type",
                   "Quadrature rule type:\n\t"
@@ -375,7 +355,8 @@ 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 = mesh->GetNodes();
+   GridFunction x(fespace);
+   mesh->SetNodalGridFunction(&x);
 
    // 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
@@ -425,10 +406,9 @@ 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;
-   // Set the perturbation of all nodes from the true nodes.
-   x->SetTrueVector();
-   x->SetFromTrueVector();
+   x -= rdm;
+   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".
@@ -439,7 +419,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;
@@ -468,15 +448,43 @@ 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;
-      default: cout << "Unknown target_id: " << target_id << endl;
-         delete metric; return 3;
+      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);
    }
-   TargetConstructor *target_c = new TargetConstructor(target_t);
    target_c->SetNodes(x0);
    TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c);
 
@@ -542,13 +550,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_metric(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
+      vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
    }
 
    // 16. Fix all boundary nodes, or fix only a given component depending on the
@@ -652,8 +660,13 @@ int main (int argc, char *argv[])
    if (tauval > 0.0)
    {
       tauval = 0.0;
-      newton = new RelaxedNewtonSolver(*ir, fespace);
-      cout << "The RelaxedNewtonSolver is used (as all det(J)>0)." << endl;
+      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";
    }
    else
    {
@@ -664,8 +677,8 @@ int main (int argc, char *argv[])
          return 3;
       }
       tauval -= 0.01 * h0.Min(); // Slightly below minJ0 to avoid div by 0.
-      newton = new DescentNewtonSolver(*ir, fespace);
-      cout << "The DescentNewtonSolver is used (as some det(J)<0)." << endl;
+      newton = new TMOPDescentNewtonSolver(*ir);
+      cout << "The TMOPDescentNewtonSolver is used (as some det(J) < 0).\n";
    }
    newton->SetPreconditioner(*S);
    newton->SetMaxIter(newton_iter);
@@ -673,8 +686,9 @@ 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."
@@ -691,12 +705,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
@@ -712,13 +726,13 @@ int main (int argc, char *argv[])
    if (visualization)
    {
       char title[] = "Final metric values";
-      vis_metric(mesh_poly_deg, *metric, *target_c, *mesh, title, 600);
+      vis_tmop_metric_s(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);
diff --git a/miniapps/meshing/pmesh-optimizer.cpp b/miniapps/meshing/pmesh-optimizer.cpp
index a8b3ffed2d..3abb7b34a8 100644
--- a/miniapps/meshing/pmesh-optimizer.cpp
+++ b/miniapps/meshing/pmesh-optimizer.cpp
@@ -32,6 +32,11 @@
 // 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:
@@ -48,223 +53,176 @@
 //     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 <fstream>
 #include <iostream>
+#include <fstream>
 
 using namespace mfem;
 using namespace std;
 
 double weight_fun(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)
+double ind_values(const Vector &x)
 {
-   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 int opt = 6;
+   const double small = 0.001, big = 0.01;
+
+   // Sine wave.
+   if (opt==1)
    {
-      sock.open("localhost", 19916);
-      sock << "solution\n";
+      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;
    }
-   pmesh.PrintAsOne(sock);
-   metric.SaveAsOne(sock);
-   if (pmesh.GetMyRank() == 0)
+
+   if (opt==2)
    {
-      sock << "window_title '"<< title << "'\n"
-           << "window_geometry "
-           << position << " " << 0 << " " << 600 << " " << 600 << "\n"
-           << "keys jRmclA" << endl;
+      // 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;
 }
 
-class RelaxedNewtonSolver : public NewtonSolver
+class HessianCoefficient : public MatrixCoefficient
 {
 private:
-   // Quadrature points that are checked for negative Jacobians etc.
-   const IntegrationRule &ir;
-   ParFiniteElementSpace *pfes;
-   mutable ParGridFunction x_gf;
+   int type;
 
 public:
-   RelaxedNewtonSolver(const IntegrationRule &irule, ParFiniteElementSpace *pf)
-      : NewtonSolver(pf->GetComm()), ir(irule), pfes(pf) { }
+   HessianCoefficient(int dim, int type_)
+      : MatrixCoefficient(dim), type(type_) { }
 
-   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      Vector pos(3);
+      T.Transform(ip, pos);
+
+      if (type == 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;
+      }
+      else
+      {
+         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;
+
+         const double tan1 = std::tanh(sf*(r-r1)),
+                      tan2 = std::tanh(sf*(r-r2));
+
+         K(0, 0) = eps + 1.0 * (tan1 - tan2);
+         K(0, 1) = 0.0;
+         K(1, 0) = 0.0;
+         K(1, 1) = 1.0;
+      }
+   }
 };
 
-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++)
-   {
-      add(x, -scale, c, x_out);
-      x_gf.SetFromTrueVector();
-
-      energy_out = nlf->GetParGridFunctionEnergy(x_gf);
-      if (energy_out > 1.2*energy_in || std::isnan(energy_out) != 0)
-      {
-         if (print_level >= 0)
-         { cout << "Scale = " << scale << " Increasing energy." << endl; }
-         scale *= 0.5; continue;
-      }
-
-      int jac_ok = 1;
-      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);
-            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());
-
-      if (jac_ok_all == 0)
-      {
-         if (print_level >= 0)
-         { cout << "Scale = " << scale << " Neg det(J) found." << endl; }
-         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)
-         { 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.
@@ -331,7 +289,9 @@ 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");
+                  "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)");
    args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
    args.AddOption(&quad_type, "-qt", "--quad-type",
                   "Quadrature rule type:\n\t"
@@ -381,6 +341,7 @@ 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(); }
 
@@ -510,20 +471,47 @@ 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;
    }
-   TargetConstructor *target_c;
-   target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
+
+   if (target_c == NULL)
+   {
+      target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
+   }
    target_c->SetNodes(x0);
-   TMOP_Integrator *he_nlf_integ;
-   he_nlf_integ = new TMOP_Integrator(metric, target_c);
+   TMOP_Integrator *he_nlf_integ= new TMOP_Integrator(metric, target_c);
 
    // 13. Setup the quadrature rule for the non-linear form integrator.
    const IntegrationRule *ir = NULL;
@@ -597,7 +585,7 @@ int main (int argc, char *argv[])
    if (visualization)
    {
       char title[] = "Initial metric values";
-      vis_metric(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
+      vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
    }
 
    // 17. Fix all boundary nodes, or fix only a given component depending on the
@@ -705,9 +693,14 @@ int main (int argc, char *argv[])
    if (tauval > 0.0)
    {
       tauval = 0.0;
-      newton = new RelaxedNewtonSolver(*ir, pfespace);
+      TMOPNewtonSolver *tns = new TMOPNewtonSolver(pfespace->GetComm(), *ir);
+      if (target_id == 5)
+      {
+         tns->SetDiscreteAdaptTC(dynamic_cast<DiscreteAdaptTC *>(target_c));
+      }
+      newton = tns;
       if (myid == 0)
-      { cout << "RelaxedNewtonSolver is used (as all det(J) > 0)." << endl; }
+      { cout << "TMOPNewtonSolver is used (as all det(J) > 0)." << endl; }
    }
    else
    {
@@ -715,15 +708,15 @@ int main (int argc, char *argv[])
            (dim == 3 && metric_id != 352) )
       {
          if (myid == 0)
-         { cout << "The mesh is inverted. Use an untangling metric." << endl; }
+         { cout << "The mesh is inverted. Use an untangling metric.\n"; }
          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 DescentNewtonSolver(*ir, pfespace);
+      newton = new TMOPDescentNewtonSolver(pfespace->GetComm(), *ir);
       if (myid == 0)
-      { cout << "DescentNewtonSolver is used (as some det(J) < 0)." << endl; }
+      { cout << "TMOPDescentNewtonSolver is used (as some det(J) < 0).\n"; }
    }
    newton->SetPreconditioner(*S);
    newton->SetMaxIter(newton_iter);
@@ -775,7 +768,7 @@ int main (int argc, char *argv[])
    if (visualization)
    {
       char title[] = "Final metric values";
-      vis_metric(mesh_poly_deg, *metric, *target_c, *pmesh, title, 600);
+      vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 600);
    }
 
    // 23. Visualize the mesh displacement.
@@ -805,6 +798,7 @@ int main (int argc, char *argv[])
    delete metric2;
    delete coeff1;
    delete target_c;
+   delete adapt_coeff;
    delete metric;
    delete pfespace;
    delete fec;
diff --git a/miniapps/meshing/square01.mesh b/miniapps/meshing/square01.mesh
new file mode 100644
index 0000000000..bc75f8555c
--- /dev/null
+++ b/miniapps/meshing/square01.mesh
@@ -0,0 +1,114 @@
+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
diff --git a/tests/unit/fem/test_pa_coeff.cpp b/tests/unit/fem/test_pa_coeff.cpp
new file mode 100644
index 0000000000..afcf19b392
--- /dev/null
+++ b/tests/unit/fem/test_pa_coeff.cpp
@@ -0,0 +1,147 @@
+// 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