www.mooseframework.org
BicubicSplineInterpolation.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
16 #include "MooseError.h"
17 
19 
21 
23  const std::vector<Real> & x2,
24  const std::vector<std::vector<Real>> & y,
25  const std::vector<Real> & yx11,
26  const std::vector<Real> & yx1n,
27  const std::vector<Real> & yx21,
28  const std::vector<Real> & yx2n)
30  _x1(x1),
31  _x2(x2),
32  _y(y),
33  _yx11(yx11),
34  _yx1n(yx1n),
35  _yx21(yx21),
36  _yx2n(yx2n)
37 {
38  errorCheck();
39  solve();
40 }
41 
42 void
43 BicubicSplineInterpolation::setData(const std::vector<Real> & x1,
44  const std::vector<Real> & x2,
45  const std::vector<std::vector<Real>> & y,
46  const std::vector<Real> & yx11,
47  const std::vector<Real> & yx1n,
48  const std::vector<Real> & yx21,
49  const std::vector<Real> & yx2n)
50 {
51  _x1 = x1;
52  _x2 = x2;
53  _y = y;
54  _yx11 = yx11;
55  _yx1n = yx1n;
56  _yx21 = yx21;
57  _yx2n = yx2n;
58  errorCheck();
59  solve();
60 }
61 
62 void
64 {
65  auto m = _x1.size(), n = _x2.size();
66 
67  if (_y.size() != m)
68  mooseError("y row dimension does not match the size of x1.");
69  else
70  for (decltype(m) i = 0; i < _y.size(); ++i)
71  if (_y[i].size() != n)
72  mooseError("y column dimension does not match the size of x2.");
73 
74  if (_yx11.empty())
75  _yx11.resize(n, _deriv_bound);
76  else if (_yx11.size() != n)
77  mooseError("The length of the vectors holding the first derivatives of y with respect to x1 "
78  "must match the length of x2.");
79 
80  if (_yx1n.empty())
81  _yx1n.resize(n, _deriv_bound);
82  else if (_yx1n.size() != n)
83  mooseError("The length of the vectors holding the first derivatives of y with respect to x1 "
84  "must match the length of x2.");
85 
86  if (_yx21.empty())
87  _yx21.resize(m, _deriv_bound);
88  else if (_yx21.size() != m)
89  mooseError("The length of the vectors holding the first derivatives of y with respect to x2 "
90  "must match the length of x1.");
91 
92  if (_yx2n.empty())
93  _yx2n.resize(m, _deriv_bound);
94  else if (_yx2n.size() != m)
95  mooseError("The length of the vectors holding the first derivatives of y with respect to x2 "
96  "must match the length of x1.");
97 }
98 
99 void
101 {
102  auto m = _x1.size();
103  _y2_rows.resize(m);
104 
105  if (_yx21.empty())
106  for (decltype(m) j = 0; j < m; ++j)
107  spline(_x2, _y[j], _y2_rows[j]);
108 
109  else
110  for (decltype(m) j = 0; j < m; ++j)
111  spline(_x2, _y[j], _y2_rows[j], _yx21[j], _yx2n[j]);
112 }
113 
114 void
116 {
117  auto m = _x1.size();
118  auto n = _x2.size();
119  _y2_columns.resize(n);
120  _y_trans.resize(n);
121 
122  for (decltype(n) j = 0; j < n; ++j)
123  _y_trans[j].resize(m);
124 
125  // transpose the _y values so the columns can be easily iterated over
126  for (decltype(n) i = 0; i < _y.size(); ++i)
127  for (decltype(n) j = 0; j < _y[0].size(); ++j)
128  _y_trans[j][i] = _y[i][j];
129 
130  if (_yx11.empty())
131  for (decltype(n) j = 0; j < n; ++j)
132  spline(_x1, _y_trans[j], _y2_columns[j]);
133 
134  else
135  for (decltype(n) j = 0; j < n; ++j)
136  spline(_x1, _y_trans[j], _y2_columns[j], _yx11[j], _yx1n[j]);
137 }
138 
139 void
141 {
144 }
145 
146 Real
148  Real x2,
149  Real yx11 /* = _deriv_bound*/,
150  Real yx1n /* = _deriv_bound*/)
151 {
152  auto m = _x1.size();
153  std::vector<Real> column_spline_second_derivs(m), row_spline_eval(m);
154 
155  // Evaluate m row-splines to get y-values for column spline construction
156  for (decltype(m) j = 0; j < m; ++j)
157  row_spline_eval[j] = SplineInterpolationBase::sample(_x2, _y[j], _y2_rows[j], x2);
158 
159  // Construct single column spline; get back the second derivatives wrt x1 coord on the x1 grid
160  // points
161  spline(_x1, row_spline_eval, column_spline_second_derivs, yx11, yx1n);
162 
163  // Evaluate newly constructed column spline
164  return SplineInterpolationBase::sample(_x1, row_spline_eval, column_spline_second_derivs, x1);
165 }
166 
167 Real
169  Real x2,
170  unsigned int deriv_var,
171  Real yp1 /* = _deriv_bound*/,
172  Real ypn /* = _deriv_bound*/)
173 {
174  auto m = _x1.size();
175 
176  // Take derivative along x1 axis
177  if (deriv_var == 1)
178  {
179  std::vector<Real> column_spline_second_derivs(m), row_spline_eval(m);
180 
181  // Evaluate m row-splines to get y-values for column spline construction
182  for (decltype(m) j = 0; j < m; ++j)
183  row_spline_eval[j] = SplineInterpolationBase::sample(_x2, _y[j], _y2_rows[j], x2);
184 
185  // Construct single column spline; get back the second derivatives wrt x1 coord on the x1 grid
186  // points
187  spline(_x1, row_spline_eval, column_spline_second_derivs, yp1, ypn);
188 
189  // Evaluate derivative wrt x1 of newly constructed column spline
191  _x1, row_spline_eval, column_spline_second_derivs, x1);
192  }
193 
194  // Take derivative along x2 axis
195  else if (deriv_var == 2)
196  {
197  auto n = _x2.size();
198  std::vector<Real> row_spline_second_derivs(n), column_spline_eval(n);
199 
200  // Evaluate n column-splines to get y-values for row spline construction
201  for (decltype(n) j = 0; j < n; ++j)
202  column_spline_eval[j] = SplineInterpolationBase::sample(_x1, _y_trans[j], _y2_columns[j], x1);
203 
204  // Construct single row spline; get back the second derivatives wrt x2 coord on the x2 grid
205  // points
206  spline(_x2, column_spline_eval, row_spline_second_derivs, yp1, ypn);
207 
208  // Evaluate derivative wrt x2 of newly constructed row spline
210  _x2, column_spline_eval, row_spline_second_derivs, x2);
211  }
212 
213  else
214  {
215  mooseError("deriv_var must be either 1 or 2");
216  return 0;
217  }
218 }
219 
220 Real
222  Real x2,
223  unsigned int deriv_var,
224  Real yp1 /* = _deriv_bound*/,
225  Real ypn /* = _deriv_bound*/)
226 {
227  auto m = _x1.size();
228 
229  // Take second derivative along x1 axis
230  if (deriv_var == 1)
231  {
232  std::vector<Real> column_spline_second_derivs(m), row_spline_eval(m);
233 
234  // Evaluate m row-splines to get y-values for column spline construction
235  for (decltype(m) j = 0; j < m; ++j)
236  row_spline_eval[j] = SplineInterpolationBase::sample(_x2, _y[j], _y2_rows[j], x2);
237 
238  // Construct single column spline; get back the second derivatives wrt x1 coord on the x1 grid
239  // points
240  spline(_x1, row_spline_eval, column_spline_second_derivs, yp1, ypn);
241 
242  // Evaluate second derivative wrt x1 of newly constructed column spline
244  _x1, row_spline_eval, column_spline_second_derivs, x1);
245  }
246 
247  // Take second derivative along x2 axis
248  else if (deriv_var == 2)
249  {
250  auto n = _x2.size();
251  std::vector<Real> row_spline_second_derivs(n), column_spline_eval(n);
252 
253  // Evaluate n column-splines to get y-values for row spline construction
254  for (decltype(n) j = 0; j < n; ++j)
255  column_spline_eval[j] = SplineInterpolationBase::sample(_x1, _y_trans[j], _y2_columns[j], x1);
256 
257  // Construct single row spline; get back the second derivatives wrt x2 coord on the x2 grid
258  // points
259  spline(_x2, column_spline_eval, row_spline_second_derivs, yp1, ypn);
260 
261  // Evaluate second derivative wrt x2 of newly constructed row spline
263  _x2, column_spline_eval, row_spline_second_derivs, x2);
264  }
265 
266  else
267  {
268  mooseError("deriv_var must be either 1 or 2");
269  return 0;
270  }
271 }
Real sample(Real x1, Real x2, Real yx11=_deriv_bound, Real yx1n=_deriv_bound)
Real sample2ndDerivative(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
Real sample(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
std::vector< std::vector< Real > > _y_trans
Transpose of _y.
std::vector< std::vector< Real > > _y
void mooseError(Args &&...args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:182
Real sampleDerivative(Real x1, Real x2, unsigned int deriv_var, Real yp1=_deriv_bound, Real ypn=_deriv_bound)
std::vector< Real > _yx11
Boundary conditions.
std::vector< std::vector< Real > > _y2_rows
Second derivatives.
PetscInt m
PetscInt n
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)
std::vector< std::vector< Real > > _y2_columns
Real sampleDerivative(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &y2, Real x_int) const
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, y- values and first derivatives.