localmatrix.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_PDELAB_GRIDOPERATOR_COMMON_LOCALMATRIX_HH
#define DUNE_PDELAB_GRIDOPERATOR_COMMON_LOCALMATRIX_HH
 
#include <dune/common/iteratorfacades.hh>
 
#include <dune/pdelab/gridfunctionspace/localvector.hh>
 
namespace Dune {
  namespace PDELab {
 
    template<typename C>
    class WeightedMatrixAccumulationView
    {
 
    public:
 
      typedef C Container;
 
 
      typedef typename Container::BaseContainer BaseContainer;
 
      typedef typename C::value_type value_type;
 
      typedef typename C::weight_type weight_type;
 
      typedef WeightedMatrixAccumulationView WeightedAccumulationView;
 
      WeightedAccumulationView weightedAccumulationView(weight_type weight)
      {
        return WeightedAccumulationView(container(),weight*this->weight());
      }
 
      typedef typename C::size_type size_type;
 
 
      weight_type weight() const
      {
        return _weight;
      }
 
 
      void setWeight(weight_type weight)
      {
        _weight = weight;
      }
 
      template<typename LFSU, typename LFSV>
      void accumulate(const LFSV& lfsv, size_type i,
                      const LFSU& lfsu, size_type j,
                      value_type v)
      {
        _modified = true;
        _container(lfsv,i,lfsu,j) += _weight * v;
      }
 
 
      template<typename LFSU, typename LFSV>
      void rawAccumulate(const LFSV& lfsv, size_type i,
                         const LFSU& lfsu, size_type j,
                         value_type v)
      {
        _modified = true;
        _container(lfsv,i,lfsu,j) += v;
      }
 
      size_type nrows() const
      {
        return _container.nrows();
      }
 
      size_type ncols() const
      {
        return _container.ncols();
      }
 
      bool modified() const
      {
        return _modified;
      }
 
 
      void resetModified()
      {
        _modified = false;
      }
 
      Container& container()
      {
        _modified = true;
        return _container;
      }
 
      const Container& container() const
      {
        return _container;
      }
 
 
      BaseContainer& base()
      {
        _modified = true;
        return _container.base();
      }
 
 
      const BaseContainer& base() const
      {
        return _container.base();
      }
 
      WeightedMatrixAccumulationView(Container& container, weight_type weight)
        : _container(container)
        , _weight(weight)
        , _modified(false)
      {}
 
    private:
      Container& _container;
      weight_type _weight;
      bool _modified;
    };
 
 
        template<typename T, typename W = T>
     class LocalMatrix
        {
        public:
 
 
          typedef std::vector<T> BaseContainer;
 
          typedef typename BaseContainer::value_type  value_type;
 
          typedef typename BaseContainer::size_type    size_type;
 
          typedef typename BaseContainer::reference    reference;
 
          typedef typename BaseContainer::const_reference const_reference;
 
 
          typedef W weight_type;
 
          typedef WeightedMatrixAccumulationView<LocalMatrix> WeightedAccumulationView;
 
          struct iterator
            : public Dune::BidirectionalIteratorFacade<iterator,value_type>
          {
 
            iterator()
              : _m(nullptr)
              , _i(0)
              , _j(0)
            {}
 
            iterator(LocalMatrix& m, size_type i, size_type j)
              : _m(&m)
              , _i(i)
              , _j(j)
            {}
 
            bool equals(const iterator& other) const
            {
              return _m == other._m && _i == other._i && _j == other._j;
            }
 
            value_type& dereference() const
            {
              return _m->getEntry(_i,_j);
            }
 
            void increment()
            {
              if (_j < _m->ncols() - 1)
                ++_j;
              else
                {
                  ++_i;
                  _j = 0;
                }
            }
 
            void decrement()
            {
              if (_j > 0)
                --_j;
              else
                {
                  --_i;
                  _j = _m->ncols() - 1;
                }
            }
 
            size_type row() const
            {
              return _i;
            }
 
            size_type col() const
            {
              return _j;
            }
 
            LocalMatrix* _m;
            size_type _i;
            size_type _j;
 
          };
 
          iterator begin()
          {
            return iterator(*this,0,0);
          }
 
          iterator end()
          {
            return ncols() > 0 ? iterator(*this, nrows(), 0) : begin();
          }
 
          LocalMatrix () {}
 
          LocalMatrix (size_type r, size_type c)
            : _container(r*c)
            , _rows(r)
            , _cols(c)
          {}
 
          LocalMatrix (size_type r, size_type c, const T& t)
            : _container(r*c,t)
            , _rows(r)
            , _cols(c)
          {}
 
          void resize (size_type r, size_type c)
          {
            _container.resize(r*c);
            _rows = r;
            _cols = c;
          }
 
          LocalMatrix& operator=(const T& t)
          {
            std::fill(_container.begin(),_container.end(),t);
            return *this;
          }
 
          void assign (size_type r, size_type c, const T& t)
          {
            _container.assign(r*c,t);
            _rows = r;
            _cols = c;
          }
 
 
          template<typename LFSU, typename LFSV>
          T& operator() (const LFSV& lfsv, size_type i, const LFSU& lfsu, size_type j)
          {
            return getEntry(lfsv.localIndex(i),lfsu.localIndex(j));
          }
 
 
          template<typename LFSU, typename LFSV>
          const T& operator() (const LFSV& lfsv, size_type i, const LFSU& lfsu, size_type j) const
          {
            return getEntry(lfsv.localIndex(i),lfsu.localIndex(j));
          }
 
          LocalMatrix& operator *= (const T& x)
          {
            using namespace std::placeholders;
            std::transform(
              _container.begin(),
              _container.end(),
              _container.begin(),
              std::bind(std::multiplies<T>(),x,_1)
              );
            return *this;
          }
 
          size_type nrows () const
          {
            return _rows;
          }
 
          size_type ncols () const
          {
            return _cols;
          }
 
          template<class X, class R>
          void umv (const X& x, R& y) const
          {
            for (size_type i=0; i<_rows; ++i)
              {
                for (size_type j=0; j<_cols; j++)
                  accessBaseContainer(y)[i] += getEntry(i,j) * accessBaseContainer(x)[j];
              }
          }
 
          template<class X, class R>
          void usmv (const value_type& alpha, const X& x, R& y) const
          {
            for (size_type i=0; i<_rows; ++i)
              {
                for (size_type j=0; j<_cols; j++)
                  accessBaseContainer(y)[i] += alpha * getEntry(i,j) * accessBaseContainer(x)[j];
              }
          }
 
          WeightedAccumulationView weightedAccumulationView(weight_type weight)
          {
            return WeightedAccumulationView(*this,weight);
          }
 
 
          BaseContainer& base()
          {
            return _container;
          }
 
 
          const BaseContainer& base() const
          {
            return _container;
          }
 
 
          value_type& getEntry(size_type i, size_type j)
          {
            return _container[j*_rows + i];
          }
 
 
          const value_type& getEntry(size_type i, size_type j) const
          {
            return _container[j*_rows + i];
          }
 
        private:
 
          std::vector<T> _container;
          size_type _rows, _cols;
     };
 
    template<class Stream, class T, class W>
    Stream &operator<<(Stream &stream, const LocalMatrix<T,W> &m) {
      for(int r = 0; r < m.nrows(); ++r) {
        if(m.ncols() >= 1)
          stream << m.getEntry(r, 0);
        for(int c = 1; c < m.ncols(); ++c)
          stream << "\t" << m.getEntry(r, c);
        stream << "\n";
      }
      return stream;
    }
 
  } // namespace PDELab
} // namespace Dune
 
#endif // DUNE_PDELAB_GRIDOPERATOR_COMMON_LOCALMATRIX_HH