www.mooseframework.org
MonotoneCubicInterpolation.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 
17 #include <fstream>
18 #include <sstream>
19 #include <string>
20 #include <stdexcept>
21 #include <cassert>
22 #include <cmath>
23 
25 
27  const std::vector<Real> & y)
28  : _x(x), _y(y)
29 {
30  errorCheck();
31  solve();
32 }
33 
34 void
35 MonotoneCubicInterpolation::setData(const std::vector<Real> & x, const std::vector<Real> & y)
36 {
37  _x = x;
38  _y = y;
39  errorCheck();
40  solve();
41 }
42 
43 void
45 {
46  if (_x.size() != _y.size())
47  throw std::domain_error("MonotoneCubicInterpolation: x and y vectors are not the same length");
48 
49  bool error = false;
50  for (unsigned i = 0; !error && i + 1 < _x.size(); ++i)
51  if (_x[i] >= _x[i + 1])
52  error = true;
53 
54  if (error)
55  throw std::domain_error("x-values are not strictly increasing");
56 }
57 
58 Real
60 {
61  if (x < 0)
62  return -1;
63  else if (x > 0)
64  return 1;
65  else
66  return 0;
67 }
68 
69 Real
70 MonotoneCubicInterpolation::phi(const Real & t) const
71 {
72  return 3. * t * t - 2. * t * t * t;
73 }
74 
75 Real
77 {
78  return 6. * t - 6. * t * t;
79 }
80 
81 Real
83 {
84  return 6. - 12. * t;
85 }
86 
87 Real
88 MonotoneCubicInterpolation::psi(const Real & t) const
89 {
90  return t * t * t - t * t;
91 }
92 
93 Real
95 {
96  return 3. * t * t - 2. * t;
97 }
98 
99 Real
101 {
102  return 6. * t - 2.;
103 }
104 
105 Real
106 MonotoneCubicInterpolation::h1(const Real & xhi, const Real & xlo, const Real & x) const
107 {
108  Real h = xhi - xlo;
109  Real t = (xhi - x) / h;
110  return phi(t);
111 }
112 
113 Real
114 MonotoneCubicInterpolation::h1Prime(const Real & xhi, const Real & xlo, const Real & x) const
115 {
116  Real h = xhi - xlo;
117  Real t = (xhi - x) / h;
118  Real tPrime = -1. / h;
119  return phiPrime(t) * tPrime;
120 }
121 
122 Real
123 MonotoneCubicInterpolation::h1DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const
124 {
125  Real h = xhi - xlo;
126  Real t = (xhi - x) / h;
127  Real tPrime = -1. / h;
128  return phiDoublePrime(t) * tPrime * tPrime;
129 }
130 
131 Real
132 MonotoneCubicInterpolation::h2(const Real & xhi, const Real & xlo, const Real & x) const
133 {
134  Real h = xhi - xlo;
135  Real t = (x - xlo) / h;
136  return phi(t);
137 }
138 
139 Real
140 MonotoneCubicInterpolation::h2Prime(const Real & xhi, const Real & xlo, const Real & x) const
141 {
142  Real h = xhi - xlo;
143  Real t = (x - xlo) / h;
144  Real tPrime = 1. / h;
145  return phiPrime(t) * tPrime;
146 }
147 
148 Real
149 MonotoneCubicInterpolation::h2DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const
150 {
151  Real h = xhi - xlo;
152  Real t = (x - xlo) / h;
153  Real tPrime = 1. / h;
154  return phiDoublePrime(t) * tPrime * tPrime;
155 }
156 
157 Real
158 MonotoneCubicInterpolation::h3(const Real & xhi, const Real & xlo, const Real & x) const
159 {
160  Real h = xhi - xlo;
161  Real t = (xhi - x) / h;
162  return -h * psi(t);
163 }
164 
165 Real
166 MonotoneCubicInterpolation::h3Prime(const Real & xhi, const Real & xlo, const Real & x) const
167 {
168  Real h = xhi - xlo;
169  Real t = (xhi - x) / h;
170  Real tPrime = -1. / h;
171  return -h * psiPrime(t) * tPrime; // psiPrime(t)
172 }
173 
174 Real
175 MonotoneCubicInterpolation::h3DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const
176 {
177  Real h = xhi - xlo;
178  Real t = (xhi - x) / h;
179  Real tPrime = -1. / h;
180  return psiDoublePrime(t) * tPrime;
181 }
182 
183 Real
184 MonotoneCubicInterpolation::h4(const Real & xhi, const Real & xlo, const Real & x) const
185 {
186  Real h = xhi - xlo;
187  Real t = (x - xlo) / h;
188  return h * psi(t);
189 }
190 
191 Real
192 MonotoneCubicInterpolation::h4Prime(const Real & xhi, const Real & xlo, const Real & x) const
193 {
194  Real h = xhi - xlo;
195  Real t = (x - xlo) / h;
196  Real tPrime = 1. / h;
197  return h * psiPrime(t) * tPrime; // psiPrime(t)
198 }
199 
200 Real
201 MonotoneCubicInterpolation::h4DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const
202 {
203  Real h = xhi - xlo;
204  Real t = (x - xlo) / h;
205  Real tPrime = 1. / h;
206  return psiDoublePrime(t) * tPrime;
207 }
208 
209 Real
211  const Real & xlo,
212  const Real & fhi,
213  const Real & flo,
214  const Real & dhi,
215  const Real & dlo,
216  const Real & x) const
217 {
218  return flo * h1(xhi, xlo, x) + fhi * h2(xhi, xlo, x) + dlo * h3(xhi, xlo, x) +
219  dhi * h4(xhi, xlo, x);
220 }
221 
222 Real
224  const Real & xlo,
225  const Real & fhi,
226  const Real & flo,
227  const Real & dhi,
228  const Real & dlo,
229  const Real & x) const
230 {
231  return flo * h1Prime(xhi, xlo, x) + fhi * h2Prime(xhi, xlo, x) + dlo * h3Prime(xhi, xlo, x) +
232  dhi * h4Prime(xhi, xlo, x);
233 }
234 
235 Real
237  const Real & xlo,
238  const Real & fhi,
239  const Real & flo,
240  const Real & dhi,
241  const Real & dlo,
242  const Real & x) const
243 {
244  return flo * h1DoublePrime(xhi, xlo, x) + fhi * h2DoublePrime(xhi, xlo, x) +
245  dlo * h3DoublePrime(xhi, xlo, x) + dhi * h4DoublePrime(xhi, xlo, x);
246 }
247 
248 void
250 {
251  for (unsigned int i = 1; i < _n_knots - 1; ++i)
252  _yp[i] = (std::pow(_h[i - 1], 2) * _y[i + 1] - std::pow(_h[i], 2) * _y[i - 1] -
253  _y[i] * (_h[i - 1] - _h[i]) * (_h[i - 1] + _h[i])) /
254  (_h[i - 1] * _h[i] * (_h[i - 1] * _h[i]));
255 
256  _yp[0] = (-std::pow(_h[0], 2) * _y[2] - _h[1] * _y[0] * (2 * _h[0] + _h[1]) +
257  _y[1] * std::pow(_h[0] + _h[1], 2)) /
258  (_h[0] * _h[1] * (_h[0] + _h[1]));
259 
260  Real hlast = _h[_n_intervals - 1];
261  Real hsecond = _h[_n_intervals - 2];
262  Real ylast = _y[_n_knots - 1];
263  Real ysecond = _y[_n_knots - 2];
264  Real ythird = _y[_n_knots - 3];
265  _yp[_n_knots - 1] = (hsecond * ylast * (hsecond + 2 * hlast) + std::pow(hlast, 2) * ythird -
266  ysecond * std::pow(hsecond + hlast, 2)) /
267  (hsecond * hlast * (hsecond + hlast));
268 }
269 
270 void
272  const Real & alpha, const Real & beta, const Real & delta, Real & yp_lo, Real & yp_hi)
273 {
274  Real tau = 3. / std::sqrt(std::pow(alpha, 2) + std::pow(beta, 2));
275  Real alpha_star = alpha * tau;
276  Real beta_star = beta * tau;
277  yp_lo = alpha_star * delta;
278  yp_hi = beta_star * delta;
279 }
280 
281 void
283 {
284  _n_knots = _x.size(), _n_intervals = _x.size() - 1, _internal_knots = _x.size() - 2;
285  _h.resize(_n_intervals);
286  _yp.resize(_n_knots);
287  _delta.resize(_n_intervals);
288  _alpha.resize(_n_intervals);
289  _beta.resize(_n_intervals);
290 
291  for (unsigned int i = 0; i < _n_intervals; ++i)
292  _h[i] = _x[i + 1] - _x[i];
293 
295  for (unsigned int i = 0; i < _n_intervals; ++i)
296  _delta[i] = (_y[i + 1] - _y[i]) / _h[i];
297  if (sign(_delta[0]) != sign(_yp[0]))
298  _yp[0] = 0;
299  if (sign(_delta[_n_intervals - 1]) != sign(_yp[_n_knots - 1]))
300  _yp[_n_knots - 1] = 0;
301  for (unsigned int i = 0; i < _internal_knots; ++i)
302  if (sign(_delta[i + 1]) == 0 || sign(_delta[i]) == 0 || sign(_delta[i + 1]) != sign(_delta[i]))
303  _yp[1 + i] = 0;
304 
305  for (unsigned int i = 0; i < _n_intervals; ++i)
306  {
307  // Test for zero slope
308  if (_yp[i] == 0)
309  _alpha[i] = 0;
310  else if (_delta[i] == 0)
311  _alpha[i] = 4;
312  else
313  _alpha[i] = _yp[i] / _delta[i];
314 
315  // Test for zero slope
316  if (_yp[i + 1] == 0)
317  _beta[i] = 0;
318  else if (_delta[i] == 0)
319  _beta[i] = 4;
320  else
321  _beta[i] = _yp[i + 1] / _delta[i];
322 
323  // check if outside region of monotonicity
324  if (std::pow(_alpha[i], 2) + std::pow(_beta[i], 2) > 9.)
325  modify_derivs(_alpha[i], _beta[i], _delta[i], _yp[i], _yp[i + 1]);
326  }
327 }
328 
329 void
331  unsigned int & klo,
332  unsigned int & khi) const
333 {
334  klo = 0;
335  khi = _n_knots - 1;
336  while (khi - klo > 1)
337  {
338  unsigned int k = (khi + klo) >> 1;
339  if (_x[k] > x)
340  khi = k;
341  else
342  klo = k;
343  }
344 }
345 
346 Real
348 {
349  // sanity check (empty MontoneCubicInterpolations are constructable
350  // so we cannot put this into the errorCheck)
351  assert(_x.size() > 0);
352 
353  unsigned int klo, khi;
354  findInterval(x, klo, khi);
355  return p(_x[khi], _x[klo], _y[khi], _y[klo], _yp[khi], _yp[klo], x);
356 }
357 
358 Real
360 {
361  unsigned int klo, khi;
362  findInterval(x, klo, khi);
363  return pPrime(_x[khi], _x[klo], _y[khi], _y[klo], _yp[khi], _yp[klo], x);
364 }
365 
366 Real
368 {
369  unsigned int klo, khi;
370  findInterval(x, klo, khi);
371  return pDoublePrime(_x[khi], _x[klo], _y[khi], _y[klo], _yp[khi], _yp[klo], x);
372 }
373 
374 void
375 MonotoneCubicInterpolation::dumpCSV(std::string filename, const std::vector<Real> & xnew)
376 {
377  unsigned int n = xnew.size();
378  std::vector<Real> ynew(n), ypnew(n), yppnew(n);
379 
380  std::ofstream out(filename.c_str());
381  for (unsigned int i = 0; i < n; ++i)
382  {
383  ynew[i] = sample(xnew[i]);
384  ypnew[i] = sampleDerivative(xnew[i]);
385  yppnew[i] = sample2ndDerivative(xnew[i]);
386  out << xnew[i] << ", " << ynew[i] << ", " << ypnew[i] << ", " << yppnew[i] << "\n";
387  }
388  out.close();
389 }
390 
391 unsigned int
393 {
394  return _x.size();
395 }
Real psiPrime(const Real &t) const
Real h2(const Real &xhi, const Real &xlo, const Real &x) const
Real phiPrime(const Real &t) const
virtual unsigned int getSampleSize()
This method returns the length of the independent variable vector.
Real h2DoublePrime(const Real &xhi, const Real &xlo, const Real &x) const
Real h3(const Real &xhi, const Real &xlo, const Real &x) const
Real phi(const Real &t) const
virtual Real pPrime(const Real &xhi, const Real &xlo, const Real &fhi, const Real &flo, const Real &dhi, const Real &dlo, const Real &x) const
static PetscErrorCode Vec x
Real psi(const Real &t) const
virtual void modify_derivs(const Real &alpha, const Real &beta, const Real &delta, Real &yp_lo, Real &yp_hi)
virtual Real sample2ndDerivative(const Real &x) const
This function will take an independent variable input and will return the second derivative of the de...
virtual void setData(const std::vector< Real > &x, const std::vector< Real > &y)
Method generally used when MonotoneCubicInterpolation object was created using the empty constructor...
Real h1DoublePrime(const Real &xhi, const Real &xlo, const Real &x) const
Real h1(const Real &xhi, const Real &xlo, const Real &x) const
MonotoneCubicInterpolation()
Empty constructor.
Real h4(const Real &xhi, const Real &xlo, const Real &x) const
virtual void findInterval(const Real &x, unsigned int &klo, unsigned int &khi) const
Real phiDoublePrime(const Real &t) const
Real h4DoublePrime(const Real &xhi, const Real &xlo, const Real &x) const
Real h1Prime(const Real &xhi, const Real &xlo, const Real &x) const
virtual Real sampleDerivative(const Real &x) const
This function will take an independent variable input and will return the derivative of the dependent...
virtual Real sample(const Real &x) const
This function will take an independent variable input and will return the dependent variable based on...
virtual Real pDoublePrime(const Real &xhi, const Real &xlo, const Real &fhi, const Real &flo, const Real &dhi, const Real &dlo, const Real &x) const
PetscInt n
virtual void dumpCSV(std::string filename, const std::vector< Real > &xnew)
This function takes an array of independent variable values and writes a CSV file with values corresp...
Real h3DoublePrime(const Real &xhi, const Real &xlo, const Real &x) const
virtual Real p(const Real &xhi, const Real &xlo, const Real &fhi, const Real &flo, const Real &dhi, const Real &dlo, const Real &x) const
Real h4Prime(const Real &xhi, const Real &xlo, const Real &x) const
Real h3Prime(const Real &xhi, const Real &xlo, const Real &x) const
Real psiDoublePrime(const Real &t) const
Real sign(const Real &x) const
Real h2Prime(const Real &xhi, const Real &xlo, const Real &x) const