Rythmos_LinearInterpolator.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_H
00030 #define Rythmos_LINEAR_INTERPOLATOR_H
00031 
00032 #include "Rythmos_InterpolatorBase.hpp"
00033 
00034 namespace Rythmos {
00035 
00036 template<class Scalar>
00037 class LinearInterpolator : virtual public InterpolatorBase<Scalar>
00038 {
00039 public:
00040 
00042   ~LinearInterpolator() {};
00043 
00045   LinearInterpolator();
00046 
00048   bool supportsCloning() const;
00049 
00051   Teuchos::RCP<InterpolatorBase<Scalar> > cloneInterpolator() const;
00052 
00054   void interpolate(
00055     const typename DataStore<Scalar>::DataStoreVector_t &data_in
00056     ,const Array<Scalar> &t_values
00057     ,typename DataStore<Scalar>::DataStoreVector_t *data_out
00058     ) const;
00059 
00061   int order() const; 
00062 
00064   std::string description() const;
00065 
00067   void describe(
00068     Teuchos::FancyOStream &out,
00069     const Teuchos::EVerbosityLevel verbLevel
00070     ) const;
00071 
00073   void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList);
00074 
00076   Teuchos::RCP<Teuchos::ParameterList> getParameterList();
00077 
00079   Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
00080 
00081 private:
00082 
00083   Teuchos::RCP<Teuchos::ParameterList> parameterList_;
00084 
00085 };
00086 
00087 template<class Scalar>
00088 LinearInterpolator<Scalar>::LinearInterpolator()
00089 {}
00090 
00091 template<class Scalar>
00092 bool LinearInterpolator<Scalar>::supportsCloning() const
00093 {
00094   return true;
00095 }
00096 
00097 template<class Scalar>
00098 Teuchos::RCP<InterpolatorBase<Scalar> >
00099 LinearInterpolator<Scalar>::cloneInterpolator() const
00100 {
00101   Teuchos::RCP<LinearInterpolator<Scalar> >
00102     interpolator = Teuchos::rcp(new LinearInterpolator<Scalar>);
00103   if (!is_null(parameterList_))
00104     interpolator->parameterList_ = parameterList(*parameterList_);
00105   return interpolator;
00106 }
00107 
00108 template<class Scalar>
00109 void LinearInterpolator<Scalar>::interpolate(
00110     const typename DataStore<Scalar>::DataStoreVector_t &data_in
00111     ,const Array<Scalar> &t_values
00112     ,typename DataStore<Scalar>::DataStoreVector_t *data_out ) const
00113 {
00114   typedef Teuchos::ScalarTraits<Scalar> ST;
00115 #ifdef TEUCHOS_DEBUG
00116   assertBaseInterpolatePreconditions(data_in,t_values,data_out);
00117 #endif // TEUCHOS_DEBUG
00118   data_out->clear();
00119   if (t_values.size() == 0) {
00120     return;
00121   }
00122   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00123   Teuchos::OSTab ostab(out,1,"LI::interpolator");
00124   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00125     *out << "data_in:" << std::endl;
00126     for (unsigned int i=0 ; i<data_in.size() ; ++i) {
00127       *out << "data_in[" << i << "] = " << std::endl;
00128       data_in[i].describe(*out,Teuchos::VERB_EXTREME);
00129     }
00130     *out << "t_values = " << std::endl;
00131     for (unsigned int i=0 ; i<t_values.size() ; ++i) {
00132       *out << "t_values[" << i << "] = " << t_values[i] << std::endl;
00133     }
00134   }
00135 
00136   if (data_in.size() == 1) {
00137     // trivial case of one node
00138     // preconditions assert that t_values[0] == data_in[0].time so we can just pass it out
00139     DataStore<Scalar> DS(data_in[0]);
00140     data_out->push_back(DS);
00141   } else {
00142     // data_in.size() >= 2
00143     int n = 0; // index into t_values
00144     for (int i=0 ; i<Teuchos::as<int>(data_in.size())-1 ; ++i) {
00145       const Scalar& ti = data_in[i].time;
00146       const Scalar& tip1 = data_in[i+1].time;
00147       const Scalar  h = tip1-ti;
00148       while ((ti <= t_values[n]) && (t_values[n] <= tip1)) {
00149         // First we check for exact node matches:
00150         if (t_values[n] == ti) {
00151           DataStore<Scalar> DS(data_in[i]);
00152           data_out->push_back(DS);
00153         } else if (t_values[n] == tip1) {
00154           DataStore<Scalar> DS(data_in[i+1]);
00155           data_out->push_back(DS);
00156         } else {
00157           const Scalar& t = t_values[n];
00158           Teuchos::RCP<const Thyra::VectorBase<Scalar> > xi =      data_in[i  ].x;
00159           Teuchos::RCP<const Thyra::VectorBase<Scalar> > xip1 =    data_in[i+1].x;
00160           Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdoti =   data_in[i  ].xdot;
00161           Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotip1 = data_in[i+1].xdot;
00162         
00163           // interpolate this point
00164           //
00165           // x(t) = (t-ti)/(tip1-ti) * xip1 + (1-(t-ti)/(tip1-ti)) * xi
00166           //
00167           // Above, it is easy to see that:
00168           //
00169           //    x(ti) = xi
00170           //    x(tip1) = xip1
00171           //
00172           DataStore<Scalar> DS;
00173           DS.time = t;
00174           const Scalar dt = t-ti;
00175           const Scalar dt_over_h = dt / h;
00176           const Scalar one_minus_dt_over_h = ST::one() - dt_over_h;
00177           // x
00178           RCP<Thyra::VectorBase<Scalar> > x = createMember(xi->space());
00179           Thyra::V_StVpStV(&*x,dt_over_h,*xip1,one_minus_dt_over_h,*xi);
00180           DS.x = x;
00181           // xdot
00182           RCP<Thyra::VectorBase<Scalar> > xdot;
00183           if ((xdoti != Teuchos::null) && (xdotip1 != Teuchos::null)) {
00184             xdot = createMember(xdoti->space());
00185             Thyra::V_StVpStV(&*xdot,dt_over_h,*xdotip1,one_minus_dt_over_h,*xdoti);
00186           }
00187           DS.xdot = xdot;
00188           // And finally we estimate our order of accuracy
00189           DS.accuracy = h;
00190           // Push DataStore object onto vector:
00191           data_out->push_back(DS);
00192           // Move to the next time to consider!
00193         }
00194         n++;
00195         if (n == Teuchos::as<int>(t_values.size())) {
00196           return;
00197         }
00198       }
00199     }
00200   } // data_in.size() == 1
00201 }
00202 
00203 template<class Scalar>
00204 int LinearInterpolator<Scalar>::order() const
00205 {
00206   return(1);
00207 }
00208 
00209 template<class Scalar>
00210 std::string LinearInterpolator<Scalar>::description() const
00211 {
00212   std::string name = "Rythmos::LinearInterpolator";
00213   return(name);
00214 }
00215 
00216 template<class Scalar>
00217 void LinearInterpolator<Scalar>::describe(
00218       Teuchos::FancyOStream       &out
00219       ,const Teuchos::EVerbosityLevel      verbLevel
00220       ) const
00221 {
00222   if ( (Teuchos::as<int>(verbLevel) == Teuchos::as<int>(Teuchos::VERB_DEFAULT) ) ||
00223        (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW)     )
00224      ) {
00225     out << description() << "::describe" << std::endl;
00226   } else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW)) {
00227   } else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_MEDIUM)) {
00228   } else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH)) {
00229   }
00230 }
00231 
00232 template <class Scalar>
00233 void LinearInterpolator<Scalar>::setParameterList(
00234   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00235   )
00236 {
00237   parameterList_ = paramList;
00238   int outputLevel = parameterList_->get( "outputLevel", int(-1) );
00239   outputLevel = std::min(std::max(outputLevel,-1),4);
00240   this->setVerbLevel(static_cast<Teuchos::EVerbosityLevel>(outputLevel));
00241   // 2007/05/18: rabartl: ToDo: Replace with standard "Verbose Object"
00242   // sublist! and validate the sublist!
00243 
00244 }
00245 
00246 template <class Scalar>
00247 Teuchos::RCP<Teuchos::ParameterList> LinearInterpolator<Scalar>::getParameterList()
00248 {
00249   return(parameterList_);
00250 }
00251 
00252 template <class Scalar>
00253 Teuchos::RCP<Teuchos::ParameterList> LinearInterpolator<Scalar>::unsetParameterList()
00254 {
00255   Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00256   parameterList_ = Teuchos::null;
00257   return(temp_param_list);
00258 }
00259 
00260 } // namespace Rythmos
00261 
00262 #endif // Rythmos_LINEAR_INTERPOLATOR_H

Generated on Tue Oct 20 12:46:08 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7