www.mooseframework.org
SplineInterpolationBase.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
16 #include "MooseError.h"
17 #include <limits>
18 
19 const Real SplineInterpolationBase::_deriv_bound = std::numeric_limits<Real>::max();
20 
22 
23 void
24 SplineInterpolationBase::spline(const std::vector<Real> & x,
25  const std::vector<Real> & y,
26  std::vector<Real> & y2,
27  Real yp1 /* = _deriv_bound*/,
28  Real ypn /* = _deriv_bound*/)
29 {
30  auto n = x.size();
31  if (n < 2)
32  mooseError("You must have at least two knots to create a spline.");
33 
34  std::vector<Real> u(n, 0.);
35  y2.assign(n, 0.);
36 
37  if (yp1 >= 1e30)
38  y2[0] = u[0] = 0.;
39  else
40  {
41  y2[0] = -0.5;
42  u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
43  }
44  // decomposition of tri-diagonal algorithm (y2 and u are used for temporary storage)
45  for (decltype(n) i = 1; i < n - 1; i++)
46  {
47  Real sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
48  Real p = sig * y2[i - 1] + 2.0;
49  y2[i] = (sig - 1.0) / p;
50  u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
51  u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
52  }
53 
54  Real qn, un;
55  if (ypn >= 1e30)
56  qn = un = 0.;
57  else
58  {
59  qn = 0.5;
60  un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
61  }
62 
63  y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.);
64  // back substitution
65  for (auto k = n - 1; k >= 1; k--)
66  y2[k - 1] = y2[k - 1] * y2[k] + u[k - 1];
67 }
68 
69 void
70 SplineInterpolationBase::findInterval(const std::vector<Real> & x,
71  Real x_int,
72  unsigned int & klo,
73  unsigned int & khi) const
74 {
75  klo = 0;
76  mooseAssert(x.size() >= 2, "You must have at least two knots to create a spline.");
77  khi = x.size() - 1;
78  while (khi - klo > 1)
79  {
80  unsigned int k = (khi + klo) >> 1;
81  if (x[k] > x_int)
82  khi = k;
83  else
84  klo = k;
85  }
86 }
87 
88 void
89 SplineInterpolationBase::computeCoeffs(const std::vector<Real> & x,
90  unsigned int klo,
91  unsigned int khi,
92  Real x_int,
93  Real & h,
94  Real & a,
95  Real & b) const
96 {
97  h = x[khi] - x[klo];
98  if (h == 0)
99  mooseError("The values of x must be distinct");
100  a = (x[khi] - x_int) / h;
101  b = (x_int - x[klo]) / h;
102 }
103 
104 Real
105 SplineInterpolationBase::sample(const std::vector<Real> & x,
106  const std::vector<Real> & y,
107  const std::vector<Real> & y2,
108  Real x_int) const
109 {
110  unsigned int klo, khi;
111  findInterval(x, x_int, klo, khi);
112 
113  Real h, a, b;
114  computeCoeffs(x, klo, khi, x_int, h, a, b);
115 
116  return a * y[klo] + b * y[khi] +
117  ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0;
118 }
119 
120 Real
122  const std::vector<Real> & y,
123  const std::vector<Real> & y2,
124  Real x_int) const
125 {
126  unsigned int klo, khi;
127  findInterval(x, x_int, klo, khi);
128 
129  Real h, a, b;
130  computeCoeffs(x, klo, khi, x_int, h, a, b);
131 
132  return (y[khi] - y[klo]) / h -
133  (((3.0 * a * a - 1.0) * y2[klo] + (3.0 * b * b - 1.0) * -y2[khi]) * h / 6.0);
134 }
135 
136 Real
138  const std::vector<Real> & /*y*/,
139  const std::vector<Real> & y2,
140  Real x_int) const
141 {
142  unsigned int klo, khi;
143  findInterval(x, x_int, klo, khi);
144 
145  Real h, a, b;
146  computeCoeffs(x, klo, khi, x_int, h, a, b);
147 
148  return a * y2[klo] + b * y2[khi];
149 }
void findInterval(const std::vector< Real > &x, Real x_int, unsigned int &klo, unsigned int &khi) const
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 mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
static PetscErrorCode Vec x
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.
Real sampleDerivative(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