SmoothnessEstimator::Fourier Namespace Reference

Functions
template<int dim, int spacedim, typename VectorType >
void	coefficient_decay (FESeries::Fourier< dim, spacedim > &fe_fourier, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &solution, Vector< float > &smoothness_indicators, const VectorTools::NormType regression_strategy=VectorTools::Linfty_norm, const double smallest_abs_coefficient=1e-10, const bool only_flagged_cells=false)
template<int dim, int spacedim, typename VectorType >
void	coefficient_decay_per_direction (FESeries::Fourier< dim, spacedim > &fe_fourier, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &solution, Vector< float > &smoothness_indicators, const ComponentMask &coefficients_predicate=ComponentMask(), const double smallest_abs_coefficient=1e-10, const bool only_flagged_cells=false)
template<int dim, int spacedim>
FESeries::Fourier< dim, spacedim >	default_fe_series (const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)

Detailed Description

Smoothness estimation strategy based on the decay of Fourier expansion coefficients.

From the definition, we can write our Fourier series expansion \(a_{\bf k}\) of the finite element solution on cell \(K\) with polynomial degree \(p\) as a matrix product

\begin{eqnarray*} u_h({\bf x}) &=& \sum_j u_j \varphi_j ({\bf x}) \\ u_{h, {\bf k}}({\bf x}) &=& \sum_{{\bf k}, \|{\bf k}\|\le p} a_{\bf k} \phi_{\bf k}({\bf x}), \quad a_{\bf k} = \sum_j {\cal F}_{{\bf k},j} u_j \end{eqnarray*}

with \(u_j\) the degrees of freedom and \(\varphi_j\) the corresponding shape functions. \(\{\phi_{\bf k}({\bf x}) = \exp(i \, 2 \pi \, {\bf k} \cdot {\bf x}) \}\) are exponential functions on cell \(K\). \(a_{\bf k}\) and \({\cal F}_{{\bf k},j}\) are coefficients and transformation matrices from the Fourier expansion of each shape function. For practical reasons, we will perform the calculation of these matrices and coefficients only on the reference cell \(\hat K\). We only have to calculate the transformation matrices once this way. However, results are only applicable if mapping from the reference cell to the actual cell is linear. We use the class FESeries::Fourier to determine all coefficients \(a_{\bf k}\).

If the finite element approximation on cell \(K\) is part of the Hilbert space \(H^s(K)\), then the following integral must exist for both the finite element and spectral representation of our solution

\begin{eqnarray*} \| \nabla^s u_h({\bf x}) \|_{L^2(K)}^2 &=& \int\limits_K \left| \nabla^s u_h({\bf x}) \right|^2 d{\bf x} < \infty \\ \| \nabla^s u_{h, {\bf k}}({\bf x}) \|_{L^2(K)}^2 &=& \int\limits_K \left| \sum\limits_{\bf k} (-i \, 2 \pi {\bf k})^s \, a_{\bf k} \, \phi_{\bf k}({\bf x}) \right|^2 d{\bf x} = (2 \pi)^{2s} \sum\limits_{\bf k} \left| a_{\bf k} \right|^2 \|{\bf k}\|_2^{2s} < \infty \end{eqnarray*}

The sum is finite only if the summands decay at least with order

\[ |a_{\bf k}|^2 \|{\bf k}\|_2^{2s} \|{\bf k}\|_2^{d - 1} = {\cal O}\left( \|{\bf k}\|_2^{-1-\epsilon} \right) \]

for all \(\epsilon > 0\). The additional factor stems from the fact that, since we sum over all multi-indices \({\bf k}\) that are located on a dim-dimensional sphere, we actually have, up to a constant, \(\|{\bf k}\|_2^{d-1}\) modes located in each increment \(\|{\bf k}\|_2 + d\|{\bf k}\|_2\) that need to be taken into account. By a comparison of exponents, we can rewrite this condition as

\[ |a_{\bf k}| = {\cal O}\left(\|{\bf k}\|_2^ {-\left(s + \frac d2 + \epsilon \right)} \right) \]

The next step is to estimate how fast these coefficients decay with \(\|{\bf k}\|_2\). Thus, we perform a least-squares fit

\[ \min_{\alpha,\sigma} \frac 12 \sum_{{\bf k}, \|{\bf k}\|_2 \le p} \left( |a_{\bf k}| - \alpha \|{\bf k}\|_2^{-\sigma}\right)^2 \]

with regression coefficients \(\alpha\) and \(\sigma\). For simplification, we apply a logarithm on our minimization problem

\[ \min_{\beta,\sigma} Q(\beta,\sigma) = \frac 12 \sum_{{\bf k}, \|{\bf k}\|_2 \le p} \left( \ln |a_{\bf k}| - \beta + \sigma \ln \|{\bf k}\|_2 \right)^2, \]

where \(\beta=\ln \alpha\). This is now a problem for which the optimality conditions \(\frac{\partial Q}{\partial\beta}=0, \frac{\partial Q}{\partial\sigma}=0\), are linear in \(\beta,\sigma\). We can write these conditions as follows:

\[ \left(\begin{array}{cc} \sum_{{\bf k}, \|{\bf k}\|_2 \le p} 1 & \sum_{{\bf k}, \|{\bf k}\|_2 \le p} \ln \|{\bf k}\|_2 \\ \sum_{{\bf k}, \|{\bf k}\|_2 \le p} \ln \|{\bf k}\|_2 & \sum_{{\bf k}, \|{\bf k}\|_2 \le p} (\ln \|{\bf k}\|_2)^2 \end{array}\right) \left(\begin{array}{c} \beta \\ -\sigma \end{array}\right) = \left(\begin{array}{c} \sum_{{\bf k}, \|{\bf k}\|_2\le p} \ln |a_{{\bf k}}| \\ \sum_{{\bf k}, \|{\bf k}\|_2\le p} \ln |a_{{\bf k}}| \ln \|{\bf k}\|_2 \end{array}\right) \]

Solving for \(\beta\) and \(\sigma\) is just a linear regression fit and to do that we will use FESeries::linear_regression().

While we are not particularly interested in the actual value of \(\beta\), the formula above gives us a means to calculate the value of the exponent \(\sigma\) that we can then use to determine that \(u(\hat{\bf x})\) is in \(H^s(K)\) with \(s=\sigma-\frac d2\). The decay rates \(\sigma\) will suffice as our smoothness indicators and will be calculated on each cell for any provided solution.

Note

An extensive demonstration of the use of these functions is provided in step-27.

Function Documentation

coefficient_decay()

template<int dim, int spacedim, typename VectorType >

void SmoothnessEstimator::Fourier::coefficient_decay	(	FESeries::Fourier< dim, spacedim > &	fe_fourier,
		const DoFHandler< dim, spacedim > &	dof_handler,
		const VectorType &	solution,
		Vector< float > &	smoothness_indicators,
		const VectorTools::NormType	regression_strategy = VectorTools::Linfty_norm,
		const double	smallest_abs_coefficient = 1e-10,
		const bool	only_flagged_cells = false
	)

In this variant of the estimation strategy for higher dimensions, we will consider all mode vectors \(\bf k\) describing Fourier polynomials \(P_{\bf k}\) and perform one least-squares fit over all coefficients at once. If there are multiple coefficients corresponding to the same absolute value of modes \(\|\bf k\|_2\), we take the maximum among those. Thus, the least-squares fit is performed on

\[ \ln \left( \max\limits_{\|{\bf k}\|_2} |a_{\bf k}| \right) \sim C - \sigma \ln \|{\bf k}\|_2 \]

for \({\bf k}=(k_1,\ldots,k_d)\) and \(k_i=0,\ldots,p\), with \(p\) the polynomial degree of the finite element. We exclude the \(\|{\bf k}\|_2=0\) modes to avoid the singularity of the logarithm.

The regression_strategy parameter determines which norm will be used on the subset of coefficients \(\bf k\) with the same absolute value \(\|{\bf k}\|_2\). Default is VectorTools::Linfty_norm for a maximum approximation.

For a provided solution vector solution defined on a DoFHandler dof_handler, this function returns a vector smoothness_indicators with as many elements as there are cells where each element contains the estimated regularity \(\sigma\).

A series expansion object fe_fourier has to be supplied, which needs to be constructed with the same FECollection object as the dof_handler.

The parameter smallest_abs_coefficient allows to ignore small (in absolute value) coefficients within the linear regression fit. In case there are less than two nonzero coefficients for a coordinate direction, this direction will be skipped. If all coefficients are zero, the returned value for this cell will be \(\sigma=\infty\).

Smoothness indicators are usually used to decide whether to perform h- or p-adaptation. So in most cases, we only need to calculate those indicators on cells flagged for refinement or coarsening. The parameter only_flagged_cells controls whether this particular subset or all cells will be considered. By default, all cells will be considered. When only flagged cells are supposed to be considered, smoothness indicators will only be set on those vector entries of flagged cells; the others will be set to a signaling NaN.

Definition at line 369 of file smoothness_estimator.cc.

coefficient_decay_per_direction()

template<int dim, int spacedim, typename VectorType >

void SmoothnessEstimator::Fourier::coefficient_decay_per_direction	(	FESeries::Fourier< dim, spacedim > &	fe_fourier,
		const DoFHandler< dim, spacedim > &	dof_handler,
		const VectorType &	solution,
		Vector< float > &	smoothness_indicators,
		const ComponentMask &	coefficients_predicate = ComponentMask(),
		const double	smallest_abs_coefficient = 1e-10,
		const bool	only_flagged_cells = false
	)

In this variant of the estimation strategy for higher dimensions, we only consider modes in each coordinate direction, i.e., only mode vectors \(\bf k\) with one nonzero entry. We perform the least-squares fit in each coordinate direction separately and take the lowest decay rate \(\sigma\) among them.

The coefficients_predicate parameter selects Fourier coefficients \(a_j\), \(j=0,\ldots,p\) for linear regression in each coordinate direction. The user is responsible for updating the vector of flags provided to this function. Note that its size is \(p+1\), where \(p\) is the polynomial degree of the FE basis on a given element. The default implementation will use all Fourier coefficients in each coordinate direction, i.e., set all the elements of the vector to true.

For a provided solution vector solution defined on a DoFHandler dof_handler, this function returns a vector smoothness_indicators with as many elements as there are cells where each element contains the estimated regularity \(\sigma\).

A series expansion object fe_fourier has to be supplied, which needs to be constructed with the same FECollection object as the dof_handler.

The parameter smallest_abs_coefficient allows to ignore small (in absolute value) coefficients within the linear regression fit. In case there are fewer than two nonzero coefficients for a coordinate direction, this direction will be skipped. If all coefficients are zero, the returned value for this cell will be \(\sigma=\infty\).

Smoothness indicators are usually used to decide whether to perform h- or p-adaptation. So in most cases, we only need to calculate those indicators on cells flagged for refinement or coarsening. The parameter only_flagged_cells controls whether this particular subset or all cells will be considered. By default, all cells will be considered. When only flagged cells are supposed to be considered, smoothness indicators will only be set on those vector entries of flagged cells; the others will be set to a signaling NaN.

Definition at line 467 of file smoothness_estimator.cc.

default_fe_series()

template<int dim, int spacedim>

FESeries::Fourier< dim, spacedim > SmoothnessEstimator::Fourier::default_fe_series	(	const hp::FECollection< dim, spacedim > &	fe_collection,
		const unsigned int	component = numbers::invalid_unsigned_int
	)

Returns a FESeries::Fourier object for Fourier series expansions with the default configuration for smoothness estimation purposes.

For each finite element of the provided fe_collection, we use as many modes as its polynomial degree plus two. Further for each element, we use a 5-point Gaussian quarature iterated in each dimension by the maximal wave number, which is the number of modes decreased by one since we start with \(k = 0\).

As the Fourier expansion can only be performed on scalar fields, this class does not operate on vector-valued finite elements and will therefore throw an assertion. However, each component of a finite element field can be treated as a scalar field, respectively, on which Fourier expansions are again possible. For this purpose, the optional parameter component defines which component of each FiniteElement will be used. The default value of component only applies to scalar FEs, in which case it indicates that the sole component is to be decomposed. For vector-valued FEs, a non-default value must be explicitly provided.

Definition at line 576 of file smoothness_estimator.cc.