Reference documentation for deal.II version GIT 6bdf9a4b6c 2022-08-12 19:20:02+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\}}\)
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 
21 #include <deal.II/fe/fe_series.h>
22 
24 
26 
30 #include <deal.II/lac/la_vector.h>
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 
49 namespace SmoothnessEstimator
50 {
51  namespace
52  {
56  template <int dim, typename CoefficientType>
57  void
58  resize(Table<dim, CoefficientType> &coeff, const unsigned int N)
59  {
60  TableIndices<dim> size;
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  {
252  TableIndices<dim> ind;
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  {
539  TableIndices<dim> ind;
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 Triangulation< dim, spacedim > & get_triangulation() const
const hp::FECollection< dim, spacedim > & get_fe_collection() 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
Definition: vector.h:110
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
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
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
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
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)