From b4d5ae4addddc3e197ed0116369f945d14fcac93 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 21 Apr 2009 11:38:27 +0000
Subject: [PATCH] fd approximation of the gradient in a dimension-independent
way
[[Imported from SVN: r4072]]
---
src/localgeodesicfefunction.hh | 52 +++++++++++-----------------------
1 file changed, 16 insertions(+), 36 deletions(-)
diff --git a/src/localgeodesicfefunction.hh b/src/localgeodesicfefunction.hh
index 06d3f944..25e27848 100644
--- a/src/localgeodesicfefunction.hh
+++ b/src/localgeodesicfefunction.hh
@@ -259,46 +259,26 @@ template <int dim, class ctype, class TargetSpace>
Dune::FieldMatrix<ctype, TargetSpace::EmbeddedTangentVector::size, dim> LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local)
{
- Dune::FieldMatrix<ctype, EmbeddedTangentVector::size, dim> result;
-
- if (dim==1) {
-
- EmbeddedTangentVector tmp = TargetSpace::interpolateDerivative(coefficients_[0], coefficients_[1], local[0]);
-
- for (int i=0; i<EmbeddedTangentVector::size; i++)
- result[i][0] = tmp[i];
-
- }
-
- if (dim==2) {
-
- double eps = 1e-6;
+ double eps = 1e-6;
- Dune::FieldVector<ctype, dim> fX = local;
- Dune::FieldVector<ctype, dim> bX = local;
- Dune::FieldVector<ctype, dim> fY = local;
- Dune::FieldVector<ctype, dim> bY = local;
-
- fX[0] += eps;
- bX[0] -= eps;
- fY[1] += eps;
- bY[1] -= eps;
-
- Quaternion<double> dX = (evaluate(fX) - evaluate(bX));
- Quaternion<double> dY = (evaluate(fY) - evaluate(bY));
-
- dX /= 2*eps;
- dY /= 2*eps;
-
- for (int i=0; i<EmbeddedTangentVector::size; i++){
- result[i][0] = dX[i];
- result[i][1] = dY[i];
- }
+ Dune::FieldMatrix<ctype, EmbeddedTangentVector::size, dim> result;
+ for (int i=0; i<dim; i++) {
+
+ Dune::FieldVector<ctype, dim> forward = local;
+ Dune::FieldVector<ctype, dim> backward = local;
+
+ forward[i] += eps;
+ backward[i] -= eps;
+
+ EmbeddedTangentVector fdDer = evaluate(forward) - evaluate(backward);
+ fdDer /= 2*eps;
+
+ for (int j=0; j<EmbeddedTangentVector::size; j++)
+ result[j][i] = fdDer[j];
+
}
- assert(dim==1 || dim==2);
-
return result;
}
--
GitLab