agrid.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ALBERTAGRID_IMP_HH
#define DUNE_ALBERTAGRID_IMP_HH

#if HAVE_ALBERTA || DOXYGEN

#include <cassert>
#include <cstddef>

#include <algorithm>
#include <iostream>
#include <fstream>
#include <memory>
#include <vector>

// Dune includes
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/stdstreams.hh>
#include <dune/common/parallel/communication.hh>

#include <dune/grid/common/grid.hh>
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/common/adaptcallback.hh>
#include <dune/grid/common/sizecache.hh>

//- Local includes
// some cpp defines and include of alberta.h
#include "albertaheader.hh"

#include <dune/grid/albertagrid/misc.hh>
#include <dune/grid/albertagrid/capabilities.hh>
#include <dune/grid/albertagrid/backuprestore.hh>

#include <dune/grid/albertagrid/coordcache.hh>
#include <dune/grid/albertagrid/gridfamily.hh>
#include <dune/grid/albertagrid/level.hh>
#include <dune/grid/albertagrid/intersection.hh>
#include <dune/grid/albertagrid/intersectioniterator.hh>
#include <dune/grid/albertagrid/datahandle.hh>
#include <dune/grid/albertagrid/entityseed.hh>

#include "indexsets.hh"
#include "geometry.hh"
#include "entity.hh"
#include "hierarchiciterator.hh"
#include "treeiterator.hh"
#include "leveliterator.hh"
#include "leafiterator.hh"

namespace Dune
{

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

  template< class Grid >
  struct DGFGridFactory;



  // AlbertaGrid
  // -----------

  template< int dim, int dimworld = Alberta::dimWorld >
  class AlbertaGrid
    : public GridDefaultImplementation
      < dim, dimworld, Alberta::Real, AlbertaGridFamily< dim, dimworld > >
  {
    typedef AlbertaGrid< dim, dimworld > This;
    typedef GridDefaultImplementation
    < dim, dimworld, Alberta::Real, AlbertaGridFamily< dim, dimworld > >
    Base;

    template< int, int, class > friend class AlbertaGridEntity;
    template< class > friend class AlbertaLevelGridView;
    template< class > friend class AlbertaLeafGridView;
    template< int, class, bool > friend class AlbertaGridTreeIterator;
    template< class > friend class AlbertaGridHierarchicIterator;

    friend class GridFactory< This >;
    friend struct DGFGridFactory< This >;

    friend class AlbertaGridIntersectionBase< const This >;
    friend class AlbertaGridLeafIntersection< const This >;

    friend class AlbertaMarkerVector< dim, dimworld >;
#if (__GNUC__ < 4) && !(defined __ICC)
    // add additional friend decls for gcc 3.4
    friend struct AlbertaMarkerVector< dim, dimworld >::MarkSubEntities<true>;
    friend struct AlbertaMarkerVector< dim, dimworld >::MarkSubEntities<false>;
#endif
    friend class AlbertaGridIndexSet< dim, dimworld >;
    friend class AlbertaGridHierarchicIndexSet< dim, dimworld >;

    template< class, class >
    friend class Alberta::AdaptRestrictProlongHandler;

  public:
    typedef AlbertaGridFamily< dim, dimworld > GridFamily;

    typedef typename GridFamily::ctype ctype;

    static const int dimension = GridFamily::dimension;
    static const int dimensionworld = GridFamily::dimensionworld;

    // the Traits
    typedef typename AlbertaGridFamily< dim, dimworld >::Traits Traits;

    typedef typename Traits::LeafIndexSet LeafIndexSet;
    typedef typename Traits::LevelIndexSet LevelIndexSet;

    typedef typename Traits::HierarchicIndexSet HierarchicIndexSet;

    typedef typename Traits::GlobalIdSet GlobalIdSet;
    typedef typename Traits::LocalIdSet LocalIdSet;

    typedef typename Traits::Communication Communication;

  private:
    typedef typename Traits::template Codim<0>::LeafIterator LeafIterator;

    typedef AlbertaGridIdSet<dim,dimworld> IdSetImp;

    struct AdaptationState
    {
      enum Phase { ComputationPhase, PreAdaptationPhase, PostAdaptationPhase };

    private:
      Phase phase_;
      int coarsenMarked_;
      int refineMarked_;

    public:
      AdaptationState ()
        : phase_( ComputationPhase ),
          coarsenMarked_( 0 ),
          refineMarked_( 0 )
      {}

      void mark ( int count )
      {
        if( count < 0 )
          ++coarsenMarked_;
        if( count > 0 )
          refineMarked_ += (2 << count);
      }

      void unmark ( int count )
      {
        if( count < 0 )
          --coarsenMarked_;
        if( count > 0 )
          refineMarked_ -= (2 << count);
      }

      bool coarsen () const
      {
        return (coarsenMarked_ > 0);
      }

      int refineMarked () const
      {
        return refineMarked_;
      }

      void preAdapt ()
      {
        if( phase_ != ComputationPhase )
          error( "preAdapt may only be called in computation phase." );
        phase_ = PreAdaptationPhase;
      }

      void adapt ()
      {
        if( phase_ != PreAdaptationPhase )
          error( "adapt may only be called in preadapdation phase." );
        phase_ = PostAdaptationPhase;
      }

      void postAdapt ()
      {
        if( phase_ != PostAdaptationPhase )
          error( "postAdapt may only be called in postadaptation phase." );
        phase_ = ComputationPhase;

        coarsenMarked_ = 0;
        refineMarked_ = 0;
      }

    private:
      void error ( const std::string &message )
      {
        DUNE_THROW( InvalidStateException, message );
      }
    };

    template< class DataHandler >
    struct AdaptationCallback;

    // max number of allowed levels is 64
    static const int MAXL = 64;

    typedef Alberta::ElementInfo< dimension > ElementInfo;
    typedef Alberta::MeshPointer< dimension > MeshPointer;
    typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
    typedef AlbertaGridLevelProvider< dimension > LevelProvider;

  public:
    AlbertaGrid ( const This & ) = delete;
    This &operator= ( const This & ) = delete;

    AlbertaGrid ();

    AlbertaGrid ( const Alberta::MacroData< dimension > &macroData,
                  const std::shared_ptr< DuneBoundaryProjection< dimensionworld > > &projection
                    = std::shared_ptr< DuneBoundaryProjection< dimensionworld > >() );

    template< class Proj, class Impl >
    AlbertaGrid ( const Alberta::MacroData< dimension > &macroData,
                  const Alberta::ProjectionFactoryInterface< Proj, Impl > &projectionFactory );

    AlbertaGrid ( const std::string &macroGridFileName );

    ~AlbertaGrid ();

    int maxLevel () const;

    template<int cd, PartitionIteratorType pitype>
    typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator
    lbegin (int level) const;

    template<int cd, PartitionIteratorType pitype>
    typename Traits::template Codim<cd>::template Partition<pitype>::LevelIterator
    lend (int level) const;

    template< int codim >
    typename Traits::template Codim< codim >::LevelIterator
    lbegin ( int level ) const;

    template< int codim >
    typename Traits::template Codim< codim >::LevelIterator
    lend ( int level ) const;

    template< int codim, PartitionIteratorType pitype >
    typename Traits
    ::template Codim< codim >::template Partition< pitype >::LeafIterator
    leafbegin () const;

    template< int codim, PartitionIteratorType pitype >
    typename Traits
    ::template Codim< codim >::template Partition< pitype >::LeafIterator
    leafend () const;

    template< int codim >
    typename Traits::template Codim< codim >::LeafIterator
    leafbegin () const;

    template< int codim >
    typename Traits::template Codim< codim >::LeafIterator
    leafend () const;

    int size (int level, int codim) const;

    int size (int level, GeometryType type) const;

    int size (int codim) const;

    int size (GeometryType type) const;

    std::size_t numBoundarySegments () const
    {
      return numBoundarySegments_;
    }

    typename Traits::LevelGridView levelGridView ( int level ) const
    {
      typedef typename Traits::LevelGridView View;
      typedef typename View::GridViewImp ViewImp;
      return View( ViewImp( *this, level ) );
    }

    typename Traits::LeafGridView leafGridView () const
    {
      typedef typename Traits::LeafGridView View;
      typedef typename View::GridViewImp ViewImp;
      return View( ViewImp( *this ) );
    }

  public:
    //***************************************************************
    //  Interface for Adaptation
    //***************************************************************
    using Base::getMark;
    using Base::mark;

    int getMark ( const typename Traits::template Codim< 0 >::Entity &e ) const;

    bool mark ( int refCount, const typename Traits::template Codim< 0 >::Entity &e );

    void globalRefine ( int refCount );

    template< class DataHandle >
    void globalRefine ( int refCount, AdaptDataHandleInterface< This, DataHandle > &handle );

    bool adapt ();

    template< class DataHandle >
    bool adapt ( AdaptDataHandleInterface< This, DataHandle > &handle );

    bool preAdapt ();

    void postAdapt();

    const Communication &comm () const
    {
      return comm_;
    }

    static std::string typeName ()
    {
      std::ostringstream s;
      s << "AlbertaGrid< " << dim << ", " << dimworld << " >";
      return s.str();
    }

    template< class EntitySeed >
    typename Traits::template Codim< EntitySeed::codimension >::Entity
    entity ( const EntitySeed &seed ) const
    {
      typedef typename Traits::template Codim< EntitySeed::codimension >::EntityImpl EntityImpl;
      return EntityImpl( *this, seed.impl().elementInfo( meshPointer() ), seed.impl().subEntity() );
    }

    //**********************************************************
    // End of Interface Methods
    //**********************************************************
    bool writeGrid( const std::string &filename, ctype time ) const;

    bool readGrid( const std::string &filename, ctype &time );

    // return hierarchic index set
    const HierarchicIndexSet & hierarchicIndexSet () const { return hIndexSet_; }

    const typename Traits :: LevelIndexSet & levelIndexSet (int level) const;

    const typename Traits :: LeafIndexSet & leafIndexSet () const;

    const GlobalIdSet &globalIdSet () const
    {
      return idSet_;
    }

    const LocalIdSet &localIdSet () const
    {
      return idSet_;
    }

    // access to mesh pointer, needed by some methods
    ALBERTA MESH* getMesh () const
    {
      return mesh_;
    };

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

    const DofNumbering &dofNumbering () const
    {
      return dofNumbering_;
    }

    const LevelProvider &levelProvider () const
    {
      return levelProvider_;
    }

    int dune2alberta ( int codim, int i ) const
    {
      return numberingMap_.dune2alberta( codim, i );
    }

    int alberta2dune ( int codim, int i ) const
    {
      return numberingMap_.alberta2dune( codim, i );
    }

    int generic2alberta ( int codim, int i ) const
    {
      return genericNumberingMap_.dune2alberta( codim, i );
    }

    int alberta2generic ( int codim, int i ) const
    {
      return genericNumberingMap_.alberta2dune( codim, i );
    }

  private:
    typedef std::vector<int> ArrayType;

    void setup ();

    // make the calculation of indexOnLevel and so on.
    // extra method because of Reihenfolge
    void calcExtras();

  private:
    // delete mesh and all vectors
    void removeMesh();

    //***********************************************************************
    //  MemoryManagement for Entities and Geometries
    //**********************************************************************
    typedef MakeableInterfaceObject< typename Traits::template Codim< 0 >::Entity >
    EntityObject;

  public:
    friend class AlbertaGridLeafIntersectionIterator< const This >;

    template< int codim >
    static int
    getTwist ( const typename Traits::template Codim< codim >::Entity &entity )
    {
      return entity.impl().twist();
    }

    template< int codim >
    static int
    getTwist ( const typename Traits::template Codim< 0 >::Entity &entity, int subEntity )
    {
      return entity.impl().template twist< codim >( subEntity );
    }

    static int
    getTwistInInside ( const typename Traits::LeafIntersection &intersection )
    {
      return intersection.impl().twistInInside();
    }

    static int
    getTwistInOutside ( const typename Traits::LeafIntersection &intersection )
    {
      return intersection.impl().twistInOutside();
    }

  public:
    // read global element number from elNumbers_
    const Alberta::GlobalVector &
    getCoord ( const ElementInfo &elementInfo, int vertex ) const;

  private:
    // pointer to an Albert Mesh, which contains the data
    MeshPointer mesh_;

    // communication
    Communication comm_;

    // maximum level of the mesh
    int maxlevel_;

    // number of boundary segments within the macro grid
    size_t numBoundarySegments_;

    // map between ALBERTA and DUNE numbering
    Alberta::NumberingMap< dimension, Alberta::Dune2AlbertaNumbering > numberingMap_;
    Alberta::NumberingMap< dimension, Alberta::Generic2AlbertaNumbering > genericNumberingMap_;

    DofNumbering dofNumbering_;

    LevelProvider levelProvider_;

    // hierarchical numbering of AlbertaGrid, unique per codim
    HierarchicIndexSet hIndexSet_;

    // the id set of this grid
    IdSetImp idSet_;

    // the level index set, is generated from the HierarchicIndexSet
    // is generated, when accessed
    mutable std::vector< typename GridFamily::LevelIndexSetImp * > levelIndexVec_;

    // the leaf index set, is generated from the HierarchicIndexSet
    // is generated, when accessed
    mutable typename GridFamily::LeafIndexSetImp* leafIndexSet_;

    SizeCache< This > sizeCache_;

    typedef AlbertaMarkerVector< dim, dimworld > MarkerVector;

    // needed for VertexIterator, mark on which element a vertex is treated
    mutable MarkerVector leafMarkerVector_;

    // needed for VertexIterator, mark on which element a vertex is treated
    mutable std::vector< MarkerVector > levelMarkerVector_;

#if DUNE_ALBERTA_CACHE_COORDINATES
    Alberta::CoordCache< dimension > coordCache_;
#endif

    // current state of adaptation
    AdaptationState adaptationState_;
  };

} // namespace Dune

#include "albertagrid.cc"

// undef all dangerous defines
#undef DIM
#undef DIM_OF_WORLD

#ifdef _ABS_NOT_DEFINED_
#undef ABS
#endif

#ifdef _MIN_NOT_DEFINED_
#undef MIN
#endif

#ifdef _MAX_NOT_DEFINED_
#undef MAX
#endif

#ifdef obstack_chunk_alloc
#undef obstack_chunk_alloc
#endif
#ifdef obstack_chunk_free
#undef obstack_chunk_free
#endif
#include <dune/grid/albertagrid/undefine-3.0.hh>

// We use MEM_ALLOC, so undefine it here.
#undef MEM_ALLOC

// We use MEM_REALLOC, so undefine it here.
#undef MEM_REALLOC

// We use MEM_CALLOC, so undefine it here.
#undef MEM_CALLOC

// We use MEM_FREE, so undefine it here.
#undef MEM_FREE

// Macro ERROR may be defined by alberta_util.h. If so, undefine it.
#ifdef ERROR
#undef ERROR
#endif // #ifdef ERROR

// Macro ERROR_EXIT may be defined by alberta_util.h. If so, undefine it.
#ifdef ERROR_EXIT
#undef ERROR_EXIT
#endif // #ifdef ERROR_EXIT

// Macro WARNING may be defined by alberta_util.h. If so, undefine it.
#ifdef WARNING
#undef WARNING
#endif // #ifdef WARNING

// Macro TEST may be defined by alberta_util.h. If so, undefine it.
#ifdef TEST
#undef TEST
#endif // #ifdef TEST

// Macro TEST_EXIT may be defined by alberta_util.h. If so, undefine it.
#ifdef TEST_EXIT
#undef TEST_EXIT
#endif // #ifdef TEST_EXIT

// Macro DEBUG_TEST may be defined by alberta_util.h. If so, undefine it.
#ifdef DEBUG_TEST
#undef DEBUG_TEST
#endif // #ifdef DEBUG_TEST

// Macro DEBUG_TEST_EXIT may be defined by alberta_util.h. If so, undefine it.
#ifdef DEBUG_TEST_EXIT
#undef DEBUG_TEST_EXIT
#endif // #ifdef DEBUG_TEST_EXIT

// Macro INFO may be defined by alberta_util.h. If so, undefine it.
#ifdef INFO
#undef INFO
#endif // #ifdef INFO

// Macro PRINT_INFO may be defined by alberta_util.h. If so, undefine it.
#ifdef PRINT_INFO
#undef PRINT_INFO
#endif // #ifdef PRINT_INFO

// Macro PRINT_INT_VEC may be defined by alberta_util.h. If so, undefine it.
#ifdef PRINT_INT_VEC
#undef PRINT_INT_VEC
#endif // #ifdef PRINT_INT_VEC

// Macro PRINT_REAL_VEC may be defined by alberta_util.h. If so, undefine it.
#ifdef PRINT_REAL_VEC
#undef PRINT_REAL_VEC
#endif // #ifdef PRINT_REAL_VEC

// Macro WAIT may be defined by alberta_util.h. If so, undefine it.
#ifdef WAIT
#undef WAIT
#endif // #ifdef WAIT

// Macro WAIT_REALLY may be defined by alberta_util.h. If so, undefine it.
#ifdef WAIT_REALLY
#undef WAIT_REALLY
#endif // #ifdef WAIT_REALLY

// Macro GET_WORKSPACE may be defined by alberta_util.h. If so, undefine it.
#ifdef GET_WORKSPACE
#undef GET_WORKSPACE
#endif // #ifdef GET_WORKSPACE

// Macro FREE_WORKSPACE may be defined by alberta_util.h. If so, undefine it.
#ifdef FREE_WORKSPACE
#undef FREE_WORKSPACE
#endif // #ifdef FREE_WORKSPACE

// Macro MAT_ALLOC may be defined by alberta_util.h. If so, undefine it.
#ifdef MAT_ALLOC
#undef MAT_ALLOC
#endif // #ifdef MAT_ALLOC

// Macro MAT_FREE may be defined by alberta_util.h. If so, undefine it.
#ifdef MAT_FREE
#undef MAT_FREE
#endif // #ifdef MAT_FREE

// Macro NAME may be defined by alberta_util.h. If so, undefine it.
#ifdef NAME
#undef NAME
#endif // #ifdef NAME

// Macro GET_STRUCT may be defined by alberta_util.h. If so, undefine it.
#ifdef GET_STRUCT
#undef GET_STRUCT
#endif // #ifdef GET_STRUCT

// Macro ADD_PARAMETER may be defined by alberta_util.h. If so, undefine it.
#ifdef ADD_PARAMETER
#undef ADD_PARAMETER
#endif // #ifdef ADD_PARAMETER

// Macro GET_PARAMETER may be defined by alberta_util.h. If so, undefine it.
#ifdef GET_PARAMETER
#undef GET_PARAMETER
#endif // #ifdef GET_PARAMETER

#define _ALBERTA_H_

#endif // HAVE_ALBERTA || DOXYGEN

#endif