diff --git a/dune/stuff/functions/expression.hh b/dune/stuff/functions/expression.hh
index 52092935e85bea9fc7982975b95fc4fdc0bea6c4..23dda75cfa127f3d6e6d5aa58998505f144b502e 100644
--- a/dune/stuff/functions/expression.hh
+++ b/dune/stuff/functions/expression.hh
@@ -210,12 +210,21 @@ public:
     const Common::Configuration cfg         = config.has_sub(sub_name) ? config.sub(sub_name) : config;
     const Common::Configuration default_cfg = default_config();
     // create
-    return Common::make_unique<ThisType>(cfg.get("variable", default_cfg.get<std::string>("variable")),
+    return cfg.has_key("gradient") ? Common::make_unique<ThisType>(
+                                         cfg.get("variable", default_cfg.get<std::string>("variable")),
                                          cfg.get("expression", default_cfg.get<std::vector<std::string>>("expression")),
                                          cfg.get("order", default_cfg.get<size_t>("order")),
-                                         cfg.get("name", default_cfg.get<std::string>("name")));
+                                         cfg.get("name", default_cfg.get<std::string>("name")),
+                                         cfg.get<Dune::FieldMatrix<std::string, dimRange, dimDomain>>("gradient"))
+                                   : Common::make_unique<ThisType>(
+                                         cfg.get("variable", default_cfg.get<std::string>("variable")),
+                                         cfg.get("expression", default_cfg.get<std::vector<std::string>>("expression")),
+                                         cfg.get("order", default_cfg.get<size_t>("order")),
+                                         cfg.get("name", default_cfg.get<std::string>("name")),
+                                         Dune::FieldMatrix<std::string, 0, 0>());
   } // ... create(...)
 
+  // constructors taking a std::vector< std::vector< std::string > > for the jacobian
   Expression(const std::string variable, const std::string expression,
              const size_t ord = default_config().get<size_t>("order"), const std::string nm = static_id(),
              const std::vector<std::vector<std::string>> gradient_expressions = std::vector<std::vector<std::string>>())
@@ -236,6 +245,29 @@ public:
     build_gradients(variable, gradient_expressions);
   }
 
+  // constructors taking a FieldMatrix for the jacobian
+  template <int rows, int cols>
+  Expression(const std::string variable, const std::string expression,
+             const size_t ord = default_config().get<size_t>("order"), const std::string nm = static_id(),
+             const Dune::FieldMatrix<std::string, rows, cols> jacobian = Dune::FieldMatrix<std::string, 0, 0>())
+    : function_(new MathExpressionFunctionType(variable, expression))
+    , order_(ord)
+    , name_(nm)
+  {
+    build_gradients(variable, jacobian);
+  }
+
+  template <int rows, int cols>
+  Expression(const std::string variable, const std::vector<std::string> expressions,
+             const size_t ord = default_config().get<size_t>("order"), const std::string nm = static_id(),
+             const Dune::FieldMatrix<std::string, rows, cols> jacobian = Dune::FieldMatrix<std::string, 0, 0>())
+    : function_(new MathExpressionFunctionType(variable, expressions))
+    , order_(ord)
+    , name_(nm)
+  {
+    build_gradients(variable, jacobian);
+  }
+
   virtual std::string name() const override
   {
     return name_;
@@ -306,6 +338,20 @@ private:
       }
   } // ... build_gradients(...)
 
+  template <int rows, int cols>
+  void build_gradients(const std::string variable, Dune::FieldMatrix<std::string, rows, cols> jacobian)
+  {
+    assert(jacobian.rows == 0 || jacobian.rows >= dimRange);
+    assert(jacobian.cols == 0 || jacobian.cols >= dimDomain);
+    if (jacobian.rows > 0 && jacobian.cols > 0)
+      for (size_t rr = 0; rr < dimRange; ++rr) {
+        std::vector<std::string> gradient_expression;
+        for (size_t cc = 0; cc < dimDomain; ++cc)
+          gradient_expression.emplace_back(jacobian[rr][cc]);
+        gradients_.emplace_back(new MathExpressionGradientType(variable, gradient_expression));
+      }
+  } // ... build_gradients(...)
+
   std::shared_ptr<const MathExpressionFunctionType> function_;
   size_t order_;
   std::string name_;