From b1fa32c3187cc16cb52c951405e974b1bf027f64 Mon Sep 17 00:00:00 2001 From: jonjenssen <69144954+jonjenssen@users.noreply.github.com> Date: Mon, 4 Dec 2023 08:36:03 +0100 Subject: [PATCH] Fault reactivation: grid info improvements (#10899) * Fix some gridding issues, add element layer information and fix K layer information * Add geomech element corner node info in result info view --- .../RigFaultReactivationModelGenerator.cpp | 77 +++++++++++++------ .../RigFaultReactivationModelGenerator.h | 6 +- .../ReservoirDataModel/RigGriddedPart3d.cpp | 63 +++++++++++++-- .../ReservoirDataModel/RigGriddedPart3d.h | 3 + .../UserInterface/RiuFemResultTextBuilder.cpp | 11 +++ 5 files changed, 129 insertions(+), 31 deletions(-) diff --git a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp index 4701019d67..3f7142df90 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp @@ -365,6 +365,33 @@ const std::vector RigFaultReactivationModelGenerator::partition( double return parts; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigFaultReactivationModelGenerator::elementKLayers( const std::vector& cellIndexColumn ) +{ + std::vector kLayers; + + size_t i, j, k; + for ( auto idx : cellIndexColumn ) + { + m_grid->ijkFromCellIndexUnguarded( idx, &i, &j, &k ); + + if ( m_activeCellInfo->isActive( idx ) ) + { + kLayers.push_back( (int)k ); + } + else + { + kLayers.push_back( -1 * (int)k ); + } + } + + std::reverse( kLayers.begin(), kLayers.end() ); + + return kLayers; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -377,8 +404,6 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t std::vector cellColumnBack; std::vector cellColumnFront; size_t i, j, k; - std::vector kLayersFront; - std::vector kLayersBack; // build column of cells behind fault m_grid->ijkFromCellIndexUnguarded( startCellIndex, &i, &j, &k ); @@ -386,15 +411,12 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t for ( size_t kLayer = 0; kLayer < m_grid->cellCountK(); kLayer++ ) { + if ( !m_grid->isCellValid( i, j, kLayer ) ) continue; + auto cellIdx = m_grid->cellIndexFromIJKUnguarded( i, j, kLayer ); if ( cellIdx != startCellIndex ) cellColumnBackSearch.push_back( cellIdx ); cellColumnBack.push_back( cellIdx ); - - if ( m_activeCellInfo->isActive( cellIdx ) ) - kLayersBack.push_back( (int)kLayer ); - else - kLayersBack.push_back( -1 ); } // build cell column of cells in front of fault, opposite to the cell column behind the fault @@ -404,21 +426,17 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t m_grid->ijkFromCellIndexUnguarded( oppositeCellIdx, &i, &j, &k ); for ( size_t kLayer = 0; kLayer < m_grid->cellCountK(); kLayer++ ) { + if ( !m_grid->isCellValid( i, j, kLayer ) ) continue; + auto cellIdx = m_grid->cellIndexFromIJKUnguarded( i, j, kLayer ); cellColumnFront.push_back( cellIdx ); - - if ( m_activeCellInfo->isActive( cellIdx ) ) - kLayersFront.push_back( (int)kLayer ); - else - kLayersFront.push_back( -1 ); } auto zPositionsBack = elementLayers( startFace, cellColumnBack ); auto zPositionsFront = elementLayers( oppositeStartFace, cellColumnFront ); - - std::reverse( kLayersBack.begin(), kLayersBack.end() ); - std::reverse( kLayersFront.begin(), kLayersFront.end() ); + auto kLayersBack = elementKLayers( cellColumnBack ); + auto kLayersFront = elementKLayers( cellColumnFront ); // add extra fault buffer below the fault, starting at the deepest bottom-most cell on either side of the fault @@ -499,17 +517,17 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t /// //-------------------------------------------------------------------------------------------------- std::map RigFaultReactivationModelGenerator::elementLayers( cvf::StructGridInterface::FaceType face, - const std::vector& cellIndexColumn ) + std::vector& cellIndexColumn ) { cvf::Plane modelPlane; modelPlane.setFromPointAndNormal( m_startPosition, m_normal ); auto cornerIndexes = faceIJCornerIndexes( face ); - std::vector klayers; - std::map zPositions; + std::vector okCells; + for ( auto cellIdx : cellIndexColumn ) { RigCell cell = m_grid->cell( cellIdx ); @@ -518,8 +536,19 @@ std::map RigFaultReactivationModelGenerator::elementLayers( cvf::Vec3d intersect1 = lineIntersect( modelPlane, corners[cornerIndexes[0]], corners[cornerIndexes[1]] ); cvf::Vec3d intersect2 = lineIntersect( modelPlane, corners[cornerIndexes[2]], corners[cornerIndexes[3]] ); - zPositions[intersect1.z()] = intersect1; - zPositions[intersect2.z()] = intersect2; + if ( intersect1.z() != intersect2.z() ) + { + zPositions[intersect1.z()] = intersect1; + zPositions[intersect2.z()] = intersect2; + okCells.push_back( cellIdx ); + } + } + + // only keep cells that have a valid height at the plane intersection + cellIndexColumn.clear(); + for ( auto idx : okCells ) + { + cellIndexColumn.push_back( idx ); } return zPositions; @@ -548,10 +577,14 @@ void RigFaultReactivationModelGenerator::splitLargeLayers( std::map& layers, std::vector& kLayers, double maxHeight ); - std::map elementLayers( cvf::StructGridInterface::FaceType face, const std::vector& cellIndexColumn ); - void addFilter( QString name, std::vector cells ); + std::map elementLayers( cvf::StructGridInterface::FaceType face, std::vector& cellIndexColumn ); + std::vector elementKLayers( const std::vector& cellIndexColumn ); + + void addFilter( QString name, std::vector cells ); size_t oppositeStartCellIndex( const std::vector cellIndexColumn, cvf::StructGridInterface::FaceType face ); diff --git a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp index 24c2a4775d..33ee6b4ccc 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp @@ -59,6 +59,7 @@ void RigGriddedPart3d::reset() m_meshLines.clear(); m_elementSets.clear(); m_elementKLayer.clear(); + m_elementLayers.clear(); } //-------------------------------------------------------------------------------------------------- @@ -144,6 +145,9 @@ std::vector RigGriddedPart3d::generateGrowingLayers( double zFrom, doubl } } + if ( std::abs( zTo - layers.back() ) < maxSize ) layers.pop_back(); + layers.push_back( zTo ); + std::sort( layers.begin(), layers.end() ); return layers; @@ -164,6 +168,37 @@ std::vector RigGriddedPart3d::extractZValues( std::vector po return layers; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigGriddedPart3d::updateReservoirElementLayers( const std::vector& reservoirLayers, const std::vector& kLayers ) +{ + const int nLayers = (int)reservoirLayers.size(); + + if ( nLayers < 2 ) return; + + int prevLayer = kLayers[0]; + double start = reservoirLayers[0].z(); + + for ( int l = 1; l < nLayers; l++ ) + { + auto currentZone = ( prevLayer >= 0 ) ? ElementSets::Reservoir : ElementSets::IntraReservoir; + + if ( l == nLayers - 1 ) + { + m_elementLayers[currentZone].push_back( std::make_pair( start, reservoirLayers[l].z() ) ); + continue; + } + + if ( ( ( prevLayer < 0 ) && ( kLayers[l] >= 0 ) ) || ( ( prevLayer >= 0 ) && ( kLayers[l] < 0 ) ) ) + { + m_elementLayers[currentZone].push_back( std::make_pair( start, reservoirLayers[l].z() ) ); + start = reservoirLayers[l].z(); + } + prevLayer = kLayers[l]; + } +} + //-------------------------------------------------------------------------------------------------- /// Point index in input /// @@ -187,7 +222,6 @@ std::vector RigGriddedPart3d::extractZValues( std::vector po /// /// //-------------------------------------------------------------------------------------------------- - void RigGriddedPart3d::generateGeometry( const std::array& inputPoints, const std::vector& reservoirLayers, const std::vector& kLayers, @@ -204,13 +238,19 @@ void RigGriddedPart3d::generateGeometry( const std::array& input std::map> layersPerRegion; layersPerRegion[Regions::LowerUnderburden] = generateGrowingLayers( inputPoints[1].z(), inputPoints[0].z(), maxCellHeight, cellSizeFactor ); - layersPerRegion[Regions::LowerUnderburden].pop_back(); layersPerRegion[Regions::UpperUnderburden] = generateConstantLayers( inputPoints[1].z(), inputPoints[2].z(), maxCellHeight ); layersPerRegion[Regions::Reservoir] = extractZValues( reservoirLayers ); - layersPerRegion[Regions::Reservoir].pop_back(); - layersPerRegion[Regions::LowerOverburden] = generateConstantLayers( inputPoints[3].z(), inputPoints[4].z(), maxCellHeight ); + layersPerRegion[Regions::LowerOverburden] = generateConstantLayers( inputPoints[3].z(), inputPoints[4].z(), maxCellHeight ); layersPerRegion[Regions::UpperOverburden] = generateGrowingLayers( inputPoints[4].z(), inputPoints[5].z(), maxCellHeight, cellSizeFactor ); + layersPerRegion[Regions::LowerUnderburden].pop_back(); // to avoid overlap with bottom of next region + layersPerRegion[Regions::Reservoir].pop_back(); // to avoid overlap with bottom of next region + + m_elementLayers[ElementSets::OverBurden] = { std::make_pair( inputPoints[3].z(), inputPoints[5].z() ) }; + m_elementLayers[ElementSets::UnderBurden] = { std::make_pair( inputPoints[0].z(), inputPoints[2].z() ) }; + + updateReservoirElementLayers( reservoirLayers, kLayers ); + size_t nVertCells = 0; size_t nHorzCells = horizontalPartition.size() - 1; @@ -353,7 +393,7 @@ void RigGriddedPart3d::generateGeometry( const std::array& input const int nextLayerIdxOff = ( (int)nHorzCells + 1 ) * ( nThicknessCells + 1 ); const int nThicknessOff = nThicknessCells + 1; - for ( int v = 0; v < (int)nVertCells - 1; v++, layer++ ) + for ( int v = 0; v < (int)nVertCells - 1; v++ ) { if ( v >= nVertCellsLower ) currentSurfaceRegion = RimFaultReactivation::BorderSurface::FaultSurface; if ( v >= nVertCellsLower + nVertCellsFault ) currentSurfaceRegion = RimFaultReactivation::BorderSurface::UpperSurface; @@ -377,7 +417,7 @@ void RigGriddedPart3d::generateGeometry( const std::array& input m_elementIndices[elementIdx].push_back( t + i + nThicknessOff + 1 ); m_elementIndices[elementIdx].push_back( t + i + 1 ); - if ( layer == 0 ) + if ( v == 0 ) { m_boundaryElements[Boundary::Bottom].push_back( elementIdx ); } @@ -401,7 +441,7 @@ void RigGriddedPart3d::generateGeometry( const std::array& input else { m_elementSets[currentElementSet].push_back( elementIdx ); - m_elementKLayer[elementIdx] = -2; + m_elementKLayer[elementIdx] = -2000; } } i += nThicknessOff; @@ -555,6 +595,15 @@ const std::vector RigGriddedPart3d::elementKLayer() const return m_elementKLayer; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const std::vector> RigGriddedPart3d::layers( RigGriddedPart3d::ElementSets elementSet ) const +{ + if ( m_elementLayers.count( elementSet ) == 0 ) return {}; + return m_elementLayers.at( elementSet ); +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h index ba4d2c42c4..410bbccab7 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h +++ b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h @@ -70,6 +70,7 @@ public: const std::map>& elementSets() const; const std::vector elementKLayer() const; const std::vector elementCorners( size_t elementIndex ) const; + const std::vector> layers( ElementSets elementSet ) const; protected: static cvf::Vec3d stepVector( cvf::Vec3d start, cvf::Vec3d stop, int nSteps ); @@ -78,6 +79,7 @@ protected: static std::vector extractZValues( std::vector ); void generateVerticalMeshlines( const std::vector& cornerPoints, const std::vector& horzPartition ); + void updateReservoirElementLayers( const std::vector& reservoirLayers, const std::vector& kLayers ); private: enum class Regions @@ -105,4 +107,5 @@ private: std::map> m_boundaryElements; std::map> m_boundaryNodes; std::map> m_elementSets; + std::map>> m_elementLayers; }; diff --git a/ApplicationLibCode/UserInterface/RiuFemResultTextBuilder.cpp b/ApplicationLibCode/UserInterface/RiuFemResultTextBuilder.cpp index 3d87efc054..8ef22bbfbe 100644 --- a/ApplicationLibCode/UserInterface/RiuFemResultTextBuilder.cpp +++ b/ApplicationLibCode/UserInterface/RiuFemResultTextBuilder.cpp @@ -181,6 +181,17 @@ QString RiuFemResultTextBuilder::geometrySelectionText( QString itemSeparator ) } } } + std::array cornerCoords; + if ( femPart->fillElementCoordinates( m_cellIndex, cornerCoords ) ) + { + text += "\n\nElement corners:\n"; + + for ( auto p : cornerCoords ) + { + text += + QString( " [E: %1, N: %2, Depth: %3]\n" ).arg( p.x(), 5, 'f', 2 ).arg( p.y(), 5, 'f', 2 ).arg( -p.z(), 5, 'f', 2 ); + } + } } }