diff --git a/ApplicationCode/ReservoirDataModel/RigMainGrid.cpp b/ApplicationCode/ReservoirDataModel/RigMainGrid.cpp index 0cdb75af06..b843be7d20 100644 --- a/ApplicationCode/ReservoirDataModel/RigMainGrid.cpp +++ b/ApplicationCode/ReservoirDataModel/RigMainGrid.cpp @@ -448,113 +448,143 @@ void RigMainGrid::calculateFaults( const RigActiveCellInfo* activeCellInfo ) const std::vector& vxs = m_mainGrid->nodes(); + std::vector& unNamedFaultFaces = unNamedFault->faultFaces(); + std::vector& unNamedFaultFacesInactive = unNamedFaultWithInactive->faultFaces(); for ( int gcIdx = 0; gcIdx < static_cast( m_cells.size() ); ++gcIdx ) { - if ( m_cells[gcIdx].isInvalid() ) + addUnNamedFaultFaces( gcIdx, + activeCellInfo, + vxs, + unNamedFaultIdx, + unNamedFaultWithInactiveIdx, + unNamedFaultFaces, + unNamedFaultFacesInactive, + m_faultsPrCellAcc.p() ); + } + + if ( computeNncs ) + { + this->nncData()->computeCompleteSetOfNncs( this, activeCellInfo, includeInactiveCells ); + } + + distributeNNCsToFaults(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigMainGrid::addUnNamedFaultFaces( int gcIdx, + const RigActiveCellInfo* activeCellInfo, + const std::vector& vxs, + int unNamedFaultIdx, + int unNamedFaultWithInactiveIdx, + std::vector& unNamedFaultFaces, + std::vector& unNamedFaultFacesInactive, + RigFaultsPrCellAccumulator* faultsPrCellAcc ) const +{ + if ( m_cells[gcIdx].isInvalid() ) + { + return; + } + + size_t neighborReservoirCellIdx; + size_t neighborGridCellIdx; + size_t i = 0; + size_t j = 0; + size_t k = 0; + + const RigGridBase* hostGrid = nullptr; + bool firstNO_FAULTFaceForCell = true; + bool isCellActive = true; + + char upperLimitForFaceType = cvf::StructGridInterface::FaceType::POS_K; + + // Compare only I and J faces + for ( char faceIdx = 0; faceIdx < upperLimitForFaceType; ++faceIdx ) + { + cvf::StructGridInterface::FaceType face = cvf::StructGridInterface::FaceType( faceIdx ); + + // For faces that has no used defined Fault assigned: + + if ( m_faultsPrCellAcc->faultIdx( gcIdx, face ) == RigFaultsPrCellAccumulator::NO_FAULT ) { - continue; - } - - size_t neighborReservoirCellIdx; - size_t neighborGridCellIdx; - size_t i = 0; - size_t j = 0; - size_t k = 0; - - RigGridBase* hostGrid = nullptr; - bool firstNO_FAULTFaceForCell = true; - bool isCellActive = true; - - char upperLimitForFaceType = cvf::StructGridInterface::FaceType::POS_K; - - // Compare only I and J faces - for ( char faceIdx = 0; faceIdx < upperLimitForFaceType; ++faceIdx ) - { - cvf::StructGridInterface::FaceType face = cvf::StructGridInterface::FaceType( faceIdx ); - - // For faces that has no used defined Fault assigned: - - if ( m_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 { - // Find neighbor cell - if ( firstNO_FAULTFaceForCell ) // To avoid doing this for every face, and only when detecting a NO_FAULT + size_t gridLocalCellIndex; + hostGrid = this->gridAndGridLocalIdxFromGlobalCellIdx( gcIdx, &gridLocalCellIndex ); + + hostGrid->ijkFromCellIndex( gridLocalCellIndex, &i, &j, &k ); + isCellActive = activeCellInfo->isActive( gcIdx ); + + firstNO_FAULTFaceForCell = false; + } + + if ( !hostGrid->cellIJKNeighbor( i, j, k, face, &neighborGridCellIdx ) ) + { + continue; + } + + neighborReservoirCellIdx = hostGrid->reservoirCellIndex( neighborGridCellIdx ); + if ( m_cells[neighborReservoirCellIdx].isInvalid() ) + { + continue; + } + + bool isNeighborCellActive = activeCellInfo->isActive( neighborReservoirCellIdx ); + + double tolerance = 1e-6; + + std::array faceIdxs; + m_cells[gcIdx].faceIndices( face, &faceIdxs ); + std::array nbFaceIdxs; + m_cells[neighborReservoirCellIdx].faceIndices( StructGridInterface::oppositeFace( face ), &nbFaceIdxs ); + + 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 + int faultIdx = unNamedFaultIdx; + if ( !( isCellActive && isNeighborCellActive ) ) faultIdx = unNamedFaultWithInactiveIdx; + + faultsPrCellAcc->setFaultIdx( gcIdx, face, faultIdx ); + faultsPrCellAcc->setFaultIdx( neighborReservoirCellIdx, StructGridInterface::oppositeFace( face ), faultIdx ); + + // 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 ); + if ( isCellActive && isNeighborCellActive ) { - size_t gridLocalCellIndex; - hostGrid = this->gridAndGridLocalIdxFromGlobalCellIdx( gcIdx, &gridLocalCellIndex ); - - hostGrid->ijkFromCellIndex( gridLocalCellIndex, &i, &j, &k ); - isCellActive = activeCellInfo->isActive( gcIdx ); - - firstNO_FAULTFaceForCell = false; - } - - if ( !hostGrid->cellIJKNeighbor( i, j, k, face, &neighborGridCellIdx ) ) - { - continue; - } - - neighborReservoirCellIdx = hostGrid->reservoirCellIndex( neighborGridCellIdx ); - if ( m_cells[neighborReservoirCellIdx].isInvalid() ) - { - continue; - } - - bool isNeighborCellActive = activeCellInfo->isActive( neighborReservoirCellIdx ); - - double tolerance = 1e-6; - - std::array faceIdxs; - m_cells[gcIdx].faceIndices( face, &faceIdxs ); - std::array nbFaceIdxs; - m_cells[neighborReservoirCellIdx].faceIndices( StructGridInterface::oppositeFace( face ), &nbFaceIdxs ); - - 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 - int faultIdx = unNamedFaultIdx; - if ( !( isCellActive && isNeighborCellActive ) ) faultIdx = unNamedFaultWithInactiveIdx; - - m_faultsPrCellAcc->setFaultIdx( gcIdx, face, faultIdx ); - m_faultsPrCellAcc->setFaultIdx( neighborReservoirCellIdx, - StructGridInterface::oppositeFace( face ), - faultIdx ); - - // 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 ); - if ( isCellActive && isNeighborCellActive ) - { - unNamedFault->faultFaces().push_back( ff ); - } - else - { - unNamedFaultWithInactive->faultFaces().push_back( ff ); - } + unNamedFaultFaces.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 + unNamedFaultFacesInactive.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 + } } } } diff --git a/ApplicationCode/ReservoirDataModel/RigMainGrid.h b/ApplicationCode/ReservoirDataModel/RigMainGrid.h index 8ba09d3340..b0d1ac25b9 100644 --- a/ApplicationCode/ReservoirDataModel/RigMainGrid.h +++ b/ApplicationCode/ReservoirDataModel/RigMainGrid.h @@ -73,6 +73,15 @@ public: cvf::Collection& faults(); void calculateFaults( const RigActiveCellInfo* activeCellInfo ); + void addUnNamedFaultFaces( int gcIdx, + const RigActiveCellInfo* activeCellInfo, + const std::vector& vxs, + int unNamedFaultIdx, + int unNamedFaultWithInactiveIdx, + std::vector& unNamedFaultFaces, + std::vector& unNamedFaultFacesInactive, + RigFaultsPrCellAccumulator* faultsPrCellAcc ) const; + void distributeNNCsToFaults(); const RigFault* findFaultFromCellIndexAndCellFace( size_t reservoirCellIndex,