diff --git a/test/Makefile.am b/test/Makefile.am
index fb417acf8ab8bc6e231e26a38a7deeb53d4c65cf..3e2a607090f9df57f7c603514f2f353df91027e3 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -10,9 +10,11 @@ check_PROGRAMS = averagedistanceassemblertest \
                  frameinvariancetest \
                  globalgfetestfunctionbasistest \
                  harmonicenergytest \
+                 interpolationtest \
                  localgeodesicfefunctiontest \
                  localgeodesicfestiffnesstest \
                  localgfetestfunctiontest \
+                 nestednesstest \
                  rotationtest \
                  svdtest \
                  targetspacetest
@@ -29,10 +31,14 @@ localgeodesicfestiffnesstest_SOURCES = localgeodesicfestiffnesstest.cc
 
 localgfetestfunctiontest_SOURCES = localgfetestfunctiontest.cc
 
+nestednesstest_SOURCES = nestednesstest.cc
+
 globalgfetestfunctionbasistest_SOURCES = globalgfetestfunctionbasistest.cc
 
 harmonicenergytest_SOURCES = harmonicenergytest.cc
 
+interpolationtest_SOURCES = interpolationtest.cc
+
 cosseratenergytest_SOURCES = cosseratenergytest.cc
 
 averagedistanceassemblertest_SOURCES = averagedistanceassemblertest.cc
diff --git a/test/averagedistanceassemblertest.cc b/test/averagedistanceassemblertest.cc
index 6edde06a820b657cecb6e4a0334057475d26196d..e1361137dfcc74203976f0071672d16c33e3cacb 100644
--- a/test/averagedistanceassemblertest.cc
+++ b/test/averagedistanceassemblertest.cc
@@ -88,7 +88,7 @@ void testWeightSet(const std::vector<TargetSpace>& corners,
 
 void testRealTuples()
 {
-    typedef RealTuple<1> TargetSpace;
+    typedef RealTuple<double,1> TargetSpace;
 
     std::vector<TargetSpace> corners = {TargetSpace(1),
                                         TargetSpace(2),
@@ -104,7 +104,7 @@ void testRealTuples()
 
 void testUnitVectors()
 {
-    typedef UnitVector<3> TargetSpace;
+    typedef UnitVector<double,3> TargetSpace;
 
     std::vector<TargetSpace> corners(dim+1);
 
@@ -126,7 +126,7 @@ void testUnitVectors()
 
 void testRotations()
 {
-    typedef Rotation<3,double> TargetSpace;
+    typedef Rotation<double,3> TargetSpace;
 
     FieldVector<double,3> xAxis(0);
     xAxis[0] = 1;
@@ -137,9 +137,9 @@ void testRotations()
 
 
     std::vector<TargetSpace> corners(dim+1);
-    corners[0] = Rotation<3,double>(xAxis,0.1);
-    corners[1] = Rotation<3,double>(yAxis,0.1);
-    corners[2] = Rotation<3,double>(zAxis,0.1);
+    corners[0] = Rotation<double,3>(xAxis,0.1);
+    corners[1] = Rotation<double,3>(yAxis,0.1);
+    corners[2] = Rotation<double,3>(zAxis,0.1);
 
     TargetSpace argument = corners[0];
     testWeightSet(corners, argument);
diff --git a/test/cosseratenergytest.cc b/test/cosseratenergytest.cc
index 44e6bad06820adb28c241b9352f6a2aa521524f4..02e8b9b350147358ec9e16a273a3b9d028ae42e8 100644
--- a/test/cosseratenergytest.cc
+++ b/test/cosseratenergytest.cc
@@ -16,7 +16,7 @@ const int dim = 2;
 
 const double eps = 1e-4;
 
-typedef RigidBodyMotion<3> TargetSpace;
+typedef RigidBodyMotion<double,3> TargetSpace;
 
 using namespace Dune;
 
@@ -138,8 +138,8 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
     // rotate the entire configuration
     std::vector<TargetSpace> rotatedCoefficients(coefficients.size());
     
-    std::vector<Rotation<3> > testRotations;
-    ValueFactory<Rotation<3> >::get(testRotations);
+    std::vector<Rotation<double,3> > testRotations;
+    ValueFactory<Rotation<double,3> >::get(testRotations);
 
     for (size_t i=0; i<testRotations.size(); i++) {
 
@@ -230,9 +230,9 @@ int main(int argc, char** argv) try
     const int domainDim = 2;
     std::cout << " --- Testing derivative of rotation matrix, domain dimension: " << domainDim << " ---" << std::endl;
 
-    std::vector<Rotation<3,double> > testPoints;
+    std::vector<Rotation<double,3> > testPoints;
     
-    ValueFactory<Rotation<3,double> >::get(testPoints);
+    ValueFactory<Rotation<double,3> >::get(testPoints);
     
     int nTestPoints = testPoints.size();
 
diff --git a/test/fdcheck.hh b/test/fdcheck.hh
index 342d8752bb9167ada1925e137e81a4e94af50c2d..00ac0cf6893e4c432e238e9385fc8e04526e8e98 100644
--- a/test/fdcheck.hh
+++ b/test/fdcheck.hh
@@ -6,7 +6,7 @@
 
 #define ABORT_ON_ERROR
 
-void infinitesimalVariation(RigidBodyMotion<3>& c, double eps, int i)
+void infinitesimalVariation(RigidBodyMotion<double,3>& c, double eps, int i)
 {
     if (i<3)
         c.r[i] += eps;
@@ -14,12 +14,12 @@ void infinitesimalVariation(RigidBodyMotion<3>& c, double eps, int i)
         Dune::FieldVector<double,3> axial(0);
         axial[i-3] = eps;
         SkewMatrix<double,3> variation(axial);
-        c.q = c.q.mult(Rotation<3,double>::exp(variation));
+        c.q = c.q.mult(Rotation<double,3>::exp(variation));
     }
 }
 
 template <class GridType>
-void strainFD(const std::vector<RigidBodyMotion<3> >& x, 
+void strainFD(const std::vector<RigidBodyMotion<double,3> >& x, 
               double pos,
               Dune::array<Dune::FieldMatrix<double,2,6>, 6>& fdStrainDerivatives,
               const RodAssembler<GridType,3>& assembler) 
@@ -32,8 +32,8 @@ void strainFD(const std::vector<RigidBodyMotion<3> >& x,
     // ///////////////////////////////////////////////////////////
     //   Compute gradient by finite-difference approximation
     // ///////////////////////////////////////////////////////////
-    std::vector<RigidBodyMotion<3> > forwardSolution = x;
-    std::vector<RigidBodyMotion<3> > backwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
     
     for (size_t i=0; i<x.size(); i++) {
         
@@ -64,7 +64,7 @@ void strainFD(const std::vector<RigidBodyMotion<3> >& x,
 
 
 template <class GridType>
-void strain2ndOrderFD(const std::vector<RigidBodyMotion<3> >& x, 
+void strain2ndOrderFD(const std::vector<RigidBodyMotion<double,3> >& x, 
                       double pos,
                       Dune::array<Dune::Matrix<Dune::FieldMatrix<double,6,6> >, 3>& translationDer,
                       Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3>& rotationDer,
@@ -88,13 +88,13 @@ void strain2ndOrderFD(const std::vector<RigidBodyMotion<3> >& x,
     // ///////////////////////////////////////////////////////////
     //   Compute gradient by finite-difference approximation
     // ///////////////////////////////////////////////////////////
-    std::vector<RigidBodyMotion<3> > forwardSolution = x;
-    std::vector<RigidBodyMotion<3> > backwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
     
-    std::vector<RigidBodyMotion<3> > forwardForwardSolution = x;
-    std::vector<RigidBodyMotion<3> > forwardBackwardSolution = x;
-    std::vector<RigidBodyMotion<3> > backwardForwardSolution = x;
-    std::vector<RigidBodyMotion<3> > backwardBackwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > forwardForwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > forwardBackwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > backwardForwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > backwardBackwardSolution = x;
 
     for (int i=0; i<2; i++) {
         
@@ -173,7 +173,7 @@ void strain2ndOrderFD(const std::vector<RigidBodyMotion<3> >& x,
 
 
 template <class GridType>
-void strain2ndOrderFD2(const std::vector<RigidBodyMotion<3> >& x, 
+void strain2ndOrderFD2(const std::vector<RigidBodyMotion<double,3> >& x, 
                        double pos,
                        Dune::FieldVector<double,1> shapeGrad[2],
                        Dune::FieldVector<double,1> shapeFunction[2],
@@ -197,8 +197,8 @@ void strain2ndOrderFD2(const std::vector<RigidBodyMotion<3> >& x,
     // ///////////////////////////////////////////////////////////
     //   Compute gradient by finite-difference approximation
     // ///////////////////////////////////////////////////////////
-    std::vector<RigidBodyMotion<3> > forwardSolution = x;
-    std::vector<RigidBodyMotion<3> > backwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
     
     for (int k=0; k<2; k++) {
         
@@ -268,9 +268,9 @@ void expHessianFD()
                 //                                            - 2*assembler.computeEnergy(x) 
                 //                                            + assembler.computeEnergy(backward)) / (eps*eps);
                 
-                hessian  = Rotation<3,double>::exp(forward);
-                hessian += Rotation<3,double>::exp(backward);
-                hessian.axpy(-2, Rotation<3,double>::exp(0,0,0));
+                hessian  = Rotation<double,3>::exp(forward);
+                hessian += Rotation<double,3>::exp(backward);
+                hessian.axpy(-2, Rotation<double,3>::exp(0,0,0));
                 hessian /= eps*eps;
                 
             } else {
@@ -288,10 +288,10 @@ void expHessianFD()
                 backwardBackward[j] -= eps;
                 
                 
-                hessian  = Rotation<3,double>::exp(forwardForward);
-                hessian += Rotation<3,double>::exp(backwardBackward);
-                hessian -= Rotation<3,double>::exp(forwardBackward);
-                hessian -= Rotation<3,double>::exp(backwardForward);
+                hessian  = Rotation<double,3>::exp(forwardForward);
+                hessian += Rotation<double,3>::exp(backwardBackward);
+                hessian -= Rotation<double,3>::exp(forwardBackward);
+                hessian -= Rotation<double,3>::exp(backwardForward);
                 hessian /= 4*eps*eps;
                 
             }
@@ -306,7 +306,7 @@ void expHessianFD()
 
 
 template <class GridType>
-void gradientFDCheck(const std::vector<RigidBodyMotion<3> >& x, 
+void gradientFDCheck(const std::vector<RigidBodyMotion<double,3> >& x, 
                      const Dune::BlockVector<Dune::FieldVector<double,6> >& gradient, 
                      const RodAssembler<GridType,3>& assembler)
 {
@@ -320,8 +320,8 @@ void gradientFDCheck(const std::vector<RigidBodyMotion<3> >& x,
     //   Compute gradient by finite-difference approximation
     // ///////////////////////////////////////////////////////////
 
-    std::vector<RigidBodyMotion<3> > forwardSolution = x;
-    std::vector<RigidBodyMotion<3> > backwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
 
     for (size_t i=0; i<x.size(); i++) {
         
@@ -364,7 +364,7 @@ void gradientFDCheck(const std::vector<RigidBodyMotion<3> >& x,
 
 
 template <class GridType>
-void hessianFDCheck(const std::vector<RigidBodyMotion<3> >& x, 
+void hessianFDCheck(const std::vector<RigidBodyMotion<double,3> >& x, 
                     const Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> >& hessian, 
                     const RodAssembler<GridType,3>& assembler)
 {
@@ -378,13 +378,13 @@ void hessianFDCheck(const std::vector<RigidBodyMotion<3> >& x,
     // ///////////////////////////////////////////////////////////
     //   Compute gradient by finite-difference approximation
     // ///////////////////////////////////////////////////////////
-    std::vector<RigidBodyMotion<3> > forwardSolution = x;
-    std::vector<RigidBodyMotion<3> > backwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > forwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > backwardSolution = x;
 
-    std::vector<RigidBodyMotion<3> > forwardForwardSolution = x;
-    std::vector<RigidBodyMotion<3> > forwardBackwardSolution = x;
-    std::vector<RigidBodyMotion<3> > backwardForwardSolution = x;
-    std::vector<RigidBodyMotion<3> > backwardBackwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > forwardForwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > forwardBackwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > backwardForwardSolution = x;
+    std::vector<RigidBodyMotion<double,3> > backwardBackwardSolution = x;
     
 
     // ///////////////////////////////////////////////////////////////
diff --git a/test/localgeodesicfefunctiontest.cc b/test/localgeodesicfefunctiontest.cc
index ca8bc5d145fdc5bb9035275f15239699497d15c9..94da4b3de1915b3d4c3ef56b34c8df21086aee1f 100644
--- a/test/localgeodesicfefunctiontest.cc
+++ b/test/localgeodesicfefunctiontest.cc
@@ -33,7 +33,7 @@ double diameter(const std::vector<TargetSpace>& v)
 
 
 template <int domainDim>
-void testDerivativeTangentiality(const RealTuple<1>& x,
+void testDerivativeTangentiality(const RealTuple<double,1>& x,
                                  const FieldMatrix<double,1,domainDim>& derivative)
 {
     // By construction, derivatives of RealTuples are always tangent
@@ -41,7 +41,7 @@ void testDerivativeTangentiality(const RealTuple<1>& x,
 
 // the columns of the derivative must be tangential to the manifold
 template <int domainDim, int vectorDim>
-void testDerivativeTangentiality(const UnitVector<vectorDim>& x,
+void testDerivativeTangentiality(const UnitVector<double,vectorDim>& x,
                                  const FieldMatrix<double,vectorDim,domainDim>& derivative)
 {
     for (int i=0; i<domainDim; i++) {
@@ -61,14 +61,14 @@ void testDerivativeTangentiality(const UnitVector<vectorDim>& x,
 
 // the columns of the derivative must be tangential to the manifold
 template <int domainDim, int vectorDim>
-void testDerivativeTangentiality(const Rotation<vectorDim-1,double>& x,
+void testDerivativeTangentiality(const Rotation<double,vectorDim-1>& x,
                                  const FieldMatrix<double,vectorDim,domainDim>& derivative)
 {
 }
 
 // the columns of the derivative must be tangential to the manifold
 template <int domainDim, int vectorDim>
-void testDerivativeTangentiality(const RigidBodyMotion<3,double>& x,
+void testDerivativeTangentiality(const RigidBodyMotion<double,3>& x,
                                  const FieldMatrix<double,vectorDim,domainDim>& derivative)
 {
 }
@@ -318,32 +318,32 @@ int main()
     ////////////////////////////////////////////////////////////////
     element.makeSimplex(1);
     
-    test<RealTuple<1>,1>(element);
-    test<UnitVector<2>,1>(element);
-    test<UnitVector<3>,1>(element);
-    test<Rotation<3,double>,1>(element);
-    test<RigidBodyMotion<3,double>,1>(element);
+    test<RealTuple<double,1>,1>(element);
+    test<UnitVector<double,2>,1>(element);
+    test<UnitVector<double,3>,1>(element);
+    test<Rotation<double,3>,1>(element);
+    test<RigidBodyMotion<double,3>,1>(element);
     
     ////////////////////////////////////////////////////////////////
     //  Test functions on 2d simplex elements
     ////////////////////////////////////////////////////////////////
     element.makeSimplex(2);
 
-    test<RealTuple<1>,2>(element);
-    test<UnitVector<2>,2>(element);
-    test<UnitVector<3>,2>(element);
-    test<Rotation<3,double>,2>(element);
-    test<RigidBodyMotion<3,double>,2>(element);
+    test<RealTuple<double,1>,2>(element);
+    test<UnitVector<double,2>,2>(element);
+    test<UnitVector<double,3>,2>(element);
+    test<Rotation<double,3>,2>(element);
+    test<RigidBodyMotion<double,3>,2>(element);
 
     ////////////////////////////////////////////////////////////////
     //  Test functions on 2d quadrilateral elements
     ////////////////////////////////////////////////////////////////
     element.makeCube(2);
 
-    test<RealTuple<1>,2>(element);
-    test<UnitVector<2>,2>(element);
-    test<UnitVector<3>,2>(element);
-    test<Rotation<3,double>,2>(element);
-    test<RigidBodyMotion<3,double>,2>(element);
+    test<RealTuple<double,1>,2>(element);
+    test<UnitVector<double,2>,2>(element);
+    test<UnitVector<double,3>,2>(element);
+    test<Rotation<double,3>,2>(element);
+    test<RigidBodyMotion<double,3>,2>(element);
 
 }
diff --git a/test/localgfetestfunctiontest.cc b/test/localgfetestfunctiontest.cc
index 3358048af0957a7eb10dbf3c5e4d4dc311da1d79..e29447675f55d61cb716f29be2e0d83b6318f46c 100644
--- a/test/localgfetestfunctiontest.cc
+++ b/test/localgfetestfunctiontest.cc
@@ -73,15 +73,15 @@ int main() try
 
     std::cout << std::setw(15) << std::setprecision(12);
     
-    test<RealTuple<1>, 1>();
+    test<RealTuple<double,1>, 1>();
     
-    test<UnitVector<2>, 1>();
-    test<UnitVector<3>, 1>();
-    test<UnitVector<2>, 2>();
-    test<UnitVector<3>, 2>();
+    test<UnitVector<double,2>, 1>();
+    test<UnitVector<double,3>, 1>();
+    test<UnitVector<double,2>, 2>();
+    test<UnitVector<double,3>, 2>();
         
-    test<Rotation<3,double>, 1>();
-    test<Rotation<3,double>, 2>();
+    test<Rotation<double,3>, 1>();
+    test<Rotation<double,3>, 2>();
 
 } catch (Exception e) {
     std::cout << e << std::endl;
diff --git a/test/targetspacetest.cc b/test/targetspacetest.cc
index 89faa5283128a64e6c3dd3bb3fc093c311ba9aed..558bff29b03488aec81b04698ac5a45eedad4077 100644
--- a/test/targetspacetest.cc
+++ b/test/targetspacetest.cc
@@ -350,13 +350,13 @@ void test()
 
 int main() try
 {
-    test<UnitVector<2> >();
-    test<UnitVector<3> >();
-    test<UnitVector<4> >();
+    test<UnitVector<double,2> >();
+    test<UnitVector<double,3> >();
+    test<UnitVector<double,4> >();
 
-    test<Rotation<3> >();
+    test<Rotation<double,3> >();
     
-    test<RigidBodyMotion<3> >();
+    test<RigidBodyMotion<double,3> >();
     
 } catch (Exception e) {
 
diff --git a/test/valuefactory.hh b/test/valuefactory.hh
index d83ce833976e33a7bfce1258ded75eb35afefb9d..db9aa269ac66a4d6bf80290e57c16ba7838514ce 100644
--- a/test/valuefactory.hh
+++ b/test/valuefactory.hh
@@ -25,10 +25,10 @@ public:
  * This is the specialization for RealTuple<1>
  */
 template <>
-class ValueFactory<RealTuple<1> >
+class ValueFactory<RealTuple<double,1> >
 {
 public:
-    static void get(std::vector<RealTuple<1> >& values) {
+    static void get(std::vector<RealTuple<double,1> >& values) {
      
         int nTestPoints = 5;
         double testPoints[5] = {-3, -1, 0, 2, 4};
@@ -37,7 +37,7 @@ public:
         
         // Set up elements of S^1
         for (int i=0; i<nTestPoints; i++)
-            values[i] = RealTuple<1>(testPoints[i]);
+            values[i] = RealTuple<double,1>(testPoints[i]);
         
     }
     
@@ -48,10 +48,10 @@ public:
  * This is the specialization for RealTuple<3>
  */
 template <>
-class ValueFactory<RealTuple<3> >
+class ValueFactory<RealTuple<double,3> >
 {
 public:
-    static void get(std::vector<RealTuple<3> >& values) {
+    static void get(std::vector<RealTuple<double,3> >& values) {
      
         int nTestPoints = 10;
         double testPoints[10][3] = {{1,0,0}, {0,1,0}, {-0.838114,0.356751,-0.412667},
@@ -67,7 +67,7 @@ public:
             Dune::FieldVector<double,3> w;
             for (int j=0; j<3; j++)
                 w[j] = testPoints[i][j];
-            values[i] = RealTuple<3>(w);
+            values[i] = RealTuple<double,3>(w);
 
         }
         
@@ -80,10 +80,10 @@ public:
  * This is the specialization for UnitVector<2>
  */
 template <>
-class ValueFactory<UnitVector<2> >
+class ValueFactory<UnitVector<double,2> >
 {
 public:
-    static void get(std::vector<UnitVector<2> >& values) {
+    static void get(std::vector<UnitVector<double,2> >& values) {
      
         int nTestPoints = 10;
         double testPoints[10][2] = {{1,0}, {0.5,0.5}, {0,1}, {-0.5,0.5}, {-1,0}, {-0.5,-0.5}, {0,-1}, {0.5,-0.5}, {0.1,1}, {1,.1}};
@@ -94,7 +94,7 @@ public:
         for (int i=0; i<nTestPoints; i++) {
         
             Dune::array<double,2> w = {{testPoints[i][0], testPoints[i][1]}};
-            values[i] = UnitVector<2>(w);
+            values[i] = UnitVector<double,2>(w);
 
         }
         
@@ -108,10 +108,10 @@ public:
  * This is the specialization for UnitVector<3>
  */
 template <>
-class ValueFactory<UnitVector<3> >
+class ValueFactory<UnitVector<double,3> >
 {
 public:
-    static void get(std::vector<UnitVector<3> >& values) {
+    static void get(std::vector<UnitVector<double,3> >& values) {
      
         int nTestPoints = 10;
         double testPoints[10][3] = {{1,0,0}, {0,1,0}, {-0.838114,0.356751,-0.412667},
@@ -125,7 +125,7 @@ public:
         for (int i=0; i<nTestPoints; i++) {
         
             Dune::array<double,3> w = {{testPoints[i][0], testPoints[i][1], testPoints[i][2]}};
-            values[i] = UnitVector<3>(w);
+            values[i] = UnitVector<double,3>(w);
 
         }
         
@@ -139,10 +139,10 @@ public:
  * This is the specialization for UnitVector<4>
  */
 template <>
-class ValueFactory<UnitVector<4> >
+class ValueFactory<UnitVector<double,4> >
 {
 public:
-    static void get(std::vector<UnitVector<4> >& values) {
+    static void get(std::vector<UnitVector<double,4> >& values) {
      
         int nTestPoints = 10;
         double testPoints[10][4] = {{1,0,0,0}, {0,1,0,0}, {-0.838114,0.356751,-0.412667,0.5},
@@ -157,7 +157,7 @@ public:
         for (int i=0; i<nTestPoints; i++) {
         
             Dune::array<double,4> w = {{testPoints[i][0], testPoints[i][1], testPoints[i][2], testPoints[i][3]}};
-            values[i] = UnitVector<4>(w);
+            values[i] = UnitVector<double,4>(w);
 
         }
         
@@ -171,10 +171,10 @@ public:
  * This is the specialization for Rotation<3>
  */
 template <>
-class ValueFactory<Rotation<3,double> >
+class ValueFactory<Rotation<double,3> >
 {
 public:
-    static void get(std::vector<Rotation<3,double> >& values) {
+    static void get(std::vector<Rotation<double,3> >& values) {
      
         int nTestPoints = 10;
         double testPoints[10][4] = {{1,0,0,0}, {0,1,0,0}, {-0.838114,0.356751,-0.412667,0.5},
@@ -189,7 +189,7 @@ public:
         for (int i=0; i<nTestPoints; i++) {
         
             Dune::array<double,4> w = {{testPoints[i][0], testPoints[i][1], testPoints[i][2], testPoints[i][3]}};
-            values[i] = Rotation<3,double>(w);
+            values[i] = Rotation<double,3>(w);
 
         }
         
@@ -202,16 +202,16 @@ public:
  * This is the specialization for RigidBodyMotion<3>
  */
 template <>
-class ValueFactory<RigidBodyMotion<3> >
+class ValueFactory<RigidBodyMotion<double,3> >
 {
 public:
-    static void get(std::vector<RigidBodyMotion<3> >& values) {
+    static void get(std::vector<RigidBodyMotion<double,3> >& values) {
      
-        std::vector<RealTuple<3> > rValues;
-        ValueFactory<RealTuple<3> >::get(rValues);
+        std::vector<RealTuple<double,3> > rValues;
+        ValueFactory<RealTuple<double,3> >::get(rValues);
         
-        std::vector<Rotation<3,double> > qValues;
-        ValueFactory<Rotation<3,double> >::get(qValues);
+        std::vector<Rotation<double,3> > qValues;
+        ValueFactory<Rotation<double,3> >::get(qValues);
                                   
         int nTestPoints = std::min(rValues.size(), qValues.size());
         
@@ -219,7 +219,7 @@ public:
         
         // Set up elements of S^1
         for (int i=0; i<nTestPoints; i++)
-            values[i] = RigidBodyMotion<3>(rValues[i].globalCoordinates(),qValues[i]);
+            values[i] = RigidBodyMotion<double,3>(rValues[i].globalCoordinates(),qValues[i]);
         
     }