Homec4science
Diffusion LibMultiScale (features/reshape_compute_python) Browse

libmultiscale/src/math/cubic_spline.hh
b05d5f2aba78features/reshape_comput…

Last Change
Raw File

cubic_spline.hh
View Options

/**
* @file cubic_spline.hh
*
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Dec 11 12:36:48 2013
*
* @brief This is a cubic spline interpolator
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __LIBMULTISCALE_CUBIC_SPLINE_HH__
#define __LIBMULTISCALE_CUBIC_SPLINE_HH__
/* -------------------------------------------------------------------------- */
#include <vector>
#include <cmath>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
class CubicSpline
{
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
CubicSpline(){
bc_first_point = 0.0;
bc_last_point = 0.0;
spacing = 1.0;
spacing_inverse = 1.0;
sixinv=1.0/6.0;
};
~CubicSpline(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
//! build the cubic spline coefficients
void buildSplineCoefficients(std::vector<Real> &x, std::vector<Real> &f, Real spacing);
//! evaluate the function at x using cubic spline interoplation
void evaluateSpline(Real x, Real &y);
//! fast evaluate the function at x using cubic spline interoplation works only for equally spaced data points
void fastEvaluateSpline(Real x, Real &y);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! set boundary conditions at first and last points
void setBCs(Real x0, Real xn);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
//private:
protected:
//! vector to store the x values
std::vector<Real> x_f;
//! vector to store the f(x) values
std::vector<Real> poly_f;
//! vector to store the coefficient b
std::vector<Real> b;
//! vector to store the coefficient c
std::vector<Real> c;
//! vector to store the coefficient d
std::vector<Real> d;
//! vector to store the second derivative of f(x) i.e. df_dyy
std::vector<Real> df_dyy;
//! boundar condition at first data point
Real bc_first_point;
//! boundar condition at last data point
Real bc_last_point;
//! equal spacing between points
Real spacing;
//! equal inverse spacing between points
Real spacing_inverse;
//! variable to store constant factor
Real sixinv;
//! number of data points
UInt n;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
//! build the cubic spline coefficients for not a knot boundary conditions
inline void CubicSpline::buildSplineCoefficients(std::vector<Real> &x, std::vector<Real> &f, Real spacing)
{
Real p,qn,sig,un;
std::vector<Real> u;
n = x.size();
// internal variables for later use
this->spacing = spacing;
spacing_inverse = 1./spacing;
x_f=x;
poly_f=f;
u.assign(n-1,0.0);
df_dyy.assign(n,0.0);
c.assign(n,0.0);
b.assign(n,0.0);
d.assign(n,0.0);
// set the boundary condition at first point
df_dyy[0] = -0.5;
u[0] = (3.0/(x[1]-x[0]))*((f[1]-f[0])/(x[1]-x[0])-bc_first_point);
for(UInt i = 1; i < n-1; i++){
sig = (x[i] - x[i-1])/(x[i+1] - x[i-1]);
p = sig*df_dyy[i-1] + 2.0;
df_dyy[i] = (sig - 1.0)/p;
u[i] = (f[i+1] - f[i])/(x[i+1] - x[i]) - (f[i] - f[i-1])/(x[i] - x[i-1]);
u[i] = (6.0*u[i]/(x[i+1] - x[i-1]) - sig*u[i-1])/p;
}
// set the boundary condition at last point
qn = 0.5;
un = (3.0/(x[n-1] - x[n-2]))*(bc_last_point - (f[n-1] - f[n-2])/(x[n-1] - x[n-2]));
// store second derivatives of f in df_dyy
df_dyy[n-1] = (un - qn*u[n-2])/(qn*df_dyy[n-2] + 1.0);
for(int k = n-2; k >= 0; k--){
df_dyy[k] = df_dyy[k]*df_dyy[k+1] + u[k];
}
// compute cubic spline polynomial coeffs b,c,d
for(UInt i = 0; i < n; ++i)
c[i] = 0.5*df_dyy[i];
for(UInt i = 0; i < n-1; ++i){
Real h = (x[i+1]-x[i]);
b[i] = (f[i+1]-f[i])/h - h*(c[i+1]+2*c[i])/3.0;
d[i] = (c[i+1]-c[i])/(h*3.0);
}
}
//! evaluate the function at x using cubic spline interoplation (useful for not equally spaced points)
inline void CubicSpline::evaluateSpline(Real x, Real &y)
{
static UInt pklo=0,pkhi=1;
Real h,b,a;
UInt klo,khi,k;
if(x_f[pklo] <= x && x_f[pkhi] > x){
klo = pklo;
khi = pkhi;
}
else{
klo = 0;
khi = n - 1;
while(khi - klo > 1){
k = (khi + klo) >> 1;
if(x_f[k] > x) khi = k;
else klo = k;
}
}
h = x_f[khi] - x_f[klo];
if(h == 0){
LM_FATAL("This should not happen. Cubic spline cannot be evaluated at this point " << x);
}
// compute the interpolated value
a = (x_f[khi] - x)/h;
b = (x - x_f[klo])/h;
y = a*poly_f[klo] + b*poly_f[khi] + ((a*a*a - a)*df_dyy[klo] + (b*b*b - b)*df_dyy[khi])*(h*h)*sixinv;
}
inline void CubicSpline::fastEvaluateSpline(Real x, Real &y)
{
UInt i=0;
// LM_ASSERT(x < x_f[0],"Cubic spline interpolation cannot be evaluated for an external point\n");
i=(x-x_f[0])*spacing_inverse;
if (i >= x_f.size()-1) {
y = 0.0;
return ;
}
Real h=x-x_f[i];
// compute the interpolated value;
y=poly_f[i]+h*(b[i]+h*(c[i]+h*d[i]));
LM_ASSERT(!std::isnan(y),"problem");
LM_ASSERT(!std::isinf(y),"problem");
}
inline void CubicSpline::setBCs(Real x0, Real xn)
{
bc_first_point = x0;
bc_last_point = xn;
}
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_CUBIC_SPLINE_HH__ */
c4science · Help