Reference documentation for deal.II version GIT relicensing-769-g16d39143b0 2024-05-28 17:50:02+00:00
Searching...
No Matches
polynomials_barycentric.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
4// Copyright (C) 2021 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16
18
20
21namespace internal
22{
27 template <int dim>
28 unsigned int
30 const std::vector<typename BarycentricPolynomials<dim>::PolyType> &polys)
31 {
32 // Since the first variable in a simplex polynomial is, e.g., in 2d,
33 //
34 // t0 = 1 - x - y
35 //
36 // (that is, it depends on the Cartesian variables), we have to compute
37 // its degree separately. An example: t0*t1*t2 has degree 1 in the affine
38 // polynomial basis but is degree 2 in the Cartesian polynomial basis.
39 std::size_t max_degree = 0;
40 for (const auto &poly : polys)
41 {
42 const TableIndices<dim + 1> degrees = poly.degrees();
43
44 const auto degree_0 = degrees[0];
45 std::size_t degree_d = 0;
46 for (unsigned int d = 1; d < dim + 1; ++d)
47 degree_d = std::max(degree_d, degrees[d]);
48
49 max_degree = std::max(max_degree, degree_d + degree_0);
50 }
51
52 return max_degree;
53 }
54} // namespace internal
55
56
57template <int dim>
60{
61 std::vector<PolyType> polys;
62
63 const auto reference_cell = ReferenceCells::get_simplex<dim>();
64
65 switch (degree)
66 {
67 case 0:
69 1);
70 break;
71 case 1:
72 {
73 for (const unsigned int v : reference_cell.vertex_indices())
75 break;
76 }
77 case 2:
78 {
79 // vertices, then lines:
80 for (const unsigned int v : reference_cell.vertex_indices())
81 polys.push_back(
84 for (const unsigned int l : reference_cell.line_indices())
85 {
86 const auto v0 = reference_cell.line_to_cell_vertices(l, 0);
87 const auto v1 = reference_cell.line_to_cell_vertices(l, 1);
88 polys.push_back(4 *
91 }
92 break;
93 }
94 case 3:
95 {
96 // vertices, then lines, then quads:
97 for (const unsigned int v : reference_cell.vertex_indices())
98 polys.push_back(
102 for (unsigned int l : reference_cell.line_indices())
103 {
104 const auto v0 = reference_cell.line_to_cell_vertices(l, 0);
105 const auto v1 = reference_cell.line_to_cell_vertices(l, 1);
106 polys.push_back(
110 polys.push_back(
114 }
115
116 if (dim == 2)
117 {
118 polys.push_back(27 *
122 }
123 else if (dim == 3)
124 {
125 polys.push_back(27 *
129 polys.push_back(27 *
133 polys.push_back(27 *
137 polys.push_back(27 *
141 }
142
143 break;
144 }
145 default:
147 }
148
149 return BarycentricPolynomials<dim>(polys);
150}
151
152
153
154template <int dim>
156 const std::vector<PolyType> &polynomials)
157 : ScalarPolynomialsBase<dim>(internal::get_degree<dim>(polynomials),
158 polynomials.size())
159{
160 polys = polynomials;
161
163 poly_hessians.resize(polynomials.size());
164 poly_third_derivatives.resize(polynomials.size());
165 poly_fourth_derivatives.resize(polynomials.size());
166
167 for (std::size_t i = 0; i < polynomials.size(); ++i)
168 {
170 for (unsigned int d = 0; d < dim; ++d)
171 poly_grads[i][d] = polynomials[i].derivative(d);
172
173 // hessians
174 for (unsigned int d0 = 0; d0 < dim; ++d0)
175 for (unsigned int d1 = 0; d1 < dim; ++d1)
176 poly_hessians[i][d0][d1] = poly_grads[i][d0].derivative(d1);
177
178 // third derivatives
179 for (unsigned int d0 = 0; d0 < dim; ++d0)
180 for (unsigned int d1 = 0; d1 < dim; ++d1)
181 for (unsigned int d2 = 0; d2 < dim; ++d2)
182 poly_third_derivatives[i][d0][d1][d2] =
183 poly_hessians[i][d0][d1].derivative(d2);
184
185 // fourth derivatives
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 for (unsigned int d3 = 0; d3 < dim; ++d3)
190 poly_fourth_derivatives[i][d0][d1][d2][d3] =
191 poly_third_derivatives[i][d0][d1][d2].derivative(d3);
192 }
193}
194
195
196
197template <int dim>
198void
200 const Point<dim> &unit_point,
201 std::vector<double> &values,
202 std::vector<Tensor<1, dim>> &grads,
204 std::vector<Tensor<3, dim>> &third_derivatives,
205 std::vector<Tensor<4, dim>> &fourth_derivatives) const
206{
207 Assert(values.size() == this->n() || values.empty(),
208 ExcDimensionMismatch2(values.size(), this->n(), 0));
210 ExcDimensionMismatch2(grads.size(), this->n(), 0));
213 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
214 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
215 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
216 ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
217
218 for (std::size_t i = 0; i < polys.size(); ++i)
219 {
220 if (values.size() == this->n())
221 values[i] = polys[i].value(unit_point);
222
224 if (grads.size() == this->n())
225 for (unsigned int d = 0; d < dim; ++d)
227
228 // hessians
230 for (unsigned int d0 = 0; d0 < dim; ++d0)
231 for (unsigned int d1 = 0; d1 < dim; ++d1)
233
234 // third derivatives
235 if (third_derivatives.size() == this->n())
236 for (unsigned int d0 = 0; d0 < dim; ++d0)
237 for (unsigned int d1 = 0; d1 < dim; ++d1)
238 for (unsigned int d2 = 0; d2 < dim; ++d2)
239 third_derivatives[i][d0][d1][d2] =
240 poly_third_derivatives[i][d0][d1][d2].value(unit_point);
241
242 // fourth derivatives
243 if (fourth_derivatives.size() == this->n())
244 for (unsigned int d0 = 0; d0 < dim; ++d0)
245 for (unsigned int d1 = 0; d1 < dim; ++d1)
246 for (unsigned int d2 = 0; d2 < dim; ++d2)
247 for (unsigned int d3 = 0; d3 < dim; ++d3)
248 fourth_derivatives[i][d0][d1][d2][d3] =
249 poly_fourth_derivatives[i][d0][d1][d2][d3].value(unit_point);
250 }
251}
252
254
255template <int dim>
256double
258 const Point<dim> &p) const
259{
260 AssertIndexRange(i, this->n());
261 return polys[i].value(p);
262}
263
264
265
266template <int dim>
269 const Point<dim> &p) const
270{
272 for (unsigned int d = 0; d < dim; ++d)
273 result[d] = poly_grads[i][d].value(p);
274 return result;
275}
276
277
278
279template <int dim>
282 const Point<dim> &p) const
283{
284 Tensor<2, dim> result;
285 for (unsigned int d0 = 0; d0 < dim; ++d0)
286 for (unsigned int d1 = 0; d1 < dim; ++d1)
287 result[d0][d1] = poly_hessians[i][d0][d1].value(p);
289 return result;
290}
291
292
293
294template <int dim>
297 const Point<dim> &p) const
298{
299 Tensor<3, dim> result;
300 for (unsigned int d0 = 0; d0 < dim; ++d0)
301 for (unsigned int d1 = 0; d1 < dim; ++d1)
302 for (unsigned int d2 = 0; d2 < dim; ++d2)
303 result[d0][d1][d2] = poly_third_derivatives[i][d0][d1][d2].value(p);
304
305 return result;
306}
307
308
310template <int dim>
313 const Point<dim> &p) const
314{
315 Tensor<4, dim> result;
316 for (unsigned int d0 = 0; d0 < dim; ++d0)
317 for (unsigned int d1 = 0; d1 < dim; ++d1)
318 for (unsigned int d2 = 0; d2 < dim; ++d2)
319 for (unsigned int d3 = 0; d3 < dim; ++d3)
320 result[d0][d1][d2][d3] =
321 poly_fourth_derivatives[i][d0][d1][d2][d3].value(p);
323 return result;
324}
325
326
327
328template <int dim>
331 const Point<dim> &p) const
332{
333 return compute_1st_derivative(i, p);
335
336
337
338template <int dim>
341 const Point<dim> &p) const
342{
343 return compute_2nd_derivative(i, p);
344}
345
346
347
348template <int dim>
349std::unique_ptr<ScalarPolynomialsBase<dim>>
351{
352 return std::make_unique<BarycentricPolynomials<dim>>(*this);
353}
354
355
356
357template <int dim>
358std::string
360{
361 return "BarycentricPolynomials<" + std::to_string(dim) + ">";
362}
363
364
365
366template <int dim>
367std::size_t
369{
370 std::size_t poly_memory = 0;
371 for (const auto &poly : polys)
372 poly_memory += poly.memory_consumption();
376 MemoryConsumption::memory_consumption(poly_third_derivatives) +
377 MemoryConsumption::memory_consumption(poly_fourth_derivatives);
378}
379
380template class BarycentricPolynomials<1>;
381template class BarycentricPolynomials<2>;
382template class BarycentricPolynomials<3>;
383
static BarycentricPolynomial< dim, Number > monomial(const unsigned int d)
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
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
Tensor< 2, dim > compute_2nd_derivative(const unsigned int i, const Point< dim > &p) const override
std::vector< PolyType > polys
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
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:111
virtual std::size_t memory_consumption() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
const unsigned int v0
const unsigned int v1
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertIndexRange(index, range)
#define DEAL_II_NOT_IMPLEMENTED()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
unsigned int get_degree(const std::vector< typename BarycentricPolynomials< dim >::PolyType > &polys)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)