diff --git a/dune/gfe/localprojectedfefunction.hh b/dune/gfe/localprojectedfefunction.hh
index 237b0ab88f948e8a19e50c9651cb156c2c86cd58..80fa763d558a2b6aad1e4191018f554caaa8a376 100644
--- a/dune/gfe/localprojectedfefunction.hh
+++ b/dune/gfe/localprojectedfefunction.hh
@@ -176,15 +176,22 @@ namespace Dune {
 
       static FieldMatrix<field_type,3,3> polarFactor(const FieldMatrix<field_type,3,3>& matrix)
       {
+        int maxIterations = 100;
+        double tol = 0.001;
         // Use Higham's method
         auto polar = matrix;
-        for (size_t i=0; i<3; i++)
+        for (size_t i=0; i<maxIterations; i++)
         {
+          auto oldPolar = polar;
           auto polarInvert = polar;
           polarInvert.invert();
           for (size_t j=0; j<polar.N(); j++)
             for (size_t k=0; k<polar.M(); k++)
               polar[j][k] = 0.5 * (polar[j][k] + polarInvert[k][j]);
+          oldPolar -= polar;
+          if (oldPolar.frobenius_norm() < tol) {
+            break;
+          }
         }
 
         return polar;
@@ -208,9 +215,11 @@ namespace Dune {
         polar = A;
 
         // Use Heron's method
-        const size_t maxIterations = 3;
+        const size_t maxIterations = 100;
+        double tol = 0.001;
         for (size_t iteration=0; iteration<maxIterations; iteration++)
         {
+          auto oldPolar = polar;
           auto polarInvert = polar;
           polarInvert.invert();
           for (size_t i=0; i<polar.N(); i++)
@@ -241,6 +250,11 @@ namespace Dune {
               for (int k=0; k<3; k++)
                 for (int l=0; l<3; l++)
                   result[i][j][k][l] = 0.5 * (result[i][j][k][l] - tmp2[i][j][k][l]);
+
+          oldPolar -= polar;
+          if (oldPolar.frobenius_norm() < tol) {
+            break;
+          }
         }
 
         return result;