blockmapper.hh

#ifndef DUNE_FEM_HPDG_SPACE_DISCONTINUOUSGALERKIN_BLOCKMAPPER_HH
#define DUNE_FEM_HPDG_SPACE_DISCONTINUOUSGALERKIN_BLOCKMAPPER_HH
 
#include <cassert>
#include <cstddef>
 
#include <algorithm>
#include <functional>
#include <limits>
#include <type_traits>
#include <utility>
#include <vector>
 
#include <dune/common/exceptions.hh>
 
#include <dune/geometry/dimension.hh>
 
#include <dune/grid/common/gridenums.hh>
#include <dune/grid/utility/persistentcontainer.hh>
 
#include <dune/fem/gridpart/common/gridpart.hh>
#include <dune/fem/io/streams/streams.hh>
#include <dune/fem/space/common/dofmanager.hh>
#include <dune/fem/space/mapper/dofmapper.hh>
 
#include <dune/fem/space/common/capabilities.hh>
 
#include "datahandle.hh"
#include "localdofstorage.hh"
 
namespace Dune
{
 
  namespace Fem
  {
 
    namespace hpDG
    {
 
      // Internal forward declaration
      // ----------------------------
 
      template< class GridPart, class LocalKeys >
      struct DiscontinuousGalerkinBlockMapper;
 
 
 
#ifndef DOXYGEN
 
      // DiscontinuousGalerkinBlockMapperTraits
      // --------------------------------------
 
      template< class GridPart, class LocalKeys >
      struct DiscontinuousGalerkinBlockMapperTraits
      {
        using DofMapperType = DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys >;
        using ElementType = typename GridPart::template Codim< 0 >::EntityType;
        using SizeType = std::size_t;
      };
 
#endif // #ifndef DOXYGEN
 
 
 
      // DiscontinuousGalerkinBlockMapper
      // --------------------------------
 
      template< class GridPart, class LocalKeys >
      struct DiscontinuousGalerkinBlockMapper final
      : public AdaptiveDofMapper< DiscontinuousGalerkinBlockMapperTraits< GridPart, LocalKeys > >
      {
        using ThisType = DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys >;
        using BaseType = AdaptiveDofMapper< DiscontinuousGalerkinBlockMapperTraits< GridPart, LocalKeys > >;
 
      public:
        using SizeType = typename BaseType::SizeType;
        using GlobalKeyType = typename BaseType::GlobalKeyType;
 
        using GridPartType = GridPart;
        using ElementType = typename BaseType::ElementType;
 
        using LocalKeysType = LocalKeys;
        using KeyType = typename LocalKeysType::KeyType;
 
      private:
        using DataHandleType = DataHandle< ThisType >;
        friend class DataHandle< DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys > >;
 
        struct Reserve;
        struct Resize;
 
        using LocalDofStorageType = LocalDofStorage< GlobalKeyType >;
 
        using GridType = typename GridPartType::GridType;
        using GridElementType = typename GridType::template Codim< 0 >::Entity;
 
      public:
        using DofManagerType = DofManager< GridType >;
 
        template< class Function >
        DiscontinuousGalerkinBlockMapper ( const GridPartType &gridPart,
                                           const LocalKeysType &localKeys,
                                           const KeyType &value, Function function )
          : gridPart_( gridPart ),
            localKeys_( localKeys ),
            key_( value ),
            keys_( gridPart_.get().grid(), 0, std::make_pair( key_, key_ ) ),
            size_( 0u ),
            dofs_( gridPart_.get().grid(), 0 ),
            dofManager_( DofManagerType::instance( gridPart.grid() ) )
        {
          auto first = gridPart.template begin< 0, All_Partition >();
          auto last = gridPart.template end< 0, All_Partition >();
          for( ; first != last; ++first )
          {
            const ElementType &element = *first;
            const KeyType key = function( element );
 
            const GridElementType &gridElement = gridEntity( element );
            keys_[ gridElement ] = std::make_pair( key, key );
 
            auto &dofs = dofs_[ gridElement ];
            resize( dofs, localKeys.blocks( gridElement.type(), key ) );
          }
 
          dofManager_.get().addIndexSet( *this );
        }
 
        DiscontinuousGalerkinBlockMapper ( const GridPartType &gridPart,
                                           const LocalKeysType &localKeys,
                                           const KeyType &value )
          : DiscontinuousGalerkinBlockMapper( gridPart, localKeys, value, [&value]( const ElementType & ){ return value; } )
        {}
 
        DiscontinuousGalerkinBlockMapper ( const ThisType & ) = delete;
 
        DiscontinuousGalerkinBlockMapper ( ThisType && ) = default;
 
        ThisType &operator= ( const ThisType & ) = delete;
 
        ThisType &operator= ( ThisType && ) = default;
 
#ifndef DOXYGEN
 
        ~DiscontinuousGalerkinBlockMapper ()
        {
          dofManager_.get().removeIndexSet( *this );
        }
 
#endif // #ifndef DOXYGEN
 
        SizeType size () const { return size_; }
 
        int maxNumDofs () const { return localKeys().maxBlocks(); }
 
        SizeType numDofs ( const ElementType &element ) const
        {
          return numDofs( element, Codim< ElementType::codimension >() );
        }
 
        template< class Entity >
        SizeType numEntityDofs ( const Entity &entity ) const
        {
          return numDofs( entity, Codim< Entity::codimension >() );
        }
 
        [[deprecated("Use onSubEntity method with char vector instead")]]
        void onSubEntity ( const ElementType &element, int i, int c, std::vector< bool > &indices ) const
        {
          std::vector< char > _idx;
          onSubEntity(element, i, c, _idx);
          indices.resize( _idx.size() );
          for (std::size_t i=0; i<_idx.size();++i)
            _idx[i] = indices[i] > 0;
        }
        // this method returns which local dofs are attached to the given subentity.
        // indices[locDofNr] =
        //  0 : not attached, not equal to 0 : attached
        // (so this method can still be used in the way the deprecated method was).
        // New: In case the dof can be associated to a component of the
        //      space, the value returned is that component+1. In other
        //      cases (normal velocity for RT for example) the value is -1).
        // So indices[i] is in [-1,dimRange+1]
        void onSubEntity ( const ElementType &element, int i, int c, std::vector< char > &indices ) const
        {
          indices.resize( numDofs(element) );
          if (c == 0)
            std::fill(indices.begin(),indices.end(),1);
          else
            std::fill(indices.begin(),indices.end(),0);
        }
 
        static constexpr bool contains ( const int codim ) { return (codim == 0); }
 
        static constexpr bool fixedDataSize ( int codim ) { return (codim != 0); }
 
        template< class Function >
        void mapEach ( const ElementType &element, Function function ) const
        {
          mapEach( element, function, Codim< ElementType::codimension >() );
        }
 
        template< class Entity, class Function >
        void mapEachEntityDof ( const Entity &entity, Function function ) const
        {
          mapEach( entity, function, Codim< Entity::codimension >() );
        }
 
        SizeType numberOfHoles ( const int block ) const
        {
          assert( block == 0 );
          return indices_.size();
        }
 
        GlobalKeyType oldIndex ( const int hole, const int block ) const
        {
          assert( block == 0 );
          GlobalKeyType dof = indices_[ hole ].first;
          assert( dof != std::numeric_limits< GlobalKeyType >::max() );
          return std::move( dof );
        }
 
        GlobalKeyType newIndex ( const int hole, const int block ) const
        {
          assert( block == 0 );
          return indices_[ hole ].second;
        }
 
        static constexpr bool consecutive () { return true; }
 
        static constexpr SizeType oldOffSet ( const int block ) { return 0u; }
 
        static constexpr SizeType offSet ( const int block ) { return 0u; }
 
        static constexpr SizeType numBlocks () { return 1u; }
 
        template< class Element >
        typename std::enable_if< (std::is_same< Element, ElementType >::value || std::is_same< Element, GridElementType >::value), const KeyType & >::type
        key ( const Element &element ) const
        {
          return keys_[ gridEntity( element ) ].first;
        }
 
        void mark ( const KeyType &key, const ElementType &element )
        {
          assert( gridPart().indexSet().contains( element ) );
          keys_[ gridEntity( element ) ].second = key;
        }
 
        KeyType getMark ( const ElementType &element ) const
        {
          assert( gridPart().indexSet().contains( element ) );
          return keys_[ gridEntity( element ) ].second;
        }
 
        template< class Function >
        bool adapt ( Function function );
 
        bool adapt ()
        {
          auto function = []( const ElementType &,
                              const KeyType &, const KeyType &,
                              const std::vector< std::size_t > &,
                              const std::vector< std::size_t > & )
          {};
          return adapt( function );
        }
 
        const DofManagerType &dofManager () const { return dofManager_.get(); }
 
        void insertEntity ( const GridElementType &gridElement )
        {
          resize( dofs_ );
          resize( keys_, std::make_pair( key_, key_ ) );
 
          resize( dofs_[ gridElement ], localKeys().blocks( gridElement.type(), key( gridElement ) ) );
        }
 
        void removeEntity ( const GridElementType &gridElement )
        {
          auto &dofs = dofs_[ gridElement ];
 
          size_ = dofs.reserve( 0u, Reserve( size_ ) );
          Resize function( holes_ );
          for( auto dof : dofs )
            function( dof );
        }
 
        void insertNewEntity ( const GridElementType &gridElement )
        {
          if( gridElement.isNew() )
          {
            resize( dofs_[ gridElement ], localKeys().blocks( gridElement.type(), key( gridElement ) ) );
            assert( gridElement.hasFather() );
            const GridElementType &father = gridElement.father();
            insertNewEntity( father );
            removeEntity( father );
          }
          else
          {
            auto &dofs = dofs_[ gridElement ];
            resize( dofs, localKeys().blocks( gridElement.type(), key( gridElement ) ) );
            for( auto dof : dofs )
            {
              auto iterator = std::lower_bound( holes_.begin(), holes_.end(), dof );
              if( iterator != holes_.end() && *iterator == dof )
                holes_.erase( iterator );
            }
          }
        }
 
        void resize ()
        {
          resize( dofs_ );
          resize( keys_, std::make_pair( key_, key_ ) );
 
          auto first = gridPart().template begin< 0, All_Partition >();
          auto last = gridPart().template end< 0, All_Partition >();
          for( ; first != last; ++first )
          {
            const ElementType &element = *first;
            insertNewEntity( gridEntity( element ) );
          }
#if 0
            auto &dofs = dofs_[ gridEntity( element ) ];
            resize( dofs, localKeys().blocks( element.type(), key( element ) ) );
#endif
        }
 
        bool compress ();
 
        template< class Traits >
        void write ( OutStreamInterface< Traits > & )
        {
          DUNE_THROW( NotImplemented, "Method write() not implemented yet" );
        }
 
        template< class Traits >
        void read ( InStreamInterface< Traits > & )
        {
          DUNE_THROW( NotImplemented, "Method read() not implemented yet" );
        }
 
        void backup () const
        {
          DUNE_THROW( NotImplemented, "Method backup() not implemented" );
        }
 
        void restore ()
        {
          DUNE_THROW( NotImplemented, "Method restore() not implemented" );
        }
 
      private:
        bool compressed () const { return holes_.empty(); }
 
        template< class Element, class Function >
        void mapEach ( const Element &element, Function function, Codim< 0 > ) const
        {
          const auto &dofs = dofs_[ gridEntity( element ) ];
          std::size_t size = dofs.size();
          for( std::size_t i = 0; i < size; ++i )
            function( i, dofs[ i ] );
        }
 
        template< class Entity, class Function, int codim >
        void mapEach ( const Entity &entity, Function function, Codim< codim > ) const
        {}
 
        template< class Element >
        SizeType numDofs ( const Element &element, Codim< 0 > ) const
        {
          return dofs_[ gridEntity( element ) ].size();
        }
 
        template< class Entity, int codim >
        SizeType numDofs ( const Entity &entity, Codim< codim > ) const
        {
          return 0u;
        }
 
        void resize ( LocalDofStorageType &dofs, SizeType n )
        {
          size_ = dofs.reserve( n, Reserve( size_ ) );
          dofs.resize( Resize( holes_ ) );
        }
 
        template< class T >
        static void resize ( PersistentContainer< GridType, T > &container, const T &value = T() )
        {
          container.resize( value );
          container.shrinkToFit();
        }
 
        const GridPartType &gridPart () const { return gridPart_.get(); }
 
        const LocalKeysType &localKeys () const { return localKeys_.get(); }
 
        std::reference_wrapper< const GridPartType > gridPart_;
        std::reference_wrapper< const LocalKeysType > localKeys_;
 
        const KeyType key_;
        PersistentContainer< GridType, std::pair< KeyType, KeyType > > keys_;
 
        SizeType size_;
        PersistentContainer< GridType, LocalDofStorageType > dofs_;
 
        std::vector< GlobalKeyType > holes_;
        std::vector< std::pair< GlobalKeyType, GlobalKeyType > > indices_;
 
        std::reference_wrapper< DofManagerType > dofManager_;
      };
 
 
 
      // DiscontinuousGalerkinBlockMapper::Reserve
      // -----------------------------------------
 
      template< class GridPart, class LocalKeys >
      struct DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys >::Reserve
      {
        explicit Reserve ( SizeType &size ) : size_( size ) {}
 
        GlobalKeyType operator() () { return size_++; }
 
        operator SizeType () const { return size_; }
 
      private:
        SizeType &size_;
      };
 
 
 
      // DiscontinuousGalerkinBlockMapper::Resize
      // ----------------------------------------
 
      template< class GridPart, class LocalKeys >
      struct DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys >::Resize
      {
        explicit Resize ( std::vector< GlobalKeyType > &holes )
          : holes_( holes )
        {
          assert( std::is_sorted( holes_.begin(), holes_.end() ) );
        }
 
        ~Resize ()
        {
          assert( std::is_sorted( holes_.begin(), holes_.end() ) );
        }
 
        void operator() ( const GlobalKeyType &dof )
        {
          auto iterator = std::lower_bound( holes_.begin(), holes_.end(), dof );
          if( iterator == holes_.end() || *iterator != dof )
            holes_.insert( iterator, dof );
        }
 
      private:
        std::vector< GlobalKeyType > &holes_;
      };
 
 
 
      // DiscontinuousGalerkinBlockMapper::adapt
      // ---------------------------------------
 
      template< class GridPart, class LocalKeys >
      template< class Function >
      bool DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys >::adapt ( Function function )
      {
        DataHandleType dataHandle( *this );
        gridPart().communicate( dataHandle, InteriorBorder_All_Interface, ForwardCommunication );
 
        std::vector< std::size_t > origin, destination;
        origin.reserve( maxNumDofs() );
        destination.reserve( maxNumDofs() );
 
        std::size_t count = 0;
 
        auto first = gridPart().template begin< 0, All_Partition >();
        auto last = gridPart().template end< 0, All_Partition >();
        for( ; first != last; ++first )
        {
          const ElementType &element = *first;
          const GridElementType &gridElement = gridEntity( element );
          auto &keys = keys_[ gridElement ];
          if( keys.first == keys.second )
            continue;
 
          // remember old state
          const KeyType prior = keys.first;
          auto &dofs = dofs_[ gridElement ];
          origin.resize( dofs.size() );
          std::copy( dofs.begin(), dofs.end(), origin.begin() );
 
          // reset to new state
          keys.first = keys.second;
          resize( dofs, localKeys().blocks( gridElement.type(), keys.first ) );
 
          // call function
          destination.resize( dofs.size() );
          std::copy( dofs.begin(), dofs.end(), destination.begin() );
          function( element, prior, keys.first, origin, destination );
 
          ++count;
        }
        return (count > 0u);
      }
 
 
 
      // DiscontinuousGalerkinBlockMapper::compress
      // ------------------------------------------
 
      template< class GridPart, class LocalKeys >
      bool DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys >::compress ()
      {
        // resize persistent containers
        resize( dofs_ );
        resize( keys_, std::make_pair( key_, key_ ) );
 
        for( auto &dofs : dofs_ )
          dofs.resize( Resize( holes_ ) );
 
        // check if compress is necessary
        if( holes_.empty() )
        {
          indices_.clear();
          return false;
        }
 
        // get number of dofs after compress
        assert( size_ >= holes_.size() );
        size_ -= holes_.size();
 
        // remove trailing indices
        auto iterator = std::lower_bound( holes_.begin(), holes_.end(), size_ );
        holes_.erase( iterator, holes_.end() );
 
        // initialize vector of old and new indices
        auto assign = []( GlobalKeyType newIndex ) -> std::pair< GlobalKeyType, GlobalKeyType >
        {
          GlobalKeyType oldIndex = std::numeric_limits< GlobalKeyType >::max();
          return std::make_pair( oldIndex, newIndex );
        };
        indices_.resize( holes_.size() );
        std::transform( holes_.begin(), holes_.end(), indices_.begin(), assign );
 
        // update local dof storages
        std::size_t hole = 0u;
        for( auto &dofs : dofs_ )
        {
          for( auto &dof : dofs )
          {
            if( dof >= size_ )
            {
              auto &indices = indices_.at( hole++ );
              indices.first = dof;
              dof = indices.second;
            }
          }
        }
 
        assert( hole == holes_.size() );
        holes_.clear();
 
        // replace all but markers currently used by default value
        PersistentContainer< GridType, std::pair< KeyType, KeyType > > tmp( gridPart().grid(), 0, std::make_pair( key_, key_ ) );
        auto first = gridPart().template begin< 0, All_Partition >();
        auto last = gridPart().template end< 0, All_Partition >();
        for( ; first != last; ++first )
        {
          const ElementType &element = *first;
          const GridElementType &gridElement = gridEntity( element );
          tmp[ gridElement ] = keys_[ gridElement ];
        }
        keys_.swap( tmp );
 
        return true;
      }
 
    } // namespace hpDG
 
    namespace Capabilities
    {
      // isConsecutiveIndexSet
      // ---------------------
 
      template< class GridPart, class LocalKeys >
      struct isConsecutiveIndexSet< hpDG::DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys > >
      {
        static const bool v = true;
      };
 
      // isAdaptiveDofMapper
      // -------------------
 
      template< class GridPart, class LocalKeys >
      struct isAdaptiveDofMapper< hpDG::DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys > >
      {
        static const bool v = true;
      };
 
    } // namespace Capabilities
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_HPDG_SPACE_DISCONTINUOUSGALERKIN_BLOCKMAPPER_HH