truncatedcompressedmgtransfer.cc

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:

template<class VectorType, class BitVectorType, class MatrixType>
void TruncatedCompressedMGTransfer<VectorType, BitVectorType, MatrixType>::prolong(const VectorType& f, VectorType& t) const
{
    // If the critical bitfield is not set, call the base class method
    if (this->critical_ == nullptr) {
        Base::prolong(f,t);
        return;
    }

    const BitVectorType& critical = *this->critical_;

    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 VectorType::Iterator      Iterator;

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

    Iterator tIt      = t.begin();

    for(size_t 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 VectorType, class BitVectorType, class MatrixType>
void TruncatedCompressedMGTransfer<VectorType, BitVectorType, MatrixType>::restrict(const VectorType& f, VectorType& t) const
{
    // If the critical bitfield is not set, call the base class method
    if (this->critical_ == nullptr) {
        Base::restrict(f,t);
        return;
    }

    const BitVectorType& critical = *this->critical_;

    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 TransferOperatorType::row_type RowType;
    typedef typename RowType::ConstIterator ColumnIterator;

    for (size_t 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 VectorType::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 VectorType, class BitVectorType, class MatrixType>
void TruncatedCompressedMGTransfer<VectorType, BitVectorType, MatrixType>::
galerkinRestrict(const MatrixType& fineMat, MatrixType& coarseMat) const
{
    // If the critical bitfield is not set, call the base class method
    if (this->recompute_ == nullptr && this->critical_ == nullptr) {
        Base::galerkinRestrict(fineMat, coarseMat);
        return;
    }

    if (this->recompute_ != nullptr && 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 MatrixType::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_ == nullptr)
      coarseMat = 0;
  else {
      for (size_t 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].any() || (*this->recompute_)[m.index()].any())
                  *m = 0;
          }
      }
  }

  // Loop over all rows of the stiffness matrix
  for (size_t 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].none() && (*this->recompute_)[jv].none())
                      continue;

                  typename MatrixType::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 (this->critical_==nullptr || (!((*this->critical_)[v][i]) && !((*this->critical_)[w][j])))
                              cm[i][j] += (*im)[0][0] * (*m)[i][j] * (*jm)[0][0];
                      }
                  }
              }
          }
      }
  }
}