diff --git a/dune/solvers/transferoperators/genericmultigridtransfer.hh b/dune/solvers/transferoperators/genericmultigridtransfer.hh
index 64352907d7f8edab35a3507727e3b5ee78be1a8e..4be44e1836409e45e795befbd3169c3da6aafcc2 100644
--- a/dune/solvers/transferoperators/genericmultigridtransfer.hh
+++ b/dune/solvers/transferoperators/genericmultigridtransfer.hh
@@ -8,6 +8,9 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/bitsetvector.hh>
+#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
#include <dune/geometry/referenceelements.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh>
@@ -64,146 +67,12 @@ public:
template<class TransferOperatorType, class GridType, class field_type>
- static void setup(TransferOperatorType& matrix, const GridType& grid, int cL, int fL)
+ static void setup(TransferOperatorType& matrix, const GridType& grid, int cL, int fL, int order = 1)
{
- typedef typename TransferOperatorType::block_type TransferMatrixBlock;
- typedef typename GridType::ctype ctype;
-
- if (fL != cL+1)
- DUNE_THROW(Dune::Exception, "The two function spaces don't belong to consecutive levels!");
-
- const int dim = GridType::dimension;
-
- const typename GridType::Traits::LevelIndexSet& coarseIndexSet = grid.levelIndexSet(cL);
- const typename GridType::Traits::LevelIndexSet& fineIndexSet = grid.levelIndexSet(fL);
-
- int rows = grid.size(fL, dim);
- int cols = grid.size(cL, dim);
-
- // A factory for the shape functions
- typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
- typedef typename P1FECache::FiniteElementType FEType;
- P1FECache p1FECache;
-
- matrix.setSize(rows,cols);
- matrix.setBuildMode(TransferOperatorType::random);
-
- // ///////////////////////////////////////////
- // Determine the indices present in the matrix
- // /////////////////////////////////////////////////
- Dune::MatrixIndexSet indices(rows, cols);
-
- for (const auto& cIt : elements(grid.levelGridView(cL))) {
-
- typedef typename GridType::template Codim<0>::Entity EntityType;
-
- // Get local finite element
- const FEType& coarseBaseSet = p1FECache.get(cIt.type());
-
- const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
-
- // preallocate vector for function evaluations
- std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
-
- for (const auto& fIt : descendantElements(cIt,fL)) {
-
- if (fIt.level()==cIt.level())
- continue;
-
- const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
-
- const typename EntityType::LocalGeometry& fGeometryInFather = fIt.geometryInFather();
-
- // Get local finite element
- const FEType& fineBaseSet = p1FECache.get(fIt.type());
-
- const size_t numFineBaseFct = fineBaseSet.localBasis().size();
-
-
- for (size_t j=0; j<numFineBaseFct; j++)
- {
- const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
-
- int globalFine = fineIndexSet.subIndex(fIt, jLocalKey.subEntity(), jLocalKey.codim());
-
- Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(jLocalKey.subEntity(), jLocalKey.codim());
- Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition);
-
- // Evaluate coarse grid base functions
- coarseBaseSet.localBasis().evaluateFunction(local, values);
-
- for (size_t 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);
- }
- }
- }
-
- }
-
- }
-
- indices.exportIdx(matrix);
-
- // /////////////////////////////////////////////
- // Compute the matrix
- // /////////////////////////////////////////////
-
- for (const auto& cIt : elements(grid.levelGridView(cL))) {
-
- // Get local finite element
- const FEType& coarseBaseSet = p1FECache.get(cIt.type());
-
- const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
-
- typedef typename GridType::template Codim<0>::Entity EntityType;
-
- // preallocate vector for function evaluations
- std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
-
- for (const auto& fIt : descendantElements(cIt,fL)) {
-
- if (fIt.level()==cIt.level())
- continue;
-
- const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
-
- const typename EntityType::LocalGeometry& fGeometryInFather = fIt.geometryInFather();
-
- // Get local finite element
- const FEType& fineBaseSet = p1FECache.get(fIt.type());
-
- const size_t numFineBaseFct = fineBaseSet.localBasis().size();
-
- for (size_t j=0; j<numFineBaseFct; j++)
- {
- const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
-
- int globalFine = fineIndexSet.subIndex(fIt, jLocalKey.subEntity(), jLocalKey.codim());
-
- Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(jLocalKey.subEntity(), jLocalKey.codim());
- Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition);
-
- // Evaluate coarse grid base functions
- coarseBaseSet.localBasis().evaluateFunction(local, values);
-
- for (size_t 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());
-
- matrix[globalFine][globalCoarse] = Dune::ScaledIdentityMatrix<ctype,TransferMatrixBlock::rows>(values[i]);
- }
- }
- }
- }
- }
+ // call function overload with two gridViews
+ auto fineGridView = grid.levelGridView(fL);
+ auto coarseGridView = grid.levelGridView(cL);
+ setup(matrix, coarseGridView, fineGridView, order);
}
@@ -213,161 +82,15 @@ public:
\param matrix The matrix object to assemble into
*/
- template<class TransferOperatorType, class GridViewCoarse, class GridViewFine, class field_type>
- static void setup(TransferOperatorType& matrix, const GridViewCoarse& gvCoarse, const GridViewFine& gvFine)
+ template<class TransferOperatorType, class GridViewCoarse, class GridViewFine, class field_type = double>
+ static void setup(TransferOperatorType& matrix, const GridViewCoarse& gvCoarse, const GridViewFine& gvFine, int order = 1)
{
- typedef typename TransferOperatorType::block_type TransferMatrixBlock;
- typedef typename GridViewCoarse::ctype ctype;
-
- const int dim = GridViewCoarse::dimension;
-
- const typename GridViewCoarse::IndexSet& coarseIndexSet = gvCoarse.indexSet();
- const typename GridViewFine::IndexSet& fineIndexSet = gvFine.indexSet();
-
- int rows = gvFine.size(dim);
- int cols = gvCoarse.size(dim);
+ // create Lagrange bases of given order for fine and coarse grid level
+ auto fineBasis = Dune::Functions::BasisFactory::makeBasis(gvFine, Dune::Functions::BasisFactory::lagrange(order));
+ auto coarseBasis = Dune::Functions::BasisFactory::makeBasis(gvCoarse,Dune::Functions::BasisFactory::lagrange(order));
- // A factory for the shape functions
- typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
- typedef typename P1FECache::FiniteElementType FEType;
- P1FECache p1FECache;
-
- matrix.setSize(rows,cols);
- matrix.setBuildMode(TransferOperatorType::random);
-
- // ///////////////////////////////////////////
- // Determine the indices present in the matrix
- // /////////////////////////////////////////////////
-
- Dune::MatrixIndexSet indices(rows, cols);
-
- for (const auto& fIt : elements(gvFine)) {
-
- const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
-
- // Get local finite element
- const FEType& fineBaseSet = p1FECache.get(fIt.type());;
-
- const size_t numFineBaseFct = fineBaseSet.localBasis().size();
-
- std::vector<Dune::FieldVector<ctype,dim> > local(numFineBaseFct);
-
- for (size_t i=0; i<numFineBaseFct; i++)
- {
- const Dune::LocalKey& iLocalKey = fineBaseSet.localCoefficients().localKey(i);
- local[i] = fineRefElement.position(iLocalKey.subEntity(), iLocalKey.codim());
- }
-
- // Get ancestor element in the coarse grid view, and the local position there.
- auto ancestor = fIt;
-
- while (not gvCoarse.contains(ancestor) && ancestor.level() != 0 ) {
-
- auto geometryInFather = ancestor.geometryInFather();
-
- for (size_t i=0; i<numFineBaseFct; i++)
- local[i] = geometryInFather.global(local[i]);
- ancestor = ancestor.father();
-
- }
-
- assert(gvCoarse.contains(ancestor));
-
- for (size_t j=0; j<numFineBaseFct; j++)
- {
- // Get local finite element
- const FEType& coarseBaseSet = p1FECache.get(ancestor.type());
-
- const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
-
- // preallocate vector for function evaluations
- std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
-
- // Evaluate coarse grid base functions
- coarseBaseSet.localBasis().evaluateFunction(local[j], values);
-
- const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
- int globalFine = fineIndexSet.subIndex(fIt, jLocalKey.subEntity(), jLocalKey.codim());
-
- for (size_t i=0; i<numCoarseBaseFct; i++)
- {
- if (values[i] > 0.001)
- {
- const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
- int globalCoarse = coarseIndexSet.subIndex(ancestor, iLocalKey.subEntity(), iLocalKey.codim());
- indices.add(globalFine, globalCoarse);
- }
- }
- }
-
- }
-
- indices.exportIdx(matrix);
-
- // /////////////////////////////////////////////
- // Compute the matrix
- // /////////////////////////////////////////////
- for (const auto& fIt : elements(gvFine)) {
-
-
- const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
-
- // Get local finite element
- const FEType& fineBaseSet = p1FECache.get(fIt.type());;
-
- const size_t numFineBaseFct = fineBaseSet.localBasis().size();
-
- std::vector<Dune::FieldVector<ctype,dim> > local(numFineBaseFct);
-
- for (size_t i=0; i<numFineBaseFct; i++)
- {
- const Dune::LocalKey& iLocalKey = fineBaseSet.localCoefficients().localKey(i);
- local[i] = fineRefElement.position(iLocalKey.subEntity(), iLocalKey.codim());
- }
-
- // Get ancestor element in the coarse grid view, and the local position there.
- auto ancestor = fIt;
-
- while (not gvCoarse.contains(ancestor) && ancestor.level() != 0 ) {
-
- auto geometryInFather = ancestor.geometryInFather();
-
- for (size_t i=0; i<numFineBaseFct; i++)
- local[i] = geometryInFather.global(local[i]);
- ancestor = ancestor.father();
-
- }
-
- assert(gvCoarse.contains(ancestor));
-
- for (size_t j=0; j<numFineBaseFct; j++)
- {
- // Get local finite element
- const FEType& coarseBaseSet = p1FECache.get(ancestor.type());
-
- const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
-
- // preallocate vector for function evaluations
- std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
-
- // Evaluate coarse grid base functions
- coarseBaseSet.localBasis().evaluateFunction(local[j], values);
-
- const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
- int globalFine = fineIndexSet.subIndex(fIt, jLocalKey.subEntity(), jLocalKey.codim());
-
- for (size_t i=0; i<numCoarseBaseFct; i++)
- {
- if (values[i] > 0.001)
- {
- const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
- int globalCoarse = coarseIndexSet.subIndex(ancestor, iLocalKey.subEntity(), iLocalKey.codim());
-
- matrix[globalFine][globalCoarse] = Dune::ScaledIdentityMatrix<ctype,TransferMatrixBlock::rows>(values[i]);
- }
- }
- }
- }
+ // call generic method to assemble the transfer matrix
+ assembleGlobalBasisTransferMatrix(matrix,coarseBasis,fineBasis);
}
//! Multiply the vector f from the right to the prolongation matrix