ResInsight/ApplicationCode/ReservoirDataModel/RigMainGrid.cpp

/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2011-2012 Statoil ASA, Ceetron AS
// 
//  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 "RigMainGrid.h"

#include "cvfAssert.h"

RigMainGrid::RigMainGrid(void)
    : RigGridBase(this)
{
    m_displayModelOffset = cvf::Vec3d::ZERO;
    
    m_gridIndex = 0;
    m_gridId = 0;
    m_gridIdToIndexMapping.push_back(0); 

    m_flipXAxis = false;
    m_flipYAxis = false;
} 


RigMainGrid::~RigMainGrid(void)
{
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigMainGrid::addLocalGrid(RigLocalGrid* localGrid)
{
    CVF_ASSERT(localGrid && localGrid->gridId() != cvf::UNDEFINED_INT); // The grid ID must be set.
    CVF_ASSERT(localGrid->gridId() >= 0); // We cant handle negative ID's if they exist.

    m_localGrids.push_back(localGrid);
    localGrid->setGridIndex(m_localGrids.size()); // Maingrid itself has grid index 0

    
    if (m_gridIdToIndexMapping.size() <= static_cast<size_t>(localGrid->gridId()))
    {
        m_gridIdToIndexMapping.resize(localGrid->gridId() + 1, cvf::UNDEFINED_SIZE_T);
    }

    m_gridIdToIndexMapping[localGrid->gridId()] = localGrid->gridIndex();
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigMainGrid::initAllSubGridsParentGridPointer()
{
    initSubGridParentPointer();
    size_t i;
    for (i = 0; i < m_localGrids.size(); ++i)
    {
       m_localGrids[i]->initSubGridParentPointer(); 
    }
}


//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigMainGrid::initAllSubCellsMainGridCellIndex()
{
    initSubCellsMainGridCellIndex();
    size_t i;
    for (i = 0; i < m_localGrids.size(); ++i)
    {
        m_localGrids[i]->initSubCellsMainGridCellIndex(); 
    }
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigMainGrid::displayModelOffset() const
{
    return m_displayModelOffset;
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigMainGrid::setDisplayModelOffset(cvf::Vec3d offset)
{
    m_displayModelOffset = offset;
}

//--------------------------------------------------------------------------------------------------
/// Initialize pointers from grid to parent grid
/// Compute cell ranges for active and valid cells
/// Compute bounding box in world coordinates based on node coordinates
//--------------------------------------------------------------------------------------------------
void RigMainGrid::computeCachedData()
{
    initAllSubGridsParentGridPointer();
    initAllSubCellsMainGridCellIndex();
}

//--------------------------------------------------------------------------------------------------
/// Returns the grid with index \a localGridIndex. Main Grid itself has index 0. First LGR starts on 1
//--------------------------------------------------------------------------------------------------
RigGridBase* RigMainGrid::gridByIndex(size_t localGridIndex)
{
    if (localGridIndex == 0) return this;
    CVF_ASSERT(localGridIndex - 1 < m_localGrids.size()) ;
    return m_localGrids[localGridIndex-1].p();
}

//--------------------------------------------------------------------------------------------------
/// Returns the grid with index \a localGridIndex. Main Grid itself has index 0. First LGR starts on 1
//--------------------------------------------------------------------------------------------------
const RigGridBase* RigMainGrid::gridByIndex(size_t localGridIndex) const
{
    if (localGridIndex == 0) return this;
    CVF_ASSERT(localGridIndex - 1 < m_localGrids.size()) ;
    return m_localGrids[localGridIndex-1].p();
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigMainGrid::setFlipAxis(bool flipXAxis, bool flipYAxis)
{
    bool needFlipX = false;
    bool needFlipY = false;

    if (m_flipXAxis != flipXAxis)
    {
        needFlipX = true;
    }

    if (m_flipYAxis != flipYAxis)
    {
        needFlipY = true;
    }

    if (needFlipX || needFlipY)
    {
        for (size_t i = 0; i < m_nodes.size(); i++)
        {
            if (needFlipX)
            {
                m_nodes[i].x() *= -1.0;
            }

            if (needFlipY)
            {
                m_nodes[i].y() *= -1.0;
            }
        }

        m_flipXAxis = flipXAxis;
        m_flipYAxis = flipYAxis;
    }
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
RigGridBase* RigMainGrid::gridById(int localGridId)
{
    CVF_ASSERT (localGridId >= 0 && static_cast<size_t>(localGridId) < m_gridIdToIndexMapping.size());
    return this->gridByIndex(m_gridIdToIndexMapping[localGridId]);
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
RigNNCData* RigMainGrid::nncData()
{
    if (m_nncData.isNull()) 
    { 
        m_nncData = new RigNNCData;
    }  
    
    return m_nncData.p();
}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigMainGrid::setFaults(const cvf::Collection<RigFault>& faults)
{
    m_faults = faults;

#pragma omp parallel for 
    for (int i = 0; i < static_cast<int>(m_faults.size()); i++)
    {
        m_faults[i]->computeFaultFacesFromCellRanges(this->mainGrid());
    }
}


//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
void RigMainGrid::calculateFaults()
{
    //RigFault::initFaultsPrCellAccumulator(m_cells.size());
    cvf::ref<RigFaultsPrCellAccumulator> faultsPrCellAcc = new RigFaultsPrCellAccumulator(m_cells.size());

    // Spread fault idx'es on the cells from the faults
    for (size_t fIdx = 0 ; fIdx < m_faults.size(); ++fIdx)
    {
        m_faults[fIdx]->accumulateFaultsPrCell(faultsPrCellAcc.p(), static_cast<int>(fIdx));
    }

    // Find the geometrical faults that is in addition

    RigFault * unNamedFault = new RigFault;
    int unNamedFaultIdx = static_cast<int>(m_faults.size());

    for (size_t gcIdx = 0 ; gcIdx < m_cells.size(); ++gcIdx)
    {
        if ( m_cells[gcIdx].isInvalid())
        {
            continue;
        }

        size_t neighborGlobalCellIdx;
        size_t neighborGridCellIdx;
        size_t i, j, k;
        RigGridBase* hostGrid = NULL; 
        bool firstNO_FAULTFaceForCell = true;

        for (char faceIdx = 0; faceIdx < 6; ++faceIdx)
        {
            cvf::StructGridInterface::FaceType face = cvf::StructGridInterface::FaceType(faceIdx);

            if (faultsPrCellAcc->faultIdx(gcIdx, face) == RigFaultsPrCellAccumulator::NO_FAULT)
            {
                // Find neighbor cell
                if (firstNO_FAULTFaceForCell) // To avoid doing this for every face, and only when detecting a NO_FAULT
                {
                    hostGrid = m_cells[gcIdx].hostGrid();
                    hostGrid->ijkFromCellIndex(m_cells[gcIdx].cellIndex(), &i,&j, &k);
                    firstNO_FAULTFaceForCell = false;
                }

                if(!hostGrid->cellIJKNeighbor(i, j, k, face, &neighborGridCellIdx))
                {
                    continue;
                }

                neighborGlobalCellIdx = hostGrid->globalGridCellIndex(neighborGridCellIdx);
                if (m_cells[neighborGlobalCellIdx].isInvalid())
                {
                    continue;
                }

                double tolerance = 1e-6;
               

                caf::SizeTArray4 faceIdxs;
                m_cells[gcIdx].faceIndices(face, &faceIdxs);
                caf::SizeTArray4 nbFaceIdxs;
                m_cells[neighborGlobalCellIdx].faceIndices(StructGridInterface::oppositeFace(face), &nbFaceIdxs);

                const std::vector<cvf::Vec3d>& vxs = m_mainGrid->nodes();

                bool sharedFaceVertices = true;
                if (sharedFaceVertices && vxs[faceIdxs[0]].pointDistance(vxs[nbFaceIdxs[0]]) > tolerance ) sharedFaceVertices = false;
                if (sharedFaceVertices && vxs[faceIdxs[1]].pointDistance(vxs[nbFaceIdxs[3]]) > tolerance ) sharedFaceVertices = false;
                if (sharedFaceVertices && vxs[faceIdxs[2]].pointDistance(vxs[nbFaceIdxs[2]]) > tolerance ) sharedFaceVertices = false;
                if (sharedFaceVertices && vxs[faceIdxs[3]].pointDistance(vxs[nbFaceIdxs[1]]) > tolerance ) sharedFaceVertices = false;

                if (sharedFaceVertices)
                {
                    continue;
                }

                // To avoid doing this calculation for the opposite face

                faultsPrCellAcc->setFaultIdx(gcIdx, face, unNamedFaultIdx);
                faultsPrCellAcc->setFaultIdx(neighborGlobalCellIdx, StructGridInterface::oppositeFace(face), unNamedFaultIdx);

                m_cells[gcIdx].setCellFaceFault(face);
                m_cells[neighborGlobalCellIdx].setCellFaceFault(StructGridInterface::oppositeFace(face));

                // Add as fault face only if the grid index is less than the neighbors

                if (gcIdx < neighborGlobalCellIdx)
                {
                    RigFault::FaultFace ff(gcIdx, cvf::StructGridInterface::FaceType(faceIdx), neighborGlobalCellIdx);
                    unNamedFault->faultFaces().push_back(ff);
                }
                else
                {
                    CVF_FAIL_MSG("Found fault with global neighbour index less than the native index. "); // Should never occur. because we flag the opposite face in the faultsPrCellAcc
                }
            }
        }
    }

    if (unNamedFault->faultFaces().size())
    {
        unNamedFault->setName("Unnamed grid faults");
        m_faults.push_back(unNamedFault);
    }

    // Distribute nnc's to the faults

    const std::vector<RigConnection>& nncs = this->nncData()->connections();
    for (size_t nncIdx = 0; nncIdx  < nncs.size(); ++nncIdx)
    {
        // Find the fault for each side of the nnc
        const RigConnection& conn = nncs[nncIdx];
        int fIdx1 =  RigFaultsPrCellAccumulator::NO_FAULT;
        int fIdx2 =  RigFaultsPrCellAccumulator::NO_FAULT;

        if (conn.m_c1Face != StructGridInterface::NO_FACE)
        {
            fIdx1 = faultsPrCellAcc->faultIdx(conn.m_c1GlobIdx, conn.m_c1Face);
            fIdx2 = faultsPrCellAcc->faultIdx(conn.m_c2GlobIdx, StructGridInterface::oppositeFace(conn.m_c1Face));
        }

        if (fIdx1 < 0 && fIdx2 < 0)
        {
            cvf::String lgrString ("Same Grid");
            if (m_cells[conn.m_c1GlobIdx].hostGrid() != m_cells[conn.m_c2GlobIdx].hostGrid() )
            {
                lgrString = "Different Grid";
            }

           //cvf::Trace::show("NNC: No Fault for NNC C1: " + cvf::String((int)conn.m_c1GlobIdx) + " C2: " + cvf::String((int)conn.m_c2GlobIdx) + " Grid: " + lgrString);
        }

        if (fIdx1 >= 0)
        {
            // Add the connection to both, if they are different.
            m_faults[fIdx1]->connectionIndices().push_back(nncIdx);
        }

        if (fIdx2 != fIdx1)
        {
            if (fIdx2 >= 0)
            {
                m_faults[fIdx2]->connectionIndices().push_back(nncIdx);
            }
        }
    }
}

//--------------------------------------------------------------------------------------------------
/// The cell is normally inverted due to Depth becoming -Z at import, 
/// but if (only) one of the flipX/Y is done, the cell is back to nomal
//--------------------------------------------------------------------------------------------------
bool RigMainGrid::faceNormalsIsOutwards() const
{
    return  m_flipXAxis ^ m_flipYAxis;
}