Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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 
118  private:
123 
128 
133 
137  std::vector<FullMatrix<CoefficientType>> fourier_transform_matrices;
138 
142  std::vector<CoefficientType> unrolled_coefficients;
143  };
144 
188  template <int dim, int spacedim = dim>
189  class Legendre : public Subscriptor
190  {
191  public:
192  using CoefficientType = double;
193 
201  Legendre(const unsigned int size_in_each_direction,
204 
210  template <typename Number>
211  void
212  calculate(const ::Vector<Number> &local_dof_values,
213  const unsigned int cell_active_fe_index,
214  Table<dim, CoefficientType> & legendre_coefficients);
215 
216  private:
220  const unsigned int N;
221 
226 
231 
235  std::vector<FullMatrix<CoefficientType>> legendre_transform_matrices;
236 
240  std::vector<CoefficientType> unrolled_coefficients;
241  };
242 
243 
258  template <int dim, typename CoefficientType>
259  std::pair<std::vector<unsigned int>, std::vector<double>>
261  const std::function<std::pair<bool, unsigned int>(
262  const TableIndices<dim> &)> & predicate,
263  const VectorTools::NormType norm);
264 
265 
266 
272  std::pair<double, double>
273  linear_regression(const std::vector<double> &x, const std::vector<double> &y);
274 
275 } // namespace FESeries
276 
279 #ifndef DOXYGEN
280 
281 // ------------------- inline and template functions ----------------
282 
283 namespace internal
284 {
285  namespace FESeriesImplementation
286  {
287  template <int dim, typename CoefficientType>
288  void
289  fill_map_index(
290  const Table<dim, CoefficientType> &coefficients,
291  const TableIndices<dim> & ind,
292  const std::function<
293  std::pair<bool, unsigned int>(const TableIndices<dim> &)> &predicate,
294  std::map<unsigned int, std::vector<CoefficientType>> &pred_to_values)
295  {
296  const std::pair<bool, unsigned int> pred_pair = predicate(ind);
297  // don't add a value if predicate is false
298  if (pred_pair.first == false)
299  return;
300 
301  const unsigned int pred_value = pred_pair.second;
302  const CoefficientType &coeff_value = coefficients(ind);
303  // If pred_value is not in the pred_to_values map, the element will be
304  // created. Otherwise a reference to the existing element is returned.
305  pred_to_values[pred_value].push_back(coeff_value);
306  }
307 
308  template <typename CoefficientType>
309  void
310  fill_map(
311  const Table<1, CoefficientType> &coefficients,
312  const std::function<
313  std::pair<bool, unsigned int>(const TableIndices<1> &)> &predicate,
314  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
315  {
316  for (unsigned int i = 0; i < coefficients.size(0); i++)
317  {
318  const TableIndices<1> ind(i);
319  fill_map_index(coefficients, ind, predicate, pred_to_values);
320  }
321  }
322 
323  template <typename CoefficientType>
324  void
325  fill_map(
326  const Table<2, CoefficientType> &coefficients,
327  const std::function<
328  std::pair<bool, unsigned int>(const TableIndices<2> &)> &predicate,
329  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
330  {
331  for (unsigned int i = 0; i < coefficients.size(0); i++)
332  for (unsigned int j = 0; j < coefficients.size(1); j++)
333  {
334  const TableIndices<2> ind(i, j);
335  fill_map_index(coefficients, ind, predicate, pred_to_values);
336  }
337  }
338 
339  template <typename CoefficientType>
340  void
341  fill_map(
342  const Table<3, CoefficientType> &coefficients,
343  const std::function<
344  std::pair<bool, unsigned int>(const TableIndices<3> &)> &predicate,
345  std::map<unsigned int, std::vector<CoefficientType>> & pred_to_values)
346  {
347  for (unsigned int i = 0; i < coefficients.size(0); i++)
348  for (unsigned int j = 0; j < coefficients.size(1); j++)
349  for (unsigned int k = 0; k < coefficients.size(2); k++)
350  {
351  const TableIndices<3> ind(i, j, k);
352  fill_map_index(coefficients, ind, predicate, pred_to_values);
353  }
354  }
355 
356 
357  template <typename Number>
358  double
359  complex_mean_value(const Number &value)
360  {
361  return value;
362  }
363 
364  template <typename Number>
365  double
366  complex_mean_value(const std::complex<Number> &value)
367  {
368  AssertThrow(false,
369  ExcMessage(
370  "FESeries::process_coefficients() can not be used with"
371  "complex-valued coefficients and VectorTools::mean norm."));
372  return std::abs(value);
373  }
374  } // namespace FESeriesImplementation
375 } // namespace internal
376 
377 
378 template <int dim, typename CoefficientType>
379 std::pair<std::vector<unsigned int>, std::vector<double>>
381  const Table<dim, CoefficientType> &coefficients,
382  const std::function<std::pair<bool, unsigned int>(const TableIndices<dim> &)>
383  & predicate,
384  const VectorTools::NormType norm)
385 {
386  std::vector<unsigned int> predicate_values;
387  std::vector<double> norm_values;
388 
389  // first, parse all table elements into a map of predicate values and
390  // coefficients. We could have stored (predicate values ->TableIndicies) map,
391  // but its processing would have been much harder later on.
392  std::map<unsigned int, std::vector<CoefficientType>> pred_to_values;
393  internal::FESeriesImplementation::fill_map(coefficients,
394  predicate,
395  pred_to_values);
396 
397  // now go through the map and populate the @p norm_values based on @p norm:
398  for (const auto &pred_to_value : pred_to_values)
399  {
400  predicate_values.push_back(pred_to_value.first);
401  Vector<CoefficientType> values(pred_to_value.second.cbegin(),
402  pred_to_value.second.cend());
403 
404  switch (norm)
405  {
407  {
408  norm_values.push_back(values.l2_norm());
409  break;
410  }
412  {
413  norm_values.push_back(values.l1_norm());
414  break;
415  }
417  {
418  norm_values.push_back(values.linfty_norm());
419  break;
420  }
421  case VectorTools::mean:
422  {
423  norm_values.push_back(
424  internal::FESeriesImplementation::complex_mean_value(
425  values.mean_value()));
426  break;
427  }
428  default:
429  AssertThrow(false, ExcNotImplemented());
430  break;
431  }
432  }
433 
434  return std::make_pair(predicate_values, norm_values);
435 }
436 
437 
438 #endif // DOXYGEN
439 
440 DEAL_II_NAMESPACE_CLOSE
441 
442 #endif // dealii_fe_series_H
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:122
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:240
SmartPointer< const hp::QCollection< dim > > q_collection
Definition: fe_series.h:230
Legendre(const unsigned int size_in_each_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection)
std::pair< double, double > linear_regression(const std::vector< double > &x, const std::vector< double > &y)
Definition: fe_series.cc:30
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
SmartPointer< const hp::QCollection< dim > > q_collection
Definition: fe_series.h:127
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition: fe_series.h:235
static ::ExceptionBase & ExcMessage(std::string arg1)
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:225
const unsigned int N
Definition: fe_series.h:220
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &legendre_coefficients)
Table< dim, Tensor< 1, dim > > k_vectors
Definition: fe_series.h:132
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:137
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:142
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:37