Rythmos_InterpolationBufferAppender.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_INTERPOLATION_BUFFER_APPENDER_H
00030 #define Rythmos_INTERPOLATION_BUFFER_APPENDER_H
00031 
00032 #include "Rythmos_InterpolationBufferBase.hpp"
00033 #include "Rythmos_Types.hpp"
00034 
00035 #include "Thyra_VectorBase.hpp"
00036 
00037 #include "Teuchos_Describable.hpp"
00038 #include "Teuchos_ParameterListAcceptor.hpp"
00039 #include "Teuchos_VerboseObject.hpp"
00040 #include "Teuchos_RCP.hpp"
00041 #include "Teuchos_implicit_cast.hpp"
00042 #include "Teuchos_Assert.hpp"
00043 
00044 
00045 namespace Rythmos {
00046 
00047 template<class Scalar>
00048 class InterpolationBufferAppenderBase
00049   : virtual public Teuchos::Describable
00050   , virtual public Teuchos::ParameterListAcceptor
00051   , virtual public Teuchos::VerboseObject<InterpolationBufferAppenderBase<Scalar> >
00052 {
00053   public:
00054 
00082     virtual void import(
00083         InterpolationBufferBase<Scalar>* IB_base, 
00084         const InterpolationBufferBase<Scalar>& IB_in, 
00085         const TimeRange<Scalar>& range
00086         ) =0;
00087   protected:
00088     void assertImportPreconditions(
00089         const InterpolationBufferBase<Scalar>& IB_base, 
00090         const InterpolationBufferBase<Scalar>& IB_in, 
00091         const TimeRange<Scalar>& range
00092         ) const;
00093 };
00094 
00095 template<class Scalar>
00096 void InterpolationBufferAppenderBase<Scalar>::assertImportPreconditions(
00097         const InterpolationBufferBase<Scalar>& IB_base, 
00098         const InterpolationBufferBase<Scalar>& IB_in, 
00099         const TimeRange<Scalar>& range
00100         ) const
00101 {
00102   TEST_FOR_EXCEPTION((range.lower() != IB_base.getTimeRange().upper()) || (range.upper() != IB_base.getTimeRange().lower()),
00103       std::logic_error, 
00104       "Error, import range is not an append nor a prepend of the base interpolation buffer.\n"
00105       );
00106   TEST_FOR_EXCEPTION(range.lower() < IB_in.getTimeRange().lower(),
00107       std::logic_error,
00108       "Error, import range's lower bound does not sit inside incoming interpolation buffer's time range.\n"
00109       );
00110   TEST_FOR_EXCEPTION(IB_in.getTimeRange().upper() < range.upper(),
00111       std::logic_error,
00112       "Error, import range's upper bound does not sit inside incoming interpolation buffer's time range.\n"
00113       );
00114 }
00115 
00116 
00117 template<class Scalar>
00118 class InterpolationBufferAppenderDefault : virtual public InterpolationBufferAppenderBase<Scalar>
00119 {
00120   public:
00121 
00122     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00123 
00127     void import(
00128         InterpolationBufferBase<Scalar>* IB_base, 
00129         const InterpolationBufferBase<Scalar>& IB_in, 
00130         const TimeRange<Scalar>& range
00131         );
00132     
00134 
00135     std::string description() const;
00136 
00138     void describe(
00139       Teuchos::FancyOStream       &out
00140       ,const Teuchos::EVerbosityLevel      verbLevel
00141       ) const;
00142 
00144 
00145     void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00146 
00148     RCP<Teuchos::ParameterList> getParameterList();
00149 
00151     RCP<Teuchos::ParameterList> unsetParameterList();
00152 
00153   private:
00154     RCP<Teuchos::ParameterList> parameterList_;
00155 
00156 };
00157 
00158 template<class Scalar>
00159 void InterpolationBufferAppenderDefault<Scalar>::import(
00160         InterpolationBufferBase<Scalar>* IB_base, 
00161         const InterpolationBufferBase<Scalar>& IB_in, 
00162         const TimeRange<Scalar>& range
00163     ) 
00164 {
00165 #ifdef TEUCHOS_DEBUG
00166   InterpolationBufferAppenderBase<Scalar>::assertImportPreconditions(*IB_base,IB_in,range);
00167 #endif // TEUCHOS_DEBUG
00168 
00169   Array<Scalar> time_vec;
00170   Array<RCP<const Thyra::VectorBase<Scalar> > > x_vec;
00171   Array<RCP<const Thyra::VectorBase<Scalar> > > xdot_vec;
00172   Array<ScalarMag> accuracy_vec;
00173 
00174   Array<Scalar> time_vec_in;
00175   IB_in.getNodes(&time_vec_in);
00176 
00177   selectPointsInTimeRange(&time_vec,time_vec_in,range);
00178 
00179   IB_in.getPoints(time_vec, &x_vec, &xdot_vec, &accuracy_vec);
00180   IB_base->addPoints(time_vec, x_vec, xdot_vec);
00181 }
00182 
00183 
00184 template<class Scalar>
00185 std::string InterpolationBufferAppenderDefault<Scalar>::description() const
00186 {
00187   std::string name = "Rythmos::InterpolationBufferAppenderDefault";
00188   return(name);
00189 }
00190 
00191 template<class Scalar>
00192 void InterpolationBufferAppenderDefault<Scalar>::describe(
00193       Teuchos::FancyOStream       &out
00194       ,const Teuchos::EVerbosityLevel      verbLevel
00195       ) const
00196 {
00197   if ( (Teuchos::as<int>(verbLevel) == Teuchos::as<int>(Teuchos::VERB_DEFAULT) ) ||
00198        (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW)     )
00199      ) {
00200     out << description() << "::describe" << std::endl;
00201   }
00202 }
00203 
00204 template<class Scalar>
00205 void InterpolationBufferAppenderDefault<Scalar>::setParameterList(RCP<Teuchos::ParameterList> const& paramList)
00206 {
00207   TEST_FOR_EXCEPTION(paramList==Teuchos::null,std::logic_error,"Error, paramList == Teuchos::null!\n");
00208   parameterList_ = paramList;
00209 }
00210 
00211 template<class Scalar>
00212 RCP<Teuchos::ParameterList> InterpolationBufferAppenderDefault<Scalar>::getParameterList()
00213 {
00214   return(parameterList_);
00215 }
00216 
00217 template<class Scalar>
00218 RCP<Teuchos::ParameterList> InterpolationBufferAppenderDefault<Scalar>::unsetParameterList()
00219 {
00220   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00221   parameterList_ = Teuchos::null;
00222   return(temp_param_list);
00223 }
00224 
00225 template<class Scalar>
00226 class InterpolationBufferAppenderSmart : virtual public InterpolationBufferAppenderBase<Scalar>
00227 {
00228   public:
00233     void import(
00234         InterpolationBufferBase<Scalar>* IB_base, 
00235         const InterpolationBufferBase<Scalar>& IB_in, 
00236         const TimeRange<Scalar>& range
00237         );
00238 };
00239 
00240 template<class Scalar>
00241 void InterpolationBufferAppenderSmart<Scalar>::import(
00242         InterpolationBufferBase<Scalar>* IB_base, 
00243         const InterpolationBufferBase<Scalar>& IB_in, 
00244         const TimeRange<Scalar>& range
00245     ) 
00246 {
00247 #ifdef TEUCHOS_DEBUG
00248   InterpolationBufferAppenderBase<Scalar>::assertImportPreconditions(*IB_base,IB_in,range);
00249 #endif // TEUCHOS_DEBUG
00250   if (IB_base->getOrder() >= IB_in.getOrder()) {
00251     // The incoming interpolation buffer's order of interpolation is lower than
00252     // the base interpolation buffer's order of interpolation.  In this case,
00253     // we just copy the data over.
00254     InterpolationBufferAppenderDefault<Scalar> defaultAppender;
00255     defaultAppender.import(IB_base,IB_in,range);
00256   } else {
00257     // In this case, the incoming interpolation buffer's order of interpolation
00258     // is higher than the base interpolation buffer's, so we'll ask it to
00259     // interpolate points before inserting into the base interpolation buffer.
00260     TEST_FOR_EXCEPTION(
00261         true,std::logic_error,
00262         "Error, the smart interpolation buffer appender is not implemented for importing interpolation buffers with higher order interpolation into interpolation buffers with a lower order of interpolation yet.\n"
00263         );
00264     // 08/09/07 tscoffe:  Note:  you can't use selectPointsInTimeRange
00265     // immediately because you may need to interpolate points before the first
00266     // node inside the range.  I.e. IB_in.getNodes = [... , range.lower(), ... , range.upper(), ... ]
00267   }
00268 }
00269 
00270 /* This is the function from InterpolationBuffer.hpp which includes some of the
00271  * import functionality I want to include above.
00272  
00273 template<class Scalar>
00274 bool InterpolationBuffer<Scalar>::setRange(
00275   const TimeRange<Scalar>& range,
00276   const InterpolationBufferBase<Scalar>& IB
00277   )
00278 {
00279   const Scalar time_lower = range.lower();
00280   const Scalar time_upper = range.upper();
00281   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00282   Teuchos::OSTab ostab(out,1,"IB::setRange");
00283   if ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_HIGH) ) {
00284     *out << "time_lower = " << time_lower << std::endl;
00285     *out << "time_upper = " << time_upper << std::endl;
00286     *out << "IB = " << IB.description() << std::endl;
00287   }
00288   Array<Scalar> input_nodes;
00289   bool status = IB.getNodes(&input_nodes);
00290   if (!status) { 
00291     return(status);
00292   }
00293   std::sort(input_nodes.begin(),input_nodes.end());
00294   if ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_HIGH) ) {
00295     *out << "input_nodes after sorting = " << std::endl;
00296     for (unsigned int i=0 ; i<input_nodes.size() ; ++i) {
00297       *out << "input_nodes[" << i << "] = " << input_nodes[i] << std::endl;
00298     }
00299   }
00300   // Remove nodes outside the range [time_lower,time_upper]
00301   typename Array<Scalar>::iterator input_it_lower = input_nodes.begin();
00302   for (; input_it_lower != input_nodes.end() ; input_it_lower++) {
00303     if (*input_it_lower >= time_lower) {
00304       break;
00305     }
00306   }
00307   if ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_HIGH) ) {
00308     int n0 = 0;
00309     int n1 = input_it_lower - input_nodes.begin();
00310     *out << "Removing input_nodes before time_lower with indices: [" << n0 << "," << n1 << ")" << std::endl;
00311     for (int i=n0 ; i<n1; ++i) {
00312       *out << "  input_nodes[" << i << "] = " << input_nodes[i] << std::endl;
00313     }
00314   }
00315   // tscoffe 10/19/06 Note:  erase removes the range [it_begin,it_end)
00316   if (input_it_lower - input_nodes.begin() >= 0) {
00317     input_nodes.erase(input_nodes.begin(),input_it_lower);
00318   }
00319   typename Array<Scalar>::iterator input_it_upper = input_nodes.end();
00320   input_it_upper--;
00321   for (; input_it_upper != input_nodes.begin() ; input_it_upper--) {
00322     if (*input_it_upper <= time_upper) {
00323       input_it_upper++;
00324       break;
00325     }
00326   }
00327   // This is to handle the case of one element in the vector
00328   if (*input_it_upper <= time_upper) {
00329     input_it_upper++;
00330   }
00331   if ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_HIGH) ) {
00332     int n0 = input_it_upper - input_nodes.begin();
00333     int n1 = input_nodes.size();
00334     *out << "Removing input_nodes after time_upper with indices [" << n0 << "," << n1 << ")" << std::endl;
00335     for (int i=n0 ; i<n1; ++i) {
00336       *out << "  input_nodes[" << i << "] = " << input_nodes[i] << std::endl;
00337     }
00338   }
00339   if (static_cast<unsigned int>(input_it_upper - input_nodes.begin()) < input_nodes.size()) {
00340     input_nodes.erase(input_it_upper,input_nodes.end());
00341   }
00342   if ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_HIGH) ) {
00343     *out << "input_nodes remaining:" << std::endl;
00344     for (unsigned int i=0 ; i<input_nodes.size() ; ++i) {
00345       *out << "input_nodes[" << i << "] = " << input_nodes[i] << std::endl;
00346     }
00347   }
00348 
00349   // Ask IB to interpolate more points if IB's order is higher than ours
00350   typedef Teuchos::ScalarTraits<Scalar> ST;
00351   Scalar h_safety = Scalar(2*ST::one());
00352   int IBOrder = IB.getOrder();
00353   if (IBOrder >= interpolator->order()) {
00354     std::list<Scalar> add_nodes;
00355     for (unsigned int i=0 ; i<input_nodes.size()-1 ; ++i) {
00356       Scalar h_0 = input_nodes[i+1] - input_nodes[i];
00357       Scalar h = pow(h_0,(IBOrder/interpolator->order())/h_safety);
00358       if ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_HIGH) ) {
00359         *out << "i = " << i << std::endl;
00360         *out << "interpolator->order() = " << interpolator->order() << std::endl;
00361         *out << "IB.getOrder() = " << IB.getOrder() << std::endl;
00362         *out << "h = " << h << std::endl;
00363       }
00364       Scalar N = ceil(h_0/h);
00365       h = Scalar(h_0/N);
00366       if ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_HIGH) ) {
00367         *out << "h_0 = " << h_0 << std::endl;
00368         *out << "N = " << N << std::endl;
00369         *out << "h = " << h << std::endl;
00370         *out << "Inserting an additional " << N-1 << " points to be interpolated:" << std::endl;
00371       }
00372       for (int j=1 ; j<N ; ++j) {
00373         if ( static_cast<int>(this->getVerbLevel()) >= static_cast<int>(Teuchos::VERB_HIGH) ) {
00374           *out << input_nodes[i]+j*h << std::endl;
00375         }
00376         add_nodes.push_back(input_nodes[i]+j*h);
00377       }
00378     }
00379     input_nodes.insert(input_nodes.end(),add_nodes.begin(),add_nodes.end());
00380     std::sort(input_nodes.begin(),input_nodes.end());
00381   }
00382   // If IB's order is lower than ours, then simply grab the node values and continue.
00383   // If IB's order is higher than ours, then grab the node values and ask IB to
00384   // interpolate extra values so that our order of accuracy is approximately
00385   // the same as the other IB's.
00386   // One approach:
00387   // Lets say IB's order is p and our order is r (p>r).
00388   // Given a particular interval with spacing h_0, the order of accuracy of IB is h_0^p
00389   // We want to find a spacing h such that h^r = h_0^p.  Clearly, this is h = h_0^{p/r}.
00390   // Given this new spacing, divide up the interval h_0 into h_0/h subintervals
00391   // and ask for the IB to interpolate points.  This will match basic order of
00392   // accuracy estimates.  Its probably a good idea to include a fudge factor in
00393   // there too.  E.g. h = h_0^{p/r}/fudge.
00394 
00395   // Don't forget to check the interval [time_lower,time_upper].
00396   // Use setPoints and check return value to make sure we observe storage_limit.
00397 
00398   Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > input_x;
00399   Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > input_xdot;
00400   Array<ScalarMag> input_accuracy;
00401   status = IB.getPoints( input_nodes, &input_x, &input_xdot, &input_accuracy );
00402   if (!status) { 
00403     return(status);
00404   }
00405   // We could check that the accuracy meets our criteria here.
00406   status = setPoints( input_nodes, input_x, input_xdot, input_accuracy );
00407   return(status);
00408 }
00409 
00410 */
00411 
00412 
00413 
00414 } // namespace Rythmos
00415 
00416 
00417 #endif //Rythmos_INTERPOLATION_BUFFER_APPENDER_H

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