diff --git a/src/amdis/gridfunctions/DiscreteFunction.hpp b/src/amdis/gridfunctions/DiscreteFunction.hpp
index 609a1af50a983425db12edd6f5f726e180cfb515..11165a97828f202e8bb30499f7e9755f1ec29072 100644
--- a/src/amdis/gridfunctions/DiscreteFunction.hpp
+++ b/src/amdis/gridfunctions/DiscreteFunction.hpp
@@ -10,7 +10,9 @@
 #include <dune/typetree/childextraction.hh>
 
 #include <amdis/GridFunctions.hpp>
+#include <amdis/LinearAlgebra.hpp>
 #include <amdis/typetree/FiniteElementType.hpp>
+#include <amdis/typetree/RangeType.hpp>
 #include <amdis/typetree/TreePath.hpp>
 
 namespace AMDiS
diff --git a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
index 52e91d28bb73b4ba84c62a65408329040e65b564..bc81d3c50fbcb63415731b9aaa533e8eed0f5f60 100644
--- a/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
+++ b/src/amdis/gridfunctions/DiscreteFunction.inc.hpp
@@ -29,6 +29,12 @@ public:
     , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
   {}
 
+  LocalFunction(LocalFunction const& other)
+    : globalFunction_(other.globalFunction_)
+    , localView_(globalFunction_.basis().localView())
+    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
+  {}
+
   /// \brief Bind the LocalView to the element
   void bind(Element const& element)
   {
@@ -100,6 +106,12 @@ public:
     , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
   {}
 
+  GradientLocalFunction(GradientLocalFunction const& other)
+    : globalFunction_(other.globalFunction_)
+    , localView_(globalFunction_.basis().localView())
+    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
+  {}
+
   void bind(Element const& element)
   {
     localView_.bind(element);
@@ -144,7 +156,7 @@ typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
 LocalFunction::operator()(Domain const& x) const
 {
   assert( bound_ );
-  auto y = Range(0);
+  Range y(0);
 
   auto&& coefficients = *globalFunction_.dofVector_;
   auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
@@ -188,9 +200,7 @@ typename DiscreteFunction<GB,VT,TP>::GradientLocalFunction::Range DiscreteFuncti
 GradientLocalFunction::operator()(Domain const& x) const
 {
   assert( bound_ );
-  Range dy;
-  for (std::size_t j = 0; j < dy.size(); ++j)
-    dy[j] = 0;
+  Range dy(0);
 
   auto&& coefficients = *globalFunction_.dofVector_;
   auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
@@ -203,6 +213,7 @@ GradientLocalFunction::operator()(Domain const& x) const
 
     auto&& fe = node.finiteElement();
     auto&& localBasis = fe.localBasis();
+    std::size_t size = localBasis.size();
 
     NodeCache<TYPEOF(node)> localBasisCache(localBasis);
     auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element().type(),x);
@@ -218,7 +229,7 @@ GradientLocalFunction::operator()(Domain const& x) const
     // Get range entry associated to this node
     auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
 
-    for (std::size_t i = 0; i < localBasis.size(); ++i) {
+    for (std::size_t i = 0; i < size; ++i) {
       auto&& multiIndex = localView_.index(node.localIndex(i));
 
       // Get coefficient associated to i-th shape function