Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
Rythmos_HermiteInterpolator_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_HERMITE_INTERPOLATOR_DEF_H
30#define Rythmos_HERMITE_INTERPOLATOR_DEF_H
31
32#include "Rythmos_HermiteInterpolator_decl.hpp"
33#include "Rythmos_InterpolatorBaseHelpers.hpp"
34#include "Thyra_VectorStdOps.hpp"
35#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
36
37namespace Rythmos {
38
39template<class Scalar>
41{
42 nodes_ = Teuchos::null;
43 parameterList_ = Teuchos::null;
44}
45
46template<class Scalar>
48 const RCP<const typename DataStore<Scalar>::DataStoreVector_t> & nodes
49 )
50{
51 nodes_ = nodes;
52}
53
54template<class Scalar>
56 const Array<Scalar> &t_values
57 ,typename DataStore<Scalar>::DataStoreVector_t *data_out ) const
58{
59
60 //TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Error, ths function is not tested!" );
61
62 typedef Teuchos::ScalarTraits<Scalar> ST;
63#ifdef HAVE_RYTHMOS_DEBUG
64 assertInterpolatePreconditions((*nodes_),t_values,data_out);
65#endif // HAVE_RYTHMOS_DEBUG
66 RCP<Teuchos::FancyOStream> out = this->getOStream();
67 Teuchos::OSTab ostab(out,1,"HI::interpolator");
68 if ( Teuchos::as<int>(this->getVerbLevel()) >= Teuchos::as<int>(Teuchos::VERB_HIGH) ) {
69 *out << "(*nodes_):" << std::endl;
70 for (Teuchos::Ordinal i=0 ; i<(*nodes_).size() ; ++i) {
71 *out << "(*nodes_)[" << i << "] = " << std::endl;
72 (*nodes_)[i].describe(*out,Teuchos::VERB_EXTREME);
73 }
74 *out << "t_values = " << std::endl;
75 for (Teuchos::Ordinal i=0 ; i<t_values.size() ; ++i) {
76 *out << "t_values[" << i << "] = " << t_values[i] << std::endl;
77 }
78 for (Teuchos::Ordinal i=0; i<data_out->size() ; ++i) {
79 *out << "data_out[" << i << "] = " << std::endl;
80 (*data_out)[i].describe(*out,Teuchos::VERB_EXTREME);
81 }
82 }
83 data_out->clear();
84 if (t_values.size() == 0) {
85 return;
86 }
87
88 if ((*nodes_).size() == 1) {
89 // trivial case of one node
90 // preconditions assert that t_values[0] == (*nodes_)[0].time so we can just pass it out
91 DataStore<Scalar> DS((*nodes_)[0]);
92 data_out->push_back(DS);
93 } else {
94 // (*nodes_).size() >= 2
95 int n = 0;
96 for (int i=0 ; i<Teuchos::as<int>((*nodes_).size())-1 ; ++i) {
97 const Scalar& t0 = (*nodes_)[i].time;
98 const Scalar& t1 = (*nodes_)[i+1].time;
99 while ((t0 <= t_values[n]) && (t_values[n] <= t1)) {
100 const Scalar& t = t_values[n];
101 // First we check for exact node matches:
102 if (t == t0) {
103 DataStore<Scalar> DS((*nodes_)[i]);
104 data_out->push_back(DS);
105 } else if (t == t1) {
106 DataStore<Scalar> DS((*nodes_)[i+1]);
107 data_out->push_back(DS);
108 } else {
109 RCP<const Thyra::VectorBase<Scalar> > x0 = (*nodes_)[i ].x;
110 RCP<const Thyra::VectorBase<Scalar> > x1 = (*nodes_)[i+1].x;
111 RCP<const Thyra::VectorBase<Scalar> > xdot0 = (*nodes_)[i ].xdot;
112 RCP<const Thyra::VectorBase<Scalar> > xdot1 = (*nodes_)[i+1].xdot;
113
114 // 10/10/06 tscoffe: this could be expensive:
115 RCP<Thyra::VectorBase<Scalar> > tmp_vec = x0->clone_v();
116 RCP<Thyra::VectorBase<Scalar> > xdot_temp = x1->clone_v();
117 Scalar dt = t1-t0;
118 Scalar dt2 = dt*dt;
119 Scalar t_t0 = t - t0;
120 Scalar t_t1 = t - t1;
121 Scalar tmp_t;
122
123 // Compute numerical divided difference:
124 Thyra::Vt_S(xdot_temp.ptr(),Scalar(ST::one()/dt));
125 Thyra::Vp_StV(xdot_temp.ptr(),Scalar(-ST::one()/dt),*x0);
126
127 // interpolate this point
128 DataStore<Scalar> DS;
129 DS.time = t;
130
131 // H_3(t) = x(t0) + xdot(t0)(t-t0) + ((x(t1)-x(t0))/(t1-t0) - xdot(t0))(t-t0)^2/(t1-t0)
132 // +(xdot(t1) - 2(x(t1)-x(t0))/(t1-t0) + xdot(t0))(t-t0)^2(t-t1)/(t1-t0)^2
133 RCP<Thyra::VectorBase<Scalar> > x_vec = x0->clone_v();
134 Thyra::Vp_StV(x_vec.ptr(),t_t0,*xdot0);
135 tmp_t = t_t0*t_t0/dt;
136 Thyra::V_StVpStV(tmp_vec.ptr(),tmp_t,*xdot_temp,Scalar(-ST::one()*tmp_t),*xdot0);
137 Thyra::Vp_V(x_vec.ptr(),*tmp_vec);
138 tmp_t = t_t0*t_t0*t_t1/dt2;
139 Thyra::V_StVpStV(tmp_vec.ptr(),tmp_t,*xdot1,Scalar(-2*tmp_t),*xdot_temp);
140 Thyra::Vp_StV(tmp_vec.ptr(),tmp_t,*xdot0);
141 Thyra::Vp_V(x_vec.ptr(),*tmp_vec);
142 DS.x = x_vec;
143
144 // H_3'(t) = xdot(t0) + 2*((x(t1)-x(t0))/(t1-t0) - xdot(t0))(t-t0)/(t1-t0)
145 // +(xdot(t1) - 2(x(t1)-x(t0))/(t1-t0) + xdot(t0))[2*(t-t0)(t-t1) + (t-t0)^2]/(t1-t0)^2
146 RCP<Thyra::VectorBase<Scalar> > xdot_vec = xdot0->clone_v();
147 tmp_t = t_t0/dt;
148 Thyra::Vp_StV(xdot_vec.ptr(),Scalar(2*tmp_t),*xdot_temp);
149 Thyra::Vp_StV(xdot_vec.ptr(),Scalar(-2*tmp_t),*xdot0);
150 tmp_t = Scalar((2*t_t0*t_t1+t_t0*t_t0)/dt2);
151 Thyra::V_StVpStV(tmp_vec.ptr(),tmp_t,*xdot1,Scalar(-2*tmp_t),*xdot_temp);
152 Thyra::Vp_StV(tmp_vec.ptr(),tmp_t,*xdot0);
153 Thyra::Vp_V(xdot_vec.ptr(),*tmp_vec);
154 DS.xdot = xdot_vec;
155
156 // Accuracy:
157 // f(x) - H_3(x) = (f^{(3)}(\xi(x))/(4!))(x-x0)^2(x-x1)^2
158 DS.accuracy = (t_t0)*(t_t0)*(t_t1)*(t_t1);
159
160 // Push DataStore object onto vector:
161 data_out->push_back(DS);
162 }
163 n++;
164 if (n == Teuchos::as<int>(t_values.size())) {
165 return;
166 }
167 }
168 }
169 } // (*nodes_).size() == 1
170}
171
172// non-member constructor
173template<class Scalar>
174RCP<HermiteInterpolator<Scalar> > hermiteInterpolator()
175{
176 RCP<HermiteInterpolator<Scalar> > hi = rcp(new HermiteInterpolator<Scalar>() );
177 return hi;
178}
179
180template<class Scalar>
182{
183 return(3);
184}
185
186template<class Scalar>
188{
189 std::string name = "Rythmos::HermiteInterpolator";
190 return(name);
191}
192
193template<class Scalar>
195 Teuchos::FancyOStream &out
196 ,const Teuchos::EVerbosityLevel verbLevel
197 ) const
198{
199 if ( (Teuchos::as<int>(verbLevel) == Teuchos::as<int>(Teuchos::VERB_DEFAULT) ) ||
200 (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW) )
201 )
202 {
203 out << description() << "::describe" << std::endl;
204 }
205 else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_LOW))
206 {}
207 else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_MEDIUM))
208 {}
209 else if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_HIGH))
210 {}
211}
212
213template <class Scalar>
214void HermiteInterpolator<Scalar>::setParameterList(RCP<ParameterList> const& paramList)
215{
216 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
217 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
218 parameterList_ = paramList;
219 Teuchos::readVerboseObjectSublist(&*parameterList_,this);
220}
221
222template <class Scalar>
224{
225 return(parameterList_);
226}
227
228template <class Scalar>
230{
231 RCP<ParameterList> temp_param_list = parameterList_;
232 parameterList_ = Teuchos::null;
233 return(temp_param_list);
234}
235
236template<class Scalar>
237RCP<const Teuchos::ParameterList> HermiteInterpolator<Scalar>::getValidParameters() const
238{
239 static RCP<Teuchos::ParameterList> validPL;
240 if (is_null(validPL)) {
241 RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
242 Teuchos::setupVerboseObjectSublist(&*pl);
243 validPL = pl;
244 }
245 return (validPL);
246}
247
248template<class Scalar>
250 const typename DataStore<Scalar>::DataStoreVector_t &data_in
251 ,const Array<Scalar> &t_values
252 ,typename DataStore<Scalar>::DataStoreVector_t *data_out
253 ) const
254{
255 assertBaseInterpolatePreconditions(data_in,t_values,data_out);
256 for (int i=0; i<Teuchos::as<int>(data_in.size()) ; ++i) {
257 TEUCHOS_TEST_FOR_EXCEPTION(
258 is_null(data_in[i].x), std::logic_error,
259 "Error, data_in[" << i << "].x == Teuchos::null.\n"
260 );
261 TEUCHOS_TEST_FOR_EXCEPTION(
262 is_null(data_in[i].xdot), std::logic_error,
263 "Error, data_in[" << i << "].xdot == Teuchos::null.\n"
264 );
265 }
266}
267
268//
269// Explicit Instantiation macro
270//
271// Must be expanded from within the Rythmos namespace!
272//
273
274#define RYTHMOS_HERMITE_INTERPOLATOR_INSTANT(SCALAR) \
275 \
276 template class HermiteInterpolator< SCALAR >; \
277 \
278 template RCP<HermiteInterpolator< SCALAR > > hermiteInterpolator();
279
280
281} // namespace Rythmos
282
283#endif // Rythmos_HERMITE_INTERPOLATOR_DEF_H