Reference documentation for deal.II version GIT c1ddcf411d 2023-11-30 16:30: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\}}\)
polynomials_barycentric.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2023 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 
17 
19 
21 
22 namespace internal
23 {
28  template <int dim>
29  unsigned int
31  const std::vector<typename BarycentricPolynomials<dim>::PolyType> &polys)
32  {
33  // Since the first variable in a simplex polynomial is, e.g., in 2d,
34  //
35  // t0 = 1 - x - y
36  //
37  // (that is, it depends on the Cartesian variables), we have to compute
38  // its degree separately. An example: t0*t1*t2 has degree 1 in the affine
39  // polynomial basis but is degree 2 in the Cartesian polynomial basis.
40  std::size_t max_degree = 0;
41  for (const auto &poly : polys)
42  {
43  const TableIndices<dim + 1> degrees = poly.degrees();
44 
45  const auto degree_0 = degrees[0];
46  std::size_t degree_d = 0;
47  for (unsigned int d = 1; d < dim + 1; ++d)
48  degree_d = std::max(degree_d, degrees[d]);
49 
50  max_degree = std::max(max_degree, degree_d + degree_0);
51  }
52 
53  return max_degree;
54  }
55 } // namespace internal
56 
57 
58 template <int dim>
61 {
62  std::vector<PolyType> polys;
63 
64  const auto reference_cell = ReferenceCells::get_simplex<dim>();
65 
66  switch (degree)
67  {
68  case 0:
69  polys.push_back(0 * BarycentricPolynomial<dim, double>::monomial(0) +
70  1);
71  break;
72  case 1:
73  {
74  for (const unsigned int v : reference_cell.vertex_indices())
76  break;
77  }
78  case 2:
79  {
80  // vertices, then lines:
81  for (const unsigned int v : reference_cell.vertex_indices())
82  polys.push_back(
85  for (const unsigned int l : reference_cell.line_indices())
86  {
87  const auto v0 = reference_cell.line_to_cell_vertices(l, 0);
88  const auto v1 = reference_cell.line_to_cell_vertices(l, 1);
89  polys.push_back(4 *
92  }
93  break;
94  }
95  case 3:
96  {
97  // vertices, then lines, then quads:
98  for (const unsigned int v : reference_cell.vertex_indices())
99  polys.push_back(
103  for (unsigned int l : reference_cell.line_indices())
104  {
105  const auto v0 = reference_cell.line_to_cell_vertices(l, 0);
106  const auto v1 = reference_cell.line_to_cell_vertices(l, 1);
107  polys.push_back(
111  polys.push_back(
115  }
116 
117  if (dim == 2)
118  {
119  polys.push_back(27 *
123  }
124  else if (dim == 3)
125  {
126  polys.push_back(27 *
130  polys.push_back(27 *
134  polys.push_back(27 *
138  polys.push_back(27 *
142  }
143 
144  break;
145  }
146  default:
147  Assert(false, ExcNotImplemented());
148  }
149 
150  return BarycentricPolynomials<dim>(polys);
151 }
152 
153 
154 
155 template <int dim>
157  const std::vector<PolyType> &polynomials)
158  : ScalarPolynomialsBase<dim>(internal::get_degree<dim>(polynomials),
159  polynomials.size())
160 {
161  polys = polynomials;
162 
163  poly_grads.resize(polynomials.size());
164  poly_hessians.resize(polynomials.size());
165  poly_third_derivatives.resize(polynomials.size());
166  poly_fourth_derivatives.resize(polynomials.size());
167 
168  for (std::size_t i = 0; i < polynomials.size(); ++i)
169  {
170  // gradients
171  for (unsigned int d = 0; d < dim; ++d)
172  poly_grads[i][d] = polynomials[i].derivative(d);
173 
174  // hessians
175  for (unsigned int d0 = 0; d0 < dim; ++d0)
176  for (unsigned int d1 = 0; d1 < dim; ++d1)
177  poly_hessians[i][d0][d1] = poly_grads[i][d0].derivative(d1);
178 
179  // third derivatives
180  for (unsigned int d0 = 0; d0 < dim; ++d0)
181  for (unsigned int d1 = 0; d1 < dim; ++d1)
182  for (unsigned int d2 = 0; d2 < dim; ++d2)
183  poly_third_derivatives[i][d0][d1][d2] =
184  poly_hessians[i][d0][d1].derivative(d2);
185 
186  // fourth derivatives
187  for (unsigned int d0 = 0; d0 < dim; ++d0)
188  for (unsigned int d1 = 0; d1 < dim; ++d1)
189  for (unsigned int d2 = 0; d2 < dim; ++d2)
190  for (unsigned int d3 = 0; d3 < dim; ++d3)
191  poly_fourth_derivatives[i][d0][d1][d2][d3] =
192  poly_third_derivatives[i][d0][d1][d2].derivative(d3);
193  }
194 }
195 
196 
197 
198 template <int dim>
199 void
201  const Point<dim> &unit_point,
202  std::vector<double> &values,
203  std::vector<Tensor<1, dim>> &grads,
204  std::vector<Tensor<2, dim>> &grad_grads,
205  std::vector<Tensor<3, dim>> &third_derivatives,
206  std::vector<Tensor<4, dim>> &fourth_derivatives) const
207 {
208  Assert(values.size() == this->n() || values.empty(),
209  ExcDimensionMismatch2(values.size(), this->n(), 0));
210  Assert(grads.size() == this->n() || grads.empty(),
211  ExcDimensionMismatch2(grads.size(), this->n(), 0));
212  Assert(grad_grads.size() == this->n() || grad_grads.empty(),
213  ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
214  Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
215  ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
216  Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
217  ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
218 
219  for (std::size_t i = 0; i < polys.size(); ++i)
220  {
221  if (values.size() == this->n())
222  values[i] = polys[i].value(unit_point);
223 
224  // gradients
225  if (grads.size() == this->n())
226  for (unsigned int d = 0; d < dim; ++d)
227  grads[i][d] = poly_grads[i][d].value(unit_point);
228 
229  // hessians
230  if (grad_grads.size() == this->n())
231  for (unsigned int d0 = 0; d0 < dim; ++d0)
232  for (unsigned int d1 = 0; d1 < dim; ++d1)
233  grad_grads[i][d0][d1] = poly_hessians[i][d0][d1].value(unit_point);
234 
235  // third derivatives
236  if (third_derivatives.size() == this->n())
237  for (unsigned int d0 = 0; d0 < dim; ++d0)
238  for (unsigned int d1 = 0; d1 < dim; ++d1)
239  for (unsigned int d2 = 0; d2 < dim; ++d2)
240  third_derivatives[i][d0][d1][d2] =
241  poly_third_derivatives[i][d0][d1][d2].value(unit_point);
242 
243  // fourth derivatives
244  if (fourth_derivatives.size() == this->n())
245  for (unsigned int d0 = 0; d0 < dim; ++d0)
246  for (unsigned int d1 = 0; d1 < dim; ++d1)
247  for (unsigned int d2 = 0; d2 < dim; ++d2)
248  for (unsigned int d3 = 0; d3 < dim; ++d3)
249  fourth_derivatives[i][d0][d1][d2][d3] =
250  poly_fourth_derivatives[i][d0][d1][d2][d3].value(unit_point);
251  }
252 }
253 
254 
255 
256 template <int dim>
257 double
259  const Point<dim> &p) const
260 {
261  AssertIndexRange(i, this->n());
262  return polys[i].value(p);
263 }
264 
265 
266 
267 template <int dim>
270  const Point<dim> &p) const
271 {
273  for (unsigned int d = 0; d < dim; ++d)
274  result[d] = poly_grads[i][d].value(p);
275  return result;
276 }
277 
278 
279 
280 template <int dim>
283  const Point<dim> &p) const
284 {
285  Tensor<2, dim> result;
286  for (unsigned int d0 = 0; d0 < dim; ++d0)
287  for (unsigned int d1 = 0; d1 < dim; ++d1)
288  result[d0][d1] = poly_hessians[i][d0][d1].value(p);
289 
290  return result;
291 }
292 
293 
294 
295 template <int dim>
298  const Point<dim> &p) const
299 {
300  Tensor<3, dim> result;
301  for (unsigned int d0 = 0; d0 < dim; ++d0)
302  for (unsigned int d1 = 0; d1 < dim; ++d1)
303  for (unsigned int d2 = 0; d2 < dim; ++d2)
304  result[d0][d1][d2] = poly_third_derivatives[i][d0][d1][d2].value(p);
305 
306  return result;
307 }
308 
309 
310 
311 template <int dim>
314  const Point<dim> &p) const
315 {
316  Tensor<4, dim> result;
317  for (unsigned int d0 = 0; d0 < dim; ++d0)
318  for (unsigned int d1 = 0; d1 < dim; ++d1)
319  for (unsigned int d2 = 0; d2 < dim; ++d2)
320  for (unsigned int d3 = 0; d3 < dim; ++d3)
321  result[d0][d1][d2][d3] =
322  poly_fourth_derivatives[i][d0][d1][d2][d3].value(p);
323 
324  return result;
325 }
326 
327 
328 
329 template <int dim>
332  const Point<dim> &p) const
333 {
334  return compute_1st_derivative(i, p);
335 }
336 
337 
338 
339 template <int dim>
342  const Point<dim> &p) const
343 {
344  return compute_2nd_derivative(i, p);
345 }
346 
347 
348 
349 template <int dim>
350 std::unique_ptr<ScalarPolynomialsBase<dim>>
352 {
353  return std::make_unique<BarycentricPolynomials<dim>>(*this);
354 }
355 
356 
357 
358 template <int dim>
359 std::string
361 {
362  return "BarycentricPolynomials<" + std::to_string(dim) + ">";
363 }
364 
365 
366 
367 template <int dim>
368 std::size_t
370 {
371  std::size_t poly_memory = 0;
372  for (const auto &poly : polys)
373  poly_memory += poly.memory_consumption();
374  return ScalarPolynomialsBase<dim>::memory_consumption() + poly_memory +
377  MemoryConsumption::memory_consumption(poly_third_derivatives) +
378  MemoryConsumption::memory_consumption(poly_fourth_derivatives);
379 }
380 
381 template class BarycentricPolynomials<1>;
382 template class BarycentricPolynomials<2>;
383 template class BarycentricPolynomials<3>;
384 
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
virtual std::size_t memory_consumption() const override
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
std::vector< GradType > poly_grads
std::string name() const override
Tensor< 3, dim > compute_3rd_derivative(const unsigned int i, const Point< dim > &p) const override
Tensor< 1, dim > compute_1st_derivative(const unsigned int i, const Point< dim > &p) 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
Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
std::vector< PolyType > polys
BarycentricPolynomials(const std::vector< BarycentricPolynomial< dim >> &polynomials)
double compute_value(const unsigned int i, const Point< dim > &p) const override
std::vector< ThirdDerivativesType > poly_third_derivatives
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() const override
Tensor< 4, dim > compute_4th_derivative(const unsigned int i, const Point< dim > &p) const override
std::vector< HessianType > poly_hessians
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
std::vector< FourthDerivativesType > poly_fourth_derivatives
Definition: point.h:112
virtual std::size_t memory_consumption() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
const unsigned int v0
const unsigned int v1
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
std::string to_string(const T &t)
Definition: patterns.h:2391
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int get_degree(const std::vector< typename BarycentricPolynomials< dim >::PolyType > &polys)