Reference documentation for deal.II version Git ca0b05a7a7 2020-12-01 11:49:21 -0500
\(\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 - 2020 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 {
86  template <int dim, int spacedim = dim>
87  class Fourier : public Subscriptor
88  {
89  public:
90  using CoefficientType = typename std::complex<double>;
91 
100  Fourier(const std::vector<unsigned int> & n_coefficients_per_direction,
103 
114  Fourier(const unsigned int n_coefficients_per_direction,
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  };
208 
209 
210 
252  template <int dim, int spacedim = dim>
253  class Legendre : public Subscriptor
254  {
255  public:
257 
266  Legendre(const std::vector<unsigned int> &n_coefficients_per_direction,
269 
279  Legendre(const unsigned int n_coefficients_per_direction,
282 
288  template <typename Number>
289  void
290  calculate(const ::Vector<Number> &local_dof_values,
291  const unsigned int cell_active_fe_index,
292  Table<dim, CoefficientType> & legendre_coefficients);
293 
298  unsigned int
299  get_n_coefficients_per_direction(const unsigned int index) const;
300 
311  void
313 
323  template <class Archive>
324  void
325  save_transformation_matrices(Archive &ar, const unsigned int version);
326 
331  template <class Archive>
332  void
333  load_transformation_matrices(Archive &ar, const unsigned int version);
334 
338  bool
339  operator==(const Legendre<dim, spacedim> &legendre) const;
340 
341  private:
346  const std::vector<unsigned int> n_coefficients_per_direction;
347 
352 
357 
361  std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
362 
366  std::vector<CoefficientType> unrolled_coefficients;
367  };
368 
369 
370 
388  template <int dim, typename CoefficientType>
389  std::pair<std::vector<unsigned int>, std::vector<double>>
391  const std::function<std::pair<bool, unsigned int>(
392  const TableIndices<dim> &)> & predicate,
393  const VectorTools::NormType norm_type,
394  const double smallest_abs_coefficient = 1e-10);
395 
401  std::pair<double, double>
402  linear_regression(const std::vector<double> &x, const std::vector<double> &y);
403 
404 } // namespace FESeries
405 
410 #ifndef DOXYGEN
411 
412 // ------------------- inline and template functions ----------------
413 
414 namespace internal
415 {
416  namespace FESeriesImplementation
417  {
418  template <int dim, typename CoefficientType>
419  void
420  fill_map_index(
421  const Table<dim, CoefficientType> &coefficients,
422  const TableIndices<dim> & ind,
423  const std::function<
424  std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
425  std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
426  {
427  const std::pair<bool, unsigned int> pred_pair = predicate(ind);
428  // don't add a value if predicate is false
429  if (pred_pair.first == false)
430  return;
431 
432  const unsigned int pred_value = pred_pair.second;
433  const CoefficientType &coeff_value = coefficients(ind);
434  // If pred_value is not in the pred_to_values map, the element will be
435  // created. Otherwise a reference to the existing element is returned.
436  pred_to_values[pred_value].push_back(coeff_value);
437  }
438 
439 
440 
441  template <typename CoefficientType>
442  void
443  fill_map(
444  const Table<1, CoefficientType> &coefficients,
445  const std::function<
446  std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
447  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
448  {
449  for (unsigned int i = 0; i < coefficients.size(0); i++)
450  {
451  const TableIndices<1> ind(i);
452  fill_map_index(coefficients, ind, predicate, pred_to_values);
453  }
454  }
455 
456 
457 
458  template <typename CoefficientType>
459  void
460  fill_map(
461  const Table<2, CoefficientType> &coefficients,
462  const std::function<
463  std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
464  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
465  {
466  for (unsigned int i = 0; i < coefficients.size(0); i++)
467  for (unsigned int j = 0; j < coefficients.size(1); j++)
468  {
469  const TableIndices<2> ind(i, j);
470  fill_map_index(coefficients, ind, predicate, pred_to_values);
471  }
472  }
473 
474 
475 
476  template <typename CoefficientType>
477  void
478  fill_map(
479  const Table<3, CoefficientType> &coefficients,
480  const std::function<
481  std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
482  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
483  {
484  for (unsigned int i = 0; i < coefficients.size(0); i++)
485  for (unsigned int j = 0; j < coefficients.size(1); j++)
486  for (unsigned int k = 0; k < coefficients.size(2); k++)
487  {
488  const TableIndices<3> ind(i, j, k);
489  fill_map_index(coefficients, ind, predicate, pred_to_values);
490  }
491  }
492 
493 
494 
495  template <typename Number>
496  double
497  complex_mean_value(const Number &value)
498  {
499  return value;
500  }
501 
502 
503 
504  template <typename Number>
505  double
506  complex_mean_value(const std::complex<Number> &value)
507  {
508  AssertThrow(false,
509  ExcMessage(
510  "FESeries::process_coefficients() can not be used with "
511  "complex-valued coefficients and VectorTools::mean norm."));
512  return std::abs(value);
513  }
514  } // namespace FESeriesImplementation
515 } // namespace internal
516 
517 
518 
519 template <int dim, typename CoefficientType>
520 std::pair<std::vector<unsigned int>, std::vector<double>>
522  const Table<dim, CoefficientType> &coefficients,
523  const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
524  & predicate,
525  const VectorTools::NormType norm_type,
526  const double smallest_abs_coefficient)
527 {
528  Assert(smallest_abs_coefficient >= 0.,
529  ExcMessage("smallest_abs_coefficient should be non-negative."));
530 
531  std::vector<unsigned int> predicate_values;
532  std::vector<double> norm_values;
533 
534  // first, parse all table elements into a map of predicate values and
535  // coefficients. We could have stored (predicate values ->TableIndicies) map,
536  // but its processing would have been much harder later on.
537  std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
538  internal::FESeriesImplementation::fill_map(coefficients,
539  predicate,
540  pred_to_values);
541 
542  // now go through the map and populate the @p norm_values based on @p norm:
543  for (const auto &pred_to_value : pred_to_values)
544  {
545  Vector<CoefficientType> values(pred_to_value.second.cbegin(),
546  pred_to_value.second.cend());
547 
548  double norm_value = 0;
549  switch (norm_type)
550  {
552  {
553  norm_value = values.l2_norm();
554  break;
555  }
557  {
558  norm_value = values.l1_norm();
559  break;
560  }
562  {
563  norm_value = values.linfty_norm();
564  break;
565  }
566  case VectorTools::mean:
567  {
568  norm_value = internal::FESeriesImplementation::complex_mean_value(
569  values.mean_value());
570  break;
571  }
572  default:
573  AssertThrow(false, ExcNotImplemented());
574  break;
575  }
576 
577  // will use all non-zero coefficients
578  if (std::abs(norm_value) > smallest_abs_coefficient)
579  {
580  predicate_values.push_back(pred_to_value.first);
581  norm_values.push_back(norm_value);
582  }
583  }
584 
585  return std::make_pair(predicate_values, norm_values);
586 }
587 
588 
589 
590 template <int dim, int spacedim>
591 template <class Archive>
592 inline void
594  Archive &ar,
595  const unsigned int /*version*/)
596 {
597  // Store information about those resources which have been used to generate
598  // the transformation matrices.
599  // mode vector
601 
602  // finite element collection
603  unsigned int size = fe_collection->size();
604  ar & size;
605  for (unsigned int i = 0; i < size; ++i)
606  ar &(*fe_collection)[i].get_name();
607 
608  // quadrature collection
609  size = q_collection.size();
610  ar &size;
611  for (unsigned int i = 0; i < size; ++i)
612  ar &q_collection[i];
613 
614  // Store the actual transform matrices.
616 }
617 
618 
619 
620 template <int dim, int spacedim>
621 template <class Archive>
622 inline void
624  Archive &ar,
625  const unsigned int /*version*/)
626 {
627  // Check whether the currently registered resources are compatible with
628  // the transformation matrices to load.
629  // mode vector
630  std::vector<unsigned int> compare_coefficients;
631  ar & compare_coefficients;
632  Assert(compare_coefficients == n_coefficients_per_direction,
633  ExcMessage("A different number of coefficients vector has been used "
634  "to generate the transformation matrices you are about "
635  "to load!"));
636 
637  // finite element collection
638  unsigned int size;
639  ar & size;
640  AssertDimension(size, fe_collection->size());
641  std::string name;
642  for (unsigned int i = 0; i < size; ++i)
643  {
644  ar &name;
645  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
646  ExcMessage("A different FECollection has been used to generate "
647  "the transformation matrices you are about to load!"));
648  }
649 
650  // quadrature collection
651  ar &size;
652  AssertDimension(size, q_collection.size());
653  Quadrature<dim> quadrature;
654  for (unsigned int i = 0; i < size; ++i)
655  {
656  ar &quadrature;
657  Assert(quadrature == q_collection[i],
658  ExcMessage("A different QCollection has been used to generate "
659  "the transformation matrices you are about to load!"));
660  }
661 
662  // Restore the transform matrices since all prerequisites are fulfilled.
664 }
665 
666 
667 
668 template <int dim, int spacedim>
669 template <class Archive>
670 inline void
672  Archive &ar,
673  const unsigned int /*version*/)
674 {
675  // Store information about those resources which have been used to generate
676  // the transformation matrices.
677  // mode vector
679 
680  // finite element collection
681  unsigned int size = fe_collection->size();
682  ar & size;
683  for (unsigned int i = 0; i < size; ++i)
684  ar &(*fe_collection)[i].get_name();
685 
686  // quadrature collection
687  size = q_collection.size();
688  ar &size;
689  for (unsigned int i = 0; i < size; ++i)
690  ar &q_collection[i];
691 
692  // Store the actual transform matrices.
693  ar &legendre_transform_matrices;
694 }
695 
696 
697 
698 template <int dim, int spacedim>
699 template <class Archive>
700 inline void
702  Archive &ar,
703  const unsigned int /*version*/)
704 {
705  // Check whether the currently registered resources are compatible with
706  // the transformation matrices to load.
707  // mode vector
708  std::vector<unsigned int> compare_coefficients;
709  ar & compare_coefficients;
710  Assert(compare_coefficients == n_coefficients_per_direction,
711  ExcMessage("A different number of coefficients vector has been used "
712  "to generate the transformation matrices you are about "
713  "to load!"));
714 
715  // finite element collection
716  unsigned int size;
717  ar & size;
718  AssertDimension(size, fe_collection->size());
719  std::string name;
720  for (unsigned int i = 0; i < size; ++i)
721  {
722  ar &name;
723  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
724  ExcMessage("A different FECollection has been used to generate "
725  "the transformation matrices you are about to load!"));
726  }
727 
728  // quadrature collection
729  ar &size;
730  AssertDimension(size, q_collection.size());
731  Quadrature<dim> quadrature;
732  for (unsigned int i = 0; i < size; ++i)
733  {
734  ar &quadrature;
735  Assert(quadrature == q_collection[i],
736  ExcMessage("A different QCollection has been used to generate "
737  "the transformation matrices you are about to load!"));
738  }
739 
740  // Restore the transform matrices since all prerequisites are fulfilled.
741  ar &legendre_transform_matrices;
742 }
743 
744 
745 #endif // DOXYGEN
746 
748 
749 #endif // dealii_fe_series_h
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:186
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
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:356
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:366
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)
Fourier(const std::vector< unsigned int > &n_coefficients_per_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition: fe_series.h:361
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:346
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1466
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:351
void precalculate_all_transformation_matrices()
typename std::complex< double > CoefficientType
Definition: fe_series.h:90
void load_transformation_matrices(Archive &ar, const unsigned int version)
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:181
Table< dim, Tensor< 1, dim > > k_vectors
Definition: fe_series.h:196
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
size_type size(const unsigned int i) const
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition: fe_series.h:201
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:206
bool operator==(const Fourier< dim, spacedim > &fourier) const
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:687
void save_transformation_matrices(Archive &ar, const unsigned int version)
const hp::QCollection< dim > q_collection
Definition: fe_series.h:191
#define DEAL_II_DEPRECATED
Definition: config.h:152
void load_transformation_matrices(Archive &ar, const unsigned int version)