https://mooseframework.inl.gov
BicubicSplineInterpolation.C
Go to the documentation of this file.
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 
11 #include "MooseError.h"
12 
14 
16 
18  const std::vector<Real> & x2,
19  const std::vector<std::vector<Real>> & y,
20  const std::vector<Real> & yx11,
21  const std::vector<Real> & yx1n,
22  const std::vector<Real> & yx21,
23  const std::vector<Real> & yx2n)
25  _x1(x1),
26  _x2(x2),
27  _y(y),
28  _yx11(yx11),
29  _yx1n(yx1n),
30  _yx21(yx21),
31  _yx2n(yx2n)
32 {
33  auto n = _x2.size();
34  _row_spline_second_derivs.resize(n);
35  _column_spline_eval.resize(n);
36 
37  auto m = _x1.size();
39  _row_spline_eval.resize(m);
40 
41  errorCheck();
42  solve();
43 }
44 
45 void
46 BicubicSplineInterpolation::setData(const std::vector<Real> & x1,
47  const std::vector<Real> & x2,
48  const std::vector<std::vector<Real>> & y,
49  const std::vector<Real> & yx11,
50  const std::vector<Real> & yx1n,
51  const std::vector<Real> & yx21,
52  const std::vector<Real> & yx2n)
53 {
54  _x1 = x1;
55  _x2 = x2;
56  _y = y;
57  _yx11 = yx11;
58  _yx1n = yx1n;
59  _yx21 = yx21;
60  _yx2n = yx2n;
61 
62  auto n = _x2.size();
63  _row_spline_second_derivs.resize(n);
64  _column_spline_eval.resize(n);
65 
66  auto m = _x1.size();
68  _row_spline_eval.resize(m);
69 
70  errorCheck();
71  solve();
72 }
73 
74 void
76 {
77  auto m = _x1.size(), n = _x2.size();
78 
79  if (_y.size() != m)
80  mooseError("y row dimension does not match the size of x1.");
81  else
82  for (decltype(m) i = 0; i < _y.size(); ++i)
83  if (_y[i].size() != n)
84  mooseError("y column dimension does not match the size of x2.");
85 
86  if (_yx11.empty())
87  _yx11.resize(n, _deriv_bound);
88  else if (_yx11.size() != n)
89  mooseError("The length of the vectors holding the first derivatives of y with respect to x1 "
90  "must match the length of x2.");
91 
92  if (_yx1n.empty())
93  _yx1n.resize(n, _deriv_bound);
94  else if (_yx1n.size() != n)
95  mooseError("The length of the vectors holding the first derivatives of y with respect to x1 "
96  "must match the length of x2.");
97 
98  if (_yx21.empty())
99  _yx21.resize(m, _deriv_bound);
100  else if (_yx21.size() != m)
101  mooseError("The length of the vectors holding the first derivatives of y with respect to x2 "
102  "must match the length of x1.");
103 
104  if (_yx2n.empty())
105  _yx2n.resize(m, _deriv_bound);
106  else if (_yx2n.size() != m)
107  mooseError("The length of the vectors holding the first derivatives of y with respect to x2 "
108  "must match the length of x1.");
109 }
110 
111 void
113 {
114  auto m = _x1.size();
115  _y2_rows.resize(m);
116 
117  if (_yx21.empty())
118  for (decltype(m) j = 0; j < m; ++j)
119  spline(_x2, _y[j], _y2_rows[j]);
120 
121  else
122  for (decltype(m) j = 0; j < m; ++j)
123  spline(_x2, _y[j], _y2_rows[j], _yx21[j], _yx2n[j]);
124 }
125 
126 void
128 {
129  auto m = _x1.size();
130  auto n = _x2.size();
131  _y2_columns.resize(n);
132  _y_trans.resize(n);
133 
134  for (decltype(n) j = 0; j < n; ++j)
135  _y_trans[j].resize(m);
136 
137  // transpose the _y values so the columns can be easily iterated over
138  for (decltype(n) i = 0; i < _y.size(); ++i)
139  for (decltype(n) j = 0; j < _y[0].size(); ++j)
140  _y_trans[j][i] = _y[i][j];
141 
142  if (_yx11.empty())
143  for (decltype(n) j = 0; j < n; ++j)
144  spline(_x1, _y_trans[j], _y2_columns[j]);
145 
146  else
147  for (decltype(n) j = 0; j < n; ++j)
148  spline(_x1, _y_trans[j], _y2_columns[j], _yx11[j], _yx1n[j]);
149 }
150 
151 void
153 {
156 }
157 
158 Real
160  Real x2,
161  Real yx11 /* = _deriv_bound*/,
162  Real yx1n /* = _deriv_bound*/)
163 {
165 
166  // Evaluate newly constructed column spline
168 }
169 
170 Real
172  Real x2,
173  unsigned int deriv_var,
174  Real yp1 /* = _deriv_bound*/,
175  Real ypn /* = _deriv_bound*/)
176 {
177  // Take derivative along x1 axis
178  if (deriv_var == 1)
179  {
181 
182  // Evaluate derivative wrt x1 of newly constructed column spline
185  }
186 
187  // Take derivative along x2 axis
188  else if (deriv_var == 2)
189  {
191 
192  // Evaluate derivative wrt x2 of newly constructed row spline
195  }
196 
197  else
198  mooseError("deriv_var must be either 1 or 2 in BicubicSplineInterpolation");
199 }
200 
201 Real
203  Real x2,
204  unsigned int deriv_var,
205  Real yp1 /* = _deriv_bound*/,
206  Real ypn /* = _deriv_bound*/)
207 {
208  // Take second derivative along x1 axis
209  if (deriv_var == 1)
210  {
212 
213  // Evaluate second derivative wrt x1 of newly constructed column spline
216  }
217 
218  // Take second derivative along x2 axis
219  else if (deriv_var == 2)
220  {
222 
223  // Evaluate second derivative wrt x2 of newly constructed row spline
226  }
227 
228  else
229  mooseError("deriv_var must be either 1 or 2 in BicubicSplineInterpolation");
230 }
231 
232 void
234  Real x2,
235  Real & y,
236  Real & dy1,
237  Real & dy2,
238  Real yx11 /* = _deriv_bound*/,
239  Real yx1n /* = _deriv_bound*/,
240  Real yx21 /* = _deriv_bound*/,
241  Real yx2n /* = _deriv_bound*/)
242 {
247 
251 }
252 
253 void
255  std::vector<Real> & column_spline_eval,
256  std::vector<Real> & row_spline_second_derivs,
257  Real yx11 /*= _deriv_bound*/,
258  Real yx1n /*= _deriv_bound*/)
259 {
260  auto n = _x2.size();
261 
262  // Find the indices that bound the point x1
263  unsigned int klo, khi;
264  findInterval(_x1, x1, klo, khi);
265 
266  // Evaluate n column-splines to get y-values for row spline construction using
267  // the indices above to avoid computing them for each j
268  for (decltype(n) j = 0; j < n; ++j)
271 
272  // Construct single row spline; get back the second derivatives wrt x2 coord
273  // on the x2 grid points
274  spline(_x2, column_spline_eval, row_spline_second_derivs, yx11, yx1n);
275 }
276 
277 void
279  std::vector<Real> & row_spline_eval,
280  std::vector<Real> & column_spline_second_derivs,
281  Real yx21 /*= _deriv_bound*/,
282  Real yx2n /*= _deriv_bound*/)
283 {
284  auto m = _x1.size();
285 
286  // Find the indices that bound the point x2
287  unsigned int klo, khi;
288  findInterval(_x2, x2, klo, khi);
289 
290  // Evaluate m row-splines to get y-values for column spline construction using
291  // the indices above to avoid computing them for each j
292  for (decltype(m) j = 0; j < m; ++j)
293  _row_spline_eval[j] = SplineInterpolationBase::sample(_x2, _y[j], _y2_rows[j], x2, klo, khi);
294 
295  // Construct single column spline; get back the second derivatives wrt x1 coord
296  // on the x1 grid points
297  spline(_x1, row_spline_eval, column_spline_second_derivs, yx21, yx2n);
298 }
Real sample(Real x1, Real x2, Real yx11=_deriv_bound, Real yx1n=_deriv_bound)
Samples value at point (x1, x2)
std::vector< Real > _x1
Independent values in the x1 direction.
std::vector< std::vector< Real > > _y_trans
Transpose of _y.
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
Real sample2ndDerivative(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
void errorCheck()
Sanity checks on input data.
std::vector< std::vector< Real > > _y
The dependent values at (x1, x2) points.
Real sample(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
Real sampleDerivative(Real x1, Real x2, unsigned int deriv_var, Real yp1=_deriv_bound, Real ypn=_deriv_bound)
Samples first derivative at point (x1, x2)
Real sampleDerivative(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
std::vector< Real > _x2
Independent values in the x2 direction.
void findInterval(const std::vector< Real > &x, Real x_int, unsigned int &klo, unsigned int &khi) const
std::vector< Real > _yx11
Boundary conditions.
std::vector< std::vector< Real > > _y2_rows
Second derivatives.
void constructColumnSpline(Real x2, std::vector< Real > &spline_eval, std::vector< Real > &spline_second_derivs, Real yx21=_deriv_bound, Real yx2n=_deriv_bound)
Helper functions to evaluate row splines and construct column spline for the given point...
void constructRowSplineSecondDerivativeTable()
Precompute tables of row (column) spline second derivatives and store them to reduce computational de...
void solve()
Calculates the tables of second derivatives.
std::vector< Real > _row_spline_second_derivs
Vectors used during sampling.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static int _file_number
File number for data dump.
void constructRowSpline(Real x1, std::vector< Real > &spline_eval, std::vector< Real > &spline_second_derivs, Real yx11=_deriv_bound, Real yx1n=_deriv_bound)
Helper functions to evaluate column splines and construct row spline for the given point...
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.
Real sample2ndDerivative(Real x1, Real x2, unsigned int deriv_var, Real yp1=_deriv_bound, Real ypn=_deriv_bound)
Samples second derivative at point (x1, x2)
void sampleValueAndDerivatives(Real x1, Real x2, Real &y, Real &dy1, Real &dy2, Real yx11=_deriv_bound, Real yx1n=_deriv_bound, Real yx21=_deriv_bound, Real yx2n=_deriv_bound)
Samples value and first derivatives at point (x1, x2) Use this function for speed when computing both...
std::vector< std::vector< Real > > _y2_columns
std::vector< Real > _column_spline_second_derivs
void setData(const std::vector< Real > &x1, const std::vector< Real > &x2, const std::vector< std::vector< Real >> &y, const std::vector< Real > &yx11=std::vector< Real >(), const std::vector< Real > &yx1n=std::vector< Real >(), const std::vector< Real > &yx21=std::vector< Real >(), const std::vector< Real > &yx2n=std::vector< Real >())
Set the x1, x2 and y values, and first derivatives at the edges.