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 : #pragma once
11 :
12 : #include <vector>
13 :
14 : #include "libmesh/libmesh_common.h"
15 :
16 : using libMesh::Real;
17 :
18 : /**
19 : * This class interpolates values given a set of data pairs and an abscissa.
20 : * The interpolation is cubic with at least C1 continuity; C2 continuity can
21 : * be violated in favor of ensuring the interpolation is monotonic.
22 : * The algorithm used is laid out in Fritsch and Carlson, SIAM J. Numer. Anal.
23 : * Vol. 17(2) April 1980
24 : */
25 : class MonotoneCubicInterpolation
26 : {
27 : public:
28 : /**
29 : * Empty constructor
30 : */
31 : MonotoneCubicInterpolation();
32 :
33 : /**
34 : * Constructor, Takes two vectors of points for which to apply the fit. One
35 : * should be of the independent variable while the other should be of the
36 : * dependent variable. These values should have a one-to-one correspondence,
37 : * e.g. the vectors must be of the same size.
38 : */
39 : MonotoneCubicInterpolation(const std::vector<Real> & x, const std::vector<Real> & y);
40 :
41 4 : virtual ~MonotoneCubicInterpolation() = default;
42 :
43 : /**
44 : * Method generally used when MonotoneCubicInterpolation object was created
45 : * using the empty constructor. Takes two vectors of points for which to apply
46 : * the fit. One should be of the independent variable while the other should
47 : * be of the dependent variable. These values should have a one-to-one
48 : * correspondence, e.g. the vectors must be of the same size.
49 : */
50 : virtual void setData(const std::vector<Real> & x, const std::vector<Real> & y);
51 :
52 : /**
53 : * This function will take an independent variable input and will return the dependent
54 : * variable based on the generated fit.
55 : */
56 : virtual Real sample(const Real & x) const;
57 :
58 : /**
59 : * This function will take an independent variable input and will return the derivative
60 : * of the dependent variable with respect to the independent variable based on the
61 : * generated fit.
62 : */
63 : virtual Real sampleDerivative(const Real & x) const;
64 :
65 : /**
66 : * This function will take an independent variable input and will return the second
67 : * derivative of the dependent variable with respect to the independent variable
68 : * based on the generated fit. Note that this can be discontinous at the knots.
69 : */
70 : virtual Real sample2ndDerivative(const Real & x) const;
71 :
72 : /**
73 : * This function takes an array of independent variable values and writes a CSV file
74 : * with values corresponding to y, y', and y''. This can be used for sanity checks
75 : * of the interpolation curve.
76 : */
77 : virtual void dumpCSV(std::string filename, const std::vector<Real> & xnew);
78 :
79 : /**
80 : * This method returns the length of the independent variable vector
81 : */
82 : virtual unsigned int getSampleSize();
83 :
84 : protected:
85 : // Error check routines run during initialization
86 : virtual void errorCheck();
87 : Real sign(const Real & x) const;
88 :
89 : // Building blocks of Hermite polynomials
90 : Real phi(const Real & t) const;
91 : Real psi(const Real & t) const;
92 : Real phiPrime(const Real & t) const;
93 : Real psiPrime(const Real & t) const;
94 : Real phiDoublePrime(const Real & t) const;
95 : Real psiDoublePrime(const Real & t) const;
96 :
97 : // Cubic Hermite polynomials
98 : Real h1(const Real & xhi, const Real & xlo, const Real & x) const;
99 : Real h2(const Real & xhi, const Real & xlo, const Real & x) const;
100 : Real h3(const Real & xhi, const Real & xlo, const Real & x) const;
101 : Real h4(const Real & xhi, const Real & xlo, const Real & x) const;
102 : Real h1Prime(const Real & xhi, const Real & xlo, const Real & x) const;
103 : Real h2Prime(const Real & xhi, const Real & xlo, const Real & x) const;
104 : Real h3Prime(const Real & xhi, const Real & xlo, const Real & x) const;
105 : Real h4Prime(const Real & xhi, const Real & xlo, const Real & x) const;
106 : Real h1DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const;
107 : Real h2DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const;
108 : Real h3DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const;
109 : Real h4DoublePrime(const Real & xhi, const Real & xlo, const Real & x) const;
110 :
111 : // Interpolating cubic polynomial and derivatives
112 : virtual Real p(const Real & xhi,
113 : const Real & xlo,
114 : const Real & fhi,
115 : const Real & flo,
116 : const Real & dhi,
117 : const Real & dlo,
118 : const Real & x) const;
119 : virtual Real pPrime(const Real & xhi,
120 : const Real & xlo,
121 : const Real & fhi,
122 : const Real & flo,
123 : const Real & dhi,
124 : const Real & dlo,
125 : const Real & x) const;
126 : virtual Real pDoublePrime(const Real & xhi,
127 : const Real & xlo,
128 : const Real & fhi,
129 : const Real & flo,
130 : const Real & dhi,
131 : const Real & dlo,
132 : const Real & x) const;
133 :
134 : // Algorithm routines
135 : virtual void initialize_derivs();
136 : virtual void modify_derivs(
137 : const Real & alpha, const Real & beta, const Real & delta, Real & yp_lo, Real & yp_hi);
138 : virtual void solve();
139 : virtual void findInterval(const Real & x, unsigned int & klo, unsigned int & khi) const;
140 :
141 : std::vector<Real> _x;
142 : std::vector<Real> _y;
143 : std::vector<Real> _h;
144 : std::vector<Real> _yp;
145 : std::vector<Real> _delta;
146 : std::vector<Real> _alpha;
147 : std::vector<Real> _beta;
148 :
149 : unsigned int _n_knots;
150 : unsigned int _n_intervals;
151 : unsigned int _internal_knots;
152 : };
|
