BicubicInterpolation Class Reference

This class interpolates tabulated data with a bicubic function. More...

#include <BicubicInterpolation.h>

Public Member Functions
	BicubicInterpolation (const std::vector< Real > &x1, const std::vector< Real > &x2, const std::vector< std::vector< Real >> &y)
virtual	~BicubicInterpolation ()=default
void	errorCheck ()
	Sanity checks on input data. More...
Real	sample (Real x1, Real x2)
	Samples value at point (x1, x2) More...
void	sampleValueAndDerivatives (Real x1, Real x2, Real &y, Real &dy1, Real &dy2)
	Samples value and first derivatives at point (x1, x2) Use this function for speed when computing both value and derivatives, as it minimizes the amount of time spent locating the point in the tabulated data. More...
Real	sampleDerivative (Real x1, Real x2, unsigned int deriv_var)
	Samples first derivative at point (x1, x2) More...
Real	sample2ndDerivative (Real x1, Real x2, unsigned int deriv_var)
	Samples second derivative at point (x1, x2) More...
void	precomputeCoefficients ()
	Precompute all of the coefficients for the bicubic interpolation to avoid calculating them repeatedly. More...

Protected Member Functions
void	findInterval (const std::vector< Real > &x, Real xi, unsigned int &klo, unsigned int &khi, Real &xs)
	Find the indices of the dependent values axis which bracket the point xi. More...
void	tableDerivatives (std::vector< std::vector< Real >> &dy_dx1, std::vector< std::vector< Real >> &dy_dx2, std::vector< std::vector< Real >> &d2y_dx1x2)
	Provides the values of the first derivatives in each direction at all points in the table of data, as well as the mixed second derivative. More...

Protected Attributes
std::vector< Real >	_x1
	Independent values in the x1 direction. More...
std::vector< Real >	_x2
	Independent values in the x2 direction. More...
std::vector< std::vector< Real > >	_y
	The dependent values at (x1, x2) points. More...
std::vector< std::vector< std::vector< std::vector< Real > > > >	_bicubic_coeffs
	Matrix of precomputed coefficients There are four coefficients in each direction at each dependent value. More...
const std::vector< std::vector< int > >	_wt
	Matrix used to calculate bicubic interpolation coefficients (from Numerical Recipes) More...

Detailed Description

This class interpolates tabulated data with a bicubic function.

In order to minimize the computational expense of each sample, the coefficients at each point in the tabulated data are computed once in advance, and then accessed during the interpolation.

As a result, this bicubic interpolation can be much faster than the bicubic spline interpolation method in BicubicSplineInterpolation.

Adapted from Numerical Recipes in C (section 3.6). The terminology used is consistent with that used in Numerical Recipes, where moving over a column corresponds to moving over the x1 coord. Likewise, moving over a row means moving over the x2 coord.

Definition at line 29 of file BicubicInterpolation.h.

Constructor & Destructor Documentation

BicubicInterpolation()

BicubicInterpolation::BicubicInterpolation	(	const std::vector< Real > &	x1,
		const std::vector< Real > &	x2,
		const std::vector< std::vector< Real >> &	y
	)

Definition at line 15 of file BicubicInterpolation.C.

  : _x1(x1), _x2(x2), _y(y)
{
  errorCheck();

  // Resize the vector of coefficients
  auto m = _x1.size();
  auto n = _x2.size();

  _bicubic_coeffs.resize(m - 1);
  for (auto i = beginIndex(_bicubic_coeffs); i < _bicubic_coeffs.size(); ++i)
  {
    _bicubic_coeffs[i].resize(n - 1);
    for (auto j = beginIndex(_bicubic_coeffs[i]); j < _bicubic_coeffs[i].size(); ++j)
    {
      _bicubic_coeffs[i][j].resize(4);
      for (unsigned int k = 0; k < 4; ++k)
        _bicubic_coeffs[i][j][k].resize(4);
    }
  }

  // Precompute the coefficients
  precomputeCoefficients();
}

~BicubicInterpolation()

virtual BicubicInterpolation::~BicubicInterpolation ( )

virtualdefault

Member Function Documentation

errorCheck()

void BicubicInterpolation::errorCheck

(

)

Sanity checks on input data.

Definition at line 347 of file BicubicInterpolation.C.

Referenced by BicubicInterpolation().

{
  auto m = _x1.size(), n = _x2.size();

  if (_y.size() != m)
    mooseError("y row dimension does not match the size of x1 in BicubicInterpolation");
  else
    for (auto i = beginIndex(_y); i < _y.size(); ++i)
      if (_y[i].size() != n)
        mooseError("y column dimension does not match the size of x2 in BicubicInterpolation");
}

findInterval()

void BicubicInterpolation::findInterval ( const std::vector< Real > &  x, Real  xi, unsigned int &  klo, unsigned int &  khi, Real &  xs  )

protected

Find the indices of the dependent values axis which bracket the point xi.

Definition at line 320 of file BicubicInterpolation.C.

Referenced by sample(), sample2ndDerivative(), sampleDerivative(), and sampleValueAndDerivatives().

{
  // Find the indices that bracket the point xi
  klo = 0;
  mooseAssert(x.size() >= 2,
              "There must be at least two points in the table in BicubicInterpolation");

  khi = x.size() - 1;
  while (khi - klo > 1)
  {
    unsigned int kmid = (khi + klo) / 2;
    if (x[kmid] > xi)
      khi = kmid;
    else
      klo = kmid;
  }

  // Now find the scaled position, normalized to [0,1]
  Real d = x[khi] - x[klo];
  xs = xi - x[klo];

  if (!MooseUtils::absoluteFuzzyEqual(d, 0.0))
    xs /= d;
}

precomputeCoefficients()

void BicubicInterpolation::precomputeCoefficients

(

)

Precompute all of the coefficients for the bicubic interpolation to avoid calculating them repeatedly.

Definition at line 191 of file BicubicInterpolation.C.

Referenced by BicubicInterpolation().

{
  // Calculate the first derivatives in each direction at each point, and the
  // mixed second derivative
  std::vector<std::vector<Real>> dy_dx1, dy_dx2, d2y_dx1x2;
  tableDerivatives(dy_dx1, dy_dx2, d2y_dx1x2);

  // Now solve for the coefficients at each point in the grid
  for (auto i = beginIndex(_bicubic_coeffs); i < _bicubic_coeffs.size(); ++i)
    for (auto j = beginIndex(_bicubic_coeffs); j < _bicubic_coeffs[i].size(); ++j)
    {
      // Distance between corner points in each direction
      const Real d1 = _x1[i + 1] - _x1[i];
      const Real d2 = _x2[j + 1] - _x2[j];

      // Values of function and derivatives of the four corner points
      std::vector<Real> y = {_y[i][j], _y[i + 1][j], _y[i + 1][j + 1], _y[i][j + 1]};
      std::vector<Real> y1 = {
          dy_dx1[i][j], dy_dx1[i + 1][j], dy_dx1[i + 1][j + 1], dy_dx1[i][j + 1]};
      std::vector<Real> y2 = {
          dy_dx2[i][j], dy_dx2[i + 1][j], dy_dx2[i + 1][j + 1], dy_dx2[i][j + 1]};
      std::vector<Real> y12 = {
          d2y_dx1x2[i][j], d2y_dx1x2[i + 1][j], d2y_dx1x2[i + 1][j + 1], d2y_dx1x2[i][j + 1]};

      std::vector<Real> cl(16), x(16);
      Real xx;
      unsigned int count = 0;

      // Temporary vector used in the matrix multiplication
      for (unsigned int k = 0; k < 4; ++k)
      {
        x[k] = y[k];
        x[k + 4] = y1[k] * d1;
        x[k + 8] = y2[k] * d2;
        x[k + 12] = y12[k] * d1 * d2;
      }

      // Multiply by the matrix of constants
      for (unsigned int k = 0; k < 16; ++k)
      {
        xx = 0.0;
        for (unsigned int q = 0; q < 16; ++q)
          xx += _wt[k][q] * x[q];

        cl[k] = xx;
      }

      // Unpack results into coefficient table
      for (unsigned int k = 0; k < 4; ++k)
        for (unsigned int q = 0; q < 4; ++q)
        {
          _bicubic_coeffs[i][j][k][q] = cl[count];
          count++;
        }
    }
}

sample()

Real BicubicInterpolation::sample	(	Real	x1,
		Real	x2
	)

Samples value at point (x1, x2)

Definition at line 43 of file BicubicInterpolation.C.

{
  unsigned int x1l, x1u, x2l, x2u;
  Real t, u;
  findInterval(_x1, x1, x1l, x1u, t);
  findInterval(_x2, x2, x2l, x2u, u);

  Real sample = 0.0;
  for (unsigned int i = 0; i < 4; ++i)
    for (unsigned int j = 0; j < 4; ++j)
      sample += _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) * MathUtils::pow(u, j);

  return sample;
}

sample2ndDerivative()

Real BicubicInterpolation::sample2ndDerivative	(	Real	x1,
		Real	x2,
		unsigned int	deriv_var
	)

Samples second derivative at point (x1, x2)

Definition at line 107 of file BicubicInterpolation.C.

{
  unsigned int x1l, x1u, x2l, x2u;
  Real t, u;
  findInterval(_x1, x1, x1l, x1u, t);
  findInterval(_x2, x2, x2l, x2u, u);

  // Take derivative along x1 axis
  // Note: sum from i = 2 as the first two terms are zero
  if (deriv_var == 1)
  {
    Real sample_deriv = 0.0;
    for (unsigned int i = 2; i < 4; ++i)
      for (unsigned int j = 0; j < 4; ++j)
        sample_deriv += i * (i - 1) * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i - 2) *
                        MathUtils::pow(u, j);

    Real d = _x1[x1u] - _x1[x1l];

    if (!MooseUtils::absoluteFuzzyEqual(d, 0.0))
      sample_deriv /= (d * d);

    return sample_deriv;
  }

  // Take derivative along x2 axis
  // Note: sum from j = 2 as the first two terms are zero
  else if (deriv_var == 2)
  {
    Real sample_deriv = 0.0;
    for (unsigned int i = 0; i < 4; ++i)
      for (unsigned int j = 2; j < 4; ++j)
        sample_deriv += j * (j - 1) * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) *
                        MathUtils::pow(u, j - 2);

    Real d = _x2[x2u] - _x2[x2l];

    if (!MooseUtils::absoluteFuzzyEqual(d, 0.0))
      sample_deriv /= (d * d);

    return sample_deriv;
  }

  else
    mooseError("deriv_var must be either 1 or 2 in BicubicInterpolation");
}

sampleDerivative()

Real BicubicInterpolation::sampleDerivative	(	Real	x1,
		Real	x2,
		unsigned int	deriv_var
	)

Samples first derivative at point (x1, x2)

Definition at line 59 of file BicubicInterpolation.C.

{
  unsigned int x1l, x1u, x2l, x2u;
  Real t, u;
  findInterval(_x1, x1, x1l, x1u, t);
  findInterval(_x2, x2, x2l, x2u, u);

  // Take derivative along x1 axis
  // Note: sum from i = 1 as the first term is zero
  if (deriv_var == 1)
  {
    Real sample_deriv = 0.0;
    for (unsigned int i = 1; i < 4; ++i)
      for (unsigned int j = 0; j < 4; ++j)
        sample_deriv +=
            i * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i - 1) * MathUtils::pow(u, j);

    Real d = _x1[x1u] - _x1[x1l];

    if (!MooseUtils::absoluteFuzzyEqual(d, 0.0))
      sample_deriv /= d;

    return sample_deriv;
  }

  // Take derivative along x2 axis
  // Note: sum from j = 1 as the first term is zero
  else if (deriv_var == 2)
  {
    Real sample_deriv = 0.0;
    for (unsigned int i = 0; i < 4; ++i)
      for (unsigned int j = 1; j < 4; ++j)
        sample_deriv +=
            j * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) * MathUtils::pow(u, j - 1);

    Real d = _x2[x2u] - _x2[x2l];

    if (!MooseUtils::absoluteFuzzyEqual(d, 0.0))
      sample_deriv /= d;

    return sample_deriv;
  }

  else
    mooseError("deriv_var must be either 1 or 2 in BicubicInterpolation");
}

sampleValueAndDerivatives()

void BicubicInterpolation::sampleValueAndDerivatives	(	Real	x1,
		Real	x2,
		Real &	y,
		Real &	dy1,
		Real &	dy2
	)

Samples value and first derivatives at point (x1, x2) Use this function for speed when computing both value and derivatives, as it minimizes the amount of time spent locating the point in the tabulated data.

Definition at line 155 of file BicubicInterpolation.C.

{
  unsigned int x1l, x1u, x2l, x2u;
  Real t, u;
  findInterval(_x1, x1, x1l, x1u, t);
  findInterval(_x2, x2, x2l, x2u, u);

  y = 0.0;
  for (unsigned int i = 0; i < 4; ++i)
    for (unsigned int j = 0; j < 4; ++j)
      y += _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) * MathUtils::pow(u, j);

  // Note: sum from i = 1 as the first term is zero
  dy1 = 0.0;
  for (unsigned int i = 1; i < 4; ++i)
    for (unsigned int j = 0; j < 4; ++j)
      dy1 += i * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i - 1) * MathUtils::pow(u, j);

  // Note: sum from j = 1 as the first term is zero
  dy2 = 0.0;
  for (unsigned int i = 0; i < 4; ++i)
    for (unsigned int j = 1; j < 4; ++j)
      dy2 += j * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) * MathUtils::pow(u, j - 1);

  Real d1 = _x1[x1u] - _x1[x1l];

  if (!MooseUtils::absoluteFuzzyEqual(d1, 0.0))
    dy1 /= d1;

  Real d2 = _x2[x2u] - _x2[x2l];

  if (!MooseUtils::absoluteFuzzyEqual(d2, 0.0))
    dy2 /= d2;
}

tableDerivatives()

void BicubicInterpolation::tableDerivatives ( std::vector< std::vector< Real >> &  dy_dx1, std::vector< std::vector< Real >> &  dy_dx2, std::vector< std::vector< Real >> &  d2y_dx1x2  )

protected

Provides the values of the first derivatives in each direction at all points in the table of data, as well as the mixed second derivative.

This is implemented using finite differencing, but could be supplied through other means (such as by sampling with cubic splines)

Definition at line 249 of file BicubicInterpolation.C.

Referenced by precomputeCoefficients().

{
  auto m = _x1.size(), n = _x2.size();
  dy_dx1.resize(m);
  dy_dx2.resize(m);
  d2y_dx1x2.resize(m);

  for (auto i = beginIndex(_x1); i < m; ++i)
  {
    dy_dx1[i].resize(n);
    dy_dx2[i].resize(n);
    d2y_dx1x2[i].resize(n);
  }

  // Central difference calculations of derivatives
  for (auto i = beginIndex(_y); i < m; ++i)
    for (auto j = beginIndex(_y[i]); j < n; ++j)
    {
      // Derivative wrt x1
      if (i == 0)
      {
        dy_dx1[i][j] = (_y[i + 1][j] - _y[i][j]) / (_x1[i + 1] - _x1[i]);
      }
      else if (i == m - 1)
        dy_dx1[i][j] = (_y[i][j] - _y[i - 1][j]) / (_x1[i] - _x1[i - 1]);
      else
        dy_dx1[i][j] = (_y[i + 1][j] - _y[i - 1][j]) / (_x1[i + 1] - _x1[i - 1]);

      // Derivative wrt x2
      if (j == 0)
        dy_dx2[i][j] = (_y[i][j + 1] - _y[i][j]) / (_x2[j + 1] - _x2[j]);
      else if (j == n - 1)
        dy_dx2[i][j] = (_y[i][j] - _y[i][j - 1]) / (_x2[j] - _x2[j - 1]);
      else
        dy_dx2[i][j] = (_y[i][j + 1] - _y[i][j - 1]) / (_x2[j + 1] - _x2[j - 1]);

      // Mixed derivative d2y/dx1x2
      if (i == 0 && j == 0)
        d2y_dx1x2[i][j] = (_y[i + 1][j + 1] - _y[i + 1][j] - _y[i][j + 1] + _y[i][j]) /
                          (_x1[i + 1] - _x1[i]) / (_x2[j + 1] - _x2[j]);
      else if (i == 0 && j == n - 1)
        d2y_dx1x2[i][j] = (_y[i + 1][j] - _y[i + 1][j - 1] - _y[i][j] + _y[i][j - 1]) /
                          (_x1[i + 1] - _x1[i]) / (_x2[j] - _x2[j - 1]);
      else if (i == m - 1 && j == 0)
        d2y_dx1x2[i][j] = (_y[i][j + 1] - _y[i][j] - _y[i - 1][j + 1] + _y[i - 1][j]) /
                          (_x1[i] - _x1[i - 1]) / (_x2[j + 1] - _x2[j]);
      else if (i == m - 1 && j == n - 1)
        d2y_dx1x2[i][j] = (_y[i][j] - _y[i][j - 1] - _y[i - 1][j] + _y[i - 1][j - 1]) /
                          (_x1[i] - _x1[i - 1]) / (_x2[j] - _x2[j - 1]);
      else if (i == 0)
        d2y_dx1x2[i][j] = (_y[i + 1][j + 1] - _y[i + 1][j - 1] - _y[i][j + 1] + _y[i][j - 1]) /
                          (_x1[i + 1] - _x1[i]) / (_x2[j + 1] - _x2[j - 1]);
      else if (i == m - 1)
        d2y_dx1x2[i][j] = (_y[i][j + 1] - _y[i][j - 1] - _y[i - 1][j + 1] + _y[i - 1][j - 1]) /
                          (_x1[i] - _x1[i - 1]) / (_x2[j + 1] - _x2[j - 1]);
      else if (j == 0)
        d2y_dx1x2[i][j] = (_y[i + 1][j + 1] - _y[i + 1][j] - _y[i - 1][j + 1] + _y[i - 1][j]) /
                          (_x1[i + 1] - _x1[i - 1]) / (_x2[j + 1] - _x2[j]);
      else if (j == n - 1)
        d2y_dx1x2[i][j] = (_y[i + 1][j] - _y[i + 1][j - 1] - _y[i - 1][j] + _y[i - 1][j - 1]) /
                          (_x1[i + 1] - _x1[i - 1]) / (_x2[j] - _x2[j - 1]);
      else
        d2y_dx1x2[i][j] =
            (_y[i + 1][j + 1] - _y[i + 1][j - 1] - _y[i - 1][j + 1] + _y[i - 1][j - 1]) /
            (_x1[i + 1] - _x1[i - 1]) / (_x2[j + 1] - _x2[j - 1]);
    }
}

Member Data Documentation

_bicubic_coeffs

std::vector<std::vector<std::vector<std::vector<Real> > > > BicubicInterpolation::_bicubic_coeffs

protected

Matrix of precomputed coefficients There are four coefficients in each direction at each dependent value.

Definition at line 98 of file BicubicInterpolation.h.

Referenced by BicubicInterpolation(), precomputeCoefficients(), sample(), sample2ndDerivative(), sampleDerivative(), and sampleValueAndDerivatives().

_wt

const std::vector<std::vector<int> > BicubicInterpolation::_wt

protected

Initial value:

{
      {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
      {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
      {-3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1, 0, 0, 0, 0},
      {2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0},
      {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
      {0, 0, 0, 0, -3, 0, 0, 3, 0, 0, 0, 0, -2, 0, 0, -1},
      {0, 0, 0, 0, 2, 0, 0, -2, 0, 0, 0, 0, 1, 0, 0, 1},
      {-3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
      {0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0},
      {9, -9, 9, -9, 6, 3, -3, -6, 6, -6, -3, 3, 4, 2, 1, 2},
      {-6, 6, -6, 6, -4, -2, 2, 4, -3, 3, 3, -3, -2, -1, -1, -2},
      {2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
      {0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0},
      {-6, 6, -6, 6, -3, -3, 3, 3, -4, 4, 2, -2, -2, -2, -1, -1},
      {4, -4, 4, -4, 2, 2, -2, -2, 2, -2, -2, 2, 1, 1, 1, 1}}

Matrix used to calculate bicubic interpolation coefficients (from Numerical Recipes)

Definition at line 102 of file BicubicInterpolation.h.

Referenced by precomputeCoefficients().

_x1

std::vector<Real> BicubicInterpolation::_x1

protected

Independent values in the x1 direction.

Definition at line 90 of file BicubicInterpolation.h.

Referenced by BicubicInterpolation(), errorCheck(), precomputeCoefficients(), sample(), sample2ndDerivative(), sampleDerivative(), sampleValueAndDerivatives(), and tableDerivatives().

_x2

std::vector<Real> BicubicInterpolation::_x2

protected

Independent values in the x2 direction.

Definition at line 92 of file BicubicInterpolation.h.

Referenced by BicubicInterpolation(), errorCheck(), precomputeCoefficients(), sample(), sample2ndDerivative(), sampleDerivative(), sampleValueAndDerivatives(), and tableDerivatives().

_y

std::vector<std::vector<Real> > BicubicInterpolation::_y

protected

The dependent values at (x1, x2) points.

Definition at line 94 of file BicubicInterpolation.h.

Referenced by errorCheck(), precomputeCoefficients(), and tableDerivatives().

---

The documentation for this class was generated from the following files:

• include/utils/BicubicInterpolation.h
• src/utils/BicubicInterpolation.C