Reference documentation for deal.II version GIT relicensing-453-g9b4d7b6a5a 2024-04-20 16:00:02+00:00
\(\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
polynomial_space.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2002 - 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
17#include <deal.II/base/table.h>
18
19#include <memory>
20
22
23
24template <int dim>
25unsigned int
27{
28 unsigned int n_pols = n;
29 for (unsigned int i = 1; i < dim; ++i)
30 {
31 n_pols *= (n + i);
32 n_pols /= (i + 1);
33 }
34 return n_pols;
35}
36
37
38template <>
39unsigned int
41{
42 return 0;
43}
44
45
46template <>
47std::array<unsigned int, 1>
48PolynomialSpace<1>::compute_index(const unsigned int i) const
49{
50 AssertIndexRange(i, index_map.size());
51 return {{index_map[i]}};
52}
53
54
55
56template <>
57std::array<unsigned int, 2>
58PolynomialSpace<2>::compute_index(const unsigned int i) const
59{
60 AssertIndexRange(i, index_map.size());
61 const unsigned int n = index_map[i];
62 // there should be a better way to
63 // write this function (not
64 // linear in n_1d), someone
65 // should think about this...
66 const unsigned int n_1d = polynomials.size();
67 unsigned int k = 0;
68 for (unsigned int iy = 0; iy < n_1d; ++iy)
69 if (n < k + n_1d - iy)
70 {
71 return {{n - k, iy}};
72 }
73 else
74 k += n_1d - iy;
75
78}
79
80
81
82template <>
83std::array<unsigned int, 3>
84PolynomialSpace<3>::compute_index(const unsigned int i) const
85{
86 AssertIndexRange(i, index_map.size());
87 const unsigned int n = index_map[i];
88 // there should be a better way to
89 // write this function (not
90 // quadratic in n_1d), someone
91 // should think about this...
92 //
93 // (ah, and yes: the original
94 // algorithm was even cubic!)
95 const unsigned int n_1d = polynomials.size();
96 unsigned int k = 0;
97 for (unsigned int iz = 0; iz < n_1d; ++iz)
98 for (unsigned int iy = 0; iy < n_1d - iz; ++iy)
99 if (n < k + n_1d - iy - iz)
100 {
101 return {{n - k, iy, iz}};
102 }
103 else
104 k += n_1d - iy - iz;
105
108}
109
110
111template <int dim>
112void
113PolynomialSpace<dim>::set_numbering(const std::vector<unsigned int> &renumber)
114{
115 Assert(renumber.size() == index_map.size(),
116 ExcDimensionMismatch(renumber.size(), index_map.size()));
117
118 index_map = renumber;
119 for (unsigned int i = 0; i < index_map.size(); ++i)
120 index_map_inverse[index_map[i]] = i;
121}
122
123
124
125template <int dim>
126double
128 const Point<dim> &p) const
129{
130 const auto ix = compute_index(i);
131 // take the product of the
132 // polynomials in the various space
133 // directions
134 double result = 1.;
135 for (unsigned int d = 0; d < dim; ++d)
136 result *= polynomials[ix[d]].value(p[d]);
137 return result;
138}
139
140
141
142template <int dim>
145 const Point<dim> &p) const
146{
147 const auto ix = compute_index(i);
148
149 Tensor<1, dim> result;
150 for (unsigned int d = 0; d < dim; ++d)
151 result[d] = 1.;
152
153 // get value and first derivative
154 std::vector<double> v(2);
155 for (unsigned int d = 0; d < dim; ++d)
156 {
157 polynomials[ix[d]].value(p[d], v);
158 result[d] *= v[1];
159 for (unsigned int d1 = 0; d1 < dim; ++d1)
160 if (d1 != d)
161 result[d1] *= v[0];
162 }
163 return result;
164}
165
166
167template <int dim>
170 const Point<dim> &p) const
171{
172 const auto ix = compute_index(i);
173
174 Tensor<2, dim> result;
175 for (unsigned int d = 0; d < dim; ++d)
176 for (unsigned int d1 = 0; d1 < dim; ++d1)
177 result[d][d1] = 1.;
178
179 // get value, first and second
180 // derivatives
181 std::vector<double> v(3);
182 for (unsigned int d = 0; d < dim; ++d)
183 {
184 polynomials[ix[d]].value(p[d], v);
185 result[d][d] *= v[2];
186 for (unsigned int d1 = 0; d1 < dim; ++d1)
187 {
188 if (d1 != d)
189 {
190 result[d][d1] *= v[1];
191 result[d1][d] *= v[1];
192 for (unsigned int d2 = 0; d2 < dim; ++d2)
193 if (d2 != d)
194 result[d1][d2] *= v[0];
195 }
196 }
197 }
198 return result;
199}
200
201
202template <int dim>
203void
205 const Point<dim> &p,
206 std::vector<double> &values,
207 std::vector<Tensor<1, dim>> &grads,
208 std::vector<Tensor<2, dim>> &grad_grads,
209 std::vector<Tensor<3, dim>> &third_derivatives,
210 std::vector<Tensor<4, dim>> &fourth_derivatives) const
211{
212 const unsigned int n_1d = polynomials.size();
213
214 Assert(values.size() == this->n() || values.empty(),
215 ExcDimensionMismatch2(values.size(), this->n(), 0));
216 Assert(grads.size() == this->n() || grads.empty(),
217 ExcDimensionMismatch2(grads.size(), this->n(), 0));
218 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
219 ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
220 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
221 ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
222 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
223 ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
224
225 unsigned int v_size = 0;
226 bool update_values = false, update_grads = false, update_grad_grads = false;
227 bool update_3rd_derivatives = false, update_4th_derivatives = false;
228 if (values.size() == this->n())
229 {
230 update_values = true;
231 v_size = 1;
232 }
233 if (grads.size() == this->n())
234 {
235 update_grads = true;
236 v_size = 2;
238 if (grad_grads.size() == this->n())
239 {
240 update_grad_grads = true;
241 v_size = 3;
242 }
243 if (third_derivatives.size() == this->n())
244 {
246 v_size = 4;
247 }
248 if (fourth_derivatives.size() == this->n())
249 {
250 update_4th_derivatives = true;
251 v_size = 5;
252 }
253
254 // Store data in a single
255 // object. Access is by
256 // v[d][n][o]
257 // d: coordinate direction
258 // n: number of 1d polynomial
259 // o: order of derivative
260 Table<2, std::vector<double>> v(dim, n_1d);
261 for (unsigned int d = 0; d < v.size()[0]; ++d)
262 for (unsigned int i = 0; i < v.size()[1]; ++i)
263 {
264 v(d, i).resize(v_size, 0.);
265 polynomials[i].value(p[d], v(d, i));
266 }
267
268 if (update_values)
269 {
270 unsigned int k = 0;
271
272 for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
273 for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
274 for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
275 values[index_map_inverse[k++]] = v[0][ix][0] *
276 ((dim > 1) ? v[1][iy][0] : 1.) *
277 ((dim > 2) ? v[2][iz][0] : 1.);
278 }
279
280 if (update_grads)
281 {
282 unsigned int k = 0;
283
284 for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
285 for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
286 for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
287 {
288 const unsigned int k2 = index_map_inverse[k++];
289 for (unsigned int d = 0; d < dim; ++d)
290 grads[k2][d] = v[0][ix][(d == 0) ? 1 : 0] *
291 ((dim > 1) ? v[1][iy][(d == 1) ? 1 : 0] : 1.) *
292 ((dim > 2) ? v[2][iz][(d == 2) ? 1 : 0] : 1.);
293 }
294 }
295
296 if (update_grad_grads)
297 {
298 unsigned int k = 0;
299
300 for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
301 for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
302 for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
303 {
304 const unsigned int k2 = index_map_inverse[k++];
305 for (unsigned int d1 = 0; d1 < dim; ++d1)
306 for (unsigned int d2 = 0; d2 < dim; ++d2)
307 {
308 // Derivative
309 // order for each
310 // direction
311 const unsigned int j0 =
312 ((d1 == 0) ? 1 : 0) + ((d2 == 0) ? 1 : 0);
313 const unsigned int j1 =
314 ((d1 == 1) ? 1 : 0) + ((d2 == 1) ? 1 : 0);
315 const unsigned int j2 =
316 ((d1 == 2) ? 1 : 0) + ((d2 == 2) ? 1 : 0);
317
318 grad_grads[k2][d1][d2] = v[0][ix][j0] *
319 ((dim > 1) ? v[1][iy][j1] : 1.) *
320 ((dim > 2) ? v[2][iz][j2] : 1.);
321 }
322 }
323 }
324
326 {
327 unsigned int k = 0;
328
329 for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
330 for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
331 for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
332 {
333 const unsigned int k2 = index_map_inverse[k++];
334 for (unsigned int d1 = 0; d1 < dim; ++d1)
335 for (unsigned int d2 = 0; d2 < dim; ++d2)
336 for (unsigned int d3 = 0; d3 < dim; ++d3)
337 {
338 // Derivative
339 // order for each
340 // direction
341 std::vector<unsigned int> deriv_order(dim, 0);
342 for (unsigned int x = 0; x < dim; ++x)
343 {
344 if (d1 == x)
345 ++deriv_order[x];
346 if (d2 == x)
347 ++deriv_order[x];
348 if (d3 == x)
349 ++deriv_order[x];
350 }
351
352 third_derivatives[k2][d1][d2][d3] =
353 v[0][ix][deriv_order[0]] *
354 ((dim > 1) ? v[1][iy][deriv_order[1]] : 1.) *
355 ((dim > 2) ? v[2][iz][deriv_order[2]] : 1.);
356 }
357 }
358 }
359
360 if (update_4th_derivatives)
361 {
362 unsigned int k = 0;
363
364 for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
365 for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
366 for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
367 {
368 const unsigned int k2 = index_map_inverse[k++];
369 for (unsigned int d1 = 0; d1 < dim; ++d1)
370 for (unsigned int d2 = 0; d2 < dim; ++d2)
371 for (unsigned int d3 = 0; d3 < dim; ++d3)
372 for (unsigned int d4 = 0; d4 < dim; ++d4)
373 {
374 // Derivative
375 // order for each
376 // direction
377 std::vector<unsigned int> deriv_order(dim, 0);
378 for (unsigned int x = 0; x < dim; ++x)
379 {
380 if (d1 == x)
381 ++deriv_order[x];
382 if (d2 == x)
383 ++deriv_order[x];
384 if (d3 == x)
385 ++deriv_order[x];
386 if (d4 == x)
387 ++deriv_order[x];
388 }
389
390 fourth_derivatives[k2][d1][d2][d3][d4] =
391 v[0][ix][deriv_order[0]] *
392 ((dim > 1) ? v[1][iy][deriv_order[1]] : 1.) *
393 ((dim > 2) ? v[2][iz][deriv_order[2]] : 1.);
394 }
395 }
396 }
397}
398
399
400
401template <int dim>
402std::unique_ptr<ScalarPolynomialsBase<dim>>
404{
405 return std::make_unique<PolynomialSpace<dim>>(*this);
406}
407
408
409template class PolynomialSpace<1>;
410template class PolynomialSpace<2>;
411template class PolynomialSpace<3>;
412
Definition point.h:111
double compute_value(const unsigned int i, const Point< dim > &p) const override
Tensor< 1, dim > compute_grad(const unsigned int i, const Point< dim > &p) const override
void set_numbering(const std::vector< unsigned int > &renumber)
static unsigned int n_polynomials(const unsigned int n)
std::array< unsigned int, dim > compute_index(const unsigned int n) const
Tensor< 2, dim > compute_grad_grad(const unsigned int i, const Point< dim > &p) const override
virtual std::unique_ptr< ScalarPolynomialsBase< dim > > clone() 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
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch2(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
@ update_values
Shape function values.
@ update_3rd_derivatives
Third derivatives of shape functions.
#define DEAL_II_ASSERT_UNREACHABLE()
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static const unsigned int invalid_unsigned_int
Definition types.h:220