From 5ac46f9897eabf99cbb9a3f168e7fd3dfb3c137e Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Fri, 13 May 2016 22:14:27 +0200 Subject: [PATCH] introduce FESeries namespace and use it in step-27 The namespace provides two classes to perform expansion into orthogonal series: Fourier and Legendre. Additionally auxiliary functions to process expansion coefficients and do linear regression are implemented. It appears that there was a bug in step-27, nammely k-vectors should have a 2\pi multiplier instead of \pi. This was fixed and step-27 was updated. --- doc/news/changes.h | 11 + examples/step-27/doc/intro.dox | 43 +- examples/step-27/doc/results.dox | 30 +- examples/step-27/step-27.cc | 332 +++++------ include/deal.II/fe/fe_series.h | 422 ++++++++++++++ source/fe/CMakeLists.txt | 1 + source/fe/fe_series.cc | 423 ++++++++++++++ tests/fe/fe_series_01.cc | 98 ++++ tests/fe/fe_series_01.output | 11 + tests/fe/fe_series_02.cc | 207 +++++++ tests/fe/fe_series_02.with_gsl=on.output | 12 + tests/fe/fe_series_03.cc | 66 +++ tests/fe/fe_series_03.output | 33 ++ tests/fe/fe_series_04.cc | 203 +++++++ tests/fe/fe_series_04.with_gsl=on.output | 27 + tests/fe/fe_series_05.cc | 268 +++++++++ tests/fe/fe_series_05.with_gsl=on.output | 29 + tests/fe/fe_series_06.cc | 78 +++ tests/fe/fe_series_06.output | 43 ++ tests/fe/fe_series_07.cc | 70 +++ tests/fe/fe_series_07.output | 21 + tests/fe/fe_series_08.cc | 60 ++ tests/fe/fe_series_08.output | 3 + tests/hp/step-27.cc | 707 +++++++++++++++++++++++ tests/hp/step-27.output | 24 + 25 files changed, 3003 insertions(+), 219 deletions(-) create mode 100644 include/deal.II/fe/fe_series.h create mode 100644 source/fe/fe_series.cc create mode 100644 tests/fe/fe_series_01.cc create mode 100644 tests/fe/fe_series_01.output create mode 100644 tests/fe/fe_series_02.cc create mode 100644 tests/fe/fe_series_02.with_gsl=on.output create mode 100644 tests/fe/fe_series_03.cc create mode 100644 tests/fe/fe_series_03.output create mode 100644 tests/fe/fe_series_04.cc create mode 100644 tests/fe/fe_series_04.with_gsl=on.output create mode 100644 tests/fe/fe_series_05.cc create mode 100644 tests/fe/fe_series_05.with_gsl=on.output create mode 100644 tests/fe/fe_series_06.cc create mode 100644 tests/fe/fe_series_06.output create mode 100644 tests/fe/fe_series_07.cc create mode 100644 tests/fe/fe_series_07.output create mode 100644 tests/fe/fe_series_08.cc create mode 100644 tests/fe/fe_series_08.output create mode 100644 tests/hp/step-27.cc create mode 100644 tests/hp/step-27.output diff --git a/doc/news/changes.h b/doc/news/changes.h index 38d4b32b26..cef4132d1a 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -139,6 +139,17 @@ inconvenience this causes.

General

    + +
  1. New: Add new classes to expand a scalar finite element solution into + the orthogonal bases FESeries::Fourier and FESeries::Legendre. Also + provide auxiliary functions to calculate norms of subsets of expansion + coefficients FESeries::process_coefficients and linear regression + FESeries::linear_regression. Update step-27 to use this namespace to drive + the hp-adaptive FEM solution process. +
    + (Denis Davydov, 2016/06/23) +
    2. +
  2. New: The tutorial step-55 shows how to solve the Stokes system in parallel with PETSc or Trilinos.
    diff --git a/examples/step-27/doc/intro.dox b/examples/step-27/doc/intro.dox index 3945e32e73..5172b791d4 100644 --- a/examples/step-27/doc/intro.dox +++ b/examples/step-27/doc/intro.dox @@ -312,27 +312,30 @@ Assuming that the cell $K$ is not degenerate, i.e. that the mapping from the unit cell to cell $K$ is sufficiently regular, above condition is of course equivalent to @f[ - \int_{\hat K} |\nabla^s \hat u(\hat{\bf x})|^2 \; d\hat{\bf x} < \infty + \int_{\hat K} |\nabla^s \hat u(\hat{\bf x})|^2 \; d\hat{\bf x} < \infty\,, @f] where $\hat u(\hat{\bf x})$ is the function $u({\bf x})$ mapped back onto the unit cell $\hat K$. From here, we can do the following: first, let us define the Fourier series of $\hat u$ as @f[ - \hat U_{\bf k} - = \frac 1{(2\pi)^{d/2}} \int_{\hat K} e^{i {\bf k}\cdot \hat{\bf x}} \hat u(\hat{\bf x}) d\hat{\bf x} + \hat u(\hat{\bf x}) + = \sum_{\bf k} \hat U_{\bf k}\,e^{-i {\bf k}\cdot \hat{\bf x}}, @f] with Fourier vectors ${\bf k}=(k_x,k_y)$ in 2d, ${\bf k}=(k_x,k_y,k_z)$ -in 3d, etc, and $k_x,k_y,k_z=0,\pi,2\pi,3\pi,\ldots$. If we re-compose $\hat u$ -from $\hat U$ using the formula +in 3d, etc, and $k_x,k_y,k_z=0,2\pi,4\pi,\ldots$. The coefficients of expansion +$\hat U_{\bf k}$ can be obtained using $L^2$-orthogonality of the exponential basis @f[ - \hat u(\hat{\bf x}) - = \frac 1{(2\pi)^{d/2}} \sum_{\bf k} e^{-i {\bf k}\cdot \hat{\bf x}} \hat U_{\bf k}, +\int_{\hat K} e^{-i {\bf m}\cdot \hat{\bf x}} e^{i {\bf n}\cdot \hat{\bf x}} d\hat{\bf x} = \delta_{\bf m \bf n}, +@f] +that leads to the following expression +@f[ + \hat U_{\bf k} + = \int_{\hat K} e^{i {\bf k}\cdot \hat{\bf x}} \hat u(\hat{\bf x}) d\hat{\bf x} \,. @f] -then it becomes clear that we can write the $H^s$ norm of $\hat u$ as +It becomes clear that we can then write the $H^s$ norm of $\hat u$ as @f[ \int_{\hat K} |\nabla^s \hat u(\hat{\bf x})|^2 \; d\hat{\bf x} = - \frac 1{(2\pi)^d} \int_{\hat K} \left| \sum_{\bf k} |{\bf k}|^s e^{-i{\bf k}\cdot \hat{\bf x}} \hat U_{\bf k} @@ -347,12 +350,12 @@ In other words, if this norm is to be finite (i.e. for $\hat u(\hat{\bf x})$ to |\hat U_{\bf k}| = {\cal O}\left(|{\bf k}|^{-\left(s+1/2+\frac{d-1}{2}+\epsilon\right)}\right). @f] Put differently: the higher regularity $s$ we want, the faster the -Fourier coefficients have to go to zero. (If you wonder where the +Fourier coefficients have to go to zero. If you wonder where the additional exponent $\frac{d-1}2$ comes from: we would like to make use of the fact that $\sum_l a_l < \infty$ if the sequence $a_l = {\cal O}(l^{-1-\epsilon})$ for any $\epsilon>0$. The problem is that we here have a summation not only over a single variable, but over all -the integer multiples of $\pi$ that are located inside the +the integer multiples of $2\pi$ that are located inside the $d$-dimensional sphere, because we have vector components $k_x, k_y, \ldots$. In the same way as we prove that the sequence $a_l$ above converges by replacing the sum by an integral over the entire line, we @@ -363,7 +366,7 @@ much the same way as we can transform the volume element $dx\;dy$ into $2\pi r\; dr$. Consequently, it is no longer $|{\bf k}|^{2s}|\hat U_{\bf k}|^2$ that has to decay as ${\cal O}(|{\bf k}|^{-1-\epsilon})$, but it is in fact $|{\bf k}|^{2s}|\hat U_{\bf k}|^2 |{\bf k}|^{d-1}$. A -comparison of exponents yields the result.) +comparison of exponents yields the result. We can turn this around: Assume we are given a function $\hat u$ of unknown smoothness. Let us compute its Fourier coefficients $\hat U_{\bf k}$ @@ -377,10 +380,10 @@ then consequently the function we had here was in $H^{\mu-d/2}$.

    What we have to do

    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 series +a cell $K$? Clearly, the first step is to compute the Fourier coefficients of our solution. Fourier series being infinite series, we simplify our task by only computing the first few terms of the series, such that -$|{\bf k}|\le N$ with a cut-off $N$. (Let us parenthetically remark +$|{\bf k}|\le 2\pi N$ with a cut-off $N$. Let us parenthetically remark that we want to choose $N$ large enough so that we capture at least the variation of those shape functions that vary the most. On the other hand, we should not choose $N$ too large: clearly, a finite @@ -388,16 +391,16 @@ element function, being a polynomial, is in $C^\infty$ on any given cell, so the coefficients will have to decay exponentially at one point; since we want to estimate the smoothness of the function this polynomial approximates, not of the polynomial itself, we need to -choose a reasonable cutoff for $N$.) Either way, computing this series +choose a reasonable cutoff for $N$. Either way, computing this series is not particularly hard: from the definition @f[ \hat U_{\bf k} - = \frac 1{(2\pi)^{d/2}} \int_{\hat K} e^{i {\bf k}\cdot \hat{\bf x}} \hat u(\hat{\bf x}) d\hat{\bf x} + = \int_{\hat K} e^{i {\bf k}\cdot \hat{\bf x}} \hat u(\hat{\bf x}) d\hat{\bf x} @f] we see that we can compute the coefficient $\hat U_{\bf k}$ as @f[ \hat U_{\bf k} - = \frac 1{(2\pi)^{d/2}} + = \sum_{i=0}^{\textrm{\tiny dofs per cell}} \left[\int_{\hat K} e^{i {\bf k}\cdot \hat{\bf x}} \hat \varphi_i(\hat{\bf x}) d\hat{\bf x} \right] u_i, @@ -411,12 +414,14 @@ cell. In other words, we can write it as a matrix-vector product with the matrix @f[ {\cal F}_{{\bf k},j} - = \frac 1{(2\pi)^{d/2}} + = \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 coefficients $\hat U_{\bf k}$ is a rather trivial job. +To simplify our life even further, we will use FESeries::Fourier class which +does exactly this. The next task is that we have to estimate how fast these coefficients decay with $|{\bf k}|$. The problem is that, of course, we have only @@ -499,6 +504,8 @@ and \right]. @f] +This is nothing else but 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 diff --git a/examples/step-27/doc/results.dox b/examples/step-27/doc/results.dox index 6f6c745a22..21d2f5c70d 100644 --- a/examples/step-27/doc/results.dox +++ b/examples/step-27/doc/results.dox @@ -15,25 +15,25 @@ Cycle 0: Number of degrees of freedom: 3264 Number of constraints : 384 Cycle 1: - Number of active cells: 996 - Number of degrees of freedom: 5327 - Number of constraints : 962 + Number of active cells: 966 + Number of degrees of freedom: 5245 + Number of constraints : 936 Cycle 2: - Number of active cells: 1335 - Number of degrees of freedom: 8947 - Number of constraints : 2056 + Number of active cells: 1143 + Number of degrees of freedom: 8441 + Number of constraints : 1929 Cycle 3: - Number of active cells: 1626 - Number of degrees of freedom: 12994 - Number of constraints : 3355 + Number of active cells: 1356 + Number of degrees of freedom: 12349 + Number of constraints : 3046 Cycle 4: - Number of active cells: 1911 - Number of degrees of freedom: 17988 - Number of constraints : 4860 + Number of active cells: 1644 + Number of degrees of freedom: 18178 + Number of constraints : 4713 Cycle 5: - Number of active cells: 2577 - Number of degrees of freedom: 26936 - Number of constraints : 7074 + Number of active cells: 1728 + Number of degrees of freedom: 22591 + Number of constraints : 6095 @endcode The first thing we learn from this is that the number of constrained degrees diff --git a/examples/step-27/step-27.cc b/examples/step-27/step-27.cc index 0adfe4283a..33a509f034 100644 --- a/examples/step-27/step-27.cc +++ b/examples/step-27/step-27.cc @@ -14,7 +14,8 @@ * --------------------------------------------------------------------- * - * Author: Wolfgang Bangerth, Texas A&M University, 2006, 2007 + * Authors: Wolfgang Bangerth, Texas A&M University, 2006, 2007; + * Denis Davydov, University of Erlangen-Nuremberg, 2016. */ @@ -47,12 +48,14 @@ // These are the new files we need. The first one provides an alternative to // the usual SparsityPattern class and the DynamicSparsityPattern class -// already discussed in step-11 and step-18. The last two provide hp +// already discussed in step-11 and step-18. The second and third provide hp // versions of the DoFHandler and FEValues classes as described in the -// introduction of this program. +// introduction of this program. The last one provides Fourier transformation +// class on the unit cell. #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. @@ -98,8 +101,9 @@ namespace Step27 void assemble_system (); void solve (); void create_coarse_grid (); - void estimate_smoothness (Vector &smoothness_indicators) const; + void estimate_smoothness (Vector &smoothness_indicators); void postprocess (const unsigned int cycle); + std::pair predicate(const TableIndices &indices); Triangulation triangulation; @@ -108,6 +112,11 @@ namespace Step27 hp::QCollection quadrature_collection; hp::QCollection face_quadrature_collection; + hp::QCollection fourier_q_collection; + std_cxx11::shared_ptr> fourier; + std::vector ln_k; + Table > fourier_coefficients; + ConstraintMatrix constraints; SparsityPattern sparsity_pattern; @@ -164,6 +173,23 @@ 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 LaplaceProblem::LaplaceProblem () : @@ -176,6 +202,48 @@ 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_cxx11::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); } @@ -612,6 +680,32 @@ namespace Step27 // @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 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 @@ -619,134 +713,18 @@ namespace Step27 template void LaplaceProblem:: - estimate_smoothness (Vector &smoothness_indicators) const + estimate_smoothness (Vector &smoothness_indicators) { - // The first thing we need to do is to define the Fourier vectors ${\bf - // k}$ for which we want to compute Fourier coefficients of the solution - // on each cell. In 2d, we pick those vectors ${\bf k}=(\pi i, \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 3d case is handled analogously. 1d and dimensions higher - // than 3 are not implemented, and we guard our implementation by making - // sure that we receive an exception in case someone tries to compile the - // program for any of these dimensions. - // - // We exclude ${\bf k}=0$ to avoid problems computing $|{\bf k}|^{-mu}$ - // and $\ln |{\bf k}|$. The other vectors are stored in the field - // k_vectors. In addition, we store the square of the - // magnitude of each of these vectors (up to a factor $\pi^2$) in the - // k_vectors_magnitude array -- we will need that when we - // attempt to find out which of those Fourier coefficients corresponding - // to Fourier vectors of the same magnitude is the largest: - 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(numbers::PI * i, - numbers::PI * j)); - k_vectors_magnitude.push_back (i*i+j*j); - } - - break; - } - - case 3: - { - for (unsigned int i=0; i(numbers::PI * i, - numbers::PI * j, - numbers::PI * k)); - k_vectors_magnitude.push_back (i*i+j*j+k*k); - } - - break; - } - - default: - Assert (false, ExcNotImplemented()); - } - - // After we have set up the Fourier vectors, we also store their total - // number for simplicity, and compute the logarithm of the magnitude of - // each of these vectors since we will need it many times over further - // down below: - const unsigned n_fourier_modes = k_vectors.size(); - std::vector ln_k (n_fourier_modes); - for (unsigned int i=0; i > > - fourier_transform_matrices (fe_collection.size()); - - // In order to compute them, we of course can't perform the Fourier - // transform analytically, but have to approximate it using quadrature. To - // this end, we use a quadrature formula 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); - - // With this, we then loop over all finite elements in use, reinitialize - // the respective matrix ${\cal F}$ to the right size, and integrate each - // entry of the matrix numerically as ${\cal F}_{{\bf k},j}=\sum_q - // e^{i{\bf k}\cdot {\bf x}}\varphi_j({\bf x}_q) w_q$, where $x_q$ are the - // quadrature points and $w_q$ are the quadrature weights. Note that the - // imaginary unit $i=\sqrt{-1}$ is obtained from the standard C++ classes - // using std::complex@(0,1). - - // Because we work on the unit cell, we can do all this work without a - // mapping from reference to real cell and consequently do not need the - // FEValues class. - for (unsigned int fe=0; fe sum = 0; - for (unsigned int q=0; q 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::pow(2*numbers::PI, 1.*dim/2); - } - } - - // The next thing is to loop over all cells and do our work there, i.e. to + // 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 two arrays as scratch arrays in the loop and - // allocate them here to avoid repeated memory allocations: - std::vector > fourier_coefficients (n_fourier_modes); + // 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: @@ -760,76 +738,58 @@ namespace Step27 // 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. We need to write out the multiplication by hand because - // the objects holding the data do not have vmult-like - // functions declared: + // 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); - for (unsigned int f=0; fget_fe().dofs_per_cell; ++i) - fourier_coefficients[f] += - fourier_transform_matrices[cell->active_fe_index()](f,i) - * - local_dof_values(i); - } + 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 create a map that for each magnitude $|{\bf k}|$ - // stores the largest $|\hat U_{{\bf k}}|$ found so far, i.e. we - // overwrite the existing value (or add it to the map) if no value for - // the current $|{\bf k}|$ exists yet, or if the current value is - // larger than the previously stored one: - std::map k_to_max_U_map; - for (unsigned int f=0; f, std::vector > res = + FESeries::process_coefficients(fourier_coefficients, + std_cxx11::bind(&LaplaceProblem::predicate, + this, + std_cxx11::_1), + 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 linear regression fit to + // obtain $\mu$. + for (unsigned int f = 0; f < res.second.size(); f++) + res.second[f] = std::log(res.second[f]); - // With these so-computed sums, we can now evaluate the formula for - // $\mu$ derived in the introduction: - 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)); + 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()) = mu - 1.*dim/2; + smoothness_indicators(cell->active_cell_index()) = -fit.first - 1.*dim/2; } } } diff --git a/include/deal.II/fe/fe_series.h b/include/deal.II/fe/fe_series.h new file mode 100644 index 0000000000..fe4c461c3c --- /dev/null +++ b/include/deal.II/fe/fe_series.h @@ -0,0 +1,422 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__fe_series_H +#define dealii__fe_series_H + + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +/*!@addtogroup feall */ +/*@{*/ + + +/** + * This namespace offers functions to calculate expansion series of the + * solution on the reference element. Coefficients of expansion are often used + * to estimate local smoothness of the underlying FiniteElement field to decide + * on h- or p-adaptive refinement strategy. + * + * @author Denis Davydov, 2016; + */ +namespace FESeries +{ + /** + * A class to calculate expansion of a scalar FE field into Fourier series + * on a reference element. The exponential form of the Fourier series is + * based on completeness and Hermitian orthogonality of the set of exponential + * functions \f$ \phi_{\bf k}({\bf x}) = \exp(2 \pi i\, {\bf k} \cdot {\bf x})\f$. + * For example in 1D the L2-orthogonality condition reads + * @f[ + * \int_0^1 \phi_k(x) \phi_l^\ast(x) dx=\delta_{kl}. + * @f] + * Note that \f$ \phi_{\bf k} = \phi_{-\bf k}^\ast \f$. + * + * The arbitrary scalar FE field on the reference element can be expanded in + * the complete orthogonal exponential basis as + * @f[ + * u({\bf x}) + * = \sum_{\bf k} c_{\bf k} \phi_{\bf k}({\bf x}). + * @f] + * From the orthogonality property of the basis, it follows that + * @f[ + * c_{\bf k} = + * \int_{[0,1]^d} u({\bf x}) \phi_{\bf k}^\ast ({\bf x}) d{\bf x}\,. + * @f] + * It is this complex-valued expansion coefficients, that are calculated by + * this class. Note that \f$ u({\bf x}) = \sum_i u_i N_i({\bf x})\f$, + * where \f$ N_i({\bf x}) \f$ are real-valued FiniteElement shape functions. + * Consequently \f$ c_{\bf k} \equiv c_{-\bf k}^\ast \f$ and + * we only need to compute \f$ c_{\bf k} \f$ for positive indices + * \f$ \bf k \f$ . + * + * @author Denis Davydov, 2016. + */ + template + class Fourier : public Subscriptor + { + public: + /** + * A non-default constructor. The @p size_in_each_direction defines the number + * of modes in each direction, @p fe_collection is the hp::FECollection + * for which expansion will be used and @p q_collection is the hp::QCollection + * used to integrate the expansion for each FiniteElement + * in @p fe_collection. + */ + Fourier(const unsigned int size_in_each_direction, + const hp::FECollection &fe_collection, + const hp::QCollection &q_collection); + + /** + * Calculate @p fourier_coefficients of the cell vector field given by + * @p local_dof_values corresponding to FiniteElement with + * @p cell_active_fe_index . + */ + void calculate(const dealii::Vector &local_dof_values, + const unsigned int cell_active_fe_index, + Table > &fourier_coefficients); + + private: + /** + * hp::FECollection for which transformation matrices will be calculated. + */ + SmartPointer > fe_collection; + + /** + * hp::QCollection used in calculation of transformation matrices. + */ + SmartPointer > q_collection; + + /** + * Ensure that the transformation matrix for FiniteElement index + * @p fe_index is calculated. If not, calculate it. + */ + void ensure_existence(const unsigned int fe_index); + + /** + * Angular frequencies \f$ 2 \pi {\bf k} \f$ . + */ + Table > k_vectors; + + /** + * Transformation matrices for each FiniteElement. + */ + std::vector > > fourier_transform_matrices; + + /** + * Auxiliary vector to store unrolled coefficients. + */ + std::vector > unrolled_coefficients; + + }; + + /** + * A class to calculate expansion of a scalar FE field into series of Legendre + * functions on a reference element. + * + * Legendre functions are solutions to Legendre's differential equation + * @f[ + * \frac{d}{dx}\left([1-x^2] \frac{d}{dx} P_n(x)\right) + + * n[n+1] P_n(x) = 0 + * @f] + * and can be expressed using Rodrigues' formula + * @f[ + * P_n(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n}[x^2-1]^n. + * @f] + * These polynomials are orthogonal with respect to the \f$ L^2 \f$ inner + * product on the interval \f$ [-1;1] \f$ + * @f[ + * \int_{-1}^1 P_m(x) P_n(x) = \frac{2}{2n + 1} \delta_{mn} + * @f] + * and are complete. + * A family of \f$ L^2 \f$-orthogonal polynomials on \f$ [0;1] \f$ can be + * constructed via + * @f[ + * \widetilde P_m = \sqrt{2} P_m(2x-1). + * @f] + * + * + * An arbitrary scalar FE field on the reference element \f$ [0;1] \f$ can be + * expanded in the complete orthogonal basis as + * @f[ + * u(x) + * = \sum_{m} c_m \widetilde P_{m}(x). + * @f] + * From the orthogonality property of the basis, it follows that + * @f[ + * c_m = \frac{2m+1}{2} + * \int_0^1 u(x) \widetilde P_m(x) dx . + * @f] + * This class calculates coefficients \f$ c_{\bf k} \f$ using + * \f$ dim \f$-dimensional Legendre polynomials constructed from + * \f$ \widetilde P_m(x) \f$ using tensor product rule. + * + * @author Denis Davydov, 2016. + */ + template + class Legendre : public Subscriptor + { + public: + /** + * A non-default constructor. The @p size_in_each_direction defines the number + * of coefficients in each direction, @p fe_collection is the hp::FECollection + * for which expansion will be used and @p q_collection is the hp::QCollection + * used to integrate the expansion for each FiniteElement + * in @p fe_collection. + */ + Legendre(const unsigned int size_in_each_direction, + const hp::FECollection &fe_collection, + const hp::QCollection &q_collection); + + /** + * Calculate @p legendre_coefficients of the cell vector field given by + * @p local_dof_values corresponding to FiniteElement with + * @p cell_active_fe_index . + */ + void calculate(const dealii::Vector &local_dof_values, + const unsigned int cell_active_fe_index, + Table &legendre_coefficients); + + private: + /** + * Number of coefficients in each direction + */ + const unsigned int N; + + /** + * hp::FECollection for which transformation matrices will be calculated. + */ + SmartPointer > fe_collection; + + /** + * hp::QCollection used in calculation of transformation matrices. + */ + SmartPointer > q_collection; + + /** + * Ensure that the transformation matrix for FiniteElement index + * @p fe_index is calculated. If not, calculate it. + */ + void ensure_existence(const unsigned int fe_index); + + /** + * Transformation matrices for each FiniteElement. + */ + std::vector > legendre_transform_matrices; + + /** + * Auxiliary vector to store unrolled coefficients. + */ + std::vector unrolled_coefficients; + + }; + + + /** + * Calculate the @p norm of subsets of @p coefficients defined by + * @p predicate being constant. Return the pair of vectors of predicate values + * and the vector of calculated subset norms. + * + * @p predicate should return a pair of bool and unsigned int. + * The former is a flag whether a given TableIndices should be used in + * calculation, whereas the latter is the unrolled value of indices according + * to which the subsets of coefficients will be formed. + * + * @note Only the following values of @p norm are implemented and make sense + * in this case: mean, L1_norm, L2_norm, Linfty_norm. The mean norm can only + * be applied to real valued coefficients. + */ + template + std::pair,std::vector > + process_coefficients(const Table &coefficients, + const std_cxx11::function(const TableIndices &)> &predicate, + const VectorTools::NormType norm); + + + + /** + * Linear regression least-square fit of $y = k \, x + b$. + * The size of the input vectors should be equal and more than 1. + * The returned pair will contain $k$ (first) and $b$ (second). + */ + std::pair linear_regression(const std::vector &x, + const std::vector &y); + +} + +/*@}*/ + +#ifndef DOXYGEN + +// ------------------- inline and template functions ---------------- + +namespace +{ + template + void fill_map_index(const Table &coefficients, + const TableIndices &ind, + const std_cxx11::function(const TableIndices &)> &predicate, + std::map > &pred_to_values) + { + const std::pair pred_pair = predicate(ind); + // don't add a value if predicate is false + if (pred_pair.first == false) + return; + + const unsigned int &pred_value = pred_pair.second; + const T &coeff_value = coefficients(ind); + // If pred_value is not in the pred_to_values map, the element will be created. + // Otherwise a reference to the existing element is returned. + pred_to_values[pred_value].push_back(coeff_value); + } + + template + void fill_map(const Table<1,T> &coefficients, + const std_cxx11::function(const TableIndices<1> &)> &predicate, + std::map > &pred_to_values) + { + for (unsigned int i = 0; i < coefficients.size(0); i++) + { + const TableIndices<1> ind(i); + fill_map_index(coefficients,ind,predicate,pred_to_values); + } + + } + + template + void fill_map(const Table<2,T> &coefficients, + const std_cxx11::function(const TableIndices<2> &)> &predicate, + std::map > &pred_to_values) + { + for (unsigned int i = 0; i < coefficients.size(0); i++) + for (unsigned int j = 0; j < coefficients.size(1); j++) + { + const TableIndices<2> ind(i,j); + fill_map_index(coefficients,ind,predicate,pred_to_values); + } + + } + + template + void fill_map(const Table<3,T> &coefficients, + const std_cxx11::function(const TableIndices<3> &)> &predicate, + std::map > &pred_to_values) + { + for (unsigned int i = 0; i < coefficients.size(0); i++) + for (unsigned int j = 0; j < coefficients.size(1); j++) + for (unsigned int k = 0; k < coefficients.size(2); k++) + { + const TableIndices<3> ind(i,j,k); + fill_map_index(coefficients,ind,predicate,pred_to_values); + } + } + + + template + double complex_mean_value(const T &value) + { + return value; + } + + template + double complex_mean_value(const std::complex &value) + { + AssertThrow(false, ExcMessage("FESeries::process_coefficients() can not be used with" + "complex-valued coefficients and VectorTools::mean norm.")); + return std::abs(value); + } + +} + + +template +std::pair,std::vector > +FESeries::process_coefficients(const Table &coefficients, + const std_cxx11::function(const TableIndices &)> &predicate, + const VectorTools::NormType norm) +{ + std::vector predicate_values; + std::vector norm_values; + + // first, parse all table elements into a map of predicate values and coefficients. + // We could have stored (predicate values ->TableIndicies) map, but its + // processing would have been much harder later on. + std::map > pred_to_values; + fill_map(coefficients,predicate,pred_to_values); + + // now go through the map and populate the @p norm_values based on @p norm: + for (typename std::map >::const_iterator it = pred_to_values.begin(); + it != pred_to_values.end(); it++) + { + predicate_values.push_back(it->first); + Vector values(it->second.begin(), + it->second.end()); + + switch (norm) + { + case VectorTools::L2_norm: + { + norm_values.push_back(values.l2_norm()); + break; + } + case VectorTools::L1_norm: + { + norm_values.push_back(values.l1_norm()); + break; + } + case VectorTools::Linfty_norm: + { + norm_values.push_back(values.linfty_norm()); + break; + } + case VectorTools::mean: + { + norm_values.push_back(complex_mean_value(values.mean_value())); + break; + } + default: + AssertThrow(false, ExcNotImplemented()); + break; + } + } + + return std::make_pair(predicate_values,norm_values); +} + + +#endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif // dealii__fe_series_H diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index 7d464fde9b..0e8f549797 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -42,6 +42,7 @@ SET(_src fe_rannacher_turek.cc fe_raviart_thomas.cc fe_raviart_thomas_nodal.cc + fe_series.cc fe_system.cc fe_tools.cc fe_tools_interpolate.cc diff --git a/source/fe/fe_series.cc b/source/fe/fe_series.cc new file mode 100644 index 0000000000..c647a7a2dc --- /dev/null +++ b/source/fe/fe_series.cc @@ -0,0 +1,423 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +#include +#include + +#include +#include + +#include +#ifdef DEAL_II_WITH_GSL +#include +#endif + + +DEAL_II_NAMESPACE_OPEN + +namespace FESeries +{ + + /*-------------- Fourier -------------------------------*/ + + void set_k_vectors(Table<1, Tensor<1,1> > &k_vectors, + const unsigned int N) + { + k_vectors.reinit(TableIndices<1>(N)); + for (unsigned int i=0; i > &k_vectors, + const unsigned int N) + { + k_vectors.reinit(TableIndices<2>(N,N)); + for (unsigned int i=0; i > &k_vectors, + const unsigned int N) + { + k_vectors.reinit(TableIndices<3>(N,N,N)); + for (unsigned int i=0; i + Fourier::Fourier(const unsigned int N, + const hp::FECollection &fe_collection, + const hp::QCollection &q_collection) + : + fe_collection(&fe_collection), + q_collection(&q_collection), + fourier_transform_matrices(fe_collection.size()) + { + set_k_vectors(k_vectors,N); + unrolled_coefficients.resize(k_vectors.n_elements()); + } + + template + void Fourier::calculate(const Vector &local_dof_values, + const unsigned int cell_active_fe_index, + Table > &fourier_coefficients) + { + ensure_existence(cell_active_fe_index); + const FullMatrix > &matrix = fourier_transform_matrices[cell_active_fe_index]; + + std::fill(unrolled_coefficients.begin(), + unrolled_coefficients.end(), + std::complex(0.)); + + Assert (unrolled_coefficients.size() == matrix.m(), + ExcInternalError()); + + Assert (local_dof_values.size() == matrix.n(), + ExcDimensionMismatch(local_dof_values.size(),matrix.n())); + + for (unsigned int i = 0; i < unrolled_coefficients.size(); i++) + for (unsigned int j = 0; j < local_dof_values.size(); j++) + unrolled_coefficients[i] += matrix[i][j] * + local_dof_values[j]; + + fourier_coefficients.fill(unrolled_coefficients.begin()); + } + + template + std::complex integrate(const FiniteElement &fe, + const Quadrature &quadrature, + const Tensor<1,dim> &k_vector, + const unsigned int j) + { + std::complex sum = 0; + for (unsigned int q=0; q &x_q = quadrature.point(q); + sum += std::exp(std::complex(0,1) * + (k_vector * x_q)) * + fe.shape_value(j,x_q) * + quadrature.weight(q); + } + return sum; + } + + + template <> + void Fourier<2>::ensure_existence(const unsigned int fe) + { + Assert (fe < fe_collection->size(), + ExcIndexRange(fe,0,fe_collection->size())) + + if (fourier_transform_matrices[fe].m() == 0) + { + fourier_transform_matrices[fe].reinit(k_vectors.n_elements(), + (*fe_collection)[fe].dofs_per_cell); + + unsigned int k = 0; + for (unsigned int k1=0; k1 + void Fourier<3>::ensure_existence(const unsigned int fe) + { + Assert (fe < fe_collection->size(), + ExcIndexRange(fe,0,fe_collection->size())) + + if (fourier_transform_matrices[fe].m() == 0) + { + fourier_transform_matrices[fe].reinit(k_vectors.n_elements(), + (*fe_collection)[fe].dofs_per_cell); + + unsigned int k = 0; + for (unsigned int k1=0; k1 + void Fourier<1>::ensure_existence(const unsigned int fe) + { + Assert (fe < fe_collection->size(), + ExcIndexRange(fe,0,fe_collection->size())) + + if (fourier_transform_matrices[fe].m() == 0) + { + fourier_transform_matrices[fe].reinit(k_vectors.n_elements(), + (*fe_collection)[fe].dofs_per_cell); + + for (unsigned int k=0; k + double Lh(const Point &x_q, + const TableIndices &indices) + { + double res = 1.0; + for (unsigned int d = 0; d < dim; d++) + { + const double x = 2.0*(x_q[d]-0.5); + Assert ( (x_q[d] <= 1.0) && (x_q[d] >= 0.), + ExcLegendre(d,x_q[d])); + const int ind = indices[d]; +#ifdef DEAL_II_WITH_GSL + res *= sqrt(2.0) * gsl_sf_legendre_Pl (ind, x); +#else + AssertThrow(false, ExcMessage("deal.II has to be configured with GSL" + "in order to use Legendre transformation.")); +#endif + } + return res; + } + + /* + * Multiplier in Legendre coefficients + */ + template + double multiplier(const TableIndices &indices) + { + double res = 1.0; + for (unsigned int d = 0; d < dim; d++) + res *= (0.5+indices[d]); + + return res; + } + + + template + Legendre::Legendre(const unsigned int size_in_each_direction, + const hp::FECollection &fe_collection, + const hp::QCollection &q_collection) + : + N(size_in_each_direction), + fe_collection(&fe_collection), + q_collection(&q_collection), + legendre_transform_matrices(fe_collection.size()), + unrolled_coefficients(Utilities::fixed_power(N), + 0.) + { + } + + template + void Legendre::calculate(const dealii::Vector &local_dof_values, + const unsigned int cell_active_fe_index, + Table &legendre_coefficients) + { + ensure_existence(cell_active_fe_index); + const FullMatrix &matrix = legendre_transform_matrices[cell_active_fe_index]; + + std::fill(unrolled_coefficients.begin(), + unrolled_coefficients.end(), + 0.); + + Assert (unrolled_coefficients.size() == matrix.m(), + ExcInternalError()); + + Assert (local_dof_values.size() == matrix.n(), + ExcDimensionMismatch(local_dof_values.size(),matrix.n())); + + for (unsigned int i = 0; i < unrolled_coefficients.size(); i++) + for (unsigned int j = 0; j < local_dof_values.size(); j++) + unrolled_coefficients[i] += matrix[i][j] * + local_dof_values[j]; + + legendre_coefficients.fill(unrolled_coefficients.begin()); + } + + + template + double integrate_Legendre(const FiniteElement &fe, + const Quadrature &quadrature, + const TableIndices &indices, + const unsigned int dof) + { + double sum = 0; + for (unsigned int q=0; q &x_q = quadrature.point(q); + sum += Lh(x_q, indices) * + fe.shape_value(dof,x_q) * + quadrature.weight(q); + } + return sum * multiplier(indices); + } + + template <> + void Legendre<1>::ensure_existence(const unsigned int fe) + { + Assert (fe < fe_collection->size(), + ExcIndexRange(fe,0,fe_collection->size())) + if (legendre_transform_matrices[fe].m() == 0) + { + legendre_transform_matrices[fe].reinit(N, + (*fe_collection)[fe].dofs_per_cell); + + for (unsigned int k=0; k(k), + j); + } + } + + + + template <> + void Legendre<2>::ensure_existence(const unsigned int fe) + { + Assert (fe < fe_collection->size(), + ExcIndexRange(fe,0,fe_collection->size())) + + if (legendre_transform_matrices[fe].m() == 0) + { + legendre_transform_matrices[fe].reinit(N*N, + (*fe_collection)[fe].dofs_per_cell); + + unsigned int k = 0; + for (unsigned int k1=0; k1(k1,k2), + j); + } + } + + template <> + void Legendre<3>::ensure_existence(const unsigned int fe) + { + Assert (fe < fe_collection->size(), + ExcIndexRange(fe,0,fe_collection->size())) + + if (legendre_transform_matrices[fe].m() == 0) + { + legendre_transform_matrices[fe].reinit(N*N*N, + (*fe_collection)[fe].dofs_per_cell); + + unsigned int k = 0; + for (unsigned int k1=0; k1(k1,k2,k3), + j); + } + } + + /*-------------- linear_regression -------------------------------*/ + std::pair linear_regression(const std::vector &x, + const std::vector &y) + { + FullMatrix K(2,2), invK(2,2); + Vector X(2), B(2); + + Assert (x.size() == y.size(), + ExcMessage("x and y are expected to have the same size")); + + Assert (x.size() >= 2, + dealii::ExcMessage("at least two points are required for linear regression fit")); + + double sum_1 = 0.0, + sum_x = 0.0, + sum_x2= 0.0, + sum_y = 0.0, + sum_xy= 0.0; + + for (unsigned int i = 0; i < x.size(); i++) + { + sum_1 += 1.0; + sum_x += x[i]; + sum_x2 += x[i]*x[i]; + sum_y += y[i]; + sum_xy += x[i]*y[i]; + } + + K(0,0) = sum_1; + K(0,1) = sum_x; + K(1,0) = sum_x; + K(1,1) = sum_x2; + + B(0) = sum_y; + B(1) = sum_xy; + + invK.invert(K); + invK.vmult(X,B,false); + + return std::make_pair(X(1),X(0)); + } + + +} // end of namespace FESeries + + + +/*-------------- Explicit Instantiations -------------------------------*/ +template class FESeries::Fourier<1>; +template class FESeries::Fourier<2>; +template class FESeries::Fourier<3>; +template class FESeries::Legendre<1>; +template class FESeries::Legendre<2>; +template class FESeries::Legendre<3>; + + + +DEAL_II_NAMESPACE_CLOSE diff --git a/tests/fe/fe_series_01.cc b/tests/fe/fe_series_01.cc new file mode 100644 index 0000000000..35854f018a --- /dev/null +++ b/tests/fe/fe_series_01.cc @@ -0,0 +1,98 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test Fourier expansion in 1D for some simple functions. +// Below is the MWE in Maxima +/********************************************************* +a:0; +b:1; +nmax:3; +define(Phi(n), exp(-2*%i*%pi*n*x/(b-a))); +f:x; +C0:integrate(f*conjugate(Phi(0)),x,a,b)/(b-a); +define(C(n),integrate(f*conjugate(Phi(n)),x,a,b)/(b-a)); +fullratsimp(map(C,makelist(i,i,1,nmax))); +fs(nmax):=C0+sum(realpart(conjugate(C(m))*Phi(-m)+C(m)*Phi(m)),m,1,nmax); +plot2d([f,fs(0),fs(1),fs(2),fs(3)],[x,0,1]); +*********************************************************/ + +#include "../tests.h" +#include +#include + +#include +#include +#include + +#include + +using namespace dealii; + +void test_1d () +{ + const unsigned int N =4; + // exact values obtained by Maxima: + std::vector > exact(N); + const std::complex I(0,1.); + exact[0] = std::complex(1./2,0.); + exact[1] = -I/(2*numbers::PI); + exact[2] = -I/(4*numbers::PI); + exact[3] = -I/(6*numbers::PI); + // + const unsigned int dim =1; + hp::FECollection fe_collection; + hp::QCollection q_collection; + + // linear FE + fe_collection.push_back(FE_Q(1)); + + QGauss<1> base_quadrature (6); + QIterated quadrature (base_quadrature, N); + q_collection.push_back(quadrature); + + FESeries::Fourier fourier(N,fe_collection,q_collection); + + Vector local_dof_values(2); + local_dof_values[0] = 0; + local_dof_values[1] = 1.; + const unsigned int cell_active_fe_index =0; + + Table > fourier_coefficients; + fourier_coefficients.reinit(TableIndices<1>(N)); + + fourier.calculate(local_dof_values, + cell_active_fe_index, + fourier_coefficients); + + deallog <<"calculated:" << std::endl; + for (unsigned int i = 0; i < N; i++) + deallog << fourier_coefficients[i].real() << " " << fourier_coefficients[i].imag() << std::endl; + deallog <<"exact:" << std::endl; + for (unsigned int i = 0; i < N; i++) + deallog << exact[i].real() << " " << exact[i].imag() << std::endl; + +} + + + +int main() +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test_1d(); +} diff --git a/tests/fe/fe_series_01.output b/tests/fe/fe_series_01.output new file mode 100644 index 0000000000..8eadf3083b --- /dev/null +++ b/tests/fe/fe_series_01.output @@ -0,0 +1,11 @@ + +DEAL::calculated: +DEAL::0.500000 0 +DEAL::0 -0.159155 +DEAL::0 -0.0795775 +DEAL::0 -0.0530516 +DEAL::exact: +DEAL::0.500000 0 +DEAL::0 -0.159155 +DEAL::0 -0.0795775 +DEAL::0 -0.0530516 diff --git a/tests/fe/fe_series_02.cc b/tests/fe/fe_series_02.cc new file mode 100644 index 0000000000..7821c31f81 --- /dev/null +++ b/tests/fe/fe_series_02.cc @@ -0,0 +1,207 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test Legendre expansion in 1D for a function given using Legendre functions. + +#include "../tests.h" +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +/** + * A 1D function given by input legendre coefficients + */ +template +class LegendreFunction : public Function +{ +public: + LegendreFunction(const std::vector coefficients) + : + dealii::Function(1), + coefficients(coefficients) + {} + + virtual double value(const dealii::Point &point, + const unsigned int component = 0 ) const; + + const std::vector &get_coefficients() const + { + return coefficients; + } + +private: + + const std::vector coefficients; +}; + +template +double LegendreFunction::value(const dealii::Point &point, + const unsigned int ) const +{ + Assert(dim==1, + dealii::ExcNotImplemented()); + + double f = 0.0; + + for (int l = 0; l < int(coefficients.size()); l++) + { + const double m = 0.5; // mid-point + const double h = 0.5; // half-length + const double x = (point[0]-m)/h; // 1D only + f += sqrt(1.0/h) * gsl_sf_legendre_Pl (l, x) * + coefficients[l]; + } + + return f; +} + +template +void test(const LegendreFunction &func, + const unsigned int poly_degree) +{ + Triangulation triangulation; + hp::DoFHandler dof_handler(triangulation); + hp::FECollection fe_collection; + fe_collection.push_back(dealii::FE_Q(poly_degree)); + + hp::QCollection quadrature_formula; + quadrature_formula.push_back(QGauss(poly_degree+6)); + + // reference cell: + GridGenerator::hyper_cube (triangulation,0.0,1.0); + + dof_handler.distribute_dofs (fe_collection); + + Vector values(dof_handler.n_dofs()); + + VectorTools::interpolate (dof_handler, + func, + values); + + const unsigned int N = poly_degree+1; + FESeries::Legendre legendre(N, + fe_collection, + quadrature_formula); + + const std::vector &coeff_in = func.get_coefficients(); + Table<1,double> coeff_out(N); + + Vector local_dof_values; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + + { + const unsigned int cell_n_dofs = cell->get_fe().dofs_per_cell; + const unsigned int cell_active_fe_index = cell->active_fe_index(); + + local_dof_values.reinit (cell_n_dofs); + cell->get_dof_values (values, local_dof_values); + + legendre.calculate(local_dof_values, + cell_active_fe_index, + coeff_out); + } + + for (unsigned int i = 0; i < coeff_in.size(); i++) + deallog << coeff_in[i] <<" "; + + deallog << std::endl; + + for (unsigned int i = 0; i < N; i++) + deallog << coeff_out[i] <<" "; + + deallog << std::endl; + + dof_handler.clear(); +} + +int main () +{ + const int dim = 1; + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + { + std::vector coeff_in(2); + coeff_in[0] = 1.0; + coeff_in[1] = 2.0; + LegendreFunction function(coeff_in); + test(function,1); + } + + { + std::vector coeff_in(3); + coeff_in[0] = 1.0; + coeff_in[1] = 2.0; + coeff_in[2] = 3.0; + LegendreFunction function(coeff_in); + test(function,2); + } + + { + std::vector coeff_in(4); + coeff_in[0] = 1.0; + coeff_in[1] = 2.0; + coeff_in[2] = 3.0; + coeff_in[3] = 4.0; + LegendreFunction function(coeff_in); + test(function,3); + } + + { + std::vector coeff_in(5); + coeff_in[0] = 1.0; + coeff_in[1] = 2.0; + coeff_in[2] = 3.0; + coeff_in[3] = 4.0; + coeff_in[4] = 5.0; + LegendreFunction function(coeff_in); + test(function,4); + } + + { + std::vector coeff_in(6); + coeff_in[0] = 1.0; + coeff_in[1] = 2.0; + coeff_in[2] = 3.0; + coeff_in[3] = 4.0; + coeff_in[4] = 5.0; + coeff_in[5] = 6.0; + LegendreFunction function(coeff_in); + test(function,5); + } + + dealii::deallog << "Ok"< +#include + +#include +#include + +using namespace dealii; + +std::pair +pred_ind(const TableIndices<2> &ind) +{ + return std::make_pair(true,ind[0]+ind[1]); +} + +void test2d (const VectorTools::NormType norm) +{ + const unsigned int dim = 2; + const unsigned int N=4; + Table coefficients(4,4); + for (unsigned int i = 0; i < N; i++) + for (unsigned int j = 0; j < N; j++) + coefficients(i,j) = i*N+j; + + std::pair,std::vector > res = + FESeries::process_coefficients<2,double>(coefficients,pred_ind,norm); + + for (unsigned int i = 0; i < res.first.size(); i++) + deallog << res.first[i] << " : " << res.second[i] << std::endl; +} + + + +int main() +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + deallog << "L2_norm" << std::endl; + test2d(VectorTools::L2_norm); + deallog << "L1_norm" << std::endl; + test2d(VectorTools::L1_norm); + deallog << "Linfty_norm" << std::endl; + test2d(VectorTools::Linfty_norm); + deallog << "mean" << std::endl; + test2d(VectorTools::mean); +} diff --git a/tests/fe/fe_series_03.output b/tests/fe/fe_series_03.output new file mode 100644 index 0000000000..0dbfbbeb1f --- /dev/null +++ b/tests/fe/fe_series_03.output @@ -0,0 +1,33 @@ + +DEAL::L2_norm +DEAL::0 : 0 +DEAL::1 : 4.12311 +DEAL::2 : 9.64365 +DEAL::3 : 16.4317 +DEAL::4 : 17.8326 +DEAL::5 : 17.8045 +DEAL::6 : 15.0000 +DEAL::L1_norm +DEAL::0 : 0 +DEAL::1 : 5.00000 +DEAL::2 : 15.0000 +DEAL::3 : 30.0000 +DEAL::4 : 30.0000 +DEAL::5 : 25.0000 +DEAL::6 : 15.0000 +DEAL::Linfty_norm +DEAL::0 : 0 +DEAL::1 : 4.00000 +DEAL::2 : 8.00000 +DEAL::3 : 12.0000 +DEAL::4 : 13.0000 +DEAL::5 : 14.0000 +DEAL::6 : 15.0000 +DEAL::mean +DEAL::0 : 0 +DEAL::1 : 2.50000 +DEAL::2 : 5.00000 +DEAL::3 : 7.50000 +DEAL::4 : 10.0000 +DEAL::5 : 12.5000 +DEAL::6 : 15.0000 diff --git a/tests/fe/fe_series_04.cc b/tests/fe/fe_series_04.cc new file mode 100644 index 0000000000..9480045a54 --- /dev/null +++ b/tests/fe/fe_series_04.cc @@ -0,0 +1,203 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test Legendre expansion in 1D for quadratic function coming from FE. +// Also test that our interpretation of GSL function is correct to have +// an orthogonal basis. + +// MWE in Maxima +/************************************** +hj : 1/2; +mj : 1/2; +define(f(x), (1.81735e-05*(1.0-x)*(0.5-x)*2 + 0.000901649*x*(x-0.5)*2 + 1.35059e-05*x*(1.0-x)*4.0)); +load("orthopoly"); +orthopoly_returns_intervals : false; +plot2d([legendre_p(0,x), legendre_p(1,x),legendre_p(2,x)], [x,-1,1])$ +define(Lh(n,h,m,x), sqrt(1/h)*legendre_p(n,(x-m)/h)); +define(C(n),integrate(f(x)*Lh(n,hj,mj,x),x,0,1)*(n+1/2)); +bfloat(C(0)), nouns; +bfloat(C(1)), nouns; +bfloat(C(2)), nouns; +bfloat(C(3)), nouns; + **************************************/ + +#include "../tests.h" +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +template +class LegendreFunction : public Function +{ +public: + LegendreFunction() + : + Function(1) + {} + + virtual double value(const dealii::Point &point, + const unsigned int component = 0 ) const; +}; + +template +double LegendreFunction::value(const Point &point, + const unsigned int ) const +{ + Assert(dim==1, + dealii::ExcNotImplemented()); + + const double &x = point[0]; + return 1.81735e-05*(1.0-x)*(0.5-x)*2 + 0.000901649*x*(x-0.5)*2 + 1.35059e-05*x*(1.0-x)*4.0; +} + + +template +void test(const LegendreFunction &func, + const unsigned int poly_degree) +{ + Triangulation triangulation; + hp::DoFHandler dof_handler(triangulation); + hp::FECollection fe_collection; + hp::QCollection quadrature_formula; + + for (unsigned int p = poly_degree; p<= poly_degree+3; p++) + { + fe_collection.push_back(dealii::FE_Q(p)); + quadrature_formula.push_back(dealii::QGauss(p+1+5)); + } + + // reference cell + GridGenerator::hyper_cube (triangulation,0.0,1.0); + + dof_handler.distribute_dofs (fe_collection); + + Vector values(dof_handler.n_dofs()); + + VectorTools::interpolate (dof_handler,func,values); + const unsigned int N = 4; + FESeries::Legendre legendre(N,fe_collection,quadrature_formula); + + Table<1,double> coeff_out(N); + Vector local_dof_values; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + { + const unsigned int cell_n_dofs = cell->get_fe().dofs_per_cell; + const unsigned int cell_active_fe_index = cell->active_fe_index(); + + local_dof_values.reinit (cell_n_dofs); + cell->get_dof_values (values, local_dof_values); + + legendre.calculate(local_dof_values, + cell_active_fe_index, + coeff_out); + + deallog << "local dofs:"; + for (unsigned int i = 0; i < cell_n_dofs; i++) + dealii::deallog << " " < coeff_exp(3); + // coeff calculated in maxima (see MWE above): + coeff_exp[0] = 1.147688635236788e-4; + coeff_exp[1] = 3.123557585310879e-4; + coeff_exp[2] = 2.104375000953028e-4; + deallog << "exact:"< quadrature (8); + deallog <<"orthogonality: " << std::endl; + for (int k1 = 0; k1 < N; k1++) + for (int k2 = 0; k2 < N; k2++) + { + double ortho = 0; + for (unsigned int q=0; q &x_q = quadrature.point(q); + const double m = 0.5; // mid-point + const double h = 0.5; // half-length + const double x = (x_q[0]-m)/h; // 1D only + Assert (std::fabs(x) < 1.0, + dealii::ExcInternalError()); + const double L1 = std::sqrt(1.0/h) * gsl_sf_legendre_Pl (k1, x); + const double L2 = std::sqrt(1.0/h) * gsl_sf_legendre_Pl (k2, x); + ortho += L1 * L2 * quadrature.weight(q); + } + ortho *=(1.0+k1+k2)/2.0; + + deallog << "("< function; + test(function,2); + + dealii::deallog << "Ok"< +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +template +class LegendreFunction : public Function +{ +public: + LegendreFunction(const Table &coefficients) + : + dealii::Function(1), + coefficients(coefficients) + { + } + + virtual double value(const Point &point, + const unsigned int component = 0 ) const; + + const Table &get_coefficients() const + { + return coefficients; + } + +private: + const Table coefficients; +}; + +// copy-paste from fe_series.cc +template +double Lh(const Point &x_q, + const TableIndices &indices) +{ + double res = 1.0; + for (unsigned int d = 0; d < dim; d++) + { + const double x = 2.0*(x_q[d]-0.5); + Assert ( (x_q[d] <= 1.0) && (x_q[d] >= 0.), + ExcMessage("x_q is not in [0,1]" + + std::to_string(x_q[d]))); + const int ind = indices[d]; + res *= sqrt(2.0) * gsl_sf_legendre_Pl (ind, x); + } + return res; +} + +template<> +double LegendreFunction<2>::value(const dealii::Point<2> &point, + const unsigned int ) const +{ + double f = 0.0; + + for (unsigned int i = 0; i < coefficients.size(0); i++) + for (unsigned int j = 0; j < coefficients.size(1); j++) + f+= Lh(point, TableIndices<2>(i,j)) * coefficients(i,j); + + return f; +} + +template<> +double LegendreFunction<3>::value(const dealii::Point<3> &point, + const unsigned int ) const +{ + double f = 0.0; + + for (unsigned int i = 0; i < coefficients.size(0); i++) + for (unsigned int j = 0; j < coefficients.size(1); j++) + for (unsigned int k = 0; k < coefficients.size(2); k++) + f+= Lh(point, TableIndices<3>(i,j,k)) * coefficients(i,j,k); + + return f; +} + +void print(const Table<2,double> &coeff) +{ + for (unsigned int i = 0; i < coeff.size(0); i++) + for (unsigned int j = 0; j < coeff.size(1); j++) + deallog << coeff(i,j) << " "; + deallog << std::endl; +} + +void print(const Table<3,double> &coeff) +{ + for (unsigned int i = 0; i < coeff.size(0); i++) + for (unsigned int j = 0; j < coeff.size(1); j++) + for (unsigned int k = 0; k < coeff.size(2); k++) + deallog << coeff(i,j,k) << " "; + deallog << std::endl; +} + +void resize(Table<2,double> &coeff, const unsigned int N) +{ + coeff.reinit(N,N); +} + +void resize(Table<3,double> &coeff, const unsigned int N) +{ + TableIndices<3> size; + for (unsigned int d=0; d<3; d++) + size[d] = N; + coeff.reinit(size); +} + + + +template +void test(const LegendreFunction &func, + const unsigned int poly_degree) +{ + const unsigned int max_poly = poly_degree+3; + deallog <<"-----------------------------------"< triangulation; + hp::DoFHandler dof_handler(triangulation); + hp::FECollection fe_collection; + hp::QCollection quadrature_formula; + + // add some extra FEs in fe_collection + for (unsigned int p = 1; p <= max_poly; p++) + { + fe_collection.push_back(FE_Q(p)); + quadrature_formula.push_back(QGauss(p+1+5)); + } + + GridGenerator::hyper_cube (triangulation,0.0,1.0); // reference cell + const unsigned int fe_index = poly_degree-1; + dof_handler.begin_active()->set_active_fe_index(fe_index); + dof_handler.distribute_dofs (fe_collection); + + Vector values(dof_handler.n_dofs()); + + VectorTools::interpolate (dof_handler,func,values); + + const unsigned int N = poly_degree+1; + FESeries::Legendre legendre(N, + fe_collection, + quadrature_formula); + + const Table &coeff_in = func.get_coefficients(); + Table coeff_out; + resize(coeff_out,N); + + Vector local_dof_values; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + { + const unsigned int cell_n_dofs = cell->get_fe().dofs_per_cell; + const unsigned int cell_active_fe_index = cell->active_fe_index(); + + local_dof_values.reinit (cell_n_dofs); + cell->get_dof_values (values, local_dof_values); + + legendre.calculate(local_dof_values, + cell_active_fe_index, + coeff_out); + } + + deallog << "calculated:" << std::endl; + print(coeff_out); + deallog <<"exact:"< coeff_in(coeff_1d,coeff_1d); + unsigned int ind = 0; + for (unsigned int i = 0; i < coeff_1d; i++) + for (unsigned int j = 0; j < coeff_1d; j++) + coeff_in(i,j) = 1.0 + ind++; + + LegendreFunction function(coeff_in); + test(function,p); + } + + { + const unsigned int dim = 2; + const unsigned int coeff_1d = 3; + const unsigned int p = 2; + Table coeff_in(coeff_1d,coeff_1d); + unsigned int ind = 0; + for (unsigned int i = 0; i < coeff_1d; i++) + for (unsigned int j = 0; j < coeff_1d; j++) + coeff_in(i,j) = 1.0 + ind++; + + LegendreFunction function(coeff_in); + test(function,p); + } + + { + const unsigned int dim = 3; + const unsigned int coeff_1d = 2; + const unsigned int p = 1; + Table coeff_in(coeff_1d,coeff_1d,coeff_1d); + unsigned int ind = 0; + for (unsigned int i = 0; i < coeff_1d; i++) + for (unsigned int j = 0; j < coeff_1d; j++) + for (unsigned int k = 0; k < coeff_1d; k++) + coeff_in(i,j,k) = 1.0 + ind++; + + LegendreFunction function(coeff_in); + test(function,p); + } + + { + const unsigned int dim = 3; + const unsigned int coeff_1d = 3; + const unsigned int p = 2; + Table coeff_in(coeff_1d,coeff_1d,coeff_1d); + unsigned int ind = 0; + for (unsigned int i = 0; i < coeff_1d; i++) + for (unsigned int j = 0; j < coeff_1d; j++) + for (unsigned int k = 0; k < coeff_1d; k++) + coeff_in(i,j,k) = 1.0 + ind++; + + LegendreFunction function(coeff_in); + test(function,p); + } + + dealii::deallog << "Ok"< +#include + +#include +#include +#include + +#include + +using namespace dealii; + +void test_2d () +{ + const unsigned int dim = 2; + const unsigned int N = 7; + hp::FECollection fe_collection; + hp::QCollection fourier_q_collection; + + for (unsigned int degree=2; degree<=N; ++degree) + { + fe_collection.push_back (FE_Q(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); + + FESeries::Fourier fourier(N,fe_collection,fourier_q_collection); + Table > fourier_coefficients; + fourier_coefficients.reinit(N,N); + + Vector local_dof_values(9); + double dofs[] = {0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 2.3801522930483391e-04, 0.0000000000000000e+00, 1.1949018981806140e-04, 0.0000000000000000e+00, 1.1949019042912971e-04, 5.9982796422221083e-05}; + for (unsigned int i = 0; i < 9; i++) + local_dof_values[i] = dofs[i]; + + const unsigned int cell_active_fe_index =0; + fourier.calculate(local_dof_values, + cell_active_fe_index, + fourier_coefficients); + + for (unsigned int i = 0; i < fourier_coefficients.size(0); i++) + for (unsigned int j = 0; j < fourier_coefficients.size(1); j++) + if ((i*i+j*j < N*N) && ( i*i+j*j>0)) + deallog << (i*i+j*j) + << " : " << fourier_coefficients(i,j) << std::endl; + +} + + + +int main() +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test_2d(); +} diff --git a/tests/fe/fe_series_06.output b/tests/fe/fe_series_06.output new file mode 100644 index 0000000000..406b5bfe0d --- /dev/null +++ b/tests/fe/fe_series_06.output @@ -0,0 +1,43 @@ + +DEAL::1 : (-4.84334e-08,-1.89889e-05) +DEAL::4 : (-1.22068e-08,-9.47054e-06) +DEAL::9 : (-5.67727e-09,-6.23164e-06) +DEAL::16 : (-3.80695e-09,-4.45243e-06) +DEAL::25 : (-4.22667e-09,-2.96277e-06) +DEAL::36 : (-1.22216e-08,-3.16485e-07) +DEAL::1 : (-4.84334e-08,-1.89889e-05) +DEAL::2 : (-6.02728e-06,3.11374e-08) +DEAL::5 : (-3.00601e-06,1.16886e-08) +DEAL::10 : (-1.97796e-06,6.93416e-09) +DEAL::17 : (-1.41323e-06,4.87421e-09) +DEAL::26 : (-9.40404e-07,3.78777e-09) +DEAL::37 : (-1.00491e-07,4.18806e-09) +DEAL::4 : (-1.22068e-08,-9.47054e-06) +DEAL::5 : (-3.00601e-06,1.16886e-08) +DEAL::8 : (-1.49921e-06,3.91395e-09) +DEAL::13 : (-9.86482e-07,2.19786e-09) +DEAL::20 : (-7.04830e-07,1.53037e-09) +DEAL::29 : (-4.69014e-07,1.28983e-09) +DEAL::40 : (-5.01095e-08,2.02474e-09) +DEAL::9 : (-5.67728e-09,-6.23164e-06) +DEAL::10 : (-1.97796e-06,6.93416e-09) +DEAL::13 : (-9.86482e-07,2.19786e-09) +DEAL::18 : (-6.49106e-07,1.19779e-09) +DEAL::25 : (-4.63779e-07,8.29497e-10) +DEAL::34 : (-3.08612e-07,7.30608e-10) +DEAL::45 : (-3.29704e-08,1.31967e-09) +DEAL::16 : (-3.80696e-09,-4.45243e-06) +DEAL::17 : (-1.41323e-06,4.87421e-09) +DEAL::20 : (-7.04830e-07,1.53036e-09) +DEAL::25 : (-4.63779e-07,8.29497e-10) +DEAL::32 : (-3.31364e-07,5.73868e-10) +DEAL::41 : (-2.20499e-07,5.09503e-10) +DEAL::25 : (-4.22667e-09,-2.96277e-06) +DEAL::26 : (-9.40404e-07,3.78777e-09) +DEAL::29 : (-4.69014e-07,1.28983e-09) +DEAL::34 : (-3.08612e-07,7.30608e-10) +DEAL::41 : (-2.20499e-07,5.09503e-10) +DEAL::36 : (-1.22216e-08,-3.16485e-07) +DEAL::37 : (-1.00491e-07,4.18806e-09) +DEAL::40 : (-5.01095e-08,2.02474e-09) +DEAL::45 : (-3.29704e-08,1.31967e-09) diff --git a/tests/fe/fe_series_07.cc b/tests/fe/fe_series_07.cc new file mode 100644 index 0000000000..8ef6c0b8f6 --- /dev/null +++ b/tests/fe/fe_series_07.cc @@ -0,0 +1,70 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// same as fe_series_03, but with ignoring some coefficients by predicate. + +#include "../tests.h" +#include +#include + +#include +#include + +using namespace dealii; + +std::pair +pred_ind(const TableIndices<2> &ind) +{ + const unsigned int val = ind[0]+ind[1]; + if (val < 4) + return std::make_pair(true,val); + else + return std::make_pair(false,val); +} + +void test2d (const VectorTools::NormType norm) +{ + const unsigned int dim = 2; + const unsigned int N=4; + Table coefficients(4,4); + for (unsigned int i = 0; i < N; i++) + for (unsigned int j = 0; j < N; j++) + coefficients(i,j) = i*N+j; + + std::pair,std::vector > res = + FESeries::process_coefficients<2,double>(coefficients,pred_ind,norm); + + for (unsigned int i = 0; i < res.first.size(); i++) + deallog << res.first[i] << " : " << res.second[i] << std::endl; +} + + + +int main() +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + deallog << "L2_norm" << std::endl; + test2d(VectorTools::L2_norm); + deallog << "L1_norm" << std::endl; + test2d(VectorTools::L1_norm); + deallog << "Linfty_norm" << std::endl; + test2d(VectorTools::Linfty_norm); + deallog << "mean" << std::endl; + test2d(VectorTools::mean); +} diff --git a/tests/fe/fe_series_07.output b/tests/fe/fe_series_07.output new file mode 100644 index 0000000000..e7955ac589 --- /dev/null +++ b/tests/fe/fe_series_07.output @@ -0,0 +1,21 @@ + +DEAL::L2_norm +DEAL::0 : 0 +DEAL::1 : 4.12311 +DEAL::2 : 9.64365 +DEAL::3 : 16.4317 +DEAL::L1_norm +DEAL::0 : 0 +DEAL::1 : 5.00000 +DEAL::2 : 15.0000 +DEAL::3 : 30.0000 +DEAL::Linfty_norm +DEAL::0 : 0 +DEAL::1 : 4.00000 +DEAL::2 : 8.00000 +DEAL::3 : 12.0000 +DEAL::mean +DEAL::0 : 0 +DEAL::1 : 2.50000 +DEAL::2 : 5.00000 +DEAL::3 : 7.50000 diff --git a/tests/fe/fe_series_08.cc b/tests/fe/fe_series_08.cc new file mode 100644 index 0000000000..b225909077 --- /dev/null +++ b/tests/fe/fe_series_08.cc @@ -0,0 +1,60 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test FESeries::linear_regression() + +#include "../tests.h" +#include +#include + +#include +#include + +using namespace dealii; + + +void test () +{ + const double k_in = numbers::PI; + const double b_in = std::sqrt(2.); + const unsigned int N = 10; + std::vector x(N), y(N); + + // fill the data + for (unsigned int i = 0; i < N; i++) + { + x[i] = 0.1*i; + y[i] = k_in * x[i] + b_in; + } + + + std::pair fit = FESeries::linear_regression(x,y); + + deallog << "exact: " << k_in << " " << b_in << std::endl; + deallog << "calculated: " << fit.first << " " << fit.second << std::endl; + +} + + + +int main() +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test(); +} diff --git a/tests/fe/fe_series_08.output b/tests/fe/fe_series_08.output new file mode 100644 index 0000000000..556fbbc3e0 --- /dev/null +++ b/tests/fe/fe_series_08.output @@ -0,0 +1,3 @@ + +DEAL::exact: 3.14159 1.41421 +DEAL::calculated: 3.14159 1.41421 diff --git a/tests/hp/step-27.cc b/tests/hp/step-27.cc new file mode 100644 index 0000000000..ec1365ddd1 --- /dev/null +++ b/tests/hp/step-27.cc @@ -0,0 +1,707 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2006 - 2015 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 at + * the top level of the deal.II distribution. + * + * -------------------------------------------------------------------- + */ + + +// A combination of step-27 from 8.4 with corrected k-vectors, that is 2\pi*k instead of \pi*k +// and a new step-27 from 8.5 which use FESeries namespace. By default, the new +// version is used, but the blessed output file is obtained using the +// modified 8.4 version. + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + + +namespace Step27 +{ + using namespace dealii; + + + + template + class LaplaceProblem + { + public: + LaplaceProblem (); + ~LaplaceProblem (); + + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void create_coarse_grid (); + void estimate_smoothness (Vector &smoothness_indicators); + void postprocess (const unsigned int cycle); + + Triangulation triangulation; + + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection quadrature_collection; + hp::QCollection face_quadrature_collection; + + hp::QCollection fourier_q_collection; + std_cxx11::shared_ptr> fourier; + std::vector ln_k; + Table > fourier_coefficients; + + ConstraintMatrix constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + const unsigned int max_degree; + }; + + + + template + class RightHandSide : public Function + { + public: + RightHandSide () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const; + }; + + + template + double + RightHandSide::value (const Point &p, + const unsigned int /*component*/) const + { + double product = 1; + for (unsigned int d=0; d + void resize(Table<2,T> &coeff, const unsigned int N) + { + coeff.reinit(N,N); + } + + + template + LaplaceProblem::LaplaceProblem () + : + dof_handler (triangulation), + max_degree (dim <= 2 ? 7 : 5) + { + for (unsigned int degree=2; degree<=max_degree; ++degree) + { + fe_collection.push_back (FE_Q(degree)); + 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_cxx11::make_shared >(N, + fe_collection, + fourier_q_collection); + resize(fourier_coefficients,N); + } + + + + template + LaplaceProblem::~LaplaceProblem () + { + dof_handler.clear (); + } + + + template + void LaplaceProblem::setup_system () + { + dof_handler.distribute_dofs (fe_collection); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + constraints); + constraints.close (); + + DynamicSparsityPattern dsp (dof_handler.n_dofs(), + dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false); + sparsity_pattern.copy_from (dsp); + + system_matrix.reinit (sparsity_pattern); + } + + + + + template + void LaplaceProblem::assemble_system () + { + hp::FEValues hp_fe_values (fe_collection, + quadrature_collection, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const RightHandSide rhs_function; + + FullMatrix cell_matrix; + Vector cell_rhs; + + std::vector local_dof_indices; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + + cell_matrix.reinit (dofs_per_cell, dofs_per_cell); + cell_matrix = 0; + + cell_rhs.reinit (dofs_per_cell); + cell_rhs = 0; + + hp_fe_values.reinit (cell); + + const FEValues &fe_values = hp_fe_values.get_present_fe_values (); + + std::vector rhs_values (fe_values.n_quadrature_points); + rhs_function.value_list (fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q_point=0; + q_pointget_dof_indices (local_dof_indices); + + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } + } + + + + + template + void LaplaceProblem::solve () + { + SolverControl solver_control (system_rhs.size(), + 1e-8*system_rhs.l2_norm()); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + constraints.distribute (solution); + } + + + + + template + void LaplaceProblem::postprocess (const unsigned int cycle) + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); + KellyErrorEstimator::estimate (dof_handler, + face_quadrature_collection, + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + + + Vector smoothness_indicators (triangulation.n_active_cells()); + estimate_smoothness (smoothness_indicators); + + // Output to VTK + if (false) + { + Vector fe_degrees (triangulation.n_active_cells()); + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + fe_degrees(cell->active_cell_index()) + = fe_collection[cell->active_fe_index()].degree; + } + + DataOut > data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.add_data_vector (estimated_error_per_cell, "error"); + data_out.add_data_vector (smoothness_indicators, "smoothness"); + data_out.add_data_vector (fe_degrees, "fe_degree"); + data_out.build_patches (); + + const std::string filename = "solution-" + + Utilities::int_to_string (cycle, 2) + + ".vtk"; + std::ofstream output (filename.c_str()); + data_out.write_vtk (output); + } + + { + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); + + float max_smoothness = *std::min_element (smoothness_indicators.begin(), + smoothness_indicators.end()), + min_smoothness = *std::max_element (smoothness_indicators.begin(), + smoothness_indicators.end()); + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->refine_flag_set()) + { + max_smoothness = std::max (max_smoothness, + smoothness_indicators(cell->active_cell_index())); + min_smoothness = std::min (min_smoothness, + smoothness_indicators(cell->active_cell_index())); + } + } + const float threshold_smoothness = (max_smoothness + min_smoothness) / 2; + + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->refine_flag_set() + && + (smoothness_indicators(cell->active_cell_index()) > threshold_smoothness) + && + (cell->active_fe_index()+1 < fe_collection.size())) + { + cell->clear_refine_flag(); + cell->set_active_fe_index (cell->active_fe_index() + 1); + } + } + + triangulation.execute_coarsening_and_refinement (); + } + } + + + + template <> + void LaplaceProblem<2>::create_coarse_grid () + { + const unsigned int dim = 2; + + static const Point<2> vertices_1[] + = { Point<2> (-1., -1.), + Point<2> (-1./2, -1.), + Point<2> (0., -1.), + Point<2> (+1./2, -1.), + Point<2> (+1, -1.), + + Point<2> (-1., -1./2.), + Point<2> (-1./2, -1./2.), + Point<2> (0., -1./2.), + Point<2> (+1./2, -1./2.), + Point<2> (+1, -1./2.), + + Point<2> (-1., 0.), + Point<2> (-1./2, 0.), + Point<2> (+1./2, 0.), + Point<2> (+1, 0.), + + Point<2> (-1., 1./2.), + Point<2> (-1./2, 1./2.), + Point<2> (0., 1./2.), + Point<2> (+1./2, 1./2.), + Point<2> (+1, 1./2.), + + Point<2> (-1., 1.), + Point<2> (-1./2, 1.), + Point<2> (0., 1.), + Point<2> (+1./2, 1.), + Point<2> (+1, 1.) + }; + const unsigned int + n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]); + const std::vector > vertices (&vertices_1[0], + &vertices_1[n_vertices]); + static const int cell_vertices[][GeometryInfo::vertices_per_cell] + = {{0, 1, 5, 6}, + {1, 2, 6, 7}, + {2, 3, 7, 8}, + {3, 4, 8, 9}, + {5, 6, 10, 11}, + {8, 9, 12, 13}, + {10, 11, 14, 15}, + {12, 13, 17, 18}, + {14, 15, 19, 20}, + {15, 16, 20, 21}, + {16, 17, 21, 22}, + {17, 18, 22, 23} + }; + const unsigned int + n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]); + + std::vector > cells (n_cells, CellData()); + for (unsigned int i=0; i::vertices_per_cell; + ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + triangulation.create_triangulation (vertices, + cells, + SubCellData()); + triangulation.refine_global (3); + } + + + + + + template + void LaplaceProblem::run () + { + for (unsigned int cycle=0; cycle<6; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + create_coarse_grid (); + + setup_system (); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl + << " Number of constraints : " + << constraints.n_constraints() + << std::endl; + + assemble_system (); + solve (); + 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(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(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 > > + fourier_transform_matrices (fe_collection.size()); + + QGauss<1> base_quadrature (2); + QIterated quadrature (base_quadrature, N); + + + for (unsigned int fe=0; fe sum = 0; + for (unsigned int q=0; q 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; fget_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; factive_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 + } +} + + + +int main () +{ + try + { + using namespace dealii; + using namespace Step27; + + LaplaceProblem<2> laplace_problem; + laplace_problem.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/hp/step-27.output b/tests/hp/step-27.output new file mode 100644 index 0000000000..c931646281 --- /dev/null +++ b/tests/hp/step-27.output @@ -0,0 +1,24 @@ +Cycle 0: + Number of active cells: 768 + Number of degrees of freedom: 3264 + Number of constraints : 384 +Cycle 1: + Number of active cells: 966 + Number of degrees of freedom: 5245 + Number of constraints : 936 +Cycle 2: + Number of active cells: 1143 + Number of degrees of freedom: 8441 + Number of constraints : 1929 +Cycle 3: + Number of active cells: 1356 + Number of degrees of freedom: 12349 + Number of constraints : 3046 +Cycle 4: + Number of active cells: 1644 + Number of degrees of freedom: 18178 + Number of constraints : 4713 +Cycle 5: + Number of active cells: 1728 + Number of degrees of freedom: 22591 + Number of constraints : 6095 -- 2.39.5