/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2025-     Equinor ASA
//
//  ResInsight 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.
//
//  ResInsight 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 at <http://www.gnu.org/licenses/gpl.html>
//  for more details.
//
/////////////////////////////////////////////////////////////////////////////////

#include "RimCornerPointCase.h"

#include "RiaApplication.h"
#include "RiaDefines.h"
#include "RiaLogging.h"
#include "RiaPreferencesGrid.h"
#include "RiaQStringFormatter.h"

#include "RifInputPropertyLoader.h"
#include "RifRoffFileTools.h"

#include "RigActiveCellInfo.h"
#include "RigCaseCellResultsData.h"
#include "RigEclipseCaseData.h"
#include "RigEclipseResultAddress.h"
#include "RigMainGrid.h"
#include "RimReloadCaseTools.h"

#include "RimEclipseInputProperty.h"
#include "RimReservoirCellResultsStorage.h"

#include "cafPdmObjectScriptingCapability.h"
#include "cafProgressInfo.h"

#include <QDir>
#include <QFileInfo>

#include <algorithm>
#include <chrono>

#ifdef USE_OPENMP
#include <omp.h>
#endif

using namespace std::chrono;

CAF_PDM_SOURCE_INIT( RimCornerPointCase, "CornerPointCase", "RimCornerPointCase" );
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimCornerPointCase::RimCornerPointCase()
    : RimEclipseCase()
{
    CAF_PDM_InitScriptableObject( "RimCornerPointCase", ":/EclipseInput48x48.png" );
    setReservoirData( new RigEclipseCaseData( this ) );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimCornerPointCase::~RimCornerPointCase()
{
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimCornerPointCase::openEclipseGridFile()
{
    return true;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<RimCornerPointCase*, QString> RimCornerPointCase::createFromCoordinatesArray( const int                 nx,
                                                                                            const int                 ny,
                                                                                            const int                 nz,
                                                                                            const std::vector<float>& coord,
                                                                                            const std::vector<float>& zcorn,
                                                                                            const std::vector<float>& actnum )
{
    CAF_ASSERT( nx > 0 );
    CAF_ASSERT( ny > 0 );
    CAF_ASSERT( nz > 0 );

    size_t ncoord = ( nx + 1 ) * ( ny + 1 ) * 2 * 3;
    size_t nzcorn = nx * ny * nz * 8;
    size_t ntot   = nx * ny * nz;

    if ( coord.size() != ncoord )
        return std::unexpected( QString( "Wrong size of coord array. Expected %1, but got %2" ).arg( ncoord ).arg( coord.size() ) );

    if ( zcorn.size() != nzcorn )
        return std::unexpected( QString( "Wrong size of zcorn array. Expected %1, but got %2" ).arg( nzcorn ).arg( zcorn.size() ) );

    if ( actnum.size() != ntot )
        return std::unexpected( QString( "Wrong size of actnum array. Expected %1, but got %2" ).arg( ntot ).arg( actnum.size() ) );

    auto cornerPointCase = new RimCornerPointCase;

    buildGrid( *cornerPointCase->eclipseCaseData(), nx, ny, nz, coord, zcorn, actnum );
    cornerPointCase->computeCachedData();
    computeDepthRelatedResults( *cornerPointCase );

    return cornerPointCase;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<void, QString> RimCornerPointCase::replaceGridFromCoordinatesArray( RimCornerPointCase&       cornerPointCase,
                                                                                  const int                 nx,
                                                                                  const int                 ny,
                                                                                  const int                 nz,
                                                                                  const std::vector<float>& coord,
                                                                                  const std::vector<float>& zcorn,
                                                                                  const std::vector<float>& actnum )
{
    RigActiveCellInfo* activeCellInfo = cornerPointCase.eclipseCaseData()->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL );
    CAF_ASSERT( activeCellInfo );
    activeCellInfo->clear();

    RigActiveCellInfo* fractureActiveCellInfo =
        cornerPointCase.eclipseCaseData()->activeCellInfo( RiaDefines::PorosityModelType::FRACTURE_MODEL );
    CAF_ASSERT( fractureActiveCellInfo );
    fractureActiveCellInfo->clear();

    RimReloadCaseTools::clearAllGridData( cornerPointCase.eclipseCaseData() );

    // Clear the existing grid geometry before building the new grid
    RigMainGrid* mainGrid = cornerPointCase.eclipseCaseData()->mainGrid();
    CAF_ASSERT( mainGrid );
    mainGrid->reservoirCells().clear();
    mainGrid->nodes().clear();

    buildGrid( *cornerPointCase.eclipseCaseData(), nx, ny, nz, coord, zcorn, actnum );
    cornerPointCase.computeCachedData();
    computeDepthRelatedResults( cornerPointCase );

    return {};
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::array<cvf::Vec3d, 8> RimCornerPointCase::getCorners( const RigMainGrid&        grid,
                                                          const std::vector<float>& coord,
                                                          const std::vector<float>& zcorn,
                                                          const size_t              cellIdx,
                                                          const cvf::Vec3d&         offset,
                                                          const cvf::Vec3d&         scale )
{
    size_t i;
    size_t j;
    size_t k;
    grid.ijkFromCellIndex( cellIdx, &i, &j, &k );

    const size_t nx = grid.cellCountI();
    const size_t ny = grid.cellCountJ();

    // Get depths from zcorn
    std::array<size_t, 8> zind;
    zind[0] = ( k * nx * ny * 8 + j * nx * 4 + i * 2 );
    zind[1] = ( zind[0] + 1 );
    zind[2] = ( zind[0] + nx * 2 );
    zind[3] = ( zind[2] + 1 );

    for ( size_t n = 0; n < 4; n++ )
        zind[n + 4] = ( zind[n] + nx * ny * 4 );

    std::array<cvf::Vec3d, 8> corners;
    for ( size_t n = 0; n < 8; n++ )
        corners[n].z() = zcorn[zind[n]];

    // calculate indices for grid pillars in COORD array
    std::array<size_t, 4> pind;
    pind[0] = j * ( nx + 1 ) * 6 + i * 6;
    pind[1] = pind[0] + 6;
    pind[2] = pind[0] + ( nx + 1 ) * 6;
    pind[3] = pind[2] + 6;

    for ( size_t n = 0; n < 4; n++ )
    {
        const double zt = coord[pind[n] + 2];
        const double zb = coord[pind[n] + 5];

        const double xt = coord[pind[n]];
        const double yt = coord[pind[n] + 1];
        const double xb = coord[pind[n] + 3];
        const double yb = coord[pind[n] + 4];

        const double diffZ = zt - zb;

        if ( diffZ == 0.0 )
        {
            corners[n].x()     = xt;
            corners[n + 4].x() = xt;
            corners[n].y()     = yt;
            corners[n + 4].y() = yt;
        }
        else
        {
            corners[n].x()     = xt + ( xb - xt ) / diffZ * ( zt - corners[n].z() );
            corners[n + 4].x() = xt + ( xb - xt ) / diffZ * ( zt - corners[n + 4].z() );
            corners[n].y()     = yt + ( yb - yt ) / diffZ * ( zt - corners[n].z() );
            corners[n + 4].y() = yt + ( yb - yt ) / diffZ * ( zt - corners[n + 4].z() );
        }
    }

    return corners;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimCornerPointCase::buildGrid( RigEclipseCaseData&       eclipseCaseData,
                                    const int                 nx,
                                    const int                 ny,
                                    const int                 nz,
                                    const std::vector<float>& coord,
                                    const std::vector<float>& zcorn,
                                    const std::vector<float>& actnum )
{
    auto startTime = high_resolution_clock::now();

    RigActiveCellInfo* activeCellInfo = eclipseCaseData.activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL );
    CAF_ASSERT( activeCellInfo );

    RigActiveCellInfo* fractureActiveCellInfo = eclipseCaseData.activeCellInfo( RiaDefines::PorosityModelType::FRACTURE_MODEL );
    CAF_ASSERT( fractureActiveCellInfo );

    RigMainGrid* mainGrid = eclipseCaseData.mainGrid();
    CAF_ASSERT( mainGrid );

    mainGrid->setCellCounts( cvf::Vec3st( nx, ny, nz ) );

    mainGrid->setGridName( "Main grid" );

    size_t totalCellCount = nx * ny * nz;

    activeCellInfo->setGridCount( 1 );
    fractureActiveCellInfo->setGridCount( 1 );

    activeCellInfo->setReservoirCellCount( totalCellCount );
    fractureActiveCellInfo->setReservoirCellCount( totalCellCount );

    // Reserve room for the cells and nodes and fill them with data
    mainGrid->reservoirCells().reserve( totalCellCount );
    mainGrid->nodes().reserve( 8 * totalCellCount );

    int               progTicks = 100;
    caf::ProgressInfo progInfo( progTicks, "" );

    int    cellCount      = static_cast<int>( totalCellCount );
    size_t cellStartIndex = mainGrid->reservoirCells().size();
    size_t nodeStartIndex = mainGrid->nodes().size();

    RigCell defaultCell;
    defaultCell.setHostGrid( mainGrid );
    mainGrid->reservoirCells().resize( cellStartIndex + cellCount, defaultCell );

    mainGrid->nodes().resize( nodeStartIndex + static_cast<size_t>( cellCount ) * 8, cvf::Vec3d( 0, 0, 0 ) );

    const size_t cellMappingECLRi[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };

    cvf::Vec3d offset( 0.0, 0.0, 0.0 );
    cvf::Vec3d scale( 1.0, 1.0, 1.0 );

    // Convert actnum to integer
    std::vector<int> activeCells( nx * ny * nz, 0 );
    CAF_ASSERT( activeCells.size() == actnum.size() );
    for ( size_t i = 0; i < activeCells.size(); i++ )
    {
        activeCells[i] = static_cast<int>( actnum[i] > 0.0 );
    }

    // Precompute the active cell matrix index
    size_t numActiveCells = RifRoffFileTools::computeActiveCellMatrixIndex( activeCells );

    // Loop over cells and fill them with data
#pragma omp parallel for
    for ( int gridLocalCellIndex = 0; gridLocalCellIndex < cellCount; ++gridLocalCellIndex )
    {
        RigCell& cell = mainGrid->cell( cellStartIndex + gridLocalCellIndex );

        cell.setGridLocalCellIndex( gridLocalCellIndex );

        // Active cell index
        int matrixActiveIndex = activeCells[gridLocalCellIndex];
        if ( matrixActiveIndex != -1 )
        {
            activeCellInfo->setCellResultIndex( ReservoirCellIndex( cellStartIndex + gridLocalCellIndex ),
                                                ActiveCellIndex( matrixActiveIndex ) );
        }

        cell.setParentCellIndex( cvf::UNDEFINED_SIZE_T );

        std::array<cvf::Vec3d, 8> corners = getCorners( *mainGrid, coord, zcorn, gridLocalCellIndex, offset, scale );

        // Corner coordinates
        for ( int cIdx = 0; cIdx < 8; ++cIdx )
        {
            double* point = mainGrid->nodes()[nodeStartIndex + (size_t)gridLocalCellIndex * 8 + cellMappingECLRi[cIdx]].ptr();

            point[0] = corners[cIdx].x();
            point[1] = corners[cIdx].y();
            point[2] = -corners[cIdx].z();

            cell.cornerIndices()[cIdx] = nodeStartIndex + (size_t)gridLocalCellIndex * 8 + cIdx;
        }

        // Mark inactive long pyramid looking cells as invalid
        cell.setInvalid( cell.isLongPyramidCell() );
    }

    activeCellInfo->setGridActiveCellCounts( 0, numActiveCells );
    fractureActiveCellInfo->setGridActiveCellCounts( 0, 0 );

    mainGrid->initAllSubGridsParentGridPointer();
    activeCellInfo->computeDerivedData();
    fractureActiveCellInfo->computeDerivedData();

    createActnumResult( eclipseCaseData );

    auto endTime = high_resolution_clock::now();

    auto totalDuration = duration_cast<milliseconds>( endTime - startTime );
    RiaLogging::info( std::format( "Total: {} ms", totalDuration.count() ) );
}

//--------------------------------------------------------------------------------------------------
/// Materialize an ACTNUM static result so the created grid exposes the same queryable property as a
/// grid loaded from an Eclipse file. The corner point case has no reader interface, so the data
/// cannot be loaded lazily and must be created here.
///
/// Why every value is 1.0: a STATIC_NATIVE result is stored in the per-active-cell layout, i.e. the
/// array has exactly one slot per active cell, indexed by the active-cell result index. Inactive
/// cells (input actnum <= 0) are already excluded upstream in buildGrid(): they receive no result
/// index and are not counted by reservoirActiveCellCount(). Inactivity is therefore encoded by the
/// absence of a slot, not by storing a 0 here, so every slot that exists corresponds to an active
/// cell and must hold 1.0. This is the same thing the native readers do
/// (RifReaderEclipseOutput/RifReaderOpmCommon::staticResult: resize(reservoirActiveCellCount(), 1.0)).
/// Walking RigActiveCellInfo to set values per cell would only be needed for a per-reservoir-cell
/// layout, where inactive cells need an explicit 0 — that is not the layout used for ACTNUM.
//--------------------------------------------------------------------------------------------------
void RimCornerPointCase::createActnumResult( RigEclipseCaseData& eclipseCaseData )
{
    RigActiveCellInfo* activeCellInfo = eclipseCaseData.activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL );
    CAF_ASSERT( activeCellInfo );

    RigCaseCellResultsData* matrixResults = eclipseCaseData.results( RiaDefines::PorosityModelType::MATRIX_MODEL );
    CAF_ASSERT( matrixResults );

    RigEclipseResultAddress resAddr( RiaDefines::ResultCatType::STATIC_NATIVE, "ACTNUM" );
    if ( matrixResults->hasResultEntry( resAddr ) ) return;

    // One slot per active cell; all active by construction (see comment above), so fill with 1.0.
    size_t activeCellCount = activeCellInfo->reservoirActiveCellCount();
    matrixResults->addStaticScalarResult( RiaDefines::ResultCatType::STATIC_NATIVE, "ACTNUM", false, activeCellCount );

    auto modifiableData = matrixResults->modifiableCellScalarResultTimesteps( resAddr );
    CAF_ASSERT( modifiableData && !modifiableData->empty() );
    std::fill( ( *modifiableData )[0].begin(), ( *modifiableData )[0].end(), 1.0 );
}

//--------------------------------------------------------------------------------------------------
/// Compute the depth related geometry properties (DEPTH, DX, DY, DZ, TOPS, BOTTOM), matching the
/// behavior when a grid is imported from file.
//--------------------------------------------------------------------------------------------------
void RimCornerPointCase::computeDepthRelatedResults( RimCornerPointCase& cornerPointCase )
{
    if ( RiaPreferencesGrid::current()->autoComputeDepthRelatedProperties() )
    {
        cornerPointCase.eclipseCaseData()->computeDepthRelatedResults();
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimCornerPointCase::defineUiOrdering( QString uiConfigName, caf::PdmUiOrdering& uiOrdering )
{
    uiOrdering.add( &m_caseUserDescription );
    uiOrdering.add( &m_displayNameOption );
    uiOrdering.add( &m_caseId );
    uiOrdering.add( &m_caseFileName );

    auto group = uiOrdering.addNewGroup( "Case Options" );
    group->add( &m_activeFormationNames );
    group->add( &m_flipXAxis );
    group->add( &m_flipYAxis );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RimCornerPointCase::locationOnDisc() const
{
    if ( gridFileName().isEmpty() ) return QString();

    QFileInfo fi( gridFileName() );
    return fi.absolutePath();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimCornerPointCase::importAsciiInputProperties( const QStringList& fileNames )
{
    return true;
}