<h3>General</h3>
<ol>
+
+ <li> 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.
+ <br>
+ (Denis Davydov, 2016/06/23)
+ </li>
+
<li> New: The tutorial step-55 shows how to solve the Stokes system
in parallel with PETSc or Trilinos.
<br>
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}
|\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
$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}$
<h4>What we have to do</h4>
So what do we have to do to estimate the local smoothness of $u({\bf x})$ on
-a cell $K$? Clearly, the first step is to compute the Fourier 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
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,
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
\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
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
* ---------------------------------------------------------------------
*
- * 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.
*/
// 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 <i>hp</i>
+// already discussed in step-11 and step-18. The second and third provide <i>hp</i>
// 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 <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_values.h>
+#include <deal.II/fe/fe_series.h>
// The last set of include files are standard C++ headers. We need support for
// complex numbers when we compute the Fourier transform.
void assemble_system ();
void solve ();
void create_coarse_grid ();
- void estimate_smoothness (Vector<float> &smoothness_indicators) const;
+ void estimate_smoothness (Vector<float> &smoothness_indicators);
void postprocess (const unsigned int cycle);
+ std::pair<bool,unsigned int> predicate(const TableIndices<dim> &indices);
Triangulation<dim> triangulation;
hp::QCollection<dim> quadrature_collection;
hp::QCollection<dim-1> face_quadrature_collection;
+ hp::QCollection<dim> fourier_q_collection;
+ std_cxx11::shared_ptr<FESeries::Fourier<dim>> fourier;
+ std::vector<double> ln_k;
+ Table<dim,std::complex<double> > fourier_coefficients;
+
ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
// face quadrature objects. We start with quadratic elements, and each
// quadrature formula is chosen so that it is appropriate for the matching
// finite element in the hp::FECollection object.
+ //
+ // Finally, we initialize FESeries::Fourier object which will be used to
+ // calculate coefficient in Fourier series as described in the introduction.
+ // In addition to the hp::FECollection, we need to provide quadrature rules
+ // hp::QCollection for integration on the reference cell.
+ //
+ // In order to resize fourier_coefficients Table, we use the following
+ // auxiliary function
+ template <int dim,typename T>
+ void resize(Table<dim,T> &coeff, const unsigned int N)
+ {
+ TableIndices<dim> size;
+ for (unsigned int d=0; d<dim; d++)
+ size[d] = N;
+ coeff.reinit(size);
+ }
+
template <int dim>
LaplaceProblem<dim>::LaplaceProblem ()
:
quadrature_collection.push_back (QGauss<dim>(degree+1));
face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
}
+
+ // As described in the introduction, we define the Fourier vectors ${\bf
+ // k}$ for which we want to compute Fourier coefficients of the solution
+ // on each cell as follows. In 2d, we will need coefficients corresponding to
+ // vectors ${\bf k}=(2 \pi i, 2\pi j)^T$
+ // for which $\sqrt{i^2+j^2}\le N$, with $i,j$ integers and $N$ being the
+ // maximal polynomial degree we use for the finite elements in this
+ // program. The FESeries::Fourier class' constructor first parameter $N$ defines
+ // the number of coefficients in 1D with the total number of coefficients
+ // being $N^{dim}$. Although we will not use coefficients corresponding to
+ // $\sqrt{i^2+j^2}> N$ and $i+j==0$, the overhead of their calculation is minimal.
+ // The transformation matrices for each FiniteElement will be calculated only
+ // once the first time they are required in the course of hp-adaptive
+ // refinement. Because we work on the unit cell, we can do all this work without a
+ // mapping from reference to real cell and consequently can precalculate
+ // these matrices. The calculation of expansion
+ // coefficients for a particular set of local degrees of freedom on a given
+ // cell then follows as a simple matrix-vector product.
+ // The 3d case is handled analogously.
+ const unsigned int N = max_degree;
+
+ // We will need to assemble the matrices that do the Fourier transforms
+ // for each of the finite elements we deal with, i.e. the matrices ${\cal
+ // F}_{{\bf k},j}$ defined in the introduction. We have to do that for
+ // each of the finite elements in use. To that end we need a quadrature
+ // rule. In this example we use the same quadrature formula for each
+ // finite element, namely that is obtained by iterating a
+ // 2-point Gauss formula as many times as the maximal exponent we use for
+ // the term $e^{i{\bf k}\cdot{\bf x}}$:
+ QGauss<1> base_quadrature (2);
+ QIterated<dim> quadrature (base_quadrature, N);
+ for (unsigned int i = 0; i < fe_collection.size(); i++)
+ fourier_q_collection.push_back(quadrature);
+
+ // Now we are ready to set-up the FESeries::Fourier object
+ fourier = std_cxx11::make_shared<FESeries::Fourier<dim> >(N,
+ fe_collection,
+ fourier_q_collection);
+
+ // We need to resize the matrix of fourier coefficients according to the
+ // number of modes N.
+ resize(fourier_coefficients,N);
}
// @sect4{LaplaceProblem::estimate_smoothness}
+ // As described in the introduction, we will need to take the maximum
+ // absolute value of fourier coefficients which correspond to $k$-vector
+ // $|{\bf k}|= const$. To filter the coefficients Table we
+ // will use the FESeries::process_coefficients() which requires a predicate
+ // to be specified. The predicate should operate on TableIndices and return
+ // a pair of <code>bool</code> and <code>unsigned int</code>. The latter
+ // is the value of the map from TableIndicies to unsigned int. It is
+ // used to define subsets of coefficients from which we search for the one
+ // with highest absolute value, i.e. $l^\infty$-norm. The <code>bool</code>
+ // parameter defines which indices should be used in processing. In the
+ // current case we are interested in coefficients which correspond to
+ // $0 < i*i+j*j < N*N$ and $0 < i*i+j*j+k*k < N*N$ in 2D and 3D, respectively.
+ template <int dim>
+ std::pair<bool,unsigned int>
+ LaplaceProblem<dim>::
+ predicate(const TableIndices<dim> &ind)
+ {
+ unsigned int v = 0;
+ for (unsigned int i = 0; i <dim; i++)
+ v += ind[i]*ind[i];
+ if (v>0 && v < max_degree*max_degree)
+ return std::make_pair(true,v);
+ else
+ return std::make_pair(false,v);
+ }
+
// This last function of significance implements the algorithm to estimate
// the smoothness exponent using the algorithms explained in detail in the
// introduction. We will therefore only comment on those points that are of
template <int dim>
void
LaplaceProblem<dim>::
- estimate_smoothness (Vector<float> &smoothness_indicators) const
+ estimate_smoothness (Vector<float> &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
- // <code>k_vectors</code>. In addition, we store the square of the
- // magnitude of each of these vectors (up to a factor $\pi^2$) in the
- // <code>k_vectors_magnitude</code> 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<Tensor<1,dim> > k_vectors;
- std::vector<unsigned int> k_vectors_magnitude;
- switch (dim)
- {
- case 2:
- {
- for (unsigned int i=0; i<N; ++i)
- for (unsigned int j=0; j<N; ++j)
- if (!((i==0) && (j==0))
- &&
- (i*i + j*j < N*N))
- {
- k_vectors.push_back (Point<dim>(numbers::PI * i,
- numbers::PI * j));
- k_vectors_magnitude.push_back (i*i+j*j);
- }
-
- break;
- }
-
- case 3:
- {
- for (unsigned int i=0; i<N; ++i)
- for (unsigned int j=0; j<N; ++j)
- for (unsigned int k=0; k<N; ++k)
- if (!((i==0) && (j==0) && (k==0))
- &&
- (i*i + j*j + k*k < N*N))
- {
- k_vectors.push_back (Point<dim>(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<double> ln_k (n_fourier_modes);
- for (unsigned int i=0; i<n_fourier_modes; ++i)
- ln_k[i] = std::log (k_vectors[i].norm());
+ // 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.
- // Next, we 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. Note that these matrices are
- // complex-valued, so we can't use the FullMatrix class. Instead, we use
- // the Table class template.
- std::vector<Table<2,std::complex<double> > >
- 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<dim> 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 <code>std::complex@<double@>(0,1)</code>.
-
- // 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<fe_collection.size(); ++fe)
- {
- fourier_transform_matrices[fe].reinit (n_fourier_modes,
- fe_collection[fe].dofs_per_cell);
-
- for (unsigned int k=0; k<n_fourier_modes; ++k)
- for (unsigned int j=0; j<fe_collection[fe].dofs_per_cell; ++j)
- {
- std::complex<double> sum = 0;
- for (unsigned int q=0; q<quadrature.size(); ++q)
- {
- const Point<dim> x_q = quadrature.point(q);
- sum += std::exp(std::complex<double>(0,1) *
- (k_vectors[k] * x_q)) *
- fe_collection[fe].shape_value(j,x_q) *
- quadrature.weight(q);
- }
- fourier_transform_matrices[fe](k,j)
- = sum / std::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<std::complex<double> > fourier_coefficients (n_fourier_modes);
+ // will use the following array as a scratch array in the loop to store
+ // local DoF values:
Vector<double> local_dof_values;
// Then here is the loop:
// <code>local_dof_values</code> array after setting it to the right
// size) and then need to compute the Fourier transform by multiplying
// this vector with the matrix ${\cal F}$ corresponding to this finite
- // element. We need to write out the multiplication by hand because
- // the objects holding the data do not have <code>vmult</code>-like
- // functions declared:
+ // element. This is done by calling FESeries::Fourier::calculate(),
+ // that has to be provided with the <code>local_dof_values</code>,
+ // <code>cell->active_fe_index()</code> and a Table to store coefficients.
local_dof_values.reinit (cell->get_fe().dofs_per_cell);
cell->get_dof_values (solution, local_dof_values);
- 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);
- }
+ 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<unsigned int, double> k_to_max_U_map;
- for (unsigned int f=0; f<n_fourier_modes; ++f)
- if ((k_to_max_U_map.find (k_vectors_magnitude[f]) ==
- k_to_max_U_map.end())
- ||
- (k_to_max_U_map[k_vectors_magnitude[f]] <
- std::abs (fourier_coefficients[f])))
- k_to_max_U_map[k_vectors_magnitude[f]]
- = std::abs (fourier_coefficients[f]);
- // Note that it comes in handy here that we have stored the magnitudes
- // of vectors as integers, since this way we do not have to deal with
- // round-off-sized differences between different values of $|{\bf
- // k}|$.
-
- // As the final task, we have to calculate the various contributions
- // to the formula for $\mu$. We'll only take those Fourier
+ // 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}|$ as explained above:
- 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];
- }
+ // k}|$ and thereby need to use VectorTools::Linfty_norm:
+ std::pair<std::vector<unsigned int>, std::vector<double> > res =
+ FESeries::process_coefficients<dim>(fourier_coefficients,
+ std_cxx11::bind(&LaplaceProblem<dim>::predicate,
+ this,
+ std_cxx11::_1),
+ VectorTools::Linfty_norm);
+
+ Assert (res.first.size() == res.second.size(),
+ ExcInternalError());
+
+ // The first vector in the <code>std::pair</code> will store values of the predicate,
+ // that is $i*i+j*j= const$ or $i*i+j*j+k*k = const$
+ // in 2D or 3D respectively. This
+ // vector will be the same for all the cells so we can calculate
+ // logarithms of the corresponding Fourier vectors $|{\bf k}|$ only once
+ // in the whole hp-refinement cycle:
+ if (ln_k.size() == 0)
+ {
+ ln_k.resize(res.first.size(),0);
+ for (unsigned int f = 0; f < ln_k.size(); f++)
+ ln_k[f] = std::log (2.0*numbers::PI*std::sqrt(1.*res.first[f]));
+ }
+
+ // We have to calculate the logarithms of absolute
+ // values of coefficients and use it in 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<double,double> fit = FESeries::linear_regression(ln_k,res.second);
// The final step is to compute the Sobolev index $s=\mu-\frac d2$ and
// store it in the vector of estimated values for each cell:
- smoothness_indicators(cell->active_cell_index()) = mu - 1.*dim/2;
+ smoothness_indicators(cell->active_cell_index()) = -fit.first - 1.*dim/2;
}
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/std_cxx11/shared_ptr.h>
+#include <deal.II/base/config.h>
+#include <deal.II/base/subscriptor.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/table.h>
+#include <deal.II/base/table_indices.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <vector>
+#include <string>
+
+
+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 <int dim>
+ 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<dim> &fe_collection,
+ const hp::QCollection<dim> &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<double> &local_dof_values,
+ const unsigned int cell_active_fe_index,
+ Table<dim,std::complex<double> > &fourier_coefficients);
+
+ private:
+ /**
+ * hp::FECollection for which transformation matrices will be calculated.
+ */
+ SmartPointer<const hp::FECollection<dim> > fe_collection;
+
+ /**
+ * hp::QCollection used in calculation of transformation matrices.
+ */
+ SmartPointer<const hp::QCollection<dim> > 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<dim, Tensor<1,dim> > k_vectors;
+
+ /**
+ * Transformation matrices for each FiniteElement.
+ */
+ std::vector<FullMatrix<std::complex<double> > > fourier_transform_matrices;
+
+ /**
+ * Auxiliary vector to store unrolled coefficients.
+ */
+ std::vector<std::complex<double> > 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 <int dim>
+ 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<dim> &fe_collection,
+ const hp::QCollection<dim> &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<double> &local_dof_values,
+ const unsigned int cell_active_fe_index,
+ Table<dim,double> &legendre_coefficients);
+
+ private:
+ /**
+ * Number of coefficients in each direction
+ */
+ const unsigned int N;
+
+ /**
+ * hp::FECollection for which transformation matrices will be calculated.
+ */
+ SmartPointer<const hp::FECollection<dim> > fe_collection;
+
+ /**
+ * hp::QCollection used in calculation of transformation matrices.
+ */
+ SmartPointer<const hp::QCollection<dim> > 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<FullMatrix<double> > legendre_transform_matrices;
+
+ /**
+ * Auxiliary vector to store unrolled coefficients.
+ */
+ std::vector<double> 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 <code>bool</code> and <code>unsigned int</code>.
+ * 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 <int dim, typename T>
+ std::pair<std::vector<unsigned int>,std::vector<double> >
+ process_coefficients(const Table<dim,T> &coefficients,
+ const std_cxx11::function<std::pair<bool,unsigned int>(const TableIndices<dim> &)> &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<double,double> linear_regression(const std::vector<double> &x,
+ const std::vector<double> &y);
+
+}
+
+/*@}*/
+
+#ifndef DOXYGEN
+
+// ------------------- inline and template functions ----------------
+
+namespace
+{
+ template <int dim,typename T>
+ void fill_map_index(const Table<dim,T> &coefficients,
+ const TableIndices<dim> &ind,
+ const std_cxx11::function<std::pair<bool,unsigned int>(const TableIndices<dim> &)> &predicate,
+ std::map<unsigned int, std::vector<T> > &pred_to_values)
+ {
+ const std::pair<bool,unsigned int> 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 <typename T>
+ void fill_map(const Table<1,T> &coefficients,
+ const std_cxx11::function<std::pair<bool,unsigned int>(const TableIndices<1> &)> &predicate,
+ std::map<unsigned int, std::vector<T> > &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 <typename T>
+ void fill_map(const Table<2,T> &coefficients,
+ const std_cxx11::function<std::pair<bool,unsigned int>(const TableIndices<2> &)> &predicate,
+ std::map<unsigned int, std::vector<T> > &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 <typename T>
+ void fill_map(const Table<3,T> &coefficients,
+ const std_cxx11::function<std::pair<bool,unsigned int>(const TableIndices<3> &)> &predicate,
+ std::map<unsigned int, std::vector<T> > &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 <typename T>
+ double complex_mean_value(const T &value)
+ {
+ return value;
+ }
+
+ template <typename T>
+ double complex_mean_value(const std::complex<T> &value)
+ {
+ AssertThrow(false, ExcMessage("FESeries::process_coefficients() can not be used with"
+ "complex-valued coefficients and VectorTools::mean norm."));
+ return std::abs(value);
+ }
+
+}
+
+
+template <int dim, typename T>
+std::pair<std::vector<unsigned int>,std::vector<double> >
+FESeries::process_coefficients(const Table<dim,T> &coefficients,
+ const std_cxx11::function<std::pair<bool,unsigned int>(const TableIndices<dim> &)> &predicate,
+ const VectorTools::NormType norm)
+{
+ std::vector<unsigned int> predicate_values;
+ std::vector<double> 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<unsigned int, std::vector<T> > 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<unsigned int, std::vector<T> >::const_iterator it = pred_to_values.begin();
+ it != pred_to_values.end(); it++)
+ {
+ predicate_values.push_back(it->first);
+ Vector<T> 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
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/fe/fe_series.h>
+#include <deal.II/base/numbers.h>
+
+#include <cctype>
+#include <iostream>
+
+#include <deal.II/base/config.h>
+#ifdef DEAL_II_WITH_GSL
+#include <gsl/gsl_sf_legendre.h>
+#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<N; ++i)
+ k_vectors(i)[0] = 2. * numbers::PI * i;
+
+ }
+
+ void set_k_vectors(Table<2, Tensor<1,2> > &k_vectors,
+ const unsigned int N)
+ {
+ k_vectors.reinit(TableIndices<2>(N,N));
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ {
+ k_vectors(i,j)[0] = 2. * numbers::PI * i;
+ k_vectors(i,j)[1] = 2. * numbers::PI * j;
+ }
+ }
+
+ void set_k_vectors(Table<3, Tensor<1,3> > &k_vectors,
+ const unsigned int N)
+ {
+ k_vectors.reinit(TableIndices<3>(N,N,N));
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ for (unsigned int k=0; k<N; ++k)
+ {
+ k_vectors(i,j,k)[0] = 2. * numbers::PI * i;
+ k_vectors(i,j,k)[0] = 2. * numbers::PI * j;
+ k_vectors(i,j,k)[0] = 2. * numbers::PI * k;
+ }
+ }
+
+
+ template <int dim>
+ Fourier<dim>::Fourier(const unsigned int N,
+ const hp::FECollection<dim> &fe_collection,
+ const hp::QCollection<dim> &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 <int dim>
+ void Fourier<dim>::calculate(const Vector<double> &local_dof_values,
+ const unsigned int cell_active_fe_index,
+ Table<dim,std::complex<double> > &fourier_coefficients)
+ {
+ ensure_existence(cell_active_fe_index);
+ const FullMatrix<std::complex<double> > &matrix = fourier_transform_matrices[cell_active_fe_index];
+
+ std::fill(unrolled_coefficients.begin(),
+ unrolled_coefficients.end(),
+ std::complex<double>(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 <int dim>
+ std::complex<double> integrate(const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const Tensor<1,dim> &k_vector,
+ const unsigned int j)
+ {
+ std::complex<double> sum = 0;
+ for (unsigned int q=0; q<quadrature.size(); ++q)
+ {
+ const Point<dim> &x_q = quadrature.point(q);
+ sum += std::exp(std::complex<double>(0,1) *
+ (k_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<k_vectors.size(0); ++k1)
+ for (unsigned int k2=0; k2<k_vectors.size(1); ++k2,k++)
+ for (unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
+ fourier_transform_matrices[fe](k,j) = integrate((*fe_collection)[fe],
+ (*q_collection)[fe],
+ k_vectors(k1,k2),
+ j);
+ }
+ }
+
+ template <>
+ 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<k_vectors.size(0); ++k1)
+ for (unsigned int k2=0; k2<k_vectors.size(1); ++k2)
+ for (unsigned int k3=0; k3<k_vectors.size(2); ++k3, k++)
+ for (unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
+ fourier_transform_matrices[fe](k,j) = integrate((*fe_collection)[fe],
+ (*q_collection)[fe],
+ k_vectors(k1,k2,k3),
+ j);
+ }
+ }
+
+ template <>
+ 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<k_vectors.size(0); ++k)
+ for (unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
+ fourier_transform_matrices[fe](k,j) = integrate((*fe_collection)[fe],
+ (*q_collection)[fe],
+ k_vectors(k),
+ j);
+ }
+ }
+
+
+ /*-------------- Legendre -------------------------------*/
+ DeclException2 (ExcLegendre, int, double,
+ << "x["<<arg1 << "] = "<<arg2 << " is not in [0,1]");
+
+ /* dim dimensional Legendre function with indices @p indices
+ * evaluated at @p x_q in [0,1]^dim.
+ */
+ template <int dim>
+ double Lh(const Point<dim> &x_q,
+ const TableIndices<dim> &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 <int dim>
+ double multiplier(const TableIndices<dim> &indices)
+ {
+ double res = 1.0;
+ for (unsigned int d = 0; d < dim; d++)
+ res *= (0.5+indices[d]);
+
+ return res;
+ }
+
+
+ template <int dim>
+ Legendre<dim>::Legendre(const unsigned int size_in_each_direction,
+ const hp::FECollection<dim> &fe_collection,
+ const hp::QCollection<dim> &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<dim>(N),
+ 0.)
+ {
+ }
+
+ template <int dim>
+ void Legendre<dim>::calculate(const dealii::Vector<double> &local_dof_values,
+ const unsigned int cell_active_fe_index,
+ Table<dim,double> &legendre_coefficients)
+ {
+ ensure_existence(cell_active_fe_index);
+ const FullMatrix<double> &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 <int dim>
+ double integrate_Legendre(const FiniteElement<dim> &fe,
+ const Quadrature<dim> &quadrature,
+ const TableIndices<dim> &indices,
+ const unsigned int dof)
+ {
+ double sum = 0;
+ for (unsigned int q=0; q<quadrature.size(); ++q)
+ {
+ const Point<dim> &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<N; ++k)
+ for (unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
+ legendre_transform_matrices[fe](k,j) =
+ integrate_Legendre((*fe_collection)[fe],
+ (*q_collection)[fe],
+ TableIndices<1>(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<N; ++k1)
+ for (unsigned int k2=0; k2<N; ++k2,k++)
+ for (unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
+ legendre_transform_matrices[fe](k,j) =
+ integrate_Legendre((*fe_collection)[fe],
+ (*q_collection)[fe],
+ TableIndices<2>(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<N; ++k1)
+ for (unsigned int k2=0; k2<N; ++k2)
+ for (unsigned int k3=0; k3<N; ++k3, k++)
+ for (unsigned int j=0; j<(*fe_collection)[fe].dofs_per_cell; ++j)
+ legendre_transform_matrices[fe](k,j) =
+ integrate_Legendre((*fe_collection)[fe],
+ (*q_collection)[fe],
+ TableIndices<3>(k1,k2,k3),
+ j);
+ }
+ }
+
+ /*-------------- linear_regression -------------------------------*/
+ std::pair<double,double> linear_regression(const std::vector<double> &x,
+ const std::vector<double> &y)
+ {
+ FullMatrix<double> K(2,2), invK(2,2);
+ Vector<double> 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <iostream>
+#include <fstream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/fe/fe_series.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_q.h>
+
+using namespace dealii;
+
+void test_1d ()
+{
+ const unsigned int N =4;
+ // exact values obtained by Maxima:
+ std::vector<std::complex<double> > exact(N);
+ const std::complex<double> I(0,1.);
+ exact[0] = std::complex<double>(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<dim> fe_collection;
+ hp::QCollection<dim> q_collection;
+
+ // linear FE
+ fe_collection.push_back(FE_Q<dim>(1));
+
+ QGauss<1> base_quadrature (6);
+ QIterated<dim> quadrature (base_quadrature, N);
+ q_collection.push_back(quadrature);
+
+ FESeries::Fourier<dim> fourier(N,fe_collection,q_collection);
+
+ Vector<double> local_dof_values(2);
+ local_dof_values[0] = 0;
+ local_dof_values[1] = 1.;
+ const unsigned int cell_active_fe_index =0;
+
+ Table<dim,std::complex<double> > 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();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <iostream>
+#include <fstream>
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/fe/fe_series.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <gsl/gsl_sf_legendre.h>
+
+using namespace dealii;
+
+/**
+ * A 1D function given by input legendre coefficients
+ */
+template<int dim>
+class LegendreFunction : public Function<dim>
+{
+public:
+ LegendreFunction(const std::vector<double> coefficients)
+ :
+ dealii::Function<dim>(1),
+ coefficients(coefficients)
+ {}
+
+ virtual double value(const dealii::Point<dim> &point,
+ const unsigned int component = 0 ) const;
+
+ const std::vector<double> &get_coefficients() const
+ {
+ return coefficients;
+ }
+
+private:
+
+ const std::vector<double> coefficients;
+};
+
+template<int dim>
+double LegendreFunction<dim>::value(const dealii::Point<dim> &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<int dim>
+void test(const LegendreFunction<dim> &func,
+ const unsigned int poly_degree)
+{
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler(triangulation);
+ hp::FECollection<dim> fe_collection;
+ fe_collection.push_back(dealii::FE_Q<dim>(poly_degree));
+
+ hp::QCollection<dim> quadrature_formula;
+ quadrature_formula.push_back(QGauss<dim>(poly_degree+6));
+
+ // reference cell:
+ GridGenerator::hyper_cube (triangulation,0.0,1.0);
+
+ dof_handler.distribute_dofs (fe_collection);
+
+ Vector<double> values(dof_handler.n_dofs());
+
+ VectorTools::interpolate (dof_handler,
+ func,
+ values);
+
+ const unsigned int N = poly_degree+1;
+ FESeries::Legendre<dim> legendre(N,
+ fe_collection,
+ quadrature_formula);
+
+ const std::vector<double> &coeff_in = func.get_coefficients();
+ Table<1,double> coeff_out(N);
+
+ Vector<double> local_dof_values;
+
+ typename hp::DoFHandler<dim>::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<double> coeff_in(2);
+ coeff_in[0] = 1.0;
+ coeff_in[1] = 2.0;
+ LegendreFunction<dim> function(coeff_in);
+ test(function,1);
+ }
+
+ {
+ std::vector<double> coeff_in(3);
+ coeff_in[0] = 1.0;
+ coeff_in[1] = 2.0;
+ coeff_in[2] = 3.0;
+ LegendreFunction<dim> function(coeff_in);
+ test(function,2);
+ }
+
+ {
+ std::vector<double> coeff_in(4);
+ coeff_in[0] = 1.0;
+ coeff_in[1] = 2.0;
+ coeff_in[2] = 3.0;
+ coeff_in[3] = 4.0;
+ LegendreFunction<dim> function(coeff_in);
+ test(function,3);
+ }
+
+ {
+ std::vector<double> 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<dim> function(coeff_in);
+ test(function,4);
+ }
+
+ {
+ std::vector<double> 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<dim> function(coeff_in);
+ test(function,5);
+ }
+
+ dealii::deallog << "Ok"<<std::endl;
+
+}
--- /dev/null
+
+DEAL::1.00000 2.00000
+DEAL::1.00000 2.00000
+DEAL::1.00000 2.00000 3.00000
+DEAL::1.00000 2.00000 3.00000
+DEAL::1.00000 2.00000 3.00000 4.00000
+DEAL::1.00000 2.00000 3.00000 4.00000
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000 6.00000
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000 6.00000
+DEAL::Ok
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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::process_coefficients()
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/fe/fe_series.h>
+
+using namespace dealii;
+
+std::pair<bool,unsigned int>
+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<dim,double> 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<unsigned int>,std::vector<double> > 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);
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <iostream>
+#include <fstream>
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/fe/fe_series.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <gsl/gsl_sf_legendre.h>
+
+using namespace dealii;
+
+template<int dim>
+class LegendreFunction : public Function<dim>
+{
+public:
+ LegendreFunction()
+ :
+ Function<dim>(1)
+ {}
+
+ virtual double value(const dealii::Point<dim> &point,
+ const unsigned int component = 0 ) const;
+};
+
+template<int dim>
+double LegendreFunction<dim>::value(const Point<dim> &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<int dim>
+void test(const LegendreFunction<dim> &func,
+ const unsigned int poly_degree)
+{
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler(triangulation);
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim> quadrature_formula;
+
+ for (unsigned int p = poly_degree; p<= poly_degree+3; p++)
+ {
+ fe_collection.push_back(dealii::FE_Q<dim>(p));
+ quadrature_formula.push_back(dealii::QGauss<dim>(p+1+5));
+ }
+
+ // reference cell
+ GridGenerator::hyper_cube (triangulation,0.0,1.0);
+
+ dof_handler.distribute_dofs (fe_collection);
+
+ Vector<double> values(dof_handler.n_dofs());
+
+ VectorTools::interpolate (dof_handler,func,values);
+ const unsigned int N = 4;
+ FESeries::Legendre<dim> legendre(N,fe_collection,quadrature_formula);
+
+ Table<1,double> coeff_out(N);
+ Vector<double> local_dof_values;
+
+ typename hp::DoFHandler<dim>::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 << " " <<local_dof_values[i];
+
+ dealii::deallog << std::endl;
+ }
+
+ deallog << "calculated:"<<std::endl;
+ for (unsigned int i = 0; i < N; i++)
+ deallog << coeff_out[i] << std::endl;
+
+ std::vector<double> 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:"<<std::endl;
+ for (unsigned int i = 0; i < coeff_exp.size(); i++)
+ deallog << coeff_exp[i] << std::endl;
+
+ dof_handler.clear();
+}
+
+
+/**
+ * Small test to first output Legendre coefficients from GSL at -1,0,1
+ * and then check that they are orthonormal
+ */
+void test_legendre_orthonormal(const unsigned int N)
+{
+ const unsigned int dim = 1;
+ deallog << "Pl @ -1;0;1"<<std::endl;
+ for (unsigned int l = 0; l < N; l++)
+ {
+ deallog << "l="<<l<<": ";
+ for (double x = -1.0; x <=1.0; x+=1.0)
+ deallog<< gsl_sf_legendre_Pl (l, x) << " ";
+
+ deallog<<std::endl;
+ }
+
+ QGauss<dim> 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<quadrature.size(); ++q)
+ {
+ const Point<dim> &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 << "("<<k1<<","<<k2<<") = " <<ortho<<std::endl;
+ }
+ deallog << std::endl;
+}
+
+int main ()
+{
+ const int dim = 1;
+
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-8);
+
+ test_legendre_orthonormal(3);
+ LegendreFunction<dim> function;
+ test(function,2);
+
+ dealii::deallog << "Ok"<<std::endl;
+}
--- /dev/null
+
+DEAL::Pl @ -1;0;1
+DEAL::l=0: 1.00000 1.00000 1.00000
+DEAL::l=1: -1.00000 0 1.00000
+DEAL::l=2: 1.00000 -0.500000 1.00000
+DEAL::orthogonality:
+DEAL::(0,0) = 1.00000
+DEAL::(0,1) = 0
+DEAL::(0,2) = 0
+DEAL::(1,0) = 0
+DEAL::(1,1) = 1.00000
+DEAL::(1,2) = 0
+DEAL::(2,0) = 0
+DEAL::(2,1) = 0
+DEAL::(2,2) = 1.00000
+DEAL::
+DEAL::local dofs: 1.81735e-05 0.000901649 1.35059e-05
+DEAL::calculated:
+DEAL::0.000114769
+DEAL::0.000312356
+DEAL::0.000210438
+DEAL::0
+DEAL::exact:
+DEAL::0.000114769
+DEAL::0.000312356
+DEAL::0.000210438
+DEAL::Ok
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 2D and 3D for a function given using Legendre
+// coefficients.
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/fe/fe_series.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <gsl/gsl_sf_legendre.h>
+
+using namespace dealii;
+
+template<int dim>
+class LegendreFunction : public Function<dim>
+{
+public:
+ LegendreFunction(const Table<dim,double> &coefficients)
+ :
+ dealii::Function<dim>(1),
+ coefficients(coefficients)
+ {
+ }
+
+ virtual double value(const Point<dim> &point,
+ const unsigned int component = 0 ) const;
+
+ const Table<dim,double> &get_coefficients() const
+ {
+ return coefficients;
+ }
+
+private:
+ const Table<dim,double> coefficients;
+};
+
+// copy-paste from fe_series.cc
+template <int dim>
+double Lh(const Point<dim> &x_q,
+ const TableIndices<dim> &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<int dim>
+void test(const LegendreFunction<dim> &func,
+ const unsigned int poly_degree)
+{
+ const unsigned int max_poly = poly_degree+3;
+ deallog <<"-----------------------------------"<<std::endl;
+ deallog << dim <<"d, p="<<poly_degree<<", max_p="<<max_poly<<std::endl;
+ deallog <<"-----------------------------------"<<std::endl;
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler(triangulation);
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim> quadrature_formula;
+
+ // add some extra FEs in fe_collection
+ for (unsigned int p = 1; p <= max_poly; p++)
+ {
+ fe_collection.push_back(FE_Q<dim>(p));
+ quadrature_formula.push_back(QGauss<dim>(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<double> values(dof_handler.n_dofs());
+
+ VectorTools::interpolate (dof_handler,func,values);
+
+ const unsigned int N = poly_degree+1;
+ FESeries::Legendre<dim> legendre(N,
+ fe_collection,
+ quadrature_formula);
+
+ const Table<dim,double> &coeff_in = func.get_coefficients();
+ Table<dim,double> coeff_out;
+ resize(coeff_out,N);
+
+ Vector<double> local_dof_values;
+
+ typename hp::DoFHandler<dim>::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:"<<std::endl;
+ print(coeff_in);
+
+ dof_handler.clear();
+}
+
+int main ()
+{
+ std::ofstream logfile("output");
+ dealii::deallog.attach(logfile,/*do not print job id*/false);
+ dealii::deallog.depth_console(0);
+ dealii::deallog.threshold_double(1e-8);
+
+ {
+ const unsigned int dim = 2;
+ const unsigned int coeff_1d = 2;
+ const unsigned int p = 1;
+ Table<dim,double> 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<dim> function(coeff_in);
+ test(function,p);
+ }
+
+ {
+ const unsigned int dim = 2;
+ const unsigned int coeff_1d = 3;
+ const unsigned int p = 2;
+ Table<dim,double> 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<dim> function(coeff_in);
+ test(function,p);
+ }
+
+ {
+ const unsigned int dim = 3;
+ const unsigned int coeff_1d = 2;
+ const unsigned int p = 1;
+ Table<dim,double> 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<dim> function(coeff_in);
+ test(function,p);
+ }
+
+ {
+ const unsigned int dim = 3;
+ const unsigned int coeff_1d = 3;
+ const unsigned int p = 2;
+ Table<dim,double> 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<dim> function(coeff_in);
+ test(function,p);
+ }
+
+ dealii::deallog << "Ok"<<std::endl;
+
+}
--- /dev/null
+DEAL::-----------------------------------
+DEAL::2d, p=1, max_p=4
+DEAL::-----------------------------------
+DEAL::calculated:
+DEAL::1.00000 2.00000 3.00000 4.00000
+DEAL::exact:
+DEAL::1.00000 2.00000 3.00000 4.00000
+DEAL::-----------------------------------
+DEAL::2d, p=2, max_p=5
+DEAL::-----------------------------------
+DEAL::calculated:
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 9.00000
+DEAL::exact:
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 9.00000
+DEAL::-----------------------------------
+DEAL::3d, p=1, max_p=4
+DEAL::-----------------------------------
+DEAL::calculated:
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000
+DEAL::exact:
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000
+DEAL::-----------------------------------
+DEAL::3d, p=2, max_p=5
+DEAL::-----------------------------------
+DEAL::calculated:
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 9.00000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000 27.0000
+DEAL::exact:
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 9.00000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000 27.0000
+DEAL::Ok
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 2D for a given vector of local DoF values.
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/fe/fe_series.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_q.h>
+
+using namespace dealii;
+
+void test_2d ()
+{
+ const unsigned int dim = 2;
+ const unsigned int N = 7;
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim> fourier_q_collection;
+
+ for (unsigned int degree=2; degree<=N; ++degree)
+ {
+ fe_collection.push_back (FE_Q<dim>(degree));
+ }
+
+ QGauss<1> base_quadrature (2);
+ QIterated<dim> quadrature (base_quadrature, N);
+ for (unsigned int i = 0; i < fe_collection.size(); i++)
+ fourier_q_collection.push_back(quadrature);
+
+ FESeries::Fourier<dim> fourier(N,fe_collection,fourier_q_collection);
+ Table<dim,std::complex<double> > fourier_coefficients;
+ fourier_coefficients.reinit(N,N);
+
+ Vector<double> 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();
+}
--- /dev/null
+
+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)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <iostream>
+#include <fstream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/fe/fe_series.h>
+
+using namespace dealii;
+
+std::pair<bool,unsigned int>
+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<dim,double> 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<unsigned int>,std::vector<double> > 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);
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <iostream>
+#include <fstream>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/fe/fe_series.h>
+
+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<double> 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<double,double> 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();
+}
--- /dev/null
+
+DEAL::exact: 3.14159 1.41421
+DEAL::calculated: 3.14159 1.41421
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_series.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+
+namespace Step27
+{
+ using namespace dealii;
+
+
+
+ template <int dim>
+ class LaplaceProblem
+ {
+ public:
+ LaplaceProblem ();
+ ~LaplaceProblem ();
+
+ void run ();
+
+ private:
+ void setup_system ();
+ void assemble_system ();
+ void solve ();
+ void create_coarse_grid ();
+ void estimate_smoothness (Vector<float> &smoothness_indicators);
+ void postprocess (const unsigned int cycle);
+
+ Triangulation<dim> triangulation;
+
+ hp::DoFHandler<dim> dof_handler;
+ hp::FECollection<dim> fe_collection;
+ hp::QCollection<dim> quadrature_collection;
+ hp::QCollection<dim-1> face_quadrature_collection;
+
+ hp::QCollection<dim> fourier_q_collection;
+ std_cxx11::shared_ptr<FESeries::Fourier<dim>> fourier;
+ std::vector<double> ln_k;
+ Table<dim,std::complex<double> > fourier_coefficients;
+
+ ConstraintMatrix constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ const unsigned int max_degree;
+ };
+
+
+
+ template <int dim>
+ class RightHandSide : public Function<dim>
+ {
+ public:
+ RightHandSide () : Function<dim> () {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const;
+ };
+
+
+ template <int dim>
+ double
+ RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+ {
+ double product = 1;
+ for (unsigned int d=0; d<dim; ++d)
+ product *= (p[d]+1);
+ return product;
+ }
+
+
+ template <typename T>
+ void resize(Table<2,T> &coeff, const unsigned int N)
+ {
+ coeff.reinit(N,N);
+ }
+
+
+ template <int dim>
+ LaplaceProblem<dim>::LaplaceProblem ()
+ :
+ dof_handler (triangulation),
+ max_degree (dim <= 2 ? 7 : 5)
+ {
+ for (unsigned int degree=2; degree<=max_degree; ++degree)
+ {
+ fe_collection.push_back (FE_Q<dim>(degree));
+ quadrature_collection.push_back (QGauss<dim>(degree+1));
+ face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+ }
+
+ const unsigned int N = max_degree;
+
+ QGauss<1> base_quadrature (2);
+ QIterated<dim> quadrature (base_quadrature, N);
+ for (unsigned int i = 0; i < fe_collection.size(); i++)
+ fourier_q_collection.push_back(quadrature);
+
+ fourier = std_cxx11::make_shared<FESeries::Fourier<dim> >(N,
+ fe_collection,
+ fourier_q_collection);
+ resize(fourier_coefficients,N);
+ }
+
+
+
+ template <int dim>
+ LaplaceProblem<dim>::~LaplaceProblem ()
+ {
+ dof_handler.clear ();
+ }
+
+
+ template <int dim>
+ void LaplaceProblem<dim>::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<dim>(),
+ 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 <int dim>
+ void LaplaceProblem<dim>::assemble_system ()
+ {
+ hp::FEValues<dim> hp_fe_values (fe_collection,
+ quadrature_collection,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const RightHandSide<dim> rhs_function;
+
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
+
+ std::vector<types::global_dof_index> local_dof_indices;
+
+ typename hp::DoFHandler<dim>::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<dim> &fe_values = hp_fe_values.get_present_fe_values ();
+
+ std::vector<double> 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_point<fe_values.n_quadrature_points;
+ ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+ cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ rhs_values[q_point] *
+ fe_values.JxW(q_point));
+ }
+
+ local_dof_indices.resize (dofs_per_cell);
+ cell->get_dof_indices (local_dof_indices);
+
+ constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
+ }
+ }
+
+
+
+
+ template <int dim>
+ void LaplaceProblem<dim>::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 <int dim>
+ void LaplaceProblem<dim>::postprocess (const unsigned int cycle)
+ {
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ face_quadrature_collection,
+ typename FunctionMap<dim>::type(),
+ solution,
+ estimated_error_per_cell);
+
+
+ Vector<float> smoothness_indicators (triangulation.n_active_cells());
+ estimate_smoothness (smoothness_indicators);
+
+ // Output to VTK
+ if (false)
+ {
+ Vector<float> fe_degrees (triangulation.n_active_cells());
+ {
+ typename hp::DoFHandler<dim>::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<dim,hp::DoFHandler<dim> > 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<dim>::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<dim>::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<Point<dim> > vertices (&vertices_1[0],
+ &vertices_1[n_vertices]);
+ static const int cell_vertices[][GeometryInfo<dim>::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<CellData<dim> > cells (n_cells, CellData<dim>());
+ for (unsigned int i=0; i<n_cells; ++i)
+ {
+ for (unsigned int j=0;
+ j<GeometryInfo<dim>::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 <int dim>
+ void LaplaceProblem<dim>::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 <int dim>
+ std::pair<bool,unsigned int>
+ predicate_ind(const TableIndices<dim> &ind);
+
+ template<>
+ std::pair<bool,unsigned int>
+ predicate_ind<2>(const TableIndices<2> &ind)
+ {
+ const unsigned int v = ind[0]*ind[0]+ind[1]*ind[1];
+ if (v>0 && v<7*7)
+ return std::make_pair(true,v);
+ else
+ return std::make_pair(false,v);
+ }
+
+ template <int dim>
+ void
+ LaplaceProblem<dim>::
+ estimate_smoothness (Vector<float> &smoothness_indicators)
+ {
+#ifdef OLD
+ const unsigned int N = max_degree;
+
+ std::vector<Tensor<1,dim> > k_vectors;
+ std::vector<unsigned int> k_vectors_magnitude;
+ switch (dim)
+ {
+ case 2:
+ {
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ if (!((i==0) && (j==0))
+ &&
+ (i*i + j*j < N*N))
+ {
+ k_vectors.push_back (Point<dim>(2.*numbers::PI * i,
+ 2.*numbers::PI * j));
+ k_vectors_magnitude.push_back (i*i+j*j);
+ }
+
+ break;
+ }
+
+ case 3:
+ {
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ for (unsigned int k=0; k<N; ++k)
+ if (!((i==0) && (j==0) && (k==0))
+ &&
+ (i*i + j*j + k*k < N*N))
+ {
+ k_vectors.push_back (Point<dim>(2.*numbers::PI * i,
+ 2.*numbers::PI * j,
+ 2.*numbers::PI * k));
+ k_vectors_magnitude.push_back (i*i+j*j+k*k);
+ }
+
+ break;
+ }
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ const unsigned n_fourier_modes = k_vectors.size();
+ std::vector<double> ln_k (n_fourier_modes);
+ for (unsigned int i=0; i<n_fourier_modes; ++i)
+ ln_k[i] = std::log (k_vectors[i].norm());
+
+
+ std::vector<Table<2,std::complex<double> > >
+ fourier_transform_matrices (fe_collection.size());
+
+ QGauss<1> base_quadrature (2);
+ QIterated<dim> quadrature (base_quadrature, N);
+
+
+ for (unsigned int fe=0; fe<fe_collection.size(); ++fe)
+ {
+ fourier_transform_matrices[fe].reinit (n_fourier_modes,
+ fe_collection[fe].dofs_per_cell);
+
+ for (unsigned int k=0; k<n_fourier_modes; ++k)
+ for (unsigned int j=0; j<fe_collection[fe].dofs_per_cell; ++j)
+ {
+ std::complex<double> sum = 0;
+ for (unsigned int q=0; q<quadrature.size(); ++q)
+ {
+ const Point<dim> x_q = quadrature.point(q);
+ sum += std::exp(std::complex<double>(0,1) *
+ (k_vectors[k] * x_q)) *
+ fe_collection[fe].shape_value(j,x_q) *
+ quadrature.weight(q);
+ }
+ fourier_transform_matrices[fe](k,j)
+ = sum;
+ }
+ }
+
+ std::vector<std::complex<double> > fourier_coefficients (n_fourier_modes);
+ Vector<double> local_dof_values;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ local_dof_values.reinit (cell->get_fe().dofs_per_cell);
+ cell->get_dof_values (solution, local_dof_values);
+
+ for (unsigned int f=0; f<n_fourier_modes; ++f)
+ {
+ fourier_coefficients[f] = 0;
+
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+ fourier_coefficients[f] +=
+ fourier_transform_matrices[cell->active_fe_index()](f,i)
+ *
+ local_dof_values(i);
+ }
+
+ std::map<unsigned int, double> k_to_max_U_map;
+ for (unsigned int f=0; f<n_fourier_modes; ++f)
+ if ((k_to_max_U_map.find (k_vectors_magnitude[f]) ==
+ k_to_max_U_map.end())
+ ||
+ (k_to_max_U_map[k_vectors_magnitude[f]] <
+ std::abs (fourier_coefficients[f])))
+ k_to_max_U_map[k_vectors_magnitude[f]]
+ = std::abs (fourier_coefficients[f]);
+
+ double sum_1 = 0,
+ sum_ln_k = 0,
+ sum_ln_k_square = 0,
+ sum_ln_U = 0,
+ sum_ln_U_ln_k = 0;
+ for (unsigned int f=0; f<n_fourier_modes; ++f)
+ if (k_to_max_U_map[k_vectors_magnitude[f]] ==
+ std::abs (fourier_coefficients[f]))
+ {
+ sum_1 += 1;
+ sum_ln_k += ln_k[f];
+ sum_ln_k_square += ln_k[f]*ln_k[f];
+ sum_ln_U += std::log (std::abs (fourier_coefficients[f]));
+ sum_ln_U_ln_k += std::log (std::abs (fourier_coefficients[f])) *
+ ln_k[f];
+ }
+
+ const double mu
+ = (1./(sum_1*sum_ln_k_square - sum_ln_k*sum_ln_k)
+ *
+ (sum_ln_k*sum_ln_U - sum_1*sum_ln_U_ln_k));
+
+ smoothness_indicators(cell->active_cell_index()) = mu - 1.*dim/2;
+ }
+#else
+ Vector<double> local_dof_values;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ local_dof_values.reinit (cell->get_fe().dofs_per_cell);
+ cell->get_dof_values (solution, local_dof_values);
+
+ fourier->calculate(local_dof_values,
+ cell->active_fe_index(),
+ fourier_coefficients);
+
+ std::pair<std::vector<unsigned int>, std::vector<double> > res =
+ FESeries::process_coefficients<dim>(fourier_coefficients,
+ predicate_ind<dim>,
+ VectorTools::Linfty_norm);
+
+ Assert (res.first.size() == res.second.size(),
+ ExcInternalError());
+
+ if (ln_k.size() == 0)
+ {
+ ln_k.resize(res.first.size(),0);
+ for (unsigned int f = 0; f < ln_k.size(); f++)
+ ln_k[f] = std::log (2.0*numbers::PI*std::sqrt(1.*res.first[f]));
+ }
+
+ for (unsigned int f = 0; f < res.second.size(); f++)
+ res.second[f] = std::log(res.second[f]);
+
+ std::pair<double,double> fit = FESeries::linear_regression(ln_k,res.second);
+ smoothness_indicators(cell->active_cell_index()) = -fit.first - 1.*dim/2;
+ }
+#endif
+ }
+}
+
+
+
+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;
+}
--- /dev/null
+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
projects / dealii.git / commitdiff