Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_InterpolationBufferHelpers.hpp
1//@HEADER
2// ***********************************************************************
3//
4// Rythmos Package
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
25//
26// ***********************************************************************
27//@HEADER
28
29#ifndef Rythmos_INTERPOLATION_BUFFER_HELPERS_HPP
30#define Rythmos_INTERPOLATION_BUFFER_HELPERS_HPP
31
32
33#include "Rythmos_InterpolationBufferBase.hpp"
34#include "Teuchos_Assert.hpp"
35#include "Teuchos_as.hpp"
36
37
38namespace Rythmos {
39
40
45template<class Scalar>
46void assertTimePointsAreSorted(const Array<Scalar>& time_vec);
47
48
66template<class Scalar>
68 const InterpolationBufferBase<Scalar> &interpBuffer,
69 const Array<Scalar>& time_vec,
70 const int &startingTimePointIndex = 0
71 );
72
73
88template<class Scalar>
90 const InterpolationBufferBase<Scalar> &interpBuffer,
91 const Array<Scalar>& time_vec
92 );
93
94
99template<class TimeType>
101 const Array<TimeType>& points_in,
102 const TimeRange<TimeType>& range,
103 const Ptr<Array<TimeType> >& points_out
104 );
105
106
111template<class TimeType>
113 Array<TimeType>* points_in,
114 const TimeRange<TimeType>& range
115 );
116
117
164template<class Scalar>
166 const InterpolationBufferBase<Scalar> &interpBuffer,
167 const Array<Scalar>& time_vec,
168 Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
169 Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
170 int *nextTimePointIndex
171 );
172
173
174} // namespace Rythmos
175
176
177//
178// Implementations
179//
180
181
182template<class Scalar>
183void Rythmos::assertTimePointsAreSorted(const Array<Scalar>& time_vec)
184{
185 const int numTimePoints = time_vec.size();
186 for ( int i = 0; i < numTimePoints-1; ++ i ) {
187 TEUCHOS_TEST_FOR_EXCEPTION(
188 time_vec[i] >= time_vec[i+1], std::logic_error,
189 "Error, the time vector points time_vec["<<i<<"] = " << time_vec[i]
190 << " >= time_vec["<<i+1<<"] = " << time_vec[i+1] << " are not [unique|sorted]!"
191 );
192 }
193}
194
195
196template<class Scalar>
197void Rythmos::assertNoTimePointsBeforeCurrentTimeRange(
198 const InterpolationBufferBase<Scalar> &interpBuffer,
199 const Array<Scalar>& time_vec,
200 const int &/* startingTimePointIndex */
201 )
202{
203 typedef ScalarTraits<Scalar> ST;
204 const int numTimePoints = time_vec.size();
205 const TimeRange<Scalar> currentTimeRange = interpBuffer.getTimeRange();
206 if (currentTimeRange.length() >= ST::zero()) {
207 for ( int i = 0; i < numTimePoints; ++i ) {
208 TEUCHOS_TEST_FOR_EXCEPTION(
209 time_vec[i] < currentTimeRange.lower(), std::out_of_range,
210 "Error, time_vec["<<i<<"] = " << time_vec[i] << " < currentTimeRange.lower() = "
211 << currentTimeRange.lower() << " for " << interpBuffer.description() << "!"
212 );
213 }
214 }
215}
216
217
218template<class Scalar>
219void Rythmos::assertNoTimePointsInsideCurrentTimeRange(
220 const InterpolationBufferBase<Scalar>& interpBuffer,
221 const Array<Scalar>& time_vec
222 )
223{
224 typedef ScalarTraits<Scalar> ST;
225 const int numTimePoints = time_vec.size();
226 const TimeRange<Scalar> currentTimeRange = interpBuffer.getTimeRange();
227 if (currentTimeRange.length() >= ST::zero()) {
228 for ( int i = 0; i < numTimePoints; ++i ) {
229 TEUCHOS_TEST_FOR_EXCEPTION(
230 currentTimeRange.isInRange(time_vec[i]), std::out_of_range,
231 "Error, time_vec["<<i<<"] = " << time_vec[i] << " is in TimeRange of "
232 << interpBuffer.description() << " = ["
233 << currentTimeRange.lower() << "," << currentTimeRange.upper() << "]!"
234 );
235 }
236 }
237}
238
239
240template<class TimeType>
241void Rythmos::selectPointsInTimeRange(
242 const Array<TimeType>& points_in,
243 const TimeRange<TimeType>& range,
244 const Ptr<Array<TimeType> >& points_out
245 )
246{
247 points_out->clear();
248 int Nt = Teuchos::as<int>(points_in.size());
249 for (int i=0; i < Nt ; ++i) {
250 if (range.isInRange(points_in[i])) {
251 points_out->push_back(points_in[i]);
252 }
253 }
254}
255
256
257template<class TimeType>
258void Rythmos::removePointsInTimeRange(
259 Array<TimeType>* points_in,
260 const TimeRange<TimeType>& range
261 )
262{
263 Array<TimeType> values_to_remove;
264 for (int i=0 ; i<Teuchos::as<int>(points_in->size()) ; ++i) {
265 if (range.isInRange((*points_in)[i])) {
266 values_to_remove.push_back((*points_in)[i]);
267 }
268 }
269 typename Array<TimeType>::iterator point_it;
270 for (int i=0 ; i< Teuchos::as<int>(values_to_remove.size()) ; ++i) {
271 point_it = std::find(points_in->begin(),points_in->end(),values_to_remove[i]);
272 TEUCHOS_TEST_FOR_EXCEPTION(
273 point_it == points_in->end(), std::logic_error,
274 "Error, point to remove = " << values_to_remove[i] << " not found with std:find!\n"
275 );
276 points_in->erase(point_it);
277 }
278}
279
280
281template<class Scalar>
282bool Rythmos::getCurrentPoints(
283 const InterpolationBufferBase<Scalar> &interpBuffer,
284 const Array<Scalar>& time_vec,
285 Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec,
286 Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec,
287 int *nextTimePointIndex_inout
288 )
289{
290
291 typedef ScalarTraits<Scalar> ST;
292 using Teuchos::as;
293
294 const int numTotalTimePoints = time_vec.size();
295
296 // Validate input
297#ifdef HAVE_RYTHMOS_DEBUG
298 TEUCHOS_TEST_FOR_EXCEPT(nextTimePointIndex_inout==0);
299 TEUCHOS_ASSERT( 0 <= *nextTimePointIndex_inout && *nextTimePointIndex_inout < numTotalTimePoints );
300 TEUCHOS_ASSERT( x_vec == 0 || as<int>(x_vec->size()) == numTotalTimePoints );
301 TEUCHOS_ASSERT( xdot_vec == 0 || as<int>(xdot_vec->size()) == numTotalTimePoints );
302#endif // HAVE_RYTHMOS_DEBUG
303
304 int &nextTimePointIndex = *nextTimePointIndex_inout;
305 const int initNextTimePointIndex = nextTimePointIndex;
306
307 const TimeRange<Scalar> currentTimeRange = interpBuffer.getTimeRange();
308
309 if (currentTimeRange.length() >= ST::zero()) {
310
311 // Load a temp array with all of the current time points that fall in the
312 // current time range.
313 Array<Scalar> current_time_vec;
314 { // scope for i to remove shadow warning.
315 int i;
316 for ( i = 0; i < numTotalTimePoints-nextTimePointIndex; ++i ) {
317 const Scalar t = time_vec[nextTimePointIndex];
318#ifdef HAVE_RYTHMOS_DEBUG
319 TEUCHOS_ASSERT( t >= currentTimeRange.lower() );
320#endif // HAVE_RYTHMOS_DEBUG
321 if ( currentTimeRange.isInRange(t) ) {
322 ++nextTimePointIndex;
323 current_time_vec.push_back(t);
324 }
325 else {
326 break;
327 }
328 }
329#ifdef HAVE_RYTHMOS_DEBUG
330 // Here I am just checking that the loop worked as expected with the data
331 // in the current time range all comming first.
332 TEUCHOS_ASSERT( nextTimePointIndex-initNextTimePointIndex == i );
333#endif
334 }
335
336 // Get points in current time range if any such points exist
337
338 const int numCurrentTimePoints = current_time_vec.size();
339
340 if ( numCurrentTimePoints > 0 ) {
341
342 // Get the state(s) for current time points from the stepper and put
343 // them into temp arrays
344 Array<RCP<const Thyra::VectorBase<Scalar> > > current_x_vec;
345 Array<RCP<const Thyra::VectorBase<Scalar> > > current_xdot_vec;
346 if (x_vec || xdot_vec) {
347 interpBuffer.getPoints(
348 current_time_vec,
349 x_vec ? &current_x_vec : 0,
350 xdot_vec ? &current_xdot_vec : 0,
351 0 // accuracy_vec
352 );
353 }
354
355 // Copy the gotten x and xdot vectors from the temp arrays to the output
356 // arrays.
357 for ( int i = initNextTimePointIndex; i < nextTimePointIndex; ++i ) {
358 if (x_vec)
359 (*x_vec)[i] = current_x_vec[i-initNextTimePointIndex];
360 if (xdot_vec)
361 (*xdot_vec)[i] = current_xdot_vec[i-initNextTimePointIndex];
362 }
363
364 }
365
366 }
367
368 return ( nextTimePointIndex == initNextTimePointIndex ? false : true );
369
370}
371
372
373#endif //Rythmos_INTERPOLATION_BUFFER_HELPERS_HPP
Base class for an interpolation buffer.
void selectPointsInTimeRange(const Array< TimeType > &points_in, const TimeRange< TimeType > &range, const Ptr< Array< TimeType > > &points_out)
Select points from an Array that sit in a TimeRange.
void assertNoTimePointsInsideCurrentTimeRange(const InterpolationBufferBase< Scalar > &interpBuffer, const Array< Scalar > &time_vec)
Assert that none of the time points fall inside the current time range for an interpolation buffer ob...
void assertTimePointsAreSorted(const Array< Scalar > &time_vec)
Assert that a time point vector is sorted.
bool getCurrentPoints(const InterpolationBufferBase< Scalar > &interpBuffer, const Array< Scalar > &time_vec, Array< RCP< const Thyra::VectorBase< Scalar > > > *x_vec, Array< RCP< const Thyra::VectorBase< Scalar > > > *xdot_vec, int *nextTimePointIndex)
Get time points in the current range of an interpolation buffer object.
void assertNoTimePointsBeforeCurrentTimeRange(const InterpolationBufferBase< Scalar > &interpBuffer, const Array< Scalar > &time_vec, const int &startingTimePointIndex=0)
Assert that none of the time points fall before the current time range for an interpolation buffer ob...
void removePointsInTimeRange(Array< TimeType > *points_in, const TimeRange< TimeType > &range)
Remove points from an Array that sit in a TimeRange.
Represent a time range.