Reference documentation for deal.II version Git a73d35cdfe 2020-05-28 10:34:58 +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_fourier.cc
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 
17 
18 #include <deal.II/base/numbers.h>
20 
21 #include <deal.II/fe/fe_series.h>
22 
23 #include <iostream>
24 
25 
27 
28 namespace
29 {
30  void set_k_vectors(Table<1, Tensor<1, 1>> &k_vectors, const unsigned int N)
31  {
32  k_vectors.reinit(TableIndices<1>(N));
33  for (unsigned int i = 0; i < N; ++i)
34  k_vectors(i)[0] = 2. * numbers::PI * i;
35  }
36 
37  void set_k_vectors(Table<2, Tensor<1, 2>> &k_vectors, const unsigned int N)
38  {
39  k_vectors.reinit(TableIndices<2>(N, N));
40  for (unsigned int i = 0; i < N; ++i)
41  for (unsigned int j = 0; j < N; ++j)
42  {
43  k_vectors(i, j)[0] = 2. * numbers::PI * i;
44  k_vectors(i, j)[1] = 2. * numbers::PI * j;
45  }
46  }
47 
48  void set_k_vectors(Table<3, Tensor<1, 3>> &k_vectors, const unsigned int N)
49  {
50  k_vectors.reinit(TableIndices<3>(N, N, N));
51  for (unsigned int i = 0; i < N; ++i)
52  for (unsigned int j = 0; j < N; ++j)
53  for (unsigned int k = 0; k < N; ++k)
54  {
55  k_vectors(i, j, k)[0] = 2. * numbers::PI * i;
56  k_vectors(i, j, k)[1] = 2. * numbers::PI * j;
57  k_vectors(i, j, k)[2] = 2. * numbers::PI * k;
58  }
59  }
60 
61 
62 
63  template <int dim, int spacedim>
64  std::complex<double>
65  integrate(const FiniteElement<dim, spacedim> &fe,
66  const Quadrature<dim> & quadrature,
67  const Tensor<1, dim> & k_vector,
68  const unsigned int j)
69  {
70  std::complex<double> sum = 0;
71  for (unsigned int q = 0; q < quadrature.size(); ++q)
72  {
73  const Point<dim> &x_q = quadrature.point(q);
74  sum += std::exp(std::complex<double>(0, 1) * (k_vector * x_q)) *
75  fe.shape_value(j, x_q) * quadrature.weight(q);
76  }
77  return sum;
78  }
79 
80 
81 
82  /*
83  * Ensure that the transformation matrix for FiniteElement index
84  * @p fe_index is calculated. If not, calculate it.
85  */
86  template <int spacedim>
87  void
88  ensure_existence(
89  const std::vector<unsigned int> & n_coefficients_per_direction,
90  const hp::FECollection<1, spacedim> & fe_collection,
91  const hp::QCollection<1> & q_collection,
92  const Table<1, Tensor<1, 1>> & k_vectors,
93  const unsigned int fe,
94  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
95  {
96  AssertIndexRange(fe, fe_collection.size());
97 
98  if (fourier_transform_matrices[fe].m() == 0)
99  {
100  fourier_transform_matrices[fe].reinit(n_coefficients_per_direction[fe],
101  fe_collection[fe].dofs_per_cell);
102 
103  for (unsigned int k = 0; k < n_coefficients_per_direction[fe]; ++k)
104  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
105  fourier_transform_matrices[fe](k, j) =
106  integrate(fe_collection[fe], q_collection[fe], k_vectors(k), j);
107  }
108  }
109 
110  template <int spacedim>
111  void
112  ensure_existence(
113  const std::vector<unsigned int> & n_coefficients_per_direction,
114  const hp::FECollection<2, spacedim> & fe_collection,
115  const hp::QCollection<2> & q_collection,
116  const Table<2, Tensor<1, 2>> & k_vectors,
117  const unsigned int fe,
118  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
119  {
120  AssertIndexRange(fe, fe_collection.size());
121 
122  if (fourier_transform_matrices[fe].m() == 0)
123  {
124  fourier_transform_matrices[fe].reinit(
125  Utilities::fixed_power<2>(n_coefficients_per_direction[fe]),
126  fe_collection[fe].dofs_per_cell);
127 
128  unsigned int k = 0;
129  for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
130  for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe];
131  ++k2, ++k)
132  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
133  fourier_transform_matrices[fe](k, j) = integrate(
134  fe_collection[fe], q_collection[fe], k_vectors(k1, k2), j);
135  }
136  }
137 
138  template <int spacedim>
139  void
140  ensure_existence(
141  const std::vector<unsigned int> & n_coefficients_per_direction,
142  const hp::FECollection<3, spacedim> & fe_collection,
143  const hp::QCollection<3> & q_collection,
144  const Table<3, Tensor<1, 3>> & k_vectors,
145  const unsigned int fe,
146  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
147  {
148  AssertIndexRange(fe, fe_collection.size());
149 
150  if (fourier_transform_matrices[fe].m() == 0)
151  {
152  fourier_transform_matrices[fe].reinit(
153  Utilities::fixed_power<3>(n_coefficients_per_direction[fe]),
154  fe_collection[fe].dofs_per_cell);
155 
156  unsigned int k = 0;
157  for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
158  for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe]; ++k2)
159  for (unsigned int k3 = 0; k3 < n_coefficients_per_direction[fe];
160  ++k3, ++k)
161  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
162  fourier_transform_matrices[fe](k, j) =
163  integrate(fe_collection[fe],
164  q_collection[fe],
165  k_vectors(k1, k2, k3),
166  j);
167  }
168  }
169 } // namespace
170 
171 
172 
173 namespace FESeries
174 {
175  template <int dim, int spacedim>
177  const std::vector<unsigned int> & n_coefficients_per_direction,
178  const hp::FECollection<dim, spacedim> &fe_collection,
179  const hp::QCollection<dim> & q_collection)
180  : n_coefficients_per_direction(n_coefficients_per_direction)
181  , fe_collection(&fe_collection)
182  , q_collection(q_collection)
183  , fourier_transform_matrices(fe_collection.size())
184  {
185  Assert(n_coefficients_per_direction.size() == fe_collection.size() &&
186  n_coefficients_per_direction.size() == q_collection.size(),
187  ExcMessage("All parameters are supposed to have the same size."));
188 
189  const unsigned int max_n_coefficients_per_direction =
190  *std::max_element(n_coefficients_per_direction.cbegin(),
191  n_coefficients_per_direction.cend());
192  set_k_vectors(k_vectors, max_n_coefficients_per_direction);
193 
194  // reserve sufficient memory
196  }
197 
198 
199 
200  template <int dim, int spacedim>
202  const unsigned int n_coefficients_per_direction,
203  const hp::FECollection<dim, spacedim> &fe_collection,
204  const hp::QCollection<dim> & q_collection)
205  : Fourier<dim, spacedim>(
206  std::vector<unsigned int>(fe_collection.size(),
207  n_coefficients_per_direction),
208  fe_collection,
209  q_collection)
210  {}
211 
212 
213 
214  template <int dim, int spacedim>
215  inline bool
217  operator==(const Fourier<dim, spacedim> &fourier) const
218  {
219  return (
220  (n_coefficients_per_direction == fourier.n_coefficients_per_direction) &&
221  (*fe_collection == *(fourier.fe_collection)) &&
222  (q_collection == fourier.q_collection) &&
223  (k_vectors == fourier.k_vectors) &&
225  }
226 
227 
228 
229  template <int dim, int spacedim>
230  void
232  {
233  Threads::TaskGroup<> task_group;
234  for (unsigned int fe = 0; fe < fe_collection->size(); ++fe)
235  task_group += Threads::new_task([&, fe]() {
236  ensure_existence(n_coefficients_per_direction,
237  *fe_collection,
238  q_collection,
239  k_vectors,
240  fe,
242  });
243 
244  task_group.join_all();
245  }
246 
247 
248 
249  template <int dim, int spacedim>
250  unsigned int
252  const unsigned int index) const
253  {
254  return n_coefficients_per_direction[index];
255  }
256 
257 
258 
259  template <int dim, int spacedim>
260  template <typename Number>
261  void
263  const Vector<Number> & local_dof_values,
264  const unsigned int cell_active_fe_index,
265  Table<dim, CoefficientType> &fourier_coefficients)
266  {
267  for (unsigned int d = 0; d < dim; ++d)
268  AssertDimension(fourier_coefficients.size(d),
269  n_coefficients_per_direction[cell_active_fe_index]);
270 
271  ensure_existence(n_coefficients_per_direction,
272  *fe_collection,
273  q_collection,
274  k_vectors,
275  cell_active_fe_index,
277 
279  fourier_transform_matrices[cell_active_fe_index];
280 
281  unrolled_coefficients.resize(Utilities::fixed_power<dim>(
282  n_coefficients_per_direction[cell_active_fe_index]));
283  std::fill(unrolled_coefficients.begin(),
284  unrolled_coefficients.end(),
285  CoefficientType(0.));
286 
287  Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError());
288 
289  Assert(local_dof_values.size() == matrix.n(),
290  ExcDimensionMismatch(local_dof_values.size(), matrix.n()));
291 
292  for (unsigned int i = 0; i < unrolled_coefficients.size(); i++)
293  for (unsigned int j = 0; j < local_dof_values.size(); j++)
294  unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j];
295 
296  fourier_coefficients.fill(unrolled_coefficients.begin());
297  }
298 } // namespace FESeries
299 
300 
301 // explicit instantiations
302 #include "fe_series_fourier.inst"
303 
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:190
size_type m() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
inline ::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &x)
Contents is actually a matrix.
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &fourier_coefficients)
Task< RT > new_task(const std::function< RT()> &function)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
const Point< dim > & point(const unsigned int i) const
STL namespace.
unsigned int size() const
Fourier(const std::vector< unsigned int > &n_coefficients_per_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection)
size_type n_elements() const
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1419
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
unsigned int size() const
Definition: q_collection.h:176
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:362
void precalculate_all_transformation_matrices()
typename std::complex< double > CoefficientType
Definition: fe_series.h:94
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:185
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
Table< dim, Tensor< 1, dim > > k_vectors
Definition: fe_series.h:200
unsigned int size() const
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
size_type size(const unsigned int i) const
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:158
Definition: tensor.h:450
static constexpr double PI
Definition: numbers.h:237
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
static const char N
bool operator==(const Fourier< dim, spacedim > &fourier) const
Definition: table.h:695
const hp::QCollection< dim > q_collection
Definition: fe_series.h:195
void fill(InputIterator entries, const bool C_style_indexing=true)
double weight(const unsigned int i) const
static ::ExceptionBase & ExcInternalError()