From 347b81e4eb5ea54a51c7335a0f406b6d365066ee Mon Sep 17 00:00:00 2001
From: Daniel Schwen <daniel@schwen.de>
Date: Sat, 18 Nov 2023 16:34:44 -0700
Subject: [PATCH] New approach (#26053)

---
 .../ProjectedStatefulMaterialStorageAction.h  | 156 +++++++++++---
 .../auxkernels/NodalPatchRecoveryAux.h        |  13 +-
 .../auxkernels/NodalPatchRecoveryAuxBase.h    |  36 ++++
 ...tedMaterialPropertyNodalPatchRecoveryAux.h |  30 +++
 .../auxkernels/ProjectedStatefulMaterialAux.h |  23 +-
 framework/include/base/Registry.h             |   6 +-
 .../materials/InterpolatedStatefulMaterial.h  |  47 ++---
 ...jectedStatefulMaterialNodalPatchRecovery.h | 104 +++++++++-
 framework/include/utils/SerialAccess.h        | 110 ++++++----
 .../include/utils_nonunity/RankFourTensor.h   |  12 --
 .../include/utils_nonunity/RankTwoTensor.h    |  15 --
 .../ProjectedStatefulMaterialStorageAction.C  | 164 +--------------
 .../src/auxkernels/NodalPatchRecoveryAux.C    |  63 +-----
 .../auxkernels/NodalPatchRecoveryAuxBase.C    |  79 +++++++
 ...tedMaterialPropertyNodalPatchRecoveryAux.C |  40 ++++
 .../auxkernels/ProjectedStatefulMaterialAux.C |  51 +++--
 .../materials/InterpolatedStatefulMaterial.C  | 112 +++-------
 ...jectedStatefulMaterialNodalPatchRecovery.C | 196 +++++++++++++++++-
 .../projected_stateful_materials/test.i       |   4 +
 unit/src/MooseServerTest.C                    |   3 +-
 unit/src/SerialAccessTest.C                   |  36 +++-
 21 files changed, 812 insertions(+), 488 deletions(-)
 create mode 100644 framework/include/auxkernels/NodalPatchRecoveryAuxBase.h
 create mode 100644 framework/include/auxkernels/ProjectedMaterialPropertyNodalPatchRecoveryAux.h
 create mode 100644 framework/src/auxkernels/NodalPatchRecoveryAuxBase.C
 create mode 100644 framework/src/auxkernels/ProjectedMaterialPropertyNodalPatchRecoveryAux.C

diff --git a/framework/include/actions/ProjectedStatefulMaterialStorageAction.h b/framework/include/actions/ProjectedStatefulMaterialStorageAction.h
index 753cdc37fa5e..21e658f1ad0d 100644
--- a/framework/include/actions/ProjectedStatefulMaterialStorageAction.h
+++ b/framework/include/actions/ProjectedStatefulMaterialStorageAction.h
@@ -10,8 +10,18 @@
 #pragma once
 
 #include "Action.h"
+#include "Registry.h"
+#include "Conversion.h"
+#include "SerialAccess.h"
+#include "RankTwoTensorForward.h"
+#include "RankFourTensorForward.h"
 #include "libmesh/fe_type.h"
 
+// created types
+#include <InterpolatedStatefulMaterial.h>
+#include <ProjectedStatefulMaterialAux.h>
+#include <ProjectedStatefulMaterialNodalPatchRecovery.h>
+
 /**
  * Set up AuxKernels and AuxVariables for projected material property storage (PoMPS).
  */
@@ -30,38 +40,26 @@ class ProjectedStatefulMaterialStorageAction : public Action
   /// List of supported types
   typedef Moose::TypeList<Real, RealVectorValue, RankTwoTensor, RankFourTensor> SupportedTypes;
 
+  /// get an enum with all supported types
+  static MooseEnum getTypeEnum();
+
+  template <typename T, bool is_ad>
+  void processProperty(const MaterialPropertyName & prop_name);
+
 protected:
-  /**
-   * Perform setup for a single scalar component of the material property prop_name, and gather the
-   * names of teh AuxVariables used to represent each component in `vars`.
-   */
-  void processComponent(const std::string & prop_name,
-                        std::vector<unsigned int> idx,
-                        std::vector<VariableName> & vars,
-                        bool is_ad);
-
-  /**
-   * Add the material object to obtain the interpolated old state (for use with
-   * InterpolatedStatefulMaterialPropertyInterface)
-   */
-  void addMaterial(const std::string & prop_type,
-                   const std::string & prop_name,
-                   std::vector<VariableName> & vars);
-
-  enum class PropertyType
+  template <typename T, int I>
+  struct ProcessProperty
   {
-    REAL,
-    REALVECTORVALUE,
-    RANKTWOTENSOR,
-    RANKFOURTENSOR
+    static void apply(ProjectedStatefulMaterialStorageAction * self,
+                      const MaterialPropertyName & prop_name)
+    {
+      self->processProperty<T, false>(prop_name);
+      self->processProperty<T, true>(prop_name);
+    }
   };
-  typedef std::pair<PropertyType, bool> PropertyInfo;
 
-  /**
-   * Return the property type and whether to use AD or not. If no property with a supported type is
-   * found, throw an error.
-   */
-  PropertyInfo checkProperty(const std::string & prop_name);
+  template <typename... Ts>
+  static MooseEnum getTypeEnum(typename Moose::TypeList<Ts...>);
 
   const std::vector<MaterialPropertyName> & _prop_names;
 
@@ -70,3 +68,105 @@ class ProjectedStatefulMaterialStorageAction : public Action
   const std::string _var_type;
   const std::string _pomps_prefix;
 };
+
+template <typename... Ts>
+MooseEnum
+ProjectedStatefulMaterialStorageAction::getTypeEnum(Moose::TypeList<Ts...>)
+{
+  return MooseEnum(Moose::stringify(std::vector<std::string>{typeid(Ts).name()...}, " "));
+}
+
+template <typename T, bool is_ad>
+void
+ProjectedStatefulMaterialStorageAction::processProperty(const MaterialPropertyName & prop_name)
+{
+  // check if a property of type T exists
+  const auto & data = _problem->getMaterialData(Moose::BLOCK_MATERIAL_DATA);
+  if (!data.haveGenericProperty<T, is_ad>(prop_name))
+    return;
+
+  // number of scalar components
+  const auto size = Moose::SerialAccess<T>::size();
+
+  // generate variable names
+  std::vector<VariableName> vars;
+  for (const auto i : make_range(size))
+    vars.push_back("_var_" + prop_name + '_' + Moose::stringify(i));
+
+  if (_current_task == "setup_projected_properties")
+  {
+    // add the AuxVars for storage
+    for (const auto & var : vars)
+    {
+      auto params = _factory.getValidParams(_var_type);
+      params.applyParameters(parameters());
+      params.set<std::vector<OutputName>>("outputs") = {"none"};
+      _problem->addAuxVariable(_var_type, var, params);
+    }
+
+    // add material
+    {
+      const auto type = Registry::getClassName<InterpolatedStatefulMaterialTempl<T>>();
+      auto params = _factory.getValidParams(type);
+      params.applySpecificParameters(parameters(), {"block"});
+      params.template set<std::vector<VariableName>>("old_state") = vars;
+      params.template set<MaterialPropertyName>("prop_name") = prop_name;
+      _problem->addMaterial(type, "_mat_" + prop_name, params);
+    }
+
+    // use nodal patch recovery for lagrange
+    if (_fe_type.family == LAGRANGE)
+    {
+      // nodal variables require patch recovery (add user object)
+      const auto & type =
+          Registry::getClassName<ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>>();
+      auto params = _factory.getValidParams(type);
+      params.applySpecificParameters(parameters(), {"block"});
+      params.template set<MaterialPropertyName>("property") = prop_name;
+      params.template set<MooseEnum>("patch_polynomial_order") = _order;
+      params.template set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
+      params.template set<bool>("force_preaux") = true;
+      _problem->addUserObject(type, "_npruo_" + prop_name, params);
+    }
+  }
+
+  if (_current_task == "add_aux_kernel")
+  {
+    // create variables
+    std::vector<std::string> auxnames;
+    for (const auto i : make_range(size))
+      auxnames.push_back("_aux_" + prop_name + '_' + Moose::stringify(i));
+
+    // use nodal patch recovery for lagrange
+    if (_fe_type.family == LAGRANGE)
+    {
+      // nodal variables require patch recovery (add aux kernel)
+      const auto & type = "ProjectedMaterialPropertyNodalPatchRecoveryAux";
+      for (const auto i : make_range(size))
+      {
+        auto params = _factory.getValidParams(type);
+        params.applySpecificParameters(parameters(), {"block"});
+        params.template set<AuxVariableName>("variable") = vars[i];
+        params.template set<unsigned int>("component") = i;
+        params.template set<UserObjectName>("nodal_patch_recovery_uo") = "_npruo_" + prop_name;
+        params.template set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
+        _problem->addAuxKernel(type, auxnames[i], params);
+      }
+    }
+    else
+    {
+      // elemental variables
+      const auto & type = Registry::getClassName<ProjectedStatefulMaterialAuxTempl<T, is_ad>>();
+      for (const auto i : make_range(size))
+      {
+        auto params = _factory.getValidParams(type);
+        params.applySpecificParameters(parameters(), {"block"});
+        params.template set<AuxVariableName>("variable") = vars[i];
+        params.template set<unsigned int>("component") = i;
+        params.template set<MaterialPropertyName>("prop") = prop_name;
+        params.template set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
+        _problem->addAuxKernel(type, auxnames[i], params);
+      }
+    }
+  }
+}
diff --git a/framework/include/auxkernels/NodalPatchRecoveryAux.h b/framework/include/auxkernels/NodalPatchRecoveryAux.h
index 4a99886e784f..5d089d2d8c84 100644
--- a/framework/include/auxkernels/NodalPatchRecoveryAux.h
+++ b/framework/include/auxkernels/NodalPatchRecoveryAux.h
@@ -9,25 +9,20 @@
 
 #pragma once
 
-#include "AuxKernel.h"
+#include "NodalPatchRecoveryAuxBase.h"
 
+/// Forward declare user object
 class NodalPatchRecoveryBase;
 
-class NodalPatchRecoveryAux : public AuxKernel
+class NodalPatchRecoveryAux : public NodalPatchRecoveryAuxBase
 {
 public:
   static InputParameters validParams();
 
   NodalPatchRecoveryAux(const InputParameters & parameters);
 
-  /**
-   * Block restrict elements on which to perform the variable/property nodal recovery.
-   */
-  void blockRestrictElements(std::vector<dof_id_type> & elem_ids,
-                             const std::vector<dof_id_type> & node_to_elem_pair_elems) const;
-
 protected:
-  virtual Real computeValue() override;
+  virtual Real nodalPatchRecovery() override;
 
 private:
   const NodalPatchRecoveryBase & _npr;
diff --git a/framework/include/auxkernels/NodalPatchRecoveryAuxBase.h b/framework/include/auxkernels/NodalPatchRecoveryAuxBase.h
new file mode 100644
index 000000000000..f24e1f12cf2b
--- /dev/null
+++ b/framework/include/auxkernels/NodalPatchRecoveryAuxBase.h
@@ -0,0 +1,36 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "AuxKernel.h"
+
+class NodalPatchRecoveryBase;
+
+class NodalPatchRecoveryAuxBase : public AuxKernel
+{
+public:
+  static InputParameters validParams();
+
+  NodalPatchRecoveryAuxBase(const InputParameters & parameters);
+
+  /**
+   * Block restrict elements on which to perform the variable/property nodal recovery.
+   */
+  void blockRestrictElements(std::vector<dof_id_type> & elem_ids,
+                             const std::vector<dof_id_type> & node_to_elem_pair_elems) const;
+
+protected:
+  virtual Real computeValue() override;
+
+  /// Override this to get the fitted value form a Nodal Patch Recovery User Object
+  virtual Real nodalPatchRecovery() = 0;
+
+  std::vector<dof_id_type> _elem_ids;
+};
diff --git a/framework/include/auxkernels/ProjectedMaterialPropertyNodalPatchRecoveryAux.h b/framework/include/auxkernels/ProjectedMaterialPropertyNodalPatchRecoveryAux.h
new file mode 100644
index 000000000000..3e3fa25f1fdc
--- /dev/null
+++ b/framework/include/auxkernels/ProjectedMaterialPropertyNodalPatchRecoveryAux.h
@@ -0,0 +1,30 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "NodalPatchRecoveryAuxBase.h"
+
+/// Forward declare user object
+class ProjectedStatefulMaterialNodalPatchRecoveryBase;
+
+class ProjectedMaterialPropertyNodalPatchRecoveryAux : public NodalPatchRecoveryAuxBase
+{
+public:
+  static InputParameters validParams();
+
+  ProjectedMaterialPropertyNodalPatchRecoveryAux(const InputParameters & parameters);
+
+protected:
+  virtual Real nodalPatchRecovery() override;
+
+private:
+  const ProjectedStatefulMaterialNodalPatchRecoveryBase & _npr;
+  const unsigned int _component;
+};
diff --git a/framework/include/auxkernels/ProjectedStatefulMaterialAux.h b/framework/include/auxkernels/ProjectedStatefulMaterialAux.h
index 5e4c56de9e8d..66b321208f3a 100644
--- a/framework/include/auxkernels/ProjectedStatefulMaterialAux.h
+++ b/framework/include/auxkernels/ProjectedStatefulMaterialAux.h
@@ -16,7 +16,7 @@
 
 class MaterialBase;
 
-template <bool is_ad>
+template <typename T, bool is_ad>
 class ProjectedStatefulMaterialAuxTempl : public AuxKernel
 {
 public:
@@ -29,11 +29,26 @@ class ProjectedStatefulMaterialAuxTempl : public AuxKernel
 protected:
   virtual Real computeValue() override;
 
-  const IndexableProperty<AuxKernel, is_ad> _prop;
   const SubdomainID & _current_subdomain_id;
 
   std::map<SubdomainID, std::vector<MaterialBase *>> _required_materials;
+
+  const GenericMaterialProperty<T, is_ad> & _prop;
+
+  const unsigned int _component;
 };
 
-typedef ProjectedStatefulMaterialAuxTempl<false> ProjectedStatefulMaterialAux;
-typedef ProjectedStatefulMaterialAuxTempl<true> ADProjectedStatefulMaterialAux;
+typedef ProjectedStatefulMaterialAuxTempl<Real, false> ProjectedStatefulMaterialRealAux;
+typedef ProjectedStatefulMaterialAuxTempl<Real, true> ADProjectedStatefulMaterialRealAux;
+typedef ProjectedStatefulMaterialAuxTempl<RealVectorValue, false>
+    ProjectedStatefulMaterialRealVectorValueAux;
+typedef ProjectedStatefulMaterialAuxTempl<RealVectorValue, true>
+    ADProjectedStatefulMaterialRealVectorValueAux;
+typedef ProjectedStatefulMaterialAuxTempl<RankTwoTensor, false>
+    ProjectedStatefulMaterialRankTwoTensorAux;
+typedef ProjectedStatefulMaterialAuxTempl<RankTwoTensor, true>
+    ADProjectedStatefulMaterialRankTwoTensorAux;
+typedef ProjectedStatefulMaterialAuxTempl<RankFourTensor, false>
+    ProjectedStatefulMaterialRankFourTensorAux;
+typedef ProjectedStatefulMaterialAuxTempl<RankFourTensor, true>
+    ADProjectedStatefulMaterialRankFourTensorAux;
diff --git a/framework/include/base/Registry.h b/framework/include/base/Registry.h
index 5ea0027c31a3..c8143ee04179 100644
--- a/framework/include/base/Registry.h
+++ b/framework/include/base/Registry.h
@@ -253,10 +253,8 @@ template <typename T>
 std::string
 Registry::getRegisteredName()
 {
-  for (const auto & [name, reg_ptr] : getRegistry()._name_to_entry)
-    if (std::dynamic_pointer_cast<T>(reg_ptr))
-      return name;
-  mooseError("The object of C++ type '", demangle(typeid(T).name()), "' has not registered.");
+  mooseDeprecated("Use Registry::getClassName() instead.");
+  return getClassName<T>();
 }
 
 template <typename T>
diff --git a/framework/include/materials/InterpolatedStatefulMaterial.h b/framework/include/materials/InterpolatedStatefulMaterial.h
index f67749ef00a7..dee50c281e81 100644
--- a/framework/include/materials/InterpolatedStatefulMaterial.h
+++ b/framework/include/materials/InterpolatedStatefulMaterial.h
@@ -14,16 +14,16 @@
 #include "RankFourTensorForward.h"
 
 /**
- * PowerLawSoftening is a smeared crack softening model that
- * uses a power law equation to soften the tensile response.
- * It is for use with ComputeSmearedCrackingStress.
+ * Reconstitute a materal property from the old and older states of projected AuxVariables. Use
+ * though the ProjectedStatefulMaterialStorageAction.
  */
-class InterpolatedStatefulMaterial : public Material
+template <typename T>
+class InterpolatedStatefulMaterialTempl : public Material
 {
 public:
   static InputParameters validParams();
 
-  InterpolatedStatefulMaterial(const InputParameters & parameters);
+  InterpolatedStatefulMaterialTempl(const InputParameters & parameters);
 
   virtual void computeQpProperties() override;
 
@@ -34,29 +34,22 @@ class InterpolatedStatefulMaterial : public Material
   /// Older projected state
   const std::vector<const VariableValue *> _older_state;
 
+  /// total number of components
+  const std::size_t _size;
+
   /// emitted property name
   const MaterialPropertyName _prop_name;
 
-  /// Property type
-  enum class PropType
-  {
-    REAL,
-    REALVECTORVALUE,
-    RANKTWOTENSOR,
-    RANKFOURTENSOR
-  } _prop_type;
-
-  ///@{ Old interpolated properties
-  MaterialProperty<Real> * _prop_old_real;
-  MaterialProperty<RealVectorValue> * _prop_old_realvectorvalue;
-  MaterialProperty<RankTwoTensor> * _prop_old_ranktwotensor;
-  MaterialProperty<RankFourTensor> * _prop_old_rankfourtensor;
-  ///@}
-
-  ///@{ Older interpolated properties
-  MaterialProperty<Real> * _prop_older_real;
-  MaterialProperty<RealVectorValue> * _prop_older_realvectorvalue;
-  MaterialProperty<RankTwoTensor> * _prop_older_ranktwotensor;
-  MaterialProperty<RankFourTensor> * _prop_older_rankfourtensor;
-  ///@}
+  /// Old interpolated property
+  MaterialProperty<T> & _prop_old;
+
+  /// Older interpolated properties
+  MaterialProperty<T> & _prop_older;
 };
+
+typedef InterpolatedStatefulMaterialTempl<Real> InterpolatedStatefulMaterialReal;
+typedef InterpolatedStatefulMaterialTempl<RealVectorValue>
+    InterpolatedStatefulMaterialRealVectorValue;
+typedef InterpolatedStatefulMaterialTempl<RankTwoTensor> InterpolatedStatefulMaterialRankTwoTensor;
+typedef InterpolatedStatefulMaterialTempl<RankFourTensor>
+    InterpolatedStatefulMaterialRankFourTensor;
diff --git a/framework/include/userobjects/ProjectedStatefulMaterialNodalPatchRecovery.h b/framework/include/userobjects/ProjectedStatefulMaterialNodalPatchRecovery.h
index 6e6688ac1c09..483c11bc4440 100644
--- a/framework/include/userobjects/ProjectedStatefulMaterialNodalPatchRecovery.h
+++ b/framework/include/userobjects/ProjectedStatefulMaterialNodalPatchRecovery.h
@@ -10,28 +10,112 @@
 #pragma once
 
 // MOOSE includes
-#include "NodalPatchRecoveryMaterialProperty.h"
+#include "ElementUserObject.h"
+#include "SerialAccess.h"
 
-#include <deque>
+class ProjectedStatefulMaterialNodalPatchRecoveryBase : public ElementUserObject
+{
+public:
+  static InputParameters validParams();
+
+  ProjectedStatefulMaterialNodalPatchRecoveryBase(const InputParameters & parameters)
+    : ElementUserObject(parameters)
+  {
+  }
 
-class MaterialBase;
+  virtual Real nodalPatchRecovery(const Point & p,
+                                  const std::vector<dof_id_type> & elem_ids,
+                                  std::size_t component) const = 0;
+};
 
-/**
- * Nodal patch recovery for material property components for projected stateful properties
- */
-class ProjectedStatefulMaterialNodalPatchRecovery : public NodalPatchRecoveryMaterialProperty
+template <typename T, bool is_ad>
+class ProjectedStatefulMaterialNodalPatchRecoveryTempl
+  : public ProjectedStatefulMaterialNodalPatchRecoveryBase
 {
 public:
   static InputParameters validParams();
 
-  ProjectedStatefulMaterialNodalPatchRecovery(const InputParameters & parameters);
+  ProjectedStatefulMaterialNodalPatchRecoveryTempl(const InputParameters & parameters);
+
+  /**
+   * Solve the least-squares problem. Use the fitted coefficients to calculate  the value at the
+   * requested point.
+   *
+   * @param p           Point at which to compute the fitted value
+   * @param elem_ids    Ids of the elements in the patch
+   * @param component   Index of the component to compute the fitted value of
+   * @return The fitted value
+   */
+  virtual Real nodalPatchRecovery(const Point & p,
+                                  const std::vector<dof_id_type> & elem_ids,
+                                  std::size_t component) const override;
 
   virtual void initialSetup() override;
 
-protected:
-  virtual Real computeValue() override;
+  virtual void initialize() override;
+  virtual void execute() override;
+  virtual void threadJoin(const UserObject &) override;
+  virtual void finalize() override;
+
+private:
+  /**
+   * Compute the P vector at a given point
+   * i.e. given dim = 2, order = 2, polynomial P has the following terms:
+   * 1
+   * x
+   * y
+   * x^2
+   * xy
+   * y^2
+   *
+   * @param q_point point at which to evaluate the polynomial basis
+   */
+  RealEigenVector evaluateBasisFunctions(const Point & q_point) const;
+
+  /// data type stored for each element
+  typedef std::pair<RealEigenMatrix, std::vector<RealEigenVector>> ElementData;
 
+  /// current quadrature point
+  unsigned int _qp;
+
+  /// number of scalar components in the recovered type
+  std::size_t _n_components;
+
+  /// The polynomial order, default is variable order
+  const unsigned int _patch_polynomial_order;
+
+  /// The multi-index table
+  const std::vector<std::vector<unsigned int>> _multi_index;
+
+  /// Number of basis functions
+  const unsigned int _q;
+
+  /// stored property
+  const GenericMaterialProperty<T, is_ad> & _prop;
+
+  /// The element-level A matrix and the element-level b vectors for each component
+  std::map<dof_id_type, ElementData> _abs;
+
+  /// Current subdomain
   const SubdomainID & _current_subdomain_id;
 
+  /// list of require materials that need to be explicityly initialized at step zero
   std::map<SubdomainID, std::vector<MaterialBase *>> _required_materials;
 };
+
+typedef ProjectedStatefulMaterialNodalPatchRecoveryTempl<Real, false>
+    ProjectedStatefulMaterialNodalPatchRecoveryReal;
+typedef ProjectedStatefulMaterialNodalPatchRecoveryTempl<Real, true>
+    ADProjectedStatefulMaterialNodalPatchRecoveryReal;
+typedef ProjectedStatefulMaterialNodalPatchRecoveryTempl<RealVectorValue, false>
+    ProjectedStatefulMaterialNodalPatchRecoveryRealVectorValue;
+typedef ProjectedStatefulMaterialNodalPatchRecoveryTempl<RealVectorValue, true>
+    ADProjectedStatefulMaterialNodalPatchRecoveryRealVectorValue;
+typedef ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankTwoTensor, false>
+    ProjectedStatefulMaterialNodalPatchRecoveryRankTwoTensor;
+typedef ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankTwoTensor, true>
+    ADProjectedStatefulMaterialNodalPatchRecoveryRankTwoTensor;
+typedef ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankFourTensor, false>
+    ProjectedStatefulMaterialNodalPatchRecoveryRankFourTensor;
+typedef ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankFourTensor, true>
+    ADProjectedStatefulMaterialNodalPatchRecoveryRankFourTensor;
diff --git a/framework/include/utils/SerialAccess.h b/framework/include/utils/SerialAccess.h
index 3f4d81561c86..c8a7241b9ca5 100644
--- a/framework/include/utils/SerialAccess.h
+++ b/framework/include/utils/SerialAccess.h
@@ -12,6 +12,10 @@
 // MOOSE includes
 #include "Moose.h"
 #include "MooseTypes.h"
+#include "RankTwoTensorForward.h"
+#include "RankFourTensorForward.h"
+
+#include <tuple>
 
 namespace Moose
 {
@@ -23,48 +27,53 @@ namespace Moose
 template <typename T>
 struct SerialAccess
 {
-  static typename T::value_type * data(T & obj)
-  {
-    static_assert(always_false<T>,
-                  "Specialize SerialAccess for this type and implement the data() method.");
-  }
-  constexpr std::size_t * size(T &)
-  {
-    static_assert(always_false<T>,
-                  "Specialize SerialAccess for this type and implement the size() method.");
-    return 0;
-  }
+  static_assert(always_false<T>, "Specialize SerialAccess for this type.");
 };
 
-/// simple Real specialization
-template <>
-struct SerialAccess<Real>
-{
-  static Real * data(Real & obj) { return &obj; }
-  static constexpr std::size_t size(Real &) { return 1u; }
-};
-template <>
-struct SerialAccess<ADReal>
-{
-  static ADReal * data(ADReal & obj) { return &obj; }
-  static constexpr std::size_t size(ADReal &) { return 1u; }
-};
+// Specializations for scalar types
+#define SERIAL_ACCESS_SCALAR(type)                                                                 \
+  template <>                                                                                      \
+  struct SerialAccess<type>                                                                        \
+  {                                                                                                \
+    static type * data(type & obj) { return &obj; }                                                \
+    static constexpr std::size_t size(type &) { return 1u; }                                       \
+    static constexpr std::size_t size() { return 1u; }                                             \
+  }
 
-/// (AD)RealVectorValue etc. specialization
-template <typename T>
-struct SerialAccess<VectorValue<T>>
-{
-  static T * data(VectorValue<T> & obj) { return &obj(0u); }
-  static constexpr std::size_t size(VectorValue<T> &) { return Moose::dim; }
-};
+SERIAL_ACCESS_SCALAR(Real);
+SERIAL_ACCESS_SCALAR(const Real);
+SERIAL_ACCESS_SCALAR(ADReal);
+SERIAL_ACCESS_SCALAR(const ADReal);
+
+// constant size containers
+#define SERIAL_ACCESS_CONST_SIZE(type, dataptr, sizeval)                                           \
+  template <typename T>                                                                            \
+  struct SerialAccess<type<T>>                                                                     \
+  {                                                                                                \
+    static auto * data(type<T> & obj) { return dataptr; }                                          \
+    static constexpr std::size_t size(type<T> &) { return sizeval; }                               \
+    static constexpr std::size_t size() { return sizeval; }                                        \
+  }
 
-/// DenseVector specialization
-template <typename T>
-struct SerialAccess<DenseVector<T>>
-{
-  static T * data(DenseVector<T> & obj) { return &obj(0u); }
-  static std::size_t size(DenseVector<T> & obj) { return obj.size(); }
-};
+SERIAL_ACCESS_CONST_SIZE(VectorValue, &obj(0u), Moose::dim);
+SERIAL_ACCESS_CONST_SIZE(const VectorValue, &obj(0u), Moose::dim);
+SERIAL_ACCESS_CONST_SIZE(RankTwoTensorTempl, &obj(0u, 0u), RankTwoTensorTempl<T>::N2);
+SERIAL_ACCESS_CONST_SIZE(const RankTwoTensorTempl, &obj(0u, 0u), RankTwoTensorTempl<T>::N2);
+SERIAL_ACCESS_CONST_SIZE(RankFourTensorTempl, &obj(0u, 0u, 0u, 0u), RankFourTensorTempl<T>::N4);
+SERIAL_ACCESS_CONST_SIZE(const RankFourTensorTempl,
+                         &obj(0u, 0u, 0u, 0u),
+                         RankFourTensorTempl<T>::N4);
+
+// dynamic size containers (determining size requires an object instance)
+#define SERIAL_ACCESS_DYNAMIC_SIZE(type, dataptr, sizeval)                                         \
+  template <typename T>                                                                            \
+  struct SerialAccess<type<T>>                                                                     \
+  {                                                                                                \
+    static auto * data(type<T> & obj) { return dataptr; }                                          \
+    static constexpr std::size_t size(type<T> & obj) { return sizeval; }                           \
+  }
+
+SERIAL_ACCESS_DYNAMIC_SIZE(DenseVector, &obj(0u), obj.size());
 
 /**
  * Value type helper (necessary for any type that does not have a value_type
@@ -89,7 +98,8 @@ struct SerialAccessVlaueTypeHelper<Real>
 template <typename T>
 class SerialAccessRange
 {
-  typedef typename SerialAccessVlaueTypeHelper<T>::value_type V;
+  typedef typename SerialAccessVlaueTypeHelper<typename std::remove_const<T>::type>::value_type R;
+  typedef typename std::conditional<std::is_const_v<T>, const R, R>::type V;
 
 public:
   class iterator
@@ -126,9 +136,10 @@ class SerialAccessRange
   }
 
   iterator begin() const { return _begin; }
-
   iterator end() const { return _end; }
 
+  V & operator[](int i) { return *(&*_begin + i); }
+
 private:
   iterator _begin, _end;
 };
@@ -144,16 +155,27 @@ serialAccess(T & obj)
 template <typename... Ts>
 struct TypeList
 {
+  typedef std::tuple<Ts...> Tuple;
+  typedef std::tuple<Ts *...> PointerTuple;
+  static constexpr std::size_t size = sizeof...(Ts);
 };
 
 /// Type loop
-template <template <typename> class L, typename T, typename... Ts, typename... As>
+template <template <typename, int> class L, int I, typename T, typename... Ts, typename... As>
 void
-typeLoop(TypeList<T, Ts...>, As... args)
+typeLoopInternal(TypeList<T, Ts...>, As... args)
 {
-  L<T>::apply(args...);
+  L<T, I>::apply(args...);
   if constexpr (sizeof...(Ts) > 0)
-    typeLoop<L>(TypeList<Ts...>{}, args...);
+    typeLoopInternal<L, I + 1>(TypeList<Ts...>{}, args...);
+}
+
+/// Type loop
+template <template <typename, int> class L, typename... Ts, typename... As>
+void
+typeLoop(TypeList<Ts...>, As... args)
+{
+  typeLoopInternal<L, 0>(TypeList<Ts...>{}, args...);
 }
 
 } // namespace Moose;
diff --git a/framework/include/utils_nonunity/RankFourTensor.h b/framework/include/utils_nonunity/RankFourTensor.h
index ed92a84f208c..2a8733d7c7f1 100644
--- a/framework/include/utils_nonunity/RankFourTensor.h
+++ b/framework/include/utils_nonunity/RankFourTensor.h
@@ -56,18 +56,6 @@ template <>
 void mooseSetToZero<ADRankFourTensor>(ADRankFourTensor & v);
 }
 
-namespace Moose
-{
-template <typename T>
-struct SerialAccess;
-template <typename T>
-struct SerialAccess<RankFourTensorTempl<T>>
-{
-  static T * data(RankFourTensorTempl<T> & obj) { return &obj(0, 0, 0, 0); }
-  static constexpr std::size_t size(RankFourTensorTempl<T> &) { return RankFourTensorTempl<T>::N4; }
-};
-}
-
 /**
  * RankFourTensorTempl is designed to handle any N-dimensional fourth order tensor, C.
  * Since N is hard-coded to 3, RankFourTensorTempl holds 81 separate C_ijkl entries,
diff --git a/framework/include/utils_nonunity/RankTwoTensor.h b/framework/include/utils_nonunity/RankTwoTensor.h
index 048d53e9802e..2e7dd558bf75 100644
--- a/framework/include/utils_nonunity/RankTwoTensor.h
+++ b/framework/include/utils_nonunity/RankTwoTensor.h
@@ -66,21 +66,6 @@ template <>
 void mooseSetToZero<ADRankTwoTensor>(ADRankTwoTensor & v);
 }
 
-namespace Moose
-{
-template <typename T>
-struct SerialAccess;
-template <typename T>
-struct SerialAccess<RankTwoTensorTempl<T>>
-{
-  static T * data(RankTwoTensorTempl<T> & obj) { return &obj(0, 0); }
-  static constexpr std::size_t size(RankTwoTensorTempl<T> &)
-  {
-    return RankTwoTensorTempl<T>::N2;
-  }
-};
-}
-
 /**
  * RankTwoTensorTempl is designed to handle the Stress or Strain Tensor for a fully anisotropic
  * material. It is designed to allow for maximum clarity of the mathematics and ease of use.
diff --git a/framework/src/actions/ProjectedStatefulMaterialStorageAction.C b/framework/src/actions/ProjectedStatefulMaterialStorageAction.C
index 8017e7f2331b..50241f680162 100644
--- a/framework/src/actions/ProjectedStatefulMaterialStorageAction.C
+++ b/framework/src/actions/ProjectedStatefulMaterialStorageAction.C
@@ -13,7 +13,6 @@
 #include "RankTwoTensorForward.h"
 #include "AddAuxVariableAction.h"
 #include "AddVariableAction.h"
-#include "Conversion.h"
 #include "FEProblemBase.h"
 #include "Factory.h"
 
@@ -66,169 +65,22 @@ void
 ProjectedStatefulMaterialStorageAction::addRelationshipManagers(
     Moose::RelationshipManagerType input_rm_type)
 {
-  auto params = _factory.getValidParams("ProjectedStatefulMaterialNodalPatchRecovery");
+  auto params = ProjectedStatefulMaterialNodalPatchRecoveryBase::validParams();
   addRelationshipManagers(input_rm_type, params);
 }
 
 void
-ProjectedStatefulMaterialStorageAction::processComponent(const std::string & prop_name,
-                                                         std::vector<unsigned int> idx,
-                                                         std::vector<VariableName> & vars,
-                                                         bool is_ad)
+ProjectedStatefulMaterialStorageAction::act()
 {
-  const std::string pomps_prefix = "_pomps_";
-  auto name = [&](const std::string & type)
-  {
-    auto name = _pomps_prefix + type + "_" + prop_name;
-    if (idx.size())
-      name += "_" + Moose::stringify(idx, "_");
-    return name;
-  };
-
-  const auto var_name = name("var");
-  vars.push_back(var_name);
-  const auto uo_name = name("uo");
-
-  if (_current_task == "setup_projected_properties")
-  {
-    // add the AuxVars for storage
-    auto params = _factory.getValidParams(_var_type);
-    params.applyParameters(parameters());
-    params.set<std::vector<OutputName>>("outputs") = {"none"};
-    _problem->addAuxVariable(_var_type, var_name, params);
-
-    // use nodal patch recovery for lagrange
-    if (_fe_type.family == LAGRANGE)
-    {
-      // nodal variables require patch recovery (add user object)
-      const auto & type_name = is_ad ? "ADProjectedStatefulMaterialNodalPatchRecovery"
-                                     : "ProjectedStatefulMaterialNodalPatchRecovery";
-      auto params = _factory.getValidParams(type_name);
-      params.applySpecificParameters(parameters(), {"block"});
-      params.set<std::vector<unsigned int>>("component") = idx;
-      params.set<MaterialPropertyName>("property") = prop_name;
-      params.set<MooseEnum>("patch_polynomial_order") = _order;
-      params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
-      params.set<bool>("force_preaux") = true;
-      _problem->addUserObject(type_name, uo_name, params);
-    }
-  }
-
-  if (_current_task == "add_aux_kernel")
+  for (const auto & prop_name : _prop_names)
   {
-    // use nodal patch recovery for lagrange
-    if (_fe_type.family == LAGRANGE)
-    {
-      // nodal variables require patch recovery (add aux kernel)
-      auto params = _factory.getValidParams("NodalPatchRecoveryAux");
-      params.applySpecificParameters(parameters(), {"block"});
-      params.set<AuxVariableName>("variable") = var_name;
-      params.set<UserObjectName>("nodal_patch_recovery_uo") = uo_name;
-      params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
-      _problem->addAuxKernel("NodalPatchRecoveryAux", name("aux"), params);
-    }
-    else
-    {
-      // elemental variables
-      {
-        const auto & type_name =
-            is_ad ? "ADProjectedStatefulMaterialAux" : "ProjectedStatefulMaterialAux";
-        auto params = _factory.getValidParams(type_name);
-        params.applySpecificParameters(parameters(), {"block"});
-        params.set<AuxVariableName>("variable") = var_name;
-        params.set<std::vector<unsigned int>>("component") = idx;
-        params.set<MaterialPropertyName>("property") = prop_name;
-        params.set<ExecFlagEnum>("execute_on") = {EXEC_INITIAL, EXEC_TIMESTEP_END};
-        _problem->addAuxKernel(type_name, name("aux"), params);
-      }
-    }
+    // loop over all supported property types
+    Moose::typeLoop<ProcessProperty>(SupportedTypes{}, this, prop_name);
   }
 }
 
-ProjectedStatefulMaterialStorageAction::PropertyInfo
-ProjectedStatefulMaterialStorageAction::checkProperty(const std::string & prop_name)
-{
-  const auto & data = _problem->getMaterialData(Moose::BLOCK_MATERIAL_DATA);
-
-  if (data.haveProperty<Real>(prop_name))
-    return {PropertyType::REAL, false};
-  if (data.haveADProperty<Real>(prop_name))
-    return {PropertyType::REAL, true};
-
-  if (data.haveProperty<RealVectorValue>(prop_name))
-    return {PropertyType::REALVECTORVALUE, false};
-  if (data.haveADProperty<RealVectorValue>(prop_name))
-    return {PropertyType::REALVECTORVALUE, true};
-
-  if (data.haveProperty<RankTwoTensor>(prop_name))
-    return {PropertyType::RANKTWOTENSOR, false};
-  if (data.haveADProperty<RankTwoTensor>(prop_name))
-    return {PropertyType::RANKTWOTENSOR, true};
-
-  if (data.haveProperty<RankFourTensor>(prop_name))
-    return {PropertyType::RANKFOURTENSOR, false};
-  if (data.haveADProperty<RankFourTensor>(prop_name))
-    return {PropertyType::RANKFOURTENSOR, true};
-
-  paramError("projected_props",
-             "Material property type of property '",
-             prop_name,
-             "' is not supported for projection.");
-}
-
-void
-ProjectedStatefulMaterialStorageAction::addMaterial(const std::string & prop_type,
-                                                    const std::string & prop_name,
-                                                    std::vector<VariableName> & vars)
-{
-  auto params = _factory.getValidParams("InterpolatedStatefulMaterial");
-  params.applySpecificParameters(parameters(), {"block"});
-  params.set<std::vector<VariableName>>("old_state") = vars;
-  params.set<MooseEnum>("prop_type") = prop_type;
-  params.set<MaterialPropertyName>("prop_name") = prop_name;
-  _problem->addMaterial("InterpolatedStatefulMaterial", _pomps_prefix + "mat_" + prop_name, params);
-}
-
-void
-ProjectedStatefulMaterialStorageAction::act()
+MooseEnum
+ProjectedStatefulMaterialStorageAction::getTypeEnum()
 {
-  for (const auto & prop_name : _prop_names)
-  {
-    const auto prop_info = checkProperty(prop_name);
-    std::vector<VariableName> vars;
-
-    switch (prop_info.first)
-    {
-      case PropertyType::REAL:
-        processComponent(prop_name, {}, vars, prop_info.second);
-        if (_current_task == "setup_projected_properties")
-          addMaterial("REAL", prop_name, vars);
-        break;
-
-      case PropertyType::REALVECTORVALUE:
-        for (unsigned int i = 0; i < Moose::dim; ++i)
-          processComponent(prop_name, {i}, vars, prop_info.second);
-        if (_current_task == "setup_projected_properties")
-          addMaterial("REALVECTORVALUE", prop_name, vars);
-        break;
-
-      case PropertyType::RANKTWOTENSOR:
-        for (unsigned int i = 0; i < Moose::dim; ++i)
-          for (unsigned int j = 0; j < Moose::dim; ++j)
-            processComponent(prop_name, {i, j}, vars, prop_info.second);
-        if (_current_task == "setup_projected_properties")
-          addMaterial("RANKTWOTENSOR", prop_name, vars);
-        break;
-
-      case PropertyType::RANKFOURTENSOR:
-        for (unsigned int i = 0; i < Moose::dim; ++i)
-          for (unsigned int j = 0; j < Moose::dim; ++j)
-            for (unsigned int k = 0; k < Moose::dim; ++k)
-              for (unsigned int l = 0; l < Moose::dim; ++l)
-                processComponent(prop_name, {i, j, k, l}, vars, prop_info.second);
-        if (_current_task == "setup_projected_properties")
-          addMaterial("RANKFOURTENSOR", prop_name, vars);
-        break;
-    }
-  }
+  return getTypeEnum(SupportedTypes{});
 }
diff --git a/framework/src/auxkernels/NodalPatchRecoveryAux.C b/framework/src/auxkernels/NodalPatchRecoveryAux.C
index 5412d67e7cb2..c8c0eceecb94 100644
--- a/framework/src/auxkernels/NodalPatchRecoveryAux.C
+++ b/framework/src/auxkernels/NodalPatchRecoveryAux.C
@@ -15,73 +15,22 @@ registerMooseObject("MooseApp", NodalPatchRecoveryAux);
 InputParameters
 NodalPatchRecoveryAux::validParams()
 {
-  InputParameters params = AuxKernel::validParams();
-  params.addClassDescription("This Auxkernel solves a least squares problem at each node to fit a "
-                             "value from quantities defined on quadrature points.");
-  params.addRequiredParam<UserObjectName>(
-      "nodal_patch_recovery_uo",
-      "The name of the userobject that sets up the least squares problem of the nodal patch.");
+  InputParameters params = NodalPatchRecoveryAuxBase::validParams();
   return params;
 }
 
 NodalPatchRecoveryAux::NodalPatchRecoveryAux(const InputParameters & parameters)
-  : AuxKernel(parameters), _npr(getUserObject<NodalPatchRecoveryBase>("nodal_patch_recovery_uo"))
+  : NodalPatchRecoveryAuxBase(parameters),
+    _npr(getUserObject<NodalPatchRecoveryBase>("nodal_patch_recovery_uo"))
 {
   // Check user object block restriction for consistency
   if (!isBlockSubset(_npr.blockIDs()))
     paramError("nodal_patch_recovery_uo",
                "Nodal patch recovery auxiliary kernel is not defined in a subset of blocks of the "
-               "associateduser object. Revise your input file.");
+               "associated user object. Revise your input file.");
 }
-
 Real
-NodalPatchRecoveryAux::computeValue()
-{
-  if (!isNodal())
-    mooseError(name(), " only runs on nodal variables.");
-
-  // get node-to-conneted-elem map
-  const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map = _mesh.nodeToElemMap();
-  auto node_to_elem_pair = node_to_elem_map.find(_current_node->id());
-  mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing entry in node to elem map");
-
-  std::vector<dof_id_type> elem_ids;
-  blockRestrictElements(elem_ids, node_to_elem_pair->second);
-
-  // consider the case for corner node
-  if (elem_ids.size() == 1)
-  {
-    dof_id_type elem_id = elem_ids[0];
-    for (auto & n : _mesh.elemPtr(elem_id)->node_ref_range())
-    {
-      node_to_elem_pair = node_to_elem_map.find(n.id());
-      std::vector<dof_id_type> elem_ids_candidate = node_to_elem_pair->second;
-      if (elem_ids_candidate.size() > elem_ids.size())
-      {
-        std::vector<dof_id_type> elem_ids_candidate_restricted;
-        blockRestrictElements(elem_ids_candidate_restricted, elem_ids_candidate);
-
-        if (elem_ids_candidate_restricted.size() > elem_ids.size())
-          elem_ids = elem_ids_candidate_restricted;
-      }
-    }
-  }
-
-  return _npr.nodalPatchRecovery(*_current_node, elem_ids);
-}
-
-void
-NodalPatchRecoveryAux::blockRestrictElements(
-    std::vector<dof_id_type> & elem_ids,
-    const std::vector<dof_id_type> & node_to_elem_pair_elems) const
+NodalPatchRecoveryAux::nodalPatchRecovery()
 {
-  if (blockRestricted())
-    for (auto elem_id : node_to_elem_pair_elems)
-    {
-      for (const auto block_id : blockIDs())
-        if (block_id == _mesh.elemPtr(elem_id)->subdomain_id())
-          elem_ids.push_back(elem_id);
-    }
-  else
-    elem_ids = node_to_elem_pair_elems;
+  return _npr.nodalPatchRecovery(*_current_node, _elem_ids);
 }
diff --git a/framework/src/auxkernels/NodalPatchRecoveryAuxBase.C b/framework/src/auxkernels/NodalPatchRecoveryAuxBase.C
new file mode 100644
index 000000000000..621b7ea73a81
--- /dev/null
+++ b/framework/src/auxkernels/NodalPatchRecoveryAuxBase.C
@@ -0,0 +1,79 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "NodalPatchRecoveryAuxBase.h"
+
+InputParameters
+NodalPatchRecoveryAuxBase::validParams()
+{
+  InputParameters params = AuxKernel::validParams();
+  params.addClassDescription("This Auxkernel solves a least squares problem at each node to fit a "
+                             "value from quantities defined on quadrature points.");
+  params.addRequiredParam<UserObjectName>(
+      "nodal_patch_recovery_uo",
+      "The name of the userobject that sets up the least squares problem of the nodal patch.");
+  return params;
+}
+
+NodalPatchRecoveryAuxBase::NodalPatchRecoveryAuxBase(const InputParameters & parameters)
+  : AuxKernel(parameters)
+{
+  if (!isNodal())
+    mooseError(name(), " only runs on nodal variables.");
+}
+
+Real
+NodalPatchRecoveryAuxBase::computeValue()
+{
+  // get node-to-conneted-elem map
+  const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map = _mesh.nodeToElemMap();
+  auto node_to_elem_pair = node_to_elem_map.find(_current_node->id());
+  mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing entry in node to elem map");
+
+  _elem_ids.clear();
+  blockRestrictElements(_elem_ids, node_to_elem_pair->second);
+
+  // consider the case for corner node
+  if (_elem_ids.size() == 1)
+  {
+    dof_id_type elem_id = _elem_ids[0];
+    for (auto & n : _mesh.elemPtr(elem_id)->node_ref_range())
+    {
+      node_to_elem_pair = node_to_elem_map.find(n.id());
+      std::vector<dof_id_type> elem_ids_candidate = node_to_elem_pair->second;
+      if (elem_ids_candidate.size() > _elem_ids.size())
+      {
+        std::vector<dof_id_type> elem_ids_candidate_restricted;
+        blockRestrictElements(elem_ids_candidate_restricted, elem_ids_candidate);
+
+        if (elem_ids_candidate_restricted.size() > _elem_ids.size())
+          _elem_ids = elem_ids_candidate_restricted;
+      }
+    }
+  }
+
+  // get the value from a userobject (overridden in derived class)
+  return nodalPatchRecovery();
+}
+
+void
+NodalPatchRecoveryAuxBase::blockRestrictElements(
+    std::vector<dof_id_type> & elem_ids,
+    const std::vector<dof_id_type> & node_to_elem_pair_elems) const
+{
+  if (blockRestricted())
+    for (auto elem_id : node_to_elem_pair_elems)
+    {
+      for (const auto block_id : blockIDs())
+        if (block_id == _mesh.elemPtr(elem_id)->subdomain_id())
+          elem_ids.push_back(elem_id);
+    }
+  else
+    elem_ids = node_to_elem_pair_elems;
+}
diff --git a/framework/src/auxkernels/ProjectedMaterialPropertyNodalPatchRecoveryAux.C b/framework/src/auxkernels/ProjectedMaterialPropertyNodalPatchRecoveryAux.C
new file mode 100644
index 000000000000..debdd959b965
--- /dev/null
+++ b/framework/src/auxkernels/ProjectedMaterialPropertyNodalPatchRecoveryAux.C
@@ -0,0 +1,40 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "ProjectedMaterialPropertyNodalPatchRecoveryAux.h"
+#include "ProjectedStatefulMaterialNodalPatchRecovery.h"
+
+registerMooseObject("MooseApp", ProjectedMaterialPropertyNodalPatchRecoveryAux);
+
+InputParameters
+ProjectedMaterialPropertyNodalPatchRecoveryAux::validParams()
+{
+  InputParameters params = NodalPatchRecoveryAuxBase::validParams();
+  params.addRequiredParam<unsigned int>("component",
+                                        "Serial access component for the recovered type");
+  return params;
+}
+
+ProjectedMaterialPropertyNodalPatchRecoveryAux::ProjectedMaterialPropertyNodalPatchRecoveryAux(
+    const InputParameters & parameters)
+  : NodalPatchRecoveryAuxBase(parameters),
+    _npr(getUserObject<ProjectedStatefulMaterialNodalPatchRecoveryBase>("nodal_patch_recovery_uo")),
+    _component(getParam<unsigned int>("component"))
+{
+  // Check user object block restriction for consistency
+  if (!isBlockSubset(_npr.blockIDs()))
+    paramError("nodal_patch_recovery_uo",
+               "Nodal patch recovery auxiliary kernel is not defined in a subset of blocks of the "
+               "associated user object. Revise your input file.");
+}
+Real
+ProjectedMaterialPropertyNodalPatchRecoveryAux::nodalPatchRecovery()
+{
+  return _npr.nodalPatchRecovery(*_current_node, _elem_ids, _component);
+}
diff --git a/framework/src/auxkernels/ProjectedStatefulMaterialAux.C b/framework/src/auxkernels/ProjectedStatefulMaterialAux.C
index 4f7ca7070732..b7677ace1846 100644
--- a/framework/src/auxkernels/ProjectedStatefulMaterialAux.C
+++ b/framework/src/auxkernels/ProjectedStatefulMaterialAux.C
@@ -9,50 +9,67 @@
 
 #include "ProjectedStatefulMaterialAux.h"
 #include "NodalPatchRecoveryBase.h"
+#include "SerialAccess.h"
 #include "MaterialBase.h"
 #include "Assembly.h"
 
-registerMooseObject("MooseApp", ProjectedStatefulMaterialAux);
-registerMooseObject("MooseApp", ADProjectedStatefulMaterialAux);
+registerMooseObject("MooseApp", ProjectedStatefulMaterialRealAux);
+registerMooseObject("MooseApp", ADProjectedStatefulMaterialRealAux);
+registerMooseObject("MooseApp", ProjectedStatefulMaterialRealVectorValueAux);
+registerMooseObject("MooseApp", ADProjectedStatefulMaterialRealVectorValueAux);
+registerMooseObject("MooseApp", ProjectedStatefulMaterialRankTwoTensorAux);
+registerMooseObject("MooseApp", ADProjectedStatefulMaterialRankTwoTensorAux);
+registerMooseObject("MooseApp", ProjectedStatefulMaterialRankFourTensorAux);
+registerMooseObject("MooseApp", ADProjectedStatefulMaterialRankFourTensorAux);
 
-template <bool is_ad>
+template <typename T, bool is_ad>
 InputParameters
-ProjectedStatefulMaterialAuxTempl<is_ad>::validParams()
+ProjectedStatefulMaterialAuxTempl<T, is_ad>::validParams()
 {
-  auto params = IndexableProperty<AuxKernel, is_ad>::validParams();
+  auto params = AuxKernel::validParams();
   params.addClassDescription(
       "Picks a component of an indexable material property to get projected on an elemental "
       "Auxvariable. For use by ProjectedStatefulMaterialStorageAction.");
+  params.addRequiredParam<MaterialPropertyName>(
+      "prop", "material property to project onto an elemental aux variable");
+  params.addRequiredParam<unsigned int>("component", "scalar component of the property to fetch");
   return params;
 }
 
-template <bool is_ad>
-ProjectedStatefulMaterialAuxTempl<is_ad>::ProjectedStatefulMaterialAuxTempl(
+template <typename T, bool is_ad>
+ProjectedStatefulMaterialAuxTempl<T, is_ad>::ProjectedStatefulMaterialAuxTempl(
     const InputParameters & parameters)
-  : AuxKernel(parameters), _prop(this), _current_subdomain_id(_assembly.currentSubdomainID())
+  : AuxKernel(parameters),
+    _current_subdomain_id(_assembly.currentSubdomainID()),
+    _prop(getGenericMaterialProperty<T, is_ad>("prop")),
+    _component(getParam<unsigned int>("component"))
 {
 }
 
-template <bool is_ad>
+template <typename T, bool is_ad>
 void
-ProjectedStatefulMaterialAuxTempl<is_ad>::initialSetup()
+ProjectedStatefulMaterialAuxTempl<T, is_ad>::initialSetup()
 {
-  _prop.check();
-
   // get all material classes that provide properties for this object
   _required_materials = buildRequiredMaterials();
 }
 
-template <bool is_ad>
+template <typename T, bool is_ad>
 Real
-ProjectedStatefulMaterialAuxTempl<is_ad>::computeValue()
+ProjectedStatefulMaterialAuxTempl<T, is_ad>::computeValue()
 {
   if (_t_step == 0)
     for (const auto & mat : _required_materials[_current_subdomain_id])
       mat->initStatefulProperties(_qrule->size());
 
-  return MetaPhysicL::raw_value(_prop[_qp]);
+  return MetaPhysicL::raw_value(Moose::serialAccess(_prop[_qp])[_component]);
 }
 
-template class ProjectedStatefulMaterialAuxTempl<false>;
-template class ProjectedStatefulMaterialAuxTempl<true>;
+template class ProjectedStatefulMaterialAuxTempl<Real, false>;
+template class ProjectedStatefulMaterialAuxTempl<Real, true>;
+template class ProjectedStatefulMaterialAuxTempl<RealVectorValue, false>;
+template class ProjectedStatefulMaterialAuxTempl<RealVectorValue, true>;
+template class ProjectedStatefulMaterialAuxTempl<RankTwoTensor, false>;
+template class ProjectedStatefulMaterialAuxTempl<RankTwoTensor, true>;
+template class ProjectedStatefulMaterialAuxTempl<RankFourTensor, false>;
+template class ProjectedStatefulMaterialAuxTempl<RankFourTensor, true>;
diff --git a/framework/src/materials/InterpolatedStatefulMaterial.C b/framework/src/materials/InterpolatedStatefulMaterial.C
index cdc30dfc7ff9..74cb95fc2c07 100644
--- a/framework/src/materials/InterpolatedStatefulMaterial.C
+++ b/framework/src/materials/InterpolatedStatefulMaterial.C
@@ -8,107 +8,53 @@
 //* https://www.gnu.org/licenses/lgpl-2.1.html
 
 #include "InterpolatedStatefulMaterial.h"
+#include "SerialAccess.h"
 
-registerMooseObject("MooseApp", InterpolatedStatefulMaterial);
+registerMooseObject("MooseApp", InterpolatedStatefulMaterialReal);
+registerMooseObject("MooseApp", InterpolatedStatefulMaterialRealVectorValue);
+registerMooseObject("MooseApp", InterpolatedStatefulMaterialRankTwoTensor);
+registerMooseObject("MooseApp", InterpolatedStatefulMaterialRankFourTensor);
 
+template <typename T>
 InputParameters
-InterpolatedStatefulMaterial::validParams()
+InterpolatedStatefulMaterialTempl<T>::validParams()
 {
   InputParameters params = Material::validParams();
   params.addClassDescription("Access old state from projected data.");
   params.addRequiredCoupledVar("old_state", "The AuxVars for the coupled components");
-  params.addParam<MooseEnum>(
-      "prop_type", MooseEnum("REAL REALVECTORVALUE RANKTWOTENSOR RANKFOURTENSOR"), "Property type");
   params.addRequiredParam<MaterialPropertyName>("prop_name", "Name to emit");
   return params;
 }
 
-InterpolatedStatefulMaterial::InterpolatedStatefulMaterial(const InputParameters & parameters)
+template <typename T>
+InterpolatedStatefulMaterialTempl<T>::InterpolatedStatefulMaterialTempl(
+    const InputParameters & parameters)
   : Material(parameters),
     _old_state(coupledValuesOld("old_state")),
-    _older_state(coupledValuesOld("old_state")),
+    _older_state(coupledValuesOld("old_state")), // older is missing
+    _size(Moose::SerialAccess<T>::size()),
     _prop_name(getParam<MaterialPropertyName>("prop_name")),
-    _prop_type(getParam<MooseEnum>("prop_type").getEnum<PropType>()),
-    _prop_old_real(_prop_type == PropType::REAL
-                       ? &declareProperty<Real>(_prop_name + "_interpolated_old")
-                       : nullptr),
-    _prop_old_realvectorvalue(
-        _prop_type == PropType::REALVECTORVALUE
-            ? &declareProperty<RealVectorValue>(_prop_name + "_interpolated_old")
-            : nullptr),
-    _prop_old_ranktwotensor(_prop_type == PropType::RANKTWOTENSOR
-                                ? &declareProperty<RankTwoTensor>(_prop_name + "_interpolated_old")
-                                : nullptr),
-    _prop_old_rankfourtensor(
-        _prop_type == PropType::RANKFOURTENSOR
-            ? &declareProperty<RankFourTensor>(_prop_name + "_interpolated_old")
-            : nullptr),
-    _prop_older_real(_prop_type == PropType::REAL
-                         ? &declareProperty<Real>(_prop_name + "_interpolated_older")
-                         : nullptr),
-    _prop_older_realvectorvalue(
-        _prop_type == PropType::REALVECTORVALUE
-            ? &declareProperty<RealVectorValue>(_prop_name + "_interpolated_older")
-            : nullptr),
-    _prop_older_ranktwotensor(
-        _prop_type == PropType::RANKTWOTENSOR
-            ? &declareProperty<RankTwoTensor>(_prop_name + "_interpolated_older")
-            : nullptr),
-    _prop_older_rankfourtensor(
-        _prop_type == PropType::RANKFOURTENSOR
-            ? &declareProperty<RankFourTensor>(_prop_name + "_interpolated_older")
-            : nullptr)
+    _prop_old(declareProperty<T>(_prop_name + "_interpolated_old")),
+    _prop_older(declareProperty<T>(_prop_name + "_interpolated_older"))
 {
+  if (_old_state.size() != _size)
+    paramError("old_state", "Wrong number of compoet AuxVariables passed in.");
 }
 
+template <typename T>
 void
-InterpolatedStatefulMaterial::computeQpProperties()
+InterpolatedStatefulMaterialTempl<T>::computeQpProperties()
 {
-  switch (_prop_type)
-  {
-    case PropType::REAL:
-      (*_prop_old_real)[_qp] = (*_old_state[0])[_qp];
-      (*_prop_older_real)[_qp] = (*_older_state[0])[_qp];
-      return;
+  std::size_t index = 0;
+  for (auto & v : Moose::serialAccess(_prop_old[_qp]))
+    v = (*_old_state[index++])[_qp];
 
-    case PropType::REALVECTORVALUE:
-    {
-      std::size_t index = 0;
-      for (const auto i : make_range(Moose::dim))
-      {
-        (*_prop_old_realvectorvalue)[_qp](i) = (*_old_state[index])[_qp];
-        (*_prop_older_realvectorvalue)[_qp](i) = (*_older_state[index])[_qp];
-        ++index;
-      }
-      return;
-    }
-
-    case PropType::RANKTWOTENSOR:
-    {
-      std::size_t index = 0;
-      for (const auto i : make_range(Moose::dim))
-        for (const auto j : make_range(Moose::dim))
-        {
-          (*_prop_old_ranktwotensor)[_qp](i, j) = (*_old_state[index])[_qp];
-          (*_prop_older_ranktwotensor)[_qp](i, j) = (*_older_state[index])[_qp];
-          ++index;
-        }
-      return;
-    }
-
-    case PropType::RANKFOURTENSOR:
-    {
-      std::size_t index = 0;
-      for (const auto i : make_range(Moose::dim))
-        for (const auto j : make_range(Moose::dim))
-          for (const auto k : make_range(Moose::dim))
-            for (const auto l : make_range(Moose::dim))
-            {
-              (*_prop_old_rankfourtensor)[_qp](i, j, k, l) = (*_old_state[index])[_qp];
-              (*_prop_older_rankfourtensor)[_qp](i, j, k, l) = (*_older_state[index])[_qp];
-              ++index;
-            }
-      return;
-    }
-  }
+  index = 0;
+  for (auto & v : Moose::serialAccess(_prop_older[_qp]))
+    v = (*_older_state[index++])[_qp];
 }
+
+template class InterpolatedStatefulMaterialTempl<Real>;
+template class InterpolatedStatefulMaterialTempl<RealVectorValue>;
+template class InterpolatedStatefulMaterialTempl<RankTwoTensor>;
+template class InterpolatedStatefulMaterialTempl<RankFourTensor>;
diff --git a/framework/src/userobjects/ProjectedStatefulMaterialNodalPatchRecovery.C b/framework/src/userobjects/ProjectedStatefulMaterialNodalPatchRecovery.C
index 4ca7efe77bee..c206c7b4488c 100644
--- a/framework/src/userobjects/ProjectedStatefulMaterialNodalPatchRecovery.C
+++ b/framework/src/userobjects/ProjectedStatefulMaterialNodalPatchRecovery.C
@@ -8,40 +8,214 @@
 //* https://www.gnu.org/licenses/lgpl-2.1.html
 
 #include "ProjectedStatefulMaterialNodalPatchRecovery.h"
+#include "ElementUserObject.h"
 #include "MaterialBase.h"
 #include "Assembly.h"
 
-registerMooseObject("MooseApp", ProjectedStatefulMaterialNodalPatchRecovery);
+// TIMPI includes
+#include "timpi/communicator.h"
+#include "timpi/parallel_sync.h"
+#include "libmesh/parallel_eigen.h"
+
+registerMooseObject("MooseApp", ProjectedStatefulMaterialNodalPatchRecoveryReal);
+registerMooseObject("MooseApp", ADProjectedStatefulMaterialNodalPatchRecoveryReal);
+registerMooseObject("MooseApp", ProjectedStatefulMaterialNodalPatchRecoveryRealVectorValue);
+registerMooseObject("MooseApp", ADProjectedStatefulMaterialNodalPatchRecoveryRealVectorValue);
+registerMooseObject("MooseApp", ProjectedStatefulMaterialNodalPatchRecoveryRankTwoTensor);
+registerMooseObject("MooseApp", ADProjectedStatefulMaterialNodalPatchRecoveryRankTwoTensor);
+registerMooseObject("MooseApp", ProjectedStatefulMaterialNodalPatchRecoveryRankFourTensor);
+registerMooseObject("MooseApp", ADProjectedStatefulMaterialNodalPatchRecoveryRankFourTensor);
 
 InputParameters
-ProjectedStatefulMaterialNodalPatchRecovery::validParams()
+ProjectedStatefulMaterialNodalPatchRecoveryBase::validParams()
 {
-  InputParameters params = NodalPatchRecoveryMaterialProperty::validParams();
+  InputParameters params = ElementUserObject::validParams();
+
+  params.addRelationshipManager("ElementSideNeighborLayers",
+                                Moose::RelationshipManagerType::ALGEBRAIC,
+                                [](const InputParameters &, InputParameters & rm_params)
+                                { rm_params.set<unsigned short>("layers") = 2; });
+
+  return params;
+}
+
+template <typename T, bool is_ad>
+InputParameters
+ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::validParams()
+{
+  InputParameters params = ProjectedStatefulMaterialNodalPatchRecoveryBase::validParams();
+  MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH");
+  params.addRequiredParam<MaterialPropertyName>(
+      "property", "Name of the material property to perform nodal patch recovery on");
+  params.addRequiredParam<MooseEnum>(
+      "patch_polynomial_order",
+      orders,
+      "Polynomial order used in least squares fitting of material property "
+      "over the local patch of elements connected to a given node");
+
+  params.addRelationshipManager("ElementSideNeighborLayers",
+                                Moose::RelationshipManagerType::ALGEBRAIC,
+                                [](const InputParameters &, InputParameters & rm_params)
+                                { rm_params.set<unsigned short>("layers") = 2; });
+
+  params.addParamNamesToGroup("patch_polynomial_order", "Advanced");
+
   return params;
 }
 
-ProjectedStatefulMaterialNodalPatchRecovery::ProjectedStatefulMaterialNodalPatchRecovery(
-    const InputParameters & parameters)
-  : NodalPatchRecoveryMaterialProperty(parameters),
+template <typename T, bool is_ad>
+ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::
+    ProjectedStatefulMaterialNodalPatchRecoveryTempl(const InputParameters & parameters)
+  : ProjectedStatefulMaterialNodalPatchRecoveryBase(parameters),
+    _qp(0),
+    _n_components(Moose::SerialAccess<T>::size()),
+    _patch_polynomial_order(
+        static_cast<unsigned int>(getParam<MooseEnum>("patch_polynomial_order"))),
+    _multi_index(MathUtils::multiIndex(_mesh.dimension(), _patch_polynomial_order)),
+    _q(_multi_index.size()),
+    _prop(getGenericMaterialProperty<T, is_ad>("property")),
     _current_subdomain_id(_assembly.currentSubdomainID())
 {
 }
 
+template <typename T, bool is_ad>
 void
-ProjectedStatefulMaterialNodalPatchRecovery::initialSetup()
+ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::initialSetup()
 {
-  NodalPatchRecoveryMaterialProperty::initialSetup();
-
   // get all material classes that provide properties for this object
   _required_materials = buildRequiredMaterials();
 }
 
+template <typename T, bool is_ad>
 Real
-ProjectedStatefulMaterialNodalPatchRecovery::computeValue()
+ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::nodalPatchRecovery(
+    const Point & x, const std::vector<dof_id_type> & elem_ids, std::size_t component) const
+{
+  // Before we go, check if we have enough sample points for solving the least square fitting
+  if (_q_point.size() * elem_ids.size() < _q)
+    mooseError("There are not enough sample points to recover the nodal value, try reducing the "
+               "polynomial order or using a higher-order quadrature scheme.");
+
+  // Assemble the least squares problem over the patch
+  RealEigenMatrix A = RealEigenMatrix::Zero(_q, _q);
+  RealEigenVector b = RealEigenVector::Zero(_q);
+  for (auto elem_id : elem_ids)
+  {
+    const auto abs = libmesh_map_find(_abs, elem_id);
+    A += abs.first;
+    b += abs.second[component];
+  }
+
+  // Solve the least squares fitting
+  RealEigenVector coef = A.completeOrthogonalDecomposition().solve(b);
+
+  // Compute the fitted nodal value
+  RealEigenVector p = evaluateBasisFunctions(x);
+  return p.dot(coef);
+}
+
+template <typename T, bool is_ad>
+RealEigenVector
+ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::evaluateBasisFunctions(
+    const Point & q_point) const
+{
+  RealEigenVector p(_q);
+  Real polynomial;
+  for (unsigned int r = 0; r < _multi_index.size(); r++)
+  {
+    polynomial = 1.0;
+    mooseAssert(_multi_index[r].size() == _mesh.dimension(), "Wrong multi-index size.");
+    for (unsigned int c = 0; c < _multi_index[r].size(); c++)
+      for (unsigned int p = 0; p < _multi_index[r][c]; p++)
+        polynomial *= q_point(c);
+    p(r) = polynomial;
+  }
+  return p;
+}
+
+template <typename T, bool is_ad>
+void
+ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::initialize()
+{
+  _abs.clear();
+}
+
+template <typename T, bool is_ad>
+void
+ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::execute()
 {
   if (_t_step == 0)
     for (const auto & mat : _required_materials[_current_subdomain_id])
       mat->initStatefulProperties(_qrule->size());
 
-  return _prop[_qp];
+  std::vector<RealEigenVector> bs(_n_components, RealEigenVector::Zero(_q));
+  RealEigenMatrix Ae = RealEigenMatrix::Zero(_q, _q);
+
+  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
+  {
+    RealEigenVector p = evaluateBasisFunctions(_q_point[_qp]);
+    Ae += p * p.transpose();
+
+    std::size_t index = 0;
+    for (const auto & v : Moose::serialAccess(_prop[_qp]))
+      bs[index++] += MetaPhysicL::raw_value(v) * p;
+  }
+
+  dof_id_type elem_id = _current_elem->id();
+  _abs[elem_id] = {Ae, bs};
 }
+
+template <typename T, bool is_ad>
+void
+ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::threadJoin(const UserObject & uo)
+{
+  const auto & npr =
+      static_cast<const ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad> &>(uo);
+  _abs.insert(npr._abs.begin(), npr._abs.end());
+}
+
+template <typename T, bool is_ad>
+void
+ProjectedStatefulMaterialNodalPatchRecoveryTempl<T, is_ad>::finalize()
+{
+  // When calling nodalPatchRecovery, we may need to know _Ae and _be on algebraically ghosted
+  // elements. However, this userobject is only run on local elements, so we need to query those
+  // information from other processors in this finalize() method.
+
+  // Populate algebraically ghosted elements to query
+  std::unordered_map<processor_id_type, std::vector<dof_id_type>> query_ids;
+  const ConstElemRange evaluable_elem_range = _fe_problem.getEvaluableElementRange();
+  for (auto elem : evaluable_elem_range)
+    if (elem->processor_id() != processor_id())
+      query_ids[elem->processor_id()].push_back(elem->id());
+
+  // Answer queries received from other processors
+  auto gather_data = [this](const processor_id_type /*pid*/,
+                            const std::vector<dof_id_type> & elem_ids,
+                            std::vector<ElementData> & abs_data)
+  {
+    for (const auto i : index_range(elem_ids))
+      abs_data.emplace_back(libmesh_map_find(_abs, elem_ids[i]));
+  };
+
+  // Gather answers received from other processors
+  auto act_on_data = [this](const processor_id_type /*pid*/,
+                            const std::vector<dof_id_type> & elem_ids,
+                            const std::vector<ElementData> & abs_data)
+  {
+    for (const auto i : index_range(elem_ids))
+      _abs[elem_ids[i]] = abs_data[i];
+  };
+
+  libMesh::Parallel::pull_parallel_vector_data<ElementData>(
+      _communicator, query_ids, gather_data, act_on_data, 0);
+}
+
+template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<Real, false>;
+template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<Real, true>;
+template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RealVectorValue, false>;
+template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RealVectorValue, true>;
+template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankTwoTensor, false>;
+template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankTwoTensor, true>;
+template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankFourTensor, false>;
+template class ProjectedStatefulMaterialNodalPatchRecoveryTempl<RankFourTensor, true>;
diff --git a/test/tests/materials/projected_stateful_materials/test.i b/test/tests/materials/projected_stateful_materials/test.i
index c548aa7f04ad..1dadfd00ff92 100644
--- a/test/tests/materials/projected_stateful_materials/test.i
+++ b/test/tests/materials/projected_stateful_materials/test.i
@@ -45,6 +45,10 @@
   num_steps = 4
 []
 
+[Debug]
+  show_material_props = true
+[]
+
 [Outputs]
   csv = true
   # in initial the freshly set up old state of the material properties differes
diff --git a/unit/src/MooseServerTest.C b/unit/src/MooseServerTest.C
index e30f8ee61ecf..a3e5ae557e24 100644
--- a/unit/src/MooseServerTest.C
+++ b/unit/src/MooseServerTest.C
@@ -933,7 +933,7 @@ TEST_F(MooseServerTest, CompletionDocumentRootLevel)
 
   EXPECT_EQ(request_id, response_id);
 
-  EXPECT_EQ(47u, completions_array.size());
+  EXPECT_EQ(48u, completions_array.size());
 
   std::ostringstream completions_actual;
 
@@ -972,6 +972,7 @@ label: InterfaceKernels                 text: [InterfaceKernels]\n  \n[]
 label: Kernels                          text: [Kernels]\n  \n[]                          desc:            application named block                 pos: [42.0]-[42.0]
 label: Materials                        text: [Materials]\n  \n[]                        desc:            application named block                 pos: [42.0]-[42.0]
 label: Mesh                             text: [Mesh]\n  \n[]                             desc:            application named block                 pos: [42.0]-[42.0]
+label: MeshDivisions                    text: [MeshDivisions]\n  \n[]                    desc:            application named block                 pos: [42.0]-[42.0]
 label: MultiApps                        text: [MultiApps]\n  \n[]                        desc:            application named block                 pos: [42.0]-[42.0]
 label: NodalKernels                     text: [NodalKernels]\n  \n[]                     desc:            application named block                 pos: [42.0]-[42.0]
 label: NodalNormals                     text: [NodalNormals]\n  \n[]                     desc:            application named block                 pos: [42.0]-[42.0]
diff --git a/unit/src/SerialAccessTest.C b/unit/src/SerialAccessTest.C
index 9a5e6554c534..b4f5ae4efad8 100644
--- a/unit/src/SerialAccessTest.C
+++ b/unit/src/SerialAccessTest.C
@@ -27,6 +27,13 @@ TEST(SerialAccess, Real)
 
   EXPECT_EQ(count, 1);
   EXPECT_EQ(r, 9.0);
+
+  // test constant type
+  const auto cr = r;
+  for (auto & a : Moose::serialAccess(cr))
+    r = a + 11;
+
+  EXPECT_EQ(r, 20.0);
 }
 
 TEST(SerialAccess, ADReal)
@@ -48,12 +55,18 @@ TEST(SerialAccess, RealVectorValue)
   RealVectorValue r;
   std::size_t count = 0;
   for (auto & a : Moose::serialAccess(r))
-  {
     a = ++count;
-  }
 
   EXPECT_EQ(count, Moose::dim);
   EXPECT_EQ(r, RealVectorValue(1, 2, 3));
+
+  // test constant type
+  const auto & cr = r;
+  Real sum = 0.0;
+  for (auto & a : Moose::serialAccess(cr))
+    sum += a;
+
+  EXPECT_EQ(sum, 6);
 }
 
 TEST(SerialAccess, ADRealVectorValue)
@@ -61,9 +74,7 @@ TEST(SerialAccess, ADRealVectorValue)
   ADRealVectorValue r;
   std::size_t count = 0;
   for (auto & a : Moose::serialAccess(r))
-  {
     a = ++count;
-  }
 
   EXPECT_EQ(count, Moose::dim);
   EXPECT_EQ(r, RealVectorValue(1, 2, 3));
@@ -74,11 +85,9 @@ TEST(SerialAccess, RankTwoTensor)
   RankTwoTensor r;
   std::size_t count = 0;
   for (auto & a : Moose::serialAccess(r))
-  {
     a = ++count;
-  }
 
-  EXPECT_EQ(count, Moose::dim * Moose::dim);
+  EXPECT_EQ(count, RankTwoTensor::N2);
   EXPECT_EQ(r, RankTwoTensor(1, 4, 7, 2, 5, 8, 3, 6, 9));
 }
 
@@ -87,10 +96,9 @@ TEST(SerialAccess, RankFourTensor)
   RankFourTensor r;
   std::size_t count = 0;
   for (auto & a : Moose::serialAccess(r))
-  {
     a = ++count;
-  }
-  EXPECT_EQ(count, Moose::dim * Moose::dim * Moose::dim * Moose::dim);
+
+  EXPECT_EQ(count, RankFourTensor::N4);
 
   count = 0;
   RankFourTensor g;
@@ -101,4 +109,12 @@ TEST(SerialAccess, RankFourTensor)
           g(i, j, k, l) = ++count;
 
   EXPECT_EQ((r - g).L2norm(), 0.0);
+
+  // test constant AD type
+  const ADRankFourTensor cr = r;
+  ADReal sum = 0.0;
+  for (auto & a : Moose::serialAccess(cr))
+    sum += a;
+
+  EXPECT_EQ(sum, 3321);
 }