16 #ifndef dealii_tensor_product_polynomials_bubbles_h
17 #define dealii_tensor_product_polynomials_bubbles_h
81 const std::vector<unsigned int> &
87 const std::vector<unsigned int> &
104 std::vector<double> & values,
185 name()
const override;
190 virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
191 clone()
const override;
220 const std::vector<Pol> &pols)
224 , index_map(tensor_polys.n() +
225 ((tensor_polys.polynomials.size() <= 2) ? 1 : dim))
226 , index_map_inverse(tensor_polys.n() +
227 ((tensor_polys.polynomials.size() <= 2) ? 1 : dim))
230 const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
232 for (
unsigned int i = 0; i <
tensor_polys.
n() + n_bubbles; ++i)
257 inline const std::vector<unsigned int> &
265 inline const std::vector<unsigned int> &
276 return "TensorProductPolynomialsBubbles";
284 const unsigned int i,
289 Assert(i < max_q_indices + ((q_degree <= 1) ? 1 : dim),
293 if (i < max_q_indices)
294 return tensor_polys.template compute_derivative<order>(i, p);
306 for (
unsigned int d = 0;
d < dim; ++
d)
308 derivative_1[
d] = 1.;
310 for (
unsigned j = 0; j < dim; ++j)
312 (
d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
314 for (
unsigned int i = 0; i < q_degree - 1; ++i)
315 derivative_1[
d] *= 2 * p(comp) - 1;
322 for (
unsigned int j = 0; j < dim; ++j)
323 value *= 4 * p(j) * (1 - p(j));
325 double tmp =
value * 2 * (q_degree - 1);
326 for (
unsigned int i = 0; i < q_degree - 2; ++i)
327 tmp *= 2 * p(comp) - 1;
328 derivative_1[comp] += tmp;
338 double v[dim + 1][3];
340 for (
unsigned int c = 0; c < dim; ++c)
342 v[c][0] = 4 * p(c) * (1 - p(c));
343 v[c][1] = 4 * (1 - 2 * p(c));
348 for (
unsigned int i = 0; i < q_degree - 1; ++i)
349 tmp *= 2 * p(comp) - 1;
354 double tmp = 2 * (q_degree - 1);
355 for (
unsigned int i = 0; i < q_degree - 2; ++i)
356 tmp *= 2 * p(comp) - 1;
364 double tmp = 4 * (q_degree - 2) * (q_degree - 1);
365 for (
unsigned int i = 0; i < q_degree - 3; ++i)
366 tmp *= 2 * p(comp) - 1;
375 for (
unsigned int d1 = 0; d1 < dim; ++d1)
376 for (
unsigned int d2 = 0; d2 < dim; ++d2)
378 grad_grad_1[d1][d2] = v[dim][0];
379 for (
unsigned int x = 0; x < dim; ++x)
381 unsigned int derivative = 0;
382 if (d1 == x || d2 == x)
389 grad_grad_1[d1][d2] *= v[x][derivative];
397 for (
unsigned int d = 0;
d < dim; ++
d)
399 grad_grad_2[
d][comp] = v[dim][1];
400 grad_grad_3[comp][
d] = v[dim][1];
401 for (
unsigned int x = 0; x < dim; ++x)
403 grad_grad_2[
d][comp] *= v[x][
d == x];
404 grad_grad_3[comp][
d] *= v[x][
d == x];
409 double psi_value = 1.;
410 for (
unsigned int x = 0; x < dim; ++x)
411 psi_value *= v[x][0];
413 for (
unsigned int d1 = 0; d1 < dim; ++d1)
414 for (
unsigned int d2 = 0; d2 < dim; ++d2)
415 derivative_2[d1][d2] =
416 grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
417 derivative_2[comp][comp] += psi_value * v[dim][2];