truncatedcompressedmgtransfer.cc

template<class DiscFuncType, class BitVectorType, class OperatorType>
void TruncatedCompressedMGTransfer<DiscFuncType, BitVectorType, OperatorType>::prolong(const DiscFuncType& f, DiscFuncType& t,
                                              const BitVectorType& critical) const
{
    if (f.size() != this->matrix_->M())
        DUNE_THROW(Dune::Exception, "Number of entries in the coarse grid vector is not equal "
                   << "to the number of columns of the prolongation matrix!");

    if (((unsigned int)this->matrix_->N()) != critical.size())
        DUNE_THROW(Dune::Exception, "Number of entries in the critical is not equal "
                   << "to the number of rows of the prolongation matrix!");

    t.resize(this->matrix_->N());

    typedef typename DiscFuncType::Iterator      Iterator;
    typedef typename DiscFuncType::ConstIterator ConstIterator;

    typedef typename TransferOperatorType::row_type RowType;
    typedef typename RowType::ConstIterator ColumnIterator;
    
    Iterator tIt      = t.begin();
    ConstIterator fIt = f.begin(); 

    for(int rowIdx=0; rowIdx<this->matrix_->N(); rowIdx++) {

        const RowType& row = (*this->matrix_)[rowIdx];

        *tIt = 0.0;
        ColumnIterator cIt    = row.begin();
        ColumnIterator cEndIt = row.end();

        for(; cIt!=cEndIt; ++cIt) {

            // The following lines are a matrix-vector loop with a multiple of the
            // identity matrix, but rows belonging to critical dofs are left out
            for (int i=0; i<blocksize; i++) {

                if (!critical[rowIdx][i])
                    (*tIt)[i] += (*cIt)[0][0] * f[cIt.index()][i];

            }

        }
        
        ++tIt;
    }

}

template<class DiscFuncType, class BitVectorType, class OperatorType>
void TruncatedCompressedMGTransfer<DiscFuncType, BitVectorType, OperatorType>::restrict(const DiscFuncType& f, DiscFuncType& t,
                                              const BitVectorType& critical) const
{
    if (f.size() != this->matrix_->N())
        DUNE_THROW(Dune::Exception, "Fine grid vector has " << f.size() << " entries "
                   << "but the interpolation matrix has " << this->matrix_->N() << " rows!");


    if (((unsigned int)this->matrix_->N()) != critical.size())
        DUNE_THROW(Dune::Exception, "Number of entries in the critical is not equal "
                   << "to the number of rows of the prolongation matrix!");

    t.resize(this->matrix_->M());
    t = 0;

    typedef typename DiscFuncType::Iterator      Iterator;
    typedef typename DiscFuncType::ConstIterator ConstIterator;

    typedef typename TransferOperatorType::row_type RowType;
    typedef typename RowType::ConstIterator ColumnIterator;

    for (int rowIdx=0; rowIdx<this->matrix_->N(); rowIdx++) {
        
        const RowType& row = (*this->matrix_)[rowIdx];

        ColumnIterator cIt    = row.begin();
        ColumnIterator cEndIt = row.end();
        
        for(; cIt!=cEndIt; ++cIt) {
            
            // The following lines are a matrix-vector loop, but rows belonging
            // to critical dofs are left out
            typename DiscFuncType::block_type& tEntry = t[cIt.index()];
            for (int i=0; i<blocksize; i++) {
                
                if (!critical[rowIdx][i])
                    tEntry[i] += (*cIt)[0][0] * f[rowIdx][i];

            }

        }

    }
 
}


template<class DiscFuncType, class BitVectorType, class OperatorType>
void TruncatedCompressedMGTransfer<DiscFuncType, BitVectorType, OperatorType>::
galerkinRestrict(const OperatorType& fineMat, OperatorType& coarseMat,
                 const BitVectorType& critical) const
{

    if (this->recompute_ != NULL && this->recompute_->size() != (unsigned int)this->matrix_->M())
        DUNE_THROW(Dune::Exception, "The size of the recompute_-bitfield doesn't match the "
                   << "size of the coarse grid space!");

  // ////////////////////////
  // Nonsymmetric case
  // ////////////////////////
  typedef typename OperatorType::row_type RowType;
  typedef typename RowType::Iterator ColumnIterator;
  typedef typename RowType::ConstIterator ConstColumnIterator;

  typedef typename TransferOperatorType::row_type::ConstIterator TransferConstColumnIterator;

  // Clear coarse matrix
  if (this->recompute_ == NULL)
      coarseMat = 0;
  else {
      for (int i=0; i<coarseMat.N(); i++) {

          RowType& row = coarseMat[i];
          
          // Loop over all columns of the stiffness matrix
          ColumnIterator m    = row.begin();
          ColumnIterator mEnd = row.end();
          
          for (; m!=mEnd; ++m) {
              
              if ((*this->recompute_)[i][0] || (*this->recompute_)[m.index()][0])
                  *m = 0;

          }

      }

  }

  // Loop over all rows of the stiffness matrix
  for (int v=0; v<fineMat.N(); v++) {
      
      const RowType& row = fineMat[v];

      // Loop over all columns of the stiffness matrix
      ConstColumnIterator m    = row.begin();
      ConstColumnIterator mEnd = row.end();
        
      for (; m!=mEnd; ++m) {

          int w = m.index();

          // Loop over all coarse grid vectors iv that have v in their support
          TransferConstColumnIterator im    = (*this->matrix_)[v].begin();
          TransferConstColumnIterator imEnd = (*this->matrix_)[v].end();
          for (; im!=imEnd; ++im) {

              int iv = im.index();

              // Loop over all coarse grid vectors jv that have w in their support
              TransferConstColumnIterator jm    = (*this->matrix_)[w].begin();
              TransferConstColumnIterator jmEnd = (*this->matrix_)[w].end();
              
              for (; jm!=jmEnd; ++jm) {
                  
                  int jv = jm.index();

                  if (this->recompute_ && !((*this->recompute_)[iv][0]) && !((*this->recompute_)[jv][0]))
                      continue;
                  
                  typename OperatorType::block_type& cm = coarseMat[iv][jv];
                  
                  // Compute im * m * jm, but omitting the critical entries
                  for (int i=0; i<blocksize; i++) {
                      
                      for (int j=0; j<blocksize; j++) {
                          
                          // Truncated Multigrid:  Omit coupling if at least
                          // one of the two vectors is critical
                          if (!critical[v][i] && !critical[w][j])
                              cm[i][j] += (*im)[0][0] * (*m)[i][j] * (*jm)[0][0];
                          
                      }
                      
                  }

              }
          }
      }
  }

}