From 264d819beb154a215e87394fa61b6e763c920bc5 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 4 Jan 2007 14:48:58 +0000
Subject: [PATCH] moved interpolation of derivatives to the quaternion class
[[Imported from SVN: r1116]]
---
src/quaternion.hh | 44 ++++++++++++++++++++++++++++++++++++++++++++
src/rodassembler.cc | 18 ++++++------------
2 files changed, 50 insertions(+), 12 deletions(-)
diff --git a/src/quaternion.hh b/src/quaternion.hh
index dec08c50..211acbc6 100644
--- a/src/quaternion.hh
+++ b/src/quaternion.hh
@@ -121,9 +121,20 @@ public:
return result;
}
+ /** \brief Invert the quaternion */
+ void invert() {
+
+ (*this)[0] *= -1;
+ (*this)[1] *= -1;
+ (*this)[2] *= -1;
+
+ (*this) /= this->two_norm2();
+ }
+
/** \brief Interpolate between two rotations */
static Quaternion<T> interpolate(const Quaternion<T>& a, const Quaternion<T>& b, double omega) {
+#if 0 // Interpolation in H plus normalization
Quaternion<T> result;
for (int i=0; i<4; i++)
@@ -131,6 +142,39 @@ public:
result.normalize();
+ return result;
+#else // Interpolation on the geodesic in SO(3) from a to b
+
+ Quaternion<T> diff = a;
+ diff.invert();
+ diff.mult(b);
+
+ T dist = 2*std::acos(diff[3]);
+
+ T invSinc = (dist < 1e-4) ? 1/(0.5+(dist*dist/48)) : dist / std::sin(dist/2);
+
+ // Compute difference on T_a SO(3)
+ Dune::FieldVector<double,3> v;
+ v[0] = diff[0] * invSinc;
+ v[1] = diff[1] * invSinc;
+ v[2] = diff[2] * invSinc;
+
+ v *= omega;
+
+ return a.mult(exp(v[0], v[1], v[2]));
+#endif
+
+ }
+
+ /** \brief Interpolate between two rotations */
+ static Quaternion<T> interpolateDerivative(const Quaternion<T>& a, const Quaternion<T>& b,
+ double omega, double intervallLength) {
+
+ Quaternion<T> result;
+
+ for (int i=0; i<4; i++)
+ result[i] = a[i] / (-intervallLength) + b[i] / intervallLength;
+
return result;
}
diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 8f15337e..1cd3c083 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -303,9 +303,8 @@ getLocalMatrix( EntityPointer &entity,
}
// Get the derivative of the rotation at the quadrature point by interpolating in $H$
- Quaternion<double> hatq_s;
- for (int i=0; i<4; i++)
- hatq_s[i] = localSolution[0].q[i]*shapeGrad[0] + localSolution[1].q[i]*shapeGrad[1];
+ Quaternion<double> hatq_s = Quaternion<double>::interpolateDerivative(localSolution[0].q, localSolution[1].q,
+ quadPos, 1/shapeGrad[1]);
// The current strain
FieldVector<double,blocksize> strain = getStrain(globalSolution, entity, quadPos);
@@ -402,7 +401,6 @@ getLocalMatrix( EntityPointer &entity,
// \partial W^2 \partial v^i_j \partial v^k_l
// All other derivatives are zero
- //double sum = duLocal_dvij[k][l][m] * (duCan_dvij[i][j] * hatq.director(m) + darbouxCan*dd_dvj[m][j]*shapeFunction[i]);
double sum = du_dvij[k][l][m] * du_dvij[i][j][m];
sum += (strain[m+3] - referenceStrain[m+3]) * du_dvij_dvkl[i][j][k][l][m];
@@ -501,9 +499,8 @@ assembleGradient(const std::vector<Configuration>& sol,
Quaternion<double> hatq = Quaternion<double>::interpolate(localSolution[0].q, localSolution[1].q,quadPos[0]);
// Get the derivative of the rotation at the quadrature point by interpolating in $H$
- Quaternion<double> hatq_s;
- for (int i=0; i<4; i++)
- hatq_s[i] = localSolution[0].q[i]*shapeGrad[0] + localSolution[1].q[i]*shapeGrad[1];
+ Quaternion<double> hatq_s = Quaternion<double>::interpolateDerivative(localSolution[0].q, localSolution[1].q,
+ quadPos, 1/shapeGrad[1]);
// The current strain
FieldVector<double,blocksize> strain = getStrain(sol, it, quadPos);
@@ -811,11 +808,8 @@ Dune::FieldVector<double, 6> Dune::RodAssembler<GridType>::getStrain(const std::
Quaternion<double> q = Quaternion<double>::interpolate(localSolution[0].q, localSolution[1].q, pos);
// Get the derivative of the rotation at the quadrature point by interpolating in $H$
- Quaternion<double> q_s;
- q_s[0] = localSolution[0].q[0]*shapeGrad[0][0] + localSolution[1].q[0]*shapeGrad[1][0];
- q_s[1] = localSolution[0].q[1]*shapeGrad[0][0] + localSolution[1].q[1]*shapeGrad[1][0];
- q_s[2] = localSolution[0].q[2]*shapeGrad[0][0] + localSolution[1].q[2]*shapeGrad[1][0];
- q_s[3] = localSolution[0].q[3]*shapeGrad[0][0] + localSolution[1].q[3]*shapeGrad[1][0];
+ Quaternion<double> q_s = Quaternion<double>::interpolateDerivative(localSolution[0].q, localSolution[1].q,
+ pos, 1/shapeGrad[1]);
// /////////////////////////////////////////////
// Sum it all up
--
GitLab