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

ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp

@@ -365,6 +365,33 @@ const std::vector<double> RigFaultReactivationModelGenerator::partition( double
     return parts;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<int> RigFaultReactivationModelGenerator::elementKLayers( const std::vector<size_t>& cellIndexColumn )
+{
+    std::vector<int> 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<size_t> cellColumnBack;
     std::vector<size_t> cellColumnFront;
     size_t              i, j, k;
-    std::vector<int>    kLayersFront;
-    std::vector<int>    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 );
-
+    auto kLayersBack     = elementKLayers( cellColumnBack );
-    std::reverse( kLayersBack.begin(), kLayersBack.end() );
+    auto kLayersFront    = elementKLayers( cellColumnFront );
-    std::reverse( kLayersFront.begin(), kLayersFront.end() );
 
     // 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<double, cvf::Vec3d> RigFaultReactivationModelGenerator::elementLayers( cvf::StructGridInterface::FaceType face,
-                                                                                const std::vector<size_t>&         cellIndexColumn )
+                                                                                std::vector<size_t>&               cellIndexColumn )
 {
     cvf::Plane modelPlane;
     modelPlane.setFromPointAndNormal( m_startPosition, m_normal );
 
     auto cornerIndexes = faceIJCornerIndexes( face );
 
-    std::vector<int> klayers;
-
     std::map<double, cvf::Vec3d> zPositions;
 
+    std::vector<size_t> okCells;
+
     for ( auto cellIdx : cellIndexColumn )
     {
         RigCell cell    = m_grid->cell( cellIdx );
@@ -518,8 +536,19 @@ std::map<double, cvf::Vec3d> RigFaultReactivationModelGenerator::elementLayers(
         cvf::Vec3d intersect1 = lineIntersect( modelPlane, corners[cornerIndexes[0]], corners[cornerIndexes[1]] );
         cvf::Vec3d intersect2 = lineIntersect( modelPlane, corners[cornerIndexes[2]], corners[cornerIndexes[3]] );
 
+        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<double, cvf:
 
     bool first = true;
     int  k     = 0;
-    kLayers.push_back( kLayers.back() );
+
+    const int nLayers = (int)layers.size();
+    int       n       = 0;
 
     for ( auto& layer : layers )
     {
+        if ( n++ == ( nLayers - 1 ) ) break;
+
         if ( first )
         {
             prevLayer = layer;
@@ -566,7 +599,7 @@ void RigFaultReactivationModelGenerator::splitLargeLayers( std::map<double, cvf:
             for ( auto& p : points )
             {
                 additionalPoints.push_back( p );
-                newKLayers.push_back( kLayers[k] );
+                newKLayers.push_back( kLayers[k - 1] );
             }
         }
 

ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.h

@@ -72,7 +72,9 @@ protected:
     static cvf::Vec3d extrapolatePoint( cvf::Vec3d startPoint, cvf::Vec3d endPoint, double stopDepth );
     static void       splitLargeLayers( std::map<double, cvf::Vec3d>& layers, std::vector<int>& kLayers, double maxHeight );
 
-    std::map<double, cvf::Vec3d> elementLayers( cvf::StructGridInterface::FaceType face, const std::vector<size_t>& cellIndexColumn );
+    std::map<double, cvf::Vec3d> elementLayers( cvf::StructGridInterface::FaceType face, std::vector<size_t>& cellIndexColumn );
+    std::vector<int>             elementKLayers( const std::vector<size_t>& cellIndexColumn );
+
     void addFilter( QString name, std::vector<size_t> cells );
 
     size_t oppositeStartCellIndex( const std::vector<size_t> cellIndexColumn, cvf::StructGridInterface::FaceType face );

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<double> 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<double> RigGriddedPart3d::extractZValues( std::vector<cvf::Vec3d> po
     return layers;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RigGriddedPart3d::updateReservoirElementLayers( const std::vector<cvf::Vec3d>& reservoirLayers, const std::vector<int>& 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<double> RigGriddedPart3d::extractZValues( std::vector<cvf::Vec3d> po
 ///
 ///
 //--------------------------------------------------------------------------------------------------
-
 void RigGriddedPart3d::generateGeometry( const std::array<cvf::Vec3d, 12>& inputPoints,
                                          const std::vector<cvf::Vec3d>&    reservoirLayers,
                                          const std::vector<int>&           kLayers,
@@ -204,13 +238,19 @@ void RigGriddedPart3d::generateGeometry( const std::array<cvf::Vec3d, 12>& input
     std::map<Regions, std::vector<double>> 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::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<cvf::Vec3d, 12>& 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<cvf::Vec3d, 12>& 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<cvf::Vec3d, 12>& 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<int> RigGriddedPart3d::elementKLayer() const
     return m_elementKLayer;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+const std::vector<std::pair<double, double>> RigGriddedPart3d::layers( RigGriddedPart3d::ElementSets elementSet ) const
+{
+    if ( m_elementLayers.count( elementSet ) == 0 ) return {};
+    return m_elementLayers.at( elementSet );
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h

@@ -70,6 +70,7 @@ public:
     const std::map<ElementSets, std::vector<unsigned int>>& elementSets() const;
     const std::vector<int>                                  elementKLayer() const;
     const std::vector<cvf::Vec3d>                           elementCorners( size_t elementIndex ) const;
+    const std::vector<std::pair<double, double>>            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<double> extractZValues( std::vector<cvf::Vec3d> );
 
     void generateVerticalMeshlines( const std::vector<cvf::Vec3d>& cornerPoints, const std::vector<double>& horzPartition );
+    void updateReservoirElementLayers( const std::vector<cvf::Vec3d>& reservoirLayers, const std::vector<int>& kLayers );
 
 private:
     enum class Regions
@@ -105,4 +107,5 @@ private:
     std::map<Boundary, std::vector<unsigned int>>                            m_boundaryElements;
     std::map<Boundary, std::vector<unsigned int>>                            m_boundaryNodes;
     std::map<ElementSets, std::vector<unsigned int>>                         m_elementSets;
+    std::map<ElementSets, std::vector<std::pair<double, double>>>            m_elementLayers;
 };

ApplicationLibCode/UserInterface/RiuFemResultTextBuilder.cpp

@@ -181,6 +181,17 @@ QString RiuFemResultTextBuilder::geometrySelectionText( QString itemSeparator )
                     }
                 }
             }
+            std::array<cvf::Vec3d, 8> 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 );
+                }
+            }
         }
     }