Some mods for compiling

Const parameter for mesh copy assignment

.gitlab/lassen.yml

@@ -30,5 +30,5 @@
 
 opt_mpi_cuda_xl_16_1_1_8:
   variables:
-    SPEC: "%xl@16.1.1.8 +mpi +cuda cuda_arch=sm_70"
+    SPEC: "%xl@16.1.1.8 +mpi +cuda cuda_arch=70"
   extends: .build_and_test_on_lassen

CHANGELOG

@@ -10,6 +10,7 @@
 
 Version 4.3.1 (development)
 ===========================
+- Added support for hr-adaptivity using TMOP-based error estimator.
 
 
 Version 4.3, released on July 29, 2021

fem/CMakeLists.txt

@@ -105,6 +105,7 @@ set(SRCS
   tmop/tmop_pa_w3.cpp
   tmop/tmop_pa_w3_c0.cpp
   tmop_tools.cpp
+  tmop_amr.cpp
   gslib.cpp
   transfer.cpp
   lor.cpp
@@ -164,6 +165,7 @@ set(HDRS
   tmop.hpp
   tmop/tmop_pa.hpp
   tmop_tools.hpp
+  tmop_amr.hpp
   gslib.hpp
   transfer.hpp
   lor.hpp

fem/fem.hpp

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

fem/gridfunc.hpp

@@ -95,6 +95,12 @@ public:
       : Vector(data, f->GetVSize())
    { fes = f; fec = NULL; fes_sequence = f->GetSequence(); UseDevice(true); }
 
+   /** @brief Construct a GridFunction using previously allocated Vector @a base
+       starting at the given offset, @a base_offset. */
+   GridFunction(FiniteElementSpace *f, Vector &base, int base_offset = 0)
+      : Vector(base, base_offset, f->GetVSize())
+   { fes = f; fec = NULL; fes_sequence = f->GetSequence(); UseDevice(true); }
+
    /// Construct a GridFunction on the given Mesh, using the data from @a input.
    /** The content of @a input should be in the format created by the method
        Save(). The reconstructed FiniteElementSpace and FiniteElementCollection

fem/gslib.cpp

@@ -610,7 +610,12 @@ void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
       {
          if (gsl_code[i] == 1) { indl2.Append(i); }
       }
-      if (indl2.Size() == 0) { return; } // no points on element borders
+      int borderPts = indl2.Size();
+#ifdef MFEM_USE_MPI
+      MPI_Allreduce(MPI_IN_PLACE, &borderPts, 1, MPI_INT, MPI_SUM, gsl_comm->c);
+#endif
+      if (borderPts == 0) { return; } // no points on element borders
+
 
       Vector field_out_l2(field_out.Size());
       VectorGridFunctionCoefficient field_in_dg(&field_in);

fem/pbilinearform.cpp

@@ -317,10 +317,6 @@ ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const double a)
 const
 {
-   MFEM_VERIFY(interior_face_integs.Size() == 0,
-               "the case of interior face integrators is not"
-               " implemented");
-
    if (X.ParFESpace() != pfes)
    {
       X.SetSpace(pfes);
@@ -334,6 +330,9 @@ const
    }
    else
    {
+      MFEM_VERIFY(interior_face_integs.Size() == 0,
+                  "the case of interior face integrators is not"
+                  " implemented");
       mat->Mult(X, Y);
    }
    pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Y, 1.0, y);

fem/pgridfunc.hpp

@@ -65,6 +65,11 @@ public:
    ParGridFunction(ParFiniteElementSpace *pf, double *data) :
       GridFunction(pf, data), pfes(pf) { }
 
+   /** @brief Construct a ParGridFunction using previously allocated Vector
+       @a base starting at the given offset, @a base_offset. */
+   ParGridFunction(ParFiniteElementSpace *pf, Vector &base, int base_offset = 0)
+      : GridFunction(pf, base, base_offset), pfes(pf) { }
+
    /// Construct a ParGridFunction using a GridFunction as external data.
    /** The parallel space @a *pf and the space used by @a *gf should match. The
        data from @a *gf is used as the local data of the ParGridFunction on each

fem/prestriction.cpp

@@ -33,6 +33,10 @@ ParL2FaceRestriction::ParL2FaceRestriction(const ParFiniteElementSpace &fes,
    // If fespace == L2
    const ParFiniteElementSpace &pfes =
       static_cast<const ParFiniteElementSpace&>(this->fes);
+
+   // Ensure the face neighbor data is constructed
+   pfes.GetParMesh()->ExchangeFaceNbrData();
+
    const FiniteElement *fe = pfes.GetFE(0);
    const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);
    MFEM_VERIFY(tfe != NULL &&

fem/restriction.cpp

@@ -15,6 +15,12 @@
 #include "../general/forall.hpp"
 #include <climits>
 
+#ifdef MFEM_USE_MPI
+
+#include "pfespace.hpp"
+
+#endif
+
 namespace mfem
 {
 
@@ -675,6 +681,19 @@ H1FaceRestriction::H1FaceRestriction(const FiniteElementSpace &fes,
      gather_indices(nf*dof)
 {
    if (nf==0) { return; }
+
+#ifdef MFEM_USE_MPI
+
+   // If the underlying finite element space is parallel, ensure the face
+   // neighbor information is generated.
+   if (const ParFiniteElementSpace *pfes
+       = dynamic_cast<const ParFiniteElementSpace*>(&fes))
+   {
+      pfes->GetParMesh()->ExchangeFaceNbrData();
+   }
+
+#endif
+
    // If fespace == H1
    const FiniteElement *fe = fes.GetFE(0);
    const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);

fem/tmop.cpp

@@ -1314,33 +1314,61 @@ static inline void device_copy(double *d_dest, const double *d_src, int size)
 } // namespace internal
 
 #ifdef MFEM_USE_MPI
-void DiscreteAdaptTC::FinalizeParDiscreteTargetSpec(const ParGridFunction
-                                                    &tspec_)
+void DiscreteAdaptTC::FinalizeParDiscreteTargetSpec(const ParGridFunction &t)
 {
    MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
    MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
 
-   ParFiniteElementSpace *ptspec_fes = tspec_.ParFESpace();
+   ParFiniteElementSpace *ptspec_fes = t.ParFESpace();
 
    adapt_eval->SetParMetaInfo(*ptspec_fes->GetParMesh(),
                               *ptspec_fes->FEColl(), ncomp);
-   adapt_eval->SetInitialField(*tspec_fes->GetMesh()->GetNodes(), tspec);
+   adapt_eval->SetInitialField(*ptspec_fes->GetMesh()->GetNodes(), tspec);
 
    tspec_sav = tspec;
 
    delete tspec_fesv;
-   tspec_fesv = new FiniteElementSpace(tspec_fes->GetMesh(),
-                                       tspec_fes->FEColl(), ncomp);
+   tspec_fesv = new FiniteElementSpace(ptspec_fes->GetMesh(),
+                                       ptspec_fes->FEColl(), ncomp);
+
+   delete ptspec_fesv;
+   ptspec_fesv = new ParFiniteElementSpace(ptspec_fes->GetParMesh(),
+                                           ptspec_fes->FEColl(), ncomp);
+
+   delete tspec_pgf;
+   tspec_pgf = new ParGridFunction(ptspec_fesv, tspec);
+   tspec_gf = tspec_pgf;
+}
+
+void DiscreteAdaptTC::ParUpdateAfterMeshTopologyChange()
+{
+   ptspec_fesv->Update();
+   if (tspec_fesv)
+   {
+      delete tspec_fesv;
+      tspec_fesv = new FiniteElementSpace(ptspec_fesv->GetMesh(),
+                                          ptspec_fesv->FEColl(), ncomp);
+   }
+   tspec_pgf->Update();
+   tspec_gf = tspec_pgf;
+   tspec.SetDataAndSize(tspec_pgf->GetData(), tspec_pgf->Size());
+   tspec_sav = tspec;
+
+   adapt_eval->SetParMetaInfo(*ptspec_fesv->GetParMesh(),
+                              *ptspec_fesv->FEColl(), ncomp);
+   adapt_eval->SetInitialField(*ptspec_fesv->GetMesh()->GetNodes(), tspec);
 }
 
 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const ParGridFunction &tspec_)
 {
-   const int vdim     = tspec_.FESpace()->GetVDim(),
-             dof_cnt  = tspec_.Size()/vdim;
+   const int vdim = tspec_.FESpace()->GetVDim(),
+             ndof = tspec_.FESpace()->GetNDofs();
+   MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTspecAtIndex.");
+
    const auto tspec__d = tspec_.Read();
    auto tspec_d = tspec.ReadWrite();
-   const int offset = idx*dof_cnt;
-   internal::device_copy(tspec_d + offset, tspec__d, dof_cnt*vdim);
+   const int offset = idx*ndof;
+   internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
    FinalizeParDiscreteTargetSpec(tspec_);
 }
 
@@ -1360,78 +1388,71 @@ void DiscreteAdaptTC::SetParDiscreteTargetSkew(const ParGridFunction &tspec_)
    FinalizeParDiscreteTargetSpec(tspec_);
 }
 
-void DiscreteAdaptTC::SetParDiscreteTargetAspectRatio(const ParGridFunction
-                                                      &tspec_)
+void DiscreteAdaptTC::SetParDiscreteTargetAspectRatio(const ParGridFunction &ar)
 {
-   if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, tspec_); return; }
+   if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
    aspectratioidx = ncomp;
-   SetDiscreteTargetBase(tspec_);
-   FinalizeParDiscreteTargetSpec(tspec_);
+   SetDiscreteTargetBase(ar);
+   FinalizeParDiscreteTargetSpec(ar);
 }
 
-void DiscreteAdaptTC::SetParDiscreteTargetOrientation(const ParGridFunction
-                                                      &tspec_)
+void DiscreteAdaptTC::SetParDiscreteTargetOrientation(const ParGridFunction &o)
 {
-   if (orientationidx > -1) { SetTspecAtIndex(orientationidx, tspec_); return; }
+   if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
    orientationidx = ncomp;
-   SetDiscreteTargetBase(tspec_);
-   FinalizeParDiscreteTargetSpec(tspec_);
+   SetDiscreteTargetBase(o);
+   FinalizeParDiscreteTargetSpec(o);
 }
 
 void DiscreteAdaptTC::SetParDiscreteTargetSpec(const ParGridFunction &tspec_)
 {
    SetParDiscreteTargetSize(tspec_);
-   FinalizeParDiscreteTargetSpec(tspec_);
 }
 #endif // MFEM_USE_MPI
 
 void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
 {
-   const int vdim     = tspec_.FESpace()->GetVDim(),
-             dof_cnt  = tspec_.Size()/vdim;
+   const int vdim = tspec_.FESpace()->GetVDim(),
+             ndof = tspec_.FESpace()->GetNDofs();
 
    ncomp += vdim;
 
-   delete tspec_fes;
-   tspec_fes = new FiniteElementSpace(tspec_.FESpace()->GetMesh(),
-                                      tspec_.FESpace()->FEColl(), 1);
-
    // need to append data to tspec
    // make a copy of tspec->tspec_temp, increase its size, and
    // copy data from tspec_temp -> tspec, then add new entries
    Vector tspec_temp = tspec;
    tspec.UseDevice(true);
    tspec_sav.UseDevice(true);
-   tspec.SetSize(ncomp*dof_cnt);
+   tspec.SetSize(ncomp*ndof);
 
    const auto tspec_temp_d = tspec_temp.Read();
    auto tspec_d = tspec.ReadWrite();
    internal::device_copy(tspec_d, tspec_temp_d, tspec_temp.Size());
 
    const auto tspec__d = tspec_.Read();
-   const int offset = (ncomp-vdim)*dof_cnt;
-   internal::device_copy(tspec_d + offset, tspec__d, dof_cnt*vdim);
+   const int offset = (ncomp-vdim)*ndof;
+   internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
 }
 
 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const GridFunction &tspec_)
 {
-   const int vdim     = tspec_.FESpace()->GetVDim(),
-             dof_cnt  = tspec_.Size()/vdim;
+   const int vdim = tspec_.FESpace()->GetVDim(),
+             ndof = tspec_.FESpace()->GetNDofs();
+   MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTargetSpec.");
 
    const auto tspec__d = tspec_.Read();
    auto tspec_d = tspec.ReadWrite();
-   const int offset = idx*dof_cnt;
-   internal::device_copy(tspec_d + offset, tspec__d, dof_cnt*vdim);
-   FinalizeSerialDiscreteTargetSpec();
+   const int offset = idx*ndof;
+   internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
+   FinalizeSerialDiscreteTargetSpec(tspec_);
 }
 
 void DiscreteAdaptTC::SetSerialDiscreteTargetSize(const GridFunction &tspec_)
 {
-
    if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
    sizeidx = ncomp;
    SetDiscreteTargetBase(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
+   FinalizeSerialDiscreteTargetSpec(tspec_);
 }
 
 void DiscreteAdaptTC::SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
@@ -1439,32 +1460,31 @@ void DiscreteAdaptTC::SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
    if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
    skewidx = ncomp;
    SetDiscreteTargetBase(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
+   FinalizeSerialDiscreteTargetSpec(tspec_);
 }
 
-void DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio(
-   const GridFunction &tspec_)
+void DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio(const GridFunction &ar)
 {
-   if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, tspec_); return; }
+   if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
    aspectratioidx = ncomp;
-   SetDiscreteTargetBase(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
+   SetDiscreteTargetBase(ar);
+   FinalizeSerialDiscreteTargetSpec(ar);
 }
 
-void DiscreteAdaptTC::SetSerialDiscreteTargetOrientation(
-   const GridFunction &tspec_)
+void DiscreteAdaptTC::SetSerialDiscreteTargetOrientation(const GridFunction &o)
 {
-   if (orientationidx > -1) { SetTspecAtIndex(orientationidx, tspec_); return; }
+   if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
    orientationidx = ncomp;
-   SetDiscreteTargetBase(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
+   SetDiscreteTargetBase(o);
+   FinalizeSerialDiscreteTargetSpec(o);
 }
 
-void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec()
+void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec(const GridFunction &t)
 {
    MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
    MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
 
+   const FiniteElementSpace *tspec_fes = t.FESpace();
    adapt_eval->SetSerialMetaInfo(*tspec_fes->GetMesh(),
                                  *tspec_fes->FEColl(), ncomp);
    adapt_eval->SetInitialField(*tspec_fes->GetMesh()->GetNodes(), tspec);
@@ -1474,12 +1494,40 @@ void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec()
    delete tspec_fesv;
    tspec_fesv = new FiniteElementSpace(tspec_fes->GetMesh(),
                                        tspec_fes->FEColl(), ncomp);
+
+   delete tspec_gf;
+   tspec_gf = new GridFunction(tspec_fesv, tspec);
+}
+
+void DiscreteAdaptTC::GetDiscreteTargetSpec(GridFunction &tspec_, int idx)
+{
+   if (idx < 0) { return; }
+   const int ndof = tspec_.FESpace()->GetNDofs(),
+             vdim = tspec_.FESpace()->GetVDim();
+   MFEM_VERIFY(ndof == tspec.Size()/ncomp,
+               "Inconsistency in GetSerialDiscreteTargetSpec.");
+
+   for (int i = 0; i < ndof*vdim; i++)
+   {
+      tspec_(i) = tspec(i + idx*ndof);
+   }
+}
+
+void DiscreteAdaptTC::UpdateAfterMeshTopologyChange()
+{
+   tspec_fesv->Update();
+   tspec_gf->Update();
+   tspec.SetDataAndSize(tspec_gf->GetData(), tspec_gf->Size());
+   tspec_sav = tspec;
+
+   adapt_eval->SetSerialMetaInfo(*tspec_fesv->GetMesh(),
+                                 *tspec_fesv->FEColl(), ncomp);
+   adapt_eval->SetInitialField(*tspec_fesv->GetMesh()->GetNodes(), tspec);
 }
 
 void DiscreteAdaptTC::SetSerialDiscreteTargetSpec(const GridFunction &tspec_)
 {
    SetSerialDiscreteTargetSize(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
 }
 
 
@@ -1509,7 +1557,7 @@ void DiscreteAdaptTC::UpdateTargetSpecificationAtNode(const FiniteElement &el,
    MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
 
    Array<int> dofs;
-   tspec_fes->GetElementDofs(T.ElementNo, dofs);
+   tspec_fesv->GetElementDofs(T.ElementNo, dofs);
    const int cnt = tspec.Size()/ncomp; // dofs per scalar-field
 
    for (int i = 0; i < ncomp; i++)
@@ -1524,7 +1572,7 @@ void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
    MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
 
    Array<int> dofs;
-   tspec_fes->GetElementDofs(T.ElementNo, dofs);
+   tspec_fesv->GetElementDofs(T.ElementNo, dofs);
    const int cnt = tspec.Size()/ncomp;
    for (int i = 0; i < ncomp; i++)
    {
@@ -1532,6 +1580,40 @@ void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
    }
 }
 
+void DiscreteAdaptTC::SetTspecFromIntRule(int e_id,
+                                          const IntegrationRule &intrule)
+{
+   switch (target_type)
+   {
+      case IDEAL_SHAPE_GIVEN_SIZE:
+      case GIVEN_SHAPE_AND_SIZE:
+      {
+         const int ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
+                   ntspec_dofs = ndofs*ncomp;
+
+         Vector tspec_vals(ntspec_dofs);
+
+         Array<int> dofs;
+         tspec_fesv->GetElementVDofs(e_id, dofs);
+         tspec.GetSubVector(dofs, tspec_vals);
+         DenseMatrix tr;
+         tspec_gf->GetVectorValues(e_id, intrule, tspec_refine, tr);
+         tspec_refine.Transpose();
+         break;
+      }
+      default:
+         MFEM_ABORT("Incompatible target type for discrete adaptation!");
+   }
+}
+
+void DiscreteAdaptTC::SetTspecDataForDerefinement(FiniteElementSpace *fes)
+{
+   coarse_tspec_fesv = fes;
+   const Operator *c_op = fes->GetUpdateOperator();
+   tspec_derefine.SetSize(c_op->Height());
+   c_op->Mult(tspec, tspec_derefine);
+}
+
 void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
                                             const IntegrationRule &ir,
                                             const Vector &elfun,
@@ -1542,6 +1624,8 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
              nqp = ir.GetNPoints();
    Jtrcomp.SetSize(dim, dim, 4*nqp);
 
+   FiniteElementSpace *src_fes = tspec_fesv;
+
    switch (target_type)
    {
       case IDEAL_SHAPE_GIVEN_SIZE:
@@ -1550,7 +1634,7 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
          const DenseMatrix &Wideal =
             Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
          const int dim = Wideal.Height(),
-                   ndofs = tspec_fes->GetFE(e_id)->GetDof(),
+                   ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
                    ntspec_dofs = ndofs*ncomp;
 
          Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
@@ -1561,11 +1645,29 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
          tspec_fesv->GetElementVDofs(e_id, dofs);
          tspec.UseDevice(true);
          tspec.GetSubVector(dofs, tspec_vals);
+         if (tspec_refine.NumCols() > 0) // Refinement
+         {
+            MFEM_VERIFY(amr_el >= 0, " Target being constructed for an AMR element.");
+            for (int i = 0; i < ncomp; i++)
+            {
+               for (int j = 0; j < ndofs; j++)
+               {
+                  tspec_vals(j + i*ndofs) = tspec_refine(j + amr_el*ndofs, i);
+               }
+            }
+         }
+         else if (tspec_derefine.Size() > 0) // Derefinement
+         {
+            dofs.SetSize(0);
+            coarse_tspec_fesv->GetElementVDofs(e_id, dofs);
+            tspec_derefine.GetSubVector(dofs, tspec_vals);
+            src_fes = coarse_tspec_fesv;
+         }
 
          for (int q = 0; q < nqp; q++)
          {
             const IntegrationPoint &ip = ir.IntPoint(q);
-            tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+            src_fes->GetFE(e_id)->CalcShape(ip, shape);
             Jtr(q) = Wideal; // Initialize to identity
             for (int d = 0; d < 4; d++)
             {
@@ -1576,9 +1678,16 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
             if (sizeidx != -1) // Set size
             {
                par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
-               const double min_size = par_vals.Min();
-               MFEM_VERIFY(min_size > 0.0,
-                           "Non-positive size propagated in the target definition.");
+               double min_size = par_vals.Min();//0.001; //
+               if (lim_min_size > 0.)
+               {
+                  min_size = lim_min_size;
+               }
+               else
+               {
+                  MFEM_VERIFY(min_size > 0.0,
+                              "Non-positive size propagated in the target definition.");
+               }
                const double size = std::max(shape * par_vals, min_size);
                Jtr(q).Set(std::pow(size, 1.0/dim), Jtr(q));
                DenseMatrix Jtrcomp_q(Jtrcomp.GetData(0 + 4*q), dim, dim);
@@ -1593,6 +1702,9 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
                {
                   par_vals.SetDataAndSize(tspec_vals.GetData()+
                                           aspectratioidx*ndofs, ndofs);
+                  const double min_size = par_vals.Min();
+                  MFEM_VERIFY(min_size > 0.0,
+                              "Non-positive aspect-ratio propagated in the target definition.");
 
                   const double aspectratio = shape * par_vals;
                   D_rho = 0.;
@@ -1777,7 +1889,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
 
                grad_phys.Mult(par_vals, grad_ptr_c1);
                Vector grad_q(dim);
-               tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+               tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                grad_e_c1.MultTranspose(shape, grad_q);
 
                const double min_size = par_vals.Min();
@@ -1810,7 +1922,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
 
                   grad_phys.Mult(par_vals, grad_ptr_c1);
                   Vector grad_q(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q);
 
                   const double aspectratio = shape * par_vals;
@@ -1841,7 +1953,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
                   grad_phys.Mult(par_vals_c2, grad_ptr_c2);
                   grad_phys.Mult(par_vals_c3, grad_ptr_c3);
                   Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q1);
                   grad_e_c2.MultTranspose(shape, grad_q2);
                   grad_e_c3.MultTranspose(shape, grad_q3);
@@ -1880,7 +1992,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
 
                   grad_phys.Mult(par_vals, grad_ptr_c1);
                   Vector grad_q(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q);
 
                   const double skew = shape * par_vals;
@@ -1913,7 +2025,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
                   grad_phys.Mult(par_vals_c2, grad_ptr_c2);
                   grad_phys.Mult(par_vals_c3, grad_ptr_c3);
                   Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q1);
                   grad_e_c2.MultTranspose(shape, grad_q2);
                   grad_e_c3.MultTranspose(shape, grad_q3);
@@ -1960,7 +2072,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
 
                   grad_phys.Mult(par_vals, grad_ptr_c1);
                   Vector grad_q(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q);
 
                   const double theta = shape * par_vals;
@@ -1991,7 +2103,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
                   grad_phys.Mult(par_vals_c2, grad_ptr_c2);
                   grad_phys.Mult(par_vals_c3, grad_ptr_c3);
                   Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q1);
                   grad_e_c2.MultTranspose(shape, grad_q2);
                   grad_e_c3.MultTranspose(shape, grad_q3);
@@ -2071,7 +2183,7 @@ void DiscreteAdaptTC::UpdateGradientTargetSpecification(const Vector &x,
 {
    if (use_flag && good_tspec_grad) { return; }
 
-   const int dim = tspec_fes->GetFE(0)->GetDim(),
+   const int dim = tspec_fesv->GetFE(0)->GetDim(),
              cnt = x.Size()/dim;
 
    tspec_pert1h.SetSize(x.Size()*ncomp);
@@ -2097,7 +2209,7 @@ void DiscreteAdaptTC::UpdateHessianTargetSpecification(const Vector &x,
 
    if (use_flag && good_tspec_hess) { return; }
 
-   const int dim    = tspec_fes->GetFE(0)->GetDim(),
+   const int dim    = tspec_fesv->GetFE(0)->GetDim(),
              cnt    = x.Size()/dim,
              totmix = 1+2*(dim-2);
 
@@ -2145,6 +2257,16 @@ void DiscreteAdaptTC::UpdateHessianTargetSpecification(const Vector &x,
    good_tspec_hess = use_flag;
 }
 
+DiscreteAdaptTC::~DiscreteAdaptTC()
+{
+   delete tspec_gf;
+   delete adapt_eval;
+   delete tspec_fesv;
+#ifdef MFEM_USE_MPI
+   delete ptspec_fesv;
+#endif
+}
+
 void AdaptivityEvaluator::SetSerialMetaInfo(const Mesh &m,
                                             const FiniteElementCollection &fec,
                                             int num_comp)
@@ -2258,6 +2380,7 @@ void TMOP_Integrator::EnableAdaptiveLimiting(const ParGridFunction &z0,
                                              AdaptivityEvaluator &ae)
 {
    zeta_0 = &z0;
+   pzeta_0 = &z0;
    delete zeta;
    zeta   = new GridFunction(z0);
    coeff_zeta = &coeff;
@@ -2270,6 +2393,33 @@ void TMOP_Integrator::EnableAdaptiveLimiting(const ParGridFunction &z0,
 }
 #endif
 
+void TMOP_Integrator::UpdateAfterMeshTopologyChange()
+{
+   if (zeta)
+   {
+      zeta->Update();
+      adapt_eval->SetSerialMetaInfo(*zeta->FESpace()->GetMesh(),
+                                    *zeta->FESpace()->FEColl(), 1);
+      adapt_eval->SetInitialField
+      (*zeta->FESpace()->GetMesh()->GetNodes(), *zeta);
+   }
+}
+
+#ifdef MFEM_USE_MPI
+void TMOP_Integrator::ParUpdateAfterMeshTopologyChange()
+{
+   if (zeta)
+   {
+      zeta->Update();
+      adapt_eval->SetParMetaInfo(*pzeta_0->ParFESpace()->GetParMesh(),
+                                 *pzeta_0->ParFESpace()->FEColl(), 1);
+      adapt_eval->SetInitialField
+      (*zeta->FESpace()->GetMesh()->GetNodes(), *zeta);
+   }
+}
+#endif
+
+
 double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
                                          ElementTransformation &T,
                                          const Vector &elfun)
@@ -2378,6 +2528,145 @@ double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
 
    return energy;
 }
+
+double TMOP_Integrator::GetRefinementElementEnergy(const FiniteElement &el,
+                                                   ElementTransformation &T,
+                                                   const Vector &elfun,
+                                                   const IntegrationRule &irule)
+{
+   int dof = el.GetDof(), dim = el.GetDim(),
+       NEsplit = elfun.Size() / (dof*dim), el_id = T.ElementNo;
+   double energy = 0.;
+
+   TargetConstructor *tc = const_cast<TargetConstructor *>(targetC);
+   DiscreteAdaptTC *dtc = dynamic_cast<DiscreteAdaptTC *>(tc);
+   // For DiscreteAdaptTC the GridFunctions used to set the targets must be
+   // mapped onto the fine elements.
+   if (dtc) { dtc->SetTspecFromIntRule(el_id, irule); }
+
+   for (int e = 0; e < NEsplit; e++)
+   {
+      DSh.SetSize(dof, dim);
+      Jrt.SetSize(dim);
+      Jpr.SetSize(dim);
+      Jpt.SetSize(dim);
+      Vector elfun_child(dof*dim);
+      for (int i = 0; i < dof; i++)
+      {
+         for (int d = 0; d < dim; d++)
+         {
+            // elfun is (xe1,xe2,...xen,ye1,ye2...yen) and has nodal coordinates
+            // for all the children element of the parent element being considered.
+            // So we must index and get (xek, yek) i.e. nodal coordinates for
+            // the fine element being considered.
+            elfun_child(i + d*dof) = elfun(i + e*dof + d*dof*NEsplit);
+         }
+      }
+      PMatI.UseExternalData(elfun_child.GetData(), dof, dim);
+
+      const IntegrationRule &ir = EnergyIntegrationRule(el);
+
+      double el_energy = 0;
+      DenseTensor Jtr(dim, dim, ir.GetNPoints());
+      if (dtc)
+      {
+         // This is used to index into the tspec vector inside DiscreteAdaptTC.
+         dtc->SetRefinementSubElement(e);
+      }
+      targetC->ComputeElementTargets(el_id, el, ir, elfun_child, Jtr);
+
+      // Define ref->physical transformation, wn a Coefficient is specified.
+      IsoparametricTransformation *Tpr = NULL;
+      if (coeff1 || coeff0)
+      {
+         Tpr = new IsoparametricTransformation;
+         Tpr->SetFE(&el);
+         Tpr->ElementNo = T.ElementNo;
+         Tpr->ElementType = ElementTransformation::ELEMENT;
+         Tpr->Attribute = T.Attribute;
+         Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
+      }
+
+      for (int i = 0; i < ir.GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(i);
+         const DenseMatrix &Jtr_i = Jtr(i);
+         h_metric->SetTargetJacobian(Jtr_i);
+         CalcInverse(Jtr_i, Jrt);
+         const double weight = ip.weight * Jtr_i.Det();
+
+         el.CalcDShape(ip, DSh);
+         MultAtB(PMatI, DSh, Jpr);
+         Mult(Jpr, Jrt, Jpt);
+
+         double val = metric_normal * h_metric->EvalW(Jpt);
+         if (coeff1) { val *= coeff1->Eval(*Tpr, ip); }
+
+         el_energy += weight * val;
+         delete Tpr;
+      }
+      energy += el_energy;
+   }
+   energy /= NEsplit;
+
+   if (dtc) { dtc->ResetRefinementTspecData(); }
+
+   return energy;
+}
+
+double TMOP_Integrator::GetDerefinementElementEnergy(const FiniteElement &el,
+                                                     ElementTransformation &T,
+                                                     const Vector &elfun)
+{
+   int dof = el.GetDof(), dim = el.GetDim();
+   double energy = 0.;
+
+   DSh.SetSize(dof, dim);
+   Jrt.SetSize(dim);
+   Jpr.SetSize(dim);
+   Jpt.SetSize(dim);
+   PMatI.UseExternalData(elfun.GetData(), dof, dim);
+
+   const IntegrationRule &ir = EnergyIntegrationRule(el);
+
+   energy = 0.0;
+   DenseTensor Jtr(dim, dim, ir.GetNPoints());
+   targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
+
+   // Define ref->physical transformation, wn a Coefficient is specified.
+   IsoparametricTransformation *Tpr = NULL;
+   if (coeff1)
+   {
+      Tpr = new IsoparametricTransformation;
+      Tpr->SetFE(&el);
+      Tpr->ElementNo = T.ElementNo;
+      Tpr->ElementType = ElementTransformation::ELEMENT;
+      Tpr->Attribute = T.Attribute;
+      Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
+   }
+
+   for (int i = 0; i < ir.GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir.IntPoint(i);
+      const DenseMatrix &Jtr_i = Jtr(i);
+      h_metric->SetTargetJacobian(Jtr_i);
+      CalcInverse(Jtr_i, Jrt);
+      const double weight = ip.weight * Jtr_i.Det();
+
+      el.CalcDShape(ip, DSh);
+      MultAtB(PMatI, DSh, Jpr);
+      Mult(Jpr, Jrt, Jpt);
+
+      double val = metric_normal * h_metric->EvalW(Jpt);
+      if (coeff1) { val *= coeff1->Eval(*Tpr, ip); }
+
+      energy += weight * val;
+   }
+
+   delete Tpr;
+   return energy;
+}
+
 void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
                                             ElementTransformation &T,
                                             const Vector &elfun, Vector &elvect)
@@ -3039,7 +3328,7 @@ void TMOP_Integrator::ComputeMinJac(const Vector &x,
    dx = detv_avg_min / dxscale;
 }
 
-void TMOP_Integrator::UpdateAfterMeshChange(const Vector &new_x)
+void TMOP_Integrator::UpdateAfterMeshPositionChange(const Vector &new_x)
 {
    if (discr_tc)
    {
@@ -3168,6 +3457,32 @@ void TMOPComboIntegrator::AssembleElementGrad(const FiniteElement &el,
    }
 }
 
+double TMOPComboIntegrator::GetRefinementElementEnergy(const FiniteElement &el,
+                                                       ElementTransformation &T,
+                                                       const Vector &elfun,
+                                                       const IntegrationRule &irule)
+{
+   double energy= 0.0;
+   for (int i = 0; i < tmopi.Size(); i++)
+   {
+      energy += tmopi[i]->GetRefinementElementEnergy(el, T, elfun, irule);
+   }
+   return energy;
+}
+
+double TMOPComboIntegrator::GetDerefinementElementEnergy(
+   const FiniteElement &el,
+   ElementTransformation &T,
+   const Vector &elfun)
+{
+   double energy= 0.0;
+   for (int i = 0; i < tmopi.Size(); i++)
+   {
+      energy += tmopi[i]->GetDerefinementElementEnergy(el, T, elfun);
+   }
+   return energy;
+}
+
 void TMOPComboIntegrator::EnableNormalization(const GridFunction &x)
 {
    const int cnt = tmopi.Size();

fem/tmop.hpp

@@ -1057,14 +1057,31 @@ protected:
    // eta1(x+h,y), eta2(x+h,y) ... etan(x+h,y), eta1(x,y+h), eta2(x,y+h) ...
    // same for tspec_pert2h and tspec_pertmix.
 
+   // DenseMatrix to hold target_spec values for the (children of the)
+   // element being refined to consider for h-refinement.
+   DenseMatrix tspec_refine;
+   // Vector to hold the target_spec values for the coarse version of the
+   // current mesh. Used for derefinement decision with hr-adaptivity.
+   Vector tspec_derefine;
+
    // Components of Target Jacobian at each quadrature point of an element. This
    // is required for computation of the derivative using chain rule.
    mutable DenseTensor Jtrcomp;
 
    // 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;
-   const FiniteElementSpace *tspec_fesv;
+   FiniteElementSpace *tspec_fesv;         //owned
+   FiniteElementSpace *coarse_tspec_fesv;  //not owned, derefinement FESpace
+   GridFunction *tspec_gf;                 //owned, uses tspec and tspec_fes
+   // discrete adaptivity
+#ifdef MFEM_USE_MPI
+   ParFiniteElementSpace *ptspec_fesv;     //owned, needed for derefinement to
+   // get update operator.
+   ParGridFunction *tspec_pgf;             // similar to tspec_gf
+#endif
+
+   int amr_el;
+   double lim_min_size;
 
    // These flags can be used by outside functions to avoid recomputing the
    // tspec and tspec_perth fields again on the same mesh.
@@ -1076,7 +1093,7 @@ protected:
 
    void SetDiscreteTargetBase(const GridFunction &tspec_);
    void SetTspecAtIndex(int idx, const GridFunction &tspec_);
-   void FinalizeSerialDiscreteTargetSpec();
+   void FinalizeSerialDiscreteTargetSpec(const GridFunction &tspec_);
 #ifdef MFEM_USE_MPI
    void SetTspecAtIndex(int idx, const ParGridFunction &tspec_);
    void FinalizeParDiscreteTargetSpec(const ParGridFunction &tspec_);
@@ -1088,16 +1105,16 @@ public:
         ncomp(0),
         sizeidx(-1), skewidx(-1), aspectratioidx(-1), orientationidx(-1),
         tspec(), tspec_sav(), tspec_pert1h(), tspec_pert2h(), tspec_pertmix(),
-        tspec_fes(NULL), tspec_fesv(NULL),
+        tspec_refine(), tspec_derefine(),
+        tspec_fesv(NULL), coarse_tspec_fesv(NULL), tspec_gf(NULL),
+#ifdef MFEM_USE_MPI
+        ptspec_fesv(NULL), tspec_pgf(NULL),
+#endif
+        amr_el(-1), lim_min_size(-0.1),
         good_tspec(false), good_tspec_grad(false), good_tspec_hess(false),
         adapt_eval(NULL) { }
 
-   virtual ~DiscreteAdaptTC()
-   {
-      delete adapt_eval;
-      delete tspec_fes;
-      delete tspec_fesv;
-   }
+   virtual ~DiscreteAdaptTC();
 
    /** @name Target specification methods.
        The following methods are used to specify geometric parameters of the
@@ -1128,6 +1145,20 @@ public:
    void ResetUpdateFlags()
    { good_tspec = good_tspec_grad = good_tspec_hess = false; }
 
+   /// Get one of the discrete fields from tspec.
+   void GetDiscreteTargetSpec(GridFunction &tspec_, int idx);
+   /// Get the FESpace associated with tspec.
+   FiniteElementSpace *GetTSpecFESpace() { return tspec_fesv; }
+   /// Get the entire tspec.
+   GridFunction *GetTSpecData() { return tspec_gf; }
+   /// Update all discrete fields based on tspec and update for AMR
+   void UpdateAfterMeshTopologyChange();
+
+#ifdef MFEM_USE_MPI
+   ParFiniteElementSpace *GetTSpecParFESpace() { return ptspec_fesv; }
+   void ParUpdateAfterMeshTopologyChange();
+#endif
+
    /** Used to update the target specification after the mesh has changed. The
        new mesh positions are given by new_x. If @a use_flags is true, repeated
        calls won't do anything until ResetUpdateFlags() is called. */
@@ -1184,6 +1215,36 @@ public:
                                               const Vector &elfun,
                                               IsoparametricTransformation &Tpr,
                                               DenseTensor &dJtr) const;
+
+   // Generates tspec_vals for target construction using intrule
+   // Used for the refinement component in hr-adaptivity.
+   void SetTspecFromIntRule(int e_id, const IntegrationRule &intrule);
+
+   // Targets based on discrete functions can result in invalid (negative)
+   // size at the quadrature points. This method can be used to set a
+   // minimum target size.
+   void SetMinSizeForTargets(double min_size_) { lim_min_size = min_size_; }
+
+   /// Computes target specification data with respect to the coarse FE space.
+   void SetTspecDataForDerefinement(FiniteElementSpace *fes);
+
+   // Reset refinement data associated with h-adaptivity component.
+   void ResetRefinementTspecData()
+   {
+      tspec_refine.Clear();
+      amr_el = -1;
+   }
+
+   // Reset derefinement data associated with h-adaptivity component.
+   void ResetDerefinementTspecData()
+   {
+      tspec_derefine.Destroy();
+      coarse_tspec_fesv = NULL;
+   }
+
+   // Used to specify the fine element for determining energy of children of a
+   // parent element.
+   void SetRefinementSubElement(int amr_el_) { amr_el = amr_el_; }
 };
 
 class TMOPNewtonSolver;
@@ -1201,6 +1262,7 @@ protected:
    friend class TMOPNewtonSolver;
    friend class TMOPComboIntegrator;
 
+   TMOP_QualityMetric *h_metric;
    TMOP_QualityMetric *metric;        // not owned
    const TargetConstructor *targetC;  // not owned
 
@@ -1227,6 +1289,9 @@ protected:
 
    // Adaptive limiting.
    const GridFunction *zeta_0;       // Not owned.
+#ifdef MFEM_USE_MPI
+   const ParGridFunction *pzeta_0;
+#endif
    GridFunction *zeta;               // Owned. Updated by adapt_eval.
    Coefficient *coeff_zeta;          // Not owned.
    AdaptivityEvaluator *adapt_eval;  // Not owned.
@@ -1337,7 +1402,7 @@ protected:
 #endif
    void ComputeMinJac(const Vector &x, const FiniteElementSpace &fes);
 
-   void UpdateAfterMeshChange(const Vector &new_x);
+   void UpdateAfterMeshPositionChange(const Vector &new_x);
 
    void DisableLimiting()
    {
@@ -1395,11 +1460,13 @@ protected:
    void ComputeAllElementTargets(const Vector &xe = Vector()) const;
 
 public:
-   /** @param[in] m  TMOP_QualityMetric that will be integrated (not owned).
-       @param[in] tc Target-matrix construction algorithm to use (not owned). */
-   TMOP_Integrator(TMOP_QualityMetric *m, TargetConstructor *tc)
-      : metric(m), targetC(tc), IntegRules(NULL), integ_order(-1),
-        coeff1(NULL), metric_normal(1.0),
+   /** @param[in] m    TMOP_QualityMetric for r-adaptivity (not owned).
+       @param[in] tc   Target-matrix construction algorithm to use (not owned).
+       @param[in] hm   TMOP_QualityMetric for h-adaptivity (not owned). */
+   TMOP_Integrator(TMOP_QualityMetric *m, TargetConstructor *tc,
+                   TMOP_QualityMetric *hm)
+      : h_metric(hm), metric(m), targetC(tc), IntegRules(NULL),
+        integ_order(-1), coeff1(NULL), metric_normal(1.0),
         nodes0(NULL), coeff0(NULL),
         lim_dist(NULL), lim_func(NULL), lim_normal(1.0),
         zeta_0(NULL), zeta(NULL), coeff_zeta(NULL), adapt_eval(NULL),
@@ -1407,6 +1474,9 @@ public:
         fdflag(false), dxscale(1.0e3), fd_call_flag(false), exact_action(false)
    { PA.enabled = false; }
 
+   TMOP_Integrator(TMOP_QualityMetric *m, TargetConstructor *tc)
+      : TMOP_Integrator(m, tc, m) { }
+
    ~TMOP_Integrator();
 
    /// Release the device memory of large PA allocations. This will copy device
@@ -1478,6 +1548,22 @@ public:
                                    ElementTransformation &T,
                                    const Vector &elfun);
 
+   /** @brief Computes the mean of the energies of the given element's children.
+
+       In addition to the inputs for GetElementEnergy, this function requires an
+       IntegrationRule to be specified that will give the decomposition of the
+       given element based on the refinement type being considered. */
+   virtual double GetRefinementElementEnergy(const FiniteElement &el,
+                                             ElementTransformation &T,
+                                             const Vector &elfun,
+                                             const IntegrationRule &irule);
+
+   /// This function is similar to GetElementEnergy, but ignores components
+   /// such as limiting etc. to compute the element energy.
+   virtual double GetDerefinementElementEnergy(const FiniteElement &el,
+                                               ElementTransformation &T,
+                                               const Vector &elfun);
+
    virtual void AssembleElementVector(const FiniteElement &el,
                                       ElementTransformation &T,
                                       const Vector &elfun, Vector &elvect);
@@ -1486,6 +1572,13 @@ public:
                                     ElementTransformation &T,
                                     const Vector &elfun, DenseMatrix &elmat);
 
+   TMOP_QualityMetric &GetAMRQualityMetric() { return *h_metric; }
+
+   void UpdateAfterMeshTopologyChange();
+#ifdef MFEM_USE_MPI
+   void ParUpdateAfterMeshTopologyChange();
+#endif
+
    // PA extension
    using NonlinearFormIntegrator::AssemblePA;
    virtual void AssemblePA(const FiniteElementSpace&);
@@ -1564,6 +1657,15 @@ public:
                                     ElementTransformation &T,
                                     const Vector &elfun, DenseMatrix &elmat);
 
+   virtual double GetRefinementElementEnergy(const FiniteElement &el,
+                                             ElementTransformation &T,
+                                             const Vector &elfun,
+                                             const IntegrationRule &irule);
+
+   virtual double GetDerefinementElementEnergy(const FiniteElement &el,
+                                               ElementTransformation &T,
+                                               const Vector &elfun);
+
    /// Normalization factor that considers all integrators in the combination.
    void EnableNormalization(const GridFunction &x);
 #ifdef MFEM_USE_MPI

fem/tmop_amr.cpp

@@ -0,0 +1,896 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "tmop_amr.hpp"
+
+namespace mfem
+{
+
+using namespace mfem;
+
+void TMOPRefinerEstimator::ComputeEstimates()
+{
+   bool iso = false;
+   bool aniso = false;
+   if (amrmetric == 1 || amrmetric == 2 || amrmetric == 58)
+   {
+      aniso = true;
+   }
+   if (amrmetric == 55 || amrmetric == 56 || amrmetric == 77 ||
+       amrmetric == 315 || amrmetric == 316 || amrmetric == 321)
+   {
+      iso = true;
+   }
+   if (amrmetric == 7 || amrmetric == 9)
+   {
+      iso = true; aniso = true;
+   }
+
+   MFEM_VERIFY(iso || aniso, "Metric type not supported in hr-adaptivity.");
+
+   const int dim = mesh->Dimension();
+   const int num_ref_types = 3 + 4*(dim-2);
+   const int NEorig = mesh->GetNE();
+
+   aniso_flags.SetSize(NEorig);
+   error_estimates.SetSize(NEorig);
+   Vector amr_base_energy(NEorig), amr_temp_energy(NEorig);
+   error_estimates = 1.*std::numeric_limits<float>::max();
+   aniso_flags = -1;
+   GetTMOPRefinementEnergy(0, amr_base_energy);
+
+   for (int i = 1; i < num_ref_types+1; i++)
+   {
+      if ( dim == 2 && i < 3  && aniso != true ) { continue; }
+      if ( dim == 2 && i == 3 && iso   != true ) { continue; }
+      if ( dim == 3 && i < 7  && aniso != true ) { continue; }
+      if ( dim == 3 && i == 7 && iso   != true ) { continue; }
+
+      GetTMOPRefinementEnergy(i, amr_temp_energy);
+
+      for (int e = 0; e < NEorig; e++)
+      {
+         if ( amr_temp_energy(e) < error_estimates(e) )
+         {
+            error_estimates(e) = amr_temp_energy(e);
+            aniso_flags[e] = i;
+         }
+      }
+   }
+   error_estimates *= energy_scaling_factor;
+
+   if (spat_gf)
+   {
+      L2_FECollection avg_fec(0, mesh->Dimension());
+      FiniteElementSpace avg_fes(spat_gf->FESpace()->GetMesh(), &avg_fec);
+      GridFunction elem_avg(&avg_fes);
+      spat_gf->GetElementAverages(elem_avg);
+      for (int i = 0; i < amr_base_energy.Size(); i++)
+      {
+         if (elem_avg(i) < spat_gf_critical) { amr_base_energy(i) = 0.; }
+      }
+   }
+
+   error_estimates -= amr_base_energy;
+   error_estimates *= -1; // error = E(parent) - scaling_factor*mean(E(children))
+   current_sequence = mesh->GetSequence();
+}
+
+void TMOPRefinerEstimator::GetTMOPRefinementEnergy(int reftype,
+                                                   Vector &el_energy_vec)
+{
+   const FiniteElementSpace *fes = mesh->GetNodalFESpace();
+   const int NE = fes->GetNE();
+   GridFunction *xdof = mesh->GetNodes();
+   xdof->SetTrueVector();
+   xdof->SetFromTrueVector();
+
+   el_energy_vec.SetSize(NE);
+   el_energy_vec = std::numeric_limits<float>::max();
+
+   for (int e = 0; e < NE; e++)
+   {
+      Geometry::Type gtype = fes->GetFE(e)->GetGeomType();
+      DenseMatrix tr, xsplit;
+      IntegrationRule *irule = NULL;
+
+      if ( (gtype == Geometry::TRIANGLE && reftype > 0 && reftype < 3) ||
+           (gtype == Geometry::CUBE && reftype > 0 && reftype < 7) ||
+           (gtype == Geometry::TETRAHEDRON && reftype > 0 && reftype < 7) )
+      {
+         continue;
+      }
+
+      switch (gtype)
+      {
+         case Geometry::TRIANGLE:
+         {
+            int ref_access = reftype == 0 ? 0 : 1;
+            xdof->GetVectorValues(e, *TriIntRule[ref_access], xsplit, tr);
+            irule = TriIntRule[ref_access];
+            break;
+         }
+         case Geometry::TETRAHEDRON:
+         {
+            int ref_access = reftype == 0 ? 0 : 1;
+            xdof->GetVectorValues(e, *TetIntRule[ref_access], xsplit, tr);
+            irule = TetIntRule[ref_access];
+            break;
+         }
+         case Geometry::SQUARE:
+         {
+            MFEM_VERIFY(QuadIntRule[reftype], " Integration rule does not exist.");
+            xdof->GetVectorValues(e, *QuadIntRule[reftype], xsplit, tr);
+            irule = QuadIntRule[reftype];
+            break;
+         }
+         case Geometry::CUBE:
+         {
+            int ref_access = reftype == 0 ? 0 : 1;
+            xdof->GetVectorValues(e, *HexIntRule[ref_access], xsplit, tr);
+            irule = HexIntRule[ref_access];
+            break;
+         }
+         default:
+            MFEM_ABORT("Incompatible geometry type!");
+      }
+      xsplit.Transpose();
+
+      el_energy_vec(e) = 0.; // Re-set to 0
+
+      // The data format is xe1,xe2,..xen,ye1,ye2..yen.
+      // We will reformat it inside GetRefinementElementEnergy
+      Vector elfun(xsplit.GetData(), xsplit.NumCols()*xsplit.NumRows());
+
+      Array<NonlinearFormIntegrator*> &integs = *(nlf->GetDNFI());
+      TMOP_Integrator *ti  = NULL;
+      TMOPComboIntegrator *co = NULL;
+      for (int i = 0; i < integs.Size(); i++)
+      {
+         ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+         if (ti)
+         {
+            el_energy_vec(e) = ti->GetRefinementElementEnergy(*fes->GetFE(e),
+                                                              *mesh->GetElementTransformation(e),
+                                                              elfun,
+                                                              *irule);
+         }
+         co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+         if (co)
+         {
+            Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
+            for (int j = 0; j < ati.Size(); j++)
+            {
+               el_energy_vec(e) += ati[j]->GetRefinementElementEnergy(*fes->GetFE(e),
+                                                                      *mesh->GetElementTransformation(e),
+                                                                      elfun,
+                                                                      *irule);
+            }
+         }
+      }
+   }
+}
+
+void TMOPRefinerEstimator::SetHexIntRules()
+{
+   HexIntRule.SetSize(1+1);
+   // Reftype = 0 -> original element
+   Mesh meshsplit = Mesh::MakeCartesian3D(1, 1, 1, Element::HEXAHEDRON);
+   Mesh base_mesh_copy(meshsplit);
+   HexIntRule[0] = SetIntRulesFromMesh(meshsplit);
+   meshsplit.Clear();
+
+   // Reftype = 7
+   for (int i = 7; i < 8; i++)
+   {
+      Array<Refinement> marked_elements;
+      Mesh mesh_ref(base_mesh_copy);
+      for (int e = 0; e < mesh_ref.GetNE(); e++)
+      {
+         marked_elements.Append(Refinement(e, i));
+      }
+      mesh_ref.GeneralRefinement(marked_elements, 1, 0);
+      HexIntRule[1] = SetIntRulesFromMesh(mesh_ref);
+      mesh_ref.Clear();
+   }
+}
+
+void TMOPRefinerEstimator::SetQuadIntRules()
+{
+   QuadIntRule.SetSize(3+1);
+
+   // Reftype = 0 -> original element
+   Mesh meshsplit = Mesh::MakeCartesian2D(1, 1, Element::QUADRILATERAL);
+   Mesh base_mesh_copy(meshsplit);
+   QuadIntRule[0] = SetIntRulesFromMesh(meshsplit);
+   meshsplit.Clear();
+
+   // Reftype = 1-3
+   for (int i = 1; i < 4; i++)
+   {
+      Array<Refinement> marked_elements;
+      Mesh mesh_ref(base_mesh_copy);
+      for (int e = 0; e < mesh_ref.GetNE(); e++)
+      {
+         marked_elements.Append(Refinement(e, i));
+      }
+      mesh_ref.GeneralRefinement(marked_elements, 1, 0);
+      QuadIntRule[i] = SetIntRulesFromMesh(mesh_ref);
+      mesh_ref.Clear();
+   }
+}
+
+void TMOPRefinerEstimator::SetTriIntRules()
+{
+   TriIntRule.SetSize(1+1);
+
+   // Reftype = 0 // original element
+   const int Nvert = 3, NEsplit = 1;
+   Mesh meshsplit(2, Nvert, NEsplit, 0 ,2);
+   const double tri_v[3][2] =
+   {
+      {0, 0}, {1, 0}, {0, 1}
+   };
+   const int tri_e[1][3] =
+   {
+      {0, 1, 2}
+   };
+
+   for (int j = 0; j < Nvert; j++)
+   {
+      meshsplit.AddVertex(tri_v[j]);
+   }
+   meshsplit.AddTriangle(tri_e[0], 1);
+   meshsplit.FinalizeTriMesh(1, 1, true);
+
+   Mesh base_mesh_copy(meshsplit);
+   TriIntRule[0] = SetIntRulesFromMesh(meshsplit);
+   meshsplit.Clear();
+
+   // no anisotropic refinements for triangle
+   // Reftype = 3
+   for (int i = 1; i < 2; i++)
+   {
+      Array<Refinement> marked_elements;
+      Mesh mesh_ref(base_mesh_copy);
+      for (int e = 0; e < mesh_ref.GetNE(); e++)
+      {
+         marked_elements.Append(Refinement(e, i));
+      }
+      mesh_ref.GeneralRefinement(marked_elements, 1, 0);
+      TriIntRule[i] = SetIntRulesFromMesh(mesh_ref);
+      mesh_ref.Clear();
+   }
+}
+
+void TMOPRefinerEstimator::SetTetIntRules()
+{
+   TetIntRule.SetSize(1+1);
+
+   // Reftype = 0 // original element
+   const int Nvert = 4, NEsplit = 1;
+   Mesh meshsplit(3, Nvert, NEsplit, 0, 3);
+   const double tet_v[4][3] =
+   {
+      {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}
+   };
+   const int tet_e[1][4] =
+   {
+      {0, 1, 2, 3}
+   };
+
+   for (int j = 0; j < Nvert; j++)
+   {
+      meshsplit.AddVertex(tet_v[j]);
+   }
+   meshsplit.AddTet(tet_e[0], 1);
+   meshsplit.FinalizeTetMesh(1, 1, true);
+
+   Mesh base_mesh_copy(meshsplit);
+   TetIntRule[0] = SetIntRulesFromMesh(meshsplit);
+   meshsplit.Clear();
+
+   // no anisotropic refinements for triangle
+   // Reftype = 7
+   for (int i = 1; i < 2; i++)
+   {
+      Array<Refinement> marked_elements;
+      Mesh mesh_ref(base_mesh_copy);
+      for (int e = 0; e < mesh_ref.GetNE(); e++)
+      {
+         marked_elements.Append(Refinement(e, i)); //ref_type will default to 7
+      }
+      mesh_ref.GeneralRefinement(marked_elements, 1, 0);
+      TetIntRule[i] = SetIntRulesFromMesh(mesh_ref);
+      mesh_ref.Clear();
+   }
+}
+
+IntegrationRule* TMOPRefinerEstimator::SetIntRulesFromMesh(Mesh &meshsplit)
+{
+   const int dim = meshsplit.Dimension();
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace nodal_fes(&meshsplit, &fec, dim);
+   meshsplit.SetNodalFESpace(&nodal_fes);
+
+   const int NEsplit = meshsplit.GetNE();
+   const int dof_cnt = nodal_fes.GetFE(0)->GetDof(),
+             pts_cnt = NEsplit * dof_cnt;
+
+   DenseMatrix pos(dof_cnt, dim);
+   Vector posV(pos.Data(), dof_cnt * dim);
+   Array<int> xdofs(dof_cnt * dim);
+
+   // Create an IntegrationRule on the nodes of the reference submesh.
+   IntegrationRule *irule = new IntegrationRule(pts_cnt);
+   GridFunction *nodesplit = meshsplit.GetNodes();
+
+   int pt_id = 0;
+   for (int i = 0; i < NEsplit; i++)
+   {
+      nodal_fes.GetElementVDofs(i, xdofs);
+      nodesplit->GetSubVector(xdofs, posV);
+      for (int j = 0; j < dof_cnt; j++)
+      {
+         if (dim == 2)
+         {
+            irule->IntPoint(pt_id).Set2(pos(j, 0), pos(j, 1));
+         }
+         else if (dim == 3)
+         {
+            irule->IntPoint(pt_id).Set3(pos(j, 0), pos(j, 1), pos(j, 2));
+         }
+         pt_id++;
+      }
+   }
+   return irule;
+}
+
+bool TMOPDeRefinerEstimator::GetDerefineEnergyForIntegrator(
+   TMOP_Integrator &tmopi,
+   Vector &fine_energy)
+{
+   DiscreteAdaptTC *tcd = tmopi.GetDiscreteAdaptTC();
+   fine_energy.SetSize(mesh->GetNE());
+
+   if (serial)
+   {
+      Mesh meshcopy(*mesh);
+      FiniteElementSpace *tcdfes = NULL;
+      if (tcd)
+      {
+         tcdfes = new FiniteElementSpace(*tcd->GetTSpecFESpace(), &meshcopy);
+      }
+
+      Vector local_err(meshcopy.GetNE());
+      local_err = 0.;
+      double threshold = std::numeric_limits<float>::max();
+      meshcopy.DerefineByError(local_err, threshold, 0, 1);
+
+      if (meshcopy.GetGlobalNE() == mesh->GetGlobalNE())
+      {
+         delete tcdfes;
+         return false;
+      }
+
+      if (tcd)
+      {
+         tcdfes->Update();
+         tcd->SetTspecDataForDerefinement(tcdfes);
+      }
+
+      Vector coarse_energy(meshcopy.GetNE());
+      GetTMOPDerefinementEnergy(meshcopy, tmopi, coarse_energy);
+      if (tcd) { tcd->ResetDerefinementTspecData(); }
+      GetTMOPDerefinementEnergy(*mesh, tmopi, fine_energy);
+
+      const CoarseFineTransformations &dtrans =
+         meshcopy.ncmesh->GetDerefinementTransforms();
+      Table coarse_to_fine;
+      dtrans.GetCoarseToFineMap(meshcopy, coarse_to_fine);
+
+      for (int pe = 0; pe < coarse_to_fine.Size(); pe++)
+      {
+         Array<int> tabrow;
+         coarse_to_fine.GetRow(pe, tabrow);
+         int nchild = tabrow.Size();
+         double parent_energy = coarse_energy(pe);
+         for (int fe = 0; fe < nchild; fe++)
+         {
+            int child = tabrow[fe];
+            MFEM_VERIFY(child < mesh->GetNE(), " invalid coarse to fine mapping");
+            fine_energy(child) -= parent_energy;
+         }
+      }
+      delete tcdfes;
+   }
+   else
+   {
+#ifdef MFEM_USE_MPI
+      ParMesh meshcopy(*pmesh);
+      ParFiniteElementSpace *tcdfes = NULL;
+      if (tcd)
+      {
+         tcdfes = new ParFiniteElementSpace(*tcd->GetTSpecParFESpace(), meshcopy);
+      }
+
+      Vector local_err(meshcopy.GetNE());
+      local_err = 0.;
+      double threshold = std::numeric_limits<float>::max();
+      meshcopy.DerefineByError(local_err, threshold, 0, 1);
+
+      if (meshcopy.GetGlobalNE() == pmesh->GetGlobalNE())
+      {
+         delete tcdfes;
+         return false;
+      }
+
+      if (tcd)
+      {
+         tcdfes->Update();
+         tcd->SetTspecDataForDerefinement(tcdfes);
+      }
+
+      Vector coarse_energy(meshcopy.GetNE());
+      GetTMOPDerefinementEnergy(meshcopy, tmopi, coarse_energy);
+      if (tcd) { tcd->ResetDerefinementTspecData(); }
+      GetTMOPDerefinementEnergy(*pmesh, tmopi, fine_energy);
+
+      const CoarseFineTransformations &dtrans =
+         meshcopy.pncmesh->GetDerefinementTransforms();
+      Table coarse_to_fine;
+      dtrans.GetCoarseToFineMap(meshcopy, coarse_to_fine);
+
+      for (int pe = 0; pe < meshcopy.GetNE(); pe++)
+      {
+         Array<int> tabrow;
+         coarse_to_fine.GetRow(pe, tabrow);
+         int nchild = tabrow.Size();
+         double parent_energy = coarse_energy(pe);
+         for (int fe = 0; fe < nchild; fe++)
+         {
+            int child = tabrow[fe];
+            MFEM_VERIFY(child < pmesh->GetNE(), " invalid coarse to fine mapping");
+            fine_energy(child) -= parent_energy;
+         }
+      }
+      delete tcdfes;
+#endif
+   }
+
+   // error_estimate(e) = energy(parent_of_e)-energy(e)
+   // Negative energy means derefinement is desirable.
+   fine_energy *= -1;
+   return true;
+}
+
+void TMOPDeRefinerEstimator::ComputeEstimates()
+{
+   Array<NonlinearFormIntegrator*> &integs = *(nlf->GetDNFI());
+   TMOP_Integrator *ti  = NULL;
+   TMOPComboIntegrator *co = NULL;
+   error_estimates.SetSize(mesh->GetNE());
+   error_estimates = 0.;
+   Vector fine_energy(mesh->GetNE());
+
+   for (int i = 0; i < integs.Size(); i++)
+   {
+      ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+      if (ti)
+      {
+         bool deref = GetDerefineEnergyForIntegrator(*ti, fine_energy);
+         if (!deref) { error_estimates = 1; return; }
+         error_estimates += fine_energy;
+      }
+      co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+      if (co)
+      {
+         Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
+         for (int j = 0; j < ati.Size(); j++)
+         {
+            bool deref = GetDerefineEnergyForIntegrator(*ati[j], fine_energy);
+            if (!deref) { error_estimates = 1; return; }
+            error_estimates += fine_energy;
+         }
+      }
+   }
+}
+
+void TMOPDeRefinerEstimator::GetTMOPDerefinementEnergy(Mesh &cmesh,
+                                                       TMOP_Integrator &tmopi,
+                                                       Vector &el_energy_vec)
+{
+   const int cNE = cmesh.GetNE();
+   el_energy_vec.SetSize(cNE);
+   const FiniteElementSpace *fespace = cmesh.GetNodalFESpace();
+
+   GridFunction *cxdof = cmesh.GetNodes();
+
+   Array<int> vdofs;
+   Vector el_x;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+
+   for (int j = 0; j < cNE; j++)
+   {
+      fe = fespace->GetFE(j);
+      fespace->GetElementVDofs(j, vdofs);
+      T = cmesh.GetElementTransformation(j);
+      cxdof->GetSubVector(vdofs, el_x);
+      el_energy_vec(j) = tmopi.GetDerefinementElementEnergy(*fe, *T, el_x);
+   }
+}
+
+
+TMOPHRSolver::TMOPHRSolver(Mesh &mesh_, NonlinearForm &nlf_,
+                           TMOPNewtonSolver &tmopns_, GridFunction &x_,
+                           bool move_bnd_, bool hradaptivity_,
+                           int mesh_poly_deg_, int amr_metric_id_,
+                           int hr_iter_, int h_per_r_iter_) :
+   mesh(&mesh_), nlf(&nlf_), tmopns(&tmopns_), x(&x_),
+   gridfuncarr(), fespacearr(),
+   move_bnd(move_bnd_), hradaptivity(hradaptivity_),
+   mesh_poly_deg(mesh_poly_deg_), amr_metric_id(amr_metric_id_),
+   serial(true), hr_iter(hr_iter_), h_per_r_iter(h_per_r_iter_)
+{
+   if (!hradaptivity) { return; }
+   tmop_r_est = new TMOPRefinerEstimator(*mesh, *nlf, mesh_poly_deg,
+                                         amr_metric_id);
+   tmop_r = new ThresholdRefiner(*tmop_r_est);
+   tmop_r->SetTotalErrorFraction(0.0);
+   tmop_r_est->SetEnergyScalingFactor(1.);
+   tmop_dr_est= new TMOPDeRefinerEstimator(*mesh, *nlf);
+   tmop_dr = new ThresholdDerefiner(*tmop_dr_est);
+   AddGridFunctionForUpdate(x);
+}
+
+#ifdef MFEM_USE_MPI
+TMOPHRSolver::TMOPHRSolver(ParMesh &pmesh_, ParNonlinearForm &pnlf_,
+                           TMOPNewtonSolver &tmopns_, ParGridFunction &px_,
+                           bool move_bnd_, bool hradaptivity_,
+                           int mesh_poly_deg_, int amr_metric_id_,
+                           int hr_iter_, int h_per_r_iter_) :
+   mesh(&pmesh_), nlf(&pnlf_), tmopns(&tmopns_), x(&px_),
+   gridfuncarr(), fespacearr(),
+   move_bnd(move_bnd_), hradaptivity(hradaptivity_),
+   mesh_poly_deg(mesh_poly_deg_), amr_metric_id(amr_metric_id_),
+   pmesh(&pmesh_), pnlf(&pnlf_), pgridfuncarr(), pfespacearr(),
+   serial(false), hr_iter(hr_iter_), h_per_r_iter(h_per_r_iter_)
+{
+   if (!hradaptivity) { return; }
+   tmop_r_est = new TMOPRefinerEstimator(*pmesh, *pnlf, mesh_poly_deg,
+                                         amr_metric_id);
+   tmop_r = new ThresholdRefiner(*tmop_r_est);
+   tmop_r->SetTotalErrorFraction(0.0);
+   tmop_r_est->SetEnergyScalingFactor(1.);
+   tmop_dr_est= new TMOPDeRefinerEstimator(*pmesh, *pnlf);
+   tmop_dr = new ThresholdDerefiner(*tmop_dr_est);
+   AddGridFunctionForUpdate(&px_);
+}
+#endif
+
+void TMOPHRSolver::Mult()
+{
+   Vector b(0);
+   int myid = 0;
+   if (serial)
+   {
+      tmopns->SetOperator(*nlf);
+   }
+   else
+   {
+#ifdef MFEM_USE_MPI
+      myid = pnlf->ParFESpace()->GetMyRank();
+      tmopns->SetOperator(*pnlf);
+#endif
+   }
+   if (!hradaptivity)
+   {
+      tmopns->Mult(b, x->GetTrueVector());
+      if (tmopns->GetConverged() == false)
+      {
+         if (myid == 0) { mfem::out << "Nonlinear solver: rtol not achieved.\n"; }
+      }
+      x->SetFromTrueVector();
+      return;
+   }
+
+   bool radaptivity = true;
+
+   tmop_dr->Reset();
+   tmop_r->Reset();
+
+   if (serial)
+   {
+      for (int i_hr = 0; i_hr < hr_iter; i_hr++)
+      {
+         if (!radaptivity)
+         {
+            break;
+         }
+         mfem::out << i_hr << " r-adaptivity iteration.\n";
+
+         tmopns->SetOperator(*nlf);
+         tmopns->Mult(b, x->GetTrueVector());
+         x->SetFromTrueVector();
+
+         mfem::out << "TMOP energy after r-adaptivity: " <<
+                   nlf->GetGridFunctionEnergy(*x)/mesh->GetNE() <<
+                   ", Elements: " << mesh->GetNE() << std::endl;
+
+         for (int i_h = 0; i_h < h_per_r_iter; i_h++)
+         {
+            // Derefinement step.
+            if (mesh->ncmesh)
+            {
+               tmop_dr->Apply(*mesh);
+               Update();
+            }
+            mfem::out << "TMOP energy after derefinement: " <<
+                      nlf->GetGridFunctionEnergy(*x)/mesh->GetNE() <<
+                      ", Elements: " << mesh->GetNE() << std::endl;
+
+            // Refinement step.
+            tmop_r->Apply(*mesh);
+            Update();
+            mfem::out << "TMOP energy after   refinement: " <<
+                      nlf->GetGridFunctionEnergy(*x)/mesh->GetNE() <<
+                      ", Elements: " << mesh->GetNE() << std::endl;
+
+            if (!tmop_dr->Derefined() && tmop_r->Stop())
+            {
+               radaptivity = false;
+               mfem::out << "AMR stopping criterion satisfied. Stop.\n";
+               break;
+            }
+         } //n_h
+      } //n_hr
+   }
+   else
+   {
+#ifdef MFEM_USE_MPI
+      int NEGlob;
+      double tmopenergy;
+      for (int i_hr = 0; i_hr < hr_iter; i_hr++)
+      {
+         if (!radaptivity)
+         {
+            break;
+         }
+         if (myid == 0) { mfem::out << i_hr << " r-adaptivity iteration.\n"; }
+         tmopns->SetOperator(*pnlf);
+         tmopns->Mult(b, x->GetTrueVector());
+         x->SetFromTrueVector();
+
+         NEGlob = pmesh->GetGlobalNE();
+         tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
+         if (myid == 0)
+         {
+            mfem::out << "TMOP energy after r-adaptivity: " << tmopenergy <<
+                      ", Elements: " << NEGlob << std::endl;
+         }
+
+         for (int i_h = 0; i_h < h_per_r_iter; i_h++)
+         {
+            // Derefinement step.
+            if (pmesh->pncmesh)
+            {
+               RebalanceParNCMesh();
+               ParUpdate();
+
+               tmop_dr->Apply(*pmesh);
+               ParUpdate();
+            }
+            NEGlob = pmesh->GetGlobalNE();
+            tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
+            if (myid == 0)
+            {
+               mfem::out << "TMOP energy after derefinement: " << tmopenergy <<
+                         ", Elements: " << NEGlob << std::endl;
+            }
+
+            // Refinement step.
+            tmop_r->Apply(*pmesh);
+            ParUpdate();
+            NEGlob = pmesh->GetGlobalNE();
+            tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
+            if (myid == 0)
+            {
+               mfem::out << "TMOP energy after   refinement: " << tmopenergy <<
+                         ", Elements: " << NEGlob << std::endl;
+            }
+
+            if (!tmop_dr->Derefined() && tmop_r->Stop())
+            {
+               radaptivity = false;
+               if (myid == 0)
+               {
+                  mfem::out << "AMR stopping criterion satisfied. Stop.\n";
+               }
+               break;
+            }
+         } //n_r limit
+      } //n_hr
+#endif
+   }
+}
+
+#ifdef MFEM_USE_MPI
+void TMOPHRSolver::RebalanceParNCMesh()
+{
+   ParNCMesh *pncmesh = pmesh->pncmesh;
+   if (pncmesh)
+   {
+      const Table &dreftable = pncmesh->GetDerefinementTable();
+      Array<int> drefs, new_ranks;
+      for (int i = 0; i < dreftable.Size(); i++)
+      {
+         drefs.Append(i);
+      }
+      pncmesh->GetFineToCoarsePartitioning(drefs, new_ranks);
+      pmesh->Rebalance(new_ranks);
+   }
+}
+#endif
+
+void TMOPHRSolver::Update()
+{
+   // Update FESpace
+   for (int i = 0; i < fespacearr.Size(); i++)
+   {
+      fespacearr[i]->Update();
+   }
+   // Update nodal GF
+   for (int i = 0; i < gridfuncarr.Size(); i++)
+   {
+      gridfuncarr[i]->Update();
+      gridfuncarr[i]->SetTrueVector();
+      gridfuncarr[i]->SetFromTrueVector();
+   }
+
+   // Update Discrete Indicator for all the TMOP_Integrators in NonLinearForm
+   Array<NonlinearFormIntegrator*> &integs = *(nlf->GetDNFI());
+   TMOP_Integrator *ti  = NULL;
+   TMOPComboIntegrator *co = NULL;
+   DiscreteAdaptTC *dtc = NULL;
+   for (int i = 0; i < integs.Size(); i++)
+   {
+      ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+      if (ti)
+      {
+         ti->UpdateAfterMeshTopologyChange();
+         dtc = ti->GetDiscreteAdaptTC();
+         if (dtc) { dtc->UpdateAfterMeshTopologyChange(); }
+      }
+      co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+      if (co)
+      {
+         Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
+         for (int j = 0; j < ati.Size(); j++)
+         {
+            ati[j]->UpdateAfterMeshTopologyChange();
+            dtc = ati[j]->GetDiscreteAdaptTC();
+            if (dtc) { dtc->UpdateAfterMeshTopologyChange(); }
+         }
+      }
+   }
+
+   // Update the Nonlinear form and set Essential BC.
+   UpdateNonlinearFormAndBC(mesh, nlf);
+}
+
+#ifdef MFEM_USE_MPI
+void TMOPHRSolver::ParUpdate()
+{
+   // Update FESpace
+   for (int i = 0; i < pfespacearr.Size(); i++)
+   {
+      pfespacearr[i]->Update();
+   }
+   // Update nodal GF
+   for (int i = 0; i < pgridfuncarr.Size(); i++)
+   {
+      pgridfuncarr[i]->Update();
+      pgridfuncarr[i]->SetTrueVector();
+      pgridfuncarr[i]->SetFromTrueVector();
+   }
+
+   // Update Discrete Indicator
+   Array<NonlinearFormIntegrator*> &integs = *(nlf->GetDNFI());
+   TMOP_Integrator *ti  = NULL;
+   TMOPComboIntegrator *co = NULL;
+   DiscreteAdaptTC *dtc = NULL;
+   for (int i = 0; i < integs.Size(); i++)
+   {
+      ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+      if (ti)
+      {
+         ti->ParUpdateAfterMeshTopologyChange();
+         dtc = ti->GetDiscreteAdaptTC();
+         if (dtc) { dtc->ParUpdateAfterMeshTopologyChange(); }
+      }
+      co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+      if (co)
+      {
+         Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
+         for (int j = 0; j < ati.Size(); j++)
+         {
+            ati[j]->ParUpdateAfterMeshTopologyChange();
+            dtc = ati[j]->GetDiscreteAdaptTC();
+            if (dtc) { dtc->ParUpdateAfterMeshTopologyChange(); }
+         }
+      }
+   }
+
+   // Update the Nonlinear form and set Essential BC.
+   UpdateNonlinearFormAndBC(pmesh, pnlf);
+}
+#endif
+
+void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
+{
+   const FiniteElementSpace &fes = *mesh->GetNodalFESpace();
+
+   // Update Nonlinear form and Set Essential BC
+   nlf->Update();
+   const int dim = fes.GetFE(0)->GetDim();
+   if (move_bnd == false)
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      nlf->SetEssentialBC(ess_bdr);
+   }
+   else
+   {
+      const int nd  = fes.GetBE(0)->GetDof();
+      int n = 0;
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         const int attr = mesh->GetBdrElement(i)->GetAttribute();
+         MFEM_VERIFY(!(dim == 2 && attr == 3),
+                     "Boundary attribute 3 must be used only for 3D meshes. "
+                     "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
+                     "components, rest for free nodes), or use -fix-bnd.");
+         if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
+         if (attr == 4) { n += nd * dim; }
+      }
+      Array<int> ess_vdofs(n), vdofs;
+      n = 0;
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         const int attr = mesh->GetBdrElement(i)->GetAttribute();
+         fes.GetBdrElementVDofs(i, vdofs);
+         if (attr == 1) // Fix x components.
+         {
+            for (int j = 0; j < nd; j++)
+            { ess_vdofs[n++] = vdofs[j]; }
+         }
+         else if (attr == 2) // Fix y components.
+         {
+            for (int j = 0; j < nd; j++)
+            { ess_vdofs[n++] = vdofs[j+nd]; }
+         }
+         else if (attr == 3) // Fix z components.
+         {
+            for (int j = 0; j < nd; j++)
+            { ess_vdofs[n++] = vdofs[j+2*nd]; }
+         }
+         else if (attr == 4) // Fix all components.
+         {
+            for (int j = 0; j < vdofs.Size(); j++)
+            { ess_vdofs[n++] = vdofs[j]; }
+         }
+      }
+      nlf->SetEssentialVDofs(ess_vdofs);
+   }
+}
+
+}

fem/tmop_amr.hpp

@@ -0,0 +1,284 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_TMOP_AMR_HPP
+#define MFEM_TMOP_AMR_HPP
+
+#include "tmop_tools.hpp"
+#include "nonlinearform.hpp"
+#include "pnonlinearform.hpp"
+#include "estimators.hpp"
+#include "../mesh/mesh_operators.hpp"
+
+namespace mfem
+{
+
+class TMOPRefinerEstimator : public AnisotropicErrorEstimator
+{
+protected:
+   Mesh *mesh; // not owned
+   NonlinearForm *nlf; // not owned
+   int order;
+   int amrmetric;
+   Array<IntegrationRule *> TriIntRule, QuadIntRule, TetIntRule, HexIntRule;
+   long current_sequence;
+   Vector error_estimates;
+   Array<int> aniso_flags;
+   // An element is refined only if
+   // [mean TMOPEnergy(children)]*energy_scaling_factor < TMOPEnergy(parent)
+   double energy_scaling_factor;
+   GridFunction *spat_gf;   // If specified, can be used to specify the
+   double spat_gf_critical; // region where hr-adaptivity is done.
+
+   /// Check if the mesh of the solution was modified.
+   bool MeshIsModified()
+   {
+      long mesh_sequence = mesh->GetSequence();
+      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
+      return (mesh_sequence > current_sequence);
+   }
+
+   /// Compute the element error estimates. For an element E in the mesh,
+   /// error(E) = TMOPEnergy(E)*energy_scaling_factor-Mean(TMOPEnergy(ChildofE)),
+   /// where TMOPEnergy of Children of E is obtained by assuming the element E
+   /// is refined using the refinement type being considered based on the TMOP
+   /// mesh quality metric.
+   void ComputeEstimates();
+
+   /// Construct the integration rules to model how each element type is split
+   /// using different refinement types. ref_type = 0 is the original element
+   /// and reftype \ in [1, 7] represent different refinement type based on
+   /// NCMesh class.
+   void SetQuadIntRules(); // supports ref_type = 1 to 3.
+   void SetTriIntRules(); // currently supports only isotropic refinement.
+   void SetHexIntRules(); // currently supports only isotropic refinement.
+   void SetTetIntRules(); // currently supports only isotropic refinement.
+
+   /// Get TMOP energy for each element corresponding to the refinement type
+   /// specified.
+   void GetTMOPRefinementEnergy(int reftype, Vector &el_energy_vec);
+
+   /// Use a mesh to setup an integration rule that will mimic the different
+   /// refinement types.
+   IntegrationRule* SetIntRulesFromMesh(Mesh &meshsplit);
+public:
+   TMOPRefinerEstimator(Mesh &mesh_, NonlinearForm &nlf_, int order_,
+                        int amrmetric_) :
+      mesh(&mesh_), nlf(&nlf_), order(order_), amrmetric(amrmetric_),
+      TriIntRule(0), QuadIntRule(0), TetIntRule(0), HexIntRule(0),
+      current_sequence(-1), error_estimates(), aniso_flags(),
+      energy_scaling_factor(1.), spat_gf(NULL), spat_gf_critical(0.)
+   {
+      if (mesh->Dimension() == 2)
+      {
+         SetQuadIntRules();
+         SetTriIntRules();
+      }
+      else
+      {
+         SetHexIntRules();
+         SetTetIntRules();
+      }
+   }
+
+   ~TMOPRefinerEstimator()
+   {
+      for (int i = 0; i < QuadIntRule.Size(); i++) { delete QuadIntRule[i]; }
+      for (int i = 0; i < TriIntRule.Size();  i++) { delete TriIntRule[i]; }
+      for (int i = 0; i < HexIntRule.Size();  i++) { delete HexIntRule[i]; }
+      for (int i = 0; i < TetIntRule.Size();  i++) { delete TetIntRule[i]; }
+   }
+
+   /// Get TMOP-based errors for each element in the mesh computed based on the
+   /// refinement types being considered.
+   virtual const Vector &GetLocalErrors()
+   {
+      if (MeshIsModified()) { ComputeEstimates(); }
+      return error_estimates;
+   }
+   /// For anisotropic refinements, get the refinement type (e.g., x or y)
+   virtual const Array<int> &GetAnisotropicFlags()
+   {
+      if (MeshIsModified()) { ComputeEstimates(); }
+      return aniso_flags;
+   }
+
+   /// Scaling factor for the TMOP refinement energy. An element is refined if
+   /// [mean TMOPEnergy(children)]*energy_scaling_factor < TMOPEnergy(parent)
+   void SetEnergyScalingFactor(double scale) { energy_scaling_factor = scale; }
+
+   /// Spatial indicator function (eta) that can be used to prevent elements
+   /// from being refined even if the energy criterion is met. Using this,
+   /// an element E is not refined if mean(@a spat_gf(E)) < @a spat_gf_critical.
+   void SetSpatialIndicator(GridFunction &spat_gf_,
+                            double spat_gf_critical_ = 0.5)
+   { spat_gf = &spat_gf_; spat_gf_critical = spat_gf_critical_; }
+   void SetSpatialIndicatorCritical(double val_) { spat_gf_critical = val_; }
+
+   /// Reset the error estimator.
+   virtual void Reset() { current_sequence = -1; }
+};
+
+class TMOPDeRefinerEstimator : public ErrorEstimator
+{
+protected:
+   Mesh *mesh;
+   NonlinearForm *nlf;
+#ifdef MFEM_USE_MPI
+   ParMesh *pmesh;
+   ParNonlinearForm *pnlf;
+#endif
+   int order;
+   int amrmetric;
+   long current_sequence;
+   Vector error_estimates;
+   bool serial;
+
+   /// Check if the mesh of the solution was modified.
+   bool MeshIsModified()
+   {
+      long mesh_sequence = mesh->GetSequence();
+      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
+      return (mesh_sequence > current_sequence);
+   }
+
+   /// Compute the element error estimates. For a given element E in the mesh,
+   /// error(E) = TMOPEnergy(parent_of_E)-TMOPEnergy(E). Children element of an
+   /// element are derefined if the mean TMOP energy of children is greated than
+   /// the TMOP energy associated with their parent.
+   void ComputeEstimates();
+
+   void GetTMOPDerefinementEnergy(Mesh &cmesh,
+                                  TMOP_Integrator &tmopi,
+                                  Vector &el_energy_vec);
+
+   bool GetDerefineEnergyForIntegrator(TMOP_Integrator &tmopi,
+                                       Vector &fine_energy);
+public:
+   TMOPDeRefinerEstimator(Mesh &mesh_, NonlinearForm &nlf_) :
+      mesh(&mesh_), nlf(&nlf_),
+      current_sequence(-1), error_estimates(), serial(true)   { }
+#ifdef MFEM_USE_MPI
+   TMOPDeRefinerEstimator(ParMesh &pmesh_, ParNonlinearForm &pnlf_) :
+      mesh(&pmesh_), nlf(&pnlf_), pmesh(&pmesh_), pnlf(&pnlf_),
+      current_sequence(-1), error_estimates(), serial(false) { }
+#endif
+
+   ~TMOPDeRefinerEstimator() { }
+
+   virtual const Vector &GetLocalErrors()
+   {
+      if (MeshIsModified()) { ComputeEstimates(); }
+      return error_estimates;
+   }
+
+   /// Reset the error estimator.
+   virtual void Reset() { current_sequence = -1; }
+};
+
+// hr-adaptivity using TMOP.
+// If hr-adaptivity is disabled, r-adaptivity is done once using the
+// TMOPNewtonSolver.
+// Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
+// "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity iteration.
+// The solver terminates early if an h-adaptivity iteration does not
+// refine/derefine any element in the mesh.
+class TMOPHRSolver
+{
+protected:
+   Mesh *mesh;
+   NonlinearForm *nlf;
+   TMOPNewtonSolver *tmopns;
+   GridFunction *x;
+   Array<GridFunction *> gridfuncarr;
+   Array<FiniteElementSpace *> fespacearr;
+   bool move_bnd, hradaptivity;
+   const int mesh_poly_deg, amr_metric_id;
+#ifdef MFEM_USE_MPI
+   ParMesh *pmesh;
+   ParNonlinearForm *pnlf;
+   Array<ParGridFunction *> pgridfuncarr;
+   Array<ParFiniteElementSpace *> pfespacearr;
+#endif
+   bool serial;
+
+   // All are owned.
+   TMOPRefinerEstimator *tmop_r_est;
+   ThresholdRefiner *tmop_r;
+   TMOPDeRefinerEstimator *tmop_dr_est;
+   ThresholdDerefiner *tmop_dr;
+
+   int hr_iter, h_per_r_iter;
+
+   void Update();
+#ifdef MFEM_USE_MPI
+   void ParUpdate();
+#endif
+   void UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf);
+
+#ifdef MFEM_USE_MPI
+   // Rebalance ParMesh such that all the children elements are moved to the same
+   // MPI rank where the parent will be if the mesh were to be derefined.
+   void RebalanceParNCMesh();
+#endif
+
+public:
+   TMOPHRSolver(Mesh &mesh_, NonlinearForm &nlf_,
+                TMOPNewtonSolver &tmopns_, GridFunction &x_,
+                bool move_bnd_, bool hradaptivity_,
+                int mesh_poly_deg_, int amr_metric_id_,
+                int hr_iter_ = 5, int h_per_r_iter_ = 1);
+#ifdef MFEM_USE_MPI
+   TMOPHRSolver(ParMesh &pmesh_, ParNonlinearForm &pnlf_,
+                TMOPNewtonSolver &tmopns_, ParGridFunction &x_,
+                bool move_bnd_, bool hradaptivity_,
+                int mesh_poly_deg_, int amr_metric_id_,
+                int hr_iter_ = 5, int h_per_r_iter_ = 1);
+#endif
+
+   void Mult();
+
+   /// These are used to update spaces and functions that are not owned by the
+   /// TMOPIntegrator or DiscreteAdaptTC. The owned ones are updated in the
+   /// functions UpdateAfterMeshTopologyChange() of both classes.
+   void AddGridFunctionForUpdate(GridFunction *gf) { gridfuncarr.Append(gf); }
+   void AddFESpaceForUpdate(FiniteElementSpace *fes) { fespacearr.Append(fes); }
+
+#ifdef MFEM_USE_MPI
+   void AddGridFunctionForUpdate(ParGridFunction *pgf_)
+   {
+      pgridfuncarr.Append(pgf_);
+   }
+   void AddFESpaceForUpdate(ParFiniteElementSpace *pfes_)
+   {
+      pfespacearr.Append(pfes_);
+   }
+#endif
+
+   ~TMOPHRSolver()
+   {
+      if (!hradaptivity) { return; }
+      delete tmop_dr;
+      delete tmop_dr_est;
+      delete tmop_r;
+      delete tmop_r_est;
+   }
+
+   /// Total number of hr-adaptivity iterations. At each iteration, we do an
+   /// r-adaptivity iteration followed by a number of h-adaptivity iterations.
+   void SetHRAdaptivityIterations(int iter) { hr_iter = iter; }
+
+   /// Total number of h-adaptivity iterations per r-adaptivity iteration.
+   void SetHAdaptivityIterations(int iter) { h_per_r_iter = iter; }
+};
+
+}
+#endif

fem/tmop_tools.cpp

@@ -407,6 +407,8 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
    {
       // Needed for the line search below. The untangling metrics see this
       // reference to detect deteriorations.
+      MFEM_VERIFY(min_det_ptr != NULL, " Initial mesh was valid, but"
+                  " intermediate mesh is invalid. Contact TMOP Developers.");
       *min_det_ptr = untangle_factor * min_detT_in;
    }
 
@@ -576,7 +578,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
          ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
          if (ti)
          {
-            ti->UpdateAfterMeshChange(x_loc);
+            ti->UpdateAfterMeshPositionChange(x_loc);
             ti->ComputeFDh(x_loc, *pfesc);
             UpdateDiscreteTC(*ti, x_loc);
          }
@@ -586,7 +588,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
             Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
             for (int j = 0; j < ati.Size(); j++)
             {
-               ati[j]->UpdateAfterMeshChange(x_loc);
+               ati[j]->UpdateAfterMeshPositionChange(x_loc);
                ati[j]->ComputeFDh(x_loc, *pfesc);
                UpdateDiscreteTC(*ati[j], x_loc);
             }
@@ -613,7 +615,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
          ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
          if (ti)
          {
-            ti->UpdateAfterMeshChange(x_loc);
+            ti->UpdateAfterMeshPositionChange(x_loc);
             ti->ComputeFDh(x_loc, *fesc);
             UpdateDiscreteTC(*ti, x_loc);
          }
@@ -623,7 +625,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
             Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
             for (int j = 0; j < ati.Size(); j++)
             {
-               ati[j]->UpdateAfterMeshChange(x_loc);
+               ati[j]->UpdateAfterMeshPositionChange(x_loc);
                ati[j]->ComputeFDh(x_loc, *fesc);
                UpdateDiscreteTC(*ati[j], x_loc);
             }

general/mem_manager.hpp

@@ -494,8 +494,7 @@ public:
    /// Copy @a size entries from @a *this to @a dest.
    /** The given @a size should not exceed the Capacity() of @a *this and the
        destination, @a dest. */
-   inline void CopyTo(Memory &dest, int size) const
-   { dest.CopyFrom(*this, size); }
+   inline void CopyTo(Memory &dest, int size) const;
 
    /// 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. */
@@ -923,6 +922,11 @@ inline void Memory<T>::Wrap(T *ptr, T *d_ptr, int size, MemoryType mt, bool own)
 template <typename T>
 inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
 {
+   MFEM_ASSERT(0 <= offset, "invalid offset = " << offset);
+   MFEM_ASSERT(0 <= size, "invalid size = " << size);
+   MFEM_ASSERT(offset + size <= base.capacity,
+               "invalid offset + size = " << offset + size
+               << " > base capacity = " << base.capacity);
    capacity = size;
    h_mt = base.h_mt;
    h_ptr = base.h_ptr + offset;
@@ -1136,6 +1140,7 @@ inline bool Memory<T>::DeviceIsValid() const
 template <typename T>
 inline void Memory<T>::CopyFrom(const Memory &src, int size)
 {
+   MFEM_VERIFY(src.capacity>=size && capacity>=size, "Incorrect size");
    if (!(flags & REGISTERED) && !(src.flags & REGISTERED))
    {
       if (h_ptr != src.h_ptr && size != 0)
@@ -1155,6 +1160,7 @@ inline void Memory<T>::CopyFrom(const Memory &src, int size)
 template <typename T>
 inline void Memory<T>::CopyFromHost(const T *src, int size)
 {
+   MFEM_VERIFY(capacity>=size, "Incorrect size");
    if (!(flags & REGISTERED))
    {
       if (h_ptr != src && size != 0)
@@ -1171,9 +1177,17 @@ inline void Memory<T>::CopyFromHost(const T *src, int size)
    }
 }
 
+template <typename T>
+inline void Memory<T>::CopyTo(Memory &dest, int size) const
+{
+   MFEM_VERIFY(capacity>=size, "Incorrect size");
+   dest.CopyFrom(*this, size);
+}
+
 template <typename T>
 inline void Memory<T>::CopyToHost(T *dest, int size) const
 {
+   MFEM_VERIFY(capacity>=size, "Incorrect size");
    if (!(flags & REGISTERED))
    {
       if (h_ptr != dest && size != 0)

linalg/imex.hpp

@@ -0,0 +1,244 @@
+#ifndef MFEM_IMEX
+#define MFEM_IMEX
+
+#include "../config/config.hpp"
+#include "operator.hpp"
+#include "ode.hpp"
+#include <vector>
+
+namespace mfem
+{
+
+/** Class for spatial discretizations of a PDE resulting in the time-dependent, 
+    nonlinear set of ODEs with implicit-explicit additive partition
+        M*du/dt = N_E(u,t) + N_I(u,t).
+    MFEM typically treats time integration as 
+        du/dt = F^{-1} G(u),
+    Here F represents what MFEM calls the “implicit” part, and G represents the
+    “explicit” part; in simpler terms, F is typically just a mass matrix. 
+
+    For BDF schemes, the ImplicitSolve function is a bit different, and it is
+    more natural to apply M and M^{-1} separate from the Mult functions, so we
+    include MassMult and MassInv as functions to be provided, and do not include
+    such actions in the Mult functions. */
+class IMEXTimeDependentOperator : public TimeDependentOperator
+{
+protected:
+    mutable Vector temp;    // Auxillary vector
+    
+public:
+    // Sets linearly implicit to false by default
+    IMEXTimeDependentOperator(int n, double t=0.0, Type type=EXPLICIT)
+        : TimeDependentOperator(n, t, type) { };
+    ~IMEXTimeDependentOperator() { };
+
+    /** Apply action of implicit part of operator y <- N_I(x,y). For fully
+        implicit schemes, this just corresponds to applying the time-dependent
+        (nonlinear) operator.
+        PREVIOUSLY CALLED ExplicitMult */
+    virtual void ImplicitMult(const Vector &x, Vector &y) const = 0;
+
+    /** Apply action of explicit part of operator y <- N_E(x,y) */
+    virtual void ExplicitMult(const Vector &x, Vector &y) const { y = 0.0; };
+
+    /** Solve k = f(x+dt*k) for stage k, where f() is the implicit part of
+        the operator. Used in Runge-Kutta methods. */
+    virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
+    { mfem::mfem_error("IMEXTimeDependentOperator::ImplicitSolve() is not overridden!"); };
+
+    /** Solve M*x - dtf(x, t) = b for solution x, where f() is the implicit
+        part of the operator. Used in BDF methods. */
+    virtual void ImplicitSolve2(const double dt, const Vector &b, Vector &x)
+    { mfem::mfem_error("IMEXTimeDependentOperator::ImplicitSolve2() is not overridden!"); };
+
+    /** Apply action mass matrix, y = M*x. 
+        If not re-implemented, this method simply generates an error.
+        PREVIOUSLY CALLED ImplictMult */
+    virtual void MassMult(const Vector &x, Vector &y) const = 0;
+
+    /** Apply action of inverse of mass matrix, y = M^{-1}*x. 
+        If not re-implemented, this method simply generates an error.
+        NOTE : only necessary for PolyIMEX methods. */
+    virtual void MassInv(const Vector &x, Vector &y) const = 0;
+};
+
+/** Class holding RK Butcher tableau, and associated data required by
+    implicit and explicit splitting. */
+class IMEXRKData 
+{
+public:
+    // Implicit Runge Kutta type. Enumeration (s, \sigma, p):
+    // - s = number of implicit stages
+    // - \sigma = number of explicit stages
+    // - p = order
+    // In this notation, when s = \sigma, we satisfy (2.3)/(2.4) in
+    // Ascher et al., and do not need to compute the final explicit
+    // stage. This is represented in the stiffly_accurate boolean.
+    enum Type { 
+        IMEX111 = 111,
+        IMEX121 = 121,
+        IMEX122 = 122,
+        IMEX222 = 222,
+        IMEX232 = 232,
+        IMEX233 = 233,
+        IMEX443 = 443,
+    // ARK ESDIRK-ERK schemes: enumeration (s,p), for total number of
+    // stages s.
+        ARK43 = -43
+    };
+
+    IMEXRKData() : s(-1) { };
+    IMEXRKData(Type ID_) : ID(ID_) { SetData(); };
+    ~IMEXRKData() { };
+    
+    /// Set explicit RK data
+    void SetExplicitData(DenseMatrix Ae_, Vector be_, Vector ce_);
+    /// Set implicit RK data
+    void SetImplicitData(DenseMatrix Ai_, Vector bi_, Vector ci_, bool esdirk_=false);
+    void SetID(Type ID_) { ID=ID_; SetData(); };
+
+    bool esdirk;
+    bool stiffly_accurate;
+    bool use_final_exp_stage;
+    int s;
+
+    DenseMatrix Ai;     // Implicit Butcher matrix
+    Vector bi;          // Implicit Butcher tableau weights
+    DenseMatrix Ae;     // Explicit Butcher matrix
+    Vector be;          // Explicit Butcher tableau weights
+    Vector c0;          // Butcher tableau nodes (same for implicit and explicit!)
+
+private:    
+    Type ID;
+    void SetData();
+    void InitData();
+};
+
+/** Class for two-part additive IMEX RK method, where explicit and implicit
+    stage vectors are stored. Assume same abscissae, {c}, for both schemes.
+    Butcher Data must be provided either in a custom IMEXRKData object, or
+    using the IMEXRKData::Type for predefined tableaux. */
+class IMEXRK : public ODESolver
+{
+protected:
+    IMEXRKData tableaux;
+    std::vector< Vector *> exp_stages;
+    std::vector< Vector *> imp_stages;
+    IMEXTimeDependentOperator *imex;    // Spatial discretization. 
+
+public:
+    IMEXRK(IMEXRKData tableaux_) : ODESolver(), tableaux(tableaux_) { };
+    IMEXRK(IMEXRKData::Type type_) : ODESolver(), tableaux(type_) { };
+    ~IMEXRK();
+    void Init(IMEXTimeDependentOperator &_imex);
+    void Step(Vector &x, double &t, double &dt) override;
+};
+
+/** Class holding BDF integrator data. Setting alpha < 0 (the default
+    constructor) defines alpha = 2/(q-1), corresponding to classical BDF
+    of order q. */
+class BDFData 
+{
+public:
+
+    enum Type {
+        BDF1 = 01, BDF2 = 02, BDF3 = 03, BDF4 = 04,
+        IMEX_BDF1 = 11, IMEX_BDF2 = 12, IMEX_BDF3 = 13,
+        IMEX_BDF4 = 14
+    };
+
+    BDFData() { };
+    BDFData(Type ID_, double alpha_=-1) : ID(ID_), alpha(alpha_) { SetData(); };
+    ~BDFData() { };
+
+    int GetID() { return static_cast<int>(ID); };
+    void SetID(Type ID_, double alpha_=-1) 
+    {
+        ID=ID_;
+        alpha = alpha_;
+        SetData();
+    };
+    void Print()
+    {
+        std::cout << "q     = " << q << "\n";
+        std::cout << "alpha = " << alpha << "\n";
+        std::cout << "A:\n";
+        A.PrintMatlab();
+        std::cout << "Be:\n";
+        Be.PrintMatlab();
+        std::cout << "Bi:\n";
+        Bi.Print();
+        std::cout << "z:\n";
+        z0.Print();
+    };
+
+    double alpha;
+    int q;              // Number of previous values stored
+    bool shifted_nodes; // false = clssical BDF, true = Polynomial BDF w/ shifted nodes
+    DenseMatrix A;      // Previous solution coefficients
+    Vector Bi;          // Implicit coefficients
+    DenseMatrix Be;     // Explicit coefficients
+    Vector z0;
+
+
+private:    
+    Type ID;
+    void SetData();
+    void InitData();
+};
+
+/** Class for IMEX-BDF methods, including classical IMEX-BDF and IMEX-
+    Polynomial-BDF (IMEX-PBDF). IMEX-PBDF methods have an additional 
+    alpha parameter, where larger alpha leads to smaller stability
+    regions and a smaller leading accuracy constant, while smaller
+    alpha leads to larger stabiltiy regions and a larger accuracy
+    constant. For classical methods, there are two implementations:
+        - ClassicalStep() stores previous solutions and the explicit
+        part of the operator evaluated on the solution, and
+        - ClassicalStepNoStore() does not store the explicit 
+        component, but must re-evaluate q times during each time
+        step.
+    This option can be set via the recompute_exp input. The type of
+    scheme must be set through the BDFData structure or BDFData::Type.
+        There is also an option to use pointwise Lagrange interpolating 
+    polynomials to provide an initial guess for the ImplicitSolve. This
+    is set via InterpolateGuess(). This option is only implemented for
+    PBDF. */
+class IMEXBDF : public ODESolver
+{
+private:
+    BDFData data;
+    bool recompute_exp;
+    bool interpolate;
+    int initialized;
+    double dt_prev;
+    std::vector< Vector*> sols;
+    std::vector< Vector*> exp_sols;
+    IMEXTimeDependentOperator *imex;    // Spatial discretization
+    IMEXRK *RKsolver;
+    std::vector<double> exp_nodes;
+
+    void AlphaStep(Vector &x, double &t, double &dt);
+    void ClassicalStep(Vector &x, double &t, double &dt);
+    void ClassicalStepNoStore(Vector &x, double &t, double &dt);
+
+public:
+    IMEXBDF(BDFData data_, bool recompute_exp_=false) :
+      ODESolver(), data(data_), recompute_exp(recompute_exp_),
+      interpolate(false) { };
+    IMEXBDF(BDFData::Type scheme, bool recompute_exp_=false) :
+      ODESolver(), recompute_exp(recompute_exp_), interpolate(false)
+      { data.SetID(scheme); };
+    IMEXBDF(BDFData::Type scheme, double alpha) :
+      ODESolver(), interpolate(false), recompute_exp(false)
+      { data.SetID(scheme, alpha); };
+    ~IMEXBDF();
+
+    void Init(IMEXTimeDependentOperator &_imex);
+    void Step(Vector &x, double &t, double &dt);
+    void InterpolateGuess() {interpolate = true; };
+};
+
+}
+
+#endif

linalg/linalg.hpp

@@ -31,6 +31,7 @@
 #include "invariants.hpp"
 #include "constraints.hpp"
 #include "auxiliary.hpp"
+#include "imex.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"

linalg/sparsemat.cpp

@@ -33,7 +33,12 @@ int SparseMatrix::SparseMatrixCount = 0;
 cusparseHandle_t SparseMatrix::handle = nullptr;
 size_t SparseMatrix::bufferSize = 0;
 void * SparseMatrix::dBuffer = nullptr;
-#endif
+#  if CUSPARSE_VERSION >=  11400
+#    define MFEM_CUSPARSE_ALG CUSPARSE_SPMV_CSR_ALG1
+#  else
+#    define MFEM_CUSPARSE_ALG CUSPARSE_CSRMV_ALG1
+#  endif // CUSPARSE_VERSION >= 11400
+#endif // MFEM_USE_CUDA
 
 void SparseMatrix::InitCuSparse()
 {
@@ -679,25 +684,16 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
          cusparseCreateMatDescr(&matA_descr);
          cusparseSetMatIndexBase(matA_descr, CUSPARSE_INDEX_BASE_ZERO);
          cusparseSetMatType(matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL);
-
 #endif
-
          initBuffers = true;
       }
       // Allocate kernel space. Buffer is shared between different sparsemats
       size_t newBufferSize = 0;
 
-#if CUDA_VERSION >= 11020
       cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
                               matA_descr,
                               vecX_descr, &beta, vecY_descr, CUDA_R_64F,
-                              CUSPARSE_SPMV_CSR_ALG1, &newBufferSize);
-#elif CUDA_VERSION >= 10010
-      cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
-                              matA_descr,
-                              vecX_descr, &beta, vecY_descr, CUDA_R_64F,
-                              CUSPARSE_CSRMV_ALG1, &newBufferSize);
-#endif
+                              MFEM_CUSPARSE_ALG, &newBufferSize);
 
       // Check if we need to resize
       if (newBufferSize > bufferSize)
@@ -707,30 +703,22 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
          CuMemAlloc(&dBuffer, bufferSize);
       }
 
-#if CUDA_VERSION >= 11020
+#if CUDA_VERSION >= 10010
       // Update input/output vectors
       cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
       cusparseDnVecSetValues(vecY_descr, d_y);
 
       // Y = alpha A * X + beta * Y
       cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA_descr,
-                   vecX_descr, &beta, vecY_descr, CUDA_R_64F, CUSPARSE_SPMV_CSR_ALG1, dBuffer);
-#elif CUDA_VERSION >= 10010
-      // Update input/output vectors
-      cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
-      cusparseDnVecSetValues(vecY_descr, d_y);
-
-      // Y = alpha A * X + beta * Y
-      cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA_descr,
-                   vecX_descr, &beta, vecY_descr, CUDA_R_64F, CUSPARSE_CSRMV_ALG1, dBuffer);
+                   vecX_descr, &beta, vecY_descr, CUDA_R_64F, MFEM_CUSPARSE_ALG, dBuffer);
 #else
       cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                      Height(), Width(), J.Capacity(),
                      &alpha, matA_descr,
                      const_cast<double *>(d_A), const_cast<int *>(d_I), const_cast<int *>(d_J),
                      const_cast<double *>(d_x), &beta, d_y);
-#endif
-#endif
+#endif // CUDA_VERSION >= 10010
+#endif // MFEM_USE_CUDA
    }
    else
    {

linalg/vector.hpp

@@ -82,6 +82,11 @@ public:
    Vector(double *data_, int size_)
    { data.Wrap(data_, size_, false); size = size_; }
 
+   /** @brief Create a Vector referencing a sub-vector of the Vector @a base
+       starting at the given offset, @a base_offset, and size @a size_. */
+   Vector(Vector &base, int base_offset, int size_)
+      : data(base.data, base_offset, size_), size(size_) { }
+
    /// Create a Vector of size @a size_ using MemoryType @a mt.
    Vector(int size_, MemoryType mt)
       : data(size_, mt), size(size_) { }

mesh/mesh.cpp

@@ -75,7 +75,7 @@ void Mesh::GetElementCenter(int i, Vector &center)
 
 double Mesh::GetElementSize(ElementTransformation *T, int type)
 {
-   DenseMatrix J(Dim);
+   DenseMatrix J(spaceDim,Dim);
 
    Geometry::Type geom = T->GetGeometryType();
    T->SetIntPoint(&Geometries.GetCenter(geom));
@@ -83,7 +83,7 @@ double Mesh::GetElementSize(ElementTransformation *T, int type)
 
    if (type == 0)
    {
-      return pow(fabs(J.Det()), 1./Dim);
+      return pow(fabs(J.Weight()), 1./Dim);
    }
    else if (type == 1)
    {
@@ -102,7 +102,7 @@ double Mesh::GetElementSize(int i, int type)
 
 double Mesh::GetElementSize(int i, const Vector &dir)
 {
-   DenseMatrix J(Dim);
+   DenseMatrix J(spaceDim,Dim);
    Vector d_hat(Dim);
    GetElementJacobian(i, J);
    J.MultTranspose(dir, d_hat);

mesh/mesh.hpp

@@ -503,7 +503,7 @@ public:
    Mesh& operator=(Mesh &&mesh);
 
    /// Explicitly delete the copy assignment operator.
-   Mesh& operator=(Mesh &mesh) = delete;
+   Mesh& operator=(const Mesh &mesh) = delete;
 
    /** @name Named mesh constructors.
 

mesh/ncmesh.cpp

@@ -2283,7 +2283,6 @@ void NCMesh::GetMeshComponents(Mesh &mesh) const
    // left uninitialized here; they will be initialized later by the Mesh from
    // Nodes -- here we just make sure mesh.vertices has the correct size.
 
-   mesh.elements.SetSize(NElements);
    mesh.elements.SetSize(0);
 
    mesh.boundary.SetSize(0);

mesh/ncmesh.hpp

@@ -131,6 +131,7 @@ public:
    int GetNVertices() const { return NVertices; }
    int GetNEdges() const { return NEdges; }
    int GetNFaces() const { return NFaces; }
+   virtual int GetNGhostElements() const { return 0; }
 
    /** Perform the given batch of refinements. Please note that in the presence
        of anisotropic splits additional refinements may be necessary to keep
@@ -160,7 +161,6 @@ public:
        derefinements may have to be skipped to preserve mesh consistency. */
    virtual void Derefine(const Array<int> &derefs);
 
-
    // master/slave lists
 
    /// Identifies a vertex/edge/face in both Mesh and NCMesh.

mesh/pncmesh.cpp

@@ -1339,6 +1339,37 @@ void ParNCMesh::LimitNCLevel(int max_nc_level)
    }
 }
 
+void ParNCMesh::GetFineToCoarsePartitioning(const Array<int> &derefs,
+                                            Array<int> &new_ranks) const
+{
+   new_ranks.SetSize(leaf_elements.Size()-GetNGhostElements());
+   for (int i = 0; i < leaf_elements.Size()-GetNGhostElements(); i++)
+   {
+      new_ranks[i] = elements[leaf_elements[i]].rank;
+   }
+
+   for (int i = 0; i < derefs.Size(); i++)
+   {
+      int row = derefs[i];
+      MFEM_VERIFY(row >= 0 && row < derefinements.Size(),
+                  "invalid derefinement number.");
+
+      const int* fine = derefinements.GetRow(row);
+      int size = derefinements.RowSize(row);
+
+      int coarse_rank = INT_MAX;
+      for (int j = 0; j < size; j++)
+      {
+         int fine_rank = elements[leaf_elements[fine[j]]].rank;
+         coarse_rank = std::min(coarse_rank, fine_rank);
+      }
+      for (int j = 0; j < size; j++)
+      {
+         new_ranks[fine[j]] = coarse_rank;
+      }
+   }
+}
+
 void ParNCMesh::Derefine(const Array<int> &derefs)
 {
    MFEM_VERIFY(Dim < 3 || Iso,

mesh/pncmesh.hpp

@@ -95,6 +95,11 @@ public:
        in sync. The interface is identical. */
    virtual void Derefine(const Array<int> &derefs);
 
+   /** Gets partitioning for the coarse mesh if the current fine mesh were to
+       be derefined. */
+   virtual void GetFineToCoarsePartitioning(const Array<int> &derefs,
+                                            Array<int> &new_ranks) const;
+
    /** Migrate leaf elements of the global refinement hierarchy (including ghost
        elements) so that each processor owns the same number of leaves (+-1).
        The default partitioning strategy is based on equal splitting of the
@@ -243,7 +248,6 @@ public:
        The debug mesh will have element attributes set to element rank + 1. */
    void GetDebugMesh(Mesh &debug_mesh) const;
 
-
 protected: // interface for ParMesh
 
    friend class ParMesh;

miniapps/gslib/pfindpts.cpp

@@ -35,7 +35,7 @@
 //    mpirun -np 2 pfindpts -m ../../data/inline-wedge.mesh -o 3
 //    mpirun -np 2 pfindpts -m ../../data/amr-quad.mesh -o 2
 //    mpirun -np 2 pfindpts -m ../../data/rt-2d-q3.mesh -o 3 -mo 4 -ft 2
-
+//    mpirun -np 2 pfindpts -m ../../data/inline-quad.mesh -ft 1 -no-vis -sr0
 
 #include "mfem.hpp"
 
@@ -74,6 +74,7 @@ int main (int argc, char *argv[])
    bool visualization    = true;
    int fieldtype         = 0;
    int ncomp             = 1;
+   bool search_on_rank_0 = false;
 
    // Parse command-line options.
    OptionsParser args(argc, argv);
@@ -94,6 +95,9 @@ int main (int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&search_on_rank_0, "-sr0", "--search-on-r0", "-no-sr0",
+                  "--no-search-on-r0",
+                  "Enable search only on rank 0 (disable to search points on all tasks).");
    args.Parse();
    if (!args.Good())
    {
@@ -211,7 +215,7 @@ int main (int argc, char *argv[])
    // Note that some points might be outside, if the mesh is not a box. Note
    // also that all tasks search the same points (not mandatory).
    const int pts_cnt_1D = 10;
-   const int pts_cnt = pow(pts_cnt_1D, dim);
+   int pts_cnt = pow(pts_cnt_1D, dim);
    Vector vxyz(pts_cnt * dim);
    if (dim == 2)
    {
@@ -237,6 +241,12 @@ int main (int argc, char *argv[])
       }
    }
 
+   if ( (myid != 0) && (search_on_rank_0) )
+   {
+      pts_cnt = 0;
+      vxyz.Destroy();
+   }
+
    // Find and Interpolate FE function values on the desired points.
    Vector interp_vals(pts_cnt*vec_dim);
    FindPointsGSLIB finder(MPI_COMM_WORLD);
@@ -246,36 +256,37 @@ int main (int argc, char *argv[])
    Array<unsigned int> task_id_out = finder.GetProc();
    Vector dist_p_out = finder.GetDist();
 
-   int face_pts = 0, not_found = 0, found_loc = 0, found_away = 0;
-   double max_err = 0.0, max_dist = 0.0;
-   Vector pos(dim);
-   int npt = 0;
-   for (int j = 0; j < vec_dim; j++)
+   // Print the results for task 0 since either 1) all tasks have the
+   // same set of points or 2) only task 0 has any points.
+   if (myid == 0 )
    {
-      for (int i = 0; i < pts_cnt; i++)
+      int face_pts = 0, not_found = 0, found_loc = 0, found_away = 0;
+      double max_err = 0.0, max_dist = 0.0;
+      Vector pos(dim);
+      int npt = 0;
+      for (int j = 0; j < vec_dim; j++)
       {
-         if (j == 0)
+         for (int i = 0; i < pts_cnt; i++)
          {
-            (task_id_out[i] == (unsigned)myid) ? found_loc++ : found_away++;
-         }
+            if (j == 0)
+            {
+               (task_id_out[i] == (unsigned)myid) ? found_loc++ : found_away++;
+            }
 
-         if (code_out[i] < 2)
-         {
-            for (int d = 0; d < dim; d++) { pos(d) = vxyz(d * pts_cnt + i); }
-            Vector exact_val(vec_dim);
-            F_exact(pos, exact_val);
-            max_err  = std::max(max_err, fabs(exact_val(j) - interp_vals(npt)));
-            max_dist = std::max(max_dist, dist_p_out(i));
-            if (code_out[i] == 1 && j == 0) { face_pts++; }
+            if (code_out[i] < 2)
+            {
+               for (int d = 0; d < dim; d++) { pos(d) = vxyz(d * pts_cnt + i); }
+               Vector exact_val(vec_dim);
+               F_exact(pos, exact_val);
+               max_err  = std::max(max_err, fabs(exact_val(j) - interp_vals(npt)));
+               max_dist = std::max(max_dist, dist_p_out(i));
+               if (code_out[i] == 1 && j == 0) { face_pts++; }
+            }
+            else { if (j == 0) { not_found++; } }
+            npt++;
          }
-         else { if (j == 0) { not_found++; } }
-         npt++;
       }
-   }
 
-   // Print the results for task 0 since all tasks have the same set of points.
-   if (myid == 0)
-   {
       cout << setprecision(16)
            << "Searched unique points: " << pts_cnt
            << "\nFound on local mesh:  " << found_loc

miniapps/meshing/mesh-optimizer.cpp

@@ -39,6 +39,11 @@
 //   Adapted analytic shape+orientation:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
 //
+//   Adapted analytic shape and/or size with hr-adaptivity:
+//     mesh-optimizer -m square01.mesh -o 2 -tid 9  -ni 50 -li 20 -hmid 55 -mid 7 -hr
+//     mesh-optimizer -m square01.mesh -o 2 -tid 10 -ni 50 -li 20 -hmid 55 -mid 7 -hr
+//     mesh-optimizer -m square01.mesh -o 2 -tid 11 -ni 50 -li 20 -hmid 58 -mid 7 -hr
+//
 //   Adapted discrete size:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor
 //   Adapted discrete size 3D with PA:
@@ -123,6 +128,8 @@ int main(int argc, char *argv[])
    int max_lin_iter      = 100;
    bool move_bnd         = true;
    int combomet          = 0;
+   bool hradaptivity     = false;
+   int h_metric_id       = -1;
    bool normalization    = false;
    bool visualization    = true;
    int verbosity_level   = 0;
@@ -131,6 +138,8 @@ int main(int argc, char *argv[])
    bool exactaction      = false;
    const char *devopt    = "cpu";
    bool pa               = false;
+   int n_hr_iter         = 5;
+   int n_h_iter          = 1;
 
    // 1. Parse command-line options.
    OptionsParser args(argc, argv);
@@ -221,6 +230,12 @@ int main(int argc, char *argv[])
                   "0: Use single metric\n\t"
                   "1: Shape + space-dependent size given analytically\n\t"
                   "2: Shape + adapted size given discretely; shared target");
+   args.AddOption(&hradaptivity, "-hr", "--hr-adaptivity", "-no-hr",
+                  "--no-hr-adaptivity",
+                  "Enable hr-adaptivity.");
+   args.AddOption(&h_metric_id, "-hmid", "--h-metric",
+                  "Same options as metric_id. Used to determine refinement"
+                  " type for each element if h-adaptivity is enabled.");
    args.AddOption(&normalization, "-nor", "--normalization", "-no-nor",
                   "--no-normalization",
                   "Make all terms in the optimization functional unitless.");
@@ -241,6 +256,11 @@ int main(int argc, char *argv[])
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&n_hr_iter, "-nhr", "--n_hr_iter",
+                  "Number of hr-adaptivity iterations.");
+   args.AddOption(&n_h_iter, "-nh", "--n_h_iter",
+                  "Number of h-adaptivity iterations per r-adaptivity"
+                  "iteration.");
    args.Parse();
    if (!args.Good())
    {
@@ -249,6 +269,13 @@ int main(int argc, char *argv[])
    }
    args.PrintOptions(cout);
 
+   if (h_metric_id < 0) { h_metric_id = metric_id; }
+
+   if (hradaptivity)
+   {
+      MFEM_VERIFY(strcmp(devopt,"cpu")==0, "HR-adaptivity is currently only"
+                  " supported on cpus.");
+   }
    Device device(devopt);
    device.Print();
 
@@ -261,7 +288,9 @@ int main(int argc, char *argv[])
    else { cout << "(NONE)"; }
    cout << endl;
 
-   // 3. Define a finite element space on the mesh. Here we use vector finite
+   if (hradaptivity) { mesh->EnsureNCMesh(); }
+
+   // 3. Define a finite element space on the mesh-> Here we use vector finite
    //    elements which are tensor products of quadratic finite elements. The
    //    number of components in the vector finite element space is specified by
    //    the last parameter of the FiniteElementSpace constructor.
@@ -392,9 +421,41 @@ int main(int argc, char *argv[])
          cout << "Unknown metric_id: " << metric_id << endl;
          return 3;
    }
+   TMOP_QualityMetric *h_metric = NULL;
+   if (hradaptivity)
+   {
+      switch (h_metric_id)
+      {
+         case 1: h_metric = new TMOP_Metric_001; break;
+         case 2: h_metric = new TMOP_Metric_002; break;
+         case 7: h_metric = new TMOP_Metric_007; break;
+         case 9: h_metric = new TMOP_Metric_009; break;
+         case 55: h_metric = new TMOP_Metric_055; break;
+         case 56: h_metric = new TMOP_Metric_056; break;
+         case 58: h_metric = new TMOP_Metric_058; break;
+         case 77: h_metric = new TMOP_Metric_077; break;
+         case 315: h_metric = new TMOP_Metric_315; break;
+         case 316: h_metric = new TMOP_Metric_316; break;
+         case 321: h_metric = new TMOP_Metric_321; break;
+         default: cout << "Metric_id not supported for h-adaptivity: " << h_metric_id <<
+                          endl;
+            return 3;
+      }
+   }
+
+   if (metric_id < 300 || h_metric_id < 300)
+   {
+      MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
+   }
+   if (metric_id >= 300 || h_metric_id >= 300)
+   {
+      MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
+   }
+
    TargetConstructor::TargetType target_t;
    TargetConstructor *target_c = NULL;
    HessianCoefficient *adapt_coeff = NULL;
+   HRHessianCoefficient *hr_adapt_coeff = NULL;
    H1_FECollection ind_fec(mesh_poly_deg, dim);
    FiniteElementSpace ind_fes(mesh, &ind_fec);
    FiniteElementSpace ind_fesv(mesh, &ind_fec, dim);
@@ -604,6 +665,18 @@ int main(int argc, char *argv[])
          target_c = tc;
          break;
       }
+      // Targets used for hr-adaptivity tests.
+      case 9:  // size target in an annular region.
+      case 10: // size+aspect-ratio in an annular region.
+      case 11: // size+aspect-ratio target for a rotate sine wave
+      {
+         target_t = TargetConstructor::GIVEN_FULL;
+         AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
+         hr_adapt_coeff = new HRHessianCoefficient(dim, target_id - 9);
+         tc->SetAnalyticTargetSpec(NULL, NULL, hr_adapt_coeff);
+         target_c = tc;
+         break;
+      }
       default: cout << "Unknown target_id: " << target_id << endl; return 3;
    }
    if (target_c == NULL)
@@ -611,7 +684,8 @@ int main(int argc, char *argv[])
       target_c = new TargetConstructor(target_t);
    }
    target_c->SetNodes(x0);
-   TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c);
+   TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c,
+                                                       h_metric);
 
    // Finite differences for computations of derivatives.
    if (fdscheme)
@@ -721,10 +795,10 @@ int main(int argc, char *argv[])
             TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE);
          target_c2->SetVolumeScale(0.01);
          target_c2->SetNodes(x0);
-         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2);
+         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
          he_nlf_integ2->SetCoefficient(coeff2);
       }
-      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c); }
+      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
       he_nlf_integ2->SetIntegrationRules(*irules, quad_order);
       if (fdscheme) { he_nlf_integ2->EnableFiniteDifferences(x); }
       he_nlf_integ2->SetExactActionFlag(exactaction);
@@ -778,7 +852,9 @@ int main(int argc, char *argv[])
       tauval -= 0.01 * h0.Min();
    }
 
-   const double init_energy = a.GetGridFunctionEnergy(x);
+   // For HR tests, the energy is normalized by the number of elements.
+   const double init_energy = a.GetGridFunctionEnergy(x) /
+                              (hradaptivity ? mesh->GetNE() : 1);
 
    // Visualize the starting mesh and metric values.
    // Note that for combinations of metrics, this only shows the first metric.
@@ -906,13 +982,25 @@ int main(int argc, char *argv[])
       solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
    }
    solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
-   solver.SetOperator(a);
-   solver.Mult(b, x.GetTrueVector());
-   x.SetFromTrueVector();
-   if (solver.GetConverged() == false)
+
+   // hr-adaptivity solver.
+   // If hr-adaptivity is disabled, r-adaptivity is done once using the
+   // TMOPNewtonSolver.
+   // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
+   // "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity.
+   // The solver terminates if an h-adaptivity iteration does not modify
+   // any element in the mesh.
+   TMOPHRSolver hr_solver(*mesh, a, solver,
+                          x, move_bnd, hradaptivity,
+                          mesh_poly_deg, h_metric_id,
+                          n_hr_iter, n_h_iter);
+   hr_solver.AddGridFunctionForUpdate(&x0);
+   if (adapt_lim_const > 0.)
    {
-      cout << "Nonlinear solver: rtol = " << solver_rtol << " not achieved.\n";
+      hr_solver.AddGridFunctionForUpdate(&zeta_0);
+      hr_solver.AddFESpaceForUpdate(&ind_fes);
    }
+   hr_solver.Mult();
 
    // 15. Save the optimized mesh to a file. This output can be viewed later
    //     using GLVis: "glvis -m optimized.mesh".
@@ -922,14 +1010,15 @@ int main(int argc, char *argv[])
       mesh->Print(mesh_ofs);
    }
 
-   // 16. Compute the amount of energy decrease.
-   const double fin_energy = a.GetGridFunctionEnergy(x);
+   const double fin_energy = a.GetGridFunctionEnergy(x) /
+                             (hradaptivity ? mesh->GetNE() : 1);
    double metric_part = fin_energy;
    if (lim_const > 0.0 || adapt_lim_const > 0.0)
    {
       lim_coeff.constant = 0.0;
       coef_zeta.constant = 0.0;
-      metric_part = a.GetGridFunctionEnergy(x);
+      metric_part = a.GetGridFunctionEnergy(x) /
+                    (hradaptivity ? mesh->GetNE() : 1);
       lim_coeff.constant = lim_const;
       coef_zeta.constant = adapt_lim_const;
    }
@@ -942,7 +1031,7 @@ int main(int argc, char *argv[])
    cout << "The strain energy decreased by: " << setprecision(12)
         << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
 
-   // 17. Visualize the final mesh and metric values.
+   // 16. Visualize the final mesh and metric values.
    if (visualization)
    {
       char title[] = "Final metric values";
@@ -956,13 +1045,13 @@ int main(int argc, char *argv[])
                              600, 600, 300, 300);
    }
 
-   // 18. Visualize the mesh displacement.
+   // 17. Visualize the mesh displacement.
    if (visualization)
    {
-      x0 -= x;
       osockstream sock(19916, "localhost");
       sock << "solution\n";
       mesh->Print(sock);
+      x0 -= x;
       x0.Save(sock);
       sock.send();
       sock << "window_title 'Displacements'\n"
@@ -971,7 +1060,6 @@ int main(int argc, char *argv[])
            << "keys jRmclA" << endl;
    }
 
-   // 19. Free the used memory.
    delete S;
    delete S_prec;
    delete target_c2;
@@ -979,7 +1067,9 @@ int main(int argc, char *argv[])
    delete coeff1;
    delete adapt_evaluator;
    delete target_c;
+   delete hr_adapt_coeff;
    delete adapt_coeff;
+   delete h_metric;
    delete metric;
    delete fespace;
    delete fec;

miniapps/meshing/mesh-optimizer.hpp

@@ -90,8 +90,6 @@ double discrete_aspr_2d(const Vector &x)
    double th = 22.5*M_PI/180.;
    double xn =  cos(th)*xc + sin(th)*yc;
    double yn = -sin(th)*xc + cos(th)*yc;
-   // double th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
-   // double stretch = 1/cos(th2);
    xc = xn; yc = yn;
 
    double tfac = 20;
@@ -223,6 +221,134 @@ public:
    }
 };
 
+class HRHessianCoefficient : public TMOPMatrixCoefficient
+{
+private:
+   int dim;
+   // 0 - size target in an annular region,
+   // 1 - size+aspect-ratio in an annular region,
+   // 2 - size+aspect-ratio target for a rotate sine wave.
+   int hr_target_type;
+
+public:
+   HRHessianCoefficient(int dim_, int hr_target_type_ = 0)
+      : TMOPMatrixCoefficient(dim_), dim(dim_),
+        hr_target_type(hr_target_type_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      Vector pos(3);
+      T.Transform(ip, pos);
+      if (hr_target_type == 0) // size only circle
+      {
+         double small = 0.001, big = 0.01;
+         if (dim == 3) { small = 0.005, big = 0.1; }
+         const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
+         double zc;
+         if (dim == 3) { zc = pos(2) - 0.5; }
+         double r = sqrt(xc*xc + yc*yc);
+         if (dim == 3) { r = sqrt(xc*xc + yc*yc + zc*zc); }
+         double r1 = 0.15; double r2 = 0.35; double sf=30.0;
+
+         const double tan1 = std::tanh(sf*(r-r1)),
+                      tan2 = std::tanh(sf*(r-r2));
+
+         double ind = (tan1 - tan2);
+         if (ind > 1.0) {ind = 1.;}
+         if (ind < 0.0) {ind = 0.;}
+         double val = ind * small + (1.0 - ind) * big;
+         K = 0.0;
+         K(0, 0) = 1.0;
+         K(0, 1) = 0.0;
+         K(1, 0) = 0.0;
+         K(1, 1) = 1.0;
+         K(0, 0) *= pow(val,0.5);
+         K(1, 1) *= pow(val,0.5);
+         if (dim == 3) { K(2, 2) = pow(val,0.5); }
+      }
+      else if (hr_target_type == 1) // circle with size and AR
+      {
+         const double small = 0.001, big = 0.01;
+         const double xc = pos(0)-0.5, yc = pos(1)-0.5;
+         const double rv = xc*xc + yc*yc;
+         double r = 0;
+         if (rv>0.) {r = sqrt(rv);}
+
+         double r1 = 0.2; double r2 = 0.3; double sf=30.0;
+         const double szfac = 1;
+         const double asfac = 4;
+         const double eps2 = szfac/asfac;
+         const double eps1 = szfac;
+
+         double tan1 = std::tanh(sf*(r-r1)+1),
+                tan2 = std::tanh(sf*(r-r2)-1);
+         double wgt = 0.5*(tan1-tan2);
+
+         tan1 = std::tanh(sf*(r-r1)),
+         tan2 = std::tanh(sf*(r-r2));
+
+         double ind = (tan1 - tan2);
+         if (ind > 1.0) {ind = 1.;}
+         if (ind < 0.0) {ind = 0.;}
+         double szval = ind * small + (1.0 - ind) * big;
+
+         double th = std::atan2(yc,xc)*180./M_PI;
+         if (wgt > 1) { wgt = 1; }
+         if (wgt < 0) { wgt = 0; }
+
+         double maxval = eps2 + eps1*(1-wgt)*(1-wgt);
+         double minval = eps1;
+         double avgval = 0.5*(maxval+minval);
+         double ampval = 0.5*(maxval-minval);
+         double val1 = avgval + ampval*sin(2.*th*M_PI/180.+90*M_PI/180.);
+         double val2 = avgval + ampval*sin(2.*th*M_PI/180.-90*M_PI/180.);
+
+         K(0,1) = 0.0;
+         K(1,0) = 0.0;
+         K(0,0) = val1;
+         K(1,1) = val2;
+
+         K(0,0) *= pow(szval,0.5);
+         K(1,1) *= pow(szval,0.5);
+      }
+      else if (hr_target_type == 2) // sharp rotated sine wave
+      {
+         double xc = pos(0)-0.5, yc = pos(1)-0.5;
+         double th = 15.5*M_PI/180.;
+         double xn =  cos(th)*xc + sin(th)*yc;
+         double yn = -sin(th)*xc + cos(th)*yc;
+         double th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
+         double stretch = 1/cos(th2);
+         xc = xn/stretch;
+         yc = yn;
+         double tfac = 20;
+         double s1 = 3;
+         double s2 = 2;
+         double yl1 = -0.025;
+         double yl2 =  0.025;
+         double wgt = std::tanh((tfac*(yc-yl1) + s2*std::sin(s1*M_PI*xc)) + 1) -
+                      std::tanh((tfac*(yc-yl2) + s2*std::sin(s1*M_PI*xc)) - 1);
+         if (wgt > 1) { wgt = 1; }
+         if (wgt < 0) { wgt = 0; }
+
+         const double eps2 = 25;
+         const double eps1 = 1;
+         K(1,1) = eps1/eps2 + eps1*(1-wgt)*(1-wgt);
+         K(0,0) = eps1;
+         K(0,1) = 0.0;
+         K(1,0) = 0.0;
+      }
+      else { MFEM_ABORT("Unsupported option / wrong input."); }
+   }
+
+   virtual void EvalGrad(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip, int comp)
+   {
+      K = 0.;
+   }
+};
+
 // Additional IntegrationRules that can be used with the --quad-type option.
 IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto);
 IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform);
@@ -252,13 +378,13 @@ double adapt_lim_fun(const Vector &x)
 
 void DiffuseField(GridFunction &field, int smooth_steps)
 {
-   //Setup the Laplacian operator
+   // Setup the Laplacian operator
    BilinearForm *Lap = new BilinearForm(field.FESpace());
    Lap->AddDomainIntegrator(new DiffusionIntegrator());
    Lap->Assemble();
    Lap->Finalize();
 
-   //Setup the smoothing operator
+   // Setup the smoothing operator
    DSmoother *S = new DSmoother(0,1.0,smooth_steps);
    S->iterative_mode = true;
    S->SetOperator(Lap->SpMat());
@@ -274,7 +400,7 @@ void DiffuseField(GridFunction &field, int smooth_steps)
 #ifdef MFEM_USE_MPI
 void DiffuseField(ParGridFunction &field, int smooth_steps)
 {
-   //Setup the Laplacian operator
+   // Setup the Laplacian operator
    ParBilinearForm *Lap = new ParBilinearForm(field.ParFESpace());
    Lap->AddDomainIntegrator(new DiffusionIntegrator());
    Lap->Assemble();

miniapps/meshing/pmesh-optimizer.cpp

@@ -39,6 +39,11 @@
 //   Adapted analytic shape+orientation:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
 //
+//   Adapted analytic shape and/or size with hr-adaptivity:
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 9  -ni 50 -li 20 -hmid 55 -mid 7 -hr
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 10 -ni 50 -li 20 -hmid 55 -mid 7 -hr
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 11 -ni 50 -li 20 -hmid 58 -mid 7 -hr
+//
 //   Adapted discrete size:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor
 //   Adapted discrete size 3D with PA:
@@ -130,6 +135,8 @@ int main (int argc, char *argv[])
    int max_lin_iter      = 100;
    bool move_bnd         = true;
    int combomet          = 0;
+   bool hradaptivity     = false;
+   int h_metric_id       = -1;
    bool normalization    = false;
    bool visualization    = true;
    int verbosity_level   = 0;
@@ -138,6 +145,8 @@ int main (int argc, char *argv[])
    bool exactaction      = false;
    const char *devopt    = "cpu";
    bool pa               = false;
+   int n_hr_iter         = 5;
+   int n_h_iter          = 1;
 
    // 2. Parse command-line options.
    OptionsParser args(argc, argv);
@@ -230,6 +239,12 @@ int main (int argc, char *argv[])
                   "0: Use single metric\n\t"
                   "1: Shape + space-dependent size given analytically\n\t"
                   "2: Shape + adapted size given discretely; shared target");
+   args.AddOption(&hradaptivity, "-hr", "--hr-adaptivity", "-no-hr",
+                  "--no-hr-adaptivity",
+                  "Enable hr-adaptivity.");
+   args.AddOption(&h_metric_id, "-hmid", "--h-metric",
+                  "Same options as metric_id. Used to determine refinement"
+                  " type for each element if h-adaptivity is enabled.");
    args.AddOption(&normalization, "-nor", "--normalization", "-no-nor",
                   "--no-normalization",
                   "Make all terms in the optimization functional unitless.");
@@ -250,6 +265,11 @@ int main (int argc, char *argv[])
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&n_hr_iter, "-nhr", "--n_hr_iter",
+                  "Number of hr-adaptivity iterations.");
+   args.AddOption(&n_h_iter, "-nh", "--n_h_iter",
+                  "Number of h-adaptivity iterations per r-adaptivity"
+                  "iteration.");
    args.Parse();
    if (!args.Good())
    {
@@ -257,7 +277,13 @@ int main (int argc, char *argv[])
       return 1;
    }
    if (myid == 0) { args.PrintOptions(cout); }
+   if (h_metric_id < 0) { h_metric_id = metric_id; }
 
+   if (hradaptivity)
+   {
+      MFEM_VERIFY(strcmp(devopt,"cpu")==0, "HR-adaptivity is currently only"
+                  " supported on cpus.");
+   }
    Device device(devopt);
    if (myid == 0) { device.Print();}
 
@@ -276,6 +302,7 @@ int main (int argc, char *argv[])
       cout << endl;
    }
 
+   if (hradaptivity) { mesh->EnsureNCMesh(); }
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
 
    delete mesh;
@@ -422,9 +449,41 @@ int main (int argc, char *argv[])
          if (myid == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
          return 3;
    }
+   TMOP_QualityMetric *h_metric = NULL;
+   if (hradaptivity)
+   {
+      switch (h_metric_id)
+      {
+         case 1: h_metric = new TMOP_Metric_001; break;
+         case 2: h_metric = new TMOP_Metric_002; break;
+         case 7: h_metric = new TMOP_Metric_007; break;
+         case 9: h_metric = new TMOP_Metric_009; break;
+         case 55: h_metric = new TMOP_Metric_055; break;
+         case 56: h_metric = new TMOP_Metric_056; break;
+         case 58: h_metric = new TMOP_Metric_058; break;
+         case 77: h_metric = new TMOP_Metric_077; break;
+         case 315: h_metric = new TMOP_Metric_315; break;
+         case 316: h_metric = new TMOP_Metric_316; break;
+         case 321: h_metric = new TMOP_Metric_321; break;
+         default: cout << "Metric_id not supported for h-adaptivity: " << h_metric_id <<
+                          endl;
+            return 3;
+      }
+   }
+
+   if (metric_id < 300 || h_metric_id < 300)
+   {
+      MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
+   }
+   if (metric_id >= 300 || h_metric_id >= 300)
+   {
+      MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
+   }
+
    TargetConstructor::TargetType target_t;
    TargetConstructor *target_c = NULL;
    HessianCoefficient *adapt_coeff = NULL;
+   HRHessianCoefficient *hr_adapt_coeff = NULL;
    H1_FECollection ind_fec(mesh_poly_deg, dim);
    ParFiniteElementSpace ind_fes(pmesh, &ind_fec);
    ParFiniteElementSpace ind_fesv(pmesh, &ind_fec, dim);
@@ -639,6 +698,18 @@ int main (int argc, char *argv[])
          target_c = tc;
          break;
       }
+      // Targets used for hr-adaptivity tests.
+      case 9:  // size target in an annular region.
+      case 10: // size+aspect-ratio in an annular region.
+      case 11: // size+aspect-ratio target for a rotate sine wave
+      {
+         target_t = TargetConstructor::GIVEN_FULL;
+         AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
+         hr_adapt_coeff = new HRHessianCoefficient(dim, target_id - 9);
+         tc->SetAnalyticTargetSpec(NULL, NULL, hr_adapt_coeff);
+         target_c = tc;
+         break;
+      }
       default:
          if (myid == 0) { cout << "Unknown target_id: " << target_id << endl; }
          return 3;
@@ -649,13 +720,13 @@ int main (int argc, char *argv[])
       target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
    }
    target_c->SetNodes(x0);
-   TMOP_Integrator *he_nlf_integ= new TMOP_Integrator(metric, target_c);
+   TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c,
+                                                       h_metric);
 
    // Finite differences for computations of derivatives.
    if (fdscheme)
    {
       MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
-
       he_nlf_integ->EnableFiniteDifferences(x);
    }
    he_nlf_integ->SetExactActionFlag(exactaction);
@@ -753,7 +824,8 @@ int main (int argc, char *argv[])
       he_nlf_integ->SetCoefficient(*coeff1);
 
       // Second metric.
-      metric2 = new TMOP_Metric_077;
+      if (dim == 2) { metric2 = new TMOP_Metric_077; }
+      else          { metric2 = new TMOP_Metric_315; }
       TMOP_Integrator *he_nlf_integ2 = NULL;
       if (combomet == 1)
       {
@@ -761,10 +833,10 @@ int main (int argc, char *argv[])
             TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE, MPI_COMM_WORLD);
          target_c2->SetVolumeScale(0.01);
          target_c2->SetNodes(x0);
-         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2);
+         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
          he_nlf_integ2->SetCoefficient(coeff2);
       }
-      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c); }
+      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
       he_nlf_integ2->SetIntegrationRules(*irules, quad_order);
       if (fdscheme) { he_nlf_integ2->EnableFiniteDifferences(x); }
       he_nlf_integ2->SetExactActionFlag(exactaction);
@@ -777,7 +849,10 @@ int main (int argc, char *argv[])
 
       a.AddDomainIntegrator(combo);
    }
-   else { a.AddDomainIntegrator(he_nlf_integ); }
+   else
+   {
+      a.AddDomainIntegrator(he_nlf_integ);
+   }
 
    if (pa) { a.Setup(); }
 
@@ -821,7 +896,9 @@ int main (int argc, char *argv[])
       tauval -= 0.01 * h0min_all;
    }
 
-   const double init_energy = a.GetParGridFunctionEnergy(x);
+   // For HR tests, the energy is normalized by the number of elements.
+   const double init_energy = a.GetParGridFunctionEnergy(x) /
+                              (hradaptivity ? pmesh->GetGlobalNE() : 1);
 
    // Visualize the starting mesh and metric values.
    // Note that for combinations of metrics, this only shows the first metric.
@@ -952,13 +1029,25 @@ int main (int argc, char *argv[])
       solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
    }
    solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
-   solver.SetOperator(a);
-   solver.Mult(b, x.GetTrueVector());
-   x.SetFromTrueVector();
-   if (myid == 0 && solver.GetConverged() == false)
+
+   // hr-adaptivity solver.
+   // If hr-adaptivity is disabled, r-adaptivity is done once using the
+   // TMOPNewtonSolver.
+   // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
+   // "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity.
+   // The solver terminates if an h-adaptivity iteration does not modify
+   // any element in the mesh.
+   TMOPHRSolver hr_solver(*pmesh, a, solver,
+                          x, move_bnd, hradaptivity,
+                          mesh_poly_deg, h_metric_id,
+                          n_hr_iter, n_h_iter);
+   hr_solver.AddGridFunctionForUpdate(&x0);
+   if (adapt_lim_const > 0.)
    {
-      cout << "Nonlinear solver: rtol = " << solver_rtol << " not achieved.\n";
+      hr_solver.AddGridFunctionForUpdate(&zeta_0);
+      hr_solver.AddFESpaceForUpdate(&ind_fes);
    }
+   hr_solver.Mult();
 
    // 16. Save the optimized mesh to a file. This output can be viewed later
    //     using GLVis: "glvis -m optimized -np num_mpi_tasks".
@@ -971,13 +1060,15 @@ int main (int argc, char *argv[])
    }
 
    // 17. Compute the amount of energy decrease.
-   const double fin_energy = a.GetParGridFunctionEnergy(x);
+   const double fin_energy = a.GetParGridFunctionEnergy(x) /
+                             (hradaptivity ? pmesh->GetGlobalNE() : 1);
    double metric_part = fin_energy;
    if (lim_const > 0.0 || adapt_lim_const > 0.0)
    {
       lim_coeff.constant = 0.0;
       coef_zeta.constant = 0.0;
-      metric_part = a.GetParGridFunctionEnergy(x);
+      metric_part = a.GetParGridFunctionEnergy(x) /
+                    (hradaptivity ? pmesh->GetGlobalNE() : 1);
       lim_coeff.constant = lim_const;
       coef_zeta.constant = adapt_lim_const;
    }
@@ -1036,7 +1127,9 @@ int main (int argc, char *argv[])
    delete coeff1;
    delete adapt_evaluator;
    delete target_c;
+   delete hr_adapt_coeff;
    delete adapt_coeff;
+   delete h_metric;
    delete metric;
    delete pfespace;
    delete fec;

tests/gitlab/build_and_test

@@ -56,8 +56,10 @@ then
         prefix_opt="--prefix=${prefix}"
     fi
 
+    echo "Fetching uberenv."
     tests/gitlab/get_mfem_uberenv || ( echo "Error fetching Uberenv" && exit 1 );
 
+    echo "Installing dependencies."
     python tests/uberenv/uberenv.py --spec="${spec}" "${upstream_opt}" "${prefix_opt}"
 
 fi

tests/gitlab/get_mfem_uberenv

@@ -18,7 +18,7 @@ set -o errexit
 set -o nounset
 
 uberenv_url="https://github.com/mfem/mfem-uberenv.git"
-uberenv_ref="dbe35eb4dbeaabc7184ad110fcb98377092baeb0"
+uberenv_ref="222e3743377d8c903e044f80ae44624e1668feb4"
 
 [[ ! -d tests/uberenv ]] && git clone ${uberenv_url} tests/uberenv
 cd tests/uberenv