Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
smoothness_estimator.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2018 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
18
20
22
25
35#include <deal.II/lac/vector.h>
36
38
39#include <algorithm>
40#include <cmath>
41#include <limits>
42#include <utility>
43
44
46
47
48namespace SmoothnessEstimator
49{
50 namespace
51 {
55 template <int dim, typename CoefficientType>
56 void
57 resize(Table<dim, CoefficientType> &coeff, const unsigned int N)
58 {
60 for (unsigned int d = 0; d < dim; ++d)
61 size[d] = N;
62 coeff.reinit(size);
63 }
64 } // namespace
65
66
67
68 namespace Legendre
69 {
70 namespace
71 {
86 template <int dim>
87 std::pair<bool, unsigned int>
88 index_sum_less_than_N(const TableIndices<dim> &ind, const unsigned int N)
89 {
90 unsigned int v = 0;
91 for (unsigned int i = 0; i < dim; ++i)
92 v += ind[i];
93
94 return std::make_pair((v < N), v);
95 }
96 } // namespace
97
98
99
100 template <int dim, int spacedim, typename VectorType>
101 void
103 const DoFHandler<dim, spacedim> & dof_handler,
104 const VectorType & solution,
105 Vector<float> & smoothness_indicators,
106 const VectorTools::NormType regression_strategy,
107 const double smallest_abs_coefficient,
108 const bool only_flagged_cells)
109 {
110 using number = typename VectorType::value_type;
111 using number_coeff =
113
114 smoothness_indicators.reinit(
115 dof_handler.get_triangulation().n_active_cells());
116
117 unsigned int n_modes;
118 Table<dim, number_coeff> expansion_coefficients;
119
120 Vector<number> local_dof_values;
121 std::vector<double> converted_indices;
122 std::pair<std::vector<unsigned int>, std::vector<double>> res;
123 for (const auto &cell : dof_handler.active_cell_iterators() |
125 {
126 if (!only_flagged_cells || cell->refine_flag_set() ||
127 cell->coarsen_flag_set())
128 {
129 n_modes = fe_legendre.get_n_coefficients_per_direction(
130 cell->active_fe_index());
131 resize(expansion_coefficients, n_modes);
132
133 local_dof_values.reinit(cell->get_fe().n_dofs_per_cell());
134 cell->get_dof_values(solution, local_dof_values);
135
136 fe_legendre.calculate(local_dof_values,
137 cell->active_fe_index(),
138 expansion_coefficients);
139
140 // We fit our exponential decay of expansion coefficients to the
141 // provided regression_strategy on each possible value of |k|.
142 // To this end, we use FESeries::process_coefficients() to
143 // rework coefficients into the desired format.
144 res = FESeries::process_coefficients<dim>(
145 expansion_coefficients,
146 [n_modes](const TableIndices<dim> &indices) {
147 return index_sum_less_than_N(indices, n_modes);
148 },
149 regression_strategy,
150 smallest_abs_coefficient);
151
152 Assert(res.first.size() == res.second.size(), ExcInternalError());
153
154 // Last, do the linear regression.
155 float regularity = std::numeric_limits<float>::infinity();
156 if (res.first.size() > 1)
157 {
158 // Prepare linear equation for the logarithmic least squares
159 // fit.
160 converted_indices.assign(res.first.begin(), res.first.end());
161
162 for (auto &residual_element : res.second)
163 residual_element = std::log(residual_element);
164
165 const std::pair<double, double> fit =
166 FESeries::linear_regression(converted_indices, res.second);
167 regularity = static_cast<float>(-fit.first);
168 }
169
170 smoothness_indicators(cell->active_cell_index()) = regularity;
171 }
172 else
173 smoothness_indicators(cell->active_cell_index()) =
174 numbers::signaling_nan<float>();
175 }
176 }
177
178
179
180 template <int dim, int spacedim, typename VectorType>
181 void
184 const DoFHandler<dim, spacedim> & dof_handler,
185 const VectorType & solution,
186 Vector<float> & smoothness_indicators,
187 const ComponentMask & coefficients_predicate,
188 const double smallest_abs_coefficient,
189 const bool only_flagged_cells)
190 {
191 Assert(smallest_abs_coefficient >= 0.,
192 ExcMessage("smallest_abs_coefficient should be non-negative."));
193
194 using number = typename VectorType::value_type;
195 using number_coeff =
197
198 smoothness_indicators.reinit(
199 dof_handler.get_triangulation().n_active_cells());
200
201 unsigned int n_modes;
202 Table<dim, number_coeff> expansion_coefficients;
203 Vector<number> local_dof_values;
204
205 // auxiliary vector to do linear regression
206 const unsigned int max_degree =
207 dof_handler.get_fe_collection().max_degree();
208
209 std::vector<double> x, y;
210 x.reserve(max_degree);
211 y.reserve(max_degree);
212
213 for (const auto &cell : dof_handler.active_cell_iterators() |
215 {
216 if (!only_flagged_cells || cell->refine_flag_set() ||
217 cell->coarsen_flag_set())
218 {
219 n_modes = fe_legendre.get_n_coefficients_per_direction(
220 cell->active_fe_index());
221 resize(expansion_coefficients, n_modes);
222
223 const unsigned int pe = cell->get_fe().degree;
224 Assert(pe > 0, ExcInternalError());
225
226 // since we use coefficients with indices [1,pe] in each
227 // direction, the number of coefficients we need to calculate is
228 // at least N=pe+1
229 AssertIndexRange(pe, n_modes);
230
231 local_dof_values.reinit(cell->get_fe().n_dofs_per_cell());
232 cell->get_dof_values(solution, local_dof_values);
233
234 fe_legendre.calculate(local_dof_values,
235 cell->active_fe_index(),
236 expansion_coefficients);
237
238 // choose the smallest decay of coefficients in each direction,
239 // i.e. the maximum decay slope k_v as in exp(-k_v)
240 double k_v = std::numeric_limits<double>::infinity();
241 for (unsigned int d = 0; d < dim; ++d)
242 {
243 x.resize(0);
244 y.resize(0);
245
246 // will use all non-zero coefficients allowed by the
247 // predicate function
248 for (unsigned int i = 0; i <= pe; ++i)
249 if (coefficients_predicate[i])
250 {
252 ind[d] = i;
253 const double coeff_abs =
254 std::abs(expansion_coefficients(ind));
255
256 if (coeff_abs > smallest_abs_coefficient)
257 {
258 x.push_back(i);
259 y.push_back(std::log(coeff_abs));
260 }
261 }
262
263 // in case we don't have enough non-zero coefficient to fit,
264 // skip this direction
265 if (x.size() < 2)
266 continue;
267
268 const std::pair<double, double> fit =
270
271 // decay corresponds to negative slope
272 // take the lesser negative slope along each direction
273 k_v = std::min(k_v, -fit.first);
274 }
275
276 smoothness_indicators(cell->active_cell_index()) =
277 static_cast<float>(k_v);
278 }
279 else
280 smoothness_indicators(cell->active_cell_index()) =
281 numbers::signaling_nan<float>();
282 }
283 }
284
285
286
287 template <int dim, int spacedim>
290 const unsigned int component)
291 {
292 // Default number of coefficients per direction.
293 //
294 // With a number of modes equal to the polynomial degree plus two for each
295 // finite element, the smoothness estimation algorithm tends to produce
296 // stable results.
297 std::vector<unsigned int> n_coefficients_per_direction;
298 for (unsigned int i = 0; i < fe_collection.size(); ++i)
299 n_coefficients_per_direction.push_back(fe_collection[i].degree + 2);
300
301 // Default quadrature collection.
302 //
303 // We initialize a FESeries::Legendre expansion object object which will
304 // be used to calculate the expansion coefficients. In addition to the
305 // hp::FECollection, we need to provide quadrature rules hp::QCollection
306 // for integration on the reference cell.
307 // We will need to assemble the expansion matrices for each of the finite
308 // elements we deal with, i.e. the matrices F_k,j. We have to do that for
309 // each of the finite elements in use. To that end we need a quadrature
310 // rule. As a default, we use the same quadrature formula for each finite
311 // element, namely a Gauss formula that yields exact results for the
312 // highest order Legendre polynomial used.
313 //
314 // We start with the zeroth Legendre polynomial which is just a constant,
315 // so the highest Legendre polynomial will be of order (n_modes - 1).
316 hp::QCollection<dim> q_collection;
317 for (unsigned int i = 0; i < fe_collection.size(); ++i)
318 {
319 const QGauss<dim> quadrature(n_coefficients_per_direction[i]);
320 const QSorted<dim> quadrature_sorted(quadrature);
321 q_collection.push_back(quadrature_sorted);
322 }
323
324 return FESeries::Legendre<dim, spacedim>(n_coefficients_per_direction,
325 fe_collection,
326 q_collection,
327 component);
328 }
329 } // namespace Legendre
330
331
332
333 namespace Fourier
334 {
335 namespace
336 {
351 template <int dim>
352 std::pair<bool, unsigned int>
353 index_norm_greater_than_zero_and_less_than_N_squared(
354 const TableIndices<dim> &ind,
355 const unsigned int N)
356 {
357 unsigned int v = 0;
358 for (unsigned int i = 0; i < dim; ++i)
359 v += ind[i] * ind[i];
360
361 return std::make_pair((v > 0 && v < N * N), v);
362 }
363 } // namespace
364
365
366
367 template <int dim, int spacedim, typename VectorType>
368 void
370 const DoFHandler<dim, spacedim> & dof_handler,
371 const VectorType & solution,
372 Vector<float> & smoothness_indicators,
373 const VectorTools::NormType regression_strategy,
374 const double smallest_abs_coefficient,
375 const bool only_flagged_cells)
376 {
377 using number = typename VectorType::value_type;
378 using number_coeff =
380
381 smoothness_indicators.reinit(
382 dof_handler.get_triangulation().n_active_cells());
383
384 unsigned int n_modes;
385 Table<dim, number_coeff> expansion_coefficients;
386
387 Vector<number> local_dof_values;
388 std::vector<double> ln_k;
389 std::pair<std::vector<unsigned int>, std::vector<double>> res;
390 for (const auto &cell : dof_handler.active_cell_iterators() |
392 {
393 if (!only_flagged_cells || cell->refine_flag_set() ||
394 cell->coarsen_flag_set())
395 {
396 n_modes = fe_fourier.get_n_coefficients_per_direction(
397 cell->active_fe_index());
398 resize(expansion_coefficients, n_modes);
399
400 // Inside the loop, we first need to get the values of the local
401 // degrees of freedom and then need to compute the series
402 // expansion by multiplying this vector with the matrix @f${\cal
403 // F}@f$ corresponding to this finite element.
404 local_dof_values.reinit(cell->get_fe().n_dofs_per_cell());
405 cell->get_dof_values(solution, local_dof_values);
406
407 fe_fourier.calculate(local_dof_values,
408 cell->active_fe_index(),
409 expansion_coefficients);
410
411 // We fit our exponential decay of expansion coefficients to the
412 // provided regression_strategy on each possible value of |k|.
413 // To this end, we use FESeries::process_coefficients() to
414 // rework coefficients into the desired format.
415 res = FESeries::process_coefficients<dim>(
416 expansion_coefficients,
417 [n_modes](const TableIndices<dim> &indices) {
418 return index_norm_greater_than_zero_and_less_than_N_squared(
419 indices, n_modes);
420 },
421 regression_strategy,
422 smallest_abs_coefficient);
423
424 Assert(res.first.size() == res.second.size(), ExcInternalError());
425
426 // Last, do the linear regression.
427 float regularity = std::numeric_limits<float>::infinity();
428 if (res.first.size() > 1)
429 {
430 // Prepare linear equation for the logarithmic least squares
431 // fit.
432 //
433 // First, calculate ln(|k|).
434 //
435 // For Fourier expansion, this translates to
436 // ln(2*pi*sqrt(predicate)) = ln(2*pi) + 0.5*ln(predicate).
437 // Since we are just interested in the slope of a linear
438 // regression later, we omit the ln(2*pi) factor.
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]));
442
443 // Second, calculate ln(U_k).
444 for (auto &residual_element : res.second)
445 residual_element = std::log(residual_element);
446
447 const std::pair<double, double> fit =
448 FESeries::linear_regression(ln_k, res.second);
449 // Compute regularity s = mu - dim/2
450 regularity = static_cast<float>(-fit.first) -
451 ((dim > 1) ? (.5 * dim) : 0);
452 }
453
454 // Store result in the vector of estimated values for each cell.
455 smoothness_indicators(cell->active_cell_index()) = regularity;
456 }
457 else
458 smoothness_indicators(cell->active_cell_index()) =
459 numbers::signaling_nan<float>();
460 }
461 }
462
463
464
465 template <int dim, int spacedim, typename VectorType>
466 void
469 const DoFHandler<dim, spacedim> & dof_handler,
470 const VectorType & solution,
471 Vector<float> & smoothness_indicators,
472 const ComponentMask & coefficients_predicate,
473 const double smallest_abs_coefficient,
474 const bool only_flagged_cells)
475 {
476 Assert(smallest_abs_coefficient >= 0.,
477 ExcMessage("smallest_abs_coefficient should be non-negative."));
478
479 using number = typename VectorType::value_type;
480 using number_coeff =
482
483 smoothness_indicators.reinit(
484 dof_handler.get_triangulation().n_active_cells());
485
486 unsigned int n_modes;
487 Table<dim, number_coeff> expansion_coefficients;
488 Vector<number> local_dof_values;
489
490 // auxiliary vector to do linear regression
491 const unsigned int max_degree =
492 dof_handler.get_fe_collection().max_degree();
493
494 std::vector<double> x, y;
495 x.reserve(max_degree);
496 y.reserve(max_degree);
497
498 for (const auto &cell : dof_handler.active_cell_iterators() |
500 {
501 if (!only_flagged_cells || cell->refine_flag_set() ||
502 cell->coarsen_flag_set())
503 {
504 n_modes = fe_fourier.get_n_coefficients_per_direction(
505 cell->active_fe_index());
506 resize(expansion_coefficients, n_modes);
507
508 const unsigned int pe = cell->get_fe().degree;
509 Assert(pe > 0, ExcInternalError());
510
511 // since we use coefficients with indices [1,pe] in each
512 // direction, the number of coefficients we need to calculate is
513 // at least N=pe+1
514 AssertIndexRange(pe, n_modes);
515
516 local_dof_values.reinit(cell->get_fe().n_dofs_per_cell());
517 cell->get_dof_values(solution, local_dof_values);
518
519 fe_fourier.calculate(local_dof_values,
520 cell->active_fe_index(),
521 expansion_coefficients);
522
523 // choose the smallest decay of coefficients in each direction,
524 // i.e. the maximum decay slope k_v as in exp(-k_v)
525 double k_v = std::numeric_limits<double>::infinity();
526 for (unsigned int d = 0; d < dim; ++d)
527 {
528 x.resize(0);
529 y.resize(0);
530
531 // will use all non-zero coefficients allowed by the
532 // predicate function
533 //
534 // skip i=0 because of logarithm
535 for (unsigned int i = 1; i <= pe; ++i)
536 if (coefficients_predicate[i])
537 {
539 ind[d] = i;
540 const double coeff_abs =
541 std::abs(expansion_coefficients(ind));
542
543 if (coeff_abs > smallest_abs_coefficient)
544 {
545 x.push_back(std::log(i));
546 y.push_back(std::log(coeff_abs));
547 }
548 }
549
550 // in case we don't have enough non-zero coefficient to fit,
551 // skip this direction
552 if (x.size() < 2)
553 continue;
554
555 const std::pair<double, double> fit =
557
558 // decay corresponds to negative slope
559 // take the lesser negative slope along each direction
560 k_v = std::min(k_v, -fit.first);
561 }
562
563 smoothness_indicators(cell->active_cell_index()) =
564 static_cast<float>(k_v);
565 }
566 else
567 smoothness_indicators(cell->active_cell_index()) =
568 numbers::signaling_nan<float>();
569 }
570 }
571
572
573
574 template <int dim, int spacedim>
577 const unsigned int component)
578 {
579 // Default number of coefficients per direction.
580 //
581 // Since we omit the zero-th mode in the Fourier decay strategy, make sure
582 // that we have at least two modes to work with per finite element. With a
583 // number of modes equal to the polynomial degree plus two for each finite
584 // element, the smoothness estimation algorithm tends to produce stable
585 // results.
586 std::vector<unsigned int> n_coefficients_per_direction;
587 for (unsigned int i = 0; i < fe_collection.size(); ++i)
588 n_coefficients_per_direction.push_back(fe_collection[i].degree + 2);
589
590 // Default quadrature collection.
591 //
592 // We initialize a series expansion object object which will be used to
593 // calculate the expansion coefficients. In addition to the
594 // hp::FECollection, we need to provide quadrature rules hp::QCollection
595 // for integration on the reference cell.
596 // We will need to assemble the expansion matrices for each of the finite
597 // elements we deal with, i.e. the matrices F_k,j. We have to do that for
598 // each of the finite elements in use. To that end we need a quadrature
599 // rule. As a default, we use the same quadrature formula for each finite
600 // element, namely one that is obtained by iterating a 5-point Gauss
601 // formula as many times as the maximal exponent we use for the term
602 // exp(ikx). Since the first mode corresponds to k = 0, the maximal wave
603 // number is k = n_modes - 1.
604 const QGauss<1> base_quadrature(5);
605 hp::QCollection<dim> q_collection;
606 for (unsigned int i = 0; i < fe_collection.size(); ++i)
607 {
608 const QIterated<dim> quadrature(base_quadrature,
609 n_coefficients_per_direction[i] - 1);
610 const QSorted<dim> quadrature_sorted(quadrature);
611 q_collection.push_back(quadrature_sorted);
612 }
613
614 return FESeries::Fourier<dim, spacedim>(n_coefficients_per_direction,
615 fe_collection,
616 q_collection,
617 component);
618 }
619 } // namespace Fourier
620} // namespace SmoothnessEstimator
621
622
623// explicit instantiations
624#include "smoothness_estimator.inst"
625
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
Definition: fe_series.h:91
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
double CoefficientType
Definition: fe_series.h:261
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
Definition: vector.h:109
unsigned int size() const
Definition: collection.h:264
unsigned int max_degree() const
void push_back(const Quadrature< dim_in > &new_quadrature)
Definition: q_collection.h:216
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
std::pair< double, double > linear_regression(const std::vector< double > &x, const std::vector< double > &y)
Definition: fe_series.cc:30
static const char N
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
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(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)
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=ComponentMask(), 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)
::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 > &)