87 const Array<Scalar> &t_values,
88 typename DataStore<Scalar>::DataStoreVector_t *data_out
93 typedef Teuchos::ScalarTraits<Scalar> ST;
95#ifdef HAVE_RYTHMOS_DEBUG
96 assertBaseInterpolatePreconditions(*nodes_,t_values,data_out);
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);
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;
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;
125 if ((*nodes_).size() == 1) {
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;
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;
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;
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);
153 while ( range_i.isInRange(t_values[n]) ) {
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;
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;
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;
185 DataStore<Scalar> DS;
186 const Scalar& t = t_values[n];
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;
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;
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);
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);
216 data_out->push_back(DS);
220 if (n == as<int>(t_values.size())) {