Merge pull request #4398 from roystgnr/fe_type_p_refinement

Add FEType::p_refinement

Author: roystgnr

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,15 @@ 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;
+
+  for (auto v : make_range(var.first_scalar_number(0),
+                           var.first_scalar_number(0) +
+                           var.n_variables()))
+    this->_variables[v].type().p_refinement = p_refine;
+
+
 #else
   libmesh_ignore(g, p_refine);
 #endif
@@ -2585,7 +2582,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 +2602,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 +2638,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

@@ -188,6 +188,18 @@ class Variable
   unsigned int            _number;
   unsigned int            _first_scalar_number;
   FEType                  _type;
+
+private:
+  /**
+   * \returns The \p FEType for this variable.  Altering this while
+   * this Variable is already in use may be dangerous!
+   */
+  FEType & type()
+  { return _type; }
+
+  // DofMap can change a VariableGroup type() to disable/enable
+  // p-refinement post-variable-addition.
+  friend class DofMap;
 };
 
 

include/fe/fe_type.h

@@ -197,23 +197,30 @@ class FEType
 {
 public:
 
+  /**
+   * The approximation order of the element (at 0 p-refinement level).
+   */
+  OrderWrapper order;
+
 #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 set_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.
-   */
-  OrderWrapper order;
-
   /**
    * The type of finite element.  Valid types are \p LAGRANGE,
    * \p HIERARCHIC, etc...
@@ -230,24 +237,26 @@ 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 set_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)
   {}
 
-  /**
-   * The approximation order in the base of the infinite element.
-   */
-  OrderWrapper order;
-
   /**
    * The approximation order in radial direction of the infinite element.
    */
@@ -276,13 +285,40 @@ class FEType
 
 #endif // ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
+  /**
+   * Whether or not the finite elements for this type increase their p
+   * refinement level on geometric elements of increased p level.
+   */
+  bool p_refinement;
+
+  /**
+   * "Fluent API" for constructing a non-default p_refinement, for
+   * easier compatibility between non-InfFE and InfFE code
+   */
+  FEType set_p_refinement(bool p) &
+  {
+    this->p_refinement = p;
+    return *this;
+  }
+
+  /**
+   * Specialization for rvalues so we don't return a dangling
+   * reference to a temporary.
+   */
+  FEType set_p_refinement(bool p) &&
+  {
+    this->p_refinement = p;
+    return std::move(*this);
+  }
+
   /**
    * 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 +344,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 +358,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 +381,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 +392,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 +436,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.