From fad74c2ccba84da622341ca6227a73e0c551bd4f Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 17 Apr 2009 13:32:23 +0000
Subject: [PATCH] a default implementation of the energy functional gradient
using finite differences
[[Imported from SVN: r4021]]
---
src/localgeodesicfestiffness.hh | 34 ++++++++++++++++++++++++++++-----
src/rodlocalstiffness.hh | 13 -------------
2 files changed, 29 insertions(+), 18 deletions(-)
diff --git a/src/localgeodesicfestiffness.hh b/src/localgeodesicfestiffness.hh
index ec31db80..1c0ed075 100644
--- a/src/localgeodesicfestiffness.hh
+++ b/src/localgeodesicfestiffness.hh
@@ -87,13 +87,37 @@ public:
};
-template <class GridType, class TargetSpace>
-void LocalGeodesicFEStiffness<GridType, TargetSpace>::
+template <class GridView, class TargetSpace>
+void LocalGeodesicFEStiffness<GridView, TargetSpace>::
assembleGradient(const Entity& element,
- const std::vector<TargetSpace>& solution,
- Dune::array<Dune::FieldVector<double,6>, 2>& gradient) const
+ const std::vector<TargetSpace>& localSolution,
+ Dune::array<Dune::FieldVector<double,6>, 2>& localGradient) const
{
-
+ // ///////////////////////////////////////////////////////////
+ // 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<6; j++) {
+
+ infinitesimalVariation(forwardSolution[i], eps, j);
+ infinitesimalVariation(backwardSolution[i], -eps, j);
+
+ localGradient[i][j] = (energy(element,forwardSolution) - energy(element,backwardSolution))
+ / (2*eps);
+
+ forwardSolution[i] = localSolution[i];
+ backwardSolution[i] = localSolution[i];
+ }
+
+ }
+
}
diff --git a/src/rodlocalstiffness.hh b/src/rodlocalstiffness.hh
index 272eaa28..8bab6c94 100644
--- a/src/rodlocalstiffness.hh
+++ b/src/rodlocalstiffness.hh
@@ -30,19 +30,6 @@ class RodLocalStiffness
enum {bendingQuadOrder = 2};
public:
- /** \brief For the fd approximations
- \todo This is public because RodAssembler uses it
- */
- static void infinitesimalVariation(RigidBodyMotion<3>& c, double eps, int i)
- {
- if (i<3)
- c.r[i] += eps;
- else
- c.q = c.q.mult(Rotation<3,double>::exp((i==3)*eps,
- (i==4)*eps,
- (i==5)*eps));
- }
-
std::vector<RigidBodyMotion<3> > localReferenceConfiguration_;
public:
--
GitLab