pqkfactory.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_PQKFACTORY_HH
#define DUNE_LOCALFUNCTIONS_LAGRANGE_PQKFACTORY_HH
 
#include <map>
 
#include <dune/geometry/type.hh>
 
#include <dune/localfunctions/common/virtualinterface.hh>
#include <dune/localfunctions/common/virtualwrappers.hh>
 
#include <dune/localfunctions/lagrange/p0.hh>
#include <dune/localfunctions/lagrange/pk.hh>
#include <dune/localfunctions/lagrange/q1.hh>
#include <dune/localfunctions/lagrange/qk.hh>
#include <dune/localfunctions/lagrange/prismp1.hh>
#include <dune/localfunctions/lagrange/prismp2.hh>
#include <dune/localfunctions/lagrange/pyramidp1.hh>
#include <dune/localfunctions/lagrange/pyramidp2.hh>
 
namespace Dune
{
 
  template<class D, class R, int d, int k>
  struct DimSpecificPQkLocalFiniteElementFactory
  {
    typedef typename P0LocalFiniteElement<D,R,d>::Traits::LocalBasisType::Traits T;
 
    static LocalFiniteElementVirtualInterface<T>* create(const GeometryType&)
    {
      return nullptr;
    }
  };
 
  template<class D, class R, int k>
  struct DimSpecificPQkLocalFiniteElementFactory<D,R,3,k>
  {
    typedef typename P0LocalFiniteElement<D,R,3>::Traits::LocalBasisType::Traits T;
    typedef PrismP1LocalFiniteElement<D,R> PrismP1;
    typedef PrismP2LocalFiniteElement<D,R> PrismP2;
    typedef PyramidP1LocalFiniteElement<D,R> PyramidP1;
    typedef PyramidP2LocalFiniteElement<D,R> PyramidP2;
 
    static LocalFiniteElementVirtualInterface<T>* create(const GeometryType& gt)
    {
      if ((gt.isPrism())and (k==1))
        return new LocalFiniteElementVirtualImp<PrismP1>(PrismP1());
      if ((gt.isPrism())and (k==2))
        return new LocalFiniteElementVirtualImp<PrismP2>(PrismP2());
      if ((gt.isPyramid())and (k==1))
        return new LocalFiniteElementVirtualImp<PyramidP1>(PyramidP1());
      if ((gt.isPyramid())and (k==2))
        return new LocalFiniteElementVirtualImp<PyramidP2>(PyramidP2());
      return nullptr;
    }
  };
 
 
  template<class D, class R, int dim, int k>
  struct PQkLocalFiniteElementFactory
  {
    typedef typename P0LocalFiniteElement<D,R,dim>::Traits::LocalBasisType::Traits T;
    typedef LocalFiniteElementVirtualInterface<T> FiniteElementType;
    typedef P0LocalFiniteElement<D,R,dim> P0;
    typedef PkLocalFiniteElement<D,R,dim,k> Pk;
    typedef QkLocalFiniteElement<D,R,dim,k> Qk;
 
 
    static FiniteElementType* create(const GeometryType& gt)
    {
      if (k==0)
        return new LocalFiniteElementVirtualImp<P0>(P0(gt));
 
      if (gt.isSimplex())
        return new LocalFiniteElementVirtualImp<Pk>(Pk());
 
      if (gt.isCube())
        return new LocalFiniteElementVirtualImp<Qk>(Qk());
 
      return DimSpecificPQkLocalFiniteElementFactory<D,R,dim,k>::create(gt);
    }
  };
 
 
 
  template<class D, class R, int dim, int k>
  class PQkLocalFiniteElementCache
  {
  protected:
    typedef typename P0LocalFiniteElement<D,R,dim>::Traits::LocalBasisType::Traits T;
    typedef LocalFiniteElementVirtualInterface<T> FE;
    typedef typename std::map<GeometryType,FE*> FEMap;
 
  public:
    typedef FE FiniteElementType;
 
    PQkLocalFiniteElementCache() {}
 
    PQkLocalFiniteElementCache(const PQkLocalFiniteElementCache& other)
    {
      typename FEMap::iterator it = other.cache_.begin();
      typename FEMap::iterator end = other.cache_.end();
      for(; it!=end; ++it)
        cache_[it->first] = (it->second)->clone();
    }
 
    ~PQkLocalFiniteElementCache()
    {
      typename FEMap::iterator it = cache_.begin();
      typename FEMap::iterator end = cache_.end();
      for(; it!=end; ++it)
        delete it->second;
    }
 
    const FiniteElementType& get(const GeometryType& gt) const
    {
      typename FEMap::const_iterator it = cache_.find(gt);
      if (it==cache_.end())
      {
        FiniteElementType* fe = PQkLocalFiniteElementFactory<D,R,dim,k>::create(gt);
        if (fe==0)
          DUNE_THROW(Dune::NotImplemented,"No Pk/Qk like local finite element available for geometry type " << gt << " and order " << k);
 
        cache_[gt] = fe;
        return *fe;
      }
      return *(it->second);
    }
 
  protected:
    mutable FEMap cache_;
 
  };
 
}
 
#endif