Reference documentation for deal.II version 9.4.0
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
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,
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) *
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 Legendre 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>
235 inline bool
238 {
239 return (
240 (n_coefficients_per_direction == legendre.n_coefficients_per_direction) &&
241 (*fe_collection == *(legendre.fe_collection)) &&
242 (q_collection == legendre.q_collection) &&
243 (legendre_transform_matrices == legendre.legendre_transform_matrices) &&
244 (component == legendre.component));
245 }
246
247
248
249 template <int dim, int spacedim>
250 void
252 {
254 for (unsigned int fe = 0; fe < fe_collection->size(); ++fe)
256 ensure_existence(n_coefficients_per_direction,
257 *fe_collection,
258 q_collection,
259 fe,
260 component,
261 legendre_transform_matrices);
262 });
263
265 }
266
267
268
269 template <int dim, int spacedim>
270 unsigned int
272 const unsigned int index) const
273 {
274 return n_coefficients_per_direction[index];
275 }
276
277
278
279 template <int dim, int spacedim>
280 template <typename Number>
281 void
283 const ::Vector<Number> &local_dof_values,
284 const unsigned int cell_active_fe_index,
285 Table<dim, CoefficientType> & legendre_coefficients)
286 {
287 for (unsigned int d = 0; d < dim; ++d)
288 AssertDimension(legendre_coefficients.size(d),
289 n_coefficients_per_direction[cell_active_fe_index]);
290
291 ensure_existence(n_coefficients_per_direction,
292 *fe_collection,
293 q_collection,
294 cell_active_fe_index,
295 component,
296 legendre_transform_matrices);
297
299 legendre_transform_matrices[cell_active_fe_index];
300
301 unrolled_coefficients.resize(Utilities::fixed_power<dim>(
302 n_coefficients_per_direction[cell_active_fe_index]));
303 std::fill(unrolled_coefficients.begin(),
304 unrolled_coefficients.end(),
305 CoefficientType(0.));
306
307 Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError());
308
309 Assert(local_dof_values.size() == matrix.n(),
310 ExcDimensionMismatch(local_dof_values.size(), matrix.n()));
311
312 for (unsigned int i = 0; i < unrolled_coefficients.size(); ++i)
313 for (unsigned int j = 0; j < local_dof_values.size(); ++j)
314 unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j];
315
316 legendre_coefficients.fill(unrolled_coefficients.begin());
317 }
318} // namespace FESeries
319
320
321// explicit instantiations
322#include "fe_series_legendre.inst"
323
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:349
double CoefficientType
Definition: fe_series.h:261
bool operator==(const Legendre< dim, spacedim > &legendre) const
const unsigned int component
Definition: fe_series.h:375
const hp::QCollection< dim > q_collection
Definition: fe_series.h:359
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:369
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:354
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:264
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)