diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index ead0fa54cd5816b2e022bc61c48b6be73e4ea2aa..aafdca5a6983154dd90f452eb894f506281dca92 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 c1d098809c9b3faee26622d99d1bb2e768f67899..16c97c4f0cd595ab06f7bd550aa61622a4a02ab4 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 a68e796f2277b2ef07acac1dafb2e55b7ed8edd4..b4541e3e722dac0d7ec67098771043d1e2ee863b 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 e767e14776e49837220842138d81b1cdde1876fd..2ca609ed8c98312c084c5a0345824db4b9323134 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 58d9f5c189b6bf19b43953b1a933de822d2b843a..2f39109d96ecc25b3553a176e2334d320b97948c 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);