BasisMapping.hpp

#pragma once
 
#include <cassert>
#include <map>
#include <vector>
 
#include <amdis/Observer.hpp>
#include <amdis/typetree/MultiIndex.hpp>
#include <dune/common/shared_ptr.hh>
 
namespace AMDiS
{
 
  template <class SourceBasis, class TargetBasis,
            class Storage = std::map<typename SourceBasis::MultiIndex, typename TargetBasis::MultiIndex>>
  class BasisMapping
      : private ObserverSequence<event::adapt,2>
  {
  public:
    using SourceBasisType = SourceBasis;
    using TargetBasisType = TargetBasis;
 
    BasisMapping(std::shared_ptr<SourceBasis> const& source, std::shared_ptr<TargetBasis> const& target)
      : ObserverSequence<event::adapt,2>(*source, *target)
      , source_(source)
      , target_(target)
    {
      // require the grids to be the same
      assert(&source_->gridView().grid() == &target_->gridView().grid());
 
      // calculate an initial mapping
      updateImpl2();
    }
 
    BasisMapping(SourceBasis& source, TargetBasis& target)
      : BasisMapping(Dune::stackobject_to_shared_ptr(source), Dune::stackobject_to_shared_ptr(target))
    {}
 
  protected:
 
    void updateImpl(event::adapt e, index_t<0> i) override { if (update_) updateImpl2(); else update_ = true; }
    void updateImpl(event::adapt e, index_t<1> i) override { if (update_) updateImpl2(); else update_ = true; }
 
  private:
 
    // Construct a mapping from the indices of the source basis to the indices of the target basis
    void updateImpl2()
    {
      mapping_.clear();
      resize(mapping_);
 
      auto sourceLocalView = source_->localView();
      auto targetLocalView = target_->localView();
      for (auto const& e : elements(source_->gridView())) {
        sourceLocalView.bind(e);
        targetLocalView.bind(e);
 
        assert(sourceLocalView.size() == targetLocalView.size());
        for (std::size_t i = 0, s = sourceLocalView.size(); i < s; ++i)
          mapping_[sourceLocalView.index(i)] = targetLocalView.index(i);
      }
 
      update_ = false;
    }
 
  public:
 
    typename TargetBasis::MultiIndex const& operator()(typename SourceBasis::MultiIndex const& mi) const
    {
      return mapping_.at(mi);
    }
 
    typename TargetBasis::MultiIndex const& operator()(std::size_t idx) const
    {
      static_assert(Traits::IsFlatIndex<typename SourceBasis::MultiIndex>::value, "");
      typename SourceBasis::MultiIndex mi{idx};
      return mapping_.at(mi);
    }
 
    TargetBasis const& sizeInfo() const
    {
      return *target_;
    }
 
 
    std::shared_ptr<SourceBasis> const& source() const
    {
      return source_;
    }
 
    std::shared_ptr<TargetBasis> const& target() const
    {
      return target_;
    }
 
  private:
    // vector must be resized before assigning indices
    template <class T, class A>
    void resize(std::vector<T,A>& mapping) const
    {
      mapping.resize(source_->dimension());
    }
 
    // map or unordered_map do not need a resize
    template <class Container>
    void resize(Container& mapping) const
    {
      /* do nothing */
    }
 
  private:
    std::shared_ptr<SourceBasis> source_;
    std::shared_ptr<TargetBasis> target_;
    Storage mapping_;
 
    bool update_ = false;
  };
 
} // end namespace AMDiS