diff --git a/dune/stuff/functions/expression.hh b/dune/stuff/functions/expression.hh
index 7f8b5f38b039ef14333a86774be393981778d6d1..11b748dd19b6a15e69965c2f692f611e434be193 100644
--- a/dune/stuff/functions/expression.hh
+++ b/dune/stuff/functions/expression.hh
@@ -16,6 +16,7 @@
 
 #include <dune/stuff/common/configuration.hh>
 #include <dune/stuff/common/exceptions.hh>
+#include <dune/stuff/common/parallel/threadstorage.hh>
 
 #include "expression/base.hh"
 #include "interfaces.hh"
@@ -216,15 +217,15 @@ public:
     bool failure = false;
     std::string error_type;
     for (size_t rr = 0; rr < dimRange; ++rr) {
-      tmp_row_ = ret[rr];
+      *tmp_row_ = ret[rr];
       for (size_t cc = 0; cc < dimRangeCols; ++cc) {
-        if (DSC::isnan(tmp_row_[cc])) {
+        if (DSC::isnan(tmp_row_->operator[](cc))) {
           failure    = true;
           error_type = "NaN";
-        } else if (DSC::isinf(tmp_row_[cc])) {
+        } else if (DSC::isinf(tmp_row_->operator[](cc))) {
           failure    = true;
           error_type = "inf";
-        } else if (std::abs(tmp_row_[cc]) > (0.9 * std::numeric_limits<double>::max())) {
+        } else if (std::abs(tmp_row_->operator[](cc)) > (0.9 * std::numeric_limits<double>::max())) {
           failure    = true;
           error_type = "an unlikely value";
         }
@@ -243,7 +244,7 @@ public:
                          << xx
                          << "\n"
                          << "The result was:                       "
-                         << tmp_row_[cc]
+                         << tmp_row_->operator[](cc)
                          << "\n\n"
                          << "You can disable this check by defining DUNE_STUFF_FUNCTIONS_EXPRESSION_DISABLE_CHECKS\n");
       }
@@ -297,11 +298,11 @@ private:
   template <size_t rC>
   void evaluate_helper(const DomainType& xx, RangeType& ret, internal::ChooseVariant<rC>) const
   {
-    function_->evaluate(xx, tmp_vector_);
+    function_->evaluate(xx, *tmp_vector_);
     for (size_t rr = 0; rr < dimRange; ++rr) {
       auto& retRow = ret[rr];
       for (size_t cc = 0; cc < dimRangeCols; ++cc)
-        retRow[cc] = tmp_vector_[rr * dimRangeCols + cc];
+        retRow[cc] = (*tmp_vector_)[rr * dimRangeCols + cc];
     }
   } // ... evaluate_helper(...)
 
@@ -380,11 +381,10 @@ private:
   } // ... get_gradient(...)
 
   std::shared_ptr<const MathExpressionFunctionType> function_;
-  size_t order_;
-  std::string name_;
-  mutable FieldVector<RangeFieldType, dimRange * dimRangeCols> tmp_vector_;
-  mutable FieldVector<RangeFieldType, dimRangeCols> tmp_row_;
-  mutable FieldVector<RangeFieldType, dimDomain> tmp_gradient_row_;
+  const size_t order_;
+  const std::string name_;
+  mutable typename DS::PerThreadValue<FieldVector<RangeFieldType, dimRange * dimRangeCols>> tmp_vector_;
+  mutable typename DS::PerThreadValue<FieldVector<RangeFieldType, dimRangeCols>> tmp_row_;
   std::vector<std::vector<std::shared_ptr<const MathExpressionGradientType>>> gradients_;
 }; // class Expression
 
diff --git a/dune/stuff/functions/expression/base.hh b/dune/stuff/functions/expression/base.hh
index 1c570e111df4b73a3a867509960ab0f0681f076a..4dadb31a22553e7425e5bb386a7b861292fcc38c 100644
--- a/dune/stuff/functions/expression/base.hh
+++ b/dune/stuff/functions/expression/base.hh
@@ -87,12 +87,14 @@ public:
   void evaluate(const Dune::FieldVector<DomainFieldType, dimDomain>& arg,
                 Dune::FieldVector<RangeFieldType, dimRange>& ret) const
   {
+    mutex_.lock();
     // copy arg
     for (typename Dune::FieldVector<DomainFieldType, dimDomain>::size_type ii = 0; ii < dimDomain; ++ii)
       *(arg_[ii]) = arg[ii];
     // copy ret
     for (typename Dune::FieldVector<RangeFieldType, dimRange>::size_type ii = 0; ii < dimRange; ++ii)
       ret[ii] = op_[ii]->Val();
+    mutex_.unlock();
   }
 
   /**
@@ -100,6 +102,7 @@ public:
    */
   void evaluate(const Dune::DynamicVector<DomainFieldType>& arg, Dune::DynamicVector<RangeFieldType>& ret) const
   {
+    mutex_.lock();
     // check for sizes
     assert(arg.size() > 0);
     if (ret.size() != dimRange)
@@ -110,11 +113,13 @@ public:
     // copy ret
     for (typename Dune::DynamicVector<RangeFieldType>::size_type ii = 0; ii < dimRange; ++ii)
       ret[ii] = op_[ii]->Val();
+    mutex_.unlock();
   }
 
   void evaluate(const Dune::FieldVector<DomainFieldType, dimDomain>& arg,
                 Dune::DynamicVector<RangeFieldType>& ret) const
   {
+    mutex_.lock();
     // check for sizes
     if (ret.size() != dimRange)
       ret = Dune::DynamicVector<RangeFieldType>(dimRange);
@@ -124,6 +129,7 @@ public:
     // copy ret
     for (typename Dune::DynamicVector<RangeFieldType>::size_type ii = 0; ii < dimRange; ++ii)
       ret[ii] = op_[ii]->Val();
+    mutex_.unlock();
   }
 
   /**
@@ -131,6 +137,7 @@ public:
    */
   void evaluate(const Dune::DynamicVector<DomainFieldType>& arg, Dune::FieldVector<RangeFieldType, dimRange>& ret) const
   {
+    mutex_.lock();
     assert(arg.size() > 0);
     // copy arg
     for (size_t ii = 0; ii < std::min(dimDomain, arg.size()); ++ii)
@@ -138,6 +145,7 @@ public:
     // copy ret
     for (size_t ii = 0; ii < dimRange; ++ii)
       ret[ii] = op_[ii]->Val();
+    mutex_.unlock();
   }
 
   void report(const std::string _name = "dune.stuff.function.mathexpressionbase", std::ostream& stream = std::cout,
@@ -209,6 +217,7 @@ private:
   RVar* var_arg_[dimDomain];
   RVar* vararray_[dimDomain];
   ROperation* op_[dimRange];
+  mutable std::mutex mutex_;
 }; // class MathExpressionBase