diff --git a/dune/solvers/transferoperators/genericmultigridtransfer.hh b/dune/solvers/transferoperators/genericmultigridtransfer.hh index c08a329ce770eeca69eaf5eaa4cfa17d002616da..6adca9fed76d2cc394339ba5e34abe3728575598 100644 --- a/dune/solvers/transferoperators/genericmultigridtransfer.hh +++ b/dune/solvers/transferoperators/genericmultigridtransfer.hh @@ -603,7 +603,7 @@ public: //! Restrict a vector valued bitfield from the fine onto the coarse grid template <class TransferOperatorType, class BitVectorType> - static void restrictBitField(const TransferOperatorType& matrix, const BitVectorType& f, BitVectorType& t) + static void restrictBitField(const TransferOperatorType& matrix, const BitVectorType& f, BitVectorType& t, bool onlyToFathers=false) { if (f.size() != matrix.N()) DUNE_THROW(Dune::Exception, "Fine grid bitfield has " << f.size() << " entries " @@ -630,7 +630,8 @@ public: { for(size_t i=0; i<fineBits.size(); ++i) if (fineBits.test(i)) - t[cIt.index()][i] = true; + if (!onlyToFathers || diagonalIsOne(*cIt, i)) + t[cIt.index()][i] = true; } } @@ -640,7 +641,7 @@ public: //! Restrict a vector valued bitfield from the fine onto the coarse grid template <class TransferOperatorType, class BitVectorType> - static void restrictBitField(const TransferOperatorType& matrix, const BitVectorType& f, BitVectorType& t, int virtualBlockSize) + static void restrictBitField(const TransferOperatorType& matrix, const BitVectorType& f, BitVectorType& t, int virtualBlockSize, bool onlyToFathers=false) { if (virtualBlockSize<0) virtualBlockSize = f.size()/matrix.N(); @@ -681,7 +682,8 @@ public: const typename BitVectorType::const_reference fineBits = f[rowIdx*virtualBlockSize+i]; for(size_t j=0; j<fineBits.size(); ++j) if (fineBits.test(j)) - t[cIt.index()*virtualBlockSize+i][j] = true; + if (!onlyToFathers || diagonalIsOne(*cIt, j)) + t[cIt.index()*virtualBlockSize+i][j] = true; } } @@ -694,35 +696,7 @@ public: template <class TransferOperatorType, class BitVectorType> static void restrictBitFieldToFathers(const TransferOperatorType& matrix, const BitVectorType& f, BitVectorType& t) { - if (f.size() != matrix.N()) - DUNE_THROW(Dune::Exception, "Fine grid bitfield has " << f.size() << " entries " - << "but the interpolation matrix has " << matrix.N() << " rows!"); - - t.resize(matrix.M()); - t.unsetAll(); - - typedef typename TransferOperatorType::row_type RowType; - typedef typename RowType::ConstIterator ColumnIterator; - - for (size_t rowIdx=0; rowIdx<matrix.N(); rowIdx++) { - const typename BitVectorType::const_reference fineBits = f[rowIdx]; - - if (fineBits.none()) - continue; - - const RowType& row = matrix[rowIdx]; - - ColumnIterator cIt = row.begin(); - ColumnIterator cEndIt = row.end(); - - for(; cIt!=cEndIt; ++cIt) - { - for(size_t i=0; i<fineBits.size(); ++i) - if (fineBits.test(i)) - if (diagonalIsOne(*cIt, i)) - t[cIt.index()][i] = true; - } - } + restrictBitField(matrix, f, t, true); } @@ -730,50 +704,7 @@ public: template <class TransferOperatorType, class BitVectorType> static void restrictBitFieldToFathers(const TransferOperatorType& matrix, const BitVectorType& f, BitVectorType& t, int virtualBlockSize) { - if (virtualBlockSize<0) - virtualBlockSize = f.size()/matrix.N(); - - if (f.size() != (unsigned int)matrix.N()*virtualBlockSize) - DUNE_THROW(Dune::Exception, "Fine grid bitfield has " << f.size() << " entries " - << "but the interpolation matrix has " << matrix.N() << " rows!"); - - t.resize(matrix.M()*virtualBlockSize); - t.unsetAll(); - - typedef typename TransferOperatorType::row_type RowType; - typedef typename RowType::ConstIterator ColumnIterator; - - for (size_t rowIdx=0; rowIdx<matrix.N(); rowIdx++) { - - bool someBitSet = false; - for(int i=0; i<virtualBlockSize; ++i) - { - if (f[rowIdx*virtualBlockSize+i].any()) - { - someBitSet = true; - continue; - } - } - if (not(someBitSet)) - continue; - - const RowType& row = matrix[rowIdx]; - - ColumnIterator cIt = row.begin(); - ColumnIterator cEndIt = row.end(); - - for(; cIt!=cEndIt; ++cIt) - { - for(int i=0; i<virtualBlockSize; ++i) - { - const typename BitVectorType::const_reference fineBits = f[rowIdx*virtualBlockSize+i]; - for(size_t j=0; j<fineBits.size(); ++j) - if (fineBits.test(j)) - if (diagonalIsOne(*cIt, j)) - t[cIt.index()*virtualBlockSize+i][j] = true; - } - } - } + restrictBitField(matrix, f, t, virtualBlockSize, true); }