Reference documentation for deal.II version GIT 574c7c8486 2023-09-22 10:20: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\}}\)
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
9 // Public License as published by the Free Software Foundation; either
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 
24 namespace 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!"));
54  AssertIndexRange(index, points.size());
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 
232 namespace 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:112
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcZero()
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcMessage(std::string arg1)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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:213