OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Refactor: extract computation for contour map projection.

Browse Source
This commit is contained in:
Kristian Bendiksen 2024-10-25 15:11:42 +02:00
parent d2cd425ad0
commit f5437e6c72
16 changed files with 1738 additions and 994 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

366
ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.cpp
View File

@ -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 <algorithm>
#include <array>
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<double, double> minmaxValAllTimeSteps = minmaxValuesAllTimeSteps();
@ -174,61 +167,16 @@ cvf::ref<cvf::UByteArray> 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<float>& 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<double>::infinity();
if ( validPorValue )
{
std::array<cvf::Vec3d, 8> 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<RigContourMapGrid>( gridBoundingBox, expandedBoundingBox, sampleSpacing() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<bool> RimGeoMechContourMapProjection::getMapCellVisibility()
{
cvf::Vec2ui nCellsIJ = numberOfElementsIJ();
std::vector<std::vector<unsigned int>> distanceImage( nCellsIJ.x(), std::vector<unsigned int>( nCellsIJ.y(), 0u ) );
std::vector<bool> mapCellVisibility;
RigFemResultAddress resAddr = view()->cellResult()->resultAddress();
if ( m_limitToPorePressureRegions )
{
resAddr = RigFemAddressDefines::elementNodalPorBarAddress();
}
std::vector<double> 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<double>::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<RigContourMapGrid>( gridBoundingBox, expandedBoundingBox, sampleSpacing() );
m_contourMapProjection = std::make_unique<RigGeoMechContourMapProjection>( *geoMechCase->geoMechData(),
*m_contourMapGrid,
m_limitToPorePressureRegions,
m_paddingAroundPorePressureRegion );
}
//--------------------------------------------------------------------------------------------------
@ -334,14 +230,30 @@ std::vector<double> RimGeoMechContourMapProjection::retrieveParameterWeights()
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RimGeoMechContourMapProjection::generateResults( int viewerStepIndex )
std::vector<double> 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<double> 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<RigGeoMechContourMapProjection*>( m_contourMapProjection.get() )
->generateAndSaveResults( resultAddress, m_resultAggregation(), timeStep );
}
}
//--------------------------------------------------------------------------------------------------
@ -349,64 +261,10 @@ std::vector<double> RimGeoMechContourMapProjection::generateResults( int viewerS
//--------------------------------------------------------------------------------------------------
std::vector<double> RimGeoMechContourMapProjection::generateResultsFromAddress( RigFemResultAddress resultAddress,
const std::vector<bool>& mapCellVisibility,
int viewerStepIndex )
int viewerStepIndex ) const
{
RigGeoMechCaseData* caseData = geoMechCase()->geoMechData();
RigFemPartResultsCollection* resultCollection = caseData->femPartResults();
size_t nCells = numberOfCells();
std::vector<double> aggregatedResults = std::vector<double>( nCells, std::numeric_limits<double>::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<float> resultValuesF = resultCollection->resultValues( resultAddress, 0, stepIdx, frameIdx );
if ( resultValuesF.empty() ) return aggregatedResults;
std::vector<double> 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<double>::infinity() )
{
value = 0.0;
}
}
}
#pragma omp parallel for
for ( int index = 0; index < static_cast<int>( 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<RigGeoMechContourMapProjection*>( m_contourMapProjection.get() )
->generateResultsFromAddress( resultAddress, mapCellVisibility, m_resultAggregation(), viewerStepIndex );
}
//--------------------------------------------------------------------------------------------------
@ -436,132 +294,6 @@ RimGridView* RimGeoMechContourMapProjection::baseView() const
return view();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<size_t> 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<cvf::Vec3d, 8> hexCorners;
m_femPartGrid->cellCornerVertices( globalCellIdx, hexCorners.data() );
cvf::BoundingBox overlapBBox;
std::array<cvf::Vec3d, 8> 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<cvf::Vec3d, 8> hexCorners;
const std::vector<cvf::Vec3f>& 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<HexIntersectionInfo> 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<double>& parameterWeights ) const
{
if ( parameterWeights.empty() ) return 1.0;
return parameterWeights[globalCellIdx];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RimGeoMechContourMapProjection::gridCellValues( RigFemResultAddress resAddr, std::vector<float>& resultValues ) const
{
std::vector<double> gridCellValues( m_femPart->elementCount(), std::numeric_limits<double>::infinity() );
for ( size_t globalCellIdx = 0; globalCellIdx < static_cast<size_t>( m_femPart->elementCount() ); ++globalCellIdx )
{
RigElementType elmType = m_femPart->elementType( globalCellIdx );
if ( !RigFemTypes::is8NodeElement( elmType ) ) continue;
if ( resAddr.resultPosType == RIG_ELEMENT )
{
gridCellValues[globalCellIdx] = static_cast<double>( resultValues[globalCellIdx] );
}
else if ( resAddr.resultPosType == RIG_ELEMENT_NODAL )
{
RiaWeightedMeanCalculator<float> cellAverage;
for ( int i = 0; i < 8; ++i )
{
size_t gridResultValueIdx =
m_femPart->resultValueIdxFromResultPosType( resAddr.resultPosType, static_cast<int>( globalCellIdx ), i );
cellAverage.addValueAndWeight( resultValues[gridResultValueIdx], 1.0 );
}
gridCellValues[globalCellIdx] = static_cast<double>( cellAverage.weightedMean() );
}
else
{
RiaWeightedMeanCalculator<float> 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<double>( cellAverage.weightedMean() );
}
}
return gridCellValues;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -672,20 +404,14 @@ std::pair<double, double> 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<double> 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 ) );
}
}
}

31
ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.h
View File

@ -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<cvf::UByteArray> getCellVisibility() const override;
cvf::BoundingBox calculateExpandedPorBarBBox( int timeStep, int frameIndex ) const;
void updateGridInformation() override;
std::vector<bool> getMapCellVisibility() override;
std::vector<double> retrieveParameterWeights() override;
std::vector<double> generateResults( int viewerStepIndex ) override;
std::vector<double> retrieveParameterWeights() override;
std::vector<double> generateResults( int timeStep ) const override;
void generateAndSaveResults( int timeStep ) override;
std::vector<double>
generateResultsFromAddress( RigFemResultAddress resultAddress, const std::vector<bool>& mapCellVisibility, int viewerStepIndex );
bool resultVariableChanged() const override;
void clearResultVariable() override;
RimGridView* baseView() const override;
std::vector<size_t> 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<double>& parameterWeights ) const override;
std::vector<double> gridCellValues( RigFemResultAddress resAddr, std::vector<float>& resultValues ) const;
generateResultsFromAddress( RigFemResultAddress resultAddress, const std::vector<bool>& 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<RigFemPart> m_femPart;
cvf::cref<RigFemPartGrid> m_femPartGrid;
RigFemResultAddress m_currentResultAddr;
size_t m_kLayers;
};

403
ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp
View File

@ -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 <algorithm>
@ -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<double>::infinity() )
, m_maxResultAllTimeSteps( -std::numeric_limits<double>::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<double>::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<double>::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<double>::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<double>::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<double>::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<size_t>( gridSize.x() ) * static_cast<size_t>( 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<double>::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<double>& gridCellValues ) const
{
const std::vector<std::pair<size_t, double>>& matchingCells = cellsAtIJ( i, j );
return RigContourMapCalculator::calculateValueInMapCell( *this, matchingCells, gridCellValues, m_resultAggregation() );
const std::vector<std::pair<size_t, double>>& 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<cvf::UByteArray> 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<double>::infinity() &&
m_maxResultAllTimeSteps != -std::numeric_limits<double>::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<double>::infinity();
m_maxResultAllTimeSteps = -std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::maxValue( const std::vector<double>& aggregatedResults ) const
{
double maxV = -std::numeric_limits<double>::infinity();
for ( size_t index = 0; index < aggregatedResults.size(); ++index )
{
if ( aggregatedResults[index] != std::numeric_limits<double>::infinity() )
{
maxV = std::max( maxV, aggregatedResults[index] );
}
}
return maxV;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::minValue( const std::vector<double>& aggregatedResults ) const
{
double minV = std::numeric_limits<double>::infinity();
for ( size_t index = 0; index < aggregatedResults.size(); ++index )
{
if ( aggregatedResults[index] != std::numeric_limits<double>::infinity() )
{
minV = std::min( minV, aggregatedResults[index] );
}
}
return minV;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<double, double> 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<double> 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<cvf::UByteArray> RimContourMapProjection::getCellVisibility() const
return baseView()->currentTotalCellVisibility();
}
//--------------------------------------------------------------------------------------------------
/// Empty default implementation
//--------------------------------------------------------------------------------------------------
std::vector<bool> RimContourMapProjection::getMapCellVisibility()
{
return std::vector<bool>( numberOfCells(), true );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::mapCellVisibilityNeedsUpdating()
bool RimContourMapProjection::mapCellVisibilityNeedsUpdating( int timestep )
{
std::vector<bool> mapCellVisiblity = getMapCellVisibility();
if ( m_contourMapProjection ) return true;
std::vector<bool> 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<double>( nVertices, std::numeric_limits<double>::infinity() );
#pragma omp parallel for
for ( int index = 0; index < static_cast<int>( 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<std::vector<std::vector<cvf::Vec4d>>> threadTriangles( numberOfThreads );
const std::vector<double>& 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<double>::infinity();
triangle[n] = vertices[vn];
uint vn = ( *faceList )[i + n];
double value = vn < aggregatedVertexResults.size() ? aggregatedVertexResults[vn] : std::numeric_limits<double>::infinity();
triangle[n] = vertices[vn];
triangleWithValues[n] = cvf::Vec4d( vertices[vn], value );
if ( value != std::numeric_limits<double>::infinity() )
{
@ -804,7 +688,7 @@ void RimContourMapProjection::generateTrianglesWithVertexValues()
}
if ( value == std::numeric_limits<double>::infinity() )
{
value = interpolateValue( cvf::Vec2d( localVertex.x(), localVertex.y() ) );
value = m_contourMapProjection->interpolateValue( cvf::Vec2d( localVertex.x(), localVertex.y() ) );
if ( value == std::numeric_limits<double>::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<double> contourLevels;
if ( resultRangeIsValid() && legendConfig()->mappingMode() != RimRegularLegendConfig::MappingType::CATEGORY_INTEGER )
{
const std::vector<double>& aggregatedVertexResults = m_contourMapProjection->aggregatedVertexResults();
legendConfig()->scalarMapper()->majorTickValues( &contourLevels );
int nContourLevels = static_cast<int>( contourLevels.size() );
@ -900,7 +786,7 @@ void RimContourMapProjection::generateContourPolygons()
}
std::vector<caf::ContourLines::ListOfLineSegments> unorderedLineSegmentsPerLevel =
caf::ContourLines::create( m_aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels );
caf::ContourLines::create( aggregatedVertexResults, xVertexPositions(), yVertexPositions(), contourLevels );
contourPolygons = std::vector<ContourPolygons>( 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<cvf::Vec4d>& 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<cvf::Vec3d, 4> 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<cvf::Vec2ui, 4> 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<double, 4> 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<double>::infinity() )
{
return std::numeric_limits<double>::infinity();
}
else
{
vertexValues[i] = vertexValue;
validBarycentricCoordsSum += baryCentricCoords[i];
}
}
if ( validBarycentricCoordsSum < 1.0e-8 )
{
return std::numeric_limits<double>::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<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::hasResultInCell( uint i, uint j ) const
{
return !cellsAtIJ( i, j ).empty();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::calculateValueAtVertex( uint vi, uint vj ) const
{
std::vector<uint> averageIs;
std::vector<uint> 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<double> 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<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<size_t, double>> 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<std::pair<size_t, double>>();
}
//--------------------------------------------------------------------------------------------------
/// Vertex positions in local coordinates (add origin2d.x() for UTM x)
//--------------------------------------------------------------------------------------------------
std::vector<double> RimContourMapProjection::xVertexPositions() const
{
return m_contourMapGrid->xVertexPositions();
if ( m_contourMapGrid ) return m_contourMapGrid->xVertexPositions();
return {};
}
//--------------------------------------------------------------------------------------------------
@ -1112,7 +864,8 @@ std::vector<double> RimContourMapProjection::xVertexPositions() const
//--------------------------------------------------------------------------------------------------
std::vector<double> 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<double>::infinity() &&
m_maxResultAllTimeSteps != -std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::clearTimeStepRange()
{
m_minResultAllTimeSteps = std::numeric_limits<double>::infinity();
m_maxResultAllTimeSteps = -std::numeric_limits<double>::infinity();
}

52
ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h
View File

@ -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<size_t, double>;
using ResultAggregation = caf::AppEnum<RigContourMapCalculator::ResultAggregationEnum>;
using ContourPolygons = std::vector<RigContourPolygonsTools::ContourPolygon>;
@ -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<double> retrieveParameterWeights() = 0;
virtual size_t kLayer( size_t globalCellIdx ) const = 0;
virtual size_t kLayers() const = 0;
virtual std::vector<size_t> 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<double>& parameterWeights ) const = 0;
virtual cvf::ref<cvf::UByteArray> getCellVisibility() const;
protected:
// Protected virtual methods to be overridden by Eclipse and Geo-mechanical contour map implementations
virtual void updateGridInformation() = 0;
virtual std::vector<double> 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<double> 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<double>& 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<double>& aggregatedResults ) const;
double minValue( const std::vector<double>& aggregatedResults ) const;
virtual std::pair<double, double> minmaxValuesAllTimeSteps() = 0;
virtual std::pair<double, double> minmaxValuesAllTimeSteps();
virtual std::vector<bool> getMapCellVisibility();
bool mapCellVisibilityNeedsUpdating();
static std::vector<std::vector<std::pair<size_t, double>>> generateGridMapping( RimContourMapProjection& contourMapProjection,
const RigContourMapGrid& contourMapGrid );
bool mapCellVisibilityNeedsUpdating( int timeStep );
void generateVertexResults();
void generateTrianglesWithVertexValues();
void generateContourPolygons();
static double sumTriangleAreas( const std::vector<cvf::Vec4d>& 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<CellIndexAndResult> cellsAtIJ( uint i, uint j ) const;
double gridEdgeOffset() const;
virtual void updateAfterResultGeneration( int timeStep ) = 0;
@ -180,11 +154,6 @@ protected:
caf::PdmField<bool> m_showContourLabels;
caf::PdmField<bool> m_smoothContourLines;
cvf::ref<cvf::UByteArray> m_cellGridIdxVisibility;
std::vector<double> m_aggregatedResults;
std::vector<double> m_aggregatedVertexResults;
std::vector<std::vector<std::pair<size_t, double>>> m_projected3dGridIndices;
cvf::Vec2d m_pickPoint;
std::vector<ContourPolygons> m_contourPolygons;
std::vector<double> m_contourLevelCumulativeAreas;
@ -195,5 +164,6 @@ protected:
double m_minResultAllTimeSteps;
double m_maxResultAllTimeSteps;
std::unique_ptr<RigContourMapGrid> m_contourMapGrid;
std::unique_ptr<RigContourMapGrid> m_contourMapGrid;
std::unique_ptr<RigContourMapProjection> m_contourMapProjection;
};

310
ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp
View File

@ -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 <algorithm>
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<double> RimEclipseContourMapProjection::generateResults( int timeStep )
std::vector<double> RimEclipseContourMapProjection::generateResults( int timeStep ) const
{
m_weightingResult->loadResult();
size_t nCells = numberOfCells();
std::vector<double> aggregatedResults = std::vector<double>( nCells, std::numeric_limits<double>::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<double> 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<int>( nCells ); ++index )
{
cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( index );
aggregatedResults[index] = calculateValueInMapCell( ij.x(), ij.y(), gridResultValues );
}
}
}
return dynamic_cast<RigEclipseContourMapProjection*>( 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<RigEclipseContourMapProjection*>( m_contourMapProjection.get() )
->generateAndSaveResults( resAddr, m_resultAggregation(), timeStep );
}
return aggregatedResults;
}
//--------------------------------------------------------------------------------------------------
@ -258,80 +222,18 @@ void RimEclipseContourMapProjection::clearResultVariable()
m_currentResultName = "";
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> 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<double>();
}
const std::vector<double>& poroResults =
resultData->cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PORO" ), 0 );
const std::vector<double>& ntgResults =
resultData->cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "NTG" ), 0 );
const std::vector<double>& 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<double> resultValues( poroResults.size(), 0.0 );
if ( resultAggregation == RigContourMapCalculator::RESULTS_OIL_COLUMN || resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN )
{
const std::vector<double>& 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<double>& 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<RigContourMapGrid>( gridBoundingBox, sampleSpacing() );
m_contourMapProjection = std::make_unique<RigEclipseContourMapProjection>( *m_contourMapGrid, *eclipseCaseData, *resultData );
}
//--------------------------------------------------------------------------------------------------
@ -376,108 +278,6 @@ RimGridView* RimEclipseContourMapProjection::baseView() const
return view();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<size_t> 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<cvf::Vec3d, 8> 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<cvf::Vec3d, 8> 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<cvf::Vec3d, 8> hexCorners;
RigCell cell = m_mainGrid->cell( globalCellIdx );
size_t localCellIdx = cell.gridLocalCellIndex();
RigGridBase* localGrid = cell.hostGrid();
localGrid->cellCornerVertices( localCellIdx, hexCorners.data() );
std::vector<HexIntersectionInfo> 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<double>& 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<double, double> RimEclipseContourMapProjection::minmaxValuesAllTimeSteps()
{
if ( !resultRangeIsValid() )
{
clearTimeStepRange();
int timeStepCount = static_cast<int>( eclipseCase()->timeStepStrings().size() );
for ( int i = 0; i < (int)timeStepCount; ++i )
{
std::vector<double> 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 );
}

24
ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.h
View File

@ -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<double> retrieveParameterWeights() override;
std::vector<double> generateResults( int timeStep ) override;
std::vector<double> generateResults( int timeStep ) const override;
void generateAndSaveResults( int timeStep ) override;
bool resultVariableChanged() const override;
void clearResultVariable() override;
RimGridView* baseView() const override;
std::vector<size_t> 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<double>& parameterWeights ) const override;
size_t gridResultIndex( size_t globalCellIdx ) const override;
// Eclipse implementation specific data generation methods
std::vector<double> calculateColumnResult( ResultAggregation resultAggregation ) const;
RimEclipseCase* eclipseCase() const;
RimEclipseContourMapView* view() const;
std::pair<double, double> minmaxValuesAllTimeSteps() override;
void updateAfterResultGeneration( int timeStep ) override;
protected:
@ -91,10 +80,5 @@ protected:
caf::PdmField<bool> m_weightByParameter;
caf::PdmChildField<RimEclipseResultDefinition*> m_weightingResult;
cvf::ref<RigMainGrid> m_mainGrid;
cvf::ref<RigActiveCellInfo> m_activeCellInfo;
size_t m_kLayers;
bool m_useActiveCellInfo;
QString m_currentResultName;
};

1
ApplicationLibCode/ProjectDataModel/RimEclipseContourMapView.cpp
View File

@ -272,6 +272,7 @@ void RimEclipseContourMapView::updateGeometry()
{
m_contourMapProjection->generateResultsIfNecessary( m_currentTimeStep() );
}
onUpdateLegends();
progress.setProgress( 30 );
}

10
ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake
View File

@ -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})

37
ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.cpp
View File

@ -24,15 +24,16 @@
#include "RigContourMapGrid.h"
#include "RimCase.h"
#include "RimContourMapProjection.h"
#include "RigContourMapProjection.h"
#include <algorithm>
#include <cmath>
#include <map>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigContourMapCalculator::calculateValueInMapCell( const RimContourMapProjection& contourMapProjection,
double RigContourMapCalculator::calculateValueInMapCell( const RigContourMapProjection& contourMapProjection,
const std::vector<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& 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<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& 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<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& 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<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& 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<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& 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<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& 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<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& 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<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& gridCellValues )
{
@ -221,14 +225,15 @@ double RigContourMapCalculator::calculateSum( const RimContourMapProjection&
///
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<std::pair<size_t, double>>>
RigContourMapCalculator::generateGridMapping( RimContourMapProjection& contourMapProjection, const RigContourMapGrid& contourMapGrid )
RigContourMapCalculator::generateGridMapping( RigContourMapProjection& contourMapProjection,
const RigContourMapGrid& contourMapGrid,
ResultAggregationEnum resultAggregation,
const std::vector<double>& weightingResultValues )
{
int nCells = contourMapGrid.numberOfCells();
std::vector<std::vector<std::pair<size_t, double>>> projected3dGridIndices( nCells );
std::vector<double> 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<std::pair<size_t, double>>>
///
//--------------------------------------------------------------------------------------------------
std::vector<RigContourMapCalculator::CellIndexAndResult>
RigContourMapCalculator::cellOverlapVolumesAndResults( const RimContourMapProjection& contourMapProjection,
RigContourMapCalculator::cellOverlapVolumesAndResults( const RigContourMapProjection& contourMapProjection,
const RigContourMapGrid& contourMapGrid,
const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues )
@ -326,7 +331,7 @@ std::vector<RigContourMapCalculator::CellIndexAndResult>
///
//--------------------------------------------------------------------------------------------------
std::vector<RigContourMapCalculator::CellIndexAndResult>
RigContourMapCalculator::cellRayIntersectionAndResults( const RimContourMapProjection& contourMapProjection,
RigContourMapCalculator::cellRayIntersectionAndResults( const RigContourMapProjection& contourMapProjection,
const RigContourMapGrid& contourMapGrid,
const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues )

28
ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.h
View File

@ -24,7 +24,7 @@
#include <vector>
class RigContourMapGrid;
class RimContourMapProjection;
class RigContourMapProjection;
//==================================================================================================
///
@ -50,20 +50,22 @@ public:
RESULTS_HC_COLUMN
};
static std::vector<std::vector<std::pair<size_t, double>>> generateGridMapping( RimContourMapProjection& contourMapProjection,
const RigContourMapGrid& contourMapGrid );
static std::vector<std::vector<std::pair<size_t, double>>> generateGridMapping( RigContourMapProjection& contourMapProjection,
const RigContourMapGrid& contourMapGrid,
ResultAggregationEnum resultAggregation,
const std::vector<double>& weightingResultValues );
static double calculateValueInMapCell( const RimContourMapProjection& contourMapProjection,
static double calculateValueInMapCell( const RigContourMapProjection& contourMapProjection,
const std::vector<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& gridCellValues,
ResultAggregationEnum resultAggregation );
static std::vector<CellIndexAndResult> cellOverlapVolumesAndResults( const RimContourMapProjection& contourMapProjection,
static std::vector<CellIndexAndResult> cellOverlapVolumesAndResults( const RigContourMapProjection& contourMapProjection,
const RigContourMapGrid& contourMapGrid,
const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues );
static std::vector<CellIndexAndResult> cellRayIntersectionAndResults( const RimContourMapProjection& contourMapProjection,
static std::vector<CellIndexAndResult> cellRayIntersectionAndResults( const RigContourMapProjection& contourMapProjection,
const RigContourMapGrid& contourMapGrid,
const cvf::Vec2d& globalPos2d,
const std::vector<double>& 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<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& gridCellValues );
static double calculateMeanValue( const RimContourMapProjection& contourMapProjection,
static double calculateMeanValue( const RigContourMapProjection& contourMapProjection,
const std::vector<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& gridCellValues );
static double calculateGeometricMeanValue( const RimContourMapProjection& contourMapProjection,
static double calculateGeometricMeanValue( const RigContourMapProjection& contourMapProjection,
const std::vector<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& gridCellValues );
static double calculateHarmonicMeanValue( const RimContourMapProjection& contourMapProjection,
static double calculateHarmonicMeanValue( const RigContourMapProjection& contourMapProjection,
const std::vector<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& gridCellValues );
static double calculateMaxValue( const RimContourMapProjection& contourMapProjection,
static double calculateMaxValue( const RigContourMapProjection& contourMapProjection,
const std::vector<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& gridCellValues );
static double calculateMinValue( const RimContourMapProjection& contourMapProjection,
static double calculateMinValue( const RigContourMapProjection& contourMapProjection,
const std::vector<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& gridCellValues );
static double calculateSum( const RimContourMapProjection& contourMapProjection,
static double calculateSum( const RigContourMapProjection& contourMapProjection,
const std::vector<std::pair<size_t, double>>& matchingCells,
const std::vector<double>& gridCellValues );
};

467
ApplicationLibCode/ReservoirDataModel/RigContourMapProjection.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "RigContourMapProjection.h"
#include "RiaWeightedMeanCalculator.h"
#include "RigContourMapCalculator.h"
#include "RigContourMapGrid.h"
#include "cvfArray.h"
#include "cvfGeometryTools.h"
#include <algorithm>
#include <array>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigContourMapProjection::RigContourMapProjection( const RigContourMapGrid& contourMapGrid )
: m_contourMapGrid( contourMapGrid )
, m_currentResultTimestep( -1 )
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<std::pair<size_t, double>>>
RigContourMapProjection::generateGridMapping( RigContourMapCalculator::ResultAggregationEnum resultAggregation,
const std::vector<double>& 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<double>::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<double>::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<size_t>( gridSize.x() ) * static_cast<size_t>( 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<double>::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<double>& gridCellValues,
RigContourMapCalculator::ResultAggregationEnum resultAggregation ) const
{
const std::vector<std::pair<size_t, double>>& matchingCells = cellsAtIJ( i, j );
return RigContourMapCalculator::calculateValueInMapCell( *this, matchingCells, gridCellValues, resultAggregation );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigContourMapProjection::maxValue( const std::vector<double>& aggregatedResults )
{
double maxV = -std::numeric_limits<double>::infinity();
for ( size_t index = 0; index < aggregatedResults.size(); ++index )
{
if ( aggregatedResults[index] != std::numeric_limits<double>::infinity() )
{
maxV = std::max( maxV, aggregatedResults[index] );
}
}
return maxV;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigContourMapProjection::minValue( const std::vector<double>& aggregatedResults )
{
double minV = std::numeric_limits<double>::infinity();
for ( size_t index = 0; index < aggregatedResults.size(); ++index )
{
if ( aggregatedResults[index] != std::numeric_limits<double>::infinity() )
{
minV = std::min( minV, aggregatedResults[index] );
}
}
return minV;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigContourMapProjection::setCellVisibility( cvf::ref<cvf::UByteArray> cellVisibility )
{
m_cellGridIdxVisibility = cellVisibility;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::ref<cvf::UByteArray> 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<double>( nVertices, std::numeric_limits<double>::infinity() );
#pragma omp parallel for
for ( int index = 0; index < static_cast<int>( nVertices ); ++index )
{
cvf::Vec2ui ij = m_contourMapGrid.ijFromVertexIndex( index );
m_aggregatedVertexResults[index] = calculateValueAtVertex( ij.x(), ij.y() );
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigContourMapProjection::sumTriangleAreas( const std::vector<cvf::Vec4d>& 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<cvf::Vec3d, 4> 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<cvf::Vec2ui, 4> 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<double, 4> 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<double>::infinity() )
{
return std::numeric_limits<double>::infinity();
}
else
{
vertexValues[i] = vertexValue;
validBarycentricCoordsSum += baryCentricCoords[i];
}
}
if ( validBarycentricCoordsSum < 1.0e-8 )
{
return std::numeric_limits<double>::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<double>::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<unsigned int> averageIs;
std::vector<unsigned int> 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<double> 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<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<size_t, double>> 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<std::pair<size_t, double>>();
}
//--------------------------------------------------------------------------------------------------
/// Vertex positions in local coordinates (add origin2d.x() for UTM x)
//--------------------------------------------------------------------------------------------------
std::vector<double> RigContourMapProjection::xVertexPositions() const
{
return m_contourMapGrid.xVertexPositions();
}
//--------------------------------------------------------------------------------------------------
/// Vertex positions in local coordinates (add origin2d.y() for UTM y)
//--------------------------------------------------------------------------------------------------
std::vector<double> RigContourMapProjection::yVertexPositions() const
{
return m_contourMapGrid.yVertexPositions();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigContourMapProjection::aggregatedResults() const
{
return m_aggregatedResults;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigContourMapProjection::aggregatedVertexResults() const
{
return m_aggregatedVertexResults;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<std::vector<std::pair<size_t, double>>>& RigContourMapProjection::projected3dGridIndices() const
{
return m_projected3dGridIndices;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigContourMapProjection::clearResults()
{
m_aggregatedResults.clear();
m_aggregatedVertexResults.clear();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigContourMapProjection::clearGridMapping()
{
m_projected3dGridIndices.clear();
}

119
ApplicationLibCode/ReservoirDataModel/RigContourMapProjection.h Normal file
View File

@ -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 <http://www.gnu.org/licenses/gpl.html>
// 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<size_t, double>;
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<double> xVertexPositions() const;
std::vector<double> 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<size_t> 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<double>& parameterWeights ) const = 0;
virtual std::vector<bool> getMapCellVisibility( int viewStepIndex, RigContourMapCalculator::ResultAggregationEnum resultAggregation ) = 0;
void setCellVisibility( cvf::ref<cvf::UByteArray> cellVisibility );
cvf::ref<cvf::UByteArray> getCellVisibility() const;
std::vector<std::vector<std::pair<size_t, double>>>
generateGridMapping( RigContourMapCalculator::ResultAggregationEnum resultAggregation, const std::vector<double>& weights );
double interpolateValue( const cvf::Vec2d& gridPosition2d ) const;
const std::vector<double>& aggregatedResults() const;
const std::vector<double>& aggregatedVertexResults() const;
const std::vector<std::vector<std::pair<size_t, double>>>& projected3dGridIndices() const;
// Cell index and position conversion
std::vector<CellIndexAndResult> cellsAtIJ( unsigned int i, unsigned int j ) const;
static double maxValue( const std::vector<double>& aggregatedResults );
static double minValue( const std::vector<double>& aggregatedResults );
static double sumTriangleAreas( const std::vector<cvf::Vec4d>& triangles );
protected:
double calculateValueInMapCell( unsigned int i,
unsigned int j,
const std::vector<double>& 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<cvf::UByteArray> m_cellGridIdxVisibility;
std::vector<double> m_aggregatedResults;
std::vector<double> m_aggregatedVertexResults;
std::vector<std::vector<std::pair<size_t, double>>> m_projected3dGridIndices;
int m_currentResultTimestep;
std::vector<bool> m_mapCellVisibility;
const RigContourMapGrid& m_contourMapGrid;
};

334
ApplicationLibCode/ReservoirDataModel/RigEclipseContourMapProjection.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/gpl.html>
// 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 <algorithm>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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<double> RigEclipseContourMapProjection::generateResults( const RigEclipseResultAddress& resultAddress,
RigContourMapCalculator::ResultAggregationEnum resultAggregation,
int timeStep ) const
{
std::pair<bool, std::vector<double>> result =
generateResults( *this, m_contourMapGrid, m_resultData, resultAddress, resultAggregation, timeStep );
return result.second;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<bool, std::vector<double>>
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<double> aggregatedResults = std::vector<double>( nCells, std::numeric_limits<double>::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<double> 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<int>( nCells ); ++index )
{
cvf::Vec2ui ij = contourMapGrid.ijFromCellIndex( index );
const std::vector<std::pair<size_t, double>>& matchingCells = contourMapProjection.cellsAtIJ( ij.x(), ij.y() );
aggregatedResults[index] =
RigContourMapCalculator::calculateValueInMapCell( contourMapProjection, matchingCells, gridResultValues, resultAggregation );
}
}
}
return { useActiveCellInfo, aggregatedResults };
}
std::vector<double> RigEclipseContourMapProjection::calculateColumnResult( RigContourMapCalculator::ResultAggregationEnum resultAggregation,
int timeStep ) const
{
return calculateColumnResult( m_resultData, resultAggregation, timeStep );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> 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<double>();
}
const std::vector<double>& poroResults =
resultData.cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "PORO" ), 0 );
const std::vector<double>& ntgResults =
resultData.cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "NTG" ), 0 );
const std::vector<double>& dzResults =
resultData.cellScalarResults( RigEclipseResultAddress( RiaDefines::ResultCatType::STATIC_NATIVE, "DZ" ), 0 );
CVF_ASSERT( poroResults.size() == ntgResults.size() && ntgResults.size() == dzResults.size() );
std::vector<double> resultValues( poroResults.size(), 0.0 );
if ( resultAggregation == RigContourMapCalculator::RESULTS_OIL_COLUMN || resultAggregation == RigContourMapCalculator::RESULTS_HC_COLUMN )
{
const std::vector<double>& 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<double>& 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<size_t> 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<cvf::Vec3d, 8> 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<cvf::Vec3d, 8> 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<cvf::Vec3d, 8> hexCorners;
const RigCell& cell = m_mainGrid->cell( globalCellIdx );
size_t localCellIdx = cell.gridLocalCellIndex();
RigGridBase* localGrid = cell.hostGrid();
localGrid->cellCornerVertices( localCellIdx, hexCorners.data() );
std::vector<HexIntersectionInfo> 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<double>& 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<bool> RigEclipseContourMapProjection::getMapCellVisibility( int viewStepIndex,
RigContourMapCalculator::ResultAggregationEnum resultAggregation )
{
return std::vector<bool>( numberOfCells(), true );
}

88
ApplicationLibCode/ReservoirDataModel/RigEclipseContourMapProjection.h Normal file
View File

@ -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 <http://www.gnu.org/licenses/gpl.html>
// 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<double> generateResults( const RigEclipseResultAddress& resultAddress,
RigContourMapCalculator::ResultAggregationEnum resultAggregation,
int timeStep ) const;
static std::pair<bool, std::vector<double>> generateResults( const RigEclipseContourMapProjection& contourMapProjection,
const RigContourMapGrid& contourMapGrid,
RigCaseCellResultsData& resultData,
const RigEclipseResultAddress& resultAddress,
RigContourMapCalculator::ResultAggregationEnum resultAggregation,
int timeStep );
std::vector<bool> getMapCellVisibility( int viewStepIndex, RigContourMapCalculator::ResultAggregationEnum resultAggregation ) override;
protected:
using CellIndexAndResult = RigContourMapProjection::CellIndexAndResult;
std::vector<size_t> 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<double>& parameterWeights ) const override;
size_t gridResultIndex( size_t globalCellIdx ) const override;
// Eclipse implementation specific data generation methods
std::vector<double> calculateColumnResult( RigContourMapCalculator::ResultAggregationEnum resultAggregation, int timeStep ) const;
static std::vector<double> calculateColumnResult( RigCaseCellResultsData& resultData,
RigContourMapCalculator::ResultAggregationEnum resultAggregation,
int timeStep );
protected:
RigEclipseCaseData& m_eclipseCaseData;
RigCaseCellResultsData& m_resultData;
cvf::ref<RigMainGrid> m_mainGrid;
cvf::ref<RigActiveCellInfo> m_activeCellInfo;
size_t m_kLayers;
bool m_useActiveCellInfo;
};

381
ApplicationLibCode/ReservoirDataModel/RigGeoMechContourMapProjection.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/gpl.html>
// 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 <algorithm>
#include <array>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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<float>& 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<double>::infinity();
if ( validPorValue )
{
std::array<cvf::Vec3d, 8> 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<bool> RigGeoMechContourMapProjection::getMapCellVisibility( int viewStepIndex,
RigContourMapCalculator::ResultAggregationEnum resultAggregation )
{
return getMapCellVisibility( m_currentResultAddr, viewStepIndex, resultAggregation );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<bool> RigGeoMechContourMapProjection::getMapCellVisibility( RigFemResultAddress resAddr,
int viewStepIndex,
RigContourMapCalculator::ResultAggregationEnum resultAggregation )
{
cvf::Vec2ui nCellsIJ = numberOfElementsIJ();
std::vector<std::vector<unsigned int>> distanceImage( nCellsIJ.x(), std::vector<unsigned int>( nCellsIJ.y(), 0u ) );
std::vector<bool> mapCellVisibility;
if ( m_limitToPorePressureRegions )
{
resAddr = RigFemAddressDefines::elementNodalPorBarAddress();
}
std::vector<double> 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<double>::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<double> RigGeoMechContourMapProjection::generateResultsFromAddress( RigFemResultAddress resultAddress,
const std::vector<bool>& mapCellVisibility,
RigContourMapCalculator::ResultAggregationEnum resultAggregation,
int viewerStepIndex ) const
{
RigFemPartResultsCollection* resultCollection = m_caseData.femPartResults();
size_t nCells = numberOfCells();
std::vector<double> aggregatedResults = std::vector<double>( nCells, std::numeric_limits<double>::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<float> resultValuesF = resultCollection->resultValues( resultAddress, 0, stepIdx, frameIdx );
if ( resultValuesF.empty() ) return aggregatedResults;
std::vector<double> 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<double>::infinity() )
{
value = 0.0;
}
}
}
#pragma omp parallel for
for ( int index = 0; index < static_cast<int>( 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<size_t> 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<cvf::Vec3d, 8> hexCorners;
m_femPartGrid->cellCornerVertices( globalCellIdx, hexCorners.data() );
cvf::BoundingBox overlapBBox;
std::array<cvf::Vec3d, 8> 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<cvf::Vec3d, 8> hexCorners;
const std::vector<cvf::Vec3f>& 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<HexIntersectionInfo> 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<double>& parameterWeights ) const
{
if ( parameterWeights.empty() ) return 1.0;
return parameterWeights[globalCellIdx];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RigGeoMechContourMapProjection::gridCellValues( RigFemResultAddress resAddr, std::vector<float>& resultValues ) const
{
std::vector<double> gridCellValues( m_femPart->elementCount(), std::numeric_limits<double>::infinity() );
for ( size_t globalCellIdx = 0; globalCellIdx < static_cast<size_t>( m_femPart->elementCount() ); ++globalCellIdx )
{
RigElementType elmType = m_femPart->elementType( globalCellIdx );
if ( !RigFemTypes::is8NodeElement( elmType ) ) continue;
if ( resAddr.resultPosType == RIG_ELEMENT )
{
gridCellValues[globalCellIdx] = static_cast<double>( resultValues[globalCellIdx] );
}
else if ( resAddr.resultPosType == RIG_ELEMENT_NODAL )
{
RiaWeightedMeanCalculator<float> cellAverage;
for ( int i = 0; i < 8; ++i )
{
size_t gridResultValueIdx =
m_femPart->resultValueIdxFromResultPosType( resAddr.resultPosType, static_cast<int>( globalCellIdx ), i );
cellAverage.addValueAndWeight( resultValues[gridResultValueIdx], 1.0 );
}
gridCellValues[globalCellIdx] = static_cast<double>( cellAverage.weightedMean() );
}
else
{
RiaWeightedMeanCalculator<float> 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<double>( cellAverage.weightedMean() );
}
}
return gridCellValues;
}

81
ApplicationLibCode/ReservoirDataModel/RigGeoMechContourMapProjection.h Normal file
View File

@ -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 <http://www.gnu.org/licenses/gpl.html>
// 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<double> generateResultsFromAddress( RigFemResultAddress resultAddress,
const std::vector<bool>& 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<bool> getMapCellVisibility( int viewStepIndex, RigContourMapCalculator::ResultAggregationEnum resultAggregation ) override;
std::vector<bool> getMapCellVisibility( RigFemResultAddress resAddr,
int viewStepIndex,
RigContourMapCalculator::ResultAggregationEnum resultAggregation );
protected:
// GeoMech implementation specific data generation methods
std::vector<size_t> 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<double>& parameterWeights ) const override;
std::vector<double> gridCellValues( RigFemResultAddress resAddr, std::vector<float>& resultValues ) const;
protected:
RigGeoMechCaseData& m_caseData;
bool m_limitToPorePressureRegions;
double m_paddingAroundPorePressureRegion;
cvf::ref<RigFemPart> m_femPart;
cvf::cref<RigFemPartGrid> m_femPartGrid;
RigFemResultAddress m_currentResultAddr;
size_t m_kLayers;
};