Reference documentation for deal.II version Git 8f0f58a98a 2021-10-21 17:29:01 -0600
\(\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 
47 
48 
55 namespace FESeries
56 {
87  template <int dim, int spacedim = dim>
88  class Fourier : public Subscriptor
89  {
90  public:
91  using CoefficientType = typename std::complex<double>;
92 
111  Fourier(const std::vector<unsigned int> & n_coefficients_per_direction,
114  const unsigned int component = numbers::invalid_unsigned_int);
115 
121  template <typename Number>
122  void
123  calculate(const ::Vector<Number> &local_dof_values,
124  const unsigned int cell_active_fe_index,
125  Table<dim, CoefficientType> & fourier_coefficients);
126 
131  unsigned int
132  get_n_coefficients_per_direction(const unsigned int index) const;
133 
144  void
146 
156  template <class Archive>
157  void
158  save_transformation_matrices(Archive &ar, const unsigned int version);
159 
164  template <class Archive>
165  void
166  load_transformation_matrices(Archive &ar, const unsigned int version);
167 
171  bool
172  operator==(const Fourier<dim, spacedim> &fourier) const;
173 
174  private:
179  const std::vector<unsigned int> n_coefficients_per_direction;
180 
185 
190 
195 
199  std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;
200 
204  std::vector<CoefficientType> unrolled_coefficients;
205 
210  const unsigned int component;
211  };
212 
213 
214 
257  template <int dim, int spacedim = dim>
258  class Legendre : public Subscriptor
259  {
260  public:
262 
281  Legendre(const std::vector<unsigned int> &n_coefficients_per_direction,
284  const unsigned int component = numbers::invalid_unsigned_int);
285 
291  template <typename Number>
292  void
293  calculate(const ::Vector<Number> &local_dof_values,
294  const unsigned int cell_active_fe_index,
295  Table<dim, CoefficientType> & legendre_coefficients);
296 
301  unsigned int
302  get_n_coefficients_per_direction(const unsigned int index) const;
303 
314  void
316 
326  template <class Archive>
327  void
328  save_transformation_matrices(Archive &ar, const unsigned int version);
329 
334  template <class Archive>
335  void
336  load_transformation_matrices(Archive &ar, const unsigned int version);
337 
341  bool
343 
344  private:
349  const std::vector<unsigned int> n_coefficients_per_direction;
350 
355 
360 
364  std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
365 
369  std::vector<CoefficientType> unrolled_coefficients;
370 
375  const unsigned int component;
376  };
377 
378 
379 
397  template <int dim, typename CoefficientType>
398  std::pair<std::vector<unsigned int>, std::vector<double>>
400  const std::function<std::pair<bool, unsigned int>(
401  const TableIndices<dim> &)> & predicate,
402  const VectorTools::NormType norm_type,
403  const double smallest_abs_coefficient = 1e-10);
404 
410  std::pair<double, double>
411  linear_regression(const std::vector<double> &x, const std::vector<double> &y);
412 
413 } // namespace FESeries
414 
419 #ifndef DOXYGEN
420 
421 // ------------------- inline and template functions ----------------
422 
423 namespace internal
424 {
425  namespace FESeriesImplementation
426  {
427  template <int dim, typename CoefficientType>
428  void
429  fill_map_index(
430  const Table<dim, CoefficientType> &coefficients,
431  const TableIndices<dim> & ind,
432  const std::function<
433  std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
434  std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
435  {
436  const std::pair<bool, unsigned int> pred_pair = predicate(ind);
437  // don't add a value if predicate is false
438  if (pred_pair.first == false)
439  return;
440 
441  const unsigned int pred_value = pred_pair.second;
442  const CoefficientType &coeff_value = coefficients(ind);
443  // If pred_value is not in the pred_to_values map, the element will be
444  // created. Otherwise a reference to the existing element is returned.
445  pred_to_values[pred_value].push_back(coeff_value);
446  }
447 
448 
449 
450  template <typename CoefficientType>
451  void
452  fill_map(
453  const Table<1, CoefficientType> &coefficients,
454  const std::function<
455  std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
456  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
457  {
458  for (unsigned int i = 0; i < coefficients.size(0); ++i)
459  {
460  const TableIndices<1> ind(i);
461  fill_map_index(coefficients, ind, predicate, pred_to_values);
462  }
463  }
464 
465 
466 
467  template <typename CoefficientType>
468  void
469  fill_map(
470  const Table<2, CoefficientType> &coefficients,
471  const std::function<
472  std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
473  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
474  {
475  for (unsigned int i = 0; i < coefficients.size(0); ++i)
476  for (unsigned int j = 0; j < coefficients.size(1); ++j)
477  {
478  const TableIndices<2> ind(i, j);
479  fill_map_index(coefficients, ind, predicate, pred_to_values);
480  }
481  }
482 
483 
484 
485  template <typename CoefficientType>
486  void
487  fill_map(
488  const Table<3, CoefficientType> &coefficients,
489  const std::function<
490  std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
491  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
492  {
493  for (unsigned int i = 0; i < coefficients.size(0); ++i)
494  for (unsigned int j = 0; j < coefficients.size(1); ++j)
495  for (unsigned int k = 0; k < coefficients.size(2); ++k)
496  {
497  const TableIndices<3> ind(i, j, k);
498  fill_map_index(coefficients, ind, predicate, pred_to_values);
499  }
500  }
501 
502 
503 
504  template <typename Number>
505  double
506  complex_mean_value(const Number &value)
507  {
508  return value;
509  }
510 
511 
512 
513  template <typename Number>
514  double
515  complex_mean_value(const std::complex<Number> &value)
516  {
517  AssertThrow(false,
518  ExcMessage(
519  "FESeries::process_coefficients() can not be used with "
520  "complex-valued coefficients and VectorTools::mean norm."));
521  return std::abs(value);
522  }
523  } // namespace FESeriesImplementation
524 } // namespace internal
525 
526 
527 
528 template <int dim, typename CoefficientType>
529 std::pair<std::vector<unsigned int>, std::vector<double>>
531  const Table<dim, CoefficientType> &coefficients,
532  const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
533  & predicate,
534  const VectorTools::NormType norm_type,
535  const double smallest_abs_coefficient)
536 {
537  Assert(smallest_abs_coefficient >= 0.,
538  ExcMessage("smallest_abs_coefficient should be non-negative."));
539 
540  std::vector<unsigned int> predicate_values;
541  std::vector<double> norm_values;
542 
543  // first, parse all table elements into a map of predicate values and
544  // coefficients. We could have stored (predicate values ->TableIndicies) map,
545  // but its processing would have been much harder later on.
546  std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
547  internal::FESeriesImplementation::fill_map(coefficients,
548  predicate,
549  pred_to_values);
550 
551  // now go through the map and populate the @p norm_values based on @p norm:
552  for (const auto &pred_to_value : pred_to_values)
553  {
554  Vector<CoefficientType> values(pred_to_value.second.cbegin(),
555  pred_to_value.second.cend());
556 
557  double norm_value = 0;
558  switch (norm_type)
559  {
561  {
562  norm_value = values.l2_norm();
563  break;
564  }
566  {
567  norm_value = values.l1_norm();
568  break;
569  }
571  {
572  norm_value = values.linfty_norm();
573  break;
574  }
575  case VectorTools::mean:
576  {
577  norm_value = internal::FESeriesImplementation::complex_mean_value(
578  values.mean_value());
579  break;
580  }
581  default:
582  AssertThrow(false, ExcNotImplemented());
583  break;
584  }
585 
586  // will use all non-zero coefficients
587  if (std::abs(norm_value) > smallest_abs_coefficient)
588  {
589  predicate_values.push_back(pred_to_value.first);
590  norm_values.push_back(norm_value);
591  }
592  }
593 
594  return std::make_pair(predicate_values, norm_values);
595 }
596 
597 
598 
599 template <int dim, int spacedim>
600 template <class Archive>
601 inline void
603  Archive &ar,
604  const unsigned int /*version*/)
605 {
606  // Store information about those resources which have been used to generate
607  // the transformation matrices.
608  // mode vector
610 
611  // finite element collection
612  unsigned int size = fe_collection->size();
613  ar & size;
614  for (unsigned int i = 0; i < size; ++i)
615  ar &(*fe_collection)[i].get_name();
616 
617  // quadrature collection
618  size = q_collection.size();
619  ar &size;
620  for (unsigned int i = 0; i < size; ++i)
621  ar &q_collection[i];
622 
623  // Store the actual transform matrices.
625 }
626 
627 
628 
629 template <int dim, int spacedim>
630 template <class Archive>
631 inline void
633  Archive &ar,
634  const unsigned int /*version*/)
635 {
636  // Check whether the currently registered resources are compatible with
637  // the transformation matrices to load.
638  // mode vector
639  std::vector<unsigned int> compare_coefficients;
640  ar & compare_coefficients;
641  Assert(compare_coefficients == n_coefficients_per_direction,
642  ExcMessage("A different number of coefficients vector has been used "
643  "to generate the transformation matrices you are about "
644  "to load!"));
645 
646  // finite element collection
647  unsigned int size;
648  ar & size;
649  AssertDimension(size, fe_collection->size());
650  std::string name;
651  for (unsigned int i = 0; i < size; ++i)
652  {
653  ar &name;
654  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
655  ExcMessage("A different FECollection has been used to generate "
656  "the transformation matrices you are about to load!"));
657  }
658 
659  // quadrature collection
660  ar &size;
661  AssertDimension(size, q_collection.size());
662  Quadrature<dim> quadrature;
663  for (unsigned int i = 0; i < size; ++i)
664  {
665  ar &quadrature;
666  Assert(quadrature == q_collection[i],
667  ExcMessage("A different QCollection has been used to generate "
668  "the transformation matrices you are about to load!"));
669  }
670 
671  // Restore the transform matrices since all prerequisites are fulfilled.
673 }
674 
675 
676 
677 template <int dim, int spacedim>
678 template <class Archive>
679 inline void
681  Archive &ar,
682  const unsigned int /*version*/)
683 {
684  // Store information about those resources which have been used to generate
685  // the transformation matrices.
686  // mode vector
688 
689  // finite element collection
690  unsigned int size = fe_collection->size();
691  ar & size;
692  for (unsigned int i = 0; i < size; ++i)
693  ar &(*fe_collection)[i].get_name();
694 
695  // quadrature collection
696  size = q_collection.size();
697  ar &size;
698  for (unsigned int i = 0; i < size; ++i)
699  ar &q_collection[i];
700 
701  // Store the actual transform matrices.
702  ar &legendre_transform_matrices;
703 }
704 
705 
706 
707 template <int dim, int spacedim>
708 template <class Archive>
709 inline void
711  Archive &ar,
712  const unsigned int /*version*/)
713 {
714  // Check whether the currently registered resources are compatible with
715  // the transformation matrices to load.
716  // mode vector
717  std::vector<unsigned int> compare_coefficients;
718  ar & compare_coefficients;
719  Assert(compare_coefficients == n_coefficients_per_direction,
720  ExcMessage("A different number of coefficients vector has been used "
721  "to generate the transformation matrices you are about "
722  "to load!"));
723 
724  // finite element collection
725  unsigned int size;
726  ar & size;
727  AssertDimension(size, fe_collection->size());
728  std::string name;
729  for (unsigned int i = 0; i < size; ++i)
730  {
731  ar &name;
732  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
733  ExcMessage("A different FECollection has been used to generate "
734  "the transformation matrices you are about to load!"));
735  }
736 
737  // quadrature collection
738  ar &size;
739  AssertDimension(size, q_collection.size());
740  Quadrature<dim> quadrature;
741  for (unsigned int i = 0; i < size; ++i)
742  {
743  ar &quadrature;
744  Assert(quadrature == q_collection[i],
745  ExcMessage("A different QCollection has been used to generate "
746  "the transformation matrices you are about to load!"));
747  }
748 
749  // Restore the transform matrices since all prerequisites are fulfilled.
750  ar &legendre_transform_matrices;
751 }
752 
753 
754 #endif // DOXYGEN
755 
757 
758 #endif // dealii_fe_series_h
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:184
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &fourier_coefficients)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
const hp::QCollection< dim > q_collection
Definition: fe_series.h:359
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:369
std::pair< double, double > linear_regression(const std::vector< double > &x, const std::vector< double > &y)
Definition: fe_series.cc:30
void save_transformation_matrices(Archive &ar, const unsigned int version)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition: fe_series.h:364
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:349
static ::ExceptionBase & ExcMessage(std::string arg1)
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)
#define Assert(cond, exc)
Definition: exceptions.h:1461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:354
void precalculate_all_transformation_matrices()
typename std::complex< double > CoefficientType
Definition: fe_series.h:91
void load_transformation_matrices(Archive &ar, const unsigned int version)
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:179
Table< dim, Tensor< 1, dim > > k_vectors
Definition: fe_series.h:194
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)
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
const unsigned int component
Definition: fe_series.h:375
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition: fe_series.h:199
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:204
bool operator==(const Fourier< dim, spacedim > &fourier) const
const unsigned int component
Definition: fe_series.h:210
static ::ExceptionBase & ExcNotImplemented()
void save_transformation_matrices(Archive &ar, const unsigned int version)
const hp::QCollection< dim > q_collection
Definition: fe_series.h:189
double legendre(unsigned int l, double x)
Definition: cmath.h:75
void load_transformation_matrices(Archive &ar, const unsigned int version)