Vector::DeleteAt on device using InclusiveScan

# Conflicts:
#	mesh/mesh.cpp

.github/workflows/builds-and-tests.yml

@@ -132,14 +132,12 @@ jobs:
             hypre-target: int32
             precision: fp64
             enzyme: true
-            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme) LDFLAGS=-L$LLVM_PREFIX/lib/c++
+            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme)
 
     name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}
 
     runs-on: ${{ matrix.os }}
 
-    continue-on-error: ${{ matrix.enzyme && true || false }}
-
     steps:
       # Fix 'No space left on device' errors for Ubuntu builds.
       - name: Run Actions Cleaner
@@ -294,12 +292,10 @@ jobs:
         run: |
           export HOMEBREW_NO_INSTALL_CLEANUP=1
           brew update
-          brew install enzyme
-          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
-          LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
-          echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
-          echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
-          echo "OMPI_CXX=$LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV
+          brew install llvm@20 enzyme
+          echo "LLVM_PREFIX=$(brew --prefix llvm@20)" >> $GITHUB_ENV
+          echo "OMPI_CC=$(brew --prefix llvm@20)/bin/clang" >> $GITHUB_ENV
+          echo "OMPI_CXX=$(brew --prefix llvm@20)/bin/clang++" >> $GITHUB_ENV
 
       # MFEM build and test
       - name: build

CHANGELOG

@@ -27,19 +27,6 @@ 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
 --------------------
 
@@ -57,10 +44,8 @@ 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
@@ -75,18 +60,10 @@ New and updated examples and miniapps
   operators as smoothers.
   These miniapps can be found in `miniapps/diag-smoothers`.
 
-- 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:
+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`.
 
@@ -94,25 +71,16 @@ 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/bilinearform.cpp

@@ -1275,6 +1275,22 @@ void BilinearForm::Update(FiniteElementSpace *nfes)
    height = width = fes->GetVSize();
 
    if (ext) { ext->Update(); }
+   for (int k = 0; k < domain_integs.Size(); ++k)
+   {
+      domain_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_integs.Size(); ++k)
+   {
+      boundary_integs[k]->Update();
+   }
+   for (int k = 0; k < interior_face_integs.Size(); ++k)
+   {
+      interior_face_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_integs.Size(); ++k)
+   {
+      boundary_face_integs[k]->Update();
+   }
 }
 
 void BilinearForm::SetDiagonalPolicy(DiagonalPolicy policy)
@@ -2337,6 +2353,31 @@ void MixedBilinearForm::Update()
    height = test_fes->GetVSize();
    width = trial_fes->GetVSize();
    if (ext) { ext->Update(); }
+
+   for (int k = 0; k < domain_integs.Size(); ++k)
+   {
+      domain_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_integs.Size(); ++k)
+   {
+      boundary_integs[k]->Update();
+   }
+   for (int k = 0; k < interior_face_integs.Size(); ++k)
+   {
+      interior_face_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_integs.Size(); ++k)
+   {
+      boundary_face_integs[k]->Update();
+   }
+   for (int k = 0; k < trace_face_integs.Size(); ++k)
+   {
+      trace_face_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_trace_face_integs.Size(); ++k)
+   {
+      boundary_trace_face_integs[k]->Update();
+   }
 }
 
 MixedBilinearForm::~MixedBilinearForm()

fem/bilininteg.cpp

@@ -21,6 +21,11 @@ using namespace std;
 namespace mfem
 {
 
+void BilinearFormIntegrator::Update()
+{
+   // default no-op
+}
+
 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&)
 {
    MFEM_ABORT("BilinearFormIntegrator::AssemblePA(fes)\n"
@@ -3460,6 +3465,12 @@ real_t ElasticityIntegrator::ComputeFluxEnergy(const FiniteElement &fluxelem,
    return energy;
 }
 
+void DGTraceIntegrator::Update()
+{
+   qspace[0].reset();
+   qspace[1].reset();
+}
+
 void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
                                            const FiniteElement &el2,
                                            FaceElementTransformations &Trans,

fem/bilininteg.hpp

@@ -46,6 +46,10 @@ public:
    // make sense for the action of the nonlinear operator (but they all make
    // sense for its Jacobian).
 
+   /// Signal this integrator that something about either the trial or test space has changed.
+   virtual void Update();
+
+
    /// Method defining partial assembly.
    /** The result of the partial assembly is stored internally so that it can be
        used later in the methods AddMultPA() and AddMultTransposePA(). */
@@ -3313,6 +3317,7 @@ protected:
    VectorCoefficient *u;
    real_t alpha, beta;
    // PA extension
+   std::unique_ptr<FaceQuadratureSpace> qspace[2];
    Vector pa_data;
    const DofToQuad *maps;             ///< Not owned
    const FaceGeometricFactors *geom;  ///< Not owned
@@ -3335,6 +3340,8 @@ public:
                      real_t a, real_t b)
    { rho = &rho_; u = &u_; alpha = a; beta = b; }
 
+   void Update() override;
+
    using BilinearFormIntegrator::AssembleFaceMatrix;
    void AssembleFaceMatrix(const FiniteElement &el1,
                            const FiniteElement &el2,

fem/coefficient.cpp

@@ -12,7 +12,6 @@
 // Implementation of Coefficient class
 
 #include "fem.hpp"
-#include "../general/forall.hpp"
 
 #include <cmath>
 #include <limits>
@@ -81,49 +80,6 @@ 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,9 +132,6 @@ 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,53 +231,22 @@ 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_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
+   /// @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
    {
       MFEM_ASSERT(!action_callbacks.empty(), "no integrators have been set");
-
-      if (mult_level == MultLevel::LVECTOR)
+      prolongation(solutions, solutions_t, solutions_l);
+      residual_l = 0.0;
+      for (auto &action : action_callbacks)
       {
-         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);
+         action(solutions_l, parameters_l, residual_l);
       }
+      prolongation_transpose(residual_l, result_t);
    }
 
    /// @brief Add a domain integrator to the operator.
@@ -376,8 +345,6 @@ 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;
@@ -387,6 +354,7 @@ 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,24 +1076,6 @@ 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

@@ -401,6 +401,9 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
    {
       MFEM_ABORT("unknown FiniteElementCollection: " << name);
    }
+   MFEM_VERIFY(!strcmp(fec->Name(), name), "input name: \"" << name
+               << "\" does not match the created collection name: \""
+               << fec->Name() << '"');
 
    return fec;
 }
@@ -2515,7 +2518,6 @@ 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,13 +464,6 @@ 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
@@ -482,13 +475,6 @@ 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/fespace.hpp

@@ -922,9 +922,6 @@ public:
    { return mesh->GetBdrElementType(i); }
 
    /// Returns ElementTransformation for the @a i-th element.
-   /// @note The returned pointer references an object owned by the associated
-   /// @a Mesh that will be modified by other calls to `GetElementTransformation`.
-   /// As such, this pointer should @b not be deleted by the caller.
    ElementTransformation *GetElementTransformation(int i) const
    { return mesh->GetElementTransformation(i); }
 

fem/gslib.cpp

@@ -85,9 +85,9 @@ namespace mfem
 {
 
 FindPointsGSLIB::FindPointsGSLIB()
-   : mesh(nullptr),
-     fec_map_lin(nullptr),
-     fdataD(nullptr), cr(nullptr), gsl_comm(nullptr),
+   : mesh(NULL),
+     fec_map_lin(NULL),
+     fdataD(NULL), cr(NULL), gsl_comm(NULL),
      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] = nullptr;
-      ir_split[i] = nullptr;
-      fes_rst_map[i] = nullptr;
-      gf_rst_map[i] = nullptr;
+      mesh_split[i] = NULL;
+      ir_split[i] = NULL;
+      fes_rst_map[i] = NULL;
+      gf_rst_map[i] = NULL;
    }
 
    gsl_comm = new gslib::comm;
@@ -117,40 +117,27 @@ 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()
 {
-   FreeData();
-#ifdef MFEM_USE_MPI
-   if (!Mpi::IsFinalized())  // currently segfaults inside gslib otherwise
-#endif
+   crystal_free(cr);
+   comm_free(gsl_comm);
+   delete gsl_comm;
+   delete cr;
+   for (int i = 0; i < 4; i++)
    {
-      crystal_free(cr);
-      comm_free(gsl_comm);
-      delete gsl_comm;
-      delete cr;
+      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; }
    }
-   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; }
+   if (fec_map_lin) { delete fec_map_lin; fec_map_lin = NULL; }
 }
 
 #ifdef MFEM_USE_MPI
 FindPointsGSLIB::FindPointsGSLIB(MPI_Comm comm_)
-   : mesh(nullptr),
-     fec_map_lin(nullptr),
-     fdataD(nullptr), cr(nullptr), gsl_comm(nullptr),
+   : mesh(NULL),
+     fec_map_lin(NULL),
+     fdataD(NULL), cr(NULL), gsl_comm(NULL),
      dim(-1), points_cnt(-1), setupflag(false), default_interp_value(0),
      avgtype(AvgType::ARITHMETIC), bdr_tol(1e-8)
 {
@@ -160,10 +147,10 @@ FindPointsGSLIB::FindPointsGSLIB(MPI_Comm comm_)
    gf_rst_map.SetSize(4);
    for (int i = 0; i < mesh_split.Size(); i++)
    {
-      mesh_split[i] = nullptr;
-      ir_split[i] = nullptr;
-      fes_rst_map[i] = nullptr;
-      gf_rst_map[i] = nullptr;
+      mesh_split[i] = NULL;
+      ir_split[i] = NULL;
+      fes_rst_map[i] = NULL;
+      gf_rst_map[i] = NULL;
    }
 
    gsl_comm = new gslib::comm;
@@ -171,21 +158,12 @@ 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
@@ -193,9 +171,37 @@ void FindPointsGSLIB::Setup(Mesh &m, const double bb_t, const double newt_tol,
 
    mesh = &m;
    dim  = mesh->Dimension();
-   const unsigned int dof1D = meshOrder+1;
+   unsigned dof1D = meshOrder + 1;
 
-   SetupSplitMeshesAndIntegrationRules(meshOrder);
+   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);
+   }
 
    GetNodalValues(mesh->GetNodes(), gsl_mesh);
 
@@ -1122,18 +1128,13 @@ void FindPointsGSLIB::Interpolate(Mesh &m, const Vector &point_pos,
 void FindPointsGSLIB::FreeData()
 {
    if (!setupflag) { return; }
-#ifdef MFEM_USE_MPI
-   if (!Mpi::IsFinalized())  // currently segfaults inside gslib otherwise
-#endif
+   if (dim == 2)
    {
-      if (dim == 2)
-      {
-         findpts_free_2((gslib::findpts_data_2 *)this->fdataD);
-      }
-      else
-      {
-         findpts_free_3((gslib::findpts_data_3 *)this->fdataD);
-      }
+      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();
@@ -1157,8 +1158,8 @@ void FindPointsGSLIB::FreeData()
 
 void FindPointsGSLIB::SetupSplitMeshes()
 {
-   if (fec_map_lin == nullptr) { fec_map_lin = new H1_FECollection(1, dim); }
-   if (dim == 2)
+   fec_map_lin = new H1_FECollection(1, dim);
+   if (mesh->Dimension() == 2)
    {
       int Nvert = 7;
       int NEsplit = 3;
@@ -1200,7 +1201,7 @@ void FindPointsGSLIB::SetupSplitMeshes()
       mesh_split[1] = new Mesh(Mesh::MakeCartesian2D(1, 1,
                                                      Element::QUADRILATERAL));
    }
-   else if (dim == 3)
+   else if (mesh->Dimension() == 3)
    {
       mesh_split[0] = new Mesh(Mesh::MakeCartesian3D(1, 1, 1,
                                                      Element::HEXAHEDRON));
@@ -1345,6 +1346,41 @@ 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,
@@ -1395,79 +1431,6 @@ 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)
 {
@@ -2118,19 +2081,6 @@ 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)
@@ -2436,10 +2386,6 @@ 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(); }
@@ -2449,7 +2395,35 @@ void OversetFindPointsGSLIB::Setup(Mesh &m, const int meshid,
    const FiniteElement *fe = mesh->GetNodalFESpace()->GetTypicalFE();
    unsigned dof1D = fe->GetOrder() + 1;
 
-   SetupSplitMeshesAndIntegrationRules(meshOrder);
+   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);
+   }
 
    GetNodalValues(mesh->GetNodes(), gsl_mesh);
 
@@ -2506,7 +2480,7 @@ void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
 {
    MFEM_VERIFY(setupflag, "Use OversetFindPointsGSLIB::Setup before "
                "finding points.");
-   MFEM_VERIFY(overset, "Please use OversetFindPoints for overlapping grids.");
+   MFEM_VERIFY(overset, "Please setup FindPoints 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,11 +13,7 @@
 #define MFEM_GSLIB
 
 #include "../config/config.hpp"
-#ifdef MFEM_USE_MPI
-#include "pgridfunc.hpp"
-#else
 #include "gridfunc.hpp"
-#endif
 
 #ifdef MFEM_USE_GSLIB
 
@@ -135,10 +131,6 @@ 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);
 
@@ -198,23 +190,14 @@ 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.
@@ -229,8 +212,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
@@ -306,12 +289,7 @@ 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. 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.
-   */
+       calls MPI_Comm_free() for internal gslib communicators. */
    virtual void FreeData();
 
    /// Return code for each point searched by FindPoints: inside element (0), on
@@ -334,9 +312,6 @@ 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 = (Tr.Elem2No >= 0)?(el2.GetDof()):(0);
+   const int dof2 = el2.GetDof();
 
 #ifdef MFEM_THREAD_SAFE
    // Local storage for element integration
@@ -219,9 +219,7 @@ void HyperbolicFormIntegrator::AssembleFaceVector(
    const IntegrationRule *ir = IntRule;
    if (!ir)
    {
-      const int max_el_order = dof2 ? std::max(el1.GetOrder(),
-                                               el2.GetOrder()) : el1.GetOrder();
-      const int order = 2*max_el_order + IntOrderOffset;
+      const int order = 2*std::max(el1.GetOrder(), el2.GetOrder()) + IntOrderOffset;
       ir = &IntRules.Get(Tr.GetGeometryType(), order);
    }
    // loop over integration points
@@ -233,22 +231,18 @@ 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);
-
-      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);
-      }
+      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.
       {
-         nor(0) = 2*Tr.GetElement1IntPoint().x - 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;
       }
       else
       {
@@ -256,18 +250,14 @@ 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 = (dof2) ? numFlux.Eval(state1, state2, nor, Tr, fluxN):
-                           fluxFunction.ComputeFluxDotN(state1, nor, Tr, fluxN);
+      const real_t speed = numFlux.Eval(state1, state2, 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);
-      if (dof2)
-      {
-         AddMult_a_VWt(+ip.weight*sign, shape2, fluxN, elvect2_mat);
-      }
+      AddMult_a_VWt(+ip.weight*sign, shape2, fluxN, elvect2_mat);
    }
 }
 
@@ -278,7 +268,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 = (Tr.Elem2No >= 0)?(el2.GetDof()):(0);
+   const int dof2 = el2.GetDof();
 
 #ifdef MFEM_THREAD_SAFE
    // Local storage for element integration
@@ -312,9 +302,7 @@ void HyperbolicFormIntegrator::AssembleFaceGrad(
    const IntegrationRule *ir = IntRule;
    if (!ir)
    {
-      const int max_el_order = dof2 ? std::max(el1.GetOrder(),
-                                               el2.GetOrder()) : el1.GetOrder();
-      const int order = 2*max_el_order + IntOrderOffset;
+      const int order = 2*std::max(el1.GetOrder(), el2.GetOrder()) + IntOrderOffset;
       ir = &IntRules.Get(Tr.GetGeometryType(), order);
    }
    // loop over integration points
@@ -324,25 +312,20 @@ void HyperbolicFormIntegrator::AssembleFaceGrad(
 
       Tr.SetAllIntPoints(&ip); // set face and element int. points
 
-      // Calculate basis functions of the first element at the face
+      // 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);
-
-      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);
-      }
+      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.
       {
-         nor(0) = 2*Tr.GetElement1IntPoint().x - 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;
       }
       else
       {
@@ -352,14 +335,7 @@ void HyperbolicFormIntegrator::AssembleFaceGrad(
       // Trial side 1
 
       // Compute hat(J) using evaluated quantities
-      if (dof2)
-      {
-         numFlux.Grad(1, state1, state2, nor, Tr, JDotN);
-      }
-      else
-      {
-         fluxFunction.ComputeFluxJacobianDotN(state1, nor, Tr, JDotN);
-      }
+      numFlux.Grad(1, state1, state2, nor, Tr, JDotN);
 
       const int ioff = fluxFunction.num_equations * dof1;
 
@@ -384,325 +360,36 @@ void HyperbolicFormIntegrator::AssembleFaceGrad(
             }
          }
 
-      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);
-      }
+      // Trial side 2
 
       // Compute hat(J) using evaluated quantities
-      numFlux.Grad(1, state_in, state_out, nor, Tr, JDotN);
+      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 < dof; j++)
-               for (int i = 0; i < dof; i++)
+            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+dof*di, j+dof*dj) += w * shape(i) * shape(j);
+                  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);
+               }
+            }
          }
    }
 }
 
-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,14 +306,12 @@ 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
@@ -333,9 +331,8 @@ private:
 
 public:
    const int num_equations;  // the number of equations
-
    /**
-    * @brief Construct a new HyperbolicFormIntegrator object
+    * @brief Construct a new Hyperbolic Form Integrator object
     *
     * @param[in] numFlux numerical flux
     * @param[in] IntOrderOffset integration order offset
@@ -346,14 +343,21 @@ public:
       const int IntOrderOffset = 0,
       const real_t sign = 1.);
 
-   /// Reset the maximum characteristic speed to zero
-   void ResetMaxCharSpeed() { max_char_speed = 0.0; }
+   /**
+    * @brief Reset the Max Char Speed 0
+    *
+    */
+   void ResetMaxCharSpeed()
+   {
+      max_char_speed = 0.0;
+   }
 
-   /// Get the maximum characteristic speed
-   real_t GetMaxCharSpeed() const { return max_char_speed; }
+   real_t GetMaxCharSpeed()
+   {
+      return max_char_speed;
+   }
 
-   /// Get the associated flux function
-   const FluxFunction &GetFluxFunction() const { return fluxFunction; }
+   const FluxFunction &GetFluxFunction() { return fluxFunction; }
 
    /**
     * @brief Implements (F(u), ∇v) with abstract F computed by
@@ -412,182 +416,6 @@ 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

@@ -147,8 +147,13 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
                                &GetRule(el.GetGeomType(), el.GetOrder(),
                                         *mesh->GetTypicalElementTransformation());
 
+   if (!qspace[static_cast<int>(type)])
+   {
+      qspace[static_cast<int>(type)].reset(
+         new FaceQuadratureSpace(*mesh, *ir, type));
+   }
 
-   FaceQuadratureSpace qs(*mesh, *ir, type);
+   FaceQuadratureSpace& qs = *qspace[static_cast<int>(type)];
    nf = qs.GetNumFaces();
    if (nf==0) { return; }
 

fem/qspace.cpp

@@ -17,12 +17,13 @@ namespace mfem
 {
 
 QuadratureSpaceBase::QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
-                                         const IntegrationRule &ir)
-   : mesh(mesh_), order(ir.GetOrder())
+                                         const IntegrationRule &ir,
+                                         QSpaceStorage storage)
+   : mesh(mesh_), order(ir.GetOrder()), storage(storage)
 {
    for (int g = 0; g < Geometry::NumGeom; g++)
    {
-      int_rule[g] = nullptr;
+      int_rule[g] = NULL;
    }
    int_rule[geom] = &ir;
 }
@@ -37,29 +38,6 @@ 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
 {
 
@@ -119,10 +97,10 @@ void QuadratureSpaceBase::Integrate(VectorCoefficient &coeff,
 
 void QuadratureSpace::ConstructOffsets()
 {
-   const int num_elem = mesh.GetNE();
-   ne = num_elem;
+   const int num_elem = ne;
 
-   if (mesh.GetNumGeometries(mesh.Dimension()) == 1)
+   if (storage == QSpaceStorage::COMPRESSED &&
+       mesh.GetNumGeometries(mesh.Dimension()) == 1)
    {
       Array<Geometry::Type> geoms;
       mesh.GetGeometries(mesh.Dimension(), geoms);
@@ -139,7 +117,7 @@ void QuadratureSpace::ConstructOffsets()
       {
          offsets[i] = offset;
          const Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-         MFEM_ASSERT(int_rule[geom] != nullptr, "Missing integration rule.");
+         MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule.");
          offset += int_rule[geom]->GetNPoints();
       }
       offsets[num_elem] = offset;
@@ -147,14 +125,9 @@ void QuadratureSpace::ConstructOffsets()
    }
 }
 
-void QuadratureSpace::Construct()
-{
-   ConstructIntRules(mesh.Dimension());
-   ConstructOffsets();
-}
-
-QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
-   : QuadratureSpaceBase(*mesh_)
+QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in,
+                                 QSpaceStorage storage)
+   : QuadratureSpaceBase(*mesh_, 0, storage)
 {
    const char *msg = "invalid input stream";
    std::string ident;
@@ -173,15 +146,24 @@ QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
       return;
    }
 
-   Construct();
+   ne = mesh.GetNE();
+   ConstructIntRules(mesh.Dimension());
 }
 
-QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir)
-   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir)
+QuadratureSpace::QuadratureSpace(Mesh *mesh_, int order_, QSpaceStorage storage)
+   : QuadratureSpaceBase(*mesh_, order_, storage)
+{
+   ne = mesh.GetNE();
+   ConstructIntRules(mesh.Dimension());
+}
+
+QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
+                                 QSpaceStorage storage)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir, storage)
 {
    MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension()) <= 1,
                "Constructor not valid for mixed meshes");
-   ConstructOffsets();
+   ne = mesh.GetNE();
 }
 
 void QuadratureSpace::Save(std::ostream &os) const
@@ -203,31 +185,32 @@ const Vector &QuadratureSpace::GetGeometricFactorWeights() const
 }
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, int order_,
-                                         FaceType face_type_)
-   : QuadratureSpaceBase(mesh_, order_), face_type(face_type_),
+                                         FaceType face_type_,
+                                         QSpaceStorage storage)
+   : QuadratureSpaceBase(mesh_, order_, storage), face_type(face_type_),
      face_indices(mesh.GetFaceIndices(face_type_)),
      face_indices_inv(mesh.GetInvFaceIndices(face_type_))
 {
-   Construct();
+   ne = face_indices.Size();
+   ConstructIntRules(mesh.Dimension() - 1);
 }
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                                         FaceType face_type_)
-   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalFaceGeometry(), ir),
-     face_type(face_type_),
-     face_indices(mesh.GetFaceIndices(face_type_)),
+                                         FaceType face_type_,
+                                         QSpaceStorage storage)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalFaceGeometry(), ir, storage),
+     face_type(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");
-   ConstructOffsets();
+   ne = face_indices.Size();
 }
 
 void FaceQuadratureSpace::ConstructOffsets()
 {
-   ne = face_indices.Size();
-
-   if (mesh.GetNumGeometries(mesh.Dimension() - 1) == 1)
+   if (storage == QSpaceStorage::COMPRESSED &&
+       mesh.GetNumGeometries(mesh.Dimension() - 1) == 1)
    {
       Array<Geometry::Type> geoms;
       mesh.GetGeometries(mesh.Dimension() - 1, geoms);
@@ -244,19 +227,13 @@ void FaceQuadratureSpace::ConstructOffsets()
       {
          offsets[i] = offset;
          Geometry::Type geom = mesh.GetFaceGeometry(face_indices[i]);
-         MFEM_ASSERT(int_rule[geom] != nullptr, "Missing integration rule");
+         MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
          offset += int_rule[geom]->GetNPoints();
       }
       offsets[face_indices.Size()] = size = offset;
    }
 }
 
-void FaceQuadratureSpace::Construct()
-{
-   ConstructIntRules(mesh.Dimension() - 1);
-   ConstructOffsets();
-}
-
 int FaceQuadratureSpace::GetPermutedIndex(int idx, int iq) const
 {
    const int f_idx = face_indices[idx];
@@ -274,16 +251,6 @@ 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,7 +19,7 @@
 namespace mfem
 {
 
-enum class QSpaceOffsetStorage
+enum class QSpaceStorage
 {
    FULL,
    COMPRESSED
@@ -31,43 +31,37 @@ enum class QSpaceOffsetStorage
 class QuadratureSpaceBase
 {
 protected:
-   friend class QuadratureFunction; // Uses the offsets.
-
    Mesh &mesh; ///< The underlying mesh.
    int order; ///< The order of integration rule.
-   int size; ///< Total number of quadrature points.
-   int ne; ///< Number of entities
+   int size = -1; ///< Total number of quadrature points. -1 indicates
+   ///< offsets/size not computed yet.
+   int ne; ///< Actual number of entities
    mutable Vector weights; ///< Integration weights.
    mutable long nodes_sequence = 0; ///< Nodes counter for cache invalidation.
 
+   QSpaceStorage storage;
+
    /// @brief Entity quadrature point offset 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.
+   /// has size num_entities + 1. Lazily constructed.
    ///
-   /// 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];
 
    /// Protected constructor. Used by derived classes.
-   QuadratureSpaceBase(Mesh &mesh_, int order_ = 0)
-      : mesh(mesh_), order(order_) { }
+   QuadratureSpaceBase(Mesh &mesh_, int order_ = 0,
+                       QSpaceStorage storage = QSpaceStorage::COMPRESSED)
+      : mesh(mesh_), order(order_), storage(storage)
+   {}
 
    /// Protected constructor. Used by derived classes.
    QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
-                       const IntegrationRule &ir);
+                       const IntegrationRule &ir,
+                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// Fill the @ref int_rule array for each geometry type using @ref order.
    void ConstructIntRules(int dim);
@@ -78,33 +72,49 @@ protected:
    /// Compute the integration weights.
    void ConstructWeights() const;
 
+   virtual void ConstructOffsets() = 0;
+
 public:
+   QSpaceStorage StorageType() const { return storage; }
+
    /// @brief Gets the offset for a given entity @a idx.
    ///
    /// The quadrature point values for entity i are stored in the indices
    /// between Offset(i) and Offset(i+1)
    int Offset(int idx) const
    {
+      if (size < 0)
+      {
+         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
+      }
       return (offsets.Size() == 1) ? (idx * offsets[0]) : offsets[idx];
    }
 
    /// @brief Entity quadrature point offset array.
    ///
-   /// If @a storage is QSpaceOffsetStorage::COMPRESSED, then the returned array
-   /// supports a constant compression scheme for meshes which have a single
+   /// 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.
    ///
-   /// 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;
+   const Array<int> &Offsets() const
+   {
+      if (size < 0)
+      {
+         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
+      }
+      return offsets;
+   }
 
    /// Return the total number of quadrature points.
-   int GetSize() const { return size; }
+   int GetSize() const
+   {
+      if (size < 0)
+      {
+         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
+      }
+      return size;
+   }
 
    /// Return the order of the quadrature rule(s) used by all elements.
    int GetOrder() const { return order; }
@@ -164,19 +174,20 @@ class QuadratureSpace : public QuadratureSpaceBase
 {
 protected:
    const Vector &GetGeometricFactorWeights() const override;
-   void ConstructOffsets();
-   void Construct();
+   void ConstructOffsets() override;
 public:
    /// Create a QuadratureSpace based on the global rules from #IntRules.
-   QuadratureSpace(Mesh *mesh_, int order_)
-      : QuadratureSpaceBase(*mesh_, order_) { Construct(); }
+   QuadratureSpace(Mesh *mesh_, int order_,
+                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// @brief Create a QuadratureSpace with an IntegrationRule, valid only when
    /// the mesh has one element type.
-   QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir);
+   QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
+                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// Read a QuadratureSpace from the stream @a in.
-   QuadratureSpace(Mesh *mesh_, std::istream &in);
+   QuadratureSpace(Mesh *mesh_, std::istream &in,
+                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// Returns number of elements in the mesh.
    inline int GetNE() const { return mesh.GetNE(); }
@@ -221,17 +232,18 @@ class FaceQuadratureSpace : public QuadratureSpaceBase
    const std::unordered_map<int,int> &face_indices_inv;
 
    const Vector &GetGeometricFactorWeights() const override;
-   void ConstructOffsets();
-   void Construct();
+   void ConstructOffsets() override;
 
 public:
    /// Create a FaceQuadratureSpace based on the global rules from #IntRules.
-   FaceQuadratureSpace(Mesh &mesh_, int order_, FaceType face_type_);
+   FaceQuadratureSpace(Mesh &mesh_, int order_, FaceType face_type_,
+                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// @brief Create a FaceQuadratureSpace with an IntegrationRule, valid only
    /// when the mesh has one type of face geometry.
    FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                       FaceType face_type_);
+                       FaceType face_type_,
+                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// Returns number of faces in the mesh.
    inline int GetNumFaces() const { return face_indices.Size(); }
@@ -240,7 +252,8 @@ public:
    FaceType GetFaceType() const { return face_type; }
 
    /// Returns the face transformation of face @a idx.
-   ElementTransformation *GetTransformation(int idx) override;
+   ElementTransformation *GetTransformation(int idx) override
+   { return mesh.GetFaceTransformation(face_indices[idx]); }
 
    /// Returns the geometry type of face @a idx.
    Geometry::Type GetGeometry(int idx) const override

fem/tmop.cpp

@@ -140,36 +140,6 @@ 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)
@@ -193,63 +163,6 @@ 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)
@@ -508,7 +421,7 @@ void TMOP_QualityMetric::DefaultAssembleH(const DenseTensor &H,
          {
             for (int cc = 0; cc < dim; cc++)
             {
-               const real_t entry_rr_cc = Hrc(rr, cc);
+               const double entry_rr_cc = Hrc(rr, cc);
 
                for (int i = 0; i < dof; i++)
                {
@@ -568,30 +481,6 @@ 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,
@@ -756,64 +645,6 @@ 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());
@@ -1019,25 +850,6 @@ 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
@@ -1159,25 +971,6 @@ 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)
@@ -1307,25 +1100,6 @@ 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.
@@ -4322,7 +4096,7 @@ real_t TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
 
          const IntegrationPoint &ip_s = ir_s->IntPoint(s);
          Tpr->SetIntPoint(&ip_s);
-         real_t w = surf_fit_coeff->Eval(*Tpr, ip_s) * surf_fit_normal *
+         double 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,14 +14,10 @@
 
 #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
@@ -73,22 +69,6 @@ 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
@@ -144,12 +124,6 @@ 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,
@@ -247,16 +221,12 @@ public:
 
    real_t EvalW(const DenseMatrix &Jpt) 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;
-
-   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override;
+   void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const override
+   { MFEM_ABORT("Not implemented"); }
 
    void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                  const real_t weight, DenseMatrix &A) const override;
+                  const real_t weight, DenseMatrix &A) const override
+   { MFEM_ABORT("Not implemented"); }
 
    // Compute mu_hat.
    real_t EvalWBarrier(const DenseMatrix &Jpt) const;
@@ -398,10 +368,6 @@ 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;
@@ -411,12 +377,6 @@ 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; }
 };
 
@@ -460,10 +420,6 @@ 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;
@@ -475,14 +431,6 @@ 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).
@@ -531,10 +479,6 @@ 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;
@@ -544,14 +488,6 @@ 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()
    {
-      if (finder) { finder->FreeData(); }
+      finder->FreeData();
       delete finder;
    }
 };

general/array.hpp

@@ -211,6 +211,9 @@ public:
    /// Delete the first entry with value == 'el'.
    inline void DeleteFirst(const T &el);
 
+   /// Delete entries at @a indices, and resize.
+   inline void DeleteAt(const Array<int> &indices);
+
    /// Delete the whole array.
    inline void DeleteAll();
 
@@ -935,6 +938,30 @@ inline void Array<T>::DeleteFirst(const T &el)
    }
 }
 
+template <class T>
+inline void Array<T>::DeleteAt(const Array<int> &indices)
+{
+   // Make a copy of the indices, sorted.
+   Array<int> sorted_indices(indices);
+   sorted_indices.Sort();
+
+   int rm_count = 0;
+   for (int i = 0; i < size; i++)
+   {
+      if (rm_count < sorted_indices.Size() && i == sorted_indices[rm_count])
+      {
+         rm_count++;
+      }
+      else
+      {
+         data[i-rm_count] = data[i]; // shift data rm_count
+      }
+   }
+
+   // Resize to remove tail
+   size -= rm_count;
+}
+
 template <class T>
 inline void Array<T>::DeleteAll()
 {

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 did not previously
+ *  return an item given the two or four indices. If the item didn't previously
  *  exist, the methods creates a new one. The Find(...) methods, on the other
- *  hand, just return NULL or -1 if the item does not exist.
+ *  hand, just return NULL or -1 if the item doesn't 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);
-   /// Deep copy
+   /// @brief Deep copy
    HashTable(const HashTable& other);
-   /// Copy assignment not supported
+   /// @brief 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 corresponds 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 correspond 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 */);
 
-   /** @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.
+   /// 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.
 
        @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 corresponds 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 correspond to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -149,8 +149,9 @@ public:
        @warning This method should only be called if T inherits from Hashed4. */
    int GetId(int p1, int p2, int p3, int p4 = -1);
 
-   /** @brief Item accessor with key (or parents) the pair p1, p2. Return
-       NULL if no value corresponds to the requested key.
+   /// 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.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -159,9 +160,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 NULL 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. Return nullptr if no value
+       correspond to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -172,8 +173,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 NULL if no value corresponds to the requested key.
+   /** @brief Item const accessor with key (or parents) the pair 'p1', 'p2'.
+       Return nullptr if no value correspond to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -182,9 +183,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 NULL if no value corresponds to the
-       requested key.
+   /** @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.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -195,12 +196,10 @@ 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;
 
-   /** @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.
+   /// 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.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -210,9 +209,8 @@ 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. Return -1 if it does not
-       exist. Default construct an item and "id" if no value corresponds to the
-       requested key.
+       the item in the underlying BlockArray<T> object. Default construct an
+       item and id if no value correspond to the requested key.
 
        @param[in] p1 First part of the key.
        @param[in] p2 Second part of the key.
@@ -223,16 +221,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;
 
-   /// Return the number of elements currently stored in the HashTable.
+   /// @brief Return the number of elements currently stored in the HashTable.
    int Size() const { return Base::Size() - unused.Size(); }
 
-   /// Return the total number of ids (used and unused) in the HashTable.
+   /// @brief Return the total number of ids (used and unused) in the HashTable.
    int NumIds() const { return Base::Size(); }
 
-   /// Return the number of free/unused ids in the HashTable.
+   /// @brief Return the number of free/unused ids in the HashTable.
    int NumFreeIds() const { return unused.Size(); }
 
-   /** @brief Return true if item @a id exists in (is used by) the container.
+   /** @brief Return true if item 'id' exists in (is used by) the container.
 
        @param[in] id Index of the item in the underlying BlockArray<T>.
 
@@ -243,13 +241,13 @@ public:
 
        @param[in] id Index of the item in the underlying BlockArray<T>.
 
-       @warning Its @a id will be reused by newly added items. */
+       @warning Its id will be reused by newly added items. */
    void Delete(int id);
 
-   /// Remove all items.
+   /// @brief Remove all items.
    void DeleteAll();
 
-   /** @brief Allocate an item at @a id. Enlarge the underlying BlockArray if
+   /** @brief Allocate an item at 'id'. Enlarge the underlying BlockArray if
        necessary.
 
        @param[in] id Index of the item in the underlying BlockArray<T>.
@@ -257,7 +255,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 @a id has already been allocated. */
+       file. Does nothing if the slot 'id' has already been allocated. */
    void Alloc(int id, int p1, int p2);
 
    /** @brief Reinitialize the internal list of unallocated items.
@@ -289,13 +287,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);
 
-   /// Return total size of allocated memory (tables plus items), in bytes.
+   /// @brief Return total size of allocated memory (tables plus items), in bytes.
    std::size_t MemoryUsage() const;
 
-   /// Write details of the memory usage to the mfem output stream.
+   /// @brief Write details of the memory usage to the mfem output stream.
    void PrintMemoryDetail() const;
 
-   /// Print a histogram of bin sizes for debugging purposes.
+   /// @brief Print a histogram of bin sizes for debugging purposes.
    void PrintStats() const;
 
    class iterator : public Base::iterator
@@ -348,7 +346,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;
 
@@ -386,11 +384,11 @@ protected:
    { return (984120265ul*p1 + 125965121ul*p2 + 495698413ul*p3) & mask; }
 
    // Delete() and Reparent() use one of these:
-   /// Hash function for items of type T that inherit from Hashed2.
+   /// @brief Hash function for items of type T that inherit from Hashed2.
    inline int Hash(const Hashed2& item) const
    { return Hash(item.p1, item.p2); }
 
-   /// Hash function for items of type T that inherit from Hashed4.
+   /// @brief 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); }
 
@@ -417,15 +415,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 @a id into bin @a idx.
+   /** @brief Insert the item 'id' into bin '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 @a 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 '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.
@@ -446,11 +444,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 @a idx.
+   /** @brief Return the size of the bin "idx".
 
        @param[in] idx The index of the bin.
        @return The size of the bin. */

general/mem_manager.cpp

@@ -1384,11 +1384,8 @@ 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 == MemoryType::DEFAULT &&
-                    d_mt == GetDualMemoryType(m.h_mt))),
+                  (m.d_mt == d_mt || (d_mt == MemoryType::DEFAULT &&
+                                      m.d_mt == GetDualMemoryType(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

@@ -12,32 +12,29 @@
 #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)
+#elif 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)
+void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
+                   Array<char> &workspace)
 {
    // 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)
       {
@@ -65,43 +62,25 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t 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,
+/// 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 workspace temporary workspace used for device scans. TODO: replace with
+/// internal temporary workspace once that's added to the memory manager.
 /// @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)
+                   Array<char> &workspace, 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)
       {
@@ -129,42 +108,25 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t 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 workspace temporary workspace used for device scans. TODO: replace with
+/// internal temporary workspace once that's added to the memory manager.
 /// @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)
+                   T init_value, Array<char> &workspace, 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)
       {
@@ -194,31 +156,17 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t 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)
+                   T init_value, Array<char> &workspace)
 {
-   ExclusiveScan(use_dev, d_in, d_out, num_items, init_value, std::plus<> {});
+   ExclusiveScan(use_dev, d_in, d_out, num_items, init_value, workspace,
+                 std::plus<> {});
 }
 
 } // namespace mfem

linalg/vector.cpp

@@ -14,6 +14,7 @@
 #include "../general/forall.hpp"
 #include "../general/reducers.hpp"
 #include "../general/hash.hpp"
+#include "../general/scan.hpp"
 #include "vector.hpp"
 
 #ifdef MFEM_USE_OPENMP
@@ -1252,4 +1253,44 @@ real_t Vector::Sum() const
    return res;
 }
 
+void Vector::DeleteAt(const Array<int> &indices)
+{
+   const bool use_dev = UseDevice();
+
+   Array<int> flag(size);
+   const auto d_flag = flag.Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   {
+      d_flag[i] = true;
+   });
+   const auto d_indices = indices.Read(use_dev);
+   mfem::forall_switch(use_dev, indices.Size(), [=] MFEM_HOST_DEVICE (int i)
+   {
+      d_flag[d_indices[i]] = false;
+   });
+
+   Array<int> out_idx(size);
+   auto d_out_idx = out_idx.Write(use_dev);
+   Array<char> workspace;
+   // Perform inclusive scan so that the last entry is the new size.
+   InclusiveScan(use_dev, d_flag, d_out_idx, size, workspace);
+
+   Vector copy(*this);
+   auto d_in = copy.Read(use_dev);
+   auto d_out = Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   {
+      if (d_flag[i])
+      {
+         // Transform inclusive scan to exclusive by shifting.
+         const int j = (i > 0) ? d_out_idx[i - 1] : 0;
+         d_out[j] = d_in[i];
+      }
+   });
+
+   // Get the new size of the vector. Copy only the last entry.
+   Memory<int> submem(out_idx.GetMemory(), out_idx.Size() - 1, 1);
+   size = submem.Read(MemoryClass::HOST, 1)[0];
+}
+
 } // namespace mfem

linalg/vector.hpp

@@ -171,6 +171,12 @@ public:
    /// Resize the vector to size @a s using the MemoryType of @a v.
    void SetSize(int s, const Vector &v) { SetSize(s, v.GetMemory().GetMemoryType()); }
 
+   /// Update \ref Capacity() to @a res (if less than current), keeping existing entries.
+   void Reserve(int res);
+
+   /// Delete entries at @a indices and resize vector accordingly.
+   void DeleteAt(const Array<int> &indices);
+
    /// Set the Vector data.
    /// @warning This method should be called only when OwnsData() is false.
    void SetData(real_t *d) { data.Wrap(d, data.Capacity(), false); }
@@ -621,6 +627,18 @@ inline void Vector::SetSize(int s, MemoryType mt)
    data.UseDevice(use_dev);
 }
 
+inline void Vector::Reserve(int res)
+{
+   if (res > Capacity())
+   {
+      Memory<real_t> p(res, data.GetMemoryType());
+      p.CopyFrom(data, size);
+      p.UseDevice(data.UseDevice());
+      data.Delete();
+      data = p;
+   }
+}
+
 inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
                                      bool own_mem)
 {

mesh/mesh.cpp

@@ -1909,8 +1909,6 @@ void Mesh::Destroy()
 
    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>();
 }
@@ -1927,11 +1925,6 @@ 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)

mesh/ncmesh.cpp

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

mesh/ncmesh.hpp

@@ -483,9 +483,6 @@ 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,13 +3897,6 @@ 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,22 +815,6 @@ 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,507 +1504,6 @@ 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)
 {
@@ -2014,14 +1513,14 @@ void ParNCMesh::Refine(const Array<Refinement> &refinements)
       return;
    }
 
-   for (int i = 0; i < refinements.Size() && Iso; i++)
+   for (int i = 0; i < refinements.Size(); i++)
    {
       const Refinement &ref = refinements[i];
-      if (ref.GetType() != Refinement::XYZ)
-      {
-         Iso = false;
-      }
+      MFEM_VERIFY(ref.GetType() == 7 || Dim < 3,
+                  "anisotropic parallel refinement not supported yet in 3D.");
    }
+   MFEM_VERIFY(Iso || Dim < 3,
+               "parallel refinement of 3D aniso meshes not supported yet.");
 
    NeighborRefinementMessage::Map send_ref;
 
@@ -2042,7 +1541,7 @@ void ParNCMesh::Refine(const Array<Refinement> &refinements)
    {
       const Refinement &ref = refinements[i];
       MFEM_ASSERT(ref.index < NElements, "");
-      const int elem = leaf_elements[ref.index];
+      int elem = leaf_elements[ref.index];
       ElementNeighborProcessors(elem, ranks);
       for (int j = 0; j < ranks.Size(); j++)
       {

mesh/pncmesh.hpp

@@ -86,12 +86,6 @@ 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;
 
@@ -217,12 +211,8 @@ 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);
@@ -596,44 +586,6 @@ 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,9 +179,12 @@ private:
 
          dres_du->Mult(z, y);
 
-         // Reuse z as a temporary vector to avoid unnecessary allocations.
-         x.GetSubVector(minsurface->ess_tdofs, z);
-         y.SetSubVector(minsurface->ess_tdofs, z);
+         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]];
+         }
       }
 
       // Pointer to the wrapped MinimalSurface operator
@@ -258,9 +261,12 @@ private:
 
          dres_du->Mult(z, y);
 
-         // Reuse z as a temporary vector to avoid unnecessary allocations.
-         x.GetSubVector(minsurface->ess_tdofs, z);
-         y.SetSubVector(minsurface->ess_tdofs, z);
+         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]];
+         }
       }
 
       const MinimalSurface *minsurface = nullptr;

miniapps/gslib/field-diff.cpp

@@ -164,14 +164,14 @@ int main (int argc, char *argv[])
       }
    }
 
-   FindPointsGSLIB finder1(mesh_1), finder2(mesh_2);
+   FindPointsGSLIB finder1, finder2;
    Vector interp_vals_1(pts_cnt), interp_vals_2(pts_cnt);
 
    // First solution.
-   finder1.Interpolate(vxyz, func_1, interp_vals_1);
+   finder1.Interpolate(mesh_1, vxyz, func_1, interp_vals_1);
 
    // Second solution.
-   finder2.Interpolate(vxyz, func_2, interp_vals_2);
+   finder2.Interpolate(mesh_2, 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,5 +245,9 @@ 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,7 +305,8 @@ int main (int argc, char *argv[])
 
    // Evaluate source grid function.
    Vector interp_vals(nodes_cnt*tar_ncomp);
-   FindPointsGSLIB finder(mesh_1);
+   FindPointsGSLIB finder;
+   finder.Setup(mesh_1);
    finder.Interpolate(vxyz, *func_source, interp_vals, point_ordering);
 
    // Project the interpolated values to the target FiniteElementSpace.
@@ -397,6 +398,9 @@ 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,7 +330,8 @@ 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(mesh);
+   FindPointsGSLIB finder;
+   finder.Setup(mesh);
    finder.SetL2AvgType(FindPointsGSLIB::NONE);
    finder.Interpolate(vxyz, field_vals, interp_vals, point_ordering);
    Array<unsigned int> code_out    = finder.GetCode();
@@ -373,6 +374,9 @@ 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,9 +359,10 @@ 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(pmesh);
+   FindPointsGSLIB finder(MPI_COMM_WORLD);
+   finder.Setup(pmesh);
    finder.SetDistanceToleranceForPointsFoundOnBoundary(10);
-   // Enable GPU to CPU fallback for GPUData only if you are using an older
+   // Enable GPU to CPU fallback for GPUData only if you must use an older
    // version of GSLIB.
    // finder.SetGPUtoCPUFallback(true);
    finder.FindPoints(vxyz, point_ordering);
@@ -442,6 +443,9 @@ 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,6 +325,9 @@ 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,6 +308,7 @@ int main(int argc, char *argv[])
    pmesh->GetNodes()->SetTrueVector();
    Vector vxyz = pmesh->GetNodes()->GetTrueVector();
 
+
    OversetFindPointsGSLIB finder(MPI_COMM_WORLD);
    finder.Setup(*pmesh, color);
 
@@ -430,6 +431,7 @@ 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 -vl 1 -btype 1
+//     mesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -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,7 +4,6 @@
 //
 // 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)
@@ -50,8 +49,6 @@ 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;
 
@@ -66,12 +63,6 @@ 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");
@@ -93,9 +84,6 @@ 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);
@@ -169,23 +157,15 @@ 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\n";
+              << hp_char[hp] << "-refinement" << endl;
 
       if (hp == 1)
       {
-         // p-refinement
+         // p-ref
          Array<pRefinement> refs;
          refs.Append(pRefinement(elem, 1));  // Increase the element order by 1
          fespace.PRefineAndUpdate(refs);
@@ -193,23 +173,9 @@ int main(int argc, char *argv[])
       }
       else
       {
-         // h-refinement
+         // h-ref
          Array<Refinement> refs;
-         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;
-            }
-         }
-
+         refs.Append(Refinement(elem));
          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 -vl 1 -btype 1
+//     mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -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,7 +20,6 @@ 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
@@ -146,11 +145,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()
 
@@ -177,7 +176,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()
 
 #-----------------------------------------------------------
@@ -185,16 +184,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()
 
 #-----------------------------------------------------------
@@ -202,15 +201,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()
 
 #-----------------------------------------------------------
@@ -226,7 +225,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})
@@ -246,25 +245,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)
 
@@ -279,25 +278,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()
 
 #-----------------------------------------------------------
@@ -305,54 +304,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)
 
 #-----------------------------------------------------------
@@ -426,8 +425,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)
 
 #-----------------------------------------------------------
@@ -436,10 +435,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

@@ -1,98 +0,0 @@
-// 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,42 +391,3 @@ 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_array.cpp

@@ -124,3 +124,19 @@ TEST_CASE("Array stl-interactions", "[Array]")
       CHECK(x[i] == y[i]);
    }
 }
+
+TEST_CASE("Array delete at indices", "[Array]")
+{
+   Array<int>       test({0,1,2,3,4,5,6,7,8});
+   Array<int> rm_indices({0,    3,4,  6,  8});
+   Array<int>     result({  1,2,    5,  7  });
+
+   test.DeleteAt(rm_indices);
+
+   REQUIRE(test.Size() == result.Size());
+
+   for (int i = 0; i < test.Size(); i++)
+   {
+      CHECK(test[i] == result[i]);
+   }
+}

tests/unit/general/test_scan.cpp

@@ -22,6 +22,7 @@ using namespace mfem;
 
 TEST_CASE("Inclusive Scan", "[Scan],[GPU]")
 {
+   Array<char> workspace;
    Array<int> a(10);
 
    for (int use_dev = 0; use_dev < 2; ++use_dev)
@@ -33,13 +34,13 @@ TEST_CASE("Inclusive Scan", "[Scan],[GPU]")
          a[i] = i;
       }
       auto dptr = a.ReadWrite(use_dev);
-      InclusiveScan(use_dev, dptr, dptr, a.Size());
+      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace);
       a.HostRead();
       for (int i = 0; i < a.Size(); ++i)
       {
          int expected = (i + 1) * i / 2;
          CAPTURE(i);
-         REQUIRE(AsConst(a)[i] == expected);
+         REQUIRE(a[i] == expected);
       }
       a.HostReadWrite();
       for (int i = 0; i < a.Size(); ++i)
@@ -47,20 +48,21 @@ TEST_CASE("Inclusive Scan", "[Scan],[GPU]")
          a[i] = i + 1;
       }
       a.ReadWrite(use_dev);
-      InclusiveScan(use_dev, dptr, dptr, a.Size(), std::multiplies<> {});
+      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace, 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);
+         REQUIRE(a[i] == expected);
       }
    }
 }
 
 TEST_CASE("Exclusive Scan", "[Scan],[GPU]")
 {
+   Array<char> workspace;
    Array<int> a(10);
 
    for (int use_dev = 0; use_dev < 2; ++use_dev)
@@ -72,13 +74,13 @@ TEST_CASE("Exclusive Scan", "[Scan],[GPU]")
          a[i] = i;
       }
       auto dptr = a.ReadWrite(use_dev);
-      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5);
+      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace);
       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);
+         REQUIRE(a[i] == expected);
       }
       a.HostReadWrite();
       for (int i = 0; i < a.Size(); ++i)
@@ -86,13 +88,14 @@ TEST_CASE("Exclusive Scan", "[Scan],[GPU]")
          a[i] = i + 1;
       }
       a.ReadWrite(use_dev);
-      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, std::multiplies<> {});
+      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace,
+                    std::multiplies<> {});
       a.HostRead();
       int expected = 5;
       for (int i = 0; i < a.Size(); ++i)
       {
          CAPTURE(i);
-         REQUIRE(AsConst(a)[i] == expected);
+         REQUIRE(a[i] == expected);
          expected *= i + 1;
       }
    }

tests/unit/linalg/test_vector.cpp

@@ -247,3 +247,21 @@ TEST_CASE("Vector Sum", "[Vector],[GPU]")
 
    REQUIRE(sum_1 == MFEM_Approx(sum_2));
 }
+
+TEST_CASE("Vector delete at indices", "[Vector][GPU]")
+{
+   Vector           test({0,1,2,3,4,5,6,7,8});
+   Array<int> rm_indices({0,    3,4,  6,  8});
+   Vector         result({  1,2,    5,  7  });
+
+   test.UseDevice(true);
+   test.DeleteAt(rm_indices);
+
+   REQUIRE(test.Size() == result.Size());
+
+   test.HostReadWrite();
+   for (int i = 0; i < test.Size(); i++)
+   {
+      CHECK(test[i] == result[i]);
+   }
+}