Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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/polynomials_piecewise.h>
17 
18 
19 DEAL_II_NAMESPACE_OPEN
20 
21 
22 
23 namespace 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  {
36  Assert(n_intervals > 0, ExcMessage("No intervals given"));
37  AssertIndexRange(interval, n_intervals);
38  }
39 
40 
41 
42  template <typename number>
43  void
45  std::vector<number> &values) const
46  {
47  Assert(values.size() > 0, ExcZero());
48 
49  value(x, values.size() - 1, values.data());
50  }
51 
52 
53 
54  template <typename number>
55  void
57  const unsigned int n_derivatives,
58  number * values) const
59  {
60  // shift polynomial if necessary
61  number y = x;
62  double derivative_change_sign = 1.;
63  if (n_intervals > 0)
64  {
65  const number step = 1. / n_intervals;
66  // polynomial spans over two intervals
67  if (spans_two_intervals)
68  {
69  const double offset = step * interval;
70  if (x < offset || x > offset + step + step)
71  {
72  for (unsigned int k = 0; k <= n_derivatives; ++k)
73  values[k] = 0;
74  return;
75  }
76  else if (x < offset + step)
77  y = x - offset;
78  else
79  {
80  derivative_change_sign = -1.;
81  y = offset + step + step - x;
82  }
83  }
84  else
85  {
86  const double offset = step * interval;
87  if (x < offset || x > offset + step)
88  {
89  for (unsigned int k = 0; k <= n_derivatives; ++k)
90  values[k] = 0;
91  return;
92  }
93  else
94  y = x - offset;
95  }
96 
97  // on subinterval boundaries, cannot evaluate derivatives properly, so
98  // set them to zero
99  if ((std::abs(y) < 1e-14 &&
100  (interval > 0 || derivative_change_sign == -1.)) ||
101  (std::abs(y - step) < 1e-14 &&
102  (interval < n_intervals - 1 || derivative_change_sign == -1.)))
103  {
104  values[0] = value(x);
105  for (unsigned int d = 1; d <= n_derivatives; ++d)
106  values[d] = 0;
107  return;
108  }
109  }
110 
111  polynomial.value(y, n_derivatives, values);
112 
113  // change sign if necessary
114  for (unsigned int j = 1; j <= n_derivatives; j += 2)
115  values[j] *= derivative_change_sign;
116  }
117 
118 
119 
120  std::vector<PiecewisePolynomial<double>>
122  const unsigned int n_subdivisions,
123  const unsigned int base_degree)
124  {
125  std::vector<Polynomial<double>> p_base =
126  LagrangeEquidistant::generate_complete_basis(base_degree);
127  for (auto &polynomial : p_base)
128  polynomial.scale(n_subdivisions);
129 
130  std::vector<PiecewisePolynomial<double>> p;
131  p.reserve(n_subdivisions * base_degree + 1);
132 
133  p.emplace_back(p_base[0], n_subdivisions, 0, false);
134  for (unsigned int s = 0; s < n_subdivisions; ++s)
135  for (unsigned int i = 0; i < base_degree; ++i)
136  p.emplace_back(p_base[i + 1],
137  n_subdivisions,
138  s,
139  i == (base_degree - 1) && s < n_subdivisions - 1);
140  return p;
141  }
142 
143 } // namespace Polynomials
144 
145 // ------------------ explicit instantiations --------------- //
146 
147 namespace Polynomials
148 {
149  template class PiecewisePolynomial<float>;
150  template class PiecewisePolynomial<double>;
151  template class PiecewisePolynomial<long double>;
152 } // namespace Polynomials
153 
154 DEAL_II_NAMESPACE_CLOSE
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
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()