1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2004 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
#include <deal.II/base/polynomials_p.h>
DEAL_II_NAMESPACE_OPEN
template <int dim>
PolynomialsP<dim>::PolynomialsP(const unsigned int p)
: PolynomialSpace<dim>(
Polynomials::Monomial<double>::generate_complete_basis(p))
, p(p)
{
std::vector<unsigned int> index_map(this->n());
create_polynomial_ordering(index_map);
this->set_numbering(index_map);
}
template <>
void
PolynomialsP<1>::create_polynomial_ordering(
std::vector<unsigned int> &index_map) const
{
Assert(index_map.size() == this->n(),
ExcDimensionMismatch(index_map.size(), this->n()));
// identity
for (unsigned int i = 0; i < this->n(); ++i)
index_map[i] = i;
}
namespace
{
const unsigned int imap2[6][21] = {
{0},
{0, 1, 2},
{0, 1, 3, 4, 2, 5},
{0, 1, 4, 5, 2, 7, 6, 8, 3, 9},
{0, 1, 5, 6, 2, 9, 7, 10, 3, 12, 11, 8, 13, 4, 14},
{0, 1, 6, 7, 2, 11, 8, 12, 3, 15, 13, 9, 16, 4, 18, 14, 17, 10, 19, 5, 20}};
}
template <>
void
PolynomialsP<2>::create_polynomial_ordering(
std::vector<unsigned int> &index_map) const
{
Assert(index_map.size() == this->n(),
ExcDimensionMismatch(index_map.size(), this->n()));
Assert(p <= 5, ExcNotImplemented());
// Given the number i of the
// polynomial in
// $1,x,y,xy,x2,y2,...$,
// index_map[i] gives the number of
// the polynomial in
// PolynomialSpace.
for (unsigned int i = 0; i < this->n(); ++i)
index_map[i] = imap2[p][i];
}
namespace
{
const unsigned int imap3[4][20] = {{0},
{0, 1, 2, 3},
{0, 1, 3, 6, 4, 7, 8, 2, 5, 9},
{0, 1, 4, 10, 5, 11, 13, 2, 7, 16,
14, 6, 12, 8, 15, 17, 18, 3, 9, 19}};
}
template <>
void
PolynomialsP<3>::create_polynomial_ordering(
std::vector<unsigned int> &index_map) const
{
Assert(index_map.size() == this->n(),
ExcDimensionMismatch(index_map.size(), this->n()));
Assert(p <= 3, ExcNotImplemented());
// Given the number i of the
// polynomial in
// $1,x,y,xy,x2,y2,...$,
// index_map[i] gives the number of
// the polynomial in
// PolynomialSpace.
for (unsigned int i = 0; i < this->n(); ++i)
index_map[i] = imap3[p][i];
}
template class PolynomialsP<1>;
template class PolynomialsP<2>;
template class PolynomialsP<3>;
DEAL_II_NAMESPACE_CLOSE
|
