Add try{}catch{} around the energy calculation and the solve() call inside the trustregion solver

In case an excpetion is thrown in either of these two calls, the trustregion solver will not accept the step and continue with a reduced radius

dune/elasticity/common/trustregionsolver.cc

@@ -442,6 +442,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
 
         CorrectionType corr(rhs.size());
         corr = 0;
+        bool solvedByInnerSolver = true;
 
         //Take the obstacles on the finest grid and give them to the multigrid solver, it will create a hierarchy for all coarser grids
         trustRegionObstacles = trustRegion.obstacles();
@@ -521,7 +522,20 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
         std::cout << "Solve quadratic problem..." << std::endl;
 
         Dune::Timer solutionTimer;
-        innerSolver_->solve();
+
+        try {
+            innerSolver_->solve();
+        } catch (Dune::Exception &e) {
+            std::cerr << "Error while solving: " << e << std::endl;
+            solvedByInnerSolver = false;
+            corr = 0;
+        }
+        double energy = 0;
+        double totalModelDecrease = 0;
+        SolutionType newIterate = x_;
+
+        if (solvedByInnerSolver) {
+
         std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
 
 #if ! HAVE_DUNE_PARMG
@@ -562,11 +576,20 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
         //   Check whether trust-region step can be accepted
         // ////////////////////////////////////////////////////
 
-        SolutionType newIterate = x_;
         for (size_t j=0; j<newIterate.size(); j++)
             newIterate[j] += corr[j];
 
-        double energy    = assembler_->computeEnergy(newIterate);
+        try {
+            energy  = assembler_->computeEnergy(newIterate);
+        } catch (Dune::Exception &e) {
+            std::cerr << "Error while computing the energy of the new Iterate: " << e << std::endl;
+            std::cerr << "Redoing trust region step with smaller radius..." << std::endl;
+            newIterate = x_;
+            solvedByInnerSolver = false;
+            energy = oldEnergy;
+        }
+        if (solvedByInnerSolver) {
+
         energy = grid_->comm().sum(energy);
 
         // compute the model decrease
@@ -576,7 +599,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
         tmp = 0;
         hessianMatrix_->umv(corr, tmp);
         double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
-        double totalModelDecrease = grid_->comm().sum(modelDecrease);
+        totalModelDecrease = grid_->comm().sum(modelDecrease);
 
         assert(totalModelDecrease >= 0);
 
@@ -607,11 +630,13 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
             x_ = newIterate;
             break;
         }
+        }
+        }
 
         // //////////////////////////////////////////////
         //   Check for acceptance of the step
         // //////////////////////////////////////////////
-        if (energy < oldEnergy && (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.9) {
+        if (solvedByInnerSolver && energy < oldEnergy && (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.9) {
             // very successful iteration
 
             x_ = newIterate;
@@ -622,7 +647,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
 
             recomputeGradientHessian = true;
 
-        } else if ( ( (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.01 || std::fabs(oldEnergy-energy) < 1e-12)) {
+        } else if (solvedByInnerSolver && ((oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.01 || std::fabs(oldEnergy-energy) < 1e-12)) {
             // successful iteration
             x_ = newIterate;