From dad60685a33bf8bc60c6c80257f98637b43ff112 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 27 Aug 2007 15:30:34 +0000
Subject: [PATCH] bugfix: third result matrix of svd needs to be transposed
[[Imported from SVN: r1603]]
---
src/averageinterface.hh | 33 +++++++++++++++++++++++++++++----
1 file changed, 29 insertions(+), 4 deletions(-)
diff --git a/src/averageinterface.hh b/src/averageinterface.hh
index 0b6c64fb..22ba084d 100644
--- a/src/averageinterface.hh
+++ b/src/averageinterface.hh
@@ -348,13 +348,38 @@ void computeAverageInterface(const BoundaryPatch<GridType>& interface,
// Get the rotational part of the deformation gradient by performing a
// polar composition.
- FieldVector<double,dim> W;
- FieldMatrix<double,dim,dim> VT;
+ FieldVector<double,dim> W, sigma;
+ FieldMatrix<double,dim,dim> V, VT, U_, VT_;
- // returns a decomposition U W VT, where U is returned in the first argument
- svdcmp<double,dim,dim>(deformationGradient, W, VT);
+ FieldMatrix<double,dim,dim> deformationGradientBackup = deformationGradient;
+ // returns a decomposition U W VT, where U is returned in the first argument
+ svdcmp<double,dim,dim>(deformationGradient, W, V);
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ VT[i][j] = V[j][i];
+#if 0
+ // Debugging
+ std::cout << "Numerical recipes:" << std::endl;
+ std::cout << deformationGradient << std::endl;
+ std::cout << W << std::endl;
+ std::cout << VT << std::endl;
+#endif
deformationGradient.rightmultiply(VT);
+ std::cout << deformationGradient << std::endl;
+
+#if 0
+ lapackSVD(deformationGradientBackup, U_, sigma, VT_);
+
+ std::cout << "Lapack:" << std::endl;
+ std::cout << U_ << std::endl;
+ std::cout << sigma << std::endl;
+ std::cout << VT_ << std::endl;
+ U_.rightmultiply(VT_);
+ std::cout << U_ << std::endl;
+#endif
+
// deformationGradient now contains the orthogonal part of the polar decomposition
assert( std::abs(1-deformationGradient.determinant()) < 1e-3);
--
GitLab