17 #include <deal.II/base/exceptions.h> 18 #include <deal.II/base/table.h> 19 #include <deal.II/base/tensor_product_polynomials_bubbles.h> 21 DEAL_II_NAMESPACE_OPEN
33 const unsigned int q_degree = this->polynomials.size() - 1;
34 const unsigned int max_q_indices = this->n_tensor_pols;
35 const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
40 if (i < max_q_indices)
43 const unsigned int comp = i - this->n_tensor_pols;
47 for (
unsigned int j = 0; j < dim; ++j)
48 value *= 4 * p(j) * (1 - p(j));
50 for (
unsigned int i = 0; i < q_degree - 1; ++i)
51 value *= 2 * p(comp) - 1;
72 const unsigned int q_degree = this->polynomials.size() - 1;
73 const unsigned int max_q_indices = this->n_tensor_pols;
74 const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
79 if (i < max_q_indices)
82 const unsigned int comp = i - this->n_tensor_pols;
85 for (
unsigned int d = 0; d < dim; ++d)
89 for (
unsigned j = 0; j < dim; ++j)
90 grad[d] *= (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
92 for (
unsigned int i = 0; i < q_degree - 1; ++i)
93 grad[d] *= 2 * p(comp) - 1;
100 for (
unsigned int j = 0; j < dim; ++j)
101 value *= 4 * p(j) * (1 - p(j));
103 double tmp = value * 2 * (q_degree - 1);
104 for (
unsigned int i = 0; i < q_degree - 2; ++i)
105 tmp *= 2 * p(comp) - 1;
117 const unsigned int i,
120 const unsigned int q_degree = this->polynomials.size() - 1;
121 const unsigned int max_q_indices = this->n_tensor_pols;
122 const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
127 if (i < max_q_indices)
130 const unsigned int comp = i - this->n_tensor_pols;
132 double v[dim + 1][3];
134 for (
unsigned int c = 0; c < dim; ++c)
136 v[c][0] = 4 * p(c) * (1 - p(c));
137 v[c][1] = 4 * (1 - 2 * p(c));
142 for (
unsigned int i = 0; i < q_degree - 1; ++i)
143 tmp *= 2 * p(comp) - 1;
148 double tmp = 2 * (q_degree - 1);
149 for (
unsigned int i = 0; i < q_degree - 2; ++i)
150 tmp *= 2 * p(comp) - 1;
158 double tmp = 4 * (q_degree - 2) * (q_degree - 1);
159 for (
unsigned int i = 0; i < q_degree - 3; ++i)
160 tmp *= 2 * p(comp) - 1;
169 for (
unsigned int d1 = 0; d1 < dim; ++d1)
170 for (
unsigned int d2 = 0; d2 < dim; ++d2)
172 grad_grad_1[d1][d2] = v[dim][0];
173 for (
unsigned int x = 0; x < dim; ++x)
175 unsigned int derivative = 0;
176 if (d1 == x || d2 == x)
183 grad_grad_1[d1][d2] *= v[x][derivative];
191 for (
unsigned int d = 0; d < dim; ++d)
193 grad_grad_2[d][comp] = v[dim][1];
194 grad_grad_3[comp][d] = v[dim][1];
195 for (
unsigned int x = 0; x < dim; ++x)
197 grad_grad_2[d][comp] *= v[x][d == x];
198 grad_grad_3[comp][d] *= v[x][d == x];
204 double psi_value = 1.;
205 for (
unsigned int x = 0; x < dim; ++x)
206 psi_value *= v[x][0];
208 for (
unsigned int d1 = 0; d1 < dim; ++d1)
209 for (
unsigned int d2 = 0; d2 < dim; ++d2)
211 grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
212 grad_grad[comp][comp] += psi_value * v[dim][2];
221 std::vector<double> & values,
227 const unsigned int q_degree = this->polynomials.size() - 1;
228 const unsigned int max_q_indices = this->n_tensor_pols;
230 const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
231 Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0,
233 Assert(grads.size() == max_q_indices + n_bubbles || grads.size() == 0,
236 grad_grads.size() == max_q_indices + n_bubbles || grad_grads.size() == 0,
238 Assert(third_derivatives.size() == max_q_indices + n_bubbles ||
239 third_derivatives.size() == 0,
241 max_q_indices + n_bubbles,
243 Assert(fourth_derivatives.size() == max_q_indices + n_bubbles ||
244 fourth_derivatives.size() == 0,
246 max_q_indices + n_bubbles,
249 bool do_values =
false, do_grads =
false, do_grad_grads =
false;
250 bool do_3rd_derivatives =
false, do_4th_derivatives =
false;
251 if (values.empty() ==
false)
253 values.resize(this->n_tensor_pols);
256 if (grads.empty() ==
false)
258 grads.resize(this->n_tensor_pols);
261 if (grad_grads.empty() ==
false)
263 grad_grads.resize(this->n_tensor_pols);
264 do_grad_grads =
true;
266 if (third_derivatives.empty() ==
false)
268 third_derivatives.resize(this->n_tensor_pols);
269 do_3rd_derivatives =
true;
271 if (fourth_derivatives.empty() ==
false)
273 fourth_derivatives.resize(this->n_tensor_pols);
274 do_4th_derivatives =
true;
278 p, values, grads, grad_grads, third_derivatives, fourth_derivatives);
280 for (
unsigned int i = this->n_tensor_pols;
281 i < this->n_tensor_pols + n_bubbles;
285 values.push_back(compute_value(i, p));
287 grads.push_back(compute_grad(i, p));
289 grad_grads.push_back(compute_grad_grad(i, p));
290 if (do_3rd_derivatives)
291 third_derivatives.push_back(compute_derivative<3>(i, p));
292 if (do_4th_derivatives)
293 fourth_derivatives.push_back(compute_derivative<4>(i, p));
303 DEAL_II_NAMESPACE_CLOSE
double compute_value(const unsigned int i, const Point< dim > &p) const
double compute_value(const unsigned int i, const Point< dim > &p) const
static::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) 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 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
#define Assert(cond, exc)
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 & ExcNotImplemented()
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
static::ExceptionBase & ExcInternalError()