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 "TrilinearInterpolation.h"
      11             : #include "MooseError.h"
      12             : 
      13           6 : TrilinearInterpolation::TrilinearInterpolation(const std::vector<Real> & x,
      14             :                                                const std::vector<Real> & y,
      15             :                                                const std::vector<Real> & z,
      16           6 :                                                const std::vector<Real> & data)
      17           6 :   : _x_axis(x), _y_axis(y), _z_axis(z), _fxyz(data)
      18             : {
      19           6 :   if (_x_axis.size() < 1)
      20           0 :     mooseError("x vector has zero elements. At least one element is required.");
      21           6 :   if (_y_axis.size() < 1)
      22           0 :     mooseError("y vector has zero elements. At least one element is required.");
      23           6 :   if (_z_axis.size() < 1)
      24           0 :     mooseError("z vector has zero elements. At least one element is required.");
      25           6 :   if (_x_axis.size() * _y_axis.size() * _z_axis.size() != data.size())
      26           1 :     mooseError("The size of data (",
      27           1 :                data.size(),
      28             :                ") does not match the supplied dimensions (",
      29           2 :                _x_axis.size(),
      30             :                ", ",
      31           2 :                _y_axis.size(),
      32             :                ", ",
      33           2 :                _z_axis.size(),
      34             :                ")");
      35           9 : }
      36             : 
      37             : void
      38         105 : TrilinearInterpolation::getCornerIndices(
      39             :     const std::vector<Real> & v, Real x, int & lower, int & upper, Real & d) const
      40             : {
      41         105 :   unsigned int N = v.size();
      42         105 :   if (x < v[0])
      43             :   {
      44           4 :     lower = 0;
      45           4 :     upper = 0;
      46             :   }
      47         101 :   else if (x >= v[N - 1])
      48             :   {
      49          47 :     lower = N - 1;
      50          47 :     upper = N - 1;
      51             :   }
      52             :   else
      53             :   {
      54          54 :     for (unsigned int i = 0; i < N - 1; i++)
      55             :     {
      56          54 :       if (x > v[i] && x < v[i + 1])
      57             :       {
      58          27 :         lower = i;
      59          27 :         upper = i + 1;
      60          27 :         d = (x - v[lower]) / (v[upper] - v[lower]);
      61          27 :         break;
      62             :       }
      63          27 :       else if (x == v[i])
      64             :       {
      65          27 :         lower = i;
      66          27 :         upper = i;
      67          27 :         break;
      68             :       }
      69             :     }
      70             :   }
      71         105 : }
      72             : 
      73             : Real
      74         280 : TrilinearInterpolation::getCornerValues(int x, int y, int z) const
      75             : {
      76         280 :   int nY = _y_axis.size();
      77         280 :   int nZ = _z_axis.size();
      78             : 
      79         280 :   return _fxyz[x * nY * nZ + y * nZ + z];
      80             : }
      81             : 
      82             : Real
      83          35 : TrilinearInterpolation::sample(Real x, Real y, Real z) const
      84             : {
      85          35 :   int x0 = 0;
      86          35 :   int y0 = 0;
      87          35 :   int z0 = 0;
      88          35 :   int x1 = 0;
      89          35 :   int y1 = 0;
      90          35 :   int z1 = 0;
      91          35 :   Real Dx = 0;
      92          35 :   Real Dy = 0;
      93          35 :   Real Dz = 0;
      94             : 
      95             :   // find the the indices of the cube, which contains the point
      96          35 :   getCornerIndices(_x_axis, x, x0, x1, Dx);
      97          35 :   getCornerIndices(_y_axis, y, y0, y1, Dy);
      98          35 :   getCornerIndices(_z_axis, z, z0, z1, Dz);
      99             : 
     100             :   // find the corresponding function values for the corner indices
     101          35 :   Real f000 = getCornerValues(x0, y0, z0);
     102          35 :   Real f001 = getCornerValues(x0, y0, z1);
     103          35 :   Real f010 = getCornerValues(x0, y1, z0);
     104          35 :   Real f011 = getCornerValues(x0, y1, z1);
     105          35 :   Real f100 = getCornerValues(x1, y0, z0);
     106          35 :   Real f101 = getCornerValues(x1, y0, z1);
     107          35 :   Real f110 = getCornerValues(x1, y1, z0);
     108          35 :   Real f111 = getCornerValues(x1, y1, z1);
     109             : 
     110             :   // interpolation
     111          35 :   Real f00 = (f100 - f000) * Dx + f000;
     112          35 :   Real f10 = (f110 - f010) * Dx + f010;
     113          35 :   Real f01 = (f101 - f001) * Dx + f001;
     114          35 :   Real f11 = (f111 - f011) * Dx + f011;
     115          35 :   Real f0 = (f10 - f00) * Dy + f00;
     116          35 :   Real f1 = (f11 - f01) * Dy + f01;
     117             : 
     118          35 :   return (f1 - f0) * Dz + f0;
     119             : }