Merge branch 'master' into minsurface-bugfix

Add GPU support for PWConstCoefficient

CHANGELOG

@@ -27,6 +27,19 @@ Discretization improvements
 - In the ParMoonolith integration, added support for variational resampling of
   H1 vector fields.
 
+- Added support for boundary integration to the hyperbolic framework. In this
+  regard, new classes `BdrHyperbolicDirichletIntegrator` and
+  `BoundaryHyperbolicFlowIntegrator` have been introduced for implementation
+  of weak Dirichlet boundary conditions with a general flux or for the linear
+  case respectively.
+
+- Added method to compute piecewise linear bounds on high-order functions on
+  tensor-product elements.
+
+- Parallel anisotropic refinement of hexahedral meshes is now supported,
+  provided that neighboring hexahedra are not refined in conflicting directions.
+  A new ParMesh method is added to check for such conflicts, before refinement.
+
 Meshing improvements
 --------------------
 
@@ -44,8 +57,10 @@ GPU computing
   set. This is most often used for setting constant essential boundary
   conditions. A new function Vector::SetSubVectorHost has been added in cases
   where host execution is always needed (e.g. when the DOFs array is small).
+
 - Introduced MFEM_FOREACH_THREAD_DIRECT, which directly maps loop tasks to GPU
   threads, assigning one task per thread.
+
 - Implemented a GPU-accelerated matrix-free AMR derefinement `GridFunction`
   update operator. This supports mixed geometry meshes and variable order
   spaces, and is the default derefinement operator constructed by
@@ -60,10 +75,18 @@ New and updated examples and miniapps
   operators as smoothers.
   These miniapps can be found in `miniapps/diag-smoothers`.
 
-API changes
+- Added a new miniapp (meshing/mesh-bounding-boxes) that computes the bounding
+  boxes for each element of a given mesh, and the bounds on the determinant of
+  the Jacobian of the transformation.
+
+- Added a new miniapp (tools/gridfunction-bounds) to compute piecewise linear
+  bounds on a given high-order grid function.
+
+API changes:
 -----------
 - mfem::internal::tensor and mfem::internal::dual have been moved to
   mfem::future::tensor and mfem::future::dual.
+
 - API addition: in class `Operator`, added virtual functions: `AbsMult`, and
   `AbsMultTranspose`; in class `Vector`, added `Abs` and `Pow`.
 
@@ -71,16 +94,25 @@ Miscellaneous
 -------------
 - Added the "gpu", "raja-gpu", and "ceed-gpu" backend aliases/shortcuts which
   automatically select between CUDA or HIP.
+
 - The CUDA-specific names used by some of the unit tests like 'cunit_tests' and
   'pcunit_tests' were replaced by names using 'gpu' instead of 'c' (short for
   CUDA) or 'cuda'. These tests automatically run the CUDA/HIP tests based on the
   MFEM build configuration.
+
 - Added the option to enable GPU-aware MPI in MFEM using the environment
   variable 'MFEM_GPU_AWARE_MPI' set to any value. Setting this environment
   variable is an alternative to calling 'Device::SetGPUAwareMPI(true)'.
+
 - Added parallel Address Sanitizer, serial and parallel Undefined Behavior
   Sanitizer and serial Memory Sanitizer GitHub actions tests on Ubuntu.
 
+- FindPointsGSLIB has a new constructor that accepts the mesh object and
+  internally calls the Setup() method so that the user does not have to.
+  The FreeData() method has also been moved to the destructor so the user does
+  not need to manually free-up the memory if the destructor is called before
+  MPI_Finalize().
+
 Version 4.8, released on Apr 9, 2025
 ====================================
 

fem/bilininteg.hpp

@@ -23,6 +23,8 @@
 
 namespace mfem
 {
+class QuadratureSpace;
+class FaceQuadratureSpace;
 
 /// Abstract base class BilinearFormIntegrator
 class BilinearFormIntegrator : public NonlinearFormIntegrator

fem/coefficient.cpp

@@ -12,6 +12,7 @@
 // Implementation of Coefficient class
 
 #include "fem.hpp"
+#include "../general/forall.hpp"
 
 #include <cmath>
 #include <limits>
@@ -80,6 +81,49 @@ real_t PWConstCoefficient::Eval(ElementTransformation & T,
    return (constants(att-1));
 }
 
+void PWConstCoefficient::Project(QuadratureFunction &qf)
+{
+   auto &qs = *qf.GetSpace();
+
+   const bool compressed =
+      qs.Offsets(QSpaceOffsetStorage::COMPRESSED).Size() == 1;
+   const int *offsets = qs.Offsets(QSpaceOffsetStorage::COMPRESSED).Read();
+   const int ne = qs.GetNE();
+
+   const int *attributes = [&]()
+   {
+      if (dynamic_cast<QuadratureSpace*>(&qs) != nullptr)
+      {
+         return qs.GetMesh()->GetElementAttributes().Read();
+      }
+      else if (auto *qs_f = dynamic_cast<FaceQuadratureSpace*>(&qs))
+      {
+         MFEM_VERIFY(qs_f->GetFaceType() == FaceType::Boundary,
+                     "Interior faces do not have attributes.");
+         return qs.GetMesh()->GetBdrFaceAttributes().Read();
+      }
+      else
+      {
+         MFEM_ABORT("Unsupported case.");
+      }
+   }();
+
+   const real_t *d_c = constants.Read();
+   real_t *d_qf = qf.Write();
+
+   mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
+   {
+      const int a = attributes[e];
+      const real_t elementConstant = d_c[a - 1];
+      const int begin = compressed ? e*offsets[0] : offsets[e];
+      const int end = compressed ? (e+1)*offsets[0] : offsets[e+1];
+      for (int i = begin; i < end; ++i)
+      {
+         d_qf[i] = elementConstant;
+      }
+   });
+}
+
 void PWCoefficient::InitMap(const Array<int> & attr,
                             const Array<Coefficient*> & coefs)
 {

fem/coefficient.hpp

@@ -132,6 +132,9 @@ public:
    /// Evaluate the coefficient.
    real_t Eval(ElementTransformation &T,
                const IntegrationPoint &ip) override;
+
+   /// Fill the QuadratureFunction @a qf with the piecewise constant values.
+   void Project(QuadratureFunction &qf) override;
 };
 
 /** @brief A piecewise coefficient with the pieces keyed off the element

fem/dfem/doperator.hpp

@@ -231,22 +231,53 @@ public:
       const std::vector<FieldDescriptor> &parameters,
       const ParMesh &mesh);
 
+   /// MultLevel enum to indicate if the T->L Operators are used in the
+   /// Mult method.
+   enum MultLevel
+   {
+      TVECTOR,
+      LVECTOR
+   };
+
+   /// @brief Set the MultLevel mode for the DifferentiableOperator.
+   /// The default is TVECTOR, which means that the Operator will use
+   /// T->L before Mult and L->T Operators after.
+   void SetMultLevel(MultLevel level)
+   {
+      mult_level = level;
+   }
+
    /// @brief Compute the action of the operator on a given vector.
    ///
-   /// @param solutions_t The solution vector in which to compute the action.
-   /// This has to be a T-dof vector.
-   /// @param result_t Result vector of the action of the operator on
-   /// solutions_t. The result is a T-dof vector.
-   void Mult(const Vector &solutions_t, Vector &result_t) const override
+   /// @param solutions_in The solution vector in which to compute the action.
+   /// This has to be a T-dof vector if MultLevel is set to TVECTOR, or L-dof
+   /// Vector if MultLevel is set to LVECTOR.
+   /// @param result_in Result vector of the action of the operator on
+   /// solutions. The result is a T-dof vector or L-dof vector depending on
+   /// the MultLevel.
+   void Mult(const Vector &solutions_in, Vector &result_in) const override
    {
       MFEM_ASSERT(!action_callbacks.empty(), "no integrators have been set");
-      prolongation(solutions, solutions_t, solutions_l);
-      residual_l = 0.0;
-      for (auto &action : action_callbacks)
+
+      if (mult_level == MultLevel::LVECTOR)
       {
-         action(solutions_l, parameters_l, residual_l);
+         get_lvectors(solutions, solutions_in, solutions_l);
+         result_in = 0.0;
+         for (auto &action : action_callbacks)
+         {
+            action(solutions_l, parameters_l, result_in);
+         }
+      }
+      else
+      {
+         prolongation(solutions, solutions_in, solutions_l);
+         residual_l = 0.0;
+         for (auto &action : action_callbacks)
+         {
+            action(solutions_l, parameters_l, residual_l);
+         }
+         prolongation_transpose(residual_l, result_in);
       }
-      prolongation_transpose(residual_l, result_t);
    }
 
    /// @brief Add a domain integrator to the operator.
@@ -345,6 +376,8 @@ public:
 private:
    const ParMesh &mesh;
 
+   MultLevel mult_level = TVECTOR;
+
    std::vector<action_t> action_callbacks;
    std::map<size_t,
        std::vector<derivative_action_t>> derivative_action_callbacks;
@@ -354,7 +387,6 @@ private:
        std::vector<assemble_derivative_hypreparmatrix_callback_t>>
        assemble_derivative_hypreparmatrix_callbacks;
 
-
    std::vector<FieldDescriptor> solutions;
    std::vector<FieldDescriptor> parameters;
    // solutions and parameters

fem/dfem/util.hpp

@@ -1076,6 +1076,24 @@ void prolongation(const std::vector<FieldDescriptor> fields,
    }
 }
 
+inline
+void get_lvectors(const std::vector<FieldDescriptor> fields,
+                  const Vector &x,
+                  std::vector<Vector> &fields_l)
+{
+   int data_offset = 0;
+   for (std::size_t i = 0; i < fields.size(); i++)
+   {
+      const int sz = GetVSize(fields[i]);
+      fields_l[i].SetSize(sz);
+
+      const Vector x_i(const_cast<Vector&>(x), data_offset, sz);
+      fields_l[i] = x_i;
+
+      data_offset += sz;
+   }
+}
+
 /// @brief Get a transpose prolongation callback for a field descriptor.
 ///
 /// In the special case of a one field operator, the transpose prolongation

fem/fe_coll.cpp

@@ -2515,6 +2515,7 @@ RT_FECollection::RT_FECollection(const int p, const int dim,
                                  const int map_type, const bool signs,
                                  const int ob_type)
    : FiniteElementCollection(p + 1)
+   , dim(dim)
    , ob_type(ob_type)
 {
    if (Quadrature1D::CheckOpen(BasisType::GetQuadrature1D(ob_type)) ==

fem/fe_coll.hpp

@@ -464,6 +464,13 @@ public:
    RT_Trace_FECollection(const int p, const int dim,
                          const int map_type = FiniteElement::INTEGRAL,
                          const int ob_type = BasisType::GaussLegendre);
+
+   FiniteElementCollection *Clone(int p) const override
+   {
+      const int map_type = (strncmp(rt_name, "RT_Trace", 8) == 0)?
+                           (FiniteElement::INTEGRAL):(FiniteElement::VALUE);
+      return new RT_Trace_FECollection(p, dim, map_type, ob_type);
+   }
 };
 
 /** Arbitrary order discontinuous finite elements defined on the interface
@@ -475,6 +482,13 @@ public:
    DG_Interface_FECollection(const int p, const int dim,
                              const int map_type = FiniteElement::VALUE,
                              const int ob_type = BasisType::GaussLegendre);
+
+   FiniteElementCollection *Clone(int p) const override
+   {
+      const int map_type = (strncmp(rt_name, "DG_Iface", 8) == 0)?
+                           (FiniteElement::VALUE):(FiniteElement::INTEGRAL);
+      return new DG_Interface_FECollection(p, dim, map_type, ob_type);
+   }
 };
 
 /// Arbitrary order H(curl)-conforming Nedelec finite elements.

fem/gslib.cpp

@@ -85,9 +85,9 @@ namespace mfem
 {
 
 FindPointsGSLIB::FindPointsGSLIB()
-   : mesh(NULL),
-     fec_map_lin(NULL),
-     fdataD(NULL), cr(NULL), gsl_comm(NULL),
+   : mesh(nullptr),
+     fec_map_lin(nullptr),
+     fdataD(nullptr), cr(nullptr), gsl_comm(nullptr),
      dim(-1), points_cnt(-1), setupflag(false), default_interp_value(0),
      avgtype(AvgType::ARITHMETIC), bdr_tol(1e-8)
 {
@@ -97,10 +97,10 @@ FindPointsGSLIB::FindPointsGSLIB()
    gf_rst_map.SetSize(4);
    for (int i = 0; i < mesh_split.Size(); i++)
    {
-      mesh_split[i] = NULL;
-      ir_split[i] = NULL;
-      fes_rst_map[i] = NULL;
-      gf_rst_map[i] = NULL;
+      mesh_split[i] = nullptr;
+      ir_split[i] = nullptr;
+      fes_rst_map[i] = nullptr;
+      gf_rst_map[i] = nullptr;
    }
 
    gsl_comm = new gslib::comm;
@@ -117,27 +117,40 @@ FindPointsGSLIB::FindPointsGSLIB()
    crystal_init(cr, gsl_comm);
 }
 
+FindPointsGSLIB::FindPointsGSLIB(Mesh &mesh_in, const double bb_t,
+                                 const double newt_tol, const int npt_max)
+   : FindPointsGSLIB()
+{
+   Setup(mesh_in, bb_t, newt_tol, npt_max);
+}
+
 FindPointsGSLIB::~FindPointsGSLIB()
 {
-   crystal_free(cr);
-   comm_free(gsl_comm);
-   delete gsl_comm;
-   delete cr;
-   for (int i = 0; i < 4; i++)
+   FreeData();
+#ifdef MFEM_USE_MPI
+   if (!Mpi::IsFinalized())  // currently segfaults inside gslib otherwise
+#endif
    {
-      if (mesh_split[i]) { delete mesh_split[i]; mesh_split[i] = NULL; }
-      if (ir_split[i]) { delete ir_split[i]; ir_split[i] = NULL; }
-      if (fes_rst_map[i]) { delete fes_rst_map[i]; fes_rst_map[i] = NULL; }
-      if (gf_rst_map[i]) { delete gf_rst_map[i]; gf_rst_map[i] = NULL; }
+      crystal_free(cr);
+      comm_free(gsl_comm);
+      delete gsl_comm;
+      delete cr;
    }
-   if (fec_map_lin) { delete fec_map_lin; fec_map_lin = NULL; }
+   for (int i = 0; i < mesh_split.Size(); i++)
+   {
+      if (mesh_split[i]) { delete mesh_split[i]; mesh_split[i] = nullptr; }
+      if (ir_split[i]) { delete ir_split[i]; ir_split[i] = nullptr; }
+      if (fes_rst_map[i]) { delete fes_rst_map[i]; fes_rst_map[i] = nullptr; }
+      if (gf_rst_map[i]) { delete gf_rst_map[i]; gf_rst_map[i] = nullptr; }
+   }
+   if (fec_map_lin) { delete fec_map_lin; fec_map_lin = nullptr; }
 }
 
 #ifdef MFEM_USE_MPI
 FindPointsGSLIB::FindPointsGSLIB(MPI_Comm comm_)
-   : mesh(NULL),
-     fec_map_lin(NULL),
-     fdataD(NULL), cr(NULL), gsl_comm(NULL),
+   : mesh(nullptr),
+     fec_map_lin(nullptr),
+     fdataD(nullptr), cr(nullptr), gsl_comm(nullptr),
      dim(-1), points_cnt(-1), setupflag(false), default_interp_value(0),
      avgtype(AvgType::ARITHMETIC), bdr_tol(1e-8)
 {
@@ -147,10 +160,10 @@ FindPointsGSLIB::FindPointsGSLIB(MPI_Comm comm_)
    gf_rst_map.SetSize(4);
    for (int i = 0; i < mesh_split.Size(); i++)
    {
-      mesh_split[i] = NULL;
-      ir_split[i] = NULL;
-      fes_rst_map[i] = NULL;
-      gf_rst_map[i] = NULL;
+      mesh_split[i] = nullptr;
+      ir_split[i] = nullptr;
+      fes_rst_map[i] = nullptr;
+      gf_rst_map[i] = nullptr;
    }
 
    gsl_comm = new gslib::comm;
@@ -158,12 +171,21 @@ FindPointsGSLIB::FindPointsGSLIB(MPI_Comm comm_)
    comm_init(gsl_comm, comm_);
    crystal_init(cr, gsl_comm);
 }
+
+FindPointsGSLIB::FindPointsGSLIB(ParMesh &mesh_in, const double bb_t,
+                                 const double newt_tol, const int npt_max)
+   : FindPointsGSLIB(mesh_in.GetComm())
+{
+   Setup(mesh_in, bb_t, newt_tol, npt_max);
+}
 #endif
 
 void FindPointsGSLIB::Setup(Mesh &m, const double bb_t, const double newt_tol,
                             const int npt_max)
 {
    MFEM_VERIFY(m.GetNodes() != NULL, "Mesh nodes are required.");
+   MFEM_VERIFY(m.SpaceDimension() == m.Dimension(),
+               "Mesh spatial dimension and reference element dimension must be the same");
    const int meshOrder = m.GetNodes()->FESpace()->GetMaxElementOrder();
 
    // call FreeData if FindPointsGSLIB::Setup has been called already
@@ -171,37 +193,9 @@ void FindPointsGSLIB::Setup(Mesh &m, const double bb_t, const double newt_tol,
 
    mesh = &m;
    dim  = mesh->Dimension();
-   unsigned dof1D = meshOrder + 1;
+   const unsigned int dof1D = meshOrder+1;
 
-   SetupSplitMeshes();
-   if (dim == 2)
-   {
-      if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
-      ir_split[0] = new IntegrationRule(3*pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[0], ir_split[0], meshOrder);
-
-      if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
-      ir_split[1] = new IntegrationRule(pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[1], ir_split[1], meshOrder);
-   }
-   else if (dim == 3)
-   {
-      if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
-      ir_split[0] = new IntegrationRule(pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[0], ir_split[0], meshOrder);
-
-      if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
-      ir_split[1] = new IntegrationRule(4*pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[1], ir_split[1], meshOrder);
-
-      if (ir_split[2]) { delete ir_split[2]; ir_split[2] = NULL; }
-      ir_split[2] = new IntegrationRule(3*pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[2], ir_split[2], meshOrder);
-
-      if (ir_split[3]) { delete ir_split[3]; ir_split[3] = NULL; }
-      ir_split[3] = new IntegrationRule(8*pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[3], ir_split[3], meshOrder);
-   }
+   SetupSplitMeshesAndIntegrationRules(meshOrder);
 
    GetNodalValues(mesh->GetNodes(), gsl_mesh);
 
@@ -1128,13 +1122,18 @@ void FindPointsGSLIB::Interpolate(Mesh &m, const Vector &point_pos,
 void FindPointsGSLIB::FreeData()
 {
    if (!setupflag) { return; }
-   if (dim == 2)
+#ifdef MFEM_USE_MPI
+   if (!Mpi::IsFinalized())  // currently segfaults inside gslib otherwise
+#endif
    {
-      findpts_free_2((gslib::findpts_data_2 *)this->fdataD);
-   }
-   else
-   {
-      findpts_free_3((gslib::findpts_data_3 *)this->fdataD);
+      if (dim == 2)
+      {
+         findpts_free_2((gslib::findpts_data_2 *)this->fdataD);
+      }
+      else
+      {
+         findpts_free_3((gslib::findpts_data_3 *)this->fdataD);
+      }
    }
    gsl_code.DeleteAll();
    gsl_proc.DeleteAll();
@@ -1158,8 +1157,8 @@ void FindPointsGSLIB::FreeData()
 
 void FindPointsGSLIB::SetupSplitMeshes()
 {
-   fec_map_lin = new H1_FECollection(1, dim);
-   if (mesh->Dimension() == 2)
+   if (fec_map_lin == nullptr) { fec_map_lin = new H1_FECollection(1, dim); }
+   if (dim == 2)
    {
       int Nvert = 7;
       int NEsplit = 3;
@@ -1201,7 +1200,7 @@ void FindPointsGSLIB::SetupSplitMeshes()
       mesh_split[1] = new Mesh(Mesh::MakeCartesian2D(1, 1,
                                                      Element::QUADRILATERAL));
    }
-   else if (mesh->Dimension() == 3)
+   else if (dim == 3)
    {
       mesh_split[0] = new Mesh(Mesh::MakeCartesian3D(1, 1, 1,
                                                      Element::HEXAHEDRON));
@@ -1346,41 +1345,6 @@ void FindPointsGSLIB::SetupSplitMeshes()
          }
       }
    }
-
-   NE_split_total = 0;
-   split_element_map.SetSize(0);
-   split_element_index.SetSize(0);
-   int NEsplit = 0;
-   for (int e = 0; e < mesh->GetNE(); e++)
-   {
-      const Geometry::Type gt   = mesh->GetElement(e)->GetGeometryType();
-      if (gt == Geometry::TRIANGLE || gt == Geometry::PRISM)
-      {
-         NEsplit = 3;
-      }
-      else if (gt == Geometry::TETRAHEDRON)
-      {
-         NEsplit = 4;
-      }
-      else if (gt == Geometry::PYRAMID)
-      {
-         NEsplit = 8;
-      }
-      else if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
-      {
-         NEsplit = 1;
-      }
-      else
-      {
-         MFEM_ABORT("Unsupported geometry type.");
-      }
-      NE_split_total += NEsplit;
-      for (int i = 0; i < NEsplit; i++)
-      {
-         split_element_map.Append(e);
-         split_element_index.Append(i);
-      }
-   }
 }
 
 void FindPointsGSLIB::SetupIntegrationRuleForSplitMesh(Mesh *meshin,
@@ -1431,6 +1395,79 @@ void FindPointsGSLIB::SetupIntegrationRuleForSplitMesh(Mesh *meshin,
    }
 }
 
+void FindPointsGSLIB::SetupSplitMeshesAndIntegrationRules(const int order)
+{
+   MFEM_VERIFY(mesh, "Setup FindPointsGSLIB with mesh first.");
+   const int dof1D = order+1;
+   const int dim = mesh->Dimension();
+
+   SetupSplitMeshes();
+   if (dim == 2)
+   {
+      if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
+      ir_split[0] = new IntegrationRule(3*pow(dof1D, dim));
+      SetupIntegrationRuleForSplitMesh(mesh_split[0], ir_split[0], order);
+
+      if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
+      ir_split[1] = new IntegrationRule(pow(dof1D, dim));
+      SetupIntegrationRuleForSplitMesh(mesh_split[1], ir_split[1], order);
+   }
+   else if (dim == 3)
+   {
+      if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
+      ir_split[0] = new IntegrationRule(pow(dof1D, dim));
+      SetupIntegrationRuleForSplitMesh(mesh_split[0], ir_split[0], order);
+
+      if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
+      ir_split[1] = new IntegrationRule(4*pow(dof1D, dim));
+      SetupIntegrationRuleForSplitMesh(mesh_split[1], ir_split[1], order);
+
+      if (ir_split[2]) { delete ir_split[2]; ir_split[2] = NULL; }
+      ir_split[2] = new IntegrationRule(3*pow(dof1D, dim));
+      SetupIntegrationRuleForSplitMesh(mesh_split[2], ir_split[2], order);
+
+      if (ir_split[3]) { delete ir_split[3]; ir_split[3] = NULL; }
+      ir_split[3] = new IntegrationRule(8*pow(dof1D, dim));
+      SetupIntegrationRuleForSplitMesh(mesh_split[3], ir_split[3], order);
+   }
+
+   // Setup map for non tensor-product elements
+   NE_split_total = 0;
+   split_element_map.SetSize(0);
+   split_element_index.SetSize(0);
+   int NEsplit = 0;
+   for (int e = 0; e < mesh->GetNE(); e++)
+   {
+      const Geometry::Type gt   = mesh->GetElement(e)->GetGeometryType();
+      if (gt == Geometry::TRIANGLE || gt == Geometry::PRISM)
+      {
+         NEsplit = 3;
+      }
+      else if (gt == Geometry::TETRAHEDRON)
+      {
+         NEsplit = 4;
+      }
+      else if (gt == Geometry::PYRAMID)
+      {
+         NEsplit = 8;
+      }
+      else if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
+      {
+         NEsplit = 1;
+      }
+      else
+      {
+         MFEM_ABORT("Unsupported geometry type.");
+      }
+      NE_split_total += NEsplit;
+      for (int i = 0; i < NEsplit; i++)
+      {
+         split_element_map.Append(e);
+         split_element_index.Append(i);
+      }
+   }
+}
+
 void FindPointsGSLIB::GetNodalValues(const GridFunction *gf_in,
                                      Vector &node_vals)
 {
@@ -2081,6 +2118,19 @@ void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,
    } // parallel
 }
 
+Array<unsigned int> FindPointsGSLIB::GetPointsNotFoundIndices() const
+{
+   Array<unsigned int> nf_idxs;
+   for (int i = 0; i < gsl_code.Size(); i++)
+   {
+      if (gsl_code[i] == 2)
+      {
+         nf_idxs.Append(i);
+      }
+   }
+   return nf_idxs;
+}
+
 void FindPointsGSLIB::DistributePointInfoToOwningMPIRanks(
    Array<unsigned int> &recv_elem, Vector &recv_ref,
    Array<unsigned int> &recv_code)
@@ -2386,6 +2436,10 @@ void OversetFindPointsGSLIB::Setup(Mesh &m, const int meshid,
 {
    MFEM_VERIFY(m.GetNodes() != NULL, "Mesh nodes are required.");
    const int meshOrder = m.GetNodes()->FESpace()->GetMaxElementOrder();
+   const int gfOrder = gfmax ? gfmax->FESpace()->GetMaxElementOrder() :
+                       meshOrder;
+   MFEM_VERIFY(meshOrder == gfOrder,
+               "Mesh order must match gfmax order in OversetFindPointsGSLIB.");
 
    // FreeData if OversetFindPointsGSLIB::Setup has been called already
    if (setupflag) { FreeData(); }
@@ -2395,35 +2449,7 @@ void OversetFindPointsGSLIB::Setup(Mesh &m, const int meshid,
    const FiniteElement *fe = mesh->GetNodalFESpace()->GetTypicalFE();
    unsigned dof1D = fe->GetOrder() + 1;
 
-   SetupSplitMeshes();
-   if (dim == 2)
-   {
-      if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
-      ir_split[0] = new IntegrationRule(3*pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[0], ir_split[0], meshOrder);
-
-      if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
-      ir_split[1] = new IntegrationRule(pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[1], ir_split[1], meshOrder);
-   }
-   else if (dim == 3)
-   {
-      if (ir_split[0]) { delete ir_split[0]; ir_split[0] = NULL; }
-      ir_split[0] = new IntegrationRule(pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[0], ir_split[0], meshOrder);
-
-      if (ir_split[1]) { delete ir_split[1]; ir_split[1] = NULL; }
-      ir_split[1] = new IntegrationRule(4*pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[1], ir_split[1], meshOrder);
-
-      if (ir_split[2]) { delete ir_split[2]; ir_split[2] = NULL; }
-      ir_split[2] = new IntegrationRule(3*pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[2], ir_split[2], meshOrder);
-
-      if (ir_split[3]) { delete ir_split[3]; ir_split[3] = NULL; }
-      ir_split[3] = new IntegrationRule(8*pow(dof1D, dim));
-      SetupIntegrationRuleForSplitMesh(mesh_split[3], ir_split[3], meshOrder);
-   }
+   SetupSplitMeshesAndIntegrationRules(meshOrder);
 
    GetNodalValues(mesh->GetNodes(), gsl_mesh);
 
@@ -2480,7 +2506,7 @@ void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
 {
    MFEM_VERIFY(setupflag, "Use OversetFindPointsGSLIB::Setup before "
                "finding points.");
-   MFEM_VERIFY(overset, "Please setup FindPoints for overlapping grids.");
+   MFEM_VERIFY(overset, "Please use OversetFindPoints for overlapping grids.");
    points_cnt = point_pos.Size() / dim;
    unsigned int match = 0; // Don't find points in the mesh if point_id=mesh_id
 

fem/gslib.hpp

@@ -13,7 +13,11 @@
 #define MFEM_GSLIB
 
 #include "../config/config.hpp"
+#ifdef MFEM_USE_MPI
+#include "pgridfunc.hpp"
+#else
 #include "gridfunc.hpp"
+#endif
 
 #ifdef MFEM_USE_GSLIB
 
@@ -131,6 +135,10 @@ protected:
                                                  IntegrationRule *irule,
                                                  int order);
 
+   /// Helper function that calls \ref SetupSplitMeshes and
+   /// \ref SetupIntegrationRuleForSplitMesh.
+   virtual void SetupSplitMeshesAndIntegrationRules(const int order);
+
    /// Get GridFunction value at the points expected by GSLIB.
    virtual void GetNodalValues(const GridFunction *gf_in, Vector &node_vals);
 
@@ -190,14 +198,23 @@ protected:
    void InterpolateOnDevice(const Vector &field_in_evec, Vector &field_out,
                             const int nel, const int ncomp,
                             const int dof1dsol, const int ordering);
+
 public:
    FindPointsGSLIB();
+   FindPointsGSLIB(Mesh &mesh_in, const double bb_t = 0.1,
+                   const double newt_tol = 1.0e-12,
+                   const int npt_max = 256);
 
 #ifdef MFEM_USE_MPI
    FindPointsGSLIB(MPI_Comm comm_);
+   FindPointsGSLIB(ParMesh &mesh_in, const double bb_t = 0.1,
+                   const double newt_tol = 1.0e-12,
+                   const int npt_max = 256);
 #endif
 
    virtual ~FindPointsGSLIB();
+   FindPointsGSLIB(const FindPointsGSLIB&) = delete;
+   FindPointsGSLIB& operator=(const FindPointsGSLIB&) = delete;
 
    /** Initializes the internal mesh in gslib, by sending the positions of the
        Gauss-Lobatto nodes of the input Mesh object \p m.
@@ -212,8 +229,8 @@ public:
        @param[in] npt_max   (Optional) Number of points for simultaneous
                             iteration. This alters performance and
                             memory footprint.*/
-   void Setup(Mesh &m, const double bb_t = 0.1,
-              const double newt_tol = 1.0e-12,
+
+   void Setup(Mesh &m, const double bb_t = 0.1, const double newt_tol = 1.0e-12,
               const int npt_max = 256);
    /** Searches positions given in physical space by \p point_pos.
        These positions can be ordered byNodes: (XXX...,YYY...,ZZZ) or
@@ -289,7 +306,12 @@ public:
 
    /** Cleans up memory allocated internally by gslib.
        Note that in parallel, this must be called before MPI_Finalize(), as it
-       calls MPI_Comm_free() for internal gslib communicators. */
+       calls MPI_Comm_free() for internal gslib communicators. FreeData is
+       also called by the class destructor and there are no memory leaks if the
+       destructor is called before MPI_Finalize(). If the destructor is called
+       after MPI_Finalize(), there will be an error because gslib will try to
+       invoke some MPI functions.
+   */
    virtual void FreeData();
 
    /// Return code for each point searched by FindPoints: inside element (0), on
@@ -312,6 +334,9 @@ public:
    /// point found by FindPoints.
    virtual const Vector &GetGSLIBReferencePosition() const { return gsl_ref; }
 
+   /// Get array of indices of not-found points.
+   Array<unsigned int> GetPointsNotFoundIndices() const;
+
    /** @name Methods to support a custom interpolation procedure.
        \brief The physical-space point that the user seeks to interpolate at
        could be located inside an element on another mpi rank.

fem/hyperbolic.cpp

@@ -181,7 +181,7 @@ void HyperbolicFormIntegrator::AssembleFaceVector(
    // current elements' the number of degrees of freedom
    // does not consider the number of equations
    const int dof1 = el1.GetDof();
-   const int dof2 = el2.GetDof();
+   const int dof2 = (Tr.Elem2No >= 0)?(el2.GetDof()):(0);
 
 #ifdef MFEM_THREAD_SAFE
    // Local storage for element integration
@@ -219,7 +219,9 @@ void HyperbolicFormIntegrator::AssembleFaceVector(
    const IntegrationRule *ir = IntRule;
    if (!ir)
    {
-      const int order = 2*std::max(el1.GetOrder(), el2.GetOrder()) + IntOrderOffset;
+      const int max_el_order = dof2 ? std::max(el1.GetOrder(),
+                                               el2.GetOrder()) : el1.GetOrder();
+      const int order = 2*max_el_order + IntOrderOffset;
       ir = &IntRules.Get(Tr.GetGeometryType(), order);
    }
    // loop over integration points
@@ -231,18 +233,22 @@ void HyperbolicFormIntegrator::AssembleFaceVector(
 
       // Calculate basis functions on both elements at the face
       el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
-      el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
 
       // Interpolate elfun at the point
       elfun1_mat.MultTranspose(shape1, state1);
-      elfun2_mat.MultTranspose(shape2, state2);
+
+      if (dof2)
+      {
+         // Calculate basis functions on both elements at the face
+         el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
+         // Interpolate elfun at the point
+         elfun2_mat.MultTranspose(shape2, state2);
+      }
 
       // Get the normal vector and the flux on the face
       if (nor.Size() == 1)  // if 1D, use 1 or -1.
       {
-         // This assume the 1D integration point is in (0,1). This may not work
-         // if this changes.
-         nor(0) = (Tr.GetElement1IntPoint().x - 0.5) * 2.0;
+         nor(0) = 2*Tr.GetElement1IntPoint().x - 1.;
       }
       else
       {
@@ -250,14 +256,18 @@ void HyperbolicFormIntegrator::AssembleFaceVector(
       }
       // Compute F(u+, x) and F(u-, x) with maximum characteristic speed
       // Compute hat(F) using evaluated quantities
-      const real_t speed = numFlux.Eval(state1, state2, nor, Tr, fluxN);
+      const real_t speed = (dof2) ? numFlux.Eval(state1, state2, nor, Tr, fluxN):
+                           fluxFunction.ComputeFluxDotN(state1, nor, Tr, fluxN);
 
       // Update the global max char speed
       max_char_speed = std::max(speed, max_char_speed);
 
       // pre-multiply integration weight to flux
       AddMult_a_VWt(-ip.weight*sign, shape1, fluxN, elvect1_mat);
-      AddMult_a_VWt(+ip.weight*sign, shape2, fluxN, elvect2_mat);
+      if (dof2)
+      {
+         AddMult_a_VWt(+ip.weight*sign, shape2, fluxN, elvect2_mat);
+      }
    }
 }
 
@@ -268,7 +278,7 @@ void HyperbolicFormIntegrator::AssembleFaceGrad(
    // current elements' the number of degrees of freedom
    // does not consider the number of equations
    const int dof1 = el1.GetDof();
-   const int dof2 = el2.GetDof();
+   const int dof2 = (Tr.Elem2No >= 0)?(el2.GetDof()):(0);
 
 #ifdef MFEM_THREAD_SAFE
    // Local storage for element integration
@@ -302,7 +312,9 @@ void HyperbolicFormIntegrator::AssembleFaceGrad(
    const IntegrationRule *ir = IntRule;
    if (!ir)
    {
-      const int order = 2*std::max(el1.GetOrder(), el2.GetOrder()) + IntOrderOffset;
+      const int max_el_order = dof2 ? std::max(el1.GetOrder(),
+                                               el2.GetOrder()) : el1.GetOrder();
+      const int order = 2*max_el_order + IntOrderOffset;
       ir = &IntRules.Get(Tr.GetGeometryType(), order);
    }
    // loop over integration points
@@ -312,20 +324,25 @@ void HyperbolicFormIntegrator::AssembleFaceGrad(
 
       Tr.SetAllIntPoints(&ip); // set face and element int. points
 
-      // Calculate basis functions on both elements at the face
+      // Calculate basis functions of the first element at the face
       el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
-      el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
 
       // Interpolate elfun at the point
       elfun1_mat.MultTranspose(shape1, state1);
-      elfun2_mat.MultTranspose(shape2, state2);
+
+      if (dof2)
+      {
+         // Calculate basis function of the second element at the face
+         el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
+
+         // Interpolate elfun at the point
+         elfun2_mat.MultTranspose(shape2, state2);
+      }
 
       // Get the normal vector and the flux on the face
       if (nor.Size() == 1)  // if 1D, use 1 or -1.
       {
-         // This assume the 1D integration point is in (0,1). This may not work
-         // if this changes.
-         nor(0) = (Tr.GetElement1IntPoint().x - 0.5) * 2.0;
+         nor(0) = 2*Tr.GetElement1IntPoint().x - 1.;
       }
       else
       {
@@ -335,7 +352,14 @@ void HyperbolicFormIntegrator::AssembleFaceGrad(
       // Trial side 1
 
       // Compute hat(J) using evaluated quantities
-      numFlux.Grad(1, state1, state2, nor, Tr, JDotN);
+      if (dof2)
+      {
+         numFlux.Grad(1, state1, state2, nor, Tr, JDotN);
+      }
+      else
+      {
+         fluxFunction.ComputeFluxJacobianDotN(state1, nor, Tr, JDotN);
+      }
 
       const int ioff = fluxFunction.num_equations * dof1;
 
@@ -360,36 +384,325 @@ void HyperbolicFormIntegrator::AssembleFaceGrad(
             }
          }
 
-      // Trial side 2
+      if (dof2)
+      {
+         // Trial side 2
+
+         // Compute hat(J) using evaluated quantities
+         numFlux.Grad(2, state1, state2, nor, Tr, JDotN);
+
+         const int joff = ioff;
+
+         for (int di = 0; di < fluxFunction.num_equations; di++)
+            for (int dj = 0; dj < fluxFunction.num_equations; dj++)
+            {
+               // pre-multiply integration weight to Jacobian
+               const real_t w = +ip.weight * sign * JDotN(di,dj);
+               for (int j = 0; j < dof2; j++)
+               {
+                  // Test side 1
+                  for (int i = 0; i < dof1; i++)
+                  {
+                     elmat(i+dof1*di, joff+j+dof2*dj) += w * shape1(i) * shape2(j);
+                  }
+
+                  // Test side 2
+                  for (int i = 0; i < dof2; i++)
+                  {
+                     elmat(ioff+i+dof2*di, joff+j+dof2*dj) -= w * shape2(i) * shape2(j);
+                  }
+               }
+            }
+      }
+   }
+}
+
+BdrHyperbolicDirichletIntegrator::BdrHyperbolicDirichletIntegrator(
+   const NumericalFlux &numFlux,
+   VectorCoefficient &bdrState,
+   const int IntOrderOffset,
+   real_t sign)
+   : NonlinearFormIntegrator(),
+     numFlux(numFlux),
+     fluxFunction(numFlux.GetFluxFunction()),
+     u_vcoeff(bdrState),
+     IntOrderOffset(IntOrderOffset),
+     sign(sign),
+     num_equations(fluxFunction.num_equations)
+{
+   MFEM_VERIFY(fluxFunction.num_equations == bdrState.GetVDim(),
+               "Flux function does not match the vector dimension of the coefficient!");
+#ifndef MFEM_THREAD_SAFE
+   state_in.SetSize(num_equations);
+   state_out.SetSize(num_equations);
+   fluxN.SetSize(num_equations);
+   JDotN.SetSize(num_equations);
+   nor.SetSize(fluxFunction.dim);
+#endif
+   ResetMaxCharSpeed();
+}
+
+void BdrHyperbolicDirichletIntegrator::AssembleFaceVector(
+   const FiniteElement &el, const FiniteElement &,
+   FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
+{
+   MFEM_ASSERT(Tr.Elem2No < 0, "Not a boundary face!");
+
+   // current elements' the number of degrees of freedom
+   // does not consider the number of equations
+   const int dof = el.GetDof();
+
+#ifdef MFEM_THREAD_SAFE
+   // Local storage for element integration
+
+   // shape function value at an integration point
+   Vector shape(dof);
+   // normal vector (usually not a unit vector)
+   Vector nor(Tr.GetSpaceDim());
+   // state value at an integration point - interior
+   Vector state_in(num_equations);
+   // state value at an integration point - boundary
+   Vector state_out(num_equations);
+   // hat(F)(u,x)
+   Vector fluxN(num_equations);
+#else
+   shape.SetSize(dof);
+#endif
+
+   elvect.SetSize(dof * num_equations);
+   elvect = 0.0;
+
+   const DenseMatrix elfun_mat(elfun.GetData(), dof, num_equations);
+
+   DenseMatrix elvect_mat(elvect.GetData(), dof, num_equations);
+
+   // Obtain integration rule. If integration is rule is given, then use it.
+   // Otherwise, get (2*p + IntOrderOffset) order integration rule
+   const IntegrationRule *ir = IntRule;
+   if (!ir)
+   {
+      const int order = 2*el.GetOrder() + IntOrderOffset;
+      ir = &IntRules.Get(Tr.GetGeometryType(), order);
+   }
+   // loop over integration points
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      Tr.SetAllIntPoints(&ip); // set face and element int. points
+
+      // Calculate basis functions at the face
+      el.CalcShape(Tr.GetElement1IntPoint(), shape);
+
+      // Interpolate elfun at the point
+      elfun_mat.MultTranspose(shape, state_in);
+
+      // Evaluate boundary state at the point
+      u_vcoeff.Eval(state_out, Tr, ip);
+
+      // Get the normal vector and the flux on the face
+      if (nor.Size() == 1)  // if 1D, use 1 or -1.
+      {
+         nor(0) = 2*Tr.GetElement1IntPoint().x - 1.;
+      }
+      else
+      {
+         CalcOrtho(Tr.Jacobian(), nor);
+      }
+      // Compute F(u+, x) and F(u_b, x) with maximum characteristic speed
+      // Compute hat(F) using evaluated quantities
+      const real_t speed = numFlux.Eval(state_in, state_out, nor, Tr, fluxN);
+
+      // Update the global max char speed
+      max_char_speed = std::max(speed, max_char_speed);
+
+      // pre-multiply integration weight to flux
+      AddMult_a_VWt(-ip.weight*sign, shape, fluxN, elvect_mat);
+   }
+}
+
+void BdrHyperbolicDirichletIntegrator::AssembleFaceGrad(
+   const FiniteElement &el, const FiniteElement &,
+   FaceElementTransformations &Tr, const Vector &elfun, DenseMatrix &elmat)
+{
+   // current elements' the number of degrees of freedom
+   // does not consider the number of equations
+   const int dof = el.GetDof();
+
+#ifdef MFEM_THREAD_SAFE
+   // Local storage for element integration
+
+   // shape function value at an integration point
+   Vector shape(dof);
+   // normal vector (usually not a unit vector)
+   Vector nor(Tr.GetSpaceDim());
+   // state value at an integration point - interior
+   Vector state_in(num_equations);
+   // state value at an integration point - boundary
+   Vector state_out(num_equations);
+   // hat(J)(u,x)
+   DenseMatrix JDotN(num_equations);
+#else
+   shape.SetSize(dof);
+#endif
+
+   elmat.SetSize(dof * num_equations);
+   elmat = 0.0;
+
+   const DenseMatrix elfun_mat(elfun.GetData(), dof, num_equations);
+
+   // Obtain integration rule. If integration is rule is given, then use it.
+   // Otherwise, get (2*p + IntOrderOffset) order integration rule
+   const IntegrationRule *ir = IntRule;
+   if (!ir)
+   {
+      const int order = 2*el.GetOrder() + IntOrderOffset;
+      ir = &IntRules.Get(Tr.GetGeometryType(), order);
+   }
+   // loop over integration points
+   for (int q = 0; q < ir->GetNPoints(); q++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(q);
+
+      Tr.SetAllIntPoints(&ip); // set face and element int. points
+
+      // Calculate basis functions at the face
+      el.CalcShape(Tr.GetElement1IntPoint(), shape);
+
+      // Interpolate elfun at the point
+      elfun_mat.MultTranspose(shape, state_in);
+
+      // Evaluate boundary state at the point
+      u_vcoeff.Eval(state_out, Tr, ip);
+
+      // Get the normal vector and the flux on the face
+      if (nor.Size() == 1)  // if 1D, use 1 or -1.
+      {
+         nor(0) = 2*Tr.GetElement1IntPoint().x - 1.;
+      }
+      else
+      {
+         CalcOrtho(Tr.Jacobian(), nor);
+      }
 
       // Compute hat(J) using evaluated quantities
-      numFlux.Grad(2, state1, state2, nor, Tr, JDotN);
-
-      const int joff = ioff;
+      numFlux.Grad(1, state_in, state_out, nor, Tr, JDotN);
 
       for (int di = 0; di < fluxFunction.num_equations; di++)
          for (int dj = 0; dj < fluxFunction.num_equations; dj++)
          {
             // pre-multiply integration weight to Jacobian
-            const real_t w = +ip.weight * sign * JDotN(di,dj);
-            for (int j = 0; j < dof2; j++)
-            {
-               // Test side 1
-               for (int i = 0; i < dof1; i++)
+            const real_t w = -ip.weight * sign * JDotN(di,dj);
+            for (int j = 0; j < dof; j++)
+               for (int i = 0; i < dof; i++)
                {
-                  elmat(i+dof1*di, joff+j+dof2*dj) += w * shape1(i) * shape2(j);
+                  elmat(i+dof*di, j+dof*dj) += w * shape(i) * shape(j);
                }
-
-               // Test side 2
-               for (int i = 0; i < dof2; i++)
-               {
-                  elmat(ioff+i+dof2*di, joff+j+dof2*dj) -= w * shape2(i) * shape2(j);
-               }
-            }
          }
    }
 }
 
+BoundaryHyperbolicFlowIntegrator::BoundaryHyperbolicFlowIntegrator(
+   const FluxFunction &flux, VectorCoefficient &u, real_t alpha_, real_t beta_,
+   const int IntOrderOffset_)
+   : fluxFunction(flux), u_vcoeff(u), alpha(alpha_), beta(beta_),
+     IntOrderOffset(IntOrderOffset_)
+{
+   MFEM_VERIFY(fluxFunction.num_equations == u_vcoeff.GetVDim(),
+               "Flux function does not match the vector dimension of the coefficient!");
+#ifndef MFEM_THREAD_SAFE
+   state.SetSize(fluxFunction.num_equations);
+   nor.SetSize(fluxFunction.dim);
+   fluxN.SetSize(fluxFunction.num_equations);
+#endif
+   ResetMaxCharSpeed();
+}
+
+void BoundaryHyperbolicFlowIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
+{
+   mfem_error("BoundaryHyperbolicFlowIntegrator::AssembleRHSElementVect\n"
+              "  is not implemented as boundary integrator!\n"
+              "  Use LinearForm::AddBdrFaceIntegrator instead of\n"
+              "  LinearForm::AddBoundaryIntegrator.");
+}
+
+void BoundaryHyperbolicFlowIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
+{
+   // current elements' the number of degrees of freedom
+   // does not consider the number of equations
+   const int dof = el.GetDof();
+
+#ifdef MFEM_THREAD_SAFE
+   // Local storage for element integration
+
+   // shape function value at an integration point
+   Vector shape(dof);
+   // state value at an integration point
+   Vector state(fluxFunction.num_equations);
+   // normal vector (usually not a unit vector)
+   Vector nor(Tr.GetSpaceDim());
+   // hat(F)(u,x)
+   Vector fluxN(fluxFunction.num_equations);
+#else
+   shape.SetSize(dof);
+#endif
+
+   elvect.SetSize(dof * fluxFunction.num_equations);
+   elvect = 0.0;
+
+   DenseMatrix elvect_mat(elvect.GetData(), dof, fluxFunction.num_equations);
+
+   // Obtain integration rule. If integration is rule is given, then use it.
+   // Otherwise, get (2*p + IntOrderOffset) order integration rule
+   const IntegrationRule *ir = IntRule;
+   if (!ir)
+   {
+      const int order = 2*el.GetOrder() + IntOrderOffset;
+      ir = &IntRules.Get(Tr.GetGeometryType(), order);
+   }
+   // loop over integration points
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      Tr.SetAllIntPoints(&ip); // set face and element int. points
+
+      // Calculate basis functions on both elements at the face
+      el.CalcShape(Tr.GetElement1IntPoint(), shape);
+
+      // Evaluate the coefficient at the point
+      u_vcoeff.Eval(state, Tr, ip);
+
+      // Get the normal vector and the flux on the face
+      if (nor.Size() == 1)  // if 1D, use 1 or -1.
+      {
+         nor(0) = 2*Tr.GetElement1IntPoint().x - 1.;
+      }
+      else
+      {
+         CalcOrtho(Tr.Jacobian(), nor);
+      }
+      // Compute F(u, x) with maximum characteristic speed
+      const real_t speed = fluxFunction.ComputeFluxDotN(state, nor, Tr, fluxN);
+
+      // Update the global max char speed
+      max_char_speed = std::max(speed, max_char_speed);
+
+      // pre-multiply integration weight to flux
+      const real_t a = 0.5 * alpha  * ip.weight;
+      const real_t b = beta * ip.weight;
+
+      for (int n = 0; n < fluxFunction.num_equations; n++)
+      {
+         fluxN(n) = a * fluxN(n) - b * fabs(fluxN(n));
+      }
+
+      AddMultVWt(shape, fluxN, elvect_mat);
+   }
+}
+
 real_t FluxFunction::ComputeFluxDotN(const Vector &U,
                                      const Vector &normal,
                                      FaceElementTransformations &Tr,

fem/hyperbolic.hpp

@@ -306,12 +306,14 @@ MFEM_DEPRECATED typedef NumericalFlux RiemannSolver;
 class HyperbolicFormIntegrator : public NonlinearFormIntegrator
 {
 private:
-   // The maximum characteristic speed, updated during element/face vector assembly
-   real_t max_char_speed;
    const NumericalFlux &numFlux;   // Numerical flux that maps F(u±,x) to F̂
    const FluxFunction &fluxFunction;
    const int IntOrderOffset; // integration order offset, 2*p + IntOrderOffset.
    const real_t sign;
+
+   // The maximum characteristic speed, updated during element/face vector assembly
+   real_t max_char_speed;
+
 #ifndef MFEM_THREAD_SAFE
    // Local storage for element integration
    Vector shape;              // shape function value at an integration point
@@ -331,8 +333,9 @@ private:
 
 public:
    const int num_equations;  // the number of equations
+
    /**
-    * @brief Construct a new Hyperbolic Form Integrator object
+    * @brief Construct a new HyperbolicFormIntegrator object
     *
     * @param[in] numFlux numerical flux
     * @param[in] IntOrderOffset integration order offset
@@ -343,21 +346,14 @@ public:
       const int IntOrderOffset = 0,
       const real_t sign = 1.);
 
-   /**
-    * @brief Reset the Max Char Speed 0
-    *
-    */
-   void ResetMaxCharSpeed()
-   {
-      max_char_speed = 0.0;
-   }
+   /// Reset the maximum characteristic speed to zero
+   void ResetMaxCharSpeed() { max_char_speed = 0.0; }
 
-   real_t GetMaxCharSpeed()
-   {
-      return max_char_speed;
-   }
+   /// Get the maximum characteristic speed
+   real_t GetMaxCharSpeed() const { return max_char_speed; }
 
-   const FluxFunction &GetFluxFunction() { return fluxFunction; }
+   /// Get the associated flux function
+   const FluxFunction &GetFluxFunction() const { return fluxFunction; }
 
    /**
     * @brief Implements (F(u), ∇v) with abstract F computed by
@@ -416,6 +412,182 @@ public:
                          const Vector &elfun, DenseMatrix &elmat) override;
 };
 
+/**
+ * @brief Abstract boundary hyperbolic form integrator, assembling
+ * <F̂(u⁻,u_b,x) n, [v]> term for scalar finite elements at the boundary.
+ *
+ * This form integrator is coupled with a NumericalFlux that implements the
+ * numerical flux F̂ at the boundary faces. The flux F is obtained from the
+ * FluxFunction assigned to the aforementioned NumericalFlux with the given
+ * boundary coefficient for the state u_b.
+ *
+ * Note the class can be used for imposing conditions on interior interfaces.
+ */
+class BdrHyperbolicDirichletIntegrator : public NonlinearFormIntegrator
+{
+private:
+   const NumericalFlux &numFlux;    // Numerical flux that maps F to F̂
+   const FluxFunction &fluxFunction;
+   VectorCoefficient &u_vcoeff;     // Boundary state vector coefficient
+   const int IntOrderOffset; // integration order offset, 2*p + IntOrderOffset.
+   const real_t sign;
+
+   // The maximum characteristic speed, updated during element/face vector assembly
+   real_t max_char_speed;
+
+#ifndef MFEM_THREAD_SAFE
+   // Local storage for element integration
+   Vector shape;  // shape function value at an integration point
+   Vector state_in;  // state value at an integration point - interior
+   Vector state_out;  // state value at an integration point - boundary
+   Vector nor;     // normal vector, see mfem::CalcOrtho()
+   Vector fluxN;   // F̂(u⁻,u_b,x) n
+   DenseMatrix JDotN;   // Ĵ(u⁻,u_b,x) n
+#endif
+
+public:
+   const int num_equations;  // the number of equations
+
+   /**
+    * @brief Construct a new BdrHyperbolicDirichletIntegrator object
+    *
+    * @param[in] numFlux numerical flux
+    * @param[in] bdrState boundary state coefficient
+    * @param[in] IntOrderOffset integration order offset
+    * @param[in] sign sign of the convection term
+    */
+   BdrHyperbolicDirichletIntegrator(
+      const NumericalFlux &numFlux,
+      VectorCoefficient &bdrState,
+      const int IntOrderOffset = 0,
+      const real_t sign = 1.);
+
+   /// Reset the maximum characteristic speed to zero
+   void ResetMaxCharSpeed() { max_char_speed = 0.0; }
+
+   /// Get the maximum characteristic speed
+   real_t GetMaxCharSpeed() const { return max_char_speed; }
+
+   /// Get the associated flux function
+   const FluxFunction &GetFluxFunction() const { return fluxFunction; }
+
+   /**
+    * @brief Implements <-F̂(u⁻,u_b,x) n, [v]> with abstract F̂ computed by
+    * NumericalFlux::Eval() of the numerical flux object
+    *
+    * @param[in] el1 finite element of the interior element
+    * @param[in] el2 not used
+    * @param[in] Tr face element transformations
+    * @param[in] elfun local coefficient of basis for the interior element
+    * @param[out] elvect evaluated dual vector <-F̂(u⁻,u_b,x) n, [v]>
+    */
+   void AssembleFaceVector(const FiniteElement &el1,
+                           const FiniteElement &el2,
+                           FaceElementTransformations &Tr,
+                           const Vector &elfun, Vector &elvect) override;
+
+   /**
+    * @brief Implements <-Ĵ(u⁻,u_b,x) n, [v]> with abstract Ĵ computed by
+    * NumericalFlux::Grad() of the numerical flux object
+    *
+    * @param[in] el1 finite element of the interior element
+    * @param[in] el2 not used
+    * @param[in] Tr face element transformations
+    * @param[in] elfun local coefficient of basis for the interior element
+    * @param[out] elmat evaluated Jacobian matrix <-Ĵ(u⁻,u_b,x) n, [v]>
+    */
+   void AssembleFaceGrad(const FiniteElement &el1,
+                         const FiniteElement &el2,
+                         FaceElementTransformations &Tr,
+                         const Vector &elfun, DenseMatrix &elmat) override;
+};
+
+/**
+ * @brief Abstract boundary hyperbolic linear form integrator, assembling
+ * <ɑ/2 F(u,x) n - β |F(u,x) n|, v> terms for scalar finite elements.
+ *
+ * This form integrator is coupled with a FluxFunction that evaluates the
+ * flux F at the boundary.
+ *
+ * Note the upwinding is performed component-wise. For general boundary
+ * integration with a numerical flux, see BdrHyperbolicDirichletIntegrator.
+ */
+class BoundaryHyperbolicFlowIntegrator : public LinearFormIntegrator
+{
+   const FluxFunction &fluxFunction;
+   VectorCoefficient &u_vcoeff;
+   const real_t alpha, beta;
+   const int IntOrderOffset; // integration order offset, 2*p + IntOrderOffset.
+
+   // The maximum characteristic speed, updated during face vector assembly
+   real_t max_char_speed;
+
+#ifndef MFEM_THREAD_SAFE
+   // Local storage for element integration
+   Vector shape;   // shape function value at an integration point
+   Vector state;   // state value at an integration point
+   Vector nor;     // normal vector, see mfem::CalcOrtho()
+   Vector fluxN;   // F(u,x) n
+#endif
+
+public:
+   /**
+    * @brief Construct a new BoundaryHyperbolicFlowIntegrator object
+    *
+    * @param[in] flux flux function
+    * @param[in] u vector state coefficient
+    * @param[in] alpha ɑ coefficient (β = ɑ/2)
+    * @param[in] IntOrderOffset integration order offset
+    */
+   BoundaryHyperbolicFlowIntegrator(const FluxFunction &flux, VectorCoefficient &u,
+                                    real_t alpha = -1., int IntOrderOffset = 0)
+      : BoundaryHyperbolicFlowIntegrator(flux, u, alpha, alpha/2., IntOrderOffset) { }
+
+   /**
+    * @brief Construct a new BoundaryHyperbolicFlowIntegrator object
+    *
+    * @param[in] flux flux function
+    * @param[in] u vector state coefficient
+    * @param[in] alpha ɑ coefficient
+    * @param[in] beta β coefficient
+    * @param[in] IntOrderOffset integration order offset
+    */
+   BoundaryHyperbolicFlowIntegrator(const FluxFunction &flux, VectorCoefficient &u,
+                                    real_t alpha, real_t beta, int IntOrderOffset = 0);
+
+   /// Reset the maximum characteristic speed to zero
+   void ResetMaxCharSpeed() { max_char_speed = 0.0; }
+
+   /// Get the maximum characteristic speed
+   real_t GetMaxCharSpeed() const { return max_char_speed; }
+
+   /// Get the associated flux function
+   const FluxFunction &GetFluxFunction() const { return fluxFunction; }
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+
+   /**
+    * @warning Boundary element integration not implemented, use
+    * AssembleRHSElementVect(const FiniteElement&,
+    * FaceElementTransformations &, Vector &) instead
+    */
+   void AssembleRHSElementVect(const FiniteElement &el,
+                               ElementTransformation &Tr,
+                               Vector &elvect) override;
+
+   /**
+    * @brief Implements <-F(u,x) n, v> with abstract F computed by
+    * FluxFunction::ComputeFluxDotN() of the flux function object
+    *
+    * @param[in] el finite element
+    * @param[in] Tr face element transformations
+    * @param[out] elvect evaluated dual vector <F(u,x) n, v>
+    */
+   void AssembleRHSElementVect(const FiniteElement &el,
+                               FaceElementTransformations &Tr,
+                               Vector &elvect) override;
+};
+
 
 /**
  * @brief Rusanov flux, also known as local Lax-Friedrichs,

fem/integ/bilininteg_dgtrace_pa.cpp

@@ -139,8 +139,6 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
                          Device::GetDeviceMemoryType() : pa_mt;
 
-   nf = fes.GetNFbyType(type);
-   if (nf==0) { return; }
    // Assumes tensor-product elements
    Mesh *mesh = fes.GetMesh();
    const FiniteElement &el = *fes.GetTypicalTraceElement();
@@ -148,6 +146,12 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
                                IntRule:
                                &GetRule(el.GetGeomType(), el.GetOrder(),
                                         *mesh->GetTypicalElementTransformation());
+
+
+   FaceQuadratureSpace qs(*mesh, *ir, type);
+   nf = qs.GetNumFaces();
+   if (nf==0) { return; }
+
    const int symmDims = 4;
    nq = ir->GetNPoints();
    dim = mesh->Dimension();
@@ -159,8 +163,6 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
    pa_data.SetSize(symmDims * nq * nf, Device::GetMemoryType());
-
-   FaceQuadratureSpace qs(*mesh, *ir, type);
    CoefficientVector vel(*u, qs, CoefficientStorage::COMPRESSED);
 
    CoefficientVector r(qs, CoefficientStorage::COMPRESSED);

fem/qspace.cpp

@@ -22,7 +22,7 @@ QuadratureSpaceBase::QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
 {
    for (int g = 0; g < Geometry::NumGeom; g++)
    {
-      int_rule[g] = NULL;
+      int_rule[g] = nullptr;
    }
    int_rule[geom] = &ir;
 }
@@ -37,6 +37,29 @@ void QuadratureSpaceBase::ConstructIntRules(int dim)
    }
 }
 
+const Array<int> &QuadratureSpaceBase::Offsets(
+   QSpaceOffsetStorage storage) const
+{
+   if (storage == QSpaceOffsetStorage::COMPRESSED || offsets.Size() > 1)
+   {
+      return offsets;
+   }
+   else
+   {
+      if (full_offset_cache.Size() == 0)
+      {
+         const int nq = size / ne;
+         full_offset_cache.SetSize(ne + 1);
+         int *d_full_offset_cache = full_offset_cache.Write();
+         mfem::forall(ne + 1, [=] MFEM_HOST_DEVICE (int e)
+         {
+            d_full_offset_cache[e] = nq * e;
+         });
+      }
+      return full_offset_cache;
+   }
+}
+
 namespace
 {
 
@@ -116,7 +139,7 @@ void QuadratureSpace::ConstructOffsets()
       {
          offsets[i] = offset;
          const Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-         MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule.");
+         MFEM_ASSERT(int_rule[geom] != nullptr, "Missing integration rule.");
          offset += int_rule[geom]->GetNPoints();
       }
       offsets[num_elem] = offset;
@@ -181,9 +204,9 @@ const Vector &QuadratureSpace::GetGeometricFactorWeights() const
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, int order_,
                                          FaceType face_type_)
-   : QuadratureSpaceBase(mesh_, order_),
-     face_type(face_type_),
-     num_faces(mesh.GetNFbyType(face_type))
+   : QuadratureSpaceBase(mesh_, order_), face_type(face_type_),
+     face_indices(mesh.GetFaceIndices(face_type_)),
+     face_indices_inv(mesh.GetInvFaceIndices(face_type_))
 {
    Construct();
 }
@@ -192,7 +215,8 @@ FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
                                          FaceType face_type_)
    : QuadratureSpaceBase(mesh_, mesh_.GetTypicalFaceGeometry(), ir),
      face_type(face_type_),
-     num_faces(mesh.GetNFbyType(face_type))
+     face_indices(mesh.GetFaceIndices(face_type_)),
+     face_indices_inv(mesh.GetInvFaceIndices(face_type_))
 {
    MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension() - 1) <= 1,
                "Constructor not valid for mixed meshes");
@@ -201,28 +225,30 @@ FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
 
 void FaceQuadratureSpace::ConstructOffsets()
 {
-   face_indices.SetSize(num_faces);
-   offsets.SetSize(num_faces + 1);
-   ne = num_faces;
-   int offset = 0;
-   int f_idx = 0;
-   for (int i = 0; i < mesh.GetNumFacesWithGhost(); i++)
-   {
-      const Mesh::FaceInformation face = mesh.GetFaceInformation(i);
-      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(face_type))
-      {
-         continue;
-      }
-      face_indices[f_idx] = i;
-      face_indices_inv[i] = f_idx;
-      offsets[f_idx] = offset;
-      Geometry::Type geom = mesh.GetFaceGeometry(i);
-      MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
-      offset += int_rule[geom]->GetNPoints();
+   ne = face_indices.Size();
 
-      f_idx++;
+   if (mesh.GetNumGeometries(mesh.Dimension() - 1) == 1)
+   {
+      Array<Geometry::Type> geoms;
+      mesh.GetGeometries(mesh.Dimension() - 1, geoms);
+      offsets.SetSize(1);
+      offsets.HostWrite();
+      offsets[0] = int_rule[geoms[0]]->GetNPoints();
+      size = ne * offsets[0];
+   }
+   else
+   {
+      offsets.SetSize(face_indices.Size() + 1);
+      int offset = 0;
+      for (int i = 0; i < mesh.GetNFbyType(face_type); ++i)
+      {
+         offsets[i] = offset;
+         Geometry::Type geom = mesh.GetFaceGeometry(face_indices[i]);
+         MFEM_ASSERT(int_rule[geom] != nullptr, "Missing integration rule");
+         offset += int_rule[geom]->GetNPoints();
+      }
+      offsets[face_indices.Size()] = size = offset;
    }
-   offsets[num_faces] = size = offset;
 }
 
 void FaceQuadratureSpace::Construct()
@@ -248,6 +274,16 @@ int FaceQuadratureSpace::GetPermutedIndex(int idx, int iq) const
    }
 }
 
+ElementTransformation *FaceQuadratureSpace::GetTransformation(int idx)
+{
+   ElementTransformation *T = mesh.GetFaceTransformation(face_indices[idx]);
+   if (face_type == FaceType::Boundary)
+   {
+      T->Attribute = mesh.GetBdrFaceAttributes()[idx];
+   }
+   return T;
+}
+
 int FaceQuadratureSpace::GetEntityIndex(const ElementTransformation &T) const
 {
    auto get_face_index = [this](const int idx)

fem/qspace.hpp

@@ -19,6 +19,12 @@
 namespace mfem
 {
 
+enum class QSpaceOffsetStorage
+{
+   FULL,
+   COMPRESSED
+};
+
 /// Abstract base class for QuadratureSpace and FaceQuadratureSpace.
 /** This class represents the storage layout for QuadratureFunction%s, that may
     be defined either on mesh elements or mesh faces. */
@@ -30,7 +36,7 @@ protected:
    Mesh &mesh; ///< The underlying mesh.
    int order; ///< The order of integration rule.
    int size; ///< Total number of quadrature points.
-   int ne; ///< Actual number of entities
+   int ne; ///< Number of entities
    mutable Vector weights; ///< Integration weights.
    mutable long nodes_sequence = 0; ///< Nodes counter for cache invalidation.
 
@@ -41,7 +47,17 @@ protected:
    /// can be computed as i * offsets[0], where i is the entity index. Otherwise
    /// has size num_entities + 1.
    ///
+   /// In the non-compressed case, the quadrature point values for entity i are
+   /// stored in the indices between offsets[i] and offsets[i+1].
    Array<int> offsets;
+
+   /// @brief Cached version of the "full" offsets, returned by Offsets() when
+   /// QSpaceOffsetStorage::FULL is provided.
+   ///
+   /// The quadrature point values for entity i are stored in the indices
+   /// between offsets[i] and offsets[i+1].
+   mutable Array<int> full_offset_cache;
+
    /// The quadrature rules used for each geometry type.
    const IntegrationRule *int_rule[Geometry::NumGeom];
 
@@ -74,12 +90,18 @@ public:
 
    /// @brief Entity quadrature point offset array.
    ///
-   /// Supports a constant compression scheme for meshes which have a single
+   /// If @a storage is QSpaceOffsetStorage::COMPRESSED, then the returned array
+   /// supports a constant compression scheme for meshes which have a single
    /// geometry type. When compressed, will have a single value. The true offset
    /// can be computed as i * offsets[0], where i is the entity index. Otherwise
    /// has size num_entities + 1.
    ///
-   const Array<int> &Offsets() const { return offsets; }
+   /// If @a storage is QSpaceOffsetStorage::FULL, then the array will never be
+   /// compressed.
+   ///
+   /// In the non-compressed case, the quadrature point values for entity i are
+   /// stored in the indices between offsets[i] and offsets[i+1].
+   const Array<int> &Offsets(QSpaceOffsetStorage storage) const;
 
    /// Return the total number of quadrature points.
    int GetSize() const { return size; }
@@ -191,13 +213,12 @@ public:
 class FaceQuadratureSpace : public QuadratureSpaceBase
 {
    FaceType face_type; ///< Is the space defined on interior or boundary faces?
-   const int num_faces; ///< Number of faces.
 
    /// Map from boundary or interior face indices to mesh face indices.
-   Array<int> face_indices;
+   const Array<int> &face_indices;
 
    /// Inverse of the map @a face_indices.
-   std::unordered_map<int,int> face_indices_inv;
+   const std::unordered_map<int,int> &face_indices_inv;
 
    const Vector &GetGeometricFactorWeights() const override;
    void ConstructOffsets();
@@ -213,14 +234,13 @@ public:
                        FaceType face_type_);
 
    /// Returns number of faces in the mesh.
-   inline int GetNumFaces() const { return num_faces; }
+   inline int GetNumFaces() const { return face_indices.Size(); }
 
    /// Returns the face type (boundary or interior).
    FaceType GetFaceType() const { return face_type; }
 
    /// Returns the face transformation of face @a idx.
-   ElementTransformation *GetTransformation(int idx) override
-   { return mesh.GetFaceTransformation(face_indices[idx]); }
+   ElementTransformation *GetTransformation(int idx) override;
 
    /// Returns the geometry type of face @a idx.
    Geometry::Type GetGeometry(int idx) const override

fem/tmop.cpp

@@ -140,6 +140,36 @@ void add_3D(const scalartype &scalar, const std::vector<type> &u,
 
 /* Metric definitions */
 
+// W = ||T||^2 - 2*det(T).
+template <typename type>
+type mu4_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   auto fnorm2 = fnorm2_2D(T);
+   auto det = det_2D(T);
+   return fnorm2 - 2*det;
+};
+
+// W = ||T-I||^2.
+template <typename type>
+type mu14_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   DenseMatrix Id(2,2); Id = 0.0;
+   Id(0,0) = 1; Id(1,1) = 1;
+
+   std::vector<type> Mat;
+   add_2D(real_t{-1.0}, T, &Id, Mat);
+
+   return fnorm2_2D(Mat);
+};
+
+// W = (det(T)-1)^2.
+template <typename type>
+type mu55_ad(const std::vector<type> &T, const std::vector<type> &W)
+{
+   auto det = det_2D(T);
+   return pow(det-1.0, 2.0);
+};
+
 // W = |T-T'|^2, where T'= |T|*I/sqrt(2).
 template <typename type>
 type mu85_ad(const std::vector<type> &T, const std::vector<type> &W)
@@ -163,6 +193,63 @@ type mu98_ad(const std::vector<type> &T, const std::vector<type> &W)
    return fnorm2_2D(Mat)/det_2D(T);
 };
 
+template <typename type>
+type make_one_type()
+{
+   return 1.0;
+}
+// add specialization for AD1Type
+template <>
+AD1Type make_one_type<AD1Type>()
+{
+   return AD1Type{1.0, 0.0};
+}
+// add specialization for AD2Type
+template <>
+AD2Type make_one_type<AD2Type>()
+{
+   return AD2Type{AD1Type{1.0, 0.0}, AD1Type{0.0, 0.0}};
+}
+
+using TWCUO = TMOP_WorstCaseUntangleOptimizer_Metric;
+template <typename type>
+type wcuo_ad(type mu,
+             const std::vector<type> &T, const std::vector<type> &W,
+             real_t alpha, real_t min_detT, real_t detT_ep,
+             int exponent, real_t max_muT, real_t muT_ep,
+             TWCUO::BarrierType bt,
+             TWCUO::WorstCaseType wct)
+{
+   type one = make_one_type<type>();
+   type zero = 0.0*one;
+   type denom = one;
+   if (bt == TWCUO::BarrierType::Shifted)
+   {
+      auto val1 = alpha*min_detT-detT_ep < 0.0 ?
+                  (alpha*min_detT-detT_ep)*one :
+                  zero;
+      denom = 2.0*(det_2D(T)-val1);
+   }
+   else if (bt == TWCUO::BarrierType::Pseudo)
+   {
+      auto detT = det_2D(T);
+      denom = detT + sqrt(detT*detT + detT_ep*detT_ep);
+   }
+   mu = mu/denom;
+
+   if (wct == TWCUO::WorstCaseType::PMean)
+   {
+      auto exp = exponent*one;
+      mu = pow(mu, exp);
+   }
+   else if (wct == TWCUO::WorstCaseType::Beta)
+   {
+      auto beta = (max_muT+muT_ep)*one;
+      mu = mu/(beta-mu);
+   }
+   return mu;
+}
+
 // W = 1/(tau^0.5) |T-I|^2.
 template <typename type>
 type mu342_ad(const std::vector<type> &T, const std::vector<type> &W)
@@ -421,7 +508,7 @@ void TMOP_QualityMetric::DefaultAssembleH(const DenseTensor &H,
          {
             for (int cc = 0; cc < dim; cc++)
             {
-               const double entry_rr_cc = Hrc(rr, cc);
+               const real_t entry_rr_cc = Hrc(rr, cc);
 
                for (int i = 0; i < dof; i++)
                {
@@ -481,6 +568,30 @@ void TMOP_Combo_QualityMetric::EvalPW(const DenseMatrix &Jpt,
    }
 }
 
+AD1Type TMOP_Combo_QualityMetric::EvalW_AD1(const std::vector<AD1Type> &T,
+                                            const std::vector<AD1Type> &W)
+const
+{
+   AD1Type metric = {0., 0.};
+   for (int i = 0; i < tmop_q_arr.Size(); i++)
+   {
+      metric += wt_arr[i]*tmop_q_arr[i]->EvalW_AD1(T, W);
+   }
+   return metric;
+}
+
+AD2Type TMOP_Combo_QualityMetric::EvalW_AD2(const std::vector<AD2Type> &T,
+                                            const std::vector<AD2Type> &W)
+const
+{
+   AD2Type metric = {{0., 0.},{0., 0.}};
+   for (int i = 0; i < tmop_q_arr.Size(); i++)
+   {
+      metric += wt_arr[i]*tmop_q_arr[i]->EvalW_AD2(T, W);
+   }
+   return metric;
+}
+
 void TMOP_Combo_QualityMetric::AssembleH(const DenseMatrix &Jpt,
                                          const DenseMatrix &DS,
                                          const real_t weight,
@@ -645,6 +756,64 @@ real_t TMOP_WorstCaseUntangleOptimizer_Metric::EvalWBarrier(
    return tmop_metric.EvalW(Jpt)/denominator;
 }
 
+AD1Type TMOP_WorstCaseUntangleOptimizer_Metric::EvalW_AD1(
+   const std::vector<AD1Type> &T,
+   const std::vector<AD1Type> &W) const
+{
+   return wcuo_ad(tmop_metric.EvalW_AD1(T,W), T, W, alpha, min_detT, detT_ep,
+                  exponent, max_muT, muT_ep, btype, wctype);
+}
+
+AD2Type TMOP_WorstCaseUntangleOptimizer_Metric::EvalW_AD2(
+   const std::vector<AD2Type> &T,
+   const std::vector<AD2Type> &W) const
+{
+   return wcuo_ad(tmop_metric.EvalW_AD2(T,W), T, W, alpha, min_detT, detT_ep,
+                  exponent, max_muT, muT_ep, btype, wctype);
+}
+
+void TMOP_WorstCaseUntangleOptimizer_Metric::EvalP(const DenseMatrix &Jpt,
+                                                   DenseMatrix &P) const
+{
+   auto mu_ad_fn = [this](std::vector<AD1Type> &T, std::vector<AD1Type> &W)
+   {
+      return EvalW_AD1(T,W);
+   };
+   if (tmop_metric.Id() == 4 || tmop_metric.Id() == 14 ||
+       tmop_metric.Id() == 66)
+   {
+      ADGrad(mu_ad_fn, P, Jpt);
+      return;
+   }
+   MFEM_ABORT("EvalW_AD1 not implemented with this metric for "
+              "TMOP_WorstCaseUntangleOptimizer_Metric. "
+              "Please use metric 4/14/66.");
+}
+
+void TMOP_WorstCaseUntangleOptimizer_Metric::AssembleH(
+   const DenseMatrix &Jpt,
+   const DenseMatrix &DS,
+   const real_t weight,
+   DenseMatrix &A) const
+{
+   DenseTensor H(Jpt.Height(), Jpt.Height(), Jpt.TotalSize());
+   H = 0.0;
+   auto mu_ad_fn = [this](std::vector<AD2Type> &T, std::vector<AD2Type> &W)
+   {
+      return EvalW_AD2(T,W);
+   };
+   if (tmop_metric.Id() == 4 || tmop_metric.Id() == 14 ||
+       tmop_metric.Id() == 66)
+   {
+      ADHessian(mu_ad_fn, H, Jpt);
+      this->DefaultAssembleH(H,DS,weight,A);
+      return;
+   }
+   MFEM_ABORT("EvalW_AD1 not implemented with this metric for "
+              "TMOP_WorstCaseUntangleOptimizer_Metric. "
+              "Please use metric 4/14/66.");
+}
+
 real_t TMOP_Metric_001::EvalW(const DenseMatrix &Jpt) const
 {
    ie.SetJacobian(Jpt.GetData());
@@ -850,6 +1019,25 @@ void TMOP_Metric_004::AssembleH(const DenseMatrix &Jpt,
    ie.Assemble_ddI2b(-2.0*weight, A.GetData());
 }
 
+template <typename type>
+type TMOP_Metric_004::EvalW_AD_impl(const std::vector<type> &T,
+                                    const std::vector<type> &W) const
+{
+   return mu4_ad(T, W);
+}
+
+AD1Type TMOP_Metric_004::EvalW_AD1(const std::vector<AD1Type> &T,
+                                   const std::vector<AD1Type> &W) const
+{
+   return EvalW_AD_impl<AD1Type>(T,W);
+}
+
+AD2Type TMOP_Metric_004::EvalW_AD2(const std::vector<AD2Type> &T,
+                                   const std::vector<AD2Type> &W) const
+{
+   return EvalW_AD_impl<AD2Type>(T,W);
+}
+
 real_t TMOP_Metric_007::EvalW(const DenseMatrix &Jpt) const
 {
    // mu_7 = |J-J^{-t}|^2 = |J|^2 + |J^{-1}|^2 - 4
@@ -971,6 +1159,25 @@ void TMOP_Metric_014::AssembleH(const DenseMatrix &Jpt,
    ie.Assemble_ddI1(weight, A.GetData());
 }
 
+template <typename type>
+type TMOP_Metric_014::EvalW_AD_impl(const std::vector<type> &T,
+                                    const std::vector<type> &W) const
+{
+   return mu14_ad(T, W);
+}
+
+AD1Type TMOP_Metric_014::EvalW_AD1(const std::vector<AD1Type> &T,
+                                   const std::vector<AD1Type> &W) const
+{
+   return EvalW_AD_impl<AD1Type>(T,W);
+}
+
+AD2Type TMOP_Metric_014::EvalW_AD2(const std::vector<AD2Type> &T,
+                                   const std::vector<AD2Type> &W) const
+{
+   return EvalW_AD_impl<AD2Type>(T,W);
+}
+
 real_t TMOP_Metric_022::EvalW(const DenseMatrix &Jpt) const
 {
    // mu_22 = (0.5*|J|^2 - det(J)) / (det(J) - tau0)
@@ -1100,6 +1307,25 @@ void TMOP_Metric_055::AssembleH(const DenseMatrix &Jpt,
    ie.Assemble_ddI2b(2*weight*(ie.Get_I2b() - 1.0), A.GetData());
 }
 
+template <typename type>
+type TMOP_Metric_055::EvalW_AD_impl(const std::vector<type> &T,
+                                    const std::vector<type> &W) const
+{
+   return mu55_ad(T, W);
+}
+
+AD1Type TMOP_Metric_055::EvalW_AD1(const std::vector<AD1Type> &T,
+                                   const std::vector<AD1Type> &W) const
+{
+   return EvalW_AD_impl<AD1Type>(T,W);
+}
+
+AD2Type TMOP_Metric_055::EvalW_AD2(const std::vector<AD2Type> &T,
+                                   const std::vector<AD2Type> &W) const
+{
+   return EvalW_AD_impl<AD2Type>(T,W);
+}
+
 real_t TMOP_Metric_056::EvalWMatrixForm(const DenseMatrix &Jpt) const
 {
    // mu_56 = 0.5 (det(J) + 1 / det(J)) - 1.
@@ -4096,7 +4322,7 @@ real_t TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
 
          const IntegrationPoint &ip_s = ir_s->IntPoint(s);
          Tpr->SetIntPoint(&ip_s);
-         double w = surf_fit_coeff->Eval(*Tpr, ip_s) * surf_fit_normal *
+         real_t w = surf_fit_coeff->Eval(*Tpr, ip_s) * surf_fit_normal *
                     1.0 / surf_fit_dof_count[scalar_dof_id];
 
          if (surf_fit_gf)

fem/tmop.hpp

@@ -14,10 +14,14 @@
 
 #include "../linalg/invariants.hpp"
 #include "nonlininteg.hpp"
+#include "../linalg/dual.hpp"
 
 namespace mfem
 {
 
+using AD1Type = future::dual<real_t, real_t>;
+using AD2Type = future::dual<AD1Type, AD1Type>;
+
 /** @brief Abstract class for local mesh quality metrics in the target-matrix
     optimization paradigm (TMOP) by P. Knupp et al. */
 class TMOP_QualityMetric : public HyperelasticModel
@@ -69,6 +73,22 @@ public:
    virtual void EvalPW(const DenseMatrix &Jpt, DenseMatrix &PW) const
    { PW = 0.0;}
 
+   /// @brief First-derivative hook for AD-based computations.
+   /// @warning Not for public use. Internal use for AD-based computations.
+   virtual AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                             const std::vector<AD1Type> &W) const
+   {
+      MFEM_ABORT("EvalW_AD1 not implemented for this metric");
+   }
+
+   /// @brief Second-derivative hook for AD-based computations.
+   /// @warning Not for public use. Internal use for AD-based computations.
+   virtual AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                             const std::vector<AD2Type> &W) const
+   {
+      MFEM_ABORT("EvalW_AD2 not implemented for this metric");
+   }
+
    /** @brief Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor
        and assemble its contribution to the local gradient matrix 'A'.
        @param[in] Jpt     Represents the target->physical transformation
@@ -124,6 +144,12 @@ public:
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
                   const real_t weight, DenseMatrix &A) const override;
 
+   AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                     const std::vector<AD1Type> &W) const override;
+
+   AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                     const std::vector<AD2Type> &W) const override;
+
    /// Computes the averages of all metrics (integral of metric / volume).
    /// Works in parallel when called with a ParGridFunction.
    void ComputeAvgMetrics(const GridFunction &nodes,
@@ -221,12 +247,16 @@ public:
 
    real_t EvalW(const DenseMatrix &Jpt) const override;
 
-   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override
-   { MFEM_ABORT("Not implemented"); }
+   AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                     const std::vector<AD1Type> &W) const override;
+
+   AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                     const std::vector<AD2Type> &W) const override;
+
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
 
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                  const real_t weight, DenseMatrix &A) const override
-   { MFEM_ABORT("Not implemented"); }
+                  const real_t weight, DenseMatrix &A) const override;
 
    // Compute mu_hat.
    real_t EvalWBarrier(const DenseMatrix &Jpt) const;
@@ -368,6 +398,10 @@ class TMOP_Metric_004 : public TMOP_QualityMetric
 protected:
    mutable InvariantsEvaluator2D<real_t> ie;
 
+   template<typename type>
+   type EvalW_AD_impl(const std::vector<type> &T,
+                      const std::vector<type> &W) const;
+
 public:
    // W = |J|^2 - 2*det(J)
    real_t EvalW(const DenseMatrix &Jpt) const override;
@@ -377,6 +411,12 @@ public:
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
                   const real_t weight, DenseMatrix &A) const override;
 
+   AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                     const std::vector<AD1Type> &W) const override;
+
+   AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                     const std::vector<AD2Type> &W) const override;
+
    int Id() const override { return 4; }
 };
 
@@ -420,6 +460,10 @@ class TMOP_Metric_014 : public TMOP_QualityMetric
 protected:
    mutable InvariantsEvaluator2D<real_t> ie;
 
+   template <typename type>
+   type EvalW_AD_impl(const std::vector<type> &T,
+                      const std::vector<type> &W) const;
+
 public:
    // W = |J - I|^2.
    real_t EvalWMatrixForm(const DenseMatrix &Jpt) const override;
@@ -431,6 +475,14 @@ public:
 
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
                   const real_t weight, DenseMatrix &A) const override;
+
+   AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                     const std::vector<AD1Type> &W) const override;
+
+   AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                     const std::vector<AD2Type> &W) const override;
+
+   int Id() const override { return 14; }
 };
 
 /// 2D Shifted barrier form of shape metric (mu_2).
@@ -479,6 +531,10 @@ class TMOP_Metric_055 : public TMOP_QualityMetric
 protected:
    mutable InvariantsEvaluator2D<real_t> ie;
 
+   template<typename type>
+   type EvalW_AD_impl(const std::vector<type> &T,
+                      const std::vector<type> &W) const;
+
 public:
    // W = (det(J) - 1)^2.
    real_t EvalW(const DenseMatrix &Jpt) const override;
@@ -488,6 +544,14 @@ public:
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
                   const real_t weight, DenseMatrix &A) const override;
 
+   AD1Type EvalW_AD1(const std::vector<AD1Type> &T,
+                     const std::vector<AD1Type> &W) const override;
+
+   AD2Type EvalW_AD2(const std::vector<AD2Type> &T,
+                     const std::vector<AD2Type> &W) const override;
+
+   int Id() const override { return 55; }
+
 };
 
 /// 2D barrier size (V) metric (polyconvex).

fem/tmop_tools.hpp

@@ -111,7 +111,7 @@ public:
 
    ~InterpolatorFP()
    {
-      finder->FreeData();
+      if (finder) { finder->FreeData(); }
       delete finder;
    }
 };

general/hash.hpp

@@ -57,9 +57,9 @@ struct Hashed4
  *  each time this class is invoked.
  *
  *  There are two main methods this class provides. The Get(...) methods always
- *  return an item given the two or four indices. If the item didn't previously
+ *  return an item given the two or four indices. If the item did not previously
  *  exist, the methods creates a new one. The Find(...) methods, on the other
- *  hand, just return NULL or -1 if the item doesn't exist.
+ *  hand, just return NULL or -1 if the item does not exist.
  *
  *  Each new item is automatically assigned a unique ID - the index of the item
  *  inside the BlockArray. The IDs may (but need not) be used as p1, p2, ... of
@@ -95,14 +95,14 @@ public:
        @param[in] init_hash_size The initial size of the hash table. Must be
                                  a power of 2. */
    HashTable(int block_size = 16*1024, int init_hash_size = 32*1024);
-   /// @brief Deep copy
+   /// Deep copy
    HashTable(const HashTable& other);
-   /// @brief Copy assignment not supported
+   /// Copy assignment not supported
    HashTable& operator=(const HashTable&) = delete;
    ~HashTable();
 
-   /** @brief Item accessor with key (or parents) the pair 'p1', 'p2'. Default
-       construct an item of type T if no value correspond to the requested key.
+   /** @brief Item accessor with key (or parents) the pair p1, p2. Default
+       construct an item of type T if no value corresponds to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -111,9 +111,9 @@ public:
        @warning This method should only be called if T inherits from Hashed2. */
    T* Get(int p1, int p2);
 
-   /** @brief Item accessor with key (or parents) the quadruplet 'p1', 'p2',
-       'p3', 'p4'. The key 'p4' is optional. Default construct an item of type T
-       if no value corresponds to the requested key.
+   /** @brief Item accessor with key (or parents) the quadruplet p1, p2, p3, p4.
+       The key p4 is optional. Default construct an item of type T if no value
+       corresponds to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -124,10 +124,10 @@ public:
       @warning This method should only be called if T inherits from Hashed4. */
    T* Get(int p1, int p2, int p3, int p4 = -1 /* p4 optional */);
 
-   /// Get id of item whose parents are p1, p2... Create it if it doesn't exist.
-   /** @brief Get the "id" of an item, this "id" corresponding to the index of the
-       item in the underlying BlockArray<T> object. Default construct an item
-       and id if no value corresponds to the requested key.
+   /** @brief Get the "id" of the item whose parents are p1, p2, this "id"
+       corresponding to the index of the item in the underlying BlockArray<T>
+       object. Default construct an item and "id" if no value corresponds to the
+       requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -136,9 +136,9 @@ public:
        @warning This method should only be called if T inherits from Hashed2. */
    int GetId(int p1, int p2);
 
-   /** @brief Get the "id" of an item, this "id" corresponding to the index of the
-       item in the underlying BlockArray<T> object. Default construct an item
-       and id if no value correspond to the requested key.
+   /** @brief Get the "id" of an item, this "id" corresponding to the index of
+       the item in the underlying BlockArray<T> object. Default construct an item
+       and "id" if no value corresponds to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -149,9 +149,8 @@ public:
        @warning This method should only be called if T inherits from Hashed4. */
    int GetId(int p1, int p2, int p3, int p4 = -1);
 
-   /// Find item whose parents are p1, p2... Return NULL if it doesn't exist.
-   /** @brief Item accessor with key (or parents) the pair 'p1', 'p2'. Return
-       nullptr if no value correspond to the requested key.
+   /** @brief Item accessor with key (or parents) the pair p1, p2. Return
+       NULL if no value corresponds to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -160,9 +159,9 @@ public:
        @warning This method should only be called if T inherits from Hashed2. */
    T* Find(int p1, int p2);
 
-   /** @brief Item accessor with key (or parents) the quadruplet 'p1', 'p2',
-       'p3', 'p4'. The key 'p4' is optional. Return nullptr if no value
-       correspond to the requested key.
+   /** @brief Item accessor with key (or parents) the quadruplet p1, p2, p3, p4.
+       The key p4 is optional. Return NULL if no value corresponds to the
+       requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -173,8 +172,8 @@ public:
        @warning This method should only be called if T inherits from Hashed4. */
    T* Find(int p1, int p2, int p3, int p4 = -1);
 
-   /** @brief Item const accessor with key (or parents) the pair 'p1', 'p2'.
-       Return nullptr if no value correspond to the requested key.
+   /** @brief Item const accessor with key (or parents) the pair p1, p2.
+       Return NULL if no value corresponds to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -183,9 +182,9 @@ public:
       @warning This method should only be called if T inherits from Hashed2. */
    const T* Find(int p1, int p2) const;
 
-   /** @brief Item const accessor with key (or parents) the quadruplet 'p1',
-       'p2', 'p3', 'p4'. The key 'p4' is optional. Return nullptr if no value
-       correspond to the requested key.
+   /** @brief Item const accessor with key (or parents) the quadruplet p1, p2,
+       p3, p4. The key p4 is optional. Return NULL if no value corresponds to the
+       requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -196,10 +195,12 @@ public:
        @warning This method should only be called if T inherits from Hashed4. */
    const T* Find(int p1, int p2, int p3, int p4 = -1) const;
 
-   /// Find id of item whose parents are p1, p2... Return -1 if it doesn't exist.
-   /** @brief Find the "id" of an item, this "id" corresponding to the index of
-       the item in the underlying BlockArray<T> object. Default construct an
-       item and id if no value correspond to the requested key.
+   /** @brief Find the "id" of an item whose parents are p1, p2. Return -1 if it
+       does not exist.
+
+       This "id" corresponds to the index of the item in the underlying
+       BlockArray<T> object. Default construct an item and "id" if no value
+       corresponds to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -209,8 +210,9 @@ public:
    int FindId(int p1, int p2) const;
 
    /** @brief Find the "id" of an item, this "id" corresponding to the index of
-       the item in the underlying BlockArray<T> object. Default construct an
-       item and id if no value correspond to the requested key.
+       the item in the underlying BlockArray<T> object. Return -1 if it does not
+       exist. Default construct an item and "id" if no value corresponds to the
+       requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -221,16 +223,16 @@ public:
        @warning This method should only be called if T inherits from Hashed4. */
    int FindId(int p1, int p2, int p3, int p4 = -1) const;
 
-   /// @brief Return the number of elements currently stored in the HashTable.
+   /// Return the number of elements currently stored in the HashTable.
    int Size() const { return Base::Size() - unused.Size(); }
 
-   /// @brief Return the total number of ids (used and unused) in the HashTable.
+   /// Return the total number of ids (used and unused) in the HashTable.
    int NumIds() const { return Base::Size(); }
 
-   /// @brief Return the number of free/unused ids in the HashTable.
+   /// Return the number of free/unused ids in the HashTable.
    int NumFreeIds() const { return unused.Size(); }
 
-   /** @brief Return true if item 'id' exists in (is used by) the container.
+   /** @brief Return true if item @a id exists in (is used by) the container.
 
        @param[in] id Index of the item in the underlying BlockArray<T>.
 
@@ -241,13 +243,13 @@ public:
 
        @param[in] id Index of the item in the underlying BlockArray<T>.
 
-       @warning Its id will be reused by newly added items. */
+       @warning Its @a id will be reused by newly added items. */
    void Delete(int id);
 
-   /// @brief Remove all items.
+   /// Remove all items.
    void DeleteAll();
 
-   /** @brief Allocate an item at 'id'. Enlarge the underlying BlockArray if
+   /** @brief Allocate an item at @a id. Enlarge the underlying BlockArray if
        necessary.
 
        @param[in] id Index of the item in the underlying BlockArray<T>.
@@ -255,7 +257,7 @@ public:
        @param[in] p2 Second part of the key.
 
        @warning This is a special purpose method used when loading data from a
-       file. Does nothing if the slot 'id' has already been allocated. */
+       file. Does nothing if the slot @a id has already been allocated. */
    void Alloc(int id, int p1, int p2);
 
    /** @brief Reinitialize the internal list of unallocated items.
@@ -287,13 +289,13 @@ public:
        @warning This method should only be called if T inherits from Hashed4. */
    void Reparent(int id, int new_p1, int new_p2, int new_p3, int new_p4 = -1);
 
-   /// @brief Return total size of allocated memory (tables plus items), in bytes.
+   /// Return total size of allocated memory (tables plus items), in bytes.
    std::size_t MemoryUsage() const;
 
-   /// @brief Write details of the memory usage to the mfem output stream.
+   /// Write details of the memory usage to the mfem output stream.
    void PrintMemoryDetail() const;
 
-   /// @brief Print a histogram of bin sizes for debugging purposes.
+   /// Print a histogram of bin sizes for debugging purposes.
    void PrintStats() const;
 
    class iterator : public Base::iterator
@@ -346,7 +348,7 @@ public:
 
 protected:
    /** The hash table: each bin is a linked list of items. For each non-empty
-       bin, this arrays stores the 'id' of the first item in the list, or -1
+       bin, this arrays stores the "id" of the first item in the list, or -1
        if the bin is empty. */
    int* table;
 
@@ -384,11 +386,11 @@ protected:
    { return (984120265ul*p1 + 125965121ul*p2 + 495698413ul*p3) & mask; }
 
    // Delete() and Reparent() use one of these:
-   /// @brief Hash function for items of type T that inherit from Hashed2.
+   /// Hash function for items of type T that inherit from Hashed2.
    inline int Hash(const Hashed2& item) const
    { return Hash(item.p1, item.p2); }
 
-   /// @brief Hash function for items of type T that inherit from Hashed4.
+   /// Hash function for items of type T that inherit from Hashed4.
    inline int Hash(const Hashed4& item) const
    { return Hash(item.p1, item.p2, item.p3); }
 
@@ -415,15 +417,15 @@ protected:
        @warning This method should only be called if T inherits from Hashed4. */
    int SearchList(int id, int p1, int p2, int p3) const;
 
-   /** @brief Insert the item 'id' into bin 'idx'.
+   /** @brief Insert the item @a id into bin @a idx.
 
        @param[in] idx The bin/bucket index.
        @param[in] id The index of the item in the BlockArray<T>.
        @param[in] item The item to insert at the beginning of the linked list.
 
-       @warning The method only works with bin 'idx' and does not check the
-                overall fill factor of the hash table. If appropriate,
-                use CheckRehash() for that. */
+       @warning The method only works with bin @a idx and does not check the
+                overall fill factor of the hash table. If appropriate, use
+                CheckRehash() for that. */
    inline void Insert(int idx, int id, T &item);
 
    /** @brief Unlink an item @a id from the linked list of bin @a idx.
@@ -444,11 +446,11 @@ protected:
        and reinsert all items into the new bins.
 
        NOTE: Rehashing is computationally expensive (O(N) in the number of items),
-       but since it is only done rarely (when the number of items doubles),
-       the amortized complexity of inserting an item is still O(1). */
+       but since it is only done rarely (when the number of items doubles), the
+       amortized complexity of inserting an item is still O(1). */
    void DoRehash();
 
-   /** @brief Return the size of the bin "idx".
+   /** @brief Return the size of the bin @a idx.
 
        @param[in] idx The index of the bin.
        @return The size of the bin. */

general/mem_manager.cpp

@@ -1384,8 +1384,11 @@ void MemoryManager::Insert(void *h_ptr, size_t bytes,
    {
       auto &m = res.first->second;
       MFEM_VERIFY(m.bytes >= bytes && m.h_mt == h_mt &&
-                  (m.d_mt == d_mt || (d_mt == MemoryType::DEFAULT &&
-                                      m.d_mt == GetDualMemoryType(h_mt))),
+                  (m.d_mt == d_mt ||
+                   (d_mt == MemoryType::DEFAULT &&
+                    m.d_mt == GetDualMemoryType(h_mt)) ||
+                   (m.d_mt == MemoryType::DEFAULT &&
+                    d_mt == GetDualMemoryType(m.h_mt))),
                   "Address already present with different attributes!");
 #ifdef MFEM_TRACK_MEM_MANAGER
       mfem::out << "[mfem memory manager]: repeated registration of h_ptr: "

general/scan.hpp

@@ -0,0 +1,228 @@
+// Copyright (c) 2010-2025, 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_SCAN_HPP
+#define MFEM_SCAN_HPP
+
+#include "backends.hpp"
+
+#ifdef MFEM_USE_CUDA
+#include <cub/device/device_scan.cuh>
+#define MFEM_CUB_NAMESPACE cub
+#elif defined(MFEM_USE_HIP)
+#include <hipcub/device/device_scan.hpp>
+#define MFEM_CUB_NAMESPACE hipcub
+#endif
+
+#include <functional>
+#include <numeric>
+#include <cstddef>
+
+namespace mfem
+{
+/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
+/// std::plus<>{})
+template <class InputIt, class OutputIt>
+void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
+{
+   // forward to InclusiveSum for potentially faster kernels
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
+   {
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
+      {
+         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
+                       workspace.Write(), bytes, d_in, d_out, num_items);
+#if defined(MFEM_USE_CUDA)
+         if (err == cudaSuccess)
+         {
+            return;
+         }
+#elif defined(MFEM_USE_HIP)
+         if (err == hipSuccess)
+         {
+            return;
+         }
+#endif
+      }
+      // try allocating a larger buffer
+      bytes = 0;
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
+                        nullptr, bytes, d_in, d_out, num_items));
+      workspace.SetSize(bytes);
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
+                        workspace.Write(), bytes, d_in, d_out, num_items));
+      return;
+   }
+#endif
+#if 0
+   std::inclusive_scan(d_in, d_in + num_items, d_out);
+#else
+   // work-around to some compilers not fully supporting C++17
+   if (num_items)
+   {
+      *d_out = *d_in;
+      auto prev = d_out;
+      ++d_in;
+      ++d_out;
+      for (size_t i = 1; i < num_items; ++i)
+      {
+         *d_out = (*prev) + (*d_in);
+         prev = d_out;
+         ++d_in;
+         ++d_out;
+      }
+   }
+#endif
+}
+
+/// @brief Performs an inclusive scan of [d_in, d_in+num_items) -> [d_out,
+/// d_out+num_items). This call is potentially asynchronous on the device.
+///
+/// @a d_in input start.
+/// @a d_out output start. Can perform in-place scans with d_out = d_in
+/// @a scan_op binary scan functor. Must be associative. If only weakly
+/// associative (i.e. floating point addition) results are not deterministic. On
+/// device this must also be commutative.
+template <class InputIt, class OutputIt, class ScanOp>
+void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
+                   ScanOp scan_op)
+{
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
+   {
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
+      {
+         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
+                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
+#if defined(MFEM_USE_CUDA)
+         if (err == cudaSuccess)
+         {
+            return;
+         }
+#elif defined(MFEM_USE_HIP)
+         if (err == hipSuccess)
+         {
+            return;
+         }
+#endif
+      }
+      // try allocating a larger buffer
+      bytes = 0;
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
+                        nullptr, bytes, d_in, d_out, scan_op, num_items));
+      workspace.SetSize(bytes);
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
+                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
+      return;
+   }
+#endif
+#if 0
+   std::inclusive_scan(d_in, d_in + num_items, d_out, scan_op);
+#else
+   // work-around to some compilers not fully supporting C++17
+   if (num_items)
+   {
+      *d_out = *d_in;
+      auto prev = d_out;
+      ++d_in;
+      ++d_out;
+      for (size_t i = 1; i < num_items; ++i)
+      {
+         *d_out = scan_op(*prev, *d_in);
+         prev = d_out;
+         ++d_in;
+         ++d_out;
+      }
+   }
+#endif
+}
+
+/// Performs an exclusive scan of [d_in, d_in+num_items) -> [d_out,
+/// d_out+num_items). This call is potentially asynchronous on the device.
+/// @a d_in input start.
+/// @a d_out output start. Can perform in-place scans with d_out = d_in
+/// @a scan_op binary scan functor. Must be associative. If only weakly
+/// associative (i.e. floating point addition) results are not deterministic. On
+/// device this must also be commutative.
+template <class InputIt, class OutputIt, class T, class ScanOp>
+void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
+                   T init_value, ScanOp scan_op)
+{
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
+   {
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
+      {
+         auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
+                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                       num_items);
+#if defined(MFEM_USE_CUDA)
+         if (err == cudaSuccess)
+         {
+            return;
+         }
+#elif defined(MFEM_USE_HIP)
+         if (err == hipSuccess)
+         {
+            return;
+         }
+#endif
+      }
+      // try allocating a larger buffer
+      bytes = 0;
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
+                        nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
+      workspace.SetSize(bytes);
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
+                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                        num_items));
+      return;
+   }
+#endif
+#if 0
+   std::exclusive_scan(d_in, d_in + num_items, d_out, init_value, scan_op);
+#else
+   // work-around to some compilers not fully supporting C++17
+   if (num_items)
+   {
+      for (size_t i = 0; i < num_items; ++i)
+      {
+         auto next = scan_op(init_value, *d_in);
+         *d_out = init_value;
+         init_value = next;
+         ++d_out;
+         ++d_in;
+      }
+   }
+#endif
+}
+
+/// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
+/// workspace, std::plus<>{})
+template <class InputIt, class OutputIt, class T>
+void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
+                   T init_value)
+{
+   ExclusiveScan(use_dev, d_in, d_out, num_items, init_value, std::plus<> {});
+}
+
+} // namespace mfem
+
+#undef MFEM_CUB_NAMESPACE
+
+#endif

mesh/mesh.cpp

@@ -979,6 +979,46 @@ const Array<int>& Mesh::GetElementAttributes() const
    return elem_attrs_cache;
 }
 
+void Mesh::ComputeFaceInfo(FaceType ftype) const
+{
+   auto &fidcs = face_indices[static_cast<int>(ftype)];
+   auto &ifidcs = inv_face_indices[static_cast<int>(ftype)];
+   fidcs.SetSize(GetNFbyType(ftype));
+   fidcs.HostWrite();
+   ifidcs.reserve(fidcs.Size());
+   int f_idx = 0;
+   for (int i = 0; i < GetNumFacesWithGhost(); ++i)
+   {
+      const FaceInformation face = GetFaceInformation(i);
+      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(ftype))
+      {
+         continue;
+      }
+      fidcs[f_idx] = i;
+      ifidcs[i] = f_idx;
+      ++f_idx;
+   }
+}
+
+const Array<int> &Mesh::GetFaceIndices(FaceType ftype) const
+{
+   if (face_indices[static_cast<int>(ftype)].Size() == 0)
+   {
+      ComputeFaceInfo(ftype);
+   }
+   return face_indices[static_cast<int>(ftype)];
+}
+
+const std::unordered_map<int, int> &
+Mesh::GetInvFaceIndices(FaceType ftype) const
+{
+   if (inv_face_indices[static_cast<int>(ftype)].empty())
+   {
+      ComputeFaceInfo(ftype);
+   }
+   return inv_face_indices[static_cast<int>(ftype)];
+}
+
 void Mesh::DeleteGeometricFactors()
 {
    for (int i = 0; i < geom_factors.Size(); i++)
@@ -1866,6 +1906,13 @@ void Mesh::Destroy()
    bdr_face_attrs_cache.DeleteAll();
    attributes.DeleteAll();
    bdr_attributes.DeleteAll();
+
+   face_indices[0].DeleteAll();
+   face_indices[1].DeleteAll();
+   // force de-allocation so after this mesh has the smallest memory footprint
+   // possible
+   inv_face_indices[0] = std::unordered_map<int, int>();
+   inv_face_indices[1] = std::unordered_map<int, int>();
 }
 
 void Mesh::ResetLazyData()
@@ -1880,6 +1927,11 @@ void Mesh::ResetLazyData()
    // set size to 0 so re-computations can potentially avoid a new allocation
    bdr_face_attrs_cache.SetSize(0);
    elem_attrs_cache.SetSize(0);
+
+   face_indices[0].SetSize(0);
+   face_indices[1].SetSize(0);
+   inv_face_indices[0].clear();
+   inv_face_indices[1].clear();
 }
 
 void Mesh::SetAttributes(bool elem_attrs_changed, bool bdr_face_attrs_changed)
@@ -8137,6 +8189,12 @@ void Mesh::GenerateFaces()
       FreeElement(f);
    }
 
+   // delete caches
+   face_indices[0].SetSize(0);
+   face_indices[1].SetSize(0);
+   inv_face_indices[0].clear();
+   inv_face_indices[1].clear();
+
    // (re)generate the interior faces and the info for them
    faces.SetSize(nfaces);
    faces_info.SetSize(nfaces);
@@ -10936,6 +10994,11 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
    // copy attribute caches
    mfem::Swap(elem_attrs_cache, other.elem_attrs_cache);
    mfem::Swap(bdr_face_attrs_cache, other.bdr_face_attrs_cache);
+
+   mfem::Swap(face_indices[0], other.face_indices[0]);
+   mfem::Swap(face_indices[1], other.face_indices[1]);
+   inv_face_indices[0].swap(other.inv_face_indices[0]);
+   inv_face_indices[1].swap(other.inv_face_indices[1]);
 }
 
 void Mesh::GetElementData(const Array<Element*> &elem_array, int geom,

mesh/mesh.hpp

@@ -278,6 +278,13 @@ protected:
 
    // used during NC mesh initialization only
    Array<Triple<int, int, int> > tmp_vertex_parents;
+   /// cache for FaceIndices(ftype)
+   mutable Array<int> face_indices[2];
+   /// cache for FaceIndices(ftype)
+   mutable std::unordered_map<int, int> inv_face_indices[2];
+
+   /// compute face_indices[ftype] and inv_face_indices[type]
+   void ComputeFaceInfo(FaceType ftype) const;
 
 public:
    typedef Geometry::Constants<Geometry::SEGMENT>     seg_t;
@@ -312,6 +319,11 @@ public:
    // (true) is set in mesh_readers.cpp.
    static bool remove_unused_vertices;
 
+   /// Map from boundary or interior face indices to mesh face indices.
+   const Array<int>& GetFaceIndices(FaceType ftype) const;
+   /// Inverse of the map FaceIndices(ftype)
+   const std::unordered_map<int, int>& GetInvFaceIndices(FaceType ftype) const;
+
 protected:
    Operation last_operation;
 

mesh/ncmesh.cpp

@@ -827,7 +827,9 @@ void NCMesh::ForceRefinement(int vn1, int vn2, int vn3, int vn4)
    Face* face = faces.Find(vn1, vn2, vn3, vn4);
    if (!face) { return; }
 
-   int elem = face->GetSingleElement();
+   MFEM_VERIFY(!IsParallel(), "ForceRefinement is supported only in serial");
+
+   const int elem = face->GetSingleElement();
    Element &el = elements[elem];
    MFEM_ASSERT(!el.ref_type, "element already refined.");
 

mesh/ncmesh.hpp

@@ -483,6 +483,9 @@ public:
 
    int PrintMemoryDetail() const;
 
+   /// Return true for ParNCMesh with more than one MPI process.
+   virtual bool IsParallel() const { return false; }
+
    using RefCoord = std::int64_t;
 
    static constexpr int MaxElemNodes =

mesh/pmesh.cpp

@@ -3897,6 +3897,13 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
 #endif
 }
 
+bool ParMesh::AnisotropicConflict(const Array<Refinement> &refinements,
+                                  std::set<int> &conflicts) const
+{
+   MFEM_VERIFY(pncmesh, "AnisotropicConflict should be called only for NCMesh");
+   return pncmesh->AnisotropicConflict(refinements, conflicts);
+}
+
 void ParMesh::NonconformingRefinement(const Array<Refinement> &refinements,
                                       int nc_limit)
 {

mesh/pmesh.hpp

@@ -815,6 +815,22 @@ public:
    /// Debugging method
    void PrintSharedEntities(const std::string &fname_prefix) const;
 
+   /** @brief Return true if the input array of refinements to be performed would
+       result in conflicting anisotropic directions on a face. Indices of
+       @a refinements entries are contained in @a conflicts, for marked elements
+       neighboring a face with a conflict.
+
+       The return value is globally MPI-reduced (true if any MPI process has a
+       conflict), whereas @a conflicts contains local indices of conflicting
+       entries of @a refinements. Conflicts are defined as anisotropic
+       refinements in different directions on a face shared by two elements.
+       Conflicts are checked for the mesh that would result from the input
+       refinements. If there are no conflicts, then the refinements can be
+       performed without forced refinements. This function is supported only for
+       3D meshes with all hexahedral elements. */
+   bool AnisotropicConflict(const Array<Refinement> &refinements,
+                            std::set<int> &conflicts) const;
+
    virtual ~ParMesh();
 };
 

mesh/pncmesh.cpp

@@ -1504,6 +1504,507 @@ void ParNCMesh::Prune()
    Update();
 }
 
+bool ParNCMesh::AnisotropicConflict(const Array<Refinement> &refinements,
+                                    std::set<int> &conflicts)
+{
+   if (Dim < 3 || NRanks == 1) { return false; }
+
+   for (int i = 0; i < refinements.Size() && Iso; i++)
+   {
+      const Refinement &ref = refinements[i];
+      if (ref.GetType() != Refinement::XYZ)
+      {
+         Iso = false;
+      }
+   }
+
+   // Reduce the Iso flag over all MPI ranks.
+   bool globalIso = false;
+   MPI_Allreduce(&Iso, &globalIso, 1, MFEM_MPI_CXX_BOOL, MPI_LAND, MyComm);
+
+   if (globalIso) { return false; }
+
+   // In the 3D parallel anisotropic case, check for conflicts on faces.
+   NeighborRefinementMessage::Map send_ref;
+
+   // Create refinement messages to all neighbors (NOTE: some may be empty).
+   Array<int> neighbors;
+   NeighborProcessors(neighbors);
+   for (int i = 0; i < neighbors.Size(); i++)
+   {
+      send_ref[neighbors[i]].SetNCMesh(this);
+   }
+
+   // Populate messages: all refinements that occur next to the processor
+   // boundary need to be sent to the adjoining neighbors so they can keep
+   // their ghost layer up to date.
+   Array<int> ranks;
+   ranks.Reserve(64);
+   for (int i = 0; i < refinements.Size(); i++)
+   {
+      const Refinement &ref = refinements[i];
+      MFEM_ASSERT(ref.index < NElements, "");
+      const int elem = leaf_elements[ref.index];
+      ElementNeighborProcessors(elem, ranks);
+      for (int j = 0; j < ranks.Size(); j++)
+      {
+         send_ref[ranks[j]].AddRefinement(elem, ref.GetType());
+      }
+   }
+
+   // Send the messages (overlap with local refinements)
+   NeighborRefinementMessage::IsendAll(send_ref, MyComm);
+
+   // Note that ghost refinements are not looked up using elemToRef. Local
+   // refinements are recorded first in elemToRef, and ghosts only need to be
+   // compared to local refinements. There is no need for ghost-to-ghost
+   // comparisons.
+   std::map<int, int> elemToRef; // Only for local refinements, not ghosts.
+   for (int i = 0; i < refinements.Size(); i++)
+   {
+      elemToRef[leaf_elements[refinements[i].index]] = i;
+   }
+
+   // Check local refinements
+   for (int i = 0; i < refinements.Size(); i++)
+   {
+      const Refinement &ref = refinements[i];
+      CheckRefinement(leaf_elements[ref.index], ref.GetType(), refinements,
+                      elemToRef, conflicts);
+   }
+
+   // Receive (ghost layer) refinements from all neighbors
+   for (int j = 0; j < neighbors.Size(); j++)
+   {
+      int rank, size;
+      NeighborRefinementMessage::Probe(rank, size, MyComm);
+
+      NeighborRefinementMessage msg;
+      msg.SetNCMesh(this);
+      msg.Recv(rank, size, MyComm);
+
+      // check the ghost refinements
+      for (int i = 0; i < msg.Size(); i++)
+      {
+         CheckRefinement(msg.elements[i], msg.values[i], refinements, elemToRef,
+                         conflicts);
+      }
+   }
+
+   // Make sure we can delete the send buffers
+   NeighborRefinementMessage::WaitAllSent(send_ref);
+
+   CheckRefinementMaster(refinements, elemToRef, conflicts);
+
+   const bool conflict = conflicts.size() > 0;
+   bool globalConflict = false;
+   MPI_Allreduce(&conflict, &globalConflict, 1, MFEM_MPI_CXX_BOOL, MPI_LOR,
+                 MyComm);
+   return globalConflict;
+}
+
+int GetHexFaceDir(int face)
+{
+   // Hexahedron face vertices
+   // From Geometry::Constants<Geometry::CUBE>::FaceVert[6][4] in fem/geom.cpp
+   // {3, 2, 1, 0}, {0, 1, 5, 4}, {1, 2, 6, 5},
+   // {2, 3, 7, 6}, {3, 0, 4, 7}, {4, 5, 6, 7}
+   constexpr std::array<int, 6> hexFaceDir = {2, 1, 0, 1, 0, 2};
+   return hexFaceDir[face];
+}
+
+char GetHexFaceRefType(const bool (&refDir)[3], int face)
+{
+   const int faceDir = GetHexFaceDir(face);
+   std::array<int, 2> faceRefDir;
+   int cnt = 0;
+   for (int d=0; d<3; ++d)
+   {
+      if (d != faceDir)
+      {
+         faceRefDir[cnt] = refDir[d] ? 1 : 0;
+         cnt++;
+      }
+   }
+
+   const char ref_type = (char)(faceRefDir[0] + (2 * faceRefDir[1]));
+   return ref_type;
+}
+
+// Assuming a vertical split of the master face with ordered vertices
+// (vn1, vn2, vn3, vn4), check whether there is a horizontal split among the
+// slave faces of this face. This recursive function is similar to
+// NCMesh::CheckAnisoFace.
+bool ParNCMesh::CheckRefAnisoFaceSplits(int vn1, int vn2, int vn3, int vn4,
+                                        int level)
+{
+   const int mid23 = FindMidEdgeNode(vn2, vn3);
+   const int mid41 = FindMidEdgeNode(vn4, vn1);
+
+   if (mid23 >= 0 && mid41 >= 0) // If horizontally split
+   {
+      const int midf = nodes.FindId(mid23, mid41);
+      if (midf >= 0)
+      {
+         if (CheckRefAnisoFaceSplits(vn1, vn2, mid23, mid41, level + 1))
+         {
+            return true;
+         }
+         if (CheckRefAnisoFaceSplits(mid41, mid23, vn3, vn4, level + 1))
+         {
+            return true;
+         }
+      }
+   }
+
+   if (level > 0) { return true; }
+
+   return false;
+}
+
+void ParNCMesh::CheckRefinementMaster(const Array<Refinement> &refinements,
+                                      const std::map<int, int> &elemToRef,
+                                      std::set<int> &conflicts)
+{
+   MFEM_VERIFY(Dim == 3, "");
+   const NCList &faceList = GetFaceList();
+
+   for (const auto &mf : faceList.masters)
+   {
+      // Check for conflicts only if the master element is marked for refinement
+      if (elemToRef.count(mf.element) == 0) { continue; }
+
+      const int refIndex = elemToRef.at(mf.element);
+      const Refinement& ref = refinements[refIndex];
+
+      bool refDir[3];
+      for (int i=0; i<3; ++i)
+         refDir[i] = ref.s[i] > real_t{0};
+
+      const char faceRefType = GetHexFaceRefType(refDir, mf.local);
+      if (faceRefType == 0) { continue; } // No refinement on this face
+
+      std::array<int, 4> fv;
+      for (int i=0; i<4; ++i)
+      {
+         fv[i] = elements[mf.element].node[
+                    Geometry::Constants<Geometry::CUBE>::FaceVert[mf.local][i]];
+      }
+
+      if (faceRefType != 2) // X or XY split w.r.t. the face.
+      {
+         // Check X face split
+         if (CheckRefAnisoFaceSplits(fv[0], fv[1], fv[2], fv[3]))
+         {
+            conflicts.insert(refIndex);
+         }
+      }
+
+      if (faceRefType != 1) // Y or XY split w.r.t. the face.
+      {
+         // Check Y face split
+         if (CheckRefAnisoFaceSplits(fv[1], fv[2], fv[3], fv[0]))
+         {
+            conflicts.insert(refIndex);
+         }
+      }
+   }
+}
+
+int FindHexFace(const int* no, int vn1, int vn2, int vn3, int vn4)
+{
+   std::set<int> v;
+   v.insert({vn1, vn2, vn3, vn4});
+
+   int face = -1;
+   for (int f=0; f<6; ++f)
+   {
+      bool allFound = true;
+      for (int i=0; i<4; ++i)
+      {
+         const int vi = no[Geometry::Constants<Geometry::CUBE>::FaceVert[f][i]];
+         if (v.count(vi) == 0)
+         {
+            allFound = false;
+         }
+      }
+
+      if (allFound)
+      {
+         MFEM_ASSERT(face == -1, "");
+         face = f;
+      }
+   }
+
+   MFEM_ASSERT(face >= 0, "");
+   return face;
+}
+
+// Assumption: v1 and v2 are indices of hex vertices connected by an edge.
+// The return value is {0,1,2} denoting split {X,Y,Z}.
+int GetHexEdgeSplit(const int* nodes, int v1, int v2)
+{
+   Array<int> v(2);
+   v[0] = v1;
+   v[1] = v2;
+   v.Sort();
+
+   // Find the edge in the hexahedron
+   int edge = -1;
+   Array<int> ev(2);
+   for (int i=0; i<12; ++i)
+   {
+      for (int j=0; j<2; ++j)
+      {
+         ev[j] = nodes[Geometry::Constants<Geometry::CUBE>::Edges[i][j]];
+      }
+      ev.Sort();
+
+      if (ev == v)
+      {
+         MFEM_ASSERT(edge == -1, "");
+         edge = i;
+      }
+   }
+
+   MFEM_ASSERT(edge >= 0, "");
+
+   constexpr int edgeDir[12] = {0, 1, 0, 1, 0, 1, 0, 1, 2, 2, 2, 2};
+   return edgeDir[edge];
+}
+
+void ParNCMesh::CheckRefAnisoFace(int elem, int vn1, int vn2, int vn3, int vn4,
+                                  const Array<Refinement> &refinements,
+                                  const std::map<int, int> &elemToRef,
+                                  std::set<int> &conflicts)
+{
+   Face* face = faces.Find(vn1, vn2, vn3, vn4);
+   if (!face) { return; }
+
+   // Find the neighbor of this face.
+   const int nghbIndex = face->elem[0] == elem ? face->elem[1] : face->elem[0];
+   if (nghbIndex < 0) { return; }
+
+   Element &nghb = elements[nghbIndex];
+   MFEM_ASSERT(nghb.ref_type == 0, "");
+
+   if (elemToRef.count(nghbIndex) > 0)
+   {
+      const int refIndex = elemToRef.at(nghbIndex);
+      const Refinement& ref = refinements[refIndex];
+
+      bool refDir[3];
+      for (int i=0; i<3; ++i)
+         refDir[i] = ref.s[i] > real_t{0};
+
+      const int localFace = FindHexFace(nghb.node, vn1, vn2, vn3, vn4);
+      const int faceDir = GetHexFaceDir(localFace);
+      const char face_ref_type = GetHexFaceRefType(refDir, localFace);
+      const bool faceAniso = face_ref_type == 1 ||
+                             face_ref_type == 2; // X or Y w.r.t. the face.
+
+      if (faceAniso)
+      {
+         // Determine whether the face is anisotropically split in the vertical
+         // direction, with respect to the vertex ordering (vn1, vn2, vn3, vn4).
+         int hexSplitOnFace = -1;
+
+         const int firstFaceDir = face_ref_type == 1 ? 0 : 1;
+
+         int cnt = 0;
+         for (int i=0; i<3; ++i)
+         {
+            if (i == faceDir) { continue; }
+
+            if (firstFaceDir == cnt)
+            {
+               MFEM_ASSERT(hexSplitOnFace == -1, "");
+               hexSplitOnFace = i;
+            }
+
+            cnt++;
+         }
+         MFEM_ASSERT(cnt == 2 && hexSplitOnFace >= 0, "");
+
+         const int edgeSplit = GetHexEdgeSplit(nghb.node, vn1, vn2);
+         if (edgeSplit != hexSplitOnFace) { conflicts.insert(refIndex); }
+      }
+   }
+   // The else case is that the neighbor is not refined, so there is no need to
+   // check for conflicts.
+}
+
+void ParNCMesh::CheckRefIsoFace(int elem, int vn1, int vn2, int vn3, int vn4,
+                                int en1, int en2, int en3, int en4,
+                                const Array<Refinement> &refinements,
+                                const std::map<int, int> &elemToRef,
+                                std::set<int> &conflicts)
+{
+   CheckRefAnisoFace(elem, vn1, vn2, en2, en4, refinements, elemToRef, conflicts);
+   CheckRefAnisoFace(elem, en4, en2, vn3, vn4, refinements, elemToRef, conflicts);
+   CheckRefAnisoFace(elem, vn4, vn1, en1, en3, refinements, elemToRef, conflicts);
+   CheckRefAnisoFace(elem, en3, en1, vn2, vn3, refinements, elemToRef, conflicts);
+}
+
+void ParNCMesh::CheckRefinement(int elem, char ref_type,
+                                const Array<Refinement> &refinements,
+                                const std::map<int, int> &elemToRef,
+                                std::set<int> &conflicts)
+{
+   const Element &el = elements[elem];
+   MFEM_ASSERT(el.geom == Geometry::CUBE && el.ref_type == 0,
+               "Element must be an unrefined hexahedron");
+
+   const int* no = el.node;
+
+   // Check the faces of this element being refined (depends on ref_type).
+   // This follows the logic of NCMesh::RefineElement().
+   if (ref_type == Refinement::X) // split along X axis
+   {
+      CheckRefAnisoFace(elem, no[0], no[1], no[5], no[4], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[2], no[3], no[7], no[6], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[4], no[5], no[6], no[7], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[3], no[2], no[1], no[0], refinements,
+                        elemToRef, conflicts);
+   }
+   else if (ref_type == Refinement::Y) // split along Y axis
+   {
+      CheckRefAnisoFace(elem, no[1], no[2], no[6], no[5], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[3], no[0], no[4], no[7], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[5], no[6], no[7], no[4], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[0], no[3], no[2], no[1], refinements,
+                        elemToRef, conflicts);
+   }
+   else if (ref_type == Refinement::Z) // split along Z axis
+   {
+      CheckRefAnisoFace(elem, no[4], no[0], no[1], no[5], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[5], no[1], no[2], no[6], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[6], no[2], no[3], no[7], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[7], no[3], no[0], no[4], refinements,
+                        elemToRef, conflicts);
+   }
+   else if (ref_type == Refinement::XY) // XY split
+   {
+      CheckRefAnisoFace(elem, no[0], no[1], no[5], no[4], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[1], no[2], no[6], no[5], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[2], no[3], no[7], no[6], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[3], no[0], no[4], no[7], refinements,
+                        elemToRef, conflicts);
+
+      const int mid01 = GetMidEdgeNode(no[0], no[1]);
+      const int mid12 = GetMidEdgeNode(no[1], no[2]);
+      const int mid23 = GetMidEdgeNode(no[2], no[3]);
+      const int mid30 = GetMidEdgeNode(no[3], no[0]);
+
+      const int mid45 = GetMidEdgeNode(no[4], no[5]);
+      const int mid56 = GetMidEdgeNode(no[5], no[6]);
+      const int mid67 = GetMidEdgeNode(no[6], no[7]);
+      const int mid74 = GetMidEdgeNode(no[7], no[4]);
+
+      CheckRefIsoFace(elem, no[3], no[2], no[1], no[0], mid23, mid12, mid01,
+                      mid30, refinements, elemToRef, conflicts);
+      CheckRefIsoFace(elem, no[4], no[5], no[6], no[7], mid45, mid56, mid67,
+                      mid74, refinements, elemToRef, conflicts);
+   }
+   else if (ref_type == Refinement::XZ) // XZ split
+   {
+      CheckRefAnisoFace(elem, no[3], no[2], no[1], no[0], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[2], no[6], no[5], no[1], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[6], no[7], no[4], no[5], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[7], no[3], no[0], no[4], refinements,
+                        elemToRef, conflicts);
+
+      const int mid01 = GetMidEdgeNode(no[0], no[1]);
+      const int mid23 = GetMidEdgeNode(no[2], no[3]);
+      const int mid45 = GetMidEdgeNode(no[4], no[5]);
+      const int mid67 = GetMidEdgeNode(no[6], no[7]);
+
+      const int mid04 = GetMidEdgeNode(no[0], no[4]);
+      const int mid15 = GetMidEdgeNode(no[1], no[5]);
+      const int mid26 = GetMidEdgeNode(no[2], no[6]);
+      const int mid37 = GetMidEdgeNode(no[3], no[7]);
+
+      CheckRefIsoFace(elem, no[0], no[1], no[5], no[4], mid01, mid15, mid45,
+                      mid04, refinements, elemToRef, conflicts);
+      CheckRefIsoFace(elem, no[2], no[3], no[7], no[6], mid23, mid37, mid67,
+                      mid26, refinements, elemToRef, conflicts);
+   }
+   else if (ref_type == Refinement::YZ) // YZ split
+   {
+      const int mid12 = GetMidEdgeNode(no[1], no[2]);
+      const int mid30 = GetMidEdgeNode(no[3], no[0]);
+      const int mid56 = GetMidEdgeNode(no[5], no[6]);
+      const int mid74 = GetMidEdgeNode(no[7], no[4]);
+
+      const int mid04 = GetMidEdgeNode(no[0], no[4]);
+      const int mid15 = GetMidEdgeNode(no[1], no[5]);
+      const int mid26 = GetMidEdgeNode(no[2], no[6]);
+      const int mid37 = GetMidEdgeNode(no[3], no[7]);
+
+      CheckRefAnisoFace(elem, no[4], no[0], no[1], no[5], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[0], no[3], no[2], no[1], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[3], no[7], no[6], no[2], refinements,
+                        elemToRef, conflicts);
+      CheckRefAnisoFace(elem, no[7], no[4], no[5], no[6], refinements,
+                        elemToRef, conflicts);
+
+      CheckRefIsoFace(elem, no[1], no[2], no[6], no[5], mid12, mid26, mid56,
+                      mid15, refinements, elemToRef, conflicts);
+      CheckRefIsoFace(elem, no[3], no[0], no[4], no[7], mid30, mid04, mid74,
+                      mid37, refinements, elemToRef, conflicts);
+   }
+   else if (ref_type == Refinement::XYZ) // XYZ split
+   {
+      const int mid01 = GetMidEdgeNode(no[0], no[1]);
+      const int mid12 = GetMidEdgeNode(no[1], no[2]);
+      const int mid23 = GetMidEdgeNode(no[2], no[3]);
+      const int mid30 = GetMidEdgeNode(no[3], no[0]);
+
+      const int mid45 = GetMidEdgeNode(no[4], no[5]);
+      const int mid56 = GetMidEdgeNode(no[5], no[6]);
+      const int mid67 = GetMidEdgeNode(no[6], no[7]);
+      const int mid74 = GetMidEdgeNode(no[7], no[4]);
+
+      const int mid04 = GetMidEdgeNode(no[0], no[4]);
+      const int mid15 = GetMidEdgeNode(no[1], no[5]);
+      const int mid26 = GetMidEdgeNode(no[2], no[6]);
+      const int mid37 = GetMidEdgeNode(no[3], no[7]);
+
+      CheckRefIsoFace(elem, no[3], no[2], no[1], no[0], mid23, mid12, mid01,
+                      mid30, refinements, elemToRef, conflicts);
+      CheckRefIsoFace(elem, no[0], no[1], no[5], no[4], mid01, mid15, mid45,
+                      mid04, refinements, elemToRef, conflicts);
+      CheckRefIsoFace(elem, no[1], no[2], no[6], no[5], mid12, mid26, mid56,
+                      mid15, refinements, elemToRef, conflicts);
+      CheckRefIsoFace(elem, no[2], no[3], no[7], no[6], mid23, mid37, mid67,
+                      mid26, refinements, elemToRef, conflicts);
+      CheckRefIsoFace(elem, no[3], no[0], no[4], no[7], mid30, mid04, mid74,
+                      mid37, refinements, elemToRef, conflicts);
+      CheckRefIsoFace(elem, no[4], no[5], no[6], no[7], mid45, mid56, mid67,
+                      mid74, refinements, elemToRef, conflicts);
+   }
+   else
+   {
+      MFEM_ABORT("Invalid refinement type.");
+   }
+}
 
 void ParNCMesh::Refine(const Array<Refinement> &refinements)
 {
@@ -1513,14 +2014,14 @@ void ParNCMesh::Refine(const Array<Refinement> &refinements)
       return;
    }
 
-   for (int i = 0; i < refinements.Size(); i++)
+   for (int i = 0; i < refinements.Size() && Iso; i++)
    {
       const Refinement &ref = refinements[i];
-      MFEM_VERIFY(ref.GetType() == 7 || Dim < 3,
-                  "anisotropic parallel refinement not supported yet in 3D.");
+      if (ref.GetType() != Refinement::XYZ)
+      {
+         Iso = false;
+      }
    }
-   MFEM_VERIFY(Iso || Dim < 3,
-               "parallel refinement of 3D aniso meshes not supported yet.");
 
    NeighborRefinementMessage::Map send_ref;
 
@@ -1541,7 +2042,7 @@ void ParNCMesh::Refine(const Array<Refinement> &refinements)
    {
       const Refinement &ref = refinements[i];
       MFEM_ASSERT(ref.index < NElements, "");
-      int elem = leaf_elements[ref.index];
+      const int elem = leaf_elements[ref.index];
       ElementNeighborProcessors(elem, ranks);
       for (int j = 0; j < ranks.Size(); j++)
       {

mesh/pncmesh.hpp

@@ -86,6 +86,12 @@ public:
        date. */
    void Refine(const Array<Refinement> &refinements) override;
 
+   /** See Mesh::AnisotropicConflict. The return value is globally MPI-reduced,
+       whereas @a conflicts contains local indices of conflicting entries of
+       @a refinements. */
+   bool AnisotropicConflict(const Array<Refinement> &refinements,
+                            std::set<int> &conflicts);
+
    /// Parallel version of NCMesh::LimitNCLevel.
    void LimitNCLevel(int max_nc_level) override;
 
@@ -211,8 +217,12 @@ public:
 
    // utility
 
+   /// Return the MPI rank for this process.
    int GetMyRank() const { return MyRank; }
 
+   /// Return true if using more than one MPI process.
+   bool IsParallel() const override { return NRanks > 1; }
+
    /// Use the communication pattern from last Rebalance() to send element DOFs.
    void SendRebalanceDofs(int old_ndofs, const Table &old_element_dofs,
                           long old_global_offset, FiniteElementSpace* space);
@@ -586,6 +596,44 @@ protected: // implementation
 
    std::size_t GroupsMemoryUsage() const;
 
+   // The following functions help with checking for anisotropic refinements in
+   // different directions on a face shared by two hexahedral elements.
+
+   /** For the face with ordered vertices vn* and neighboring element @a elem,
+       check whether the other neighboring element (if it exists) is marked for
+       a horizontal refinement conflicting with a vertical split. */
+   void CheckRefAnisoFace(int elem, int vn1, int vn2, int vn3, int vn4,
+                          const Array<Refinement> &refinements,
+                          const std::map<int, int> &elemToRef,
+                          std::set<int> &conflicts);
+
+   /** For the face with ordered vertices vn*, edge midpoints en*, and
+       neighboring element @a elem, check whether the other neighboring element
+       (if it exists) is marked for a refinement conflicting with an isotropic
+       refinement of the face. */
+   void CheckRefIsoFace(int elem, int vn1, int vn2, int vn3, int vn4,
+                        int en1, int en2, int en3, int en4,
+                        const Array<Refinement> &refinements,
+                        const std::map<int, int> &elemToRef,
+                        std::set<int> &conflicts);
+
+   /// Check whether any master face is marked for a conflicting refinement.
+   void CheckRefinementMaster(const Array<Refinement> &refinements,
+                              const std::map<int, int> &elemToRef,
+                              std::set<int> &conflicts);
+
+   /** Check whether the refinement of the element with index @a elem and type
+       @a ref_type would cause a conflict. */
+   void CheckRefinement(int elem, char ref_type,
+                        const Array<Refinement> &refinements,
+                        const std::map<int, int> &elemToRef,
+                        std::set<int> &conflicts);
+
+   /** For a vertical split of the master face with ordered vertices
+       (vn1, vn2, vn3, vn4), check whether there is a horizontal split among the
+       slave faces. */
+   bool CheckRefAnisoFaceSplits(int vn1, int vn2, int vn3, int vn4,
+                                int level = 0);
    friend class NeighborRowMessage;
    friend class NeighborOrderMessage;
 };

miniapps/dfem/dfem-minimal-surface.cpp

@@ -179,12 +179,9 @@ private:
 
          dres_du->Mult(z, y);
 
-         auto d_y = y.HostReadWrite();
-         const auto d_x = x.HostRead();
-         for (int i = 0; i < minsurface->ess_tdofs.Size(); i++)
-         {
-            d_y[minsurface->ess_tdofs[i]] = d_x[minsurface->ess_tdofs[i]];
-         }
+         // Reuse z as a temporary vector to avoid unnecessary allocations.
+         x.GetSubVector(minsurface->ess_tdofs, z);
+         y.SetSubVector(minsurface->ess_tdofs, z);
       }
 
       // Pointer to the wrapped MinimalSurface operator
@@ -261,12 +258,9 @@ private:
 
          dres_du->Mult(z, y);
 
-         auto d_y = y.HostReadWrite();
-         const auto d_x = x.HostRead();
-         for (int i = 0; i < minsurface->ess_tdofs.Size(); i++)
-         {
-            d_y[minsurface->ess_tdofs[i]] = d_x[minsurface->ess_tdofs[i]];
-         }
+         // Reuse z as a temporary vector to avoid unnecessary allocations.
+         x.GetSubVector(minsurface->ess_tdofs, z);
+         y.SetSubVector(minsurface->ess_tdofs, z);
       }
 
       const MinimalSurface *minsurface = nullptr;

miniapps/gslib/field-diff.cpp

@@ -164,14 +164,14 @@ int main (int argc, char *argv[])
       }
    }
 
-   FindPointsGSLIB finder1, finder2;
+   FindPointsGSLIB finder1(mesh_1), finder2(mesh_2);
    Vector interp_vals_1(pts_cnt), interp_vals_2(pts_cnt);
 
    // First solution.
-   finder1.Interpolate(mesh_1, vxyz, func_1, interp_vals_1);
+   finder1.Interpolate(vxyz, func_1, interp_vals_1);
 
    // Second solution.
-   finder2.Interpolate(mesh_2, vxyz, func_2, interp_vals_2);
+   finder2.Interpolate(vxyz, func_2, interp_vals_2);
 
    // Compute differences between the two sets of values.
    double avg_diff = 0.0, max_diff = 0.0, diff_p;
@@ -245,9 +245,5 @@ int main (int argc, char *argv[])
    const double vol_diff = diff * lf;
    std::cout << "Vol diff: " << vol_diff << std::endl;
 
-   // Free the internal gslib data.
-   finder1.FreeData();
-   finder2.FreeData();
-
    return 0;
 }

miniapps/gslib/field-interp.cpp

@@ -305,8 +305,7 @@ int main (int argc, char *argv[])
 
    // Evaluate source grid function.
    Vector interp_vals(nodes_cnt*tar_ncomp);
-   FindPointsGSLIB finder;
-   finder.Setup(mesh_1);
+   FindPointsGSLIB finder(mesh_1);
    finder.Interpolate(vxyz, *func_source, interp_vals, point_ordering);
 
    // Project the interpolated values to the target FiniteElementSpace.
@@ -398,9 +397,6 @@ int main (int argc, char *argv[])
    func_target.Save(rho_ofs);
    rho_ofs.close();
 
-   // Free the internal gslib data.
-   finder.FreeData();
-
    // Delete remaining memory.
    if (func_source->OwnFEC())
    {

miniapps/gslib/findpts.cpp

@@ -330,8 +330,7 @@ int main (int argc, char *argv[])
 
    // Find and Interpolate FE function values on the desired points.
    Vector interp_vals(pts_cnt*vec_dim);
-   FindPointsGSLIB finder;
-   finder.Setup(mesh);
+   FindPointsGSLIB finder(mesh);
    finder.SetL2AvgType(FindPointsGSLIB::NONE);
    finder.Interpolate(vxyz, field_vals, interp_vals, point_ordering);
    Array<unsigned int> code_out    = finder.GetCode();
@@ -374,9 +373,6 @@ int main (int argc, char *argv[])
         << "\nPoints not found:    " << not_found
         << "\nPoints on faces:     " << face_pts << endl;
 
-   // Free the internal gslib data.
-   finder.FreeData();
-
    delete fec;
 
    return 0;

miniapps/gslib/pfindpts.cpp

@@ -359,10 +359,9 @@ int main (int argc, char *argv[])
 
    // Find and Interpolate FE function values on the desired points.
    Vector interp_vals(pts_cnt*vec_dim);
-   FindPointsGSLIB finder(MPI_COMM_WORLD);
-   finder.Setup(pmesh);
+   FindPointsGSLIB finder(pmesh);
    finder.SetDistanceToleranceForPointsFoundOnBoundary(10);
-   // Enable GPU to CPU fallback for GPUData only if you must use an older
+   // Enable GPU to CPU fallback for GPUData only if you are using an older
    // version of GSLIB.
    // finder.SetGPUtoCPUFallback(true);
    finder.FindPoints(vxyz, point_ordering);
@@ -443,9 +442,6 @@ int main (int argc, char *argv[])
            << endl;
    }
 
-   // Free the internal gslib data.
-   finder.FreeData();
-
    delete fec;
    if (randomization != 0) { delete mesh; }
 

miniapps/gslib/schwarz_ex1.cpp

@@ -325,9 +325,6 @@ int main(int argc, char *argv[])
       }
    }
 
-   // Free the used memory.
-   finder1.FreeData();
-   finder2.FreeData();
    for (int i = 0; i < nmeshes; i++)
    {
       delete a_ar[i];

miniapps/gslib/schwarz_ex1p.cpp

@@ -308,7 +308,6 @@ int main(int argc, char *argv[])
    pmesh->GetNodes()->SetTrueVector();
    Vector vxyz = pmesh->GetNodes()->GetTrueVector();
 
-
    OversetFindPointsGSLIB finder(MPI_COMM_WORLD);
    finder.Setup(*pmesh, color);
 
@@ -431,7 +430,6 @@ int main(int argc, char *argv[])
    }
 
    // 15. Free the used memory.
-   finder.FreeData();
    delete a;
    delete fespace;
    if (order > 0) { delete fec; }

miniapps/meshing/mesh-optimizer.cpp

@@ -105,7 +105,7 @@
 //   2D untangling:
 //     mesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
 //   2D untangling with shifted barrier metric:
-//     mesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -vl 1 -btype 1
+//     mesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -vl 1 -btype 1
 //   3D untangling (the mesh is in the mfem/data GitHub repository):
 //   * mesh-optimizer -m ../../../mfem_data/cube-holes-inv.mesh -o 3 -mid 313 -tid 1 -rtol 1e-5 -li 50 -qo 4 -fd -vl 1
 

miniapps/meshing/phpref.cpp

@@ -4,6 +4,7 @@
 //
 // Sample runs:  mpirun -np 4 phpref -dim 2 -n 1000
 //               mpirun -np 8 phpref -dim 3 -n 200
+//               mpirun -np 8 phpref -dim 3 -n 20 --anisotropic --fixed-order
 //
 // Description:  This example demonstrates h- and p-refinement in a parallel
 //               finite element discretization of the Poisson problem (cf. ex1p)
@@ -49,6 +50,8 @@ int main(int argc, char *argv[])
    bool visualization = true;
    int numIter = 0;
    int dim = 2;
+   bool anisotropic = false;
+   bool fixedOrder = false;
    bool deterministic = true;
    bool projectSolution = false;
 
@@ -63,6 +66,12 @@ int main(int argc, char *argv[])
                   "Enable or disable GLVis visualization.");
    args.AddOption(&numIter, "-n", "--num-iter", "Number of hp-ref iterations");
    args.AddOption(&dim, "-dim", "--dim", "Mesh dimension (2 or 3)");
+   args.AddOption(&anisotropic, "-aniso", "--anisotropic", "-iso",
+                  "--isotropic",
+                  "Whether to use anisotropic refinements");
+   args.AddOption(&fixedOrder, "-fo", "--fixed-order", "-vo",
+                  "--variable-order",
+                  "Whether to fix the finite element order on all elements");
    args.AddOption(&deterministic, "-det", "--deterministic", "-not-det",
                   "--not-deterministic",
                   "Use deterministic random refinements");
@@ -84,6 +93,9 @@ int main(int argc, char *argv[])
       args.PrintOptions(cout);
    }
 
+   MFEM_VERIFY(!anisotropic || fixedOrder,
+               "Variable-order is not supported with anisotropic refinement");
+
    // 3. Enable hardware devices such as GPUs, and programming models such as
    //    CUDA, OCCA, RAJA and OpenMP based on command line options.
    Device device(device_config);
@@ -157,15 +169,23 @@ int main(int argc, char *argv[])
       const int r2 = deterministic ? DetRand(seed) : rand();
       const int elem = r1 % pmesh.GetNE();
       int hp = r2 % 2;
+      char htype = 7;
       MPI_Bcast(&hp, 1, MPI_INT, 0, MPI_COMM_WORLD);
 
+      if (fixedOrder) { hp = 0; } // Only perform h-refinement
+      if (anisotropic)
+      {
+         const int r3 = deterministic ? DetRand(seed) : rand();
+         htype = (r3 % 7) + 1;
+      }
+
       if (myid == 0)
          cout << "hp-refinement iteration " << iter << ": "
-              << hp_char[hp] << "-refinement" << endl;
+              << hp_char[hp] << "-refinement\n";
 
       if (hp == 1)
       {
-         // p-ref
+         // p-refinement
          Array<pRefinement> refs;
          refs.Append(pRefinement(elem, 1));  // Increase the element order by 1
          fespace.PRefineAndUpdate(refs);
@@ -173,9 +193,23 @@ int main(int argc, char *argv[])
       }
       else
       {
-         // h-ref
+         // h-refinement
          Array<Refinement> refs;
-         refs.Append(Refinement(elem));
+         refs.Append(Refinement(elem, htype));
+         if (anisotropic)
+         {
+            std::set<int> conflicts; // Indices in refs of conflicting elements
+            const bool conflict = pmesh.AnisotropicConflict(refs, conflicts);
+            if (conflict)
+            {
+               if (myid == 0)
+                  cout << "Anisotropic conflict on iteration " << iter
+                       << ", retrying\n";
+               iter--;
+               continue;
+            }
+         }
+
          pmesh.GeneralRefinement(refs);
          fespace.Update(false);
          numH++;

miniapps/meshing/pmesh-optimizer.cpp

@@ -107,7 +107,7 @@
 //   2D untangling:
 //     mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
 //   2D untangling with shifted barrier metric:
-//     mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -vl 1 -btype 1
+//     mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -vl 1 -btype 1
 //   3D untangling (the mesh is in the mfem/data GitHub repository):
 //   * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/cube-holes-inv.mesh -o 3 -mid 313 -tid 1 -rtol 1e-5 -li 50 -qo 4 -fd -vl 1
 //   Shape optimization for a Kershaw transformed mesh using partial assembly:

tests/unit/CMakeLists.txt

@@ -20,8 +20,10 @@ set(UNIT_TESTS_SRCS
   dfem/test_diffusion.cpp
   dfem/test_divergence.cpp
   dfem/test_mass.cpp
+  dfem/test_lvector_interface.cpp
   general/test_array.cpp
   general/test_reduction.cpp
+  general/test_scan.cpp
   general/test_arrays_by_name.cpp
   general/test_error.cpp
   general/test_mem.cpp
@@ -144,11 +146,11 @@ set(UNIT_TESTS_SRCS
 # SERIAL CPU TESTS: unit_tests
 #-----------------------------------------------------------
 if (MFEM_USE_CUDA)
-   set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
+    set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
                 PROPERTY LANGUAGE CUDA)
 endif()
 if (MFEM_USE_HIP)
-  set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
+    set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
                PROPERTY HIP_SOURCE_PROPERTY_FORMAT TRUE)
 endif()
 
@@ -175,7 +177,7 @@ COMMAND ${CMAKE_COMMAND} -E copy_directory
 #   make unit_tests
 #   ctest -R unit_tests [-V]
 if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-  add_test(NAME unit_tests COMMAND unit_tests)
+    add_test(NAME unit_tests COMMAND unit_tests)
 endif()
 
 #-----------------------------------------------------------
@@ -183,16 +185,16 @@ endif()
 #-----------------------------------------------------------
 # Create CUDA executable and test
 if (MFEM_USE_CUDA)
-   # gpu_unit_tests
-   set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
-   set_property(SOURCE ${GPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
-   mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-   target_link_libraries(gpu_unit_tests mfem)
-   add_dependencies(gpu_unit_tests copy_data)
-   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
-   if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-      add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
-   endif()
+    # gpu_unit_tests
+    set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
+    set_property(SOURCE ${GPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+    mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+    target_link_libraries(gpu_unit_tests mfem)
+    add_dependencies(gpu_unit_tests copy_data)
+    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
+    if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+        add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
+    endif()
 endif()
 
 #-----------------------------------------------------------
@@ -200,15 +202,15 @@ endif()
 #-----------------------------------------------------------
 # Create HIP 'gpu_unit_tests' executable and test
 if (MFEM_USE_HIP)
-   # gpu_unit_tests
-   set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
-   mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-   target_link_libraries(gpu_unit_tests mfem)
-   add_dependencies(gpu_unit_tests copy_data)
-   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
-   if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-      add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
-   endif()
+    # gpu_unit_tests
+    set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
+    mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+    target_link_libraries(gpu_unit_tests mfem)
+    add_dependencies(gpu_unit_tests copy_data)
+    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
+    if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+        add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
+    endif()
 endif()
 
 #-----------------------------------------------------------
@@ -224,7 +226,7 @@ function(add_serial_miniapp_test name test_uvm)
 
     set(${NAME}_TESTS_SRCS miniapps/test_${name}.cpp)
     if (MFEM_USE_CUDA)
-       set_property(SOURCE ${${NAME}_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+        set_property(SOURCE ${${NAME}_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
     endif(MFEM_USE_CUDA)
 
     mfem_add_executable(${name}_tests_cpu ${${NAME}_TESTS_SRCS})
@@ -244,25 +246,25 @@ function(add_serial_miniapp_test name test_uvm)
     endif()
 
     if (MFEM_USE_CUDA OR MFEM_USE_HIP)
-       mfem_add_executable(${name}_tests_gpu ${${NAME}_TESTS_SRCS})
-       target_compile_definitions(${name}_tests_gpu PUBLIC MFEM_${NAME}_DEVICE="gpu")
-       target_link_libraries(${name}_tests_gpu mfem)
-       add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_gpu)
-       if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-           add_test(NAME ${name}_tests_gpu COMMAND ${name}_tests_gpu)
-       endif()
+        mfem_add_executable(${name}_tests_gpu ${${NAME}_TESTS_SRCS})
+        target_compile_definitions(${name}_tests_gpu PUBLIC MFEM_${NAME}_DEVICE="gpu")
+        target_link_libraries(${name}_tests_gpu mfem)
+        add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_gpu)
+        if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+            add_test(NAME ${name}_tests_gpu COMMAND ${name}_tests_gpu)
+        endif()
 
-       if (test_uvm)
-           mfem_add_executable(${name}_tests_gpu_uvm ${${NAME}_TESTS_SRCS})
-           target_compile_definitions(${name}_tests_gpu_uvm PUBLIC
+        if (test_uvm)
+            mfem_add_executable(${name}_tests_gpu_uvm ${${NAME}_TESTS_SRCS})
+            target_compile_definitions(${name}_tests_gpu_uvm PUBLIC
                MFEM_${NAME}_DEVICE="gpu:uvm")
-           target_link_libraries(${name}_tests_gpu_uvm mfem)
-           add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME}
+            target_link_libraries(${name}_tests_gpu_uvm mfem)
+            add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME}
                ${name}_tests_gpu_uvm)
-           if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-               add_test(NAME ${name}_tests_gpu_uvm COMMAND ${name}_tests_gpu_uvm)
-           endif()
-       endif()
+            if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+                add_test(NAME ${name}_tests_gpu_uvm COMMAND ${name}_tests_gpu_uvm)
+            endif()
+        endif()
     endif()
 endfunction(add_serial_miniapp_test)
 
@@ -277,25 +279,25 @@ add_dependencies(tmop_pa_tests_cpu copy_miniapps_meshing_data)
 #-----------------------------------------------------------
 # Add 'ceed_tests' executable and test; add extra tests 'ceed_test_*'
 if (MFEM_USE_CEED)
-   set(CEED_TESTS_SRCS
+    set(CEED_TESTS_SRCS
       ceed/test_ceed.cpp
       ceed/test_ceed_main.cpp)
-   if (MFEM_USE_CUDA)
-      set_property(SOURCE ${CEED_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
-   endif(MFEM_USE_CUDA)
-   mfem_add_executable(ceed_tests ${CEED_TESTS_SRCS})
-   target_link_libraries(ceed_tests mfem)
-   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ceed_tests)
-   # Add CEED tests
-   add_test(NAME ceed_tests COMMAND ceed_tests)
-   if (MFEM_USE_CUDA)
-      add_test(NAME ceed_tests_cuda_ref
+    if (MFEM_USE_CUDA)
+        set_property(SOURCE ${CEED_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+    endif(MFEM_USE_CUDA)
+    mfem_add_executable(ceed_tests ${CEED_TESTS_SRCS})
+    target_link_libraries(ceed_tests mfem)
+    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ceed_tests)
+    # Add CEED tests
+    add_test(NAME ceed_tests COMMAND ceed_tests)
+    if (MFEM_USE_CUDA)
+        add_test(NAME ceed_tests_cuda_ref
                COMMAND ceed_tests --device ceed-cuda:/gpu/cuda/ref)
-      add_test(NAME ceed_tests_cuda_shared
+        add_test(NAME ceed_tests_cuda_shared
                COMMAND ceed_tests --device ceed-cuda:/gpu/cuda/shared)
-      add_test(NAME ceed_tests_cuda_gen
+        add_test(NAME ceed_tests_cuda_gen
                COMMAND ceed_tests --device ceed-cuda:/gpu/cuda/gen)
-   endif()
+    endif()
 endif()
 
 #-----------------------------------------------------------
@@ -303,54 +305,54 @@ endif()
 #-----------------------------------------------------------
 # Define executables and tests
 if (MFEM_USE_MPI)
-   # punit_tests
-   if (MFEM_USE_CUDA)
-      set_property(SOURCE punit_test_main.cpp PROPERTY LANGUAGE CUDA)
-   endif()
-   mfem_add_executable(punit_tests punit_test_main.cpp ${UNIT_TESTS_SRCS})
-   target_link_libraries(punit_tests mfem)
-   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} punit_tests)
-   foreach(np 1 ${MFEM_MPI_NP})
-      if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-         add_test(NAME punit_tests_np=${np}
+    # punit_tests
+    if (MFEM_USE_CUDA)
+        set_property(SOURCE punit_test_main.cpp PROPERTY LANGUAGE CUDA)
+    endif()
+    mfem_add_executable(punit_tests punit_test_main.cpp ${UNIT_TESTS_SRCS})
+    target_link_libraries(punit_tests mfem)
+    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} punit_tests)
+    foreach(np 1 ${MFEM_MPI_NP})
+        if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+            add_test(NAME punit_tests_np=${np}
                   COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${np}
                   ${MPIEXEC_PREFLAGS} $<TARGET_FILE:punit_tests>
                   ${MPIEXEC_POSTFLAGS})
-      endif()
-   endforeach()
-   if (MFEM_USE_CUDA)
-       # pgpu_unit_tests
-       set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
-       set_property(SOURCE ${PGPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
-       mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-       add_dependencies(pgpu_unit_tests copy_data)
-       target_link_libraries(pgpu_unit_tests mfem)
-       add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
-       foreach(np 1 ${MFEM_MPI_NP})
-          if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-             add_test(NAME pgpu_unit_tests_np=${np}
+        endif()
+    endforeach()
+    if (MFEM_USE_CUDA)
+        # pgpu_unit_tests
+        set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
+        set_property(SOURCE ${PGPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+        mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+        add_dependencies(pgpu_unit_tests copy_data)
+        target_link_libraries(pgpu_unit_tests mfem)
+        add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
+        foreach(np 1 ${MFEM_MPI_NP})
+            if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+                add_test(NAME pgpu_unit_tests_np=${np}
                       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${np}
                       ${MPIEXEC_PREFLAGS} $<TARGET_FILE:pgpu_unit_tests>
                       ${MPIEXEC_POSTFLAGS})
-          endif()
-       endforeach()
+            endif()
+        endforeach()
     endif()
     if (MFEM_USE_HIP)
-       # pgpu_unit_tests
-       set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
-       mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-       add_dependencies(pgpu_unit_tests copy_data)
-       target_link_libraries(pgpu_unit_tests mfem)
-       add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
-       foreach(np 1 ${MFEM_MPI_NP})
-          if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-             add_test(NAME pgpu_unit_tests_np=${np}
+        # pgpu_unit_tests
+        set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
+        mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+        add_dependencies(pgpu_unit_tests copy_data)
+        target_link_libraries(pgpu_unit_tests mfem)
+        add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
+        foreach(np 1 ${MFEM_MPI_NP})
+            if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+                add_test(NAME pgpu_unit_tests_np=${np}
                       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${np}
                       ${MPIEXEC_PREFLAGS} $<TARGET_FILE:pgpu_unit_tests>
                       ${MPIEXEC_POSTFLAGS})
-          endif()
-       endforeach()
-   endif()
+            endif()
+        endforeach()
+    endif()
 endif(MFEM_USE_MPI)
 
 #-----------------------------------------------------------
@@ -424,8 +426,8 @@ endfunction(add_parallel_miniapp_test)
 
 # Additional MPI unit tests
 if (MFEM_USE_MPI)
-   add_parallel_miniapp_test(sedov TRUE)
-   add_parallel_miniapp_test(tmop_pa FALSE)
+    add_parallel_miniapp_test(sedov TRUE)
+    add_parallel_miniapp_test(tmop_pa FALSE)
 endif(MFEM_USE_MPI)
 
 #-----------------------------------------------------------
@@ -434,10 +436,10 @@ endif(MFEM_USE_MPI)
 #-----------------------------------------------------------
 set(DEBUG_DEVICE_SRCS miniapps/test_debug_device.cpp)
 if (MFEM_USE_CUDA)
-   set_property(SOURCE ${DEBUG_DEVICE_SRCS} PROPERTY LANGUAGE CUDA)
+    set_property(SOURCE ${DEBUG_DEVICE_SRCS} PROPERTY LANGUAGE CUDA)
 endif()
 if (MFEM_USE_HIP)
-   set_property(SOURCE ${DEBUG_DEVICE_SRCS}
+    set_property(SOURCE ${DEBUG_DEVICE_SRCS}
                 PROPERTY HIP_SOURCE_PROPERTY_FORMAT TRUE)
 endif()
 mfem_add_executable(debug_device_tests ${DEBUG_DEVICE_SRCS})

tests/unit/dfem/test_diffusion.cpp

@@ -255,8 +255,8 @@ void DFemDiffusion(const char *filename, int p, const int r)
          DOperator dop_mf(vsol, {{Coords, mfes}}, pmesh);
          const auto mf_vector_diffusion_qf =
             [] MFEM_HOST_DEVICE (const tensor<dscalar_t, DIM, DIM> &dudxi,
-               const tensor<real_t, DIM, DIM> &J,
-               const real_t &w)
+                                 const tensor<real_t, DIM, DIM> &J,
+                                 const real_t &w)
          {
             const auto invJ = inv(J), TinJ = transpose(invJ);
             return tuple{ (dudxi * invJ) * TinJ * det(J) * w };

tests/unit/dfem/test_lvector_interface.cpp

@@ -0,0 +1,98 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../unit_tests.hpp"
+#include "mfem.hpp"
+#include <utility>
+
+#ifdef MFEM_USE_MPI
+
+using namespace mfem;
+using namespace mfem::future;
+using mfem::future::tensor;
+
+constexpr int DIM = 3;
+
+namespace kernels
+{
+struct MFApply
+{
+   MFEM_HOST_DEVICE inline auto operator()(const tensor<real_t, DIM> &dudxi,
+                                           const tensor<real_t, DIM, DIM> &J,
+                                           const real_t &w) const
+   {
+      const auto invJ = inv(J);
+      return tuple{ (dudxi * invJ) * transpose(invJ) * det(J) * w };
+   }
+};
+}
+
+TEST_CASE("DFEM L-Vector interface", "[Parallel][DFEM]")
+{
+   constexpr int p = 2; // Polynomial order
+   constexpr int r = 1;
+   constexpr int q = 2 * p + r;
+
+   const auto filename = GENERATE("../../data/fichera.mesh");
+   Mesh smesh(filename);
+   ParMesh pmesh(MPI_COMM_WORLD, smesh);
+   MFEM_VERIFY(pmesh.Dimension() == DIM, "Mesh dimension mismatch");
+
+   pmesh.EnsureNodes();
+   auto *nodes = static_cast<ParGridFunction *>(pmesh.GetNodes());
+   smesh.Clear();
+
+   Array<int> all_domain_attr;
+   if (pmesh.attributes.Size() > 0)
+   {
+      all_domain_attr.SetSize(pmesh.attributes.Max());
+      all_domain_attr = 1;
+   }
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace pfes(&pmesh, &fec);
+   ParFiniteElementSpace *mfes = nodes->ParFESpace();
+
+   const auto *ir = &IntRules.Get(pmesh.GetTypicalElementGeometry(), q);
+
+   ParGridFunction x(&pfes), y(&pfes), z(&pfes);
+   Vector X(pfes.GetTrueVSize()), Y(pfes.GetTrueVSize()), Z(pfes.GetTrueVSize());
+
+   X.Randomize(1);
+   x.SetFromTrueDofs(X);
+
+   ParBilinearForm blf_fa(&pfes);
+   blf_fa.AddDomainIntegrator(new DiffusionIntegrator(ir));
+   blf_fa.Assemble();
+   blf_fa.Finalize();
+
+   static constexpr int U = 0, Coords = 1;
+
+   const auto solution = std::vector{FieldDescriptor{U, &pfes}};
+   DifferentiableOperator dop(solution, {{Coords, mfes}}, pmesh);
+
+   kernels::MFApply mf_apply_qf;
+   dop.AddDomainIntegrator(mf_apply_qf,
+                           tuple{Gradient<U>{}, Gradient<Coords>{}, Weight{}},
+                           tuple{Gradient<U>{}}, *ir, all_domain_attr);
+
+   // Use the L-vector interface to multiply
+   dop.SetMultLevel(DifferentiableOperator::MultLevel::LVECTOR);
+   dop.SetParameters({nodes});
+   dop.Mult(x, z);
+
+   blf_fa.Mult(x, y);
+
+   z -= y;
+   REQUIRE(z.Normlinf() == MFEM_Approx(0.0));
+}
+
+#endif

tests/unit/fem/test_coefficient.cpp

@@ -391,3 +391,42 @@ TEST_CASE("Symmetric Matrix Coefficient", "[Coefficient]")
    // Require equality
    REQUIRE(qf.DistanceTo(values) == MFEM_Approx(0.0));
 }
+
+TEST_CASE("Piecewise Constant Coefficient", "[Coefficient]")
+{
+   Mesh mesh("../../data/beam-quad.mesh");
+
+   QuadratureSpace qs(&mesh, 2);
+   FaceQuadratureSpace qs_f(mesh, 2, FaceType::Boundary);
+   QuadratureFunction qf(qs);
+   QuadratureFunction qf_f(qs_f);
+
+   Vector values({1.0, 2.0, 3.0});
+   PWConstCoefficient coeff(values);
+
+   coeff.Project(qf);
+   for (int e = 0; e < mesh.GetNE(); ++e)
+   {
+      Vector vals;
+      qf.GetValues(e, vals);
+      const int a = mesh.GetAttribute(e);
+      for (const real_t val : vals)
+      {
+         REQUIRE(val == a);
+      }
+   }
+
+   coeff.Project(qf_f);
+   for (int be = 0; be < mesh.GetNBE(); ++be)
+   {
+      const int f = mesh.GetBdrElementFaceIndex(be);
+      const int bf = mesh.GetInvFaceIndices(FaceType::Boundary).at(f);
+      Vector vals;
+      qf_f.GetValues(bf, vals);
+      const int a = mesh.GetBdrAttribute(be);
+      for (const real_t val : vals)
+      {
+         REQUIRE(val == a);
+      }
+   }
+}

tests/unit/general/test_scan.cpp

@@ -0,0 +1,99 @@
+// Copyright (c) 2010-2025, 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 <algorithm>
+#include <limits>
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+// must be included after mfem.hpp
+#include "general/scan.hpp"
+
+using namespace mfem;
+
+TEST_CASE("Inclusive Scan", "[Scan],[GPU]")
+{
+   Array<int> a(10);
+
+   for (int use_dev = 0; use_dev < 2; ++use_dev)
+   {
+      CAPTURE(use_dev);
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i;
+      }
+      auto dptr = a.ReadWrite(use_dev);
+      InclusiveScan(use_dev, dptr, dptr, a.Size());
+      a.HostRead();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         int expected = (i + 1) * i / 2;
+         CAPTURE(i);
+         REQUIRE(AsConst(a)[i] == expected);
+      }
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i + 1;
+      }
+      a.ReadWrite(use_dev);
+      InclusiveScan(use_dev, dptr, dptr, a.Size(), std::multiplies<> {});
+      a.HostRead();
+      int expected = 1;
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         expected *= i + 1;
+         CAPTURE(i);
+         REQUIRE(AsConst(a)[i] == expected);
+      }
+   }
+}
+
+TEST_CASE("Exclusive Scan", "[Scan],[GPU]")
+{
+   Array<int> a(10);
+
+   for (int use_dev = 0; use_dev < 2; ++use_dev)
+   {
+      CAPTURE(use_dev);
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i;
+      }
+      auto dptr = a.ReadWrite(use_dev);
+      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5);
+      a.HostRead();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         int expected = (i + 1) * i / 2 - i + 5;
+         CAPTURE(i);
+         REQUIRE(AsConst(a)[i] == expected);
+      }
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i + 1;
+      }
+      a.ReadWrite(use_dev);
+      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, std::multiplies<> {});
+      a.HostRead();
+      int expected = 5;
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         CAPTURE(i);
+         REQUIRE(AsConst(a)[i] == expected);
+         expected *= i + 1;
+      }
+   }
+}