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

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