BTSpline1D Class Reference

#include <BTSpline1D.hh>

List of all members.

Public Member Functions

double operator() (double x) const
 BTSpline1D (const BTSpline1D &s)
 BTSpline1D ()
 BTSpline1D (int extent, const Data_t *nodes, const Data_t *values)
 BTSpline1D (int extent, const Data_t *nodes, double values_function(double x))
 BTSpline1D (int extent, const Data_t *nodes, const Data_t *values, Data_t slope0, Data_t slopeN)
 BTSpline1D (int extent, const Data_t *nodes, double values_function(double x), Data_t slope0, Data_t slopeN)
 BTSpline1D (int extent, const Data_t *nodes, const Data_t *values, const Data_t *gradients)
 BTSpline1D (int extent, const Data_t *nodes, void values_grads_function(double x, double *val, double *grads))
 ~BTSpline1D ()

Protected Member Functions

void captureGrid (int extent, const Data_t *nodes)
void captureValues (int extent, const Data_t *values)
void fillValues (int extent, double values_function(double x))
void computeSecondDerivs ()

Protected Attributes

int extent_
Data_t * nodePoints_
Data_t * distances_
Data_t * values_
Data_t * grads_
Data_t * secondDerivs_
bool naturalBoundaryConditions
Data_t slope0_
Data_t slopeN_

Private Types

typedef double Data_t

Detailed Description

Definition at line 107 of file BTSpline1D.hh.

Member Typedef Documentation

typedef double BTSpline1D::Data_t [private]

Definition at line 113 of file BTSpline1D.hh.

Constructor & Destructor Documentation

BTSpline1D::BTSpline1D ( const BTSpline1D &  s  ) 

Definition at line 46 of file BTSpline1D.cc.

References distances_, extent_, grads_, nodePoints_, secondDerivs_, and values_.

00046                                           : grads_(0), secondDerivs_(0) {
00047     extent_ = s.extent_;
00048     int i;
00049     nodePoints_ = new Data_t[extent_];
00050     distances_ = new Data_t[extent_];
00051     values_ = new Data_t[extent_];
00052     for (i=0; i<extent_; ++i) {
00053         nodePoints_[i] = s.nodePoints_[i];
00054     }
00055     for (i=0; i<extent_; ++i) {
00056         distances_[i] = s.distances_[i];
00057     }
00058     for (i=0; i<extent_; ++i) {
00059         values_[i] = s.values_[i];
00060     }
00061     if (s.grads_) {
00062         grads_ = new Data_t[extent_];
00063         for (i=0; i<extent_; ++i) {
00064             grads_[i] = s.grads_[i];
00065         }
00066     }
00067     if (s.secondDerivs_) {
00068         secondDerivs_ = new Data_t[extent_];
00069         for (i=0; i<extent_; ++i) {
00070             secondDerivs_[i] = s.secondDerivs_[i];
00071         }
00072     }
00073 
00074 }

BTSpline1D::BTSpline1D (  ) 

Definition at line 77 of file BTSpline1D.cc.

00077                       : extent_(0), nodePoints_(0),
00078 values_(0), grads_(0), secondDerivs_(0){}

BTSpline1D::BTSpline1D ( int  extent,
const Data_t *  nodes,
const Data_t *  values 
)

Definition at line 256 of file BTSpline1D.cc.

References captureGrid(), captureValues(), and computeSecondDerivs().

00260   : grads_(0), naturalBoundaryConditions (true) {
00261 
00262     captureGrid ( extent, nodes );
00263     captureValues ( extent, values );
00264     computeSecondDerivs ();
00265 
00266 }                                                     // BTSpline1D

BTSpline1D::BTSpline1D ( int  extent,
const Data_t *  nodes,
double   values_functiondouble x 
)

Definition at line 269 of file BTSpline1D.cc.

References captureGrid(), computeSecondDerivs(), and fillValues().

00273   : grads_(0), naturalBoundaryConditions (true) {
00274 
00275     captureGrid ( extent, nodes );
00276     fillValues (extent, values_function);
00277     computeSecondDerivs ();
00278 
00279 }                                                     // BTSpline1D

BTSpline1D::BTSpline1D ( int  extent,
const Data_t *  nodes,
const Data_t *  values,
Data_t  slope0,
Data_t  slopeN 
)

Definition at line 282 of file BTSpline1D.cc.

References captureGrid(), captureValues(), and computeSecondDerivs().

00287   : grads_(0), naturalBoundaryConditions (false),
00288 slope0_ (slope0),  slopeN_ (slopeN) {
00289 
00290     captureGrid ( extent, nodes );
00291     captureValues ( extent, values );
00292     computeSecondDerivs ();
00293 
00294 }                                                     // BTSpline1D

BTSpline1D::BTSpline1D ( int  extent,
const Data_t *  nodes,
double   values_functiondouble x,
Data_t  slope0,
Data_t  slopeN 
)

Definition at line 297 of file BTSpline1D.cc.

References captureGrid(), computeSecondDerivs(), and fillValues().

00302   : grads_(0), naturalBoundaryConditions (false),
00303 slope0_ (slope0),  slopeN_ (slopeN) {
00304 
00305     captureGrid ( extent, nodes );
00306     fillValues (extent, values_function);
00307     computeSecondDerivs ();
00308 
00309 }                                                     // BTSpline1D

BTSpline1D::BTSpline1D ( int  extent,
const Data_t *  nodes,
const Data_t *  values,
const Data_t *  gradients 
)

Definition at line 318 of file BTSpline1D.cc.

References captureGrid(), extent_, grads_, and values_.

00325   : secondDerivs_(0) {
00326 
00327     captureGrid ( extent, nodes );
00328 
00329     // Capture values and grads by copying the arrays.
00330     //-| This could be sped up by the use of memcopy,
00331     //-| but is a one-shot deal so I opt for maximal clarity.
00332 
00333     const Data_t* valp  = values;
00334     const Data_t* gradp = gradients;
00335     Data_t* vp = values_;
00336     Data_t* gp = grads_;
00337     for ( int i=0; i < extent_; i++ ) {
00338         *vp++ = *valp++;
00339         *gp++ = *gradp++;
00340     }
00341 
00342 }                                                     // BTSpline1D - separate tables for values, grads

BTSpline1D::BTSpline1D ( int  extent,
const Data_t *  nodes,
void   values_grads_functiondouble x, double *val, double *grads 
)
BTSpline1D::~BTSpline1D (  ) 

Definition at line 375 of file BTSpline1D.cc.

References distances_, grads_, nodePoints_, secondDerivs_, and values_.

00375                          {
00376 
00377     delete[] nodePoints_;
00378     delete[] distances_;
00379     delete[] values_;
00380     if (grads_) delete[] grads_;
00381     if (secondDerivs_) delete[] secondDerivs_;
00382 
00383 }                                                     // Destructor

Member Function Documentation

void BTSpline1D::captureGrid ( int  extent,
const Data_t *  nodes 
) [protected]

Definition at line 80 of file BTSpline1D.cc.

References distances_, extent_, and nodePoints_.

Referenced by BTSpline1D().

00080                                                                {
00081 
00082     // Operations common to all the constructors of BTSpline1D:  Capture the
00083     // shape and nodes of the grid, and compute differences.
00084 
00085     int i;
00086 
00087     // Capture N_, extents, and nodes
00088 
00089     extent_ = extent;
00090 
00091     nodePoints_ = new Data_t[extent];
00092     for ( i = 0; i < extent_; i++ ) {
00093         nodePoints_[i] = nodes[i];
00094     }
00095 
00096     // Compute h_ -- the differences between consecutive grid points
00097 
00098     distances_ = new Data_t[extent];
00099     for ( i = 0; i < extent_-1; i++ ) {
00100         distances_[i] = nodePoints_[i+1] - nodePoints_[i];
00101         if ( distances_[i] <= 0 ) {
00102             std::cout <<
00103                 "Spline1D initialized with nodePoints not monotonic increasing\n";
00104             std::abort();
00105         }
00106     }
00107 
00108 }                                                     /* captureGrid() */

void BTSpline1D::captureValues ( int  extent,
const Data_t *  values 
) [protected]

Definition at line 111 of file BTSpline1D.cc.

References extent_, and values_.

Referenced by BTSpline1D().

00111                                                                   {
00112     values_ = new Data_t[extent_];
00113     const Data_t* valp  = values;
00114     Data_t* vp = values_;
00115     for ( int i=0; i < extent_; i++ ) {
00116         *vp++ = *valp++;
00117     }
00118 }                                                     // captureValues

void BTSpline1D::computeSecondDerivs (  )  [protected]

Definition at line 133 of file BTSpline1D.cc.

References distances_, extent_, naturalBoundaryConditions, secondDerivs_, slope0_, slopeN_, and values_.

Referenced by BTSpline1D().

00133                                       {
00134 
00135     secondDerivs_ = new Data_t[extent_];
00136 
00137     // Based on Numerical Recipes page 114.
00138 
00139 #ifdef MATH
00140 
00141     For 1 <= j <= N-2, we have equations matching first derivatives at the
00142         node points.  These are
00143 
00144         distance[j-1]/6              * y[j-1]
00145         + (distance[j]+distance[j-1])/3 * y[j]
00146         +  distance[j]/6                * y[j+1]
00147         = dyRight/dxRight - dyLeft/dxLeft
00148 
00149         We also have, for "natural" boundary conditions, y2[0] = y2[N-1] = 0.
00150             (Or for specified slopes, a different relation, reflected in the elses.)
00151 
00152     So our equations are
00153 
00154                 M Y2 = S
00155 
00156                 where Y2 are the desired second derivatives,
00157 
00158                 S is a vector {0, (dyRight/dxRight - dyLeft/dxLeft) from 1...N-2, 0}
00159     or for slope boundaries,
00160             S[0]   = 3 * ( (y[1]-y[0])/(x[1]-x[0]) - v )         / (x[1]-x[0])      = S(v)
00161             S[N-1] = 3 * ( w - (y[N-1]-y[N-2])/(x[N-1]-x[N-2]) ) / (x[N-1]-x[N-2])  = S(w)
00162 
00163             and M is a matrix of the form
00164 
00165             1      0 or .5
00166             d0/6  (d0+d1)/3  d1/6
00167             d1/6   (d1+d2)/3  d2/6
00168             d2/6   (d2+d3)/3  d3/6
00169             ...
00170             dN-3/6 (dN-3+DN-2)/3 DN-2/3
00171             0 or .5       1
00172 
00173             Step 1 in solving this tridiagonal system is to eliminate the lower part.
00174             That only affects the first rank above the diagonal; we will call the
00175         affected rank u.  (u is a vector with indices ranging from 0 to N-2 (not N-1)).
00176         By inspection of M, u[0] starts out as 0 for our natural condition.
00177 
00178             We will also hold S in the place where we want y2 to end up.  S[0] starts out
00179             as 0 as well, for the natural condition.
00180 
00181                 We will not form u in advance, nor form the below-diagonal, because they are
00182                 simple expressions in terns of the distances.  For each row, we get a new
00183                 pivot from the original diagonal element modified by the u above it when
00184                 eliminating the below-diagonal term, and a new y2 similarly, and then we divide
00185                 thru by the pivot.
00186 #endif
00187 
00188     // Decompose the tridiagonal matrix representing the conditions of
00189     // continuous second derivative.
00190 
00191                 double* u = new double[extent_];
00192 
00193     Data_t* y = values_;                              // aliases just to shorten the expressions.
00194     Data_t* y2 = secondDerivs_;
00195     int N = extent_;
00196 
00197     if (naturalBoundaryConditions) {
00198         y2[0] = u[0] = 0;
00199     }
00200     else {
00201         y2[0] = 3 * ( (y[1]-y[0])/distances_[0] - slope0_ ) / distances_[0];
00202     }
00203 
00204     double dxLeft;
00205     double dyLeft;
00206     double dxRight;
00207     double dyRight;
00208     double pivot;
00209     for ( int i=1; i < N-1; ++i ) {
00210         dxLeft  = distances_[i-1];
00211         dxRight = distances_[i];
00212         dyLeft  = y[i]  - y[i-1];
00213         dyRight = y[i+1]- y[i];
00214         pivot   = (dxLeft+dxRight)/3.0 - dxLeft * u[i-1] / 6.0;
00215         y2[i]   = dyRight/dxRight - dyLeft/dxLeft - dxLeft * y2[i-1] / 6.0;
00216         u[i]    = dxRight / (6.0 * pivot);
00217         y2[i]   = y2[i] / pivot;
00218     }
00219 
00220     double lower;                                     // The element one to the left of the bottom right of M
00221     if (naturalBoundaryConditions) {
00222         y2[N-1] = 0;
00223         lower = 0;
00224     }
00225     else {
00226         y2[N-1] =
00227             3 * ( slopeN_ - (y[N-1]-y[N-2])/distances_[N-2] ) / distances_[N-2];
00228         lower = .5;
00229     }
00230     pivot    = 1 - lower*u[N-2];
00231     y2[N-1] -= lower * y2[N-2];
00232     y2[N-1]  = y2[N-1] / pivot;
00233 
00234     // Backsolve that decomposition to get the second dervatives.
00235 
00236     for ( int k=N-2; k >= 0; --k) {
00237         y2[k] -= u[k] * y2[k+1];
00238     }
00239 
00240     // secondDerivs_ is already y2, and we no longer need u
00241 
00242     delete[] u;
00243 
00244     return;
00245 }

void BTSpline1D::fillValues ( int  extent,
double   values_functiondouble x 
) [protected]

Definition at line 121 of file BTSpline1D.cc.

References extent_, nodePoints_, and values_.

Referenced by BTSpline1D().

00121                                                                           {
00122     values_ = new Data_t[extent_];
00123     double x;
00124     Data_t* vp = values_;
00125     for ( int i = 0; i < extent_; i++ ) {
00126         x = nodePoints_[i];
00127         *vp = values_function ( x );
00128         vp++;
00129     }
00130 }                                                     // fillValues

double BTSpline1D::operator() ( double  x  )  const

Definition at line 386 of file BTSpline1D.cc.

References distances_, extent_, grads_, nodePoints_, secondDerivs_, and values_.

00386                                                {
00387 
00388     // Find the interval that x lies within.  This has a lower coordinate
00389     // of origin.  Also compute the f0, f1, g0 and g1 functions.
00390     // (See mathematics comment at end of this file):
00391     //-| This step will point to the exactly matching point if there is one,
00392     //-| and to the first or next-to-last point if out of range.
00393 
00394     double f0;
00395     double f1;
00396     double g0;
00397     double g1;
00398 
00399     int a = 0;                                        // Highest node value known NOT to exceed x.
00400     //   Will end up as highest node which
00401     //   DOES not exceed x.
00402     int b = extent_ - 1;                              // Lowest node value which must exceed xd,
00403     //   assuming x is not outside the range.
00404 
00405     while (b != (a+1) ) {
00406         int c = (a+b+1)>>1;
00407         if (x > nodePoints_[c]) {
00408             a = c;
00409         }
00410         else {
00411             b = c;
00412         }
00413     }
00414 
00415     // Now use either the f0,f1, g0, g1 algorithm if gradients are knowen,
00416     // or the method in Numerical Recipes if second derivs have been computed:
00417 
00418     double spline;
00419 
00420     double dx;
00421     double oneMinusDx;
00422 
00423     dx = (x - nodePoints_[a])/distances_[a];
00424     oneMinusDx  = 1 - dx;
00425 
00426     if ( grads_ ) {
00427 
00428         double oneMinusX2;
00429         double x2;
00430 
00431         x2 = dx * dx;
00432         oneMinusX2 = oneMinusDx * oneMinusDx;
00433 
00434         f0 = (2. * dx + 1.) * oneMinusX2;
00435         f1 = (3. - 2. * dx) * x2;
00436         g0 =  distances_[a] * dx * oneMinusX2;
00437         g1 = -distances_[a] * oneMinusDx * x2;
00438 
00439         // Sum these F and G elements times the corresponding value (for F) and
00440         // gradient in each direction (for G) to get the answer.
00441 
00442         spline = f0 * values_[a] + f1 * values_[a+1]
00443             + g0 * grads_[a]  + g1 * grads_[a+1] ;
00444 
00445     }
00446     else {
00447         double spacingFactor = distances_[a] * distances_[a] / 6.0;
00448         double leftCubic  = (dx*dx*dx - dx) * secondDerivs_[b];
00449         double rightCubic =
00450             (oneMinusDx*oneMinusDx*oneMinusDx-oneMinusDx) * secondDerivs_[a];
00451         spline =   dx * values_[b] + oneMinusDx * values_[a]
00452             + (leftCubic + rightCubic) * spacingFactor;
00453     }
00454 
00455     return spline;
00456 
00457 }                                                     /* operator() */

Member Data Documentation

Data_t* BTSpline1D::distances_ [protected]

Definition at line 213 of file BTSpline1D.hh.

Referenced by BTSpline1D(), captureGrid(), computeSecondDerivs(), operator()(), and ~BTSpline1D().

int BTSpline1D::extent_ [protected]

Definition at line 211 of file BTSpline1D.hh.

Referenced by BTSpline1D(), captureGrid(), captureValues(), computeSecondDerivs(), fillValues(), and operator()().

Data_t* BTSpline1D::grads_ [protected]

Definition at line 215 of file BTSpline1D.hh.

Referenced by BTSpline1D(), operator()(), and ~BTSpline1D().

bool BTSpline1D::naturalBoundaryConditions [protected]

Definition at line 218 of file BTSpline1D.hh.

Referenced by computeSecondDerivs().

Data_t* BTSpline1D::nodePoints_ [protected]

Definition at line 212 of file BTSpline1D.hh.

Referenced by BTSpline1D(), captureGrid(), fillValues(), operator()(), and ~BTSpline1D().

Data_t* BTSpline1D::secondDerivs_ [protected]

Definition at line 216 of file BTSpline1D.hh.

Referenced by BTSpline1D(), computeSecondDerivs(), operator()(), and ~BTSpline1D().

Data_t BTSpline1D::slope0_ [protected]

Definition at line 219 of file BTSpline1D.hh.

Referenced by computeSecondDerivs().

Data_t BTSpline1D::slopeN_ [protected]

Definition at line 220 of file BTSpline1D.hh.

Referenced by computeSecondDerivs().

Data_t* BTSpline1D::values_ [protected]

Definition at line 214 of file BTSpline1D.hh.

Referenced by BTSpline1D(), captureValues(), computeSecondDerivs(), fillValues(), operator()(), and ~BTSpline1D().

The documentation for this class was generated from the following files:
Generated on 15 Jun 2016 for g4sim by  doxygen 1.6.1