Reference documentation for deal.II version GIT ffb4c3937f 2023-03-31 14:25:02+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.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2012 - 2021 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*x_i*(1-x_i)
84  double value = 1.;
85  for (unsigned int j = 0; j < dim; ++j)
86  value *= 4 * p(j) * (1 - p(j));
87 
88  // Then multiply with (2x_i-1)^{r-1}. Since q_degree is generally a
89  // small integer, using a loop is likely faster than using std::pow.
90  for (unsigned int i = 0; i < q_degree - 1; ++i)
91  value *= (2 * p(comp) - 1);
92  return value;
93 }
94 
95 
96 
97 template <>
98 double
100  const Point<0> &) const
101 {
102  Assert(false, ExcNotImplemented());
103  return 0.;
104 }
105 
106 
107 
108 template <int dim>
111  const Point<dim> & p) const
112 {
113  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
114  const unsigned int max_q_indices = tensor_polys.n();
115  Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
116  ExcInternalError());
117 
118  // treat the regular basis functions
119  if (i < max_q_indices)
120  return tensor_polys.compute_grad(i, p);
121 
122  const unsigned int comp = i - tensor_polys.n();
123  Tensor<1, dim> grad;
124 
125  for (unsigned int d = 0; d < dim; ++d)
126  {
127  grad[d] = 1.;
128  // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
129  for (unsigned j = 0; j < dim; ++j)
130  grad[d] *= (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
131  // and multiply with (2*x_i-1)^{r-1}
132  for (unsigned int i = 0; i < q_degree - 1; ++i)
133  grad[d] *= 2 * p(comp) - 1;
134  }
135 
136  if (q_degree >= 2)
137  {
138  // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
139  double value = 1.;
140  for (unsigned int j = 0; j < dim; ++j)
141  value *= 4 * p(j) * (1 - p(j));
142  // and multiply with grad(2*x_i-1)^{r-1}
143  double tmp = value * 2 * (q_degree - 1);
144  for (unsigned int i = 0; i < q_degree - 2; ++i)
145  tmp *= 2 * p(comp) - 1;
146  grad[comp] += tmp;
147  }
148 
149  return grad;
150 }
151 
152 
153 
154 template <int dim>
157  const unsigned int i,
158  const Point<dim> & p) const
159 {
160  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
161  const unsigned int max_q_indices = tensor_polys.n();
162  Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
163  ExcInternalError());
164 
165  // treat the regular basis functions
166  if (i < max_q_indices)
167  return tensor_polys.compute_grad_grad(i, p);
168 
169  const unsigned int comp = i - tensor_polys.n();
170 
171  double v[dim + 1][3];
172  {
173  for (unsigned int c = 0; c < dim; ++c)
174  {
175  v[c][0] = 4 * p(c) * (1 - p(c));
176  v[c][1] = 4 * (1 - 2 * p(c));
177  v[c][2] = -8;
178  }
179 
180  double tmp = 1.;
181  for (unsigned int i = 0; i < q_degree - 1; ++i)
182  tmp *= 2 * p(comp) - 1;
183  v[dim][0] = tmp;
184 
185  if (q_degree >= 2)
186  {
187  double tmp = 2 * (q_degree - 1);
188  for (unsigned int i = 0; i < q_degree - 2; ++i)
189  tmp *= 2 * p(comp) - 1;
190  v[dim][1] = tmp;
191  }
192  else
193  v[dim][1] = 0.;
194 
195  if (q_degree >= 3)
196  {
197  double tmp = 4 * (q_degree - 2) * (q_degree - 1);
198  for (unsigned int i = 0; i < q_degree - 3; ++i)
199  tmp *= 2 * p(comp) - 1;
200  v[dim][2] = tmp;
201  }
202  else
203  v[dim][2] = 0.;
204  }
205 
206  // calculate (\partial_j \partial_k \psi) * monomial
207  Tensor<2, dim> grad_grad_1;
208  for (unsigned int d1 = 0; d1 < dim; ++d1)
209  for (unsigned int d2 = 0; d2 < dim; ++d2)
210  {
211  grad_grad_1[d1][d2] = v[dim][0];
212  for (unsigned int x = 0; x < dim; ++x)
213  {
214  unsigned int derivative = 0;
215  if (d1 == x || d2 == x)
216  {
217  if (d1 == d2)
218  derivative = 2;
219  else
220  derivative = 1;
221  }
222  grad_grad_1[d1][d2] *= v[x][derivative];
223  }
224  }
225 
226  // calculate (\partial_j \psi) *(\partial_k monomial)
227  // and (\partial_k \psi) *(\partial_j monomial)
228  Tensor<2, dim> grad_grad_2;
229  Tensor<2, dim> grad_grad_3;
230  for (unsigned int d = 0; d < dim; ++d)
231  {
232  grad_grad_2[d][comp] = v[dim][1];
233  grad_grad_3[comp][d] = v[dim][1];
234  for (unsigned int x = 0; x < dim; ++x)
235  {
236  grad_grad_2[d][comp] *= v[x][d == x];
237  grad_grad_3[comp][d] *= v[x][d == x];
238  }
239  }
240 
241  // calculate \psi *(\partial j \partial_k monomial) and sum
242  Tensor<2, dim> grad_grad;
243  double psi_value = 1.;
244  for (unsigned int x = 0; x < dim; ++x)
245  psi_value *= v[x][0];
246 
247  for (unsigned int d1 = 0; d1 < dim; ++d1)
248  for (unsigned int d2 = 0; d2 < dim; ++d2)
249  grad_grad[d1][d2] =
250  grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
251  grad_grad[comp][comp] += psi_value * v[dim][2];
252 
253  return grad_grad;
254 }
255 
256 
257 
258 template <int dim>
259 void
261  const Point<dim> & p,
262  std::vector<double> & values,
263  std::vector<Tensor<1, dim>> &grads,
264  std::vector<Tensor<2, dim>> &grad_grads,
265  std::vector<Tensor<3, dim>> &third_derivatives,
266  std::vector<Tensor<4, dim>> &fourth_derivatives) const
267 {
268  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
269  const unsigned int max_q_indices = tensor_polys.n();
270  (void)max_q_indices;
271  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
272  Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0,
273  ExcDimensionMismatch2(values.size(), max_q_indices + n_bubbles, 0));
274  Assert(grads.size() == max_q_indices + n_bubbles || grads.size() == 0,
275  ExcDimensionMismatch2(grads.size(), max_q_indices + n_bubbles, 0));
276  Assert(
277  grad_grads.size() == max_q_indices + n_bubbles || grad_grads.size() == 0,
278  ExcDimensionMismatch2(grad_grads.size(), max_q_indices + n_bubbles, 0));
279  Assert(third_derivatives.size() == max_q_indices + n_bubbles ||
280  third_derivatives.size() == 0,
281  ExcDimensionMismatch2(third_derivatives.size(),
282  max_q_indices + n_bubbles,
283  0));
284  Assert(fourth_derivatives.size() == max_q_indices + n_bubbles ||
285  fourth_derivatives.size() == 0,
286  ExcDimensionMismatch2(fourth_derivatives.size(),
287  max_q_indices + n_bubbles,
288  0));
289 
290  bool do_values = false, do_grads = false, do_grad_grads = false;
291  bool do_3rd_derivatives = false, do_4th_derivatives = false;
292  if (values.empty() == false)
293  {
294  values.resize(tensor_polys.n());
295  do_values = true;
296  }
297  if (grads.empty() == false)
298  {
299  grads.resize(tensor_polys.n());
300  do_grads = true;
301  }
302  if (grad_grads.empty() == false)
303  {
304  grad_grads.resize(tensor_polys.n());
305  do_grad_grads = true;
306  }
307  if (third_derivatives.empty() == false)
308  {
309  third_derivatives.resize(tensor_polys.n());
310  do_3rd_derivatives = true;
311  }
312  if (fourth_derivatives.empty() == false)
313  {
314  fourth_derivatives.resize(tensor_polys.n());
315  do_4th_derivatives = true;
316  }
317 
318  tensor_polys.evaluate(
319  p, values, grads, grad_grads, third_derivatives, fourth_derivatives);
320 
321  for (unsigned int i = tensor_polys.n(); i < tensor_polys.n() + n_bubbles; ++i)
322  {
323  if (do_values)
324  values.push_back(compute_value(i, p));
325  if (do_grads)
326  grads.push_back(compute_grad(i, p));
327  if (do_grad_grads)
328  grad_grads.push_back(compute_grad_grad(i, p));
329  if (do_3rd_derivatives)
330  third_derivatives.push_back(compute_derivative<3>(i, p));
331  if (do_4th_derivatives)
332  fourth_derivatives.push_back(compute_derivative<4>(i, p));
333  }
334 }
335 
336 
337 
338 template <int dim>
339 std::unique_ptr<ScalarPolynomialsBase<dim>>
341 {
342  return std::make_unique<TensorProductPolynomialsBubbles<dim>>(*this);
343 }
344 
345 
346 /* ------------------- explicit instantiations -------------- */
350 
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< 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
void set_numbering(const std::vector< unsigned int > &renumber)
double compute_value(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)