Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20: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\}}\)
Loading...
Searching...
No Matches
fe_series_fourier.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16
19
21
22#include <iostream>
23
24
26
27namespace
28{
29 void
30 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
38 set_k_vectors(Table<2, Tensor<1, 2>> &k_vectors, const unsigned int N)
39 {
40 k_vectors.reinit(TableIndices<2>(N, N));
41 for (unsigned int i = 0; i < N; ++i)
42 for (unsigned int j = 0; j < N; ++j)
43 {
44 k_vectors(i, j)[0] = 2. * numbers::PI * i;
45 k_vectors(i, j)[1] = 2. * numbers::PI * j;
46 }
47 }
48
49 void
50 set_k_vectors(Table<3, Tensor<1, 3>> &k_vectors, const unsigned int N)
51 {
52 k_vectors.reinit(TableIndices<3>(N, N, N));
53 for (unsigned int i = 0; i < N; ++i)
54 for (unsigned int j = 0; j < N; ++j)
55 for (unsigned int k = 0; k < N; ++k)
56 {
57 k_vectors(i, j, k)[0] = 2. * numbers::PI * i;
58 k_vectors(i, j, k)[1] = 2. * numbers::PI * j;
59 k_vectors(i, j, k)[2] = 2. * numbers::PI * k;
60 }
61 }
62
63
64
65 template <int dim, int spacedim>
66 std::complex<double>
67 integrate(const FiniteElement<dim, spacedim> &fe,
68 const Quadrature<dim> &quadrature,
69 const Tensor<1, dim> &k_vector,
70 const unsigned int j,
71 const unsigned int component)
72 {
73 std::complex<double> sum = 0;
74 for (unsigned int q = 0; q < quadrature.size(); ++q)
75 {
76 const Point<dim> &x_q = quadrature.point(q);
77 sum += std::exp(std::complex<double>(0, 1) * (k_vector * x_q)) *
78 fe.shape_value_component(j, x_q, component) *
79 quadrature.weight(q);
80 }
81 return sum;
82 }
83
84
85
86 /*
87 * Ensure that the transformation matrix for FiniteElement index
88 * @p fe_index is calculated. If not, calculate it.
89 */
90 template <int spacedim>
91 void
92 ensure_existence(
93 const std::vector<unsigned int> &n_coefficients_per_direction,
94 const hp::FECollection<1, spacedim> &fe_collection,
95 const hp::QCollection<1> &q_collection,
96 const Table<1, Tensor<1, 1>> &k_vectors,
97 const unsigned int fe,
98 const unsigned int component,
99 std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
100 {
101 AssertIndexRange(fe, fe_collection.size());
102
103 if (fourier_transform_matrices[fe].m() == 0)
104 {
105 fourier_transform_matrices[fe].reinit(
106 n_coefficients_per_direction[fe],
107 fe_collection[fe].n_dofs_per_cell());
108
109 for (unsigned int k = 0; k < n_coefficients_per_direction[fe]; ++k)
110 for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell(); ++j)
111 fourier_transform_matrices[fe](k, j) = integrate(
112 fe_collection[fe], q_collection[fe], k_vectors(k), j, component);
113 }
114 }
115
116 template <int spacedim>
117 void
118 ensure_existence(
119 const std::vector<unsigned int> &n_coefficients_per_direction,
120 const hp::FECollection<2, spacedim> &fe_collection,
121 const hp::QCollection<2> &q_collection,
122 const Table<2, Tensor<1, 2>> &k_vectors,
123 const unsigned int fe,
124 const unsigned int component,
125 std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
126 {
127 AssertIndexRange(fe, fe_collection.size());
128
129 if (fourier_transform_matrices[fe].m() == 0)
130 {
131 fourier_transform_matrices[fe].reinit(
132 Utilities::fixed_power<2>(n_coefficients_per_direction[fe]),
133 fe_collection[fe].n_dofs_per_cell());
134
135 unsigned int k = 0;
136 for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
137 for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe];
138 ++k2, ++k)
139 for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
140 ++j)
141 fourier_transform_matrices[fe](k, j) =
142 integrate(fe_collection[fe],
143 q_collection[fe],
144 k_vectors(k1, k2),
145 j,
146 component);
147 }
148 }
149
150 template <int spacedim>
151 void
152 ensure_existence(
153 const std::vector<unsigned int> &n_coefficients_per_direction,
154 const hp::FECollection<3, spacedim> &fe_collection,
155 const hp::QCollection<3> &q_collection,
156 const Table<3, Tensor<1, 3>> &k_vectors,
157 const unsigned int fe,
158 const unsigned int component,
159 std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
160 {
161 AssertIndexRange(fe, fe_collection.size());
162
163 if (fourier_transform_matrices[fe].m() == 0)
164 {
165 fourier_transform_matrices[fe].reinit(
166 Utilities::fixed_power<3>(n_coefficients_per_direction[fe]),
167 fe_collection[fe].n_dofs_per_cell());
168
169 unsigned int k = 0;
170 for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
171 for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe]; ++k2)
172 for (unsigned int k3 = 0; k3 < n_coefficients_per_direction[fe];
173 ++k3, ++k)
174 for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
175 ++j)
176 fourier_transform_matrices[fe](k, j) =
177 integrate(fe_collection[fe],
178 q_collection[fe],
179 k_vectors(k1, k2, k3),
180 j,
181 component);
182 }
183 }
184} // namespace
185
186
187
188namespace FESeries
189{
190 template <int dim, int spacedim>
192 const std::vector<unsigned int> &n_coefficients_per_direction,
193 const hp::FECollection<dim, spacedim> &fe_collection,
194 const hp::QCollection<dim> &q_collection,
195 const unsigned int component_)
196 : n_coefficients_per_direction(n_coefficients_per_direction)
197 , fe_collection(&fe_collection)
198 , q_collection(q_collection)
199 , fourier_transform_matrices(fe_collection.size())
200 , component(component_ != numbers::invalid_unsigned_int ? component_ : 0)
201 {
204 ExcMessage("All parameters are supposed to have the same size."));
205
206 if (fe_collection[0].n_components() > 1)
207 Assert(
208 component_ != numbers::invalid_unsigned_int,
210 "For vector-valued problems, you need to explicitly specify for "
211 "which vector component you will want to do a Fourier decomposition "
212 "by setting the 'component' argument of this constructor."));
213
214 AssertIndexRange(component, fe_collection[0].n_components());
215
216 const unsigned int max_n_coefficients_per_direction =
217 *std::max_element(n_coefficients_per_direction.cbegin(),
219 set_k_vectors(k_vectors, max_n_coefficients_per_direction);
220
221 // reserve sufficient memory
222 unrolled_coefficients.reserve(k_vectors.n_elements());
223 }
224
225
226
227 template <int dim, int spacedim>
228 inline bool
230 const Fourier<dim, spacedim> &fourier) const
231 {
232 return (
233 (n_coefficients_per_direction == fourier.n_coefficients_per_direction) &&
234 (*fe_collection == *(fourier.fe_collection)) &&
235 (q_collection == fourier.q_collection) &&
236 (k_vectors == fourier.k_vectors) &&
237 (fourier_transform_matrices == fourier.fourier_transform_matrices) &&
238 (component == fourier.component));
239 }
240
241
242
243 template <int dim, int spacedim>
244 void
246 {
247 Threads::TaskGroup<> task_group;
248 for (unsigned int fe = 0; fe < fe_collection->size(); ++fe)
249 task_group += Threads::new_task([&, fe]() {
250 ensure_existence(n_coefficients_per_direction,
251 *fe_collection,
252 q_collection,
253 k_vectors,
254 fe,
255 component,
256 fourier_transform_matrices);
257 });
258
259 task_group.join_all();
260 }
261
262
263
264 template <int dim, int spacedim>
265 unsigned int
267 const unsigned int index) const
268 {
269 return n_coefficients_per_direction[index];
270 }
271
272
273
274 template <int dim, int spacedim>
275 template <typename Number>
276 void
278 const Vector<Number> &local_dof_values,
279 const unsigned int cell_active_fe_index,
280 Table<dim, CoefficientType> &fourier_coefficients)
281 {
282 for (unsigned int d = 0; d < dim; ++d)
283 AssertDimension(fourier_coefficients.size(d),
284 n_coefficients_per_direction[cell_active_fe_index]);
285
286 ensure_existence(n_coefficients_per_direction,
287 *fe_collection,
288 q_collection,
289 k_vectors,
290 cell_active_fe_index,
291 component,
292 fourier_transform_matrices);
293
294 const FullMatrix<CoefficientType> &matrix =
295 fourier_transform_matrices[cell_active_fe_index];
296
297 unrolled_coefficients.resize(Utilities::fixed_power<dim>(
298 n_coefficients_per_direction[cell_active_fe_index]));
299 std::fill(unrolled_coefficients.begin(),
300 unrolled_coefficients.end(),
301 CoefficientType(0.));
302
303 Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError());
304
305 Assert(local_dof_values.size() == matrix.n(),
306 ExcDimensionMismatch(local_dof_values.size(), matrix.n()));
307
308 for (unsigned int i = 0; i < unrolled_coefficients.size(); ++i)
309 for (unsigned int j = 0; j < local_dof_values.size(); ++j)
310 unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j];
311
312 fourier_coefficients.fill(unrolled_coefficients.begin());
313 }
314} // namespace FESeries
315
316
317// explicit instantiations
318#include "fe_series_fourier.inst"
319
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition fe_series.h:185
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:195
std::vector< CoefficientType > unrolled_coefficients
Definition fe_series.h:205
const unsigned int component
Definition fe_series.h:211
typename std::complex< double > CoefficientType
Definition fe_series.h:92
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:180
bool operator==(const Fourier< dim, spacedim > &fourier) const
std::vector< FullMatrix< CoefficientType > > fourier_transform_matrices
Definition fe_series.h:200
const hp::QCollection< dim > q_collection
Definition fe_series.h:190
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
virtual size_type size() const override
unsigned int size() const
Definition collection.h:264
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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:220
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)