www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes | List of all members
SplineInterpolation Class Reference

This class interpolates tabulated functions with cubic splines. More...

#include <SplineInterpolation.h>

Inheritance diagram for SplineInterpolation:
[legend]

Public Member Functions

 SplineInterpolation ()
 
 SplineInterpolation (const std::vector< Real > &x, const std::vector< Real > &y, Real yp1=_deriv_bound, Real ypn=_deriv_bound)
 Construct the object. More...
 
virtual ~SplineInterpolation ()=default
 
void setData (const std::vector< Real > &x, const std::vector< Real > &y, Real yp1=_deriv_bound, Real ypn=_deriv_bound)
 Set the x-, y- values and first derivatives. More...
 
void errorCheck ()
 
Real sample (Real x) const
 This function will take an independent variable input and will return the dependent variable based on the generated fit. More...
 
Real sampleDerivative (Real x) const
 
Real sample2ndDerivative (Real x) const
 
void dumpSampleFile (std::string base_name, std::string x_label="X", std::string y_label="Y", float xmin=0, float xmax=0, float ymin=0, float ymax=0)
 This function will dump GNUPLOT input files that can be run to show the data points and function fits. More...
 
unsigned int getSampleSize ()
 This function returns the size of the array holding the points, i.e. More...
 
Real domain (int i) const
 
Real range (int i) const
 
Real sample (const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
 
Real sampleDerivative (const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
 
Real sample2ndDerivative (const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
 

Protected Member Functions

void solve ()
 
Real sample (const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int, unsigned int klo, unsigned int khi) const
 Sample value at point x_int given the indices of the vector of dependent values that bound the point. More...
 
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. More...
 
void findInterval (const std::vector< Real > &x, Real x_int, unsigned int &klo, unsigned int &khi) const
 
void computeCoeffs (const std::vector< Real > &x, unsigned int klo, unsigned int khi, Real x_int, Real &h, Real &a, Real &b) const
 

Protected Attributes

std::vector< Real > _x
 
std::vector< Real > _y
 
Real _yp1
 boundary conditions More...
 
Real _ypn
 
std::vector< Real > _y2
 Second derivatives. More...
 

Static Protected Attributes

static int _file_number = 0
 
static const Real _deriv_bound = std::numeric_limits<Real>::max()
 

Detailed Description

This class interpolates tabulated functions with cubic splines.

Adopted from Numerical Recipes in C (section 3.3).

Definition at line 21 of file SplineInterpolation.h.

Constructor & Destructor Documentation

◆ SplineInterpolation() [1/2]

SplineInterpolation::SplineInterpolation ( )

Definition at line 20 of file SplineInterpolation.C.

20 {}

◆ SplineInterpolation() [2/2]

SplineInterpolation::SplineInterpolation ( const std::vector< Real > &  x,
const std::vector< Real > &  y,
Real  yp1 = _deriv_bound,
Real  ypn = _deriv_bound 
)

Construct the object.

Parameters
xTabulated function (x-positions)
yTabulated function (y-positions)
yp1First derivative of the interpolating function at point 1
ypnFirst derivative of the interpolating function at point n

If yp1, ypn are not specified or greater or equal that _deriv_bound, we use natural spline

Definition at line 22 of file SplineInterpolation.C.

26  : SplineInterpolationBase(), _x(x), _y(y), _yp1(yp1), _ypn(ypn)
27 {
28  errorCheck();
29  solve();
30 }
Real _yp1
boundary conditions
static PetscErrorCode Vec x
std::vector< Real > _x
std::vector< Real > _y

◆ ~SplineInterpolation()

virtual SplineInterpolation::~SplineInterpolation ( )
virtualdefault

Member Function Documentation

◆ computeCoeffs()

void SplineInterpolationBase::computeCoeffs ( const std::vector< Real > &  x,
unsigned int  klo,
unsigned int  khi,
Real  x_int,
Real &  h,
Real &  a,
Real &  b 
) const
protectedinherited

Definition at line 84 of file SplineInterpolationBase.C.

Referenced by SplineInterpolationBase::sample(), SplineInterpolationBase::sample2ndDerivative(), and SplineInterpolationBase::sampleDerivative().

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 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:208
static PetscErrorCode Vec x

◆ domain()

Real SplineInterpolation::domain ( int  i) const

Definition at line 86 of file SplineInterpolation.C.

87 {
88  return _x[i];
89 }
std::vector< Real > _x

◆ dumpSampleFile()

void SplineInterpolation::dumpSampleFile ( std::string  base_name,
std::string  x_label = "X",
std::string  y_label = "Y",
float  xmin = 0,
float  xmax = 0,
float  ymin = 0,
float  ymax = 0 
)

This function will dump GNUPLOT input files that can be run to show the data points and function fits.

Definition at line 98 of file SplineInterpolation.C.

105 {
106  std::stringstream filename, filename_pts;
107  const unsigned char fill_character = '0';
108  const unsigned int field_width = 4;
109 
110  filename.fill(fill_character);
111  filename << base_name;
112  filename.width(field_width);
113  filename << _file_number << ".plt";
114 
115  filename_pts.fill(fill_character);
116  filename_pts << base_name << "_pts";
117  filename_pts.width(field_width);
118  filename_pts << _file_number << ".dat";
119 
120  /* First dump the GNUPLOT file with the Piecewise Linear Equations */
121  std::ofstream out(filename.str().c_str());
122  out.precision(15);
123  out.fill(fill_character);
124 
125  out << "set terminal postscript color enhanced\n"
126  << "set output \"" << base_name;
127  out.width(field_width);
128  out << _file_number << ".eps\"\n"
129  << "set xlabel \"" << x_label << "\"\n"
130  << "set ylabel \"" << y_label << "\"\n";
131  if (xmin != 0 && xmax != 0)
132  out << "set xrange [" << xmin << ":" << xmax << "]\n";
133  if (ymin != 0 && ymax != 0)
134  out << "set yrange [" << ymin << ":" << ymax << "]\n";
135  out << "set key left top\n"
136  << "f(x)=";
137 
138  for (unsigned int i = 1; i < _x.size(); ++i)
139  {
140  Real m = (_y[i] - _y[i - 1]) / (_x[i] - _x[i - 1]);
141  Real b = (_y[i] - m * _x[i]);
142 
143  out << _x[i - 1] << "<=x && x<" << _x[i] << " ? " << m << "*x+(" << b << ") : ";
144  }
145  out << " 1/0\n";
146 
147  out << "\nplot f(x) with lines, '" << filename_pts.str() << "' using 1:2 title \"Points\"\n";
148  out.close();
149 
150  libmesh_assert(_x.size() == _y.size());
151 
152  out.open(filename_pts.str().c_str());
153  if (out.fail())
154  throw std::runtime_error(std::string("Unable to open file ") + filename_pts.str());
155 
156  /* Next dump the data points into a seperate file */
157  for (unsigned int i = 0; i < _x.size(); ++i)
158  out << _x[i] << " " << _y[i] << "\n";
159  out << std::endl;
160 
161  ++_file_number;
162  out.close();
163 }
PetscInt m
std::vector< Real > _x
std::vector< Real > _y

◆ errorCheck()

void SplineInterpolation::errorCheck ( )

Definition at line 47 of file SplineInterpolation.C.

Referenced by setData(), and SplineInterpolation().

48 {
49  if (_x.size() != _y.size())
50  mooseError("SplineInterpolation: vectors are not the same length");
51 
52  bool error = false;
53  for (unsigned i = 0; !error && i + 1 < _x.size(); ++i)
54  if (_x[i] >= _x[i + 1])
55  error = true;
56 
57  if (error)
58  mooseError("x-values are not strictly increasing");
59 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:208
std::vector< Real > _x
std::vector< Real > _y

◆ findInterval()

void SplineInterpolationBase::findInterval ( const std::vector< Real > &  x,
Real  x_int,
unsigned int &  klo,
unsigned int &  khi 
) const
protectedinherited

Definition at line 65 of file SplineInterpolationBase.C.

Referenced by BicubicSplineInterpolation::constructColumnSpline(), BicubicSplineInterpolation::constructRowSpline(), SplineInterpolationBase::sample(), SplineInterpolationBase::sample2ndDerivative(), and SplineInterpolationBase::sampleDerivative().

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 }
static PetscErrorCode Vec x

◆ getSampleSize()

unsigned int SplineInterpolation::getSampleSize ( )

This function returns the size of the array holding the points, i.e.

the number of sample points

Definition at line 166 of file SplineInterpolation.C.

167 {
168  return _x.size();
169 }
std::vector< Real > _x

◆ range()

Real SplineInterpolation::range ( int  i) const

Definition at line 92 of file SplineInterpolation.C.

93 {
94  return _y[i];
95 }
std::vector< Real > _y

◆ sample() [1/3]

Real SplineInterpolationBase::sample ( const std::vector< Real > &  x,
const std::vector< Real > &  y,
const std::vector< Real > &  y2,
Real  x_int 
) const
inherited

Definition at line 100 of file SplineInterpolationBase.C.

Referenced by BicubicSplineInterpolation::constructColumnSpline(), BicubicSplineInterpolation::constructRowSpline(), sample(), BicubicSplineInterpolation::sample(), and BicubicSplineInterpolation::sampleValueAndDerivatives().

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 }
Real sample(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

◆ sample() [2/3]

Real SplineInterpolation::sample ( Real  x) const

This function will take an independent variable input and will return the dependent variable based on the generated fit.

Definition at line 68 of file SplineInterpolation.C.

Referenced by SplineFunction::value().

69 {
71 }
Real sample(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
static PetscErrorCode Vec x
std::vector< Real > _y2
Second derivatives.
std::vector< Real > _x
std::vector< Real > _y

◆ sample() [3/3]

Real SplineInterpolationBase::sample ( const std::vector< Real > &  x,
const std::vector< Real > &  y,
const std::vector< Real > &  y2,
Real  x_int,
unsigned int  klo,
unsigned int  khi 
) const
protectedinherited

Sample value at point x_int given the indices of the vector of dependent values that bound the point.

This method is useful in bicubic spline interpolation, where several spline evaluations are needed to sample from a 2D point.

Definition at line 143 of file SplineInterpolationBase.C.

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 computeCoeffs(const std::vector< Real > &x, unsigned int klo, unsigned int khi, Real x_int, Real &h, Real &a, Real &b) const
static PetscErrorCode Vec x

◆ sample2ndDerivative() [1/2]

Real SplineInterpolationBase::sample2ndDerivative ( const std::vector< Real > &  x,
const std::vector< Real > &  y,
const std::vector< Real > &  y2,
Real  x_int 
) const
inherited

Definition at line 128 of file SplineInterpolationBase.C.

Referenced by sample2ndDerivative(), and BicubicSplineInterpolation::sample2ndDerivative().

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 }
void computeCoeffs(const std::vector< Real > &x, unsigned int klo, unsigned int khi, Real x_int, Real &h, Real &a, Real &b) const
static PetscErrorCode Vec x
void findInterval(const std::vector< Real > &x, Real x_int, unsigned int &klo, unsigned int &khi) const

◆ sample2ndDerivative() [2/2]

Real SplineInterpolation::sample2ndDerivative ( Real  x) const

Definition at line 80 of file SplineInterpolation.C.

Referenced by SplineFunction::secondDerivative().

81 {
83 }
Real sample2ndDerivative(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
static PetscErrorCode Vec x
std::vector< Real > _y2
Second derivatives.
std::vector< Real > _x
std::vector< Real > _y

◆ sampleDerivative() [1/2]

Real SplineInterpolationBase::sampleDerivative ( const std::vector< Real > &  x,
const std::vector< Real > &  y,
const std::vector< Real > &  y2,
Real  x_int 
) const
inherited

Definition at line 112 of file SplineInterpolationBase.C.

Referenced by sampleDerivative(), BicubicSplineInterpolation::sampleDerivative(), and BicubicSplineInterpolation::sampleValueAndDerivatives().

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 }
void computeCoeffs(const std::vector< Real > &x, unsigned int klo, unsigned int khi, Real x_int, Real &h, Real &a, Real &b) const
static PetscErrorCode Vec x
void findInterval(const std::vector< Real > &x, Real x_int, unsigned int &klo, unsigned int &khi) const

◆ sampleDerivative() [2/2]

Real SplineInterpolation::sampleDerivative ( Real  x) const

Definition at line 74 of file SplineInterpolation.C.

Referenced by SplineFunction::derivative().

75 {
77 }
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
std::vector< Real > _y2
Second derivatives.
std::vector< Real > _x
std::vector< Real > _y

◆ setData()

void SplineInterpolation::setData ( const std::vector< Real > &  x,
const std::vector< Real > &  y,
Real  yp1 = _deriv_bound,
Real  ypn = _deriv_bound 
)

Set the x-, y- values and first derivatives.

Definition at line 33 of file SplineInterpolation.C.

37 {
38  _x = x;
39  _y = y;
40  _yp1 = yp1;
41  _ypn = ypn;
42  errorCheck();
43  solve();
44 }
Real _yp1
boundary conditions
static PetscErrorCode Vec x
std::vector< Real > _x
std::vector< Real > _y

◆ solve()

void SplineInterpolation::solve ( )
protected

Definition at line 62 of file SplineInterpolation.C.

Referenced by setData(), and SplineInterpolation().

63 {
64  spline(_x, _y, _y2, _yp1, _ypn);
65 }
Real _yp1
boundary conditions
std::vector< Real > _y2
Second derivatives.
std::vector< Real > _x
std::vector< Real > _y
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.

◆ spline()

void SplineInterpolationBase::spline ( const std::vector< Real > &  x,
const std::vector< Real > &  y,
std::vector< Real > &  y2,
Real  yp1 = _deriv_bound,
Real  ypn = _deriv_bound 
)
protectedinherited

This function calculates the second derivatives based on supplied x and y-vectors.

Definition at line 19 of file SplineInterpolationBase.C.

Referenced by BicubicSplineInterpolation::constructColumnSpline(), BicubicSplineInterpolation::constructColumnSplineSecondDerivativeTable(), BicubicSplineInterpolation::constructRowSpline(), BicubicSplineInterpolation::constructRowSplineSecondDerivativeTable(), and solve().

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 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:208
static PetscErrorCode Vec x
PetscInt n

Member Data Documentation

◆ _deriv_bound

const Real SplineInterpolationBase::_deriv_bound = std::numeric_limits<Real>::max()
staticprotectedinherited

Definition at line 75 of file SplineInterpolationBase.h.

Referenced by BicubicSplineInterpolation::errorCheck().

◆ _file_number

int SplineInterpolation::_file_number = 0
staticprotected

Definition at line 92 of file SplineInterpolation.h.

Referenced by dumpSampleFile().

◆ _x

std::vector<Real> SplineInterpolation::_x
protected

◆ _y

std::vector<Real> SplineInterpolation::_y
protected

◆ _y2

std::vector<Real> SplineInterpolation::_y2
protected

Second derivatives.

Definition at line 88 of file SplineInterpolation.h.

Referenced by sample(), sample2ndDerivative(), sampleDerivative(), and solve().

◆ _yp1

Real SplineInterpolation::_yp1
protected

boundary conditions

Definition at line 86 of file SplineInterpolation.h.

Referenced by setData(), and solve().

◆ _ypn

Real SplineInterpolation::_ypn
protected

Definition at line 86 of file SplineInterpolation.h.

Referenced by setData(), and solve().


The documentation for this class was generated from the following files: