From 032e865077325aae887894c0fa19cebde378e15f Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 13 Jun 2011 16:13:23 +0000
Subject: [PATCH] new implementation of the fd gradient. Hopefully less buggy
[[Imported from SVN: r7424]]
---
dune/gfe/localgeodesicfestiffness.hh | 43 +++++++++++++++++++++++-----
1 file changed, 36 insertions(+), 7 deletions(-)
diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index dcfa715e..3c16f578 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -103,6 +103,8 @@ class LocalGeodesicFEStiffnessImp
// to a neighborhood around S^n
forwardSolution[i] = localSolution[i];
backwardSolution[i] = localSolution[i];
+
+
LocalGeodesicFEStiffnessImp<GridView,TargetSpace>::infinitesimalVariation(forwardSolution[i], eps, j);
LocalGeodesicFEStiffnessImp<GridView,TargetSpace>::infinitesimalVariation(backwardSolution[i], -eps, j);
@@ -184,16 +186,43 @@ assembleGradient(const Entity& element,
const std::vector<TargetSpace>& localSolution,
std::vector<typename TargetSpace::TangentVector>& localGradient) const
{
- std::vector<typename TargetSpace::EmbeddedTangentVector> embeddedLocalGradient;
- // first compute the gradient in embedded coordinates
- LocalGeodesicFEStiffnessImp<GridView,TargetSpace>::assembleEmbeddedGradient(element, localSolution, embeddedLocalGradient, this);
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
- // transform to coordinates on the tangent space
- localGradient.resize(embeddedLocalGradient.size());
+ double eps = 1e-6;
+
+ localGradient.resize(localSolution.size());
+
+ std::vector<TargetSpace> forwardSolution = localSolution;
+ std::vector<TargetSpace> backwardSolution = localSolution;
+
+ for (size_t i=0; i<localSolution.size(); i++) {
+
+ // basis vectors of the tangent space of the i-th entry of localSolution
+ const Dune::FieldMatrix<double,blocksize,embeddedBlocksize> B = localSolution[i].orthonormalFrame();
+
+ for (int j=0; j<blocksize; j++) {
+
+ typename TargetSpace::EmbeddedTangentVector forwardCorrection = B[j];
+ forwardCorrection *= eps;
+
+ typename TargetSpace::EmbeddedTangentVector backwardCorrection = B[j];
+ backwardCorrection *= -eps;
+
+ forwardSolution[i] = TargetSpace::exp(localSolution[i], forwardCorrection);
+ backwardSolution[i] = TargetSpace::exp(localSolution[i], backwardCorrection);
+
+ localGradient[i][j] = (energy(element,forwardSolution) - energy(element,backwardSolution)) / (2*eps);
+
+ }
+
+ forwardSolution[i] = localSolution[i];
+ backwardSolution[i] = localSolution[i];
+
+ }
- for (size_t i=0; i<localGradient.size(); i++)
- localSolution[i].orthonormalFrame().mv(embeddedLocalGradient[i], localGradient[i]);
}
template <class GridView, class TargetSpace>
--
GitLab