Reference documentation for deal.II version GIT b6bf1e606d 2022-08-11 15:25: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\}}\)
fe_series.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 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 
16 #ifndef dealii_fe_series_h
17 #define dealii_fe_series_h
18 
19 
20 
21 #include <deal.II/base/config.h>
22 
25 #include <deal.II/base/table.h>
27 #include <deal.II/base/tensor.h>
28 
31 
33 #include <deal.II/lac/vector.h>
34 
36 
37 #include <memory>
38 #include <string>
39 #include <vector>
40 
41 
43 
44 
57 namespace FESeries
58 {
89  template <int dim, int spacedim = dim>
90  class Fourier : public Subscriptor
91  {
92  public:
93  using CoefficientType = typename std::complex<double>;
94 
113  Fourier(const std::vector<unsigned int> & n_coefficients_per_direction,
116  const unsigned int component = numbers::invalid_unsigned_int);
117 
123  template <typename Number>
124  void
125  calculate(const ::Vector<Number> &local_dof_values,
126  const unsigned int cell_active_fe_index,
127  Table<dim, CoefficientType> & fourier_coefficients);
128 
133  unsigned int
134  get_n_coefficients_per_direction(const unsigned int index) const;
135 
146  void
148 
158  template <class Archive>
159  void
160  save_transformation_matrices(Archive &ar, const unsigned int version);
161 
166  template <class Archive>
167  void
168  load_transformation_matrices(Archive &ar, const unsigned int version);
169 
173  bool
174  operator==(const Fourier<dim, spacedim> &fourier) const;
175 
176  private:
181  const std::vector<unsigned int> n_coefficients_per_direction;
182 
187 
192 
197 
201  std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;
202 
206  std::vector<CoefficientType> unrolled_coefficients;
207 
212  const unsigned int component;
213  };
214 
215 
216 
259  template <int dim, int spacedim = dim>
260  class Legendre : public Subscriptor
261  {
262  public:
264 
283  Legendre(const std::vector<unsigned int> &n_coefficients_per_direction,
286  const unsigned int component = numbers::invalid_unsigned_int);
287 
293  template <typename Number>
294  void
295  calculate(const ::Vector<Number> &local_dof_values,
296  const unsigned int cell_active_fe_index,
297  Table<dim, CoefficientType> & legendre_coefficients);
298 
303  unsigned int
304  get_n_coefficients_per_direction(const unsigned int index) const;
305 
316  void
318 
328  template <class Archive>
329  void
330  save_transformation_matrices(Archive &ar, const unsigned int version);
331 
336  template <class Archive>
337  void
338  load_transformation_matrices(Archive &ar, const unsigned int version);
339 
343  bool
345 
346  private:
351  const std::vector<unsigned int> n_coefficients_per_direction;
352 
357 
362 
366  std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
367 
371  std::vector<CoefficientType> unrolled_coefficients;
372 
377  const unsigned int component;
378  };
379 
380 
381 
399  template <int dim, typename CoefficientType>
400  std::pair<std::vector<unsigned int>, std::vector<double>>
402  const std::function<std::pair<bool, unsigned int>(
403  const TableIndices<dim> &)> & predicate,
404  const VectorTools::NormType norm_type,
405  const double smallest_abs_coefficient = 1e-10);
406 
412  std::pair<double, double>
413  linear_regression(const std::vector<double> &x, const std::vector<double> &y);
414 
415 } // namespace FESeries
416 
421 #ifndef DOXYGEN
422 
423 // ------------------- inline and template functions ----------------
424 
425 namespace internal
426 {
427  namespace FESeriesImplementation
428  {
429  template <int dim, typename CoefficientType>
430  void
431  fill_map_index(
432  const Table<dim, CoefficientType> &coefficients,
433  const TableIndices<dim> & ind,
434  const std::function<
435  std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
436  std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
437  {
438  const std::pair<bool, unsigned int> pred_pair = predicate(ind);
439  // don't add a value if predicate is false
440  if (pred_pair.first == false)
441  return;
442 
443  const unsigned int pred_value = pred_pair.second;
444  const CoefficientType &coeff_value = coefficients(ind);
445  // If pred_value is not in the pred_to_values map, the element will be
446  // created. Otherwise a reference to the existing element is returned.
447  pred_to_values[pred_value].push_back(coeff_value);
448  }
449 
450 
451 
452  template <typename CoefficientType>
453  void
454  fill_map(
455  const Table<1, CoefficientType> &coefficients,
456  const std::function<
457  std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
458  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
459  {
460  for (unsigned int i = 0; i < coefficients.size(0); ++i)
461  {
462  const TableIndices<1> ind(i);
463  fill_map_index(coefficients, ind, predicate, pred_to_values);
464  }
465  }
466 
467 
468 
469  template <typename CoefficientType>
470  void
471  fill_map(
472  const Table<2, CoefficientType> &coefficients,
473  const std::function<
474  std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
475  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
476  {
477  for (unsigned int i = 0; i < coefficients.size(0); ++i)
478  for (unsigned int j = 0; j < coefficients.size(1); ++j)
479  {
480  const TableIndices<2> ind(i, j);
481  fill_map_index(coefficients, ind, predicate, pred_to_values);
482  }
483  }
484 
485 
486 
487  template <typename CoefficientType>
488  void
489  fill_map(
490  const Table<3, CoefficientType> &coefficients,
491  const std::function<
492  std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
493  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
494  {
495  for (unsigned int i = 0; i < coefficients.size(0); ++i)
496  for (unsigned int j = 0; j < coefficients.size(1); ++j)
497  for (unsigned int k = 0; k < coefficients.size(2); ++k)
498  {
499  const TableIndices<3> ind(i, j, k);
500  fill_map_index(coefficients, ind, predicate, pred_to_values);
501  }
502  }
503 
504 
505 
506  template <typename Number>
507  double
508  complex_mean_value(const Number &value)
509  {
510  return value;
511  }
512 
513 
514 
515  template <typename Number>
516  double
517  complex_mean_value(const std::complex<Number> &value)
518  {
519  AssertThrow(false,
520  ExcMessage(
521  "FESeries::process_coefficients() can not be used with "
522  "complex-valued coefficients and VectorTools::mean norm."));
523  return std::abs(value);
524  }
525  } // namespace FESeriesImplementation
526 } // namespace internal
527 
528 
529 
530 template <int dim, typename CoefficientType>
531 std::pair<std::vector<unsigned int>, std::vector<double>>
533  const Table<dim, CoefficientType> &coefficients,
534  const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
535  & predicate,
536  const VectorTools::NormType norm_type,
537  const double smallest_abs_coefficient)
538 {
539  Assert(smallest_abs_coefficient >= 0.,
540  ExcMessage("smallest_abs_coefficient should be non-negative."));
541 
542  std::vector<unsigned int> predicate_values;
543  std::vector<double> norm_values;
544 
545  // first, parse all table elements into a map of predicate values and
546  // coefficients. We could have stored (predicate values ->TableIndicies) map,
547  // but its processing would have been much harder later on.
548  std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
549  internal::FESeriesImplementation::fill_map(coefficients,
550  predicate,
551  pred_to_values);
552 
553  // now go through the map and populate the @p norm_values based on @p norm:
554  for (const auto &pred_to_value : pred_to_values)
555  {
556  Vector<CoefficientType> values(pred_to_value.second.cbegin(),
557  pred_to_value.second.cend());
558 
559  double norm_value = 0;
560  switch (norm_type)
561  {
563  {
564  norm_value = values.l2_norm();
565  break;
566  }
568  {
569  norm_value = values.l1_norm();
570  break;
571  }
573  {
574  norm_value = values.linfty_norm();
575  break;
576  }
577  case VectorTools::mean:
578  {
579  norm_value = internal::FESeriesImplementation::complex_mean_value(
580  values.mean_value());
581  break;
582  }
583  default:
584  AssertThrow(false, ExcNotImplemented());
585  break;
586  }
587 
588  // will use all non-zero coefficients
589  if (std::abs(norm_value) > smallest_abs_coefficient)
590  {
591  predicate_values.push_back(pred_to_value.first);
592  norm_values.push_back(norm_value);
593  }
594  }
595 
596  return std::make_pair(predicate_values, norm_values);
597 }
598 
599 
600 
601 template <int dim, int spacedim>
602 template <class Archive>
603 inline void
605  Archive &ar,
606  const unsigned int /*version*/)
607 {
608  // Store information about those resources which have been used to generate
609  // the transformation matrices.
610  // mode vector
611  ar &n_coefficients_per_direction;
612 
613  // finite element collection
614  unsigned int size = fe_collection->size();
615  ar & size;
616  for (unsigned int i = 0; i < size; ++i)
617  ar &(*fe_collection)[i].get_name();
618 
619  // quadrature collection
620  size = q_collection.size();
621  ar &size;
622  for (unsigned int i = 0; i < size; ++i)
623  ar &q_collection[i];
624 
625  // Store the actual transform matrices.
626  ar &fourier_transform_matrices;
627 }
628 
629 
630 
631 template <int dim, int spacedim>
632 template <class Archive>
633 inline void
635  Archive &ar,
636  const unsigned int /*version*/)
637 {
638  // Check whether the currently registered resources are compatible with
639  // the transformation matrices to load.
640  // mode vector
641  std::vector<unsigned int> compare_coefficients;
642  ar & compare_coefficients;
643  Assert(compare_coefficients == n_coefficients_per_direction,
644  ExcMessage("A different number of coefficients vector has been used "
645  "to generate the transformation matrices you are about "
646  "to load!"));
647 
648  // finite element collection
649  unsigned int size;
650  ar & size;
651  AssertDimension(size, fe_collection->size());
652  std::string name;
653  for (unsigned int i = 0; i < size; ++i)
654  {
655  ar &name;
656  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
657  ExcMessage("A different FECollection has been used to generate "
658  "the transformation matrices you are about to load!"));
659  }
660 
661  // quadrature collection
662  ar &size;
663  AssertDimension(size, q_collection.size());
664  Quadrature<dim> quadrature;
665  for (unsigned int i = 0; i < size; ++i)
666  {
667  ar &quadrature;
668  Assert(quadrature == q_collection[i],
669  ExcMessage("A different QCollection has been used to generate "
670  "the transformation matrices you are about to load!"));
671  }
672 
673  // Restore the transform matrices since all prerequisites are fulfilled.
674  ar &fourier_transform_matrices;
675 }
676 
677 
678 
679 template <int dim, int spacedim>
680 template <class Archive>
681 inline void
683  Archive &ar,
684  const unsigned int /*version*/)
685 {
686  // Store information about those resources which have been used to generate
687  // the transformation matrices.
688  // mode vector
689  ar &n_coefficients_per_direction;
690 
691  // finite element collection
692  unsigned int size = fe_collection->size();
693  ar & size;
694  for (unsigned int i = 0; i < size; ++i)
695  ar &(*fe_collection)[i].get_name();
696 
697  // quadrature collection
698  size = q_collection.size();
699  ar &size;
700  for (unsigned int i = 0; i < size; ++i)
701  ar &q_collection[i];
702 
703  // Store the actual transform matrices.
704  ar &legendre_transform_matrices;
705 }
706 
707 
708 
709 template <int dim, int spacedim>
710 template <class Archive>
711 inline void
713  Archive &ar,
714  const unsigned int /*version*/)
715 {
716  // Check whether the currently registered resources are compatible with
717  // the transformation matrices to load.
718  // mode vector
719  std::vector<unsigned int> compare_coefficients;
720  ar & compare_coefficients;
721  Assert(compare_coefficients == n_coefficients_per_direction,
722  ExcMessage("A different number of coefficients vector has been used "
723  "to generate the transformation matrices you are about "
724  "to load!"));
725 
726  // finite element collection
727  unsigned int size;
728  ar & size;
729  AssertDimension(size, fe_collection->size());
730  std::string name;
731  for (unsigned int i = 0; i < size; ++i)
732  {
733  ar &name;
734  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
735  ExcMessage("A different FECollection has been used to generate "
736  "the transformation matrices you are about to load!"));
737  }
738 
739  // quadrature collection
740  ar &size;
741  AssertDimension(size, q_collection.size());
742  Quadrature<dim> quadrature;
743  for (unsigned int i = 0; i < size; ++i)
744  {
745  ar &quadrature;
746  Assert(quadrature == q_collection[i],
747  ExcMessage("A different QCollection has been used to generate "
748  "the transformation matrices you are about to load!"));
749  }
750 
751  // Restore the transform matrices since all prerequisites are fulfilled.
752  ar &legendre_transform_matrices;
753 }
754 
755 
756 #endif // DOXYGEN
757 
759 
760 #endif // dealii_fe_series_h
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:186
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &fourier_coefficients)
Table< dim, Tensor< 1, dim > > k_vectors
Definition: fe_series.h:196
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:206
const unsigned int component
Definition: fe_series.h:212
typename std::complex< double > CoefficientType
Definition: fe_series.h:93
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:181
bool operator==(const Fourier< dim, spacedim > &fourier) const
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition: fe_series.h:201
const hp::QCollection< dim > q_collection
Definition: fe_series.h:191
Fourier(const std::vector< unsigned int > &n_coefficients_per_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection, const unsigned int component=numbers::invalid_unsigned_int)
void precalculate_all_transformation_matrices()
void load_transformation_matrices(Archive &ar, const unsigned int version)
void save_transformation_matrices(Archive &ar, const unsigned int version)
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:351
void save_transformation_matrices(Archive &ar, const unsigned int version)
bool operator==(const Legendre< dim, spacedim > &legendre) const
const unsigned int component
Definition: fe_series.h:377
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition: fe_series.h:366
const hp::QCollection< dim > q_collection
Definition: fe_series.h:361
Legendre(const std::vector< unsigned int > &n_coefficients_per_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection, const unsigned int component=numbers::invalid_unsigned_int)
void precalculate_all_transformation_matrices()
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:371
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:356
void load_transformation_matrices(Archive &ar, const unsigned int version)
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &legendre_coefficients)
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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:30
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
static const unsigned int invalid_unsigned_int
Definition: types.h:201
double legendre(unsigned int l, double x)
Definition: cmath.h:75
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)