Merge branch 'master' into dfem-lvector-interface

.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/dfem/doperator.hpp

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

fem/dfem/util.hpp

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

fem/fe_coll.cpp

@@ -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); }
 

tests/unit/CMakeLists.txt

@@ -20,6 +20,7 @@ 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_arrays_by_name.cpp
@@ -144,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()
 
@@ -175,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()
 
 #-----------------------------------------------------------
@@ -183,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()
 
 #-----------------------------------------------------------
@@ -200,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()
 
 #-----------------------------------------------------------
@@ -224,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})
@@ -244,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)
 
@@ -277,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()
 
 #-----------------------------------------------------------
@@ -303,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)
 
 #-----------------------------------------------------------
@@ -424,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)
 
 #-----------------------------------------------------------
@@ -434,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

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