hierarchy.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_AMGHIERARCHY_HH
#define DUNE_AMGHIERARCHY_HH

#include <list>
#include <memory>
#include <limits>
#include <dune/common/stdstreams.hh>
#include <dune/common/timer.hh>
#include <dune/common/bigunsignedint.hh>
#include <dune/istl/paamg/construction.hh>

namespace Dune
{
  namespace Amg
  {
    template<typename T, typename A=std::allocator<T> >
    class Hierarchy
    {
    public:
      typedef T MemberType;

      template<typename T1, typename T2>
      class LevelIterator;

    private:
      struct Element
      {
        friend class LevelIterator<Hierarchy<T,A>, T>;
        friend class LevelIterator<const Hierarchy<T,A>, const T>;

        std::weak_ptr<Element> coarser_;

        std::shared_ptr<Element> finer_;

        std::shared_ptr<MemberType> element_;

        std::shared_ptr<MemberType> redistributed_;
      };
    public:

      using Allocator = typename std::allocator_traits<A>::template rebind_alloc<Element>;

      typedef typename ConstructionTraits<T>::Arguments Arguments;

      Hierarchy(const std::shared_ptr<MemberType> & first);

      Hierarchy() : levels_(0)
      {}

      Hierarchy(const Hierarchy& other);

      void addCoarser(Arguments& args);

      void addRedistributedOnCoarsest(Arguments& args);

      void addFiner(Arguments& args);

      template<class C, class T1>
      class LevelIterator
        : public BidirectionalIteratorFacade<LevelIterator<C,T1>,T1,T1&>
      {
        friend class LevelIterator<typename std::remove_const<C>::type,
            typename std::remove_const<T1>::type >;
        friend class LevelIterator<const typename std::remove_const<C>::type,
            const typename std::remove_const<T1>::type >;

      public:
        LevelIterator() = default;

        LevelIterator(std::shared_ptr<Element> element)
          : element_(std::move(element))
        {}

        template <class C_, class T1_,
          std::enable_if_t<std::is_same_v<std::remove_const_t<C>, std::remove_const_t<C_>>, int> = 0,
          std::enable_if_t<std::is_same_v<std::remove_const_t<T1>, std::remove_const_t<T1_>>, int> = 0>
        LevelIterator(const LevelIterator<C_,T1_>& other)
          : element_(other.element_)
        {}

        template <class C_, class T1_,
          std::enable_if_t<std::is_same_v<std::remove_const_t<C>, std::remove_const_t<C_>>, int> = 0,
          std::enable_if_t<std::is_same_v<std::remove_const_t<T1>, std::remove_const_t<T1_>>, int> = 0>
        LevelIterator& operator=(const LevelIterator<C_,T1_>& other)
        {
          element_ = other.element_;
          return *this;
        }

        bool equals(const LevelIterator<typename std::remove_const<C>::type,
                        typename std::remove_const<T1>::type>& other) const
        {
          return element_ == other.element_;
        }

        bool equals(const LevelIterator<const typename std::remove_const<C>::type,
                        const typename std::remove_const<T1>::type>& other) const
        {
          return element_ == other.element_;
        }

        T1& dereference() const
        {
          return *(element_->element_);
        }

        void increment()
        {
          element_ = element_->coarser_.lock();
        }

        void decrement()
        {
          element_ = element_->finer_;
        }

        bool isRedistributed() const
        {
          return (bool)element_->redistributed_;
        }

        T1& getRedistributed() const
        {
          assert(element_->redistributed_);
          return *element_->redistributed_;
        }
        void addRedistributed(std::shared_ptr<T1> t)
        {
          element_->redistributed_ = t;
        }

        void deleteRedistributed()
        {
          element_->redistributed_ = nullptr;
        }

      private:
        std::shared_ptr<Element> element_ = {};
      };

      typedef LevelIterator<Hierarchy<T,A>,T> Iterator;

      typedef LevelIterator<const Hierarchy<T,A>, const T> ConstIterator;

      Iterator finest();

      Iterator coarsest();


      ConstIterator finest() const;

      ConstIterator coarsest() const;

      std::size_t levels() const;

    private:
      std::shared_ptr<MemberType> originalFinest_;
      std::shared_ptr<Element> finest_;
      std::shared_ptr<Element> coarsest_;
      Allocator allocator_;
      int levels_;
    };

    template<class T, class A>
    Hierarchy<T,A>::Hierarchy(const std::shared_ptr<MemberType> & first)
      : originalFinest_(first)
    {
      finest_ = std::allocate_shared<Element>(allocator_);
      finest_->element_ = originalFinest_;
      coarsest_ = finest_;
      levels_ = 1;
    }

    //TODO: do we actually want to support this? This might be very expensive?!
    template<class T, class A>
    Hierarchy<T,A>::Hierarchy(const Hierarchy& other)
    : allocator_(other.allocator_),
      levels_(other.levels_)
    {
      if(!other.finest_)
      {
        finest_=coarsest_=nullptr;
        return;
      }
      finest_ = std::allocate_shared<Element>(allocator_);
      std::shared_ptr<Element> finer_;
      std::shared_ptr<Element> current_ = finest_;
      std::weak_ptr<Element> otherWeak_ = other.finest_;

      while(! otherWeak_.expired())
      {
        // create shared_ptr from weak_ptr, we just checked that this is safe
        std::shared_ptr<Element> otherCurrent_ = std::shared_ptr<Element>(otherWeak_);
        // clone current level
        //TODO: should we use the allocator?
        current_->element_ =
          std::make_shared<MemberType>(*(otherCurrent_->element_));
        current_->finer_=finer_;
        if(otherCurrent_->redistributed_)
          current_->redistributed_ =
            std::make_shared<MemberType>(*(otherCurrent_->redistributed_));
        finer_=current_;
        if(not otherCurrent_->coarser_.expired())
        {
          auto c = std::allocate_shared<Element>(allocator_);
          current_->coarser_ = c;
          current_ = c;
        }
        // go to coarser level
        otherWeak_ = otherCurrent_->coarser_;
      }
      coarsest_=current_;
    }

    template<class T, class A>
    std::size_t Hierarchy<T,A>::levels() const
    {
      return levels_;
    }

    template<class T, class A>
    void Hierarchy<T,A>::addRedistributedOnCoarsest(Arguments& args)
    {
      coarsest_->redistributed_ = ConstructionTraits<MemberType>::construct(args);
    }

    template<class T, class A>
    void Hierarchy<T,A>::addCoarser(Arguments& args)
    {
      if(!coarsest_) {
        // we have no levels at all...
        assert(!finest_);
        // allocate into the shared_ptr
        originalFinest_ = ConstructionTraits<MemberType>::construct(args);
        coarsest_ = std::allocate_shared<Element>(allocator_);
        coarsest_->element_ = originalFinest_;
        finest_ = coarsest_;
      }else{
        auto old_coarsest = coarsest_;
        coarsest_ = std::allocate_shared<Element>(allocator_);
        coarsest_->finer_ = old_coarsest;
        coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
        old_coarsest->coarser_ = coarsest_;
      }
      ++levels_;
    }


    template<class T, class A>
    void Hierarchy<T,A>::addFiner(Arguments& args)
    {
      //TODO: wouldn't it be better to do this in the constructor?'
      if(!finest_) {
        // we have no levels at all...
        assert(!coarsest_);
        // allocate into the shared_ptr
        originalFinest_ = ConstructionTraits<MemberType>::construct(args);
        finest_ = std::allocate_shared<Element>(allocator_);
        finest_->element = originalFinest_;
        coarsest_ = finest_;
      }else{
        finest_->finer_ = std::allocate_shared<Element>(allocator_);
        finest_->finer_->coarser_ = finest_;
        finest_ = finest_->finer_;
        finest_->element = ConstructionTraits<T>::construct(args);
      }
      ++levels_;
    }

    template<class T, class A>
    typename Hierarchy<T,A>::Iterator Hierarchy<T,A>::finest()
    {
      return Iterator(finest_);
    }

    template<class T, class A>
    typename Hierarchy<T,A>::Iterator Hierarchy<T,A>::coarsest()
    {
      return Iterator(coarsest_);
    }

    template<class T, class A>
    typename Hierarchy<T,A>::ConstIterator Hierarchy<T,A>::finest() const
    {
      return ConstIterator(finest_);
    }

    template<class T, class A>
    typename Hierarchy<T,A>::ConstIterator Hierarchy<T,A>::coarsest() const
    {
      return ConstIterator(coarsest_);
    }
  } // namespace Amg
} // namespace Dune

#endif