Reference documentation for deal.II version 9.3.3
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
21namespace 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
56template <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
103template <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
146template <int dim>
147void
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
205template <int dim>
206double
208 const Point<dim> & p) const
209{
210 AssertIndexRange(i, this->n());
211 return polys[i].value(p);
212}
213
214
215
216template <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}
227
228
229template <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
244template <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);
255 return result;
256}
257
258
259
260template <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;
275
276
277
278template <int dim>
281 const Point<dim> & p) const
282{
283 return compute_1st_derivative(i, p);
284}
285
286
287
288template <int dim>
291 const Point<dim> & p) const
292{
293 return compute_2nd_derivative(i, p);
295
296
297
298template <int dim>
299std::unique_ptr<ScalarPolynomialsBase<dim>>
301{
302 return std::make_unique<BarycentricPolynomials<dim>>(*this);
303}
304
305
307template <int dim>
308std::string
310{
311 return "BarycentricPolynomials<" + std::to_string(dim) + ">";
313
314
315
316template <int dim>
317std::size_t
319{
320 std::size_t poly_memory = 0;
321 for (const auto &poly : polys)
322 poly_memory += poly.memory_consumption();
325 poly_third_derivatives.memory_consumption() +
326 poly_fourth_derivatives.memory_consumption();
327}
328
329template class BarycentricPolynomials<1>;
330template class BarycentricPolynomials<2>;
331template class BarycentricPolynomials<3>;
332
static BarycentricPolynomial< dim, Number > monomial(const unsigned int d)
Table< 2, BarycentricPolynomial< dim > > poly_grads
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
BarycentricPolynomials(const std::vector< BarycentricPolynomial< dim > > &polynomials)
virtual std::size_t memory_consumption() const override
std::vector< BarycentricPolynomial< dim > > polys
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
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
Table< 3, BarycentricPolynomial< dim > > poly_hessians
Tensor< 2, dim > compute_2nd_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
double compute_value(const unsigned int i, const Point< dim > &p) const override
Table< 4, BarycentricPolynomial< dim > > poly_third_derivatives
Table< 5, BarycentricPolynomial< dim > > poly_fourth_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
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Definition: point.h:111
virtual std::size_t memory_consumption() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
std::string to_string(const T &t)
Definition: patterns.h:2329
static ::ExceptionBase & ExcDimensionMismatch2(int arg1, int arg2, int arg3)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int get_degree(const std::vector< BarycentricPolynomial< dim > > &polys)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)