diff --git a/src/trivialnonlinearity.hh b/src/trivialnonlinearity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8b30a169d283771d26d059b8043c4fdf9d132014
--- /dev/null
+++ b/src/trivialnonlinearity.hh
@@ -0,0 +1,137 @@
+/* -*- mode:c++; mode:semantic -*- */
+
+#ifndef TRIVIALNONLINEARITY_HH
+#define TRIVIALNONLINEARITY_HH
+
+#include <dune/tnnmg/problem-classes/nonlinearity.hh>
+
+template <class LocalVectorTypeTEMPLATE = Dune::FieldVector<double, 1>,
+          class LocalMatrixTypeTEMPLATE = Dune::FieldMatrix<double, 1, 1>>
+class TrivialNonlinearity
+    : Nonlinearity<LocalVectorTypeTEMPLATE, LocalMatrixTypeTEMPLATE> {
+public:
+  typedef typename Nonlinearity<LocalVectorTypeTEMPLATE,
+                                LocalMatrixTypeTEMPLATE>::IndexSet IndexSet;
+  typedef typename Nonlinearity<LocalVectorTypeTEMPLATE,
+                                LocalMatrixTypeTEMPLATE>::LocalVectorType
+  LocalVectorType;
+  typedef typename Nonlinearity<LocalVectorTypeTEMPLATE,
+                                LocalMatrixTypeTEMPLATE>::MatrixType MatrixType;
+  typedef typename Nonlinearity<LocalVectorTypeTEMPLATE,
+                                LocalMatrixTypeTEMPLATE>::VectorType VectorType;
+
+  enum {
+    block_size = Nonlinearity<LocalVectorTypeTEMPLATE,
+                              LocalMatrixTypeTEMPLATE>::block_size
+  };
+
+  //! Evaluate the nonlinearity at v.
+  virtual double operator()(const VectorType& v) const { return 0; }
+
+  //! Add the gradient of the nonlinearity at v to the vector gradient.
+  virtual void addGradient(const VectorType& v, VectorType& gradient) const {}
+
+  //! Add the Hessian matrix of the nonlinearity at v to the matrix Hessian.
+  virtual void addHessian(const VectorType& v, MatrixType& hessian) const {}
+
+  //! Add the indices of the Hessian matrix to the index set.
+  virtual void addHessianIndices(IndexSet& indices) const {}
+
+  /** \brief Set the internal position vector u_pos to v.
+   *
+   * This is only needed if the nonlinearity does not decouple in the Euclidean
+   *directions.
+   * If the nonlinearity decouples in the Euclidean directions this can be
+   *empty.
+   */
+  virtual void setVector(const VectorType& v) {}
+
+  /** \brief Update the (i,j)-th entry of the internal position vector u_pos to
+   *x.
+   *
+   * \param i global index
+   * \param j local index
+   * \param x new value of the entry (u_pos)_(i,j)
+   */
+  virtual void updateEntry(int i, double x, int j) {}
+
+  /** \brief Compute the subdifferential of the nonlinearity restricted to the
+   * line u_pos' +t e_(i,j) at t=0.
+   *
+   * Here e_(i,j) is the (i,j)-th Euclidean unit vector,
+   * and u_pos' is the internal position vector u_pos with the (i,j)-the entry
+   *replaced by x.
+   * If the nonlinearity decouples in the Euclidean directions this is simply
+   *the (i,j)-th
+   * component of the subdifferential.
+   *
+   * \param i global index
+   * \param j local index
+   * \param x value of the (i,j)-th entry of position to evaluate the
+   *nonlinearity at
+   * \param[out] D the subdifferential
+   */
+  virtual void subDiff(int i, double x, Interval<double>& D, int j) const {
+    D[0] = D[1] = 0.0;
+  }
+
+  /** \brief Return the regularity of the nonlinearity restricted to the
+   * line u_pos' +t e_(i,j) at t=0.
+   *
+   * Here e_(i,j) is the (i,j)-th Euclidean unit vector,
+   * and u_pos' is the internal position vector u_pos with the (i,j)-the entry
+   *replaced by x.
+   * Usually this will be the third derivative or a local Lipschitz constant of
+   *the second
+   * derivative. Note that if the subdifferential is set-valued at this
+   *position, this
+   * value will normally be infinity.
+   *
+   * \param i global index
+   * \param j local index
+   * \param x value of the (i,j)-th entry of position to evaluate the
+   *nonlinearity at
+   * \returns a value measuring the regularity
+   */
+  virtual double regularity(int i, double x, int j) const { return 0; }
+
+  /** \brief Compute the domain of the nonlinearity restricted to the
+   * line t e_(i,j) at t=0 where e_(i,j) is the (i,j)-th Euclidean unit vector.
+   *
+   * Notice that this does not depend on the position since the nonsmooth
+   * part of the nonlinearity must decouple in the Euclidean directions.
+   *
+   * \param i global index
+   * \param j local index
+   * \param[out] dom the domain
+   */
+  virtual void domain(int i, Interval<double>& dom, int j) const {
+    dom[1] = std::numeric_limits<double>::max();
+    dom[0] = std::numeric_limits<double>::min();
+  }
+
+  /** \brief Compute a local domain around t=0 where the nonlinearity restricted
+   *to the
+   * line u_pos' +t e_(i,j) is smooth.
+   *
+   * Here e_(i,j) is the (i,j)-th Euclidean unit vector,
+   * and u_pos' is the internal position vector u_pos with the (i,j)-the entry
+   *replaced by x.
+   * This defaults to the whole domain. Usually there is no need to use
+   *different values here.
+   *
+   * \param i global index
+   * \param j local index
+   * \param x value of the (i,j)-th entry of position to evaluate the
+   *nonlinearity at
+   * \param[out] dom the domain
+   */
+  virtual void smoothDomain(int i, double x, double regularity,
+                            Interval<double>& dom, int j) const {
+    dom[1] = std::numeric_limits<double>::max();
+    dom[0] = std::numeric_limits<double>::min();
+    return;
+  }
+};
+
+#endif