diff --git a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
index 45a56cc1fa991964de28abdb61a843a417ea2560..49937c2c5c07afa854581b5ebd08b47bc4e6c311 100644
--- a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
@@ -112,14 +112,7 @@ energy(const typename Basis::LocalView& localView,
energy >>= pureEnergy;
trace_off();
-#if 0
- size_t tape_stats[STAT_SIZE];
- tapestats(rank,tape_stats); // reading of tape statistics
- cout<<"maxlive "<<tape_stats[NUM_MAX_LIVES]<<"\n";
- cout<<"tay_stack_size "<<tape_stats[TAY_STACK_SIZE]<<"\n";
- cout<<"total number of operations "<<tape_stats[NUM_OPERATIONS]<<"\n";
- // ..... print other tape stats
-#endif
+
return pureEnergy;
}
@@ -156,10 +149,6 @@ assembleGradient(const typename Basis::LocalView& localView,
for (size_t j=0; j<embeddedBlocksize; j++)
localEmbeddedGradient[i][j] = g[idx++];
- // std::cout << "localEmbeddedGradient:\n";
- // for (size_t i=0; i<nDofs; i++)
- // std::cout << localEmbeddedGradient[i] << std::endl;
-
// Express gradient in local coordinate system
for (size_t i=0; i<nDofs; i++) {
Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> orthonormalFrame = localSolution[i].orthonormalFrame();
@@ -358,9 +347,6 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
}
}
-
- // std::cout << "ADOL-C stiffness:\n";
- // printmatrix(std::cout, this->A_, "foo", "--");
}
#endif
diff --git a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
index 6ad5ad3e1832628c6cc09c379e1da99112b4d462..fd2a36c4b21c489c4cc61921309f27015f73f95d 100644
--- a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
@@ -52,14 +52,6 @@ public:
/** \brief Assemble the local tangent matrix and gradient at the current position
This implementation uses finite-difference approximations
-
- The formula for the Riemannian Hessian has been taken from Absil, Mahony, Sepulchre:
- 'Optimization algorithms on matrix manifolds', page 107. There it says that
- \f[
- \langle Hess f(x)[\xi], \eta \rangle
- = \frac 12 \frac{d^2}{dt^2} \Big(f(\exp_x(t(\xi + \eta))) - f(\exp_x(t\xi)) - f(\exp_x(t\eta))\Big)\Big|_{t=0}.
- \f]
- We compute that using a finite difference approximation.
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
diff --git a/dune/gfe/assemblers/localgeodesicfestiffness.hh b/dune/gfe/assemblers/localgeodesicfestiffness.hh
index 9dd98819b8f430c7647aaef85241b24c25a01672..4b3ccafef695070fc1d36f97fb1ee9457bfa6a91 100644
--- a/dune/gfe/assemblers/localgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfestiffness.hh
@@ -39,8 +39,7 @@ public:
const std::vector<TargetSpace>& localSolution) const = 0;
/** \brief Assemble the element gradient of the energy functional
-
- The default implementation in this class uses a finite difference approximation */
+ */
virtual void assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& solution,
std::vector<typename TargetSpace::TangentVector>& gradient) const = 0;
diff --git a/dune/gfe/assemblers/mixedgfeassembler.hh b/dune/gfe/assemblers/mixedgfeassembler.hh
index a57b86d754a4351a519f733e3c8b23738c91a9e2..29c61a41aca5f0e690cd870788bb65f4ece0d114 100644
--- a/dune/gfe/assemblers/mixedgfeassembler.hh
+++ b/dune/gfe/assemblers/mixedgfeassembler.hh
@@ -60,11 +60,7 @@ public:
Dune::BlockVector<Dune::FieldVector<double, blocksize1> >& gradient1,
MatrixType& hessian,
bool computeOccupationPattern=true) const;
-#if 0
- /** \brief Assemble the gradient */
- virtual void assembleGradient(const std::vector<TargetSpace>& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
-#endif
+
/** \brief Compute the energy of a deformation state */
virtual double computeEnergy(const std::vector<TargetSpace0>& configuration0,
const std::vector<TargetSpace1>& configuration1) const;
@@ -255,47 +251,6 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
}
}
-#if 0
-template <class Basis, class TargetSpace>
-void GeodesicFEAssembler<Basis,TargetSpace>::
-assembleGradient(const std::vector<TargetSpace>& sol,
- Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
-{
- if (sol.size()!=basis_.size())
- DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
-
- grad.resize(sol.size());
- grad = 0;
-
- ElementIterator it = basis_.getGridView().template begin<0,Dune::Interior_Partition>();
- ElementIterator endIt = basis_.getGridView().template end<0,Dune::Interior_Partition>();
-
- // Loop over all elements
- for (; it!=endIt; ++it) {
-
- // A 1d grid has two vertices
- const int nDofs = basis_.getLocalFiniteElement(*it).localBasis().size();
-
- // Extract local solution
- std::vector<TargetSpace> localSolution(nDofs);
-
- for (int i=0; i<nDofs; i++)
- localSolution[i] = sol[basis_.index(*it,i)];
-
- // Assemble local gradient
- std::vector<Dune::FieldVector<double,blocksize> > localGradient(nDofs);
-
- localStiffness_->assembleGradient(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient);
-
- // Add to global gradient
- for (int i=0; i<nDofs; i++)
- grad[basis_.index(*it,i)] += localGradient[i];
-
- }
-
-}
-#endif
-
template <class Basis, class TargetSpace0, class TargetSpace1>
double MixedGFEAssembler<Basis, TargetSpace0, TargetSpace1>::
computeEnergy(const std::vector<TargetSpace0>& configuration0,
diff --git a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh b/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
index 59c946a5d138496762955ab858465a4ea9c6b54b..83d5212b57058b2e60899892c4b5645384bcac50 100644
--- a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
@@ -24,17 +24,6 @@ public:
constexpr static int orientationBlocksize = OrientationTargetSpace::TangentVector::dimension;
/** \brief Assemble the local stiffness matrix at the current position
-
- This default implementation used finite-difference approximations to compute the second derivatives
-
- The formula for the Riemannian Hessian has been taken from Absil, Mahony, Sepulchre:
- 'Optimization algorithms on matrix manifolds', page 107. There it says that
- \f[
- \langle Hess f(x)[\xi], \eta \rangle
- = \frac 12 \frac{d^2}{dt^2} \Big(f(\exp_x(t(\xi + \eta))) - f(\exp_x(t\xi)) - f(\exp_x(t\eta))\Big)\Big|_{t=0}.
- \f]
- We compute that using a finite difference approximation.
-
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<DeformationTargetSpace>& localDisplacementConfiguration,
@@ -49,18 +38,7 @@ public:
virtual RT energy (const typename Basis::LocalView& localView,
const std::vector<DeformationTargetSpace>& localDeformationConfiguration,
const std::vector<OrientationTargetSpace>& localOrientationConfiguration) const = 0;
-#if 0
- /** \brief Assemble the element gradient of the energy functional
- The default implementation in this class uses a finite difference approximation */
- virtual void assembleGradient(const Entity& element,
- const LocalFiniteElement& localFiniteElement,
- const std::vector<TargetSpace>& solution,
- std::vector<typename TargetSpace::TangentVector>& gradient) const;
- {
- DUNE_THROW(Dune::NotImplemented, "!");
- }
-#endif
// assembled data
Dune::Matrix<Dune::FieldMatrix<RT, deformationBlocksize, deformationBlocksize> > A00_;
Dune::Matrix<Dune::FieldMatrix<RT, deformationBlocksize, orientationBlocksize> > A01_;
diff --git a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
index b8399ea451fb095f7b69562c0b4dde3053143909..18a120a5f5806ebb3d0e2686c8247609116dbe66 100644
--- a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
@@ -53,16 +53,7 @@ public:
virtual RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace0>& localConfiguration0,
const std::vector<TargetSpace1>& localConfiguration1) const;
-#if 0
- /** \brief Assemble the element gradient of the energy functional
- This uses the automatic differentiation toolbox ADOL_C.
- */
- virtual void assembleGradient(const Entity& element,
- const LocalFiniteElement& localFiniteElement,
- const std::vector<TargetSpace>& solution,
- std::vector<typename TargetSpace::TangentVector>& gradient) const;
-#endif
/** \brief Assemble the local stiffness matrix at the current position
This uses the automatic differentiation toolbox ADOL_C.
@@ -133,60 +124,9 @@ energy(const typename Basis::LocalView& localView,
energy >>= pureEnergy;
trace_off();
-#if 0
- size_t tape_stats[STAT_SIZE];
- tapestats(1,tape_stats); // reading of tape statistics
- cout<<"maxlive "<<tape_stats[NUM_MAX_LIVES]<<"\n";
- cout<<"tay_stack_size "<<tape_stats[TAY_STACK_SIZE]<<"\n";
- cout<<"total number of operations "<<tape_stats[NUM_OPERATIONS]<<"\n";
- // ..... print other tape stats
-#endif
- return pureEnergy;
-}
-
-#if 0
-template <class GridView, class LocalFiniteElement, class TargetSpace>
-void LocalGeodesicFEADOLCStiffness<GridView, LocalFiniteElement, TargetSpace>::
-assembleGradient(const Entity& element,
- const LocalFiniteElement& localFiniteElement,
- const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient) const
-{
- // Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
- energy(element, localFiniteElement, localSolution);
-
- // Compute the actual gradient
- size_t nDofs = localSolution.size();
- size_t nDoubles = nDofs*embeddedBlocksize;
- std::vector<double> xp(nDoubles);
- int idx=0;
- for (size_t i=0; i<nDofs; i++)
- for (size_t j=0; j<embeddedBlocksize; j++)
- xp[idx++] = localSolution[i].globalCoordinates()[j];
-
- // Compute gradient
- std::vector<double> g(nDoubles);
- gradient(1,nDoubles,xp.data(),g.data()); // gradient evaluation
- // Copy into Dune type
- std::vector<typename TargetSpace::EmbeddedTangentVector> localEmbeddedGradient(localSolution.size());
-
- idx=0;
- for (size_t i=0; i<nDofs; i++)
- for (size_t j=0; j<embeddedBlocksize; j++)
- localEmbeddedGradient[i][j] = g[idx++];
-
- // std::cout << "localEmbeddedGradient:\n";
- // for (size_t i=0; i<nDofs; i++)
- // std::cout << localEmbeddedGradient[i] << std::endl;
-
- // 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]);
- }
+ return pureEnergy;
}
-#endif
// ///////////////////////////////////////////////////////////
// Compute gradient and Hessian together