Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_LinearInterpolator_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_LINEAR_INTERPOLATOR_DEF_H
00030 #define Rythmos_LINEAR_INTERPOLATOR_DEF_H
00031 
00032 #include "Rythmos_LinearInterpolator_decl.hpp"
00033 #include "Rythmos_InterpolatorBaseHelpers.hpp"
00034 #include "Thyra_VectorStdOps.hpp"
00035 #include "Thyra_VectorSpaceBase.hpp"
00036 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00037 
00038 
00039 namespace Rythmos {
00040 
00041 
00042 // non-member constructor
00043 template<class Scalar>
00044 RCP<LinearInterpolator<Scalar> > linearInterpolator()
00045 {
00046   RCP<LinearInterpolator<Scalar> > li = rcp(new LinearInterpolator<Scalar>() );
00047   return li;
00048 }
00049 
00050 template<class Scalar>
00051 LinearInterpolator<Scalar>::LinearInterpolator()
00052 {
00053   nodes_ = Teuchos::null;
00054   parameterList_ = Teuchos::null;
00055 }
00056 
00057 
00058 template<class Scalar>
00059 bool LinearInterpolator<Scalar>::supportsCloning() const
00060 {
00061   return true;
00062 }
00063 
00064 
00065 template<class Scalar>
00066 RCP<InterpolatorBase<Scalar> >
00067 LinearInterpolator<Scalar>::cloneInterpolator() const
00068 {
00069   RCP<LinearInterpolator<Scalar> >
00070     interpolator = Teuchos::rcp(new LinearInterpolator<Scalar>);
00071   if (!is_null(parameterList_))
00072     interpolator->parameterList_ = parameterList(*parameterList_);
00073   return interpolator;
00074 }
00075 
00076 template<class Scalar>
00077 void LinearInterpolator<Scalar>::setNodes(
00078     const RCP<const typename DataStore<Scalar>::DataStoreVector_t> & nodes
00079     )
00080 {
00081   nodes_ = nodes;
00082 }
00083 
00084 
00085 template<class Scalar>
00086 void LinearInterpolator<Scalar>::interpolate(
00087   const Array<Scalar> &t_values,
00088   typename DataStore<Scalar>::DataStoreVector_t *data_out
00089   ) const
00090 {
00091 
00092   using Teuchos::as;
00093   typedef Teuchos::ScalarTraits<Scalar> ST;
00094 
00095 #ifdef RYTHMOS_DEBUG
00096   assertBaseInterpolatePreconditions(*nodes_,t_values,data_out);
00097 #endif // RYTHMOS_DEBUG
00098   
00099   // Output info
00100   const RCP<FancyOStream> out = this->getOStream();
00101   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel(); 
00102   Teuchos::OSTab ostab(out,1,"LI::interpolator");
00103   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00104     *out << "nodes_:" << std::endl;
00105     for (Teuchos::Ordinal i=0 ; i<(*nodes_).size() ; ++i) {
00106       *out << "nodes_[" << i << "] = " << std::endl;
00107       (*nodes_)[i].describe(*out,Teuchos::VERB_EXTREME);
00108     }
00109     *out << "t_values = " << std::endl;
00110     for (Teuchos::Ordinal i=0 ; i<t_values.size() ; ++i) {
00111       *out << "t_values[" << i << "] = " << t_values[i] << std::endl;
00112     }
00113   }
00114 
00115   data_out->clear();
00116 
00117   // Return immediatly if not time points are requested ...
00118   if (t_values.size() == 0) {
00119     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00120       *out << "Returning because no time points were requested" << std::endl;
00121     }
00122     return;
00123   }
00124 
00125   if ((*nodes_).size() == 1) {
00126     // trivial case of one node.  Preconditions assert that t_values[0] ==
00127     // (*nodes_)[0].time so we can just pass it out
00128     DataStore<Scalar> DS((*nodes_)[0]);
00129     data_out->push_back(DS);
00130     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00131       *out << "Only a single node is in the buffer, so preconditions assert that this must be the point requested" << std::endl;
00132     }
00133   }
00134   else { // (*nodes_).size() >= 2
00135     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00136       *out << "More than two nodes, looping through the intervals looking for points to interpolate" << std::endl;
00137     }
00138     int n = 0; // index into t_values
00139     // Loop through all of the time interpolation points in the buffer and
00140     // satisfiy all of the requested time points that you find.  NOTE: The
00141     // loop will be existed once all of the time points are satisified (see
00142     // return below).
00143     for (int i=0 ; i < as<int>((*nodes_).size())-1; ++i) {
00144       if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00145         *out << "Looking for interval containing: t_values["<<n<<"] = " << t_values[n] << std::endl;
00146       }
00147       const Scalar& ti = (*nodes_)[i].time;
00148       const Scalar& tip1 = (*nodes_)[i+1].time;
00149       const Scalar  h = tip1-ti;
00150       const TimeRange<Scalar> range_i(ti,tip1);
00151       // For the interpolation range of [ti,tip1], satisify all of the
00152       // requested points in this range.
00153       while ( range_i.isInRange(t_values[n]) ) {
00154         // First we check for exact node matches:
00155         if (compareTimeValues(t_values[n],ti)==0) {
00156           DataStore<Scalar> DS((*nodes_)[i]);
00157           data_out->push_back(DS);
00158           if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00159             *out << "Found an exact node match (on left), shallow copying." << std::endl;
00160             *out << "Found t_values["<<n<<"] = " << t_values[n] << " on boundary of interval ["<<ti<<","<<tip1<<"]" << std::endl;
00161           }
00162         }
00163         else if (compareTimeValues(t_values[n],tip1)==0) {
00164           DataStore<Scalar> DS((*nodes_)[i+1]);
00165           data_out->push_back(DS);
00166           if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00167             *out << "Found an exact node match (on right), shallow copying." << std::endl;
00168             *out << "Found t_values["<<n<<"] = " << t_values[n] << " on boundary of interval ["<<ti<<","<<tip1<<"]" << std::endl;
00169           }
00170         }
00171         else {
00172           if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00173             *out << "Interpolating a point (creating a new vector)..." << std::endl;
00174             *out << "Found t_values["<<n<<"] = " << t_values[n] << " in interior of interval ["<<ti<<","<<tip1<<"]" << std::endl;
00175           }
00176           // interpolate this point
00177           //
00178           // x(t) = (t-ti)/(tip1-ti) * xip1 + (1-(t-ti)/(tip1-ti)) * xi
00179           //
00180           // Above, it is easy to see that:
00181           //
00182           //    x(ti) = xi
00183           //    x(tip1) = xip1
00184           //
00185           DataStore<Scalar> DS;
00186           const Scalar& t = t_values[n];
00187           DS.time = t;
00188           // Get the time and interpolation node points
00189           RCP<const Thyra::VectorBase<Scalar> > xi = (*nodes_)[i].x;
00190           RCP<const Thyra::VectorBase<Scalar> > xip1 = (*nodes_)[i+1].x;
00191           RCP<const Thyra::VectorBase<Scalar> > xdoti = (*nodes_)[i].xdot;
00192           RCP<const Thyra::VectorBase<Scalar> > xdotip1 = (*nodes_)[i+1].xdot;
00193           // Get constants used in interplation
00194           const Scalar dt = t-ti;
00195           const Scalar dt_over_h = dt / h;
00196           const Scalar one_minus_dt_over_h = ST::one() - dt_over_h;
00197           // x = dt/h * xip1 + (1-dt/h) * xi
00198           RCP<Thyra::VectorBase<Scalar> > x;
00199           if (!is_null(xi) && !is_null(xip1)) {
00200             x = createMember(xi->space());
00201             Thyra::V_StVpStV(x.ptr(),dt_over_h,*xip1,one_minus_dt_over_h,*xi);
00202           }
00203           DS.x = x;
00204           // x = dt/h * xdotip1 + (1-dt/h) * xdoti
00205           RCP<Thyra::VectorBase<Scalar> > xdot;
00206           if (!is_null(xdoti) && !is_null(xdotip1)) {
00207             xdot = createMember(xdoti->space());
00208             Thyra::V_StVpStV(xdot.ptr(),dt_over_h,*xdotip1,one_minus_dt_over_h,*xdoti);
00209           }
00210           DS.xdot = xdot;
00211           // Estimate our accuracy ???
00212           DS.accuracy = h;
00213           // 2007/12/06: rabartl: Above, should the be a relative value of
00214           // some type.  What does this really mean?
00215           // Store this interplation
00216           data_out->push_back(DS);
00217         }
00218         // Move to the next user time point to consider!
00219         n++;
00220         if (n == as<int>(t_values.size())) {
00221           // WE ARE ALL DONE!  MOVE OUT!
00222           return;
00223         }
00224       }
00225       // Move on the the next interpolation time range
00226     }
00227   } // (*nodes_).size() == 1
00228 }
00229 
00230 
00231 template<class Scalar>
00232 int LinearInterpolator<Scalar>::order() const
00233 {
00234   return(1);
00235 }
00236 
00237 
00238 template<class Scalar>
00239 std::string LinearInterpolator<Scalar>::description() const
00240 {
00241   std::string name = "Rythmos::LinearInterpolator";
00242   return(name);
00243 }
00244 
00245 
00246 template<class Scalar>
00247 void LinearInterpolator<Scalar>::describe(
00248   FancyOStream &out,
00249   const Teuchos::EVerbosityLevel verbLevel
00250   ) const
00251 {
00252   using Teuchos::as;
00253   if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
00254        (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)     )
00255      )
00256   {
00257     out << description() << "::describe" << std::endl;
00258   }
00259   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW))
00260   {}
00261   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM))
00262   {}
00263   else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH))
00264   {}
00265 }
00266 
00267 
00268 template <class Scalar>
00269 void LinearInterpolator<Scalar>::setParameterList(
00270   RCP<ParameterList> const& paramList
00271   )
00272 {
00273   TEST_FOR_EXCEPT(is_null(paramList));
00274   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00275   parameterList_ = paramList;
00276   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00277 }
00278 
00279 
00280 template <class Scalar>
00281 RCP<ParameterList>
00282 LinearInterpolator<Scalar>::getNonconstParameterList()
00283 {
00284   return(parameterList_);
00285 }
00286 
00287 
00288 template <class Scalar>
00289 RCP<ParameterList>
00290 LinearInterpolator<Scalar>::unsetParameterList()
00291 {
00292   RCP<ParameterList> temp_param_list;
00293   std::swap( temp_param_list, parameterList_ );
00294   return(temp_param_list);
00295 }
00296 
00297 template<class Scalar>
00298 RCP<const Teuchos::ParameterList> LinearInterpolator<Scalar>::getValidParameters() const
00299 {
00300   static RCP<Teuchos::ParameterList> validPL;
00301   if (is_null(validPL)) {
00302     RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00303     Teuchos::setupVerboseObjectSublist(&*pl);
00304     validPL = pl;
00305   }
00306   return (validPL);
00307 }
00308 
00309 // 
00310 // Explicit Instantiation macro
00311 //
00312 // Must be expanded from within the Rythmos namespace!
00313 //
00314 
00315 #define RYTHMOS_LINEAR_INTERPOLATOR_INSTANT(SCALAR) \
00316   \
00317   template class LinearInterpolator< SCALAR >; \
00318   \
00319   template RCP<LinearInterpolator< SCALAR > > linearInterpolator(); 
00320 
00321 } // namespace Rythmos
00322 
00323 
00324 #endif // Rythmos_LINEAR_INTERPOLATOR_DEF_H
 All Classes Functions Variables Typedefs Friends