From 2b0449b0312b884328b30cb0f254c856e5638ec5 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 27 May 2013 13:37:05 +0000
Subject: [PATCH] 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]]
---
.../genericmultigridtransfer.hh | 173 ++++++++++++++++++
1 file changed, 173 insertions(+)
diff --git a/dune/solvers/transferoperators/genericmultigridtransfer.hh b/dune/solvers/transferoperators/genericmultigridtransfer.hh
index fdc9522..133817d 100644
--- a/dune/solvers/transferoperators/genericmultigridtransfer.hh
+++ b/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)
--
GitLab