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   Array<TimeType>* points_out, 
00102   const Array<TimeType>& points_in,
00103   const TimeRange<TimeType>& range 
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 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         );
00234     }
00235   }
00236 }
00237 
00238 
00239 template<class TimeType>
00240 void Rythmos::selectPointsInTimeRange(
00241     Array<TimeType>* points_out, 
00242     const Array<TimeType>& points_in,
00243     const TimeRange<TimeType>& range 
00244     )
00245 {
00246   points_out->clear();
00247   int Nt = Teuchos::as<int>(points_in.size());
00248   for (int i=0; i < Nt ; ++i) {
00249     if (range.isInRange(points_in[i])) {
00250       points_out->push_back(points_in[i]);
00251     }
00252   }
00253 }
00254 
00255 
00256 template<class TimeType>
00257 void Rythmos::removePointsInTimeRange(
00258     Array<TimeType>* points_in, 
00259     const TimeRange<TimeType>& range 
00260     )
00261 {
00262   Array<TimeType> values_to_remove;
00263   for (int i=0 ; i<Teuchos::as<int>(points_in->size()) ; ++i) {
00264     if (range.isInRange((*points_in)[i])) {
00265       values_to_remove.push_back((*points_in)[i]);
00266     }
00267   }
00268   typename Array<TimeType>::iterator point_it;
00269   for (int i=0 ; i< Teuchos::as<int>(values_to_remove.size()) ; ++i) {
00270     point_it = std::find(points_in->begin(),points_in->end(),values_to_remove[i]);
00271     TEST_FOR_EXCEPTION(
00272         point_it == points_in->end(), std::logic_error,
00273         "Error, point to remove = " << values_to_remove[i] << " not found with std:find!\n"
00274         );
00275     points_in->erase(point_it);
00276   }
00277 }
00278 
00279 
00280 template<class Scalar>
00281 bool Rythmos::getCurrentPoints(
00282   const InterpolationBufferBase<Scalar> &interpBuffer,
00283   const Array<Scalar>& time_vec,
00284   Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
00285   Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
00286   int *nextTimePointIndex_inout
00287   )
00288 {
00289 
00290   typedef ScalarTraits<Scalar> ST;
00291   using Teuchos::as;
00292 
00293   const int numTotalTimePoints = time_vec.size();
00294 
00295   // Validate input
00296 #ifdef TEUCHOS_DEBUG
00297   TEST_FOR_EXCEPT(nextTimePointIndex_inout==0);
00298   TEUCHOS_ASSERT( 0 <= *nextTimePointIndex_inout && *nextTimePointIndex_inout < numTotalTimePoints );
00299   TEUCHOS_ASSERT( x_vec == 0 || as<int>(x_vec->size()) == numTotalTimePoints );
00300   TEUCHOS_ASSERT( xdot_vec == 0 || as<int>(xdot_vec->size()) == numTotalTimePoints );
00301 #endif // TEUCHOS_DEBUG
00302 
00303   int &nextTimePointIndex = *nextTimePointIndex_inout;
00304   const int initNextTimePointIndex = nextTimePointIndex;
00305 
00306   const TimeRange<Scalar> currentTimeRange = interpBuffer.getTimeRange();
00307   
00308   if (currentTimeRange.length() >= ST::zero()) {
00309 
00310     // Load a temp array with all of the current time points that fall in the
00311     // current time range.
00312     Array<Scalar> current_time_vec;
00313     int i;
00314     for ( i = 0; i < numTotalTimePoints-nextTimePointIndex; ++i ) {
00315       const Scalar t = time_vec[nextTimePointIndex];
00316 #ifdef TEUCHOS_DEBUG
00317       TEUCHOS_ASSERT( t >= currentTimeRange.lower() );
00318 #endif // TEUCHOS_DEBUG
00319       if ( currentTimeRange.isInRange(t) ) {
00320         ++nextTimePointIndex;
00321         current_time_vec.push_back(t);
00322       }
00323       else {
00324         break;
00325       }
00326     }
00327 #ifdef TEUCHOS_DEBUG
00328     // Here I am just checking that the loop worked as expected with the data
00329     // in the current time range all comming first.
00330     TEUCHOS_ASSERT( nextTimePointIndex-initNextTimePointIndex == i );
00331 #endif
00332 
00333     // Get points in current time range if any such points exist
00334 
00335     const int numCurrentTimePoints = current_time_vec.size();
00336 
00337     if ( numCurrentTimePoints > 0 ) {
00338 
00339       // Get the state(s) for current time points from the stepper and put
00340       // them into temp arrays
00341       Array<RCP<const Thyra::VectorBase<Scalar> > > current_x_vec;
00342       Array<RCP<const Thyra::VectorBase<Scalar> > > current_xdot_vec;
00343       interpBuffer.getPoints(
00344         current_time_vec,
00345         x_vec ? &current_x_vec : 0,
00346         xdot_vec ? &current_xdot_vec : 0,
00347         0 // accuracy_vec
00348         );
00349 
00350       // Copy the gotten x and xdot vectors from the temp arrays to the output
00351       // arrays.
00352       for ( int i = initNextTimePointIndex; i < nextTimePointIndex; ++i ) {
00353         if (x_vec)
00354           (*x_vec)[i] = current_x_vec[i-initNextTimePointIndex];
00355         if (xdot_vec)
00356           (*xdot_vec)[i] = current_xdot_vec[i-initNextTimePointIndex];
00357       }
00358 
00359     }
00360 
00361   }
00362 
00363   return ( nextTimePointIndex == initNextTimePointIndex ? false : true );
00364 
00365 }
00366 
00367 
00368 #endif //Rythmos_INTERPOLATION_BUFFER_HELPERS_HPP

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