Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor_product_polynomials_bubbles.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2012 - 2019 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 
20 
22 
23 
24 
25 /* ------------------- TensorProductPolynomialsBubbles -------------- */
26 
27 
28 
29 template <int dim>
30 void
32 {
33  std::array<unsigned int, dim> ix;
34  for (unsigned int i = 0; i < tensor_polys.n(); ++i)
35  {
37  out << i << "\t";
38  for (unsigned int d = 0; d < dim; ++d)
39  out << ix[d] << " ";
40  out << std::endl;
41  }
42 }
43 
44 
45 
46 template <int dim>
47 void
49  const std::vector<unsigned int> &renumber)
50 {
51  Assert(renumber.size() == index_map.size(),
52  ExcDimensionMismatch(renumber.size(), index_map.size()));
53 
54  index_map = renumber;
55  for (unsigned int i = 0; i < index_map.size(); ++i)
57 
58  std::vector<unsigned int> renumber_base;
59  for (unsigned int i = 0; i < tensor_polys.n(); ++i)
60  renumber_base.push_back(renumber[i]);
61 
62  tensor_polys.set_numbering(renumber_base);
63 }
64 
65 
66 template <int dim>
67 double
69  const Point<dim> & p) const
70 {
71  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
72  const unsigned int max_q_indices = tensor_polys.n();
73  Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
75 
76  // treat the regular basis functions
77  if (i < max_q_indices)
78  return tensor_polys.compute_value(i, p);
79 
80  const unsigned int comp = i - tensor_polys.n();
81 
82  // compute \prod_{i=1}^d 4*(1-x_i^2)(p)
83  double value = 1.;
84  for (unsigned int j = 0; j < dim; ++j)
85  value *= 4 * p(j) * (1 - p(j));
86  // and multiply with (2x_i-1)^{r-1}
87  for (unsigned int i = 0; i < q_degree - 1; ++i)
88  value *= 2 * p(comp) - 1;
89  return value;
90 }
91 
92 
93 
94 template <>
95 double
97  const Point<0> &) const
98 {
99  Assert(false, ExcNotImplemented());
100  return 0.;
101 }
102 
103 
104 
105 template <int dim>
108  const Point<dim> & p) const
109 {
110  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
111  const unsigned int max_q_indices = tensor_polys.n();
112  Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
113  ExcInternalError());
114 
115  // treat the regular basis functions
116  if (i < max_q_indices)
117  return tensor_polys.compute_grad(i, p);
118 
119  const unsigned int comp = i - tensor_polys.n();
120  Tensor<1, dim> grad;
121 
122  for (unsigned int d = 0; d < dim; ++d)
123  {
124  grad[d] = 1.;
125  // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
126  for (unsigned j = 0; j < dim; ++j)
127  grad[d] *= (d == j ? 4 * (1 - 2 * p(j)) : 4 * p(j) * (1 - p(j)));
128  // and multiply with (2*x_i-1)^{r-1}
129  for (unsigned int i = 0; i < q_degree - 1; ++i)
130  grad[d] *= 2 * p(comp) - 1;
131  }
132 
133  if (q_degree >= 2)
134  {
135  // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
136  double value = 1.;
137  for (unsigned int j = 0; j < dim; ++j)
138  value *= 4 * p(j) * (1 - p(j));
139  // and multiply with grad(2*x_i-1)^{r-1}
140  double tmp = value * 2 * (q_degree - 1);
141  for (unsigned int i = 0; i < q_degree - 2; ++i)
142  tmp *= 2 * p(comp) - 1;
143  grad[comp] += tmp;
144  }
145 
146  return grad;
147 }
148 
149 
150 
151 template <int dim>
154  const unsigned int i,
155  const Point<dim> & p) const
156 {
157  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
158  const unsigned int max_q_indices = tensor_polys.n();
159  Assert(i < max_q_indices + /* n_bubbles= */ ((q_degree <= 1) ? 1 : dim),
160  ExcInternalError());
161 
162  // treat the regular basis functions
163  if (i < max_q_indices)
164  return tensor_polys.compute_grad_grad(i, p);
165 
166  const unsigned int comp = i - tensor_polys.n();
167 
168  double v[dim + 1][3];
169  {
170  for (unsigned int c = 0; c < dim; ++c)
171  {
172  v[c][0] = 4 * p(c) * (1 - p(c));
173  v[c][1] = 4 * (1 - 2 * p(c));
174  v[c][2] = -8;
175  }
176 
177  double tmp = 1.;
178  for (unsigned int i = 0; i < q_degree - 1; ++i)
179  tmp *= 2 * p(comp) - 1;
180  v[dim][0] = tmp;
181 
182  if (q_degree >= 2)
183  {
184  double tmp = 2 * (q_degree - 1);
185  for (unsigned int i = 0; i < q_degree - 2; ++i)
186  tmp *= 2 * p(comp) - 1;
187  v[dim][1] = tmp;
188  }
189  else
190  v[dim][1] = 0.;
191 
192  if (q_degree >= 3)
193  {
194  double tmp = 4 * (q_degree - 2) * (q_degree - 1);
195  for (unsigned int i = 0; i < q_degree - 3; ++i)
196  tmp *= 2 * p(comp) - 1;
197  v[dim][2] = tmp;
198  }
199  else
200  v[dim][2] = 0.;
201  }
202 
203  // calculate (\partial_j \partial_k \psi) * monomial
204  Tensor<2, dim> grad_grad_1;
205  for (unsigned int d1 = 0; d1 < dim; ++d1)
206  for (unsigned int d2 = 0; d2 < dim; ++d2)
207  {
208  grad_grad_1[d1][d2] = v[dim][0];
209  for (unsigned int x = 0; x < dim; ++x)
210  {
211  unsigned int derivative = 0;
212  if (d1 == x || d2 == x)
213  {
214  if (d1 == d2)
215  derivative = 2;
216  else
217  derivative = 1;
218  }
219  grad_grad_1[d1][d2] *= v[x][derivative];
220  }
221  }
222 
223  // calculate (\partial_j \psi) *(\partial_k monomial)
224  // and (\partial_k \psi) *(\partial_j monomial)
225  Tensor<2, dim> grad_grad_2;
226  Tensor<2, dim> grad_grad_3;
227  for (unsigned int d = 0; d < dim; ++d)
228  {
229  grad_grad_2[d][comp] = v[dim][1];
230  grad_grad_3[comp][d] = v[dim][1];
231  for (unsigned int x = 0; x < dim; ++x)
232  {
233  grad_grad_2[d][comp] *= v[x][d == x];
234  grad_grad_3[comp][d] *= v[x][d == x];
235  }
236  }
237 
238  // calculate \psi *(\partial j \partial_k monomial) and sum
239  Tensor<2, dim> grad_grad;
240  double psi_value = 1.;
241  for (unsigned int x = 0; x < dim; ++x)
242  psi_value *= v[x][0];
243 
244  for (unsigned int d1 = 0; d1 < dim; ++d1)
245  for (unsigned int d2 = 0; d2 < dim; ++d2)
246  grad_grad[d1][d2] =
247  grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
248  grad_grad[comp][comp] += psi_value * v[dim][2];
249 
250  return grad_grad;
251 }
252 
253 
254 
255 template <int dim>
256 void
258  const Point<dim> & p,
259  std::vector<double> & values,
260  std::vector<Tensor<1, dim>> &grads,
261  std::vector<Tensor<2, dim>> &grad_grads,
262  std::vector<Tensor<3, dim>> &third_derivatives,
263  std::vector<Tensor<4, dim>> &fourth_derivatives) const
264 {
265  const unsigned int q_degree = tensor_polys.polynomials.size() - 1;
266  const unsigned int max_q_indices = tensor_polys.n();
267  (void)max_q_indices;
268  const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim);
269  Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0,
270  ExcDimensionMismatch2(values.size(), max_q_indices + n_bubbles, 0));
271  Assert(grads.size() == max_q_indices + n_bubbles || grads.size() == 0,
272  ExcDimensionMismatch2(grads.size(), max_q_indices + n_bubbles, 0));
273  Assert(
274  grad_grads.size() == max_q_indices + n_bubbles || grad_grads.size() == 0,
275  ExcDimensionMismatch2(grad_grads.size(), max_q_indices + n_bubbles, 0));
276  Assert(third_derivatives.size() == max_q_indices + n_bubbles ||
277  third_derivatives.size() == 0,
278  ExcDimensionMismatch2(third_derivatives.size(),
279  max_q_indices + n_bubbles,
280  0));
281  Assert(fourth_derivatives.size() == max_q_indices + n_bubbles ||
282  fourth_derivatives.size() == 0,
283  ExcDimensionMismatch2(fourth_derivatives.size(),
284  max_q_indices + n_bubbles,
285  0));
286 
287  bool do_values = false, do_grads = false, do_grad_grads = false;
288  bool do_3rd_derivatives = false, do_4th_derivatives = false;
289  if (values.empty() == false)
290  {
291  values.resize(tensor_polys.n());
292  do_values = true;
293  }
294  if (grads.empty() == false)
295  {
296  grads.resize(tensor_polys.n());
297  do_grads = true;
298  }
299  if (grad_grads.empty() == false)
300  {
301  grad_grads.resize(tensor_polys.n());
302  do_grad_grads = true;
303  }
304  if (third_derivatives.empty() == false)
305  {
306  third_derivatives.resize(tensor_polys.n());
307  do_3rd_derivatives = true;
308  }
309  if (fourth_derivatives.empty() == false)
310  {
311  fourth_derivatives.resize(tensor_polys.n());
312  do_4th_derivatives = true;
313  }
314 
316  p, values, grads, grad_grads, third_derivatives, fourth_derivatives);
317 
318  for (unsigned int i = tensor_polys.n(); i < tensor_polys.n() + n_bubbles; ++i)
319  {
320  if (do_values)
321  values.push_back(compute_value(i, p));
322  if (do_grads)
323  grads.push_back(compute_grad(i, p));
324  if (do_grad_grads)
325  grad_grads.push_back(compute_grad_grad(i, p));
326  if (do_3rd_derivatives)
327  third_derivatives.push_back(compute_derivative<3>(i, p));
328  if (do_4th_derivatives)
329  fourth_derivatives.push_back(compute_derivative<4>(i, p));
330  }
331 }
332 
333 
334 
335 template <int dim>
336 std::unique_ptr<ScalarPolynomialsBase<dim>>
338 {
339  return std_cxx14::make_unique<TensorProductPolynomialsBubbles<dim>>(*this);
340 }
341 
342 
343 /* ------------------- explicit instantiations -------------- */
347 
TensorProductPolynomialsBubbles::compute_grad_grad
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
Definition: tensor_product_polynomials_bubbles.cc:153
TensorProductPolynomialsBubbles::set_numbering
void set_numbering(const std::vector< unsigned int > &renumber)
Definition: tensor_product_polynomials_bubbles.cc:48
tensor_product_polynomials_bubbles.h
TensorProductPolynomialsBubbles::index_map
std::vector< unsigned int > index_map
Definition: tensor_product_polynomials_bubbles.h:202
TensorProductPolynomialsBubbles::clone
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Definition: tensor_product_polynomials_bubbles.cc:337
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
TensorProductPolynomialsBubbles
Definition: tensor_product_polynomials.h:37
TensorProductPolynomials::compute_grad
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
Definition: tensor_product_polynomials.cc:164
TensorProductPolynomials::compute_grad_grad
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const
Definition: tensor_product_polynomials.cc:212
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
TensorProductPolynomialsBubbles::output_indices
void output_indices(std::ostream &out) const
Definition: tensor_product_polynomials_bubbles.cc:31
StandardExceptions::ExcDimensionMismatch2
static ::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
TensorProductPolynomials::polynomials
std::vector< PolynomialType > polynomials
Definition: tensor_product_polynomials.h:222
ScalarPolynomialsBase::n
unsigned int n() const
Definition: scalar_polynomials_base.h:164
TensorProductPolynomialsBubbles::tensor_polys
TensorProductPolynomials< dim > tensor_polys
Definition: tensor_product_polynomials_bubbles.h:197
TensorProductPolynomialsBubbles::compute_value
double compute_value(const unsigned int i, const Point< dim > &p) const
Definition: tensor_product_polynomials_bubbles.cc:68
Tensor< 1, dim >
TensorProductPolynomialsBubbles::index_map_inverse
std::vector< unsigned int > index_map_inverse
Definition: tensor_product_polynomials_bubbles.h:207
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
TensorProductPolynomials::compute_value
double compute_value(const unsigned int i, const Point< dim > &p) const
Definition: tensor_product_polynomials.cc:144
exceptions.h
TensorProductPolynomials::evaluate
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
Definition: tensor_product_polynomials.cc:268
value
static const bool value
Definition: dof_tools_constraints.cc:433
TensorProductPolynomialsBubbles::evaluate
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
Definition: tensor_product_polynomials_bubbles.cc:257
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point< dim >
TensorProductPolynomials::compute_index
void compute_index(const unsigned int i, std::array< unsigned int, dim > &indices) const
Definition: tensor_product_polynomials.cc:83
memory.h
TensorProductPolynomialsBubbles::compute_grad
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const
Definition: tensor_product_polynomials_bubbles.cc:107
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
TensorProductPolynomials::set_numbering
void set_numbering(const std::vector< unsigned int > &renumber)
Definition: tensor_product_polynomials.cc:117