Merge pull request #24180 from YaqiWang/enable_p_refinement_24141

Enable p refinement

framework/include/base/Adaptivity.h

@@ -50,7 +50,7 @@ class Adaptivity : public ConsoleStreamInterface,
                    public libMesh::ParallelObject
 {
 public:
-  Adaptivity(FEProblemBase & subproblem);
+  Adaptivity(FEProblemBase & fe_problem);
   virtual ~Adaptivity();
 
   /**
@@ -128,6 +128,11 @@ class Adaptivity : public ConsoleStreamInterface,
    */
   void setRecomputeMarkersFlag(const bool flag) { _recompute_markers_during_cycles = flag; }
 
+  /**
+   * Switch from h-refinement to p-refinement
+   */
+  void switchHToPRefinement();
+
   /**
    * Adapts the mesh based on the error estimator used
    *
@@ -245,7 +250,7 @@ class Adaptivity : public ConsoleStreamInterface,
   bool isAdaptivityDue();
 
 protected:
-  FEProblemBase & _subproblem;
+  FEProblemBase & _fe_problem;
   MooseMesh & _mesh;
 
   /// on/off flag reporting if the adaptivity is being used
@@ -302,6 +307,8 @@ class Adaptivity : public ConsoleStreamInterface,
 
   /// Stores pointers to ErrorVectors associated with indicator field names
   std::map<std::string, std::unique_ptr<ErrorVector>> _indicator_field_to_error_vector;
+
+  bool _p_refinement_flag = false;
 };
 
 template <typename T>

framework/include/base/Assembly.h

@@ -568,6 +568,11 @@ class Assembly
    */
   bool needDual() const { return _need_dual; }
 
+  /**
+   * Set the cached quadrature rules to nullptr
+   */
+  void clearCachedQRules();
+
 private:
   /**
    * Set the qrule to be used for lower dimensional integration.
@@ -595,6 +600,11 @@ class Assembly
    */
   void reinit(const Elem * elem);
 
+  /**
+   * Set the volumetric quadrature rule based on the provided element
+   */
+  void setVolumeQRule(const Elem * elem);
+
   /**
    * Reinitialize FE data for the given element on the given side, optionally
    * with a given set of reference points
@@ -664,6 +674,11 @@ class Assembly
    */
   void reinit(const Elem * elem, const std::vector<Point> & reference_points);
 
+  /**
+   * Set the face quadrature rule based on the provided element and side
+   */
+  void setFaceQRule(const Elem * const elem, const unsigned int side);
+
   /**
    * Reinitialize the assembly data on an side of an element
    */
@@ -1817,6 +1832,11 @@ class Assembly
    */
   const Elem * const & msmElem() const { return _msm_elem; }
 
+  /**
+   * Indicate that we have p-refinement
+   */
+  void havePRefinement();
+
 private:
   /**
    * Just an internal helper function to reinit the volume FE objects.
@@ -2150,6 +2170,11 @@ class Assembly
     return _jacobian_block_nonlocal_used[tag][ivar][_block_diagonal_matrix ? 0 : jvar];
   }
 
+  /**
+   * request phi, dphi, xyz, JxW, etc. data through the FE helper functions
+   */
+  void helpersRequestData();
+
   SystemBase & _sys;
   SubProblem & _subproblem;
 
@@ -2194,6 +2219,28 @@ class Assembly
 
   unsigned int _mesh_dimension;
 
+  /// The finite element type of the FE helper classes. The helper class gives us data like JxW, the
+  /// physical quadrature point locations, etc.
+  const FEType _helper_type;
+
+  /// Whether user code requested a \p FEType the same as our \p _helper_type
+  mutable bool _user_added_fe_of_helper_type;
+  mutable bool _user_added_fe_face_of_helper_type;
+  mutable bool _user_added_fe_face_neighbor_of_helper_type;
+  mutable bool _user_added_fe_neighbor_of_helper_type;
+  mutable bool _user_added_fe_lower_of_helper_type;
+
+  /// Containers for holding unique FE helper types if we are doing p-refinement. If we are not
+  /// doing p-refinement then the helper data is owned by the \p _fe data members
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_helper;
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_face_helper;
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_face_neighbor_helper;
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_neighbor_helper;
+  std::vector<std::unique_ptr<FEBase>> _unique_fe_lower_helper;
+
+  /// Whether we are currently building the FE classes for the helpers
+  bool _building_helpers;
+
   /// The XFEM controller
   std::shared_ptr<XFEMInterface> _xfem;
 
@@ -2222,7 +2269,7 @@ class Assembly
   /// Each dimension's actual vector fe objects indexed on type
   mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe;
   /// Each dimension's helper objects
-  std::map<unsigned int, FEBase **> _holder_fe_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_helper;
   /// The current helper object for transforming coordinates
   FEBase * _current_fe_helper;
   /// The current current quadrature rule being used (could be either volumetric or arbitrary - for dirac kernels)
@@ -2334,7 +2381,7 @@ class Assembly
   /// types of vector finite elements
   mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_face;
   /// Each dimension's helper objects
-  std::map<unsigned int, FEBase **> _holder_fe_face_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_face_helper;
   /// helper object for transforming coordinates
   FEBase * _current_fe_face_helper;
   /// quadrature rule used on faces
@@ -2370,15 +2417,15 @@ class Assembly
   mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_face_neighbor;
 
   /// Each dimension's helper objects
-  std::map<unsigned int, FEBase **> _holder_fe_neighbor_helper;
-  std::map<unsigned int, FEBase **> _holder_fe_face_neighbor_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_neighbor_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_face_neighbor_helper;
 
   /// FE objects for lower dimensional elements
   mutable std::map<unsigned int, std::map<FEType, FEBase *>> _fe_lower;
   /// Vector FE objects for lower dimensional elements
   mutable std::map<unsigned int, std::map<FEType, FEVectorBase *>> _vector_fe_lower;
   /// helper object for transforming coordinates for lower dimensional element quadrature points
-  std::map<unsigned int, FEBase **> _holder_fe_lower_helper;
+  std::map<unsigned int, FEBase *> _holder_fe_lower_helper;
 
   /// quadrature rule used on neighbors
   QBase * _current_qrule_neighbor;

framework/include/loops/MaxQpsThread.h

@@ -38,16 +38,11 @@ class MaxQpsThread
 
   unsigned int max() const { return _max; }
 
-  unsigned int max_shape_funcs() const { return _max_shape_funcs; }
-
 protected:
   FEProblemBase & _fe_problem;
 
   THREAD_ID _tid;
 
   /// Maximum number of qps encountered
   unsigned int _max;
-
-  /// Maximum number of shape functions encountered
-  unsigned int _max_shape_funcs;
 };

framework/include/problems/DisplacedProblem.h

@@ -354,6 +354,11 @@ class DisplacedProblem : public SubProblem
 
   virtual const std::vector<VectorTag> & currentResidualVectorTags() const override;
 
+  /**
+   * Indicate that we have p-refinement
+   */
+  void havePRefinement();
+
 protected:
   FEProblemBase & _mproblem;
   MooseMesh & _mesh;

framework/include/problems/FEProblemBase.h

@@ -339,11 +339,6 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   unsigned int getMaxQps() const;
 
-  /**
-   * @return The maximum number of quadrature points in use on any element in this problem.
-   */
-  unsigned int getMaxShapeFunctions() const;
-
   /**
    * @return The maximum order for all scalar variables in this problem's systems.
    */
@@ -2102,6 +2097,11 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   void clearCurrentResidualVectorTags();
 
+  /**
+   * Indicate that we have p-refinement
+   */
+  void havePRefinement();
+
 protected:
   /// Create extra tagged vectors and matrices
   void createTagVectors();
@@ -2384,9 +2384,6 @@ class FEProblemBase : public SubProblem, public Restartable
   /// Maximum number of quadrature points used in the problem
   unsigned int _max_qps;
 
-  /// Maximum number of shape functions on any element in the problem
-  unsigned int _max_shape_funcs;
-
   /// Maximum scalar variable order
   Order _max_scalar_order;
 

framework/src/actions/AdaptivityAction.C

@@ -79,6 +79,7 @@ AdaptivityAction::validParams()
       "show_initial_progress", true, "Show the progress of the initial adaptivity");
   params.addParam<bool>(
       "recompute_markers_during_cycles", false, "Recompute markers during adaptivity cycles");
+  params.addParam<bool>("switch_h_to_p_refinement", false, "True to perform p-refinement");
   return params;
 }
 
@@ -202,6 +203,8 @@ AdaptivityAction::act()
 
     adapt.setTimeActive(getParam<Real>("start_time"), getParam<Real>("stop_time"));
     adapt.setInterval(getParam<unsigned int>("interval"));
+    if (getParam<bool>("switch_h_to_p_refinement"))
+      adapt.switchHToPRefinement();
   }
 }
 

framework/src/actions/SetAdaptivityOptionsAction.C

@@ -49,6 +49,7 @@ SetAdaptivityOptionsAction::validParams()
       "The number of adaptive steps to use when on each timestep during a Transient simulation.");
   params.addParam<bool>(
       "recompute_markers_during_cycles", false, "Recompute markers during adaptivity cycles");
+  params.addParam<bool>("switch_h_to_p_refinement", false, "True to perform p-refinement");
   return params;
 }
 
@@ -137,5 +138,7 @@ SetAdaptivityOptionsAction::act()
     adapt.setInterval(getParam<unsigned int>("interval"));
 
     adapt.setRecomputeMarkersFlag(getParam<bool>("recompute_markers_during_cycles"));
+    if (getParam<bool>("switch_h_to_p_refinement"))
+      adapt.switchHToPRefinement();
   }
 }