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 "MooseTypes.h"
13 : #include "BidimensionalInterpolation.h"
14 :
15 : /**
16 : * This class interpolates tabulated data with a bicubic function. In order to
17 : * minimize the computational expense of each sample, the coefficients at each
18 : * point in the tabulated data are computed once in advance, and then accessed
19 : * during the interpolation.
20 : *
21 : * As a result, this bicubic interpolation can be much faster than the bicubic
22 : * spline interpolation method in BicubicSplineInterpolation.
23 : *
24 : * Adapted from Numerical Recipes in C (section 3.6). The terminology used is
25 : * consistent with that used in Numerical Recipes, where moving over a column
26 : * corresponds to moving over the x1 coord. Likewise, moving over a row means
27 : * moving over the x2 coord.
28 : */
29 : class BicubicInterpolation : public BidimensionalInterpolation
30 : {
31 : public:
32 : BicubicInterpolation(const std::vector<Real> & x1,
33 : const std::vector<Real> & x2,
34 : const std::vector<std::vector<Real>> & y);
35 :
36 1 : virtual ~BicubicInterpolation() = default;
37 :
38 : /**
39 : * Sanity checks on input data
40 : */
41 : void errorCheck();
42 :
43 : /**
44 : * Samples value at point (x1, x2)
45 : */
46 : Real sample(const Real x1, const Real x2) const override;
47 : ADReal sample(const ADReal & x1, const ADReal & x2) const override;
48 : ChainedReal sample(const ChainedReal & x1, const ChainedReal & x2) const override;
49 :
50 : /**
51 : * Samples value and first derivatives at point (x1, x2)
52 : * Use this function for speed when computing both value and derivatives,
53 : * as it minimizes the amount of time spent locating the point in the
54 : * tabulated data
55 : */
56 : using BidimensionalInterpolation::sampleValueAndDerivatives;
57 : virtual void
58 : sampleValueAndDerivatives(Real x1, Real x2, Real & y, Real & dy1, Real & dy2) const override;
59 : virtual void sampleValueAndDerivatives(
60 : const ADReal & x1, const ADReal & x2, ADReal & y, ADReal & dy1, ADReal & dy2) const override;
61 : virtual void sampleValueAndDerivatives(const ChainedReal & x1,
62 : const ChainedReal & x2,
63 : ChainedReal & y,
64 : ChainedReal & dy1,
65 : ChainedReal & dy2) const override;
66 :
67 : /**
68 : * Samples first derivative at point (x1, x2)
69 : */
70 : using BidimensionalInterpolation::sampleDerivative;
71 : Real sampleDerivative(Real x1, Real x2, unsigned int deriv_var) const override;
72 :
73 : /**
74 : * Samples second derivative at point (x1, x2)
75 : */
76 : Real sample2ndDerivative(Real x1, Real x2, unsigned int deriv_var) const override;
77 :
78 : /**
79 : * Precompute all of the coefficients for the bicubic interpolation to avoid
80 : * calculating them repeatedly
81 : */
82 : void precomputeCoefficients();
83 :
84 : protected:
85 : /**
86 : * Find the indices of the dependent values axis which bracket the point xi
87 : */
88 : template <typename T>
89 : void findInterval(const std::vector<Real> & x,
90 : const T & xi,
91 : unsigned int & klo,
92 : unsigned int & khi,
93 : T & xs) const;
94 :
95 : template <typename T>
96 : T sampleInternal(const T & x1, const T & x2) const;
97 :
98 : template <typename T>
99 : void sampleValueAndDerivativesInternal(T x1, T x2, T & y, T & dy1, T & dy2) const;
100 :
101 : /**
102 : * Provides the values of the first derivatives in each direction at all
103 : * points in the table of data, as well as the mixed second derivative.
104 : * This is implemented using finite differencing, but could be supplied through
105 : * other means (such as by sampling with cubic splines)
106 : */
107 : void tableDerivatives(std::vector<std::vector<Real>> & dy_dx1,
108 : std::vector<std::vector<Real>> & dy_dx2,
109 : std::vector<std::vector<Real>> & d2y_dx1x2);
110 :
111 : /// The dependent values at (x1, x2) points
112 : std::vector<std::vector<Real>> _y;
113 :
114 : /// Matrix of precomputed coefficients
115 : /// There are four coefficients in each direction at each dependent value
116 : std::vector<std::vector<std::vector<std::vector<Real>>>> _bicubic_coeffs;
117 :
118 : /// Matrix used to calculate bicubic interpolation coefficients
119 : /// (from Numerical Recipes)
120 : const std::vector<std::vector<int>> _wt{
121 : {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
122 : {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
123 : {-3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1, 0, 0, 0, 0},
124 : {2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0},
125 : {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
126 : {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
127 : {0, 0, 0, 0, -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1},
128 : {0, 0, 0, 0, 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1},
129 : {-3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
130 : {0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0},
131 : {9, -9, 9, -9, 6, 3, -3, -6, 6, -6, -3, 3, 4, 2, 1, 2},
132 : {-6, 6, -6, 6, -4, -2, 2, 4, -3, 3, 3, -3, -2, -1, -1, -2},
133 : {2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
134 : {0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0},
135 : {-6, 6, -6, 6, -3, -3, 3, 3, -4, 4, 2, -2, -2, -2, -1, -1},
136 : {4, -4, 4, -4, 2, 2, -2, -2, 2, -2, -2, 2, 1, 1, 1, 1}};
137 : };
|
