Correct the iterative calculation of the polar decompositon and its derivative

Before, the calculation was terminated after exactly three iterations, now the
stopping criterion is fulfilled if the change in the norm of the polar factor is
small enough.

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;