Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_LinearInterpolator_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_LINEAR_INTERPOLATOR_DEF_H
30#define Rythmos_LINEAR_INTERPOLATOR_DEF_H
31
32#include "Rythmos_LinearInterpolator_decl.hpp"
33#include "Rythmos_InterpolatorBaseHelpers.hpp"
34#include "Thyra_VectorStdOps.hpp"
35#include "Thyra_VectorSpaceBase.hpp"
36#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
37
38
39namespace Rythmos {
40
41
42// non-member constructor
43template<class Scalar>
44RCP<LinearInterpolator<Scalar> > linearInterpolator()
45{
46 RCP<LinearInterpolator<Scalar> > li = rcp(new LinearInterpolator<Scalar>() );
47 return li;
48}
49
50template<class Scalar>
52{
53 nodes_ = Teuchos::null;
54 parameterList_ = Teuchos::null;
55}
56
57
58template<class Scalar>
60{
61 return true;
62}
63
64
65template<class Scalar>
66RCP<InterpolatorBase<Scalar> >
68{
69 RCP<LinearInterpolator<Scalar> >
70 interpolator = Teuchos::rcp(new LinearInterpolator<Scalar>);
71 if (!is_null(parameterList_))
72 interpolator->parameterList_ = parameterList(*parameterList_);
73 return interpolator;
74}
75
76template<class Scalar>
78 const RCP<const typename DataStore<Scalar>::DataStoreVector_t> & nodes
79 )
80{
81 nodes_ = nodes;
82}
83
84
85template<class Scalar>
87 const Array<Scalar> &t_values,
88 typename DataStore<Scalar>::DataStoreVector_t *data_out
89 ) const
90{
91
92 using Teuchos::as;
93 typedef Teuchos::ScalarTraits<Scalar> ST;
94
95#ifdef HAVE_RYTHMOS_DEBUG
96 assertBaseInterpolatePreconditions(*nodes_,t_values,data_out);
97#endif // HAVE_RYTHMOS_DEBUG
98
99 // Output info
100 const RCP<FancyOStream> out = this->getOStream();
101 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
102 Teuchos::OSTab ostab(out,1,"LI::interpolator");
103 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
104 *out << "nodes_:" << std::endl;
105 for (Teuchos::Ordinal i=0 ; i<(*nodes_).size() ; ++i) {
106 *out << "nodes_[" << i << "] = " << std::endl;
107 (*nodes_)[i].describe(*out,Teuchos::VERB_EXTREME);
108 }
109 *out << "t_values = " << std::endl;
110 for (Teuchos::Ordinal i=0 ; i<t_values.size() ; ++i) {
111 *out << "t_values[" << i << "] = " << t_values[i] << std::endl;
112 }
113 }
114
115 data_out->clear();
116
117 // Return immediatly if not time points are requested ...
118 if (t_values.size() == 0) {
119 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
120 *out << "Returning because no time points were requested" << std::endl;
121 }
122 return;
123 }
124
125 if ((*nodes_).size() == 1) {
126 // trivial case of one node. Preconditions assert that t_values[0] ==
127 // (*nodes_)[0].time so we can just pass it out
128 DataStore<Scalar> DS((*nodes_)[0]);
129 data_out->push_back(DS);
130 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
131 *out << "Only a single node is in the buffer, so preconditions assert that this must be the point requested" << std::endl;
132 }
133 }
134 else { // (*nodes_).size() >= 2
135 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
136 *out << "More than two nodes, looping through the intervals looking for points to interpolate" << std::endl;
137 }
138 int n = 0; // index into t_values
139 // Loop through all of the time interpolation points in the buffer and
140 // satisfiy all of the requested time points that you find. NOTE: The
141 // loop will be existed once all of the time points are satisified (see
142 // return below).
143 for (int i=0 ; i < as<int>((*nodes_).size())-1; ++i) {
144 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
145 *out << "Looking for interval containing: t_values["<<n<<"] = " << t_values[n] << std::endl;
146 }
147 const Scalar& ti = (*nodes_)[i].time;
148 const Scalar& tip1 = (*nodes_)[i+1].time;
149 const Scalar h = tip1-ti;
150 const TimeRange<Scalar> range_i(ti,tip1);
151 // For the interpolation range of [ti,tip1], satisify all of the
152 // requested points in this range.
153 while ( range_i.isInRange(t_values[n]) ) {
154 // First we check for exact node matches:
155 if (compareTimeValues(t_values[n],ti)==0) {
156 DataStore<Scalar> DS((*nodes_)[i]);
157 data_out->push_back(DS);
158 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
159 *out << "Found an exact node match (on left), shallow copying." << std::endl;
160 *out << "Found t_values["<<n<<"] = " << t_values[n] << " on boundary of interval ["<<ti<<","<<tip1<<"]" << std::endl;
161 }
162 }
163 else if (compareTimeValues(t_values[n],tip1)==0) {
164 DataStore<Scalar> DS((*nodes_)[i+1]);
165 data_out->push_back(DS);
166 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
167 *out << "Found an exact node match (on right), shallow copying." << std::endl;
168 *out << "Found t_values["<<n<<"] = " << t_values[n] << " on boundary of interval ["<<ti<<","<<tip1<<"]" << std::endl;
169 }
170 }
171 else {
172 if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
173 *out << "Interpolating a point (creating a new vector)..." << std::endl;
174 *out << "Found t_values["<<n<<"] = " << t_values[n] << " in interior of interval ["<<ti<<","<<tip1<<"]" << std::endl;
175 }
176 // interpolate this point
177 //
178 // x(t) = (t-ti)/(tip1-ti) * xip1 + (1-(t-ti)/(tip1-ti)) * xi
179 //
180 // Above, it is easy to see that:
181 //
182 // x(ti) = xi
183 // x(tip1) = xip1
184 //
185 DataStore<Scalar> DS;
186 const Scalar& t = t_values[n];
187 DS.time = t;
188 // Get the time and interpolation node points
189 RCP<const Thyra::VectorBase<Scalar> > xi = (*nodes_)[i].x;
190 RCP<const Thyra::VectorBase<Scalar> > xip1 = (*nodes_)[i+1].x;
191 RCP<const Thyra::VectorBase<Scalar> > xdoti = (*nodes_)[i].xdot;
192 RCP<const Thyra::VectorBase<Scalar> > xdotip1 = (*nodes_)[i+1].xdot;
193 // Get constants used in interplation
194 const Scalar dt = t-ti;
195 const Scalar dt_over_h = dt / h;
196 const Scalar one_minus_dt_over_h = ST::one() - dt_over_h;
197 // x = dt/h * xip1 + (1-dt/h) * xi
198 RCP<Thyra::VectorBase<Scalar> > x;
199 if (!is_null(xi) && !is_null(xip1)) {
200 x = createMember(xi->space());
201 Thyra::V_StVpStV(x.ptr(),dt_over_h,*xip1,one_minus_dt_over_h,*xi);
202 }
203 DS.x = x;
204 // x = dt/h * xdotip1 + (1-dt/h) * xdoti
205 RCP<Thyra::VectorBase<Scalar> > xdot;
206 if (!is_null(xdoti) && !is_null(xdotip1)) {
207 xdot = createMember(xdoti->space());
208 Thyra::V_StVpStV(xdot.ptr(),dt_over_h,*xdotip1,one_minus_dt_over_h,*xdoti);
209 }
210 DS.xdot = xdot;
211 // Estimate our accuracy ???
212 DS.accuracy = h;
213 // 2007/12/06: rabartl: Above, should the be a relative value of
214 // some type. What does this really mean?
215 // Store this interplation
216 data_out->push_back(DS);
217 }
218 // Move to the next user time point to consider!
219 n++;
220 if (n == as<int>(t_values.size())) {
221 // WE ARE ALL DONE! MOVE OUT!
222 return;
223 }
224 }
225 // Move on the the next interpolation time range
226 }
227 } // (*nodes_).size() == 1
228}
229
230
231template<class Scalar>
233{
234 return(1);
235}
236
237
238template<class Scalar>
240{
241 std::string name = "Rythmos::LinearInterpolator";
242 return(name);
243}
244
245
246template<class Scalar>
248 FancyOStream &out,
249 const Teuchos::EVerbosityLevel verbLevel
250 ) const
251{
252 using Teuchos::as;
253 if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
254 (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) )
255 )
256 {
257 out << description() << "::describe" << std::endl;
258 }
259 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW))
260 {}
261 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM))
262 {}
263 else if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH))
264 {}
265}
266
267
268template <class Scalar>
270 RCP<ParameterList> const& paramList
271 )
272{
273 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
274 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
275 parameterList_ = paramList;
276 Teuchos::readVerboseObjectSublist(&*parameterList_,this);
277}
278
279
280template <class Scalar>
281RCP<ParameterList>
283{
284 return(parameterList_);
285}
286
287
288template <class Scalar>
289RCP<ParameterList>
291{
292 RCP<ParameterList> temp_param_list;
293 std::swap( temp_param_list, parameterList_ );
294 return(temp_param_list);
295}
296
297template<class Scalar>
298RCP<const Teuchos::ParameterList> LinearInterpolator<Scalar>::getValidParameters() const
299{
300 static RCP<Teuchos::ParameterList> validPL;
301 if (is_null(validPL)) {
302 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
303 Teuchos::setupVerboseObjectSublist(&*pl);
304 validPL = pl;
305 }
306 return (validPL);
307}
308
309//
310// Explicit Instantiation macro
311//
312// Must be expanded from within the Rythmos namespace!
313//
314
315#define RYTHMOS_LINEAR_INTERPOLATOR_INSTANT(SCALAR) \
316 \
317 template class LinearInterpolator< SCALAR >; \
318 \
319 template RCP<LinearInterpolator< SCALAR > > linearInterpolator();
320
321} // namespace Rythmos
322
323
324#endif // Rythmos_LINEAR_INTERPOLATOR_DEF_H
Concrete implemenation of InterpolatorBase just just does simple linear interploation.