Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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 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
186namespace 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 {
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,
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(),
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) &&
249 (fourier_transform_matrices == fourier.fourier_transform_matrices) &&
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,
268 fourier_transform_matrices);
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,
304 fourier_transform_matrices);
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
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:208
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:218
const unsigned int component
Definition: fe_series.h:224
typename std::complex< double > CoefficientType
Definition: fe_series.h:91
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:193
bool operator==(const Fourier< dim, spacedim > &fourier) const
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition: fe_series.h:213
const hp::QCollection< dim > q_collection
Definition: fe_series.h:203
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:111
const Point< dim > & point(const unsigned int i) const
double weight(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:472
Definition: vector.h:110
unsigned int size() const
Definition: collection.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
size_type size() const
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:231
static const unsigned int invalid_unsigned_int
Definition: types.h:196
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)