diff --git a/ApplicationLibCode/Application/Tools/RiaCellDividingTools.cpp b/ApplicationLibCode/Application/Tools/RiaCellDividingTools.cpp index 0838a7086b..f93d307054 100644 --- a/ApplicationLibCode/Application/Tools/RiaCellDividingTools.cpp +++ b/ApplicationLibCode/Application/Tools/RiaCellDividingTools.cpp @@ -32,6 +32,11 @@ //-------------------------------------------------------------------------------------------------- std::vector splitLine( cvf::Vec3d ptStart, cvf::Vec3d ptEnd, size_t partCount ); +//-------------------------------------------------------------------------------------------------- +/// Splits a line using cumulative fractions instead of equal parts +//-------------------------------------------------------------------------------------------------- +std::vector splitLineNonUniform( cvf::Vec3d ptStart, cvf::Vec3d ptEnd, const std::vector& cumulativeFractions ); + //-------------------------------------------------------------------------------------------------- /// Calculates all points on a face described by edge points from all four edges. /// The result is a grid of points including the given edge points @@ -50,6 +55,15 @@ std::vector> calcFacePoints( const std::vector edgeYPtsLow, const std::vector edgeYPtsHigh ); +//-------------------------------------------------------------------------------------------------- +/// Calculates face points using non-uniform X-direction fractions for interior rows +//-------------------------------------------------------------------------------------------------- +std::vector> calcFacePointsNonUniform( const std::vector& edgeXPtsLow, + const std::vector& edgeXPtsHigh, + const std::vector& edgeYPtsLow, + const std::vector& edgeYPtsHigh, + const std::vector& xCumulativeFractions ); + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -204,3 +218,135 @@ std::vector> calcFacePoints( const std::vector splitLineNonUniform( cvf::Vec3d ptStart, cvf::Vec3d ptEnd, const std::vector& cumulativeFractions ) +{ + std::vector pts = { ptStart }; + + for ( size_t i = 0; i < cumulativeFractions.size() - 1; ++i ) + { + double f = cumulativeFractions[i]; + pts.push_back( cvf::Vec3d( ptStart.x() + ( ptEnd.x() - ptStart.x() ) * f, + ptStart.y() + ( ptEnd.y() - ptStart.y() ) * f, + ptStart.z() + ( ptEnd.z() - ptStart.z() ) * f ) ); + } + pts.push_back( ptEnd ); + return pts; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector> calcFacePointsNonUniform( const std::vector& edgeXPtsLow, + const std::vector& edgeXPtsHigh, + const std::vector& edgeYPtsLow, + const std::vector& edgeYPtsHigh, + const std::vector& xCumulativeFractions ) +{ + CVF_ASSERT( edgeYPtsLow.size() == edgeYPtsHigh.size() ); + + size_t ySize = edgeYPtsLow.size(); + + std::vector> pts; + + pts.push_back( edgeXPtsLow ); + + for ( size_t y = 1; y < ySize - 1; y++ ) + { + auto interiorPts = splitLineNonUniform( edgeYPtsLow[y], edgeYPtsHigh[y], xCumulativeFractions ); + pts.push_back( interiorPts ); + } + + pts.push_back( edgeXPtsHigh ); + return pts; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RiaCellDividingTools::createHexCornerCoords( const std::array& mainCellCorners, + const std::vector& cumulativeFractionsX, + const std::vector& cumulativeFractionsY, + const std::vector& cumulativeFractionsZ ) +{ + size_t nx = cumulativeFractionsX.size(); + size_t ny = cumulativeFractionsY.size(); + size_t nz = cumulativeFractionsZ.size(); + + std::array, 12> edgeCorners = { + std::make_pair( 0, 1 ), + std::make_pair( 3, 2 ), + std::make_pair( 4, 5 ), + std::make_pair( 7, 6 ), // X + std::make_pair( 0, 3 ), + std::make_pair( 4, 7 ), + std::make_pair( 1, 2 ), + std::make_pair( 5, 6 ), // Y + std::make_pair( 0, 4 ), + std::make_pair( 1, 5 ), + std::make_pair( 3, 7 ), + std::make_pair( 2, 6 ), // Z + }; + + // Split edges using appropriate cumulative fractions per dimension + const std::vector* dimFractions[3] = { &cumulativeFractionsX, &cumulativeFractionsY, &cumulativeFractionsZ }; + + std::array, 12> edgePoints; + for ( int i = 0; i < 12; i++ ) + { + int dimIndex = i / 4; + edgePoints[i] = + splitLineNonUniform( mainCellCorners[edgeCorners[i].first], mainCellCorners[edgeCorners[i].second], *dimFractions[dimIndex] ); + } + + std::vector>> nodes; + nodes.reserve( ( nx + 1 ) * ( ny + 1 ) * ( nz + 1 ) ); + + auto xyFacePtsLow = calcFacePointsNonUniform( edgePoints[0], edgePoints[1], edgePoints[4], edgePoints[6], cumulativeFractionsX ); + auto xyFacePtsHigh = calcFacePointsNonUniform( edgePoints[2], edgePoints[3], edgePoints[5], edgePoints[7], cumulativeFractionsX ); + auto yzFacePtsLow = calcFacePointsNonUniform( edgePoints[4], edgePoints[5], edgePoints[8], edgePoints[10], cumulativeFractionsY ); + auto yzFacePtsHigh = calcFacePointsNonUniform( edgePoints[6], edgePoints[7], edgePoints[9], edgePoints[11], cumulativeFractionsY ); + auto xzFacePtsLow = calcFacePointsNonUniform( edgePoints[0], edgePoints[2], edgePoints[8], edgePoints[9], cumulativeFractionsX ); + auto xzFacePtsHigh = calcFacePointsNonUniform( edgePoints[1], edgePoints[3], edgePoints[10], edgePoints[11], cumulativeFractionsX ); + + nodes.push_back( xyFacePtsLow ); + + for ( size_t z = 1; z < nz; z++ ) + { + auto xyFacePoints = + calcFacePointsNonUniform( xzFacePtsLow[z], xzFacePtsHigh[z], yzFacePtsLow[z], yzFacePtsHigh[z], cumulativeFractionsX ); + nodes.push_back( xyFacePoints ); + } + + nodes.push_back( xyFacePtsHigh ); + + std::vector coords; + coords.reserve( nx * ny * nz * 8 ); + + for ( size_t z = 1; z < nz + 1; z++ ) + { + for ( size_t y = 1; y < ny + 1; y++ ) + { + for ( size_t x = 1; x < nx + 1; x++ ) + { + std::array cs; + + cs[0] = nodes[z - 1][y - 1][x - 1]; + cs[1] = nodes[z - 1][y - 1][x]; + cs[2] = nodes[z - 1][y][x]; + cs[3] = nodes[z - 1][y][x - 1]; + + cs[4] = nodes[z][y - 1][x - 1]; + cs[5] = nodes[z][y - 1][x]; + cs[6] = nodes[z][y][x]; + cs[7] = nodes[z][y][x - 1]; + + coords.insert( coords.end(), cs.begin(), cs.end() ); + } + } + } + return coords; +} diff --git a/ApplicationLibCode/Application/Tools/RiaCellDividingTools.h b/ApplicationLibCode/Application/Tools/RiaCellDividingTools.h index 7e63d18c03..d8ea4276a5 100644 --- a/ApplicationLibCode/Application/Tools/RiaCellDividingTools.h +++ b/ApplicationLibCode/Application/Tools/RiaCellDividingTools.h @@ -31,5 +31,11 @@ class RiaCellDividingTools public: static std::vector createHexCornerCoords( std::array mainCellCorners, size_t nx, size_t ny, size_t nz ); + // Non-uniform overload: takes cumulative fractions per dimension instead of integer counts + static std::vector createHexCornerCoords( const std::array& mainCellCorners, + const std::vector& cumulativeFractionsX, + const std::vector& cumulativeFractionsY, + const std::vector& cumulativeFractionsZ ); + static double computeFlowDistance( const std::array& cellVertices, const cvf::Vec3d& areaCenter ); }; diff --git a/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelFeature.cpp b/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelFeature.cpp index 1512d6d5da..975e819540 100644 --- a/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelFeature.cpp +++ b/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelFeature.cpp @@ -101,6 +101,11 @@ void RicExportSectorModelFeature::doExport( RicExportSectorModelUi* exportSettin settings.setMax( exportSettings->max() ); settings.setRefinement( exportSettings->refinement() ); + if ( exportSettings->hasNonUniformRefinement() ) + { + settings.setNonUniformRefinement( exportSettings->nonUniformRefinement() ); + } + std::vector keywordsToRemove; for ( auto& s : exportSettings->keywordsToRemove() ) { diff --git a/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelUi.cpp b/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelUi.cpp index f391b419a6..4a7b30f17b 100644 --- a/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelUi.cpp +++ b/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelUi.cpp @@ -44,6 +44,17 @@ #include +namespace caf +{ +template <> +void AppEnum::setUp() +{ + addItem( RicExportSectorModelUi::UNIFORM, "UNIFORM", "Uniform" ); + addItem( RicExportSectorModelUi::NON_UNIFORM, "NON_UNIFORM", "Non-Uniform" ); + setDefault( RicExportSectorModelUi::UNIFORM ); +} +} // namespace caf + CAF_PDM_SOURCE_INIT( RicExportSectorModelUi, "RicExportSectorModelUi" ); //-------------------------------------------------------------------------------------------------- @@ -99,6 +110,20 @@ RicExportSectorModelUi::RicExportSectorModelUi() m_refinementCountJ.setRange( 1, 10 ); m_refinementCountK.setRange( 1, 10 ); + CAF_PDM_InitFieldNoDefault( &m_refinementMode, "RefinementMode", "Refinement Mode" ); + m_refinementMode.uiCapability()->setUiEditorTypeName( caf::PdmUiRadioButtonEditor::uiEditorTypeName() ); + + CAF_PDM_InitField( &m_nonUniformDirection, "NonUniformDirection", 0, "Direction" ); + m_nonUniformDirection.setRange( 0, 2 ); + + CAF_PDM_InitField( &m_nonUniformRangeStart, "NonUniformRangeStart", 1, "Cell Range Start" ); + m_nonUniformRangeStart.setMinValue( 1 ); + + CAF_PDM_InitField( &m_nonUniformRangeEnd, "NonUniformRangeEnd", 1, "Cell Range End" ); + m_nonUniformRangeEnd.setMinValue( 1 ); + + CAF_PDM_InitField( &m_nonUniformIntervals, "NonUniformIntervals", QString( "0.5, 0.5" ), "Fractional Widths" ); + CAF_PDM_InitFieldNoDefault( &m_bcpropKeywords, "BcpropKeywords", "BCPROP Keywords" ); m_bcpropKeywords.uiCapability()->setUiEditorTypeName( caf::PdmUiTableViewEditor::uiEditorTypeName() ); m_bcpropKeywords.uiCapability()->setUiLabelPosition( caf::PdmUiItemInfo::LabelPosition::HIDDEN ); @@ -241,13 +266,35 @@ void RicExportSectorModelUi::defineUiOrdering( QString uiConfigName, caf::PdmUiO uiOrdering.add( &m_refineGrid ); uiOrdering.addNewLabel( "" ); - uiOrdering.add( &m_refinementCountI, { .newRow = true, .totalColumnSpan = 2, .leftLabelColumnSpan = 1 } ); - uiOrdering.appendToRow( &m_refinementCountJ ); - uiOrdering.appendToRow( &m_refinementCountK ); + uiOrdering.add( &m_refinementMode ); + uiOrdering.addNewLabel( "" ); - m_refinementCountI.uiCapability()->setUiReadOnly( !m_refineGrid() ); - m_refinementCountJ.uiCapability()->setUiReadOnly( !m_refineGrid() ); - m_refinementCountK.uiCapability()->setUiReadOnly( !m_refineGrid() ); + bool isEnabled = m_refineGrid(); + + if ( m_refinementMode() == UNIFORM ) + { + uiOrdering.add( &m_refinementCountI, { .newRow = true, .totalColumnSpan = 2, .leftLabelColumnSpan = 1 } ); + uiOrdering.appendToRow( &m_refinementCountJ ); + uiOrdering.appendToRow( &m_refinementCountK ); + + m_refinementCountI.uiCapability()->setUiReadOnly( !isEnabled ); + m_refinementCountJ.uiCapability()->setUiReadOnly( !isEnabled ); + m_refinementCountK.uiCapability()->setUiReadOnly( !isEnabled ); + } + else + { + uiOrdering.add( &m_nonUniformDirection ); + uiOrdering.add( &m_nonUniformRangeStart ); + uiOrdering.add( &m_nonUniformRangeEnd ); + uiOrdering.add( &m_nonUniformIntervals ); + + m_nonUniformDirection.uiCapability()->setUiReadOnly( !isEnabled ); + m_nonUniformRangeStart.uiCapability()->setUiReadOnly( !isEnabled ); + m_nonUniformRangeEnd.uiCapability()->setUiReadOnly( !isEnabled ); + m_nonUniformIntervals.uiCapability()->setUiReadOnly( !isEnabled ); + } + + m_refinementMode.uiCapability()->setUiReadOnly( !isEnabled ); } else if ( uiConfigName == m_pageNames[WizardPageEnum::BoundaryConditions] ) { @@ -466,13 +513,65 @@ void RicExportSectorModelUi::setMax( const caf::VecIjk0& max ) //-------------------------------------------------------------------------------------------------- cvf::Vec3st RicExportSectorModelUi::refinement() const { - if ( !m_refineGrid() ) + if ( !m_refineGrid() || m_refinementMode() == NON_UNIFORM ) { return cvf::Vec3st( 1, 1, 1 ); } return cvf::Vec3st( m_refinementCountI(), m_refinementCountJ(), m_refinementCountK() ); } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RicExportSectorModelUi::RefinementMode RicExportSectorModelUi::refinementMode() const +{ + if ( !m_refineGrid() ) return UNIFORM; + return m_refinementMode(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigNonUniformRefinement RicExportSectorModelUi::nonUniformRefinement() const +{ + cvf::Vec3st sectorSize( max().x() - min().x() + 1, max().y() - min().y() + 1, max().z() - min().z() + 1 ); + + RigNonUniformRefinement result( sectorSize ); + + if ( !m_refineGrid() || m_refinementMode() != NON_UNIFORM ) return result; + + // Parse the fractional widths from the text field + QStringList parts = m_nonUniformIntervals().split( ",", Qt::SkipEmptyParts ); + std::vector widths; + for ( const auto& part : parts ) + { + bool ok = false; + double value = part.trimmed().toDouble( &ok ); + if ( ok && value > 0.0 ) widths.push_back( value ); + } + + if ( widths.empty() ) return result; + + auto dim = static_cast( m_nonUniformDirection() ); + int start = std::max( 0, m_nonUniformRangeStart() - 1 ); // Convert to 0-based, clamp + int end = std::min( static_cast( result.sectorSize( dim ) ) - 1, m_nonUniformRangeEnd() - 1 ); + + if ( start > end ) return result; + + // Distribute the widths across the combined range of cells + result.distributeWidthsAcrossCells( dim, start, end, widths ); + + return result; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RicExportSectorModelUi::hasNonUniformRefinement() const +{ + return m_refineGrid() && m_refinementMode() == NON_UNIFORM; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelUi.h b/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelUi.h index 927dfba60b..d8daa8c887 100644 --- a/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelUi.h +++ b/ApplicationLibCode/Commands/ExportCommands/RicExportSectorModelUi.h @@ -25,6 +25,7 @@ #include "cafVecIjk.h" #include "RiaModelExportDefines.h" +#include "RigNonUniformRefinement.h" #include #include @@ -43,8 +44,15 @@ class RicExportSectorModelUi : public caf::PdmObject { CAF_PDM_HEADER_INIT; + enum RefinementMode + { + UNIFORM, + NON_UNIFORM + }; + using GridBoxSelectionEnum = caf::AppEnum; using BoundaryConditionEnum = caf::AppEnum; + using RefinementModeEnum = caf::AppEnum; public: RicExportSectorModelUi(); @@ -63,6 +71,10 @@ public: cvf::Vec3st refinement() const; + RefinementMode refinementMode() const; + RigNonUniformRefinement nonUniformRefinement() const; + bool hasNonUniformRefinement() const; + std::vector keywordsToRemove() const; std::vector bcpropKeywords() const; @@ -126,6 +138,12 @@ private: caf::PdmField m_refinementCountJ; caf::PdmField m_refinementCountK; + caf::PdmField m_refinementMode; + caf::PdmField m_nonUniformDirection; // 0=I, 1=J, 2=K + caf::PdmField m_nonUniformRangeStart; + caf::PdmField m_nonUniformRangeEnd; + caf::PdmField m_nonUniformIntervals; + caf::PdmField m_boundaryCondition; caf::PdmChildArrayField m_bcpropKeywords; caf::PdmField m_porvMultiplier; diff --git a/ApplicationLibCode/FileInterface/RifEclipseInputFileTools.cpp b/ApplicationLibCode/FileInterface/RifEclipseInputFileTools.cpp index 0a7e5f99b4..074887a7da 100644 --- a/ApplicationLibCode/FileInterface/RifEclipseInputFileTools.cpp +++ b/ApplicationLibCode/FileInterface/RifEclipseInputFileTools.cpp @@ -335,6 +335,93 @@ std::expected, std::string> RifEclipseInputFileTools::extrac return filteredResults; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::expected, std::string> RifEclipseInputFileTools::extractKeywordData( RigEclipseCaseData* eclipseCase, + const QString& keyword, + const cvf::Vec3st& min, + const cvf::Vec3st& maxIn, + const RigNonUniformRefinement& nonUniformRefinement ) +{ + RigCaseCellResultsData* cellResultsData = eclipseCase->results( RiaDefines::PorosityModelType::MATRIX_MODEL ); + RigActiveCellInfo* activeCells = cellResultsData->activeCellInfo(); + RigMainGrid* mainGrid = eclipseCase->mainGrid(); + + cvf::Vec3st max = maxIn; + if ( max.isUndefined() ) + { + max = cvf::Vec3st( mainGrid->cellCountI() - 1, mainGrid->cellCountJ() - 1, mainGrid->cellCountK() - 1 ); + } + + auto allResultAddresses = cellResultsData->existingResults(); + auto findResultAddress = [&allResultAddresses]( const QString& keyword ) -> RigEclipseResultAddress + { + for ( auto it = allResultAddresses.rbegin(); it != allResultAddresses.rend(); ++it ) + { + if ( it->resultName() == keyword ) return *it; + } + return {}; + }; + + RigEclipseResultAddress resAddr = findResultAddress( keyword ); + if ( !cellResultsData->hasResultEntry( resAddr ) ) + return std::unexpected( "Keyword '" + keyword.toStdString() + "' not found in results" ); + + if ( !cellResultsData->ensureKnownResultLoaded( resAddr ) ) + return std::unexpected( "Keyword '" + keyword.toStdString() + "' result not loaded." ); + + if ( cellResultsData->cellScalarResults( resAddr ).empty() ) + return std::unexpected( "Keyword '" + keyword.toStdString() + "' result empty." ); + + std::vector resultValues = cellResultsData->cellScalarResults( resAddr )[0]; + if ( resultValues.empty() ) return std::unexpected( "Keyword '" + keyword.toStdString() + "' has no data" ); + + double defaultExportValue = RiaResultNames::isCategoryResult( keyword ) ? 1.0 : 0.0; + + RiaDefines::PorosityModelType porosityModel = RiaDefines::PorosityModelType::MATRIX_MODEL; + auto resultAcc = RigResultAccessorFactory::createFromResultAddress( eclipseCase, 0, porosityModel, 0, resAddr ); + + std::vector filteredResults; + size_t sectorNi = max.x() - min.x() + 1; + size_t sectorNj = max.y() - min.y() + 1; + size_t sectorNk = max.z() - min.z() + 1; + + for ( size_t origK = 0; origK < sectorNk; ++origK ) + { + size_t mainK = min.z() + origK; + size_t subCountK = nonUniformRefinement.subcellCount( RigNonUniformRefinement::DimK, origK ); + for ( size_t rk = 0; rk < subCountK; ++rk ) + { + for ( size_t origJ = 0; origJ < sectorNj; ++origJ ) + { + size_t mainJ = min.y() + origJ; + size_t subCountJ = nonUniformRefinement.subcellCount( RigNonUniformRefinement::DimJ, origJ ); + for ( size_t rj = 0; rj < subCountJ; ++rj ) + { + for ( size_t origI = 0; origI < sectorNi; ++origI ) + { + size_t mainI = min.x() + origI; + size_t reservoirCellIndex = mainGrid->cellIndexFromIJK( mainI, mainJ, mainK ); + size_t resIndex = activeCells->cellResultIndex( reservoirCellIndex ); + size_t subCountI = nonUniformRefinement.subcellCount( RigNonUniformRefinement::DimI, origI ); + + double value = ( resIndex != cvf::UNDEFINED_SIZE_T ) ? resultAcc->cellScalarGlobIdx( reservoirCellIndex ) + : defaultExportValue; + + for ( size_t ri = 0; ri < subCountI; ++ri ) + { + filteredResults.push_back( value ); + } + } + } + } + } + } + + return filteredResults; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/FileInterface/RifEclipseInputFileTools.h b/ApplicationLibCode/FileInterface/RifEclipseInputFileTools.h index c0b0df4d41..b6291a4753 100644 --- a/ApplicationLibCode/FileInterface/RifEclipseInputFileTools.h +++ b/ApplicationLibCode/FileInterface/RifEclipseInputFileTools.h @@ -25,6 +25,7 @@ #include "RiaDefines.h" #include "RigFault.h" +#include "RigNonUniformRefinement.h" #include "ert/ecl/ecl_kw.h" @@ -82,6 +83,12 @@ public: const cvf::Vec3st& max = cvf::Vec3st::UNDEFINED, const cvf::Vec3st& refinement = cvf::Vec3st( 1, 1, 1 ) ); + static std::expected, std::string> extractKeywordData( RigEclipseCaseData* eclipseCase, + const QString& keyword, + const cvf::Vec3st& min, + const cvf::Vec3st& max, + const RigNonUniformRefinement& nonUniformRefinement ); + static void saveFault( QString completeFilename, const RigMainGrid* mainGrid, const std::vector& faultFaces, diff --git a/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake b/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake index 126fb93136..76dbfa6897 100644 --- a/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake +++ b/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake @@ -88,6 +88,7 @@ set(SOURCE_GROUP_HEADER_FILES ${CMAKE_CURRENT_LIST_DIR}/RigDoglegTools.h ${CMAKE_CURRENT_LIST_DIR}/RigResdataGridConverter.h ${CMAKE_CURRENT_LIST_DIR}/RigGridExportAdapter.h + ${CMAKE_CURRENT_LIST_DIR}/RigNonUniformRefinement.h ) set(SOURCE_GROUP_SOURCE_FILES @@ -175,6 +176,7 @@ set(SOURCE_GROUP_SOURCE_FILES ${CMAKE_CURRENT_LIST_DIR}/RigDoglegTools.cpp ${CMAKE_CURRENT_LIST_DIR}/RigResdataGridConverter.cpp ${CMAKE_CURRENT_LIST_DIR}/RigGridExportAdapter.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigNonUniformRefinement.cpp ) list(APPEND CODE_HEADER_FILES ${SOURCE_GROUP_HEADER_FILES}) diff --git a/ApplicationLibCode/ReservoirDataModel/RigGridExportAdapter.cpp b/ApplicationLibCode/ReservoirDataModel/RigGridExportAdapter.cpp index 110e2a8faa..23ccf327e8 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGridExportAdapter.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigGridExportAdapter.cpp @@ -26,6 +26,7 @@ #include "RigCell.h" #include "RigEclipseCaseData.h" #include "RigMainGrid.h" +#include "RigNonUniformRefinement.h" #include "cafVecIjk.h" #include "cvfAssert.h" @@ -64,6 +65,35 @@ static std::array generateRefinedCellCorners( const std::array generateRefinedCellCornersNonUniform( const std::array& originalCorners, + const std::vector& cumulativeFractionsX, + const std::vector& cumulativeFractionsY, + const std::vector& cumulativeFractionsZ, + size_t subI, + size_t subJ, + size_t subK ) +{ + auto allRefinedCorners = + RiaCellDividingTools::createHexCornerCoords( originalCorners, cumulativeFractionsX, cumulativeFractionsY, cumulativeFractionsZ ); + + size_t nx = cumulativeFractionsX.size(); + size_t ny = cumulativeFractionsY.size(); + + size_t subcellIndex = subK * ( nx * ny ) + subJ * nx + subI; + size_t cornerStartIndex = subcellIndex * 8; + + std::array refinedCorners; + for ( size_t i = 0; i < 8; ++i ) + { + refinedCorners[i] = allRefinedCorners[cornerStartIndex + i]; + } + + return refinedCorners; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -99,6 +129,46 @@ RigGridExportAdapter::RigGridExportAdapter( RigEclipseCaseData* eclipseCase, m_refinedNI = ( maxActual.x() - min.x() + 1 ) * refinement.x(); m_refinedNJ = ( maxActual.y() - min.y() + 1 ) * refinement.y(); m_refinedNK = ( maxActual.z() - min.z() + 1 ) * refinement.z(); + + // Create equivalent non-uniform refinement for internal use + cvf::Vec3st sectorSize( maxActual.x() - min.x() + 1, maxActual.y() - min.y() + 1, maxActual.z() - min.z() + 1 ); + m_nonUniformRefinement = RigNonUniformRefinement::fromUniform( refinement, sectorSize ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigGridExportAdapter::RigGridExportAdapter( RigEclipseCaseData* eclipseCase, + const cvf::Vec3st& min, + const cvf::Vec3st& max, + const RigNonUniformRefinement& nonUniformRefinement, + const cvf::UByteArray* cellVisibilityOverrideForActnum ) + : m_mainGrid( nullptr ) + , m_activeCellInfo( nullptr ) + , m_cellVisibilityOverride( cellVisibilityOverrideForActnum ) + , m_min( min ) + , m_max( max ) + , m_refinement( 1, 1, 1 ) + , m_nonUniformRefinement( nonUniformRefinement ) + , m_refinedNI( 0 ) + , m_refinedNJ( 0 ) + , m_refinedNK( 0 ) +{ + CVF_ASSERT( eclipseCase ); + + m_mainGrid = eclipseCase->mainGrid(); + m_activeCellInfo = eclipseCase->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL ); + + CVF_ASSERT( m_mainGrid ); + CVF_ASSERT( m_activeCellInfo ); + + cvf::Vec3st maxActual = + max.isUndefined() ? cvf::Vec3st( m_mainGrid->cellCountI() - 1, m_mainGrid->cellCountJ() - 1, m_mainGrid->cellCountK() - 1 ) : max; + m_max = maxActual; + + m_refinedNI = m_nonUniformRefinement.totalRefinedCount( RigNonUniformRefinement::DimI ); + m_refinedNJ = m_nonUniformRefinement.totalRefinedCount( RigNonUniformRefinement::DimJ ); + m_refinedNK = m_nonUniformRefinement.totalRefinedCount( RigNonUniformRefinement::DimK ); } //-------------------------------------------------------------------------------------------------- @@ -157,14 +227,33 @@ std::array RigGridExportAdapter::getFaceCorners( size_t i, size_t // Apply coordinate transformations if needed applyCoordinateTransformation( originalCorners ); - // Generate refined cell corners using ResInsight's refinement algorithm - auto refinedCorners = generateRefinedCellCorners( originalCorners, - m_refinement.x(), - m_refinement.y(), - m_refinement.z(), - mapping.subI, - mapping.subJ, - mapping.subK ); + // Generate refined cell corners using the appropriate refinement algorithm + std::array refinedCorners; + if ( m_nonUniformRefinement.hasNonUniformRefinement() ) + { + size_t sectorI = mapping.originalI - m_min.x(); + size_t sectorJ = mapping.originalJ - m_min.y(); + size_t sectorK = mapping.originalK - m_min.z(); + + refinedCorners = + generateRefinedCellCornersNonUniform( originalCorners, + m_nonUniformRefinement.cumulativeFractions( RigNonUniformRefinement::DimI, sectorI ), + m_nonUniformRefinement.cumulativeFractions( RigNonUniformRefinement::DimJ, sectorJ ), + m_nonUniformRefinement.cumulativeFractions( RigNonUniformRefinement::DimK, sectorK ), + mapping.subI, + mapping.subJ, + mapping.subK ); + } + else + { + refinedCorners = generateRefinedCellCorners( originalCorners, + m_refinement.x(), + m_refinement.y(), + m_refinement.z(), + mapping.subI, + mapping.subJ, + mapping.subK ); + } // Extract the 4 corners for the requested face from the refined corners std::array faceCorners; @@ -299,10 +388,23 @@ std::array RigGridExportAdapter::getOriginalCellCorners( size_t o std::array RigGridExportAdapter::getRefinedCellCorners( size_t origI, size_t origJ, size_t origK, size_t subI, size_t subJ, size_t subK ) const { - // Get original cell corners first std::array originalCorners = getOriginalCellCorners( origI, origJ, origK ); - // Generate refined subcell corners + if ( m_nonUniformRefinement.hasNonUniformRefinement() ) + { + size_t sectorI = origI - m_min.x(); + size_t sectorJ = origJ - m_min.y(); + size_t sectorK = origK - m_min.z(); + + return generateRefinedCellCornersNonUniform( originalCorners, + m_nonUniformRefinement.cumulativeFractions( RigNonUniformRefinement::DimI, sectorI ), + m_nonUniformRefinement.cumulativeFractions( RigNonUniformRefinement::DimJ, sectorJ ), + m_nonUniformRefinement.cumulativeFractions( RigNonUniformRefinement::DimK, sectorK ), + subI, + subJ, + subK ); + } + return generateRefinedCellCorners( originalCorners, m_refinement.x(), m_refinement.y(), m_refinement.z(), subI, subJ, subK ); } @@ -343,15 +445,29 @@ RigGridExportAdapter::CellMapping RigGridExportAdapter::mapRefinedToOriginal( si { CellMapping mapping; - // Calculate which original cell this refined cell belongs to - mapping.originalI = m_min.x() + refinedI / m_refinement.x(); - mapping.originalJ = m_min.y() + refinedJ / m_refinement.y(); - mapping.originalK = m_min.z() + refinedK / m_refinement.z(); + if ( m_nonUniformRefinement.hasNonUniformRefinement() ) + { + auto [origI, subI] = m_nonUniformRefinement.mapRefinedToOriginal( RigNonUniformRefinement::DimI, refinedI ); + auto [origJ, subJ] = m_nonUniformRefinement.mapRefinedToOriginal( RigNonUniformRefinement::DimJ, refinedJ ); + auto [origK, subK] = m_nonUniformRefinement.mapRefinedToOriginal( RigNonUniformRefinement::DimK, refinedK ); - // Calculate subcell indices within the original cell - mapping.subI = refinedI % m_refinement.x(); - mapping.subJ = refinedJ % m_refinement.y(); - mapping.subK = refinedK % m_refinement.z(); + mapping.originalI = m_min.x() + origI; + mapping.originalJ = m_min.y() + origJ; + mapping.originalK = m_min.z() + origK; + mapping.subI = subI; + mapping.subJ = subJ; + mapping.subK = subK; + } + else + { + mapping.originalI = m_min.x() + refinedI / m_refinement.x(); + mapping.originalJ = m_min.y() + refinedJ / m_refinement.y(); + mapping.originalK = m_min.z() + refinedK / m_refinement.z(); + + mapping.subI = refinedI % m_refinement.x(); + mapping.subJ = refinedJ % m_refinement.y(); + mapping.subK = refinedK % m_refinement.z(); + } return mapping; } @@ -385,7 +501,7 @@ cvf::Vec3st RigGridExportAdapter::refinement() const //-------------------------------------------------------------------------------------------------- bool RigGridExportAdapter::hasRefinement() const { - return m_refinement.x() > 1 || m_refinement.y() > 1 || m_refinement.z() > 1; + return m_refinement.x() > 1 || m_refinement.y() > 1 || m_refinement.z() > 1 || m_nonUniformRefinement.hasNonUniformRefinement(); } //-------------------------------------------------------------------------------------------------- @@ -467,3 +583,57 @@ std::expected RigGridExportAdapter::transformIjkToSectorC ( ijk.z() - min.z() ) * refinement.z() ); } } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const RigNonUniformRefinement& RigGridExportAdapter::nonUniformRefinement() const +{ + return m_nonUniformRefinement; +} + +//-------------------------------------------------------------------------------------------------- +/// Transform IJK coordinates from global grid space to sector-relative space with non-uniform refinement +//-------------------------------------------------------------------------------------------------- +std::expected RigGridExportAdapter::transformIjkToSectorCoordinates( const caf::VecIjk0& ijk, + const caf::VecIjk0& min, + const caf::VecIjk0& max, + const RigNonUniformRefinement& nonUniformRefinement, + bool applyRefinementCentering, + bool isBoxMaxCoordinate ) +{ + RigBoundingBoxIjk sectorBox( min, max ); + if ( !sectorBox.contains( ijk ) ) + { + return std::unexpected( QString( "IJK coordinates (%1) are outside sector bounds [(%2), (%3)]" ) + .arg( QString::fromStdString( ijk.toString() ) ) + .arg( QString::fromStdString( min.toString() ) ) + .arg( QString::fromStdString( max.toString() ) ) ); + } + + size_t sectorI = ijk.x() - min.x(); + size_t sectorJ = ijk.y() - min.y(); + size_t sectorK = ijk.z() - min.z(); + + if ( applyRefinementCentering ) + { + return caf::VecIjk0( nonUniformRefinement.cumulativeOffset( RigNonUniformRefinement::DimI, sectorI ) + + nonUniformRefinement.subcellCount( RigNonUniformRefinement::DimI, sectorI ) / 2, + nonUniformRefinement.cumulativeOffset( RigNonUniformRefinement::DimJ, sectorJ ) + + nonUniformRefinement.subcellCount( RigNonUniformRefinement::DimJ, sectorJ ) / 2, + nonUniformRefinement.cumulativeOffset( RigNonUniformRefinement::DimK, sectorK ) + + nonUniformRefinement.subcellCount( RigNonUniformRefinement::DimK, sectorK ) / 2 ); + } + else if ( isBoxMaxCoordinate ) + { + return caf::VecIjk0( nonUniformRefinement.cumulativeOffset( RigNonUniformRefinement::DimI, sectorI + 1 ) - 1, + nonUniformRefinement.cumulativeOffset( RigNonUniformRefinement::DimJ, sectorJ + 1 ) - 1, + nonUniformRefinement.cumulativeOffset( RigNonUniformRefinement::DimK, sectorK + 1 ) - 1 ); + } + else + { + return caf::VecIjk0( nonUniformRefinement.cumulativeOffset( RigNonUniformRefinement::DimI, sectorI ), + nonUniformRefinement.cumulativeOffset( RigNonUniformRefinement::DimJ, sectorJ ), + nonUniformRefinement.cumulativeOffset( RigNonUniformRefinement::DimK, sectorK ) ); + } +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigGridExportAdapter.h b/ApplicationLibCode/ReservoirDataModel/RigGridExportAdapter.h index 448f7fb804..1d65964234 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGridExportAdapter.h +++ b/ApplicationLibCode/ReservoirDataModel/RigGridExportAdapter.h @@ -18,6 +18,8 @@ #pragma once +#include "RigNonUniformRefinement.h" + #include "cafVecIjk.h" #include "cvfArray.h" @@ -66,6 +68,12 @@ public: const cvf::Vec3st& refinement, const cvf::UByteArray* cellVisibilityOverrideForActnum = nullptr ); + RigGridExportAdapter( RigEclipseCaseData* eclipseCase, + const cvf::Vec3st& min, + const cvf::Vec3st& max, + const RigNonUniformRefinement& nonUniformRefinement, + const cvf::UByteArray* cellVisibilityOverrideForActnum = nullptr ); + // Grid dimensions (after refinement) size_t cellCountI() const; size_t cellCountJ() const; @@ -88,6 +96,9 @@ public: cvf::Vec3st refinement() const; bool hasRefinement() const; + // Non-uniform refinement accessor + const RigNonUniformRefinement& nonUniformRefinement() const; + // Coordinate transformation utilities static std::expected transformIjkToSectorCoordinates( const caf::VecIjk0& originalIjk, const caf::VecIjk0& min, @@ -96,6 +107,13 @@ public: bool applyRefinementCentering = false, bool isBoxMaxCoordinate = false ); + static std::expected transformIjkToSectorCoordinates( const caf::VecIjk0& originalIjk, + const caf::VecIjk0& min, + const caf::VecIjk0& max, + const RigNonUniformRefinement& nonUniformRefinement, + bool applyRefinementCentering = false, + bool isBoxMaxCoordinate = false ); + private: // Internal methods to handle original vs refined cell access std::array getOriginalCellCorners( size_t origI, size_t origJ, size_t origK ) const; @@ -119,6 +137,8 @@ private: cvf::Vec3st m_max; cvf::Vec3st m_refinement; + RigNonUniformRefinement m_nonUniformRefinement; + // Cached refined grid dimensions size_t m_refinedNI; size_t m_refinedNJ; diff --git a/ApplicationLibCode/ReservoirDataModel/RigNonUniformRefinement.cpp b/ApplicationLibCode/ReservoirDataModel/RigNonUniformRefinement.cpp new file mode 100644 index 0000000000..485272996f --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigNonUniformRefinement.cpp @@ -0,0 +1,248 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2026- Equinor ASA +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "RigNonUniformRefinement.h" + +#include +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigNonUniformRefinement::RigNonUniformRefinement() + : m_sectorSize( 0, 0, 0 ) +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigNonUniformRefinement::RigNonUniformRefinement( const cvf::Vec3st& sectorSize ) + : m_sectorSize( sectorSize ) +{ + // Initialize each cell with default {1.0} (no refinement) + for ( auto dim : { DimI, DimJ, DimK } ) + { + m_fractions[dim].resize( m_sectorSize[dim], std::vector{ 1.0 } ); + rebuildOffsets( dim ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigNonUniformRefinement::setCumulativeFractions( Dimension dim, size_t origIndex, const std::vector& cumulativeFractions ) +{ + if ( origIndex >= m_sectorSize[dim] ) return; + + m_fractions[dim][origIndex] = cumulativeFractions; + rebuildOffsets( dim ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigNonUniformRefinement::widthsToCumulativeFractions( const std::vector& widths ) +{ + if ( widths.empty() ) return { 1.0 }; + + double sum = std::accumulate( widths.begin(), widths.end(), 0.0 ); + if ( sum <= 0.0 ) return { 1.0 }; + + std::vector cumulative; + cumulative.reserve( widths.size() ); + + double runningSum = 0.0; + for ( size_t i = 0; i < widths.size(); ++i ) + { + runningSum += widths[i] / sum; + cumulative.push_back( runningSum ); + } + + // Ensure the last value is exactly 1.0 + cumulative.back() = 1.0; + + return cumulative; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigNonUniformRefinement::subcellCount( Dimension dim, size_t origIndex ) const +{ + if ( origIndex >= m_sectorSize[dim] ) return 1; + return m_fractions[dim][origIndex].size(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigNonUniformRefinement::cumulativeOffset( Dimension dim, size_t origIndex ) const +{ + // Allow origIndex == m_sectorSize[dim] to get the total count (end sentinel) + if ( origIndex > m_sectorSize[dim] ) return 0; + return m_offsets[dim][origIndex]; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigNonUniformRefinement::totalRefinedCount( Dimension dim ) const +{ + if ( m_sectorSize[dim] == 0 ) return 0; + return m_offsets[dim][m_sectorSize[dim]]; +} + +//-------------------------------------------------------------------------------------------------- +/// Maps a refined index back to {originalIndex, subIndex} using binary search on the offset table +//-------------------------------------------------------------------------------------------------- +std::pair RigNonUniformRefinement::mapRefinedToOriginal( Dimension dim, size_t refinedIndex ) const +{ + if ( m_sectorSize[dim] == 0 ) return { 0, 0 }; + + // Binary search: find the largest origIndex such that m_offsets[dim][origIndex] <= refinedIndex + // m_offsets[dim] has size m_sectorSize[dim]+1 + // We want upper_bound - 1 + auto it = std::upper_bound( m_offsets[dim].begin(), m_offsets[dim].end(), refinedIndex ); + size_t origIndex = static_cast( std::distance( m_offsets[dim].begin(), it ) ) - 1; + + size_t subIndex = refinedIndex - m_offsets[dim][origIndex]; + return { origIndex, subIndex }; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const std::vector& RigNonUniformRefinement::cumulativeFractions( Dimension dim, size_t origIndex ) const +{ + static const std::vector defaultFraction = { 1.0 }; + if ( origIndex >= m_sectorSize[dim] ) return defaultFraction; + return m_fractions[dim][origIndex]; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigNonUniformRefinement::sectorSize( Dimension dim ) const +{ + return m_sectorSize[dim]; +} + +//-------------------------------------------------------------------------------------------------- +/// Distribute fractional widths across a range of cells. +/// +/// The widths define subdivision boundaries in the combined range of cells [startIndex, endIndex]. +/// Each original cell occupies 1/numCells of the global range (in index space). +/// Global subdivision boundaries that fall within a cell are converted to cell-local fractions. +/// +/// Example: 3 cells, widths [0.2, 0.3, 0.5] -> cumulative [0.2, 0.5, 1.0] +/// Cell 0 [0, 0.333]: boundary 0.2 falls here -> local fracs [0.6, 1.0] -> 2 subcells +/// Cell 1 [0.333, 0.667]: boundary 0.5 falls here -> local fracs [0.5, 1.0] -> 2 subcells +/// Cell 2 [0.667, 1.0]: no interior boundaries -> local fracs [1.0] -> 1 subcell +//-------------------------------------------------------------------------------------------------- +void RigNonUniformRefinement::distributeWidthsAcrossCells( Dimension dim, size_t startIndex, size_t endIndex, const std::vector& widths ) +{ + if ( startIndex > endIndex || endIndex >= m_sectorSize[dim] || widths.empty() ) return; + + // Convert widths to global cumulative fractions + auto globalFractions = widthsToCumulativeFractions( widths ); + + size_t numCells = endIndex - startIndex + 1; + + for ( size_t c = 0; c < numCells; ++c ) + { + double cellStart = static_cast( c ) / static_cast( numCells ); + double cellEnd = static_cast( c + 1 ) / static_cast( numCells ); + double cellWidth = cellEnd - cellStart; + + std::vector localFractions; + + // Find global boundaries that fall strictly within this cell + for ( double gf : globalFractions ) + { + if ( gf > cellStart && gf < cellEnd ) + { + // Convert to cell-local fraction + localFractions.push_back( ( gf - cellStart ) / cellWidth ); + } + } + + // Always add 1.0 as the final boundary (the cell end) + localFractions.push_back( 1.0 ); + + setCumulativeFractions( dim, startIndex + c, localFractions ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RigNonUniformRefinement::hasNonUniformRefinement() const +{ + for ( auto dim : { DimI, DimJ, DimK } ) + { + if ( totalRefinedCount( dim ) != m_sectorSize[dim] ) return true; + } + return false; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigNonUniformRefinement RigNonUniformRefinement::fromUniform( const cvf::Vec3st& refinement, const cvf::Vec3st& sectorSize ) +{ + RigNonUniformRefinement result( sectorSize ); + + size_t refCount[3] = { refinement.x(), refinement.y(), refinement.z() }; + + for ( auto dim : { DimI, DimJ, DimK } ) + { + if ( refCount[dim] <= 1 ) continue; + + // Create uniform cumulative fractions + std::vector uniformFractions; + uniformFractions.reserve( refCount[dim] ); + for ( size_t i = 1; i <= refCount[dim]; ++i ) + { + uniformFractions.push_back( static_cast( i ) / static_cast( refCount[dim] ) ); + } + uniformFractions.back() = 1.0; + + for ( size_t idx = 0; idx < result.m_sectorSize[dim]; ++idx ) + { + result.m_fractions[dim][idx] = uniformFractions; + } + result.rebuildOffsets( dim ); + } + + return result; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigNonUniformRefinement::rebuildOffsets( Dimension dim ) +{ + m_offsets[dim].resize( m_sectorSize[dim] + 1 ); + m_offsets[dim][0] = 0; + + for ( size_t i = 0; i < m_sectorSize[dim]; ++i ) + { + m_offsets[dim][i + 1] = m_offsets[dim][i] + m_fractions[dim][i].size(); + } +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigNonUniformRefinement.h b/ApplicationLibCode/ReservoirDataModel/RigNonUniformRefinement.h new file mode 100644 index 0000000000..9331ad1ce1 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigNonUniformRefinement.h @@ -0,0 +1,93 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2026- Equinor ASA +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "cvfVector3.h" + +#include +#include +#include + +//================================================================================================== +// +// Non-uniform per-dimension refinement specification for sector model export. +// +// Since refinement "crosses over" the model (all cells at a given I index share the same +// I-refinement), each dimension is independent. +// +// For each dimension: a vector-of-vectors indexed by original cell index (relative to sector), +// containing cumulative fractions in (0, 1]. Cells with no explicit refinement default to {1.0} +// (1 subcell = no refinement). +// +// A cumulative offset table per dimension enables O(1) original->refined lookup +// and O(log n) reverse lookup via binary search. +// +//================================================================================================== +class RigNonUniformRefinement +{ +public: + enum Dimension : size_t + { + DimI = 0, + DimJ = 1, + DimK = 2 + }; + + RigNonUniformRefinement(); + explicit RigNonUniformRefinement( const cvf::Vec3st& sectorSize ); + + // Set cumulative fractions for a specific cell index in a dimension. + // The fractions must be strictly increasing and end at 1.0. + void setCumulativeFractions( Dimension dim, size_t origIndex, const std::vector& cumulativeFractions ); + + // Set fractional widths for a specific cell index in a dimension. + // Widths are normalized so they sum to 1.0, then converted to cumulative fractions. + static std::vector widthsToCumulativeFractions( const std::vector& widths ); + + // Query methods + size_t subcellCount( Dimension dim, size_t origIndex ) const; + size_t cumulativeOffset( Dimension dim, size_t origIndex ) const; + size_t totalRefinedCount( Dimension dim ) const; + std::pair mapRefinedToOriginal( Dimension dim, size_t refinedIndex ) const; + const std::vector& cumulativeFractions( Dimension dim, size_t origIndex ) const; + size_t sectorSize( Dimension dim ) const; + + bool hasNonUniformRefinement() const; + + // Distribute fractional widths across a range of cells in a given dimension. + // The widths span the combined range [startIndex, endIndex] — subdivision boundaries + // are placed in the global fraction space and then assigned to individual cells. + // Each cell in the range occupies an equal share (1/numCells) of the global range. + void distributeWidthsAcrossCells( Dimension dim, size_t startIndex, size_t endIndex, const std::vector& widths ); + + // Create a non-uniform refinement that is equivalent to uniform refinement + static RigNonUniformRefinement fromUniform( const cvf::Vec3st& refinement, const cvf::Vec3st& sectorSize ); + +private: + void rebuildOffsets( Dimension dim ); + + // Per-dimension data: m_fractions[dim][origIndex] = cumulative fractions + std::vector> m_fractions[3]; + + // Per-dimension cumulative offset table: m_offsets[dim][n] = total refined cells for cells 0..n-1 + std::vector m_offsets[3]; + + // Sector dimensions (number of original cells per dimension) + cvf::Vec3st m_sectorSize; +}; diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputSettings.cpp b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputSettings.cpp index 8e32b035c1..8dffdb589d 100644 --- a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputSettings.cpp +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputSettings.cpp @@ -93,6 +93,30 @@ void RigSimulationInputSettings::setRefinement( const cvf::Vec3st& refinement ) m_refinement = refinement; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const RigNonUniformRefinement& RigSimulationInputSettings::nonUniformRefinement() const +{ + return m_nonUniformRefinement; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigSimulationInputSettings::setNonUniformRefinement( const RigNonUniformRefinement& refinement ) +{ + m_nonUniformRefinement = refinement; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RigSimulationInputSettings::hasNonUniformRefinement() const +{ + return m_nonUniformRefinement.hasNonUniformRefinement(); +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputSettings.h b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputSettings.h index 655b6d72c7..819477488c 100644 --- a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputSettings.h +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputSettings.h @@ -20,6 +20,7 @@ #include "RiaModelExportDefines.h" #include "RigModelPaddingSettings.h" +#include "RigNonUniformRefinement.h" #include "cafAppEnum.h" #include "cafVecIjk.h" @@ -60,6 +61,11 @@ public: cvf::Vec3st refinement() const; void setRefinement( const cvf::Vec3st& refinement ); + // Non-uniform refinement + const RigNonUniformRefinement& nonUniformRefinement() const; + void setNonUniformRefinement( const RigNonUniformRefinement& refinement ); + bool hasNonUniformRefinement() const; + // Keyword handling const std::vector& keywordsToRemove() const; void setKeywordsToRemove( const std::vector& keywords ); @@ -89,6 +95,7 @@ private: caf::VecIjk0 m_min; caf::VecIjk0 m_max; cvf::Vec3st m_refinement; + RigNonUniformRefinement m_nonUniformRefinement; std::vector m_keywordsToRemove; std::vector m_bcpropKeywords; BoundaryConditionEnum m_boundaryCondition; diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.cpp b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.cpp index 621eed8ad3..46041634ed 100644 --- a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.cpp +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.cpp @@ -34,6 +34,7 @@ #include "RigGridExportAdapter.h" #include "RigMainGrid.h" #include "RigModelPaddingSettings.h" +#include "RigNonUniformRefinement.h" #include "RigPadModel.h" #include "RigResdataGridConverter.h" #include "RigSimulationInputSettings.h" @@ -193,8 +194,16 @@ std::expected RigSimulationInputTool::updateCornerPointGridInDeck RifOpmFlowDeckFile& deckFile ) { // Get grid bounds for extraction + auto createGridAdapter = [&]() + { + if ( settings.hasNonUniformRefinement() ) + { + return RigGridExportAdapter( eclipseCase->eclipseCaseData(), settings.min(), settings.max(), settings.nonUniformRefinement() ); + } + return RigGridExportAdapter( eclipseCase->eclipseCaseData(), settings.min(), settings.max(), settings.refinement() ); + }; - RigGridExportAdapter gridAdapter( eclipseCase->eclipseCaseData(), settings.min(), settings.max(), settings.refinement() ); + RigGridExportAdapter gridAdapter = createGridAdapter(); std::vector coordArray; std::vector zcornArray; @@ -293,11 +302,16 @@ std::expected RigSimulationInputTool::replaceKeywordValuesInDeckF } // Try to extract keyword data - auto result = RifEclipseInputFileTools::extractKeywordData( eclipseCase->eclipseCaseData(), - keyword, - settings.min(), - settings.max(), - settings.refinement() ); + auto result = settings.hasNonUniformRefinement() ? RifEclipseInputFileTools::extractKeywordData( eclipseCase->eclipseCaseData(), + keyword, + settings.min(), + settings.max(), + settings.nonUniformRefinement() ) + : RifEclipseInputFileTools::extractKeywordData( eclipseCase->eclipseCaseData(), + keyword, + settings.min(), + settings.max(), + settings.refinement() ); if ( result ) { // Replace keyword values in deck with extracted data @@ -337,6 +351,14 @@ std::set RigSimulationInputTool::cropDataKeywordsInDeckFile( RimEcl const auto max = settings.max(); const auto refinement = settings.refinement(); + // Build non-uniform refinement for cropping (handles both uniform and non-uniform cases) + const bool hasNonUniform = settings.hasNonUniformRefinement(); + const auto nonUniformRef = + hasNonUniform + ? settings.nonUniformRefinement() + : RigNonUniformRefinement::fromUniform( refinement, + cvf::Vec3st( max.x() - min.x() + 1, max.y() - min.y() + 1, max.z() - min.z() + 1 ) ); + auto keywords = deckFile.keywords( false ); for ( const auto& name : keywords ) @@ -357,20 +379,35 @@ std::set RigSimulationInputTool::cropDataKeywordsInDeckFile( RimEcl auto cropAndReplace = [&]( const std::vector& fullData ) { std::vector croppedData; - cvf::Vec3st refinedMin( min.x() * refinement.x(), min.y() * refinement.y(), min.z() * refinement.z() ); - cvf::Vec3st refinedMax( ( max.x() + 1 ) * refinement.x(), ( max.y() + 1 ) * refinement.y(), ( max.z() + 1 ) * refinement.z() ); - for ( size_t k = refinedMin.z(); k < refinedMax.z(); ++k ) + size_t sectorNk = max.z() - min.z() + 1; + size_t sectorNj = max.y() - min.y() + 1; + size_t sectorNi = max.x() - min.x() + 1; + + for ( size_t origK = 0; origK < sectorNk; ++origK ) { - size_t mainK = k / refinement.z(); - for ( size_t j = refinedMin.y(); j < refinedMax.y(); ++j ) + size_t mainK = min.z() + origK; + size_t subCountK = nonUniformRef.subcellCount( RigNonUniformRefinement::DimK, origK ); + for ( size_t rk = 0; rk < subCountK; ++rk ) { - size_t mainJ = j / refinement.y(); - for ( size_t i = refinedMin.x(); i < refinedMax.x(); ++i ) + for ( size_t origJ = 0; origJ < sectorNj; ++origJ ) { - size_t mainI = i / refinement.x(); - size_t sourceIdx = mainI + mainJ * fullNx + mainK * fullNx * fullNy; - croppedData.push_back( fullData[sourceIdx] ); + size_t mainJ = min.y() + origJ; + size_t subCountJ = nonUniformRef.subcellCount( RigNonUniformRefinement::DimJ, origJ ); + for ( size_t rj = 0; rj < subCountJ; ++rj ) + { + for ( size_t origI = 0; origI < sectorNi; ++origI ) + { + size_t mainI = min.x() + origI; + size_t sourceIdx = mainI + mainJ * fullNx + mainK * fullNx * fullNy; + size_t subCountI = nonUniformRef.subcellCount( RigNonUniformRefinement::DimI, origI ); + + for ( size_t ri = 0; ri < subCountI; ++ri ) + { + croppedData.push_back( fullData[sourceIdx] ); + } + } + } } } } @@ -419,7 +456,19 @@ std::expected RigSimulationInputTool::addBorderBoundaryConditions } // Create grid adapter for refined grid operations - RigGridExportAdapter gridAdapter( eclipseCase->eclipseCaseData(), settings.min(), settings.max(), settings.refinement(), visibility.p() ); + auto createBorderGridAdapter = [&]() + { + if ( settings.hasNonUniformRefinement() ) + { + return RigGridExportAdapter( eclipseCase->eclipseCaseData(), + settings.min(), + settings.max(), + settings.nonUniformRefinement(), + visibility.p() ); + } + return RigGridExportAdapter( eclipseCase->eclipseCaseData(), settings.min(), settings.max(), settings.refinement(), visibility.p() ); + }; + RigGridExportAdapter gridAdapter = createBorderGridAdapter(); // Create refined visibility array matching refined grid dimensions std::vector refinedVisibility = createRefinedVisibility( gridAdapter ); @@ -657,13 +706,14 @@ static std::expected processBoxKeywordEntry( const Opm::D /// Crop a data keyword to the intersection of its BOX context and the sector. /// Returns nullopt if no modification is needed (e.g. keyword has unexpected structure). //-------------------------------------------------------------------------------------------------- -static std::optional cropDataKeywordInBoxContext( const Opm::DeckKeyword& kw, - const Opm::FileDeck::Index& index, - const caf::VecIjk0& boxMin, - const caf::VecIjk0& boxMax, - const caf::VecIjk0& sectorMin, - const caf::VecIjk0& sectorMax, - const cvf::Vec3st& refinement ) +static std::optional cropDataKeywordInBoxContext( const Opm::DeckKeyword& kw, + const Opm::FileDeck::Index& index, + const caf::VecIjk0& boxMin, + const caf::VecIjk0& boxMax, + const caf::VecIjk0& sectorMin, + const caf::VecIjk0& sectorMax, + const cvf::Vec3st& refinement, + const RigNonUniformRefinement* nonUniformRef = nullptr ) { const std::string& name = kw.name(); @@ -691,11 +741,15 @@ static std::optional cropDataKeywordInBoxContext( const Opm::DeckK for ( size_t k = intMin.z(); k <= static_cast( intMax.z() ); ++k ) { - for ( size_t rk = 0; rk < refinement.z(); ++rk ) + size_t subCountK = nonUniformRef ? nonUniformRef->subcellCount( RigNonUniformRefinement::DimK, k - sectorMin.z() ) + : refinement.z(); + for ( size_t rk = 0; rk < subCountK; ++rk ) { for ( size_t j = intMin.y(); j <= static_cast( intMax.y() ); ++j ) { - for ( size_t rj = 0; rj < refinement.y(); ++rj ) + size_t subCountJ = nonUniformRef ? nonUniformRef->subcellCount( RigNonUniformRefinement::DimJ, j - sectorMin.y() ) + : refinement.y(); + for ( size_t rj = 0; rj < subCountJ; ++rj ) { for ( size_t i = intMin.x(); i <= static_cast( intMax.x() ); ++i ) { @@ -704,7 +758,9 @@ static std::optional cropDataKeywordInBoxContext( const Opm::DeckK size_t boxRelK = k - boxMin.z(); size_t srcIdx = boxRelI + boxRelJ * boxNi + boxRelK * boxNi * boxNj; - for ( size_t ri = 0; ri < refinement.x(); ++ri ) + size_t subCountI = nonUniformRef ? nonUniformRef->subcellCount( RigNonUniformRefinement::DimI, i - sectorMin.x() ) + : refinement.x(); + for ( size_t ri = 0; ri < subCountI; ++ri ) { croppedData.push_back( sourceData[srcIdx] ); } @@ -843,6 +899,9 @@ std::expected RigSimulationInputTool::transformKeywordsInDeckFile const auto sectorMax = settings.max(); const auto refinement = settings.refinement(); + // Build non-uniform refinement pointer for functions that support it + const RigNonUniformRefinement* nonUniformRefPtr = settings.hasNonUniformRefinement() ? &settings.nonUniformRefinement() : nullptr; + // Map keyword names to their record processor functions const std::map recordProcessors = { { "EQUALS", processEqualsRecord }, @@ -878,7 +937,7 @@ std::expected RigSimulationInputTool::transformKeywordsInDeckFile } else if ( insideBox && kw.isDataKeyword() ) { - auto result = cropDataKeywordInBoxContext( kw, it, boxMin, boxMax, sectorMin, sectorMax, refinement ); + auto result = cropDataKeywordInBoxContext( kw, it, boxMin, boxMax, sectorMin, sectorMax, refinement, nonUniformRefPtr ); if ( result ) modifications.push_back( std::move( result.value() ) ); } else if ( auto procIt = recordProcessors.find( name ); procIt != recordProcessors.end() ) @@ -1013,7 +1072,28 @@ std::expected RigSimulationInputTool::updateWelldimsKeyword( RimE if ( wellDims.empty() ) return std::unexpected( QString( "Missing WELLDIMS keyword" ) ); // Scale the max connections with the refinement - int newMaxConnections = wellDims[1] * static_cast( settings.refinement().x() * settings.refinement().y() * settings.refinement().z() ); + int refinementProduct = 1; + if ( settings.hasNonUniformRefinement() ) + { + const auto& nuRef = settings.nonUniformRefinement(); + // Use the maximum subcell product across all cell combinations + size_t maxProduct = 1; + for ( auto dim : { RigNonUniformRefinement::DimI, RigNonUniformRefinement::DimJ, RigNonUniformRefinement::DimK } ) + { + size_t maxSubcells = 1; + for ( size_t idx = 0; idx < nuRef.sectorSize( dim ); ++idx ) + { + maxSubcells = std::max( maxSubcells, nuRef.subcellCount( dim, idx ) ); + } + maxProduct *= maxSubcells; + } + refinementProduct = static_cast( maxProduct ); + } + else + { + refinementProduct = static_cast( settings.refinement().x() * settings.refinement().y() * settings.refinement().z() ); + } + int newMaxConnections = wellDims[1] * refinementProduct; deckFile.setWelldims( wellDims[0], newMaxConnections, wellDims[2], wellDims[3] ); @@ -1133,8 +1213,20 @@ std::expected RigSimulationInputTool::processWelspecsR // Transform with clamped coordinates (this will always succeed) // Center the coordinate in the refined cell block - size_t sectorI = ( clampedI - settings.min().x() ) * settings.refinement().x() + ( settings.refinement().x() + 1 ) / 2; - size_t sectorJ = ( clampedJ - settings.min().y() ) * settings.refinement().y() + ( settings.refinement().y() + 1 ) / 2; + size_t sectorI, sectorJ; + if ( settings.hasNonUniformRefinement() ) + { + const auto& nuRef = settings.nonUniformRefinement(); + size_t relI = clampedI - settings.min().x(); + size_t relJ = clampedJ - settings.min().y(); + sectorI = nuRef.cumulativeOffset( RigNonUniformRefinement::DimI, relI ) + nuRef.subcellCount( RigNonUniformRefinement::DimI, relI ) / 2; + sectorJ = nuRef.cumulativeOffset( RigNonUniformRefinement::DimJ, relJ ) + nuRef.subcellCount( RigNonUniformRefinement::DimJ, relJ ) / 2; + } + else + { + sectorI = ( clampedI - settings.min().x() ) * settings.refinement().x() + ( settings.refinement().x() + 1 ) / 2; + sectorJ = ( clampedJ - settings.min().y() ) * settings.refinement().y() + ( settings.refinement().y() + 1 ) / 2; + } if ( origI < static_cast( settings.min().x() ) || origI > static_cast( settings.max().x() ) || origJ < static_cast( settings.min().y() ) || origJ > static_cast( settings.max().y() ) ) @@ -1190,12 +1282,16 @@ std::expected RigSimulationInputTool::processCompdatRe // Transform K1 caf::VecIjk0 origIjkK1( origI, origJ, origK1 ); auto transformResultK1 = - RigGridExportAdapter::transformIjkToSectorCoordinates( origIjkK1, settings.min(), settings.max(), settings.refinement(), true ); + settings.hasNonUniformRefinement() + ? RigGridExportAdapter::transformIjkToSectorCoordinates( origIjkK1, settings.min(), settings.max(), settings.nonUniformRefinement(), true ) + : RigGridExportAdapter::transformIjkToSectorCoordinates( origIjkK1, settings.min(), settings.max(), settings.refinement(), true ); // Transform K2 caf::VecIjk0 origIjkK2( origI, origJ, origK2 ); auto transformResultK2 = - RigGridExportAdapter::transformIjkToSectorCoordinates( origIjkK2, settings.min(), settings.max(), settings.refinement(), true ); + settings.hasNonUniformRefinement() + ? RigGridExportAdapter::transformIjkToSectorCoordinates( origIjkK2, settings.min(), settings.max(), settings.nonUniformRefinement(), true ) + : RigGridExportAdapter::transformIjkToSectorCoordinates( origIjkK2, settings.min(), settings.max(), settings.refinement(), true ); if ( !transformResultK1 ) { @@ -1253,8 +1349,17 @@ std::expected RigSimulationInputTool::processCompsegsR // Transform I, J, K (first three items) caf::VecIjk1 origIjk = extractIjk( record, 0, 1, 2 ); - auto transformResult = - RigGridExportAdapter::transformIjkToSectorCoordinates( origIjk.toZeroBased(), settings.min(), settings.max(), settings.refinement(), true ); + auto transformResult = settings.hasNonUniformRefinement() + ? RigGridExportAdapter::transformIjkToSectorCoordinates( origIjk.toZeroBased(), + settings.min(), + settings.max(), + settings.nonUniformRefinement(), + true ) + : RigGridExportAdapter::transformIjkToSectorCoordinates( origIjk.toZeroBased(), + settings.min(), + settings.max(), + settings.refinement(), + true ); if ( !transformResult ) { @@ -2440,11 +2545,12 @@ std::expected RigSimulationInputTool::exportEditNncKeyword( RimEc QString( "Found %1 internal EDITNNC connections in sector (out of %2 total)" ).arg( sectorConnections.size() ).arg( allConnections.size() ) ); // Step 3: Transform and refine - bool hasRefinement = ( settings.refinement().x() != 1 || settings.refinement().y() != 1 || settings.refinement().z() != 1 ); + bool hasRef = ( settings.refinement().x() != 1 || settings.refinement().y() != 1 || settings.refinement().z() != 1 ) || + settings.hasNonUniformRefinement(); std::vector transformedConnections; - if ( hasRefinement ) + if ( hasRef ) { // Refine to subcells with corresponding positions for ( const auto& conn : sectorConnections ) diff --git a/ApplicationLibCode/UnitTests/CMakeLists.txt b/ApplicationLibCode/UnitTests/CMakeLists.txt index 2f4c5817af..0a5e6d3614 100644 --- a/ApplicationLibCode/UnitTests/CMakeLists.txt +++ b/ApplicationLibCode/UnitTests/CMakeLists.txt @@ -23,6 +23,7 @@ set(SOURCE_UNITTEST_FILES ${CMAKE_CURRENT_LIST_DIR}/RigReservoir-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigResdataGridConverter-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigGridExportAdapter-Test.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigNonUniformRefinement-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigStatisticsMath-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigWellPathIntersectionTools-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RimWellLogExtractionCurveImpl-Test.cpp diff --git a/ApplicationLibCode/UnitTests/RigNonUniformRefinement-Test.cpp b/ApplicationLibCode/UnitTests/RigNonUniformRefinement-Test.cpp new file mode 100644 index 0000000000..2cff83304d --- /dev/null +++ b/ApplicationLibCode/UnitTests/RigNonUniformRefinement-Test.cpp @@ -0,0 +1,394 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2026- Equinor ASA +// +// ResInsight is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. +// +// See the GNU General Public License at +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "gtest/gtest.h" + +#include "RigNonUniformRefinement.h" + +#include "RiaCellDividingTools.h" +#include "RigGridExportAdapter.h" + +#include "cafVecIjk.h" + +//-------------------------------------------------------------------------------------------------- +/// Test fromUniform() matches expected uniform counts +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, FromUniform ) +{ + auto nuRef = RigNonUniformRefinement::fromUniform( cvf::Vec3st( 3, 2, 1 ), cvf::Vec3st( 4, 5, 6 ) ); + + // Total refined count should match uniform: sectorSize * refinement + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 4u * 3 ); + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimJ ), 5u * 2 ); + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimK ), 6u * 1 ); + + // Each cell in dim 0 should have 3 subcells + for ( size_t i = 0; i < 4; ++i ) + { + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, i ), 3u ); + } + + // Each cell in dim 1 should have 2 subcells + for ( size_t j = 0; j < 5; ++j ) + { + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimJ, j ), 2u ); + } + + // Each cell in dim 2 should have 1 subcell + for ( size_t k = 0; k < 6; ++k ) + { + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimK, k ), 1u ); + } + + EXPECT_TRUE( nuRef.hasNonUniformRefinement() ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test no refinement case +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, NoRefinement ) +{ + auto nuRef = RigNonUniformRefinement::fromUniform( cvf::Vec3st( 1, 1, 1 ), cvf::Vec3st( 3, 3, 3 ) ); + + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 3u ); + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimJ ), 3u ); + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimK ), 3u ); + EXPECT_FALSE( nuRef.hasNonUniformRefinement() ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test subcellCount and cumulativeOffset for mixed refinement +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, MixedRefinement ) +{ + RigNonUniformRefinement nuRef( cvf::Vec3st( 3, 1, 1 ) ); + + // Cell 0: 2 subcells, Cell 1: 5 subcells, Cell 2: 1 subcell (default) + auto fractions2 = RigNonUniformRefinement::widthsToCumulativeFractions( { 0.5, 0.5 } ); + auto fractions5 = RigNonUniformRefinement::widthsToCumulativeFractions( { 0.2, 0.4, 0.5, 0.6, 0.8 } ); + + nuRef.setCumulativeFractions( RigNonUniformRefinement::DimI, 0, fractions2 ); + nuRef.setCumulativeFractions( RigNonUniformRefinement::DimI, 1, fractions5 ); + + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 0 ), 2u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 1 ), 5u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 2 ), 1u ); + + EXPECT_EQ( nuRef.cumulativeOffset( RigNonUniformRefinement::DimI, 0 ), 0u ); + EXPECT_EQ( nuRef.cumulativeOffset( RigNonUniformRefinement::DimI, 1 ), 2u ); + EXPECT_EQ( nuRef.cumulativeOffset( RigNonUniformRefinement::DimI, 2 ), 7u ); + + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 8u ); // 2 + 5 + 1 + EXPECT_TRUE( nuRef.hasNonUniformRefinement() ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test mapRefinedToOriginal round-trips with cumulativeOffset +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, MapRefinedToOriginalRoundTrip ) +{ + RigNonUniformRefinement nuRef( cvf::Vec3st( 4, 1, 1 ) ); + + // Cell 0: 3 subcells, Cell 1: 1 subcell, Cell 2: 2 subcells, Cell 3: 4 subcells + nuRef.setCumulativeFractions( RigNonUniformRefinement::DimI, 0, RigNonUniformRefinement::widthsToCumulativeFractions( { 1.0, 1.0, 1.0 } ) ); + // Cell 1: default (1 subcell) + nuRef.setCumulativeFractions( RigNonUniformRefinement::DimI, 2, RigNonUniformRefinement::widthsToCumulativeFractions( { 1.0, 1.0 } ) ); + nuRef.setCumulativeFractions( RigNonUniformRefinement::DimI, + 3, + RigNonUniformRefinement::widthsToCumulativeFractions( { 1.0, 1.0, 1.0, 1.0 } ) ); + + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 10u ); // 3 + 1 + 2 + 4 + + // Verify round-trip: for each original cell, cumulativeOffset should map back correctly + for ( size_t origIdx = 0; origIdx < 4; ++origIdx ) + { + size_t offset = nuRef.cumulativeOffset( RigNonUniformRefinement::DimI, origIdx ); + auto [mappedOrig, mappedSub] = nuRef.mapRefinedToOriginal( RigNonUniformRefinement::DimI, offset ); + EXPECT_EQ( mappedOrig, origIdx ); + EXPECT_EQ( mappedSub, 0u ); + } + + // Check specific refined indices + auto [orig0, sub0] = nuRef.mapRefinedToOriginal( RigNonUniformRefinement::DimI, 0 ); // First subcell of cell 0 + EXPECT_EQ( orig0, 0u ); + EXPECT_EQ( sub0, 0u ); + + auto [orig1, sub1] = nuRef.mapRefinedToOriginal( RigNonUniformRefinement::DimI, 2 ); // Last subcell of cell 0 + EXPECT_EQ( orig1, 0u ); + EXPECT_EQ( sub1, 2u ); + + auto [orig2, sub2] = nuRef.mapRefinedToOriginal( RigNonUniformRefinement::DimI, 3 ); // Cell 1 (only subcell) + EXPECT_EQ( orig2, 1u ); + EXPECT_EQ( sub2, 0u ); + + auto [orig3, sub3] = nuRef.mapRefinedToOriginal( RigNonUniformRefinement::DimI, 4 ); // First subcell of cell 2 + EXPECT_EQ( orig3, 2u ); + EXPECT_EQ( sub3, 0u ); + + auto [orig4, sub4] = nuRef.mapRefinedToOriginal( RigNonUniformRefinement::DimI, 9 ); // Last subcell of cell 3 + EXPECT_EQ( orig4, 3u ); + EXPECT_EQ( sub4, 3u ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test widthsToCumulativeFractions normalization +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, WidthsToCumulativeFractions ) +{ + // Equal widths + auto fractions = RigNonUniformRefinement::widthsToCumulativeFractions( { 1.0, 1.0, 1.0 } ); + EXPECT_EQ( fractions.size(), 3u ); + EXPECT_NEAR( fractions[0], 1.0 / 3.0, 1e-10 ); + EXPECT_NEAR( fractions[1], 2.0 / 3.0, 1e-10 ); + EXPECT_DOUBLE_EQ( fractions[2], 1.0 ); + + // Non-equal widths + auto fractions2 = RigNonUniformRefinement::widthsToCumulativeFractions( { 1.0, 3.0 } ); + EXPECT_EQ( fractions2.size(), 2u ); + EXPECT_NEAR( fractions2[0], 0.25, 1e-10 ); + EXPECT_DOUBLE_EQ( fractions2[1], 1.0 ); + + // Empty input + auto fractions3 = RigNonUniformRefinement::widthsToCumulativeFractions( {} ); + EXPECT_EQ( fractions3.size(), 1u ); + EXPECT_DOUBLE_EQ( fractions3[0], 1.0 ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test single cell edge case +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, SingleCell ) +{ + RigNonUniformRefinement nuRef( cvf::Vec3st( 1, 1, 1 ) ); + + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 1u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 0 ), 1u ); + EXPECT_EQ( nuRef.cumulativeOffset( RigNonUniformRefinement::DimI, 0 ), 0u ); + + auto [orig, sub] = nuRef.mapRefinedToOriginal( RigNonUniformRefinement::DimI, 0 ); + EXPECT_EQ( orig, 0u ); + EXPECT_EQ( sub, 0u ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test distributeWidthsAcrossCells: 3 equal widths across 3 cells → 1 subcell each +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, DistributeEqualWidthsAcrossEqualCells ) +{ + RigNonUniformRefinement nuRef( cvf::Vec3st( 3, 1, 1 ) ); + + // 3 equal widths across 3 cells: each boundary lands on a cell boundary + // so each cell gets exactly 1 subcell + nuRef.distributeWidthsAcrossCells( RigNonUniformRefinement::DimI, 0, 2, { 1.0, 1.0, 1.0 } ); + + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 0 ), 1u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 1 ), 1u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 2 ), 1u ); + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 3u ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test distributeWidthsAcrossCells: 6 equal widths across 3 cells → 2 subcells each +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, DistributeDoubleWidthsAcrossCells ) +{ + RigNonUniformRefinement nuRef( cvf::Vec3st( 3, 1, 1 ) ); + + // 6 equal widths across 3 cells: 2 subcell boundaries per cell + nuRef.distributeWidthsAcrossCells( RigNonUniformRefinement::DimI, 0, 2, { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 } ); + + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 0 ), 2u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 1 ), 2u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 2 ), 2u ); + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 6u ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test distributeWidthsAcrossCells: uneven widths distributed across cells +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, DistributeUnevenWidthsAcrossCells ) +{ + RigNonUniformRefinement nuRef( cvf::Vec3st( 3, 1, 1 ) ); + + // 5 widths [0.1, 0.1, 0.1, 0.1, 0.6] across 3 cells + // Cumulative fractions: [0.1, 0.2, 0.3, 0.4, 1.0] + // Cell 0 [0, 0.333]: boundaries 0.1, 0.2, 0.3 fall inside → 4 subcells + // Cell 1 [0.333, 0.667]: boundary 0.4 falls inside → 2 subcells + // Cell 2 [0.667, 1.0]: no interior boundaries → 1 subcell + nuRef.distributeWidthsAcrossCells( RigNonUniformRefinement::DimI, 0, 2, { 0.1, 0.1, 0.1, 0.1, 0.6 } ); + + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 0 ), 4u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 1 ), 2u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 2 ), 1u ); + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 7u ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test distributeWidthsAcrossCells: single cell range +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, DistributeWidthsSingleCell ) +{ + RigNonUniformRefinement nuRef( cvf::Vec3st( 5, 1, 1 ) ); + + // Apply 3 widths to just cell 2 + nuRef.distributeWidthsAcrossCells( RigNonUniformRefinement::DimI, 2, 2, { 0.2, 0.3, 0.5 } ); + + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 0 ), 1u ); // Unaffected + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 1 ), 1u ); // Unaffected + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 2 ), 3u ); // 3 subcells + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 3 ), 1u ); // Unaffected + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 4 ), 1u ); // Unaffected + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 7u ); // 1+1+3+1+1 +} + +//-------------------------------------------------------------------------------------------------- +/// Test distributeWidthsAcrossCells: many widths concentrate in first cell +//-------------------------------------------------------------------------------------------------- +TEST( RigNonUniformRefinement, DistributeConcentratedWidths ) +{ + RigNonUniformRefinement nuRef( cvf::Vec3st( 2, 1, 1 ) ); + + // 4 widths [0.1, 0.1, 0.1, 0.7] across 2 cells + // Cumulative fractions: [0.1, 0.2, 0.3, 1.0] + // Cell 0 [0, 0.5]: boundaries 0.1, 0.2, 0.3 fall inside → 4 subcells + // Cell 1 [0.5, 1.0]: no interior boundaries → 1 subcell + nuRef.distributeWidthsAcrossCells( RigNonUniformRefinement::DimI, 0, 1, { 0.1, 0.1, 0.1, 0.7 } ); + + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 0 ), 4u ); + EXPECT_EQ( nuRef.subcellCount( RigNonUniformRefinement::DimI, 1 ), 1u ); + EXPECT_EQ( nuRef.totalRefinedCount( RigNonUniformRefinement::DimI ), 5u ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test non-uniform geometry: equal fractions should match uniform createHexCornerCoords +//-------------------------------------------------------------------------------------------------- +TEST( RiaCellDividingTools, NonUniformEqualFractionsMatchUniform ) +{ + std::array corners = { cvf::Vec3d( 0, 0, 0 ), + cvf::Vec3d( 1, 0, 0 ), + cvf::Vec3d( 1, 1, 0 ), + cvf::Vec3d( 0, 1, 0 ), + cvf::Vec3d( 0, 0, 1 ), + cvf::Vec3d( 1, 0, 1 ), + cvf::Vec3d( 1, 1, 1 ), + cvf::Vec3d( 0, 1, 1 ) }; + + size_t nx = 3, ny = 2, nz = 2; + + // Uniform result + auto uniformResult = RiaCellDividingTools::createHexCornerCoords( corners, nx, ny, nz ); + + // Non-uniform with equal fractions + std::vector fracX = { 1.0 / 3, 2.0 / 3, 1.0 }; + std::vector fracY = { 0.5, 1.0 }; + std::vector fracZ = { 0.5, 1.0 }; + + auto nonUniformResult = RiaCellDividingTools::createHexCornerCoords( corners, fracX, fracY, fracZ ); + + ASSERT_EQ( uniformResult.size(), nonUniformResult.size() ); + + for ( size_t i = 0; i < uniformResult.size(); ++i ) + { + EXPECT_NEAR( uniformResult[i].x(), nonUniformResult[i].x(), 1e-10 ) << "Mismatch at index " << i << " (x)"; + EXPECT_NEAR( uniformResult[i].y(), nonUniformResult[i].y(), 1e-10 ) << "Mismatch at index " << i << " (y)"; + EXPECT_NEAR( uniformResult[i].z(), nonUniformResult[i].z(), 1e-10 ) << "Mismatch at index " << i << " (z)"; + } +} + +//-------------------------------------------------------------------------------------------------- +/// Test non-uniform fractions produce geometrically correct subcells +//-------------------------------------------------------------------------------------------------- +TEST( RiaCellDividingTools, NonUniformFractionsGeometry ) +{ + std::array corners = { cvf::Vec3d( 0, 0, 0 ), + cvf::Vec3d( 10, 0, 0 ), + cvf::Vec3d( 10, 10, 0 ), + cvf::Vec3d( 0, 10, 0 ), + cvf::Vec3d( 0, 0, 10 ), + cvf::Vec3d( 10, 0, 10 ), + cvf::Vec3d( 10, 10, 10 ), + cvf::Vec3d( 0, 10, 10 ) }; + + // Non-uniform X: 20% and 80% + std::vector fracX = { 0.2, 1.0 }; + std::vector fracY = { 1.0 }; // No refinement + std::vector fracZ = { 1.0 }; // No refinement + + auto result = RiaCellDividingTools::createHexCornerCoords( corners, fracX, fracY, fracZ ); + + // Should produce 2*1*1 = 2 subcells, each with 8 corners + ASSERT_EQ( result.size(), 16u ); + + // First subcell: x range [0, 2] + EXPECT_NEAR( result[0].x(), 0.0, 1e-10 ); // corner 0 x + EXPECT_NEAR( result[1].x(), 2.0, 1e-10 ); // corner 1 x + + // Second subcell: x range [2, 10] + EXPECT_NEAR( result[8].x(), 2.0, 1e-10 ); // corner 0 x + EXPECT_NEAR( result[9].x(), 10.0, 1e-10 ); // corner 1 x +} + +//-------------------------------------------------------------------------------------------------- +/// Test transformIjkToSectorCoordinates with non-uniform refinement +//-------------------------------------------------------------------------------------------------- +TEST( RigGridExportAdapter, TransformIjkNonUniform ) +{ + RigNonUniformRefinement nuRef( cvf::Vec3st( 3, 2, 1 ) ); + + // I: cell 0 -> 2 subcells, cell 1 -> 3 subcells, cell 2 -> 1 subcell + nuRef.setCumulativeFractions( RigNonUniformRefinement::DimI, 0, RigNonUniformRefinement::widthsToCumulativeFractions( { 1.0, 1.0 } ) ); + nuRef.setCumulativeFractions( RigNonUniformRefinement::DimI, 1, RigNonUniformRefinement::widthsToCumulativeFractions( { 1.0, 1.0, 1.0 } ) ); + // Cell 2: default (1 subcell) + + // J: cell 0 -> 1 subcell, cell 1 -> 2 subcells + nuRef.setCumulativeFractions( RigNonUniformRefinement::DimJ, 1, RigNonUniformRefinement::widthsToCumulativeFractions( { 1.0, 1.0 } ) ); + + caf::VecIjk0 sectorMin( 10, 20, 30 ); + caf::VecIjk0 sectorMax( 12, 21, 30 ); // 3 cells I, 2 cells J, 1 cell K + + // Test min coordinate (default): cell (10, 20, 30) -> offset(0), offset(0), offset(0) = (0, 0, 0) + auto result = RigGridExportAdapter::transformIjkToSectorCoordinates( caf::VecIjk0( 10, 20, 30 ), sectorMin, sectorMax, nuRef ); + ASSERT_TRUE( result.has_value() ); + EXPECT_EQ( result->x(), 0u ); + EXPECT_EQ( result->y(), 0u ); + EXPECT_EQ( result->z(), 0u ); + + // Test cell (11, 20, 30) -> offset(1) = 2, offset(0) = 0, offset(0) = 0 + auto result2 = RigGridExportAdapter::transformIjkToSectorCoordinates( caf::VecIjk0( 11, 20, 30 ), sectorMin, sectorMax, nuRef ); + ASSERT_TRUE( result2.has_value() ); + EXPECT_EQ( result2->x(), 2u ); + EXPECT_EQ( result2->y(), 0u ); + + // Test max coordinate: cell (12, 21, 30) -> offset(3)-1 = 5, offset(2)-1 = 2, offset(1)-1 = 0 + auto result3 = + RigGridExportAdapter::transformIjkToSectorCoordinates( caf::VecIjk0( 12, 21, 30 ), sectorMin, sectorMax, nuRef, false, true ); + ASSERT_TRUE( result3.has_value() ); + EXPECT_EQ( result3->x(), 5u ); // 2 + 3 + 1 - 1 = 5 + EXPECT_EQ( result3->y(), 2u ); // 1 + 2 - 1 = 2 + EXPECT_EQ( result3->z(), 0u ); + + // Test centering: cell (11, 20, 30) -> offset(1) + subcellCount(1)/2 = 2 + 1 = 3 + auto result4 = RigGridExportAdapter::transformIjkToSectorCoordinates( caf::VecIjk0( 11, 20, 30 ), sectorMin, sectorMax, nuRef, true ); + ASSERT_TRUE( result4.has_value() ); + EXPECT_EQ( result4->x(), 3u ); // offset(1)=2, subcellCount(1)=3, center=2+3/2=3 + EXPECT_EQ( result4->y(), 0u ); // offset(0)=0, subcellCount(0)=1, center=0+1/2=0 + + // Test out of bounds + auto result5 = RigGridExportAdapter::transformIjkToSectorCoordinates( caf::VecIjk0( 9, 20, 30 ), sectorMin, sectorMax, nuRef ); + EXPECT_FALSE( result5.has_value() ); +}