From d18526907a179b4216c717bea1f37db067cf6c3d Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 21 Apr 2009 10:06:39 +0000
Subject: [PATCH] temporary: use the fd approximation of the gradient, and add
 some debugging code

[[Imported from SVN: r4069]]
---
 src/harmonicenergystiffness.hh | 23 +++++++++++++++++++++--
 1 file changed, 21 insertions(+), 2 deletions(-)

diff --git a/src/harmonicenergystiffness.hh b/src/harmonicenergystiffness.hh
index 42ea57fd..c0306a84 100644
--- a/src/harmonicenergystiffness.hh
+++ b/src/harmonicenergystiffness.hh
@@ -69,7 +69,7 @@ energy(const Entity& element,
         double weight = quad[pt].weight() * integrationElement;
 
         // The derivative of the local function defined on the reference element
-        Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos);
+        Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> referenceDerivative = localGeodesicFEFunction.evaluateDerivativeFD(quadPos);
 
         // The derivative of the function defined on the actual element
         Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::size, gridDim> derivative(0);
@@ -77,6 +77,25 @@ energy(const Entity& element,
         for (int comp=0; comp<4; comp++)
             jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
 
+        TargetSpace value = localGeodesicFEFunction.evaluate(quadPos);
+
+        for (int i=0; i<gridDim; i++) {
+            double dotproduct = 0;
+            for (int j=0; j<4; j++)
+                dotproduct += value[j]*derivative[j][i];
+            assert(std::fabs(dotproduct) < 1e-6);
+        }
+
+#if 0
+        double derivativenorm0 = 0;
+        double derivativenorm1 = 0;
+        for (int i=0; i<4; i++) {
+            derivativenorm0 += derivative[i][0]*derivative[i][0];
+            derivativenorm1 += derivative[i][1]*derivative[i][1];
+        }
+        std::cout << "Derivative norm: " << derivativenorm0 << ",  " << derivativenorm1 << std::endl;
+#endif
+
         for (int comp=0; comp<4; comp++) {
 
             for (int dir=0; dir<gridDim; dir++)
@@ -85,7 +104,7 @@ energy(const Entity& element,
         }
 
     }
-    
+
     return energy;
 }
 
-- 
GitLab