tension_spline_blueprints.cpp
Go to the documentation of this file.
00001 
00010 /*****************************************************************************
00011 ** Includes
00012 *****************************************************************************/
00013 
00014 #include <cmath>
00015 #include <ecl/containers/array.hpp>
00016 #include <ecl/exceptions/standard_exception.hpp>
00017 #include "../../include/ecl/geometry/tension_spline.hpp"
00018 
00019 /*****************************************************************************
00020 ** Namespaces
00021 *****************************************************************************/
00022 
00023 namespace ecl {
00024 namespace blueprints {
00025 
00026 /*****************************************************************************
00027 ** Using
00028 *****************************************************************************/
00029 
00030 using ecl::Array;
00031 using ecl::OutOfRangeError;
00032 using ecl::StandardException;
00033 
00034 /*****************************************************************************
00035 ** Implementation [C2CubicSpline]
00036 *****************************************************************************/
00037 /*
00038  * This algorithm can be found on page 323 of "Numerical Analysis" [Kincaid &
00039  * Cheney].
00040  */
00041 C2TensionSpline::C2TensionSpline(const ecl::Array<double>& x_set, const ecl::Array<double>& y_set,
00042                 const double &tau) ecl_assert_throw_decl(ecl::StandardException) :
00043     x_data(x_set),
00044     y_data(y_set)
00045 {
00046     ecl_assert_throw( x_data.size() == y_data.size(), StandardException(LOC,OutOfRangeError) );
00047 
00048     /*********************
00049     ** Basic Parameters
00050     **********************/
00051     unsigned int n = x_data.size()-1; // No. of data points (x_0...x_n)
00052     yddot_data.resize(n+1);
00053     yddot_data[0] = 0;
00054     tension = tau;
00055 
00056     /******************************************
00057     ** Tridiagonal parameters
00058     *******************************************/
00059     // These equations are derived from the c2 constraints - it results in a tridiagonal system
00060     // of equations. See p.319 of kincaid and cheney for the tridiagonal solution (note this
00061     // is only valid for diagonally dominant tridiagonal matrices! Is this a problem for us?
00062     Array<double> h(n), a(n), beta(n), gamma(n), u(n), v(n);
00063     h[0] = x_set[1] - x_set[0];
00064     for ( unsigned int i = 0; i < n; ++i ) {
00065         h[i] = x_set[i+1]-x_set[i];
00066         a[i] = 1/h[i] - tension/sinh(tension*h[i]);
00067         beta[i] = tension*(cosh(tension*h[i])/sinh(tension*h[i])) - 1/h[i];
00068         gamma[i] = tension*tension*(y_data[i+1]-y_data[i])/h[i];
00069     }
00070     /*
00071      * [ u1 a1 0   0   0 ] [z1]    [ v1 ]
00072      * [ a1 u2 a2  0   0 ] [z2]    [ v2 ]
00073      * [  0 a2 u3 a3   0 ] [z3]    [ v3 ]
00074      *
00075      * but with back substitution here,
00076      *
00077      * [ u1      a1       0   0   0 ] [z1]    [     v1      ]
00078      * [  0 (u2-a1^2/u1) a2   0   0 ] [z2]    [ v2-a1*v1/u1 ]
00079      *
00080      * etc
00081      */
00082     u[1] = beta[1] + beta[0];
00083     v[1] = gamma[1]-gamma[0];
00084     for ( unsigned int i = 2; i < n; ++i ) {
00085         u[i] = beta[i] + beta[i-1] - a[i-1]*a[i-1]/u[i-1];
00086         v[i] = gamma[i]-gamma[i-1] - a[i-1]*v[i-1]/u[i-1];
00087     }
00088     /*
00089      * Now we can roll the upper triangular matrix backwards
00090      */
00091     yddot_data[n] = 0;
00092     for ( unsigned int i = n-1; i > 0; --i ) {
00093         yddot_data[i] = (v[i] - a[i]*yddot_data[i+1])/u[i];
00094     }
00095     yddot_data[0] = 0;
00096 }
00097 
00098 TensionSpline C2TensionSpline::instantiate() {
00099     TensionSpline cubic;
00100     apply(cubic);
00101     return cubic;
00102 };
00103 
00104 void C2TensionSpline::apply(TensionSpline& spline) const {
00105 
00106     spline.discretised_domain = x_data;
00107     spline.tension = tension;
00108     spline.functions.resize(x_data.size()-1); // One less polynomials than there is points
00109     for (unsigned int i = 0; i < spline.functions.size(); ++i ) {
00110         spline.functions[i] = TensionFunction::Interpolation(
00111                         x_data[i],   y_data[i],   yddot_data[i],
00112                         x_data[i+1], y_data[i+1], yddot_data[i+1]  );
00113     }
00114 };
00115 
00116 } // namespace blueprints
00117 
00118 /*****************************************************************************
00119 ** Implementation [BluePrintFactory][TensionSpline]
00120 *****************************************************************************/
00121 
00122 using blueprints::C2TensionSpline;
00123 
00124 
00125 C2TensionSpline BluePrintFactory< TensionSpline >::Natural(const Array<double>& x_set,
00126                 const Array<double>& y_set, const double &tau) ecl_assert_throw_decl(StandardException) {
00127     return C2TensionSpline(x_set, y_set, tau);
00128 };
00129 
00130 } // namespace ecl


ecl_geometry
Author(s): Daniel Stonier
autogenerated on Sun Oct 5 2014 23:36:08