hp-Adaptivity: smoothness percentage control flexibility

Author: skvarjun

examples/adaptivity/adaptivity_ex3/adaptivity_ex3.C

@@ -150,6 +150,7 @@ int main(int argc, char ** argv)
   const Real coarsen_percentage     = input_file("coarsen_percentage", 0.5);
   const unsigned int uniform_refine = input_file("uniform_refine",0);
   const std::string refine_type     = input_file("refinement_type", "h");
+  const Real smoothness_threshold   = input_file("smoothness_threshold", 0.9);
   const std::string approx_type     = input_file("approx_type", "LAGRANGE");
   const unsigned int approx_order   = input_file("approx_order", 1);
   n_vars                            = input_file("n_vars", n_vars);
@@ -523,7 +524,7 @@ int main(int argc, char ** argv)
               else if (refine_type == "hp")
                 {
                   HPCoarsenTest hpselector;
-                  hpselector.select_refinement(system);
+                  hpselector.select_refinement(system, smoothness_threshold);
                 }
               // If we are doing "singular hp" refinement, we
               // try switching most elements from h to p

examples/adaptivity/adaptivity_ex3/adaptivity_ex3.in

@@ -39,6 +39,12 @@ indicator_type = kelly
 # on every element; hp refinement tries to do one or the other
 refinement_type = h
 
+# Fraction of smoothness values above which p-refinement is
+# favored in hp-adaptivity. For instance, default value is 0.9;
+# means p refining favored (for hp adaptivity) on an element
+# only if local smoothness is above 90%.
+smoothness_threshold = 0.9
+
 # Mesh dimensionality
 # We generate 1D meshes on [-1, 0].
 # We load lshaped.xda for 2 dimensional problems, lshaped3D.xda for 3

include/error_estimation/hp_coarsentest.h

@@ -97,7 +97,7 @@ class HPCoarsenTest : public HPSelector
    * refinement and potentially change the desired
    * refinement type.
    */
-  virtual void select_refinement (System & system) override;
+  virtual void select_refinement (System & system, Real smoothness_threshold = 0.9) override;
 
   /**
    * Because the coarsening test seems to always choose p refinement, we're

include/error_estimation/hp_selector.h

@@ -70,7 +70,7 @@ class HPSelector
    * refinement and potentially change the desired
    * refinement type.
    */
-  virtual void select_refinement (System & system) = 0;
+  virtual void select_refinement (System & system, Real smoothness_threshold = 0.9) = 0;
 
   /**
    * This vector can be used to "scale" certain

src/error_estimation/hp_coarsentest.C

@@ -135,7 +135,7 @@ void HPCoarsenTest::add_projection(const System & system,
     }
 }
 
-void HPCoarsenTest::select_refinement (System & system)
+void HPCoarsenTest::select_refinement (System & system, Real smoothness_threshold)
 {
   LOG_SCOPE("select_refinement()", "HPCoarsenTest");
 
@@ -541,7 +541,11 @@ void HPCoarsenTest::select_refinement (System & system)
       const Real h_value =
         std::sqrt(h_error_per_cell[e_id]) /
         static_cast<Real>(new_h_dofs);
-      if (smoothness[elem->id()] > smoothness.mean() && p_value > h_value)
+
+      const Real _threshold = smoothness.minimum() +
+        smoothness_threshold * (smoothness.maximum() - smoothness.minimum());
+
+      if (smoothness[elem->id()] > _threshold && p_value > h_value)
         {
           elem->set_p_refinement_flag(Elem::REFINE);
           elem->set_refinement_flag(Elem::DO_NOTHING);