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/storage/entitygeometry.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/cachingquadrature.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
  {
    template< class Entity, class ShapeFunctionSet >
    class BasisFunctionSetStorage : public EntityGeometryStorage< Entity >
    {
      typedef EntityGeometryStorage< Entity >  BaseType;
      typedef BasisFunctionSetStorage< Entity, ShapeFunctionSet > ThisType;
 
    public:
      typedef typename BaseType :: EntityType  EntityType;
      typedef typename BaseType :: Geometry    Geometry;
 
      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;
 
      static const int pointSetId = detail::SelectPointSetId< ShapeFunctionSetType >::value;
 
      BasisFunctionSetStorage () {}
 
      explicit BasisFunctionSetStorage( const EntityType &entity, const ShapeFunctionSet &shapeFunctionSet = ShapeFunctionSet() )
        : BaseType( entity ),
          shapeFunctionSet_( shapeFunctionSet )
      {
      }
 
      BasisFunctionSetStorage ( const BasisFunctionSetStorage &other )
        : BaseType( other ),
          shapeFunctionSet_( other.shapeFunctionSet_ )
      {
      }
 
      BasisFunctionSetStorage &operator= ( const BasisFunctionSetStorage &other )
      {
        BaseType::operator=(other);
        shapeFunctionSet_ = other.shapeFunctionSet_;
        return *this;
      }
 
      using BaseType :: entity;
      using BaseType :: valid;
      using BaseType :: geometry;
      using BaseType :: type;
      using BaseType :: referenceElement;
 
      // Non-interface methods
      // ---------------------
 
      const ShapeFunctionSetType &shapeFunctionSet () const { return shapeFunctionSet_; }
 
      // Basis Function Set Interface Methods
      // ------------------------------------
 
      int order () const { return shapeFunctionSet().order(); }
 
      std::size_t size () const { return shapeFunctionSet().size(); }
 
    protected:
      ShapeFunctionSetType shapeFunctionSet_;
    };
 
 
    // DefaultBasisFunctionSet
    // -----------------------
 
    template< class Entity, class ShapeFunctionSet >
    class DefaultBasisFunctionSet
      : public BasisFunctionSetStorage< Entity, ShapeFunctionSet >
    {
      typedef BasisFunctionSetStorage< Entity, ShapeFunctionSet > BaseType;
      typedef DefaultBasisFunctionSet< Entity, ShapeFunctionSet > ThisType;
 
    public:
      typedef typename BaseType::EntityType  EntityType;
 
      typedef typename BaseType::Geometry    Geometry ;
 
      typedef typename Geometry::ctype       ctype;
 
      typedef typename BaseType::ShapeFunctionSetType  ShapeFunctionSetType;
 
      // if underlying shape function set was created with storage CodegenStorage
      // then this value should be true (see selectcaching.hh)
      using BaseType :: codegenShapeFunctionSet;
 
    protected:
      typedef typename ShapeFunctionSetType::FunctionSpaceType   LocalFunctionSpaceType;
      typedef typename LocalFunctionSpaceType::JacobianRangeType LocalJacobianRangeType;
      typedef typename LocalFunctionSpaceType::HessianRangeType  LocalHessianRangeType;
 
      typedef typename LocalFunctionSpaceType::RangeFieldType RangeFieldType;
 
    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() )
        : BaseType( entity, shapeFunctionSet )
      {
      }
 
      DefaultBasisFunctionSet ( const DefaultBasisFunctionSet &other ) = default;
      DefaultBasisFunctionSet &operator= ( const DefaultBasisFunctionSet &other ) = default;
 
      using BaseType :: entity;
      using BaseType :: geometry;
      using BaseType :: type;
      using BaseType :: shapeFunctionSet;
      using BaseType :: order;
      using BaseType :: size;
      using BaseType :: referenceElement;
 
      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
            baseEval->evaluateJacobians( quad, geometry(), 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;
        Transformation transformation( geometry(), 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;
 
        Transformation transformation( geometry(), 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;
        Transformation transformation( geometry(), 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;
        Transformation transformation( geometry(), coordinate( x ) );
        AssignFunctor< HessianRangeArray, Transformation > f( hessians, transformation );
        shapeFunctionSet().hessianEach( x, f );
      }
    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
            baseEval->axpyJacobians( quad, geometry(), 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 );
      }
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_BASISFUNCTIONSET_DEFAULT_HH