diff --git a/dune/tectonic/ellipticenergy.hh b/dune/tectonic/ellipticenergy.hh
index 11924f548c2b88bc5e5845ce1f0eac9f83aa3741..876df28cc6e8d4cfe70a35ef8e6685503284fb28 100644
--- a/dune/tectonic/ellipticenergy.hh
+++ b/dune/tectonic/ellipticenergy.hh
@@ -5,7 +5,6 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/stdstreams.hh>
 
-#include <dune/fufem/interval.hh>
 #include <dune/fufem/arithmetic.hh>
 
 #include "localfriction.hh"
diff --git a/dune/tectonic/globalnonlinearity.hh b/dune/tectonic/globalnonlinearity.hh
index 5531162231a7f7038ccd6f92329bce9b61c639de..0a899a56db4cd7d81d3472920c06f9a4ff1b25e9 100644
--- a/dune/tectonic/globalnonlinearity.hh
+++ b/dune/tectonic/globalnonlinearity.hh
@@ -7,6 +7,8 @@
 #include <dune/istl/bvector.hh>
 #include <dune/istl/matrixindexset.hh>
 
+#include <dune/solvers/common/interval.hh>
+
 #include "localfriction.hh"
 
 template <class Matrix, class Vector> class GlobalNonlinearity {
@@ -45,16 +47,17 @@ template <class Matrix, class Vector> class GlobalNonlinearity {
   }
 
   void directionalDomain(Vector const &, Vector const &,
-                         Interval<double> &dom) const {
+                         Dune::Solvers::Interval<double> &dom) const {
     dom[0] = -std::numeric_limits<double>::max();
     dom[1] = std::numeric_limits<double>::max();
   }
 
   void directionalSubDiff(Vector const &u, Vector const &v,
-                          Interval<double> &subdifferential) const {
+                          Dune::Solvers::Interval<double> &subdifferential)
+      const {
     subdifferential[0] = subdifferential[1] = 0;
     for (size_t i = 0; i < u.size(); ++i) {
-      Interval<double> D;
+      Dune::Solvers::Interval<double> D;
       auto const res = restriction(i);
       res->directionalSubDiff(u[i], v[i], D);
       subdifferential[0] += D[0];
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index baac090430c2d67e28830c39267c64c072bfe120..04b8121a9b3d06c82a15c44bed2fa4f2471ef67c 100644
--- a/dune/tectonic/localfriction.hh
+++ b/dune/tectonic/localfriction.hh
@@ -8,7 +8,7 @@
 #include <dune/common/fmatrix.hh>
 
 #include <dune/fufem/arithmetic.hh>
-#include <dune/fufem/interval.hh>
+#include <dune/solvers/common/interval.hh>
 
 #include "frictionpotential.hh"
 
@@ -54,7 +54,7 @@ template <size_t dimension> class LocalFriction {
   // directional subdifferential: at u on the line u + t*v
   // u and v are assumed to be non-zero
   void directionalSubDiff(VectorType const &x, VectorType const &v,
-                          Interval<double> &D) const {
+                          Dune::Solvers::Interval<double> &D) const {
     double const xnorm = x.two_norm();
     if (xnorm <= 0.0)
       D[0] = D[1] = func->differential(0.0) * v.two_norm();
@@ -113,7 +113,7 @@ template <size_t dimension> class LocalFriction {
   }
 
   void directionalDomain(VectorType const &, VectorType const &,
-                         Interval<double> &dom) const {
+                         Dune::Solvers::Interval<double> &dom) const {
     dom[0] = -std::numeric_limits<double>::max();
     dom[1] = std::numeric_limits<double>::max();
   }
diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index ee5388cea7cfb101dcffa8e5e77e4fdcd590cd80..d19db04f38774cfa81dbd9fef7fde77ce58d5ee9 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -11,6 +11,7 @@
 #include "dune/solvers/computeenergy.hh"
 
 #include <dune/fufem/arithmetic.hh>
+#include <dune/solvers/common/interval.hh>
 #include <dune/tnnmg/problem-classes/bisection.hh>
 
 #include "globalnonlinearity.hh"
@@ -98,7 +99,7 @@ template <class ConvexProblem> class MyBlockProblem {
     psi(computeDirectionalA(problem.A, v),
         computeDirectionalb(problem.A, problem.f, u, v), problem.phi, u, v);
 
-    Interval<double> D;
+    Dune::Solvers::Interval<double> D;
     psi.subDiff(0, D);
     // NOTE: Numerical instability can actually get us here
     if (D[1] > 0)
diff --git a/dune/tectonic/mydirectionalconvexfunction.hh b/dune/tectonic/mydirectionalconvexfunction.hh
index 9c792a2d92358412256d1aa483b4067890fe184e..83cd3b40522969ee8f68d1a852cc7e879a283bed 100644
--- a/dune/tectonic/mydirectionalconvexfunction.hh
+++ b/dune/tectonic/mydirectionalconvexfunction.hh
@@ -5,7 +5,7 @@
 #define MY_DIRECTIONAL_CONVEX_FUNCTION_HH
 
 #include <dune/fufem/arithmetic.hh>
-#include <dune/fufem/interval.hh>
+#include <dune/solvers/common/interval.hh>
 
 /*
   1/2 <A(u + hv),u + hv> - <b, u + hv>
@@ -41,7 +41,7 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
 
   double linearPart() const { return b; }
 
-  void subDiff(double x, Interval<double> &D) const {
+  void subDiff(double x, Dune::Solvers::Interval<double> &D) const {
     Vector uxv = u;
     Arithmetic::addProduct(uxv, x, v);
     phi.directionalSubDiff(uxv, v, D);
@@ -49,7 +49,7 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
     D[1] += A * x - b;
   }
 
-  void domain(Interval<double> &domain) const {
+  void domain(Dune::Solvers::Interval<double> &domain) const {
     domain[0] = this->dom[0];
     domain[1] = this->dom[1];
   }
@@ -62,7 +62,7 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
   Vector const &u;
   Vector const &v;
 
-  Interval<double> dom;
+  Dune::Solvers::Interval<double> dom;
 };
 
 #endif
diff --git a/src/state/compute_state_dieterich_euler.cc b/src/state/compute_state_dieterich_euler.cc
index 87f657244c9e7b4d7ea1982ac753a857807e9519..103c56f2c3d514c086ad91cfdf58bf71a4191220 100644
--- a/src/state/compute_state_dieterich_euler.cc
+++ b/src/state/compute_state_dieterich_euler.cc
@@ -8,7 +8,7 @@
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
-#include <dune/fufem/interval.hh>
+#include <dune/solvers/common/interval.hh>
 
 /* Dieterich's law:
 
@@ -49,12 +49,12 @@ class DieterichNonlinearity {
   DieterichNonlinearity() {}
 
   void directionalSubDiff(VectorType const &u, VectorType const &v,
-                          Interval<double> &D) const {
+                          Dune::Solvers::Interval<double> &D) const {
     D[0] = D[1] = v[0] * (-std::exp(-u[0]));
   }
 
   void directionalDomain(VectorType const &, VectorType const &,
-                         Interval<double> &dom) const {
+                         Dune::Solvers::Interval<double> &dom) const {
     dom[0] = -std::numeric_limits<double>::max();
     dom[1] = std::numeric_limits<double>::max();
   }