From 3250a985f29f1939a87e21d156da92f501823414 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sat, 22 Jan 2011 11:05:11 +0000
Subject: [PATCH] bugfix: properly implement setRotation

[[Imported from SVN: r6819]]
---
 dune/gfe/averageinterface.hh                  | 21 ++++++++++++-------
 .../rodcontinuumsteklovpoincarestep.hh        |  8 +------
 2 files changed, 15 insertions(+), 14 deletions(-)

diff --git a/dune/gfe/averageinterface.hh b/dune/gfe/averageinterface.hh
index 41e9f637..3897d5c5 100644
--- a/dune/gfe/averageinterface.hh
+++ b/dune/gfe/averageinterface.hh
@@ -817,12 +817,22 @@ void computeAverageInterface(const BoundaryPatchBase<GridView>& interface,
 template <class GridView>
 void setRotation(const BoundaryPatchBase<GridView>& dirichletBoundary,
                  Dune::BlockVector<Dune::FieldVector<double,GridView::dimension> >& deformation,
-                 const RigidBodyMotion<3>& relativeMovement)
+                 const RigidBodyMotion<3>& referenceInterface,
+                 const RigidBodyMotion<3>& lambda)
 {
     const typename GridView::IndexSet& indexSet = dirichletBoundary.gridView().indexSet();
     const int dim        = GridView::dimension;
     const int dimworld   = GridView::dimensionworld;
 
+    // Get the relative rotation, first as a quaternion...
+    Rotation<3,double> relativeRotation;
+    relativeRotation = referenceInterface.q.inverse();
+    relativeRotation = lambda.q.mult(relativeRotation);
+
+    // ... then as a matrix
+    Dune::FieldMatrix<double,3,3> rotation;
+    relativeRotation.matrix(rotation);
+
     // ///////////////////////////////////////////
     //   Loop over all vertices
     // ///////////////////////////////////////////
@@ -840,15 +850,12 @@ void setRotation(const BoundaryPatchBase<GridView>& dirichletBoundary,
             int globalIdx = indexSet.subIndex(*it->inside(), cornerIdx, dim);
 
             // Get vertex position
-            Dune::FieldVector<double,dimworld> pos = it->inside()->geometry().corner(cornerIdx);
+            const Dune::FieldVector<double,dimworld> pos = it->inside()->geometry().corner(cornerIdx);
             
             // Action of the rigid body motion
-            Dune::FieldMatrix<double,3,3> rotation;
-            relativeMovement.q.matrix(rotation);
-            
             Dune::FieldVector<double,dimworld> rpos;
-            rotation.mv(pos, rpos);
-            rpos += relativeMovement.r;
+            rotation.mv(pos-referenceInterface.r, rpos);
+            rpos += lambda.r;
             
             // We compute _displacements_, not positions
             rpos -= pos;
diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index a3beeeba..fe7fd337 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -442,14 +442,8 @@ continuumDirichletToNeumannMap(const RigidBodyMotion<3>& lambda) const
     x3d = 0;
 
     // Turn \lambda \in TSE(3) into a Dirichlet value for the continuum
-    RigidBodyMotion<3> relativeMovement;
-    relativeMovement.r = lambda.r - referenceInterface_.r;
-    relativeMovement.q = referenceInterface_.q;
-    relativeMovement.q.invert();
-    relativeMovement.q = lambda.q.mult(relativeMovement.q);
-
     const LeafBoundaryPatch<ContinuumGridType>& foo = complex_.coupling(couplingName).continuumInterfaceBoundary_;
-    setRotation(foo, x3d, relativeMovement);
+    setRotation(foo, x3d, referenceInterface_, lambda);
     
     // Set the correct Dirichlet nodes
     dynamic_cast<IterationStep<VectorType>* >(solver_->iterationStep_)->ignoreNodes_ = &dirichletAndCouplingNodes_;
-- 
GitLab