twistprovider.hh

#ifndef DUNE_FEM_TWISTPROVIDER_HH
#define DUNE_FEM_TWISTPROVIDER_HH
 
//- System includes
#include <cassert>
 
#include <map>
#include <memory>
#include <vector>
 
//- Dune includes
#include <dune/geometry/referenceelements.hh>
 
#include <dune/fem/quadrature/quadrature.hh>
 
#include <dune/fem/storage/singleton.hh>
 
//- Local includes
#include "pointmapper.hh"
#include "topology.hh"
 
namespace Dune
{
 
  namespace Fem
  {
 
    // Forward declaration
    template <class ct, int dim>
    class TwistProvider;
    template <class ct, int dim>
    class TwistMapperCreator;
    template <class ct, int dim>
    class TwistStorage;
 
 
 
 
    template <class ct, int dim>
    class TwistStorage
    {
      typedef CachingTraits<ct, dim> Traits;
    public:
      typedef typename Traits::PointType PointType;
      typedef typename Traits::PointVectorType PointVectorType;
      typedef typename Traits::MapperType MapperType;
 
    public:
      explicit TwistStorage(int minTwist, int maxTwist);
 
      void addMapper(const MapperType& mapper, int twist);
 
      size_t addPoint(const PointType& points);
 
      const MapperType& getMapper(int twist) const;
 
      const PointVectorType& getPoints() const;
 
      int minTwist() const { return minTwist_; }
 
      int maxTwist() const { return maxTwist_; }
 
    private:
      typedef typename Traits::MapperVectorType MapperVectorType;
 
    private:
      MapperVectorType mappers_;
      PointVectorType points_;
 
      int minTwist_;
      int maxTwist_;
    };
 
 
 
    // TiwstProvider
    // -------------
 
    template <class ct, int dim>
    class TwistProvider
    {
      typedef CachingTraits<ct, dim> Traits;
    public:
      typedef typename Traits::QuadratureType QuadratureType;
      typedef typename TwistMapperCreator<ct, dim>::TwistStorageType TwistStorageType;
 
    public:
      static const TwistStorageType& getTwistStorage(const QuadratureType& quad);
 
    private:
      typedef std::vector< const TwistStorageType* > MapperContainerType;
      typedef typename MapperContainerType::iterator IteratorType;
 
    private:
      // singleton class holding map with storages
      class MapperContainer
      {
        // instance of map
        MapperContainerType mappers_;
 
      public:
        MapperContainer() : mappers_(100, (TwistStorageType*) 0)
        {}
 
        ~MapperContainer()
        {
          IteratorType endit = mappers_.end();
          for(IteratorType it = mappers_.begin(); it != endit; ++it)
          {
            delete (*it);
          }
        }
 
        friend class Dune::Fem::Singleton< MapperContainerType >;
        static MapperContainerType& instance()
        {
          return Singleton< MapperContainer > :: instance().mappers_;
        }
      };
    };
 
 
 
    // TwistMapperStrategy
    // -------------------
 
    template <class ct, int dim>
    struct TwistMapperStrategy
    {
      typedef FieldMatrix<ct, dim+1, dim> MatrixType;
 
    public:
      TwistMapperStrategy(int minTwist, int maxTwist) :
        minTwist_(minTwist),
        maxTwist_(maxTwist)
      {}
 
      virtual ~TwistMapperStrategy () = default;
 
      virtual MatrixType buildTransformationMatrix ( int twist ) const = 0;
 
      int minTwist() const { return minTwist_; }
      int maxTwist() const { return maxTwist_; }
 
    private:
      int minTwist_;
      int maxTwist_;
    };
 
 
 
    // TwistMapperCreator
    // ------------------
 
    template <class ct, int dim>
    class TwistMapperCreator
    {
      typedef CachingTraits<ct, dim> Traits;
    public:
      typedef typename Traits::QuadratureType QuadratureType;
      typedef typename Traits::PointType PointType;
      typedef typename Traits::MapperType MapperType;
      typedef FieldVector<ct, dim+1> CoordinateType;
      typedef TwistStorage<ct, dim> TwistStorageType;
 
    public:
      TwistMapperCreator(const QuadratureType& quad);
 
      ~TwistMapperCreator();
 
      const TwistStorageType* createStorage() const;
 
      int minTwist() const {
        return helper_->minTwist();
      }
 
      int maxTwist() const {
        return helper_->maxTwist();
      }
 
    private:
      typedef typename TwistMapperStrategy<ct, dim>::MatrixType MatrixType;
 
    private:
      TwistMapperCreator(const TwistMapperCreator&);
      TwistMapperCreator& operator=(const TwistMapperCreator&);
 
    private:
      const QuadratureType& quad_;
      TwistMapperStrategy<ct, dim>* helper_;
 
      static const ct eps_;
    };
 
 
 
    // PointTwistMapperStrategy
    // ------------------------
 
    template <class ct, int dim>
    class PointTwistMapperStrategy final
      : public TwistMapperStrategy<ct, dim>
    {
      typedef TwistMapperStrategy<ct, dim> BaseType;
      typedef typename BaseType::MatrixType MatrixType;
 
    public:
      PointTwistMapperStrategy(GeometryType geo)
        : BaseType( 0, 1 )
      {
        assert(dim == 0);
      }
 
      virtual MatrixType buildTransformationMatrix ( int twist ) const
      {
        const auto &refElement = Dune::ReferenceElements< ct, dim >::cube();
 
        assert( twist == 0 );
        MatrixType mat;
        mat[ 0 ] = refElement.position( 0, dim );
        return mat;
      }
    };
 
 
 
    // LineTwistMapperStrategy
    // -----------------------
 
    template <class ct, int dim>
    class LineTwistMapperStrategy final
      : public TwistMapperStrategy<ct, dim>
    {
      typedef TwistMapperStrategy<ct, dim> BaseType;
      typedef typename BaseType::MatrixType MatrixType;
 
    public:
      LineTwistMapperStrategy(GeometryType geo)
        : BaseType( 0, 2 )
      {
        assert(dim == 1);
      }
 
      MatrixType buildTransformationMatrix ( int twist ) const override
      {
        const auto &refElement = Dune::ReferenceElements< ct, dim >::cube();
 
        assert( (twist == 0) || (twist == 1) );
        MatrixType mat;
        mat[ twist   ] = refElement.position( 0, dim );
        mat[ 1-twist ] = refElement.position( 1, dim );
        return mat;
      }
    };
 
 
 
    // TriangleTwistMapperStrategy
    // ---------------------------
 
    template <class ct, int dim>
    class TriangleTwistMapperStrategy final
      : public TwistMapperStrategy<ct, dim>
    {
      typedef TwistMapperStrategy<ct, dim> BaseType;
      typedef typename BaseType::MatrixType MatrixType;
 
    public:
      TriangleTwistMapperStrategy(GeometryType geo)
        : BaseType( -3, 3 )
      {
        assert(dim == 2);
      }
 
      virtual MatrixType buildTransformationMatrix ( int twist ) const override
      {
        typedef Dune::Fem::FaceTopologyMapping<tetra> FaceTopo;
 
        const auto &refElement = Dune::ReferenceElements< ct, dim >::simplex();
 
        MatrixType mat( ct( 0 ) );
        for (int idx = 0; idx < dim+1; ++idx)
          mat[idx] = refElement.position( FaceTopo::twistedDuneIndex(idx, twist), dim); // dim == codim here
        return mat;
      }
    };
 
 
 
    // QuadrilateralTwistMapperStrategy
    // --------------------------------
 
    template <class ct, int dim>
    class QuadrilateralTwistMapperStrategy final
      : public TwistMapperStrategy<ct, dim>
    {
      typedef TwistMapperStrategy<ct, dim> BaseType;
      typedef typename BaseType::MatrixType MatrixType;
 
    public:
      QuadrilateralTwistMapperStrategy(GeometryType geo)
        : BaseType( -4, 4 )
      {
        assert(dim == 2);
      }
 
      MatrixType buildTransformationMatrix ( int twist ) const override
      {
        typedef Dune::Fem::FaceTopologyMapping<hexa> FaceTopo;
 
        const auto &refElement = Dune::ReferenceElements< ct, dim >::cube();
 
        MatrixType mat( ct( 0 ) );
        for (int idx = 0; idx < dim+1; ++idx)
          mat[idx] = refElement.position( FaceTopo::twistedDuneIndex(idx, twist), dim); // dim == codim here
        return mat;
      }
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#include "twistprovider.cc"
#endif // #ifndef DUNE_FEM_TWISTPROVIDER_HH