From 207ce563d8369ea15ddcde921ee1476dc2a7ac2d Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 22 Oct 2019 14:52:45 +0200
Subject: [PATCH] Move the computeDR method to the Rotation class
Because it really is a functionality of rotation matrices in quaternion
representation. Only the Cosserat energies use it at the moment, but
this is coincidental.
In the same process, the method is renamed to quaternionTangentToMatrixTangent.
That seems to be a more telling name.
---
dune/gfe/cosseratenergystiffness.hh | 67 +-----------------------
dune/gfe/nonplanarcosseratshellenergy.hh | 33 +-----------
dune/gfe/rigidbodymotion.hh | 37 +++++++++++++
dune/gfe/rotation.hh | 37 +++++++++++++
test/rigidbodymotiontest.cc | 12 ++---
5 files changed, 82 insertions(+), 104 deletions(-)
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index ead0fa54..aafdca5a 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -121,67 +121,6 @@ class CosseratEnergyLocalStiffness
return result;
}
-public: // for testing
- /** \brief Compute the derivative of the rotation (with respect to x), but wrt matrix coordinates
- \param value Value of the gfe function at a certain point
- \param derivative First derivative of the gfe function wrt x at that point, in quaternion coordinates
- \param DR First derivative of the gfe function wrt x at that point, in matrix coordinates
- */
- static void computeDR(const RigidBodyMotion<field_type,3>& value,
- const Dune::FieldMatrix<field_type,7,gridDim>& derivative,
- Tensor3<field_type,3,3,gridDim>& DR)
- {
- // The LocalGFEFunction class gives us the derivatives of the orientation variable,
- // but as a map into quaternion space. To obtain matrix coordinates we use the
- // chain rule, which means that we have to multiply the given derivative with
- // the derivative of the embedding of the unit quaternion into the space of 3x3 matrices.
- // This second derivative is almost given by the method getFirstDerivativesOfDirectors.
- // However, 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
- //
- // So, if I am not mistaken, DR[i][j][k] contains \partial R_ij / \partial k
- Tensor3<field_type,3 , 3, 4> dd_dq;
- value.q.getFirstDerivativesOfDirectors(dd_dq);
-
- DR = field_type(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] * derivative[l+3][k];
-
- }
-
- /** \brief Compute the derivative of the rotation (with respect to x), but wrt matrix coordinates
- \param value Value of the gfe function at a certain point
- \param derivative First derivative of the gfe function wrt x at that point, in quaternion coordinates
- \param DR First derivative of the gfe function wrt x at that point, in matrix coordinates
- */
- static void computeDR(const Rotation<field_type,3>& value,
- const Dune::FieldMatrix<field_type,4,gridDim>& derivative,
- Tensor3<field_type,3,3,gridDim>& DR)
- {
- // The LocalGFEFunction class gives us the derivatives of the orientation variable,
- // but as a map into quaternion space. To obtain matrix coordinates we use the
- // chain rule, which means that we have to multiply the given derivative with
- // the derivative of the embedding of the unit quaternion into the space of 3x3 matrices.
- // This second derivative is almost given by the method getFirstDerivativesOfDirectors.
- // However, 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
- //
- // So, if I am not mistaken, DR[i][j][k] contains \partial R_ij / \partial k
- Tensor3<field_type,3 , 3, 4> dd_dq;
- value.getFirstDerivativesOfDirectors(dd_dq);
-
- DR = field_type(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] * derivative[l][k];
-
- }
-
public:
/** \brief Constructor with a set of material parameters
@@ -432,8 +371,7 @@ energy(const typename Basis::LocalView& localView,
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
- Tensor3<field_type,3,3,gridDim> DR;
- computeDR(value, derivative, DR);
+ Tensor3<field_type,3,3,gridDim> DR = value.quaternionTangentToMatrixTangent(derivative);
// Add the local energy density
if (gridDim==2) {
@@ -598,8 +536,7 @@ energy(const typename Basis::LocalView& localView,
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
- Tensor3<field_type,3,3,gridDim> DR;
- computeDR(orientationValue, orientationDerivative, DR);
+ Tensor3<field_type,3,3,gridDim> DR = orientationValue.quaternionTangentToMatrixTangent(orientationDerivative);
// Add the local energy density
if (gridDim==2) {
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index c1d09880..16c97c4f 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -147,36 +147,6 @@ class NonplanarCosseratShellEnergy
return result;
}
-
- /** \brief Compute the derivative of the rotation (with respect to x), but wrt matrix coordinates
- \param value Value of the gfe function at a certain point
- \param derivative First derivative of the gfe function wrt x at that point, in quaternion coordinates
- \param DR First derivative of the gfe function wrt x at that point, in matrix coordinates
- */
- static void computeDR(const RigidBodyMotion<field_type,3>& value,
- const Dune::FieldMatrix<field_type,7,gridDim>& derivative,
- Tensor3<field_type,3,3,gridDim>& DR)
- {
- // The LocalGFEFunction class gives us the derivatives of the orientation variable,
- // but as a map into quaternion space. To obtain matrix coordinates we use the
- // chain rule, which means that we have to multiply the given derivative with
- // the derivative of the embedding of the unit quaternion into the space of 3x3 matrices.
- // This second derivative is almost given by the method getFirstDerivativesOfDirectors.
- // However, 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
- //
- // So, DR[i][j][k] contains \partial R_ij / \partial k
- Tensor3<field_type,3 , 3, 4> dd_dq;
- value.q.getFirstDerivativesOfDirectors(dd_dq);
-
- DR = field_type(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] * derivative[l+3][k];
- }
-
public:
/** \brief Constructor with a set of material parameters
@@ -329,8 +299,7 @@ energy(const typename Basis::LocalView& localView,
value.q.matrix(R);
auto RT = transpose(R);
- Tensor3<field_type,3,3,gridDim> DR;
- computeDR(value, derivative, DR);
+ Tensor3<field_type,3,3,gridDim> DR = value.quaternionTangentToMatrixTangent(derivative);
//////////////////////////////////////////////////////////
// Fundamental forms and curvature
diff --git a/dune/gfe/rigidbodymotion.hh b/dune/gfe/rigidbodymotion.hh
index a68e796f..b4541e3e 100644
--- a/dune/gfe/rigidbodymotion.hh
+++ b/dune/gfe/rigidbodymotion.hh
@@ -284,6 +284,43 @@ public:
return result;
}
+ /** \brief Transform tangent vectors from quaternion representation to matrix representation
+ *
+ * This class represents rotations as unit quaternions, and therefore variations of rotations
+ * (i.e., tangent vector) are represented in quaternion space, too. However, some applications
+ * require the tangent vectors as matrices. To obtain matrix coordinates we use the
+ * chain rule, which means that we have to multiply the given derivative with
+ * the derivative of the embedding of the unit quaternion into the space of 3x3 matrices.
+ * This second derivative is almost given by the method getFirstDerivativesOfDirectors.
+ * However, 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
+ *
+ * As a typical GFE assembler will require this for several tangent vectors at once,
+ * the implementation here is a vector one: It allows to treat several tangent vectors
+ * together, which have to come as the columns of the `derivative` parameter.
+ *
+ * So, if I am not mistaken, result[i][j][k] contains \partial R_ij / \partial k
+ *
+ * \param derivative Tangent vector in quaternion coordinates
+ * \returns DR Tangent vector in matrix coordinates
+ */
+ template <int blocksize>
+ Tensor3<T,3,3,blocksize> quaternionTangentToMatrixTangent(const Dune::FieldMatrix<T,7,blocksize>& derivative) const
+ {
+ Tensor3<T,3,3,blocksize> result = T(0);
+
+ Tensor3<T,3 , 3, 4> dd_dq;
+ q.getFirstDerivativesOfDirectors(dd_dq);
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<blocksize; k++)
+ for (int l=0; l<4; l++)
+ result[i][j][k] += dd_dq[j][i][l] * derivative[l+3][k];
+
+ return result;
+ }
+
/** \brief The Weingarten map */
EmbeddedTangentVector weingarten(const EmbeddedTangentVector& z, const EmbeddedTangentVector& v) const {
diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh
index e767e147..2ca609ed 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -671,6 +671,43 @@ public:
}
+ /** \brief Transform tangent vectors from quaternion representation to matrix representation
+ *
+ * This class represents rotations as unit quaternions, and therefore variations of rotations
+ * (i.e., tangent vector) are represented in quaternion space, too. However, some applications
+ * require the tangent vectors as matrices. To obtain matrix coordinates we use the
+ * chain rule, which means that we have to multiply the given derivative with
+ * the derivative of the embedding of the unit quaternion into the space of 3x3 matrices.
+ * This second derivative is almost given by the method getFirstDerivativesOfDirectors.
+ * However, 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
+ *
+ * As a typical GFE assembler will require this for several tangent vectors at once,
+ * the implementation here is a vector one: It allows to treat several tangent vectors
+ * together, which have to come as the columns of the `derivative` parameter.
+ *
+ * So, if I am not mistaken, result[i][j][k] contains \partial R_ij / \partial k
+ *
+ * \param derivative Tangent vector in quaternion coordinates
+ * \returns DR Tangent vector in matrix coordinates
+ */
+ template <int blocksize>
+ Tensor3<T,3,3,blocksize> quaternionTangentToMatrixTangent(const Dune::FieldMatrix<T,4,blocksize>& derivative) const
+ {
+ Tensor3<T,3,3,blocksize> result = T(0);
+
+ Tensor3<T,3 , 3, 4> dd_dq;
+ getFirstDerivativesOfDirectors(dd_dq);
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<blocksize; k++)
+ for (int l=0; l<4; l++)
+ result[i][j][k] += dd_dq[j][i][l] * derivative[l][k];
+
+ return result;
+ }
+
/** \brief Compute the derivative of the squared distance function with respect to the second argument
*
* The squared distance function is
diff --git a/test/rigidbodymotiontest.cc b/test/rigidbodymotiontest.cc
index 58d9f5c1..2f39109d 100644
--- a/test/rigidbodymotiontest.cc
+++ b/test/rigidbodymotiontest.cc
@@ -1,14 +1,14 @@
#include "config.h"
+#include <dune/common/parallel/mpihelper.hh>
+
#include <dune/geometry/type.hh>
#include <dune/geometry/quadraturerules.hh>
-#include <dune/grid/uggrid.hh>
-
-#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/cosseratenergystiffness.hh>
#include "multiindex.hh"
#include "valuefactory.hh"
@@ -75,9 +75,7 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
// evaluate actual derivative
Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::dimension, domainDim> derivative = f.evaluateDerivative(quadPos);
- Tensor3<double,3,3,domainDim> DR;
- typedef Dune::Functions::LagrangeBasis<typename UGGrid<domainDim>::LeafGridView,1> FEBasis;
- CosseratEnergyLocalStiffness<FEBasis,3>::computeDR(f.evaluate(quadPos),derivative, DR);
+ Tensor3<double,3,3,domainDim> DR = f.evaluate(quadPos).quaternionTangentToMatrixTangent(derivative);
// evaluate fd approximation of derivative
Tensor3<double,3,3,domainDim> DR_fd = evaluateDerivativeFD(f,quadPos);
--
GitLab