lagrangesimplex.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGESIMPLEX_HH
#define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGESIMPLEX_HH
 
#include <array>
#include <cstddef>
#include <numeric>
#include <algorithm>
#include <span>
#include <type_traits>
 
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/math.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/std/mdarray.hh>
#include <dune/common/std/mdspan.hh>
 
#include <dune/geometry/referenceelements.hh>
 
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/localfunctions/common/localkey.hh>
 
namespace Dune { namespace Impl
{
  // The traits provide static or dynamic order and size information
  template<unsigned int dim, int compileTimeOrder>
  struct LagrangeSimplexTraits
  {
    static constexpr bool is_static_order = (compileTimeOrder >= 0);
 
    constexpr LagrangeSimplexTraits(int /*runTimeOrder*/ = compileTimeOrder)
    {
      static_assert(compileTimeOrder >= 0, "LagrangeSimplex: order must be non-negative");
    }
 
    static constexpr unsigned int order() { return compileTimeOrder; }
 
    static constexpr unsigned int size() { return binomial(compileTimeOrder+dim,dim); }
  };
 
 
  // this specialization is for the dynamic order case
  template<unsigned int dim>
  struct LagrangeSimplexTraits<dim,-1>
  {
    const unsigned int k_;
    const unsigned int size_;
 
    static constexpr bool is_static_order = false;
 
    constexpr explicit LagrangeSimplexTraits(int runTimeOrder)
      : k_(runTimeOrder >= 0 ? (unsigned int)(runTimeOrder) : 0u)
      , size_(binomial(k_+dim,dim))
    {
      if (runTimeOrder < 0)
        DUNE_THROW(Dune::InvalidStateException, "LagrangeSimplex: order must be non-negative");
    }
 
    constexpr unsigned int order() const
    {
      return k_;
    }
 
    constexpr unsigned int size() const
    {
      return size_;
    }
  };
 
 
  template<class R, unsigned int dim, int k>
  struct LagrangeSimplexLocalBasisBuffers
  {
    LagrangeSimplexLocalBasisBuffers(int /*order*/) {}
 
    // Buffer for the evaluation of shape function, used in evaluateFunction()
    static auto makeValueBuffer(int /*order*/)
    {
      using E = Std::extents<int,dim+1,k+1>;
      using C = std::array<R,(dim+1)*(k+1)>;
      return Std::mdarray<R,E,Std::layout_right,C>{};
    }
 
    // Buffer for the evaluation of shape function Jacobians, used in evaluateJacobian()
    static auto makeJacobianBuffer(int /*order*/)
    {
      using E = Std::extents<int,dim+1,2,k+1>;
      using C = std::array<R,(dim+1)*2*(k+1)>;
      return Std::mdarray<R,E,Std::layout_right,C>{};
    }
 
    // Buffer for the evaluation of shape function partial derivatives, used in partial()
    template <class T, T derivativeOrder>
    static auto makePartialBuffer(std::integral_constant<T,derivativeOrder>, int /*order*/)
    {
      using E = Std::extents<int,dim+1,derivativeOrder+1,k+1>;
      using C = std::array<R,(dim+1)*(derivativeOrder+1)*(k+1)>;
      return Std::mdarray<R,E,Std::layout_right,C>{};
    }
  };
 
 
  template<class R, unsigned int dim>
  struct LagrangeSimplexLocalBasisBuffers<R,dim,-1>
  {
    LagrangeSimplexLocalBasisBuffers(int order)
      : buffer_((dim+1)*std::max<int>(2,order+1)*(order+1))
    {}
 
    mutable std::vector<R> buffer_{};
 
    // Buffer for the evaluation of shape function, used in evaluateFunction()
    auto makeValueBuffer(int order) const
    {
      using E = Std::extents<int,dim+1,std::dynamic_extent>;
      return Std::mdspan<R,E>{buffer_.data(), E(order+1)};
    }
 
    // Buffer for the evaluation of shape function Jacobians, used in evaluateJacobian()
    auto makeJacobianBuffer(int order) const
    {
      using E = Std::extents<int,dim+1,2,std::dynamic_extent>;
      return Std::mdspan<R,E>{buffer_.data(), E(order+1)};
    }
 
    // Buffer for the evaluation of shape function partial derivatives, used in partial()
    auto makePartialBuffer(int derivativeOrder, int order) const
    {
      using E = Std::extents<int,dim+1,std::dynamic_extent,std::dynamic_extent>;
      return Std::mdspan<R,E>{buffer_.data(), E(derivativeOrder+1, order+1)};
    }
  };
 
 
  template<class D, class R, unsigned int dim, int compileTimeOrder>
  class LagrangeSimplexLocalBasis
    : public LagrangeSimplexTraits<dim,compileTimeOrder>
    , private LagrangeSimplexLocalBasisBuffers<R,dim,compileTimeOrder>
  {
    template <class> friend class LagrangeSimplexLocalInterpolation;
 
  public:
 
    using Base = LagrangeSimplexTraits<dim,compileTimeOrder>;
    using Buffers = LagrangeSimplexLocalBasisBuffers<R,dim,compileTimeOrder>;
    static constexpr bool is_static_order = Base::is_static_order;
 
    constexpr LagrangeSimplexLocalBasis(int runTimeOrder = compileTimeOrder)
      : Base(runTimeOrder)
      , Buffers(runTimeOrder)
    {}
 
    using Base::order;
    using Base::size;
 
  private:
 
    // Fix the first index in a multi-dimensional array to `i`. For matrices this
    // corresponds to the ith row of the matrix.
    // Equivalent to std::submdspan(L, i, std::full_extent...)
    template <class I, std::size_t e0, std::size_t... ee, class A>
    static auto slice(Std::mdspan<R,Std::extents<I,e0,ee...>,Std::layout_right,A> L, int i)
    {
      return Dune::unpackIntegerSequence([&](auto... ii) {
        return Std::mdspan<R,Std::extents<I,ee...>,Std::layout_right>{
          L.data_handle() + i * L.stride(0), L.extent(ii+1)...};
      }, std::make_index_sequence<sizeof...(ee)>{});
    }
 
    // Overload for mdarray type: forward to the mdspan implementation
    template <class I, std::size_t e0, std::size_t... ee, class C>
    static auto slice(Std::mdarray<R,Std::extents<I,e0,ee...>,Std::layout_right,C>& L, int i)
    {
      return slice(L.to_mdspan(), i);
    }
 
    // Compute the rescaled barycentric coordinates of x.
    // We rescale the simplex by k and then compute the
    // barycentric coordinates with respect to the points
    // p_i = e_i (for i=0,...,dim-1) and p_dim=0.
    // Notice that then the Lagrange points have the barycentric
    // coordinates (i_0,...,i_d) where i_j are all non-negative
    // integers satisfying the constraint sum i_j = k.
    constexpr auto barycentric(const auto& x) const
    {
      auto b = std::array<R,dim+1>{};
      b[dim] = order();
      for (auto i : Dune::range(dim))
      {
        b[i] = order() * x[i];
        b[dim] -= b[i];
      }
      return b;
    }
 
    // Evaluate the univariate Lagrange polynomials L_i(t) for i=0,...,k where
    //
    // L_i(t) = (t-0)/(i-0) * ... * (t-(i-1))/(i-(i-1))
    //        = (t-0)*...*(t-(i-1))/(i!);
    static constexpr void evaluateLagrangePolynomials(const R& t, auto&& L, int k)
    {
      L[0] = 1;
      for (auto i : Dune::range(k))
        L[i+1] = L[i] * (t - i) / (i+1);
    }
 
    // Evaluate the univariate Lagrange polynomial derivatives L_i(t) for i=0,...,k
    // up to given maxDerivativeOrder.
    static constexpr void evaluateLagrangePolynomialDerivative(const R& t, auto&& LL, unsigned int maxDerivativeOrder, int k)
    {
      auto L = slice(LL,0);
      L[0] = 1;
      for (auto i : Dune::range(k))
        L[i+1] = L[i] * (t - i) / (i+1);
      for(auto j : Dune::range(maxDerivativeOrder))
      {
        auto F = slice(LL,j);
        auto DF = slice(LL,j+1);
        DF[0] = 0;
        for (auto i : Dune::range(k))
          DF[i+1] = (DF[i] * (t - i) + (j+1)*F[i]) / (i+1);
      }
    }
 
    using BarycentricMultiIndex = std::array<int,dim+1>;
 
 
    // This computed the required partial derivative given by the multi-index
    // beta of a product of a function given as a product of dim+1 derivatives
    // of univariate Lagrange polynomials of the dim+1 barycentric coordinates.
    // The polynomials in the product are specified as follows:
    //
    // The table L contains all required derivatives of all univariate
    // polynomials evaluated at all barycentric coordinates.  The two
    // multi-indices i and alpha encode that the polynomial for the
    // j-th barycentric coordinate is the alpha-j-th derivative of
    // the i_j-the Lagrange polynomial.
    //
    // Hence this method computes D_beta f(x) where f(x) is the product
    // \f$f(x) = \prod_{j=0}^{d} L_{i_j}^{(alpha_j)}(x_j) \f$.
    static constexpr R barycentricDerivative(
        BarycentricMultiIndex beta,
        const auto&L,
        int k,
        const BarycentricMultiIndex& i,
        const BarycentricMultiIndex& alpha = {})
    {
      // If there are unprocessed derivatives left we search the first unprocessed
      // partial derivative direction j and compute it using the product and chain rule.
      // The remaining partial derivatives are forwarded to the recursion.
      // Notice that the partial derivative of the last barycentric component
      // wrt to the j-th component is -1 which is responsible for the second
      // term in the sum. Furthermore we get the factor k due to the scaling of
      // the simplex.
      for(auto j : Dune::range(dim))
      {
        if (beta[j] > 0)
        {
          auto leftDerivatives = alpha;
          auto rightDerivatives = alpha;
          leftDerivatives[j]++;
          rightDerivatives.back()++;
          beta[j]--;
          return (barycentricDerivative(beta, L, k, i, leftDerivatives) -
                  barycentricDerivative(beta, L, k, i, rightDerivatives)) * k;
        }
      }
      // If there are no unprocessed derivatives we can simply evaluate
      // the product of the derivatives of the Lagrange polynomials with
      // given indices i_j and derivative orders alpha_j
      // Evaluate the product of the univariate Lagrange polynomials
      // with given indices and orders.
      auto y = R(1);
      for(auto j : Dune::range(dim+1))
        y *= L(j,alpha[j],i[j]);
      return y;
    }
 
 
  public:
    using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
 
    void evaluateFunction(const typename Traits::DomainType& x,
                          std::vector<typename Traits::RangeType>& out) const
    {
      const int k = order();
      out.resize(size());
 
      // Specialization for zero-order case
      if (k==0)
      {
        out[0] = 1;
        return;
      }
 
      // Specialization for first-order case
      if (k==1)
      {
        out[0] = 1.0;
        for (size_t i=0; i<dim; i++)
        {
          out[0]  -= x[i];
          out[i+1] = x[i];
        }
        return;
      }
 
      // Compute rescaled barycentric coordinates of x
      auto z = barycentric(x);
 
      auto L = Buffers::makeValueBuffer(k);
      for (auto j : Dune::range(dim+1))
        evaluateLagrangePolynomials(z[j], slice(L,j), k);
 
      if (dim==1)
      {
        unsigned int n = 0;
        for (auto i0 : Dune::range(k + 1))
          for (auto i1 : std::array{k - i0})
            out[n++] = L(0,i0) * L(1,i1);
        return;
      }
      if (dim==2)
      {
        unsigned int n=0;
        for (auto i1 : Dune::range(k + 1))
          for (auto i0 : Dune::range(k - i1 + 1))
            for (auto i2 : std::array{k - i1 - i0})
              out[n++] = L(0,i0) * L(1,i1) * L(2,i2);
        return;
      }
      if (dim==3)
      {
        unsigned int n = 0;
        for (auto i2 : Dune::range(k + 1))
          for (auto i1 : Dune::range(k - i2 + 1))
            for (auto i0 : Dune::range(k - i2 - i1 + 1))
              for (auto i3 : std::array{k - i2 - i1 - i0})
                out[n++] = L(0,i0) * L(1,i1)  * L(2,i2) * L(3,i3);
        return;
      }
 
      DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis for k>=2 only implemented for dim<=3");
    }
 
    void evaluateJacobian(const typename Traits::DomainType& x,
                          std::vector<typename Traits::JacobianType>& out) const
    {
      const int k = order();
      out.resize(size());
 
      // Specialization for k==0
      if (k==0)
      {
        std::fill(out[0][0].begin(), out[0][0].end(), 0);
        return;
      }
 
      // Specialization for k==1
      if (k==1)
      {
        std::fill(out[0][0].begin(), out[0][0].end(), -1);
 
        for (unsigned int i=0; i<dim; i++)
          for (unsigned int j=0; j<dim; j++)
            out[i+1][0][j] = (i==j);
 
        return;
      }
 
      // Compute rescaled barycentric coordinates of x
      auto z = barycentric(x);
 
      // L(j,m,i) is the m-th derivative of the i-th Lagrange polynomial at z[j]
      auto L = Buffers::makeJacobianBuffer(k);
      for (auto j : Dune::range(dim+1))
        evaluateLagrangePolynomialDerivative(z[j], slice(L,j), 1, k);
 
      if (dim==1)
      {
        unsigned int n = 0;
        for (auto i0 : Dune::range(k + 1))
        {
          for (auto i1 : std::array{k-i0})
          {
            out[n][0][0] = (L(0,1,i0) * L(1,0,i1) - L(0,0,i0) * L(1,1,i1))*k;
            ++n;
          }
        }
        return;
      }
      if (dim==2)
      {
        unsigned int n=0;
        for (auto i1 : Dune::range(k + 1))
        {
          for (auto i0 : Dune::range(k - i1 + 1))
          {
            for (auto i2 : std::array{k - i1 - i0})
            {
              out[n][0][0] = (L(0,1,i0) * L(1,0,i1) * L(2,0,i2) - L(0,0,i0) * L(1,0,i1) * L(2,1,i2))*k;
              out[n][0][1] = (L(0,0,i0) * L(1,1,i1) * L(2,0,i2) - L(0,0,i0) * L(1,0,i1) * L(2,1,i2))*k;
              ++n;
            }
          }
        }
        return;
      }
      if (dim==3)
      {
        unsigned int n = 0;
        for (auto i2 : Dune::range(k + 1))
        {
          for (auto i1 : Dune::range(k - i2 + 1))
          {
            for (auto i0 : Dune::range(k - i2 - i1 + 1))
            {
              for (auto i3 : std::array{k - i2 - i1 - i0})
              {
                out[n][0][0] = (L(0,1,i0) * L(1,0,i1) * L(2,0,i2) * L(3,0,i3) - L(0,0,i0) * L(1,0,i1) * L(2,0,i2) * L(3,1,i3))*k;
                out[n][0][1] = (L(0,0,i0) * L(1,1,i1) * L(2,0,i2) * L(3,0,i3) - L(0,0,i0) * L(1,0,i1) * L(2,0,i2) * L(3,1,i3))*k;
                out[n][0][2] = (L(0,0,i0) * L(1,0,i1) * L(2,1,i2) * L(3,0,i3) - L(0,0,i0) * L(1,0,i1) * L(2,0,i2) * L(3,1,i3))*k;
                ++n;
              }
            }
          }
        }
 
        return;
      }
 
      DUNE_THROW(NotImplemented, "LagrangeSimplexLocalBasis for k>=2 only implemented for dim<=3");
    }
 
    void partial(const std::array<unsigned int,dim>& derivativeOrder,
                 const typename Traits::DomainType& in,
                 std::vector<typename Traits::RangeType>& out) const
    {
      const int k = order();
      int totalOrder = std::accumulate(derivativeOrder.begin(), derivativeOrder.end(), 0u);
 
      out.resize(size());
 
      // Derivative order zero corresponds to the function evaluation.
      if (totalOrder == 0)
      {
        evaluateFunction(in, out);
        return;
      }
 
      // Derivatives of order >k are all zero.
      if (totalOrder > k)
      {
        for(auto& out_i : out)
          out_i = 0;
        return;
      }
 
      // It remains to cover the cases 0 < totalOrder<= k.
 
      if (k==1)
      {
        if (totalOrder==1)
        {
          auto direction = std::find(derivativeOrder.begin(), derivativeOrder.end(), 1);
          out[0] = -1;
          for (unsigned int i=0; i<dim; i++)
            out[i+1] = (i==(direction-derivativeOrder.begin()));
        }
        return;
      }
 
      // Since the required stack storage depends on the dynamic total order,
      // we need to do a dynamic to static dispatch by enumerating all supported
      // static orders.
      auto supportedOrders = [&]{
        if constexpr(is_static_order)
          return Dune::range(Dune::index_constant<1>{}, Dune::index_constant<compileTimeOrder+1>{});
        else
          return Dune::range(1, k+1);
      }();
      return Dune::Hybrid::switchCases(supportedOrders, totalOrder, [&](auto staticTotalOrder) {
 
        // Compute rescaled barycentric coordinates of x
        auto z = barycentric(in);
 
        // L(j,m,i) is the m-th derivative of the i-th Lagrange polynomial at z[j]
        auto L = Buffers::makePartialBuffer(staticTotalOrder,k);
        for (auto j : Dune::range(dim))
          evaluateLagrangePolynomialDerivative(z[j], slice(L,j), derivativeOrder[j], k);
        evaluateLagrangePolynomialDerivative(z[dim], slice(L,dim), totalOrder, k);
 
        auto barycentricOrder = BarycentricMultiIndex{};
        for (auto j : Dune::range(dim))
          barycentricOrder[j] = derivativeOrder[j];
        barycentricOrder[dim] = 0;
 
        if constexpr (dim==1)
        {
          unsigned int n = 0;
          for (auto i0 : Dune::range(k + 1))
            for (auto i1 : std::array{k - i0})
              out[n++] = barycentricDerivative(barycentricOrder, L, k, BarycentricMultiIndex{i0, i1});
        }
        if constexpr (dim==2)
        {
          unsigned int n=0;
          for (auto i1 : Dune::range(k + 1))
            for (auto i0 : Dune::range(k - i1 + 1))
              for (auto i2 : std::array{k - i1 - i0})
                out[n++] = barycentricDerivative(barycentricOrder, L, k, BarycentricMultiIndex{i0, i1, i2});
        }
        if constexpr (dim==3)
        {
          unsigned int n = 0;
          for (auto i2 : Dune::range(k + 1))
            for (auto i1 : Dune::range(k - i2 + 1))
              for (auto i0 : Dune::range(k - i2 - i1 + 1))
                for (auto i3 : std::array{k - i2 - i1 - i0})
                  out[n++] = barycentricDerivative(barycentricOrder, L, k, BarycentricMultiIndex{i0, i1, i2, i3});
        }
      });
    }
  };
 
  template<unsigned int dim, int compileTimeOrder>
  class LagrangeSimplexLocalCoefficients
    : public LagrangeSimplexTraits<dim,compileTimeOrder>
  {
    using Base = LagrangeSimplexTraits<dim,compileTimeOrder>;
 
  public:
    LagrangeSimplexLocalCoefficients (int runTimeOrder = compileTimeOrder)
      : Base(runTimeOrder)
      , localKeys_(Base::size())
    {
      const int k = Base::order();
      if (k==0)
      {
        localKeys_[0] = LocalKey(0,0,0);
        return;
      }
 
      if (k==1)
      {
        for (std::size_t i=0; i<Base::size(); i++)
          localKeys_[i] = LocalKey(i,dim,0);
        return;
      }
 
      if (dim==1)
      {
        // Order is at least 2 here
        localKeys_[0] = LocalKey(0,1,0);          // vertex dof
        for (int i=1; i<k; i++)
          localKeys_[i] = LocalKey(0,0,i-1);      // element dofs
        localKeys_.back() = LocalKey(1,1,0);      // vertex dof
        return;
      }
 
      if (dim==2)
      {
        int n=0;
        int c=0;
        for (int j=0; j<=k; j++)
          for (int i=0; i<=k-j; i++)
          {
            if (i==0 && j==0)
            {
              localKeys_[n++] = LocalKey(0,2,0);
              continue;
            }
            if (i==k && j==0)
            {
              localKeys_[n++] = LocalKey(1,2,0);
              continue;
            }
            if (i==0 && j==k)
            {
              localKeys_[n++] = LocalKey(2,2,0);
              continue;
            }
            if (j==0)
            {
              localKeys_[n++] = LocalKey(0,1,i-1);
              continue;
            }
            if (i==0)
            {
              localKeys_[n++] = LocalKey(1,1,j-1);
              continue;
            }
            if (i+j==k)
            {
              localKeys_[n++] = LocalKey(2,1,j-1);
              continue;
            }
            localKeys_[n++] = LocalKey(0,0,c++);
          }
        return;
      }
 
      if (dim==3)
      {
        std::array<unsigned int, dim+1> vertexMap;
        for (unsigned int i=0; i<=dim; i++)
          vertexMap[i] = i;
        generateLocalKeys(vertexMap);
        return;
      }
      DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for k<=1 or dim<=3!");
    }
 
    constexpr LagrangeSimplexLocalCoefficients (const std::array<unsigned int, dim+1> vertexMap)
      : localKeys_(Base::size())
    {
      if (dim!=2 && dim!=3)
        DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for dim==2 and dim==3!");
 
      generateLocalKeys(vertexMap);
    }
 
 
    template<class VertexMap>
    constexpr LagrangeSimplexLocalCoefficients(const VertexMap &vertexmap)
      : localKeys_(Base::size())
    {
      if (dim!=2 && dim!=3)
        DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for dim==2 and dim==3!");
 
      std::array<unsigned int, dim+1> vertexmap_array;
      std::copy(vertexmap, vertexmap + dim + 1, vertexmap_array.begin());
      generateLocalKeys(vertexmap_array);
    }
 
    constexpr const LocalKey& localKey (std::size_t i) const
    {
      return localKeys_[i];
    }
 
    using Base::order;
    using Base::size;
 
  private:
    std::vector<LocalKey> localKeys_;
 
    constexpr void generateLocalKeys(const std::array<unsigned int, dim+1> vertexMap)
    {
      const int k = Base::order();
      if (k==0)
      {
        localKeys_[0] = LocalKey(0,0,0);
        return;
      }
 
      if (dim==2)
      {
        // Create default assignment
        int n=0;
        int c=0;
        for (int j=0; j<=k; j++)
          for (int i=0; i<=k-j; i++)
          {
            if (i==0 && j==0)
            {
              localKeys_[n++] = LocalKey(0,2,0);
              continue;
            }
            if (i==k && j==0)
            {
              localKeys_[n++] = LocalKey(1,2,0);
              continue;
            }
            if (i==0 && j==k)
            {
              localKeys_[n++] = LocalKey(2,2,0);
              continue;
            }
            if (j==0)
            {
              localKeys_[n++] = LocalKey(0,1,i-1);
              continue;
            }
            if (i==0)
            {
              localKeys_[n++] = LocalKey(1,1,j-1);
              continue;
            }
            if (i+j==k)
            {
              localKeys_[n++] = LocalKey(2,1,j-1);
              continue;
            }
            localKeys_[n++] = LocalKey(0,0,c++);
          }
 
        // Flip edge orientations, if requested
        bool flip[3];
        flip[0] = vertexMap[0] > vertexMap[1];
        flip[1] = vertexMap[0] > vertexMap[2];
        flip[2] = vertexMap[1] > vertexMap[2];
        for (std::size_t i=0; i<Base::size(); i++)
          if (localKeys_[i].codim()==1 && flip[localKeys_[i].subEntity()])
            localKeys_[i].index(k-2-localKeys_[i].index());
 
        return;
      }
 
      if (dim!=3)
        DUNE_THROW(NotImplemented, "LagrangeSimplexLocalCoefficients only implemented for dim==3!");
 
      unsigned int subindex[16];
      unsigned int codim_count[4] = {0};
      for (unsigned int m = 1; m < 16; ++m)
      {
        unsigned int codim = !(m&1) + !(m&2) + !(m&4) + !(m&8);
        subindex[m] = codim_count[codim]++;
      }
 
      int a1 = (3*k + 12)*k + 11;
      int a2 = -3*k - 6;
      unsigned int dof_count[16] = {0};
      int i[4];
      for (i[3] = 0; i[3] <= k; ++i[3])
        for (i[2] = 0; i[2] <= k - i[3]; ++i[2])
          for (i[1] = 0; i[1] <= k - i[2] - i[3]; ++i[1])
          {
            i[0] = k - i[1] - i[2] - i[3];
            int j[4];
            unsigned int entity = 0;
            unsigned int codim = 0;
            for (unsigned int m = 0; m < 4; ++m)
            {
              j[m] = i[vertexMap[m]];
              entity += !!j[m] << m;
              codim += !j[m];
            }
            int local_index = j[3]*(a1 + (a2 + j[3])*j[3])/6
                              + j[2]*(2*(k - j[3]) + 3 - j[2])/2 + j[1];
            localKeys_[local_index] = LocalKey(subindex[entity], codim, dof_count[entity]++);
          }
    }
  };
 
  template<class LocalBasis>
  class LagrangeSimplexLocalInterpolation
    : public LocalBasis::Base
  {
    using Base = typename LocalBasis::Base;
 
  public:
    constexpr LagrangeSimplexLocalInterpolation () = default;
 
    explicit constexpr LagrangeSimplexLocalInterpolation (const LocalBasis& localBasis)
      : Base(localBasis)
    {}
 
    template<typename F, typename C>
    constexpr void interpolate (const F& f, std::vector<C>& out) const
    {
      constexpr auto dim = LocalBasis::Traits::dimDomain;
      const int k = Base::order();
      using D = typename LocalBasis::Traits::DomainFieldType;
 
      typename LocalBasis::Traits::DomainType x;
 
      out.resize(Base::size());
 
      // Specialization for zero-order case
      if (k==0)
      {
        auto center = ReferenceElements<D,dim>::simplex().position(0,0);
        out[0] = f(center);
        return;
      }
 
      // Specialization for first-order case
      if (k==1)
      {
        // vertex 0
        std::fill(x.begin(), x.end(), 0);
        out[0] = f(x);
 
        // remaining vertices
        for (int i=0; i<dim; i++)
        {
          for (int j=0; j<dim; j++)
            x[j] = (i==j);
 
          out[i+1] = f(x);
        }
        return;
      }
 
      if (dim==1)
      {
        for (int i=0; i<k+1; i++)
        {
          x[0] = ((D)i)/k;
          out[i] = f(x);
        }
        return;
      }
 
      if (dim==2)
      {
        int n=0;
        for (int j=0; j<=k; j++)
          for (int i=0; i<=k-j; i++)
          {
            x = { ((D)i)/k, ((D)j)/k };
            out[n] = f(x);
            n++;
          }
        return;
      }
 
      if (dim!=3)
        DUNE_THROW(NotImplemented, "LagrangeSimplexLocalInterpolation only implemented for dim<=3!");
 
      const int kdiv = (k==0 ? 1 : k);
 
      int n=0;
      for (int i2 = 0; i2 <= k; i2++)
        for (int i1 = 0; i1 <= k-i2; i1++)
          for (int i0 = 0; i0 <= k-i1-i2; i0++)
          {
            x[0] = ((D)i0)/((D)kdiv);
            x[1] = ((D)i1)/((D)kdiv);
            x[2] = ((D)i2)/((D)kdiv);
            out[n] = f(x);
            n++;
          }
    }
 
  };
 
} }    // namespace Dune::Impl
 
namespace Dune
{
  template<class D, class R, int d, int compileTimeOrder = -1>
  class LagrangeSimplexLocalFiniteElement
  {
  public:
    using Traits = LocalFiniteElementTraits<
      Impl::LagrangeSimplexLocalBasis<D,R,d,compileTimeOrder>,
      Impl::LagrangeSimplexLocalCoefficients<d,compileTimeOrder>,
      Impl::LagrangeSimplexLocalInterpolation<Impl::LagrangeSimplexLocalBasis<D,R,d,compileTimeOrder> > >;
 
    constexpr LagrangeSimplexLocalFiniteElement()
      : basis_()
      , coefficients_()
      , interpolation_(basis_)
    {
      static_assert(compileTimeOrder >= 0, "Default constructor only allowed for compile-time >= 0");
    }
 
    explicit constexpr LagrangeSimplexLocalFiniteElement(int runTimeOrder)
      : basis_(runTimeOrder)
      , coefficients_(runTimeOrder)
      , interpolation_(basis_)
    {
      if (runTimeOrder < 0)
        DUNE_THROW(Dune::InvalidStateException, "LagrangeSimplexLocalFiniteElement: run-time order must be non-negative!");
    }
 
    template<typename VertexMap>
    explicit constexpr LagrangeSimplexLocalFiniteElement(const VertexMap& vertexmap) requires(compileTimeOrder >= 0)
      : basis_()
      , coefficients_(vertexmap)
      , interpolation_(basis_)
    {}
 
    constexpr const typename Traits::LocalBasisType& localBasis () const
    {
      return basis_;
    }
 
    constexpr const typename Traits::LocalCoefficientsType& localCoefficients () const
    {
      return coefficients_;
    }
 
    constexpr const typename Traits::LocalInterpolationType& localInterpolation () const
    {
      return interpolation_;
    }
 
    constexpr std::size_t size () const
    {
      return basis_.size();
    }
 
    static constexpr GeometryType type ()
    {
      return GeometryTypes::simplex(d);
    }
 
  private:
    typename Traits::LocalBasisType basis_;
    typename Traits::LocalCoefficientsType coefficients_;
    typename Traits::LocalInterpolationType interpolation_;
  };
 
}        // namespace Dune
 
#endif   // DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGESIMPLEX_HH