Merge pull request #3652 from lindsayad/more-add-p-level

More add_p_level API additions

Author: lindsayad

include/base/dof_map.h

@@ -1621,6 +1621,11 @@ class DofMap : public ReferenceCountedObject<DofMap>,
    */
   bool should_p_refine(unsigned int g) const;
 
+  /**
+   * Whether the given variable should be p-refined
+   */
+  bool should_p_refine_var(unsigned int var) const;
+
   // Prevent bad user implicit conversions
   void should_p_refine(FEFamily, bool) = delete;
   void should_p_refine(Order, bool) = delete;
@@ -2344,6 +2349,18 @@ unsigned int DofMap::var_group_from_var_number(const unsigned int var_num) const
   libmesh_assert(var_num < n_variables());
   return libmesh_map_find(_var_to_vg, var_num);
 }
+
+inline
+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);
+#else
+  libmesh_ignore(var);
+  return false;
+#endif
+}
 } // namespace libMesh
 
 #endif // LIBMESH_DOF_MAP_H

include/fe/fe.h

@@ -364,7 +364,8 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
    */
   static void nodal_soln(const Elem * elem, const Order o,
                          const std::vector<Number> & elem_soln,
-                         std::vector<Number> & nodal_soln);
+                         std::vector<Number> & nodal_soln,
+                         bool add_p_level = true);
 
   /**
    * Build the nodal soln on one side from the (full) element soln.
@@ -375,7 +376,8 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
   static void side_nodal_soln(const Elem * elem, const Order o,
                               const unsigned int side,
                               const std::vector<Number> & elem_soln,
-                              std::vector<Number> & nodal_soln_on_side);
+                              std::vector<Number> & nodal_soln_on_side,
+                              bool add_p_level = true);
 
   /**
    * \returns The number of shape functions associated with
@@ -716,7 +718,8 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
   static void default_side_nodal_soln(const Elem * elem, const Order o,
                                       const unsigned int side,
                                       const std::vector<Number> & elem_soln,
-                                      std::vector<Number> & nodal_soln_on_side);
+                                      std::vector<Number> & nodal_soln_on_side,
+                                      bool add_p_level = true);
 
   /**
    * An array of the node locations on the last
@@ -1401,7 +1404,8 @@ OutputShape fe_fdm_second_deriv(const ElemType type,
 void lagrange_nodal_soln(const Elem * elem,
                          const Order order,
                          const std::vector<Number> & elem_soln,
-                         std::vector<Number> &       nodal_soln);
+                         std::vector<Number> &       nodal_soln,
+                         bool add_p_level = true);
 
 /**
  * Helper functions for Discontinuous-Pn type basis functions.

include/fe/fe_abstract.h

@@ -609,6 +609,11 @@ class FEAbstract : public ReferenceCountedObject<FEAbstract>
    */
   void add_p_level_in_reinit(bool value) { _add_p_level_in_reinit = value; }
 
+  /**
+   * Whether to add p-refinement levels in init/reinit methods
+   */
+  bool add_p_level_in_reinit() const { return _add_p_level_in_reinit; }
+
 protected:
 
   /**

include/fe/fe_interface.h

@@ -95,7 +95,8 @@ class FEInterface
    * Non-deprecated version of function above.
    */
   static unsigned int n_shape_functions(const FEType & fe_t,
-                                        const Elem * elem);
+                                        const Elem * elem,
+                                        const bool add_p_level = true);
 
   /**
    * Same as above, but ignores the elem->p_level() and uses the
@@ -226,7 +227,8 @@ class FEInterface
    * The non-deprecated version of the function above.
    */
   static unsigned int n_dofs_per_elem(const FEType & fe_t,
-                                      const Elem * elem);
+                                      const Elem * elem,
+                                      bool add_p_level = true);
 
   /**
    * Same thing but internally elem->p_level() is ignored and extra_order is used instead.
@@ -289,7 +291,8 @@ class FEInterface
                          const FEType & fe_t,
                          const Elem * elem,
                          const std::vector<Number> & elem_soln,
-                         std::vector<Number> & nodal_soln);
+                         std::vector<Number> & nodal_soln,
+                         bool add_p_level = true);
 
 
   /**
@@ -304,7 +307,8 @@ class FEInterface
                               const Elem * elem,
                               const unsigned int side,
                               const std::vector<Number> & elem_soln,
-                              std::vector<Number> & nodal_soln);
+                              std::vector<Number> & nodal_soln,
+                              bool add_p_level = true);
 
   /**
    * This is now deprecated; use FEMap::map instead.
@@ -386,12 +390,13 @@ class FEInterface
 
   /**
    * Non-deprecated version of the shape() function. The
-   * Elem::p_level() is accounted for internally.
+   * Elem::p_level() is accounted for internally if \p add_p_level
    */
   static Real shape(const FEType & fe_t,
                     const Elem * elem,
                     const unsigned int i,
-                    const Point & p);
+                    const Point & p,
+                    const bool add_p_level = true);
 
   /**
    * Non-deprecated version of the shape() function. The

include/fe/fe_macro.h

@@ -51,7 +51,7 @@
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::init_shape_functions(const std::vector<Point> &, const Elem *); \
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::init_dual_shape_functions(unsigned int, unsigned int); \
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_all_shape_derivs (const Elem * elem, const Order o, const std::vector<Point> & p, std::vector<std::vector<Real>> * comps[3], const bool add_p_level); \
-  template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_side_nodal_soln(const Elem * elem, const Order o, const unsigned int side, const std::vector<Number> & elem_soln, std::vector<Number> & nodal_soln_on_side)
+  template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_side_nodal_soln(const Elem * elem, const Order o, const unsigned int side, const std::vector<Number> & elem_soln, std::vector<Number> & nodal_soln_on_side, bool add_p_level)
 
 #else // LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -64,7 +64,7 @@
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::init_shape_functions(const std::vector<Point> &, const Elem *); \
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::init_dual_shape_functions(unsigned int, unsigned int); \
   template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_all_shape_derivs (const Elem * elem, const Order o, const std::vector<Point> & p, std::vector<std::vector<Real>> * comps[3], const bool add_p_level); \
-  template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_side_nodal_soln(const Elem * elem, const Order o, const unsigned int side, const std::vector<Number> & elem_soln, std::vector<Number> & nodal_soln_on_side)
+  template LIBMESH_EXPORT void         FE<2,SUBDIVISION>::default_side_nodal_soln(const Elem * elem, const Order o, const unsigned int side, const std::vector<Number> & elem_soln, std::vector<Number> & nodal_soln_on_side, bool add_p_level)
 
 #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -116,8 +116,9 @@ template <>                                                             \
 void FE<_dim,_fetype>::nodal_soln(const Elem * elem,                    \
                                   const Order order,                    \
                                   const std::vector<Number> & elem_soln,\
-                                  std::vector<Number> & nodal_soln)     \
-{ _funcname(elem, order, elem_soln, nodal_soln); }
+                                  std::vector<Number> & nodal_soln,     \
+                                  const bool add_p_level)               \
+{ _funcname(elem, order, elem_soln, nodal_soln, add_p_level); }
 
 #define LIBMESH_FE_NODAL_SOLN(fetype, _funcname)                        \
 LIBMESH_FE_NODAL_SOLN_DIM(fetype, _funcname, 0)                         \
@@ -132,8 +133,9 @@ void FE<_dim,_fetype>::side_nodal_soln(const Elem * elem,               \
                                        const Order order,               \
                                        const unsigned int side,         \
                                        const std::vector<Number> & elem_soln,\
-                                       std::vector<Number> & nodal_soln)\
-{ default_side_nodal_soln(elem, order, side, elem_soln, nodal_soln); }
+                                       std::vector<Number> & nodal_soln, \
+                                       const bool add_p_level)          \
+{ default_side_nodal_soln(elem, order, side, elem_soln, nodal_soln, add_p_level); }
 
 #define LIBMESH_FE_SIDE_NODAL_SOLN(fetype)                              \
 LIBMESH_FE_SIDE_NODAL_SOLN_DIM(fetype, 0)                               \

include/systems/generic_projector.h

@@ -1656,7 +1656,7 @@ void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::SortAndCopy
                       }
             }
 
-          if (FEInterface::n_dofs_per_elem(fe_type, elem) ||
+          if (FEInterface::n_dofs_per_elem(fe_type, elem, add_p_level) ||
               (has_interior_nodes &&
                FEInterface::n_dofs_at_node(fe_type, elem, n_nodes-1, add_p_level)))
             {
@@ -2549,10 +2549,12 @@ 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);
 
           // This may be a p refined element
           fe_type.order =
-            libMesh::Order (fe_type.order + elem.p_level());
+            libMesh::Order (fe_type.order + add_p_level*elem.p_level());
 
           // If this is a Lagrange element with DoFs on sides then by
           // convention we interpolate at the node rather than project
@@ -2664,7 +2666,7 @@ void GenericProjector<FFunctor, GFunctor, FValue, ProjectionAction>::ProjectSide
 
           std::vector<unsigned int> side_dofs;
           FEInterface::dofs_on_side(&elem, dim, base_fe_type,
-                                    context.side, side_dofs);
+                                    context.side, side_dofs, add_p_level);
 
           this->construct_projection
             (dof_indices, side_dofs, var_component,
@@ -2708,10 +2710,13 @@ 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);
 
           // This may be a p refined element
           fe_type.order =
-            libMesh::Order (fe_type.order + elem->p_level());
+            libMesh::Order (fe_type.order + add_p_level * elem->p_level());
 
           const unsigned int var_component =
             system.variable_scalar_number(var, 0);

src/base/dof_map.C

@@ -813,7 +813,7 @@ void DofMap::reinit(MeshBase & mesh)
             }
           // Allocate the element DOFs
           const unsigned int dofs_per_elem =
-            FEInterface::n_dofs_per_elem(base_fe_type, elem);
+            FEInterface::n_dofs_per_elem(base_fe_type, elem, add_p_level);
 
           elem->set_n_comp_group(sys_num, vg, dofs_per_elem);
 
@@ -2509,7 +2509,7 @@ void DofMap::_dof_indices (const Elem & elem,
         }
 
       // If there are any element-based DOF numbers, get them
-      const unsigned int nc = FEInterface::n_dofs_per_elem(fe_type, p_level, &elem);
+      const unsigned int nc = FEInterface::n_dofs_per_elem(fe_type, add_p_level*p_level, &elem);
 
       // We should never have fewer dofs than necessary on an
       // element unless we're getting indices on a parent element

src/fe/fe.C

@@ -706,10 +706,11 @@ void
 FE<Dim,T>::default_side_nodal_soln(const Elem * elem, const Order o,
                                    const unsigned int side,
                                    const std::vector<Number> & elem_soln,
-                                   std::vector<Number> & nodal_soln_on_side)
+                                   std::vector<Number> & nodal_soln_on_side,
+                                   const bool add_p_level)
 {
   std::vector<Number> full_nodal_soln;
-  nodal_soln(elem, o, elem_soln, full_nodal_soln);
+  nodal_soln(elem, o, elem_soln, full_nodal_soln, add_p_level);
   const std::vector<unsigned int> side_nodes =
     elem->nodes_on_side(side);
 

src/fe/fe_base.C

@@ -1549,10 +1549,12 @@ 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);
 
   // Construct FE objects for this element and its neighbors.
   std::unique_ptr<FEGenericBase<OutputShape>> my_fe
     (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
+  my_fe->add_p_level_in_reinit(add_p_level);
   const FEContinuity cont = my_fe->get_continuity();
 
   // We don't need to constrain discontinuous elements
@@ -1570,6 +1572,7 @@ FEGenericBase<OutputType>::compute_proj_constraints (DofConstraints & constraint
 
   std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
     (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
+  neigh_fe->add_p_level_in_reinit(add_p_level);
 
   QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
   my_fe->attach_quadrature_rule (&my_qface);
@@ -1630,15 +1633,14 @@ FEGenericBase<OutputType>::compute_proj_constraints (DofConstraints & constraint
           // matrix with this and then constrain away all higher p
           // DoFs.
           libmesh_assert(neigh->active());
-          const unsigned int min_p_level =
+          const unsigned int min_p_level = add_p_level *
             std::min(elem->p_level(), neigh->p_level());
-
           // we may need to make the FE objects reinit with the
           // minimum shared p_level
-          const unsigned int old_elem_level = elem->p_level();
-          if (elem->p_level() != min_p_level)
+          const unsigned int old_elem_level = add_p_level * elem->p_level();
+          if (old_elem_level != min_p_level)
             my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + min_p_level - old_elem_level);
-          const unsigned int old_neigh_level = neigh->p_level();
+          const unsigned int old_neigh_level = add_p_level * neigh->p_level();
           if (old_neigh_level != min_p_level)
             neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + min_p_level - old_neigh_level);
 
@@ -1808,20 +1810,23 @@ FEGenericBase<OutputType>::compute_proj_constraints (DofConstraints & constraint
           neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + old_neigh_level - min_p_level);
         }
 
-      // p refinement constraints:
-      // constrain dofs shared between
-      // active elements and neighbors with
-      // lower polynomial degrees
-      const unsigned int min_p_level =
-        neigh->min_p_level_by_neighbor(elem, elem->p_level());
-      if (min_p_level < elem->p_level())
+      if (add_p_level)
+      {
+        // p refinement constraints:
+        // constrain dofs shared between
+        // active elements and neighbors with
+        // lower polynomial degrees
+        const unsigned int min_p_level =
+          neigh->min_p_level_by_neighbor(elem, elem->p_level());
+        if (min_p_level < elem->p_level())
         {
           // Adaptive p refinement of non-hierarchic bases will
           // require more coding
           libmesh_assert(my_fe->is_hierarchic());
           dof_map.constrain_p_dofs(variable_number, elem,
                                    s, min_p_level);
         }
+      }
     }
 }
 

src/fe/fe_bernstein.C

@@ -40,15 +40,16 @@ namespace {
 void bernstein_nodal_soln(const Elem * elem,
                           const Order order,
                           const std::vector<Number> & elem_soln,
-                          std::vector<Number> & nodal_soln)
+                          std::vector<Number> & nodal_soln,
+                          const bool add_p_level)
 {
   const unsigned int n_nodes = elem->n_nodes();
 
   const ElemType elem_type = elem->type();
 
   nodal_soln.resize(n_nodes);
 
-  const Order totalorder = static_cast<Order>(order + elem->p_level());
+  const Order totalorder = static_cast<Order>(order + add_p_level*elem->p_level());
 
   // FEType object to be passed to various FEInterface functions below.
   FEType fe_type(order, BERNSTEIN);