From ac21922a66fb9624640327f8c64c0489e09acc31 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 12 Oct 2007 08:40:59 +0000
Subject: [PATCH] proper derivative of the geodesic interpolation
[[Imported from SVN: r1677]]
---
src/rodassembler.cc | 79 ++++++++++++++++++++++++++++++++++++++++++---
src/rodassembler.hh | 4 +++
2 files changed, 79 insertions(+), 4 deletions(-)
diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index c088ae23..4f7a4347 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -17,10 +17,6 @@ energy(const EntityPointer& element,
{
RT energy = 0;
- // Extract local solution on this element
- const Dune::LagrangeShapeFunctionSet<double, double, 1> & baseSet
- = Dune::LagrangeShapeFunctions<double, double, 1>::general(element->type(), k);
-
// ///////////////////////////////////////////////////////////////////////////////
// The following two loops are a reduced integration scheme. We integrate
// the transverse shear and extensional energy with a first-order quadrature
@@ -75,6 +71,81 @@ energy(const EntityPointer& element,
return energy;
}
+template <class GridType, class RT>
+template <class T>
+void RodLocalStiffness<GridType, RT>::
+interpolationDerivative(const Quaternion<T>& q0, const Quaternion<T>& q1, double s,
+ Dune::array<Quaternion<double>,6>& grad)
+{
+ // Clear output array
+ for (int i=0; i<6; i++)
+ grad[i] = 0;
+
+ // Compute q_0^{-1}
+ Quaternion<T> q0Inv = q0;
+ q0Inv.invert();
+
+ // Compute v = s \exp^{-1} ( q_0^{-1} q_1)
+ Dune::FieldVector<T,3> v = Quaternion<T>::expInv(q0Inv.mult(q1));
+ v *= s;
+
+ Dune::FieldMatrix<T,4,3> dExp_v = Quaternion<T>::Dexp(v);
+
+ Dune::FieldMatrix<T,3,4> dExpInv = Quaternion<T>::DexpInv(q0Inv.mult(q1));
+
+ /** \todo Compute this once and for all */
+ Dune::FieldMatrix<T,4,3> dExp_v_0 = Quaternion<T>::Dexp(Dune::FieldVector<T,3>(0));
+
+
+
+ Dune::FieldMatrix<T,4,4> mat(0);
+ for (int i=0; i<4; i++)
+ for (int j=0; j<4; j++)
+ for (int k=0; k<3; k++)
+ mat[i][j] += s * dExp_v[i][k] * dExpInv[k][j];
+
+
+ // The derivatives with respect to w^0
+ for (int i=0; i<3; i++) {
+
+ Quaternion<T> dw;
+ for (int j=0; j<4; j++)
+ dw[j] = dExp_v_0[j][i];
+
+ // First addend
+ Quaternion<T> grad0 = q0.mult(dw.mult(Quaternion<T>::exp(v)));
+
+ // Second addend
+ dw.conjugate();
+
+ dw[3] -= 2 * dExp_v_0[3][i];
+
+ dw = dw.mult(q0Inv.mult(q1));
+
+ mat.umv(dw,grad[i]);
+ grad[i] = q0.mult(grad[i]);
+
+ // Add the two addends
+ for (int j=0; j<4; j++)
+ grad[i][j] = grad0[j] + grad[i][j];
+
+ }
+
+ // The derivatives with respect to w^1
+ for (int i=3; i<6; i++) {
+
+ Quaternion<T> dw;
+ for (int j=0; j<4; j++)
+ dw[j] = dExp_v_0[j][i-3];
+
+ dw = q0Inv.mult(q1.mult(dw));
+
+ mat.umv(dw,grad[i]);
+ grad[i] = q0.mult(grad[i]);
+
+ }
+}
+
template <class GridType, class RT>
Dune::FieldVector<double, 6> RodLocalStiffness<GridType, RT>::
getStrain(const std::vector<Configuration>& localSolution,
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index eaecc939..da2246ca 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -92,6 +92,10 @@ public:
const std::vector<Configuration>& localReferenceConfiguration,
int k=1);
+ template <class T>
+ static void interpolationDerivative(const Quaternion<T>& q0, const Quaternion<T>& q1, double s,
+ Dune::array<Quaternion<double>,6>& grad);
+
Dune::FieldVector<double, 6> getStrain(const std::vector<Configuration>& localSolution,
const EntityPointer& element,
const Dune::FieldVector<double,1>& pos) const;
--
GitLab