variableqkdgfem.hh

// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_FINITEELEMENTMAP_VARIABLEQKDGFEM_HH
#define DUNE_PDELAB_FINITEELEMENTMAP_VARIABLEQKDGFEM_HH
 
#include <array>
#include <memory>
 
#include <dune/geometry/type.hh>
 
#include <dune/localfunctions/common/virtualwrappers.hh>
#include "finiteelementmap.hh"
#include "qkdg.hh"
 
namespace Dune {
  namespace PDELab {
 
    namespace {
      template<class D, class R, int d, int p>
      struct InitVariableQkDGLocalFiniteElementMap
      {
        template<typename C>
        static void init(C & c)
        {
          typedef Dune::QkDGLagrangeLocalFiniteElement<D,R,p,d> LFE;
          typedef typename C::value_type ptr;
          c[p] = ptr(new LocalFiniteElementVirtualImp<LFE>(LFE()));
 
          InitVariableQkDGLocalFiniteElementMap<D,R,d,p-1>::init(c);
        }
      };
      template<class D, class R, int d>
      struct InitVariableQkDGLocalFiniteElementMap<D,R,d,-1>
      {
        template<typename C>
        static void init(C & c) {}
      };
    }
 
    template<class M, class D, class R, int d, int maxP=6>
    class VariableQkDGLocalFiniteElementMap
    {
      typedef typename QkDGLagrangeLocalFiniteElement<D,R,0,d>::Traits::LocalBasisType::Traits T;
      typedef LocalFiniteElementVirtualInterface<T> FiniteElementType;
    public:
      typedef FiniteElementMapTraits<FiniteElementType> Traits;
 
      static constexpr int dimension = d;
 
      VariableQkDGLocalFiniteElementMap (const M & m, unsigned int defaultP) :
        mapper_(m), polOrder_(mapper_.size(), defaultP), defaultP_(defaultP)
      {
        InitVariableQkDGLocalFiniteElementMap<D,R,d,maxP>::init(finiteElements_);
      }
 
      template<class EntityType>
      const typename Traits::FiniteElementType& find (const EntityType& e) const
      {
        return getFEM(getOrder(e));
      }
 
      const typename Traits::FiniteElementType& getFEM (unsigned int p) const
      {
        return *(finiteElements_[p]);
      }
 
      const typename Traits::FiniteElementType& getFEM () const
      {
        return *(finiteElements_[defaultP_]);
      }
 
      template<class EntityType>
      void setOrder (const EntityType& e, unsigned int p)
      {
        assert(p <= maxP);
        unsigned int i = mapper_.map(e);
        polOrder_[i] = p;
      }
 
      template<class EntityType>
      unsigned int getOrder (const EntityType& e) const
      {
        unsigned int i = mapper_.map(e);
        unsigned int p = polOrder_[i];
        assert(p <= maxP);
        return p;
      }
 
      static constexpr bool fixedSize()
      {
        return false;
      }
 
      static constexpr bool hasDOFs(int codim)
      {
        return codim == 0;
      }
 
      std::size_t size(GeometryType gt) const
      {
        DUNE_THROW(Dune::Exception,"This should not be called!");
      }
 
      std::size_t maxLocalSize() const
      {
        return getFEM(maxP).localCoefficients().size();
      }
 
    private:
      const M & mapper_;
      std::vector<unsigned char> polOrder_;
      unsigned int defaultP_;
      std::array< std::shared_ptr<FiniteElementType>, maxP+1 > finiteElements_;
    };
 
 
  }
}
 
#endif // DUNE_PDELAB_FINITEELEMENTMAP_VARIABLEQKDGFEM_HH