diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index d74091c095554195da330ad5b5a242751dc38f39..5adf4b8c0d6edce6c74f03c61b9e5c43dc78ffe6 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,4 +1,5 @@
-set(programs cosserat-continuum
+set(programs compute-disc-error
+ cosserat-continuum
harmonicmaps
mixed-cosserat-continuum
rod-eoc
diff --git a/src/compute-disc-error.cc b/src/compute-disc-error.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f08b7cb97ac6e785f07beaf981fa3db6bbc3d2e7
--- /dev/null
+++ b/src/compute-disc-error.cc
@@ -0,0 +1,260 @@
+#include <config.h>
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/basisinterpolator.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
+#include <dune/fufem/discretizationerror.hh>
+#include <dune/fufem/dunepython.hh>
+
+#include <dune/gfe/rotation.hh>
+#include <dune/gfe/unitvector.hh>
+#include <dune/gfe/realtuple.hh>
+#include <dune/gfe/embeddedglobalgfefunction.hh>
+#include <dune/gfe/cosseratvtkwriter.hh>
+
+// grid dimension
+const int dim = 2;
+const int dimworld = 2;
+
+// Image space of the geodesic fe functions
+const int targetDim = 3;
+// typedef Rotation<double,targetDim> TargetSpace;
+// typedef UnitVector<double,targetDim> TargetSpace;
+// typedef RealTuple<double,targetDim> TargetSpace;
+using TargetSpace = RigidBodyMotion<double,targetDim>;
+
+// Tangent vector of the image space
+const int blocksize = TargetSpace::TangentVector::dimension;
+
+using namespace Dune;
+
+template <class GridView, int order>
+void measureEOC(const GridView gridView,
+ const GridView referenceGridView,
+ const ParameterTree& parameterSet)
+{
+ typedef std::vector<TargetSpace> SolutionType;
+
+ //////////////////////////////////////////////////////////////////////////////////
+ // Construct the scalar function space bases corresponding to the GFE space
+ //////////////////////////////////////////////////////////////////////////////////
+
+ typedef Dune::Functions::PQkNodalBasis<GridView, order> FEBasis;
+ FEBasis feBasis(gridView);
+ FEBasis referenceFEBasis(referenceGridView);
+
+ //////////////////////////////////////////////////////////////////////////////////
+ // Read the data whose error is to be measured, and the reference data
+ //////////////////////////////////////////////////////////////////////////////////
+
+ // Input data
+ typedef BlockVector<TargetSpace::CoordinateType> EmbeddedVectorType;
+
+ EmbeddedVectorType embeddedX(feBasis.size());
+ std::ifstream inFile(parameterSet.get<std::string>("simulationData"), std::ios_base::binary);
+ if (not inFile)
+ DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("simulationData") << " could not be opened.");
+ GenericVector::readBinary(inFile, embeddedX);
+ inFile.close();
+
+ SolutionType x(embeddedX.size());
+ for (size_t i=0; i<x.size(); i++)
+ x[i] = TargetSpace(embeddedX[i]);
+
+ // Reference configuration
+ EmbeddedVectorType embeddedReferenceX(referenceFEBasis.size());
+ inFile.open(parameterSet.get<std::string>("referenceData"), std::ios_base::binary);
+ if (not inFile)
+ DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("referenceData") << " could not be opened.");
+ GenericVector::readBinary(inFile, embeddedReferenceX);
+
+ SolutionType referenceX(embeddedReferenceX.size());
+ for (size_t i=0; i<referenceX.size(); i++)
+ referenceX[i] = TargetSpace(embeddedReferenceX[i]);
+
+ //CosseratVTKWriter<typename GridView::Grid>::template write<FEBasis>(feBasis,x, "cosserat_infile_" + std::to_string(gridView.grid().maxLevel()+1));
+
+
+ /////////////////////////////////////////////////////////////////
+ // Measure the discretization error
+ /////////////////////////////////////////////////////////////////
+#if 0
+ if (parameterSet.get<std::string>("discretizationErrorMode")=="analytical")
+ {
+ // Read reference solution and its derivative into a PythonFunction
+ typedef VirtualDifferentiableFunction<FieldVector<double, dim>, TargetSpace::CoordinateType> FBase;
+
+ Python::Module module = Python::import(parameterSet.get<std::string>("referenceSolution"));
+ auto referenceSolution = module.get("fdf").toC<std::shared_ptr<FBase>>();
+
+ // The numerical solution, as a grid function
+ GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(feBasis, x);
+
+ // QuadratureRule for the integral of the L^2 error
+ QuadratureRuleKey quadKey(dim,6);
+
+ // Compute the embedded L^2 error
+ double l2Error = DiscretizationError<GridType::LeafGridView>::computeL2Error(&numericalSolution,
+ referenceSolution.get(),
+ quadKey);
+
+ // Compute the embedded H^1 error
+ double h1Error = DiscretizationError<GridType::LeafGridView>::computeH1HalfNormDifferenceSquared(grid->leafGridView(),
+ &numericalSolution,
+ referenceSolution.get(),
+ quadKey);
+
+ std::cout << "elements: " << grid->leafGridView().size(0)
+ << " "
+ << "L^2 error: " << l2Error
+ << " ";
+ std::cout << "H^1 error: " << std::sqrt(l2Error*l2Error + h1Error) << std::endl;
+ }
+#endif
+
+ if (parameterSet.get<std::string>("discretizationErrorMode")=="discrete")
+ {
+ std::cout << "FOO" << std::endl;
+ using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
+ FufemFEBasis fufemReferenceFEBasis(referenceFEBasis);
+ FufemFEBasis fufemFEBasis(feBasis);
+
+ // The numerical solution, as a grid function
+ GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> referenceSolution(fufemReferenceFEBasis, referenceX);
+
+ // The numerical solution, as a grid function
+ GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(fufemFEBasis, x);
+
+ double l2ErrorSquared = 0;
+ double h1ErrorSquared = 0;
+
+ HierarchicSearch<typename GridView::Grid,typename GridView::IndexSet> hierarchicSearch(gridView.grid(), gridView.indexSet());
+
+ std::cout << "l2: " << l2ErrorSquared << std::endl;
+ for (const auto& rElement : elements(referenceGridView))
+ {
+ const auto& quadRule = QuadratureRules<double, dim>::rule(rElement.type(), 6);
+
+ for (const auto& qp : quadRule)
+ {
+ auto integrationElement = rElement.geometry().integrationElement(qp.position());
+
+ auto globalPos = rElement.geometry().global(qp.position());
+
+ auto element = hierarchicSearch.findEntity(globalPos);
+ auto localPos = element.geometry().local(globalPos);
+
+ auto diff = referenceSolution(rElement, qp.position()) - numericalSolution(element, localPos);
+
+ l2ErrorSquared += integrationElement * qp.weight() * diff.two_norm2();
+
+ auto derDiff = referenceSolution.derivative(rElement, qp.position()) - numericalSolution.derivative(element, localPos);
+
+ h1ErrorSquared += integrationElement * qp.weight() * derDiff.frobenius_norm2();
+
+ }
+ }
+ std::cout << "l2: " << l2ErrorSquared << std::endl;
+
+ std::cout << "levels: " << gridView.grid().maxLevel()+1
+ << " "
+ << "L^2 error: " << std::sqrt(l2ErrorSquared)
+ << " "
+ << "H^1 error: " << std::sqrt(l2ErrorSquared + h1ErrorSquared)
+ << std::endl;
+ }
+}
+
+
+int main (int argc, char *argv[]) try
+{
+ // Start Python interpreter
+ Python::start();
+ Python::Reference main = Python::import("__main__");
+ Python::run("import math");
+
+ Python::runStream()
+ << std::endl << "import sys"
+ << std::endl << "sys.path.append('/home/sander/dune/dune-gfe/problems')"
+ << std::endl;
+
+ // parse data file
+ ParameterTree parameterSet;
+ if (argc < 2)
+ DUNE_THROW(Exception, "Usage: ./compute-disc-error <parameter file>");
+
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
+
+ ParameterTreeParser::readOptions(argc, argv, parameterSet);
+
+ // Print all parameters, to have them in the log file
+ parameterSet.report();
+
+ /////////////////////////////////////////
+ // Create the grids
+ /////////////////////////////////////////
+ typedef UGGrid<dim> GridType;
+
+ const int numLevels = parameterSet.get<int>("numLevels");
+
+ shared_ptr<GridType> grid, referenceGrid;
+
+ FieldVector<double,dimworld> lower(0), upper(1);
+
+ if (parameterSet.get<bool>("structuredGrid"))
+ {
+ lower = parameterSet.get<FieldVector<double,dimworld> >("lower");
+ upper = parameterSet.get<FieldVector<double,dimworld> >("upper");
+
+ array<unsigned int,dim> elements = parameterSet.get<array<unsigned int,dim> >("elements");
+ grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+ referenceGrid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+ }
+ else
+ {
+ std::string path = parameterSet.get<std::string>("path");
+ std::string gridFile = parameterSet.get<std::string>("gridFile");
+ grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+ referenceGrid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+ }
+
+ grid->globalRefine(numLevels-1);
+ referenceGrid->globalRefine(parameterSet.get<int>("numReferenceLevels")-1);
+
+ // Do the actual measurement
+ const int order = parameterSet.get<int>("order");
+ switch (order)
+ {
+ case 1:
+ measureEOC<GridType::LeafGridView,1>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
+ break;
+
+ case 2:
+ measureEOC<GridType::LeafGridView,2>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
+ break;
+
+ case 3:
+ measureEOC<GridType::LeafGridView,3>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
+ break;
+
+ default:
+ DUNE_THROW(NotImplemented, "Order '" << order << "' is not implemented");
+ }
+
+ return 0;
+}
+catch (Exception e)
+{
+ std::cout << e << std::endl;
+ return 1;
+}