From 9c3ff811eeb9f341eb6da1dcda0a3cb843fd27bc Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 13 Sep 2012 15:37:20 +0000
Subject: [PATCH] properly implement getting the Lagrange nodes. Works without
a priori information for all orders and dimensions
[[Imported from SVN: r8869]]
---
test/nestednesstest.cc | 54 +++++++++++++++++++++++-------------------
1 file changed, 30 insertions(+), 24 deletions(-)
diff --git a/test/nestednesstest.cc b/test/nestednesstest.cc
index 96b488e2..8aa40b89 100644
--- a/test/nestednesstest.cc
+++ b/test/nestednesstest.cc
@@ -31,34 +31,40 @@ double diameter(const std::vector<TargetSpace>& v)
return d;
}
-template <int domainDim, int elementOrder>
-std::vector<FieldVector<double,domainDim> > lagrangeNodes(const GeometryType& type)
+// The function f(x_1,x_2,x_3) := x_i, needed to get the Lagrange nodes
+// from a Dune LocalFunction
+template <int dim>
+struct CoordinateFunction
+ : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,1> >
{
- std::vector<FieldVector<double,domainDim> > result;
+ CoordinateFunction(int dir)
+ : dir_(dir)
+ {}
- // Special case d=1, q=anything, because that one is easy
- if (domainDim==1) {
-
- result.resize(elementOrder+1);
-
- for (size_t i=0; i<result.size(); i++)
- result[i] = double(i)/double(elementOrder);
-
- return result;
+ void evaluate(const FieldVector<double, dim>& x, FieldVector<double,1>& out) const {
+ out = x[dir_];
}
-
- assert(elementOrder==2);
+
+ int dir_;
+};
+
+
+
+template <int domainDim, int elementOrder>
+std::vector<FieldVector<double,domainDim> > lagrangeNodes(const GeometryType& type)
+{
PQkLocalFiniteElementCache<double,double,domainDim,elementOrder> feCache;
- result.resize(feCache.get(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(type);
-
- for (size_t i=0; i<feCache.get(type).localBasis().size(); i++) {
+
+ std::vector<FieldVector<double,domainDim> > result(feCache.get(type).localBasis().size());
+ std::vector<FieldVector<double,1> > resultCoeff(feCache.get(type).localBasis().size());
+
+ for (size_t i=0; i<domainDim; i++) {
- result[i] = refElement.position(feCache.get(type).localCoefficients().localKey(i).subEntity(),
- feCache.get(type).localCoefficients().localKey(i).codim());
+ CoordinateFunction<domainDim> cF(i);
+ feCache.get(type).localInterpolation().interpolate(cF, resultCoeff);
+
+ for (size_t j=0; j<result.size(); j++)
+ result[j][i] = resultCoeff[j];
}
@@ -110,7 +116,7 @@ void testNestedness(const LocalGeodesicFEFunction<domainDim,double,typename PQkL
assert(false);
if ( diff.infinity_norm() > eps ) {
- std::cout << className<TargetSpace>() << ": Values doe not match." << std::endl;
+ std::cout << className<TargetSpace>() << ": Values does not match." << std::endl;
std::cout << "Low order : " << loValue << std::endl;
std::cout << "High order: " << hoValue << std::endl;
assert(false);
--
GitLab