Reference documentation for deal.II version 9.5.0
Searching...
No Matches
polynomials_hermite.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2023 - 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
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
20
21namespace Polynomials
22{
23 namespace
24 {
25 std::vector<double>
26 hermite_poly_coeffs(const unsigned int regularity, const unsigned int index)
27 {
28 AssertIndexRange(index, 2 * regularity + 2);
29
30 const unsigned int curr_index = index % (regularity + 1);
31 const unsigned int side = (index > regularity) ? 1 : 0;
32
33 // Signed ints are used here to protect against underflow errors
34 const int loop_control_1 = static_cast<int>(regularity + 1 - curr_index);
35 const int loop_control_2 = (side == 1) ?
36 static_cast<int>(curr_index + 1) :
37 static_cast<int>(regularity + 2);
38
39 std::vector<double> poly_coeffs(2 * regularity + 2, 0.0);
40
41 if (side == 1) // right side: g polynomials
42 {
43 int binomial_1 = (curr_index % 2) ? -1 : 1;
44
45 for (int i = 0; i < loop_control_2; ++i)
46 {
47 int inv_binomial = 1;
48
49 for (int j = 0; j < loop_control_1; ++j)
50 {
51 int binomial_2 = 1;
52
53 for (int k = 0; k < j + 1; ++k)
54 {
55 poly_coeffs[regularity + i + k + 1] +=
56 binomial_1 * inv_binomial * binomial_2;
57 binomial_2 *= k - j;
58 binomial_2 /= k + 1;
59 }
60 inv_binomial *= regularity + j + 1;
61 inv_binomial /= j + 1;
62 }
63 // ints used here to protect against underflow errors
64 binomial_1 *= -static_cast<int>(curr_index - i);
65 binomial_1 /= i + 1;
66 }
67 }
68 else // left side: f polynomials
69 {
70 int binomial = 1;
71
72 for (int i = 0; i < loop_control_2; ++i)
73 {
74 int inv_binomial = 1;
75
76 for (int j = 0; j < loop_control_1; ++j)
77 {
78 poly_coeffs[curr_index + i + j] += binomial * inv_binomial;
79 inv_binomial *= regularity + j + 1;
80 inv_binomial /= j + 1;
81 }
82 // Protection needed here against underflow errors
83 binomial *= -static_cast<int>(regularity + 1 - i);
84 binomial /= i + 1;
85 }
86 }
87
88 // rescale coefficients by a factor of 4^curr_index to account for reduced
89 // L2-norms
90 double precond_factor = Utilities::pow(4, curr_index);
91 for (auto &it : poly_coeffs)
92 it *= precond_factor;
93
94 return poly_coeffs;
95 }
96 } // namespace
97
98
99
100 PolynomialsHermite::PolynomialsHermite(const unsigned int regularity,
101 const unsigned int index)
102 : Polynomial<double>(hermite_poly_coeffs(regularity, index))
103 , degree(2 * regularity + 1)
104 , regularity(regularity)
105 , side_index(index % (regularity + 1))
106 , side((index >= regularity + 1) ? 1 : 0)
107 {
108 AssertIndexRange(index, 2 * (regularity + 1));
109 }
110
111
112
113 std::vector<Polynomial<double>>
114 PolynomialsHermite::generate_complete_basis(const unsigned int regularity)
115 {
116 std::vector<Polynomial<double>> polys;
117 const unsigned int sz = 2 * regularity + 2;
118 polys.reserve(sz);
119
120 for (unsigned int i = 0; i < sz; ++i)
121 polys.emplace_back(PolynomialsHermite(regularity, i));
122
123 return polys;
124 }
125} // namespace Polynomials
126
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int regularity)
PolynomialsHermite(const unsigned int regularity, const unsigned int index)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define AssertIndexRange(index, range)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:447