diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index f7ed1ab0ecc2f9f53961b07892e7798d1e4a97ed..1440bde74f8ced098538d1dc66efa012161e952e 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/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,97 +395,116 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
std::cout << "Solve quadratic problem..." << std::endl;
Dune::Timer solutionTimer;
- innerSolver_->solve();
- std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+ try {
+ innerSolver_->solve();
+ } 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
- if (mgStep)
- corr = mgStep->getSol();
+ std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+ if (mgStep)
+ corr = mgStep->getSol();
- auto correctionInfinityNorm = corr.infinity_norm();
+ auto correctionInfinityNorm = corr.infinity_norm();
#else
- corr = iterationStep->getSol();
+ std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds and " << step << " steps." << std::endl;
+ corr = iterationStep->getSol();
- auto correctionInfinityNorm = parallelInfinityNorm(corr);
+ 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;
-
- 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 (this->verbosity_ != NumProc::QUIET and rank==0)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
- break;
- }
-
- // ////////////////////////////////////////////////////
- // 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);
- energy = grid_->comm().sum(energy);
-
- // compute the model decrease
- // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
- // Note that rhs = -g
- CorrectionType tmp(corr.size());
- tmp = 0;
- hessianMatrix_->umv(corr, tmp);
- double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
- double totalModelDecrease = grid_->comm().sum(modelDecrease);
-
- assert(totalModelDecrease >= 0);
-
- double relativeModelDecrease = totalModelDecrease / std::fabs(energy);
- double relativeFunctionalDecrease = (oldEnergy - energy)/std::fabs(energy);
-
- if (this->verbosity_ == NumProc::FULL and rank==0) {
- std::cout << "Absolute model decrease: " << totalModelDecrease << std::endl;
- std::cout << "Functional decrease: " << oldEnergy - energy << std::endl;
- std::cout << "oldEnergy = " << oldEnergy << " and new energy = " << energy << std::endl;
- std::cout << "Relative model decrease: " << relativeModelDecrease << std::endl;
- std::cout << "Relative functional decrease: " << relativeFunctionalDecrease << std::endl;
- }
-
-
- if (energy >= oldEnergy) {
- if (this->verbosity_ == NumProc::FULL and rank==0)
- printf("Richtung ist keine Abstiegsrichtung!\n");
- }
-
- if (energy >= oldEnergy && (std::fabs(relativeFunctionalDecrease) < 1e-9 || std::fabs(relativeModelDecrease) < 1e-9)) {
- if (this->verbosity_ == NumProc::FULL and rank==0)
- std::cout << "Suspecting rounding problems" << std::endl;
-
- if (this->verbosity_ != NumProc::QUIET and rank==0)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
-
- x_ = newIterate;
- break;
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Infinity norm of the correction: " << correctionInfinityNorm << std::endl;
+
+ 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 (this->verbosity_ != NumProc::QUIET and rank==0)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+ break;
+ }
+
+ // ////////////////////////////////////////////////////
+ // Check whether trust-region step can be accepted
+ // ////////////////////////////////////////////////////
+
+ for (size_t j=0; j<newIterate.size(); j++)
+ newIterate[j] += corr[j];
+
+ 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
+ // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
+ // Note that rhs = -g
+ CorrectionType tmp(corr.size());
+ tmp = 0;
+ hessianMatrix_->umv(corr, tmp);
+ double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
+ double totalModelDecrease = grid_->comm().sum(modelDecrease);
+
+ assert(totalModelDecrease >= 0);
+
+ double relativeModelDecrease = totalModelDecrease / std::fabs(energy);
+ double relativeFunctionalDecrease = (oldEnergy - energy)/std::fabs(energy);
+
+ if (this->verbosity_ == NumProc::FULL and rank==0) {
+ std::cout << "Absolute model decrease: " << totalModelDecrease << std::endl;
+ std::cout << "Functional decrease: " << oldEnergy - energy << std::endl;
+ std::cout << "oldEnergy = " << oldEnergy << " and new energy = " << energy << std::endl;
+ std::cout << "Relative model decrease: " << relativeModelDecrease << std::endl;
+ std::cout << "Relative functional decrease: " << relativeFunctionalDecrease << std::endl;
+ }
+
+
+ if (energy >= oldEnergy) {
+ if (this->verbosity_ == NumProc::FULL and rank==0)
+ printf("Richtung ist keine Abstiegsrichtung!\n");
+ }
+
+ if (energy >= oldEnergy && (std::fabs(relativeFunctionalDecrease) < 1e-9 || std::fabs(relativeModelDecrease) < 1e-9)) {
+ if (this->verbosity_ == NumProc::FULL and rank==0)
+ std::cout << "Suspecting rounding problems" << std::endl;
+
+ if (this->verbosity_ != NumProc::QUIET and rank==0)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+
+ 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;
diff --git a/dune/elasticity/materials/mooneyrivlinenergy.hh b/dune/elasticity/materials/mooneyrivlinenergy.hh
index 508112e0b8f59644d05fa60f8eb1fd61dd95f690..6f589619b3f19967297f7164038b94bc2c48537d 100644
--- a/dune/elasticity/materials/mooneyrivlinenergy.hh
+++ b/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)