Set up a transfer operator between two arbitrary grid views of one grid

This is needed, e.g., to implement two-level Schwarz methods.
It includes the previous 'setup' method for two adjacent level views
as a special case.  I haven't touched the old method, but we may
consider removing it.

[[Imported from SVN: r8759]]

dune/solvers/transferoperators/genericmultigridtransfer.hh

@@ -223,6 +223,179 @@ public:
     }
 
 
+    /** \brief Set up transfer operator between two arbitrary grid views of the same grid 
+     
+       The 'fine' grid view gvFine is expected to be a refinement of the 'coarse' one gvCoarse.
+     
+       \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)
+    {
+        typedef typename TransferOperatorType::block_type TransferMatrixBlock;
+        typedef typename GridViewCoarse::ctype ctype;
+        typedef typename GridViewCoarse::Grid GridType;
+
+        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);
+
+        // 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);
+
+        typedef typename GridViewFine::template Codim<0>::Iterator ElementIterator;
+
+        ElementIterator fIt    = gvFine.template begin<0>();
+        ElementIterator fEndIt = gvFine.template end<0>();
+
+
+        // ///////////////////////////////////////////
+        // Determine the indices present in the matrix
+        // /////////////////////////////////////////////////
+        Dune::MatrixIndexSet indices(rows, cols);
+
+        for (; fIt != fEndIt; ++fIt) {
+
+            const Dune::ReferenceElement<ctype,dim>& 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.
+            typename GridViewFine::template Codim<0>::EntityPointer ancestor(fIt);
+            
+            while (not gvCoarse.contains(ancestor) && ancestor->level() != 0 ) {
+                
+                typename GridType::template Codim<0>::LocalGeometry 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 (fIt = gvFine.template begin<0>(); fIt != fEndIt; ++fIt) {
+
+
+            const Dune::ReferenceElement<ctype,dim>& 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.
+            typename GridViewFine::template Codim<0>::EntityPointer ancestor(fIt);
+            
+            while (not gvCoarse.contains(ancestor) && ancestor->level() != 0 ) {
+                
+                typename GridType::template Codim<0>::LocalGeometry 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]);
+                    }
+                }
+            }
+
+        }
+
+    }
+
     //! Multiply the vector f from the right to the prolongation matrix
     template<class TransferOperatorType, class DiscFuncType>
     static void prolong(const TransferOperatorType& matrix, const DiscFuncType& f, DiscFuncType &t)