legendre.hh

#ifndef DUNE_FEM_SPACE_SHAPEFUNCTIONSET_LEGENDRE_HH
#define DUNE_FEM_SPACE_SHAPEFUNCTIONSET_LEGENDRE_HH
 
#include <cassert>
#include <cstddef>
 
#include <algorithm>
#include <array>
#include <iterator>
#include <vector>
 
#include <dune/fem/common/coordinate.hh>
#include <dune/fem/space/shapefunctionset/legendrepolynomials.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    // LegendreShapeFunction
    // ---------------------
 
    template< class FunctionSpace >
    class LegendreShapeFunction
    {
      static_assert( FunctionSpace::dimRange == 1, "FunctionSpace must be scalar (i.e., dimRange = 1)." );
 
    public:
      typedef FunctionSpace FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
      typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
    private:
      static const int dimDomain = FunctionSpaceType::dimDomain;
 
    public:
      LegendreShapeFunction () = default;
 
      template< class MultiIndex >
      explicit LegendreShapeFunction ( const MultiIndex &multiIndex )
      {
        using std::begin;
        using std::end;
 
        auto first = begin( multiIndex );
        auto last = end( multiIndex );
        assert( std::distance( first, last ) == dimDomain );
        std::copy( first, last, multiIndex_.begin() );
      }
 
      int order () const noexcept
      {
        return *std::max_element( multiIndex_.begin(), multiIndex_.end() );
      }
 
      const std::array< int, FunctionSpaceType::dimDomain > &orders () const noexcept
      {
        return multiIndex_;
      }
 
      void evaluate ( const DomainType &x, RangeType &value ) const noexcept
      {
        value[ 0 ] = RangeFieldType( 1 );
        for( int i = 0; i < dimDomain; ++i )
          value[ 0 ] *= LegendrePolynomials::evaluate( multiIndex_[ i ], x[ i ] );
      }
 
      void jacobian ( const DomainType &x, JacobianRangeType &jacobian ) const noexcept
      {
        jacobian = JacobianRangeType( 1 );
        for( int k = 0; k < dimDomain; ++k )
        {
          const RangeFieldType phi = LegendrePolynomials::evaluate( multiIndex_[ k ], x[ k ] );
          const RangeFieldType dphi = LegendrePolynomials::jacobian( multiIndex_[ k ], x[ k ]);
          for( int i = 0; i < dimDomain; ++i )
            jacobian[ 0 ][ i ] *= ( k == i ) ? dphi : phi;
        }
      }
 
      void hessian ( const DomainType &x, HessianRangeType &hessian ) const noexcept
      {
        hessian = HessianRangeType( typename HessianRangeType::value_type( 1 ) );
        for( int k = 0; k < dimDomain; ++k )
        {
          const RangeFieldType phi = LegendrePolynomials::evaluate( multiIndex_[ k ], x[ k ] );
          const RangeFieldType dphi = LegendrePolynomials::jacobian( multiIndex_[ k ], x[ k ] );
          for( int i = 0; i < dimDomain; ++i )
          {
            hessian[ 0 ][ i ][ i ] *= ( k == i ) ? LegendrePolynomials::hessian( multiIndex_[ i ], x[ i ]) : phi;
            for( int j = i+1; j < dimDomain; ++j )
            {
              RangeFieldType tmp = ( k == i || k == j ) ? dphi : phi;
              hessian[ 0 ][ i ][ j ] *= tmp;
              hessian[ 0 ][ j ][ i ] *= tmp;
            }
          }
        }
      }
 
    private:
      std::array< int, dimDomain > multiIndex_;
    };
 
 
 
#ifndef DOXYGEN
 
    namespace __LegendreShapeFunctionSet
    {
 
      // DefaultFactory
      // --------------
 
      template< class FunctionSpace >
      class DefaultFactory
      {
        typedef LegendreShapeFunction< FunctionSpace > ShapeFunctionType;
 
        static const int dimDomain = FunctionSpace::dimDomain;
        typedef std::array< int, dimDomain > MultiIndex;
 
      public:
        explicit DefaultFactory ( int order )
          : order_( order )
        {}
 
        int order () const noexcept { return order_; }
 
        std::size_t size () const noexcept
        {
          std::size_t size = 1;
          for( int i = 0; i < dimDomain; ++i )
            size *= order() + 1;
          return size;
        }
 
        template< class InputIterator >
        void operator() ( InputIterator first ) const noexcept
        {
          auto function = [&first]( const MultiIndex &multiIndex )
          {
            *first = ShapeFunctionType( multiIndex );
            ++first;
          };
 
          MultiIndex multiIndex;
          fill( multiIndex, function, &multiIndex[ 0 ], dimDomain, order() );
        }
 
      private:
        template< class Function >
        static void fill ( MultiIndex &multiIndex, Function function,
                           int *begin, std::size_t n, int order )
        {
          if( n > 0u )
          {
            for( *begin = 0; *begin <= order; *begin += 1 )
              fill( multiIndex, function, begin+1, n-1, order );
          }
          else
            function( multiIndex );
        }
 
        int order_;
      };
 
    } // namespace __LegendreShapeFunctionSet
 
#endif // #ifndef DOXYGEN
 
 
 
    // LegendreShapeFunctionSet
    // ------------------------
 
    template< class FunctionSpace, bool hierarchicalOrdering = false >
    class LegendreShapeFunctionSet
    {
    protected:
      typedef LegendreShapeFunction< FunctionSpace > ShapeFunctionType;
 
    public:
      typedef FunctionSpace FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
    protected:
      struct Compare
      {
        bool operator() ( const ShapeFunctionType &lhs, const ShapeFunctionType &rhs ) noexcept
        {
          if( lhs.order() != rhs.order() )
            return lhs.order() < rhs.order();
          else
          {
            const auto &a = lhs.orders();
            const auto &b = rhs.orders();
            return std::lexicographical_compare( a.begin(), a.end(), b.begin(), b.end() );
          }
        }
      };
 
    public:
      LegendreShapeFunctionSet () = default;
 
      explicit LegendreShapeFunctionSet ( int order )
        : LegendreShapeFunctionSet( __LegendreShapeFunctionSet::DefaultFactory< FunctionSpaceType >( order ) )
      {}
 
      template< class Factory >
      LegendreShapeFunctionSet ( const Factory &factory )
        : shapeFunctions_( factory.size() ),
          order_( factory.order() )
      {
        factory( shapeFunctions_.begin() );
 
        // if hierarchical ordering was selected sort shape functions
        // according to polynomial order
        if( hierarchicalOrdering )
        {
          std::sort( shapeFunctions_.begin(), shapeFunctions_.end(), Compare() );
        }
      }
 
      int order () const noexcept { return order_; }
 
      std::size_t size () const noexcept { return shapeFunctions_.size(); }
 
      template< class Point, class Functor >
      void evaluateEach ( const Point &x, Functor functor ) const noexcept
      {
        RangeType value;
        std::size_t i = 0u;
        for( const ShapeFunctionType &shapeFunction : shapeFunctions_ )
        {
          shapeFunction.evaluate( coordinate( x ), value );
          functor( i++, value );
        }
      }
 
      template< class Point, class Functor >
      void jacobianEach ( const Point &x, Functor functor ) const noexcept
      {
        JacobianRangeType jacobian;
        std::size_t i = 0u;
        for( const ShapeFunctionType &shapeFunction : shapeFunctions_ )
        {
          shapeFunction.jacobian( coordinate( x ), jacobian );
          functor( i++, jacobian );
        }
      }
 
      template< class Point, class Functor >
      void hessianEach ( const Point &x, Functor functor ) const noexcept
      {
        HessianRangeType hessian;
        std::size_t i = 0u;
        for( const ShapeFunctionType &shapeFunction : shapeFunctions_ )
        {
          shapeFunction.hessian( coordinate( x ), hessian );
          functor( i++, hessian );
        }
      }
 
    protected:
      std::vector< ShapeFunctionType > shapeFunctions_;
      int order_;
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_SPACE_SHAPEFUNCTIONSET_LEGENDRE_HH