Reference documentation for deal.II version Git 423cd11810 2020-05-28 01:07:25 +0200
\(\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 
57 namespace FESeries
58 {
90  template <int dim, int spacedim = dim>
91  class Fourier : public Subscriptor
92  {
93  public:
94  using CoefficientType = typename std::complex<double>;
95 
104  Fourier(const std::vector<unsigned int> & n_coefficients_per_direction,
107 
118  Fourier(const unsigned int n_coefficients_per_direction,
121 
127  template <typename Number>
128  void
129  calculate(const ::Vector<Number> &local_dof_values,
130  const unsigned int cell_active_fe_index,
131  Table<dim, CoefficientType> & fourier_coefficients);
132 
137  unsigned int
138  get_n_coefficients_per_direction(const unsigned int index) const;
139 
150  void
152 
162  template <class Archive>
163  void
164  save_transformation_matrices(Archive &ar, const unsigned int version);
165 
170  template <class Archive>
171  void
172  load_transformation_matrices(Archive &ar, const unsigned int version);
173 
177  bool
178  operator==(const Fourier<dim, spacedim> &fourier) const;
179 
180  private:
185  const std::vector<unsigned int> n_coefficients_per_direction;
186 
191 
196 
201 
205  std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;
206 
210  std::vector<CoefficientType> unrolled_coefficients;
211  };
212 
213 
214 
258  template <int dim, int spacedim = dim>
259  class Legendre : public Subscriptor
260  {
261  public:
263 
272  Legendre(const std::vector<unsigned int> &n_coefficients_per_direction,
275 
285  Legendre(const unsigned int n_coefficients_per_direction,
288 
294  template <typename Number>
295  void
296  calculate(const ::Vector<Number> &local_dof_values,
297  const unsigned int cell_active_fe_index,
298  Table<dim, CoefficientType> & legendre_coefficients);
299 
304  unsigned int
305  get_n_coefficients_per_direction(const unsigned int index) const;
306 
317  void
319 
329  template <class Archive>
330  void
331  save_transformation_matrices(Archive &ar, const unsigned int version);
332 
337  template <class Archive>
338  void
339  load_transformation_matrices(Archive &ar, const unsigned int version);
340 
344  bool
345  operator==(const Legendre<dim, spacedim> &legendre) const;
346 
347  private:
352  const std::vector<unsigned int> n_coefficients_per_direction;
353 
358 
363 
367  std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
368 
372  std::vector<CoefficientType> unrolled_coefficients;
373  };
374 
375 
376 
394  template <int dim, typename CoefficientType>
395  std::pair<std::vector<unsigned int>, std::vector<double>>
397  const std::function<std::pair<bool, unsigned int>(
398  const TableIndices<dim> &)> & predicate,
399  const VectorTools::NormType norm_type,
400  const double smallest_abs_coefficient = 1e-10);
401 
407  std::pair<double, double>
408  linear_regression(const std::vector<double> &x, const std::vector<double> &y);
409 
410 } // namespace FESeries
411 
416 #ifndef DOXYGEN
417 
418 // ------------------- inline and template functions ----------------
419 
420 namespace internal
421 {
422  namespace FESeriesImplementation
423  {
424  template <int dim, typename CoefficientType>
425  void
426  fill_map_index(
427  const Table<dim, CoefficientType> &coefficients,
428  const TableIndices<dim> & ind,
429  const std::function<
430  std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
431  std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
432  {
433  const std::pair<bool, unsigned int> pred_pair = predicate(ind);
434  // don't add a value if predicate is false
435  if (pred_pair.first == false)
436  return;
437 
438  const unsigned int pred_value = pred_pair.second;
439  const CoefficientType &coeff_value = coefficients(ind);
440  // If pred_value is not in the pred_to_values map, the element will be
441  // created. Otherwise a reference to the existing element is returned.
442  pred_to_values[pred_value].push_back(coeff_value);
443  }
444 
445 
446 
447  template <typename CoefficientType>
448  void
449  fill_map(
450  const Table<1, CoefficientType> &coefficients,
451  const std::function<
452  std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
453  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
454  {
455  for (unsigned int i = 0; i < coefficients.size(0); i++)
456  {
457  const TableIndices<1> ind(i);
458  fill_map_index(coefficients, ind, predicate, pred_to_values);
459  }
460  }
461 
462 
463 
464  template <typename CoefficientType>
465  void
466  fill_map(
467  const Table<2, CoefficientType> &coefficients,
468  const std::function<
469  std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
470  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
471  {
472  for (unsigned int i = 0; i < coefficients.size(0); i++)
473  for (unsigned int j = 0; j < coefficients.size(1); j++)
474  {
475  const TableIndices<2> ind(i, j);
476  fill_map_index(coefficients, ind, predicate, pred_to_values);
477  }
478  }
479 
480 
481 
482  template <typename CoefficientType>
483  void
484  fill_map(
485  const Table<3, CoefficientType> &coefficients,
486  const std::function<
487  std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
488  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
489  {
490  for (unsigned int i = 0; i < coefficients.size(0); i++)
491  for (unsigned int j = 0; j < coefficients.size(1); j++)
492  for (unsigned int k = 0; k < coefficients.size(2); k++)
493  {
494  const TableIndices<3> ind(i, j, k);
495  fill_map_index(coefficients, ind, predicate, pred_to_values);
496  }
497  }
498 
499 
500 
501  template <typename Number>
502  double
503  complex_mean_value(const Number &value)
504  {
505  return value;
506  }
507 
508 
509 
510  template <typename Number>
511  double
512  complex_mean_value(const std::complex<Number> &value)
513  {
514  AssertThrow(false,
515  ExcMessage(
516  "FESeries::process_coefficients() can not be used with "
517  "complex-valued coefficients and VectorTools::mean norm."));
518  return std::abs(value);
519  }
520  } // namespace FESeriesImplementation
521 } // namespace internal
522 
523 
524 
525 template <int dim, typename CoefficientType>
526 std::pair<std::vector<unsigned int>, std::vector<double>>
528  const Table<dim, CoefficientType> &coefficients,
529  const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
530  & predicate,
531  const VectorTools::NormType norm_type,
532  const double smallest_abs_coefficient)
533 {
534  Assert(smallest_abs_coefficient >= 0.,
535  ExcMessage("smallest_abs_coefficient should be non-negative."));
536 
537  std::vector<unsigned int> predicate_values;
538  std::vector<double> norm_values;
539 
540  // first, parse all table elements into a map of predicate values and
541  // coefficients. We could have stored (predicate values ->TableIndicies) map,
542  // but its processing would have been much harder later on.
543  std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
544  internal::FESeriesImplementation::fill_map(coefficients,
545  predicate,
546  pred_to_values);
547 
548  // now go through the map and populate the @p norm_values based on @p norm:
549  for (const auto &pred_to_value : pred_to_values)
550  {
551  Vector<CoefficientType> values(pred_to_value.second.cbegin(),
552  pred_to_value.second.cend());
553 
554  double norm_value = 0;
555  switch (norm_type)
556  {
558  {
559  norm_value = values.l2_norm();
560  break;
561  }
563  {
564  norm_value = values.l1_norm();
565  break;
566  }
568  {
569  norm_value = values.linfty_norm();
570  break;
571  }
572  case VectorTools::mean:
573  {
574  norm_value = internal::FESeriesImplementation::complex_mean_value(
575  values.mean_value());
576  break;
577  }
578  default:
579  AssertThrow(false, ExcNotImplemented());
580  break;
581  }
582 
583  // will use all non-zero coefficients
584  if (std::abs(norm_value) > smallest_abs_coefficient)
585  {
586  predicate_values.push_back(pred_to_value.first);
587  norm_values.push_back(norm_value);
588  }
589  }
590 
591  return std::make_pair(predicate_values, norm_values);
592 }
593 
594 
595 
596 template <int dim, int spacedim>
597 template <class Archive>
598 inline void
600  Archive &ar,
601  const unsigned int /*version*/)
602 {
603  // Store information about those resources which have been used to generate
604  // the transformation matrices.
605  // mode vector
607 
608  // finite element collection
609  unsigned int size = fe_collection->size();
610  ar & size;
611  for (unsigned int i = 0; i < size; ++i)
612  ar &(*fe_collection)[i].get_name();
613 
614  // quadrature collection
615  size = q_collection.size();
616  ar &size;
617  for (unsigned int i = 0; i < size; ++i)
618  ar &q_collection[i];
619 
620  // Store the actual transform matrices.
622 }
623 
624 
625 
626 template <int dim, int spacedim>
627 template <class Archive>
628 inline void
630  Archive &ar,
631  const unsigned int /*version*/)
632 {
633  // Check whether the currently registered resources are compatible with
634  // the transformation matrices to load.
635  // mode vector
636  std::vector<unsigned int> compare_coefficients;
637  ar & compare_coefficients;
638  Assert(compare_coefficients == n_coefficients_per_direction,
639  ExcMessage("A different number of coefficients vector has been used "
640  "to generate the transformation matrices you are about "
641  "to load!"));
642 
643  // finite element collection
644  unsigned int size;
645  ar & size;
646  AssertDimension(size, fe_collection->size());
647  std::string name;
648  for (unsigned int i = 0; i < size; ++i)
649  {
650  ar &name;
651  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
652  ExcMessage("A different FECollection has been used to generate "
653  "the transformation matrices you are about to load!"));
654  }
655 
656  // quadrature collection
657  ar &size;
658  AssertDimension(size, q_collection.size());
659  Quadrature<dim> quadrature;
660  for (unsigned int i = 0; i < size; ++i)
661  {
662  ar &quadrature;
663  Assert(quadrature == q_collection[i],
664  ExcMessage("A different QCollection has been used to generate "
665  "the transformation matrices you are about to load!"));
666  }
667 
668  // Restore the transform matrices since all prerequisites are fulfilled.
670 }
671 
672 
673 
674 template <int dim, int spacedim>
675 template <class Archive>
676 inline void
678  Archive &ar,
679  const unsigned int /*version*/)
680 {
681  // Store information about those resources which have been used to generate
682  // the transformation matrices.
683  // mode vector
685 
686  // finite element collection
687  unsigned int size = fe_collection->size();
688  ar & size;
689  for (unsigned int i = 0; i < size; ++i)
690  ar &(*fe_collection)[i].get_name();
691 
692  // quadrature collection
693  size = q_collection.size();
694  ar &size;
695  for (unsigned int i = 0; i < size; ++i)
696  ar &q_collection[i];
697 
698  // Store the actual transform matrices.
699  ar &legendre_transform_matrices;
700 }
701 
702 
703 
704 template <int dim, int spacedim>
705 template <class Archive>
706 inline void
708  Archive &ar,
709  const unsigned int /*version*/)
710 {
711  // Check whether the currently registered resources are compatible with
712  // the transformation matrices to load.
713  // mode vector
714  std::vector<unsigned int> compare_coefficients;
715  ar & compare_coefficients;
716  Assert(compare_coefficients == n_coefficients_per_direction,
717  ExcMessage("A different number of coefficients vector has been used "
718  "to generate the transformation matrices you are about "
719  "to load!"));
720 
721  // finite element collection
722  unsigned int size;
723  ar & size;
724  AssertDimension(size, fe_collection->size());
725  std::string name;
726  for (unsigned int i = 0; i < size; ++i)
727  {
728  ar &name;
729  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
730  ExcMessage("A different FECollection has been used to generate "
731  "the transformation matrices you are about to load!"));
732  }
733 
734  // quadrature collection
735  ar &size;
736  AssertDimension(size, q_collection.size());
737  Quadrature<dim> quadrature;
738  for (unsigned int i = 0; i < size; ++i)
739  {
740  ar &quadrature;
741  Assert(quadrature == q_collection[i],
742  ExcMessage("A different QCollection has been used to generate "
743  "the transformation matrices you are about to load!"));
744  }
745 
746  // Restore the transform matrices since all prerequisites are fulfilled.
747  ar &legendre_transform_matrices;
748 }
749 
750 
751 #endif // DOXYGEN
752 
754 
755 #endif // dealii_fe_series_h
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:190
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
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:362
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:372
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:1531
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition: fe_series.h:367
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:352
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1419
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:362
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:357
void precalculate_all_transformation_matrices()
typename std::complex< double > CoefficientType
Definition: fe_series.h:94
void load_transformation_matrices(Archive &ar, const unsigned int version)
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:185
Table< dim, Tensor< 1, dim > > k_vectors
Definition: fe_series.h:200
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:205
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:361
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:210
bool operator==(const Fourier< dim, spacedim > &fourier) const
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:695
void save_transformation_matrices(Archive &ar, const unsigned int version)
const hp::QCollection< dim > q_collection
Definition: fe_series.h:195
static const bool value
#define DEAL_II_DEPRECATED
Definition: config.h:151
void load_transformation_matrices(Archive &ar, const unsigned int version)