Reference documentation for deal.II version Git 8596a7cd07 2020-12-04 07:30:43 +0100
\(\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(
101  n_coefficients_per_direction[fe],
102  fe_collection[fe].n_dofs_per_cell());
103 
104  for (unsigned int k = 0; k < n_coefficients_per_direction[fe]; ++k)
105  for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell(); ++j)
106  fourier_transform_matrices[fe](k, j) =
107  integrate(fe_collection[fe], q_collection[fe], k_vectors(k), j);
108  }
109  }
110 
111  template <int spacedim>
112  void
113  ensure_existence(
114  const std::vector<unsigned int> & n_coefficients_per_direction,
115  const hp::FECollection<2, spacedim> & fe_collection,
116  const hp::QCollection<2> & q_collection,
117  const Table<2, Tensor<1, 2>> & k_vectors,
118  const unsigned int fe,
119  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
120  {
121  AssertIndexRange(fe, fe_collection.size());
122 
123  if (fourier_transform_matrices[fe].m() == 0)
124  {
125  fourier_transform_matrices[fe].reinit(
126  Utilities::fixed_power<2>(n_coefficients_per_direction[fe]),
127  fe_collection[fe].n_dofs_per_cell());
128 
129  unsigned int k = 0;
130  for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
131  for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe];
132  ++k2, ++k)
133  for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
134  ++j)
135  fourier_transform_matrices[fe](k, j) = integrate(
136  fe_collection[fe], q_collection[fe], k_vectors(k1, k2), j);
137  }
138  }
139 
140  template <int spacedim>
141  void
142  ensure_existence(
143  const std::vector<unsigned int> & n_coefficients_per_direction,
144  const hp::FECollection<3, spacedim> & fe_collection,
145  const hp::QCollection<3> & q_collection,
146  const Table<3, Tensor<1, 3>> & k_vectors,
147  const unsigned int fe,
148  std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
149  {
150  AssertIndexRange(fe, fe_collection.size());
151 
152  if (fourier_transform_matrices[fe].m() == 0)
153  {
154  fourier_transform_matrices[fe].reinit(
155  Utilities::fixed_power<3>(n_coefficients_per_direction[fe]),
156  fe_collection[fe].n_dofs_per_cell());
157 
158  unsigned int k = 0;
159  for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
160  for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe]; ++k2)
161  for (unsigned int k3 = 0; k3 < n_coefficients_per_direction[fe];
162  ++k3, ++k)
163  for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
164  ++j)
165  fourier_transform_matrices[fe](k, j) =
166  integrate(fe_collection[fe],
167  q_collection[fe],
168  k_vectors(k1, k2, k3),
169  j);
170  }
171  }
172 } // namespace
173 
174 
175 
176 namespace FESeries
177 {
178  template <int dim, int spacedim>
180  const std::vector<unsigned int> & n_coefficients_per_direction,
181  const hp::FECollection<dim, spacedim> &fe_collection,
182  const hp::QCollection<dim> & q_collection)
183  : n_coefficients_per_direction(n_coefficients_per_direction)
184  , fe_collection(&fe_collection)
185  , q_collection(q_collection)
186  , fourier_transform_matrices(fe_collection.size())
187  {
188  Assert(n_coefficients_per_direction.size() == fe_collection.size() &&
189  n_coefficients_per_direction.size() == q_collection.size(),
190  ExcMessage("All parameters are supposed to have the same size."));
191 
192  const unsigned int max_n_coefficients_per_direction =
193  *std::max_element(n_coefficients_per_direction.cbegin(),
194  n_coefficients_per_direction.cend());
195  set_k_vectors(k_vectors, max_n_coefficients_per_direction);
196 
197  // reserve sufficient memory
199  }
200 
201 
202 
203  template <int dim, int spacedim>
205  const unsigned int n_coefficients_per_direction,
206  const hp::FECollection<dim, spacedim> &fe_collection,
207  const hp::QCollection<dim> & q_collection)
208  : Fourier<dim, spacedim>(
209  std::vector<unsigned int>(fe_collection.size(),
210  n_coefficients_per_direction),
211  fe_collection,
212  q_collection)
213  {}
214 
215 
216 
217  template <int dim, int spacedim>
218  inline bool
220  operator==(const Fourier<dim, spacedim> &fourier) const
221  {
222  return (
223  (n_coefficients_per_direction == fourier.n_coefficients_per_direction) &&
224  (*fe_collection == *(fourier.fe_collection)) &&
225  (q_collection == fourier.q_collection) &&
226  (k_vectors == fourier.k_vectors) &&
228  }
229 
230 
231 
232  template <int dim, int spacedim>
233  void
235  {
236  Threads::TaskGroup<> task_group;
237  for (unsigned int fe = 0; fe < fe_collection->size(); ++fe)
238  task_group += Threads::new_task([&, fe]() {
239  ensure_existence(n_coefficients_per_direction,
240  *fe_collection,
241  q_collection,
242  k_vectors,
243  fe,
245  });
246 
247  task_group.join_all();
248  }
249 
250 
251 
252  template <int dim, int spacedim>
253  unsigned int
255  const unsigned int index) const
256  {
257  return n_coefficients_per_direction[index];
258  }
259 
260 
261 
262  template <int dim, int spacedim>
263  template <typename Number>
264  void
266  const Vector<Number> & local_dof_values,
267  const unsigned int cell_active_fe_index,
268  Table<dim, CoefficientType> &fourier_coefficients)
269  {
270  for (unsigned int d = 0; d < dim; ++d)
271  AssertDimension(fourier_coefficients.size(d),
272  n_coefficients_per_direction[cell_active_fe_index]);
273 
274  ensure_existence(n_coefficients_per_direction,
275  *fe_collection,
276  q_collection,
277  k_vectors,
278  cell_active_fe_index,
280 
282  fourier_transform_matrices[cell_active_fe_index];
283 
284  unrolled_coefficients.resize(Utilities::fixed_power<dim>(
285  n_coefficients_per_direction[cell_active_fe_index]));
286  std::fill(unrolled_coefficients.begin(),
287  unrolled_coefficients.end(),
288  CoefficientType(0.));
289 
290  Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError());
291 
292  Assert(local_dof_values.size() == matrix.n(),
293  ExcDimensionMismatch(local_dof_values.size(), matrix.n()));
294 
295  for (unsigned int i = 0; i < unrolled_coefficients.size(); i++)
296  for (unsigned int j = 0; j < local_dof_values.size(); j++)
297  unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j];
298 
299  fourier_coefficients.fill(unrolled_coefficients.begin());
300  }
301 } // namespace FESeries
302 
303 
304 // explicit instantiations
305 #include "fe_series_fourier.inst"
306 
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:186
size_type m() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
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:1691
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:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
unsigned int size() const
Definition: q_collection.h:202
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
void precalculate_all_transformation_matrices()
typename std::complex< double > CoefficientType
Definition: fe_series.h:90
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:181
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:196
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:186
Definition: tensor.h:448
static constexpr double PI
Definition: numbers.h:231
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition: fe_series.h:201
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:206
static const char N
bool operator==(const Fourier< dim, spacedim > &fourier) const
Definition: table.h:687
const hp::QCollection< dim > q_collection
Definition: fe_series.h:191
void fill(InputIterator entries, const bool C_style_indexing=true)
double weight(const unsigned int i) const
static ::ExceptionBase & ExcInternalError()