FlatPreBasis.hpp

#pragma once

#include <dune/common/version.hh>
#include <dune/functions/functionspacebases/basistags.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>

#include <amdis/Output.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/typetree/FiniteElementType.hpp>

namespace AMDiS
{
  // Convert the index-merging strategy to FlatLexicographic or FlatInterleaved for
  // composite or power bases, respectively.
  template <class PreBasis,
            class MultiIndex = Dune::Functions::FlatMultiIndex<std::size_t>>
  struct FlatPreBasis
  {
    using type = PreBasis;

    template <class PB>
    static type create(PB const& preBasis)
    {
      return {preBasis.gridView()};
    }

    static PreBasis const& create(PreBasis const& preBasis)
    {
      return preBasis;
    }
  };

  template <class PreBasis,
            class MultiIndex = Dune::Functions::FlatMultiIndex<std::size_t>>
  using FlatPreBasis_t = typename FlatPreBasis<PreBasis, MultiIndex>::type;

  template <class PreBasis>
  auto flatPreBasis(PreBasis const& preBasis)
  {
    return FlatPreBasis<PreBasis>::create(preBasis);
  }

  // specialization for Lagrange basis that needs an additional `order` parameter.
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
  template <class GV, class MI, class MultiIndex>
  struct FlatPreBasis<Dune::Functions::LagrangePreBasis<GV,-1,MI>, MultiIndex>
  {
    using type = Dune::Functions::LagrangePreBasis<GV,-1,MultiIndex>;
#else
  template <class GV, class MI, class R, class MultiIndex>
  struct FlatPreBasis<Dune::Functions::LagrangePreBasis<GV,-1,MI,R>, MultiIndex>
  {
    using type = Dune::Functions::LagrangePreBasis<GV,-1,MultiIndex,R>;
#endif

    template <class PB>
    static type create(PB const& preBasis)
    {
      auto node = preBasis.makeNode();
      node.bind(*preBasis.gridView().template begin<0>());
      return {preBasis.gridView(), (unsigned int)(node.finiteElement().localBasis().order())};
    }

    static type const& create(type const& preBasis)
    {
      return preBasis;
    }
  };

  namespace Impl
  {
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
    // NOTE: dirty hack to get access to protected member variable, due to missing
    // function in dune-functions-2.6
    template <class PreBasis>
    class DerivedPreBasis {
    public:
      // specialization for PowerPreBasis
      template <class PB, class SPB = typename PB::SubPreBasis>
      static SPB const& subPreBasis(PB const& pb)
      {
        return access(pb).subPreBasis_;
      }

      // specialization for CompositePreBasis
      template <std::size_t I, class PB, class SPB = typename PB::template SubPreBasis<I>>
      static SPB const& subPreBasis(PB const& pb, Dune::index_constant<I>)
      {
        return std::get<I>(access(pb).subPreBases_);
      }

      template <class T>
      class Accessor : public T { friend class DerivedPreBasis<PreBasis>; };

      template <class T>
      static Accessor<T> const& access(T const& obj) { return static_cast<Accessor<T> const&>(obj); }
    };

    template <class PreBasis, class... Index>
    static auto const& subPreBasis(PreBasis const& preBasis, Index... ii)
    {
      return DerivedPreBasis<PreBasis>::subPreBasis(preBasis,ii...);
    }
#else
    template <class PreBasis, class... Index>
    static auto const& subPreBasis(PreBasis const& preBasis, Index... ii)
    {
      return preBasis.subPreBasis(ii...);
    }
#endif
  } // end namespace Impl


  // specialization for composite bases
  template <class MI, class IMS, class... SPB, class MultiIndex>
  struct FlatPreBasis<Dune::Functions::CompositePreBasis<MI, IMS, SPB...>, MultiIndex>
  {
    using FIMS = Dune::Functions::BasisFactory::FlatLexicographic;
    using type = Dune::Functions::CompositePreBasis<MultiIndex, FIMS, FlatPreBasis_t<SPB, MultiIndex>...>;

    template <class PreBasis>
    static type create(PreBasis const& preBasis)
    {
      return create(preBasis, std::index_sequence_for<SPB...>{});
    }

    template <class PreBasis, std::size_t... I>
    static type create(PreBasis const& preBasis, std::index_sequence<I...>)
    {
      test_warning(std::is_same_v<IMS,FIMS>, "Basis converted into flat index-merging strategy.");
      return {FlatPreBasis<SPB, MultiIndex>::create(Impl::subPreBasis(preBasis,Dune::index_constant<I>{}))...};
    }
  };

  // specialization for flat power bases
  template <class MI, class IMS, class SPB, std::size_t C, class MultiIndex>
  struct FlatPreBasis<Dune::Functions::PowerPreBasis<MI, IMS, SPB, C>, MultiIndex>
  {
    using type = Dune::Functions::PowerPreBasis<MultiIndex, IMS, SPB, C>;

    template <class PreBasis>
    static type create(PreBasis const& preBasis)
    {
      return {FlatPreBasis<SPB, MultiIndex>::create(Impl::subPreBasis(preBasis))};
    }
  };

  // specialization for blocked power bases
  template <class MI, class SPB, std::size_t C, class MultiIndex>
  struct FlatPreBasis<Dune::Functions::PowerPreBasis
    <MI, Dune::Functions::BasisFactory::BlockedInterleaved, SPB, C>, MultiIndex>
  {
    using FIMS = Dune::Functions::BasisFactory::FlatInterleaved;
    using type = Dune::Functions::PowerPreBasis<MultiIndex, FIMS, FlatPreBasis_t<SPB, MultiIndex>, C>;

    template <class PreBasis>
    static type create(PreBasis const& preBasis)
    {
      warning("Basis converted into flat index-merging strategy.");
      return {FlatPreBasis<SPB, MultiIndex>::create(Impl::subPreBasis(preBasis))};
    }
  };

  // specialization for blocked power bases
  template <class MI, class SPB, std::size_t C, class MultiIndex>
  struct FlatPreBasis<Dune::Functions::PowerPreBasis
    <MI, Dune::Functions::BasisFactory::BlockedLexicographic, SPB, C>, MultiIndex>
  {
    using FIMS = Dune::Functions::BasisFactory::FlatLexicographic;
    using type = Dune::Functions::PowerPreBasis<MultiIndex, FIMS, FlatPreBasis_t<SPB, MultiIndex>, C>;

    template <class PreBasis>
    static type create(PreBasis const& preBasis)
    {
      warning("Basis converted into flat index-merging strategy.");
      return {FlatPreBasis<SPB, MultiIndex>::create(Impl::subPreBasis(preBasis))};
    }
  };

} // end namespace AMDiS