localprojectedfefunction.hh

#ifndef DUNE_GFE_FUNCTIONS_LOCALPROJECTEDFEFUNCTION_HH
#define DUNE_GFE_FUNCTIONS_LOCALPROJECTEDFEFUNCTION_HH

#include <vector>

#include <dune/common/fvector.hh>

#include <dune/geometry/type.hh>

#include <dune/gfe/linearalgebra.hh>
#include <dune/gfe/polardecomposition.hh>
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/productmanifold.hh>
#include <dune/gfe/spaces/rotation.hh>

namespace Dune::GFE
{

  /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold
   *
   * \tparam dim Dimension of the reference element
   * \tparam ctype Type used for coordinates on the reference element
   * \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
   * \tparam TargetSpace The manifold that the function takes its values in
   * \tparam conforming Geometrical conformity of the functions (omits projections when set to false)
   */
  template <int dim, class ctype, class LocalFiniteElement, class TS, bool conforming=true>
  class LocalProjectedFEFunction
  {
  public:
    using TargetSpace=TS;
  private:
    typedef typename TargetSpace::ctype RT;

    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
    static const int embeddedDim = EmbeddedTangentVector::dimension;

  public:

    /** \brief The type used for derivatives */
    typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;

    /** \brief Constructor
     * \param localFiniteElement A Lagrangian finite element that provides the interpolation points
     * \param coefficients Values of the function at the Lagrange points
     */
    LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
                             const std::vector<TargetSpace>& coefficients)
      : localFiniteElement_(localFiniteElement),
      coefficients_(coefficients)
    {
      assert(localFiniteElement_.localBasis().size() == coefficients_.size());
    }

    /** \brief Rebind the FEFunction to another TargetSpace */
    template<class U>
    struct rebind
    {
      using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U,conforming>;
    };

    /** \brief The number of Lagrange points */
    unsigned int size() const
    {
      return localFiniteElement_.localBasis().size();
    }

    /** \brief The type of the reference element */
    Dune::GeometryType type() const
    {
      return localFiniteElement_.type();
    }

    /** \brief The scalar finite element used as the interpolation weights
     *
     * \note This method was added for InterpolationDerivatives, which needs it
     * to construct a copy of a LocalGeodesicFEFunction with ADOL-C's adouble
     * number type.  This is not optimal, because the localFiniteElement
     * really is an implementation detail of LocalGeodesicFEFunction and
     * should not be needed just to copy an entire object.  Other non-Euclidean
     * interpolation rules may not have such a finite element at all.
     * Therefore, this method may disappear again eventually.
     */
    const LocalFiniteElement& localFiniteElement() const
    {
      return localFiniteElement_;
    }

    /** \brief Evaluate the function */
    auto evaluate(const Dune::FieldVector<ctype, dim>& local) const;

    /** \brief Evaluate the derivative of the function */
    DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const;

    /** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
     *        \param local Local coordinates in the reference element where to evaluate the derivative
     *        \param q Value of the local gfe function at 'local'.  If you provide something wrong here the result will be wrong, too!
     *
     * \note This method is only usable in the conforming setting, because it requires the caller
     * to hand over the interpolation value as a TargetSpace object.
     */
    DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
                                      const TargetSpace& q) const;

    /** \brief Evaluate the value and the derivative of the interpolation function
     *
     * \return A std::pair containing the value and the first derivative of the interpolation function.
     * If the interpolation is conforming then the first member of the pair will be a TargetSpace.
     * Otherwise it will be a RealTuple.
     */
    auto evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const;

    /** \brief Get the i'th base coefficient. */
    const TargetSpace& coefficient(int i) const
    {
      return coefficients_[i];
    }
  private:

    /** \brief The scalar local finite element, which provides the weighting factors
     *        \todo We really only need the local basis
     */
    const LocalFiniteElement& localFiniteElement_;

    /** \brief The coefficient vector */
    std::vector<TargetSpace> coefficients_;

  };

  template <int dim, class ctype, class LocalFiniteElement, class TargetSpace, bool conforming>
  auto LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
  evaluate(const Dune::FieldVector<ctype, dim>& local) const
  {
    // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
    std::vector<Dune::FieldVector<ctype,1> > w;
    localFiniteElement_.localBasis().evaluateFunction(local,w);

    typename TargetSpace::CoordinateType c(0);
    for (size_t i=0; i<coefficients_.size(); i++)
      c.axpy(w[i][0], coefficients_[i].globalCoordinates());

    if constexpr (conforming)
      return TargetSpace::projectOnto(c);
    else
      return (RealTuple<RT, TargetSpace::CoordinateType::dimension>)c;

  }

  template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
  typename LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::DerivativeType
  LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
  evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
  {
    if constexpr(conforming)
    {
      // the function value at the point where we are evaluating the derivative
      TargetSpace q = evaluate(local);

      // Actually compute the derivative
      return evaluateDerivative(local, q);
    }
    else {
      std::vector<Dune::FieldMatrix<ctype, 1, dim> > wDer;
      localFiniteElement_.localBasis().evaluateJacobian(local, wDer);

      Dune::FieldMatrix<RT, embeddedDim, dim> derivative(0);
      for (size_t i = 0; i < embeddedDim; i++)
        for (size_t j = 0; j < dim; j++)
          for (size_t k = 0; k < coefficients_.size(); k++)
            derivative[i][j] += wDer[k][0][j] * coefficients_[k].globalCoordinates()[i];

      return derivative;
    }
  }

  template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
  typename LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::DerivativeType
  LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
  evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace& q) const
  {
    // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
    std::vector<Dune::FieldVector<ctype,1> > w;
    localFiniteElement_.localBasis().evaluateFunction(local,w);

    std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
    localFiniteElement_.localBasis().evaluateJacobian(local,wDer);

    typename TargetSpace::CoordinateType embeddedInterpolation(0);
    for (size_t i=0; i<coefficients_.size(); i++)
      embeddedInterpolation.axpy(w[i][0], coefficients_[i].globalCoordinates());

    Dune::FieldMatrix<RT,embeddedDim,dim> derivative(0);
    for (size_t i=0; i<embeddedDim; i++)
      for (size_t j=0; j<dim; j++)
        for (size_t k=0; k<coefficients_.size(); k++)
          derivative[i][j] += wDer[k][0][j] * coefficients_[k].globalCoordinates()[i];

    auto derivativeOfProjection = TargetSpace::derivativeOfProjection(embeddedInterpolation);

    return derivativeOfProjection*derivative;
  }

  template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
  auto
  LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
  evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const
  {
    // Construct the type of the result -- it depends on whether the interpolation
    // is conforming or not.
    using Value = std::conditional_t<conforming,TargetSpace,RealTuple<RT,embeddedDim> >;

    std::pair<Value,DerivativeType> result;

    ///////////////////////////////////////////////////////////
    //  Compute the value of the interpolation function
    ///////////////////////////////////////////////////////////

    std::vector<Dune::FieldVector<ctype,1> > w;
    localFiniteElement_.localBasis().evaluateFunction(local,w);

    typename TargetSpace::CoordinateType embeddedInterpolation(0);
    for (size_t i=0; i<coefficients_.size(); i++)
      embeddedInterpolation.axpy(w[i][0], coefficients_[i].globalCoordinates());

    if constexpr (conforming)
      result.first = TargetSpace::projectOnto(embeddedInterpolation);
    else
      result.first = (RealTuple<RT, TargetSpace::CoordinateType::dimension>)embeddedInterpolation;

    ///////////////////////////////////////////////////////////
    //  Compute the derivative of the interpolation function
    ///////////////////////////////////////////////////////////

    // Compute the interpolation in the surrounding space
    std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
    localFiniteElement_.localBasis().evaluateJacobian(local,wDer);

    result.second = 0.0;
    for (size_t i=0; i<embeddedDim; i++)
      for (size_t j=0; j<dim; j++)
        for (size_t k=0; k<coefficients_.size(); k++)
          result.second[i][j] += wDer[k][0][j] * coefficients_[k].globalCoordinates()[i];

    // The derivative of the projection onto the manifold
    if constexpr(conforming)
      result.second = TargetSpace::derivativeOfProjection(embeddedInterpolation) * result.second;

    return result;
  }


  /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for SO(3)
   *
   * \tparam dim Dimension of the reference element
   * \tparam ctype Type used for coordinates on the reference element
   * \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
   */
  template <int dim, class ctype, class LocalFiniteElement, class field_type, bool conforming>
  class LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3>,conforming>
  {
  public:
    typedef Rotation<field_type,3> TargetSpace;
  private:
    typedef typename TargetSpace::ctype RT;

    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
    static const int embeddedDim = EmbeddedTangentVector::dimension;

    /**
     * \param A The argument of the projection
     * \param polar The image of the projection, i.e., the polar factor of A
     */
    static std::array<std::array<FieldMatrix<field_type,3,3>, 3>, 3> derivativeOfProjection(const FieldMatrix<field_type,3,3>& A,
                                                                                            FieldMatrix<field_type,3,3>& polar)
    {
      std::array<std::array<FieldMatrix<field_type,3,3>, 3>, 3> result;

      for (int i=0; i<3; i++)
        for (int j=0; j<3; j++)
          for (int k=0; k<3; k++)
            for (int l=0; l<3; l++)
              result[i][j][k][l] = (i==k) and (j==l);

      polar = A;

      // Use Heron's method
      const size_t maxIterations = 100;
      double tol = 0.001;
      for (size_t iteration=0; iteration<maxIterations; iteration++)
      {
        auto oldPolar = polar;
        auto polarInvert = polar;
        polarInvert.invert();
        for (size_t i=0; i<polar.N(); i++)
          for (size_t j=0; j<polar.M(); j++)
            polar[i][j] = 0.5 * (polar[i][j] + polarInvert[j][i]);

        // Alternative name to align the code better with a description in a math text
        const auto& dQT = result;

        // Multiply from the right with Q^{-T}
        decltype(result) tmp2;
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++)
                tmp2[i][j][k][l] = 0.0;

        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++)
                for (int m=0; m<3; m++)
                  for (int o=0; o<3; o++)
                    tmp2[i][j][k][l] += polarInvert[m][i] * dQT[o][m][k][l] * polarInvert[j][o];

        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++)
                result[i][j][k][l] = 0.5 * (result[i][j][k][l] - tmp2[i][j][k][l]);

        oldPolar -= polar;
        if (oldPolar.frobenius_norm() < tol) {
          break;
        }
      }

      return result;
    }


  public:

    /** \brief The type used for derivatives */
    typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;

    /** \brief Constructor
     * \param localFiniteElement A Lagrangian finite element that provides the interpolation points
     * \param coefficients Values of the function at the Lagrange points
     */
    LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
                             const std::vector<TargetSpace>& coefficients)
      : localFiniteElement_(localFiniteElement),
      coefficients_(coefficients)
    {
      assert(localFiniteElement_.localBasis().size() == coefficients_.size());
    }

    /** \brief Rebind the FEFunction to another TargetSpace */
    template<class U>
    struct rebind
    {
      using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U>;
    };

    /** \brief The number of Lagrange points */
    unsigned int size() const
    {
      return localFiniteElement_.size();
    }

    /** \brief The type of the reference element */
    Dune::GeometryType type() const
    {
      return localFiniteElement_.type();
    }

    /** \brief The scalar finite element used as the interpolation weights
     *
     * \note This method was added for InterpolationDerivatives, which needs it
     * to construct a copy of a LocalGeodesicFEFunction with ADOL-C's adouble
     * number type.  This is not optimal, because the localFiniteElement
     * really is an implementation detail of LocalGeodesicFEFunction and
     * should not be needed just to copy an entire object.  Other non-Euclidean
     * interpolation rules may not have such a finite element at all.
     * Therefore, this method may disappear again eventually.
     */
    const LocalFiniteElement& localFiniteElement() const
    {
      return localFiniteElement_;
    }

    /** \brief Evaluate the function */
    TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const
    {
      Rotation<field_type,3> result;

      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
      std::vector<Dune::FieldVector<ctype,1> > w;
      localFiniteElement_.localBasis().evaluateFunction(local,w);

      // Interpolate in R^{3x3}
      FieldMatrix<field_type,3,3> interpolatedMatrix(0);
      for (size_t i=0; i<coefficients_.size(); i++)
      {
        FieldMatrix<field_type,3,3> coefficientAsMatrix;
        coefficients_[i].matrix(coefficientAsMatrix);
        interpolatedMatrix.axpy(w[i][0], coefficientAsMatrix);
      }

      // Project back onto SO(3)
      result.set(Dune::GFE::PolarDecomposition()(interpolatedMatrix));

      return result;
    }

    /** \brief Evaluate the derivative of the function */
    DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
    {
      // the function value at the point where we are evaluating the derivative
      TargetSpace q = evaluate(local);

      // Actually compute the derivative
      return evaluateDerivative(local,q);
    }

    /** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
     *        \param local Local coordinates in the reference element where to evaluate the derivative
     *        \param q Value of the local function at 'local'.  If you provide something wrong here the result will be wrong, too!
     */
    DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
                                      const TargetSpace& q) const
    {
      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
      std::vector<Dune::FieldVector<ctype,1> > w;
      localFiniteElement_.localBasis().evaluateFunction(local,w);

      std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
      localFiniteElement_.localBasis().evaluateJacobian(local,wDer);

      // Compute matrix representations for all coefficients (we only have them in quaternion representation)
      std::vector<Dune::FieldMatrix<field_type,3,3> > coefficientsAsMatrix(coefficients_.size());
      for (size_t i=0; i<coefficients_.size(); i++)
        coefficients_[i].matrix(coefficientsAsMatrix[i]);

      // Interpolate in R^{3x3}
      FieldMatrix<field_type,3,3> interpolatedMatrix(0);
      for (size_t i=0; i<coefficients_.size(); i++)
        interpolatedMatrix.axpy(w[i][0], coefficientsAsMatrix[i]);

      Tensor3<RT,dim,3,3> derivative(0);

      for (size_t dir=0; dir<dim; dir++)
        for (size_t i=0; i<3; i++)
          for (size_t j=0; j<3; j++)
            for (size_t k=0; k<coefficients_.size(); k++)
              derivative[dir][i][j] += wDer[k][0][dir] * coefficientsAsMatrix[k][i][j];

      FieldMatrix<field_type,3,3> polarFactor;
      auto derivativeOfProjection = this->derivativeOfProjection(interpolatedMatrix,polarFactor);

      Tensor3<field_type,dim,3,3> intermediateResult(0);

      for (size_t dir=0; dir<dim; dir++)
        for (size_t i=0; i<3; i++)
          for (size_t j=0; j<3; j++)
            for (size_t k=0; k<3; k++)
              for (size_t l=0; l<3; l++)
                intermediateResult[dir][i][j] += derivativeOfProjection[i][j][k][l]*derivative[dir][k][l];

      // One more application of the chain rule: we need to go from orthogonal matrices to quaternions
      Tensor3<field_type,4,3,3> derivativeOfMatrixToQuaternion = Rotation<field_type,3>::derivativeOfMatrixToQuaternion(polarFactor);

      DerivativeType result(0);

      for (size_t dir0=0; dir0<4; dir0++)
        for (size_t dir1=0; dir1<dim; dir1++)
          for (size_t i=0; i<3; i++)
            for (size_t j=0; j<3; j++)
              result[dir0][dir1] += derivativeOfMatrixToQuaternion[dir0][i][j] * intermediateResult[dir1][i][j];

      return result;
    }

    /** \brief Evaluate the value and the derivative of the interpolation function
     *
     * \return A std::pair containing the value and the first derivative of the interpolation function.
     * If the interpolation is conforming then the first member of the pair will be a TargetSpace.
     * Otherwise it will be a RealTuple.
     */
    auto evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const
    {
      // Construct the type of the result -- it depends on whether the interpolation
      // is conforming or not.
      using Value = std::conditional_t<conforming,TargetSpace,RealTuple<RT,embeddedDim> >;

      std::pair<Value,DerivativeType> result;

      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
      std::vector<Dune::FieldVector<ctype,1> > w;
      localFiniteElement_.localBasis().evaluateFunction(local,w);

      // Interpolate in R^{3x3}
      FieldMatrix<field_type,3,3> interpolatedMatrix(0);
      for (size_t i=0; i<coefficients_.size(); i++)
      {
        FieldMatrix<field_type,3,3> coefficientAsMatrix;
        coefficients_[i].matrix(coefficientAsMatrix);
        interpolatedMatrix.axpy(w[i][0], coefficientAsMatrix);
      }

      // Project back onto SO(3)
      static_assert(conforming, "Nonconforming interpolation into SO(3) is not implemented!");
      result.first.set(Dune::GFE::PolarDecomposition()(interpolatedMatrix));

      result.second = evaluateDerivative(local, result.first);

      return result;
    }

    /** \brief Get the i'th base coefficient. */
    const TargetSpace& coefficient(int i) const
    {
      return coefficients_[i];
    }
  private:

    /** \brief The scalar local finite element, which provides the weighting factors
     *        \todo We really only need the local basis
     */
    const LocalFiniteElement& localFiniteElement_;

    /** \brief The coefficient vector */
    std::vector<TargetSpace> coefficients_;

  };


  /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for R^3 x SO(3)
   *
   * \tparam dim Dimension of the reference element
   * \tparam ctype Type used for coordinates on the reference element
   * \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
   */
  template <int dim, class ctype, class LocalFiniteElement, class field_type>
  class LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
  {
  public:
    using TargetSpace = ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >;
  private:
    typedef typename TargetSpace::ctype RT;

    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
    static const int embeddedDim = EmbeddedTangentVector::dimension;

  public:

    /** \brief The type used for derivatives */
    typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;

    /** \brief Constructor
     * \param localFiniteElement A Lagrangian finite element that provides the interpolation points
     * \param coefficients Values of the function at the Lagrange points
     */
    LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
                             const std::vector<TargetSpace>& coefficients)
      : localFiniteElement_(localFiniteElement),
      coefficients_(coefficients),
      translationCoefficients_(coefficients.size())
    {
      assert(localFiniteElement.localBasis().size() == coefficients.size());
      using namespace Dune::Indices;

      for (size_t i=0; i<coefficients.size(); i++)
        translationCoefficients_[i] = coefficients[i][_0].globalCoordinates();

      std::vector<Rotation<field_type,3> > orientationCoefficients(coefficients.size());
      for (size_t i=0; i<coefficients.size(); i++)
        orientationCoefficients[i] = coefficients[i][_1];

      orientationFunction_ = std::make_unique<LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > (localFiniteElement,orientationCoefficients);
    }

    /** \brief Rebind the FEFunction to another TargetSpace */
    template<class U>
    struct rebind
    {
      using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U>;
    };

    /** \brief The number of Lagrange points */
    unsigned int size() const
    {
      return localFiniteElement_.size();
    }

    /** \brief The type of the reference element */
    Dune::GeometryType type() const
    {
      return localFiniteElement_.type();
    }

    /** \brief Evaluate the function */
    TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const
    {
      using namespace Dune::Indices;

      TargetSpace result;

      // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
      std::vector<Dune::FieldVector<ctype,1> > w;
      localFiniteElement_.localBasis().evaluateFunction(local,w);

      result[_0] = FieldVector<field_type,3>(0.0);
      for (size_t i=0; i<w.size(); i++)
        result[_0].globalCoordinates().axpy(w[i][0], translationCoefficients_[i]);

      result[_1] = orientationFunction_->evaluate(local);

      return result;
    }

    /** \brief Evaluate the derivative of the function */
    DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
    {
      // the function value at the point where we are evaluating the derivative
      TargetSpace q = evaluate(local);

      // Actually compute the derivative
      return evaluateDerivative(local,q);
    }

    /** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
     *        \param local Local coordinates in the reference element where to evaluate the derivative
     *        \param q Value of the local function at 'local'.  If you provide something wrong here the result will be wrong, too!
     */
    DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
                                      const TargetSpace& q) const
    {
      using namespace Dune::Indices;

      DerivativeType result(0);

      // get translation part
      std::vector<Dune::FieldMatrix<ctype,1,dim> > sfDer(translationCoefficients_.size());
      localFiniteElement_.localBasis().evaluateJacobian(local, sfDer);

      for (size_t i=0; i<translationCoefficients_.size(); i++)
        for (int j=0; j<3; j++)
          result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);

      // get orientation part
      Dune::FieldMatrix<field_type,4,dim> qResult = orientationFunction_->evaluateDerivative(local,q[_1]);

      for (int i=0; i<4; i++)
        for (int j=0; j<dim; j++)
          result[3+i][j] = qResult[i][j];

      return result;
    }

    /** \brief Get the i'th base coefficient. */
    const TargetSpace& coefficient(int i) const
    {
      return coefficients_[i];
    }
  private:

    /** \brief The scalar local finite element, which provides the weighting factors
     *        \todo We really only need the local basis
     */
    const LocalFiniteElement& localFiniteElement_;

    // The coefficients of this interpolation rule
    std::vector<TargetSpace> coefficients_;

    // The coefficients again.  Yes, we store it twice:  To evaluate the interpolation rule efficiently
    // we need access to the coefficients for the various factor spaces separately.
    std::vector<Dune::FieldVector<field_type,3> > translationCoefficients_;

    std::unique_ptr<LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type, 3> > > orientationFunction_;


  };

}  // namespace Dune::GFE

#endif