diff --git a/fem/transfer/mgtransfer.cc b/fem/transfer/mgtransfer.cc index 084f4bd8c2a91afe9be16d5aca1d3c6ec48b67ed..1c42cf17d1cbcdf88ea4b7b5d4822d77c78c3b1a 100644 --- a/fem/transfer/mgtransfer.cc +++ b/fem/transfer/mgtransfer.cc @@ -15,8 +15,8 @@ void Dune::MGTransfer<DiscFuncType>::setup(const FunctionSpaceType& fS, int cL, typedef typename FunctionSpaceType::GridType GridType; - int rows = fS.size(coarseLevel); - int cols = fS.size(fineLevel); + int rows = fS.size(fineLevel); + int cols = fS.size(coarseLevel); GridType& grid = fS.getGrid(); @@ -78,7 +78,7 @@ void Dune::MGTransfer<DiscFuncType>::setup(const FunctionSpaceType& fS, int cL, //std::cout << "value: " << value << "\n"; - matrix_.set(globalCoarse, globalFine, value[0]); + matrix_.set(globalFine, globalCoarse, value[0]); } @@ -88,22 +88,20 @@ void Dune::MGTransfer<DiscFuncType>::setup(const FunctionSpaceType& fS, int cL, } - //matrix_.print(std::cout); } template<class DiscFuncType> -void Dune::MGTransfer<DiscFuncType>::restrict (const DiscFuncType & f, DiscFuncType& t) const +void Dune::MGTransfer<DiscFuncType>::prolong(const DiscFuncType& f, DiscFuncType& t) const { - assert(t.getFunctionSpace().size(coarseLevel) == matrix_.rows()); - assert(f.getFunctionSpace().size(fineLevel) == matrix_.cols()); + assert(t.getFunctionSpace().size(fineLevel) == matrix_.rows()); + assert(f.getFunctionSpace().size(coarseLevel) == matrix_.cols()); typedef typename DiscFuncType::DofIteratorType DofIteratorType; typedef typename SparseRowMatrix<double>::ColumnIterator ColumnIterator; - - DofIteratorType tIt = t.dbegin( coarseLevel ); - DofIteratorType fIt = f.dbegin( fineLevel ); + DofIteratorType tIt = t.dbegin( fineLevel ); + DofIteratorType fIt = f.dbegin( coarseLevel ); for(int row=0; row<matrix_.rows(); row++) { @@ -121,18 +119,18 @@ void Dune::MGTransfer<DiscFuncType>::restrict (const DiscFuncType & f, DiscFuncT } template<class DiscFuncType> -void Dune::MGTransfer<DiscFuncType>::prolong(const DiscFuncType& f, DiscFuncType& t) const +void Dune::MGTransfer<DiscFuncType>::restrict (const DiscFuncType & f, DiscFuncType& t) const { - assert(f.getFunctionSpace().size(coarseLevel) == matrix_.rows()); - assert(t.getFunctionSpace().size(fineLevel) == matrix_.cols()); + assert(f.getFunctionSpace().size(fineLevel) == matrix_.rows()); + assert(t.getFunctionSpace().size(coarseLevel) == matrix_.cols()); typedef typename DiscFuncType::DofIteratorType DofIteratorType; typedef typename SparseRowMatrix<double>::ColumnIterator ColumnIterator; - t.clearLevel(fineLevel); + t.clearLevel(coarseLevel); - DofIteratorType tIt = t.dbegin( fineLevel ); - DofIteratorType fIt = f.dbegin( coarseLevel ); + DofIteratorType tIt = t.dbegin( coarseLevel ); + DofIteratorType fIt = f.dbegin( fineLevel ); for (int row=0; row<matrix_.rows(); row++) { @@ -145,6 +143,7 @@ void Dune::MGTransfer<DiscFuncType>::prolong(const DiscFuncType& f, DiscFuncType tIt[cIt.col()] += fIt[row] * (*cIt); } + } @@ -152,59 +151,35 @@ template<class DiscFuncType> Dune::SparseRowMatrix<double> Dune::MGTransfer<DiscFuncType>:: galerkinRestrict(const Dune::SparseRowMatrix<double>& fineMat) const { - typedef typename SparseRowMatrix<double>::ColumnIterator ColumnIterator; - //return matrix_.applyFromLeftAndRightTo(fineMat); - - // Hack: We need the transposed matrix - - SparseRowMatrix<double> transpose(matrix_.cols(), matrix_.rows(), matrix_.numNonZeros()); - - for (int row=0; row<matrix_.rows(); row++) { - - /** \todo Remove casts */ - ColumnIterator cIt = const_cast<SparseRowMatrix<double>&>(matrix_).template rbegin(row); - ColumnIterator cEndIt = const_cast<SparseRowMatrix<double>&>(matrix_).template rend(row); - - for(; cIt!=cEndIt; ++cIt) - transpose.set(cIt.col(), row, *cIt); - } + typedef typename SparseRowMatrix<double>::ColumnIterator ColumnIterator; - SparseRowMatrix<double> result(matrix_.rows() ,matrix_.rows() , fineMat.numNonZeros()); + SparseRowMatrix<double> result(matrix_.cols() ,matrix_.cols() , fineMat.numNonZeros()); - //for (v=FIRSTVECTOR(FineGrid); v!=NULL; v=SUCCVC(v)) for (int row=0; row<fineMat.rows(); row++) { - // for (m=VSTART(v); m!=NULL; m=MNEXT(m)) { - // w = MDEST(m); ColumnIterator cIt = const_cast<SparseRowMatrix<double>&>(fineMat).template rbegin(row); ColumnIterator cEndIt = const_cast<SparseRowMatrix<double>&>(fineMat).template rend(row); for(; cIt!=cEndIt; ++cIt) { - //mvalue = MVALUE(m,mc); double mvalue = *cIt; - ColumnIterator tciIt = transpose.template rbegin(row); - ColumnIterator tciEndIt = transpose.template rend(row); + ColumnIterator tciIt = const_cast<SparseRowMatrix<double>&>(matrix_).template rbegin(row); + ColumnIterator tciEndIt = const_cast<SparseRowMatrix<double>&>(matrix_).template rend(row); - //for (im=VISTART(v); im!= NULL; im = NEXT(im)) { for (; tciIt!=tciEndIt; ++tciIt) { - //fac = mvalue*MVALUE(im,0); double fac = mvalue* (*tciIt); - ColumnIterator tcjIt = transpose.template rbegin(cIt.col()); - ColumnIterator tcjEndIt = transpose.template rend(cIt.col()); + ColumnIterator tcjIt = const_cast<SparseRowMatrix<double>&>(matrix_).template rbegin(cIt.col()); + ColumnIterator tcjEndIt = const_cast<SparseRowMatrix<double>&>(matrix_).template rend(cIt.col()); - //for (jm=VISTART(v); jm!= NULL; jm = NEXT(jm)) { for (; tcjIt!=tcjEndIt; ++tcjIt) { - // jv = MDEST(jm); - // cm = GetMatrix(iv,jv); - // MVALUE(cm,mc) += - // fac * MVALUE(jm,0); - result.add( tciIt.col(), tcjIt.col(), fac* (*tcjIt)); + //result.add( tciIt.col(), tcjIt.col(), fac* (*tcjIt)); + result.add( tcjIt.col(), tciIt.col(), fac* (*tcjIt)); + } }