Reference documentation for deal.II version Git fe553f7db3 2020-12-03 08:11:48 +0100
\(\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\}}\)
polynomial_space.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2002 - 2020 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 
18 #include <deal.II/base/table.h>
19 
20 #include <memory>
21 
23 
24 
25 template <int dim>
26 unsigned int
28 {
29  unsigned int n_pols = n;
30  for (unsigned int i = 1; i < dim; ++i)
31  {
32  n_pols *= (n + i);
33  n_pols /= (i + 1);
34  }
35  return n_pols;
36 }
37 
38 
39 template <>
40 unsigned int
42 {
43  return 0;
44 }
45 
46 
47 template <>
48 std::array<unsigned int, 1>
49 PolynomialSpace<1>::compute_index(const unsigned int i) const
50 {
51  AssertIndexRange(i, index_map.size());
52  return {{index_map[i]}};
53 }
54 
55 
56 
57 template <>
58 std::array<unsigned int, 2>
59 PolynomialSpace<2>::compute_index(const unsigned int i) const
60 {
61  AssertIndexRange(i, index_map.size());
62  const unsigned int n = index_map[i];
63  // there should be a better way to
64  // write this function (not
65  // linear in n_1d), someone
66  // should think about this...
67  const unsigned int n_1d = polynomials.size();
68  unsigned int k = 0;
69  for (unsigned int iy = 0; iy < n_1d; ++iy)
70  if (n < k + n_1d - iy)
71  {
72  return {{n - k, iy}};
73  }
74  else
75  k += n_1d - iy;
76 
77  Assert(false, ExcInternalError());
79 }
80 
81 
82 
83 template <>
84 std::array<unsigned int, 3>
85 PolynomialSpace<3>::compute_index(const unsigned int i) const
86 {
87  AssertIndexRange(i, index_map.size());
88  const unsigned int n = index_map[i];
89  // there should be a better way to
90  // write this function (not
91  // quadratic in n_1d), someone
92  // should think about this...
93  //
94  // (ah, and yes: the original
95  // algorithm was even cubic!)
96  const unsigned int n_1d = polynomials.size();
97  unsigned int k = 0;
98  for (unsigned int iz = 0; iz < n_1d; ++iz)
99  for (unsigned int iy = 0; iy < n_1d - iz; ++iy)
100  if (n < k + n_1d - iy - iz)
101  {
102  return {{n - k, iy, iz}};
103  }
104  else
105  k += n_1d - iy - iz;
106 
107  Assert(false, ExcInternalError());
109 }
110 
111 
112 template <int dim>
113 void
114 PolynomialSpace<dim>::set_numbering(const std::vector<unsigned int> &renumber)
115 {
116  Assert(renumber.size() == index_map.size(),
117  ExcDimensionMismatch(renumber.size(), index_map.size()));
118 
119  index_map = renumber;
120  for (unsigned int i = 0; i < index_map.size(); ++i)
121  index_map_inverse[index_map[i]] = i;
122 }
123 
124 
125 
126 template <int dim>
127 double
129  const Point<dim> & p) const
130 {
131  const auto ix = compute_index(i);
132  // take the product of the
133  // polynomials in the various space
134  // directions
135  double result = 1.;
136  for (unsigned int d = 0; d < dim; ++d)
137  result *= polynomials[ix[d]].value(p(d));
138  return result;
139 }
140 
141 
142 
143 template <int dim>
146  const Point<dim> & p) const
147 {
148  const auto ix = compute_index(i);
149 
150  Tensor<1, dim> result;
151  for (unsigned int d = 0; d < dim; ++d)
152  result[d] = 1.;
153 
154  // get value and first derivative
155  std::vector<double> v(2);
156  for (unsigned int d = 0; d < dim; ++d)
157  {
158  polynomials[ix[d]].value(p(d), v);
159  result[d] *= v[1];
160  for (unsigned int d1 = 0; d1 < dim; ++d1)
161  if (d1 != d)
162  result[d1] *= v[0];
163  }
164  return result;
165 }
166 
167 
168 template <int dim>
171  const Point<dim> & p) const
172 {
173  const auto ix = compute_index(i);
174 
175  Tensor<2, dim> result;
176  for (unsigned int d = 0; d < dim; ++d)
177  for (unsigned int d1 = 0; d1 < dim; ++d1)
178  result[d][d1] = 1.;
179 
180  // get value, first and second
181  // derivatives
182  std::vector<double> v(3);
183  for (unsigned int d = 0; d < dim; ++d)
184  {
185  polynomials[ix[d]].value(p(d), v);
186  result[d][d] *= v[2];
187  for (unsigned int d1 = 0; d1 < dim; ++d1)
188  {
189  if (d1 != d)
190  {
191  result[d][d1] *= v[1];
192  result[d1][d] *= v[1];
193  for (unsigned int d2 = 0; d2 < dim; ++d2)
194  if (d2 != d)
195  result[d1][d2] *= v[0];
196  }
197  }
198  }
199  return result;
200 }
201 
202 
203 template <int dim>
204 void
206  const Point<dim> & p,
207  std::vector<double> & values,
208  std::vector<Tensor<1, dim>> &grads,
209  std::vector<Tensor<2, dim>> &grad_grads,
210  std::vector<Tensor<3, dim>> &third_derivatives,
211  std::vector<Tensor<4, dim>> &fourth_derivatives) const
212 {
213  const unsigned int n_1d = polynomials.size();
214 
215  Assert(values.size() == this->n() || values.size() == 0,
216  ExcDimensionMismatch2(values.size(), this->n(), 0));
217  Assert(grads.size() == this->n() || grads.size() == 0,
218  ExcDimensionMismatch2(grads.size(), this->n(), 0));
219  Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
220  ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
221  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
222  ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
223  Assert(fourth_derivatives.size() == this->n() ||
224  fourth_derivatives.size() == 0,
225  ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
226 
227  unsigned int v_size = 0;
228  bool update_values = false, update_grads = false, update_grad_grads = false;
229  bool update_3rd_derivatives = false, update_4th_derivatives = false;
230  if (values.size() == this->n())
231  {
232  update_values = true;
233  v_size = 1;
234  }
235  if (grads.size() == this->n())
236  {
237  update_grads = true;
238  v_size = 2;
239  }
240  if (grad_grads.size() == this->n())
241  {
242  update_grad_grads = true;
243  v_size = 3;
244  }
245  if (third_derivatives.size() == this->n())
246  {
247  update_3rd_derivatives = true;
248  v_size = 4;
249  }
250  if (fourth_derivatives.size() == this->n())
251  {
252  update_4th_derivatives = true;
253  v_size = 5;
254  }
255 
256  // Store data in a single
257  // object. Access is by
258  // v[d][n][o]
259  // d: coordinate direction
260  // n: number of 1d polynomial
261  // o: order of derivative
262  Table<2, std::vector<double>> v(dim, n_1d);
263  for (unsigned int d = 0; d < v.size()[0]; ++d)
264  for (unsigned int i = 0; i < v.size()[1]; ++i)
265  {
266  v(d, i).resize(v_size, 0.);
267  polynomials[i].value(p(d), v(d, i));
268  }
269 
270  if (update_values)
271  {
272  unsigned int k = 0;
273 
274  for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
275  for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
276  for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
277  values[index_map_inverse[k++]] = v[0][ix][0] *
278  ((dim > 1) ? v[1][iy][0] : 1.) *
279  ((dim > 2) ? v[2][iz][0] : 1.);
280  }
281 
282  if (update_grads)
283  {
284  unsigned int k = 0;
285 
286  for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
287  for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
288  for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
289  {
290  const unsigned int k2 = index_map_inverse[k++];
291  for (unsigned int d = 0; d < dim; ++d)
292  grads[k2][d] = v[0][ix][(d == 0) ? 1 : 0] *
293  ((dim > 1) ? v[1][iy][(d == 1) ? 1 : 0] : 1.) *
294  ((dim > 2) ? v[2][iz][(d == 2) ? 1 : 0] : 1.);
295  }
296  }
297 
298  if (update_grad_grads)
299  {
300  unsigned int k = 0;
301 
302  for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
303  for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
304  for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
305  {
306  const unsigned int k2 = index_map_inverse[k++];
307  for (unsigned int d1 = 0; d1 < dim; ++d1)
308  for (unsigned int d2 = 0; d2 < dim; ++d2)
309  {
310  // Derivative
311  // order for each
312  // direction
313  const unsigned int j0 =
314  ((d1 == 0) ? 1 : 0) + ((d2 == 0) ? 1 : 0);
315  const unsigned int j1 =
316  ((d1 == 1) ? 1 : 0) + ((d2 == 1) ? 1 : 0);
317  const unsigned int j2 =
318  ((d1 == 2) ? 1 : 0) + ((d2 == 2) ? 1 : 0);
319 
320  grad_grads[k2][d1][d2] = v[0][ix][j0] *
321  ((dim > 1) ? v[1][iy][j1] : 1.) *
322  ((dim > 2) ? v[2][iz][j2] : 1.);
323  }
324  }
325  }
326 
327  if (update_3rd_derivatives)
328  {
329  unsigned int k = 0;
330 
331  for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
332  for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
333  for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
334  {
335  const unsigned int k2 = index_map_inverse[k++];
336  for (unsigned int d1 = 0; d1 < dim; ++d1)
337  for (unsigned int d2 = 0; d2 < dim; ++d2)
338  for (unsigned int d3 = 0; d3 < dim; ++d3)
339  {
340  // Derivative
341  // order for each
342  // direction
343  std::vector<unsigned int> deriv_order(dim, 0);
344  for (unsigned int x = 0; x < dim; ++x)
345  {
346  if (d1 == x)
347  ++deriv_order[x];
348  if (d2 == x)
349  ++deriv_order[x];
350  if (d3 == x)
351  ++deriv_order[x];
352  }
353 
354  third_derivatives[k2][d1][d2][d3] =
355  v[0][ix][deriv_order[0]] *
356  ((dim > 1) ? v[1][iy][deriv_order[1]] : 1.) *
357  ((dim > 2) ? v[2][iz][deriv_order[2]] : 1.);
358  }
359  }
360  }
361 
362  if (update_4th_derivatives)
363  {
364  unsigned int k = 0;
365 
366  for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
367  for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
368  for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
369  {
370  const unsigned int k2 = index_map_inverse[k++];
371  for (unsigned int d1 = 0; d1 < dim; ++d1)
372  for (unsigned int d2 = 0; d2 < dim; ++d2)
373  for (unsigned int d3 = 0; d3 < dim; ++d3)
374  for (unsigned int d4 = 0; d4 < dim; ++d4)
375  {
376  // Derivative
377  // order for each
378  // direction
379  std::vector<unsigned int> deriv_order(dim, 0);
380  for (unsigned int x = 0; x < dim; ++x)
381  {
382  if (d1 == x)
383  ++deriv_order[x];
384  if (d2 == x)
385  ++deriv_order[x];
386  if (d3 == x)
387  ++deriv_order[x];
388  if (d4 == x)
389  ++deriv_order[x];
390  }
391 
392  fourth_derivatives[k2][d1][d2][d3][d4] =
393  v[0][ix][deriv_order[0]] *
394  ((dim > 1) ? v[1][iy][deriv_order[1]] : 1.) *
395  ((dim > 2) ? v[2][iz][deriv_order[2]] : 1.);
396  }
397  }
398  }
399 }
400 
401 
402 
403 template <int dim>
404 std::unique_ptr<ScalarPolynomialsBase<dim>>
406 {
407  return std::make_unique<PolynomialSpace<dim>>(*this);
408 }
409 
410 
411 template class PolynomialSpace<1>;
412 template class PolynomialSpace<2>;
413 template class PolynomialSpace<3>;
414 
Shape function values.
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
static const unsigned int invalid_unsigned_int
Definition: types.h:196
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
static ::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
void set_numbering(const std::vector< unsigned int > &renumber)
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Third derivatives of shape functions.
#define Assert(cond, exc)
Definition: exceptions.h:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
static unsigned int n_polynomials(const unsigned int n)
std::array< unsigned int, dim > compute_index(const unsigned int n) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
size_type size(const unsigned int i) const
double compute_value(const unsigned int i, const Point< dim > &p) const override
unsigned int compute_index()
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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
Definition: table.h:687
static ::ExceptionBase & ExcInternalError()