diff --git a/dune/gfe/geodesicfeassembler.hh b/dune/gfe/geodesicfeassembler.hh
index 0be023a8a04559414b5385cc003f8cdde14e8245..b533168abc531ead8175cf8ea2606ac72f8951bc 100644
--- a/dune/gfe/geodesicfeassembler.hh
+++ b/dune/gfe/geodesicfeassembler.hh
@@ -43,11 +43,15 @@ public:
localStiffness_(localStiffness)
{}
- /** \brief Assemble the tangent stiffness matrix
+ /** \brief Assemble the tangent stiffness matrix and the functional gradient together
+ *
+ * This is more efficient than computing them separately, because you need the gradient
+ * anyway to compute the Riemannian Hessian.
*/
- virtual void assembleMatrix(const std::vector<TargetSpace>& sol,
- Dune::BCRSMatrix<MatrixBlock>& matrix,
- bool computeOccupationPattern=true) const;
+ virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+ Dune::BCRSMatrix<MatrixBlock>& hessian,
+ bool computeOccupationPattern=true) const;
/** \brief Assemble the gradient */
virtual void assembleGradient(const std::vector<TargetSpace>& sol,
@@ -97,19 +101,23 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
-assembleMatrix(const std::vector<TargetSpace>& sol,
- Dune::BCRSMatrix<MatrixBlock>& matrix,
- bool computeOccupationPattern) const
+assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
+ Dune::BCRSMatrix<MatrixBlock>& hessian,
+ bool computeOccupationPattern) const
{
if (computeOccupationPattern) {
Dune::MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex);
- neighborsPerVertex.exportIdx(matrix);
+ neighborsPerVertex.exportIdx(hessian);
}
- matrix = 0;
+ hessian = 0;
+
+ gradient.resize(sol.size());
+ gradient = 0;
ElementIterator it = basis_.getGridView().template begin<0>();
ElementIterator endit = basis_.getGridView().template end<0> ();
@@ -124,8 +132,10 @@ assembleMatrix(const std::vector<TargetSpace>& sol,
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[basis_.index(*it,i)];
- // setup matrix
- localStiffness_->assembleHessian(*it, basis_.getLocalFiniteElement(*it), localSolution);
+ std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct);
+
+ // setup local matrix and gradient
+ localStiffness_->assembleGradientAndHessian(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
@@ -135,11 +145,15 @@ assembleMatrix(const std::vector<TargetSpace>& sol,
for (int j=0; j<numOfBaseFct; j++ ) {
int col = basis_.index(*it,j);
- matrix[row][col] += localStiffness_->A_[i][j];
+ hessian[row][col] += localStiffness_->A_[i][j];
}
}
+ // Add local gradient to global gradient
+ for (int i=0; i<numOfBaseFct; i++)
+ gradient[basis_.index(*it,i)] += localGradient[i];
+
}
}
diff --git a/dune/gfe/localgeodesicfeadolcstiffness.hh b/dune/gfe/localgeodesicfeadolcstiffness.hh
index 9547060884b5400a6f59e0d020025d324d4da459..c1208508934f51203cfd71798beaca6cb481b0d5 100644
--- a/dune/gfe/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/localgeodesicfeadolcstiffness.hh
@@ -61,9 +61,10 @@ public:
This uses the automatic differentiation toolbox ADOL_C.
*/
- virtual void assembleHessian(const Entity& e,
+ virtual void assembleGradientAndHessian(const Entity& e,
const LocalFiniteElement& localFiniteElement,
- const std::vector<TargetSpace>& localSolution);
+ const std::vector<TargetSpace>& localSolution,
+ std::vector<typename TargetSpace::TangentVector>& localGradient);
const LocalGeodesicFEStiffness<GridView, LocalFiniteElement, ATargetSpace>* localEnergy_;
@@ -162,13 +163,16 @@ assembleGradient(const Entity& element,
// ///////////////////////////////////////////////////////////
-// Compute gradient by finite-difference approximation
+// Compute gradient and Hessian together
+// To compute the Hessian we need to compute the gradient anyway, so we may
+// as well return it. This saves assembly time.
// ///////////////////////////////////////////////////////////
template <class GridType, class LocalFiniteElement, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<GridType, LocalFiniteElement, TargetSpace>::
-assembleHessian(const Entity& element,
+assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement& localFiniteElement,
- const std::vector<TargetSpace>& localSolution)
+ const std::vector<TargetSpace>& localSolution,
+ std::vector<typename TargetSpace::TangentVector>& localGradient)
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(element, localFiniteElement, localSolution);
@@ -199,6 +203,12 @@ assembleHessian(const Entity& element,
for (size_t j=0; j<embeddedBlocksize; j++)
localEmbeddedGradient[i][j] = g[idx++];
+ // Express gradient in local coordinate system
+ for (size_t i=0; i<nDofs; i++) {
+ Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> orthonormalFrame = localSolution[i].orthonormalFrame();
+ orthonormalFrame.mv(localEmbeddedGradient[i],localGradient[i]);
+ }
+
/////////////////////////////////////////////////////////////////
// Compute Hessian
/////////////////////////////////////////////////////////////////
diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index 75afca5caa9c23b0d1b37df2ec0233b36f30bbd5..b223ce7710c237c24816501cec76fed7f8151bd4 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -39,9 +39,10 @@ public:
We compute that using a finite difference approximation.
*/
- virtual void assembleHessian(const Entity& e,
+ virtual void assembleGradientAndHessian(const Entity& e,
const LocalFiniteElement& localFiniteElement,
- const std::vector<TargetSpace>& localSolution);
+ const std::vector<TargetSpace>& localSolution,
+ std::vector<typename TargetSpace::TangentVector>& localGradient);
/** \brief Compute the energy at the current configuration */
virtual RT energy (const Entity& e,
@@ -114,9 +115,10 @@ assembleGradient(const Entity& element,
// ///////////////////////////////////////////////////////////
template <class GridType, class LocalFiniteElement, class TargetSpace>
void LocalGeodesicFEStiffness<GridType, LocalFiniteElement, TargetSpace>::
-assembleHessian(const Entity& element,
+assembleGradientAndHessian(const Entity& element,
const LocalFiniteElement& localFiniteElement,
- const std::vector<TargetSpace>& localSolution)
+ const std::vector<TargetSpace>& localSolution,
+ std::vector<typename TargetSpace::TangentVector>& localGradient)
{
// Number of degrees of freedom for this element
size_t nDofs = localSolution.size();
@@ -161,9 +163,22 @@ assembleHessian(const Entity& element,
}
- // finite-difference approximation
- // we loop over the lower left triangular half of the matrix.
- // The other half follows from symmetry
+ //////////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ //////////////////////////////////////////////////////////////
+
+ localGradient.resize(localSolution.size());
+
+ for (size_t i=0; i<localSolution.size(); i++)
+ for (int j=0; j<blocksize; j++)
+ localGradient[i][j] = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps);
+
+
+ ///////////////////////////////////////////////////////////////////////////
+ // Compute Riemannian Hesse matrix by finite-difference approximation.
+ // We loop over the lower left triangular half of the matrix.
+ // The other half follows from symmetry.
+ ///////////////////////////////////////////////////////////////////////////
//#pragma omp parallel for schedule (dynamic)
for (size_t i=0; i<localSolution.size(); i++) {
for (size_t i2=0; i2<blocksize; i2++) {
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index b3d329d87d37794299a8fea6ed335f6e86448e2d..9b977f39cccaa7e05b491e87861f748c2b21970c 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -240,19 +240,19 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
if (recomputeGradientHessian) {
- assembler_->assembleGradient(x_, rhs);
+ assembler_->assembleGradientAndHessian(x_,
+ rhs,
+ *hessianMatrix_,
+ i==0 // assemble occupation pattern only for the first call
+ );
+
rhs *= -1; // The right hand side is the _negative_ gradient
- if (this->verbosity_ == Solver::FULL)
- std::cout << "gradient assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
- gradientTimer.reset();
- assembler_->assembleMatrix(x_,
- *hessianMatrix_,
- i==0 // assemble occupation pattern only for the first call
- );
if (this->verbosity_ == Solver::FULL)
- std::cout << "hessian assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+ std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+
recomputeGradientHessian = false;
+
}
/* std::cout << "rhs:\n" << rhs << std::endl;