16 #include <deal.II/base/exceptions.h> 17 #include <deal.II/base/polynomials_piecewise.h> 18 #include <deal.II/base/table.h> 19 #include <deal.II/base/tensor_product_polynomials.h> 21 #include <boost/container/small_vector.hpp> 25 DEAL_II_NAMESPACE_OPEN
38 compute_tensor_index(
const unsigned int,
41 unsigned int (&)[dim])
47 compute_tensor_index(
const unsigned int n,
50 unsigned int (&indices)[1])
56 compute_tensor_index(
const unsigned int n,
57 const unsigned int n_pols_0,
59 unsigned int (&indices)[2])
61 indices[0] = n % n_pols_0;
62 indices[1] = n / n_pols_0;
66 compute_tensor_index(
const unsigned int n,
67 const unsigned int n_pols_0,
68 const unsigned int n_pols_1,
69 unsigned int (&indices)[3])
71 indices[0] = n % n_pols_0;
72 indices[1] = (n / n_pols_0) % n_pols_1;
73 indices[2] = n / (n_pols_0 * n_pols_1);
80 template <
int dim,
typename PolynomialType>
84 unsigned int (&indices)[(dim > 0 ? dim : 1)])
const 86 Assert(i < Utilities::fixed_power<dim>(polynomials.size()),
88 internal::compute_tensor_index(index_map[i],
96 template <
int dim,
typename PolynomialType>
99 std::ostream &out)
const 101 unsigned int ix[dim];
102 for (
unsigned int i = 0; i < n_tensor_pols; ++i)
104 compute_index(i, ix);
106 for (
unsigned int d = 0; d < dim; ++d)
114 template <
int dim,
typename PolynomialType>
117 const std::vector<unsigned int> &renumber)
119 Assert(renumber.size() == index_map.size(),
122 index_map = renumber;
123 for (
unsigned int i = 0; i < index_map.size(); ++i)
124 index_map_inverse[index_map[i]] = i;
141 template <
int dim,
typename PolynomialType>
144 const unsigned int i,
149 unsigned int indices[dim];
150 compute_index(i, indices);
153 for (
unsigned int d = 0; d < dim; ++d)
154 value *= polynomials[indices[d]].value(p(d));
161 template <
int dim,
typename PolynomialType>
164 const unsigned int i,
167 unsigned int indices[dim];
168 compute_index(i, indices);
176 std::vector<double> tmp(2);
177 for (
unsigned int d = 0; d < dim; ++d)
179 polynomials[indices[d]].value(p(d), tmp);
186 for (
unsigned int d = 0; d < dim; ++d)
189 for (
unsigned int x = 0; x < dim; ++x)
190 grad[d] *= v[x][d == x];
198 template <
int dim,
typename PolynomialType>
201 const unsigned int i,
204 unsigned int indices[dim];
205 compute_index(i, indices);
209 std::vector<double> tmp(3);
210 for (
unsigned int d = 0; d < dim; ++d)
212 polynomials[indices[d]].value(p(d), tmp);
220 for (
unsigned int d1 = 0; d1 < dim; ++d1)
221 for (
unsigned int d2 = 0; d2 < dim; ++d2)
223 grad_grad[d1][d2] = 1.;
224 for (
unsigned int x = 0; x < dim; ++x)
226 unsigned int derivative = 0;
227 if (d1 == x || d2 == x)
234 grad_grad[d1][d2] *= v[x][derivative];
243 template <
int dim,
typename PolynomialType>
247 std::vector<double> & values,
254 Assert(values.size() == n_tensor_pols || values.size() == 0,
256 Assert(grads.size() == n_tensor_pols || grads.size() == 0,
258 Assert(grad_grads.size() == n_tensor_pols || grad_grads.size() == 0,
260 Assert(third_derivatives.size() == n_tensor_pols ||
261 third_derivatives.size() == 0,
263 Assert(fourth_derivatives.size() == n_tensor_pols ||
264 fourth_derivatives.size() == 0,
268 update_grads = (grads.size() == n_tensor_pols),
269 update_grad_grads = (grad_grads.size() == n_tensor_pols),
271 (third_derivatives.size() == n_tensor_pols),
272 update_4th_derivatives =
273 (fourth_derivatives.size() == n_tensor_pols);
276 unsigned int n_values_and_derivatives = 0;
278 n_values_and_derivatives = 1;
280 n_values_and_derivatives = 2;
281 if (update_grad_grads)
282 n_values_and_derivatives = 3;
284 n_values_and_derivatives = 4;
285 if (update_4th_derivatives)
286 n_values_and_derivatives = 5;
293 const unsigned int n_polynomials = polynomials.size();
294 boost::container::small_vector<std::array<std::array<double, 5>, dim>, 20>
295 values_1d(n_polynomials);
296 if (n_values_and_derivatives == 1)
297 for (
unsigned int i = 0; i < n_polynomials; ++i)
298 for (
unsigned int d = 0; d < dim; ++d)
299 values_1d[i][d][0] = polynomials[i].value(p(d));
301 for (
unsigned int i = 0; i < n_polynomials; ++i)
302 for (
unsigned d = 0; d < dim; ++d)
303 polynomials[i].value(p(d),
304 n_values_and_derivatives,
305 &values_1d[i][d][0]);
307 unsigned int indices[3];
308 unsigned int ind = 0;
309 for (indices[2] = 0; indices[2] < (dim > 2 ? n_polynomials : 1); ++indices[2])
310 for (indices[1] = 0; indices[1] < (dim > 1 ? n_polynomials : 1);
312 if (n_values_and_derivatives == 1)
313 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
315 double value = values_1d[indices[0]][0][0];
316 for (
unsigned int d = 1; d < dim; ++d)
317 value *= values_1d[indices[d]][d][0];
318 values[index_map_inverse[ind]] = value;
321 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
323 unsigned int i = index_map_inverse[ind];
327 double value = values_1d[indices[0]][0][0];
328 for (
unsigned int x = 1; x < dim; ++x)
329 value *= values_1d[indices[x]][x][0];
334 for (
unsigned int d = 0; d < dim; ++d)
337 for (
unsigned int x = 0; x < dim; ++x)
338 grad *= values_1d[indices[x]][x][(d == x) ? 1 : 0];
342 if (update_grad_grads)
343 for (
unsigned int d1 = 0; d1 < dim; ++d1)
344 for (
unsigned int d2 = 0; d2 < dim; ++d2)
347 for (
unsigned int x = 0; x < dim; ++x)
349 unsigned int derivative = 0;
355 der2 *= values_1d[indices[x]][x][derivative];
357 grad_grads[i][d1][d2] = der2;
361 for (
unsigned int d1 = 0; d1 < dim; ++d1)
362 for (
unsigned int d2 = 0; d2 < dim; ++d2)
363 for (
unsigned int d3 = 0; d3 < dim; ++d3)
366 for (
unsigned int x = 0; x < dim; ++x)
368 unsigned int derivative = 0;
376 der3 *= values_1d[indices[x]][x][derivative];
378 third_derivatives[i][d1][d2][d3] = der3;
381 if (update_4th_derivatives)
382 for (
unsigned int d1 = 0; d1 < dim; ++d1)
383 for (
unsigned int d2 = 0; d2 < dim; ++d2)
384 for (
unsigned int d3 = 0; d3 < dim; ++d3)
385 for (
unsigned int d4 = 0; d4 < dim; ++d4)
388 for (
unsigned int x = 0; x < dim; ++x)
390 unsigned int derivative = 0;
400 der4 *= values_1d[indices[x]][x][derivative];
402 fourth_derivatives[i][d1][d2][d3][d4] = der4;
416 , n_tensor_pols(get_n_tensor_pols(pols))
419 for (
unsigned int d = 0; d < dim; ++d)
420 Assert(pols[d].size() > 0,
421 ExcMessage(
"The number of polynomials must be larger than zero " 422 "for all coordinate directions."));
430 unsigned int (&indices)[dim])
const 433 unsigned int n_poly = 1;
434 for (
unsigned int d = 0; d < dim; ++d)
440 internal::compute_tensor_index(i,
445 internal::compute_tensor_index(i,
458 unsigned int indices[dim];
462 for (
unsigned int d = 0; d < dim; ++d)
474 unsigned int indices[dim];
481 std::vector<std::vector<double>> v(dim, std::vector<double>(2));
482 for (
unsigned int d = 0; d < dim; ++d)
486 for (
unsigned int d = 0; d < dim; ++d)
489 for (
unsigned int x = 0; x < dim; ++x)
490 grad[d] *= v[x][d == x];
502 unsigned int indices[dim];
505 std::vector<std::vector<double>> v(dim, std::vector<double>(3));
506 for (
unsigned int d = 0; d < dim; ++d)
510 for (
unsigned int d1 = 0; d1 < dim; ++d1)
511 for (
unsigned int d2 = 0; d2 < dim; ++d2)
513 grad_grad[d1][d2] = 1.;
514 for (
unsigned int x = 0; x < dim; ++x)
516 unsigned int derivative = 0;
517 if (d1 == x || d2 == x)
524 grad_grad[d1][d2] *= v[x][derivative];
537 std::vector<double> & values,
550 third_derivatives.size() == 0,
553 fourth_derivatives.size() == 0,
561 update_4th_derivatives =
567 unsigned int n_values_and_derivatives = 0;
569 n_values_and_derivatives = 1;
571 n_values_and_derivatives = 2;
572 if (update_grad_grads)
573 n_values_and_derivatives = 3;
575 n_values_and_derivatives = 4;
576 if (update_4th_derivatives)
577 n_values_and_derivatives = 5;
583 std::vector<std::vector<std::vector<double>>> v(dim);
584 for (
unsigned int d = 0; d < dim; ++d)
587 for (
unsigned int i = 0; i <
polynomials[d].size(); ++i)
589 v[d][i].resize(n_values_and_derivatives, 0.);
600 unsigned int indices[dim];
606 for (
unsigned int x = 0; x < dim; ++x)
607 values[i] *= v[x][indices[x]][0];
611 for (
unsigned int d = 0; d < dim; ++d)
614 for (
unsigned int x = 0; x < dim; ++x)
615 grads[i][d] *= v[x][indices[x]][d == x ? 1 : 0];
618 if (update_grad_grads)
619 for (
unsigned int d1 = 0; d1 < dim; ++d1)
620 for (
unsigned int d2 = 0; d2 < dim; ++d2)
622 grad_grads[i][d1][d2] = 1.;
623 for (
unsigned int x = 0; x < dim; ++x)
625 unsigned int derivative = 0;
631 grad_grads[i][d1][d2] *= v[x][indices[x]][derivative];
636 for (
unsigned int d1 = 0; d1 < dim; ++d1)
637 for (
unsigned int d2 = 0; d2 < dim; ++d2)
638 for (
unsigned int d3 = 0; d3 < dim; ++d3)
640 third_derivatives[i][d1][d2][d3] = 1.;
641 for (
unsigned int x = 0; x < dim; ++x)
643 unsigned int derivative = 0;
651 third_derivatives[i][d1][d2][d3] *=
652 v[x][indices[x]][derivative];
656 if (update_4th_derivatives)
657 for (
unsigned int d1 = 0; d1 < dim; ++d1)
658 for (
unsigned int d2 = 0; d2 < dim; ++d2)
659 for (
unsigned int d3 = 0; d3 < dim; ++d3)
660 for (
unsigned int d4 = 0; d4 < dim; ++d4)
662 fourth_derivatives[i][d1][d2][d3][d4] = 1.;
663 for (
unsigned int x = 0; x < dim; ++x)
665 unsigned int derivative = 0;
675 fourth_derivatives[i][d1][d2][d3][d4] *=
676 v[x][indices[x]][derivative];
698 for (
unsigned int d = 0; d < dim; ++d)
715 Polynomials::PiecewisePolynomial<double>>;
718 Polynomials::PiecewisePolynomial<double>>;
724 DEAL_II_NAMESPACE_CLOSE
double compute_value(const unsigned int i, const Point< dim > &p) const
static::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
void compute_index(const unsigned int i, unsigned int(&indices)[(dim > 0?dim:1)]) const
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
void set_numbering(const std::vector< unsigned int > &renumber)
static unsigned int get_n_tensor_pols(const std::vector< std::vector< Polynomials::Polynomial< double >>> &pols)
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
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
static::ExceptionBase & ExcMessage(std::string arg1)
Third derivatives of shape functions.
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const std::vector< std::vector< Polynomials::Polynomial< double > > > polynomials
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
void output_indices(std::ostream &out) const
static::ExceptionBase & ExcNotImplemented()
AnisotropicPolynomials(const std::vector< std::vector< Polynomials::Polynomial< double >>> &base_polynomials)
const unsigned int n_tensor_pols
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
void compute_index(const unsigned int i, unsigned int(&indices)[dim]) const
double compute_value(const unsigned int i, const Point< dim > &p) const
static::ExceptionBase & ExcInternalError()