Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_series_legendre.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2018 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 
18 #include <deal.II/base/config.h>
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 
28 DEAL_II_NAMESPACE_OPEN
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(const hp::FECollection<1, spacedim> &fe_collection,
111  const hp::QCollection<1> & q_collection,
112  const unsigned int N,
113  const unsigned int fe,
114  std::vector<FullMatrix<double>> &legendre_transform_matrices)
115  {
116  AssertIndexRange(fe, fe_collection.size());
117 
118  if (legendre_transform_matrices[fe].m() == 0)
119  {
120  legendre_transform_matrices[fe].reinit(N,
121  fe_collection[fe].dofs_per_cell);
122 
123  for (unsigned int k = 0; k < N; ++k)
124  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
125  legendre_transform_matrices[fe](k, j) = integrate(
126  fe_collection[fe], q_collection[fe], TableIndices<1>(k), j);
127  }
128  }
129 
130  template <int spacedim>
131  void
132  ensure_existence(const hp::FECollection<2, spacedim> &fe_collection,
133  const hp::QCollection<2> & q_collection,
134  const unsigned int N,
135  const unsigned int fe,
136  std::vector<FullMatrix<double>> &legendre_transform_matrices)
137  {
138  AssertIndexRange(fe, fe_collection.size());
139 
140  if (legendre_transform_matrices[fe].m() == 0)
141  {
142  legendre_transform_matrices[fe].reinit(N * N,
143  fe_collection[fe].dofs_per_cell);
144 
145  unsigned int k = 0;
146  for (unsigned int k1 = 0; k1 < N; ++k1)
147  for (unsigned int k2 = 0; k2 < N; ++k2, k++)
148  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++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  }
155  }
156 
157  template <int spacedim>
158  void
159  ensure_existence(const hp::FECollection<3, spacedim> &fe_collection,
160  const hp::QCollection<3> & q_collection,
161  const unsigned int N,
162  const unsigned int fe,
163  std::vector<FullMatrix<double>> &legendre_transform_matrices)
164  {
165  AssertIndexRange(fe, fe_collection.size());
166 
167  if (legendre_transform_matrices[fe].m() == 0)
168  {
169  legendre_transform_matrices[fe].reinit(N * N * N,
170  fe_collection[fe].dofs_per_cell);
171 
172  unsigned int k = 0;
173  for (unsigned int k1 = 0; k1 < N; ++k1)
174  for (unsigned int k2 = 0; k2 < N; ++k2)
175  for (unsigned int k3 = 0; k3 < N; ++k3, k++)
176  for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j)
177  legendre_transform_matrices[fe](k, j) =
178  integrate(fe_collection[fe],
179  q_collection[fe],
180  TableIndices<3>(k1, k2, k3),
181  j);
182  }
183  }
184 } // namespace
185 
186 
187 
188 namespace FESeries
189 {
190  template <int dim, int spacedim>
192  const unsigned int size_in_each_direction,
193  const hp::FECollection<dim, spacedim> &fe_collection,
194  const hp::QCollection<dim> & q_collection)
195  : N(size_in_each_direction)
196  , fe_collection(&fe_collection)
197  , q_collection(&q_collection)
198  , legendre_transform_matrices(fe_collection.size())
199  , unrolled_coefficients(Utilities::fixed_power<dim>(N), 0.)
200  {}
201 
202 
203 
204  template <int dim, int spacedim>
205  template <typename Number>
206  void
208  const ::Vector<Number> &local_dof_values,
209  const unsigned int cell_active_fe_index,
210  Table<dim, CoefficientType> & legendre_coefficients)
211  {
212  ensure_existence(*fe_collection,
213  *q_collection,
214  N,
215  cell_active_fe_index,
216  legendre_transform_matrices);
217 
218  const FullMatrix<CoefficientType> &matrix =
219  legendre_transform_matrices[cell_active_fe_index];
220 
221  std::fill(unrolled_coefficients.begin(),
222  unrolled_coefficients.end(),
223  CoefficientType(0.));
224 
225  Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError());
226 
227  Assert(local_dof_values.size() == matrix.n(),
228  ExcDimensionMismatch(local_dof_values.size(), matrix.n()));
229 
230  for (unsigned int i = 0; i < unrolled_coefficients.size(); i++)
231  for (unsigned int j = 0; j < local_dof_values.size(); j++)
232  unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j];
233 
234  legendre_coefficients.fill(unrolled_coefficients.begin());
235  }
236 } // namespace FESeries
237 
238 
239 // explicit instantiations
240 #include "fe_series_legendre.inst"
241 
242 DEAL_II_NAMESPACE_CLOSE
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
Legendre(const unsigned int size_in_each_direction, const hp::FECollection< dim, spacedim > &fe_collection, const hp::QCollection< dim > &q_collection)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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
Definition: cuda.h:31
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:159
Definition: table.h:37
void fill(InputIterator entries, const bool C_style_indexing=true)
double weight(const unsigned int i) const
static ::ExceptionBase & ExcInternalError()