Reference documentation for deal.II version GIT relicensing-249-g48dc7357c7 2024-03-29 12:30: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
polynomials_piecewise.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2013 - 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
18
20
21
22
23namespace Polynomials
24{
25 template <typename number>
27 const Polynomial<number> &coefficients_on_interval,
28 const unsigned int n_intervals,
29 const unsigned int interval,
30 const bool spans_next_interval)
31 : polynomial(coefficients_on_interval)
32 , n_intervals(n_intervals)
33 , interval(interval)
34 , spans_two_intervals(spans_next_interval)
35 , index(numbers::invalid_unsigned_int)
36 {
37 Assert(n_intervals > 0, ExcMessage("No intervals given"));
39 }
40
41
42
43 template <typename number>
45 const std::vector<Point<1, number>> &points,
46 const unsigned int index)
47 : n_intervals(numbers::invalid_unsigned_int)
48 , interval(numbers::invalid_unsigned_int)
49 , spans_two_intervals(false)
50 , index(index)
51 {
52 Assert(points.size() > 1, ExcMessage("No enough points given!"));
54
55 this->points.resize(points.size());
56 for (unsigned int i = 0; i < points.size(); ++i)
57 this->points[i] = points[i][0];
58
59 this->one_over_lengths.resize(points.size() - 1);
60 for (unsigned int i = 0; i < points.size() - 1; ++i)
61 this->one_over_lengths[i] =
62 number(1.0) / (points[i + 1][0] - points[i][0]);
63 }
64
65
66
67 template <typename number>
68 void
70 std::vector<number> &values) const
71 {
72 Assert(values.size() > 0, ExcZero());
73
74 value(x, values.size() - 1, values.data());
75 }
76
77
78
79 template <typename number>
80 void
82 const unsigned int n_derivatives,
83 number *values) const
84 {
85 if (points.size() > 0)
86 {
87 if (x > points[index])
88 values[0] = std::max<number>(0.0,
89 1.0 - (x - points[index]) *
90 one_over_lengths[index]);
91 else if (x < points[index])
92 values[0] = std::max<number>(0.0,
93 0.0 + (x - points[index - 1]) *
94 one_over_lengths[index - 1]);
95 else
96 values[0] = 1.0;
97
98 if (n_derivatives >= 1)
99 {
100 if ((x > points[index]) && (points[index + 1] >= x))
101 values[1] = -1.0 * one_over_lengths[index];
102 else if ((x < points[index]) && (points[index - 1] <= x))
103 values[1] = +1.0 * one_over_lengths[index - 1];
104 else
105 values[1] = 0.0;
106 }
107
108 // all other derivatives are zero
109 for (unsigned int i = 2; i <= n_derivatives; ++i)
110 values[i] = 0.0;
111
112 return;
113 }
114
115 // shift polynomial if necessary
116 number y = x;
117 double derivative_change_sign = 1.;
118 if (n_intervals > 0)
119 {
120 const number step = 1. / n_intervals;
121 // polynomial spans over two intervals
122 if (spans_two_intervals)
123 {
124 const double offset = step * interval;
125 if (x < offset || x > offset + step + step)
126 {
127 for (unsigned int k = 0; k <= n_derivatives; ++k)
128 values[k] = 0;
129 return;
130 }
131 else if (x < offset + step)
132 y = x - offset;
133 else
134 {
135 derivative_change_sign = -1.;
136 y = offset + step + step - x;
137 }
138 }
139 else
140 {
141 const double offset = step * interval;
142 if (x < offset || x > offset + step)
143 {
144 for (unsigned int k = 0; k <= n_derivatives; ++k)
145 values[k] = 0;
146 return;
147 }
148 else
149 y = x - offset;
150 }
151
152 // on subinterval boundaries, cannot evaluate derivatives properly, so
153 // set them to zero
154 if ((std::abs(y) < 1e-14 &&
155 (interval > 0 || derivative_change_sign == -1.)) ||
156 (std::abs(y - step) < 1e-14 &&
157 (interval < n_intervals - 1 || derivative_change_sign == -1.)))
158 {
159 values[0] = value(x);
160 for (unsigned int d = 1; d <= n_derivatives; ++d)
161 values[d] = 0;
162 return;
163 }
164 }
165
166 polynomial.value(y, n_derivatives, values);
167
168 // change sign if necessary
169 for (unsigned int j = 1; j <= n_derivatives; j += 2)
170 values[j] *= derivative_change_sign;
171 }
172
173
174
175 template <typename number>
176 std::size_t
178 {
179 return (polynomial.memory_consumption() +
182 MemoryConsumption::memory_consumption(spans_two_intervals) +
185 }
186
187
188
189 std::vector<PiecewisePolynomial<double>>
191 const unsigned int n_subdivisions,
192 const unsigned int base_degree)
193 {
194 std::vector<Polynomial<double>> p_base =
196 for (auto &polynomial : p_base)
197 polynomial.scale(n_subdivisions);
198
199 std::vector<PiecewisePolynomial<double>> p;
200 p.reserve(n_subdivisions * base_degree + 1);
201
202 p.emplace_back(p_base[0], n_subdivisions, 0, false);
203 for (unsigned int s = 0; s < n_subdivisions; ++s)
204 for (unsigned int i = 0; i < base_degree; ++i)
205 p.emplace_back(p_base[i + 1],
206 n_subdivisions,
207 s,
208 i == (base_degree - 1) && s < n_subdivisions - 1);
209 return p;
210 }
211
212
213
214 std::vector<PiecewisePolynomial<double>>
216 const std::vector<Point<1>> &points)
217 {
218 std::vector<PiecewisePolynomial<double>> p;
219 p.reserve(points.size());
220
221 for (unsigned int s = 0; s < points.size(); ++s)
222 p.emplace_back(points, s);
223
224 return p;
225 }
226
227} // namespace Polynomials
228
229// ------------------ explicit instantiations --------------- //
230
231namespace Polynomials
232{
233 template class PiecewisePolynomial<float>;
234 template class PiecewisePolynomial<double>;
235 template class PiecewisePolynomial<long double>;
236} // namespace Polynomials
237
Definition point.h:111
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
PiecewisePolynomial(const Polynomial< number > &coefficients_on_interval, const unsigned int n_intervals, const unsigned int interval, const bool spans_next_interval)
number value(const number x) const
virtual std::size_t memory_consumption() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcZero()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
std::vector< PiecewisePolynomial< double > > generate_complete_linear_basis_on_subdivisions(const std::vector< Point< 1 > > &points)
std::vector< PiecewisePolynomial< double > > generate_complete_Lagrange_basis_on_subdivisions(const unsigned int n_subdivisions, const unsigned int base_degree)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)