www.mooseframework.org
BicubicSplineInterpolation.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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();
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();
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:207
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...
PetscInt m
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.
static int _file_number
File number for data dump.
PetscInt n
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.