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