CpGrid.hpp

Go to the documentation of this file.

//===========================================================================
//
// File: CpGrid.hpp
//
// Created: Fri May 29 20:26:36 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//            B�rd Skaflestad     <bard.skaflestad@sintef.no>
//            Antonella Ritorto   <antonella.ritorto@opm-op.com>
//
// $Date$
//
// $Revision$
//
//===========================================================================
 
/*
  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010, 2014, 2022-2023 Equinor ASA.
  Copyright 2014, 2015 Dr. Blatt - HPC-Simulartion-Software & Services
  Copyright 2015       NTNU
 
  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_CPGRID_HEADER
#define OPM_CPGRID_HEADER
 
// Warning suppression for Dune includes.
#include <opm/grid/utility/platform_dependent/disable_warnings.h>  // Not really needed it seems, but alas.
 
#include <dune/grid/common/grid.hh>
#include <opm/grid/cpgrid/CpGridDataTraits.hpp>
#include <opm/grid/cpgrid/OrientedEntityTable.hpp>
#include <opm/grid/cpgpreprocess/preprocess.h>
#include <opm/grid/utility/platform_dependent/reenable_warnings.h> //  Not really needed it seems, but alas.
#include "common/GridEnums.hpp"   
#include <opm/grid/utility/OpmWellType.hpp>  
 
#include <iostream>
#if ! HAVE_MPI
#include <list>
#endif
 
namespace Opm
{
class EclipseGrid;
class EclipseState;
template<typename Grid, typename GridView> class LookUpData;
template<typename Grid, typename GridView> class LookUpCartesianData;
class NNC;
}
 
namespace Dune
{
    class CpGrid;
 
    namespace cpgrid
    {
    class CpGridData;
    template <int> class Entity;
    template<int,int> class Geometry;
    class HierarchicIterator;
    class IntersectionIterator;
    template<int, PartitionIteratorType> class Iterator;
    class LevelGlobalIdSet;
    class GlobalIdSet;
    class Intersection;
    class IntersectionIterator;
    class IndexSet;
    class IdSet;
    
    }
}
 
void noNNC_check(Dune::CpGrid&,
                 const std::vector<std::array<int,3>>&,
                 const std::vector<std::array<int,3>>&,
                 const std::vector<std::array<int,3>>&,
                 const std::vector<std::string>&);
 
void testCase(const std::string&,
              const Opm::NNC&,
              const std::vector<std::array<int,3>>&,
              const std::vector<std::array<int,3>>&,
              const std::vector<std::array<int,3>>&,
              const std::vector<std::string>&,
              bool);
 
void disjointPatches_check(Dune::CpGrid&,
                           const std::vector<std::array<int,3>>&,
                           const std::vector<std::array<int,3>>&);
 
void lookup_check(const Dune::CpGrid&);
 
void refine_and_check(const Dune::cpgrid::Geometry<3, 3>&,
                      const std::array<int, 3>&,
                      bool);
 
void refinePatch_and_check(Dune::CpGrid&,
                           const std::vector<std::array<int,3>>&,
                           const std::vector<std::array<int,3>>&,
                           const std::vector<std::array<int,3>>&,
                           const std::vector<std::string>&);
 
void refinePatch_and_check(const std::array<int,3>&,
                           const std::array<int,3>&,
                           const std::array<int,3>&);
 
void check_global_refine(const Dune::CpGrid&,
                         const Dune::CpGrid&);
namespace Dune
{
 
    //
    //   CpGridTraits
    //
 
    struct CpGridTraits
    {
        typedef CpGrid Grid;
 
        typedef cpgrid::Intersection LeafIntersection;
        typedef cpgrid::Intersection LevelIntersection;
        typedef cpgrid::IntersectionIterator LeafIntersectionIterator;
        typedef cpgrid::IntersectionIterator LevelIntersectionIterator;
 
        typedef cpgrid::HierarchicIterator HierarchicIterator;
 
        template <int cd>
        struct Codim
        {
             typedef cpgrid::Geometry<3-cd, 3> Geometry;
            //typedef Dune::Geometry<3-cd, 3, CpGrid, cpgrid::Geometry> Geometry;
             typedef cpgrid::Geometry<3-cd, 3> LocalGeometry;
            //typedef Dune::Geometry<3-cd, 3, CpGrid, cpgrid::Geometry> LocalGeometry;
            typedef cpgrid::Entity<cd> Entity;
 
            typedef cpgrid::Iterator<cd, All_Partition> LevelIterator;
 
            typedef cpgrid::Iterator<cd, All_Partition> LeafIterator;
 
            typedef cpgrid::Entity<cd> EntitySeed;
 
            template <PartitionIteratorType pitype>
            struct Partition
            {
                typedef cpgrid::Iterator<cd, pitype> LevelIterator;
                typedef cpgrid::Iterator<cd, pitype> LeafIterator;
            };
        };
 
        template <PartitionIteratorType pitype>
        struct Partition
        {
            typedef Dune::GridView<DefaultLevelGridViewTraits<CpGrid> > LevelGridView;
            typedef Dune::GridView<DefaultLeafGridViewTraits<CpGrid> > LeafGridView;
 
        };
 
        typedef Dune::GridView<DefaultLevelGridViewTraits<CpGrid>> LevelGridView;
        typedef Dune::GridView<DefaultLeafGridViewTraits<CpGrid>> LeafGridView;
 
        typedef cpgrid::IndexSet LevelIndexSet;
        typedef cpgrid::IndexSet LeafIndexSet;
        typedef cpgrid::GlobalIdSet GlobalIdSet;
        typedef GlobalIdSet LocalIdSet;
 
        using Communication = cpgrid::CpGridDataTraits::Communication;
        using CollectiveCommunication = cpgrid::CpGridDataTraits::CollectiveCommunication;
    };
 
    //
    //   CpGridFamily
    //
 
    struct CpGridFamily
    {
        typedef CpGridTraits Traits;
    };
 
    //
    //   CpGrid
    //
 
    class CpGrid
        : public GridDefaultImplementation<3, 3, double, CpGridFamily>
    {
        friend class cpgrid::CpGridData;
        friend class cpgrid::Entity<0>;
        friend class cpgrid::Entity<1>;
        friend class cpgrid::Entity<2>;
        friend class cpgrid::Entity<3>;
        template<typename Grid, typename GridView> friend class Opm::LookUpData;
        template<typename Grid, typename GridView> friend class Opm::LookUpCartesianData;
        template<int dim>
        friend cpgrid::Entity<dim> createEntity(const CpGrid&,int,bool);
        friend void ::noNNC_check(Dune::CpGrid&,
                                  const std::vector<std::array<int,3>>&,
                                  const std::vector<std::array<int,3>>&,
                                  const std::vector<std::array<int,3>>&,
                                  const std::vector<std::string>&);
        friend void ::testCase(const std::string&,
                               const Opm::NNC&,
                               const std::vector<std::array<int,3>>&,
                               const std::vector<std::array<int,3>>&,
                               const std::vector<std::array<int,3>>&,
                               const std::vector<std::string>&,
                               bool);
        friend void ::disjointPatches_check(Dune::CpGrid&,
                                            const std::vector<std::array<int,3>>&,
                                            const std::vector<std::array<int,3>>&);
        friend void ::lookup_check(const Dune::CpGrid&);
        friend
        void ::refine_and_check(const Dune::cpgrid::Geometry<3,3>&,
                                const std::array<int,3>&,
                                bool);
        friend
        void ::refinePatch_and_check(Dune::CpGrid&,
                                     const std::vector<std::array<int,3>>&,
                                     const std::vector<std::array<int,3>>&,
                                     const std::vector<std::array<int,3>>&,
                                     const std::vector<std::string>&);
        friend
        void ::refinePatch_and_check(const std::array<int,3>&,
                                     const std::array<int,3>&,
                                     const std::array<int,3>&);
        friend
        void ::check_global_refine(const Dune::CpGrid&,
                                   const Dune::CpGrid&);
        
    public:
 
        // --- Typedefs ---
 
 
        typedef CpGridFamily GridFamily;
 
 
        // --- Methods ---
 
 
        CpGrid();
 
        CpGrid(MPIHelper::MPICommunicator comm);
 
 
 
        void readSintefLegacyFormat(const std::string& grid_prefix);
 
 
        void writeSintefLegacyFormat(const std::string& grid_prefix) const;
 
 
#if HAVE_ECL_INPUT
        std::vector<std::size_t> processEclipseFormat(const Opm::EclipseGrid* ecl_grid,
                                                      Opm::EclipseState* ecl_state,
                                                      bool periodic_extension, bool turn_normals, bool clip_z,
                                                      bool pinchActive);
 
        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);
 
#endif
 
        void processEclipseFormat(const grdecl& input_data, bool remove_ij_boundary, bool turn_normals = false);
 
 
 
 
        void createCartesian(const std::array<int, 3>& dims,
                             const std::array<double, 3>& cellsize,
                             const std::array<int, 3>& shift = {0,0,0});
 
        const std::array<int, 3>& logicalCartesianSize() const;
 
        const std::vector<int>& globalCell() const;
 
        const std::vector<std::shared_ptr<Dune::cpgrid::CpGridData>>& chooseData() const;
 
        void getIJK(const int c, std::array<int,3>& ijk) const;
 
        bool uniqueBoundaryIds() const;
 
        void setUniqueBoundaryIds(bool uids);
       
 
        // --- Dune interface below ---
 
        // \@{
        std::string name() const;
        
        int maxLevel() 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>::LevelIterator lbegin (int level) const;
        template<int codim, PartitionIteratorType PiType>
        typename Traits::template Codim<codim>::template Partition<PiType>::LevelIterator lend (int level) const;
 
        template<int codim>
        typename Traits::template Codim<codim>::LeafIterator leafbegin() const;
        template<int codim>
        typename Traits::template Codim<codim>::LeafIterator leafend() 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;
 
        int size (int level, int codim) const;
 
        int size (int codim) const;
 
        int size (int level, GeometryType type) const;
 
        int size (GeometryType type) const;
 
        const Traits::GlobalIdSet& globalIdSet() const;
 
        const Traits::LocalIdSet& localIdSet() const;
        
        const Traits::LevelIndexSet& levelIndexSet(int level) const;
 
        const Traits::LeafIndexSet& leafIndexSet() const;
 
        void globalRefine (int);
 
        const std::vector<Dune::GeometryType>& geomTypes(const int) const;
 
        template <int codim>
        cpgrid::Entity<codim> entity(const cpgrid::Entity<codim>& seed) const;
 
        void addLgrUpdateLeafView(const std::array<int,3>& cells_per_dim, const std::array<int,3>& startIJK,
                                  const std::array<int,3>& endIJK,  const std::string& lgr_name);
 
        void addLgrsUpdateLeafView(const std::vector<std::array<int,3>>& cells_per_dim_vec,
                                   const std::vector<std::array<int,3>>& startIJK_vec,
                                   const std::vector<std::array<int,3>>& endIJK_vec,
                                   const std::vector<std::string>& lgr_name_vec);
 
        // @brief TO BE DONE
        const std::map<std::string,int>& getLgrNameToLevel() 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> getEclCentroid(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> getEclCentroid(const cpgrid::Entity<0>& elem) const;
 
        // @brief Return parent (coarse) intersection (face) of a refined face on the leaf grid view, whose neighboring cells
        //        are two: one coarse cell (equivalent to its origin cell from level 0), and one refined cell
        //        from certain LGR
        Dune::cpgrid::Intersection getParentIntersectionFromLgrBoundaryFace(const Dune::cpgrid::Intersection& intersection) const;
 
        /*  No refinement implemented. GridDefaultImplementation's methods will be used.
 
 
        bool mark(int refCount, const typename Traits::template Codim<0>::EntityPointer & e)
        {
            return hostgrid_->mark(refCount, getHostEntity<0>(*e));
        }
 
 
        int getMark(const typename Traits::template Codim<0>::EntityPointer & e) const
        {
            return hostgrid_->getMark(getHostEntity<0>(*e));
        }
 
        bool preAdapt() {
            return hostgrid_->preAdapt();
        }
 
 
        bool adapt()
        {
            return hostgrid_->adapt();
        }
 
        void postAdapt() {
            return hostgrid_->postAdapt();
        }
 
        end of refinement section */
 
        unsigned int overlapSize(int) const;
 
 
        unsigned int ghostSize(int) const;
        
        unsigned int overlapSize(int, int) const;
 
        unsigned int ghostSize(int, int) const;
 
        unsigned int numBoundarySegments() const;
 
        void setZoltanParams(const std::map<std::string,std::string>& params); 
      
 
        // loadbalance is not part of the grid interface therefore we skip it.
 
        bool loadBalance(int overlapLayers=1, bool useZoltan=true)
        {
            using std::get;
            return get<0>(scatterGrid(defaultTransEdgeWgt, false, nullptr, false, nullptr, true, overlapLayers, useZoltan ));
        }
 
        // loadbalance is not part of the grid interface therefore we skip it.
 
        std::pair<bool,std::vector<std::pair<std::string,bool>>>
        loadBalance(const std::vector<cpgrid::OpmWellType> * wells,
                    const double* transmissibilities = nullptr,
                    int overlapLayers=1, bool useZoltan=true)
        {
            return scatterGrid(defaultTransEdgeWgt, false, wells, false, transmissibilities, false, overlapLayers, useZoltan);
        }
 
        // loadbalance is not part of the grid interface therefore we skip it.
 
        std::pair<bool,std::vector<std::pair<std::string,bool>>>
        loadBalance(EdgeWeightMethod method, const std::vector<cpgrid::OpmWellType> * wells,
                    const double* transmissibilities = nullptr, bool ownersFirst=false,
                    bool addCornerCells=false, int overlapLayers=1,
                    bool useZoltan = true)
        {
            return scatterGrid(method, ownersFirst, wells, false, transmissibilities, addCornerCells, overlapLayers, useZoltan);
        }
 
        template<class DataHandle>
        std::pair<bool, std::vector<std::pair<std::string,bool> > >
        loadBalance(DataHandle& data,
                    const std::vector<cpgrid::OpmWellType> * wells,
                    const double* transmissibilities = nullptr,
                    int overlapLayers=1, bool useZoltan = true)
        {
            auto ret = loadBalance(wells, transmissibilities, overlapLayers, useZoltan);
            using std::get;
            if (get<0>(ret))
            {
                scatterData(data);
            }
            return ret;
        }
 
        template<class DataHandle>
        std::pair<bool, std::vector<std::pair<std::string,bool> > >
        loadBalance(DataHandle& data, EdgeWeightMethod method,
                    const std::vector<cpgrid::OpmWellType> * wells,
                    bool serialPartitioning,
                    const double* transmissibilities = nullptr, bool ownersFirst=false,
                    bool addCornerCells=false, int overlapLayers=1, bool useZoltan = true,
                    double zoltanImbalanceTol = 1.1,
                    bool allowDistributedWells = false)
        {
            auto ret = scatterGrid(method, ownersFirst, wells, serialPartitioning, transmissibilities,
                                   addCornerCells, overlapLayers, useZoltan, zoltanImbalanceTol, allowDistributedWells);
            using std::get;
            if (get<0>(ret))
            {
                scatterData(data);
            }
            return ret;
        }
 
        template<class DataHandle>
        std::pair<bool, std::vector<std::pair<std::string,bool> > >
        loadBalance(DataHandle& data, const std::vector<int>& parts,
                    const std::vector<cpgrid::OpmWellType> * wells,
                    bool ownersFirst=false,
                    bool addCornerCells=false, int overlapLayers=1)
        {
            using std::get;
            auto ret = scatterGrid(defaultTransEdgeWgt,  ownersFirst, wells,
                                   /* serialPartitioning = */ false,
                                   /* transmissibilities = */ {},
                                   addCornerCells, overlapLayers, /* useZoltan =*/ false,
                                   /* zoltanImbalanceTol (ignored) = */ 0.0,
                                   /* allowDistributedWells = */ true, parts);
            using std::get;
            if (get<0>(ret))
            {
                scatterData(data);
            }
            return ret;
        }
 
        template<class DataHandle>
        bool loadBalance(DataHandle& data,
                         decltype(data.fixedSize(0,0)) overlapLayers=1, bool useZoltan = true)
        {
            // decltype usage needed to tell the compiler not to use this function if first
            // argument is std::vector but rather loadbalance by parts
            bool ret = loadBalance(overlapLayers, useZoltan);
            if (ret)
            {
                scatterData(data);
            }
            return ret;
        }
 
        bool loadBalance(const std::vector<int>& parts, bool ownersFirst=false,
                         bool addCornerCells=false, int overlapLayers=1)
        {
            using std::get;
            return get<0>(scatterGrid(defaultTransEdgeWgt,  ownersFirst, /* wells = */ {},
                                      /* serialPartitioning = */ false,
                                      /* trabsmissibilities = */ {},
                                      addCornerCells, overlapLayers, /* useZoltan =*/ false,
                                      /* zoltanImbalanceTol (ignored) = */ 0.0,
                                      /* allowDistributedWells = */ true, parts));
        }
 
        template<class DataHandle>
        bool loadBalance(DataHandle& data, const std::vector<int>& parts, bool ownersFirst=false,
                         bool addCornerCells=false, int overlapLayers=1)
        {
            bool ret = loadBalance(parts, ownersFirst, addCornerCells, overlapLayers);
            if (ret)
            {
                scatterData(data);
            }
            return ret;
        }
 
         std::vector<int>
         zoltanPartitionWithoutScatter(const std::vector<cpgrid::OpmWellType>* wells,
                                       const double* transmissibilities,
                                       const int     numParts,
                                       const double  zoltanImbalanceTol) const;
 
        template<class DataHandle>
        void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir, int /*level*/) const
        {
            communicate(data, iftype, dir);
        }
 
        template<class DataHandle>
        void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const;
 
        const typename CpGridTraits::Communication& comm () const;
 
        // ------------ End of Dune interface, start of simplified interface --------------
 
 
        // enum { dimension = 3 }; // already defined
 
        typedef Dune::FieldVector<double, 3> Vector;
 
 
        const std::vector<double>& zcornData() const;
 
 
        // Topology
        int numCells() const;
 
        int numFaces() const;
 
        int numVertices() const;
 
 
        int numCellFaces(int cell) const;
 
        int cellFace(int cell, int local_index) const;
 
        const cpgrid::OrientedEntityTable<0,1>::row_type cellFaceRow(int cell) const;
 
        int faceCell(int face, int local_index) const;
      
        int numCellFaces() const;
 
        int numFaceVertices(int face) const;
 
        int faceVertex(int face, int local_index) const;
 
        double cellCenterDepth(int cell_index) const;
 
 
        const Vector faceCenterEcl(int cell_index, int face, const Dune::cpgrid::Intersection& intersection) const;
 
        const Vector faceAreaNormalEcl(int face) const;
 
 
        // Geometry
        const Vector& vertexPosition(int vertex) const;
 
        double faceArea(int face) const;
 
        const Vector& faceCentroid(int face) const;
 
        const Vector& faceNormal(int face) const;
 
        double cellVolume(int cell) const;
 
        const Vector& cellCentroid(int cell) const;
 
        template<int codim>
        class CentroidIterator
            : public RandomAccessIteratorFacade<CentroidIterator<codim>,
                                                FieldVector<double, 3>,
                                                const FieldVector<double, 3>&, int>
        {
        public:
            typedef typename std::vector<cpgrid::Geometry<3-codim, 3> >::const_iterator
            GeometryIterator;
            CentroidIterator(GeometryIterator iter)
            : iter_(iter)
            {}
 
            const FieldVector<double,3>& dereference() const
            {
                return iter_->center();
            }
            void increment()
            {
                ++iter_;
            }
            const FieldVector<double,3>& elementAt(int n)
            {
                return iter_[n]->center();
            }
            void advance(int n){
                iter_+=n;
            }
            void decrement()
            {
                --iter_;
            }
            int distanceTo(const CentroidIterator& o)
            {
                return o-iter_;
            }
            bool equals(const CentroidIterator& o) const{
                return o==iter_;
            }
        private:
            GeometryIterator iter_;
        };
 
        CentroidIterator<0> beginCellCentroids() const;
 
        CentroidIterator<1> beginFaceCentroids() const;
 
        // Extra
        int boundaryId(int face) const;
 
        template<class Cell2FacesRowIterator>
        int
        faceTag(const Cell2FacesRowIterator& cell_face) const;
 
 
        // ------------ End of simplified interface --------------
 
        //------------- methods not in the DUNE grid interface.
 
 
 
        template<class DataHandle>
        void scatterData(DataHandle& handle) const;
 
        template<class DataHandle>
        void gatherData(DataHandle& handle) const;
 
        using InterfaceMap = cpgrid::CpGridDataTraits::InterfaceMap;
 
        const InterfaceMap& cellScatterGatherInterface() const;
 
        const InterfaceMap& pointScatterGatherInterface() const;
 
        void switchToGlobalView();
 
        void switchToDistributedView();
 
#if HAVE_MPI
        using ParallelIndexSet = cpgrid::CpGridDataTraits::ParallelIndexSet;
        using RemoteIndices = cpgrid::CpGridDataTraits::RemoteIndices;
 
        using CommunicationType = cpgrid::CpGridDataTraits::CommunicationType;
 
        const CommunicationType& cellCommunication() const;
 
        ParallelIndexSet& getCellIndexSet();
 
        RemoteIndices& getCellRemoteIndices();
 
        const ParallelIndexSet& getCellIndexSet() const;
 
        const RemoteIndices& getCellRemoteIndices() const;
#endif
 
        const std::vector<int>& sortedNumAquiferCells() const;
 
    private:
        std::pair<bool, std::vector<std::pair<std::string,bool> > >
        scatterGrid(EdgeWeightMethod method,
                    bool ownersFirst,
                    const std::vector<cpgrid::OpmWellType> * wells,
                    bool serialPartitioning,
                    const double* transmissibilities,
                    bool addCornerCells,
                    int overlapLayers,
                    bool useZoltan = true,
                    double zoltanImbalanceTol = 1.1,
                    bool allowDistributedWells = true,
                    const std::vector<int>& input_cell_part = {});
 
        std::vector<std::shared_ptr<cpgrid::CpGridData>> data_;
        cpgrid::CpGridData* current_view_data_;
        std::vector<std::shared_ptr<cpgrid::CpGridData>> distributed_data_;
        std::map<std::string,int> lgr_names_ = {{"GLOBAL", 0}};
        std::shared_ptr<InterfaceMap> cell_scatter_gather_interfaces_;
        /*
         * @brief Interface for scattering and gathering point data.
         *
         * @warning Will only update owner cells
         */
        std::shared_ptr<InterfaceMap> point_scatter_gather_interfaces_;
        std::shared_ptr<cpgrid::GlobalIdSet> global_id_set_ptr_;
 
 
        std::map<std::string,std::string> zoltanParams;
 
    }; // end Class CpGrid
 
} // end namespace Dune
 
#include <opm/grid/cpgrid/Entity.hpp>
#include <opm/grid/cpgrid/Iterators.hpp>
#include <opm/grid/cpgrid/CpGridData.hpp>
 
 
namespace Dune
{
 
    namespace Capabilities
    {
        template <>
        struct hasEntity<CpGrid, 0>
        {
            static const bool v = true;
        };
 
        template <>
        struct hasEntity<CpGrid, 3>
        {
            static const bool v = true;
        };
 
        template<>
        struct canCommunicate<CpGrid,0>
        {
            static const bool v = true;
        };
 
        template<>
        struct canCommunicate<CpGrid,3>
        {
            static const bool v = true;
        };
 
        template <>
        struct hasBackupRestoreFacilities<CpGrid>
        {
            static const bool v = false;
        };
 
    }
 
    template<class DataHandle>
    void CpGrid::communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const
    {
        current_view_data_->communicate(data, iftype, dir);
    }
 
 
    template<class DataHandle>
    void CpGrid::scatterData(DataHandle& handle) const
    {
#if HAVE_MPI
        if(distributed_data_.empty())
            OPM_THROW(std::runtime_error, "Moving Data only allowed with a load balanced grid!");
        distributed_data_[0]->scatterData(handle, data_[0].get(), distributed_data_[0].get(), cellScatterGatherInterface(),
                                          pointScatterGatherInterface());
#else
        // Suppress warnings for unused argument.
        (void) handle;
#endif
    }
 
    template<class DataHandle>
    void CpGrid::gatherData(DataHandle& handle) const
    {
#if HAVE_MPI
        if(distributed_data_.empty())
            OPM_THROW(std::runtime_error, "Moving Data only allowed with a load balance grid!");
        distributed_data_[0]->gatherData(handle, data_[0].get(), distributed_data_[0].get());
#else
        // Suppress warnings for unused argument.
        (void) handle;
#endif
    }
 
 
    template<class Cell2FacesRowIterator>
    int
    CpGrid::faceTag(const Cell2FacesRowIterator& cell_face) const
    {
        // Note that this relies on the following implementation detail:
        // The grid is always constructed such that the interior faces constructed
        // with orientation set to true are
        // oriented along the positive IJK direction. Oriented means that
        // the first cell attached to face has the lower index.
        // For faces along the boundary (only one cell, always  attached at index 0)
        // the orientation has to be determined by the orientation of the cell.
        // If it is true then in UnstructuredGrid it would be stored at index 0,
        // otherwise at index 1.
        const int cell = cell_face.getCellIndex();
        const int face = *cell_face;
        assert (0 <= cell);  assert (cell < numCells());
        assert (0 <= face);  assert (face < numFaces());
 
        typedef cpgrid::OrientedEntityTable<1,0>::row_type F2C;
 
        const cpgrid::EntityRep<1> f(face, true);
        const F2C&     f2c = current_view_data_->face_to_cell_[f];
        const face_tag tag = current_view_data_->face_tag_[f];
 
        assert ((f2c.size() == 1) || (f2c.size() == 2));
 
        int inside_cell = 0;
 
        if ( f2c.size() == 2 ) // Two cells => interior
        {
            if ( f2c[1].index() == cell )
            {
                inside_cell = 1;
            }
        }
        const bool normal_is_in = ! f2c[inside_cell].orientation();
 
        switch (tag) {
        case I_FACE:
            //                    LEFT : RIGHT
            return normal_is_in ? 0    : 1; // min(I) : max(I)
        case J_FACE:
            //                    BACK : FRONT
            return normal_is_in ? 2    : 3; // min(J) : max(J)
        case K_FACE:
            // Note: TOP at min(K) as 'z' measures *depth*.
            //                    TOP  : BOTTOM
            return normal_is_in ? 4    : 5; // min(K) : max(K)
        case NNC_FACE:
            // For nnc faces we return the otherwise unused value -1.
            return -1;
        default:
            OPM_THROW(std::logic_error, "Unhandled face tag. This should never happen!");
        }
    }
 
    template<int dim>
    cpgrid::Entity<dim> createEntity(const CpGrid&, int, bool);
 
} // namespace Dune
 
#include <opm/grid/cpgrid/PersistentContainer.hpp>
#include <opm/grid/cpgrid/CartesianIndexMapper.hpp>
#include "cpgrid/Intersection.hpp"
#include "cpgrid/Geometry.hpp"
#include "cpgrid/Indexsets.hpp"
 
#endif // OPM_CPGRID_HEADER