From 55164cfc320bb005329d9c644a45190846c85bd2 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 2 Sep 2010 08:17:48 +0000
Subject: [PATCH] implement finite difference test for first and second
derivative
[[Imported from SVN: r6310]]
---
test/unitvectortest.cc | 102 +++++++++++++++++++++++++++++++++++++++++
1 file changed, 102 insertions(+)
diff --git a/test/unitvectortest.cc b/test/unitvectortest.cc
index b2b35ef6..8e7019a2 100644
--- a/test/unitvectortest.cc
+++ b/test/unitvectortest.cc
@@ -3,11 +3,112 @@
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/rotation.hh>
+/** \file
+ \brief Unit tests for the UnitVector class
+*/
using namespace Dune;
+const double eps = 1e-6;
+
+double square(double a)
+{
+ return a*a;
+}
+
+template <int dim>
+void testDerivativesOfSquaredDistance(const UnitVector<dim>& a, const UnitVector<dim>& b)
+{
+
+ ///////////////////////////////////////////////////////////////////
+ // Test derivative with respect to second argument
+ ///////////////////////////////////////////////////////////////////
+ typename UnitVector<dim>::EmbeddedTangentVector d2 = UnitVector<dim>::derivativeOfDistanceSquaredWRTSecondArgument(a, b);
+
+ // finite-difference approximation
+ typename UnitVector<dim>::EmbeddedTangentVector d2_fd;
+ for (size_t i=0; i<dim; i++) {
+ FieldVector<double,dim> bPlus = b.globalCoordinates();
+ FieldVector<double,dim> bMinus = b.globalCoordinates();
+ bPlus[i] += eps;
+ bMinus[i] -= eps;
+ d2_fd[i] = (square(UnitVector<dim>::distance(a,bPlus)) - square(UnitVector<dim>::distance(a,bMinus))) / (2*eps);
+ }
+
+ std::cout << "Analytical: " << d2 << std::endl;
+ std::cout << "FD : " << d2_fd << std::endl;
+
+
+ ///////////////////////////////////////////////////////////////////
+ // Test second derivative with respect to second argument
+ ///////////////////////////////////////////////////////////////////
+ FieldMatrix<double,dim,dim> d2d2 = UnitVector<dim>::secondDerivativeOfDistanceSquaredWRTSecondArgument(a, b);
+
+ // finite-difference approximation
+ FieldMatrix<double,dim,dim> d2d2_fd;
+
+ for (size_t i=0; i<dim; i++) {
+ for (size_t j=0; j<dim; j++) {
+
+ if (i==j) {
+
+ FieldVector<double,dim> bPlus = b.globalCoordinates();
+ FieldVector<double,dim> bMinus = b.globalCoordinates();
+ bPlus[i] += eps;
+ bMinus[i] -= eps;
+
+ d2d2_fd[i][i] = (square(UnitVector<dim>::distance(a,bPlus)) - 2*square(UnitVector<dim>::distance(a,b)) + square(UnitVector<dim>::distance(a,bMinus))) / (eps*eps);
+ } else {
+
+ FieldVector<double,dim> bPlusPlus = b.globalCoordinates();
+ FieldVector<double,dim> bPlusMinus = b.globalCoordinates();
+ FieldVector<double,dim> bMinusPlus = b.globalCoordinates();
+ FieldVector<double,dim> bMinusMinus = b.globalCoordinates();
+
+ bPlusPlus[i] += eps; bPlusPlus[j] += eps;
+ bPlusMinus[i] += eps; bPlusMinus[j] -= eps;
+ bMinusPlus[i] -= eps; bMinusPlus[j] += eps;
+ bMinusMinus[i] -= eps; bMinusMinus[j] -= eps;
+
+ d2d2_fd[i][j] = (square(UnitVector<dim>::distance(a,bPlusPlus)) + square(UnitVector<dim>::distance(a,bMinusMinus))
+ - square(UnitVector<dim>::distance(a,bPlusMinus)) - square(UnitVector<dim>::distance(a,bMinusPlus))) / (4*eps*eps);
+ }
+ }
+ }
+
+ std::cout << "Analytical:\n" << d2d2 << std::endl;
+ std::cout << "FD :\n" << d2d2_fd << std::endl;
+
+ //////////////////////////////////////////////////////////////////////////////
+ // Test mixed second derivative with respect to first and second argument
+ //////////////////////////////////////////////////////////////////////////////
+
+ FieldMatrix<double,dim,dim> d1d2 = UnitVector<dim>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(a, b);
+
+
+ /////////////////////////////////////////////////////////////////////////////////////////////
+ // Test mixed third derivative with respect to first (once) and second (twice) argument
+ /////////////////////////////////////////////////////////////////////////////////////////////
+
+ Tensor3<double,dim,dim,dim> d1d2d2 = UnitVector<dim>::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(a, b);
+}
+
int main()
{
+ // Set up elements of S^1
+ Dune::array<double,2> w0 = {{0,1}};
+ UnitVector<2> uv0(w0);
+ Dune::array<double,2> w1 = {{1,1}};
+ UnitVector<2> uv1(w1);
+ testDerivativesOfSquaredDistance<2>(uv0, uv1);
+
+ Dune::array<double,3> w3_0 = {{0,1,0}};
+ UnitVector<3> v3_0(w3_0);
+ Dune::array<double,3> w3_1 = {{1,1,0}};
+ UnitVector<3> v3_1(w3_1);
+ testDerivativesOfSquaredDistance<3>(v3_0, v3_1);
+
+#if 0
// Set up elements of S^1
FieldVector<double,2> v;
v[0] = 1; v[1] = 1;
@@ -27,5 +128,6 @@ int main()
std::cout << UnitVector<2>::secondDerivativeOfDistanceSquaredWRTSecondArgument(uv0, uv1) << std::endl;
std::cout << Rotation<2,double>::secondDerivativeOfDistanceSquaredWRTSecondArgument(ro0, ro1) << std::endl;
+#endif
}
--
GitLab