diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index d1274f81b49190a0133a5f2e6552e1091300a7b8..6e2fb19f32550a872a1c361c3ddba7fadfc063d3 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -82,8 +82,6 @@ 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);
@@ -138,8 +136,7 @@ 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_ * DR.frobenius_norm(),q_);
+ return mu_ * std::pow(L_c_ * curl(DR).frobenius_norm(),q_);
}
RT bendingEnergy(const Dune::FieldMatrix<double,dim,dim>& R, const Tensor3<double,3,3,3>& DR) const
@@ -183,7 +180,7 @@ public:
/** \brief Curvature exponent */
double q_;
- const BoundaryPatch<GridView>* neumannBoundary_;
+
};
template <class GridView, class LocalFiniteElement, int dim>
@@ -198,7 +195,7 @@ energy(const Entity& element,
LocalGeodesicFEFunction<gridDim, double, LocalFiniteElement, TargetSpace> localGeodesicFEFunction(localFiniteElement,
localSolution);
- int quadOrder = 2*gridDim;
+ int quadOrder = 1;//gridDim;
const Dune::QuadratureRule<double, gridDim>& quad
= Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
@@ -279,37 +276,7 @@ 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 95be5a8693671c54acf130294b7b6e31105b08f9..e3f02bf43dd26360d29b0cfe1c2d8c44fde21a05 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -62,9 +62,7 @@ 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);
@@ -74,11 +72,6 @@ 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>;
@@ -112,21 +105,12 @@ 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()))
@@ -206,7 +190,7 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
std::vector<std::vector<BoxConstraint<field_type,blocksize> > > trustRegionObstacles((mgStep)
? mgStep->numLevels_
- : 1);
+ : 0);
// /////////////////////////////////////////////////////
// Set up the log file, if requested
@@ -258,6 +242,8 @@ 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
@@ -291,17 +277,10 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
}
- //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);
+ mgStep->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1);
trustRegionObstacles.back() = trustRegion.obstacles();
- //mgStep->obstacles_ = &trustRegionObstacles;
- std::dynamic_pointer_cast<QuadraticIPOptSolver<MatrixType,CorrectionType> >(innerSolver_)->obstacles_ = &trustRegionObstacles.back();
+ mgStep->obstacles_ = &trustRegionObstacles;
innerSolver_->preprocess();
@@ -441,7 +420,7 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
}
if (energy >= oldEnergy &&
- (std::abs((oldEnergy-energy)/energy) < 1e-9 || std::abs(modelDecrease/energy) < 1e-9)) {
+ (std::abs(oldEnergy-energy)/energy < 1e-9 || modelDecrease/energy < 1e-9)) {
if (this->verbosity_ == NumProc::FULL)
std::cout << "Suspecting rounding problems" << std::endl;