simple.hh

#ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_SIMPLE_HH
#define DUNE_FEM_SPACE_BASISFUNCTIONSET_SIMPLE_HH
 
#include <cassert>
#include <cstddef>
 
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
 
#include <dune/fem/space/basisfunctionset/functor.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    // SimpleBasisFunctionSet
    // ----------------------
 
    template< class LocalFunctionSet >
    class SimpleBasisFunctionSet
    {
      typedef SimpleBasisFunctionSet< LocalFunctionSet > ThisType;
 
    public:
      typedef LocalFunctionSet LocalFunctionSetType;
 
      typedef typename LocalFunctionSetType::EntityType EntityType;
      typedef typename LocalFunctionSetType::Geometry   Geometry;
      typedef typename LocalFunctionSetType::FunctionSpaceType FunctionSpaceType;
 
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename FunctionSpaceType::RangeFieldType  RangeFieldType;
      typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
 
      typedef Dune::ReferenceElement< typename EntityType::Geometry > ReferenceElementType;
 
      /* default constructor
       *
       * Note: we require a local function set to have a default constructor;
       *       eventually use LocalFunctionSetProxy.
       */
      SimpleBasisFunctionSet () {}
 
      explicit SimpleBasisFunctionSet ( const LocalFunctionSetType &localFunctionSet )
      : localFunctionSet_( localFunctionSet )
      {}
 
      int order () const { return localFunctionSet().order(); }
 
      std::size_t size () const { return localFunctionSet().size(); }
 
      decltype(auto) referenceElement () const
      {
        return Dune::referenceElement< typename EntityType::Geometry::ctype, EntityType::Geometry::coorddimension >( entity().type() );
      }
 
      template< class Quadrature, class Vector, class DofVector >
      void axpy ( const Quadrature &quad, const Vector &values, DofVector &dofs ) const
      {
        const unsigned int nop = quad.nop();
        for( unsigned int qp = 0; qp < nop; ++qp )
          axpy( quad[ qp ], values[ qp ], dofs );
      }
 
      template< class Quadrature, class VectorA, class VectorB, class DofVector >
      void axpy ( const Quadrature &quad, const VectorA &valuesA, const VectorB &valuesB, DofVector &dofs ) const
      {
        const unsigned int nop = quad.nop();
        for( unsigned int qp = 0; qp < nop; ++qp )
        {
          axpy( quad[ qp ], valuesA[ qp ], dofs );
          axpy( quad[ qp ], valuesB[ qp ], dofs );
        }
      }
 
      template< class Point, class DofVector >
      void axpy ( const Point &x, const RangeType &valueFactor, DofVector &dofs ) const
      {
        FunctionalAxpyFunctor< RangeType, DofVector > functor( valueFactor, dofs );
        localFunctionSet().evaluateEach( x, functor );
      }
 
      template< class Point, class DofVector >
      void axpy ( const Point &x, const JacobianRangeType &jacobianFactor, DofVector &dofs ) const
      {
        FunctionalAxpyFunctor< JacobianRangeType, DofVector > functor( jacobianFactor, dofs );
        localFunctionSet().jacobianEach( x, functor );
      }
      template< class Point, class DofVector >
      void axpy ( const Point &x, const HessianRangeType &hessianFactor, DofVector &dofs ) const
      {
        FunctionalAxpyFunctor< HessianRangeType, DofVector > functor( hessianFactor, dofs );
        localFunctionSet().hessianEach( x, functor );
      }
 
      template< class Point, class DofVector >
      void axpy ( const Point &x, const RangeType &valueFactor,
                  const JacobianRangeType &jacobianFactor,
                  DofVector &dofs ) const
      {
        axpy( x, valueFactor, dofs );
        axpy( x, jacobianFactor, dofs );
      }
 
      template< class Quadrature, class DofVector, class RangeArray >
      void evaluateAll ( const Quadrature &quad, const DofVector &dofs, RangeArray &ranges ) const
      {
        const unsigned int nop = quad.nop();
        for( unsigned int qp = 0; qp < nop; ++qp )
          evaluateAll( quad[ qp ], dofs, ranges[ qp ] );
      }
 
      template< class Point, class DofVector >
      void evaluateAll ( const Point &x, const DofVector &dofs, RangeType &value ) const
      {
        value = RangeType( 0 );
        AxpyFunctor< DofVector, RangeType > functor( dofs, value );
        localFunctionSet().evaluateEach( x, functor );
      }
 
      template< class Point, class RangeArray >
      void evaluateAll ( const Point &x, RangeArray &values ) const
      {
        AssignFunctor< RangeArray > functor( values );
        localFunctionSet().evaluateEach( x, functor );
      }
 
      template< class Quadrature, class DofVector, class JacobianRangeArray >
      void jacobianAll ( const Quadrature &quad, const DofVector &dofs, JacobianRangeArray &jacobians ) const
      {
        const unsigned int nop = quad.nop();
        for( unsigned int qp = 0; qp < nop; ++qp )
          jacobianAll( quad[ qp ], dofs, jacobians[ qp ] );
      }
 
      template< class Point, class DofVector >
      void jacobianAll ( const Point &x, const DofVector &dofs, JacobianRangeType &jacobian ) const
      {
        jacobian = JacobianRangeType( 0 );
        AxpyFunctor< DofVector, JacobianRangeType > functor( dofs, jacobian );
        localFunctionSet().jacobianEach( x, functor );
      }
 
      template< class Point, class JacobianRangeArray >
      void jacobianAll ( const Point &x, JacobianRangeArray &jacobians ) const
      {
        AssignFunctor< JacobianRangeArray > functor( jacobians );
        localFunctionSet().jacobianEach( x, functor );
      }
 
      template< class Point, class DofVector >
      void hessianAll ( const Point &x, const DofVector &dofs, HessianRangeType &hessian ) const
      {
        hessian = HessianRangeType( typename HessianRangeType::value_type( typename RangeType::value_type( 0 ) ) );
        AxpyFunctor< DofVector, HessianRangeType > functor( dofs, hessian );
        localFunctionSet().hessianEach( x, functor );
      }
 
      template< class Point, class HessianRangeArray >
      void hessianAll ( const Point &x, HessianRangeArray &hessians ) const
      {
        AssignFunctor< HessianRangeArray > functor( hessians );
        localFunctionSet().hessianEach( x, functor );
      }
 
      const EntityType &entity () const { return localFunctionSet().entity(); }
 
      const Geometry &geometry () const { return localFunctionSet().geometry(); }
 
      // Non-interface methods
      // ---------------------
 
      const LocalFunctionSetType localFunctionSet () const { return localFunctionSet_; }
 
    private:
      LocalFunctionSetType localFunctionSet_;
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_SIMPLE_HH