ElementGridFunction.hpp

#pragma once
 
#include <functional>
#include <type_traits>
 
#include <dune/common/exceptions.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
 
#include <amdis/ElementVector.hpp>
#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/gridfunctions/ConstantGridFunction.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>
 
namespace AMDiS
{
#ifndef DOXYGEN
  template <class Signature, class IndexSet, class Container>
  class ElementLocalFunction;
#endif
 
  template <class R, class D, class IndexSet, class Container>
  class ElementLocalFunction<R(D), IndexSet, Container>
  {
  public:
    using Domain = D;
 
    using Range = R;
 
    enum { hasDerivative = true };
 
    using Element = typename IndexSet::template Codim<0>::Entity;
 
  public:
    ElementLocalFunction(IndexSet const& indexSet, Container const& container)
      : indexSet_(&indexSet)
      , container_(&container)
    {}
 
    Element const& localContext() const
    {
      assert(!!element_);
      return *element_;
    }
 
    void bind(Element const& element)
    {
      element_.emplace(element);
      value_ = (*container_)[indexSet_->index(element)];
    }
 
    void unbind() {}
 
    bool bound() const
    {
      return !!element_;
    }
 
    Range const& operator()(Domain const& /*local*/) const
    {
      return value_;
    }
 
    template <class Type>
    auto makeDerivative(Type const& /*type*/) const
    {
      using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
      using DerivativeRange = typename DerivativeTraits<RawSignature,Type>::Range;
      DerivativeRange diff(0);
      ConstantLocalFunction<DerivativeRange(D), Element, DerivativeRange> df{diff};
      if (bound())
        df.bind(*element_);
      return df;
    }
 
    int order() const
    {
      return 0;
    }
 
  private:
    IndexSet const* indexSet_;
    Container const* container_;
    std::optional<Element> element_;
    Range value_ = 0;
  };
 
 
 
  template <class GridView, class Container>
  class ElementGridFunction
  {
  public:
    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
    using Domain = typename EntitySet::GlobalCoordinate;
    using Range = typename Container::value_type;
 
    enum { hasDerivative = false };
 
  private:
    using LocalDomain = typename EntitySet::LocalCoordinate;
    using LocalFunction = ElementLocalFunction<Range(LocalDomain), typename GridView::IndexSet, Container>;
 
  public:
    ElementGridFunction(GridView const& gridView, Container const& container)
      : entitySet_(gridView)
      , container_(&container)
    {}
 
    template <class G, class T,
      std::enable_if_t<std::is_same_v<typename G::LeafGridView,GridView>, int> = 0,
      std::enable_if_t<std::is_same_v<typename ElementVector<G,T>::Data, Container>, int> = 0>
    ElementGridFunction(ElementVector<G,T> const& elementVector)
      : entitySet_(elementVector.gridView())
      , container_(&elementVector.data())
    {}
 
    Range operator()(Domain const& /*x*/) const
    {
      DUNE_THROW(Dune::NotImplemented, "Global evaluation not yet implemented");
      return Range(0);
    }
 
    EntitySet const& entitySet() const
    {
      return entitySet_;
    }
 
    LocalFunction makeLocalFunction() const
    {
      return {entitySet_.gridView().indexSet(), *container_};
    }
 
  private:
    EntitySet entitySet_;
    Container const* container_;
  };
 
  // deduction guide
  template <class G, class T>
  ElementGridFunction(ElementVector<G,T> const&)
    -> ElementGridFunction<typename G::LeafGridView, typename ElementVector<G,T>::Data>;
 
  template <class G, class T>
  auto valueOf(ElementVector<G,T> const& ev)
  {
    return ElementGridFunction{ev.gridView(), ev.data()};
  }
 
} // end namespace AMDiS