Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 :
10 : #include "MonotoneCubicInterpolation.h"
11 : #include "Conversion.h"
12 :
13 : #include <fstream>
14 : #include <sstream>
15 : #include <string>
16 : #include <stdexcept>
17 : #include <cassert>
18 : #include <cmath>
19 :
20 0 : MonotoneCubicInterpolation::MonotoneCubicInterpolation() {}
21 :
22 7 : MonotoneCubicInterpolation::MonotoneCubicInterpolation(const std::vector<Real> & x,
23 7 : const std::vector<Real> & y)
24 7 : : _x(x), _y(y)
25 : {
26 7 : errorCheck();
27 4 : solve();
28 25 : }
29 :
30 : void
31 0 : MonotoneCubicInterpolation::setData(const std::vector<Real> & x, const std::vector<Real> & y)
32 : {
33 0 : _x = x;
34 0 : _y = y;
35 0 : errorCheck();
36 0 : solve();
37 0 : }
38 :
39 : void
40 7 : MonotoneCubicInterpolation::errorCheck()
41 : {
42 7 : if (_x.size() != _y.size())
43 1 : throw std::domain_error("MonotoneCubicInterpolation: x and y vectors are not the same length");
44 :
45 6 : if (_x.size() < 3)
46 1 : throw std::domain_error("MonotoneCubicInterpolation: " + Moose::stringify(_x.size()) +
47 2 : " points is not enough data for a cubic interpolation");
48 :
49 5 : bool error = false;
50 36 : for (unsigned i = 0; !error && i + 1 < _x.size(); ++i)
51 31 : if (_x[i] >= _x[i + 1])
52 1 : error = true;
53 :
54 5 : if (error)
55 1 : throw std::domain_error("x-values are not strictly increasing");
56 4 : }
57 :
58 : Real
59 98 : MonotoneCubicInterpolation::sign(const Real & x) const
60 : {
61 98 : if (x < 0)
62 12 : return -1;
63 86 : else if (x > 0)
64 75 : return 1;
65 : else
66 11 : return 0;
67 : }
68 :
69 : Real
70 44 : MonotoneCubicInterpolation::phi(const Real & t) const
71 : {
72 44 : return 3. * t * t - 2. * t * t * t;
73 : }
74 :
75 : Real
76 788 : MonotoneCubicInterpolation::phiPrime(const Real & t) const
77 : {
78 788 : return 6. * t - 6. * t * t;
79 : }
80 :
81 : Real
82 22 : MonotoneCubicInterpolation::phiDoublePrime(const Real & t) const
83 : {
84 22 : return 6. - 12. * t;
85 : }
86 :
87 : Real
88 44 : MonotoneCubicInterpolation::psi(const Real & t) const
89 : {
90 44 : return t * t * t - t * t;
91 : }
92 :
93 : Real
94 788 : MonotoneCubicInterpolation::psiPrime(const Real & t) const
95 : {
96 788 : return 3. * t * t - 2. * t;
97 : }
98 :
99 : Real
100 22 : MonotoneCubicInterpolation::psiDoublePrime(const Real & t) const
101 : {
102 22 : return 6. * t - 2.;
103 : }
104 :
105 : Real
106 22 : MonotoneCubicInterpolation::h1(const Real & xhi, const Real & xlo, const Real & x) const
107 : {
108 22 : Real h = xhi - xlo;
109 22 : Real t = (xhi - x) / h;
110 44 : return phi(t);
111 : }
112 :
113 : Real
114 394 : MonotoneCubicInterpolation::h1Prime(const Real & xhi, const Real & xlo, const Real & x) const
115 : {
116 394 : Real h = xhi - xlo;
117 394 : Real t = (xhi - x) / h;
118 394 : Real tPrime = -1. / h;
119 394 : return phiPrime(t) * tPrime;
120 : }
121 :
122 : Real
123 11 : MonotoneCubicInterpolation::h1DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const
124 : {
125 11 : Real h = xhi - xlo;
126 11 : Real t = (xhi - x) / h;
127 11 : Real tPrime = -1. / h;
128 11 : return phiDoublePrime(t) * tPrime * tPrime;
129 : }
130 :
131 : Real
132 22 : MonotoneCubicInterpolation::h2(const Real & xhi, const Real & xlo, const Real & x) const
133 : {
134 22 : Real h = xhi - xlo;
135 22 : Real t = (x - xlo) / h;
136 44 : return phi(t);
137 : }
138 :
139 : Real
140 394 : MonotoneCubicInterpolation::h2Prime(const Real & xhi, const Real & xlo, const Real & x) const
141 : {
142 394 : Real h = xhi - xlo;
143 394 : Real t = (x - xlo) / h;
144 394 : Real tPrime = 1. / h;
145 394 : return phiPrime(t) * tPrime;
146 : }
147 :
148 : Real
149 11 : MonotoneCubicInterpolation::h2DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const
150 : {
151 11 : Real h = xhi - xlo;
152 11 : Real t = (x - xlo) / h;
153 11 : Real tPrime = 1. / h;
154 11 : return phiDoublePrime(t) * tPrime * tPrime;
155 : }
156 :
157 : Real
158 22 : MonotoneCubicInterpolation::h3(const Real & xhi, const Real & xlo, const Real & x) const
159 : {
160 22 : Real h = xhi - xlo;
161 22 : Real t = (xhi - x) / h;
162 22 : return -h * psi(t);
163 : }
164 :
165 : Real
166 394 : MonotoneCubicInterpolation::h3Prime(const Real & xhi, const Real & xlo, const Real & x) const
167 : {
168 394 : Real h = xhi - xlo;
169 394 : Real t = (xhi - x) / h;
170 394 : Real tPrime = -1. / h;
171 394 : return -h * psiPrime(t) * tPrime; // psiPrime(t)
172 : }
173 :
174 : Real
175 11 : MonotoneCubicInterpolation::h3DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const
176 : {
177 11 : Real h = xhi - xlo;
178 11 : Real t = (xhi - x) / h;
179 11 : Real tPrime = -1. / h;
180 11 : return psiDoublePrime(t) * tPrime;
181 : }
182 :
183 : Real
184 22 : MonotoneCubicInterpolation::h4(const Real & xhi, const Real & xlo, const Real & x) const
185 : {
186 22 : Real h = xhi - xlo;
187 22 : Real t = (x - xlo) / h;
188 22 : return h * psi(t);
189 : }
190 :
191 : Real
192 394 : MonotoneCubicInterpolation::h4Prime(const Real & xhi, const Real & xlo, const Real & x) const
193 : {
194 394 : Real h = xhi - xlo;
195 394 : Real t = (x - xlo) / h;
196 394 : Real tPrime = 1. / h;
197 394 : return h * psiPrime(t) * tPrime; // psiPrime(t)
198 : }
199 :
200 : Real
201 11 : MonotoneCubicInterpolation::h4DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const
202 : {
203 11 : Real h = xhi - xlo;
204 11 : Real t = (x - xlo) / h;
205 11 : Real tPrime = 1. / h;
206 11 : return psiDoublePrime(t) * tPrime;
207 : }
208 :
209 : Real
210 22 : MonotoneCubicInterpolation::p(const Real & xhi,
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 22 : return flo * h1(xhi, xlo, x) + fhi * h2(xhi, xlo, x) + dlo * h3(xhi, xlo, x) +
219 22 : dhi * h4(xhi, xlo, x);
220 : }
221 :
222 : Real
223 394 : MonotoneCubicInterpolation::pPrime(const Real & xhi,
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 394 : return flo * h1Prime(xhi, xlo, x) + fhi * h2Prime(xhi, xlo, x) + dlo * h3Prime(xhi, xlo, x) +
232 394 : dhi * h4Prime(xhi, xlo, x);
233 : }
234 :
235 : Real
236 11 : MonotoneCubicInterpolation::pDoublePrime(const Real & xhi,
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 11 : return flo * h1DoublePrime(xhi, xlo, x) + fhi * h2DoublePrime(xhi, xlo, x) +
245 11 : dlo * h3DoublePrime(xhi, xlo, x) + dhi * h4DoublePrime(xhi, xlo, x);
246 : }
247 :
248 : void
249 4 : MonotoneCubicInterpolation::initialize_derivs()
250 : {
251 30 : for (unsigned int i = 1; i < _n_knots - 1; ++i)
252 52 : _yp[i] = (std::pow(_h[i - 1], 2) * _y[i + 1] - std::pow(_h[i], 2) * _y[i - 1] -
253 26 : _y[i] * (_h[i - 1] - _h[i]) * (_h[i - 1] + _h[i])) /
254 26 : (_h[i - 1] * _h[i] * (_h[i - 1] * _h[i]));
255 :
256 4 : _yp[0] = (-std::pow(_h[0], 2) * _y[2] - _h[1] * _y[0] * (2 * _h[0] + _h[1]) +
257 4 : _y[1] * std::pow(_h[0] + _h[1], 2)) /
258 4 : (_h[0] * _h[1] * (_h[0] + _h[1]));
259 :
260 4 : Real hlast = _h[_n_intervals - 1];
261 4 : Real hsecond = _h[_n_intervals - 2];
262 4 : Real ylast = _y[_n_knots - 1];
263 4 : Real ysecond = _y[_n_knots - 2];
264 4 : Real ythird = _y[_n_knots - 3];
265 4 : _yp[_n_knots - 1] = (hsecond * ylast * (hsecond + 2 * hlast) + std::pow(hlast, 2) * ythird -
266 4 : ysecond * std::pow(hsecond + hlast, 2)) /
267 4 : (hsecond * hlast * (hsecond + hlast));
268 4 : }
269 :
270 : void
271 3 : MonotoneCubicInterpolation::modify_derivs(
272 : const Real & alpha, const Real & beta, const Real & delta, Real & yp_lo, Real & yp_hi)
273 : {
274 3 : Real tau = 3. / std::sqrt(std::pow(alpha, 2) + std::pow(beta, 2));
275 3 : Real alpha_star = alpha * tau;
276 3 : Real beta_star = beta * tau;
277 3 : yp_lo = alpha_star * delta;
278 3 : yp_hi = beta_star * delta;
279 3 : }
280 :
281 : void
282 4 : MonotoneCubicInterpolation::solve()
283 : {
284 4 : _n_knots = _x.size(), _n_intervals = _x.size() - 1, _internal_knots = _x.size() - 2;
285 4 : _h.resize(_n_intervals);
286 4 : _yp.resize(_n_knots);
287 4 : _delta.resize(_n_intervals);
288 4 : _alpha.resize(_n_intervals);
289 4 : _beta.resize(_n_intervals);
290 :
291 34 : for (unsigned int i = 0; i < _n_intervals; ++i)
292 30 : _h[i] = _x[i + 1] - _x[i];
293 :
294 4 : initialize_derivs();
295 34 : for (unsigned int i = 0; i < _n_intervals; ++i)
296 30 : _delta[i] = (_y[i + 1] - _y[i]) / _h[i];
297 4 : if (sign(_delta[0]) != sign(_yp[0]))
298 1 : _yp[0] = 0;
299 4 : if (sign(_delta[_n_intervals - 1]) != sign(_yp[_n_knots - 1]))
300 0 : _yp[_n_knots - 1] = 0;
301 30 : for (unsigned int i = 0; i < _internal_knots; ++i)
302 26 : if (sign(_delta[i + 1]) == 0 || sign(_delta[i]) == 0 || sign(_delta[i + 1]) != sign(_delta[i]))
303 12 : _yp[1 + i] = 0;
304 :
305 34 : for (unsigned int i = 0; i < _n_intervals; ++i)
306 : {
307 : // Test for zero slope
308 30 : if (_yp[i] == 0)
309 14 : _alpha[i] = 0;
310 16 : else if (_delta[i] == 0)
311 0 : _alpha[i] = 4;
312 : else
313 16 : _alpha[i] = _yp[i] / _delta[i];
314 :
315 : // Test for zero slope
316 30 : if (_yp[i + 1] == 0)
317 12 : _beta[i] = 0;
318 18 : else if (_delta[i] == 0)
319 0 : _beta[i] = 4;
320 : else
321 18 : _beta[i] = _yp[i + 1] / _delta[i];
322 :
323 : // check if outside region of monotonicity
324 30 : if (std::pow(_alpha[i], 2) + std::pow(_beta[i], 2) > 9.)
325 3 : modify_derivs(_alpha[i], _beta[i], _delta[i], _yp[i], _yp[i + 1]);
326 : }
327 4 : }
328 :
329 : void
330 427 : MonotoneCubicInterpolation::findInterval(const Real & x,
331 : unsigned int & klo,
332 : unsigned int & khi) const
333 : {
334 427 : klo = 0;
335 427 : khi = _n_knots - 1;
336 1865 : while (khi - klo > 1)
337 : {
338 1438 : unsigned int k = (khi + klo) >> 1;
339 1438 : if (_x[k] > x)
340 601 : khi = k;
341 : else
342 837 : klo = k;
343 : }
344 427 : }
345 :
346 : Real
347 22 : MonotoneCubicInterpolation::sample(const Real & x) const
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 22 : findInterval(x, klo, khi);
355 44 : return p(_x[khi], _x[klo], _y[khi], _y[klo], _yp[khi], _yp[klo], x);
356 : }
357 :
358 : Real
359 394 : MonotoneCubicInterpolation::sampleDerivative(const Real & x) const
360 : {
361 : unsigned int klo, khi;
362 394 : findInterval(x, klo, khi);
363 788 : return pPrime(_x[khi], _x[klo], _y[khi], _y[klo], _yp[khi], _yp[klo], x);
364 : }
365 :
366 : Real
367 11 : MonotoneCubicInterpolation::sample2ndDerivative(const Real & x) const
368 : {
369 : unsigned int klo, khi;
370 11 : findInterval(x, klo, khi);
371 22 : return pDoublePrime(_x[khi], _x[klo], _y[khi], _y[klo], _yp[khi], _yp[klo], x);
372 : }
373 :
374 : void
375 0 : MonotoneCubicInterpolation::dumpCSV(std::string filename, const std::vector<Real> & xnew)
376 : {
377 0 : unsigned int n = xnew.size();
378 0 : std::vector<Real> ynew(n), ypnew(n), yppnew(n);
379 :
380 0 : std::ofstream out(filename.c_str());
381 0 : for (unsigned int i = 0; i < n; ++i)
382 : {
383 0 : ynew[i] = sample(xnew[i]);
384 0 : ypnew[i] = sampleDerivative(xnew[i]);
385 0 : yppnew[i] = sample2ndDerivative(xnew[i]);
386 0 : out << xnew[i] << ", " << ynew[i] << ", " << ypnew[i] << ", " << yppnew[i] << "\n";
387 : }
388 :
389 0 : out << std::flush;
390 0 : out.close();
391 0 : }
392 :
393 : unsigned int
394 1 : MonotoneCubicInterpolation::getSampleSize()
395 : {
396 1 : return _x.size();
397 : }
|
