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

Throw an exception in the Mooney-Rivlin-Energy-Calculation if det(∇φ) < 0
Then:     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

@@ -316,6 +316,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
 
         CorrectionType corr(rhs.size());
         corr = 0;
+        bool solved = 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();
@@ -394,22 +395,31 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
         std::cout << "Solve quadratic problem..." << std::endl;
 
         Dune::Timer solutionTimer;
+        try {
             innerSolver_->solve();
-        std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+        } catch (Dune::Exception &e) {
+            std::cerr << "Error while solving: " << e << std::endl;
+            solved = false;
+            corr = 0;
+        }
+        double energy = 0;
+        double totalModelDecrease = 0;
+        SolutionType newIterate = x_;
 
+        if (solved) {
 #if ! HAVE_DUNE_PARMG
+            std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
             if (mgStep)
                 corr = mgStep->getSol();
 
             auto correctionInfinityNorm = corr.infinity_norm();
 #else
+            std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds and " << step << " steps." << std::endl;
             corr = iterationStep->getSol();
 
             auto correctionInfinityNorm = parallelInfinityNorm(corr);
 #endif
 
-        //std::cout << "Correction: " << std::endl << corr << std::endl;
-
             if (this->verbosity_ == NumProc::FULL)
                 std::cout << "Infinity norm of the correction: " << correctionInfinityNorm << std::endl;
 
@@ -435,11 +445,19 @@ 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_;
+                solved = false;
+                energy = oldEnergy;
+            }
+            if (solved) {
                 energy = grid_->comm().sum(energy);
 
                 // compute the model decrease
@@ -480,11 +498,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 (solved && energy < oldEnergy && (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.9) {
             // very successful iteration
 
             x_ = newIterate;
@@ -495,7 +515,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 (solved && ((oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.01 || std::fabs(oldEnergy-energy) < 1e-12)) {
             // successful iteration
             x_ = newIterate;
 

dune/elasticity/materials/mooneyrivlinenergy.hh

@@ -152,6 +152,8 @@ energy(const Entity& element,
       for (int j = i+1; j < dim; j++)
         c2Tilde += sigmaSquared[j]*sigmaSquared[i];
     }
+    if (detF < 0)
+      DUNE_THROW( Dune::Exception, "det(∇φ) < 0, so it is not possible to calculate the MooneyRivlinEnergy.");
 
     field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3;
     // \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)

[oliver.sander_at_tu-dresden.de]

Can this variable name be a bit less generic? It is about the inner solver, right?

[lisa_julia.nebel_at_tu-dresden.de]

Changed it to solvedByInnerSolver, see commit 8fb2b257

[oliver.sander_at_tu-dresden.de]

Can you write "det(F)"? I am not sure whether C++ allows unicode strings nowadays.

[lisa_julia.nebel_at_tu-dresden.de]

Changed it to det(F), see commit 8fb2b257