CpGrid.hpp

//===========================================================================
//
// 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>

#include <dune/grid/common/grid.hh>

#include <opm/grid/utility/platform_dependent/reenable_warnings.h>

#include <opm/grid/common/GridEnums.hpp>

#include <opm/grid/cpgrid/CpGridDataTraits.hpp>
#include <opm/grid/cpgrid/DefaultGeometryPolicy.hpp>
#include <opm/grid/cpgrid/OrientedEntityTable.hpp>

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

#include <opm/grid/utility/OpmWellType.hpp>

#include <set>

namespace Opm
{
struct NNCdata;
class EclipseGrid;
class EclipseState;
}

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;

    }
}

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<int dim>
        friend cpgrid::Entity<dim> createEntity(const CpGrid&,int,bool);

    public:

        // --- Typedefs ---


        typedef CpGridFamily GridFamily;


        // --- Methods ---


        CpGrid();

        explicit CpGrid(MPIHelper::MPICommunicator comm);

#if HAVE_OPM_COMMON
        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,
                             bool edge_conformal);

        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 edge_conformal = false);
#endif // HAVE_OPM_COMMON

        void processEclipseFormat(const grdecl& input_data,
                                  bool remove_ij_boundary,
                                  bool turn_normals = false,
                                  bool edge_conformal = 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>>& currentData() const;

        std::vector<std::shared_ptr<Dune::cpgrid::CpGridData>>& currentData();

        const Dune::cpgrid::CpGridData& currentLeafData() const;

        Dune::cpgrid::CpGridData& currentLeafData();

        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 refCount, bool throwOnFailure = false);

        const std::vector<Dune::GeometryType>& geomTypes(const int) const;

        template <int codim>
        cpgrid::Entity<codim> entity(const cpgrid::Entity<codim>& seed) const;

        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,
                                   const std::vector<std::string>& lgr_parent_grid_name_vec = std::vector<std::string>{});

        void autoRefine(const std::array<int,3>& nxnynz);

        // @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;

        Dune::cpgrid::Intersection getParentIntersectionFromLgrBoundaryFace(const Dune::cpgrid::Intersection& intersection) 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();

        bool refineAndUpdateGrid(bool throwOnFailure,
                                 const std::vector<std::array<int,3>>& cells_per_dim_vec,
                                 const std::vector<int>& assignRefinedLevel,
                                 const std::vector<std::string>& lgr_name_vec,
                                 const std::vector<std::array<int,3>>& startIJK_vec = std::vector<std::array<int,3>>{},
                                 const std::vector<std::array<int,3>>& endIJK_vec = std::vector<std::array<int,3>>{});

        void postAdapt();

    private:
        void updateCornerHistoryLevels(const std::vector<std::vector<std::array<int,2>>>& cornerInMarkedElemWithEquivRefinedCorner,
                                       const std::map<std::array<int,2>,std::array<int,2>>& elemLgrAndElemLgrCorner_to_refinedLevelAndRefinedCorner,
                                       const std::unordered_map<int,std::array<int,2>>& adaptedCorner_to_elemLgrAndElemLgrCorner,
                                       const int& corner_count,
                                       const std::vector<std::array<int,2>>& preAdaptGrid_corner_history,
                                       const int& preAdaptMaxLevel,
                                       const int& newLevels);

        void globalIdsPartitionTypesLgrAndLeafGrids(const std::vector<int>& assignRefinedLevel,
                                                    const std::vector<std::array<int,3>>& cells_per_dim_vec,
                                                    const std::vector<int>& lgr_with_at_least_one_active_cell);

        void getFirstChildGlobalIds([[maybe_unused]] std::vector<int>& parentToFirstChildGlobalIds);
    public:
        void syncDistributedGlobalCellIds();
    private:

        void computeGlobalCellLgr(const int& level, const std::array<int,3>& startIJK, std::vector<int>& global_cell_lgr);

        void computeGlobalCellLeafGridViewWithLgrs(std::vector<int>& global_cell_leaf);

    private:
        void markElemAssignLevelDetectActiveLgrs(const std::vector<std::array<int,3>>& startIJK_vec,
                                                 const std::vector<std::array<int,3>>& endIJK_vec,
                                                 std::vector<int>& assignRefinedLevel,
                                                 std::vector<int>& lgr_with_at_least_one_active_cell);

        void populateCellIndexSetRefinedGrid(int level);

        void populateCellIndexSetLeafGridView();

        void populateLeafGlobalIdSet();

    public:

        std::vector<std::unordered_map<std::size_t, std::size_t>> mapLocalCartesianIndexSetsToLeafIndexSet() const;

        std::vector<std::array<int,2>> mapLeafIndexSetToLocalCartesianIndexSets() const;

        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 setPartitioningParams(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,
                         int partitionMethod = Dune::PartitionMethod::zoltan,
                         double imbalanceTol = 1.1,
                         int level =-1)
        {
            using std::get;
            return get<0>(scatterGrid(/* edgeWeightMethod = */ defaultTransEdgeWgt,
                                      /* ownersFirst = */ false,
                                      /* wells = */ nullptr,
                                      /* possibleFutureConnections = */ {},
                                      /* serialPartitioning = */ false,
                                      /* transmissibilities = */ nullptr,
                                      /* addCornerCells = */ true,
                                      overlapLayers,
                                      partitionMethod,
                                      imbalanceTol,
                                      /* allowDistributedWells = */ false,
                                      /* input_cell_part = */ {},
                                      level));
        }

        // loadbalance is not part of the grid interface therefore we skip it.

        bool loadBalanceSerial(int overlapLayers=1,
                               int partitionMethod = Dune::PartitionMethod::zoltan,
                               int edgeWeightMethod = Dune::EdgeWeightMethod::defaultTransEdgeWgt,
                               double imbalanceTol = 1.1,
                               int level = -1)
        {
            using std::get;
            return get<0>(scatterGrid(EdgeWeightMethod(edgeWeightMethod),
                                      /* ownersFirst = */ false,
                                      /* wells = */ nullptr,
                                      /* possibleFutureConnections = */ {},
                                      /* serialPartitioning = */ true,
                                      /* transmissibilities = */ nullptr,
                                      /* addCornerCells = */ true,
                                      overlapLayers,
                                      partitionMethod,
                                      imbalanceTol,
                                      /* allowDistributedWells = */ false,
                                      /* input_cell_part = */ {},
                                      level));
        }

        // 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 std::unordered_map<std::string, std::set<int>>& possibleFutureConnections = {},
                    const double* transmissibilities = nullptr,
                    int overlapLayers=1, int partitionMethod=Dune::PartitionMethod::zoltanGoG,
                    int level = -1)
        {
            return scatterGrid(defaultTransEdgeWgt, /* ownersFirst = */ false, wells, possibleFutureConnections,
                               /* serialPartitioning = */ false, transmissibilities, /* addCornerCells = */ false,
                               overlapLayers, partitionMethod, /* imbalanceTol = */ 1.1, /* allowDistributeWells = */ false,
                               /* input_cell_part = */ {}, level);
        }

        // 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 std::unordered_map<std::string, std::set<int>>& possibleFutureConnections = {},
                    const double* transmissibilities = nullptr, bool ownersFirst=false,
                    bool addCornerCells=false, int overlapLayers=1,
                    int partitionMethod = Dune::PartitionMethod::zoltanGoG,
                    double imbalanceTol = 1.1,
                    int level = -1)
        {
            return scatterGrid(method, ownersFirst, wells, possibleFutureConnections,  /* serialPartitioning = */ false,
                               transmissibilities, addCornerCells, overlapLayers, partitionMethod, imbalanceTol,
                               /* allowDistributeWells = */ false, /* input_cell_part = */ {}, level);
        }

        template<class DataHandle>
        std::pair<bool, std::vector<std::pair<std::string,bool> > >
        loadBalance(DataHandle& data,
                    const std::vector<cpgrid::OpmWellType> * wells,
                    const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections = {},
                    const double* transmissibilities = nullptr,
                    int overlapLayers=1, int partitionMethod = 1, int level =-1)
        {
            auto ret = loadBalance(wells, possibleFutureConnections, transmissibilities, overlapLayers, partitionMethod, level);
            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,
                    const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                    bool serialPartitioning,
                    const double* transmissibilities = nullptr, bool ownersFirst=false,
                    bool addCornerCells=false, int overlapLayers=1, int partitionMethod = Dune::PartitionMethod::zoltanGoG,
                    double imbalanceTol = 1.1,
                    bool allowDistributedWells = false,
                    bool useTransToFilterOverlap = true)
        {
            auto ret = scatterGrid(method, ownersFirst, wells, possibleFutureConnections, serialPartitioning, transmissibilities,
                                   addCornerCells, overlapLayers, partitionMethod, imbalanceTol, allowDistributedWells,
                                   /* input_cell_parts = */ std::vector<int>{}, /* level = */ 0,
                                   useTransToFilterOverlap);
            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,
                    const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections = {},
                    bool ownersFirst=false,
                    bool addCornerCells=false, int overlapLayers=1)
        {
            using std::get;
            auto ret = scatterGrid(defaultTransEdgeWgt,  ownersFirst, wells,
                                   possibleFutureConnections,
                                   /* serialPartitioning = */ false,
                                   /* transmissibilities = */ {},
                                   addCornerCells, overlapLayers, /* partitionMethod =*/ Dune::PartitionMethod::simple,
                                   /* imbalanceTol (ignored) = */ 0.0,
                                   /* allowDistributedWells = */ true, parts, /* level = */ 0);
            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, int partitionMethod = Dune::PartitionMethod::zoltan)
        {
            // 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, partitionMethod);
            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, /* partitionMethod =*/ Dune::PartitionMethod::simple,
                                      /* imbalanceTol (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 std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                                       const double* transmissibilities,
                                       const int     numParts,
                                       const double  imbalanceTol) 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(int level = -1) const;

        int numFaces(int level = -1) const;

        int numVertices() const;


        int numCellFaces(int cell, int level = -1) const;

        int cellFace(int cell, int local_index, int level = -1) const;

        const cpgrid::OrientedEntityTable<0,1>::row_type cellFaceRow(int cell) const;

        int faceCell(int face, int local_index, int level = -1) 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;
            explicit 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,
                    const std::unordered_map<std::string, std::set<int>>& possibleFutureConnections,
                    bool serialPartitioning,
                    const double* transmissibilities,
                    bool addCornerCells,
                    int overlapLayers,
                    int partitionMethod = Dune::PartitionMethod::zoltanGoG,
                    double imbalanceTol = 1.1,
                    bool allowDistributedWells = true,
                    const std::vector<int>& input_cell_part = {},
                    int level = -1,
                    bool useTransToFilterOverlap = true);

        std::vector<std::shared_ptr<cpgrid::CpGridData>> data_;
        std::vector<std::shared_ptr<cpgrid::CpGridData>> distributed_data_;
        std::vector<std::shared_ptr<cpgrid::CpGridData>>* current_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> partitioningParams;

    }; // 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_data_->back()->communicate(data, iftype, dir);
    }


    template<class DataHandle>
    void CpGrid::scatterData([[maybe_unused]] DataHandle& handle) const
    {
#if HAVE_MPI
        if (distributed_data_.empty()) {
            OPM_THROW(std::runtime_error, "Moving Data only allowed with a load balanced grid!");
        } else {
            distributed_data_[0]->scatterData(handle, data_[0].get(),
                                              distributed_data_[0].get(),
                                              cellScatterGatherInterface(),
                                              pointScatterGatherInterface());
        }
#endif
    }

    template<class DataHandle>
    void CpGrid::gatherData([[maybe_unused]] DataHandle& handle) const
    {
#if HAVE_MPI
        if (distributed_data_.empty()) {
            OPM_THROW(std::runtime_error, "Moving Data only allowed with a load balance grid!");
        } else {
            distributed_data_[0]->gatherData(handle, data_[0].get(), distributed_data_[0].get());
        }
#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_data_->back()->face_to_cell_[f];
        const face_tag tag = current_data_->back()->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