From 54d3bc81afa0d4e349c25f445bb6b96feb9c4e7e Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 15 Jul 2011 09:45:52 +0000
Subject: [PATCH] Move the computation of the derivatives of the rotations in
matrix coordinates into a separate method, and fix an index bug on the way.
Now this derivative passes an fd check.
[[Imported from SVN: r7565]]
---
dune/gfe/cosseratenergystiffness.hh | 50 +++++++++++++++++------------
1 file changed, 30 insertions(+), 20 deletions(-)
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index c179b51d..793256f3 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -68,7 +68,34 @@ class CosseratEnergyLocalStiffness
return result;
}
-
+
+public: // for testing
+ /** \brief Compute the derivative of the rotation, but wrt matrix coordinates */
+ static void computeDR(const RigidBodyMotion<3>& value,
+ const Dune::FieldMatrix<double,7,gridDim>& derivative,
+ Tensor3<double,3,3,3>& DR)
+ {
+ // derivative of the rotation part in quaternion coordinates
+ Dune::FieldMatrix<double,4,gridDim> DR_quat;
+ for (int i=0; i<4; i++)
+ for (int j=0; j<gridDim; j++)
+ DR_quat[i][j] = derivative[i+3][j];
+
+ // first get the derivative of the embedding of H_1 into R^{3\times3}
+ // Since the directors of a given unit quaternion are the _columns_ of the
+ // corresponding orthogonal matrix, we need to invert the i and j indices
+ Tensor3<double,3 , 3, 4> dd_dq;
+ value.q.getFirstDerivativesOfDirectors(dd_dq);
+
+ //
+ DR = 0;
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<gridDim; k++)
+ for (int l=0; l<4; l++)
+ DR[i][j][k] += dd_dq[j][i][l] * DR_quat[l][k];
+
+ }
public:
@@ -235,25 +262,8 @@ energy(const Entity& element,
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
- // derivative of the rotation part in quaternion coordinates
- Dune::FieldMatrix<double,4,gridDim> DR_quat;
- for (int i=0; i<4; i++)
- for (int j=0; j<gridDim; j++)
- DR_quat[i][j] = derivative[i+3][j];
-
- // transform to matrix coordinates:
- // first get the derivative of the embedding of H_1 into R^{3\times3}
- Tensor3<double,3 , 3, 4> dd_dq;
- value.q.getFirstDerivativesOfDirectors(dd_dq);
-
- //
- Tensor3<double,3,3,3> DR(0);
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- for (int k=0; k<gridDim; k++) {
- for (int l=0; l<4; l++)
- DR[i][j][k] += dd_dq[i][j][l] * DR_quat[l][k];
- }
+ Tensor3<double,3,3,3> DR;
+ computeDR(derivative, DR);
// Add the local energy density
if (gridDim==2) {
--
GitLab