libdune-grid-dev 2.3.1-1 / usr / include / dune / grid / albertagrid / dofadmin.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ALBERTA_DOFADMIN_HH
#define DUNE_ALBERTA_DOFADMIN_HH

#include <dune/grid/albertagrid/misc.hh>
#include <dune/grid/albertagrid/elementinfo.hh>

#if HAVE_ALBERTA

namespace Dune
{

  namespace Alberta
  {

    // External Forward Declarations
    // -----------------------------

    template< int dim >
    class MeshPointer;



    // DofAccess
    // ---------

    template< int dim, int codim >
    class DofAccess
    {
      static const int codimtype = CodimType< dim, codim >::value;

    public:
      static const int numSubEntities = NumSubEntities< dim, codim >::value;

      static const int dimension = dim;
      static const int codimension = codim;

      typedef Alberta::ElementInfo< dimension > ElementInfo;

      DofAccess ()
        : node_( -1 )
      {}

      explicit DofAccess ( const DofSpace *dofSpace )
      {
        assert( dofSpace );
        node_ = dofSpace->admin->mesh->node[ codimtype ];
        index_ = dofSpace->admin->n0_dof[ codimtype ];
      }

      int operator() ( const Element *element, int subEntity, int i ) const
      {
        assert( element );
        assert( node_ != -1 );
        assert( subEntity < numSubEntities );
        return element->dof[ node_ + subEntity ][ index_ + i ];
      }

      int operator() ( const Element *element, int subEntity ) const
      {
        return (*this)( element, subEntity, 0 );
      }

      int operator() ( const ElementInfo &elementInfo, int subEntity, int i ) const
      {
        return (*this)( elementInfo.el(), subEntity, i );
      }

      int operator() ( const ElementInfo &elementInfo, int subEntity ) const
      {
        return (*this)( elementInfo.el(), subEntity );
      }

    private:
      int node_;
      int index_;
    };



    // HierarchyDofNumbering
    // ---------------------

    template< int dim >
    class HierarchyDofNumbering
    {
      typedef HierarchyDofNumbering< dim > This;

    public:
      static const int dimension = dim;

      typedef Alberta::MeshPointer< dimension > MeshPointer;
      typedef Alberta::ElementInfo< dimension > ElementInfo;

    private:
      static const int nNodeTypes = N_NODE_TYPES;

      template< int codim >
      struct CreateDofSpace;

      template< int codim >
      struct CacheDofSpace;

      typedef std::pair< int, int > Cache;

    public:
      HierarchyDofNumbering ()
      {}

    private:
      HierarchyDofNumbering ( const This & );
      This &operator= ( const This & );

    public:
      ~HierarchyDofNumbering ()
      {
        release();
      }

      int operator() ( const Element *element, int codim, unsigned int subEntity ) const
      {
        assert( !(*this) == false );
        assert( (codim >= 0) && (codim <= dimension) );
        const Cache &cache = cache_[ codim ];
        return element->dof[ cache.first + subEntity ][ cache.second ];
      }

      int operator() ( const ElementInfo &element, int codim, unsigned int subEntity ) const
      {
        return (*this)( element.el(), codim, subEntity );
      }

      operator bool () const
      {
        return (bool)mesh_;
      }

      const DofSpace *dofSpace ( int codim ) const
      {
        assert( *this );
        assert( (codim >= 0) && (codim <= dimension) );
        return dofSpace_[ codim ];
      }

      const DofSpace *emptyDofSpace () const
      {
        assert( *this );
        return emptySpace_;
      }

      const MeshPointer &mesh () const
      {
        return mesh_;
      }

      int size ( int codim ) const
      {
        return dofSpace( codim )->admin->size;
      }

      void create ( const MeshPointer &mesh );

      void release ()
      {
        if( *this )
        {
          for( int codim = 0; codim <= dimension; ++codim )
            freeDofSpace( dofSpace_[ codim ] );
          freeDofSpace( emptySpace_ );
          mesh_ = MeshPointer();
        }
      }

    private:
      static const DofSpace *createEmptyDofSpace ( const MeshPointer &mesh );
      static const DofSpace *createDofSpace ( const MeshPointer &mesh,
                                              const std::string &name,
                                              const int (&ndof)[ nNodeTypes ],
                                              const bool periodic = false );
      static void freeDofSpace ( const DofSpace *dofSpace );

      MeshPointer mesh_;
      const DofSpace *emptySpace_;
      const DofSpace *dofSpace_[ dimension+1 ];
      Cache cache_[ dimension+1 ];
    };



    template< int dim >
    inline void
    HierarchyDofNumbering< dim >::create ( const MeshPointer &mesh )
    {
      release();

      if( !mesh )
        return;

      mesh_ = mesh;
      ForLoop< CreateDofSpace, 0, dimension >::apply( mesh_, dofSpace_ );
      ForLoop< CacheDofSpace, 0, dimension >::apply( dofSpace_, cache_ );

      emptySpace_ = createEmptyDofSpace( mesh_ );
      for( int i = 0; i < nNodeTypes; ++i )
        assert( emptySpace_->admin->n_dof[ i ] == 0 );
    }



    template< int dim >
    inline const DofSpace *
    HierarchyDofNumbering< dim >::createEmptyDofSpace ( const MeshPointer &mesh )
    {
      int ndof[ nNodeTypes ];
      for( int i = 0; i < nNodeTypes; ++i )
        ndof[ i ] = 0;
      std::string name = "Empty";
      return createDofSpace( mesh, name, ndof );
    }


#if DUNE_ALBERTA_VERSION >= 0x300
    template< int dim >
    inline const DofSpace *
    HierarchyDofNumbering< dim >::createDofSpace ( const MeshPointer &mesh,
                                                   const std::string &name,
                                                   const int (&ndof)[ nNodeTypes ],
                                                   const bool periodic )
    {
      const ALBERTA FLAGS flags
        = ADM_PRESERVE_COARSE_DOFS | (periodic ? ADM_PERIODIC : 0);
      return ALBERTA get_dof_space ( mesh, name.c_str(), ndof, flags );
    }
#endif // #if DUNE_ALBERTA_VERSION >= 0x300

#if DUNE_ALBERTA_VERSION == 0x200
    template< int dim >
    inline const DofSpace *
    HierarchyDofNumbering< dim >::createDofSpace ( const MeshPointer &mesh,
                                                   const std::string &name,
                                                   const int (&ndof)[ nNodeTypes ],
                                                   const bool periodic )
    {
      return ALBERTA get_fe_space( mesh, name.c_str(), ndof, NULL, 1 );
    }
#endif // #if DUNE_ALBERTA_VERSION == 0x200


#if DUNE_ALBERTA_VERSION >= 0x300
    template< int dim >
    inline void
    HierarchyDofNumbering< dim >::freeDofSpace ( const DofSpace *dofSpace )
    {
      ALBERTA free_fe_space( dofSpace );
    }
#endif // #if DUNE_ALBERTA_VERSION >= 0x300

#if DUNE_ALBERTA_VERSION == 0x200
    template< int dim >
    inline void
    HierarchyDofNumbering< dim >::freeDofSpace ( const DofSpace *dofSpace )
    {
      // the const cast is needed due to a bug in ALBERTA 2.0
      ALBERTA free_fe_space( const_cast< DofSpace * >( dofSpace ) );
    }
#endif // #if DUNE_ALBERTA_VERSION == 0x200



    // HierarchyDofNumbering::CreateDofSpace
    // -------------------------------------

    template< int dim >
    template< int codim >
    struct HierarchyDofNumbering< dim >::CreateDofSpace
    {
      static void apply ( const MeshPointer &mesh, const DofSpace *(&dofSpace)[ dim+1 ] )
      {
        int ndof[ nNodeTypes ];
        for( int i = 0; i < nNodeTypes; ++i )
          ndof[ i ] = 0;
        ndof[ CodimType< dim, codim >::value ] = 1;

        std::string name = "Codimension ";
        name += (char)(codim + '0');

        dofSpace[ codim ] = createDofSpace( mesh, name, ndof );
        assert( dofSpace[ codim ] );
      }
    };



    // HierarchyDofNumbering::CacheDofSpace
    // ------------------------------------

    template< int dim >
    template< int codim >
    struct HierarchyDofNumbering< dim >::CacheDofSpace
    {
      static void apply ( const DofSpace *(&dofSpace)[ dim+1 ], Cache (&cache)[ dim+1 ] )
      {
        assert( dofSpace[ codim ] );
        const int codimtype = CodimType< dim, codim >::value;
        cache[ codim ].first = dofSpace[ codim ]->mesh->node[ codimtype ];
        cache[ codim ].second = dofSpace[ codim ]->admin->n0_dof[ codimtype ];
      }
    };

  } // namespace Alberta

} // namespace Dune

#endif // #if HAVE_ALBERTA

#endif // #ifndef DUNE_ALBERTA_DOFADMIN_HH