diff --git a/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.cpp b/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.cpp index 8821377cd1..e4f5e1b2ea 100644 --- a/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.cpp +++ b/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.cpp @@ -17,19 +17,16 @@ ///////////////////////////////////////////////////////////////////////////////// #include "RimGeoMechContourMapProjection.h" -#include "RiaImageTools.h" -#include "RiaWeightedMeanCalculator.h" - -#include "RigCellGeometryTools.h" #include "RigContourMapCalculator.h" #include "RigContourMapGrid.h" +#include "RigContourMapProjection.h" #include "RigFemAddressDefines.h" #include "RigFemPart.h" #include "RigFemPartCollection.h" #include "RigFemPartGrid.h" #include "RigFemPartResultsCollection.h" #include "RigGeoMechCaseData.h" -#include "RigHexIntersectionTools.h" +#include "RigGeoMechContourMapProjection.h" #include "RimCellFilterCollection.h" #include "RimGeoMechCase.h" @@ -46,16 +43,12 @@ #include "cvfStructGridGeometryGenerator.h" #include "cvfVector3.h" -#include -#include - CAF_PDM_SOURCE_INIT( RimGeoMechContourMapProjection, "RimGeoMechContourMapProjection" ); //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RimGeoMechContourMapProjection::RimGeoMechContourMapProjection() - : m_kLayers( 0u ) { CAF_PDM_InitObject( "RimContourMapProjection", ":/2DMapProjection16x16.png" ); CAF_PDM_InitField( &m_limitToPorePressureRegions, "LimitToPorRegion", true, "Limit to Pore Pressure regions" ); @@ -97,8 +90,8 @@ void RimGeoMechContourMapProjection::updateLegend() { RimGeoMechCellColors* cellColors = view()->cellResult(); - double minVal = minValue( m_aggregatedResults ); - double maxVal = maxValue( m_aggregatedResults ); + double minVal = m_contourMapProjection->minValue(); + double maxVal = m_contourMapProjection->maxValue(); std::pair minmaxValAllTimeSteps = minmaxValuesAllTimeSteps(); @@ -174,61 +167,16 @@ cvf::ref RimGeoMechContourMapProjection::getCellVisibility() co return cellGridIdxVisibility; } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -cvf::BoundingBox RimGeoMechContourMapProjection::calculateExpandedPorBarBBox( int timeStep, int frameIndex ) const -{ - RigFemResultAddress porBarAddr( RigFemResultPosEnum::RIG_ELEMENT_NODAL, - RigFemAddressDefines::porBar(), - view()->cellResult()->resultComponentName().toStdString() ); - RigGeoMechCaseData* caseData = geoMechCase()->geoMechData(); - RigFemPartResultsCollection* resultCollection = caseData->femPartResults(); - - const std::vector& resultValues = resultCollection->resultValues( porBarAddr, 0, timeStep, frameIndex ); - cvf::BoundingBox boundingBox; - - if ( resultValues.empty() ) - { - return boundingBox; - } - - for ( int i = 0; i < m_femPart->elementCount(); ++i ) - { - size_t resValueIdx = m_femPart->elementNodeResultIdx( (int)i, 0 ); - CVF_ASSERT( resValueIdx < resultValues.size() ); - double scalarValue = resultValues[resValueIdx]; - bool validPorValue = scalarValue != std::numeric_limits::infinity(); - - if ( validPorValue ) - { - std::array hexCorners; - m_femPartGrid->cellCornerVertices( i, hexCorners.data() ); - for ( size_t c = 0; c < 8; ++c ) - { - boundingBox.add( hexCorners[c] ); - } - } - } - cvf::Vec3d boxMin = boundingBox.min(); - cvf::Vec3d boxMax = boundingBox.max(); - cvf::Vec3d boxExtent = boundingBox.extent(); - boxMin.x() -= boxExtent.x() * 0.5 * m_paddingAroundPorePressureRegion(); - boxMin.y() -= boxExtent.y() * 0.5 * m_paddingAroundPorePressureRegion(); - boxMax.x() += boxExtent.x() * 0.5 * m_paddingAroundPorePressureRegion(); - boxMax.y() += boxExtent.y() * 0.5 * m_paddingAroundPorePressureRegion(); - return cvf::BoundingBox( boxMin, boxMax ); -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimGeoMechContourMapProjection::updateGridInformation() { RimGeoMechCase* geoMechCase = this->geoMechCase(); - m_femPart = geoMechCase->geoMechData()->femParts()->part( 0 ); - m_femPartGrid = m_femPart->getOrCreateStructGrid(); - m_kLayers = m_femPartGrid->cellCountK(); + if ( !geoMechCase ) return; + + m_femPart = geoMechCase->geoMechData()->femParts()->part( 0 ); + m_femPartGrid = m_femPart->getOrCreateStructGrid(); m_femPart->ensureIntersectionSearchTreeIsBuilt(); cvf::BoundingBox gridBoundingBox = geoMechCase->activeCellsBoundingBox(); @@ -238,7 +186,12 @@ void RimGeoMechContourMapProjection::updateGridInformation() { auto [stepIdx, frameIdx] = view()->currentStepAndDataFrame(); - expandedBoundingBox = calculateExpandedPorBarBBox( stepIdx, frameIdx ); + expandedBoundingBox = + RigGeoMechContourMapProjection::calculateExpandedPorBarBBox( *geoMechCase->geoMechData(), + view()->cellResult()->resultComponentName().toStdString(), + stepIdx, + frameIdx, + m_paddingAroundPorePressureRegion() ); if ( !expandedBoundingBox.isValid() ) { m_limitToPorePressureRegions = false; @@ -259,68 +212,11 @@ void RimGeoMechContourMapProjection::updateGridInformation() } expandedBoundingBox = cvf::BoundingBox( minExpandedPoint, maxExpandedPoint ); - m_contourMapGrid = std::make_unique( gridBoundingBox, expandedBoundingBox, sampleSpacing() ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector RimGeoMechContourMapProjection::getMapCellVisibility() -{ - cvf::Vec2ui nCellsIJ = numberOfElementsIJ(); - std::vector> distanceImage( nCellsIJ.x(), std::vector( nCellsIJ.y(), 0u ) ); - - std::vector mapCellVisibility; - RigFemResultAddress resAddr = view()->cellResult()->resultAddress(); - - if ( m_limitToPorePressureRegions ) - { - resAddr = RigFemAddressDefines::elementNodalPorBarAddress(); - } - - std::vector cellResults = generateResultsFromAddress( resAddr, mapCellVisibility, view()->currentTimeStep() ); - - mapCellVisibility.resize( numberOfCells(), true ); - CVF_ASSERT( mapCellVisibility.size() == cellResults.size() ); - - { - cvf::BoundingBox validResBoundingBox; - for ( size_t cellIndex = 0; cellIndex < cellResults.size(); ++cellIndex ) - { - cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( cellIndex ); - if ( cellResults[cellIndex] != std::numeric_limits::infinity() ) - { - distanceImage[ij.x()][ij.y()] = 1u; - validResBoundingBox.add( cvf::Vec3d( m_contourMapGrid->cellCenterPosition( ij.x(), ij.y() ), 0.0 ) ); - } - else - { - mapCellVisibility[cellIndex] = false; - } - } - - if ( m_limitToPorePressureRegions && m_paddingAroundPorePressureRegion > 0.0 ) - { - RiaImageTools::distanceTransform2d( distanceImage ); - - cvf::Vec3d porExtent = validResBoundingBox.extent(); - double radius = std::max( porExtent.x(), porExtent.y() ) * 0.25; - double expansion = m_paddingAroundPorePressureRegion * radius; - size_t cellPadding = std::ceil( expansion / sampleSpacing() ); - for ( size_t cellIndex = 0; cellIndex < cellResults.size(); ++cellIndex ) - { - if ( !mapCellVisibility[cellIndex] ) - { - cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( cellIndex ); - if ( distanceImage[ij.x()][ij.y()] < cellPadding * cellPadding ) - { - mapCellVisibility[cellIndex] = true; - } - } - } - } - } - return mapCellVisibility; + m_contourMapGrid = std::make_unique( gridBoundingBox, expandedBoundingBox, sampleSpacing() ); + m_contourMapProjection = std::make_unique( *geoMechCase->geoMechData(), + *m_contourMapGrid, + m_limitToPorePressureRegions, + m_paddingAroundPorePressureRegion ); } //-------------------------------------------------------------------------------------------------- @@ -334,14 +230,30 @@ std::vector RimGeoMechContourMapProjection::retrieveParameterWeights() //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::vector RimGeoMechContourMapProjection::generateResults( int viewerStepIndex ) +std::vector RimGeoMechContourMapProjection::generateResults( int viewerStepIndex ) const { - RimGeoMechCellColors* cellColors = view()->cellResult(); - RigFemResultAddress resultAddress = cellColors->resultAddress(); + if ( m_contourMapProjection ) + { + RimGeoMechCellColors* cellColors = view()->cellResult(); + RigFemResultAddress resultAddress = cellColors->resultAddress(); + return generateResultsFromAddress( resultAddress, m_mapCellVisibility, viewerStepIndex ); + } - std::vector aggregatedResults = generateResultsFromAddress( resultAddress, m_mapCellVisibility, viewerStepIndex ); + return {}; +} - return aggregatedResults; +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimGeoMechContourMapProjection::generateAndSaveResults( int timeStep ) +{ + if ( m_contourMapProjection ) + { + RimGeoMechCellColors* cellColors = view()->cellResult(); + RigFemResultAddress resultAddress = cellColors->resultAddress(); + dynamic_cast( m_contourMapProjection.get() ) + ->generateAndSaveResults( resultAddress, m_resultAggregation(), timeStep ); + } } //-------------------------------------------------------------------------------------------------- @@ -349,64 +261,10 @@ std::vector RimGeoMechContourMapProjection::generateResults( int viewerS //-------------------------------------------------------------------------------------------------- std::vector RimGeoMechContourMapProjection::generateResultsFromAddress( RigFemResultAddress resultAddress, const std::vector& mapCellVisibility, - int viewerStepIndex ) + int viewerStepIndex ) const { - RigGeoMechCaseData* caseData = geoMechCase()->geoMechData(); - RigFemPartResultsCollection* resultCollection = caseData->femPartResults(); - size_t nCells = numberOfCells(); - std::vector aggregatedResults = std::vector( nCells, std::numeric_limits::infinity() ); - - auto [stepIdx, frameIdx] = caseData->femPartResults()->stepListIndexToTimeStepAndDataFrameIndex( viewerStepIndex ); - - bool wasInvalid = false; - if ( !resultAddress.isValid() ) - { - wasInvalid = true; - resultAddress = RigFemAddressDefines::elementNodalPorBarAddress(); - } - - if ( resultAddress.fieldName == "PP" ) - { - resultAddress.fieldName = RigFemAddressDefines::porBar(); // More likely to be in memory than POR - } - if ( resultAddress.fieldName == RigFemAddressDefines::porBar() ) - { - resultAddress.resultPosType = RIG_ELEMENT_NODAL; - } - else if ( resultAddress.resultPosType == RIG_FORMATION_NAMES ) - { - resultAddress.resultPosType = RIG_ELEMENT_NODAL; // formation indices are stored per element node result. - } - - std::vector resultValuesF = resultCollection->resultValues( resultAddress, 0, stepIdx, frameIdx ); - if ( resultValuesF.empty() ) return aggregatedResults; - - std::vector resultValues = gridCellValues( resultAddress, resultValuesF ); - - if ( wasInvalid ) - { - // For invalid result addresses we just use the POR-Bar result to get the reservoir region - // And display a dummy 0-result in the region. - for ( double& value : resultValues ) - { - if ( value != std::numeric_limits::infinity() ) - { - value = 0.0; - } - } - } - -#pragma omp parallel for - for ( int index = 0; index < static_cast( nCells ); ++index ) - { - if ( mapCellVisibility.empty() || mapCellVisibility[index] ) - { - cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( index ); - aggregatedResults[index] = calculateValueInMapCell( ij.x(), ij.y(), resultValues ); - } - } - - return aggregatedResults; + return dynamic_cast( m_contourMapProjection.get() ) + ->generateResultsFromAddress( resultAddress, mapCellVisibility, m_resultAggregation(), viewerStepIndex ); } //-------------------------------------------------------------------------------------------------- @@ -436,132 +294,6 @@ RimGridView* RimGeoMechContourMapProjection::baseView() const return view(); } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector RimGeoMechContourMapProjection::findIntersectingCells( const cvf::BoundingBox& bbox ) const -{ - return m_femPart->findIntersectingElementsWithExistingSearchTree( bbox ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -size_t RimGeoMechContourMapProjection::kLayer( size_t globalCellIdx ) const -{ - size_t i, j, k; - m_femPartGrid->ijkFromCellIndex( globalCellIdx, &i, &j, &k ); - return k; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -size_t RimGeoMechContourMapProjection::kLayers() const -{ - return m_kLayers; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimGeoMechContourMapProjection::calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const -{ - std::array hexCorners; - m_femPartGrid->cellCornerVertices( globalCellIdx, hexCorners.data() ); - - cvf::BoundingBox overlapBBox; - std::array overlapCorners; - if ( RigCellGeometryTools::estimateHexOverlapWithBoundingBox( hexCorners, bbox, &overlapCorners, &overlapBBox ) ) - { - double overlapVolume = RigCellGeometryTools::calculateCellVolume( overlapCorners ); - return overlapVolume; - } - return 0.0; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimGeoMechContourMapProjection::calculateRayLengthInCell( size_t globalCellIdx, - const cvf::Vec3d& highestPoint, - const cvf::Vec3d& lowestPoint ) const -{ - std::array hexCorners; - - const std::vector& nodeCoords = m_femPart->nodes().coordinates; - const int* cornerIndices = m_femPart->connectivities( globalCellIdx ); - - hexCorners[0] = cvf::Vec3d( nodeCoords[cornerIndices[0]] ); - hexCorners[1] = cvf::Vec3d( nodeCoords[cornerIndices[1]] ); - hexCorners[2] = cvf::Vec3d( nodeCoords[cornerIndices[2]] ); - hexCorners[3] = cvf::Vec3d( nodeCoords[cornerIndices[3]] ); - hexCorners[4] = cvf::Vec3d( nodeCoords[cornerIndices[4]] ); - hexCorners[5] = cvf::Vec3d( nodeCoords[cornerIndices[5]] ); - hexCorners[6] = cvf::Vec3d( nodeCoords[cornerIndices[6]] ); - hexCorners[7] = cvf::Vec3d( nodeCoords[cornerIndices[7]] ); - - std::vector intersections; - - if ( RigHexIntersectionTools::lineHexCellIntersection( highestPoint, lowestPoint, hexCorners.data(), 0, &intersections ) ) - { - double lengthInCell = ( intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint ).length(); - return lengthInCell; - } - return 0.0; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimGeoMechContourMapProjection::getParameterWeightForCell( size_t globalCellIdx, const std::vector& parameterWeights ) const -{ - if ( parameterWeights.empty() ) return 1.0; - - return parameterWeights[globalCellIdx]; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector RimGeoMechContourMapProjection::gridCellValues( RigFemResultAddress resAddr, std::vector& resultValues ) const -{ - std::vector gridCellValues( m_femPart->elementCount(), std::numeric_limits::infinity() ); - for ( size_t globalCellIdx = 0; globalCellIdx < static_cast( m_femPart->elementCount() ); ++globalCellIdx ) - { - RigElementType elmType = m_femPart->elementType( globalCellIdx ); - if ( !RigFemTypes::is8NodeElement( elmType ) ) continue; - - if ( resAddr.resultPosType == RIG_ELEMENT ) - { - gridCellValues[globalCellIdx] = static_cast( resultValues[globalCellIdx] ); - } - else if ( resAddr.resultPosType == RIG_ELEMENT_NODAL ) - { - RiaWeightedMeanCalculator cellAverage; - for ( int i = 0; i < 8; ++i ) - { - size_t gridResultValueIdx = - m_femPart->resultValueIdxFromResultPosType( resAddr.resultPosType, static_cast( globalCellIdx ), i ); - cellAverage.addValueAndWeight( resultValues[gridResultValueIdx], 1.0 ); - } - - gridCellValues[globalCellIdx] = static_cast( cellAverage.weightedMean() ); - } - else - { - RiaWeightedMeanCalculator cellAverage; - const int* elmNodeIndices = m_femPart->connectivities( globalCellIdx ); - for ( int i = 0; i < 8; ++i ) - { - cellAverage.addValueAndWeight( resultValues[elmNodeIndices[i]], 1.0 ); - } - gridCellValues[globalCellIdx] = static_cast( cellAverage.weightedMean() ); - } - } - return gridCellValues; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -672,20 +404,14 @@ std::pair RimGeoMechContourMapProjection::minmaxValuesAllTimeSte { clearTimeStepRange(); - m_minResultAllTimeSteps = std::min( m_minResultAllTimeSteps, minValue( m_aggregatedResults ) ); - m_maxResultAllTimeSteps = std::max( m_maxResultAllTimeSteps, maxValue( m_aggregatedResults ) ); - - if ( geoMechCase() && geoMechCase()->geoMechData() && geoMechCase()->geoMechData()->femPartResults() ) + if ( geoMechCase()->geoMechData()->femPartResults() ) { int steps = geoMechCase()->geoMechData()->femPartResults()->totalSteps(); - for ( int stepIdx = 0; stepIdx < steps; stepIdx++ ) { - if ( stepIdx == m_currentResultTimestep ) continue; - std::vector aggregatedResults = generateResults( stepIdx ); - m_minResultAllTimeSteps = std::min( m_minResultAllTimeSteps, minValue( aggregatedResults ) ); - m_maxResultAllTimeSteps = std::max( m_maxResultAllTimeSteps, maxValue( aggregatedResults ) ); + m_minResultAllTimeSteps = std::min( m_minResultAllTimeSteps, RigContourMapProjection::minValue( aggregatedResults ) ); + m_maxResultAllTimeSteps = std::max( m_maxResultAllTimeSteps, RigContourMapProjection::maxValue( aggregatedResults ) ); } } } diff --git a/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.h b/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.h index ea6926e229..f78d0079e0 100644 --- a/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.h +++ b/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.h @@ -21,19 +21,13 @@ #include "RigFemPart.h" #include "RigFemResultAddress.h" -#include "RimCheckableNamedObject.h" #include "RimContourMapProjection.h" -#include "cafDisplayCoordTransform.h" -#include "cafPdmChildField.h" #include "cafPdmField.h" #include "cafPdmObject.h" #include "cvfArray.h" #include "cvfBoundingBox.h" -#include "cvfGeometryBuilderFaceList.h" -#include "cvfString.h" -#include "cvfVector2.h" class RimGeoMechContourMapView; class RimGeoMechCase; @@ -58,27 +52,19 @@ public: double sampleSpacing() const override; protected: - using CellIndexAndResult = RimContourMapProjection::CellIndexAndResult; - // GeoMech implementation specific data generation methods cvf::ref getCellVisibility() const override; cvf::BoundingBox calculateExpandedPorBarBBox( int timeStep, int frameIndex ) const; void updateGridInformation() override; - std::vector getMapCellVisibility() override; - std::vector retrieveParameterWeights() override; - std::vector generateResults( int viewerStepIndex ) override; + + std::vector retrieveParameterWeights() override; + std::vector generateResults( int timeStep ) const override; + void generateAndSaveResults( int timeStep ) override; std::vector - generateResultsFromAddress( RigFemResultAddress resultAddress, const std::vector& mapCellVisibility, int viewerStepIndex ); - bool resultVariableChanged() const override; - void clearResultVariable() override; - RimGridView* baseView() const override; - std::vector findIntersectingCells( const cvf::BoundingBox& bbox ) const override; - size_t kLayer( size_t globalCellIdx ) const override; - size_t kLayers() const override; - double calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const override; - double calculateRayLengthInCell( size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint ) const override; - double getParameterWeightForCell( size_t globalCellIdx, const std::vector& parameterWeights ) const override; - std::vector gridCellValues( RigFemResultAddress resAddr, std::vector& resultValues ) const; + generateResultsFromAddress( RigFemResultAddress resultAddress, const std::vector& mapCellVisibility, int viewerStepIndex ) const; + bool resultVariableChanged() const override; + void clearResultVariable() override; + RimGridView* baseView() const override; RimGeoMechCase* geoMechCase() const; RimGeoMechContourMapView* view() const; @@ -100,5 +86,4 @@ protected: cvf::ref m_femPart; cvf::cref m_femPartGrid; RigFemResultAddress m_currentResultAddr; - size_t m_kLayers; }; diff --git a/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp b/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp index d95a581789..3e0121ded4 100644 --- a/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp +++ b/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp @@ -19,11 +19,11 @@ #include "RimContourMapProjection.h" #include "RiaOpenMPTools.h" -#include "RiaWeightedMeanCalculator.h" #include "RigCellGeometryTools.h" #include "RigContourMapCalculator.h" #include "RigContourMapGrid.h" +#include "RigContourMapProjection.h" #include "RimCase.h" #include "RimGridView.h" @@ -37,10 +37,10 @@ #include "cafProgressInfo.h" #include "cvfArray.h" -#include "cvfGeometryTools.h" #include "cvfGeometryUtils.h" #include "cvfScalarMapper.h" #include "cvfStructGridGeometryGenerator.h" +#include "cvfVector2.h" #include @@ -66,6 +66,7 @@ void RimContourMapProjection::ResultAggregation::setUp() setDefault( RigContourMapCalculator::RESULTS_MEAN_VALUE ); } } // namespace caf + CAF_PDM_ABSTRACT_SOURCE_INIT( RimContourMapProjection, "RimContourMapProjection" ); //-------------------------------------------------------------------------------------------------- @@ -76,6 +77,7 @@ RimContourMapProjection::RimContourMapProjection() , m_currentResultTimestep( -1 ) , m_minResultAllTimeSteps( std::numeric_limits::infinity() ) , m_maxResultAllTimeSteps( -std::numeric_limits::infinity() ) + { CAF_PDM_InitObject( "RimContourMapProjection", ":/2DMapProjection16x16.png" ); @@ -106,22 +108,27 @@ void RimContourMapProjection::generateResultsIfNecessary( int timeStep ) { caf::ProgressInfo progress( 100, "Generate Results", true ); - updateGridInformation(); + if ( !m_contourMapGrid || !m_contourMapProjection ) updateGridInformation(); + + auto cellVisibility = getCellVisibility(); + m_contourMapProjection->setCellVisibility( cellVisibility ); + progress.setProgress( 10 ); - if ( gridMappingNeedsUpdating() || mapCellVisibilityNeedsUpdating() || resultVariableChanged() ) + if ( gridMappingNeedsUpdating() || mapCellVisibilityNeedsUpdating( timeStep ) || resultVariableChanged() ) { clearResults(); clearTimeStepRange(); - m_cellGridIdxVisibility = getCellVisibility(); + auto cellVisibility = getCellVisibility(); + m_contourMapProjection->setCellVisibility( cellVisibility ); if ( gridMappingNeedsUpdating() ) { - m_projected3dGridIndices = RigContourMapCalculator::generateGridMapping( *this, *m_contourMapGrid ); + m_contourMapProjection->generateGridMapping( m_resultAggregation(), retrieveParameterWeights() ); } progress.setProgress( 20 ); - m_mapCellVisibility = getMapCellVisibility(); + m_mapCellVisibility = m_contourMapProjection->getMapCellVisibility( timeStep, m_resultAggregation() ); progress.setProgress( 30 ); } else @@ -132,7 +139,8 @@ void RimContourMapProjection::generateResultsIfNecessary( int timeStep ) if ( resultsNeedsUpdating( timeStep ) ) { clearGeometry(); - m_aggregatedResults = generateResults( timeStep ); + generateAndSaveResults( timeStep ); + progress.setProgress( 80 ); generateVertexResults(); } @@ -279,7 +287,8 @@ QString RimContourMapProjection::currentTimeStepName() const //-------------------------------------------------------------------------------------------------- double RimContourMapProjection::maxValue() const { - return maxValue( m_aggregatedResults ); + if ( m_contourMapProjection ) return m_contourMapProjection->maxValue(); + return -std::numeric_limits::infinity(); } //-------------------------------------------------------------------------------------------------- @@ -287,7 +296,8 @@ double RimContourMapProjection::maxValue() const //-------------------------------------------------------------------------------------------------- double RimContourMapProjection::minValue() const { - return minValue( m_aggregatedResults ); + if ( m_contourMapProjection ) return m_contourMapProjection->minValue(); + return std::numeric_limits::infinity(); } //-------------------------------------------------------------------------------------------------- @@ -295,7 +305,8 @@ double RimContourMapProjection::minValue() const //-------------------------------------------------------------------------------------------------- double RimContourMapProjection::meanValue() const { - return sumAllValues() / numberOfValidCells(); + if ( m_contourMapProjection ) return m_contourMapProjection->meanValue(); + return std::numeric_limits::infinity(); } //-------------------------------------------------------------------------------------------------- @@ -303,24 +314,8 @@ double RimContourMapProjection::meanValue() const //-------------------------------------------------------------------------------------------------- double RimContourMapProjection::sumAllValues() const { - double sum = 0.0; - - for ( size_t index = 0; index < m_aggregatedResults.size(); ++index ) - { - if ( m_aggregatedResults[index] != std::numeric_limits::infinity() ) - { - sum += m_aggregatedResults[index]; - } - } - return sum; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -cvf::Vec2ui RimContourMapProjection::numberOfElementsIJ() const -{ - return m_contourMapGrid->numberOfElementsIJ(); + if ( m_contourMapProjection ) return m_contourMapProjection->sumAllValues(); + return 0.0; } //-------------------------------------------------------------------------------------------------- @@ -328,7 +323,8 @@ cvf::Vec2ui RimContourMapProjection::numberOfElementsIJ() const //-------------------------------------------------------------------------------------------------- cvf::Vec2ui RimContourMapProjection::numberOfVerticesIJ() const { - return m_contourMapGrid->numberOfVerticesIJ(); + if ( m_contourMapGrid ) return m_contourMapGrid->numberOfVerticesIJ(); + return cvf::Vec2ui( 0, 0 ); } //-------------------------------------------------------------------------------------------------- @@ -344,11 +340,7 @@ bool RimContourMapProjection::isColumnResult() const //-------------------------------------------------------------------------------------------------- double RimContourMapProjection::valueAtVertex( uint i, uint j ) const { - size_t index = m_contourMapGrid->vertexIndexFromIJ( i, j ); - if ( index < numberOfVertices() ) - { - return m_aggregatedVertexResults.at( index ); - } + if ( m_contourMapProjection ) return m_contourMapProjection->valueAtVertex( i, j ); return std::numeric_limits::infinity(); } @@ -357,7 +349,8 @@ double RimContourMapProjection::valueAtVertex( uint i, uint j ) const //-------------------------------------------------------------------------------------------------- uint RimContourMapProjection::numberOfCells() const { - return m_contourMapGrid->numberOfCells(); + if ( m_contourMapGrid ) return m_contourMapGrid->numberOfCells(); + return 0u; } //-------------------------------------------------------------------------------------------------- @@ -365,16 +358,8 @@ uint RimContourMapProjection::numberOfCells() const //-------------------------------------------------------------------------------------------------- uint RimContourMapProjection::numberOfValidCells() const { - uint validCount = 0u; - for ( uint i = 0; i < numberOfCells(); ++i ) - { - cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( i ); - if ( hasResultInCell( ij.x(), ij.y() ) ) - { - validCount++; - } - } - return validCount; + if ( m_contourMapProjection ) return m_contourMapProjection->numberOfValidCells(); + return 0u; } //-------------------------------------------------------------------------------------------------- @@ -382,8 +367,8 @@ uint RimContourMapProjection::numberOfValidCells() const //-------------------------------------------------------------------------------------------------- size_t RimContourMapProjection::numberOfVertices() const { - cvf::Vec2ui gridSize = numberOfVerticesIJ(); - return static_cast( gridSize.x() ) * static_cast( gridSize.y() ); + if ( m_contourMapGrid ) return m_contourMapGrid->numberOfVertices(); + return 0; } //-------------------------------------------------------------------------------------------------- @@ -391,22 +376,7 @@ size_t RimContourMapProjection::numberOfVertices() const //-------------------------------------------------------------------------------------------------- bool RimContourMapProjection::checkForMapIntersection( const cvf::Vec3d& domainPoint3d, cvf::Vec2d* contourMapPoint, double* valueAtPoint ) const { - CVF_TIGHT_ASSERT( contourMapPoint ); - CVF_TIGHT_ASSERT( valueAtPoint ); - - const cvf::Vec3d& minPoint = m_contourMapGrid->expandedBoundingBox().min(); - cvf::Vec3d mapPos3d = domainPoint3d - minPoint; - cvf::Vec2d mapPos2d( mapPos3d.x(), mapPos3d.y() ); - cvf::Vec2d gridorigin( minPoint.x(), minPoint.y() ); - - double value = interpolateValue( mapPos2d ); - if ( value != std::numeric_limits::infinity() ) - { - *valueAtPoint = value; - *contourMapPoint = mapPos2d + gridorigin; - - return true; - } + if ( m_contourMapProjection ) return m_contourMapProjection->checkForMapIntersection( domainPoint3d, contourMapPoint, valueAtPoint ); return false; } @@ -439,8 +409,8 @@ size_t RimContourMapProjection::gridResultIndex( size_t globalCellIdx ) const //-------------------------------------------------------------------------------------------------- double RimContourMapProjection::calculateValueInMapCell( uint i, uint j, const std::vector& gridCellValues ) const { - const std::vector>& matchingCells = cellsAtIJ( i, j ); - return RigContourMapCalculator::calculateValueInMapCell( *this, matchingCells, gridCellValues, m_resultAggregation() ); + const std::vector>& matchingCells = m_contourMapProjection->cellsAtIJ( i, j ); + return RigContourMapCalculator::calculateValueInMapCell( *m_contourMapProjection, matchingCells, gridCellValues, m_resultAggregation() ); } //-------------------------------------------------------------------------------------------------- @@ -448,16 +418,19 @@ double RimContourMapProjection::calculateValueInMapCell( uint i, uint j, const s //-------------------------------------------------------------------------------------------------- bool RimContourMapProjection::gridMappingNeedsUpdating() const { - if ( m_projected3dGridIndices.size() != numberOfCells() ) return true; + if ( !m_contourMapProjection ) return true; - if ( m_cellGridIdxVisibility.isNull() ) return true; + if ( m_contourMapProjection->projected3dGridIndices().size() != numberOfCells() ) return true; + + auto cellGridIdxVisibility = m_contourMapProjection->getCellVisibility(); + if ( cellGridIdxVisibility.isNull() ) return true; cvf::ref currentVisibility = getCellVisibility(); - CVF_ASSERT( currentVisibility->size() == m_cellGridIdxVisibility->size() ); + CVF_ASSERT( currentVisibility->size() == cellGridIdxVisibility->size() ); for ( size_t i = 0; i < currentVisibility->size(); ++i ) { - if ( ( *currentVisibility )[i] != ( *m_cellGridIdxVisibility )[i] ) return true; + if ( ( *currentVisibility )[i] != ( *cellGridIdxVisibility )[i] ) return true; } return false; @@ -468,8 +441,10 @@ bool RimContourMapProjection::gridMappingNeedsUpdating() const //-------------------------------------------------------------------------------------------------- bool RimContourMapProjection::resultsNeedsUpdating( int timeStep ) const { - return ( m_aggregatedResults.size() != numberOfCells() || m_aggregatedVertexResults.size() != numberOfVertices() || - timeStep != m_currentResultTimestep ); + if ( !m_contourMapProjection ) return true; + + return ( m_contourMapProjection->aggregatedResults().size() != numberOfCells() || + m_contourMapProjection->aggregatedVertexResults().size() != numberOfVertices() || timeStep != m_currentResultTimestep ); } //-------------------------------------------------------------------------------------------------- @@ -480,15 +455,6 @@ bool RimContourMapProjection::geometryNeedsUpdating() const return m_contourPolygons.empty() || m_trianglesWithVertexValues.empty(); } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -bool RimContourMapProjection::resultRangeIsValid() const -{ - return m_minResultAllTimeSteps != std::numeric_limits::infinity() && - m_maxResultAllTimeSteps != -std::numeric_limits::infinity(); -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -496,7 +462,8 @@ void RimContourMapProjection::clearGridMapping() { clearResults(); clearTimeStepRange(); - m_projected3dGridIndices.clear(); + + if ( m_contourMapProjection ) m_contourMapProjection->clearGridMapping(); m_mapCellVisibility.clear(); } @@ -507,81 +474,12 @@ void RimContourMapProjection::clearResults() { clearGeometry(); - m_aggregatedResults.clear(); - m_aggregatedVertexResults.clear(); + if ( m_contourMapProjection ) m_contourMapProjection->clearResults(); m_currentResultTimestep = -1; clearResultVariable(); } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RimContourMapProjection::clearTimeStepRange() -{ - m_minResultAllTimeSteps = std::numeric_limits::infinity(); - m_maxResultAllTimeSteps = -std::numeric_limits::infinity(); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimContourMapProjection::maxValue( const std::vector& aggregatedResults ) const -{ - double maxV = -std::numeric_limits::infinity(); - - for ( size_t index = 0; index < aggregatedResults.size(); ++index ) - { - if ( aggregatedResults[index] != std::numeric_limits::infinity() ) - { - maxV = std::max( maxV, aggregatedResults[index] ); - } - } - return maxV; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimContourMapProjection::minValue( const std::vector& aggregatedResults ) const -{ - double minV = std::numeric_limits::infinity(); - - for ( size_t index = 0; index < aggregatedResults.size(); ++index ) - { - if ( aggregatedResults[index] != std::numeric_limits::infinity() ) - { - minV = std::min( minV, aggregatedResults[index] ); - } - } - return minV; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::pair RimContourMapProjection::minmaxValuesAllTimeSteps() -{ - if ( !resultRangeIsValid() ) - { - clearTimeStepRange(); - - m_minResultAllTimeSteps = std::min( m_minResultAllTimeSteps, minValue( m_aggregatedResults ) ); - m_maxResultAllTimeSteps = std::max( m_maxResultAllTimeSteps, maxValue( m_aggregatedResults ) ); - - for ( int i = 0; i < (int)baseView()->ownerCase()->timeStepStrings().size() - 1; ++i ) - { - if ( i != m_currentResultTimestep ) - { - std::vector aggregatedResults = generateResults( i ); - m_minResultAllTimeSteps = std::min( m_minResultAllTimeSteps, minValue( aggregatedResults ) ); - m_maxResultAllTimeSteps = std::max( m_maxResultAllTimeSteps, maxValue( aggregatedResults ) ); - } - } - } - return std::make_pair( m_minResultAllTimeSteps, m_maxResultAllTimeSteps ); -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -590,20 +488,14 @@ cvf::ref RimContourMapProjection::getCellVisibility() const return baseView()->currentTotalCellVisibility(); } -//-------------------------------------------------------------------------------------------------- -/// Empty default implementation -//-------------------------------------------------------------------------------------------------- -std::vector RimContourMapProjection::getMapCellVisibility() -{ - return std::vector( numberOfCells(), true ); -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -bool RimContourMapProjection::mapCellVisibilityNeedsUpdating() +bool RimContourMapProjection::mapCellVisibilityNeedsUpdating( int timestep ) { - std::vector mapCellVisiblity = getMapCellVisibility(); + if ( m_contourMapProjection ) return true; + + std::vector mapCellVisiblity = m_contourMapProjection->getMapCellVisibility( timestep, m_resultAggregation() ); return !( mapCellVisiblity == m_mapCellVisibility ); } @@ -612,17 +504,7 @@ bool RimContourMapProjection::mapCellVisibilityNeedsUpdating() //-------------------------------------------------------------------------------------------------- void RimContourMapProjection::generateVertexResults() { - size_t nCells = numberOfCells(); - if ( nCells != m_aggregatedResults.size() ) return; - - size_t nVertices = numberOfVertices(); - m_aggregatedVertexResults = std::vector( nVertices, std::numeric_limits::infinity() ); -#pragma omp parallel for - for ( int index = 0; index < static_cast( nVertices ); ++index ) - { - cvf::Vec2ui ij = m_contourMapGrid->ijFromVertexIndex( index ); - m_aggregatedVertexResults[index] = calculateValueAtVertex( ij.x(), ij.y() ); - } + if ( m_contourMapProjection ) m_contourMapProjection->generateVertexResults(); } //-------------------------------------------------------------------------------------------------- @@ -668,6 +550,8 @@ void RimContourMapProjection::generateTrianglesWithVertexValues() std::vector>> threadTriangles( numberOfThreads ); + const std::vector& aggregatedVertexResults = m_contourMapProjection->aggregatedVertexResults(); + #pragma omp parallel { int myThread = RiaOpenMPTools::currentThreadIndex(); @@ -681,9 +565,9 @@ void RimContourMapProjection::generateTrianglesWithVertexValues() bool anyValidVertex = false; for ( size_t n = 0; n < 3; ++n ) { - uint vn = ( *faceList )[i + n]; - double value = vn < m_aggregatedVertexResults.size() ? m_aggregatedVertexResults[vn] : std::numeric_limits::infinity(); - triangle[n] = vertices[vn]; + uint vn = ( *faceList )[i + n]; + double value = vn < aggregatedVertexResults.size() ? aggregatedVertexResults[vn] : std::numeric_limits::infinity(); + triangle[n] = vertices[vn]; triangleWithValues[n] = cvf::Vec4d( vertices[vn], value ); if ( value != std::numeric_limits::infinity() ) { @@ -804,7 +688,7 @@ void RimContourMapProjection::generateTrianglesWithVertexValues() } if ( value == std::numeric_limits::infinity() ) { - value = interpolateValue( cvf::Vec2d( localVertex.x(), localVertex.y() ) ); + value = m_contourMapProjection->interpolateValue( cvf::Vec2d( localVertex.x(), localVertex.y() ) ); if ( value == std::numeric_limits::infinity() ) { value = contourLevels[c]; @@ -841,7 +725,7 @@ void RimContourMapProjection::generateTrianglesWithVertexValues() allTrianglesThisLevel.insert( allTrianglesThisLevel.end(), threadTriangles[i][c].begin(), threadTriangles[i][c].end() ); } - double triangleAreasThisLevel = sumTriangleAreas( allTrianglesThisLevel ); + double triangleAreasThisLevel = RigContourMapProjection::sumTriangleAreas( allTrianglesThisLevel ); if ( c >= m_contourLevelCumulativeAreas.size() || triangleAreasThisLevel > 1.0e-3 * m_contourLevelCumulativeAreas[c] ) { trianglesPerLevel[c] = allTrianglesThisLevel; @@ -867,6 +751,8 @@ void RimContourMapProjection::generateContourPolygons() std::vector contourLevels; if ( resultRangeIsValid() && legendConfig()->mappingMode() != RimRegularLegendConfig::MappingType::CATEGORY_INTEGER ) { + const std::vector& aggregatedVertexResults = m_contourMapProjection->aggregatedVertexResults(); + legendConfig()->scalarMapper()->majorTickValues( &contourLevels ); int nContourLevels = static_cast( contourLevels.size() ); @@ -900,7 +786,7 @@ void RimContourMapProjection::generateContourPolygons() } std::vector unorderedLineSegmentsPerLevel = - caf::ContourLines::create( m_aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels ); + caf::ContourLines::create( aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels ); contourPolygons = std::vector( unorderedLineSegmentsPerLevel.size() ); const double areaThreshold = 1.5 * ( sampleSpacing() * sampleSpacing() ) / ( sampleSpacingFactor() * sampleSpacingFactor() ); @@ -948,23 +834,6 @@ void RimContourMapProjection::generateContourPolygons() m_contourPolygons = contourPolygons; } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimContourMapProjection::sumTriangleAreas( const std::vector& triangles ) -{ - double sumArea = 0.0; - for ( size_t i = 0; i < triangles.size(); i += 3 ) - { - cvf::Vec3d v1( triangles[i].x(), triangles[i].y(), triangles[i].z() ); - cvf::Vec3d v2( triangles[i + 1].x(), triangles[i + 1].y(), triangles[i + 1].z() ); - cvf::Vec3d v3( triangles[i + 2].x(), triangles[i + 2].y(), triangles[i + 2].z() ); - double area = 0.5 * ( ( v3 - v1 ) ^ ( v2 - v1 ) ).length(); - sumArea += area; - } - return sumArea; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -981,130 +850,13 @@ bool RimContourMapProjection::isStraightSummationResult() const return RigContourMapCalculator::isStraightSummationResult( m_resultAggregation() ); } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimContourMapProjection::interpolateValue( const cvf::Vec2d& gridPos2d ) const -{ - cvf::Vec2ui cellContainingPoint = m_contourMapGrid->ijFromLocalPos( gridPos2d ); - cvf::Vec2d cellCenter = m_contourMapGrid->cellCenterPosition( cellContainingPoint.x(), cellContainingPoint.y() ); - - std::array x; - x[0] = cvf::Vec3d( cellCenter + cvf::Vec2d( -sampleSpacing() * 0.5, -sampleSpacing() * 0.5 ), 0.0 ); - x[1] = cvf::Vec3d( cellCenter + cvf::Vec2d( sampleSpacing() * 0.5, -sampleSpacing() * 0.5 ), 0.0 ); - x[2] = cvf::Vec3d( cellCenter + cvf::Vec2d( sampleSpacing() * 0.5, sampleSpacing() * 0.5 ), 0.0 ); - x[3] = cvf::Vec3d( cellCenter + cvf::Vec2d( -sampleSpacing() * 0.5, sampleSpacing() * 0.5 ), 0.0 ); - - cvf::Vec4d baryCentricCoords = cvf::GeometryTools::barycentricCoords( x[0], x[1], x[2], x[3], cvf::Vec3d( gridPos2d, 0.0 ) ); - - std::array v; - v[0] = cellContainingPoint; - v[1] = cvf::Vec2ui( cellContainingPoint.x() + 1u, cellContainingPoint.y() ); - v[2] = cvf::Vec2ui( cellContainingPoint.x() + 1u, cellContainingPoint.y() + 1u ); - v[3] = cvf::Vec2ui( cellContainingPoint.x(), cellContainingPoint.y() + 1u ); - - std::array vertexValues; - double validBarycentricCoordsSum = 0.0; - for ( int i = 0; i < 4; ++i ) - { - double vertexValue = valueAtVertex( v[i].x(), v[i].y() ); - if ( vertexValue == std::numeric_limits::infinity() ) - { - return std::numeric_limits::infinity(); - } - else - { - vertexValues[i] = vertexValue; - validBarycentricCoordsSum += baryCentricCoords[i]; - } - } - - if ( validBarycentricCoordsSum < 1.0e-8 ) - { - return std::numeric_limits::infinity(); - } - - // Calculate final value - double value = 0.0; - for ( int i = 0; i < 4; ++i ) - { - value += baryCentricCoords[i] / validBarycentricCoordsSum * vertexValues[i]; - } - - return value; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimContourMapProjection::valueInCell( uint i, uint j ) const -{ - size_t index = m_contourMapGrid->cellIndexFromIJ( i, j ); - if ( index < numberOfCells() ) - { - return m_aggregatedResults.at( index ); - } - return std::numeric_limits::infinity(); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -bool RimContourMapProjection::hasResultInCell( uint i, uint j ) const -{ - return !cellsAtIJ( i, j ).empty(); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimContourMapProjection::calculateValueAtVertex( uint vi, uint vj ) const -{ - std::vector averageIs; - std::vector averageJs; - - if ( vi > 0u ) averageIs.push_back( vi - 1 ); - if ( vj > 0u ) averageJs.push_back( vj - 1 ); - if ( vi < m_contourMapGrid->mapSize().x() ) averageIs.push_back( vi ); - if ( vj < m_contourMapGrid->mapSize().y() ) averageJs.push_back( vj ); - - RiaWeightedMeanCalculator calc; - for ( uint j : averageJs ) - { - for ( uint i : averageIs ) - { - if ( hasResultInCell( i, j ) ) - { - calc.addValueAndWeight( valueInCell( i, j ), 1.0 ); - } - } - } - if ( calc.validAggregatedWeight() ) - { - return calc.weightedMean(); - } - return std::numeric_limits::infinity(); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector> RimContourMapProjection::cellsAtIJ( uint i, uint j ) const -{ - size_t cellIndex = m_contourMapGrid->cellIndexFromIJ( i, j ); - if ( cellIndex < m_projected3dGridIndices.size() ) - { - return m_projected3dGridIndices[cellIndex]; - } - return std::vector>(); -} - //-------------------------------------------------------------------------------------------------- /// Vertex positions in local coordinates (add origin2d.x() for UTM x) //-------------------------------------------------------------------------------------------------- std::vector RimContourMapProjection::xVertexPositions() const { - return m_contourMapGrid->xVertexPositions(); + if ( m_contourMapGrid ) return m_contourMapGrid->xVertexPositions(); + return {}; } //-------------------------------------------------------------------------------------------------- @@ -1112,7 +864,8 @@ std::vector RimContourMapProjection::xVertexPositions() const //-------------------------------------------------------------------------------------------------- std::vector RimContourMapProjection::yVertexPositions() const { - return m_contourMapGrid->yVertexPositions(); + if ( m_contourMapGrid ) return m_contourMapGrid->yVertexPositions(); + return {}; } //-------------------------------------------------------------------------------------------------- @@ -1207,3 +960,21 @@ void RimContourMapProjection::defineUiTreeOrdering( caf::PdmUiTreeOrdering& uiTr void RimContourMapProjection::initAfterRead() { } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RimContourMapProjection::resultRangeIsValid() const +{ + return m_minResultAllTimeSteps != std::numeric_limits::infinity() && + m_maxResultAllTimeSteps != -std::numeric_limits::infinity(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimContourMapProjection::clearTimeStepRange() +{ + m_minResultAllTimeSteps = std::numeric_limits::infinity(); + m_maxResultAllTimeSteps = -std::numeric_limits::infinity(); +} diff --git a/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h b/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h index 9ecc867255..ac97552d4e 100644 --- a/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h +++ b/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h @@ -27,12 +27,12 @@ #include "cafPdmObject.h" #include "cvfArray.h" -#include "cvfBoundingBox.h" #include "cvfVector2.h" class RigContourMapGrid; class RimGridView; class RimRegularLegendConfig; +class RigContourMapProjection; //================================================================================================== /// @@ -43,8 +43,6 @@ class RimContourMapProjection : public RimCheckableNamedObject CAF_PDM_HEADER_INIT; public: - using CellIndexAndResult = std::pair; - using ResultAggregation = caf::AppEnum; using ContourPolygons = std::vector; @@ -78,7 +76,6 @@ public: double meanValue() const; double sumAllValues() const; - cvf::Vec2ui numberOfElementsIJ() const; cvf::Vec2ui numberOfVerticesIJ() const; bool isColumnResult() const; @@ -108,26 +105,19 @@ public: virtual std::vector retrieveParameterWeights() = 0; - virtual size_t kLayer( size_t globalCellIdx ) const = 0; - virtual size_t kLayers() const = 0; - virtual std::vector findIntersectingCells( const cvf::BoundingBox& bbox ) const = 0; - virtual double calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const = 0; - virtual double calculateRayLengthInCell( size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint ) const = 0; - virtual double getParameterWeightForCell( size_t globalCellIdx, const std::vector& parameterWeights ) const = 0; - virtual cvf::ref getCellVisibility() const; protected: // Protected virtual methods to be overridden by Eclipse and Geo-mechanical contour map implementations - virtual void updateGridInformation() = 0; - virtual std::vector generateResults( int timeStep ) = 0; - virtual bool resultVariableChanged() const = 0; - virtual void clearResultVariable() = 0; - virtual RimGridView* baseView() const = 0; + virtual void updateGridInformation() = 0; + virtual std::vector generateResults( int timeStep ) const = 0; + virtual void generateAndSaveResults( int timeStep ) = 0; + virtual bool resultVariableChanged() const = 0; + virtual void clearResultVariable() = 0; + virtual RimGridView* baseView() const = 0; double calculateValueInMapCell( uint i, uint j, const std::vector& gridCellValues ) const; -protected: // Keep track of whether cached data needs updating bool gridMappingNeedsUpdating() const; bool resultsNeedsUpdating( int timeStep ) const; @@ -137,30 +127,14 @@ protected: void clearResults(); void clearTimeStepRange(); - double maxValue( const std::vector& aggregatedResults ) const; - double minValue( const std::vector& aggregatedResults ) const; + virtual std::pair minmaxValuesAllTimeSteps() = 0; - virtual std::pair minmaxValuesAllTimeSteps(); - - virtual std::vector getMapCellVisibility(); - bool mapCellVisibilityNeedsUpdating(); - static std::vector>> generateGridMapping( RimContourMapProjection& contourMapProjection, - const RigContourMapGrid& contourMapGrid ); + bool mapCellVisibilityNeedsUpdating( int timeStep ); void generateVertexResults(); void generateTrianglesWithVertexValues(); void generateContourPolygons(); - static double sumTriangleAreas( const std::vector& triangles ); - - double interpolateValue( const cvf::Vec2d& gridPosition2d ) const; - double valueInCell( uint i, uint j ) const; - bool hasResultInCell( uint i, uint j ) const; - double calculateValueAtVertex( uint i, uint j ) const; - - // Cell index and position conversion - std::vector cellsAtIJ( uint i, uint j ) const; - double gridEdgeOffset() const; virtual void updateAfterResultGeneration( int timeStep ) = 0; @@ -180,11 +154,6 @@ protected: caf::PdmField m_showContourLabels; caf::PdmField m_smoothContourLines; - cvf::ref m_cellGridIdxVisibility; - std::vector m_aggregatedResults; - std::vector m_aggregatedVertexResults; - std::vector>> m_projected3dGridIndices; - cvf::Vec2d m_pickPoint; std::vector m_contourPolygons; std::vector m_contourLevelCumulativeAreas; @@ -195,5 +164,6 @@ protected: double m_minResultAllTimeSteps; double m_maxResultAllTimeSteps; - std::unique_ptr m_contourMapGrid; + std::unique_ptr m_contourMapGrid; + std::unique_ptr m_contourMapProjection; }; diff --git a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp index 4346ab9c18..d23a540f8a 100644 --- a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp +++ b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp @@ -18,15 +18,13 @@ #include "RimEclipseContourMapProjection.h" -#include "RigActiveCellInfo.h" +#include "RiaPorosityModel.h" #include "RigCaseCellResultsData.h" -#include "RigCell.h" -#include "RigCellGeometryTools.h" #include "RigContourMapCalculator.h" #include "RigContourMapGrid.h" #include "RigEclipseCaseData.h" +#include "RigEclipseContourMapProjection.h" #include "RigEclipseResultAddress.h" -#include "RigHexIntersectionTools.h" #include "RigMainGrid.h" #include "Rim3dView.h" @@ -37,8 +35,6 @@ #include "RimEclipseView.h" #include "RimRegularLegendConfig.h" -#include - CAF_PDM_SOURCE_INIT( RimEclipseContourMapProjection, "RimEclipseContourMapProjection" ); //-------------------------------------------------------------------------------------------------- @@ -46,8 +42,6 @@ CAF_PDM_SOURCE_INIT( RimEclipseContourMapProjection, "RimEclipseContourMapProjec //-------------------------------------------------------------------------------------------------- RimEclipseContourMapProjection::RimEclipseContourMapProjection() : RimContourMapProjection() - , m_kLayers( 0u ) - , m_useActiveCellInfo( true ) { CAF_PDM_InitObject( "RimEclipseContourMapProjection", ":/2DMapProjection16x16.png" ); @@ -109,8 +103,8 @@ void RimEclipseContourMapProjection::updateLegend() { RimEclipseCellColors* cellColors = view()->cellResult(); - double minVal = minValue( m_aggregatedResults ); - double maxVal = maxValue( m_aggregatedResults ); + double minVal = minValue(); + double maxVal = maxValue(); auto [minValAllTimeSteps, maxValAllTimeSteps] = minmaxValuesAllTimeSteps(); @@ -138,13 +132,16 @@ void RimEclipseContourMapProjection::updateLegend() //-------------------------------------------------------------------------------------------------- double RimEclipseContourMapProjection::sampleSpacing() const { - if ( m_mainGrid.notNull() ) + if ( auto ec = eclipseCase() ) { - return m_relativeSampleSpacing * m_mainGrid->characteristicIJCellSize(); + if ( auto mainGrid = ec->mainGrid() ) + { + return m_relativeSampleSpacing * mainGrid->characteristicIJCellSize(); + } } + return 0.0; } - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -162,76 +159,43 @@ void RimEclipseContourMapProjection::clearGridMappingAndRedraw() //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::vector RimEclipseContourMapProjection::generateResults( int timeStep ) +std::vector RimEclipseContourMapProjection::generateResults( int timeStep ) const { m_weightingResult->loadResult(); - size_t nCells = numberOfCells(); - - std::vector aggregatedResults = std::vector( nCells, std::numeric_limits::infinity() ); - - RimEclipseCellColors* cellColors = view()->cellResult(); - auto gridCellResult = view()->currentGridCellResults(); + if ( m_contourMapProjection ) { - { - auto resultAdr = cellColors->eclipseResultAddress(); - if ( resultAdr.isValid() && gridCellResult->hasResultEntry( resultAdr ) ) - m_useActiveCellInfo = gridCellResult->isUsingGlobalActiveIndex( resultAdr ); - } + RimEclipseCellColors* cellColors = view()->cellResult(); + RigEclipseResultAddress resAddr( cellColors->resultType(), + cellColors->resultVariable(), + cellColors->timeLapseBaseTimeStep(), + cellColors->caseDiffIndex() ); - if ( !cellColors->isTernarySaturationSelected() ) - { - std::vector gridResultValues; - if ( isColumnResult() ) - { - m_currentResultName = ""; - gridCellResult->ensureKnownResultLoaded( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PORO" ) ); - gridCellResult->ensureKnownResultLoaded( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "NTG" ) ); - gridCellResult->ensureKnownResultLoaded( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "DZ" ) ); - if ( m_resultAggregation == RigContourMapCalculator::RESULTS_OIL_COLUMN || - m_resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN ) - { - gridCellResult->ensureKnownResultLoaded( - RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::soil() ) ); - } - if ( m_resultAggregation == RigContourMapCalculator::RESULTS_GAS_COLUMN || - m_resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN ) - { - gridCellResult->ensureKnownResultLoaded( - RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() ) ); - } - gridResultValues = calculateColumnResult( m_resultAggregation() ); - } - else - { - if ( cellColors->hasStaticResult() && timeStep > 0 ) timeStep = 0; - - m_currentResultName = cellColors->resultVariable(); - RigEclipseResultAddress resAddr( cellColors->resultType(), - cellColors->resultVariable(), - cellColors->timeLapseBaseTimeStep(), - cellColors->caseDiffIndex() ); - - // When loading a project file, grid calculator results are not computed the first time this function is - // called. Must check if result is loaded. See RimReloadCaseTools::updateAll3dViews() - if ( resAddr.isValid() && gridCellResult->hasResultEntry( resAddr ) && gridCellResult->isResultLoaded( resAddr ) ) - { - gridResultValues = gridCellResult->cellScalarResults( resAddr, timeStep ); - } - } - - if ( !gridResultValues.empty() ) - { -#pragma omp parallel for - for ( int index = 0; index < static_cast( nCells ); ++index ) - { - cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( index ); - aggregatedResults[index] = calculateValueInMapCell( ij.x(), ij.y(), gridResultValues ); - } - } - } + return dynamic_cast( m_contourMapProjection.get() ) + ->generateResults( resAddr, m_resultAggregation(), timeStep ); + } + + return {}; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RimEclipseContourMapProjection::generateAndSaveResults( int timeStep ) +{ + m_weightingResult->loadResult(); + + if ( m_contourMapProjection ) + { + RimEclipseCellColors* cellColors = view()->cellResult(); + RigEclipseResultAddress resAddr( cellColors->resultType(), + cellColors->resultVariable(), + cellColors->timeLapseBaseTimeStep(), + cellColors->caseDiffIndex() ); + + dynamic_cast( m_contourMapProjection.get() ) + ->generateAndSaveResults( resAddr, m_resultAggregation(), timeStep ); } - return aggregatedResults; } //-------------------------------------------------------------------------------------------------- @@ -258,80 +222,18 @@ void RimEclipseContourMapProjection::clearResultVariable() m_currentResultName = ""; } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector RimEclipseContourMapProjection::calculateColumnResult( ResultAggregation resultAggregation ) const -{ - const RigCaseCellResultsData* resultData = eclipseCase()->results( RiaDefines::PorosityModelType::MATRIX_MODEL ); - bool hasPoroResult = resultData->hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PORO" ) ); - bool hasNtgResult = resultData->hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "NTG" ) ); - bool hasDzResult = resultData->hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "DZ" ) ); - - if ( !( hasPoroResult && hasNtgResult && hasDzResult ) ) - { - return std::vector(); - } - - const std::vector& poroResults = - resultData->cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PORO" ), 0 ); - const std::vector& ntgResults = - resultData->cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "NTG" ), 0 ); - const std::vector& dzResults = - resultData->cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "DZ" ), 0 ); - - CVF_ASSERT( poroResults.size() == ntgResults.size() && ntgResults.size() == dzResults.size() ); - - int timeStep = view()->currentTimeStep(); - - std::vector resultValues( poroResults.size(), 0.0 ); - - if ( resultAggregation == RigContourMapCalculator::RESULTS_OIL_COLUMN || resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN ) - { - const std::vector& soilResults = - resultData->cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::soil() ), - timeStep ); - for ( size_t cellResultIdx = 0; cellResultIdx < resultValues.size(); ++cellResultIdx ) - { - resultValues[cellResultIdx] = soilResults[cellResultIdx]; - } - } - - if ( resultAggregation == RigContourMapCalculator::RESULTS_GAS_COLUMN || resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN ) - { - bool hasGasResult = - resultData->hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() ) ); - if ( hasGasResult ) - { - const std::vector& sgasResults = - resultData->cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() ), - timeStep ); - for ( size_t cellResultIdx = 0; cellResultIdx < resultValues.size(); ++cellResultIdx ) - { - resultValues[cellResultIdx] += sgasResults[cellResultIdx]; - } - } - } - - for ( size_t cellResultIdx = 0; cellResultIdx < resultValues.size(); ++cellResultIdx ) - { - resultValues[cellResultIdx] *= poroResults[cellResultIdx] * ntgResults[cellResultIdx] * dzResults[cellResultIdx]; - } - return resultValues; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimEclipseContourMapProjection::updateGridInformation() { - auto eclipseCase = this->eclipseCase(); - m_mainGrid = eclipseCase->eclipseCaseData()->mainGrid(); - m_activeCellInfo = eclipseCase->eclipseCaseData()->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL ); - m_kLayers = m_mainGrid->cellCountK(); + auto eclipseCase = this->eclipseCase(); + auto eclipseCaseData = eclipseCase->eclipseCaseData(); + auto resultData = eclipseCaseData->results( RiaDefines::PorosityModelType::MATRIX_MODEL ); cvf::BoundingBox gridBoundingBox = eclipseCase->activeCellsBoundingBox(); m_contourMapGrid = std::make_unique( gridBoundingBox, sampleSpacing() ); + m_contourMapProjection = std::make_unique( *m_contourMapGrid, *eclipseCaseData, *resultData ); } //-------------------------------------------------------------------------------------------------- @@ -376,108 +278,6 @@ RimGridView* RimEclipseContourMapProjection::baseView() const return view(); } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector RimEclipseContourMapProjection::findIntersectingCells( const cvf::BoundingBox& bbox ) const -{ - return m_mainGrid->findIntersectingCells( bbox ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -size_t RimEclipseContourMapProjection::kLayer( size_t globalCellIdx ) const -{ - const RigCell& cell = m_mainGrid->cell( globalCellIdx ); - size_t mainGridCellIdx = cell.mainGridCellIndex(); - size_t i, j, k; - m_mainGrid->ijkFromCellIndex( mainGridCellIdx, &i, &j, &k ); - return k; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -size_t RimEclipseContourMapProjection::kLayers() const -{ - return m_kLayers; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimEclipseContourMapProjection::calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const -{ - std::array hexCorners; - - const RigCell& cell = m_mainGrid->cell( globalCellIdx ); - - size_t localCellIdx = cell.gridLocalCellIndex(); - RigGridBase* localGrid = cell.hostGrid(); - - localGrid->cellCornerVertices( localCellIdx, hexCorners.data() ); - - cvf::BoundingBox overlapBBox; - std::array overlapCorners; - if ( RigCellGeometryTools::estimateHexOverlapWithBoundingBox( hexCorners, bbox, &overlapCorners, &overlapBBox ) ) - { - double overlapVolume = RigCellGeometryTools::calculateCellVolume( overlapCorners ); - return overlapVolume; - } - return 0.0; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimEclipseContourMapProjection::calculateRayLengthInCell( size_t globalCellIdx, - const cvf::Vec3d& highestPoint, - const cvf::Vec3d& lowestPoint ) const -{ - std::array hexCorners; - - RigCell cell = m_mainGrid->cell( globalCellIdx ); - - size_t localCellIdx = cell.gridLocalCellIndex(); - RigGridBase* localGrid = cell.hostGrid(); - - localGrid->cellCornerVertices( localCellIdx, hexCorners.data() ); - std::vector intersections; - - if ( RigHexIntersectionTools::lineHexCellIntersection( highestPoint, lowestPoint, hexCorners.data(), 0, &intersections ) ) - { - double lengthInCell = ( intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint ).length(); - return lengthInCell; - } - return 0.0; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RimEclipseContourMapProjection::getParameterWeightForCell( size_t cellResultIdx, const std::vector& cellWeights ) const -{ - if ( cellWeights.empty() ) return 1.0; - - double result = std::max( cellWeights[cellResultIdx], 0.0 ); - if ( result < 1.0e-6 ) - { - result = 0.0; - } - return result; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -size_t RimEclipseContourMapProjection::gridResultIndex( size_t globalCellIdx ) const -{ - if ( m_useActiveCellInfo ) return m_activeCellInfo->cellResultIndex( globalCellIdx ); - - return globalCellIdx; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -545,3 +345,23 @@ void RimEclipseContourMapProjection::initAfterRead() m_weightingResult->setEclipseCase( eclipseCase() ); } } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RimEclipseContourMapProjection::minmaxValuesAllTimeSteps() +{ + if ( !resultRangeIsValid() ) + { + clearTimeStepRange(); + + int timeStepCount = static_cast( eclipseCase()->timeStepStrings().size() ); + for ( int i = 0; i < (int)timeStepCount; ++i ) + { + std::vector aggregatedResults = generateResults( i ); + m_minResultAllTimeSteps = std::min( m_minResultAllTimeSteps, RigContourMapProjection::minValue( aggregatedResults ) ); + m_maxResultAllTimeSteps = std::max( m_maxResultAllTimeSteps, RigContourMapProjection::maxValue( aggregatedResults ) ); + } + } + return std::make_pair( m_minResultAllTimeSteps, m_maxResultAllTimeSteps ); +} diff --git a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.h b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.h index 0a5b975ee9..6129329231 100644 --- a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.h +++ b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.h @@ -24,8 +24,6 @@ #include "cafPdmField.h" #include "cafPdmObject.h" -#include "cvfBoundingBox.h" - class RigActiveCellInfo; class RigMainGrid; class RigContourMapGrid; @@ -57,28 +55,19 @@ public: double sampleSpacing() const override; protected: - using CellIndexAndResult = RimContourMapProjection::CellIndexAndResult; - void updateGridInformation() override; std::vector retrieveParameterWeights() override; - std::vector generateResults( int timeStep ) override; + std::vector generateResults( int timeStep ) const override; + void generateAndSaveResults( int timeStep ) override; bool resultVariableChanged() const override; void clearResultVariable() override; RimGridView* baseView() const override; - std::vector findIntersectingCells( const cvf::BoundingBox& bbox ) const override; - size_t kLayer( size_t globalCellIdx ) const override; - size_t kLayers() const override; - double calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const override; - double calculateRayLengthInCell( size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint ) const override; - double getParameterWeightForCell( size_t cellResultIdx, const std::vector& parameterWeights ) const override; - size_t gridResultIndex( size_t globalCellIdx ) const override; - - // Eclipse implementation specific data generation methods - std::vector calculateColumnResult( ResultAggregation resultAggregation ) const; RimEclipseCase* eclipseCase() const; RimEclipseContourMapView* view() const; + std::pair minmaxValuesAllTimeSteps() override; + void updateAfterResultGeneration( int timeStep ) override; protected: @@ -91,10 +80,5 @@ protected: caf::PdmField m_weightByParameter; caf::PdmChildField m_weightingResult; - cvf::ref m_mainGrid; - cvf::ref m_activeCellInfo; - size_t m_kLayers; - bool m_useActiveCellInfo; - QString m_currentResultName; }; diff --git a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapView.cpp b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapView.cpp index f830879552..abab36c1e7 100644 --- a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapView.cpp +++ b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapView.cpp @@ -272,6 +272,7 @@ void RimEclipseContourMapView::updateGeometry() { m_contourMapProjection->generateResultsIfNecessary( m_currentTimeStep() ); } + onUpdateLegends(); progress.setProgress( 30 ); } diff --git a/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake b/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake index 7c6b11ec86..1b775c4969 100644 --- a/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake +++ b/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake @@ -83,6 +83,10 @@ set(SOURCE_GROUP_HEADER_FILES ${CMAKE_CURRENT_LIST_DIR}/RigVisibleCategoriesCalculator.h ${CMAKE_CURRENT_LIST_DIR}/RigWbsParameter.h ${CMAKE_CURRENT_LIST_DIR}/RigWeightedMeanCalc.h + ${CMAKE_CURRENT_LIST_DIR}/RigContourMapCalculator.h + ${CMAKE_CURRENT_LIST_DIR}/RigContourMapProjection.h + ${CMAKE_CURRENT_LIST_DIR}/RigEclipseContourMapProjection.h + ${CMAKE_CURRENT_LIST_DIR}/RigGeoMechContourMapProjection.h ) set(SOURCE_GROUP_SOURCE_FILES @@ -166,6 +170,12 @@ set(SOURCE_GROUP_SOURCE_FILES ${CMAKE_CURRENT_LIST_DIR}/RigVisibleCategoriesCalculator.cpp ${CMAKE_CURRENT_LIST_DIR}/RigWbsParameter.cpp ${CMAKE_CURRENT_LIST_DIR}/RigWeightedMeanCalc.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigContourMapGrid.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigContourPolygonsTools.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigContourMapCalculator.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigContourMapProjection.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigEclipseContourMapProjection.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigGeoMechContourMapProjection.cpp ) list(APPEND CODE_HEADER_FILES ${SOURCE_GROUP_HEADER_FILES}) diff --git a/ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.cpp b/ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.cpp index d7a55db529..16ca635f72 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.cpp @@ -24,15 +24,16 @@ #include "RigContourMapGrid.h" -#include "RimCase.h" -#include "RimContourMapProjection.h" +#include "RigContourMapProjection.h" #include +#include +#include //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigContourMapCalculator::calculateValueInMapCell( const RimContourMapProjection& contourMapProjection, +double RigContourMapCalculator::calculateValueInMapCell( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues, ResultAggregationEnum resultAggregation ) @@ -70,7 +71,7 @@ double RigContourMapCalculator::calculateValueInMapCell( const RimContourMapProj //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigContourMapCalculator::calculateTopValue( const RimContourMapProjection& contourMapProjection, +double RigContourMapCalculator::calculateTopValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ) { @@ -88,7 +89,7 @@ double RigContourMapCalculator::calculateTopValue( const RimContourMapProjection //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigContourMapCalculator::calculateMeanValue( const RimContourMapProjection& contourMapProjection, +double RigContourMapCalculator::calculateMeanValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ) { @@ -111,7 +112,7 @@ double RigContourMapCalculator::calculateMeanValue( const RimContourMapProjectio //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigContourMapCalculator::calculateGeometricMeanValue( const RimContourMapProjection& contourMapProjection, +double RigContourMapCalculator::calculateGeometricMeanValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ) { @@ -138,7 +139,7 @@ double RigContourMapCalculator::calculateGeometricMeanValue( const RimContourMap //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigContourMapCalculator::calculateHarmonicMeanValue( const RimContourMapProjection& contourMapProjection, +double RigContourMapCalculator::calculateHarmonicMeanValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ) { @@ -165,7 +166,7 @@ double RigContourMapCalculator::calculateHarmonicMeanValue( const RimContourMapP //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigContourMapCalculator::calculateMaxValue( const RimContourMapProjection& contourMapProjection, +double RigContourMapCalculator::calculateMaxValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ) { @@ -188,7 +189,7 @@ double RigContourMapCalculator::calculateMaxValue( const RimContourMapProjection //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigContourMapCalculator::calculateMinValue( const RimContourMapProjection& contourMapProjection, +double RigContourMapCalculator::calculateMinValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ) { @@ -201,7 +202,10 @@ double RigContourMapCalculator::calculateMinValue( const RimContourMapProjection return minValue; } -double RigContourMapCalculator::calculateSum( const RimContourMapProjection& contourMapProjection, +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapCalculator::calculateSum( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ) { @@ -221,14 +225,15 @@ double RigContourMapCalculator::calculateSum( const RimContourMapProjection& /// //-------------------------------------------------------------------------------------------------- std::vector>> - RigContourMapCalculator::generateGridMapping( RimContourMapProjection& contourMapProjection, const RigContourMapGrid& contourMapGrid ) + RigContourMapCalculator::generateGridMapping( RigContourMapProjection& contourMapProjection, + const RigContourMapGrid& contourMapGrid, + ResultAggregationEnum resultAggregation, + const std::vector& weightingResultValues ) { int nCells = contourMapGrid.numberOfCells(); std::vector>> projected3dGridIndices( nCells ); - std::vector weightingResultValues = contourMapProjection.retrieveParameterWeights(); - - if ( contourMapProjection.isStraightSummationResult() ) + if ( RigContourMapCalculator::isStraightSummationResult( resultAggregation ) ) { #pragma omp parallel for for ( int index = 0; index < nCells; ++index ) @@ -260,7 +265,7 @@ std::vector>> /// //-------------------------------------------------------------------------------------------------- std::vector - RigContourMapCalculator::cellOverlapVolumesAndResults( const RimContourMapProjection& contourMapProjection, + RigContourMapCalculator::cellOverlapVolumesAndResults( const RigContourMapProjection& contourMapProjection, const RigContourMapGrid& contourMapGrid, const cvf::Vec2d& globalPos2d, const std::vector& weightingResultValues ) @@ -326,7 +331,7 @@ std::vector /// //-------------------------------------------------------------------------------------------------- std::vector - RigContourMapCalculator::cellRayIntersectionAndResults( const RimContourMapProjection& contourMapProjection, + RigContourMapCalculator::cellRayIntersectionAndResults( const RigContourMapProjection& contourMapProjection, const RigContourMapGrid& contourMapGrid, const cvf::Vec2d& globalPos2d, const std::vector& weightingResultValues ) diff --git a/ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.h b/ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.h index 19257c1f56..199f578ea7 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.h +++ b/ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.h @@ -24,7 +24,7 @@ #include class RigContourMapGrid; -class RimContourMapProjection; +class RigContourMapProjection; //================================================================================================== /// @@ -50,20 +50,22 @@ public: RESULTS_HC_COLUMN }; - static std::vector>> generateGridMapping( RimContourMapProjection& contourMapProjection, - const RigContourMapGrid& contourMapGrid ); + static std::vector>> generateGridMapping( RigContourMapProjection& contourMapProjection, + const RigContourMapGrid& contourMapGrid, + ResultAggregationEnum resultAggregation, + const std::vector& weightingResultValues ); - static double calculateValueInMapCell( const RimContourMapProjection& contourMapProjection, + static double calculateValueInMapCell( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues, ResultAggregationEnum resultAggregation ); - static std::vector cellOverlapVolumesAndResults( const RimContourMapProjection& contourMapProjection, + static std::vector cellOverlapVolumesAndResults( const RigContourMapProjection& contourMapProjection, const RigContourMapGrid& contourMapGrid, const cvf::Vec2d& globalPos2d, const std::vector& weightingResultValues ); - static std::vector cellRayIntersectionAndResults( const RimContourMapProjection& contourMapProjection, + static std::vector cellRayIntersectionAndResults( const RigContourMapProjection& contourMapProjection, const RigContourMapGrid& contourMapGrid, const cvf::Vec2d& globalPos2d, const std::vector& weightingResultValues ); @@ -73,26 +75,26 @@ public: static bool isStraightSummationResult( ResultAggregationEnum aggregationType ); private: - static double calculateTopValue( const RimContourMapProjection& contourMapProjection, + static double calculateTopValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ); - static double calculateMeanValue( const RimContourMapProjection& contourMapProjection, + static double calculateMeanValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ); - static double calculateGeometricMeanValue( const RimContourMapProjection& contourMapProjection, + static double calculateGeometricMeanValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ); - static double calculateHarmonicMeanValue( const RimContourMapProjection& contourMapProjection, + static double calculateHarmonicMeanValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ); - static double calculateMaxValue( const RimContourMapProjection& contourMapProjection, + static double calculateMaxValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ); - static double calculateMinValue( const RimContourMapProjection& contourMapProjection, + static double calculateMinValue( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ); - static double calculateSum( const RimContourMapProjection& contourMapProjection, + static double calculateSum( const RigContourMapProjection& contourMapProjection, const std::vector>& matchingCells, const std::vector& gridCellValues ); }; diff --git a/ApplicationLibCode/ReservoirDataModel/RigContourMapProjection.cpp b/ApplicationLibCode/ReservoirDataModel/RigContourMapProjection.cpp new file mode 100644 index 0000000000..809dd47529 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigContourMapProjection.cpp @@ -0,0 +1,467 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2018- 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 "RigContourMapProjection.h" + +#include "RiaWeightedMeanCalculator.h" + +#include "RigContourMapCalculator.h" +#include "RigContourMapGrid.h" + +#include "cvfArray.h" +#include "cvfGeometryTools.h" + +#include +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigContourMapProjection::RigContourMapProjection( const RigContourMapGrid& contourMapGrid ) + : m_contourMapGrid( contourMapGrid ) + , m_currentResultTimestep( -1 ) +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector>> + RigContourMapProjection::generateGridMapping( RigContourMapCalculator::ResultAggregationEnum resultAggregation, + const std::vector& weights ) +{ + m_projected3dGridIndices = RigContourMapCalculator::generateGridMapping( *this, m_contourMapGrid, resultAggregation, weights ); + return m_projected3dGridIndices; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::maxValue() const +{ + return maxValue( m_aggregatedResults ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::minValue() const +{ + return minValue( m_aggregatedResults ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::meanValue() const +{ + return sumAllValues() / numberOfValidCells(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::sumAllValues() const +{ + double sum = 0.0; + + for ( size_t index = 0; index < m_aggregatedResults.size(); ++index ) + { + if ( m_aggregatedResults[index] != std::numeric_limits::infinity() ) + { + sum += m_aggregatedResults[index]; + } + } + return sum; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2ui RigContourMapProjection::numberOfElementsIJ() const +{ + return m_contourMapGrid.numberOfElementsIJ(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2ui RigContourMapProjection::numberOfVerticesIJ() const +{ + return m_contourMapGrid.numberOfVerticesIJ(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::valueAtVertex( unsigned int i, unsigned int j ) const +{ + size_t index = m_contourMapGrid.vertexIndexFromIJ( i, j ); + if ( index < numberOfVertices() ) + { + return m_aggregatedVertexResults.at( index ); + } + return std::numeric_limits::infinity(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +unsigned int RigContourMapProjection::numberOfCells() const +{ + return m_contourMapGrid.numberOfCells(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +unsigned int RigContourMapProjection::numberOfValidCells() const +{ + unsigned int validCount = 0u; + for ( unsigned int i = 0; i < numberOfCells(); ++i ) + { + cvf::Vec2ui ij = m_contourMapGrid.ijFromCellIndex( i ); + if ( hasResultInCell( ij.x(), ij.y() ) ) + { + validCount++; + } + } + return validCount; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigContourMapProjection::numberOfVertices() const +{ + cvf::Vec2ui gridSize = numberOfVerticesIJ(); + return static_cast( gridSize.x() ) * static_cast( gridSize.y() ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RigContourMapProjection::checkForMapIntersection( const cvf::Vec3d& domainPoint3d, cvf::Vec2d* contourMapPoint, double* valueAtPoint ) const +{ + CVF_TIGHT_ASSERT( contourMapPoint ); + CVF_TIGHT_ASSERT( valueAtPoint ); + + const cvf::Vec3d& minPoint = m_contourMapGrid.expandedBoundingBox().min(); + cvf::Vec3d mapPos3d = domainPoint3d - minPoint; + cvf::Vec2d mapPos2d( mapPos3d.x(), mapPos3d.y() ); + cvf::Vec2d gridorigin( minPoint.x(), minPoint.y() ); + + double value = interpolateValue( mapPos2d ); + if ( value != std::numeric_limits::infinity() ) + { + *valueAtPoint = value; + *contourMapPoint = mapPos2d + gridorigin; + + return true; + } + return false; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec3d RigContourMapProjection::origin3d() const +{ + return m_contourMapGrid.expandedBoundingBox().min(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigContourMapProjection::gridResultIndex( size_t globalCellIdx ) const +{ + return globalCellIdx; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::calculateValueInMapCell( unsigned int i, + unsigned int j, + const std::vector& gridCellValues, + RigContourMapCalculator::ResultAggregationEnum resultAggregation ) const +{ + const std::vector>& matchingCells = cellsAtIJ( i, j ); + return RigContourMapCalculator::calculateValueInMapCell( *this, matchingCells, gridCellValues, resultAggregation ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::maxValue( const std::vector& aggregatedResults ) +{ + double maxV = -std::numeric_limits::infinity(); + + for ( size_t index = 0; index < aggregatedResults.size(); ++index ) + { + if ( aggregatedResults[index] != std::numeric_limits::infinity() ) + { + maxV = std::max( maxV, aggregatedResults[index] ); + } + } + return maxV; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::minValue( const std::vector& aggregatedResults ) +{ + double minV = std::numeric_limits::infinity(); + + for ( size_t index = 0; index < aggregatedResults.size(); ++index ) + { + if ( aggregatedResults[index] != std::numeric_limits::infinity() ) + { + minV = std::min( minV, aggregatedResults[index] ); + } + } + return minV; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigContourMapProjection::setCellVisibility( cvf::ref cellVisibility ) +{ + m_cellGridIdxVisibility = cellVisibility; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::ref RigContourMapProjection::getCellVisibility() const +{ + return m_cellGridIdxVisibility; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigContourMapProjection::generateVertexResults() +{ + size_t nCells = numberOfCells(); + if ( nCells != m_aggregatedResults.size() ) return; + + size_t nVertices = numberOfVertices(); + m_aggregatedVertexResults = std::vector( nVertices, std::numeric_limits::infinity() ); +#pragma omp parallel for + for ( int index = 0; index < static_cast( nVertices ); ++index ) + { + cvf::Vec2ui ij = m_contourMapGrid.ijFromVertexIndex( index ); + m_aggregatedVertexResults[index] = calculateValueAtVertex( ij.x(), ij.y() ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::sumTriangleAreas( const std::vector& triangles ) +{ + double sumArea = 0.0; + for ( size_t i = 0; i < triangles.size(); i += 3 ) + { + cvf::Vec3d v1( triangles[i].x(), triangles[i].y(), triangles[i].z() ); + cvf::Vec3d v2( triangles[i + 1].x(), triangles[i + 1].y(), triangles[i + 1].z() ); + cvf::Vec3d v3( triangles[i + 2].x(), triangles[i + 2].y(), triangles[i + 2].z() ); + double area = 0.5 * ( ( v3 - v1 ) ^ ( v2 - v1 ) ).length(); + sumArea += area; + } + return sumArea; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::interpolateValue( const cvf::Vec2d& gridPos2d ) const +{ + cvf::Vec2ui cellContainingPoint = m_contourMapGrid.ijFromLocalPos( gridPos2d ); + cvf::Vec2d cellCenter = m_contourMapGrid.cellCenterPosition( cellContainingPoint.x(), cellContainingPoint.y() ); + + double sampleSpacing = m_contourMapGrid.sampleSpacing(); + std::array x; + x[0] = cvf::Vec3d( cellCenter + cvf::Vec2d( -sampleSpacing * 0.5, -sampleSpacing * 0.5 ), 0.0 ); + x[1] = cvf::Vec3d( cellCenter + cvf::Vec2d( sampleSpacing * 0.5, -sampleSpacing * 0.5 ), 0.0 ); + x[2] = cvf::Vec3d( cellCenter + cvf::Vec2d( sampleSpacing * 0.5, sampleSpacing * 0.5 ), 0.0 ); + x[3] = cvf::Vec3d( cellCenter + cvf::Vec2d( -sampleSpacing * 0.5, sampleSpacing * 0.5 ), 0.0 ); + + cvf::Vec4d baryCentricCoords = cvf::GeometryTools::barycentricCoords( x[0], x[1], x[2], x[3], cvf::Vec3d( gridPos2d, 0.0 ) ); + + std::array v; + v[0] = cellContainingPoint; + v[1] = cvf::Vec2ui( cellContainingPoint.x() + 1u, cellContainingPoint.y() ); + v[2] = cvf::Vec2ui( cellContainingPoint.x() + 1u, cellContainingPoint.y() + 1u ); + v[3] = cvf::Vec2ui( cellContainingPoint.x(), cellContainingPoint.y() + 1u ); + + std::array vertexValues; + double validBarycentricCoordsSum = 0.0; + for ( int i = 0; i < 4; ++i ) + { + double vertexValue = valueAtVertex( v[i].x(), v[i].y() ); + if ( vertexValue == std::numeric_limits::infinity() ) + { + return std::numeric_limits::infinity(); + } + else + { + vertexValues[i] = vertexValue; + validBarycentricCoordsSum += baryCentricCoords[i]; + } + } + + if ( validBarycentricCoordsSum < 1.0e-8 ) + { + return std::numeric_limits::infinity(); + } + + // Calculate final value + double value = 0.0; + for ( int i = 0; i < 4; ++i ) + { + value += baryCentricCoords[i] / validBarycentricCoordsSum * vertexValues[i]; + } + + return value; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::valueInCell( unsigned int i, unsigned int j ) const +{ + size_t index = m_contourMapGrid.cellIndexFromIJ( i, j ); + if ( index < numberOfCells() ) + { + return m_aggregatedResults.at( index ); + } + return std::numeric_limits::infinity(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RigContourMapProjection::hasResultInCell( unsigned int i, unsigned int j ) const +{ + return !cellsAtIJ( i, j ).empty(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapProjection::calculateValueAtVertex( unsigned int vi, unsigned int vj ) const +{ + std::vector averageIs; + std::vector averageJs; + + if ( vi > 0u ) averageIs.push_back( vi - 1 ); + if ( vj > 0u ) averageJs.push_back( vj - 1 ); + if ( vi < m_contourMapGrid.mapSize().x() ) averageIs.push_back( vi ); + if ( vj < m_contourMapGrid.mapSize().y() ) averageJs.push_back( vj ); + + RiaWeightedMeanCalculator calc; + for ( unsigned int j : averageJs ) + { + for ( unsigned int i : averageIs ) + { + if ( hasResultInCell( i, j ) ) + { + calc.addValueAndWeight( valueInCell( i, j ), 1.0 ); + } + } + } + if ( calc.validAggregatedWeight() ) + { + return calc.weightedMean(); + } + return std::numeric_limits::infinity(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector> RigContourMapProjection::cellsAtIJ( unsigned int i, unsigned int j ) const +{ + size_t cellIndex = m_contourMapGrid.cellIndexFromIJ( i, j ); + if ( cellIndex < m_projected3dGridIndices.size() ) + { + return m_projected3dGridIndices[cellIndex]; + } + return std::vector>(); +} + +//-------------------------------------------------------------------------------------------------- +/// Vertex positions in local coordinates (add origin2d.x() for UTM x) +//-------------------------------------------------------------------------------------------------- +std::vector RigContourMapProjection::xVertexPositions() const +{ + return m_contourMapGrid.xVertexPositions(); +} + +//-------------------------------------------------------------------------------------------------- +/// Vertex positions in local coordinates (add origin2d.y() for UTM y) +//-------------------------------------------------------------------------------------------------- +std::vector RigContourMapProjection::yVertexPositions() const +{ + return m_contourMapGrid.yVertexPositions(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const std::vector& RigContourMapProjection::aggregatedResults() const +{ + return m_aggregatedResults; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const std::vector& RigContourMapProjection::aggregatedVertexResults() const +{ + return m_aggregatedVertexResults; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const std::vector>>& RigContourMapProjection::projected3dGridIndices() const +{ + return m_projected3dGridIndices; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigContourMapProjection::clearResults() +{ + m_aggregatedResults.clear(); + m_aggregatedVertexResults.clear(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigContourMapProjection::clearGridMapping() +{ + m_projected3dGridIndices.clear(); +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigContourMapProjection.h b/ApplicationLibCode/ReservoirDataModel/RigContourMapProjection.h new file mode 100644 index 0000000000..0b68ac2e9e --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigContourMapProjection.h @@ -0,0 +1,119 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2018- 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 "RigContourMapCalculator.h" +#include "RigContourPolygonsTools.h" + +#include "cvfArray.h" +#include "cvfBoundingBox.h" +#include "cvfVector2.h" + +class RigContourMapGrid; +class RimGridView; +class RimRegularLegendConfig; + +//================================================================================================== +/// +/// +//================================================================================================== +class RigContourMapProjection +{ +public: + using CellIndexAndResult = std::pair; + + RigContourMapProjection( const RigContourMapGrid& ); + + void clearResults(); + void clearGridMapping(); + + void generateVertexResults(); + + double maxValue() const; + double minValue() const; + + double meanValue() const; + double sumAllValues() const; + + cvf::Vec2ui numberOfElementsIJ() const; + cvf::Vec2ui numberOfVerticesIJ() const; + + double valueAtVertex( unsigned int i, unsigned int j ) const; + + unsigned int numberOfCells() const; + unsigned int numberOfValidCells() const; + size_t numberOfVertices() const; + + bool checkForMapIntersection( const cvf::Vec3d& domainPoint3d, cvf::Vec2d* contourMapPoint, double* valueAtPoint ) const; + cvf::Vec3d origin3d() const; + + std::vector xVertexPositions() const; + std::vector yVertexPositions() const; + + // Use this function to get the result index into grid cell results. The index will differ if we have active cells + virtual size_t gridResultIndex( size_t globalCellIdx ) const; + + virtual size_t kLayer( size_t globalCellIdx ) const = 0; + virtual size_t kLayers() const = 0; + virtual std::vector findIntersectingCells( const cvf::BoundingBox& bbox ) const = 0; + virtual double calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const = 0; + virtual double calculateRayLengthInCell( size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint ) const = 0; + virtual double getParameterWeightForCell( size_t globalCellIdx, const std::vector& parameterWeights ) const = 0; + virtual std::vector getMapCellVisibility( int viewStepIndex, RigContourMapCalculator::ResultAggregationEnum resultAggregation ) = 0; + + void setCellVisibility( cvf::ref cellVisibility ); + cvf::ref getCellVisibility() const; + + std::vector>> + generateGridMapping( RigContourMapCalculator::ResultAggregationEnum resultAggregation, const std::vector& weights ); + + double interpolateValue( const cvf::Vec2d& gridPosition2d ) const; + + const std::vector& aggregatedResults() const; + const std::vector& aggregatedVertexResults() const; + const std::vector>>& projected3dGridIndices() const; + + // Cell index and position conversion + std::vector cellsAtIJ( unsigned int i, unsigned int j ) const; + + static double maxValue( const std::vector& aggregatedResults ); + static double minValue( const std::vector& aggregatedResults ); + static double sumTriangleAreas( const std::vector& triangles ); + +protected: + double calculateValueInMapCell( unsigned int i, + unsigned int j, + const std::vector& gridCellValues, + RigContourMapCalculator::ResultAggregationEnum resultAggregation ) const; + + double valueInCell( unsigned int i, unsigned int j ) const; + bool hasResultInCell( unsigned int i, unsigned int j ) const; + double calculateValueAtVertex( unsigned int i, unsigned int j ) const; + +protected: + cvf::ref m_cellGridIdxVisibility; + std::vector m_aggregatedResults; + std::vector m_aggregatedVertexResults; + std::vector>> m_projected3dGridIndices; + + int m_currentResultTimestep; + std::vector m_mapCellVisibility; + + const RigContourMapGrid& m_contourMapGrid; +}; diff --git a/ApplicationLibCode/ReservoirDataModel/RigEclipseContourMapProjection.cpp b/ApplicationLibCode/ReservoirDataModel/RigEclipseContourMapProjection.cpp new file mode 100644 index 0000000000..cee5e82926 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigEclipseContourMapProjection.cpp @@ -0,0 +1,334 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2018- 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 "RigEclipseContourMapProjection.h" + +#include "RigActiveCellInfo.h" +#include "RigCaseCellResultsData.h" +#include "RigCell.h" +#include "RigCellGeometryTools.h" +#include "RigContourMapCalculator.h" +#include "RigContourMapGrid.h" +#include "RigEclipseCaseData.h" +#include "RigEclipseResultAddress.h" +#include "RigHexIntersectionTools.h" +#include "RigMainGrid.h" + +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigEclipseContourMapProjection::RigEclipseContourMapProjection( const RigContourMapGrid& contourMapGrid, + RigEclipseCaseData& eclipseCaseData, + RigCaseCellResultsData& resultData ) + : RigContourMapProjection( contourMapGrid ) + , m_eclipseCaseData( eclipseCaseData ) + , m_resultData( resultData ) + , m_kLayers( 0u ) + , m_useActiveCellInfo( true ) +{ + m_mainGrid = m_eclipseCaseData.mainGrid(); + m_activeCellInfo = m_eclipseCaseData.activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL ); + m_kLayers = m_mainGrid->cellCountK(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigEclipseContourMapProjection::~RigEclipseContourMapProjection() +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigEclipseContourMapProjection::generateAndSaveResults( const RigEclipseResultAddress& resultAddress, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int timeStep ) +{ + std::tie( m_useActiveCellInfo, m_aggregatedResults ) = + generateResults( *this, m_contourMapGrid, m_resultData, resultAddress, resultAggregation, timeStep ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigEclipseContourMapProjection::generateResults( const RigEclipseResultAddress& resultAddress, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int timeStep ) const +{ + std::pair> result = + generateResults( *this, m_contourMapGrid, m_resultData, resultAddress, resultAggregation, timeStep ); + return result.second; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair> + RigEclipseContourMapProjection::generateResults( const RigEclipseContourMapProjection& contourMapProjection, + const RigContourMapGrid& contourMapGrid, + RigCaseCellResultsData& resultData, + const RigEclipseResultAddress& resultAddress, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int timeStep ) +{ + size_t nCells = contourMapProjection.numberOfCells(); + + std::vector aggregatedResults = std::vector( nCells, std::numeric_limits::infinity() ); + + bool useActiveCellInfo = resultAddress.isValid() && resultData.hasResultEntry( resultAddress ) && + resultData.isUsingGlobalActiveIndex( resultAddress ); + + auto isTernaryResult = []( const RigEclipseResultAddress& address ) -> bool + { + return address.resultCatType() == RiaDefines::ResultCatType::DYNAMIC_NATIVE && + ( address.resultName().compare( RiaResultNames::ternarySaturationResultName(), Qt::CaseInsensitive ) == 0 ); + }; + + if ( !isTernaryResult( resultAddress ) ) + { + std::vector gridResultValues; + if ( RigContourMapCalculator::isColumnResult( resultAggregation ) ) + { + resultData.ensureKnownResultLoaded( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PORO" ) ); + resultData.ensureKnownResultLoaded( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "NTG" ) ); + resultData.ensureKnownResultLoaded( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "DZ" ) ); + if ( resultAggregation == RigContourMapCalculator::RESULTS_OIL_COLUMN || + resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN ) + { + resultData.ensureKnownResultLoaded( + RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::soil() ) ); + } + if ( resultAggregation == RigContourMapCalculator::RESULTS_GAS_COLUMN || + resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN ) + { + resultData.ensureKnownResultLoaded( + RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() ) ); + } + gridResultValues = calculateColumnResult( resultData, resultAggregation, timeStep ); + } + else + { + // TODO: this was RimEclipseCellColors->hasStaticResult() + if ( resultAddress.resultCatType() == RiaDefines::ResultCatType::STATIC_NATIVE && timeStep > 0 ) timeStep = 0; + + // When loading a project file, grid calculator results are not computed the first time this function is + // called. Must check if result is loaded. See RimReloadCaseTools::updateAll3dViews() + if ( resultAddress.isValid() && resultData.hasResultEntry( resultAddress ) && resultData.isResultLoaded( resultAddress ) ) + { + gridResultValues = resultData.cellScalarResults( resultAddress, timeStep ); + } + } + + if ( !gridResultValues.empty() ) + { +#pragma omp parallel for + for ( int index = 0; index < static_cast( nCells ); ++index ) + { + cvf::Vec2ui ij = contourMapGrid.ijFromCellIndex( index ); + const std::vector>& matchingCells = contourMapProjection.cellsAtIJ( ij.x(), ij.y() ); + aggregatedResults[index] = + RigContourMapCalculator::calculateValueInMapCell( contourMapProjection, matchingCells, gridResultValues, resultAggregation ); + } + } + } + + return { useActiveCellInfo, aggregatedResults }; +} + +std::vector RigEclipseContourMapProjection::calculateColumnResult( RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int timeStep ) const +{ + return calculateColumnResult( m_resultData, resultAggregation, timeStep ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigEclipseContourMapProjection::calculateColumnResult( RigCaseCellResultsData& resultData, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int timeStep ) +{ + bool hasPoroResult = resultData.hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PORO" ) ); + bool hasNtgResult = resultData.hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "NTG" ) ); + bool hasDzResult = resultData.hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "DZ" ) ); + + if ( !( hasPoroResult && hasNtgResult && hasDzResult ) ) + { + return std::vector(); + } + + const std::vector& poroResults = + resultData.cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PORO" ), 0 ); + const std::vector& ntgResults = + resultData.cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "NTG" ), 0 ); + const std::vector& dzResults = + resultData.cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "DZ" ), 0 ); + + CVF_ASSERT( poroResults.size() == ntgResults.size() && ntgResults.size() == dzResults.size() ); + + std::vector resultValues( poroResults.size(), 0.0 ); + + if ( resultAggregation == RigContourMapCalculator::RESULTS_OIL_COLUMN || resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN ) + { + const std::vector& soilResults = + resultData.cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::soil() ), + timeStep ); + for ( size_t cellResultIdx = 0; cellResultIdx < resultValues.size(); ++cellResultIdx ) + { + resultValues[cellResultIdx] = soilResults[cellResultIdx]; + } + } + + if ( resultAggregation == RigContourMapCalculator::RESULTS_GAS_COLUMN || resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN ) + { + bool hasGasResult = + resultData.hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() ) ); + if ( hasGasResult ) + { + const std::vector& sgasResults = + resultData.cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() ), + timeStep ); + for ( size_t cellResultIdx = 0; cellResultIdx < resultValues.size(); ++cellResultIdx ) + { + resultValues[cellResultIdx] += sgasResults[cellResultIdx]; + } + } + } + + for ( size_t cellResultIdx = 0; cellResultIdx < resultValues.size(); ++cellResultIdx ) + { + resultValues[cellResultIdx] *= poroResults[cellResultIdx] * ntgResults[cellResultIdx] * dzResults[cellResultIdx]; + } + + return resultValues; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigEclipseContourMapProjection::findIntersectingCells( const cvf::BoundingBox& bbox ) const +{ + return m_mainGrid->findIntersectingCells( bbox ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigEclipseContourMapProjection::kLayer( size_t globalCellIdx ) const +{ + const RigCell& cell = m_mainGrid->cell( globalCellIdx ); + size_t mainGridCellIdx = cell.mainGridCellIndex(); + size_t i, j, k; + m_mainGrid->ijkFromCellIndex( mainGridCellIdx, &i, &j, &k ); + return k; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigEclipseContourMapProjection::kLayers() const +{ + return m_kLayers; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigEclipseContourMapProjection::calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const +{ + std::array hexCorners; + + const RigCell& cell = m_mainGrid->cell( globalCellIdx ); + + size_t localCellIdx = cell.gridLocalCellIndex(); + RigGridBase* localGrid = cell.hostGrid(); + + localGrid->cellCornerVertices( localCellIdx, hexCorners.data() ); + + cvf::BoundingBox overlapBBox; + std::array overlapCorners; + if ( RigCellGeometryTools::estimateHexOverlapWithBoundingBox( hexCorners, bbox, &overlapCorners, &overlapBBox ) ) + { + double overlapVolume = RigCellGeometryTools::calculateCellVolume( overlapCorners ); + return overlapVolume; + } + return 0.0; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigEclipseContourMapProjection::calculateRayLengthInCell( size_t globalCellIdx, + const cvf::Vec3d& highestPoint, + const cvf::Vec3d& lowestPoint ) const +{ + std::array hexCorners; + + const RigCell& cell = m_mainGrid->cell( globalCellIdx ); + + size_t localCellIdx = cell.gridLocalCellIndex(); + RigGridBase* localGrid = cell.hostGrid(); + + localGrid->cellCornerVertices( localCellIdx, hexCorners.data() ); + std::vector intersections; + + if ( RigHexIntersectionTools::lineHexCellIntersection( highestPoint, lowestPoint, hexCorners.data(), 0, &intersections ) ) + { + double lengthInCell = ( intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint ).length(); + return lengthInCell; + } + return 0.0; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigEclipseContourMapProjection::getParameterWeightForCell( size_t cellResultIdx, const std::vector& cellWeights ) const +{ + if ( cellWeights.empty() ) return 1.0; + + double result = std::max( cellWeights[cellResultIdx], 0.0 ); + if ( result < 1.0e-6 ) + { + result = 0.0; + } + return result; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigEclipseContourMapProjection::gridResultIndex( size_t globalCellIdx ) const +{ + if ( m_useActiveCellInfo ) return m_activeCellInfo->cellResultIndex( globalCellIdx ); + + return globalCellIdx; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigEclipseContourMapProjection::getMapCellVisibility( int viewStepIndex, + RigContourMapCalculator::ResultAggregationEnum resultAggregation ) + +{ + return std::vector( numberOfCells(), true ); +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigEclipseContourMapProjection.h b/ApplicationLibCode/ReservoirDataModel/RigEclipseContourMapProjection.h new file mode 100644 index 0000000000..18cc04c141 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigEclipseContourMapProjection.h @@ -0,0 +1,88 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2018- 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 "RigContourMapCalculator.h" +#include "RigContourMapProjection.h" + +#include "cvfBoundingBox.h" + +class RigActiveCellInfo; +class RigMainGrid; +class RigContourMapGrid; +class RigResultAccessor; +class RigEclipseCaseData; +class RigCaseCellResultsData; +class RigEclipseResultAddress; + +//================================================================================================== +/// +/// +//================================================================================================== +class RigEclipseContourMapProjection : public RigContourMapProjection +{ +public: + RigEclipseContourMapProjection( const RigContourMapGrid& contourMapGrid, + RigEclipseCaseData& eclipseCaseData, + RigCaseCellResultsData& resultData ); + virtual ~RigEclipseContourMapProjection(); + + void generateAndSaveResults( const RigEclipseResultAddress& resultAddress, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int timeStep ); + + std::vector generateResults( const RigEclipseResultAddress& resultAddress, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int timeStep ) const; + + static std::pair> generateResults( const RigEclipseContourMapProjection& contourMapProjection, + const RigContourMapGrid& contourMapGrid, + RigCaseCellResultsData& resultData, + const RigEclipseResultAddress& resultAddress, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int timeStep ); + + std::vector getMapCellVisibility( int viewStepIndex, RigContourMapCalculator::ResultAggregationEnum resultAggregation ) override; + +protected: + using CellIndexAndResult = RigContourMapProjection::CellIndexAndResult; + + std::vector findIntersectingCells( const cvf::BoundingBox& bbox ) const override; + size_t kLayer( size_t globalCellIdx ) const override; + size_t kLayers() const override; + double calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const override; + double calculateRayLengthInCell( size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint ) const override; + double getParameterWeightForCell( size_t cellResultIdx, const std::vector& parameterWeights ) const override; + size_t gridResultIndex( size_t globalCellIdx ) const override; + + // Eclipse implementation specific data generation methods + std::vector calculateColumnResult( RigContourMapCalculator::ResultAggregationEnum resultAggregation, int timeStep ) const; + + static std::vector calculateColumnResult( RigCaseCellResultsData& resultData, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int timeStep ); + +protected: + RigEclipseCaseData& m_eclipseCaseData; + RigCaseCellResultsData& m_resultData; + cvf::ref m_mainGrid; + cvf::ref m_activeCellInfo; + size_t m_kLayers; + bool m_useActiveCellInfo; +}; diff --git a/ApplicationLibCode/ReservoirDataModel/RigGeoMechContourMapProjection.cpp b/ApplicationLibCode/ReservoirDataModel/RigGeoMechContourMapProjection.cpp new file mode 100644 index 0000000000..58c9fe6108 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigGeoMechContourMapProjection.cpp @@ -0,0 +1,381 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2018- 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 "RigGeoMechContourMapProjection.h" + +#include "RiaImageTools.h" +#include "RiaWeightedMeanCalculator.h" + +#include "RigCellGeometryTools.h" +#include "RigContourMapCalculator.h" +#include "RigContourMapGrid.h" +#include "RigFemAddressDefines.h" +#include "RigFemPart.h" +#include "RigFemPartCollection.h" +#include "RigFemPartGrid.h" +#include "RigFemPartResultsCollection.h" +#include "RigGeoMechCaseData.h" +#include "RigHexIntersectionTools.h" + +#include "RivFemElmVisibilityCalculator.h" + +#include "cvfVector3.h" + +#include +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigGeoMechContourMapProjection::RigGeoMechContourMapProjection( RigGeoMechCaseData& caseData, + const RigContourMapGrid& contourMapGrid, + bool limitToPorePressureRegions, + double paddingAroundPorePressureRegion ) + : RigContourMapProjection( contourMapGrid ) + , m_caseData( caseData ) + , m_limitToPorePressureRegions( limitToPorePressureRegions ) + , m_paddingAroundPorePressureRegion( paddingAroundPorePressureRegion ) + , m_kLayers( 0u ) +{ + m_femPart = m_caseData.femParts()->part( 0 ); + m_femPartGrid = m_femPart->getOrCreateStructGrid(); + m_kLayers = m_femPartGrid->cellCountK(); + m_femPart->ensureIntersectionSearchTreeIsBuilt(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::BoundingBox RigGeoMechContourMapProjection::calculateExpandedPorBarBBox( RigGeoMechCaseData& caseData, + const std::string& resultComponentName, + int timeStep, + int frameIndex, + double paddingAroundPorePressureRegion ) +{ + RigFemResultAddress porBarAddr( RigFemResultPosEnum::RIG_ELEMENT_NODAL, RigFemAddressDefines::porBar(), resultComponentName ); + RigFemPartResultsCollection* resultCollection = caseData.femPartResults(); + + const std::vector& resultValues = resultCollection->resultValues( porBarAddr, 0, timeStep, frameIndex ); + cvf::BoundingBox boundingBox; + + if ( resultValues.empty() ) + { + return boundingBox; + } + + auto femPart = caseData.femParts()->part( 0 ); + auto femPartGrid = femPart->getOrCreateStructGrid(); + for ( int i = 0; i < femPart->elementCount(); ++i ) + { + size_t resValueIdx = femPart->elementNodeResultIdx( (int)i, 0 ); + CVF_ASSERT( resValueIdx < resultValues.size() ); + double scalarValue = resultValues[resValueIdx]; + bool validPorValue = scalarValue != std::numeric_limits::infinity(); + + if ( validPorValue ) + { + std::array hexCorners; + femPartGrid->cellCornerVertices( i, hexCorners.data() ); + for ( size_t c = 0; c < 8; ++c ) + { + boundingBox.add( hexCorners[c] ); + } + } + } + + cvf::Vec3d boxMin = boundingBox.min(); + cvf::Vec3d boxMax = boundingBox.max(); + cvf::Vec3d boxExtent = boundingBox.extent(); + boxMin.x() -= boxExtent.x() * 0.5 * paddingAroundPorePressureRegion; + boxMin.y() -= boxExtent.y() * 0.5 * paddingAroundPorePressureRegion; + boxMax.x() += boxExtent.x() * 0.5 * paddingAroundPorePressureRegion; + boxMax.y() += boxExtent.y() * 0.5 * paddingAroundPorePressureRegion; + return cvf::BoundingBox( boxMin, boxMax ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigGeoMechContourMapProjection::getMapCellVisibility( int viewStepIndex, + RigContourMapCalculator::ResultAggregationEnum resultAggregation ) +{ + return getMapCellVisibility( m_currentResultAddr, viewStepIndex, resultAggregation ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigGeoMechContourMapProjection::getMapCellVisibility( RigFemResultAddress resAddr, + int viewStepIndex, + RigContourMapCalculator::ResultAggregationEnum resultAggregation ) +{ + cvf::Vec2ui nCellsIJ = numberOfElementsIJ(); + std::vector> distanceImage( nCellsIJ.x(), std::vector( nCellsIJ.y(), 0u ) ); + + std::vector mapCellVisibility; + + if ( m_limitToPorePressureRegions ) + { + resAddr = RigFemAddressDefines::elementNodalPorBarAddress(); + } + + std::vector cellResults = generateResultsFromAddress( resAddr, mapCellVisibility, resultAggregation, viewStepIndex ); + + mapCellVisibility.resize( numberOfCells(), true ); + CVF_ASSERT( mapCellVisibility.size() == cellResults.size() ); + + { + cvf::BoundingBox validResBoundingBox; + for ( size_t cellIndex = 0; cellIndex < cellResults.size(); ++cellIndex ) + { + cvf::Vec2ui ij = m_contourMapGrid.ijFromCellIndex( cellIndex ); + if ( cellResults[cellIndex] != std::numeric_limits::infinity() ) + { + distanceImage[ij.x()][ij.y()] = 1u; + validResBoundingBox.add( cvf::Vec3d( m_contourMapGrid.cellCenterPosition( ij.x(), ij.y() ), 0.0 ) ); + } + else + { + mapCellVisibility[cellIndex] = false; + } + } + + if ( m_limitToPorePressureRegions && m_paddingAroundPorePressureRegion > 0.0 ) + { + RiaImageTools::distanceTransform2d( distanceImage ); + + cvf::Vec3d porExtent = validResBoundingBox.extent(); + double radius = std::max( porExtent.x(), porExtent.y() ) * 0.25; + double expansion = m_paddingAroundPorePressureRegion * radius; + size_t cellPadding = std::ceil( expansion / m_contourMapGrid.sampleSpacing() ); + for ( size_t cellIndex = 0; cellIndex < cellResults.size(); ++cellIndex ) + { + if ( !mapCellVisibility[cellIndex] ) + { + cvf::Vec2ui ij = m_contourMapGrid.ijFromCellIndex( cellIndex ); + if ( distanceImage[ij.x()][ij.y()] < cellPadding * cellPadding ) + { + mapCellVisibility[cellIndex] = true; + } + } + } + } + } + return mapCellVisibility; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigGeoMechContourMapProjection::generateAndSaveResults( RigFemResultAddress resultAddress, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int viewerStepIndex ) +{ + m_aggregatedResults = generateResultsFromAddress( resultAddress, m_mapCellVisibility, resultAggregation, viewerStepIndex ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigGeoMechContourMapProjection::generateResultsFromAddress( RigFemResultAddress resultAddress, + const std::vector& mapCellVisibility, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int viewerStepIndex ) const +{ + RigFemPartResultsCollection* resultCollection = m_caseData.femPartResults(); + size_t nCells = numberOfCells(); + std::vector aggregatedResults = std::vector( nCells, std::numeric_limits::infinity() ); + + auto [stepIdx, frameIdx] = m_caseData.femPartResults()->stepListIndexToTimeStepAndDataFrameIndex( viewerStepIndex ); + + bool wasInvalid = false; + if ( !resultAddress.isValid() ) + { + wasInvalid = true; + resultAddress = RigFemAddressDefines::elementNodalPorBarAddress(); + } + + if ( resultAddress.fieldName == "PP" ) + { + resultAddress.fieldName = RigFemAddressDefines::porBar(); // More likely to be in memory than POR + } + if ( resultAddress.fieldName == RigFemAddressDefines::porBar() ) + { + resultAddress.resultPosType = RIG_ELEMENT_NODAL; + } + else if ( resultAddress.resultPosType == RIG_FORMATION_NAMES ) + { + resultAddress.resultPosType = RIG_ELEMENT_NODAL; // formation indices are stored per element node result. + } + + std::vector resultValuesF = resultCollection->resultValues( resultAddress, 0, stepIdx, frameIdx ); + if ( resultValuesF.empty() ) return aggregatedResults; + + std::vector resultValues = gridCellValues( resultAddress, resultValuesF ); + + if ( wasInvalid ) + { + // For invalid result addresses we just use the POR-Bar result to get the reservoir region + // And display a dummy 0-result in the region. + for ( double& value : resultValues ) + { + if ( value != std::numeric_limits::infinity() ) + { + value = 0.0; + } + } + } + +#pragma omp parallel for + for ( int index = 0; index < static_cast( nCells ); ++index ) + { + if ( mapCellVisibility.empty() || mapCellVisibility[index] ) + { + cvf::Vec2ui ij = m_contourMapGrid.ijFromCellIndex( index ); + aggregatedResults[index] = calculateValueInMapCell( ij.x(), ij.y(), resultValues, resultAggregation ); + } + } + + return aggregatedResults; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigGeoMechContourMapProjection::findIntersectingCells( const cvf::BoundingBox& bbox ) const +{ + return m_femPart->findIntersectingElementsWithExistingSearchTree( bbox ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigGeoMechContourMapProjection::kLayer( size_t globalCellIdx ) const +{ + size_t i, j, k; + m_femPartGrid->ijkFromCellIndex( globalCellIdx, &i, &j, &k ); + return k; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigGeoMechContourMapProjection::kLayers() const +{ + return m_kLayers; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigGeoMechContourMapProjection::calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const +{ + std::array hexCorners; + m_femPartGrid->cellCornerVertices( globalCellIdx, hexCorners.data() ); + + cvf::BoundingBox overlapBBox; + std::array overlapCorners; + if ( RigCellGeometryTools::estimateHexOverlapWithBoundingBox( hexCorners, bbox, &overlapCorners, &overlapBBox ) ) + { + double overlapVolume = RigCellGeometryTools::calculateCellVolume( overlapCorners ); + return overlapVolume; + } + return 0.0; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigGeoMechContourMapProjection::calculateRayLengthInCell( size_t globalCellIdx, + const cvf::Vec3d& highestPoint, + const cvf::Vec3d& lowestPoint ) const +{ + std::array hexCorners; + + const std::vector& nodeCoords = m_femPart->nodes().coordinates; + const int* cornerIndices = m_femPart->connectivities( globalCellIdx ); + + hexCorners[0] = cvf::Vec3d( nodeCoords[cornerIndices[0]] ); + hexCorners[1] = cvf::Vec3d( nodeCoords[cornerIndices[1]] ); + hexCorners[2] = cvf::Vec3d( nodeCoords[cornerIndices[2]] ); + hexCorners[3] = cvf::Vec3d( nodeCoords[cornerIndices[3]] ); + hexCorners[4] = cvf::Vec3d( nodeCoords[cornerIndices[4]] ); + hexCorners[5] = cvf::Vec3d( nodeCoords[cornerIndices[5]] ); + hexCorners[6] = cvf::Vec3d( nodeCoords[cornerIndices[6]] ); + hexCorners[7] = cvf::Vec3d( nodeCoords[cornerIndices[7]] ); + + std::vector intersections; + + if ( RigHexIntersectionTools::lineHexCellIntersection( highestPoint, lowestPoint, hexCorners.data(), 0, &intersections ) ) + { + double lengthInCell = ( intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint ).length(); + return lengthInCell; + } + return 0.0; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigGeoMechContourMapProjection::getParameterWeightForCell( size_t globalCellIdx, const std::vector& parameterWeights ) const +{ + if ( parameterWeights.empty() ) return 1.0; + + return parameterWeights[globalCellIdx]; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigGeoMechContourMapProjection::gridCellValues( RigFemResultAddress resAddr, std::vector& resultValues ) const +{ + std::vector gridCellValues( m_femPart->elementCount(), std::numeric_limits::infinity() ); + for ( size_t globalCellIdx = 0; globalCellIdx < static_cast( m_femPart->elementCount() ); ++globalCellIdx ) + { + RigElementType elmType = m_femPart->elementType( globalCellIdx ); + if ( !RigFemTypes::is8NodeElement( elmType ) ) continue; + + if ( resAddr.resultPosType == RIG_ELEMENT ) + { + gridCellValues[globalCellIdx] = static_cast( resultValues[globalCellIdx] ); + } + else if ( resAddr.resultPosType == RIG_ELEMENT_NODAL ) + { + RiaWeightedMeanCalculator cellAverage; + for ( int i = 0; i < 8; ++i ) + { + size_t gridResultValueIdx = + m_femPart->resultValueIdxFromResultPosType( resAddr.resultPosType, static_cast( globalCellIdx ), i ); + cellAverage.addValueAndWeight( resultValues[gridResultValueIdx], 1.0 ); + } + + gridCellValues[globalCellIdx] = static_cast( cellAverage.weightedMean() ); + } + else + { + RiaWeightedMeanCalculator cellAverage; + const int* elmNodeIndices = m_femPart->connectivities( globalCellIdx ); + for ( int i = 0; i < 8; ++i ) + { + cellAverage.addValueAndWeight( resultValues[elmNodeIndices[i]], 1.0 ); + } + gridCellValues[globalCellIdx] = static_cast( cellAverage.weightedMean() ); + } + } + return gridCellValues; +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigGeoMechContourMapProjection.h b/ApplicationLibCode/ReservoirDataModel/RigGeoMechContourMapProjection.h new file mode 100644 index 0000000000..bba506e53b --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigGeoMechContourMapProjection.h @@ -0,0 +1,81 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2018- 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 "RigContourMapCalculator.h" +#include "RigContourMapProjection.h" +#include "RigFemPart.h" +#include "RigFemResultAddress.h" + +#include "cvfBoundingBox.h" + +class RigGeoMechCaseData; + +//================================================================================================== +/// +/// +//================================================================================================== +class RigGeoMechContourMapProjection : public RigContourMapProjection +{ +public: + RigGeoMechContourMapProjection( RigGeoMechCaseData& caseData, + const RigContourMapGrid&, + bool limitToPorePressureRegions, + double paddingAroundPorePressureRegion ); + + void generateAndSaveResults( RigFemResultAddress resultAddress, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int viewerStepIndex ); + + std::vector generateResultsFromAddress( RigFemResultAddress resultAddress, + const std::vector& mapCellVisibility, + RigContourMapCalculator::ResultAggregationEnum resultAggregation, + int viewerStepIndex ) const; + + static cvf::BoundingBox calculateExpandedPorBarBBox( RigGeoMechCaseData& caseData, + const std::string& resultComponentName, + int timeStep, + int frameIndex, + double padding ); + + std::vector getMapCellVisibility( int viewStepIndex, RigContourMapCalculator::ResultAggregationEnum resultAggregation ) override; + + std::vector getMapCellVisibility( RigFemResultAddress resAddr, + int viewStepIndex, + RigContourMapCalculator::ResultAggregationEnum resultAggregation ); + +protected: + // GeoMech implementation specific data generation methods + std::vector findIntersectingCells( const cvf::BoundingBox& bbox ) const override; + size_t kLayer( size_t globalCellIdx ) const override; + size_t kLayers() const override; + double calculateOverlapVolume( size_t globalCellIdx, const cvf::BoundingBox& bbox ) const override; + double calculateRayLengthInCell( size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint ) const override; + double getParameterWeightForCell( size_t globalCellIdx, const std::vector& parameterWeights ) const override; + std::vector gridCellValues( RigFemResultAddress resAddr, std::vector& resultValues ) const; + +protected: + RigGeoMechCaseData& m_caseData; + bool m_limitToPorePressureRegions; + double m_paddingAroundPorePressureRegion; + cvf::ref m_femPart; + cvf::cref m_femPartGrid; + RigFemResultAddress m_currentResultAddr; + size_t m_kLayers; +};