From: Marc Fehling Date: Wed, 21 Nov 2018 10:07:00 +0000 (+0100) Subject: SmoothnessEstimator: Initial implementation with FESeries template. X-Git-Tag: v9.2.0-rc1~12^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b649028d317e41361558033b0213179c179ce6f4;p=dealii.git SmoothnessEstimator: Initial implementation with FESeries template. --- diff --git a/cmake/config/template-arguments.in b/cmake/config/template-arguments.in index efe2f5111e..836fb7fc33 100644 --- a/cmake/config/template-arguments.in +++ b/cmake/config/template-arguments.in @@ -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 } diff --git a/examples/step-27/doc/intro.dox b/examples/step-27/doc/intro.dox index 0597da6789..e30a733ca0 100644 --- a/examples/step-27/doc/intro.dox +++ b/examples/step-27/doc/intro.dox @@ -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}$. -

What we have to do

+

What has to be done

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.

Compensating for anisotropy

diff --git a/examples/step-27/step-27.cc b/examples/step-27/step-27.cc index 54fbca5271..3d32f9638f 100644 --- a/examples/step-27/step-27.cc +++ b/examples/step-27/step-27.cc @@ -49,11 +49,11 @@ // These are the new files we need. The first and second provide hp // 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 #include -#include +#include // 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 &smoothness_indicators); void postprocess(const unsigned int cycle); - std::pair predicate(const TableIndices &indices); Triangulation triangulation; @@ -110,11 +108,6 @@ namespace Step27 hp::QCollection quadrature_collection; hp::QCollection face_quadrature_collection; - hp::QCollection fourier_q_collection; - std::shared_ptr> fourier; - std::vector ln_k; - Table> fourier_coefficients; - AffineConstraints 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 - void resize(Table &coeff, const unsigned int N) - { - TableIndices size; - for (unsigned int d = 0; d < dim; d++) - size[d] = N; - coeff.reinit(size); - } - template LaplaceProblem::LaplaceProblem() : dof_handler(triangulation) @@ -195,48 +171,6 @@ namespace Step27 quadrature_collection.push_back(QGauss(degree + 1)); face_quadrature_collection.push_back(QGauss(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 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>(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 estimated_error_per_cell(triangulation.n_active_cells()); KellyErrorEstimator::estimate( dof_handler, @@ -413,9 +346,11 @@ namespace Step27 solution, estimated_error_per_cell); - - Vector 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 smoothness_indicators; + SmoothnessEstimator::estimate_by_coeff_decay>( + 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 bool and unsigned int. 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 bool - // 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 - std::pair - LaplaceProblem::predicate(const TableIndices &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 - void - LaplaceProblem::estimate_smoothness(Vector &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 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 - // local_dof_values 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 local_dof_values, - // cell->active_fe_index() 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> res = - FESeries::process_coefficients( - fourier_coefficients, - [this](const TableIndices &indices) { - return this->predicate(indices); - }, - VectorTools::Linfty_norm); - - Assert(res.first.size() == res.second.size(), ExcInternalError()); - - // The first vector in the std::pair 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 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 diff --git a/include/deal.II/hp/fe_collection.h b/include/deal.II/hp/fe_collection.h index 1764bffd5d..746fea8972 100644 --- a/include/deal.II/hp/fe_collection.h +++ b/include/deal.II/hp/fe_collection.h @@ -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 + unsigned int + FECollection::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 unsigned int FECollection::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 index 0000000000..9965cd8444 --- /dev/null +++ b/include/deal.II/numerics/smoothness_estimator.h @@ -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 + +#include + +#include + +#include + + +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 + void + estimate_by_coeff_decay( + FESeriesType & fe_series, + const DoFHandlerType & dof_handler, + const std::vector &all_solutions, + const std::vector *> & all_smoothness_indicators, + const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm); + + /** + * Same as the function above, only for one @p solution vector. + */ + template + void + estimate_by_coeff_decay( + FESeriesType & fe_series, + const DoFHandlerType & dof_handler, + const VectorType & solution, + Vector & 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>( + * dof_handler, all_solutions, all_smoothness_indicators); + * @endcode + */ + template + void + estimate_by_coeff_decay( + const DoFHandlerType & dof_handler, + const std::vector &all_solutions, + const std::vector *> & all_smoothness_indicators, + const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm); + + /** + * Same as the function above, only for one @p solution vector. + */ + template + void + estimate_by_coeff_decay( + const DoFHandlerType & dof_handler, + const VectorType & solution, + Vector & 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 index 0000000000..2570ddcbc6 --- /dev/null +++ b/include/deal.II/numerics/smoothness_estimator.templates.h @@ -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 + +#include + +#include + +#include + +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace SmoothnessEstimator +{ + namespace + { + /** + * Resizes @p coeff to @p N in each dimension. + */ + template + void + resize(Table &coeff, const unsigned int N) + { + TableIndices 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 + std::pair + predicate(const TableIndices &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 + void + estimate_by_coeff_decay( + FESeriesType & fe_series, + const DoFHandlerType & dof_handler, + const std::vector &all_solutions, + const std::vector *> & 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 expansion_coefficients; + resize(expansion_coefficients, max_degree); + + Vector local_dof_values; + std::vector ln_k; + std::pair, std::vector> 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( + expansion_coefficients, + std::bind(&predicate, 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 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 + void + estimate_by_coeff_decay(FESeriesType & fe_series, + const DoFHandlerType & dof_handler, + const VectorType & solution, + Vector & smoothness_indicators, + const VectorTools::NormType regression_strategy) + { + const std::vector all_solutions(1, &solution); + const std::vector *> all_smoothness_indicators( + 1, &smoothness_indicators); + + estimate_by_coeff_decay(fe_series, + dof_handler, + all_solutions, + all_smoothness_indicators, + regression_strategy); + } + + + + template + void + estimate_by_coeff_decay( + const DoFHandlerType & dof_handler, + const std::vector &all_solutions, + const std::vector *> & 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 quadrature(base_quadrature, max_degree); + + hp::QCollection 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 + void + estimate_by_coeff_decay(const DoFHandlerType & dof_handler, + const VectorType & solution, + Vector & smoothness_indicators, + const VectorTools::NormType regression_strategy) + { + const std::vector all_solutions(1, &solution); + const std::vector *> all_smoothness_indicators( + 1, &smoothness_indicators); + + estimate_by_coeff_decay(dof_handler, + all_solutions, + all_smoothness_indicators, + regression_strategy); + } +} // namespace SmoothnessEstimator + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/numerics/CMakeLists.txt b/source/numerics/CMakeLists.txt index 9a9494c928..028dab93e7 100644 --- a/source/numerics/CMakeLists.txt +++ b/source/numerics/CMakeLists.txt @@ -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 index 0000000000..efba58d335 --- /dev/null +++ b/source/numerics/smoothness_estimator.cc @@ -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 +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +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 index 0000000000..25db4d0a05 --- /dev/null +++ b/source/numerics/smoothness_estimator.inst.in @@ -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, + DH, + VEC>(EXP &, + const DH &, + const std::vector &, + const std::vector *> &, + const VectorTools::NormType); + template void SmoothnessEstimator::estimate_by_coeff_decay< + EXP, + DH, + VEC>(EXP &, + const DH &, + const VEC &, + Vector &, + const VectorTools::NormType); + + template void SmoothnessEstimator::estimate_by_coeff_decay< + EXP, + DH, + VEC>(const DH &, + const std::vector &, + const std::vector *> &, + const VectorTools::NormType); + template void SmoothnessEstimator::estimate_by_coeff_decay< + EXP, + DH, + VEC>(const DH &, + const VEC &, + Vector &, + const VectorTools::NormType); +#endif + } diff --git a/tests/hp/step-27.cc b/tests/hp/step-27.cc index 2356e18509..446008e619 100644 --- a/tests/hp/step-27.cc +++ b/tests/hp/step-27.cc @@ -29,7 +29,6 @@ #include #include -#include #include #include @@ -51,6 +50,7 @@ #include #include #include +#include #include #include @@ -84,8 +84,6 @@ namespace Step27 void create_coarse_grid(); void - estimate_smoothness(Vector &smoothness_indicators); - void postprocess(const unsigned int cycle); Triangulation triangulation; @@ -95,11 +93,6 @@ namespace Step27 hp::QCollection quadrature_collection; hp::QCollection face_quadrature_collection; - hp::QCollection fourier_q_collection; - std::shared_ptr> fourier; - std::vector ln_k; - Table> fourier_coefficients; - AffineConstraints constraints; SparsityPattern sparsity_pattern; @@ -138,13 +131,6 @@ namespace Step27 } - template - void resize(Table<2, T> &coeff, const unsigned int N) - { - coeff.reinit(N, N); - } - - template LaplaceProblem::LaplaceProblem() : dof_handler(triangulation) @@ -156,18 +142,6 @@ namespace Step27 quadrature_collection.push_back(QGauss(degree + 1)); face_quadrature_collection.push_back(QGauss(degree + 1)); } - - const unsigned int N = max_degree; - - QGauss<1> base_quadrature(2); - QIterated quadrature(base_quadrature, N); - for (unsigned int i = 0; i < fe_collection.size(); i++) - fourier_q_collection.push_back(quadrature); - - fourier = std::make_shared>(N, - fe_collection, - fourier_q_collection); - resize(fourier_coefficients, N); } @@ -295,9 +269,9 @@ namespace Step27 solution, estimated_error_per_cell); - - Vector smoothness_indicators(triangulation.n_active_cells()); - estimate_smoothness(smoothness_indicators); + Vector smoothness_indicators; + SmoothnessEstimator::estimate_by_coeff_decay>( + dof_handler, solution, smoothness_indicators); // Output to VTK if (false) @@ -458,192 +432,6 @@ namespace Step27 postprocess(cycle); } } - - template - std::pair - predicate_ind(const TableIndices &ind); - - template <> - std::pair - 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 - void - LaplaceProblem::estimate_smoothness(Vector &smoothness_indicators) - { -#ifdef OLD - const unsigned int N = max_degree; - - std::vector> k_vectors; - std::vector 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(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(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 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>> fourier_transform_matrices( - fe_collection.size()); - - QGauss<1> base_quadrature(2); - QIterated 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 sum = 0; - for (unsigned int q = 0; q < quadrature.size(); ++q) - { - const Point x_q = quadrature.point(q); - sum += std::exp(std::complex(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> fourier_coefficients(n_fourier_modes); - Vector local_dof_values; - - typename hp::DoFHandler::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 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 local_dof_values; - - typename hp::DoFHandler::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> res = - FESeries::process_coefficients(fourier_coefficients, - predicate_ind, - 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 fit = - FESeries::linear_regression(ln_k, res.second); - smoothness_indicators(cell->active_cell_index()) = - -fit.first - 1. * dim / 2; - } -#endif - } } // namespace Step27