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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/subscriptor.h>
#include <deal.II/base/polynomial.h>
#include <deal.II/base/point.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);

}


/* -------------------------- 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;
  }

}

DEAL_II_NAMESPACE_CLOSE

#endif