fetransfer.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 FETRANSFER_HH
#define FETRANSFER_HH
 
 
#include <limits>
#include <memory>       // std::unique_ptr
#include <vector>
#include <map>
#include <set>
#include <stack>
#include <iterator>
#include <utility>      // std::move
 
 
#include "dune/grid/common/grid.hh"
#include "dune/istl/bcrsmatrix.hh"
#include "dune/common/fvector.hh"
#include "dune/common/fmatrix.hh"
#include "dune/grid/io/file/dgfparser/dgfparser.hh"
//#include "dune/grid/geometrygrid.hh" // why should this be needed?
 
#include "fem/fixdune.hh"
#include "fem/gridmanager.hh"
#include "fem/mllgeometry.hh"
#include "linalg/dynamicMatrix.hh"
 
namespace Kaskade
{
  template<class G, class T> class CellData;
  template <class Space, class CP> class TransferData;
  struct AdaptationCoarseningPolicy ;
 
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
  template <class Space, class ShapeFunctionSet>
  void evaluateGlobalShapeFunctions(Space const& space,
                                    Cell<typename Space::GridView> const& cell,
                                    std::vector<LocalPosition<typename Space::GridView>> const& xi,
                                    DynamicMatrix<Dune::FieldMatrix<typename Space::Scalar,Space::sfComponents,1>>& sfValues,
                                    ShapeFunctionSet const& shapeFunctionSet)
  {
 
    sfValues.setSize(xi.size(),shapeFunctionSet.size());
 
    // Evaluation and transformation of shape function values. This
    // realizes \Phi.
    shapeFunctionSet.evaluate(xi,sfValues);
 
    // Convert the local to global values. This realizes \Psi
    typename Space::Mapper::Converter psi(cell);
    for (int i=0; i<sfValues.N(); ++i) 
    {
      psi.setLocalPosition(xi[i]);
      for (int j=0; j<sfValues.M(); ++j)
        sfValues[i][j] = psi.global(sfValues[i][j]);
    }
  }
 
  template <class Space>
  void evaluateGlobalShapeFunctions(Space const& space,
                                    Cell<typename Space::GridView> const& cell,
                                    std::vector<LocalPosition<typename Space::GridView>> const& x,
                                    DynamicMatrix<Dune::FieldMatrix<typename Space::Scalar,Space::sfComponents,1> >& sfValues)
  {
    evaluateGlobalShapeFunctions(space,cell,x,sfValues,space.mapper().shapefunctions(cell));
  }
  
  template <class Space, class ShapeFunctionSet>
  void approximateGlobalValues(Space const& space,
                               typename Space::Grid::template Codim<0>::Entity const& cell,
                               DynamicMatrix<Dune::FieldMatrix<typename Space::Scalar,Space::sfComponents,1>>& globalValues,
                               DynamicMatrix<Dune::FieldMatrix<typename Space::Scalar,1,1>>& coeff,
                               ShapeFunctionSet const& sfs)
  {
    // Transform global values to local values. This is \Psi^{-1}.
    typename Space::Mapper::Converter psi(cell);
    
    auto const& iNodes = sfs.interpolationNodes();
 
    for (int i=0; i<globalValues.N(); ++i) 
    {
      psi.setLocalPosition(iNodes[i]);
      for (int j=0; j<globalValues.M(); ++j)
        globalValues[i][j] = psi.local(globalValues[i][j]);
    }
    // Interpolate local values, giving shape function coefficients. This is \Phi^+.
    sfs.interpolate(globalValues,coeff);
  }
 
  template <class Space>
  void approximateGlobalValues(Space const& space,
                               typename Space::Grid::template Codim<0>::Entity const& cell,
                               DynamicMatrix<Dune::FieldMatrix<typename Space::Scalar,Space::sfComponents,1>>& globalValues,
                               DynamicMatrix<Dune::FieldMatrix<typename Space::Scalar,1,1>>& coeff)
  {
    approximateGlobalValues(space,cell,globalValues,coeff,space.mapper().shapefunctions(cell));
  }
 
  template <class ChildSpace, class FatherSpace>
  void localProlongationMatrix(ChildSpace const& childSpace,
                               Cell<typename ChildSpace::Grid> const& child,
                               FatherSpace const& fatherSpace,
                               Cell<typename FatherSpace::Grid> const& father,
                               DynamicMatrix<Dune::FieldMatrix<typename ChildSpace::Scalar,1,1>>& prolongation,
                               typename ChildSpace::Mapper::ShapeFunctionSet const& childSfs,
                               typename FatherSpace::Mapper::ShapeFunctionSet const& fatherSfs)
  {
    using GridView = typename ChildSpace::GridView;
    typedef typename ChildSpace::Grid Grid;
    typedef typename ChildSpace::Scalar   Scalar;
    int const sfComponents = ChildSpace::sfComponents;
 
    // Obtain the child interpolation nodes and map them to the father.
    std::vector<LocalPosition<GridView>> iNodes(childSfs.interpolationNodes());
 
    std::vector<int> dummy(iNodes.size());
    MultiLevelLocalGeometry<Grid> const mlGeometry(childSpace.grid(),child,father,MultiLevelLocalGeometry<Grid>::ChildIsGlobal);
    mlGeometry.local(iNodes,dummy);
    assert(dummy.size()==iNodes.size()); // \todo: we assume the child is contained in the father. Is this guaranteed?
 
    // Evaluate global shape functions on father cell.
    DynamicMatrix<Dune::FieldMatrix<Scalar,sfComponents,1>> afValues;
    evaluateGlobalShapeFunctions(fatherSpace,father,iNodes,afValues,fatherSfs);
 
    // Interpolate global values on the child cell.
    approximateGlobalValues(childSpace,child,afValues,prolongation,childSfs);
  }
 
  template <class ChildSpace, class FatherSpace>
  void localProlongationMatrix(ChildSpace const& childSpace,
                               Cell<typename ChildSpace::GridView> const& child,
                               FatherSpace const& fatherSpace,
                               Cell<typename FatherSpace::GridView> const& father,
                               DynamicMatrix<Dune::FieldMatrix<typename ChildSpace::Scalar,1,1> >& prolongation)
  {
    localProlongationMatrix(childSpace,child,fatherSpace,father,prolongation,
                            childSpace.mapper().shapefunctions(child),fatherSpace.mapper().shapefunctions(father));
  }
  
  struct AdaptationCoarseningPolicy
  {
    void update (int) { }
 
    template <class Cell>
    bool accept(Cell const& cell) const 
    { 
      return cell.isLeaf(); 
    }
 
    template <class Cell>
    bool operator()(Cell const& cell) const 
    { 
      return cell.mightBeCoarsened(); 
    }
  };
 
  struct MultilevelCoarseningPolicy
  {
    MultilevelCoarseningPolicy(int level_):
      level(level_)
    {}
 
    void update( int level_ ) { level = level_ ; }
 
    template <class Cell>
    bool accept(Cell const& cell) const 
    { 
      return cell.level()==level; 
    }
 
    template <class Cell>
    bool operator()(Cell const& cell) const 
    { 
      return false; 
    }
 
  private:
    int level;
  };
 
  template <class Space, class CoarseningPolicy=AdaptationCoarseningPolicy>
  class LocalTransfer
  {
  public:
    typedef typename Space::Grid            Grid;
    typedef typename Space::IndexSet        IndexSet;
    typedef typename Space::Scalar          Scalar;
    typedef typename Grid::GlobalIdSet::IdType Id;
    typedef Dune::FieldVector<Scalar,Space::sfComponents>   SfValue;
    using Cell = Kaskade::Cell<Grid>;
 
    typedef DynamicMatrix<Dune::FieldMatrix<Scalar,1,1>>    LocalTransferMatrix;
 
 
  private:
    using ShapeFunctionSet = typename Space::Mapper::ShapeFunctionSet;
 
  public:
    LocalTransfer() = default;
 
    LocalTransfer(LocalTransfer&& other)
    : globalIndices(std::forward<std::vector<size_t>>(other.globalIndices))
    , restrictionMatrix(std::forward<LocalTransferMatrix>(other.restrictionMatrix))
    , shapeFunctionSet(other.shapeFunctionSet)
    {    }
 
    LocalTransfer<Space,CoarseningPolicy>&
    operator =(LocalTransfer&& other)
    {
      globalIndices = std::forward<std::vector<size_t>>(other.globalIndices);
      restrictionMatrix = std::forward<LocalTransferMatrix>(other.restrictionMatrix);
      shapeFunctionSet = other.shapeFunctionSet;
      return *this;
    }
 
    LocalTransfer(LocalTransfer const& other) = default;
 
    LocalTransfer(Space const& space,
                  Cell father,
                  CoarseningPolicy const& relevancePolicy = CoarseningPolicy())
    {
      // TODO: if father.isLeaf() is true, the local transfer matrix is essentially the identity
      //       (ok, it's the combiner). Check and implement a shortcut. After all, at least half of
      //       the cells visited here are leafs, and extracting the identity is tremendously cheap.
      // TODO: Higher performance solution (for special case, though) has been developed in coarsening.hh.
      //       See whether this can be used to improve the performance here.
      Grid const& grid = space.grid();
      assert(&grid == &space.grid());
 
      // Collect all relevant children descending from the father cell. "Relevant" meaning that 
      // projection to the father shall be done from that child cell. For normal grid transfer, 
      // these are the leaf cells, for multilevel projections, these are children on a certain 
      // level.
      std::vector<Cell> children;
      if ( relevancePolicy.accept(father) )
        children.push_back(father);
      else
      {
        auto end = father.hend(grid.maxLevel());
        for (auto hi=father.hbegin(grid.maxLevel()); hi!=end; ++hi)
          if (relevancePolicy.accept(*hi))
            children.push_back(Cell(*hi));
      }
 
      // Get the shape function set of the father cell.
      auto const& fatherShapeFunctions = space.mapper().shapefunctions(father);
 
      // Compute the global indices associated to the relevant child cells 
      // (possibly including the father itself).
      for (Cell const& child: children)
      {
        auto gi = space.mapper().globalIndices(child);
        globalIndices.insert(globalIndices.end(),gi.begin(),gi.end());
      }
      std::sort(globalIndices.begin(),globalIndices.end());
      globalIndices.erase(std::unique(globalIndices.begin(),globalIndices.end()),globalIndices.end());
 
      // Initialize local restriction data.
      // The restriction data stores for each interpolation point of the father cell (corresponding
      // to the rows of the restrictionData matrix) the value of all child's global ansatz functions
      // (corresponding to the rows of the restrictionData matrix) - considering only those ansatz
      // functions with a support intersecting the father cell.
      typedef typename Space::Mapper::ShapeFunctionSet::SfValueArray SfValueArray;
      SfValueArray restrictionData;
      restrictionData.setSize(fatherShapeFunctions.interpolationNodes().size(),globalIndices.size());
      restrictionData = 0;
 
      // Compute local restriction matrices.
      // First declare variables that allocate memory in front of the loop to prevent frequent reallocations.
      std::vector<int> fCount(restrictionData.N(),0);
      // LocalTransferMatrix p;
      SfValueArray afValues; 
      std::vector<size_t> globIdx;
      using Position = LocalPosition<Grid>;
      std::vector<Position> iNodes;
      std::vector<int> nodesInside;
      
      for (Cell const& child: children)
      {
        // Compute the mapping "globIdx" of degrees of freedom on the child to 
        // the set of dofs on all relevant children.
        auto gi = space.mapper().globalIndices(child);
        globIdx.resize(gi.size());
        for (int j=0; j<gi.size(); ++j)
        {
          globIdx[j] = std::lower_bound(globalIndices.begin(),globalIndices.end(),gi[j]) - globalIndices.begin();
          assert(globIdx[j]<globalIndices.size());
        }
 
        // Compute and scatter prolongation. First compute "p" as the matrix with p_kj the contribution
        // of the k-th child shape function to the j-th father shape function.
        // localProlongationMatrix(space,child,space,father,p,
        //                         space.mapper().shapefunctions(child),
        //                         space.mapper().shapefunctions(child)); // this makes only sense if the father and child shape functions are the same... TODO: replace this by father - however, as the father cell is not guaranteed to be a leaf cell, there need not be a father shape function set associated with the father cell. Maybe we can extract a "generic" shape function set from the father space?
 
 
        // Compute and scatter restriction
        // Create a geometry mapping from the child coordinate system (local) to father (global).
        MultiLevelLocalGeometry<Grid> mlGeometry(grid,child,father,MultiLevelLocalGeometry<Grid>::FatherIsGlobal,false);
        
        shapeFunctionSet = &fatherShapeFunctions;
        
        iNodes = fatherShapeFunctions.interpolationNodes();
        nodesInside.resize(iNodes.size());
        std::iota(begin(nodesInside),end(nodesInside),0);
        mlGeometry.local(iNodes,nodesInside);
 
        // Evaluate the global shape functions at the interpolation nodes...
        evaluateGlobalShapeFunctions(space,child,iNodes,afValues);
 
        // ... and translate them into ansatz function values.
        space.mapper().combiner(child,space.indexSet().index(child)).rightTransform(afValues);
        assert(afValues.M() == globIdx.size()); // make sure the number of ansatz functions is consistent
 
        for (int k=0; k<afValues.N(); ++k)
        {
          int rIdx = nodesInside[k];
          assert(0<=rIdx && rIdx<fCount.size());
          ++fCount[rIdx];
          for (int i=0; i<afValues.M(); ++i)
          {
            assert(!std::isnan(restrictionData[rIdx][globIdx[i]]));
            restrictionData[rIdx][globIdx[i]] += afValues[k][i];
            assert(!std::isnan(restrictionData[rIdx][globIdx[i]]));
          }
        }
      }
 
      // Do averaging of restriction of discontinuous FE functions.
      for (int i=0; i<restrictionData.N(); ++i)
      {
        assert(fCount[i]>0);
        double avg = 1.0 / fCount[i];
 
        for (int j=0; j<restrictionData.M(); ++j)
        {
          assert(!std::isnan(restrictionData[i][j]));
          restrictionData[i][j] *= avg;
        }
      }
 
      // Interpolate father shape functions.
      approximateGlobalValues(space,father,restrictionData,restrictionMatrix);
    }
 
    std::vector<size_t> globalIndices;
    LocalTransferMatrix restrictionMatrix;
    ShapeFunctionSet const* shapeFunctionSet;
  };
 
 
  //---------------------------------------------------------------------
 
 
  template <class Space, class CoarseningPolicy=AdaptationCoarseningPolicy>
  class TransferData 
  {
  public:
    typedef typename Space::Grid            Grid;
    typedef typename Space::GridView        GridView;
    typedef typename Space::IndexSet        IndexSet;
    typedef typename Space::Scalar          Scalar;
    typedef typename Grid::GlobalIdSet::IdType Id;
    typedef Dune::FieldVector<Scalar,Space::sfComponents>             SfValue;
    typedef typename Grid::template Codim<0>::Entity                  Cell;
 
  private:
    typedef DynamicMatrix<Dune::FieldMatrix<Scalar,1,1>>    RestrictionMatrix;
    typedef LocalTransfer<Space,CoarseningPolicy>           RestrictionData;
    typedef std::map<Id, RestrictionData>                   HistoryMap;
 
  public:
 
    class TransferMatrix
    {
    public:
      template <class StorageValue>
      std::unique_ptr<Dune::BlockVector<StorageValue>> apply(Dune::BlockVector<StorageValue> const& oldCoeff) const
      {
        std::unique_ptr<Dune::BlockVector<StorageValue>> newCoeff(new Dune::BlockVector<StorageValue>(tm.size()));
        for (int i=0; i<tm.size(); ++i)
        {
          StorageValue tmp(0.0);
          for (int j=0; j<tm[i].first.size(); ++j)
            tmp += tm[i].first[j].second * oldCoeff[tm[i].first[j].first];
          (*newCoeff)[i] = tmp;
        }
        return newCoeff;
      }
 
      TransferMatrix(size_t rows, size_t cols):
        tm(rows), nrows(rows), ncols(cols)
      {
        for (int i=0; i<tm.size(); ++i)
          tm[i].second = std::numeric_limits<int>::max(); // sentinel: restriction level
      }
 
    private:
 
      friend class TransferData;
 
      // Enters transfer matrix entries for a specific row.
      // This is done only if there have not yet been data entered with a smaller
      // restriction level. Transfer coefficients are potentially stored for a whole
      // cell hierarchy - we make sure that transfer is done only from the topmost
      // grid level.
      template <class ColIter, class EntryIter>
      void storeRow(int row, int restrictionLevel, ColIter firstCol, ColIter lastCol,
                    EntryIter firstEntry)
      {
        if (tm[row].second > restrictionLevel)
        {
          tm[row].second = restrictionLevel;
          tm[row].first.clear();
          tm[row].first.reserve(std::distance(firstCol,lastCol));
          while (firstCol != lastCol)
          {
            if (std::abs(*firstEntry)>1e-14)
              tm[row].first.push_back(std::make_pair(*firstCol,*firstEntry));
            ++firstCol;
            ++firstEntry;
          }
        }
      }
 
    public:
      void out(std::ostream& o) const
      {
        o << "shape=[" << nrows << ' ' << ncols << "];\n";
        o << "data=[\n";
        for (int i=0; i<tm.size(); ++i)
          for (int j=0; j<tm[i].first.size(); ++j)
            o << i+1 << ' ' << tm[i].first[j].first+1 << ' ' << tm[i].first[j].second << '\n';
        o << "];\n";
 
      }
 
      int nRows() const { return nrows; }
      int nCols() const { return ncols; }
 
 
    private:
      // Data for TransferMatrix - a vector with one entry for each row
      // tm.first: sparse column-vector <index,value>  tm.second: globalIdx
      std::vector<std::pair<std::vector<std::pair<int,Scalar> >,int>> tm;
      int nrows,ncols;
    }; // end of class TransferMatrix
 
 
 
    TransferData(Space const& space_, CoarseningPolicy const& mightBeCoarsened=CoarseningPolicy())
    : space(space_), DOFcoarseSpace(space.degreesOfFreedom())
    {
      auto const& idSet = space.grid().globalIdSet();
 
      typename Space::Evaluator evaluator(space);
 
      std::set<Id> processedAncestors;
 
 
      // Gather all information that is necessary to create local restriction matrices.
      // We reach down from the leaf cells to father cells. Since a father cell has
      // multiple children, father cells are visited several times, but need to be
      // processed only once. We thus check whether the father has already been processed.
      // This requires synchronization of the differen threads running in parallel on
      // leaf cell ranges.
      std::mutex mutex;
      auto const& cellRanges = space.gridManager().cellRanges(space.gridView());
      int nTasks = std::min(NumaThreadPool::instance().cpus(),cellRanges.maxRanges());
      parallelFor([&](int i, int n)
      {
        HistoryMap myHistoryData;                   // thread-local result data
 
        std::vector<std::pair<Id,Cell>> ancestors;  // for storing the ancestors to be processed in this thread
                                                    // (defined here to prevent frequent reallocations)
        for (auto const& cell: cellRanges.range(n,i))
        {
          // obtain and store only global indices of ansatz functions for all leaf cells.
          // This does not conflict with other leaf cells in different threads and need
          // not synchronized
          myHistoryData.insert(std::pair(idSet.id(cell),
                                         RestrictionData(space,cell,mightBeCoarsened)));
 
          // lower level cells might be important for grid transfer
          // todo@ !(entity->mightBeCoarsened) does not mean that this entity will exist after a grid change.
          // e.g. if an entity with the same father is marked for coarsening, then there are problems.
          {
            Cell ancestor = cell;
            std::lock_guard<std::mutex> lock(mutex);
 
            while (ancestor.mightVanish() && ancestor.level()>0)
            {
              ancestor = ancestor.father();
              Id id = idSet.id(ancestor);
 
              // If no data for the current ancestor exists, create some. Mark taking over responsibility
              // to prevent concurrent threads from processing the same ancestor.
              if (processedAncestors.find(id) == processedAncestors.end())
              {
                processedAncestors.insert(id);
                ancestors.push_back(std::pair(id,ancestor));
              }
            }
          }
 
          // Now process all the ancestors in our responsibility.
          for (auto const& [id,ancestor]: ancestors)
            myHistoryData.insert(std::pair(id,RestrictionData(space,ancestor,
                                                              mightBeCoarsened)));
          ancestors.clear();
        }
 
        // Now we need to move our local leaf-level history data over to the global history data map.
        // This requires synchronization.
        std::lock_guard<std::mutex> lock(mutex);
        historyData.merge(myHistoryData);
      },nTasks);
    }
 
 
    std::unique_ptr<TransferMatrix> transferMatrix() const
    {
      typename Grid::GlobalIdSet const& idSet = space.grid().globalIdSet();
      auto transfer = std::make_unique<TransferMatrix>(space.degreesOfFreedom(),DOFcoarseSpace);
 
 
      for (auto cell: elements(space.gridView()))
      {
        // Among the ancestors of the current cell that have a
        // restriction matrix stored, we select the one with the highest
        // refinement level (i.e. the most exact representation of the
        // FE functions before mesh adaptation.
        Cell ancestor(cell);
        typename HistoryMap::const_iterator restriction;
        int colevel = 0;  // The restriction level, i.e. the distance from the target grid view.
        while((restriction=historyData.find(idSet.id(ancestor))) == historyData.end() && ancestor.level()> 0)
        {
          colevel++;
          ancestor = ancestor.father();
        }
        assert(restriction != historyData.end());
        RestrictionMatrix const& rm = restriction->second.restrictionMatrix;
 
        auto index = space.indexSet().index(cell);
 
        RestrictionMatrix prolongation;
        // WARNING: having the same shape function set for child and father makes only sense if these are the same...
        //localProlongationMatrix(space,Cell(ci),space,ancestor,prolongation,*(restriction->second.shapeFunctionSet),*(restriction->second.shapeFunctionSet));
        // Then let's try the following.
        localProlongationMatrix(space,cell,space,ancestor,prolongation,
                                space.mapper().shapefunctions(index),*(restriction->second.shapeFunctionSet));
 
        // compute matrix product prolongation * restriction
        RestrictionMatrix localTransfer;
        MatMult(localTransfer,prolongation,rm);
 
 
        // apply K^+
        space.mapper().combiner(cell,index).leftPseudoInverse(localTransfer);
 
 
 
        // scatter the local transfer
        auto rowIdx = space.mapper().globalIndices(index);
        std::vector<size_t> const& columnIdx = restriction->second.globalIndices;
        for (int i=0; i<localTransfer.N(); ++i)
        {
          transfer->storeRow(rowIdx[i], colevel,
                             columnIdx.begin(), columnIdx.end(),
                             localTransfer[i].begin());
        }
      }
 
      return transfer;
    }
 
  private:
    Space const&  space;
    HistoryMap    historyData;
    int           DOFcoarseSpace;
  };
 
 
  // ----------------------------------------------------------------------------------------------
  
  struct InverseVolume
  {
    template<class Entity>
    double operator()(Entity const& e) {return 1.0/e.geometry().volume();}
    double scalingFactor() {return 1.0;}
    double offset() {return 0.0;}
  };
 
  struct Volume
  {
    template<class Entity>
    double operator()(Entity const& e) {return e.geometry().volume();}
    double scalingFactor() {return 1.0;}
    double offset() {return 0.0;}
  };
 
  struct PlainAverage
  {
    template<class Entity>
    double operator()(Entity const& e) {return 1.0;}
    double scalingFactor() {return 1.0;}
    double offset() {return 0.0;}
  };
 
  template<int scl=1>
  struct SumUp
  {
    template<class Entity>
    double operator()(Entity const& e) {return 1.0;}
    double scalingFactor() {return 0.0;}
    double offset() {return scl;}
  };
 
 
  template<typename Weighing=PlainAverage, typename FSElement, typename Function>
  void interpolateGlobally(FSElement& fse, Function const& fu)
  {
    typedef typename FSElement::Space ImageSpace;
    typedef typename ImageSpace::Grid Grid;
 
    DynamicMatrix< Dune::FieldMatrix<typename ImageSpace::Scalar,ImageSpace::sfComponents,1>> globalValues;
    DynamicMatrix< Dune::FieldMatrix<typename ImageSpace::Scalar, 1, 1 >> localCoefficients;
    std::vector<typename Grid::ctype> sumOfWeights(fse.space().degreesOfFreedom(),0.0);
 
    // Initialize the target to zero.
    fse.coefficients() = typename FSElement::Scalar(0.0);
 
    typename ImageSpace::Evaluator isfs(fse.space());
    typename Function::Space::Evaluator dsfs(fu.space());
    Weighing w;
    
    using DirectValueType = decltype(fu.value(dsfs));
    using ValueType = std::conditional_t<std::is_arithmetic<DirectValueType>::value,
                                         Dune::FieldVector<DirectValueType,1>,
                                         DirectValueType>;
    std::vector<ValueType> fuvalue; // function values - declare here to prevent reallocations
 
    // Step through all the cells and perform local interpolation on each cell.
    for (auto const& cell: elements(fse.space().gridView())) 
    {
      auto index = fse.space().indexSet().index(cell);
      
      isfs.moveTo(cell);
      dsfs.moveTo(cell);
 
      typename Grid::ctype myWeight = w(cell);
      for (int i=0; i<isfs.size(); ++i) 
        sumOfWeights[isfs.globalIndices()[i]] += myWeight;
 
      // Obtain interpolation nodes of target on this cell
      auto const& iNodes(isfs.shapeFunctions().interpolationNodes());
 
      // Evaluate source
      globalValues.setSize(iNodes.size(),1);
      localCoefficients.setSize(iNodes.size(),1);
      fuvalue.resize(iNodes.size());
      for (int i=0; i<iNodes.size(); ++i) 
      {
        dsfs.evaluateAt(iNodes[i]);
        fuvalue[i] = fu.value(dsfs);
      }
 
      // Perform local interpolation
      for(int k=0; k<FSElement::components/ImageSpace::sfComponents; ++k) 
      {
        for (int i=0; i<iNodes.size(); ++i)
          for(int j=0; j<ImageSpace::sfComponents; ++j)
            globalValues[i][0][j][0] = myWeight*fuvalue[i][ImageSpace::sfComponents*k+j];
 
        approximateGlobalValues(fse.space(),cell,globalValues,localCoefficients);
 
        assert(localCoefficients.N() == isfs.globalIndices().size());
 
        // transform global values onto shape function coefficients
        fse.space().mapper().combiner(cell,index).leftPseudoInverse(localCoefficients);
 
        // add up the shape function coefficients to the ansatz function coefficients
        for (int i=0; i<localCoefficients.N(); ++i) 
          fse.coefficients()[isfs.globalIndices()[i]][k] += localCoefficients[i][0];
      }
    }
 
    // On vertices, edges, and faces, degrees of freedom may have received
    // updates from multiple cells. Average this out.
    for(int i=0; i<sumOfWeights.size(); ++i)
      fse.coefficients()[i] /= (sumOfWeights[i]*w.scalingFactor()+w.offset());
  }
 
 
 
  template<typename Weighing=PlainAverage, typename Target, typename Source>
  void interpolateGloballyWeak(Target& target, Source const& fu)
  {
    typedef typename Target::Space ImageSpace;
    typedef typename ImageSpace::Grid Grid;
    using Scalar = typename ImageSpace::Scalar;
 
    auto const& space = target.space();
 
    DynamicMatrix< Dune::FieldMatrix<Scalar, ImageSpace::sfComponents, 1>> globalValues;
    DynamicMatrix< Dune::FieldMatrix<Scalar, 1, 1>> localCoefficients;
    std::vector<typename Grid::ctype> sumOfWeights(space.degreesOfFreedom(),0.0);
 
    // Initialize the target to zero.
    target = 0.0;
 
    typename ImageSpace::Evaluator isfs(space);
    Weighing w;
 
 
    using ValueType = decltype(fu.value(*space.gridView().template end<0>(),
                                        Dune::FieldVector<typename Grid::ctype, ImageSpace::dim>()));
    static_assert(std::is_same_v<ValueType,typename Target::ValueType>,
                  "return type of the weak function view must agree with the value type of the assignment target");
 
    std::vector<ValueType> fuvalue; // function values - declare here to prevent reallocations
 
    for (auto const& cell: elements(space.gridView()))
    {
      auto index = space.indexSet().index(cell);
      isfs.moveTo(cell);
 
      typename Grid::ctype myWeight = w(cell);
      for (int i=0; i<isfs.size(); ++i) 
        sumOfWeights[isfs.globalIndices()[i]] += myWeight;
 
      // Obtain interpolation nodes of target on this cell
      auto const& iNodes(isfs.shapeFunctions().interpolationNodes());
 
      // Evaluate source
      globalValues.setSize(iNodes.size(),1);
      localCoefficients.setSize(iNodes.size(),1);
      fuvalue.resize(iNodes.size());
      for (int i=0; i<iNodes.size(); ++i)
        fuvalue[i] = fu.value(cell,iNodes[i]);
 
      // Perform local interpolation
      for(int k=0; k<Target::components/ImageSpace::sfComponents; ++k)
      {
        for (int i=0; i<iNodes.size(); ++i)
          for(int j=0; j<ImageSpace::sfComponents; ++j)
            globalValues[i][0][j][0] = myWeight*fuvalue[i][ImageSpace::sfComponents*k+j];
 
        approximateGlobalValues(space,cell,globalValues,localCoefficients);
     
        if (isfs.globalIndices().size()!=0)
        {
          assert(localCoefficients.N() == isfs.globalIndices().size());
 
          // transform global values onto shape function coefficients
          space.mapper().combiner(cell,index).leftPseudoInverse(localCoefficients);
 
          // add up the shape function coefficients to the ansatz function coefficients
          for (int i=0; i<localCoefficients.N(); ++i)
            target.coefficients()[isfs.globalIndices()[i]][k]+= localCoefficients[i][0];
        }
      }
    }
    
    // On vertices, edges, and faces, degrees of freedom may have received
    // updates from multiple cells. Average this out.
    for(int i=0; i<sumOfWeights.size(); ++i)
      target.coefficients()[i] /= sumOfWeights[i]*w.scalingFactor()+w.offset();
  }
 
 
  
  // ----------------------------------------------------------------------------------------------
 
 
 
  template <class Functor>
  class WeakFunctionViewAdaptor
  {
  public:
    WeakFunctionViewAdaptor(Functor const& g_)
    : g(g_)
    {}
 
    template <class ...Args>
    auto value(Args... args) const
    {
      return g(std::forward<Args>(args)...);
    }
 
  private:
    Functor g;
  };
 
  template <class Functor>
  auto makeWeakFunctionView(Functor const& g)
  {
    return WeakFunctionViewAdaptor<Functor>(g);
  }
 
  template <class Space_, class Functor>
  class FunctionViewAdaptor
  {
  public:
    using Space = Space_;
    
    FunctionViewAdaptor(Space const& space, Functor const& g_): g(g_), sp(space)
    {}
    
    Space const& space() const
    {
      return sp;
    }
    
    template <class ...Args>
    auto value(Args... args) const 
    {
      return g(std::forward<Args>(args)...);
    }
    
    
  private:
    Functor g;
    Space const&  sp;
  };
  
  template <class Space, class Functor>
  auto makeFunctionView(Space const& space, Functor const& g)
  {
    return FunctionViewAdaptor<Space,Functor>(space,g);
  }
  
  // ----------------------------------------------------------------------------------------------
  
  template <class Fu1, class Fu2>
  [[deprecated]]
  void spaceTransfer(Fu1& f1, Fu2 const& f2)
  {
    // typedef TransferData<typename Fu2::Space> TrData;
    // 
    // TrData trdata(f2.space());
    // std::unique_ptr<typename TrData::TransferMatrix> trMat=trdata.transferMatrix(f1.space());
    // 
    // 
    // (*f1) = *(trMat->apply(*f2));
    interpolateGlobally(*f1,*f2);
  }
  
  
  
  
 
} /* end of namespace Kaskade */
#endif