diff --git a/dune/gfe/localprojectedfefunction.hh b/dune/gfe/localprojectedfefunction.hh
index b408ab44df6c31efbef13cf1f60341aff6eee2a2..b04c92974888e1b24645d276e00ae8f1ab1a09dd 100644
--- a/dune/gfe/localprojectedfefunction.hh
+++ b/dune/gfe/localprojectedfefunction.hh
@@ -241,8 +241,8 @@ namespace Dune {
      * \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,Rotation<field_type,3> >
+    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;
@@ -459,6 +459,42 @@ namespace Dune {
         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
       {