Reference documentation for deal.II version GIT 8ad845fa53 2023-12-10 20:50:02+00:00
\(\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.empty(),
216  ExcDimensionMismatch2(values.size(), this->n(), 0));
217  Assert(grads.size() == this->n() || grads.empty(),
218  ExcDimensionMismatch2(grads.size(), this->n(), 0));
219  Assert(grad_grads.size() == this->n() || grad_grads.empty(),
220  ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
221  Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
222  ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
223  Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
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 
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::make_unique<PolynomialSpace<dim>>(*this);
407 }
408 
409 
410 template class PolynomialSpace<1>;
411 template class PolynomialSpace<2>;
412 template class PolynomialSpace<3>;
413 
Definition: point.h:112
std::array< unsigned int, dim > compute_index(const unsigned int n) const
double compute_value(const unsigned int i, const Point< dim > &p) const override
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
void set_numbering(const std::vector< unsigned int > &renumber)
static unsigned int n_polynomials(const unsigned int n)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static const unsigned int invalid_unsigned_int
Definition: types.h:221