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