]> https://gitweb.dealii.org/ - dealii.git/commitdiff
SmoothnessEstimator: Initial implementation with FESeries template.
authorMarc Fehling <marc.fehling@gmx.net>
Wed, 21 Nov 2018 10:07:00 +0000 (11:07 +0100)
committerMarc Fehling <marc.fehling@gmx.net>
Tue, 12 May 2020 14:01:47 +0000 (16:01 +0200)
cmake/config/template-arguments.in
examples/step-27/doc/intro.dox
examples/step-27/step-27.cc
include/deal.II/hp/fe_collection.h
include/deal.II/numerics/smoothness_estimator.h [new file with mode: 0644]
include/deal.II/numerics/smoothness_estimator.templates.h [new file with mode: 0644]
source/numerics/CMakeLists.txt
source/numerics/smoothness_estimator.cc [new file with mode: 0644]
source/numerics/smoothness_estimator.inst.in [new file with mode: 0644]
tests/hp/step-27.cc

index efe2f5111ee9fb2d7ada4e0c3551dcc27ed085a8..836fb7fc33fc425440a6e1792e08b0ef55751764 100644 (file)
@@ -258,6 +258,11 @@ AFFINE_CONSTRAINTS_SP_BLOCK := { BlockSparsityPattern;
                                  @DEAL_II_EXPAND_TRILINOS_BLOCK_SPARSITY_PATTERN@;
                                }
 
+// Series expansion templates
+SERIES_EXPANSION_TEMPLATES := { FESeries::Fourier;
+                                FESeries::Legendre;
+                              }
+
 // all supported logical dimensions
 DIMENSIONS := { 1; 2; 3 }
 
index 0597da67891a27f40f6045b82c60b51977435c43..e30a733ca081f5bcf89bd6f9f1432601703cb846 100644 (file)
@@ -367,7 +367,7 @@ and see how fast they decay. If they decay as
 then consequently the function we had here was in $H^{\mu-d/2}$.
 
 
-<h4>What we have to do</h4>
+<h4>What has to be done</h4>
 
 So what do we have to do to estimate the local smoothness of $u({\bf x})$ on
 a cell $K$? Clearly, the first step is to compute the Fourier coefficients
@@ -403,9 +403,9 @@ cell. In other words, we can write it as a matrix-vector product
 @f]
 with the matrix
 @f[
-        {\cal F}_{{\bf k},j}
-        =
-        \int_{\hat K} e^{i {\bf k}\cdot \hat{\bf x}} \hat \varphi_j(\hat{\bf x}) d\hat{\bf x}.
+   {\cal F}_{{\bf k},j}
+   =
+   \int_{\hat K} e^{i {\bf k}\cdot \hat{\bf x}} \hat \varphi_j(\hat{\bf x}) d\hat{\bf x}.
 @f]
 This matrix is easily computed for a given number of shape functions
 $\varphi_j$ and Fourier modes $N$. Consequently, finding the
@@ -509,6 +509,10 @@ $\beta$, the formula above gives us a mean to calculate the value of
 the exponent $\mu$ that we can then use to determine that
 $\hat u(\hat{\bf x})$ is in $H^s(\hat K)$ with $s=\mu-\frac d2$.
 
+These steps outlined above are applicable on many different scenarios, which
+motivated the introduction of a generic function
+SmoothnessEstimator::estimate_by_coeff_decay() in deal.II, that combines all
+the tasks described in this section in one simple function call.
 
 <h4>Compensating for anisotropy</h4>
 
index 54fbca5271969e9621c9e8642da347406b22b924..3d32f9638f3ae4808f41e40262e64385ea17597d 100644 (file)
 
 // These are the new files we need. The first and second provide <i>hp</i>
 // versions of the DoFHandler and FEValues classes as described in the
-// introduction of this program. The last one provides Fourier transformation
-// class on the unit cell.
+// introduction of this program. The last one provides the smoothness estimation
+// algorithms on decaying series expansion coefficients.
 #include <deal.II/hp/dof_handler.h>
 #include <deal.II/hp/fe_values.h>
-#include <deal.II/fe/fe_series.h>
+#include <deal.II/numerics/smoothness_estimator.h>
 
 // The last set of include files are standard C++ headers. We need support for
 // complex numbers when we compute the Fourier transform.
@@ -99,9 +99,7 @@ namespace Step27
     void assemble_system();
     void solve();
     void create_coarse_grid();
-    void estimate_smoothness(Vector<float> &smoothness_indicators);
     void postprocess(const unsigned int cycle);
-    std::pair<bool, unsigned int> predicate(const TableIndices<dim> &indices);
 
     Triangulation<dim> triangulation;
 
@@ -110,11 +108,6 @@ namespace Step27
     hp::QCollection<dim>     quadrature_collection;
     hp::QCollection<dim - 1> face_quadrature_collection;
 
-    hp::QCollection<dim>                    fourier_q_collection;
-    std::shared_ptr<FESeries::Fourier<dim>> fourier;
-    std::vector<double>                     ln_k;
-    Table<dim, std::complex<double>>        fourier_coefficients;
-
     AffineConstraints<double> constraints;
 
     SparsityPattern      sparsity_pattern;
@@ -167,23 +160,6 @@ namespace Step27
   // face quadrature objects. We start with quadratic elements, and each
   // quadrature formula is chosen so that it is appropriate for the matching
   // finite element in the hp::FECollection object.
-  //
-  // Finally, we initialize FESeries::Fourier object which will be used to
-  // calculate coefficient in Fourier series as described in the introduction.
-  // In addition to the hp::FECollection, we need to provide quadrature rules
-  // hp::QCollection for integration on the reference cell.
-  //
-  // In order to resize fourier_coefficients Table, we use the following
-  // auxiliary function
-  template <int dim, typename T>
-  void resize(Table<dim, T> &coeff, const unsigned int N)
-  {
-    TableIndices<dim> size;
-    for (unsigned int d = 0; d < dim; d++)
-      size[d] = N;
-    coeff.reinit(size);
-  }
-
   template <int dim>
   LaplaceProblem<dim>::LaplaceProblem()
     : dof_handler(triangulation)
@@ -195,48 +171,6 @@ namespace Step27
         quadrature_collection.push_back(QGauss<dim>(degree + 1));
         face_quadrature_collection.push_back(QGauss<dim - 1>(degree + 1));
       }
-
-    // As described in the introduction, we define the Fourier vectors ${\bf
-    // k}$ for which we want to compute Fourier coefficients of the solution
-    // on each cell as follows. In 2d, we will need coefficients corresponding
-    // to vectors ${\bf k}=(2 \pi i, 2\pi j)^T$ for which $\sqrt{i^2+j^2}\le N$,
-    // with $i,j$ integers and $N$ being the maximal polynomial degree we use
-    // for the finite elements in this program. The FESeries::Fourier class'
-    // constructor first parameter $N$ defines the number of coefficients in 1D
-    // with the total number of coefficients being $N^{dim}$. Although we will
-    // not use coefficients corresponding to
-    // $\sqrt{i^2+j^2}> N$ and $i+j==0$, the overhead of their calculation is
-    // minimal. The transformation matrices for each FiniteElement will be
-    // calculated only once the first time they are required in the course of
-    // hp-adaptive refinement. Because we work on the unit cell, we can do all
-    // this work without a mapping from reference to real cell and consequently
-    // can precalculate these matrices. The calculation of expansion
-    // coefficients for a particular set of local degrees of freedom on a given
-    // cell then follows as a simple matrix-vector product.
-    // The 3d case is handled analogously.
-    const unsigned int N = max_degree;
-
-    // We will need to assemble the matrices that do the Fourier transforms
-    // for each of the finite elements we deal with, i.e. the matrices ${\cal
-    // F}_{{\bf k},j}$ defined in the introduction. We have to do that for
-    // each of the finite elements in use. To that end we need a quadrature
-    // rule. In this example we use the same quadrature formula for each
-    // finite element, namely that is obtained by iterating a
-    // 2-point Gauss formula as many times as the maximal exponent we use for
-    // the term $e^{i{\bf k}\cdot{\bf x}}$:
-    QGauss<1>      base_quadrature(2);
-    QIterated<dim> quadrature(base_quadrature, N);
-    for (unsigned int i = 0; i < fe_collection.size(); i++)
-      fourier_q_collection.push_back(quadrature);
-
-    // Now we are ready to set-up the FESeries::Fourier object
-    fourier = std::make_shared<FESeries::Fourier<dim>>(N,
-                                                       fe_collection,
-                                                       fourier_q_collection);
-
-    // We need to resize the matrix of fourier coefficients according to the
-    // number of modes N.
-    resize(fourier_coefficients, N);
   }
 
 
@@ -403,8 +337,7 @@ namespace Step27
     // Let us start with computing estimated error and smoothness indicators,
     // which each are one number for each active cell of our
     // triangulation. For the error indicator, we use the KellyErrorEstimator
-    // class as always. Estimating the smoothness is done in the respective
-    // function of this class; that function is discussed further down below:
+    // class as always.
     Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
     KellyErrorEstimator<dim>::estimate(
       dof_handler,
@@ -413,9 +346,11 @@ namespace Step27
       solution,
       estimated_error_per_cell);
 
-
-    Vector<float> smoothness_indicators(triangulation.n_active_cells());
-    estimate_smoothness(smoothness_indicators);
+    // Estimating the smoothness is performed with the method of decaing
+    // expansion coefficients as outlined in the introduction.
+    Vector<float> smoothness_indicators;
+    SmoothnessEstimator::estimate_by_coeff_decay<FESeries::Fourier<dim>>(
+      dof_handler, solution, smoothness_indicators);
 
     // Next we want to generate graphical output. In addition to the two
     // estimated quantities derived above, we would also like to output the
@@ -613,118 +548,6 @@ namespace Step27
         postprocess(cycle);
       }
   }
-
-
-  // @sect4{LaplaceProblem::estimate_smoothness}
-
-  // As described in the introduction, we will need to take the maximum
-  // absolute value of fourier coefficients which correspond to $k$-vector
-  // $|{\bf k}|= const$. To filter the coefficients Table we
-  // will use the FESeries::process_coefficients() which requires a predicate
-  // to be specified. The predicate should operate on TableIndices and return
-  // a pair of <code>bool</code> and <code>unsigned int</code>. The latter
-  // is the value of the map from TableIndicies to unsigned int.  It is
-  // used to define subsets of coefficients from which we search for the one
-  // with highest absolute value, i.e. $l^\infty$-norm. The <code>bool</code>
-  // parameter defines which indices should be used in processing. In the
-  // current case we are interested in coefficients which correspond to
-  // $0 < i*i+j*j < N*N$ and $0 < i*i+j*j+k*k < N*N$ in 2D and 3D, respectively.
-  template <int dim>
-  std::pair<bool, unsigned int>
-  LaplaceProblem<dim>::predicate(const TableIndices<dim> &ind)
-  {
-    unsigned int v = 0;
-    for (unsigned int i = 0; i < dim; i++)
-      v += ind[i] * ind[i];
-    if (v > 0 && v < max_degree * max_degree)
-      return std::make_pair(true, v);
-    else
-      return std::make_pair(false, v);
-  }
-
-  // This last function of significance implements the algorithm to estimate
-  // the smoothness exponent using the algorithms explained in detail in the
-  // introduction. We will therefore only comment on those points that are of
-  // implementational importance.
-  template <int dim>
-  void
-  LaplaceProblem<dim>::estimate_smoothness(Vector<float> &smoothness_indicators)
-  {
-    // Since most of the hard work is done for us in FESeries::Fourier and
-    // we set up the object of this class in the constructor, what we are left
-    // to do here is apply this class to calculate coefficients and then
-    // perform linear regression to fit their decay slope.
-
-
-    // First thing to do is to loop over all cells and do our work there, i.e.
-    // to locally do the Fourier transform and estimate the decay coefficient.
-    // We will use the following array as a scratch array in the loop to store
-    // local DoF values:
-    Vector<double> local_dof_values;
-
-    // Then here is the loop:
-    for (const auto &cell : dof_handler.active_cell_iterators())
-      {
-        // Inside the loop, we first need to get the values of the local
-        // degrees of freedom (which we put into the
-        // <code>local_dof_values</code> array after setting it to the right
-        // size) and then need to compute the Fourier transform by multiplying
-        // this vector with the matrix ${\cal F}$ corresponding to this finite
-        // element. This is done by calling FESeries::Fourier::calculate(),
-        // that has to be provided with the <code>local_dof_values</code>,
-        // <code>cell->active_fe_index()</code> and a Table to store
-        // coefficients.
-        local_dof_values.reinit(cell->get_fe().dofs_per_cell);
-        cell->get_dof_values(solution, local_dof_values);
-
-        fourier->calculate(local_dof_values,
-                           cell->active_fe_index(),
-                           fourier_coefficients);
-
-        // The next thing, as explained in the introduction, is that we wanted
-        // to only fit our exponential decay of Fourier coefficients to the
-        // largest coefficients for each possible value of $|{\bf k}|$. To
-        // this end, we use FESeries::process_coefficients() to rework
-        // coefficients into the desired format. We'll only take those Fourier
-        // coefficients with the largest magnitude for a given value of $|{\bf
-        // k}|$ and thereby need to use VectorTools::Linfty_norm:
-        std::pair<std::vector<unsigned int>, std::vector<double>> res =
-          FESeries::process_coefficients<dim>(
-            fourier_coefficients,
-            [this](const TableIndices<dim> &indices) {
-              return this->predicate(indices);
-            },
-            VectorTools::Linfty_norm);
-
-        Assert(res.first.size() == res.second.size(), ExcInternalError());
-
-        // The first vector in the <code>std::pair</code> will store values of
-        // the predicate, that is $i*i+j*j= const$ or $i*i+j*j+k*k = const$ in
-        // 2D or 3D respectively. This vector will be the same for all the cells
-        // so we can calculate logarithms of the corresponding Fourier vectors
-        // $|{\bf k}|$ only once in the whole hp-refinement cycle:
-        if (ln_k.size() == 0)
-          {
-            ln_k.resize(res.first.size(), 0);
-            for (unsigned int f = 0; f < ln_k.size(); f++)
-              ln_k[f] =
-                std::log(2.0 * numbers::PI * std::sqrt(1. * res.first[f]));
-          }
-
-        // We have to calculate the logarithms of absolute values of
-        // coefficients and use it in a linear regression fit to obtain $\mu$.
-        for (double &residual_element : res.second)
-          residual_element = std::log(residual_element);
-
-        std::pair<double, double> fit =
-          FESeries::linear_regression(ln_k, res.second);
-
-        // The final step is to compute the Sobolev index $s=\mu-\frac d2$ and
-        // store it in the vector of estimated values for each cell:
-        smoothness_indicators(cell->active_cell_index()) =
-          -fit.first - 1. * dim / 2;
-      }
-  }
 } // namespace Step27
 
 
index 1764bffd5dc60027940ef5e4e8b3beb77df72ddd..746fea8972ca47835ad079945d1da6047308f65b 100644 (file)
@@ -227,6 +227,12 @@ namespace hp
     unsigned int
     n_blocks() const;
 
+    /**
+     * Return the maximal degree over all elements of this collection.
+     */
+    unsigned int
+    max_degree() const;
+
     /**
      * Return the maximal number of degrees of freedom per vertex over all
      * elements of this collection.
@@ -906,6 +912,22 @@ namespace hp
 
 
 
+  template <int dim, int spacedim>
+  unsigned int
+  FECollection<dim, spacedim>::max_degree() const
+  {
+    Assert(finite_elements.size() > 0, ExcNoFiniteElements());
+
+    unsigned int max = 0;
+    for (unsigned int i = 0; i < finite_elements.size(); ++i)
+      if (finite_elements[i]->degree > max)
+        max = finite_elements[i]->degree;
+
+    return max;
+  }
+
+
+
   template <int dim, int spacedim>
   unsigned int
   FECollection<dim, spacedim>::max_dofs_per_vertex() const
diff --git a/include/deal.II/numerics/smoothness_estimator.h b/include/deal.II/numerics/smoothness_estimator.h
new file mode 100644 (file)
index 0000000..9965cd8
--- /dev/null
@@ -0,0 +1,169 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_smoothness_estimator_h
+#define dealii_smoothness_estimator_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <vector>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/**
+ * Estimate the smoothness of a solution based on the decay of coefficients from
+ * a series expansion.
+ *
+ * From the definition, we can write our series expansion $\hat U_{\bf k}$ as a
+ * matrix product
+ * @f[
+ *    \hat U_{\bf k}
+ *    = {\cal F}_{{\bf k},j} u_j,
+ * @f]
+ * with $u_j$ the coefficients and ${\cal F}_{{\bf k},j}$ the transformation
+ * matrix from the expansion. We use the classes FESeries::Fourier and
+ * FESeries::Legendre to determine all coefficients $u_j$.
+ *
+ * The next step is that we have to estimate how fast these coefficients
+ * decay with $|{\bf k}|$. Thus, we perform a least-squares fit
+ * @f[
+ *    \min_{\alpha,\mu}
+ *    \frac 12 \sum_{{\bf k}, |{\bf k}|\le N}
+ *    \left( |\hat U_{\bf k}| - \alpha |{\bf k}|^{-\mu}\right)^2
+ * @f]
+ * with linear regressions coefficients $\alpha$ and $\mu$. For simplification,
+ * we apply a logarithm on our minimization problem
+ * @f[
+ *    \min_{\beta,\mu}
+ *    Q(\beta,\mu) =
+ *    \frac 12 \sum_{{\bf k}, |{\bf k}|\le N}
+ *    \left( \ln |\hat U_{\bf k}| - \beta + \mu \ln |{\bf k}|\right)^2,
+ * @f]
+ * where $\beta=\ln \alpha$. This is now a problem for which the
+ * optimality conditions $\frac{\partial Q}{\partial\beta}=0,
+ * \frac{\partial Q}{\partial\mu}=0$, are linear in $\beta,\mu$. We can
+ * write these conditions as follows:
+ * @f[
+ *    \left(\begin{array}{cc}
+ *    \sum_{{\bf k}, |{\bf k}|\le N} 1 &
+ *    \sum_{{\bf k}, |{\bf k}|\le N} \ln |{\bf k}|
+ *    \\
+ *    \sum_{{\bf k}, |{\bf k}|\le N} \ln |{\bf k}| &
+ *    \sum_{{\bf k}, |{\bf k}|\le N} (\ln |{\bf k}|)^2
+ *    \end{array}\right)
+ *    \left(\begin{array}{c}
+ *    \beta \\ -\mu
+ *    \end{array}\right)
+ *    =
+ *    \left(\begin{array}{c}
+ *    \sum_{{\bf k}, |{\bf k}|\le N} \ln |\hat U_{{\bf k}}|
+ *    \\
+ *    \sum_{{\bf k}, |{\bf k}|\le N} \ln |\hat U_{{\bf k}}| \ln |{\bf k}|
+ *    \end{array}\right)
+ * @f]
+ * Solving for $\beta$ and $\mu$ is nothing else but a linear regression fit and
+ * to do that we will use FESeries::linear_regression().
+ *
+ * While we are not particularly interested in the actual value of
+ * $\beta$, the formula above gives us a mean to calculate the value of
+ * the exponent $\mu$ that we can then use to determine that
+ * $\hat u(\hat{\bf x})$ is in $H^s(\hat K)$ with $s=\mu-\frac d2$. These
+ * Sobolev indices $s$ will suffice as our smoothness estimators and will be
+ * calculated on each cell for any provided solution.
+ *
+ * @note An extensive demonstration of the use of these functions is provided in step-27.
+ *
+ * @ingroup numerics
+ * @author Denis Davydov, 2016, Marc Fehling, 2018
+ */
+namespace SmoothnessEstimator
+{
+  /**
+   * Estimates the smoothness of the provided solutions using the method of
+   * decaying coefficents as outlined above.
+   *
+   * The @p regression_strategy parameter determines which norm will be used on the subset of
+   * coeffiecients $\mathbf{k}$ with the same absolute value $|\mathbf{k}|$.
+   * Default is VectorTools::Linfty_norm for a maximum approximation.
+   *
+   * Smoothness indicators will be calculated for each solution in @p all_solutions
+   * and stored in @p all_smoothness_indicators in the same order.
+   *
+   * An individual @p fe_series object can be supplied, which has to be constructed with the
+   * same FECollection object as the @p dof_handler.
+   */
+  template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+  void
+  estimate_by_coeff_decay(
+    FESeriesType &                         fe_series,
+    const DoFHandlerType &                 dof_handler,
+    const std::vector<const VectorType *> &all_solutions,
+    const std::vector<Vector<float> *> &   all_smoothness_indicators,
+    const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm);
+
+  /**
+   * Same as the function above, only for one @p solution vector.
+   */
+  template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+  void
+  estimate_by_coeff_decay(
+    FESeriesType &              fe_series,
+    const DoFHandlerType &      dof_handler,
+    const VectorType &          solution,
+    Vector<float> &             smoothness_indicators,
+    const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm);
+
+  /**
+   * Same as the function above, but with a default configuration for the chosen
+   * series expansion, using 2-point Gaussian quadrature for each finite
+   * element.
+   *
+   * Provide the desired series expansion as a template argument, i.e.
+   * @code
+   *    SmoothnessEstimator::estimate_by_coeff_decay<FESeries::Fourier<dim>>(
+   *      dof_handler, all_solutions, all_smoothness_indicators);
+   * @endcode
+   */
+  template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+  void
+  estimate_by_coeff_decay(
+    const DoFHandlerType &                 dof_handler,
+    const std::vector<const VectorType *> &all_solutions,
+    const std::vector<Vector<float> *> &   all_smoothness_indicators,
+    const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm);
+
+  /**
+   * Same as the function above, only for one @p solution vector.
+   */
+  template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+  void
+  estimate_by_coeff_decay(
+    const DoFHandlerType &      dof_handler,
+    const VectorType &          solution,
+    Vector<float> &             smoothness_indicators,
+    const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm);
+} // namespace SmoothnessEstimator
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/numerics/smoothness_estimator.templates.h b/include/deal.II/numerics/smoothness_estimator.templates.h
new file mode 100644 (file)
index 0000000..2570ddc
--- /dev/null
@@ -0,0 +1,262 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_smoothness_estimator_templates_h
+#define dealii_smoothness_estimator_templates_h
+
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/fe/fe_series.h>
+
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/numerics/smoothness_estimator.h>
+
+#include <algorithm>
+#include <cmath>
+#include <functional>
+#include <utility>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace SmoothnessEstimator
+{
+  namespace
+  {
+    /**
+     * Resizes @p coeff to @p N in each dimension.
+     */
+    template <int dim, typename CoefficientType>
+    void
+    resize(Table<dim, CoefficientType> &coeff, const unsigned int N)
+    {
+      TableIndices<dim> size;
+      for (unsigned int d = 0; d < dim; d++)
+        size[d] = N;
+      coeff.reinit(size);
+    }
+
+
+
+    /**
+     * Calculates predicates of @p ind in the form
+     * \f$
+     *    v = \sum\limits_{d=0}^{dim} ind[d]^2
+     * \f$.
+     *
+     * We flag the predicate whether it fulfills the criterion
+     * \f$
+     *    0 < v < max_degree^2
+     * \f$
+     * using @p max_degree.
+     */
+    template <int dim>
+    std::pair<bool, unsigned int>
+    predicate(const TableIndices<dim> &ind, const unsigned int max_degree)
+    {
+      unsigned int v = 0;
+      for (unsigned int i = 0; i < dim; i++)
+        v += ind[i] * ind[i];
+      if (v > 0 && v < max_degree * max_degree)
+        return std::make_pair(true, v);
+      else
+        return std::make_pair(false, v);
+    }
+  } // namespace
+
+
+
+  template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+  void
+  estimate_by_coeff_decay(
+    FESeriesType &                         fe_series,
+    const DoFHandlerType &                 dof_handler,
+    const std::vector<const VectorType *> &all_solutions,
+    const std::vector<Vector<float> *> &   all_smoothness_indicators,
+    const VectorTools::NormType            regression_strategy)
+  {
+    AssertDimension(all_solutions.size(), all_smoothness_indicators.size());
+
+    for (auto &smoothness_indicator : all_smoothness_indicators)
+      smoothness_indicator->reinit(
+        dof_handler.get_triangulation().n_active_cells());
+
+    const unsigned int dim = DoFHandlerType::dimension;
+    const unsigned int max_degree =
+      dof_handler.get_fe_collection().max_degree();
+
+    Table<dim, typename FESeriesType::CoefficientType> expansion_coefficients;
+    resize(expansion_coefficients, max_degree);
+
+    Vector<typename VectorType::value_type>                   local_dof_values;
+    std::vector<double>                                       ln_k;
+    std::pair<std::vector<unsigned int>, std::vector<double>> res;
+    for (auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          local_dof_values.reinit(cell->get_fe().dofs_per_cell);
+
+          auto solution_it              = all_solutions.cbegin();
+          auto smoothness_indicators_it = all_smoothness_indicators.begin();
+          for (; solution_it != all_solutions.cend();
+               ++solution_it, ++smoothness_indicators_it)
+            {
+              // Inside the loop, we first need to get the values of the local
+              // degrees of freedom and then need to compute the series
+              // expansion by multiplying this vector with the matrix ${\cal F}$
+              // corresponding to this finite element.
+              cell->get_dof_values(*(*solution_it), local_dof_values);
+
+              fe_series.calculate(local_dof_values,
+                                  cell->active_fe_index(),
+                                  expansion_coefficients);
+
+              // We fit our exponential decay of expansion coefficients to the
+              // provided regression_strategy on each possible value of |k|. To
+              // this end, we use FESeries::process_coefficients() to rework
+              // coefficients into the desired format.
+              res = FESeries::process_coefficients<dim>(
+                expansion_coefficients,
+                std::bind(&predicate<dim>, std::placeholders::_1, max_degree),
+                regression_strategy);
+
+              Assert(res.first.size() == res.second.size(), ExcInternalError());
+
+              // Prepare linear equation for the logarithmic least squares fit.
+              //
+              // First, calculate ln(|k|). This vector will be the same for all
+              // the cells so we can calculate ln(|k|) only once.
+              //
+              // For Fourier expansion, this translates to
+              // ln(2*pi*sqrt(predicate)) = ln(2*pi) + 0.5*ln(predicate). Since
+              // we are just interested in a linear regression later, we omit
+              // the ln(2*pi) factor.
+              // For Legendre expansion, this translates to
+              // 0.5*ln(predicate) as well, without the pi factor.
+              if (ln_k.empty())
+                {
+                  ln_k.resize(res.first.size());
+                  for (unsigned int f = 0; f < res.first.size(); ++f)
+                    ln_k[f] = 0.5 * std::log((double)res.first[f]);
+                }
+
+              // Second, calculate ln(U_k).
+              for (double &residual_element : res.second)
+                residual_element = std::log(residual_element);
+
+              // Last, do the linear regression.
+              std::pair<double, double> fit =
+                FESeries::linear_regression(ln_k, res.second);
+
+              // Compute the Sobolev index s=mu-dim/2 and store it in the vector
+              // of estimated values for each cell.
+              (*(*smoothness_indicators_it))(cell->active_cell_index()) =
+                (float)(-fit.first - .5 * dim);
+            }
+        }
+  }
+
+
+
+  template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+  void
+  estimate_by_coeff_decay(FESeriesType &              fe_series,
+                          const DoFHandlerType &      dof_handler,
+                          const VectorType &          solution,
+                          Vector<float> &             smoothness_indicators,
+                          const VectorTools::NormType regression_strategy)
+  {
+    const std::vector<const VectorType *> all_solutions(1, &solution);
+    const std::vector<Vector<float> *>    all_smoothness_indicators(
+      1, &smoothness_indicators);
+
+    estimate_by_coeff_decay(fe_series,
+                            dof_handler,
+                            all_solutions,
+                            all_smoothness_indicators,
+                            regression_strategy);
+  }
+
+
+
+  template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+  void
+  estimate_by_coeff_decay(
+    const DoFHandlerType &                 dof_handler,
+    const std::vector<const VectorType *> &all_solutions,
+    const std::vector<Vector<float> *> &   all_smoothness_indicators,
+    const VectorTools::NormType            regression_strategy)
+  {
+    const unsigned int dim = DoFHandlerType::dimension;
+    const unsigned int max_degree =
+      dof_handler.get_fe_collection().max_degree();
+
+    // We initialize a series expansion object object which will be used to
+    // calculate the expansion coefficients. In addition to the
+    // hp::FECollection, we need to provide quadrature rules hp::QCollection for
+    // integration on the reference cell.
+    // We will need to assemble the expansion matrices for each of the finite
+    // 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 2-point Gauss formula
+    // as many times as the maximal exponent we use for the term exp(ikx).
+    QGauss<1>      base_quadrature(2);
+    QIterated<dim> quadrature(base_quadrature, max_degree);
+
+    hp::QCollection<dim> expansion_q_collection;
+    for (unsigned int i = 0; i < dof_handler.get_fe_collection().size(); ++i)
+      expansion_q_collection.push_back(quadrature);
+
+    // The FESeries::Fourier class' constructor first parameter $N$ defines the
+    // number of coefficients in 1D with the total number of coefficients being
+    // $N^{dim}$.
+    FESeriesType fe_series(max_degree,
+                           dof_handler.get_fe_collection(),
+                           expansion_q_collection);
+
+    estimate_by_coeff_decay(fe_series,
+                            dof_handler,
+                            all_solutions,
+                            all_smoothness_indicators,
+                            regression_strategy);
+  }
+
+
+
+  template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+  void
+  estimate_by_coeff_decay(const DoFHandlerType &      dof_handler,
+                          const VectorType &          solution,
+                          Vector<float> &             smoothness_indicators,
+                          const VectorTools::NormType regression_strategy)
+  {
+    const std::vector<const VectorType *> all_solutions(1, &solution);
+    const std::vector<Vector<float> *>    all_smoothness_indicators(
+      1, &smoothness_indicators);
+
+    estimate_by_coeff_decay<FESeriesType>(dof_handler,
+                                          all_solutions,
+                                          all_smoothness_indicators,
+                                          regression_strategy);
+  }
+} // namespace SmoothnessEstimator
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 9a9494c928f31fcfe9ba91c5c565f09b73daefd8..028dab93e728f7e3da2605666480561b51f55f4e 100644 (file)
@@ -48,6 +48,7 @@ SET(_separate_src
   matrix_creator_inst2.cc
   matrix_creator_inst3.cc
   point_value_history.cc
+  smoothness_estimator.cc
   solution_transfer.cc
   solution_transfer_inst2.cc
   solution_transfer_inst3.cc
@@ -90,6 +91,7 @@ SET(_inst
   matrix_creator.inst.in
   matrix_tools.inst.in
   point_value_history.inst.in
+  smoothness_estimator.inst.in
   solution_transfer.inst.in
   time_dependent.inst.in
   vector_tools_boundary.inst.in
diff --git a/source/numerics/smoothness_estimator.cc b/source/numerics/smoothness_estimator.cc
new file mode 100644 (file)
index 0000000..efba58d
--- /dev/null
@@ -0,0 +1,38 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/hp/dof_handler.h>
+
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/la_vector.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/smoothness_estimator.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+#include "smoothness_estimator.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/numerics/smoothness_estimator.inst.in b/source/numerics/smoothness_estimator.inst.in
new file mode 100644 (file)
index 0000000..25db4d0
--- /dev/null
@@ -0,0 +1,55 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2010 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS;
+     VEC : REAL_VECTOR_TYPES;
+     DH : DOFHANDLER_TEMPLATES;
+     EXP : SERIES_EXPANSION_TEMPLATES)
+  {
+#if deal_II_dimension != 1 && deal_II_dimension <= deal_II_space_dimension
+    template void SmoothnessEstimator::estimate_by_coeff_decay<
+      EXP<deal_II_dimension, deal_II_space_dimension>,
+      DH<deal_II_dimension, deal_II_space_dimension>,
+      VEC>(EXP<deal_II_dimension, deal_II_space_dimension> &,
+           const DH<deal_II_dimension, deal_II_space_dimension> &,
+           const std::vector<const VEC *> &,
+           const std::vector<Vector<float> *> &,
+           const VectorTools::NormType);
+    template void SmoothnessEstimator::estimate_by_coeff_decay<
+      EXP<deal_II_dimension, deal_II_space_dimension>,
+      DH<deal_II_dimension, deal_II_space_dimension>,
+      VEC>(EXP<deal_II_dimension, deal_II_space_dimension> &,
+           const DH<deal_II_dimension, deal_II_space_dimension> &,
+           const VEC &,
+           Vector<float> &,
+           const VectorTools::NormType);
+
+    template void SmoothnessEstimator::estimate_by_coeff_decay<
+      EXP<deal_II_dimension, deal_II_space_dimension>,
+      DH<deal_II_dimension, deal_II_space_dimension>,
+      VEC>(const DH<deal_II_dimension, deal_II_space_dimension> &,
+           const std::vector<const VEC *> &,
+           const std::vector<Vector<float> *> &,
+           const VectorTools::NormType);
+    template void SmoothnessEstimator::estimate_by_coeff_decay<
+      EXP<deal_II_dimension, deal_II_space_dimension>,
+      DH<deal_II_dimension, deal_II_space_dimension>,
+      VEC>(const DH<deal_II_dimension, deal_II_space_dimension> &,
+           const VEC &,
+           Vector<float> &,
+           const VectorTools::NormType);
+#endif
+  }
index 2356e185094873d7053ddfedbb9aa937eef690f0..446008e61957bbf41350bf4a39be21692cc5f444 100644 (file)
@@ -29,7 +29,6 @@
 #include <deal.II/dofs/dof_tools.h>
 
 #include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_series.h>
 
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_refinement.h>
@@ -51,6 +50,7 @@
 #include <deal.II/numerics/data_out.h>
 #include <deal.II/numerics/error_estimator.h>
 #include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/smoothness_estimator.h>
 #include <deal.II/numerics/vector_tools.h>
 
 #include <complex>
@@ -84,8 +84,6 @@ namespace Step27
     void
     create_coarse_grid();
     void
-    estimate_smoothness(Vector<float> &smoothness_indicators);
-    void
     postprocess(const unsigned int cycle);
 
     Triangulation<dim> triangulation;
@@ -95,11 +93,6 @@ namespace Step27
     hp::QCollection<dim>     quadrature_collection;
     hp::QCollection<dim - 1> face_quadrature_collection;
 
-    hp::QCollection<dim>                    fourier_q_collection;
-    std::shared_ptr<FESeries::Fourier<dim>> fourier;
-    std::vector<double>                     ln_k;
-    Table<dim, std::complex<double>>        fourier_coefficients;
-
     AffineConstraints<double> constraints;
 
     SparsityPattern      sparsity_pattern;
@@ -138,13 +131,6 @@ namespace Step27
   }
 
 
-  template <typename T>
-  void resize(Table<2, T> &coeff, const unsigned int N)
-  {
-    coeff.reinit(N, N);
-  }
-
-
   template <int dim>
   LaplaceProblem<dim>::LaplaceProblem()
     : dof_handler(triangulation)
@@ -156,18 +142,6 @@ namespace Step27
         quadrature_collection.push_back(QGauss<dim>(degree + 1));
         face_quadrature_collection.push_back(QGauss<dim - 1>(degree + 1));
       }
-
-    const unsigned int N = max_degree;
-
-    QGauss<1>      base_quadrature(2);
-    QIterated<dim> quadrature(base_quadrature, N);
-    for (unsigned int i = 0; i < fe_collection.size(); i++)
-      fourier_q_collection.push_back(quadrature);
-
-    fourier = std::make_shared<FESeries::Fourier<dim>>(N,
-                                                       fe_collection,
-                                                       fourier_q_collection);
-    resize(fourier_coefficients, N);
   }
 
 
@@ -295,9 +269,9 @@ namespace Step27
       solution,
       estimated_error_per_cell);
 
-
-    Vector<float> smoothness_indicators(triangulation.n_active_cells());
-    estimate_smoothness(smoothness_indicators);
+    Vector<float> smoothness_indicators;
+    SmoothnessEstimator::estimate_by_coeff_decay<FESeries::Fourier<dim>>(
+      dof_handler, solution, smoothness_indicators);
 
     // Output to VTK
     if (false)
@@ -458,192 +432,6 @@ namespace Step27
         postprocess(cycle);
       }
   }
-
-  template <int dim>
-  std::pair<bool, unsigned int>
-  predicate_ind(const TableIndices<dim> &ind);
-
-  template <>
-  std::pair<bool, unsigned int>
-  predicate_ind<2>(const TableIndices<2> &ind)
-  {
-    const unsigned int v = ind[0] * ind[0] + ind[1] * ind[1];
-    if (v > 0 && v < 7 * 7)
-      return std::make_pair(true, v);
-    else
-      return std::make_pair(false, v);
-  }
-
-  template <int dim>
-  void
-  LaplaceProblem<dim>::estimate_smoothness(Vector<float> &smoothness_indicators)
-  {
-#ifdef OLD
-    const unsigned int N = max_degree;
-
-    std::vector<Tensor<1, dim>> k_vectors;
-    std::vector<unsigned int>   k_vectors_magnitude;
-    switch (dim)
-      {
-        case 2:
-          {
-            for (unsigned int i = 0; i < N; ++i)
-              for (unsigned int j = 0; j < N; ++j)
-                if (!((i == 0) && (j == 0)) && (i * i + j * j < N * N))
-                  {
-                    k_vectors.push_back(
-                      Point<dim>(2. * numbers::PI * i, 2. * numbers::PI * j));
-                    k_vectors_magnitude.push_back(i * i + j * j);
-                  }
-
-            break;
-          }
-
-        case 3:
-          {
-            for (unsigned int i = 0; i < N; ++i)
-              for (unsigned int j = 0; j < N; ++j)
-                for (unsigned int k = 0; k < N; ++k)
-                  if (!((i == 0) && (j == 0) && (k == 0)) &&
-                      (i * i + j * j + k * k < N * N))
-                    {
-                      k_vectors.push_back(Point<dim>(2. * numbers::PI * i,
-                                                     2. * numbers::PI * j,
-                                                     2. * numbers::PI * k));
-                      k_vectors_magnitude.push_back(i * i + j * j + k * k);
-                    }
-
-            break;
-          }
-
-        default:
-          Assert(false, ExcNotImplemented());
-      }
-
-    const unsigned      n_fourier_modes = k_vectors.size();
-    std::vector<double> ln_k(n_fourier_modes);
-    for (unsigned int i = 0; i < n_fourier_modes; ++i)
-      ln_k[i] = std::log(k_vectors[i].norm());
-
-
-    std::vector<Table<2, std::complex<double>>> fourier_transform_matrices(
-      fe_collection.size());
-
-    QGauss<1>      base_quadrature(2);
-    QIterated<dim> quadrature(base_quadrature, N);
-
-
-    for (unsigned int fe = 0; fe < fe_collection.size(); ++fe)
-      {
-        fourier_transform_matrices[fe].reinit(n_fourier_modes,
-                                              fe_collection[fe].dofs_per_cell);
-
-        for (unsigned int k = 0; k < n_fourier_modes; ++k)
-          for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
-            {
-              std::complex<double> sum = 0;
-              for (unsigned int q = 0; q < quadrature.size(); ++q)
-                {
-                  const Point<dim> x_q = quadrature.point(q);
-                  sum += std::exp(std::complex<double>(0, 1) *
-                                  (k_vectors[k] * x_q)) *
-                         fe_collection[fe].shape_value(j, x_q) *
-                         quadrature.weight(q);
-                }
-              fourier_transform_matrices[fe](k, j) = sum;
-            }
-      }
-
-    std::vector<std::complex<double>> fourier_coefficients(n_fourier_modes);
-    Vector<double>                    local_dof_values;
-
-    typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler
-                                                                .begin_active(),
-                                                       endc = dof_handler.end();
-    for (; cell != endc; ++cell)
-      {
-        local_dof_values.reinit(cell->get_fe().dofs_per_cell);
-        cell->get_dof_values(solution, local_dof_values);
-
-        for (unsigned int f = 0; f < n_fourier_modes; ++f)
-          {
-            fourier_coefficients[f] = 0;
-
-            for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i)
-              fourier_coefficients[f] +=
-                fourier_transform_matrices[cell->active_fe_index()](f, i) *
-                local_dof_values(i);
-          }
-
-        std::map<unsigned int, double> k_to_max_U_map;
-        for (unsigned int f = 0; f < n_fourier_modes; ++f)
-          if ((k_to_max_U_map.find(k_vectors_magnitude[f]) ==
-               k_to_max_U_map.end()) ||
-              (k_to_max_U_map[k_vectors_magnitude[f]] <
-               std::abs(fourier_coefficients[f])))
-            k_to_max_U_map[k_vectors_magnitude[f]] =
-              std::abs(fourier_coefficients[f]);
-
-        double sum_1 = 0, sum_ln_k = 0, sum_ln_k_square = 0, sum_ln_U = 0,
-               sum_ln_U_ln_k = 0;
-        for (unsigned int f = 0; f < n_fourier_modes; ++f)
-          if (k_to_max_U_map[k_vectors_magnitude[f]] ==
-              std::abs(fourier_coefficients[f]))
-            {
-              sum_1 += 1;
-              sum_ln_k += ln_k[f];
-              sum_ln_k_square += ln_k[f] * ln_k[f];
-              sum_ln_U += std::log(std::abs(fourier_coefficients[f]));
-              sum_ln_U_ln_k +=
-                std::log(std::abs(fourier_coefficients[f])) * ln_k[f];
-            }
-
-        const double mu =
-          (1. / (sum_1 * sum_ln_k_square - sum_ln_k * sum_ln_k) *
-           (sum_ln_k * sum_ln_U - sum_1 * sum_ln_U_ln_k));
-
-        smoothness_indicators(cell->active_cell_index()) = mu - 1. * dim / 2;
-      }
-#else
-    Vector<double> local_dof_values;
-
-    typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler
-                                                                .begin_active(),
-                                                       endc = dof_handler.end();
-    for (; cell != endc; ++cell)
-      {
-        local_dof_values.reinit(cell->get_fe().dofs_per_cell);
-        cell->get_dof_values(solution, local_dof_values);
-
-        fourier->calculate(local_dof_values,
-                           cell->active_fe_index(),
-                           fourier_coefficients);
-
-        std::pair<std::vector<unsigned int>, std::vector<double>> res =
-          FESeries::process_coefficients<dim>(fourier_coefficients,
-                                              predicate_ind<dim>,
-                                              VectorTools::Linfty_norm);
-
-        Assert(res.first.size() == res.second.size(), ExcInternalError());
-
-        if (ln_k.size() == 0)
-          {
-            ln_k.resize(res.first.size(), 0);
-            for (unsigned int f = 0; f < ln_k.size(); f++)
-              ln_k[f] =
-                std::log(2.0 * numbers::PI * std::sqrt(1. * res.first[f]));
-          }
-
-        for (unsigned int f = 0; f < res.second.size(); f++)
-          res.second[f] = std::log(res.second[f]);
-
-        std::pair<double, double> fit =
-          FESeries::linear_regression(ln_k, res.second);
-        smoothness_indicators(cell->active_cell_index()) =
-          -fit.first - 1. * dim / 2;
-      }
-#endif
-  }
 } // namespace Step27
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.