Reference documentation for deal.II version GIT 5983d193e2 2023-05-27 16: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\}}\)
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
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 (unsigned int v : reference_cell.vertex_indices())
76  break;
77  }
78  case 2:
79  {
80  // vertices, then lines:
81  for (unsigned int v : reference_cell.vertex_indices())
82  polys.push_back(
85  for (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  default:
96  Assert(false, ExcNotImplemented());
97  }
98 
99  return BarycentricPolynomials<dim>(polys);
100 }
101 
102 
103 
104 template <int dim>
106  const std::vector<PolyType> &polynomials)
107  : ScalarPolynomialsBase<dim>(internal::get_degree<dim>(polynomials),
108  polynomials.size())
109 {
110  polys = polynomials;
111 
112  poly_grads.resize(polynomials.size());
113  poly_hessians.resize(polynomials.size());
114  poly_third_derivatives.resize(polynomials.size());
115  poly_fourth_derivatives.resize(polynomials.size());
116 
117  for (std::size_t i = 0; i < polynomials.size(); ++i)
118  {
119  // gradients
120  for (unsigned int d = 0; d < dim; ++d)
121  poly_grads[i][d] = polynomials[i].derivative(d);
122 
123  // hessians
124  for (unsigned int d0 = 0; d0 < dim; ++d0)
125  for (unsigned int d1 = 0; d1 < dim; ++d1)
126  poly_hessians[i][d0][d1] = poly_grads[i][d0].derivative(d1);
127 
128  // third derivatives
129  for (unsigned int d0 = 0; d0 < dim; ++d0)
130  for (unsigned int d1 = 0; d1 < dim; ++d1)
131  for (unsigned int d2 = 0; d2 < dim; ++d2)
132  poly_third_derivatives[i][d0][d1][d2] =
133  poly_hessians[i][d0][d1].derivative(d2);
134 
135  // fourth derivatives
136  for (unsigned int d0 = 0; d0 < dim; ++d0)
137  for (unsigned int d1 = 0; d1 < dim; ++d1)
138  for (unsigned int d2 = 0; d2 < dim; ++d2)
139  for (unsigned int d3 = 0; d3 < dim; ++d3)
140  poly_fourth_derivatives[i][d0][d1][d2][d3] =
141  poly_third_derivatives[i][d0][d1][d2].derivative(d3);
142  }
143 }
144 
145 
146 
147 template <int dim>
148 void
150  const Point<dim> & unit_point,
151  std::vector<double> & values,
152  std::vector<Tensor<1, dim>> &grads,
153  std::vector<Tensor<2, dim>> &grad_grads,
154  std::vector<Tensor<3, dim>> &third_derivatives,
155  std::vector<Tensor<4, dim>> &fourth_derivatives) const
156 {
157  Assert(values.size() == this->n() || values.size() == 0,
158  ExcDimensionMismatch2(values.size(), this->n(), 0));
159  Assert(grads.size() == this->n() || grads.size() == 0,
160  ExcDimensionMismatch2(grads.size(), this->n(), 0));
161  Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
162  ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
163  Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
164  ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
165  Assert(fourth_derivatives.size() == this->n() ||
166  fourth_derivatives.size() == 0,
167  ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
168 
169  for (std::size_t i = 0; i < polys.size(); ++i)
170  {
171  if (values.size() == this->n())
172  values[i] = polys[i].value(unit_point);
173 
174  // gradients
175  if (grads.size() == this->n())
176  for (unsigned int d = 0; d < dim; ++d)
177  grads[i][d] = poly_grads[i][d].value(unit_point);
178 
179  // hessians
180  if (grad_grads.size() == this->n())
181  for (unsigned int d0 = 0; d0 < dim; ++d0)
182  for (unsigned int d1 = 0; d1 < dim; ++d1)
183  grad_grads[i][d0][d1] = poly_hessians[i][d0][d1].value(unit_point);
184 
185  // third derivatives
186  if (third_derivatives.size() == this->n())
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  third_derivatives[i][d0][d1][d2] =
191  poly_third_derivatives[i][d0][d1][d2].value(unit_point);
192 
193  // fourth derivatives
194  if (fourth_derivatives.size() == this->n())
195  for (unsigned int d0 = 0; d0 < dim; ++d0)
196  for (unsigned int d1 = 0; d1 < dim; ++d1)
197  for (unsigned int d2 = 0; d2 < dim; ++d2)
198  for (unsigned int d3 = 0; d3 < dim; ++d3)
199  fourth_derivatives[i][d0][d1][d2][d3] =
200  poly_fourth_derivatives[i][d0][d1][d2][d3].value(unit_point);
201  }
202 }
203 
204 
205 
206 template <int dim>
207 double
209  const Point<dim> & p) const
210 {
211  AssertIndexRange(i, this->n());
212  return polys[i].value(p);
213 }
214 
215 
216 
217 template <int dim>
220  const Point<dim> & p) const
221 {
222  Tensor<1, dim> result;
223  for (unsigned int d = 0; d < dim; ++d)
224  result[d] = poly_grads[i][d].value(p);
225  return result;
226 }
227 
228 
229 
230 template <int dim>
233  const Point<dim> & p) const
234 {
235  Tensor<2, dim> result;
236  for (unsigned int d0 = 0; d0 < dim; ++d0)
237  for (unsigned int d1 = 0; d1 < dim; ++d1)
238  result[d0][d1] = poly_hessians[i][d0][d1].value(p);
239 
240  return result;
241 }
242 
243 
244 
245 template <int dim>
248  const Point<dim> & p) const
249 {
250  Tensor<3, dim> result;
251  for (unsigned int d0 = 0; d0 < dim; ++d0)
252  for (unsigned int d1 = 0; d1 < dim; ++d1)
253  for (unsigned int d2 = 0; d2 < dim; ++d2)
254  result[d0][d1][d2] = poly_third_derivatives[i][d0][d1][d2].value(p);
255 
256  return result;
257 }
258 
259 
260 
261 template <int dim>
264  const Point<dim> & p) const
265 {
266  Tensor<4, dim> result;
267  for (unsigned int d0 = 0; d0 < dim; ++d0)
268  for (unsigned int d1 = 0; d1 < dim; ++d1)
269  for (unsigned int d2 = 0; d2 < dim; ++d2)
270  for (unsigned int d3 = 0; d3 < dim; ++d3)
271  result[d0][d1][d2][d3] =
272  poly_fourth_derivatives[i][d0][d1][d2][d3].value(p);
273 
274  return result;
275 }
276 
277 
278 
279 template <int dim>
282  const Point<dim> & p) const
283 {
284  return compute_1st_derivative(i, p);
285 }
286 
287 
288 
289 template <int dim>
292  const Point<dim> & p) const
293 {
294  return compute_2nd_derivative(i, p);
295 }
296 
297 
298 
299 template <int dim>
300 std::unique_ptr<ScalarPolynomialsBase<dim>>
302 {
303  return std::make_unique<BarycentricPolynomials<dim>>(*this);
304 }
305 
306 
307 
308 template <int dim>
309 std::string
311 {
312  return "BarycentricPolynomials<" + std::to_string(dim) + ">";
313 }
314 
315 
316 
317 template <int dim>
318 std::size_t
320 {
321  std::size_t poly_memory = 0;
322  for (const auto &poly : polys)
323  poly_memory += poly.memory_consumption();
324  return ScalarPolynomialsBase<dim>::memory_consumption() + poly_memory +
327  MemoryConsumption::memory_consumption(poly_third_derivatives) +
328  MemoryConsumption::memory_consumption(poly_fourth_derivatives);
329 }
330 
331 template class BarycentricPolynomials<1>;
332 template class BarycentricPolynomials<2>;
333 template class BarycentricPolynomials<3>;
334 
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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
const unsigned int v0
Definition: grid_tools.cc:1062
const unsigned int v1
Definition: grid_tools.cc:1062
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
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< T >::value, std::size_t > memory_consumption(const T &t)
std::string to_string(const T &t)
Definition: patterns.h:2392
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)