Rythmos_LinearInterpolationBuffer.hpp

Go to the documentation of this file.
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_INTERPOLATION_BUFFER_H
00030 #define Rythmos_LINEAR_INTERPOLATION_BUFFER_H
00031 
00032 #include "Rythmos_InterpolationBuffer.hpp"
00033 #include "Thyra_VectorBase.hpp"
00034 
00035 namespace Rythmos {
00036 
00038 template<class Scalar> 
00039 class LinearInterpolationBuffer : virtual public InterpolationBuffer<Scalar>
00040 {
00041   public:
00042 
00043     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00044     
00046     LinearInterpolationBuffer();
00047     LinearInterpolationBuffer( int storage );
00048 
00049     SetStorage( int storage );
00050         
00052     ~LinearInterpolationBuffer() {};
00053 
00055     bool SetPoints(
00056       const std::vector<Scalar>& time_list
00057       ,const std::vector<Thyra::VectorBase<Scalar> >& x_list
00058       ,const std::vector<Thyra::VectorBase<Scalar> >& xdot_list);
00059 
00061     bool GetPoints(
00062       const std::vector<Scalar>& time_list
00063       ,std::vector<Thyra::VectorBase<Scalar> >* x_list
00064       ,std::vector<Thyra::VectorBase<Scalar> >* xdot_list
00065       ,std::vector<ScalarMag>* accuracy_list) const;
00066 
00068     bool SetRange(
00069       const Scalar& time_lower
00070       ,const Scalar& time_upper
00071       ,const InterpolationBuffer<Scalar>& IB);
00072 
00074     bool GetNodes(std::vector<Scalar>* time_list) const;
00075 
00077     bool RemoveNodes(std::vector<Scalar>* time_list) const;
00078 
00080     int GetOrder() const;
00081 
00082   private:
00083 
00084     int storage_limit;
00085     Thyra::VectorBase<Scalar> tmp_vec;
00086     std::list<Teuchos::RefCountPtr<DataStore<Scalar> > > node_list;
00087     int order;
00088 
00089     template<class Scalar>
00090     class DataStore
00091     {
00092       public:
00093         ~DataStore();
00094         DataStore();
00095         DataStore(ScalarMag &time
00096           ,Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > &x
00097           ,Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > &xdot);
00098         ScalarMag time;
00099         Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > x;
00100         Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > xdot;
00101         bool operator< (const DataStore<Scalar>& d1, const DataStore<Scalar>& d2) const
00102         { return( d1.time < d2.time ); }
00103     }
00104     void VectorToDataStoreList(
00105       std::list<Teuchos::RefCountPtr<DataStore<Scalar> > > *list_ds
00106       ,const std::vector<ScalarMag> &time_vec
00107       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > &x_vec
00108       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > &xdot_vec) const;
00109     void DataStoreListToVector(
00110       std::vector<ScalarMag> *time_vec
00111       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > *x_vec
00112       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > *xdot_vec
00113       ,const std::list<Teuchos::RefCountPtr<DataStore<Scalar> > > &list_ds) const;
00114 
00115 
00116 };
00117 
00118 // ////////////////////////////
00119 // Defintions
00120 template<class Scalar>
00121 LinearInterpolationBuffer<Scalar>::LinearInterpolationBuffer()
00122 {
00123   order = 1;
00124   SetStorage(2);
00125 }
00126 
00127 template<class Scalar>
00128 LinearInterpolationBuffer<Scalar>::LinearInterpolationBuffer( int storage_ )
00129 {
00130   order = 1;
00131   SetStorage(storage_);
00132 }
00133 
00134 template<class Scalar>
00135 LinearInterpolationBuffer<Scalar>::SetStorage( int storage_ )
00136 {
00137   storage_limit = min(2,storage_); // Minimum of two points so interpolation is possible
00138 }
00139 
00140 template<class Scalar>
00141 bool LinearInterpolationBuffer<Scalar>::SetPoints( 
00142     const std::vector<Scalar>& time_vec
00143     ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& x_vec
00144     ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >& xdot_vec );
00145 {
00146   std::list<DataStore<Scalar> > input_list;
00147   VectorToDataStoreList(&input_list,time_vec,x_vec,xdot_vec);
00148   input_list.sort();
00149   // Determine if time is already in list and if so, replace existing data with new data.
00150   std::list<DataStore<Scalar> >::iterator node_it;
00151   std::list<DataStore<Scalar> >::iterator input_it = input_list.begin();
00152   for (; input_it != input_list.end() ; input_it++)
00153   {
00154     node_it = find(node_list.begin(),node_list.end(),*input_it);
00155     if (node_it != node_list.end())
00156     {
00157       node_list.splice(node_it,input_list,input_it);
00158       node_list.erase(node_it);
00159     }
00160   }
00161   // Check that we're not going to exceed our storage limit:
00162   if ((node_list.size()+input_list.size()) > storage_limit)
00163     return(false);
00164   // Now add all the remaining points to node_list
00165   node_list.merge(input_list);
00166   return(true);
00167 }
00168 
00169 template<class Scalar>
00170 bool LinearInterpolationBuffer<Scalar>::GetPoints(
00171     const std::vector<Scalar>& time_vec
00172     ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > >* x_vec
00173     ,std::vector<Tuechos::RefCountPtr<Thyra::VectorBase<Scalar> > >* xdot_vec
00174     ,std::vector<ScalarMag>* accuracy_vec) const
00175 {
00176   // Copy the const time_vec to a local sorted time_vec
00177   std::vector<Scalar> local_time_vec = time_vec;
00178   std::sort(local_time_vec.begin(),local_time_vec.end());
00179   // If there are fewer than 2 points in node_list, then return failure
00180   if (node_list.size() < 2)
00181     return(false);
00182   // If time is outside range of t_values, then return failure
00183   if ( (*local_time_vec.begin() < node_list.begin().t) || (*local_time_vec.end() > node_list.end().t) )
00184     return(false);
00185   // Find t values on either side of time
00186   std::vector<Scalar>::iterator input_it = local_time_vec.begin();
00187   std::vector<DataStore<Scalar> >::iterator node_it = node_list.begin();
00188   for ( ; node_it != node_list.end() ; node_it++ )
00189   {
00190     while ((*input_it >= node_it->t) && (*input_it <= (node_it+1)->t))
00191     {
00192       Scalar& t = *input_it;
00193       Scalar& ti = node_it->t;
00194       Scalar& tip1 = (node_it+1)->t;
00195       Thyra::VectorBase<Scalar>& xi = *(node_it->x);
00196       Thyra::VectorBase<Scalar>& xip1 = *((node_it+1)->x);
00197       Thyra::VectorBase<Scalar>& xdoti = *(node_it->xdot);
00198       Thyra::VectorBase<Scalar>& xdotip1 = *((node_it+1)->xdot);
00199 
00200       // interpolate this point
00201       Scalar h = tip1-ti;
00202       // First we work on x.
00203       V_StVpStV(&*tmp_vec,Scalar(ST::one()/h),xi,Scalar(-ST::one()/h),xip1);
00204       Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > x = (node_it->x).clone_v();
00205       V_StVpStV(&*x, ST::one(), xi, t-ti, tmp_vec);
00206       x_vec.pushback(x);
00207       // Then we work on xdot.
00208       V_StVpStV(&*tmp_vec,Scalar(ST::one()/h),xdoti,Scalar(-ST::one()/h),xdotip1);
00209       Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > xdot = (node_it->xdot).clone_v();
00210       V_StVpStV(&*xdot, ST::one(), xdoti, t-ti, tmp_vec);
00211       xdot_vec.pushback(xdot);
00212       // And finally we estimate our order of accuracy
00213       accuracy_vec.pushback(h); 
00214       // Now we increment iterator for local_time_vec
00215       input_it++;
00216     }
00217   }
00218   return(true);
00219 }
00220 
00221 template<class Scalar>
00222 bool LinearInterpolationBuffer<Scalar>::SetRange(
00223     const Scalar& time_lower
00224     ,const Scalar& time_upper
00225     ,const InterpolationBuffer<Scalar>& IB )
00226 {
00227   std::vector<ScalarMag> input_nodes;
00228   bool status = IB.GetNodes(&input_nodes);
00229   if (status == false) return(status);
00230   std::sort(input_nodes.begin(),input_nodes.end());
00231   // Remove nodes outside the range [time_lower,time_upper]
00232   std::vector<ScalarMag>::iterator input_it = input_nodes.begin();
00233   for (; input_it != input_nodes.end() ; input_it++)
00234   {
00235     if (*input_it >= time_lower)
00236     {
00237       input_it--;
00238       break;
00239     }
00240   }
00241   input_nodes.erase(input_nodes.begin(),input_it);
00242   input_it = input_nodes.end();
00243   for (; input_it != input_nodes.begin() ; input_it--)
00244   {
00245     if (*input_it <= time_upper)
00246     {
00247       input_it++;
00248       break;
00249     }
00250   }
00251   input_nodes.erase(input_it,input_nodes.end());
00252 
00253   // Ask IB to interpolate more points if IB's order is higher than ours
00254   ScalarMag h_safety = ScalarMag(2*ST::one());
00255   int IBOrder = IB.GetOrder();
00256   if (IBOrder >= order)
00257   {
00258     input_it = input_nodes.begin();
00259     for (; input_it != input_nodes.end() ; input_it++)
00260     {
00261       std::vector<ScalarMag>::iterator input_it_next = input_it++;
00262       if (input_it_next == input_nodes.end())
00263         break;
00264       ScalarMag h_0 = *input_it_next - *input_it;
00265       ScalarMag h = h_0^(IBOrder/order)/h_safety;
00266       int N = ceil(h_0/h);
00267       h = ScalarMag(h_0/N);
00268       for (int i=1 ; i<N ; ++i)
00269       {
00270         input_nodes.insert(input_it_next,*input_it+i*h);
00271       }
00272     }
00273   }
00274   // If IB's order is lower than ours, then simply grab the node values and continue.
00275   // If IB's order is higher than ours, then grab the node values and ask IB to
00276   // interpolate extra values so that our order of accuracy is approximately
00277   // the same as the other IB's.
00278   // One approach:
00279   // Lets say IB's order is p and our order is r (p>r).
00280   // Given a particular interval with spacing h_0, the order of accuracy of IB is h_0^p
00281   // We want to find a spacing h such that h^r = h_0^p.  Clearly, this is h = h_0^{p/r}.
00282   // Given this new spacing, divide up the interval h_0 into h_0/h subintervals
00283   // and ask for the IB to interpolate points.  This will match basic order of
00284   // accuracy estimates.  Its probably a good idea to include a fudge factor in
00285   // there too.  E.g. h = h_0^{p/r}/fudge.
00286 
00287   // Don't forget to check the interval [time_lower,time_upper].
00288   // Use SetPoints and check return value to make sure we observe storage_limit.
00289 
00290   std::vector<Teuchos::RefCountPtr<Scalar> > input_x;
00291   std::vector<Teuchos::RefCountPtr<Scalar> > input_xdot;
00292   std::vector<ScalarMag> input_accuracy;
00293   status = IB.GetPoints( input_nodes, &input_x, &input_xdot &input_accuracy );
00294   if (status == false) return(status);
00295   // We could check that the accuracy meets our criteria here.
00296   status = SetPoints( input_nodes, input_x, input_xdot );
00297   return(status);
00298 }
00299 
00300 template<class Scalar>
00301 bool LinearInterpolationBuffer<Scalar>::GetNodes( std::vector<Scalar>* time_list ) const
00302 {
00303   std::list<DataStore<Scalar> >::iterator list_it = node_list.begin();
00304   for (; list_it != node_list.end() ; list_it++)
00305   {
00306     time_list->push_back(list_it->t);
00307   }
00308   return(true);
00309 }
00310 
00311 template<class Scalar>
00312 bool LinearInterpolationBuffer<Scalar>::RemoveNodes( std::vector<Scalar>& time_vec ) const
00313 {
00314   int N = time_vec.size();
00315   for (int i=0; i<N ; ++i)
00316   {
00317     node_list.remove(time_vec[i]);
00318   }
00319   return(true);
00320 }
00321 
00322 template<class Scalar>
00323 int LinearInterpolationBuffer<Scalar>::GetOrder() const
00324 {
00325   return(order);
00326 }
00327 
00328 // DataStore definitions:
00329 template<class Scalar>
00330 DataStore<Scalar>::DataStore(ScalarMag &time_
00331   ,Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > &x_
00332   ,Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > &xdot_)
00333 {
00334   time = time_;
00335   x = x_;
00336   xdot = xdot_;
00337 }
00338 
00339 template<class Scalar>
00340 void LinearInterpolationBuffer<Scalar>::VectorToDataStoreList(
00341       std::list<Teuchos::RefCountPtr<DataStore<Scalar> > > *list_ds
00342       ,const std::vector<ScalarMag> &time_vec
00343       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > &x_vec
00344       ,const std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > &xdot_vec) const
00345 {
00346   int N = time_vec.size();
00347   int Nx = x_vec.size();
00348   int Nxdot = xdot_vec.size();
00349   if ((N != Nx) || (N != Nxdot))
00350   {
00351     list_ds = NULL;
00352     return;
00353   }
00354   list_ds->clear();
00355   for (int i=0; i<N ; ++i)
00356   {
00357     list_ds->push_back(Teuchos::rcp(new DataStore(time_list[i],x_list[i],xdot_list[i])));
00358   }
00359 }
00360 
00361 template<class Scalar>
00362 void LinearInterpolationBuffer<Scalar>::DataStoreListToVector(
00363       std::vector<ScalarMag> *time_vec
00364       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > *x_vec
00365       ,std::vector<Teuchos::RefCountPtr<Thyra::VectorBase<Scalar> > > *xdot_vec
00366       ,const std::list<Teuchos::RefCountPtr<DataStore<Scalar> > > &list_ds) const
00367 {
00368   int N = list_ds.size();
00369   time_list->reserve(N); time_list->clear();
00370   x_list->reserve(N);    x_list->clear();
00371   xdot_list->reserve(N); xdot_list->clear();
00372   std::list<DataStore<Scalar> >::iterator list_it = list_ds.begin();
00373   for (; list_it != list_ds.end() ; list_it++)
00374   {
00375     time_list->push_back(list_it->time);
00376     x_list->push_back(list_it->x);
00377     x_dot_list->push_back(list_it->xdot);
00378   }
00379 }
00380 
00381 
00382 } // namespace Rythmos
00383 
00384 #endif // Rythmos_LINEAR_INTERPOLATION_BUFFER_H

Generated on Thu Sep 18 12:30:05 2008 for Rythmos - Transient Integration for Differential Equations by doxygen 1.3.9.1