CpGridData.hpp

//===========================================================================
//
// File: CpGridData.hpp
//
// Created: Sep 17 21:11:41 2013
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//            Bård Skaflestad     <bard.skaflestad@sintef.no>
//            Markus Blatt        <markus@dr-blatt.de>
//            Antonella Ritorto   <antonella.ritorto@opm-op.com>
//
// Comment: Major parts of this file originated in dune/grid/CpGrid.hpp
//          and got transfered here during refactoring for the parallelization.
//
// $Date$
//
// $Revision$
//
//===========================================================================

/*
  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010, 2013, 2022-2023 Equinor ASA.
  Copyright 2013 Dr. Blatt - HPC-Simulation-Software & Services

  This file is part of The Open Porous Media project  (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_CPGRIDDATA_HEADER
#define OPM_CPGRIDDATA_HEADER


#include <dune/common/parallel/mpihelper.hh>
#ifdef HAVE_DUNE_ISTL
#include <dune/istl/owneroverlapcopy.hh>
#endif

#include <dune/common/parallel/communication.hh>
#include <dune/common/parallel/variablesizecommunicator.hh>
#include <dune/grid/common/gridenums.hh>

#if HAVE_OPM_COMMON
#include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/input/eclipse/EclipseState/Grid/NNC.hpp>
#endif

#include <opm/grid/cpgpreprocess/preprocess.h>

#include "Entity2IndexDataHandle.hpp"
#include "CpGridDataTraits.hpp"
//#include "DataHandleWrappers.hpp"
//#include "GlobalIdMapping.hpp"
#include "Geometry.hpp"

#include <array>
#include <initializer_list>
#include <set>
#include <vector>

namespace Opm
{
class EclipseState;
}
namespace Dune
{
class CpGrid;

namespace cpgrid
{

class IndexSet;
class IdSet;
class LevelGlobalIdSet;
class PartitionTypeIndicator;
template<int,int> class Geometry;
template<int> class Entity;
template<int> class EntityRep;
}
}

void refine_and_check(const Dune::cpgrid::Geometry<3, 3>&,
                      const std::array<int, 3>&,
                      bool);

namespace Dune
{
namespace cpgrid
{
namespace mover
{
template<class T, int i> struct Mover;
}

class CpGridData
{
    template<class T, int i> friend struct mover::Mover;
    friend class GlobalIdSet;
    friend class HierarchicIterator;
    friend class Dune::cpgrid::IndexSet;
    friend class Dune::cpgrid::IdSet;
    friend class Dune::cpgrid::LevelGlobalIdSet;

    friend
    void ::refine_and_check(const Dune::cpgrid::Geometry<3, 3>&,
                            const std::array<int, 3>&,
                            bool);

private:
    CpGridData(const CpGridData& g);

public:
    enum{
#ifndef MAX_DATA_COMMUNICATED_PER_ENTITY
        MAX_DATA_PER_CELL = 2
#else
        MAX_DATA_PER_CELL = MAX_DATA_COMMUNICATED_PER_ENTITY
#endif
    };

    CpGridData() = delete;

    explicit CpGridData(MPIHelper::MPICommunicator comm,  std::vector<std::shared_ptr<CpGridData>>& data);



    explicit CpGridData(std::vector<std::shared_ptr<CpGridData>>& data);
    ~CpGridData();




    int size(int codim) const;

    int size (GeometryType type) const
    {
        if (type.isCube()) {
            return size(3 - type.dim());
        } else {
            return 0;
        }
    }

    void readEclipseFormat(const std::string& filename,
                           bool periodic_extension,
                           bool turn_normals = false,
                           bool edge_conformal = false);

#if HAVE_OPM_COMMON
    void processEclipseFormat(const Opm::Deck& deck,
                              bool periodic_extension,
                              bool turn_normals = false,
                              bool clip_z = false,
                              const std::vector<double>& poreVolume = std::vector<double>{},
                              bool edge_conformal = false);

    std::vector<std::size_t>
    processEclipseFormat(const Opm::EclipseGrid* ecl_grid,
                         Opm::EclipseState* ecl_state,
                         bool periodic_extension,
                         bool turn_normals = false,
                         bool clip_z = false,
                         bool pinchActive = true,
                         bool edge_conformal = false);
#endif

    void processEclipseFormat(const grdecl& input_data,
#if HAVE_OPM_COMMON
                              Opm::EclipseState* ecl_state,
#endif
                              std::array<std::set<std::pair<int, int>>, 2>& nnc,
                              bool remove_ij_boundary,
                              bool turn_normals,
                              bool pinchActive,
                              double tolerance_unique_points,
                              bool edge_conformal);

    void getIJK(int c, std::array<int,3>& ijk) const;

    int cellFace(int cell, int local_index) const
    {
        return cell_to_face_[cpgrid::EntityRep<0>(cell, true)][local_index].index();
    }

    auto cellToFace(int cellIdx) const
    {
        return cell_to_face_[cpgrid::EntityRep<0>(cellIdx, true)];
    }

    const auto& cellToPoint() const
    {
        return cell_to_point_;
    }

    const auto& cellToPoint(int cellIdx) const
    {
        return cell_to_point_[cellIdx];
    }

    int faceToCellSize(int face) const {
        Dune::cpgrid::EntityRep<1> faceRep(face, true);
        return face_to_cell_[faceRep].size();
    }

    auto faceTag(int faceIdx) const
    {
        Dune::cpgrid::EntityRep<1> faceRep(faceIdx, true);
        return face_tag_[faceRep];
    }

    auto faceNormals(int faceIdx) const
    {
        Dune::cpgrid::EntityRep<1> faceRep(faceIdx, true);
        return face_normals_[faceRep];
    }

    auto faceToPoint(int faceIdx) const
    {
        return face_to_point_[faceIdx];
    }

    int numFaces() const
    {
        return face_to_cell_.size();
    }

    auto cornerHistorySize() const
    {
        return corner_history_.size();
    }

    const auto& getCornerHistory(int cornerIdx) const
    {
        if(cornerHistorySize()) {
            return corner_history_[cornerIdx];
        }
        else {
            OPM_THROW(std::logic_error, "Vertex has no history record.\n");
        }
    }

    const std::vector<int>& globalCell() const
    {
        return  global_cell_;
    }

    bool hasNNCs(const std::vector<int>& cellIndices) const;

    bool mark(int refCount, const cpgrid::Entity<0>& element, bool throwOnFailure = false);

    int getMark(const cpgrid::Entity<0>& element) const;

    bool preAdapt();

    bool adapt();

    void postAdapt();

private:
    std::array<Dune::FieldVector<double,3>,8> getReferenceRefinedCorners(int idx_in_parent_cell, const std::array<int,3>& cells_per_dim) const;

public:
    int getGridIdx() const {
        // Not the nicest way of checking if "this" points at the leaf grid view of a mixed grid (with coarse and refined cells).
        // 1. When the grid has been refined at least onece, level_data_ptr_ ->size() >1. Therefore, there is a chance of "this" pointing at the leaf grid view.
        // 2. Unfortunately, level_ is default initialized by 0. This implies, in particular, that if someone wants to check the value of
        //    "this->level_" when "this" points at the leaf grid view of a grid that has been refined, this value is - unfortunately - equal to 0.
        // 3. Due to 2. we need an extra bool value to distinguish between the actual level 0 grid and such a leaf grid view (with incorrect level_ == 0). For this
        //    reason we check if child_to_parent_cells_.empty() [true for actual level 0 grid, false for the leaf grid view].
        // --- TO BE IMPROVED ---
        if ((level_data_ptr_ ->size() >1) && (level_ == 0) && (!child_to_parent_cells_.empty())) {
            return level_data_ptr_->size() -1;
        }
        return level_;
    }
    const std::vector<std::shared_ptr<Dune::cpgrid::CpGridData>>& levelData() const
    {
        if (level_data_ptr_->empty()) {
            OPM_THROW(std::logic_error, "Level data has not been initialized\n");
        }
        return *level_data_ptr_;
    }

    const std::tuple<int,std::vector<int>>& getChildrenLevelAndIndexList(int elemIdx) const {
        return parent_to_children_cells_[elemIdx];
    }

    const std::vector<std::tuple<int,std::vector<int>>>& getParentToChildren() const {
        return parent_to_children_cells_;
    }

    const cpgrid::DefaultGeometryPolicy getGeometry() const
    {
        return geometry_;
    }

    int getLeafIdxFromLevelIdx(int level_cell_idx) const
    {
        if (level_to_leaf_cells_.empty()) {
            OPM_THROW(std::logic_error, "Grid has no LGRs. No mapping to the leaf.\n");
        }
        return level_to_leaf_cells_[level_cell_idx];
    }

    std::tuple< const std::shared_ptr<CpGridData>,
                const std::vector<std::array<int,2>>>                  // parent_to_refined_corners(~boundary_old_to_new_corners)
    refineSingleCell(const std::array<int,3>& cells_per_dim,
                     const int& parent_idx,
                     std::vector<std::vector<std::pair<int, std::vector<int>>>>& faceInMarkedElemAndRefinedFaces) const;

    // @breif Compute center of an entity/element/cell in the Eclipse way:
    //        - Average of the 4 corners of the bottom face.
    //        - Average of the 4 corners of the top face.
    //        Return average of the previous computations.
    // @param [in]   int   Index of a cell.
    // @return            'eclipse centroid'
    std::array<double,3> computeEclCentroid(const int idx) const;

    // @breif Compute center of an entity/element/cell in the Eclipse way:
    //        - Average of the 4 corners of the bottom face.
    //        - Average of the 4 corners of the top face.
    //        Return average of the previous computations.
    // @param [in]   Entity<0>   Entity
    // @return                   'eclipse centroid'
    std::array<double,3> computeEclCentroid(const Entity<0>& elem) const;

    // Make unique boundary ids for all intersections.
    void computeUniqueBoundaryIds();

    bool uniqueBoundaryIds() const
    {
        return use_unique_boundary_ids_;
    }

    void setUniqueBoundaryIds(bool uids)
    {
        use_unique_boundary_ids_ = uids;
        if (use_unique_boundary_ids_ && unique_boundary_ids_.empty()) {
            computeUniqueBoundaryIds();
        }
    }

    const std::vector<double>& zcornData() const {
        return zcorn;
    }


    const IndexSet& indexSet() const
    {
        return *index_set_;
    }

    const cpgrid::IdSet& localIdSet() const
    {
        return *local_id_set_;
    }

    const cpgrid::LevelGlobalIdSet& globalIdSet() const
    {
        return *global_id_set_;
    }

    const std::array<int, 3>& logicalCartesianSize() const
    {
        return logical_cartesian_size_;
    }

    void distributeGlobalGrid(CpGrid& grid,
                              const CpGridData& view_data,
                              const std::vector<int>& cell_part);

    template<class DataHandle>
    void communicate(DataHandle& data, InterfaceType iftype, CommunicationDirection dir);

    void computeCellPartitionType();

    void computePointPartitionType();

    void computeCommunicationInterfaces(int noexistingPoints);

    using MPICommunicator = CpGridDataTraits::MPICommunicator ;
    using Communication = CpGridDataTraits::Communication;
    using CollectiveCommunication = CpGridDataTraits::CollectiveCommunication;

    using AttributeSet = CpGridDataTraits::AttributeSet;
#if HAVE_MPI
    using Communicator = CpGridDataTraits::Communicator;

    using InterfaceMap = CpGridDataTraits::InterfaceMap;

    using CommunicationType = CpGridDataTraits::CommunicationType;

    using  ParallelIndexSet = CpGridDataTraits::ParallelIndexSet;

    using RemoteIndices = CpGridDataTraits::RemoteIndices;

    CommunicationType& cellCommunication()
    {
        return cell_comm_;
    }

    const CommunicationType& cellCommunication() const
    {
        return cell_comm_;
    }

    ParallelIndexSet& cellIndexSet()
    {
        return cellCommunication().indexSet();
    }

    const ParallelIndexSet& cellIndexSet() const
    {
        return cellCommunication().indexSet();
    }

    RemoteIndices& cellRemoteIndices()
    {
        return cellCommunication().remoteIndices();
    }

    const RemoteIndices& cellRemoteIndices() const
    {
            return cellCommunication().remoteIndices();
    }
#endif

    const std::vector<int>& sortedNumAquiferCells() const
    {
        return aquifer_cells_;
    }

private:

    void populateGlobalCellIndexSet();

#if HAVE_MPI

    template<class DataHandle>
    void gatherData(DataHandle& data, CpGridData* global_view,
                    CpGridData* distributed_view);


    template<int codim, class DataHandle>
    void gatherCodimData(DataHandle& data, CpGridData* global_data,
                         CpGridData* distributed_data);

    template<class DataHandle>
    void scatterData(DataHandle& data, const CpGridData* global_data,
                     const CpGridData* distributed_data, const InterfaceMap& cell_inf,
                     const InterfaceMap& point_inf);

    template<int codim, class DataHandle>
    void scatterCodimData(DataHandle& data, CpGridData* global_data,
                          CpGridData* distributed_data);

    template<int codim, class DataHandle>
    void communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
                          const Interface& interface);

    template<int codim, class DataHandle>
    void communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
                          const InterfaceMap& interface);

#endif

    void computeGeometry(const CpGrid& grid,
                         const DefaultGeometryPolicy&  globalGeometry,
                         const std::vector<int>& globalAquiferCells,
                         const OrientedEntityTable<0, 1>& globalCell2Faces,
                         DefaultGeometryPolicy& geometry,
                         std::vector<int>& aquiferCells,
                         const OrientedEntityTable<0, 1>& cell2Faces,
                         const std::vector< std::array<int,8> >& cell2Points);

    // Representing the topology
    cpgrid::OrientedEntityTable<0, 1> cell_to_face_;
    cpgrid::OrientedEntityTable<1, 0> face_to_cell_;
    Opm::SparseTable<int>             face_to_point_;
    std::vector< std::array<int,8> >       cell_to_point_;
    std::array<int, 3>                logical_cartesian_size_{};
    std::vector<int>                  global_cell_;
    cpgrid::EntityVariable<enum face_tag, 1> face_tag_;
    cpgrid::DefaultGeometryPolicy geometry_;
    typedef FieldVector<double, 3> PointType;
    cpgrid::SignedEntityVariable<PointType, 1> face_normals_;
    cpgrid::EntityVariable<int, 1> unique_boundary_ids_;
    std::unique_ptr<cpgrid::IndexSet> index_set_;
    std::shared_ptr<const cpgrid::IdSet> local_id_set_;
    std::shared_ptr<LevelGlobalIdSet> global_id_set_;
    std::shared_ptr<PartitionTypeIndicator> partition_type_indicator_;
    std::vector<int> mark_;
    int level_{0};
    std::vector<std::shared_ptr<CpGridData>>* level_data_ptr_;
    // SUITABLE FOR ALL LEVELS EXCEPT FOR LEAFVIEW
    std::vector<int> level_to_leaf_cells_; // In entry 'level cell index', we store 'leafview cell index' // {level LGR, {child0, child1, ...}}
    std::vector<std::tuple<int,std::vector<int>>> parent_to_children_cells_; // {# children in x-direction, ... y-, ... z-}
    std::array<int,3> cells_per_dim_;
    // SUITABLE ONLY FOR LEAFVIEW // {level, cell index in that level}
    std::vector<std::array<int,2>> leaf_to_level_cells_;
    std::vector<std::array<int,2>> corner_history_;
    // SUITABLE FOR ALL LEVELS INCLUDING LEAFVIEW // {level parent cell, parent cell index}
    std::vector<std::array<int,2>> child_to_parent_cells_;
    std::vector<int> cell_to_idxInParentCell_;
    int refinement_max_level_{0};

    Communication ccobj_;

    // Boundary information (optional).
    bool use_unique_boundary_ids_;

    std::vector<double> zcorn;

    std::vector<int> aquifer_cells_;

#if HAVE_MPI

    CommunicationType cell_comm_;

    std::tuple<Interface,Interface,Interface,Interface,Interface> cell_interfaces_;
    /*
    // code deactivated, because users cannot access face indices and therefore
    // communication on faces makes no sense!
    std::tuple<InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap>
    face_interfaces_;
    */
    std::tuple<InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap>
    point_interfaces_;

#endif

    // Return the geometry vector corresponding to the given codim.
    template <int codim>
    const EntityVariable<Geometry<3 - codim, 3>, codim>& geomVector() const
    {
        return geometry_.geomVector<codim>();
    }

    friend class Dune::CpGrid;
    template<int> friend class Entity;
    template<int> friend class EntityRep;
    friend class Intersection;
    friend class PartitionTypeIndicator;
};



#if HAVE_MPI

namespace
{
template<class T>
T& getInterface(InterfaceType iftype,
                std::tuple<T,T,T,T,T>& interfaces)
{
    switch(iftype)
    {
    case 0:
        return std::get<0>(interfaces);
    case 1:
        return std::get<1>(interfaces);
    case 2:
        return std::get<2>(interfaces);
    case 3:
        return std::get<3>(interfaces);
    case 4:
        return std::get<4>(interfaces);
    }
    OPM_THROW(std::runtime_error, "Invalid Interface type was used during communication");
}

} // end unnamed namespace

template<int codim, class DataHandle>
void CpGridData::communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
                                  const Interface& interface)
{
    this->template communicateCodim<codim>(data, dir, interface.interfaces());
}

template<int codim, class DataHandle>
void CpGridData::communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data_wrapper, CommunicationDirection dir,
                                  const InterfaceMap& interface)
{
    Communicator comm(ccobj_, interface);

    if(dir==ForwardCommunication)
        comm.forward(data_wrapper);
    else
        comm.backward(data_wrapper);
}
#endif

template<class DataHandle>
void CpGridData::communicate(DataHandle& data, InterfaceType iftype,
                             CommunicationDirection dir)
{
#if HAVE_MPI
    if(data.contains(3,0))
    {
        Entity2IndexDataHandle<DataHandle, 0> data_wrapper(*this, data);
        communicateCodim<0>(data_wrapper, dir, getInterface(iftype, cell_interfaces_));
    }
    if(data.contains(3,3))
    {
        Entity2IndexDataHandle<DataHandle, 3> data_wrapper(*this, data);
        communicateCodim<3>(data_wrapper, dir, getInterface(iftype, point_interfaces_));
    }
#else
    // Suppress warnings for unused arguments.
    (void) data;
    (void) iftype;
    (void) dir;
#endif
}
}}

#if HAVE_MPI
#include <opm/grid/cpgrid/Entity.hpp>
#include <opm/grid/cpgrid/Indexsets.hpp>

namespace Dune {
namespace cpgrid {

namespace mover
{
template<class T>
class MoveBuffer
{
    friend class Dune::cpgrid::CpGridData;
public:
    void read(T& data)
    {
        data=buffer_[index_++];
    }
    void write(const T& data)
    {
        buffer_[index_++]=data;
    }
    void reset()
    {
        index_=0;
    }
    void resize(std::size_t size)
    {
        buffer_.resize(size);
        index_=0;
    }
private:
    std::vector<T> buffer_;
    typename std::vector<T>::size_type index_;
};
template<class DataHandle,int codim>
struct Mover
{
};

template<class DataHandle>
struct BaseMover
{
    explicit BaseMover(DataHandle& data)
    : data_(data)
    {}
    template<class E>
    void moveData(const E& from, const E& to)
    {
        std::size_t size=data_.size(from);
        buffer.resize(size);
        data_.gather(buffer, from);
        buffer.reset();
        data_.scatter(buffer, to, size);
    }
    DataHandle& data_;
    MoveBuffer<typename DataHandle::DataType> buffer;
};


template<class DataHandle>
struct Mover<DataHandle,0> : public BaseMover<DataHandle>
{
    Mover(DataHandle& data, CpGridData* gatherView,
          CpGridData* scatterView)
    : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
    {}

    void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
    {
        Entity<0> from_entity=Entity<0>(*gatherView_, from_cell_index, true);
        Entity<0> to_entity=Entity<0>(*scatterView_, to_cell_index, true);
        this->moveData(from_entity, to_entity);
    }
    CpGridData* gatherView_;
    CpGridData* scatterView_;
};

template<class DataHandle>
struct Mover<DataHandle,1> : public BaseMover<DataHandle>
{
    Mover(DataHandle& data, CpGridData* gatherView,
          CpGridData* scatterView)
    : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
    {}

    void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
    {
        typedef typename OrientedEntityTable<0,1>::row_type row_type;
        EntityRep<0> from_cell=EntityRep<0>(from_cell_index, true);
        EntityRep<0> to_cell=EntityRep<0>(to_cell_index, true);
        const OrientedEntityTable<0,1>& table = gatherView_->cell_to_face_;
        row_type from_faces=table.operator[](from_cell);
        row_type to_faces=scatterView_->cell_to_face_[to_cell];

        for(int i=0; i<from_faces.size(); ++i)
            this->moveData(from_faces[i], to_faces[i]);
    }
    CpGridData *gatherView_;
    CpGridData *scatterView_;
};

template<class DataHandle>
struct Mover<DataHandle,3> : public BaseMover<DataHandle>
{
    Mover(DataHandle& data, CpGridData* gatherView,
          CpGridData* scatterView)
    : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
    {}
    void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
    {
        const std::array<int,8>& from_cell_points=
            gatherView_->cell_to_point_[from_cell_index];
        const std::array<int,8>& to_cell_points=
            scatterView_->cell_to_point_[to_cell_index];
        for(std::size_t i=0; i<8; ++i)
        {
            this->moveData(Entity<3>(*gatherView_, from_cell_points[i], true),
                           Entity<3>(*scatterView_, to_cell_points[i], true));
        }
    }
    CpGridData* gatherView_;
    CpGridData* scatterView_;
};

} // end mover namespace

template<class DataHandle>
void CpGridData::scatterData(DataHandle& data, const CpGridData* global_data,
                             const CpGridData* distributed_data, const InterfaceMap& cell_inf,
                             const InterfaceMap& point_inf)
{
#if HAVE_MPI
    if(data.contains(3,0))
    {
        Entity2IndexDataHandle<DataHandle, 0> data_wrapper(*global_data, *distributed_data, data);
        communicateCodim<0>(data_wrapper, ForwardCommunication, cell_inf);
    }
    if(data.contains(3,3))
    {
        Entity2IndexDataHandle<DataHandle, 3> data_wrapper(*global_data, *distributed_data, data);
        communicateCodim<3>(data_wrapper, ForwardCommunication, point_inf);
    }
#endif
}

template<int codim, class DataHandle>
void CpGridData::scatterCodimData(DataHandle& data, CpGridData* global_data,
                                  CpGridData* distributed_data)
{
    CpGridData *gather_view, *scatter_view;
    gather_view=global_data;
    scatter_view=distributed_data;

    mover::Mover<DataHandle,codim> mover(data, gather_view, scatter_view);


    for(auto index=distributed_data->cellIndexSet().begin(),
            end = distributed_data->cellIndexSet().end();
        index!=end; ++index)
    {
        std::size_t from=index->global();
        std::size_t to=index->local();
        mover(from,to);
    }
}

namespace
{

template<int codim, class T, class F>
void visitInterior(CpGridData& distributed_data, T begin, T endit, F& func)
{
    for(T it=begin; it!=endit; ++it)
    {
        Entity<codim> entity(distributed_data, it-begin, true);
        PartitionType pt = entity.partitionType();
        if(pt==Dune::InteriorEntity)
        {
            func(*it, entity);
        }
    }
}

template<class DataHandle>
struct GlobalIndexSizeGatherer
{
    GlobalIndexSizeGatherer(DataHandle& data_,
                            std::vector<int>& ownedGlobalIndices_,
                            std::vector<int>& ownedSizes_)
        : data(data_), ownedGlobalIndices(ownedGlobalIndices_), ownedSizes(ownedSizes_)
    {}

    template<class T, class E>
    void operator()(T& i, E& entity)
    {
            ownedGlobalIndices.push_back(i);
            ownedSizes.push_back(data.size(entity));
    }
    DataHandle& data;
    std::vector<int>& ownedGlobalIndices;
    std::vector<int>& ownedSizes;
};

template<class DataHandle>
struct DataGatherer
{
    DataGatherer(mover::MoveBuffer<typename DataHandle::DataType>& buffer_,
                 DataHandle& data_)
        : buffer(buffer_), data(data_)
    {}

    template<class T, class E>
    void operator()(T& /* it */, E& entity)
    {
        data.gather(buffer, entity);
    }
    mover::MoveBuffer<typename DataHandle::DataType>& buffer;
    DataHandle& data;
};

}

template<class DataHandle>
void CpGridData::gatherData(DataHandle& data, CpGridData* global_data,
                            CpGridData* distributed_data)
{
#if HAVE_MPI
    if(data.contains(3,0))
       gatherCodimData<0>(data, global_data, distributed_data);
    if(data.contains(3,3))
       gatherCodimData<3>(data, global_data, distributed_data);
#endif
}

template<int codim, class DataHandle>
void CpGridData::gatherCodimData(DataHandle& data, CpGridData* global_data,
                                 CpGridData* distributed_data)
{
#if HAVE_MPI
    // Get the mapping to global index from  the global id set
    const std::vector<int>& mapping =
        distributed_data->global_id_set_->getMapping<codim>();

    // Get the global indices and data size for the entities whose data is
    // to be sent, i.e. the ones that we own.
    std::vector<int>         owned_global_indices;
    std::vector<int> owned_sizes;
    owned_global_indices.reserve(mapping.size());
    owned_sizes.reserve(mapping.size());

    GlobalIndexSizeGatherer<DataHandle> gisg(data, owned_global_indices, owned_sizes);
    visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gisg);

    // communicate the number of indices that each processor sends
    int no_indices=owned_sizes.size();
    // We will take the address of the first elemet for MPI_Allgather below.
    // Make sure the containers have such an element.
    if ( owned_global_indices.empty() )
        owned_global_indices.resize(1);
    if ( owned_sizes.empty() )
        owned_sizes.resize(1);
    std::vector<int> no_indices_to_recv(distributed_data->ccobj_.size());
    distributed_data->ccobj_.allgather(&no_indices, 1, &(no_indices_to_recv[0]));
    // compute size of the vector capable for receiving all indices
    // and allgather the global indices and the sizes.
    // calculate displacements
    std::vector<int> displ(distributed_data->ccobj_.size()+1, 0);
    std::transform(displ.begin(), displ.end()-1, no_indices_to_recv.begin(), displ.begin()+1,
                   std::plus<int>());
    int global_size=displ[displ.size()-1];//+no_indices_to_recv[displ.size()-1];
    std::vector<int>         global_indices(global_size);
    std::vector<int> global_sizes(global_size);
    MPI_Allgatherv(&(owned_global_indices[0]), no_indices, MPITraits<int>::getType(),
                   &(global_indices[0]), &(no_indices_to_recv[0]), &(displ[0]),
                   MPITraits<int>::getType(),
                   distributed_data->ccobj_);
    MPI_Allgatherv(&(owned_sizes[0]), no_indices, MPITraits<int>::getType(),
                   &(global_sizes[0]), &(no_indices_to_recv[0]), &(displ[0]),
                   MPITraits<int>::getType(),
                   distributed_data->ccobj_);
    std::vector<int>().swap(owned_global_indices); // free data for reuse.
    // Compute the number of data items to send
    std::vector<int> no_data_send(distributed_data->ccobj_.size());
    for(typename std::vector<int>::iterator begin=no_data_send.begin(),
            i=begin, end=no_data_send.end(); i!=end; ++i)
        *i = std::accumulate(global_sizes.begin()+displ[i-begin],
                            global_sizes.begin()+displ[i-begin+1], std::size_t());
    // free at least some memory that can be reused.
    std::vector<int>().swap(owned_sizes);
    // compute the displacements for receiving with allgatherv
    displ[0]=0;
    std::transform(displ.begin(), displ.end()-1, no_data_send.begin(), displ.begin()+1,
                   std::plus<std::size_t>());
    // Compute the number of data items we will receive
    int no_data_recv = displ[displ.size()-1];//+global_sizes[displ.size()-1];

    // Collect the data to send, gather it
    mover::MoveBuffer<typename DataHandle::DataType> local_data_buffer, global_data_buffer;
    if ( no_data_send[distributed_data->ccobj_.rank()] )
    {
        local_data_buffer.resize(no_data_send[distributed_data->ccobj_.rank()]);
    }
    else
    {
        local_data_buffer.resize(1);
    }
    global_data_buffer.resize(no_data_recv);

    DataGatherer<DataHandle> gatherer(local_data_buffer, data);
    visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gatherer);
    MPI_Allgatherv(&(local_data_buffer.buffer_[0]), no_data_send[distributed_data->ccobj_.rank()],
                   MPITraits<typename DataHandle::DataType>::getType(),
                   &(global_data_buffer.buffer_[0]), &(no_data_send[0]), &(displ[0]),
                   MPITraits<typename DataHandle::DataType>::getType(),
                   distributed_data->ccobj_);
    Entity2IndexDataHandle<DataHandle, codim> edata(*global_data, data);
    int offset=0;
    for(int i=0; i< codim; ++i)
        offset+=global_data->size(i);

    typename std::vector<int>::const_iterator s=global_sizes.begin();
    for(typename std::vector<int>::const_iterator i=global_indices.begin(),
            end=global_indices.end();
        i!=end; ++s, ++i)
    {
        edata.scatter(global_data_buffer, *i-offset, *s);
    }
#endif
}

} // end namespace cpgrid
} // end namespace Dune

#endif

#endif