Suppress Doxygen warnings

Introduce some workarounds for MSVC:

- Need to expand variadic macros explicitly
- Adjacent parameter packs in function templates not supported; use nested
  structs instead
- Can't use default parameters in template definition if not present in
  declaration

Some improvements to the kernel dispatch:

- Run is now a static member function (don't need to call Get)
- DIM is not treated differently from the other parameters

.gitignore

@@ -15,6 +15,9 @@
 CMakeCache.txt
 CMakeFiles/
 
+# Clangd server cache
+*.cache*
+
 # Backup files
 *~
 

fem/CMakeLists.txt

@@ -112,8 +112,6 @@ set(SRCS
   qinterp/eval_by_vdim.cpp
   qinterp/grad_by_nodes.cpp
   qinterp/grad_by_vdim.cpp
-  qinterp/grad_phys_by_nodes.cpp
-  qinterp/grad_phys_by_vdim.cpp
   qspace.cpp
   quadinterpolator.cpp
   quadinterpolator_face.cpp
@@ -192,6 +190,9 @@ set(HDRS
   hybridization.hpp
   intrules.hpp
   intrules_cut.hpp
+  kernel_dispatch.hpp
+  kernel_reporter.hpp
+  kernels.hpp
   ceed/interface/basis.hpp
   ceed/interface/integrator.hpp
   ceed/interface/interface.hpp
@@ -223,7 +224,6 @@ set(HDRS
   nonlinearform_ext.hpp
   nonlininteg.hpp
   qfunction.hpp
-  qinterp/dispatch.hpp
   qinterp/eval.hpp
   qinterp/grad.hpp
   qspace.hpp

fem/bilininteg.hpp

@@ -19,6 +19,8 @@
 #include "qfunction.hpp"
 #include <memory>
 
+#include "kernel_dispatch.hpp"
+
 namespace mfem
 {
 
@@ -2127,6 +2129,22 @@ public:
     can be a scalar or a matrix coefficient. */
 class DiffusionIntegrator: public BilinearFormIntegrator
 {
+public:
+
+   using ApplyKernelType = void(*)(const int, const bool, const Array<real_t>&,
+                                   const Array<real_t>&, const Array<real_t>&,
+                                   const Array<real_t>&,
+                                   const Vector&, const Vector&,
+                                   Vector&, const int, const int);
+
+   using DiagonalKernelType = void(*)(const int, const bool, const Array<real_t>&,
+                                      const Array<real_t>&, const Vector&, Vector&,
+                                      const int, const int);
+
+   MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType, (int, int, int));
+   MFEM_REGISTER_KERNELS(DiagonalPAKernels, DiagonalKernelType, (int, int, int));
+   static struct Kernels { Kernels(); } kernels;
+
 protected:
    Coefficient *Q;
    VectorCoefficient *VQ;
@@ -2287,6 +2305,13 @@ public:
    bool SupportsCeed() const { return DeviceCanUseCeed(); }
 
    Coefficient *GetCoefficient() const { return Q; }
+
+   template <int DIM, int D1D, int Q1D>
+   static void AddSpecialization()
+   {
+      ApplyPAKernels::Specialization<DIM,D1D,Q1D>::Add();
+      DiagonalPAKernels::Specialization<DIM,D1D,Q1D>::Add();
+   }
 };
 
 /** Class for local mass matrix assembling $a(u,v) := (Q u, v)$ */
@@ -2306,6 +2331,20 @@ protected:
    const FaceGeometricFactors *face_geom; ///< Not owned
    int dim, ne, nq, dofs1D, quad1D;
 
+public:
+
+   using ApplyKernelType = void(*)(const int, const Array<real_t>&,
+                                   const Array<real_t>&, const Vector&,
+                                   const Vector&, Vector&, const int, const int);
+
+   using DiagonalKernelType =  void(*)(const int, const Array<real_t>&,
+                                       const Vector&, Vector&, const int,
+                                       const int);
+
+   MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType, (int, int, int));
+   MFEM_REGISTER_KERNELS(DiagonalPAKernels, DiagonalKernelType, (int, int, int));
+   static struct Kernels { Kernels(); } kernels;
+
 public:
    MassIntegrator(const IntegrationRule *ir = NULL)
       : BilinearFormIntegrator(ir), Q(NULL), maps(NULL), geom(NULL) { }
@@ -2351,6 +2390,13 @@ public:
    bool SupportsCeed() const { return DeviceCanUseCeed(); }
 
    const Coefficient *GetCoefficient() const { return Q; }
+
+   template <int DIM, int D1D, int Q1D>
+   static void AddSpecialization()
+   {
+      ApplyPAKernels::Specialization<DIM,D1D,Q1D>::Add();
+      DiagonalPAKernels::Specialization<DIM,D1D,Q1D>::Add();
+   }
 };
 
 /** Mass integrator $(u, v)$ restricted to the boundary of a domain */

fem/integ/bilininteg_diffusion_kernels.cpp

@@ -14,6 +14,33 @@
 namespace mfem
 {
 
+// PA Diffusion Integrator
+
+DiffusionIntegrator::Kernels DiffusionIntegrator::kernels;
+DiffusionIntegrator::Kernels::Kernels()
+{
+   // 2D
+   DiffusionIntegrator::AddSpecialization<2,2,2>();
+   DiffusionIntegrator::AddSpecialization<2,3,3>();
+   DiffusionIntegrator::AddSpecialization<2,4,4>();
+   DiffusionIntegrator::AddSpecialization<2,5,5>();
+   DiffusionIntegrator::AddSpecialization<2,6,6>();
+   DiffusionIntegrator::AddSpecialization<2,7,7>();
+   DiffusionIntegrator::AddSpecialization<2,8,8>();
+   DiffusionIntegrator::AddSpecialization<2,9,9>();
+   // 3D
+   DiffusionIntegrator::AddSpecialization<3,2,2>();
+   DiffusionIntegrator::AddSpecialization<3,2,3>();
+   DiffusionIntegrator::AddSpecialization<3,3,4>();
+   DiffusionIntegrator::AddSpecialization<3,4,5>();
+   DiffusionIntegrator::AddSpecialization<3,4,6>();
+   DiffusionIntegrator::AddSpecialization<3,5,6>();
+   DiffusionIntegrator::AddSpecialization<3,5,8>();
+   DiffusionIntegrator::AddSpecialization<3,6,7>();
+   DiffusionIntegrator::AddSpecialization<3,7,8>();
+   DiffusionIntegrator::AddSpecialization<3,8,9>();
+}
+
 namespace internal
 {
 
@@ -363,118 +390,6 @@ void OccaPADiffusionSetup3D(const int D1D,
 }
 #endif // MFEM_USE_OCCA
 
-void PADiffusionAssembleDiagonal(const int dim,
-                                 const int D1D,
-                                 const int Q1D,
-                                 const int NE,
-                                 const bool symm,
-                                 const Array<real_t> &B,
-                                 const Array<real_t> &G,
-                                 const Vector &D,
-                                 Vector &Y)
-{
-   if (dim == 2)
-   {
-      switch ((D1D << 4 ) | Q1D)
-      {
-         case 0x22: return SmemPADiffusionDiagonal2D<2,2,8>(NE,symm,B,G,D,Y);
-         case 0x33: return SmemPADiffusionDiagonal2D<3,3,8>(NE,symm,B,G,D,Y);
-         case 0x44: return SmemPADiffusionDiagonal2D<4,4,4>(NE,symm,B,G,D,Y);
-         case 0x55: return SmemPADiffusionDiagonal2D<5,5,4>(NE,symm,B,G,D,Y);
-         case 0x66: return SmemPADiffusionDiagonal2D<6,6,2>(NE,symm,B,G,D,Y);
-         case 0x77: return SmemPADiffusionDiagonal2D<7,7,2>(NE,symm,B,G,D,Y);
-         case 0x88: return SmemPADiffusionDiagonal2D<8,8,1>(NE,symm,B,G,D,Y);
-         case 0x99: return SmemPADiffusionDiagonal2D<9,9,1>(NE,symm,B,G,D,Y);
-         default: return PADiffusionDiagonal2D(NE,symm,B,G,D,Y,D1D,Q1D);
-      }
-   }
-   else if (dim == 3)
-   {
-      switch ((D1D << 4 ) | Q1D)
-      {
-         case 0x22: return SmemPADiffusionDiagonal3D<2,2>(NE,symm,B,G,D,Y);
-         case 0x23: return SmemPADiffusionDiagonal3D<2,3>(NE,symm,B,G,D,Y);
-         case 0x34: return SmemPADiffusionDiagonal3D<3,4>(NE,symm,B,G,D,Y);
-         case 0x45: return SmemPADiffusionDiagonal3D<4,5>(NE,symm,B,G,D,Y);
-         case 0x46: return SmemPADiffusionDiagonal3D<4,6>(NE,symm,B,G,D,Y);
-         case 0x56: return SmemPADiffusionDiagonal3D<5,6>(NE,symm,B,G,D,Y);
-         case 0x67: return SmemPADiffusionDiagonal3D<6,7>(NE,symm,B,G,D,Y);
-         case 0x78: return SmemPADiffusionDiagonal3D<7,8>(NE,symm,B,G,D,Y);
-         case 0x89: return SmemPADiffusionDiagonal3D<8,9>(NE,symm,B,G,D,Y);
-         case 0x9A: return SmemPADiffusionDiagonal3D<9,10>(NE,symm,B,G,D,Y);
-         default: return PADiffusionDiagonal3D(NE,symm,B,G,D,Y,D1D,Q1D);
-      }
-   }
-   MFEM_ABORT("Unknown kernel.");
-}
-
-void PADiffusionApply(const int dim,
-                      const int D1D,
-                      const int Q1D,
-                      const int NE,
-                      const bool symm,
-                      const Array<real_t> &B,
-                      const Array<real_t> &G,
-                      const Array<real_t> &Bt,
-                      const Array<real_t> &Gt,
-                      const Vector &D,
-                      const Vector &X,
-                      Vector &Y)
-{
-#ifdef MFEM_USE_OCCA
-   if (DeviceCanUseOcca())
-   {
-      if (dim == 2)
-      {
-         OccaPADiffusionApply2D(D1D,Q1D,NE,B,G,Bt,Gt,D,X,Y);
-         return;
-      }
-      if (dim == 3)
-      {
-         OccaPADiffusionApply3D(D1D,Q1D,NE,B,G,Bt,Gt,D,X,Y);
-         return;
-      }
-      MFEM_ABORT("OCCA PADiffusionApply unknown kernel!");
-   }
-#endif // MFEM_USE_OCCA
-   const int id = (D1D << 4) | Q1D;
-
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x22: return SmemPADiffusionApply2D<2,2,16>(NE,symm,B,G,D,X,Y);
-         case 0x33: return SmemPADiffusionApply2D<3,3,16>(NE,symm,B,G,D,X,Y);
-         case 0x44: return SmemPADiffusionApply2D<4,4,8>(NE,symm,B,G,D,X,Y);
-         case 0x55: return SmemPADiffusionApply2D<5,5,8>(NE,symm,B,G,D,X,Y);
-         case 0x66: return SmemPADiffusionApply2D<6,6,4>(NE,symm,B,G,D,X,Y);
-         case 0x77: return SmemPADiffusionApply2D<7,7,4>(NE,symm,B,G,D,X,Y);
-         case 0x88: return SmemPADiffusionApply2D<8,8,2>(NE,symm,B,G,D,X,Y);
-         case 0x99: return SmemPADiffusionApply2D<9,9,2>(NE,symm,B,G,D,X,Y);
-         default:   return PADiffusionApply2D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
-      }
-   }
-
-   if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x22: return SmemPADiffusionApply3D<2,2>(NE,symm,B,G,D,X,Y);
-         case 0x23: return SmemPADiffusionApply3D<2,3>(NE,symm,B,G,D,X,Y);
-         case 0x34: return SmemPADiffusionApply3D<3,4>(NE,symm,B,G,D,X,Y);
-         case 0x45: return SmemPADiffusionApply3D<4,5>(NE,symm,B,G,D,X,Y);
-         case 0x46: return SmemPADiffusionApply3D<4,6>(NE,symm,B,G,D,X,Y);
-         case 0x56: return SmemPADiffusionApply3D<5,6>(NE,symm,B,G,D,X,Y);
-         case 0x58: return SmemPADiffusionApply3D<5,8>(NE,symm,B,G,D,X,Y);
-         case 0x67: return SmemPADiffusionApply3D<6,7>(NE,symm,B,G,D,X,Y);
-         case 0x78: return SmemPADiffusionApply3D<7,8>(NE,symm,B,G,D,X,Y);
-         case 0x89: return SmemPADiffusionApply3D<8,9>(NE,symm,B,G,D,X,Y);
-         default:   return PADiffusionApply3D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
-      }
-   }
-   MFEM_ABORT("Unknown kernel: 0x"<<std::hex << id << std::dec);
-}
-
 #ifdef MFEM_USE_OCCA
 void OccaPADiffusionApply2D(const int D1D,
                             const int Q1D,

fem/integ/bilininteg_diffusion_kernels.hpp

@@ -12,6 +12,7 @@
 #ifndef MFEM_BILININTEG_DIFFUSION_KERNELS_HPP
 #define MFEM_BILININTEG_DIFFUSION_KERNELS_HPP
 
+#include "../kernel_dispatch.hpp"
 #include "../../config/config.hpp"
 #include "../../general/array.hpp"
 #include "../../general/forall.hpp"
@@ -36,7 +37,7 @@ void PADiffusionSetup(const int dim,
                       const Vector &C,
                       Vector &D);
 
-// PA Diffusion Assemble 2D kernel
+// PA Diffusion Assemble 2D f
 template<int T_SDIM>
 void PADiffusionSetup2D(const int Q1D,
                         const int coeffDim,
@@ -151,8 +152,23 @@ inline void PADiffusionDiagonal2D(const int NE,
    });
 }
 
+namespace diffusion
+{
+constexpr int ipow(int x, int p) { return p == 0 ? 1 : x*ipow(x, p-1); }
+constexpr int D11(int x) { return (11 - x)/2; }
+constexpr int D10(int x) { return (10 - x)/2; }
+constexpr int NBZApply(int D1D)
+{
+   return ipow(2, D11(D1D) >= 0 ? D11(D1D) : 0);
+}
+constexpr int NBZDiagonal(int D1D)
+{
+   return ipow(2, D10(D1D) >= 0 ? D10(D1D) : 0);
+}
+}
+
 // Shared memory PA Diffusion Diagonal 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
+template<int T_D1D = 0, int T_Q1D = 0>
 inline void SmemPADiffusionDiagonal2D(const int NE,
                                       const bool symmetric,
                                       const Array<real_t> &b_,
@@ -162,9 +178,10 @@ inline void SmemPADiffusionDiagonal2D(const int NE,
                                       const int d1d = 0,
                                       const int q1d = 0)
 {
+   static constexpr int T_NBZ = diffusion::NBZDiagonal(T_D1D);
+   static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
-   constexpr int NBZ = T_NBZ ? T_NBZ : 1;
    const int max_q1d = T_Q1D ? T_Q1D : DeviceDofQuadLimits::Get().MAX_Q1D;
    const int max_d1d = T_D1D ? T_D1D : DeviceDofQuadLimits::Get().MAX_D1D;
    MFEM_VERIFY(D1D <= max_d1d, "");
@@ -178,7 +195,6 @@ inline void SmemPADiffusionDiagonal2D(const int NE,
       const int tidz = MFEM_THREAD_ID(z);
       const int D1D = T_D1D ? T_D1D : d1d;
       const int Q1D = T_Q1D ? T_Q1D : q1d;
-      constexpr int NBZ = T_NBZ ? T_NBZ : 1;
       constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
       constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
       MFEM_SHARED real_t BG[2][MQ1*MD1];
@@ -628,20 +644,23 @@ inline void PADiffusionApply2D(const int NE,
 }
 
 // Shared memory PA Diffusion Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
+template<int T_D1D = 0, int T_Q1D = 0>
 inline void SmemPADiffusionApply2D(const int NE,
                                    const bool symmetric,
                                    const Array<real_t> &b_,
                                    const Array<real_t> &g_,
+                                   const Array<real_t> &bt_,
+                                   const Array<real_t> &gt_,
                                    const Vector &d_,
                                    const Vector &x_,
                                    Vector &y_,
                                    const int d1d = 0,
                                    const int q1d = 0)
 {
+   static constexpr int T_NBZ = diffusion::NBZApply(T_D1D);
+   static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
-   constexpr int NBZ = T_NBZ ? T_NBZ : 1;
    const int max_q1d = T_Q1D ? T_Q1D : DeviceDofQuadLimits::Get().MAX_Q1D;
    const int max_d1d = T_D1D ? T_D1D : DeviceDofQuadLimits::Get().MAX_D1D;
    MFEM_VERIFY(D1D <= max_d1d, "");
@@ -656,7 +675,6 @@ inline void SmemPADiffusionApply2D(const int NE,
       const int tidz = MFEM_THREAD_ID(z);
       const int D1D = T_D1D ? T_D1D : d1d;
       const int Q1D = T_Q1D ? T_Q1D : q1d;
-      constexpr int NBZ = T_NBZ ? T_NBZ : 1;
       constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
       constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
       MFEM_SHARED real_t sBG[2][MQ1*MD1];
@@ -984,6 +1002,8 @@ inline void SmemPADiffusionApply3D(const int NE,
                                    const bool symmetric,
                                    const Array<real_t> &b_,
                                    const Array<real_t> &g_,
+                                   const Array<real_t> &,
+                                   const Array<real_t> &,
                                    const Vector &d_,
                                    const Vector &x_,
                                    Vector &y_,
@@ -1203,6 +1223,44 @@ inline void SmemPADiffusionApply3D(const int NE,
 
 } // namespace internal
 
+namespace
+{
+using ApplyKernelType = DiffusionIntegrator::ApplyKernelType;
+using DiagonalKernelType = DiffusionIntegrator::DiagonalKernelType;
+}
+
+template<int DIM, int T_D1D, int T_Q1D>
+ApplyKernelType DiffusionIntegrator::ApplyPAKernels::Kernel()
+{
+   if (DIM == 2) { return internal::SmemPADiffusionApply2D<T_D1D,T_Q1D>; }
+   else if (DIM == 3) { return internal::SmemPADiffusionApply3D<T_D1D, T_Q1D>; }
+   else { MFEM_ABORT(""); }
+}
+
+inline
+ApplyKernelType DiffusionIntegrator::ApplyPAKernels::Fallback(int DIM, int, int)
+{
+   if (DIM == 2) { return internal::PADiffusionApply2D; }
+   else if (DIM == 3) { return internal::PADiffusionApply3D; }
+   else { MFEM_ABORT(""); }
+}
+
+template<int DIM, int D1D, int Q1D>
+DiagonalKernelType DiffusionIntegrator::DiagonalPAKernels::Kernel()
+{
+   if (DIM == 2) { return internal::SmemPADiffusionDiagonal2D<D1D,Q1D>; }
+   else if (DIM == 3) { return internal::SmemPADiffusionDiagonal3D<D1D, Q1D>; }
+   else { MFEM_ABORT(""); }
+}
+
+inline DiagonalKernelType
+DiffusionIntegrator::DiagonalPAKernels::Fallback(int DIM, int, int)
+{
+   if (DIM == 2) { return internal::PADiffusionDiagonal2D; }
+   else if (DIM == 3) { return internal::PADiffusionDiagonal3D; }
+   else { MFEM_ABORT(""); }
+}
+
 } // namespace mfem
 
 #endif

fem/integ/bilininteg_diffusion_pa.cpp

@@ -19,6 +19,73 @@
 namespace mfem
 {
 
+void DiffusionIntegrator::AssembleDiagonalPA(Vector &diag)
+{
+   if (DeviceCanUseCeed())
+   {
+      ceedOp->GetDiagonal(diag);
+   }
+   else
+   {
+      if (pa_data.Size() == 0) { AssemblePA(*fespace); }
+      const Array<real_t> &B = maps->B;
+      const Array<real_t> &G = maps->G;
+      const Vector &Dv = pa_data;
+      DiagonalPAKernels::Run(dim, dofs1D, quad1D, ne, symmetric, B, G, Dv,
+                             diag, dofs1D, quad1D);
+   }
+}
+
+// PA Diffusion Apply kernel
+void DiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
+{
+   if (DeviceCanUseCeed())
+   {
+      ceedOp->AddMult(x, y);
+   }
+   else
+   {
+      const Array<real_t> &B = maps->B;
+      const Array<real_t> &G = maps->G;
+      const Array<real_t> &Bt = maps->Bt;
+      const Array<real_t> &Gt = maps->Gt;
+      const Vector &Dv = pa_data;
+
+#ifdef MFEM_USE_OCCA
+      if (DeviceCanUseOcca())
+      {
+         if (dim == 2)
+         {
+            OccaPADiffusionApply2D(dofs1D,quad1D,ne,B,G,Bt,Gt,Dv,x,y);
+            return;
+         }
+         if (dim == 3)
+         {
+            OccaPADiffusionApply3D(dofs1D,quad1D,ne,B,G,Bt,Gt,Dv,x,y);
+            return;
+         }
+         MFEM_ABORT("OCCA PADiffusionApply unknown kernel!");
+      }
+#endif // MFEM_USE_OCCA
+
+      ApplyPAKernels::Run(dim, dofs1D, quad1D, ne, symmetric, B, G, Bt,
+                          Gt, Dv, x, y, dofs1D, quad1D);
+   }
+}
+
+void DiffusionIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
+{
+   if (symmetric)
+   {
+      AddMultPA(x, y);
+   }
+   else
+   {
+      MFEM_ABORT("DiffusionIntegrator::AddMultTransposePA only implemented in "
+                 "the symmetric case.")
+   }
+}
+
 void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
@@ -98,47 +165,6 @@ void DiffusionIntegrator::AssemblePatchPA(const int patch,
    SetupPatchPA(patch, mesh);  // For full quadrature, unitWeights = false
 }
 
-void DiffusionIntegrator::AssembleDiagonalPA(Vector &diag)
-{
-   if (DeviceCanUseCeed())
-   {
-      ceedOp->GetDiagonal(diag);
-   }
-   else
-   {
-      if (pa_data.Size()==0) { AssemblePA(*fespace); }
-      internal::PADiffusionAssembleDiagonal(dim, dofs1D, quad1D, ne, symmetric,
-                                            maps->B, maps->G, pa_data, diag);
-   }
-}
-
-void DiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
-{
-   if (DeviceCanUseCeed())
-   {
-      ceedOp->AddMult(x, y);
-   }
-   else
-   {
-      internal::PADiffusionApply(dim, dofs1D, quad1D, ne, symmetric,
-                                 maps->B, maps->G, maps->Bt, maps->Gt,
-                                 pa_data, x, y);
-   }
-}
-
-void DiffusionIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
-{
-   if (symmetric)
-   {
-      AddMultPA(x, y);
-   }
-   else
-   {
-      MFEM_ABORT("DiffusionIntegrator::AddMultTransposePA only implemented in "
-                 "the symmetric case.")
-   }
-}
-
 // This version uses full 1D quadrature rules, taking into account the
 // minimum interaction between basis functions and integration points.
 void DiffusionIntegrator::AddMultPatchPA(const int patch, const Vector &x,

fem/integ/bilininteg_mass_kernels.cpp

@@ -14,78 +14,34 @@
 namespace mfem
 {
 
+MassIntegrator::Kernels MassIntegrator::kernels;
+MassIntegrator::Kernels::Kernels()
+{
+   // 2D
+   MassIntegrator::AddSpecialization<2,2,2>();
+   MassIntegrator::AddSpecialization<2,3,3>();
+   MassIntegrator::AddSpecialization<2,4,4>();
+   MassIntegrator::AddSpecialization<2,5,5>();
+   MassIntegrator::AddSpecialization<2,6,6>();
+   MassIntegrator::AddSpecialization<2,7,7>();
+   MassIntegrator::AddSpecialization<2,8,8>();
+   MassIntegrator::AddSpecialization<2,9,9>();
+   // 3D
+   MassIntegrator::AddSpecialization<3,2,2>();
+   MassIntegrator::AddSpecialization<3,2,3>();
+   MassIntegrator::AddSpecialization<3,3,4>();
+   MassIntegrator::AddSpecialization<3,4,5>();
+   MassIntegrator::AddSpecialization<3,4,6>();
+   MassIntegrator::AddSpecialization<3,5,6>();
+   MassIntegrator::AddSpecialization<3,5,8>();
+   MassIntegrator::AddSpecialization<3,6,7>();
+   MassIntegrator::AddSpecialization<3,7,8>();
+   MassIntegrator::AddSpecialization<3,8,9>();
+}
+
 namespace internal
 {
 
-// PA Mass Diagonal 1D kernel
-static void PAMassAssembleDiagonal1D(const int NE,
-                                     const Array<real_t> &b,
-                                     const Vector &d,
-                                     Vector &y,
-                                     const int D1D,
-                                     const int Q1D)
-{
-   auto B = Reshape(b.Read(), Q1D, D1D);
-   auto D = Reshape(d.Read(), Q1D, NE);
-   auto Y = Reshape(y.ReadWrite(), D1D, NE);
-   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
-   {
-      for (int dx = 0; dx < D1D; ++dx)
-      {
-         for (int qx = 0; qx < Q1D; ++qx)
-         {
-            Y(dx, e) += B(qx, dx) * B(qx, dx) * D(qx, e);
-         }
-      }
-   });
-}
-
-void PAMassAssembleDiagonal(const int dim, const int D1D,
-                            const int Q1D, const int NE,
-                            const Array<real_t> &B,
-                            const Vector &D,
-                            Vector &Y)
-{
-   if (dim == 1)
-   {
-      return PAMassAssembleDiagonal1D(NE,B,D,Y,D1D,Q1D);
-   }
-   else if (dim == 2)
-   {
-      switch ((D1D << 4 ) | Q1D)
-      {
-         case 0x22: return SmemPAMassAssembleDiagonal2D<2,2,16>(NE,B,D,Y);
-         case 0x33: return SmemPAMassAssembleDiagonal2D<3,3,16>(NE,B,D,Y);
-         case 0x44: return SmemPAMassAssembleDiagonal2D<4,4,8>(NE,B,D,Y);
-         case 0x55: return SmemPAMassAssembleDiagonal2D<5,5,8>(NE,B,D,Y);
-         case 0x66: return SmemPAMassAssembleDiagonal2D<6,6,4>(NE,B,D,Y);
-         case 0x77: return SmemPAMassAssembleDiagonal2D<7,7,4>(NE,B,D,Y);
-         case 0x88: return SmemPAMassAssembleDiagonal2D<8,8,2>(NE,B,D,Y);
-         case 0x99: return SmemPAMassAssembleDiagonal2D<9,9,2>(NE,B,D,Y);
-         default:   return PAMassAssembleDiagonal2D(NE,B,D,Y,D1D,Q1D);
-      }
-   }
-   else if (dim == 3)
-   {
-      switch ((D1D << 4 ) | Q1D)
-      {
-         case 0x23: return SmemPAMassAssembleDiagonal3D<2,3>(NE,B,D,Y);
-         case 0x24: return SmemPAMassAssembleDiagonal3D<2,4>(NE,B,D,Y);
-         case 0x26: return SmemPAMassAssembleDiagonal3D<2,6>(NE,B,D,Y);
-         case 0x34: return SmemPAMassAssembleDiagonal3D<3,4>(NE,B,D,Y);
-         case 0x35: return SmemPAMassAssembleDiagonal3D<3,5>(NE,B,D,Y);
-         case 0x45: return SmemPAMassAssembleDiagonal3D<4,5>(NE,B,D,Y);
-         case 0x48: return SmemPAMassAssembleDiagonal3D<4,8>(NE,B,D,Y);
-         case 0x56: return SmemPAMassAssembleDiagonal3D<5,6>(NE,B,D,Y);
-         case 0x67: return SmemPAMassAssembleDiagonal3D<6,7>(NE,B,D,Y);
-         case 0x78: return SmemPAMassAssembleDiagonal3D<7,8>(NE,B,D,Y);
-         case 0x89: return SmemPAMassAssembleDiagonal3D<8,9>(NE,B,D,Y);
-         default:   return PAMassAssembleDiagonal3D(NE,B,D,Y,D1D,Q1D);
-      }
-   }
-   MFEM_ABORT("Unknown kernel.");
-}
-
 #ifdef MFEM_USE_OCCA
 void OccaPAMassApply2D(const int D1D,
                        const int Q1D,
@@ -176,154 +132,6 @@ void OccaPAMassApply3D(const int D1D,
 }
 #endif // MFEM_USE_OCCA
 
-MFEM_HOST_DEVICE inline
-void PAMassApply1D_Element(const int e,
-                           const int NE,
-                           const real_t *b_,
-                           const real_t *bt_,
-                           const real_t *d_,
-                           const real_t *x_,
-                           real_t *y_,
-                           const int d1d = 0,
-                           const int q1d = 0)
-{
-   const int D1D = d1d;
-   const int Q1D = q1d;
-   auto B = ConstDeviceMatrix(b_, Q1D, D1D);
-   auto Bt = ConstDeviceMatrix(bt_, D1D, Q1D);
-   auto D = ConstDeviceMatrix(d_, Q1D, NE);
-   auto X = ConstDeviceMatrix(x_, D1D, NE);
-   auto Y = DeviceMatrix(y_, D1D, NE);
-
-   real_t XQ[DofQuadLimits::MAX_Q1D];
-   for (int qx = 0; qx < Q1D; ++qx)
-   {
-      XQ[qx] = 0.0;
-   }
-   for (int dx = 0; dx < D1D; ++dx)
-   {
-      const real_t s = X(dx,e);
-      for (int qx = 0; qx < Q1D; ++qx)
-      {
-         XQ[qx] += B(qx,dx)*s;
-      }
-   }
-   for (int qx = 0; qx < Q1D; ++qx)
-   {
-      const real_t q = XQ[qx]*D(qx,e);
-      for (int dx = 0; dx < D1D; ++dx)
-      {
-         Y(dx,e) += Bt(dx,qx) * q;
-      }
-   }
-}
-
-// PA Mass Apply 1D kernel
-static void PAMassApply1D(const int NE,
-                          const Array<real_t> &b_,
-                          const Array<real_t> &bt_,
-                          const Vector &d_,
-                          const Vector &x_,
-                          Vector &y_,
-                          const int d1d = 0,
-                          const int q1d = 0)
-{
-   MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
-   MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
-
-   const auto B = b_.Read();
-   const auto Bt = bt_.Read();
-   const auto D = d_.Read();
-   const auto X = x_.Read();
-   auto Y = y_.ReadWrite();
-
-   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
-   {
-      internal::PAMassApply1D_Element(e, NE, B, Bt, D, X, Y, d1d, q1d);
-   });
-}
-
-void PAMassApply(const int dim,
-                 const int D1D,
-                 const int Q1D,
-                 const int NE,
-                 const Array<real_t> &B,
-                 const Array<real_t> &Bt,
-                 const Vector &D,
-                 const Vector &X,
-                 Vector &Y)
-{
-#ifdef MFEM_USE_OCCA
-   if (DeviceCanUseOcca())
-   {
-      if (dim == 2)
-      {
-         return OccaPAMassApply2D(D1D,Q1D,NE,B,Bt,D,X,Y);
-      }
-      if (dim == 3)
-      {
-         return OccaPAMassApply3D(D1D,Q1D,NE,B,Bt,D,X,Y);
-      }
-      MFEM_ABORT("OCCA PA Mass Apply unknown kernel!");
-   }
-#endif // MFEM_USE_OCCA
-   const int id = (D1D << 4) | Q1D;
-
-   if (dim == 1)
-   {
-      return PAMassApply1D(NE,B,Bt,D,X,Y,D1D,Q1D);
-   }
-   else if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x22: return SmemPAMassApply2D<2,2,16>(NE,B,Bt,D,X,Y);
-         case 0x24: return SmemPAMassApply2D<2,4,16>(NE,B,Bt,D,X,Y);
-         case 0x33: return SmemPAMassApply2D<3,3,16>(NE,B,Bt,D,X,Y);
-         case 0x34: return SmemPAMassApply2D<3,4,16>(NE,B,Bt,D,X,Y);
-         case 0x35: return SmemPAMassApply2D<3,5,16>(NE,B,Bt,D,X,Y);
-         case 0x36: return SmemPAMassApply2D<3,6,16>(NE,B,Bt,D,X,Y);
-         case 0x44: return SmemPAMassApply2D<4,4,8>(NE,B,Bt,D,X,Y);
-         case 0x46: return SmemPAMassApply2D<4,6,8>(NE,B,Bt,D,X,Y);
-         case 0x48: return SmemPAMassApply2D<4,8,4>(NE,B,Bt,D,X,Y);
-         case 0x55: return SmemPAMassApply2D<5,5,8>(NE,B,Bt,D,X,Y);
-         case 0x57: return SmemPAMassApply2D<5,7,8>(NE,B,Bt,D,X,Y);
-         case 0x58: return SmemPAMassApply2D<5,8,2>(NE,B,Bt,D,X,Y);
-         case 0x66: return SmemPAMassApply2D<6,6,4>(NE,B,Bt,D,X,Y);
-         case 0x77: return SmemPAMassApply2D<7,7,4>(NE,B,Bt,D,X,Y);
-         case 0x88: return SmemPAMassApply2D<8,8,2>(NE,B,Bt,D,X,Y);
-         case 0x99: return SmemPAMassApply2D<9,9,2>(NE,B,Bt,D,X,Y);
-         default:   return PAMassApply2D(NE,B,Bt,D,X,Y,D1D,Q1D);
-      }
-   }
-   else if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x22: return SmemPAMassApply3D<2,2>(NE,B,Bt,D,X,Y);
-         case 0x23: return SmemPAMassApply3D<2,3>(NE,B,Bt,D,X,Y);
-         case 0x24: return SmemPAMassApply3D<2,4>(NE,B,Bt,D,X,Y);
-         case 0x26: return SmemPAMassApply3D<2,6>(NE,B,Bt,D,X,Y);
-         case 0x34: return SmemPAMassApply3D<3,4>(NE,B,Bt,D,X,Y);
-         case 0x35: return SmemPAMassApply3D<3,5>(NE,B,Bt,D,X,Y);
-         case 0x36: return SmemPAMassApply3D<3,6>(NE,B,Bt,D,X,Y);
-         case 0x37: return SmemPAMassApply3D<3,7>(NE,B,Bt,D,X,Y);
-         case 0x45: return SmemPAMassApply3D<4,5>(NE,B,Bt,D,X,Y);
-         case 0x46: return SmemPAMassApply3D<4,6>(NE,B,Bt,D,X,Y);
-         case 0x48: return SmemPAMassApply3D<4,8>(NE,B,Bt,D,X,Y);
-         case 0x56: return SmemPAMassApply3D<5,6>(NE,B,Bt,D,X,Y);
-         case 0x58: return SmemPAMassApply3D<5,8>(NE,B,Bt,D,X,Y);
-         case 0x67: return SmemPAMassApply3D<6,7>(NE,B,Bt,D,X,Y);
-         case 0x78: return SmemPAMassApply3D<7,8>(NE,B,Bt,D,X,Y);
-         case 0x89: return SmemPAMassApply3D<8,9>(NE,B,Bt,D,X,Y);
-         case 0x9A: return SmemPAMassApply3D<9,10>(NE,B,Bt,D,X,Y);
-         default:   return PAMassApply3D(NE,B,Bt,D,X,Y,D1D,Q1D);
-      }
-   }
-   mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
-   MFEM_ABORT("Unknown kernel.");
-}
-
 } // namespace internal
 
 } // namespace mfem

fem/integ/bilininteg_mass_kernels.hpp

@@ -25,11 +25,95 @@ namespace mfem
 namespace internal
 {
 
-void PAMassAssembleDiagonal(const int dim, const int D1D,
-                            const int Q1D, const int NE,
-                            const Array<real_t> &B,
-                            const Vector &D,
-                            Vector &Y);
+// PA Mass Diagonal 1D kernel
+static void PAMassAssembleDiagonal1D(const int NE,
+                                     const Array<real_t> &b,
+                                     const Vector &d,
+                                     Vector &y,
+                                     const int D1D,
+                                     const int Q1D)
+{
+   auto B = Reshape(b.Read(), Q1D, D1D);
+   auto D = Reshape(d.Read(), Q1D, NE);
+   auto Y = Reshape(y.ReadWrite(), D1D, NE);
+   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   {
+      for (int dx = 0; dx < D1D; ++dx)
+      {
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            Y(dx, e) += B(qx, dx) * B(qx, dx) * D(qx, e);
+         }
+      }
+   });
+}
+
+MFEM_HOST_DEVICE inline
+void PAMassApply1D_Element(const int e,
+                           const int NE,
+                           const real_t *b_,
+                           const real_t *bt_,
+                           const real_t *d_,
+                           const real_t *x_,
+                           real_t *y_,
+                           const int d1d = 0,
+                           const int q1d = 0)
+{
+   const int D1D = d1d;
+   const int Q1D = q1d;
+   auto B = ConstDeviceMatrix(b_, Q1D, D1D);
+   auto Bt = ConstDeviceMatrix(bt_, D1D, Q1D);
+   auto D = ConstDeviceMatrix(d_, Q1D, NE);
+   auto X = ConstDeviceMatrix(x_, D1D, NE);
+   auto Y = DeviceMatrix(y_, D1D, NE);
+
+   real_t XQ[DofQuadLimits::MAX_Q1D];
+   for (int qx = 0; qx < Q1D; ++qx)
+   {
+      XQ[qx] = 0.0;
+   }
+   for (int dx = 0; dx < D1D; ++dx)
+   {
+      const real_t s = X(dx,e);
+      for (int qx = 0; qx < Q1D; ++qx)
+      {
+         XQ[qx] += B(qx,dx)*s;
+      }
+   }
+   for (int qx = 0; qx < Q1D; ++qx)
+   {
+      const double q = XQ[qx]*D(qx,e);
+      for (int dx = 0; dx < D1D; ++dx)
+      {
+         Y(dx,e) += Bt(dx,qx) * q;
+      }
+   }
+}
+
+// PA Mass Apply 1D kernel
+static void PAMassApply1D(const int NE,
+                          const Array<real_t> &b_,
+                          const Array<real_t> &bt_,
+                          const Vector &d_,
+                          const Vector &x_,
+                          Vector &y_,
+                          const int d1d = 0,
+                          const int q1d = 0)
+{
+   MFEM_VERIFY(d1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
+
+   const auto B = b_.Read();
+   const auto Bt = bt_.Read();
+   const auto D = d_.Read();
+   const auto X = x_.Read();
+   auto Y = y_.ReadWrite();
+
+   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   {
+      internal::PAMassApply1D_Element(e, NE, B, Bt, D, X, Y, d1d, q1d);
+   });
+}
 
 // PA Mass Diagonal 2D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
@@ -78,8 +162,18 @@ inline void PAMassAssembleDiagonal2D(const int NE,
    });
 }
 
+namespace mass
+{
+constexpr int ipow(int x, int p) { return p == 0 ? 1 : x*ipow(x, p-1); }
+constexpr int D(int D1D) { return (11 - D1D) / 2; }
+constexpr int NBZ(int D1D)
+{
+   return ipow(2, D(D1D) >= 0 ? D(D1D) : 0);
+}
+}
+
 // Shared memory PA Mass Diagonal 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
+template<int T_D1D = 0, int T_Q1D = 0>
 inline void SmemPAMassAssembleDiagonal2D(const int NE,
                                          const Array<real_t> &b_,
                                          const Vector &d_,
@@ -87,9 +181,10 @@ inline void SmemPAMassAssembleDiagonal2D(const int NE,
                                          const int d1d = 0,
                                          const int q1d = 0)
 {
+   static constexpr int T_NBZ = mass::NBZ(T_D1D);
+   static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
-   constexpr int NBZ = T_NBZ ? T_NBZ : 1;
    const int max_q1d = T_Q1D ? T_Q1D : DeviceDofQuadLimits::Get().MAX_Q1D;
    const int max_d1d = T_D1D ? T_D1D : DeviceDofQuadLimits::Get().MAX_D1D;
    MFEM_VERIFY(D1D <= max_d1d, "");
@@ -102,7 +197,6 @@ inline void SmemPAMassAssembleDiagonal2D(const int NE,
       const int tidz = MFEM_THREAD_ID(z);
       const int D1D = T_D1D ? T_D1D : d1d;
       const int Q1D = T_Q1D ? T_Q1D : q1d;
-      constexpr int NBZ = T_NBZ ? T_NBZ : 1;
       constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
       constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
       MFEM_SHARED real_t B[MQ1][MD1];
@@ -302,16 +396,6 @@ inline void SmemPAMassAssembleDiagonal3D(const int NE,
    });
 }
 
-void PAMassApply(const int dim,
-                 const int D1D,
-                 const int Q1D,
-                 const int NE,
-                 const Array<real_t> &B,
-                 const Array<real_t> &Bt,
-                 const Vector &D,
-                 const Vector &X,
-                 Vector &Y);
-
 #ifdef MFEM_USE_OCCA
 // OCCA PA Mass Apply 2D kernel
 void OccaPAMassApply2D(const int D1D,
@@ -964,7 +1048,7 @@ inline void PAMassApply2D(const int NE,
 }
 
 // Shared memory PA Mass Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
+template<int T_D1D = 0, int T_Q1D = 0>
 inline void SmemPAMassApply2D(const int NE,
                               const Array<real_t> &b_,
                               const Array<real_t> &bt_,
@@ -975,9 +1059,10 @@ inline void SmemPAMassApply2D(const int NE,
                               const int q1d = 0)
 {
    MFEM_CONTRACT_VAR(bt_);
+   static constexpr int T_NBZ = mass::NBZ(T_D1D);
+   static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
-   constexpr int NBZ = T_NBZ ? T_NBZ : 1;
    const int max_q1d = T_Q1D ? T_Q1D : DeviceDofQuadLimits::Get().MAX_Q1D;
    const int max_d1d = T_D1D ? T_D1D : DeviceDofQuadLimits::Get().MAX_D1D;
    MFEM_VERIFY(D1D <= max_d1d, "");
@@ -988,8 +1073,8 @@ inline void SmemPAMassApply2D(const int NE,
    auto Y = y_.ReadWrite();
    mfem::forall_2D_batch(NE, Q1D, Q1D, NBZ, [=] MFEM_HOST_DEVICE (int e)
    {
-      internal::SmemPAMassApply2D_Element<T_D1D,T_Q1D,T_NBZ>(e, NE, b, D, x, Y, d1d,
-                                                             q1d);
+      internal::SmemPAMassApply2D_Element<T_D1D,T_Q1D,T_NBZ>(
+         e, NE, b, D, x, Y, d1d, q1d);
    });
 }
 
@@ -1049,6 +1134,48 @@ inline void SmemPAMassApply3D(const int NE,
 
 } // namespace internal
 
+namespace
+{
+using ApplyKernelType = MassIntegrator::ApplyKernelType;
+using DiagonalKernelType = MassIntegrator::DiagonalKernelType;
+}
+
+template<int DIM, int T_D1D, int T_Q1D>
+ApplyKernelType MassIntegrator::ApplyPAKernels::Kernel()
+{
+   if (DIM == 1) { return internal::PAMassApply1D; }
+   else if (DIM == 2) { return internal::SmemPAMassApply2D<T_D1D,T_Q1D>; }
+   else if (DIM == 3) { return internal::SmemPAMassApply3D<T_D1D, T_Q1D>; }
+   else { MFEM_ABORT(""); }
+}
+
+inline ApplyKernelType MassIntegrator::ApplyPAKernels::Fallback(
+   int DIM, int, int)
+{
+   if (DIM == 1) { return internal::PAMassApply1D; }
+   else if (DIM == 2) { return internal::PAMassApply2D; }
+   else if (DIM == 3) { return internal::PAMassApply3D; }
+   else { MFEM_ABORT(""); }
+}
+
+template<int DIM, int T_D1D, int T_Q1D>
+DiagonalKernelType MassIntegrator::DiagonalPAKernels::Kernel()
+{
+   if (DIM == 1) { return internal::PAMassAssembleDiagonal1D; }
+   else if (DIM == 2) { return internal::SmemPAMassAssembleDiagonal2D<T_D1D,T_Q1D>; }
+   else if (DIM == 3) { return internal::SmemPAMassAssembleDiagonal3D<T_D1D, T_Q1D>; }
+   else { MFEM_ABORT(""); }
+}
+
+inline DiagonalKernelType MassIntegrator::DiagonalPAKernels::Fallback(
+   int DIM, int, int)
+{
+   if (DIM == 1) { return internal::PAMassAssembleDiagonal1D; }
+   else if (DIM == 2) { return internal::PAMassAssembleDiagonal2D; }
+   else if (DIM == 3) { return internal::PAMassAssembleDiagonal3D; }
+   else { MFEM_ABORT(""); }
+}
+
 } // namespace mfem
 
 #endif

fem/integ/bilininteg_mass_pa.cpp

@@ -19,6 +19,8 @@
 namespace mfem
 {
 
+// PA Mass Integrator
+
 void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
@@ -195,8 +197,8 @@ void MassIntegrator::AssembleDiagonalPA(Vector &diag)
    }
    else
    {
-      internal::PAMassAssembleDiagonal(dim, dofs1D, quad1D, ne, maps->B, pa_data,
-                                       diag);
+      DiagonalPAKernels::Run(dim, dofs1D, quad1D, ne, maps->B, pa_data,
+                             diag, dofs1D, quad1D);
    }
 }
 
@@ -208,8 +210,26 @@ void MassIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
    else
    {
-      internal::PAMassApply(dim, dofs1D, quad1D, ne, maps->B, maps->Bt, pa_data, x,
-                            y);
+      const int D1D = dofs1D;
+      const int Q1D = quad1D;
+      const Array<real_t> &B = maps->B;
+      const Array<real_t> &Bt = maps->Bt;
+      const Vector &D = pa_data;
+#ifdef MFEM_USE_OCCA
+      if (DeviceCanUseOcca())
+      {
+         if (dim == 2)
+         {
+            return OccaPAMassApply2D(D1D,Q1D,ne,B,Bt,D,x,y);
+         }
+         if (dim == 3)
+         {
+            return OccaPAMassApply3D(D1D,Q1D,ne,B,Bt,D,x,y);
+         }
+         MFEM_ABORT("OCCA PA Mass Apply unknown kernel!");
+      }
+#endif // MFEM_USE_OCCA
+      ApplyPAKernels::Run(dim, D1D, Q1D, ne, B, Bt, D, x, y, D1D, Q1D);
    }
 }
 

fem/kernel_dispatch.hpp

@@ -0,0 +1,183 @@
+// Copyright (c) 2010-2024, 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_KERNEL_DISPATCH_HPP
+#define MFEM_KERNEL_DISPATCH_HPP
+
+#include "../config/config.hpp"
+#include "kernel_reporter.hpp"
+#include <unordered_map>
+#include <tuple>
+#include <cstddef>
+
+namespace mfem
+{
+
+// The MFEM_REGISTER_KERNELS macro registers kernels for runtime dispatch using
+// a dispatch map.
+//
+// This creates a dispatch table (a static member variable) named @a KernelName
+// containing function points of type @a KernelType. These are followed by one
+// or two sets of parenthesized argument types.
+//
+// The first set of argument types contains the types that are used to dispatch
+// to either specialized or fallback kernels. The second set of argument types
+// can be used to further specialize the kernel without participating in
+// dispatch (a canonical example is NBZ, determining the size of the thread
+// blocks; this is required to specialize kernels for optimal performance, but
+// is not relevant for dispatch).
+//
+// After calling this macro, the user must implement the Kernel and Fallback
+// static member functions, which return pointers to the appropriate kernel
+// functions depending on the parameters.
+//
+// Specialized functions can be registered using the static AddSpecialization
+// member function.
+
+#define MFEM_EXPAND(X) X // Workaround needed for MSVC compiler
+
+#define MFEM_REGISTER_KERNELS(KernelName, KernelType, ...)                     \
+   MFEM_EXPAND(MFEM_EXPAND(MFEM_REGISTER_KERNELS_N(__VA_ARGS__,2,1,))          \
+      (KernelName,KernelType,__VA_ARGS__))
+
+#define MFEM_REGISTER_KERNELS_N(_1, _2, N, ...) MFEM_REGISTER_KERNELS_##N
+
+// Expands a variable length macro parameter so that multiple variable length
+// parameters can be passed to the same macro.
+#define MFEM_PARAM_LIST(...) __VA_ARGS__
+
+// Version of MFEM_REGISTER_KERNELS without any "optional" (non-dispatch)
+// parameters.
+#define MFEM_REGISTER_KERNELS_1(KernelName, KernelType, Params)                \
+   MFEM_REGISTER_KERNELS_(KernelName, KernelType, Params, (), Params)
+
+// Version of MFEM_REGISTER_KERNELS without any optional (non-dispatch)
+// parameters (e.g. NBZ).
+#define MFEM_REGISTER_KERNELS_2(KernelName, KernelType, Params, OptParams)     \
+   MFEM_REGISTER_KERNELS_(KernelName, KernelType, Params, OptParams,           \
+                          (MFEM_PARAM_LIST Params, MFEM_PARAM_LIST OptParams))
+
+// P1 are the parameters, P2 are the optional (non-dispatch parameters), and P3
+// is the concatenation of P1 and P2. We need to pass it as a separate argument
+// to avoid a trailing comma in the case that P2 is empty.
+#define MFEM_REGISTER_KERNELS_(KernelName, KernelType, P1, P2, P3)             \
+   class KernelName : public                                                   \
+   KernelDispatchTable<KernelName, KernelType,                                 \
+      internal::KernelTypeList<MFEM_PARAM_LIST P1>,                            \
+      internal::KernelTypeList<MFEM_PARAM_LIST P2>>                            \
+   {                                                                           \
+   public:                                                                     \
+      const char *kernel_name = MFEM_KERNEL_NAME(KernelName);                  \
+      using KernelSignature = KernelType;                                      \
+      template <MFEM_PARAM_LIST P3>                                            \
+      static KernelSignature Kernel();                                         \
+      static KernelSignature Fallback(MFEM_PARAM_LIST P1);                     \
+      static KernelName &Get()                                                 \
+      { static KernelName table; return table;}                                \
+   }
+
+/// @brief Hashes variadic packs for which each type contained in the variadic
+/// pack has a specialization of `std::hash` available.
+///
+/// For example, packs containing int, bool, enum values, etc.
+template<typename ...KernelParameters>
+struct KernelDispatchKeyHash
+{
+private:
+   template<int N>
+   size_t operator()(std::tuple<KernelParameters...> value) const { return 0; }
+
+   // The hashing formula here is taken directly from the Boost library, with
+   // the magic number 0x9e3779b9 chosen to minimize hashing collisions.
+   template<std::size_t N, typename THead, typename... TTail>
+   size_t operator()(std::tuple<KernelParameters...> value) const
+   {
+      constexpr int Index = N - sizeof...(TTail) - 1;
+      auto lhs_hash = std::hash<THead>()(std::get<Index>(value));
+      auto rhs_hash = operator()<N, TTail...>(value);
+      return lhs_hash^(rhs_hash + 0x9e3779b9 + (lhs_hash<<6) + (lhs_hash>>2));
+   }
+public:
+   /// Returns the hash of the given @a value.
+   size_t operator()(std::tuple<KernelParameters...> value) const
+   {
+      return operator()<sizeof...(KernelParameters),KernelParameters...>(value);
+   }
+};
+
+namespace internal { template<typename... Types> struct KernelTypeList { }; }
+
+template<typename... T> class KernelDispatchTable { };
+
+template <typename Kernels,
+          typename Signature,
+          typename... Params,
+          typename... OptParams>
+class KernelDispatchTable<Kernels,
+         Signature,
+         internal::KernelTypeList<Params...>,
+         internal::KernelTypeList<OptParams...>>
+{
+   std::unordered_map<std::tuple<Params...>,
+       Signature,
+       KernelDispatchKeyHash<Params...>> table;
+
+public:
+   /// @brief Run the kernel with the given dispatch parameters and arguments.
+   ///
+   /// If a compile-time specialized version of the kernel with the given
+   /// parameters has been registered, it will be called. Otherwise, the
+   /// fallback kernel will be called.
+   template<typename... Args>
+   static void Run(Params... params, Args&&... args)
+   {
+      const auto &table = Kernels::Get().table;
+      const std::tuple<Params...> key = std::make_tuple(params...);
+      const auto it = table.find(key);
+      if (it != table.end())
+      {
+         it->second(std::forward<Args>(args)...);
+      }
+      else
+      {
+         ReportFallback(Kernels::Get().kernel_name, params...);
+         Kernels::Fallback(params...)(std::forward<Args>(args)...);
+      }
+   }
+
+   /// Register a specialized kernel for dispatch.
+   template <Params... PARAMS>
+   struct Specialization
+   {
+      // Version without optional parameters
+      static void Add()
+      {
+         std::tuple<Params...> param_tuple(PARAMS...);
+         Kernels::Get().table[param_tuple] =
+            Kernels:: template Kernel<PARAMS...>();
+      };
+      // Version with optional parameters
+      template <OptParams... OPT_PARAMS>
+      struct Opt
+      {
+         static void Add()
+         {
+            std::tuple<Params...> param_tuple(PARAMS...);
+            Kernels::Get().table[param_tuple] =
+               Kernels:: template Kernel<PARAMS..., OPT_PARAMS...>();
+         }
+      };
+   };
+};
+
+}
+
+#endif

fem/kernel_reporter.hpp

@@ -0,0 +1,82 @@
+// Copyright (c) 2010-2024, 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_KERNEL_REPORTER_HPP
+#define MFEM_KERNEL_REPORTER_HPP
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_REPORT_KERNELS
+
+#include "../general/globals.hpp"
+#include <set>
+#include <sstream>
+#include <string>
+
+#define MFEM_STR_(X) #X
+#define MFEM_STR(X) MFEM_STR_(X)
+#define MFEM_KERNEL_NAME(KernelName)  \
+   __FILE__ ":" MFEM_STR(__LINE__) " : " #KernelName
+
+namespace mfem
+{
+
+namespace
+{
+
+template <typename Last>
+static void Stringify_(std::ostream &o, Last &&arg)
+{
+   o << arg;
+}
+
+template <typename T1, typename T2, typename... Rest>
+static void Stringify_(std::ostream &o, T1 &&a1, T2 &&a2, Rest&&... rest)
+{
+   o << int(a1) << ",";
+   Stringify_(o, a2, rest...);
+}
+
+template <typename... Args>
+static std::string Stringify(Args&&... args)
+{
+   std::stringstream o;
+   Stringify_(o, args...);
+   return o.str();
+}
+
+} // namespace
+
+template <typename... Params>
+void ReportFallback(const std::string &kernel_name, Params&&... params)
+{
+   static std::set<std::string> reported_fallbacks;
+   const std::string requested_kernel =
+      kernel_name + "<" + Stringify(params...) + ">";
+   if (reported_fallbacks.find(requested_kernel) == reported_fallbacks.end())
+   {
+      reported_fallbacks.insert(requested_kernel);
+      mfem::err << "Fallback kernel. Requested "
+                << requested_kernel << std::endl;
+   }
+}
+
+} // namespace mfem
+
+#else // #ifdef MFEM_REPORT_KERNELS
+
+// No-op
+#define MFEM_KERNEL_NAME(KernelName) ""
+template <typename... T> void ReportFallback(T&&...) { }
+
+#endif
+
+#endif

fem/qinterp/det.cpp

@@ -27,12 +27,16 @@ namespace quadrature_interpolator
 {
 
 static void Det1D(const int NE,
+                  const real_t *b,
                   const real_t *g,
                   const real_t *x,
                   real_t *y,
                   const int d1d,
-                  const int q1d)
+                  const int q1d,
+                  Vector *d_buff = nullptr)
 {
+   MFEM_CONTRACT_VAR(b);
+   MFEM_CONTRACT_VAR(d_buff);
    const auto G = Reshape(g, q1d, d1d);
    const auto X = Reshape(x, d1d, NE);
 
@@ -59,8 +63,10 @@ static void Det2D(const int NE,
                   const real_t *x,
                   real_t *y,
                   const int d1d = 0,
-                  const int q1d = 0)
+                  const int q1d = 0,
+                  Vector *d_buff = nullptr)
 {
+   MFEM_CONTRACT_VAR(d_buff);
    static constexpr int SDIM = 2;
    static constexpr int NBZ = 1;
 
@@ -109,8 +115,11 @@ static void Det2DSurface(const int NE,
                          const real_t *x,
                          real_t *y,
                          const int d1d = 0,
-                         const int q1d = 0)
+                         const int q1d = 0,
+                         Vector *d_buff = nullptr)
 {
+   MFEM_CONTRACT_VAR(d_buff);
+
    static constexpr int SDIM = 3;
    static constexpr int NBZ = 1;
 
@@ -272,91 +281,63 @@ static void Det3D(const int NE,
    });
 }
 
-// Tensor-product evaluation of quadrature point determinants: dispatch
-// function.
-void TensorDeterminants(const int NE,
-                        const int vdim,
-                        const DofToQuad &maps,
-                        const Vector &e_vec,
-                        Vector &q_det,
-                        Vector &d_buff)
+void InitDetKernels()
 {
-   if (NE == 0) { return; }
-   const int dim = maps.FE->GetDim();
-   const int D1D = maps.ndof;
-   const int Q1D = maps.nqpt;
-   const real_t *B = maps.B.Read();
-   const real_t *G = maps.G.Read();
-   const real_t *X = e_vec.Read();
-   real_t *Y = q_det.Write();
-
-   const int id = (vdim<<8) | (D1D<<4) | Q1D;
-
-   if (dim == 1)
-   {
-      MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D,
-                  "Orders higher than " << DeviceDofQuadLimits::Get().MAX_D1D-1
-                  << " are not supported!");
-      MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D,
-                  "Quadrature rules with more than "
-                  << DeviceDofQuadLimits::Get().MAX_Q1D << " 1D points are not supported!");
-      Det1D(NE, G, X, Y, D1D, Q1D);
-      return;
-   }
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x222: return Det2D<2,2>(NE,B,G,X,Y);
-         case 0x223: return Det2D<2,3>(NE,B,G,X,Y);
-         case 0x224: return Det2D<2,4>(NE,B,G,X,Y);
-         case 0x226: return Det2D<2,6>(NE,B,G,X,Y);
-         case 0x234: return Det2D<3,4>(NE,B,G,X,Y);
-         case 0x236: return Det2D<3,6>(NE,B,G,X,Y);
-         case 0x244: return Det2D<4,4>(NE,B,G,X,Y);
-         case 0x246: return Det2D<4,6>(NE,B,G,X,Y);
-         case 0x256: return Det2D<5,6>(NE,B,G,X,Y);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_D1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_Q1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            if (vdim == 2) { Det2D(NE,B,G,X,Y,D1D,Q1D); }
-            else if (vdim == 3) { Det2DSurface(NE,B,G,X,Y,D1D,Q1D); }
-            else { MFEM_ABORT("Invalid space dimension."); }
-            return;
-         }
-      }
-   }
-   if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x324: return Det3D<2,4>(NE,B,G,X,Y);
-         case 0x333: return Det3D<3,3>(NE,B,G,X,Y);
-         case 0x335: return Det3D<3,5>(NE,B,G,X,Y);
-         case 0x336: return Det3D<3,6>(NE,B,G,X,Y);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_DET_1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_DET_1D;
-            // Highest orders that fit in shared memory
-            if (D1D <= MD && Q1D <= MQ)
-            { return Det3D<0,0,true>(NE,B,G,X,Y,D1D,Q1D); }
-            // Last fall-back will use global memory
-            return Det3D<0,0,false>(
-                      NE,B,G,X,Y,D1D,Q1D,&d_buff);
-         }
-      }
-   }
-   MFEM_ABORT("Kernel " << std::hex << id << std::dec << " not supported yet");
+   using k = QuadratureInterpolator::DetKernels;
+   // 2D
+   k::Specialization<2,2,2,2>::Add();
+   k::Specialization<2,2,2,3>::Add();
+   k::Specialization<2,2,2,4>::Add();
+   k::Specialization<2,2,2,6>::Add();
+   k::Specialization<2,2,3,4>::Add();
+   k::Specialization<2,2,3,6>::Add();
+   k::Specialization<2,2,4,4>::Add();
+   k::Specialization<2,2,4,6>::Add();
+   k::Specialization<2,2,5,6>::Add();
+   // 3D
+   k::Specialization<3,3,2,4>::Add();
+   k::Specialization<3,3,3,3>::Add();
+   k::Specialization<3,3,3,5>::Add();
+   k::Specialization<3,3,3,6>::Add();
 }
 
 } // namespace quadrature_interpolator
 
 } // namespace internal
 
+/// @cond Suppress_Doxygen_warnings
+
+namespace
+{
+using DetKernel = QuadratureInterpolator::DetKernelType;
+}
+
+template<int DIM, int SDIM, int D1D, int Q1D>
+DetKernel QuadratureInterpolator::DetKernels::Kernel()
+{
+   if (DIM == 1) { return internal::quadrature_interpolator::Det1D; }
+   else if (DIM == 2 && SDIM == 2) { return internal::quadrature_interpolator::Det2D<D1D, Q1D>; }
+   else if (DIM == 2 && SDIM == 3) { return internal::quadrature_interpolator::Det2DSurface<D1D, Q1D>; }
+   else if (DIM == 3) { return internal::quadrature_interpolator::Det3D<D1D, Q1D>; }
+   else { MFEM_ABORT(""); }
+}
+
+DetKernel QuadratureInterpolator::DetKernels::Fallback(
+   int DIM, int SDIM, int D1D, int Q1D)
+{
+   if (DIM == 1) { return internal::quadrature_interpolator::Det1D; }
+   else if (DIM == 2 && SDIM == 2) { return internal::quadrature_interpolator::Det2D; }
+   else if (DIM == 2 && SDIM == 3) { return internal::quadrature_interpolator::Det2DSurface; }
+   else if (DIM == 3)
+   {
+      const int MD = DeviceDofQuadLimits::Get().MAX_DET_1D;
+      const int MQ = DeviceDofQuadLimits::Get().MAX_DET_1D;
+      if (D1D <= MD && Q1D <= MQ) { return internal::quadrature_interpolator::Det3D<0,0,true>; }
+      else { return internal::quadrature_interpolator::Det3D<0,0,false>; }
+   }
+   else { MFEM_ABORT(""); }
+}
+
+/// @endcond
+
 } // namespace mfem

fem/qinterp/dispatch.hpp

@@ -1,64 +0,0 @@
-// Copyright (c) 2010-2024, 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.
-
-// Internal header, included only by .cpp files
-
-#include "../quadinterpolator.hpp"
-
-namespace mfem
-{
-
-namespace internal
-{
-
-namespace quadrature_interpolator
-{
-
-// Tensor-product evaluation of quadrature point values: dispatch function.
-template<QVectorLayout VL>
-void TensorValues(const int NE,
-                  const int vdim,
-                  const DofToQuad &maps,
-                  const Vector &e_vec,
-                  Vector &q_val);
-
-// Tensor-product evaluation of quadrature point derivatives: dispatch function.
-template<QVectorLayout VL>
-void TensorDerivatives(const int NE,
-                       const int vdim,
-                       const DofToQuad &maps,
-                       const Vector &e_vec,
-                       Vector &q_der);
-
-// Tensor-product evaluation of quadrature point physical derivatives: dispatch
-// function.
-template<QVectorLayout VL>
-void TensorPhysDerivatives(const int NE,
-                           const int vdim,
-                           const DofToQuad &maps,
-                           const GeometricFactors &geom,
-                           const Vector &e_vec,
-                           Vector &q_der);
-
-// Tensor-product evaluation of quadrature point determinants: dispatch
-// function.
-void TensorDeterminants(const int NE,
-                        const int vdim,
-                        const DofToQuad &maps,
-                        const Vector &e_vec,
-                        Vector &q_det,
-                        Vector &d_buff);
-
-} // namespace quadrature_interpolator
-
-} // namespace internal
-
-} // namespace mfem

fem/qinterp/eval.hpp

@@ -12,6 +12,9 @@
 // Internal header, included only by .cpp files.
 // Template function implementations.
 
+#ifndef MFEM_QUADINTERP_EVAL
+#define MFEM_QUADINTERP_EVAL
+
 #include "../quadinterpolator.hpp"
 #include "../../general/forall.hpp"
 #include "../../linalg/dtensor.hpp"
@@ -63,7 +66,7 @@ static void Values1D(const int NE,
 // Template compute kernel for Values in 2D: tensor product version.
 template<QVectorLayout Q_LAYOUT,
          int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
-         int T_NBZ = 1, int MAX_D1D = 0, int MAX_Q1D = 0>
+         int T_NBZ = 1>
 static void Values2D(const int NE,
                      const real_t *b_,
                      const real_t *x_,
@@ -193,4 +196,21 @@ static void Values3D(const int NE,
 
 } // namespace internal
 
+/// @cond Suppress_Doxygen_warnings
+
+template<int DIM, QVectorLayout Q_LAYOUT,
+         int VDIM, int D1D, int Q1D, int NBZ>
+QuadratureInterpolator::TensorEvalKernelType
+QuadratureInterpolator::TensorEvalKernels::Kernel()
+{
+   if (DIM == 1) { return internal::quadrature_interpolator::Values1D<Q_LAYOUT>; }
+   else if (DIM == 2) { return internal::quadrature_interpolator::Values2D<Q_LAYOUT, VDIM, D1D, Q1D, NBZ>; }
+   else if (DIM == 3) { return internal::quadrature_interpolator::Values3D<Q_LAYOUT, VDIM, D1D, Q1D>; }
+   else { MFEM_ABORT(""); }
+}
+
+/// @endcond
+
 } // namespace mfem
+
+#endif

fem/qinterp/eval_by_nodes.cpp

@@ -10,143 +10,75 @@
 // CONTRIBUTING.md for details.
 
 #include "../quadinterpolator.hpp"
-#include "dispatch.hpp"
 #include "eval.hpp"
 
 namespace mfem
 {
-
 namespace internal
 {
-
 namespace quadrature_interpolator
 {
 
-// Tensor-product evaluation of quadrature point values: dispatch function.
-// Instantiation for the case QVectorLayout::byNODES.
-template<>
-void TensorValues<QVectorLayout::byNODES>(const int NE,
-                                          const int vdim,
-                                          const DofToQuad &maps,
-                                          const Vector &e_vec,
-                                          Vector &q_val)
+void InitEvalByNodesKernels()
 {
-   if (NE == 0) { return; }
-   const int dim = maps.FE->GetDim();
-   const int D1D = maps.ndof;
-   const int Q1D = maps.nqpt;
-   const real_t *B = maps.B.Read();
-   const real_t *X = e_vec.Read();
-   real_t *Y = q_val.Write();
+   using k = QuadratureInterpolator::TensorEvalKernels;
 
-   constexpr QVectorLayout L = QVectorLayout::byNODES;
+   // 2D
+   k::Specialization<2,QVectorLayout::byNODES,1,3,3>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,1,2,4>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,1,3,2>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,1,3,4>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,1,4,3>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,1,4,4>::Opt<1>::Add();
 
-   const int id = (vdim<<8) | (D1D<<4) | Q1D;
+   k::Specialization<2,QVectorLayout::byNODES,2,2,2>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,2,3>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,2,4>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,2,5>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,2,6>::Opt<1>::Add();
 
-   if (dim == 1)
-   {
-      MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D,
-                  "Orders higher than " << DeviceDofQuadLimits::Get().MAX_D1D-1
-                  << " are not supported!");
-      MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D,
-                  "Quadrature rules with more than "
-                  << DeviceDofQuadLimits::Get().MAX_Q1D << " 1D points are not supported!");
-      Values1D<L>(NE, B, X, Y, vdim, D1D, Q1D);
-      return;
-   }
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x133: return Values2D<L,1,3,3>(NE,B,X,Y);
-         case 0x124: return Values2D<L,1,2,4>(NE,B,X,Y);
-         case 0x132: return Values2D<L,1,3,2>(NE,B,X,Y);
-         case 0x134: return Values2D<L,1,3,4>(NE,B,X,Y);
-         case 0x143: return Values2D<L,1,4,3>(NE,B,X,Y);
-         case 0x144: return Values2D<L,1,4,4>(NE,B,X,Y);
+   k::Specialization<2,QVectorLayout::byNODES,2,3,3>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,3,4>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,3,6>::Opt<1>::Add();
 
-         case 0x222: return Values2D<L,2,2,2>(NE,B,X,Y);
-         case 0x223: return Values2D<L,2,2,3>(NE,B,X,Y);
-         case 0x224: return Values2D<L,2,2,4>(NE,B,X,Y);
-         case 0x225: return Values2D<L,2,2,5>(NE,B,X,Y);
-         case 0x226: return Values2D<L,2,2,6>(NE,B,X,Y);
+   k::Specialization<2,QVectorLayout::byNODES,2,4,3>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,4,4>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,4,5>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,4,6>::Opt<1>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,2,4,7>::Opt<1>::Add();
 
-         case 0x233: return Values2D<L,2,3,3>(NE,B,X,Y);
-         case 0x234: return Values2D<L,2,3,4>(NE,B,X,Y);
-         case 0x236: return Values2D<L,2,3,6>(NE,B,X,Y);
+   k::Specialization<2,QVectorLayout::byNODES,2,5,6>::Opt<1>::Add();
 
-         case 0x243: return Values2D<L,2,4,3>(NE,B,X,Y);
-         case 0x244: return Values2D<L,2,4,4>(NE,B,X,Y);
-         case 0x245: return Values2D<L,2,4,5>(NE,B,X,Y);
-         case 0x246: return Values2D<L,2,4,6>(NE,B,X,Y);
-         case 0x247: return Values2D<L,2,4,7>(NE,B,X,Y);
+   // 3D
+   k::Specialization<3,QVectorLayout::byNODES,1,2,4>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,1,3,3>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,1,3,4>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,1,3,6>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,1,4,3>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,1,4,4>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,1,4,8>::Opt<1>::Add();
 
-         case 0x256: return Values2D<L,2,5,6>(NE,B,X,Y);
+   k::Specialization<3,QVectorLayout::byNODES,2,2,2>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,2,2,3>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,2,3,4>::Opt<1>::Add();
 
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_D1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_Q1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Values2D<L>(NE,B,X,Y,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x124: return Values3D<L,1,2,4>(NE,B,X,Y);
-         case 0x133: return Values3D<L,1,3,3>(NE,B,X,Y);
-         case 0x134: return Values3D<L,1,3,4>(NE,B,X,Y);
-         case 0x136: return Values3D<L,1,3,6>(NE,B,X,Y);
-         case 0x143: return Values3D<L,1,4,3>(NE,B,X,Y);
-         case 0x144: return Values3D<L,1,4,4>(NE,B,X,Y);
-         case 0x148: return Values3D<L,1,4,8>(NE,B,X,Y);
+   k::Specialization<3,QVectorLayout::byNODES,3,2,3>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,2,4>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,2,5>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,2,6>::Opt<1>::Add();
 
-         case 0x222: return Values3D<L,2,2,2>(NE,B,X,Y);
-         case 0x223: return Values3D<L,2,2,3>(NE,B,X,Y);
-         case 0x234: return Values3D<L,2,3,4>(NE,B,X,Y);
+   k::Specialization<3,QVectorLayout::byNODES,3,3,3>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,3,4>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,3,5>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,3,6>::Opt<1>::Add();
 
-         case 0x323: return Values3D<L,3,2,3>(NE,B,X,Y);
-         case 0x324: return Values3D<L,3,2,4>(NE,B,X,Y);
-         case 0x325: return Values3D<L,3,2,5>(NE,B,X,Y);
-         case 0x326: return Values3D<L,3,2,6>(NE,B,X,Y);
-
-         case 0x333: return Values3D<L,3,3,3>(NE,B,X,Y);
-         case 0x334: return Values3D<L,3,3,4>(NE,B,X,Y);
-         case 0x335: return Values3D<L,3,3,5>(NE,B,X,Y);
-         case 0x336: return Values3D<L,3,3,6>(NE,B,X,Y);
-
-         case 0x343: return Values3D<L,3,4,3>(NE,B,X,Y);
-         case 0x344: return Values3D<L,3,4,4>(NE,B,X,Y);
-         case 0x346: return Values3D<L,3,4,6>(NE,B,X,Y);
-         case 0x347: return Values3D<L,3,4,7>(NE,B,X,Y);
-         case 0x348: return Values3D<L,3,4,8>(NE,B,X,Y);
-
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Values3D<L>(NE,B,X,Y,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
-   MFEM_ABORT("Kernel not supported yet");
+   k::Specialization<3,QVectorLayout::byNODES,3,4,3>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,4,4>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,4,6>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,4,7>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,3,4,8>::Opt<1>::Add();
 }
 
 } // namespace quadrature_interpolator
-
 } // namespace internal
-
 } // namespace mfem

fem/qinterp/eval_by_vdim.cpp

@@ -10,117 +10,45 @@
 // CONTRIBUTING.md for details.
 
 #include "../quadinterpolator.hpp"
-#include "dispatch.hpp"
 #include "eval.hpp"
 
 namespace mfem
 {
-
 namespace internal
 {
-
 namespace quadrature_interpolator
 {
 
-// Tensor-product evaluation of quadrature point values: dispatch function.
-// Instantiation for the case QVectorLayout::byVDIM.
-template<>
-void TensorValues<QVectorLayout::byVDIM>(const int NE,
-                                         const int vdim,
-                                         const DofToQuad &maps,
-                                         const Vector &e_vec,
-                                         Vector &q_val)
+void InitEvalByVDimKernels()
 {
-   if (NE == 0) { return; }
-   const int dim = maps.FE->GetDim();
-   const int D1D = maps.ndof;
-   const int Q1D = maps.nqpt;
-   const real_t *B = maps.B.Read();
-   const real_t *X = e_vec.Read();
-   real_t *Y = q_val.Write();
+   using k = QuadratureInterpolator::TensorEvalKernels;
+   // 2D
+   k::Specialization<2,QVectorLayout::byVDIM,1,2,4>::Opt<8>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,1,3,6>::Opt<4>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,1,4,8>::Opt<2>::Add();
 
-   constexpr QVectorLayout L = QVectorLayout::byVDIM;
+   k::Specialization<2,QVectorLayout::byVDIM,2,2,4>::Opt<8>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,2,3,4>::Opt<8>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,2,3,6>::Opt<4>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,2,4,8>::Opt<2>::Add();
+   // 3D
+   k::Specialization<3,QVectorLayout::byVDIM,1,2,4>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,1,3,6>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,1,4,8>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,2,4>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,3,6>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,4,8>::Opt<1>::Add();
 
-   const int id = (vdim<<8) | (D1D<<4) | Q1D;
-
-   if (dim == 1)
-   {
-      MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D,
-                  "Orders higher than " << DeviceDofQuadLimits::Get().MAX_D1D-1
-                  << " are not supported!");
-      MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D,
-                  "Quadrature rules with more than "
-                  << DeviceDofQuadLimits::Get().MAX_Q1D << " 1D points are not supported!");
-      Values1D<L>(NE, B, X, Y, vdim, D1D, Q1D);
-      return;
-   }
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x124: return Values2D<L,1,2,4,8>(NE,B,X,Y);
-         case 0x136: return Values2D<L,1,3,6,4>(NE,B,X,Y);
-         case 0x148: return Values2D<L,1,4,8,2>(NE,B,X,Y);
-
-         case 0x224: return Values2D<L,2,2,4,8>(NE,B,X,Y);
-         case 0x234: return Values2D<L,2,3,4,8>(NE,B,X,Y);
-         case 0x236: return Values2D<L,2,3,6,4>(NE,B,X,Y);
-         case 0x248: return Values2D<L,2,4,8,2>(NE,B,X,Y);
-
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_D1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_Q1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Values2D<L>(NE,B,X,Y,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x124: return Values3D<L,1,2,4>(NE,B,X,Y);
-         case 0x136: return Values3D<L,1,3,6>(NE,B,X,Y);
-         case 0x148: return Values3D<L,1,4,8>(NE,B,X,Y);
-
-         case 0x324: return Values3D<L,3,2,4>(NE,B,X,Y);
-         case 0x336: return Values3D<L,3,3,6>(NE,B,X,Y);
-         case 0x348: return Values3D<L,3,4,8>(NE,B,X,Y);
-
-         // Used for LOR batched assembly
-         case 0x322: return Values3D<L,3,2,2>(NE,B,X,Y);
-         case 0x333: return Values3D<L,3,3,3>(NE,B,X,Y);
-         case 0x344: return Values3D<L,3,4,4>(NE,B,X,Y);
-         case 0x355: return Values3D<L,3,5,5>(NE,B,X,Y);
-         case 0x366: return Values3D<L,3,6,6>(NE,B,X,Y);
-         case 0x377: return Values3D<L,3,7,7>(NE,B,X,Y);
-         case 0x388: return Values3D<L,3,8,8>(NE,B,X,Y);
-         case 0x399: return Values3D<L,3,9,9>(NE,B,X,Y);
-
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Values3D<L>(NE,B,X,Y,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
-   MFEM_ABORT("Kernel not supported yet");
+   k::Specialization<3,QVectorLayout::byVDIM,3,2,2>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,3,3>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,4,4>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,5,5>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,6,6>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,7,7>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,8,8>::Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,3,9,9>::Opt<1>::Add();
 }
 
 } // namespace quadrature_interpolator
-
 } // namespace internal
-
 } // namespace mfem

fem/qinterp/grad.hpp

@@ -12,6 +12,9 @@
 // Internal header, included only by .cpp files.
 // Template function implementations.
 
+#ifndef MFEM_QUADINTERP_GRAD
+#define MFEM_QUADINTERP_GRAD
+
 #include "../quadinterpolator.hpp"
 #include "../../general/forall.hpp"
 #include "../../linalg/dtensor.hpp"
@@ -29,6 +32,7 @@ namespace quadrature_interpolator
 
 template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS>
 static void Derivatives1D(const int NE,
+                          const real_t *b_,
                           const real_t *g_,
                           const real_t *j_,
                           const real_t *x_,
@@ -38,6 +42,7 @@ static void Derivatives1D(const int NE,
                           const int d1d,
                           const int q1d)
 {
+   MFEM_CONTRACT_VAR(b_);
    const auto g = Reshape(g_, q1d, d1d);
    const auto j = Reshape(j_, q1d, sdim, NE);
    const auto x = Reshape(x_, d1d, vdim, NE);
@@ -232,6 +237,7 @@ static void Derivatives3D(const int NE,
                           const real_t *j_,
                           const real_t *x_,
                           real_t *y_,
+                          const int sdim = 3,
                           const int vdim = 0,
                           const int d1d = 0,
                           const int q1d = 0)
@@ -370,4 +376,21 @@ static void Derivatives3D(const int NE,
 
 } // namespace internal
 
+/// @cond Suppress_Doxygen_warnings
+
+template<int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
+         int VDIM, int D1D, int Q1D, int NBZ>
+QuadratureInterpolator::GradKernelType
+QuadratureInterpolator::GradKernels::Kernel()
+{
+   if (DIM == 1) { return internal::quadrature_interpolator::Derivatives1D<Q_LAYOUT, GRAD_PHYS>; }
+   else if (DIM == 2) { return internal::quadrature_interpolator::Derivatives2D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, Q1D, NBZ>; }
+   else if (DIM == 3) { return internal::quadrature_interpolator::Derivatives3D<Q_LAYOUT, GRAD_PHYS, VDIM, D1D, Q1D>; }
+   else { MFEM_ABORT(""); }
+}
+
+/// @endcond
+
 } // namespace mfem
+
+#endif

fem/qinterp/grad_by_nodes.cpp

@@ -9,128 +9,69 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#include "dispatch.hpp"
+#include "../quadinterpolator.hpp"
 #include "grad.hpp"
 
 namespace mfem
 {
-
 namespace internal
 {
-
 namespace quadrature_interpolator
 {
 
-// Tensor-product evaluation of quadrature point derivatives: dispatch function.
-// Instantiation for the case QVectorLayout::byNODES.
-template<>
-void TensorDerivatives<QVectorLayout::byNODES>(const int NE,
-                                               const int vdim,
-                                               const DofToQuad &maps,
-                                               const Vector &e_vec,
-                                               Vector &q_der)
+template <bool P>
+void InitGradByNodesKernels()
 {
-   if (NE == 0) { return; }
-   const int dim = maps.FE->GetDim();
-   const int D1D = maps.ndof;
-   const int Q1D = maps.nqpt;
-   const real_t *B = maps.B.Read();
-   const real_t *G = maps.G.Read();
-   const real_t *J = nullptr; // not used in DERIVATIVES (non-GRAD_PHYS) mode
-   const real_t *X = e_vec.Read();
-   real_t *Y = q_der.Write();
+   using k = QuadratureInterpolator::GradKernels;
+   // 2D
+   k::Specialization<2,QVectorLayout::byNODES,P,1,3,3>::template Opt<16>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,1,3,4>::template Opt<16>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,1,4,3>::template Opt<16>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,1,4,4>::template Opt<16>::Add();
 
-   constexpr QVectorLayout L = QVectorLayout::byNODES;
-   constexpr bool P = false; // GRAD_PHYS
+   k::Specialization<2,QVectorLayout::byNODES,P,2,2,2>::template Opt<16>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,2,3>::template Opt<8>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,2,4>::template Opt<4>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,2,5>::template Opt<4>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,2,6>::template Opt<2>::Add();
 
-   const int id = (vdim<<8) | (D1D<<4) | Q1D;
+   k::Specialization<2,QVectorLayout::byNODES,P,2,3,3>::template Opt<2>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,3,4>::template Opt<4>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,4,3>::template Opt<4>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,3,6>::template Opt<2>::Add();
 
-   if (dim == 1)
-   {
-      return Derivatives1D<L,P>(NE,G,J,X,Y,dim,vdim,D1D,Q1D);
-   }
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x133: return Derivatives2D<L,P,1,3,3,16>(NE,B,G,J,X,Y);
-         case 0x134: return Derivatives2D<L,P,1,3,4,16>(NE,B,G,J,X,Y);
-         case 0x143: return Derivatives2D<L,P,1,4,3,16>(NE,B,G,J,X,Y);
-         case 0x144: return Derivatives2D<L,P,1,4,4,16>(NE,B,G,J,X,Y);
+   k::Specialization<2,QVectorLayout::byNODES,P,2,4,4>::template Opt<2>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,4,5>::template Opt<2>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,4,6>::template Opt<2>::Add();
+   k::Specialization<2,QVectorLayout::byNODES,P,2,4,7>::template Opt<2>::Add();
 
-         case 0x222: return Derivatives2D<L,P,2,2,2,16>(NE,B,G,J,X,Y);
-         case 0x223: return Derivatives2D<L,P,2,2,3,8>(NE,B,G,J,X,Y);
-         case 0x224: return Derivatives2D<L,P,2,2,4,4>(NE,B,G,J,X,Y);
-         case 0x225: return Derivatives2D<L,P,2,2,5,4>(NE,B,G,J,X,Y);
-         case 0x226: return Derivatives2D<L,P,2,2,6,2>(NE,B,G,J,X,Y);
+   k::Specialization<2,QVectorLayout::byNODES,P,2,5,6>::template Opt<2>::Add();
+   // 3D
+   k::Specialization<3,QVectorLayout::byNODES,P,1,2,4>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,1,3,3>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,1,3,4>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,1,3,6>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,1,4,4>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,1,4,8>::template Opt<1>::Add();
 
-         case 0x233: return Derivatives2D<L,P,2,3,3,2>(NE,B,G,J,X,Y);
-         case 0x234: return Derivatives2D<L,P,2,3,4,4>(NE,B,G,J,X,Y);
-         case 0x243: return Derivatives2D<L,P,2,4,3,4>(NE,B,G,J,X,Y);
-         case 0x236: return Derivatives2D<L,P,2,3,6,2>(NE,B,G,J,X,Y);
+   k::Specialization<3,QVectorLayout::byNODES,P,3,2,3>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,2,4>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,2,5>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,2,6>::template Opt<1>::Add();
 
-         case 0x244: return Derivatives2D<L,P,2,4,4,2>(NE,B,G,J,X,Y);
-         case 0x245: return Derivatives2D<L,P,2,4,5,2>(NE,B,G,J,X,Y);
-         case 0x246: return Derivatives2D<L,P,2,4,6,2>(NE,B,G,J,X,Y);
-         case 0x247: return Derivatives2D<L,P,2,4,7,2>(NE,B,G,J,X,Y);
-
-         case 0x256: return Derivatives2D<L,P,2,5,6,2>(NE,B,G,J,X,Y);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_D1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_Q1D;
-            if (D1D > MD || Q1D > MQ)
-            {
-               MFEM_ABORT("");
-            }
-            Derivatives2D<L,P>(NE,B,G,J,X,Y,dim,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x124: return Derivatives3D<L,P,1,2,4>(NE,B,G,J,X,Y);
-         case 0x133: return Derivatives3D<L,P,1,3,3>(NE,B,G,J,X,Y);
-         case 0x134: return Derivatives3D<L,P,1,3,4>(NE,B,G,J,X,Y);
-         case 0x136: return Derivatives3D<L,P,1,3,6>(NE,B,G,J,X,Y);
-         case 0x144: return Derivatives3D<L,P,1,4,4>(NE,B,G,J,X,Y);
-         case 0x148: return Derivatives3D<L,P,1,4,8>(NE,B,G,J,X,Y);
-
-         case 0x323: return Derivatives3D<L,P,3,2,3>(NE,B,G,J,X,Y);
-         case 0x324: return Derivatives3D<L,P,3,2,4>(NE,B,G,J,X,Y);
-         case 0x325: return Derivatives3D<L,P,3,2,5>(NE,B,G,J,X,Y);
-         case 0x326: return Derivatives3D<L,P,3,2,6>(NE,B,G,J,X,Y);
-
-         case 0x333: return Derivatives3D<L,P,3,3,3>(NE,B,G,J,X,Y);
-         case 0x334: return Derivatives3D<L,P,3,3,4>(NE,B,G,J,X,Y);
-         case 0x335: return Derivatives3D<L,P,3,3,5>(NE,B,G,J,X,Y);
-         case 0x336: return Derivatives3D<L,P,3,3,6>(NE,B,G,J,X,Y);
-         case 0x344: return Derivatives3D<L,P,3,4,4>(NE,B,G,J,X,Y);
-         case 0x346: return Derivatives3D<L,P,3,4,6>(NE,B,G,J,X,Y);
-         case 0x347: return Derivatives3D<L,P,3,4,7>(NE,B,G,J,X,Y);
-         case 0x348: return Derivatives3D<L,P,3,4,8>(NE,B,G,J,X,Y);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Derivatives3D<L,P>(NE,B,G,J,X,Y,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
-   MFEM_ABORT("Kernel not supported yet");
+   k::Specialization<3,QVectorLayout::byNODES,P,3,3,3>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,3,4>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,3,5>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,3,6>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,4,4>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,4,6>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,4,7>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byNODES,P,3,4,8>::template Opt<1>::Add();
 }
 
+template void InitGradByNodesKernels<true>();
+template void InitGradByNodesKernels<false>();
+
 } // namespace quadrature_interpolator
-
 } // namespace internal
-
 } // namespace mfem

fem/qinterp/grad_by_vdim.cpp

@@ -9,100 +9,41 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#include "dispatch.hpp"
+#include "../quadinterpolator.hpp"
 #include "grad.hpp"
 
 namespace mfem
 {
-
 namespace internal
 {
-
 namespace quadrature_interpolator
 {
 
-// Tensor-product evaluation of quadrature point derivatives: dispatch function.
-// Instantiation for the case QVectorLayout::byVDIM.
-template<>
-void TensorDerivatives<QVectorLayout::byVDIM>(const int NE,
-                                              const int vdim,
-                                              const DofToQuad &maps,
-                                              const Vector &e_vec,
-                                              Vector &q_der)
+template <bool P>
+void InitGradByVDimKernels()
 {
-   if (NE == 0) { return; }
-   const int dim = maps.FE->GetDim();
-   const int D1D = maps.ndof;
-   const int Q1D = maps.nqpt;
-   const real_t *B = maps.B.Read();
-   const real_t *G = maps.G.Read();
-   const real_t *J = nullptr; // not used in DERIVATIVES (non-GRAD_PHYS) mode
-   const real_t *X = e_vec.Read();
-   real_t *Y = q_der.Write();
+   using k = QuadratureInterpolator::GradKernels;
+   // 2D
+   k::Specialization<2,QVectorLayout::byVDIM,P,1,3,4>::template Opt<8>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,P,1,4,6>::template Opt<4>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,P,1,5,8>::template Opt<2>::Add();
 
-   constexpr QVectorLayout L = QVectorLayout::byVDIM;
-   constexpr bool P = false; // GRAD_PHYS
-
-   const int id = (vdim<<8) | (D1D<<4) | Q1D;
-
-   if (dim == 1)
-   {
-      return Derivatives1D<L,P>(NE,G,J,X,Y,dim,vdim,D1D,Q1D);
-   }
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x134: return Derivatives2D<L,P,1,3,4,8>(NE,B,G,J,X,Y);
-         case 0x146: return Derivatives2D<L,P,1,4,6,4>(NE,B,G,J,X,Y);
-         case 0x158: return Derivatives2D<L,P,1,5,8,2>(NE,B,G,J,X,Y);
-
-         case 0x234: return Derivatives2D<L,P,2,3,4,8>(NE,B,G,J,X,Y);
-         case 0x246: return Derivatives2D<L,P,2,4,6,4>(NE,B,G,J,X,Y);
-         case 0x258: return Derivatives2D<L,P,2,5,8,2>(NE,B,G,J,X,Y);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_D1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_Q1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Derivatives2D<L,P>(NE,B,G,J,X,Y,dim,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x134: return Derivatives3D<L,P,1,3,4>(NE,B,G,J,X,Y);
-         case 0x146: return Derivatives3D<L,P,1,4,6>(NE,B,G,J,X,Y);
-         case 0x158: return Derivatives3D<L,P,1,5,8>(NE,B,G,J,X,Y);
-
-         case 0x334: return Derivatives3D<L,P,3,3,4>(NE,B,G,J,X,Y);
-         case 0x346: return Derivatives3D<L,P,3,4,6>(NE,B,G,J,X,Y);
-         case 0x358: return Derivatives3D<L,P,3,5,8>(NE,B,G,J,X,Y);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Derivatives3D<L,P>(NE,B,G,J,X,Y,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
-   MFEM_ABORT("Kernel not supported yet");
+   k::Specialization<2,QVectorLayout::byVDIM,P,2,3,3>::template Opt<8>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,P,2,3,4>::template Opt<8>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,P,2,4,6>::template Opt<4>::Add();
+   k::Specialization<2,QVectorLayout::byVDIM,P,2,5,8>::template Opt<2>::Add();
+   // 3D
+   k::Specialization<3,QVectorLayout::byVDIM,P,1,3,4>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,P,1,4,6>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,P,1,5,8>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,P,3,3,4>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,P,3,4,6>::template Opt<1>::Add();
+   k::Specialization<3,QVectorLayout::byVDIM,P,3,5,8>::template Opt<1>::Add();
 }
 
+template void InitGradByVDimKernels<true>();
+template void InitGradByVDimKernels<false>();
+
 } // namespace quadrature_interpolator
-
 } // namespace internal
-
 } // namespace mfem

fem/qinterp/grad_phys_by_nodes.cpp

@@ -1,123 +0,0 @@
-// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "dispatch.hpp"
-#include "grad.hpp"
-
-namespace mfem
-{
-
-namespace internal
-{
-
-namespace quadrature_interpolator
-{
-
-// Tensor-product evaluation of quadrature point physical derivatives: dispatch
-// function.
-// Instantiation for the case QVectorLayout::byNODES.
-template<>
-void TensorPhysDerivatives<QVectorLayout::byNODES>(const int NE,
-                                                   const int vdim,
-                                                   const DofToQuad &maps,
-                                                   const GeometricFactors &geom,
-                                                   const Vector &e_vec,
-                                                   Vector &q_der)
-{
-   if (NE == 0) { return; }
-   const int dim = maps.FE->GetDim();
-   const int D1D = maps.ndof;
-   const int Q1D = maps.nqpt;
-
-   const int sdim = geom.mesh->SpaceDimension();
-
-   const real_t *B = maps.B.Read();
-   const real_t *G = maps.G.Read();
-   const real_t *J = geom.J.Read();
-   const real_t *X = e_vec.Read();
-   real_t *Y = q_der.Write();
-
-   constexpr QVectorLayout L = QVectorLayout::byNODES;
-   constexpr bool P = true; // GRAD_PHYS
-
-   const int id = (vdim<<8) | (D1D<<4) | Q1D;
-
-   if (dim == 1)
-   {
-      return Derivatives1D<L,P>(NE,G,J,X,Y,sdim,vdim,D1D,Q1D);
-   }
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x133: return Derivatives2D<L,P,1,3,3,8>(NE,B,G,J,X,Y,sdim);
-         case 0x134: return Derivatives2D<L,P,1,3,4,8>(NE,B,G,J,X,Y,sdim);
-         case 0x143: return Derivatives2D<L,P,1,4,3,4>(NE,B,G,J,X,Y,sdim);
-         case 0x144: return Derivatives2D<L,P,1,4,4,4>(NE,B,G,J,X,Y,sdim);
-         case 0x146: return Derivatives2D<L,P,1,4,6,4>(NE,B,G,J,X,Y,sdim);
-         case 0x158: return Derivatives2D<L,P,1,5,8,2>(NE,B,G,J,X,Y,sdim);
-
-         case 0x233: return Derivatives2D<L,P,2,3,3,8>(NE,B,G,J,X,Y,sdim);
-         case 0x234: return Derivatives2D<L,P,2,3,4,8>(NE,B,G,J,X,Y,sdim);
-         case 0x243: return Derivatives2D<L,P,2,4,3,4>(NE,B,G,J,X,Y,sdim);
-         case 0x244: return Derivatives2D<L,P,2,4,4,4>(NE,B,G,J,X,Y,sdim);
-         case 0x246: return Derivatives2D<L,P,2,4,6,4>(NE,B,G,J,X,Y,sdim);
-         case 0x258: return Derivatives2D<L,P,2,5,8,2>(NE,B,G,J,X,Y,sdim);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_D1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_Q1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Derivatives2D<L,P>(NE,B,G,J,X,Y,sdim,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x133: return Derivatives3D<L,P,1,3,3>(NE,B,G,J,X,Y);
-         case 0x134: return Derivatives3D<L,P,1,3,4>(NE,B,G,J,X,Y);
-         case 0x144: return Derivatives3D<L,P,1,4,4>(NE,B,G,J,X,Y);
-         case 0x146: return Derivatives3D<L,P,1,4,6>(NE,B,G,J,X,Y);
-         case 0x158: return Derivatives3D<L,P,1,5,8>(NE,B,G,J,X,Y);
-
-         case 0x333: return Derivatives3D<L,P,3,3,3>(NE,B,G,J,X,Y);
-         case 0x334: return Derivatives3D<L,P,3,3,4>(NE,B,G,J,X,Y);
-         case 0x344: return Derivatives3D<L,P,3,4,4>(NE,B,G,J,X,Y);
-         case 0x346: return Derivatives3D<L,P,3,4,6>(NE,B,G,J,X,Y);
-         case 0x358: return Derivatives3D<L,P,3,5,8>(NE,B,G,J,X,Y);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Derivatives3D<L,P>(NE,B,G,J,X,Y,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
-   MFEM_ABORT("Unknown kernel");
-}
-
-} // namespace quadrature_interpolator
-
-} // namespace internal
-
-} // namespace mfem

fem/qinterp/grad_phys_by_vdim.cpp

@@ -1,114 +0,0 @@
-// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "dispatch.hpp"
-#include "grad.hpp"
-
-namespace mfem
-{
-
-namespace internal
-{
-
-namespace quadrature_interpolator
-{
-
-// Tensor-product evaluation of quadrature point physical derivatives: dispatch
-// function.
-// Instantiation for the case QVectorLayout::byVDIM.
-template<>
-void TensorPhysDerivatives<QVectorLayout::byVDIM>(const int NE,
-                                                  const int vdim,
-                                                  const DofToQuad &maps,
-                                                  const GeometricFactors &geom,
-                                                  const Vector &e_vec,
-                                                  Vector &q_der)
-{
-   if (NE == 0) { return; }
-   const int dim = maps.FE->GetDim();
-   const int D1D = maps.ndof;
-   const int Q1D = maps.nqpt;
-
-   const int sdim = geom.mesh->SpaceDimension();
-
-   const real_t *B = maps.B.Read();
-   const real_t *G = maps.G.Read();
-   const real_t *J = geom.J.Read();
-   const real_t *X = e_vec.Read();
-   real_t *Y = q_der.Write();
-
-   constexpr QVectorLayout L = QVectorLayout::byVDIM;
-   constexpr bool P = true; // GRAD_PHYS
-
-   const int id = (vdim<<8) | (D1D<<4) | Q1D;
-
-   if (dim == 1)
-   {
-      return Derivatives1D<L,P>(NE,G,J,X,Y,sdim,vdim,D1D,Q1D);
-   }
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x134: return Derivatives2D<L,P,1,3,4,8>(NE,B,G,J,X,Y,sdim);
-         case 0x146: return Derivatives2D<L,P,1,4,6,4>(NE,B,G,J,X,Y,sdim);
-         case 0x158: return Derivatives2D<L,P,1,5,8,2>(NE,B,G,J,X,Y,sdim);
-
-         case 0x233: return Derivatives2D<L,P,2,3,3,8>(NE,B,G,J,X,Y,sdim);
-         case 0x234: return Derivatives2D<L,P,2,3,4,8>(NE,B,G,J,X,Y,sdim);
-         case 0x246: return Derivatives2D<L,P,2,4,6,4>(NE,B,G,J,X,Y,sdim);
-         case 0x258: return Derivatives2D<L,P,2,5,8,2>(NE,B,G,J,X,Y,sdim);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_D1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_Q1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Derivatives2D<L,P>(NE,B,G,J,X,Y,sdim,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x134: return Derivatives3D<L,P,1,3,4>(NE,B,G,J,X,Y);
-         case 0x146: return Derivatives3D<L,P,1,4,6>(NE,B,G,J,X,Y);
-         case 0x158: return Derivatives3D<L,P,1,5,8>(NE,B,G,J,X,Y);
-
-         case 0x334: return Derivatives3D<L,P,3,3,4>(NE,B,G,J,X,Y);
-         case 0x346: return Derivatives3D<L,P,3,4,6>(NE,B,G,J,X,Y);
-         case 0x358: return Derivatives3D<L,P,3,5,8>(NE,B,G,J,X,Y);
-         default:
-         {
-            const int MD = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            const int MQ = DeviceDofQuadLimits::Get().MAX_INTERP_1D;
-            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
-                        << " are not supported!");
-            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
-                        << MQ << " 1D points are not supported!");
-            Derivatives3D<L,P>(NE,B,G,J,X,Y,vdim,D1D,Q1D);
-            return;
-         }
-      }
-   }
-   mfem::out << "Unknown kernel 0x" << std::hex << id << std::endl;
-   MFEM_ABORT("Unknown kernel");
-}
-
-} // namespace quadrature_interpolator
-
-} // namespace internal
-
-} // namespace mfem

fem/quadinterpolator.cpp

@@ -10,7 +10,8 @@
 // CONTRIBUTING.md for details.
 
 #include "quadinterpolator.hpp"
-#include "qinterp/dispatch.hpp"
+#include "qinterp/grad.hpp"
+#include "qinterp/eval.hpp"
 #include "qspace.hpp"
 #include "../general/forall.hpp"
 #include "../linalg/dtensor.hpp"
@@ -19,6 +20,38 @@
 namespace mfem
 {
 
+namespace internal
+{
+namespace quadrature_interpolator
+{
+void InitEvalByNodesKernels();
+void InitEvalByVDimKernels();
+void InitEvalKernels();
+void InitDetKernels();
+template <bool P> void InitGradByNodesKernels();
+template <bool P> void InitGradByVDimKernels();
+}
+}
+
+QuadratureInterpolator::Kernels QuadratureInterpolator::kernels;
+QuadratureInterpolator::Kernels::Kernels()
+{
+   using namespace internal::quadrature_interpolator;
+
+   InitEvalByNodesKernels();
+   InitEvalByVDimKernels();
+   // Non-phys grad kernels
+   InitGradByNodesKernels<false>();
+   InitGradByVDimKernels<false>();
+   // Phys grad kernels
+   InitGradByNodesKernels<true>();
+   InitGradByVDimKernels<true>();
+   // Determinants
+   InitDetKernels();
+   // Non-tensor
+   InitEvalKernels();
+}
+
 QuadratureInterpolator::QuadratureInterpolator(const FiniteElementSpace &fes,
                                                const IntegrationRule &ir):
 
@@ -467,6 +500,7 @@ void QuadratureInterpolator::Mult(const Vector &e_vec,
    const int ne = fespace->GetNE();
    if (ne == 0) { return; }
    const int vdim = fespace->GetVDim();
+   const int sdim = fespace->GetMesh()->SpaceDimension();
    const FiniteElement *fe = fespace->GetFE(0);
    const bool use_tensor_eval =
       use_tensor_products &&
@@ -477,6 +511,8 @@ void QuadratureInterpolator::Mult(const Vector &e_vec,
       use_tensor_eval ? DofToQuad::TENSOR : DofToQuad::FULL;
    const DofToQuad &maps = fe->GetDofToQuad(*ir, mode);
    const int dim = maps.FE->GetDim();
+   const int nd = maps.ndof;
+   const int nq = maps.nqpt;
    const GeometricFactors *geom = nullptr;
    if (eval_flags & PHYSICAL_DERIVATIVES)
    {
@@ -492,202 +528,31 @@ void QuadratureInterpolator::Mult(const Vector &e_vec,
 
    if (use_tensor_eval)
    {
-      // TODO: use fused kernels
-      if (q_layout == QVectorLayout::byNODES)
+      if (eval_flags & VALUES)
       {
-         if (eval_flags & VALUES)
-         {
-            TensorValues<QVectorLayout::byNODES>(ne, vdim, maps, e_vec, q_val);
-         }
-         if (eval_flags & DERIVATIVES)
-         {
-            TensorDerivatives<QVectorLayout::byNODES>(
-               ne, vdim, maps, e_vec, q_der);
-         }
-         if (eval_flags & PHYSICAL_DERIVATIVES)
-         {
-            TensorPhysDerivatives<QVectorLayout::byNODES>(
-               ne, vdim, maps, *geom, e_vec, q_der);
-         }
+         TensorEvalKernels::Run(dim, q_layout, vdim, nd, nq, ne, maps.B.Read(),
+                                e_vec.Read(), q_val.Write(), vdim, nd, nq);
       }
-
-      if (q_layout == QVectorLayout::byVDIM)
+      if (eval_flags & (DERIVATIVES | PHYSICAL_DERIVATIVES))
       {
-         if (eval_flags & VALUES)
-         {
-            TensorValues<QVectorLayout::byVDIM>(ne, vdim, maps, e_vec, q_val);
-         }
-         if (eval_flags & DERIVATIVES)
-         {
-            TensorDerivatives<QVectorLayout::byVDIM>(
-               ne, vdim, maps, e_vec, q_der);
-         }
-         if (eval_flags & PHYSICAL_DERIVATIVES)
-         {
-            TensorPhysDerivatives<QVectorLayout::byVDIM>(
-               ne, vdim, maps, *geom, e_vec, q_der);
-         }
+         const bool phys = (eval_flags & PHYSICAL_DERIVATIVES);
+         const real_t *J = phys ? geom->J.Read() : nullptr;
+         const int s_dim = phys ? sdim : dim;
+         GradKernels::Run(dim, q_layout, phys, vdim, nd, nq, ne,
+                          maps.B.Read(), maps.G.Read(), J, e_vec.Read(),
+                          q_der.Write(), s_dim, vdim, nd, nq);
       }
       if (eval_flags & DETERMINANTS)
       {
-         TensorDeterminants(ne, vdim, maps, e_vec, q_det, d_buffer);
+         DetKernels::Run(dim, vdim, nd, nq, ne, maps.B.Read(),
+                         maps.G.Read(), e_vec.Read(), q_det.Write(), nd,
+                         nq, &d_buffer);
       }
    }
    else // use_tensor_eval == false
    {
-      const int nd = maps.ndof;
-      const int nq = maps.nqpt;
-
-      void (*mult)(const int NE,
-                   const int vdim,
-                   const QVectorLayout q_layout,
-                   const GeometricFactors *geom,
-                   const DofToQuad &maps,
-                   const Vector &e_vec,
-                   Vector &q_val,
-                   Vector &q_der,
-                   Vector &q_det,
-                   const int eval_flags) = NULL;
-
-      if (dim == 1)
-      {
-         mult = &Eval1D;
-      }
-      else if (vdim == 1) // dim == 2 || dim == 3
-      {
-         if (dim == 2)
-         {
-            switch (100*nd + nq)
-            {
-               // Q0
-               case 101: mult = &Eval2D<1,1,1>; break;
-               case 104: mult = &Eval2D<1,1,4>; break;
-               // Q1
-               case 404: mult = &Eval2D<1,4,4>; break;
-               case 409: mult = &Eval2D<1,4,9>; break;
-               // Q2
-               case 909: mult = &Eval2D<1,9,9>; break;
-               case 916: mult = &Eval2D<1,9,16>; break;
-               // Q3
-               case 1616: mult = &Eval2D<1,16,16>; break;
-               case 1625: mult = &Eval2D<1,16,25>; break;
-               case 1636: mult = &Eval2D<1,16,36>; break;
-               // Q4
-               case 2525: mult = &Eval2D<1,25,25>; break;
-               case 2536: mult = &Eval2D<1,25,36>; break;
-               case 2549: mult = &Eval2D<1,25,49>; break;
-               case 2564: mult = &Eval2D<1,25,64>; break;
-            }
-            if (nq >= 100 || !mult)
-            {
-               mult = &Eval2D<1,0,0>;
-            }
-         }
-         else if (dim == 3)
-         {
-            switch (1000*nd + nq)
-            {
-               // Q0
-               case 1001: mult = &Eval3D<1,1,1>; break;
-               case 1008: mult = &Eval3D<1,1,8>; break;
-               // Q1
-               case 8008: mult = &Eval3D<1,8,8>; break;
-               case 8027: mult = &Eval3D<1,8,27>; break;
-               // Q2
-               case 27027: mult = &Eval3D<1,27,27>; break;
-               case 27064: mult = &Eval3D<1,27,64>; break;
-               // Q3
-               case 64064: mult = &Eval3D<1,64,64>; break;
-               case 64125: mult = &Eval3D<1,64,125>; break;
-               case 64216: mult = &Eval3D<1,64,216>; break;
-               // Q4
-               case 125125: mult = &Eval3D<1,125,125>; break;
-               case 125216: mult = &Eval3D<1,125,216>; break;
-            }
-            if (nq >= 1000 || !mult)
-            {
-               mult = &Eval3D<1,0,0>;
-            }
-         }
-      }
-      else if (vdim == 3 && dim == 2)
-      {
-         switch (100*nd + nq)
-         {
-            // Q0
-            case 101: mult = &Eval2D<3,1,1>; break;
-            case 104: mult = &Eval2D<3,1,4>; break;
-            // Q1
-            case 404: mult = &Eval2D<3,4,4>; break;
-            case 409: mult = &Eval2D<3,4,9>; break;
-            // Q2
-            case 904: mult = &Eval2D<3,9,4>; break;
-            case 909: mult = &Eval2D<3,9,9>; break;
-            case 916: mult = &Eval2D<3,9,16>; break;
-            case 925: mult = &Eval2D<3,9,25>; break;
-            // Q3
-            case 1616: mult = &Eval2D<3,16,16>; break;
-            case 1625: mult = &Eval2D<3,16,25>; break;
-            case 1636: mult = &Eval2D<3,16,36>; break;
-            // Q4
-            case 2525: mult = &Eval2D<3,25,25>; break;
-            case 2536: mult = &Eval2D<3,25,36>; break;
-            case 2549: mult = &Eval2D<3,25,49>; break;
-            case 2564: mult = &Eval2D<3,25,64>; break;
-            default:   mult = &Eval2D<3,0,0>;
-         }
-      }
-      else if (vdim == dim)
-      {
-         if (dim == 2)
-         {
-            switch (100*nd + nq)
-            {
-               // Q1
-               case 404: mult = &Eval2D<2,4,4>; break;
-               case 409: mult = &Eval2D<2,4,9>; break;
-               // Q2
-               case 909: mult = &Eval2D<2,9,9>; break;
-               case 916: mult = &Eval2D<2,9,16>; break;
-               // Q3
-               case 1616: mult = &Eval2D<2,16,16>; break;
-               case 1625: mult = &Eval2D<2,16,25>; break;
-               case 1636: mult = &Eval2D<2,16,36>; break;
-               // Q4
-               case 2525: mult = &Eval2D<2,25,25>; break;
-               case 2536: mult = &Eval2D<2,25,36>; break;
-               case 2549: mult = &Eval2D<2,25,49>; break;
-               case 2564: mult = &Eval2D<2,25,64>; break;
-            }
-            if (nq >= 100 || !mult) { mult = &Eval2D<2,0,0>; }
-         }
-         else if (dim == 3)
-         {
-            switch (1000*nd + nq)
-            {
-               // Q1
-               case 8008: mult = &Eval3D<3,8,8>; break;
-               case 8027: mult = &Eval3D<3,8,27>; break;
-               // Q2
-               case 27027: mult = &Eval3D<3,27,27>; break;
-               case 27064: mult = &Eval3D<3,27,64>; break;
-               case 27125: mult = &Eval3D<3,27,125>; break;
-               // Q3
-               case 64064: mult = &Eval3D<3,64,64>; break;
-               case 64125: mult = &Eval3D<3,64,125>; break;
-               case 64216: mult = &Eval3D<3,64,216>; break;
-               // Q4
-               case 125125: mult = &Eval3D<3,125,125>; break;
-               case 125216: mult = &Eval3D<3,125,216>; break;
-            }
-            if (nq >= 1000 || !mult) {  mult = &Eval3D<3,0,0>; }
-         }
-      }
-      if (mult)
-      {
-         mult(ne,vdim,q_layout,geom,maps,e_vec,q_val,q_der,q_det,eval_flags);
-      }
-      else { MFEM_ABORT("case not supported yet"); }
+      EvalKernels::Run(dim, vdim, maps.ndof, maps.nqpt, ne,vdim,q_layout,
+                       geom, maps,e_vec, q_val,q_der,q_det,eval_flags);
    }
 }
 
@@ -731,4 +596,186 @@ void QuadratureInterpolator::Determinants(const Vector &e_vec,
    Mult(e_vec, DETERMINANTS, empty, empty, q_det);
 }
 
+/// @cond Suppress_Doxygen_warnings
+
+namespace
+{
+using EvalKernel = QuadratureInterpolator::EvalKernelType;
+using TensorEvalKernel = QuadratureInterpolator::TensorEvalKernelType;
+using GradKernel = QuadratureInterpolator::GradKernelType;
+
+template <QVectorLayout Q_LAYOUT>
+TensorEvalKernel FallbackTensorEvalKernel(int DIM)
+{
+   if (DIM == 1) { return internal::quadrature_interpolator::Values1D<Q_LAYOUT>; }
+   else if (DIM == 2) { return internal::quadrature_interpolator::Values2D<Q_LAYOUT>; }
+   else if (DIM == 3) { return internal::quadrature_interpolator::Values3D<Q_LAYOUT>; }
+   else { MFEM_ABORT(""); }
+}
+
+template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS>
+GradKernel GetGradKernel(int DIM)
+{
+   if (DIM == 1) { return internal::quadrature_interpolator::Derivatives1D<Q_LAYOUT, GRAD_PHYS>; }
+   else if (DIM == 2) { return internal::quadrature_interpolator::Derivatives2D<Q_LAYOUT, GRAD_PHYS>; }
+   else if (DIM == 3) { return internal::quadrature_interpolator::Derivatives3D<Q_LAYOUT, GRAD_PHYS>; }
+   else { MFEM_ABORT(""); }
+}
+
+template<QVectorLayout Q_LAYOUT>
+GradKernel GetGradKernel(int DIM, bool GRAD_PHYS)
+{
+   if (GRAD_PHYS) { return GetGradKernel<Q_LAYOUT, true>(DIM); }
+   else { return GetGradKernel<Q_LAYOUT, false>(DIM); }
+}
+} // namespace
+
+template <int DIM, int VDIM, int ND, int NQ>
+EvalKernel QuadratureInterpolator::EvalKernels::Kernel()
+{
+   using namespace internal::quadrature_interpolator;
+   if (DIM == 1) { return Eval1D; }
+   else if (DIM == 2) { return Eval2D<VDIM,ND,NQ>; }
+   else if (DIM == 3) { return Eval3D<VDIM,ND,NQ>; }
+   else { MFEM_ABORT(""); }
+}
+
+template <int DIM>
+EvalKernel GetEvalKernelVDimFallback(int VDIM)
+{
+   using EvalKernels = QuadratureInterpolator::EvalKernels;
+   if (VDIM == 1) { return EvalKernels::Kernel<DIM,1,0,0>(); }
+   else if (VDIM == 2) { return EvalKernels::Kernel<DIM,2,0,0>(); }
+   else if (VDIM == 3) { return EvalKernels::Kernel<DIM,3,0,0>(); }
+   else { MFEM_ABORT(""); }
+}
+
+EvalKernel QuadratureInterpolator::EvalKernels::Fallback(
+   int DIM, int VDIM, int ND, int NQ)
+{
+   if (DIM == 1) { return GetEvalKernelVDimFallback<1>(VDIM); }
+   else if (DIM == 2) { return GetEvalKernelVDimFallback<2>(VDIM); }
+   else if (DIM == 3) { return GetEvalKernelVDimFallback<3>(VDIM); }
+   else { MFEM_ABORT(""); }
+}
+
+TensorEvalKernel QuadratureInterpolator::TensorEvalKernels::Fallback(
+   int DIM, QVectorLayout Q_LAYOUT, int, int, int)
+{
+   if (Q_LAYOUT == QVectorLayout::byNODES) { return FallbackTensorEvalKernel<QVectorLayout::byNODES>(DIM); }
+   else { return FallbackTensorEvalKernel<QVectorLayout::byVDIM>(DIM); }
+}
+
+GradKernel QuadratureInterpolator::GradKernels::Fallback(
+   int DIM, QVectorLayout Q_LAYOUT, bool GRAD_PHYS, int, int, int)
+{
+   if (Q_LAYOUT == QVectorLayout::byNODES) { return GetGradKernel<QVectorLayout::byNODES>(DIM, GRAD_PHYS); }
+   else { return GetGradKernel<QVectorLayout::byVDIM>(DIM, GRAD_PHYS); }
+}
+
+/// @endcond
+
+namespace internal
+{
+namespace quadrature_interpolator
+{
+void InitEvalKernels()
+{
+   using k = QuadratureInterpolator::EvalKernels;
+   // 2D, VDIM = 1
+   k::Specialization<2,1,1,1>::Add();
+   k::Specialization<2,1,1,4>::Add();
+   // Q1
+   k::Specialization<2,1,4,4>::Add();
+   k::Specialization<2,1,4,9>::Add();
+   // Q2
+   k::Specialization<2,1,9,9>::Add();
+   k::Specialization<2,1,9,16>::Add();
+   // Q3
+   k::Specialization<2,1,16,16>::Add();
+   k::Specialization<2,1,16,25>::Add();
+   k::Specialization<2,1,16,36>::Add();
+   // Q4
+   k::Specialization<2,1,25,25>::Add();
+   k::Specialization<2,1,25,36>::Add();
+   k::Specialization<2,1,25,49>::Add();
+   k::Specialization<2,1,25,64>::Add();
+
+   // 3D, VDIM = 1
+   // Q0
+   k::Specialization<3,1,1,1>::Add();
+   k::Specialization<3,1,1,8>::Add();
+   // Q1
+   k::Specialization<3,1,8,8>::Add();
+   k::Specialization<3,1,8,27>::Add();
+   // Q2
+   k::Specialization<3,1,27,27>::Add();
+   k::Specialization<3,1,27,64>::Add();
+   // Q3
+   k::Specialization<3,1,64,64>::Add();
+   k::Specialization<3,1,64,125>::Add();
+   k::Specialization<3,1,64,216>::Add();
+   // Q4
+   k::Specialization<3,1,125,125>::Add();
+   k::Specialization<3,1,125,216>::Add();
+
+   // 2D, VDIM = 3
+   // Q0
+   k::Specialization<2,3,1,1>::Add();
+   k::Specialization<2,3,1,4>::Add();
+   // Q1
+   k::Specialization<2,3,4,4>::Add();
+   k::Specialization<2,3,4,9>::Add();
+   // Q2
+   k::Specialization<2,3,9,4>::Add();
+   k::Specialization<2,3,9,9>::Add();
+   k::Specialization<2,3,9,16>::Add();
+   k::Specialization<2,3,9,25>::Add();
+   // Q3
+   k::Specialization<2,3,16,16>::Add();
+   k::Specialization<2,3,16,25>::Add();
+   k::Specialization<2,3,16,36>::Add();
+   // Q4
+   k::Specialization<2,3,25,25>::Add();
+   k::Specialization<2,3,25,36>::Add();
+   k::Specialization<2,3,25,49>::Add();
+   k::Specialization<2,3,25,64>::Add();
+
+   // 2D, VDIM = 2
+   // Q1
+   k::Specialization<2,2,4,4>::Add();
+   k::Specialization<2,2,4,9>::Add();
+   // Q2
+   k::Specialization<2,2,9,9>::Add();
+   k::Specialization<2,2,9,16>::Add();
+   // Q3
+   k::Specialization<2,2,16,16>::Add();
+   k::Specialization<2,2,16,25>::Add();
+   k::Specialization<2,2,16,36>::Add();
+   // Q4
+   k::Specialization<2,2,25,25>::Add();
+   k::Specialization<2,2,25,36>::Add();
+   k::Specialization<2,2,25,49>::Add();
+   k::Specialization<2,2,25,64>::Add();
+
+   // 3D, VDIM = 3
+   // Q1
+   k::Specialization<3,3,8,8>::Add();
+   k::Specialization<3,3,8,27>::Add();
+   // Q2
+   k::Specialization<3,3,27,27>::Add();
+   k::Specialization<3,3,27,64>::Add();
+   k::Specialization<3,3,27,125>::Add();
+   // Q3
+   k::Specialization<3,3,64,64>::Add();
+   k::Specialization<3,3,64,125>::Add();
+   k::Specialization<3,3,64,216>::Add();
+   // Q4
+   k::Specialization<3,3,125,125>::Add();
+   k::Specialization<3,3,125,216>::Add();
+}
+
+} // namespace quadrature_Interpolator
+} // namespace internal
+
 } // namespace mfem

fem/quadinterpolator.hpp

@@ -13,6 +13,7 @@
 #define MFEM_QUADINTERP
 
 #include "fespace.hpp"
+#include "kernel_dispatch.hpp"
 
 namespace mfem
 {
@@ -130,6 +131,29 @@ public:
    /// Perform the transpose operation of Mult(). (TODO)
    void MultTranspose(unsigned eval_flags, const Vector &q_val,
                       const Vector &q_der, Vector &e_vec) const;
+
+
+   using TensorEvalKernelType = void(*)(const int, const real_t *, const real_t *,
+                                        real_t *, const int, const int, const int);
+   using GradKernelType = void(*)(const int, const real_t *, const real_t *,
+                                  const real_t *, const real_t *, real_t *,
+                                  const int, const int, const int, const int);
+   using DetKernelType = void(*)(const int NE, const real_t *, const real_t *,
+                                 const real_t *, real_t *, const int, const int,
+                                 Vector *);
+   using EvalKernelType = void(*)(const int, const int, const QVectorLayout,
+                                  const GeometricFactors *, const DofToQuad &,
+                                  const Vector &, Vector &, Vector &, Vector &,
+                                  const int);
+
+   MFEM_REGISTER_KERNELS(TensorEvalKernels, TensorEvalKernelType,
+                         (int, QVectorLayout, int, int, int), (int));
+   MFEM_REGISTER_KERNELS(GradKernels, GradKernelType,
+                         (int, QVectorLayout, bool, int, int, int), (int));
+   MFEM_REGISTER_KERNELS(DetKernels, DetKernelType, (int, int, int, int));
+   MFEM_REGISTER_KERNELS(EvalKernels, EvalKernelType, (int, int, int, int));
+
+   static struct Kernels { Kernels(); } kernels;
 };
 
 }

miniapps/hdiv-linear-solver/change_basis.cpp

@@ -10,7 +10,6 @@
 // CONTRIBUTING.md for details.
 
 #include "change_basis.hpp"
-#include "../../fem/qinterp/dispatch.hpp"
 #include "../../general/forall.hpp"
 #include "../../linalg/dtensor.hpp"
 
@@ -105,17 +104,25 @@ ChangeOfBasis_L2::ChangeOfBasis_L2(FiniteElementSpace &fes)
 void ChangeOfBasis_L2::Mult(const Vector &x, Vector &y) const
 {
    if (no_op) { y = x; return; }
-   using namespace internal::quadrature_interpolator;
    dof2quad.B.MakeRef(B_1d);
-   TensorValues<QVectorLayout::byVDIM>(ne, 1, dof2quad, x, y);
+   const int dim = dof2quad.FE->GetDim();
+   const int nd = dof2quad.ndof;
+   const int nq = dof2quad.nqpt;
+   QuadratureInterpolator::TensorEvalKernels::Run(
+      dim, QVectorLayout::byVDIM, 1, nd, nq, ne, dof2quad.B.Read(), x.Read(),
+      y.Write(), 1, nd, nq);
 }
 
 void ChangeOfBasis_L2::MultTranspose(const Vector &x, Vector &y) const
 {
    if (no_op) { y = x; return; }
-   using namespace internal::quadrature_interpolator;
    dof2quad.B.MakeRef(Bt_1d);
-   TensorValues<QVectorLayout::byVDIM>(ne, 1, dof2quad, x, y);
+   const int dim = dof2quad.FE->GetDim();
+   const int nd = dof2quad.ndof;
+   const int nq = dof2quad.nqpt;
+   QuadratureInterpolator::TensorEvalKernels::Run(
+      dim, QVectorLayout::byVDIM, 1, nd, nq, ne, dof2quad.B.Read(), x.Read(),
+      y.Write(), 1, nd, nq);
 }
 
 ChangeOfBasis_RT::ChangeOfBasis_RT(FiniteElementSpace &fes)