Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_CubicSplineInterpolator_def.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_CUBIC_SPLINE_INTERPOLATOR_DEF_H
30#define Rythmos_CUBIC_SPLINE_INTERPOLATOR_DEF_H
31
32#include "Rythmos_CubicSplineInterpolator_decl.hpp"
33#include "Rythmos_InterpolatorBaseHelpers.hpp"
34#include "Thyra_VectorStdOps.hpp"
35#include "Thyra_VectorSpaceBase.hpp"
36#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
37
38namespace Rythmos {
39
40template<class Scalar>
41Teuchos::RCP<Rythmos::CubicSplineInterpolator<Scalar> >
42cubicSplineInterpolator()
43{
44 RCP<CubicSplineInterpolator<Scalar> > csi = Teuchos::rcp(new CubicSplineInterpolator<Scalar>() );
45 return csi;
46}
47
48template<class Scalar>
49void computeCubicSplineCoeff(
50 const typename DataStore<Scalar>::DataStoreVector_t & data,
51 const Ptr<CubicSplineCoeff<Scalar> > & coeffPtr
52 )
53{
54 using Teuchos::outArg;
55 typedef Teuchos::ScalarTraits<Scalar> ST;
56 using Thyra::createMember;
57 TEUCHOS_TEST_FOR_EXCEPTION(
58 (data.size() < 2), std::logic_error,
59 "Error! A minimum of two data points is required for this cubic spline."
60 );
61 // time data in the DataStoreVector should be unique and sorted
62 Array<Scalar> t;
63 Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
64 dataStoreVectorToVector<Scalar>( data, &t, &x_vec, &xdot_vec, NULL );
65#ifdef HAVE_RYTHMOS_DEBUG
66 assertTimePointsAreSorted<Scalar>( t );
67#endif // HAVE_RYTHMOS_DEBUG
68 // 11/18/08 tscoffe: Question: Should I erase everything in coeffPtr or
69 // re-use what I can? For now, I'll erase and create new each time.
70 CubicSplineCoeff<Scalar>& coeff = *coeffPtr;
71 // If there are only two points, then we do something special and just create
72 // a linear polynomial between the points and return.
73 if (t.size() == 2) {
74 coeff.t.clear();
75 coeff.a.clear(); coeff.b.clear(); coeff.c.clear(); coeff.d.clear();
76 coeff.t.reserve(2);
77 coeff.a.reserve(1); coeff.b.reserve(1); coeff.c.reserve(1); coeff.d.reserve(1);
78 coeff.t.push_back(t[0]);
79 coeff.t.push_back(t[1]);
80 coeff.a.push_back(x_vec[0]->clone_v());
81 coeff.b.push_back(createMember(x_vec[0]->space()));
82 coeff.c.push_back(createMember(x_vec[0]->space()));
83 coeff.d.push_back(createMember(x_vec[0]->space()));
84 Scalar h = coeff.t[1] - coeff.t[0];
85 V_StVpStV(outArg(*coeff.b[0]),ST::one()/h,*x_vec[1],-ST::one()/h,*x_vec[0]);
86 V_S(outArg(*coeff.c[0]),ST::zero());
87 V_S(outArg(*coeff.d[0]),ST::zero());
88 return;
89 }
90 // Data objects we'll need:
91 int n = t.length()-1; // Number of intervals
92 coeff.t.clear(); coeff.t.reserve(n+1);
93 coeff.a.clear(); coeff.a.reserve(n+1);
94 coeff.b.clear(); coeff.b.reserve(n);
95 coeff.c.clear(); coeff.c.reserve(n+1);
96 coeff.d.clear(); coeff.d.reserve(n);
97
98 Array<Scalar> h(n);
99 Array<RCP<Thyra::VectorBase<Scalar> > > alpha(n), z(n+1);
100 Array<Scalar> l(n+1), mu(n);
101 for (int i=0 ; i<n ; ++i) {
102 coeff.t.push_back(t[i]);
103 coeff.a.push_back(x_vec[i]->clone_v());
104 coeff.b.push_back(Thyra::createMember(x_vec[0]->space()));
105 coeff.c.push_back(Thyra::createMember(x_vec[0]->space()));
106 coeff.d.push_back(Thyra::createMember(x_vec[0]->space()));
107 alpha[i] = Thyra::createMember(x_vec[0]->space());
108 z[i] = Thyra::createMember(x_vec[0]->space());
109 }
110 coeff.a.push_back(x_vec[n]->clone_v());
111 coeff.t.push_back(t[n]);
112 coeff.c.push_back(Thyra::createMember(x_vec[0]->space()));
113 z[n] = Thyra::createMember(x_vec[0]->space());
114 Scalar zero = ST::zero();
115 Scalar one = ST::one();
116 Scalar two = Scalar(2*ST::one());
117 Scalar three = Scalar(3*ST::one());
118
119 // Algorithm starts here:
120 for (int i=0 ; i<n ; ++i) {
121 h[i] = coeff.t[i+1]-coeff.t[i];
122 }
123 for (int i=1 ; i<n ; ++i) {
124 V_StVpStV(outArg(*(alpha[i])),three/h[i],*coeff.a[i+1],-3/h[i],*coeff.a[i]);
125 Vp_StV(outArg(*(alpha[i])),-three/h[i-1],*coeff.a[i]);
126 Vp_StV(outArg(*(alpha[i])),+three/h[i-1],*coeff.a[i-1]);
127 }
128 l[0] = one;
129 mu[0] = zero;
130 V_S(outArg(*(z[0])),zero);
131 for (int i=1 ; i<n ; ++i) {
132 l[i] = 2*(coeff.t[i+1]-coeff.t[i-1])-h[i-1]*mu[i-1];
133 mu[i] = h[i]/l[i];
134 V_StVpStV(outArg(*(z[i])),one/l[i],*alpha[i],-h[i-1]/l[i],*z[i-1]);
135 }
136 l[n] = one;
137 V_S(outArg(*(z[n])),zero);
138 V_S(outArg(*(coeff.c[n])),zero);
139 for (int j=n-1 ; j >= 0 ; --j) {
140 V_StVpStV(outArg(*(coeff.c[j])),one,*z[j],-mu[j],*coeff.c[j+1]);
141 V_StVpStV(outArg(*(coeff.b[j])),one/h[j],*coeff.a[j+1],-one/h[j],*coeff.a[j]);
142 Vp_StV(outArg(*(coeff.b[j])),-h[j]/three,*coeff.c[j+1]);
143 Vp_StV(outArg(*(coeff.b[j])),-h[j]*two/three,*coeff.c[j]);
144 V_StVpStV(outArg(*(coeff.d[j])),one/(three*h[j]),*coeff.c[j+1],-one/(three*h[j]),*coeff.c[j]);
145 }
146 // Pop the last entry off of a and c to make them the right size.
147 coeff.a.erase(coeff.a.end()-1);
148 coeff.c.erase(coeff.c.end()-1);
149}
150
151
152template<class Scalar>
153void validateCubicSplineCoeff(const CubicSplineCoeff<Scalar>& coeff)
154{
155 int t_n = coeff.t.size();
156 int a_n = coeff.a.size();
157 int b_n = coeff.b.size();
158 int c_n = coeff.c.size();
159 int d_n = coeff.d.size();
160 TEUCHOS_TEST_FOR_EXCEPTION(
161 ((a_n != t_n-1) || (a_n != b_n) || (a_n != c_n) || (a_n != d_n)),
162 std::logic_error,
163 "Error! The sizes of the data structures in the CubicSplineCoeff object do not match"
164 );
165}
166
167
168template<class Scalar>
169void evaluateCubicSpline(
170 const CubicSplineCoeff<Scalar>& coeff,
171 Teuchos::Ordinal j,
172 const Scalar& t,
173 const Ptr<Thyra::VectorBase<Scalar> >& S,
174 const Ptr<Thyra::VectorBase<Scalar> >& Sp,
175 const Ptr<Thyra::VectorBase<Scalar> >& Spp
176 )
177{
178 using Teuchos::outArg;
179 using Teuchos::as;
180 typedef Teuchos::ScalarTraits<Scalar> ST;
181 // Assert preconditions:
182 validateCubicSplineCoeff<Scalar>(coeff);
183 TEUCHOS_TEST_FOR_EXCEPTION( as<Teuchos::Ordinal>(j) >= coeff.a.size(),
184 std::out_of_range, "Error!, j is out of range" );
185
186 Scalar dt = t-coeff.t[j];
187 const Thyra::VectorBase<Scalar>& a = *(coeff.a[j]);
188 const Thyra::VectorBase<Scalar>& b = *(coeff.b[j]);
189 const Thyra::VectorBase<Scalar>& c = *(coeff.c[j]);
190 const Thyra::VectorBase<Scalar>& d = *(coeff.d[j]);
191
192 if (!Teuchos::is_null(S)) {
193 // Evaluate S:
194 //*S = (a) + (b)*dt + (c)*dt*dt + (d)*dt*dt*dt;
195 V_StVpStV(outArg(*S),dt*dt*dt,d,dt*dt,c);
196 Vp_StV(outArg(*S),dt,b);
197 Vp_StV(outArg(*S),ST::one(),a);
198 }
199 if (!Teuchos::is_null(Sp)) {
200 // Evaluate S':
201 //*Sp = (b) + (c)*2*dt + (d)*3*dt*dt;
202 V_StVpStV(outArg(*Sp),Scalar(3*ST::one())*dt*dt,d,Scalar(2*ST::one())*dt,c);
203 Vp_StV(outArg(*Sp),ST::one(),b);
204 }
205 if (!Teuchos::is_null(Spp)) {
206 // Evaluate S'':
207 //*Spp = (c)*2 + (d)*6*dt;
208 V_StVpStV(outArg(*Spp),Scalar(6*ST::one())*dt,d,Scalar(2*ST::one()),c);
209 }
210}
211
212
213
214
215template<class Scalar>
217{
218 splineCoeffComputed_ = false;
219 nodesSet_ = false;
220}
221
222
223template<class Scalar>
225{
226 return true;
227}
228
229
230template<class Scalar>
231RCP<InterpolatorBase<Scalar> >
233{
234 RCP<CubicSplineInterpolator<Scalar> >
235 interpolator = Teuchos::rcp(new CubicSplineInterpolator<Scalar>);
236 if (!is_null(parameterList_))
237 interpolator->parameterList_ = parameterList(*parameterList_);
238 return interpolator;
239}
240
241template<class Scalar>
243 const RCP<const typename DataStore<Scalar>::DataStoreVector_t> & nodesPtr
244 )
245{
246 nodes_ = nodesPtr;
247 nodesSet_ = true;
248 splineCoeffComputed_ = false;
249#ifdef HAVE_RYTHMOS_DEBUG
250 const typename DataStore<Scalar>::DataStoreVector_t & nodes = *nodesPtr;
251 // Copy nodes to internal data structure for verification upon calls to interpolate
252 nodes_copy_ = Teuchos::rcp(new typename DataStore<Scalar>::DataStoreVector_t);
253 nodes_copy_->reserve(nodes.size());
254 for (int i=0 ; i<Teuchos::as<int>(nodes.size()) ; ++i) {
255 nodes_copy_->push_back(*nodes[i].clone());
256 }
257#endif // HAVE_RYTHMOS_DEBUG
258}
259
260template<class Scalar>
262 const Array<Scalar> &t_values,
263 typename DataStore<Scalar>::DataStoreVector_t *data_out
264 ) const
265{
266 using Teuchos::as;
267 using Teuchos::outArg;
268 typedef Teuchos::ScalarTraits<Scalar> ST;
269
270 TEUCHOS_TEST_FOR_EXCEPTION( nodesSet_ == false, std::logic_error,
271 "Error!, setNodes must be called before interpolate"
272 );
273#ifdef HAVE_RYTHMOS_DEBUG
274 // Check that our nodes_ have not changed between the call to setNodes and interpolate
275 assertNodesUnChanged<Scalar>(*nodes_,*nodes_copy_);
276 // Assert that the base interpolator preconditions are satisfied
277 assertBaseInterpolatePreconditions(*nodes_,t_values,data_out);
278#endif // HAVE_RYTHMOS_DEBUG
279
280 // Output info
281 const RCP<FancyOStream> out = this->getOStream();
282 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
283 Teuchos::OSTab ostab(out,1,"CSI::interpolator");
284 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
285 *out << "nodes_:" << std::endl;
286 for (Teuchos::Ordinal i=0 ; i<(*nodes_).size() ; ++i) {
287 *out << "nodes_[" << i << "] = " << std::endl;
288 (*nodes_)[i].describe(*out,Teuchos::VERB_EXTREME);
289 }
290 *out << "t_values = " << std::endl;
291 for (Teuchos::Ordinal i=0 ; i<t_values.size() ; ++i) {
292 *out << "t_values[" << i << "] = " << t_values[i] << std::endl;
293 }
294 }
295
296 data_out->clear();
297
298 // Return immediately if no time points are requested ...
299 if (t_values.size() == 0) {
300 return;
301 }
302
303 if ((*nodes_).size() == 1) {
304 // trivial case of one node. Preconditions assert that t_values[0] ==
305 // (*nodes_)[0].time so we can just pass it out
306 DataStore<Scalar> DS((*nodes_)[0]);
307 data_out->push_back(DS);
308 }
309 else { // (*nodes_).size() >= 2
310 int n = 0; // index into t_values
311 // Loop through all of the time interpolation points in the buffer and
312 // satisfiy all of the requested time points that you find. NOTE: The
313 // loop will be existed once all of the time points are satisified (see
314 // return below).
315 for (Teuchos::Ordinal i=0 ; i < (*nodes_).size()-1; ++i) {
316 const Scalar& ti = (*nodes_)[i].time;
317 const Scalar& tip1 = (*nodes_)[i+1].time;
318 const TimeRange<Scalar> range_i(ti,tip1);
319 // For the interpolation range of [ti,tip1], satisify all of the
320 // requested points in this range.
321 while ( range_i.isInRange(t_values[n]) ) {
322 // First we check for exact node matches:
323 if (compareTimeValues(t_values[n],ti)==0) {
324 DataStore<Scalar> DS((*nodes_)[i]);
325 data_out->push_back(DS);
326 }
327 else if (compareTimeValues(t_values[n],tip1)==0) {
328 DataStore<Scalar> DS((*nodes_)[i+1]);
329 data_out->push_back(DS);
330 } else {
331 if (!splineCoeffComputed_) {
332 computeCubicSplineCoeff<Scalar>(*nodes_,outArg(splineCoeff_));
333 splineCoeffComputed_ = true;
334 }
335 DataStore<Scalar> DS;
336 RCP<Thyra::VectorBase<Scalar> > x = createMember((*nodes_)[i].x->space());
337 evaluateCubicSpline<Scalar>( splineCoeff_, i, t_values[n], outArg(*x) );
338 DS.time = t_values[n];
339 DS.x = x;
340 DS.accuracy = ST::zero();
341 data_out->push_back(DS);
342 }
343 // Move to the next user time point to consider!
344 n++;
345 if (n == as<int>(t_values.size())) {
346 // WE ARE ALL DONE! MOVE OUT!
347 return;
348 }
349 }
350 // Move on the the next interpolation time range
351 }
352 } // (*nodes_).size() == 1
353}
354
355
356template<class Scalar>
358{
359 return(1);
360}
361
362
363template<class Scalar>
365{
366 std::string name = "Rythmos::CubicSplineInterpolator";
367 return(name);
368}
369
370
371template<class Scalar>
373 FancyOStream &out,
374 const Teuchos::EVerbosityLevel verbLevel
375 ) const
376{
377 using Teuchos::as;
378 if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
379 (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
380 )
381 {
382 out << description() << "::describe" << std::endl;
383 }
384 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW))
385 {}
386 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM))
387 {}
388 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH))
389 {}
390}
391
392
393template <class Scalar>
395 RCP<ParameterList> const& paramList
396 )
397{
398 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
399 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
400 parameterList_ = paramList;
401 Teuchos::readVerboseObjectSublist(&*parameterList_,this);
402}
403
404
405template <class Scalar>
406RCP<ParameterList>
408{
409 return(parameterList_);
410}
411
412
413template <class Scalar>
414RCP<ParameterList>
416{
417 RCP<ParameterList> temp_param_list;
418 std::swap( temp_param_list, parameterList_ );
419 return(temp_param_list);
420}
421
422template<class Scalar>
423RCP<const Teuchos::ParameterList>
425{
426 static RCP<Teuchos::ParameterList> validPL;
427 if (is_null(validPL)) {
428 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
429 Teuchos::setupVerboseObjectSublist(&*pl);
430 validPL = pl;
431 }
432 return (validPL);
433}
434
435
436//
437// Explicit Instantiation macro
438//
439// Must be expanded from within the Rythmos namespace!
440//
441
442
443#define RYTHMOS_CUBIC_SPLINE_INTERPOLATOR_INSTANT(SCALAR) \
444 \
445 template class CubicSplineInterpolator< SCALAR >; \
446 \
447 template class CubicSplineCoeff< SCALAR >; \
448 template RCP<CubicSplineInterpolator< SCALAR > > cubicSplineInterpolator(); \
449 template void computeCubicSplineCoeff( \
450 const DataStore< SCALAR >::DataStoreVector_t & data, \
451 const Ptr<CubicSplineCoeff< SCALAR > > & coeffPtr \
452 ); \
453 template void validateCubicSplineCoeff(const CubicSplineCoeff< SCALAR >& coeff); \
454 template void evaluateCubicSpline( \
455 const CubicSplineCoeff< SCALAR >& coeff, \
456 Teuchos::Ordinal j, \
457 const SCALAR & t, \
458 const Ptr<Thyra::VectorBase< SCALAR > >& S, \
459 const Ptr<Thyra::VectorBase< SCALAR > >& Sp, \
460 const Ptr<Thyra::VectorBase< SCALAR > >& Spp \
461 );
462
463
464
465} // namespace Rythmos
466
467
468#endif // Rythmos_CUBIC_SPLINE_INTERPOLATOR_DEF_H
Concrete implemenation of InterpolatorBase that implements cubic spline interpolation.
void setParameterList(RCP< ParameterList > const &paramList)
void interpolate(const Array< Scalar > &t_values, typename DataStore< Scalar >::DataStoreVector_t *data_out) const
RCP< InterpolatorBase< Scalar > > cloneInterpolator() const
void describe(FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
RCP< const Teuchos::ParameterList > getValidParameters() const
void setNodes(const RCP< const typename DataStore< Scalar >::DataStoreVector_t > &nodes)
Represent a time range.
virtual bool isInRange(const TimeType &t) const