diff --git a/dirneucoupling.cc b/dirneucoupling.cc
index a6bca1bb34f09f653ab11a95def6148fb9ebfa2e..d4e086f5a72212fd5e7c83236692695b87786d98 100644
--- a/dirneucoupling.cc
+++ b/dirneucoupling.cc
@@ -52,7 +52,7 @@ using std::vector;
 // Some types that I need
 //typedef BCRSMatrix<FieldMatrix<double, dim, dim> > OperatorType;
 //typedef BlockVector<FieldVector<double, dim> >     VectorType;
-typedef vector<RigidBodyMotion<dim> >              RodSolutionType;
+typedef vector<RigidBodyMotion<double,dim> >       RodSolutionType;
 typedef BlockVector<FieldVector<double, 6> >       RodDifferenceType;
 
 
@@ -163,7 +163,7 @@ int main (int argc, char *argv[]) try
     //   Determine rod Dirichlet values
     // /////////////////////////////////////////
     rodFactory.create(complex.rods_["rod"].dirichletValues_,
-                      RigidBodyMotion<3>(FieldVector<double,3>(0), Rotation<3,double>::identity()));
+                      RigidBodyMotion<double,3>(FieldVector<double,3>(0), Rotation<double,3>::identity()));
     BitSetVector<1> rodDNodes(complex.rods_["rod"].dirichletValues_.size(), false);
     
     // we need at least one Dirichlet side
@@ -171,11 +171,11 @@ int main (int argc, char *argv[]) try
 
     if (parameterSet.hasKey("dirichletValue0")){
         
-        RigidBodyMotion<3> dirichletValue;
+        RigidBodyMotion<double,3> dirichletValue;
         dirichletValue.r = parameterSet.get("dirichletValue0", rodRestEndPoint[0]);
         FieldVector<double,3> axis = parameterSet.get("dirichletAxis0", FieldVector<double,3>(0));
         double angle = parameterSet.get("dirichletAngle0", double(0));
-        dirichletValue.q = Rotation<3,double>(axis, M_PI*angle/180);
+        dirichletValue.q = Rotation<double,3>(axis, M_PI*angle/180);
         
         rodX[0] = dirichletValue;
 
@@ -187,11 +187,11 @@ int main (int argc, char *argv[]) try
 
     if (parameterSet.hasKey("dirichletValue1")) {
         
-        RigidBodyMotion<3> dirichletValue;
+        RigidBodyMotion<double,3> dirichletValue;
         dirichletValue.r = parameterSet.get("dirichletValue1", rodRestEndPoint[1]);
         FieldVector<double,3> axis = parameterSet.get("dirichletAxis1", FieldVector<double,3>(0));
         double angle = parameterSet.get("dirichletAngle1", double(0));
-        dirichletValue.q = Rotation<3,double>(axis, M_PI*angle/180);
+        dirichletValue.q = Rotation<double,3>(axis, M_PI*angle/180);
         
         rodX.back() = dirichletValue;
     
@@ -302,7 +302,7 @@ int main (int argc, char *argv[]) try
 
     RodAssembler<RodGridType::LeafGridView,3> rodAssembler(complex.rods_["rod"].grid_->leafView(), &rodLocalStiffness);
 
-    RiemannianTrustRegionSolver<RodGridType,RigidBodyMotion<3> > rodSolver;
+    RiemannianTrustRegionSolver<RodGridType,RigidBodyMotion<double,3> > rodSolver;
     rodSolver.setup(*complex.rods_["rod"].grid_, 
                     &rodAssembler,
                     rodX,
@@ -379,14 +379,14 @@ int main (int argc, char *argv[]) try
     //   Dirichlet-Neumann Solver
     // /////////////////////////////////////////////////////
 
-    RigidBodyMotion<3> referenceInterface = (parameterSet.hasKey("dirichletValue1"))
+    RigidBodyMotion<double,3> referenceInterface = (parameterSet.hasKey("dirichletValue1"))
                                           ? rodX[0]
                                           : rodX.back();
                                           
     complex.couplings_[interfaceName].referenceInterface_ = referenceInterface;
 
     // Init interface value
-    std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> > lambda;
+    std::map<std::pair<std::string,std::string>, RigidBodyMotion<double,3> > lambda;
     lambda[interfaceName] = referenceInterface;
     
     ///////////////////////////////////////////////////////////////////////////7
@@ -412,7 +412,7 @@ int main (int argc, char *argv[]) try
         std::cout << "----------------------------------------------------" << std::endl;
         
         // Backup of the current iterate for the error computation later on
-        std::map<std::pair<std::string,std::string>, RigidBodyMotion<3> > oldLambda  = lambda;
+        std::map<std::pair<std::string,std::string>, RigidBodyMotion<double,3> > oldLambda  = lambda;
         
         if (ddType=="FixedPointIteration") {
 
@@ -481,7 +481,7 @@ int main (int argc, char *argv[]) try
         //   Compute error in the energy norm
         // ////////////////////////////////////////////
 
-        double lengthOfCorrection = RigidBodyMotion<3>::distance(oldLambda[interfaceName], lambda[interfaceName]);
+        double lengthOfCorrection = RigidBodyMotion<double,3>::distance(oldLambda[interfaceName], lambda[interfaceName]);
 
         double convRate = lengthOfCorrection / normOfOldCorrection;
 
diff --git a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
index c08341540157f0d841dd6bae65e230a683e56f26..bf877835d269defee857acafc8712f4db8304bc7 100644
--- a/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
+++ b/dune/gfe/coupling/rodcontinuumsteklovpoincarestep.hh
@@ -503,7 +503,7 @@ continuumDirichletToNeumannMap(const std::string& continuumName,
     //  Extract the residual stresses
     //////////////////////////////////////////////////////////////////////////////
             
-    std::map<std::pair<std::string,std::string>, RigidBodyMotion<3>::TangentVector > result;
+    std::map<std::pair<std::string,std::string>, RigidBodyMotion<double,3>::TangentVector > result;
 
     for (it = lambda.begin(); it!=lambda.end(); ++it) {
         
diff --git a/harmonicmaps-eoc.cc b/harmonicmaps-eoc.cc
index b2c3896cfa5d151604183ea35d9b15b8f555e29b..6c5d5d6c0fe91e9e52ef443beb1247a691331dc9 100644
--- a/harmonicmaps-eoc.cc
+++ b/harmonicmaps-eoc.cc
@@ -35,7 +35,7 @@
 // grid dimension
 const int dim = 3;
 
-typedef UnitVector<3> TargetSpace;
+typedef UnitVector<double,3> TargetSpace;
 typedef std::vector<TargetSpace> SolutionType;
 
 const int blocksize = TargetSpace::TangentVector::dimension;
diff --git a/rod-eoc.cc b/rod-eoc.cc
index f6a9f28e2fc41316e274c0b0a2116ba752ba878c..bc4aa11e3d6ce1e196073031b198669606bfb531 100644
--- a/rod-eoc.cc
+++ b/rod-eoc.cc
@@ -26,8 +26,8 @@
 
 typedef Dune::OneDGrid GridType;
 
-typedef RigidBodyMotion<3> TargetSpace;
-typedef std::vector<RigidBodyMotion<3> > SolutionType;
+typedef RigidBodyMotion<double,3> TargetSpace;
+typedef std::vector<RigidBodyMotion<double,3> > SolutionType;
 
 const int blocksize = TargetSpace::TangentVector::dimension;
 
@@ -71,7 +71,7 @@ void solve (const GridType& grid,
         x[i].r[0] = 0;
         x[i].r[1] = 0;
         x[i].r[2] = double(i)/(x.size()-1);
-        x[i].q    = Rotation<3,double>::identity();
+        x[i].q    = Rotation<double,3>::identity();
     }
 
     //  set Dirichlet value
@@ -89,7 +89,7 @@ void solve (const GridType& grid,
 
     RodAssembler<GridType::LeafGridView,3> rodAssembler(grid.leafView(), &localStiffness);
 
-    RiemannianTrustRegionSolver<GridType,RigidBodyMotion<3> > rodSolver;
+    RiemannianTrustRegionSolver<GridType,RigidBodyMotion<double,3> > rodSolver;
 #if 1
     rodSolver.setup(grid, 
                     &rodAssembler,
@@ -150,7 +150,7 @@ int main (int argc, char *argv[]) try
     //   Read Dirichlet values
     // /////////////////////////////////////////
 
-    RigidBodyMotion<3> dirichletValue;
+    RigidBodyMotion<double,3> dirichletValue;
     dirichletValue.r[0] = parameterSet.get<double>("dirichletValueX");
     dirichletValue.r[1] = parameterSet.get<double>("dirichletValueY");
     dirichletValue.r[2] = parameterSet.get<double>("dirichletValueZ");
@@ -161,7 +161,7 @@ int main (int argc, char *argv[]) try
     axis[2] = parameterSet.get<double>("dirichletAxisZ");
     double angle = parameterSet.get<double>("dirichletAngle");
 
-    dirichletValue.q = Rotation<3,double>(axis, M_PI*angle/180);
+    dirichletValue.q = Rotation<double,3>(axis, M_PI*angle/180);
 
     // ///////////////////////////////////////////////////////////
     //   First compute the 'exact' solution on a very fine grid