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