@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 }
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
@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
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>
// 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.
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;
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;
// 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)
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);
}
// 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,
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
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
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.
+ 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
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
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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
+ }
#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>
#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>
void
create_coarse_grid();
void
- estimate_smoothness(Vector<float> &smoothness_indicators);
- void
postprocess(const unsigned int cycle);
Triangulation<dim> triangulation;
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;
}
- 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)
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);
}
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)
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
projects / dealii.git / commitdiff