53 template <
int dim,
typename CoefficientType>
58 for (
unsigned int d = 0;
d < dim;
d++)
85 std::pair<bool, unsigned int>
89 for (
unsigned int i = 0; i < dim; ++i)
92 return std::make_pair((v <
N), v);
98 template <
int dim,
int spacedim,
typename VectorType>
103 Vector<float> & smoothness_indicators,
105 const double smallest_abs_coefficient,
106 const bool only_flagged_cells)
108 using number =
typename VectorType::value_type;
112 smoothness_indicators.reinit(
115 unsigned int n_modes;
119 std::vector<double> converted_indices;
120 std::pair<std::vector<unsigned int>, std::vector<double>> res;
122 if (cell->is_locally_owned())
124 if (!only_flagged_cells || cell->refine_flag_set() ||
125 cell->coarsen_flag_set())
128 cell->active_fe_index());
129 resize(expansion_coefficients, n_modes);
131 local_dof_values.
reinit(cell->get_fe().dofs_per_cell);
132 cell->get_dof_values(solution, local_dof_values);
135 cell->active_fe_index(),
136 expansion_coefficients);
142 res = FESeries::process_coefficients<dim>(
143 expansion_coefficients,
145 return index_sum_less_than_N(indices, n_modes);
148 smallest_abs_coefficient);
150 Assert(res.first.size() == res.second.size(),
154 float regularity = std::numeric_limits<float>::infinity();
155 if (res.first.size() > 1)
159 converted_indices.assign(res.first.begin(),
162 for (
auto &residual_element : res.second)
163 residual_element = std::log(residual_element);
165 const std::pair<double, double> fit =
168 regularity =
static_cast<float>(-fit.first);
171 smoothness_indicators(cell->active_cell_index()) = regularity;
174 smoothness_indicators(cell->active_cell_index()) =
175 numbers::signaling_nan<float>();
181 template <
int dim,
int spacedim,
typename VectorType>
187 Vector<float> & smoothness_indicators,
189 const double smallest_abs_coefficient,
190 const bool only_flagged_cells)
192 Assert(smallest_abs_coefficient >= 0.,
193 ExcMessage(
"smallest_abs_coefficient should be non-negative."));
195 using number =
typename VectorType::value_type;
199 smoothness_indicators.reinit(
202 unsigned int n_modes;
207 const unsigned int max_degree =
210 std::vector<double> x, y;
211 x.reserve(max_degree);
212 y.reserve(max_degree);
215 if (cell->is_locally_owned())
217 if (!only_flagged_cells || cell->refine_flag_set() ||
218 cell->coarsen_flag_set())
221 cell->active_fe_index());
222 resize(expansion_coefficients, n_modes);
224 const unsigned int pe = cell->get_fe().degree;
232 local_dof_values.
reinit(cell->get_fe().dofs_per_cell);
233 cell->get_dof_values(solution, local_dof_values);
236 cell->active_fe_index(),
237 expansion_coefficients);
241 double k_v = std::numeric_limits<double>::infinity();
242 for (
unsigned int d = 0;
d < dim; ++
d)
249 for (
unsigned int i = 0; i <= pe; ++i)
250 if (coefficients_predicate[i])
254 const double coeff_abs =
255 std::abs(expansion_coefficients(ind));
257 if (coeff_abs > smallest_abs_coefficient)
260 y.push_back(std::log(coeff_abs));
269 const std::pair<double, double> fit =
277 smoothness_indicators(cell->active_cell_index()) =
278 static_cast<float>(k_v);
281 smoothness_indicators(cell->active_cell_index()) =
282 numbers::signaling_nan<float>();
288 template <
int dim,
int spacedim>
293 std::vector<unsigned int> n_coefficients_per_direction;
294 for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
295 n_coefficients_per_direction.push_back(fe_collection[i].degree + 1);
313 for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
315 const QGauss<dim> quadrature(n_coefficients_per_direction[i]);
317 q_collection.
push_back(quadrature_sorted);
347 std::pair<bool, unsigned int>
348 index_norm_greater_than_zero_and_less_than_N_squared(
350 const unsigned int N)
353 for (
unsigned int i = 0; i < dim; ++i)
354 v += ind[i] * ind[i];
356 return std::make_pair((v > 0 && v <
N *
N), v);
362 template <
int dim,
int spacedim,
typename VectorType>
367 Vector<float> & smoothness_indicators,
369 const double smallest_abs_coefficient,
370 const bool only_flagged_cells)
372 using number =
typename VectorType::value_type;
376 smoothness_indicators.reinit(
379 unsigned int n_modes;
383 std::vector<double> ln_k;
384 std::pair<std::vector<unsigned int>, std::vector<double>> res;
386 if (cell->is_locally_owned())
388 if (!only_flagged_cells || cell->refine_flag_set() ||
389 cell->coarsen_flag_set())
392 cell->active_fe_index());
393 resize(expansion_coefficients, n_modes);
399 local_dof_values.
reinit(cell->get_fe().dofs_per_cell);
400 cell->get_dof_values(solution, local_dof_values);
403 cell->active_fe_index(),
404 expansion_coefficients);
410 res = FESeries::process_coefficients<dim>(
411 expansion_coefficients,
413 return index_norm_greater_than_zero_and_less_than_N_squared(
417 smallest_abs_coefficient);
419 Assert(res.first.size() == res.second.size(),
423 float regularity = std::numeric_limits<float>::infinity();
424 if (res.first.size() > 1)
435 ln_k.resize(res.first.size());
436 for (
unsigned int f = 0; f < res.first.size(); ++f)
438 0.5 * std::log(
static_cast<double>(res.first[f]));
441 for (
auto &residual_element : res.second)
442 residual_element = std::log(residual_element);
444 const std::pair<double, double> fit =
447 regularity =
static_cast<float>(-fit.first) -
448 ((dim > 1) ? (.5 * dim) : 0);
452 smoothness_indicators(cell->active_cell_index()) = regularity;
455 smoothness_indicators(cell->active_cell_index()) =
456 numbers::signaling_nan<float>();
462 template <
int dim,
int spacedim,
typename VectorType>
468 Vector<float> & smoothness_indicators,
470 const double smallest_abs_coefficient,
471 const bool only_flagged_cells)
473 Assert(smallest_abs_coefficient >= 0.,
474 ExcMessage(
"smallest_abs_coefficient should be non-negative."));
476 using number =
typename VectorType::value_type;
480 smoothness_indicators.reinit(
483 unsigned int n_modes;
488 const unsigned int max_degree =
491 std::vector<double> x, y;
492 x.reserve(max_degree);
493 y.reserve(max_degree);
496 if (cell->is_locally_owned())
498 if (!only_flagged_cells || cell->refine_flag_set() ||
499 cell->coarsen_flag_set())
502 cell->active_fe_index());
503 resize(expansion_coefficients, n_modes);
505 const unsigned int pe = cell->get_fe().degree;
513 local_dof_values.
reinit(cell->get_fe().dofs_per_cell);
514 cell->get_dof_values(solution, local_dof_values);
517 cell->active_fe_index(),
518 expansion_coefficients);
522 double k_v = std::numeric_limits<double>::infinity();
523 for (
unsigned int d = 0;
d < dim; ++
d)
532 for (
unsigned int i = 1; i <= pe; ++i)
533 if (coefficients_predicate[i])
537 const double coeff_abs =
538 std::abs(expansion_coefficients(ind));
540 if (coeff_abs > smallest_abs_coefficient)
542 x.push_back(std::log(i));
543 y.push_back(std::log(coeff_abs));
552 const std::pair<double, double> fit =
560 smoothness_indicators(cell->active_cell_index()) =
561 static_cast<float>(k_v);
564 smoothness_indicators(cell->active_cell_index()) =
565 numbers::signaling_nan<float>();
571 template <
int dim,
int spacedim>
576 std::vector<unsigned int> n_coefficients_per_direction;
577 for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
578 n_coefficients_per_direction.push_back(
579 std::max<unsigned int>(3, fe_collection[i].degree + 1));
597 for (
unsigned int i = 0; i < fe_collection.
size(); ++i)
600 n_coefficients_per_direction[i] - 1);
602 q_collection.
push_back(quadrature_sorted);
614 #include "smoothness_estimator.inst"