diff --git a/Makefile.am b/Makefile.am
index ea8fba908986b09fb6ab846d92619a5f4ff5f383..1507982c6c504c60e2fa436d5e6c2afe238e4038 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -9,7 +9,7 @@ SUBDIRS = m4 src test
LDADD = $(IPOPT_LDFLAGS) $(IPOPT_LIBS)
AM_CPPFLAGS += $(IPOPT_CPPFLAGS)
-noinst_PROGRAMS = rodobstacle rod3d harmonicmaps dirneucoupling rod-eoc
+noinst_PROGRAMS = rodobstacle rod3d harmonicmaps harmonicmaps-eoc dirneucoupling rod-eoc
rodobstacle_SOURCES = rodobstacle.cc
rod3d_SOURCES = rod3d.cc
@@ -21,6 +21,11 @@ harmonicmaps_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $
harmonicmaps_LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \
$(IPOPT_LDFLAGS) $(IPOPT_LIBS) $(PSURFACE_LDFLAGS) $(PSURFACE_LIBS)
+harmonicmaps_eoc_SOURCES = harmonicmaps-eoc.cc
+harmonicmaps_eoc_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $(PSURFACE_CPPFLAGS)
+harmonicmaps_eoc_LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \
+ $(IPOPT_LDFLAGS) $(IPOPT_LIBS) $(PSURFACE_LDFLAGS) $(PSURFACE_LIBS)
+
dirneucoupling_SOURCES = dirneucoupling.cc
dirneucoupling_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $(PSURFACE_CPPFLAGS)
dirneucoupling_LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \
diff --git a/harmonicmaps-eoc.cc b/harmonicmaps-eoc.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5d73e84ed5e0ad91f9316a33f203f92c1679df44
--- /dev/null
+++ b/harmonicmaps-eoc.cc
@@ -0,0 +1,253 @@
+#include <config.h>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/configparser.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/onedgrid.hh>
+#include <dune/grid/../../doc/grids/gridfactory/structuredgridfactory.hh>
+
+#include <dune/ag-common/functionspacebases/p1nodalbasis.hh>
+#include <dune/ag-common/assemblers/operatorassembler.hh>
+#include <dune/ag-common/assemblers/localassemblers/laplaceassembler.hh>
+#include <dune/ag-common/assemblers/localassemblers/massassembler.hh>
+
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+
+#include "src/unitvector.hh"
+#include "src/harmonicenergystiffness.hh"
+#include "src/geodesicfeassembler.hh"
+#include "src/riemanniantrsolver.hh"
+#include "src/rodrefine.hh"
+#include "src/rodwriter.hh"
+
+// grid dimension
+const int dim = 2;
+
+typedef UnitVector<3> TargetSpace;
+typedef std::vector<TargetSpace> SolutionType;
+
+const int blocksize = TargetSpace::TangentVector::size;
+
+using namespace Dune;
+using std::string;
+
+template <class GridType>
+void solve (const shared_ptr<GridType>& grid,
+ SolutionType& x,
+ int numLevels,
+ 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 Dirichlet values
+ // /////////////////////////////////////////
+
+ BitSetVector<1> allNodes(grid->size(dim));
+ allNodes.setAll();
+ LeafBoundaryPatch<GridType> dirichletBoundary(*grid, allNodes);
+
+ BitSetVector<blocksize> dirichletNodes(grid->size(dim));
+ for (int i=0; i<dirichletNodes.size(); i++)
+ dirichletNodes[i] = dirichletBoundary.containsVertex(i);
+
+ // //////////////////////////
+ // Initial solution
+ // //////////////////////////
+
+ x.resize(grid->size(dim));
+
+ FieldVector<double,3> yAxis(0);
+ yAxis[1] = 1;
+
+ typename GridType::LeafGridView::template Codim<dim>::Iterator vIt = grid->template leafbegin<dim>();
+ typename GridType::LeafGridView::template Codim<dim>::Iterator vEndIt = grid->template leafend<dim>();
+
+ for (; vIt!=vEndIt; ++vIt) {
+ int idx = grid->leafIndexSet().index(*vIt);
+
+ FieldVector<double,3> v;
+ FieldVector<double,2> pos = vIt->geometry().corner(0);
+ FieldVector<double,3> axis;
+ axis[0] = pos[0]; axis[1] = pos[1]; axis[2] = 1;
+ Rotation<3,double> rotation(axis, pos.two_norm()*M_PI*1.5);
+
+ if (dirichletNodes[idx][0]) {
+// FieldMatrix<double,3,3> rMat;
+// rotation.matrix(rMat);
+// v = rMat[2];
+ v[0] = std::sin(pos[0]*M_PI);
+ v[1] = 0;
+ v[2] = std::cos(pos[0]*M_PI);
+ } else {
+ v[0] = 1;
+ v[1] = 0;
+ v[2] = 0;
+ }
+
+ x[idx] = v;
+
+ }
+
+
+ // ////////////////////////////////////////////////////////////
+ // Create an assembler for the Harmonic Energy Functional
+ // ////////////////////////////////////////////////////////////
+
+ HarmonicEnergyLocalStiffness<typename GridType::LeafGridView,TargetSpace> harmonicEnergyLocalStiffness;
+
+ GeodesicFEAssembler<typename GridType::LeafGridView,TargetSpace> assembler(grid->leafView(),
+ &harmonicEnergyLocalStiffness);
+
+ // ///////////////////////////////////////////
+ // Create a solver for the rod problem
+ // ///////////////////////////////////////////
+
+ RiemannianTrustRegionSolver<GridType,TargetSpace> solver;
+
+ solver.setup(*grid,
+ &assembler,
+ x,
+ dirichletNodes,
+ tolerance,
+ maxTrustRegionSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ innerTolerance,
+ 1, 3, 3,
+ 100, // iterations of the base solver
+ 1e-8, // base tolerance
+ false); // instrumentation
+
+ // /////////////////////////////////////////////////////
+ // Solve!
+ // /////////////////////////////////////////////////////
+
+ solver.setInitialSolution(x);
+ solver.solve();
+
+ x = solver.getSol();
+}
+
+int main (int argc, char *argv[]) try
+{
+ // parse data file
+ ConfigParser parameterSet;
+ if (argc==2)
+ parameterSet.parseFile(argv[1]);
+ else
+ parameterSet.parseFile("harmonicmaps-eoc.parset");
+
+ // read solver settings
+ const int numLevels = parameterSet.get<int>("numLevels");
+ const int baseIterations = parameterSet.get<int>("baseIt");
+ const double baseTolerance = parameterSet.get<double>("baseTolerance");
+
+ const int numBaseElements = parameterSet.get<int>("numBaseElements");
+
+ // /////////////////////////////////////////
+ // Read Dirichlet values
+ // /////////////////////////////////////////
+
+
+
+
+ // ///////////////////////////////////////////////////////////
+ // First compute the 'exact' solution on a very fine grid
+ // ///////////////////////////////////////////////////////////
+
+ typedef std::conditional<dim==1,OneDGrid,UGGrid<dim> >::type GridType;
+
+ // Create the reference grid
+
+ array<unsigned int,dim> elements;
+ elements.fill(4);
+ shared_ptr<GridType> referenceGrid = StructuredGridFactory<GridType>::createSimplexGrid(FieldVector<double,dim>(0),
+ FieldVector<double,dim>(1),
+ elements);
+ referenceGrid->globalRefine(numLevels-1);
+
+ // Solve the rod Dirichlet problem
+ SolutionType referenceSolution;
+ solve(referenceGrid, referenceSolution, numLevels, parameterSet);
+
+
+ // //////////////////////////////////////////////////////////////////////
+ // Compute mass matrix and laplace matrix to emulate L2 and H1 norms
+ // //////////////////////////////////////////////////////////////////////
+#if 0
+ typedef P1NodalBasis<GridType::LeafGridView,double> FEBasis;
+ FEBasis basis(referenceGrid->leafView());
+ OperatorAssembler<FEBasis,FEBasis> operatorAssembler(basis, basis);
+
+ LaplaceAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> laplaceLocalAssembler;
+ MassAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> massMatrixLocalAssembler;
+
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
+ ScalarMatrixType laplace, massMatrix;
+
+ operatorAssembler.assemble(laplaceLocalAssembler, laplace);
+ operatorAssembler.assemble(massMatrixLocalAssembler, massMatrix);
+#endif
+ // ///////////////////////////////////////////////////////////
+ // Compute on all coarser levels, and compare
+ // ///////////////////////////////////////////////////////////
+
+ for (int i=1; i<=numLevels; i++) {
+
+ array<unsigned int,dim> elements;
+ elements.fill(numBaseElements);
+ shared_ptr<GridType> grid = StructuredGridFactory<GridType>::createSimplexGrid(FieldVector<double,dim>(0),
+ FieldVector<double,dim>(1),
+ elements);
+
+ grid->globalRefine(i-1);
+
+ // compute again
+ SolutionType solution;
+ solve(grid, solution, i, parameterSet);
+
+#if 0
+ // Prolong solution to the very finest grid
+ for (int j=i; j<numLevels; j++)
+ globalRodRefine(grid, solution);
+
+ std::stringstream numberAsAscii;
+ numberAsAscii << i;
+ writeRod(solution, "rodGrid_" + numberAsAscii.str());
+
+ assert(referenceSolution.size() == solution.size());
+#endif
+#if 0
+ BlockVector<TargetSpace::TangentVector> difference = computeGeodesicDifference(solution,referenceSolution);
+
+ H1SemiNorm< BlockVector<TargetSpace::TangentVector> > h1Norm(laplace);
+ H1SemiNorm< BlockVector<TargetSpace::TangentVector> > l2Norm(massMatrix);
+
+ // Compute max-norm difference
+ std::cout << "Level: " << i-1
+ << ", max-norm error: " << difference.infinity_norm()
+ << std::endl;
+
+ std::cout << "Level: " << i-1
+ << ", L2 error: " << l2Norm(difference)
+ << std::endl;
+
+ std::cout << "Level: " << i-1
+ << ", H1 error: " << h1Norm(difference)
+ << std::endl;
+#endif
+ }
+
+ } catch (Exception e) {
+
+ std::cout << e << std::endl;
+
+ }
diff --git a/harmonicmaps-eoc.parset b/harmonicmaps-eoc.parset
new file mode 100644
index 0000000000000000000000000000000000000000..ce10259d114217dc5b3c176b83f63099b8d008fe
--- /dev/null
+++ b/harmonicmaps-eoc.parset
@@ -0,0 +1,31 @@
+# Number of grid levels
+numLevels = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-12
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 200
+
+# Initial trust-region radius
+initialTrustRegionRadius = 1
+
+# Number of multigrid iterations per trust-region step
+numIt = 200
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the inner solver
+innerTolerance = 1e-5
+
+# Tolerance of the base grid solver
+baseTolerance = -1
+
+############################
+# Problem specifications
+############################
+
+# Number of elements of the rod grid
+numBaseElements = 4
+