From 10890b074aa5de3a51e5ba126f2b5bec3239a720 Mon Sep 17 00:00:00 2001 From: Oliver Sander <sander@igpm.rwth-aachen.de> Date: Sun, 7 Mar 2010 09:36:06 +0000 Subject: [PATCH] new test averagedistanceassemblertest [[Imported from SVN: r5669]] --- test/Makefile.am | 4 +- test/averagedistanceassemblertest.cc | 155 +++++++++++++++++++++++++++ 2 files changed, 158 insertions(+), 1 deletion(-) create mode 100644 test/averagedistanceassemblertest.cc diff --git a/test/Makefile.am b/test/Makefile.am index 9cd17092..d11177b9 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -5,7 +5,7 @@ LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) AM_CPPFLAGS += $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) -Wall check_PROGRAMS = frameinvariancetest rotationtest fdcheck localgeodesicfefunctiontest \ - harmonicenergytest + harmonicenergytest averagedistanceassemblertest frameinvariancetest_SOURCES = frameinvariancetest.cc @@ -17,6 +17,8 @@ localgeodesicfefunctiontest_SOURCES = localgeodesicfefunctiontest.cc harmonicenergytest_SOURCES = harmonicenergytest.cc +averagedistanceassemblertest_SOURCES = averagedistanceassemblertest.cc + # don't follow the full GNU-standard # we need automake 1.5 AUTOMAKE_OPTIONS = foreign 1.5 diff --git a/test/averagedistanceassemblertest.cc b/test/averagedistanceassemblertest.cc new file mode 100644 index 00000000..62c4f1d7 --- /dev/null +++ b/test/averagedistanceassemblertest.cc @@ -0,0 +1,155 @@ +#include <config.h> + +#include <iostream> + +#include <dune/common/fvector.hh> +#include <dune/grid/common/quadraturerules.hh> + +#include <dune/src/rotation.hh> +#include <dune/src/realtuple.hh> +#include <dune/src/unitvector.hh> + +#include <dune/src/localgeodesicfefunction.hh> + +// Domain dimension +const int dim = 2; + +using namespace Dune; + +/** \brief Test whether interpolation is invariant under permutation of the simplex vertices + */ +template <class TargetSpace> +void testPoint(const std::vector<TargetSpace>& corners, + const std::vector<double>& weights, + const TargetSpace& argument) +{ + // create the assembler + AverageDistanceAssembler<TargetSpace> assembler(corners, weights); + + // test the functional + double value = assembler.value(argument); + assert(!std::isnan(value)); + assert(value >= 0); + + // test the gradient + typename TargetSpace::TangentVector gradient; + assembler.assembleGradient(argument, gradient); + + // test the hessian + FieldMatrix<double, TargetSpace::TangentVector::size, TargetSpace::TangentVector::size> hessian; + FieldMatrix<double, TargetSpace::TangentVector::size, TargetSpace::TangentVector::size> hessianApproximation; + + assembler.assembleHessian(argument, hessian); + assembler.assembleHessianApproximation(argument, hessianApproximation); + + for (size_t i=0; i<hessian.N(); i++) + for (size_t j=0; j<hessian.M(); j++) { + assert(!std::isnan(hessian[i][j])); + assert(!std::isnan(hessianApproximation[i][j])); + } + + FieldMatrix<double, TargetSpace::TangentVector::size, TargetSpace::TangentVector::size> diff = hessian; + diff -= hessianApproximation; + std::cout << "Matrix inf diff: " << diff.infinity_norm() << std::endl; +} + + +template <class TargetSpace> +void testWeightSet(const std::vector<TargetSpace>& corners, + const TargetSpace& argument) +{ + // A quadrature rule as a set of test points + int quadOrder = 3; + + const Dune::QuadratureRule<double, dim>& quad + = Dune::QuadratureRules<double, dim>::rule(GeometryType(GeometryType::simplex,dim), quadOrder); + + for (size_t pt=0; pt<quad.size(); pt++) { + + const Dune::FieldVector<double,dim>& quadPos = quad[pt].position(); + + // local to barycentric coordinates + std::vector<double> weights(dim+1); + weights[0] = 1; + for (int i=0; i<dim; i++) { + weights[0] -= quadPos[i]; + weights[i+1] = quadPos[i]; + } + + testPoint(corners, weights, argument); + + } + +} + + +void testRealTuples() +{ + typedef RealTuple<1> TargetSpace; + + std::vector<TargetSpace> corners = {TargetSpace(1), + TargetSpace(2), + TargetSpace(3)}; + + TargetSpace argument = corners[0]; + testWeightSet(corners, argument); + argument = corners[1]; + testWeightSet(corners, argument); + argument = corners[2]; + testWeightSet(corners, argument); +} + +void testUnitVectors() +{ + typedef UnitVector<3> TargetSpace; + + std::vector<TargetSpace> corners(dim+1); + + FieldVector<double,3> input; + input[0] = 1; input[1] = 0; input[2] = 0; + corners[0] = input; + input[0] = 0; input[1] = 1; input[2] = 0; + corners[1] = input; + input[0] = 0; input[1] = 0; input[2] = 1; + corners[2] = input; + + TargetSpace argument = corners[0]; + testWeightSet(corners, argument); + argument = corners[1]; + testWeightSet(corners, argument); + argument = corners[2]; + testWeightSet(corners, argument); +} + +void testRotations() +{ + typedef Rotation<3,double> TargetSpace; + + FieldVector<double,3> xAxis(0); + xAxis[0] = 1; + FieldVector<double,3> yAxis(0); + yAxis[1] = 1; + FieldVector<double,3> zAxis(0); + zAxis[2] = 1; + + + std::vector<TargetSpace> corners(dim+1); + corners[0] = Rotation<3,double>(xAxis,0.1); + corners[1] = Rotation<3,double>(yAxis,0.1); + corners[2] = Rotation<3,double>(zAxis,0.1); + + TargetSpace argument = corners[0]; + testWeightSet(corners, argument); + argument = corners[1]; + testWeightSet(corners, argument); + argument = corners[2]; + testWeightSet(corners, argument); +} + + +int main() +{ + testRealTuples(); + testUnitVectors(); + //testRotations(); +} -- GitLab