diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index b00b59e1b3d4e79c3a1c56786ae26ea0b7a9e883..7af60d3251dcf6c3e8f16daae82ba7a26cb2b9c5 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -55,6 +55,15 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
         // /////////////////////////////
         //    Solve !
         // /////////////////////////////
+        // TODO: The GramSchmidtSolver may not be the smartest choice here, because it needs
+        // a matrix that is positive definite on the complement of its kernel.  This can limit
+        // the radius of convergence of the Newton solver.  As an example consider interpolation
+        // between two points on the unit sphere that are more than pi/2 apart.  There is a
+        // well-defined unique shortest geodesic between two such points, but depending on the
+        // weights, the weighted sum of squared distances does not everywhere have a positive definite
+        // second derivative.
+        // Maybe the Newton solver has to be replaced by a trust-region or line-search solver
+        // for such situations.
         Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> basis = x_.orthonormalFrame();
         GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(hesseMatrix, corr, rhs, basis);
 
diff --git a/test/harmonicenergytest.cc b/test/harmonicenergytest.cc
index a6e3b3de4620ddfaf3fc2ff28e0a156955f6c506..21e418a1be84032d4463023e16f720af96de20f5 100644
--- a/test/harmonicenergytest.cc
+++ b/test/harmonicenergytest.cc
@@ -17,6 +17,17 @@ typedef UnitVector<double,3> TargetSpace;
 
 using namespace Dune;
 
+/** \brief Computes the diameter of a set */
+template <class TargetSpace>
+double diameter(const std::vector<TargetSpace>& v)
+{
+    double d = 0;
+    for (size_t i=0; i<v.size(); i++)
+        for (size_t j=0; j<v.size(); j++)
+            d = std::max(d, TargetSpace::distance(v[i],v[j]));
+    return d;
+}
+
 template <class Basis>
 void testEnergy(const Basis& basis, const std::vector<TargetSpace>& coefficients)
 {
@@ -102,6 +113,10 @@ void testUnitVector3d()
         for (int j=0; j<dim+1; j++)
             coefficients[j] = testPoints[index[j]];
 
+        // This may be overly restrictive, but see the TODO in targetspacetrsolver.cc
+        if (diameter(coefficients) > TargetSpace::convexityRadius)
+            continue;
+
         testEnergy<Basis>(basis, coefficients);
         
     }