diff --git a/test/Makefile.am b/test/Makefile.am
index 9cd1709286580cdc76f30b717a33f2bd6b42c8fc..d11177b97a7d31685d75940215baee9e1ea1b489 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -5,7 +5,7 @@ LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS)
 AM_CPPFLAGS += $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) -Wall
 
 check_PROGRAMS = frameinvariancetest rotationtest fdcheck localgeodesicfefunctiontest \
-                 harmonicenergytest
+                 harmonicenergytest averagedistanceassemblertest
 
 frameinvariancetest_SOURCES = frameinvariancetest.cc
 
@@ -17,6 +17,8 @@ localgeodesicfefunctiontest_SOURCES = localgeodesicfefunctiontest.cc
 
 harmonicenergytest_SOURCES = harmonicenergytest.cc
 
+averagedistanceassemblertest_SOURCES = averagedistanceassemblertest.cc
+
 # don't follow the full GNU-standard
 # we need automake 1.5
 AUTOMAKE_OPTIONS = foreign 1.5
diff --git a/test/averagedistanceassemblertest.cc b/test/averagedistanceassemblertest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..62c4f1d796eb04f13dbbf6d3fcb5bf19f9d9279c
--- /dev/null
+++ b/test/averagedistanceassemblertest.cc
@@ -0,0 +1,155 @@
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/fvector.hh>
+#include <dune/grid/common/quadraturerules.hh>
+
+#include <dune/src/rotation.hh>
+#include <dune/src/realtuple.hh>
+#include <dune/src/unitvector.hh>
+
+#include <dune/src/localgeodesicfefunction.hh>
+
+// Domain dimension
+const int dim = 2;
+
+using namespace Dune;
+
+/** \brief Test whether interpolation is invariant under permutation of the simplex vertices
+ */
+template <class TargetSpace>
+void testPoint(const std::vector<TargetSpace>& corners, 
+               const std::vector<double>& weights,
+               const TargetSpace& argument)
+{
+    // create the assembler
+    AverageDistanceAssembler<TargetSpace> assembler(corners, weights);
+
+    // test the functional
+    double value = assembler.value(argument);
+    assert(!std::isnan(value));
+    assert(value >= 0);
+    
+    // test the gradient
+    typename TargetSpace::TangentVector gradient;
+    assembler.assembleGradient(argument, gradient);
+
+    // test the hessian
+    FieldMatrix<double, TargetSpace::TangentVector::size, TargetSpace::TangentVector::size> hessian;
+    FieldMatrix<double, TargetSpace::TangentVector::size, TargetSpace::TangentVector::size> hessianApproximation;
+
+    assembler.assembleHessian(argument, hessian);
+    assembler.assembleHessianApproximation(argument, hessianApproximation);
+
+    for (size_t i=0; i<hessian.N(); i++)
+        for (size_t j=0; j<hessian.M(); j++) {
+            assert(!std::isnan(hessian[i][j]));
+            assert(!std::isnan(hessianApproximation[i][j]));
+        }
+
+    FieldMatrix<double, TargetSpace::TangentVector::size, TargetSpace::TangentVector::size> diff = hessian;
+    diff -= hessianApproximation;
+    std::cout << "Matrix inf diff: " << diff.infinity_norm() << std::endl;
+}
+
+
+template <class TargetSpace>
+void testWeightSet(const std::vector<TargetSpace>& corners, 
+                   const TargetSpace& argument)
+{
+    // A quadrature rule as a set of test points
+    int quadOrder = 3;
+    
+    const Dune::QuadratureRule<double, dim>& quad 
+        = Dune::QuadratureRules<double, dim>::rule(GeometryType(GeometryType::simplex,dim), quadOrder);
+    
+    for (size_t pt=0; pt<quad.size(); pt++) {
+        
+        const Dune::FieldVector<double,dim>& quadPos = quad[pt].position();
+        
+        // local to barycentric coordinates
+        std::vector<double> weights(dim+1);
+        weights[0] = 1;
+        for (int i=0; i<dim; i++) {
+            weights[0] -= quadPos[i];
+            weights[i+1] = quadPos[i];
+        }
+
+        testPoint(corners, weights, argument);
+
+    }
+
+}
+
+
+void testRealTuples()
+{
+    typedef RealTuple<1> TargetSpace;
+
+    std::vector<TargetSpace> corners = {TargetSpace(1),
+                                        TargetSpace(2),
+                                        TargetSpace(3)};
+
+    TargetSpace argument = corners[0];
+    testWeightSet(corners, argument);
+    argument = corners[1];
+    testWeightSet(corners, argument);
+    argument = corners[2];
+    testWeightSet(corners, argument);
+}
+
+void testUnitVectors()
+{
+    typedef UnitVector<3> TargetSpace;
+
+    std::vector<TargetSpace> corners(dim+1);
+
+    FieldVector<double,3> input;
+    input[0] = 1;  input[1] = 0;  input[2] = 0;
+    corners[0] = input;
+    input[0] = 0;  input[1] = 1;  input[2] = 0;
+    corners[1] = input;
+    input[0] = 0;  input[1] = 0;  input[2] = 1;
+    corners[2] = input;
+
+    TargetSpace argument = corners[0];
+    testWeightSet(corners, argument);
+    argument = corners[1];
+    testWeightSet(corners, argument);
+    argument = corners[2];
+    testWeightSet(corners, argument);
+}
+
+void testRotations()
+{
+    typedef Rotation<3,double> TargetSpace;
+
+    FieldVector<double,3> xAxis(0);
+    xAxis[0] = 1;
+    FieldVector<double,3> yAxis(0);
+    yAxis[1] = 1;
+    FieldVector<double,3> zAxis(0);
+    zAxis[2] = 1;
+
+
+    std::vector<TargetSpace> corners(dim+1);
+    corners[0] = Rotation<3,double>(xAxis,0.1);
+    corners[1] = Rotation<3,double>(yAxis,0.1);
+    corners[2] = Rotation<3,double>(zAxis,0.1);
+
+    TargetSpace argument = corners[0];
+    testWeightSet(corners, argument);
+    argument = corners[1];
+    testWeightSet(corners, argument);
+    argument = corners[2];
+    testWeightSet(corners, argument);
+}
+
+
+int main()
+{
+    testRealTuples();
+    testUnitVectors();
+    //testRotations();
+}