Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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
20
22
23#include <iostream>
24
25
27
28namespace
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
189namespace 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,
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
295 const FullMatrix<CoefficientType> &matrix =
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
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
size_type size() const
unsigned int size() const
Definition collection.h:265
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
Task< RT > new_task(const std::function< RT()> &function)
T sum(const T &t, const MPI_Comm mpi_communicator)
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 > &)