From 15279bd3c221272e50dcb8f671b6fc1cc49a3702 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 31 Jan 2011 08:36:45 +0000
Subject: [PATCH] bugfix: the linearized rod NtD map gets its Neumann values in
 local coordinates (wrt the directors)

[[Imported from SVN: r6901]]
---
 .../rodcontinuumsteklovpoincarestep.hh        | 22 ++++++++++++++++++-
 1 file changed, 21 insertions(+), 1 deletion(-)

diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index 89eb4f7a..152d592f 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -651,6 +651,26 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
         if (couplingName.first != rodName)
             continue;
         
+        // Transform force from canonical coordinates to coordinates given by lambda
+        Dune::FieldVector<double,3> canonicalForce, canonicalTorque;
+        for (int i=0; i<3; i++) {
+            canonicalForce[i]  = it->second[i];
+            canonicalTorque[i] = it->second[3+i];
+        }
+        
+        Dune::FieldMatrix<double,3,3> orientationMatrix;
+        centerOfTorque.find(couplingName)->second.q.matrix(orientationMatrix);
+        
+        /** \todo Why don't we have to transform the force as well? */
+        Dune::FieldVector<double,3> localForce, localTorque;
+        orientationMatrix.mv(canonicalTorque,localTorque);
+        
+        RigidBodyMotion<3>::TangentVector localForceTorque;
+        for (int i=0; i<3; i++) {
+            localForceTorque[i] = localForce[i];
+            localForceTorque[3+i] = localTorque[i];
+        }
+        
         // Use 'forceTorque' as a Neumann value for the rod
         const LeafBoundaryPatch<RodGridType>& interfaceBoundary = complex_.coupling(couplingName).rodInterfaceBoundary_;
 
@@ -659,7 +679,7 @@ linearizedRodNeumannToDirichletMap(const std::string& rodName,
          */
         for (size_t i=0; i<rhs.size(); i++)
             if (interfaceBoundary.containsVertex(i))
-                rhs[i] += it->second;
+                rhs[i] += localForceTorque;
     }
         
     ///////////////////////////////////////////////////////////
-- 
GitLab