transformed.hh

#ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_TRANSFORMED_HH
#define DUNE_FEM_SPACE_BASISFUNCTIONSET_TRANSFORMED_HH
 
// C++ includes
#include <cassert>
#include <cstddef>
 
// dune-geometry includes
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
 
// dune-fem includes
#include <dune/fem/common/coordinate.hh>
#include <dune/fem/space/basisfunctionset/functor.hh>
#include <dune/fem/space/basisfunctionset/transformation.hh>
#include <dune/fem/space/common/functionspace.hh>
 
 
namespace Dune
{
 
  namespace Fem
  {
 
    // TransformedBasisFunctionSet
    // ---------------------------
 
    template< class Entity, class ShapeFunctionSet, class Transformation >
    class TransformedBasisFunctionSet
    {
      typedef TransformedBasisFunctionSet< Entity, ShapeFunctionSet, Transformation > ThisType;
 
    public:
      typedef Entity EntityType;
      typedef ShapeFunctionSet ShapeFunctionSetType;
 
    protected:
      typedef typename EntityType::Geometry GeometryType;
      typedef typename GeometryType::ctype ctype;
      typedef typename GeometryType::JacobianTransposed JacobianTransposed;
 
    public:
      typedef typename ShapeFunctionSetType::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 std::decay_t< decltype( Dune::ReferenceElements< ctype, GeometryType::coorddimension >::general( std::declval< const Dune::GeometryType & >() ) ) > ReferenceElementType;
 
      TransformedBasisFunctionSet ()
        : entity_( nullptr )
      {}
 
      explicit TransformedBasisFunctionSet ( const EntityType &entity, const ShapeFunctionSet &shapeFunctionSet = ShapeFunctionSet() )
        : entity_( &entity ),
          shapeFunctionSet_( shapeFunctionSet )
      {}
 
 
      // Basis Function Set Interface Methods
      // ------------------------------------
 
      int order () const { return shapeFunctionSet().order(); }
 
      std::size_t size () const { return shapeFunctionSet().size(); }
 
      auto referenceElement () const
        -> decltype( Dune::ReferenceElements< ctype, GeometryType::coorddimension >::general( std::declval< const Dune::GeometryType & >() ) )
      {
        return Dune::ReferenceElements< ctype, GeometryType::coorddimension >::general( type() );
      }
 
      template< class QuadratureType, class Vector, class DofVector >
      void axpy ( const QuadratureType &quad, const Vector &values, DofVector &dofs ) const
      {
        // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
        const unsigned int nop = quad.nop();
        for( unsigned int qp = 0; qp < nop; ++qp )
          axpy( quad[ qp ], values[ qp ], dofs );
      }
 
      template< class QuadratureType, class VectorA, class VectorB, class DofVector >
      void axpy ( const QuadratureType &quad, const VectorA &valuesA, const VectorB &valuesB, DofVector &dofs ) const
      {
        // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
        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
      {
        RangeType factor = transformation( coordinate( x ) ).apply_t( valueFactor );
        FunctionalAxpyFunctor< RangeType, DofVector > f( factor, dofs );
        shapeFunctionSet().evaluateEach( x, f );
      }
 
      template< class Point, class DofVector >
      void axpy ( const Point &x, const JacobianRangeType &jacobianFactor, DofVector &dofs ) const
      {
        typedef typename GeometryType::JacobianInverseTransposed GeometryJacobianInverseTransposedType;
        const GeometryType &geo = geometry();
        const GeometryJacobianInverseTransposedType &gjit = geo.jacobianInverseTransposed( coordinate( x ) );
        JacobianRangeType tmpJacobianFactor( RangeFieldType( 0 ) );
        for( int r = 0; r < FunctionSpaceType::dimRange; ++r )
          gjit.mtv( jacobianFactor[ r ], tmpJacobianFactor[ r ] );
 
        tmpJacobianFactor = transformation( coordinate( x ) ).apply_t( tmpJacobianFactor );
        FunctionalAxpyFunctor< JacobianRangeType, DofVector > f( tmpJacobianFactor, dofs );
        shapeFunctionSet().jacobianEach( x, f );
      }
 
      template< class Point, class DofVector >
      void axpy( const Point &x, const HessianRangeType &hessianFactor, DofVector &dofs ) const
      {
        DUNE_THROW( NotImplemented, "hessian axpy for TransformedBasisFunctionSet not implemented." );
      }
 
      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 QuadratureType, class DofVector, class RangeArray >
      void evaluateAll ( const QuadratureType &quad, const DofVector &dofs, RangeArray &ranges ) const
      {
        // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
        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 = 0;
        AxpyFunctor< DofVector, RangeType > f( dofs, value );
        shapeFunctionSet().evaluateEach( x, f );
        value = transformation( coordinate( x ) ).apply( value );
      }
 
      template< class Point, class RangeArray >
      void evaluateAll ( const Point &x, RangeArray &values ) const
      {
        assert( values.size() >= size() );
        AssignFunctor< RangeArray > f( values );
        shapeFunctionSet().evaluateEach( x, f );
        auto transform = transformation( coordinate( x ) );
        for( auto &e : values )
          e = transform.apply( e );
      }
 
      template< class QuadratureType, class DofVector, class JacobianArray >
      void jacobianAll ( const QuadratureType &quad, const DofVector &dofs, JacobianArray &jacobians ) const
      {
        // call axpy method for each entry of the given vector, e.g. rangeVector or jacobianVector
        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
      {
        JacobianRangeType localJacobian( RangeFieldType( 0 ) );
        AxpyFunctor< DofVector, JacobianRangeType > f( dofs, localJacobian );
        shapeFunctionSet().jacobianEach( x, f );
        const GeometryType &geo = geometry();
        JacobianTransformation< GeometryType >( geo, coordinate( x ) )( localJacobian, jacobian );
        jacobian = transformation( coordinate( x ) ).apply( jacobian );
      }
 
      template< class Point, class JacobianRangeArray >
      void jacobianAll ( const Point &x, JacobianRangeArray &jacobians ) const
      {
        assert( jacobians.size() >= size() );
        const GeometryType &geo = geometry();
        typedef JacobianTransformation< GeometryType > JacobianTransformation;
        JacobianTransformation jacobianTransformation( geo, coordinate( x ) );
        AssignFunctor< JacobianRangeArray, JacobianTransformation > f( jacobians, jacobianTransformation );
        shapeFunctionSet().jacobianEach( x, f );
        auto transform = transformation( coordinate( x ) );
        for( auto &jacobian : jacobians )
          jacobian = transform.apply( jacobian );
      }
 
      template< class Point, class DofVector, class HessianRange >
      void hessianAll ( const Point &x, const DofVector &dofs, HessianRange& hessian ) const
      {
        if( ! std::is_same< HessianRange, HessianRangeType >:: value )
          DUNE_THROW( NotImplemented, "hessianAll: HessianRange mismatch!" );
 
        DUNE_THROW( NotImplemented, "hessianAll for TransformedBasisFunctionSet not implemented." );
      }
 
      template< class Point, class HessianRangeArray >
      void hessianAll ( const Point &x, HessianRangeArray &hessians ) const
      {
        DUNE_THROW( NotImplemented, "hessianAll for TransformedBasisFunctionSet not implemented." );
      }
 
      const Entity &entity () const
      {
        assert( entity_ );
        return *entity_;
      }
 
      bool valid () const { return bool(entity_); }
 
      Dune::GeometryType type () const { return entity().type(); }
 
 
      // Non-interface methods
      // ---------------------
 
      const ShapeFunctionSetType &shapeFunctionSet () const { return shapeFunctionSet_; }
 
      Transformation transformation ( const DomainType& x ) const
      {
        return Transformation( geometry(), x );
      }
 
    protected:
      GeometryType geometry () const { return entity().geometry(); }
 
    private:
      const EntityType *entity_;
      ShapeFunctionSetType shapeFunctionSet_;
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_TRANSFORMED_HH