hierarchicallagrangebasis.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_HIERARCHICALLAGRANGEBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_HIERARCHICALLAGRANGEBASIS_HH
 
#include <dune/common/exceptions.hh>
#include <dune/localfunctions/hierarchical/hierarchicalp2.hh>
 
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/grid/common/mcmgmapper.hh>
 
namespace Dune {
  namespace Functions {
 
    // *****************************************************************************
    // Implementation for Hierarchical Lagrange Basis
    //
    // -- only order k=2 is implemented up to now --
    // -- currently only supports simplex grids --
    //
    // This is the reusable part of the HierarchicalLagrangeBasis. It contains
    //
    //  HierarchicalLagrangePreBasis
    //  HierarchicalLagrangeNode
    //
    // 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 can be
    // used without a global basis.
    // *****************************************************************************
 
    template<typename GV, int k, typename R=double>
    class HierarchicalLagrangeNode;
 
    template<typename GV, int k, class MI, typename R=double>
    class HierarchicalLagrangePreBasis;
 
    template<typename GV, int k, class MI, typename R>
    class HierarchicalLagrangePreBasis
    {
      static const int dim = GV::dimension;
 
    public:
 
      using GridView = GV;
 
      using size_type = std::size_t;
 
      using Node = HierarchicalLagrangeNode<GV, k, R>;
 
      using IndexSet = Impl::DefaultNodeIndexSet<HierarchicalLagrangePreBasis>;
 
      using MultiIndex = MI;
 
      using SizePrefix = Dune::ReservedVector<size_type, 1>;
 
      HierarchicalLagrangePreBasis(const GridView& gv) : gridView_(gv) , mcmgMapper_(gv,p2Layout())
      {}
 
      void initializeIndices()
      {}
 
      const GridView& gridView() const
      {
        return gridView_;
      }
 
      void update (const GridView& gv)
      {
        gridView_ = gv;
        mcmgMapper_.update(gv);
      }
 
      Node makeNode() const
      {
        return Node{};
      }
 
      IndexSet makeIndexSet() const
      {
        return IndexSet{*this};
      }
 
      size_type size() const
      {
        return mcmgMapper_.size();
      }
 
      size_type size(const SizePrefix prefix) const
      {
        assert(prefix.size() == 0 || prefix.size() == 1);
        return (prefix.size() == 0) ? size() : 0;
      }
 
      size_type dimension() const
      {
        return size();
      }
 
      size_type maxNodeSize() const
      {
        // That cast to unsigned int is necessary because GV::dimension is an enum
        return Dune::binomial(std::size_t(order() + (unsigned int)GV::dimension),std::size_t(order()));
      }
 
      template<typename It>
      It indices(const Node& node, It it) const
      {
        for (size_type i = 0, end = node.finiteElement().size() ; i < end ; ++it, ++i)
        {
          Dune::LocalKey localKey = node.finiteElement().localCoefficients().localKey(i);
          const auto& element = node.element();
 
          *it = {{ (size_type)(mcmgMapper_.subIndex(element,localKey.subEntity(),localKey.codim())) }};
        }
        return it;
      }
 
    protected:
      GridView gridView_;
 
      unsigned int order() const
      {
        return 2;
      }
 
      MultipleCodimMultipleGeomTypeMapper<GridView> mcmgMapper_;
 
    private:
      static auto p2Layout()
      {
        return [](Dune::GeometryType type, int gridDim)
               {
                 if (type.isVertex())
                   return 1;
                 if (type.isLine())
                   return 1;
                 if (type.isTriangle())
                   return 0;
                 assert(type.isTetrahedron());
                 return 0;
               };
      }
    };
 
 
 
    template<typename GV, int k, typename R>
    class HierarchicalLagrangeNode :
      public LeafBasisNode
    {
      static const int dim = GV::dimension;
 
    public:
 
      using size_type = std::size_t;
      using Element = typename GV::template Codim<0>::Entity;
      using FiniteElement = HierarchicalP2LocalFiniteElement<typename GV::ctype,R,dim>;
 
      HierarchicalLagrangeNode() :
        finiteElement_(),
        element_(nullptr)
      {}
 
      const Element& element() const
      {
        return *element_;
      }
 
      const FiniteElement& finiteElement() const
      {
        return finiteElement_;
      }
 
      void bind(const Element& e)
      {
        element_ = &e;
 
        if (e.type() != finiteElement_.type())
          DUNE_THROW(Dune::Exception,
                     "HierarchicalLagrange-elements do not exist for elements of type " << e.type());
 
        this->setSize(finiteElement_.size());
      }
 
    protected:
 
      unsigned int order() const
      {
        return 2;
      }
 
      const FiniteElement finiteElement_;
      const Element* element_;
    };
 
 
 
    namespace BasisFactory {
 
      namespace Impl {
 
        template<int k, typename R=double>
        class HierarchicalLagrangePreBasisFactory
        {
 
        public:
          static const std::size_t requiredMultiIndexSize = 1;
 
          template<class MultiIndex, class GridView>
          auto makePreBasis(const GridView& gridView) const
          {
            return HierarchicalLagrangePreBasis<GridView, k, MultiIndex, R>(gridView);
          }
 
        };
      } // end namespace BasisFactory::Impl
 
 
      template<std::size_t k, typename R=double>
      auto hierarchicalLagrange()
      {
        return Impl::HierarchicalLagrangePreBasisFactory<k,R>();
      }
 
    } // end namespace BasisFactory
 
    template<typename GV, int k, typename R=double>
    using HierarchicalLagrangeBasis = DefaultGlobalBasis<HierarchicalLagrangePreBasis<GV, k, FlatMultiIndex<std::size_t>, R> >;
 
  } // end namespace Functions
} // end namespace Dune
 
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_HIERARCHICALLAGRANGEBASIS_HH