///////////////////////////////////////////////////////////////////////////////// // // 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 // for more details. // ///////////////////////////////////////////////////////////////////////////////// #include "RigMainGrid.h" #include "cvfAssert.h" #include "RimDefines.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(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(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& faults) { m_faults = faults; #pragma omp parallel for for (int i = 0; i < static_cast(m_faults.size()); i++) { m_faults[i]->computeFaultFacesFromCellRanges(this->mainGrid()); } } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RigMainGrid::calculateFaults() { //RigFault::initFaultsPrCellAccumulator(m_cells.size()); cvf::ref 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(fIdx)); } // Find the geometrical faults that is in addition RigFault * unNamedFault = new RigFault; int unNamedFaultIdx = static_cast(m_faults.size()); for (int gcIdx = 0 ; gcIdx < static_cast(m_cells.size()); ++gcIdx) { if ( m_cells[gcIdx].isInvalid()) { continue; } size_t neighborReservoirCellIdx; 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; } neighborReservoirCellIdx = hostGrid->reservoirCellIndex(neighborGridCellIdx); if (m_cells[neighborReservoirCellIdx].isInvalid()) { continue; } double tolerance = 1e-6; caf::SizeTArray4 faceIdxs; m_cells[gcIdx].faceIndices(face, &faceIdxs); caf::SizeTArray4 nbFaceIdxs; m_cells[neighborReservoirCellIdx].faceIndices(StructGridInterface::oppositeFace(face), &nbFaceIdxs); const std::vector& 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(neighborReservoirCellIdx, StructGridInterface::oppositeFace(face), unNamedFaultIdx); //m_cells[gcIdx].setCellFaceFault(face); //m_cells[neighborReservoirCellIdx].setCellFaceFault(StructGridInterface::oppositeFace(face)); // Add as fault face only if the grid index is less than the neighbors if (static_cast(gcIdx) < neighborReservoirCellIdx) { { RigFault::FaultFace ff(gcIdx, cvf::StructGridInterface::FaceType(faceIdx), neighborReservoirCellIdx); unNamedFault->faultFaces().push_back(ff); } } else { CVF_FAIL_MSG("Found fault with global neighbor index less than the native index. "); // Should never occur. because we flag the opposite face in the faultsPrCellAcc } } } } if (unNamedFault->faultFaces().size()) { unNamedFault->setName(RimDefines::undefinedGridFaultName()); m_faults.push_back(unNamedFault); } // Distribute nnc's to the faults const std::vector& 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; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const RigFault* RigMainGrid::findFaultFromCellIndexAndCellFace(size_t cellIndex, cvf::StructGridInterface::FaceType face) const { for (size_t i = 0; i < m_faults.size(); i++) { const RigFault* rigFault = m_faults.at(i); const std::vector& faultFaces = rigFault->faultFaces(); for (size_t fIdx = 0; fIdx < faultFaces.size(); fIdx++) { if (faultFaces[fIdx].m_nativeReservoirCellIndex == cellIndex) { if (face == faultFaces[fIdx].m_nativeFace ) { return rigFault; } } if (faultFaces[fIdx].m_oppositeReservoirCellIndex == cellIndex) { if (face == cvf::StructGridInterface::oppositeFace(faultFaces[fIdx].m_nativeFace)) { return rigFault; } } } } return NULL; }