subdomain.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUNCTIONS_COMMON_SUBDOMAIN_HH
#define DUNE_FUNCTIONS_COMMON_SUBDOMAIN_HH
 
#include <array>
#include <cstddef>
#include <limits>
#include <optional>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>
 
#include <dune/common/exceptions.hh>
#include <dune/common/iteratorfacades.hh>
#include <dune/common/iteratorrange.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/std/no_unique_address.hh>
 
#include <dune/geometry/type.hh>
#include <dune/geometry/typeindex.hh>
 
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/partitionset.hh>
 
namespace Dune::Functions::Experimental {
 
  namespace Impl {
 
    template<class GV, class ContainsCallback>
    class GlobalIntersectionIt
      : public Dune::IteratorFacade<GlobalIntersectionIt<GV, ContainsCallback>, std::forward_iterator_tag, const typename GV::Intersection>
    {
      using Facade = Dune::IteratorFacade<GlobalIntersectionIt<GV, ContainsCallback>, std::forward_iterator_tag, const typename GV::Intersection>;
 
    public:
 
      using GridView = GV;
      using Element = typename GridView::template Codim<0>::Entity;
      using ElementIterator = typename GridView::template Codim<0>::Iterator;
      using IntersectionIterator = typename GridView::IntersectionIterator;
      using Intersection = typename GridView::Intersection;
 
    private:
      static constexpr bool cacheIntersection = not std::is_lvalue_reference_v<decltype(std::declval<IntersectionIterator>().operator*())>;
      static constexpr bool cacheElement = not std::is_lvalue_reference_v<decltype(std::declval<ElementIterator>().operator*())>;
 
      using ElementStorage = std::conditional_t<cacheElement, std::optional<Element>, std::monostate>;
      using IntersectionStorage = std::conditional_t<cacheIntersection, std::optional<Intersection>, std::monostate>;
 
      const Element& element() const
      {
        if constexpr (cacheElement)
          return *element_;
        else
          return *elementIt_;
      }
 
      void updateElement()
      {
        if constexpr (cacheElement)
          element_ = *elementIt_;
      }
 
      const Intersection& intersection() const
      {
        if constexpr (cacheIntersection)
          return *intersection_;
        else
          return **iIt_;
      }
 
      void updateIntersection()
      {
        if constexpr (cacheIntersection)
          intersection_ = **iIt_;
      }
 
    public:
 
      class SentinelIterator
      {};
 
      GlobalIntersectionIt(const GridView& gridView, const ContainsCallback& contains, ElementIterator elementIt, ElementIterator elementEnd)
        : gridView_(gridView)
        , contains_(contains)
        , elementIt_(std::move(elementIt))
        , elementEnd_(std::move(elementEnd))
      {
        if (elementIt_ != elementEnd_)
        {
          updateElement();
          iIt_ = gridView_.ibegin(element());
          updateIntersection();
          if (not contains_(intersection()))
            ++(*this);
        }
      }
 
      GlobalIntersectionIt(const GridView& gridView, const ContainsCallback& contains, ElementIterator elementIt)
        : GlobalIntersectionIt(gridView, contains, elementIt, gridView.template end<0>())
      {}
 
      GlobalIntersectionIt(const GridView& gridView, const ContainsCallback& contains)
        : GlobalIntersectionIt(gridView, contains, gridView.template begin<0>())
      {}
 
      using reference = typename Facade::reference;
 
      reference operator*() const
      {
        return intersection();
      }
 
      GlobalIntersectionIt& operator++()
      {
        while(true)
        {
          ++(*iIt_);
          if (*iIt_ == gridView_.iend(element()))
          {
            ++elementIt_;
            if (elementIt_ == elementEnd_)
              return *this;
            updateElement();
            iIt_ = gridView_.ibegin(element());
          }
          updateIntersection();
          if (contains_(intersection()))
            return *this;
        }
        return *this;
      }
 
      friend bool operator==(const GlobalIntersectionIt& it1, const GlobalIntersectionIt& it2)
      {
        if (it1.elementIt_ != it2.elementIt_)
          return false;
        if (it1.elementIt_ == it1.elementEnd_)
          return true;
        return (*(it1.iIt_) == *(it2.iIt_));
      }
 
      friend bool operator==(const GlobalIntersectionIt& it1, const SentinelIterator& it2)
      {
        return it1.elementIt_ == it1.elementEnd_;
      }
 
    private:
      GridView gridView_;
      ContainsCallback contains_;
      ElementIterator elementIt_;
      ElementIterator elementEnd_;
      std::optional<IntersectionIterator> iIt_;
      DUNE_NO_UNIQUE_ADDRESS ElementStorage element_;
      DUNE_NO_UNIQUE_ADDRESS IntersectionStorage intersection_;
    };
 
  }
 
 
 
  template<class HGV>
  class SubDomainIndexSet
  {
    using HostGridView = HGV;
 
  public:
 
    using Grid = typename HostGridView::Grid;
    using Types = std::vector<Dune::GeometryType>;
    using IndexType = std::size_t;
 
    template<int codim>
    struct Codim
    {
      using Entity = typename Grid::template Codim<codim>::Entity;
      using EntitySeed = typename Grid::template Codim<codim>::EntitySeed;
      using Geometry = typename Grid::template Codim<codim>::Geometry;
      using LocalGeometry = typename Grid::template Codim<codim>::LocalGeometry;
    };
 
    enum {dimension = Grid::dimension};
 
  private:
 
    using AllEntityMapper = Dune::MultipleCodimMultipleGeomTypeMapper<HostGridView>;
 
    static auto allCodimLayout()
    {
      return [](Dune::GeometryType, int) { return true; };
    }
 
    static constexpr auto typeIndexSize = Dune::GlobalGeometryTypeIndex::size(dimension);
    static constexpr auto unusesIndex = std::numeric_limits<IndexType>::max();
 
  public:
 
    SubDomainIndexSet(const HostGridView& hostGridView)
      : hostGridView_(hostGridView)
      , allEntityMapper_(hostGridView_, allCodimLayout())
    {
      clear();
    }
 
    // *********************************
    // IndexSet interface methods
    // *********************************
 
    IndexType size(Dune::GeometryType gt) const
    {
      if (gt.dim() <= dimension)
        return sizePerGT_[Dune::GlobalGeometryTypeIndex::index(gt)];
      else
        return 0;
    }
 
    IndexType size(int codim) const
    {
      return sizePerCodim_[codim];
    }
 
    template<class Entity>
    IndexType index(const Entity& entity) const
    {
      auto index = indices_[allEntityMapper_.index(entity)];
      if (index==unusesIndex)
        DUNE_THROW(Dune::InvalidStateException, "Accessing nonexisting entry using SubDomainIndexSet::index()!");
      return index;
    }
 
    template<int cc>
    IndexType index(const typename Codim<cc>::Entity& entity) const
    {
      return index<typename Codim<cc>::Entity>(entity);
    }
 
    template<class Entity>
    IndexType subIndex(const Entity& entity, int subEntity, unsigned int codim) const
    {
      auto index = indices_[allEntityMapper_.subIndex(entity, subEntity, codim)];
      if (index==unusesIndex)
        DUNE_THROW(Dune::InvalidStateException, "Accessing nonexisting entry using SubDomainIndexSet::subIndex()!");
      return index;
    }
 
    template<int cc>
    IndexType subIndex(const typename Codim<cc>::Entity& entity, int subEntity, unsigned int codim) const
    {
      return subIndex<typename Codim<cc>::Entity>(entity, subEntity, codim);
    }
 
    template<class Entity >
    bool contains(const Entity& entity) const
    {
      return (indices_[allEntityMapper_.index(entity)] != unusesIndex);
    }
 
    Types types(int codim) const
    {
      return typesPerCodim_[codim];
    }
 
    // *********************************
    // Extended methods
    // *********************************
 
    const HostGridView& hostGridView() const
    {
      return hostGridView_;
    }
 
    void insertElement(const typename Codim<0>::Entity& element)
    {
      const auto& re = referenceElement(element);
      for (auto codim : Dune::range(0, dimension+1))
      {
        for (auto subEntity : Dune::range(re.size(codim)))
        {
          auto& index = indices_[allEntityMapper_.subIndex(element, subEntity, codim)];
          if (index==unusesIndex)
          {
            const auto& type = re.type(subEntity, codim);
            const auto typeIndex = Dune::GlobalGeometryTypeIndex::index(type);
            index = sizePerGT_[typeIndex];
            if (sizePerGT_[typeIndex]==0)
              typesPerCodim_[codim].push_back(type);
            sizePerGT_[typeIndex]++;
            sizePerCodim_[codim]++;
          }
        }
      }
    }
 
  protected:
 
    // Clear all data
    void clear()
    {
      for(auto& size : sizePerGT_)
        size = 0;
      for(auto& size : sizePerCodim_)
        size = 0;
      for(auto& types : typesPerCodim_)
        types.clear();
      indices_.clear();
      indices_.resize(allEntityMapper_.size(), unusesIndex);
    }
 
    HostGridView hostGridView_;
 
    // Global size information
    std::array<std::size_t, typeIndexSize> sizePerGT_;
    std::array<std::size_t, dimension+1> sizePerCodim_;
    std::array<Types, dimension+1> typesPerCodim_;
 
    AllEntityMapper allEntityMapper_;
 
    // Index map
    std::vector<IndexType> indices_;
  };
 
 
 
  template<class HGV>
  class SubDomainGridView
  {
 
    template<int codim>
    class NonImplementedIterator
    {
    public:
      NonImplementedIterator()
      {
        static_assert(codim==0, "SubDomainGridView::Codim::Iterator<codim> is only implemented for codim=0");
      }
    };
 
    template<PartitionIteratorType pit>
    class ElementIterator
    {
      using Element = typename HGV::template Codim<0>::Entity;
    public:
 
      using HostElementIterator = typename HGV::template Codim<0>::template Partition<pit>::Iterator;
 
      ElementIterator(const SubDomainIndexSet<HGV>& indexSet, HostElementIterator&& it, HostElementIterator&& endIt)
        : indexSet_(&indexSet)
        , hostIt_(std::move(it))
        , hostEndIt_(std::move(endIt))
      {
        while ((hostIt_!= hostEndIt_) and (not indexSet_->contains(*hostIt_)))
          ++hostIt_;
      }
 
      ElementIterator& operator++()
      {
        ++hostIt_;
        while ((hostIt_!= hostEndIt_) and (not indexSet_->contains(*hostIt_)))
          ++hostIt_;
        return *this;
      }
 
      const Element& operator*() const
      {
        return *hostIt_;
      }
 
      friend bool operator==(const ElementIterator& a, const ElementIterator& b)
      {
        return a.hostIt_==b.hostIt_;
      }
 
    private:
      HostElementIterator hostIt_;
      HostElementIterator hostEndIt_;
      const SubDomainIndexSet<HGV>* indexSet_;
    };
 
  public:
 
    using HostGridView = HGV;
 
    using Grid = typename HostGridView::Grid;
    using ctype = typename Grid::ctype;
    using IndexSet = SubDomainIndexSet<HostGridView>;
    using Intersection = typename HostGridView::Intersection;
    using IntersectionIterator = typename HostGridView::IntersectionIterator;
 
    template<int codim>
    struct Codim
    {
      using Entity = typename Grid::template Codim<codim>::Entity;
      using EntitySeed = typename Grid::template Codim<codim>::EntitySeed;
      using Geometry = typename Grid::template Codim<codim>::Geometry;
      using LocalGeometry = typename Grid::template Codim<codim>::LocalGeometry;
      using Iterator = std::conditional_t<codim==0, ElementIterator<All_Partition>, NonImplementedIterator<codim>>;
 
      template<PartitionIteratorType pit>
      struct Partition
      {
        using Iterator = std::conditional_t<codim==0, ElementIterator<pit>, NonImplementedIterator<codim>>;
      };
    };
 
    enum {dimension = Grid::dimension};
    enum {dimensionworld = Grid::dimensionworld};
 
    SubDomainGridView(const IndexSet& indexSet)
      : indexSet_(&indexSet)
    {}
 
    SubDomainGridView(const SubDomainGridView& other) = default;
 
    const Grid& grid() const
    {
      return indexSet_->hostGridView().grid();
    }
 
    const IndexSet& indexSet() const
    {
      return *indexSet_;
    }
 
    int size(int codim) const
    {
      return indexSet_->size(codim);
    }
 
    int size(Dune::GeometryType gt) const
    {
      return indexSet_->size(gt);
    }
 
    template<class Entity>
    bool contains(const Entity& entity) const
    {
      return indexSet_->contains(entity);
    }
 
    template<int codim, PartitionIteratorType pit = All_Partition>
    typename Codim<codim>::template Partition<pit>::Iterator begin() const
    {
      static_assert(codim==0, "SubDomainGridView::begin<codim> is only implemented for codim=0");
      return {indexSet(), hostGridView().template begin<codim, pit>(), hostGridView().template end<codim, pit>()};
    }
 
    template<int codim, PartitionIteratorType pit = All_Partition>
    typename Codim<codim>::template Partition<pit>::Iterator end() const
    {
      static_assert(codim==0, "SubDomainGridView::end<codim> is only implemented for codim=0");
      return {indexSet(), hostGridView().template end<codim, pit>(), hostGridView().template end<codim, pit>()};
    }
 
    decltype(auto) comm() const
    {
      return hostGridView().comm();
    }
 
    decltype(auto) ibegin(const typename Codim<0>::Entity& element) const
    {
      return hostGridView().ibegin(element);
    }
 
    decltype(auto) iend(const typename Codim<0>::Entity& element) const
    {
      return hostGridView().iend(element);
    }
 
    const HostGridView& hostGridView() const
    {
      return indexSet_->hostGridView();
    }
 
  protected:
    const IndexSet* indexSet_;
  };
 
 
 
  template<class HostGridView>
  auto elements(const SubDomainGridView<HostGridView>& subDomainGridView)
  {
    return Dune::IteratorRange(subDomainGridView.template begin<0>(), subDomainGridView.template end<0>());
  }
 
  template<class HostGridView, unsigned int partitions>
  auto elements(const SubDomainGridView<HostGridView>& subDomainGridView, Dune::PartitionSet<partitions> partitionSet)
  {
    constexpr auto pit = partitionSet.partitionIterator();
    return Dune::IteratorRange(subDomainGridView.template begin<0, pit>(), subDomainGridView.template end<0, pit>());
  }
 
  template<class HostGridView, class Element>
  auto intersections(const SubDomainGridView<HostGridView>& subDomainGridView, const Element& element)
  {
    return Dune::IteratorRange(subDomainGridView.ibegin(element), subDomainGridView.iend(element));
  }
 
 
 
  template<class HGV>
  class SubDomain
  {
  public:
 
    using HostGridView = HGV;
    using Grid = typename HostGridView::Grid;
    using IndexSet = SubDomainIndexSet<HostGridView>;
    using GridView = SubDomainGridView<HostGridView>;
 
    template<int codim>
    struct Codim
    {
      using Entity = typename Grid::template Codim<codim>::Entity;
      using EntitySeed = typename Grid::template Codim<codim>::EntitySeed;
      using Geometry = typename Grid::template Codim<codim>::Geometry;
      using LocalGeometry = typename Grid::template Codim<codim>::LocalGeometry;
    };
 
    enum {dimension = Grid::dimension};
 
    SubDomain(const HostGridView& hostGridView)
      : indexSet_(hostGridView)
    {}
 
    const IndexSet& indexSet() const
    {
      return indexSet_;
    }
 
    GridView gridView() const
    {
      return GridView(indexSet_);
    }
 
    HostGridView hostGridView() const
    {
      return indexSet_.hostGridView();
    }
 
    void insertElement(const typename Codim<0>::Entity& element)
    {
      indexSet_.insertElement(element);
    }
 
    bool contains(const typename Codim<0>::Entity& element) const
    {
      return indexSet_.contains(element);
    }
 
  private:
    IndexSet indexSet_;
  };
 
 
 
  template<class SubDomainA, class SubDomainB>
  class SubDomainInterface
  {
    static_assert(
      std::is_same_v<typename SubDomainA::GridView::Intersection, typename SubDomainB::GridView::Intersection>,
      "SubDomainInterface requires that both SubDomain types have the same Intersection type");
 
  public:
 
    using Intersection = typename SubDomainA::GridView::Intersection;
 
    SubDomainInterface(const SubDomainA& subDomainA, const SubDomainB& subDomainB)
      : subDomainA_(subDomainA)
      , subDomainB_(subDomainB)
    {}
 
    bool contains(const Intersection& is) const
    {
      if (is.boundary() or not(is.neighbor()))
        return false;
      return (subDomainA_.contains(is.inside()) && subDomainB_.contains(is.outside()))
          || (subDomainA_.contains(is.outside()) && subDomainB_.contains(is.inside()));
    }
 
    bool isOriented(const Intersection& is) const
    {
      if (is.boundary() or not(is.neighbor()))
        return false;
      return (subDomainA_.contains(is.inside()) && subDomainB_.contains(is.outside()));
    }
 
    const auto begin() const
    {
      return Impl::GlobalIntersectionIt(subDomainA_.gridView(), [&](const auto& is) {
        if (is.boundary() or not(is.neighbor()))
          return false;
        return subDomainB_.indexSet().contains(is.outside());
      });
    }
 
    const auto end() const
    {
      return typename decltype(begin())::SentinelIterator();
    }
 
  private:
    const SubDomainA& subDomainA_;
    const SubDomainB& subDomainB_;
  };
 
 
 
  template<class SubDomain>
  class SubDomainSkeleton
  {
  public:
 
    using Intersection = typename SubDomain::GridView::Intersection;
 
    SubDomainSkeleton(const SubDomain& subDomain)
      : subDomain_(subDomain)
    {}
 
    bool contains(const Intersection& is) const
    {
      if (is.boundary() or not(is.neighbor()))
        return false;
      return subDomain_.contains(is.inside()) and subDomain_.contains(is.outside());
    }
 
  private:
    const SubDomain& subDomain_;
  };
 
 
 
} // namespace Dune::Functions::Experimental
 
#endif// DUNE_FUNCTIONS_COMMON_SUBDOMAIN_HH