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