16 #include <deal.II/base/exceptions.h>
17 #include <deal.II/base/polynomials_piecewise.h>
18 #include <deal.II/base/tensor_product_polynomials.h>
20 #include <boost/container/small_vector.hpp>
24 DEAL_II_NAMESPACE_OPEN
37 compute_tensor_index(
const unsigned int,
40 unsigned int (&)[dim])
46 compute_tensor_index(
const unsigned int n,
49 unsigned int (&indices)[1])
55 compute_tensor_index(
const unsigned int n,
56 const unsigned int n_pols_0,
58 unsigned int (&indices)[2])
60 indices[0] = n % n_pols_0;
61 indices[1] = n / n_pols_0;
65 compute_tensor_index(
const unsigned int n,
66 const unsigned int n_pols_0,
67 const unsigned int n_pols_1,
68 unsigned int (&indices)[3])
70 indices[0] = n % n_pols_0;
71 indices[1] = (n / n_pols_0) % n_pols_1;
72 indices[2] = n / (n_pols_0 * n_pols_1);
79 template <
int dim,
typename PolynomialType>
83 unsigned int (&indices)[(dim > 0 ? dim : 1)])
const
85 Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
87 internal::compute_tensor_index(index_map[i],
95 template <
int dim,
typename PolynomialType>
98 std::ostream &out)
const
100 unsigned int ix[dim];
101 for (
unsigned int i = 0; i < n_tensor_pols; ++i)
103 compute_index(i, ix);
105 for (
unsigned int d = 0; d < dim; ++d)
113 template <
int dim,
typename PolynomialType>
116 const std::vector<unsigned int> &renumber)
118 Assert(renumber.size() == index_map.size(),
121 index_map = renumber;
122 for (
unsigned int i = 0; i < index_map.size(); ++i)
123 index_map_inverse[index_map[i]] = i;
140 template <
int dim,
typename PolynomialType>
143 const unsigned int i,
148 unsigned int indices[dim];
149 compute_index(i, indices);
152 for (
unsigned int d = 0; d < dim; ++d)
153 value *= polynomials[indices[d]].value(p(d));
160 template <
int dim,
typename PolynomialType>
163 const unsigned int i,
166 unsigned int indices[dim];
167 compute_index(i, indices);
175 std::vector<double> tmp(2);
176 for (
unsigned int d = 0; d < dim; ++d)
178 polynomials[indices[d]].value(p(d), tmp);
185 for (
unsigned int d = 0; d < dim; ++d)
188 for (
unsigned int x = 0; x < dim; ++x)
189 grad[d] *= v[x][d == x];
197 template <
int dim,
typename PolynomialType>
200 const unsigned int i,
203 unsigned int indices[dim];
204 compute_index(i, indices);
208 std::vector<double> tmp(3);
209 for (
unsigned int d = 0; d < dim; ++d)
211 polynomials[indices[d]].value(p(d), tmp);
219 for (
unsigned int d1 = 0; d1 < dim; ++d1)
220 for (
unsigned int d2 = 0; d2 < dim; ++d2)
222 grad_grad[d1][d2] = 1.;
223 for (
unsigned int x = 0; x < dim; ++x)
225 unsigned int derivative = 0;
226 if (d1 == x || d2 == x)
233 grad_grad[d1][d2] *= v[x][derivative];
242 template <
int dim,
typename PolynomialType>
246 std::vector<double> & values,
253 Assert(values.size() == n_tensor_pols || values.size() == 0,
255 Assert(grads.size() == n_tensor_pols || grads.size() == 0,
257 Assert(grad_grads.size() == n_tensor_pols || grad_grads.size() == 0,
259 Assert(third_derivatives.size() == n_tensor_pols ||
260 third_derivatives.size() == 0,
262 Assert(fourth_derivatives.size() == n_tensor_pols ||
263 fourth_derivatives.size() == 0,
266 const bool update_values = (values.size() == n_tensor_pols),
267 update_grads = (grads.size() == n_tensor_pols),
268 update_grad_grads = (grad_grads.size() == n_tensor_pols),
269 update_3rd_derivatives =
270 (third_derivatives.size() == n_tensor_pols),
271 update_4th_derivatives =
272 (fourth_derivatives.size() == n_tensor_pols);
275 unsigned int n_values_and_derivatives = 0;
277 n_values_and_derivatives = 1;
279 n_values_and_derivatives = 2;
280 if (update_grad_grads)
281 n_values_and_derivatives = 3;
282 if (update_3rd_derivatives)
283 n_values_and_derivatives = 4;
284 if (update_4th_derivatives)
285 n_values_and_derivatives = 5;
292 const unsigned int n_polynomials = polynomials.size();
293 boost::container::small_vector<std::array<std::array<double, 5>, dim>, 20>
294 values_1d(n_polynomials);
295 if (n_values_and_derivatives == 1)
296 for (
unsigned int i = 0; i < n_polynomials; ++i)
297 for (
unsigned int d = 0; d < dim; ++d)
298 values_1d[i][d][0] = polynomials[i].value(p(d));
300 for (
unsigned int i = 0; i < n_polynomials; ++i)
301 for (
unsigned d = 0; d < dim; ++d)
302 polynomials[i].value(p(d),
303 n_values_and_derivatives,
304 values_1d[i][d].data());
306 unsigned int indices[3];
307 unsigned int ind = 0;
308 for (indices[2] = 0; indices[2] < (dim > 2 ? n_polynomials : 1); ++indices[2])
309 for (indices[1] = 0; indices[1] < (dim > 1 ? n_polynomials : 1);
311 if (n_values_and_derivatives == 1)
312 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
314 double value = values_1d[indices[0]][0][0];
315 for (
unsigned int d = 1; d < dim; ++d)
316 value *= values_1d[indices[d]][d][0];
317 values[index_map_inverse[ind]] = value;
320 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
322 unsigned int i = index_map_inverse[ind];
326 double value = values_1d[indices[0]][0][0];
327 for (
unsigned int x = 1; x < dim; ++x)
328 value *= values_1d[indices[x]][x][0];
333 for (
unsigned int d = 0; d < dim; ++d)
336 for (
unsigned int x = 0; x < dim; ++x)
337 grad *= values_1d[indices[x]][x][(d == x) ? 1 : 0];
341 if (update_grad_grads)
342 for (
unsigned int d1 = 0; d1 < dim; ++d1)
343 for (
unsigned int d2 = 0; d2 < dim; ++d2)
346 for (
unsigned int x = 0; x < dim; ++x)
348 unsigned int derivative = 0;
354 der2 *= values_1d[indices[x]][x][derivative];
356 grad_grads[i][d1][d2] = der2;
359 if (update_3rd_derivatives)
360 for (
unsigned int d1 = 0; d1 < dim; ++d1)
361 for (
unsigned int d2 = 0; d2 < dim; ++d2)
362 for (
unsigned int d3 = 0; d3 < dim; ++d3)
365 for (
unsigned int x = 0; x < dim; ++x)
367 unsigned int derivative = 0;
375 der3 *= values_1d[indices[x]][x][derivative];
377 third_derivatives[i][d1][d2][d3] = der3;
380 if (update_4th_derivatives)
381 for (
unsigned int d1 = 0; d1 < dim; ++d1)
382 for (
unsigned int d2 = 0; d2 < dim; ++d2)
383 for (
unsigned int d3 = 0; d3 < dim; ++d3)
384 for (
unsigned int d4 = 0; d4 < dim; ++d4)
387 for (
unsigned int x = 0; x < dim; ++x)
389 unsigned int derivative = 0;
399 der4 *= values_1d[indices[x]][x][derivative];
401 fourth_derivatives[i][d1][d2][d3][d4] = der4;
415 , n_tensor_pols(get_n_tensor_pols(pols))
418 for (
unsigned int d = 0; d < dim; ++d)
419 Assert(pols[d].size() > 0,
420 ExcMessage(
"The number of polynomials must be larger than zero "
421 "for all coordinate directions."));
429 unsigned int (&indices)[dim])
const
432 unsigned int n_poly = 1;
433 for (
unsigned int d = 0; d < dim; ++d)
434 n_poly *= polynomials[d].size();
439 internal::compute_tensor_index(i,
440 polynomials[0].size(),
444 internal::compute_tensor_index(i,
445 polynomials[0].size(),
446 polynomials[1].size(),
457 unsigned int indices[dim];
458 compute_index(i, indices);
461 for (
unsigned int d = 0; d < dim; ++d)
462 value *= polynomials[d][indices[d]].value(p(d));
473 unsigned int indices[dim];
474 compute_index(i, indices);
480 std::vector<std::vector<double>> v(dim, std::vector<double>(2));
481 for (
unsigned int d = 0; d < dim; ++d)
482 polynomials[d][indices[d]].value(p(d), v[d]);
485 for (
unsigned int d = 0; d < dim; ++d)
488 for (
unsigned int x = 0; x < dim; ++x)
489 grad[d] *= v[x][d == x];
501 unsigned int indices[dim];
502 compute_index(i, indices);
504 std::vector<std::vector<double>> v(dim, std::vector<double>(3));
505 for (
unsigned int d = 0; d < dim; ++d)
506 polynomials[d][indices[d]].value(p(d), v[d]);
509 for (
unsigned int d1 = 0; d1 < dim; ++d1)
510 for (
unsigned int d2 = 0; d2 < dim; ++d2)
512 grad_grad[d1][d2] = 1.;
513 for (
unsigned int x = 0; x < dim; ++x)
515 unsigned int derivative = 0;
516 if (d1 == x || d2 == x)
523 grad_grad[d1][d2] *= v[x][derivative];
536 std::vector<double> & values,
542 Assert(values.size() == n_tensor_pols || values.size() == 0,
544 Assert(grads.size() == n_tensor_pols || grads.size() == 0,
546 Assert(grad_grads.size() == n_tensor_pols || grad_grads.size() == 0,
548 Assert(third_derivatives.size() == n_tensor_pols ||
549 third_derivatives.size() == 0,
551 Assert(fourth_derivatives.size() == n_tensor_pols ||
552 fourth_derivatives.size() == 0,
555 const bool update_values = (values.size() == n_tensor_pols),
556 update_grads = (grads.size() == n_tensor_pols),
557 update_grad_grads = (grad_grads.size() == n_tensor_pols),
558 update_3rd_derivatives =
559 (third_derivatives.size() == n_tensor_pols),
560 update_4th_derivatives =
561 (fourth_derivatives.size() == n_tensor_pols);
566 unsigned int n_values_and_derivatives = 0;
568 n_values_and_derivatives = 1;
570 n_values_and_derivatives = 2;
571 if (update_grad_grads)
572 n_values_and_derivatives = 3;
573 if (update_3rd_derivatives)
574 n_values_and_derivatives = 4;
575 if (update_4th_derivatives)
576 n_values_and_derivatives = 5;
582 std::vector<std::vector<std::vector<double>>> v(dim);
583 for (
unsigned int d = 0; d < dim; ++d)
585 v[d].resize(polynomials[d].size());
586 for (
unsigned int i = 0; i < polynomials[d].size(); ++i)
588 v[d][i].resize(n_values_and_derivatives, 0.);
589 polynomials[d][i].value(p(d), v[d][i]);
593 for (
unsigned int i = 0; i < n_tensor_pols; ++i)
599 unsigned int indices[dim];
600 compute_index(i, indices);
605 for (
unsigned int x = 0; x < dim; ++x)
606 values[i] *= v[x][indices[x]][0];
610 for (
unsigned int d = 0; d < dim; ++d)
613 for (
unsigned int x = 0; x < dim; ++x)
614 grads[i][d] *= v[x][indices[x]][d == x ? 1 : 0];
617 if (update_grad_grads)
618 for (
unsigned int d1 = 0; d1 < dim; ++d1)
619 for (
unsigned int d2 = 0; d2 < dim; ++d2)
621 grad_grads[i][d1][d2] = 1.;
622 for (
unsigned int x = 0; x < dim; ++x)
624 unsigned int derivative = 0;
630 grad_grads[i][d1][d2] *= v[x][indices[x]][derivative];
634 if (update_3rd_derivatives)
635 for (
unsigned int d1 = 0; d1 < dim; ++d1)
636 for (
unsigned int d2 = 0; d2 < dim; ++d2)
637 for (
unsigned int d3 = 0; d3 < dim; ++d3)
639 third_derivatives[i][d1][d2][d3] = 1.;
640 for (
unsigned int x = 0; x < dim; ++x)
642 unsigned int derivative = 0;
650 third_derivatives[i][d1][d2][d3] *=
651 v[x][indices[x]][derivative];
655 if (update_4th_derivatives)
656 for (
unsigned int d1 = 0; d1 < dim; ++d1)
657 for (
unsigned int d2 = 0; d2 < dim; ++d2)
658 for (
unsigned int d3 = 0; d3 < dim; ++d3)
659 for (
unsigned int d4 = 0; d4 < dim; ++d4)
661 fourth_derivatives[i][d1][d2][d3][d4] = 1.;
662 for (
unsigned int x = 0; x < dim; ++x)
664 unsigned int derivative = 0;
674 fourth_derivatives[i][d1][d2][d3][d4] *=
675 v[x][indices[x]][derivative];
687 return n_tensor_pols;
697 for (
unsigned int d = 0; d < dim; ++d)
723 DEAL_II_NAMESPACE_CLOSE