diff --git a/framework/include/base/Adaptivity.h b/framework/include/base/Adaptivity.h
index f87ad8a7dcbc..2e8fc651fe72 100644
--- a/framework/include/base/Adaptivity.h
+++ b/framework/include/base/Adaptivity.h
@@ -50,7 +50,7 @@ class Adaptivity : public ConsoleStreamInterface,
                    public libMesh::ParallelObject
 {
 public:
-  Adaptivity(FEProblemBase & subproblem);
+  Adaptivity(FEProblemBase & fe_problem);
   virtual ~Adaptivity();
 
   /**
@@ -128,6 +128,11 @@ class Adaptivity : public ConsoleStreamInterface,
    */
   void setRecomputeMarkersFlag(const bool flag) { _recompute_markers_during_cycles = flag; }
 
+  /**
+   * Switch from h-refinement to p-refinement
+   */
+  void switchHToPRefinement();
+
   /**
    * Adapts the mesh based on the error estimator used
    *
@@ -245,7 +250,7 @@ class Adaptivity : public ConsoleStreamInterface,
   bool isAdaptivityDue();
 
 protected:
-  FEProblemBase & _subproblem;
+  FEProblemBase & _fe_problem;
   MooseMesh & _mesh;
 
   /// on/off flag reporting if the adaptivity is being used
@@ -302,6 +307,8 @@ class Adaptivity : public ConsoleStreamInterface,
 
   /// Stores pointers to ErrorVectors associated with indicator field names
   std::map<std::string, std::unique_ptr<ErrorVector>> _indicator_field_to_error_vector;
+
+  bool _p_refinement_flag = false;
 };
 
 template <typename T>
diff --git a/framework/include/base/Assembly.h b/framework/include/base/Assembly.h
index c96b61d8c847..0e610a2283f4 100644
--- a/framework/include/base/Assembly.h
+++ b/framework/include/base/Assembly.h
@@ -568,6 +568,11 @@ class Assembly
    */
   bool needDual() const { return _need_dual; }
 
+  /**
+   * Set the cached quadrature rules to nullptr
+   */
+  void clearCachedQRules();
+
 private:
   /**
    * Set the qrule to be used for lower dimensional integration.
@@ -595,6 +600,11 @@ class Assembly
    */
   void reinit(const Elem * elem);
 
+  /**
+   * Set the volumetric quadrature rule based on the provided element
+   */
+  void setVolumeQRule(const Elem * elem);
+
   /**
    * Reinitialize FE data for the given element on the given side, optionally
    * with a given set of reference points
@@ -664,6 +674,11 @@ class Assembly
    */
   void reinit(const Elem * elem, const std::vector<Point> & reference_points);
 
+  /**
+   * Set the face quadrature rule based on the provided element and side
+   */
+  void setFaceQRule(const Elem * const elem, const unsigned int side);
+
   /**
    * Reinitialize the assembly data on an side of an element
    */
@@ -1817,6 +1832,11 @@ class Assembly
    */
   const Elem * const & msmElem() const { return _msm_elem; }
 
+  /**
+   * Indicate that we have p-refinement
+   */
+  void havePRefinement();
+
 private:
   /**
    * Just an internal helper function to reinit the volume FE objects.
@@ -2150,6 +2170,11 @@ class Assembly
     return _jacobian_block_nonlocal_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
   }
 
+  /**
+   * request phi, dphi, xyz, JxW, etc. data through the FE helper functions
+   */
+  void helpersRequestData();
+
   SystemBase & _sys;
   SubProblem & _subproblem;
 
@@ -2194,6 +2219,28 @@ class Assembly
 
   unsigned int _mesh_dimension;
 
+  /// The finite element type of the FE helper classes. The helper class gives us data like JxW, the
+  /// physical quadrature point locations, etc.
+  const FEType _helper_type;
+
+  /// Whether user code requested a \p FEType the same as our \p _helper_type
+  mutable bool _user_added_fe_of_helper_type;
+  mutable bool _user_added_fe_face_of_helper_type;
+  mutable bool _user_added_fe_face_neighbor_of_helper_type;
+  mutable bool _user_added_fe_neighbor_of_helper_type;
+  mutable bool _user_added_fe_lower_of_helper_type;
+
+  /// Containers for holding unique FE helper types if we are doing p-refinement. If we are not
+  /// doing p-refinement then the helper data is owned by the \p _fe data members
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_helper;
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_face_helper;
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_face_neighbor_helper;
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_neighbor_helper;
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_lower_helper;
+
+  /// Whether we are currently building the FE classes for the helpers
+  bool _building_helpers;
+
   /// The XFEM controller
   std::shared_ptr<XFEMInterface> _xfem;
 
@@ -2222,7 +2269,7 @@ class Assembly
   /// Each dimension's actual vector fe objects indexed on type
   mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe;
   /// Each dimension's helper objects
-  std::map<unsigned int, FEBase **> _holder_fe_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_helper;
   /// The current helper object for transforming coordinates
   FEBase * _current_fe_helper;
   /// The current current quadrature rule being used (could be either volumetric or arbitrary - for dirac kernels)
@@ -2334,7 +2381,7 @@ class Assembly
   /// types of vector finite elements
   mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_face;
   /// Each dimension's helper objects
-  std::map<unsigned int, FEBase **> _holder_fe_face_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_face_helper;
   /// helper object for transforming coordinates
   FEBase * _current_fe_face_helper;
   /// quadrature rule used on faces
@@ -2370,15 +2417,15 @@ class Assembly
   mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_face_neighbor;
 
   /// Each dimension's helper objects
-  std::map<unsigned int, FEBase **> _holder_fe_neighbor_helper;
-  std::map<unsigned int, FEBase **> _holder_fe_face_neighbor_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_neighbor_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_face_neighbor_helper;
 
   /// FE objects for lower dimensional elements
   mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_lower;
   /// Vector FE objects for lower dimensional elements
   mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_lower;
   /// helper object for transforming coordinates for lower dimensional element quadrature points
-  std::map<unsigned int, FEBase **> _holder_fe_lower_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_lower_helper;
 
   /// quadrature rule used on neighbors
   QBase * _current_qrule_neighbor;
diff --git a/framework/include/loops/MaxQpsThread.h b/framework/include/loops/MaxQpsThread.h
index 11a0c1f146b6..9b3c003c5633 100644
--- a/framework/include/loops/MaxQpsThread.h
+++ b/framework/include/loops/MaxQpsThread.h
@@ -38,8 +38,6 @@ class MaxQpsThread
 
   unsigned int max() const { return _max; }
 
-  unsigned int max_shape_funcs() const { return _max_shape_funcs; }
-
 protected:
   FEProblemBase & _fe_problem;
 
@@ -47,7 +45,4 @@ class MaxQpsThread
 
   /// Maximum number of qps encountered
   unsigned int _max;
-
-  /// Maximum number of shape functions encountered
-  unsigned int _max_shape_funcs;
 };
diff --git a/framework/include/problems/DisplacedProblem.h b/framework/include/problems/DisplacedProblem.h
index 874640e74512..56dcc638a2fa 100644
--- a/framework/include/problems/DisplacedProblem.h
+++ b/framework/include/problems/DisplacedProblem.h
@@ -354,6 +354,11 @@ class DisplacedProblem : public SubProblem
 
   virtual const std::vector<VectorTag> & currentResidualVectorTags() const override;
 
+  /**
+   * Indicate that we have p-refinement
+   */
+  void havePRefinement();
+
 protected:
   FEProblemBase & _mproblem;
   MooseMesh & _mesh;
diff --git a/framework/include/problems/FEProblemBase.h b/framework/include/problems/FEProblemBase.h
index 25fe5514e39d..4e71b3065f47 100644
--- a/framework/include/problems/FEProblemBase.h
+++ b/framework/include/problems/FEProblemBase.h
@@ -339,11 +339,6 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   unsigned int getMaxQps() const;
 
-  /**
-   * @return The maximum number of quadrature points in use on any element in this problem.
-   */
-  unsigned int getMaxShapeFunctions() const;
-
   /**
    * @return The maximum order for all scalar variables in this problem's systems.
    */
@@ -2102,6 +2097,11 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   void clearCurrentResidualVectorTags();
 
+  /**
+   * Indicate that we have p-refinement
+   */
+  void havePRefinement();
+
 protected:
   /// Create extra tagged vectors and matrices
   void createTagVectors();
@@ -2384,9 +2384,6 @@ class FEProblemBase : public SubProblem, public Restartable
   /// Maximum number of quadrature points used in the problem
   unsigned int _max_qps;
 
-  /// Maximum number of shape functions on any element in the problem
-  unsigned int _max_shape_funcs;
-
   /// Maximum scalar variable order
   Order _max_scalar_order;
 
diff --git a/framework/src/actions/AdaptivityAction.C b/framework/src/actions/AdaptivityAction.C
index 5cb1c2dfae65..77a92a5af009 100644
--- a/framework/src/actions/AdaptivityAction.C
+++ b/framework/src/actions/AdaptivityAction.C
@@ -79,6 +79,7 @@ AdaptivityAction::validParams()
       "show_initial_progress", true, "Show the progress of the initial adaptivity");
   params.addParam<bool>(
       "recompute_markers_during_cycles", false, "Recompute markers during adaptivity cycles");
+  params.addParam<bool>("switch_h_to_p_refinement", false, "True to perform p-refinement");
   return params;
 }
 
@@ -202,6 +203,8 @@ AdaptivityAction::act()
 
     adapt.setTimeActive(getParam<Real>("start_time"), getParam<Real>("stop_time"));
     adapt.setInterval(getParam<unsigned int>("interval"));
+    if (getParam<bool>("switch_h_to_p_refinement"))
+      adapt.switchHToPRefinement();
   }
 }
 
diff --git a/framework/src/actions/SetAdaptivityOptionsAction.C b/framework/src/actions/SetAdaptivityOptionsAction.C
index 8bb4a840d70f..261d641a8fe3 100644
--- a/framework/src/actions/SetAdaptivityOptionsAction.C
+++ b/framework/src/actions/SetAdaptivityOptionsAction.C
@@ -49,6 +49,7 @@ SetAdaptivityOptionsAction::validParams()
       "The number of adaptive steps to use when on each timestep during a Transient simulation.");
   params.addParam<bool>(
       "recompute_markers_during_cycles", false, "Recompute markers during adaptivity cycles");
+  params.addParam<bool>("switch_h_to_p_refinement", false, "True to perform p-refinement");
   return params;
 }
 
@@ -137,5 +138,7 @@ SetAdaptivityOptionsAction::act()
     adapt.setInterval(getParam<unsigned int>("interval"));
 
     adapt.setRecomputeMarkersFlag(getParam<bool>("recompute_markers_during_cycles"));
+    if (getParam<bool>("switch_h_to_p_refinement"))
+      adapt.switchHToPRefinement();
   }
 }
diff --git a/framework/src/base/Adaptivity.C b/framework/src/base/Adaptivity.C
index d1c30fadf1ca..064a401879cd 100644
--- a/framework/src/base/Adaptivity.C
+++ b/framework/src/base/Adaptivity.C
@@ -28,19 +28,19 @@
 
 #ifdef LIBMESH_ENABLE_AMR
 
-Adaptivity::Adaptivity(FEProblemBase & subproblem)
-  : ConsoleStreamInterface(subproblem.getMooseApp()),
-    PerfGraphInterface(subproblem.getMooseApp().perfGraph(), "Adaptivity"),
-    ParallelObject(subproblem.getMooseApp()),
-    _subproblem(subproblem),
-    _mesh(_subproblem.mesh()),
+Adaptivity::Adaptivity(FEProblemBase & fe_problem)
+  : ConsoleStreamInterface(fe_problem.getMooseApp()),
+    PerfGraphInterface(fe_problem.getMooseApp().perfGraph(), "Adaptivity"),
+    ParallelObject(fe_problem.getMooseApp()),
+    _fe_problem(fe_problem),
+    _mesh(_fe_problem.mesh()),
     _mesh_refinement_on(false),
     _initialized(false),
     _initial_steps(0),
     _steps(0),
     _print_mesh_changed(false),
-    _t(_subproblem.time()),
-    _step(_subproblem.timeStep()),
+    _t(_fe_problem.time()),
+    _step(_fe_problem.timeStep()),
     _interval(1),
     _start_time(-std::numeric_limits<Real>::max()),
     _stop_time(std::numeric_limits<Real>::max()),
@@ -59,12 +59,12 @@ Adaptivity::init(unsigned int steps, unsigned int initial_steps)
   // Get the pointer to the DisplacedProblem, this cannot be done at construction because
   // DisplacedProblem
   // does not exist at that point.
-  _displaced_problem = _subproblem.getDisplacedProblem();
+  _displaced_problem = _fe_problem.getDisplacedProblem();
 
   _mesh_refinement = std::make_unique<MeshRefinement>(_mesh);
   _error = std::make_unique<ErrorVector>();
 
-  EquationSystems & es = _subproblem.es();
+  EquationSystems & es = _fe_problem.es();
   es.parameters.set<bool>("adaptivity") = true;
 
   _initial_steps = initial_steps;
@@ -72,7 +72,7 @@ Adaptivity::init(unsigned int steps, unsigned int initial_steps)
   _mesh_refinement_on = true;
 
   _mesh_refinement->set_periodic_boundaries_ptr(
-      _subproblem.getNonlinearSystemBase().dofMap().get_periodic_boundaries());
+      _fe_problem.getNonlinearSystemBase().dofMap().get_periodic_boundaries());
 
   // displaced problem
   if (_displaced_problem != nullptr)
@@ -88,7 +88,7 @@ Adaptivity::init(unsigned int steps, unsigned int initial_steps)
     // i.e. neighbors across periodic boundaries, for the purposes of
     // refinement.
     _displaced_mesh_refinement->set_periodic_boundaries_ptr(
-        _subproblem.getNonlinearSystemBase().dofMap().get_periodic_boundaries());
+        _fe_problem.getNonlinearSystemBase().dofMap().get_periodic_boundaries());
 
     // TODO: This is currently an empty function on the DisplacedProblem... could it be removed?
     _displaced_problem->initAdaptivity();
@@ -141,7 +141,7 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
 
       std::vector<Number> serialized_solution;
 
-      auto distributed_mesh = dynamic_cast<DistributedMesh *>(&_subproblem.mesh().getMesh());
+      auto distributed_mesh = dynamic_cast<DistributedMesh *>(&_fe_problem.mesh().getMesh());
 
       // Element range
       std::unique_ptr<ConstElemRange> all_elems;
@@ -152,7 +152,7 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
       if (distributed_mesh && !distributed_mesh->is_serial_on_zero())
       {
         // We update here to make sure local solution is up-to-date
-        _subproblem.getAuxiliarySystem().update();
+        _fe_problem.getAuxiliarySystem().update();
         distributed_adaptivity = true;
 
         // We can not assume that geometric and algebraic ghosting functors cover
@@ -166,13 +166,13 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
         // After we set markers for all local elements, we will do a global
         // communication to sync markers for ghosted elements from their owners.
         all_elems = std::make_unique<ConstElemRange>(
-            _subproblem.mesh().getMesh().active_local_elements_begin(),
-            _subproblem.mesh().getMesh().active_local_elements_end());
+            _fe_problem.mesh().getMesh().active_local_elements_begin(),
+            _fe_problem.mesh().getMesh().active_local_elements_end());
       }
       else // This is not scalable but it might be useful for small-size problems
       {
-        _subproblem.getAuxiliarySystem().solution().close();
-        _subproblem.getAuxiliarySystem().solution().localize(serialized_solution);
+        _fe_problem.getAuxiliarySystem().solution().close();
+        _fe_problem.getAuxiliarySystem().solution().localize(serialized_solution);
         distributed_adaptivity = false;
 
         // For a replicated mesh or a serialized distributed mesh, the solution
@@ -182,20 +182,20 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
         // We might not care about much since a replicated mesh
         // or a serialized distributed mesh is not scalable anyway.
         all_elems =
-            std::make_unique<ConstElemRange>(_subproblem.mesh().getMesh().active_elements_begin(),
-                                             _subproblem.mesh().getMesh().active_elements_end());
+            std::make_unique<ConstElemRange>(_fe_problem.mesh().getMesh().active_elements_begin(),
+                                             _fe_problem.mesh().getMesh().active_elements_end());
       }
 
       FlagElementsThread fet(
-          _subproblem, serialized_solution, _max_h_level, marker_name, !distributed_adaptivity);
+          _fe_problem, serialized_solution, _max_h_level, marker_name, !distributed_adaptivity);
       Threads::parallel_reduce(*all_elems, fet);
-      _subproblem.getAuxiliarySystem().solution().close();
+      _fe_problem.getAuxiliarySystem().solution().close();
     }
   }
   else
   {
     // Compute the error for each active element
-    _error_estimator->estimate_error(_subproblem.getNonlinearSystemBase().system(), *_error);
+    _error_estimator->estimate_error(_fe_problem.getNonlinearSystemBase().system(), *_error);
 
     // Flag elements to be refined and coarsened
     _mesh_refinement->flag_elements_by_error_fraction(*_error);
@@ -220,6 +220,9 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
   if (distributed_adaptivity)
     _mesh_refinement->make_flags_parallel_consistent();
 
+  if (_p_refinement_flag)
+    _mesh_refinement->switch_h_to_p_refinement();
+
   // Perform refinement and coarsening
   mesh_changed = _mesh_refinement->refine_and_coarsen_elements();
 
@@ -229,6 +232,10 @@ Adaptivity::adaptMesh(std::string marker_name /*=std::string()*/)
     // we sync them here.
     if (distributed_adaptivity)
       _displaced_mesh_refinement->make_flags_parallel_consistent();
+
+    if (_p_refinement_flag)
+      _displaced_mesh_refinement->switch_h_to_p_refinement();
+
 #ifndef NDEBUG
     bool displaced_mesh_changed =
 #endif
@@ -303,7 +310,7 @@ Adaptivity::uniformRefineWithProjection()
 
     if (_displaced_problem)
       displaced_mesh_refinement.uniformly_refine(1);
-    _subproblem.meshChanged();
+    _fe_problem.meshChanged();
   }
 }
 
@@ -364,12 +371,12 @@ Adaptivity::updateErrorVectors()
   }
 
   // Fill the vectors with the local contributions
-  UpdateErrorVectorsThread uevt(_subproblem, _indicator_field_to_error_vector);
+  UpdateErrorVectorsThread uevt(_fe_problem, _indicator_field_to_error_vector);
   Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), uevt);
 
   // Now sum across all processors
   for (const auto & it : _indicator_field_to_error_vector)
-    _subproblem.comm().sum((std::vector<float> &)*(it.second));
+    _fe_problem.comm().sum((std::vector<float> &)*(it.second));
 }
 
 bool
@@ -378,4 +385,11 @@ Adaptivity::isAdaptivityDue()
   return _mesh_refinement_on && (_start_time <= _t && _t < _stop_time) && _step % _interval == 0;
 }
 
+void
+Adaptivity::switchHToPRefinement()
+{
+  _p_refinement_flag = true;
+  _fe_problem.havePRefinement();
+}
+
 #endif // LIBMESH_ENABLE_AMR
diff --git a/framework/src/base/Assembly.C b/framework/src/base/Assembly.C
index 6bb20c57e486..e142eda8649b 100644
--- a/framework/src/base/Assembly.C
+++ b/framework/src/base/Assembly.C
@@ -86,6 +86,13 @@ Assembly::Assembly(SystemBase & sys, THREAD_ID tid)
     _tid(tid),
     _mesh(sys.mesh()),
     _mesh_dimension(_mesh.dimension()),
+    _helper_type(_mesh.hasSecondOrderElements() ? SECOND : FIRST, LAGRANGE),
+    _user_added_fe_of_helper_type(false),
+    _user_added_fe_face_of_helper_type(false),
+    _user_added_fe_face_neighbor_of_helper_type(false),
+    _user_added_fe_neighbor_of_helper_type(false),
+    _user_added_fe_lower_of_helper_type(false),
+    _building_helpers(false),
     _current_qrule(nullptr),
     _current_qrule_volume(nullptr),
     _current_qrule_arbitrary(nullptr),
@@ -131,55 +138,39 @@ Assembly::Assembly(SystemBase & sys, THREAD_ID tid)
     _calculate_curvatures(false),
     _calculate_ad_coord(false)
 {
-  Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
+  const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
+  _building_helpers = true;
   // Build fe's for the helpers
   buildFE(FEType(helper_order, LAGRANGE));
   buildFaceFE(FEType(helper_order, LAGRANGE));
   buildNeighborFE(FEType(helper_order, LAGRANGE));
   buildFaceNeighborFE(FEType(helper_order, LAGRANGE));
   buildLowerDFE(FEType(helper_order, LAGRANGE));
+  _building_helpers = false;
 
   // Build an FE helper object for this type for each dimension up to the dimension of the current
   // mesh
   for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
   {
-    _holder_fe_helper[dim] = &_fe[dim][FEType(helper_order, LAGRANGE)];
-    (*_holder_fe_helper[dim])->get_phi();
-    (*_holder_fe_helper[dim])->get_dphi();
-    (*_holder_fe_helper[dim])->get_xyz();
-    (*_holder_fe_helper[dim])->get_JxW();
-
-    _holder_fe_face_helper[dim] = &_fe_face[dim][FEType(helper_order, LAGRANGE)];
-    (*_holder_fe_face_helper[dim])->get_phi();
-    (*_holder_fe_face_helper[dim])->get_dphi();
-    (*_holder_fe_face_helper[dim])->get_xyz();
-    (*_holder_fe_face_helper[dim])->get_JxW();
-    (*_holder_fe_face_helper[dim])->get_normals();
-
-    _holder_fe_face_neighbor_helper[dim] = &_fe_face_neighbor[dim][FEType(helper_order, LAGRANGE)];
-    (*_holder_fe_face_neighbor_helper[dim])->get_xyz();
-    (*_holder_fe_face_neighbor_helper[dim])->get_JxW();
-    (*_holder_fe_face_neighbor_helper[dim])->get_normals();
-
-    _holder_fe_neighbor_helper[dim] = &_fe_neighbor[dim][FEType(helper_order, LAGRANGE)];
-    (*_holder_fe_neighbor_helper[dim])->get_xyz();
-    (*_holder_fe_neighbor_helper[dim])->get_JxW();
+    _holder_fe_helper[dim] = _fe[dim][FEType(helper_order, LAGRANGE)];
+    _holder_fe_face_helper[dim] = _fe_face[dim][FEType(helper_order, LAGRANGE)];
+    _holder_fe_face_neighbor_helper[dim] = _fe_face_neighbor[dim][FEType(helper_order, LAGRANGE)];
+    _holder_fe_neighbor_helper[dim] = _fe_neighbor[dim][FEType(helper_order, LAGRANGE)];
   }
 
-  for (unsigned int dim = 0; dim <= _mesh_dimension - 1; dim++)
-  {
-    _holder_fe_lower_helper[dim] = &_fe_lower[dim][FEType(helper_order, LAGRANGE)];
-    // We need these computations in order to compute correct lower-d element volumes in curvilinear
-    // coordinates
-    (*_holder_fe_lower_helper[dim])->get_xyz();
-    (*_holder_fe_lower_helper[dim])->get_JxW();
-  }
+  for (unsigned int dim = 0; dim < _mesh_dimension; dim++)
+    _holder_fe_lower_helper[dim] = _fe_lower[dim][FEType(helper_order, LAGRANGE)];
+
+  // request phi, dphi, xyz, JxW, etc. data
+  helpersRequestData();
 
   // For 3D mortar, mortar segments are always TRI3 elements so we want FIRST LAGRANGE regardless
   // of discretization
   _fe_msm = (_mesh_dimension == 2)
                 ? FEGenericBase<Real>::build(_mesh_dimension - 1, FEType(helper_order, LAGRANGE))
                 : FEGenericBase<Real>::build(_mesh_dimension - 1, FEType(FIRST, LAGRANGE));
+  // This FE object should not take part in p-refinement
+  _fe_msm->add_p_level_in_reinit(false);
   _JxW_msm = &_fe_msm->get_JxW();
   // Prerequest xyz so that it is computed for _fe_msm so that it can be used for calculating
   // _coord_msm
@@ -267,6 +258,9 @@ Assembly::JxWNeighbor() const
 void
 Assembly::buildFE(FEType type) const
 {
+  if (!_building_helpers && type == _helper_type)
+    _user_added_fe_of_helper_type = true;
+
   if (!_fe_shape_data[type])
     _fe_shape_data[type] = std::make_unique<FEShapeData>();
 
@@ -290,6 +284,9 @@ Assembly::buildFE(FEType type) const
 void
 Assembly::buildFaceFE(FEType type) const
 {
+  if (!_building_helpers && type == _helper_type)
+    _user_added_fe_face_of_helper_type = true;
+
   if (!_fe_shape_data_face[type])
     _fe_shape_data_face[type] = std::make_unique<FEShapeData>();
 
@@ -309,6 +306,9 @@ Assembly::buildFaceFE(FEType type) const
 void
 Assembly::buildNeighborFE(FEType type) const
 {
+  if (!_building_helpers && type == _helper_type)
+    _user_added_fe_neighbor_of_helper_type = true;
+
   if (!_fe_shape_data_neighbor[type])
     _fe_shape_data_neighbor[type] = std::make_unique<FEShapeData>();
 
@@ -328,6 +328,9 @@ Assembly::buildNeighborFE(FEType type) const
 void
 Assembly::buildFaceNeighborFE(FEType type) const
 {
+  if (!_building_helpers && type == _helper_type)
+    _user_added_fe_face_neighbor_of_helper_type = true;
+
   if (!_fe_shape_data_face_neighbor[type])
     _fe_shape_data_face_neighbor[type] = std::make_unique<FEShapeData>();
 
@@ -347,6 +350,9 @@ Assembly::buildFaceNeighborFE(FEType type) const
 void
 Assembly::buildLowerDFE(FEType type) const
 {
+  if (!_building_helpers && type == _helper_type)
+    _user_added_fe_lower_of_helper_type = true;
+
   if (!_fe_shape_data_lower[type])
     _fe_shape_data_lower[type] = std::make_unique<FEShapeData>();
 
@@ -632,6 +638,11 @@ Assembly::setVolumeQRule(QBase * qrule, unsigned int dim)
       it.second->attach_quadrature_rule(qrule);
     for (auto & it : _vector_fe[dim])
       it.second->attach_quadrature_rule(qrule);
+    if (!_unique_fe_helper.empty())
+    {
+      mooseAssert(dim < _unique_fe_helper.size(), "We should not be indexing out of bounds");
+      _unique_fe_helper[dim]->attach_quadrature_rule(qrule);
+    }
   }
 }
 
@@ -644,6 +655,11 @@ Assembly::setFaceQRule(QBase * qrule, unsigned int dim)
     it.second->attach_quadrature_rule(qrule);
   for (auto & it : _vector_fe_face[dim])
     it.second->attach_quadrature_rule(qrule);
+  if (!_unique_fe_face_helper.empty())
+  {
+    mooseAssert(dim < _unique_fe_face_helper.size(), "We should not be indexing out of bounds");
+    _unique_fe_face_helper[dim]->attach_quadrature_rule(qrule);
+  }
 }
 
 void
@@ -658,6 +674,11 @@ Assembly::setLowerQRule(QBase * qrule, unsigned int dim)
     it.second->attach_quadrature_rule(qrule);
   for (auto & it : _vector_fe_lower[dim])
     it.second->attach_quadrature_rule(qrule);
+  if (!_unique_fe_lower_helper.empty())
+  {
+    mooseAssert(dim < _unique_fe_lower_helper.size(), "We should not be indexing out of bounds");
+    _unique_fe_lower_helper[dim]->attach_quadrature_rule(qrule);
+  }
 }
 
 void
@@ -669,6 +690,21 @@ Assembly::setNeighborQRule(QBase * qrule, unsigned int dim)
     it.second->attach_quadrature_rule(qrule);
   for (auto & it : _vector_fe_face_neighbor[dim])
     it.second->attach_quadrature_rule(qrule);
+  if (!_unique_fe_face_neighbor_helper.empty())
+  {
+    mooseAssert(dim < _unique_fe_face_neighbor_helper.size(),
+                "We should not be indexing out of bounds");
+    _unique_fe_face_neighbor_helper[dim]->attach_quadrature_rule(qrule);
+  }
+}
+
+void
+Assembly::clearCachedQRules()
+{
+  _current_qrule = nullptr;
+  _current_qrule_face = nullptr;
+  _current_qrule_lower = nullptr;
+  _current_qrule_neighbor = nullptr;
 }
 
 void
@@ -740,12 +776,17 @@ Assembly::reinitFE(const Elem * elem)
       fesd._curl_phi.shallowCopy(
           const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe.get_curl_phi()));
   }
+  if (!_unique_fe_helper.empty())
+  {
+    mooseAssert(dim < _unique_fe_helper.size(), "We should be in bounds here");
+    _unique_fe_helper[dim]->reinit(elem);
+  }
 
   // During that last loop the helper objects will have been reinitialized as well
   // We need to dig out the q_points and JxW from it.
   _current_q_points.shallowCopy(
-      const_cast<std::vector<Point> &>((*_holder_fe_helper[dim])->get_xyz()));
-  _current_JxW.shallowCopy(const_cast<std::vector<Real> &>((*_holder_fe_helper[dim])->get_JxW()));
+      const_cast<std::vector<Point> &>(_holder_fe_helper[dim]->get_xyz()));
+  _current_JxW.shallowCopy(const_cast<std::vector<Real> &>(_holder_fe_helper[dim]->get_JxW()));
 
   if (_subproblem.haveADObjects())
   {
@@ -755,7 +796,7 @@ Assembly::reinitFE(const Elem * elem)
     {
       const auto & qw = _current_qrule->get_weights();
       for (unsigned int qp = 0; qp != n_qp; qp++)
-        computeSinglePointMapAD(elem, qw, qp, *_holder_fe_helper[dim]);
+        computeSinglePointMapAD(elem, qw, qp, _holder_fe_helper[dim]);
     }
     else
       for (unsigned qp = 0; qp < n_qp; ++qp)
@@ -1235,15 +1276,20 @@ Assembly::reinitFEFace(const Elem * elem, unsigned int side)
       fesd._curl_phi.shallowCopy(
           const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
   }
+  if (!_unique_fe_face_helper.empty())
+  {
+    mooseAssert(dim < _unique_fe_face_helper.size(), "We should be in bounds here");
+    _unique_fe_face_helper[dim]->reinit(elem, side);
+  }
 
   // During that last loop the helper objects will have been reinitialized as well
   // We need to dig out the q_points and JxW from it.
   _current_q_points_face.shallowCopy(
-      const_cast<std::vector<Point> &>((*_holder_fe_face_helper[dim])->get_xyz()));
+      const_cast<std::vector<Point> &>(_holder_fe_face_helper[dim]->get_xyz()));
   _current_JxW_face.shallowCopy(
-      const_cast<std::vector<Real> &>((*_holder_fe_face_helper[dim])->get_JxW()));
+      const_cast<std::vector<Real> &>(_holder_fe_face_helper[dim]->get_JxW()));
   _current_normals.shallowCopy(
-      const_cast<std::vector<Point> &>((*_holder_fe_face_helper[dim])->get_normals()));
+      const_cast<std::vector<Point> &>(_holder_fe_face_helper[dim]->get_normals()));
 
   _mapped_normals.resize(_current_normals.size(), Eigen::Map<RealDIMValue>(nullptr));
   for (unsigned int i = 0; i < _current_normals.size(); i++)
@@ -1252,7 +1298,7 @@ Assembly::reinitFEFace(const Elem * elem, unsigned int side)
 
   if (_calculate_curvatures)
     _curvatures.shallowCopy(
-        const_cast<std::vector<Real> &>((*_holder_fe_face_helper[dim])->get_curvatures()));
+        const_cast<std::vector<Real> &>(_holder_fe_face_helper[dim]->get_curvatures()));
 
   computeADFace(*elem, side);
 
@@ -1277,9 +1323,9 @@ Assembly::computeFaceMap(const Elem & elem, const unsigned int side, const std::
   const Elem & side_elem = _compute_face_map_side_elem_builder(elem, side);
   const auto dim = elem.dim();
   const auto n_qp = qw.size();
-  const auto & dpsidxi_map = (*_holder_fe_face_helper[dim])->get_fe_map().get_dpsidxi();
-  const auto & dpsideta_map = (*_holder_fe_face_helper[dim])->get_fe_map().get_dpsideta();
-  const auto & psi_map = (*_holder_fe_face_helper[dim])->get_fe_map().get_psi();
+  const auto & dpsidxi_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsidxi();
+  const auto & dpsideta_map = _holder_fe_face_helper[dim]->get_fe_map().get_dpsideta();
+  const auto & psi_map = _holder_fe_face_helper[dim]->get_fe_map().get_psi();
   std::vector<std::vector<Real>> const * d2psidxi2_map = nullptr;
   std::vector<std::vector<Real>> const * d2psidxideta_map = nullptr;
   std::vector<std::vector<Real>> const * d2psideta2_map = nullptr;
@@ -1288,9 +1334,9 @@ Assembly::computeFaceMap(const Elem & elem, const unsigned int side, const std::
 
   if (_calculate_curvatures)
   {
-    d2psidxi2_map = &(*_holder_fe_face_helper[dim])->get_fe_map().get_d2psidxi2();
-    d2psidxideta_map = &(*_holder_fe_face_helper[dim])->get_fe_map().get_d2psidxideta();
-    d2psideta2_map = &(*_holder_fe_face_helper[dim])->get_fe_map().get_d2psideta2();
+    d2psidxi2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxi2();
+    d2psidxideta_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psidxideta();
+    d2psideta2_map = &_holder_fe_face_helper[dim]->get_fe_map().get_d2psideta2();
   }
 
   switch (dim)
@@ -1530,6 +1576,12 @@ Assembly::reinitFEFaceNeighbor(const Elem * neighbor, const std::vector<Point> &
       fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
           fe_face_neighbor.get_curl_phi()));
   }
+  if (!_unique_fe_face_neighbor_helper.empty())
+  {
+    mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(),
+                "We should be in bounds here");
+    _unique_fe_face_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points);
+  }
 }
 
 void
@@ -1577,6 +1629,11 @@ Assembly::reinitFENeighbor(const Elem * neighbor, const std::vector<Point> & ref
       fesd._curl_phi.shallowCopy(
           const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_neighbor.get_curl_phi()));
   }
+  if (!_unique_fe_neighbor_helper.empty())
+  {
+    mooseAssert(neighbor_dim < _unique_fe_neighbor_helper.size(), "We should be in bounds here");
+    _unique_fe_neighbor_helper[neighbor_dim]->reinit(neighbor, &reference_points);
+  }
 }
 
 void
@@ -1596,7 +1653,7 @@ Assembly::reinitNeighbor(const Elem * neighbor, const std::vector<Point> & refer
   if (_need_neighbor_elem_volume)
   {
     unsigned int dim = neighbor->dim();
-    FEBase & fe = **_holder_fe_neighbor_helper[dim];
+    FEBase & fe = *_holder_fe_neighbor_helper[dim];
     QBase * qrule = qrules(dim).vol.get();
 
     fe.attach_quadrature_rule(qrule);
@@ -1691,7 +1748,7 @@ Assembly::reinitAtPhysical(const Elem * elem, const std::vector<Point> & physica
   _current_elem_volume_computed = false;
 
   FEInterface::inverse_map(elem->dim(),
-                           (*_holder_fe_helper[elem->dim()])->get_fe_type(),
+                           _holder_fe_helper[elem->dim()]->get_fe_type(),
                            elem,
                            physical_points,
                            _temp_reference_points);
@@ -1702,6 +1759,16 @@ Assembly::reinitAtPhysical(const Elem * elem, const std::vector<Point> & physica
   _current_physical_points = physical_points;
 }
 
+void
+Assembly::setVolumeQRule(const Elem * const elem)
+{
+  unsigned int elem_dimension = elem->dim();
+  _current_qrule_volume = qrules(elem_dimension).vol.get();
+  // Make sure the qrule is the right one
+  if (_current_qrule != _current_qrule_volume)
+    setVolumeQRule(_current_qrule_volume, elem_dimension);
+}
+
 void
 Assembly::reinit(const Elem * elem)
 {
@@ -1711,14 +1778,7 @@ Assembly::reinit(const Elem * elem)
               "current subdomain has been set incorrectly");
   _current_elem_volume_computed = false;
 
-  unsigned int elem_dimension = elem->dim();
-
-  _current_qrule_volume = qrules(elem_dimension).vol.get();
-
-  // Make sure the qrule is the right one
-  if (_current_qrule != _current_qrule_volume)
-    setVolumeQRule(_current_qrule_volume, elem_dimension);
-
+  setVolumeQRule(elem);
   reinitFE(elem);
 
   computeCurrentElemVolume();
@@ -1805,7 +1865,17 @@ Assembly::qruleArbitraryFace(const Elem * elem, unsigned int side)
 }
 
 void
-Assembly::reinit(const Elem * elem, unsigned int side)
+Assembly::setFaceQRule(const Elem * const elem, const unsigned int side)
+{
+  const auto elem_dimension = elem->dim();
+  //// Make sure the qrule is the right one
+  auto rule = qruleFace(elem, side);
+  if (_current_qrule_face != rule)
+    setFaceQRule(rule, elem_dimension);
+}
+
+void
+Assembly::reinit(const Elem * const elem, const unsigned int side)
 {
   _current_elem = elem;
   _current_neighbor_elem = nullptr;
@@ -1815,15 +1885,9 @@ Assembly::reinit(const Elem * elem, unsigned int side)
   _current_elem_volume_computed = false;
   _current_side_volume_computed = false;
 
-  unsigned int elem_dimension = elem->dim();
-
   _current_side_elem = &_current_side_elem_builder(*elem, side);
 
-  //// Make sure the qrule is the right one
-  auto rule = qruleFace(elem, side);
-  if (_current_qrule_face != rule)
-    setFaceQRule(rule, elem_dimension);
-
+  setFaceQRule(elem, side);
   reinitFEFace(elem, side);
 
   computeCurrentFaceVolume();
@@ -1899,7 +1963,7 @@ Assembly::reinitElemAndNeighbor(const Elem * elem,
 
     // compute JxW on the neighbor's face
     _current_JxW_neighbor.shallowCopy(const_cast<std::vector<Real> &>(
-        (*_holder_fe_face_neighbor_helper[_current_neighbor_side_elem->dim()])->get_JxW()));
+        _holder_fe_face_neighbor_helper[_current_neighbor_side_elem->dim()]->get_JxW()));
   }
 
   reinitFEFaceNeighbor(neighbor, *reference_points_ptr);
@@ -1971,20 +2035,26 @@ Assembly::reinitElemFaceRef(const Elem * elem,
       fesd._curl_phi.shallowCopy(
           const_cast<std::vector<std::vector<VectorValue<Real>>> &>(fe_face.get_curl_phi()));
   }
+  if (!_unique_fe_face_helper.empty())
+  {
+    mooseAssert(elem_dim < _unique_fe_face_helper.size(), "We should be in bounds here");
+    _unique_fe_face_helper[elem_dim]->reinit(elem, elem_side, tolerance, pts, weights);
+  }
+
   // During that last loop the helper objects will have been reinitialized
   _current_q_points_face.shallowCopy(
-      const_cast<std::vector<Point> &>((*_holder_fe_face_helper[elem_dim])->get_xyz()));
+      const_cast<std::vector<Point> &>(_holder_fe_face_helper[elem_dim]->get_xyz()));
   _current_normals.shallowCopy(
-      const_cast<std::vector<Point> &>((*_holder_fe_face_helper[elem_dim])->get_normals()));
+      const_cast<std::vector<Point> &>(_holder_fe_face_helper[elem_dim]->get_normals()));
   _current_tangents.shallowCopy(const_cast<std::vector<std::vector<Point>> &>(
-      (*_holder_fe_face_helper[elem_dim])->get_tangents()));
+      _holder_fe_face_helper[elem_dim]->get_tangents()));
   // Note that if the user did pass in points and not weights to this method, JxW will be garbage
   // and should not be used
   _current_JxW_face.shallowCopy(
-      const_cast<std::vector<Real> &>((*_holder_fe_face_helper[elem_dim])->get_JxW()));
+      const_cast<std::vector<Real> &>(_holder_fe_face_helper[elem_dim]->get_JxW()));
   if (_calculate_curvatures)
     _curvatures.shallowCopy(
-        const_cast<std::vector<Real> &>((*_holder_fe_face_helper[elem_dim])->get_curvatures()));
+        const_cast<std::vector<Real> &>(_holder_fe_face_helper[elem_dim]->get_curvatures()));
 
   computeADFace(*elem, elem_side);
 }
@@ -2012,7 +2082,7 @@ Assembly::computeADFace(const Elem & elem, const unsigned int side)
       const std::vector<Real> dummy_qw(n_qp, 1.);
 
       for (unsigned int qp = 0; qp != n_qp; qp++)
-        computeSinglePointMapAD(&elem, dummy_qw, qp, *_holder_fe_face_helper[dim]);
+        computeSinglePointMapAD(&elem, dummy_qw, qp, _holder_fe_face_helper[dim]);
     }
     else
       for (unsigned qp = 0; qp < n_qp; ++qp)
@@ -2129,12 +2199,19 @@ Assembly::reinitNeighborFaceRef(const Elem * neighbor,
       fesd._curl_phi.shallowCopy(const_cast<std::vector<std::vector<VectorValue<Real>>> &>(
           fe_face_neighbor.get_curl_phi()));
   }
+  if (!_unique_fe_face_neighbor_helper.empty())
+  {
+    mooseAssert(neighbor_dim < _unique_fe_face_neighbor_helper.size(),
+                "We should be in bounds here");
+    _unique_fe_face_neighbor_helper[neighbor_dim]->reinit(
+        neighbor, neighbor_side, tolerance, pts, weights);
+  }
   // During that last loop the helper objects will have been reinitialized as well
   // We need to dig out the q_points from it
-  _current_q_points_face_neighbor.shallowCopy(const_cast<std::vector<Point> &>(
-      (*_holder_fe_face_neighbor_helper[neighbor_dim])->get_xyz()));
+  _current_q_points_face_neighbor.shallowCopy(
+      const_cast<std::vector<Point> &>(_holder_fe_face_neighbor_helper[neighbor_dim]->get_xyz()));
   _current_neighbor_normals.shallowCopy(const_cast<std::vector<Point> &>(
-      (*_holder_fe_face_neighbor_helper[neighbor_dim])->get_normals()));
+      _holder_fe_face_neighbor_helper[neighbor_dim]->get_normals()));
 }
 
 void
@@ -2210,6 +2287,11 @@ Assembly::reinitLowerDElem(const Elem * elem,
             const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower.get_dual_d2phi()));
     }
   }
+  if (!_unique_fe_lower_helper.empty())
+  {
+    mooseAssert(elem_dim < _unique_fe_lower_helper.size(), "We should be in bounds here");
+    _unique_fe_lower_helper[elem_dim]->reinit(elem);
+  }
 
   if (!_need_lower_d_elem_volume)
     return;
@@ -2231,7 +2313,7 @@ Assembly::reinitLowerDElem(const Elem * elem,
   else
   {
     // During that last loop the helper objects will have been reinitialized as well
-    FEBase & helper_fe = **_holder_fe_lower_helper[elem_dim];
+    FEBase & helper_fe = *_holder_fe_lower_helper[elem_dim];
     const auto & physical_q_points = helper_fe.get_xyz();
     const auto & JxW = helper_fe.get_JxW();
     MooseArray<Real> coord;
@@ -2296,7 +2378,7 @@ Assembly::reinitNeighborAtPhysical(const Elem * neighbor,
   // compute JxW on the neighbor's face
   unsigned int neighbor_side_dim = _current_neighbor_side_elem->dim();
   _current_JxW_neighbor.shallowCopy(const_cast<std::vector<Real> &>(
-      (*_holder_fe_face_neighbor_helper[neighbor_side_dim])->get_JxW()));
+      _holder_fe_face_neighbor_helper[neighbor_side_dim]->get_JxW()));
 
   reinitFEFaceNeighbor(neighbor, reference_points);
   reinitNeighbor(neighbor, reference_points);
@@ -4605,6 +4687,102 @@ Assembly::adCurvatures() const
   return _ad_curvatures;
 }
 
+void
+Assembly::helpersRequestData()
+{
+  for (unsigned int dim = 0; dim <= _mesh_dimension; dim++)
+  {
+    _holder_fe_helper[dim]->get_phi();
+    _holder_fe_helper[dim]->get_dphi();
+    _holder_fe_helper[dim]->get_xyz();
+    _holder_fe_helper[dim]->get_JxW();
+
+    _holder_fe_face_helper[dim]->get_phi();
+    _holder_fe_face_helper[dim]->get_dphi();
+    _holder_fe_face_helper[dim]->get_xyz();
+    _holder_fe_face_helper[dim]->get_JxW();
+    _holder_fe_face_helper[dim]->get_normals();
+
+    _holder_fe_face_neighbor_helper[dim]->get_xyz();
+    _holder_fe_face_neighbor_helper[dim]->get_JxW();
+    _holder_fe_face_neighbor_helper[dim]->get_normals();
+
+    _holder_fe_neighbor_helper[dim]->get_xyz();
+    _holder_fe_neighbor_helper[dim]->get_JxW();
+  }
+
+  for (unsigned int dim = 0; dim < _mesh_dimension; dim++)
+  {
+    // We need these computations in order to compute correct lower-d element volumes in curvilinear
+    // coordinates
+    _holder_fe_lower_helper[dim]->get_xyz();
+    _holder_fe_lower_helper[dim]->get_JxW();
+  }
+}
+
+void
+Assembly::havePRefinement()
+{
+  const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
+  const FEType helper_type(helper_order, LAGRANGE);
+  auto setup_helpers = [&helper_type](auto & unique_helper_container,
+                                      auto & helper_container,
+                                      const unsigned int num_dimensionalities,
+                                      const bool user_added_helper_type,
+                                      auto & fe_container)
+  {
+    unique_helper_container.resize(num_dimensionalities);
+    for (const auto dim : make_range(num_dimensionalities))
+    {
+      auto & unique_helper = unique_helper_container[dim];
+      unique_helper = FEGenericBase<Real>::build(dim, helper_type);
+      // don't participate in p-refinement
+      unique_helper->add_p_level_in_reinit(false);
+      helper_container[dim] = unique_helper.get();
+
+      // If the user did not request the helper type then we should erase it from our FE container
+      // so that they're not penalized (in the "we should be able to do p-refinement sense") for our
+      // perhaps silly helpers
+      if (!user_added_helper_type)
+      {
+        auto & fe_container_dim = fe_container[dim];
+        auto fe_it = fe_container_dim.find(helper_type);
+        mooseAssert(fe_it != fe_container_dim.end(), "We should have the helper type");
+        delete fe_it->second;
+        fe_container_dim.erase(fe_it);
+      }
+    }
+  };
+
+  setup_helpers(_unique_fe_helper,
+                _holder_fe_helper,
+                _mesh_dimension + 1,
+                _user_added_fe_of_helper_type,
+                _fe);
+  setup_helpers(_unique_fe_face_helper,
+                _holder_fe_face_helper,
+                _mesh_dimension + 1,
+                _user_added_fe_face_of_helper_type,
+                _fe_face);
+  setup_helpers(_unique_fe_face_neighbor_helper,
+                _holder_fe_face_neighbor_helper,
+                _mesh_dimension + 1,
+                _user_added_fe_face_neighbor_of_helper_type,
+                _fe_face_neighbor);
+  setup_helpers(_unique_fe_neighbor_helper,
+                _holder_fe_neighbor_helper,
+                _mesh_dimension + 1,
+                _user_added_fe_neighbor_of_helper_type,
+                _fe_neighbor);
+  setup_helpers(_unique_fe_lower_helper,
+                _holder_fe_lower_helper,
+                _mesh_dimension,
+                _user_added_fe_lower_of_helper_type,
+                _fe_lower);
+
+  helpersRequestData();
+}
+
 template void coordTransformFactor<Point, Real>(const SubProblem & s,
                                                 SubdomainID sub_id,
                                                 const Point & point,
diff --git a/framework/src/loops/MaxQpsThread.C b/framework/src/loops/MaxQpsThread.C
index 841523baacd9..2517a54dc10a 100644
--- a/framework/src/loops/MaxQpsThread.C
+++ b/framework/src/loops/MaxQpsThread.C
@@ -11,20 +11,17 @@
 #include "FEProblem.h"
 #include "Assembly.h"
 
-#include "libmesh/fe_base.h"
+#include "libmesh/fe.h"
 #include "libmesh/threads.h"
 #include LIBMESH_INCLUDE_UNORDERED_SET
 LIBMESH_DEFINE_HASH_POINTERS
 #include "libmesh/quadrature.h"
 
-MaxQpsThread::MaxQpsThread(FEProblemBase & fe_problem)
-  : _fe_problem(fe_problem), _max(0), _max_shape_funcs(0)
-{
-}
+MaxQpsThread::MaxQpsThread(FEProblemBase & fe_problem) : _fe_problem(fe_problem), _max(0) {}
 
 // Splitting Constructor
 MaxQpsThread::MaxQpsThread(MaxQpsThread & x, Threads::split /*split*/)
-  : _fe_problem(x._fe_problem), _max(x._max), _max_shape_funcs(x._max_shape_funcs)
+  : _fe_problem(x._fe_problem), _max(x._max)
 {
 }
 
@@ -50,45 +47,22 @@ MaxQpsThread::operator()(const ConstElemRange & range)
     // This ensures we can access the correct qrules if any block-specific
     // qrules have been created.
     assem.setCurrentSubdomainID(elem->subdomain_id());
+    assem.setVolumeQRule(elem);
 
-    FEType fe_type(FIRST, LAGRANGE);
-    unsigned int dim = elem->dim();
-    unsigned int side = 0; // we assume that any element will have at least one side ;)
-
-    // We cannot mess with the FE objects in Assembly, because we might need to request second
-    // derivatives
-    // later on. If we used them, we'd call reinit on them, thus making the call to request second
-    // derivatives harmful (i.e. leading to segfaults/asserts). Thus, we have to use a locally
-    // allocated object here.
-    //
-    // We'll use one for element interiors, which calculates nothing
-    std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
-    fe->get_nothing();
-
-    // And another for element sides, which calculates the minimum
-    // libMesh currently allows for that
-    std::unique_ptr<FEBase> side_fe(FEBase::build(dim, fe_type));
-    side_fe->get_xyz();
-
-    // figure out the number of qps for volume
-    auto qrule = assem.attachQRuleElem(dim, *fe);
-    fe->reinit(elem);
+    auto & qrule = assem.writeableQRule();
+    qrule->init(elem->type(), elem->p_level());
     if (qrule->n_points() > _max)
       _max = qrule->n_points();
 
-    unsigned int n_shape_funcs = fe->n_shape_functions();
-    if (n_shape_funcs > _max_shape_funcs)
-      _max_shape_funcs = n_shape_funcs;
-
-    // figure out the number of qps for the face
-    // NOTE: user might specify higher order rule for faces, thus possibly ending up with more qps
-    // than in the volume
-    auto qrule_face = assem.attachQRuleFace(dim, *side_fe);
-    if (dim > 0) // side reinit in 0D makes no sense, but we may have NodeElems
+    if (elem->dim())
     {
-      side_fe->reinit(elem, side);
+      // We are assuming here that side 0 ends up with at least as many quadrature points as any
+      // other side
+      assem.setFaceQRule(elem, /*side=*/0);
+      auto & qrule_face = assem.writeableQRuleFace();
+      qrule_face->init(elem->side_type(0), elem->p_level());
       if (qrule_face->n_points() > _max)
-        _max = qrule_face->n_points();
+        _max = qrule->n_points();
     }
 
     // In initial conditions nodes are enumerated as pretend quadrature points
@@ -97,6 +71,11 @@ MaxQpsThread::operator()(const ConstElemRange & range)
     if (elem->n_nodes() > _max)
       _max = elem->n_nodes();
   }
+
+  // Clear the cached quadrature rules because we may add FE objects in between now and simulation
+  // start and we only ensure we set all the FE quadrature rules if a quadrature rule is different
+  // from the cached quadrature rule
+  assem.clearCachedQRules();
 }
 
 void
@@ -104,7 +83,4 @@ MaxQpsThread::join(const MaxQpsThread & y)
 {
   if (y._max > _max)
     _max = y._max;
-
-  if (y._max_shape_funcs > _max_shape_funcs)
-    _max_shape_funcs = y._max_shape_funcs;
 }
diff --git a/framework/src/problems/DisplacedProblem.C b/framework/src/problems/DisplacedProblem.C
index 2c7a64ce91d1..bc9d3a4c9722 100644
--- a/framework/src/problems/DisplacedProblem.C
+++ b/framework/src/problems/DisplacedProblem.C
@@ -1303,3 +1303,11 @@ DisplacedProblem::currentResidualVectorTags() const
 {
   return _mproblem.currentResidualVectorTags();
 }
+
+void
+DisplacedProblem::havePRefinement()
+{
+  for (auto & assembly_vecs : _assembly)
+    for (auto & assembly : assembly_vecs)
+      assembly->havePRefinement();
+}
diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C
index f91c69f13a02..110b10344489 100644
--- a/framework/src/problems/FEProblemBase.C
+++ b/framework/src/problems/FEProblemBase.C
@@ -362,7 +362,6 @@ FEProblemBase::FEProblemBase(const InputParameters & parameters)
     _material_dependency_check(getParam<bool>("material_dependency_check")),
     _uo_aux_state_check(getParam<bool>("check_uo_aux_state")),
     _max_qps(std::numeric_limits<unsigned int>::max()),
-    _max_shape_funcs(std::numeric_limits<unsigned int>::max()),
     _max_scalar_order(INVALID_ORDER),
     _has_time_integrator(false),
     _has_exception(false),
@@ -1362,14 +1361,6 @@ FEProblemBase::getMaxQps() const
   return _max_qps;
 }
 
-unsigned int
-FEProblemBase::getMaxShapeFunctions() const
-{
-  if (_max_shape_funcs == std::numeric_limits<unsigned int>::max())
-    mooseError("Max shape functions uninitialized");
-  return _max_shape_funcs;
-}
-
 Order
 FEProblemBase::getMaxScalarOrder() const
 {
@@ -5241,13 +5232,11 @@ FEProblemBase::updateMaxQps()
     MaxQpsThread mqt(*this);
     Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), mqt);
     _max_qps = mqt.max();
-    _max_shape_funcs = mqt.max_shape_funcs();
 
     // If we have more shape functions or more quadrature points on
     // another processor, then we may need to handle those elements
     // ourselves later after repartitioning.
     _communicator.max(_max_qps);
-    _communicator.max(_max_shape_funcs);
   }
 
   unsigned int max_qpts = getMaxQps();
@@ -8095,3 +8084,14 @@ FEProblemBase::reinitMortarUserObjects(const BoundaryID primary_boundary_id,
     mortar_uo->reinit();
   }
 }
+
+void
+FEProblemBase::havePRefinement()
+{
+  for (auto & assembly_vecs : _assembly)
+    for (auto & assembly : assembly_vecs)
+      assembly->havePRefinement();
+
+  if (_displaced_problem)
+    _displaced_problem->havePRefinement();
+}
diff --git a/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg_test.i b/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg_test.i
index ba5eda6628e8..7c2238c39e42 100644
--- a/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg_test.i
+++ b/test/tests/dgkernels/2d_diffusion_dg/2d_diffusion_dg_test.i
@@ -28,30 +28,30 @@
 [Variables]
   active = 'u'
 
-  [./u]
+  [u]
     order = FIRST
     family = MONOMIAL
 
-    [./InitialCondition]
+    [InitialCondition]
       type = ConstantIC
       value = 1
-    [../]
-  [../]
+    []
+  []
 []
 
 [Functions]
   active = 'forcing_fn exact_fn'
 
-  [./forcing_fn]
+  [forcing_fn]
     type = ParsedFunction
 #    function = -4.0+(x*x)+(y*y)
 #    function = x
 #    function = (x*x)-2.0
     expression = 2*pow(e,-x-(y*y))*(1-2*y*y)
 #    function = (x*x*x)-6.0*x
-  [../]
+  []
 
-  [./exact_fn]
+  [exact_fn]
     type = ParsedGradFunction
 #    function = x
 #    grad_x = 1
@@ -65,58 +65,58 @@
 #    grad_x = 2*x
 #    grad_y = 0
 
-    value = pow(e,-x-(y*y))
+    expression = pow(e,-x-(y*y))
     grad_x = -pow(e,-x-(y*y))
     grad_y = -2*y*pow(e,-x-(y*y))
 
 #    function = (x*x*x)
 #    grad_x = 3*x*x
 #    grad_y = 0
-  [../]
+  []
 []
 
 [Kernels]
   active = 'diff abs forcing'
 
-  [./diff]
+  [diff]
     type = Diffusion
     variable = u
-  [../]
+  []
 
-  [./abs]          # u * v
+  [abs]          # u * v
     type = Reaction
     variable = u
-  [../]
+  []
 
-  [./forcing]
+  [forcing]
     type = BodyForce
     variable = u
     function = forcing_fn
-  [../]
+  []
 []
 
 [DGKernels]
   active = 'dg_diff'
 
-  [./dg_diff]
+  [dg_diff]
     type = DGDiffusion
     variable = u
     epsilon = -1
     sigma = 6
-  [../]
+  []
 []
 
 [BCs]
   active = 'all'
 
-  [./all]
+  [all]
     type = DGFunctionDiffusionDirichletBC
     variable = u
     boundary = '0 1 2 3'
     function = exact_fn
     epsilon = -1
     sigma = 6
-  [../]
+  []
 []
 
 [Executioner]
@@ -129,12 +129,12 @@
 
 #  petsc_options = '-snes_mf'
 #  max_r_steps = 2
-  [./Adaptivity]
+  [Adaptivity]
     steps = 2
     refine_fraction = 1.0
     coarsen_fraction = 0
     max_h_level = 8
-  [../]
+  []
 
   nl_rel_tol = 1e-10
 
@@ -144,19 +144,19 @@
 [Postprocessors]
   active = 'h dofs l2_err'
 
-  [./h]
+  [h]
     type = AverageElementSize
-  [../]
+  []
 
-  [./dofs]
+  [dofs]
     type = NumDOFs
-  [../]
+  []
 
-  [./l2_err]
+  [l2_err]
     type = ElementL2Error
     variable = u
     function = exact_fn
-  [../]
+  []
 []
 
 [Outputs]
diff --git a/test/tests/dgkernels/2d_diffusion_dg/gold/p_refinement.e-s003 b/test/tests/dgkernels/2d_diffusion_dg/gold/p_refinement.e-s003
new file mode 100644
index 000000000000..f9fb8f769d70
Binary files /dev/null and b/test/tests/dgkernels/2d_diffusion_dg/gold/p_refinement.e-s003 differ
diff --git a/test/tests/dgkernels/2d_diffusion_dg/tests b/test/tests/dgkernels/2d_diffusion_dg/tests
index 8c0b8c9d85f7..4db589f741d7 100644
--- a/test/tests/dgkernels/2d_diffusion_dg/tests
+++ b/test/tests/dgkernels/2d_diffusion_dg/tests
@@ -5,16 +5,23 @@
     type = 'Exodiff'
     input = '2d_diffusion_dg_test.i'
     exodiff = 'out.e-s003'
-    max_parallel = 1
     requirement = "The system shall support solving 2D diffusion using the discontinuous Galerkin "
                   "method."
   []
+  [p_refinement_test]
+    type = 'Exodiff'
+    input = '2d_diffusion_dg_test.i'
+    cli_args = 'Outputs/file_base=p_refinement Executioner/Adaptivity/switch_h_to_p_refinement=true'
+    exodiff = 'p_refinement.e-s003'
+    requirement = "The system shall support solving 2D diffusion using the discontinuous Galerkin "
+                  "method with p-refinement."
+    issues = '#869 #24141'
+  []
   [constMonomial]
     type = 'Exodiff'
     input = '2d_diffusion_dg_test.i'
     exodiff = 'constMonomial.e-s003'
     cli_args = 'Outputs/file_base=constMonomial Variables/u/order=CONSTANT'
-    max_parallel = 1
     issues = '#21314'
     requirement = "The system shall support solving 2D diffusion using the discontinuous Galerkin "
                   "method with constant monomial shape functions."
diff --git a/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_dg_test.i b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_dg_test.i
index a0bdff3c19d8..bc34bf5f90e3 100644
--- a/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_dg_test.i
+++ b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_dg_test.i
@@ -16,76 +16,76 @@
 [Variables]
   active = 'u'
 
-  [./u]
+  [u]
     order = FIRST
     family = MONOMIAL
 
-    [./InitialCondition]
+    [InitialCondition]
       type = ConstantIC
       value = 0.5
-    [../]
-  [../]
+    []
+  []
 []
 
 [Functions]
   active = 'forcing_fn exact_fn'
 
-  [./forcing_fn]
+  [forcing_fn]
     type = ParsedFunction
     expression = 2*pow(e,-x-(y*y))*(1-2*y*y)
-  [../]
+  []
 
-  [./exact_fn]
+  [exact_fn]
     type = ParsedGradFunction
 
-    value = pow(e,-x-(y*y))
+    expression = pow(e,-x-(y*y))
     grad_x = -pow(e,-x-(y*y))
     grad_y = -2*y*pow(e,-x-(y*y))
-  [../]
+  []
 []
 
 [Kernels]
   active = 'diff abs forcing'
 
-  [./diff]
+  [diff]
     type = Diffusion
     variable = u
-  [../]
+  []
 
-  [./abs]          # u * v
+  [abs]          # u * v
     type = Reaction
     variable = u
-  [../]
+  []
 
-  [./forcing]
+  [forcing]
     type = BodyForce
     variable = u
     function = forcing_fn
-  [../]
+  []
 []
 
 [DGKernels]
   active = 'dg_diff'
 
-  [./dg_diff]
+  [dg_diff]
     type = DGDiffusion
     variable = u
     epsilon = -1
     sigma = 6
-  [../]
+  []
 []
 
 [BCs]
   active = 'all'
 
-  [./all]
+  [all]
     type = DGFunctionDiffusionDirichletBC
     variable = u
     boundary = '0 1 2 3 4 5'
     function = exact_fn
     epsilon = -1
     sigma = 6
-  [../]
+  []
 []
 
 [Executioner]
@@ -97,22 +97,22 @@
 [Postprocessors]
   active = 'h dofs l2_err'
 
-  [./h]
+  [h]
     type = AverageElementSize
     execute_on = 'initial timestep_end'
-  [../]
+  []
 
-  [./dofs]
+  [dofs]
     type = NumDOFs
     execute_on = 'initial timestep_end'
-  [../]
+  []
 
-  [./l2_err]
+  [l2_err]
     type = ElementL2Error
     variable = u
     function = exact_fn
     execute_on = 'initial timestep_end'
-  [../]
+  []
 []
 
 [Outputs]
diff --git a/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i
new file mode 100644
index 000000000000..83e3bd21ea5f
--- /dev/null
+++ b/test/tests/dgkernels/3d_diffusion_dg/3d_diffusion_p_refinement.i
@@ -0,0 +1,107 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 3
+  nx = 1
+  ny = 1
+  nz = 1
+  xmin = 0
+  xmax = 1
+  ymin = 0
+  ymax = 1
+  zmin = 0
+  zmax = 1
+  elem_type = HEX8
+[]
+
+[Variables]
+  [u]
+    order = FIRST
+    family = MONOMIAL
+    [InitialCondition]
+      type = ConstantIC
+      value = 0.5
+    []
+  []
+[]
+
+[Functions]
+  [forcing_fn]
+    type = ParsedFunction
+    expression = 2*pow(e,-x-(y*y))*(1-2*y*y)
+  []
+  [exact_fn]
+    type = ParsedGradFunction
+    expression = pow(e,-x-(y*y))
+    grad_x = -pow(e,-x-(y*y))
+    grad_y = -2*y*pow(e,-x-(y*y))
+  []
+[]
+
+[Kernels]
+  [diff]
+    type = Diffusion
+    variable = u
+  []
+  [abs]
+    type = Reaction
+    variable = u
+  []
+  [forcing]
+    type = BodyForce
+    variable = u
+    function = forcing_fn
+  []
+[]
+
+[DGKernels]
+  [dg_diff]
+    type = DGDiffusion
+    variable = u
+    epsilon = -1
+    sigma = 6
+  []
+[]
+
+[BCs]
+  [all]
+    type = DGFunctionDiffusionDirichletBC
+    variable = u
+    boundary = '0 1 2 3 4 5'
+    function = exact_fn
+    epsilon = -1
+    sigma = 6
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = 'PJFNK'
+  [Adaptivity]
+    switch_h_to_p_refinement = true
+    steps = 2
+    refine_fraction = 1.0
+    coarsen_fraction = 0
+    max_h_level = 8
+  []
+[]
+
+[Postprocessors]
+  [h]
+    type = AverageElementSize
+    execute_on = 'initial timestep_end'
+  []
+  [dofs]
+    type = NumDOFs
+    execute_on = 'initial timestep_end'
+  []
+  [l2_err]
+    type = ElementL2Error
+    variable = u
+    function = exact_fn
+    execute_on = 'initial timestep_end'
+  []
+[]
+
+[Outputs]
+  exodus = true
+[]
diff --git a/test/tests/dgkernels/3d_diffusion_dg/gold/3d_diffusion_p_refinement_out.e-s003 b/test/tests/dgkernels/3d_diffusion_dg/gold/3d_diffusion_p_refinement_out.e-s003
new file mode 100644
index 000000000000..cfd26d76da41
Binary files /dev/null and b/test/tests/dgkernels/3d_diffusion_dg/gold/3d_diffusion_p_refinement_out.e-s003 differ
diff --git a/test/tests/dgkernels/3d_diffusion_dg/tests b/test/tests/dgkernels/3d_diffusion_dg/tests
index 7a1206169c62..1c7a1b9c80b9 100644
--- a/test/tests/dgkernels/3d_diffusion_dg/tests
+++ b/test/tests/dgkernels/3d_diffusion_dg/tests
@@ -1,13 +1,20 @@
 [Tests]
-  issues = '#869'
   design = 'source/bcs/DGFunctionDiffusionDirichletBC.md source/dgkernels/DGDiffusion.md'
-  [./test]
+  [test]
     type = 'Exodiff'
     input = '3d_diffusion_dg_test.i'
     exodiff = 'out.e'
     group = 'adaptive'
-    max_parallel = 1
     valgrind = 'HEAVY'
     requirement = "The system shall support solving 3D diffusion using the discontinuous Galerkin method."
-  [../]
+    issues = '#869'
+  []
+  [p_refinement_test]
+    type = 'Exodiff'
+    input = '3d_diffusion_p_refinement.i'
+    exodiff = '3d_diffusion_p_refinement_out.e-s003'
+    requirement = "The system shall support solving 3D diffusion using the discontinuous Galerkin "
+                  "method with p-refinement."
+    issues = '#869 #24141'
+  []
 []