diff --git a/Makefile.am b/Makefile.am
index 2c0fae94bdefff84c4c0869bf043b7d74e554558..684a8ea4499f2c563ae4a8b5b303f0dec3a37105 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -4,11 +4,14 @@
LDADD = $(IPOPT_LDFLAGS) $(IPOPT_LIBS)
AM_CPPFLAGS += $(IPOPT_CPPFLAGS)
-noinst_PROGRAMS = staticrod staticrod2 rod3d harmonicmaps dirneucoupling
+noinst_PROGRAMS = staticrod staticrod2 rod3d harmonicmaps dirneucoupling rod-eoc
staticrod_SOURCES = staticrod.cc
staticrod2_SOURCES = staticrod2.cc
rod3d_SOURCES = rod3d.cc
+
+rod_eoc_SOURCES = rod-eoc.cc
+
harmonicmaps_SOURCES = harmonicmaps.cc
harmonicmaps_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $(PSURFACE_CPPFLAGS)
harmonicmaps_LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \
diff --git a/rod-eoc.cc b/rod-eoc.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e7f0f0f808a20da6220e1495e6f4f1ef0b27e708
--- /dev/null
+++ b/rod-eoc.cc
@@ -0,0 +1,207 @@
+#include <config.h>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/configparser.hh>
+
+#include <dune/grid/onedgrid.hh>
+
+#include <dune/istl/io.hh>
+
+
+#include <dune-solvers/solvers/iterativesolver.hh>
+#include <dune-solvers/norms/energynorm.hh>
+
+#include "src/rigidbodymotion.hh"
+#include "src/geodesicdifference.hh"
+#include "src/rotation.hh"
+#include "src/rodassembler.hh"
+#include "src/riemanniantrsolver.hh"
+#include "src/rodrefine.hh"
+
+typedef Dune::OneDGrid GridType;
+
+typedef RigidBodyMotion<3> TargetSpace;
+typedef std::vector<RigidBodyMotion<3> > SolutionType;
+
+const int blocksize = TargetSpace::TangentVector::size;
+
+using namespace Dune;
+using std::string;
+
+void solve (const GridType& grid,
+ SolutionType& x, int numLevels,
+ const TargetSpace& dirichletValue,
+ ConfigParser& parameters)
+{
+ // read solver setting
+ const double innerTolerance = parameters.get<double>("innerTolerance");
+ const double tolerance = parameters.get<double>("tolerance");
+ const int maxTrustRegionSteps = parameters.get<int>("maxTrustRegionSteps");
+ const double initialTrustRegionRadius = parameters.get<double>("initialTrustRegionRadius");
+ const int multigridIterations = parameters.get<int>("numIt");
+
+ // read rod parameter settings
+ const double A = parameters.get<double>("A");
+ const double J1 = parameters.get<double>("J1");
+ const double J2 = parameters.get<double>("J2");
+ const double E = parameters.get<double>("E");
+ const double nu = parameters.get<double>("nu");
+
+ // Create a local assembler
+ RodLocalStiffness<OneDGrid::LeafGridView,double> localStiffness(grid.leafView(),
+ A, J1, J2, E, nu);
+
+
+
+
+ x.resize(grid.size(1));
+
+ // //////////////////////////
+ // Initial solution
+ // //////////////////////////
+
+ for (int i=0; i<x.size(); i++) {
+ 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();
+ }
+
+ // set Dirichlet value
+ x.back() = dirichletValue;
+
+ // Both ends are Dirichlet
+ BitSetVector<blocksize> dirichletNodes(grid.size(1));
+ dirichletNodes.unsetAll();
+
+ dirichletNodes[0] = dirichletNodes.back() = true;
+
+ // ///////////////////////////////////////////
+ // Create a solver for the rod problem
+ // ///////////////////////////////////////////
+
+ RodAssembler<GridType> rodAssembler(grid, &localStiffness);
+
+ RiemannianTrustRegionSolver<GridType,RigidBodyMotion<3> > rodSolver;
+#if 0
+ rodSolver.setup(grid,
+ &rodAssembler,
+ x,
+ dirichletNodes,
+ tolerance,
+ maxTrustRegionSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ innerTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ false); // instrumentation
+#else
+ rodSolver.setupTCG(grid,
+ &rodAssembler,
+ x,
+ dirichletNodes,
+ tolerance,
+ maxTrustRegionSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ innerTolerance,
+ false); // instrumentation
+#endif
+
+ // /////////////////////////////////////////////////////
+ // Solve!
+ // /////////////////////////////////////////////////////
+
+ rodSolver.setInitialSolution(x);
+ rodSolver.solve();
+
+ x = rodSolver.getSol();
+}
+
+int main (int argc, char *argv[]) try
+{
+ // parse data file
+ ConfigParser parameterSet;
+ if (argc==2)
+ parameterSet.parseFile(argv[1]);
+ else
+ parameterSet.parseFile("rod-eoc.parset");
+
+ // read solver settings
+ const int numLevels = parameterSet.get<int>("numLevels");
+ const int nu1 = parameterSet.get<int>("nu1");
+ const int nu2 = parameterSet.get<int>("nu2");
+ const int mu = parameterSet.get<int>("mu");
+ const int baseIterations = parameterSet.get<int>("baseIt");
+ const double baseTolerance = parameterSet.get<double>("baseTolerance");
+
+ const int numRodBaseElements = parameterSet.get<int>("numRodBaseElements");
+
+ // /////////////////////////////////////////
+ // Read Dirichlet values
+ // /////////////////////////////////////////
+
+ RigidBodyMotion<3> dirichletValue;
+ dirichletValue.r[0] = parameterSet.get<double>("dirichletValueX");
+ dirichletValue.r[1] = parameterSet.get<double>("dirichletValueY");
+ dirichletValue.r[2] = parameterSet.get<double>("dirichletValueZ");
+
+ FieldVector<double,3> axis;
+ axis[0] = parameterSet.get<double>("dirichletAxisX");
+ axis[1] = parameterSet.get<double>("dirichletAxisY");
+ axis[2] = parameterSet.get<double>("dirichletAxisZ");
+ double angle = parameterSet.get<double>("dirichletAngle");
+
+ dirichletValue.q = Rotation<3,double>(axis, M_PI*angle/180);
+
+ // ///////////////////////////////////////////////////////////
+ // First compute the 'exact' solution on a very fine grid
+ // ///////////////////////////////////////////////////////////
+
+ // Create the reference grid
+
+ GridType referenceGrid(numRodBaseElements, 0, 1);
+
+ referenceGrid.globalRefine(numLevels-1);
+
+ // Solve the rod Dirichlet problem
+ SolutionType referenceSolution;
+ solve(referenceGrid, referenceSolution, numLevels, dirichletValue, parameterSet);
+
+
+
+ // ///////////////////////////////////////////////////////////
+ // Compute on all coarser levels, and compare
+ // ///////////////////////////////////////////////////////////
+
+
+
+ for (int i=1; i<=numLevels; i++) {
+
+ GridType grid(numRodBaseElements, 0, 1);
+ grid.globalRefine(i-1);
+
+ // compute again
+ SolutionType solution;
+ solve(grid, solution, i, dirichletValue, parameterSet);
+
+ // Prolong solution to the very finest grid
+ for (int j=i; j<numLevels; j++)
+ globalRodRefine(grid, solution);
+
+ assert(referenceSolution.size() == solution.size());
+
+ // Compute max-norm difference
+ std::cout << "Level: " << i-1
+ << ", max-norm error: " << computeGeodesicDifference(solution,referenceSolution).infinity_norm()
+ << std::endl;
+
+ }
+
+ } catch (Exception e) {
+
+ std::cout << e << std::endl;
+
+ }