p2hierarchicalbasis.hh

#ifndef P2_HIERARCHICAL_BASIS_HH
#define P2_HIERARCHICAL_BASIS_HH

/**
   @file
   @brief

   @author
 */

#include <dune/common/version.hh>

#include <dune/localfunctions/hierarchical/hierarchicalp2.hh>
#include <dune/localfunctions/hierarchical/hierarchicalprismp2.hh>

#include <dune/fufem/functionspacebases/functionspacebasis.hh>
#include <dune/fufem/functionspacebases/p2nodalbasis.hh>


template <class GV, class RT=double>
class P2HierarchicalBasis :
    public FunctionSpaceBasis<GV, RT,
        Dune::LocalFiniteElementVirtualInterface<
#if DUNE_VERSION_GTE(DUNE_LOCALFUNCTIONS, 2, 7)
            typename Dune::Impl::LagrangeSimplexLocalBasis<typename GV::Grid::ctype, RT, GV::dimension, 1>::Traits
#else
            typename Dune::P1LocalBasis<typename GV::Grid::ctype, RT, GV::dimension>::Traits
#endif
          >
>
{
    protected:
        typedef typename GV::Grid::ctype ctype;

#if DUNE_VERSION_GTE(DUNE_LOCALFUNCTIONS, 2, 7)
        using P1Traits = typename Dune::Impl::LagrangeSimplexLocalBasis<typename GV::Grid::ctype, RT, GV::dimension, 1>::Traits;
#else
        using P1Traits = typename Dune::P1LocalBasis<typename GV::Grid::ctype, RT, GV::dimension>::Traits;
#endif

        typedef Dune::LocalFiniteElementVirtualInterface<P1Traits> LFE;

        typedef FunctionSpaceBasis<GV, RT, LFE> Base;
        typedef typename Base::Element Element;

        using Base::dim;
        using Base::gridview_;


    public:
        typedef typename Base::GridView GridView;
        typedef typename Base::ReturnType ReturnType;
        typedef typename Base::LocalFiniteElement LocalFiniteElement;
        typedef typename Base::LinearCombination LinearCombination;

        P2HierarchicalBasis(const GridView& gridview) :
            Base(gridview),
            mapper_(gridview)
        {}

        size_t size() const
        {
            return mapper_.size();
        }

        void update(const GridView& gridview)
        {
            Base::update(gridview);
            mapper_.update();
        }

        const LocalFiniteElement& getLocalFiniteElement(const Element& e) const
        {
            if (e.type().isSimplex())
                return localSimplexFE_;
            else if (e.type().isPrism())
                // without the cast this doesn't compile for dim!=3
                return dynamic_cast<const LocalFiniteElement&>(localPrismFE_);
            else
                DUNE_THROW(Dune::NotImplemented, "There is no shape function set for a " << e.type());
        }

        int index(const Element& e, const int i) const
        {
            return mapper_.index(e, getLocalFiniteElement(e).localCoefficients().localKey(i));
        }

    protected:
        const Dune::LocalFiniteElementVirtualImp<Dune::HierarchicalPrismP2LocalFiniteElement<ctype, RT> > localPrismFE_;
        const Dune::LocalFiniteElementVirtualImp<Dune::HierarchicalP2LocalFiniteElement<ctype, RT, dim> > localSimplexFE_;
        P2BasisMapper<GV> mapper_;
};

#endif