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