finiteelementmap.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
#ifndef DUNE_PDELAB_FINITEELEMENTMAP_FINITEELEMENTMAP_HH
#define DUNE_PDELAB_FINITEELEMENTMAP_FINITEELEMENTMAP_HH
 
#include <dune/common/deprecated.hh>
#include <dune/pdelab/common/exceptions.hh>
 
#include <dune/geometry/referenceelements.hh>
 
namespace Dune {
  namespace PDELab {
 
 
    class FiniteElementMapError : public Exception {};
    class VariableElementSize : public FiniteElementMapError {};
    class InvalidGeometryType : public FiniteElementMapError {};
 
    template<class T>
    struct FiniteElementMapTraits
    {
      typedef T FiniteElementType;
 
      typedef T FiniteElement;
    };
 
    template<class T>
    struct LocalFiniteElementMapTraits : FiniteElementMapTraits<T> {};
 
    template<class T, class Imp>
    class LocalFiniteElementMapInterface
    {
    public:
      typedef T Traits;
 
      template<class EntityType>
      const typename Traits::FiniteElementType&
      find (const EntityType& e) const = delete;
 
      static constexpr int dimension = -1;
      static constexpr bool fixedSize() = delete;
      static constexpr std::size_t size(GeometryType gt) = delete;
      static constexpr bool hasDOFs(int codim) = delete;
      static constexpr std::size_t maxLocalSize() = delete;
    };
 
    template<class Imp, int dim_>
    class SimpleLocalFiniteElementMap :
      public LocalFiniteElementMapInterface<LocalFiniteElementMapTraits<Imp>,
                                            SimpleLocalFiniteElementMap<Imp,dim_> >
    {
    public:
      typedef LocalFiniteElementMapTraits<Imp> Traits;
 
      SimpleLocalFiniteElementMap ()
      {}
 
      SimpleLocalFiniteElementMap (const Imp& imp_) : imp(imp_)
      {}
 
      template<class EntityType>
      const typename Traits::FiniteElementType& find (const EntityType& e) const
      {
        return imp;
      }
 
      static constexpr int dimension = dim_;
 
    private:
      Imp imp; // create once
    };
 
    template<typename GV, typename FE, typename Imp>
    class EdgeS0LocalFiniteElementMap
      : public LocalFiniteElementMapInterface<
          LocalFiniteElementMapTraits<FE>,
          Imp
        >
    {
      typedef typename GV::IndexSet IndexSet;
      static const int dim = GV::dimension;
 
    public:
      typedef LocalFiniteElementMapTraits<FE> Traits;
 
      static constexpr int dimension = GV::dimension;
 
      EdgeS0LocalFiniteElementMap (const GV& gv_)
        : gv(gv_), orient(gv_.size(0))
      {
        using ct = typename GV::Grid::ctype;
        auto& refElem =
          ReferenceElements<ct, dim>::general(FE().type());
 
        auto &idSet = gv.grid().globalIdSet();
 
        // create all variants
        variant.resize(1 << refElem.size(dim-1));
        for (std::size_t i=0; i<variant.size(); i++)
          variant[i] = FE(i);
 
        // compute orientation for all elements
        auto& indexSet = gv.indexSet();
 
        // loop once over the grid
        for (const auto& element : elements(gv))
          {
            auto elemid = indexSet.index(element);
            orient[elemid] = 0;
 
            std::vector<typename GV::Grid::GlobalIdSet::IdType> vid(refElem.size(dim));
            for(std::size_t i = 0; i < vid.size(); ++i)
              vid[i] = idSet.subId(element, i, dim);
 
            // loop over all edges of the element
            for(int i = 0; i < refElem.size(dim-1); ++i) {
              auto v0 = refElem.subEntity(i, dim-1, 0, dim);
              auto v1 = refElem.subEntity(i, dim-1, 1, dim);
              // if (edge orientation in refelement) != (edge orientation in indexset)
              if((v0 > v1) != (vid[v0] > vid[v1]))
                orient[elemid] |= 1 << i;
            }
          }
      }
 
      template<class EntityType>
      const typename Traits::FiniteElementType& find (const EntityType& e) const
      {
        return variant[orient[gv.indexSet().index(e)]];
      }
 
    private:
      GV gv;
      std::vector<FE> variant;
      std::vector<unsigned char> orient;
    };
 
    template<typename GV, typename FE, typename Imp, std::size_t Variants>
    class RTLocalFiniteElementMap :
      public LocalFiniteElementMapInterface<
        LocalFiniteElementMapTraits<FE>,
        Imp >
    {
      typedef FE FiniteElement;
      typedef typename GV::IndexSet IndexSet;
 
    public:
      typedef LocalFiniteElementMapTraits<FE> Traits;
 
      static constexpr int dimension = GV::dimension;
 
      RTLocalFiniteElementMap(const GV& gv_)
        : gv(gv_), is(gv_.indexSet()), orient(gv_.size(0))
      {
        // create all variants
        for (std::size_t i = 0; i < Variants; i++)
        {
          variant[i] = FiniteElement(i);
        }
 
        // compute orientation for all elements
        // loop once over the grid
        for(const auto& cell : elements(gv))
        {
          auto myId = is.index(cell);
          orient[myId] = 0;
 
          for (const auto& intersection : intersections(gv,cell))
          {
            if (intersection.neighbor()
                && is.index(intersection.outside()) > myId)
            {
              orient[myId] |= 1 << intersection.indexInInside();
            }
          }
        }
      }
 
      template<class EntityType>
      const typename Traits::FiniteElementType& find(const EntityType& e) const
      {
        return variant[orient[is.index(e)]];
      }
 
    private:
      GV gv;
      std::array<FiniteElement,Variants> variant;
      const IndexSet& is;
      std::vector<unsigned char> orient;
    };
 
 
  } // namespace PDELab
} // namespace Dune
 
#endif // DUNE_PDELAB_FINITEELEMENTMAP_FINITEELEMENTMAP_HH