Reference documentation for deal.II version GIT f9a1fd2718 2022-08-19 02:45:02+00:00
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
21 #include <deal.II/fe/fe_series.h>
22
23 #include <iostream>
24
25
27
28 namespace
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],
120  TableIndices<1>(k),
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
196 namespace 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,
217  ExcMessage(
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
237  const Legendre<dim, spacedim> &legendre) const
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  {
253  Threads::TaskGroup<> task_group;
254  for (unsigned int fe = 0; fe < fe_collection->size(); ++fe)
255  task_group += Threads::new_task([&, fe]() {
256  ensure_existence(n_coefficients_per_direction,
257  *fe_collection,
258  q_collection,
259  fe,
260  component,
261  legendre_transform_matrices);
262  });
263
264  task_group.join_all();
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:351
bool operator==(const Legendre< dim, spacedim > &legendre) const
const unsigned int component
Definition: fe_series.h:377
const hp::QCollection< dim > q_collection
Definition: fe_series.h:361
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:371
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:356
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
double weight(const unsigned int i) const
const Point< dim > & point(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
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#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 & 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:201
double legendre(unsigned int l, double x)
Definition: cmath.h:75
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)