Twist.hpp

#pragma once

#include <vector>

#include <dune/geometry/referenceelements.hh>
#include <dune/localfunctions/common/localkey.hh>

#include <amdis/Output.hpp>

namespace AMDiS
{
  template <class IdSet, int dim>
  class Twist
  {
    using IdType = typename IdSet::IdType;

    using RefElements = typename Dune::ReferenceElements<double,dim>;
    using RefElement = typename RefElements::ReferenceElement;

  public:
    Twist(IdSet const& idSet)
      : idSet_(idSet)
    {}

    template <class Element>
    void bind(Element const& element)
    {
      static_assert(dim == Element::mydimension, "");
      refElem_ = &RefElements::general(element.type());

      ids_.resize(refElem_->size(dim));
      for (int i = 0; i < refElem_->size(dim); ++i)
        ids_[i] = idSet_.subId(element, i, dim);
    }


    unsigned int get(Dune::LocalKey const& localKey, unsigned int size) const
    {
      int subDim = dim - localKey.codim();
      if (subDim == 1 && 1 < dim) // facet == edge
        return id(localKey,0) < id(localKey,1) ? localKey.index() : size - localKey.index() - 1;
      else if (subDim == 2 && 2 < dim) { // facet == face
        test_exit(size == 1, "Bases with more then one DoF per subentity are not yet supported.");
        return localKey.index();
      }

      return localKey.index();
    }

  private:
    IdType const& id(Dune::LocalKey const& localKey, int ii) const
    {
      return ids_[refElem_->subEntity(localKey.subEntity(), localKey.codim(), ii, dim)];
    }

  private:
    IdSet const& idSet_;
    RefElement const* refElem_ = nullptr;
    std::vector<IdType> ids_;
  };

} // end namespace AMDiS