Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
18 namespace Tempus {
19 
20 template<class Scalar>
22 interpolate(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 
32 template<class Scalar>
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 
45 template<class Scalar>
49 {
50  return pl_;
51 }
52 
53 template<class Scalar>
57 {
59  pl_ = Teuchos::null;
60  return tmp;
61 }
62 
63 template<class Scalar>
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 
74 template<class Scalar>
75 void
77 lagrange(const int p, const Scalar& t, SolutionState<Scalar>* state_out) const
78 {
79  using Teuchos::RCP;
80  using Teuchos::is_null;
81  using Thyra::VectorBase;
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
bool is_null(const boost::shared_ptr< T > &p)
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getXDot()
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool is_null(const std::shared_ptr< T > &p)
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getXDotDot()
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void lagrange(const int p, const Scalar &t, SolutionState< Scalar > *state_out) const
virtual Teuchos::RCP< const SolutionStateMetaData< Scalar > > getMetaData() const
bool floating_compare_equals(const Scalar &a, const Scalar &b, const Scalar &scale)
Helper function for comparing times.
void interpolate(const Scalar &t, SolutionState< Scalar > *state_out) const
Perform an interpolation.
#define TEUCHOS_ASSERT(assertion_test)
virtual Teuchos::RCP< Thyra::VectorBase< Scalar > > getX()
Solution state for integrators and steppers.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
virtual void copy(const Teuchos::RCP< const SolutionState< Scalar > > &ss)
This is a deep copy.