monomial.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_LOCALFUNCTIONS_MONOMIAL_HH
#define DUNE_LOCALFUNCTIONS_MONOMIAL_HH
 
#include <cassert>
#include <cstddef>
#include <cstdlib>
#include <memory>
#include <vector>
 
#include <dune/geometry/type.hh>
 
#include "common/localfiniteelementtraits.hh"
#include "common/localtoglobaladaptors.hh"
#include "monomial/monomiallocalbasis.hh"
#include "monomial/monomiallocalcoefficients.hh"
#include "monomial/monomiallocalinterpolation.hh"
 
namespace Dune
{
 
 
  template<class D, class R, int d, int p>
  class MonomialLocalFiniteElement
  {
    enum { static_size = MonomImp::Size<d,p>::val };
 
  public:
    typedef LocalFiniteElementTraits<
        MonomialLocalBasis<D,R,d,p>,
        MonomialLocalCoefficients<static_size>,
        MonomialLocalInterpolation<MonomialLocalBasis<D,R,d,p>,static_size>
        > Traits;
 
    MonomialLocalFiniteElement (const GeometryType &gt_)
      : basis(), interpolation(gt_, basis), gt(gt_)
    {}
 
    const typename Traits::LocalBasisType& localBasis () const
    {
      return basis;
    }
 
    const typename Traits::LocalCoefficientsType& localCoefficients () const
    {
      return coefficients;
    }
 
    const typename Traits::LocalInterpolationType& localInterpolation () const
    {
      return interpolation;
    }
 
    unsigned int size () const
    {
      return basis.size();
    }
 
    GeometryType type () const
    {
      return gt;
    }
 
  private:
    MonomialLocalBasis<D,R,d,p> basis;
    MonomialLocalCoefficients<static_size> coefficients;
    MonomialLocalInterpolation<MonomialLocalBasis<D,R,d,p>,static_size> interpolation;
    GeometryType gt;
  };
 
 
  template<class Geometry, class RF, std::size_t p>
  class MonomialFiniteElementFactory {
    typedef typename Geometry::ctype DF;
    static const std::size_t dim = Geometry::mydimension;
 
    typedef MonomialLocalFiniteElement<DF, RF, dim, p> LocalFE;
 
    std::vector<std::shared_ptr<const LocalFE> > localFEs;
 
    void init(const GeometryType &gt) {
      std::size_t index = gt.id() >> 1;
      if(localFEs.size() <= index)
        localFEs.resize(index+1);
      localFEs[index].reset(new LocalFE(gt));
    }
 
  public:
    typedef ScalarLocalToGlobalFiniteElementAdaptor<LocalFE, Geometry>
    FiniteElement;
 
 
    template<class ForwardIterator>
    MonomialFiniteElementFactory(const ForwardIterator &begin,
                              const ForwardIterator &end)
    {
      for(ForwardIterator it = begin; it != end; ++it)
        init(*it);
    }
 
 
    MonomialFiniteElementFactory(const GeometryType &gt)
    { init(gt); }
 
 
    MonomialFiniteElementFactory() {
      static_assert(dim <= 3, "MonomFiniteElementFactory knows the "
                    "available geometry types only up to dimension 3");
 
      GeometryType gt;
      switch(dim) {
      case 0 :
        gt = Dune::GeometryTypes::vertex;        init(gt);
        break;
      case 1 :
        gt = Dune::GeometryTypes::line;          init(gt);
        break;
      case 2 :
        gt = Dune::GeometryTypes::triangle;      init(gt);
        gt = Dune::GeometryTypes::quadrilateral; init(gt);
        break;
      case 3 :
        gt = Dune::GeometryTypes::tetrahedron;   init(gt);
        gt = Dune::GeometryTypes::pyramid;       init(gt);
        gt = Dune::GeometryTypes::prism;         init(gt);
        gt = Dune::GeometryTypes::hexahedron;    init(gt);
        break;
      default :
        // this should never happen -- it should be caught by the static
        // assert above.
        std::abort();
      };
    }
 
 
    const FiniteElement make(const Geometry& geometry) {
      std::size_t index = geometry.type().id() >> 1;
      assert(localFEs.size() > index && localFEs[index]);
      return FiniteElement(*localFEs[index], geometry);
    }
  };
}
 
#endif // DUNE_LOCALFUNCTIONS_MONOMIAL_HH