tensor_product_polynomials.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2016 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.
//
// ---------------------------------------------------------------------

#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/base/polynomials_piecewise.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/table.h>
#include <deal.II/base/std_cxx11/array.h>

DEAL_II_NAMESPACE_OPEN



/* ------------------- TensorProductPolynomials -------------- */


namespace internal
{
  namespace
  {
    template <int dim>
    inline
    void compute_tensor_index(const unsigned int,
                              const unsigned int,
                              const unsigned int,
                              unsigned int      ( &)[dim])
    {
      Assert(false, ExcNotImplemented());
    }

    inline
    void compute_tensor_index(const unsigned int n,
                              const unsigned int ,
                              const unsigned int ,
                              unsigned int       (&indices)[1])
    {
      indices[0] = n;
    }

    inline
    void compute_tensor_index(const unsigned int n,
                              const unsigned int n_pols_0,
                              const unsigned int ,
                              unsigned int       (&indices)[2])
    {
      indices[0] = n % n_pols_0;
      indices[1] = n / n_pols_0;
    }

    inline
    void compute_tensor_index(const unsigned int n,
                              const unsigned int n_pols_0,
                              const unsigned int n_pols_1,
                              unsigned int       (&indices)[3])
    {
      indices[0] = n % n_pols_0;
      indices[1] = (n/n_pols_0) % n_pols_1;
      indices[2] = n / (n_pols_0*n_pols_1);
    }
  }
}



template <int dim, typename PolynomialType>
inline
void
TensorProductPolynomials<dim,PolynomialType>::
compute_index (const unsigned int i,
               unsigned int       (&indices)[(dim > 0 ? dim : 1)]) const
{
  Assert (i<Utilities::fixed_power<dim>(polynomials.size()),ExcInternalError());
  internal::compute_tensor_index(index_map[i], polynomials.size(),
                                 polynomials.size(), indices);
}



template <int dim, typename PolynomialType>
void
TensorProductPolynomials<dim,PolynomialType>::output_indices(std::ostream &out) const
{
  unsigned int ix[dim];
  for (unsigned int i=0; i<n_tensor_pols; ++i)
    {
      compute_index(i,ix);
      out << i << "\t";
      for (unsigned int d=0; d<dim; ++d)
        out << ix[d] << " ";
      out << std::endl;
    }
}



template <int dim, typename PolynomialType>
void
TensorProductPolynomials<dim,PolynomialType>::set_numbering
(const std::vector<unsigned int> &renumber)
{
  Assert(renumber.size()==index_map.size(),
         ExcDimensionMismatch(renumber.size(), index_map.size()));

  index_map=renumber;
  for (unsigned int i=0; i<index_map.size(); ++i)
    index_map_inverse[index_map[i]]=i;
}



template <>
double
TensorProductPolynomials<0,Polynomials::Polynomial<double> >
::compute_value(const unsigned int,
                const Point<0> &) const
{
  Assert (false, ExcNotImplemented());
  return 0;
}



template <int dim, typename PolynomialType>
double
TensorProductPolynomials<dim,PolynomialType>::compute_value
(const unsigned int  i,
 const Point<dim>   &p) const
{
  Assert(dim>0, ExcNotImplemented());

  unsigned int indices[dim];
  compute_index (i, indices);

  double value=1.;
  for (unsigned int d=0; d<dim; ++d)
    value *= polynomials[indices[d]].value(p(d));

  return value;
}



template <int dim, typename PolynomialType>
Tensor<1,dim>
TensorProductPolynomials<dim,PolynomialType>::compute_grad (const unsigned int  i,
                                                            const Point<dim>   &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);

  // compute values and
  // uni-directional derivatives at
  // the given point in each
  // co-ordinate direction
  double v [dim][2];
  {
    std::vector<double> tmp (2);
    for (unsigned int d=0; d<dim; ++d)
      {
        polynomials[indices[d]].value (p(d), tmp);
        v[d][0] = tmp[0];
        v[d][1] = tmp[1];
      }
  }

  Tensor<1,dim> grad;
  for (unsigned int d=0; d<dim; ++d)
    {
      grad[d] = 1.;
      for (unsigned int x=0; x<dim; ++x)
        grad[d] *= v[x][d==x];
    }

  return grad;
}



template <int dim, typename PolynomialType>
Tensor<2,dim>
TensorProductPolynomials<dim,PolynomialType>::compute_grad_grad
(const unsigned int  i,
 const Point<dim>   &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);

  double v [dim][3];
  {
    std::vector<double> tmp (3);
    for (unsigned int d=0; d<dim; ++d)
      {
        polynomials[indices[d]].value (p(d), tmp);
        v[d][0] = tmp[0];
        v[d][1] = tmp[1];
        v[d][2] = tmp[2];
      }
  }

  Tensor<2,dim> grad_grad;
  for (unsigned int d1=0; d1<dim; ++d1)
    for (unsigned int d2=0; d2<dim; ++d2)
      {
        grad_grad[d1][d2] = 1.;
        for (unsigned int x=0; x<dim; ++x)
          {
            unsigned int derivative=0;
            if (d1==x || d2==x)
              {
                if (d1==d2)
                  derivative=2;
                else
                  derivative=1;
              }
            grad_grad[d1][d2] *= v[x][derivative];
          }
      }

  return grad_grad;
}




template <int dim, typename PolynomialType>
void
TensorProductPolynomials<dim,PolynomialType>::
compute (const Point<dim>            &p,
         std::vector<double>         &values,
         std::vector<Tensor<1,dim> > &grads,
         std::vector<Tensor<2,dim> > &grad_grads,
         std::vector<Tensor<3,dim> > &third_derivatives,
         std::vector<Tensor<4,dim> > &fourth_derivatives) const
{
  Assert (values.size()==n_tensor_pols    || values.size()==0,
          ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
  Assert (grads.size()==n_tensor_pols     || grads.size()==0,
          ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
  Assert (grad_grads.size()==n_tensor_pols|| grad_grads.size()==0,
          ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
  Assert (third_derivatives.size()==n_tensor_pols|| third_derivatives.size()==0,
          ExcDimensionMismatch2(third_derivatives.size(), n_tensor_pols, 0));
  Assert (fourth_derivatives.size()==n_tensor_pols|| fourth_derivatives.size()==0,
          ExcDimensionMismatch2(fourth_derivatives.size(), n_tensor_pols, 0));

  const bool update_values     = (values.size()==n_tensor_pols),
             update_grads      = (grads.size()==n_tensor_pols),
             update_grad_grads = (grad_grads.size()==n_tensor_pols),
             update_3rd_derivatives      = (third_derivatives.size()==n_tensor_pols),
             update_4th_derivatives = (fourth_derivatives.size()==n_tensor_pols);

  // check how many
  // values/derivatives we have to
  // compute
  unsigned int n_values_and_derivatives = 0;
  if (update_values)
    n_values_and_derivatives = 1;
  if (update_grads)
    n_values_and_derivatives = 2;
  if (update_grad_grads)
    n_values_and_derivatives = 3;
  if (update_3rd_derivatives)
    n_values_and_derivatives = 4;
  if (update_4th_derivatives)
    n_values_and_derivatives = 5;

  // Provide a shortcut if only values are requested. For this case usually the
  // temporary memory allocation below does not pay off. Also we loop over all
  // tensor polynomials in a peculiar way to avoid all dynamic memory allocation.
  // We need to evaluate the polynomial value n_polynomials*dim times, and
  // need to compute n_polynomials^dim tensor polynomials. Therefore we can split
  // the loop over the tensor polynomials into one loop over n_polynomials,
  // evaluate this polynomial for each dimension and multiply it with the
  // associated n_polynomials^(dim-1) tensor polynomials before moving to the
  // next polynomial.
  if (n_values_and_derivatives == 1)
    {
      const unsigned int n_polynomials = polynomials.size();
      const unsigned int factor = n_tensor_pols/n_polynomials;

      std::fill(values.begin(),values.end(),1.0);

      for (unsigned int polynomial=0; polynomial<n_polynomials; ++polynomial)
        {
          for (unsigned int d=0; d<dim; ++d)
            {
              const double value = polynomials[polynomial].value(p(d));

              for (unsigned int f=0; f<factor; ++f)
                {
                  unsigned int index = 0;
                  switch (d)
                    {
                    case 0:
                      index = polynomial+f*n_polynomials;
                      break;
                    case 1:
                      index = polynomial*n_polynomials + (f/n_polynomials)*n_polynomials*n_polynomials + f%n_polynomials;
                      break;
                    case 2:
                      index = polynomial*factor + f;
                      break;
                    default:
                      AssertThrow(false, ExcNotImplemented());
                    }

                  const unsigned int i = index_map_inverse[index];
                  values[i] *= value;
                }
            }
        }
      return;
    }


  // Compute the values (and derivatives, if
  // necessary) of all polynomials at this
  // evaluation point. To avoid expensive memory allocation,
  // we use a small amount of memory on the stack, and store the
  // result in an array (that has enough
  // fields for any evaluation of values and
  // derivatives, up to the 4th derivative, for up to 20 polynomials).
  // If someone uses a larger number of
  // polynomials, we need to allocate more memory on the heap.
  std_cxx11::array<std_cxx11::array<double,5>, dim> *v;
  std_cxx11::array<std_cxx11::array<std_cxx11::array<double,5>, dim>, 20> small_array;
  std::vector<std_cxx11::array<std_cxx11::array<double,5>, dim> > large_array;

  const unsigned int n_polynomials = polynomials.size();
  if (n_polynomials > 20)
    {
      large_array.resize(n_polynomials);
      v = &large_array[0];
    }
  else
    v = &small_array[0];

  for (unsigned int i=0; i<n_polynomials; ++i)
    for (unsigned int d=0; d<dim; ++d)
      polynomials[i].value(p(d), n_values_and_derivatives, &v[i][d][0]);

  for (unsigned int i=0; i<n_tensor_pols; ++i)
    {
      // first get the
      // one-dimensional indices of
      // this particular tensor
      // product polynomial
      unsigned int indices[dim];
      compute_index (i, indices);

      if (update_values)
        {
          values[i] = 1;
          for (unsigned int x=0; x<dim; ++x)
            values[i] *= v[indices[x]][x][0];
        }

      if (update_grads)
        for (unsigned int d=0; d<dim; ++d)
          {
            grads[i][d] = 1.;
            for (unsigned int x=0; x<dim; ++x)
              grads[i][d] *= v[indices[x]][x][d==x];
          }

      if (update_grad_grads)
        for (unsigned int d1=0; d1<dim; ++d1)
          for (unsigned int d2=0; d2<dim; ++d2)
            {
              grad_grads[i][d1][d2] = 1.;
              for (unsigned int x=0; x<dim; ++x)
                {
                  unsigned int derivative=0;
                  if (d1==x) ++derivative;
                  if (d2==x) ++derivative;

                  grad_grads[i][d1][d2]
                  *= v[indices[x]][x][derivative];
                }
            }

      if (update_3rd_derivatives)
        for (unsigned int d1=0; d1<dim; ++d1)
          for (unsigned int d2=0; d2<dim; ++d2)
            for (unsigned int d3=0; d3<dim; ++d3)
              {
                third_derivatives[i][d1][d2][d3] = 1.;
                for (unsigned int x=0; x<dim; ++x)
                  {
                    unsigned int derivative=0;
                    if (d1==x) ++derivative;
                    if (d2==x) ++derivative;
                    if (d3==x) ++derivative;

                    third_derivatives[i][d1][d2][d3]
                    *= v[indices[x]][x][derivative];
                  }
              }

      if (update_4th_derivatives)
        for (unsigned int d1=0; d1<dim; ++d1)
          for (unsigned int d2=0; d2<dim; ++d2)
            for (unsigned int d3=0; d3<dim; ++d3)
              for (unsigned int d4=0; d4<dim; ++d4)
                {
                  fourth_derivatives[i][d1][d2][d3][d4] = 1.;
                  for (unsigned int x=0; x<dim; ++x)
                    {
                      unsigned int derivative=0;
                      if (d1==x) ++derivative;
                      if (d2==x) ++derivative;
                      if (d3==x) ++derivative;
                      if (d4==x) ++derivative;

                      fourth_derivatives[i][d1][d2][d3][d4]
                      *= v[indices[x]][x][derivative];
                    }
                }
    }
}




/* ------------------- AnisotropicPolynomials -------------- */


template <int dim>
AnisotropicPolynomials<dim>::
AnisotropicPolynomials(const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols)
  :
  polynomials (pols),
  n_tensor_pols(get_n_tensor_pols(pols))
{
  Assert (pols.size() == dim, ExcDimensionMismatch(pols.size(), dim));
  for (unsigned int d=0; d<dim; ++d)
    Assert (pols[d].size() > 0,
            ExcMessage ("The number of polynomials must be larger than zero "
                        "for all coordinate directions."));
}




template <int dim>
void
AnisotropicPolynomials<dim>::
compute_index (const unsigned int i,
               unsigned int       (&indices)[dim]) const
{
#ifdef DEBUG
  unsigned int n_poly = 1;
  for (unsigned int d=0; d<dim; ++d)
    n_poly *= polynomials[d].size();
  Assert (i < n_poly, ExcInternalError());
#endif

  if (dim==1)
    internal::compute_tensor_index(i, polynomials[0].size(),
                                   0 /*not used*/, indices);
  else
    internal::compute_tensor_index(i, polynomials[0].size(),
                                   polynomials[1].size(), indices);
}



template <int dim>
double
AnisotropicPolynomials<dim>::compute_value (const unsigned int i,
                                            const Point<dim> &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);

  double value=1.;
  for (unsigned int d=0; d<dim; ++d)
    value *= polynomials[d][indices[d]].value(p(d));

  return value;
}


template <int dim>
Tensor<1,dim>
AnisotropicPolynomials<dim>::compute_grad (const unsigned int i,
                                           const Point<dim> &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);

  // compute values and
  // uni-directional derivatives at
  // the given point in each
  // co-ordinate direction
  std::vector<std::vector<double> > v(dim, std::vector<double> (2));
  for (unsigned int d=0; d<dim; ++d)
    polynomials[d][indices[d]].value(p(d), v[d]);

  Tensor<1,dim> grad;
  for (unsigned int d=0; d<dim; ++d)
    {
      grad[d] = 1.;
      for (unsigned int x=0; x<dim; ++x)
        grad[d] *= v[x][d==x];
    }

  return grad;
}


template <int dim>
Tensor<2,dim>
AnisotropicPolynomials<dim>::compute_grad_grad (const unsigned int i,
                                                const Point<dim> &p) const
{
  unsigned int indices[dim];
  compute_index (i, indices);

  std::vector<std::vector<double> > v(dim, std::vector<double> (3));
  for (unsigned int d=0; d<dim; ++d)
    polynomials[d][indices[d]].value(p(d), v[d]);

  Tensor<2,dim> grad_grad;
  for (unsigned int d1=0; d1<dim; ++d1)
    for (unsigned int d2=0; d2<dim; ++d2)
      {
        grad_grad[d1][d2] = 1.;
        for (unsigned int x=0; x<dim; ++x)
          {
            unsigned int derivative=0;
            if (d1==x || d2==x)
              {
                if (d1==d2)
                  derivative=2;
                else
                  derivative=1;
              }
            grad_grad[d1][d2] *= v[x][derivative];
          }
      }

  return grad_grad;
}




template <int dim>
void
AnisotropicPolynomials<dim>::
compute (const Point<dim>            &p,
         std::vector<double>         &values,
         std::vector<Tensor<1,dim> > &grads,
         std::vector<Tensor<2,dim> > &grad_grads,
         std::vector<Tensor<3,dim> > &third_derivatives,
         std::vector<Tensor<4,dim> > &fourth_derivatives) const
{
  Assert (values.size()==n_tensor_pols || values.size()==0,
          ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
  Assert (grads.size()==n_tensor_pols|| grads.size()==0,
          ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
  Assert (grad_grads.size()==n_tensor_pols|| grad_grads.size()==0,
          ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
  Assert (third_derivatives.size()==n_tensor_pols|| third_derivatives.size()==0,
          ExcDimensionMismatch2(third_derivatives.size(), n_tensor_pols, 0));
  Assert (fourth_derivatives.size()==n_tensor_pols|| fourth_derivatives.size()==0,
          ExcDimensionMismatch2(fourth_derivatives.size(), n_tensor_pols, 0));

  const bool update_values     = (values.size() == n_tensor_pols),
             update_grads      = (grads.size()==n_tensor_pols),
             update_grad_grads = (grad_grads.size()==n_tensor_pols),
             update_3rd_derivatives      = (third_derivatives.size()==n_tensor_pols),
             update_4th_derivatives = (fourth_derivatives.size()==n_tensor_pols);

  // check how many
  // values/derivatives we have to
  // compute
  unsigned int n_values_and_derivatives = 0;
  if (update_values)
    n_values_and_derivatives = 1;
  if (update_grads)
    n_values_and_derivatives = 2;
  if (update_grad_grads)
    n_values_and_derivatives = 3;
  if (update_3rd_derivatives)
    n_values_and_derivatives = 4;
  if (update_4th_derivatives)
    n_values_and_derivatives = 5;

  // compute the values (and
  // derivatives, if necessary) of
  // all polynomials at this
  // evaluation point
  std::vector<std::vector<std::vector<double> > > v(dim);
  for (unsigned int d=0; d<dim; ++d)
    {
      v[d].resize (polynomials[d].size());
      for (unsigned int i=0; i<polynomials[d].size(); ++i)
        {
          v[d][i].resize (n_values_and_derivatives, 0.);
          polynomials[d][i].value(p(d), v[d][i]);
        };
    }

  for (unsigned int i=0; i<n_tensor_pols; ++i)
    {
      // first get the
      // one-dimensional indices of
      // this particular tensor
      // product polynomial
      unsigned int indices[dim];
      compute_index (i, indices);

      if (update_values)
        {
          values[i] = 1;
          for (unsigned int x=0; x<dim; ++x)
            values[i] *= v[x][indices[x]][0];
        }

      if (update_grads)
        for (unsigned int d=0; d<dim; ++d)
          {
            grads[i][d] = 1.;
            for (unsigned int x=0; x<dim; ++x)
              grads[i][d] *= v[x][indices[x]][d==x ? 1 : 0];
          }

      if (update_grad_grads)
        for (unsigned int d1=0; d1<dim; ++d1)
          for (unsigned int d2=0; d2<dim; ++d2)
            {
              grad_grads[i][d1][d2] = 1.;
              for (unsigned int x=0; x<dim; ++x)
                {
                  unsigned int derivative=0;
                  if (d1==x) ++derivative;
                  if (d2==x) ++derivative;

                  grad_grads[i][d1][d2]
                  *= v[x][indices[x]][derivative];
                }
            }

      if (update_3rd_derivatives)
        for (unsigned int d1=0; d1<dim; ++d1)
          for (unsigned int d2=0; d2<dim; ++d2)
            for (unsigned int d3=0; d3<dim; ++d3)
              {
                third_derivatives[i][d1][d2][d3] = 1.;
                for (unsigned int x=0; x<dim; ++x)
                  {
                    unsigned int derivative=0;
                    if (d1==x) ++derivative;
                    if (d2==x) ++derivative;
                    if (d3==x) ++derivative;

                    third_derivatives[i][d1][d2][d3]
                    *= v[x][indices[x]][derivative];
                  }
              }

      if (update_4th_derivatives)
        for (unsigned int d1=0; d1<dim; ++d1)
          for (unsigned int d2=0; d2<dim; ++d2)
            for (unsigned int d3=0; d3<dim; ++d3)
              for (unsigned int d4=0; d4<dim; ++d4)
                {
                  fourth_derivatives[i][d1][d2][d3][d4] = 1.;
                  for (unsigned int x=0; x<dim; ++x)
                    {
                      unsigned int derivative=0;
                      if (d1==x) ++derivative;
                      if (d2==x) ++derivative;
                      if (d3==x) ++derivative;
                      if (d4==x) ++derivative;

                      fourth_derivatives[i][d1][d2][d3][d4]
                      *= v[x][indices[x]][derivative];
                    }
                }
    }
}



template<int dim>
unsigned int
AnisotropicPolynomials<dim>::n() const
{
  return n_tensor_pols;
}


template <int dim>
unsigned int
AnisotropicPolynomials<dim>::
get_n_tensor_pols (const std::vector<std::vector<Polynomials::Polynomial<double> > > &pols)
{
  unsigned int y = 1;
  for (unsigned int d=0; d<dim; ++d)
    y *= pols[d].size();
  return y;
}



/* ------------------- explicit instantiations -------------- */
template class TensorProductPolynomials<1,Polynomials::Polynomial<double> >;
template class TensorProductPolynomials<2,Polynomials::Polynomial<double> >;
template class TensorProductPolynomials<3,Polynomials::Polynomial<double> >;

template class TensorProductPolynomials<1,Polynomials::PiecewisePolynomial<double> >;
template class TensorProductPolynomials<2,Polynomials::PiecewisePolynomial<double> >;
template class TensorProductPolynomials<3,Polynomials::PiecewisePolynomial<double> >;

template class AnisotropicPolynomials<1>;
template class AnisotropicPolynomials<2>;
template class AnisotropicPolynomials<3>;

DEAL_II_NAMESPACE_CLOSE