diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 8b277d92c6304409248ba0704756d6d2a16a7141..c842b351e252e2453be925ea32ca9dd39d7c12aa 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -13,7 +13,7 @@
 
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/solver.hh>
-#include <dune/istl/umfpack.hh>
+#include <dune/istl/cholmod.hh>
 
 #include <dune/mmdg/nonconformingp1vtkfunction.hh>
 #include <dune/mmdg/eigenvaluehelper.hh>
@@ -57,7 +57,8 @@ public:
     A->compress();
 
     Dune::InverseOperatorResult result;
-    Dune::UMFPack<Matrix> solver(*A);
+    Dune::Cholmod<Vector> solver;
+    solver.setMatrix(*A);
     solver.apply(d, b, result);
 
     //write solution to a vtk file
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 3511adf96acb239aa34eec93afa1b68ee0bbea0b..a83541dbf9d2d190a9b2c948f63fc04c0adfa8df 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -41,7 +41,8 @@ public:
     Base::A->compress();
 
     Dune::InverseOperatorResult result;
-    Dune::UMFPack<Matrix> solver(*Base::A); //TODO: ilu0bicgSTAB, cg (ilo0, ssor), gmres
+    Dune::Cholmod<Vector> solver;
+    solver.setMatrix(*Base::A);
     solver.apply(Base::d, Base::b, result);
 
     //write solution to a vtk file
diff --git a/dune/mmdg/problems/mmdgproblem3.hh b/dune/mmdg/problems/mmdgproblem3.hh
index d3f284c1c2a69d257bdf3acd9a11583477eb876d..cf7f896a9aa6b61d19dc18a5cecf7d20c9b1b0fa 100644
--- a/dune/mmdg/problems/mmdgproblem3.hh
+++ b/dune/mmdg/problems/mmdgproblem3.hh
@@ -20,7 +20,7 @@ public:
   //the exact bulk solution at position pos
   Scalar exactSolution (const Coordinate& pos) const
   {
-    Scalar xPlusy = pos[0] + pos[1];// pos * Coordinate(1.0);
+    Scalar xPlusy = pos[0] + pos[1];
     Scalar solution = std::exp(xPlusy);
     return (xPlusy < 1.0) ? solution :
       0.5 * solution + std::exp(1.0) * (0.5 + 1.5 / sqrt(2.0) * d_);
@@ -41,7 +41,7 @@ public:
   //bulk source term at position pos
   Scalar q (const Coordinate& pos) const
   {
-    Scalar xPlusy = pos[0] + pos[1]; //pos * Coordinate(1.0);
+    Scalar xPlusy = pos[0] + pos[1];
     Scalar source = -std::exp(xPlusy);
     return (xPlusy < 1.0) ? 2.0 * source : source;
   }
@@ -59,7 +59,8 @@ public:
     return d_;
   }
 
-  //permeability per aperture of the fracture in normal direction at position pos
+  //permeability per aperture of the fracture in normal direction
+  //at position pos
   Scalar Kperp (const Coordinate& pos) const
   {
     return 1.0 / d_;