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