Rythmos_HermiteInterpolator.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_HERMITE_INTERPOLATOR_H
00030 #define Rythmos_HERMITE_INTERPOLATOR_H
00031 
00032 #include "Rythmos_InterpolatorBase.hpp"
00033 
00034 namespace Rythmos {
00052 template<class Scalar>
00053 class HermiteInterpolator : virtual public InterpolatorBase<Scalar>
00054 {
00055   public:
00056 
00058     ~HermiteInterpolator() {};
00059 
00061     HermiteInterpolator();
00062 
00064 
00071     void interpolate(
00072         const typename DataStore<Scalar>::DataStoreVector_t &data_in
00073         ,const Array<Scalar> &t_values
00074         ,typename DataStore<Scalar>::DataStoreVector_t *data_out
00075         ) const;
00076 
00078     int order() const; 
00079 
00081 
00082     std::string description() const;
00083 
00085     void describe(
00086       Teuchos::FancyOStream       &out
00087       ,const Teuchos::EVerbosityLevel      verbLevel
00088       ) const;
00089 
00091 
00092     void setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList);
00093 
00095     Teuchos::RCP<Teuchos::ParameterList> getParameterList();
00096 
00098     Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
00099 
00100     void assertInterpolatePreconditions(
00101         const typename DataStore<Scalar>::DataStoreVector_t &data_in
00102         ,const Array<Scalar> &t_values
00103         ,typename DataStore<Scalar>::DataStoreVector_t *data_out
00104         ) const;
00105 
00106   private:
00107 
00108     Teuchos::RCP<Teuchos::ParameterList> parameterList;
00109 
00110 };
00111 
00112 template<class Scalar>
00113 HermiteInterpolator<Scalar>::HermiteInterpolator()
00114 {
00115 }
00116 
00117 template<class Scalar>
00118 void HermiteInterpolator<Scalar>::interpolate(
00119     const typename DataStore<Scalar>::DataStoreVector_t &data_in
00120     ,const Array<Scalar> &t_values
00121     ,typename DataStore<Scalar>::DataStoreVector_t *data_out ) const
00122 {
00123   typedef Teuchos::ScalarTraits<Scalar> ST;
00124 #ifdef TEUCHOS_DEBUG
00125   assertInterpolatePreconditions(data_in,t_values,data_out);
00126 #endif // TEUCHOS_DEBUG
00127   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00128   Teuchos::OSTab ostab(out,1,"HI::interpolator");
00129   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00130     *out << "data_in:" << std::endl;
00131     for (unsigned int i=0 ; i<data_in.size() ; ++i) {
00132       *out << "data_in[" << i << "] = " << std::endl;
00133       data_in[i].describe(*out,Teuchos::VERB_EXTREME);
00134     }
00135     *out << "t_values = " << std::endl;
00136     for (unsigned int i=0 ; i<t_values.size() ; ++i) {
00137       *out << "t_values[" << i << "] = " << t_values[i] << std::endl;
00138     }
00139     for (unsigned int i=0; i<data_out->size() ; ++i) {
00140       *out << "data_out[" << i << "] = " << std::endl;
00141       (*data_out)[i].describe(*out,Teuchos::VERB_EXTREME);
00142     }
00143   }
00144   data_out->clear();
00145   if (t_values.size() == 0) {
00146     return;
00147   }
00148   
00149   if (data_in.size() == 1) {
00150     // trivial case of one node
00151     // preconditions assert that t_values[0] == data_in[0].time so we can just pass it out
00152     DataStore<Scalar> DS(data_in[0]);
00153     data_out->push_back(DS);
00154   } else {
00155     // data_in.size() >= 2
00156     int n = 0;
00157     for (int i=0 ; i<Teuchos::as<int>(data_in.size())-1 ; ++i) {
00158       const Scalar& ti = data_in[i].time;
00159       const Scalar& tip1 = data_in[i+1].time;
00160       while ((ti <= t_values[n]) && (t_values[n] <= tip1)) {
00161         const Scalar& t = t_values[n];
00162         // First we check for exact node matches:
00163         if (t == ti) {
00164           DataStore<Scalar> DS(data_in[i]);
00165           data_out->push_back(DS);
00166         } else if (t == tip1) {
00167           DataStore<Scalar> DS(data_in[i+1]);
00168           data_out->push_back(DS);
00169         } else {
00170           Teuchos::RCP<const Thyra::VectorBase<Scalar> > xi =      data_in[i  ].x;
00171           Teuchos::RCP<const Thyra::VectorBase<Scalar> > xip1 =    data_in[i+1].x;
00172           Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdoti =   data_in[i  ].xdot;
00173           Teuchos::RCP<const Thyra::VectorBase<Scalar> > xdotip1 = data_in[i+1].xdot;
00174           
00175           // 10/10/06 tscoffe:  this could be expensive:
00176           Teuchos::RCP<Thyra::VectorBase<Scalar> > tmp_vec = xi->clone_v(); 
00177           Teuchos::RCP<Thyra::VectorBase<Scalar> > xdot_temp = xip1->clone_v(); 
00178           Scalar dt = tip1-ti;
00179           Scalar dt2 = dt*dt;
00180           Scalar t_t0 = t - ti;
00181           Scalar t_t1 = t - tip1;
00182           Scalar tmp_t;
00183 
00184           // Compute numerical divided difference:
00185           Thyra::Vt_S(&*xdot_temp,Scalar(ST::one()/dt));
00186           Thyra::Vp_StV(&*xdot_temp,Scalar(-ST::one()/dt),*xi);
00187 
00188           // interpolate this point
00189           DataStore<Scalar> DS;
00190           DS.time = t;
00191 
00192           //  H_3(x) = f(x0) + f'(x0)(x-x0) + ((f(x1)-f(x0))/(x1-x0) - f'(x0))(x-x0)^2/(x1-x0)
00193           //           +(f'(x1) - 2(f(x1)-f(x0))/(x1-x0) + f'(x0))(x-x0)^2(x-x1)/(x1-x0)^2
00194           Teuchos::RCP<Thyra::VectorBase<Scalar> > x_vec = xi->clone_v(); 
00195           Thyra::Vp_StV(&*x_vec,t_t0,*xdoti);
00196           tmp_t = t_t0*t_t0/dt;
00197           Thyra::V_StVpStV(&*tmp_vec,tmp_t,*xdot_temp,Scalar(-ST::one()*tmp_t),*xdoti);
00198           Thyra::Vp_V(&*x_vec,*tmp_vec);
00199           tmp_t = t_t0*t_t0*t_t1/dt2;
00200           Thyra::V_StVpStV(&*tmp_vec,tmp_t,*xdotip1,Scalar(-2*tmp_t),*xdot_temp);
00201           Thyra::Vp_StV(&*tmp_vec,tmp_t,*xdoti);
00202           Thyra::Vp_V(&*x_vec,*tmp_vec);
00203           DS.x = x_vec;
00204 
00205           //  H_3'(x) =        f'(x0) + 2*((f(x1)-f(x0))/(x1-x0) - f'(x0))(x-x0)/(x1-x0)
00206           //           +(f'(x1) - 2(f(x1)-f(x0))/(x1-x0) + f'(x0))[2*(x-x0)(x-x1) + (x-x0)^2]/(x1-x0)^2
00207           Teuchos::RCP<Thyra::VectorBase<Scalar> > xdot_vec = xdoti->clone_v(); 
00208           tmp_t = t_t0/dt;
00209           Thyra::Vp_StV(&*xdot_vec,Scalar(2*tmp_t),*xdot_temp);
00210           Thyra::Vp_StV(&*xdot_vec,Scalar(-ST::one()*tmp_t),*xdoti);
00211           tmp_t = Scalar((2*t_t0*t_t1+t_t0*t_t0)/dt2);
00212           Thyra::V_StVpStV(&*tmp_vec,tmp_t,*xdotip1,Scalar(-2*tmp_t),*xdot_temp);
00213           Thyra::Vp_StV(&*tmp_vec,tmp_t,*xdoti);
00214           Thyra::Vp_V(&*xdot_vec,*tmp_vec);
00215           DS.xdot = xdot_vec;
00216           
00217           // Accuracy:
00218           // f(x) - H_3(x) = (f^{(3)}(\xi(x))/(4!))(x-x0)^2(x-x1)^2
00219           DS.accuracy = (t_t0)*(t_t0)*(t_t1)*(t_t1);
00220 
00221           // Push DataStore object onto vector:
00222           data_out->push_back(DS);
00223         }
00224         n++;
00225         if (n == Teuchos::as<int>(t_values.size())) {
00226           return;
00227         }
00228       }
00229     }
00230   } // data_in.size() == 1
00231 }
00232 
00233 template<class Scalar>
00234 int HermiteInterpolator<Scalar>::order() const
00235 {
00236   return(2);
00237 }
00238 
00239 template<class Scalar>
00240 std::string HermiteInterpolator<Scalar>::description() const
00241 {
00242   std::string name = "Rythmos::HermiteInterpolator";
00243   return(name);
00244 }
00245 
00246 template<class Scalar>
00247 void HermiteInterpolator<Scalar>::describe(
00248       Teuchos::FancyOStream       &out
00249       ,const Teuchos::EVerbosityLevel      verbLevel
00250       ) const
00251 {
00252   if ( (Teuchos::as<int>(verbLevel) == Teuchos::as<int>(Teuchos::VERB_DEFAULT) ) ||
00253        (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW)     )
00254      ) {
00255     out << description() << "::describe" << std::endl;
00256   } else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW)) {
00257   } else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_MEDIUM)) {
00258   } else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH)) {
00259   }
00260 }
00261 
00262 template <class Scalar>
00263 void HermiteInterpolator<Scalar>::setParameterList(Teuchos::RCP<Teuchos::ParameterList> const& paramList)
00264 {
00265   parameterList = paramList;
00266   int outputLevel = parameterList->get( "outputLevel", int(-1) );
00267   outputLevel = std::min(std::max(outputLevel,-1),4);
00268   this->setVerbLevel(static_cast<Teuchos::EVerbosityLevel>(outputLevel));
00269 }
00270 
00271 template <class Scalar>
00272 Teuchos::RCP<Teuchos::ParameterList> HermiteInterpolator<Scalar>::getParameterList()
00273 {
00274   return(parameterList);
00275 }
00276 
00277 template <class Scalar>
00278 Teuchos::RCP<Teuchos::ParameterList> HermiteInterpolator<Scalar>::unsetParameterList()
00279 {
00280   Teuchos::RCP<Teuchos::ParameterList> temp_param_list = parameterList;
00281   parameterList = Teuchos::null;
00282   return(temp_param_list);
00283 }
00284 
00285 template<class Scalar>
00286 void HermiteInterpolator<Scalar>::assertInterpolatePreconditions(
00287         const typename DataStore<Scalar>::DataStoreVector_t &data_in
00288         ,const Array<Scalar> &t_values
00289         ,typename DataStore<Scalar>::DataStoreVector_t *data_out
00290         ) const
00291 {
00292   assertBaseInterpolatePreconditions(data_in,t_values,data_out);
00293   for (int i=0; i<Teuchos::as<int>(data_in.size()) ; ++i) {
00294     TEST_FOR_EXCEPTION(
00295         data_in[i].xdot == Teuchos::null, std::logic_error,
00296         "Error, data_in[" << i << "].xdot == Teuchos::null.\n"
00297         );
00298   }
00299 }
00300 
00301 } // namespace Rythmos
00302 
00303 #endif // Rythmos_HERMITE_INTERPOLATOR_H

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