@DEAL_II_EXPAND_TRILINOS_BLOCK_SPARSITY_PATTERN@;
}
-// Series expansion templates
-SERIES_EXPANSION_TEMPLATES := { FESeries::Fourier;
- FESeries::Legendre;
- }
-
// all supported logical dimensions
DIMENSIONS := { 1; 2; 3 }
--- /dev/null
+New: Add FESeries::Legendre::get_size_in_each_direction() to retrieve
+the number of coefficients in each direction. Also add an Assert in
+FESeries::Legendre::calculate() to check the dimension of the table to store
+coefficients.
+<br>
+(Denis Davydov, 2018/12/27)
These steps outlined above are applicable on many different scenarios, which
motivated the introduction of a generic function
-SmoothnessEstimator::estimate_by_coeff_decay() in deal.II, that combines all
+SmoothnessEstimator::estimate_by_coefficient_decay() in deal.II, that combines all
the tasks described in this section in one simple function call.
<h4>Compensating for anisotropy</h4>
// Estimating the smoothness is performed with the method of decaing
// expansion coefficients as outlined in the introduction.
Vector<float> smoothness_indicators;
- SmoothnessEstimator::estimate_by_coeff_decay<FESeries::Fourier<dim>>(
- dof_handler, solution, smoothness_indicators);
+ SmoothnessEstimator::fourier_coefficient_decay(dof_handler,
+ solution,
+ smoothness_indicators);
// Next we want to generate graphical output. In addition to the two
// estimated quantities derived above, we would also like to output the
const unsigned int cell_active_fe_index,
Table<dim, CoefficientType> & legendre_coefficients);
+ /**
+ * Return number of coefficients in each coordinate direction.
+ */
+ unsigned int
+ get_size_in_each_direction() const;
+
/**
* Calculate all transformation matrices to transfer the finite element
* solution to the series expansion representation.
#include <deal.II/base/config.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_series.h>
+
+#include <deal.II/hp/dof_handler.h>
+
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/vector_tools.h>
/**
- * Estimate the smoothness of a solution based on the decay of coefficients from
- * a series expansion.
- *
- * From the definition, we can write our series expansion $\hat U_{\bf k}$ as a
- * matrix product
- * @f[
- * \hat U_{\bf k}
- * = {\cal F}_{{\bf k},j} u_j,
- * @f]
- * with $u_j$ the coefficients and ${\cal F}_{{\bf k},j}$ the transformation
- * matrix from the expansion. We use the classes FESeries::Fourier and
- * FESeries::Legendre to determine all coefficients $u_j$.
- *
- * The next step is that we have to estimate how fast these coefficients
- * decay with $|{\bf k}|$. Thus, we perform a least-squares fit
- * @f[
- * \min_{\alpha,\mu}
- * \frac 12 \sum_{{\bf k}, |{\bf k}|\le N}
- * \left( |\hat U_{\bf k}| - \alpha |{\bf k}|^{-\mu}\right)^2
- * @f]
- * with linear regressions coefficients $\alpha$ and $\mu$. For simplification,
- * we apply a logarithm on our minimization problem
- * @f[
- * \min_{\beta,\mu}
- * Q(\beta,\mu) =
- * \frac 12 \sum_{{\bf k}, |{\bf k}|\le N}
- * \left( \ln |\hat U_{\bf k}| - \beta + \mu \ln |{\bf k}|\right)^2,
- * @f]
- * where $\beta=\ln \alpha$. This is now a problem for which the
- * optimality conditions $\frac{\partial Q}{\partial\beta}=0,
- * \frac{\partial Q}{\partial\mu}=0$, are linear in $\beta,\mu$. We can
- * write these conditions as follows:
- * @f[
- * \left(\begin{array}{cc}
- * \sum_{{\bf k}, |{\bf k}|\le N} 1 &
- * \sum_{{\bf k}, |{\bf k}|\le N} \ln |{\bf k}|
- * \\
- * \sum_{{\bf k}, |{\bf k}|\le N} \ln |{\bf k}| &
- * \sum_{{\bf k}, |{\bf k}|\le N} (\ln |{\bf k}|)^2
- * \end{array}\right)
- * \left(\begin{array}{c}
- * \beta \\ -\mu
- * \end{array}\right)
- * =
- * \left(\begin{array}{c}
- * \sum_{{\bf k}, |{\bf k}|\le N} \ln |\hat U_{{\bf k}}|
- * \\
- * \sum_{{\bf k}, |{\bf k}|\le N} \ln |\hat U_{{\bf k}}| \ln |{\bf k}|
- * \end{array}\right)
- * @f]
- * Solving for $\beta$ and $\mu$ is nothing else but a linear regression fit and
- * to do that we will use FESeries::linear_regression().
- *
- * While we are not particularly interested in the actual value of
- * $\beta$, the formula above gives us a mean to calculate the value of
- * the exponent $\mu$ that we can then use to determine that
- * $\hat u(\hat{\bf x})$ is in $H^s(\hat K)$ with $s=\mu-\frac d2$. These
- * Sobolev indices $s$ will suffice as our smoothness estimators and will be
- * calculated on each cell for any provided solution.
- *
- * @note An extensive demonstration of the use of these functions is provided in step-27.
- *
- * @ingroup numerics
- * @author Denis Davydov, 2016, Marc Fehling, 2018
+ * A namespace for various smoothness estimation strategies for hp-adaptive FEM.
*/
namespace SmoothnessEstimator
{
/**
- * Estimates the smoothness of the provided solutions using the method of
- * decaying coefficents as outlined above.
+ * Estimate smoothness from decay of Legendre absolute values of coefficients
+ * on the reference cell.
+ *
+ * In one dimension, the finite element solution on the reference element with
+ * polynomial degree $p$ can be written as
+ * @f[
+ * u_h(\hat x) = \sum_{j=0}^{p} a_j P_j(\hat x)
+ * @f]
+ * where $\{P_j(x)\}$ are Legendre polynomials. The decay of the coefficients
+ * is estimated by performing the linear regression fit of
+ * @f[
+ * \ln |a_j| \sim C - \sigma j
+ * @f]
+ * for $j=0,..,p$. The rate of the decay $\sigma$ can be used to estimate the
+ * smoothness. For example, one strategy to implement hp-refinement
+ * criteria is to perform p-refinement if $\sigma>1$.
+ *
+ * Extension to higher dimension is done by performing the fit in each
+ * coordinate direction separately and then taking the lowest value of
+ * $\sigma$.
+ *
+ * For each input vector of degrees of freedom defined on a DoFHandler,
+ * this function returns a vector with as many elements as there are cells
+ * where each element contains $\exp(-\sigma)$, which is a so-called
+ * analyticity (see references below).
+ *
+ * @param [in] fe_series FESeries::Legendre object to calculate coefficients.
+ * This object needs to be initialized to have at least $p+1$ coefficients in
+ * each direction, where $p$ is the maximum polynomial degree to be used.
+ * @param [in] dof_hander An hp::DoFHandler
+ * @param [in] all_solutions A vector of pointers to the solution vectors
+ * @param [out] all_smoothness_indicators A vector of pointers to the smoothness indicators for each @p all_solutions.
+ * @param [in] coefficients_predicate A predicate to select Legendre
+ * coefficients $a_j \;\; j=0\dots 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. Default
+ * implementation will use all Legendre coefficients in each coordinate
+ * direction, i.e. set all elements of the vector to `true`.
+ * @param [in] smallest_abs_coefficient The smallest absolute value of the
+ * coefficient to be used in linear regression in each coordinate direction.
+ * Note that Legendre coefficients of some functions may have a repeating
+ * pattern of zero coefficients (i.e. for functions that are locally symmetric
+ * or antisymmetric about the midpoint of the element in any coordinate
+ * direction). Thus this parameters allows to ingore small (in absolute value)
+ * coefficients within the linear regression fit. In case there are less than
+ * two non-zero coefficients for a coordinate direction, this direction will
+ * be skipped. If all coefficients are zero, the returned value for this cell
+ * will be zero (i.e. corresponding to the $\sigma=\infty$).
+ *
+ * For more theoretical details see
+ * @code{.bib}
+ * @Article{Mavriplis1994,
+ * author = {Mavriplis, Catherine},
+ * title = {Adaptive mesh strategies for the spectral element method},
+ * journal = {{Computer Methods in Applied Mechanics and Engineering}},
+ * year = {1994},
+ * volume = {116},
+ * number = {1},
+ * pages = {77--86},
+ * publisher = {Elsevier},
+ * }
+ * @article{Houston2005,
+ * author = {Houston, Paul and S{\"u}li, Endre},
+ * title = {A note on the design of hp-adaptive finite element
+ * methods for elliptic partial differential equations},
+ * journal = {{Computer Methods in Applied Mechanics and Engineering}},
+ * number = {2},
+ * pages = {229--243},
+ * publisher = {Elsevier},
+ * volume = {194},
+ * year = {2005}
+ * }
+ * @article{Eibner2007,
+ * author = {Eibner, Tino and Melenk, Jens Markus},
+ * title = {An adaptive strategy for hp-FEM based on testing for
+ * analyticity},
+ * journal = {{Computational Mechanics}},
+ * year = {2007},
+ * volume = {39},
+ * number = {5},
+ * pages = {575--595},
+ * publisher = {Springer},
+ * }
+ * @endcode
+ * and for the application within the deal.II:
+ * @code{.bib}
+ * @article{Davydov2017,
+ * author = {Denis Davydov and Tymofiy Gerasimov and Jean-Paul Pelteret and
+ * Paul Steinmann},
+ * title = {Convergence study of the h-adaptive PUM and the hp-adaptive
+ * FEM applied to eigenvalue problems in quantum mechanics},
+ * journal = {{Advanced Modeling and Simulation in Engineering Sciences}},
+ * year = {2017},
+ * volume = {4},
+ * number = {1},
+ * pages = {7},
+ * issn = {2213-7467},
+ * doi = {10.1186/s40323-017-0093-0},
+ * }
+ * @endcode
+ *
+ * @ingroup numerics
+ * @author Denis Davydov, 2018
+ */
+ template <int dim, int spacedim, typename VectorType>
+ void
+ legendre_coefficient_decay(
+ FESeries::Legendre<dim, spacedim> & fe_series,
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
+ const std::vector<const VectorType *> &all_solutions,
+ const std::vector<Vector<float> *> & all_smoothness_indicators,
+ const std::function<void(std::vector<bool> &flags)> coefficients_predicate =
+ [](std::vector<bool> &flags) -> void {
+ std::fill(flags.begin(), flags.end(), true);
+ },
+ const double smallest_abs_coefficient = 1e-10);
+
+ /**
+ * Same as above, but for a single solution vector.
+ */
+ template <int dim, int spacedim, typename VectorType>
+ void
+ legendre_coefficient_decay(
+ FESeries::Legendre<dim, spacedim> & fe_series,
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
+ const VectorType & solution,
+ Vector<float> & smoothness_indicators,
+ const std::function<void(std::vector<bool> &flags)> coefficients_predicate =
+ [](std::vector<bool> &flags) -> void {
+ std::fill(flags.begin(), flags.end(), true);
+ },
+ const double smallest_abs_coefficient = 1e-10);
+
+ /**
+ * Same as above, but for a single solution vector and with the default
+ * FESeries::Legendre.
+ */
+ template <int dim, int spacedim, typename VectorType>
+ void
+ legendre_coefficient_decay(
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
+ const VectorType & solution,
+ Vector<float> & smoothness_indicators,
+ const std::function<void(std::vector<bool> &flags)> coefficients_predicate =
+ [](std::vector<bool> &flags) -> void {
+ std::fill(flags.begin(), flags.end(), true);
+ },
+ const double smallest_abs_coefficient = 1e-10);
+
+ /**
+ * Estimate the smoothness of a solution based on the decay of coefficients
+ * from a series expansion.
+ *
+ * From the definition, we can write our series expansion $\hat U_{\bf k}$ as
+ * a matrix product
+ * @f[
+ * \hat U_{\bf k}
+ * = {\cal F}_{{\bf k},j} u_j,
+ * @f]
+ * with $u_j$ the coefficients and ${\cal F}_{{\bf k},j}$ the transformation
+ * matrix from the expansion. We use the classes FESeries::Fourier and
+ * FESeries::Legendre to determine all coefficients $u_j$.
+ *
+ * The next step is that we have to estimate how fast these coefficients
+ * decay with $|{\bf k}|$. Thus, we perform a least-squares fit
+ * @f[
+ * \min_{\alpha,\mu}
+ * \frac 12 \sum_{{\bf k}, |{\bf k}|\le N}
+ * \left( |\hat U_{\bf k}| - \alpha |{\bf k}|^{-\mu}\right)^2
+ * @f]
+ * with linear regressions coefficients $\alpha$ and $\mu$. For
+ * simplification, we apply a logarithm on our minimization problem
+ * @f[
+ * \min_{\beta,\mu}
+ * Q(\beta,\mu) =
+ * \frac 12 \sum_{{\bf k}, |{\bf k}|\le N}
+ * \left( \ln |\hat U_{\bf k}| - \beta + \mu \ln |{\bf k}|\right)^2,
+ * @f]
+ * where $\beta=\ln \alpha$. This is now a problem for which the
+ * optimality conditions $\frac{\partial Q}{\partial\beta}=0,
+ * \frac{\partial Q}{\partial\mu}=0$, are linear in $\beta,\mu$. We can
+ * write these conditions as follows:
+ * @f[
+ * \left(\begin{array}{cc}
+ * \sum_{{\bf k}, |{\bf k}|\le N} 1 &
+ * \sum_{{\bf k}, |{\bf k}|\le N} \ln |{\bf k}|
+ * \\
+ * \sum_{{\bf k}, |{\bf k}|\le N} \ln |{\bf k}| &
+ * \sum_{{\bf k}, |{\bf k}|\le N} (\ln |{\bf k}|)^2
+ * \end{array}\right)
+ * \left(\begin{array}{c}
+ * \beta \\ -\mu
+ * \end{array}\right)
+ * =
+ * \left(\begin{array}{c}
+ * \sum_{{\bf k}, |{\bf k}|\le N} \ln |\hat U_{{\bf k}}|
+ * \\
+ * \sum_{{\bf k}, |{\bf k}|\le N} \ln |\hat U_{{\bf k}}| \ln |{\bf k}|
+ * \end{array}\right)
+ * @f]
+ * Solving for $\beta$ and $\mu$ is nothing else but a linear regression fit
+ * and to do that we will use FESeries::linear_regression().
+ *
+ * While we are not particularly interested in the actual value of
+ * $\beta$, the formula above gives us a mean to calculate the value of
+ * the exponent $\mu$ that we can then use to determine that
+ * $\hat u(\hat{\bf x})$ is in $H^s(\hat K)$ with $s=\mu-\frac d2$. These
+ * Sobolev indices $s$ will suffice as our smoothness estimators and will be
+ * calculated on each cell for any provided solution.
+ *
+ * @note An extensive demonstration of the use of these functions is provided in step-27.
*
* The @p regression_strategy parameter determines which norm will be used on the subset of
* coeffiecients $\mathbf{k}$ with the same absolute value $|\mathbf{k}|$.
*
* An individual @p fe_series object can be supplied, which has to be constructed with the
* same FECollection object as the @p dof_handler.
+ *
+ * @ingroup numerics
+ * @author Denis Davydov, 2016, Marc Fehling, 2018
+ *
*/
- template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+ template <int dim, int spacedim, typename VectorType>
void
- estimate_by_coeff_decay(
- FESeriesType & fe_series,
- const DoFHandlerType & dof_handler,
+ fourier_coefficient_decay(
+ FESeries::Fourier<dim, spacedim> & fe_series,
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
const std::vector<const VectorType *> &all_solutions,
const std::vector<Vector<float> *> & all_smoothness_indicators,
const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm);
/**
* Same as the function above, only for one @p solution vector.
*/
- template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+ template <int dim, int spacedim, typename VectorType>
void
- estimate_by_coeff_decay(
- FESeriesType & fe_series,
- const DoFHandlerType & dof_handler,
- const VectorType & solution,
- Vector<float> & smoothness_indicators,
+ fourier_coefficient_decay(
+ FESeries::Fourier<dim, spacedim> & fe_series,
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const VectorType & solution,
+ Vector<float> & smoothness_indicators,
const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm);
/**
*
* Provide the desired series expansion as a template argument, i.e.
* @code
- * SmoothnessEstimator::estimate_by_coeff_decay<FESeries::Fourier<dim>>(
+ * SmoothnessEstimator::estimate_by_coefficient_decay<FESeries::Fourier<dim>>(
* dof_handler, all_solutions, all_smoothness_indicators);
* @endcode
*/
- template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+ template <int dim, int spacedim, typename VectorType>
void
- estimate_by_coeff_decay(
- const DoFHandlerType & dof_handler,
+ fourier_coefficient_decay(
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
const std::vector<const VectorType *> &all_solutions,
const std::vector<Vector<float> *> & all_smoothness_indicators,
const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm);
/**
* Same as the function above, only for one @p solution vector.
*/
- template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+ template <int dim, int spacedim, typename VectorType>
void
- estimate_by_coeff_decay(
- const DoFHandlerType & dof_handler,
- const VectorType & solution,
- Vector<float> & smoothness_indicators,
+ fourier_coefficient_decay(
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const VectorType & solution,
+ Vector<float> & smoothness_indicators,
const VectorTools::NormType regression_strategy = VectorTools::Linfty_norm);
} // namespace SmoothnessEstimator
/**
- * Calculates predicates of @p ind in the form
- * \f$
- * v = \sum\limits_{d=0}^{dim} ind[d]^2
- * \f$.
- *
- * We flag the predicate whether it fulfills the criterion
- * \f$
- * 0 < v < max_degree^2
- * \f$
- * using @p max_degree.
+ * 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>
} // namespace
+ template <int dim, int spacedim, typename VectorType>
+ void
+ legendre_coefficient_decay(
+ FESeries::Legendre<dim, spacedim> & fe_legendre,
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
+ const std::vector<const VectorType *> &all_solutions,
+ const std::vector<Vector<float> *> & all_smoothness_indicators,
+ const std::function<void(std::vector<bool> &flags)> coefficients_predicate,
+ const double smallest_abs_coefficient)
+ {
+ Assert(smallest_abs_coefficient >= 0.,
+ ExcMessage("smallest_abs_coefficient should be non-negative."));
+
+ using number = typename VectorType::value_type;
+ using number_coeff =
+ typename FESeries::Legendre<dim, spacedim>::CoefficientType;
+
+ AssertDimension(all_solutions.size(), all_smoothness_indicators.size());
+
+ for (auto &smoothness_indicator : all_smoothness_indicators)
+ smoothness_indicator->reinit(
+ dof_handler.get_triangulation().n_active_cells());
+
+ Table<dim, number_coeff> expansion_coefficients;
+ resize(expansion_coefficients, fe_legendre.get_size_in_each_direction());
+
+ Vector<number> local_dof_values;
+
+ // auxiliary vector to do linear regression
+ std::vector<number_coeff> x;
+ std::vector<number_coeff> y;
+
+ x.reserve(dof_handler.get_fe_collection().max_degree());
+ y.reserve(dof_handler.get_fe_collection().max_degree());
+
+ // precalculate predicates for each degree:
+ std::vector<std::vector<bool>> predicates(
+ dof_handler.get_fe_collection().max_degree());
+ for (unsigned int p = 1; p <= dof_handler.get_fe_collection().max_degree();
+ ++p)
+ {
+ auto &pred = predicates[p - 1];
+ // we have p+1 coefficients for degree p
+ pred.resize(p + 1);
+ coefficients_predicate(pred);
+ }
+
+ for (auto &cell : dof_handler.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ local_dof_values.reinit(cell->get_fe().dofs_per_cell);
+
+ const unsigned int pe = cell->get_fe().degree;
+
+ Assert(pe > 0, ExcInternalError());
+ const auto &pred = predicates[pe - 1];
+
+ // since we use coefficients with indices [1,pe] in each direction,
+ // the number of coefficients we need to calculate is at least N=pe+1
+ AssertIndexRange(pe, fe_legendre.get_size_in_each_direction());
+
+ auto solution_it = all_solutions.cbegin();
+ auto smoothness_indicators_it = all_smoothness_indicators.begin();
+ for (; solution_it != all_solutions.cend();
+ ++solution_it, ++smoothness_indicators_it)
+ {
+ cell->get_dof_values(*(*solution_it), local_dof_values);
+ fe_legendre.calculate(local_dof_values,
+ cell->active_fe_index(),
+ expansion_coefficients);
+
+ // choose the smallest decay of coefficients in each direction,
+ // i.e. the maximum decay slope k_v
+ number_coeff k_v = -std::numeric_limits<number_coeff>::max();
+ for (unsigned int d = 0; d < dim; d++)
+ {
+ x.resize(0);
+ y.resize(0);
+
+ // will use all non-zero coefficients allowed by the predicate
+ // function
+ Assert(pred.size() == pe + 1, ExcInternalError());
+ for (unsigned int i = 0; i <= pe; i++)
+ if (pred[i])
+ {
+ TableIndices<dim> ind;
+ ind[d] = i;
+ const number_coeff coeff_abs =
+ std::abs(expansion_coefficients(ind));
+
+ if (coeff_abs > smallest_abs_coefficient)
+ {
+ y.push_back(std::log(coeff_abs));
+ x.push_back(i);
+ }
+ }
+
+ // in case we don't have enough non-zero coefficient to fit,
+ // skip this direction
+ if (x.size() < 2)
+ continue;
+
+ const std::pair<number_coeff, number_coeff> fit =
+ FESeries::linear_regression(x, y);
+
+ // decay corresponds to negative slope
+ // take the lesser negative slope along each direction
+ k_v = std::max(k_v, fit.first);
+ }
+
+ (*(*smoothness_indicators_it))(cell->active_cell_index()) =
+ std::exp(k_v);
+ }
+ }
+ }
- template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+
+
+ template <int dim, int spacedim, typename VectorType>
void
- estimate_by_coeff_decay(
- FESeriesType & fe_series,
- const DoFHandlerType & dof_handler,
+ legendre_coefficient_decay(
+ FESeries::Legendre<dim, spacedim> & fe_legendre,
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
+ const VectorType & solution,
+ Vector<float> & smoothness_indicators,
+ const std::function<void(std::vector<bool> &flags)> coefficients_predicate,
+ const double smallest_abs_coefficient)
+ {
+ const std::vector<const VectorType *> all_solutions(1, &solution);
+ const std::vector<Vector<float> *> all_smoothness_indicators(
+ 1, &smoothness_indicators);
+
+ legendre_coefficient_decay(fe_legendre,
+ dof_handler,
+ all_solutions,
+ all_smoothness_indicators,
+ coefficients_predicate,
+ smallest_abs_coefficient);
+ }
+
+
+
+ template <int dim, int spacedim, typename VectorType>
+ void
+ legendre_coefficient_decay(
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
+ const VectorType & solution,
+ Vector<float> & smoothness_indicators,
+ const std::function<void(std::vector<bool> &flags)> coefficients_predicate,
+ const double smallest_abs_coefficient)
+ {
+ const unsigned int max_degree =
+ dof_handler.get_fe_collection().max_degree();
+
+ // We initialize a FESeries::Legendre expansion object object which will be
+ // used to calculate the expansion coefficients. In addition to the
+ // hp::FECollection, we need to provide quadrature rules hp::QCollection for
+ // integration on the reference cell.
+ // We will need to assemble the expansion matrices for each of the finite
+ // elements we deal with, i.e. the matrices F_k,j. We have to do that for
+ // each of the finite elements in use. To that end we need a quadrature
+ // rule. As a default, we use the same quadrature formula for each finite
+ // element, namely one that is obtained by iterating a 2-point Gauss formula
+ // as many times as the maximal polynomial degree.
+ QGauss<1> base_quadrature(2);
+ QIterated<dim> quadrature(base_quadrature, max_degree);
+
+ hp::QCollection<dim> expansion_q_collection;
+ for (unsigned int i = 0; i < dof_handler.get_fe_collection().size(); ++i)
+ expansion_q_collection.push_back(quadrature);
+
+ FESeries::Legendre<dim, spacedim> legendre(max_degree + 1,
+ dof_handler.get_fe_collection(),
+ expansion_q_collection);
+
+ legendre_coefficient_decay(legendre,
+ dof_handler,
+ solution,
+ smoothness_indicators,
+ coefficients_predicate,
+ smallest_abs_coefficient);
+ }
+
+
+
+ template <int dim, int spacedim, typename VectorType>
+ void
+ fourier_coefficient_decay(
+ FESeries::Fourier<dim, spacedim> & fe_series,
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
const std::vector<const VectorType *> &all_solutions,
const std::vector<Vector<float> *> & all_smoothness_indicators,
const VectorTools::NormType regression_strategy)
{
+ using number = typename VectorType::value_type;
+ using number_coeff =
+ typename FESeries::Fourier<dim, spacedim>::CoefficientType;
+
AssertDimension(all_solutions.size(), all_smoothness_indicators.size());
for (auto &smoothness_indicator : all_smoothness_indicators)
smoothness_indicator->reinit(
dof_handler.get_triangulation().n_active_cells());
- const unsigned int dim = DoFHandlerType::dimension;
const unsigned int max_degree =
dof_handler.get_fe_collection().max_degree();
- Table<dim, typename FESeriesType::CoefficientType> expansion_coefficients;
+ Table<dim, number_coeff> expansion_coefficients;
resize(expansion_coefficients, max_degree);
- Vector<typename VectorType::value_type> local_dof_values;
+ Vector<number> local_dof_values;
std::vector<double> ln_k;
std::pair<std::vector<unsigned int>, std::vector<double>> res;
for (auto &cell : dof_handler.active_cell_iterators())
}
// Second, calculate ln(U_k).
- for (double &residual_element : res.second)
+ for (auto &residual_element : res.second)
residual_element = std::log(residual_element);
// Last, do the linear regression.
- std::pair<double, double> fit =
- FESeries::linear_regression(ln_k, res.second);
+ const auto fit = FESeries::linear_regression(ln_k, res.second);
// Compute the Sobolev index s=mu-dim/2 and store it in the vector
// of estimated values for each cell.
- template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+ template <int dim, int spacedim, typename VectorType>
void
- estimate_by_coeff_decay(FESeriesType & fe_series,
- const DoFHandlerType & dof_handler,
- const VectorType & solution,
- Vector<float> & smoothness_indicators,
- const VectorTools::NormType regression_strategy)
+ fourier_coefficient_decay(FESeries::Fourier<dim, spacedim> & fe_series,
+ const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const VectorType & solution,
+ Vector<float> & smoothness_indicators,
+ const VectorTools::NormType regression_strategy)
{
const std::vector<const VectorType *> all_solutions(1, &solution);
const std::vector<Vector<float> *> all_smoothness_indicators(
1, &smoothness_indicators);
- estimate_by_coeff_decay(fe_series,
- dof_handler,
- all_solutions,
- all_smoothness_indicators,
- regression_strategy);
+ fourier_coefficient_decay(fe_series,
+ dof_handler,
+ all_solutions,
+ all_smoothness_indicators,
+ regression_strategy);
}
- template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+ template <int dim, int spacedim, typename VectorType>
void
- estimate_by_coeff_decay(
- const DoFHandlerType & dof_handler,
+ fourier_coefficient_decay(
+ const hp::DoFHandler<dim, spacedim> & dof_handler,
const std::vector<const VectorType *> &all_solutions,
const std::vector<Vector<float> *> & all_smoothness_indicators,
const VectorTools::NormType regression_strategy)
{
- const unsigned int dim = DoFHandlerType::dimension;
const unsigned int max_degree =
dof_handler.get_fe_collection().max_degree();
// The FESeries::Fourier class' constructor first parameter $N$ defines the
// number of coefficients in 1D with the total number of coefficients being
// $N^{dim}$.
- FESeriesType fe_series(max_degree,
- dof_handler.get_fe_collection(),
- expansion_q_collection);
-
- estimate_by_coeff_decay(fe_series,
- dof_handler,
- all_solutions,
- all_smoothness_indicators,
- regression_strategy);
+ FESeries::Fourier<dim, spacedim> fe_series(max_degree,
+ dof_handler.get_fe_collection(),
+ expansion_q_collection);
+
+ fourier_coefficient_decay(fe_series,
+ dof_handler,
+ all_solutions,
+ all_smoothness_indicators,
+ regression_strategy);
}
- template <typename FESeriesType, typename DoFHandlerType, typename VectorType>
+ template <int dim, int spacedim, typename VectorType>
void
- estimate_by_coeff_decay(const DoFHandlerType & dof_handler,
- const VectorType & solution,
- Vector<float> & smoothness_indicators,
- const VectorTools::NormType regression_strategy)
+ fourier_coefficient_decay(const hp::DoFHandler<dim, spacedim> &dof_handler,
+ const VectorType & solution,
+ Vector<float> & smoothness_indicators,
+ const VectorTools::NormType regression_strategy)
{
const std::vector<const VectorType *> all_solutions(1, &solution);
const std::vector<Vector<float> *> all_smoothness_indicators(
1, &smoothness_indicators);
- estimate_by_coeff_decay<FESeriesType>(dof_handler,
- all_solutions,
- all_smoothness_indicators,
- regression_strategy);
+ fourier_coefficient_decay(dof_handler,
+ all_solutions,
+ all_smoothness_indicators,
+ regression_strategy);
}
} // namespace SmoothnessEstimator
+ template <int dim, int spacedim>
+ unsigned int
+ Legendre<dim, spacedim>::get_size_in_each_direction() const
+ {
+ return N;
+ }
+
+
+
template <int dim, int spacedim>
template <typename Number>
void
const unsigned int cell_active_fe_index,
Table<dim, CoefficientType> & legendre_coefficients)
{
+ for (unsigned int d = 0; d < dim; ++d)
+ AssertDimension(legendre_coefficients.size(d), N);
+
ensure_existence(*fe_collection,
*q_collection,
N,
//
// ---------------------------------------------------------------------
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_handler.h>
-
-#include <deal.II/hp/dof_handler.h>
-
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS;
- VEC : REAL_VECTOR_TYPES;
- DH : DOFHANDLER_TEMPLATES;
- EXP : SERIES_EXPANSION_TEMPLATES)
+ VEC : REAL_VECTOR_TYPES)
{
#if deal_II_dimension != 1 && deal_II_dimension <= deal_II_space_dimension
- template void SmoothnessEstimator::estimate_by_coeff_decay<
- EXP<deal_II_dimension, deal_II_space_dimension>,
- DH<deal_II_dimension, deal_II_space_dimension>,
- VEC>(EXP<deal_II_dimension, deal_II_space_dimension> &,
- const DH<deal_II_dimension, deal_II_space_dimension> &,
+ template void SmoothnessEstimator::legendre_coefficient_decay<
+ deal_II_dimension,
+ deal_II_space_dimension,
+ VEC>(FESeries::Legendre<deal_II_dimension, deal_II_space_dimension> &,
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+ const std::vector<const VEC *> &,
+ const std::vector<Vector<float> *> &,
+ const std::function<void(std::vector<bool> & flags)>,
+ const double);
+
+ template void SmoothnessEstimator::legendre_coefficient_decay<
+ deal_II_dimension,
+ deal_II_space_dimension,
+ VEC>(FESeries::Legendre<deal_II_dimension, deal_II_space_dimension> &,
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+ const VEC &,
+ Vector<float> &,
+ const std::function<void(std::vector<bool> & flags)>,
+ const double);
+
+ template void SmoothnessEstimator::legendre_coefficient_decay<
+ deal_II_dimension,
+ deal_II_space_dimension,
+ VEC>(const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+ const VEC &,
+ Vector<float> &,
+ const std::function<void(std::vector<bool> & flags)>,
+ const double);
+
+ template void SmoothnessEstimator::fourier_coefficient_decay<
+ deal_II_dimension,
+ deal_II_space_dimension,
+ VEC>(FESeries::Fourier<deal_II_dimension, deal_II_space_dimension> &,
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
const std::vector<const VEC *> &,
const std::vector<Vector<float> *> &,
const VectorTools::NormType);
- template void SmoothnessEstimator::estimate_by_coeff_decay<
- EXP<deal_II_dimension, deal_II_space_dimension>,
- DH<deal_II_dimension, deal_II_space_dimension>,
- VEC>(EXP<deal_II_dimension, deal_II_space_dimension> &,
- const DH<deal_II_dimension, deal_II_space_dimension> &,
+ template void SmoothnessEstimator::fourier_coefficient_decay<
+ deal_II_dimension,
+ deal_II_space_dimension,
+ VEC>(FESeries::Fourier<deal_II_dimension, deal_II_space_dimension> &,
+ const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
const VEC &,
Vector<float> &,
const VectorTools::NormType);
- template void SmoothnessEstimator::estimate_by_coeff_decay<
- EXP<deal_II_dimension, deal_II_space_dimension>,
- DH<deal_II_dimension, deal_II_space_dimension>,
- VEC>(const DH<deal_II_dimension, deal_II_space_dimension> &,
+ template void SmoothnessEstimator::fourier_coefficient_decay<
+ deal_II_dimension,
+ deal_II_space_dimension,
+ VEC>(const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
const std::vector<const VEC *> &,
const std::vector<Vector<float> *> &,
const VectorTools::NormType);
- template void SmoothnessEstimator::estimate_by_coeff_decay<
- EXP<deal_II_dimension, deal_II_space_dimension>,
- DH<deal_II_dimension, deal_II_space_dimension>,
- VEC>(const DH<deal_II_dimension, deal_II_space_dimension> &,
+ template void SmoothnessEstimator::fourier_coefficient_decay<
+ deal_II_dimension,
+ deal_II_space_dimension,
+ VEC>(const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
const VEC &,
Vector<float> &,
const VectorTools::NormType);
+
#endif
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// base header for hp-FEM test on Laplace equation.
+
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_hierarchical.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/derivative_approximation.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/smoothness_estimator.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <cmath>
+#include <fstream>
+
+#include "../tests.h"
+
+/**
+ * Basic class for Laplace problem
+ */
+template <int dim>
+class Laplace
+{
+public:
+ Laplace(const Function<dim> &force_function,
+ const Function<dim> &exact_solution,
+ const Function<dim> &boundary_conditions,
+ const unsigned int n_cycles,
+ const std::string txt_file_name);
+
+ virtual ~Laplace();
+
+ void
+ run();
+
+ hp::DoFHandler<dim> &
+ get_dof_handler();
+
+ void
+ setup_solve_estimate(Vector<float> &output_estimate);
+
+protected:
+ void
+ setup_system();
+
+ virtual void
+ setup_geometry() = 0;
+
+ void
+ assemble();
+
+ virtual void
+ solve();
+
+ /**
+ * estimate error
+ */
+ virtual void
+ estimate_error() = 0;
+
+ /**
+ * mark cells for h-refinement based on error estimation only
+ */
+ virtual void
+ mark_h_cells() = 0;
+
+ /**
+ * remove h-refinement flag from some cells and populate @p p_cells with
+ * iterators to those cells, that shall be p-refined.
+ */
+ virtual std::pair<unsigned int, unsigned int>
+ substitute_h_for_p(
+ std::vector<typename Triangulation<dim>::active_cell_iterator>
+ &p_cells) = 0;
+
+ void
+ refine_grid(const unsigned int cycle);
+
+ void
+ calculate_error();
+
+ void
+ output_results(int cycle);
+
+ void
+ print_errors();
+
+ const Function<dim> &force_function;
+ const Function<dim> &exact_solution;
+ const Function<dim> &boundary_conditions;
+
+ Triangulation<dim> triangulation;
+ hp::FECollection<dim> fe;
+ hp::DoFHandler<dim> dof_handler;
+ hp::QCollection<dim> quadrature;
+ hp::QCollection<dim> quadrature_infty;
+
+ AffineConstraints<double> constraints;
+ SparsityPattern sparsity_pattern;
+ TrilinosWrappers::SparseMatrix system_matrix;
+ TrilinosWrappers::MPI::Vector system_rhs;
+ TrilinosWrappers::MPI::Vector solution;
+ TrilinosWrappers::MPI::Vector solution_locally_relevant;
+
+ Vector<float> estimated_error_per_cell;
+ double total_error;
+
+ std::pair<unsigned int, unsigned int> hp_number;
+
+ double L2_error;
+ double H1_error;
+ double Linfty_error;
+
+ const unsigned int n_cycles;
+
+ std::string sp;
+ ConvergenceTable error_table;
+ std::string output_name;
+
+ MPI_Comm mpi_communicator;
+ const unsigned int n_mpi_processes;
+ const unsigned int this_mpi_process;
+ ConditionalOStream pcout;
+
+ IndexSet locally_owned_dofs;
+ IndexSet locally_relevant_dofs;
+};
+
+
+
+// implementatoin
+template <int dim>
+Laplace<dim>::Laplace(const Function<dim> &force_function,
+ const Function<dim> &exact_solution,
+ const Function<dim> &boundary_conditions,
+ const unsigned int n_cycles,
+ const std::string output_name)
+ : force_function(force_function)
+ , exact_solution(exact_solution)
+ , boundary_conditions(boundary_conditions)
+ , dof_handler(triangulation)
+ , n_cycles(n_cycles)
+ , sp(" ")
+ , output_name(output_name)
+ , mpi_communicator(MPI_COMM_WORLD)
+ , n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
+ , this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
+ , pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
+{
+ hp_number.first = 0;
+ hp_number.second = 0;
+
+ deallog << std::endl;
+}
+
+template <int dim>
+Laplace<dim>::~Laplace()
+{
+ dof_handler.clear();
+}
+
+template <int dim>
+hp::DoFHandler<dim> &
+Laplace<dim>::get_dof_handler()
+{
+ return dof_handler;
+}
+
+template <int dim>
+void
+Laplace<dim>::setup_solve_estimate(Vector<float> &output_estimate)
+{
+ setup_system();
+ assemble();
+ solve();
+ estimate_error();
+ output_estimate = estimated_error_per_cell;
+}
+
+template <int dim>
+void
+Laplace<dim>::setup_system()
+{
+ GridTools::partition_triangulation(n_mpi_processes, triangulation);
+
+ dof_handler.distribute_dofs(fe);
+
+ locally_owned_dofs = dof_handler.locally_owned_dofs();
+
+ locally_relevant_dofs.clear();
+ DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+ AssertThrow(locally_relevant_dofs.n_elements() == dof_handler.n_dofs(),
+ ExcInternalError());
+
+ // init vectors
+ solution_locally_relevant.reinit(locally_owned_dofs,
+ locally_relevant_dofs,
+ mpi_communicator);
+ solution.reinit(locally_owned_dofs, mpi_communicator);
+ solution = 0;
+ solution_locally_relevant = solution;
+ system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+ system_rhs = 0;
+
+ constraints.clear();
+ // constraints.reinit(locally_relevant_dofs);
+ DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+ VectorTools::interpolate_boundary_values(dof_handler,
+ 0,
+ boundary_conditions,
+ constraints);
+ if (dim == 1)
+ VectorTools::interpolate_boundary_values(dof_handler,
+ 1,
+ boundary_conditions,
+ constraints);
+ constraints.close();
+
+ TrilinosWrappers::SparsityPattern sp(locally_owned_dofs, mpi_communicator);
+ DoFTools::make_sparsity_pattern(
+ dof_handler, sp, constraints, false, this_mpi_process);
+ sp.compress();
+
+ system_matrix.reinit(sp);
+
+ estimated_error_per_cell.reinit(triangulation.n_active_cells());
+
+ // print out some info:
+ pcout << "Number of active cells: " << triangulation.n_active_cells()
+ << std::endl;
+ pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
+ << std::endl;
+}
+
+template <int dim>
+void
+Laplace<dim>::assemble()
+{
+ pcout << "Assembling...";
+
+ hp::FEValues<dim> hp_fe_values(fe,
+ quadrature,
+ update_values | update_gradients |
+ update_quadrature_points |
+ update_JxW_values);
+
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
+
+ std::vector<types::global_dof_index> local_dof_indices;
+
+ for (auto &cell : dof_handler.active_cell_iterators())
+ if (cell->subdomain_id() == this_mpi_process)
+ {
+ hp_fe_values.reinit(cell);
+ const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
+ const unsigned int & dofs_per_cell = fe_values.dofs_per_cell;
+
+ local_dof_indices.resize(dofs_per_cell);
+ cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+ cell_rhs.reinit(dofs_per_cell);
+
+ cell_matrix = 0.;
+ cell_rhs = 0.;
+
+ const unsigned int n_q_points =
+ hp_fe_values.get_present_fe_values().n_quadrature_points;
+ for (unsigned int q_index = 0; q_index < n_q_points; q_index++)
+ {
+ for (unsigned int i = 0; i < dofs_per_cell; i++)
+ {
+ for (unsigned int j = i; j < dofs_per_cell; j++)
+ {
+ cell_matrix(i, j) += (fe_values.shape_grad(i, q_index) *
+ fe_values.shape_grad(j, q_index)) *
+ fe_values.JxW(q_index);
+ }
+
+ cell_rhs(i) +=
+ force_function.value(fe_values.quadrature_point(q_index)) *
+ fe_values.shape_value(i, q_index) * fe_values.JxW(q_index);
+ }
+ }
+
+ // exploit symmetry
+ for (unsigned int i = 0; i < dofs_per_cell; i++)
+ for (unsigned int j = i; j < dofs_per_cell; j++)
+ cell_matrix(j, i) = cell_matrix(i, j);
+
+
+ cell->get_dof_indices(local_dof_indices);
+ constraints.distribute_local_to_global(
+ cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+ }
+
+ system_matrix.compress(VectorOperation::add);
+ system_rhs.compress(VectorOperation::add);
+
+ pcout << " done." << std::endl;
+}
+
+//#define DIRECT
+
+template <int dim>
+void
+Laplace<dim>::solve()
+{
+ pcout << "Solving...";
+
+ SolverControl solver_control(system_rhs.size(),
+ 1e-8 * system_rhs.l2_norm(),
+ /*log_history*/ false,
+ /*log_result*/ false);
+
+ constraints.set_zero(solution);
+ constraints.set_zero(system_rhs);
+#ifdef DIRECT
+ std::string solver_name =
+ "Amesos_Superludist"; //"Amesos_Mumps" || "Amesos_Klu"
+
+ TrilinosWrappers::SolverDirect::AdditionalData additional_data(false,
+ solver_name);
+
+ TrilinosWrappers::SolverDirect solver(solver_control, additional_data);
+
+ solver.solve(system_matrix, solution, system_rhs);
+
+ TrilinosWrappers::MPI::Vector tmp(solution);
+ const double l2 = system_matrix.residual(tmp, solution, system_rhs);
+ solver_control.check(1, l2);
+#else
+ TrilinosWrappers::SolverCG cg(solver_control);
+
+ TrilinosWrappers::PreconditionSSOR preconditioner;
+ TrilinosWrappers::PreconditionSSOR::AdditionalData data(1.2);
+ preconditioner.initialize(system_matrix, data);
+
+ cg.solve(system_matrix, solution, system_rhs, preconditioner);
+#endif
+
+ constraints.distribute(solution);
+ solution_locally_relevant = solution;
+
+ pcout << " done." << std::endl;
+}
+
+template <int dim>
+void
+Laplace<dim>::refine_grid(const unsigned int cycle)
+{
+ pcout << "Refining mesh..." << std::endl;
+
+ // 3.2. Mark cells for h-refinement
+ mark_h_cells();
+
+ // 3.3. Substitute h for p refinement
+ std::vector<typename Triangulation<dim>::active_cell_iterator> p_cells;
+ hp_number = substitute_h_for_p(p_cells);
+
+ triangulation.prepare_coarsening_and_refinement();
+
+ // 3.4. Solution Transfer
+ SolutionTransfer<dim, TrilinosWrappers::MPI::Vector, hp::DoFHandler<dim>>
+ soltrans(dof_handler);
+
+ // copy current functions
+ TrilinosWrappers::MPI::Vector solution_coarse;
+ solution_coarse.reinit(locally_owned_dofs,
+ locally_relevant_dofs,
+ mpi_communicator);
+ solution_coarse = solution;
+ soltrans.prepare_for_coarsening_and_refinement(solution_coarse);
+
+ // increase fe_index()
+ for (unsigned int i = 0; i < p_cells.size(); i++)
+ {
+ typename hp::DoFHandler<dim>::active_cell_iterator cell(
+ &triangulation, p_cells[i]->level(), p_cells[i]->index(), &dof_handler);
+
+ const unsigned int incremented_index = cell->active_fe_index() + 1;
+ Assert(incremented_index < fe.size(), ExcInternalError());
+ cell->set_active_fe_index(incremented_index);
+ }
+
+ // 3.5. Refinement
+ triangulation.execute_coarsening_and_refinement();
+
+ // FIXME: some hp strategies might need:
+ // post_execute_coarsening_and_refinement();
+
+ // 3.6. Setup
+ setup_system();
+
+ // 3.7. Solution Transfer finish
+ soltrans.interpolate(solution_coarse, solution);
+}
+
+template <int dim>
+void
+Laplace<dim>::calculate_error()
+{
+ L2_error = 0.0;
+ H1_error = 0.0;
+ Linfty_error = 0.0;
+
+ hp::FEValues<dim> hp_fe_values_linf(fe,
+ quadrature_infty,
+ update_values | update_quadrature_points);
+ hp::FEValues<dim> hp_fe_values(fe,
+ quadrature,
+ update_values | update_gradients |
+ update_quadrature_points |
+ update_JxW_values);
+
+ std::vector<double> values, exact_values;
+ std::vector<double> values_linf, exact_values_linf;
+ std::vector<Tensor<1, dim>> gradients, exact_gradients;
+
+ for (auto &cell : dof_handler.active_cell_iterators())
+ if (cell->subdomain_id() == this_mpi_process)
+ {
+ hp_fe_values.reinit(cell);
+ hp_fe_values_linf.reinit(cell);
+ const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
+ const unsigned int n_q_points = fe_values.n_quadrature_points;
+
+ const FEValues<dim> &fe_values_linf =
+ hp_fe_values_linf.get_present_fe_values();
+ const unsigned int n_q_points_linf = fe_values_linf.n_quadrature_points;
+
+ values_linf.resize(n_q_points_linf);
+ exact_values_linf.resize(n_q_points_linf);
+
+ values.resize(n_q_points);
+ exact_values.resize(n_q_points);
+ gradients.resize(n_q_points);
+ exact_gradients.resize(n_q_points);
+
+ fe_values.get_function_values(solution_locally_relevant, values);
+ fe_values_linf.get_function_values(solution_locally_relevant,
+ values_linf);
+ fe_values.get_function_gradients(solution_locally_relevant, gradients);
+
+ exact_solution.value_list(fe_values.get_quadrature_points(),
+ exact_values);
+
+ exact_solution.value_list(fe_values_linf.get_quadrature_points(),
+ exact_values_linf);
+
+ exact_solution.gradient_list(fe_values.get_quadrature_points(),
+ exact_gradients);
+
+ double cell_L2 = 0.0, cell_Linf = 0.0, cell_H1 = 0.0;
+
+ for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+ {
+ const double diff_values = exact_values[q_point] - values[q_point];
+ const Tensor<1, dim> diff_grad =
+ exact_gradients[q_point] - gradients[q_point];
+ cell_L2 += diff_values * diff_values * fe_values.JxW(q_point);
+ cell_H1 += (diff_grad * diff_grad) * fe_values.JxW(q_point);
+ }
+
+ for (unsigned int q_point = 0; q_point < n_q_points_linf; ++q_point)
+ {
+ cell_Linf = std::max(cell_Linf,
+ std::abs(exact_values_linf[q_point] -
+ values_linf[q_point]));
+ }
+
+
+ // calculate l2_norm() for cell-vectors for L2 and H1
+ // and linfty_norm() for Linf:
+ L2_error += cell_L2;
+ H1_error += cell_H1;
+ Linfty_error = std::max(Linfty_error, cell_Linf);
+ } // end of loop over cells
+
+ // finish l2_norm() / linfty_norm() calculation:
+ L2_error = sqrt(Utilities::MPI::sum(L2_error, mpi_communicator));
+ H1_error = sqrt(Utilities::MPI::sum(H1_error, mpi_communicator));
+ Linfty_error = Utilities::MPI::max(Linfty_error, mpi_communicator);
+}
+
+template <int dim>
+void
+Laplace<dim>::output_results(int cycle)
+{
+ // log:
+ error_table.add_value("cycle", cycle);
+ error_table.add_value("cells", triangulation.n_active_cells());
+ error_table.add_value("h-cells", hp_number.first);
+ error_table.add_value("p-cells", hp_number.second);
+ error_table.add_value("dofs", dof_handler.n_dofs());
+ error_table.add_value("L2", L2_error);
+ error_table.add_value("H1", H1_error);
+ error_table.add_value("Linfty", Linfty_error);
+ error_table.add_value("estimated", total_error);
+
+ if (this_mpi_process == 0)
+ deallog << cycle << sp << triangulation.n_active_cells() << sp
+ << hp_number.first << sp << hp_number.second << sp
+ << dof_handler.n_dofs() << sp << L2_error << sp << H1_error << sp
+ << Linfty_error << sp << total_error << sp << std::endl;
+}
+
+template <int dim>
+void
+Laplace<dim>::print_errors()
+{
+ error_table.set_precision("L2", 3);
+ error_table.set_precision("H1", 3);
+ error_table.set_precision("Linfty", 3);
+ error_table.set_precision("estimated", 3);
+ error_table.set_scientific("L2", true);
+ error_table.set_scientific("H1", true);
+ error_table.set_scientific("Linfty", true);
+ error_table.set_scientific("estimated", true);
+
+ pcout << std::endl << "Error analysis:" << std::endl;
+ if (this_mpi_process == 0)
+ {
+ error_table.write_text(std::cout);
+
+ const std::string fname = output_name + ".gp";
+ std::ofstream output(fname.c_str(), std::ios::out | std::ios::trunc);
+
+ // use Gnuplot datablocks:
+ output << "$data << EOD" << std::endl;
+ error_table.write_text(output);
+ output << "EOD" << std::endl << std::endl;
+
+ output
+ << "set terminal postscript eps enhanced color dashed \"Helvetica\" 22"
+ << std::endl
+ << "set style line 1 linetype 1 linecolor rgb \"#e41a1c\" linewidth 2.000 pointtype 4 pointsize 2.0"
+ << std::endl
+ << "set style line 2 linetype 1 linecolor rgb \"#377eb8\" linewidth 2.000 pointtype 6 pointsize 2.0"
+ << std::endl
+ << "set xlabel \"DoF\"" << std::endl
+ << "set ylabel \"L2+H1\"" << std::endl
+ << "set logscale xy" << std::endl
+ << "set format x \"10^{%T}\"" << std::endl
+ << "set format y \"10^{%T}\"" << std::endl
+ << "set output \'" << output_name << ".eps\'" << std::endl
+ << "plot \"$data\" using ($5):($6+$7) axis x1y1 with lp ls 1 title \"error\", \\"
+ << std::endl
+ << " \"$data\" using ($5):($9) axis x1y1 with lp ls 2 title \"{/Symbol h}_{/Symbol W}\""
+ << std::endl;
+ }
+}
+
+template <int dim>
+void
+Laplace<dim>::run()
+{
+ // 1. Define problem
+ setup_geometry();
+ setup_system();
+
+ for (unsigned int cycle = 0; cycle <= n_cycles; cycle++)
+ {
+ pcout << std::endl << "Cycle " << cycle << std::endl;
+
+ // 2. Solve Problem
+ assemble();
+ solve();
+ calculate_error();
+
+ estimate_error();
+ total_error = estimated_error_per_cell.l2_norm();
+
+ output_results(cycle);
+
+ // Do refinement (Yes/No) ?
+ if (cycle < n_cycles)
+ {
+ refine_grid(cycle);
+ }
+ }
+ print_errors();
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// test SmoothnessEstimator::legendre_coefficient_decay() on
+// problem 4 (peak) in Mitchel 2014.
+
+#include "laplace.h"
+
+
+template <int dim>
+class ForcingFunction : public Function<dim>
+{
+public:
+ ForcingFunction(const double alpha, const Point<dim> center)
+ : Function<dim>(1)
+ , alpha(alpha)
+ , center(center)
+ {}
+
+ virtual double
+ value(const Point<dim> &point, const unsigned int component = 0) const;
+
+private:
+ const double alpha;
+ const Point<dim> center;
+};
+
+template <int dim>
+double
+ForcingFunction<dim>::value(const Point<dim> &point, const unsigned int) const
+{
+ const double x = point[0];
+ const double y = point[1];
+
+ return -exp(-alpha * (point - center).norm_square()) * 4 * alpha *
+ (alpha * (point - center).norm_square() - 1);
+}
+
+template <int dim>
+class ExactSolution : public Function<dim>
+{
+public:
+ ExactSolution(const double alpha, const Point<dim> center)
+ : Function<dim>(1)
+ , alpha(alpha)
+ , center(center){};
+
+ virtual double
+ value(const Point<dim> &point, const unsigned int component = 0) const;
+
+ virtual Tensor<1, dim>
+ gradient(const Point<dim> &point, const unsigned int component = 0) const;
+
+private:
+ const double alpha;
+ const Point<dim> center;
+};
+
+template <int dim>
+double
+ExactSolution<dim>::value(const Point<dim> &point, const unsigned int) const
+{
+ return exp(-alpha * ((point - center).norm_square()));
+}
+
+template <int dim>
+Tensor<1, dim>
+ExactSolution<dim>::gradient(const Point<dim> &point, const unsigned int) const
+{
+ Tensor<1, dim> grad_u = point - center;
+ grad_u *= -2 * alpha * exp(-alpha * ((point - center).norm_square()));
+ return grad_u;
+}
+
+
+template <int dim>
+class Problem4 : public Laplace<dim>
+{
+public:
+ Problem4(const Function<dim> &force_function,
+ const Function<dim> &exact_solution,
+ const Function<dim> &boundary_conditions,
+ const unsigned int n_cycles,
+ const std::string output_name);
+
+
+private:
+ void
+ setup_geometry();
+ void
+ estimate_error();
+ void
+ mark_h_cells();
+
+ std::pair<unsigned int, unsigned int>
+ substitute_h_for_p(
+ std::vector<typename Triangulation<dim>::active_cell_iterator> &p_cells);
+
+ hp::QCollection<dim - 1> quadrature_face;
+};
+
+template <int dim>
+Problem4<dim>::Problem4(const Function<dim> &force_function,
+ const Function<dim> &exact_solution,
+ const Function<dim> &boundary_conditions,
+ const unsigned int n_cycles,
+ const std::string output_name)
+ : Laplace<dim>(force_function,
+ exact_solution,
+ boundary_conditions,
+ n_cycles,
+ output_name)
+{
+ for (unsigned int p = 1; p <= n_cycles; p++)
+ {
+ // Laplace<dim>::fe.push_back(FE_Q_Hierarchical<dim>(p));
+ Laplace<dim>::fe.push_back(FE_Q<dim>(p));
+ Laplace<dim>::quadrature.push_back(QSorted<dim>(QGauss<dim>(p + 1)));
+
+ quadrature_face.push_back(QSorted<dim - 1>(QGauss<dim - 1>(p + 1)));
+
+ const QTrapez<1> q_trapez;
+ const QIterated<dim> q_iterated(q_trapez, p + 3);
+ Laplace<dim>::quadrature_infty.push_back(QSorted<dim>(q_iterated));
+ }
+}
+
+
+
+template <int dim>
+std::pair<unsigned int, unsigned int>
+Problem4<dim>::substitute_h_for_p(
+ std::vector<typename Triangulation<dim>::active_cell_iterator> &p_cells)
+{
+ Vector<float> smoothness_indicators(
+ Laplace<dim>::triangulation.n_active_cells());
+ SmoothnessEstimator::legendre_coefficient_decay(Laplace<dim>::dof_handler,
+ Laplace<dim>::solution,
+ smoothness_indicators);
+
+ unsigned int num_p_cells = 0;
+ unsigned int num_h_cells = 0;
+ for (auto &cell : Laplace<dim>::dof_handler.active_cell_iterators())
+ if (cell->refine_flag_set())
+ {
+ typename Triangulation<dim>::active_cell_iterator tria_cell(
+ &(Laplace<dim>::triangulation), cell->level(), cell->index());
+
+ const unsigned int cur_fe_index = cell->active_fe_index();
+ const bool p_ref = smoothness_indicators(cell->index()) < exp(-1.);
+
+ if (cur_fe_index < Laplace<dim>::fe.size() - 1 && p_ref)
+ {
+ ++num_p_cells;
+ cell->clear_refine_flag();
+ p_cells.push_back(tria_cell);
+ }
+ else
+ {
+ ++num_h_cells;
+ }
+ }
+
+ return std::make_pair(num_h_cells, num_p_cells);
+}
+
+
+
+template <int dim>
+void
+Problem4<dim>::setup_geometry()
+{
+ std::vector<unsigned int> number_elements(2);
+ number_elements[0] = 16;
+ number_elements[1] = 16;
+
+ GridGenerator::subdivided_hyper_rectangle(Laplace<dim>::triangulation,
+ number_elements,
+ Point<dim>(0, 0),
+ Point<dim>(1, 1),
+ false);
+}
+
+
+
+template <int dim>
+void
+Problem4<dim>::estimate_error()
+{
+ KellyErrorEstimator<dim>::estimate(
+ Laplace<dim>::dof_handler,
+ quadrature_face,
+ std::map<types::boundary_id, const Function<dim> *>(),
+ Laplace<dim>::solution,
+ Laplace<dim>::estimated_error_per_cell);
+}
+
+template <int dim>
+void
+Problem4<dim>::mark_h_cells()
+{
+ GridRefinement::refine_and_coarsen_fixed_number(
+ Laplace<dim>::triangulation,
+ Laplace<dim>::estimated_error_per_cell,
+ 0.2,
+ 0.0);
+}
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+ const int dim = 2;
+
+ initlog();
+
+ // peak strength
+ const double alpha = 1000;
+ // peak position:
+ const double xc = 0.5;
+ const double yc = 0.5;
+
+ const Point<dim> center(xc, yc);
+
+ ForcingFunction<dim> ff(alpha, center);
+ ExactSolution<dim> ex(alpha, center);
+
+ Problem4<dim> problem(ff, ex, ex, 10, "convergence");
+ problem.run();
+}
--- /dev/null
+
+DEAL::
+DEAL::0 256 0 0 289 0.267081 1.78414 1.23386 1.29971
+DEAL::1 256 0 56 477 0.0236257 0.743531 0.249786 1.18708
+DEAL::2 268 4 48 801 0.00108376 0.293919 0.0466411 0.361176
+DEAL::3 349 22 32 1730 0.000121568 0.0610935 0.00469929 0.0464975
+DEAL::4 421 18 52 2818 3.08015e-05 0.0204013 0.00315297 0.0197872
+DEAL::5 496 11 74 4632 6.14178e-06 0.00604078 0.000714754 0.00530244
+DEAL::6 622 28 72 7303 2.03947e-06 0.00210333 0.000148863 0.00185557
+DEAL::7 712 17 108 10086 5.30920e-07 0.000744017 7.24474e-05 0.000775517
+DEAL::8 811 20 123 13127 2.97635e-07 0.000379131 3.62580e-05 0.000400767
+DEAL::9 901 15 148 17913 1.18544e-07 0.000194730 1.48736e-05 0.000216587
+DEAL::10 1030 20 161 23692 4.78884e-08 0.000106719 9.54505e-06 0.000117773
estimated_error_per_cell);
Vector<float> smoothness_indicators;
- SmoothnessEstimator::estimate_by_coeff_decay<FESeries::Fourier<dim>>(
- dof_handler, solution, smoothness_indicators);
+ SmoothnessEstimator::fourier_coefficient_decay(dof_handler,
+ solution,
+ smoothness_indicators);
// Output to VTK
if (false)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// essentially similar to fe/fe_series_05.cc but test smoothness estimation.
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_series.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/smoothness_estimator.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <gsl/gsl_sf_legendre.h>
+
+#include <iostream>
+
+#include "../tests.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]" + Utilities::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
+compare(const Table<2, double> &coeff1, const Table<2, double> &coeff2)
+{
+ double linf = 0.;
+ for (unsigned int i = 0; i < coeff1.size(0); i++)
+ for (unsigned int j = 0; j < coeff1.size(1); j++)
+ linf = std::max(linf, std::abs(coeff1(i, j) - coeff2(i, j)));
+
+ deallog << "Linf norm in exact and calculate Legendre coefficients:"
+ << std::endl
+ << linf << std::endl;
+}
+
+void
+compare(const Table<3, double> &coeff1, const Table<3, double> &coeff2)
+{
+ double linf = 0.;
+ for (unsigned int i = 0; i < coeff1.size(0); i++)
+ for (unsigned int j = 0; j < coeff1.size(1); j++)
+ for (unsigned int k = 0; k < coeff1.size(2); k++)
+ linf = std::max(linf, std::abs(coeff1(i, j, k) - coeff2(i, j, k)));
+
+ deallog << "Linf norm in exact and calculate Legendre coefficients:"
+ << std::endl
+ << linf << 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)
+{
+ // custom predicate:
+ // p-ref for linear elements and use j=1,...,pe otherwise.
+ const auto coefficients_predicate = [](std::vector<bool> &flags) -> void {
+ std::fill(flags.begin(), flags.end(), flags.size() > 2);
+ flags[0] = false;
+ };
+
+ 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);
+ }
+
+ compare(coeff_in, coeff_out);
+
+ // finally test smoothness estimator:
+ Vector<float> smoothness(1);
+ SmoothnessEstimator::legendre_coefficient_decay(
+ legendre, dof_handler, values, smoothness, coefficients_predicate);
+
+ deallog << "smoothness:" << std::endl << smoothness[0] << std::endl;
+
+ dof_handler.clear();
+}
+
+int
+main()
+{
+ std::ofstream logfile("output");
+ dealii::deallog.attach(logfile, /*do not print job id*/ false);
+ dealii::deallog.depth_console(0);
+
+ // for linear elements we expect p-refinement by convention
+ {
+ 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);
+ deallog << "expected smoothness:" << std::endl << 0. << std::endl;
+ }
+
+ // for quadratic we can already assign exponential decay: a_i = C exp ( -k
+ // i) set one with different k's
+ {
+ const double k1 = 1.;
+ const double k2 = 2.;
+
+ 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++)
+ coeff_in(i, 0) = exp(-k1 * i);
+
+ for (unsigned int i = 0; i < coeff_1d; i++)
+ coeff_in(0, i) = exp(-k2 * i);
+
+ // make sure predicate skips 0-th:
+ coeff_in(0, 0) = 12345;
+
+ LegendreFunction<dim> function(coeff_in);
+ test(function, p);
+
+ deallog << "expected smoothness:" << std::endl
+ << exp(-std::min(k1, k2)) << std::endl;
+ }
+
+ // linear elements in 3D (expect zero output)
+ {
+ 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);
+ deallog << "expected smoothness:" << std::endl << 0. << std::endl;
+ }
+
+ // cubic in 3D
+ {
+ const double k1 = 2.;
+ const double k2 = 3.;
+ const double k3 = 4.;
+ const unsigned int dim = 3;
+ const unsigned int coeff_1d = 4;
+ const unsigned int p = 3;
+ Table<dim, double> coeff_in(coeff_1d, coeff_1d, coeff_1d);
+ for (unsigned int i = 0; i < coeff_1d; i++)
+ coeff_in(i, 0, 0) = exp(-k1 * i);
+
+ for (unsigned int j = 0; j < coeff_1d; j++)
+ coeff_in(0, j, 0) = exp(-k2 * j);
+
+ for (unsigned int k = 0; k < coeff_1d; k++)
+ coeff_in(0, 0, k) = exp(-k3 * k);
+
+ // make sure predicate skips 0-th:
+ coeff_in(0, 0, 0) = 12345;
+
+ LegendreFunction<dim> function(coeff_in);
+ test(function, p);
+
+ deallog << "expected smoothness:" << std::endl
+ << exp(-std::min(k1, std::min(k2, k3))) << std::endl;
+ }
+
+
+ // 4-th order in 3D but with some coefficients being zero
+ {
+ const double k1 = 2.;
+ const double k2 = k1 + 1.;
+ const unsigned int dim = 3;
+ const unsigned int coeff_1d = 5;
+ const unsigned int p = 4;
+ Table<dim, double> coeff_in(coeff_1d, coeff_1d, coeff_1d);
+ // all non-zero:
+ for (unsigned int i = 0; i < coeff_1d; i++)
+ coeff_in(i, 0, 0) = exp(-k2 * i);
+
+ // some non-zero (2nd and 4th), the slowest decay will be from this
+ // direction
+ for (unsigned int j = 2; j < coeff_1d; j = j + 2)
+ coeff_in(0, j, 0) = exp(-k1 * j);
+
+ // all but one zero:
+ for (unsigned int k = 3; k < coeff_1d; k = k + 10)
+ coeff_in(0, 0, k) = exp(-k2 * k);
+
+ // make sure predicate skips 0-th:
+ coeff_in(0, 0, 0) = 12345;
+
+ LegendreFunction<dim> function(coeff_in);
+ test(function, p);
+
+ deallog << "expected smoothness:" << std::endl << exp(-k1) << std::endl;
+ }
+
+ // cubic in 3D (zero)
+ {
+ const unsigned int dim = 3;
+ const unsigned int coeff_1d = 4;
+ const unsigned int p = 3;
+ Table<dim, double> coeff_in(coeff_1d, coeff_1d, coeff_1d);
+
+ LegendreFunction<dim> function(coeff_in);
+ test(function, p);
+
+ deallog << "expected smoothness:" << std::endl << 0. << std::endl;
+ }
+
+ dealii::deallog << "Ok" << std::endl;
+}
--- /dev/null
+DEAL::-----------------------------------
+DEAL::2d, p=1, max_p=4
+DEAL::-----------------------------------
+DEAL::Linf norm in exact and calculate Legendre coefficients:
+DEAL::1.77636e-15
+DEAL::smoothness:
+DEAL::0.00000
+DEAL::expected smoothness:
+DEAL::0.00000
+DEAL::-----------------------------------
+DEAL::2d, p=2, max_p=5
+DEAL::-----------------------------------
+DEAL::Linf norm in exact and calculate Legendre coefficients:
+DEAL::7.21645e-16
+DEAL::smoothness:
+DEAL::0.367879
+DEAL::expected smoothness:
+DEAL::0.367879
+DEAL::-----------------------------------
+DEAL::3d, p=1, max_p=4
+DEAL::-----------------------------------
+DEAL::Linf norm in exact and calculate Legendre coefficients:
+DEAL::3.55271e-15
+DEAL::smoothness:
+DEAL::0.00000
+DEAL::expected smoothness:
+DEAL::0.00000
+DEAL::-----------------------------------
+DEAL::3d, p=3, max_p=6
+DEAL::-----------------------------------
+DEAL::Linf norm in exact and calculate Legendre coefficients:
+DEAL::3.19189e-15
+DEAL::smoothness:
+DEAL::0.135335
+DEAL::expected smoothness:
+DEAL::0.135335
+DEAL::-----------------------------------
+DEAL::3d, p=4, max_p=7
+DEAL::-----------------------------------
+DEAL::Linf norm in exact and calculate Legendre coefficients:
+DEAL::4.31599e-15
+DEAL::smoothness:
+DEAL::0.135335
+DEAL::expected smoothness:
+DEAL::0.135335
+DEAL::-----------------------------------
+DEAL::3d, p=3, max_p=6
+DEAL::-----------------------------------
+DEAL::Linf norm in exact and calculate Legendre coefficients:
+DEAL::0.00000
+DEAL::smoothness:
+DEAL::0.00000
+DEAL::expected smoothness:
+DEAL::0.00000
+DEAL::Ok
projects / dealii.git / commitdiff