diff --git a/Makefile.am b/Makefile.am
index 91b69b52530cd434ed2b769f7a43b19934193d08..2c0fae94bdefff84c4c0869bf043b7d74e554558 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -4,11 +4,15 @@
LDADD = $(IPOPT_LDFLAGS) $(IPOPT_LIBS)
AM_CPPFLAGS += $(IPOPT_CPPFLAGS)
-noinst_PROGRAMS = staticrod staticrod2 rod3d dirneucoupling
+noinst_PROGRAMS = staticrod staticrod2 rod3d harmonicmaps dirneucoupling
staticrod_SOURCES = staticrod.cc
staticrod2_SOURCES = staticrod2.cc
rod3d_SOURCES = rod3d.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) \
+ $(IPOPT_LDFLAGS) $(IPOPT_LIBS) $(PSURFACE_LDFLAGS) $(PSURFACE_LIBS)
dirneucoupling_SOURCES = dirneucoupling.cc
dirneucoupling_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $(PSURFACE_CPPFLAGS)
diff --git a/harmonicmaps.cc b/harmonicmaps.cc
new file mode 100644
index 0000000000000000000000000000000000000000..06356b9b30ecc03182922bb95a9f1fe81b7b59d8
--- /dev/null
+++ b/harmonicmaps.cc
@@ -0,0 +1,215 @@
+#include <config.h>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/configparser.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/amirameshreader.hh>
+
+#include <dune-solvers/solvers/iterativesolver.hh>
+#include <dune-solvers/norms/energynorm.hh>
+
+#include "src/roddifference.hh"
+#include "src/rodwriter.hh"
+#include "src/rotation.hh"
+#include "src/geodesicfeassembler.hh"
+#include "src/rodsolver.hh"
+
+// grid dimension
+const int dim = 3;
+
+// Image space of the geodesic fe functions
+typedef Rotation<3,double> TargetSpace;
+
+// Tangent vector of the image space
+const int blocksize = TargetSpace::TangentVector::size;
+
+using namespace Dune;
+
+int main (int argc, char *argv[]) try
+{
+ typedef std::vector<TargetSpace> SolutionType;
+
+ // parse data file
+ ConfigParser parameterSet;
+ if (argc==2)
+ parameterSet.parseFile(argv[1]);
+ else
+ parameterSet.parseFile("harmonicmaps.parset");
+
+ // read solver settings
+ const int numLevels = parameterSet.get<int>("numLevels");
+ const double tolerance = parameterSet.get<double>("tolerance");
+ const int maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
+ const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+ const int multigridIterations = parameterSet.get<int>("numIt");
+ 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 mgTolerance = parameterSet.get<double>("mgTolerance");
+ const double baseTolerance = parameterSet.get<double>("baseTolerance");
+ const bool instrumented = parameterSet.get<bool>("instrumented");
+ std::string resultPath = parameterSet.get("resultPath", "");
+
+ // read problem settings
+ std::string path = parameterSet.get<std::string>("path");
+ std::string gridFile = parameterSet.get<std::string>("gridFile");
+
+ // ///////////////////////////////////////
+ // Create the grid
+ // ///////////////////////////////////////
+ typedef UGGrid<dim> GridType;
+ GridType grid;
+ grid.setRefinementType(GridType::COPY);
+ AmiraMeshReader<GridType>::read(grid, path + gridFile);
+
+ grid.globalRefine(numLevels-1);
+
+ SolutionType x(grid.size(dim));
+
+ // //////////////////////////
+ // Initial solution
+ // //////////////////////////
+
+ for (int i=0; i<x.size(); i++) {
+ x[i] = Rotation<3,double>::identity();
+ }
+
+ // backup for error measurement later
+ SolutionType initialIterate = x;
+
+ // /////////////////////////////////////////
+ // Read Dirichlet values
+ // /////////////////////////////////////////
+
+ BitSetVector<1> allNodes(grid.size(dim));
+ allNodes.setAll();
+ LeafBoundaryPatch<GridType> dirichletBoundary(grid, allNodes);
+
+ BitSetVector<blocksize> dirichletNodes(grid.size(1));
+ for (int i=0; i<dirichletNodes.size(); i++)
+ dirichletNodes[i] = dirichletBoundary.containsVertex(i);
+
+ // ////////////////////////////////////////////////////////////
+ // Create an assembler for the Harmonic Energy Functional
+ // ////////////////////////////////////////////////////////////
+ GeodesicFEAssembler<GridType::LeafGridView,TargetSpace> assembler(grid.leafView());
+
+ // /////////////////////////////////////////////////
+ // Create a Riemannian trust-region solver
+ // /////////////////////////////////////////////////
+#if 0
+ RodSolver<GridType> rodSolver;
+ rodSolver.setup(grid,
+ &assembler,
+ x,
+ dirichletNodes,
+ tolerance,
+ maxTrustRegionSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance,
+ instrumented);
+
+ // /////////////////////////////////////////////////////
+ // Solve!
+ // /////////////////////////////////////////////////////
+
+ std::cout << "Energy: " << assembler.computeEnergy(x) << std::endl;
+
+ rodSolver.setInitialSolution(x);
+ rodSolver.solve();
+
+ x = rodSolver.getSol();
+#endif
+ // //////////////////////////////
+ // Output result
+ // //////////////////////////////
+#if 0
+ writeRod(x, resultPath + "rod3d.result");
+#endif
+
+#if 0
+ // //////////////////////////////////////////////////////////
+ // Recompute and compare against exact solution
+ // //////////////////////////////////////////////////////////
+
+ SolutionType exactSolution = x;
+
+ // //////////////////////////////////////////////////////////
+ // Compute hessian of the rod functional at the exact solution
+ // for use of the energy norm it creates.
+ // //////////////////////////////////////////////////////////
+
+ BCRSMatrix<FieldMatrix<double, 6, 6> > hessian;
+ MatrixIndexSet indices(exactSolution.size(), exactSolution.size());
+ rodAssembler.getNeighborsPerVertex(indices);
+ indices.exportIdx(hessian);
+ rodAssembler.assembleMatrixFD(exactSolution, hessian);
+
+
+ double error = std::numeric_limits<double>::max();
+
+ SolutionType intermediateSolution(x.size());
+
+ // Create statistics file
+ std::ofstream statisticsFile((resultPath + "trStatistics").c_str());
+
+ // Compute error of the initial iterate
+ typedef BlockVector<FieldVector<double,6> > RodDifferenceType;
+ RodDifferenceType rodDifference = computeRodDifference(exactSolution, initialIterate);
+ double oldError = std::sqrt(computeEnergyNormSquared(rodDifference, hessian));
+
+ int i;
+ for (i=0; i<maxTrustRegionSteps; i++) {
+
+ // /////////////////////////////////////////////////////
+ // Read iteration from file
+ // /////////////////////////////////////////////////////
+ char iSolFilename[100];
+ sprintf(iSolFilename, "tmp/intermediateSolution_%04d", i);
+
+ FILE* fp = fopen(iSolFilename, "rb");
+ if (!fp)
+ DUNE_THROW(IOError, "Couldn't open intermediate solution '" << iSolFilename << "'");
+ for (int j=0; j<intermediateSolution.size(); j++) {
+ fread(&intermediateSolution[j].r, sizeof(double), 3, fp);
+ fread(&intermediateSolution[j].q, sizeof(double), 4, fp);
+ }
+
+ fclose(fp);
+
+ // /////////////////////////////////////////////////////
+ // Compute error
+ // /////////////////////////////////////////////////////
+
+ rodDifference = computeRodDifference(exactSolution, intermediateSolution);
+
+ error = std::sqrt(computeEnergyNormSquared(rodDifference, hessian));
+
+
+ double convRate = error / oldError;
+
+ // Output
+ std::cout << "Trust-region iteration: " << i << " error : " << error << ", "
+ << "convrate " << convRate << std::endl;
+ statisticsFile << i << " " << error << " " << convRate << std::endl;
+
+ if (error < 1e-12)
+ break;
+
+ oldError = error;
+
+ }
+#endif
+
+ // //////////////////////////////
+ } catch (Exception e) {
+
+ std::cout << e << std::endl;
+
+ }