Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
Tempus_InterpolatorLagrange_impl.hpp
Go to the documentation of this file.
1// @HEADER
2// ****************************************************************************
3// Tempus: Copyright (2017) Sandia Corporation
4//
5// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6// ****************************************************************************
7// @HEADER
8
9#ifndef Tempus_InterpolatorLagrange_impl_hpp
10#define Tempus_InterpolatorLagrange_impl_hpp
11
12#include <algorithm>
13#include <iterator>
14
15#include "Teuchos_ScalarTraits.hpp"
16#include "Thyra_VectorStdOps.hpp"
17
18namespace Tempus {
19
20template<class Scalar>
22interpolate(const Scalar& t, SolutionState<Scalar>* state_out) const
23{
24 int n = (*nodes_).size();
25 TEUCHOS_ASSERT(n > 0);
26 if (n >= order_+1)
27 lagrange(order_, t, state_out);
28 else
29 lagrange(n-1, t, state_out);
30}
31
32template<class Scalar>
34setParameterList(const Teuchos::RCP<Teuchos::ParameterList>& paramList)
35{
36 pl_ = Teuchos::parameterList();
37 if (paramList == Teuchos::null)
38 *pl_ = *(this->getValidParameters());
39 else
40 *pl_ = *paramList;
41 pl_->validateParametersAndSetDefaults(*this->getValidParameters());
42 order_ = pl_->get<int>("Order");
43}
44
45template<class Scalar>
46Teuchos::RCP<Teuchos::ParameterList>
52
53template<class Scalar>
54Teuchos::RCP<Teuchos::ParameterList>
57{
58 Teuchos::RCP<Teuchos::ParameterList> tmp = pl_;
59 pl_ = Teuchos::null;
60 return tmp;
61}
62
63template<class Scalar>
64Teuchos::RCP<const Teuchos::ParameterList>
67{
68 Teuchos::RCP<Teuchos::ParameterList> tmp = Teuchos::parameterList();
69 tmp->set<std::string>("Interpolator Type", "Lagrange");
70 tmp->set<int>("Order",0);
71 return tmp;
72}
73
74template<class Scalar>
75void
77lagrange(const int p, const Scalar& t, SolutionState<Scalar>* state_out) const
78{
79 using Teuchos::RCP;
80 using Teuchos::is_null;
82 using Thyra::V_StVpStV;
83
84 // Here we assume we have at least p nodes
85 int n = nodes_->size();
86 TEUCHOS_ASSERT(n >= p);
87
88 const Scalar t_begin = (*nodes_)[0]->getTime();
89 const Scalar t_final = (*nodes_)[n-1]->getTime();
90
91 // Find largest index i such that
92 // (*nodes)[i]->getTime() <= t
93 int i;
94 if (t <= t_begin)
95 i = 0;
96 else if (t >= t_final)
97 i = n-1;
98 else {
99 auto it = std::find_if(nodes_->begin(), nodes_->end(),
100 [=](const RCP<SolutionState<Scalar> >& s)
101 {
102 return t <= s->getTime();
103 });
104 i = std::distance(nodes_->begin(), it)-1;
105 }
106 TEUCHOS_ASSERT(i >= 0 && i<=n-1);
107
108 // First we check for exact node matches:
109 if (floating_compare_equals((*nodes_)[i]->getTime(),t,t_final)) {
110 state_out->copy((*nodes_)[i]);
111 return;
112 }
113 if (i < n-1 && floating_compare_equals((*nodes_)[i+1]->getTime(),t,t_final)) {
114 state_out->copy((*nodes_)[i+1]);
115 return;
116 }
117
118 // Put t roughly in the middle of the interpolation window
119 int k = i - p/2;
120 if (k < 0)
121 k = 0;
122 if (k+p+1 > n)
123 k = n-p-1;
124 TEUCHOS_ASSERT(k >= 0 && k+p+1 <= n);
125
126 RCP<VectorBase<Scalar> > x = state_out->getX();
127 RCP<VectorBase<Scalar> > xdot = state_out->getXDot();
128 RCP<VectorBase<Scalar> > xdotdot = state_out->getXDotDot();
129
130 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
131 Thyra::assign(x.ptr(), zero);
132 if (!is_null(xdot))
133 Thyra::assign(xdot.ptr(), zero);
134 if (!is_null(xdotdot))
135 Thyra::assign(xdotdot.ptr(), zero);
136
137 for (int j=0; j<=p; ++j) {
138 const Scalar tj = (*nodes_)[k+j]->getTime();
139 RCP<const VectorBase<Scalar> > xj = (*nodes_)[k+j]->getX();
140 RCP<const VectorBase<Scalar> > xdotj = (*nodes_)[k+j]->getXDot();
141 RCP<const VectorBase<Scalar> > xdotdotj = (*nodes_)[k+j]->getXDotDot();
142
143 Scalar num = 1.0;
144 Scalar den = 1.0;
145 for (int l=0; l<=p; ++l) {
146 if (l != j) {
147 const Scalar tl = (*nodes_)[k+l]->getTime();
148 num *= t-tl;
149 den *= tj-tl;
150 }
151 }
152 const Scalar a = num / den;
153 Thyra::Vp_StV(x.ptr(), a, *xj);
154 if (!is_null(xdot))
155 Thyra::Vp_StV(xdot.ptr(), a, *xdotj);
156 if (!is_null(xdotdot))
157 Thyra::Vp_StV(xdotdot.ptr(), a, *xdotdotj);
158 }
159
160 // Set meta-data for interpolated state
161 state_out->getMetaData()->copy((*nodes_)[i]->getMetaData());
162 state_out->getMetaData()->setTime(t);
163 state_out->getMetaData()->setDt(t-(*nodes_)[i]->getTime());
164 state_out->getMetaData()->setIsInterpolated(true);
165 // What else should we set??
166}
167
168} // namespace Tempus
169
170#endif // Tempus_InterpolatorLagrange_impl_hpp
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
void lagrange(const int p, const Scalar &t, SolutionState< Scalar > *state_out) const
void interpolate(const Scalar &t, SolutionState< Scalar > *state_out) const
Perform an interpolation.
bool floating_compare_equals(const Scalar &a, const Scalar &b, const Scalar &scale)
Helper function for comparing times.