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
 *~
 

CHANGELOG

@@ -11,7 +11,19 @@
 Version 4.7.1 (development)
 ===========================
 
-- Added NURBS-based H(div) and H(curl) elements in 2D and 3D. Only on single 
+- Refactored the `ARKStepSolver` class (ARKODE interface) to use
+  `TimeDependentOperator::Mult` only when the associated ODE operator is
+  expressed in explicit form (i.e., `TimeDependentOperator::isExplicit()`),
+  otherwise `TimeDependentOperator::ExplicitMult` is used. A check has been
+  added to `ARKStepSolver` to verify that the associated ODE operator is not in
+  explicit form when a mass matrix solver is enabled via a call to either the
+  `UseMFEMMassLinearSolver` or `UseSundialsMassLinearSolver` methods. This is
+  because enabling a mass matrix solver assumes that F(u,k,t) = M k in the
+  associated ODE operator.
+
+- Added support for custom interpolation procedure in FindPointsGSLIB.
+
+- Added NURBS-based H(div) and H(curl) elements in 2D and 3D. Only on single
   patch meshes. Only implemented for serial computations.
 
 - Added miniapps to demonstrate the H(div) and H(curl) NURBS elements.
@@ -45,6 +57,9 @@ Meshing improvements
 
 - Added support for internal boundary elements in nonconforming meshes.
 
+- Added ExodusII output capability. The writer can handle first-order (Pyramid5,
+  Wedge6, Hex8, Tet4) and second-order FE types (Pyramid14, Wedge18, Hex27, Tet10).
+
 - The ReadCubit Genesis mesh importer has been rewritten to improve readability.
 
 Discretization improvements

examples/hiop/ex9p.cpp

@@ -96,6 +96,7 @@ public:
    {
       Vector w_glob(width);
       pfes.Dof_TrueDof_Matrix()->MultTranspose(w, w_glob);
+      w_glob.HostReadWrite(); // read+write -> can use w_glob(i) (non-const)
       for (int i = 0; i < width; i++) { grad(0, i) = w_glob(i); }
    }
 

examples/sundials/CMakeLists.txt

@@ -31,11 +31,21 @@ include_directories(BEFORE ${PROJECT_BINARY_DIR})
 add_custom_target(test_sundials
   ${CMAKE_CTEST_COMMAND} -R sundials USES_TERMINAL)
 
-# Add one executable per cpp file, adding "sundials_" as prefix. Sets
-# "test_sundials" as a target that depends on the given examples.
+# Add one executable per cpp file, adding "sundials_" as prefix so the CMake
+# target is unique from those in the non-SUNDIALS examples. Also sets
+# "test_sundials" as a target that depends on the given SUNDIALS examples.
 set(PFX sundials_)
 add_mfem_examples(SUNDIALS_EXAMPLES_SRCS ${PFX} "" test_sundials)
 
+# Remove "sundials_" prefix from exectuable name for consistency with GNU build
+# system.
+foreach(SRC_FILE ${SUNDIALS_EXAMPLES_SRCS})
+  get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
+  string(REPLACE ".cpp" "" TARGET_NAME "${PFX}${SRC_FILENAME}")
+  string(REPLACE ${PFX} "" EXE_NAME ${TARGET_NAME})
+  set_target_properties(${TARGET_NAME} PROPERTIES OUTPUT_NAME ${EXE_NAME})
+endforeach()
+
 # Testing.
 # The SUNDIALS tests can be run separately using the target "test_sundials"
 # which builds the examples and runs:
@@ -51,7 +61,10 @@ if (MFEM_ENABLE_TESTING)
   set(EX10_COMMON_OPTS -m ../../data/beam-quad.mesh -o 2 -s 5 -dt 0.15 -tf 6 -vs 10)
   set(EX10_TEST_OPTS ${EX10_COMMON_OPTS}  -r  2)
   set(EX10P_TEST_OPTS ${EX10_COMMON_OPTS} -rp 1)
-  # Example 16: use the default options
+  # Example 16: test ARKODE with implicit time stepping using mass form
+  set(EX16_COMMON_OPTS -s 15)
+  set(EX16_TEST_OPTS ${EX16_COMMON_OPTS})
+  set(EX16P_TEST_OPTS ${EX16_COMMON_OPTS})
 
   # Add the tests: one test per source file.
   foreach(SRC_FILE ${SUNDIALS_EXAMPLES_SRCS})

examples/sundials/ex10.cpp

@@ -1,7 +1,9 @@
 //                                MFEM Example 10
 //                             SUNDIALS Modification
 //
-// Compile with: make ex10
+// Compile with:
+//    make ex10             (GNU make)
+//    make sundials_ex10    (CMake)
 //
 // Sample runs:
 //    ex10 -m ../../data/beam-quad.mesh -r 2 -o 2 -s 12 -dt 0.15 -vs 10

examples/sundials/ex10p.cpp

@@ -1,7 +1,9 @@
 //                       MFEM Example 10 - Parallel Version
 //                             SUNDIALS Modification
 //
-// Compile with: make ex10p
+// Compile with:
+//    make ex10p            (GNU make)
+//    make sundials_ex10p   (CMake)
 //
 // Sample runs:
 //    mpirun -np 4 ex10p -m ../../data/beam-quad.mesh -rp 1 -o 2 -s 12 -dt 0.15 -vs 10

examples/sundials/ex16.cpp

@@ -1,15 +1,21 @@
 //                                MFEM Example 16
 //                             SUNDIALS Modification
 //
-// Compile with: make ex16
+// Compile with:
+//    make ex16             (GNU make)
+//    make sundials_ex16    (CMake)
 //
 // Sample runs:  ex16
 //               ex16 -m ../../data/inline-tri.mesh
 //               ex16 -m ../../data/disc-nurbs.mesh -tf 2
 //               ex16 -s 12 -a 0.0 -k 1.0
+//               ex16 -s 15 -a 0.0 -k 1.0
 //               ex16 -s 8 -a 1.0 -k 0.0 -dt 1e-4 -tf 5e-2 -vs 25
+//               ex16 -s 11 -a 1.0 -k 0.0 -dt 1e-4 -tf 5e-2 -vs 25
 //               ex16 -s 9 -a 0.5 -k 0.5 -o 4 -dt 1e-4 -tf 2e-2 -vs 25
+//               ex16 -s 12 -a 0.5 -k 0.5 -o 4 -dt 1e-4 -tf 2e-2 -vs 25
 //               ex16 -s 10 -dt 1.0e-4 -tf 4.0e-2 -vs 40
+//               ex16 -s 13 -dt 1.0e-4 -tf 4.0e-2 -vs 40
 //               ex16 -m ../../data/fichera-q2.mesh
 //               ex16 -m ../../data/escher.mesh
 //               ex16 -m ../../data/beam-tet.mesh -tf 10 -dt 0.1
@@ -37,75 +43,102 @@
 using namespace std;
 using namespace mfem;
 
-/** After spatial discretization, the conduction model can be written as:
+/** After spatial discretization, the conduction model is expressed as
  *
- *     du/dt = M^{-1}(-Ku)
+ *   M du/dt = - K(u) u
  *
  *  where u is the vector representing the temperature, M is the mass matrix,
- *  and K is the diffusion operator with diffusivity depending on u:
+ *  and K(u) is the diffusion operator with diffusivity depending on u:
  *  (\kappa + \alpha u).
  *
- *  Class ConductionOperator represents the right-hand side of the above ODE.
+ *  Class ConductionOperatorOperator represents the above ODE operator in the
+ *  general form F(u, k, t) = G(u, t) where
+ *
+ *    1. F(u, du/dt, t) = du/dt            (ODE is expressed in EXPLICIT form)
+ *       G(u, t) = - inv(M) K(u) u
+ *    2. F(u, du/dt, t) = M du/dt          (ODE is expressed in IMPLICIT form)
+ *       G(u, t) = - K(u) u
  */
 class ConductionOperator : public TimeDependentOperator
 {
-protected:
    FiniteElementSpace &fespace;
    Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
 
-   BilinearForm *M;
-   BilinearForm *K;
+   BilinearForm M;
+   SparseMatrix Mmat;
 
-   SparseMatrix Mmat, Kmat;
-   SparseMatrix *T; // T = M + dt K
+   const real_t alpha, kappa;
+   std::unique_ptr<BilinearForm> K;
+   SparseMatrix Kmat;
+
+   std::unique_ptr<SparseMatrix> T; // T = M + gam K(u)
 
    CGSolver M_solver; // Krylov solver for inverting the mass matrix M
    DSmoother M_prec;  // Preconditioner for the mass matrix M
 
-   CGSolver T_solver; // Implicit solver for T = M + dt K
+   CGSolver T_solver; // Implicit solver for T = M + gam K(u)
    DSmoother T_prec;  // Preconditioner for the implicit solver
 
-   double alpha, kappa;
-
    mutable Vector z; // auxiliary vector
 
 public:
-   ConductionOperator(FiniteElementSpace &f, double alpha, double kappa,
-                      const Vector &u);
 
-   virtual void Mult(const Vector &u, Vector &du_dt) const;
+   ConductionOperator(FiniteElementSpace &f, const real_t alpha,
+                      const real_t kappa, const Vector &u,
+                      const Type &ode_expression_type);
 
-   /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
-       This is the only requirement for high-order SDIRK implicit integration.*/
-   virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
+   // Compute K(u_n) for use as an approximation in - K(u) u
+   void SetConductionTensor(const Vector &u);
 
-   /// Custom Jacobian system solver for the SUNDIALS time integrators.
-   /** For the ODE system represented by ConductionOperator
+   /** Compute G(u, t) as defined in the IMPLICIT expression form of the ODE
+    operator, i.e., @a v = - K(u_n) @a u. Note that K(u_n) is an
+      approximation to K(u). */
+   void ExplicitMult(const Vector &u, Vector &v) const override;
 
-           M du/dt = -K(u),
+   /** Solve for k in F(u, k, t) = G(u, t) for either EXPLICIT or IMPLICIT
+       expression forms of the ODE operator, i.e., @a k = - inv(M) K(u_n) @a u.
+       Note that K(u_n) is an approximation to K(u). */
+   void Mult(const Vector &u, Vector &k) const override;
 
-       this class facilitates the solution of linear systems of the form
+   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t) for either EXPLICIT
+       or IMPLICIT expression forms of the ODE operator, i.e.,
+       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u . Note that K(u_n) is an
+       approximation to K(u). */
+   void ImplicitSolve(const real_t gam, const Vector &u, Vector &k) override;
 
-           (M + γK) y = M b,
+   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F for
+       either EXPLICIT or IMPLICIT expression forms of the ODE operator, i.e.,
+       [M - @a gam Jf(u)] dk = G - F, where Jf(u) is an approximation of the
+       Jacobian of -K(u) u. The approximation chosen here is Jf(u) = -K(u_n). */
+   int SUNImplicitSetup(const Vector &u, const Vector &fu, int jok, int *jcur,
+                        real_t gam) override;
 
-       for given b, u (not used), and γ = GetTimeStep(). */
+   /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
+       with the residual @a r providing either
+        1. @a r = G - F = inv(M) f(u) - k       (EXPLICIT expression form)
+        1. @a r = G - F = f(u) - M k            (IMPLICIT expression form)
+       */
+   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
 
-   /** Setup the system (M + dt K) x = M b. This method is used by the implicit
-       SUNDIALS solvers. */
-   virtual int SUNImplicitSetup(const Vector &x, const Vector &fx,
-                                int jok, int *jcur, double gamma);
+   int SUNMassSetup() override;
 
-   /** Solve the system (M + dt K) x = M b. This method is used by the implicit
-       SUNDIALS solvers. */
-   virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol);
+   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override;
 
-   /// Update the diffusion BilinearForm K using the given true-dof vector `u`.
-   void SetParameters(const Vector &u);
-
-   virtual ~ConductionOperator();
+   int SUNMassMult(const Vector &x, Vector &v) override;
 };
 
-double InitialTemperature(const Vector &x);
+real_t InitialTemperature(const Vector &x)
+{
+   if (x.Norml2() < 0.5)
+   {
+      return 2.0;
+   }
+   else
+   {
+      return 1.0;
+   }
+}
+
 
 int main(int argc, char *argv[])
 {
@@ -117,16 +150,16 @@ int main(int argc, char *argv[])
    int ref_levels = 2;
    int order = 2;
    int ode_solver_type = 9; // CVODE implicit BDF
-   double t_final = 0.5;
-   double dt = 1.0e-2;
-   double alpha = 1.0e-2;
-   double kappa = 0.5;
+   real_t t_final = 0.5;
+   real_t dt = 1.0e-2;
+   real_t alpha = 1.0e-2;
+   real_t kappa = 0.5;
    bool visualization = true;
    bool visit = false;
    int vis_steps = 5;
 
    // Relative and absolute tolerances for CVODE and ARKODE.
-   const double reltol = 1e-4, abstol = 1e-4;
+   const real_t reltol = 1e-4, abstol = 1e-4;
 
    int precision = 8;
    cout.precision(precision);
@@ -151,7 +184,10 @@ int main(int argc, char *argv[])
                   "9  - CVODE (implicit BDF),\n\t"
                   "10 - ARKODE (default explicit),\n\t"
                   "11 - ARKODE (explicit Fehlberg-6-4-5),\n\t"
-                  "12 - ARKODE (default impicit).");
+                  "12 - ARKODE (default implicit),\n\t"
+                  "13 - ARKODE (default explicit with MFEM mass solve),\n\t"
+                  "14 - ARKODE (explicit Fehlberg-6-4-5 with MFEM mass solve),\n\t"
+                  "15 - ARKODE (default implicit with MFEM mass solve).");
    args.AddOption(&t_final, "-tf", "--t-final",
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
@@ -174,16 +210,13 @@ int main(int argc, char *argv[])
       args.PrintUsage(cout);
       return 1;
    }
-   if (ode_solver_type < 1 || ode_solver_type > 12)
-   {
-      cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
-      return 3;
-   }
    args.PrintOptions(cout);
 
+   bool use_mass_solver = ode_solver_type >= 13;
+
    // 2. Read the mesh from the given mesh file. We can handle triangular,
    //    quadrilateral, tetrahedral and hexahedral meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   std::unique_ptr<Mesh> mesh(new Mesh(mesh_file, 1, 1));
    int dim = mesh->Dimension();
 
    // 3. Refine the mesh to increase the resolution. In this example we do
@@ -197,7 +230,7 @@ int main(int argc, char *argv[])
    // 4. Define the vector finite element space representing the current and the
    //    initial temperature, u_ref.
    H1_FECollection fe_coll(order, dim);
-   FiniteElementSpace fespace(mesh, &fe_coll);
+   FiniteElementSpace fespace(mesh.get(), &fe_coll);
 
    int fe_size = fespace.GetTrueVSize();
    cout << "Number of temperature unknowns: " << fe_size << endl;
@@ -211,8 +244,17 @@ int main(int argc, char *argv[])
    Vector u;
    u_gf.GetTrueDofs(u);
 
-   // 6. Initialize the conduction operator and the visualization.
-   ConductionOperator oper(fespace, alpha, kappa, u);
+   // 6. Initialize the conduction ODE operator and the visualization.
+   ConductionOperator::Type ode_expression_type;
+   if (use_mass_solver)
+   {
+      ode_expression_type = ConductionOperator::Type::IMPLICIT;
+   }
+   else
+   {
+      ode_expression_type = ConductionOperator::Type::EXPLICIT;
+   }
+   ConductionOperator oper(fespace, alpha, kappa, u, ode_expression_type);
 
    u_gf.SetFromTrueDofs(u);
    {
@@ -224,7 +266,7 @@ int main(int argc, char *argv[])
       u_gf.Save(osol);
    }
 
-   VisItDataCollection visit_dc("Example16", mesh);
+   VisItDataCollection visit_dc("Example16", mesh.get());
    visit_dc.RegisterField("temperature", &u_gf);
    if (visit)
    {
@@ -258,52 +300,75 @@ int main(int argc, char *argv[])
    }
 
    // 7. Define the ODE solver used for time integration.
-   double t = 0.0;
-   ODESolver *ode_solver = NULL;
-   CVODESolver *cvode = NULL;
-   ARKStepSolver *arkode = NULL;
+   real_t t = 0.0;
+   std::unique_ptr<ODESolver> ode_solver;
    switch (ode_solver_type)
    {
       // MFEM explicit methods
-      case 1: ode_solver = new ForwardEulerSolver; break;
-      case 2: ode_solver = new RK2Solver(0.5); break; // midpoint method
-      case 3: ode_solver = new RK3SSPSolver; break;
-      case 4: ode_solver = new RK4Solver; break;
+      case 1: ode_solver = std::make_unique<ForwardEulerSolver>(); break;
+      case 2: ode_solver = std::make_unique<RK2Solver>(0.5); break; // midpoint method
+      case 3: ode_solver = std::make_unique<RK3SSPSolver>(); break;
+      case 4: ode_solver = std::make_unique<RK4Solver>(); break;
       // MFEM implicit L-stable methods
-      case 5: ode_solver = new BackwardEulerSolver; break;
-      case 6: ode_solver = new SDIRK23Solver(2); break;
-      case 7: ode_solver = new SDIRK33Solver; break;
+      case 5: ode_solver = std::make_unique<BackwardEulerSolver>(); break;
+      case 6: ode_solver = std::make_unique<SDIRK23Solver>(2); break;
+      case 7: ode_solver = std::make_unique<SDIRK33Solver>(); break;
       // CVODE
       case 8:
-         cvode = new CVODESolver(CV_ADAMS);
-         cvode->Init(oper);
-         cvode->SetSStolerances(reltol, abstol);
-         cvode->SetMaxStep(dt);
-         ode_solver = cvode; break;
       case 9:
-         cvode = new CVODESolver(CV_BDF);
+      {
+         int cvode_solver_type;
+         if (ode_solver_type == 8)
+         {
+            cvode_solver_type = CV_ADAMS;
+         }
+         else
+         {
+            cvode_solver_type = CV_BDF;
+         }
+         std::unique_ptr<CVODESolver> cvode(new CVODESolver(cvode_solver_type));
          cvode->Init(oper);
          cvode->SetSStolerances(reltol, abstol);
          cvode->SetMaxStep(dt);
-         ode_solver = cvode; break;
+         ode_solver = std::move(cvode);
+         break;
+      }
       // ARKODE
       case 10:
       case 11:
-         arkode = new ARKStepSolver(ARKStepSolver::EXPLICIT);
+      case 12:
+      case 13:
+      case 14:
+      case 15:
+      {
+         ARKStepSolver::Type arkode_solver_type;
+         if (ode_solver_type == 12 || ode_solver_type == 15)
+         {
+            arkode_solver_type = ARKStepSolver::IMPLICIT;
+         }
+         else
+         {
+            arkode_solver_type = ARKStepSolver::EXPLICIT;
+         }
+         std::unique_ptr<ARKStepSolver> arkode(
+            new ARKStepSolver(arkode_solver_type));
          arkode->Init(oper);
          arkode->SetSStolerances(reltol, abstol);
          arkode->SetMaxStep(dt);
-         if (ode_solver_type == 11)
+         if (ode_solver_type == 11 || ode_solver_type == 14)
          {
             arkode->SetERKTableNum(ARKODE_FEHLBERG_13_7_8);
          }
-         ode_solver = arkode; break;
-      case 12:
-         arkode = new ARKStepSolver(ARKStepSolver::IMPLICIT);
-         arkode->Init(oper);
-         arkode->SetSStolerances(reltol, abstol);
-         arkode->SetMaxStep(dt);
-         ode_solver = arkode; break;
+         if (use_mass_solver)
+         {
+            arkode->UseMFEMMassLinearSolver(SUNFALSE);
+         }
+         ode_solver = std::move(arkode);
+         break;
+      }
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
    }
 
    // Initialize MFEM integrators, SUNDIALS integrators are initialized above
@@ -311,8 +376,14 @@ int main(int argc, char *argv[])
 
    // Since we want to update the diffusion coefficient after every time step,
    // we need to use the "one-step" mode of the SUNDIALS solvers.
-   if (cvode) { cvode->SetStepMode(CV_ONE_STEP); }
-   if (arkode) { arkode->SetStepMode(ARK_ONE_STEP); }
+   if (CVODESolver* cvode = dynamic_cast<CVODESolver*>(ode_solver.get()))
+   {
+      cvode->SetStepMode(CV_ONE_STEP);
+   }
+   else if (ARKStepSolver* arkode = dynamic_cast<ARKStepSolver*>(ode_solver.get()))
+   {
+      arkode->SetStepMode(ARK_ONE_STEP);
+   }
 
    // 8. Perform time-integration (looping over the time iterations, ti, with a
    //    time-step dt).
@@ -323,7 +394,7 @@ int main(int argc, char *argv[])
    bool last_step = false;
    for (int ti = 1; !last_step; ti++)
    {
-      double dt_real = min(dt, t_final - t);
+      real_t dt_real = min(dt, t_final - t);
 
       // Note that since we are using the "one-step" mode of the SUNDIALS
       // solvers, they will, generally, step over the final time and will not
@@ -337,8 +408,14 @@ int main(int argc, char *argv[])
       if (last_step || (ti % vis_steps) == 0)
       {
          cout << "step " << ti << ", t = " << t << endl;
-         if (cvode) { cvode->PrintInfo(); }
-         if (arkode) { arkode->PrintInfo(); }
+         if (CVODESolver* cvode = dynamic_cast<CVODESolver*>(ode_solver.get()))
+         {
+            cvode->PrintInfo();
+         }
+         else if (ARKStepSolver* arkode = dynamic_cast<ARKStepSolver*>(ode_solver.get()))
+         {
+            arkode->PrintInfo();
+         }
 
          u_gf.SetFromTrueDofs(u);
          if (visualization)
@@ -353,137 +430,153 @@ int main(int argc, char *argv[])
             visit_dc.Save();
          }
       }
-      oper.SetParameters(u);
+      oper.SetConductionTensor(u);
    }
    tic_toc.Stop();
    cout << "Done, " << tic_toc.RealTime() << "s." << endl;
 
    // 9. Save the final solution. This output can be viewed later using GLVis:
    //    "glvis -m ex16.mesh -g ex16-final.gf".
-   {
-      ofstream osol("ex16-final.gf");
-      osol.precision(precision);
-      u_gf.Save(osol);
-   }
-
-   // 10. Free the used memory.
-   delete ode_solver;
-   delete mesh;
+   u_gf.Save("ex16-final.gf", precision);
 
    return 0;
 }
 
-ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al,
-                                       double kap, const Vector &u)
-   : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
-     T(NULL), z(height)
+ConductionOperator::ConductionOperator(FiniteElementSpace &fes,
+                                       const real_t alpha, const real_t kappa,
+                                       const Vector &u,
+                                       const Type &ode_expression_type)
+   : TimeDependentOperator(fes.GetTrueVSize(), 0.0, ode_expression_type),
+     fespace(fes), alpha(alpha), kappa(kappa), M(&fespace), z(height)
 {
-   const double rel_tol = 1e-8;
+   // specify a relative tolerance for all solves with MFEM integrators
+   const real_t rel_tol = 1e-8;
 
-   M = new BilinearForm(&fespace);
-   M->AddDomainIntegrator(new MassIntegrator());
-   M->Assemble();
-   M->FormSystemMatrix(ess_tdof_list, Mmat);
+   M.AddDomainIntegrator(new MassIntegrator());
+   M.Assemble();
+   M.FormSystemMatrix(ess_tdof_list, Mmat);
 
    M_solver.iterative_mode = false;
-   M_solver.SetRelTol(rel_tol);
+   M_solver.SetRelTol(rel_tol); // will be overwritten with SUNDIALS integrators
    M_solver.SetAbsTol(0.0);
    M_solver.SetMaxIter(50);
    M_solver.SetPrintLevel(0);
    M_solver.SetPreconditioner(M_prec);
    M_solver.SetOperator(Mmat);
 
-   alpha = al;
-   kappa = kap;
-
    T_solver.iterative_mode = false;
-   T_solver.SetRelTol(rel_tol);
+   T_solver.SetRelTol(rel_tol); // will be overwritten with SUNDIALS integrators
    T_solver.SetAbsTol(0.0);
    T_solver.SetMaxIter(100);
    T_solver.SetPrintLevel(0);
    T_solver.SetPreconditioner(T_prec);
 
-   SetParameters(u);
+   SetConductionTensor(u);
 }
 
-void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
-{
-   // Compute:
-   //    du_dt = M^{-1}*-K(u)
-   // for du_dt
-   Kmat.Mult(u, z);
-   z.Neg(); // z = -z
-   M_solver.Mult(z, du_dt);
-}
-
-void ConductionOperator::ImplicitSolve(const double dt,
-                                       const Vector &u, Vector &du_dt)
-{
-   // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt
-   if (T) { delete T; }
-   T = Add(1.0, Mmat, dt, Kmat);
-   T_solver.SetOperator(*T);
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
-}
-
-void ConductionOperator::SetParameters(const Vector &u)
+void ConductionOperator::SetConductionTensor(const Vector &u)
 {
+   // Compute K(u_n).
    GridFunction u_alpha_gf(&fespace);
    u_alpha_gf.SetFromTrueDofs(u);
    for (int i = 0; i < u_alpha_gf.Size(); i++)
    {
       u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
    }
-
-   delete K;
-   K = new BilinearForm(&fespace);
-
    GridFunctionCoefficient u_coeff(&u_alpha_gf);
 
+   K = std::make_unique<BilinearForm>(&fespace);
    K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
    K->Assemble();
    K->FormSystemMatrix(ess_tdof_list, Kmat);
 }
 
-int ConductionOperator::SUNImplicitSetup(const Vector &x,
-                                         const Vector &fx, int jok, int *jcur,
-                                         double gamma)
+void ConductionOperator::ExplicitMult(const Vector &u, Vector &v) const
 {
-   // Setup the ODE Jacobian T = M + gamma K.
-   if (T) { delete T; }
-   T = Add(1.0, Mmat, gamma, Kmat);
+   // Compute - K(u_n) u.
+   Kmat.Mult(u, v);
+   v.Neg();
+}
+
+void ConductionOperator::Mult(const Vector &u, Vector &k) const
+{
+   // Compute - inv(M) K(u_n) u.
+   ExplicitMult(u, z);
+   M_solver.Mult(z, k);
+}
+
+void ConductionOperator::ImplicitSolve(const real_t gam, const Vector &u,
+                                       Vector &k)
+{
+   // Solve for k in M k = - K(u_n) [u + gam*k].
+   ExplicitMult(u, z);
+   T = std::unique_ptr<SparseMatrix>(Add(1.0, Mmat, gam, Kmat));
    T_solver.SetOperator(*T);
-   *jcur = 1;
-   return (0);
+   T_solver.Mult(z, k);
 }
 
-int ConductionOperator::SUNImplicitSolve(const Vector &b, Vector &x, double tol)
+int ConductionOperator::SUNImplicitSetup(const Vector &u, const Vector &fu,
+                                         int jok, int *jcur, real_t gam)
 {
-   // Solve the system A x = z => (M - gamma K) x = M b.
-   Mmat.Mult(b, z);
-   T_solver.Mult(z, x);
-   return (0);
+   // Compute T = M + gamma K(u_n).
+   T = std::unique_ptr<SparseMatrix>(Add(1.0, Mmat, gam, Kmat));
+   T_solver.SetOperator(*T);
+   *jcur = SUNTRUE; // this should eventually only be set true if K(u) is used
+   return SUNLS_SUCCESS;
 }
 
-ConductionOperator::~ConductionOperator()
+int ConductionOperator::SUNImplicitSolve(const Vector &r, Vector &dk,
+                                         real_t tol)
 {
-   delete T;
-   delete M;
-   delete K;
-}
-
-double InitialTemperature(const Vector &x)
-{
-   if (x.Norml2() < 0.5)
+   // Solve the system [M + gamma K(u_n)] dk = - K(u_n) u - M k.
+   // What value r is providing depends on the ODE expression form:
+   //   EXPLICIT form: r = -inv(M) K(u_n) u - k
+   //   IMPLICIT form: r = -K(u_n) u - M k
+   T_solver.SetRelTol(tol);
+   if (isExplicit())
    {
-      return 2.0;
+      Mmat.Mult(r, z);
+      T_solver.Mult(z, dk);
    }
    else
    {
-      return 1.0;
+      T_solver.Mult(r, dk);
+   }
+   if (T_solver.GetConverged())
+   {
+      return SUNLS_SUCCESS;
+   }
+   else
+   {
+      return SUNLS_CONV_FAIL;
    }
 }
+
+int ConductionOperator::SUNMassSetup()
+{
+   // Do nothing b/c mass solver was setup in constructor.
+   return SUNLS_SUCCESS;
+}
+
+int ConductionOperator::SUNMassSolve(const Vector &b, Vector &x, real_t tol)
+{
+   // Solve the system M x = b.
+   M_solver.SetRelTol(tol);
+   M_solver.Mult(b, x);
+   if (M_solver.GetConverged())
+   {
+      return SUNLS_SUCCESS;
+   }
+   else
+   {
+      return SUNLS_CONV_FAIL;
+   }
+}
+
+int ConductionOperator::SUNMassMult(const Vector &x, Vector &v)
+{
+   // Compute M x.
+   Mmat.Mult(x, v);
+   return SUNLS_SUCCESS;
+}
+

examples/sundials/ex16p.cpp

@@ -1,16 +1,22 @@
 //                         MFEM Example 16 - Parallel Version
 //                              SUNDIALS Modification
 //
-// Compile with: make ex16p
+// Compile with:
+//    make ex16p            (GNU make)
+//    make sundials_ex16p   (CMake)
 //
 // Sample runs:
 //     mpirun -np 4 ex16p
 //     mpirun -np 4 ex16p -m ../../data/inline-tri.mesh
 //     mpirun -np 4 ex16p -m ../../data/disc-nurbs.mesh -tf 2
 //     mpirun -np 4 ex16p -s 12 -a 0.0 -k 1.0
+//     mpirun -np 4 ex16p -s 15 -a 0.0 -k 1.0
 //     mpirun -np 4 ex16p -s 8 -a 1.0 -k 0.0 -dt 4e-6 -tf 2e-2 -vs 50
+//     mpirun -np 4 ex16p -s 11 -a 1.0 -k 0.0 -dt 4e-6 -tf 2e-2 -vs 50
 //     mpirun -np 8 ex16p -s 9 -a 0.5 -k 0.5 -o 4 -dt 8e-6 -tf 2e-2 -vs 50
+//     mpirun -np 8 ex16p -s 12 -a 0.5 -k 0.5 -o 4 -dt 8e-6 -tf 2e-2 -vs 50
 //     mpirun -np 4 ex16p -s 10 -dt 2.0e-4 -tf 4.0e-2
+//     mpirun -np 4 ex16p -s 13 -dt 2.0e-4 -tf 4.0e-2
 //     mpirun -np 16 ex16p -m ../../data/fichera-q2.mesh
 //     mpirun -np 16 ex16p -m ../../data/escher-p2.mesh
 //     mpirun -np 8 ex16p -m ../../data/beam-tet.mesh -tf 10 -dt 0.1
@@ -38,66 +44,102 @@
 using namespace std;
 using namespace mfem;
 
-/** After spatial discretization, the conduction model can be written as:
+/** After spatial discretization, the conduction model is expressed as
  *
- *     du/dt = M^{-1}(-Ku)
+ *   M du/dt = - K(u) u
  *
  *  where u is the vector representing the temperature, M is the mass matrix,
- *  and K is the diffusion operator with diffusivity depending on u:
+ *  and K(u) is the diffusion operator with diffusivity depending on u:
  *  (\kappa + \alpha u).
  *
- *  Class ConductionOperator represents the right-hand side of the above ODE.
+ *  Class ConductionOperatorOperator represents the above ODE operator in the
+ *  general form F(u, k, t) = G(u, t) where either
+ *
+ *    1. F(u, du/dt, t) = du/dt            (ODE is expressed in EXPLICIT form)
+ *       G(u, t) = - inv(M) K(u) u
+ *    2. F(u, du/dt, t) = M du/dt          (ODE is expressed in IMPLICIT form)
+ *       G(u, t) = - K(u) u
  */
 class ConductionOperator : public TimeDependentOperator
 {
-protected:
    ParFiniteElementSpace &fespace;
    Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
 
-   ParBilinearForm *M;
-   ParBilinearForm *K;
-
+   ParBilinearForm M;
    HypreParMatrix Mmat;
+
+   const real_t alpha, kappa;
+   std::unique_ptr<BilinearForm> K;
    HypreParMatrix Kmat;
-   HypreParMatrix *T; // T = M + dt K
-   double current_dt;
 
-   CGSolver M_solver;    // Krylov solver for inverting the mass matrix M
-   HypreSmoother M_prec; // Preconditioner for the mass matrix M
+   std::unique_ptr<HypreParMatrix> T; // T = M + gam K(u)
 
-   CGSolver T_solver;    // Implicit solver for T = M + dt K
-   HypreSmoother T_prec; // Preconditioner for the implicit solver
+   CGSolver M_solver; // Krylov solver for inverting the mass matrix M
+   HypreSmoother M_prec;  // Preconditioner for the mass matrix M
 
-   double alpha, kappa;
+   CGSolver T_solver; // Implicit solver for T = M + gam K(u)
+   HypreSmoother T_prec;  // Preconditioner for the implicit solver
 
    mutable Vector z; // auxiliary vector
 
 public:
-   ConductionOperator(ParFiniteElementSpace &f, double alpha, double kappa,
-                      const Vector &u);
 
-   virtual void Mult(const Vector &u, Vector &du_dt) const;
+   ConductionOperator(ParFiniteElementSpace &f, const real_t alpha,
+                      const real_t kappa, const Vector &u,
+                      const Type &ode_expression_type);
 
-   /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
-       This is the only requirement for high-order SDIRK implicit integration.*/
-   virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
+   // Compute K(u_n) for use as an approximation in - K(u) u
+   void SetConductionTensor(const Vector &u);
 
-   /** Setup the system (M + dt K) x = M b. This method is used by the implicit
-       SUNDIALS solvers. */
-   virtual int SUNImplicitSetup(const Vector &x, const Vector &fx,
-                                int jok, int *jcur, double gamma);
+   /** Compute G(u, t) as defined in the IMPLICIT expression form of the ODE
+       operator, i.e., @a v = - K(u_n) @a u. Note that K(u_n) is an
+       approximation to K(u). */
+   void ExplicitMult(const Vector &u, Vector &v) const override;
 
-   /** Solve the system (M + dt K) x = M b. This method is used by the implicit
-       SUNDIALS solvers. */
-   virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol);
+   /** Solve for k in F(u, k, t) = G(u, t) for either EXPLICIT or IMPLICIT
+       expression forms of the ODE operator, i.e., @a k = - inv(M) K(u_n) @a u.
+       Note that K(u_n) is an approximation to K(u). */
+   void Mult(const Vector &u, Vector &k) const override;
 
-   /// Update the diffusion BilinearForm K using the given true-dof vector `u`.
-   void SetParameters(const Vector &u);
+   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t) for either EXPLICIT
+       or IMPLICIT expression forms of the ODE operator, i.e.,
+       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u . Note that K(u_n) is an
+       approximation to K(u). */
+   void ImplicitSolve(const real_t gam, const Vector &u, Vector &k) override;
 
-   virtual ~ConductionOperator();
+   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F for
+       either EXPLICIT or IMPLICIT expression forms of the ODE operator, i.e.,
+       [M - @a gam Jf(u)] dk = G - F, where Jf(u) is an approximation of the
+       Jacobian of -K(u) u. The approximation chosen here is Jf(u) = -K(u_n). */
+   int SUNImplicitSetup(const Vector &u, const Vector &fu, int jok, int *jcur,
+                        real_t gam) override;
+
+   /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
+       with the residual @a r providing either
+        1. @a r = G - F = inv(M) f(u) - k       (EXPLICIT expression form)
+        1. @a r = G - F = f(u) - M k            (IMPLICIT expression form)
+       */
+   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
+
+   int SUNMassSetup() override;
+
+   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override;
+
+   int SUNMassMult(const Vector &x, Vector &v) override;
 };
 
-double InitialTemperature(const Vector &x);
+real_t InitialTemperature(const Vector &x)
+{
+   if (x.Norml2() < 0.5)
+   {
+      return 2.0;
+   }
+   else
+   {
+      return 1.0;
+   }
+}
+
 
 int main(int argc, char *argv[])
 {
@@ -114,16 +156,16 @@ int main(int argc, char *argv[])
    int par_ref_levels = 1;
    int order = 2;
    int ode_solver_type = 9; // CVODE implicit BDF
-   double t_final = 0.5;
-   double dt = 1.0e-2;
-   double alpha = 1.0e-2;
-   double kappa = 0.5;
+   real_t t_final = 0.5;
+   real_t dt = 1.0e-2;
+   real_t alpha = 1.0e-2;
+   real_t kappa = 0.5;
    bool visualization = true;
    bool visit = false;
    int vis_steps = 5;
 
    // Relative and absolute tolerances for CVODE and ARKODE.
-   const double reltol = 1e-4, abstol = 1e-4;
+   const real_t reltol = 1e-4, abstol = 1e-4;
 
    int precision = 8;
    cout.precision(precision);
@@ -150,7 +192,10 @@ int main(int argc, char *argv[])
                   "9  - CVODE (implicit BDF),\n\t"
                   "10 - ARKODE (default explicit),\n\t"
                   "11 - ARKODE (explicit Fehlberg-6-4-5),\n\t"
-                  "12 - ARKODE (default impicit).");
+                  "12 - ARKODE (default implicit),\n\t"
+                  "13 - ARKODE (default explicit with MFEM mass solve),\n\t"
+                  "14 - ARKODE (explicit Fehlberg-6-4-5 with MFEM mass solve),\n\t"
+                  "15 - ARKODE (default implicit with MFEM mass solve).");
    args.AddOption(&t_final, "-tf", "--t-final",
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
@@ -174,40 +219,33 @@ int main(int argc, char *argv[])
       return 1;
    }
 
-   if (myid == 0)
+   if (Mpi::Root())
    {
       args.PrintOptions(cout);
    }
 
-   // check for valid ODE solver option
-   if (ode_solver_type < 1 || ode_solver_type > 12)
-   {
-      if (myid == 0)
-      {
-         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
-      }
-      return 1;
-   }
+   bool use_mass_solver = ode_solver_type >= 13;
 
-   // 3. Read the serial mesh from the given mesh file on all processors. We can
+   // 3. Define a parallel mesh by a partitioning of a serial mesh. Read the
+   //    serial mesh from the given mesh file on all processors. We can
    //    handle triangular, quadrilateral, tetrahedral and hexahedral meshes
    //    with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-
-   // 4. Refine the mesh in serial to increase the resolution. In this example
-   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
-   //    a command-line parameter.
-   for (int lev = 0; lev < ser_ref_levels; lev++)
+   std::unique_ptr<ParMesh> pmesh;
    {
-      mesh->UniformRefinement();
-   }
+      std::unique_ptr<Mesh> mesh(new Mesh(mesh_file, 1, 1));
 
-   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted.
-   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
+      // 4. Refine the mesh in serial to increase the resolution. In this example
+      //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+      //    a command-line parameter.
+      for (int lev = 0; lev < ser_ref_levels; lev++)
+      {
+         mesh->UniformRefinement();
+      }
+
+      // 5. Refine this mesh further in parallel to increase the resolution.
+      //    Once the parallel mesh is defined, the serial mesh can be deleted.
+      pmesh = std::make_unique<ParMesh>(MPI_COMM_WORLD, *mesh);
+   }
    for (int lev = 0; lev < par_ref_levels; lev++)
    {
       pmesh->UniformRefinement();
@@ -215,8 +253,9 @@ int main(int argc, char *argv[])
 
    // 6. Define the vector finite element space representing the current and the
    //    initial temperature, u_ref.
+   int dim = pmesh->Dimension();
    H1_FECollection fe_coll(order, dim);
-   ParFiniteElementSpace fespace(pmesh, &fe_coll);
+   ParFiniteElementSpace fespace(pmesh.get(), &fe_coll);
 
    int fe_size = fespace.GlobalTrueVSize();
    if (myid == 0)
@@ -233,8 +272,17 @@ int main(int argc, char *argv[])
    Vector u;
    u_gf.GetTrueDofs(u);
 
-   // 8. Initialize the conduction operator and the VisIt visualization.
-   ConductionOperator oper(fespace, alpha, kappa, u);
+   // 8. Initialize the conduction ODE operator and the visualization.
+   ConductionOperator::Type ode_expression_type;
+   if (use_mass_solver)
+   {
+      ode_expression_type = ConductionOperator::Type::IMPLICIT;
+   }
+   else
+   {
+      ode_expression_type = ConductionOperator::Type::EXPLICIT;
+   }
+   ConductionOperator oper(fespace, alpha, kappa, u, ode_expression_type);
 
    u_gf.SetFromTrueDofs(u);
    {
@@ -249,7 +297,7 @@ int main(int argc, char *argv[])
       u_gf.Save(osol);
    }
 
-   VisItDataCollection visit_dc("Example16-Parallel", pmesh);
+   VisItDataCollection visit_dc("Example16-Parallel", pmesh.get());
    visit_dc.RegisterField("temperature", &u_gf);
    if (visit)
    {
@@ -293,52 +341,76 @@ int main(int argc, char *argv[])
    }
 
    // 9. Define the ODE solver used for time integration.
-   double t = 0.0;
-   ODESolver *ode_solver = NULL;
-   CVODESolver *cvode = NULL;
-   ARKStepSolver *arkode = NULL;
+   real_t t = 0.0;
+   std::unique_ptr<ODESolver> ode_solver;
    switch (ode_solver_type)
    {
       // MFEM explicit methods
-      case 1: ode_solver = new ForwardEulerSolver; break;
-      case 2: ode_solver = new RK2Solver(0.5); break; // midpoint method
-      case 3: ode_solver = new RK3SSPSolver; break;
-      case 4: ode_solver = new RK4Solver; break;
+      case 1: ode_solver = std::make_unique<ForwardEulerSolver>(); break;
+      case 2: ode_solver = std::make_unique<RK2Solver>(0.5); break; // midpoint method
+      case 3: ode_solver = std::make_unique<RK3SSPSolver>(); break;
+      case 4: ode_solver = std::make_unique<RK4Solver>(); break;
       // MFEM implicit L-stable methods
-      case 5: ode_solver = new BackwardEulerSolver; break;
-      case 6: ode_solver = new SDIRK23Solver(2); break;
-      case 7: ode_solver = new SDIRK33Solver; break;
+      case 5: ode_solver = std::make_unique<BackwardEulerSolver>(); break;
+      case 6: ode_solver = std::make_unique<SDIRK23Solver>(2); break;
+      case 7: ode_solver = std::make_unique<SDIRK33Solver>(); break;
       // CVODE
       case 8:
-         cvode = new CVODESolver(MPI_COMM_WORLD, CV_ADAMS);
-         cvode->Init(oper);
-         cvode->SetSStolerances(reltol, abstol);
-         cvode->SetMaxStep(dt);
-         ode_solver = cvode; break;
       case 9:
-         cvode = new CVODESolver(MPI_COMM_WORLD, CV_BDF);
+      {
+         int cvode_solver_type;
+         if (ode_solver_type == 8)
+         {
+            cvode_solver_type = CV_ADAMS;
+         }
+         else
+         {
+            cvode_solver_type = CV_BDF;
+         }
+         std::unique_ptr<CVODESolver> cvode(
+            new CVODESolver(MPI_COMM_WORLD, cvode_solver_type));
          cvode->Init(oper);
          cvode->SetSStolerances(reltol, abstol);
          cvode->SetMaxStep(dt);
-         ode_solver = cvode; break;
+         ode_solver = std::move(cvode);
+         break;
+      }
       // ARKODE
       case 10:
       case 11:
-         arkode = new ARKStepSolver(MPI_COMM_WORLD, ARKStepSolver::EXPLICIT);
+      case 12:
+      case 13:
+      case 14:
+      case 15:
+      {
+         ARKStepSolver::Type arkode_solver_type;
+         if (ode_solver_type == 12 || ode_solver_type == 15)
+         {
+            arkode_solver_type = ARKStepSolver::IMPLICIT;
+         }
+         else
+         {
+            arkode_solver_type = ARKStepSolver::EXPLICIT;
+         }
+         std::unique_ptr<ARKStepSolver> arkode(
+            new ARKStepSolver(MPI_COMM_WORLD, arkode_solver_type));
          arkode->Init(oper);
          arkode->SetSStolerances(reltol, abstol);
          arkode->SetMaxStep(dt);
-         if (ode_solver_type == 11)
+         if (ode_solver_type == 11 || ode_solver_type == 14)
          {
             arkode->SetERKTableNum(ARKODE_FEHLBERG_13_7_8);
          }
-         ode_solver = arkode; break;
-      case 12:
-         arkode = new ARKStepSolver(MPI_COMM_WORLD, ARKStepSolver::IMPLICIT);
-         arkode->Init(oper);
-         arkode->SetSStolerances(reltol, abstol);
-         arkode->SetMaxStep(dt);
-         ode_solver = arkode; break;
+         if (use_mass_solver)
+         {
+            arkode->UseMFEMMassLinearSolver(SUNFALSE);
+         }
+         ode_solver = std::move(arkode);
+         break;
+      }
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
    }
 
    // Initialize MFEM integrators, SUNDIALS integrators are initialized above
@@ -346,12 +418,18 @@ int main(int argc, char *argv[])
 
    // Since we want to update the diffusion coefficient after every time step,
    // we need to use the "one-step" mode of the SUNDIALS solvers.
-   if (cvode) { cvode->SetStepMode(CV_ONE_STEP); }
-   if (arkode) { arkode->SetStepMode(ARK_ONE_STEP); }
+   if (CVODESolver* cvode = dynamic_cast<CVODESolver*>(ode_solver.get()))
+   {
+      cvode->SetStepMode(CV_ONE_STEP);
+   }
+   else if (ARKStepSolver* arkode = dynamic_cast<ARKStepSolver*>(ode_solver.get()))
+   {
+      arkode->SetStepMode(ARK_ONE_STEP);
+   }
 
    // 10. Perform time-integration (looping over the time iterations, ti, with a
    //     time-step dt).
-   if (myid == 0)
+   if (Mpi::Root())
    {
       cout << "Integrating the ODE ..." << endl;
    }
@@ -361,7 +439,7 @@ int main(int argc, char *argv[])
    bool last_step = false;
    for (int ti = 1; !last_step; ti++)
    {
-      double dt_real = min(dt, t_final - t);
+      real_t dt_real = min(dt, t_final - t);
 
       // Note that since we are using the "one-step" mode of the SUNDIALS
       // solvers, they will, generally, step over the final time and will not
@@ -377,8 +455,14 @@ int main(int argc, char *argv[])
          if (myid == 0)
          {
             cout << "step " << ti << ", t = " << t << endl;
-            if (cvode) { cvode->PrintInfo(); }
-            if (arkode) { arkode->PrintInfo(); }
+            if (CVODESolver* cvode = dynamic_cast<CVODESolver*>(ode_solver.get()))
+            {
+               cvode->PrintInfo();
+            }
+            else if (ARKStepSolver* arkode = dynamic_cast<ARKStepSolver*>(ode_solver.get()))
+            {
+               arkode->PrintInfo();
+            }
          }
 
          u_gf.SetFromTrueDofs(u);
@@ -395,46 +479,38 @@ int main(int argc, char *argv[])
             visit_dc.Save();
          }
       }
-      oper.SetParameters(u);
+      oper.SetConductionTensor(u);
    }
    tic_toc.Stop();
-   if (myid == 0)
+   if (Mpi::Root())
    {
       cout << "Done, " << tic_toc.RealTime() << "s." << endl;
    }
 
    // 11. Save the final solution in parallel. This output can be viewed later
    //     using GLVis: "glvis -np <np> -m ex16-mesh -g ex16-final".
-   {
-      ostringstream sol_name;
-      sol_name << "ex16-final." << setfill('0') << setw(6) << myid;
-      ofstream osol(sol_name.str().c_str());
-      osol.precision(precision);
-      u_gf.Save(osol);
-   }
-
-   // 12. Free the used memory.
-   delete ode_solver;
-   delete pmesh;
+   u_gf.Save("ex16-final", precision);
 
    return 0;
 }
 
-ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
-                                       double kap, const Vector &u)
-   : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
-     T(NULL),
-     M_solver(f.GetComm()), T_solver(f.GetComm()), z(height)
+ConductionOperator::ConductionOperator(ParFiniteElementSpace &fes,
+                                       const real_t alpha, const real_t kappa,
+                                       const Vector &u,
+                                       const Type &ode_expression_type)
+   : TimeDependentOperator(fes.GetTrueVSize(), 0.0, ode_expression_type),
+     fespace(fes), alpha(alpha), kappa(kappa), M(&fespace),
+     M_solver(fes.GetComm()), T_solver(fes.GetComm()), z(height)
 {
-   const double rel_tol = 1e-8;
+   // specify a relative tolerance for all solves with MFEM integrators
+   const real_t rel_tol = 1e-8;
 
-   M = new ParBilinearForm(&fespace);
-   M->AddDomainIntegrator(new MassIntegrator());
-   M->Assemble(0); // keep sparsity pattern of M and K the same
-   M->FormSystemMatrix(ess_tdof_list, Mmat);
+   M.AddDomainIntegrator(new MassIntegrator());
+   M.Assemble(0); // keep zeros to keep sparsity pattern of M and K the same
+   M.FormSystemMatrix(ess_tdof_list, Mmat);
 
    M_solver.iterative_mode = false;
-   M_solver.SetRelTol(rel_tol);
+   M_solver.SetRelTol(rel_tol); // will be overwritten with SUNDIALS integrators
    M_solver.SetAbsTol(0.0);
    M_solver.SetMaxIter(100);
    M_solver.SetPrintLevel(0);
@@ -442,97 +518,118 @@ ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
    M_solver.SetPreconditioner(M_prec);
    M_solver.SetOperator(Mmat);
 
-   alpha = al;
-   kappa = kap;
-
    T_solver.iterative_mode = false;
-   T_solver.SetRelTol(rel_tol);
+   T_solver.SetRelTol(rel_tol); // will be overwritten with SUNDIALS integrators
    T_solver.SetAbsTol(0.0);
    T_solver.SetMaxIter(100);
    T_solver.SetPrintLevel(0);
    T_solver.SetPreconditioner(T_prec);
 
-   SetParameters(u);
+   SetConductionTensor(u);
 }
 
-void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
-{
-   // Compute:
-   //    du_dt = M^{-1}*-K(u)
-   // for du_dt
-   Kmat.Mult(u, z);
-   z.Neg(); // z = -z
-   M_solver.Mult(z, du_dt);
-}
-
-void ConductionOperator::ImplicitSolve(const double dt,
-                                       const Vector &u, Vector &du_dt)
-{
-   // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt
-   if (T) { delete T; }
-   T = Add(1.0, Mmat, dt, Kmat);
-   T_solver.SetOperator(*T);
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
-}
-
-int ConductionOperator::SUNImplicitSetup(const Vector &x,
-                                         const Vector &fx, int jok, int *jcur,
-                                         double gamma)
-{
-   // Setup the ODE Jacobian T = M + gamma K.
-   if (T) { delete T; }
-   T = Add(1.0, Mmat, gamma, Kmat);
-   T_solver.SetOperator(*T);
-   *jcur = 1;
-   return (0);
-}
-
-int ConductionOperator::SUNImplicitSolve(const Vector &b, Vector &x, double tol)
-{
-   // Solve the system A x = z => (M - gamma K) x = M b.
-   Mmat.Mult(b, z);
-   T_solver.Mult(z, x);
-   return (0);
-}
-
-void ConductionOperator::SetParameters(const Vector &u)
+void ConductionOperator::SetConductionTensor(const Vector &u)
 {
+   // Compute K(u_n).
    ParGridFunction u_alpha_gf(&fespace);
    u_alpha_gf.SetFromTrueDofs(u);
    for (int i = 0; i < u_alpha_gf.Size(); i++)
    {
       u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
    }
-
-   delete K;
-   K = new ParBilinearForm(&fespace);
-
    GridFunctionCoefficient u_coeff(&u_alpha_gf);
 
+   K = std::make_unique<ParBilinearForm>(&fespace);
    K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
-   K->Assemble(0); // keep sparsity pattern of M and K the same
+   K->Assemble(0); // keep zeros to keep sparsity pattern of M and K the same
    K->FormSystemMatrix(ess_tdof_list, Kmat);
 }
 
-ConductionOperator::~ConductionOperator()
+void ConductionOperator::ExplicitMult(const Vector &u, Vector &v) const
 {
-   delete T;
-   delete M;
-   delete K;
+   // Compute - K(u_n) u.
+   Kmat.Mult(u, v);
+   v.Neg();
 }
 
-double InitialTemperature(const Vector &x)
+void ConductionOperator::Mult(const Vector &u, Vector &k) const
 {
-   if (x.Norml2() < 0.5)
+   // Compute - inv(M) K(u_n) u.
+   ExplicitMult(u, z);
+   M_solver.Mult(z, k);
+}
+
+void ConductionOperator::ImplicitSolve(const real_t gam, const Vector &u,
+                                       Vector &k)
+{
+   // Solve for k in M k = - K(u_n) [u + gam*k].
+   ExplicitMult(u, z);
+   T = std::unique_ptr<HypreParMatrix>(Add(1.0, Mmat, gam, Kmat));
+   T_solver.SetOperator(*T);
+   T_solver.Mult(z, k);
+}
+
+int ConductionOperator::SUNImplicitSetup(const Vector &u, const Vector &fu,
+                                         int jok, int *jcur, real_t gam)
+{
+   // Compute T = M + gamma K(u_n).
+   T = std::unique_ptr<HypreParMatrix>(Add(1.0, Mmat, gam, Kmat));
+   T_solver.SetOperator(*T);
+   *jcur = SUNTRUE; // this should eventually only be set true if K(u) is used
+   return SUNLS_SUCCESS;
+}
+
+int ConductionOperator::SUNImplicitSolve(const Vector &r, Vector &dk,
+                                         real_t tol)
+{
+   // Solve the system [M + gamma K(u_n)] dk = - K(u_n) u - M k.
+   // What value r is providing depends on the ODE expression form:
+   //   EXPLICIT form: r = -inv(M) K(u_n) u - k
+   //   IMPLICIT form: r = -K(u_n) u - M k
+   T_solver.SetRelTol(tol);
+   if (isExplicit())
    {
-      return 2.0;
+      Mmat.Mult(r, z);
+      T_solver.Mult(z, dk);
    }
    else
    {
-      return 1.0;
+      T_solver.Mult(r, dk);
+   }
+   if (T_solver.GetConverged())
+   {
+      return SUNLS_SUCCESS;
+   }
+   else
+   {
+      return SUNLS_CONV_FAIL;
    }
 }
+
+int ConductionOperator::SUNMassSetup()
+{
+   // Do nothing b/c mass solver was setup in constructor.
+   return SUNLS_SUCCESS;
+}
+
+int ConductionOperator::SUNMassSolve(const Vector &b, Vector &x, real_t tol)
+{
+   // Solve the system M x = b.
+   M_solver.SetRelTol(tol);
+   M_solver.Mult(b, x);
+   if (M_solver.GetConverged())
+   {
+      return SUNLS_SUCCESS;
+   }
+   else
+   {
+      return SUNLS_CONV_FAIL;
+   }
+}
+
+int ConductionOperator::SUNMassMult(const Vector &x, Vector &v)
+{
+   // Compute M x.
+   Mmat.Mult(x, v);
+   return SUNLS_SUCCESS;
+}

examples/sundials/ex9.cpp

@@ -1,7 +1,9 @@
 //                                MFEM Example 9
 //                             SUNDIALS Modification
 //
-// Compile with: make ex9
+// Compile with:
+//    make ex9              (GNU make)
+//    make sundials_ex9     (CMake)
 //
 // Sample runs:
 //    ex9 -m ../../data/periodic-segment.mesh -p 0 -r 2 -s 7 -dt 0.005

examples/sundials/ex9p.cpp

@@ -1,7 +1,9 @@
 //                       MFEM Example 9 - Parallel Version
 //                             SUNDIALS Modification
 //
-// Compile with: make ex9p
+// Compile with:
+//    make ex9p             (GNU make)
+//    make sundials_ex9p    (CMake)
 //
 // Sample runs:
 //    mpirun -np 4 ex9p -m ../../data/periodic-segment.mesh -p 1 -rp 1 -s 7 -dt 0.0025

examples/sundials/makefile

@@ -100,6 +100,12 @@ ex10-test-seq: ex10
 	@$(call mfem-test,$<,, $(SERIAL_NAME),$(EX10_ARGS))
 ex10p-test-par: ex10p
 	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_NAME),$(EX10P_ARGS))
+# Example 16: test ARKODE with implicit time stepping using mass form
+EX16_COMMON_ARGS := -s 15
+ex16-test-seq: ex16
+	@$(call mfem-test,$<,, $(SERIAL_NAME),$(EX16_COMMON_ARGS))
+ex16p-test-par: ex16p
+	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_NAME),$(EX16_COMMON_ARGS))
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 

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/gslib.cpp

@@ -1168,6 +1168,120 @@ void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,
    } // parallel
 }
 
+void FindPointsGSLIB::DistributePointInfoToOwningMPIRanks(
+   Array<unsigned int> &recv_elem, Vector &recv_ref,
+   Array<unsigned int> &recv_code)
+{
+   MFEM_VERIFY(points_cnt,
+               "Invalid size. Please make sure to call FindPoints method "
+               "before calling this function.");
+
+   // Pack data to send via crystal router
+   struct gslib::array *outpt = new gslib::array;
+
+   struct out_pt { double rst[3]; uint index, elem, proc, code; };
+   struct out_pt *pt;
+   array_init(struct out_pt, outpt, points_cnt);
+   outpt->n=points_cnt;
+   pt = (struct out_pt *)outpt->ptr;
+
+   for (int index = 0; index < points_cnt; index++)
+   {
+      pt->index = index;
+      pt->elem = gsl_mfem_elem[index];
+      pt->proc  = gsl_proc[index];
+      pt->code = gsl_code[index];
+      for (int d = 0; d < dim; ++d)
+      {
+         pt->rst[d]= gsl_mfem_ref(index*dim + d);
+      }
+      ++pt;
+   }
+
+   // Transfer data to target MPI ranks
+   sarray_transfer(struct out_pt, outpt, proc, 1, cr);
+
+   // Store received data
+   const int points_recv = outpt->n;
+   recv_proc.SetSize(points_recv);
+   recv_elem.SetSize(points_recv);
+   recv_index.SetSize(points_recv);
+   recv_code.SetSize(points_recv);
+   recv_ref.SetSize(points_recv*dim);
+
+   pt = (struct out_pt *)outpt->ptr;
+   for (int index = 0; index < points_recv; index++)
+   {
+      recv_index[index] = pt->index;
+      recv_elem[index] = pt->elem;
+      recv_proc[index] = pt->proc;
+      recv_code[index] = pt->code;
+      for (int d = 0; d < dim; ++d)
+      {
+         recv_ref(index*dim + d)= pt->rst[d];
+      }
+      ++pt;
+   }
+
+   array_free(outpt);
+   delete outpt;
+}
+
+void FindPointsGSLIB::DistributeInterpolatedValues(const Vector &int_vals,
+                                                   const int vdim,
+                                                   const int ordering,
+                                                   Vector &field_out) const
+{
+   const int points_recv = recv_index.Size();;
+   MFEM_VERIFY(points_recv == 0 ||
+               int_vals.Size() % points_recv == 0,
+               "Incompatible size. Please return interpolated values"
+               "corresponding to points received using"
+               "SendCoordinatesToOwningProcessors.");
+   field_out.SetSize(points_cnt*vdim);
+
+   for (int v = 0; v < vdim; v++)
+   {
+      // Pack data to send via crystal router
+      struct gslib::array *outpt = new gslib::array;
+      struct out_pt { double val; uint index, proc; };
+      struct out_pt *pt;
+      array_init(struct out_pt, outpt, points_recv);
+      outpt->n=points_recv;
+      pt = (struct out_pt *)outpt->ptr;
+      for (int index = 0; index < points_recv; index++)
+      {
+         pt->index = recv_index[index];
+         pt->proc  = recv_proc[index];
+         pt->val = ordering == Ordering::byNODES ?
+                   int_vals(index + v*points_recv) :
+                   int_vals(index*vdim + v);
+         ++pt;
+      }
+
+      // Transfer data to target MPI ranks
+      sarray_transfer(struct out_pt, outpt, proc, 1, cr);
+
+      // Store received data
+      MFEM_VERIFY(outpt->n == points_cnt, "Incompatible size. Number of points "
+                  "received does not match the number of points originally "
+                  "found using FindPoints.");
+
+      pt = (struct out_pt *)outpt->ptr;
+      for (int index = 0; index < points_cnt; index++)
+      {
+         int idx = ordering == Ordering::byNODES ?
+                   pt->index + v*points_cnt :
+                   pt->index*vdim + v;
+         field_out(idx) = pt->val;
+         ++pt;
+      }
+
+      array_free(outpt);
+      delete outpt;
+   }
+}
+
 void OversetFindPointsGSLIB::Setup(Mesh &m, const int meshid,
                                    GridFunction *gfmax,
                                    const double bb_t, const double newt_tol,

fem/gslib.hpp

@@ -34,7 +34,7 @@ namespace mfem
  *
  *  There are three key functions in FindPointsGSLIB:
  *
- *  1. Setup - constructs the internal data structures of gslib.
+ *  1. Setup - constructs the internal data structures of gslib. See \ref Setup.
  *
  *  2. FindPoints - for any given arbitrary set of points in physical space,
  *     gslib finds the element number, MPI rank, and the reference space
@@ -45,12 +45,23 @@ namespace mfem
  *     on an element edge/face or near the domain boundary, and gslib also
  *     returns a distance to the border. Points near (but outside) the domain
  *     boundary must then be marked as not found using the distance returned
- *     by gslib.
+ *     by gslib. See \ref FindPoints.
  *
  *  3. Interpolate - Interpolates any grid function at the points found using 2.
+ *     For functions in L2 finite element space, use \ref SetL2AvgType to
+ *     specify how to interpolate values at points located at element boundaries
+ *     where the function might be multi-valued. See \ref Interpolate.
  *
- *  FindPointsGSLIB provides interface to use these functions individually or
- *  using a single call.
+ *  FindPointsGSLIB also provides interface to use these functions through a
+ *  single call.
+ *
+ *  For custom interpolation (e.g., evaluating strain rate tensor), we provide
+ *  functions that use gslib to send element index and corresponding
+ *  reference-space coordinates for each point to the mpi rank that the element
+ *  is located on. Then, custom interpolation can be defined locally by the user
+ *  before sending the values back to mpi ranks where the query originated from.
+ *  See \ref DistributePointInfoToOwningMPIRanks and
+ *  \ref DistributeInterpolatedValues.
  */
 class FindPointsGSLIB
 {
@@ -74,7 +85,8 @@ protected:
    int dim, points_cnt;
    Array<unsigned int> gsl_code, gsl_proc, gsl_elem, gsl_mfem_elem;
    Vector gsl_mesh, gsl_ref, gsl_dist, gsl_mfem_ref;
-   bool setupflag;              // flag to indicate whether gslib data has been setup
+   Array<unsigned int> recv_proc, recv_index; // data for custom interpolation
+   bool setupflag;              // flag to indicate if gslib data has been setup
    double default_interp_value; // used for points that are not found in the mesh
    AvgType avgtype;             // average type used for L2 functions
    Array<int> split_element_map;
@@ -118,9 +130,9 @@ public:
    virtual ~FindPointsGSLIB();
 
    /** Initializes the internal mesh in gslib, by sending the positions of the
-       Gauss-Lobatto nodes of the input Mesh object @a m.
+       Gauss-Lobatto nodes of the input Mesh object \p m.
        Note: not tested with periodic (L2).
-       Note: the input mesh @a m must have Nodes set.
+       Note: the input mesh \p m must have Nodes set.
 
        @param[in] m         Input mesh.
        @param[in] bb_t      (Optional) Relative size of bounding box around
@@ -133,9 +145,9 @@ public:
    void Setup(Mesh &m, const double bb_t = 0.1,
               const double newt_tol = 1.0e-12,
               const int npt_max = 256);
-   /** Searches positions given in physical space by @a point_pos.
+   /** Searches positions given in physical space by \p point_pos.
        These positions can be ordered byNodes: (XXX...,YYY...,ZZZ) or
-       byVDim: (XYZ,XYZ,....XYZ) specified by @a point_pos_ordering.
+       byVDim: (XYZ,XYZ,....XYZ) specified by \p point_pos_ordering.
        This function populates the following member variables:
        #gsl_code        Return codes for each point: inside element (0),
                         element boundary (1), not found (2).
@@ -164,20 +176,20 @@ public:
    /** Interpolation of field values at prescribed reference space positions.
        @param[in] field_in    Function values that will be interpolated on the
                               reference positions. Note: it is assumed that
-                              @a field_in is in H1 and in the same space as the
+                              \p field_in is in H1 and in the same space as the
                               mesh that was given to Setup().
        @param[out] field_out  Interpolated values. For points that are not found
                               the value is set to #default_interp_value. */
    virtual void Interpolate(const GridFunction &field_in, Vector &field_out);
    /** Search positions and interpolate. The ordering (byNODES or byVDIM) of
-       the output values in @a field_out corresponds to the ordering used
-       in the input GridFunction @a field_in. */
+       the output values in \p field_out corresponds to the ordering used
+       in the input GridFunction \p field_in. */
    void Interpolate(const Vector &point_pos, const GridFunction &field_in,
                     Vector &field_out,
                     int point_pos_ordering = Ordering::byNODES);
    /** Setup FindPoints, search positions and interpolate. The ordering (byNODES
-       or byVDIM) of the output values in @a field_out corresponds to the
-       ordering used in the input GridFunction @a field_in. */
+       or byVDIM) of the output values in \p field_out corresponds to the
+       ordering used in the input GridFunction \p field_in. */
    void Interpolate(Mesh &m, const Vector &point_pos,
                     const GridFunction &field_in, Vector &field_out,
                     int point_pos_ordering = Ordering::byNODES);
@@ -225,6 +237,41 @@ public:
    /// Return reference coordinates in [-1,1] (internal range in GSLIB) for each
    /// point found by FindPoints.
    virtual const Vector &GetGSLIBReferencePosition() const { return gsl_ref; }
+
+   /** @name Methods to support a custom interpolation procedure.
+       \brief The physical-space point that the user seeks to interpolate at
+       could be located inside an element on another mpi rank.
+       To enable a custom interpolation procedure (e.g., strain tensor computation)
+       we need a mechanism to first send element indices and reference-space
+       coordinates to the mpi-ranks where each point is found. Then the custom
+       interpolation can be done locally by the user before sending the
+       interpolated values back to the mpi-ranks that the query originated from.
+       Example usage looks something like this:
+
+       FindPoints() -> DistributePointInfoToOwningMPIRanks() -> Computation by
+       user -> DistributeInterpolatedValues().
+   */
+   ///@{
+   /// Distribute element indices in #gsl_mfem_elem, the reference coordinates
+   /// #gsl_mfem_ref, and the code #gsl_code to the corresponding mpi-rank
+   /// #gsl_proc for each point. The received information is provided locally
+   /// in \p recv_elem, \p recv_ref (ordered by vdim), and \p recv_code.
+   /// Note: The user can send empty Array/Vectors to the method as they are
+   /// appropriately sized and filled internally.
+   virtual void DistributePointInfoToOwningMPIRanks(
+      Array<unsigned int> &recv_elem, Vector &recv_ref,
+      Array<unsigned int> &recv_code);
+   /// Return interpolated values back to the mpi-ranks #recv_proc that had
+   /// sent the element indices and corresponding reference-space coordinates.
+   /// Specify \p vdim and \p ordering (by nodes or by vdim) based on how the
+   /// \p int_vals are structured. The received values are filled in
+   /// \p field_out consistent with the original ordering of the points that
+   /// were used in \ref FindPoints.
+   virtual void DistributeInterpolatedValues(const Vector &int_vals,
+                                             const int vdim,
+                                             const int ordering,
+                                             Vector &field_out) const;
+   ///@}
 };
 
 /** \brief OversetFindPointsGSLIB enables use of findpts for arbitrary number of
@@ -249,9 +296,9 @@ public:
 #endif
 
    /** Initializes the internal mesh in gslib, by sending the positions of the
-       Gauss-Lobatto nodes of the input Mesh object @a m.
+       Gauss-Lobatto nodes of the input Mesh object \p m.
        Note: not tested with periodic meshes (L2).
-       Note: the input mesh @a m must have Nodes set.
+       Note: the input mesh \p m must have Nodes set.
 
        @param[in] m         Input mesh.
        @param[in] meshid    A unique # for each overlapping mesh. This id is
@@ -274,12 +321,12 @@ public:
               const double bb_t = 0.1, const double newt_tol = 1.0e-12,
               const int npt_max = 256);
 
-   /** Searches positions given in physical space by @a point_pos. All output
+   /** Searches positions given in physical space by \p point_pos. All output
        Arrays and Vectors are expected to have the correct size.
 
        @param[in]  point_pos           Positions to be found.
        @param[in]  point_id            Index of the mesh that the point belongs
-                                       to (corresponding to @a meshid in Setup).
+                                       to (corresponding to \p meshid in Setup).
        @param[in]  point_pos_ordering  Ordering of the points:
                                        byNodes: (XXX...,YYY...,ZZZ) or
                                        byVDim: (XYZ,XYZ,....XYZ) */
@@ -342,7 +389,7 @@ public:
    enum GSOp {ADD, MUL, MIN, MAX};
 
    /// Update the identifiers used for the gather-scatter operator.
-   /// Same @a ids get grouped together and id == 0 does not participate.
+   /// Same \p ids get grouped together and id == 0 does not participate.
    /// See class description.
    void UpdateIdentifiers(const Array<long long> &ids);
 

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

fem/restriction.cpp

@@ -855,7 +855,8 @@ void ConformingFaceRestriction::SetFaceDofsScatterIndices(
                "This method should not be used on nonconforming coarse faces.");
    MFEM_ASSERT(face.element[0].orientation==0,
                "FaceRestriction used on degenerated mesh.");
-   MFEM_CONTRACT_VAR(f_ordering); // not supported yet
+   MFEM_VERIFY(f_ordering == ElementDofOrdering::LEXICOGRAPHIC,
+               "NATIVE ordering is not supported yet");
 
    fes.GetFE(0)->GetFaceMap(face.element[0].local_face_id, face_map);
 
@@ -883,7 +884,8 @@ void ConformingFaceRestriction::SetFaceDofsGatherIndices(
 {
    MFEM_ASSERT(!(face.IsNonconformingCoarse()),
                "This method should not be used on nonconforming coarse faces.");
-   MFEM_CONTRACT_VAR(f_ordering); // not supported yet
+   MFEM_VERIFY(f_ordering == ElementDofOrdering::LEXICOGRAPHIC,
+               "NATIVE ordering is not supported yet");
 
    fes.GetFE(0)->GetFaceMap(face.element[0].local_face_id, face_map);
 

linalg/operator.hpp

@@ -420,8 +420,9 @@ public:
    /** @brief Perform the action of the explicit part of the operator, G:
        @a v = G(@a u, t) where t is the current time.
 
-       Presently, this method is used by some PETSc ODE solvers, for more
-       details, see the PETSc Manual. */
+       Presently, this method is used by some PETSc ODE solvers and the
+       SUNDIALS ARKStep integrator, for more details, see either the PETSc
+       Manual or the ARKode User Guide, respectively. */
    virtual void ExplicitMult(const Vector &u, Vector &v) const;
 
    /** @brief Perform the action of the implicit part of the operator, F:
@@ -445,7 +446,7 @@ public:
 
        Regardless of the choice of F and G, this function should always compute
        @a k = inv(M) g(@a u, t). */
-   virtual void Mult(const Vector &u, Vector &v) const override;
+   virtual void Mult(const Vector &u, Vector &k) const override;
 
    /** @brief Solve for the unknown @a k, at the current time t, the following
        equation:
@@ -496,7 +497,17 @@ public:
        details, see the PETSc Manual. */
    virtual Operator& GetExplicitGradient(const Vector &u) const;
 
-   /** @brief Setup a linear system as needed by some SUNDIALS ODE solvers.
+   /** @brief Setup a linear system as needed by some SUNDIALS ODE solvers to
+       perform a similar action to ImplicitSolve, i.e., solve for k, at the
+       current time t, in F(u + gamma k, k, t) = G(u + gamma k, t).
+
+       The SUNDIALS ODE solvers iteratively solve for k, as knew = kold + dk.
+       The linear system here is for dk, obtained by linearizing the nonlinear
+       system F(u + gamma knew, knew, t) = G(u + gamma knew, t) about dk = 0:
+          F(u + gamma (kold + dk), kold + dk, t) = G(u + gamma (kold + dk), t)
+          => [dF/dk + gamma (dF/du - dG/du)] dk = G - F + O(dk^2)
+       In other words, the linear system to be setup here is A dk = r, where
+       A = [dF/dk + gamma (dF/du - dG/du)] and r = G - F.
 
        For solving an ordinary differential equation of the form
        $ M \frac{dy}{dt} = g(y,t) $, recall that F and G can be defined as one
@@ -506,7 +517,7 @@ public:
          2. F(u,k,t) = M k and G(u,t) = g(u,t)
          3. F(u,k,t) = M k - g(u,t) and G(u,t) = 0
 
-       This function performs setup to solve $ A x = b $ where A is either
+       This function performs setup to solve $ A dk = r $ where A is either
 
          1. A(@a y,t) = I - @a gamma inv(M) J(@a y,t)
          2. A(@a y,t) = M - @a gamma J(@a y,t)
@@ -527,18 +538,26 @@ public:
    virtual int SUNImplicitSetup(const Vector &y, const Vector &v,
                                 int jok, int *jcur, real_t gamma);
 
-   /** @brief Solve the ODE linear system A @a x = @a b, where A is defined by
-       the method SUNImplicitSetup().
+   /** @brief Solve the ODE linear system A @a dk = @a r , where A and r are
+       defined by the method SUNImplicitSetup().
 
-       @param[in]      b   The linear system right-hand side.
-       @param[in,out]  x   On input, the initial guess. On output, the solution.
+       For solving an ordinary differential equation of the form
+       $ M \frac{dy}{dt} = g(y,t) $, recall that F and G can be defined as one
+       of the following:
+
+         1. F(u,k,t) = k and G(u,t) = inv(M) g(u,t)
+         2. F(u,k,t) = M k and G(u,t) = g(u,t)
+         3. F(u,k,t) = M k - g(u,t) and G(u,t) = 0
+
+       @param[in]      r   inv(M) g(y,t) - k for 1 or g(y,t) - M k for 2 & 3.
+       @param[in,out]  dk  On input, the initial guess. On output, the solution.
        @param[in]      tol Linear solve tolerance.
 
        If not re-implemented, this method simply generates an error.
 
        Presently, this method is used by SUNDIALS ODE solvers, for more
        details, see the SUNDIALS User Guides. */
-   virtual int SUNImplicitSolve(const Vector &b, Vector &x, real_t tol);
+   virtual int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol);
 
    /** @brief Setup the mass matrix in the ODE system
        $ M \frac{dy}{dt} = g(y,t) $ .

linalg/sparsemat.cpp

@@ -1267,24 +1267,32 @@ real_t SparseMatrix::InnerProduct(const Vector &x, const Vector &y) const
 
 void SparseMatrix::GetRowSums(Vector &x) const
 {
-   for (int i = 0; i < height; i++)
+   if (Finalized())
    {
-      real_t a = 0.0;
-      if (A)
+      auto d_I = ReadI();
+      auto d_A = ReadData();
+      auto d_x = x.Write();
+      mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
       {
-         for (int j = I[i], end = I[i+1]; j < end; j++)
+         real_t sum = 0.0;
+         for (int j = d_I[i], end = d_I[i+1]; j < end; j++)
          {
-            a += A[j];
+            sum += d_A[j];
          }
-      }
-      else
+         d_x[i] = sum;
+      });
+   }
+   else
+   {
+      for (int i = 0; i < height; i++)
       {
+         real_t a = 0.0;
          for (RowNode *np = Rows[i]; np != NULL; np = np->Prev)
          {
             a += np->Value;
          }
+         x(i) = a;
       }
-      x(i) = a;
    }
 }
 
@@ -3300,7 +3308,7 @@ void SparseMatrix::Print(std::ostream & os, int width_) const
 {
    int i, j;
 
-   if (A == NULL)
+   if (A.Empty())
    {
       RowNode *nd;
       for (i = 0; i < height; i++)
@@ -3354,7 +3362,7 @@ void SparseMatrix::PrintMatlab(std::ostream & os) const
    os.setf(ios::scientific);
    std::streamsize old_prec = os.precision(14);
 
-   if (A == NULL)
+   if (A.Empty())
    {
       RowNode *nd;
       for (i = 0; i < height; i++)
@@ -3397,7 +3405,7 @@ void SparseMatrix::PrintMM(std::ostream & os) const
 
    os << height << " " << width << " " << NumNonZeroElems() << '\n';
 
-   if (A == NULL)
+   if (A.Empty())
    {
       RowNode *nd;
       for (i = 0; i < height; i++)

linalg/sparsemat.hpp

@@ -216,7 +216,7 @@ public:
    void ClearCuSparse() { ClearGPUSparse(); }
 
    /// Check if the SparseMatrix is empty.
-   bool Empty() const { return (A == NULL) && (Rows == NULL); }
+   bool Empty() const { return A.Empty() && (Rows == NULL); }
 
    /// Return the array #I.
    inline int *GetI() { return I; }

linalg/sundials.cpp

@@ -1341,38 +1341,59 @@ CVODESSolver::~CVODESSolver()
 // ARKStep interface
 // ---------------------------------------------------------------------------
 
-int ARKStepSolver::RHS1(realtype t, const N_Vector y, N_Vector ydot,
+int ARKStepSolver::RHS1(realtype t, const N_Vector y, N_Vector result,
                         void *user_data)
 {
    // Get data from N_Vectors
    const SundialsNVector mfem_y(y);
-   SundialsNVector mfem_ydot(ydot);
+   SundialsNVector mfem_result(result);
    ARKStepSolver *self = static_cast<ARKStepSolver*>(user_data);
 
-   // Compute f(t, y) in y' = f(t, y) or fe(t, y) in y' = fe(t, y) + fi(t, y)
+   // Compute either f(t, y) in one of
+   //   1. y' = f(t, y)
+   //   2. M y' = f(t, y)
+   // or fe(t, y) in one of
+   //   1. y' = fe(t, y) + fi(t, y)
+   //   2. M y' = fe(t, y) + fi(t, y)
    self->f->SetTime(t);
    if (self->rk_type == IMEX)
    {
       self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_1);
    }
-   self->f->Mult(mfem_y, mfem_ydot);
+   if (self->f->isExplicit()) // ODE is in form 1
+   {
+      self->f->Mult(mfem_y, mfem_result);
+   }
+   else // ODE is in form 2
+   {
+      self->f->ExplicitMult(mfem_y, mfem_result);
+   }
 
    // Return success
    return (0);
 }
 
-int ARKStepSolver::RHS2(realtype t, const N_Vector y, N_Vector ydot,
+int ARKStepSolver::RHS2(realtype t, const N_Vector y, N_Vector result,
                         void *user_data)
 {
    // Get data from N_Vectors
    const SundialsNVector mfem_y(y);
-   SundialsNVector mfem_ydot(ydot);
+   SundialsNVector mfem_result(result);
    ARKStepSolver *self = static_cast<ARKStepSolver*>(user_data);
 
-   // Compute fi(t, y) in y' = fe(t, y) + fi(t, y)
+   // Compute fi(t, y) in one of
+   //   1. y' = fe(t, y) + fi(t, y)       (ODE is expressed in EXPLICIT form)
+   //   2. M y' = fe(t, y) + fi(y, t)     (ODE is expressed in IMPLICIT form)
    self->f->SetTime(t);
    self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_2);
-   self->f->Mult(mfem_y, mfem_ydot);
+   if (self->f->isExplicit())
+   {
+      self->f->Mult(mfem_y, mfem_result);
+   }
+   else
+   {
+      self->f->ExplicitMult(mfem_y, mfem_result);
+   }
 
    // Return success
    return (0);
@@ -1567,7 +1588,7 @@ void ARKStepSolver::Init(TimeDependentOperator &f_)
    reinit = true;
 }
 
-void ARKStepSolver::Step(Vector &x, double &t, double &dt)
+void ARKStepSolver::Step(Vector &x, real_t &t, real_t &dt)
 {
    Y->MakeRef(x, 0, x.Size());
    MFEM_VERIFY(Y->Size() == x.Size(), "size mismatch");
@@ -1666,7 +1687,7 @@ void ARKStepSolver::UseMFEMMassLinearSolver(int tdep)
    LSM->content      = this;
    LSM->ops->gettype = LSGetType;
    LSM->ops->solve   = ARKStepSolver::MassSysSolve;
-   LSA->ops->free    = LSFree;
+   LSM->ops->free    = LSFree;
 
    M = SUNMatNewEmpty(Sundials::GetContext());
    MFEM_VERIFY(M, "error in SUNMatNewEmpty()");
@@ -1683,6 +1704,9 @@ void ARKStepSolver::UseMFEMMassLinearSolver(int tdep)
    // Set the linear system function
    flag = ARKStepSetMassFn(sundials_mem, ARKStepSolver::MassSysSetup);
    MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetMassFn()");
+
+   // Check that the ODE is not expressed in EXPLICIT form
+   MFEM_VERIFY(!f->isExplicit(), "ODE operator is expressed in EXPLICIT form")
 }
 
 void ARKStepSolver::UseSundialsMassLinearSolver(int tdep)
@@ -1703,6 +1727,9 @@ void ARKStepSolver::UseSundialsMassLinearSolver(int tdep)
    flag = ARKStepSetMassTimes(sundials_mem, NULL, ARKStepSolver::MassMult2,
                               this);
    MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetMassTimes()");
+
+   // Check that the ODE is not expressed in EXPLICIT form
+   MFEM_VERIFY(!f->isExplicit(), "ODE operator is expressed in EXPLICIT form")
 }
 
 void ARKStepSolver::SetStepMode(int itask)

linalg/sundials.hpp

@@ -763,7 +763,7 @@ public:
        @note On input, the values of @a t and @a dt are used to compute desired
        output time for the integration, tout = @a t + @a dt.
    */
-   virtual void Step(Vector &x, double &t, double &dt);
+   virtual void Step(Vector &x, real_t &t, real_t &dt) override;
 
    /** @brief Attach the linear system setup and solve methods from the
        TimeDependentOperator i.e., SUNImplicitSetup() and SUNImplicitSolve() to

mesh/CMakeLists.txt

@@ -12,6 +12,7 @@
 set(SRCS
   attribute_sets.cpp
   element.cpp
+  exodus_writer.cpp
   face_nbr_geom.cpp
   gmsh.cpp
   hexahedron.cpp

mesh/mesh.hpp

@@ -2334,6 +2334,11 @@ public:
                     bool high_order_output=false,
                     int compression_level=0);
 
+#ifdef MFEM_USE_NETCDF
+   /// @brief Export a mesh to an Exodus II file.
+   void PrintExodusII(const std::string fpath);
+#endif
+
    /** @brief Prints the mesh with boundary elements given by the boundary of
        the subdomains, so that the boundary of subdomain i has boundary
        attribute i+1. */

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)

miniapps/meshing/mesh-explorer.cpp

@@ -386,6 +386,9 @@ int main (int argc, char *argv[])
            "S) Save in MFEM serial format\n"
            "T) Save in MFEM parallel format using the current partitioning\n"
            "V) Save in VTK format (only linear and quadratic meshes)\n"
+#ifdef MFEM_USE_NETCDF
+           "X) Save in Exodus II format (only linear and quadratic meshes)\n"
+#endif
            "D) Save as a DataCollection\n"
            "q) Quit\n"
 #ifdef MFEM_USE_ZLIB
@@ -1288,6 +1291,15 @@ int main (int argc, char *argv[])
          cout << "New VTK mesh file: " << omesh_file << endl;
       }
 
+#ifdef MFEM_USE_NETCDF
+      if (mk == 'X')
+      {
+         const char omesh_file[] = "mesh-explorer.e";
+         mesh->PrintExodusII(omesh_file);
+         cout << "New Exodus II mesh file: " << omesh_file << endl;
+      }
+#endif
+
       if (mk == 'D')
       {
          cout << "What type of DataCollection?\n"

tests/unit/CMakeLists.txt

@@ -56,6 +56,7 @@ set(UNIT_TESTS_SRCS
   mesh/test_submesh.cpp
   mesh/test_vtu.cpp
   mesh/test_nurbs.cpp
+  mesh/test_exodus_writer.cpp
   fem/test_1d_bilininteg.cpp
   fem/test_2d_bilininteg.cpp
   fem/test_3d_bilininteg.cpp
@@ -85,6 +86,7 @@ set(UNIT_TESTS_SRCS
   fem/test_get_value.cpp
   fem/test_getderivative.cpp
   fem/test_getgradient.cpp
+  fem/test_gslib.cpp
   fem/test_intrules.cpp
   fem/test_intruletypes.cpp
   fem/test_inversetransform.cpp

tests/unit/fem/test_gslib.cpp

@@ -369,6 +369,124 @@ TEST_CASE("GSLIBInterpolateL2ElementBoundary",
 }
 
 #ifdef MFEM_USE_MPI
+// Custom interpolation procedure with gslib
+TEST_CASE("GSLIBCustomInterpolation",
+          "[GSLIBCustomInterpolation][Parallel][GSLIB]")
+{
+   int myid;
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   int dim      = GENERATE(2, 3);
+   bool simplex = GENERATE(true, false);
+
+   CAPTURE(dim, simplex);
+
+   int nex = 4;
+   int mesh_order = 2;
+   Mesh mesh;
+   if (dim == 2)
+   {
+      Element::Type type = simplex ? Element::TRIANGLE : Element::QUADRILATERAL;
+      mesh = Mesh::MakeCartesian2D(nex, nex, type);
+   }
+   else
+   {
+      Element::Type type = simplex ? Element::TETRAHEDRON : Element::HEXAHEDRON;
+      mesh = Mesh::MakeCartesian3D(nex, nex, nex, type);
+   }
+
+   mesh.SetCurvature(mesh_order);
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+
+   // f(x,y,z) = x^2 + y^2 + z^2
+   auto func = [](const Vector &x)
+   {
+      const int dim = x.Size();
+      double res = 0.0;
+      for (int d = 0; d < dim; d++) { res += std::pow(x(d), 2); }
+      return res;
+   };
+
+   // \nabla f(x,y,z) = [2*x,2*y,2*z]
+   auto func_grad = [](const Vector &x, Vector &p)
+   {
+      const int dim = x.Size();
+      p.SetSize(dim);
+      for (int d = 0; d < dim; d++) { p(d) = 2.0*x(d); }
+   };
+
+   // Set GridFunction to be interpolated
+   int func_order = 3;
+   H1_FECollection c_fec(func_order, dim);
+   FiniteElementSpace c_fespace(&pmesh, &c_fec, 1);
+   GridFunction field_vals(&c_fespace);
+
+   FunctionCoefficient f(func);
+   field_vals.ProjectCoefficient(f);
+
+   // Generate randomized points in [0, 1]^D. Assume ordering by VDIM.
+   int npt = 101;
+   Vector xyz(npt*dim);
+   xyz.Randomize(myid + 1);
+
+   // Find points on the ParMesh
+   Vector interp_vals(npt);
+   FindPointsGSLIB finder;
+   finder.Setup(pmesh);
+   finder.FindPoints(xyz, Ordering::byVDIM);
+
+   /** Interpolate gradient using custom interpolation procedure. */
+   // We first send information to MPI ranks that own the element corresponding
+   // to each point.
+   Array<unsigned int> recv_elem, recv_code;
+   Vector recv_rst;
+   finder.DistributePointInfoToOwningMPIRanks(recv_elem, recv_rst, recv_code);
+   int npt_recv = recv_elem.Size();
+   // Compute gradient locally
+   Vector grad(npt_recv*dim);
+   for (int i = 0; i < npt_recv; i++)
+   {
+      const int e = recv_elem[i];
+
+      IntegrationPoint ip;
+      if (dim == 2)
+      {
+         ip.Set2(recv_rst(dim*i + 0),recv_rst(dim*i + 1));
+      }
+      else
+      {
+         ip.Set3(recv_rst(dim*i + 0),recv_rst(dim*i + 1),
+                 recv_rst(dim*i + 2));
+      }
+      ElementTransformation *Tr = c_fespace.GetElementTransformation(e);
+      Tr->SetIntPoint(&ip);
+
+      Vector gradloc(grad.GetData()+i*dim,dim);
+      field_vals.GetGradient(*Tr, gradloc);
+   }
+
+   // Send the computed gradient back to the ranks that requested it.
+   Vector recv_grad;
+   finder.DistributeInterpolatedValues(grad, dim, Ordering::byVDIM, recv_grad);
+
+   // Check if the received gradient matched analytic gradient.
+   for (int i = 0; i < npt && myid == 0; i++)
+   {
+      Vector x(xyz.GetData()+i*dim,dim);
+      Vector grad_exact(dim);
+      func_grad(x, grad_exact);
+
+      Vector recv_grad_i(recv_grad.GetData()+i*dim,dim);
+
+      for (int d = 0; d < dim; d++)
+      {
+         REQUIRE(grad_exact(d) == Approx(recv_grad(i*dim + d)));
+      }
+   }
+
+   finder.FreeData();
+}
+
 TEST_CASE("GSLIBGSOP", "[GSLIBGSOP][Parallel][GSLIB]")
 {
    int myid;
@@ -434,7 +552,7 @@ TEST_CASE("GSLIBGSOP", "[GSLIBGSOP][Parallel][GSLIB]")
       REQUIRE(vals(i) < 0);
    }
 }
-#endif
+#endif // MFEM_USE_MPI
 
 } //namespace_gslib
 #endif

tests/unit/general/test_umpire_mem.cpp

@@ -18,12 +18,13 @@
 #include <unistd.h>
 #include <stdio.h>
 #include "umpire/Umpire.hpp"
+#include <umpire/strategy/QuickPool.hpp>
 
 #ifdef MFEM_USE_CUDA
 #include <cuda.h>
 constexpr const char * device_name = "cuda";
 #elif defined(MFEM_USE_HIP)
-constexpr const char * device_name = "raja-hip";
+constexpr const char * device_name = "hip";
 #endif
 
 using namespace mfem;
@@ -45,10 +46,12 @@ static bool is_pinned_host(void * h_p)
    unsigned flags;
 #ifdef MFEM_USE_CUDA
    auto err = cudaHostGetFlags(&flags, h_p);
+   cudaGetLastError(); // also resets last error
    if (err == cudaSuccess) { return true; }
    else if (err == cudaErrorInvalidValue) { return false; }
 #elif defined(MFEM_USE_HIP)
    auto err = hipHostGetFlags(&flags, h_p);
+   hipGetLastError(); // also resets last error
    if (err == hipSuccess) { return true; }
    else if (err == hipErrorInvalidValue) { return false; }
 #endif

tests/unit/mesh/test_exodus_writer.cpp

@@ -0,0 +1,129 @@
+// 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 "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+#ifdef MFEM_USE_NETCDF
+static void CompareMeshes(Mesh &mesh1, Mesh &mesh2)
+{
+   REQUIRE(mesh1.GetNE() == mesh2.GetNE());
+   REQUIRE(mesh1.GetNV() == mesh2.GetNV());
+   REQUIRE(mesh1.GetNBE() == mesh2.GetNBE());
+   REQUIRE(mesh1.GetNFaces() == mesh2.GetNFaces());
+
+   const FiniteElementSpace *fespace1 = mesh1.GetNodalFESpace();
+   const FiniteElementSpace *fespace2 = mesh2.GetNodalFESpace();
+
+   // Check elements.
+   Array<int> element_faces1, element_faces2;
+   Array<int> element_orient1, element_orient2;
+   Array<int> dofs1, dofs2;
+
+   for (int ielement = 0; ielement < mesh1.GetNE(); ielement++)
+   {
+      int attr1 = mesh1.GetAttribute(ielement);
+      int attr2 = mesh2.GetAttribute(ielement);
+
+      REQUIRE(attr1 == attr2);
+
+      Element::Type type1 = mesh1.GetElementType(ielement);
+      Element::Type type2 = mesh2.GetElementType(ielement);
+
+      REQUIRE(type1 == type2);
+
+      mesh1.GetElementFaces(ielement, element_faces1, element_orient1);
+      mesh2.GetElementFaces(ielement, element_faces2, element_orient2);
+
+      REQUIRE(element_faces1 == element_faces2);
+      REQUIRE(element_orient1 == element_orient2);
+
+      if (fespace1 && fespace2)
+      {
+         fespace1->GetElementDofs(ielement, dofs1);
+         fespace2->GetElementDofs(ielement, dofs2);
+      }
+      else
+      {
+         mesh1.GetElementVertices(ielement, dofs1);
+         mesh2.GetElementVertices(ielement, dofs2);
+      }
+
+      REQUIRE(dofs1 == dofs2);
+   }
+
+   // Check bdr elements.
+   for (int ibdr_element = 0; ibdr_element < mesh1.GetNBE(); ibdr_element++)
+   {
+      int attr1 = mesh1.GetBdrAttribute(ibdr_element);
+      int attr2 = mesh2.GetBdrAttribute(ibdr_element);
+
+      REQUIRE(attr1 == attr2);
+
+      Element::Type type1 = mesh1.GetBdrElementType(ibdr_element);
+      Element::Type type2 = mesh2.GetBdrElementType(ibdr_element);
+
+      REQUIRE(type1 == type2);
+
+      int face_index1 = mesh1.GetBdrElementFaceIndex(ibdr_element);
+      int face_index2 = mesh2.GetBdrElementFaceIndex(ibdr_element);
+
+      REQUIRE(face_index1 == face_index2);
+   }
+
+   // Check face vertices.
+   Array<int> face_vertices1, face_vertices2;
+   for (int iface_index = 0; iface_index < mesh1.GetNFaces(); iface_index++)
+   {
+      mesh1.GetFaceVertices(iface_index, face_vertices1);
+      mesh2.GetFaceVertices(iface_index, face_vertices2);
+
+      REQUIRE(face_vertices1 == face_vertices2);
+   }
+}
+#endif
+
+TEST_CASE("ExodusII Writer", "[Mesh][ExodusII][MFEMData]")
+{
+#ifdef MFEM_USE_NETCDF
+   // NB: wedge, pyramid and mixed mesh tests require the ExodusII reader PR
+   // to be merged. Pyramid14 tests require the pyramid-dev branch to be merged.
+   auto filename = GENERATE("simple-cube-hex8.e",
+                            "simple-cube-hex27.e",
+                            "simple-cube-tet4.e",
+                            "simple-cube-tet10.e"//,
+                            //  "simple-cube-wedge6.e",
+                            //  "simple-cube-wedge18.e",
+                            //  "simple-cube-pyramid5.e",
+                            //  "simple-cube-pyramid14.e",
+                            //  "simple-cube-multi-element-order1.e",
+                            //  "simple-cube-multi-element-order2.e"
+                           );
+
+   // Load Exodus II mesh from file. NB: do NOT refine as this changes vertex ordering!
+   Mesh original_mesh = Mesh::LoadFromFile(mfem_data_dir + "/exodusii/" + filename,
+                                           0, 0, true);
+
+   // Write generated Exodus II mesh to file.
+   std::string filename_generated = "generated-mesh.e";
+   original_mesh.PrintExodusII(filename_generated);
+
+   // Load generated Exodus II mesh.
+   Mesh generated_mesh = Mesh::LoadFromFile(filename_generated, 0, 0, true);
+
+   CompareMeshes(original_mesh, generated_mesh);
+
+   // Remove temporary file.
+   REQUIRE(remove(filename_generated.c_str()) == 0);
+#endif
+}