Merge branch 'gpu_em' into gpu-maxwell

[BUG] Fixed type narrowing in ProjectBdrCoefficientNormal unit test

CHANGELOG

@@ -66,6 +66,8 @@ Linear and nonlinear solvers
 
 GPU computing
 -------------
+- Added device assembly support for 3D H(curl) VectorFEDomainLFIntegrator.
+
 - Added NVIDIA cuDSS library interface. Implementation examples have been
   added to ex1 and ex1p. See https://developer.nvidia.com/cudss for more
   details. Supported versions >= 0.6.0.

config/cmake/modules/FindGSLIB.cmake

@@ -18,19 +18,17 @@
 
 if (MFEM_FETCH_GSLIB OR MFEM_FETCH_TPLS)
   enable_language(C)
-  string(TOUPPER "${CMAKE_BUILD_TYPE}" BUILD_TYPE)
   set(GSLIB_FETCH_VERSION 1.0.9)
-  set(GSLIB_C_FLAGS ${CMAKE_C_FLAGS_${BUILD_TYPE}})
-  if (CMAKE_C_FLAGS)
-    set(GSLIB_C_FLAGS "${CMAKE_C_FLAGS} ${CMAKE_C_FLAGS_${BUILD_TYPE}}")
-  endif()
-  if (BUILD_SHARED_LIBS)
-    set(GSLIB_C_FLAGS "${GSLIB_C_FLAGS} -fPIC")
-  endif()
   add_library(GSLIB STATIC IMPORTED)
+  # set options (technically flags because GSLIB does not use cmake)
+  string(TOUPPER "${CMAKE_BUILD_TYPE}" BUILD_TYPE)
+  set(GSLIB_FLAGS "${CMAKE_C_FLAGS} ${CMAKE_C_FLAGS_${BUILD_TYPE}}")
+  if (BUILD_SHARED_LIBS)
+    set(GSLIB_FLAGS "${GSLIB_FLAGS} -fPIC")
+  endif()
   # define external project and create future include directory so it is present
   # to pass CMake checks at end of MFEM configuration step
-  message(STATUS "Will fetch GSLIB ${GSLIB_FETCH_VERSION} to be built with ${GSLIB_C_FLAGS}")
+  message(STATUS "Will fetch GSLIB ${GSLIB_FETCH_VERSION} to be built with ${GSLIB_FLAGS}")
   set(PREFIX ${CMAKE_BINARY_DIR}/fetch/gslib)
   include(ExternalProject)
   ExternalProject_Add(gslib
@@ -40,7 +38,7 @@ if (MFEM_FETCH_GSLIB OR MFEM_FETCH_TPLS)
     UPDATE_DISCONNECTED TRUE
     PREFIX ${PREFIX}
     CONFIGURE_COMMAND ""
-    BUILD_COMMAND cd ${PREFIX}/src/gslib && $(MAKE) clean && $(MAKE) DESTDIR=${PREFIX} MPI=$<BOOL:${MFEM_USE_MPI}> "CFLAGS= ${GSLIB_C_FLAGS}"
+    BUILD_COMMAND cd ${PREFIX}/src/gslib && $(MAKE) clean && $(MAKE) DESTDIR=${PREFIX} MPI=$<BOOL:${MFEM_USE_MPI}> "CFLAGS=${GSLIB_FLAGS}"
     INSTALL_COMMAND "")
   file(MAKE_DIRECTORY ${PREFIX}/include)
   # set imported library target properties

config/cmake/modules/FindHYPRE.cmake

@@ -44,6 +44,9 @@ if (MFEM_FETCH_HYPRE OR MFEM_FETCH_TPLS)
   # set options and associated dependencies
   set(HYPRE_CMAKE_OPTIONS "")
   list(APPEND HYPRE_CMAKE_OPTIONS -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE})
+  if (BUILD_SHARED_LIBS)
+    list(APPEND HYPRE_CMAKE_OPTIONS -DCMAKE_POSITION_INDEPENDENT_CODE:BOOL=ON)
+  endif()
   # collect all HYPRE_ENABLE variables and pass them to hypre, assuming they are BOOL.
   get_cmake_property(all_vars VARIABLES)
   foreach(var ${all_vars})
@@ -95,7 +98,6 @@ if (MFEM_FETCH_HYPRE OR MFEM_FETCH_TPLS)
     UPDATE_DISCONNECTED TRUE
     SOURCE_SUBDIR src
     PREFIX ${HYPRE_INSTALL}
-    BUILD_COMMAND ${CMAKE_COMMAND} --build . -- -j${CMAKE_BUILD_PARALLEL_LEVEL}
     CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${HYPRE_INSTALL} -DCMAKE_INSTALL_LIBDIR:PATH=lib ${HYPRE_CMAKE_OPTIONS})
   file(MAKE_DIRECTORY ${HYPRE_INSTALL}/include)
   # set imported library target properties

config/cmake/modules/FindMETIS.cmake

@@ -19,10 +19,18 @@
 #   - METIS_VERSION_5 (cache variable)
 
 if (MFEM_FETCH_METIS OR MFEM_FETCH_TPLS)
+  enable_language(C)
   set(METIS_FETCH_VERSION 4.0.3)
   add_library(METIS STATIC IMPORTED)
+  # set options (technically flags because METIS does not use cmake)
+  set(METIS_FLAGS "-Wno-implicit-int -Wno-incompatible-pointer-types")
+  string(TOUPPER "${CMAKE_BUILD_TYPE}" BUILD_TYPE)
+  set(METIS_FLAGS "${METIS_FLAGS} ${CMAKE_C_FLAGS} ${CMAKE_C_FLAGS_${BUILD_TYPE}}")
+  if (BUILD_SHARED_LIBS)
+    set(METIS_FLAGS "${METIS_FLAGS} -fPIC")
+  endif()
   # define external project
-  message(STATUS "Will fetch METIS ${METIS_FETCH_VERSION} to be built with default options")
+  message(STATUS "Will fetch METIS ${METIS_FETCH_VERSION} to be built with ${METIS_FLAGS}")
   set(PREFIX ${CMAKE_BINARY_DIR}/fetch/metis)
   include(ExternalProject)
   ExternalProject_Add(metis
@@ -32,7 +40,7 @@ if (MFEM_FETCH_METIS OR MFEM_FETCH_TPLS)
     UPDATE_DISCONNECTED TRUE
     PREFIX ${PREFIX}
     CONFIGURE_COMMAND tar -xzf ../metis/metis-${METIS_FETCH_VERSION}-mac.tgz --strip=1
-    BUILD_COMMAND $(MAKE) COPTIONS=-Wno-incompatible-pointer-types
+    BUILD_COMMAND $(MAKE) clean && $(MAKE) "OPTFLAGS=${METIS_FLAGS}"
     INSTALL_COMMAND mkdir -p ${PREFIX}/lib && cp libmetis.a ${PREFIX}/lib/)
   # set imported library target properties
   add_dependencies(METIS metis)

fem/bilinearform.cpp

@@ -1255,6 +1255,31 @@ void BilinearForm::Mult(const Vector &x, Vector &y) const
    }
 }
 
+void BilinearForm::AddMult(const Vector &x, Vector &y, const real_t a) const
+{
+   if (ext)
+   {
+      ext->AddMult(x, y, a);
+   }
+   else
+   {
+      mat->AddMult(x, y, a);
+   }
+}
+
+void BilinearForm::AddMultTranspose(const Vector &x, Vector &y,
+                                    const real_t a) const
+{
+   if (ext)
+   {
+      ext->AddMultTranspose(x, y, a);
+   }
+   else
+   {
+      mat->AddMultTranspose(x, y, a);
+   }
+}
+
 void BilinearForm::MultTranspose(const Vector & x, Vector & y) const
 {
    if (ext)

fem/bilinearform.hpp

@@ -307,8 +307,8 @@ public:
    { mat->Mult(x, y); mat_e->AddMult(x, y); }
 
    /// Add the matrix vector multiple to a vector:  $ y += a M x $
-   void AddMult(const Vector &x, Vector &y, const real_t a = 1.0) const override
-   { mat -> AddMult (x, y, a); }
+   void AddMult(const Vector &x, Vector &y,
+                const real_t a = 1.0) const override;
 
    /** @brief Add the original uneliminated matrix vector multiple to a vector.
        The original matrix is $ M + Me $ so we have:
@@ -318,8 +318,7 @@ public:
 
    /// Add the matrix transpose vector multiplication:  $ y += a M^T x $
    void AddMultTranspose(const Vector & x, Vector & y,
-                         const real_t a = 1.0) const override
-   { mat->AddMultTranspose(x, y, a); }
+                         const real_t a = 1.0) const override;
 
    /** @brief Add the original uneliminated matrix transpose vector
        multiple to a vector. The original matrix is $ M + M_e $

fem/bilininteg.hpp

@@ -1110,20 +1110,36 @@ public:
    MixedDotProductIntegrator(VectorCoefficient &vq)
       : MixedScalarVectorIntegrator(vq, true) {}
 
-   inline virtual bool VerifyFiniteElementTypes(
+   inline bool VerifyFiniteElementTypes(
       const FiniteElement & trial_fe,
-      const FiniteElement & test_fe) const
+      const FiniteElement & test_fe) const override
    {
       return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
    }
 
-   inline virtual const char * FiniteElementTypeFailureMessage() const
+   inline const char * FiniteElementTypeFailureMessage() const override
    {
       return "MixedDotProductIntegrator:  "
              "Trial space must be a vector field "
              "and the test space must be a scalar field";
    }
+
+   using BilinearFormIntegrator::AssemblePA;
+   void AssemblePA(const FiniteElementSpace &trial_fes,
+                   const FiniteElementSpace &test_fes) override;
+
+   void AddMultPA(const Vector&, Vector&) const override;
+   void AddMultTransposePA(const Vector&, Vector&) const override;
+
+private:
+   Vector pa_data;
+   const DofToQuad *mapsO = nullptr;      ///< Not owned. Trial open map.
+   const DofToQuad *mapsC = nullptr;      ///< Not owned. Trial closed map.
+   const DofToQuad *mapsTest = nullptr;   ///< Not owned. Scalar test map.
+   const GeometricFactors *geom = nullptr;///< Not owned.
+   int dim = 0, ne = 0, dofs1D = 0, dofs1Dtest = 0, quad1D = 0;
+   bool test_map_integral = false;
 };
 
 /** Class for integrating the bilinear form $a(u,v) := (-\vec{V} \cdot u, \nabla \cdot v)$ in 2D or
@@ -2995,11 +3011,10 @@ public:
     vector (diagonal matrix), or matrix), trial function $u$ is in $H(curl$ or
     $H(div)$, and test function $v$ is in $H(curl$, $H(div)$, or $v=(v_1,\dots,v_n)$, where
     $v_i$ are in $H^1$. */
-class VectorFEMassIntegrator: public BilinearFormIntegrator
+class VectorFEMassIntegrator : public BilinearFormIntegrator
 {
 private:
-   void Init(Coefficient *q, DiagonalMatrixCoefficient *dq, MatrixCoefficient *mq)
-   { Q = q; DQ = dq; MQ = mq; }
+   void Init(Coefficient *q, DiagonalMatrixCoefficient *dq, MatrixCoefficient *mq);
 
 #ifndef MFEM_THREAD_SAFE
    Vector shape;
@@ -3022,7 +3037,8 @@ protected:
    const DofToQuad *mapsOtest;     ///< Not owned. DOF-to-quad map, open.
    const DofToQuad *mapsCtest;     ///< Not owned. DOF-to-quad map, closed.
    const GeometricFactors *geom;   ///< Not owned
-   int dim, ne, nq, dofs1D, dofs1Dtest, quad1D, trial_fetype, test_fetype;
+   int dim, ne, nq, dofs1D, dofs1Dtest, quad1D;
+   FiniteElement::DerivType trial_fetype, test_fetype;
    bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient
 
 public:
@@ -3053,6 +3069,29 @@ public:
                    const bool add) override;
 
    const Coefficient *GetCoefficient() const { return Q; }
+
+   using ApplyKernelType =
+      void (*)(const int NE, bool symmetric, const bool scalar_coeff,
+               const Array<real_t> &trialBO, const Array<real_t> &trialBC,
+               const Array<real_t> &testBOt, const Array<real_t> &testBCt,
+               const Vector &pa_data, const Vector &x, Vector &y,
+               const int triald1d, const int testd1d, const int q1d);
+
+   /// parameters: trial_fetype, test_fetype, ndims, trial_d1d, test_d1d, q1d
+   MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType,
+                         (FiniteElement::DerivType, FiniteElement::DerivType,
+                          int, int, int, int));
+
+   struct Kernels { Kernels(); };
+
+   template <FiniteElement::DerivType TrialType,
+             FiniteElement::DerivType TestType, int DIM, int TRIAL_D1D,
+             int TEST_D1D, int Q1D>
+   static void AddSpecialization()
+   {
+      ApplyPAKernels::Specialization<TrialType, TestType, DIM, TRIAL_D1D,
+                     TEST_D1D, Q1D>::Add();
+   }
 };
 
 /** Integrator for $(Q \nabla \cdot u, v)$ where $u=(u_1,\cdots,u_n)$ and all $u_i$ are in the same
@@ -3906,7 +3945,7 @@ class DiscreteInterpolator : public BilinearFormIntegrator { };
 
 
 /** Class for constructing the gradient as a DiscreteLinearOperator from an
-    $H^1$-conforming space to an $H(curl$-conforming space. The range space can be
+    $H^1$-conforming space to an $H(curl)$-conforming space. The range space can be
     vector $L_2$ space as well. */
 class GradientInterpolator : public DiscreteInterpolator
 {
@@ -4015,12 +4054,48 @@ public:
     discrete curl matrix. */
 class CurlInterpolator : public DiscreteInterpolator
 {
+   int dim, ne;
+   // "dof" are the domain fespace dof counts
+   int ndof_o;
+   // "quads" are the range fespace dof counts
+   int nquad_o;
+
+   Vector pa_data;
+
 public:
+   CurlInterpolator();
+
    void AssembleElementMatrix2(const FiniteElement &dom_fe,
                                const FiniteElement &ran_fe,
                                ElementTransformation &Trans,
                                DenseMatrix &elmat) override
    { ran_fe.ProjectCurl(dom_fe, Trans, elmat); }
+
+   void AssemblePA(const FiniteElementSpace &dom_fes,
+                   const FiniteElementSpace &ran_fes) override;
+   void AssemblePA(const FiniteElementSpace &fes) override
+   {
+      AssemblePA(fes, fes);
+   }
+   void AddMultPA(const Vector &x, Vector &y) const override;
+   void AddMultTransposePA(const Vector &x, Vector &y) const override;
+
+   using ApplyKernelType = void (*)(const int ne, const int ndof_o,
+                                    const int nquad_o, const Vector &pa,
+                                    const Vector &x, Vector &y);
+
+   /// arguments: DIM, ndof_o, nquad_o
+   MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType, (int, int, int));
+   /// arguments: DIM, ndof_o, nquad_o
+   MFEM_REGISTER_KERNELS(ApplyTPAKernels, ApplyKernelType, (int, int, int));
+
+   template <int DIM, int NDOF_O, int NQUAD_O> static void AddSpecialization()
+   {
+      ApplyPAKernels::Specialization<DIM, NDOF_O, NQUAD_O>::Add();
+      ApplyTPAKernels::Specialization<DIM, NDOF_O, NQUAD_O>::Add();
+   }
+
+   struct Kernels { Kernels(); };
 };
 
 

fem/coefficient.cpp

@@ -54,6 +54,8 @@ void Coefficient::Project(QuadratureFunction &qf)
    QuadratureSpaceBase &qspace = *qf.GetSpace();
    const int ne = qspace.GetNE();
    Vector values;
+   // GetValues makes a reference, but we need it to be valid on Host
+   qf.HostWrite();
    for (int iel = 0; iel < ne; ++iel)
    {
       qf.GetValues(iel, values);
@@ -327,6 +329,8 @@ void VectorCoefficient::Project(QuadratureFunction &qf)
    const int ne = qspace.GetNE();
    DenseMatrix values;
    Vector col;
+   // GetValues makes a reference, but we need it to be valid on Host
+   qf.HostWrite();
    for (int iel = 0; iel < ne; ++iel)
    {
       qf.GetValues(iel, values);
@@ -695,6 +699,8 @@ void MatrixCoefficient::Project(QuadratureFunction &qf, bool transpose)
    QuadratureSpaceBase &qspace = *qf.GetSpace();
    const int ne = qspace.GetNE();
    DenseMatrix values, matrix;
+   // GetValues makes a reference, but we need it to be valid on Host
+   qf.HostWrite();
    for (int iel = 0; iel < ne; ++iel)
    {
       qf.GetValues(iel, values);

fem/fe/fe_base.hpp

@@ -1446,6 +1446,8 @@ public:
                                 dof2quad_array_open);
    }
 
+   const Poly_1D::Basis &GetOpenBasis1D() const { return obasis1d; }
+
    virtual ~VectorTensorFiniteElement();
 };
 

fem/integ/bilininteg_diffusion_kernels.hpp

@@ -19,6 +19,8 @@
 #include "../../linalg/vector.hpp"
 #include "../bilininteg.hpp"
 
+#include "bilininteg_diffusion_pa_simplices.hpp"
+
 namespace mfem
 {
 

fem/integ/bilininteg_hcurl_kernels.cpp

@@ -147,18 +147,15 @@ void PAHcurlMassAssembleDiagonal3D(const int D1D,
    }); // end of element loop
 }
 
-void PAHcurlMassApply2D(const int D1D,
-                        const int Q1D,
-                        const int NE,
-                        const bool symmetric,
-                        const Array<real_t> &bo,
-                        const Array<real_t> &bc,
-                        const Array<real_t> &bot,
-                        const Array<real_t> &bct,
-                        const Vector &pa_data,
-                        const Vector &x,
-                        Vector &y)
+void PAHcurlMassApply2D(const int NE, const bool symmetric,
+                        const bool scalar_coeff, const Array<real_t> &bo,
+                        const Array<real_t> &bc, const Array<real_t> &bot,
+                        const Array<real_t> &bct, const Vector &pa_data,
+                        const Vector &x, Vector &y, const int D1D,
+                        const int TestD1D, const int Q1D)
 {
+   MFEM_ASSERT(D1D == TestD1D,
+               "Trial and Test space must have the same number of dofs");
    auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
    auto Bc = Reshape(bc.Read(), Q1D, D1D);
    auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
@@ -277,18 +274,15 @@ void PAHcurlMassApply2D(const int D1D,
    }); // end of element loop
 }
 
-void PAHcurlMassApply3D(const int D1D,
-                        const int Q1D,
-                        const int NE,
-                        const bool symmetric,
-                        const Array<real_t> &bo,
-                        const Array<real_t> &bc,
-                        const Array<real_t> &bot,
-                        const Array<real_t> &bct,
-                        const Vector &pa_data,
-                        const Vector &x,
-                        Vector &y)
+void PAHcurlMassApply3D(const int NE, const bool symmetric,
+                        const bool scalar_coeff, const Array<real_t> &bo,
+                        const Array<real_t> &bc, const Array<real_t> &bot,
+                        const Array<real_t> &bct, const Vector &pa_data,
+                        const Vector &x, Vector &y, const int D1D,
+                        const int TestD1D, const int Q1D)
 {
+   MFEM_VERIFY(D1D == TestD1D,
+               "Trial and test spaces must have same number of dofs");
    MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
                "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,

fem/integ/bilininteg_hcurl_kernels.hpp

@@ -181,228 +181,309 @@ inline void SmemPAHcurlMassAssembleDiagonal3D(const int d1d,
 }
 
 // PA H(curl) Mass Apply 2D kernel
-void PAHcurlMassApply2D(const int D1D,
-                        const int Q1D,
-                        const int NE,
-                        const bool symmetric,
-                        const Array<real_t> &bo,
-                        const Array<real_t> &bc,
-                        const Array<real_t> &bot,
-                        const Array<real_t> &bct,
-                        const Vector &pa_data,
-                        const Vector &x,
-                        Vector &y);
+void PAHcurlMassApply2D(const int NE, const bool symmetric,
+                        const bool scalar_coeff, const Array<real_t> &bo,
+                        const Array<real_t> &bc, const Array<real_t> &bot,
+                        const Array<real_t> &bct, const Vector &pa_data,
+                        const Vector &x, Vector &y, const int TrialD1D,
+                        const int TestD1D, const int Q1D);
 
 // PA H(curl) Mass Apply 3D kernel
-void PAHcurlMassApply3D(const int D1D,
-                        const int Q1D,
-                        const int NE,
-                        const bool symmetric,
-                        const Array<real_t> &bo,
-                        const Array<real_t> &bc,
-                        const Array<real_t> &bot,
-                        const Array<real_t> &bct,
-                        const Vector &pa_data,
-                        const Vector &x,
-                        Vector &y);
+void PAHcurlMassApply3D(const int NE, const bool symmetric,
+                        const bool scalar_coeff, const Array<real_t> &bo,
+                        const Array<real_t> &bc, const Array<real_t> &bot,
+                        const Array<real_t> &bct, const Vector &pa_data,
+                        const Vector &x, Vector &y, const int TrialD1D,
+                        const int TestD1D, const int Q1D);
 
 // Shared memory PA H(curl) Mass Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
-inline void SmemPAHcurlMassApply3D(const int d1d,
-                                   const int q1d,
-                                   const int NE,
-                                   const bool symmetric,
-                                   const Array<real_t> &bo,
-                                   const Array<real_t> &bc,
-                                   const Array<real_t> &bot,
-                                   const Array<real_t> &bct,
-                                   const Vector &pa_data,
-                                   const Vector &x,
-                                   Vector &y)
+template <int T_D1D = 0, int T_Q1D = 0, int TBATCH = 0, bool ACCUMULATE = true>
+inline void SmemPAHcurlMassApply3D(
+   const int NE, const bool symmetric, const bool scalar_coeff,
+   const Array<real_t> &bo, const Array<real_t> &bc, const Array<real_t> &bot,
+   const Array<real_t> &bct, const Vector &pa_data, const Vector &x, Vector &y,
+   const int d1d = 0, const int = 0, const int q1d = 0)
 {
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+
    MFEM_VERIFY(T_D1D || d1d <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
                "Error: d1d > HCURL_MAX_D1D");
    MFEM_VERIFY(T_Q1D || q1d <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,
                "Error: q1d > HCURL_MAX_Q1D");
-   const int D1D = T_D1D ? T_D1D : d1d;
-   const int Q1D = T_Q1D ? T_Q1D : q1d;
-
+   MFEM_ASSERT(Q1D >= D1D, "Expected Q1D >= D1D");
    const int dataSize = symmetric ? 6 : 9;
 
-   auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
-   auto Bc = Reshape(bc.Read(), Q1D, D1D);
-   auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, dataSize, NE);
-   auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
-   auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
+   // assume trial space == test space
+   auto Bo = bo.Read();
+   auto Bc = bc.Read();
+   auto op =
+      Reshape(pa_data.Read(), Q1D, Q1D, Q1D, dataSize, NE);
+   auto X_ = Reshape(x.Read(), 3 * (D1D - 1) * D1D * D1D, NE);
+   auto y_ = y.ReadWrite();
 
-   mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+   constexpr int MD_ = T_D1D ? T_D1D : DofQuadLimits::HCURL_MAX_D1D;
+   constexpr int MQ_ = T_Q1D ? T_Q1D : DofQuadLimits::HCURL_MAX_Q1D;
+   constexpr int MDQ_ = std::max(MD_, MQ_);
+   constexpr int MB_ = TBATCH ? TBATCH : 1;
+
+   mfem::forall_2D_batch<MDQ_ * MDQ_ * MDQ_ * MB_>(
+      NE, MDQ_ * MDQ_ * MDQ_, 1, MB_, [=] MFEM_HOST_DEVICE(int e)
    {
+#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
+      constexpr int nbz = TBATCH ? TBATCH : 1;
+      int tidz = MFEM_THREAD_ID(z);
+#else
+      constexpr int nbz = 1;
+      constexpr int tidz = 0;
+#endif
+
       constexpr int VDIM = 3;
       constexpr int MD1D = T_D1D ? T_D1D : DofQuadLimits::HCURL_MAX_D1D;
       constexpr int MQ1D = T_Q1D ? T_Q1D : DofQuadLimits::HCURL_MAX_Q1D;
-      const int D1D = T_D1D ? T_D1D : d1d;
-      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      constexpr int MDQ = std::max(MD1D, MQ1D);
 
-      MFEM_SHARED real_t sBo[MQ1D][MD1D];
-      MFEM_SHARED real_t sBc[MQ1D][MD1D];
+      // nvcc limit work-around: can't have Y_ be captured first in
+      // if constexpr, so capture y_ and construct Y_ locally
+      // only works on GPU
+      auto Y = Reshape(y_, VDIM * (D1D - 1) * D1D * D1D, NE);
 
-      real_t op9[9];
-      MFEM_SHARED real_t sop[9*MQ1D*MQ1D];
-      MFEM_SHARED real_t mass[MQ1D][MQ1D][3];
+      MFEM_SHARED real_t sBo[MDQ * (MD1D - 1)];
+      MFEM_SHARED real_t sBc[MDQ * MD1D];
+      auto BO = Reshape(sBo, Q1D, D1D - 1);
+      auto BC = Reshape(sBc, Q1D, D1D);
 
-      MFEM_SHARED real_t sX[MD1D][MD1D][MD1D];
+      MFEM_SHARED real_t sX[nbz * VDIM * (MD1D - 1) * MD1D * MD1D];
+      MFEM_SHARED real_t sm0[nbz * VDIM * MDQ * MDQ * MDQ];
+      MFEM_SHARED real_t sm1[nbz * VDIM * MDQ * MDQ * MDQ];
 
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      real_t(*X)[nbz][(MD1D - 1) * MD1D * MD1D] =
+         (real_t(*)[nbz][(MD1D - 1) * MD1D * MD1D])(sX);
+      // shapes of buffers always use MQ1D to mitigate shared memory bank
+      // conflicts
+      real_t(*DDQ)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm0);
+      real_t(*DQQ)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm1);
+      real_t(*QQQ)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm0);
+      real_t(*QQD)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm1);
+      real_t(*QDD)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm0);
+
+      // load dofs into smem
+      const int offset = (D1D - 1) * D1D * D1D;
+      MFEM_FOREACH_THREAD_DIRECT(ix, x, offset)
       {
-         MFEM_FOREACH_THREAD(qy,y,Q1D)
+         for (int dim = 0; dim < VDIM; ++dim)
          {
-            MFEM_FOREACH_THREAD(qz,z,Q1D)
-            {
-               for (int i=0; i<dataSize; ++i)
-               {
-                  op9[i] = op(qx,qy,qz,i,e);
-               }
-            }
+            X[dim][tidz][ix] = X_(ix + dim * offset, e);
          }
       }
-
-      const int tidx = MFEM_THREAD_ID(x);
-      const int tidy = MFEM_THREAD_ID(y);
-      const int tidz = MFEM_THREAD_ID(z);
-
+      // load basis functions data
       if (tidz == 0)
       {
-         MFEM_FOREACH_THREAD(d,y,D1D)
+         MFEM_FOREACH_THREAD_DIRECT(ix, x, D1D * Q1D) { sBc[ix] = Bc[ix]; }
+         MFEM_FOREACH_THREAD_DIRECT(ix, x, (D1D - 1) * Q1D)
          {
-            MFEM_FOREACH_THREAD(q,x,Q1D)
+            sBo[ix] = Bo[ix];
+         }
+      }
+
+      for (int dim0 = 0; dim0 < VDIM; ++dim0)
+      {
+         MFEM_SYNC_THREAD;
+         // sum factor to QQQ = Q_{dim0,dim1} B X_{dim1}
+         for (int dim1 = 0; dim1 < VDIM; ++dim1)
+         {
+            const int D1Dz = (dim1 == 2) ? D1D - 1 : D1D;
+            const int D1Dy = (dim1 == 1) ? D1D - 1 : D1D;
+            const int D1Dx = (dim1 == 0) ? D1D - 1 : D1D;
+
+            // threads assigned to mitigate bank conflicts
+            MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, Q1D, D1Dy, D1Dz,
+                                                 Q1D, Q1D, Q1D)
             {
-               sBc[q][d] = Bc(q,d);
-               if (d < D1D-1)
+               real_t u = 0;
+               for (int dx = 0; dx < D1Dx; ++dx)
                {
-                  sBo[q][d] = Bo(q,d);
+                  real_t b;
+                  if (dim1 == 0)
+                  {
+                     b = BO(qx, dx);
+                  }
+                  else
+                  {
+                     b = BC(qx, dx);
+                  }
+                  u += X[dim1][tidz][dx + (dy + dz * D1Dy) * D1Dx] * b;
+               }
+               DDQ[dim1][tidz][dz][dy][qx] = u;
+            }
+         }
+         MFEM_SYNC_THREAD;
+         for (int dim1 = 0; dim1 < VDIM; ++dim1)
+         {
+            const int D1Dz = (dim1 == 2) ? D1D - 1 : D1D;
+            const int D1Dy = (dim1 == 1) ? D1D - 1 : D1D;
+            // const int D1Dx = (dim1 == 0) ? D1D - 1 : D1D;
+            // threads assigned to mitigate bank conflicts
+            MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, Q1D, Q1D, D1Dz,
+                                                 Q1D, Q1D, Q1D)
+            {
+               real_t u = 0;
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  real_t b;
+                  if (dim1 == 1)
+                  {
+                     b = BO(qy, dy);
+                  }
+                  else
+                  {
+                     b = BC(qy, dy);
+                  }
+                  u += DDQ[dim1][tidz][dz][dy][qx] * b;
+               }
+               DQQ[dim1][tidz][dz][qy][qx] = u;
+            }
+         }
+         MFEM_SYNC_THREAD;
+         for (int dim1 = 0; dim1 < VDIM; ++dim1)
+         {
+            const int D1Dz = (dim1 == 2) ? D1D - 1 : D1D;
+            // const int D1Dy = (dim1 == 1) ? D1D - 1 : D1D;
+            // const int D1Dx = (dim1 == 0) ? D1D - 1 : D1D;
+            MFEM_FOREACH_THREAD_DIRECT_3D(qx, qy, qz, x, Q1D, Q1D, Q1D)
+            {
+               real_t u = 0;
+               for (int dz = 0; dz < D1Dz; ++dz)
+               {
+                  real_t b;
+                  if (dim1 == 2)
+                  {
+                     b = BO(qz, dz);
+                  }
+                  else
+                  {
+                     b = BC(qz, dz);
+                  }
+                  u += DQQ[dim1][tidz][dz][qy][qx] * b;
+               }
+               // pa_data is row major
+               int idx;
+               if (symmetric)
+               {
+                  int row;
+                  int col;
+                  if (dim0 > dim1)
+                  {
+                     row = dim1;
+                     col = dim0;
+                  }
+                  else
+                  {
+                     row = dim0;
+                     col = dim1;
+                  }
+                  idx = col + VDIM * row - row * (row + 1) / 2;
+               }
+               else
+               {
+                  idx = dim0 * VDIM + dim1;
+               }
+               QQQ[dim1][tidz][qz][qy][qx] = op(qx, qy, qz, idx, e) * u;
+            }
+         }
+         MFEM_SYNC_THREAD;
+         // sum factor back to Y
+         // Assume bot and bct == bo^t and bc^t respectively (i.e. test ==
+         // trial functions), skip loading them again.
+         {
+            const int D1Dz = (dim0 == 2) ? D1D - 1 : D1D;
+            const int D1Dy = (dim0 == 1) ? D1D - 1 : D1D;
+            const int D1Dx = (dim0 == 0) ? D1D - 1 : D1D;
+            // threads assigned to mitigate bank conflicts
+            MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, D1Dz, Q1D, Q1D,
+                                                 Q1D, Q1D, Q1D)
+            {
+               for (int dim1 = 0; dim1 < VDIM; ++dim1)
+               {
+                  real_t u = 0;
+                  for (int qz = 0; qz < Q1D; ++qz)
+                  {
+                     real_t b = 0;
+                     if (dim0 == 2)
+                     {
+                        b = BO(qz, dz);
+                     }
+                     else
+                     {
+                        b = BC(qz, dz);
+                     }
+                     u += QQQ[dim1][tidz][qz][qy][qx] * b;
+                  }
+                  QQD[dim1][tidz][qy][qx][dz] = u;
+               }
+            }
+            MFEM_SYNC_THREAD;
+            // threads assigned to mitigate bank conflicts
+            MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, D1Dy, D1Dz, Q1D,
+                                                 Q1D, Q1D, Q1D)
+            {
+               for (int dim1 = 0; dim1 < VDIM; ++dim1)
+               {
+                  real_t u = 0;
+                  for (int qy = 0; qy < Q1D; ++qy)
+                  {
+                     real_t b;
+                     if (dim0 == 1)
+                     {
+                        b = BO(qy, dy);
+                     }
+                     else
+                     {
+                        b = BC(qy, dy);
+                     }
+                     u += QQD[dim1][tidz][qy][qx][dz] * b;
+                  }
+                  QDD[dim1][tidz][qx][dz][dy] = u;
+               }
+            }
+            MFEM_SYNC_THREAD;
+
+            MFEM_FOREACH_THREAD_DIRECT_3D(dx, dy, dz, x, D1Dx, D1Dy, D1Dz)
+            {
+               int ix = dx + D1Dx * (dy + D1Dy * dz);
+               real_t u = 0;
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  real_t b;
+                  if (dim0 == 0)
+                  {
+                     b = BO(qx, dx);
+                  }
+                  else
+                  {
+                     b = BC(qx, dx);
+                  }
+                  for (int dim1 = 0; dim1 < VDIM; ++dim1)
+                  {
+                     u += QDD[dim1][tidz][qx][dz][dy] * b;
+                  }
+               }
+               if constexpr (ACCUMULATE)
+               {
+                  Y(ix + dim0 * offset, e) += u;
+               }
+               else
+               {
+                  Y(ix + dim0 * offset, e) = u;
                }
             }
          }
       }
-      MFEM_SYNC_THREAD;
-
-      for (int qz=0; qz < Q1D; ++qz)
-      {
-         int osc = 0;
-         for (int c = 0; c < VDIM; ++c)  // loop over x, y, z components
-         {
-            const int D1Dz = (c == 2) ? D1D - 1 : D1D;
-            const int D1Dy = (c == 1) ? D1D - 1 : D1D;
-            const int D1Dx = (c == 0) ? D1D - 1 : D1D;
-
-            MFEM_FOREACH_THREAD(dz,z,D1Dz)
-            {
-               MFEM_FOREACH_THREAD(dy,y,D1Dy)
-               {
-                  MFEM_FOREACH_THREAD(dx,x,D1Dx)
-                  {
-                     sX[dz][dy][dx] = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
-                  }
-               }
-            }
-            MFEM_SYNC_THREAD;
-
-            if (tidz == qz)
-            {
-               for (int i=0; i<dataSize; ++i)
-               {
-                  sop[i + (dataSize*tidx) + (dataSize*Q1D*tidy)] = op9[i];
-               }
-
-               MFEM_FOREACH_THREAD(qy,y,Q1D)
-               {
-                  MFEM_FOREACH_THREAD(qx,x,Q1D)
-                  {
-                     real_t u = 0.0;
-
-                     for (int dz = 0; dz < D1Dz; ++dz)
-                     {
-                        const real_t wz = (c == 2) ? sBo[qz][dz] : sBc[qz][dz];
-                        for (int dy = 0; dy < D1Dy; ++dy)
-                        {
-                           const real_t wy = (c == 1) ? sBo[qy][dy] : sBc[qy][dy];
-                           for (int dx = 0; dx < D1Dx; ++dx)
-                           {
-                              const real_t t = sX[dz][dy][dx];
-                              const real_t wx = (c == 0) ? sBo[qx][dx] : sBc[qx][dx];
-                              u += t * wx * wy * wz;
-                           }
-                        }
-                     }
-
-                     mass[qy][qx][c] = u;
-                  } // qx
-               } // qy
-            } // tidz == qz
-
-            osc += D1Dx * D1Dy * D1Dz;
-            MFEM_SYNC_THREAD;
-         } // c
-
-         MFEM_SYNC_THREAD;  // Sync mass[qy][qx][d] and sop
-
-         osc = 0;
-         for (int c = 0; c < VDIM; ++c)  // loop over x, y, z components
-         {
-            const int D1Dz = (c == 2) ? D1D - 1 : D1D;
-            const int D1Dy = (c == 1) ? D1D - 1 : D1D;
-            const int D1Dx = (c == 0) ? D1D - 1 : D1D;
-
-            real_t dxyz = 0.0;
-
-            MFEM_FOREACH_THREAD(dz,z,D1Dz)
-            {
-               const real_t wz = (c == 2) ? sBo[qz][dz] : sBc[qz][dz];
-
-               MFEM_FOREACH_THREAD(dy,y,D1Dy)
-               {
-                  MFEM_FOREACH_THREAD(dx,x,D1Dx)
-                  {
-                     for (int qy = 0; qy < Q1D; ++qy)
-                     {
-                        const real_t wy = (c == 1) ? sBo[qy][dy] : sBc[qy][dy];
-                        for (int qx = 0; qx < Q1D; ++qx)
-                        {
-                           const int os = (dataSize*qx) + (dataSize*Q1D*qy);
-                           const int id1 = os + ((c == 0) ? 0 : ((c == 1) ? (symmetric ? 1 : 3) :
-                                                                 (symmetric ? 2 : 6))); // O11, O21, O31
-                           const int id2 = os + ((c == 0) ? 1 : ((c == 1) ? (symmetric ? 3 : 4) :
-                                                                 (symmetric ? 4 : 7))); // O12, O22, O32
-                           const int id3 = os + ((c == 0) ? 2 : ((c == 1) ? (symmetric ? 4 : 5) :
-                                                                 (symmetric ? 5 : 8))); // O13, O23, O33
-
-                           const real_t m_c = (sop[id1] * mass[qy][qx][0]) + (sop[id2] * mass[qy][qx][1]) +
-                                              (sop[id3] * mass[qy][qx][2]);
-
-                           const real_t wx = (c == 0) ? sBo[qx][dx] : sBc[qx][dx];
-                           dxyz += m_c * wx * wy * wz;
-                        }
-                     }
-                  }
-               }
-            }
-
-            MFEM_SYNC_THREAD;
-
-            MFEM_FOREACH_THREAD(dz,z,D1Dz)
-            {
-               MFEM_FOREACH_THREAD(dy,y,D1Dy)
-               {
-                  MFEM_FOREACH_THREAD(dx,x,D1Dx)
-                  {
-                     Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += dxyz;
-                  }
-               }
-            }
-
-            osc += D1Dx * D1Dy * D1Dz;
-         } // c loop
-      } // qz
    }); // end of element loop
 }
 

fem/integ/bilininteg_hcurlhdiv_kernels.hpp

@@ -62,6 +62,30 @@ void PAHcurlHdivMassApply2D(const int D1D,
                             const Vector &x_,
                             Vector &y_);
 
+/// H(curl) test, H(div) trial
+inline void
+PAHcurlHdivMassApply2D(const int NE, const bool, const bool scalarCoeff,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int D1Dtest, const int Q1D)
+{
+   return PAHcurlHdivMassApply2D(D1D, D1Dtest, Q1D, NE, scalarCoeff, false,
+                                 false, Bo_, Bc_, Bot_, Bct_, op_, x_, y_);
+}
+
+/// H(div) test, H(curl) trial
+inline void
+PAHdivHcurlMassApply2D(const int NE, const bool, const bool scalarCoeff,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int D1Dtest, const int Q1D)
+{
+   return PAHcurlHdivMassApply2D(D1D, D1Dtest, Q1D, NE, scalarCoeff, true,
+                                 false, Bo_, Bc_, Bot_, Bct_, op_, x_, y_);
+}
+
 // PA H(curl)-H(div) Mass Apply 3D kernel
 void PAHcurlHdivMassApply3D(const int D1D,
                             const int D1Dtest,
@@ -78,6 +102,30 @@ void PAHcurlHdivMassApply3D(const int D1D,
                             const Vector &x_,
                             Vector &y_);
 
+/// H(curl) test, H(div) trial
+inline void
+PAHcurlHdivMassApply3D(const int NE, const bool, const bool scalarCoeff,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int D1Dtest, const int Q1D)
+{
+   PAHcurlHdivMassApply3D(D1D, D1Dtest, Q1D, NE, scalarCoeff, false, false, Bo_,
+                          Bc_, Bot_, Bct_, op_, x_, y_);
+}
+
+/// H(div) test, H(curl) trial
+inline void
+PAHdivHcurlMassApply3D(const int NE, const bool, const bool scalarCoeff,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int D1Dtest, const int Q1D)
+{
+   PAHcurlHdivMassApply3D(D1D, D1Dtest, Q1D, NE, scalarCoeff, true, false, Bo_,
+                          Bc_, Bot_, Bct_, op_, x_, y_);
+}
+
 // PA H(curl)-H(div) Curl Apply 3D kernel
 template<int T_D1D = 0, int T_D1D_TEST = 0, int T_Q1D = 0>
 inline void PAHcurlHdivApply3D(const int d1d,
@@ -816,8 +864,656 @@ inline void PAHcurlHdivApplyTranspose3D(const int d1d,
    }); // end of element loop
 }
 
+namespace curlinterp
+{
+constexpr int NBZ3D(int ndof_o, int nquad_o, int mdq)
+{
+   if (ndof_o <= 0 || nquad_o <= 0)
+   {
+      return 1;
+   }
+   int ndof_c = ndof_o + 1;
+   int nquad_c = nquad_o + 1;
+   // z dimension is capped at 64 on nvidia and amd gpus
+   int tmp =
+      std::min((128 + mdq * mdq * (mdq - 1) - 1) / (mdq * mdq * (mdq - 1)), 64);
+   int smem_req =
+      sizeof(mfem::real_t) *
+      ((3 * ndof_c * ndof_c * ndof_o + 2 * 2 * mdq * mdq * mdq) * tmp +
+       ndof_c * nquad_o + ndof_c * nquad_c + ndof_o * nquad_o);
+   // assume GPU has at least 48k shared memory
+   return std::max(std::min(tmp, (48 * 1024 + smem_req - 1) / smem_req), 1);
+}
+}
+
+template <int T_NDOF_O, int T_NQUAD_O>
+void CurlInterpolatorApply3DSmem(const int ne, const int ndof_o,
+                                 const int nquad_o, const Vector &pa,
+                                 const Vector &x_, Vector &y_)
+{
+   constexpr int mnd_o = T_NDOF_O ? T_NDOF_O : DofQuadLimits::HCURL_MAX_D1D - 1;
+   constexpr int mnq_o =
+      T_NQUAD_O ? T_NQUAD_O : DofQuadLimits::HDIV_MAX_D1D - 1;
+   constexpr int mndq = std::max(mnd_o + 1, mnq_o + 1);
+   constexpr int tbatch = curlinterp::NBZ3D(T_NDOF_O, T_NQUAD_O, mndq);
+   MFEM_VERIFY(ndof_o <= mnd_o, "Error: H(curl) order larger than supported");
+   MFEM_VERIFY(nquad_o <= mnq_o, "Error: H(div) order larger than supported");
+   int mnq = std::max(ndof_o + 1, nquad_o + 1);
+   auto pa_data = pa.Read();
+   auto x_d = x_.Read();
+   auto y_d = y_.ReadWrite();
+   mfem::forall_2D_batch<mndq * mndq * (mndq - 1) * tbatch>(
+      ne, mnq * mnq * (mnq - 1), 1, tbatch, [=] MFEM_HOST_DEVICE(int e)
+   {
+      constexpr int MND_O =
+         T_NDOF_O ? T_NDOF_O : DofQuadLimits::HCURL_MAX_D1D - 1;
+      constexpr int MNQ_O =
+         T_NQUAD_O ? T_NQUAD_O : DofQuadLimits::HDIV_MAX_D1D - 1;
+      constexpr int MNDQ = std::max(MND_O + 1, MNQ_O + 1);
+#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
+      constexpr int nbz = curlinterp::NBZ3D(T_NDOF_O, T_NQUAD_O, MNDQ);
+      int tidz = MFEM_THREAD_ID(z);
+      // Make mnq a local variable since capturing would result in different
+      // captures between host/device versions, and spuriously fails
+      int mnq = std::max(ndof_o + 1, nquad_o + 1);
+#else
+      constexpr int nbz = 1;
+      constexpr int tidz = 0;
+#endif
+      const int NDOF_O = T_NDOF_O ? T_NDOF_O : ndof_o;
+      const int NQUAD_O = T_NQUAD_O ? T_NQUAD_O : nquad_o;
+      const int NDOF_C = NDOF_O + 1;
+      const int NQUAD_C = NQUAD_O + 1;
+      MFEM_SHARED real_t
+      sBG[(MND_O + 1) * MNQ_O + (MND_O + 1) * (MNQ_O + 1) + MND_O * MNQ_O];
+      auto X_ = Reshape(x_d, 3 * NDOF_C * NDOF_C * NDOF_O, ne);
+      auto Y = Reshape(y_d, 3 * NQUAD_C * NQUAD_O * NQUAD_O, ne);
+      auto Gco = Reshape(sBG, NQUAD_O, NDOF_C);
+      auto Bcc = Reshape(sBG + NDOF_C * NQUAD_O, NQUAD_C, NDOF_C);
+      auto Boo =
+         Reshape(sBG + NDOF_C * NQUAD_O + NDOF_C * NQUAD_C, NQUAD_O, NDOF_O);
+      MFEM_SHARED real_t X[3][nbz][MND_O * (MND_O + 1) * (MND_O + 1)];
+      MFEM_SHARED real_t sm0[nbz * 2 * MNDQ * MNDQ * MNDQ];
+      MFEM_SHARED real_t sm1[nbz * 2 * MNDQ * MNDQ * MNDQ];
+
+      // shapes of buffers always use MNDQ to mitigate shared memory bank
+      // conflicts
+      real_t(*DDQ)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm0);
+      real_t(*DQQ)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm1);
+      real_t(*QQQ)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm0);
+      const int offset = NDOF_O * NDOF_C * NDOF_C;
+      const int offsetq = NQUAD_C * NQUAD_O * NQUAD_O;
+      MFEM_FOREACH_THREAD_DIRECT(ix, x, offset)
+      {
+         for (int dim = 0; dim < 3; ++dim)
+         {
+            X[dim][tidz][ix] = X_(ix + dim * offset, e);
+         }
+      }
+      // load basis functions data
+      if (tidz == 0)
+      {
+         auto npts = NDOF_C * NQUAD_O + NDOF_C * NQUAD_C + NDOF_O * NQUAD_O;
+         MFEM_FOREACH_THREAD(ix, x, npts) { sBG[ix] = pa_data[ix]; }
+      }
+      MFEM_SYNC_THREAD;
+
+      // x: Vz Bcc Gco Boo - Vy Bcc Boo Gco
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_C, NDOF_C,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_C; ++dx)
+         {
+            u += X[2][tidz][dx + (dy + dz * NDOF_C) * NDOF_C] * Bcc(qx, dx);
+         }
+         DDQ[0][tidz][dz][dy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_C, NDOF_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_C; ++dx)
+         {
+            u += X[1][tidz][dx + (dy + dz * NDOF_O) * NDOF_C] * Bcc(qx, dx);
+         }
+         DDQ[1][tidz][dz][dy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_C, NQUAD_O,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_C; ++dy)
+         {
+            u += DDQ[0][tidz][dz][dy][qx] * Gco(qy, dy);
+         }
+         DQQ[0][tidz][dz][qy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_C, NQUAD_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_O; ++dy)
+         {
+            u += DDQ[1][tidz][dz][dy][qx] * Boo(qy, dy);
+         }
+         DQQ[1][tidz][dz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_C, NQUAD_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_O; ++dz)
+         {
+            u += DQQ[0][tidz][dz][qy][qx] * Boo(qz, dz);
+         }
+         QQQ[0][tidz][qz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_C, NQUAD_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_C; ++dz)
+         {
+            u += DQQ[1][tidz][dz][qy][qx] * Gco(qz, dz);
+         }
+         Y(qx + (qy + qz * NQUAD_O) * NQUAD_C, e) =
+            QQQ[0][tidz][qz][qy][qx] - u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // y: Vx Boo Bcc Gco - Vz Gco Bcc Boo
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_O, NDOF_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_O; ++dx)
+         {
+            u += X[0][tidz][dx + (dy + dz * NDOF_C) * NDOF_O] * Boo(qx, dx);
+         }
+         DDQ[0][tidz][dz][dy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_O, NDOF_C,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_C; ++dx)
+         {
+            u += X[2][tidz][dx + (dy + dz * NDOF_C) * NDOF_C] * Gco(qx, dx);
+         }
+         DDQ[1][tidz][dz][dy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_O, NQUAD_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_C; ++dy)
+         {
+            u += DDQ[0][tidz][dz][dy][qx] * Bcc(qy, dy);
+         }
+         DQQ[0][tidz][dz][qy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_O, NQUAD_C,
+                                           NDOF_O, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_C; ++dy)
+         {
+            u += DDQ[1][tidz][dz][dy][qx] * Bcc(qy, dy);
+         }
+         DQQ[1][tidz][dz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_O, NQUAD_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_C; ++dz)
+         {
+            u += DQQ[0][tidz][dz][qy][qx] * Gco(qz, dz);
+         }
+         QQQ[0][tidz][qz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_O, NQUAD_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_O; ++dz)
+         {
+            u += DQQ[1][tidz][dz][qy][qx] * Boo(qz, dz);
+         }
+         Y(qx + (qy + qz * NQUAD_C) * NQUAD_O + offsetq, e) =
+            QQQ[0][tidz][qz][qy][qx] - u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // z: Vy Gco Boo Bcc - Vx Boo Gco Bcc
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_O, NDOF_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_C; ++dx)
+         {
+            u += X[1][tidz][dx + (dy + dz * NDOF_O) * NDOF_C] * Gco(qx, dx);
+         }
+         DDQ[0][tidz][dz][dy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_O, NDOF_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_O; ++dx)
+         {
+            u += X[0][tidz][dx + (dy + dz * NDOF_C) * NDOF_O] * Boo(qx, dx);
+         }
+         DDQ[1][tidz][dz][dy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_O, NQUAD_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_O; ++dy)
+         {
+            u += DDQ[0][tidz][dz][dy][qx] * Boo(qy, dy);
+         }
+         DQQ[0][tidz][dz][qy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_O, NQUAD_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_C; ++dy)
+         {
+            u += DDQ[1][tidz][dz][dy][qx] * Gco(qy, dy);
+         }
+         DQQ[1][tidz][dz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_O, NQUAD_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_C; ++dz)
+         {
+            u += DQQ[0][tidz][dz][qy][qx] * Bcc(qz, dz);
+         }
+         QQQ[0][tidz][qz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_O, NQUAD_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_C; ++dz)
+         {
+            u += DQQ[1][tidz][dz][qy][qx] * Bcc(qz, dz);
+         }
+         Y(qx + (qy + qz * NQUAD_O) * NQUAD_O + 2 * offsetq, e) =
+            QQQ[0][tidz][qz][qy][qx] - u;
+      }
+      MFEM_SYNC_THREAD;
+   });
+}
+
+template <int T_NDOF_O, int T_NQUAD_O>
+void CurlInterpolatorTApply3DSmem(const int ne, const int ndof_o,
+                                  const int nquad_o, const Vector &pa,
+                                  const Vector &x_, Vector &y_)
+{
+   constexpr int mnd_o = T_NDOF_O ? T_NDOF_O : DofQuadLimits::HCURL_MAX_D1D - 1;
+   constexpr int mnq_o =
+      T_NQUAD_O ? T_NQUAD_O : DofQuadLimits::HDIV_MAX_D1D - 1;
+   constexpr int mndq = std::max(mnd_o + 1, mnq_o + 1);
+   constexpr int tbatch = curlinterp::NBZ3D(T_NDOF_O, T_NQUAD_O, mndq);
+   MFEM_VERIFY(ndof_o <= mnd_o, "Error: H(curl) order larger than supported");
+   MFEM_VERIFY(nquad_o <= mnq_o, "Error: H(div) order larger than supported");
+   int mnq = std::max(ndof_o + 1, nquad_o + 1);
+   auto pa_data = pa.Read();
+   auto x_d = x_.Read();
+   auto y_d = y_.ReadWrite();
+   mfem::forall_2D_batch<mndq * mndq * (mndq - 1) * tbatch>(
+      ne, mnq * mnq * (mnq - 1), 1, tbatch, [=] MFEM_HOST_DEVICE(int e)
+   {
+      constexpr int MND_O =
+         T_NDOF_O ? T_NDOF_O : DofQuadLimits::HCURL_MAX_D1D - 1;
+      constexpr int MNQ_O =
+         T_NQUAD_O ? T_NQUAD_O : DofQuadLimits::HDIV_MAX_D1D - 1;
+      constexpr int MNDQ = std::max(MND_O + 1, MNQ_O + 1);
+#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
+      constexpr int nbz = curlinterp::NBZ3D(T_NDOF_O, T_NQUAD_O, MNDQ);
+      int tidz = MFEM_THREAD_ID(z);
+      // Make mnq a local variable since capturing would result in different
+      // captures between host/device versions, and spuriously fails
+      int mnq = std::max(ndof_o + 1, nquad_o + 1);
+#else
+      constexpr int nbz = 1;
+      constexpr int tidz = 0;
+#endif
+      const int NDOF_O = T_NDOF_O ? T_NDOF_O : ndof_o;
+      const int NQUAD_O = T_NQUAD_O ? T_NQUAD_O : nquad_o;
+      const int NDOF_C = NDOF_O + 1;
+      const int NQUAD_C = NQUAD_O + 1;
+      MFEM_SHARED real_t
+      sBG[(MND_O + 1) * MNQ_O + (MND_O + 1) * (MNQ_O + 1) + MND_O * MNQ_O];
+      auto X_ = Reshape(x_d, 3 * NQUAD_C * NQUAD_O * NQUAD_O, ne);
+      auto Y = Reshape(y_d, 3 * NDOF_C * NDOF_C * NDOF_O, ne);
+      auto Gco = Reshape(sBG, NQUAD_O, NDOF_C);
+      auto Bcc = Reshape(sBG + NDOF_C * NQUAD_O, NQUAD_C, NDOF_C);
+      auto Boo =
+         Reshape(sBG + NDOF_C * NQUAD_O + NDOF_C * NQUAD_C, NQUAD_O, NDOF_O);
+      MFEM_SHARED real_t X[3][nbz][MNQ_O * MNQ_O * (MNQ_O + 1)];
+      MFEM_SHARED real_t sm0[nbz * 2 * MNDQ * MNDQ * MNDQ];
+      MFEM_SHARED real_t sm1[nbz * 2 * MNDQ * MNDQ * MNDQ];
+
+      // shapes of buffers always use MNDQ to mitigate shared memory bank
+      // conflicts
+      real_t(*QQD)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm0);
+      real_t(*QDD)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm1);
+      real_t(*DDD)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm0);
+      const int offset = NDOF_O * NDOF_C * NDOF_C;
+      const int offsetq = NQUAD_C * NQUAD_O * NQUAD_O;
+      MFEM_FOREACH_THREAD_DIRECT(ix, x, offsetq)
+      {
+         for (int dim = 0; dim < 3; ++dim)
+         {
+            X[dim][tidz][ix] = X_(ix + dim * offsetq, e);
+         }
+      }
+      // load basis functions data
+      if (tidz == 0)
+      {
+         auto npts = NDOF_C * NQUAD_O + NDOF_C * NQUAD_C + NDOF_O * NQUAD_O;
+         MFEM_FOREACH_THREAD(ix, x, npts) { sBG[ix] = pa_data[ix]; }
+      }
+      MFEM_SYNC_THREAD;
+
+      // x: Vy Boo Bcc Gco - Vz Boo Gco Bcc
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_C, NQUAD_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_O; ++qz)
+         {
+            u += X[1][tidz][qx + (qy + qz * NQUAD_C) * NQUAD_O] * Gco(qz, dz);
+         }
+         QQD[0][tidz][qy][qx][dz] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_C, NQUAD_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_C; ++qz)
+         {
+            u += X[2][tidz][qx + (qy + qz * NQUAD_O) * NQUAD_O] * Bcc(qz, dz);
+         }
+         QQD[1][tidz][qy][qx][dz] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_C, NDOF_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_C; ++qy)
+         {
+            u += QQD[0][tidz][qy][qx][dz] * Bcc(qy, dy);
+         }
+         QDD[0][tidz][qx][dz][dy] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_C, NDOF_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_O; ++qy)
+         {
+            u += QQD[1][tidz][qy][qx][dz] * Gco(qy, dy);
+         }
+         QDD[1][tidz][qx][dz][dy] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_O, NDOF_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_O; ++qx)
+         {
+            u += QDD[0][tidz][qx][dz][dy] * Boo(qx, dx);
+         }
+         DDD[0][tidz][dz][dy][dx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_O, NDOF_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_O; ++qx)
+         {
+            u += QDD[1][tidz][qx][dz][dy] * Boo(qx, dx);
+         }
+         Y(dx + (dy + dz * NDOF_C) * NDOF_O, e) =
+            DDD[0][tidz][dz][dy][dx] - u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // y: Vz Gco Boo Bcc - Vx Bcc Boo Gco
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_C, NQUAD_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_C; ++qz)
+         {
+            u += X[2][tidz][qx + (qy + qz * NQUAD_O) * NQUAD_O] * Bcc(qz, dz);
+         }
+         QQD[0][tidz][qy][qx][dz] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_C, NQUAD_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_O; ++qz)
+         {
+            u += X[0][tidz][qx + (qy + qz * NQUAD_O) * NQUAD_C] * Gco(qz, dz);
+         }
+         QQD[1][tidz][qy][qx][dz] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_O, NDOF_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_O; ++qy)
+         {
+            u += QQD[0][tidz][qy][qx][dz] * Boo(qy, dy);
+         }
+         QDD[0][tidz][qx][dz][dy] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_O, NDOF_C,
+                                           NQUAD_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_O; ++qy)
+         {
+            u += QQD[1][tidz][qy][qx][dz] * Boo(qy, dy);
+         }
+         QDD[1][tidz][qx][dz][dy] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_C, NDOF_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_O; ++qx)
+         {
+            u += QDD[0][tidz][qx][dz][dy] * Gco(qx, dx);
+         }
+         DDD[0][tidz][dz][dy][dx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_C, NDOF_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_C; ++qx)
+         {
+            u += QDD[1][tidz][qx][dz][dy] * Bcc(qx, dx);
+         }
+         Y(dx + (dy + dz * NDOF_O) * NDOF_C + offset, e) =
+            DDD[0][tidz][dz][dy][dx] - u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // z: Vx Bcc Gco Boo - Vy Gco Bcc Boo
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_O, NQUAD_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_O; ++qz)
+         {
+            u += X[0][tidz][qx + (qy + qz * NQUAD_O) * NQUAD_C] * Boo(qz, dz);
+         }
+         QQD[0][tidz][qy][qx][dz] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_O, NQUAD_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_O; ++qz)
+         {
+            u += X[1][tidz][qx + (qy + qz * NQUAD_C) * NQUAD_O] * Boo(qz, dz);
+         }
+         QQD[1][tidz][qy][qx][dz] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_C, NDOF_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_O; ++qy)
+         {
+            u += QQD[0][tidz][qy][qx][dz] * Gco(qy, dy);
+         }
+         QDD[0][tidz][qx][dz][dy] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_C, NDOF_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_C; ++qy)
+         {
+            u += QQD[1][tidz][qy][qx][dz] * Bcc(qy, dy);
+         }
+         QDD[1][tidz][qx][dz][dy] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_C, NDOF_C,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_C; ++qx)
+         {
+            u += QDD[0][tidz][qx][dz][dy] * Bcc(qx, dx);
+         }
+         DDD[0][tidz][dz][dy][dx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_C, NDOF_C,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_O; ++qx)
+         {
+            u += QDD[1][tidz][qx][dz][dy] * Gco(qx, dx);
+         }
+         Y(dx + (dy + dz * NDOF_C) * NDOF_C + 2 * offset, e) =
+            DDD[0][tidz][dz][dy][dx] - u;
+      }
+      MFEM_SYNC_THREAD;
+   });
+}
+
 } // namespace internal
 
+template <int DIM, int NDOF_O, int NQUAD_O>
+CurlInterpolator::ApplyKernelType
+CurlInterpolator::ApplyPAKernels::Kernel()
+{
+   if constexpr (DIM == 3)
+   {
+      return internal::CurlInterpolatorApply3DSmem<NDOF_O, NQUAD_O>;
+   }
+   MFEM_ABORT("Bad dimension!");
+}
+
+template <int DIM, int NDOF_O, int NQUAD_O>
+CurlInterpolator::ApplyKernelType
+CurlInterpolator::ApplyTPAKernels::Kernel()
+{
+   if constexpr (DIM == 3)
+   {
+      return internal::CurlInterpolatorTApply3DSmem<NDOF_O, NQUAD_O>;
+   }
+   MFEM_ABORT("Bad dimension!");
+}
+
 } // namespace mfem
 
 /// \endcond DO_NOT_DOCUMENT

fem/integ/bilininteg_hdiv_kernels.cpp

@@ -294,61 +294,14 @@ void PAHdivMassAssembleDiagonal3D(const int D1D,
    }); // end of element loop
 }
 
-void PAHdivMassApply(const int dim,
-                     const int D1D,
-                     const int Q1D,
-                     const int NE,
-                     const bool symmetric,
-                     const Array<real_t> &Bo,
-                     const Array<real_t> &Bc,
-                     const Array<real_t> &Bot,
-                     const Array<real_t> &Bct,
-                     const Vector &op,
-                     const Vector &x,
-                     Vector &y)
-{
-   const int id = (D1D << 4) | Q1D;
-
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x22: return SmemPAHdivMassApply2D<2,2>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x33: return SmemPAHdivMassApply2D<3,3>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x44: return SmemPAHdivMassApply2D<4,4>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x55: return SmemPAHdivMassApply2D<5,5>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         default: // fallback
-            return PAHdivMassApply2D(D1D,Q1D,NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-      }
-   }
-   else if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x23: return SmemPAHdivMassApply3D<2,3>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x34: return SmemPAHdivMassApply3D<3,4>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x45: return SmemPAHdivMassApply3D<4,5>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x56: return SmemPAHdivMassApply3D<5,6>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x67: return SmemPAHdivMassApply3D<6,7>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x78: return SmemPAHdivMassApply3D<7,8>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         default: // fallback
-            return PAHdivMassApply3D(D1D,Q1D,NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-      }
-   }
-}
-
-void PAHdivMassApply2D(const int D1D,
-                       const int Q1D,
-                       const int NE,
-                       const bool symmetric,
-                       const Array<real_t> &Bo_,
-                       const Array<real_t> &Bc_,
-                       const Array<real_t> &Bot_,
-                       const Array<real_t> &Bct_,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_)
+void PAHdivMassApply2D(const int NE, const bool symmetric, const bool,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int TestD1D, const int Q1D)
 {
+   MFEM_VERIFY(D1D == TestD1D,
+               "Trial and test spaces must have same number of dofs");
    auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
    auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
    auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
@@ -468,18 +421,14 @@ void PAHdivMassApply2D(const int D1D,
    }); // end of element loop
 }
 
-void PAHdivMassApply3D(const int D1D,
-                       const int Q1D,
-                       const int NE,
-                       const bool symmetric,
-                       const Array<real_t> &Bo_,
-                       const Array<real_t> &Bc_,
-                       const Array<real_t> &Bot_,
-                       const Array<real_t> &Bct_,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_)
+void PAHdivMassApply3D(const int NE, const bool symmetric, const bool,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int TestD1D, const int Q1D)
 {
+   MFEM_VERIFY(D1D == TestD1D,
+               "Trial and test spaces must have same number of dofs");
    MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
                "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,

fem/integ/bilininteg_hdiv_kernels.hpp

@@ -66,58 +66,29 @@ void PAHdivMassAssembleDiagonal3D(const int D1D,
                                   const Vector &op_,
                                   Vector &diag_);
 
-void PAHdivMassApply(const int dim,
-                     const int D1D,
-                     const int Q1D,
-                     const int NE,
-                     const bool symmetric,
-                     const Array<real_t> &Bo,
-                     const Array<real_t> &Bc,
-                     const Array<real_t> &Bot,
-                     const Array<real_t> &Bct,
-                     const Vector &op,
-                     const Vector &x,
-                     Vector &y);
-
 // PA H(div) Mass Apply 2D kernel
-void PAHdivMassApply2D(const int D1D,
-                       const int Q1D,
-                       const int NE,
-                       const bool symmetric,
-                       const Array<real_t> &Bo_,
-                       const Array<real_t> &Bc_,
-                       const Array<real_t> &Bot_,
-                       const Array<real_t> &Bct_,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_);
+void PAHdivMassApply2D(const int NE, const bool symmetric,
+                       const bool scalar_coeff, const Array<real_t> &Bo_,
+                       const Array<real_t> &Bc_, const Array<real_t> &Bot_,
+                       const Array<real_t> &Bct_, const Vector &op_,
+                       const Vector &x_, Vector &y_, const int D1D,
+                       const int TestD1D, const int Q1D);
 
 // PA H(div) Mass Apply 3D kernel
-void PAHdivMassApply3D(const int D1D,
-                       const int Q1D,
-                       const int NE,
-                       const bool symmetric,
-                       const Array<real_t> &Bo_,
-                       const Array<real_t> &Bc_,
-                       const Array<real_t> &Bot_,
-                       const Array<real_t> &Bct_,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_);
+void PAHdivMassApply3D(const int NE, const bool symmetric,
+                       const bool scalar_coeff, const Array<real_t> &Bo_,
+                       const Array<real_t> &Bc_, const Array<real_t> &Bot_,
+                       const Array<real_t> &Bct_, const Vector &op_,
+                       const Vector &x_, Vector &y_, const int D1D,
+                       const int TestD1D, const int Q1D);
 
 // Shared memory PA H(div) Mass Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
-inline void SmemPAHdivMassApply2D(const int NE,
-                                  const bool symmetric,
-                                  const Array<real_t> &Bo_,
-                                  const Array<real_t> &Bc_,
-                                  const Array<real_t> &Bot_,
-                                  const Array<real_t> &Bct_,
-                                  const Vector &op_,
-                                  const Vector &x_,
-                                  Vector &y_,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+template <int T_D1D = 0, int T_Q1D = 0>
+inline void SmemPAHdivMassApply2D(
+   const int NE, const bool symmetric, const bool, const Array<real_t> &Bo_,
+   const Array<real_t> &Bc_, const Array<real_t> &Bot_,
+   const Array<real_t> &Bct_, const Vector &op_, const Vector &x_, Vector &y_,
+   const int d1d = 0, const int = 0, const int q1d = 0)
 {
    MFEM_CONTRACT_VAR(Bot_);
    MFEM_CONTRACT_VAR(Bct_);
@@ -280,18 +251,13 @@ inline void SmemPAHdivMassApply2D(const int NE,
 }
 
 // Shared memory PA H(div) Mass Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
-inline void SmemPAHdivMassApply3D(const int NE,
-                                  const bool symmetric,
-                                  const Array<real_t> &Bo_,
-                                  const Array<real_t> &Bc_,
-                                  const Array<real_t> &Bot_,
-                                  const Array<real_t> &Bct_,
-                                  const Vector &op_,
-                                  const Vector &x_,
-                                  Vector &y_,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+template <int T_D1D = 0, int T_Q1D = 0>
+inline void
+SmemPAHdivMassApply3D(const int NE, const bool symmetric, const bool,
+                      const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                      const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                      const Vector &op_, const Vector &x_, Vector &y_,
+                      const int d1d = 0, const int = 0, const int q1d = 0)
 {
    MFEM_CONTRACT_VAR(Bot_);
    MFEM_CONTRACT_VAR(Bct_);

fem/integ/bilininteg_interp_pa.cpp

@@ -14,6 +14,8 @@
 #include "../gridfunc.hpp"
 #include "../qfunction.hpp"
 
+#include "bilininteg_hcurlhdiv_kernels.hpp"
+
 namespace mfem
 {
 
@@ -1950,4 +1952,138 @@ void IdentityInterpolator::AddMultTransposePA(const Vector &x, Vector &y) const
    }
 }
 
+void CurlInterpolator::AssemblePA(const FiniteElementSpace &dom_fes,
+                                  const FiniteElementSpace &ran_fes)
+{
+   // TODO: 1D and 2D meshes
+   Mesh *mesh = dom_fes.GetMesh();
+   const VectorTensorFiniteElement *dom_el =
+      dynamic_cast<const VectorTensorFiniteElement *>(dom_fes.GetTypicalFE());
+   const VectorTensorFiniteElement *ran_el =
+      dynamic_cast<const VectorTensorFiniteElement *>(ran_fes.GetTypicalFE());
+   MFEM_VERIFY(dom_el != NULL, "Only VectorTensorFiniteElement is supported!");
+   MFEM_VERIFY(ran_el != NULL, "Only VectorTensorFiniteElement is supported!");
+   // only supports H(curl) -> H(div) because of discontinuity requirements
+   MFEM_VERIFY(dom_el->GetDerivType() == FiniteElement::CURL,
+               "Domain space must be H(curl)");
+   MFEM_VERIFY(ran_el->GetDerivType() == FiniteElement::DIV,
+               "Range space must be H(div)");
+
+   const int dims = dom_el->GetDim();
+   MFEM_VERIFY(dims == 3, "");
+   dim = mesh->Dimension();
+
+   ne = dom_fes.GetNE();
+   ndof_o = dom_el->GetOrder();
+   int ndof_c = ndof_o + 1;
+   nquad_o = ran_el->GetOrder();
+   int nquad_c = nquad_o + 1;
+
+   // extract the tensor product range dof locations
+   std::vector<real_t> qc(nquad_c);
+   std::vector<real_t> qo(nquad_o);
+   {
+      const IntegrationRule &ran_nodes = ran_el->GetNodes();
+      const Array<int> &quad_map = ran_el->GetDofMap();
+      for (int i = 0; i < nquad_c; ++i)
+      {
+         int idx = UnsignIndex(quad_map[i]);
+         qc[i] = ran_nodes.IntPoint(idx).x;
+      }
+      int offset = ndof_c * ndof_o * ndof_o;
+      for (int i = 0; i < nquad_o; ++i)
+      {
+         int idx = UnsignIndex(quad_map[i + offset]);
+         qo[i] = ran_nodes.IntPoint(idx).x;
+      }
+   }
+
+   // evaluate closed/open 1D basis (and their derivatives) at closed and
+   // open quads
+   // storage order: GCO, BCC, BOO
+   pa_data.SetSize(ndof_c * nquad_o + ndof_c * nquad_c + ndof_o * nquad_o);
+   auto ptr = pa_data.HostWrite();
+   auto &cbasis1d = dom_el->GetBasis1D();
+   auto &obasis1d = dom_el->GetOpenBasis1D();
+   Vector b, g;
+   b.SetSize(ndof_c);
+   g.SetSize(ndof_c);
+   for (int j = 0; j < nquad_o; ++j)
+   {
+      cbasis1d.Eval(qo[j], b, g);
+      for (int i = 0; i < ndof_c; ++i)
+      {
+         ptr[j + i * nquad_o] = g[i];
+      }
+   }
+   ptr += nquad_o * ndof_c;
+
+   for (int j = 0; j < nquad_c; ++j)
+   {
+      cbasis1d.Eval(qc[j], b);
+      for (int i = 0; i < ndof_c; ++i)
+      {
+         ptr[j + i * nquad_c] = b[i];
+      }
+   }
+   ptr += ndof_c * nquad_c;
+
+   b.SetSize(ndof_o);
+   for (int j = 0; j < nquad_o; ++j)
+   {
+      obasis1d.Eval(qo[j], b);
+      for (int i = 0; i < ndof_o; ++i)
+      {
+         ptr[j + i * nquad_o] = b[i];
+      }
+   }
+}
+
+CurlInterpolator::Kernels::Kernels()
+{
+   CurlInterpolator::AddSpecialization<3, 1, 1>();
+   CurlInterpolator::AddSpecialization<3, 2, 2>();
+   CurlInterpolator::AddSpecialization<3, 3, 3>();
+   CurlInterpolator::AddSpecialization<3, 4, 4>();
+   CurlInterpolator::AddSpecialization<3, 5, 5>();
+}
+
+CurlInterpolator::CurlInterpolator() { static Kernels kernels{}; }
+
+void CurlInterpolator::AddMultPA(const Vector &x, Vector &y) const
+{
+   ApplyPAKernels::Run(dim, ndof_o, nquad_o, ne, ndof_o, nquad_o, pa_data, x,
+                       y);
+}
+
+void CurlInterpolator::AddMultTransposePA(const Vector &x, Vector &y) const
+{
+   ApplyTPAKernels::Run(dim, ndof_o, nquad_o, ne, ndof_o, nquad_o, pa_data, x,
+                        y);
+}
+
+/// \cond DO_NOT_DOCUMENT
+
+CurlInterpolator::ApplyKernelType
+CurlInterpolator::ApplyPAKernels::Fallback(int DIM, int, int)
+{
+   if (DIM == 3)
+   {
+      return internal::CurlInterpolatorApply3DSmem<0, 0>;
+   }
+   MFEM_ABORT("Bad dimension!");
+}
+
+CurlInterpolator::ApplyKernelType
+CurlInterpolator::ApplyTPAKernels::Fallback(int DIM, int, int)
+{
+   if (DIM == 3)
+   {
+      return internal::CurlInterpolatorTApply3DSmem<0, 0>;
+   }
+   MFEM_ABORT("Bad dimension!");
+}
+
+/// \endcond DO_NOT_DOCUMENT
+
 } // namespace mfem

fem/integ/bilininteg_mass_kernels.hpp

@@ -19,6 +19,8 @@
 #include "../../linalg/vector.hpp"
 #include "../bilininteg.hpp"
 
+#include "bilininteg_mass_pa_simplices.hpp"
+
 namespace mfem
 {
 

fem/integ/bilininteg_mixedcurl_pa.cpp

@@ -18,6 +18,779 @@
 namespace mfem
 {
 
+namespace
+{
+
+void PAHcurlDotSetup2D(const int q1d,
+                       const int ne,
+                       const bool test_map_integral,
+                       const Array<real_t> &w,
+                       const Vector &jacobians,
+                       const Vector &coeff,
+                       Vector &op)
+{
+   auto W = Reshape(w.Read(), q1d, q1d);
+   auto J = Reshape(jacobians.Read(), q1d, q1d, 2, 2, ne);
+   auto C = Reshape(coeff.Read(), 2, q1d, q1d, ne);
+   auto O = Reshape(op.Write(), 2, q1d, q1d, ne);
+
+   mfem::forall_2D(ne, q1d, q1d, [=] MFEM_HOST_DEVICE (int e)
+   {
+      MFEM_FOREACH_THREAD(qy, y, q1d)
+      {
+         MFEM_FOREACH_THREAD(qx, x, q1d)
+         {
+            const real_t J11 = J(qx, qy, 0, 0, e);
+            const real_t J12 = J(qx, qy, 1, 0, e);
+            const real_t J21 = J(qx, qy, 0, 1, e);
+            const real_t J22 = J(qx, qy, 1, 1, e);
+            const real_t detJ = (J11 * J22) - (J21 * J12);
+            const real_t scale = W(qx, qy) * (test_map_integral ? 1.0 / detJ : 1.0);
+            const real_t Vx = C(0, qx, qy, e);
+            const real_t Vy = C(1, qx, qy, e);
+
+            O(0, qx, qy, e) = scale * ( J22 * Vx - J12 * Vy);
+            O(1, qx, qy, e) = scale * (-J21 * Vx + J11 * Vy);
+         }
+      }
+   });
+}
+
+void PAHcurlDotSetup3D(const int q1d,
+                       const int ne,
+                       const bool test_map_integral,
+                       const Array<real_t> &w,
+                       const Vector &jacobians,
+                       const Vector &coeff,
+                       Vector &op)
+{
+   auto W = Reshape(w.Read(), q1d, q1d, q1d);
+   auto J = Reshape(jacobians.Read(), q1d, q1d, q1d, 3, 3, ne);
+   auto C = Reshape(coeff.Read(), 3, q1d, q1d, q1d, ne);
+   auto O = Reshape(op.Write(), 3, q1d, q1d, q1d, ne);
+
+   mfem::forall_3D(ne, q1d, q1d, q1d, [=] MFEM_HOST_DEVICE (int e)
+   {
+      MFEM_FOREACH_THREAD(qz, z, q1d)
+      {
+         MFEM_FOREACH_THREAD(qy, y, q1d)
+         {
+            MFEM_FOREACH_THREAD(qx, x, q1d)
+            {
+               const real_t J11 = J(qx, qy, qz, 0, 0, e);
+               const real_t J12 = J(qx, qy, qz, 0, 1, e);
+               const real_t J13 = J(qx, qy, qz, 0, 2, e);
+               const real_t J21 = J(qx, qy, qz, 1, 0, e);
+               const real_t J22 = J(qx, qy, qz, 1, 1, e);
+               const real_t J23 = J(qx, qy, qz, 1, 2, e);
+               const real_t J31 = J(qx, qy, qz, 2, 0, e);
+               const real_t J32 = J(qx, qy, qz, 2, 1, e);
+               const real_t J33 = J(qx, qy, qz, 2, 2, e);
+               const real_t detJ = J11 * (J22 * J33 - J32 * J23)
+                                   - J21 * (J12 * J33 - J32 * J13)
+                                   + J31 * (J12 * J23 - J22 * J13);
+               const real_t scale = W(qx, qy, qz) *
+                                    (test_map_integral ? 1.0 / detJ : 1.0);
+               const real_t Vx = C(0, qx, qy, qz, e);
+               const real_t Vy = C(1, qx, qy, qz, e);
+               const real_t Vz = C(2, qx, qy, qz, e);
+
+               O(0, qx, qy, qz, e) = scale *
+                                     ((J22 * J33 - J23 * J32) * Vx +
+                                      (J13 * J32 - J12 * J33) * Vy +
+                                      (J12 * J23 - J13 * J22) * Vz);
+               O(1, qx, qy, qz, e) = scale *
+                                     ((J23 * J31 - J21 * J33) * Vx +
+                                      (J11 * J33 - J13 * J31) * Vy +
+                                      (J13 * J21 - J11 * J23) * Vz);
+               O(2, qx, qy, qz, e) = scale *
+                                     ((J21 * J32 - J22 * J31) * Vx +
+                                      (J12 * J31 - J11 * J32) * Vy +
+                                      (J11 * J22 - J12 * J21) * Vz);
+            }
+         }
+      }
+   });
+}
+
+void PAHcurlDotApply2D(const int d1d,
+                       const int d1d_test,
+                       const int q1d,
+                       const int ne,
+                       const Array<real_t> &bo,
+                       const Array<real_t> &bc,
+                       const Array<real_t> &bt,
+                       const Vector &pa_data,
+                       const Vector &x,
+                       Vector &y)
+{
+   MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D, "");
+   MFEM_VERIFY(d1d_test <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D, "");
+
+   auto Bo = Reshape(bo.Read(), q1d, d1d - 1);
+   auto Bc = Reshape(bc.Read(), q1d, d1d);
+   auto Bt = Reshape(bt.Read(), d1d_test, q1d);
+   auto O = Reshape(pa_data.Read(), 2, q1d, q1d, ne);
+   auto X = Reshape(x.Read(), 2 * (d1d - 1) * d1d, ne);
+   auto Y = Reshape(y.ReadWrite(), d1d_test, d1d_test, ne);
+
+   mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
+   {
+      constexpr int MAX_D1D = DofQuadLimits::MAX_D1D;
+      constexpr int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
+
+      real_t u0[MAX_Q1D][MAX_Q1D];
+      real_t u1[MAX_Q1D][MAX_Q1D];
+
+      for (int qy = 0; qy < q1d; ++qy)
+      {
+         for (int qx = 0; qx < q1d; ++qx)
+         {
+            u0[qy][qx] = 0.0;
+            u1[qy][qx] = 0.0;
+         }
+      }
+
+      int osc = 0;
+      for (int dy = 0; dy < d1d; ++dy)
+      {
+         real_t mass_x[MAX_Q1D];
+         for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] = 0.0; }
+         for (int dx = 0; dx < d1d - 1; ++dx)
+         {
+            const real_t t = X(dx + (dy * (d1d - 1)) + osc, e);
+            for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] += t * Bo(qx, dx); }
+         }
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            const real_t wy = Bc(qy, dy);
+            for (int qx = 0; qx < q1d; ++qx) { u0[qy][qx] += mass_x[qx] * wy; }
+         }
+      }
+
+      osc += (d1d - 1) * d1d;
+      for (int dy = 0; dy < d1d - 1; ++dy)
+      {
+         real_t mass_x[MAX_Q1D];
+         for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] = 0.0; }
+         for (int dx = 0; dx < d1d; ++dx)
+         {
+            const real_t t = X(dx + (dy * d1d) + osc, e);
+            for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] += t * Bc(qx, dx); }
+         }
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            const real_t wy = Bo(qy, dy);
+            for (int qx = 0; qx < q1d; ++qx) { u1[qy][qx] += mass_x[qx] * wy; }
+         }
+      }
+
+      for (int qy = 0; qy < q1d; ++qy)
+      {
+         real_t sol_x[MAX_D1D];
+         for (int dx = 0; dx < d1d_test; ++dx) { sol_x[dx] = 0.0; }
+         for (int qx = 0; qx < q1d; ++qx)
+         {
+            const real_t s = O(0, qx, qy, e) * u0[qy][qx]
+                             + O(1, qx, qy, e) * u1[qy][qx];
+            for (int dx = 0; dx < d1d_test; ++dx)
+            {
+               sol_x[dx] += s * Bt(dx, qx);
+            }
+         }
+         for (int dy = 0; dy < d1d_test; ++dy)
+         {
+            const real_t wy = Bt(dy, qy);
+            for (int dx = 0; dx < d1d_test; ++dx)
+            {
+               Y(dx, dy, e) += sol_x[dx] * wy;
+            }
+         }
+      }
+   });
+}
+
+void PAHcurlDotApplyTranspose2D(const int d1d,
+                                const int d1d_test,
+                                const int q1d,
+                                const int ne,
+                                const Array<real_t> &bo,
+                                const Array<real_t> &bc,
+                                const Array<real_t> &b,
+                                const Vector &pa_data,
+                                const Vector &x,
+                                Vector &y)
+{
+   MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D, "");
+   MFEM_VERIFY(d1d_test <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D, "");
+
+   auto Bo = Reshape(bo.Read(), q1d, d1d - 1);
+   auto Bc = Reshape(bc.Read(), q1d, d1d);
+   auto B = Reshape(b.Read(), q1d, d1d_test);
+   auto O = Reshape(pa_data.Read(), 2, q1d, q1d, ne);
+   auto X = Reshape(x.Read(), d1d_test, d1d_test, ne);
+   auto Y = Reshape(y.ReadWrite(), 2 * (d1d - 1) * d1d, ne);
+
+   mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
+   {
+      constexpr int MAX_D1D = DofQuadLimits::MAX_D1D;
+      constexpr int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
+
+      real_t mass[MAX_Q1D][MAX_Q1D];
+      for (int qy = 0; qy < q1d; ++qy)
+      {
+         for (int qx = 0; qx < q1d; ++qx)
+         {
+            mass[qy][qx] = 0.0;
+         }
+      }
+
+      for (int dy = 0; dy < d1d_test; ++dy)
+      {
+         real_t sol_x[MAX_Q1D];
+         for (int qx = 0; qx < q1d; ++qx) { sol_x[qx] = 0.0; }
+         for (int dx = 0; dx < d1d_test; ++dx)
+         {
+            const real_t t = X(dx, dy, e);
+            for (int qx = 0; qx < q1d; ++qx) { sol_x[qx] += t * B(qx, dx); }
+         }
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            const real_t wy = B(qy, dy);
+            for (int qx = 0; qx < q1d; ++qx) { mass[qy][qx] += sol_x[qx] * wy; }
+         }
+      }
+
+      int osc = 0;
+      for (int qy = 0; qy < q1d; ++qy)
+      {
+         real_t mass_x[MAX_D1D];
+         for (int dx = 0; dx < d1d - 1; ++dx) { mass_x[dx] = 0.0; }
+         for (int qx = 0; qx < q1d; ++qx)
+         {
+            const real_t s = O(0, qx, qy, e) * mass[qy][qx];
+            for (int dx = 0; dx < d1d - 1; ++dx) { mass_x[dx] += s * Bo(qx, dx); }
+         }
+         for (int dy = 0; dy < d1d; ++dy)
+         {
+            const real_t wy = Bc(qy, dy);
+            for (int dx = 0; dx < d1d - 1; ++dx)
+            {
+               Y(dx + (dy * (d1d - 1)) + osc, e) += mass_x[dx] * wy;
+            }
+         }
+      }
+
+      osc += (d1d - 1) * d1d;
+      for (int qy = 0; qy < q1d; ++qy)
+      {
+         real_t mass_x[MAX_D1D];
+         for (int dx = 0; dx < d1d; ++dx) { mass_x[dx] = 0.0; }
+         for (int qx = 0; qx < q1d; ++qx)
+         {
+            const real_t s = O(1, qx, qy, e) * mass[qy][qx];
+            for (int dx = 0; dx < d1d; ++dx) { mass_x[dx] += s * Bc(qx, dx); }
+         }
+         for (int dy = 0; dy < d1d - 1; ++dy)
+         {
+            const real_t wy = Bo(qy, dy);
+            for (int dx = 0; dx < d1d; ++dx)
+            {
+               Y(dx + (dy * d1d) + osc, e) += mass_x[dx] * wy;
+            }
+         }
+      }
+   });
+}
+
+void PAHcurlDotApply3D(const int d1d,
+                       const int d1d_test,
+                       const int q1d,
+                       const int ne,
+                       const Array<real_t> &bo,
+                       const Array<real_t> &bc,
+                       const Array<real_t> &bt,
+                       const Vector &pa_data,
+                       const Vector &x,
+                       Vector &y)
+{
+   MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D, "");
+   MFEM_VERIFY(d1d_test <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D, "");
+
+   auto Bo = Reshape(bo.Read(), q1d, d1d - 1);
+   auto Bc = Reshape(bc.Read(), q1d, d1d);
+   auto Bt = Reshape(bt.Read(), d1d_test, q1d);
+   auto O = Reshape(pa_data.Read(), 3, q1d, q1d, q1d, ne);
+   auto X = Reshape(x.Read(), 3 * (d1d - 1) * d1d * d1d, ne);
+   auto Y = Reshape(y.ReadWrite(), d1d_test, d1d_test, d1d_test, ne);
+
+   mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
+   {
+      constexpr int MAX_D1D = DofQuadLimits::MAX_D1D;
+      constexpr int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
+
+      real_t u[MAX_Q1D][MAX_Q1D][MAX_Q1D][3];
+      for (int qz = 0; qz < q1d; ++qz)
+      {
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            for (int qx = 0; qx < q1d; ++qx)
+            {
+               for (int c = 0; c < 3; ++c) { u[qz][qy][qx][c] = 0.0; }
+            }
+         }
+      }
+
+      int osc = 0;
+      for (int dz = 0; dz < d1d; ++dz)
+      {
+         real_t mass_xy[MAX_Q1D][MAX_Q1D];
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            for (int qx = 0; qx < q1d; ++qx) { mass_xy[qy][qx] = 0.0; }
+         }
+
+         for (int dy = 0; dy < d1d; ++dy)
+         {
+            real_t mass_x[MAX_Q1D];
+            for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] = 0.0; }
+            for (int dx = 0; dx < d1d - 1; ++dx)
+            {
+               const real_t t = X(dx + ((dy + (dz * d1d)) * (d1d - 1)) + osc, e);
+               for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] += t * Bo(qx, dx); }
+            }
+            for (int qy = 0; qy < q1d; ++qy)
+            {
+               const real_t wy = Bc(qy, dy);
+               for (int qx = 0; qx < q1d; ++qx) { mass_xy[qy][qx] += mass_x[qx] * wy; }
+            }
+         }
+
+         for (int qz = 0; qz < q1d; ++qz)
+         {
+            const real_t wz = Bc(qz, dz);
+            for (int qy = 0; qy < q1d; ++qy)
+            {
+               for (int qx = 0; qx < q1d; ++qx) { u[qz][qy][qx][0] += mass_xy[qy][qx] * wz; }
+            }
+         }
+      }
+
+      osc += (d1d - 1) * d1d * d1d;
+      for (int dz = 0; dz < d1d; ++dz)
+      {
+         real_t mass_xy[MAX_Q1D][MAX_Q1D];
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            for (int qx = 0; qx < q1d; ++qx) { mass_xy[qy][qx] = 0.0; }
+         }
+
+         for (int dy = 0; dy < d1d - 1; ++dy)
+         {
+            real_t mass_x[MAX_Q1D];
+            for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] = 0.0; }
+            for (int dx = 0; dx < d1d; ++dx)
+            {
+               const real_t t = X(dx + ((dy + (dz * (d1d - 1))) * d1d) + osc, e);
+               for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] += t * Bc(qx, dx); }
+            }
+            for (int qy = 0; qy < q1d; ++qy)
+            {
+               const real_t wy = Bo(qy, dy);
+               for (int qx = 0; qx < q1d; ++qx) { mass_xy[qy][qx] += mass_x[qx] * wy; }
+            }
+         }
+
+         for (int qz = 0; qz < q1d; ++qz)
+         {
+            const real_t wz = Bc(qz, dz);
+            for (int qy = 0; qy < q1d; ++qy)
+            {
+               for (int qx = 0; qx < q1d; ++qx) { u[qz][qy][qx][1] += mass_xy[qy][qx] * wz; }
+            }
+         }
+      }
+
+      osc += (d1d - 1) * d1d * d1d;
+      for (int dz = 0; dz < d1d - 1; ++dz)
+      {
+         real_t mass_xy[MAX_Q1D][MAX_Q1D];
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            for (int qx = 0; qx < q1d; ++qx) { mass_xy[qy][qx] = 0.0; }
+         }
+
+         for (int dy = 0; dy < d1d; ++dy)
+         {
+            real_t mass_x[MAX_Q1D];
+            for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] = 0.0; }
+            for (int dx = 0; dx < d1d; ++dx)
+            {
+               const real_t t = X(dx + ((dy + (dz * d1d)) * d1d) + osc, e);
+               for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] += t * Bc(qx, dx); }
+            }
+            for (int qy = 0; qy < q1d; ++qy)
+            {
+               const real_t wy = Bc(qy, dy);
+               for (int qx = 0; qx < q1d; ++qx) { mass_xy[qy][qx] += mass_x[qx] * wy; }
+            }
+         }
+
+         for (int qz = 0; qz < q1d; ++qz)
+         {
+            const real_t wz = Bo(qz, dz);
+            for (int qy = 0; qy < q1d; ++qy)
+            {
+               for (int qx = 0; qx < q1d; ++qx) { u[qz][qy][qx][2] += mass_xy[qy][qx] * wz; }
+            }
+         }
+      }
+
+      for (int qz = 0; qz < q1d; ++qz)
+      {
+         real_t mass_xy[MAX_D1D][MAX_D1D];
+         for (int dy = 0; dy < d1d_test; ++dy)
+         {
+            for (int dx = 0; dx < d1d_test; ++dx) { mass_xy[dy][dx] = 0.0; }
+         }
+
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            real_t mass_x[MAX_D1D];
+            for (int dx = 0; dx < d1d_test; ++dx) { mass_x[dx] = 0.0; }
+            for (int qx = 0; qx < q1d; ++qx)
+            {
+               const real_t s = O(0, qx, qy, qz, e) * u[qz][qy][qx][0]
+                                + O(1, qx, qy, qz, e) * u[qz][qy][qx][1]
+                                + O(2, qx, qy, qz, e) * u[qz][qy][qx][2];
+               for (int dx = 0; dx < d1d_test; ++dx) { mass_x[dx] += s * Bt(dx, qx); }
+            }
+            for (int dy = 0; dy < d1d_test; ++dy)
+            {
+               const real_t wy = Bt(dy, qy);
+               for (int dx = 0; dx < d1d_test; ++dx) { mass_xy[dy][dx] += mass_x[dx] * wy; }
+            }
+         }
+
+         for (int dz = 0; dz < d1d_test; ++dz)
+         {
+            const real_t wz = Bt(dz, qz);
+            for (int dy = 0; dy < d1d_test; ++dy)
+            {
+               for (int dx = 0; dx < d1d_test; ++dx)
+               {
+                  Y(dx, dy, dz, e) += mass_xy[dy][dx] * wz;
+               }
+            }
+         }
+      }
+   });
+}
+
+void PAHcurlDotApplyTranspose3D(const int d1d,
+                                const int d1d_test,
+                                const int q1d,
+                                const int ne,
+                                const Array<real_t> &bo,
+                                const Array<real_t> &bc,
+                                const Array<real_t> &b,
+                                const Vector &pa_data,
+                                const Vector &x,
+                                Vector &y)
+{
+   MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D, "");
+   MFEM_VERIFY(d1d_test <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D, "");
+
+   auto Bo = Reshape(bo.Read(), q1d, d1d - 1);
+   auto Bc = Reshape(bc.Read(), q1d, d1d);
+   auto B = Reshape(b.Read(), q1d, d1d_test);
+   auto O = Reshape(pa_data.Read(), 3, q1d, q1d, q1d, ne);
+   auto X = Reshape(x.Read(), d1d_test, d1d_test, d1d_test, ne);
+   auto Y = Reshape(y.ReadWrite(), 3 * (d1d - 1) * d1d * d1d, ne);
+
+   mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
+   {
+      constexpr int MAX_D1D = DofQuadLimits::MAX_D1D;
+      constexpr int MAX_Q1D = DofQuadLimits::HCURL_MAX_Q1D;
+
+      real_t mass[MAX_Q1D][MAX_Q1D][MAX_Q1D];
+      for (int qz = 0; qz < q1d; ++qz)
+      {
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            for (int qx = 0; qx < q1d; ++qx) { mass[qz][qy][qx] = 0.0; }
+         }
+      }
+
+      for (int dz = 0; dz < d1d_test; ++dz)
+      {
+         real_t mass_xy[MAX_Q1D][MAX_Q1D];
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            for (int qx = 0; qx < q1d; ++qx) { mass_xy[qy][qx] = 0.0; }
+         }
+
+         for (int dy = 0; dy < d1d_test; ++dy)
+         {
+            real_t mass_x[MAX_Q1D];
+            for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] = 0.0; }
+            for (int dx = 0; dx < d1d_test; ++dx)
+            {
+               const real_t t = X(dx, dy, dz, e);
+               for (int qx = 0; qx < q1d; ++qx) { mass_x[qx] += t * B(qx, dx); }
+            }
+            for (int qy = 0; qy < q1d; ++qy)
+            {
+               const real_t wy = B(qy, dy);
+               for (int qx = 0; qx < q1d; ++qx) { mass_xy[qy][qx] += mass_x[qx] * wy; }
+            }
+         }
+
+         for (int qz = 0; qz < q1d; ++qz)
+         {
+            const real_t wz = B(qz, dz);
+            for (int qy = 0; qy < q1d; ++qy)
+            {
+               for (int qx = 0; qx < q1d; ++qx) { mass[qz][qy][qx] += mass_xy[qy][qx] * wz; }
+            }
+         }
+      }
+
+      int osc = 0;
+      for (int qz = 0; qz < q1d; ++qz)
+      {
+         real_t mass_xy[MAX_D1D][MAX_D1D];
+         for (int dy = 0; dy < d1d; ++dy)
+         {
+            for (int dx = 0; dx < d1d - 1; ++dx) { mass_xy[dy][dx] = 0.0; }
+         }
+
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            real_t mass_x[MAX_D1D];
+            for (int dx = 0; dx < d1d - 1; ++dx) { mass_x[dx] = 0.0; }
+            for (int qx = 0; qx < q1d; ++qx)
+            {
+               const real_t s = O(0, qx, qy, qz, e) * mass[qz][qy][qx];
+               for (int dx = 0; dx < d1d - 1; ++dx) { mass_x[dx] += s * Bo(qx, dx); }
+            }
+            for (int dy = 0; dy < d1d; ++dy)
+            {
+               const real_t wy = Bc(qy, dy);
+               for (int dx = 0; dx < d1d - 1; ++dx) { mass_xy[dy][dx] += mass_x[dx] * wy; }
+            }
+         }
+
+         for (int dz = 0; dz < d1d; ++dz)
+         {
+            const real_t wz = Bc(qz, dz);
+            for (int dy = 0; dy < d1d; ++dy)
+            {
+               for (int dx = 0; dx < d1d - 1; ++dx)
+               {
+                  Y(dx + ((dy + (dz * d1d)) * (d1d - 1)) + osc, e) += mass_xy[dy][dx] * wz;
+               }
+            }
+         }
+      }
+
+      osc += (d1d - 1) * d1d * d1d;
+      for (int qz = 0; qz < q1d; ++qz)
+      {
+         real_t mass_xy[MAX_D1D][MAX_D1D];
+         for (int dy = 0; dy < d1d - 1; ++dy)
+         {
+            for (int dx = 0; dx < d1d; ++dx) { mass_xy[dy][dx] = 0.0; }
+         }
+
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            real_t mass_x[MAX_D1D];
+            for (int dx = 0; dx < d1d; ++dx) { mass_x[dx] = 0.0; }
+            for (int qx = 0; qx < q1d; ++qx)
+            {
+               const real_t s = O(1, qx, qy, qz, e) * mass[qz][qy][qx];
+               for (int dx = 0; dx < d1d; ++dx) { mass_x[dx] += s * Bc(qx, dx); }
+            }
+            for (int dy = 0; dy < d1d - 1; ++dy)
+            {
+               const real_t wy = Bo(qy, dy);
+               for (int dx = 0; dx < d1d; ++dx) { mass_xy[dy][dx] += mass_x[dx] * wy; }
+            }
+         }
+
+         for (int dz = 0; dz < d1d; ++dz)
+         {
+            const real_t wz = Bc(qz, dz);
+            for (int dy = 0; dy < d1d - 1; ++dy)
+            {
+               for (int dx = 0; dx < d1d; ++dx)
+               {
+                  Y(dx + ((dy + (dz * (d1d - 1))) * d1d) + osc, e) += mass_xy[dy][dx] * wz;
+               }
+            }
+         }
+      }
+
+      osc += (d1d - 1) * d1d * d1d;
+      for (int qz = 0; qz < q1d; ++qz)
+      {
+         real_t mass_xy[MAX_D1D][MAX_D1D];
+         for (int dy = 0; dy < d1d; ++dy)
+         {
+            for (int dx = 0; dx < d1d; ++dx) { mass_xy[dy][dx] = 0.0; }
+         }
+
+         for (int qy = 0; qy < q1d; ++qy)
+         {
+            real_t mass_x[MAX_D1D];
+            for (int dx = 0; dx < d1d; ++dx) { mass_x[dx] = 0.0; }
+            for (int qx = 0; qx < q1d; ++qx)
+            {
+               const real_t s = O(2, qx, qy, qz, e) * mass[qz][qy][qx];
+               for (int dx = 0; dx < d1d; ++dx) { mass_x[dx] += s * Bc(qx, dx); }
+            }
+            for (int dy = 0; dy < d1d; ++dy)
+            {
+               const real_t wy = Bc(qy, dy);
+               for (int dx = 0; dx < d1d; ++dx) { mass_xy[dy][dx] += mass_x[dx] * wy; }
+            }
+         }
+
+         for (int dz = 0; dz < d1d - 1; ++dz)
+         {
+            const real_t wz = Bo(qz, dz);
+            for (int dy = 0; dy < d1d; ++dy)
+            {
+               for (int dx = 0; dx < d1d; ++dx)
+               {
+                  Y(dx + ((dy + (dz * d1d)) * d1d) + osc, e) += mass_xy[dy][dx] * wz;
+               }
+            }
+         }
+      }
+   });
+}
+
+} // namespace
+
+void MixedDotProductIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
+                                           const FiniteElementSpace &test_fes)
+{
+   Mesh *mesh = trial_fes.GetMesh();
+   const FiniteElement *trial_fel = trial_fes.GetTypicalFE();
+   const FiniteElement *test_fel = test_fes.GetTypicalFE();
+
+   const VectorTensorFiniteElement *trial_el =
+      dynamic_cast<const VectorTensorFiniteElement*>(trial_fel);
+   MFEM_VERIFY(trial_el != NULL, "Only VectorTensorFiniteElement is supported!");
+
+   const TensorBasisElement *test_tensor_el =
+      dynamic_cast<const TensorBasisElement*>(test_fel);
+   MFEM_VERIFY(test_tensor_el != NULL,
+               "Only tensor-product scalar test elements are supported!");
+
+   MFEM_VERIFY(trial_el->GetDerivType() == mfem::FiniteElement::CURL,
+               "Only H(curl) trial spaces are supported!");
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == nullptr)
+   {
+      const int order = trial_fel->GetOrder() + test_fel->GetOrder()
+                        + mesh->GetTypicalElementTransformation()->OrderW();
+      ir = &IntRules.Get(trial_fel->GetGeomType(), order);
+   }
+
+   dim = mesh->Dimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Unsupported dimension!");
+   MFEM_VERIFY(trial_el->GetDim() == dim && test_fel->GetDim() == dim,
+               "Trial/test dimension mismatch.");
+
+   ne = trial_fes.GetNE();
+   MFEM_VERIFY(ne == test_fes.GetNE(),
+               "Different meshes for test and trial spaces");
+
+   geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
+   mapsC = &trial_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
+   mapsO = &trial_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
+   mapsTest = &test_fel->GetDofToQuad(*ir, DofToQuad::TENSOR);
+
+   dofs1D = mapsC->ndof;
+   dofs1Dtest = mapsTest->ndof;
+   quad1D = mapsC->nqpt;
+   test_map_integral = (test_fel->GetMapType() == FiniteElement::INTEGRAL);
+
+   MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
+   MFEM_VERIFY(quad1D == mapsTest->nqpt, "Trial/test quadrature mismatch");
+   MFEM_VERIFY(dofs1D <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D, "");
+   MFEM_VERIFY(dofs1Dtest <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(quad1D <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D, "");
+
+   const int nq = ir->GetNPoints();
+   if (dim == 2) { MFEM_VERIFY(nq == quad1D * quad1D, ""); }
+   else { MFEM_VERIFY(nq == quad1D * quad1D * quad1D, ""); }
+
+   QuadratureSpace qs(*mesh, *ir);
+   CoefficientVector coeff(*VQ, qs, CoefficientStorage::FULL);
+   MFEM_VERIFY(coeff.GetVDim() == dim, "Vector coefficient dimension mismatch.");
+
+   pa_data.SetSize(dim * nq * ne, Device::GetMemoryType());
+
+   if (dim == 2)
+   {
+      PAHcurlDotSetup2D(quad1D, ne, test_map_integral, ir->GetWeights(),
+                        geom->J, coeff, pa_data);
+   }
+   else
+   {
+      PAHcurlDotSetup3D(quad1D, ne, test_map_integral, ir->GetWeights(),
+                        geom->J, coeff, pa_data);
+   }
+}
+
+void MixedDotProductIntegrator::AddMultPA(const Vector &x, Vector &y) const
+{
+   if (dim == 2)
+   {
+      PAHcurlDotApply2D(dofs1D, dofs1Dtest, quad1D, ne,
+                        mapsO->B, mapsC->B, mapsTest->Bt, pa_data, x, y);
+   }
+   else if (dim == 3)
+   {
+      PAHcurlDotApply3D(dofs1D, dofs1Dtest, quad1D, ne,
+                        mapsO->B, mapsC->B, mapsTest->Bt, pa_data, x, y);
+   }
+   else
+   {
+      MFEM_ABORT("Unsupported dimension!");
+   }
+}
+
+void MixedDotProductIntegrator::AddMultTransposePA(const Vector &x,
+                                                   Vector &y) const
+{
+   if (dim == 2)
+   {
+      PAHcurlDotApplyTranspose2D(dofs1D, dofs1Dtest, quad1D, ne,
+                                 mapsO->B, mapsC->B, mapsTest->B,
+                                 pa_data, x, y);
+   }
+   else if (dim == 3)
+   {
+      PAHcurlDotApplyTranspose3D(dofs1D, dofs1Dtest, quad1D, ne,
+                                 mapsO->B, mapsC->B, mapsTest->B,
+                                 pa_data, x, y);
+   }
+   else
+   {
+      MFEM_ABORT("Unsupported dimension!");
+   }
+}
+
 void MixedScalarCurlIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
                                            const FiniteElementSpace &test_fes)
 {
@@ -320,6 +1093,7 @@ void MixedVectorWeakCurlIntegrator::AssemblePA(const FiniteElementSpace
    CoefficientVector coeff(qs, CoefficientStorage::FULL);
    if (Q) { coeff.Project(*Q); }
    else if (DQ) { coeff.Project(*DQ); }
+   else if (MQ) { MFEM_ABORT("Not implemented."); }
    else { coeff.SetConstant(1.0); }
 
    if (trialType == mfem::FiniteElement::CURL && dim == 3)

fem/integ/bilininteg_vectorfemass_kernels.hpp

@@ -0,0 +1,131 @@
+// 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_BILININTEG_VECTORFEMASS_KERNELS_HPP
+#define MFEM_BILININTEG_VECTORFEMASS_KERNELS_HPP
+
+#include "../../config/config.hpp"
+#include "../../general/array.hpp"
+#include "../../general/forall.hpp"
+#include "../../linalg/dtensor.hpp"
+#include "../../linalg/vector.hpp"
+#include "../bilininteg.hpp"
+
+#include "bilininteg_diffusion_kernels.hpp"
+#include "bilininteg_hcurl_kernels.hpp"
+#include "bilininteg_hdiv_kernels.hpp"
+#include "bilininteg_hcurlhdiv_kernels.hpp"
+
+namespace mfem
+{
+/// \cond DO_NOT_DOCUMENT
+namespace internal
+{
+namespace hcurlmass
+{
+constexpr int NBZ3D(int d1d, int q1d)
+{
+   if (d1d <= 1 || q1d <= 0)
+   {
+      return 1;
+   }
+   // assume q1d >= d1d
+   // z dimension is capped at 64 on nvidia and amd gpus
+   int tmp = std::min((128 + q1d * q1d * q1d - 1) / (q1d * q1d * q1d), 64);
+   int smem_req =
+      sizeof(mfem::real_t) *
+      (3 * ((d1d - 1) * d1d * d1d + 2 * q1d * q1d * q1d) * tmp +
+       q1d * (d1d - 1) + q1d * d1d);
+   // assume GPU has at least 48k shared memory
+   return std::max(std::min(tmp, (48 * 1024 + smem_req - 1) / smem_req), 1);
+}
+} // namespace hcurlmass
+} // namespace internal
+
+template <FiniteElement::DerivType TrialType, FiniteElement::DerivType TestType,
+          int DIM, int TrialD1D, int TestD1D, int Q1D>
+VectorFEMassIntegrator::ApplyKernelType
+VectorFEMassIntegrator::ApplyPAKernels::Kernel()
+{
+   constexpr bool trial_curl = (TrialType == mfem::FiniteElement::CURL);
+   constexpr bool trial_div = (TrialType == mfem::FiniteElement::DIV);
+   constexpr bool test_curl = (TestType == mfem::FiniteElement::CURL);
+   constexpr bool test_div = (TestType == mfem::FiniteElement::DIV);
+
+   if constexpr (DIM == 3)
+   {
+      if constexpr (trial_curl && test_curl)
+      {
+         if (Device::Allows(Backend::DEVICE_MASK))
+         {
+            // assume TrialD1D == TestD1D
+            return internal::SmemPAHcurlMassApply3D<
+                   TrialD1D, Q1D, internal::hcurlmass::NBZ3D(TrialD1D, Q1D)>;
+         }
+         else
+         {
+            return internal::PAHcurlMassApply3D;
+         }
+      }
+      else if constexpr (trial_div && test_div)
+      {
+         if (Device::Allows(Backend::DEVICE_MASK))
+         {
+            // assumes TrialD1D == TestD1D
+            return internal::SmemPAHdivMassApply3D<TrialD1D, Q1D>;
+         }
+         else
+         {
+            return internal::PAHdivMassApply3D;
+         }
+      }
+      else if constexpr (trial_curl && test_div)
+      {
+         return internal::PAHdivHcurlMassApply3D;
+      }
+      else if constexpr (trial_div && test_curl)
+      {
+         return internal::PAHcurlHdivMassApply3D;
+      }
+   }
+   else if constexpr (DIM == 2) // 2D
+   {
+      if constexpr (trial_curl && test_curl)
+      {
+         return internal::PAHcurlMassApply2D;
+      }
+      else if constexpr (trial_div && test_div)
+      {
+         if (Device::Allows(Backend::DEVICE_MASK))
+         {
+            // assumes TrialD1D == TestD1D
+            return internal::SmemPAHdivMassApply2D<TrialD1D, Q1D>;
+         }
+         else
+         {
+            return internal::PAHdivMassApply2D;
+         }
+      }
+      else if constexpr (trial_curl && test_div)
+      {
+         return internal::PAHdivHcurlMassApply2D;
+      }
+      else if constexpr (trial_div && test_curl)
+      {
+         return internal::PAHcurlHdivMassApply2D;
+      }
+   }
+   MFEM_ABORT("Unknown kernel.");
+}
+/// \endcond DO_NOT_DOCUMENT
+}
+
+#endif

fem/integ/bilininteg_vectorfemass_pa.cpp

@@ -12,13 +12,115 @@
 #include "../bilininteg.hpp"
 #include "../gridfunc.hpp"
 #include "../qfunction.hpp"
-#include "bilininteg_diffusion_kernels.hpp"
-#include "bilininteg_hcurl_kernels.hpp"
-#include "bilininteg_hdiv_kernels.hpp"
-#include "bilininteg_hcurlhdiv_kernels.hpp"
+#include "bilininteg_vectorfemass_kernels.hpp"
 
 namespace mfem
 {
+/// \cond DO_NOT_DOCUMENT
+VectorFEMassIntegrator::ApplyKernelType
+VectorFEMassIntegrator::ApplyPAKernels::Fallback(
+   FiniteElement::DerivType TrialType, FiniteElement::DerivType TestType,
+   int dim, int, int, int)
+{
+   const bool trial_curl = (TrialType == mfem::FiniteElement::CURL);
+   const bool trial_div = (TrialType == mfem::FiniteElement::DIV);
+   const bool test_curl = (TestType == mfem::FiniteElement::CURL);
+   const bool test_div = (TestType == mfem::FiniteElement::DIV);
+
+   if (dim == 3)
+   {
+      if (trial_curl && test_curl)
+      {
+         return internal::PAHcurlMassApply3D;
+      }
+      else if (trial_div && test_div)
+      {
+         return internal::PAHdivMassApply3D;
+      }
+      else if (trial_curl && test_div)
+      {
+         return internal::PAHdivHcurlMassApply3D;
+      }
+      else if (trial_div && test_curl)
+      {
+         return internal::PAHcurlHdivMassApply3D;
+      }
+   }
+   else if (dim == 2) // 2D
+   {
+      if (trial_curl && test_curl)
+      {
+         return internal::PAHcurlMassApply2D;
+      }
+      else if (trial_div && test_div)
+      {
+         return internal::PAHdivMassApply2D;
+      }
+      else if (trial_curl && test_div)
+      {
+         return internal::PAHdivHcurlMassApply2D;
+      }
+      else if (trial_div && test_curl)
+      {
+         return internal::PAHcurlHdivMassApply2D;
+      }
+   }
+   MFEM_ABORT("Unknown kernel.");
+}
+/// \endcond DO_NOT_DOCUMENT
+
+VectorFEMassIntegrator::Kernels::Kernels()
+{
+   // h(curl), h(curl)
+   // P = Q (3D)
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 2, 2, 2>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 3, 3, 3>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 4, 4, 4>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 5, 5, 5>();
+
+   // P = Q + 1 (3D)
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 2, 2, 3>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 3, 3, 4>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 4, 4, 5>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 5, 5, 6>();
+   // h(div), h(div)
+   // P = Q (2D)
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 2, 2, 2, 2>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 2, 3, 3, 3>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 2, 4, 4, 4>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 2, 5, 5, 5>();
+
+   // P = Q + 1 (3D)
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 3, 2, 2, 3>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 3, 3, 3, 4>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 3, 4, 4, 5>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 3, 5, 5, 6>();
+}
+
+void VectorFEMassIntegrator::Init(Coefficient *q, DiagonalMatrixCoefficient *dq,
+                                  MatrixCoefficient *mq)
+{
+   static Kernels kernels{};
+   Q = q;
+   DQ = dq;
+   MQ = mq;
+}
 
 void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
@@ -67,8 +169,8 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 
    MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
 
-   trial_fetype = trial_el->GetDerivType();
-   test_fetype = test_el->GetDerivType();
+   trial_fetype = static_cast<FiniteElement::DerivType>(trial_el->GetDerivType());
+   test_fetype = static_cast<FiniteElement::DerivType>(test_el->GetDerivType());
 
    const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
    const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
@@ -215,225 +317,36 @@ void VectorFEMassIntegrator::AssembleDiagonalPA(Vector& diag)
 
 void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
-   const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
-   const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
-   const bool test_curl = (test_fetype == mfem::FiniteElement::CURL);
-   const bool test_div = (test_fetype == mfem::FiniteElement::DIV);
-
-   if (dim == 3)
-   {
-      if (trial_curl && test_curl)
-      {
-         if (Device::Allows(Backend::DEVICE_MASK))
-         {
-            const int ID = (dofs1D << 4) | quad1D;
-            switch (ID)
-            {
-               case 0x23:
-                  return internal::SmemPAHcurlMassApply3D<2,3>(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-               case 0x34:
-                  return internal::SmemPAHcurlMassApply3D<3,4>(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-               case 0x45:
-                  return internal::SmemPAHcurlMassApply3D<4,5>(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-               case 0x56:
-                  return internal::SmemPAHcurlMassApply3D<5,6>(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-               default:
-                  return internal::SmemPAHcurlMassApply3D(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-            }
-         }
-         else
-         {
-            internal::PAHcurlMassApply3D(dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
-                                         mapsO->Bt, mapsC->Bt, pa_data, x, y);
-         }
-      }
-      else if (trial_div && test_div)
-      {
-         internal::PAHdivMassApply(3, dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
-                                   mapsO->Bt, mapsC->Bt, pa_data, x, y);
-      }
-      else if (trial_curl && test_div)
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
-                                          true, false, mapsO->B, mapsC->B, mapsOtest->Bt,
-                                          mapsCtest->Bt, pa_data, x, y);
-      }
-      else if (trial_div && test_curl)
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
-                                          false, false, mapsO->B, mapsC->B, mapsOtest->Bt,
-                                          mapsCtest->Bt, pa_data, x, y);
-      }
-      else
-      {
-         MFEM_ABORT("Unknown kernel.");
-      }
-   }
-   else // 2D
-   {
-      if (trial_curl && test_curl)
-      {
-         internal::PAHcurlMassApply2D(dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
-                                      mapsO->Bt, mapsC->Bt, pa_data, x, y);
-      }
-      else if (trial_div && test_div)
-      {
-         internal::PAHdivMassApply(2, dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
-                                   mapsO->Bt,
-                                   mapsC->Bt, pa_data, x, y);
-      }
-      else if ((trial_curl && test_div) || (trial_div && test_curl))
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
-                                          trial_curl, false, mapsO->B, mapsC->B,
-                                          mapsOtest->Bt, mapsCtest->Bt, pa_data, x, y);
-      }
-      else
-      {
-         MFEM_ABORT("Unknown kernel.");
-      }
-   }
+   const bool scalar_coeff = !(DQ || MQ);
+   ApplyPAKernels::Run(trial_fetype, test_fetype, dim, dofs1D, dofs1Dtest,
+                       quad1D, ne, symmetric, scalar_coeff, mapsO->B, mapsC->B,
+                       mapsOtest->Bt, mapsCtest->Bt, pa_data, x, y, dofs1D,
+                       dofs1Dtest, quad1D);
 }
 
 void VectorFEMassIntegrator::AddAbsMultPA(const Vector &x, Vector &y) const
 {
-   const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
-   const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
-   const bool test_curl = (test_fetype == mfem::FiniteElement::CURL);
-   const bool test_div = (test_fetype == mfem::FiniteElement::DIV);
+   const bool scalar_coeff = !(DQ || MQ);
 
    Vector abs_pa_data(pa_data);
    abs_pa_data.Abs();
 
    Array<real_t> absBo(mapsO->B);
    Array<real_t> absBc(mapsC->B);
-   Array<real_t> absBto(mapsO->Bt);
-   Array<real_t> absBtc(mapsC->Bt);
+   // Array<real_t> absBto(mapsO->Bt);
+   // Array<real_t> absBtc(mapsC->Bt);
    Array<real_t> absBto_t(mapsOtest->Bt);
    Array<real_t> absBtc_t(mapsCtest->Bt);
 
    absBo.Abs();
    absBc.Abs();
-   absBto.Abs();
-   absBtc.Abs();
    absBto_t.Abs();
    absBtc_t.Abs();
 
-   if (dim == 3)
-   {
-      if (trial_curl && test_curl)
-      {
-         if (Device::Allows(Backend::DEVICE_MASK))
-         {
-            const int ID = (dofs1D << 4) | quad1D;
-            switch (ID)
-            {
-               case 0x23:
-                  return internal::SmemPAHcurlMassApply3D<2,3>(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-               case 0x34:
-                  return internal::SmemPAHcurlMassApply3D<3,4>(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-               case 0x45:
-                  return internal::SmemPAHcurlMassApply3D<4,5>(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-               case 0x56:
-                  return internal::SmemPAHcurlMassApply3D<5,6>(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-               default:
-                  return internal::SmemPAHcurlMassApply3D(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-            }
-         }
-         else
-         {
-            internal::PAHcurlMassApply3D(dofs1D, quad1D, ne, symmetric,
-                                         absBo, absBc, absBto, absBtc,
-                                         abs_pa_data, x, y);
-         }
-      }
-      else if (trial_div && test_div)
-      {
-         internal::PAHdivMassApply(3, dofs1D, quad1D, ne, symmetric,
-                                   absBo, absBc, absBto, absBtc,
-                                   abs_pa_data, x, y);
-      }
-      else if (trial_curl && test_div)
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne,
-                                          scalarCoeff, true, false,
-                                          absBo, absBc, absBto_t, absBtc_t,
-                                          abs_pa_data, x, y);
-      }
-      else if (trial_div && test_curl)
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne,
-                                          scalarCoeff, false, false,
-                                          absBo, absBc, absBto_t, absBtc_t,
-                                          abs_pa_data, x, y);
-      }
-      else
-      {
-         MFEM_ABORT("Unknown kernel.");
-      }
-   }
-   else // 2D
-   {
-      if (trial_curl && test_curl)
-      {
-         internal::PAHcurlMassApply2D(dofs1D, quad1D, ne, symmetric,
-                                      absBo, absBc, absBto, absBtc,
-                                      abs_pa_data, x, y);
-      }
-      else if (trial_div && test_div)
-      {
-         internal::PAHdivMassApply(2, dofs1D, quad1D, ne, symmetric,
-                                   absBo, absBc, absBto, absBtc,
-                                   abs_pa_data, x, y);
-      }
-      else if ((trial_curl && test_div) || (trial_div && test_curl))
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne,
-                                          scalarCoeff, trial_curl, false,
-                                          absBo, absBc, absBto_t, absBtc_t,
-                                          abs_pa_data, x, y);
-      }
-      else
-      {
-         MFEM_ABORT("Unknown kernel.");
-      }
-   }
+   ApplyPAKernels::Run(trial_fetype, test_fetype, dim, dofs1D, dofs1Dtest,
+                       quad1D, ne, symmetric, scalar_coeff, absBo, absBc,
+                       absBto_t, absBtc_t, abs_pa_data, x, y, dofs1D,
+                       dofs1Dtest, quad1D);
 }
 
 void VectorFEMassIntegrator::AddMultTransposePA(const Vector &x,

fem/integ/lininteg_domain_kernels.hpp

@@ -307,6 +307,506 @@ DomainLFIntegrator::AssembleKernels::Kernel()
    MFEM_ABORT("");
 }
 
+template <int T_D1D = 0, int T_Q1D = 0>
+static void HdivDLFAssemble2D(const int ne, const Array<int> &markers,
+                              const Vector &jac, const Array<real_t> &weights,
+                              const Array<real_t> &testBO,
+                              const Array<real_t> &testBC, const Vector &coeff,
+                              Vector &y, const int d, const int q)
+{
+   MFEM_VERIFY(T_D1D || d <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
+               "Problem size too large.");
+   MFEM_VERIFY(T_Q1D || q <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
+               "Problem size too large.");
+   MFEM_VERIFY(y.Size() == 2 * (d - 1) * d * ne, "");
+
+   constexpr int vdim = 2;
+   const auto F = coeff.Read();
+   const auto M = markers.Read();
+   const auto BO = Reshape(testBO.Read(), q, d-1);
+   const auto BC = Reshape(testBC.Read(), q, d);
+   const auto J = Reshape(jac.Read(), q, q, vdim, vdim, ne);
+   const auto W = Reshape(weights.Read(), q, q);
+   const bool cst = coeff.Size() == vdim;
+   const auto C = cst ? Reshape(F,vdim,1,1,1) : Reshape(F,vdim,q,q,ne);
+   auto Y = y.ReadWrite();
+
+   mfem::forall_3D(ne, q, q, vdim, [=] MFEM_HOST_DEVICE (int e)
+   {
+      constexpr int vdim = 2;
+      if (M[e] == 0) { return; } // ignore
+
+      constexpr int Q = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
+      constexpr int D = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
+
+      MFEM_SHARED real_t sBot[Q*D];
+      MFEM_SHARED real_t sBct[Q*D];
+      MFEM_SHARED real_t sQQ[vdim*Q*Q];
+      MFEM_SHARED real_t sQD[vdim*Q*D];
+
+      // Bo and Bc into shared memory
+      const DeviceMatrix Bot(sBot, d-1, q);
+      kernels::internal::LoadB<D,Q>(d-1, q, BO, sBot);
+      const DeviceMatrix Bct(sBct, d, q);
+      kernels::internal::LoadB<D,Q>(d, q, BC, sBct);
+
+      const DeviceCube QQ(sQQ, q, q, vdim);
+      const DeviceCube QD(sQD, q, d, vdim);
+
+      MFEM_FOREACH_THREAD(vd,z,vdim)
+      {
+         const real_t cst_val_0 = C(0,0,0,0);
+         const real_t cst_val_1 = C(1,0,0,0);
+         MFEM_FOREACH_THREAD(y,y,q)
+         {
+            MFEM_FOREACH_THREAD(x,x,q)
+            {
+               const real_t J0 = J(x,y,0,vd,e);
+               const real_t J1 = J(x,y,1,vd,e);
+               const real_t C0 = cst ? cst_val_0 : C(0,x,y,e);
+               const real_t C1 = cst ? cst_val_1 : C(1,x,y,e);
+               QQ(x,y,vd) = W(x,y)*(J0*C0 + J1*C1);
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(vd,z,vdim)
+      {
+         const int nx = (vd == 0) ? d : d-1;
+         DeviceMatrix Btx = (vd == 0) ? Bct : Bot;
+         MFEM_FOREACH_THREAD(qy,y,q)
+         {
+            MFEM_FOREACH_THREAD(dx,x,nx)
+            {
+               real_t qd = 0.0;
+               for (int qx = 0; qx < q; ++qx)
+               {
+                  qd += QQ(qx,qy,vd) * Btx(dx,qx);
+               }
+               QD(dx,qy,vd) = qd;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(vd,z,vdim)
+      {
+         const int nx = (vd == 0) ? d : d-1;
+         const int ny = (vd == 1) ? d : d-1;
+         DeviceMatrix Bty = (vd == 1) ? Bct : Bot;
+         DeviceTensor<4> Yxy(Y, nx, ny, vdim, ne);
+         MFEM_FOREACH_THREAD(dy,y,ny)
+         {
+            MFEM_FOREACH_THREAD(dx,x,nx)
+            {
+               real_t dd = 0.0;
+               for (int qy = 0; qy < q; ++qy)
+               {
+                  dd += QD(dx,qy,vd) * Bty(dy,qy);
+               }
+               Yxy(dx,dy,vd,e) += dd;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+   });
+}
+
+template <int T_D1D = 0, int T_Q1D = 0>
+static void HdivDLFAssemble3D(const int ne, const Array<int> &markers,
+                              const Vector &jac, const Array<real_t> &weights,
+                              const Array<real_t> &testBO,
+                              const Array<real_t> &testBC, const Vector &coeff,
+                              Vector &y, const int d, const int q)
+{
+   MFEM_VERIFY(T_D1D || d <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
+               "Problem size too large.");
+   MFEM_VERIFY(T_Q1D || q <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
+               "Problem size too large.");
+   MFEM_VERIFY(y.Size() == 3 * (d - 1) * (d - 1) * d * ne, "y wrong length");
+
+   constexpr int vdim = 3;
+   const auto F = coeff.Read();
+   const auto M = markers.Read();
+   const auto BO = Reshape(testBO.Read(), q, d-1);
+   const auto BC = Reshape(testBC.Read(), q, d);
+   const auto J = Reshape(jac.Read(), q, q, q, vdim, vdim, ne);
+   const auto W = Reshape(weights.Read(), q, q, q);
+   const bool cst = coeff.Size() == vdim;
+   const auto C = cst ? Reshape(F,vdim,1,1,1,1) : Reshape(F,vdim,q,q,q,ne);
+   auto Y = y.ReadWrite();
+
+   mfem::forall_3D(ne, q, q, vdim, [=] MFEM_HOST_DEVICE (int e)
+   {
+      constexpr int vdim = 3;
+      if (M[e] == 0) { return; } // ignore
+
+      constexpr int Q = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
+      constexpr int D = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
+
+      MFEM_SHARED real_t sBot[Q*D];
+      MFEM_SHARED real_t sBct[Q*D];
+
+      // Bo and Bc into shared memory
+      const DeviceMatrix Bot(sBot, d-1, q);
+      kernels::internal::LoadB<D,Q>(d-1, q, BO, sBot);
+      const DeviceMatrix Bct(sBct, d, q);
+      kernels::internal::LoadB<D,Q>(d, q, BC, sBct);
+
+      MFEM_SHARED real_t sm0[vdim*Q*Q*Q];
+      MFEM_SHARED real_t sm1[vdim*Q*Q*Q];
+      DeviceTensor<4> QQQ(sm1, q, q, q, vdim);
+      DeviceTensor<4> DQQ(sm0, d, q, q, vdim);
+      DeviceTensor<4> DDQ(sm1, d, d, q, vdim);
+
+      MFEM_FOREACH_THREAD(vd,z,vdim)
+      {
+         const real_t cst_val_0 = C(0,0,0,0,0);
+         const real_t cst_val_1 = C(1,0,0,0,0);
+         const real_t cst_val_2 = C(2,0,0,0,0);
+         MFEM_FOREACH_THREAD(y,y,q)
+         {
+            MFEM_FOREACH_THREAD(x,x,q)
+            {
+               for (int z = 0; z < q; ++z)
+               {
+                  const real_t J0 = J(x,y,z,0,vd,e);
+                  const real_t J1 = J(x,y,z,1,vd,e);
+                  const real_t J2 = J(x,y,z,2,vd,e);
+                  const real_t C0 = cst ? cst_val_0 : C(0,x,y,z,e);
+                  const real_t C1 = cst ? cst_val_1 : C(1,x,y,z,e);
+                  const real_t C2 = cst ? cst_val_2 : C(2,x,y,z,e);
+                  QQQ(x,y,z,vd) = W(x,y,z)*(J0*C0 + J1*C1 + J2*C2);
+               }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      // Apply Bt operator
+      MFEM_FOREACH_THREAD(vd,z,vdim)
+      {
+         const int nx = (vd == 0) ? d : d-1;
+         DeviceMatrix Btx = (vd == 0) ? Bct : Bot;
+         MFEM_FOREACH_THREAD(qy,y,q)
+         {
+            MFEM_FOREACH_THREAD(dx,x,nx)
+            {
+               real_t u[Q];
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz) { u[qz] = 0.0; }
+               MFEM_UNROLL(Q)
+               for (int qx = 0; qx < q; ++qx)
+               {
+                  MFEM_UNROLL(Q)
+                  for (int qz = 0; qz < q; ++qz)
+                  {
+                     u[qz] += QQQ(qx,qy,qz,vd) * Btx(dx,qx);
+                  }
+               }
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz) { DQQ(dx,qy,qz,vd) = u[qz]; }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(vd,z,vdim)
+      {
+         const int nx = (vd == 0) ? d : d-1;
+         const int ny = (vd == 1) ? d : d-1;
+         DeviceMatrix Bty = (vd == 1) ? Bct : Bot;
+         MFEM_FOREACH_THREAD(dy,y,ny)
+         {
+            MFEM_FOREACH_THREAD(dx,x,nx)
+            {
+               real_t u[Q];
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz) { u[qz] = 0.0; }
+               MFEM_UNROLL(Q)
+               for (int qy = 0; qy < q; ++qy)
+               {
+                  MFEM_UNROLL(Q)
+                  for (int qz = 0; qz < q; ++qz)
+                  {
+                     u[qz] += DQQ(dx,qy,qz,vd) * Bty(dy,qy);
+                  }
+               }
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz) { DDQ(dx,dy,qz,vd) = u[qz]; }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(vd,z,vdim)
+      {
+         const int nx = (vd == 0) ? d : d-1;
+         const int ny = (vd == 1) ? d : d-1;
+         const int nz = (vd == 2) ? d : d-1;
+         DeviceTensor<5> Yxyz(Y, nx, ny, nz, vdim, ne);
+         DeviceMatrix Btz = (vd == 2) ? Bct : Bot;
+         MFEM_FOREACH_THREAD(dy,y,ny)
+         {
+            MFEM_FOREACH_THREAD(dx,x,nx)
+            {
+               real_t u[D];
+               MFEM_UNROLL(D)
+               for (int dz = 0; dz < nz; ++dz) { u[dz] = 0.0; }
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz)
+               {
+                  MFEM_UNROLL(D)
+                  for (int dz = 0; dz < nz; ++dz)
+                  {
+                     u[dz] += DDQ(dx,dy,qz,vd) * Btz(dz,qz);
+                  }
+               }
+               MFEM_UNROLL(D)
+               for (int dz = 0; dz < nz; ++dz) { Yxyz(dx,dy,dz,vd,e) += u[dz]; }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+   });
+}
+
+/// @param ne number of elements
+/// @param markers array where entry markers[e] == 0 to skip assembly over
+/// element e element
+/// @param jac Spatial Jacobians evaluated at all quadrature points
+/// @param weights 1D quadrature weights
+/// @param testBO 1D open basis test functions
+/// @param testBC 1D closed basis test functions
+/// @param coeff coefficient values evaluated at quadrature points, possibly
+/// compressed.
+/// @param d number of 1D closed dofs
+/// @param q number of 1D quadrature points
+/// @tparam T_D1D maximum number of dofs along any direction, or 0
+/// @tparam T_Q1D maximum number of quadrature points along any direction, or 0
+template <int T_D1D = 0, int T_Q1D = 0>
+static void HcurlDLFAssemble3D(const int ne, const Array<int> &markers,
+                               const Vector &jac, const Array<real_t> &weights,
+                               const Array<real_t> &testBO,
+                               const Array<real_t> &testBC, const Vector &coeff,
+                               Vector &y, const int d, const int q)
+{
+   MFEM_VERIFY(T_D1D || d <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
+               "Problem size too large.");
+   MFEM_VERIFY(T_Q1D || q <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,
+               "Problem size too large.");
+   MFEM_VERIFY(y.Size() == 3 * (d - 1) * d * d * ne, "y wrong length");
+
+   constexpr int vdim = 3;
+   const auto F = coeff.Read();
+   const auto M = markers.Read();
+   const auto BO = Reshape(testBO.Read(), q, d-1);
+   const auto BC = Reshape(testBC.Read(), q, d);
+   const auto J = Reshape(jac.Read(), q, q, q, vdim, vdim, ne);
+   const auto W = Reshape(weights.Read(), q, q, q);
+   const bool cst = coeff.Size() == vdim;
+   const auto C = cst ? Reshape(F,vdim,1,1,1,1) : Reshape(F,vdim,q,q,q,ne);
+   auto Y = y.ReadWrite();
+
+   mfem::forall_3D(ne, q, q, vdim, [=] MFEM_HOST_DEVICE(int e)
+   {
+      if (M[e] == 0)
+      {
+         // ignore
+         return;
+      }
+
+      constexpr int vdim = 3;
+      constexpr int Q = T_Q1D ? T_Q1D : DofQuadLimits::HCURL_MAX_Q1D;
+      constexpr int D = T_D1D ? T_D1D : DofQuadLimits::HCURL_MAX_D1D;
+
+      MFEM_SHARED real_t sBot[Q * D];
+      MFEM_SHARED real_t sBct[Q * D];
+
+      // Bo and Bc into shared memory
+      const DeviceMatrix Bot(sBot, d - 1, q);
+      kernels::internal::LoadB<D, Q>(d - 1, q, BO, sBot);
+      const DeviceMatrix Bct(sBct, d, q);
+      kernels::internal::LoadB<D, Q>(d, q, BC, sBct);
+
+      MFEM_SHARED real_t sm0[vdim * Q * Q * Q];
+      MFEM_SHARED real_t sm1[vdim * Q * Q * Q];
+      DeviceTensor<4> QQQ(sm1, q, q, q, vdim);
+      DeviceTensor<4> DQQ(sm0, d, q, q, vdim);
+      DeviceTensor<4> DDQ(sm1, d, d, q, vdim);
+
+      const real_t cst_val_0 = C(0, 0, 0, 0, 0);
+      const real_t cst_val_1 = C(1, 0, 0, 0, 0);
+      const real_t cst_val_2 = C(2, 0, 0, 0, 0);
+
+      MFEM_FOREACH_THREAD(vd, z, vdim)
+      {
+         MFEM_FOREACH_THREAD(y, y, q)
+         {
+            MFEM_FOREACH_THREAD(x, x, q)
+            {
+               for (int z = 0; z < q; ++z)
+               {
+                  real_t curr[3];
+                  curr[0] = cst ? cst_val_0 : C(0, x, y, z, e);
+                  curr[1] = cst ? cst_val_1 : C(1, x, y, z, e);
+                  curr[2] = cst ? cst_val_2 : C(2, x, y, z, e);
+
+                  const real_t J11 = J(x, y, z, 0, 0, e);
+                  const real_t J21 = J(x, y, z, 1, 0, e);
+                  const real_t J31 = J(x, y, z, 2, 0, e);
+                  const real_t J12 = J(x, y, z, 0, 1, e);
+                  const real_t J22 = J(x, y, z, 1, 1, e);
+                  const real_t J32 = J(x, y, z, 2, 1, e);
+                  const real_t J13 = J(x, y, z, 0, 2, e);
+                  const real_t J23 = J(x, y, z, 1, 2, e);
+                  const real_t J33 = J(x, y, z, 2, 2, e);
+                  // adj(J)
+                  const real_t A11 = (J22 * J33) - (J23 * J32);
+                  const real_t A12 = (J32 * J13) - (J12 * J33);
+                  const real_t A13 = (J12 * J23) - (J22 * J13);
+                  const real_t A21 = (J31 * J23) - (J21 * J33);
+                  const real_t A22 = (J11 * J33) - (J13 * J31);
+                  const real_t A23 = (J21 * J13) - (J11 * J23);
+                  const real_t A31 = (J21 * J32) - (J31 * J22);
+                  const real_t A32 = (J31 * J12) - (J11 * J32);
+                  const real_t A33 = (J11 * J22) - (J12 * J21);
+                  const real_t A[9] = {A11, A12, A13, A21, A22,
+                                       A23, A31, A32, A33
+                                      };
+                  QQQ(x, y, z, vd) = W(x, y, z) * (A[vd * vdim] * curr[0] +
+                                                   A[vd * vdim + 1] * curr[1] +
+                                                   A[vd * vdim + 2] * curr[2]);
+               }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      // Apply Bt operator
+      MFEM_FOREACH_THREAD(vd, z, vdim)
+      {
+         const int nx = (vd == 0) ? d - 1 : d;
+         DeviceMatrix Btx = (vd == 0) ? Bot : Bct;
+         MFEM_FOREACH_THREAD(qy, y, q)
+         {
+            MFEM_FOREACH_THREAD(dx, x, nx)
+            {
+               real_t u[Q];
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz)
+               {
+                  u[qz] = 0.0;
+               }
+               MFEM_UNROLL(Q)
+               for (int qx = 0; qx < q; ++qx)
+               {
+                  MFEM_UNROLL(Q)
+                  for (int qz = 0; qz < q; ++qz)
+                  {
+                     u[qz] += QQQ(qx, qy, qz, vd) * Btx(dx, qx);
+                  }
+               }
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz)
+               {
+                  DQQ(dx, qy, qz, vd) = u[qz];
+               }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(vd, z, vdim)
+      {
+         const int nx = (vd == 0) ? d - 1 : d;
+         const int ny = (vd == 1) ? d - 1 : d;
+         DeviceMatrix Bty = (vd == 1) ? Bot : Bct;
+         MFEM_FOREACH_THREAD(dy, y, ny)
+         {
+            MFEM_FOREACH_THREAD(dx, x, nx)
+            {
+               real_t u[Q];
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz)
+               {
+                  u[qz] = 0.0;
+               }
+               MFEM_UNROLL(Q)
+               for (int qy = 0; qy < q; ++qy)
+               {
+                  MFEM_UNROLL(Q)
+                  for (int qz = 0; qz < q; ++qz)
+                  {
+                     u[qz] += DQQ(dx, qy, qz, vd) * Bty(dy, qy);
+                  }
+               }
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz)
+               {
+                  DDQ(dx, dy, qz, vd) = u[qz];
+               }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(vd, z, vdim)
+      {
+         const int nx = (vd == 0) ? d - 1 : d;
+         const int ny = (vd == 1) ? d - 1 : d;
+         const int nz = (vd == 2) ? d - 1 : d;
+         DeviceTensor<5> Yxyz(Y, nx, ny, nz, vdim, ne);
+         DeviceMatrix Btz = (vd == 2) ? Bot : Bct;
+         MFEM_FOREACH_THREAD(dy, y, ny)
+         {
+            MFEM_FOREACH_THREAD(dx, x, nx)
+            {
+               real_t u[D];
+               MFEM_UNROLL(D)
+               for (int dz = 0; dz < nz; ++dz)
+               {
+                  u[dz] = 0.0;
+               }
+               MFEM_UNROLL(Q)
+               for (int qz = 0; qz < q; ++qz)
+               {
+                  MFEM_UNROLL(D)
+                  for (int dz = 0; dz < nz; ++dz)
+                  {
+                     u[dz] += DDQ(dx, dy, qz, vd) * Btz(dz, qz);
+                  }
+               }
+               MFEM_UNROLL(D)
+               for (int dz = 0; dz < nz; ++dz)
+               {
+                  Yxyz(dx, dy, dz, vd, e) += u[dz];
+               }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+   });
+}
+
+template <FiniteElement::DerivType TestType, int DIM, int TEST_D1D, int Q1D>
+VectorFEDomainLFIntegrator::AssembleKernelType
+VectorFEDomainLFIntegrator::AssembleKernels::Kernel()
+{
+   if constexpr (TestType == FiniteElement::DIV)
+   {
+      if constexpr (DIM == 2)
+      {
+         return HdivDLFAssemble2D<TEST_D1D, Q1D>;
+      }
+      if constexpr (DIM == 3)
+      {
+         return HdivDLFAssemble3D<TEST_D1D, Q1D>;
+      }
+   }
+   if constexpr (TestType == FiniteElement::CURL)
+   {
+      if constexpr (DIM == 3)
+      {
+         return HcurlDLFAssemble3D<TEST_D1D, Q1D>;
+      }
+   }
+   MFEM_ABORT("");
+}
+
 /// \endcond DO_NOT_DOCUMENT
 
 } // namespace mfem

fem/integ/lininteg_domain_vectorfe.cpp

@@ -13,317 +13,76 @@
 #include "../../fem/kernels.hpp"
 #include "../fem.hpp"
 
+#include "lininteg_domain_kernels.hpp"
+
 namespace mfem
 {
 
-template<int T_D1D = 0, int T_Q1D = 0>
-static void HdivDLFAssemble2D(
-   const int ne, const int d, const int q, const int *markers, const real_t *bo,
-   const real_t *bc, const real_t *j, const real_t *weights,
-   const Vector &coeff, real_t *y)
+VectorFEDomainLFIntegrator::Kernels::Kernels()
 {
-   MFEM_VERIFY(T_D1D || d <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
-               "Problem size too large.");
-   MFEM_VERIFY(T_Q1D || q <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
-               "Problem size too large.");
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 2, 1, 1>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 2, 2, 2>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 2, 3, 3>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 2, 4, 4>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 2, 5, 5>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 2, 6, 6>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 2, 7, 7>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 2, 8, 8>();
 
-   static constexpr int vdim = 2;
-   const auto F = coeff.Read();
-   const auto M = Reshape(markers, ne);
-   const auto BO = Reshape(bo, q, d-1);
-   const auto BC = Reshape(bc, q, d);
-   const auto J = Reshape(j, q, q, vdim, vdim, ne);
-   const auto W = Reshape(weights, q, q);
-   const bool cst = coeff.Size() == vdim;
-   const auto C = cst ? Reshape(F,vdim,1,1,1) : Reshape(F,vdim,q,q,ne);
-   auto Y = Reshape(y, 2*(d-1)*d, ne);
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 3, 1, 1>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 3, 2, 2>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 3, 3, 3>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 3, 4, 4>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 3, 5, 5>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 3, 6, 6>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 3, 7, 7>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::DIV, 3, 8, 8>();
 
-   mfem::forall_3D(ne, q, q, vdim, [=] MFEM_HOST_DEVICE (int e)
-   {
-      if (M(e) == 0) { return; } // ignore
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 1, 1>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 2, 2>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 3, 3>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 4, 4>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 5, 5>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 6, 6>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 7, 7>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 8, 8>();
 
-      constexpr int Q = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
-      constexpr int D = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
-
-      MFEM_SHARED real_t sBot[Q*D];
-      MFEM_SHARED real_t sBct[Q*D];
-      MFEM_SHARED real_t sQQ[vdim*Q*Q];
-      MFEM_SHARED real_t sQD[vdim*Q*D];
-
-      // Bo and Bc into shared memory
-      const DeviceMatrix Bot(sBot, d-1, q);
-      kernels::internal::LoadB<D,Q>(d-1, q, BO, sBot);
-      const DeviceMatrix Bct(sBct, d, q);
-      kernels::internal::LoadB<D,Q>(d, q, BC, sBct);
-
-      const DeviceCube QQ(sQQ, q, q, vdim);
-      const DeviceCube QD(sQD, q, d, vdim);
-
-      MFEM_FOREACH_THREAD(vd,z,vdim)
-      {
-         const real_t cst_val_0 = C(0,0,0,0);
-         const real_t cst_val_1 = C(1,0,0,0);
-         MFEM_FOREACH_THREAD(y,y,q)
-         {
-            MFEM_FOREACH_THREAD(x,x,q)
-            {
-               const real_t J0 = J(x,y,0,vd,e);
-               const real_t J1 = J(x,y,1,vd,e);
-               const real_t C0 = cst ? cst_val_0 : C(0,x,y,e);
-               const real_t C1 = cst ? cst_val_1 : C(1,x,y,e);
-               QQ(x,y,vd) = W(x,y)*(J0*C0 + J1*C1);
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(vd,z,vdim)
-      {
-         const int nx = (vd == 0) ? d : d-1;
-         DeviceMatrix Btx = (vd == 0) ? Bct : Bot;
-         MFEM_FOREACH_THREAD(qy,y,q)
-         {
-            MFEM_FOREACH_THREAD(dx,x,nx)
-            {
-               real_t qd = 0.0;
-               for (int qx = 0; qx < q; ++qx)
-               {
-                  qd += QQ(qx,qy,vd) * Btx(dx,qx);
-               }
-               QD(dx,qy,vd) = qd;
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(vd,z,vdim)
-      {
-         const int nx = (vd == 0) ? d : d-1;
-         const int ny = (vd == 1) ? d : d-1;
-         DeviceMatrix Bty = (vd == 1) ? Bct : Bot;
-         DeviceTensor<4> Yxy(Y, nx, ny, vdim, ne);
-         MFEM_FOREACH_THREAD(dy,y,ny)
-         {
-            MFEM_FOREACH_THREAD(dx,x,nx)
-            {
-               real_t dd = 0.0;
-               for (int qy = 0; qy < q; ++qy)
-               {
-                  dd += QD(dx,qy,vd) * Bty(dy,qy);
-               }
-               Yxy(dx,dy,vd,e) += dd;
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-   });
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 1, 2>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 2, 3>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 3, 4>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 4, 5>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 5, 6>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 6, 7>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 7, 8>();
+   VectorFEDomainLFIntegrator::AddSpecialization<FiniteElement::CURL, 3, 8, 9>();
 }
 
-template<int T_D1D = 0, int T_Q1D = 0>
-static void HdivDLFAssemble3D(
-   const int ne, const int d, const int q, const int *markers, const real_t *bo,
-   const real_t *bc, const real_t *j, const real_t *weights,
-   const Vector &coeff, real_t *y)
+/// \cond DO_NOT_DOCUMENT
+VectorFEDomainLFIntegrator::AssembleKernelType
+VectorFEDomainLFIntegrator::AssembleKernels::Fallback(
+   FiniteElement::DerivType TestType, int DIM, int, int)
 {
-   MFEM_VERIFY(T_D1D || d <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
-               "Problem size too large.");
-   MFEM_VERIFY(T_Q1D || q <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,
-               "Problem size too large.");
-
-   static constexpr int vdim = 3;
-   const auto F = coeff.Read();
-   const auto M = Reshape(markers, ne);
-   const auto BO = Reshape(bo, q, d-1);
-   const auto BC = Reshape(bc, q, d);
-   const auto J = Reshape(j, q, q, q, vdim, vdim, ne);
-   const auto W = Reshape(weights, q, q, q);
-   const bool cst = coeff.Size() == vdim;
-   const auto C = cst ? Reshape(F,vdim,1,1,1,1) : Reshape(F,vdim,q,q,q,ne);
-   auto Y = Reshape(y, 2*(d-1)*(d-1)*d, ne);
-
-   mfem::forall_3D(ne, q, q, vdim, [=] MFEM_HOST_DEVICE (int e)
+   if (TestType == FiniteElement::DIV)
    {
-      if (M(e) == 0) { return; } // ignore
-
-      constexpr int Q = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
-      constexpr int D = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
-
-      MFEM_SHARED real_t sBot[Q*D];
-      MFEM_SHARED real_t sBct[Q*D];
-
-      // Bo and Bc into shared memory
-      const DeviceMatrix Bot(sBot, d-1, q);
-      kernels::internal::LoadB<D,Q>(d-1, q, BO, sBot);
-      const DeviceMatrix Bct(sBct, d, q);
-      kernels::internal::LoadB<D,Q>(d, q, BC, sBct);
-
-      MFEM_SHARED real_t sm0[vdim*Q*Q*Q];
-      MFEM_SHARED real_t sm1[vdim*Q*Q*Q];
-      DeviceTensor<4> QQQ(sm1, q, q, q, vdim);
-      DeviceTensor<4> DQQ(sm0, d, q, q, vdim);
-      DeviceTensor<4> DDQ(sm1, d, d, q, vdim);
-
-      MFEM_FOREACH_THREAD(vd,z,vdim)
+      if (DIM == 2)
       {
-         const real_t cst_val_0 = C(0,0,0,0,0);
-         const real_t cst_val_1 = C(1,0,0,0,0);
-         const real_t cst_val_2 = C(2,0,0,0,0);
-         MFEM_FOREACH_THREAD(y,y,q)
-         {
-            MFEM_FOREACH_THREAD(x,x,q)
-            {
-               for (int z = 0; z < q; ++z)
-               {
-                  const real_t J0 = J(x,y,z,0,vd,e);
-                  const real_t J1 = J(x,y,z,1,vd,e);
-                  const real_t J2 = J(x,y,z,2,vd,e);
-                  const real_t C0 = cst ? cst_val_0 : C(0,x,y,z,e);
-                  const real_t C1 = cst ? cst_val_1 : C(1,x,y,z,e);
-                  const real_t C2 = cst ? cst_val_2 : C(2,x,y,z,e);
-                  QQQ(x,y,z,vd) = W(x,y,z)*(J0*C0 + J1*C1 + J2*C2);
-               }
-            }
-         }
+         return HdivDLFAssemble2D<0, 0>;
       }
-      MFEM_SYNC_THREAD;
-      // Apply Bt operator
-      MFEM_FOREACH_THREAD(vd,z,vdim)
+      if (DIM == 3)
       {
-         const int nx = (vd == 0) ? d : d-1;
-         DeviceMatrix Btx = (vd == 0) ? Bct : Bot;
-         MFEM_FOREACH_THREAD(qy,y,q)
-         {
-            MFEM_FOREACH_THREAD(dx,x,nx)
-            {
-               real_t u[Q];
-               MFEM_UNROLL(Q)
-               for (int qz = 0; qz < q; ++qz) { u[qz] = 0.0; }
-               MFEM_UNROLL(Q)
-               for (int qx = 0; qx < q; ++qx)
-               {
-                  MFEM_UNROLL(Q)
-                  for (int qz = 0; qz < q; ++qz)
-                  {
-                     u[qz] += QQQ(qx,qy,qz,vd) * Btx(dx,qx);
-                  }
-               }
-               MFEM_UNROLL(Q)
-               for (int qz = 0; qz < q; ++qz) { DQQ(dx,qy,qz,vd) = u[qz]; }
-            }
-         }
+         return HdivDLFAssemble3D<0, 0>;
       }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(vd,z,vdim)
-      {
-         const int nx = (vd == 0) ? d : d-1;
-         const int ny = (vd == 1) ? d : d-1;
-         DeviceMatrix Bty = (vd == 1) ? Bct : Bot;
-         MFEM_FOREACH_THREAD(dy,y,ny)
-         {
-            MFEM_FOREACH_THREAD(dx,x,nx)
-            {
-               real_t u[Q];
-               MFEM_UNROLL(Q)
-               for (int qz = 0; qz < q; ++qz) { u[qz] = 0.0; }
-               MFEM_UNROLL(Q)
-               for (int qy = 0; qy < q; ++qy)
-               {
-                  MFEM_UNROLL(Q)
-                  for (int qz = 0; qz < q; ++qz)
-                  {
-                     u[qz] += DQQ(dx,qy,qz,vd) * Bty(dy,qy);
-                  }
-               }
-               MFEM_UNROLL(Q)
-               for (int qz = 0; qz < q; ++qz) { DDQ(dx,dy,qz,vd) = u[qz]; }
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(vd,z,vdim)
-      {
-         const int nx = (vd == 0) ? d : d-1;
-         const int ny = (vd == 1) ? d : d-1;
-         const int nz = (vd == 2) ? d : d-1;
-         DeviceTensor<5> Yxyz(Y, nx, ny, nz, vdim, ne);
-         DeviceMatrix Btz = (vd == 2) ? Bct : Bot;
-         MFEM_FOREACH_THREAD(dy,y,ny)
-         {
-            MFEM_FOREACH_THREAD(dx,x,nx)
-            {
-               real_t u[D];
-               MFEM_UNROLL(D)
-               for (int dz = 0; dz < nz; ++dz) { u[dz] = 0.0; }
-               MFEM_UNROLL(Q)
-               for (int qz = 0; qz < q; ++qz)
-               {
-                  MFEM_UNROLL(D)
-                  for (int dz = 0; dz < nz; ++dz)
-                  {
-                     u[dz] += DDQ(dx,dy,qz,vd) * Btz(dz,qz);
-                  }
-               }
-               MFEM_UNROLL(D)
-               for (int dz = 0; dz < nz; ++dz) { Yxyz(dx,dy,dz,vd,e) += u[dz]; }
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-   });
-}
-
-static void HdivDLFAssemble(const FiniteElementSpace &fes,
-                            const IntegrationRule *ir,
-                            const Array<int> &markers,
-                            const Vector &coeff,
-                            Vector &y)
-{
-   Mesh &mesh = *fes.GetMesh();
-   const int dim = mesh.Dimension();
-   const FiniteElement *el = fes.GetTypicalFE();
-   const auto *vel = dynamic_cast<const VectorTensorFiniteElement *>(el);
-   MFEM_VERIFY(vel != nullptr, "Must be VectorTensorFiniteElement");
-   const MemoryType mt = Device::GetDeviceMemoryType();
-   const DofToQuad &maps_o = vel->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
-   const DofToQuad &maps_c = vel->GetDofToQuad(*ir, DofToQuad::TENSOR);
-   const int d = maps_c.ndof, q = maps_c.nqpt;
-   constexpr int flags = GeometricFactors::JACOBIANS;
-   const GeometricFactors *geom = mesh.GetGeometricFactors(*ir, flags, mt);
-   decltype(&HdivDLFAssemble2D<>) ker =
-      dim == 2 ? HdivDLFAssemble2D<> : HdivDLFAssemble3D<>;
-
-   if (dim==2)
-   {
-      if (d==1 && q==1) { ker=HdivDLFAssemble2D<1,1>; }
-      if (d==2 && q==2) { ker=HdivDLFAssemble2D<2,2>; }
-      if (d==3 && q==3) { ker=HdivDLFAssemble2D<3,3>; }
-      if (d==4 && q==4) { ker=HdivDLFAssemble2D<4,4>; }
-      if (d==5 && q==5) { ker=HdivDLFAssemble2D<5,5>; }
-      if (d==6 && q==6) { ker=HdivDLFAssemble2D<6,6>; }
-      if (d==7 && q==7) { ker=HdivDLFAssemble2D<7,7>; }
-      if (d==8 && q==8) { ker=HdivDLFAssemble2D<8,8>; }
    }
-
-   if (dim==3)
+   else if (TestType == FiniteElement::CURL)
    {
-      if (d==2 && q==2) { ker=HdivDLFAssemble3D<2,2>; }
-      if (d==3 && q==3) { ker=HdivDLFAssemble3D<3,3>; }
-      if (d==4 && q==4) { ker=HdivDLFAssemble3D<4,4>; }
-      if (d==5 && q==5) { ker=HdivDLFAssemble3D<5,5>; }
-      if (d==6 && q==6) { ker=HdivDLFAssemble3D<6,6>; }
-      if (d==7 && q==7) { ker=HdivDLFAssemble3D<7,7>; }
-      if (d==8 && q==8) { ker=HdivDLFAssemble3D<8,8>; }
+      if (DIM == 3)
+      {
+         return HcurlDLFAssemble3D<0, 0>;
+      }
    }
-
-   MFEM_VERIFY(ker, "No kernel ndof " << d << " nqpt " << q);
-
-   const int ne = mesh.GetNE();
-   const int *M = markers.Read();
-   const real_t *Bo = maps_o.B.Read();
-   const real_t *Bc = maps_c.B.Read();
-   const real_t *J = geom->J.Read();
-   const real_t *W = ir->GetWeights().Read();
-   real_t *Y = y.ReadWrite();
-   ker(ne, d, q, M, Bo, Bc, J, W, coeff, Y);
+   MFEM_ABORT("");
 }
+/// \endcond DO_NOT_DOCUMENT
 
 void VectorFEDomainLFIntegrator::AssembleDevice(const FiniteElementSpace &fes,
                                                 const Array<int> &markers,
@@ -337,15 +96,23 @@ void VectorFEDomainLFIntegrator::AssembleDevice(const FiniteElementSpace &fes,
    QuadratureSpace qs(*fes.GetMesh(), *ir);
    CoefficientVector coeff(QF, qs, CoefficientStorage::COMPRESSED);
 
-   const int fe_type = fe.GetDerivType();
-   if (fe_type == FiniteElement::DIV)
-   {
-      HdivDLFAssemble(fes, ir, markers, coeff, b);
-   }
-   else
-   {
-      MFEM_ABORT("Not implemented.");
-   }
+   const FiniteElement::DerivType fe_type =
+      static_cast<FiniteElement::DerivType>(fe.GetDerivType());
+
+   Mesh &mesh = *fes.GetMesh();
+   const int dim = mesh.Dimension();
+   const FiniteElement *el = fes.GetTypicalFE();
+   const auto *vel = dynamic_cast<const VectorTensorFiniteElement *>(el);
+   MFEM_VERIFY(vel != nullptr, "Must be VectorTensorFiniteElement");
+   const MemoryType mt = Device::GetDeviceMemoryType();
+   const DofToQuad &maps_o = vel->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
+   const DofToQuad &maps_c = vel->GetDofToQuad(*ir, DofToQuad::TENSOR);
+   const int d = maps_c.ndof, q = maps_c.nqpt;
+   constexpr int flags = GeometricFactors::JACOBIANS;
+   const GeometricFactors *geom = mesh.GetGeometricFactors(*ir, flags, mt);
+
+   AssembleKernels::Run(fe_type, dim, d, q, mesh.GetNE(), markers, geom->J,
+                        ir->GetWeights(), maps_o.B, maps_c.B, coeff, b, d, q);
 }
 
 } // namespace mfem

fem/lininteg.cpp

@@ -471,6 +471,13 @@ void VectorBoundaryLFIntegrator::AssembleRHSElementVect(
    }
 }
 
+VectorFEDomainLFIntegrator::VectorFEDomainLFIntegrator(
+   VectorCoefficient &F, const IntegrationRule *ir)
+   : DeltaLFIntegrator(F, ir), QF(F)
+{
+   static Kernels kernels{};
+}
+
 void VectorFEDomainLFIntegrator::AssembleRHSElementVect(
    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {

fem/lininteg.hpp

@@ -369,8 +369,8 @@ private:
    Vector vec;
 
 public:
-   VectorFEDomainLFIntegrator(VectorCoefficient &F)
-      : DeltaLFIntegrator(F), QF(F) { }
+   VectorFEDomainLFIntegrator(VectorCoefficient &F,
+                              const IntegrationRule *ir = nullptr);
 
    void AssembleRHSElementVect(const FiniteElement &el,
                                ElementTransformation &Tr,
@@ -387,6 +387,40 @@ public:
                        Vector &b) override;
 
    using LinearFormIntegrator::AssembleRHSElementVect;
+
+   /// @param ne number of elements
+   /// @param markers array where entry markers[e] == 0 to skip assembly over
+   /// element e element
+   /// @param jac Spatial Jacobians evaluated at all quadrature points
+   /// @param weights 1D quadrature weights
+   /// @param testBO 1D open basis test functions
+   /// @param testBC 1D closed basis test functions
+   /// @param coeff coefficient values evaluated at quadrature points, possibly
+   /// compressed.
+   /// @param d number of 1D closed dofs
+   /// @param q number of 1D quadrature points
+   using AssembleKernelType = void (*)(const int NE, const Array<int> &markers,
+                                       const Vector &jac,
+                                       const Array<real_t> &weights,
+                                       const Array<real_t> &testBO,
+                                       const Array<real_t> &testBC,
+                                       const Vector &coeff, Vector &y,
+                                       const int testd1d, const int q1d);
+
+   /// parameters: test_fetype, ndims, test_d1d, q1d
+   MFEM_REGISTER_KERNELS(AssembleKernels, AssembleKernelType,
+                         (FiniteElement::DerivType, int, int, int));
+
+   struct Kernels
+   {
+      Kernels();
+   };
+
+   template <FiniteElement::DerivType TestType, int DIM, int TEST_D1D, int Q1D>
+   static void AddSpecialization()
+   {
+      AssembleKernels::Specialization<TestType, DIM, TEST_D1D, Q1D>::Add();
+   }
 };
 
 /// $ (Q, \mathrm{curl}(v))_{\Omega} $ for Nedelec Elements

fem/lor/lor_dg_impl.hpp

@@ -258,8 +258,8 @@ void BatchedLOR_DG::Assemble2D()
                   const int w_idx = (n_idx == 0) ? iy : ix;
                   const int x_idx = (n_idx == 0) ? i_0 : j_0;
 
-                  const real_t J1 = J(ix, iy, n_idx, !n_idx, iel_ho);
-                  const real_t J2 = J(ix, iy, !n_idx, !n_idx, iel_ho);
+                  const real_t J1 = J(ix, iy, n_idx, (int)!n_idx, iel_ho);
+                  const real_t J2 = J(ix, iy, (int)!n_idx, (int)!n_idx, iel_ho);
                   const real_t Jh = (J1*J1 + J2*J2) / detJ(ix, iy, iel_ho);
 
                   V(v_idx, ix, iy, iel_ho) =

fem/qinterp/eval_hdiv.hpp

@@ -136,8 +136,8 @@ inline void EvalHDiv2D(const int NE,
                {
                   qq += QD(qx,dy,vd) * By(dy,qy);
                }
-               if (FLAGS & (QuadratureInterpolator::PHYSICAL_VALUES |
-                            QuadratureInterpolator::PHYSICAL_MAGNITUDES))
+               if ((FLAGS & (QuadratureInterpolator::PHYSICAL_VALUES |
+                             QuadratureInterpolator::PHYSICAL_MAGNITUDES)) != 0)
                {
                   QQ(qx,qy,vd) = qq;
                }
@@ -153,8 +153,8 @@ inline void EvalHDiv2D(const int NE,
          }
       }
       MFEM_SYNC_THREAD;
-      if (FLAGS & (QuadratureInterpolator::PHYSICAL_VALUES |
-                   QuadratureInterpolator::PHYSICAL_MAGNITUDES))
+      if ((FLAGS & (QuadratureInterpolator::PHYSICAL_VALUES |
+                    QuadratureInterpolator::PHYSICAL_MAGNITUDES)) != 0)
       {
          if (tidz == 0)
          {
@@ -178,7 +178,7 @@ inline void EvalHDiv2D(const int NE,
                   const real_t detJ = kernels::Det<DIM>(J_loc);
                   kernels::Mult(DIM, DIM, J_loc, u_ref, u_phys);
                   kernels::Set(DIM, 1, 1_r/detJ, u_phys, u_phys);
-                  if (FLAGS & QuadratureInterpolator::PHYSICAL_VALUES)
+                  if ((FLAGS & QuadratureInterpolator::PHYSICAL_VALUES) != 0)
                   {
                      MFEM_UNROLL(DIM)
                      for (int sd = 0; sd < DIM; sd++)
@@ -193,7 +193,8 @@ inline void EvalHDiv2D(const int NE,
                         }
                      }
                   }
-                  else if (FLAGS & QuadratureInterpolator::PHYSICAL_MAGNITUDES)
+                  else if ((FLAGS &
+                            QuadratureInterpolator::PHYSICAL_MAGNITUDES) != 0)
                   {
                      y(qx,qy,0,e) = kernels::Norml2(DIM, u_phys);
                   }
@@ -372,8 +373,9 @@ inline void EvalHDiv3D(const int NE,
                MFEM_UNROLL(MQ1)
                for (int qz = 0; qz < Q1D; ++qz)
                {
-                  if (FLAGS & (QuadratureInterpolator::PHYSICAL_VALUES |
-                               QuadratureInterpolator::PHYSICAL_MAGNITUDES))
+                  if ((FLAGS & (QuadratureInterpolator::PHYSICAL_VALUES |
+                                QuadratureInterpolator::PHYSICAL_MAGNITUDES)) !=
+                      0)
                   {
                      QQQ(qx,qy,qz,vd) = u[qz];
                   }
@@ -390,8 +392,8 @@ inline void EvalHDiv3D(const int NE,
          }
       }
       MFEM_SYNC_THREAD;
-      if (FLAGS & (QuadratureInterpolator::PHYSICAL_VALUES |
-                   QuadratureInterpolator::PHYSICAL_MAGNITUDES))
+      if ((FLAGS & (QuadratureInterpolator::PHYSICAL_VALUES |
+                    QuadratureInterpolator::PHYSICAL_MAGNITUDES)) != 0)
       {
          MFEM_FOREACH_THREAD(qz,z,Q1D)
          {
@@ -415,7 +417,7 @@ inline void EvalHDiv3D(const int NE,
                   const real_t detJ = kernels::Det<DIM>(J_loc);
                   kernels::Mult(DIM, DIM, J_loc, u_ref, u_phys);
                   kernels::Set(DIM, 1, 1_r/detJ, u_phys, u_phys);
-                  if (FLAGS & QuadratureInterpolator::PHYSICAL_VALUES)
+                  if ((FLAGS & QuadratureInterpolator::PHYSICAL_VALUES) != 0)
                   {
                      MFEM_UNROLL(DIM)
                      for (int sd = 0; sd < DIM; sd++)
@@ -430,7 +432,8 @@ inline void EvalHDiv3D(const int NE,
                         }
                      }
                   }
-                  else if (FLAGS & QuadratureInterpolator::PHYSICAL_MAGNITUDES)
+                  else if ((FLAGS &
+                            QuadratureInterpolator::PHYSICAL_MAGNITUDES) != 0)
                   {
                      y(qx,qy,qz,0,e) = kernels::Norml2(DIM, u_phys);
                   }

general/backends.hpp

@@ -75,6 +75,16 @@ constexpr bool mfem_use_gpu = false;
 #define MFEM_THREAD_SIZE(k) 1
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=0; i<N; i++)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) MFEM_FOREACH_THREAD(i,k,N)
+// Assigns a thread block shaped (SX,SY,SZ) contiguous in x.
+#define MFEM_FOREACH_THREAD_DIRECT_3D(ix, iy, iz, k, SX, SY, SZ)               \
+   for (int iz = 0; iz < SZ; ++iz)                                             \
+      for (int iy = 0; iy < SY; ++iy)                                          \
+         for (int ix = 0; ix < SX; ++ix)
+// Assigns a thread block shaped (OX,OY,OZ) to work on items (SX,SY,SZ),
+// contiguous in x. This intentionally offsets threads
+#define MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(ix, iy, iz, k, SX, SY, SZ, OX,    \
+                                             OY, OZ)                           \
+   MFEM_FOREACH_THREAD_DIRECT_3D(ix, iy, iz, k, SX, SY, SZ)
 #endif
 
 // 'double' and 'float' atomicAdd implementation for previous versions of CUDA

general/cuda.hpp

@@ -49,6 +49,16 @@ constexpr bool mfem_use_gpu = true;
 #define MFEM_THREAD_SIZE(k) blockDim.k
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=threadIdx.k; i<N; i+=blockDim.k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) if(const int i=threadIdx.k; i<N)
+// Assigns a thread block shaped (SX,SY,SZ) contiguous in x.
+#define MFEM_FOREACH_THREAD_DIRECT_3D(ix, iy, iz, k, SX, SY, SZ)               \
+   if (int ix = threadIdx.k % (SX), iy = threadIdx.k / (SX), iz = iy / (SY);   \
+       (iy %= (SY)), (threadIdx.k < (SX) * (SY) * (SZ)))
+// Assigns a thread block shaped (OX,OY,OZ) to work on items (SX,SY,SZ),
+// contiguous in x. This intentionally offsets threads
+#define MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(ix, iy, iz, k, SX, SY, SZ, OX,    \
+                                             OY, OZ)                           \
+   if (int ix = threadIdx.k % (OX), iy = threadIdx.k / (OX), iz = iy / (OY);   \
+       (ix < (SX)) && ((iy %= (OY)) < (SY)) && (iz < (SZ)))
 #endif // defined(__CUDA_ARCH__)
 #endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 

general/error.cpp

@@ -171,7 +171,14 @@ void mfem_error(const char *msg)
 #ifdef MFEM_USE_EXCEPTIONS
    if (mfem_error_action == MFEM_ERROR_THROW)
    {
-      throw ErrorException(msg);
+      if (msg)
+      {
+         throw ErrorException(msg);
+      }
+      else
+      {
+         throw ErrorException("");
+      }
    }
 #endif
 

general/hip.hpp

@@ -51,6 +51,18 @@ constexpr bool mfem_use_gpu = true;
    for(int i=hipThreadIdx_ ##k; i<N; i+=hipBlockDim_ ##k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) \
    if(const int i=hipThreadIdx_ ##k; i<N)
+// Assigns a thread block shaped (SX,SY,SZ) contiguous in x.
+#define MFEM_FOREACH_THREAD_DIRECT_3D(ix, iy, iz, k, SX, SY, SZ)               \
+   if (int ix = hipThreadIdx_##k % (SX), iy = hipThreadIdx_##k / (SX),         \
+       iz = iy / (SY);                                                         \
+       (iy %= (SY)), (hipThreadIdx_##k < (SX) * (SY) * (SZ)))
+// Assigns a thread block shaped (OX,OY,OZ) to work on items (SX,SY,SZ),
+// contiguous in x. This intentionally offsets threads
+#define MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(ix, iy, iz, k, SX, SY, SZ, OX,    \
+                                             OY, OZ)                           \
+   if (int ix = hipThreadIdx_##k % (OX), iy = hipThreadIdx_##k / (OX),         \
+       iz = iy / (OY);                                                         \
+       (ix < (SX)) && ((iy %= (OY)) < (SY)) && (iz < (SZ)))
 #endif // defined(__HIP_DEVICE_COMPILE__)
 #endif // defined(MFEM_USE_HIP) && defined(__HIP__)
 

general/optparser.cpp

@@ -148,9 +148,26 @@ void parseVector(char * str, Vector & var)
    }
 }
 
+bool parseEnumOption(const char *str, size_t &var,
+                     const std::vector<std::string> &options)
+{
+   for (size_t i = 0; i < options.size(); ++i)
+   {
+      auto &v = options[i];
+      if (v == str)
+      {
+         var = i;
+         // success
+         return true;
+      }
+   }
+   // failure
+   return false;
+}
+
 void OptionsParser::Parse()
 {
-   option_check.SetSize(options.Size());
+   option_check.SetSize(options.size());
    option_check = 0;
    for (int i = 1; i < argc; )
    {
@@ -161,9 +178,9 @@ void OptionsParser::Parse()
          return;
       }
 
-      for (int j = 0; true; j++)
+      for (size_t j = 0; true; j++)
       {
-         if (j >= options.Size())
+         if (j >= options.size())
          {
             // unrecognized option
             error_type = 2;
@@ -171,8 +188,8 @@ void OptionsParser::Parse()
             return;
          }
 
-         if (strcmp(argv[i], options[j].short_name) == 0 ||
-             strcmp(argv[i], options[j].long_name) == 0)
+         if (strcmp(argv[i], options[j].short_name.c_str()) == 0 ||
+             strcmp(argv[i], options[j].long_name.c_str()) == 0)
          {
             OptionType type = options[j].type;
 
@@ -224,6 +241,11 @@ void OptionsParser::Parse()
                case VECTOR:
                   parseVector(argv[i++], *(Vector*)(options[j].var_ptr) );
                   break;
+               case ENUM_OPTION:
+                  isValid =
+                     parseEnumOption(argv[i++], *(size_t *)(options[j].var_ptr),
+                                     options[j].options);
+                  break;
             }
 
             if (!isValid)
@@ -239,7 +261,7 @@ void OptionsParser::Parse()
    }
 
    // check for missing required options
-   for (int i = 0; i < options.Size(); i++)
+   for (size_t i = 0; i < options.size(); i++)
       if (options[i].required &&
           (option_check[i] == 0 ||
            (options[i].type == ENABLE && option_check[++i] == 0)))
@@ -323,6 +345,12 @@ void OptionsParser::WriteValue(const Option &opt, std::ostream &os)
          break;
       }
 
+      case ENUM_OPTION:
+      {
+         os << opt.options.at(*((size_t *)opt.var_ptr));
+         break;
+      }
+
       default: // provide a default to suppress warning
          break;
    }
@@ -333,7 +361,7 @@ void OptionsParser::PrintOptions(ostream &os) const
    static const char *indent = "   ";
 
    os << "Options used:\n";
-   for (int j = 0; j < options.Size(); j++)
+   for (size_t j = 0; j < options.size(); j++)
    {
       OptionType type = options[j].type;
 
@@ -355,7 +383,7 @@ void OptionsParser::PrintOptions(ostream &os) const
          os << options[j].long_name << " ";
          WriteValue(options[j], os);
       }
-      os << '\n';
+      os << std::endl;
    }
 }
 
@@ -410,12 +438,13 @@ void OptionsParser::PrintHelp(ostream &os) const
    static const char *line_sep = "";
    static const char *types[] = { " <int>", " <double>", " <string>",
                                   " <string>", "", "", " '<int>...'",
-                                  " '<double>...'"
+                                  " '<double>...'",
+                                  " <string>",
                                 };
 
    os << indent << "-h" << seprtr << "--help" << descr_sep
       << "Print this help message and exit.\n" << line_sep;
-   for (int j = 0; j < options.Size(); j++)
+   for (size_t j = 0; j < options.size(); j++)
    {
       OptionType type = options[j].type;
 
@@ -451,9 +480,22 @@ void OptionsParser::PrintHelp(ostream &os) const
       }
       os << descr_sep;
 
-      if (options[j].description)
+      if (options[j].description.size())
       {
-         os << options[j].description << '\n';
+         os << options[j].description << std::endl;
+      }
+      if (options[j].type == ENUM_OPTION)
+      {
+         os << "\tavailable choices: ";
+         for (size_t i = 0; i < options[j].options.size(); ++i)
+         {
+            os << '"' << options[j].options[i] << '"';
+            if (i + 1 < options[j].options.size())
+            {
+               os << ", ";
+            }
+         }
+         os << std::endl;
       }
       os << line_sep;
    }

general/optparser.hpp

@@ -15,6 +15,9 @@
 #include "../config/config.hpp"
 #include "array.hpp"
 
+#include <string>
+#include <vector>
+
 namespace mfem
 {
 
@@ -31,17 +34,18 @@ class Vector;
 class OptionsParser
 {
 public:
-   enum OptionType { INT, DOUBLE, STRING, STD_STRING, ENABLE, DISABLE, ARRAY, VECTOR };
+   enum OptionType { INT, DOUBLE, STRING, STD_STRING, ENABLE, DISABLE, ARRAY, VECTOR, ENUM_OPTION };
 
 private:
    struct Option
    {
       OptionType type;
       void *var_ptr;
-      const char *short_name;
-      const char *long_name;
-      const char *description;
+      std::string short_name;
+      std::string long_name;
+      std::string description;
       bool required;
+      std::vector<std::string> options;
 
       Option() = default;
 
@@ -49,11 +53,19 @@ private:
              const char *long_name_, const char *description_, bool req)
          : type(type_), var_ptr(var_ptr_), short_name(short_name_),
            long_name(long_name_), description(description_), required(req) { }
+
+      Option(size_t *var_ptr_, const char *short_name_, const char *long_name_,
+             const char *description_, std::vector<std::string> &&options_,
+             bool req)
+         : type(ENUM_OPTION), var_ptr(var_ptr_), short_name(short_name_),
+           long_name(long_name_), description(description_),
+           required(req), options(std::move(options_))
+      {}
    };
 
    int argc;
    char **argv;
-   Array<Option> options;
+   std::vector<Option> options;
    Array<int> option_check;
    // error_type can be:
    //  0 - no error
@@ -84,26 +96,26 @@ public:
                   const char *disable_long_name, const char *description,
                   bool required = false)
    {
-      options.Append(Option(ENABLE, var, enable_short_name, enable_long_name,
-                            description, required));
-      options.Append(Option(DISABLE, var, disable_short_name, disable_long_name,
-                            description, required));
+      options.emplace_back(ENABLE, var, enable_short_name, enable_long_name,
+                           description, required);
+      options.emplace_back(DISABLE, var, disable_short_name, disable_long_name,
+                           description, required);
    }
 
    /// Add an integer option and set 'var' to receive the value.
    void AddOption(int *var, const char *short_name, const char *long_name,
                   const char *description, bool required = false)
    {
-      options.Append(Option(INT, var, short_name, long_name, description,
-                            required));
+      options.emplace_back(INT, var, short_name, long_name, description,
+                           required);
    }
 
    /// Add a double option and set 'var' to receive the value.
    void AddOption(real_t *var, const char *short_name, const char *long_name,
                   const char *description, bool required = false)
    {
-      options.Append(Option(DOUBLE, var, short_name, long_name, description,
-                            required));
+      options.emplace_back(DOUBLE, var, short_name, long_name, description,
+                           required);
    }
 
    /// Add a string (char*) option and set 'var' to receive the value.
@@ -111,8 +123,8 @@ public:
                   const char *long_name, const char *description,
                   bool required = false)
    {
-      options.Append(Option(STRING, var, short_name, long_name, description,
-                            required));
+      options.emplace_back(STRING, var, short_name, long_name, description,
+                           required);
    }
 
    /// Add a string (std::string) option and set 'var' to receive the value.
@@ -120,8 +132,8 @@ public:
                   const char *long_name, const char *description,
                   bool required = false)
    {
-      options.Append(Option(STD_STRING, var, short_name, long_name, description,
-                            required));
+      options.emplace_back(STD_STRING, var, short_name, long_name, description,
+                           required);
    }
 
    /** Add an integer array (separated by spaces) option and set 'var' to
@@ -130,8 +142,8 @@ public:
                   const char *long_name, const char *description,
                   bool required = false)
    {
-      options.Append(Option(ARRAY, var, short_name, long_name, description,
-                            required));
+      options.emplace_back(ARRAY, var, short_name, long_name, description,
+                           required);
    }
 
    /** Add a vector (doubles separated by spaces) option and set 'var' to
@@ -140,8 +152,20 @@ public:
                   const char *long_name, const char *description,
                   bool required = false)
    {
-      options.Append(Option(VECTOR, var, short_name, long_name, description,
-                            required));
+      options.emplace_back(VECTOR, var, short_name, long_name, description,
+                           required);
+   }
+
+   /**
+    * Add an option which must be one of a given list of options
+    */
+   void AddOptionChoice(size_t *index, const char *short_name,
+                        const char *long_name, const char *description,
+                        std::vector<std::string> options_,
+                        bool required = false)
+   {
+      options.emplace_back(index, short_name, long_name, description,
+                           std::move(options_), required);
    }
 
    /** @brief Parse the command-line options.

miniapps/electromagnetics/CMakeLists.txt

@@ -34,6 +34,11 @@ if (MFEM_USE_MPI)
    EXTRA_HEADERS maxwell_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
    LIBRARIES mfem-common)
 
+ add_mfem_miniapp(maxwell-gpu
+   MAIN maxwell_gpu.cpp
+   EXTRA_HEADERS ${MFEM_MINIAPPS_COMMON_HEADERS}
+   LIBRARIES mfem-common)
+
  if (MFEM_USE_GSLIB)
   add_mfem_miniapp(lorentz
     MAIN lorentz.cpp

miniapps/electromagnetics/electromagnetics.hpp

@@ -12,6 +12,8 @@
 #ifndef MFEM_ELECTROMAGNETICS_HPP
 #define MFEM_ELECTROMAGNETICS_HPP
 
+#include <cmath>
+
 namespace mfem
 {
 

miniapps/electromagnetics/maxwell_solver.cpp

@@ -308,7 +308,7 @@ MaxwellSolver::~MaxwellSolver()
    {
       int i = mit1->first;
       delete pcg_[i];
-      delete diagScale_[i];
+      // delete diagScale_[i];
       delete A1_[i];
       delete a1_[i];
    }
@@ -377,7 +377,7 @@ MaxwellSolver::Mult(const Vector &B, Vector &dEdt) const
 void
 MaxwellSolver::ImplicitSolve(real_t dt, const Vector &B, Vector &dEdt)
 {
-   implicitSolve(dt, B, dEdt);
+   const_cast<const MaxwellSolver *>(this)->implicitSolve(dt, B, dEdt);
 }
 
 void
@@ -420,12 +420,12 @@ MaxwellSolver::setupSolver(const int idt, const real_t dt) const
       a1_[idt]->Finalize();
       A1_[idt] = a1_[idt]->ParallelAssemble();
 
-      diagScale_[idt] = new HypreDiagScale(*A1_[idt]);
+      // diagScale_[idt] = new HypreDiagScale(*A1_[idt]);
       pcg_[idt] = new HyprePCG(*A1_[idt]);
       pcg_[idt]->SetTol(1.0e-12);
       pcg_[idt]->SetMaxIter(200);
       pcg_[idt]->SetPrintLevel(0);
-      pcg_[idt]->SetPreconditioner(*diagScale_[idt]);
+      // pcg_[idt]->SetPreconditioner(*diagScale_[idt]);
    }
 }
 

tests/unit/fem/test_3d_bilininteg.cpp

@@ -5920,4 +5920,206 @@ TEST_CASE("3D Bilinear Div Div Integrators",
    }
 }
 
+TEST_CASE("3D Bilinear VectorFE Integrators PartialAssembly",
+          "[BilinearFormIntegrator]"
+          "[GPU]")
+{
+   auto order = GENERATE(1, 3);
+   CAPTURE(order);
+   int dim = 3;
+
+   FunctionCoefficient q3_coeff(q3);
+   VectorFunctionCoefficient F3_coeff(dim, F3);
+   MatrixFunctionCoefficient M3_coeff(dim, M3);
+
+   auto mesh_fname =
+      GENERATE("../../data/fichera-amr.mesh", "../../data/ball-nurbs.mesh");
+   CAPTURE(mesh_fname);
+   Mesh mesh(mesh_fname);
+   REQUIRE(mesh.Dimension() == dim);
+   REQUIRE(mesh.SpaceDimension() == dim);
+
+   // convert nurbs into piecewise-quadratic curved mesh
+   if (mesh.NURBSext)
+   {
+      mesh.UniformRefinement();
+      mesh.SetCurvature(2);
+   }
+
+   SECTION("H(curl) H(curl) Scalar Coeff")
+   {
+      ND_FECollection fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      BilinearForm bfa(&fespace_nd);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(q3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      BilinearForm bpa(&fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(q3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_nd), y_fa(&fespace_nd), y_pa(&fespace_nd);
+      x.Randomize(1234);
+      bfa.Mult(x, y_fa);
+      bpa.Mult(x, y_pa);
+      y_pa -= y_fa;
+      REQUIRE( y_pa.Normlinf() == MFEM_Approx(0_r) );
+   }
+
+   SECTION("H(curl) H(curl) Diagonal Matrix Coeff")
+   {
+      ND_FECollection fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      BilinearForm bfa(&fespace_nd);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(F3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      BilinearForm bpa(&fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(F3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_nd), y_fa(&fespace_nd), y_pa(&fespace_nd);
+      x.Randomize(1234);
+      bfa.Mult(x, y_fa);
+      bpa.Mult(x, y_pa);
+      y_pa -= y_fa;
+      REQUIRE( y_pa.Normlinf() == MFEM_Approx(0_r) );
+   }
+
+   SECTION("H(curl) H(curl) Matrix Coeff")
+   {
+      ND_FECollection fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      BilinearForm bfa(&fespace_nd);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      BilinearForm bpa(&fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_nd), y_fa(&fespace_nd), y_pa(&fespace_nd);
+      x.Randomize(1234);
+      bfa.Mult(x, y_fa);
+      bpa.Mult(x, y_pa);
+      y_pa -= y_fa;
+      REQUIRE( y_pa.Normlinf() == MFEM_Approx(0_r) );
+   }
+}
+
+TEST_CASE("3D Bilinear Weak Curl Integrators Partial Assembly",
+          "[MixedVectorWeakCurlIntegrator]"
+          "[BilinearFormIntegrator]"
+          "[GPU]")
+{
+   auto order = GENERATE(1, 3);
+   CAPTURE(order);
+   int dim = 3;
+
+   FunctionCoefficient q3_coeff(q3);
+   VectorFunctionCoefficient F3_coeff(dim, F3);
+
+   auto mesh_fname =
+      GENERATE("../../data/fichera-amr.mesh", "../../data/ball-nurbs.mesh");
+   CAPTURE(mesh_fname);
+   Mesh mesh(mesh_fname);
+   REQUIRE(mesh.Dimension() == dim);
+   REQUIRE(mesh.SpaceDimension() == dim);
+
+   // convert nurbs into piecewise-quadratic curved mesh
+   if (mesh.NURBSext)
+   {
+      mesh.UniformRefinement();
+      mesh.SetCurvature(2);
+   }
+
+   SECTION("H(div) H(curl) No Coeff")
+   {
+      ND_FECollection fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+      RT_FECollection fec_rt(order - 1, dim);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+      MixedBilinearForm bfa(&fespace_rt, &fespace_nd);
+      bfa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator);
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_rt, &fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator);
+      bpa.Assemble();
+
+      GridFunction x(&fespace_rt), y_fa(&fespace_nd), y_pa(&fespace_nd);
+      x.Randomize(1234);
+      REQUIRE(bfa.Height() == y_fa.Size());
+      REQUIRE(bfa.Width() == x.Size());
+      REQUIRE(bpa.Height() == y_fa.Size());
+      REQUIRE(bpa.Width() == x.Size());
+      bfa.Mult(x, y_fa);
+      bpa.Mult(x, y_pa);
+      y_pa -= y_fa;
+      REQUIRE( y_pa.Normlinf() == MFEM_Approx(0_r) );
+   }
+
+   SECTION("H(div) H(curl) Scalar Coeff")
+   {
+      ND_FECollection fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+      RT_FECollection fec_rt(order - 1, dim);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+      MixedBilinearForm bfa(&fespace_rt, &fespace_nd);
+      bfa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(q3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_rt, &fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(q3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_rt), y_fa(&fespace_nd), y_pa(&fespace_nd);
+      x.Randomize(1234);
+      bfa.Mult(x, y_fa);
+      bpa.Mult(x, y_pa);
+      y_pa -= y_fa;
+      REQUIRE( y_pa.Normlinf() == MFEM_Approx(0_r) );
+   }
+
+   SECTION("H(div) H(curl) Diagonal Matrix Coeff")
+   {
+      ND_FECollection fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+      RT_FECollection fec_rt(order - 1, dim);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+      MixedBilinearForm bfa(&fespace_rt, &fespace_nd);
+      bfa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(F3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_rt, &fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(F3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_rt), y_fa(&fespace_nd), y_pa(&fespace_nd);
+      x.Randomize(1234);
+      bfa.Mult(x, y_fa);
+      bpa.Mult(x, y_pa);
+      y_pa -= y_fa;
+      REQUIRE( y_pa.Normlinf() == MFEM_Approx(0_r) );
+   }
+}
+
 } // namespace bilininteg_3d

tests/unit/fem/test_lin_interp.cpp

@@ -1069,4 +1069,75 @@ TEST_CASE("Exact Sequence Properties: d(df)=0",
    }
 }
 
+TEST_CASE("Partial Assemble Linear Interpolator",
+          "[CurlInterpolator]"
+          "[GPU]")
+{
+   const int maxOrder = 3;
+   auto order = GENERATE_COPY(range(1, maxOrder + 1));
+   CAPTURE(order);
+
+   int n = 3, dim = -1;
+   real_t tol = 1e-10;
+
+   auto type = Element::HEXAHEDRON;
+   CAPTURE(type);
+
+   Mesh mesh;
+
+   {
+      dim = 3;
+      mesh = Mesh::MakeCartesian3D(n, n, n, (Element::Type)type,
+                                   2.0, 3.0, 5.0);
+   }
+
+   ND_FECollection    fec_nd(order, dim);
+   RT_FECollection    fec_rt(order - 1, dim);
+
+   FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+   FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+   // 3D
+   {
+      DiscreteLinearOperator CurlFA(&fespace_nd, &fespace_rt);
+      CurlFA.AddDomainInterpolator(new CurlInterpolator());
+      CurlFA.Assemble();
+      CurlFA.Finalize();
+      DiscreteLinearOperator CurlPA(&fespace_nd, &fespace_rt);
+      CurlPA.AddDomainInterpolator(new CurlInterpolator());
+      CurlPA.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      CurlPA.Assemble();
+
+      SECTION("Curl (3D)")
+      {
+         SparseMatrix &Curl = CurlFA.SpMat();
+         GridFunction x(&fespace_nd), y_fa(&fespace_rt), y_pa(&fespace_rt);
+         VectorFunctionCoefficient coeff(3, [](const Vector &x, Vector &y)
+         {
+            y.SetSize(3);
+            y[0] = sin(2 * M_PI * x[2] / 5) - cos(2 * M_PI * x[1] / 3);
+            y[1] = sin(2 * M_PI * x[0] / 2) - cos(2 * M_PI * x[2] / 5);
+            y[2] = sin(2 * M_PI * x[1] / 3) - cos(2 * M_PI * x[0] / 2);
+         });
+         x.ProjectCoefficient(coeff);
+         REQUIRE(x.Size() == Curl.Width());
+         REQUIRE(y_fa.Size() == Curl.Height());
+         REQUIRE(x.Size() == CurlPA.Width());
+         REQUIRE(y_pa.Size() == CurlPA.Height());
+         Curl.Mult(x, y_fa);
+         CurlPA.Mult(x, y_pa);
+         y_pa -= y_fa;
+         REQUIRE(y_pa.Normlinf() < tol);
+
+         // transpose
+         y_fa.ProjectCoefficient(coeff);
+         GridFunction x_fa(&fespace_nd), x_pa(&fespace_nd);
+         Curl.MultTranspose(y_fa, x_fa);
+         CurlPA.MultTranspose(y_fa, x_pa);
+         x_pa -= x_fa;
+         REQUIRE(x_pa.Normlinf() < tol);
+      }
+   }
+}
+
 } // namespace lin_interp

tests/unit/fem/test_linearform_ext.cpp

@@ -328,7 +328,7 @@ TEST_CASE("Linear Form Extension", "[LinearFormExtension], [GPU]")
    }
 }
 
-TEST_CASE("H(div) Linear Form Extension", "[LinearFormExtension], [GPU]")
+TEST_CASE("Vector FE Linear Form Extension", "[LinearFormExtension], [GPU]")
 {
    const bool all = launch_all_non_regression_tests;
 
@@ -341,26 +341,44 @@ TEST_CASE("H(div) Linear Form Extension", "[LinearFormExtension], [GPU]")
    Mesh mesh(mesh_file);
    const int dim = mesh.Dimension();
 
-   CAPTURE(mesh_file, dim, p);
+   {
+      const auto space_type =
+         dim == 3 ? GENERATE(FiniteElement::DIV, FiniteElement::CURL)
+         : FiniteElement::DIV;
 
-   RT_FECollection fec(p, dim);
-   FiniteElementSpace fes(&mesh, &fec);
+      CAPTURE(mesh_file, dim, p, space_type);
 
-   VectorFunctionCoefficient coeff(dim, fvec_dim);
+      std::unique_ptr<FiniteElementCollection> fec;
 
-   LinearForm d1(&fes);
-   d1.AddDomainIntegrator(new VectorFEDomainLFIntegrator(coeff));
-   d1.UseFastAssembly(true);
-   d1.Assemble();
+      switch (space_type)
+      {
+         case FiniteElement::DIV:
+            fec.reset(new RT_FECollection(p, dim));
+            break;
+         case FiniteElement::CURL:
+            fec.reset(new ND_FECollection(p, dim));
+            break;
+         default:
+            MFEM_ABORT("unsupported space type");
+      }
+      FiniteElementSpace fes(&mesh, fec.get());
 
-   LinearForm d2(&fes);
-   d2.AddDomainIntegrator(new VectorFEDomainLFIntegrator(coeff));
-   d2.UseFastAssembly(false);
-   d2.Assemble();
+      VectorFunctionCoefficient coeff(dim, fvec_dim);
 
-   CAPTURE(d1.Norml2(), d2.Norml2());
-   d1 -= d2;
-   REQUIRE(d1.Norml2() == MFEM_Approx(0.0));
+      LinearForm d1(&fes);
+      d1.AddDomainIntegrator(new VectorFEDomainLFIntegrator(coeff));
+      d1.UseFastAssembly(true);
+      d1.Assemble();
+
+      LinearForm d2(&fes);
+      d2.AddDomainIntegrator(new VectorFEDomainLFIntegrator(coeff));
+      d2.UseFastAssembly(false);
+      d2.Assemble();
+
+      CAPTURE(d1.Norml2(), d2.Norml2());
+      d1 -= d2;
+      REQUIRE(d1.Norml2() == MFEM_Approx(0.0));
+   }
 }
 
 #ifdef MFEM_USE_MPI

tests/unit/fem/test_project_bdr.cpp

@@ -124,8 +124,8 @@ TEST_CASE("3D ProjectBdrCoefficientNormal Scalar",
 
    const double tol = 1e-6;
 
-   const char bdrs_axis[] = {2, 1, 0, 1, 0, 2};
-   const char bdrs_sign[] = {-1, -1, +1, +1, -1, +1};
+   const int bdrs_axis[] = {2, 1, 0, 1, 0, 2};
+   const int bdrs_sign[] = {-1, -1, +1, +1, -1, +1};
 
    for (int type = (int)Element::TETRAHEDRON;
         type <= (int)Element::HEXAHEDRON; type++)