Reference documentation for deal.II version 9.4.0
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
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 auto M = [](const unsigned int d) {
65 };
66 switch (degree)
67 {
68 case 0:
69 polys.push_back(0 * M(0) + 1);
70 break;
71 case 1:
72 {
73 for (unsigned int d = 0; d < dim + 1; ++d)
74 polys.push_back(M(d));
75 break;
76 }
77 case 2:
78 {
79 for (unsigned int d = 0; d < dim + 1; ++d)
80 polys.push_back(M(d) * (2 * M(d) - 1));
81 polys.push_back(4 * M(1) * M(0));
82 if (dim >= 2)
83 {
84 polys.push_back(4 * M(1) * M(2));
85 polys.push_back(4 * M(2) * M(0));
86 }
87 if (dim == 3)
88 {
89 polys.push_back(4 * M(3) * M(0));
90 polys.push_back(4 * M(1) * M(3));
91 polys.push_back(4 * M(2) * M(3));
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
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 {
120 for (unsigned int d = 0; d < dim; ++d)
122
123 // hessians
124 for (unsigned int d0 = 0; d0 < dim; ++d0)
125 for (unsigned int d1 = 0; d1 < dim; ++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,
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));
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
176 for (unsigned int d = 0; d < dim; ++d)
178
179 // hessians
181 for (unsigned int d0 = 0; d0 < dim; ++d0)
182 for (unsigned int d1 = 0; d1 < dim; ++d1)
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)
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;
258
259
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
283{
284 return compute_1st_derivative(i, p);
285}
286
288
289template <int dim>
292 const Point<dim> & p) const
293{
295}
296
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
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();
328 MemoryConsumption::memory_consumption(poly_fourth_derivatives);
329}
330
331template class BarycentricPolynomials<1>;
332template class BarycentricPolynomials<2>;
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::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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::string to_string(const T &t)
Definition: patterns.h:2403
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
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< typename BarycentricPolynomials< dim >::PolyType > &polys)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)