diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 6e2fb19f32550a872a1c361c3ddba7fadfc063d3..d1274f81b49190a0133a5f2e6552e1091300a7b8 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -82,6 +82,8 @@ public:  // for testing
         // This second derivative is almost given by the method getFirstDerivativesOfDirectors.
         // However, since the directors of a given unit quaternion are the _columns_ of the
         // corresponding orthogonal matrix, we need to invert the i and j indices
+        //
+        // So, if I am not mistaken, DR[i][j][k] contains \partial M_ij / \partial k
         Tensor3<double,3 , 3, 4> dd_dq;
         value.q.getFirstDerivativesOfDirectors(dd_dq);
         
@@ -136,7 +138,8 @@ public:
 
     RT curvatureEnergy(const Tensor3<double,3,3,3>& DR) const
     {
-        return mu_ * std::pow(L_c_ * curl(DR).frobenius_norm(),q_);
+        //return mu_ * std::pow(L_c_ * curl(DR).frobenius_norm(),q_);
+        return mu_ * std::pow(L_c_ * DR.frobenius_norm(),q_);
     }
 
     RT bendingEnergy(const Dune::FieldMatrix<double,dim,dim>& R, const Tensor3<double,3,3,3>& DR) const
@@ -180,7 +183,7 @@ public:
     /** \brief Curvature exponent */
     double q_;
     
-    
+    const BoundaryPatch<GridView>* neumannBoundary_;
 };
 
 template <class GridView, class LocalFiniteElement, int dim>
@@ -195,7 +198,7 @@ energy(const Entity& element,
     LocalGeodesicFEFunction<gridDim, double, LocalFiniteElement, TargetSpace> localGeodesicFEFunction(localFiniteElement,
                                                                                                       localSolution);
 
-    int quadOrder = 1;//gridDim;
+    int quadOrder = 2*gridDim;
 
     const Dune::QuadratureRule<double, gridDim>& quad 
         = Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
@@ -276,7 +279,37 @@ energy(const Entity& element,
             DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
 
     }
+    
+    
+    //////////////////////////////////////////////////////////////////////////////
+    //   Assemble boundary contributions
+    //////////////////////////////////////////////////////////////////////////////
+
+    for (typename Entity::LeafIntersectionIterator it = element.ileafbegin(); it != element.ileafend(); ++it) {
+        
+        if (not neumannBoundary_ or not neumannBoundary_->contains(*it))
+            continue;
+        
+        const Dune::QuadratureRule<double, gridDim-1>& quad 
+            = Dune::QuadratureRules<double, gridDim-1>::rule(it->type(), quadOrder);
+    
+        for (size_t pt=0; pt<quad.size(); pt++) {
+        
+            // Local position of the quadrature point
+            const Dune::FieldVector<double,gridDim>& quadPos = it->geometryInInside().global(quad[pt].position());
+        
+            const double integrationElement = it->geometry().integrationElement(quad[pt].position());
 
+            // The value of the local function
+            RigidBodyMotion<dim> value = localGeodesicFEFunction.evaluate(quadPos);
+
+            // Constant Neumann force in z-direction
+            energy += thickness_ * 40 * value.r[2] * quad[pt].weight() * integrationElement;
+            
+        }
+        
+    }
+    
     return energy;
 }
 
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index e3f02bf43dd26360d29b0cfe1c2d8c44fde21a05..95be5a8693671c54acf130294b7b6e31105b08f9 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -62,7 +62,9 @@ setup(const GridType& grid,
     //   Create a multigrid solver
     // ////////////////////////////////
 
+#if 0
     // First create a Gauss-seidel base solver
+    
     TrustRegionGSStep<MatrixType, CorrectionType>* baseSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
 
     EnergyNorm<MatrixType, CorrectionType>* baseEnergyNorm = new EnergyNorm<MatrixType, CorrectionType>(*baseSolverStep);
@@ -72,6 +74,11 @@ setup(const GridType& grid,
                                                                             baseTolerance,
                                                                             baseEnergyNorm,
                                                                             Solver::QUIET);
+#else
+    QuadraticIPOptSolver<MatrixType,CorrectionType>* baseSolver = new QuadraticIPOptSolver<MatrixType,CorrectionType>();
+    baseSolver->tolerance_ = baseTolerance;
+    baseSolver->verbosity_ = Solver::QUIET;
+#endif
 
     // Make pre and postsmoothers
     TrustRegionGSStep<MatrixType, CorrectionType>* presmoother  = new TrustRegionGSStep<MatrixType, CorrectionType>;
@@ -105,12 +112,21 @@ setup(const GridType& grid,
 
     h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
 
+#if 0
     innerSolver_ = std::shared_ptr<LoopSolver<CorrectionType> >(new ::LoopSolver<CorrectionType>(mmgStep,
                                                                                                    innerIterations_,
                                                                                                    innerTolerance_,
                                                                                                    h1SemiNorm_,
                                                                                                  Solver::FULL));
 
+#else
+    QuadraticIPOptSolver<MatrixType,CorrectionType>* solver = new QuadraticIPOptSolver<MatrixType,CorrectionType>();
+    solver->tolerance_ = innerTolerance_;
+    solver->verbosity_ = Solver::REDUCED;
+    innerSolver_ = std::shared_ptr<QuadraticIPOptSolver<MatrixType,CorrectionType> >(solver);
+#endif
+
+    
     // Write all intermediate solutions, if requested
     if (instrumented_
         && dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
@@ -190,7 +206,7 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
 
     std::vector<std::vector<BoxConstraint<field_type,blocksize> > > trustRegionObstacles((mgStep) 
                                                                                          ? mgStep->numLevels_
-                                                                                         : 0);
+                                                                                         : 1);
 
    // /////////////////////////////////////////////////////
     //   Set up the log file, if requested
@@ -242,8 +258,6 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
         // The right hand side is the _negative_ gradient
         rhs *= -1;
 
-/*        std::cout << "rhs:\n" << rhs << std::endl;
-        std::cout << "matrix[0][0]:\n" << (*hessianMatrix_)[0][0] << std::endl;*/
         
         // //////////////////////////////////////////////////////////////////////
         //   Modify matrix and right-hand side to account for Dirichlet values
@@ -277,10 +291,17 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
 
         }
 
-        mgStep->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1);
+        //std::cout << "rhs:\n" << rhs << std::endl;
+        //std::cout << "matrix[0][0]:\n" << (*hessianMatrix_)[0][0] << std::endl;
+        //printmatrix(std::cout, (*hessianMatrix_), "hessian", "--");
+        //writeMatrixToMatlab(*hessianMatrix_, "cosserat_hessian");
+        
+        //mgStep->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1);
+        std::dynamic_pointer_cast<QuadraticIPOptSolver<MatrixType,CorrectionType> >(innerSolver_)->setProblem(*hessianMatrix_, corr, rhs);
         
         trustRegionObstacles.back() = trustRegion.obstacles();
-        mgStep->obstacles_ = &trustRegionObstacles;
+        //mgStep->obstacles_ = &trustRegionObstacles;
+        std::dynamic_pointer_cast<QuadraticIPOptSolver<MatrixType,CorrectionType> >(innerSolver_)->obstacles_ = &trustRegionObstacles.back();
 
         innerSolver_->preprocess();
         
@@ -420,7 +441,7 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
         }
 
         if (energy >= oldEnergy &&
-            (std::abs(oldEnergy-energy)/energy < 1e-9 || modelDecrease/energy < 1e-9)) {
+            (std::abs((oldEnergy-energy)/energy) < 1e-9 || std::abs(modelDecrease/energy) < 1e-9)) {
             if (this->verbosity_ == NumProc::FULL)
                 std::cout << "Suspecting rounding problems" << std::endl;
 
diff --git a/dune/gfe/rodassembler.cc b/dune/gfe/rodassembler.cc
index 6ab87a730b25587e38821837790817e25b5faf3f..9fa2e96e847942862b0cbd937969cc69940cbcd6 100644
--- a/dune/gfe/rodassembler.cc
+++ b/dune/gfe/rodassembler.cc
@@ -18,9 +18,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
 {
     using namespace Dune;
 
-    const typename GridView::Traits::IndexSet& indexSet = this->basis_.getGridView().indexSet();
-
-    if (sol.size()!=indexSet.size(gridDim))
+    if (sol.size()!=this->basis_.size())
         DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
 
     grad.resize(sol.size());
@@ -39,7 +37,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
         std::vector<RigidBodyMotion<3> > localSolution(nDofs);
         
         for (int i=0; i<nDofs; i++)
-            localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+            localSolution[i] = sol[this->basis_.index(*it,i)];
 
         // Assemble local gradient
         std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
@@ -51,7 +49,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
 
         // Add to global gradient
         for (int i=0; i<nDofs; i++)
-            grad[indexSet.subIndex(*it,i,gridDim)] += localGradient[i];
+            grad[this->basis_.index(*it,i)] += localGradient[i];
 
     }
 
@@ -65,17 +63,17 @@ getStrain(const std::vector<RigidBodyMotion<3> >& sol,
 {
     using namespace Dune;
 
-    const typename GridView::Traits::IndexSet& indexSet = this->gridView_.indexSet();
+    const typename GridView::Traits::IndexSet& indexSet = this->basis_.getGridView().indexSet();
 
-    if (sol.size()!=indexSet.size(gridDim))
+    if (sol.size()!=this->basis_.size())
         DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
 
     // Strain defined on each element
     strain.resize(indexSet.size(0));
     strain = 0;
 
-    ElementIterator it    = this->gridView_.template begin<0>();
-    ElementIterator endIt = this->gridView_.template end<0>();
+    ElementIterator it    = this->basis_.getGridView().template begin<0>();
+    ElementIterator endIt = this->basis_.getGridView().template end<0>();
 
     // Loop over all elements
     for (; it!=endIt; ++it) {
@@ -227,15 +225,13 @@ void RodAssembler<GridView,2>::
 assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
                Dune::BCRSMatrix<MatrixBlock>& matrix)
 {
-    const typename GridView::IndexSet& indexSet = this->gridView_.indexSet();
-
     Dune::MatrixIndexSet neighborsPerVertex;
     this->getNeighborsPerVertex(neighborsPerVertex);
     
     matrix = 0;
     
-    ElementIterator it    = this->gridView_.template begin<0>();
-    ElementIterator endit = this->gridView_.template end<0>  ();
+    ElementIterator it    = this->basis_.getGridView().template begin<0>();
+    ElementIterator endit = this->basis_.getGridView().template end<0>  ();
 
     Dune::Matrix<MatrixBlock> mat;
     
@@ -247,7 +243,7 @@ assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
         std::vector<RigidBodyMotion<2> > localSolution(numOfBaseFct);
         
         for (int i=0; i<numOfBaseFct; i++)
-            localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+            localSolution[i] = sol[this->basis_.index(*it,i)];
 
         // setup matrix 
         getLocalMatrix( *it, localSolution, mat);
@@ -255,11 +251,11 @@ assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
         // Add element matrix to global stiffness matrix
         for(int i=0; i<numOfBaseFct; i++) { 
             
-            int row = indexSet.subIndex(*it,i,gridDim);
+            int row = this->basis_.index(*it,i);
 
             for (int j=0; j<numOfBaseFct; j++ ) {
                 
-                int col = indexSet.subIndex(*it,j,gridDim);
+                int col = this->basis_.index(*it,j);
                 matrix[row][col] += mat[i][j];
                 
             }
@@ -280,8 +276,6 @@ getLocalMatrix( EntityType &entity,
                 const std::vector<RigidBodyMotion<2> >& localSolution,
                 Dune::Matrix<MatrixBlock>& localMat) const
 {
-    const typename GridView::IndexSet& indexSet = this->gridView_.indexSet();
-
     /* ndof is the number of vectors of the element */
     int ndof = localSolution.size();
 
@@ -429,16 +423,14 @@ void RodAssembler<GridView,2>::
 assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
                  Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
 {
-    const typename GridView::IndexSet& indexSet = this->gridView_.indexSet();
-
-    if (sol.size()!=this->gridView_.size(gridDim))
+    if (sol.size()!=this->basis_.size())
         DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
 
     grad.resize(sol.size());
     grad = 0;
 
-    ElementIterator it    = this->gridView_.template begin<0>();
-    ElementIterator endIt = this->gridView_.template end<0>();
+    ElementIterator it    = this->basis_.getGridView().template begin<0>();
+    ElementIterator endIt = this->basis_.getGridView().template end<0>();
 
     // Loop over all elements
     for (; it!=endIt; ++it) {
@@ -450,7 +442,7 @@ assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
         RigidBodyMotion<2> localSolution[numOfBaseFct];
         
         for (int i=0; i<numOfBaseFct; i++)
-            localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+            localSolution[i] = sol[this->basis_.index(*it,i)];
 
         // Get quadrature rule
         const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), 2);
@@ -500,7 +492,7 @@ assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
 
             for (int dof=0; dof<numOfBaseFct; dof++) {
 
-                int globalDof = indexSet.subIndex(*it,dof,gridDim);
+                int globalDof = this->basis_.index(*it,dof);
 
                 //printf("globalDof: %d   partA1: %g   partA3: %g\n", globalDof, partA1, partA3);
 
@@ -531,13 +523,11 @@ computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const
 {
     double energy = 0;
 
-    const typename GridView::IndexSet& indexSet = this->gridView_.indexSet();
-
-    if (sol.size()!=this->gridView_.size(gridDim))
+    if (sol.size()!=this->basis_.size())
         DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
 
-    ElementIterator it    = this->gridView_.template begin<0>();
-    ElementIterator endIt = this->gridView_.template end<0>();
+    ElementIterator it    = this->basis_.getGridView().template begin<0>();
+    ElementIterator endIt = this->basis_.getGridView().template end<0>();
 
     // Loop over all elements
     for (; it!=endIt; ++it) {
@@ -550,7 +540,7 @@ computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const
         RigidBodyMotion<2> localSolution[numOfBaseFct];
         
         for (int i=0; i<numOfBaseFct; i++)
-            localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+            localSolution[i] = sol[this->basis_.index(*it,i)];
 
         // Get quadrature rule
         const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), 2);