pkqkfem.hh

// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_FINITEELEMENTMAP_PKQKFEM_HH
#define DUNE_PDELAB_FINITEELEMENTMAP_PKQKFEM_HH
 
#include <array>
#include <memory>
 
#include <dune/geometry/type.hh>
 
#include <dune/localfunctions/common/virtualwrappers.hh>
#include <dune/localfunctions/common/virtualinterface.hh>
#include "finiteelementmap.hh"
#include "qkfem.hh"
#include "pkfem.hh"
 
namespace Dune {
    namespace PDELab {
 
        namespace {
 
            template<class D, class R, int d, int k>
            struct InitPkQkLocalFiniteElementMap
            {
                template<typename C>
                static void init(C & c, unsigned int order)
                {
                    if (order == k)
                    {
                        typedef Dune::LagrangeCubeLocalFiniteElement<D,R,d,k> QkLFE;
                        typedef Dune::LagrangeSimplexLocalFiniteElement<D,R,d,k> PkLFE;
                        typedef typename C::value_type ptr;
                        c[0] = ptr(new LocalFiniteElementVirtualImp<QkLFE>(QkLFE()));
                        c[1] = ptr(new LocalFiniteElementVirtualImp<PkLFE>(PkLFE()));
                    }
                    else
                        InitPkQkLocalFiniteElementMap<D,R,d,k-1>::init(c,order);
                }
            };
            template<class D, class R, int d>
            struct InitPkQkLocalFiniteElementMap<D,R,d,-1>
            {
                template<typename C>
                    static void init(C & c, unsigned int order)
                {
                    DUNE_THROW(Exception, "Sorry, but we failed to initialize a QkPk FiniteElementMap of order " << order);
                }
            };
        }
 
        template<class D, class R, int d, int maxP=6>
        class PkQkLocalFiniteElementMap
        {
            typedef LocalFiniteElementVirtualInterface<Dune::LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> > > FiniteElementType;
        public:
            typedef FiniteElementMapTraits<FiniteElementType> Traits;
 
            static constexpr int dimension = d;
 
            PkQkLocalFiniteElementMap ()
              : order_(maxP)
            {
                InitPkQkLocalFiniteElementMap<D,R,d,maxP>::init(finiteElements_,maxP);
            }
 
            PkQkLocalFiniteElementMap (unsigned int order)
              : order_(order)
            {
                InitPkQkLocalFiniteElementMap<D,R,d,maxP>::init(finiteElements_,order);
            }
 
            template<class EntityType>
            const typename Traits::FiniteElementType& find (const EntityType& e) const
            {
                typename Dune::GeometryType geoType=e.type();
                return getFEM(geoType);
            }
 
            const typename Traits::FiniteElementType& getFEM (Dune::GeometryType gt) const
            {
                if (gt.isCube())
                {
                    return *(finiteElements_[0]);
                }
                if (gt.isSimplex())
                {
                    return *(finiteElements_[1]);
                }
                DUNE_THROW(Exception, "We can only handle cubes and simplices");
            }
 
            static constexpr bool fixedSize()
            {
                return false;
            }
 
            bool hasDOFs(int codim) const
            {
              switch (codim)
                {
                case 0:
                  return order_ == 0 || order_ > 1;
                case d:
                  return order_ > 0;
                default:
                  return order_ > 1;
                }
            }
 
            std::size_t size(GeometryType gt) const
            {
                DUNE_THROW(NotImplemented, "PkQkLocalFiniteElement is not fixed-size!");
            }
 
            static constexpr std::size_t maxLocalSize()
            {
                return (1<<d);
            }
 
        private:
            std::array< std::shared_ptr<FiniteElementType>, 2 > finiteElements_;
            const std::size_t order_;
        };
    }
}
 
#endif // DUNE_PDELAB_FINITEELEMENTMAP_PKQKFEM_HH