From bbad25b92cbe36cdb23036dcab03e5748afc842f Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Thu, 11 Feb 2021 18:39:26 +0100
Subject: [PATCH] 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 | 18 ++++++++++++++++--
1 file changed, 16 insertions(+), 2 deletions(-)
diff --git a/dune/gfe/localprojectedfefunction.hh b/dune/gfe/localprojectedfefunction.hh
index 237b0ab8..80fa763d 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;
--
GitLab