Reference documentation for deal.II version Git 8596a7cd07 2020-12-04 07:30:43 +0100
\(\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.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2012 - 2019 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 
19 
20 #include <memory>
21 
23 
24 
25 
26 /* ------------------- TensorProductPolynomialsBubbles -------------- */
27 
28 
29 
30 template <int dim>
31 void
33 {
34  std::array<unsigned int, dim> ix;
35  for (unsigned int i = 0; i < tensor_polys.n(); ++i)
36  {
37  tensor_polys.compute_index(i, ix);
38  out << i << "\t";
39  for (unsigned int d = 0; d < dim; ++d)
40  out << ix[d] << " ";
41  out << std::endl;
42  }
43 }
44 
45 
46 
47 template <int dim>
48 void
50  const std::vector<unsigned int> &renumber)
51 {
52  Assert(renumber.size() == index_map.size(),
53  ExcDimensionMismatch(renumber.size(), index_map.size()));
54 
55  index_map = renumber;
56  for (unsigned int i = 0; i < index_map.size(); ++i)
57  index_map_inverse[index_map[i]] = i;
58 
59  std::vector<unsigned int> renumber_base;
60  for (unsigned int i = 0; i < tensor_polys.n(); ++i)
61  renumber_base.push_back(renumber[i]);
62 
63  tensor_polys.set_numbering(renumber_base);
64 }
65 
66 
67 template <int dim>
68 double
70  const Point<dim> & p) const
71 {
72  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
73  const unsigned int max_q_indices = tensor_polys.n();
74  Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
76 
77  // treat the regular basis functions
78  if (i < max_q_indices)
79  return tensor_polys.compute_value(i, p);
80 
81  const unsigned int comp = i - tensor_polys.n();
82 
83  // compute \prod_{i=1}^d 4*(1-x_i^2)(p)
84  double value = 1.;
85  for (unsigned int j = 0; j < dim; ++j)
86  value *= 4 * p(j) * (1 - p(j));
87  // and multiply with (2x_i-1)^{r-1}
88  for (unsigned int i = 0; i < q_degree - 1; ++i)
89  value *= 2 * p(comp) - 1;
90  return value;
91 }
92 
93 
94 
95 template <>
96 double
98  const Point<0> &) const
99 {
100  Assert(false, ExcNotImplemented());
101  return 0.;
102 }
103 
104 
105 
106 template <int dim>
109  const Point<dim> & p) const
110 {
111  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
112  const unsigned int max_q_indices = tensor_polys.n();
113  Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
114  ExcInternalError());
115 
116  // treat the regular basis functions
117  if (i < max_q_indices)
118  return tensor_polys.compute_grad(i, p);
119 
120  const unsigned int comp = i - tensor_polys.n();
121  Tensor<1, dim> grad;
122 
123  for (unsigned int d = 0; d < dim; ++d)
124  {
125  grad[d] = 1.;
126  // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
127  for (unsigned j = 0; j < dim; ++j)
128  grad[d] *= (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
129  // and multiply with (2*x_i-1)^{r-1}
130  for (unsigned int i = 0; i < q_degree - 1; ++i)
131  grad[d] *= 2 * p(comp) - 1;
132  }
133 
134  if (q_degree >= 2)
135  {
136  // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
137  double value = 1.;
138  for (unsigned int j = 0; j < dim; ++j)
139  value *= 4 * p(j) * (1 - p(j));
140  // and multiply with grad(2*x_i-1)^{r-1}
141  double tmp = value * 2 * (q_degree - 1);
142  for (unsigned int i = 0; i < q_degree - 2; ++i)
143  tmp *= 2 * p(comp) - 1;
144  grad[comp] += tmp;
145  }
146 
147  return grad;
148 }
149 
150 
151 
152 template <int dim>
155  const unsigned int i,
156  const Point<dim> & p) const
157 {
158  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
159  const unsigned int max_q_indices = tensor_polys.n();
160  Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
161  ExcInternalError());
162 
163  // treat the regular basis functions
164  if (i < max_q_indices)
165  return tensor_polys.compute_grad_grad(i, p);
166 
167  const unsigned int comp = i - tensor_polys.n();
168 
169  double v[dim + 1][3];
170  {
171  for (unsigned int c = 0; c < dim; ++c)
172  {
173  v[c][0] = 4 * p(c) * (1 - p(c));
174  v[c][1] = 4 * (1 - 2 * p(c));
175  v[c][2] = -8;
176  }
177 
178  double tmp = 1.;
179  for (unsigned int i = 0; i < q_degree - 1; ++i)
180  tmp *= 2 * p(comp) - 1;
181  v[dim][0] = tmp;
182 
183  if (q_degree >= 2)
184  {
185  double tmp = 2 * (q_degree - 1);
186  for (unsigned int i = 0; i < q_degree - 2; ++i)
187  tmp *= 2 * p(comp) - 1;
188  v[dim][1] = tmp;
189  }
190  else
191  v[dim][1] = 0.;
192 
193  if (q_degree >= 3)
194  {
195  double tmp = 4 * (q_degree - 2) * (q_degree - 1);
196  for (unsigned int i = 0; i < q_degree - 3; ++i)
197  tmp *= 2 * p(comp) - 1;
198  v[dim][2] = tmp;
199  }
200  else
201  v[dim][2] = 0.;
202  }
203 
204  // calculate (\partial_j \partial_k \psi) * monomial
205  Tensor<2, dim> grad_grad_1;
206  for (unsigned int d1 = 0; d1 < dim; ++d1)
207  for (unsigned int d2 = 0; d2 < dim; ++d2)
208  {
209  grad_grad_1[d1][d2] = v[dim][0];
210  for (unsigned int x = 0; x < dim; ++x)
211  {
212  unsigned int derivative = 0;
213  if (d1 == x || d2 == x)
214  {
215  if (d1 == d2)
216  derivative = 2;
217  else
218  derivative = 1;
219  }
220  grad_grad_1[d1][d2] *= v[x][derivative];
221  }
222  }
223 
224  // calculate (\partial_j \psi) *(\partial_k monomial)
225  // and (\partial_k \psi) *(\partial_j monomial)
226  Tensor<2, dim> grad_grad_2;
227  Tensor<2, dim> grad_grad_3;
228  for (unsigned int d = 0; d < dim; ++d)
229  {
230  grad_grad_2[d][comp] = v[dim][1];
231  grad_grad_3[comp][d] = v[dim][1];
232  for (unsigned int x = 0; x < dim; ++x)
233  {
234  grad_grad_2[d][comp] *= v[x][d == x];
235  grad_grad_3[comp][d] *= v[x][d == x];
236  }
237  }
238 
239  // calculate \psi *(\partial j \partial_k monomial) and sum
240  Tensor<2, dim> grad_grad;
241  double psi_value = 1.;
242  for (unsigned int x = 0; x < dim; ++x)
243  psi_value *= v[x][0];
244 
245  for (unsigned int d1 = 0; d1 < dim; ++d1)
246  for (unsigned int d2 = 0; d2 < dim; ++d2)
247  grad_grad[d1][d2] =
248  grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
249  grad_grad[comp][comp] += psi_value * v[dim][2];
250 
251  return grad_grad;
252 }
253 
254 
255 
256 template <int dim>
257 void
259  const Point<dim> & p,
260  std::vector<double> & values,
261  std::vector<Tensor<1, dim>> &grads,
262  std::vector<Tensor<2, dim>> &grad_grads,
263  std::vector<Tensor<3, dim>> &third_derivatives,
264  std::vector<Tensor<4, dim>> &fourth_derivatives) const
265 {
266  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
267  const unsigned int max_q_indices = tensor_polys.n();
268  (void)max_q_indices;
269  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
270  Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0,
271  ExcDimensionMismatch2(values.size(), max_q_indices + n_bubbles, 0));
272  Assert(grads.size() == max_q_indices + n_bubbles || grads.size() == 0,
273  ExcDimensionMismatch2(grads.size(), max_q_indices + n_bubbles, 0));
274  Assert(
275  grad_grads.size() == max_q_indices + n_bubbles || grad_grads.size() == 0,
276  ExcDimensionMismatch2(grad_grads.size(), max_q_indices + n_bubbles, 0));
277  Assert(third_derivatives.size() == max_q_indices + n_bubbles ||
278  third_derivatives.size() == 0,
279  ExcDimensionMismatch2(third_derivatives.size(),
280  max_q_indices + n_bubbles,
281  0));
282  Assert(fourth_derivatives.size() == max_q_indices + n_bubbles ||
283  fourth_derivatives.size() == 0,
284  ExcDimensionMismatch2(fourth_derivatives.size(),
285  max_q_indices + n_bubbles,
286  0));
287 
288  bool do_values = false, do_grads = false, do_grad_grads = false;
289  bool do_3rd_derivatives = false, do_4th_derivatives = false;
290  if (values.empty() == false)
291  {
292  values.resize(tensor_polys.n());
293  do_values = true;
294  }
295  if (grads.empty() == false)
296  {
297  grads.resize(tensor_polys.n());
298  do_grads = true;
299  }
300  if (grad_grads.empty() == false)
301  {
302  grad_grads.resize(tensor_polys.n());
303  do_grad_grads = true;
304  }
305  if (third_derivatives.empty() == false)
306  {
307  third_derivatives.resize(tensor_polys.n());
308  do_3rd_derivatives = true;
309  }
310  if (fourth_derivatives.empty() == false)
311  {
312  fourth_derivatives.resize(tensor_polys.n());
313  do_4th_derivatives = true;
314  }
315 
316  tensor_polys.evaluate(
317  p, values, grads, grad_grads, third_derivatives, fourth_derivatives);
318 
319  for (unsigned int i = tensor_polys.n(); i < tensor_polys.n() + n_bubbles; ++i)
320  {
321  if (do_values)
322  values.push_back(compute_value(i, p));
323  if (do_grads)
324  grads.push_back(compute_grad(i, p));
325  if (do_grad_grads)
326  grad_grads.push_back(compute_grad_grad(i, p));
327  if (do_3rd_derivatives)
328  third_derivatives.push_back(compute_derivative<3>(i, p));
329  if (do_4th_derivatives)
330  fourth_derivatives.push_back(compute_derivative<4>(i, p));
331  }
332 }
333 
334 
335 
336 template <int dim>
337 std::unique_ptr<ScalarPolynomialsBase<dim>>
339 {
340  return std::make_unique<TensorProductPolynomialsBubbles<dim>>(*this);
341 }
342 
343 
344 /* ------------------- explicit instantiations -------------- */
348 
static ::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
double compute_value(const unsigned int i, const Point< dim > &p) const override
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
void set_numbering(const std::vector< unsigned int > &renumber)
#define Assert(cond, exc)
Definition: exceptions.h:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
static ::ExceptionBase & ExcNotImplemented()
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
static ::ExceptionBase & ExcInternalError()