diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index 2ba48011010a4cea99c62dc9e478c7972e62202e..92e45badc309eeecd4e8d128103f6fbb6913dede 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -62,7 +62,6 @@ setup(const typename BasisType::GridView::Grid& grid,
     baseTolerance_            = baseTolerance;
     damping_                  = damping;
 
-    int numLevels = grid_->maxLevel()+1;
     const auto dim = VectorType::value_type::dimension;
 
 #if HAVE_DUNE_PARMG
@@ -289,6 +288,8 @@ setup(const typename BasisType::GridView::Grid& grid,
       isP1Basis = std::is_same<Basis,Dune::Functions::LagrangeBasis<typename Basis::GridView, 1> >::value;
     }
 
+    int numLevels = grid_->maxLevel()+1;
+
     using TransferOperatorType = typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType;
     std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<CorrectionType>>> transferOperators(isP1Basis ? numLevels-1 : numLevels);
 
@@ -550,16 +551,12 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
 
             if (correctionInfinityNorm < this->tolerance_) {
                 if (this->verbosity_ == NumProc::FULL and rank==0)
-                    std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
-
-                if (this->verbosity_ != NumProc::QUIET and rank==0)
-                    std::cout << i+1 << " trust-region steps were taken." << std::endl;
-                break;
-            }
-
-            if (trustRegion.radius() < this->tolerance_) {
-              if (this->verbosity_ == NumProc::FULL and rank==0)
-                    std::cout << "TRUST REGION RADIUS IS TOO SMALL, SMALLER THAN TOLERANCE, STOPPING NOW" << std::endl;
+                {
+                    if (correctionInfinityNorm < trustRegion.radius())
+                        std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+                    else
+                        std::cout << "TRUST-REGION UNDERFLOW!" << std::endl;
+                }
 
                 if (this->verbosity_ != NumProc::QUIET and rank==0)
                     std::cout << i+1 << " trust-region steps were taken." << std::endl;
diff --git a/dune/elasticity/materials/localintegralenergy.hh b/dune/elasticity/materials/localintegralenergy.hh
index e19d9c831b84f4a580253eb52d58aa6dc4bfd0a1..6ef381539d8163f6409229c65dcd72631591587d 100644
--- a/dune/elasticity/materials/localintegralenergy.hh
+++ b/dune/elasticity/materials/localintegralenergy.hh
@@ -2,6 +2,7 @@
 #define DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
 
 #include <dune/common/fmatrix.hh>
+#include <dune/common/transpose.hh>
 
 #include <dune/geometry/quadraturerules.hh>
 
@@ -53,8 +54,9 @@ energy(const LocalView& localView,
   field_type energy = 0;
 
   // store gradients of shape functions and base functions
-  std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
-  std::vector<FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+  using FiniteElement = typename LocalView::Tree::ChildType::FiniteElement;
+  using Jacobian = typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
+  std::vector<Jacobian> jacobians(localFiniteElement.size());
 
   int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
                                                : localFiniteElement.localBasis().order() * gridDim;
@@ -65,23 +67,23 @@ energy(const LocalView& localView,
   {
     const DT integrationElement = element.geometry().integrationElement(qp.position());
 
-    const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
+    const auto geometryJacobianIT = element.geometry().jacobianInverseTransposed(qp.position());
 
     // Global position
     auto x = element.geometry().global(qp.position());
 
     // Get gradients of shape functions
-    localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
+    localFiniteElement.localBasis().evaluateJacobian(qp.position(), jacobians);
 
     // compute gradients of base functions
-    for (size_t i=0; i<gradients.size(); ++i)
-      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+    for (size_t i=0; i<jacobians.size(); ++i)
+      jacobians[i] = jacobians[i] * transpose(geometryJacobianIT);
 
     // Deformation gradient
     FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
-    for (size_t i=0; i<gradients.size(); i++)
+    for (size_t i=0; i<jacobians.size(); i++)
       for (int j=0; j<gridDim; j++)
-        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], gradients[i]);
+        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], jacobians[i][0]);
     // Integrate the energy density
     energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
   }