/*
  Copyright 2013 Statoil ASA.

  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/>.
 */

#include <cstdio>
#include <ctime>
#include <iostream>
#include <math.h>
#include <memory>
#include <optional>
#include <numeric>
#include <stdexcept>
#include <unistd.h>

#include <opm/common/utility/FileSystem.hpp>
#include <opm/common/utility/OpmInputError.hpp>

#define BOOST_TEST_MODULE EclipseGridTests
#include <boost/test/unit_test.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>

#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
#include <opm/parser/eclipse/Deck/DeckSection.hpp>

#include <opm/parser/eclipse/Units/UnitSystem.hpp>

#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/GridDims.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/PinchMode.hpp>

#include <opm/parser/eclipse/Parser/Parser.hpp>

#include <opm/io/eclipse/EclFile.hpp>

BOOST_AUTO_TEST_CASE(CreateMissingDIMENS_throws) {
    Opm::Deck deck;
    Opm::Parser parser;
    deck.addKeyword( Opm::DeckKeyword( parser.getKeyword("RUNSPEC" )));
    deck.addKeyword( Opm::DeckKeyword( parser.getKeyword("GRID" )));
    deck.addKeyword( Opm::DeckKeyword( parser.getKeyword("EDIT" )));

    BOOST_CHECK_THROW(Opm::EclipseGrid{ deck } , std::invalid_argument);
}

static Opm::Deck createDeckHeaders() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

static Opm::Deck createDeckDIMENS() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 13 17 19/\n"
        "GRID\n"
        "EDIT\n"
        "\n";
    Opm::Parser parser;
    return parser.parseString( deckData);
}

static Opm::Deck createDeckSPECGRID() {
    const char* deckData =
        "GRID\n"
        "SPECGRID \n"
        "  13 17 19 / \n"
        "COORD\n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "EDIT\n"
        "\n";
    Opm::Parser parser;
    return parser.parseString( deckData);
}

static Opm::Deck createDeckMissingDIMS() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "GRID\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

BOOST_AUTO_TEST_CASE(MissingDimsThrows) {
    Opm::Deck deck = createDeckMissingDIMS();
    BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument);
}

BOOST_AUTO_TEST_CASE(HasGridKeywords) {
    Opm::Deck deck = createDeckHeaders();
    BOOST_CHECK( !Opm::EclipseGrid::hasCornerPointKeywords( deck ));
    BOOST_CHECK( !Opm::EclipseGrid::hasCartesianKeywords( deck ));
}

BOOST_AUTO_TEST_CASE(CreateGridNoCells) {
    Opm::Deck deck = createDeckHeaders();
    BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument);

    const Opm::GridDims grid( deck);
    BOOST_CHECK_EQUAL( 10U , grid.getNX());
    BOOST_CHECK_EQUAL( 10U , grid.getNY());
    BOOST_CHECK_EQUAL( 10U , grid.getNZ());

    BOOST_CHECK_EQUAL(10U, grid[0]);
    BOOST_CHECK_EQUAL(10U, grid[2]);
    BOOST_CHECK_THROW( grid[10], std::invalid_argument);

    BOOST_CHECK_EQUAL( 1000U , grid.getCartesianSize());
}

BOOST_AUTO_TEST_CASE(CheckGridIndex) {
    Opm::EclipseGrid grid(17, 19, 41); // prime time

    auto v_start = grid.getIJK(0);
    BOOST_CHECK_EQUAL(v_start[0], 0);
    BOOST_CHECK_EQUAL(v_start[1], 0);
    BOOST_CHECK_EQUAL(v_start[2], 0);

    auto v_end = grid.getIJK(17*19*41 - 1);
    BOOST_CHECK_EQUAL(v_end[0], 16);
    BOOST_CHECK_EQUAL(v_end[1], 18);
    BOOST_CHECK_EQUAL(v_end[2], 40);

    auto v167 = grid.getIJK(167);
    BOOST_CHECK_EQUAL(v167[0], 14);
    BOOST_CHECK_EQUAL(v167[1], 9);
    BOOST_CHECK_EQUAL(v167[2], 0);
    BOOST_CHECK_EQUAL(grid.getGlobalIndex(14, 9, 0), 167U);

    auto v5723 = grid.getIJK(5723);
    BOOST_CHECK_EQUAL(v5723[0], 11);
    BOOST_CHECK_EQUAL(v5723[1], 13);
    BOOST_CHECK_EQUAL(v5723[2], 17);
    BOOST_CHECK_EQUAL(grid.getGlobalIndex(11, 13, 17), 5723U);

    BOOST_CHECK_EQUAL(std::size_t(17 * 19 * 41), grid.getCartesianSize());
}

static Opm::Deck createCPDeck() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "COORD\n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

static Opm::Deck createPinchedCPDeck() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "COORD\n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "PINCH \n"
        "  0.2 / \n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

static Opm::Deck createPinchedNOGAPCPDeck() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "COORD\n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "PINCH \n"
        "  0.2 NOGAP / \n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

static Opm::Deck createMinpvDefaultCPDeck() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "COORD\n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "MINPV \n"
        "  / \n"
        "MINPVFIL \n"
        "  / \n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}


static Opm::Deck createMinpvCPDeck() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "COORD\n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "MINPV \n"
        "  10 / \n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}



static Opm::Deck createCARTDeck() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "DX\n"
        "1000*0.25 /\n"
        "DYV\n"
        "10*0.25 /\n"
        "DZ\n"
        "1000*0.25 /\n"
        "TOPS\n"
        "100*0.25 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}


static Opm::Deck createCARTDeckDEPTHZ() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "DXV\n"
        "10*0.25 /\n"
        "DYV\n"
        "10*0.25 /\n"
        "DZV\n"
        "10*0.25 /\n"
        "DEPTHZ\n"
        "121*0.25 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}


static Opm::Deck createCARTInvalidDeck() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "DX\n"
        "1000*0.25 /\n"
        "DYV\n"
        "1000*0.25 /\n"
        "DZ\n"
        "1000*0.25 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

BOOST_AUTO_TEST_CASE(CREATE_SIMPLE) {
    Opm::EclipseGrid grid(10,20,30);

    BOOST_CHECK_EQUAL( grid.getNX() , 10U );
    BOOST_CHECK_EQUAL( grid.getNY() , 20U );
    BOOST_CHECK_EQUAL( grid.getNZ() , 30U );
    BOOST_CHECK_EQUAL( grid.getCartesianSize() , 6000U );
}

BOOST_AUTO_TEST_CASE(DEPTHZ_EQUAL_TOPS) {
    Opm::Deck deck1 = createCARTDeck();
    Opm::Deck deck2 = createCARTDeckDEPTHZ();

    Opm::EclipseGrid grid1( deck1 );
    Opm::EclipseGrid grid2( deck2 );

    BOOST_CHECK( grid1.equal( grid2 ) );

    {
        BOOST_CHECK_THROW( grid1.getCellVolume(1000) , std::invalid_argument);
        BOOST_CHECK_THROW( grid1.getCellVolume(10,0,0) , std::invalid_argument);
        BOOST_CHECK_THROW( grid1.getCellVolume(0,10,0) , std::invalid_argument);
        BOOST_CHECK_THROW( grid1.getCellVolume(0,0,10) , std::invalid_argument);

        for (size_t g=0; g < 1000; g++)
            BOOST_CHECK_CLOSE( grid1.getCellVolume(g) , 0.25*0.25*0.25 , 0.001);


        for (size_t k= 0; k < 10; k++)
            for (size_t j= 0; j < 10; j++)
                for (size_t i= 0; i < 10; i++)
                    BOOST_CHECK_CLOSE( grid1.getCellVolume(i, j, k) , 0.25*0.25*0.25 , 0.001 );

    }

    {
        BOOST_CHECK_THROW( grid1.getCellCenter(1000) , std::invalid_argument);
        BOOST_CHECK_THROW( grid1.getCellCenter(10,0,0) , std::invalid_argument);
        BOOST_CHECK_THROW( grid1.getCellCenter(0,10,0) , std::invalid_argument);
        BOOST_CHECK_THROW( grid1.getCellCenter(0,0,10) , std::invalid_argument);

        for (size_t k= 0; k < 10; k++) {
            for (size_t j= 0; j < 10; j++) {
                for (size_t i= 0; i < 10; i++) {
                    auto pos = grid1.getCellCenter(i, j, k);

                    BOOST_CHECK_CLOSE( std::get<0>(pos) , i*0.25 + 0.125, 0.001);
                    BOOST_CHECK_CLOSE( std::get<1>(pos) , j*0.25 + 0.125, 0.001);
                    BOOST_CHECK_CLOSE( std::get<2>(pos) , k*0.25 + 0.125 + 0.25, 0.001);
                }
            }
        }
    }
}

BOOST_AUTO_TEST_CASE(HasCPKeywords) {
    Opm::Deck deck = createCPDeck();
    BOOST_CHECK(  Opm::EclipseGrid::hasCornerPointKeywords( deck ));
    BOOST_CHECK( !Opm::EclipseGrid::hasCartesianKeywords( deck ));
}

BOOST_AUTO_TEST_CASE(HasCartKeywords) {
    Opm::Deck deck = createCARTDeck();
    BOOST_CHECK( !Opm::EclipseGrid::hasCornerPointKeywords( deck ));
    BOOST_CHECK(  Opm::EclipseGrid::hasCartesianKeywords( deck ));
}

BOOST_AUTO_TEST_CASE(HasCartKeywordsDEPTHZ) {
    Opm::Deck deck = createCARTDeckDEPTHZ();
    BOOST_CHECK( !Opm::EclipseGrid::hasCornerPointKeywords( deck ));
    BOOST_CHECK(  Opm::EclipseGrid::hasCartesianKeywords( deck ));
}

BOOST_AUTO_TEST_CASE(HasINVALIDCartKeywords) {
    Opm::Deck deck = createCARTInvalidDeck();
    BOOST_CHECK( !Opm::EclipseGrid::hasCornerPointKeywords( deck ));
    BOOST_CHECK( !Opm::EclipseGrid::hasCartesianKeywords( deck ));
}

BOOST_AUTO_TEST_CASE(CreateMissingGRID_throws) {
    auto deck= createDeckHeaders();
    BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument);
}

static Opm::Deck createInvalidDXYZCARTDeck() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "DX\n"
        "99*0.25 /\n"
        "DY\n"
        "1000*0.25 /\n"
        "DZ\n"
        "1000*0.25 /\n"
        "TOPS\n"
        "1000*0.25 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

BOOST_AUTO_TEST_CASE(CreateCartesianGRID) {
    auto deck = createInvalidDXYZCARTDeck();
    BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument);
}

static Opm::Deck createInvalidDXYZCARTDeckDEPTHZ() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "DX\n"
        "100*0.25 /\n"
        "DY\n"
        "1000*0.25 /\n"
        "DZ\n"
        "1000*0.25 /\n"
        "DEPTHZ\n"
        "101*0.25 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

BOOST_AUTO_TEST_CASE(CreateCartesianGRIDDEPTHZ) {
    auto deck = createInvalidDXYZCARTDeckDEPTHZ();
    BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument);
}

static Opm::Deck createOnlyTopDZCartGrid() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 5 20 /\n"
        "GRID\n"
        "DX\n"
        "1000*0.25 /\n"
        "DY\n"
        "1000*0.25 /\n"
        "DZ\n"
        "101*0.25 /\n"
        "TOPS\n"
        "110*0.25 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}


static Opm::Deck createInvalidDEPTHZDeck1 () {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 5 20 /\n"
        "GRID\n"
        "DXV\n"
        "1000*0.25 /\n"
        "DYV\n"
        "5*0.25 /\n"
        "DZV\n"
        "20*0.25 /\n"
        "DEPTHZ\n"
        "66*0.25 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

BOOST_AUTO_TEST_CASE(CreateCartesianGRIDInvalidDEPTHZ1) {
    auto deck = createInvalidDEPTHZDeck1();
    BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument);
}

static Opm::Deck createInvalidDEPTHZDeck2 () {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 5 20 /\n"
        "GRID\n"
        "DXV\n"
        "10*0.25 /\n"
        "DYV\n"
        "5*0.25 /\n"
        "DZV\n"
        "20*0.25 /\n"
        "DEPTHZ\n"
        "67*0.25 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

BOOST_AUTO_TEST_CASE(CreateCartesianGRIDInvalidDEPTHZ2) {
    auto deck = createInvalidDEPTHZDeck2();
    BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument);
}

BOOST_AUTO_TEST_CASE(CreateCartesianGRIDOnlyTopLayerDZ) {
    Opm::Deck deck = createOnlyTopDZCartGrid();
    Opm::EclipseGrid grid( deck );
    BOOST_CHECK_EQUAL( 10U , grid.getNX( ));
    BOOST_CHECK_EQUAL(  5U , grid.getNY( ));
    BOOST_CHECK_EQUAL( 20U , grid.getNZ( ));
    BOOST_CHECK_EQUAL( 1000U , grid.getNumActive());
}

BOOST_AUTO_TEST_CASE(AllActiveExportActnum) {
    Opm::Deck deck = createOnlyTopDZCartGrid();
    Opm::EclipseGrid grid( deck );

    std::vector<int> actnum = grid.getACTNUM();

    BOOST_CHECK_EQUAL( 1000U , actnum.size());
}

BOOST_AUTO_TEST_CASE(CornerPointSizeMismatchCOORD) {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "COORD\n"
        "  725*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    auto deck = parser.parseString( deckData) ;
    const auto& zcorn = deck.getKeyword("ZCORN");
    BOOST_CHECK_EQUAL( 8000U , zcorn.getDataSize( ));

    BOOST_CHECK_THROW(Opm::EclipseGrid{ deck }, std::invalid_argument);
}


BOOST_AUTO_TEST_CASE(CornerPointSizeMismatchZCORN) {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "COORD\n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8001*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    auto deck = parser.parseString( deckData) ;
    BOOST_CHECK_THROW(Opm::EclipseGrid{ deck }, std::invalid_argument);
}

BOOST_AUTO_TEST_CASE(ResetACTNUM) {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "GRID\n"
        "COORD\n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    auto deck = parser.parseString( deckData) ;

    Opm::EclipseGrid grid( deck);
    BOOST_CHECK_EQUAL( 1000U , grid.getNumActive());
    std::vector<int> actnum(1000);
    actnum[0] = 1;
    actnum[2] = 1;
    actnum[4] = 1;
    actnum[6] = 1;

    grid.resetACTNUM( actnum );

    BOOST_CHECK_EQUAL( 4U , grid.getNumActive() );
    {
        std::vector<int> full(grid.getCartesianSize());
        std::iota(full.begin(), full.end(), 0);

        auto compressed = grid.compressedVector( full );

        BOOST_CHECK_EQUAL( compressed.size() , 4U );
        BOOST_CHECK_EQUAL( compressed[0] , 0 );
        BOOST_CHECK_EQUAL( compressed[1] , 2 );
        BOOST_CHECK_EQUAL( compressed[2] , 4 );
        BOOST_CHECK_EQUAL( compressed[3] , 6 );
    }

    {
        const auto& activeMap = grid.getActiveMap( );
        BOOST_CHECK_EQUAL( 4U , activeMap.size() );
        BOOST_CHECK_EQUAL( 0 , activeMap[0] );
        BOOST_CHECK_EQUAL( 2 , activeMap[1] );
        BOOST_CHECK_EQUAL( 4 , activeMap[2] );
        BOOST_CHECK_EQUAL( 6 , activeMap[3] );
    }

    grid.resetACTNUM();

    BOOST_CHECK_EQUAL( 1000U , grid.getNumActive() );
    {
        const auto&  activeMap = grid.getActiveMap( );
        BOOST_CHECK_EQUAL( 1000U , activeMap.size() );
        BOOST_CHECK_EQUAL( 0 , activeMap[0] );
        BOOST_CHECK_EQUAL( 1 , activeMap[1] );
        BOOST_CHECK_EQUAL( 2 , activeMap[2] );
        BOOST_CHECK_EQUAL( 999 , activeMap[999] );
    }

    actnum.assign(1000, 1);

    actnum[0] = 0;
    actnum[1] = 0;
    actnum[2] = 0;
    actnum[11] = 0;
    actnum[21] = 0;
    actnum[430] = 0;
    actnum[431] = 0;

    grid.resetACTNUM( actnum );

    std::vector<int> actMap = grid.getActiveMap();

    BOOST_CHECK_EQUAL(actMap.size(), 993U);
    BOOST_CHECK_THROW(grid.getGlobalIndex(993), std::out_of_range);
    BOOST_CHECK_EQUAL(grid.getGlobalIndex(0), 3U);
    BOOST_CHECK_EQUAL(grid.getGlobalIndex(33), 38U);
    BOOST_CHECK_EQUAL(grid.getGlobalIndex(450), 457U);
    BOOST_CHECK_EQUAL(grid.getGlobalIndex(1,2,3), 321U);
}

BOOST_AUTO_TEST_CASE(TestCP_example) {
    const char* deckData =

    "RUNSPEC\n"
    "\n"
    "DIMENS\n"
    " 3 2 1 /\n"
    "GRID\n"
    "COORD\n"
    " 2000.0000  2000.0000  2000.0000   1999.9476  2000.0000  2002.9995\n"
    " 2049.9924  2000.0000  2000.8726   2049.9400  2000.0000  2003.8722 \n"
    " 2099.9848  2000.0000  2001.7452   2099.9324  2000.0000  2004.7448 \n"
    " 2149.9772  2000.0000  2002.6179   2149.9248  2000.0000  2005.6174 \n"
    " 2000.0000  2050.0000  2000.0000   1999.9476  2050.0000  2002.9995 \n"
    " 2049.9924  2050.0000  2000.8726   2049.9400  2050.0000  2003.8722 \n"
    " 2099.9848  2050.0000  2001.7452   2099.9324  2050.0000  2004.7448 \n"
    " 2149.9772  2050.0000  2002.6179   2149.9248  2050.0000  2005.6174 \n"
    " 2000.0000  2100.0000  2000.0000   1999.9476  2100.0000  2002.9995 \n"
    " 2049.9924  2100.0000  2000.8726   2049.9400  2100.0000  2003.8722 \n"
    " 2099.9848  2100.0000  2001.7452   2099.9324  2100.0000  2004.7448 \n"
    " 2149.9772  2100.0000  2002.6179   2149.9248  2100.0000  2005.6174 / \n"
    "ZCORN\n"
    " 2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179 \n"
    " 2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179 \n"
    " 2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179 \n"
    " 2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179 \n"
    " 2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174 \n"
    " 2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174 \n"
    " 2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174 \n"
    " 2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174 / \n"
    "EDIT\n"
    "\n";

    Opm::Parser parser;
    auto deck = parser.parseString( deckData) ;

    Opm::EclipseGrid grid( deck);
    BOOST_CHECK_EQUAL( 6U , grid.getNumActive());

    std::vector<int> actnum(6, 0);
    actnum[0] = 1;
    actnum[2] = 1;
    actnum[4] = 1;

    grid.resetACTNUM( actnum );

    BOOST_CHECK_EQUAL( 3U , grid.getNumActive() );
}



BOOST_AUTO_TEST_CASE(ConstructorNORUNSPEC) {
    const char* deckData =
        "GRID\n"
        "SPECGRID \n"
        "  10 10 10 / \n"
        "COORD\n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "PORO\n"
        "  1000*0.15 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    auto deck1 = parser.parseString( deckData) ;
    auto deck2 = createCPDeck();

    Opm::EclipseGrid grid1(deck1);
    Opm::EclipseGrid grid2(deck2);

    BOOST_CHECK(grid1.equal( grid2 ));
}

BOOST_AUTO_TEST_CASE(ConstructorNoSections) {
    const char* deckData =
        "DIMENS \n"
        "  10 10 10 / \n"
        "COORD \n"
        "  726*1 / \n"
        "ZCORN \n"
        "  8000*1 / \n"
        "ACTNUM \n"
        "  1000*1 / \n"
        "PORO\n"
        "  1000*0.15 /\n"
        "\n";

    Opm::Parser parser;
    auto deck1 = parser.parseString( deckData) ;
    auto deck2 = createCPDeck();

    Opm::EclipseGrid grid1(deck1);
    Opm::EclipseGrid grid2(deck2);

    BOOST_CHECK(grid1.equal( grid2 ));
}

BOOST_AUTO_TEST_CASE(ConstructorNORUNSPEC_PINCH) {
    auto deck1 = createCPDeck();
    auto deck2 = createPinchedCPDeck();
    auto deck3 = createPinchedNOGAPCPDeck();

    Opm::EclipseGrid grid1(deck1);
    Opm::EclipseGrid grid2(deck2);
    Opm::EclipseGrid grid3(deck3);

    BOOST_CHECK(!grid1.equal( grid2 ));

    BOOST_CHECK(!grid1.isPinchActive());
    BOOST_CHECK_THROW(grid1.getPinchThresholdThickness(), std::bad_optional_access);
    BOOST_CHECK(grid2.isPinchActive());
    BOOST_CHECK_EQUAL(grid2.getPinchThresholdThickness(), 0.2);
    BOOST_CHECK_EQUAL(grid2.getPinchGapMode(), Opm::PinchMode::ModeEnum::GAP);
    BOOST_CHECK_EQUAL(grid3.getPinchGapMode(), Opm::PinchMode::ModeEnum::NOGAP);
}

BOOST_AUTO_TEST_CASE(ConstructorMINPV) {
    auto deck1 = createCPDeck();
    auto deck2 = createMinpvDefaultCPDeck();
    auto deck3 = createMinpvCPDeck();

    Opm::EclipseGrid grid1(deck1);
    BOOST_CHECK_THROW(Opm::EclipseGrid grid2(deck2), std::invalid_argument);
    Opm::EclipseGrid grid3(deck3);

    BOOST_CHECK(!grid1.equal( grid3 ));
    BOOST_CHECK_EQUAL(grid1.getMinpvMode(), Opm::MinpvMode::ModeEnum::Inactive);
    BOOST_CHECK_EQUAL(grid3.getMinpvMode(), Opm::MinpvMode::ModeEnum::EclSTD);
    BOOST_CHECK_EQUAL(grid3.getMinpvVector()[0], 10.0);
}

static Opm::Deck createActnumDeck() {
    const char* deckData = "RUNSPEC\n"
            "\n"
            "DIMENS \n"
            "  2 2 2 / \n"
            "GRID\n"
            "DXV\n"
            "  2*0.25 /\n"
            "DYV\n"
            "  2*0.25 /\n"
            "DZV\n"
            "  2*0.25 /\n"
            "DEPTHZ\n"
            "  9*0.25 /\n"
            "EQUALS\n"
            " ACTNUM 0 1 1 1 1 1 1 /\n"
            "/ \n"
            "PORO\n"
            "  8*0.15 /\n"
            "FLUXNUM\n"
            "8*0 /\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}


/// creates a deck where the top-layer has ACTNUM = 0 and two partially
/// overlapping 2*2*2 boxes in the center, one [5,7]^3 and one [6,8]^3
/// have ACTNUM = 0
static Opm::Deck createActnumBoxDeck() {
    const char* deckData = "RUNSPEC\n"
            "\n"
            "DIMENS \n"
            "  10 10 10 / \n"
            "GRID\n"
            "DXV\n"
            "  10*0.25 /\n"
            "DYV\n"
            "  10*0.25 /\n"
            "DZV\n"
            "  10*0.25 /\n"
            "DEPTHZ\n"
            "  121*0.25 /\n"
            "PORO \n"
            "  1000*0.15 /\n"
            "EQUALS\n"
            " ACTNUM 0 1 10 1 10 1 1 /\n" // disable top layer
            "/ \n"
            // start box
            "BOX\n"
            "  5 7 5 7 5 7 /\n"
            "ACTNUM \n"
            "    0 0 0 0 0 0 0 0 0\n"
            "    0 0 0 0 0 0 0 0 0\n"
            "    0 0 0 0 0 0 0 0 0\n"
            "/\n"
            "BOX\n" // don't need ENDBOX
            "  6 8 6 8 6 8 /\n"
            "ACTNUM \n"
            "    27*0\n"
            "/\n"
            "ENDBOX\n"
            // end   box
            "FLUXNUM\n"
            "1000*0 /\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

BOOST_AUTO_TEST_CASE(GridBoxActnum) {
    auto deck = createActnumBoxDeck();
    Opm::EclipseState es( deck);
    const auto& fp = es.fieldProps();
    const auto& grid = es.getInputGrid();

    BOOST_CHECK_NO_THROW(fp.get_int("ACTNUM"));

    size_t active = 10 * 10 * 10     // 1000
                    - (10 * 10 * 1)  // - top layer
                    - ( 3 *  3 * 3)  // - [5,7]^3 box
                    - ( 3 *  3 * 3)  // - [6,8]^3 box
                    + ( 2 *  2 * 2); // + inclusion/exclusion

    BOOST_CHECK_NO_THROW(grid.getNumActive());
    BOOST_CHECK_EQUAL(grid.getNumActive(), active);

    BOOST_CHECK_EQUAL(es.getInputGrid().getNumActive(), active);

    {
        size_t active_index = 0;
        // NB: The implementation of this test actually assumes that
        //     the loops are running with z as the outer and x as the
        //     inner direction.
        for (size_t z = 0; z < grid.getNZ(); z++) {
            for (size_t y = 0; y < grid.getNY(); y++) {
                for (size_t x = 0; x < grid.getNX(); x++) {
                    if (z == 0)
                        BOOST_CHECK(!grid.cellActive(x, y, z));
                    else if (x >= 4 && x <= 6 && y >= 4 && y <= 6 && z >= 4 && z <= 6)
                        BOOST_CHECK(!grid.cellActive(x, y, z));
                    else if (x >= 5 && x <= 7 && y >= 5 && y <= 7 && z >= 5 && z <= 7)
                        BOOST_CHECK(!grid.cellActive(x, y, z));
                    else {
                        size_t g = grid.getGlobalIndex( x,y,z );

                        BOOST_CHECK(grid.cellActive(x, y, z));
                        BOOST_CHECK_EQUAL( grid.activeIndex(x,y,z) , active_index );
                        BOOST_CHECK_EQUAL( grid.activeIndex(g) , active_index );

                        active_index++;
                    }
                }
            }
        }

        BOOST_CHECK_THROW( grid.activeIndex(0,0,0) , std::invalid_argument );
    }
}

BOOST_AUTO_TEST_CASE(GridActnumVia3D) {
    auto deck = createActnumDeck();

    Opm::EclipseState es( deck);
    const auto& fp = es.fieldProps();
    const auto& grid = es.getInputGrid();
    Opm::EclipseGrid grid2( grid );

    std::vector<int> actnum = {1, 1, 0, 1, 1, 0, 1, 1};
    Opm::EclipseGrid grid3( grid , actnum);

    BOOST_CHECK_NO_THROW(fp.get_int("ACTNUM"));
    BOOST_CHECK_NO_THROW(grid.getNumActive());
    BOOST_CHECK_EQUAL(grid.getNumActive(), std::size_t(2 * 2 * 2 - 1));

    BOOST_CHECK_NO_THROW(grid2.getNumActive());
    BOOST_CHECK_EQUAL(grid2.getNumActive(), std::size_t(2 * 2 * 2 - 1));

    BOOST_CHECK_EQUAL(grid3.getNumActive(), 6U);
}


BOOST_AUTO_TEST_CASE(GridActnumViaState) {
    auto deck = createActnumDeck();

    BOOST_CHECK_NO_THROW( std::unique_ptr<Opm::EclipseState>(new Opm::EclipseState( deck)));
    Opm::EclipseState es( deck);
    BOOST_CHECK_EQUAL(es.getInputGrid().getNumActive(), std::size_t(2 * 2 * 2 - 1));
}


BOOST_AUTO_TEST_CASE(GridDimsSPECGRID) {
    auto deck =  createDeckSPECGRID();
    auto gd = Opm::GridDims( deck );
    BOOST_CHECK_EQUAL(gd.getNX(), 13U);
    BOOST_CHECK_EQUAL(gd.getNY(), 17U);
    BOOST_CHECK_EQUAL(gd.getNZ(), 19U);
}


BOOST_AUTO_TEST_CASE(GridDimsDIMENS) {
    auto deck =  createDeckDIMENS();
    auto gd = Opm::GridDims( deck );
    BOOST_CHECK_EQUAL(gd.getNX(), 13U);
    BOOST_CHECK_EQUAL(gd.getNY(), 17U);
    BOOST_CHECK_EQUAL(gd.getNZ(), 19U);
}


BOOST_AUTO_TEST_CASE(ProcessedCopy) {
    Opm::EclipseGrid gd(10,10,10);
    std::vector<double> zcorn;
    std::vector<int> actnum;

    zcorn = gd.getZCORN();
    actnum = gd.getACTNUM();

    Opm::EclipseGrid gd1(gd , actnum );
    BOOST_CHECK( gd.equal( gd1 ));
    {
        Opm::EclipseGrid gd2(gd , zcorn.data() , actnum );
        BOOST_CHECK( gd.equal( gd2 ));
    }

    zcorn[0] -= 1.0;
    {
        Opm::EclipseGrid gd2(gd , zcorn.data() , actnum );
        BOOST_CHECK( !gd.equal( gd2 ));
    }

    {
        Opm::EclipseGrid gd2(gd , actnum );
        BOOST_CHECK( gd.equal( gd2 ));
    }

    actnum.assign( gd.getCartesianSize() , 1U);
    actnum[0] = 0;
    {
        Opm::EclipseGrid gd2(gd , actnum );
        BOOST_CHECK( !gd.equal( gd2 ));
        BOOST_CHECK( !gd2.cellActive( 0 ));
    }
}

BOOST_AUTO_TEST_CASE(regularCartGrid) {

    int nx = 3;
    int ny = 4;
    int nz = 5;

    double dx = 25;
    double dy = 35;
    double dz = 2;

    double ref_volume = dx* dy* dz;

    Opm::EclipseGrid grid(nx, ny, nz, dx, dy, dz);

    std::array<int, 3> dims = grid.getNXYZ();

    int nCells = dims[0]*dims[1]*dims[2];

    for (int n=0; n<nCells; n++){
        BOOST_CHECK_CLOSE(grid.getCellVolume(n), ref_volume, 1e-12);
        BOOST_CHECK_CLOSE(grid.getCellThickness(n), 2.0, 1e-12);
    }

    for (int k=0; k< dims[2]; k++){
        double ref_depth = k*dz + dz/2.0;
        for (int j=0; j< dims[1]; j++){
            for (int i=0; i< dims[0]; i++){
                BOOST_CHECK_CLOSE(grid.getCellDepth(i,j,k), ref_depth, 1e-12);
            }
        }
    }

    for (int k=0; k< dims[2]; k++){
        double ref_z = k*dz + dz/2.0;
        for (int j=0; j< dims[1]; j++){
            double ref_y = j*dy + dy/2.0;
            for (int i=0; i< dims[0]; i++){
                double ref_x = i*dx + dx/2.0;
                std::array<double, 3> cc = grid.getCellCenter(i, j, k);
                BOOST_CHECK_CLOSE(cc[0], ref_x, 1e-12);
                BOOST_CHECK_CLOSE(cc[1], ref_y, 1e-12);
                BOOST_CHECK_CLOSE(cc[2], ref_z, 1e-12);
            }
        }
    }
}

BOOST_AUTO_TEST_CASE(ZcornMapper) {

    int nx = 3;
    int ny = 4;
    int nz = 5;

    Opm::EclipseGrid grid(nx, ny, nz);
    Opm::ZcornMapper zmp = grid.zcornMapper( );

    BOOST_CHECK_THROW(zmp.index(nx,1,1,0) , std::invalid_argument);
    BOOST_CHECK_THROW(zmp.index(0,ny,1,0) , std::invalid_argument);
    BOOST_CHECK_THROW(zmp.index(0,1,nz,0) , std::invalid_argument);
    BOOST_CHECK_THROW(zmp.index(0,1,2,8) , std::invalid_argument);

    auto points_adjusted = grid.fixupZCORN();

    std::vector<int> actnum = grid.getACTNUM();
    std::vector<double> zcorn = grid.getZCORN();

    zcorn[42] = zcorn[42] + 2.0;
    zcorn[96] = zcorn[96] + 2.0;

    Opm::EclipseGrid grid2(grid , zcorn.data() , actnum );
    points_adjusted = grid2.getZcornFixed();
    BOOST_CHECK_EQUAL( points_adjusted , 4U );

    points_adjusted = grid2.fixupZCORN();
    BOOST_CHECK_EQUAL( points_adjusted , 0U );

    zcorn = grid.getZCORN();

    BOOST_CHECK( zmp.validZCORN( zcorn ));

    // Manually destroy it - cell internal
    zcorn[ zmp.index(0,0,0,4) ] = zcorn[ zmp.index(0,0,0,0) ] - 0.1;
    BOOST_CHECK( !zmp.validZCORN( zcorn ));
    points_adjusted = zmp.fixupZCORN( zcorn );
    BOOST_CHECK_EQUAL( points_adjusted , 1U );
    BOOST_CHECK( zmp.validZCORN( zcorn ));

    // Manually destroy it - cell 2 cell
    zcorn[ zmp.index(0,0,0,4) ] = zcorn[ zmp.index(0,0,1,0) ] + 0.1;
    BOOST_CHECK( !zmp.validZCORN( zcorn ));
    points_adjusted = zmp.fixupZCORN( zcorn );
    BOOST_CHECK_EQUAL( points_adjusted , 1U );
    BOOST_CHECK( zmp.validZCORN( zcorn ));

    // Manually destroy it - cell 2 cell and cell internal
    zcorn[ zmp.index(0,0,0,4) ] = zcorn[ zmp.index(0,0,1,0) ] + 0.1;
    zcorn[ zmp.index(0,0,0,0) ] = zcorn[ zmp.index(0,0,0,4) ] + 0.1;
    BOOST_CHECK( !zmp.validZCORN( zcorn ));
    points_adjusted = zmp.fixupZCORN( zcorn );
    BOOST_CHECK_EQUAL( points_adjusted , 2U );
    BOOST_CHECK( zmp.validZCORN( zcorn ));
}

BOOST_AUTO_TEST_CASE(MoveTest) {
    int nx = 3;
    int ny = 4;
    int nz = 5;
    Opm::EclipseGrid grid1(nx,ny,nz);
    Opm::EclipseGrid grid2( std::move( grid1 )); // grid2 should be move constructed from grid1

    BOOST_CHECK( !grid1.circle( ));
}

static Opm::Deck radial_missing_INRAD() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 10 /\n"
        "RADIAL\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}




BOOST_AUTO_TEST_CASE(RadialTest) {
    Opm::Deck deck = radial_missing_INRAD();
    BOOST_CHECK_THROW( Opm::EclipseGrid{ deck }, std::invalid_argument);
}

static Opm::Deck radial_keywords_DRV_size_mismatch() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "10 6 12 /\n"
        "RADIAL\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "9*1 /\n"
        "DTHETAV\n"
        "6*60 /\n"
        "DZV\n"
        "12*0.25 /\n"
        "TOPS\n"
        "60*0.0 /\n"
        "PORO \n"
        "  720*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}


static Opm::Deck radial_keywords_DZV_size_mismatch() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "10 6 12 /\n"
        "RADIAL\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "10*1 /\n"
        "DTHETAV\n"
        "6*60 /\n"
        "DZV\n"
        "11*0.25 /\n"
        "TOPS\n"
        "60*0.0 /\n"
        "PORO \n"
        "  720*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

static Opm::Deck radial_keywords_DZ_size_mismatch() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "10 6 12 /\n"
        "RADIAL\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "10*1 /\n"
        "DTHETAV\n"
        "6*60 /\n"
        "DZ\n"
        "660*0.25 /\n"
        "TOPS\n"
        "60*0.0 /\n"
        "PORO \n"
        "  720*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

static Opm::Deck radial_keywords_DTHETAV_size_mismatch() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "10 6 12 /\n"
        "RADIAL\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "10*1 /\n"
        "DTHETAV\n"
        "5*60 /\n"
        "DZV\n"
        "12*0.25 /\n"
        "TOPS\n"
        "60*0.0 /\n"
        "PORO \n"
        "  720*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

//  This is stricter than the ECLIPSE implementation; we assume that
//  *only* the top layer is explicitly given.

static Opm::Deck radial_keywords_TOPS_size_mismatch() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "10 6 12 /\n"
        "RADIAL\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "10*1 /\n"
        "DTHETAV\n"
        "6*60 /\n"
        "DZV\n"
        "12*0.25 /\n"
        "TOPS\n"
        "65*0.0 /\n"
        "PORO \n"
        "  720*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

static Opm::Deck radial_keywords_ANGLE_OVERFLOW() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "10 6 12 /\n"
        "RADIAL\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "10*1 /\n"
        "DTHETAV\n"
        "6*70 /\n"
        "DZV\n"
        "12*0.25 /\n"
        "TOPS\n"
        "60*0.0 /\n"
        "PORO \n"
        "  720*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}


BOOST_AUTO_TEST_CASE(RadialKeywords_SIZE_ERROR) {
    BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_DRV_size_mismatch() } , std::invalid_argument);
    BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_DZV_size_mismatch() } , std::invalid_argument);
    BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_DZ_size_mismatch() } , std::invalid_argument);
    BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_TOPS_size_mismatch() } , std::invalid_argument);
    BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_DTHETAV_size_mismatch() } , std::invalid_argument);
    BOOST_CHECK_THROW( Opm::EclipseGrid{ radial_keywords_ANGLE_OVERFLOW() } , std::invalid_argument);
}


static Opm::Deck spider_details() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "1 5 2 /\n"
        "SPIDER\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "1 /\n"
        "DTHETAV\n"
        "3*90 60 30/\n"
        "DZV\n"
        "2*1 /\n"
        "TOPS\n"
        "5*1.0 /\n"
        "PORO \n"
        "  10*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

BOOST_AUTO_TEST_CASE(SpiderDetails) {
    Opm::Deck deck = spider_details();
    Opm::EclipseGrid grid( deck );

    BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 0 , 0 ) , 0.5*(2*2 - 1)*1, 0.0001);
    BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 3 , 0 ) , sqrt(3.0)*0.25*( 4 - 1 ) , 0.0001);
    auto pos0 = grid.getCellCenter(0,0,0);
    auto pos2 = grid.getCellCenter(0,2,0);

    BOOST_CHECK_CLOSE( std::get<0>(pos0) , 0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<1>(pos0) , 0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<2>(pos0) , 1.50 , 0.0001);

    BOOST_CHECK_CLOSE( std::get<0>(pos2) , -0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<1>(pos2) , -0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<2>(pos2) , 1.50 , 0.0001);

    {
        const auto& p0 = grid.getCornerPos( 0,0,0 , 0 );
        const auto& p6 = grid.getCornerPos( 0,0,0 , 6 );
        BOOST_CHECK_CLOSE( p0[0]*p0[0] + p0[1]*p0[1] , 1.0, 0.0001);
        BOOST_CHECK_CLOSE( p6[0]*p6[0] + p6[1]*p6[1] , 1.0, 0.0001);

        BOOST_CHECK_THROW( grid.getCornerPos( 0,0,0 , 8 ) , std::invalid_argument);
    }
}

static Opm::Deck spider_details_dz() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "1 5 2 /\n"
        "SPIDER\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "1 /\n"
        "DTHETAV\n"
        "3*90 60 30/\n"
        "DZ\n"
        "10*1 /\n"
        "TOPS\n"
        "5*1.0 /\n"
        "PORO \n"
        "  10*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

BOOST_AUTO_TEST_CASE(SpiderDetailsDZ) {
    Opm::Deck deck = spider_details_dz();
    Opm::EclipseGrid grid( deck );

    BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 0 , 0 ) , 0.5*(2*2 - 1)*1, 0.0001);
    BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 3 , 0 ) , sqrt(3.0)*0.25*( 4 - 1 ) , 0.0001);
    auto pos0 = grid.getCellCenter(0,0,0);
    auto pos2 = grid.getCellCenter(0,2,0);

    BOOST_CHECK_CLOSE( std::get<0>(pos0) , 0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<1>(pos0) , 0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<2>(pos0) , 1.50 , 0.0001);

    BOOST_CHECK_CLOSE( std::get<0>(pos2) , -0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<1>(pos2) , -0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<2>(pos2) , 1.50 , 0.0001);

    {
        const auto& p0 = grid.getCornerPos( 0,0,0 , 0 );
        const auto& p6 = grid.getCornerPos( 0,0,0 , 6 );
        BOOST_CHECK_CLOSE( p0[0]*p0[0] + p0[1]*p0[1] , 1.0, 0.0001);
        BOOST_CHECK_CLOSE( p6[0]*p6[0] + p6[1]*p6[1] , 1.0, 0.0001);

        BOOST_CHECK_THROW( grid.getCornerPos( 0,0,0 , 8 ) , std::invalid_argument);
    }
}

static Opm::Deck radial_details() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "1 5 2 /\n"
        "RADIAL\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "1 /\n"
        "DTHETAV\n"
        "3*90 60 30/\n"
        "DZV\n"
        "2*1 /\n"
        "TOPS\n"
        "5*1.0 /\n"
        "PORO \n"
        "  10*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

BOOST_AUTO_TEST_CASE(RadialDetails) {
    Opm::Deck deck = radial_details();
    Opm::EclipseGrid grid( deck );

    BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 0 , 0 ) , 0.75 * M_PI, 0.0001);
    BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 3 , 0 ) , 0.5 * M_PI , 0.0001);
    auto pos0 = grid.getCellCenter(0,0,0);
    auto pos2 = grid.getCellCenter(0,2,0);

    BOOST_CHECK_CLOSE( std::get<0>(pos0) , 0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<1>(pos0) , 0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<2>(pos0) , 1.50 , 0.0001);

    BOOST_CHECK_CLOSE( std::get<0>(pos2) , -0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<1>(pos2) , -0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<2>(pos2) , 1.50 , 0.0001);

    {
        const auto& p0 = grid.getCornerPos( 0,0,0 , 0 );
        const auto& p6 = grid.getCornerPos( 0,0,0 , 6 );
        BOOST_CHECK_CLOSE( p0[0]*p0[0] + p0[1]*p0[1] , 1.0, 0.0001);
        BOOST_CHECK_CLOSE( p6[0]*p6[0] + p6[1]*p6[1] , 1.0, 0.0001);

        BOOST_CHECK_THROW( grid.getCornerPos( 0,0,0 , 8 ) , std::invalid_argument);
    }
}

static Opm::Deck radial_details_dz() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "1 5 2 /\n"
        "RADIAL\n"
        "GRID\n"
        "INRAD\n"
        "1 /\n"
        "DRV\n"
        "1 /\n"
        "DTHETAV\n"
        "3*90 60 30/\n"
        "DZ\n"
        "10*1 /\n"
        "TOPS\n"
        "5*1.0 /\n"
        "PORO \n"
        "  10*0.15 /"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

BOOST_AUTO_TEST_CASE(RadialDetailsDZ) {
    Opm::Deck deck = radial_details_dz();
    Opm::EclipseGrid grid( deck );

    BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 0 , 0 ) , 0.75 * M_PI, 0.0001);
    BOOST_CHECK_CLOSE( grid.getCellVolume( 0 , 3 , 0 ) , 0.5 * M_PI , 0.0001);
    auto pos0 = grid.getCellCenter(0,0,0);
    auto pos2 = grid.getCellCenter(0,2,0);

    BOOST_CHECK_CLOSE( std::get<0>(pos0) , 0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<1>(pos0) , 0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<2>(pos0) , 1.50 , 0.0001);

    BOOST_CHECK_CLOSE( std::get<0>(pos2) , -0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<1>(pos2) , -0.75 , 0.0001);
    BOOST_CHECK_CLOSE( std::get<2>(pos2) , 1.50 , 0.0001);

    {
        const auto& p0 = grid.getCornerPos( 0,0,0 , 0 );
        const auto& p6 = grid.getCornerPos( 0,0,0 , 6 );
        BOOST_CHECK_CLOSE( p0[0]*p0[0] + p0[1]*p0[1] , 1.0, 0.0001);
        BOOST_CHECK_CLOSE( p6[0]*p6[0] + p6[1]*p6[1] , 1.0, 0.0001);

        BOOST_CHECK_THROW( grid.getCornerPos( 0,0,0 , 8 ) , std::invalid_argument);
    }
}

BOOST_AUTO_TEST_CASE(CoordMapper) {
    size_t nx = 10;
    size_t ny = 7;
    Opm::CoordMapper cmp = Opm::CoordMapper( nx , ny );
    BOOST_CHECK_THROW( cmp.index(12,6,0,0), std::invalid_argument );
    BOOST_CHECK_THROW( cmp.index(10,8,0,0), std::invalid_argument );
    BOOST_CHECK_THROW( cmp.index(10,7,5,0), std::invalid_argument );
    BOOST_CHECK_THROW( cmp.index(10,5,1,2), std::invalid_argument );

    BOOST_CHECK_EQUAL( cmp.index(10,7,2,1) + 1 , cmp.size( ));
}

static Opm::Deck createCARTDeckTest3x4x2() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 3 4 2 /\n"
        "GRID\n"
        "DX\n"
        "100 120 110 100 120 110 100 120 110 100 120 110 /\n"
        "DY\n"
        "70  80  85  80  70  80  85  80  70  80  85  80 /\n"
        "DZ\n"
        "12*25 12*35 /\n"
        "TOPS\n"
        "2500 2510 2520  2520 2530 2540  2540 2550 2560  2560 2570 2580 /\n"
        "PORO \n"
        "  24*0.15 /"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

BOOST_AUTO_TEST_CASE(CART_Deck_3x4x2) {

    Opm::Deck deck1 = createCARTDeckTest3x4x2();
    Opm::EclipseGrid grid1( deck1 );

    const std::vector<double> t_coord = grid1.getCOORD();
    const std::vector<double> t_zcorn = grid1.getZCORN();

    std::vector<float> ref_coord = {0, 0, 2500, 0, 0, 2560, 100, 0, 2510, 100, 0, 2570, 220, 0, 2520, 220, 0,
        2580, 330, 0, 2520, 330, 0, 2580, 0, 70, 2520, 0, 70, 2580, 100, 80, 2530, 100, 80, 2590, 220, 85,
        2540, 220, 85, 2600, 330, 85, 2540, 330, 85, 2600, 0, 150, 2540, 0, 150, 2600, 100, 150, 2550, 100,
        150, 2610, 220, 165, 2560, 220, 165, 2620, 330, 165, 2560, 330, 165, 2620, 0, 235, 2560, 0, 235,
        2620, 100, 230, 2570, 100, 230, 2630, 220, 235, 2580, 220, 235, 2640, 330, 235, 2580, 330, 235,
        2640, 0, 315, 2560, 0, 315, 2620, 100, 315, 2570, 100, 315, 2630, 220, 315, 2580, 220, 315, 2640,
        330, 315, 2580, 330, 315, 2640};

    std::vector<float> ref_zcorn = {2500, 2500, 2510, 2510, 2520, 2520, 2500, 2500, 2510, 2510, 2520, 2520,
        2520, 2520, 2530, 2530, 2540, 2540, 2520, 2520, 2530, 2530, 2540, 2540, 2540, 2540, 2550, 2550, 2560,
        2560, 2540, 2540, 2550, 2550, 2560, 2560, 2560, 2560, 2570, 2570, 2580, 2580, 2560, 2560, 2570, 2570,
        2580, 2580, 2525, 2525, 2535, 2535, 2545, 2545, 2525, 2525, 2535, 2535, 2545, 2545, 2545, 2545, 2555,
        2555, 2565, 2565, 2545, 2545, 2555, 2555, 2565, 2565, 2565, 2565, 2575, 2575, 2585, 2585, 2565, 2565,
        2575, 2575, 2585, 2585, 2585, 2585, 2595, 2595, 2605, 2605, 2585, 2585, 2595, 2595, 2605, 2605, 2525,
        2525, 2535, 2535, 2545, 2545, 2525, 2525, 2535, 2535, 2545, 2545, 2545, 2545, 2555, 2555, 2565, 2565,
        2545, 2545, 2555, 2555, 2565, 2565, 2565, 2565, 2575, 2575, 2585, 2585, 2565, 2565, 2575, 2575, 2585,
        2585, 2585, 2585, 2595, 2595, 2605, 2605, 2585, 2585, 2595, 2595, 2605, 2605, 2560, 2560, 2570, 2570,
        2580, 2580, 2560, 2560, 2570, 2570, 2580, 2580, 2580, 2580, 2590, 2590, 2600, 2600, 2580, 2580, 2590,
        2590, 2600, 2600, 2600, 2600, 2610, 2610, 2620, 2620, 2600, 2600, 2610, 2610, 2620, 2620, 2620, 2620,
        2630, 2630, 2640, 2640, 2620, 2620, 2630, 2630, 2640, 2640};

    BOOST_CHECK_EQUAL( t_coord.size() , ref_coord.size());

    for (size_t i=0; i< t_coord.size(); i++) {
        BOOST_CHECK_CLOSE( t_coord[i] , ref_coord[i], 1.0e-5);
    }

    BOOST_CHECK_EQUAL( t_zcorn.size() , ref_zcorn.size());

    for (size_t i=0; i< t_zcorn.size(); i++) {
        BOOST_CHECK_CLOSE( t_zcorn[i] , ref_zcorn[i], 1.0e-5);
    }
}

static Opm::Deck createCARTDeckDEPTHZ_2x3x2() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 2 3 2 /\n"
        "GRID\n"
        "DXV\n"
        "100 120 /\n"
        "DYV\n"
        "70  80  85 /\n"
        "DZV\n"
        "25 35 /\n"
        "DEPTHZ\n"
        "2500 2510 2520  2502 2512 2522  2504 2514 2524  2505 2515 2525 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    return parser.parseString( deckData) ;
}

BOOST_AUTO_TEST_CASE(CART_Deck_DEPTHZ_2x3x2) {

    Opm::Deck deck1 = createCARTDeckDEPTHZ_2x3x2();
    Opm::EclipseGrid grid1( deck1 );

    std::vector<double> ref_coord = { 0, 0, 2500, 0, 0, 2560, 100, 0, 2510, 100, 0, 2570, 220, 0, 2520, 220, 0,
        2580, 0, 70, 2502, 0, 70, 2562, 100, 70, 2512, 100, 70, 2572, 220, 70, 2522, 220, 70, 2582, 0, 150,
        2504, 0, 150, 2564, 100, 150, 2514, 100, 150, 2574, 220, 150, 2524, 220, 150, 2584, 0, 235, 2505, 0,
        235, 2565, 100, 235, 2515, 100, 235, 2575, 220, 235, 2525, 220, 235, 2585 };

    std::vector<double> ref_zcorn = { 2500, 2510, 2510, 2520, 2502, 2512, 2512, 2522, 2502, 2512, 2512, 2522,
        2504, 2514, 2514, 2524, 2504, 2514, 2514, 2524, 2505, 2515, 2515, 2525, 2525, 2535, 2535, 2545, 2527,
        2537, 2537, 2547, 2527, 2537, 2537, 2547, 2529, 2539, 2539, 2549, 2529, 2539, 2539, 2549, 2530, 2540,
        2540, 2550, 2525, 2535, 2535, 2545, 2527, 2537, 2537, 2547, 2527, 2537, 2537, 2547, 2529, 2539, 2539,
        2549, 2529, 2539, 2539, 2549, 2530, 2540, 2540, 2550, 2560, 2570, 2570, 2580, 2562, 2572, 2572, 2582,
        2562, 2572, 2572, 2582, 2564, 2574, 2574, 2584, 2564, 2574, 2574, 2584, 2565, 2575, 2575, 2585 };

    const std::vector<double> t_coord = grid1.getCOORD();
    const std::vector<double> t_zcorn = grid1.getZCORN();

    BOOST_CHECK_EQUAL( t_coord.size() , ref_coord.size());

    for (size_t i=0; i< t_coord.size(); i++) {
        BOOST_CHECK_CLOSE( t_coord[i] , ref_coord[i], 1.0e-5);
    }

    BOOST_CHECK_EQUAL( t_zcorn.size() , ref_zcorn.size());

    for (size_t i=0; i< t_zcorn.size(); i++) {
        BOOST_CHECK_CLOSE( t_zcorn[i] , ref_zcorn[i], 1.0e-5);
    }
}

static Opm::Deck BAD_CP_GRID() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "GRID\n"
        "SPECGRID\n"
        " 2 2 2 1 F /\n"
        "COORD\n"
        "  2002.0000  2002.0000   100.0000   1999.8255  1999.9127   108.4935\n"
        "  2011.9939  2000.0000   100.3490   2009.8194  1999.9127   108.8425\n"
        "  2015.9878  2000.0000   100.6980   2019.8133  1999.9127   109.1915\n"
        "  2000.0000  2009.9985   100.1745   1999.8255  2009.9112   108.6681 \n"
        "  2010.9939  2011.9985   100.5235   2009.8194  2009.9112   109.0170\n"
        "  2019.9878  2009.9985   100.8725   2019.8133  2009.9112   109.3660\n"
        "  2005.0000  2019.9970   100.3490   1999.8255  2019.9097   108.8426\n"
        "  2009.9939  2019.9970   100.6980   2009.8194  2019.9097   109.1916\n"
        "  2016.9878  2019.9970   101.0470   2019.8133  2019.9097   109.5406 /\n"
        "ZCORN\n"
        "    98.0000   100.3490    97.3490   100.6980   100.1745   100.5235\n"
        "   100.5235   100.8725   100.1745   100.5235   100.5235   100.8725\n"
        "   100.3490   101.6980   101.6980   102.5470   102.4973   102.1463\n"
        "   103.2463   104.1953   103.6719   104.0209   104.0209   104.3698\n"
        "   103.6719   104.0209   104.0209   104.3698   103.8464   104.1954\n"
        "   104.1954   104.5444   103.4973   103.8463   103.8463   104.1953\n"
        "   103.6719   104.0209   104.0209   104.3698   103.6719   104.0209\n"
        "   104.0209   104.3698   103.8464   104.1954   104.1954   104.5444\n"
        "   108.4935   108.8425   108.8425   109.1915   108.6681   109.0170\n"
        "   109.0170   109.3660   108.6681   109.0170   109.0170   109.3660\n"
        "   108.8426   109.1916   109.1916   109.5406  /\n"
        "\n"
        "PORO\n"
        "  8*0.15 /\n"
        "EDIT\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

static Opm::Deck BAD_CP_GRID_MAPAXES() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "GRID\n"
        "MAPAXES\n"
        " 0.  100.  0.  0.  100.  0.  /\n"
        "\n"
        "SPECGRID\n"
        " 2 2 2 1 F /\n"
        "COORD\n"
        "  2002.0000  2002.0000   100.0000   1999.8255  1999.9127   108.4935\n"
        "  2011.9939  2000.0000   100.3490   2009.8194  1999.9127   108.8425\n"
        "  2015.9878  2000.0000   100.6980   2019.8133  1999.9127   109.1915\n"
        "  2000.0000  2009.9985   100.1745   1999.8255  2009.9112   108.6681 \n"
        "  2010.9939  2011.9985   100.5235   2009.8194  2009.9112   109.0170\n"
        "  2019.9878  2009.9985   100.8725   2019.8133  2009.9112   109.3660\n"
        "  2005.0000  2019.9970   100.3490   1999.8255  2019.9097   108.8426\n"
        "  2009.9939  2019.9970   100.6980   2009.8194  2019.9097   109.1916\n"
        "  2016.9878  2019.9970   101.0470   2019.8133  2019.9097   109.5406 /\n"
        "ZCORN\n"
        "    98.0000   100.3490    97.3490   100.6980   100.1745   100.5235\n"
        "   100.5235   100.8725   100.1745   100.5235   100.5235   100.8725\n"
        "   100.3490   101.6980   101.6980   102.5470   102.4973   102.1463\n"
        "   103.2463   104.1953   103.6719   104.0209   104.0209   104.3698\n"
        "   103.6719   104.0209   104.0209   104.3698   103.8464   104.1954\n"
        "   104.1954   104.5444   103.4973   103.8463   103.8463   104.1953\n"
        "   103.6719   104.0209   104.0209   104.3698   103.6719   104.0209\n"
        "   104.0209   104.3698   103.8464   104.1954   104.1954   104.5444\n"
        "   108.4935   108.8425   108.8425   109.1915   108.6681   109.0170\n"
        "   109.0170   109.3660   108.6681   109.0170   109.0170   109.3660\n"
        "   108.8426   109.1916   109.1916   109.5406  /\n"
        "\n"
        "PORO\n"
        "  8*0.15 /\n"
        "EDIT\n";

    Opm::Parser parser;
    return parser.parseString( deckData);
}

BOOST_AUTO_TEST_CASE(SAVE_FIELD_UNITS) {

    const char* deckData =

        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 4 4 3 /\n"
        "FIELD\n"
        "GRID\n"
        "DX\n"
        " 48*300 /\n"
        "DY\n"
        " 48*300 /\n"
        "DZ\n"
        " 16*20 16*30 16*50 / \n"
        "TOPS\n"
        " 16*8325 / \n"
        "PORO\n"
        "  48*0.15 /\n"
        "EDIT\n"
        "\n";

    const char* deckData2 =

        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 4 4 3 /\n"
        "FIELD\n"
        "GRID\n"
        "MAPUNITS\n"
        " METRES /\n"
        "MAPAXES\n"
        " 0.0  101.1  0.0  0.0  102.2  0.0  /\n"
        "DX\n"
        " 48*300 /\n"
        "DY\n"
        " 48*300 /\n"
        "DZ\n"
        " 16*20 16*30 16*50 / \n"
        "TOPS\n"
        " 16*8325 / \n"
        "PORO\n"
        "  48*0.15 /\n"
        "EDIT\n"
        "\n";

    const char* deckData3 =

        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 4 4 3 /\n"
        "FIELD\n"
        "GRID\n"
        "MAPUNITS\n"
        " FEET /\n"
        "MAPAXES\n"
        " 0.0  102.2  0.0  0.0  103.3  0.0  /\n"
        "DX\n"
        " 48*300 /\n"
        "DY\n"
        " 48*300 /\n"
        "DZ\n"
        " 16*20 16*30 16*50 / \n"
        "TOPS\n"
        " 16*8325 / \n"
        "PORO\n"
        "  48*0.15 /\n"
        "EDIT\n"
        "\n";

    std::vector<float> ref2_mapaxes = {0.0, 101.1, 0.0, 0.0, 102.2, 0.0 };
    std::vector<float> ref3_mapaxes = {0.0, 102.2, 0.0, 0.0, 103.3, 0.0 };

    Opm::Parser parser;
    auto deck = parser.parseString( deckData) ;

    Opm::EclipseState es(deck);
    Opm::UnitSystem units = es.getDeckUnitSystem();
    const auto length = ::Opm::UnitSystem::measure::length;

    const auto& grid1 = es.getInputGrid();

    Opm::NNC nnc(grid1,  deck);
    bool formatted = false;

    time_t timer;
    time(&timer);

    std::string cwd = Opm::filesystem::current_path().c_str();
    std::string testDir = cwd + "/tmp_dir_" + std::to_string(timer);

    if ( Opm::filesystem::exists( testDir )) {
        Opm::filesystem::remove_all(testDir);
    }

    Opm::filesystem::create_directory(testDir);

    std::string fileName = testDir + "/" + "TMP.EGRID";
    grid1.save(fileName, formatted, nnc.input(), units);

    Opm::EclIO::EclFile file1(fileName);

    // Values getZCORNed from the grid needs to be converted from SI to Field units
    // and then converted from double to single precissions before comparing with values saved to
    // the EGRID file

    // check coord
    const std::vector<float> coord_egrid = file1.get<float>("COORD");
    std::vector<double> coord_input_si = grid1.getCOORD();

    BOOST_CHECK( coord_egrid.size() == coord_input_si.size());

    std::vector<float> coord_input_f;
    coord_input_f.reserve(coord_input_si.size());

    for (size_t n =0; n< coord_egrid.size(); n++) {
        coord_input_f.push_back( static_cast<float>(units.from_si(length, coord_input_si[n])));
        BOOST_CHECK_CLOSE( coord_input_f[n] , coord_egrid[n], 1e-6 );
    }

    // check zcorn
    const std::vector<float> zcorn_egrid = file1.get<float>("ZCORN");
    std::vector<double> zcorn_input_si = grid1.getZCORN();

    BOOST_CHECK( zcorn_egrid.size() == zcorn_input_si.size());

    std::vector<float> zcorn_input_f;
    zcorn_input_f.reserve(zcorn_input_si.size());

    for (size_t n =0; n< zcorn_egrid.size(); n++) {
        zcorn_input_f.push_back( static_cast<float>(units.from_si(length, zcorn_input_si[n])));
        BOOST_CHECK_CLOSE( zcorn_input_f[n] , zcorn_egrid[n], 1e-6 );
    }

    BOOST_CHECK( file1.hasKey("GRIDUNIT"));
    const std::vector<std::string> gridunits = file1.get<std::string>("GRIDUNIT");

    BOOST_CHECK( gridunits[0]=="FEET");

    // input deck do not hold MAPAXES or MAPUNITS entries. Below keywords should not be written to EGRID file
    BOOST_CHECK( !file1.hasKey("MAPAXES"));
    BOOST_CHECK( !file1.hasKey("MAPUNITS"));

    // this deck do not have any nnc. Below keywords should not be written to EGRID file
    BOOST_CHECK( !file1.hasKey("NNCHEAD"));
    BOOST_CHECK( !file1.hasKey("NNC1"));
    BOOST_CHECK( !file1.hasKey("NNC2"));

    // testing deck in field units and MAPUNITS in METRES
    auto deck2 = parser.parseString( deckData2) ;

    Opm::EclipseState es2(deck2);
    Opm::UnitSystem units2 = es.getDeckUnitSystem();
    const auto& grid2 = es2.getInputGrid();
    Opm::NNC nnc2(grid2,  deck2 );


    std::string fileName2 = testDir + "/" + "TMP2.FEGRID";

    grid2.save(fileName2, true, nnc2.input(), units);

    Opm::EclIO::EclFile file2(fileName2);

    const std::vector<std::string>& test_mapunits2 = file2.get<std::string>("MAPUNITS");
    BOOST_CHECK( test_mapunits2[0] == "METRES");

    const std::vector<float>& test_mapaxes2 = file2.get<float>("MAPAXES");

    BOOST_CHECK( test_mapaxes2.size() == ref2_mapaxes.size());

    for (size_t n =0; n< ref2_mapaxes.size(); n++) {
        BOOST_CHECK( ref2_mapaxes[n] == test_mapaxes2[n]);
    }

    // testing deck in field units and MAPUNITS in FEET
    auto deck3 = parser.parseString( deckData3) ;

    Opm::EclipseState es3(deck3);
    Opm::UnitSystem units3 = es.getDeckUnitSystem();
    const auto& grid3 = es3.getInputGrid();
    Opm::NNC nnc3(grid3, deck3);


    std::string fileName3 = testDir + "/" + "TMP3.FEGRID";

    grid3.save(fileName3, true, nnc3.input(), units3);

    Opm::EclIO::EclFile file3(fileName3);

    const std::vector<std::string>& test_mapunits3 = file3.get<std::string>("MAPUNITS");
    BOOST_CHECK( test_mapunits3[0] == "FEET");

    const std::vector<float>& test_mapaxes3 = file3.get<float>("MAPAXES");

    BOOST_CHECK( test_mapaxes3.size() == ref3_mapaxes.size());

    for (size_t n =0; n< ref3_mapaxes.size(); n++) {
        BOOST_CHECK( ref3_mapaxes[n] == test_mapaxes3[n]);
    }

    Opm::filesystem::remove_all(testDir);
}

BOOST_AUTO_TEST_CASE(SAVE_METRIC_UNITS) {

    const char* deckData1 =

        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 4 4 3 /\n"
        "GRID\n"
        "MAPAXES\n"
        " 0.0 45000.0 0.0 0.0 720000.0 0.0 / \n"
        "MAPUNITS\n"
        " METRES / \n"
        "DX\n"
        " 48*300 /\n"
        "DY\n"
        " 48*300 /\n"
        "DZ\n"
        " 16*20 16*30 16*50 / \n"
        "TOPS\n"
        " 16*8325 / \n"
        "NNC\n"
        " 2 2 1  2 3 2   0.95 / \n"
        " 3 2 1  3 3 2   1.05 / \n"
        " 4 2 1  4 3 2   1.15 / \n"
        "/ \n"
        "PORO\n"
        "  48*0.15 /\n"
        "EDIT\n"
        "\n";

    const char* deckData2 =

        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 4 4 3 /\n"
        "GRID\n"
        "MAPAXES\n"
        " 0.0 450.0 0.0 0.0 7200.0 0.0 / \n"
        "MAPUNITS\n"
        " FEET / \n"
        "DX\n"
        " 48*300 /\n"
        "DY\n"
        " 48*300 /\n"
        "DZ\n"
        " 16*20 16*30 16*50 / \n"
        "TOPS\n"
        " 16*8325 / \n"
        "NNC\n"
        " 2 2 1  2 3 2   0.95 / \n"
        " 3 2 1  3 3 2   1.05 / \n"
        " 4 2 1  4 3 2   1.15 / \n"
        "/ \n"
        "PORO\n"
        "  48*0.15 /\n"
        "EDIT\n"
        "\n";

    std::vector<float> ref_mapaxes1 = { 0.0, 45000.0, 0.0, 0.0, 720000.0, 0.0 };
    std::vector<float> ref_mapaxes2 = { 0.0, 450.0, 0.0, 0.0, 7200.0, 0.0 };

    Opm::Parser parser;
    auto deck1 = parser.parseString( deckData1) ;

    Opm::EclipseState es1(deck1);
    Opm::UnitSystem units1 = es1.getDeckUnitSystem();
    const auto length = ::Opm::UnitSystem::measure::length;

    const auto& grid1 = es1.getInputGrid();
    Opm::NNC nnc(grid1, deck1);

    bool formatted = true;

    time_t timer;
    time(&timer);

    std::string cwd = Opm::filesystem::current_path().c_str();
    std::string testDir = cwd + "/tmp_dir_" + std::to_string(timer);

    if ( Opm::filesystem::exists( testDir )) {
        Opm::filesystem::remove_all(testDir);
    }

    Opm::filesystem::create_directory(testDir);

    std::string fileName = testDir + "/" + "TMP.FEGRID";
    grid1.save(fileName, formatted, nnc.input(), units1);

    Opm::EclIO::EclFile file1(fileName);

    // Values getZCORNed from the grid have same units as input deck (metric), however these needs to be
    // converted from double to single precissions before comparing with values saved to the EGRID file

    // check coord
    const std::vector<float> coord_egrid = file1.get<float>("COORD");
    std::vector<double> coord_input_si = grid1.getCOORD();

    BOOST_CHECK( coord_egrid.size() == coord_input_si.size());

    std::vector<float> coord_input_f;
    coord_input_f.reserve(coord_input_si.size());

    for (size_t n =0; n< coord_egrid.size(); n++) {
        coord_input_f.push_back( static_cast<float>(units1.from_si(length, coord_input_si[n])));
        BOOST_CHECK_CLOSE( coord_input_f[n] , coord_egrid[n], 1e-6 );
    }

    // check zcorn
    const std::vector<float> zcorn_egrid = file1.get<float>("ZCORN");
    std::vector<double> zcorn_input_si = grid1.getZCORN();

    BOOST_CHECK( zcorn_egrid.size() == zcorn_input_si.size());

    std::vector<float> zcorn_input_f;
    zcorn_input_f.reserve(zcorn_input_si.size());

    for (size_t n =0; n< zcorn_egrid.size(); n++) {
        zcorn_input_f.push_back( static_cast<float>(units1.from_si(length, zcorn_input_si[n])));
        BOOST_CHECK_CLOSE( zcorn_input_f[n] , zcorn_egrid[n], 1e-6 );
    }

    BOOST_CHECK( file1.hasKey("GRIDUNIT"));
    const std::vector<std::string> gridunits = file1.get<std::string>("GRIDUNIT");

    BOOST_CHECK( gridunits[0]=="METRES");

    BOOST_CHECK( file1.hasKey("MAPAXES"));
    std::vector<float> mapaxes = file1.get<float>("MAPAXES");

    for (size_t n = 0; n < 6; n++) {
        BOOST_CHECK_CLOSE( mapaxes[n] , ref_mapaxes1[n], 1e-6 );
    }

    BOOST_CHECK( file1.hasKey("MAPUNITS"));
    const std::vector<std::string> mapunits = file1.get<std::string>("MAPUNITS");
    BOOST_CHECK( gridunits[0]=="METRES");

    BOOST_CHECK( file1.hasKey("NNCHEAD"));
    const std::vector<int> nnchead = file1.get<int>("NNCHEAD");

    BOOST_CHECK( nnchead[0] == static_cast<int>(nnc.input().size()) );

    std::vector<int> ref_nnc1 = { 6, 7, 8 };
    std::vector<int> ref_nnc2 = { 26, 27, 28 };

    BOOST_CHECK( file1.hasKey("NNC1"));
    BOOST_CHECK( file1.hasKey("NNC2"));

    const std::vector<int> nnc1 = file1.get<int>("NNC1");
    const std::vector<int> nnc2 = file1.get<int>("NNC2");

    BOOST_CHECK( nnc1.size() == nnc2.size() );

    for (size_t n =0; n< nnc1.size(); n++) {
        BOOST_CHECK( nnc1[n] == ref_nnc1[n] );
    }

    for (size_t n =0; n< nnc2.size(); n++) {
        BOOST_CHECK( nnc2[n] == ref_nnc2[n] );
    }

    // testing deck in metric units with mapaxes in field units
    auto deck2 = parser.parseString( deckData2) ;

    Opm::EclipseState es2(deck2);
    Opm::UnitSystem units2 = es2.getDeckUnitSystem();

    const auto& grid2 = es2.getInputGrid();
    //Opm::NNC nnc( deck2 );

    std::string fileName2 = testDir + "/" + "TMP2.FEGRID";

    grid2.save(fileName2, true, nnc.input(), units2);

    Opm::EclIO::EclFile file2(fileName2);

    const std::vector<std::string>& test_mapunits2 = file2.get<std::string>("MAPUNITS");
    BOOST_CHECK( test_mapunits2[0] == "FEET");

    const std::vector<float>& test_mapaxes2 = file2.get<float>("MAPAXES");

    BOOST_CHECK( test_mapaxes2.size() == ref_mapaxes2.size());

    for (size_t n =0; n< ref_mapaxes2.size(); n++) {
        BOOST_CHECK( ref_mapaxes2[n] == test_mapaxes2[n]);
    }


    Opm::filesystem::remove_all(testDir);
}

BOOST_AUTO_TEST_CASE(CalcCellDims) {

    Opm::Deck deck = BAD_CP_GRID();
    Opm::EclipseGrid grid( deck );

    std::array<int, 3> dims = grid.getNXYZ();

    size_t nCells = dims[0]*dims[1]*dims[2];

    std::vector<double> dz_ref = { 0.33223500E+01, 0.40973248E+01, 0.32474000E+01, 0.28723750E+01, 0.49961748E+01, 0.49961748E+01,
                                   0.49961748E+01, 0.49961748E+01
                                 };

    std::vector<double> dx_ref = { 0.10309320E+02, 0.70301223E+01, 0.84377403E+01, 0.85725355E+01, 0.10140956E+02, 0.89693098E+01,
                                   0.94102650E+01, 0.94102678E+01
                                 };

    std::vector<double> dy_ref = { 0.99226236E+01, 0.10826077E+02, 0.93370037E+01, 0.93144703E+01, 0.10008223E+02, 0.10302064E+02,
                                   0.97221985E+01, 0.97221985E+01
                                 };

    std::vector<double> depth_ref = { 0.10142293E+03, 0.10190942E+03, 0.10230995E+03, 0.10284644E+03, 0.10625719E+03,
                                      0.10660616E+03, 0.10643174E+03, 0.10678072E+03
                                    };

    for (size_t n=0; n<nCells; n++) {

        BOOST_CHECK_CLOSE( grid.getCellThickness(n) , dz_ref[n], 1e-5 );

        std::array<double, 3> cellDims = grid.getCellDims(n);

        BOOST_CHECK_CLOSE( cellDims[0] , dx_ref[n], 1e-5 );
        BOOST_CHECK_CLOSE( cellDims[1] , dy_ref[n], 1e-5 );
        BOOST_CHECK_CLOSE( cellDims[2] , dz_ref[n], 1e-5 );

        BOOST_CHECK_CLOSE( grid.getCellDepth(n) , depth_ref[n], 1e-5 );
    }

    for (int k = 0; k < dims[2]; k++) {
        for (int j = 0; j < dims[1]; j++) {
            for (int i = 0; i < dims[0]; i++) {
                size_t globInd = i + j*dims[0] + k*dims[0]*dims[1];
                BOOST_CHECK_CLOSE( grid.getCellThickness(i, j, k) , dz_ref[globInd], 1e-5 );

                std::array<double, 3> cellDims = grid.getCellDims(i, j, k);

                BOOST_CHECK_CLOSE( cellDims[0] , dx_ref[globInd], 1e-5 );
                BOOST_CHECK_CLOSE( cellDims[1] , dy_ref[globInd], 1e-5 );
                BOOST_CHECK_CLOSE( cellDims[2] , dz_ref[globInd], 1e-5 );

                BOOST_CHECK_CLOSE( grid.getCellDepth(i, j, k) , depth_ref[globInd], 1e-5 );
            }
        }
    }
}

BOOST_AUTO_TEST_CASE(ExportMAPAXES_TEST) {

    Opm::Deck deck1 = BAD_CP_GRID_MAPAXES();
    Opm::EclipseGrid grid1( deck1 );

    std::vector<double> ref_mapaxes = { 0.0, 100.0, 0.0, 0.0, 100.0, 0.0 };

    std::vector<double> mapaxes = grid1.getMAPAXES();

    for (size_t n=0; n< mapaxes.size(); n++ ) {
        BOOST_CHECK_EQUAL( ref_mapaxes[n] , mapaxes[n]);
    }

    Opm::Deck deck2 = BAD_CP_GRID();
    Opm::EclipseGrid grid2( deck2 );

    BOOST_CHECK( !grid1.equal( grid2 ));

    std::vector<double> coord = grid1.getCOORD();
    std::vector<double> zcorn = grid1.getZCORN();
    std::vector<int> actnum = grid1.getACTNUM();

    std::array<int, 3> dims = grid1.getNXYZ();

    Opm::EclipseGrid grid3(dims, coord, zcorn, actnum.data(), mapaxes.data());

    BOOST_CHECK( grid3.equal( grid1 ));

    mapaxes[1] = 101;
    Opm::EclipseGrid grid4(dims, coord, zcorn, actnum.data(), mapaxes.data());

    BOOST_CHECK( !grid4.equal( grid1 ));
}

BOOST_AUTO_TEST_CASE(TESTCP_ACTNUM_UPDATE) {
    const char* deckData =

        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 3 2 1 /\n"
        "GRID\n"
        "COORD\n"
        " 2000.0000  2000.0000  2000.0000   1999.9476  2000.0000  2002.9995\n"
        " 2049.9924  2000.0000  2000.8726   2049.9400  2000.0000  2003.8722 \n"
        " 2099.9848  2000.0000  2001.7452   2099.9324  2000.0000  2004.7448 \n"
        " 2149.9772  2000.0000  2002.6179   2149.9248  2000.0000  2005.6174 \n"
        " 2000.0000  2050.0000  2000.0000   1999.9476  2050.0000  2002.9995 \n"
        " 2049.9924  2050.0000  2000.8726   2049.9400  2050.0000  2003.8722 \n"
        " 2099.9848  2050.0000  2001.7452   2099.9324  2050.0000  2004.7448 \n"
        " 2149.9772  2050.0000  2002.6179   2149.9248  2050.0000  2005.6174 \n"
        " 2000.0000  2100.0000  2000.0000   1999.9476  2100.0000  2002.9995 \n"
        " 2049.9924  2100.0000  2000.8726   2049.9400  2100.0000  2003.8722 \n"
        " 2099.9848  2100.0000  2001.7452   2099.9324  2100.0000  2004.7448 \n"
        " 2149.9772  2100.0000  2002.6179   2149.9248  2100.0000  2005.6174 / \n"
        "ZCORN\n"
        " 2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179 \n"
        " 2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179 \n"
        " 2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179 \n"
        " 2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179 \n"
        " 2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174 \n"
        " 2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174 \n"
        " 2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174 \n"
        " 2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174 / \n"
        "ACTNUM\n"
        " 0 1 1 1 0 1 / \n"
        "PORO\n"
        "  6*0.15 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    auto deck = parser.parseString( deckData) ;

    std::vector<int> actInDeck = {0, 1, 1, 1, 0, 1};
    std::vector<int> newAct = {1, 0, 0, 0, 1, 0};

    Opm::EclipseGrid grid1( deck);
    Opm::EclipseGrid grid2( deck, newAct.data());

    std::vector<int> actGrid1 = grid1.getACTNUM();
    std::vector<int> actGrid2 = grid2.getACTNUM();

    BOOST_CHECK( actGrid1.size() == actGrid2.size());

    for (size_t n=0; n< actGrid1.size(); n++) {
        BOOST_CHECK_EQUAL( actGrid1[n], actInDeck[n]);
        BOOST_CHECK_EQUAL( actGrid2[n], newAct[n]);
    }
}


BOOST_AUTO_TEST_CASE(TESTCP_ACTNUM_AQUNUM) {
    const std::string deckData = R"(
         RUNSPEC
         DIMENS
         3 2 1 /
         GRID
         COORD
         2000.0000  2000.0000  2000.0000   1999.9476  2000.0000  2002.9995
         2049.9924  2000.0000  2000.8726   2049.9400  2000.0000  2003.8722
         2099.9848  2000.0000  2001.7452   2099.9324  2000.0000  2004.7448
         2149.9772  2000.0000  2002.6179   2149.9248  2000.0000  2005.6174
         2000.0000  2050.0000  2000.0000   1999.9476  2050.0000  2002.9995
         2049.9924  2050.0000  2000.8726   2049.9400  2050.0000  2003.8722
         2099.9848  2050.0000  2001.7452   2099.9324  2050.0000  2004.7448
         2149.9772  2050.0000  2002.6179   2149.9248  2050.0000  2005.6174
         2000.0000  2100.0000  2000.0000   1999.9476  2100.0000  2002.9995
         2049.9924  2100.0000  2000.8726   2049.9400  2100.0000  2003.8722
         2099.9848  2100.0000  2001.7452   2099.9324  2100.0000  2004.7448
         2149.9772  2100.0000  2002.6179   2149.9248  2100.0000  2005.6174 /
         ZCORN
         2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179
         2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179
         2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179
         2000.0000  2000.8726  2000.8726  2001.7452  2001.7452  2002.6179
         2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174
         2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174
         2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174
         2002.9995  2003.8722  2003.8722  2004.7448  2004.7448  2005.6174 /
         ACTNUM
          0 0 1 1 0 1 /
         PORO
          6*0.15 /
         AQUNUM
            1     1  1  1   1000000.0   10000   0.25   400    2585.00   285.00    1   1  /
            1     2  1  1   1000000.0   10000   0.25   400    2585.00   285.00    1   1  /
         /
         AQUCON
           1     3    3      2    2     1   1   'J+'      1.00      1  /
         /
         EDIT)";

    Opm::Parser parser;
    const auto deck = parser.parseString( deckData) ;
    Opm::EclipseGrid grid( deck);

    const std::vector<int>& grid_actnum = grid.getACTNUM();
    const std::vector<int> desired_actnum = {1, 1, 1, 1, 0, 1};

    BOOST_CHECK( grid_actnum.size() == desired_actnum.size() );

    for (size_t n=0; n< grid_actnum.size(); n++) {
        BOOST_CHECK_EQUAL( grid_actnum[n], desired_actnum[n] );
    }

    Opm::EclipseState es(deck);
    const auto& grid_actnum2 = es.getInputGrid().getACTNUM();

    BOOST_CHECK( desired_actnum.size() == grid_actnum2.size());
    for (size_t n=0; n< grid_actnum.size(); n++) {
        BOOST_CHECK_EQUAL( grid_actnum2[n], desired_actnum[n] );
    }
}

BOOST_AUTO_TEST_CASE(TEST_altGridConstructors) {

    const char* deckData =

        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 2 2 3 /\n"
        "GRID\n"
        "SPECGRID\n"
        " 2 2 3 1 F / \n"
        "\n"
        "COORD\n"
        "  2000.0000  2000.0000  1497.0000   1999.9127  2000.0000  1511.9977\n"
        "  2049.9924  2000.0000  1500.8726   2049.9051  2000.0000  1515.8703\n"
        "  2099.9848  2000.0000  1501.7452   2099.8975  2000.0000  1516.7430 \n"
        "  2000.0000  2050.0000  1497.0000   1999.9127  2050.0000  1511.9977\n"
        "  2049.9924  2050.0000  1500.8726   2049.9051  2050.0000  1515.8703\n"
        "  2099.9848  2050.0000  1501.7452   2099.8975  2050.0000  1516.7430\n"
        "  2000.0000  2100.0000  1497.0000   1999.9127  2100.0000  1511.9977\n"
        "  2049.9924  2100.0000  1500.8726   2049.9051  2100.0000  1515.8703\n"
        "  2099.9848  2100.0000  1501.7452   2099.8975  2100.0000  1516.7430 /\n"
        "\n"
        "ZCORN\n"
        "  1497.0000  1497.8726  1500.8726  1501.7452  1497.0000  1497.8726\n"
        "  1500.8726  1501.7452  1497.0000  1497.8726  1500.8726  1501.7452\n"
        "  1497.0000  1497.8726  1500.8726  1501.7452  1501.9992  1502.8719\n"
        "  1505.8719  1506.7445  1501.9992  1502.8719  1505.8719  1506.7445\n"
        "  1501.9992  1502.8719  1505.8719  1506.7445  1501.9992  1502.8719\n"
        "  1505.8719  1506.7445  1501.9992  1502.8719  1505.8719  1506.7445\n"
        "  1501.9992  1502.8719  1505.8719  1506.7445  1501.9992  1502.8719\n"
        "  1505.8719  1506.7445  1501.9992  1502.8719  1505.8719  1506.7445\n"
        "  1506.9985  1507.8711  1510.8711  1511.7437  1506.9985  1507.8711\n"
        "  1510.8711  1511.7437  1506.9985  1507.8711  1510.8711  1511.7437\n"
        "  1506.9985  1507.8711  1510.8711  1511.7437  1506.9985  1507.8711\n"
        "  1510.8711  1511.7437  1506.9985  1507.8711  1510.8711  1511.7437\n"
        "  1506.9985  1507.8711  1510.8711  1511.7437  1506.9985  1507.8711\n"
        "  1510.8711  1511.7437  1511.9977  1512.8703  1515.8703  1516.7430\n"
        "  1511.9977  1512.8703  1515.8703  1516.7430  1511.9977  1512.8703\n"
        "  1515.8703  1516.7430  1511.9977  1512.8703  1515.8703  1516.7430 /\n"
        "\n"
        "ACTNUM\n"
        " 1 1 1 1 1 0 1 1 1 0 1 1  /\n"
        "PORO\n"
        "  24*0.15 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    auto deck = parser.parseString( deckData) ;

    Opm::EclipseGrid grid1( deck);

    std::vector<int> actnum = grid1.getACTNUM();
    std::vector<double> coord = grid1.getCOORD();
    std::vector<double> zcorn = grid1.getZCORN();

    Opm::EclipseGrid grid2( grid1 , zcorn.data(), actnum);
    //Opm::EclipseGrid grid2( grid1 , zcorn, actnum);

    BOOST_CHECK( grid1.equal( grid2) );

    std::vector<double> emptyZcorn;

    Opm::EclipseGrid grid3( grid1 , emptyZcorn.data(), actnum);
    BOOST_CHECK( grid1.equal( grid3) );
}

static Opm::Deck BAD_CP_GRID_ACTNUM() {
    const char* deckData =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "GRID\n"
        "SPECGRID\n"
        " 2 2 2 1 F /\n"
        "COORD\n"
        "  2002.0000  2002.0000   100.0000   1999.8255  1999.9127   108.4935\n"
        "  2011.9939  2000.0000   100.3490   2009.8194  1999.9127   108.8425\n"
        "  2015.9878  2000.0000   100.6980   2019.8133  1999.9127   109.1915\n"
        "  2000.0000  2009.9985   100.1745   1999.8255  2009.9112   108.6681 \n"
        "  2010.9939  2011.9985   100.5235   2009.8194  2009.9112   109.0170\n"
        "  2019.9878  2009.9985   100.8725   2019.8133  2009.9112   109.3660\n"
        "  2005.0000  2019.9970   100.3490   1999.8255  2019.9097   108.8426\n"
        "  2009.9939  2019.9970   100.6980   2009.8194  2019.9097   109.1916\n"
        "  2016.9878  2019.9970   101.0470   2019.8133  2019.9097   109.5406 /\n"
        "ZCORN\n"
        "    98.0000   100.3490    97.3490   100.6980   100.1745   100.5235\n"
        "   100.5235   100.8725   100.1745   100.5235   100.5235   100.8725\n"
        "   100.3490   101.6980   101.6980   102.5470   102.4973   102.1463\n"
        "   103.2463   104.1953   103.6719   104.0209   104.0209   104.3698\n"
        "   103.6719   104.0209   104.0209   104.3698   103.8464   104.1954\n"
        "   104.1954   104.5444   103.4973   103.8463   103.8463   104.1953\n"
        "   103.6719   104.0209   104.0209   104.3698   103.6719   104.0209\n"
        "   104.0209   104.3698   103.8464   104.1954   104.1954   104.5444\n"
        "   108.4935   108.8425   108.8425   109.1915   108.6681   109.0170\n"
        "   109.0170   109.3660   108.6681   109.0170   109.0170   109.3660\n"
        "   108.8426   109.1916   109.1916   109.5406  /\n"
        "\n"
        "ACTNUM\n"
        "   1 1 1 1 0 1 0 1  /\n"
        "PORO\n"
        "  8*0.15 /\n"
        "EDIT\n";

    Opm::Parser parser;
    return parser.parseString( deckData );
}

BOOST_AUTO_TEST_CASE(TEST_getCellCenters) {

    Opm::Deck deck1 = BAD_CP_GRID_ACTNUM();
    Opm::EclipseGrid grid1( deck1 );

    std::vector<std::array<double, 3>> ref_centers_prev = {
        { 2.006104082026e+03, 2.005869733884e+03, 1.014229250000e+02 },
        { 2.014871600974e+03, 2.005382958337e+03, 1.019094125000e+02 },
        { 2.006149627699e+03, 2.015375548018e+03, 1.023099500000e+02 },
        { 2.014617670881e+03, 2.015373620907e+03, 1.028464375000e+02 },
        { 2.014794142285e+03, 2.005084683036e+03, 1.066061625000e+02 },
        { 2.014720614821e+03, 2.015083182885e+03, 1.067807125000e+02 }
    };

    std::vector<std::array<double, 3>> ref_dims_prev = {
        { 1.030932076324e+01, 9.922624077310e+00, 3.322350000000e+00 },
        { 7.030122275667e+00, 1.082607766242e+01, 4.097325000000e+00 },
        { 8.437740561450e+00, 9.337003646693e+00, 3.247400000000e+00 },
        { 8.572535163904e+00, 9.314470420138e+00, 2.872375000000e+00 },
        { 8.969310273018e+00, 1.030206365401e+01, 4.996175000000e+00 },
        { 9.410267880374e+00, 9.722198202980e+00, 4.996175000000e+00 }
    };

    std::vector<int> actMap = grid1.getActiveMap();

    int n = 0;
    for (auto ind : actMap) {
        std::array<double, 3> cellC = grid1.getCellCenter(ind);
        std::array<double, 3> cellD = grid1.getCellDims(ind);

        for (size_t i = 0; i < 3; i++) {
            BOOST_CHECK_CLOSE( ref_centers_prev[n][i], cellC[i], 1e-10);
            BOOST_CHECK_CLOSE( ref_dims_prev[n][i], cellD[i], 1e-10);
        }

        n++;
    }
}

BOOST_AUTO_TEST_CASE(LoadFromBinary) {
    BOOST_CHECK_THROW(Opm::EclipseGrid( "No/does/not/exist" ) , std::runtime_error);
}


BOOST_AUTO_TEST_CASE(TEST_constructFromEgrid) {

    const char* deckData =

        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        " 10 10 3 /\n"
        "FIELD\n"
        "GRID\n"
        "DX\n"
        " 300*1000 /\n"
        "DY\n"
        " 300*1000 /\n"
        "DZ\n"
        " 100*20 100*30 100*50 / \n"
        "TOPS\n"
        " 100*8325 / \n"
        "ACTNUM\n"
        " 44*1 3*0 7*1 3*0  243*1/\n"
        "PORO\n"
        "  300*0.15 /\n"
        "EDIT\n"
        "\n";

    Opm::Parser parser;
    auto deck = parser.parseString( deckData) ;

    Opm::EclipseGrid grid1( deck);
    Opm::EclipseGrid grid2( "SPE1CASE1.EGRID");

    // compare actnum
    std::vector<int> actGrid1 = grid1.getACTNUM();
    std::vector<int> actGrid2 = grid2.getACTNUM();

    BOOST_CHECK( actGrid1.size() == actGrid2.size() );

    for (size_t n= 0; n< actGrid1.size(); n++) {
        BOOST_CHECK( actGrid1[n] == actGrid2[n] );
    }

    // compare coord
    std::vector<double> coordGrid1 = grid1.getCOORD();
    std::vector<double> coordGrid2 = grid2.getCOORD();

    BOOST_CHECK( coordGrid1.size() == coordGrid2.size() );

    std::vector<double> zcornGrid1 = grid1.getZCORN();
    std::vector<double> zcornGrid2 = grid2.getZCORN();

    BOOST_CHECK( zcornGrid1.size() == zcornGrid2.size() );

    for (size_t n = 0; n < zcornGrid1.size(); n++){
        BOOST_CHECK_CLOSE( zcornGrid1[n], zcornGrid2[n], 1.0e-6 );
    }

    BOOST_CHECK( grid1.getCartesianSize() == grid2.getCartesianSize() );

    for (size_t n=0; n < grid1.getCartesianSize(); n++){
        BOOST_CHECK_CLOSE( grid1.getCellVolume(n), grid2.getCellVolume(n), 1e-6);

    }
}



BOOST_AUTO_TEST_CASE(TEST_GDFILE_1) {

    const char* deckData1 =
        "RUNSPEC\n"
        "DIMENS\n"
        "1 1 2 /\n"
        "GRID\n"
        "COORD\n"
        "10.0000    10.0000  2000.0000      9.8255    10.0000  2014.9977\n"
        "109.9848   10.0000  2001.7452    109.8102    10.0000  2016.7430\n"
        "10.0000   110.0000  2000.0000      9.8255   110.0000  2014.9977\n"
        "109.9848   110.0000  2001.7452    109.8102   110.0000  2016.7430 /\n"
        "ZCORN\n"
        "2000.0000  2001.7452  2000.0000  2001.7452  2004.9992  2006.7445\n"
        "2004.9992  2006.7445  2004.9992  2006.7445  2004.9992  2006.7445\n"
        "2014.9977  2016.7430  2014.9977  2016.7430 /\n"
        "PORO\n"
        "   2*0.15 /\n";

    Opm::Parser parser;
    auto deck1 = parser.parseString( deckData1) ;
    BOOST_CHECK_NO_THROW( Opm::EclipseGrid grid1(deck1) );
}


BOOST_AUTO_TEST_CASE(TEST_GDFILE_2) {

    const char* deckData1 =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "GRID\n"
        "SPECGRID\n"
        " 2 2 2 1 F /\n"
        "COORD\n"
        "  2002.0000  2002.0000   100.0000   1999.8255  1999.9127   108.4935\n"
        "  2011.9939  2000.0000   100.3490   2009.8194  1999.9127   108.8425\n"
        "  2015.9878  2000.0000   100.6980   2019.8133  1999.9127   109.1915\n"
        "  2000.0000  2009.9985   100.1745   1999.8255  2009.9112   108.6681 \n"
        "  2010.9939  2011.9985   100.5235   2009.8194  2009.9112   109.0170\n"
        "  2019.9878  2009.9985   100.8725   2019.8133  2009.9112   109.3660\n"
        "  2005.0000  2019.9970   100.3490   1999.8255  2019.9097   108.8426\n"
        "  2009.9939  2019.9970   100.6980   2009.8194  2019.9097   109.1916\n"
        "  2016.9878  2019.9970   101.0470   2019.8133  2019.9097   109.5406 /\n"
        "ZCORN\n"
        "    98.0000   100.3490    97.3490   100.6980   100.1745   100.5235\n"
        "   100.5235   100.8725   100.1745   100.5235   100.5235   100.8725\n"
        "   100.3490   101.6980   101.6980   102.5470   102.4973   102.1463\n"
        "   103.2463   104.1953   103.6719   104.0209   104.0209   104.3698\n"
        "   103.6719   104.0209   104.0209   104.3698   103.8464   104.1954\n"
        "   104.1954   104.5444   103.4973   103.8463   103.8463   104.1953\n"
        "   103.6719   104.0209   104.0209   104.3698   103.6719   104.0209\n"
        "   104.0209   104.3698   103.8464   104.1954   104.1954   104.5444\n"
        "   108.4935   108.8425   108.8425   109.1915   108.6681   109.0170\n"
        "   109.0170   109.3660   108.6681   109.0170   109.0170   109.3660\n"
        "   108.8426   109.1916   109.1916   109.5406  /\n"
        "\n"
        "PORO\n"
        "   8*0.15 /\n"
        "EDIT\n";

    const char* deckData1a =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "GRID\n"
        "SPECGRID\n"
        " 2 2 2 1 F /\n"
        "COORD\n"
        "  2002.0000  2002.0000   100.0000   1999.8255  1999.9127   108.4935\n"
        "  2011.9939  2000.0000   100.3490   2009.8194  1999.9127   108.8425\n"
        "  2015.9878  2000.0000   100.6980   2019.8133  1999.9127   109.1915\n"
        "  2000.0000  2009.9985   100.1745   1999.8255  2009.9112   108.6681 \n"
        "  2010.9939  2011.9985   100.5235   2009.8194  2009.9112   109.0170\n"
        "  2019.9878  2009.9985   100.8725   2019.8133  2009.9112   109.3660\n"
        "  2005.0000  2019.9970   100.3490   1999.8255  2019.9097   108.8426\n"
        "  2009.9939  2019.9970   100.6980   2009.8194  2019.9097   109.1916\n"
        "  2016.9878  2019.9970   101.0470   2019.8133  2019.9097   109.5406 /\n"
        "ZCORN\n"
        "    98.0000   100.3490    97.3490   100.6980   100.1745   100.5235\n"
        "   100.5235   100.8725   100.1745   100.5235   100.5235   100.8725\n"
        "   100.3490   101.6980   101.6980   102.5470   102.4973   102.1463\n"
        "   103.2463   104.1953   103.6719   104.0209   104.0209   104.3698\n"
        "   103.6719   104.0209   104.0209   104.3698   103.8464   104.1954\n"
        "   104.1954   104.5444   103.4973   103.8463   103.8463   104.1953\n"
        "   103.6719   104.0209   104.0209   104.3698   103.6719   104.0209\n"
        "   104.0209   104.3698   103.8464   104.1954   104.1954   104.5444\n"
        "   108.4935   108.8425   108.8425   109.1915   108.6681   109.0170\n"
        "   109.0170   109.3660   108.6681   109.0170   109.0170   109.3660\n"
        "   108.8426   109.1916   109.1916   109.5406  /\n"
        "\n"
        "ACTNUM\n"
        " 1 1 1 1 0 1 0 1 /\n"
        "PORO\n"
        "   8*0.15 /\n"
        "EDIT\n";

    const char* deckData1b =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "FIELD\n"
        "GRID\n"
        "MAPUNITS\n"
        " METRES /\n"
        "MAPAXES\n"
        " 0.  100.  0.  0.  100.  0.  /\n"
        "SPECGRID\n"
        " 2 2 2 1 F /\n"
        "COORD\n"
        "  2002.0000  2002.0000   100.0000   1999.8255  1999.9127   108.4935\n"
        "  2011.9939  2000.0000   100.3490   2009.8194  1999.9127   108.8425\n"
        "  2015.9878  2000.0000   100.6980   2019.8133  1999.9127   109.1915\n"
        "  2000.0000  2009.9985   100.1745   1999.8255  2009.9112   108.6681 \n"
        "  2010.9939  2011.9985   100.5235   2009.8194  2009.9112   109.0170\n"
        "  2019.9878  2009.9985   100.8725   2019.8133  2009.9112   109.3660\n"
        "  2005.0000  2019.9970   100.3490   1999.8255  2019.9097   108.8426\n"
        "  2009.9939  2019.9970   100.6980   2009.8194  2019.9097   109.1916\n"
        "  2016.9878  2019.9970   101.0470   2019.8133  2019.9097   109.5406 /\n"
        "ZCORN\n"
        "    98.0000   100.3490    97.3490   100.6980   100.1745   100.5235\n"
        "   100.5235   100.8725   100.1745   100.5235   100.5235   100.8725\n"
        "   100.3490   101.6980   101.6980   102.5470   102.4973   102.1463\n"
        "   103.2463   104.1953   103.6719   104.0209   104.0209   104.3698\n"
        "   103.6719   104.0209   104.0209   104.3698   103.8464   104.1954\n"
        "   104.1954   104.5444   103.4973   103.8463   103.8463   104.1953\n"
        "   103.6719   104.0209   104.0209   104.3698   103.6719   104.0209\n"
        "   104.0209   104.3698   103.8464   104.1954   104.1954   104.5444\n"
        "   108.4935   108.8425   108.8425   109.1915   108.6681   109.0170\n"
        "   109.0170   109.3660   108.6681   109.0170   109.0170   109.3660\n"
        "   108.8426   109.1916   109.1916   109.5406  /\n"
        "\n"
        "ACTNUM\n"
        " 1 1 1 1 0 1 0 1 /\n"
        "PORO\n"
        "   8*0.15 /\n"
        "EDIT\n";

    const char* deckData2 =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "GRID\n"
        "GDFILE\n"
        " 'BAD_CP_M.EGRID' /\n"
        "EDIT\n";

    const char* deckData3a =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "GRID\n"
        "ACTNUM\n"
        " 1 0 1 0 1 1 1 1 /\n"
        "MAPUNITS\n"
        " FEET /\n"
        "MAPAXES\n"
        " 0.  200.  0.  0.  200.  0.  /\n"
        "GDFILE\n"
        " 'BAD_CP_M.EGRID' /\n"
        "EDIT\n";

    const char* deckData3b =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "GRID\n"
        "MAPUNITS\n"
        " FEET /\n"
        "MAPAXES\n"
        " 0.  200.  0.  0.  200.  0.  /\n"
        "GDFILE\n"
        " 'BAD_CP_F.EGRID' /\n"
        "ACTNUM\n"
        " 1 0 1 0 1 1 1 1 /\n"
        "EDIT\n";

    const char* deckData3c =
        "RUNSPEC\n"
        "\n"
        "DIMENS\n"
        "2 2 2 /\n"
        "GRID\n"
        "GDFILE\n"
        " 'BAD_CP_F.EGRID' /\n"
        "MAPUNITS\n"
        " FEET /\n"
        "MAPAXES\n"
        " 0.  200.  0.  0.  200.  0.  /\n"
        "ACTNUM\n"
        " 1 0 1 0 1 1 1 1 /\n"
        "EDIT\n";

    Opm::Parser parser;

    std::vector<int> ref_act_egrid = {1, 1, 1, 1, 0, 1, 0, 1};
    std::vector<int> ref_act_deck3 = {1, 0, 1, 0, 1, 1, 1, 1};

    std::vector<double> ref_mapaxes_egrid = { 0.0, 100.0, 0.0, 0.0, 100.0, 0.0 };
    std::vector<double> ref_mapaxes_deck = { 0.0, 200.0, 0.0, 0.0, 200.0, 0.0 };


    // egrid file in si units, no conversion requied by grid constructor
    std::vector<double> refDepthGrid3a = {101.42292, 101.90941, 102.30995, 102.84644, 106.25719, 106.60616, 106.43174, 106.78071 };

    // egrid file in field units, depths converted to SI units when loaded by grid constructor
    std::vector<double> refDepthGrid3b = {30.913707, 31.061988, 31.184072, 31.347594, 32.38719, 32.493558, 32.440393 };

    auto deck1a = parser.parseString( deckData1a) ;

    Opm::EclipseState es1a( deck1a );
    Opm::UnitSystem units1a = es1a.getDeckUnitSystem();

    const auto& grid1a = es1a.getInputGrid();
    Opm::NNC nnc(grid1a, deck1a);

    grid1a.save("BAD_CP_M.EGRID", false, nnc.input(), units1a);

    auto deck1b = parser.parseString( deckData1b) ;
    Opm::EclipseState es1b( deck1b );
    Opm::UnitSystem units1b = es1b.getDeckUnitSystem();
    const auto& grid1b = es1b.getInputGrid();

    grid1b.save("BAD_CP_F.EGRID", false, nnc.input(), units1b);

    auto deck1 = parser.parseString( deckData1) ;
    Opm::EclipseGrid grid1( deck1);

    Opm::EclIO::EclFile file1("BAD_CP_M.EGRID");

    // actnum not defined in deck. keyword GDFILE not present in the DECK
    // check that coord and zcorn from deck-grid identical to coord and zcorn
    // from egrid - grid

    std::vector<double> coordGrid1 = grid1.getCOORD();
    std::vector<double> zcornGrid1 = grid1.getZCORN();
    std::vector<float> coordGrid1_f(coordGrid1.begin(), coordGrid1.end() );
    std::vector<float> zcornGrid1_f(zcornGrid1.begin(), zcornGrid1.end() );

    const std::vector<float> coord_egrid_f = file1.get<float>("COORD");
    const std::vector<float> zcorn_egrid_f = file1.get<float>("ZCORN");

    BOOST_CHECK( coordGrid1.size() == coord_egrid_f.size() );
    BOOST_CHECK( zcornGrid1.size() == zcorn_egrid_f.size() );

    for (size_t n = 0; n < coordGrid1.size(); n++){
        BOOST_CHECK( coordGrid1_f[n] == coord_egrid_f[n] );
    }

    for (size_t n = 0; n < zcornGrid1.size(); n++){
        BOOST_CHECK( zcornGrid1_f[n] == zcorn_egrid_f[n] );
    }

    // all cells are active, since ACTNUM not present
    std::vector<int> actGrid1 = grid1.getACTNUM();
    for (size_t n = 0; n < actGrid1.size(); n++){
        BOOST_CHECK( actGrid1[n] == 1 );
    }

    BOOST_CHECK( grid1.getMAPUNITS() == "" );

    std::vector<double> mapaxes = grid1.getMAPAXES();
    BOOST_CHECK( mapaxes.size() == 0 );


    auto deck2 = parser.parseString( deckData2) ;
    Opm::EclipseGrid grid2( deck2);

    std::vector<int> actGrid2 = grid2.getACTNUM();

    // check that actnum is reset from gdfile

    for (size_t n = 0; n < actGrid2.size(); n++){
        BOOST_CHECK( actGrid2[n] == ref_act_egrid[n] );
    }

    BOOST_CHECK( grid2.getMAPUNITS() == "" );

    mapaxes = grid2.getMAPAXES();
    BOOST_CHECK( mapaxes.size() == 0 );


    auto deck3a = parser.parseString( deckData3a) ;
    Opm::EclipseGrid grid3a( deck3a);

    // mapunits and mapaxes define in deck (only)

    BOOST_CHECK( grid3a.getMAPUNITS() == "FEET" );

    mapaxes = grid3a.getMAPAXES();
    BOOST_CHECK( mapaxes.size() == 6 );

    for (size_t n = 0; n < mapaxes.size(); n++){
        BOOST_CHECK( mapaxes[n] == ref_mapaxes_deck[n] );
    }

    std::vector<int> actGrid3 = grid3a.getACTNUM();

    // check that actnum is reset from gdfile, ACTNUM input in deck
    // but before keyword GDFILE

    for (size_t n = 0; n < actGrid3.size(); n++){
        BOOST_CHECK( actGrid3[n] == ref_act_egrid[n] );
    }

    // check that depth values are in SI units
    for (size_t n = 0; n < refDepthGrid3a.size(); n++){
        BOOST_CHECK_CLOSE( grid3a.getCellDepth(n), refDepthGrid3a[n], 1e-3 );
    }

    auto deck3b = parser.parseString( deckData3b) ;
    Opm::EclipseGrid grid3b( deck3b);

    // mapunits and mapaxes both in egrid and deck. Uses properties
    // from the egrid keyword gdfile input after MAPUNITS and MAPAXES

    BOOST_CHECK( grid3b.getMAPUNITS() == "METRES" );

    mapaxes = grid3b.getMAPAXES();
    BOOST_CHECK( mapaxes.size() == 6 );

    for (size_t n = 0; n < mapaxes.size(); n++){
        BOOST_CHECK( mapaxes[n] == ref_mapaxes_egrid[n] );
    }

    actGrid3 = grid3b.getACTNUM();

    // check that actnum is reset from deck since input after keyword GDFILE
    for (size_t n = 0; n < actGrid3.size(); n++){
        BOOST_CHECK( actGrid3[n] == ref_act_deck3[n] );
    }

    // check that depth values are converted from Field to SI units
    for (size_t n = 0; n < refDepthGrid3b.size(); n++){
        BOOST_CHECK_CLOSE( grid3b.getCellDepth(n), refDepthGrid3b[n], 1e-3 );
    }

    // mapunits and mapaxes both in egrid and deck. Uses properties
    // from the deck sinze these are input after GDfile

    auto deck3c = parser.parseString( deckData3c) ;
    Opm::EclipseGrid grid3c( deck3c);

    BOOST_CHECK( grid3c.getMAPUNITS() == "FEET" );

    mapaxes = grid3c.getMAPAXES();
    BOOST_CHECK( mapaxes.size() == 6 );

    for (size_t n = 0; n < mapaxes.size(); n++){
        BOOST_CHECK( mapaxes[n] == ref_mapaxes_deck[n] );
    }
}

BOOST_AUTO_TEST_CASE(TEST_COLLAPSED_CELL) {
    Opm::EclipseGrid grid(2,2,2,1,1,0);
    for (std::size_t g = 0; g < grid.getCartesianSize(); g++)
        BOOST_CHECK_EQUAL(grid.getCellVolume(g), 0);
}