Reference documentation for deal.II version 9.3.0
polynomials_barycentric.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2021 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
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
18
20
21 namespace internal
22 {
27  template <int dim>
28  unsigned int
29  get_degree(const std::vector<BarycentricPolynomial<dim>> &polys)
30  {
31  // Since the first variable in a simplex polynomial is, e.g., in 2D,
32  //
33  // t0 = 1 - x - y
34  //
35  // (that is, it depends on the Cartesian variables), we have to compute
36  // its degree separately. An example: t0*t1*t2 has degree 1 in the affine
37  // polynomial basis but is degree 2 in the Cartesian polynomial basis.
38  std::size_t max_degree = 0;
39  for (const auto &poly : polys)
40  {
41  const TableIndices<dim + 1> degrees = poly.degrees();
42
43  const auto degree_0 = degrees[0];
44  std::size_t degree_d = 0;
45  for (unsigned int d = 1; d < dim + 1; ++d)
46  degree_d = std::max(degree_d, degrees[d]);
47
48  max_degree = std::max(max_degree, degree_d + degree_0);
49  }
50
51  return max_degree;
52  }
53 } // namespace internal
54
55
56 template <int dim>
59 {
60  std::vector<BarycentricPolynomial<dim>> polys;
61
62  auto M = [](const unsigned int d) {
64  };
65  switch (degree)
66  {
67  case 0:
68  polys.push_back(0 * M(0) + 1);
69  break;
70  case 1:
71  {
72  for (unsigned int d = 0; d < dim + 1; ++d)
73  polys.push_back(M(d));
74  break;
75  }
76  case 2:
77  {
78  for (unsigned int d = 0; d < dim + 1; ++d)
79  polys.push_back(M(d) * (2 * M(d) - 1));
80  polys.push_back(4 * M(1) * M(0));
81  if (dim >= 2)
82  {
83  polys.push_back(4 * M(1) * M(2));
84  polys.push_back(4 * M(2) * M(0));
85  }
86  if (dim == 3)
87  {
88  polys.push_back(4 * M(3) * M(0));
89  polys.push_back(4 * M(1) * M(3));
90  polys.push_back(4 * M(2) * M(3));
91  }
92  break;
93  }
94  default:
95  Assert(false, ExcNotImplemented());
96  }
97
98  return BarycentricPolynomials<dim>(polys);
99 }
100
101
102
103 template <int dim>
105  const std::vector<BarycentricPolynomial<dim>> &polynomials)
106  : ScalarPolynomialsBase<dim>(internal::get_degree(polynomials),
107  polynomials.size())
108 {
109  polys = polynomials;
110
112  poly_hessians.reinit({polynomials.size(), dim, dim});
113  poly_third_derivatives.reinit({polynomials.size(), dim, dim, dim});
114  poly_fourth_derivatives.reinit({polynomials.size(), dim, dim, dim, dim});
115
116  for (std::size_t i = 0; i < polynomials.size(); ++i)
117  {
119  for (unsigned int d = 0; d < dim; ++d)
121
122  // hessians
123  for (unsigned int d0 = 0; d0 < dim; ++d0)
124  for (unsigned int d1 = 0; d1 < dim; ++d1)
126
127  // third derivatives
128  for (unsigned int d0 = 0; d0 < dim; ++d0)
129  for (unsigned int d1 = 0; d1 < dim; ++d1)
130  for (unsigned int d2 = 0; d2 < dim; ++d2)
131  poly_third_derivatives[i][d0][d1][d2] =
132  poly_hessians[i][d0][d1].derivative(d2);
133
134  // fourth derivatives
135  for (unsigned int d0 = 0; d0 < dim; ++d0)
136  for (unsigned int d1 = 0; d1 < dim; ++d1)
137  for (unsigned int d2 = 0; d2 < dim; ++d2)
138  for (unsigned int d3 = 0; d3 < dim; ++d3)
139  poly_fourth_derivatives[i][d0][d1][d2][d3] =
140  poly_third_derivatives[i][d0][d1][d2].derivative(d3);
141  }
142 }
143
144
145
146 template <int dim>
147 void
149  const Point<dim> & unit_point,
150  std::vector<double> & values,
153  std::vector<Tensor<3, dim>> &third_derivatives,
154  std::vector<Tensor<4, dim>> &fourth_derivatives) const
155 {
156  Assert(values.size() == this->n() || values.size() == 0,
157  ExcDimensionMismatch2(values.size(), this->n(), 0));
162  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
163  ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
164  Assert(fourth_derivatives.size() == this->n() ||
165  fourth_derivatives.size() == 0,
166  ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
167
168  for (std::size_t i = 0; i < polys.size(); ++i)
169  {
170  if (values.size() == this->n())
171  values[i] = polys[i].value(unit_point);
172
175  for (unsigned int d = 0; d < dim; ++d)
177
178  // hessians
180  for (unsigned int d0 = 0; d0 < dim; ++d0)
181  for (unsigned int d1 = 0; d1 < dim; ++d1)
183
184  // third derivatives
185  if (third_derivatives.size() == this->n())
186  for (unsigned int d0 = 0; d0 < dim; ++d0)
187  for (unsigned int d1 = 0; d1 < dim; ++d1)
188  for (unsigned int d2 = 0; d2 < dim; ++d2)
189  third_derivatives[i][d0][d1][d2] =
190  poly_third_derivatives[i][d0][d1][d2].value(unit_point);
191
192  // fourth derivatives
193  if (fourth_derivatives.size() == this->n())
194  for (unsigned int d0 = 0; d0 < dim; ++d0)
195  for (unsigned int d1 = 0; d1 < dim; ++d1)
196  for (unsigned int d2 = 0; d2 < dim; ++d2)
197  for (unsigned int d3 = 0; d3 < dim; ++d3)
198  fourth_derivatives[i][d0][d1][d2][d3] =
199  poly_fourth_derivatives[i][d0][d1][d2][d3].value(unit_point);
200  }
201 }
202
203
204
205 template <int dim>
206 double
208  const Point<dim> & p) const
209 {
210  AssertIndexRange(i, this->n());
211  return polys[i].value(p);
212 }
213
214
215
216 template <int dim>
219  const Point<dim> & p) const
220 {
221  Tensor<1, dim> result;
222  for (unsigned int d = 0; d < dim; ++d)
224  return result;
225 }
226
227
228
229 template <int dim>
232  const Point<dim> & p) const
233 {
234  Tensor<2, dim> result;
235  for (unsigned int d0 = 0; d0 < dim; ++d0)
236  for (unsigned int d1 = 0; d1 < dim; ++d1)
237  result[d0][d1] = poly_hessians[i][d0][d1].value(p);
238
239  return result;
240 }
241
242
243
244 template <int dim>
247  const Point<dim> & p) const
248 {
249  Tensor<3, dim> result;
250  for (unsigned int d0 = 0; d0 < dim; ++d0)
251  for (unsigned int d1 = 0; d1 < dim; ++d1)
252  for (unsigned int d2 = 0; d2 < dim; ++d2)
253  result[d0][d1][d2] = poly_third_derivatives[i][d0][d1][d2].value(p);
254
255  return result;
256 }
257
258
259
260 template <int dim>
263  const Point<dim> & p) const
264 {
265  Tensor<4, dim> result;
266  for (unsigned int d0 = 0; d0 < dim; ++d0)
267  for (unsigned int d1 = 0; d1 < dim; ++d1)
268  for (unsigned int d2 = 0; d2 < dim; ++d2)
269  for (unsigned int d3 = 0; d3 < dim; ++d3)
270  result[d0][d1][d2][d3] =
271  poly_fourth_derivatives[i][d0][d1][d2][d3].value(p);
272
273  return result;
274 }
275
276
277
278 template <int dim>
281  const Point<dim> & p) const
282 {
283  return compute_1st_derivative(i, p);
284 }
285
286
287
288 template <int dim>
291  const Point<dim> & p) const
292 {
293  return compute_2nd_derivative(i, p);
294 }
295
296
297
298 template <int dim>
299 std::unique_ptr<ScalarPolynomialsBase<dim>>
301 {
302  return std::make_unique<BarycentricPolynomials<dim>>(*this);
303 }
304
305
306
307 template <int dim>
308 std::string
310 {
311  return "BarycentricPolynomials<" + std::to_string(dim) + ">";
312 }
313
314
315
316 template <int dim>
317 std::size_t
319 {
320  std::size_t poly_memory = 0;
321  for (const auto &poly : polys)
322  poly_memory += poly.memory_consumption();
323  return ScalarPolynomialsBase<dim>::memory_consumption() + poly_memory +
325  poly_third_derivatives.memory_consumption() +
326  poly_fourth_derivatives.memory_consumption();
327 }
328
329 template class BarycentricPolynomials<1>;
330 template class BarycentricPolynomials<2>;
331 template class BarycentricPolynomials<3>;
332
Table< 2, BarycentricPolynomial< dim > > poly_grads
Table< 5, BarycentricPolynomial< dim > > poly_fourth_derivatives
double compute_value(const unsigned int i, const Point< dim > &p) const override
static BarycentricPolynomial< dim, Number > monomial(const unsigned int d)
std::string name() const override
static ::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
Table< 4, BarycentricPolynomial< dim > > poly_third_derivatives
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
unsigned int get_degree(const std::vector< BarycentricPolynomial< dim >> &polys)
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
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
std::string to_string(const T &t)
Definition: patterns.h:2329
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) const override
Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
virtual std::size_t memory_consumption() const
static ::ExceptionBase & ExcNotImplemented()
BarycentricPolynomials(const std::vector< BarycentricPolynomial< dim >> &polynomials)
Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
Table< 3, BarycentricPolynomial< dim > > poly_hessians
std::vector< BarycentricPolynomial< dim > > polys
virtual std::size_t memory_consumption() const override