functionspace.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2023 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef FUNCTIONSPACE_HH
#define FUNCTIONSPACE_HH
 
#include <cassert>
#include <map>
#include <sstream>
 
#include <boost/mpl/int.hpp>
#include <boost/mpl/range_c.hpp>
#include <boost/multi_array.hpp>
#include <boost/signals2.hpp>
#include <boost/type_traits/remove_pointer.hpp>
#include <boost/type_traits/remove_reference.hpp>
 
#include <boost/fusion/algorithm.hpp>
#include <boost/fusion/sequence.hpp>
 
#include "dune/common/config.h"
// #include "dune/common/dynmatrixev.hh"
#include "dune/common/fmatrixev.hh"
#include "dune/geometry/quadraturerules.hh"
#include "dune/istl/bvector.hh"
 
#include "fem/gridmanager.hh"
#include "fem/pshapefunctions.hh"
#include "fem/shapefunctioncache.hh"
#include "linalg/crsutil.hh"
#include "utilities/detailed_exception.hh"
 
namespace Kaskade
{
 
  // forward declarations
  template <class LocalToGlobalMapper> class FEFunctionSpace;
 
  template <template <class, class> class, class, class> class BoundaryMapper;
 
  namespace FunctionSpace_Detail
  {
    template <class Mapper>
    struct IsBoundaryMapperStruct : std::false_type {};
    template<template <class, class> class DomainMapper, class Scalar, class GridView>
    struct IsBoundaryMapperStruct<BoundaryMapper<DomainMapper,Scalar,GridView>> : std::true_type {};
    // if Mapper is a BoundaryMapper this returns true, otherwise false
    template <class Mapper>
    constexpr bool isBoundaryMapper() {
      return IsBoundaryMapperStruct<Mapper>::value;
    }
 
    template <class Mapper>
    struct ChooseDomainMapperStruct {
      using type = void;
    };
    template<template <class, class> class DomainMapper, class Scalar, class GridView>
    struct ChooseDomainMapperStruct<BoundaryMapper<DomainMapper,Scalar,GridView>> {
      using type = DomainMapper<Scalar,GridView>;
    };
    // if Mapper is a BoundaryMapper this gives the underlying FE mapper, otherwise it gives the void type
    template<class Mapper>
    using ChooseDomainMapper = typename ChooseDomainMapperStruct<Mapper>::type;
 
    // checks whether two LeafGridViews are the same (i.e. whether they belong to the same grid (tested by memory address of grid))
    template<class Grid>
    bool gridViewsAreSame(typename Grid::LeafGridView gv1, typename Grid::LeafGridView gv2) {
      return &(gv1.grid()) == &(gv2.grid());
    }
    // checks whether two LevelGridViews are the same (i.e. whether they belong to the same grid and are of the same level)
    template<class Grid>
    bool gridViewsAreSame(typename Grid::LevelGridView gv1, typename Grid::LevelGridView gv2) {
      return (gv1.template begin<0>()->level() == gv2.template begin<0>()->level()) && (&(gv1.grid()) == &(gv2.grid()));
    }
 
    // --------------------------------------------------------------------------------------------
 
    enum CoefficientTransfer { START, FINISH };
  }
 
  enum ProblemType { VariationalFunctional, WeakFormulation };
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
  template <class Spaces, int Idx>
  struct SpaceType
  {
    using type = std::remove_pointer_t<typename boost::fusion::result_of::value_at_c<Spaces,Idx>::type>;
  };
 
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
  template <class Cell, class ctype>
  std::pair<Cell,Dune::FieldVector<ctype,Cell::dimension>>
  findCellOnLevel(Cell cell, Dune::FieldVector<ctype,Cell::dimension> xi, int level=std::numeric_limits<int>::max())
  {
    ctype const tol = 100*std::numeric_limits<ctype>::epsilon();
    assert(checkInside(cell.type(),xi) < tol);
    assert(level>=0);
 
    using namespace Dune;
 
    if (level <= cell.level()) // go downwards - that's easy
    {
      while (level < cell.level())
      {
        xi = cell.geometryInFather().global(xi);
        cell = cell.father();
      }
    }
    else // go upwards - we have to check in which child the position is 
    {
      // Move up if we have not yet blundered into the leaf or the max level 
      while (cell.level()<level && !cell.isLeaf())
        // Check all the direct children
        for (Cell const& c: descendantElements(cell,cell.level()+1))
        {
          auto xiC = c.geometryInFather().local(xi);
          // and if the point is inside one child, that's it
          if (checkInside(c.type(),xiC) < tol)
          {
            cell = c;
            xi = xiC;
            break;
          }
        }
    } 
    
    return std::make_pair(cell,xi);
  }
  
  
  //---------------------------------------------------------------------
  
    
  template <class GridView>
  Cell<GridView> findCell(GridView const& gv, GlobalPosition<GridView> global)
  {
    const double tol = 1e-12 ;
 
    using namespace Dune;
 
    // Do linear search on coarse grid (level 0).
    auto coarseIterator = std::find_if(gv.grid().levelGridView(0).template begin<0>(), gv.grid().levelGridView(0).template end<0>(),
                                       [&](Cell<GridView> const& cell)
                                       {
                                         return checkInside(cell.type(),cell.geometry().local(global)) < tol;
                                       });
 
    if (coarseIterator == gv.grid().levelGridView(0).template end<0>()) // not found
    {
      std::ostringstream msg;
      msg << "findCell(): global point << " << global << " is not contained in any level 0 cell in the grid!";
      msg.flush();
      throw LookupException(msg.str(),__FILE__,__LINE__);
    }
 
    // Do a hierarchical search for a leaf cell containing the point.
    Cell<GridView> ci(*coarseIterator);
 
    int lastFound = 0;
    
    for( int level = 0; level <= gv.grid().maxLevel(); level++ )
    {
      if (ci.isLeaf())
        return ci;
      
      for (auto hi = ci.hbegin(level+1); hi != ci.hend(level+1) ; ++hi)
        if( checkInside(hi->type(),hi->geometry().local(global)) < tol )
        {
          ci = typename GridView::Grid::template Codim<0>::Entity(*hi);
          lastFound = level;
          break ;
        }
    }
 
    // We should never get here... For a position outside the domain, no coarse grid cell should have been
    // found, leading to an exception being raised above. If, on the other hand, a coarse grid cell has been
    // found, there should be a leaf cell containing the given point as well.
    std::cerr << "findCell() at " << __FILE__ << ':' << __LINE__ << '\n'
        << "Global point << " << global
        << " is not contained in any cell in the grid!\n"
        << " last found on " << lastFound << " of " << gv.grid().maxLevel() << "\n";
    abort();
    return ci; // never get here! (just for formal correctness return something)
  }
 
 
  //---------------------------------------------------------------------
 
 
  template <class FunctionSpace, int m>
  class FunctionSpaceElement
  {
  public:
    using Space = FunctionSpace;
 
    using GridView = typename Space::GridView;
 
    using ctype = typename GridView::ctype;
 
    typedef FunctionSpaceElement <Space, m>  Self;
    typedef typename Space::Scalar Scalar;
    
    using field_type = typename Space::field_type;
    
    typedef Dune::FieldVector<Scalar,m>                 StorageValueType;
    
    using SfValueType = typename Space::SfValueType;
 
    using block_type = StorageValueType;
 
    typedef Dune::BlockVector<StorageValueType>         StorageType;
 
    static int const components = Space::sfComponents * m;
 
    typedef Dune::FieldVector<Scalar,components>        ValueType;
    
    typedef Dune::FieldMatrix<Scalar,components,Space::dimworld> DerivativeType;
 
    using HessianType = Tensor<Scalar,components,Space::dimworld,Space::dimworld>;
 
 
    FunctionSpaceElement(Space const& fs):
      data(fs.degreesOfFreedom()), sp(&fs), transferMe(*this)
    {
      data = 0;
      transferConnection = sp->requestToRefine.connect(transferMe);
    }
 
    FunctionSpaceElement(Self const& fse):
      data(fse.data), sp(fse.sp), transferMe(*this)
    {
      transferConnection = sp->requestToRefine.connect(transferMe);
    }
 
    ~FunctionSpaceElement() 
    {
      transferConnection.disconnect();
      sp = 0;
    }
    
    Self& operator=(Self const& fse) 
    {
      if(this != &fse)     // nothing to do on self-assignment
      {
        if (sp == fse.sp)  // the same function space
          data = fse.data;
        else  // different function spaces
          interpolateGlobally<PlainAverage>(*this,fse);
      }
      return *this;
    }
 
    Self& operator=(Self&& fse) 
    {
      if(this != &fse) { // nothing to do on self-assignment
        if (sp == fse.sp)  // the same function space
        {
          data = std::move(fse.data);
          fse.transferConnection.disconnect();
        }
        else  // different function spaces
          interpolateGlobally<PlainAverage>(*this,fse);
      }
      return *this;
    }
 
    template <class OtherSpace>
    Self& operator=(FunctionSpaceElement<OtherSpace,m> const& f)
    {
      interpolateGlobally<PlainAverage>(*this,f);
      return *this;
    }
 
    template <template <class, class> class DomainMapper>
    Self& operator=(FunctionSpaceElement<FEFunctionSpace<BoundaryMapper<DomainMapper,Scalar,typename Space::GridView>>,m> const& bf)
    {
      bool hasEquivalentUnderlyingMapper = !FunctionSpace_Detail::isBoundaryMapper<typename Space::Mapper>();
      hasEquivalentUnderlyingMapper = (hasEquivalentUnderlyingMapper && std::is_same<DomainMapper<Scalar,typename Space::GridView>, typename Space::Mapper>::value);
      hasEquivalentUnderlyingMapper = (hasEquivalentUnderlyingMapper && space().mapper().maxOrder()==bf.space().mapper().maxOrder());
      hasEquivalentUnderlyingMapper = (hasEquivalentUnderlyingMapper && FunctionSpace_Detail::gridViewsAreSame<typename Space::GridView::Grid>(space().gridView(),bf.space().gridView()));
      if(hasEquivalentUnderlyingMapper) {
        BoundaryMapper<DomainMapper, Scalar, typename Space::GridView> const& boundaryMapper = bf.space().mapper();
        *this = 0.0;
        for(typename StorageType::size_type vi=0;vi<size();++vi) {
          std::pair<bool,typename StorageType::size_type> restrictedIndex = boundaryMapper.unrestrictedToRestrictedIndex(vi);
          if(restrictedIndex.first) {
            data[vi] = bf[restrictedIndex.second];
          }
        }
      } else {
        // TODO: When interpolateGlobally works correctly (better said: as one would expect) for boundary functions as second argument remove the following three lines
        FEFunctionSpace<DomainMapper<Scalar,typename Space::GridView>> domainSpace(bf.space().gridManager(),bf.space().gridView(),bf.space().mapper().maxOrder());
        typename FEFunctionSpace<DomainMapper<Scalar,typename Space::GridView>>::template Element_t<m> df(domainSpace);
        df = bf;
        interpolateGlobally<PlainAverage>(*this,df);
      }
      return *this;
    }
 
    template <class Mapper>
    std::enable_if_t<std::is_same<Mapper,FunctionSpace_Detail::ChooseDomainMapper<typename Space::Mapper>>::value,Self&>
    operator=(FunctionSpaceElement<FEFunctionSpace<Mapper>,m> const& f)
    {
      bool hasEquivalentUnderlyingMapper = FunctionSpace_Detail::isBoundaryMapper<typename Space::Mapper>();
      hasEquivalentUnderlyingMapper = (hasEquivalentUnderlyingMapper && space().mapper().maxOrder()==f.space().mapper().maxOrder());
      hasEquivalentUnderlyingMapper = (hasEquivalentUnderlyingMapper && FunctionSpace_Detail::gridViewsAreSame<typename Space::GridView::Grid>(space().gridView(),f.space().gridView()));
      if (hasEquivalentUnderlyingMapper)
      {
        typename Space::Mapper const& boundaryMapper = space().mapper();
        for(typename StorageType::size_type vi=0;vi<f.size();++vi)
        {
          std::pair<bool,typename StorageType::size_type> restrictedIndex = boundaryMapper.unrestrictedToRestrictedIndex(vi);
          if (restrictedIndex.first)
            data[restrictedIndex.second] = f[vi];
        }
      }
      else
        interpolateGlobally<PlainAverage>(*this,f);
 
      return *this;
    }
 
    template <class Functor>
    Self& operator=(WeakFunctionViewAdaptor<Functor> const& f)
    {
      interpolateGloballyWeak(*this,f);
      return *this;
    }
 
    template <class Space, class Functor>
    Self& operator=(FunctionViewAdaptor<Space,Functor> const& f)
    {
      interpolateGlobally(*this,f);
      return *this;
    }
 
    Self& operator=(StorageType const& v)
    {
      assert(data.size()==v.size());
      data = v;
      return *this;
    }
 
    Self& operator=(StorageValueType const& a)  { for (int i=0; i<data.N(); ++i) data[i] = a; return *this; }
 
    Self& operator=(Scalar val) { for(size_t i=0; i<data.N(); ++i) data[i] = val; return *this; }
    
    Self& operator+=(Self const& f) 
    {
      if (sp==f.sp) // same function space
        coefficients() += f.coefficients();
      else {
        Self tmp(space()); tmp = f; *this += tmp;
      }
      return *this;
    }
 
    template <class OtherSpace>
    Self& operator+=(FunctionSpaceElement<OtherSpace,m> const& f) 
    {
      Self tmp(space()); tmp = f; *this += tmp;
      return *this;
    }
 
    Self& operator+=(StorageType const& v) { assert(data.size()==v.size()); data += v; return *this; }
    
    Self& operator+=(Scalar val)
    {
      for (int i=0; i<data.N(); ++i)
        data[i] += val;
      return *this;
    }
 
    Self& operator-=(Self const& f)
    {
      if (sp==f.sp) // same function space
        coefficients() -= f.coefficients();
      else {
        Self tmp(space()); tmp = f; *this -= tmp;
      }
      return *this;
    }
 
    template <class OtherSpace>
    Self& operator-=(FunctionSpaceElement<OtherSpace,m> const& f) 
    {
      Self tmp(space()); tmp = f; *this -= tmp;
      return *this;
    }
 
    Self& operator-=(StorageType const& v) 
    { 
      assert(data.size()==v.size()); 
      data -= v; 
      return *this; 
    }
 
    Self& axpy(Scalar a, Self const& f)
    {
      if (sp==f.sp) // same function space
        coefficients().axpy(a,f.coefficients());
      else {
        Self tmp(space()); tmp = f; this->axpy(a,tmp);
      }
      return *this;
    }
 
    template <class OtherSpace>
    Self& axpy(Scalar a, FunctionSpaceElement<OtherSpace,m> const& f)
    {
      Self tmp(space()); tmp = f; this->axpy(a,tmp);
      return *this;
    }
 
    Self& axpy(Scalar a, StorageType const& v) { assert(data.size()==v.size()); data.axpy(a,v); return *this; }
 
    Self& operator*=(Scalar a) { data *= a; return *this; }
    
    ValueType value(typename Space::Evaluator const& evaluator) const
    {
      ValueType y(0);
 
      size_t esize = evaluator.size();
 
      for (size_t i=0; i<esize; ++i)
      {
        StorageValueType const& coeff = data[evaluator.globalIndices()[i]];
        auto const& sfVal = evaluator.evalData[i].value;
        for (int n=0; n<m; ++n)
          for (int l=0; l<Space::sfComponents; ++l)
            y[n*Space::sfComponents+l] += coeff[n]*sfVal[l];
      }
 
      return y;
    }
 
    ValueType value(Cell<GridView> const& cell,
                    LocalPosition<GridView> const& local) const
    {
      typename Space::Evaluator eval(space());
      eval.moveTo(cell);
      eval.evaluateAt(local);
 
      return value(eval);
    }
 
    ValueType value(GlobalPosition<GridView> const& x) const
    {
      auto cell = findCell(space().gridView(),x);     // throws LookupException if no cell is found
      auto xi = cell.geometry().local(x);             // may end up slightly outside of reference element
                                                      // due to rounding errors...
      return value(cell,xi);                          // ... but this should have a negligible impact.
    }
 
 
    DerivativeType derivative(typename FunctionSpace::Evaluator const& evaluator) const
    {
      DerivativeType dy(0);
 
      size_t esize = evaluator.size();
      for (size_t i=0; i<esize; ++i)                              // step through all restricted ansatz functions
      {
        auto const& coeff = data[evaluator.globalIndices()[i]];
        auto const& sfGrad = evaluator.evalData[i].derivative;
        for (int n=0; n<m; ++n)                                   // step through all FE function components
          for (int l=0; l<Space::sfComponents; ++l)               // step through all shape function components
            dy[n*Space::sfComponents+l] += coeff[n] * sfGrad[l];
      }
 
      return dy;
    }
    
    DerivativeType derivative(Cell<GridView> const& cell,
                              LocalPosition<GridView> const& local) const
    {
      typename Space::Evaluator eval(space());
      eval.moveTo(cell);
      eval.evaluateAt(local);
 
      return derivative(eval);
    }
 
    DerivativeType derivative(GlobalPosition<GridView> const& global) const
    {
      auto ci = findCell(space().gridView(),global) ;
 
      return derivative(ci,ci.geometry().local(global));
    }
 
 
    HessianType hessian(typename FunctionSpace::Evaluator const& evaluator) const
    {
      HessianType ddy(0);
 
      assert(evaluator.derivatives() >= 2);
 
      size_t esize = evaluator.size();
      for (size_t i=0; i<esize; ++i) {
        StorageValueType const& coeff = data[evaluator.globalIndices()[i]];
        auto const& sfHess = evaluator.evalData[i].hessian;
        for (int n=0; n<m; ++n)
          for (int l=0; l<Space::sfComponents; ++l) {
            auto hess = sfHess[l];
            hess *= coeff[n];
            ddy[n*Space::sfComponents+l] += hess;
          }
      }
 
      return ddy;
    }
 
    HessianType hessian(Cell<GridView> const& cell,
                        LocalPosition<GridView> const& local) const
    {
      typename Space::Evaluator eval(space(),nullptr,2);
      eval.moveTo(cell);
      eval.evaluateAt(local);
 
      return hessian(eval);
    }
 
    HessianType hessian(GlobalPosition<GridView> const& global) const
    {
      auto ci = findCell(space().gridView(),global);
      return hessian(ci,ci.geometry().local(global));
    }
    
    StorageType const& coefficients() const { return data; }
 
    StorageType& coefficients() { return data; }
    
    StorageValueType& operator[](size_t i) { return data[i]; }
 
    StorageValueType const& operator[](size_t i) const { return data[i]; }
    
    int dim() const { return data.dim(); }
 
    int size() const { assert(data.size()*m==dim()); return data.size(); }
 
    Space const& space() const { return *sp; }
 
    int order(Cell<GridView> const& cell) const
    {
      return space().mapper().shapefunctions(cell).order();
    }
 
    // Why not simply calling sfs.order directly? Can we delete this method?
    template<class SFS>
    [[deprecated]] int order(SFS const& sfs) const
    {
      return sfs.order();
    }
 
    void transfer(typename TransferData<Space>::TransferMatrix const& transferMatrix)
    {
      std::unique_ptr<Dune::BlockVector<StorageValueType>> newCoeff = transferMatrix.apply(data);
      data = *newCoeff;
    }
 
  private:
 
    struct TransferMe
    {
      TransferMe(Self& fe) : myFSE(fe) {}
 
      TransferMe(TransferMe const& other)
      : myFSE(other.myFSE)
      {}
 
      void operator()(typename TransferData<Space>::TransferMatrix const* transfer)
      {
        if (transfer)
        {
          // auto f = [&]() { myFSE.transfer(*transfer); };
          // f();
          future = std::async(std::launch::async,[transfer,f=&myFSE]() { f->transfer(*transfer); });
        }
        else
        {
          assert(future.valid());
          future.get();
        }
      }
 
    private:
      Self& myFSE;
      std::future<void> future;
    };
 
 
 
  protected:
    Dune::BlockVector<StorageValueType> data;
    Space const* sp;
 
  private:
    TransferMe transferMe;
    boost::signals2::connection transferConnection;
  };
  
  template <class Space, int m, class OutIter>
  OutIter vectorToSequence(FunctionSpaceElement<Space,m> const& v, OutIter iter)
  {
    return vectorToSequence(v.coefficients(),iter);
  }
  
  template <class Space, int m, class InIter>
  InIter vectorFromSequence(FunctionSpaceElement<Space,m>& v, InIter in) 
  {
    return vectorFromSequence(v.coefficients(),in);
  }
 
  namespace FunctionSpace_Detail
  {
    // if an FEFunction is restricted to the boundary transform it in to a new FEFunction over the whole domain,
    // otherwise leave it unchanged. Used by writeVTK.
    template <class FEFunction>
    class ToDomainRepresentation {
    public:
      ToDomainRepresentation(FEFunction const& f_) : f(f_){}
 
      FEFunction const& get() const {
        return f;
      }
 
    private:
      FEFunction const& f;
    };
 
    template <template <class, class> class DomainMapper, typename Scalar, typename GridView, int m>
    class ToDomainRepresentation<FunctionSpaceElement<FEFunctionSpace<BoundaryMapper<DomainMapper, Scalar, GridView>>, m>>
    {
 
      using BoundaryFEFunction = FunctionSpaceElement<FEFunctionSpace<BoundaryMapper<DomainMapper, Scalar, GridView>>, m>;
      using DomainFEFunction = FunctionSpaceElement<FEFunctionSpace<DomainMapper<Scalar, GridView>>, m>;
      using DomainFESpace = typename DomainFEFunction::Space;
 
    public:
      ToDomainRepresentation(BoundaryFEFunction const& f_) : space(f_.space().gridManager(),f_.space().gridView(),f_.space().mapper().maxOrder()), f(space) {
        f = f_;
      }
 
      DomainFEFunction const& get() const {
        return f;
      }
 
    private:
      DomainFESpace space;
      DomainFEFunction f;
    };
  }
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
  template <class LocalToGlobalMapper>
  class FEFunctionSpace
  {
    typedef FEFunctionSpace<LocalToGlobalMapper> Self;
 
  public:
    using GridView = typename LocalToGlobalMapper::GridView;
    typedef typename GridView::Grid                    Grid;
    typedef  typename GridView::IndexSet               IndexSet;
 
    typedef typename LocalToGlobalMapper::Scalar    Scalar;
    
    using field_type = Scalar;
 
    typedef typename Grid::ctype ctype;
    
    typedef LocalToGlobalMapper Mapper;
 
    static int const dim = Grid::dimension;
 
    static int const dimworld = Grid::dimensionworld;
 
    static int const sfComponents = Mapper::ShapeFunctionSet::comps;
 
    static int const continuity = Mapper::continuity;
 
    template <int m>
    struct [[deprecated]] Element
    {
      typedef FunctionSpaceElement<Self,m> type;
    };
    
    template <int m>
    using Element_t = FunctionSpaceElement<Self,m>;
 
    struct Evaluator
    {
      typedef typename IndexSet::IndexType IndexType;
 
      typedef Self Space;
      
      using Cell = Kaskade::Cell<GridView>;
 
      static int const sfComponents = Mapper::ShapeFunctionSet::comps;
 
      typedef VariationalArg<Scalar,dimworld,sfComponents> VarArg;
 
      typedef typename Mapper::GlobalIndexRange GlobalIndexRange;
    
      typedef typename Mapper::SortedIndexRange SortedIndexRange;
 
      Evaluator(Self const& space_, ShapeFunctionCache<Grid,Scalar>* cache_=nullptr, int deriv_ = 1) :
        space(space_), sf(nullptr), cache(cache_), sfValues(nullptr), combiner_(nullptr),
        globIdx(space.mapper().initGlobalIndexRange()), sortedIdx(space.mapper().initSortedIndexRange()),
        deriv(deriv_), currentCell(nullptr)
      { }
 
      Evaluator(Evaluator const& eval):
        evalData(eval.evalData),
        space(eval.space), sf(eval.sf), sfSize(eval.sfSize), cache(eval.cache), sfValues(eval.sfValues), lastQr(eval.lastQr), lastSubId(eval.lastSubId), combiner_(nullptr),
        globIdx(eval.globIdx), sortedIdx(eval.sortedIdx), deriv(eval.deriv), currentXloc(eval.currentXloc), currentCell(nullptr)
      {
        if (eval.combiner_) 
          combiner_ = std::make_unique<typename Mapper::Combiner>(*eval.combiner_);
        if (eval.currentCell) 
          currentCell = std::make_unique<Cell>(*eval.currentCell);
      }
 
      Evaluator& operator=(Evaluator const& eval)
      {
        assert(&space==&eval.space);
        sf = eval.sf;
        sfSize = eval.sfSize();
        cache = eval.cache;
        sfValues = eval.sfValues;
        lastQr = nullptr;
        lastSubId = -1;
        converter = eval.converter;
        if (eval.combiner_) {
          combiner_ = std::make_unique<typename Mapper::Combiner>(*eval.combiner_);
        } else {
          combiner_.reset();
        }
        globIdx = eval.globIdx;
        sortedIdx = eval.sortedIdx;
        evalData = eval.evalData;
        deriv = eval.deriv;
        currentXloc = eval.currentXloc;
        if (eval.currentCell) {
          currentCell = std::make_unique<Cell>(*eval.currentCell);
        } else {
          currentCell.reset();
        }
      
        return *this;
      }
 
      void moveTo(Cell const& cell, IndexType const index)
      {
        sf = &space.mapper().shapefunctions(index);
        sfSize = sf->size();
        assert(sfSize>=0);
 
        currentCell = std::make_unique<Cell>(cell);
 
        converter.moveTo(*currentCell);
 
        combiner_ = std::make_unique<typename Mapper::Combiner>(space.mapper().combiner(*currentCell,index));
 
        globIdx = space.mapper().globalIndices(index);
        sortedIdx = space.mapper().sortedIndices(index);
      }
      
      void moveTo(Cell const& cell)
      {
        moveTo(cell,space.indexSet().index(cell));
      }
    
      template <int subDim>
      void useQuadratureRule(Dune::QuadratureRule<typename Grid::ctype,subDim> const& qr, int subId)
      {
        // reusing data does not work for boundary mapper in its current (non optimal) implementation,
        // where every cell has its own restricted shape function set.
        // TODO: When the implementation of boundary mapper is improved (i.e. all cells use the
        // same shape function set and restriction is handled by the Combiner) remove the following line.
        if(FunctionSpace_Detail::isBoundaryMapper<Mapper>())
          return;
        // check whether we have the data corresponding to the lookup key (qr,subId) already available
        if (sfValues && lastQr==&qr && lastSubId==subId) 
          return;
        
        if (cache)
        {
          sfValues = &cache->evaluate(*sf,qr,subId); // extract data from cache
          lastQr = &qr;                              // remember extraction key
          lastSubId = subId;
        }
      }
 
      void evaluateAt(Dune::FieldVector<ctype,dim> const& xi)
      {
        currentXloc = xi;
        evalData.clear();
 
        converter.setLocalPosition(xi);
 
        evalData.reserve(sfSize);
        for (int i=0; i<sfSize; ++i) 
        {
          typedef VariationalArg<Scalar,dim,sfComponents> SfArg;
          
          typename Mapper::ShapeFunctionSet::value_type const& sfun = (*sf)[i];
          evalData.push_back( derivatives()<=1?
                                converter.global(SfArg(sfun.evaluateFunction(xi),sfun.evaluateDerivative(xi)),1) :
                                converter.global(SfArg(sfun.evaluateFunction(xi),sfun.evaluateDerivative(xi),sfun.evaluate2ndDerivative(xi)),2) );
        }
        combiner_->rightTransform(evalData);
      }
 
      template <int subDim>
      void evaluateAt(Dune::FieldVector<ctype,dim> const& x,
                      Dune::QuadratureRule<typename Grid::ctype,subDim> const& qr, int ip, int subId)
      {
        if (cache) 
        {
          currentXloc = x;
          evalData.clear();
          converter.setLocalPosition(x);
          
          if (sfValues) 
            for (int i=0; i<sfSize; ++i)
              evalData.push_back( converter.global((*sfValues)[ip][i],derivatives()) );
          else
          {
            auto const& localData = cache->evaluate(*sf,qr,ip,subId);
            
            for (int i=0; i<sfSize; ++i)
              evalData.push_back( converter.global(localData[i],derivatives()) );
          }
 
          combiner_->rightTransform(evalData);
        } 
        else 
          evaluateAt(x);
      }
 
      GlobalIndexRange globalIndices() const { return globIdx; }
    
      SortedIndexRange sortedIndices() const { return sortedIdx; }
 
      int size() const 
      { 
        return globalIndices().size(); 
      }
 
 
      int order() const { return sf->order(); }
 
      typename Mapper::ShapeFunctionSet const& shapeFunctions() const { return *sf; }
 
      typename Mapper::Combiner const& combiner() const { return *combiner_; }
 
      std::vector<VarArg> evalData;
 
      GridView const& gridView() const
      {
        return space.gridView();
      }
 
      LocalToGlobalMapper const& mapper() const
      {
        return space.mapper();
      }
      
      Cell const& cell() const
      {
        assert(currentCell);
        return *currentCell;
      }
      
      Dune::FieldVector<ctype,dim> const& xloc() const { return currentXloc; }
      
      int derivatives() const { return deriv; }
 
    private:
      Self const&                                    space;
      typename Mapper::ShapeFunctionSet const*       sf;
      int                                            sfSize;
      ShapeFunctionCache<Grid,Scalar>*               cache;
      typename ShapeFunctionCache<Grid,Scalar>::template DataType<sfComponents>::type const* sfValues;
      void const*                                    lastQr;    // lookup key for cached shape function values
      int                                            lastSubId; // lookup key for cached shape function values
      typename Mapper::Converter                     converter;
      std::unique_ptr<typename Mapper::Combiner>     combiner_; // TODO: why held by pointer? Wouldn't a member variable make more sense and simplify things? (For this, something like a default constructor would be needed)
      GlobalIndexRange                               globIdx;
      SortedIndexRange                               sortedIdx;
      int                                            deriv;         // number of derivatives to evaluate
      Dune::FieldVector<ctype,dim>                   currentXloc;   // local coordinates of current evaluation point
      std::unique_ptr<Cell>                          currentCell;   // cell containing current evaluation point (want Evaluator to store its own Entity object
                                                                    // (previously it only stored a raw pointer),
                                                                    // unfortunately not all grid implementations (especially AluGrid) provide a 
                                                                    // default constructor for Entity, so we cannot hold it by value.
    }; // end of class Evaluator
 
    using SfValueType = typename Evaluator::VarArg::ValueType;
 
    template <typename... Args>
    FEFunctionSpace(GridManagerBase<Grid>& gridMan,
                    GridView const& gridView_,
                    Args... args)
    : gridman(gridMan)
    , m_gridView(gridView_)
    , mapper_(m_gridView,args...)
    , rtf(*this,gridman)
    {
      refConnection = gridman.signals.informAboutRefinement.connect(int(GridSignals::allocResources),rtf);
    }
 
    FEFunctionSpace(GridManagerBase<Grid>& gridman_, Mapper&& m)
    : gridman(gridman_)
    , m_gridView(m.gridView())
    , mapper_(m)
    , rtf(*this,gridman)
    {
      refConnection = gridman.signals.informAboutRefinement.connect(int(GridSignals::allocResources),rtf);
    }
 
    FEFunctionSpace(FEFunctionSpace<LocalToGlobalMapper> const& fs) :
      gridman(fs.gridman),
      m_gridView(fs.m_gridView),
      mapper_(fs.mapper()),
      rtf(*this,gridman)
    {
      refConnection = gridman.signals.informAboutRefinement.connect(int(GridSignals::allocResources),rtf);
    }
 
    ~FEFunctionSpace() 
    {
      refConnection.disconnect();
    }
 
    size_t degreesOfFreedom() const { return mapper_.size(); }
 
 
    Grid const& grid() const { return gridman.grid(); }
 
    GridManagerBase<Grid>&  gridManager() const { return gridman; }
 
    IndexSet const& indexSet() const { return m_gridView.indexSet(); }
 
    GridView const& gridView() const { return m_gridView; }
 
    LocalToGlobalMapper const & mapper() const { return mapper_; }
    
    LocalToGlobalMapper /* const */ & mapper() /* const */ { return mapper_; }
 
    Evaluator evaluator(ShapeFunctionCache<Grid,Scalar>* cache=nullptr, int deriv = 1) const { return Evaluator(*this,cache,deriv); }
    
    template <int m>
    Element_t<m> element() const
    {
      return Element_t<m>(*this);
    }
 
    std::vector<std::pair<size_t,SfValueType>> linearCombination(Evaluator const& evaluator) const
    {
      std::vector<std::pair<size_t,SfValueType>> ic;
 
      size_t esize = evaluator.size();
      ic.reserve(esize);
 
      for (size_t j=0; j<esize; ++j)
      {
        size_t i = evaluator.globalIndices()[j];
        auto const& c = evaluator.evalData[j].value;
        ic.push_back(std::make_pair(i,c));
      }
 
      return ic;
    }
 
    std::vector<std::pair<size_t,SfValueType>>
    linearCombination(Cell<GridView> const& cell, LocalPosition<GridView> const& xi) const
    {
      Evaluator eval(*this);
      eval.moveTo(cell);
      eval.evaluateAt(xi);
 
      return linearCombination(eval);
    }
 
    std::vector<std::pair<size_t,SfValueType>>
    linearCombination(GlobalPosition<GridView> const& x) const
    {
      auto cell = findCell(gridView(),x);             // throws LookupException if no cell is found
      auto xi = cell.geometry().local(x);             // may end up slightly outside of reference element
 
      return linearCombination(cell,xi);
    }
 
  private:
 
    // The slot class that connects to the grid manager signal to be informed about the 
    // state of the refinement process.
    struct ReactToRefinement
    {
      ReactToRefinement(Self & space, GridManagerBase<Grid>& gm_) 
      : mySpace(space), gm(gm_)  {};
 
      ReactToRefinement(ReactToRefinement const& rtr)
      : mySpace(rtr.mySpace)
      , gm(rtr.gm)
      {}
 
      void operator()(typename GridSignals::Status const ref)
      {
        switch (ref)
        {
          case GridSignals::BeforeRefinement:
            assert(!future.valid());
            if (!mySpace.requestToRefine.empty())               // no functions registered -> nothing to do
              future = std::async(std::launch::async,           // otherwise create transfer matrix data
                                  [&]() { transferData.reset(new TransferData<Self>(mySpace)); });
            break;
 
          case GridSignals::BeforeRefinementComplete:           // Now we need to have the data ready
            if (future.valid())                                 // (if there is any).
              future.get();                                     // Block execution until the data is complete.
            break;
 
          case GridSignals::AfterRefinement:
            mySpace.mapper().update();                          // new grid available -> new indexing
            if(!mySpace.requestToRefine.empty())                // no functions registered -> nothing to do
            {
              assert(transferData!=0);
              assert(!future.valid());
              future = std::async(std::launch::async,
                                  [&]() {
                                    transferMatrix = transferData->transferMatrix();
                                    transferData.reset();
                                    mySpace.requestToRefine(transferMatrix.get());
                                  });
            }
            break;
 
          case GridSignals::TransferCompleted:
            if (future.valid())
              future.get();
            mySpace.requestToRefine(nullptr);
            transferMatrix.reset();
        }
      }
      
    private:
      Self&                               mySpace;
      GridManagerBase<Grid>&              gm;
      std::unique_ptr<TransferData<Self>> transferData;
      std::unique_ptr<typename TransferData<Self>::TransferMatrix> transferMatrix;
      std::future<void>                   future;
    };
 
    GridManagerBase<Grid>& gridman;
    GridView               m_gridView;
    LocalToGlobalMapper    mapper_;
    ReactToRefinement      rtf;
    boost::signals2::connection refConnection;
 
  public:
    mutable typename boost::signals2::signal<void (typename TransferData<Self>::TransferMatrix const*) > requestToRefine;
 
  };
 
 
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
  
  template <class Cache>
  class GetEvaluators
  {
  public:
    GetEvaluators(Cache* cache_, std::map<void const*,int> const& deriv_ = std::map<void const*,int>()):
      cache(cache_), deriv(deriv_)
    {}
 
    template <class Space>
    typename Space::Evaluator operator()(Space const* space) const
    {
      auto it = deriv.find(space);
      int derivatives = it==deriv.end()? 1: it->second; // defaults to first derivatives
      return typename Space::Evaluator(*space,cache,derivatives);
    }
 
  private:
    Cache* cache;
    std::map<void const*,int> const& deriv;
  };
 
  class GetEvaluatorTypes
  {
  public:
    template <class Space>
    typename Space::Evaluator operator()(Space const*);
  };
 
  template <class Spaces, class ShapeFunctionCache>
  auto getEvaluators(Spaces const& spaces, ShapeFunctionCache* cache,
                     std::map<void const*,int> const& deriv = std::map<void const*,int>())
  {
//     return boost::fusion::as_vector(boost::fusion::transform(spaces,GetEvaluators<ShapeFunctionCache>(cache,deriv)));
    return boost::fusion::as_vector(boost::fusion::transform(spaces,[&](auto const* space) 
    { 
      auto it = deriv.find(space);
      int derivatives = it==deriv.end()? 1: it->second; // defaults to first derivatives
      return space->evaluator(cache,derivatives); 
    }));
  }
  
  template <class Spaces>
  using Evaluators = typename boost::fusion::result_of::as_vector<
                          typename boost::fusion::result_of::transform<Spaces,GetEvaluatorTypes>::type
                       >::type;
 
  //---------------------------------------------------------------------
  template<class Evaluators>
  class [[deprecated]] GetVarArgs
  {
  public:
    template <class T> struct result {};
 
    template <class VD>
    struct result<GetVarArgs(VD)>
    {
      static int const sidx = VD::spaceIndex;
      typedef typename boost::fusion::result_of::value_at_c<Evaluators, sidx>::type::VarArg type;
    };
  };
 
  //---------------------------------------------------------------------
  
  template <class Evaluators, class Cell>
  void moveEvaluatorsToCell(Evaluators& evals, Cell const& cell)
  {
    boost::fusion::for_each(evals,[&cell](auto& eval) { eval.moveTo(cell); });
  }
 
  template <class Evaluators, class Cell, class Index>
  void moveEvaluatorsToCell(Evaluators& evals, Cell const& cell, Index idx)
  {
    boost::fusion::for_each(evals,[&cell,idx](auto& eval) { eval.moveTo(cell,idx); });
  }
  
//---------------------------------------------------------------------
 
  template <class Evaluator, class QuadratureRule>
  void useQuadratureRuleInEvaluators(Evaluator& evals, QuadratureRule const& qr, int subId)
  {
    boost::fusion::for_each(evals,[&qr,subId](auto& eval) { eval.useQuadratureRule(qr,subId); });
  }
 
  //---------------------------------------------------------------------
 
  template <class Evaluators, class CoordType, int dim, int subDim>
  void moveEvaluatorsToIntegrationPoint(Evaluators& evals, Dune::FieldVector<CoordType,dim> const& x,
                                        Dune::QuadratureRule<CoordType,subDim> const& qr, int ip, int subId)
  {
    boost::fusion::for_each(evals,[&x,&qr,ip,subId](auto& eval) { eval.evaluateAt(x,qr,ip,subId); });
  }
 
  template <class Evaluators, class CoordType, int dim>
  void moveEvaluatorsToIntegrationPoint(Evaluators& evals, Dune::FieldVector<CoordType,dim> const& x)
  {
    boost::fusion::for_each(evals,[&x](auto& eval) { eval.evaluateAt(x); });
  }
 
  //---------------------------------------------------------------------
  
  template <class Evaluators>
  int maxOrder(Evaluators const& evals)
  {
    auto getOrder = [](auto const& eval) -> int { return eval.order(); };
    auto max = [](auto a, auto b) -> int { return std::max(a,b); };
    return boost::fusion::fold(boost::fusion::transform(evals,getOrder),0,max);
  }
 
} /* end of namespace Kaskade */
 
#endif