Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2972-g3b1da50de6 2025-03-29 20:10:00+00:00
\(\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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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.
143 res = FESeries::process_coefficients<dim>(
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()) =
173 numbers::signaling_nan<float>();
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>::max();
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()) =
280 numbers::signaling_nan<float>();
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 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);
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>::max();
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 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);
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 "numerics/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: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:316
unsigned int max_degree() const
void push_back(const Quadrature< dim_in > &new_quadrature)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
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::size_t size
Definition mpi.cc:745
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)
::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 > &)