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