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