deal.II version GIT relicensing-1809-g201751ccb9 2024-09-08 15:20:00+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 , polys(polynomials)
160{
162 poly_hessians.resize(polynomials.size());
163 poly_third_derivatives.resize(polynomials.size());
164 poly_fourth_derivatives.resize(polynomials.size());
165
166 for (std::size_t i = 0; i < polynomials.size(); ++i)
167 {
169 for (unsigned int d = 0; d < dim; ++d)
171
172 // hessians
173 for (unsigned int d0 = 0; d0 < dim; ++d0)
174 for (unsigned int d1 = 0; d1 < dim; ++d1)
176
177 // third derivatives
178 for (unsigned int d0 = 0; d0 < dim; ++d0)
179 for (unsigned int d1 = 0; d1 < dim; ++d1)
180 for (unsigned int d2 = 0; d2 < dim; ++d2)
181 poly_third_derivatives[i][d0][d1][d2] =
182 poly_hessians[i][d0][d1].derivative(d2);
183
184 // fourth derivatives
185 for (unsigned int d0 = 0; d0 < dim; ++d0)
186 for (unsigned int d1 = 0; d1 < dim; ++d1)
187 for (unsigned int d2 = 0; d2 < dim; ++d2)
188 for (unsigned int d3 = 0; d3 < dim; ++d3)
189 poly_fourth_derivatives[i][d0][d1][d2][d3] =
190 poly_third_derivatives[i][d0][d1][d2].derivative(d3);
191 }
192}
193
194
195
196template <int dim>
197void
199 const Point<dim> &unit_point,
200 std::vector<double> &values,
203 std::vector<Tensor<3, dim>> &third_derivatives,
204 std::vector<Tensor<4, dim>> &fourth_derivatives) const
205{
206 Assert(values.size() == this->n() || values.empty(),
207 ExcDimensionMismatch2(values.size(), this->n(), 0));
212 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
213 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
214 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
215 ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
216
217 for (std::size_t i = 0; i < polys.size(); ++i)
218 {
219 if (values.size() == this->n())
220 values[i] = polys[i].value(unit_point);
221
224 for (unsigned int d = 0; d < dim; ++d)
226
227 // hessians
229 for (unsigned int d0 = 0; d0 < dim; ++d0)
230 for (unsigned int d1 = 0; d1 < dim; ++d1)
232
233 // third derivatives
234 if (third_derivatives.size() == this->n())
235 for (unsigned int d0 = 0; d0 < dim; ++d0)
236 for (unsigned int d1 = 0; d1 < dim; ++d1)
237 for (unsigned int d2 = 0; d2 < dim; ++d2)
238 third_derivatives[i][d0][d1][d2] =
239 poly_third_derivatives[i][d0][d1][d2].value(unit_point);
240
241 // fourth derivatives
242 if (fourth_derivatives.size() == this->n())
243 for (unsigned int d0 = 0; d0 < dim; ++d0)
244 for (unsigned int d1 = 0; d1 < dim; ++d1)
245 for (unsigned int d2 = 0; d2 < dim; ++d2)
246 for (unsigned int d3 = 0; d3 < dim; ++d3)
247 fourth_derivatives[i][d0][d1][d2][d3] =
248 poly_fourth_derivatives[i][d0][d1][d2][d3].value(unit_point);
249 }
250}
251
252
254template <int dim>
255double
257 const Point<dim> &p) const
259 AssertIndexRange(i, this->n());
260 return polys[i].value(p);
261}
262
263
264
265template <int dim>
268 const Point<dim> &p) const
269{
270 Tensor<1, dim> result;
271 for (unsigned int d = 0; d < dim; ++d)
273 return result;
274}
275
276
277
278template <int dim>
281 const Point<dim> &p) const
283 Tensor<2, dim> result;
284 for (unsigned int d0 = 0; d0 < dim; ++d0)
285 for (unsigned int d1 = 0; d1 < dim; ++d1)
286 result[d0][d1] = poly_hessians[i][d0][d1].value(p);
287
288 return result;
289}
290
291
292
293template <int dim>
296 const Point<dim> &p) const
297{
298 Tensor<3, dim> result;
299 for (unsigned int d0 = 0; d0 < dim; ++d0)
300 for (unsigned int d1 = 0; d1 < dim; ++d1)
301 for (unsigned int d2 = 0; d2 < dim; ++d2)
302 result[d0][d1][d2] = poly_third_derivatives[i][d0][d1][d2].value(p);
303
304 return result;
305}
306
307
308
309template <int dim>
312 const Point<dim> &p) const
313{
314 Tensor<4, dim> result;
315 for (unsigned int d0 = 0; d0 < dim; ++d0)
316 for (unsigned int d1 = 0; d1 < dim; ++d1)
317 for (unsigned int d2 = 0; d2 < dim; ++d2)
318 for (unsigned int d3 = 0; d3 < dim; ++d3)
319 result[d0][d1][d2][d3] =
320 poly_fourth_derivatives[i][d0][d1][d2][d3].value(p);
321
322 return result;
323}
324
325
326
327template <int dim>
330 const Point<dim> &p) const
331{
332 return compute_1st_derivative(i, p);
333}
335
336
337template <int dim>
340 const Point<dim> &p) const
341{
342 return compute_2nd_derivative(i, p);
343}
344
345
346
347template <int dim>
348std::unique_ptr<ScalarPolynomialsBase<dim>>
350{
351 return std::make_unique<BarycentricPolynomials<dim>>(*this);
352}
353
354
355
356template <int dim>
357std::string
359{
360 return "BarycentricPolynomials<" + std::to_string(dim) + ">";
361}
362
363
364
365template <int dim>
366std::size_t
368{
369 std::size_t poly_memory = 0;
370 for (const auto &poly : polys)
371 poly_memory += poly.memory_consumption();
375 MemoryConsumption::memory_consumption(poly_third_derivatives) +
376 MemoryConsumption::memory_consumption(poly_fourth_derivatives);
377}
378
379template class BarycentricPolynomials<1>;
380template class BarycentricPolynomials<2>;
381template class BarycentricPolynomials<3>;
382
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
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:500
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:501
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 > &)