/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2018 Statoil 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 "RimFractureContainmentTools.h"

#include "RigEclipseCaseData.h"
#include "RigFault.h"
#include "RigHexIntersectionTools.h"
#include "RigMainGrid.h"
#include "RigReservoirGridTools.h"

#include "RimEclipseCase.h"
#include "RimEclipseView.h"
#include "RimFracture.h"
#include "RimFractureContainment.h"
#include "RimFractureTemplate.h"

#include "cvfStructGrid.h"

#include <array>

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t findNeighborReservoirCellIndex(const RigMainGrid*                 mainGrid,
                                      cvf::StructGridInterface::FaceType face,
                                      size_t                             globalReservoirCellIndex)
{
    size_t neighborGlobalReservoirCellIndex = cvf::UNDEFINED_SIZE_T;

    if (mainGrid)
    {
        size_t             gridLocalCellIndex = cvf::UNDEFINED_SIZE_T;
        const RigGridBase* hostGrid =
            mainGrid->gridAndGridLocalIdxFromGlobalCellIdx(globalReservoirCellIndex, &gridLocalCellIndex);

        if (hostGrid && gridLocalCellIndex != cvf::UNDEFINED_SIZE_T)
        {
            size_t i, j, k;
            hostGrid->ijkFromCellIndex(gridLocalCellIndex, &i, &j, &k);

            size_t neighborGridLocalCellIndex;

            bool foundCell = hostGrid->cellIJKNeighbor(i, j, k, face, &neighborGridLocalCellIndex);
            if (foundCell)
            {
                neighborGlobalReservoirCellIndex = hostGrid->reservoirCellIndex(neighborGridLocalCellIndex);
            }
        }
    }

    return neighborGlobalReservoirCellIndex;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimFractureContainmentTools::appendNeighborCellForFace(const std::set<size_t>&            allFracturedCells,
                                                            const RigMainGrid*                 mainGrid,
                                                            size_t                             currentCell,
                                                            cvf::StructGridInterface::FaceType face,
                                                            std::set<size_t>&                  connectedCells,
                                                            double                             minimumFaultThrow)
{
    size_t candidate = findNeighborReservoirCellIndex(mainGrid, face, currentCell);
    if (candidate != cvf::UNDEFINED_SIZE_T)
    {
        appendNeighborCells(allFracturedCells, mainGrid, candidate, connectedCells, minimumFaultThrow);
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double computeAverageZForTwoDeepestZ(const RigMainGrid*                 mainGrid,
                                     size_t                             globalReservoirCellIndex,
                                     cvf::StructGridInterface::FaceType face)
{
    const RigCell& currentCell = mainGrid->globalCellArray()[globalReservoirCellIndex];

    cvf::Vec3d hexCorners[8];
    mainGrid->cellCornerVertices(globalReservoirCellIndex, hexCorners);

    double avgZ = 0.0;

    cvf::ubyte faceVertexIndices[4];
    cvf::StructGridInterface::cellFaceVertexIndices(face, faceVertexIndices);

    for (const auto& faceIdx : faceVertexIndices)
    {
        // Face indices 0-3 are defined to have deepest Z
        // See void StructGridInterface::cellFaceVertexIndices(FaceType face, cvf::ubyte vertexIndices[4])

        if (faceIdx < 4)
        {
            avgZ += hexCorners[faceIdx].z();
        }
    }

    avgZ /= 2.0;

    return avgZ;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimFractureContainmentTools::checkFaultAndAppendNeighborCell(const std::set<size_t>&            allFracturedCells,
                                                                  const RigMainGrid*                 mainGrid,
                                                                  size_t                             globalReservoirCellIndex,
                                                                  cvf::StructGridInterface::FaceType face,
                                                                  std::set<size_t>&                  connectedCells,
                                                                  double                             minimumFaultThrow)
{
    const RigFault* fault = mainGrid->findFaultFromCellIndexAndCellFace(globalReservoirCellIndex, face);
    if (fault)
    {
        {
            // See RigMainGrid::calculateFaults() for reference

            // This function is intended to support fractures in LGR-grids
            // Currently, only faults in main grid is supported when reading fault specifications from input text files
            // Eclipse 300 supports faults in LGR
            // https://github.com/OPM/ResInsight/issues/3019

            size_t neighborGlobalReservoirCellIndex = findNeighborReservoirCellIndex(mainGrid, face, globalReservoirCellIndex);
            if (neighborGlobalReservoirCellIndex == cvf::UNDEFINED_SIZE_T)
            {
                // This is probably an assert condition, but we return directly to ensure we are robust
                return;
            }

            double currentCellAvgZ  = computeAverageZForTwoDeepestZ(mainGrid, globalReservoirCellIndex, face);
            double neighborCellAvgZ = computeAverageZForTwoDeepestZ(
                mainGrid, neighborGlobalReservoirCellIndex, cvf::StructGridInterface::oppositeFace(face));

            double faultThrow = fabs(currentCellAvgZ - neighborCellAvgZ);
            if (faultThrow > minimumFaultThrow)
            {
                return;
            }
        }
    }

    appendNeighborCellForFace(allFracturedCells, mainGrid, globalReservoirCellIndex, face, connectedCells, minimumFaultThrow);
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimFractureContainmentTools::appendNeighborCells(const std::set<size_t>& allFracturedCells,
                                                      const RigMainGrid*      mainGrid,
                                                      size_t                  currentCell,
                                                      std::set<size_t>&       connectedCells,
                                                      double                  minimumFaultThrow)
{
    if (std::find(connectedCells.begin(), connectedCells.end(), currentCell) != connectedCells.end())
    {
        // currentCell is already handled
        return;
    }

    if (std::find(allFracturedCells.begin(), allFracturedCells.end(), currentCell) == allFracturedCells.end())
    {
        // currentCell is not found among the set of fracture cells
        return;
    }

    connectedCells.insert(currentCell);

    // Check faults in IJ directions
    checkFaultAndAppendNeighborCell(
        allFracturedCells, mainGrid, currentCell, cvf::StructGridInterface::NEG_I, connectedCells, minimumFaultThrow);
    checkFaultAndAppendNeighborCell(
        allFracturedCells, mainGrid, currentCell, cvf::StructGridInterface::POS_I, connectedCells, minimumFaultThrow);
    checkFaultAndAppendNeighborCell(
        allFracturedCells, mainGrid, currentCell, cvf::StructGridInterface::NEG_J, connectedCells, minimumFaultThrow);
    checkFaultAndAppendNeighborCell(
        allFracturedCells, mainGrid, currentCell, cvf::StructGridInterface::POS_J, connectedCells, minimumFaultThrow);

    // Append cells without fault check in K direction
    appendNeighborCellForFace(
        allFracturedCells, mainGrid, currentCell, cvf::StructGridInterface::NEG_K, connectedCells, minimumFaultThrow);
    appendNeighborCellForFace(
        allFracturedCells, mainGrid, currentCell, cvf::StructGridInterface::POS_K, connectedCells, minimumFaultThrow);
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::set<size_t> RimFractureContainmentTools::fracturedCellsTruncatedByFaults(const RimEclipseCase* eclipseCase,
                                                                              const RimFracture*    fracture)
{
    std::set<size_t> fracturedCellsContainedByFaults;

    if (eclipseCase && fracture)
    {
        auto eclipseCaseData = eclipseCase->eclipseCaseData();
        if (eclipseCaseData)
        {
            auto mainGrid = eclipseCaseData->mainGrid();

            if (mainGrid)
            {
                std::set<size_t> cellsIntersectingFracturePlane = getCellsIntersectingFracturePlane(mainGrid, fracture);

                // Negative faultThrow disables test on faultThrow
                double maximumFaultThrow = -1.0;
                if (fracture->fractureTemplate() && fracture->fractureTemplate()->fractureContainment())
                {
                    maximumFaultThrow = fracture->fractureTemplate()->fractureContainment()->minimumFaultThrow();
                }

                if (maximumFaultThrow > -1.0)
                {
                    size_t anchorCellGlobalIndex = mainGrid->findReservoirCellIndexFromPoint(fracture->anchorPosition());
                    appendNeighborCells(cellsIntersectingFracturePlane,
                                        mainGrid,
                                        anchorCellGlobalIndex,
                                        fracturedCellsContainedByFaults,
                                        maximumFaultThrow);
                }
                else
                {
                    fracturedCellsContainedByFaults = cellsIntersectingFracturePlane;
                }
            }

            /*
            NB : Please do not delete this code, used to create input to cell based range filter to see the computed fracture
            cells

                        qDebug() << "FracturedCells - Truncated";
                        qDebug() << RigReservoirGridTools::globalCellIndicesToOneBasedIJKText(
                            fracturedCellsContainedByFaults.begin(), fracturedCellsContainedByFaults.end(), mainGrid);
            */
        }
    }

    return fracturedCellsContainedByFaults;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::set<size_t> RimFractureContainmentTools::getCellsIntersectingFracturePlane(const RigMainGrid* mainGrid,
                                                                                const RimFracture* fracture)
{
    std::set<size_t> eclipseCellIndices;

    {
        const auto indicesToPotentiallyFracturedCells = fracture->getPotentiallyFracturedCells(mainGrid);

        for (const auto& globalCellIndex : indicesToPotentiallyFracturedCells)
        {
            std::array<cvf::Vec3d, 8> hexCorners;
            mainGrid->cellCornerVertices(globalCellIndex, hexCorners.data());
            std::vector<std::vector<cvf::Vec3d>> planeCellPolygons;

            bool isPlanIntersected =
                RigHexIntersectionTools::planeHexIntersectionPolygons(hexCorners, fracture->transformMatrix(), planeCellPolygons);
            if (isPlanIntersected || !planeCellPolygons.empty())
            {
                eclipseCellIndices.insert(globalCellIndex);
            }
        }
    }

    return eclipseCellIndices;
}