diff --git a/dune/solvers/transferoperators/mandelobsrestrictor.cc b/dune/solvers/transferoperators/mandelobsrestrictor.cc index 69aebce58dee0121f85348234802eeedc9991f1b..e63ed9353190d842e6e19a23465eb7bd06615b28 100644 --- a/dune/solvers/transferoperators/mandelobsrestrictor.cc +++ b/dune/solvers/transferoperators/mandelobsrestrictor.cc @@ -59,3 +59,63 @@ restrict(const std::vector<BoxConstraint<typename DiscFuncType::field_type,block } } + +template <class DiscFuncType> +void MandelObstacleRestrictor<DiscFuncType>:: +restrict(const std::vector<BoxConstraint<typename DiscFuncType::field_type,blocksize> >&f, + std::vector<BoxConstraint<typename DiscFuncType::field_type,blocksize> > &t, + const Dune::BitSetVector<blocksize>& fHasObstacle, + const Dune::BitSetVector<blocksize>& tHasObstacle, + const MultigridTransfer<DiscFuncType>& transfer, + const Dune::BitSetVector<blocksize>& critical) +{ + if (!(dynamic_cast<const CompressedMultigridTransfer<DiscFuncType>*>(&transfer))) + DUNE_THROW(Dune::NotImplemented, "Implementation assumes transfer to be of type CompressedMultigridTransfer!"); + + const Dune::BCRSMatrix<Dune::FieldMatrix<field_type, 1, 1> >& transferMat = + dynamic_cast<const CompressedMultigridTransfer<DiscFuncType>*>(&transfer)->getMatrix(); + + t.resize(transferMat.M()); + for (size_t i=0; i<t.size(); i++) + t[i].clear(); + + // The Mandel restriction computes the maximum $m$ of all nodal base coefficients + // in the support of a given coarse grid base function \f$ \phi_c \f$. If the + // coefficient of \f$ \phi_c \f$ is smaller than \f$ m \f$ , it is set to \f$ m \f$. + + // We use the multigrid transfer matrix to find all relevant fine grid dofs + typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<field_type,1,1> >::row_type RowType; + typedef typename RowType::ConstIterator ColumnIterator; + + // Loop over all coarse grid dofs + for (size_t i=0; i<transferMat.N(); i++) { + + if (fHasObstacle[i].none()) + continue; + + const RowType& row = transferMat[i]; + + ColumnIterator cIt = row.begin(); + ColumnIterator cEndIt = row.end(); + + // Loop over all fine grid base functions in the support of tObsIt + for(; cIt!=cEndIt; ++cIt) { + + if (tHasObstacle[cIt.index()].none()) + continue; + + // Sort out spurious entries due to numerical dirt + if ((*cIt)[0][0] < 1e-3) + continue; + + for (int j=0; j<blocksize; j++) + if (!critical[i][j] && fHasObstacle[i][j] && tHasObstacle[cIt.index()][j]) { + t[cIt.index()].lower(j) = std::max(t[cIt.index()].lower(j), f[i].lower(j)); + t[cIt.index()].upper(j) = std::min(t[cIt.index()].upper(j), f[i].upper(j)); + } + + } + + } + +} diff --git a/dune/solvers/transferoperators/mandelobsrestrictor.hh b/dune/solvers/transferoperators/mandelobsrestrictor.hh index c22621da09d614bced389e8562836b6033d711fa..4ab2d0558faa72c5d1a73c96144f857ae18e93d2 100644 --- a/dune/solvers/transferoperators/mandelobsrestrictor.hh +++ b/dune/solvers/transferoperators/mandelobsrestrictor.hh @@ -21,6 +21,13 @@ public: const Dune::BitSetVector<1>& tHasObstacle, const MultigridTransfer<DiscFuncType>& transfer, const Dune::BitSetVector<blocksize>& critical); + + virtual void restrict(const std::vector<BoxConstraint<field_type,blocksize> >& f, + std::vector<BoxConstraint<field_type,blocksize> >& t, + const Dune::BitSetVector<blocksize>& fHasObstacle, + const Dune::BitSetVector<blocksize>& tHasObstacle, + const MultigridTransfer<DiscFuncType>& transfer, + const Dune::BitSetVector<blocksize>& critical); }; // include implementation