Reference documentation for deal.II version GIT 3c37cfefb2 2023-03-24 04:00:03+00:00
\(\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\}}\)
tensor_product_polynomials_bubbles.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2012 - 2022 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 
23 #include <deal.II/base/point.h>
25 #include <deal.II/base/tensor.h>
27 #include <deal.II/base/utilities.h>
28 
29 #include <vector>
30 
32 
52 template <int dim>
54 {
55 public:
60  static constexpr unsigned int dimension = dim;
61 
67  template <class Pol>
68  TensorProductPolynomialsBubbles(const std::vector<Pol> &pols);
69 
73  void
74  output_indices(std::ostream &out) const;
75 
81  void
82  set_numbering(const std::vector<unsigned int> &renumber);
83 
87  const std::vector<unsigned int> &
88  get_numbering() const;
89 
93  const std::vector<unsigned int> &
95 
108  void
109  evaluate(const Point<dim> & unit_point,
110  std::vector<double> & values,
111  std::vector<Tensor<1, dim>> &grads,
112  std::vector<Tensor<2, dim>> &grad_grads,
113  std::vector<Tensor<3, dim>> &third_derivatives,
114  std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
115 
128  double
129  compute_value(const unsigned int i, const Point<dim> &p) const override;
130 
143  template <int order>
145  compute_derivative(const unsigned int i, const Point<dim> &p) const;
146 
150  virtual Tensor<1, dim>
151  compute_1st_derivative(const unsigned int i,
152  const Point<dim> & p) const override;
153 
157  virtual Tensor<2, dim>
158  compute_2nd_derivative(const unsigned int i,
159  const Point<dim> & p) const override;
160 
164  virtual Tensor<3, dim>
165  compute_3rd_derivative(const unsigned int i,
166  const Point<dim> & p) const override;
167 
171  virtual Tensor<4, dim>
172  compute_4th_derivative(const unsigned int i,
173  const Point<dim> & p) const override;
174 
188  compute_grad(const unsigned int i, const Point<dim> &p) const override;
189 
203  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
204 
211  unsigned int
212  n() const;
213 
218  std::string
219  name() const override;
220 
224  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
225  clone() const override;
226 
227 private:
232 
236  std::vector<unsigned int> index_map;
237 
241  std::vector<unsigned int> index_map_inverse;
242 };
243 
247 /* ---------------- template and inline functions ---------- */
248 
249 #ifndef DOXYGEN
250 
251 template <int dim>
252 template <class Pol>
254  const std::vector<Pol> &pols)
255  : ScalarPolynomialsBase<dim>(1,
256  Utilities::fixed_power<dim>(pols.size()) + dim)
257  , tensor_polys(pols)
258  , index_map(tensor_polys.n() +
259  ((tensor_polys.polynomials.size() <= 2) ? 1 : dim))
260  , index_map_inverse(tensor_polys.n() +
261  ((tensor_polys.polynomials.size() <= 2) ? 1 : dim))
262 {
263  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
264  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
265  // append index for renumbering
266  for (unsigned int i = 0; i < tensor_polys.n() + n_bubbles; ++i)
267  {
268  index_map[i] = i;
269  index_map_inverse[i] = i;
270  }
271 }
272 
273 
274 template <int dim>
275 inline unsigned int
277 {
278  return tensor_polys.n() + dim;
279 }
280 
281 
282 template <>
283 inline unsigned int
285 {
287 }
288 
289 
290 template <int dim>
291 inline const std::vector<unsigned int> &
293 {
294  return index_map;
295 }
296 
297 
298 template <int dim>
299 inline const std::vector<unsigned int> &
301 {
302  return index_map_inverse;
303 }
304 
305 
306 template <int dim>
307 inline std::string
309 {
310  return "TensorProductPolynomialsBubbles";
311 }
312 
313 
314 template <int dim>
315 template <int order>
318  const unsigned int i,
319  const Point<dim> & p) const
320 {
321  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
322  const unsigned int max_q_indices = tensor_polys.n();
323  Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
324  ExcInternalError());
325 
326  // treat the regular basis functions
327  if (i < max_q_indices)
328  return tensor_polys.template compute_derivative<order>(i, p);
329 
330  const unsigned int comp = i - tensor_polys.n();
331 
332  Tensor<order, dim> derivative;
333  switch (order)
334  {
335  case 1:
336  {
337  Tensor<1, dim> &derivative_1 =
338  *reinterpret_cast<Tensor<1, dim> *>(&derivative);
339 
340  for (unsigned int d = 0; d < dim; ++d)
341  {
342  derivative_1[d] = 1.;
343  // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
344  for (unsigned j = 0; j < dim; ++j)
345  derivative_1[d] *=
346  (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
347  // and multiply with (2*x_i-1)^{r-1}
348  for (unsigned int i = 0; i < q_degree - 1; ++i)
349  derivative_1[d] *= 2 * p(comp) - 1;
350  }
351 
352  if (q_degree >= 2)
353  {
354  // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
355  double value = 1.;
356  for (unsigned int j = 0; j < dim; ++j)
357  value *= 4 * p(j) * (1 - p(j));
358  // and multiply with grad(2*x_i-1)^{r-1}
359  double tmp = value * 2 * (q_degree - 1);
360  for (unsigned int i = 0; i < q_degree - 2; ++i)
361  tmp *= 2 * p(comp) - 1;
362  derivative_1[comp] += tmp;
363  }
364 
365  return derivative;
366  }
367  case 2:
368  {
369  Tensor<2, dim> &derivative_2 =
370  *reinterpret_cast<Tensor<2, dim> *>(&derivative);
371 
372  double v[dim + 1][3];
373  {
374  for (unsigned int c = 0; c < dim; ++c)
375  {
376  v[c][0] = 4 * p(c) * (1 - p(c));
377  v[c][1] = 4 * (1 - 2 * p(c));
378  v[c][2] = -8;
379  }
380 
381  double tmp = 1.;
382  for (unsigned int i = 0; i < q_degree - 1; ++i)
383  tmp *= 2 * p(comp) - 1;
384  v[dim][0] = tmp;
385 
386  if (q_degree >= 2)
387  {
388  double tmp = 2 * (q_degree - 1);
389  for (unsigned int i = 0; i < q_degree - 2; ++i)
390  tmp *= 2 * p(comp) - 1;
391  v[dim][1] = tmp;
392  }
393  else
394  v[dim][1] = 0.;
395 
396  if (q_degree >= 3)
397  {
398  double tmp = 4 * (q_degree - 2) * (q_degree - 1);
399  for (unsigned int i = 0; i < q_degree - 3; ++i)
400  tmp *= 2 * p(comp) - 1;
401  v[dim][2] = tmp;
402  }
403  else
404  v[dim][2] = 0.;
405  }
406 
407  // calculate (\partial_j \partial_k \psi) * monomial
408  Tensor<2, dim> grad_grad_1;
409  for (unsigned int d1 = 0; d1 < dim; ++d1)
410  for (unsigned int d2 = 0; d2 < dim; ++d2)
411  {
412  grad_grad_1[d1][d2] = v[dim][0];
413  for (unsigned int x = 0; x < dim; ++x)
414  {
415  unsigned int derivative = 0;
416  if (d1 == x || d2 == x)
417  {
418  if (d1 == d2)
419  derivative = 2;
420  else
421  derivative = 1;
422  }
423  grad_grad_1[d1][d2] *= v[x][derivative];
424  }
425  }
426 
427  // calculate (\partial_j \psi) *(\partial_k monomial)
428  // and (\partial_k \psi) *(\partial_j monomial)
429  Tensor<2, dim> grad_grad_2;
430  Tensor<2, dim> grad_grad_3;
431  for (unsigned int d = 0; d < dim; ++d)
432  {
433  grad_grad_2[d][comp] = v[dim][1];
434  grad_grad_3[comp][d] = v[dim][1];
435  for (unsigned int x = 0; x < dim; ++x)
436  {
437  grad_grad_2[d][comp] *= v[x][d == x];
438  grad_grad_3[comp][d] *= v[x][d == x];
439  }
440  }
441 
442  // calculate \psi *(\partial j \partial_k monomial) and sum
443  double psi_value = 1.;
444  for (unsigned int x = 0; x < dim; ++x)
445  psi_value *= v[x][0];
446 
447  for (unsigned int d1 = 0; d1 < dim; ++d1)
448  for (unsigned int d2 = 0; d2 < dim; ++d2)
449  derivative_2[d1][d2] =
450  grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
451  derivative_2[comp][comp] += psi_value * v[dim][2];
452 
453  return derivative;
454  }
455  default:
456  {
457  Assert(false, ExcNotImplemented());
458  return derivative;
459  }
460  }
461 }
462 
463 
464 
465 template <int dim>
466 inline Tensor<1, dim>
468  const unsigned int i,
469  const Point<dim> & p) const
470 {
471  return compute_derivative<1>(i, p);
472 }
473 
474 
475 
476 template <int dim>
477 inline Tensor<2, dim>
479  const unsigned int i,
480  const Point<dim> & p) const
481 {
482  return compute_derivative<2>(i, p);
483 }
484 
485 
486 
487 template <int dim>
488 inline Tensor<3, dim>
490  const unsigned int i,
491  const Point<dim> & p) const
492 {
493  return compute_derivative<3>(i, p);
494 }
495 
496 
497 
498 template <int dim>
499 inline Tensor<4, dim>
501  const unsigned int i,
502  const Point<dim> & p) const
503 {
504  return compute_derivative<4>(i, p);
505 }
506 
507 #endif // DOXYGEN
509 
510 #endif
Definition: point.h:110
void evaluate(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 override
Tensor< order, dim > compute_derivative(const unsigned int i, const Point< dim > &p) const
const std::vector< unsigned int > & get_numbering_inverse() const
std::string name() const override
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
TensorProductPolynomialsBubbles(const std::vector< Pol > &pols)
virtual Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
void set_numbering(const std::vector< unsigned int > &renumber)
virtual Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
double compute_value(const unsigned int i, const Point< dim > &p) const override
const std::vector< unsigned int > & get_numbering() const
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
virtual Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
virtual Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T fixed_power(const T t)
Definition: utilities.h:966
static const unsigned int invalid_unsigned_int
Definition: types.h:213