Rythmos_HermiteInterpolator_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_HERMITE_INTERPOLATOR_DEF_H
00030 #define Rythmos_HERMITE_INTERPOLATOR_DEF_H
00031 
00032 #include "Rythmos_HermiteInterpolator_decl.hpp"
00033 #include "Rythmos_InterpolatorBaseHelpers.hpp"
00034 #include "Thyra_VectorStdOps.hpp"
00035 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00036 
00037 namespace Rythmos {
00038 
00039 template<class Scalar>
00040 HermiteInterpolator<Scalar>::HermiteInterpolator()
00041 {
00042   nodes_ = Teuchos::null;
00043   parameterList_ = Teuchos::null;
00044 }
00045 
00046 template<class Scalar>
00047 void HermiteInterpolator<Scalar>::setNodes(
00048   const RCP<const typename DataStore<Scalar>::DataStoreVector_t> & nodes
00049   )
00050 {
00051   nodes_ = nodes;
00052 }
00053 
00054 template<class Scalar>
00055 void HermiteInterpolator<Scalar>::interpolate(
00056     const Array<Scalar> &t_values
00057     ,typename DataStore<Scalar>::DataStoreVector_t *data_out ) const
00058 {
00059 
00060   //TEST_FOR_EXCEPT_MSG(true, "Error, ths function is not tested!" );
00061 
00062   typedef Teuchos::ScalarTraits<Scalar> ST;
00063 #ifdef RYTHMOS_DEBUG
00064   assertInterpolatePreconditions((*nodes_),t_values,data_out);
00065 #endif // RYTHMOS_DEBUG
00066   RCP<Teuchos::FancyOStream> out = this->getOStream();
00067   Teuchos::OSTab ostab(out,1,"HI::interpolator");
00068   if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
00069     *out << "(*nodes_):" << std::endl;
00070     for (unsigned int i=0 ; i<(*nodes_).size() ; ++i) {
00071       *out << "(*nodes_)[" << i << "] = " << std::endl;
00072       (*nodes_)[i].describe(*out,Teuchos::VERB_EXTREME);
00073     }
00074     *out << "t_values = " << std::endl;
00075     for (unsigned int i=0 ; i<t_values.size() ; ++i) {
00076       *out << "t_values[" << i << "] = " << t_values[i] << std::endl;
00077     }
00078     for (unsigned int i=0; i<data_out->size() ; ++i) {
00079       *out << "data_out[" << i << "] = " << std::endl;
00080       (*data_out)[i].describe(*out,Teuchos::VERB_EXTREME);
00081     }
00082   }
00083   data_out->clear();
00084   if (t_values.size() == 0) {
00085     return;
00086   }
00087   
00088   if ((*nodes_).size() == 1) {
00089     // trivial case of one node
00090     // preconditions assert that t_values[0] == (*nodes_)[0].time so we can just pass it out
00091     DataStore<Scalar> DS((*nodes_)[0]);
00092     data_out->push_back(DS);
00093   } else {
00094     // (*nodes_).size() >= 2
00095     int n = 0;
00096     for (int i=0 ; i<Teuchos::as<int>((*nodes_).size())-1 ; ++i) {
00097       const Scalar& t0 = (*nodes_)[i].time;
00098       const Scalar& t1 = (*nodes_)[i+1].time;
00099       while ((t0 <= t_values[n]) && (t_values[n] <= t1)) {
00100         const Scalar& t = t_values[n];
00101         // First we check for exact node matches:
00102         if (t == t0) {
00103           DataStore<Scalar> DS((*nodes_)[i]);
00104           data_out->push_back(DS);
00105         } else if (t == t1) {
00106           DataStore<Scalar> DS((*nodes_)[i+1]);
00107           data_out->push_back(DS);
00108         } else {
00109           RCP<const Thyra::VectorBase<Scalar> > x0    = (*nodes_)[i  ].x;
00110           RCP<const Thyra::VectorBase<Scalar> > x1    = (*nodes_)[i+1].x;
00111           RCP<const Thyra::VectorBase<Scalar> > xdot0 = (*nodes_)[i  ].xdot;
00112           RCP<const Thyra::VectorBase<Scalar> > xdot1 = (*nodes_)[i+1].xdot;
00113           
00114           // 10/10/06 tscoffe:  this could be expensive:
00115           RCP<Thyra::VectorBase<Scalar> > tmp_vec = x0->clone_v(); 
00116           RCP<Thyra::VectorBase<Scalar> > xdot_temp = x1->clone_v(); 
00117           Scalar dt = t1-t0;
00118           Scalar dt2 = dt*dt;
00119           Scalar t_t0 = t - t0;
00120           Scalar t_t1 = t - t1;
00121           Scalar tmp_t;
00122 
00123           // Compute numerical divided difference:
00124           Thyra::Vt_S(&*xdot_temp,Scalar(ST::one()/dt));
00125           Thyra::Vp_StV(&*xdot_temp,Scalar(-ST::one()/dt),*x0);
00126 
00127           // interpolate this point
00128           DataStore<Scalar> DS;
00129           DS.time = t;
00130 
00131           //  H_3(t) = x(t0) + xdot(t0)(t-t0) + ((x(t1)-x(t0))/(t1-t0) - xdot(t0))(t-t0)^2/(t1-t0)
00132           //           +(xdot(t1) - 2(x(t1)-x(t0))/(t1-t0) + xdot(t0))(t-t0)^2(t-t1)/(t1-t0)^2
00133           RCP<Thyra::VectorBase<Scalar> > x_vec = x0->clone_v(); 
00134           Thyra::Vp_StV(&*x_vec,t_t0,*xdot0);
00135           tmp_t = t_t0*t_t0/dt;
00136           Thyra::V_StVpStV(&*tmp_vec,tmp_t,*xdot_temp,Scalar(-ST::one()*tmp_t),*xdot0);
00137           Thyra::Vp_V(&*x_vec,*tmp_vec);
00138           tmp_t = t_t0*t_t0*t_t1/dt2;
00139           Thyra::V_StVpStV(&*tmp_vec,tmp_t,*xdot1,Scalar(-2*tmp_t),*xdot_temp);
00140           Thyra::Vp_StV(&*tmp_vec,tmp_t,*xdot0);
00141           Thyra::Vp_V(&*x_vec,*tmp_vec);
00142           DS.x = x_vec;
00143 
00144           //  H_3'(t) =        xdot(t0) + 2*((x(t1)-x(t0))/(t1-t0) - xdot(t0))(t-t0)/(t1-t0)
00145           //           +(xdot(t1) - 2(x(t1)-x(t0))/(t1-t0) + xdot(t0))[2*(t-t0)(t-t1) + (t-t0)^2]/(t1-t0)^2
00146           RCP<Thyra::VectorBase<Scalar> > xdot_vec = xdot0->clone_v(); 
00147           tmp_t = t_t0/dt;
00148           Thyra::Vp_StV(&*xdot_vec,Scalar(2*tmp_t),*xdot_temp);
00149           Thyra::Vp_StV(&*xdot_vec,Scalar(-2*tmp_t),*xdot0);
00150           tmp_t = Scalar((2*t_t0*t_t1+t_t0*t_t0)/dt2);
00151           Thyra::V_StVpStV(&*tmp_vec,tmp_t,*xdot1,Scalar(-2*tmp_t),*xdot_temp);
00152           Thyra::Vp_StV(&*tmp_vec,tmp_t,*xdot0);
00153           Thyra::Vp_V(&*xdot_vec,*tmp_vec);
00154           DS.xdot = xdot_vec;
00155           
00156           // Accuracy:
00157           // f(x) - H_3(x) = (f^{(3)}(\xi(x))/(4!))(x-x0)^2(x-x1)^2
00158           DS.accuracy = (t_t0)*(t_t0)*(t_t1)*(t_t1);
00159 
00160           // Push DataStore object onto vector:
00161           data_out->push_back(DS);
00162         }
00163         n++;
00164         if (n == Teuchos::as<int>(t_values.size())) {
00165           return;
00166         }
00167       }
00168     }
00169   } // (*nodes_).size() == 1
00170 }
00171 
00172 // non-member constructor
00173 template<class Scalar>
00174 RCP<HermiteInterpolator<Scalar> > hermiteInterpolator()
00175 {
00176   RCP<HermiteInterpolator<Scalar> > hi = rcp(new HermiteInterpolator<Scalar>() );
00177   return hi;
00178 }
00179 
00180 template<class Scalar>
00181 int HermiteInterpolator<Scalar>::order() const
00182 {
00183   return(3);
00184 }
00185 
00186 template<class Scalar>
00187 std::string HermiteInterpolator<Scalar>::description() const
00188 {
00189   std::string name = "Rythmos::HermiteInterpolator";
00190   return(name);
00191 }
00192 
00193 template<class Scalar>
00194 void HermiteInterpolator<Scalar>::describe(
00195       Teuchos::FancyOStream       &out
00196       ,const Teuchos::EVerbosityLevel      verbLevel
00197       ) const
00198 {
00199   if ( (Teuchos::as<int>(verbLevel) == Teuchos::as<int>(Teuchos::VERB_DEFAULT) ) ||
00200        (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW)     )
00201      )
00202   {
00203     out << description() << "::describe" << std::endl;
00204   }
00205   else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW))
00206   {}
00207   else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_MEDIUM))
00208   {}
00209   else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH))
00210   {}
00211 }
00212 
00213 template <class Scalar>
00214 void HermiteInterpolator<Scalar>::setParameterList(RCP<ParameterList> const& paramList)
00215 {
00216   TEST_FOR_EXCEPT(is_null(paramList));
00217   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00218   parameterList_ = paramList;
00219   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00220 }
00221 
00222 template <class Scalar>
00223 RCP<ParameterList> HermiteInterpolator<Scalar>::getNonconstParameterList()
00224 {
00225   return(parameterList_);
00226 }
00227 
00228 template <class Scalar>
00229 RCP<ParameterList> HermiteInterpolator<Scalar>::unsetParameterList()
00230 {
00231   RCP<ParameterList> temp_param_list = parameterList_;
00232   parameterList_ = Teuchos::null;
00233   return(temp_param_list);
00234 }
00235 
00236 template<class Scalar>
00237 RCP<const Teuchos::ParameterList> HermiteInterpolator<Scalar>::getValidParameters() const
00238 {
00239   static RCP<Teuchos::ParameterList> validPL;
00240   if (is_null(validPL)) {
00241     RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00242     Teuchos::setupVerboseObjectSublist(&*pl);
00243     validPL = pl;
00244   }
00245   return (validPL);
00246 }
00247 
00248 template<class Scalar>
00249 void HermiteInterpolator<Scalar>::assertInterpolatePreconditions(
00250         const typename DataStore<Scalar>::DataStoreVector_t &data_in
00251         ,const Array<Scalar> &t_values
00252         ,typename DataStore<Scalar>::DataStoreVector_t *data_out
00253         ) const
00254 {
00255   assertBaseInterpolatePreconditions(data_in,t_values,data_out);
00256   for (int i=0; i<Teuchos::as<int>(data_in.size()) ; ++i) {
00257     TEST_FOR_EXCEPTION(
00258         is_null(data_in[i].x), std::logic_error,
00259         "Error, data_in[" << i << "].x == Teuchos::null.\n"
00260         );
00261     TEST_FOR_EXCEPTION(
00262         is_null(data_in[i].xdot), std::logic_error,
00263         "Error, data_in[" << i << "].xdot == Teuchos::null.\n"
00264         );
00265   }
00266 }
00267 
00268 // 
00269 // Explicit Instantiation macro
00270 //
00271 // Must be expanded from within the Rythmos namespace!
00272 //
00273 
00274 #define RYTHMOS_HERMITE_INTERPOLATOR_INSTANT(SCALAR) \
00275   \
00276   template class HermiteInterpolator< SCALAR >; \
00277   \
00278   template RCP<HermiteInterpolator< SCALAR > > hermiteInterpolator(); 
00279 
00280 
00281 } // namespace Rythmos
00282 
00283 #endif // Rythmos_HERMITE_INTERPOLATOR_DEF_H

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