Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_InterpolationBufferHelpers.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_HELPERS_HPP
00030 #define Rythmos_INTERPOLATION_BUFFER_HELPERS_HPP
00031 
00032 
00033 #include "Rythmos_InterpolationBufferBase.hpp"
00034 #include "Teuchos_Assert.hpp"
00035 #include "Teuchos_as.hpp"
00036 
00037 
00038 namespace Rythmos {
00039 
00040 
00045 template<class Scalar>
00046 void assertTimePointsAreSorted(const Array<Scalar>& time_vec);
00047 
00048 
00066 template<class Scalar>
00067 void assertNoTimePointsBeforeCurrentTimeRange(
00068   const InterpolationBufferBase<Scalar> &interpBuffer,
00069   const Array<Scalar>& time_vec,
00070   const int &startingTimePointIndex = 0
00071   );
00072 
00073   
00088 template<class Scalar>
00089 void assertNoTimePointsInsideCurrentTimeRange(
00090   const InterpolationBufferBase<Scalar> &interpBuffer,
00091   const Array<Scalar>& time_vec
00092   );
00093 
00094 
00099 template<class TimeType>
00100 void selectPointsInTimeRange(
00101   const Array<TimeType>& points_in,
00102   const TimeRange<TimeType>& range,
00103   const Ptr<Array<TimeType> >& points_out 
00104   );
00105 
00106 
00111 template<class TimeType>
00112 void removePointsInTimeRange(
00113   Array<TimeType>* points_in, 
00114   const TimeRange<TimeType>& range 
00115   );
00116 
00117 
00164 template<class Scalar>
00165 bool getCurrentPoints(
00166   const InterpolationBufferBase<Scalar> &interpBuffer,
00167   const Array<Scalar>& time_vec,
00168   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00169   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00170   int *nextTimePointIndex
00171   );
00172 
00173 
00174 } // namespace Rythmos
00175 
00176 
00177 //
00178 // Implementations
00179 //
00180 
00181 
00182 template<class Scalar>
00183 void Rythmos::assertTimePointsAreSorted(const Array<Scalar>& time_vec)
00184 {
00185   const int numTimePoints = time_vec.size();
00186   for ( int i = 0; i < numTimePoints-1; ++ i ) {
00187     TEST_FOR_EXCEPTION(
00188       time_vec[i] >= time_vec[i+1], std::logic_error,
00189       "Error, the time vector points time_vec["<<i<<"] = " << time_vec[i]
00190       << " >= time_vec["<<i+1<<"] = " << time_vec[i+1] << " are not [unique|sorted]!"
00191       );
00192   }
00193 }
00194 
00195 
00196 template<class Scalar>
00197 void Rythmos::assertNoTimePointsBeforeCurrentTimeRange(
00198   const InterpolationBufferBase<Scalar> &interpBuffer,
00199   const Array<Scalar>& time_vec,
00200   const int &startingTimePointIndex
00201   )
00202 {
00203   typedef ScalarTraits<Scalar> ST;
00204   const int numTimePoints = time_vec.size();
00205   const TimeRange<Scalar> currentTimeRange = interpBuffer.getTimeRange();
00206   if (currentTimeRange.length() >= ST::zero()) {
00207     for ( int i = 0; i < numTimePoints; ++i ) {
00208       TEST_FOR_EXCEPTION(
00209         time_vec[i] < currentTimeRange.lower(), std::out_of_range,
00210         "Error, time_vec["<<i<<"] = " << time_vec[i] << " < currentTimeRange.lower() = "
00211         << currentTimeRange.lower() << " for " << interpBuffer.description() << "!"
00212         );
00213     }
00214   }
00215 }
00216 
00217 
00218 template<class Scalar>
00219 void Rythmos::assertNoTimePointsInsideCurrentTimeRange(
00220   const InterpolationBufferBase<Scalar>& interpBuffer,
00221   const Array<Scalar>& time_vec
00222   )
00223 {
00224   typedef ScalarTraits<Scalar> ST;
00225   const int numTimePoints = time_vec.size();
00226   const TimeRange<Scalar> currentTimeRange = interpBuffer.getTimeRange();
00227   if (currentTimeRange.length() >= ST::zero()) {
00228     for ( int i = 0; i < numTimePoints; ++i ) {
00229       TEST_FOR_EXCEPTION(
00230         currentTimeRange.isInRange(time_vec[i]), std::out_of_range,
00231         "Error, time_vec["<<i<<"] = " << time_vec[i] << " is in TimeRange of " 
00232         << interpBuffer.description() << " = ["
00233         << currentTimeRange.lower() << "," << currentTimeRange.upper() << "]!"
00234         );
00235     }
00236   }
00237 }
00238 
00239 
00240 template<class TimeType>
00241 void Rythmos::selectPointsInTimeRange(
00242     const Array<TimeType>& points_in,
00243     const TimeRange<TimeType>& range,
00244     const Ptr<Array<TimeType> >& points_out 
00245     )
00246 {
00247   points_out->clear();
00248   int Nt = Teuchos::as<int>(points_in.size());
00249   for (int i=0; i < Nt ; ++i) {
00250     if (range.isInRange(points_in[i])) {
00251       points_out->push_back(points_in[i]);
00252     }
00253   }
00254 }
00255 
00256 
00257 template<class TimeType>
00258 void Rythmos::removePointsInTimeRange(
00259     Array<TimeType>* points_in, 
00260     const TimeRange<TimeType>& range 
00261     )
00262 {
00263   Array<TimeType> values_to_remove;
00264   for (int i=0 ; i<Teuchos::as<int>(points_in->size()) ; ++i) {
00265     if (range.isInRange((*points_in)[i])) {
00266       values_to_remove.push_back((*points_in)[i]);
00267     }
00268   }
00269   typename Array<TimeType>::iterator point_it;
00270   for (int i=0 ; i< Teuchos::as<int>(values_to_remove.size()) ; ++i) {
00271     point_it = std::find(points_in->begin(),points_in->end(),values_to_remove[i]);
00272     TEST_FOR_EXCEPTION(
00273         point_it == points_in->end(), std::logic_error,
00274         "Error, point to remove = " << values_to_remove[i] << " not found with std:find!\n"
00275         );
00276     points_in->erase(point_it);
00277   }
00278 }
00279 
00280 
00281 template<class Scalar>
00282 bool Rythmos::getCurrentPoints(
00283   const InterpolationBufferBase<Scalar> &interpBuffer,
00284   const Array<Scalar>& time_vec,
00285   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00286   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00287   int *nextTimePointIndex_inout
00288   )
00289 {
00290 
00291   typedef ScalarTraits<Scalar> ST;
00292   using Teuchos::as;
00293 
00294   const int numTotalTimePoints = time_vec.size();
00295 
00296   // Validate input
00297 #ifdef RYTHMOS_DEBUG
00298   TEST_FOR_EXCEPT(nextTimePointIndex_inout==0);
00299   TEUCHOS_ASSERT( 0 <= *nextTimePointIndex_inout && *nextTimePointIndex_inout < numTotalTimePoints );
00300   TEUCHOS_ASSERT( x_vec == 0 || as<int>(x_vec->size()) == numTotalTimePoints );
00301   TEUCHOS_ASSERT( xdot_vec == 0 || as<int>(xdot_vec->size()) == numTotalTimePoints );
00302 #endif // RYTHMOS_DEBUG
00303 
00304   int &nextTimePointIndex = *nextTimePointIndex_inout;
00305   const int initNextTimePointIndex = nextTimePointIndex;
00306 
00307   const TimeRange<Scalar> currentTimeRange = interpBuffer.getTimeRange();
00308   
00309   if (currentTimeRange.length() >= ST::zero()) {
00310 
00311     // Load a temp array with all of the current time points that fall in the
00312     // current time range.
00313     Array<Scalar> current_time_vec;
00314     { // scope for i to remove shadow warning.
00315       int i;
00316       for ( i = 0; i < numTotalTimePoints-nextTimePointIndex; ++i ) {
00317         const Scalar t = time_vec[nextTimePointIndex];
00318 #ifdef RYTHMOS_DEBUG
00319         TEUCHOS_ASSERT( t >= currentTimeRange.lower() );
00320 #endif // RYTHMOS_DEBUG
00321         if ( currentTimeRange.isInRange(t) ) {
00322           ++nextTimePointIndex;
00323           current_time_vec.push_back(t);
00324         }
00325         else {
00326           break;
00327         }
00328       }
00329 #ifdef RYTHMOS_DEBUG
00330       // Here I am just checking that the loop worked as expected with the data
00331       // in the current time range all comming first.
00332       TEUCHOS_ASSERT( nextTimePointIndex-initNextTimePointIndex == i );
00333 #endif
00334     }
00335 
00336     // Get points in current time range if any such points exist
00337 
00338     const int numCurrentTimePoints = current_time_vec.size();
00339 
00340     if ( numCurrentTimePoints > 0 ) {
00341 
00342       // Get the state(s) for current time points from the stepper and put
00343       // them into temp arrays
00344       Array<RCP<const Thyra::VectorBase<Scalar> > > current_x_vec;
00345       Array<RCP<const Thyra::VectorBase<Scalar> > > current_xdot_vec;
00346       if (x_vec || xdot_vec) {
00347         interpBuffer.getPoints(
00348           current_time_vec,
00349           x_vec ? &current_x_vec : 0,
00350           xdot_vec ? &current_xdot_vec : 0,
00351           0 // accuracy_vec
00352           );
00353       }
00354 
00355       // Copy the gotten x and xdot vectors from the temp arrays to the output
00356       // arrays.
00357       for ( int i = initNextTimePointIndex; i < nextTimePointIndex; ++i ) {
00358         if (x_vec)
00359           (*x_vec)[i] = current_x_vec[i-initNextTimePointIndex];
00360         if (xdot_vec)
00361           (*xdot_vec)[i] = current_xdot_vec[i-initNextTimePointIndex];
00362       }
00363 
00364     }
00365 
00366   }
00367 
00368   return ( nextTimePointIndex == initNextTimePointIndex ? false : true );
00369 
00370 }
00371 
00372 
00373 #endif //Rythmos_INTERPOLATION_BUFFER_HELPERS_HPP
 All Classes Functions Variables Typedefs Friends