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

.gitignore

@@ -63,7 +63,6 @@ examples/ex3[0-9]
 examples/ex3[0-9]p
 examples/ex4[0-9]
 examples/ex4[0-9]p
-examples/vector-dg-diffusion
 
 examples/refined.mesh
 examples/displaced.mesh

examples/vector-dg-diffusion.cpp

@@ -1,164 +0,0 @@
-#include "mfem.hpp"
-#include "vector-dg-diffusion.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-class RepeatedCoefficient : public VectorCoefficient
-{
-   Coefficient &coeff;
-public:
-   RepeatedCoefficient(int dim, Coefficient &coeff_)
-      : VectorCoefficient(dim), coeff(coeff_)
-   { }
-   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
-   {
-      V.SetSize(vdim);
-      V = coeff.Eval(T, ip);
-   }
-};
-
-real_t u_fn(const Vector &xvec);
-real_t f_fn(const Vector &xvec);
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int ref_levels = 0;
-   int order = 1;
-   real_t sigma = -1.0;
-   real_t kappa = -1.0;
-   const char *device_config = "cpu";
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ref_levels, "-r", "--refine",
-                  "Number of times to refine the mesh uniformly, -1 for auto.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) >= 0.");
-   args.AddOption(&sigma, "-s", "--sigma",
-                  "One of the three DG penalty parameters, typically +1/-1."
-                  " See the documentation of class DGDiffusionIntegrator.");
-   args.AddOption(&kappa, "-k", "--kappa",
-                  "One of the three DG penalty parameters, should be positive."
-                  " Negative values are replaced with (order+1)^2.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.ParseCheck();
-
-   if (kappa < 0)
-   {
-      kappa = (order+1)*(order+1);
-   }
-
-   Device device(device_config);
-   device.Print();
-
-   Mesh mesh(mesh_file);
-   const int dim = mesh.Dimension();
-
-   {
-      if (ref_levels < 0)
-      {
-         ref_levels = (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
-      }
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   DG_FECollection fec(order, dim);
-   FiniteElementSpace fespace(&mesh, &fec, dim);
-   cout << "Number of unknowns: " << fespace.GetVSize() << endl;
-
-
-   FunctionCoefficient scalar_f_coeff(f_fn);
-   FunctionCoefficient scalar_u_coeff(u_fn);
-   RepeatedCoefficient f_coeff(dim, scalar_f_coeff);
-   RepeatedCoefficient u_coeff(dim, scalar_u_coeff);
-
-   ConstantCoefficient one(1.0);
-   ConstantCoefficient zero(5.0);
-   RepeatedCoefficient zero_vec(dim, zero);
-
-   LinearForm b(&fespace);
-   b.AddDomainIntegrator(new VectorDomainLFIntegrator(f_coeff));
-   b.AddBdrFaceIntegrator(
-      new VectorDGDirichletLFIntegrator(u_coeff, one, sigma, kappa));
-   b.Assemble();
-
-   GridFunction x(&fespace);
-   x = 0.0;
-
-   BilinearForm a(&fespace);
-   a.AddDomainIntegrator(new VectorDiffusionIntegrator(one));
-   a.AddInteriorFaceIntegrator(new VectorDGDiffusionIntegrator(
-                                  one, sigma, kappa, dim));
-   a.AddBdrFaceIntegrator(new VectorDGDiffusionIntegrator(
-                             one, sigma, kappa, dim));
-   a.Assemble();
-   a.Finalize();
-
-   const SparseMatrix &A = a.SpMat();
-#ifndef MFEM_USE_SUITESPARSE
-   GSSmoother M(A);
-   if (sigma == -1.0)
-   {
-      PCG(A, M, b, x, 1, 500, 1e-12, 0.0);
-   }
-   else
-   {
-      GMRES(A, M, b, x, 1, 500, 10, 1e-12, 0.0);
-   }
-#else
-   UMFPackSolver umf_solver;
-   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
-   umf_solver.SetOperator(A);
-   umf_solver.Mult(b, x);
-#endif
-
-   ParaViewDataCollection pv("DGDiffusion", &mesh);
-   pv.SetPrefixPath("ParaView");
-   pv.SetHighOrderOutput(true);
-   pv.SetLevelsOfDetail(order);
-   pv.RegisterField("u", &x);
-   pv.SetCycle(0);
-   pv.SetTime(0.0);
-   pv.Save();
-
-   cout << "L2 error: " << x.ComputeL2Error(u_coeff) << '\n';
-
-   return 0;
-}
-
-constexpr real_t pi = M_PI;
-constexpr real_t pi2 = pi*pi;
-
-real_t u_fn(const Vector &xvec)
-{
-   int dim = xvec.Size();
-   real_t x = pi*xvec[0], y = pi*xvec[1];
-   if (dim == 2) { return sin(x)*sin(y); }
-   else { real_t z = pi*xvec[2]; return sin(x)*sin(y)*sin(z); }
-}
-
-real_t f_fn(const Vector &xvec)
-{
-   int dim = xvec.Size();
-   real_t x = pi*xvec[0], y = pi*xvec[1];
-
-   if (dim == 2)
-   {
-      return 2*pi2*sin(x)*sin(y);
-   }
-   else // dim == 3
-   {
-      real_t z = pi*xvec[2];
-      return 3*pi2*sin(x)*sin(y)*sin(z);
-   }
-}

examples/vector-dg-diffusion.hpp

@@ -1,433 +0,0 @@
-#include "mfem.hpp"
-
-namespace mfem
-{
-
-class VectorDGDiffusionIntegrator : public BilinearFormIntegrator
-{
-protected:
-   Coefficient *Q = nullptr;
-   MatrixCoefficient *MQ = nullptr;
-   real_t sigma, kappa;
-   int vdim;
-
-   // these are not thread-safe!
-   Vector shape1, shape2, dshape1dn, dshape2dn, nor, nh, ni;
-   DenseMatrix jmat, dshape1, dshape2, mq, adjJ;
-
-public:
-   VectorDGDiffusionIntegrator(real_t s, real_t k, int vd=-1)
-      : sigma(s), kappa(k), vdim(vd) { }
-   VectorDGDiffusionIntegrator(Coefficient &q, real_t s, real_t k, int vd=-1)
-      : Q(&q), sigma(s), kappa(k), vdim(vd) { }
-   VectorDGDiffusionIntegrator(MatrixCoefficient &mq, real_t s, real_t k,
-                               int vd=-1)
-      : MQ(&mq), sigma(s), kappa(k), vdim(vd) { }
-
-   using BilinearFormIntegrator::AssembleFaceMatrix;
-   virtual void AssembleFaceMatrix(const FiniteElement &el1,
-                                   const FiniteElement &el2,
-                                   FaceElementTransformations &Trans,
-                                   DenseMatrix &full_elmat);
-};
-
-class VectorDGDirichletLFIntegrator : public LinearFormIntegrator
-{
-protected:
-   VectorCoefficient &uD;
-   Coefficient *Q = nullptr;
-   MatrixCoefficient *MQ = nullptr;
-   real_t sigma, kappa;
-   int vdim;
-
-   // these are not thread-safe!
-   Vector shape, dshape_dn, nor, nh, ni, uD_vec;
-   DenseMatrix dshape, mq, adjJ;
-
-public:
-   VectorDGDirichletLFIntegrator(VectorCoefficient &u, real_t s, real_t k,
-                                 int vd=-1)
-      : uD(u), sigma(s), kappa(k), vdim(vd) { }
-   VectorDGDirichletLFIntegrator(VectorCoefficient &u, Coefficient &q, real_t s,
-                                 real_t k, int vd=-1)
-      : uD(u), Q(&q), sigma(s), kappa(k), vdim(vd) { }
-   VectorDGDirichletLFIntegrator(VectorCoefficient &u, MatrixCoefficient &mq,
-                                 real_t s, real_t k, int vd=-1)
-      : uD(u), MQ(&mq), sigma(s), kappa(k), vdim(vd) { }
-
-   using LinearFormIntegrator::AssembleRHSElementVect;
-
-   void AssembleRHSElementVect(const FiniteElement &el,
-                               ElementTransformation &Tr,
-                               Vector &elvect) override
-   { MFEM_ABORT("Not implemented."); }
-
-   void AssembleRHSElementVect(const FiniteElement &el,
-                               FaceElementTransformations &Tr,
-                               Vector &elvect) override;
-};
-
-void VectorDGDiffusionIntegrator::AssembleFaceMatrix(
-   const FiniteElement &el1, const FiniteElement &el2,
-   FaceElementTransformations &Trans, DenseMatrix &full_elmat)
-{
-   int dim, ndof1, ndof2, ndofs;
-   bool kappa_is_nonzero = (kappa != 0.);
-   real_t w, wq = 0.0;
-
-   const int sdim = Trans.GetSpaceDim();
-   if (vdim < 0) { vdim = sdim; }
-
-   dim = el1.GetDim();
-   ndof1 = el1.GetDof();
-
-   nor.SetSize(dim);
-   nh.SetSize(dim);
-   ni.SetSize(dim);
-   adjJ.SetSize(dim);
-   if (MQ)
-   {
-      mq.SetSize(dim);
-   }
-
-   shape1.SetSize(ndof1);
-   dshape1.SetSize(ndof1, dim);
-   dshape1dn.SetSize(ndof1);
-   if (Trans.Elem2No >= 0)
-   {
-      ndof2 = el2.GetDof();
-      shape2.SetSize(ndof2);
-      dshape2.SetSize(ndof2, dim);
-      dshape2dn.SetSize(ndof2);
-   }
-   else
-   {
-      ndof2 = 0;
-   }
-
-   ndofs = ndof1 + ndof2;
-   DenseMatrix elmat;
-   elmat.SetSize(ndofs);
-   elmat = 0.0;
-   if (kappa_is_nonzero)
-   {
-      jmat.SetSize(ndofs);
-      jmat = 0.;
-   }
-
-   const IntegrationRule *ir = IntRule;
-   if (ir == nullptr)
-   {
-      // a simple choice for the integration order
-      int order;
-      if (ndof2)
-      {
-         order = 2 * std::max(el1.GetOrder(), el2.GetOrder());
-      }
-      else
-      {
-         order = 2 * el1.GetOrder();
-      }
-      ir = &IntRules.Get(Trans.GetGeometryType(), order);
-   }
-
-   // assemble: < {(Q \nabla u).n},[v] >      --> elmat
-   //           kappa < {h^{-1} Q} [u],[v] >  --> jmat
-   for (int p = 0; p < ir->GetNPoints(); p++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(p);
-
-      // Set the integration point in the face and the neighboring elements
-      Trans.SetAllIntPoints(&ip);
-
-      // Access the neighboring elements' integration points
-      // Note: eip2 will only contain valid data if Elem2 exists
-      const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
-      const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
-
-      if (dim == 1)
-      {
-         nor(0) = 2 * eip1.x - 1.0;
-      }
-      else
-      {
-         CalcOrtho(Trans.Jacobian(), nor);
-      }
-
-      el1.CalcShape(eip1, shape1);
-      el1.CalcDShape(eip1, dshape1);
-      w = ip.weight / Trans.Elem1->Weight();
-      if (ndof2)
-      {
-         w /= 2;
-      }
-      if (!MQ)
-      {
-         if (Q)
-         {
-            w *= Q->Eval(*Trans.Elem1, eip1);
-         }
-         ni.Set(w, nor);
-      }
-      else
-      {
-         nh.Set(w, nor);
-         MQ->Eval(mq, *Trans.Elem1, eip1);
-         mq.MultTranspose(nh, ni);
-      }
-      CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
-      adjJ.Mult(ni, nh);
-      if (kappa_is_nonzero)
-      {
-         wq = ni * nor;
-      }
-
-      // Note: in the jump term, we use 1/h1 = |nor|/det(J1) which is
-      // independent of Loc1 and always gives the size of element 1 in
-      // direction perpendicular to the face. Indeed, for linear transformation
-      //
-      //      |nor|=measure(face)/measure(ref. face),
-      //
-      //      det(J1)=measure(element)/measure(ref. element),
-      //
-      //      and the ratios measure(ref. element)/measure(ref. face)
-      //      are compatible for all element/face pairs.
-      //
-      // For example: meas(ref. tetrahedron)/meas(ref. triangle) = 1/3, and
-      // for any tetrahedron vol(tet)=(1/3)*height*area(base).
-      //
-      // For interior faces: q_e/h_e=(q1/h1+q2/h2)/2.
-
-      dshape1.Mult(nh, dshape1dn);
-      for (int i = 0; i < ndof1; i++)
-         for (int j = 0; j < ndof1; j++)
-         {
-            elmat(i, j) += shape1(i) * dshape1dn(j);
-         }
-
-      if (ndof2)
-      {
-         el2.CalcShape(eip2, shape2);
-         el2.CalcDShape(eip2, dshape2);
-         w = ip.weight / 2 / Trans.Elem2->Weight();
-         if (!MQ)
-         {
-            if (Q)
-            {
-               w *= Q->Eval(*Trans.Elem2, eip2);
-            }
-            ni.Set(w, nor);
-         }
-         else
-         {
-            nh.Set(w, nor);
-            MQ->Eval(mq, *Trans.Elem2, eip2);
-            mq.MultTranspose(nh, ni);
-         }
-         CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
-         adjJ.Mult(ni, nh);
-         if (kappa_is_nonzero)
-         {
-            wq += ni * nor;
-         }
-
-         dshape2.Mult(nh, dshape2dn);
-
-         for (int i = 0; i < ndof1; i++)
-            for (int j = 0; j < ndof2; j++)
-            {
-               elmat(i, ndof1 + j) += shape1(i) * dshape2dn(j);
-            }
-
-         for (int i = 0; i < ndof2; i++)
-            for (int j = 0; j < ndof1; j++)
-            {
-               elmat(ndof1 + i, j) -= shape2(i) * dshape1dn(j);
-            }
-
-         for (int i = 0; i < ndof2; i++)
-            for (int j = 0; j < ndof2; j++)
-            {
-               elmat(ndof1 + i, ndof1 + j) -= shape2(i) * dshape2dn(j);
-            }
-      }
-
-      if (kappa_is_nonzero)
-      {
-         // only assemble the lower triangular part of jmat
-         wq *= kappa;
-         for (int i = 0; i < ndof1; i++)
-         {
-            const real_t wsi = wq * shape1(i);
-            for (int j = 0; j <= i; j++)
-            {
-               jmat(i, j) += wsi * shape1(j);
-            }
-         }
-         if (ndof2)
-         {
-            for (int i = 0; i < ndof2; i++)
-            {
-               const int i2 = ndof1 + i;
-               const real_t wsi = wq * shape2(i);
-               for (int j = 0; j < ndof1; j++)
-               {
-                  jmat(i2, j) -= wsi * shape1(j);
-               }
-               for (int j = 0; j <= i; j++)
-               {
-                  jmat(i2, ndof1 + j) += wsi * shape2(j);
-               }
-            }
-         }
-      }
-   }
-
-   // elmat := -elmat + sigma*elmat^t + jmat
-   if (kappa_is_nonzero)
-   {
-      for (int i = 0; i < ndofs; i++)
-      {
-         for (int j = 0; j < i; j++)
-         {
-            real_t aij = elmat(i, j), aji = elmat(j, i), mij = jmat(i, j);
-            elmat(i, j) = sigma * aji - aij + mij;
-            elmat(j, i) = sigma * aij - aji + mij;
-         }
-         elmat(i, i) = (sigma - 1.) * elmat(i, i) + jmat(i, i);
-      }
-   }
-   else
-   {
-      for (int i = 0; i < ndofs; i++)
-      {
-         for (int j = 0; j < i; j++)
-         {
-            real_t aij = elmat(i, j), aji = elmat(j, i);
-            elmat(i, j) = sigma * aji - aij;
-            elmat(j, i) = sigma * aij - aji;
-         }
-         elmat(i, i) *= (sigma - 1.);
-      }
-   }
-
-   // populate full matrix following github issue #2909
-   full_elmat.SetSize(vdim*(ndof1 + ndof2));
-   full_elmat = 0.0;
-   for (int d=0; d<vdim; ++d)
-   {
-      for (int j=0; j<ndofs; ++j)
-      {
-         int jj = (j < ndof1) ? j + d*ndof1 : j - ndof1 + d*ndof2 + vdim*ndof1;
-         for (int i=0; i<ndofs; ++i)
-         {
-            int ii = (i < ndof1) ? i + d*ndof1 : i - ndof1 + d*ndof2 + vdim*ndof1;
-            full_elmat(ii, jj) += elmat(i, j);
-         }
-      }
-   }
-};
-
-void VectorDGDirichletLFIntegrator::AssembleRHSElementVect(
-   const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
-{
-   const int dim = el.GetDim();
-   const int sdim = Tr.GetSpaceDim();
-
-   if (vdim < 0) { vdim = sdim; }
-
-   const int ndof = el.GetDof();
-
-   bool kappa_is_nonzero = (kappa != 0.);
-   real_t w;
-
-   nor.SetSize(dim);
-   nh.SetSize(dim);
-   ni.SetSize(dim);
-   adjJ.SetSize(dim);
-   if (MQ)
-   {
-      mq.SetSize(dim);
-   }
-
-   shape.SetSize(ndof);
-   dshape.SetSize(ndof, dim);
-   dshape_dn.SetSize(ndof);
-
-   elvect.SetSize(vdim * ndof);
-   elvect = 0.0;
-
-   const IntegrationRule *ir = IntRule;
-   if (ir == NULL)
-   {
-      // a simple choice for the integration order; is this OK?
-      int order = 2*el.GetOrder();
-      ir = &IntRules.Get(Tr.GetGeometryType(), order);
-   }
-
-   for (int p = 0; p < ir->GetNPoints(); p++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(p);
-
-      // Set the integration point in the face and the neighboring element
-      Tr.SetAllIntPoints(&ip);
-
-      // Access the neighboring element's integration point
-      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
-
-      uD.Eval(uD_vec, Tr, ip);
-
-      if (dim == 1)
-      {
-         nor(0) = 2*eip.x - 1.0;
-      }
-      else
-      {
-         CalcOrtho(Tr.Jacobian(), nor);
-      }
-
-      el.CalcShape(eip, shape);
-      el.CalcDShape(eip, dshape);
-
-      // compute uD through the face transformation
-      w = ip.weight / Tr.Elem1->Weight();
-      if (!MQ)
-      {
-         if (Q)
-         {
-            w *= Q->Eval(*Tr.Elem1, eip);
-         }
-         ni.Set(w, nor);
-      }
-      else
-      {
-         nh.Set(w, nor);
-         MQ->Eval(mq, *Tr.Elem1, eip);
-         mq.MultTranspose(nh, ni);
-      }
-      CalcAdjugate(Tr.Elem1->Jacobian(), adjJ);
-      adjJ.Mult(ni, nh);
-
-      dshape.Mult(nh, dshape_dn);
-
-      for (int vd = 0; vd < vdim; ++vd)
-      {
-         for (int i = 0; i < ndof; ++i)
-         {
-            elvect[i + vd*ndof] += sigma * uD_vec[vd] * dshape_dn[i];
-         }
-      }
-      if (kappa_is_nonzero)
-      {
-         for (int vd = 0; vd < vdim; ++vd)
-         {
-            for (int i = 0; i < ndof; ++i)
-            {
-               elvect[i + vd*ndof] += kappa*(ni*nor) * uD_vec[vd] * shape[i];
-            }
-         }
-      }
-   }
-}
-
-} // namespace mfem

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

@@ -23,6 +23,8 @@
 
 namespace mfem
 {
+class QuadratureSpace;
+class FaceQuadratureSpace;
 
 /// Abstract base class BilinearFormIntegrator
 class BilinearFormIntegrator : public NonlinearFormIntegrator
@@ -44,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(). */
@@ -3311,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
@@ -3333,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/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;
 }

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/integ/bilininteg_dgtrace_pa.cpp

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

fem/qspace.cpp

@@ -17,8 +17,9 @@ 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++)
    {
@@ -96,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);
@@ -124,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;
@@ -150,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
@@ -180,55 +185,53 @@ const Vector &QuadratureSpace::GetGeometricFactorWeights() const
 }
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, int order_,
-                                         FaceType face_type_)
-   : QuadratureSpaceBase(mesh_, order_),
-     face_type(face_type_),
-     num_faces(mesh.GetNFbyType(face_type))
+                                         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_),
-     num_faces(mesh.GetNFbyType(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()
 {
-   face_indices.SetSize(num_faces);
-   offsets.SetSize(num_faces + 1);
-   ne = num_faces;
-   int offset = 0;
-   int f_idx = 0;
-   for (int i = 0; i < mesh.GetNumFacesWithGhost(); i++)
+   if (storage == QSpaceStorage::COMPRESSED &&
+       mesh.GetNumGeometries(mesh.Dimension() - 1) == 1)
    {
-      const Mesh::FaceInformation face = mesh.GetFaceInformation(i);
-      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(face_type))
-      {
-         continue;
-      }
-      face_indices[f_idx] = i;
-      face_indices_inv[i] = f_idx;
-      offsets[f_idx] = offset;
-      Geometry::Type geom = mesh.GetFaceGeometry(i);
-      MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
-      offset += int_rule[geom]->GetNPoints();
-
-      f_idx++;
+      Array<Geometry::Type> geoms;
+      mesh.GetGeometries(mesh.Dimension() - 1, geoms);
+      offsets.SetSize(1);
+      offsets.HostWrite();
+      offsets[0] = int_rule[geoms[0]]->GetNPoints();
+      size = ne * offsets[0];
+   }
+   else
+   {
+      offsets.SetSize(face_indices.Size() + 1);
+      int offset = 0;
+      for (int i = 0; i < mesh.GetNFbyType(face_type); ++i)
+      {
+         offsets[i] = offset;
+         Geometry::Type geom = mesh.GetFaceGeometry(face_indices[i]);
+         MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
+         offset += int_rule[geom]->GetNPoints();
+      }
+      offsets[face_indices.Size()] = size = offset;
    }
-   offsets[num_faces] = size = offset;
-}
-
-void FaceQuadratureSpace::Construct()
-{
-   ConstructIntRules(mesh.Dimension() - 1);
-   ConstructOffsets();
 }
 
 int FaceQuadratureSpace::GetPermutedIndex(int idx, int iq) const

fem/qspace.hpp

@@ -19,39 +19,49 @@
 namespace mfem
 {
 
+enum class QSpaceStorage
+{
+   FULL,
+   COMPRESSED
+};
+
 /// Abstract base class for QuadratureSpace and FaceQuadratureSpace.
 /** This class represents the storage layout for QuadratureFunction%s, that may
     be defined either on mesh elements or mesh faces. */
 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 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.
    ///
    Array<int> offsets;
    /// 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);
@@ -62,13 +72,21 @@ 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];
    }
 
@@ -79,10 +97,24 @@ public:
    /// can be computed as i * offsets[0], where i is the entity index. Otherwise
    /// has size num_entities + 1.
    ///
-   const Array<int> &Offsets() const { return offsets; }
+   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; }
@@ -142,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(); }
@@ -191,29 +224,29 @@ public:
 class FaceQuadratureSpace : public QuadratureSpaceBase
 {
    FaceType face_type; ///< Is the space defined on interior or boundary faces?
-   const int num_faces; ///< Number of faces.
 
    /// Map from boundary or interior face indices to mesh face indices.
-   Array<int> face_indices;
+   const Array<int> &face_indices;
 
    /// Inverse of the map @a face_indices.
-   std::unordered_map<int,int> face_indices_inv;
+   const std::unordered_map<int,int> &face_indices_inv;
 
    const Vector &GetGeometricFactorWeights() const override;
-   void ConstructOffsets();
-   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 num_faces; }
+   inline int GetNumFaces() const { return face_indices.Size(); }
 
    /// Returns the face type (boundary or interior).
    FaceType GetFaceType() const { return face_type; }

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/scan.hpp

@@ -0,0 +1,176 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_SCAN_HPP
+#define MFEM_SCAN_HPP
+
+#ifdef MFEM_USE_CUDA
+#include <cub/device/device_scan.cuh>
+#define MFEM_CUB_NAMESPACE cub
+#elif MFEM_USE_HIP
+#include <hipcub/device/device_scan.hpp>
+#define MFEM_CUB_NAMESPACE hipcub
+#endif
+
+#include <functional>
+#include <numeric>
+
+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,
+                   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))
+   {
+      size_t bytes = workspace.Size();
+      if (bytes)
+      {
+         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
+                       workspace.Write(), bytes, d_in, d_out, num_items);
+#if defined(MFEM_USE_CUDA)
+         if (err == cudaSuccess)
+         {
+            return;
+         }
+#elif defined(MFEM_USE_HIP)
+         if (err == hipSuccess)
+         {
+            return;
+         }
+#endif
+      }
+      // try allocating a larger buffer
+      bytes = 0;
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
+                        nullptr, bytes, d_in, d_out, num_items));
+      workspace.SetSize(bytes);
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
+                        workspace.Write(), bytes, d_in, d_out, num_items));
+      return;
+   }
+#endif
+   std::inclusive_scan(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,
+                   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))
+   {
+      size_t bytes = workspace.Size();
+      if (bytes)
+      {
+         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
+                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
+#if defined(MFEM_USE_CUDA)
+         if (err == cudaSuccess)
+         {
+            return;
+         }
+#elif defined(MFEM_USE_HIP)
+         if (err == hipSuccess)
+         {
+            return;
+         }
+#endif
+      }
+      // try allocating a larger buffer
+      bytes = 0;
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
+                        nullptr, bytes, d_in, d_out, scan_op, num_items));
+      workspace.SetSize(bytes);
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
+                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
+      return;
+   }
+#endif
+   std::inclusive_scan(d_in, d_in + num_items, d_out, scan_op);
+}
+
+/// 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, 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))
+   {
+      size_t bytes = workspace.Size();
+      if (bytes)
+      {
+         auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
+                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                       num_items);
+#if defined(MFEM_USE_CUDA)
+         if (err == cudaSuccess)
+         {
+            return;
+         }
+#elif defined(MFEM_USE_HIP)
+         if (err == hipSuccess)
+         {
+            return;
+         }
+#endif
+      }
+      // try allocating a larger buffer
+      bytes = 0;
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
+                        nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
+      workspace.SetSize(bytes);
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
+                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                        num_items));
+      return;
+   }
+#endif
+   std::exclusive_scan(d_in, d_in + num_items, d_out, init_value, scan_op);
+}
+
+/// 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, Array<char> &workspace)
+{
+   ExclusiveScan(use_dev, d_in, d_out, num_items, init_value, workspace,
+                 std::plus<> {});
+}
+
+} // namespace mfem
+
+#undef MFEM_CUB_NAMESPACE
+
+#endif

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

@@ -979,6 +979,46 @@ const Array<int>& Mesh::GetElementAttributes() const
    return elem_attrs_cache;
 }
 
+void Mesh::ComputeFaceInfo(FaceType ftype) const
+{
+   auto &fidcs = face_indices[static_cast<int>(ftype)];
+   auto &ifidcs = inv_face_indices[static_cast<int>(ftype)];
+   fidcs.SetSize(GetNFbyType(ftype));
+   fidcs.HostWrite();
+   ifidcs.reserve(fidcs.Size());
+   int f_idx = 0;
+   for (int i = 0; i < GetNumFacesWithGhost(); ++i)
+   {
+      const FaceInformation face = GetFaceInformation(i);
+      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(ftype))
+      {
+         continue;
+      }
+      fidcs[f_idx] = i;
+      ifidcs[i] = f_idx;
+      ++f_idx;
+   }
+}
+
+const Array<int> &Mesh::GetFaceIndices(FaceType ftype) const
+{
+   if (face_indices[static_cast<int>(ftype)].Size() == 0)
+   {
+      ComputeFaceInfo(ftype);
+   }
+   return face_indices[static_cast<int>(ftype)];
+}
+
+const std::unordered_map<int, int> &
+Mesh::GetInvFaceIndices(FaceType ftype) const
+{
+   if (inv_face_indices[static_cast<int>(ftype)].empty())
+   {
+      ComputeFaceInfo(ftype);
+   }
+   return inv_face_indices[static_cast<int>(ftype)];
+}
+
 void Mesh::DeleteGeometricFactors()
 {
    for (int i = 0; i < geom_factors.Size(); i++)
@@ -1866,6 +1906,11 @@ void Mesh::Destroy()
    bdr_face_attrs_cache.DeleteAll();
    attributes.DeleteAll();
    bdr_attributes.DeleteAll();
+
+   face_indices[0].DeleteAll();
+   face_indices[1].DeleteAll();
+   inv_face_indices[0] = std::unordered_map<int, int>();
+   inv_face_indices[1] = std::unordered_map<int, int>();
 }
 
 void Mesh::ResetLazyData()
@@ -8137,6 +8182,12 @@ void Mesh::GenerateFaces()
       FreeElement(f);
    }
 
+   // delete caches
+   face_indices[0].SetSize(0);
+   face_indices[1].SetSize(0);
+   inv_face_indices[0].clear();
+   inv_face_indices[1].clear();
+
    // (re)generate the interior faces and the info for them
    faces.SetSize(nfaces);
    faces_info.SetSize(nfaces);
@@ -10936,6 +10987,11 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
    // copy attribute caches
    mfem::Swap(elem_attrs_cache, other.elem_attrs_cache);
    mfem::Swap(bdr_face_attrs_cache, other.bdr_face_attrs_cache);
+
+   mfem::Swap(face_indices[0], other.face_indices[0]);
+   mfem::Swap(face_indices[1], other.face_indices[1]);
+   inv_face_indices[0].swap(other.inv_face_indices[0]);
+   inv_face_indices[1].swap(other.inv_face_indices[1]);
 }
 
 void Mesh::GetElementData(const Array<Element*> &elem_array, int geom,

mesh/mesh.hpp

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

tests/unit/CMakeLists.txt

@@ -22,6 +22,7 @@ set(UNIT_TESTS_SRCS
   dfem/test_mass.cpp
   general/test_array.cpp
   general/test_reduction.cpp
+  general/test_scan.cpp
   general/test_arrays_by_name.cpp
   general/test_error.cpp
   general/test_mem.cpp

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

@@ -0,0 +1,102 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include <algorithm>
+#include <limits>
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+// must be included after mfem.hpp
+#include "general/scan.hpp"
+
+using namespace mfem;
+
+TEST_CASE("Inclusive Scan", "[Scan],[GPU]")
+{
+   Array<char> workspace;
+   Array<int> a(10);
+
+   for (int use_dev = 0; use_dev < 2; ++use_dev)
+   {
+      CAPTURE(use_dev);
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i;
+      }
+      auto dptr = a.ReadWrite(use_dev);
+      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace);
+      a.HostRead();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         int expected = (i + 1) * i / 2;
+         CAPTURE(i);
+         REQUIRE(a[i] == expected);
+      }
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i + 1;
+      }
+      a.ReadWrite(use_dev);
+      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace, std::multiplies<> {});
+      a.HostRead();
+      int expected = 1;
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         expected *= i + 1;
+         CAPTURE(i);
+         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)
+   {
+      CAPTURE(use_dev);
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i;
+      }
+      auto dptr = a.ReadWrite(use_dev);
+      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace);
+      a.HostRead();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         int expected = (i + 1) * i / 2 - i + 5;
+         CAPTURE(i);
+         REQUIRE(a[i] == expected);
+      }
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i + 1;
+      }
+      a.ReadWrite(use_dev);
+      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace,
+                    std::multiplies<> {});
+      a.HostRead();
+      int expected = 5;
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         CAPTURE(i);
+         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]);
+   }
+}