Adjusted settings for default FESeries object.

diff --git a/include/deal.II/numerics/smoothness_estimator.h b/include/deal.II/numerics/smoothness_estimator.h
index 4150e6c6172a8eb6aba30a2efdc858369041f5db..7aa7200d62a156bfff60a0704ce2fb8206d5a7f1 100644 (file)
--- a/include/deal.II/numerics/smoothness_estimator.h
+++ b/include/deal.II/numerics/smoothness_estimator.h
@@ -454,10 +454,10 @@ namespace SmoothnessEstimator
      * the default configuration for smoothness estimation purposes.
      *
      * For each finite element of the provided @p fe_collection, we use as many
-     * modes as its polynomial degree plus one, and at least three modes.
-     * Further for each element, we use a 4-point Gaussian quarature iterated in
-     * each dimension by the maximal wave number, which is the number of modes
-     * decreased by one since we start with $k = 0$.
+     * modes as its polynomial degree plus two. Further for each element, we use
+     * a 5-point Gaussian quarature iterated in each dimension by the maximal
+     * wave number, which is the number of modes decreased by one since we start
+     * with $k = 0$.
      */
     template <int dim, int spacedim>
     FESeries::Fourier<dim, spacedim>

diff --git a/source/numerics/smoothness_estimator.cc b/source/numerics/smoothness_estimator.cc
index 39b8f88cf3de8e3baab1dfa21595ed855f997593..82c651701c9f480e9f7e069bbdd057e0ccc22310 100644 (file)
--- a/source/numerics/smoothness_estimator.cc
+++ b/source/numerics/smoothness_estimator.cc
@@ -573,10 +573,15 @@ namespace SmoothnessEstimator
     default_fe_series(const hp::FECollection<dim, spacedim> &fe_collection)
     {
       // Default number of coefficients per direction.
+      //
+      // Since we omit the zero-th mode in the Fourier decay strategy, make sure
+      // that we have at least two modes to work with per finite element. With a
+      // number of modes equal to the polynomial degree plus two for each finite
+      // element, the smoothness estimation algorithm tends to produce stable
+      // results.
       std::vector<unsigned int> n_coefficients_per_direction;
       for (unsigned int i = 0; i < fe_collection.size(); ++i)
-        n_coefficients_per_direction.push_back(
-          std::max<unsigned int>(3, fe_collection[i].degree + 1));
+        n_coefficients_per_direction.push_back(fe_collection[i].degree + 2);
 
       // Default quadrature collection.
       //
@@ -588,11 +593,11 @@ namespace SmoothnessEstimator
       // elements we deal with, i.e. the matrices F_k,j. We have to do that for
       // each of the finite elements in use. To that end we need a quadrature
       // rule. As a default, we use the same quadrature formula for each finite
-      // element, namely one that is obtained by iterating a 4-point Gauss
+      // element, namely one that is obtained by iterating a 5-point Gauss
       // formula as many times as the maximal exponent we use for the term
       // exp(ikx). Since the first mode corresponds to k = 0, the maximal wave
       // number is k = n_modes - 1.
-      const QGauss<1>      base_quadrature(4);
+      const QGauss<1>      base_quadrature(5);
       hp::QCollection<dim> q_collection;
       for (unsigned int i = 0; i < fe_collection.size(); ++i)
         {

diff --git a/tests/numerics/smoothness_estimator_02.cc b/tests/numerics/smoothness_estimator_02.cc
index 9f048f17b62563aa0a5e955814c9878ce874d4a5..996db5739e93c1e42e0602943085719844a9f8c9 100644 (file)
--- a/tests/numerics/smoothness_estimator_02.cc
+++ b/tests/numerics/smoothness_estimator_02.cc
@@ -189,7 +189,7 @@ template <int dim>
 void
 test(const unsigned int poly_degree)
 {
-  const unsigned int    max_poly = 3;
+  const unsigned int    max_poly = 2;
   hp::FECollection<dim> fe_collection;
   for (unsigned int p = 1; p <= max_poly; ++p)
     fe_collection.push_back(FE_Q<dim>(p));
@@ -202,7 +202,7 @@ test(const unsigned int poly_degree)
     fourier.get_n_coefficients_per_direction(fe_index);
 
   Assert((poly_degree >= 1) && (poly_degree <= max_poly), ExcInternalError());
-  Assert((n_modes >= 3) && (n_modes <= max_poly + 1), ExcInternalError());
+  Assert((n_modes >= 3) && (n_modes <= max_poly + 2), ExcInternalError());
 
   deallog << "-----------------------------------" << std::endl;
   deallog << dim << "d, p=" << poly_degree << ", max_p=" << max_poly
@@ -314,7 +314,7 @@ main()
   dealii::deallog.attach(logfile, /*do not print job id*/ false);
   dealii::deallog.depth_console(0);
 
-  for (unsigned int poly_degree = 1; poly_degree <= 3; ++poly_degree)
+  for (unsigned int poly_degree = 1; poly_degree <= 2; ++poly_degree)
     {
       test<1>(poly_degree);
       test<2>(poly_degree);

diff --git a/tests/numerics/smoothness_estimator_02.output b/tests/numerics/smoothness_estimator_02.output
index ff88be2087f76f8472da61305f7b4ccc0fef15e4..1d9308dc07ba6ed83aafd8d3ff210b994bbd4238 100644 (file)
--- a/tests/numerics/smoothness_estimator_02.output
+++ b/tests/numerics/smoothness_estimator_02.output
@@ -1,64 +1,43 @@
 DEAL::-----------------------------------
-DEAL::1d, p=1, max_p=3, n_modes=3
+DEAL::1d, p=1, max_p=2, n_modes=3
 DEAL::-----------------------------------
 DEAL::Linf norm in exact and calculate Fourier coefficients:
-DEAL::0.000831516
+DEAL::2.80948e-05
 DEAL::estimated regularity:
-DEAL::0.988446
+DEAL::1.00043
 DEAL::-----------------------------------
-DEAL::2d, p=1, max_p=3, n_modes=3
+DEAL::2d, p=1, max_p=2, n_modes=3
 DEAL::-----------------------------------
 DEAL::Linf norm in exact and calculate Fourier coefficients:
-DEAL::0.000415758
+DEAL::1.40474e-05
 DEAL::estimated regularity:
-DEAL::1.09966
+DEAL::1.11550
 DEAL::-----------------------------------
-DEAL::3d, p=1, max_p=3, n_modes=3
+DEAL::3d, p=1, max_p=2, n_modes=3
 DEAL::-----------------------------------
 DEAL::Linf norm in exact and calculate Fourier coefficients:
-DEAL::0.000207879
+DEAL::7.02369e-06
 DEAL::estimated regularity:
-DEAL::0.890459
+DEAL::0.905620
 DEAL::-----------------------------------
-DEAL::1d, p=2, max_p=3, n_modes=3
+DEAL::1d, p=2, max_p=2, n_modes=4
 DEAL::-----------------------------------
 DEAL::Linf norm in exact and calculate Fourier coefficients:
-DEAL::0.000637547
+DEAL::1.59345e-05
 DEAL::estimated regularity:
-DEAL::1.04031
+DEAL::1.04044
 DEAL::-----------------------------------
-DEAL::2d, p=2, max_p=3, n_modes=3
+DEAL::2d, p=2, max_p=2, n_modes=4
 DEAL::-----------------------------------
 DEAL::Linf norm in exact and calculate Fourier coefficients:
-DEAL::0.000212516
+DEAL::5.31151e-06
 DEAL::estimated regularity:
-DEAL::0.854283
+DEAL::0.712977
 DEAL::-----------------------------------
-DEAL::3d, p=2, max_p=3, n_modes=3
+DEAL::3d, p=2, max_p=2, n_modes=4
 DEAL::-----------------------------------
 DEAL::Linf norm in exact and calculate Fourier coefficients:
-DEAL::7.08386e-05
+DEAL::1.77050e-06
 DEAL::estimated regularity:
-DEAL::0.504717
-DEAL::-----------------------------------
-DEAL::1d, p=3, max_p=3, n_modes=4
-DEAL::-----------------------------------
-DEAL::Linf norm in exact and calculate Fourier coefficients:
-DEAL::0.000375709
-DEAL::estimated regularity:
-DEAL::0.972651
-DEAL::-----------------------------------
-DEAL::2d, p=3, max_p=3, n_modes=4
-DEAL::-----------------------------------
-DEAL::Linf norm in exact and calculate Fourier coefficients:
-DEAL::9.39272e-05
-DEAL::estimated regularity:
-DEAL::0.522683
-DEAL::-----------------------------------
-DEAL::3d, p=3, max_p=3, n_modes=4
-DEAL::-----------------------------------
-DEAL::Linf norm in exact and calculate Fourier coefficients:
-DEAL::2.34818e-05
-DEAL::estimated regularity:
-DEAL::0.417735
+DEAL::0.706966
 DEAL::Ok