54 template <
int dim,
typename CoefficientType>
59 for (
unsigned int d = 0;
d < dim; ++
d)
86 std::pair<bool, unsigned int>
90 for (
unsigned int i = 0; i < dim; ++i)
93 return std::make_pair((v < N), v);
99 template <
int dim,
int spacedim,
typename VectorType>
103 const VectorType &solution,
106 const double smallest_abs_coefficient,
107 const bool only_flagged_cells)
109 using number =
typename VectorType::value_type;
113 smoothness_indicators.
reinit(
116 unsigned int n_modes;
120 std::vector<double> converted_indices;
121 std::pair<std::vector<unsigned int>, std::vector<double>> res;
125 if (!only_flagged_cells || cell->refine_flag_set() ||
126 cell->coarsen_flag_set())
129 cell->active_fe_index());
130 resize(expansion_coefficients, n_modes);
132 local_dof_values.
reinit(cell->get_fe().n_dofs_per_cell());
133 cell->get_dof_values(solution, local_dof_values);
136 cell->active_fe_index(),
137 expansion_coefficients);
143 res = FESeries::process_coefficients<dim>(
144 expansion_coefficients,
146 return index_sum_less_than_N(indices, n_modes);
149 smallest_abs_coefficient);
154 float regularity = std::numeric_limits<float>::infinity();
155 if (res.first.size() > 1)
159 converted_indices.assign(res.first.begin(), res.first.end());
161 for (
auto &residual_element : res.second)
162 residual_element =
std::log(residual_element);
164 const std::pair<double, double> fit =
166 regularity =
static_cast<float>(-fit.first);
169 smoothness_indicators(cell->active_cell_index()) = regularity;
172 smoothness_indicators(cell->active_cell_index()) =
173 numbers::signaling_nan<float>();
179 template <
int dim,
int spacedim,
typename VectorType>
184 const VectorType &solution,
187 const double smallest_abs_coefficient,
188 const bool only_flagged_cells)
190 Assert(smallest_abs_coefficient >= 0.,
191 ExcMessage(
"smallest_abs_coefficient should be non-negative."));
193 using number =
typename VectorType::value_type;
197 smoothness_indicators.
reinit(
200 unsigned int n_modes;
205 const unsigned int max_degree =
208 std::vector<double> x, y;
209 x.reserve(max_degree);
210 y.reserve(max_degree);
215 if (!only_flagged_cells || cell->refine_flag_set() ||
216 cell->coarsen_flag_set())
219 cell->active_fe_index());
220 resize(expansion_coefficients, n_modes);
222 const unsigned int pe = cell->get_fe().degree;
230 local_dof_values.
reinit(cell->get_fe().n_dofs_per_cell());
231 cell->get_dof_values(solution, local_dof_values);
234 cell->active_fe_index(),
235 expansion_coefficients);
239 double k_v = std::numeric_limits<double>::max();
240 for (
unsigned int d = 0; d < dim; ++d)
247 for (
unsigned int i = 0; i <= pe; ++i)
248 if (coefficients_predicate[i])
252 const double coeff_abs =
253 std::abs(expansion_coefficients(ind));
255 if (coeff_abs > smallest_abs_coefficient)
267 const std::pair<double, double> fit =
275 smoothness_indicators(cell->active_cell_index()) =
276 static_cast<float>(k_v);
279 smoothness_indicators(cell->active_cell_index()) =
280 numbers::signaling_nan<float>();
286 template <
int dim,
int spacedim>
289 const unsigned int component)
296 std::vector<unsigned int> n_coefficients_per_direction;
297 n_coefficients_per_direction.reserve(fe_collection.
size());
298 for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
299 n_coefficients_per_direction.push_back(fe_collection[i].degree + 2);
317 for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
319 const QGauss<dim> quadrature(n_coefficients_per_direction[i]);
321 q_collection.
push_back(quadrature_sorted);
352 std::pair<bool, unsigned int>
353 index_norm_greater_than_zero_and_less_than_N_squared(
355 const unsigned int N)
358 for (
unsigned int i = 0; i < dim; ++i)
359 v += ind[i] * ind[i];
361 return std::make_pair((v > 0 && v < N * N), v);
367 template <
int dim,
int spacedim,
typename VectorType>
371 const VectorType &solution,
374 const double smallest_abs_coefficient,
375 const bool only_flagged_cells)
377 using number =
typename VectorType::value_type;
381 smoothness_indicators.
reinit(
384 unsigned int n_modes;
388 std::vector<double> ln_k;
389 std::pair<std::vector<unsigned int>, std::vector<double>> res;
393 if (!only_flagged_cells || cell->refine_flag_set() ||
394 cell->coarsen_flag_set())
397 cell->active_fe_index());
398 resize(expansion_coefficients, n_modes);
404 local_dof_values.
reinit(cell->get_fe().n_dofs_per_cell());
405 cell->get_dof_values(solution, local_dof_values);
408 cell->active_fe_index(),
409 expansion_coefficients);
415 res = FESeries::process_coefficients<dim>(
416 expansion_coefficients,
418 return index_norm_greater_than_zero_and_less_than_N_squared(
422 smallest_abs_coefficient);
427 float regularity = std::numeric_limits<float>::infinity();
428 if (res.first.size() > 1)
439 ln_k.resize(res.first.size());
440 for (
unsigned int f = 0; f < res.first.size(); ++f)
441 ln_k[f] = 0.5 *
std::log(
static_cast<double>(res.first[f]));
444 for (
auto &residual_element : res.second)
445 residual_element =
std::log(residual_element);
447 const std::pair<double, double> fit =
450 regularity =
static_cast<float>(-fit.first) -
451 ((dim > 1) ? (.5 * dim) : 0);
455 smoothness_indicators(cell->active_cell_index()) = regularity;
458 smoothness_indicators(cell->active_cell_index()) =
459 numbers::signaling_nan<float>();
465 template <
int dim,
int spacedim,
typename VectorType>
470 const VectorType &solution,
473 const double smallest_abs_coefficient,
474 const bool only_flagged_cells)
476 Assert(smallest_abs_coefficient >= 0.,
477 ExcMessage(
"smallest_abs_coefficient should be non-negative."));
479 using number =
typename VectorType::value_type;
483 smoothness_indicators.
reinit(
486 unsigned int n_modes;
491 const unsigned int max_degree =
494 std::vector<double> x, y;
495 x.reserve(max_degree);
496 y.reserve(max_degree);
501 if (!only_flagged_cells || cell->refine_flag_set() ||
502 cell->coarsen_flag_set())
505 cell->active_fe_index());
506 resize(expansion_coefficients, n_modes);
508 const unsigned int pe = cell->get_fe().degree;
516 local_dof_values.
reinit(cell->get_fe().n_dofs_per_cell());
517 cell->get_dof_values(solution, local_dof_values);
520 cell->active_fe_index(),
521 expansion_coefficients);
525 double k_v = std::numeric_limits<double>::max();
526 for (
unsigned int d = 0; d < dim; ++d)
535 for (
unsigned int i = 1; i <= pe; ++i)
536 if (coefficients_predicate[i])
540 const double coeff_abs =
541 std::abs(expansion_coefficients(ind));
543 if (coeff_abs > smallest_abs_coefficient)
555 const std::pair<double, double> fit =
563 smoothness_indicators(cell->active_cell_index()) =
564 static_cast<float>(k_v);
567 smoothness_indicators(cell->active_cell_index()) =
568 numbers::signaling_nan<float>();
574 template <
int dim,
int spacedim>
577 const unsigned int component)
586 std::vector<unsigned int> n_coefficients_per_direction;
587 n_coefficients_per_direction.reserve(fe_collection.
size());
588 for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
589 n_coefficients_per_direction.push_back(fe_collection[i].degree + 2);
607 for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
610 n_coefficients_per_direction[i] - 1);
612 q_collection.
push_back(quadrature_sorted);
625#include "numerics/smoothness_estimator.inst"
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &fourier_coefficients)
typename std::complex< double > CoefficientType
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &legendre_coefficients)
unsigned int n_active_cells() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
unsigned int size() const
unsigned int max_degree() const
void push_back(const Quadrature< dim_in > &new_quadrature)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::pair< double, double > linear_regression(const std::vector< double > &x, const std::vector< double > &y)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
FESeries::Fourier< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)
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={}, const double smallest_abs_coefficient=1e-10, const bool only_flagged_cells=false)
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)
FESeries::Legendre< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)
void coefficient_decay(FESeries::Legendre< dim, spacedim > &fe_legendre, 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)
void coefficient_decay_per_direction(FESeries::Legendre< dim, spacedim > &fe_legendre, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &solution, Vector< float > &smoothness_indicators, const ComponentMask &coefficients_predicate={}, const double smallest_abs_coefficient=1e-10, const bool only_flagged_cells=false)
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)