getRecoveryGradient vor higher polynomial degrees

AMDiS/src/DOFVector.cc

@@ -412,7 +412,7 @@ namespace AMDiS {
     volume.set(0.0);
 
     const BasisFunction *basFcts = feSpace->getBasisFcts();
-    int nPreDofs = feSpace->getAdmin()->getNumberOfPreDofs(VERTEX);
+    int nBasisFcts = basFcts->getNumber();
 
     DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
     WorldVector<double> grd;
@@ -421,22 +421,22 @@ namespace AMDiS {
     Mesh *mesh = feSpace->getMesh();
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(mesh, -1,
-					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET | 
+					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
 					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
 
-    ElementVector localUh(basFcts->getNumber());
+    mtl::dense_vector<double> localUh(basFcts->getNumber());
+    std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
 
     while (elInfo) {
       double det = elInfo->getDet();
-      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
       const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
       getLocalVector(elInfo->getElement(), localUh);
       basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
+      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
 
-      for (int i = 0; i < dim + 1; i++) {
-	DegreeOfFreedom dofIndex = dof[i][nPreDofs];
-	(*result)[dofIndex] += grd * det;
-	volume[dofIndex] += det;
+      for (int i = 0; i < nBasisFcts; i++) {
+	(*result)[localIndices[i]] += grd * det;
+	volume[localIndices[i]] += det;
       }
 
       elInfo = stack.traverseNext(elInfo);

AMDiS/src/io/PngReader.cc

@@ -64,6 +64,8 @@ namespace AMDiS {
     cout << "Bytes per pixel: " << bytesPerPixel << endl;
 
     double value = 0;
+    const DOFAdmin* admin = vec->getFeSpace()->getAdmin();
+
     const BasisFunction *basFcts = vec->getFeSpace()->getBasisFcts();
     int numBasFcts = basFcts->getNumber();
     DegreeOfFreedom *localIndices = new DegreeOfFreedom[numBasFcts];