From 8b343d922a40c6b4f531ba9b0efd1a5b9972822e Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 2 Sep 2010 08:36:09 +0000
Subject: [PATCH] test mixed second derivative

[[Imported from SVN: r6312]]
---
 test/unitvectortest.cc | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/test/unitvectortest.cc b/test/unitvectortest.cc
index e13ac291..26cfe8ac 100644
--- a/test/unitvectortest.cc
+++ b/test/unitvectortest.cc
@@ -23,6 +23,12 @@ double energy(const UnitVector<dim>& a, const FieldVector<double,dim>& b)
     return UnitVector<dim>::distance(a,b) * UnitVector<dim>::distance(a,b);
 }
 
+template <int dim>
+double energy(const FieldVector<double,dim>& a, const FieldVector<double,dim>& b)
+{
+    return UnitVector<dim>::distance(a,b) * UnitVector<dim>::distance(a,b);
+}
+
 template <int dim>
 void testDerivativesOfSquaredDistance(const UnitVector<dim>& a, const UnitVector<dim>& b)
 {
@@ -92,6 +98,31 @@ void testDerivativesOfSquaredDistance(const UnitVector<dim>& a, const UnitVector
 
     FieldMatrix<double,dim,dim> d1d2 = UnitVector<dim>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(a, b);
     
+    // finite-difference approximation
+    FieldMatrix<double,dim,dim> d1d2_fd;
+    
+    for (size_t i=0; i<dim; i++) {
+        for (size_t j=0; j<dim; j++) {
+            
+            FieldVector<double,dim> aPlus  = a.globalCoordinates();
+            FieldVector<double,dim> aMinus = a.globalCoordinates();
+            aPlus[i]  += eps;
+            aMinus[i] -= eps;
+
+            FieldVector<double,dim> bPlus  = b.globalCoordinates();
+            FieldVector<double,dim> bMinus = b.globalCoordinates();
+            bPlus[i]  += eps;
+            bMinus[i] -= eps;
+                
+            d1d2_fd[i][j] = (energy(aPlus,bPlus) + energy(aMinus,bMinus)
+                            - energy(aPlus,bMinus) - energy(aMinus,bPlus)) / (4*eps*eps);
+
+        }
+    }
+    
+    std::cout << "Analytical:\n" << d1d2 << std::endl;
+    std::cout << "FD        :\n" << d1d2_fd << std::endl;
+
     
     /////////////////////////////////////////////////////////////////////////////////////////////
     //  Test mixed third derivative with respect to first (once) and second (twice) argument
-- 
GitLab