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_legendre.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 DeclException2(ExcLegendre,
31 int,
32 double,
33 << "x[" << arg1 << "] = " << arg2 << " is not in [0,1]");
34
35 /*
36 * dim dimensional Legendre function with indices @p indices
37 * evaluated at @p x_q in [0,1]^dim.
38 */
39 template <int dim>
40 double
41 Lh(const Point<dim> &x_q, const TableIndices<dim> &indices)
42 {
43 double res = 1.0;
44 for (unsigned int d = 0; d < dim; d++)
45 {
46 const double x = 2.0 * (x_q[d] - 0.5);
47 Assert((x_q[d] <= 1.0) && (x_q[d] >= 0.), ExcLegendre(d, x_q[d]));
48 const unsigned int ind = indices[d];
49 res *= std::sqrt(2.0) * std_cxx17::legendre(ind, x);
50 }
51 return res;
52 }
53
54
55
56 /*
57 * Multiplier in Legendre coefficients
58 */
59 template <int dim>
60 double
61 multiplier(const TableIndices<dim> &indices)
62 {
63 double res = 1.0;
64 for (unsigned int d = 0; d < dim; d++)
65 res *= (0.5 + indices[d]);
66
67 return res;
68 }
69
70
71
72 template <int dim, int spacedim>
73 double
74 integrate(const FiniteElement<dim, spacedim> &fe,
75 const Quadrature<dim> & quadrature,
76 const TableIndices<dim> & indices,
77 const unsigned int dof,
78 const unsigned int component)
79 {
80 double sum = 0;
81 for (unsigned int q = 0; q < quadrature.size(); ++q)
82 {
83 const Point<dim> &x_q = quadrature.point(q);
84 sum += Lh(x_q, indices) *
85 fe.shape_value_component(dof, x_q, component) *
86 quadrature.weight(q);
87 }
88 return sum * multiplier(indices);
89 }
90
91
92
97 template <int spacedim>
98 void
99 ensure_existence(
100 const std::vector<unsigned int> & n_coefficients_per_direction,
101 const hp::FECollection<1, spacedim> &fe_collection,
102 const hp::QCollection<1> & q_collection,
103 const unsigned int fe,
104 const unsigned int component,
105 std::vector<FullMatrix<double>> & legendre_transform_matrices)
106 {
107 AssertIndexRange(fe, fe_collection.size());
108
109 if (legendre_transform_matrices[fe].m() == 0)
110 {
111 legendre_transform_matrices[fe].reinit(
112 n_coefficients_per_direction[fe],
113 fe_collection[fe].n_dofs_per_cell());
114
115 for (unsigned int k = 0; k < n_coefficients_per_direction[fe]; ++k)
116 for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell(); ++j)
117 legendre_transform_matrices[fe](k, j) =
118 integrate(fe_collection[fe],
119 q_collection[fe],
121 j,
122 component);
123 }
124 }
125
126 template <int spacedim>
127 void
128 ensure_existence(
129 const std::vector<unsigned int> & n_coefficients_per_direction,
130 const hp::FECollection<2, spacedim> &fe_collection,
131 const hp::QCollection<2> & q_collection,
132 const unsigned int fe,
133 const unsigned int component,
134 std::vector<FullMatrix<double>> & legendre_transform_matrices)
135 {
136 AssertIndexRange(fe, fe_collection.size());
137
138 if (legendre_transform_matrices[fe].m() == 0)
139 {
140 legendre_transform_matrices[fe].reinit(
141 Utilities::fixed_power<2>(n_coefficients_per_direction[fe]),
142 fe_collection[fe].n_dofs_per_cell());
143
144 unsigned int k = 0;
145 for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
146 for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe];
147 ++k2, k++)
148 for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
149 ++j)
150 legendre_transform_matrices[fe](k, j) =
151 integrate(fe_collection[fe],
152 q_collection[fe],
153 TableIndices<2>(k1, k2),
154 j,
155 component);
156 }
157 }
158
159 template <int spacedim>
160 void
161 ensure_existence(
162 const std::vector<unsigned int> & n_coefficients_per_direction,
163 const hp::FECollection<3, spacedim> &fe_collection,
164 const hp::QCollection<3> & q_collection,
165 const unsigned int fe,
166 const unsigned int component,
167 std::vector<FullMatrix<double>> & legendre_transform_matrices)
168 {
169 AssertIndexRange(fe, fe_collection.size());
170
171 if (legendre_transform_matrices[fe].m() == 0)
172 {
173 legendre_transform_matrices[fe].reinit(
174 Utilities::fixed_power<3>(n_coefficients_per_direction[fe]),
175 fe_collection[fe].n_dofs_per_cell());
176
177 unsigned int k = 0;
178 for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
179 for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe]; ++k2)
180 for (unsigned int k3 = 0; k3 < n_coefficients_per_direction[fe];
181 ++k3, k++)
182 for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell();
183 ++j)
184 legendre_transform_matrices[fe](k, j) =
185 integrate(fe_collection[fe],
186 q_collection[fe],
187 TableIndices<3>(k1, k2, k3),
188 j,
189 component);
190 }
191 }
192} // namespace
193
194
195
196namespace FESeries
197{
198 template <int dim, int spacedim>
200 const std::vector<unsigned int> & n_coefficients_per_direction,
201 const hp::FECollection<dim, spacedim> &fe_collection,
202 const hp::QCollection<dim> & q_collection,
203 const unsigned int component_)
204 : n_coefficients_per_direction(n_coefficients_per_direction)
205 , fe_collection(&fe_collection)
206 , q_collection(q_collection)
207 , legendre_transform_matrices(fe_collection.size())
208 , component(component_ != numbers::invalid_unsigned_int ? component_ : 0)
209 {
212 ExcMessage("All parameters are supposed to have the same size."));
213
214 if (fe_collection[0].n_components() > 1)
215 Assert(
216 component_ != numbers::invalid_unsigned_int,
218 "For vector-valued problems, you need to explicitly specify for "
219 "which vector component you will want to do a Fourier decomposition "
220 "by setting the 'component' argument of this constructor."));
221
222 AssertIndexRange(component, fe_collection[0].n_components());
223
224 // reserve sufficient memory
225 const unsigned int max_n_coefficients_per_direction =
226 *std::max_element(n_coefficients_per_direction.cbegin(),
228 unrolled_coefficients.reserve(
229 Utilities::fixed_power<dim>(max_n_coefficients_per_direction));
230 }
231
232
233
234 template <int dim, int spacedim>
236 const unsigned int n_coefficients_per_direction,
237 const hp::FECollection<dim, spacedim> &fe_collection,
238 const hp::QCollection<dim> & q_collection)
239 : Legendre<dim, spacedim>(
240 std::vector<unsigned int>(fe_collection.size(),
241 n_coefficients_per_direction),
242 fe_collection,
243 q_collection)
244 {}
245
246
247
248 template <int dim, int spacedim>
249 inline bool
252 {
253 return (
254 (n_coefficients_per_direction == legendre.n_coefficients_per_direction) &&
255 (*fe_collection == *(legendre.fe_collection)) &&
256 (q_collection == legendre.q_collection) &&
257 (legendre_transform_matrices == legendre.legendre_transform_matrices) &&
258 (component == legendre.component));
259 }
260
261
262
263 template <int dim, int spacedim>
264 void
266 {
267 Threads::TaskGroup<> task_group;
268 for (unsigned int fe = 0; fe < fe_collection->size(); ++fe)
269 task_group += Threads::new_task([&, fe]() {
270 ensure_existence(n_coefficients_per_direction,
271 *fe_collection,
272 q_collection,
273 fe,
274 component,
275 legendre_transform_matrices);
276 });
277
278 task_group.join_all();
279 }
280
281
282
283 template <int dim, int spacedim>
284 unsigned int
286 const unsigned int index) const
287 {
288 return n_coefficients_per_direction[index];
289 }
290
291
292
293 template <int dim, int spacedim>
294 template <typename Number>
295 void
297 const ::Vector<Number> &local_dof_values,
298 const unsigned int cell_active_fe_index,
299 Table<dim, CoefficientType> & legendre_coefficients)
300 {
301 for (unsigned int d = 0; d < dim; ++d)
302 AssertDimension(legendre_coefficients.size(d),
303 n_coefficients_per_direction[cell_active_fe_index]);
304
305 ensure_existence(n_coefficients_per_direction,
306 *fe_collection,
307 q_collection,
308 cell_active_fe_index,
309 component,
310 legendre_transform_matrices);
311
313 legendre_transform_matrices[cell_active_fe_index];
314
315 unrolled_coefficients.resize(Utilities::fixed_power<dim>(
316 n_coefficients_per_direction[cell_active_fe_index]));
317 std::fill(unrolled_coefficients.begin(),
318 unrolled_coefficients.end(),
319 CoefficientType(0.));
320
321 Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError());
322
323 Assert(local_dof_values.size() == matrix.n(),
324 ExcDimensionMismatch(local_dof_values.size(), matrix.n()));
325
326 for (unsigned int i = 0; i < unrolled_coefficients.size(); i++)
327 for (unsigned int j = 0; j < local_dof_values.size(); j++)
328 unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j];
329
330 legendre_coefficients.fill(unrolled_coefficients.begin());
331 }
332} // namespace FESeries
333
334
335// explicit instantiations
336#include "fe_series_legendre.inst"
337
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:376
double CoefficientType
Definition: fe_series.h:275
bool operator==(const Legendre< dim, spacedim > &legendre) const
const unsigned int component
Definition: fe_series.h:402
const hp::QCollection< dim > q_collection
Definition: fe_series.h:386
Legendre(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()
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:396
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:381
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &legendre_coefficients)
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)
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 DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
#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)
Task< RT > new_task(const std::function< RT()> &function)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
double legendre(unsigned int l, double x)
Definition: cmath.h:75
STL namespace.
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)