Reference documentation for deal.II version Git f40be01994 2020-04-09 07:13:12 +0200
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_series.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2018 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 
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/subscriptor.h>
25 #include <deal.II/base/table.h>
26 #include <deal.II/base/table_indices.h>
27 #include <deal.II/base/tensor.h>
28 
29 #include <deal.II/hp/fe_collection.h>
30 #include <deal.II/hp/q_collection.h>
31 
32 #include <deal.II/lac/full_matrix.h>
33 #include <deal.II/lac/vector.h>
34 
35 #include <deal.II/numerics/vector_tools.h>
36 
37 #include <memory>
38 #include <string>
39 #include <vector>
40 
41 
42 DEAL_II_NAMESPACE_OPEN
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 
103  Fourier(const unsigned int size_in_each_direction,
106 
112  template <typename Number>
113  void
114  calculate(const ::Vector<Number> &local_dof_values,
115  const unsigned int cell_active_fe_index,
116  Table<dim, CoefficientType> & fourier_coefficients);
117 
128  void
130 
140  template <class Archive>
141  void
142  save_transformation_matrices(Archive &ar, const unsigned int version);
143 
148  template <class Archive>
149  void
150  load_transformation_matrices(Archive &ar, const unsigned int version);
151 
155  bool
156  operator==(const Fourier<dim, spacedim> &fourier) const;
157 
158  private:
163 
168 
173 
177  std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;
178 
182  std::vector<CoefficientType> unrolled_coefficients;
183  };
184 
228  template <int dim, int spacedim = dim>
229  class Legendre : public Subscriptor
230  {
231  public:
232  using CoefficientType = double;
233 
241  Legendre(const unsigned int size_in_each_direction,
244 
250  template <typename Number>
251  void
252  calculate(const ::Vector<Number> &local_dof_values,
253  const unsigned int cell_active_fe_index,
254  Table<dim, CoefficientType> & legendre_coefficients);
255 
266  void
268 
278  template <class Archive>
279  void
280  save_transformation_matrices(Archive &ar, const unsigned int version);
281 
286  template <class Archive>
287  void
288  load_transformation_matrices(Archive &ar, const unsigned int version);
289 
293  bool
294  operator==(const Legendre<dim, spacedim> &legendre) const;
295 
296  private:
300  const unsigned int N;
301 
306 
311 
315  std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
316 
320  std::vector<CoefficientType> unrolled_coefficients;
321  };
322 
323 
338  template <int dim, typename CoefficientType>
339  std::pair<std::vector<unsigned int>, std::vector<double>>
341  const std::function<std::pair<bool, unsigned int>(
342  const TableIndices<dim> &)> & predicate,
343  const VectorTools::NormType norm);
344 
345 
346 
352  std::pair<double, double>
353  linear_regression(const std::vector<double> &x, const std::vector<double> &y);
354 
355 } // namespace FESeries
356 
359 #ifndef DOXYGEN
360 
361 // ------------------- inline and template functions ----------------
362 
363 namespace internal
364 {
365  namespace FESeriesImplementation
366  {
367  template <int dim, typename CoefficientType>
368  void
369  fill_map_index(
370  const Table<dim, CoefficientType> &coefficients,
371  const TableIndices<dim> & ind,
372  const std::function<
373  std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
374  std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
375  {
376  const std::pair<bool, unsigned int> pred_pair = predicate(ind);
377  // don't add a value if predicate is false
378  if (pred_pair.first == false)
379  return;
380 
381  const unsigned int pred_value = pred_pair.second;
382  const CoefficientType &coeff_value = coefficients(ind);
383  // If pred_value is not in the pred_to_values map, the element will be
384  // created. Otherwise a reference to the existing element is returned.
385  pred_to_values[pred_value].push_back(coeff_value);
386  }
387 
388  template <typename CoefficientType>
389  void
390  fill_map(
391  const Table<1, CoefficientType> &coefficients,
392  const std::function<
393  std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
394  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
395  {
396  for (unsigned int i = 0; i < coefficients.size(0); i++)
397  {
398  const TableIndices<1> ind(i);
399  fill_map_index(coefficients, ind, predicate, pred_to_values);
400  }
401  }
402 
403  template <typename CoefficientType>
404  void
405  fill_map(
406  const Table<2, CoefficientType> &coefficients,
407  const std::function<
408  std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
409  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
410  {
411  for (unsigned int i = 0; i < coefficients.size(0); i++)
412  for (unsigned int j = 0; j < coefficients.size(1); j++)
413  {
414  const TableIndices<2> ind(i, j);
415  fill_map_index(coefficients, ind, predicate, pred_to_values);
416  }
417  }
418 
419  template <typename CoefficientType>
420  void
421  fill_map(
422  const Table<3, CoefficientType> &coefficients,
423  const std::function<
424  std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
425  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
426  {
427  for (unsigned int i = 0; i < coefficients.size(0); i++)
428  for (unsigned int j = 0; j < coefficients.size(1); j++)
429  for (unsigned int k = 0; k < coefficients.size(2); k++)
430  {
431  const TableIndices<3> ind(i, j, k);
432  fill_map_index(coefficients, ind, predicate, pred_to_values);
433  }
434  }
435 
436 
437  template <typename Number>
438  double
439  complex_mean_value(const Number &value)
440  {
441  return value;
442  }
443 
444  template <typename Number>
445  double
446  complex_mean_value(const std::complex<Number> &value)
447  {
448  AssertThrow(false,
449  ExcMessage(
450  "FESeries::process_coefficients() can not be used with "
451  "complex-valued coefficients and VectorTools::mean norm."));
452  return std::abs(value);
453  }
454  } // namespace FESeriesImplementation
455 } // namespace internal
456 
457 
458 template <int dim, typename CoefficientType>
459 std::pair<std::vector<unsigned int>, std::vector<double>>
461  const Table<dim, CoefficientType> &coefficients,
462  const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
463  & predicate,
464  const VectorTools::NormType norm)
465 {
466  std::vector<unsigned int> predicate_values;
467  std::vector<double> norm_values;
468 
469  // first, parse all table elements into a map of predicate values and
470  // coefficients. We could have stored (predicate values ->TableIndicies) map,
471  // but its processing would have been much harder later on.
472  std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
473  internal::FESeriesImplementation::fill_map(coefficients,
474  predicate,
475  pred_to_values);
476 
477  // now go through the map and populate the @p norm_values based on @p norm:
478  for (const auto &pred_to_value : pred_to_values)
479  {
480  predicate_values.push_back(pred_to_value.first);
481  Vector<CoefficientType> values(pred_to_value.second.cbegin(),
482  pred_to_value.second.cend());
483 
484  switch (norm)
485  {
487  {
488  norm_values.push_back(values.l2_norm());
489  break;
490  }
492  {
493  norm_values.push_back(values.l1_norm());
494  break;
495  }
497  {
498  norm_values.push_back(values.linfty_norm());
499  break;
500  }
501  case VectorTools::mean:
502  {
503  norm_values.push_back(
504  internal::FESeriesImplementation::complex_mean_value(
505  values.mean_value()));
506  break;
507  }
508  default:
509  AssertThrow(false, ExcNotImplemented());
510  break;
511  }
512  }
513 
514  return std::make_pair(predicate_values, norm_values);
515 }
516 
517 
518 
519 template <int dim, int spacedim>
520 template <class Archive>
521 inline void
523  Archive &ar,
524  const unsigned int /*version*/)
525 {
526  // Store information about those resources which have been used to generate
527  // the transformation matrices.
528  // mode vector
529  unsigned int size = k_vectors.n_elements();
530  ar & size;
531 
532  // finite element collection
533  size = fe_collection->size();
534  ar &size;
535  for (unsigned int i = 0; i < size; ++i)
536  ar &(*fe_collection)[i].get_name();
537 
538  // quadrature collection
539  size = q_collection->size();
540  ar &size;
541  for (unsigned int i = 0; i < size; ++i)
542  ar &(*q_collection)[i];
543 
544  // Store the actual transform matrices.
546 }
547 
548 
549 
550 template <int dim, int spacedim>
551 template <class Archive>
552 inline void
554  Archive &ar,
555  const unsigned int /*version*/)
556 {
557  // Check whether the currently registered resources are compatible with
558  // the transformation matrices to load.
559  // mode vector
560  unsigned int size;
561  ar & size;
563 
564  // finite element collection
565  ar &size;
566  AssertDimension(size, fe_collection->size());
567  std::string name;
568  for (unsigned int i = 0; i < size; ++i)
569  {
570  ar &name;
571  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
572  ExcMessage("A different FECollection has been used to generate "
573  "the transformation matrices you are about to load!"));
574  }
575 
576  // quadrature collection
577  ar &size;
578  AssertDimension(size, q_collection->size());
579  Quadrature<dim> quadrature;
580  for (unsigned int i = 0; i < size; ++i)
581  {
582  ar &quadrature;
583  Assert(quadrature == (*q_collection)[i],
584  ExcMessage("A different QCollection has been used to generate "
585  "the transformation matrices you are about to load!"));
586  }
587 
588  // Restore the transform matrices since all prerequisites are fulfilled.
590 }
591 
592 
593 
594 template <int dim, int spacedim>
595 template <class Archive>
596 inline void
598  Archive &ar,
599  const unsigned int /*version*/)
600 {
601  // Store information about those resources which have been used to generate
602  // the transformation matrices.
603  // mode vector
604  unsigned int size = N;
605  ar & size;
606 
607  // finite element collection
608  size = fe_collection->size();
609  ar &size;
610  for (unsigned int i = 0; i < size; ++i)
611  ar &(*fe_collection)[i].get_name();
612 
613  // quadrature collection
614  size = q_collection->size();
615  ar &size;
616  for (unsigned int i = 0; i < size; ++i)
617  ar &(*q_collection)[i];
618 
619  // Store the actual transform matrices.
620  ar &legendre_transform_matrices;
621 }
622 
623 
624 
625 template <int dim, int spacedim>
626 template <class Archive>
627 inline void
629  Archive &ar,
630  const unsigned int /*version*/)
631 {
632  // Check whether the currently registered resources are compatible with
633  // the transformation matrices to load.
634  // mode vector
635  unsigned int size;
636  ar & size;
637  AssertDimension(size, N);
638 
639  // finite element collection
640  ar &size;
641  AssertDimension(size, fe_collection->size());
642  std::string name;
643  for (unsigned int i = 0; i < size; ++i)
644  {
645  ar &name;
646  Assert(name.compare((*fe_collection)[i].get_name()) == 0,
647  ExcMessage("A different FECollection has been used to generate "
648  "the transformation matrices you are about to load!"));
649  }
650 
651  // quadrature collection
652  ar &size;
653  AssertDimension(size, q_collection->size());
654  Quadrature<dim> quadrature;
655  for (unsigned int i = 0; i < size; ++i)
656  {
657  ar &quadrature;
658  Assert(quadrature == (*q_collection)[i],
659  ExcMessage("A different QCollection has been used to generate "
660  "the transformation matrices you are about to load!"));
661  }
662 
663  // Restore the transform matrices since all prerequisites are fulfilled.
664  ar &legendre_transform_matrices;
665 }
666 
667 
668 #endif // DOXYGEN
669 
670 DEAL_II_NAMESPACE_CLOSE
671 
672 #endif // dealii_fe_series_H
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:162
#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)
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:320
SmartPointer< const hp::QCollection< dim > > q_collection
Definition: fe_series.h:310
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:1531
SmartPointer< const hp::QCollection< dim > > q_collection
Definition: fe_series.h:167
size_type n_elements() const
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition: fe_series.h:315
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1419
Fourier(const unsigned int size_in_each_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection)
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:305
void precalculate_all_transformation_matrices()
const unsigned int N
Definition: fe_series.h:300
void load_transformation_matrices(Archive &ar, const unsigned int version)
Table< dim, Tensor< 1, dim > > k_vectors
Definition: fe_series.h:172
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)
size_type size(const unsigned int i) const
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition: fe_series.h:177
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:182
bool operator==(const Fourier< dim, spacedim > &fourier) const
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:699
void save_transformation_matrices(Archive &ar, const unsigned int version)
void load_transformation_matrices(Archive &ar, const unsigned int version)