Reference documentation for deal.II version Git 3ff8e288d2 2020-02-20 17:30:00 +0100
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
polynomials_piecewise.cc
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
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 
16 #include <deal.II/base/memory_consumption.h>
17 #include <deal.II/base/polynomials_piecewise.h>
18 
19 
20 DEAL_II_NAMESPACE_OPEN
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)
36  {
37  Assert(n_intervals > 0, ExcMessage("No intervals given"));
38  AssertIndexRange(interval, n_intervals);
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 =
139  LagrangeEquidistant::generate_complete_basis(base_degree);
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 
160 namespace Polynomials
161 {
162  template class PiecewisePolynomial<float>;
163  template class PiecewisePolynomial<double>;
164  template class PiecewisePolynomial<long double>;
165 } // namespace Polynomials
166 
167 DEAL_II_NAMESPACE_CLOSE
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1411
std::vector< PiecewisePolynomial< double > > generate_complete_Lagrange_basis_on_subdivisions(const unsigned int n_subdivisions, const unsigned int base_degree)
PiecewisePolynomial(const Polynomial< number > &coefficients_on_interval, const unsigned int n_intervals, const unsigned int interval, const bool spans_next_interval)
static ::ExceptionBase & ExcZero()
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)