int main(int argc, char** argv)
try
{
    Opm::parameter::ParameterGroup param(argc, argv);
    Dune::CpGrid grid;
    grid.init(param);
    grid.setUniqueBoundaryIds(true);
    GridInterfaceEuler<Dune::CpGrid> gi(grid);
    typedef FlowBC FBC;
    std::array<FBC, 6> fcond = {{ FBC(FBC::Periodic, 1.0e5),
                                   FBC(FBC::Periodic, -1.0e5),
                                   FBC(FBC::Periodic, 0.0),
                                   FBC(FBC::Periodic, 0.0),
                                   FBC(FBC::Neumann, 0.0),
                                   FBC(FBC::Neumann, 0.0) }};
    typedef SatBC SBC;
    std::array<SBC, 6> scond = {{ SBC(SBC::Periodic, 0.0),
                                   SBC(SBC::Periodic, 0.0),
                                   SBC(SBC::Periodic, 0.0),
                                   SBC(SBC::Periodic, 0.0),
                                   SBC(SBC::Dirichlet, 0.0),
                                   SBC(SBC::Dirichlet, 0.0) }};
    BasicBoundaryConditions<true, true> bcs;
    createPeriodic(bcs, gi, fcond, scond);
    std::cout << bcs;
}
catch (const std::exception &e) {
    std::cerr << "Program threw an exception: " << e.what() << "\n";
    throw;
}

void extractInternalFaces(const Dune::CpGrid& grid,
                          Eigen::Array<int, Eigen::Dynamic, 1>& internal_faces,
                          Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>& nbi)
{
    // Extracts the internal faces of the grid.
    // These are stored in internal_faces.
    int nf=numFaces(grid);
    int num_internal=0;
    for(int f=0; f<nf; ++f)
    {
        if(grid.faceCell(f, 0)<0 || grid.faceCell(f, 1)<0)
            continue;
        ++num_internal;
    }
    // std::cout << num_internal << " internal faces." << std::endl;
    nbi.resize(num_internal, 2);
    internal_faces.resize(num_internal);
    int fi = 0;

    for (int f = 0; f < nf; ++f) {
        if(grid.faceCell(f, 0)>=0 && grid.faceCell(f, 1)>=0) {
            internal_faces[fi] = f;
            nbi(fi,0) = grid.faceCell(f, 0);
            nbi(fi,1) = grid.faceCell(f, 1);
            ++fi;
        }
    }
}

int main(int argc, char** argv)
{
    Opm::parameter::ParameterGroup param(argc, argv);
    Dune::MPIHelper::instance(argc,argv);

    // Make a grid
    Dune::CpGrid grid;

    Opm::EclipseGridParser parser(param.get<std::string>("filename"));
    double z_tol = param.getDefault<double>("z_tolerance", 0.0);
    std::tr1::array<int,3> cartDims;
    build_grid(parser, z_tol, grid, cartDims);

    // Make the grid interface
    Opm::GridInterfaceEuler<Dune::CpGrid> g(grid);

    // Reservoir properties.
    typedef ReservoirPropertyCapillaryAnisotropicRelperm<3> RP;
    RP res_prop;
    res_prop.init(parser, grid.globalCell());

    assign_permeability<3>(res_prop, g.numberOfCells(), 0.1*Opm::unit::darcy);
    test_flowsolver<3>(g, res_prop);

#if 0
    // Make flow equation boundary conditions.
    // Pressure 1.0e5 on the left, 0.0 on the right.
    // Recall that the boundary ids range from 1 to 6 for the cartesian edges,
    // and that boundary id 0 means interiour face/intersection.
    typedef FlowBoundaryCondition BC;
    FlowBoundaryConditions flow_bcond(7);
    flow_bcond[1] = BC(BC::Dirichlet, 1.0e5);
    flow_bcond[2] = BC(BC::Dirichlet, 0.0);

    // Make transport equation boundary conditions.
    // The default one is fine (sat = 1.0 on inflow).
    SaturationBoundaryConditions sat_bcond(7);

    // No injection or production.
    SparseVector<double> injection_rates(g.numberOfCells());

    // Make a solver.
    typedef EulerUpstream<GridInterface,
                          RP,
                          SaturationBoundaryConditions> TransportSolver;
    TransportSolver transport_solver(g, res_prop, sat_bcond, injection_rates);

    // Define a flow field with constant velocity
    Dune::FieldVector<double, 3> vel(0.0);
    vel[0] = 2.0;
    vel[1] = 1.0;
    TestSolution<GridInterface> flow_solution(g, vel);
    // Solve a step.
    double time = 1.0;
    std::vector<double> sat(g.numberOfCells(), 0.0);
    Dune::FieldVector<double, 3> gravity(0.0);
    gravity[2] = -9.81;
    transport_solver.transportSolve(sat, time, gravity, flow_solution);
#endif
}

void Opm::initCPGrid(Dune::CpGrid& grid, const Opm::ParameterGroup& param) {
    std::string fileformat = param.get<std::string>("fileformat");
    if (fileformat == "sintef_legacy") {
        std::string grid_prefix = param.get<std::string>("grid_prefix");
        grid.readSintefLegacyFormat(grid_prefix);
    } else if (fileformat == "eclipse") {
        std::string filename = param.get<std::string>("filename");
        if (param.has("z_tolerance")) {
            std::cerr << "****** Warning: z_tolerance parameter is obsolete, use PINCH in deck input instead\n";
        }
        bool periodic_extension = param.getDefault<bool>("periodic_extension", false);
        bool turn_normals = param.getDefault<bool>("turn_normals", false);
        Opm::Parser parser;
        auto deck = parser.parseFile(filename);
        Opm::EclipseGrid inputGrid(deck);
        grid.processEclipseFormat(inputGrid, periodic_extension , turn_normals );
    } else if (fileformat == "cartesian") {
        std::array<int, 3> dims = {{ param.getDefault<int>("nx", 1),
                                     param.getDefault<int>("ny", 1),
                                     param.getDefault<int>("nz", 1) }};
        std::array<double, 3> cellsz = {{ param.getDefault<double>("dx", 1.0),
                                          param.getDefault<double>("dy", 1.0),
                                          param.getDefault<double>("dz", 1.0) }};
        grid.createCartesian(dims, cellsz);
    } else {
        OPM_THROW(std::runtime_error, "Unknown file format string: " << fileformat);
    }
}

 /// \brief Constructor.
 /// \param sendGrid   The grid that the data is attached to when sending.
 /// \param recvGrid   The grid that the data is attached to when receiving.
 /// \param sendState  The state where we will retieve the values to be sent.
 /// \param recvState The state where we will store the received values.
 BlackoilStateDataHandle(const Dune::CpGrid& sendGrid,
                         const Dune::CpGrid& recvGrid,
                         const BlackoilState& sendState,
                         BlackoilState& recvState)
     : sendGrid_(sendGrid), recvGrid_(recvGrid), sendState_(sendState), recvState_(recvState)
 {
     // construction does not resize surfacevol and hydroCarbonState. Do it manually.
     recvState.surfacevol().resize(recvGrid.numCells()*sendState.numPhases(),
                                   std::numeric_limits<double>::max());
     recvState.hydroCarbonState().resize(recvGrid.numCells());
 }

// Extracts the information about the data decomposition from the grid for dune-istl
void extractParallelGridInformationToISTL(const Dune::CpGrid& grid, boost::any& anyComm)
{
    if(grid.comm().size()>1)
    {
        // this is a parallel run with distributed data.
        Dune::CpGrid& mgrid=const_cast<Dune::CpGrid&>(grid);
        Dune::CpGrid::ParallelIndexSet& idx=mgrid.getCellIndexSet();
        Dune::CpGrid::RemoteIndices& ridx=mgrid.getCellRemoteIndices();
        anyComm=boost::any(Opm::ParallelISTLInformation(Dune::stackobject_to_shared_ptr(idx),
                                                        Dune::stackobject_to_shared_ptr(ridx),
                                                        grid.comm()));
    }
}

void check_cpgrid()
//-----------------------------------------------------------------------------
{
    std::cout << '\n' << "CpGrid<" << refinement << ">\n" << std::endl;

    Dune::CpGrid grid;
    Dune::array<int   , 3> dims;
    std::fill(dims.begin(), dims.end(), 1 << refinement);
    Dune::array<double, 3> cell_sz;
    std::fill(cell_sz.begin(), cell_sz.end(), 1.0 / (1 << refinement));

    grid.createCartesian(dims, cell_sz);

    // Test the interface
    Dune::GridInterfaceEuler<Dune::CpGrid> gie(grid);
    test_evaluator<3>(gie);
#if 0
    test_flowsolver<3>(gie);
#endif
}

 inline void setupGridAndPropsEclipse(Opm::DeckConstPtr deck,
                                      bool periodic_extension,
                                      bool turn_normals,
                                      bool clip_z,
                                      bool unique_bids,
                                      double perm_threshold,
                                      const std::string& rock_list,
                                      bool use_jfunction_scaling,
                                      double sigma,
                                      double theta,
                                      Dune::CpGrid& grid,
                                      ResProp<3>& res_prop)
 {
     grid.processEclipseFormat(deck, periodic_extension, turn_normals, clip_z);
     const std::string* rl_ptr = (rock_list == "no_list") ? 0 : &rock_list;
     res_prop.init(deck, grid.globalCell(), perm_threshold, rl_ptr, use_jfunction_scaling, sigma, theta);
     if (unique_bids) {
         grid.setUniqueBoundaryIds(true);
     }
 }

    void outputStateVtk(const Dune::CpGrid& grid,
                        const Opm::BlackoilState& state,
                        const int step,
                        const std::string& output_dir)
    {
        // Write data in VTK format.
        std::ostringstream vtkfilename;
        std::ostringstream vtkpath;
        vtkpath << output_dir << "/vtk_files";
        vtkpath << "/output-" << std::setw(3) << std::setfill('0') << step;
        boost::filesystem::path fpath(vtkpath.str());
        try {
            create_directories(fpath);
        }
        catch (...) {
            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
        }
        vtkfilename << "output-" << std::setw(3) << std::setfill('0') << step;
#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
        Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafGridView(), Dune::VTK::nonconforming);
#else
        Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafView(), Dune::VTK::nonconforming);
#endif
        writer.addCellData(state.saturation(), "saturation", state.numPhases());
        writer.addCellData(state.pressure(), "pressure", 1);
        
        std::vector<double> cell_velocity;
        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
                                  AutoDiffGrid::numFaces(grid),
                                  AutoDiffGrid::beginFaceCentroids(grid),
                                  AutoDiffGrid::faceCells(grid),
                                  AutoDiffGrid::beginCellCentroids(grid),
                                  AutoDiffGrid::beginCellVolumes(grid),
                                  AutoDiffGrid::dimensions(grid),
                                  state.faceflux(), cell_velocity);
        writer.addCellData(cell_velocity, "velocity", Dune::CpGrid::dimension);
        writer.pwrite(vtkfilename.str(), vtkpath.str(), std::string("."), Dune::VTK::ascii);
    }

void build_grid(const Opm::EclipseGridParser& parser,
                const double z_tol, Dune::CpGrid& grid,
                std::tr1::array<int,3>& cartDims)
{
    Opm::EclipseGridInspector insp(parser);

    grdecl g;
    cartDims[0] = g.dims[0] = insp.gridSize()[0];
    cartDims[1] = g.dims[1] = insp.gridSize()[1];
    cartDims[2] = g.dims[2] = insp.gridSize()[2];

    g.coord = &parser.getFloatingPointValue("COORD")[0];
    g.zcorn = &parser.getFloatingPointValue("ZCORN")[0];

    if (parser.hasField("ACTNUM")) {
        g.actnum = &parser.getIntegerValue("ACTNUM")[0];
        grid.processEclipseFormat(g, z_tol, false, false);
    } else {
        std::vector<int> dflt_actnum(g.dims[0] * g.dims[1] * g.dims[2], 1);

        g.actnum = &dflt_actnum[0];
        grid.processEclipseFormat(g, z_tol, false, false);
    }
}

int numCellFaces(const Dune::CpGrid& grid)
{
    return grid.numCellFaces();    
}

    inline void setupGridAndProps(const Opm::parameter::ParameterGroup& param,
                                  Dune::CpGrid& grid,
                                  ResProp<3>& res_prop)
    {
        // Initialize grid and reservoir properties.
        // Parts copied from Dune::CpGrid::init().
        std::string fileformat = param.getDefault<std::string>("fileformat", "cartesian");
        if (fileformat == "sintef_legacy") {
            std::string grid_prefix = param.get<std::string>("grid_prefix");
            grid.readSintefLegacyFormat(grid_prefix);
            OPM_MESSAGE("Warning: We do not yet read legacy reservoir properties. Using defaults.");
            res_prop.init(grid.size(0));
        } else if (fileformat == "eclipse") {
            std::string ecl_file = param.get<std::string>("filename");

            Opm::ParseContext parseContext;
            Opm::ParserPtr parser(new Opm::Parser());
            Opm::DeckConstPtr deck(parser->parseFile(ecl_file , parseContext));
            if (param.has("z_tolerance")) {
                std::cerr << "****** Warning: z_tolerance parameter is obsolete, use PINCH in deck input instead\n";
            }
            bool periodic_extension = param.getDefault<bool>("periodic_extension", false);
            bool turn_normals = param.getDefault<bool>("turn_normals", false);
            grid.processEclipseFormat(deck, periodic_extension, turn_normals);
            // Save EGRID file in case we are writing ECL output.
            if (param.getDefault("output_ecl", false)) {
                OPM_THROW(std::runtime_error, "Saving to EGRID files is not yet implemented");
                /*
                boost::filesystem::path ecl_path(ecl_file);
                const std::vector<int>& globalCell = grid.globalCell();
                ecl_path.replace_extension(".EGRID");
                parser.saveEGRID(ecl_path.string() , (int) globalCell.size() , &globalCell[0]);
                */
            }
            double perm_threshold_md = param.getDefault("perm_threshold_md", 0.0);
            double perm_threshold = Opm::unit::convert::from(perm_threshold_md, Opm::prefix::milli*Opm::unit::darcy);
            std::string rock_list = param.getDefault<std::string>("rock_list", "no_list");
            std::string* rl_ptr = (rock_list == "no_list") ? 0 : &rock_list;
            bool use_j = param.getDefault("use_jfunction_scaling", useJ<ResProp<3> >());
            double sigma = 1.0;
            double theta = 0.0;
            if (use_j) {
                sigma = param.getDefault("sigma", sigma);
                theta = param.getDefault("theta", theta);
            }
            if (param.has("viscosity1") || param.has("viscosity2")) {
                double v1 = param.getDefault("viscosity1", 0.001);
                double v2 = param.getDefault("viscosity2", 0.003);
                res_prop.setViscosities(v1, v2);
            }
            res_prop.init(deck, grid.globalCell(), perm_threshold, rl_ptr,
                          use_j, sigma, theta);
        } else if (fileformat == "cartesian") {
            std::array<int, 3> dims = {{ param.getDefault<int>("nx", 1),
                                    param.getDefault<int>("ny", 1),
                                    param.getDefault<int>("nz", 1) }};
            std::array<double, 3> cellsz = {{ param.getDefault<double>("dx", 1.0),
                                         param.getDefault<double>("dy", 1.0),
                                         param.getDefault<double>("dz", 1.0) }};
            grid.createCartesian(dims, cellsz);
            double default_poro = param.getDefault("default_poro", 0.2);
            double default_perm_md = param.getDefault("default_perm_md", 100.0);
            double default_perm = Opm::unit::convert::from(default_perm_md, Opm::prefix::milli*Opm::unit::darcy);
            OPM_MESSAGE("Warning: For generated cartesian grids, we use uniform reservoir properties.");
            res_prop.init(grid.size(0), default_poro, default_perm);
        } else {
            OPM_THROW(std::runtime_error, "Unknown file format string: " << fileformat);
        }
        if (param.getDefault("use_unique_boundary_ids", false)) {
            grid.setUniqueBoundaryIds(true);
        }
    }

double cellCentroidCoordinate(const Dune::CpGrid& grid, int cell_index,
                              int coordinate)
{
    return grid.cellCentroid(cell_index)[coordinate];
}

inline
void distributeGridAndData( Dune::CpGrid& grid,
                            Opm::DeckConstPtr deck,
                            EclipseStateConstPtr eclipseState,
                            BlackoilState& state,
                            BlackoilPropsAdFromDeck& properties,
                            DerivedGeology& geology,
                            std::shared_ptr<BlackoilPropsAdFromDeck::MaterialLawManager>& material_law_manager,
                            std::vector<double>& threshold_pressures,
                            boost::any& parallelInformation,
                            const bool useLocalPerm)
{
    Dune::CpGrid global_grid ( grid );
    global_grid.switchToGlobalView();

    // distribute the grid and switch to the distributed view
    grid.loadBalance(eclipseState, geology.transmissibility().data());
    grid.switchToDistributedView();
    std::vector<int> compressedToCartesianIdx;
    Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
    typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
    auto distributed_material_law_manager = std::make_shared<MaterialLawManager>();
    distributed_material_law_manager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
    // copy the values from the global to the local MaterialLawManager
    // We should actually communicate these to be future proof. But that is
    // really, really cumbersome for the underlying vector<shared_ptr>
    // where the classes pointed to even have more shared_ptr stored in them.
    typedef Dune::CpGrid::ParallelIndexSet IndexSet;
    const IndexSet& local_indices  = grid.getCellIndexSet();
    for ( auto index : local_indices )
    {
        distributed_material_law_manager->materialLawParamsPointerReferenceHack(index.local()) =
            material_law_manager->materialLawParamsPointerReferenceHack(index.global());

        distributed_material_law_manager->oilWaterScaledEpsInfoDrainagePointerReferenceHack(index.local()) =
            material_law_manager->oilWaterScaledEpsInfoDrainagePointerReferenceHack(index.global());
    }
    BlackoilPropsAdFromDeck distributed_props(properties,
                                              distributed_material_law_manager,
                                              grid.numCells());
    BlackoilState distributed_state(grid.numCells(), grid.numFaces(), state.numPhases());
    BlackoilStateDataHandle state_handle(global_grid, grid,
                                         state, distributed_state);
    BlackoilPropsDataHandle props_handle(properties,
                                         distributed_props);
    grid.scatterData(state_handle);
    grid.scatterData(props_handle);
    // Create a distributed Geology. Some values will be updated using communication
    // below
    DerivedGeology distributed_geology(grid,
                                       distributed_props, eclipseState,
                                       useLocalPerm, geology.gravity());
    GeologyDataHandle geo_handle(global_grid, grid,
                                 geology, distributed_geology);
    grid.scatterData(geo_handle);

    std::vector<double> distributed_pressures;

    if( !threshold_pressures.empty() ) // Might be empty if not specified
    {
        if( threshold_pressures.size() !=
            static_cast<std::size_t>(UgGridHelpers::numFaces(global_grid)) )
        {
            OPM_THROW(std::runtime_error, "NNCs not yet supported for parallel runs. "
                      << UgGridHelpers::numFaces(grid) << " faces but " <<
                      threshold_pressures.size()<<" threshold pressure values");
        }
        distributed_pressures.resize(UgGridHelpers::numFaces(grid));
        ThresholdPressureDataHandle press_handle(global_grid, grid,
                                                 threshold_pressures,
                                                 distributed_pressures);
        grid.scatterData(press_handle);
    }

    // copy states
    properties           = distributed_props;
    geology              = distributed_geology;
    state                = distributed_state;
    material_law_manager = distributed_material_law_manager;
    threshold_pressures   = distributed_pressures;
    extractParallelGridInformationToISTL(grid, parallelInformation);
}

FaceCentroidTraits<Dune::CpGrid>::ValueType
faceCentroid(const Dune::CpGrid& grid, int face_index)
{
    return grid.faceCentroid(face_index);
}

double cellVolume(const  Dune::CpGrid& grid, int cell_index)
{
    return grid.cellVolume(cell_index);
}

double faceArea(const Dune::CpGrid& grid, int face_index)
{
    return grid.faceArea(face_index);
}