Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/math/cubic_spline.hh
8ae23a91b711features/coupler_compon…

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 <cmath>
#include <vector>
/* -------------------------------------------------------------------------- */
__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