Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
tensor_product_polynomials_bubbles.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2012 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/exceptions.h>
18 #include <deal.II/base/tensor_product_polynomials_bubbles.h>
19 
20 DEAL_II_NAMESPACE_OPEN
21 
22 
23 
24 /* ------------------- TensorProductPolynomialsBubbles -------------- */
25 
26 
27 template <int dim>
28 double
30  const Point<dim> & p) const
31 {
32  const unsigned int q_degree = this->polynomials.size() - 1;
33  const unsigned int max_q_indices = this->n_tensor_pols;
34  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
35  (void)n_bubbles;
36  Assert(i < max_q_indices + n_bubbles, ExcInternalError());
37 
38  // treat the regular basis functions
39  if (i < max_q_indices)
41 
42  const unsigned int comp = i - this->n_tensor_pols;
43 
44  // compute \prod_{i=1}^d 4*(1-x_i^2)(p)
45  double value = 1.;
46  for (unsigned int j = 0; j < dim; ++j)
47  value *= 4 * p(j) * (1 - p(j));
48  // and multiply with (2x_i-1)^{r-1}
49  for (unsigned int i = 0; i < q_degree - 1; ++i)
50  value *= 2 * p(comp) - 1;
51  return value;
52 }
53 
54 
55 
56 template <>
57 double
59  const Point<0> &) const
60 {
61  Assert(false, ExcNotImplemented());
62  return 0.;
63 }
64 
65 
66 template <int dim>
69  const Point<dim> & p) const
70 {
71  const unsigned int q_degree = this->polynomials.size() - 1;
72  const unsigned int max_q_indices = this->n_tensor_pols;
73  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
74  (void)n_bubbles;
75  Assert(i < max_q_indices + n_bubbles, ExcInternalError());
76 
77  // treat the regular basis functions
78  if (i < max_q_indices)
80 
81  const unsigned int comp = i - this->n_tensor_pols;
82  Tensor<1, dim> grad;
83 
84  for (unsigned int d = 0; d < dim; ++d)
85  {
86  grad[d] = 1.;
87  // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
88  for (unsigned j = 0; j < dim; ++j)
89  grad[d] *= (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
90  // and multiply with (2*x_i-1)^{r-1}
91  for (unsigned int i = 0; i < q_degree - 1; ++i)
92  grad[d] *= 2 * p(comp) - 1;
93  }
94 
95  if (q_degree >= 2)
96  {
97  // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
98  double value = 1.;
99  for (unsigned int j = 0; j < dim; ++j)
100  value *= 4 * p(j) * (1 - p(j));
101  // and multiply with grad(2*x_i-1)^{r-1}
102  double tmp = value * 2 * (q_degree - 1);
103  for (unsigned int i = 0; i < q_degree - 2; ++i)
104  tmp *= 2 * p(comp) - 1;
105  grad[comp] += tmp;
106  }
107 
108  return grad;
109 }
110 
111 
112 
113 template <int dim>
116  const unsigned int i,
117  const Point<dim> & p) const
118 {
119  const unsigned int q_degree = this->polynomials.size() - 1;
120  const unsigned int max_q_indices = this->n_tensor_pols;
121  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
122  (void)n_bubbles;
123  Assert(i < max_q_indices + n_bubbles, ExcInternalError());
124 
125  // treat the regular basis functions
126  if (i < max_q_indices)
128 
129  const unsigned int comp = i - this->n_tensor_pols;
130 
131  double v[dim + 1][3];
132  {
133  for (unsigned int c = 0; c < dim; ++c)
134  {
135  v[c][0] = 4 * p(c) * (1 - p(c));
136  v[c][1] = 4 * (1 - 2 * p(c));
137  v[c][2] = -8;
138  }
139 
140  double tmp = 1.;
141  for (unsigned int i = 0; i < q_degree - 1; ++i)
142  tmp *= 2 * p(comp) - 1;
143  v[dim][0] = tmp;
144 
145  if (q_degree >= 2)
146  {
147  double tmp = 2 * (q_degree - 1);
148  for (unsigned int i = 0; i < q_degree - 2; ++i)
149  tmp *= 2 * p(comp) - 1;
150  v[dim][1] = tmp;
151  }
152  else
153  v[dim][1] = 0.;
154 
155  if (q_degree >= 3)
156  {
157  double tmp = 4 * (q_degree - 2) * (q_degree - 1);
158  for (unsigned int i = 0; i < q_degree - 3; ++i)
159  tmp *= 2 * p(comp) - 1;
160  v[dim][2] = tmp;
161  }
162  else
163  v[dim][2] = 0.;
164  }
165 
166  // calculate (\partial_j \partial_k \psi) * monomial
167  Tensor<2, dim> grad_grad_1;
168  for (unsigned int d1 = 0; d1 < dim; ++d1)
169  for (unsigned int d2 = 0; d2 < dim; ++d2)
170  {
171  grad_grad_1[d1][d2] = v[dim][0];
172  for (unsigned int x = 0; x < dim; ++x)
173  {
174  unsigned int derivative = 0;
175  if (d1 == x || d2 == x)
176  {
177  if (d1 == d2)
178  derivative = 2;
179  else
180  derivative = 1;
181  }
182  grad_grad_1[d1][d2] *= v[x][derivative];
183  }
184  }
185 
186  // calculate (\partial_j \psi) *(\partial_k monomial)
187  // and (\partial_k \psi) *(\partial_j monomial)
188  Tensor<2, dim> grad_grad_2;
189  Tensor<2, dim> grad_grad_3;
190  for (unsigned int d = 0; d < dim; ++d)
191  {
192  grad_grad_2[d][comp] = v[dim][1];
193  grad_grad_3[comp][d] = v[dim][1];
194  for (unsigned int x = 0; x < dim; ++x)
195  {
196  grad_grad_2[d][comp] *= v[x][d == x];
197  grad_grad_3[comp][d] *= v[x][d == x];
198  }
199  }
200 
201  // calculate \psi *(\partial j \partial_k monomial) and sum
202  Tensor<2, dim> grad_grad;
203  double psi_value = 1.;
204  for (unsigned int x = 0; x < dim; ++x)
205  psi_value *= v[x][0];
206 
207  for (unsigned int d1 = 0; d1 < dim; ++d1)
208  for (unsigned int d2 = 0; d2 < dim; ++d2)
209  grad_grad[d1][d2] =
210  grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
211  grad_grad[comp][comp] += psi_value * v[dim][2];
212 
213  return grad_grad;
214 }
215 
216 template <int dim>
217 void
219  const Point<dim> & p,
220  std::vector<double> & values,
221  std::vector<Tensor<1, dim>> &grads,
222  std::vector<Tensor<2, dim>> &grad_grads,
223  std::vector<Tensor<3, dim>> &third_derivatives,
224  std::vector<Tensor<4, dim>> &fourth_derivatives) const
225 {
226  const unsigned int q_degree = this->polynomials.size() - 1;
227  const unsigned int max_q_indices = this->n_tensor_pols;
228  (void)max_q_indices;
229  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
230  Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0,
231  ExcDimensionMismatch2(values.size(), max_q_indices + n_bubbles, 0));
232  Assert(grads.size() == max_q_indices + n_bubbles || grads.size() == 0,
233  ExcDimensionMismatch2(grads.size(), max_q_indices + n_bubbles, 0));
234  Assert(
235  grad_grads.size() == max_q_indices + n_bubbles || grad_grads.size() == 0,
236  ExcDimensionMismatch2(grad_grads.size(), max_q_indices + n_bubbles, 0));
237  Assert(third_derivatives.size() == max_q_indices + n_bubbles ||
238  third_derivatives.size() == 0,
239  ExcDimensionMismatch2(third_derivatives.size(),
240  max_q_indices + n_bubbles,
241  0));
242  Assert(fourth_derivatives.size() == max_q_indices + n_bubbles ||
243  fourth_derivatives.size() == 0,
244  ExcDimensionMismatch2(fourth_derivatives.size(),
245  max_q_indices + n_bubbles,
246  0));
247 
248  bool do_values = false, do_grads = false, do_grad_grads = false;
249  bool do_3rd_derivatives = false, do_4th_derivatives = false;
250  if (values.empty() == false)
251  {
252  values.resize(this->n_tensor_pols);
253  do_values = true;
254  }
255  if (grads.empty() == false)
256  {
257  grads.resize(this->n_tensor_pols);
258  do_grads = true;
259  }
260  if (grad_grads.empty() == false)
261  {
262  grad_grads.resize(this->n_tensor_pols);
263  do_grad_grads = true;
264  }
265  if (third_derivatives.empty() == false)
266  {
267  third_derivatives.resize(this->n_tensor_pols);
268  do_3rd_derivatives = true;
269  }
270  if (fourth_derivatives.empty() == false)
271  {
272  fourth_derivatives.resize(this->n_tensor_pols);
273  do_4th_derivatives = true;
274  }
275 
277  p, values, grads, grad_grads, third_derivatives, fourth_derivatives);
278 
279  for (unsigned int i = this->n_tensor_pols;
280  i < this->n_tensor_pols + n_bubbles;
281  ++i)
282  {
283  if (do_values)
284  values.push_back(compute_value(i, p));
285  if (do_grads)
286  grads.push_back(compute_grad(i, p));
287  if (do_grad_grads)
288  grad_grads.push_back(compute_grad_grad(i, p));
289  if (do_3rd_derivatives)
290  third_derivatives.push_back(compute_derivative<3>(i, p));
291  if (do_4th_derivatives)
292  fourth_derivatives.push_back(compute_derivative<4>(i, p));
293  }
294 }
295 
296 
297 /* ------------------- explicit instantiations -------------- */
301 
302 DEAL_II_NAMESPACE_CLOSE
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
static ::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
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
#define Assert(cond, exc)
Definition: exceptions.h:1407
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
double compute_value(const unsigned int i, const Point< dim > &p) const
static ::ExceptionBase & ExcNotImplemented()
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
double compute_value(const unsigned int i, const Point< dim > &p) const
static ::ExceptionBase & ExcInternalError()