tension_spline_blueprints.cpp
Go to the documentation of this file.
1 
10 /*****************************************************************************
11 ** Includes
12 *****************************************************************************/
13 
14 #include <cmath>
15 #include <ecl/containers/array.hpp>
17 #include "../../include/ecl/geometry/tension_spline.hpp"
18 
19 /*****************************************************************************
20 ** Namespaces
21 *****************************************************************************/
22 
23 namespace ecl {
24 namespace blueprints {
25 
26 /*****************************************************************************
27 ** Using
28 *****************************************************************************/
29 
30 using ecl::Array;
33 
34 /*****************************************************************************
35 ** Implementation [C2CubicSpline]
36 *****************************************************************************/
37 /*
38  * This algorithm can be found on page 323 of "Numerical Analysis" [Kincaid &
39  * Cheney].
40  */
42  const double &tau) ecl_assert_throw_decl(ecl::StandardException) :
43  x_data(x_set),
44  y_data(y_set)
45 {
46  ecl_assert_throw( x_data.size() == y_data.size(), StandardException(LOC,OutOfRangeError) );
47 
48  /*********************
49  ** Basic Parameters
50  **********************/
51  unsigned int n = x_data.size()-1; // No. of data points (x_0...x_n)
52  yddot_data.resize(n+1);
53  yddot_data[0] = 0;
54  tension = tau;
55 
56  /******************************************
57  ** Tridiagonal parameters
58  *******************************************/
59  // These equations are derived from the c2 constraints - it results in a tridiagonal system
60  // of equations. See p.319 of kincaid and cheney for the tridiagonal solution (note this
61  // is only valid for diagonally dominant tridiagonal matrices! Is this a problem for us?
62  Array<double> h(n), a(n), beta(n), gamma(n), u(n), v(n);
63  h[0] = x_set[1] - x_set[0];
64  for ( unsigned int i = 0; i < n; ++i ) {
65  h[i] = x_set[i+1]-x_set[i];
66  a[i] = 1/h[i] - tension/sinh(tension*h[i]);
67  beta[i] = tension*(cosh(tension*h[i])/sinh(tension*h[i])) - 1/h[i];
68  gamma[i] = tension*tension*(y_data[i+1]-y_data[i])/h[i];
69  }
70  /*
71  * [ u1 a1 0 0 0 ] [z1] [ v1 ]
72  * [ a1 u2 a2 0 0 ] [z2] [ v2 ]
73  * [ 0 a2 u3 a3 0 ] [z3] [ v3 ]
74  *
75  * but with back substitution here,
76  *
77  * [ u1 a1 0 0 0 ] [z1] [ v1 ]
78  * [ 0 (u2-a1^2/u1) a2 0 0 ] [z2] [ v2-a1*v1/u1 ]
79  *
80  * etc
81  */
82  u[1] = beta[1] + beta[0];
83  v[1] = gamma[1]-gamma[0];
84  for ( unsigned int i = 2; i < n; ++i ) {
85  u[i] = beta[i] + beta[i-1] - a[i-1]*a[i-1]/u[i-1];
86  v[i] = gamma[i]-gamma[i-1] - a[i-1]*v[i-1]/u[i-1];
87  }
88  /*
89  * Now we can roll the upper triangular matrix backwards
90  */
91  yddot_data[n] = 0;
92  for ( unsigned int i = n-1; i > 0; --i ) {
93  yddot_data[i] = (v[i] - a[i]*yddot_data[i+1])/u[i];
94  }
95  yddot_data[0] = 0;
96 }
97 
98 TensionSpline C2TensionSpline::instantiate() {
99  TensionSpline cubic;
100  apply(cubic);
101  return cubic;
102 };
103 
105 
106  spline.discretised_domain = x_data;
107  spline.tension = tension;
108  spline.functions.resize(x_data.size()-1); // One less polynomials than there is points
109  for (unsigned int i = 0; i < spline.functions.size(); ++i ) {
111  x_data[i], y_data[i], yddot_data[i],
112  x_data[i+1], y_data[i+1], yddot_data[i+1] );
113  }
114 };
115 
116 } // namespace blueprints
117 
118 /*****************************************************************************
119 ** Implementation [BluePrintFactory][TensionSpline]
120 *****************************************************************************/
121 
123 
124 
126  const Array<double>& y_set, const double &tau) ecl_assert_throw_decl(StandardException) {
127  return C2TensionSpline(x_set, y_set, tau);
128 };
129 
130 } // namespace ecl
Embedded control libraries.
Storage container for a tension spline interpolation.
Array< double > discretised_domain
Blueprint for generating a tension spline satisfying C2 constraints.
C2TensionSpline()
Default constructor.
void apply(base_type &spline) const
Apply the blueprint to configure an existing object.
static blueprints::TensionSecondDerivativeInterpolation Interpolation(const double x_i, const double y_i, const double yddot_i, const double x_f, const double y_f, const double yddot_f)
Blueprint for interpolating a tension function between two end points using second derivatives...
C2TensionSpline(const ecl::Array< double > &x_set, const ecl::Array< double > &y_set, const double &tau) ecl_assert_throw_decl(ecl base_type instantiate()
Constructor that properly configures/initialises the blueprint.
#define ecl_assert_throw(expression, exception)
OutOfRangeError
static size_type size()
#define ecl_assert_throw_decl(exception)
Array< TensionFunction > functions


ecl_geometry
Author(s): Daniel Stonier
autogenerated on Mon Jun 10 2019 13:08:37