Reference documentation for deal.II version 9.3.3
polynomials_piecewise.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 2019 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)
36 {
37 Assert(n_intervals > 0, ExcMessage("No intervals given"));
39 }
40
41
42
43 template <typename number>
44 void
46 std::vector<number> &values) const
47 {
48 Assert(values.size() > 0, ExcZero());
49
50 value(x, values.size() - 1, values.data());
51 }
52
53
54
55 template <typename number>
56 void
58 const unsigned int n_derivatives,
59 number * values) const
60 {
61 // shift polynomial if necessary
62 number y = x;
63 double derivative_change_sign = 1.;
64 if (n_intervals > 0)
65 {
66 const number step = 1. / n_intervals;
67 // polynomial spans over two intervals
68 if (spans_two_intervals)
69 {
70 const double offset = step * interval;
71 if (x < offset || x > offset + step + step)
72 {
73 for (unsigned int k = 0; k <= n_derivatives; ++k)
74 values[k] = 0;
75 return;
76 }
77 else if (x < offset + step)
78 y = x - offset;
79 else
80 {
81 derivative_change_sign = -1.;
82 y = offset + step + step - x;
83 }
84 }
85 else
86 {
87 const double offset = step * interval;
88 if (x < offset || x > offset + step)
89 {
90 for (unsigned int k = 0; k <= n_derivatives; ++k)
91 values[k] = 0;
92 return;
93 }
94 else
95 y = x - offset;
96 }
97
98 // on subinterval boundaries, cannot evaluate derivatives properly, so
99 // set them to zero
100 if ((std::abs(y) < 1e-14 &&
101 (interval > 0 || derivative_change_sign == -1.)) ||
102 (std::abs(y - step) < 1e-14 &&
103 (interval < n_intervals - 1 || derivative_change_sign == -1.)))
104 {
105 values[0] = value(x);
106 for (unsigned int d = 1; d <= n_derivatives; ++d)
107 values[d] = 0;
108 return;
109 }
110 }
111
112 polynomial.value(y, n_derivatives, values);
113
114 // change sign if necessary
115 for (unsigned int j = 1; j <= n_derivatives; j += 2)
116 values[j] *= derivative_change_sign;
117 }
118
119
120
121 template <typename number>
122 std::size_t
124 {
125 return (polynomial.memory_consumption() +
128 MemoryConsumption::memory_consumption(spans_two_intervals));
129 }
130
131
132
133 std::vector<PiecewisePolynomial<double>>
135 const unsigned int n_subdivisions,
136 const unsigned int base_degree)
137 {
138 std::vector<Polynomial<double>> p_base =
140 for (auto &polynomial : p_base)
141 polynomial.scale(n_subdivisions);
142
143 std::vector<PiecewisePolynomial<double>> p;
144 p.reserve(n_subdivisions * base_degree + 1);
145
146 p.emplace_back(p_base[0], n_subdivisions, 0, false);
147 for (unsigned int s = 0; s < n_subdivisions; ++s)
148 for (unsigned int i = 0; i < base_degree; ++i)
149 p.emplace_back(p_base[i + 1],
150 n_subdivisions,
151 s,
152 i == (base_degree - 1) && s < n_subdivisions - 1);
153 return p;
154 }
155
156} // namespace Polynomials
157
158// ------------------ explicit instantiations --------------- //
159
160namespace Polynomials
161{
162 template class PiecewisePolynomial<float>;
163 template class PiecewisePolynomial<double>;
164 template class PiecewisePolynomial<long double>;
165} // namespace Polynomials
166
static std::vector< Polynomial< double > > generate_complete_basis(const unsigned int degree)
Definition: polynomial.cc:678
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:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcZero()
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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_Lagrange_basis_on_subdivisions(const unsigned int n_subdivisions, const unsigned int base_degree)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)