diff --git a/test/unitvectortest.cc b/test/unitvectortest.cc index e13ac291bf795dda11d3c3b46dd93ce085c7e599..26cfe8ac37c32ec3140dfdc041d2f4dc0ca1b96a 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