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