Work in progress

ApplicationLibCode/FileInterface/RifReaderOpmCommon.cpp

@@ -330,7 +330,8 @@ void RifReaderOpmCommon::transferStaticNNCData( Opm::EclIO::EGrid& opmMainGrid,
             RigGridBase* grid1 = mainGrid->gridByIndex( c.grid1_Id );
             RigGridBase* grid2 = mainGrid->gridByIndex( c.grid2_Id );
 
-            RigConnection nncConnection( grid1->reservoirCellIndex( c.grid1_CellIdx - 1 ), grid2->reservoirCellIndex( c.grid2_CellIdx - 1 ) );
+            RigConnection nncConnection( grid1->localCellIndexToNative( c.grid1_CellIdx - 1 ),
+                                         grid2->localCellIndexToNative( c.grid2_CellIdx - 1 ) );
 
             nncConnections.push_back( nncConnection );
 

ApplicationLibCode/FileInterface/RifReaderOpmCommonActive.cpp

@@ -122,8 +122,6 @@ bool RifReaderOpmCommonActive::importGrid( RigMainGrid* /* mainGrid*/, RigEclips
             updateActiveCellInfo( eclipseCaseData, opmGrid, lgrGrids, activeGrid );
         }
 
-        // TODO - loop over all grids
-
         size_t anInactiveCellIndex = activeGrid->transferActiveInformation( 0,
                                                                             eclipseCaseData,
                                                                             opmGrid.totalActiveCells(),
@@ -156,7 +154,7 @@ bool RifReaderOpmCommonActive::importGrid( RigMainGrid* /* mainGrid*/, RigEclips
                               .arg( QString::fromStdString( RiaStdStringTools::formatThousandGrouping( opmGrid.totalNumberOfCells() ) ) ) );
 
         auto task = progInfo.task( "Loading Active Cell Main Grid Geometry", 1 );
-        transferActiveGeometry( opmGrid, activeGrid, eclipseCaseData );
+        transferActiveGeometry( opmGrid, opmGrid, activeGrid, activeGrid, eclipseCaseData );
 
         bool hasParentInfo = ( lgr_parent_names.size() >= (size_t)numLGRs );
 
@@ -169,7 +167,7 @@ bool RifReaderOpmCommonActive::importGrid( RigMainGrid* /* mainGrid*/, RigEclips
             RigLocalGrid* localGrid = static_cast<RigLocalGrid*>( activeGrid->gridById( lgrIdx + 1 ) );
             localGrid->setParentGrid( parentGrid );
 
-            transferActiveGeometry( opmGrid, lgrGrids[lgrIdx], mainGrid, localGrid, eclipseCaseData );
+            transferActiveGeometry( opmGrid, lgrGrids[lgrIdx], activeGrid, localGrid, eclipseCaseData );
         }
     }
 
@@ -211,51 +209,148 @@ bool RifReaderOpmCommonActive::importGrid( RigMainGrid* /* mainGrid*/, RigEclips
 
     return true;
 }
+//
+////--------------------------------------------------------------------------------------------------
+/////
+////--------------------------------------------------------------------------------------------------
+// void RifReaderOpmCommonActive::transferActiveGeometry( Opm::EclIO::EGrid&  opmMainGrid,
+//                                                        RigActiveCellGrid*  activeGrid,
+//                                                        RigEclipseCaseData* eclipseCaseData )
+//{
+//     int cellCount = opmMainGrid.totalActiveCells();
+//
+//     RigCell defaultCell;
+//     defaultCell.setHostGrid( activeGrid );
+//     for ( size_t i = 0; i < 8; i++ )
+//         defaultCell.cornerIndices()[i] = 0;
+//
+//     activeGrid->reservoirCells().resize( cellCount + 1, defaultCell );
+//     activeGrid->reservoirCells()[cellCount].setInvalid( true );
+//
+//     activeGrid->nodes().resize( ( cellCount + 1 ) * 8, cvf::Vec3d( 0, 0, 0 ) );
+//
+//     auto& riNodes = activeGrid->nodes();
+//
+//     opmMainGrid.loadData();
+//     opmMainGrid.load_grid_data();
+//
+//     const bool  isRadialGrid      = opmMainGrid.is_radial();
+//     const auto& activeMatIndexes  = opmMainGrid.active_indexes();
+//     const auto& activeFracIndexes = opmMainGrid.active_frac_indexes();
+//
+//     // Compute the center of the LGR radial grid cells for each K layer
+//     auto radialGridCenterTopLayerOpm = isRadialGrid
+//                                            ? RifOpmRadialGridTools::computeXyCenterForTopOfCells( opmMainGrid, opmMainGrid, activeGrid )
+//                                            : std::map<int, std::pair<double, double>>();
+//
+//     const bool invalidateLongPyramidCells = invalidateLongThinCells();
+//
+//     // use same mapping as resdata
+//     const size_t cellMappingECLRi[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };
+//
+// #pragma omp parallel for
+//     for ( int opmCellIndex = 0; opmCellIndex < static_cast<int>( opmMainGrid.totalNumberOfCells() ); opmCellIndex++ )
+//     {
+//         if ( ( activeMatIndexes[opmCellIndex] < 0 ) && ( activeFracIndexes[opmCellIndex] < 0 ) ) continue;
+//
+//         auto opmIJK = opmMainGrid.ijk_from_global_index( opmCellIndex );
+//
+//         double xCenterCoordOpm = 0.0;
+//         double yCenterCoordOpm = 0.0;
+//
+//         if ( isRadialGrid && radialGridCenterTopLayerOpm.contains( opmIJK[2] ) )
+//         {
+//             const auto& [xCenter, yCenter] = radialGridCenterTopLayerOpm[opmIJK[2]];
+//             xCenterCoordOpm                = xCenter;
+//             yCenterCoordOpm                = yCenter;
+//         }
+//
+//         auto     nativeIndex = activeGrid->cellIndexFromIJK( opmIJK[0], opmIJK[1], opmIJK[2] );
+//         RigCell& cell        = activeGrid->nativeCell( nativeIndex );
+//         // auto     globalIndex = activeGrid->nativeCellIndexToGlobal( nativeIndex );
+//         cell.setGridLocalCellIndex( nativeIndex );
+//         cell.setParentCellIndex( cvf::UNDEFINED_SIZE_T );
+//
+//         // corner coordinates
+//         std::array<double, 8> opmX{};
+//         std::array<double, 8> opmY{};
+//         std::array<double, 8> opmZ{};
+//         opmMainGrid.getCellCorners( opmCellIndex, opmX, opmY, opmZ );
+//
+//         // Each cell has 8 nodes, use active cell index and multiply to find first node index for cell
+//         auto riNodeStartIndex = nativeIndex * 8;
+//
+//         for ( size_t opmNodeIndex = 0; opmNodeIndex < 8; opmNodeIndex++ )
+//         {
+//             auto   riCornerIndex = cellMappingECLRi[opmNodeIndex];
+//             size_t riNodeIndex   = riNodeStartIndex + riCornerIndex;
+//
+//             auto& riNode = riNodes[riNodeIndex];
+//             riNode.x()   = opmX[opmNodeIndex] + xCenterCoordOpm;
+//             riNode.y()   = opmY[opmNodeIndex] + yCenterCoordOpm;
+//             riNode.z()   = -opmZ[opmNodeIndex];
+//
+//             cell.cornerIndices()[riCornerIndex] = riNodeIndex;
+//         }
+//
+//         if ( invalidateLongPyramidCells )
+//         {
+//             cell.setInvalid( cell.isLongPyramidCell() );
+//         }
+//     }
+//
+//     if ( riNodes.size() > 1 ) riNodes[riNodes.size() - 1] = riNodes[0];
+// }
 
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
 void RifReaderOpmCommonActive::transferActiveGeometry( Opm::EclIO::EGrid&  opmMainGrid,
+                                                       Opm::EclIO::EGrid&  opmGrid,
                                                        RigActiveCellGrid*  activeGrid,
+                                                       RigGridBase*        localGrid,
                                                        RigEclipseCaseData* eclipseCaseData )
 {
-    int cellCount = opmMainGrid.totalActiveCells();
+    int    cellCount      = opmGrid.totalActiveCells();
+    size_t cellStartIndex = activeGrid->reservoirCells().size();
+    size_t nodeStartIndex = activeGrid->nodes().size();
+
+    const bool invalidateLongPyramidCells = invalidateLongThinCells();
 
     RigCell defaultCell;
-    defaultCell.setHostGrid( activeGrid );
+    defaultCell.setHostGrid( localGrid );
     for ( size_t i = 0; i < 8; i++ )
         defaultCell.cornerIndices()[i] = 0;
 
-    activeGrid->reservoirCells().resize( cellCount + 1, defaultCell );
+    const auto newCellCount = cellStartIndex + cellCount + 1;
-    activeGrid->reservoirCells()[cellCount].setInvalid( true );
+    activeGrid->reservoirCells().resize( newCellCount, defaultCell );
-
+    activeGrid->reservoirCells()[newCellCount - 1].setInvalid( true );
-    activeGrid->nodes().resize( ( cellCount + 1 ) * 8, cvf::Vec3d( 0, 0, 0 ) );
+    activeGrid->nodes().resize( ( newCellCount ) * 8, cvf::Vec3d( 0, 0, 0 ) );
 
     auto& riNodes = activeGrid->nodes();
 
-    opmMainGrid.loadData();
+    opmGrid.loadData();
-    opmMainGrid.load_grid_data();
+    opmGrid.load_grid_data();
 
-    const bool  isRadialGrid      = opmMainGrid.is_radial();
+    const bool  isRadialGrid          = opmGrid.is_radial();
-    const auto& activeMatIndexes  = opmMainGrid.active_indexes();
+    const auto& activeMatIndexes      = opmGrid.active_indexes();
-    const auto& activeFracIndexes = opmMainGrid.active_frac_indexes();
+    const auto& activeFracIndexes     = opmGrid.active_frac_indexes();
+    const auto& gridDimension         = opmGrid.dimension();
+    const auto& hostCellGlobalIndices = opmGrid.hostCellsGlobalIndex();
 
     // Compute the center of the LGR radial grid cells for each K layer
-    auto radialGridCenterTopLayerOpm = isRadialGrid
+    auto radialGridCenterTopLayerOpm = isRadialGrid ? RifOpmRadialGridTools::computeXyCenterForTopOfCells( opmMainGrid, opmGrid, localGrid )
-                                           ? RifOpmRadialGridTools::computeXyCenterForTopOfCells( opmMainGrid, opmMainGrid, activeGrid )
+                                                    : std::map<int, std::pair<double, double>>();
-                                           : std::map<int, std::pair<double, double>>();
-
-    const bool invalidateLongPyramidCells = invalidateLongThinCells();
 
     // use same mapping as resdata
     const size_t cellMappingECLRi[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };
 
 #pragma omp parallel for
-    for ( int opmCellIndex = 0; opmCellIndex < static_cast<int>( opmMainGrid.totalNumberOfCells() ); opmCellIndex++ )
+    for ( int opmCellIndex = 0; opmCellIndex < static_cast<int>( opmGrid.totalNumberOfCells() ); opmCellIndex++ )
     {
         if ( ( activeMatIndexes[opmCellIndex] < 0 ) && ( activeFracIndexes[opmCellIndex] < 0 ) ) continue;
 
-        auto opmIJK = opmMainGrid.ijk_from_global_index( opmCellIndex );
+        auto opmIJK = opmGrid.ijk_from_global_index( opmCellIndex );
 
         double xCenterCoordOpm = 0.0;
         double yCenterCoordOpm = 0.0;
@@ -268,19 +363,27 @@ void RifReaderOpmCommonActive::transferActiveGeometry( Opm::EclIO::EGrid&  opmMa
         }
 
         auto     nativeIndex = activeGrid->cellIndexFromIJK( opmIJK[0], opmIJK[1], opmIJK[2] );
-        RigCell& cell        = activeGrid->nativeCell( nativeIndex );
+        RigCell& cell        = activeGrid->nativeCell( cellStartIndex + nativeIndex );
-        // auto     globalIndex = activeGrid->nativeCellIndexToGlobal( nativeIndex );
         cell.setGridLocalCellIndex( nativeIndex );
-        cell.setParentCellIndex( cvf::UNDEFINED_SIZE_T );
+
+        // parent cell index
+        if ( ( hostCellGlobalIndices.size() > (size_t)opmCellIndex ) && hostCellGlobalIndices[opmCellIndex] >= 0 )
+        {
+            cell.setParentCellIndex( hostCellGlobalIndices[opmCellIndex] );
+        }
+        else
+        {
+            cell.setParentCellIndex( cvf::UNDEFINED_SIZE_T );
+        }
 
         // corner coordinates
         std::array<double, 8> opmX{};
         std::array<double, 8> opmY{};
         std::array<double, 8> opmZ{};
-        opmMainGrid.getCellCorners( opmCellIndex, opmX, opmY, opmZ );
+        opmGrid.getCellCorners( opmCellIndex, opmX, opmY, opmZ );
 
         // Each cell has 8 nodes, use active cell index and multiply to find first node index for cell
-        auto riNodeStartIndex = nativeIndex * 8;
+        auto riNodeStartIndex = nodeStartIndex + nativeIndex * 8;
 
         for ( size_t opmNodeIndex = 0; opmNodeIndex < 8; opmNodeIndex++ )
         {
@@ -293,6 +396,22 @@ void RifReaderOpmCommonActive::transferActiveGeometry( Opm::EclIO::EGrid&  opmMa
             riNode.z()   = -opmZ[opmNodeIndex];
 
             cell.cornerIndices()[riCornerIndex] = riNodeIndex;
+
+            // First grid dimension is radius, check if cell are at the outer-most slice
+            if ( isRadialGrid && !hostCellGlobalIndices.empty() && ( gridDimension[0] - 1 == opmIJK[0] ) )
+            {
+                auto hostCellIndex = hostCellGlobalIndices[opmCellIndex];
+
+                RifOpmRadialGridTools::lockToHostPillars( riNode,
+                                                          opmMainGrid,
+                                                          opmGrid,
+                                                          opmIJK,
+                                                          hostCellIndex,
+                                                          opmCellIndex,
+                                                          opmNodeIndex,
+                                                          xCenterCoordOpm,
+                                                          yCenterCoordOpm );
+            }
         }
 
         if ( invalidateLongPyramidCells )

ApplicationLibCode/FileInterface/RifReaderOpmCommonActive.h

@@ -34,10 +34,10 @@ public:
 
 protected:
     bool importGrid( RigMainGrid* mainGrid, RigEclipseCaseData* caseData ) override;
-    void transferActiveGeometry( Opm::EclIO::EGrid& opmMainGrid, RigActiveCellGrid* riMainGrid, RigEclipseCaseData* caseData );
+
     void transferActiveGeometry( Opm::EclIO::EGrid&  opmMainGrid,
                                  Opm::EclIO::EGrid&  opmGrid,
-                                 RigMainGrid*        mainGrid,
+                                 RigActiveCellGrid*  activeGrid,
                                  RigGridBase*        localGrid,
                                  RigEclipseCaseData* eclipseCaseData );
 };

ApplicationLibCode/ReservoirDataModel/ResultAccessors/RigActiveCellsResultAccessor.cpp

@@ -40,7 +40,7 @@ double RigActiveCellsResultAccessor::cellScalar( size_t gridLocalCellIndex ) con
 {
     if ( m_reservoirResultValues == nullptr || m_reservoirResultValues->empty() ) return HUGE_VAL;
 
-    size_t reservoirCellIndex = m_grid->reservoirCellIndex( gridLocalCellIndex );
+    size_t reservoirCellIndex = m_grid->localCellIndexToNative( gridLocalCellIndex );
     size_t resultValueIndex   = m_activeCellInfo->cellResultIndex( reservoirCellIndex );
     if ( resultValueIndex == cvf::UNDEFINED_SIZE_T ) return HUGE_VAL;
 

ApplicationLibCode/ReservoirDataModel/ResultAccessors/RigAllGridCellsResultAccessor.cpp

@@ -39,7 +39,7 @@ double RigAllGridCellsResultAccessor::cellScalar( size_t gridLocalCellIndex ) co
 {
     if ( m_reservoirResultValues->empty() ) return HUGE_VAL;
 
-    size_t reservoirCellIndex = m_grid->reservoirCellIndex( gridLocalCellIndex );
+    size_t reservoirCellIndex = m_grid->localCellIndexToNative( gridLocalCellIndex );
     CVF_TIGHT_ASSERT( reservoirCellIndex < m_reservoirResultValues->size() );
 
     return m_reservoirResultValues->at( reservoirCellIndex );

ApplicationLibCode/ReservoirDataModel/RigCaseCellResultCalculator.cpp

@@ -185,7 +185,7 @@ bool RigCaseCellResultCalculator::computeDifference( RigEclipseCaseData*
 #pragma omp parallel for
             for ( long localGridCellIdx = 0; localGridCellIdx < static_cast<long>( grid->cellCount() ); localGridCellIdx++ )
             {
-                size_t reservoirCellIndex = grid->reservoirCellIndex( localGridCellIdx );
+                size_t reservoirCellIndex = grid->localNativeToNative( localGridCellIdx );
                 if ( activeCellInfo->isActive( reservoirCellIndex ) )
                 {
                     double sourceVal = sourceResultAccessor->cellScalar( localGridCellIdx );
@@ -290,7 +290,7 @@ bool RigCaseCellResultCalculator::computeDivideByCellFaceArea( RigMainGrid*
 #pragma omp parallel for
             for ( int localGridCellIdx = 0; localGridCellIdx < static_cast<int>( grid->cellCount() ); localGridCellIdx++ )
             {
-                const size_t reservoirCellIndex = grid->reservoirCellIndex( localGridCellIdx );
+                const size_t reservoirCellIndex = grid->localNativeToNative( localGridCellIdx );
                 if ( activeCellInfo->isActive( reservoirCellIndex ) )
                 {
                     double sourceVal = sourceResultAccessor->cellScalar( localGridCellIdx );

ApplicationLibCode/ReservoirDataModel/RigGridBase.cpp

@@ -603,9 +603,7 @@ bool RigGridCellFaceVisibilityFilter::isFaceVisible( size_t
 
     // Do not show cell geometry if a fault is present to avoid z fighting between surfaces
     // It will always be a better solution to avoid geometry creation instead of part priority and polygon offset
-    const auto      cellIndex     = m_grid->nativeCellIndexToGlobal( nativeCellIndex );
+    const RigFault* fault = m_grid->mainGrid()->findFaultFromCellIndexAndCellFace( nativeCellIndex, face );
-    size_t          resvCellIndex = m_grid->reservoirCellIndex( cellIndex );
-    const RigFault* fault         = m_grid->mainGrid()->findFaultFromCellIndexAndCellFace( resvCellIndex, face );
     if ( fault )
     {
         return false;

ApplicationLibCode/ReservoirDataModel/RigReservoirBuilder.cpp

@@ -109,7 +109,7 @@ void RigReservoirBuilder::createGridsAndCells( RigEclipseCaseData* eclipseCase )
         eclipseCase->mainGrid()->addLocalGrid( localGrid );
         localGrid->setParentGrid( eclipseCase->mainGrid() );
 
-        localGrid->setIndexToStartOfCells( mainGridNodes.size() / 8 );
+        localGrid->setIndexToGlobalStartOfCells( mainGridNodes.size() / 8 );
         cvf::Vec3st gridPointDimensions( lgr.m_singleCellRefinementFactors.x() *
                                                  ( lgr.m_mainGridMaxCellPosition.x() - lgr.m_mainGridMinCellPosition.x() + 1 ) +
                                              1,

ApplicationLibCode/ReservoirDataModel/RigResultModifier.h

@@ -45,10 +45,10 @@ public:
 
     void setCellScalar( size_t gridLocalCellIndex, double scalarValue ) override
     {
-        size_t reservoirCellIndex = m_grid->reservoirCellIndex( gridLocalCellIndex );
+        size_t nativeCellIndex = m_grid->localCellIndexToNative( gridLocalCellIndex );
-        CVF_TIGHT_ASSERT( reservoirCellIndex < m_reservoirResultValues->size() );
+        CVF_TIGHT_ASSERT( nativeCellIndex < m_reservoirResultValues->size() );
 
-        ( *m_reservoirResultValues )[reservoirCellIndex] = scalarValue;
+        ( *m_reservoirResultValues )[nativeCellIndex] = scalarValue;
     }
 
 private:
@@ -71,8 +71,8 @@ public:
 
     void setCellScalar( size_t gridLocalCellIndex, double scalarValue ) override
     {
-        size_t reservoirCellIndex = m_grid->reservoirCellIndex( gridLocalCellIndex );
+        size_t nativeCellIndex  = m_grid->localCellIndexToNative( gridLocalCellIndex );
-        size_t resultValueIndex   = m_activeCellInfo->cellResultIndex( reservoirCellIndex );
+        size_t resultValueIndex = m_activeCellInfo->cellResultIndex( nativeCellIndex );
 
         CVF_TIGHT_ASSERT( m_reservoirResultValues != nullptr && resultValueIndex < m_reservoirResultValues->size() );