34 std::array<unsigned int, dim> ix;
35 for (
unsigned int i = 0; i < tensor_polys.n(); ++i)
37 tensor_polys.compute_index(i, ix);
39 for (
unsigned int d = 0; d < dim; ++d)
50 const std::vector<unsigned int> &renumber)
52 Assert(renumber.size() == index_map.size(),
56 for (
unsigned int i = 0; i < index_map.size(); ++i)
57 index_map_inverse[index_map[i]] = i;
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]);
63 tensor_polys.set_numbering(renumber_base);
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 + ((q_degree <= 1) ? 1 : dim),
78 if (i < max_q_indices)
79 return tensor_polys.compute_value(i, p);
81 const unsigned int comp = i - tensor_polys.n();
85 for (
unsigned int j = 0; j < dim; ++j)
86 value *= 4 * p(j) * (1 - p(j));
90 for (
unsigned int i = 0; i < q_degree - 1; ++i)
91 value *= (2 * p(comp) - 1);
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 + ((q_degree <= 1) ? 1 : dim),
119 if (i < max_q_indices)
120 return tensor_polys.compute_grad(i, p);
122 const unsigned int comp = i - tensor_polys.n();
125 for (
unsigned int d = 0; d < dim; ++d)
129 for (
unsigned j = 0; j < dim; ++j)
130 grad[d] *= (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
132 for (
unsigned int i = 0; i < q_degree - 1; ++i)
133 grad[d] *= 2 * p(comp) - 1;
140 for (
unsigned int j = 0; j < dim; ++j)
141 value *= 4 * p(j) * (1 - p(j));
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;
157 const unsigned int i,
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 + ((q_degree <= 1) ? 1 : dim),
166 if (i < max_q_indices)
167 return tensor_polys.compute_grad_grad(i, p);
169 const unsigned int comp = i - tensor_polys.n();
171 double v[dim + 1][3];
173 for (
unsigned int c = 0; c < dim; ++c)
175 v[c][0] = 4 * p(c) * (1 - p(c));
176 v[c][1] = 4 * (1 - 2 * p(c));
181 for (
unsigned int i = 0; i < q_degree - 1; ++i)
182 tmp *= 2 * p(comp) - 1;
187 double tmp = 2 * (q_degree - 1);
188 for (
unsigned int i = 0; i < q_degree - 2; ++i)
189 tmp *= 2 * p(comp) - 1;
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;
208 for (
unsigned int d1 = 0; d1 < dim; ++d1)
209 for (
unsigned int d2 = 0; d2 < dim; ++d2)
211 grad_grad_1[d1][d2] = v[dim][0];
212 for (
unsigned int x = 0; x < dim; ++x)
214 unsigned int derivative = 0;
215 if (d1 == x || d2 == x)
222 grad_grad_1[d1][d2] *= v[x][derivative];
230 for (
unsigned int d = 0; d < dim; ++d)
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)
236 grad_grad_2[d][comp] *= v[x][d == x];
237 grad_grad_3[comp][d] *= v[x][d == x];
243 double psi_value = 1.;
244 for (
unsigned int x = 0; x < dim; ++x)
245 psi_value *= v[x][0];
247 for (
unsigned int d1 = 0; d1 < dim; ++d1)
248 for (
unsigned int d2 = 0; d2 < dim; ++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];
262 std::vector<double> & values,
268 const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
269 const unsigned int max_q_indices = tensor_polys.n();
271 const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
272 Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0,
274 Assert(grads.size() == max_q_indices + n_bubbles || grads.size() == 0,
277 grad_grads.size() == max_q_indices + n_bubbles || grad_grads.size() == 0,
279 Assert(third_derivatives.size() == max_q_indices + n_bubbles ||
280 third_derivatives.size() == 0,
282 max_q_indices + n_bubbles,
284 Assert(fourth_derivatives.size() == max_q_indices + n_bubbles ||
285 fourth_derivatives.size() == 0,
287 max_q_indices + n_bubbles,
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)
294 values.resize(tensor_polys.n());
297 if (grads.empty() ==
false)
299 grads.resize(tensor_polys.n());
302 if (grad_grads.empty() ==
false)
304 grad_grads.resize(tensor_polys.n());
305 do_grad_grads =
true;
307 if (third_derivatives.empty() ==
false)
309 third_derivatives.resize(tensor_polys.n());
310 do_3rd_derivatives =
true;
312 if (fourth_derivatives.empty() ==
false)
314 fourth_derivatives.resize(tensor_polys.n());
315 do_4th_derivatives =
true;
318 tensor_polys.evaluate(
319 p, values, grads, grad_grads, third_derivatives, fourth_derivatives);
321 for (
unsigned int i = tensor_polys.n(); i < tensor_polys.n() + n_bubbles; ++i)
324 values.push_back(compute_value(i, p));
326 grads.push_back(compute_grad(i, p));
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));
339std::unique_ptr<ScalarPolynomialsBase<dim>>
342 return std::make_unique<TensorProductPolynomialsBubbles<dim>>(*this);
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)
void output_indices(std::ostream &out) const
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
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)