compositebasis.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_COMPOSITEBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_COMPOSITEBASIS_HH
 
#include <tuple>
#include <utility>
 
#include <dune/common/hybridutilities.hh>
#include <dune/common/reservedvector.hh>
#include <dune/common/typeutilities.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/tupleutility.hh>
#include <dune/common/tuplevector.hh>
 
#include <dune/functions/common/staticforloop.hh>
#include <dune/functions/common/type_traits.hh>
#include <dune/functions/common/utility.hh>
#include <dune/functions/functionspacebases/basistags.hh>
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/concepts.hh>
#include <dune/functions/functionspacebases/containerdescriptors.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
 
 
namespace Dune {
namespace Functions {
 
// *****************************************************************************
// This is the reusable part of the composite bases. It contains
//
//   CompositePreBasis
//
// The pre-basis allows to create the others and is the owner of possible shared
// state. These components do _not_ depend on the global basis and local view
// and can be used without a global basis.
// *****************************************************************************
 
 
template<class IMS, class... SPB>
class CompositePreBasis
{
  static const bool isBlocked = std::is_same_v<IMS,BasisFactory::BlockedLexicographic> or std::is_same_v<IMS,BasisFactory::BlockedInterleaved>;
public:
 
  using SubPreBases = std::tuple<SPB...>;
 
  template<std::size_t i>
  using SubPreBasis = std::tuple_element_t<i, SubPreBases>;
 
  using GridView = typename std::tuple_element_t<0, SubPreBases>::GridView;
 
  using size_type = std::size_t;
 
  using IndexMergingStrategy = IMS;
 
protected:
  static const std::size_t children = sizeof...(SPB);
 
  using ChildIndices = std::make_index_sequence<children>;
 
public:
 
  using Node = CompositeBasisNode<typename SPB::Node...>;
 
  static constexpr size_type maxMultiIndexSize = std::max({SPB::maxMultiIndexSize...}) + isBlocked;
  static constexpr size_type minMultiIndexSize = std::min({SPB::minMultiIndexSize...}) + isBlocked;
  static constexpr size_type multiIndexBufferSize = std::max({SPB::multiIndexBufferSize...}) + isBlocked;
 
  template<class... SFArgs,
    disableCopyMove<CompositePreBasis, SFArgs...> = 0,
    enableIfConstructible<std::tuple<SPB...>, SFArgs...> = 0>
  CompositePreBasis(SFArgs&&... sfArgs) :
    subPreBases_(std::forward<SFArgs>(sfArgs)...)
  {
    Hybrid::forEach(subPreBases_, [&](const auto& subPreBasis){
      static_assert(models<Concept::PreBasis<GridView>, std::decay_t<decltype(subPreBasis)>>(), "Subprebases passed to CompositePreBasis does not model the PreBasis concept.");
    });
  }
 
  template<class GV,
    std::enable_if_t<std::conjunction_v<
      std::bool_constant<(children > 1)>,    // Avoid ambiguous constructor if there's only one child
      std::is_same<GV, GridView>,
      std::is_constructible<SPB, GridView>...
    >, int> = 0>
  CompositePreBasis(const GV& gv) :
    subPreBases_(SPB(gv)...)
  {
    Hybrid::forEach(subPreBases_, [&](const auto& subPreBasis){
      static_assert(models<Concept::PreBasis<GridView>, std::decay_t<decltype(subPreBasis)>>(), "Subprebases passed to CompositePreBasis does not model the PreBasis concept.");
    });
  }
 
  void initializeIndices()
  {
    Hybrid::forEach(ChildIndices(), [&](auto i) {
      this->subPreBasis(i).initializeIndices();
    });
  }
 
  const GridView& gridView() const
  {
    return std::get<0>(subPreBases_).gridView();
  }
 
  void update(const GridView& gv)
  {
    Hybrid::forEach(ChildIndices(), [&](auto i) {
      this->subPreBasis(i).update(gv);
    });
  }
 
  Node makeNode() const
  {
    auto node = Node{};
    Hybrid::forEach(ChildIndices(), [&](auto i) {
      node.setChild(this->subPreBasis(i).makeNode(), i);
    });
    return node;
  }
 
  size_type size() const
  {
    return size(Dune::ReservedVector<size_type, multiIndexBufferSize>{});
  }
 
  template<class SizePrefix>
  size_type size(const SizePrefix& prefix) const
  {
    return size(prefix, IndexMergingStrategy{});
  }
 
private:
 
  template<class MultiIndex>
  static void multiIndexPopFront(MultiIndex& M)
  {
    for(std::size_t i=0; i<M.size()-1; ++i)
      M[i] = M[i+1];
    M.resize(M.size()-1);
  }
 
  template<class SizePrefix>
  size_type size(SizePrefix prefix, BasisFactory::BlockedLexicographic) const
  {
    if (prefix.size() == 0)
      return children;
 
    auto front = prefix.front();
    multiIndexPopFront(prefix);
    return Hybrid::switchCases(ChildIndices(), front, [&] (auto i) {
      return this->subPreBasis(i).size(prefix);
    }, []() {
      return size_type(0);
    });
  }
 
  template<class SizePrefix>
  size_type size(SizePrefix prefix, BasisFactory::FlatLexicographic) const
  {
    size_type result = 0;
    if (prefix.size() == 0)
      Hybrid::forEach(ChildIndices(), [&](auto i) {
        result += this->subPreBasis(i).size();
      });
    else {
      staticFindInRange<0, children>([&](auto i) {
          auto firstDigitSize = this->subPreBasis(i).size();
          if (prefix[0] < firstDigitSize)
          {
            result = this->subPreBasis(i).size(prefix);
            return true;
          }
          prefix[0] -= firstDigitSize;
          return false;
        });
    }
    return result;
  }
 
public:
 
  size_type dimension() const
  {
    size_type r=0;
    // Accumulate dimension() for all subprebases
    Hybrid::forEach(ChildIndices(), [&](auto i) {
      r += this->subPreBasis(i).dimension();
    });
    return r;
  }
 
  size_type maxNodeSize() const
  {
    size_type r=0;
    // Accumulate maxNodeSize() for all subprebases
    Hybrid::forEach(ChildIndices(), [&](auto i) {
      r += this->subPreBasis(i).maxNodeSize();
    });
    return r;
  }
 
  template<std::size_t i>
  const SubPreBasis<i>& subPreBasis(Dune::index_constant<i> = {}) const
  {
    return std::get<i>(subPreBases_);
  }
 
  template<std::size_t i>
  SubPreBasis<i>& subPreBasis(Dune::index_constant<i> = {})
  {
    return std::get<i>(subPreBases_);
  }
 
  const auto& subPreBases() const
  {
    return subPreBases_;
  }
 
  template<typename It>
  It indices(const Node& node, It it) const
  {
    return indices(node, it, IndexMergingStrategy{});
  }
 
  auto containerDescriptor() const
  {
    namespace CD = Dune::Functions::ContainerDescriptors;
    if constexpr(std::is_same_v<IMS, BasisFactory::BlockedLexicographic>) {
      return std::apply([&](auto const&... spb) {
        return CD::makeDescriptor(Dune::Functions::containerDescriptor(spb)...);
      }, subPreBases_);
    }
    else if constexpr(std::is_same_v<IMS, BasisFactory::FlatLexicographic>) {
      return CD::Unknown{}; // Not yet implemented
    }
    else
      return CD::Unknown{};
  }
 
private:
 
  template<typename It>
  It indices(const Node& node, It multiIndices, BasisFactory::FlatLexicographic) const
  {
    size_type firstComponentOffset = 0;
    // Loop over all children
    Hybrid::forEach(ChildIndices(), [&](auto child){
      size_type subTreeSize = node.child(child).size();
      // Fill indices for current child into index buffer starting from current
      // buffer position and shift first index component of any index for current
      // child by suitable offset to get lexicographic indices.
      subPreBasis(child).indices(node.child(child), multiIndices);
      for (std::size_t i = 0; i<subTreeSize; ++i)
        multiIndices[i][0] += firstComponentOffset;
      // Increment offset by the size for first index component of the current child
      firstComponentOffset += subPreBasis(child).size();
      // Increment buffer iterator by the number of indices processed for current child
      multiIndices += subTreeSize;
    });
    return multiIndices;
  }
 
  template<class MultiIndex>
  static void multiIndexPushFront(MultiIndex& M, size_type M0)
  {
    M.resize(M.size()+1);
    for(std::size_t i=M.size()-1; i>0; --i)
      M[i] = M[i-1];
    M[0] = M0;
  }
 
  template<typename It>
  It indices(const Node& node, It multiIndices, BasisFactory::BlockedLexicographic) const
  {
    // Loop over all children
    Hybrid::forEach(ChildIndices(), [&](auto child){
      size_type subTreeSize = node.child(child).size();
      // Fill indices for current child into index buffer starting from current position
      subPreBasis(child).indices(node.child(child), multiIndices);
      // Insert child index before first component of all indices of current child.
      for (std::size_t i = 0; i<subTreeSize; ++i)
        this->multiIndexPushFront(multiIndices[i], child);
      // Increment buffer iterator by the number of indices processed for current child
      multiIndices += subTreeSize;
    });
    return multiIndices;
  }
 
  std::tuple<SPB...> subPreBases_;
};
 
 
 
namespace BasisFactory {
 
namespace Imp {
 
template<class IndexMergingStrategy, class... ChildPreBasisFactory>
class CompositePreBasisFactory
{
 
  template<class GridView, class... ChildPreBasis>
  auto makePreBasisFromChildPreBases(const GridView&, ChildPreBasis&&... childPreBasis) const
  {
    return CompositePreBasis<IndexMergingStrategy, std::decay_t<ChildPreBasis>...>(std::forward<ChildPreBasis>(childPreBasis)...);
  }
 
public:
 
  CompositePreBasisFactory(const ChildPreBasisFactory&... childPreBasisFactory) :
    childPreBasisFactories_(childPreBasisFactory...)
  {}
 
  CompositePreBasisFactory(ChildPreBasisFactory&&... childPreBasisFactory) :
    childPreBasisFactories_(std::move(childPreBasisFactory)...)
  {}
 
  template<class GridView>
  auto operator()(const GridView& gridView) const
  {
    // Use std::apply to unpack the tuple childPreBasisFactories_
    return std::apply([&](const auto&... childPreBasisFactory) {
        return this->makePreBasisFromChildPreBases(gridView, childPreBasisFactory(gridView)...);
      }, childPreBasisFactories_);
  }
 
private:
  std::tuple<ChildPreBasisFactory...> childPreBasisFactories_;
};
 
} // end namespace BasisFactory::Imp
 
 
 
template<
  typename... Args,
  std::enable_if_t<Concept::isIndexMergingStrategy<typename LastType<Args...>::type>(),int> = 0>
auto composite(Args&&... args)
{
  // We have to separate the last entry which is the IndexMergingStrategy
  // and the preceding ones, which are the ChildPreBasisFactories
 
  using ArgTuple = std::tuple<std::decay_t<Args>...>;
 
  // Compute number of children and index of the IndexMergingStrategy argument
  constexpr std::size_t children = sizeof...(Args) - 1;
 
  // Use last type as IndexMergingStrategy
  using IndexMergingStrategy = std::tuple_element_t<children, ArgTuple>;
 
  // Index sequence for all but the last entry for partial tuple unpacking
  auto childIndices = std::make_index_sequence<children>{};
 
  // Unpack tuple only for those entries related to children
  return applyPartial([](auto&&... childPreBasisFactory){
    return Imp::CompositePreBasisFactory<IndexMergingStrategy, std::decay_t<decltype(childPreBasisFactory)>...>(std::forward<decltype(childPreBasisFactory)>(childPreBasisFactory)...);
  },
  std::forward_as_tuple(std::forward<Args>(args)...),
  childIndices);
}
 
template<
  typename... Args,
  std::enable_if_t<not Concept::isIndexMergingStrategy<typename LastType<Args...>::type>(),int> = 0>
auto composite(Args&&... args)
{
  return Imp::CompositePreBasisFactory<BasisFactory::BlockedLexicographic, std::decay_t<Args>...>(std::forward<Args>(args)...);
}
 
} // end namespace BasisFactory
 
// Backward compatibility
namespace [[deprecated("Will be removed after Dune 2.10")]] BasisBuilder {
 
  using namespace BasisFactory;
 
}
 
 
 
} // end namespace Functions
} // end namespace Dune
 
 
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_COMPOSITEBASIS_HH