evaluatecaller.hh

#ifndef DUNE_FEM_EVALUATECALLER_HH
#define DUNE_FEM_EVALUATECALLER_HH
 
#include <cstdlib>
#include <iostream>
#include <memory>
#include <vector>
 
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
 
#include <dune/fem/space/basisfunctionset/evaluatecallerdeclaration.hh>
 
#include <dune/fem/misc/threads/threadsafevalue.hh>
#include <dune/fem/common/utility.hh>
#include <dune/fem/quadrature/quadrature.hh>
 
//#ifdef USE_BASEFUNCTIONSET_CODEGEN
//#ifndef DUNE_FEM_INCLUDE_AUTOGENERATEDCODE_FILENAME_SPEC
// default filename for autogenerated code is simply autogeneratedcode.hh
// this filename is overloaded by the codegeneration for the python modules
//#include <autogeneratedcode.hh>
//#endif
 
#ifdef DUNE_FEM_INCLUDE_AUTOGENERATEDCODE_FILENAME_SPEC
#define CODEGEN_INCLUDEMAXNUMS
// include max number definitions
#include DUNE_FEM_INCLUDE_AUTOGENERATEDCODE_FILENAME_SPEC
#undef CODEGEN_INCLUDEMAXNUMS
#endif
 
//
// pre-define these values for faster compilation
//
#ifndef MAX_NUMBER_OF_QUAD_POINTS
#define MAX_NUMBER_OF_QUAD_POINTS 20
#endif
 
#ifndef MAX_NUMBER_OF_BASE_FCT
#define MAX_NUMBER_OF_BASE_FCT 20
#endif
 
#ifndef MIN_NUMBER_OF_QUAD_POINTS
#define MIN_NUMBER_OF_QUAD_POINTS 1
#endif
 
#ifndef MIN_NUMBER_OF_BASE_FCT
#define MIN_NUMBER_OF_BASE_FCT 1
#endif
 
#include <dune/fem/space/basisfunctionset/evaluatecallerdefaultimpl.hh>
 
#ifdef DUNE_FEM_INCLUDE_AUTOGENERATEDCODE_FILENAME_SPEC
// defined in evaluatecaller.hh
#include DUNE_FEM_INCLUDE_AUTOGENERATEDCODE_FILENAME_SPEC
#endif // DUNE_FEM_INCLUDE_AUTOGENERATEDCODE_FILENAME_SPEC
 
namespace Dune
{
  namespace Fem
  {
    namespace Codegen
    {
 
    // forward declaration
    template <class Traits,
              int quadNop,
              int numBaseFct >
    class EvaluateCaller;
 
    template< class QuadratureImp,
              class FactorImp,
              class LocalDofVectorImp,
              class GeometryImp = EmptyGeometry >
    struct EvaluateCallerInterfaceTraits
    {
      typedef QuadratureImp     QuadratureType;
      typedef FactorImp         FactorType;
      typedef LocalDofVectorImp LocalDofVectorType;
      typedef GeometryImp       Geometry;
    };
 
    template <class Traits,
              class BaseFunctionSet,
              class RangeVectorImp>
    struct EvaluateCallerTraits
    {
      typedef Traits  BaseTraits;
      typedef typename Traits :: QuadratureType      QuadratureType ;
      typedef typename Traits :: FactorType          FactorType ;
      typedef typename Traits :: LocalDofVectorType  LocalDofVectorType ;
      typedef typename Traits :: Geometry            Geometry ;
 
      typedef BaseFunctionSet   BaseFunctionSetType;
      typedef RangeVectorImp    RangeVectorType;
    };
 
 
    //- base function evaluation interface
    template <class Traits>
    class EvaluateCallerInterface
    {
      typedef EvaluateCallerInterface< Traits >  ThisType;
    public:
      typedef std::unique_ptr< ThisType > StoragePointerType;
 
    protected:
      static const int maxNumBaseFunctions = MAX_NUMBER_OF_BASE_FCT;
      static const int minNumBaseFunctions = MIN_NUMBER_OF_BASE_FCT;
 
      static const int maxQuadNop = MAX_NUMBER_OF_QUAD_POINTS;
      static const int minQuadNop = MIN_NUMBER_OF_QUAD_POINTS;
 
      // maximal number of different quadratures we can use here
      static const int maxQuadratures = 25;
 
      class EvaluatorStorage
      {
        class Item : public std::pair< bool, StoragePointerType >
        {
          typedef std::pair< bool, StoragePointerType > BaseType;
 
        public:
          Item() : BaseType(false, StoragePointerType() ) {}
        };
 
      protected:
        std::vector< Item > storage_;
      public:
        EvaluatorStorage() :
          storage_( maxQuadratures ) {}
 
        Item& get( const size_t id )
        {
          if( id >= storage_.size() )
          {
            storage_.resize( id + 10 );
          }
          return storage_[ id ];
        }
      };
 
 
      EvaluateCallerInterface() {}
 
    public:
      typedef typename Traits :: QuadratureType      QuadratureType ;
      typedef typename Traits :: FactorType          FactorType ;
      typedef typename Traits :: LocalDofVectorType  LocalDofVectorType ;
      typedef typename Traits :: Geometry            Geometry ;
 
      virtual ~EvaluateCallerInterface() {}
 
      virtual void* storageAddress () const = 0;
      virtual size_t storageSize () const = 0;
 
      virtual void axpyRanges( const QuadratureType&,
                               const FactorType& ,
                               LocalDofVectorType & ) const = 0;
 
      virtual void evaluateRanges( const QuadratureType& quad,
                                   const LocalDofVectorType & dofs,
                                   FactorType& factors) const = 0;
 
      virtual void axpyJacobians( const QuadratureType&,
                                  const Geometry&,
                                  const FactorType& ,
                                  LocalDofVectorType & ) const = 0;
 
      virtual void evaluateJacobians( const QuadratureType&,
                                      const Geometry&,
                                      const LocalDofVectorType&,
                                      FactorType&) const = 0;
 
      template < class BaseFunctionSet, class Storage >
      static const StoragePointerType& storage(const BaseFunctionSet& baseSet,
                                               const Storage& dataCache,
                                               const QuadratureType& quad)
      {
        // static vector holding all evaluator instances
        static ThreadSafeValue< EvaluatorStorage > evaluatorStorage;
        EvaluatorStorage& evaluators = *evaluatorStorage;
 
        // check if object already created
        const size_t quadId = quad.id();
        auto& item = evaluators.get( quadId );
        if( ! item.first )
        {
          const int nop = quad.nop();
          static const int dimRange = BaseFunctionSet :: FunctionSpaceType:: dimRange;
          const int numBaseFct = baseSet.size() / dimRange;
 
          typedef EvaluateCallerTraits< Traits, BaseFunctionSet, Storage> NewTraits;
 
#if 0 // NDEBUG
          if( quad.isInterpolationQuadrature( numBaseFct ) )
            std::cout << "EvaluateCallerInterface::storage: Not creating implementation because of interpolation feature!" <<std::endl;
#endif
 
          // if quadrature points or number of basis functions are not within
          // the range of generated code snippets then don't search for
          // a matching combination
          // do not create an evaluation if the quadrature is an interpolation
          // quadrature and matches the number of shape functions
          if( (nop >= minQuadNop && nop <= maxQuadNop) &&
              (numBaseFct >= minNumBaseFunctions && numBaseFct <= maxNumBaseFunctions) &&
              ! quad.isInterpolationQuadrature( numBaseFct ) )
          {
            item.second.reset(
                EvaluateCaller< NewTraits, maxQuadNop, maxNumBaseFunctions >
                   :: create( dataCache , nop, numBaseFct ) );
          }
 
          // if pointer was checked, set flag to true, pointer may still be a nullptr
          item.first = true;
        }
 
        // this can be a nullptr (in this case the default implementation is used)
        return item.second;
      }
    };
 
    template <class Traits,
              int dimRange,
              int quadNop,
              int numBaseFct >
    class EvaluateRealImplementation
      : public EvaluateCallerInterface< typename Traits :: BaseTraits >
    {
    protected:
      typedef typename Traits :: BaseFunctionSetType BaseFunctionSetType;
      typedef typename Traits :: QuadratureType      QuadratureType ;
      typedef typename Traits :: FactorType          FactorType ;
      typedef typename Traits :: LocalDofVectorType  LocalDofVectorType ;
      typedef typename Traits :: Geometry            Geometry ;
      typedef typename Traits :: RangeVectorType     RangeVectorType ;
      typedef typename RangeVectorType :: value_type :: field_type FieldType;
 
      typedef EvaluateRealImplementation< Traits, dimRange, quadNop, numBaseFct > ThisType;
      typedef EvaluateCallerInterface< typename Traits :: BaseTraits >   BaseType;
 
      // A copy is made of the storage here, otherwise strange memory corruptions happen
      // TODO: needs further investigation
      RangeVectorType rangeStorage_; // basis function evaluations
 
      std::vector< std::vector< FieldType > > rangeStorageTransposed_;
      std::vector< std::vector< FieldType > > rangeStorageFlat_;
      mutable std::vector< std::vector< std::vector< FieldType > > > rangeStorageTwisted_;
 
      template <class K>
      int getDim( const DenseVector< K >& vec) const
      {
        return vec.size();
      }
 
      template <class K>
      int getDim( const DenseMatrix< K >& mat) const
      {
        // dimRange == rows which is 1 for basis storage
        return getDim( mat[ 0 ] );
      }
 
      // initialize storage for ranges (i.e. scalar)
      void initRangeStorageTransposed( const std::integral_constant< bool, true > )
      {
        assert( rangeStorage_[ 0 ].size() == 1 );
        {
          const int quadPoints = rangeStorage_.size() / numBaseFct;
          const int faces = quadPoints / quadNop;
          rangeStorageTransposed_.resize( faces );
          for( int f=0; f<faces; ++f )
          {
            auto& rangeStorageTransposed = rangeStorageTransposed_[ f ];
 
            // rearrange such that we store for one basis functions all
            // evaluations for all quadrature points, i.e. numBaseFct * quadNop
            rangeStorageTransposed.resize( numBaseFct * quadNop );
            for( int i=0; i<numBaseFct; ++i )
            {
              const int idx  = i * quadNop;
              for( int j=0; j<quadNop; ++j )
              {
                int qp = f * quadNop + j ;
                assert( j*numBaseFct + i < int(rangeStorage_.size()) );
                // copy and transpose
                rangeStorageTransposed[ idx + j ] = rangeStorage_[ qp*numBaseFct + i ][ 0 ];
              }
            }
          }
        }
      }
 
      // initialize storage for jacobians (i.e. vectors)
      void initRangeStorageTransposed( const std::integral_constant< bool, false > )
      {
        const int dim = rangeStorage_[ 0 ][ 0 ].size();
        {
          const int quadPoints = rangeStorage_.size() / numBaseFct;
          const int faces = quadPoints / quadNop;
          rangeStorageTransposed_.resize( faces );
          rangeStorageFlat_.resize( faces );
          for( int f=0; f<faces; ++f )
          {
            auto& rangeStorageTransposed = rangeStorageTransposed_[ f ];
            auto& rangeStorageFlat = rangeStorageFlat_[ f ];
 
            // rearrange such that we store for one basis functions all
            // evaluations for all quadrature points, i.e. numBaseFct * quadNop
            rangeStorageTransposed.resize( numBaseFct * quadNop * dim );
            rangeStorageFlat.resize( numBaseFct * quadNop * dim );
            for( int i=0; i<numBaseFct; ++i )
            {
              const int idx  = i * (quadNop * dim);
              for( int j=0; j<quadNop; ++j )
              {
                int qp = f * quadNop + j ;
                for( int d=0; d<dim; ++d )
                {
                  rangeStorageFlat[ j*numBaseFct*dim + (i * dim) + d ] = rangeStorage_[ qp*numBaseFct + i ][ 0 ][ d ];
                  rangeStorageTransposed[ idx + (j * dim) + d ] = rangeStorage_[ qp*numBaseFct + i ][ 0 ][ d ];
                }
              }
            }
          }
        }
      }
 
      template <class Quadrature>
      //const DynamicArray< FieldType >&
      const std::vector< FieldType >&
      getTwistedStorage( const Quadrature& quad ) const
      {
        // for evaluation the range storage dimension should be 1 and therefore
        // rangeStorageTransposed should have been filled
        assert( ! rangeStorageTransposed_.empty() );
 
        // if we are in the ranges cases then basis can be stored transposed
        // quadrature points is the outer loop
        if constexpr ( Quadrature::twisted() )
        {
          if( quad.twistId() == 4 ) // no twist
            return rangeStorageTransposed_[ quad.localFaceIndex() ];
 
          auto& rangeStorageTwisted = rangeStorageTwisted_[ quad.twistId() ];
          if( rangeStorageTwisted.empty() )
          {
            const int quadPoints = rangeStorage_.size() / numBaseFct;
            const int faces = quadPoints / quadNop;
            rangeStorageTwisted.resize( faces );
          }
 
          const int f = quad.localFaceIndex() ;
          auto& rangeStorageFace = rangeStorageTwisted[ f ];
          if( rangeStorageFace.empty() )
          {
            // either 1 or dim of grid
            const int dim = getDim( rangeStorage_[ 0 ] );
 
            const auto& rangeStorageTransposed = rangeStorageTransposed_[ f ];
 
            // rearrange such that we store for one basis functions all
            // evaluations for all quadrature points including the twisted mapping
            rangeStorageFace.resize( rangeStorageTransposed.size() );
            for( int i=0; i<numBaseFct; ++i )
            {
              const int idx  = i * (quadNop * dim);
              for( int j=0; j<quadNop; ++j )
              {
                const int qp = quad.localCachingPoint( j );
                for( int d=0; d<dim; ++d )
                {
                  rangeStorageFace[ idx + (j * dim) + d ] = rangeStorageTransposed[ idx + (qp * dim) + d ];
                }
              }
            }
          } // end if( rangeStorageTwisted.empty() )
          return rangeStorageFace;
        }
        else // no twist (i.e. twist = 0 and twistId == 4 (-4 is mapped to 0))
        {
          return rangeStorageTransposed_[ quad.localFaceIndex() ];
        }
      }
    public:
      // type of interface class
      typedef BaseType InterfaceType;
 
      EvaluateRealImplementation( const RangeVectorType& rangeStorage )
        : rangeStorage_( rangeStorage ), rangeStorageTwisted_( 8 ) // 8 different twists
      {
        initRangeStorageTransposed( std::integral_constant< bool,
                                    std::is_same< typename RangeVectorType::value_type,
                                                  Dune::FieldVector< double, 1 > > :: value > () );
      }
 
      virtual void* storageAddress() const { return (void *) &rangeStorage_ ; }
      virtual size_t storageSize() const { return rangeStorage_.size() ; }
 
      virtual void axpyRanges( const QuadratureType& quad,
                               const FactorType& rangeFactors,
                               LocalDofVectorType & dofs ) const
      {
        Codegen::AxpyRanges< BaseFunctionSetType, Geometry, dimRange, quadNop, numBaseFct > :: axpy
          ( quad, rangeStorage_, rangeFactors, dofs );
      }
 
      virtual void evaluateRanges( const QuadratureType& quad,
                                   const LocalDofVectorType & dofs,
                                   FactorType& rangeFactors) const
      {
        Codegen::EvaluateRanges< BaseFunctionSetType, Geometry, dimRange, quadNop, numBaseFct >
          :: eval ( quad, getTwistedStorage( quad ), dofs, rangeFactors );
      }
 
      virtual void axpyJacobians( const QuadratureType& quad,
                                  const Geometry& geometry,
                                  const FactorType& jacFactors,
                                  LocalDofVectorType& dofs) const
      {
        Codegen::AxpyJacobians< BaseFunctionSetType, Geometry, dimRange, quadNop, numBaseFct > :: axpy
          ( quad, geometry, rangeStorageFlat_[ quad.localFaceIndex() ], jacFactors, dofs );
      }
 
      virtual void evaluateJacobians( const QuadratureType& quad,
                                      const Geometry& geometry,
                                      const LocalDofVectorType& dofs,
                                      FactorType& jacFactors) const
      {
        Codegen::EvaluateJacobians< BaseFunctionSetType, Geometry, dimRange, quadNop, numBaseFct > :: eval
          ( quad, geometry, getTwistedStorage( quad ), dofs, jacFactors );
      }
 
      static InterfaceType* create( const RangeVectorType& rangeStorage )
      {
        return new ThisType( rangeStorage );
      }
    };
 
    // The default EvaluateImplementation is empty
    // to create this has to be specified and derived from EvaluateCallerDefault
    template <class Traits,
              int dimRange,
              int quadNop,
              int numBaseFct >
    class EvaluateImplementation
      : public EvaluateCallerInterface< typename Traits :: BaseTraits >
    {
    protected:
      typedef typename Traits :: BaseFunctionSetType BaseFunctionSetType;
      typedef typename Traits :: QuadratureType      QuadratureType ;
      typedef typename Traits :: FactorType          FactorType ;
      typedef typename Traits :: LocalDofVectorType  LocalDofVectorType ;
      typedef typename Traits :: Geometry            Geometry ;
      typedef typename Traits :: RangeVectorType     RangeVectorType ;
 
      typedef EvaluateImplementation< Traits, dimRange, quadNop, numBaseFct > ThisType;
 
      typedef EvaluateCallerInterface< typename Traits :: BaseTraits >   BaseType;
    public:
      // type of interface class
      typedef BaseType InterfaceType;
 
      EvaluateImplementation( const RangeVectorType& rangeStorage )
      {}
 
      virtual void axpyRanges( const QuadratureType& quad,
                               const FactorType& rangeFactors,
                               LocalDofVectorType & dofs ) const
      {
        std::cerr << "ERROR: EvaluateImplementation::axpyRanges not overloaded!" << std::endl;
        std::abort();
      }
 
      virtual void axpyJacobians( const QuadratureType& quad,
                                  const Geometry& geometry,
                                  const FactorType& jacFactors,
                                  LocalDofVectorType& dofs) const
      {
        std::cerr << "ERROR: EvaluateImplementation::axpyJacobians not overloaded!" << std::endl;
        std::abort();
      }
 
      virtual void evaluateRanges( const QuadratureType& quad,
                                   const LocalDofVectorType & dofs,
                                   FactorType& rangeFactors) const
      {
        std::cerr << "ERROR: EvaluateImplementation::evaluateRanges not overloaded!" << std::endl;
        std::abort();
      }
 
      virtual void evaluateJacobians( const QuadratureType& quad,
                                      const Geometry& geometry,
                                      const LocalDofVectorType& dofs,
                                      FactorType& jacFactors) const
      {
        std::cerr << "ERROR: EvaluateImplementation::evaluateJacobians not overloaded!" << std::endl;
        std::abort();
      }
 
      static InterfaceType* create( const RangeVectorType& )
      {
      #ifndef NDEBUG
        std::cout << "Optimized EvaluateImplementation for < dimR="<<dimRange<< ", qp=" << quadNop << ", bases=" << numBaseFct << " > not created, falling back to default!" << std::endl;
        //DUNE_THROW(NotImplemented,"EvaluateImplementation for < " << quadNop << " , " << numBaseFct << " > not created!");
        //return (InterfaceType*) 0;
      #endif
        return nullptr;
      }
    };
 
    template <class Traits,
              int quadNop,
              int numBaseFct >
    class EvaluateCaller
    {
    protected:
      typedef typename Traits :: BaseFunctionSetType BaseFunctionSetType;
      static const int dimRange = BaseFunctionSetType :: FunctionSpaceType:: dimRange;
      typedef typename Traits :: RangeVectorType     RangeVectorType ;
      typedef EvaluateCallerInterface< typename Traits :: BaseTraits >  InterfaceType;
    public:
      static InterfaceType* createObj( const RangeVectorType& rangeStorage,
                                       const size_t numbase )
      {
        if( numBaseFct == numbase )
          return EvaluateImplementation< Traits, dimRange, quadNop, numBaseFct > :: create( rangeStorage );
        else
          return EvaluateCaller< Traits, quadNop, numBaseFct - 1 > :: createObj( rangeStorage, numbase );
      }
 
      static InterfaceType* create( const RangeVectorType& rangeStorage,
                                    const size_t quadnop, const size_t numbase )
      {
        if( quadNop == quadnop )
          return EvaluateCaller< Traits, quadNop, numBaseFct > :: createObj( rangeStorage, numbase );
        else
          return EvaluateCaller< Traits, quadNop - 1, numBaseFct > :: create( rangeStorage, quadnop, numbase );
      }
    };
 
    template <class Traits,
              int numBaseFct >
    class EvaluateCaller< Traits, MIN_NUMBER_OF_QUAD_POINTS, numBaseFct >
    {
    protected:
      typedef typename Traits :: BaseFunctionSetType BaseFunctionSetType;
      static const int dimRange = BaseFunctionSetType :: FunctionSpaceType:: dimRange;
      enum { quadNop = MIN_NUMBER_OF_QUAD_POINTS };
      typedef typename Traits :: RangeVectorType     RangeVectorType ;
      typedef EvaluateCallerInterface< typename Traits :: BaseTraits >  InterfaceType;
    public:
      static InterfaceType* createObj( const RangeVectorType& rangeStorage,
                                       const size_t numbase )
      {
        if( numBaseFct == numbase )
          return EvaluateImplementation< Traits, dimRange, quadNop, numBaseFct > :: create( rangeStorage );
        else
          return EvaluateCaller< Traits, quadNop, numBaseFct - 1 > :: createObj( rangeStorage, numbase );
      }
 
      static InterfaceType* create( const RangeVectorType& rangeStorage,
                                    const size_t quadnop, const size_t numbase )
      {
        if( quadNop == quadnop )
          return EvaluateCaller< Traits, quadNop, numBaseFct > :: createObj( rangeStorage, numbase );
        else
        {
          std::cerr << "ERROR: EvaluateCaller< "<< quadNop << ", " << numBaseFct << " >::createObj: no working combination!" << std::endl;
          std::abort();
        }
      }
    };
 
    template <class Traits,
              int quadNop>
    class EvaluateCaller< Traits, quadNop, MIN_NUMBER_OF_BASE_FCT >
    {
    protected:
      typedef typename Traits :: BaseFunctionSetType BaseFunctionSetType;
      static const int dimRange = BaseFunctionSetType :: FunctionSpaceType:: dimRange;
      enum { numBaseFct = MIN_NUMBER_OF_BASE_FCT };
      typedef typename Traits :: RangeVectorType     RangeVectorType ;
      typedef EvaluateCallerInterface< typename Traits :: BaseTraits >  InterfaceType;
    public:
      static InterfaceType* createObj( const RangeVectorType& rangeStorage,
                                       const size_t numbase )
      {
        if( numBaseFct == numbase )
          return EvaluateImplementation< Traits, dimRange, quadNop, numBaseFct > :: create( rangeStorage );
        else
        {
          std::cerr << "ERROR: EvaluateCaller< "<< quadNop << ", " << numBaseFct << " >::createObj: no working combination!" << std::endl;
          std::abort();
        }
      }
 
      static InterfaceType* create( const RangeVectorType& rangeStorage,
                                    const size_t quadnop, const size_t numbase )
      {
        if( quadNop == quadnop )
          return EvaluateCaller< Traits, quadNop, numBaseFct > :: createObj( rangeStorage, numbase );
        else
        {
          return EvaluateCaller< Traits, quadNop - 1, numBaseFct > :: create( rangeStorage, quadnop, numbase );
        }
      }
    };
 
    template <class Traits>
    class EvaluateCaller< Traits, MIN_NUMBER_OF_QUAD_POINTS, MIN_NUMBER_OF_BASE_FCT>
    {
    protected:
      typedef typename Traits :: BaseFunctionSetType BaseFunctionSetType;
      static const int dimRange = BaseFunctionSetType :: FunctionSpaceType:: dimRange;
      enum { quadNop = MIN_NUMBER_OF_QUAD_POINTS };
      enum { numBaseFct = MIN_NUMBER_OF_BASE_FCT };
      typedef typename Traits :: RangeVectorType     RangeVectorType ;
      typedef EvaluateCallerInterface< typename Traits :: BaseTraits >  InterfaceType;
    public:
      static InterfaceType* createObj( const RangeVectorType& rangeStorage,
                                       const size_t numbase )
      {
        if( numBaseFct == numbase )
          return EvaluateImplementation< Traits, dimRange, quadNop, numBaseFct > :: create( rangeStorage );
        else
        {
          std::cerr << "ERROR: EvaluateCaller< "<< quadNop << ", " << numBaseFct << " >::createObj: no working combination!" << std::endl;
          std::abort();
        }
      }
 
      static InterfaceType* create( const RangeVectorType& rangeStorage,
                                    const size_t quadnop, const size_t numbase )
      {
        if( quadNop == quadnop )
          return EvaluateCaller< Traits, quadNop, numBaseFct > :: createObj( rangeStorage, numbase );
        else
        {
          std::cerr << "ERROR: EvaluateCaller< "<< quadNop << ", " << numBaseFct << " >::create: no working combination!" << std::endl;
          std::abort();
        }
      }
    };
 
    } // namespace Codegen
 
  } // namespace Fem
 
} // namespace Dune
 
#ifdef DUNE_FEM_INCLUDE_AUTOGENERATEDCODE_FILENAME_SPEC
#define CODEGEN_INCLUDEEVALCALLERS
// include specializations of EvaluateImplementation
#include DUNE_FEM_INCLUDE_AUTOGENERATEDCODE_FILENAME_SPEC
#undef CODEGEN_INCLUDEEVALCALLERS
#endif
#endif // #ifndef DUNE_FEM_EVALUATECALLER_HH