From 8e596f7f32ba71a1f6a747d51f4bb2799664c589 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sun, 28 Mar 2010 20:50:18 +0000
Subject: [PATCH] add a buggy but somewhat helpful implementation of the
Hessian
[[Imported from SVN: r5808]]
---
src/localgeodesicfestiffness.hh | 58 +++++++++++++++++++++++++++++++--
1 file changed, 55 insertions(+), 3 deletions(-)
diff --git a/src/localgeodesicfestiffness.hh b/src/localgeodesicfestiffness.hh
index 93eb6495..044bc52b 100644
--- a/src/localgeodesicfestiffness.hh
+++ b/src/localgeodesicfestiffness.hh
@@ -425,8 +425,6 @@ void LocalGeodesicFEStiffness<GridType,UnitVector<dim> >::
assemble(const Entity& element,
const std::vector<TargetSpace>& localSolution)
{
-#warning Dummy Hessian implementation
-
// 1 degree of freedom per element vertex
int nDofs = element.template count<gridDim>();
@@ -435,10 +433,63 @@ assemble(const Entity& element,
this->A = 0;
+#if 1
+#warning Dummy Hessian implementation
for (int i=0; i<nDofs; i++)
for (int j=0; j<blocksize; j++)
this->A[i][i][j][j] = 1;
-#if 0
+#else
+
+#if 1
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+
+ double eps = 1e-6;
+
+ std::vector<TargetSpace> forwardSolution = localSolution;
+ std::vector<TargetSpace> backwardSolution = localSolution;
+
+ for (size_t i=0; i<localSolution.size(); i++) {
+
+ for (int j=0; j<blocksize; j++) {
+
+ // The return value does not have unit norm. But assigning it to a UnitVector object
+ // will normalize it. This amounts to an extension of the energy functional
+ // to a neighborhood around S^n
+ forwardSolution[i] = infinitesimalVariation(localSolution[i], eps, j);
+ backwardSolution[i] = infinitesimalVariation(localSolution[i], -eps, j);
+
+ std::vector<Dune::FieldVector<double,blocksize> > forwardGradient;
+ std::vector<Dune::FieldVector<double,blocksize> > backwardGradient;
+ assembleGradient(element, forwardSolution, forwardGradient);
+ assembleGradient(element, backwardSolution, backwardGradient);
+
+ for (int k=0; k<localSolution.size(); k++)
+ for (int l=0; l<blocksize; l++)
+ this->A[i][k][j][l] = (forwardGradient[k][l] - backwardGradient[k][l]) / (2*eps);
+
+ forwardSolution[i] = localSolution[i];
+ backwardSolution[i] = localSolution[i];
+
+ }
+
+ }
+
+ // Project gradient in embedding space onto the tangent space
+ for (size_t i=0; i<localSolution.size(); i++)
+ for (size_t j=0; j<localSolution.size(); j++)
+ for (int k=0; k<blocksize; k++)
+ this->A[i][j][k] = localSolution[j].projectOntoTangentSpace(this->A[i][j][k]);
+
+
+
+
+
+#else
+
+
double eps = 1e-4;
typedef typename Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> >::row_type::iterator ColumnIterator;
@@ -578,6 +629,7 @@ assemble(const Entity& element,
}
#endif
+#endif
}
--
GitLab