diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index adfb45e8237acfe78afa5f47be624517213cd94a..8b277d92c6304409248ba0704756d6d2a16a7141 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -16,6 +16,7 @@
 #include <dune/istl/umfpack.hh>
 
 #include <dune/mmdg/nonconformingp1vtkfunction.hh>
+#include <dune/mmdg/eigenvaluehelper.hh>
 
 template<class GridView, class Mapper, class Problem>
 class DG
@@ -532,7 +533,7 @@ protected:
   const Scalar getLargestEigenvalue (const auto& matrix) const
   {
     Coordinate eigenvalues;
-    Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
+    EigenvalueHelper::eigenValues(matrix, eigenvalues);
     return eigenvalues.infinity_norm();
   }
 
diff --git a/dune/mmdg/eigenvaluehelper.hh b/dune/mmdg/eigenvaluehelper.hh
new file mode 100644
index 0000000000000000000000000000000000000000..63bde0fdbac7fd10c587c741c768be6247d42d66
--- /dev/null
+++ b/dune/mmdg/eigenvaluehelper.hh
@@ -0,0 +1,123 @@
+#ifndef DUNE_MMDG_EIGENVALUEHELPER_HH
+#define DUNE_MMDG_EIGENVALUEHELPER_HH
+
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/math.hh>
+
+/** \file
+ * \brief Eigenvalue computations for the Dune::FieldMatrix class
+ * NOTE: taken from DUNE-common, fmatrixev.hh, with adapted tolerance
+ */
+class EigenvalueHelper {
+
+public:
+  /** \brief calculates the eigenvalues of a symmetric field matrix
+      \param[in]  matrix matrix eigenvalues are calculated for
+      \param[out] eigenvalues Dune::FieldVector that contains eigenvalues in
+                  ascending order
+   */
+  template <typename K>
+  static void eigenValues(const Dune::FieldMatrix<K, 1, 1>& matrix,
+                          Dune::FieldVector<K, 1>& eigenvalues)
+  {
+    eigenvalues[0] = matrix[0][0];
+  }
+
+  /** \brief calculates the eigenvalues of a symmetric field matrix
+      \param[in]  matrix matrix eigenvalues are calculated for
+      \param[out] eigenvalues Dune::FieldVector that contains eigenvalues in
+                  ascending order
+   */
+  template <typename K>
+  static void eigenValues(const Dune::FieldMatrix<K, 2, 2>& matrix,
+                          Dune::FieldVector<K, 2>& eigenvalues)
+  {
+    using std::sqrt;
+    const K detM = matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1];
+    const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
+    K q = p * p - detM;
+    if( q < 0 && q > -1e-13/*-1e-14*/ ) q = 0;
+    if (q < 0)
+    {
+      std::cout << std::numeric_limits<K>::epsilon() << "\n" << q << "\n" << matrix << std::endl;
+      // Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
+      DUNE_THROW(Dune::MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
+    }
+
+    // get square root
+    q = sqrt(q);
+
+    // store eigenvalues in ascending order
+    eigenvalues[0] = p - q;
+    eigenvalues[1] = p + q;
+  }
+
+  /** \brief Calculates the eigenvalues of a symmetric 3x3 field matrix
+      \param[in]  matrix matrix eigenvalues are calculated for
+      \param[out] eigenvalues Eigenvalues in ascending order
+
+      \note If the input matrix is not symmetric the behavior of this method is undefined.
+
+      This implementation was adapted from the pseudo-code (Python?) implementation found on
+      http://en.wikipedia.org/wiki/Eigenvalue_algorithm  (retrieved late August 2014).
+      Wikipedia claims to have taken it from
+        Smith, Oliver K. (April 1961), Eigenvalues of a symmetric 3 × 3 matrix.,
+        Communications of the ACM 4 (4): 168, doi:10.1145/355578.366316
+   */
+  template <typename K>
+  static void eigenValues(const Dune::FieldMatrix<K, 3, 3>& matrix,
+                          Dune::FieldVector<K, 3>& eigenvalues)
+  {
+    using std::sqrt;
+    using std::acos;
+    using real_type = typename Dune::FieldTraits<K>::real_type;
+    const K pi = Dune::MathematicalConstants<K>::pi();
+    K p1 = matrix[0][1]*matrix[0][1] + matrix[0][2]*matrix[0][2] + matrix[1][2]*matrix[1][2];
+
+    if (p1 <= 1e-8)
+    {
+      // A is diagonal.
+      eigenvalues[0] = matrix[0][0];
+      eigenvalues[1] = matrix[1][1];
+      eigenvalues[2] = matrix[2][2];
+    }
+    else
+    {
+      // q = trace(A)/3
+      K q = 0;
+      for (int i=0; i<3; i++)
+        q += matrix[i][i]/3.0;
+
+      K p2 = (matrix[0][0] - q)*(matrix[0][0] - q) + (matrix[1][1] - q)*(matrix[1][1] - q) + (matrix[2][2] - q)*(matrix[2][2] - q) + 2 * p1;
+      K p = sqrt(p2 / 6);
+      // B = (1 / p) * (A - q * I);       // I is the identity matrix
+      Dune::FieldMatrix<K,3,3> B;
+      for (int i=0; i<3; i++)
+        for (int j=0; j<3; j++)
+          B[i][j] = (real_type(1.0)/p) * (matrix[i][j] - q*(i==j));
+
+      K r = B.determinant() / 2.0;
+
+      // In exact arithmetic for a symmetric matrix  -1 <= r <= 1
+      // but computation error can leave it slightly outside this range.
+      K phi;
+      if (r <= -1)
+        phi = pi / 3.0;
+      else if (r >= 1)
+        phi = 0;
+      else
+        phi = acos(r) / 3;
+
+      // the eigenvalues satisfy eig[2] <= eig[1] <= eig[0]
+      eigenvalues[2] = q + 2 * p * cos(phi);
+      eigenvalues[0] = q + 2 * p * cos(phi + (2*pi/3));
+      eigenvalues[1] = 3 * q - eigenvalues[0] - eigenvalues[2];     // since trace(matrix) = eig1 + eig2 + eig3
+    }
+  }
+};
+
+#endif