Rythmos_CubicSplineInterpolator_def.hpp

00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_CUBIC_SPLINE_INTERPOLATOR_DEF_H
00030 #define Rythmos_CUBIC_SPLINE_INTERPOLATOR_DEF_H
00031 
00032 #include "Rythmos_CubicSplineInterpolator_decl.hpp"
00033 #include "Rythmos_InterpolatorBaseHelpers.hpp"
00034 #include "Thyra_VectorStdOps.hpp"
00035 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00036 
00037 namespace Rythmos {
00038 
00039 // non-member constructor
00040 template<class Scalar>
00041 RCP<CubicSplineInterpolator<Scalar> > cubicSplineInterpolator()
00042 {
00043   RCP<CubicSplineInterpolator<Scalar> > csi = Teuchos::rcp(new CubicSplineInterpolator<Scalar>() );
00044   return csi;
00045 }
00046 
00047 template<class Scalar>
00048 void computeCubicSplineCoeff(
00049     const typename DataStore<Scalar>::DataStoreVector_t & data,
00050     const Ptr<CubicSplineCoeff<Scalar> > & coeffPtr
00051     )
00052 {
00053   using Teuchos::outArg;
00054   typedef Teuchos::ScalarTraits<Scalar> ST;
00055   TEST_FOR_EXCEPTION( 
00056       (data.size() < 2), std::logic_error,
00057       "Error!  A minimum of two data points is required for this cubic spline."
00058       );
00059   // time data in the DataStoreVector should be unique and sorted 
00060   Array<Scalar> t;
00061   Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
00062   dataStoreVectorToVector<Scalar>( data, &t, &x_vec, &xdot_vec, NULL );
00063 #ifdef RYTHMOS_DEBUG
00064   assertTimePointsAreSorted<Scalar>( t );
00065 #endif // RYTHMOS_DEBUG
00066   // 11/18/08 tscoffe:  Question:  Should I erase everything in coeffPtr or
00067   // re-use what I can?  For now, I'll erase and create new each time.
00068   CubicSplineCoeff<Scalar>& coeff = *coeffPtr;
00069   // If there are only two points, then we do something special and just create
00070   // a linear polynomial between the points and return.
00071   if (t.size() == 2) {
00072     coeff.t.clear(); 
00073     coeff.a.clear(); coeff.b.clear(); coeff.c.clear(); coeff.d.clear();
00074     coeff.t.reserve(2); 
00075     coeff.a.reserve(1); coeff.b.reserve(1); coeff.c.reserve(1); coeff.d.reserve(1);
00076     coeff.t.push_back(t[0]);
00077     coeff.t.push_back(t[1]);
00078     coeff.a.push_back(x_vec[0]->clone_v());
00079     coeff.b.push_back(Thyra::createMember(x_vec[0]->space()));
00080     coeff.c.push_back(Thyra::createMember(x_vec[0]->space()));
00081     coeff.d.push_back(Thyra::createMember(x_vec[0]->space()));
00082     Scalar h = coeff.t[1] - coeff.t[0];
00083     V_StVpStV(outArg(*coeff.b[0]),ST::one()/h,*x_vec[1],-ST::one()/h,*x_vec[0]);
00084     V_S(outArg(*coeff.c[0]),ST::zero());
00085     V_S(outArg(*coeff.d[0]),ST::zero());
00086     return;
00087   }
00088   // Data objects we'll need:
00089   int n = t.length()-1; // Number of intervals
00090   coeff.t.clear(); coeff.t.reserve(n+1);
00091   coeff.a.clear(); coeff.a.reserve(n+1);
00092   coeff.b.clear(); coeff.b.reserve(n);
00093   coeff.c.clear(); coeff.c.reserve(n+1);
00094   coeff.d.clear(); coeff.d.reserve(n);
00095 
00096   Array<Scalar> h(n);
00097   Array<RCP<Thyra::VectorBase<Scalar> > > alpha(n), z(n+1);
00098   Array<Scalar> l(n+1), mu(n);
00099   for (int i=0 ; i<n ; ++i) {
00100     coeff.t.push_back(t[i]);
00101     coeff.a.push_back(x_vec[i]->clone_v());
00102     coeff.b.push_back(Thyra::createMember(x_vec[0]->space()));
00103     coeff.c.push_back(Thyra::createMember(x_vec[0]->space()));
00104     coeff.d.push_back(Thyra::createMember(x_vec[0]->space()));
00105     alpha[i] = Thyra::createMember(x_vec[0]->space());
00106     z[i] = Thyra::createMember(x_vec[0]->space());
00107   }
00108   coeff.a.push_back(x_vec[n]->clone_v());
00109   coeff.t.push_back(t[n]);
00110   coeff.c.push_back(Thyra::createMember(x_vec[0]->space()));
00111   z[n] = Thyra::createMember(x_vec[0]->space());
00112   Scalar zero = ST::zero();
00113   Scalar one = ST::one();
00114   Scalar two = Scalar(2*ST::one());
00115   Scalar three = Scalar(3*ST::one());
00116 
00117   // Algorithm starts here:
00118   for (int i=0 ; i<n ; ++i) {
00119     h[i] = coeff.t[i+1]-coeff.t[i];
00120   }
00121   for (int i=1 ; i<n ; ++i) {
00122     V_StVpStV(outArg(*(alpha[i])),three/h[i],*coeff.a[i+1],-3/h[i],*coeff.a[i]);
00123     Vp_StV(outArg(*(alpha[i])),-three/h[i-1],*coeff.a[i]);
00124     Vp_StV(outArg(*(alpha[i])),+three/h[i-1],*coeff.a[i-1]);
00125   }
00126   l[0] = one;
00127   mu[0] = zero;
00128   V_S(outArg(*(z[0])),zero);
00129   for (int i=1 ; i<n ; ++i) {
00130     l[i] = 2*(coeff.t[i+1]-coeff.t[i-1])-h[i-1]*mu[i-1];
00131     mu[i] = h[i]/l[i];
00132     V_StVpStV(outArg(*(z[i])),one/l[i],*alpha[i],-h[i-1]/l[i],*z[i-1]);
00133   }
00134   l[n] = one;
00135   V_S(outArg(*(z[n])),zero);
00136   V_S(outArg(*(coeff.c[n])),zero);
00137   for (int j=n-1 ; j >= 0 ; --j) {
00138     V_StVpStV(outArg(*(coeff.c[j])),one,*z[j],-mu[j],*coeff.c[j+1]);
00139     V_StVpStV(outArg(*(coeff.b[j])),one/h[j],*coeff.a[j+1],-one/h[j],*coeff.a[j]);
00140     Vp_StV(outArg(*(coeff.b[j])),-h[j]/three,*coeff.c[j+1]);
00141     Vp_StV(outArg(*(coeff.b[j])),-h[j]*two/three,*coeff.c[j]);
00142     V_StVpStV(outArg(*(coeff.d[j])),one/(three*h[j]),*coeff.c[j+1],-one/(three*h[j]),*coeff.c[j]);
00143   }
00144   // Pop the last entry off of a and c to make them the right size.
00145   coeff.a.erase(coeff.a.end()-1);
00146   coeff.c.erase(coeff.c.end()-1);
00147 }
00148 
00149 template<class Scalar>
00150 void validateCubicSplineCoeff(const CubicSplineCoeff<Scalar>& coeff) 
00151 {
00152   int t_n = coeff.t.size();
00153   int a_n = coeff.a.size();
00154   int b_n = coeff.b.size();
00155   int c_n = coeff.c.size();
00156   int d_n = coeff.d.size();
00157   TEST_FOR_EXCEPTION( 
00158       ((a_n != t_n-1) || (a_n != b_n) || (a_n != c_n) || (a_n != d_n)),
00159       std::logic_error,
00160       "Error!  The sizes of the data structures in the CubicSplineCoeff object do not match"
00161       );
00162 }
00163 
00164 // Non-member helper function
00165 // Evaluate cubic spline 
00166 template<class Scalar>
00167 void evaluateCubicSpline(
00168     const CubicSplineCoeff<Scalar>& coeff,
00169     unsigned int j, 
00170     const Scalar& t,
00171     const Ptr<Thyra::VectorBase<Scalar> >& S,
00172     const Ptr<Thyra::VectorBase<Scalar> >& Sp, 
00173     const Ptr<Thyra::VectorBase<Scalar> >& Spp
00174     )
00175 {
00176   using Teuchos::outArg;
00177   typedef Teuchos::ScalarTraits<Scalar> ST;
00178   // Assert preconditions:
00179   validateCubicSplineCoeff<Scalar>(coeff);
00180   TEST_FOR_EXCEPTION( ((j < 0) || (j >= coeff.a.size())), std::out_of_range, "Error!, j is out of range" );
00181 
00182   Scalar dt = t-coeff.t[j];
00183   const Thyra::VectorBase<Scalar>& a = *(coeff.a[j]);
00184   const Thyra::VectorBase<Scalar>& b = *(coeff.b[j]);
00185   const Thyra::VectorBase<Scalar>& c = *(coeff.c[j]);
00186   const Thyra::VectorBase<Scalar>& d = *(coeff.d[j]);
00187   
00188   if (!Teuchos::is_null(S)) {
00189     // Evaluate S:
00190     //*S = (a) + (b)*dt + (c)*dt*dt + (d)*dt*dt*dt;
00191     V_StVpStV(outArg(*S),dt*dt*dt,d,dt*dt,c);
00192     Vp_StV(outArg(*S),dt,b);
00193     Vp_StV(outArg(*S),ST::one(),a);
00194   } 
00195   if (!Teuchos::is_null(Sp)) {
00196     // Evaluate S':
00197     //*Sp = (b) + (c)*2*dt + (d)*3*dt*dt;
00198     V_StVpStV(outArg(*Sp),Scalar(3*ST::one())*dt*dt,d,Scalar(2*ST::one())*dt,c);
00199     Vp_StV(outArg(*Sp),ST::one(),b);
00200   }
00201   if (!Teuchos::is_null(Spp)) {
00202     // Evaluate S'':
00203     //*Spp = (c)*2 + (d)*6*dt;
00204     V_StVpStV(outArg(*Spp),Scalar(6*ST::one())*dt,d,Scalar(2*ST::one()),c);
00205   }
00206 }
00207 
00208 template<class Scalar>
00209 CubicSplineInterpolator<Scalar>::CubicSplineInterpolator()
00210 {
00211   splineCoeffComputed_ = false;
00212   nodesSet_ = false;
00213 }
00214 
00215 
00216 template<class Scalar>
00217 bool CubicSplineInterpolator<Scalar>::supportsCloning() const
00218 {
00219   return true;
00220 }
00221 
00222 
00223 template<class Scalar>
00224 RCP<InterpolatorBase<Scalar> >
00225 CubicSplineInterpolator<Scalar>::cloneInterpolator() const
00226 {
00227   RCP<CubicSplineInterpolator<Scalar> >
00228     interpolator = Teuchos::rcp(new CubicSplineInterpolator<Scalar>);
00229   if (!is_null(parameterList_))
00230     interpolator->parameterList_ = parameterList(*parameterList_);
00231   return interpolator;
00232 }
00233 
00234 template<class Scalar>
00235 void CubicSplineInterpolator<Scalar>::setNodes(
00236     const RCP<const typename DataStore<Scalar>::DataStoreVector_t> & nodesPtr
00237     )
00238 {
00239   nodes_ = nodesPtr;
00240   nodesSet_ = true;
00241   splineCoeffComputed_ = false;
00242 #ifdef RYTHMOS_DEBUG
00243   const typename DataStore<Scalar>::DataStoreVector_t & nodes = *nodesPtr;
00244   // Copy nodes to internal data structure for verification upon calls to interpolate
00245   nodes_copy_ = Teuchos::rcp(new typename DataStore<Scalar>::DataStoreVector_t);
00246   nodes_copy_->reserve(nodes.size());
00247   for (int i=0 ; i<Teuchos::as<int>(nodes.size()) ; ++i) {
00248     nodes_copy_->push_back(*nodes[i].clone());
00249   }
00250 #endif // RYTHMOS_DEBUG
00251 }
00252 
00253 template<class Scalar>
00254 void CubicSplineInterpolator<Scalar>::interpolate(
00255   const Array<Scalar> &t_values,
00256   typename DataStore<Scalar>::DataStoreVector_t *data_out
00257   ) const
00258 {
00259   using Teuchos::as;
00260   using Teuchos::outArg;
00261   typedef Teuchos::ScalarTraits<Scalar> ST;
00262 
00263   TEST_FOR_EXCEPTION( nodesSet_ == false, std::logic_error,
00264       "Error!, setNodes must be called before interpolate"
00265       );
00266 #ifdef RYTHMOS_DEBUG
00267   // Check that our nodes_ have not changed between the call to setNodes and interpolate
00268   assertNodesUnChanged<Scalar>(*nodes_,*nodes_copy_);
00269   // Assert that the base interpolator preconditions are satisfied
00270   assertBaseInterpolatePreconditions(*nodes_,t_values,data_out);
00271 #endif // RYTHMOS_DEBUG
00272   
00273   // Output info
00274   const RCP<FancyOStream> out = this->getOStream();
00275   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel(); 
00276   Teuchos::OSTab ostab(out,1,"CSI::interpolator");
00277   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00278     *out << "nodes_:" << std::endl;
00279     for (unsigned int i=0 ; i<(*nodes_).size() ; ++i) {
00280       *out << "nodes_[" << i << "] = " << std::endl;
00281       (*nodes_)[i].describe(*out,Teuchos::VERB_EXTREME);
00282     }
00283     *out << "t_values = " << std::endl;
00284     for (unsigned int i=0 ; i<t_values.size() ; ++i) {
00285       *out << "t_values[" << i << "] = " << t_values[i] << std::endl;
00286     }
00287   }
00288 
00289   data_out->clear();
00290 
00291   // Return immediately if no time points are requested ...
00292   if (t_values.size() == 0) {
00293     return;
00294   }
00295 
00296   if ((*nodes_).size() == 1) {
00297     // trivial case of one node.  Preconditions assert that t_values[0] ==
00298     // (*nodes_)[0].time so we can just pass it out
00299     DataStore<Scalar> DS((*nodes_)[0]);
00300     data_out->push_back(DS);
00301   }
00302   else { // (*nodes_).size() >= 2
00303     int n = 0; // index into t_values
00304     // Loop through all of the time interpolation points in the buffer and
00305     // satisfiy all of the requested time points that you find.  NOTE: The
00306     // loop will be existed once all of the time points are satisified (see
00307     // return below).
00308     for (int i=0 ; i < as<int>((*nodes_).size())-1; ++i) {
00309       const Scalar& ti = (*nodes_)[i].time;
00310       const Scalar& tip1 = (*nodes_)[i+1].time;
00311       const TimeRange<Scalar> range_i(ti,tip1);
00312       // For the interpolation range of [ti,tip1], satisify all of the
00313       // requested points in this range.
00314       while ( range_i.isInRange(t_values[n]) ) {
00315         // First we check for exact node matches:
00316         if (compareTimeValues(t_values[n],ti)==0) {
00317           DataStore<Scalar> DS((*nodes_)[i]);
00318           data_out->push_back(DS);
00319         }
00320         else if (compareTimeValues(t_values[n],tip1)==0) {
00321           DataStore<Scalar> DS((*nodes_)[i+1]);
00322           data_out->push_back(DS);
00323         } else {
00324           if (!splineCoeffComputed_) {
00325             computeCubicSplineCoeff<Scalar>(*nodes_,outArg(splineCoeff_));
00326             splineCoeffComputed_ = true;
00327           }
00328           DataStore<Scalar> DS;
00329           RCP<Thyra::VectorBase<Scalar> > x = createMember((*nodes_)[i].x->space());
00330           evaluateCubicSpline<Scalar>( splineCoeff_, i, t_values[n], outArg(*x) );
00331           DS.time = t_values[n];
00332           DS.x = x;
00333           DS.accuracy = ST::zero();
00334           data_out->push_back(DS);
00335         }
00336         // Move to the next user time point to consider!
00337         n++;
00338         if (n == as<int>(t_values.size())) {
00339           // WE ARE ALL DONE!  MOVE OUT!
00340           return;
00341         }
00342       }
00343       // Move on the the next interpolation time range
00344     }
00345   } // (*nodes_).size() == 1
00346 }
00347 
00348 
00349 template<class Scalar>
00350 int CubicSplineInterpolator<Scalar>::order() const
00351 {
00352   return(1);
00353 }
00354 
00355 
00356 template<class Scalar>
00357 std::string CubicSplineInterpolator<Scalar>::description() const
00358 {
00359   std::string name = "Rythmos::CubicSplineInterpolator";
00360   return(name);
00361 }
00362 
00363 
00364 template<class Scalar>
00365 void CubicSplineInterpolator<Scalar>::describe(
00366   FancyOStream &out,
00367   const Teuchos::EVerbosityLevel verbLevel
00368   ) const
00369 {
00370   using Teuchos::as;
00371   if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
00372        (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)     )
00373      )
00374   {
00375     out << description() << "::describe" << std::endl;
00376   }
00377   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW))
00378   {}
00379   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM))
00380   {}
00381   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH))
00382   {}
00383 }
00384 
00385 
00386 template <class Scalar>
00387 void CubicSplineInterpolator<Scalar>::setParameterList(
00388   RCP<ParameterList> const& paramList
00389   )
00390 {
00391   TEST_FOR_EXCEPT(is_null(paramList));
00392   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00393   parameterList_ = paramList;
00394   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00395 }
00396 
00397 
00398 template <class Scalar>
00399 RCP<ParameterList>
00400 CubicSplineInterpolator<Scalar>::getNonconstParameterList()
00401 {
00402   return(parameterList_);
00403 }
00404 
00405 
00406 template <class Scalar>
00407 RCP<ParameterList>
00408 CubicSplineInterpolator<Scalar>::unsetParameterList()
00409 {
00410   RCP<ParameterList> temp_param_list;
00411   std::swap( temp_param_list, parameterList_ );
00412   return(temp_param_list);
00413 }
00414 
00415 template<class Scalar>
00416 RCP<const Teuchos::ParameterList> CubicSplineInterpolator<Scalar>::getValidParameters() const
00417 {
00418   static RCP<Teuchos::ParameterList> validPL;
00419   if (is_null(validPL)) {
00420     RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00421     Teuchos::setupVerboseObjectSublist(&*pl);
00422     validPL = pl;
00423   }
00424   return (validPL);
00425 }
00426 
00427 // 
00428 // Explicit Instantiation macro
00429 //
00430 // Must be expanded from within the Rythmos namespace!
00431 //
00432 
00433 #define RYTHMOS_CUBIC_SPLINE_INTERPOLATOR_INSTANT(SCALAR) \
00434   \
00435   template class CubicSplineInterpolator< SCALAR >; \
00436   \
00437   template class CubicSplineCoeff< SCALAR >; \
00438   template RCP<CubicSplineInterpolator< SCALAR > > cubicSplineInterpolator(); \
00439   template void computeCubicSplineCoeff( \
00440       const DataStore< SCALAR >::DataStoreVector_t & data, \
00441       const Ptr<CubicSplineCoeff< SCALAR > > & coeffPtr \
00442       ); \
00443   template void validateCubicSplineCoeff(const CubicSplineCoeff< SCALAR >& coeff); \
00444   template void evaluateCubicSpline( \
00445       const CubicSplineCoeff< SCALAR >& coeff, \
00446       unsigned int j,  \
00447       const  SCALAR & t, \
00448       const Ptr<Thyra::VectorBase< SCALAR > >& S, \
00449       const Ptr<Thyra::VectorBase< SCALAR > >& Sp,  \
00450       const Ptr<Thyra::VectorBase< SCALAR > >& Spp \
00451       ); 
00452 
00453 
00454 
00455 } // namespace Rythmos
00456 
00457 
00458 #endif // Rythmos_CUBIC_SPLINE_INTERPOLATOR_DEF_H

Generated on Wed May 12 21:25:42 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7