Reference documentation for deal.II version Git 6cf3cbfd16 2020-07-12 08:18:15 -0400
\(\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\}}\)
fe_series_legendre.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2020 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 
19 
20 #include <deal.II/fe/fe_series.h>
21 
22 #include <iostream>
23 #ifdef DEAL_II_WITH_GSL
24 # include <gsl/gsl_sf_legendre.h>
25 #endif
26 
27 
29 
30 namespace
31 {
32  DeclException2(ExcLegendre,
33  int,
34  double,
35  << "x[" << arg1 << "] = " << arg2 << " is not in [0,1]");
36 
37  /*
38  * dim dimensional Legendre function with indices @p indices
39  * evaluated at @p x_q in [0,1]^dim.
40  */
41  template <int dim>
42  double
43  Lh(const Point<dim> &x_q, const TableIndices<dim> &indices)
44  {
45 #ifdef DEAL_II_WITH_GSL
46  double res = 1.0;
47  for (unsigned int d = 0; d < dim; d++)
48  {
49  const double x = 2.0 * (x_q[d] - 0.5);
50  Assert((x_q[d] <= 1.0) && (x_q[d] >= 0.), ExcLegendre(d, x_q[d]));
51  const int ind = indices[d];
52  res *= std::sqrt(2.0) * gsl_sf_legendre_Pl(ind, x);
53  }
54  return res;
55 
56 #else
57 
58  (void)x_q;
59  (void)indices;
60  AssertThrow(false,
61  ExcMessage("deal.II has to be configured with GSL "
62  "in order to use Legendre transformation."));
63  return 0;
64 #endif
65  }
66 
67 
68 
69  /*
70  * Multiplier in Legendre coefficients
71  */
72  template <int dim>
73  double
74  multiplier(const TableIndices<dim> &indices)
75  {
76  double res = 1.0;
77  for (unsigned int d = 0; d < dim; d++)
78  res *= (0.5 + indices[d]);
79 
80  return res;
81  }
82 
83 
84 
85  template <int dim, int spacedim>
86  double
87  integrate(const FiniteElement<dim, spacedim> &fe,
88  const Quadrature<dim> & quadrature,
89  const TableIndices<dim> & indices,
90  const unsigned int dof)
91  {
92  double sum = 0;
93  for (unsigned int q = 0; q < quadrature.size(); ++q)
94  {
95  const Point<dim> &x_q = quadrature.point(q);
96  sum +=
97  Lh(x_q, indices) * fe.shape_value(dof, x_q) * quadrature.weight(q);
98  }
99  return sum * multiplier(indices);
100  }
101 
102 
103 
108  template <int spacedim>
109  void
110  ensure_existence(
111  const std::vector<unsigned int> & n_coefficients_per_direction,
112  const hp::FECollection<1, spacedim> &fe_collection,
113  const hp::QCollection<1> & q_collection,
114  const unsigned int fe,
115  std::vector<FullMatrix<double>> & legendre_transform_matrices)
116  {
117  AssertIndexRange(fe, fe_collection.size());
118 
119  if (legendre_transform_matrices[fe].m() == 0)
120  {
121  legendre_transform_matrices[fe].reinit(n_coefficients_per_direction[fe],
122  fe_collection[fe].dofs_per_cell);
123 
124  for (unsigned int k = 0; k < n_coefficients_per_direction[fe]; ++k)
125  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
126  legendre_transform_matrices[fe](k, j) = integrate(
127  fe_collection[fe], q_collection[fe], TableIndices<1>(k), j);
128  }
129  }
130 
131  template <int spacedim>
132  void
133  ensure_existence(
134  const std::vector<unsigned int> & n_coefficients_per_direction,
135  const hp::FECollection<2, spacedim> &fe_collection,
136  const hp::QCollection<2> & q_collection,
137  const unsigned int fe,
138  std::vector<FullMatrix<double>> & legendre_transform_matrices)
139  {
140  AssertIndexRange(fe, fe_collection.size());
141 
142  if (legendre_transform_matrices[fe].m() == 0)
143  {
144  legendre_transform_matrices[fe].reinit(
145  Utilities::fixed_power<2>(n_coefficients_per_direction[fe]),
146  fe_collection[fe].dofs_per_cell);
147 
148  unsigned int k = 0;
149  for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
150  for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe];
151  ++k2, k++)
152  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
153  legendre_transform_matrices[fe](k, j) =
154  integrate(fe_collection[fe],
155  q_collection[fe],
156  TableIndices<2>(k1, k2),
157  j);
158  }
159  }
160 
161  template <int spacedim>
162  void
163  ensure_existence(
164  const std::vector<unsigned int> & n_coefficients_per_direction,
165  const hp::FECollection<3, spacedim> &fe_collection,
166  const hp::QCollection<3> & q_collection,
167  const unsigned int fe,
168  std::vector<FullMatrix<double>> & legendre_transform_matrices)
169  {
170  AssertIndexRange(fe, fe_collection.size());
171 
172  if (legendre_transform_matrices[fe].m() == 0)
173  {
174  legendre_transform_matrices[fe].reinit(
175  Utilities::fixed_power<3>(n_coefficients_per_direction[fe]),
176  fe_collection[fe].dofs_per_cell);
177 
178  unsigned int k = 0;
179  for (unsigned int k1 = 0; k1 < n_coefficients_per_direction[fe]; ++k1)
180  for (unsigned int k2 = 0; k2 < n_coefficients_per_direction[fe]; ++k2)
181  for (unsigned int k3 = 0; k3 < n_coefficients_per_direction[fe];
182  ++k3, k++)
183  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++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  }
190  }
191 } // namespace
192 
193 
194 
195 namespace 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  : n_coefficients_per_direction(n_coefficients_per_direction)
203  , fe_collection(&fe_collection)
204  , q_collection(q_collection)
205  , legendre_transform_matrices(fe_collection.size())
206  {
207  Assert(n_coefficients_per_direction.size() == fe_collection.size() &&
208  n_coefficients_per_direction.size() == q_collection.size(),
209  ExcMessage("All parameters are supposed to have the same size."));
210 
211  // reserve sufficient memory
212  const unsigned int max_n_coefficients_per_direction =
213  *std::max_element(n_coefficients_per_direction.cbegin(),
214  n_coefficients_per_direction.cend());
215  unrolled_coefficients.reserve(
216  Utilities::fixed_power<dim>(max_n_coefficients_per_direction));
217  }
218 
219 
220 
221  template <int dim, int spacedim>
223  const unsigned int n_coefficients_per_direction,
224  const hp::FECollection<dim, spacedim> &fe_collection,
225  const hp::QCollection<dim> & q_collection)
226  : Legendre<dim, spacedim>(
227  std::vector<unsigned int>(fe_collection.size(),
228  n_coefficients_per_direction),
229  fe_collection,
230  q_collection)
231  {}
232 
233 
234 
235  template <int dim, int spacedim>
236  inline bool
238  operator==(const Legendre<dim, spacedim> &legendre) const
239  {
240  return (
241  (n_coefficients_per_direction == legendre.n_coefficients_per_direction) &&
242  (*fe_collection == *(legendre.fe_collection)) &&
243  (q_collection == legendre.q_collection) &&
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,
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,
295 
297  legendre_transform_matrices[cell_active_fe_index];
298 
299  unrolled_coefficients.resize(Utilities::fixed_power<dim>(
300  n_coefficients_per_direction[cell_active_fe_index]));
301  std::fill(unrolled_coefficients.begin(),
302  unrolled_coefficients.end(),
303  CoefficientType(0.));
304 
305  Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError());
306 
307  Assert(local_dof_values.size() == matrix.n(),
308  ExcDimensionMismatch(local_dof_values.size(), matrix.n()));
309 
310  for (unsigned int i = 0; i < unrolled_coefficients.size(); i++)
311  for (unsigned int j = 0; j < local_dof_values.size(); j++)
312  unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j];
313 
314  legendre_coefficients.fill(unrolled_coefficients.begin());
315  }
316 } // namespace FESeries
317 
318 
319 // explicit instantiations
320 #include "fe_series_legendre.inst"
321 
size_type m() const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1560
Contents is actually a matrix.
Task< RT > new_task(const std::function< RT()> &function)
const hp::QCollection< dim > q_collection
Definition: fe_series.h:356
std::vector< CoefficientType > unrolled_coefficients
Definition: fe_series.h:366
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
unsigned int get_n_coefficients_per_direction(const unsigned int index) const
const Point< dim > & point(const unsigned int i) const
STL namespace.
unsigned int size() const
Legendre(const std::vector< unsigned int > &n_coefficients_per_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
std::vector< FullMatrix< CoefficientType > > legendre_transform_matrices
Definition: fe_series.h:361
const std::vector< unsigned int > n_coefficients_per_direction
Definition: fe_series.h:346
static ::ExceptionBase & ExcMessage(std::string arg1)
bool operator==(const Legendre< dim, spacedim > &legendre) const
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1403
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
unsigned int size() const
Definition: q_collection.h:174
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
SmartPointer< const hp::FECollection< dim, spacedim > > fe_collection
Definition: fe_series.h:351
void calculate(const ::Vector< Number > &local_dof_values, const unsigned int cell_active_fe_index, Table< dim, CoefficientType > &legendre_coefficients)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
unsigned int size() const
void precalculate_all_transformation_matrices()
size_type size(const unsigned int i) const
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:158
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
double CoefficientType
Definition: fe_series.h:256
Definition: table.h:687
void fill(InputIterator entries, const bool C_style_indexing=true)
double weight(const unsigned int i) const
inline ::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
static ::ExceptionBase & ExcInternalError()