lagrangebasis.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH
 
#include <array>
#include <type_traits>
#include <vector>
 
#include <dune/common/exceptions.hh>
#include <dune/common/typelist.hh>
#include <dune/common/indices.hh>
#include <dune/common/math.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/std/mdarray.hh>
 
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>
 
#include <dune/localfunctions/common/localkey.hh>
#include <dune/localfunctions/lagrange/lagrangecube.hh>
#include <dune/localfunctions/lagrange/lagrangeprism.hh>
#include <dune/localfunctions/lagrange/lagrangepyramid.hh>
#include <dune/localfunctions/lagrange/lagrangesimplex.hh>
#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
 
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/mcmgmapper.hh>
 
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/leafprebasismixin.hh>
 
#include <dune/grid/common/capabilities.hh>
 
namespace Dune {
namespace Functions {
 
 
 
namespace Experimental {
 
template<std::size_t dim>
class FaceOrientations
{
  static constexpr std::size_t maxEdges = 12;
  static constexpr std::size_t maxFacets = 6;
  static constexpr std::size_t bitsPerFacet = 4;
 
  using Data = std::bitset<maxEdges + maxFacets*bitsPerFacet>;
 
  static constexpr Data singleFacetMask = (1<<bitsPerFacet) - 1;
 
  // Offset of the 2d facet within the bitfield
  static constexpr auto facetBitOffset(int subEntity)
  {
    return maxEdges + bitsPerFacet*subEntity;
  }
 
  template<unsigned int blockSize, class T>
  static constexpr T getBitBlock(const T& bits, unsigned int i)
  {
    constexpr T blockMask = (1 << (blockSize)) - 1;
    return (bits >> (blockSize*i)) & blockMask;
  }
 
  constexpr void computeEdgeOrientation(const auto& vertexIndices, const auto& re, int subEntity)
  {
    data_[subEntity] = vertexIndices[0] > vertexIndices[1];
  }
 
  constexpr void computeTriangleOrientation(const auto& vertexIndices, const auto& re, int subEntity)
  {
    // All possible triangle orientations are stored in a binary encoded lookup table.
    // We map the orientations of the three adjacent edges to a flat
    // index in 0,...,7 and use it to access a 4 bit block in the
    // table encoding the triangle orientation.
    // The leading bit is always zero to indicate that the face is a triangle.
    constexpr uint32_t edgeOrientationsToTriangleOrientation = 0b0101'0001'0011'0100'0010'0011'0110'0000;
    auto edgeOrientations = (data_[re.subEntity(subEntity, 1, 0, 2)] << 0)
                          | (data_[re.subEntity(subEntity, 1, 1, 2)] << 1)
                          | (data_[re.subEntity(subEntity, 1, 2, 2)] << 2);
    auto faceOrientation = getBitBlock<4>(edgeOrientationsToTriangleOrientation, edgeOrientations);
    // Set bits associated to this facet.
    data_ |= (Data(faceOrientation) << facetBitOffset(subEntity));
  }
 
  constexpr void computeQuadrilateralOrientation(const auto& vertexIndices, const auto& re, auto subEntity)
  {
    // The following computes the local index i_min of the quadrilateral
    // vertex with smallest index. I.e. it is equivalent to
    //
    //   std::size_t i_min = 0;
    //   for(auto i: Dune::range(1, 4))
    //     if (vertexIndices[i] < vertexIndices[i_min])
    //       i_min = i;
    //
    // This code would always do 3 index comparisons although several of them have
    // already been done for the edges. Hence we rather want to reuse the latter
    // to save comparisons (for large index/id types). The following code maps
    // the 16 possible outcomes of the edge comparisons to the index i which
    // is then used to lookup i_min from a precomputed index table. Since i_min
    // is from 0,..,3 we can represent the value by two bits and encode the
    // lookup table as 16*2 bits. There are two cases (5 and 10), where i_min
    // cannot be deduced from edge comparisons because the two vertices with the
    // smallest index are not neighbors. These cases require an additional comparison
    // of the indices of those two vertices.
    //
    // In contrast to the linear search above this reduces the total number of comparisons
    // for a hexahedron from 36=12+6*4 to a number between 18=12+6 and 24=12+6*2,
    // depending on the number of required non-neighbor comparisons.
    //
    // For index/id types that are cheap to compare like e.g. int, the code does
    // not make a measurable difference. For larger types it is slightly faster.
    constexpr uint32_t edgeOrientationsToMinVertex = 0b11'11'01'01'11'00'00'00'10'00'00'01'10'00'10'00;
    auto edgeOrientations = (data_[re.subEntity(subEntity, 1, 0, 2)] << 0)
                          | (data_[re.subEntity(subEntity, 1, 1, 2)] << 1)
                          | (data_[re.subEntity(subEntity, 1, 2, 2)] << 2)
                          | (data_[re.subEntity(subEntity, 1, 3, 2)] << 3);
    std::size_t i_min = 0;
    // In case 5 and 10 we need to do an extra comparison, because
    // we cannot deduce i_min from edge orientations. In all other
    // cases, we can lookup i_min by extracting two bits from the table.
    if (edgeOrientations==5)
      i_min = (vertexIndices[1] < vertexIndices[2]) ? 1 : 2;
    else if (edgeOrientations==10)
      i_min = (vertexIndices[0] < vertexIndices[3]) ? 0 : 3;
    else
      i_min = getBitBlock<2>(edgeOrientationsToMinVertex, edgeOrientations);
    // The lowest two bits of the faceOrientation indicate the number of
    // counter clockwise rotations needed to map vertex i_min to vertex 0.
    // The third bit indicates if we need a reflection to revert the order.
    // The leading bit is always one to indicate that the face is a quadrilateral.
    unsigned long faceOrientation = 0;
    if (i_min==0)
      faceOrientation = 0b1000 | ((vertexIndices[2] < vertexIndices[1]) << 2);
    else if (i_min==1)
      faceOrientation = 0b1011 | ((vertexIndices[0] < vertexIndices[3]) << 2);
    else if (i_min==2)
      faceOrientation = 0b1001 | ((vertexIndices[3] < vertexIndices[0]) << 2);
    else if (i_min==3)
      faceOrientation = 0b1010 | ((vertexIndices[1] < vertexIndices[2]) << 2);
    // Set bits associated to this facet.
    data_ |= (Data(faceOrientation) << facetBitOffset(subEntity));
  }
 
public:
 
  FaceOrientations() = default;
 
  template <class VertexIds, std::size_t... codims>
  FaceOrientations(const Dune::GeometryType& geometryType, const VertexIds& vertexIds, const Dune::index_constant<codims>&... codimensions)
    : data_(0)
  {
    if constexpr ((dim == 1) or (sizeof...(codims)==0))
      return;
 
    // Get reference element
    const auto& re = Dune::referenceElement<double,dim>(geometryType);
 
    // Access and store vertex ids only once
    using Index = std::decay_t<decltype(vertexIds[0])>;
    auto cachedVertexIds = std::array<Index, Dune::power(2, dim)>{};
    for(auto k : Dune::range(re.size(dim)))
      cachedVertexIds[k] = vertexIds[k];
 
    // Subrange of face vertex ids
    auto faceVertexIds = [&](auto faceIndex, auto faceCodim) {
      return Dune::transformedRangeView(re.subEntities(faceIndex, faceCodim, dim), [&](auto localVertexIndex) {
        return cachedVertexIds[localVertexIndex];
      });
    };
 
    if constexpr ((dim > 1) and ((codims<=(dim-1)) || ...))
    {
      for(auto subEntity : Dune::range(re.size(dim-1)))
        computeEdgeOrientation(faceVertexIds(subEntity, dim-1), re, subEntity);
    }
    if constexpr ((dim == 3) and ((codims==(dim-2)) || ...))
    {
      for(auto subEntity : Dune::range(re.size(1)))
      {
        if (re.type(subEntity, 1).isTriangle())
          computeTriangleOrientation(faceVertexIds(subEntity, 1), re, subEntity);
        if (re.type(subEntity, 1).isQuadrilateral())
          computeQuadrilateralOrientation(faceVertexIds(subEntity, 1), re, subEntity);
      }
    }
  }
 
  std::size_t faceOrientationIndex(std::size_t subEntity, std::size_t codim) const
  {
    if (codim == dim-1)
      return data_[subEntity];
    if (codim == 1)
      return ((data_ >> facetBitOffset(subEntity)) & singleFacetMask).to_ulong();
    return 0;
  }
 
private:
 
  Data data_;
};
 
 
 
template<class IndexIdSet>
class LagrangeFaceDOFPermutation
{
  static constexpr int dim = IndexIdSet::dimension;
 
  // The following methods pre-compute the tables for renumbering the face DOFs.
  // The rows of a table correspond to the local orientations which are indexed
  // using the face orientation index as provided by the FaceOrientations class.
  // Each row stores the local indices of the face DOF indices wrt the global orientation.
  // I.e. given the table T, a local face DOF index i with respect to the
  // local orientation j, the globally unique index of this DOF with respect to
  // the global face orientation is T[j][i].
  //
  // The computation of the table is based on the following observation:
  // Assume that F_j is the transformation mapping the global orientation
  // of the face to its local orientation with index j. Now let p_i an
  // interpolation point in the face with local index i. Then the
  // globally oriented index of the basis function associated to p_i
  // is obtained as the local index of the point F_j(p_i). Hence we can
  // compute the renumbering table T[j][i] as follows:
  // For given local orientation j and locally oriented DOF index i
  // we first compute the point p_i, then we apply the transformation
  // F_j (not its inverse as one might guess at first sight) which
  // consists of rotations and zero or one reflection and use the index
  // of the obtained point.
 
  // Compute table of all possible permuted triangle DOF indices
  static auto globallyOrientedTriangleDOFTable(int order)
  {
    auto dofPermutation = std::array<std::vector<std::size_t>, 8>();
    // Number of DOFs on the base line
    std::size_t m = order-2;
    // Total number of DOFs
    std::size_t n = m*(m+1)/2;
    if (order<3)
      return dofPermutation;
    // In order to reflect we need to apply the transformation
    // from global to local orientation as documented in the
    // FaceOrientations class for each face orientation index.
    // To do this we first map flat DOF indices to barycentric
    // multi-indices wrt. the vertices and summing up to m-1.
    // Then all transformations can be represented by flips
    // of the components for the barycentric multi-index.
    using BarycentricIndex = std::array<std::size_t, 3>;
    auto barycentricIndices = std::vector<BarycentricIndex>();
    for(auto i : Dune::range(m))
      for(auto j : Dune::range(m-i))
        barycentricIndices.push_back(BarycentricIndex{m-1-(i+j), j, i});
    auto flatToBarycentric = [&](auto i) { return barycentricIndices[i]; };
    auto barycentricToFlat = [&](auto b) { return b[1] + m*b[2] - b[2]*(b[2]-1)/2; };
    // Loop of all 8 orientations
    for(auto j : Dune::range(8))
    {
      dofPermutation[j].resize(n);
      for(auto i : Dune::range(n))
      {
        auto b = flatToBarycentric(i);
        // Bit 0 and 1: Number of counter-clockwise vertex-wise rotations
        if ((j & 0b011) == 1)
        {
          // The counter-clockwise rotation by one vertex can be decomposed
          // into two reflections.
          std::swap(b[0], b[1]);
          std::swap(b[0], b[2]);
        }
        if ((j & 0b011) == 2)
        {
          // The counter-clockwise rotation by two vertices is the inverse
          // of the rotation by one vertex and can thus be decomposed into
          // the same two reflections in opposite order.
          std::swap(b[0], b[2]);
          std::swap(b[0], b[1]);
        }
        // Bit 2: Reflect across xy-diagonal
        if (j & 0b100)
          std::swap(b[1], b[2]);
        dofPermutation[j][i] = barycentricToFlat(b);
      }
    }
    return dofPermutation;
  }
 
  // Compute table of all possible permuted quadrilateral DOF indices
  static auto globallyOrientedQuadrilateralDOFTable(int order)
  {
    auto dofPermutation = std::array<std::vector<std::size_t>, 8>();
    if (order<2)
      return dofPermutation;
    // Number of DOFs on the base line
    std::size_t m = order-1;
    // Total number of DOFs
    std::size_t n = m*m;
    // In order to reflect we map the flat DOF index into a cartesian
    // multi-index wrt. the x- and y-axis.
    using CartesianIndex = std::array<std::size_t, 2>;
    auto flatToCartesian = [&](auto i) { return CartesianIndex{i % m, i / m}; };
    auto cartesianToFlat = [&](auto c) { return c[0] + m*c[1]; };
    for(auto j : Dune::range(8))
    {
      dofPermutation[j].resize(n);
      for(auto i : Dune::range(n))
      {
        auto c = flatToCartesian(i);
        // Bit 0 and 1: Number of counter-clockwise vertex-wise rotations
        if ((j & 0b011) == 1)
        {
          std::swap(c[0], c[1]);
          c[0] = m-1-c[0];
        }
        if ((j & 0b011) == 2)
        {
          c[0] = m-1-c[0];
          c[1] = m-1-c[1];
        }
        if ((j & 0b011) == 3)
        {
          c[0] = m-1-c[0];
          std::swap(c[0], c[1]);
        }
        // Bit 2: Reflect across xy-diagonal
        if (j & 0b100)
          std::swap(c[0], c[1]);
        dofPermutation[j][i] = cartesianToFlat(c);
      }
    }
    return dofPermutation;
  }
 
  std::size_t order() const
  {
    return order_;
  }
 
  /*
   * Type used to store local per-face indices.
   *
   * Since the maximal number of face DOFs is n=(order-1)*(order-1),
   * this type must be able to represent the range 0,...,n-1.
   */
  using FaceIndexType = unsigned char;
 
public:
 
  static constexpr unsigned int maxOrder3d = 17;
 
  LagrangeFaceDOFPermutation(const IndexIdSet& indexSet, int order)
    : indexIdSet_(&indexSet)
    , order_(order)
    , facetFlipTable_((dim < 3) ? 0 : std::max(0, (order-1)*(order-1)))
  {
    if constexpr (dim == 3)
    {
      // Fill table with all possible permuted 2d DOF indices.
      // This is indexed with the faceOrientationIndex which first
      // enumerates the orientations of triangles and then the ones
      // of quadrilateral.
      auto triangleFlipTable = globallyOrientedTriangleDOFTable(order);
      auto quadrilateralFlipTable = globallyOrientedQuadrilateralDOFTable(order);
      for(std::size_t k : Dune::range(8))
        for(std::size_t i : Dune::range(triangleFlipTable[k].size()))
          facetFlipTable_(k,i) = triangleFlipTable[k][i];
      for(std::size_t k : Dune::range(8))
        for(std::size_t i : Dune::range(quadrilateralFlipTable[k].size()))
          facetFlipTable_(8+k,i) = quadrilateralFlipTable[k][i];
    }
  }
 
  template <class Element>
  FaceOrientations<dim> computeFaceOrientations(const Element& element) const
  {
    const auto& re = Dune::referenceElement<double,dim>(element.type());
    auto vertexIds = Dune::transformedRangeView(re.subEntities(0,0,dim), [&](auto localVertexIndex) {
      if constexpr (requires { indexIdSet_->subIndex(element, 0, dim); })
        return indexIdSet_->subIndex(element, localVertexIndex, dim);
      else if constexpr (requires { indexIdSet_->subId(element, 0, dim); })
        return indexIdSet_->subId(element, localVertexIndex, dim);
    });
    using namespace Dune::Indices;
    auto k = order();
    if ((dim==1) or (k<3))
      return FaceOrientations<dim>(element.type(), vertexIds);
    else if ((k==3) and element.type().isTetrahedron())
      return FaceOrientations<dim>(element.type(), vertexIds, _2);
    else
      return FaceOrientations<dim>(element.type(), vertexIds, _2, _1);
  }
 
  unsigned int permuteFaceDOF(unsigned int subEntity, unsigned int codim, unsigned int index, const FaceOrientations<dim>& orientations) const
  {
    auto k = order();
    if constexpr(dim == 1)
      return index;
    if constexpr(dim == 2)
    {
      // Flip edge DOFs if reorientation is needed
      if ((k>2) and (codim == 1) and orientations.faceOrientationIndex(subEntity, codim))
        return k-2-index;
      return index;
    }
    if constexpr(dim == 3)
    {
      if (k>2)
      {
        // Flip edge DOFs if reorientation is needed
        if ((codim == 2) and orientations.faceOrientationIndex(subEntity, codim))
          return k-2-index;
        // Permute 2d facet DOFs according to lookup table
        if (codim == 1)
          return facetFlipTable_(orientations.faceOrientationIndex(subEntity, codim), index);
      }
      return index;
    }
  }
 
  unsigned int permuteFaceDOF(const Dune::LocalKey& localKey, const FaceOrientations<dim>& orientations) const
  {
    return permuteFaceDOF(localKey.subEntity(), localKey.codim(), localKey.index(), orientations);
  }
 
private:
  const IndexIdSet* indexIdSet_;
  std::size_t order_;
  Std::mdarray<FaceIndexType, Std::extents<std::size_t, 16, std::dynamic_extent>> facetFlipTable_;
};
 
} // namespace Experimental
 
 
 
// *****************************************************************************
// This is the reusable part of the LagrangeBasis. It contains
//
//   LagrangePreBasis
//   LagrangeNode
//
// The pre-basis allows to create the others and is the owner of possible shared
// state. These components do _not_ depend on the global basis and local view
// and can be used without a global basis.
// *****************************************************************************
 
template<typename GV, int k, typename R=double>
class LagrangeNode;
 
template<typename GV, int k, typename R=double>
class LagrangePreBasis;
 
 
 
template<typename GV, int k, typename R>
class LagrangePreBasis :
  public LeafPreBasisMixin< LagrangePreBasis<GV,k,R> >
{
  static const int dim = GV::dimension;
  static const bool useDynamicOrder = (k<0);
 
  static Dune::MCMGLayout dofLayout(unsigned int order)
  {
    return [order](Dune::GeometryType type, int dim) -> std::size_t {
      if (order==0)
        return type.dim() == dim;
      if (type.isSimplex())
        return Dune::binomial(order-1, type.dim());
      if (type.isCube())
        return Dune::power(order-1, type.dim());
      if (type.isPyramid())
        return 0;
      if (type.isPrism() and (order>1))
        return (order-1)*(order-1)*(order-2)/2;
      return 0;
    };
  }
 
  using FaceDOFPermutation = Experimental::LagrangeFaceDOFPermutation<typename GV::Grid::GlobalIdSet>;
 
public:
 
  using GridView = GV;
 
  using size_type = std::size_t;
 
  using Node = LagrangeNode<GV, k, R>;
 
  LagrangePreBasis(const GridView& gv)
    : LagrangePreBasis(gv, std::numeric_limits<unsigned int>::max())
  {}
 
  LagrangePreBasis(const GridView& gv, unsigned int runTimeOrder)
    : gridView_(gv)
    , order_(runTimeOrder)
    , mapper_(gridView_, dofLayout(useDynamicOrder ? runTimeOrder : k))
    , faceDOFPermutation_(gridView_.grid().globalIdSet(), useDynamicOrder ? runTimeOrder : k)
  {
    if (!useDynamicOrder && runTimeOrder!=std::numeric_limits<unsigned int>::max())
      DUNE_THROW(RangeError, "Template argument k has to be -1 when supplying a run-time order!");
    if ((order() > 2) and (gridView_.indexSet().size(Dune::GeometryTypes::pyramid) > 0))
      DUNE_THROW(RangeError, "Polynomial order >2 is not supported for grids containing pyramid elements.");
    if ((order() > 2) and (gridView_.indexSet().size(Dune::GeometryTypes::prism) > 0))
      DUNE_THROW(RangeError, "Polynomial order >2 is not supported for grids containing prism elements.");
    if ((dim == 3) and (order() > FaceDOFPermutation::maxOrder3d))
      DUNE_THROW(RangeError, "Polynomial order >" << FaceDOFPermutation::maxOrder3d << " is not supported in 3d");
  }
 
  void initializeIndices()
  {}
 
  const GridView& gridView() const
  {
    return gridView_;
  }
 
  void update (const GridView& gv)
  {
    if ((order() > 2) and (gv.indexSet().size(Dune::GeometryTypes::pyramid) > 0))
      DUNE_THROW(RangeError, "Polynomial order >2 is not supported for grids containing pyramid elements.");
    if ((order() > 2) and (gv.indexSet().size(Dune::GeometryTypes::prism) > 0))
      DUNE_THROW(RangeError, "Polynomial order >2 is not supported for grids containing prism elements.");
    gridView_ = gv;
    mapper_.update(gridView_);
    faceDOFPermutation_ = FaceDOFPermutation(gridView_.grid().globalIdSet(), order());
  }
 
  Node makeNode() const
  {
    return Node{order()};
  }
 
  size_type dimension() const
  {
    return mapper_.size();
  }
 
  size_type maxNodeSize() const
  {
    // That cast to unsigned int is necessary because GV::dimension is an enum,
    // which is not recognized by the power method as an integer type...
    return power(order()+1, (unsigned int)GV::dimension);
  }
 
  template<class Node, class It>
  It indices(const Node& node, It it) const
  {
    const auto& element = node.element();
    const auto& localCoefficients = node.finiteElement().localCoefficients();
 
    // Short-cut for order 1 since directly using the IndexSet is cheaper than the Mapper.
    if (order() == 1)
    {
      for(auto localIndex : Dune::range(localCoefficients.size()))
      {
        auto localKey = localCoefficients.localKey(localIndex);
        auto globalIndex = mapper_.gridView().indexSet().subIndex(element, localKey.subEntity(), localKey.codim());
        *it = {{ (size_type)(globalIndex) }};
        ++it;
      }
      return it;
    }
 
    // Precompute orientations of all faces
    auto faceOrientations = faceDOFPermutation_.computeFaceOrientations(element);
    for(auto localIndex : Dune::range(localCoefficients.size()))
    {
      auto localKey = localCoefficients.localKey(localIndex);
      auto globalIndex = mapper_.subIndex(element, localKey.subEntity(), localKey.codim());
      globalIndex += faceDOFPermutation_.permuteFaceDOF(localKey, faceOrientations);
      *it = {{ (size_type)(globalIndex) }};
      ++it;
    }
    return it;
  }
 
  unsigned int order() const
  {
    return (useDynamicOrder) ? order_ : k;
  }
 
protected:
  GridView gridView_;
 
  // Run-time order, only valid if k<0
  unsigned int order_;
 
  Dune::MultipleCodimMultipleGeomTypeMapper<GridView> mapper_;
  FaceDOFPermutation faceDOFPermutation_;
};
 
 
 
template<typename GV, int k, typename R>
class LagrangeNode :
  public LeafBasisNode
{
  static constexpr int dim = GV::dimension;
 
  // A simple cache handing storing exactly one LFE.  This can be
  // used for grids with only a single element type. In contrast to
  // StaticLagrangeLocalFiniteElementCache this also supports
  // Lagrange*LocalFiniteElement with run-time order.
  template <class LFE>
  class SingleLocalFiniteElementCache
  {
    LFE lfe_;
  public:
    using FiniteElementType = LFE;
 
    template<class... Args>
    SingleLocalFiniteElementCache(Args&&... args)
      : lfe_(std::forward<Args>(args)...)
    {}
 
    const FiniteElementType& get ([[maybe_unused]] Dune::GeometryType type) const
    {
      return lfe_;
    }
  };
 
  // Utility function to construct the FiniteElementCache type.
  // Since the function is just a helper to generate a type,
  // it hands out a MetaType<T> instead of a raw T.
  static constexpr auto makeCacheType()
  {
    using D = typename GV::ctype;
    if constexpr (Dune::Capabilities::hasSingleGeometryType<typename GV::Grid>::v)
    {
      constexpr auto type = Dune::GeometryType(Dune::Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId, GV::dimension);
      if constexpr (type.isSimplex())
        return Dune::MetaType<SingleLocalFiniteElementCache<Dune::LagrangeSimplexLocalFiniteElement<D,R,dim,k>>>{};
      else if constexpr (type.isCube())
        return Dune::MetaType<SingleLocalFiniteElementCache<Dune::LagrangeCubeLocalFiniteElement<D,R,dim,k>>>{};
      else if constexpr (type.isPrism())
        return Dune::MetaType<SingleLocalFiniteElementCache<Dune::LagrangePrismLocalFiniteElement<D,R,k>>>{};
      else if constexpr (type.isPyramid())
        return Dune::MetaType<SingleLocalFiniteElementCache<Dune::LagrangePyramidLocalFiniteElement<D,R,k>>>{};
    }
    else
      return Dune::MetaType<Dune::LagrangeLocalFiniteElementCache<D,R,dim,k>>{};
  }
 
  using FiniteElementCache = typename decltype(makeCacheType())::type;
 
public:
 
  using size_type = std::size_t;
  using Element = typename GV::template Codim<0>::Entity;
  using FiniteElement = typename FiniteElementCache::FiniteElementType;
 
  LagrangeNode() :
    LagrangeNode(k)
  {}
 
  LagrangeNode(unsigned int order) :
    cache_(order),
    finiteElement_(nullptr),
    element_(nullptr)
  {}
 
  const Element& element() const
  {
    return *element_;
  }
 
  const FiniteElement& finiteElement() const
  {
    return *finiteElement_;
  }
 
  void bind(const Element& e)
  {
    element_ = &e;
    finiteElement_ = &(cache_.get(element_->type()));
    this->setSize(finiteElement_->size());
  }
 
protected:
 
  FiniteElementCache cache_;
  const FiniteElement* finiteElement_;
  const Element* element_;
};
 
 
 
namespace BasisFactory {
 
template<std::size_t k, typename R=double>
auto lagrange()
{
  return [](const auto& gridView) {
    return LagrangePreBasis<std::decay_t<decltype(gridView)>, k, R>(gridView);
  };
}
 
template<typename R=double>
auto lagrange(int order)
{
  return [=](const auto& gridView) {
    return LagrangePreBasis<std::decay_t<decltype(gridView)>, -1, R>(gridView, order);
  };
}
 
} // end namespace BasisFactory
 
 
 
template<typename GV, int k=-1, typename R=double>
using LagrangeBasis = DefaultGlobalBasis<LagrangePreBasis<GV, k, R> >;
 
 
 
 
 
} // end namespace Functions
} // end namespace Dune
 
 
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEBASIS_HH