From 4c936881669a30b2816e4dbf4453b4b7008bd1fd Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 13 Sep 2012 09:25:33 +0000
Subject: [PATCH] actually add the test for nestedness of GFE functions of
different order
[[Imported from SVN: r8863]]
---
test/nestednesstest.cc | 181 +++++++++++++++++++++++++++++++++++++++++
1 file changed, 181 insertions(+)
create mode 100644 test/nestednesstest.cc
diff --git a/test/nestednesstest.cc b/test/nestednesstest.cc
new file mode 100644
index 00000000..2a82f10d
--- /dev/null
+++ b/test/nestednesstest.cc
@@ -0,0 +1,181 @@
+#include <config.h>
+
+#include <fenv.h>
+#include <iostream>
+
+#include <dune/common/fvector.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/gfe/rotation.hh>
+#include <dune/gfe/realtuple.hh>
+#include <dune/gfe/unitvector.hh>
+
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include "multiindex.hh"
+#include "valuefactory.hh"
+
+const double eps = 1e-6;
+
+using namespace Dune;
+
+/** \brief Computes the diameter of a set */
+template <class TargetSpace>
+double diameter(const std::vector<TargetSpace>& v)
+{
+ double d = 0;
+ for (size_t i=0; i<v.size(); i++)
+ for (size_t j=0; j<v.size(); j++)
+ d = std::max(d, TargetSpace::distance(v[i],v[j]));
+ return d;
+}
+
+
+
+template <int domainDim, class TargetSpace>
+void testNestedness(const LocalGeodesicFEFunction<domainDim,double,typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType, TargetSpace>& loF)
+{
+ // Make higher order local gfe function
+ PQkLocalFiniteElementCache<double,double,domainDim,2> feCache;
+ typedef typename PQkLocalFiniteElementCache<double,double,domainDim,2>::FiniteElementType LocalFiniteElement;
+
+ std::vector<TargetSpace> nodalValues(feCache.get(loF.type()).localBasis().size());
+
+ // evaluate loF at the Lagrange points of the second-order function
+ const Dune::GenericReferenceElement<double,domainDim>& refElement
+ = Dune::GenericReferenceElements<double,domainDim>::general(loF.type());
+
+ for (size_t i=0; i<feCache.get(loF.type()).localBasis().size(); i++) {
+
+ nodalValues[i] = loF.evaluate(refElement.position(feCache.get(loF.type()).localCoefficients().localKey(i).subEntity(),
+ feCache.get(loF.type()).localCoefficients().localKey(i).codim()));
+
+ }
+
+
+ LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> hoF(feCache.get(loF.type()),nodalValues);
+
+
+ // A quadrature rule as a set of test points
+ int quadOrder = 10; // high, I want many many points
+
+ const Dune::QuadratureRule<double, domainDim>& quad
+ = Dune::QuadratureRules<double, domainDim>::rule(loF.type(), quadOrder);
+
+ for (size_t pt=0; pt<quad.size(); pt++) {
+
+ const Dune::FieldVector<double,domainDim>& quadPos = quad[pt].position();
+
+ // evaluate value of low-order function
+ TargetSpace loValue = loF.evaluate(quadPos);
+
+ // evaluate value of high-order function
+ TargetSpace hoValue = hoF.evaluate(quadPos);
+
+ typename TargetSpace::CoordinateType diff = loValue.globalCoordinates();
+ diff -= hoValue.globalCoordinates();
+
+ static double maxDiff = 0;
+ maxDiff = std::max(maxDiff, diff.infinity_norm());
+ std::cout << "maxDiff: " << maxDiff << std::endl;
+
+ if (maxDiff > 0.2)
+ assert(false);
+
+ if ( false and diff.infinity_norm() > eps ) {
+ std::cout << className<TargetSpace>() << ": Values doe not match." << std::endl;
+ std::cout << "Low order : " << loValue << std::endl;
+ std::cout << "High order: " << hoValue << std::endl;
+ //assert(false);
+ }
+
+ }
+
+}
+
+
+
+
+template <class TargetSpace, int domainDim>
+void test(const GeometryType& element)
+{
+ std::cout << " --- Testing " << className<TargetSpace>() << ", domain dimension: " << element.dim() << " ---" << std::endl;
+
+ std::vector<TargetSpace> testPoints;
+ ValueFactory<TargetSpace>::get(testPoints);
+
+ int nTestPoints = testPoints.size();
+ size_t nVertices = Dune::GenericReferenceElements<double,domainDim>::general(element).size(domainDim);
+
+ // Set up elements of the target space
+ std::vector<TargetSpace> corners(nVertices);
+
+ MultiIndex index(nVertices, nTestPoints);
+ int numIndices = index.cycle();
+
+ for (int i=0; i<numIndices; i++, ++index) {
+
+ for (size_t j=0; j<nVertices; j++)
+ corners[j] = testPoints[index[j]];
+
+ if (diameter(corners) > 0.5*M_PI)
+ continue;
+
+ // Make local gfe function to be tested
+ PQkLocalFiniteElementCache<double,double,domainDim,1> feCache;
+ typedef typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
+
+ LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> f(feCache.get(element),corners);
+
+ testNestedness<domainDim>(f);
+
+ }
+
+}
+
+
+int main()
+{
+ // choke on NaN -- don't enable this by default, as there are
+ // a few harmless NaN in the loopsolver
+ //feenableexcept(FE_INVALID);
+
+ std::cout << std::setw(15) << std::setprecision(12);
+
+ GeometryType element;
+
+ ////////////////////////////////////////////////////////////////
+ // Test functions on 1d elements
+ ////////////////////////////////////////////////////////////////
+ element.makeSimplex(1);
+
+ test<RealTuple<double,1>,1>(element);
+ test<UnitVector<double,2>,1>(element);
+ test<UnitVector<double,3>,1>(element);
+ test<Rotation<double,3>,1>(element);
+ test<RigidBodyMotion<double,3>,1>(element);
+ exit(0);
+ ////////////////////////////////////////////////////////////////
+ // Test functions on 2d simplex elements
+ ////////////////////////////////////////////////////////////////
+ element.makeSimplex(2);
+
+ test<RealTuple<double,1>,2>(element);
+ test<UnitVector<double,2>,2>(element);
+ test<UnitVector<double,3>,2>(element);
+ test<Rotation<double,3>,2>(element);
+ test<RigidBodyMotion<double,3>,2>(element);
+
+ ////////////////////////////////////////////////////////////////
+ // Test functions on 2d quadrilateral elements
+ ////////////////////////////////////////////////////////////////
+ element.makeCube(2);
+
+ test<RealTuple<double,1>,2>(element);
+ test<UnitVector<double,2>,2>(element);
+ test<UnitVector<double,3>,2>(element);
+ test<Rotation<double,3>,2>(element);
+ test<RigidBodyMotion<double,3>,2>(element);
+
+}
--
GitLab