Reference documentation for deal.II version 9.3.0
\(\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 - 2021 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  const unsigned int component)
70  {
71  std::complex<double> sum = 0;
72  for (unsigned int q = 0; q < quadrature.size(); ++q)
73  {
74  const Point<dim> &x_q = quadrature.point(q);
75  sum += std::exp(std::complex<double>(0, 1) * (k_vector * x_q)) *
76  fe.shape_value_component(j, x_q, component) *
77  quadrature.weight(q);
78  }
79  return sum;
80  }
81 
82 
83 
84  /*
85  * Ensure that the transformation matrix for FiniteElement index
86  * @p fe_index is calculated. If not, calculate it.
87  */
88  template <int spacedim>
89  void
90  ensure_existence(
91  const std::vector<unsigned int> & n_coefficients_per_direction,
92  const hp::FECollection<1, spacedim> & fe_collection,
93  const hp::QCollection<1> & q_collection,
94  const Table<1, Tensor<1, 1>> & k_vectors,
95  const unsigned int fe,
96  const unsigned int component,
97  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
98  {
99  AssertIndexRange(fe, fe_collection.size());
100 
101  if (fourier_transform_matrices[fe].m() == 0)
102  {
103  fourier_transform_matrices[fe].reinit(
104  n_coefficients_per_direction[fe],
105  fe_collection[fe].n_dofs_per_cell());
106 
107  for (unsigned int k = 0; k < n_coefficients_per_direction[fe]; ++k)
108  for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell(); ++j)
109  fourier_transform_matrices[fe](k, j) = integrate(
110  fe_collection[fe], q_collection[fe], k_vectors(k), j, component);
111  }
112  }
113 
114  template <int spacedim>
115  void
116  ensure_existence(
117  const std::vector<unsigned int> & n_coefficients_per_direction,
118  const hp::FECollection<2, spacedim> & fe_collection,
119  const hp::QCollection<2> & q_collection,
120  const Table<2, Tensor<1, 2>> & k_vectors,
121  const unsigned int fe,
122  const unsigned int component,
123  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
124  {
125  AssertIndexRange(fe, fe_collection.size());
126 
127  if (fourier_transform_matrices[fe].m() == 0)
128  {
129  fourier_transform_matrices[fe].reinit(
130  Utilities::fixed_power<2>(n_coefficients_per_direction[fe]),
131  fe_collection[fe].n_dofs_per_cell());
132 
133  unsigned int k = 0;
134  for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
135  for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe];
136  ++k2, ++k)
137  for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
138  ++j)
139  fourier_transform_matrices[fe](k, j) =
140  integrate(fe_collection[fe],
141  q_collection[fe],
142  k_vectors(k1, k2),
143  j,
144  component);
145  }
146  }
147 
148  template <int spacedim>
149  void
150  ensure_existence(
151  const std::vector<unsigned int> & n_coefficients_per_direction,
152  const hp::FECollection<3, spacedim> & fe_collection,
153  const hp::QCollection<3> & q_collection,
154  const Table<3, Tensor<1, 3>> & k_vectors,
155  const unsigned int fe,
156  const unsigned int component,
157  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
158  {
159  AssertIndexRange(fe, fe_collection.size());
160 
161  if (fourier_transform_matrices[fe].m() == 0)
162  {
163  fourier_transform_matrices[fe].reinit(
164  Utilities::fixed_power<3>(n_coefficients_per_direction[fe]),
165  fe_collection[fe].n_dofs_per_cell());
166 
167  unsigned int k = 0;
168  for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
169  for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe]; ++k2)
170  for (unsigned int k3 = 0; k3 < n_coefficients_per_direction[fe];
171  ++k3, ++k)
172  for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
173  ++j)
174  fourier_transform_matrices[fe](k, j) =
175  integrate(fe_collection[fe],
176  q_collection[fe],
177  k_vectors(k1, k2, k3),
178  j,
179  component);
180  }
181  }
182 } // namespace
183 
184 
185 
186 namespace FESeries
187 {
188  template <int dim, int spacedim>
190  const std::vector<unsigned int> & n_coefficients_per_direction,
191  const hp::FECollection<dim, spacedim> &fe_collection,
192  const hp::QCollection<dim> & q_collection,
193  const unsigned int component_)
194  : n_coefficients_per_direction(n_coefficients_per_direction)
195  , fe_collection(&fe_collection)
196  , q_collection(q_collection)
197  , fourier_transform_matrices(fe_collection.size())
198  , component(component_ != numbers::invalid_unsigned_int ? component_ : 0)
199  {
200  Assert(n_coefficients_per_direction.size() == fe_collection.size() &&
201  n_coefficients_per_direction.size() == q_collection.size(),
202  ExcMessage("All parameters are supposed to have the same size."));
203 
204  if (fe_collection[0].n_components() > 1)
205  Assert(
206  component_ != numbers::invalid_unsigned_int,
207  ExcMessage(
208  "For vector-valued problems, you need to explicitly specify for "
209  "which vector component you will want to do a Fourier decomposition "
210  "by setting the 'component' argument of this constructor."));
211 
212  AssertIndexRange(component, fe_collection[0].n_components());
213 
214  const unsigned int max_n_coefficients_per_direction =
215  *std::max_element(n_coefficients_per_direction.cbegin(),
216  n_coefficients_per_direction.cend());
217  set_k_vectors(k_vectors, max_n_coefficients_per_direction);
218 
219  // reserve sufficient memory
220  unrolled_coefficients.reserve(k_vectors.n_elements());
221  }
222 
223 
224 
225  template <int dim, int spacedim>
227  const unsigned int n_coefficients_per_direction,
228  const hp::FECollection<dim, spacedim> &fe_collection,
229  const hp::QCollection<dim> & q_collection)
230  : Fourier<dim, spacedim>(
231  std::vector<unsigned int>(fe_collection.size(),
232  n_coefficients_per_direction),
233  fe_collection,
234  q_collection)
235  {}
236 
237 
238 
239  template <int dim, int spacedim>
240  inline bool
242  operator==(const Fourier<dim, spacedim> &fourier) const
243  {
244  return (
245  (n_coefficients_per_direction == fourier.n_coefficients_per_direction) &&
246  (*fe_collection == *(fourier.fe_collection)) &&
247  (q_collection == fourier.q_collection) &&
248  (k_vectors == fourier.k_vectors) &&
250  (component == fourier.component));
251  }
252 
253 
254 
255  template <int dim, int spacedim>
256  void
258  {
259  Threads::TaskGroup<> task_group;
260  for (unsigned int fe = 0; fe < fe_collection->size(); ++fe)
261  task_group += Threads::new_task([&, fe]() {
262  ensure_existence(n_coefficients_per_direction,
263  *fe_collection,
264  q_collection,
265  k_vectors,
266  fe,
267  component,
269  });
270 
271  task_group.join_all();
272  }
273 
274 
275 
276  template <int dim, int spacedim>
277  unsigned int
279  const unsigned int index) const
280  {
281  return n_coefficients_per_direction[index];
282  }
283 
284 
285 
286  template <int dim, int spacedim>
287  template <typename Number>
288  void
290  const Vector<Number> & local_dof_values,
291  const unsigned int cell_active_fe_index,
292  Table<dim, CoefficientType> &fourier_coefficients)
293  {
294  for (unsigned int d = 0; d < dim; ++d)
295  AssertDimension(fourier_coefficients.size(d),
296  n_coefficients_per_direction[cell_active_fe_index]);
297 
298  ensure_existence(n_coefficients_per_direction,
299  *fe_collection,
300  q_collection,
301  k_vectors,
302  cell_active_fe_index,
303  component,
305 
307  fourier_transform_matrices[cell_active_fe_index];
308 
309  unrolled_coefficients.resize(Utilities::fixed_power<dim>(
310  n_coefficients_per_direction[cell_active_fe_index]));
311  std::fill(unrolled_coefficients.begin(),
312  unrolled_coefficients.end(),
313  CoefficientType(0.));
314 
315  Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError());
316 
317  Assert(local_dof_values.size() == matrix.n(),
318  ExcDimensionMismatch(local_dof_values.size(), matrix.n()));
319 
320  for (unsigned int i = 0; i < unrolled_coefficients.size(); i++)
321  for (unsigned int j = 0; j < local_dof_values.size(); j++)
322  unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j];
323 
324  fourier_coefficients.fill(unrolled_coefficients.begin());
325  }
326 } // namespace FESeries
327 
328 
329 // explicit instantiations
330 #include "fe_series_fourier.inst"
331 
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:198
size_type m() const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
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)
virtual double shape_value_component(const unsigned int i, const Point< dim > &p, const unsigned int component) const
Definition: fe.cc:197
Task< RT > new_task(const std::function< RT()> &function)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
const Point< dim > & point(const unsigned int i) const
STL namespace.
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
Fourier(const std::vector< unsigned int > &n_coefficients_per_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection, const unsigned int component=numbers::invalid_unsigned_int)
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
void precalculate_all_transformation_matrices()
typename std::complex< double > CoefficientType
Definition: fe_series.h:91
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:193
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:208
unsigned int size() const
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
Definition: tensor.h:462
static constexpr double PI
Definition: numbers.h:231
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition: fe_series.h:213
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:218
static const char N
bool operator==(const Fourier< dim, spacedim > &fourier) const
const unsigned int component
Definition: fe_series.h:224
const hp::QCollection< dim > q_collection
Definition: fe_series.h:203
double weight(const unsigned int i) const
static ::ExceptionBase & ExcInternalError()