Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
tensor_product_polynomials_bubbles.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2012 - 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 #ifndef dealii_tensor_product_polynomials_bubbles_h
17 #define dealii_tensor_product_polynomials_bubbles_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/point.h>
24 #include <deal.II/base/polynomial.h>
25 #include <deal.II/base/tensor.h>
26 #include <deal.II/base/tensor_product_polynomials.h>
27 #include <deal.II/base/utilities.h>
28 
29 #include <vector>
30 
31 DEAL_II_NAMESPACE_OPEN
32 
33 
47 template <int dim>
49 {
50 public:
56  template <class Pol>
57  TensorProductPolynomialsBubbles(const std::vector<Pol> &pols);
58 
71  void
72  compute(const Point<dim> & unit_point,
73  std::vector<double> & values,
74  std::vector<Tensor<1, dim>> &grads,
75  std::vector<Tensor<2, dim>> &grad_grads,
76  std::vector<Tensor<3, dim>> &third_derivatives,
77  std::vector<Tensor<4, dim>> &fourth_derivatives) const;
78 
91  double
92  compute_value(const unsigned int i, const Point<dim> &p) const;
93 
106  template <int order>
108  compute_derivative(const unsigned int i, const Point<dim> &p) const;
109 
123  compute_grad(const unsigned int i, const Point<dim> &p) const;
124 
138  compute_grad_grad(const unsigned int i, const Point<dim> &p) const;
139 
146  unsigned int
147  n() const;
148 };
149 
153 /* ---------------- template and inline functions ---------- */
154 
155 #ifndef DOXYGEN
156 
157 template <int dim>
158 template <class Pol>
160  const std::vector<Pol> &pols)
161  : TensorProductPolynomials<dim>(pols)
162 {
163  const unsigned int q_degree = this->polynomials.size() - 1;
164  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
165  // append index for renumbering
166  for (unsigned int i = 0; i < n_bubbles; ++i)
167  {
168  this->index_map.push_back(i + this->n_tensor_pols);
169  this->index_map_inverse.push_back(i + this->n_tensor_pols);
170  }
171 }
172 
173 
174 
175 template <int dim>
176 inline unsigned int
178 {
179  return this->n_tensor_pols + dim;
180 }
181 
182 
183 
184 template <>
185 inline unsigned int
187 {
189 }
190 
191 template <int dim>
192 template <int order>
195  const unsigned int i,
196  const Point<dim> & p) const
197 {
198  const unsigned int q_degree = this->polynomials.size() - 1;
199  const unsigned int max_q_indices = this->n_tensor_pols;
200  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
201  (void)n_bubbles;
202  Assert(i < max_q_indices + n_bubbles, ExcInternalError());
203 
204  // treat the regular basis functions
205  if (i < max_q_indices)
206  return this
207  ->TensorProductPolynomials<dim>::template compute_derivative<order>(i, p);
208 
209  const unsigned int comp = i - this->n_tensor_pols;
210 
211  Tensor<order, dim> derivative;
212  switch (order)
213  {
214  case 1:
215  {
216  Tensor<1, dim> &derivative_1 =
217  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
218 
219  for (unsigned int d = 0; d < dim; ++d)
220  {
221  derivative_1[d] = 1.;
222  // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
223  for (unsigned j = 0; j < dim; ++j)
224  derivative_1[d] *=
225  (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
226  // and multiply with (2*x_i-1)^{r-1}
227  for (unsigned int i = 0; i < q_degree - 1; ++i)
228  derivative_1[d] *= 2 * p(comp) - 1;
229  }
230 
231  if (q_degree >= 2)
232  {
233  // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
234  double value = 1.;
235  for (unsigned int j = 0; j < dim; ++j)
236  value *= 4 * p(j) * (1 - p(j));
237  // and multiply with grad(2*x_i-1)^{r-1}
238  double tmp = value * 2 * (q_degree - 1);
239  for (unsigned int i = 0; i < q_degree - 2; ++i)
240  tmp *= 2 * p(comp) - 1;
241  derivative_1[comp] += tmp;
242  }
243 
244  return derivative;
245  }
246  case 2:
247  {
248  Tensor<2, dim> &derivative_2 =
249  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
250 
251  double v[dim + 1][3];
252  {
253  for (unsigned int c = 0; c < dim; ++c)
254  {
255  v[c][0] = 4 * p(c) * (1 - p(c));
256  v[c][1] = 4 * (1 - 2 * p(c));
257  v[c][2] = -8;
258  }
259 
260  double tmp = 1.;
261  for (unsigned int i = 0; i < q_degree - 1; ++i)
262  tmp *= 2 * p(comp) - 1;
263  v[dim][0] = tmp;
264 
265  if (q_degree >= 2)
266  {
267  double tmp = 2 * (q_degree - 1);
268  for (unsigned int i = 0; i < q_degree - 2; ++i)
269  tmp *= 2 * p(comp) - 1;
270  v[dim][1] = tmp;
271  }
272  else
273  v[dim][1] = 0.;
274 
275  if (q_degree >= 3)
276  {
277  double tmp = 4 * (q_degree - 2) * (q_degree - 1);
278  for (unsigned int i = 0; i < q_degree - 3; ++i)
279  tmp *= 2 * p(comp) - 1;
280  v[dim][2] = tmp;
281  }
282  else
283  v[dim][2] = 0.;
284  }
285 
286  // calculate (\partial_j \partial_k \psi) * monomial
287  Tensor<2, dim> grad_grad_1;
288  for (unsigned int d1 = 0; d1 < dim; ++d1)
289  for (unsigned int d2 = 0; d2 < dim; ++d2)
290  {
291  grad_grad_1[d1][d2] = v[dim][0];
292  for (unsigned int x = 0; x < dim; ++x)
293  {
294  unsigned int derivative = 0;
295  if (d1 == x || d2 == x)
296  {
297  if (d1 == d2)
298  derivative = 2;
299  else
300  derivative = 1;
301  }
302  grad_grad_1[d1][d2] *= v[x][derivative];
303  }
304  }
305 
306  // calculate (\partial_j \psi) *(\partial_k monomial)
307  // and (\partial_k \psi) *(\partial_j monomial)
308  Tensor<2, dim> grad_grad_2;
309  Tensor<2, dim> grad_grad_3;
310  for (unsigned int d = 0; d < dim; ++d)
311  {
312  grad_grad_2[d][comp] = v[dim][1];
313  grad_grad_3[comp][d] = v[dim][1];
314  for (unsigned int x = 0; x < dim; ++x)
315  {
316  grad_grad_2[d][comp] *= v[x][d == x];
317  grad_grad_3[comp][d] *= v[x][d == x];
318  }
319  }
320 
321  // calculate \psi *(\partial j \partial_k monomial) and sum
322  double psi_value = 1.;
323  for (unsigned int x = 0; x < dim; ++x)
324  psi_value *= v[x][0];
325 
326  for (unsigned int d1 = 0; d1 < dim; ++d1)
327  for (unsigned int d2 = 0; d2 < dim; ++d2)
328  derivative_2[d1][d2] =
329  grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
330  derivative_2[comp][comp] += psi_value * v[dim][2];
331 
332  return derivative;
333  }
334  default:
335  {
336  Assert(false, ExcNotImplemented());
337  return derivative;
338  }
339  }
340 }
341 
342 
343 #endif // DOXYGEN
344 DEAL_II_NAMESPACE_CLOSE
345 
346 #endif
static const unsigned int invalid_unsigned_int
Definition: types.h:173
std::vector< Polynomials::Polynomial< double > > polynomials
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
void compute(const Point< dim > &unit_point, std::vector< double > &values, std::vector< Tensor< 1, dim >> &grads, std::vector< Tensor< 2, dim >> &grad_grads, std::vector< Tensor< 3, dim >> &third_derivatives, std::vector< Tensor< 4, dim >> &fourth_derivatives) const
#define Assert(cond, exc)
Definition: exceptions.h:1407
std::vector< unsigned int > index_map_inverse
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TensorProductPolynomialsBubbles(const std::vector< Pol > &pols)
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
double compute_value(const unsigned int i, const Point< dim > &p) const
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInternalError()