BoundaryManager.hpp

#pragma once

#include <memory>
#include <vector>

#include <dune/common/hybridutilities.hh>
#include <dune/common/std/type_traits.hh>

#include <amdis/Boundary.hpp>
#include <amdis/common/Concepts.hpp>

namespace AMDiS
{
  class BoundaryManagerBase
  {
  public:
    BoundaryManagerBase(std::size_t numBoundarySegments)
      : boundaryIds_(numBoundarySegments, BoundaryType{1})
    {}

    template <class Intersection>
    BoundaryType boundaryId(Intersection const& intersection) const
    {
      return boundaryIds_[intersection.boundarySegmentIndex()];
    }

    std::vector<BoundaryType> const& boundaryIds() const
    {
      return boundaryIds_;
    }

  protected:
    std::vector<BoundaryType> boundaryIds_; // maps a boundarySegmentIndex to an ID
  };



  template <class G>
  class BoundaryManager
      : public BoundaryManagerBase
  {
    using Super = BoundaryManagerBase;

  public:
    using Grid = G;
    enum { dim = Grid::dimension };
    enum { dow = Grid::dimensionworld };

    using Segment = typename Grid::LevelGridView::Intersection;
    using Domain = typename Grid::template Codim<0>::Geometry::GlobalCoordinate;

  public:
    BoundaryManager(std::shared_ptr<G> const& grid)
      : Super(grid->numBoundarySegments())
      , grid_(grid)
    {
      setBoundaryId();
    }

    BoundaryManager(G& grid)
      : BoundaryManager{Dune::stackobject_to_shared_ptr(grid)}
    {}


    void setBoxBoundary(std::array<BoundaryType, 2*dow> const& ids)
    {
      auto gv = grid_->leafGridView();
      for (auto const& e : elements(gv))
      {
        for (auto const& segment : intersections(gv,e)) {
          if (!segment.boundary())
            continue;

          auto n = segment.centerUnitOuterNormal();
          auto index = segment.boundarySegmentIndex();

          for (int i = 0; i < dow; ++i) {
            if (n[i] < -0.5)
              boundaryIds_[index] = ids[2*i];
            else if (n[i] > 0.5)
              boundaryIds_[index] = ids[2*i+1];
          }
        }
      }
    }


    template <class Indicator,
      REQUIRES(Concepts::Functor<Indicator, int(Domain)>) >
    void setIndicator(Indicator const& indicator)
    {
      auto gv = grid_->leafGridView();
      for (auto const& e : elements(gv))
      {
        for (auto const& segment : intersections(gv,e)) {
          if (!segment.boundary())
            continue;

          auto index = segment.boundarySegmentIndex();
          boundaryIds_[index] = indicator(segment.geometry().center());
        }
      }
    }


    template <class Predicate,
      REQUIRES(Concepts::Functor<Predicate, bool(Domain)>) >
    void setPredicate(Predicate const& pred, BoundaryType id)
    {
      auto gv = grid_->leafGridView();
      for (auto const& e : elements(gv))
      {
        for (auto const& segment : intersections(gv,e)) {
          if (!segment.boundary())
            continue;

          auto index = segment.boundarySegmentIndex();
          if (pred(segment.geometry().center()))
            boundaryIds_[index] = id;
        }
      }
    }


    template <class I>
    using HasBoundaryId = decltype(std::declval<I>().boundaryId());

    void setBoundaryId()
    {
      if (!Dune::Std::is_detected<HasBoundaryId, Segment>::value)
        return;

      auto gv = grid_->leafGridView();
      for (auto const& e : elements(gv))
      {
        for (auto const& segment : intersections(gv,e)) {
          if (!segment.boundary())
            continue;

          if constexpr (Dune::Std::is_detected<HasBoundaryId, Segment>::value) {
            auto index = segment.boundarySegmentIndex();
            boundaryIds_[index] = segment.boundaryId();
          }
        }
      }
    }

    template <class I>
    void setBoundaryIds(std::vector<I> const& ids)
    {
      test_exit(ids.size() == boundaryIds_.size(), "Number of boundary IDs does not match!");
      std::copy(ids.begin(), ids.end(), boundaryIds_.begin());
    }

  private:
    std::shared_ptr<Grid> grid_;
    using Super::boundaryIds_;
  };

} // end namespace AMDiS