Reference documentation for deal.II version 9.3.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 
19 #include <deal.II/fe/fe_series.h>
20 
21 #include <deal.II/hp/dof_handler.h>
23 
27 #include <deal.II/lac/la_vector.h>
33 #include <deal.II/lac/vector.h>
34 
36 
37 #include <algorithm>
38 #include <cmath>
39 #include <limits>
40 #include <utility>
41 
42 
44 
45 
46 namespace SmoothnessEstimator
47 {
48  namespace
49  {
53  template <int dim, typename CoefficientType>
54  void
55  resize(Table<dim, CoefficientType> &coeff, const unsigned int N)
56  {
57  TableIndices<dim> size;
58  for (unsigned int d = 0; d < dim; d++)
59  size[d] = N;
60  coeff.reinit(size);
61  }
62  } // namespace
63 
64 
65 
66  namespace Legendre
67  {
68  namespace
69  {
84  template <int dim>
85  std::pair<bool, unsigned int>
86  index_sum_less_than_N(const TableIndices<dim> &ind, const unsigned int N)
87  {
88  unsigned int v = 0;
89  for (unsigned int i = 0; i < dim; ++i)
90  v += ind[i];
91 
92  return std::make_pair((v < N), v);
93  }
94  } // namespace
95 
96 
97 
98  template <int dim, int spacedim, typename VectorType>
99  void
101  const DoFHandler<dim, spacedim> & dof_handler,
102  const VectorType & solution,
103  Vector<float> & smoothness_indicators,
104  const VectorTools::NormType regression_strategy,
105  const double smallest_abs_coefficient,
106  const bool only_flagged_cells)
107  {
108  using number = typename VectorType::value_type;
109  using number_coeff =
111 
112  smoothness_indicators.reinit(
113  dof_handler.get_triangulation().n_active_cells());
114 
115  unsigned int n_modes;
116  Table<dim, number_coeff> expansion_coefficients;
117 
118  Vector<number> local_dof_values;
119  std::vector<double> converted_indices;
120  std::pair<std::vector<unsigned int>, std::vector<double>> res;
121  for (const auto &cell : dof_handler.active_cell_iterators())
122  if (cell->is_locally_owned())
123  {
124  if (!only_flagged_cells || cell->refine_flag_set() ||
125  cell->coarsen_flag_set())
126  {
127  n_modes = fe_legendre.get_n_coefficients_per_direction(
128  cell->active_fe_index());
129  resize(expansion_coefficients, n_modes);
130 
131  local_dof_values.reinit(cell->get_fe().n_dofs_per_cell());
132  cell->get_dof_values(solution, local_dof_values);
133 
134  fe_legendre.calculate(local_dof_values,
135  cell->active_fe_index(),
136  expansion_coefficients);
137 
138  // We fit our exponential decay of expansion coefficients to the
139  // provided regression_strategy on each possible value of |k|.
140  // To this end, we use FESeries::process_coefficients() to
141  // rework coefficients into the desired format.
142  res = FESeries::process_coefficients<dim>(
143  expansion_coefficients,
144  [n_modes](const TableIndices<dim> &indices) {
145  return index_sum_less_than_N(indices, n_modes);
146  },
147  regression_strategy,
148  smallest_abs_coefficient);
149 
150  Assert(res.first.size() == res.second.size(),
151  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(),
160  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,
167  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())
215  if (cell->is_locally_owned())
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())
392  if (cell->is_locally_owned())
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(),
426  ExcInternalError());
427 
428  // Last, do the linear regression.
429  float regularity = std::numeric_limits<float>::infinity();
430  if (res.first.size() > 1)
431  {
432  // Prepare linear equation for the logarithmic least squares
433  // fit.
434  //
435  // First, calculate ln(|k|).
436  //
437  // For Fourier expansion, this translates to
438  // ln(2*pi*sqrt(predicate)) = ln(2*pi) + 0.5*ln(predicate).
439  // Since we are just interested in the slope of a linear
440  // regression later, we omit the ln(2*pi) factor.
441  ln_k.resize(res.first.size());
442  for (unsigned int f = 0; f < res.first.size(); ++f)
443  ln_k[f] =
444  0.5 * std::log(static_cast<double>(res.first[f]));
445 
446  // Second, calculate ln(U_k).
447  for (auto &residual_element : res.second)
448  residual_element = std::log(residual_element);
449 
450  const std::pair<double, double> fit =
451  FESeries::linear_regression(ln_k, res.second);
452  // Compute regularity s = mu - dim/2
453  regularity = static_cast<float>(-fit.first) -
454  ((dim > 1) ? (.5 * dim) : 0);
455  }
456 
457  // Store result in the vector of estimated values for each cell.
458  smoothness_indicators(cell->active_cell_index()) = regularity;
459  }
460  else
461  smoothness_indicators(cell->active_cell_index()) =
462  numbers::signaling_nan<float>();
463  }
464  }
465 
466 
467 
468  template <int dim, int spacedim, typename VectorType>
469  void
472  const DoFHandler<dim, spacedim> & dof_handler,
473  const VectorType & solution,
474  Vector<float> & smoothness_indicators,
475  const ComponentMask & coefficients_predicate,
476  const double smallest_abs_coefficient,
477  const bool only_flagged_cells)
478  {
479  Assert(smallest_abs_coefficient >= 0.,
480  ExcMessage("smallest_abs_coefficient should be non-negative."));
481 
482  using number = typename VectorType::value_type;
483  using number_coeff =
485 
486  smoothness_indicators.reinit(
487  dof_handler.get_triangulation().n_active_cells());
488 
489  unsigned int n_modes;
490  Table<dim, number_coeff> expansion_coefficients;
491  Vector<number> local_dof_values;
492 
493  // auxiliary vector to do linear regression
494  const unsigned int max_degree =
495  dof_handler.get_fe_collection().max_degree();
496 
497  std::vector<double> x, y;
498  x.reserve(max_degree);
499  y.reserve(max_degree);
500 
501  for (const auto &cell : dof_handler.active_cell_iterators())
502  if (cell->is_locally_owned())
503  {
504  if (!only_flagged_cells || cell->refine_flag_set() ||
505  cell->coarsen_flag_set())
506  {
507  n_modes = fe_fourier.get_n_coefficients_per_direction(
508  cell->active_fe_index());
509  resize(expansion_coefficients, n_modes);
510 
511  const unsigned int pe = cell->get_fe().degree;
512  Assert(pe > 0, ExcInternalError());
513 
514  // since we use coefficients with indices [1,pe] in each
515  // direction, the number of coefficients we need to calculate is
516  // at least N=pe+1
517  AssertIndexRange(pe, n_modes);
518 
519  local_dof_values.reinit(cell->get_fe().n_dofs_per_cell());
520  cell->get_dof_values(solution, local_dof_values);
521 
522  fe_fourier.calculate(local_dof_values,
523  cell->active_fe_index(),
524  expansion_coefficients);
525 
526  // choose the smallest decay of coefficients in each direction,
527  // i.e. the maximum decay slope k_v as in exp(-k_v)
528  double k_v = std::numeric_limits<double>::infinity();
529  for (unsigned int d = 0; d < dim; ++d)
530  {
531  x.resize(0);
532  y.resize(0);
533 
534  // will use all non-zero coefficients allowed by the
535  // predicate function
536  //
537  // skip i=0 because of logarithm
538  for (unsigned int i = 1; i <= pe; ++i)
539  if (coefficients_predicate[i])
540  {
541  TableIndices<dim> ind;
542  ind[d] = i;
543  const double coeff_abs =
544  std::abs(expansion_coefficients(ind));
545 
546  if (coeff_abs > smallest_abs_coefficient)
547  {
548  x.push_back(std::log(i));
549  y.push_back(std::log(coeff_abs));
550  }
551  }
552 
553  // in case we don't have enough non-zero coefficient to fit,
554  // skip this direction
555  if (x.size() < 2)
556  continue;
557 
558  const std::pair<double, double> fit =
560 
561  // decay corresponds to negative slope
562  // take the lesser negative slope along each direction
563  k_v = std::min(k_v, -fit.first);
564  }
565 
566  smoothness_indicators(cell->active_cell_index()) =
567  static_cast<float>(k_v);
568  }
569  else
570  smoothness_indicators(cell->active_cell_index()) =
571  numbers::signaling_nan<float>();
572  }
573  }
574 
575 
576 
577  template <int dim, int spacedim>
580  const unsigned int component)
581  {
582  // Default number of coefficients per direction.
583  //
584  // Since we omit the zero-th mode in the Fourier decay strategy, make sure
585  // that we have at least two modes to work with per finite element. With a
586  // number of modes equal to the polynomial degree plus two for each finite
587  // element, the smoothness estimation algorithm tends to produce stable
588  // results.
589  std::vector<unsigned int> n_coefficients_per_direction;
590  for (unsigned int i = 0; i < fe_collection.size(); ++i)
591  n_coefficients_per_direction.push_back(fe_collection[i].degree + 2);
592 
593  // Default quadrature collection.
594  //
595  // We initialize a series expansion object object which will be used to
596  // calculate the expansion coefficients. In addition to the
597  // hp::FECollection, we need to provide quadrature rules hp::QCollection
598  // for integration on the reference cell.
599  // We will need to assemble the expansion matrices for each of the finite
600  // elements we deal with, i.e. the matrices F_k,j. We have to do that for
601  // each of the finite elements in use. To that end we need a quadrature
602  // rule. As a default, we use the same quadrature formula for each finite
603  // element, namely one that is obtained by iterating a 5-point Gauss
604  // formula as many times as the maximal exponent we use for the term
605  // exp(ikx). Since the first mode corresponds to k = 0, the maximal wave
606  // number is k = n_modes - 1.
607  const QGauss<1> base_quadrature(5);
608  hp::QCollection<dim> q_collection;
609  for (unsigned int i = 0; i < fe_collection.size(); ++i)
610  {
611  const QIterated<dim> quadrature(base_quadrature,
612  n_coefficients_per_direction[i] - 1);
613  const QSorted<dim> quadrature_sorted(quadrature);
614  q_collection.push_back(quadrature_sorted);
615  }
616 
617  return FESeries::Fourier<dim, spacedim>(n_coefficients_per_direction,
618  fe_collection,
619  q_collection,
620  component);
621  }
622  } // namespace Fourier
623 } // namespace SmoothnessEstimator
624 
625 
626 // explicit instantiations
627 #include "smoothness_estimator.inst"
628 
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 calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &fourier_coefficients)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
std::pair< double, double > linear_regression(const std::vector< double > &x, const std::vector< double > &y)
Definition: fe_series.cc:30
const hp::FECollection< dim, spacedim > & get_fe_collection() const
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1465
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
typename std::complex< double > CoefficientType
Definition: fe_series.h:91
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &legendre_coefficients)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void push_back(const Quadrature< dim_in > &new_quadrature)
Definition: q_collection.h:223
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
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)
const Triangulation< dim, spacedim > & get_triangulation() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
static const char N
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
FESeries::Legendre< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)
static ::ExceptionBase & ExcInternalError()