www.mooseframework.org
SplineInterpolationBase.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 #include "MooseError.h"
12 #include <limits>
13 
15 
17 
18 void
19 SplineInterpolationBase::spline(const std::vector<Real> & x,
20  const std::vector<Real> & y,
21  std::vector<Real> & y2,
22  Real yp1 /* = _deriv_bound*/,
23  Real ypn /* = _deriv_bound*/)
24 {
25  auto n = x.size();
26  if (n < 2)
27  mooseError("You must have at least two knots to create a spline.");
28 
29  std::vector<Real> u(n, 0.);
30  y2.assign(n, 0.);
31 
32  if (yp1 >= 1e30)
33  y2[0] = u[0] = 0.;
34  else
35  {
36  y2[0] = -0.5;
37  u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
38  }
39  // decomposition of tri-diagonal algorithm (y2 and u are used for temporary storage)
40  for (decltype(n) i = 1; i < n - 1; i++)
41  {
42  Real sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
43  Real p = sig * y2[i - 1] + 2.0;
44  y2[i] = (sig - 1.0) / p;
45  u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
46  u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
47  }
48 
49  Real qn, un;
50  if (ypn >= 1e30)
51  qn = un = 0.;
52  else
53  {
54  qn = 0.5;
55  un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
56  }
57 
58  y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.);
59  // back substitution
60  for (auto k = n - 1; k >= 1; k--)
61  y2[k - 1] = y2[k - 1] * y2[k] + u[k - 1];
62 }
63 
64 void
65 SplineInterpolationBase::findInterval(const std::vector<Real> & x,
66  Real x_int,
67  unsigned int & klo,
68  unsigned int & khi) const
69 {
70  klo = 0;
71  mooseAssert(x.size() >= 2, "You must have at least two knots to create a spline.");
72  khi = x.size() - 1;
73  while (khi - klo > 1)
74  {
75  unsigned int k = (khi + klo) >> 1;
76  if (x[k] > x_int)
77  khi = k;
78  else
79  klo = k;
80  }
81 }
82 
83 template <typename T>
84 void
85 SplineInterpolationBase::computeCoeffs(const std::vector<Real> & x,
86  unsigned int klo,
87  unsigned int khi,
88  const T & x_int,
89  Real & h,
90  T & a,
91  T & b) const
92 {
93  h = x[khi] - x[klo];
94  if (h == 0)
95  mooseError("The values of x must be distinct");
96  a = (x[khi] - x_int) / h;
97  b = (x_int - x[klo]) / h;
98 }
99 
100 Real
101 SplineInterpolationBase::sample(const std::vector<Real> & x,
102  const std::vector<Real> & y,
103  const std::vector<Real> & y2,
104  Real x_int) const
105 {
106  unsigned int klo, khi;
107  findInterval(x, x_int, klo, khi);
108 
109  return sample(x, y, y2, x_int, klo, khi);
110 }
111 
112 ADReal
113 SplineInterpolationBase::sample(const std::vector<Real> & x,
114  const std::vector<Real> & y,
115  const std::vector<Real> & y2,
116  const ADReal & x_int) const
117 {
118  unsigned int klo, khi;
119  findInterval(x, MetaPhysicL::raw_value(x_int), klo, khi);
120 
121  return sample(x, y, y2, x_int, klo, khi);
122 }
123 
124 Real
125 SplineInterpolationBase::sampleDerivative(const std::vector<Real> & x,
126  const std::vector<Real> & y,
127  const std::vector<Real> & y2,
128  Real x_int) const
129 {
130  unsigned int klo, khi;
131  findInterval(x, x_int, klo, khi);
132 
133  Real h, a, b;
134  computeCoeffs(x, klo, khi, x_int, h, a, b);
135 
136  return (y[khi] - y[klo]) / h -
137  (((3.0 * a * a - 1.0) * y2[klo] + (3.0 * b * b - 1.0) * -y2[khi]) * h / 6.0);
138 }
139 
140 Real
142  const std::vector<Real> & /*y*/,
143  const std::vector<Real> & y2,
144  Real x_int) const
145 {
146  unsigned int klo, khi;
147  findInterval(x, x_int, klo, khi);
148 
149  Real h, a, b;
150  computeCoeffs(x, klo, khi, x_int, h, a, b);
151 
152  return a * y2[klo] + b * y2[khi];
153 }
154 
155 template <typename T>
156 T
157 SplineInterpolationBase::sample(const std::vector<Real> & x,
158  const std::vector<Real> & y,
159  const std::vector<Real> & y2,
160  const T & x_int,
161  unsigned int klo,
162  unsigned int khi) const
163 {
164  Real h;
165  T a, b;
166  computeCoeffs(x, klo, khi, x_int, h, a, b);
167 
168  return a * y[klo] + b * y[khi] +
169  ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0;
170 }
171 
172 template Real SplineInterpolationBase::sample<Real>(const std::vector<Real> & x,
173  const std::vector<Real> & y,
174  const std::vector<Real> & y2,
175  const Real & x_int,
176  unsigned int klo,
177  unsigned int khi) const;
178 template ADReal SplineInterpolationBase::sample<ADReal>(const std::vector<Real> & x,
179  const std::vector<Real> & y,
180  const std::vector<Real> & y2,
181  const ADReal & x_int,
182  unsigned int klo,
183  unsigned int khi) const;
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:284
Real sample2ndDerivative(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
Real sample(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
auto raw_value(const Eigen::Map< T > &in)
Definition: ADReal.h:73
void computeCoeffs(const std::vector< Real > &x, unsigned int klo, unsigned int khi, const T &x_int, Real &h, T &a, T &b) const
Real sampleDerivative(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
auto max(const L &left, const R &right)
void findInterval(const std::vector< Real > &x, Real x_int, unsigned int &klo, unsigned int &khi) const
DualReal ADReal
Definition: ADRealForward.h:14
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void spline(const std::vector< Real > &x, const std::vector< Real > &y, std::vector< Real > &y2, Real yp1=_deriv_bound, Real ypn=_deriv_bound)
This function calculates the second derivatives based on supplied x and y-vectors.