Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
22namespace 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
58template <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:
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 default:
96 Assert(false, ExcNotImplemented());
97 }
98
99 return BarycentricPolynomials<dim>(polys);
100}
101
102
103
104template <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
147template <int dim>
148void
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
206template <int dim>
207double
209 const Point<dim> & p) const
210{
211 AssertIndexRange(i, this->n());
212 return polys[i].value(p);
213}
214
215
216
217template <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
230template <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
245template <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
260
261template <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
279template <int dim>
282 const Point<dim> & p) const
284 return compute_1st_derivative(i, p);
285}
286
287
288
289template <int dim>
292 const Point<dim> & p) const
293{
294 return compute_2nd_derivative(i, p);
295}
297
298
299template <int dim>
300std::unique_ptr<ScalarPolynomialsBase<dim>>
302{
303 return std::make_unique<BarycentricPolynomials<dim>>(*this);
304}
305
306
307
308template <int dim>
309std::string
311{
312 return "BarycentricPolynomials<" + std::to_string(dim) + ">";
313}
314
315
316
317template <int dim>
318std::size_t
320{
321 std::size_t poly_memory = 0;
322 for (const auto &poly : polys)
323 poly_memory += poly.memory_consumption();
327 MemoryConsumption::memory_consumption(poly_third_derivatives) +
328 MemoryConsumption::memory_consumption(poly_fourth_derivatives);
330
331template class BarycentricPolynomials<1>;
332template class BarycentricPolynomials<2>;
333template class BarycentricPolynomials<3>;
334
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::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
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:112
virtual std::size_t memory_consumption() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
const unsigned int v0
const unsigned int v1
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertIndexRange(index, range)
std::enable_if_t< std::is_fundamental< T >::value, 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 > &)