Add FEType::p_refinement

include/base/dof_map.h

@@ -1764,8 +1764,9 @@ class DofMap : public DofMapBase,
                                                          bool use_condensed_system = false) const;
 
   /**
-   * Describe whether the given variable group should be p-refined. If this API is not called with
-   * \p false, the default is to p-refine
+   * Set whether the given variable group should be p-refined on a
+   * p-refined Elem.  This changes the FEType of the variable group to
+   * enable or disable p-refinement.
    */
   void should_p_refine(unsigned int g, bool p_refine);
 
@@ -2257,11 +2258,6 @@ class DofMap : public DofMapBase,
    * non-SCALAR variable v
    */
   std::vector<dof_id_type> _first_old_scalar_df;
-
-  /**
-   * A container of variable groups that we should not p-refine
-   */
-  std::unordered_set<unsigned int> _dont_p_refine;
 #endif
 
 #ifdef LIBMESH_ENABLE_CONSTRAINTS
@@ -2568,14 +2564,9 @@ inline
 void DofMap::should_p_refine(const unsigned int g, const bool p_refine)
 {
 #ifdef LIBMESH_ENABLE_AMR
-  if (p_refine)
-  {
-    auto it = _dont_p_refine.find(g);
-    if (it != _dont_p_refine.end())
-      _dont_p_refine.erase(it);
-  }
-  else
-    _dont_p_refine.insert(g);
+  VariableGroup & var = _variable_groups[g];
+  var.type().p_refinement = p_refine;
+
 #else
   libmesh_ignore(g, p_refine);
 #endif
@@ -2585,7 +2576,8 @@ inline
 bool DofMap::should_p_refine(const unsigned int g) const
 {
 #ifdef LIBMESH_ENABLE_AMR
-  return !_dont_p_refine.count(g);
+  const VariableGroup & var = this->variable_group(g);
+  return var.type().p_refinement;
 #else
   libmesh_ignore(g);
   return false;
@@ -2604,7 +2596,7 @@ bool DofMap::should_p_refine_var(const unsigned int var) const
 {
 #ifdef LIBMESH_ENABLE_AMR
   const auto vg = this->var_group_from_var_number(var);
-  return !_dont_p_refine.count(vg);
+  return this->should_p_refine(vg);
 #else
   libmesh_ignore(var);
   return false;
@@ -2640,12 +2632,7 @@ void DofMap::_dof_indices (const Elem & elem,
 
       FEType fe_type = var.type();
 
-      const bool add_p_level =
-#ifdef LIBMESH_ENABLE_AMR
-          !_dont_p_refine.count(vg);
-#else
-          false;
-#endif
+      const bool add_p_level = fe_type.p_refinement;
 
 #ifdef DEBUG
       // The number of dofs per element is non-static for subdivision FE

include/base/variable.h

@@ -144,6 +144,13 @@ class Variable
   const FEType & type() const
   { return _type; }
 
+  /**
+   * \returns The \p FEType for this variable.  Altering this while
+   * this Variable is already in use may be dangerous!
+   */
+  FEType & type()
+  { return _type; }
+
   /**
    * \returns The number of components of this variable if the \p FEFamily is \p SCALAR or if the
    * associated \p FEFieldType is \p TYPE_SCALAR. Otherwise this will error because determination of

include/fe/fe_type.h

@@ -200,17 +200,24 @@ class FEType
 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
   /**
-   * Constructor.  Optionally takes the approximation \p Order
-   * and the finite element family \p FEFamily
+   * Constructor.  Optionally takes the approximation \p Order,
+   * the finite element family \p FEFamily.
+   *
+   * We can't set p_refinement in this argument list without
+   * conflicting with the \p ro parameter in the InfFE-compatible
+   * constructor version below, so we use the with_p_refinement API
+   * below to potentially convert an FEType with p-refinement (the
+   * default) to one without.
    */
   FEType(const int      o = 1,
          const FEFamily f = LAGRANGE) :
     order(o),
-    family(f)
+    family(f),
+    p_refinement(true)
   {}
 
   /**
-   * The approximation order of the element.
+   * The approximation order of the element (at 0 p-refinement level).
    */
   OrderWrapper order;
 
@@ -220,6 +227,12 @@ class FEType
    */
   FEFamily family;
 
+  /**
+   * Whether or not the finite elements for this type increase their p
+   * refinement level on geometric elements of increased p level.
+   */
+  bool p_refinement;
+
 #else
 
   /**
@@ -230,17 +243,24 @@ class FEType
    * are the same, as with the \p family and \p base_family.  It must
    * be so, otherwise what we switch on would change when infinite
    * elements are not compiled in.
+   *
+   * We can't set p_refinement in this argument list in a way
+   * that matches the non-InfFE-enabled constructor version above, so
+   * we use the with_p_refinement API below to potentially convert an
+   * FEType with p-refinement (the default) to one without.
    */
   FEType(const int        o  = 1,
          const FEFamily   f  = LAGRANGE,
          const int        ro = THIRD,
          const FEFamily   rf = JACOBI_20_00,
-         const InfMapType im = CARTESIAN) :
+         const InfMapType im = CARTESIAN
+         ) :
     order(o),
     radial_order(ro),
     family(f),
     radial_family(rf),
-    inf_map(im)
+    inf_map(im),
+    p_refinement(true)
   {}
 
   /**
@@ -274,15 +294,33 @@ class FEType
    */
   InfMapType inf_map;
 
+  /**
+   * Whether or not the finite elements for this type increase their p
+   * refinement level on geometric elements of increased p level.
+   */
+  bool p_refinement;
+
 #endif // ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
+  /**
+   * "Fluent API" for constructing a non-default p_refinement, for
+   * easier compatibility between non-InfFE and InfFE code
+   */
+  FEType with_p_refinement(bool p)
+  {
+    FEType returnval = *this;
+    returnval.p_refinement = p;
+    return returnval;
+  }
+
   /**
    * Tests equality
    */
   bool operator== (const FEType & f2) const
   {
     return (order == f2.order
             && family == f2.family
+            && p_refinement == f2.p_refinement
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
             && radial_order == f2.radial_order
             && radial_family == f2.radial_family
@@ -308,7 +346,8 @@ class FEType
       return (order < f2.order);
     if (family != f2.family)
       return (family < f2.family);
-
+    if (p_refinement != f2.p_refinement)
+      return (p_refinement < f2.p_refinement);
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
     if (radial_order != f2.radial_order)
       return (radial_order < f2.radial_order);
@@ -321,16 +360,19 @@ class FEType
   }
 
   /**
-   * \returns The default quadrature order for this \p FEType.  The
-   * default quadrature order is calculated assuming a polynomial of
-   * degree \p order and is based on integrating the mass matrix for
-   * such an element exactly on affine elements.
+   * \returns The default quadrature order for this \p FEType, on an
+   * element without p refinement.  The default quadrature order is
+   * calculated assuming a polynomial of degree \p order and is based
+   * on integrating the mass matrix for such an element exactly on
+   * affine elements.
    */
   Order default_quadrature_order () const;
 
   /**
    * \returns A quadrature rule of appropriate type and order for this \p
-   * FEType.  The default quadrature rule is based on integrating the mass
+   * FEType, on an element without p refinement.
+   *
+   * The default quadrature rule is based on integrating the mass
    * matrix for such an element exactly, with an additional power on
    * the basis order to help account for nonlinearities and/or
    * nonuniform coefficients.  Higher or lower degree rules can be
@@ -341,7 +383,9 @@ class FEType
 
   /**
    * \returns The default quadrature order for integrating unweighted
-   * basis functions of this \p FEType.
+   * basis functions of this \p FEType, on an element without p
+   * refinement.
+   *
    * The unweighted quadrature order is calculated assuming a
    * polynomial of degree \p order and is based on integrating the
    * shape functions for such an element exactly on affine elements.
@@ -350,10 +394,12 @@ class FEType
 
   /**
    * \returns A quadrature rule of appropriate type and order for
-   * unweighted integration of this \p FEType.  The default quadrature
-   * rule is based on integrating the shape functions on an affine
-   * element exactly.  Higher or lower degree rules can be chosen by
-   * changing the extraorder parameter.
+   * unweighted integration of this \p FEType, on an element without p
+   * refinement
+   *
+   * The default quadrature rule is based on integrating the shape
+   * functions on an affine element exactly.  Higher or lower degree
+   * rules can be chosen by changing the extraorder parameter.
    */
   std::unique_ptr<QBase> unweighted_quadrature_rule (const unsigned int dim,
                                                      const int extraorder=0) const;
@@ -392,7 +438,13 @@ struct hash<libMesh::FEType>
       // Old compiler versions seem to need the static_cast
       libMesh::boostcopy::hash_combine(seed, static_cast<int>(fe_type.family));
       libMesh::boostcopy::hash_combine(seed, fe_type.order._order);
-      return seed;
+      libMesh::boostcopy::hash_combine(seed, static_cast<int>(fe_type.p_refinement));
+#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
+      libMesh::boostcopy::hash_combine(seed, static_cast<int>(fe_type.radial_family));
+      libMesh::boostcopy::hash_combine(seed, fe_type.radial_order._order);
+      libMesh::boostcopy::hash_combine(seed, static_cast<int>(fe_type.inf_map));
+#endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
+     return seed;
     }
 };
 }

include/systems/generic_projector.h

@@ -1628,9 +1628,7 @@ void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy
           if (!var.active_on_subdomain(elem->subdomain_id()))
             continue;
           const FEType fe_type = var.type();
-          const auto & dof_map = this->projector.system.get_dof_map();
-          const auto vg = dof_map.var_group_from_var_number(v_num);
-          const bool add_p_level = dof_map.should_p_refine(vg);
+          const bool add_p_level = fe_type.p_refinement;
 
           // If we're trying to do projections on an isogeometric
           // analysis mesh, only the finite element nodes constrained
@@ -2305,8 +2303,7 @@ void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectVert
                           base_fe_type,
                           &elem,
                           elem.get_node_index(&vertex),
-                          this->projector.system.get_dof_map().should_p_refine(
-                              this->projector.system.get_dof_map().var_group_from_var_number(var))),
+                          base_fe_type.p_refinement),
                       (unsigned int)(1 + dim));
 
                   const FValue val =
@@ -2561,8 +2558,7 @@ void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectSide
             continue;
 
           FEType fe_type = base_fe_type;
-          const auto & dof_map = system.get_dof_map();
-          const bool add_p_level = dof_map.should_p_refine_var(var);
+          const bool add_p_level = base_fe_type.p_refinement;
 
           // This may be a p refined element
           fe_type.order = fe_type.order + add_p_level*elem.p_level();
@@ -2721,9 +2717,7 @@ void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectInte
           context.get_element_fe( var, fe, dim );
 
           FEType fe_type = base_fe_type;
-          const auto & dof_map = system.get_dof_map();
-          const auto vg = dof_map.var_group_from_var_number(var);
-          const bool add_p_level = dof_map.should_p_refine(vg);
+          const bool add_p_level = fe_type.p_refinement;
 
           // This may be a p refined element
           fe_type.order = fe_type.order + add_p_level * elem->p_level();

include/systems/system.h

@@ -1189,12 +1189,15 @@ class System : public ReferenceCountedObject<System>,
    * for this system.  Same as before, but assumes \p LAGRANGE
    * as default value for \p FEType.family. If \p active_subdomains is either
    * \p nullptr (the default) or points to an empty set, then it will be assumed
-   * that \p var has no subdomain restrictions
+   * that \p var has no subdomain restrictions. If \p p_refinement is
+   * false, then even when on an \p Elem with non-zero \p p_level()
+   * this variable will not be p-refined.
    */
   unsigned int add_variable (std::string_view var,
                              const Order order = FIRST,
                              const FEFamily = LAGRANGE,
-                             const std::set<subdomain_id_type> * const active_subdomains = nullptr);
+                             const std::set<subdomain_id_type> * const active_subdomains = nullptr,
+                             const bool p_refinement = true);
 
   /**
    * Adds the variables \p vars to the list of variables
@@ -1235,12 +1238,15 @@ class System : public ReferenceCountedObject<System>,
    * for this system.  Same as before, but assumes \p LAGRANGE
    * as default value for \p FEType.family. If \p active_subdomains is either
    * \p nullptr (the default) or points to an empty set, then it will be assumed that
-   * \p var has no subdomain restrictions
+   * \p var has no subdomain restrictions. If \p p_refinement is
+   * false, then even when on an \p Elem with non-zero \p p_level()
+   * this variable will not be p-refined.
    */
   unsigned int add_variables (const std::vector<std::string> & vars,
                               const Order order = FIRST,
                               const FEFamily = LAGRANGE,
-                              const std::set<subdomain_id_type> * const active_subdomains = nullptr);
+                              const std::set<subdomain_id_type> * const active_subdomains = nullptr,
+                              const bool p_refinement = true);
 
   /**
    * Return a constant reference to \p Variable \p var.

src/base/dof_map.C

@@ -579,12 +579,7 @@ void DofMap::reinit
       const unsigned int n_var_in_group = vg_description.n_variables();
       const FEType & base_fe_type        = vg_description.type();
 
-      const bool add_p_level =
-#ifdef LIBMESH_ENABLE_AMR
-          !_dont_p_refine.count(vg);
-#else
-          false;
-#endif
+      const bool add_p_level = base_fe_type.p_refinement;
 
       // Don't need to loop over elements for a SCALAR variable
       // Just increment _n_SCALAR_dofs
@@ -2510,12 +2505,7 @@ void DofMap::_node_dof_indices (const Elem & elem,
   const bool extra_hanging_dofs =
     FEInterface::extra_hanging_dofs(fe_type);
 
-  const bool add_p_level =
-#ifdef LIBMESH_ENABLE_AMR
-      !_dont_p_refine.count(vg);
-#else
-      false;
-#endif
+  const bool add_p_level = fe_type.p_refinement;
 
   // There is a potential problem with h refinement.  Imagine a
   // quad9 that has a linear FE on it.  Then, on the hanging side,
@@ -2765,12 +2755,8 @@ void DofMap::old_dof_indices (const Elem * const elem,
                     // or just for the specified variable
 
                     FEType fe_type = var.type();
-                    const bool add_p_level =
-#ifdef LIBMESH_ENABLE_AMR
-                        !_dont_p_refine.count(vg);
-#else
-                        false;
-#endif
+                    const bool add_p_level = fe_type.p_refinement;
+
                     // Increase the polynomial order on p refined elements,
                     // but make sure you get the right polynomial order for
                     // the OLD degrees of freedom

src/fe/fe_base.C

@@ -1550,7 +1550,7 @@ FEGenericBase<OutputType>::compute_proj_constraints (DofConstraints & constraint
 
   const Variable & var = dof_map.variable(variable_number);
   const FEType & base_fe_type = var.type();
-  const bool add_p_level = dof_map.should_p_refine_var(variable_number);
+  const bool add_p_level = base_fe_type.p_refinement;
 
   // Construct FE objects for this element and its neighbors.
   std::unique_ptr<FEGenericBase<OutputShape>> my_fe