cubic_spline_blueprints.cpp
Go to the documentation of this file.
1 
12 /*****************************************************************************
13 ** Includes
14 *****************************************************************************/
15 
16 #include <ecl/containers/array.hpp>
18 #include "../../include/ecl/geometry/cubic_spline.hpp"
19 
20 /*****************************************************************************
21 ** Namespaces
22 *****************************************************************************/
23 
24 namespace ecl {
25 namespace blueprints {
26 
27 /*****************************************************************************
28 ** Using
29 *****************************************************************************/
30 
31 using ecl::Array;
34 
35 /*****************************************************************************
36 ** Implementation [C2CubicSpline]
37 *****************************************************************************/
38 /*
39  * This algorithm can be found in page 115 of "Numerical Recipes in C"
40  * [The art of scientific computing] by Press, Vetterling, Tuekolsky, Flannery.
41  */
43  const double ydot_0, const double ydot_f) ecl_assert_throw_decl(StandardException) :
44  x_data(x_set),
45  y_data(y_set)
46 {
47  if (x_data.size() < 2 || y_data.size() < 2)
48  {
49  ecl_throw(StandardException(LOC,OutOfRangeError) );
50  }
51  ecl_assert_throw( x_data.size() == y_data.size(), StandardException(LOC,OutOfRangeError) );
52  unsigned int n = x_data.size();
53  yddot_data.resize(n);
54  Array<double> u(n); // u is a temporary used in the algorithm
55 
56  // Initial boundary conditions
57  yddot_data[0] = -0.5;
58  u[0] = (3.0/(x_data[1]-x_data[0])) * ((y_data[1]-y_data[0])/(x_data[1]-x_data[0])-ydot_0);
59  // Set up the tridiagonal matrix to solve for the accelerations
60  for (unsigned int i = 1; i <= n-2; ++i){
61  double sig = (x_data[i]-x_data[i-1]) / (x_data[i+1]-x_data[i-1]);
62  double p = sig*yddot_data[i-1]+2.0;
63  yddot_data[i] = (sig-1.0)/p;
64  u[i] = (y_data[i+1]-y_data[i])/(x_data[i+1]-x_data[i]) -
65  (y_data[i]-y_data[i-1])/(x_data[i]-x_data[i-1]);
66  u[i] = (6.0*u[i]/(x_data[i+1]-x_data[i-1]) - sig*u[i-1])/p;
67  }
68  // Final boundary conditions
69  double qn = 0.5;
70  u[n-1] = (3.0/(x_data[n-1]-x_data[n-2])) * (ydot_f - (y_data[n-1]-y_data[n-2])/(x_data[n-1]-x_data[n-2]));
71 
72  // Back substitution loop of the tridiagonal algorithm
73  yddot_data[n-1] = ( u[n-1] - qn*u[n-2]) / ( qn*yddot_data[n-2] + 1.0 );
74  for ( int k = n-2; k >= 0; --k ) {
75  yddot_data[k] = yddot_data[k]*yddot_data[k+1] + u[k];
76  }
77 }
78 
79 /*
80  * This is a special case of the above. You'll note that only initial and
81  * boundary conditions change here.
82  */
83 C2CubicSpline::C2CubicSpline(const ecl::Array<double>& x_set, const ecl::Array<double>& y_set) ecl_assert_throw_decl(StandardException) :
84  x_data(x_set),
85  y_data(y_set)
86 {
87  if (x_data.size() < 2 || y_data.size() < 2)
88  {
89  ecl_throw(StandardException(LOC,OutOfRangeError) );
90  }
91  ecl_assert_throw( x_data.size() == y_data.size(), StandardException(LOC,OutOfRangeError) );
92  unsigned int n = x_data.size();
93  yddot_data.resize(n);
94  Array<double> u(n); // u is a temporary used in the algorithm
95 
96  // Initial boundary conditions
97  yddot_data[0] = 0.0;
98  u[0] = 0.0;
99  // Set up the tridiagonal matrix to solve for the accelerations
100  for (unsigned int i = 1; i <= n-2; ++i){
101  double sig = (x_data[i]-x_data[i-1]) / (x_data[i+1]-x_data[i-1]);
102  double p = sig*yddot_data[i-1]+2.0;
103  yddot_data[i] = (sig-1.0)/p;
104  u[i] = (y_data[i+1]-y_data[i])/(x_data[i+1]-x_data[i]) -
105  (y_data[i]-y_data[i-1])/(x_data[i]-x_data[i-1]);
106  u[i] = (6.0*u[i]/(x_data[i+1]-x_data[i-1]) - sig*u[i-1])/p;
107  }
108  // Final boundary conditions
109  double qn = 0.0;
110  u[n-1] = 0.0;
111 
112  // Back substitution loop of the tridiagonal algorithm
113  yddot_data[n-1] = ( u[n-1] - qn*u[n-2]) / ( qn*yddot_data[n-2] + 1.0 );
114  for ( int k = n-2; k >= 0; --k ) {
115  yddot_data[k] = yddot_data[k]*yddot_data[k+1] + u[k];
116  }
117 }
118 
120  CubicSpline cubic;
121  apply(cubic);
122  return cubic;
123 };
124 
126 
127  spline.discretised_domain = x_data;
128  spline.cubic_polynomials.resize(x_data.size()-1); // One less polynomials than there is points
129  for (unsigned int i = 0; i < spline.cubic_polynomials.size(); ++i ) {
130  spline.cubic_polynomials[i] = CubicPolynomial::SecondDerivativeInterpolation(
131  x_data[i], y_data[i], yddot_data[i],
132  x_data[i+1], y_data[i+1], yddot_data[i+1] );
133  }
134 };
135 
136 
137 /*****************************************************************************
138 ** Implementation [DerivativeHeuristicCubicSpline]
139 *****************************************************************************/
140 
142  const double ydot_0, const double ydot_f) ecl_assert_throw_decl(StandardException) :
143  x_data(x_set),
144  y_data(y_set)
145 {
147  /*********************
148  ** Via Velocities
149  **********************/
150  ydot_data.resize(x_data.size());
151  ydot_data[0] = ydot_0;
152  for (unsigned int i = 1; i < (x_data.size()-1); ++i ) { // Skip first and last.
153  double ydot_before, ydot_after;
154  ydot_before = (y_data[i]-y_data[i-1])/(x_data[i]-x_data[i-1]);
155  ydot_after = (y_data[i+1]-y_data[i])/(x_data[i+1]-x_data[i]);
156  ydot_data[i] = (ydot_before + ydot_after)/2;
157  }
158  ydot_data[x_data.size()-1] = ydot_f;
159 }
160 
162  CubicSpline cubic;
163  apply(cubic);
164  return cubic;
165 };
166 
168 
169  spline.discretised_domain = x_data;
170  spline.cubic_polynomials.resize(x_data.size()-1); // One less polynomials than there is points
171  for (unsigned int i = 0; i < spline.cubic_polynomials.size(); ++i ) {
172  spline.cubic_polynomials[i] = CubicPolynomial::DerivativeInterpolation(
173  x_data[i], y_data[i], ydot_data[i],
174  x_data[i+1], y_data[i+1], ydot_data[i+1] );
175  }
176 };
177 
178 } // namespace blueprints
179 
180 /*****************************************************************************
181 ** Implementation [BluePrintFactory][CubicSpline]
182 *****************************************************************************/
183 
186 
188  return C2CubicSpline(x_set, y_set);
189 };
190 
192  const Array<double>& x_set, const Array<double>& y_set, const double ydot_0, const double ydot_f) {
193  return C2CubicSpline(x_set, y_set, ydot_0, ydot_f);
194 };
195 
197  const Array<double>& x_set, const Array<double>& y_set, const double ydot_0, const double ydot_f) {
198  return DerivativeHeuristicCubicSpline(x_set, y_set, ydot_0, ydot_f);
199 }
200 
201 } // namespace ecl
Array< CubicPolynomial > cubic_polynomials
Embedded control libraries.
Blueprint for generating a cubic spline satisfying C2 constraints.
#define ecl_throw(exception)
#define LOC
#define ecl_assert_throw(expression, exception)
OutOfRangeError
void apply(ecl::CubicSpline &spline) const
Apply the blueprint to configure an existing object.
ecl::CubicSpline instantiate()
Instantiate a copy of the object that is blueprinted.
Storage container for a cubic spline interpolation.
Array< double > discretised_domain
ecl::Array< double > yddot_data
static size_type size()
#define ecl_assert_throw_decl(exception)
ecl::Array< double > x_data
Blueprint for generating a cubic spline satisfying C2 constraints.
void apply(ecl::CubicSpline &spline) const
Apply the blueprint to configure an existing object.
C2CubicSpline(const ecl::Array< double > &x_set, const ecl::Array< double > &y_set) ecl_assert_throw_decl(ecl ecl::CubicSpline instantiate()
Constructor that properly configures/initialises the blueprint.
C2CubicSpline()
Default constructor.
ecl::Array< double > y_data


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