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