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