diff --git a/CHANGELOG.md b/CHANGELOG.md
index 996c8f04e0..953933b369 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -45,6 +45,10 @@ Specifically, directories set with `CeedAddJitSourceRoot(ceed, "foo/bar")` will
 - OCCA backends were retired.
 - Use clang-tidy to automatically fix if-else blocks which are missing braces.
 
+### Backends
+
+- CUDA: Change the maximum number of points from a compile-time parameter to a run-time parameter to reduce the need for recompilation of at-points kernels.
+
 (v0-12)=
 
 ## v0.12 (Oct 31, 2023)
diff --git a/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp b/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp
index 950949b2ee..2cab29ce47 100644
--- a/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp
+++ b/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp
@@ -731,7 +731,7 @@ static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedO
 // QFunction
 //------------------------------------------------------------------------------
 static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, Tab &tab, CeedInt max_dim,
-                                                     CeedInt max_num_points, CeedInt num_input_fields, CeedOperatorField *op_input_fields,
+                                                     CeedInt num_input_fields, CeedOperatorField *op_input_fields,
                                                      CeedQFunctionField *qf_input_fields, CeedInt num_output_fields,
                                                      CeedOperatorField *op_output_fields, CeedQFunctionField *qf_output_fields,
                                                      std::string qfunction_name, CeedInt Q_1d, bool is_all_tensor, bool is_at_points,
@@ -806,7 +806,7 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C
 
     code << tab << "// -- Coordinates\n";
     code << tab << "CeedScalar r_x[max_dim];\n";
-    code << tab << "ReadPoint<max_dim, coords_comp_stride, max_num_points>(data, elem, p, max_num_points, points.indices, points.coords, r_x);\n\n";
+    code << tab << "ReadPoint<max_dim>(data, elem, p, max_num_points, max_num_points, coords_comp_stride, points.indices, points.coords, r_x);\n\n";
 
     code << tab << "// -- Input fields\n";
     for (CeedInt i = 0; i < num_input_fields; i++) {
@@ -822,18 +822,18 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C
       switch (eval_mode) {
         case CEED_EVAL_NONE:
           code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
-          code << tab << "ReadPoint<num_comp" << var_suffix << ", comp_stride" << var_suffix
-               << ", max_num_points>(data, elem, p, max_num_points, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n";
+          code << tab << "ReadPoint<num_comp" << var_suffix << ">(data, elem, p, max_num_points, max_num_points, comp_stride" << var_suffix
+               << ", indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n";
           break;
         case CEED_EVAL_INTERP:
           code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
-          code << tab << "InterpAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name
-               << ">(data, i, r_c" << var_suffix << ", r_x, r_s" << var_suffix << ");\n";
+          code << tab << "InterpAtPoints" << max_dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ">(data, i, r_c"
+               << var_suffix << ", r_x, r_s" << var_suffix << ");\n";
           break;
         case CEED_EVAL_GRAD:
           code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n";
-          code << tab << "GradAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name
-               << ">(data, i, r_c" << var_suffix << ", r_x, r_s" << var_suffix << ");\n";
+          code << tab << "GradAtPoints" << max_dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ">(data, i, r_c" << var_suffix
+               << ", r_x, r_s" << var_suffix << ");\n";
           break;
         case CEED_EVAL_WEIGHT:
           code << tab << "CeedScalar r_s" << var_suffix << "[1];\n";
@@ -1068,8 +1068,8 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C
           code << tab << "{\n";
           tab.push();
           code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n\n";
-          code << tab << "WritePoint<num_comp" << var_suffix << ", comp_stride" << var_suffix
-               << ", max_num_points>(data, elem, i, points.num_per_elem[elem], indices.outputs[" << i << "]"
+          code << tab << "WritePoint<num_comp" << var_suffix << ">(data, elem, i, max_num_points, points.num_per_elem[elem], comp_stride"
+               << var_suffix << ", indices.outputs[" << i << "]"
                << ", r_s" << var_suffix << ", d" << var_suffix << ");\n";
           tab.pop();
           code << tab << "}\n";
@@ -1081,8 +1081,8 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C
           code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << "; j++) r_s" << var_suffix << "[j] = 0.0;\n";
           tab.pop();
           code << tab << "}\n";
-          code << tab << "InterpTransposeAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name
-               << ">(data, i, r_s" << var_suffix << ", r_x, r_c" << var_suffix << ");\n";
+          code << tab << "InterpTransposeAtPoints" << max_dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name
+               << ">(data, i, max_num_points, r_s" << var_suffix << ", r_x, r_c" << var_suffix << ");\n";
           break;
         case CEED_EVAL_GRAD:
           code << tab << "if (i >= points.num_per_elem[elem]) {\n";
@@ -1090,8 +1090,8 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C
           code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << "*dim" << var_suffix << "; j++) r_s" << var_suffix << "[j] = 0.0;\n";
           tab.pop();
           code << tab << "}\n";
-          code << tab << "GradTransposeAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name
-               << ">(data, i, r_s" << var_suffix << ", r_x, r_c" << var_suffix << ");\n";
+          code << tab << "GradTransposeAtPoints" << max_dim << "d<num_comp" << var_suffix << ", " << P_name << ", " << Q_name
+               << ">(data, i, max_num_points, r_s" << var_suffix << ", r_x, r_c" << var_suffix << ");\n";
           break;
           // LCOV_EXCL_START
         case CEED_EVAL_WEIGHT:
@@ -1158,7 +1158,7 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C
 extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op, bool *is_good_build) {
   bool                    is_all_tensor = true, is_all_nontensor = true, is_at_points = false, use_3d_slices = false;
   Ceed                    ceed;
-  CeedInt                 Q = 0, Q_1d = 0, num_input_fields, num_output_fields, max_dim = 1, max_num_points = 0, coords_comp_stride = 0;
+  CeedInt                 Q = 0, Q_1d = 0, num_input_fields, num_output_fields, max_dim = 1;
   CeedQFunctionField     *qf_input_fields, *qf_output_fields;
   CeedQFunction_Cuda_gen *qf_data;
   CeedQFunction           qf;
@@ -1263,21 +1263,25 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op, bool *is_good_b
     CeedElemRestriction       rstr_points = NULL;
 
     CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL));
-    CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points));
-    CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &coords_comp_stride));
     CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &coords_dim));
     CeedCallBackend(CeedElemRestrictionGetData(rstr_points, &rstr_data));
     data->points.indices = (CeedInt *)rstr_data->d_offsets;
-    CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
     if (max_dim == 0) max_dim = coords_dim;
-    if (Q_1d == 0) max_num_points = ceil(pow(max_num_points, 1.0 / max_dim));
+    if (Q_1d == 0) {
+      CeedInt elem_size, num_comp;
+
+      CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_points, &elem_size));
+      CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &num_comp));
+      Q_1d = CeedIntPow(elem_size / num_comp, 1 / max_dim);
+    }
+    CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
   }
   if (max_dim == 0) max_dim = 1;
   data->dim = max_dim;
   if (is_at_points) use_3d_slices = false;
   if (Q_1d == 0) {
     if (is_at_points) {
-      Q_1d = max_num_points;
+      Q_1d = 1;
     } else {
       CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q_1d));
     }
@@ -1398,8 +1402,8 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op, bool *is_good_b
     code << tab << "const CeedInt Q_1d = " << Q_1d << ";\n";
   }
   if (is_at_points) {
-    code << tab << "const CeedInt max_num_points = " << max_num_points << ";\n";
-    code << tab << "const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n";
+    code << tab << "const CeedInt max_num_points = points.max_num_points;\n";
+    code << tab << "const CeedInt coords_comp_stride = points.coords_comp_stride;\n";
   }
 
   // Shared data
@@ -1630,9 +1634,9 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op, bool *is_good_b
   }
 
   // -- Q function
-  CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, max_num_points, num_input_fields, op_input_fields,
-                                                            qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, qfunction_name,
-                                                            Q_1d, is_all_tensor, is_at_points, use_3d_slices, false));
+  CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, num_input_fields, op_input_fields, qf_input_fields,
+                                                            num_output_fields, op_output_fields, qf_output_fields, qfunction_name, Q_1d,
+                                                            is_all_tensor, is_at_points, use_3d_slices, false));
 
   // -- Output basis and restriction
   code << "\n" << tab << "// -- Output field basis action and restrictions\n";
@@ -1686,7 +1690,7 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op, bool *is_good_b
 static int CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(CeedOperator op, bool is_full, bool *is_good_build) {
   bool                    is_all_tensor = true, is_at_points = false, use_3d_slices = false;
   Ceed                    ceed;
-  CeedInt                 Q, Q_1d, num_input_fields, num_output_fields, max_dim = 1, max_num_points = 0, coords_comp_stride = 0;
+  CeedInt                 Q, Q_1d, num_input_fields, num_output_fields, max_dim = 1;
   CeedQFunctionField     *qf_input_fields, *qf_output_fields;
   CeedQFunction_Cuda_gen *qf_data;
   CeedQFunction           qf;
@@ -1705,14 +1709,6 @@ static int CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(CeedOperator op, boo
   Q       = data->Q;
   Q_1d    = data->Q_1d;
   max_dim = data->dim;
-  {
-    CeedElemRestriction rstr_points = NULL;
-
-    CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL));
-    CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points));
-    CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &coords_comp_stride));
-    CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
-  }
   CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
   CeedCallBackend(CeedQFunctionGetData(qf, &qf_data));
   CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
@@ -1797,8 +1793,8 @@ static int CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(CeedOperator op, boo
 
   code << tab << "const CeedInt max_dim = " << max_dim << ";\n";
   code << tab << "const CeedInt Q_1d = " << Q_1d << ";\n";
-  code << tab << "const CeedInt max_num_points = " << max_num_points << ";\n";
-  code << tab << "const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n";
+  code << tab << "const CeedInt max_num_points = points.max_num_points;\n";
+  code << tab << "const CeedInt coords_comp_stride = points.coords_comp_stride;\n";
 
   // Shared data
   code << tab << "extern __shared__ CeedScalar slice[];\n";
@@ -2056,9 +2052,9 @@ static int CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(CeedOperator op, boo
   }
 
   // -- Q function
-  CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, max_num_points, num_input_fields, op_input_fields,
-                                                            qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, qfunction_name,
-                                                            Q_1d, is_all_tensor, is_at_points, use_3d_slices, true));
+  CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, num_input_fields, op_input_fields, qf_input_fields,
+                                                            num_output_fields, op_output_fields, qf_output_fields, qfunction_name, Q_1d,
+                                                            is_all_tensor, is_at_points, use_3d_slices, true));
 
   // -- Output basis and restriction
   code << "\n" << tab << "// -- Output field basis action and restrictions\n";
@@ -2175,7 +2171,7 @@ extern "C" int CeedOperatorBuildKernelFullAssemblyAtPoints_Cuda_gen(CeedOperator
 extern "C" int CeedOperatorBuildKernelLinearAssembleQFunction_Cuda_gen(CeedOperator op, bool *is_good_build) {
   bool                    is_all_tensor = true, is_all_nontensor = true, is_at_points = false, use_3d_slices = false;
   Ceed                    ceed;
-  CeedInt                 Q, Q_1d, num_input_fields, num_output_fields, max_dim = 1, max_num_points = 0;
+  CeedInt                 Q, Q_1d, num_input_fields, num_output_fields, max_dim = 1;
   CeedQFunctionField     *qf_input_fields, *qf_output_fields;
   CeedQFunction_Cuda_gen *qf_data;
   CeedQFunction           qf;
@@ -2673,9 +2669,9 @@ extern "C" int CeedOperatorBuildKernelLinearAssembleQFunction_Cuda_gen(CeedOpera
   }
 
   // -- Q function
-  CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, max_num_points, num_input_fields, op_input_fields,
-                                                            qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, qfunction_name,
-                                                            Q_1d, is_all_tensor, is_at_points, use_3d_slices, true));
+  CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, num_input_fields, op_input_fields, qf_input_fields,
+                                                            num_output_fields, op_output_fields, qf_output_fields, qfunction_name, Q_1d,
+                                                            is_all_tensor, is_at_points, use_3d_slices, true));
 
   // -- Output basis and restriction
   code << "\n" << tab << "// -- Output field basis action and restrictions\n";
diff --git a/backends/cuda-gen/ceed-cuda-gen-operator.c b/backends/cuda-gen/ceed-cuda-gen-operator.c
index f8228fbfbd..a28bcff3f4 100644
--- a/backends/cuda-gen/ceed-cuda-gen-operator.c
+++ b/backends/cuda-gen/ceed-cuda-gen-operator.c
@@ -174,17 +174,20 @@ static int CeedOperatorApplyAddCore_Cuda_gen(CeedOperator op, CUstream stream, c
   CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points));
   if (is_at_points) {
     // Coords
-    CeedVector vec;
+    CeedVector          vec;
+    CeedElemRestriction rstr_points = NULL;
 
-    CeedCallBackend(CeedOperatorAtPointsGetPoints(op, NULL, &vec));
+    CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &vec));
     CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->points.coords));
+    CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &data->points.max_num_points));
+    CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &data->points.coords_comp_stride));
     CeedCallBackend(CeedVectorDestroy(&vec));
+    CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
 
     // Points per elem
     if (num_elem != data->points.num_elem) {
-      CeedInt            *points_per_elem;
-      const CeedInt       num_bytes   = num_elem * sizeof(CeedInt);
-      CeedElemRestriction rstr_points = NULL;
+      CeedInt      *points_per_elem;
+      const CeedInt num_bytes = num_elem * sizeof(CeedInt);
 
       data->points.num_elem = num_elem;
       CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL));
@@ -608,11 +611,15 @@ static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Cuda_gen(CeedOperator o
 
     // Point coordinates
     {
-      CeedVector vec;
+      CeedVector          vec;
+      CeedElemRestriction rstr_points = NULL;
 
-      CeedCallBackend(CeedOperatorAtPointsGetPoints(op, NULL, &vec));
+      CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &vec));
       CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->points.coords));
+      CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &data->points.max_num_points));
+      CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &data->points.coords_comp_stride));
       CeedCallBackend(CeedVectorDestroy(&vec));
+      CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
 
       // Points per elem
       if (num_elem != data->points.num_elem) {
@@ -774,11 +781,15 @@ static int CeedOperatorAssembleSingleAtPoints_Cuda_gen(CeedOperator op, CeedInt
 
     // Point coordinates
     {
-      CeedVector vec;
+      CeedVector          vec;
+      CeedElemRestriction rstr_points = NULL;
 
-      CeedCallBackend(CeedOperatorAtPointsGetPoints(op, NULL, &vec));
+      CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &vec));
       CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->points.coords));
+      CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &data->points.max_num_points));
+      CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &data->points.coords_comp_stride));
       CeedCallBackend(CeedVectorDestroy(&vec));
+      CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
 
       // Points per elem
       if (num_elem != data->points.num_elem) {
diff --git a/backends/cuda-ref/ceed-cuda-ref-basis.c b/backends/cuda-ref/ceed-cuda-ref-basis.c
index 684f329e89..555c15bce9 100644
--- a/backends/cuda-ref/ceed-cuda-ref-basis.c
+++ b/backends/cuda-ref/ceed-cuda-ref-basis.c
@@ -105,7 +105,7 @@ static int CeedBasisApplyAdd_Cuda(CeedBasis basis, const CeedInt num_elem, CeedT
 static int CeedBasisApplyAtPointsCore_Cuda(CeedBasis basis, bool apply_add, const CeedInt num_elem, const CeedInt *num_points,
                                            CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) {
   Ceed              ceed;
-  CeedInt           Q_1d, dim, max_num_points = num_points[0];
+  CeedInt           P_1d, Q_1d, dim, max_num_points = num_points[0];
   const CeedInt     is_transpose   = t_mode == CEED_TRANSPOSE;
   const int         max_block_size = 32;
   const CeedScalar *d_x, *d_u;
@@ -114,6 +114,7 @@ static int CeedBasisApplyAtPointsCore_Cuda(CeedBasis basis, bool apply_add, cons
 
   CeedCallBackend(CeedBasisGetData(basis, &data));
   CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
+  CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
   CeedCallBackend(CeedBasisGetDimension(basis, &dim));
 
   // Weight handled separately
@@ -126,6 +127,7 @@ static int CeedBasisApplyAtPointsCore_Cuda(CeedBasis basis, bool apply_add, cons
 
   // Check padded to uniform number of points per elem
   for (CeedInt i = 1; i < num_elem; i++) max_num_points = CeedIntMax(max_num_points, num_points[i]);
+  data->num_points = max_num_points;
   {
     CeedInt  num_comp, q_comp;
     CeedSize len, len_required;
@@ -158,36 +160,30 @@ static int CeedBasisApplyAtPointsCore_Cuda(CeedBasis basis, bool apply_add, cons
     }
   }
 
-  // Build kernels if needed
-  if (data->num_points != max_num_points) {
-    CeedInt P_1d;
-
-    CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
-    data->num_points = max_num_points;
-
-    // -- Create interp matrix to Chebyshev coefficients
-    if (!data->d_chebyshev_interp_1d) {
-      CeedSize    interp_bytes;
-      CeedScalar *chebyshev_interp_1d;
-
-      interp_bytes = P_1d * Q_1d * sizeof(CeedScalar);
-      CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d));
-      CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d));
-      CeedCallCuda(ceed, cudaMalloc((void **)&data->d_chebyshev_interp_1d, interp_bytes));
-      CeedCallCuda(ceed, cudaMemcpy(data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice));
-      CeedCallBackend(CeedFree(&chebyshev_interp_1d));
-    }
+  // Create interp matrix to Chebyshev coefficients
+  if (!data->d_chebyshev_interp_1d) {
+    CeedSize    interp_bytes;
+    CeedScalar *chebyshev_interp_1d;
+
+    interp_bytes = P_1d * Q_1d * sizeof(CeedScalar);
+    CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d));
+    CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d));
+    CeedCallCuda(ceed, cudaMalloc((void **)&data->d_chebyshev_interp_1d, interp_bytes));
+    CeedCallCuda(ceed, cudaMemcpy(data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice));
+    CeedCallBackend(CeedFree(&chebyshev_interp_1d));
+  }
 
+  // Build kernels if needed
+  if (!data->moduleAtPoints) {
     // -- Compile kernels
     const char basis_kernel_source[] = "// AtPoints basis source\n#include <ceed/jit-source/cuda/cuda-ref-basis-tensor-at-points.h>\n";
     CeedInt    num_comp;
 
-    if (data->moduleAtPoints) CeedCallCuda(ceed, cuModuleUnload(data->moduleAtPoints));
     CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
-    CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, "basis_at_points", &data->moduleAtPoints, 9, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d,
+    CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, "basis_at_points", &data->moduleAtPoints, 8, "BASIS_Q_1D", Q_1d, "BASIS_P_1D", P_1d,
                                      "BASIS_BUF_LEN", Q_1d * CeedIntPow(Q_1d > P_1d ? Q_1d : P_1d, dim - 1), "BASIS_DIM", dim, "BASIS_NUM_COMP",
-                                     num_comp, "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_NUM_PTS",
-                                     max_num_points, "POINTS_BUFF_LEN", CeedIntPow(Q_1d, dim - 1)));
+                                     num_comp, "BASIS_NUM_NODES", CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "POINTS_BUFF_LEN",
+                                     CeedIntPow(Q_1d, dim - 1)));
     CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "InterpAtPoints", &data->InterpAtPoints));
     CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "InterpTransposeAtPoints", &data->InterpTransposeAtPoints));
     CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "GradAtPoints", &data->GradAtPoints));
@@ -212,15 +208,15 @@ static int CeedBasisApplyAtPointsCore_Cuda(CeedBasis basis, bool apply_add, cons
   // Basis action
   switch (eval_mode) {
     case CEED_EVAL_INTERP: {
-      void         *interp_args[] = {(void *)&num_elem, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v};
-      const CeedInt block_size    = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size);
+      void *interp_args[] = {(void *)&num_elem, (void *)&data->num_points, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v};
+      const CeedInt block_size = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size);
 
       CeedCallBackend(CeedRunKernel_Cuda(ceed, is_transpose ? data->InterpTransposeAtPoints : data->InterpAtPoints, num_elem, block_size,
                                          interp_args));
     } break;
     case CEED_EVAL_GRAD: {
-      void         *grad_args[] = {(void *)&num_elem, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v};
-      const CeedInt block_size  = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size);
+      void *grad_args[] = {(void *)&num_elem, (void *)&data->num_points, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v};
+      const CeedInt block_size = CeedIntMin(CeedIntPow(Q_1d, dim), max_block_size);
 
       CeedCallBackend(CeedRunKernel_Cuda(ceed, is_transpose ? data->GradTransposeAtPoints : data->GradAtPoints, num_elem, block_size, grad_args));
     } break;
diff --git a/backends/cuda-ref/ceed-cuda-ref-restriction.c b/backends/cuda-ref/ceed-cuda-ref-restriction.c
index ae53a9628f..d94421f94e 100644
--- a/backends/cuda-ref/ceed-cuda-ref-restriction.c
+++ b/backends/cuda-ref/ceed-cuda-ref-restriction.c
@@ -68,14 +68,12 @@ static inline int CeedElemRestrictionSetupCompile_Cuda(CeedElemRestriction rstr)
       CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OffsetTranspose", &impl->ApplyTranspose));
     } break;
     case CEED_RESTRICTION_POINTS: {
-      const char restriction_kernel_source[] =
-          "// AtPoints restriction source\n#include <ceed/jit-source/cuda/cuda-ref-restriction-at-points.h>\n\n"
-          "// Standard restriction source\n#include <ceed/jit-source/cuda/cuda-ref-restriction-offset.h>\n";
+      const char restriction_kernel_source[] = "// AtPoints restriction source\n#include <ceed/jit-source/cuda/cuda-ref-restriction-at-points.h>\n";
 
-      CeedCallBackend(CeedCompile_Cuda(ceed, restriction_kernel_source, "restriction_at_points", &impl->module, 6, "RSTR_ELEM_SIZE", elem_size,
-                                       "RSTR_NUM_ELEM", num_elem, "RSTR_NUM_COMP", num_comp, "RSTR_NUM_NODES", impl->num_nodes, "RSTR_COMP_STRIDE",
-                                       comp_stride, "USE_DETERMINISTIC", is_deterministic ? 1 : 0));
-      CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OffsetNoTranspose", &impl->ApplyNoTranspose));
+      CeedCallBackend(CeedCompile_Cuda(ceed, restriction_kernel_source, "restriction_at_points", &impl->module, 5, "RSTR_NUM_ELEM", num_elem,
+                                       "RSTR_NUM_COMP", num_comp, "RSTR_NUM_NODES", impl->num_nodes, "RSTR_COMP_STRIDE", comp_stride,
+                                       "USE_DETERMINISTIC", is_deterministic ? 1 : 0));
+      CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "AtPointsNoTranspose", &impl->ApplyNoTranspose));
       CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "AtPointsTranspose", &impl->ApplyTranspose));
     } break;
     case CEED_RESTRICTION_ORIENTED: {
@@ -155,7 +153,15 @@ static inline int CeedElemRestrictionApply_Cuda_Core(CeedElemRestriction rstr, C
 
         CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyNoTranspose, grid, block_size, args));
       } break;
-      case CEED_RESTRICTION_POINTS:
+      case CEED_RESTRICTION_POINTS: {
+        CeedInt max_num_points;
+
+        CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr, &max_num_points));
+
+        void *args[] = {(void *)&max_num_points, &impl->d_offsets, &d_u, &d_v};
+
+        CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyNoTranspose, grid, block_size, args));
+      } break;
       case CEED_RESTRICTION_STANDARD: {
         void *args[] = {&impl->d_offsets, &d_u, &d_v};
 
@@ -201,12 +207,17 @@ static inline int CeedElemRestrictionApply_Cuda_Core(CeedElemRestriction rstr, C
         CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args));
       } break;
       case CEED_RESTRICTION_POINTS: {
+        CeedInt max_num_points;
+
+        CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr, &max_num_points));
         if (!is_deterministic) {
-          void *args[] = {&impl->d_offsets, &impl->d_points_per_elem, &d_u, &d_v};
+          void *args[] = {(void *)&max_num_points, &impl->d_offsets, &impl->d_points_per_elem, &d_u, &d_v};
 
           CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args));
         } else {
-          void *args[] = {&impl->d_l_vec_indices, &impl->d_t_indices, &impl->d_points_per_elem, &impl->d_t_offsets, &d_u, &d_v};
+          void *args[] = {
+              (void *)&max_num_points, &impl->d_l_vec_indices, &impl->d_t_indices, &impl->d_points_per_elem, &impl->d_t_offsets, &d_u, &d_v
+          };
 
           CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args));
         }
diff --git a/backends/cuda-shared/ceed-cuda-shared-basis.c b/backends/cuda-shared/ceed-cuda-shared-basis.c
index a3f8d74d19..61305391ed 100644
--- a/backends/cuda-shared/ceed-cuda-shared-basis.c
+++ b/backends/cuda-shared/ceed-cuda-shared-basis.c
@@ -214,13 +214,14 @@ static int CeedBasisApplyAtPointsCore_Cuda_shared(CeedBasis basis, bool apply_ad
                                                   CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) {
   Ceed                   ceed;
   Ceed_Cuda             *ceed_Cuda;
-  CeedInt                Q_1d, dim, num_comp, max_num_points = num_points[0];
+  CeedInt                P_1d, Q_1d, dim, num_comp, max_num_points = num_points[0];
   const CeedInt          is_transpose = t_mode == CEED_TRANSPOSE;
   const CeedScalar      *d_x, *d_u;
   CeedScalar            *d_v;
   CeedBasis_Cuda_shared *data;
 
   CeedCallBackend(CeedBasisGetData(basis, &data));
+  CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
   CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
   CeedCallBackend(CeedBasisGetDimension(basis, &dim));
   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
@@ -236,6 +237,7 @@ static int CeedBasisApplyAtPointsCore_Cuda_shared(CeedBasis basis, bool apply_ad
 
   // Check padded to uniform number of points per elem
   for (CeedInt i = 1; i < num_elem; i++) max_num_points = CeedIntMax(max_num_points, num_points[i]);
+  data->num_points = max_num_points;
   {
     CeedInt  q_comp;
     CeedSize len, len_required;
@@ -266,35 +268,29 @@ static int CeedBasisApplyAtPointsCore_Cuda_shared(CeedBasis basis, bool apply_ad
     }
   }
 
-  // Build kernels if needed
-  if (data->num_points != max_num_points) {
-    CeedInt P_1d;
-
-    CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
-    data->num_points = max_num_points;
-
-    // -- Create interp matrix to Chebyshev coefficients
-    if (!data->d_chebyshev_interp_1d) {
-      CeedSize    interp_bytes;
-      CeedScalar *chebyshev_interp_1d;
-
-      interp_bytes = P_1d * Q_1d * sizeof(CeedScalar);
-      CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d));
-      CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d));
-      CeedCallCuda(ceed, cudaMalloc((void **)&data->d_chebyshev_interp_1d, interp_bytes));
-      CeedCallCuda(ceed, cudaMemcpy(data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice));
-      CeedCallBackend(CeedFree(&chebyshev_interp_1d));
-    }
+  // Create interp matrix to Chebyshev coefficients, if needed
+  if (!data->d_chebyshev_interp_1d) {
+    CeedSize    interp_bytes;
+    CeedScalar *chebyshev_interp_1d;
+
+    interp_bytes = P_1d * Q_1d * sizeof(CeedScalar);
+    CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d));
+    CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d));
+    CeedCallCuda(ceed, cudaMalloc((void **)&data->d_chebyshev_interp_1d, interp_bytes));
+    CeedCallCuda(ceed, cudaMemcpy(data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice));
+    CeedCallBackend(CeedFree(&chebyshev_interp_1d));
+  }
 
+  // Build kernels if needed
+  if (!data->moduleAtPoints) {
     // -- Compile kernels
     const char basis_kernel_source[] = "// AtPoints basis source\n#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points.h>\n";
     CeedInt    num_comp;
 
-    if (data->moduleAtPoints) CeedCallCuda(ceed, cuModuleUnload(data->moduleAtPoints));
     CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
-    CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, "basis_at_points_shared", &data->moduleAtPoints, 8, "BASIS_Q_1D", Q_1d, "BASIS_P_1D",
+    CeedCallBackend(CeedCompile_Cuda(ceed, basis_kernel_source, "basis_at_points_shared", &data->moduleAtPoints, 7, "BASIS_Q_1D", Q_1d, "BASIS_P_1D",
                                      P_1d, "BASIS_T_1D", CeedIntMax(Q_1d, P_1d), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp, "BASIS_NUM_NODES",
-                                     CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim), "BASIS_NUM_PTS", max_num_points));
+                                     CeedIntPow(P_1d, dim), "BASIS_NUM_QPTS", CeedIntPow(Q_1d, dim)));
     CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "InterpAtPoints", &data->InterpAtPoints));
     CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "InterpTransposeAtPoints", &data->InterpTransposeAtPoints));
     CeedCallBackend(CeedGetKernel_Cuda(ceed, data->moduleAtPoints, "InterpTransposeAddAtPoints", &data->InterpTransposeAddAtPoints));
@@ -325,7 +321,7 @@ static int CeedBasisApplyAtPointsCore_Cuda_shared(CeedBasis basis, bool apply_ad
       CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
       CeedInt thread_1d = CeedIntMax(Q_1d, P_1d);
 
-      void *interp_args[] = {(void *)&num_elem, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v};
+      void *interp_args[] = {(void *)&num_elem, (void *)&data->num_points, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v};
 
       if (dim == 1) {
         // avoid >512 total threads
@@ -375,7 +371,7 @@ static int CeedBasisApplyAtPointsCore_Cuda_shared(CeedBasis basis, bool apply_ad
       CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
       CeedInt thread_1d = CeedIntMax(Q_1d, P_1d);
 
-      void *grad_args[] = {(void *)&num_elem, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v};
+      void *grad_args[] = {(void *)&num_elem, (void *)&data->num_points, &data->d_chebyshev_interp_1d, &data->d_points_per_elem, &d_x, &d_u, &d_v};
 
       if (dim == 1) {
         // avoid >512 total threads
diff --git a/include/ceed/jit-source/cuda/cuda-gen-templates.h b/include/ceed/jit-source/cuda/cuda-gen-templates.h
index 5fd998d9e9..c8cf8984d0 100644
--- a/include/ceed/jit-source/cuda/cuda-gen-templates.h
+++ b/include/ceed/jit-source/cuda/cuda-gen-templates.h
@@ -24,27 +24,29 @@ inline __device__ void LoadMatrix(SharedData_Cuda &data, const CeedScalar *__res
 //------------------------------------------------------------------------------
 // L-vector -> single point
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int COMP_STRIDE, int NUM_PTS>
-inline __device__ void ReadPoint(SharedData_Cuda &data, const CeedInt elem, const CeedInt p, const CeedInt points_in_elem,
-                                 const CeedInt *__restrict__ indices, const CeedScalar *__restrict__ d_u, CeedScalar *r_u) {
-  const CeedInt ind = indices[p + elem * NUM_PTS];
+template <int NUM_COMP>
+inline __device__ void ReadPoint(SharedData_Cuda &data, const CeedInt elem, const CeedInt p, const CeedInt max_num_points,
+                                 const CeedInt points_in_elem, const CeedInt     comp_stride, const CeedInt *__restrict__ indices,
+                                 const CeedScalar *__restrict__ d_u, CeedScalar *r_u) {
+  const CeedInt ind = indices[p + elem * max_num_points];
 
   for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
-    r_u[comp] = d_u[ind + comp * COMP_STRIDE];
+    r_u[comp] = d_u[ind + comp * comp_stride];
   }
 }
 
 //------------------------------------------------------------------------------
 // Single point -> L-vector
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int COMP_STRIDE, int NUM_PTS>
-inline __device__ void WritePoint(SharedData_Cuda &data, const CeedInt elem, const CeedInt p, const CeedInt points_in_elem,
-                                  const CeedInt *__restrict__ indices, const CeedScalar *__restrict__ r_u, CeedScalar *d_u) {
+template <int NUM_COMP>
+inline __device__ void WritePoint(SharedData_Cuda &data, const CeedInt elem, const CeedInt p, const CeedInt max_num_points,
+                                  const CeedInt points_in_elem, const CeedInt     comp_stride, const CeedInt *__restrict__ indices,
+                                  const CeedScalar *__restrict__ r_u, CeedScalar *d_u) {
   if (p < points_in_elem) {
-    const CeedInt ind = indices[p + elem * NUM_PTS];
+    const CeedInt ind = indices[p + elem * max_num_points];
 
     for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
-      d_u[ind + comp * COMP_STRIDE] += r_u[comp];
+      d_u[ind + comp * comp_stride] += r_u[comp];
     }
   }
 }
diff --git a/include/ceed/jit-source/cuda/cuda-ref-basis-tensor-at-points.h b/include/ceed/jit-source/cuda/cuda-ref-basis-tensor-at-points.h
index 1ac83beb8d..9f87ecc5ef 100644
--- a/include/ceed/jit-source/cuda/cuda-ref-basis-tensor-at-points.h
+++ b/include/ceed/jit-source/cuda/cuda-ref-basis-tensor-at-points.h
@@ -41,7 +41,7 @@ inline __device__ void ChebyshevDerivativeAtPoint(const CeedScalar x, CeedScalar
 //------------------------------------------------------------------------------
 // Interp
 //------------------------------------------------------------------------------
-extern "C" __global__ void InterpAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ chebyshev_interp_1d,
+extern "C" __global__ void InterpAtPoints(const CeedInt num_elem, const CeedInt max_num_points, const CeedScalar *__restrict__ chebyshev_interp_1d,
                                           const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ coords,
                                           const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v) {
   const CeedInt i = threadIdx.x;
@@ -59,9 +59,9 @@ extern "C" __global__ void InterpAtPoints(const CeedInt num_elem, const CeedScal
   const CeedInt P             = BASIS_P_1D;
   const CeedInt Q             = BASIS_Q_1D;
   const CeedInt u_stride      = BASIS_NUM_NODES;
-  const CeedInt v_stride      = BASIS_NUM_PTS;
+  const CeedInt v_stride      = max_num_points;
   const CeedInt u_comp_stride = num_elem * BASIS_NUM_NODES;
-  const CeedInt v_comp_stride = num_elem * BASIS_NUM_PTS;
+  const CeedInt v_comp_stride = num_elem * max_num_points;
   const CeedInt u_size        = BASIS_NUM_NODES;
 
   // Apply basis element by element
@@ -96,7 +96,7 @@ extern "C" __global__ void InterpAtPoints(const CeedInt num_elem, const CeedScal
 
       // Map to point
       __syncthreads();
-      for (CeedInt p = threadIdx.x; p < BASIS_NUM_PTS; p += blockDim.x) {
+      for (CeedInt p = threadIdx.x; p < max_num_points; p += blockDim.x) {
         pre  = BASIS_NUM_QPTS;
         post = 1;
         for (CeedInt d = 0; d < BASIS_DIM; d++) {
@@ -124,9 +124,10 @@ extern "C" __global__ void InterpAtPoints(const CeedInt num_elem, const CeedScal
   }
 }
 
-extern "C" __global__ void InterpTransposeAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ chebyshev_interp_1d,
-                                                   const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ coords,
-                                                   const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v) {
+extern "C" __global__ void InterpTransposeAtPoints(const CeedInt num_elem, const CeedInt max_num_points,
+                                                   const CeedScalar *__restrict__ chebyshev_interp_1d, const CeedInt *__restrict__ points_per_elem,
+                                                   const CeedScalar *__restrict__ coords, const CeedScalar *__restrict__ u,
+                                                   CeedScalar *__restrict__ v) {
   const CeedInt i = threadIdx.x;
 
   __shared__ CeedScalar s_mem[BASIS_Q_1D * BASIS_P_1D + 2 * BASIS_BUF_LEN + POINTS_BUFF_LEN * BASIS_Q_1D];
@@ -141,11 +142,11 @@ extern "C" __global__ void InterpTransposeAtPoints(const CeedInt num_elem, const
 
   const CeedInt P             = BASIS_P_1D;
   const CeedInt Q             = BASIS_Q_1D;
-  const CeedInt u_stride      = BASIS_NUM_PTS;
+  const CeedInt u_stride      = max_num_points;
   const CeedInt v_stride      = BASIS_NUM_NODES;
-  const CeedInt u_comp_stride = num_elem * BASIS_NUM_PTS;
+  const CeedInt u_comp_stride = num_elem * max_num_points;
   const CeedInt v_comp_stride = num_elem * BASIS_NUM_NODES;
-  const CeedInt u_size        = BASIS_NUM_PTS;
+  const CeedInt u_size        = max_num_points;
 
   // Apply basis element by element
   for (CeedInt elem = blockIdx.x; elem < num_elem; elem += gridDim.x) {
@@ -162,7 +163,7 @@ extern "C" __global__ void InterpTransposeAtPoints(const CeedInt num_elem, const
 
       // Map from point
       __syncthreads();
-      for (CeedInt p = threadIdx.x; p < BASIS_NUM_PTS; p += blockDim.x) {
+      for (CeedInt p = threadIdx.x; p < max_num_points; p += blockDim.x) {
         if (p >= points_per_elem[elem]) continue;
         pre  = 1;
         post = 1;
@@ -223,7 +224,7 @@ extern "C" __global__ void InterpTransposeAtPoints(const CeedInt num_elem, const
 //------------------------------------------------------------------------------
 // Grad
 //------------------------------------------------------------------------------
-extern "C" __global__ void GradAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ chebyshev_interp_1d,
+extern "C" __global__ void GradAtPoints(const CeedInt num_elem, const CeedInt max_num_points, const CeedScalar *__restrict__ chebyshev_interp_1d,
                                         const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ coords,
                                         const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v) {
   const CeedInt i = threadIdx.x;
@@ -241,12 +242,12 @@ extern "C" __global__ void GradAtPoints(const CeedInt num_elem, const CeedScalar
   const CeedInt P             = BASIS_P_1D;
   const CeedInt Q             = BASIS_Q_1D;
   const CeedInt u_stride      = BASIS_NUM_NODES;
-  const CeedInt v_stride      = BASIS_NUM_PTS;
+  const CeedInt v_stride      = max_num_points;
   const CeedInt u_comp_stride = num_elem * BASIS_NUM_NODES;
-  const CeedInt v_comp_stride = num_elem * BASIS_NUM_PTS;
+  const CeedInt v_comp_stride = num_elem * max_num_points;
   const CeedInt u_size        = BASIS_NUM_NODES;
   const CeedInt u_dim_stride  = 0;
-  const CeedInt v_dim_stride  = num_elem * BASIS_NUM_PTS * BASIS_NUM_COMP;
+  const CeedInt v_dim_stride  = num_elem * max_num_points * BASIS_NUM_COMP;
 
   // Apply basis element by element
   for (CeedInt elem = blockIdx.x; elem < num_elem; elem += gridDim.x) {
@@ -279,7 +280,7 @@ extern "C" __global__ void GradAtPoints(const CeedInt num_elem, const CeedScalar
 
       // Map to point
       __syncthreads();
-      for (CeedInt p = threadIdx.x; p < BASIS_NUM_PTS; p += blockDim.x) {
+      for (CeedInt p = threadIdx.x; p < max_num_points; p += blockDim.x) {
         for (CeedInt dim_1 = 0; dim_1 < BASIS_DIM; dim_1++) {
           CeedScalar *cur_v = &v[elem * v_stride + dim_1 * v_dim_stride + comp * v_comp_stride];
 
@@ -315,9 +316,10 @@ extern "C" __global__ void GradAtPoints(const CeedInt num_elem, const CeedScalar
   }
 }
 
-extern "C" __global__ void GradTransposeAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ chebyshev_interp_1d,
-                                                 const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ coords,
-                                                 const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v) {
+extern "C" __global__ void GradTransposeAtPoints(const CeedInt num_elem, const CeedInt max_num_points,
+                                                 const CeedScalar *__restrict__ chebyshev_interp_1d, const CeedInt *__restrict__ points_per_elem,
+                                                 const CeedScalar *__restrict__ coords, const CeedScalar *__restrict__ u,
+                                                 CeedScalar *__restrict__ v) {
   const CeedInt i = threadIdx.x;
 
   __shared__ CeedScalar s_mem[BASIS_Q_1D * BASIS_P_1D + 2 * BASIS_BUF_LEN + POINTS_BUFF_LEN * BASIS_Q_1D];
@@ -332,12 +334,12 @@ extern "C" __global__ void GradTransposeAtPoints(const CeedInt num_elem, const C
 
   const CeedInt P             = BASIS_P_1D;
   const CeedInt Q             = BASIS_Q_1D;
-  const CeedInt u_stride      = BASIS_NUM_PTS;
+  const CeedInt u_stride      = max_num_points;
   const CeedInt v_stride      = BASIS_NUM_NODES;
-  const CeedInt u_comp_stride = num_elem * BASIS_NUM_PTS;
+  const CeedInt u_comp_stride = num_elem * max_num_points;
   const CeedInt v_comp_stride = num_elem * BASIS_NUM_NODES;
-  const CeedInt u_size        = BASIS_NUM_PTS;
-  const CeedInt u_dim_stride  = num_elem * BASIS_NUM_PTS * BASIS_NUM_COMP;
+  const CeedInt u_size        = max_num_points;
+  const CeedInt u_dim_stride  = num_elem * max_num_points * BASIS_NUM_COMP;
   const CeedInt v_dim_stride  = 0;
 
   // Apply basis element by element
@@ -354,7 +356,7 @@ extern "C" __global__ void GradTransposeAtPoints(const CeedInt num_elem, const C
 
       // Map from point
       __syncthreads();
-      for (CeedInt p = threadIdx.x; p < BASIS_NUM_PTS; p += blockDim.x) {
+      for (CeedInt p = threadIdx.x; p < max_num_points; p += blockDim.x) {
         if (p >= points_per_elem[elem]) continue;
         for (CeedInt dim_1 = 0; dim_1 < BASIS_DIM; dim_1++) {
           const CeedScalar *cur_u = &u[elem * u_stride + dim_1 * u_dim_stride + comp * u_comp_stride];
diff --git a/include/ceed/jit-source/cuda/cuda-ref-restriction-at-points.h b/include/ceed/jit-source/cuda/cuda-ref-restriction-at-points.h
index 73ecc3bb25..605f2da2cf 100644
--- a/include/ceed/jit-source/cuda/cuda-ref-restriction-at-points.h
+++ b/include/ceed/jit-source/cuda/cuda-ref-restriction-at-points.h
@@ -9,27 +9,44 @@
 /// Internal header for CUDA offset element restriction kernels
 #include <ceed/types.h>
 
+//------------------------------------------------------------------------------
+// L-vector -> E-vector, standard (with offsets)
+//------------------------------------------------------------------------------
+extern "C" __global__ void AtPointsNoTranspose(const CeedInt max_num_points, const CeedInt *__restrict__ indices, const CeedScalar *__restrict__ u,
+                                               CeedScalar *__restrict__ v) {
+  for (CeedInt node = blockIdx.x * blockDim.x + threadIdx.x; node < RSTR_NUM_ELEM * max_num_points; node += blockDim.x * gridDim.x) {
+    const CeedInt ind      = indices[node];
+    const CeedInt loc_node = node % max_num_points;
+    const CeedInt elem     = node / max_num_points;
+
+    for (CeedInt comp = 0; comp < RSTR_NUM_COMP; comp++) {
+      v[loc_node + comp * max_num_points * RSTR_NUM_ELEM + elem * max_num_points] = u[ind + comp * RSTR_COMP_STRIDE];
+    }
+  }
+}
+
 //------------------------------------------------------------------------------
 // E-vector -> L-vector, standard (with offsets)
 //------------------------------------------------------------------------------
 #if !USE_DETERMINISTIC
-extern "C" __global__ void AtPointsTranspose(const CeedInt *__restrict__ indices, const CeedInt *__restrict__ points_per_elem,
-                                             const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v) {
-  for (CeedInt node = blockIdx.x * blockDim.x + threadIdx.x; node < RSTR_NUM_ELEM * RSTR_ELEM_SIZE; node += blockDim.x * gridDim.x) {
+extern "C" __global__ void AtPointsTranspose(const CeedInt max_num_points, const CeedInt *__restrict__ indices,
+                                             const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ u,
+                                             CeedScalar *__restrict__ v) {
+  for (CeedInt node = blockIdx.x * blockDim.x + threadIdx.x; node < RSTR_NUM_ELEM * max_num_points; node += blockDim.x * gridDim.x) {
     const CeedInt ind      = indices[node];
-    const CeedInt loc_node = node % RSTR_ELEM_SIZE;
-    const CeedInt elem     = node / RSTR_ELEM_SIZE;
+    const CeedInt loc_node = node % max_num_points;
+    const CeedInt elem     = node / max_num_points;
 
     if (loc_node >= points_per_elem[elem]) continue;
     for (CeedInt comp = 0; comp < RSTR_NUM_COMP; comp++) {
-      atomicAdd(&v[ind + comp * RSTR_COMP_STRIDE], u[loc_node + comp * RSTR_ELEM_SIZE * RSTR_NUM_ELEM + elem * RSTR_ELEM_SIZE]);
+      atomicAdd(&v[ind + comp * RSTR_COMP_STRIDE], u[loc_node + comp * max_num_points * RSTR_NUM_ELEM + elem * max_num_points]);
     }
   }
 }
 #else
-extern "C" __global__ void AtPointsTranspose(const CeedInt *__restrict__ l_vec_indices, const CeedInt *__restrict__ t_indices,
-                                             const CeedInt *__restrict__ points_per_elem, const CeedInt *__restrict__ t_offsets,
-                                             const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v) {
+extern "C" __global__ void AtPointsTranspose(const CeedInt max_num_points, const CeedInt *__restrict__ l_vec_indices,
+                                             const CeedInt *__restrict__ t_indices, const CeedInt *__restrict__ points_per_elem,
+                                             const CeedInt *__restrict__ t_offsets, const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v) {
   CeedScalar value[RSTR_NUM_COMP];
 
   for (CeedInt i = blockIdx.x * blockDim.x + threadIdx.x; i < RSTR_NUM_NODES; i += blockDim.x * gridDim.x) {
@@ -41,12 +58,12 @@ extern "C" __global__ void AtPointsTranspose(const CeedInt *__restrict__ l_vec_i
 
     for (CeedInt j = range_1; j < range_N; j++) {
       const CeedInt t_ind    = t_indices[j];
-      const CeedInt loc_node = t_ind % RSTR_ELEM_SIZE;
-      const CeedInt elem     = t_ind / RSTR_ELEM_SIZE;
+      const CeedInt loc_node = t_ind % max_num_points;
+      const CeedInt elem     = t_ind / max_num_points;
 
       if (loc_node >= points_per_elem[elem]) continue;
       for (CeedInt comp = 0; comp < RSTR_NUM_COMP; comp++) {
-        value[comp] += u[loc_node + comp * RSTR_ELEM_SIZE * RSTR_NUM_ELEM + elem * RSTR_ELEM_SIZE];
+        value[comp] += u[loc_node + comp * max_num_points * RSTR_NUM_ELEM + elem * max_num_points];
       }
     }
 
diff --git a/include/ceed/jit-source/cuda/cuda-shared-basis-read-write-templates.h b/include/ceed/jit-source/cuda/cuda-shared-basis-read-write-templates.h
index ececd93ae6..f6183d5af7 100644
--- a/include/ceed/jit-source/cuda/cuda-shared-basis-read-write-templates.h
+++ b/include/ceed/jit-source/cuda/cuda-shared-basis-read-write-templates.h
@@ -177,11 +177,11 @@ inline __device__ void SumElementStrided3d(SharedData_Cuda &data, const CeedInt
 //------------------------------------------------------------------------------
 // E-vector -> single point
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_PTS>
-inline __device__ void ReadPoint(SharedData_Cuda &data, const CeedInt elem, const CeedInt p, const CeedInt points_in_elem,
-                                 const CeedInt strides_point, const CeedInt strides_comp, const CeedInt strides_elem,
+template <int NUM_COMP>
+inline __device__ void ReadPoint(SharedData_Cuda &data, const CeedInt elem, const CeedInt p, const CeedInt max_num_points,
+                                 const CeedInt points_in_elem, const CeedInt strides_point, const CeedInt strides_comp, const CeedInt strides_elem,
                                  const CeedScalar *__restrict__ d_u, CeedScalar *r_u) {
-  const CeedInt ind = (p % NUM_PTS) * strides_point + elem * strides_elem;
+  const CeedInt ind = (p % max_num_points) * strides_point + elem * strides_elem;
 
   if (p < points_in_elem) {
     for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
@@ -197,12 +197,12 @@ inline __device__ void ReadPoint(SharedData_Cuda &data, const CeedInt elem, cons
 //------------------------------------------------------------------------------
 // Single point -> E-vector
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_PTS>
-inline __device__ void WritePoint(SharedData_Cuda &data, const CeedInt elem, const CeedInt p, const CeedInt points_in_elem,
-                                  const CeedInt strides_point, const CeedInt strides_comp, const CeedInt strides_elem, const CeedScalar *r_v,
-                                  CeedScalar *d_v) {
+template <int NUM_COMP>
+inline __device__ void WritePoint(SharedData_Cuda &data, const CeedInt elem, const CeedInt p, const CeedInt max_num_points,
+                                  const CeedInt points_in_elem, const CeedInt strides_point, const CeedInt strides_comp, const CeedInt strides_elem,
+                                  const CeedScalar *r_v, CeedScalar *d_v) {
   if (p < points_in_elem) {
-    const CeedInt ind = (p % NUM_PTS) * strides_point + elem * strides_elem;
+    const CeedInt ind = (p % max_num_points) * strides_point + elem * strides_elem;
 
     for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
       d_v[ind + comp * strides_comp] = r_v[comp];
diff --git a/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h b/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h
index bad37bfb81..e805b0404a 100644
--- a/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h
+++ b/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h
@@ -41,7 +41,7 @@ inline __device__ void ChebyshevDerivativeAtPoint(const CeedScalar x, CeedScalar
 //------------------------------------------------------------------------------
 // 1D interpolate to points
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
+template <int NUM_COMP, int P_1D, int Q_1D>
 inline __device__ void InterpAtPoints1d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_C, const CeedScalar *r_X,
                                         CeedScalar *__restrict__ r_V) {
   CeedScalar chebyshev_x[Q_1D];
@@ -62,9 +62,9 @@ inline __device__ void InterpAtPoints1d(SharedData_Cuda &data, const CeedInt p,
 //------------------------------------------------------------------------------
 // 1D interpolate transpose
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
-inline __device__ void InterpTransposeAtPoints1d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_U, const CeedScalar *r_X,
-                                                 CeedScalar *__restrict__ r_C) {
+template <int NUM_COMP, int P_1D, int Q_1D>
+inline __device__ void InterpTransposeAtPoints1d(SharedData_Cuda &data, const CeedInt p, const CeedInt max_num_points,
+                                                 const CeedScalar *__restrict__ r_U, const CeedScalar *r_X, CeedScalar *__restrict__ r_C) {
   CeedScalar chebyshev_x[Q_1D];
 
   ChebyshevPolynomialsAtPoint<Q_1D>(r_X[0], chebyshev_x);
@@ -73,7 +73,7 @@ inline __device__ void InterpTransposeAtPoints1d(SharedData_Cuda &data, const Ce
     if (data.t_id_x < Q_1D) data.slice[data.t_id_x] = 0.0;
     __syncthreads();
     // Contract x direction
-    if (p < NUM_POINTS) {
+    if (p < max_num_points) {
       for (CeedInt i = 0; i < Q_1D; i++) {
         atomicAdd_block(&data.slice[comp * Q_1D + (i + data.t_id_x) % Q_1D], chebyshev_x[(i + data.t_id_x) % Q_1D] * r_U[comp]);
       }
@@ -87,7 +87,7 @@ inline __device__ void InterpTransposeAtPoints1d(SharedData_Cuda &data, const Ce
 //------------------------------------------------------------------------------
 // 1D derivatives at points
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
+template <int NUM_COMP, int P_1D, int Q_1D>
 inline __device__ void GradAtPoints1d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_C, const CeedScalar *r_X,
                                       CeedScalar *__restrict__ r_V) {
   CeedScalar chebyshev_x[Q_1D];
@@ -109,9 +109,9 @@ inline __device__ void GradAtPoints1d(SharedData_Cuda &data, const CeedInt p, co
 //------------------------------------------------------------------------------
 // 1D derivatives transpose
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
-inline __device__ void GradTransposeAtPoints1d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_U, const CeedScalar *r_X,
-                                               CeedScalar *__restrict__ r_C) {
+template <int NUM_COMP, int P_1D, int Q_1D>
+inline __device__ void GradTransposeAtPoints1d(SharedData_Cuda &data, const CeedInt p, const CeedInt max_num_points,
+                                               const CeedScalar *__restrict__ r_U, const CeedScalar *r_X, CeedScalar *__restrict__ r_C) {
   CeedScalar chebyshev_x[Q_1D];
 
   ChebyshevDerivativeAtPoint<Q_1D>(r_X[0], chebyshev_x);
@@ -120,7 +120,7 @@ inline __device__ void GradTransposeAtPoints1d(SharedData_Cuda &data, const Ceed
     if (data.t_id_x < Q_1D) data.slice[data.t_id_x] = 0.0;
     __syncthreads();
     // Contract x direction
-    if (p < NUM_POINTS) {
+    if (p < max_num_points) {
       for (CeedInt i = 0; i < Q_1D; i++) {
         atomicAdd_block(&data.slice[comp * Q_1D + (i + data.t_id_x) % Q_1D], chebyshev_x[(i + data.t_id_x) % Q_1D] * r_U[comp]);
       }
@@ -138,7 +138,7 @@ inline __device__ void GradTransposeAtPoints1d(SharedData_Cuda &data, const Ceed
 //------------------------------------------------------------------------------
 // 2D interpolate to points
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
+template <int NUM_COMP, int P_1D, int Q_1D>
 inline __device__ void InterpAtPoints2d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_C, const CeedScalar *r_X,
                                         CeedScalar *__restrict__ r_V) {
   for (CeedInt i = 0; i < NUM_COMP; i++) r_V[i] = 0.0;
@@ -169,9 +169,9 @@ inline __device__ void InterpAtPoints2d(SharedData_Cuda &data, const CeedInt p,
 //------------------------------------------------------------------------------
 // 2D interpolate transpose
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
-inline __device__ void InterpTransposeAtPoints2d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_U, const CeedScalar *r_X,
-                                                 CeedScalar *__restrict__ r_C) {
+template <int NUM_COMP, int P_1D, int Q_1D>
+inline __device__ void InterpTransposeAtPoints2d(SharedData_Cuda &data, const CeedInt p, const CeedInt max_num_points,
+                                                 const CeedScalar *__restrict__ r_U, const CeedScalar *r_X, CeedScalar *__restrict__ r_C) {
   for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
     CeedScalar buffer[Q_1D];
     CeedScalar chebyshev_x[Q_1D];
@@ -181,7 +181,7 @@ inline __device__ void InterpTransposeAtPoints2d(SharedData_Cuda &data, const Ce
     __syncthreads();
     // Contract y direction
     ChebyshevPolynomialsAtPoint<Q_1D>(r_X[1], chebyshev_x);
-    const CeedScalar r_u = p < NUM_POINTS ? r_U[comp] : 0.0;
+    const CeedScalar r_u = p < max_num_points ? r_U[comp] : 0.0;
 
     for (CeedInt i = 0; i < Q_1D; i++) {
       buffer[i] = chebyshev_x[i] * r_u;
@@ -207,7 +207,7 @@ inline __device__ void InterpTransposeAtPoints2d(SharedData_Cuda &data, const Ce
 //------------------------------------------------------------------------------
 // 2D derivatives at points
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
+template <int NUM_COMP, int P_1D, int Q_1D>
 inline __device__ void GradAtPoints2d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_C, const CeedScalar *r_X,
                                       CeedScalar *__restrict__ r_V) {
   for (CeedInt i = 0; i < NUM_COMP * 2; i++) r_V[i] = 0.0;
@@ -248,9 +248,9 @@ inline __device__ void GradAtPoints2d(SharedData_Cuda &data, const CeedInt p, co
 //------------------------------------------------------------------------------
 // 2D derivatives transpose
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
-inline __device__ void GradTransposeAtPoints2d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_U, const CeedScalar *r_X,
-                                               CeedScalar *__restrict__ r_C) {
+template <int NUM_COMP, int P_1D, int Q_1D>
+inline __device__ void GradTransposeAtPoints2d(SharedData_Cuda &data, const CeedInt p, const CeedInt max_num_points,
+                                               const CeedScalar *__restrict__ r_U, const CeedScalar *r_X, CeedScalar *__restrict__ r_C) {
   for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
     CeedScalar buffer[Q_1D];
     CeedScalar chebyshev_x[Q_1D];
@@ -265,7 +265,7 @@ inline __device__ void GradTransposeAtPoints2d(SharedData_Cuda &data, const Ceed
       } else {
         ChebyshevPolynomialsAtPoint<Q_1D>(r_X[1], chebyshev_x);
       }
-      const CeedScalar r_u = p < NUM_POINTS ? r_U[comp + dim * NUM_COMP] : 0.0;
+      const CeedScalar r_u = p < max_num_points ? r_U[comp + dim * NUM_COMP] : 0.0;
 
       for (CeedInt i = 0; i < Q_1D; i++) {
         buffer[i] = chebyshev_x[i] * r_u;
@@ -300,7 +300,7 @@ inline __device__ void GradTransposeAtPoints2d(SharedData_Cuda &data, const Ceed
 //------------------------------------------------------------------------------
 // 3D interpolate to points
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
+template <int NUM_COMP, int P_1D, int Q_1D>
 inline __device__ void InterpAtPoints3d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_C, const CeedScalar *r_X,
                                         CeedScalar *__restrict__ r_V) {
   for (CeedInt i = 0; i < NUM_COMP; i++) r_V[i] = 0.0;
@@ -337,9 +337,9 @@ inline __device__ void InterpAtPoints3d(SharedData_Cuda &data, const CeedInt p,
 //------------------------------------------------------------------------------
 // 3D interpolate transpose
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
-inline __device__ void InterpTransposeAtPoints3d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_U, const CeedScalar *r_X,
-                                                 CeedScalar *__restrict__ r_C) {
+template <int NUM_COMP, int P_1D, int Q_1D>
+inline __device__ void InterpTransposeAtPoints3d(SharedData_Cuda &data, const CeedInt p, const CeedInt max_num_points,
+                                                 const CeedScalar *__restrict__ r_U, const CeedScalar *r_X, CeedScalar *__restrict__ r_C) {
   for (CeedInt k = 0; k < Q_1D; k++) {
     CeedScalar buffer[Q_1D];
     CeedScalar chebyshev_x[Q_1D];
@@ -354,7 +354,7 @@ inline __device__ void InterpTransposeAtPoints3d(SharedData_Cuda &data, const Ce
       __syncthreads();
       // Contract y and z direction
       ChebyshevPolynomialsAtPoint<Q_1D>(r_X[1], chebyshev_x);
-      const CeedScalar r_u = p < NUM_POINTS ? r_U[comp] : 0.0;
+      const CeedScalar r_u = p < max_num_points ? r_U[comp] : 0.0;
 
       for (CeedInt i = 0; i < Q_1D; i++) {
         buffer[i] = chebyshev_x[i] * r_u * z;
@@ -381,7 +381,7 @@ inline __device__ void InterpTransposeAtPoints3d(SharedData_Cuda &data, const Ce
 //------------------------------------------------------------------------------
 // 3D derivatives at points
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
+template <int NUM_COMP, int P_1D, int Q_1D>
 inline __device__ void GradAtPoints3d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_C, const CeedScalar *r_X,
                                       CeedScalar *__restrict__ r_V) {
   for (CeedInt i = 0; i < NUM_COMP * 3; i++) r_V[i] = 0.0;
@@ -434,9 +434,9 @@ inline __device__ void GradAtPoints3d(SharedData_Cuda &data, const CeedInt p, co
 //------------------------------------------------------------------------------
 // 3D derivatives transpose
 //------------------------------------------------------------------------------
-template <int NUM_COMP, int NUM_POINTS, int P_1D, int Q_1D>
-inline __device__ void GradTransposeAtPoints3d(SharedData_Cuda &data, const CeedInt p, const CeedScalar *__restrict__ r_U, const CeedScalar *r_X,
-                                               CeedScalar *__restrict__ r_C) {
+template <int NUM_COMP, int P_1D, int Q_1D>
+inline __device__ void GradTransposeAtPoints3d(SharedData_Cuda &data, const CeedInt p, const CeedInt max_num_points,
+                                               const CeedScalar *__restrict__ r_U, const CeedScalar *r_X, CeedScalar *__restrict__ r_C) {
   for (CeedInt k = 0; k < Q_1D; k++) {
     CeedScalar buffer[Q_1D];
     CeedScalar chebyshev_x[Q_1D];
@@ -461,7 +461,7 @@ inline __device__ void GradTransposeAtPoints3d(SharedData_Cuda &data, const Ceed
           ChebyshevPolynomialsAtPoint<Q_1D>(r_X[1], chebyshev_x);
         }
         const CeedScalar zz  = dim == 2 ? dz : z;
-        const CeedScalar r_u = (p < NUM_POINTS) ? r_U[comp + dim * NUM_COMP] : 0.0;
+        const CeedScalar r_u = (p < max_num_points) ? r_U[comp + dim * NUM_COMP] : 0.0;
 
         for (CeedInt i = 0; i < Q_1D; i++) {
           buffer[i] = chebyshev_x[i] * r_u * zz;
diff --git a/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points.h b/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points.h
index fc812792e4..7c2ecfe22b 100644
--- a/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points.h
+++ b/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points.h
@@ -20,8 +20,9 @@
 //------------------------------------------------------------------------------
 // Interp
 //------------------------------------------------------------------------------
-extern "C" __global__ void InterpAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ c_B, const CeedInt *__restrict__ points_per_elem,
-                                          const CeedScalar *__restrict__ d_X, const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
+extern "C" __global__ void InterpAtPoints(const CeedInt num_elem, const CeedInt max_num_points, const CeedScalar *__restrict__ c_B,
+                                          const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ d_X,
+                                          const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
   extern __shared__ CeedScalar slice[];
 
   SharedData_Cuda data;
@@ -57,25 +58,25 @@ extern "C" __global__ void InterpAtPoints(const CeedInt num_elem, const CeedScal
     }
 
     // Map to points
-    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * BASIS_NUM_PTS / (blockDim.x * blockDim.y));
+    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * max_num_points / (blockDim.x * blockDim.y));
 
     for (CeedInt i = threadIdx.x + threadIdx.y * blockDim.x; i < point_loop_bound; i += blockDim.x * blockDim.y) {
-      const CeedInt p = i % BASIS_NUM_PTS;
+      const CeedInt p = i % max_num_points;
 
-      ReadPoint<BASIS_DIM, BASIS_NUM_PTS>(data, elem, p, BASIS_NUM_PTS, 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_X, r_X);
+      ReadPoint<BASIS_DIM>(data, elem, p, max_num_points, max_num_points, 1, num_elem * max_num_points, max_num_points, d_X, r_X);
       if (BASIS_DIM == 1) {
-        InterpAtPoints1d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
+        InterpAtPoints1d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
       } else if (BASIS_DIM == 2) {
-        InterpAtPoints2d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
+        InterpAtPoints2d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
       } else if (BASIS_DIM == 3) {
-        InterpAtPoints3d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
+        InterpAtPoints3d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
       }
-      WritePoint<BASIS_NUM_COMP, BASIS_NUM_PTS>(data, elem, p, BASIS_NUM_PTS, 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, r_V, d_V);
+      WritePoint<BASIS_NUM_COMP>(data, elem, p, max_num_points, max_num_points, 1, num_elem * max_num_points, max_num_points, r_V, d_V);
     }
   }
 }
 
-extern "C" __global__ void InterpTransposeAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ c_B,
+extern "C" __global__ void InterpTransposeAtPoints(const CeedInt num_elem, const CeedInt max_num_points, const CeedScalar *__restrict__ c_B,
                                                    const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ d_X,
                                                    const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
   extern __shared__ CeedScalar slice[];
@@ -114,19 +115,19 @@ extern "C" __global__ void InterpTransposeAtPoints(const CeedInt num_elem, const
     }
 
     // Map from points
-    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * BASIS_NUM_PTS / (blockDim.x * blockDim.y));
+    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * max_num_points / (blockDim.x * blockDim.y));
 
     for (CeedInt i = threadIdx.x + threadIdx.y * blockDim.x; i < point_loop_bound; i += blockDim.x * blockDim.y) {
-      const CeedInt p = i % BASIS_NUM_PTS;
+      const CeedInt p = i % max_num_points;
 
-      ReadPoint<BASIS_DIM, BASIS_NUM_PTS>(data, elem, p, BASIS_NUM_PTS, 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_X, r_X);
-      ReadPoint<BASIS_NUM_COMP, BASIS_NUM_PTS>(data, elem, i, points_per_elem[elem], 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_U, r_U);
+      ReadPoint<BASIS_DIM>(data, elem, p, max_num_points, max_num_points, 1, num_elem * max_num_points, max_num_points, d_X, r_X);
+      ReadPoint<BASIS_NUM_COMP>(data, elem, i, max_num_points, points_per_elem[elem], 1, num_elem * max_num_points, max_num_points, d_U, r_U);
       if (BASIS_DIM == 1) {
-        InterpTransposeAtPoints1d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        InterpTransposeAtPoints1d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       } else if (BASIS_DIM == 2) {
-        InterpTransposeAtPoints2d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        InterpTransposeAtPoints2d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       } else if (BASIS_DIM == 3) {
-        InterpTransposeAtPoints3d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        InterpTransposeAtPoints3d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       }
     }
 
@@ -145,7 +146,7 @@ extern "C" __global__ void InterpTransposeAtPoints(const CeedInt num_elem, const
   }
 }
 
-extern "C" __global__ void InterpTransposeAddAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ c_B,
+extern "C" __global__ void InterpTransposeAddAtPoints(const CeedInt num_elem, const CeedInt max_num_points, const CeedScalar *__restrict__ c_B,
                                                       const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ d_X,
                                                       const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
   extern __shared__ CeedScalar slice[];
@@ -173,19 +174,19 @@ extern "C" __global__ void InterpTransposeAddAtPoints(const CeedInt num_elem, co
     for (CeedInt i = 0; i < BASIS_NUM_COMP * (BASIS_DIM > 2 ? BASIS_Q_1D : 1); i++) r_C[i] = 0.0;
 
     // Map from points
-    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * BASIS_NUM_PTS / (blockDim.x * blockDim.y));
+    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * max_num_points / (blockDim.x * blockDim.y));
 
     for (CeedInt i = threadIdx.x + threadIdx.y * blockDim.x; i < point_loop_bound; i += blockDim.x * blockDim.y) {
-      const CeedInt p = i % BASIS_NUM_PTS;
+      const CeedInt p = i % max_num_points;
 
-      ReadPoint<BASIS_DIM, BASIS_NUM_PTS>(data, elem, p, BASIS_NUM_PTS, 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_X, r_X);
-      ReadPoint<BASIS_NUM_COMP, BASIS_NUM_PTS>(data, elem, i, points_per_elem[elem], 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_U, r_U);
+      ReadPoint<BASIS_DIM>(data, elem, p, max_num_points, max_num_points, 1, num_elem * max_num_points, max_num_points, d_X, r_X);
+      ReadPoint<BASIS_NUM_COMP>(data, elem, i, max_num_points, points_per_elem[elem], 1, num_elem * max_num_points, max_num_points, d_U, r_U);
       if (BASIS_DIM == 1) {
-        InterpTransposeAtPoints1d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        InterpTransposeAtPoints1d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       } else if (BASIS_DIM == 2) {
-        InterpTransposeAtPoints2d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        InterpTransposeAtPoints2d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       } else if (BASIS_DIM == 3) {
-        InterpTransposeAtPoints3d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        InterpTransposeAtPoints3d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       }
     }
 
@@ -207,8 +208,9 @@ extern "C" __global__ void InterpTransposeAddAtPoints(const CeedInt num_elem, co
 //------------------------------------------------------------------------------
 // Grad
 //------------------------------------------------------------------------------
-extern "C" __global__ void GradAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ c_B, const CeedInt *__restrict__ points_per_elem,
-                                        const CeedScalar *__restrict__ d_X, const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
+extern "C" __global__ void GradAtPoints(const CeedInt num_elem, const CeedInt max_num_points, const CeedScalar *__restrict__ c_B,
+                                        const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ d_X,
+                                        const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
   extern __shared__ CeedScalar slice[];
 
   SharedData_Cuda data;
@@ -244,25 +246,25 @@ extern "C" __global__ void GradAtPoints(const CeedInt num_elem, const CeedScalar
     }
 
     // Map to points
-    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * BASIS_NUM_PTS / (blockDim.x * blockDim.y));
+    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * max_num_points / (blockDim.x * blockDim.y));
 
     for (CeedInt i = threadIdx.x + threadIdx.y * blockDim.x; i < point_loop_bound; i += blockDim.x * blockDim.y) {
-      const CeedInt p = i % BASIS_NUM_PTS;
+      const CeedInt p = i % max_num_points;
 
-      ReadPoint<BASIS_DIM, BASIS_NUM_PTS>(data, elem, p, BASIS_NUM_PTS, 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_X, r_X);
+      ReadPoint<BASIS_DIM>(data, elem, p, max_num_points, max_num_points, 1, num_elem * max_num_points, max_num_points, d_X, r_X);
       if (BASIS_DIM == 1) {
-        GradAtPoints1d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
+        GradAtPoints1d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
       } else if (BASIS_DIM == 2) {
-        GradAtPoints2d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
+        GradAtPoints2d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
       } else if (BASIS_DIM == 3) {
-        GradAtPoints3d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
+        GradAtPoints3d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, r_C, r_X, r_V);
       }
-      WritePoint<BASIS_NUM_COMP * BASIS_DIM, BASIS_NUM_PTS>(data, elem, p, BASIS_NUM_PTS, 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, r_V, d_V);
+      WritePoint<BASIS_NUM_COMP * BASIS_DIM>(data, elem, p, max_num_points, max_num_points, 1, num_elem * max_num_points, max_num_points, r_V, d_V);
     }
   }
 }
 
-extern "C" __global__ void GradTransposeAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ c_B,
+extern "C" __global__ void GradTransposeAtPoints(const CeedInt num_elem, const CeedInt max_num_points, const CeedScalar *__restrict__ c_B,
                                                  const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ d_X,
                                                  const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
   extern __shared__ CeedScalar slice[];
@@ -301,20 +303,20 @@ extern "C" __global__ void GradTransposeAtPoints(const CeedInt num_elem, const C
     }
 
     // Map from points
-    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * BASIS_NUM_PTS / (blockDim.x * blockDim.y));
+    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * max_num_points / (blockDim.x * blockDim.y));
 
     for (CeedInt i = threadIdx.x + threadIdx.y * blockDim.x; i < point_loop_bound; i += blockDim.x * blockDim.y) {
-      const CeedInt p = i % BASIS_NUM_PTS;
+      const CeedInt p = i % max_num_points;
 
-      ReadPoint<BASIS_DIM, BASIS_NUM_PTS>(data, elem, p, BASIS_NUM_PTS, 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_X, r_X);
-      ReadPoint<BASIS_NUM_COMP * BASIS_DIM, BASIS_NUM_PTS>(data, elem, i, points_per_elem[elem], 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_U,
-                                                           r_U);
+      ReadPoint<BASIS_DIM>(data, elem, p, max_num_points, max_num_points, 1, num_elem * max_num_points, max_num_points, d_X, r_X);
+      ReadPoint<BASIS_NUM_COMP * BASIS_DIM>(data, elem, i, max_num_points, points_per_elem[elem], 1, num_elem * max_num_points, max_num_points, d_U,
+                                            r_U);
       if (BASIS_DIM == 1) {
-        GradTransposeAtPoints1d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        GradTransposeAtPoints1d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       } else if (BASIS_DIM == 2) {
-        GradTransposeAtPoints2d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        GradTransposeAtPoints2d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       } else if (BASIS_DIM == 3) {
-        GradTransposeAtPoints3d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        GradTransposeAtPoints3d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       }
     }
 
@@ -333,7 +335,7 @@ extern "C" __global__ void GradTransposeAtPoints(const CeedInt num_elem, const C
   }
 }
 
-extern "C" __global__ void GradTransposeAddAtPoints(const CeedInt num_elem, const CeedScalar *__restrict__ c_B,
+extern "C" __global__ void GradTransposeAddAtPoints(const CeedInt num_elem, const CeedInt max_num_points, const CeedScalar *__restrict__ c_B,
                                                     const CeedInt *__restrict__ points_per_elem, const CeedScalar *__restrict__ d_X,
                                                     const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
   extern __shared__ CeedScalar slice[];
@@ -361,20 +363,20 @@ extern "C" __global__ void GradTransposeAddAtPoints(const CeedInt num_elem, cons
     for (CeedInt i = 0; i < BASIS_NUM_COMP * (BASIS_DIM > 2 ? BASIS_Q_1D : 1); i++) r_C[i] = 0.0;
 
     // Map from points
-    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * BASIS_NUM_PTS / (blockDim.x * blockDim.y));
+    const CeedInt point_loop_bound = (blockDim.x * blockDim.y) * ceil(1.0 * max_num_points / (blockDim.x * blockDim.y));
 
     for (CeedInt i = threadIdx.x + threadIdx.y * blockDim.x; i < point_loop_bound; i += blockDim.x * blockDim.y) {
-      const CeedInt p = i % BASIS_NUM_PTS;
+      const CeedInt p = i % max_num_points;
 
-      ReadPoint<BASIS_DIM, BASIS_NUM_PTS>(data, elem, p, BASIS_NUM_PTS, 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_X, r_X);
-      ReadPoint<BASIS_NUM_COMP * BASIS_DIM, BASIS_NUM_PTS>(data, elem, i, points_per_elem[elem], 1, num_elem * BASIS_NUM_PTS, BASIS_NUM_PTS, d_U,
-                                                           r_U);
+      ReadPoint<BASIS_DIM>(data, elem, p, max_num_points, max_num_points, 1, num_elem * max_num_points, max_num_points, d_X, r_X);
+      ReadPoint<BASIS_NUM_COMP * BASIS_DIM>(data, elem, i, max_num_points, points_per_elem[elem], 1, num_elem * max_num_points, max_num_points, d_U,
+                                            r_U);
       if (BASIS_DIM == 1) {
-        GradTransposeAtPoints1d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        GradTransposeAtPoints1d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       } else if (BASIS_DIM == 2) {
-        GradTransposeAtPoints2d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        GradTransposeAtPoints2d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       } else if (BASIS_DIM == 3) {
-        GradTransposeAtPoints3d<BASIS_NUM_COMP, BASIS_NUM_PTS, BASIS_P_1D, BASIS_Q_1D>(data, i, r_U, r_X, r_C);
+        GradTransposeAtPoints3d<BASIS_NUM_COMP, BASIS_P_1D, BASIS_Q_1D>(data, i, max_num_points, r_U, r_X, r_C);
       }
     }
 
diff --git a/include/ceed/jit-source/cuda/cuda-types.h b/include/ceed/jit-source/cuda/cuda-types.h
index 58b2961246..ab8f5234f2 100644
--- a/include/ceed/jit-source/cuda/cuda-types.h
+++ b/include/ceed/jit-source/cuda/cuda-types.h
@@ -25,6 +25,8 @@ typedef struct {
 
 typedef struct {
   CeedInt           num_elem;
+  CeedInt           max_num_points;
+  CeedInt           coords_comp_stride;
   const CeedInt    *num_per_elem;
   const CeedInt    *indices;
   const CeedScalar *coords;