default.hh

#ifndef DUNE_FEM_BASISFUNCTIONSET_DEFAULT_HH
#define DUNE_FEM_BASISFUNCTIONSET_DEFAULT_HH
 
// C++ includes
#include <cassert>
#include <cstddef>
#include <memory>
#include <type_traits>
#include <utility>
#include <optional>
 
// dune-geometry includes
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
 
#include <dune/fem/space/shapefunctionset/caching.hh>
 
// dune-fem includes
#include <dune/fem/space/basisfunctionset/functor.hh>
#include <dune/fem/space/basisfunctionset/transformation.hh>
#include <dune/fem/space/shapefunctionset/caching.hh>
#include <dune/fem/quadrature/cachingpointlist.hh>
#include <dune/fem/space/common/functionspace.hh>
#include <dune/fem/version.hh>
 
#include <dune/fem/space/basisfunctionset/codegen.hh>
#include <dune/fem/space/basisfunctionset/evaluatecaller.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    // DefaultBasisFunctionSet
    // -----------------------
 
    template< class Entity, class ShapeFunctionSet >
    class DefaultBasisFunctionSet
    {
      typedef DefaultBasisFunctionSet< Entity, ShapeFunctionSet > ThisType;
 
    public:
      typedef Entity EntityType;
      typedef ShapeFunctionSet ShapeFunctionSetType;
 
      // if underlying shape function set was created with storage CodegenStorage
      // then this value should be true (see selectcaching.hh)
      static constexpr bool codegenShapeFunctionSet = detail::IsCodegenShapeFunctionSet< ShapeFunctionSetType >::value;
 
    protected:
      typedef typename ShapeFunctionSetType::FunctionSpaceType   LocalFunctionSpaceType;
      typedef typename LocalFunctionSpaceType::JacobianRangeType LocalJacobianRangeType;
      typedef typename LocalFunctionSpaceType::HessianRangeType  LocalHessianRangeType;
 
      typedef typename LocalFunctionSpaceType::RangeFieldType RangeFieldType;
 
      typedef typename EntityType::Geometry Geometry ;
 
      typedef typename Geometry::ctype ctype;
    public:
      //  slight misuse of struct ToLocalFunctionSpace!!!
      typedef typename ToNewDimDomainFunctionSpace< LocalFunctionSpaceType, Geometry::coorddimension > :: Type  FunctionSpaceType;
 
      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;
 
      typedef typename FunctionSpaceType::ScalarFunctionSpaceType ScalarFunctionSpaceType;
 
      typedef typename ScalarFunctionSpaceType::RangeType         ScalarRangeType;
      typedef typename ScalarFunctionSpaceType::JacobianRangeType ScalarJacobianRangeType;
 
      typedef std::decay_t< decltype( Dune::ReferenceElements< ctype, Geometry::coorddimension >::general( std::declval< const Dune::GeometryType & >() ) ) > ReferenceElementType;
 
      static const int pointSetId = detail::SelectPointSetId< ShapeFunctionSetType >::value;
 
      DefaultBasisFunctionSet () {}
 
      explicit DefaultBasisFunctionSet ( const EntityType &entity, const ShapeFunctionSet &shapeFunctionSet = ShapeFunctionSet() )
        : entity_( &entity ),
          shapeFunctionSet_( shapeFunctionSet )
      {
        // Note that this should be geometry_ = entity.geometry()
        // But Dune::Geometries are not assignable ...
        geometry_.reset();
        geometry_.emplace( entity.geometry() );
      }
 
      DefaultBasisFunctionSet ( const DefaultBasisFunctionSet &other )
        : entity_( other.entity_ ),
          shapeFunctionSet_( other.shapeFunctionSet_ )
      {
        // Note that this should be geometry_ = entity.geometry()
        // But Dune::Geometries are not assignable ...
        geometry_.reset();
        if( other.geometry_ )
          geometry_.emplace( other.geometry_.value() );
      }
 
      DefaultBasisFunctionSet &operator= ( const DefaultBasisFunctionSet &other )
      {
        entity_ = other.entity_;
        shapeFunctionSet_ = other.shapeFunctionSet_;
 
        // Note that this should be geometry_ = entity.geometry()
        // But Dune::Geometries are not assignable ...
        geometry_.reset();
        if( other.geometry_ )
          geometry_.emplace( other.geometry_.value() );
        return *this;
      }
 
      // 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, Geometry::coorddimension >::general( std::declval< const Dune::GeometryType & >() ) )
      {
        return Dune::ReferenceElements< ctype, Geometry::coorddimension >::general( type() );
      }
 
      template< class QuadratureType, class Vector, class DofVector >
      void axpy ( const QuadratureType &quad, const Vector &values, DofVector &dofs ) const
      {
        axpyImpl( quad, values, dofs, values[ 0 ] );
      }
 
      template< class QuadratureType, class VectorA, class VectorB, class DofVector >
      void axpy ( const QuadratureType &quad, const VectorA &valuesA, const VectorB &valuesB, DofVector &dofs ) const
      {
        assert( valuesA.size() > 0 );
        assert( valuesB.size() > 0 );
 
        axpyImpl( quad, valuesA, dofs, valuesA[ 0 ] );
        axpyImpl( quad, valuesB, dofs, valuesB[ 0 ] );
      }
 
      template< class Point, class DofVector >
      void axpy ( const Point &x, const RangeType &valueFactor, DofVector &dofs ) const
      {
        FunctionalAxpyFunctor< RangeType, DofVector > f( valueFactor, dofs );
        shapeFunctionSet().evaluateEach( x, f );
      }
 
      template< class Point, class DofVector >
      void axpy ( const Point &x, const JacobianRangeType &jacobianFactor, DofVector &dofs ) const
      {
        typedef typename Geometry::JacobianInverseTransposed GeometryJacobianInverseTransposedType;
        const Geometry &geo = geometry();
        const GeometryJacobianInverseTransposedType &gjit = geo.jacobianInverseTransposed( coordinate( x ) );
        LocalJacobianRangeType tmpJacobianFactor;
        for( int r = 0; r < FunctionSpaceType::dimRange; ++r )
          gjit.mtv( jacobianFactor[ r ], tmpJacobianFactor[ r ] );
 
        FunctionalAxpyFunctor< LocalJacobianRangeType, DofVector > f( tmpJacobianFactor, dofs );
        shapeFunctionSet().jacobianEach( x, f );
      }
 
      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 Point, class DofVector >
      void axpy ( const Point &x, const HessianRangeType &hessianFactor, DofVector &dofs ) const
      {
        typedef typename Geometry::JacobianInverseTransposed GeometryJacobianInverseTransposedType;
        const Geometry &geo = geometry();
        const GeometryJacobianInverseTransposedType &gjit = geo.jacobianInverseTransposed( coordinate( x ) );
        LocalHessianRangeType tmpHessianFactor( RangeFieldType(0) );
        hessianTransformation(gjit.transposed(),hessianFactor,tmpHessianFactor);
 
        FunctionalAxpyFunctor< LocalHessianRangeType, DofVector > f( tmpHessianFactor, dofs );
        shapeFunctionSet().hessianEach( x, f );
      }
 
 
      template< class QuadratureType, class DofVector, class RangeArray >
      void evaluateAll ( const QuadratureType &quad, const DofVector &dofs, RangeArray &ranges ) const
      {
        assert( ranges.size() >= quad.nop() );
 
        // if shape function set supports codegen and quadrature supports caching
        if constexpr ( codegenShapeFunctionSet && std::is_base_of< CachingInterface, QuadratureType > :: value)
        {
          typedef Codegen :: EvaluateCallerInterfaceTraits< QuadratureType, RangeArray, DofVector > Traits;
          typedef Codegen :: EvaluateCallerInterface< Traits > BaseEvaluationType;
 
          // get base function evaluate caller
          const auto& baseEval = BaseEvaluationType::storage( *this, rangeCache( quad ), quad );
 
          // true if implementation exists, otherwise this is a nullptr
          if( baseEval )
          {
            baseEval->evaluateRanges( quad, dofs, ranges );
            return ;
          }
        }
 
        {
          // 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 = RangeType( 0 );
        AxpyFunctor< DofVector, RangeType > f( dofs, value );
        shapeFunctionSet().evaluateEach( x, f );
      }
 
      template< class Point, class RangeArray >
      void evaluateAll ( const Point &x, RangeArray &values ) const
      {
        assert( values.size() >= size() );
        std::fill( values.begin(), values.end(), RangeType(0) );
        AssignFunctor< RangeArray > f( values );
        shapeFunctionSet().evaluateEach( x, f );
      }
 
      template< class QuadratureType, class DofVector, class JacobianArray >
      void jacobianAll ( const QuadratureType &quad, const DofVector &dofs, JacobianArray &jacobians ) const
      {
        assert( jacobians.size() >= quad.nop() );
 
        // if shape function set supports codegen and quadrature supports caching
        if constexpr ( codegenShapeFunctionSet && std::is_base_of< CachingInterface, QuadratureType > :: value)
        {
          typedef Codegen :: EvaluateCallerInterfaceTraits< QuadratureType, JacobianArray, DofVector, Geometry >  Traits;
          typedef Codegen :: EvaluateCallerInterface< Traits > BaseEvaluationType;
 
          // get base function evaluate caller (calls axpyRanges)
          const auto& baseEval = BaseEvaluationType::storage( *this, jacobianCache( quad ), quad );
 
          // true if implementation exists
          if( baseEval )
          {
            // call appropriate axpyJacobian method
            const Geometry &geo = geometry();
            baseEval->evaluateJacobians( quad, geo, dofs, jacobians );
            return ;
          }
        }
 
        {
          // 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
      {
        LocalJacobianRangeType localJacobian( RangeFieldType( 0 ) );
        AxpyFunctor< DofVector, LocalJacobianRangeType > f( dofs, localJacobian );
        shapeFunctionSet().jacobianEach( x, f );
 
        typedef JacobianTransformation< Geometry > Transformation;
        const Geometry &geo = geometry();
        Transformation transformation( geo, coordinate( x ) );
        transformation( localJacobian, jacobian );
      }
 
      template< class Point, class JacobianRangeArray >
      void jacobianAll ( const Point &x, JacobianRangeArray &jacobians ) const
      {
        assert( jacobians.size() >= size() );
        typedef JacobianTransformation< Geometry > Transformation;
        const Geometry &geo = geometry();
 
        Transformation transformation( geo, coordinate( x ) );
        AssignFunctor< JacobianRangeArray, Transformation > f( jacobians, transformation );
        shapeFunctionSet().jacobianEach( x, f );
      }
 
      template< class QuadratureType, class DofVector, class HessianArray >
      void hessianAll ( const QuadratureType &quad, const DofVector &dofs, HessianArray &hessians ) const
      {
        assert( hessians.size() >= quad.nop() );
        // 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 )
        {
          hessianAll( quad[ qp ], dofs, hessians[ qp ] );
        }
      }
 
      template< class Point, class DofVector >
      void hessianAll ( const Point &x, const DofVector &dofs, HessianRangeType &hessian ) const
      {
        LocalHessianRangeType localHessian( typename LocalHessianRangeType::value_type( RangeFieldType( 0 ) ) );
        AxpyFunctor< DofVector, LocalHessianRangeType > f( dofs, localHessian );
        shapeFunctionSet().hessianEach( x, f );
 
        typedef HessianTransformation< Geometry > Transformation;
        const Geometry &geo = geometry();
        Transformation transformation( geo, coordinate( x ) );
        transformation( localHessian, hessian );
      }
 
      template< class Point, class HessianRangeArray >
      void hessianAll ( const Point &x, HessianRangeArray &hessians ) const
      {
        assert( hessians.size() >= size() );
        typedef HessianTransformation< Geometry > Transformation;
        const Geometry &geo = geometry();
        Transformation transformation( geo, coordinate( x ) );
        AssignFunctor< HessianRangeArray, Transformation > f( hessians, transformation );
        shapeFunctionSet().hessianEach( x, f );
      }
 
      const Entity &entity () const
      {
        assert( valid() );
        return *entity_;
      }
 
      bool valid () const { return bool(entity_); }
 
      Dune::GeometryType type () const { return entity().type(); }
 
      // Non-interface methods
      // ---------------------
 
      const ShapeFunctionSetType &shapeFunctionSet () const { return shapeFunctionSet_; }
 
    protected:
      template< class QuadratureType, class RangeArray, class DofVector >
      void axpyImpl ( const QuadratureType &quad, const RangeArray &rangeFactors, DofVector &dofs, const RangeType& ) const
      {
        assert( rangeFactors.size() >= quad.nop() );
 
        // if shape function set supports codegen and quadrature supports caching
        if constexpr ( codegenShapeFunctionSet && std::is_base_of< CachingInterface, QuadratureType > :: value)
        {
          typedef Codegen :: EvaluateCallerInterfaceTraits< QuadratureType, RangeArray, DofVector > Traits;
          typedef Codegen :: EvaluateCallerInterface< Traits > BaseEvaluationType;
 
          // get base function evaluate caller (calls axpyRanges)
          const auto& baseEval = BaseEvaluationType::storage( *this, rangeCache( quad ), quad );
 
          // true if implementation exists
          if( baseEval )
          {
            baseEval->axpyRanges( quad, rangeFactors, dofs );
            return ;
          }
        }
 
        {
          // 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 ], rangeFactors[ qp ], dofs );
          }
        }
      }
 
      template< class QuadratureType, class JacobianArray, class DofVector >
      void axpyImpl ( const QuadratureType &quad, const JacobianArray &jacobianFactors, DofVector &dofs, const JacobianRangeType& ) const
      {
        assert( jacobianFactors.size() >= quad.nop() );
        // if shape function set supports codegen and quadrature supports caching
        if constexpr ( codegenShapeFunctionSet && std::is_base_of< CachingInterface, QuadratureType > :: value)
        {
          typedef Codegen :: EvaluateCallerInterfaceTraits< QuadratureType, JacobianArray, DofVector, Geometry >  Traits;
          typedef Codegen :: EvaluateCallerInterface< Traits > BaseEvaluationType;
 
          // get base function evaluate caller (calls axpyRanges)
          const auto& baseEval = BaseEvaluationType::storage( *this, jacobianCache( quad ), quad );
 
          // true if implementation exists
          if( baseEval )
          {
            // call appropriate axpyRanges method
            const Geometry &geo = geometry();
            baseEval->axpyJacobians( quad, geo, jacobianFactors, dofs );
            return ;
          }
        }
 
        {
          // 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 ], jacobianFactors[ qp ], dofs );
          }
        }
      }
 
      template< class QuadratureType, class HessianArray, class DofVector >
      void axpyImpl ( const QuadratureType &quad, const HessianArray &hessianFactors, DofVector &dofs, const HessianRangeType& ) const
      {
        assert( hessianFactors.size() >= quad.nop() );
        /* TODO add code generation for hessians
        // if shape function set supports codegen and quadrature supports caching
        if constexpr ( codegenShapeFunctionSet && std::is_base_of< CachingInterface, QuadratureType > :: value)
        {
          typedef Codegen :: EvaluateCallerInterfaceTraits< QuadratureType, HessianArray, DofVector, Geometry >  Traits;
          typedef Codegen :: EvaluateCallerInterface< Traits > BaseEvaluationType;
 
          // get base function evaluate caller (calls axpyRanges)
          const auto& baseEval = BaseEvaluationType::storage( *this, hessianCache( quad ), quad );
 
          // true if implementation exists
          if( baseEval )
          {
            // call appropriate axpyRanges method
            const Geometry &geo = geometry();
            baseEval->axpyHessian( quad, geo, hessianFactors, dofs );
            return ;
          }
        }
        */
        {
          const unsigned int nop = quad.nop();
          for( unsigned int qp = 0; qp < nop; ++qp )
          {
            axpy( quad[ qp ], hessianFactors[ qp ], dofs );
          }
        }
      }
 
      template <class QuadratureType>
      const auto& rangeCache( const QuadratureType& quad ) const
      {
        return shapeFunctionSet().scalarShapeFunctionSet().impl().rangeCache( quad );
      }
 
      template <class QuadratureType>
      const auto& jacobianCache( const QuadratureType& quad ) const
      {
        return shapeFunctionSet().scalarShapeFunctionSet().impl().jacobianCache( quad );
      }
 
    protected:
      Geometry geometry () const { return geometry_.value(); }
 
      const EntityType *entity_ = nullptr;
      ShapeFunctionSetType shapeFunctionSet_;
      std::optional< Geometry > geometry_;
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_BASISFUNCTIONSET_DEFAULT_HH