polynomials_piecewise.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_polynomials_piecewise_h
#define dealii_polynomials_piecewise_h



#include <deal.II/base/config.h>

#include <deal.II/base/exceptions.h>
#include <deal.II/base/point.h>
#include <deal.II/base/polynomial.h>
#include <deal.II/base/subscriptor.h>

#include <vector>

DEAL_II_NAMESPACE_OPEN

namespace Polynomials
{
  template <typename number>
  class PiecewisePolynomial : public Subscriptor
  {
  public:
    PiecewisePolynomial(const Polynomial<number> &coefficients_on_interval,
                        const unsigned int        n_intervals,
                        const unsigned int        interval,
                        const bool                spans_next_interval);

    number
    value(const number x) const;

    void
    value(const number x, std::vector<number> &values) const;

    void
    value(const number       x,
          const unsigned int n_derivatives,
          number *           values) const;

    unsigned int
    degree() const;

    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int version);

  protected:
    Polynomial<number> polynomial;

    unsigned int n_intervals;

    unsigned int interval;

    bool spans_two_intervals;
  };



  std::vector<PiecewisePolynomial<double>>
  generate_complete_Lagrange_basis_on_subdivisions(
    const unsigned int n_subdivisions,
    const unsigned int base_degree);

} // namespace Polynomials


/* -------------------------- inline functions --------------------- */

namespace Polynomials
{
  template <typename number>
  inline unsigned int
  PiecewisePolynomial<number>::degree() const
  {
    return polynomial.degree();
  }



  template <typename number>
  inline number
  PiecewisePolynomial<number>::value(const number x) const
  {
    AssertIndexRange(interval, n_intervals);
    number y = x;
    // shift polynomial if necessary
    if (n_intervals > 1)
      {
        const number step = 1. / n_intervals;

        // polynomial spans over two intervals
        if (spans_two_intervals == true)
          {
            const number offset = step * interval;
            if (x < offset)
              return 0;
            else if (x > offset + step + step)
              return 0;
            else if (x < offset + step)
              y = x - offset;
            else
              y = offset + step + step - x;
          }
        else
          {
            const number offset = step * interval;
            if (x < offset || x > offset + step)
              return 0;
            else
              y = x - offset;
          }

        return polynomial.value(y);
      }
    else
      return polynomial.value(x);
  }



  template <typename number>
  template <class Archive>
  inline void
  PiecewisePolynomial<number>::serialize(Archive &ar, const unsigned int)
  {
    // forward to serialization function in the base class.
    ar &static_cast<Subscriptor &>(*this);
    ar &polynomial;
    ar &n_intervals;
    ar &interval;
    ar &spans_two_intervals;
  }

} // namespace Polynomials

DEAL_II_NAMESPACE_CLOSE

#endif