Tempus_InterpolatorLagrange_impl.hpp

Go to the documentation of this file.

// @HEADER
// ****************************************************************************
//                Tempus: Copyright (2017) Sandia Corporation
//
// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
// ****************************************************************************
// @HEADER

#ifndef Tempus_InterpolatorLagrange_impl_hpp
#define Tempus_InterpolatorLagrange_impl_hpp

#include <algorithm>
#include <iterator>

#include "Teuchos_ScalarTraits.hpp"
#include "Thyra_VectorStdOps.hpp"

namespace Tempus {

template<class Scalar>
void InterpolatorLagrange<Scalar>::
interpolate(const Scalar& t, SolutionState<Scalar>* state_out) const
{
  int n = (*nodes_).size();
  TEUCHOS_ASSERT(n > 0);
  if (n >= order_+1)
    lagrange(order_, t, state_out);
  else
    lagrange(n-1, t, state_out);
}

template<class Scalar>
void InterpolatorLagrange<Scalar>::
setParameterList(const Teuchos::RCP<Teuchos::ParameterList>& paramList)
{
  pl_ = Teuchos::parameterList();
  if (paramList == Teuchos::null)
    *pl_ = *(this->getValidParameters());
  else
    *pl_ = *paramList;
  pl_->validateParametersAndSetDefaults(*this->getValidParameters());
  order_ = pl_->get<int>("Order");
}

template<class Scalar>
Teuchos::RCP<Teuchos::ParameterList>
InterpolatorLagrange<Scalar>::
getNonconstParameterList()
{
  return pl_;
}

template<class Scalar>
Teuchos::RCP<Teuchos::ParameterList>
InterpolatorLagrange<Scalar>::
unsetParameterList()
{
  Teuchos::RCP<Teuchos::ParameterList> tmp = pl_;
  pl_ = Teuchos::null;
  return tmp;
}

template<class Scalar>
Teuchos::RCP<const Teuchos::ParameterList>
InterpolatorLagrange<Scalar>::
getValidParameters() const
{
  Teuchos::RCP<Teuchos::ParameterList> tmp = Teuchos::parameterList();
  tmp->set<std::string>("Interpolator Type", "Lagrange");
  tmp->set<int>("Order",0);
  return tmp;
}

template<class Scalar>
void
InterpolatorLagrange<Scalar>::
lagrange(const int p, const Scalar& t, SolutionState<Scalar>* state_out) const
{
  using Teuchos::RCP;
  using Teuchos::is_null;
  using Thyra::VectorBase;
  using Thyra::V_StVpStV;

  // Here we assume we have at least p nodes
  int n = nodes_->size();
  TEUCHOS_ASSERT(n >= p);

  const Scalar t_begin = (*nodes_)[0]->getTime();
  const Scalar t_final = (*nodes_)[n-1]->getTime();

  // Find largest index i such that
  //      (*nodes)[i]->getTime() <= t
  int i;
  if (t <= t_begin)
    i = 0;
  else if (t >= t_final)
    i = n-1;
  else {
    auto it = std::find_if(nodes_->begin(), nodes_->end(),
                           [=](const RCP<SolutionState<Scalar> >& s)
    {
      return t <= s->getTime();
    });
    i = std::distance(nodes_->begin(), it)-1;
  }
  TEUCHOS_ASSERT(i >= 0 && i<=n-1);

  // First we check for exact node matches:
  if (floating_compare_equals((*nodes_)[i]->getTime(),t,t_final)) {
    state_out->copy((*nodes_)[i]);
    return;
  }
  if (i < n-1 && floating_compare_equals((*nodes_)[i+1]->getTime(),t,t_final)) {
    state_out->copy((*nodes_)[i+1]);
    return;
  }

  // Put t roughly in the middle of the interpolation window
  int k = i - p/2;
  if (k < 0)
    k = 0;
  if (k+p+1 > n)
    k = n-p-1;
  TEUCHOS_ASSERT(k >= 0 && k+p+1 <= n);

  RCP<VectorBase<Scalar> > x = state_out->getX();
  RCP<VectorBase<Scalar> > xdot = state_out->getXDot();
  RCP<VectorBase<Scalar> > xdotdot = state_out->getXDotDot();

  Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
  Thyra::assign(x.ptr(), zero);
  if (!is_null(xdot))
    Thyra::assign(xdot.ptr(), zero);
  if (!is_null(xdotdot))
    Thyra::assign(xdotdot.ptr(), zero);

  for (int j=0; j<=p; ++j) {
    const Scalar tj = (*nodes_)[k+j]->getTime();
    RCP<const VectorBase<Scalar> > xj = (*nodes_)[k+j]->getX();
    RCP<const VectorBase<Scalar> > xdotj = (*nodes_)[k+j]->getXDot();
    RCP<const VectorBase<Scalar> > xdotdotj = (*nodes_)[k+j]->getXDotDot();

    Scalar num = 1.0;
    Scalar den = 1.0;
    for (int l=0; l<=p; ++l) {
      if (l != j) {
        const Scalar tl = (*nodes_)[k+l]->getTime();
        num *=  t-tl;
        den *= tj-tl;
      }
    }
    const Scalar a = num / den;
    Thyra::Vp_StV(x.ptr(), a, *xj);
    if (!is_null(xdot))
      Thyra::Vp_StV(xdot.ptr(), a, *xdotj);
    if (!is_null(xdotdot))
      Thyra::Vp_StV(xdotdot.ptr(), a, *xdotdotj);
  }

  // Set meta-data for interpolated state
  state_out->getMetaData()->copy((*nodes_)[i]->getMetaData());
  state_out->getMetaData()->setTime(t);
  state_out->getMetaData()->setDt(t-(*nodes_)[i]->getTime());
  state_out->getMetaData()->setIsInterpolated(true);
  // What else should we set??
}

} // namespace Tempus

#endif // Tempus_InterpolatorLagrange_impl_hpp