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 
14 const Real SplineInterpolationBase::_deriv_bound = std::numeric_limits<Real>::max();
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 void
84 SplineInterpolationBase::computeCoeffs(const std::vector<Real> & x,
85  unsigned int klo,
86  unsigned int khi,
87  Real x_int,
88  Real & h,
89  Real & a,
90  Real & b) const
91 {
92  h = x[khi] - x[klo];
93  if (h == 0)
94  mooseError("The values of x must be distinct");
95  a = (x[khi] - x_int) / h;
96  b = (x_int - x[klo]) / h;
97 }
98 
99 Real
100 SplineInterpolationBase::sample(const std::vector<Real> & x,
101  const std::vector<Real> & y,
102  const std::vector<Real> & y2,
103  Real x_int) const
104 {
105  unsigned int klo, khi;
106  findInterval(x, x_int, klo, khi);
107 
108  return sample(x, y, y2, x_int, klo, khi);
109 }
110 
111 Real
113  const std::vector<Real> & y,
114  const std::vector<Real> & y2,
115  Real x_int) const
116 {
117  unsigned int klo, khi;
118  findInterval(x, x_int, klo, khi);
119 
120  Real h, a, b;
121  computeCoeffs(x, klo, khi, x_int, h, a, b);
122 
123  return (y[khi] - y[klo]) / h -
124  (((3.0 * a * a - 1.0) * y2[klo] + (3.0 * b * b - 1.0) * -y2[khi]) * h / 6.0);
125 }
126 
127 Real
129  const std::vector<Real> & /*y*/,
130  const std::vector<Real> & y2,
131  Real x_int) const
132 {
133  unsigned int klo, khi;
134  findInterval(x, x_int, klo, khi);
135 
136  Real h, a, b;
137  computeCoeffs(x, klo, khi, x_int, h, a, b);
138 
139  return a * y2[klo] + b * y2[khi];
140 }
141 
142 Real
143 SplineInterpolationBase::sample(const std::vector<Real> & x,
144  const std::vector<Real> & y,
145  const std::vector<Real> & y2,
146  Real x_int,
147  unsigned int klo,
148  unsigned int khi) const
149 {
150  Real h, a, b;
151  computeCoeffs(x, klo, khi, x_int, h, a, b);
152 
153  return a * y[klo] + b * y[khi] +
154  ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0;
155 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:207
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
void computeCoeffs(const std::vector< Real > &x, unsigned int klo, unsigned int khi, Real x_int, Real &h, Real &a, Real &b) const
Real sampleDerivative(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
static PetscErrorCode Vec x
void findInterval(const std::vector< Real > &x, Real x_int, unsigned int &klo, unsigned int &khi) const
PetscInt n
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.