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 "BicubicSplineFunction.h"
11 :
12 : registerMooseObject("MooseApp", BicubicSplineFunction);
13 :
14 : InputParameters
15 14346 : BicubicSplineFunction::validParams()
16 : {
17 14346 : InputParameters params = Function::validParams();
18 14346 : params.addClassDescription(
19 : "Define a bicubic spline function from interpolated data defined by input parameters.");
20 :
21 14346 : MooseEnum normal_component("x=0 y=1 z=2", "z");
22 14346 : params.addParam<MooseEnum>(
23 : "normal_component",
24 : normal_component,
25 : "The component of the geometry that is normal to the spline x1/x2 values");
26 14346 : params.addRequiredParam<std::vector<Real>>("x1", "The first independent coordinate.");
27 14346 : params.addRequiredParam<std::vector<Real>>("x2", "The second independent coordinate.");
28 14346 : params.addRequiredParam<std::vector<Real>>("y", "The dependent values");
29 14346 : params.addParam<std::vector<Real>>(
30 : "yx11", "The values of the derivative wrt x1 on the lower interpolation grid points.");
31 14346 : params.addParam<std::vector<Real>>(
32 : "yx1n", "The values of the derivative wrt x1 on the upper interpolation grid points.");
33 14346 : params.addParam<std::vector<Real>>(
34 : "yx21", "The values of the derivative wrt x2 on the lower interpolation grid points.");
35 14346 : params.addParam<std::vector<Real>>(
36 : "yx2n", "The values of the derivative wrt x2 on the upper interpolation grid points.");
37 14346 : params.addParam<FunctionName>(
38 : "yx1", "1e30", "The functional form of the derivative with respect to x1.");
39 14346 : params.addParam<FunctionName>(
40 : "yx2", "1e30", "The functional form of the derivative with respect to x2.");
41 :
42 28692 : return params;
43 14346 : }
44 :
45 42 : BicubicSplineFunction::BicubicSplineFunction(const InputParameters & parameters)
46 : : Function(parameters),
47 : FunctionInterface(this),
48 42 : _normal_component(getParam<MooseEnum>("normal_component")),
49 42 : _yx1(getFunction("yx1")),
50 126 : _yx2(getFunction("yx2"))
51 : {
52 42 : _x1 = getParam<std::vector<Real>>("x1");
53 42 : _x2 = getParam<std::vector<Real>>("x2");
54 42 : std::vector<Real> yvec = getParam<std::vector<Real>>("y");
55 42 : if (isParamValid("yx11"))
56 42 : _yx11 = getParam<std::vector<Real>>("yx11");
57 42 : if (isParamValid("yx1n"))
58 42 : _yx1n = getParam<std::vector<Real>>("yx1n");
59 42 : if (isParamValid("yx21"))
60 42 : _yx21 = getParam<std::vector<Real>>("yx21");
61 42 : if (isParamValid("yx2n"))
62 42 : _yx2n = getParam<std::vector<Real>>("yx2n");
63 :
64 42 : unsigned int m = _x1.size(), n = _x2.size(), mn = yvec.size();
65 42 : if (m * n != mn)
66 0 : mooseError("The length of the supplied y must be equal to the lengths of x1 and x2 multiplied "
67 : "together");
68 :
69 42 : std::vector<std::vector<Real>> y(m, std::vector<Real>(n));
70 42 : unsigned int k = 0;
71 168 : for (unsigned int i = 0; i < m; ++i)
72 630 : for (unsigned int j = 0; j < n; ++j)
73 504 : y[i][j] = yvec[k++];
74 :
75 42 : if (_yx11.empty())
76 0 : _yx11.resize(n, 1e30);
77 42 : else if (_yx11.size() != n)
78 0 : mooseError("The length of the vectors holding the first derivatives of y with respect to x1 "
79 : "must match the length of x2.");
80 :
81 42 : if (_yx1n.empty())
82 0 : _yx1n.resize(n, 1e30);
83 42 : else if (_yx1n.size() != n)
84 0 : mooseError("The length of the vectors holding the first derivatives of y with respect to x1 "
85 : "must match the length of x2.");
86 :
87 42 : if (_yx21.empty())
88 0 : _yx21.resize(m, 1e30);
89 42 : else if (_yx21.size() != m)
90 0 : mooseError("The length of the vectors holding the first derivatives of y with respect to x2 "
91 : "must match the length of x1.");
92 :
93 42 : if (_yx2n.empty())
94 0 : _yx2n.resize(m, 1e30);
95 42 : else if (_yx2n.size() != m)
96 0 : mooseError("The length of the vectors holding the first derivatives of y with respect to x2 "
97 : "must match the length of x1.");
98 :
99 42 : _ipol.setData(_x1, _x2, y, _yx11, _yx1n, _yx21, _yx2n);
100 :
101 42 : if (_normal_component == 0)
102 : {
103 : // YZ plane
104 14 : _x1_index = 1;
105 14 : _x2_index = 2;
106 : }
107 28 : else if (_normal_component == 1)
108 : {
109 : // ZX plane
110 14 : _x1_index = 2;
111 14 : _x2_index = 0;
112 : }
113 : else
114 : {
115 : // XY plane
116 14 : _x1_index = 0;
117 14 : _x2_index = 1;
118 : }
119 42 : }
120 :
121 : Real
122 7450 : BicubicSplineFunction::value(Real /*t*/, const Point & p) const
123 : {
124 : // Call yx11/yx1n with the correctly oriented points
125 7450 : Real x1_begin = _x1[0];
126 7450 : Real x1_end = p(_x2_index);
127 7450 : Real xn_begin = _x1.back();
128 7450 : Real xn_end = p(_x2_index);
129 :
130 7450 : Point x1(0, 0, 0);
131 7450 : Point xn(0, 0, 0);
132 :
133 7450 : x1(_x1_index) = x1_begin;
134 7450 : x1(_x2_index) = x1_end;
135 7450 : xn(_x1_index) = xn_begin;
136 7450 : xn(_x2_index) = xn_end;
137 :
138 7450 : Real yx11 = _yx1.value(0, x1);
139 7450 : Real yx1n = _yx1.value(0, xn);
140 :
141 14900 : return _ipol.sample(p(_x1_index), p(_x2_index), yx11, yx1n);
142 : }
143 :
144 : Real
145 19000 : BicubicSplineFunction::derivative(const Point & p, unsigned int deriv_var) const
146 : {
147 : Real yp1, ypn;
148 19000 : Point x1(0, 0, 0);
149 19000 : Point xn(0, 0, 0);
150 19000 : if (deriv_var == 1)
151 : {
152 : // Call yx11/yx1n with the correctly oriented points
153 9500 : Real x1_begin = _x1[0];
154 9500 : Real x1_end = p(_x2_index);
155 9500 : Real xn_begin = _x1.back();
156 9500 : Real xn_end = p(_x2_index);
157 :
158 9500 : x1(_x1_index) = x1_begin;
159 9500 : x1(_x2_index) = x1_end;
160 9500 : xn(_x1_index) = xn_begin;
161 9500 : xn(_x2_index) = xn_end;
162 :
163 9500 : yp1 = _yx1.value(0, x1);
164 9500 : ypn = _yx1.value(0, xn);
165 : }
166 9500 : else if (deriv_var == 2)
167 : {
168 : // Call yx11/yx1n with the correctly oriented points
169 9500 : Real x1_begin = p(_x1_index);
170 9500 : Real x1_end = _x2[0];
171 9500 : Real xn_begin = p(_x1_index);
172 9500 : Real xn_end = _x2.back();
173 :
174 9500 : x1(_x1_index) = x1_begin;
175 9500 : x1(_x2_index) = x1_end;
176 9500 : xn(_x1_index) = xn_begin;
177 9500 : xn(_x2_index) = xn_end;
178 :
179 9500 : yp1 = _yx2.value(0, x1);
180 9500 : ypn = _yx2.value(0, xn);
181 : }
182 : else
183 0 : mooseError("deriv_var must equal 1 or 2");
184 :
185 38000 : return _ipol.sampleDerivative(p(_x1_index), p(_x2_index), deriv_var, yp1, ypn);
186 : }
187 :
188 : RealGradient
189 9500 : BicubicSplineFunction::gradient(Real /*t*/, const Point & p) const
190 : {
191 9500 : RealGradient grad = RealGradient(0, 0, 0);
192 :
193 9500 : Real dF_dx1 = derivative(p, 1);
194 9500 : Real dF_dx2 = derivative(p, 2);
195 :
196 9500 : grad(_x1_index) = dF_dx1;
197 9500 : grad(_x2_index) = dF_dx2;
198 :
199 9500 : return grad;
200 : }
201 :
202 : Real
203 0 : BicubicSplineFunction::secondDerivative(const Point & p, unsigned int deriv_var) const
204 : {
205 : Real yp1, ypn;
206 0 : Point x1(0, 0, 0);
207 0 : Point xn(0, 0, 0);
208 0 : if (deriv_var == 1)
209 : {
210 : // Call yx11/yx1n with the correctly oriented points
211 0 : Real x1_begin = _x1[0];
212 0 : Real x1_end = p(_x2_index);
213 0 : Real xn_begin = _x1.back();
214 0 : Real xn_end = p(_x2_index);
215 :
216 0 : x1(_x1_index) = x1_begin;
217 0 : x1(_x2_index) = x1_end;
218 0 : xn(_x1_index) = xn_begin;
219 0 : xn(_x2_index) = xn_end;
220 :
221 0 : yp1 = _yx1.value(0, x1);
222 0 : ypn = _yx1.value(0, xn);
223 : }
224 0 : else if (deriv_var == 2)
225 : {
226 : // Call yx11/yx1n with the correctly oriented points
227 0 : Real x1_begin = p(_x1_index);
228 0 : Real x1_end = _x2[0];
229 0 : Real xn_begin = p(_x1_index);
230 0 : Real xn_end = _x2.back();
231 :
232 0 : x1(_x1_index) = x1_begin;
233 0 : x1(_x2_index) = x1_end;
234 0 : xn(_x1_index) = xn_begin;
235 0 : xn(_x2_index) = xn_end;
236 :
237 0 : yp1 = _yx2.value(0, x1);
238 0 : ypn = _yx2.value(0, xn);
239 : }
240 : else
241 0 : mooseError("deriv_var must equal 1 or 2");
242 :
243 0 : return _ipol.sample2ndDerivative(p(_x1_index), p(_x2_index), deriv_var, yp1, ypn);
244 : }
|
