diff --git a/dune/gfe/rodreader.hh b/dune/gfe/rodreader.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cc61c80d9c6c99aedac21944585b91410575414e
--- /dev/null
+++ b/dune/gfe/rodreader.hh
@@ -0,0 +1,118 @@
+#ifndef ROD_READER_HH
+#define ROD_READER_HH
+
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+
+#include "rigidbodymotion.hh"
+
+class RodReader
+{
+public:
+
+    /** \brief Read a spatial rod from file
+     */
+    static void readRod(std::vector<RigidBodyMotion<double,3> >& rod, 
+            const std::string& filename)
+    {
+
+        // /////////////////////////////////////////////////////
+        // Load the AmiraMesh file
+        // ////////////////////////////////////////////////
+
+        AmiraMesh* am = AmiraMesh::read(filename.c_str());
+
+        if(!am)
+            DUNE_THROW(Dune::IOError, "Could not open AmiraMesh file: " << filename);
+
+        int i, j;
+
+        bool datafound=false;
+
+        // get the translation
+        AmiraMesh::Data* am_ValueData =  am->findData("Nodes", HxFLOAT, 3, "Coordinates");
+        if (am_ValueData) {
+            datafound = true;
+
+            if (rod.size()<am->nElements("Nodes"))
+                DUNE_THROW(Dune::IOError, "When reading data from a surface field the "
+                        << "array you provide has to have at least the size of the surface!");
+
+            float* am_values_float = (float*) am_ValueData->dataPtr();
+
+            for (i=0; i<am->nElements("Nodes"); i++) {
+                for (j=0; j<3; j++)
+                    rod[i].r[j] = am_values_float[i*3+j];
+            }
+
+        } else
+            DUNE_THROW(Dune::IOError, "Couldn't find Coordinates: " << filename);
+
+        // get the rotation
+        Dune::BlockVector<Dune::FieldVector<double,3> > dir0(rod.size()),dir1(rod.size());
+        am_ValueData =  am->findData("Nodes", HxFLOAT, 3, "Directors0");
+        if (am_ValueData) {
+            datafound = true;
+
+            if (rod.size()<am->nElements("Nodes"))
+                DUNE_THROW(Dune::IOError, "When reading data from a surface field the "
+                        << "array you provide has to have at least the size of the surface!");
+
+            float* am_values_float = (float*) am_ValueData->dataPtr();
+
+            for (i=0; i<am->nElements("Nodes"); i++) {
+                for (j=0; j<3; j++)
+                    dir0[i][j] = am_values_float[i*3+j];
+            }
+
+        } else 
+            DUNE_THROW(Dune::IOError, "Couldn't find Directors0: " << filename);
+
+        am_ValueData =  am->findData("Nodes", HxFLOAT, 3, "Directors1");
+        if (am_ValueData) {
+            datafound = true;
+
+            if (rod.size()<am->nElements("Nodes"))
+                DUNE_THROW(Dune::IOError, "When reading data from a surface field the "
+                        << "array you provide has to have at least the size of the surface!");
+
+            float* am_values_float = (float*) am_ValueData->dataPtr();
+
+            for (i=0; i<am->nElements("Nodes"); i++) {
+                for (j=0; j<3; j++)
+                    dir1[i][j] = am_values_float[i*3+j];
+            }
+
+        } else 
+            DUNE_THROW(Dune::IOError, "Couldn't find Directors1: " << filename);
+
+
+        for (i=0; i<rod.size(); i++) {
+            dir0[i] /= dir0[i].two_norm(); 
+            dir1[i] /= dir1[i].two_norm(); 
+            
+            // compute third director as the outer normal of the cross-section
+            Dune::FieldVector<double,3> dir2;
+            dir2[0] = dir0[i][1]*dir1[i][2]-dir0[i][2]*dir1[i][1];    
+            dir2[1] = dir0[i][2]*dir1[i][0]-dir0[i][0]*dir1[i][2];    
+            dir2[2] = dir0[i][0]*dir1[i][1]-dir0[i][1]*dir1[i][0];    
+            dir2 /= dir2.two_norm();
+
+            // the rotation matrix corresponding to the directors
+            Dune::FieldMatrix<double,3,3> rot(0);
+            for (j=0;j<3;j++) {
+                rot[j][0] = dir0[i][j];
+                rot[j][1] = dir1[i][j];
+                rot[j][2] = dir2[j];
+            }
+            rod[i].q.set(rot);
+        }
+
+        std::cout << "Rod successfully read from: " << filename << std::endl;
+
+    }
+};
+#endif