diff --git a/dune-solvers/transferoperators/genericmultigridtransfer.hh b/dune-solvers/transferoperators/genericmultigridtransfer.hh index ea4878181cadc874b6eb86e20ca4cd6607688fbb..526c844eb7674d75f2f0c05225a4d936745da05b 100644 --- a/dune-solvers/transferoperators/genericmultigridtransfer.hh +++ b/dune-solvers/transferoperators/genericmultigridtransfer.hh @@ -135,6 +135,9 @@ public: const int numCoarseBaseFct = coarseBaseSet->localBasis().size(); + // preallocate vector for function evaluations + std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct); + HierarchicIterator fIt = cIt->hbegin(fL); HierarchicIterator fEndIt = cIt->hend(fL); @@ -153,28 +156,28 @@ public: const int numFineBaseFct = fineBaseSet->localBasis().size(); - for (int i=0; i<numCoarseBaseFct; i++) { - - int globalCoarse = coarseIndexSet.subIndex(*cIt,i,dim); - for (int j=0; j<numFineBaseFct; j++) { + for (int j=0; j<numFineBaseFct; j++) + { + const Dune::LocalKey& jLocalKey = fineBaseSet->localCoefficients().localKey(j); - int globalFine = fineIndexSet.subIndex(*fIt,j,dim); + int globalFine = fineIndexSet.subIndex(*fIt, jLocalKey.subEntity(), jLocalKey.codim()); - Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(j,dim); - Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition); + Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(jLocalKey.subEntity(), jLocalKey.codim()); + Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition); - // Evaluate coarse grid base function - std::vector<Dune::FieldVector<field_type,1> > values; - coarseBaseSet->localBasis().evaluateFunction(local, values); - double value = values[i]; + // Evaluate coarse grid base functions + coarseBaseSet->localBasis().evaluateFunction(local, values); - if (value > 0.001) + for (int i=0; i<numCoarseBaseFct; i++) + { + if (values[i] > 0.001) + { + const Dune::LocalKey& iLocalKey = coarseBaseSet->localCoefficients().localKey(i); + int globalCoarse = coarseIndexSet.subIndex(*cIt, iLocalKey.subEntity(), iLocalKey.codim()); indices.add(globalFine, globalCoarse); - + } } - - } } @@ -197,6 +200,9 @@ public: typedef typename GridType::template Codim<0>::Entity EntityType; typedef typename EntityType::HierarchicIterator HierarchicIterator; + // preallocate vector for function evaluations + std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct); + HierarchicIterator fIt = cIt->hbegin(fL); HierarchicIterator fEndIt = cIt->hend(fL); @@ -215,42 +221,31 @@ public: const int numFineBaseFct = fineBaseSet->localBasis().size(); - for (int i=0; i<numCoarseBaseFct; i++) { - - int globalCoarse = coarseIndexSet.subIndex(*cIt,i,dim); - - for (int j=0; j<numFineBaseFct; j++) { + for (int j=0; j<numFineBaseFct; j++) + { + const Dune::LocalKey& jLocalKey = fineBaseSet->localCoefficients().localKey(j); - int globalFine = fineIndexSet.subIndex(*fIt,j,dim); + int globalFine = fineIndexSet.subIndex(*fIt, jLocalKey.subEntity(), jLocalKey.codim()); - // Evaluate coarse grid base function at the location of the fine grid dof + Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(jLocalKey.subEntity(), jLocalKey.codim()); + Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition); - // first determine local fine grid dof position - Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(j,dim); - Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition); + // Evaluate coarse grid base functions + coarseBaseSet->localBasis().evaluateFunction(local, values); - // Evaluate coarse grid base function - std::vector<Dune::FieldVector<field_type,1> > values; - coarseBaseSet->localBasis().evaluateFunction(local, values); - double value = values[i]; + for (int i=0; i<numCoarseBaseFct; i++) + { + if (values[i] > 0.001) + { + const Dune::LocalKey& iLocalKey = coarseBaseSet->localCoefficients().localKey(i); + int globalCoarse = coarseIndexSet.subIndex(*cIt, iLocalKey.subEntity(), iLocalKey.codim()); - // The following conditional is a hack: evaluating the coarse - // grid base function will often return 0. However, we don't - // want explicit zero entries in our prolongation matrix. Since - // the whole code works for P1-elements only anyways, we know - // that value can only be 0, 0.5, or 1. Thus testing for nonzero - // by testing > 0.001 is safe. - if (value > 0.001) { TransferMatrixBlock matValue = identity; - matValue *= value; + matValue *= values[i]; matrix[globalFine][globalCoarse] = matValue; } - } - - } - } }