From 336dc575db64a97dd429bf641e6d7f9bfbb8a890 Mon Sep 17 00:00:00 2001 From: jonjenssen <69144954+jonjenssen@users.noreply.github.com> Date: Mon, 12 Feb 2024 14:52:15 +0100 Subject: [PATCH] Fault Reactivation: Update element sets (#11191) * Update element sets based on active cells * Bonus: Fix viewer crash --- .../RigFaultReactivationModelGenerator.cpp | 69 +----- .../RigFaultReactivationModelGenerator.h | 6 +- .../ReservoirDataModel/RigGriddedPart3d.cpp | 223 ++++++++---------- .../ReservoirDataModel/RigGriddedPart3d.h | 23 +- .../UserInterface/RiuViewer.cpp | 21 +- 5 files changed, 138 insertions(+), 204 deletions(-) diff --git a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp index 7881484a63..a2b6d31524 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp @@ -361,33 +361,6 @@ 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; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -543,11 +516,6 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t sta std::map layersFront; std::vector cellColumnFront = buildCellColumn( oppositeStartCellIdx, oppositeStartFace, layersFront ); - // identify k layers used to classify reservoir/inter-reservoir - - 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 double front_bottom = layersFront.begin()->first; @@ -588,11 +556,11 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t sta // make sure layers aren't too small or too thick - mergeTinyLayers( layersFront, kLayersFront, m_minCellHeight ); - mergeTinyLayers( layersBack, kLayersBack, m_minCellHeight ); + mergeTinyLayers( layersFront, m_minCellHeight ); + mergeTinyLayers( layersBack, m_minCellHeight ); - splitLargeLayers( layersFront, kLayersFront, m_maxCellHeight ); - splitLargeLayers( layersBack, kLayersBack, m_maxCellHeight ); + splitLargeLayers( layersFront, m_maxCellHeight ); + splitLargeLayers( layersBack, m_maxCellHeight ); std::vector frontReservoirLayers; for ( auto& kvp : layersFront ) @@ -608,7 +576,6 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t sta frontPart->generateGeometry( m_frontPoints, frontReservoirLayers, - kLayersFront, m_maxCellHeight, m_cellSizeHeightFactor, m_horizontalPartition, @@ -618,7 +585,6 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t sta m_faultZoneCells ); backPart->generateGeometry( m_backPoints, backReservoirLayers, - kLayersBack, m_maxCellHeight, m_cellSizeHeightFactor, m_horizontalPartition, @@ -674,9 +640,8 @@ cvf::Vec3d RigFaultReactivationModelGenerator::extrapolatePoint( cvf::Vec3d star //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RigFaultReactivationModelGenerator::mergeTinyLayers( std::map& layers, std::vector& kLayers, double minHeight ) +void RigFaultReactivationModelGenerator::mergeTinyLayers( std::map& layers, double minHeight ) { - std::vector newKLayers; std::vector newLayers; const int nLayers = (int)layers.size(); @@ -692,7 +657,6 @@ void RigFaultReactivationModelGenerator::mergeTinyLayers( std::map& layers, std::vector& kLayers, double maxHeight ) +void RigFaultReactivationModelGenerator::splitLargeLayers( std::map& layers, double maxHeight ) { std::vector additionalPoints; - std::vector newKLayers; - - const int nLayers = (int)layers.size(); - const int nKLayers = (int)kLayers.size(); + const int nLayers = (int)layers.size(); std::vector keys; std::vector vals; @@ -764,23 +717,15 @@ void RigFaultReactivationModelGenerator::splitLargeLayers( std::map sideFacesIJ( FaceType face ); static cvf::Vec3d extrapolatePoint( cvf::Vec3d startPoint, cvf::Vec3d endPoint, double stopDepth ); - static void splitLargeLayers( std::map& layers, std::vector& kLayers, double maxHeight ); - static void mergeTinyLayers( std::map& layers, std::vector& kLayers, double minHeight ); - - std::vector elementKLayers( const std::vector& cellIndexColumn ); + static void splitLargeLayers( std::map& layers, double maxHeight ); + static void mergeTinyLayers( std::map& layers, double minHeight ); std::vector buildCellColumn( size_t startCell, FaceType startFace, std::map& layers ); diff --git a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp index 02d0581d19..07bc7d89b0 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp @@ -37,6 +37,8 @@ RigGriddedPart3d::RigGriddedPart3d() : m_useLocalCoordinates( false ) , m_topHeight( 0.0 ) , m_faultSafetyDistance( 1.0 ) + , m_nVertElements( 0 ) + , m_nHorzElements( 0 ) { } @@ -61,8 +63,9 @@ void RigGriddedPart3d::reset() m_elementIndices.clear(); m_meshLines.clear(); m_elementSets.clear(); - m_elementKLayer.clear(); - m_elementLayers.clear(); + m_nVertElements = 0; + m_nHorzElements = 0; + m_topHeight = 0.0; } //-------------------------------------------------------------------------------------------------- @@ -174,37 +177,6 @@ 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 /// @@ -226,11 +198,9 @@ void RigGriddedPart3d::updateReservoirElementLayers( const std::vector& inputPoints, const std::vector& reservoirLayers, - const std::vector& kLayers, const double maxCellHeight, double cellSizeFactor, const std::vector& horizontalPartition, @@ -254,17 +224,12 @@ void RigGriddedPart3d::generateGeometry( const std::array& input 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() ) }; - m_boundaryNodes[Boundary::Bottom] = {}; m_boundaryNodes[Boundary::FarSide] = {}; m_boundaryNodes[Boundary::Fault] = {}; - updateReservoirElementLayers( reservoirLayers, kLayers ); - - size_t nVertCells = 0; - size_t nHorzCells = horizontalPartition.size() - 1; + size_t nVertCells = 0; + const size_t nHorzCells = horizontalPartition.size() - 1; for ( auto region : allRegions() ) { @@ -279,6 +244,9 @@ void RigGriddedPart3d::generateGeometry( const std::array& input m_nodes.reserve( reserveSize ); m_dataNodes.reserve( reserveSize ); + m_nHorzElements = (int)nHorzCells; + m_nVertElements = (int)nVertCells - 1; + unsigned int nodeIndex = 0; unsigned int layer = 0; @@ -389,7 +357,6 @@ void RigGriddedPart3d::generateGeometry( const std::array& input // ** generate elements of type hex8 m_elementIndices.resize( (size_t)( ( nVertCells - 1 ) * nHorzCells * nThicknessCells ) ); - m_elementKLayer.resize( (size_t)( ( nVertCells - 1 ) * nHorzCells * nThicknessCells ) ); m_borderSurfaceElements[RimFaultReactivation::BorderSurface::Seabed] = {}; m_borderSurfaceElements[RimFaultReactivation::BorderSurface::UpperSurface] = {}; @@ -409,18 +376,12 @@ void RigGriddedPart3d::generateGeometry( const std::array& input int layerIndexOffset = 0; int elementIdx = 0; layer = 0; - int kLayer = 0; const int nVertCellsLower = (int)layersPerRegion[Regions::LowerUnderburden].size(); const int nVertCellsFault = (int)( layersPerRegion[Regions::UpperUnderburden].size() + layersPerRegion[Regions::Reservoir].size() + layersPerRegion[Regions::LowerOverburden].size() ); - const int nVertCellsUnderburden = - (int)( layersPerRegion[Regions::LowerUnderburden].size() + layersPerRegion[Regions::UpperUnderburden].size() ); - const int nVertCellsReservoir = nVertCellsUnderburden + (int)( layersPerRegion[Regions::Reservoir].size() ); - RimFaultReactivation::BorderSurface currentSurfaceRegion = RimFaultReactivation::BorderSurface::LowerSurface; - RimFaultReactivation::ElementSets currentElementSet = RimFaultReactivation::ElementSets::UnderBurden; const int nextLayerIdxOff = ( (int)nHorzCells + 1 ) * ( nThicknessCells + 1 ); const int nThicknessOff = nThicknessCells + 1; @@ -433,9 +394,6 @@ void RigGriddedPart3d::generateGeometry( const std::array& input if ( v >= nVertCellsLower ) currentSurfaceRegion = RimFaultReactivation::BorderSurface::FaultSurface; if ( v >= nVertCellsLower + nVertCellsFault ) currentSurfaceRegion = RimFaultReactivation::BorderSurface::UpperSurface; - if ( v >= nVertCellsUnderburden ) currentElementSet = RimFaultReactivation::ElementSets::Reservoir; - if ( v >= nVertCellsReservoir ) currentElementSet = RimFaultReactivation::ElementSets::OverBurden; - int i = layerIndexOffset; for ( int h = 0; h < (int)nHorzCells; h++ ) @@ -468,27 +426,6 @@ void RigGriddedPart3d::generateGeometry( const std::array& input bool inFaultZone = ( currentSurfaceRegion == RimFaultReactivation::BorderSurface::FaultSurface ) && ( h > nFaultZoneStart ); if ( inFaultZone ) m_elementSets[RimFaultReactivation::ElementSets::FaultZone].push_back( elementIdx ); - - if ( currentElementSet == RimFaultReactivation::ElementSets::Reservoir ) - { - m_elementKLayer[elementIdx] = kLayers[kLayer]; - if ( !inFaultZone ) - { - if ( kLayers[kLayer] < 0 ) - { - m_elementSets[RimFaultReactivation::ElementSets::IntraReservoir].push_back( elementIdx ); - } - else - { - m_elementSets[currentElementSet].push_back( elementIdx ); - } - } - } - else - { - if ( !inFaultZone ) m_elementSets[currentElementSet].push_back( elementIdx ); - m_elementKLayer[elementIdx] = -2000; - } } i += nThicknessOff; } @@ -497,11 +434,6 @@ void RigGriddedPart3d::generateGeometry( const std::array& input m_borderSurfaceElements[currentSurfaceRegion].push_back( elementIdx - 2 ); m_borderSurfaceElements[currentSurfaceRegion].push_back( elementIdx - 1 ); - if ( currentElementSet == RimFaultReactivation::ElementSets::Reservoir ) - { - kLayer++; - } - layerIndexOffset += nextLayerIdxOff; } @@ -628,7 +560,7 @@ const std::vector>& RigGriddedPart3d::elementIndices() //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::vector RigGriddedPart3d::elementCorners( size_t elementIndex ) const +const std::vector RigGriddedPart3d::elementCorners( size_t elementIndex ) const { return extractCornersForElement( m_elementIndices, m_nodes, elementIndex ); } @@ -636,11 +568,19 @@ std::vector RigGriddedPart3d::elementCorners( size_t elementIndex ) //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::vector RigGriddedPart3d::elementDataCorners( size_t elementIndex ) const +const std::vector RigGriddedPart3d::elementDataCorners( size_t elementIndex ) const { return extractCornersForElement( m_elementIndices, m_dataNodes, elementIndex ); } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const std::pair RigGriddedPart3d::elementCountHorzVert() const +{ + return { m_nHorzElements, m_nVertElements }; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -677,23 +617,6 @@ const std::vector>& RigGriddedPart3d::meshLines() const return m_meshLines; } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -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 ); -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -739,52 +662,108 @@ void RigGriddedPart3d::generateLocalNodes( const cvf::Mat4d transform ) } //-------------------------------------------------------------------------------------------------- -/// Make sure any active cells outside the flat reservoir zone is added to the reservoir element set +/// //-------------------------------------------------------------------------------------------------- void RigGriddedPart3d::postProcessElementSets( const RigMainGrid* mainGrid, const RigActiveCellInfo* cellInfo ) { - std::map> newElementSets; + std::set usedElements; + + // fault zone elements are already assigned + for ( auto elIdx : m_elementSets[ElementSets::FaultZone] ) + { + usedElements.insert( elIdx ); + } - for ( auto elSet : { ElementSets::OverBurden, ElementSets::UnderBurden, ElementSets::IntraReservoir } ) + // look for overburden, starting at top going down + updateElementSet( ElementSets::OverBurden, usedElements, mainGrid, cellInfo, m_nVertElements - 1, -1, -1 ); + + // look for underburden, starting at bottom going up + updateElementSet( ElementSets::UnderBurden, usedElements, mainGrid, cellInfo, 0, m_nVertElements, 1 ); + + // remaining elements are in the reservoir + m_elementSets[ElementSets::IntraReservoir] = {}; + m_elementSets[ElementSets::Reservoir] = {}; + + for ( unsigned int element = 0; element < m_elementIndices.size(); element++ ) { - newElementSets[elSet] = {}; + if ( usedElements.contains( element ) ) continue; + + auto corners = elementCorners( element ); + bool bActive = false; + + size_t cellIdx = 0; + for ( const auto& p : corners ) + { + cellIdx = mainGrid->findReservoirCellIndexFromPoint( p ); + + bActive = ( cellIdx != cvf::UNDEFINED_SIZE_T ) && ( cellInfo->isActive( cellIdx ) ); + if ( bActive ) break; + } + + if ( bActive ) + { + m_elementSets[ElementSets::Reservoir].push_back( element ); + } + else + { + m_elementSets[ElementSets::IntraReservoir].push_back( element ); + } + } +} - for ( auto element : m_elementSets[elSet] ) +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigGriddedPart3d::updateElementSet( ElementSets elSet, + std::set& usedElements, + const RigMainGrid* mainGrid, + const RigActiveCellInfo* cellInfo, + int rowStart, + int rowEnd, + int rowInc ) +{ + for ( int col = 0; col < m_nHorzElements; col++ ) + { + for ( int row = rowStart; row != rowEnd; row += rowInc ) { - auto corners = elementCorners( element ); - int nFound = 0; + const unsigned int elIdx = (unsigned int)( 2 * ( ( row * m_nHorzElements ) + col ) ); - size_t cellIdx = 0; - for ( const auto& p : corners ) + bool bStop = false; + + for ( unsigned int t = 0; t < 2; t++ ) { - cellIdx = mainGrid->findReservoirCellIndexFromPoint( p ); + if ( usedElements.contains( elIdx + t ) ) + { + bStop = true; + break; + } - if ( cellIdx != cvf::UNDEFINED_SIZE_T ) + auto corners = elementCorners( elIdx + t ); + + size_t cellIdx = 0; + for ( const auto& p : corners ) { - if ( cellInfo->isActive( cellIdx ) ) + cellIdx = mainGrid->findReservoirCellIndexFromPoint( p ); + + if ( ( cellIdx != cvf::UNDEFINED_SIZE_T ) && ( cellInfo->isActive( cellIdx ) ) ) { - nFound++; + bStop = true; + break; } } } - if ( nFound > 0 ) + if ( bStop ) { - m_elementSets[ElementSets::Reservoir].push_back( element ); + break; } else { - newElementSets[elSet].push_back( element ); + m_elementSets[elSet].push_back( elIdx ); + m_elementSets[elSet].push_back( elIdx + 1 ); + usedElements.insert( elIdx ); + usedElements.insert( elIdx + 1 ); } } } - - for ( auto elSet : { ElementSets::OverBurden, ElementSets::UnderBurden, ElementSets::IntraReservoir } ) - { - m_elementSets[elSet].clear(); - for ( auto element : newElementSets[elSet] ) - { - m_elementSets[elSet].push_back( element ); - } - } } diff --git a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h index 3ce27e5544..49f44269c8 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h +++ b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h @@ -26,6 +26,7 @@ #include #include +#include #include class RigMainGrid; @@ -48,7 +49,6 @@ class RigGriddedPart3d : public cvf::Object void generateGeometry( const std::array& inputPoints, const std::vector& reservoirLayers, - const std::vector& kLayers, double maxCellHeight, double cellSizeFactor, const std::vector& horizontalPartition, @@ -79,10 +79,9 @@ class RigGriddedPart3d : public cvf::Object const std::map>& boundaryElements() const; const std::map>& boundaryNodes() const; const std::map>& elementSets() const; - const std::vector elementKLayer() const; - std::vector elementCorners( size_t elementIndex ) const; - std::vector elementDataCorners( size_t elementIndex ) const; - const std::vector> layers( ElementSets elementSet ) const; + const std::vector elementCorners( size_t elementIndex ) const; + const std::vector elementDataCorners( size_t elementIndex ) const; + const std::pair elementCountHorzVert() const; protected: static cvf::Vec3d stepVector( cvf::Vec3d start, cvf::Vec3d stop, int nSteps ); @@ -91,7 +90,14 @@ class RigGriddedPart3d : public cvf::Object 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 ); + + void updateElementSet( ElementSets elSet, + std::set& usedElements, + const RigMainGrid* mainGrid, + const RigActiveCellInfo* cellInfo, + int rowStart, + int rowEnd, + int rowInc ); private: enum class Regions @@ -114,15 +120,16 @@ class RigGriddedPart3d : public cvf::Object double m_topHeight; double m_faultSafetyDistance; + int m_nVertElements; + int m_nHorzElements; + std::vector m_nodes; std::vector m_dataNodes; std::vector m_localNodes; std::vector> m_elementIndices; - std::vector m_elementKLayer; std::map> m_borderSurfaceElements; std::vector> m_meshLines; std::map> m_boundaryElements; std::map> m_boundaryNodes; std::map> m_elementSets; - std::map>> m_elementLayers; }; diff --git a/ApplicationLibCode/UserInterface/RiuViewer.cpp b/ApplicationLibCode/UserInterface/RiuViewer.cpp index ea1b8d0ecf..c929010946 100644 --- a/ApplicationLibCode/UserInterface/RiuViewer.cpp +++ b/ApplicationLibCode/UserInterface/RiuViewer.cpp @@ -482,17 +482,20 @@ void RiuViewer::paintOverlayItems( QPainter* painter ) { Rim3dView* view = dynamic_cast( m_rimView.p() ); - QString stepName = view->timeStepName( view->currentTimeStep() ); + if ( view ) + { + QString stepName = view->timeStepName( view->currentTimeStep() ); - m_animationProgress->setFormat( "Time Step: %v/%m " + stepName ); - m_animationProgress->setMinimum( 0 ); - m_animationProgress->setMaximum( static_cast( view->timeStepCount() ) - 1 ); - m_animationProgress->setValue( view->currentTimeStep() ); + m_animationProgress->setFormat( "Time Step: %v/%m " + stepName ); + m_animationProgress->setMinimum( 0 ); + m_animationProgress->setMaximum( static_cast( view->timeStepCount() ) - 1 ); + m_animationProgress->setValue( view->currentTimeStep() ); - m_animationProgress->resize( columnWidth, m_animationProgress->sizeHint().height() ); - m_animationProgress->render( painter, QPoint( columnPos, yPos ) ); + m_animationProgress->resize( columnWidth, m_animationProgress->sizeHint().height() ); + m_animationProgress->render( painter, QPoint( columnPos, yPos ) ); - yPos += m_animationProgress->height() + margin; + yPos += m_animationProgress->height() + margin; + } } if ( m_showInfoText && !isComparisonViewActive() ) @@ -1037,6 +1040,8 @@ RimViewWindow* RiuViewer::ownerViewWindow() const //-------------------------------------------------------------------------------------------------- void RiuViewer::optimizeClippingPlanes() { + if ( m_rimView == nullptr ) return; + if ( m_showWindowEdgeAxes ) { m_windowEdgeAxisOverlay->setDisplayCoordTransform( m_rimView->displayCoordTransform().p() );