SplineInterpolatorOffline.cc

Go to the documentation of this file.

#include "SplineInterpolator.h"

#include <vector>
#include <boost/multi_array.hpp>

using namespace utl::Spline;

Interpolator1D::Interpolator1D(const double xStart,
                               const double xStop, /* inclusive! */
                               const std::vector<double>& values,
                               const BoundaryCondition& bcX)
{
  const size_t xn = values.size();

  fSpline.Configure(0, KnotVector(xn, xStart, xStop));

  fData.resize(xn+2);
  for (size_t i = 0; i < xn; ++i)
    fData[i] = values[i];

  Solver<1, KnotVector, BasisFunction> solver;
  solver.Configure(0, fSpline.GetKnotVector(), fSpline.GetBasisVector(), bcX);
  solver(&fData.front());
}

double
Interpolator1D::operator()(const double x)
  const
{
  return fSpline(&fData.front(), x);
}

Interpolator2D::Interpolator2D(const double xStart,
                               const double xStop, /* inclusive! */
                               const double yStart,
                               const double yStop, /* inclusive! */
                               const boost::multi_array<double,2>& values,
                               const BoundaryCondition bcX,
                               const BoundaryCondition bcY)
{
  const size_t xn = values.shape()[0];
  const size_t yn = values.shape()[1];

  fSpline.Configure(0, KnotVector(xn, xStart, xStop));
  fSpline.Configure(1, KnotVector(yn, yStart, yStop));

  fData.resize((xn+2)*(yn+2));
  for (size_t i = 0; i < xn; ++i)
    for (size_t j = 0; j < yn; ++j)
      fData[j*(xn+2) + i] = values[i][j];

  Solver<2, KnotVector, BasisFunction> solver;
  solver.Configure(0, fSpline.GetKnotVector(0), fSpline.GetBasisVector(0), bcX);
  solver.Configure(1, fSpline.GetKnotVector(1), fSpline.GetBasisVector(1), bcY);
  solver(&fData.front());
}

double
Interpolator2D::operator()(const double x, const double y)
  const
{
  fX[0] = x;
  fX[1] = y;
  return fSpline(&fData.front(), fX);
}

Interpolator3D::Interpolator3D(const double xStart,
                               const double xStop, /* inclusive! */
                               const double yStart,
                               const double yStop, /* inclusive! */
                               const double zStart,
                               const double zStop, /* inclusive! */
                               const boost::multi_array<double,3>& values,
                               const BoundaryCondition bcX,
                               const BoundaryCondition bcY,
                               const BoundaryCondition bcZ)
{
  const size_t xn = values.shape()[0];
  const size_t yn = values.shape()[1];
  const size_t zn = values.shape()[2];

  fSpline.Configure(0, KnotVector(xn, xStart, xStop));
  fSpline.Configure(1, KnotVector(yn, yStart, yStop));
  fSpline.Configure(2, KnotVector(zn, zStart, zStop));

  fData.resize((xn+2)*(yn+2)*(zn+2));
  for (size_t i = 0; i < xn; ++i)
    for (size_t j = 0; j < yn; ++j)
      for (size_t k = 0; k < zn; ++k)
        fData[k*(yn+2)*(xn+2) + j*(xn+2) + i] = values[i][j][k];

  Solver<3, KnotVector, BasisFunction> solver;
  solver.Configure(0, fSpline.GetKnotVector(0), fSpline.GetBasisVector(0), bcX);
  solver.Configure(1, fSpline.GetKnotVector(1), fSpline.GetBasisVector(1), bcY);
  solver.Configure(2, fSpline.GetKnotVector(2), fSpline.GetBasisVector(2), bcZ);
  solver(&fData.front());
}

double
Interpolator3D::operator()(const double x, const double y, const double z)
  const
{
  fX[0] = x;
  fX[1] = y;
  fX[2] = z;
  return fSpline(&fData.front(),fX);
}

VectorInterpolator1D::VectorInterpolator1D(const double xStart,
                                           const double xStop, /* inclusive! */
                                           const boost::multi_array<double,2>& fValue, // first index is function vector index!
                                           const BoundaryCondition& bcX)
{
  const size_t fn = fValue.shape()[0];
  const size_t xn = fValue.shape()[1];

  fSpline.Configure(0, KnotVector(xn, xStart, xStop));

  fData.resize(fn, std::vector<double>(xn+2));

  Solver<1, KnotVector, BasisFunction> solver;
  solver.Configure(0, fSpline.GetKnotVector(), fSpline.GetBasisVector(), bcX);
  for (size_t m = 0; m < fn; ++m)
  {
    for (size_t i = 0; i < xn; ++i)
      fData[m][i] = fValue[m][i];
    solver(&fData[m].front());
  }
}

void
VectorInterpolator1D::operator()(std::vector<double>& result, const double x)
  const
{
  for (size_t m = 0; m < fData.size(); ++m)
    result[m] = fSpline(&fData[m].front(),x);
}

VectorInterpolator2D::VectorInterpolator2D(const double xStart,
                                           const double xStop, /* inclusive! */
                                           const double yStart,
                                           const double yStop, /* inclusive! */
                                           const boost::multi_array<double,3>& values, // first index is function vector index!
                                           const BoundaryCondition& bcX,
                                           const BoundaryCondition& bcY)
{
  const size_t fn = values.shape()[0];
  const size_t xn = values.shape()[1];
  const size_t yn = values.shape()[2];

  fSpline.Configure(0, KnotVector(xn, xStart, xStop));
  fSpline.Configure(1, KnotVector(yn, yStart, yStop));

  fData.resize(fn, std::vector<double>((xn+2)*(yn+2)));

  Solver<2, KnotVector, BasisFunction> solver;
  solver.Configure(0, fSpline.GetKnotVector(0), fSpline.GetBasisVector(0), bcX);
  solver.Configure(1, fSpline.GetKnotVector(1), fSpline.GetBasisVector(1), bcY);
  for (size_t m = 0; m < fn; ++m)
  {
    for (size_t i = 0; i < xn; ++i)
      for (size_t j = 0; j < yn; ++j)
        fData[m][j*(xn+2) + i] = values[m][i][j];
    solver(&fData[m].front());
  }
}

void
VectorInterpolator2D::operator()(std::vector<double>& result, const double x, const double y)
  const
{
  fX[0] = x;
  fX[1] = y;
  for (size_t m = 0; m < fData.size(); ++m)
    result[m] = fSpline(&fData[m].front(), fX);
}

VectorInterpolator3D::VectorInterpolator3D(const double xStart,
                                           const double xStop, /* inclusive! */
                                           const double yStart,
                                           const double yStop, /* inclusive! */
                                           const double zStart,
                                           const double zStop, /* inclusive! */
                                           const boost::multi_array<double,4>& values, // first index is function vector index!
                                           const BoundaryCondition& bcX,
                                           const BoundaryCondition& bcY,
                                           const BoundaryCondition& bcZ)
{
  const size_t fn = values.shape()[0];
  const size_t xn = values.shape()[1];
  const size_t yn = values.shape()[2];
  const size_t zn = values.shape()[3];

  fSpline.Configure(0, KnotVector(xn, xStart, xStop));
  fSpline.Configure(1, KnotVector(yn, yStart, yStop));
  fSpline.Configure(2, KnotVector(zn, zStart, zStop));

  fData.resize(fn, std::vector<double>((xn+2)*(yn+2)*(zn+2)));

  Solver<3, KnotVector, BasisFunction> solver;
  solver.Configure(0, fSpline.GetKnotVector(0), fSpline.GetBasisVector(0), bcX);
  solver.Configure(1, fSpline.GetKnotVector(1), fSpline.GetBasisVector(1), bcY);
  solver.Configure(2, fSpline.GetKnotVector(2), fSpline.GetBasisVector(2), bcZ);
  for (size_t m = 0; m < fn; ++m)
  {
    for (size_t i = 0; i < xn; ++i)
      for (size_t j = 0; j < yn; ++j)
        for (size_t k = 0; k < zn; ++k)
          fData[m][k*(yn+2)*(xn+2) + j*(xn+2) + i] = values[m][i][j][k];
    solver(&fData[m].front());
  }
}

void
VectorInterpolator3D::operator()(std::vector<double>& result, const double x, const double y, const double z)
  const
{
  fX[0] = x;
  fX[1] = y;
  fX[2] = z;
  for (size_t m = 0; m < fData.size(); ++m)
    result[m] = fSpline(&fData[m].front(), fX);
}