ConstantGridFunction.hpp

#pragma once

#include <functional>
#include <type_traits>

#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>

#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/gridfunctions/AnalyticGridFunction.hpp>
#include <amdis/gridfunctions/GridFunction.hpp>

namespace AMDiS
{
#ifndef DOXYGEN
  template <class Signature, class LocalContext, class Function>
  class ConstantLocalFunction;
#endif

  template <class R, class D, class LC, class T>
  class ConstantLocalFunction<R(D), LC, T>
  {
  public:
    using Domain = D;

    using Range = R;

    enum { hasDerivative = true };

  private:
    using LocalContext = LC;
    using Geometry = typename LC::Geometry;

  public:
    ConstantLocalFunction(T const& value)
      : value_(value)
    {}

    void bind(LocalContext const& element)
    {
      localContext_.emplace(element);
    }

    void unbind() { /* do nothing */ }

    bool bound() const 
    {
      return !!localContext_;
    }
    
    LocalContext const& localContext() const
    {
      assert( !!localContext_ );
      return *localContext_;
    }

    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);
      auto df = ConstantLocalFunction<DerivativeRange(D), LC, DerivativeRange>{diff};
      // bind derivative if this is bound
      if (bound())
        df.bind(*localContext_);
      return df;
    }

    int order() const
    {
      return 0;
    }

  private:
    T value_;
    std::optional<LocalContext> localContext_;
  };



  template <class T, class GridView>
  class ConstantGridFunction
  {
  public:
    using EntitySet = Dune::Functions::GridViewEntitySet<GridView, 0>;
    using Domain = typename EntitySet::GlobalCoordinate;
    using Range = Underlying_t<T>;

    enum { hasDerivative = false };

  private:
    using Element = typename EntitySet::Element;
    using LocalDomain = typename EntitySet::LocalCoordinate;
    using LocalFunction = ConstantLocalFunction<Range(LocalDomain), Element, T>;

  public:
    ConstantGridFunction(T const& value, GridView const& gridView)
      : value_(value)
      , entitySet_(gridView)
    {}

    Range const& operator()(Domain const& /*x*/) const
    {
      return value_;
    }

    EntitySet const& entitySet() const
    {
      return entitySet_;
    }

    LocalFunction makeLocalFunction() const
    {
      return {value_};
    }

  private:
    T value_;
    EntitySet entitySet_;
  };


  namespace Concepts
  {
    namespace Definition
    {
      template <class T>
      struct ConstantToGridFunction
        : std::is_arithmetic<T> {};

      template <class T>
      struct ConstantToGridFunction<std::reference_wrapper<T>>
        : ConstantToGridFunction<T> {};

      template <class T, int N>
      struct ConstantToGridFunction<Dune::FieldVector<T, N>>
        : ConstantToGridFunction<T> {};

      template <class T, int N, int M>
      struct ConstantToGridFunction<Dune::FieldMatrix<T, N, M>>
        : ConstantToGridFunction<T> {};

      template <class T, int N>
      struct ConstantToGridFunction<Dune::DiagonalMatrix<T, N>>
        : ConstantToGridFunction<T> {};

    } // end namespace Definition


    template <class T>
    constexpr bool ConstantToGridFunction =
      Definition::ConstantToGridFunction<T>::value;

  } // end namespace Concepts


  template <class Value>
  struct GridFunctionCreator<Value, std::enable_if_t<Concepts::ConstantToGridFunction<Value>>>
  {
    template <class GridView>
    static auto create(Value const& value, GridView const& gridView)
    {
      return ConstantGridFunction<Value,GridView>{value, gridView};
    }
  };

} // end namespace AMDiS