mirror of
https://github.com/OPM/ResInsight.git
synced 2024-12-28 18:01:08 -06:00
415 lines
19 KiB
C++
415 lines
19 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// 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 "RigContourMapCalculator.h"
|
|
|
|
#include "RiaWeightedGeometricMeanCalculator.h"
|
|
#include "RiaWeightedHarmonicMeanCalculator.h"
|
|
#include "RiaWeightedMeanCalculator.h"
|
|
|
|
#include "RigContourMapGrid.h"
|
|
|
|
#include "RigContourMapProjection.h"
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <map>
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapCalculator::calculateValueInMapCell( const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<std::pair<size_t, double>>& matchingCells,
|
|
const std::vector<double>& gridCellValues,
|
|
ResultAggregationType resultAggregation )
|
|
{
|
|
if ( matchingCells.empty() ) return std::numeric_limits<double>::infinity();
|
|
|
|
switch ( resultAggregation )
|
|
{
|
|
case TOP_VALUE:
|
|
return calculateTopValue( contourMapProjection, matchingCells, gridCellValues );
|
|
case MEAN:
|
|
return calculateMeanValue( contourMapProjection, matchingCells, gridCellValues );
|
|
case GEOMETRIC_MEAN:
|
|
return calculateGeometricMeanValue( contourMapProjection, matchingCells, gridCellValues );
|
|
case HARMONIC_MEAN:
|
|
return calculateHarmonicMeanValue( contourMapProjection, matchingCells, gridCellValues );
|
|
case MAX_VALUE:
|
|
return calculateMaxValue( contourMapProjection, matchingCells, gridCellValues );
|
|
case MIN_VALUE:
|
|
return calculateMinValue( contourMapProjection, matchingCells, gridCellValues );
|
|
case VOLUME_SUM:
|
|
case SUM:
|
|
case OIL_COLUMN:
|
|
case GAS_COLUMN:
|
|
case HYDROCARBON_COLUMN:
|
|
return calculateSum( contourMapProjection, matchingCells, gridCellValues );
|
|
default:
|
|
{
|
|
CVF_TIGHT_ASSERT( false );
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapCalculator::calculateTopValue( const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<std::pair<size_t, double>>& matchingCells,
|
|
const std::vector<double>& gridCellValues )
|
|
{
|
|
for ( auto [cellIdx, weight] : matchingCells )
|
|
{
|
|
double cellValue = gridCellValues[contourMapProjection.gridResultIndex( cellIdx )];
|
|
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
|
|
{
|
|
return cellValue;
|
|
}
|
|
}
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapCalculator::calculateMeanValue( const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<std::pair<size_t, double>>& matchingCells,
|
|
const std::vector<double>& gridCellValues )
|
|
{
|
|
RiaWeightedMeanCalculator<double> calculator;
|
|
for ( auto [cellIdx, weight] : matchingCells )
|
|
{
|
|
double cellValue = gridCellValues[contourMapProjection.gridResultIndex( cellIdx )];
|
|
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
|
|
{
|
|
calculator.addValueAndWeight( cellValue, weight );
|
|
}
|
|
}
|
|
if ( calculator.validAggregatedWeight() )
|
|
{
|
|
return calculator.weightedMean();
|
|
}
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapCalculator::calculateGeometricMeanValue( const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<std::pair<size_t, double>>& matchingCells,
|
|
const std::vector<double>& gridCellValues )
|
|
{
|
|
RiaWeightedGeometricMeanCalculator calculator;
|
|
for ( auto [cellIdx, weight] : matchingCells )
|
|
{
|
|
double cellValue = gridCellValues[contourMapProjection.gridResultIndex( cellIdx )];
|
|
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
|
|
{
|
|
if ( cellValue < 1.0e-8 )
|
|
{
|
|
return 0.0;
|
|
}
|
|
calculator.addValueAndWeight( cellValue, weight );
|
|
}
|
|
}
|
|
if ( calculator.validAggregatedWeight() )
|
|
{
|
|
return calculator.weightedMean();
|
|
}
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapCalculator::calculateHarmonicMeanValue( const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<std::pair<size_t, double>>& matchingCells,
|
|
const std::vector<double>& gridCellValues )
|
|
{
|
|
RiaWeightedHarmonicMeanCalculator calculator;
|
|
for ( auto [cellIdx, weight] : matchingCells )
|
|
{
|
|
double cellValue = gridCellValues[contourMapProjection.gridResultIndex( cellIdx )];
|
|
if ( std::fabs( cellValue ) < 1.0e-8 )
|
|
{
|
|
return 0.0;
|
|
}
|
|
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
|
|
{
|
|
calculator.addValueAndWeight( cellValue, weight );
|
|
}
|
|
}
|
|
if ( calculator.validAggregatedWeight() )
|
|
{
|
|
return calculator.weightedMean();
|
|
}
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapCalculator::calculateMaxValue( const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<std::pair<size_t, double>>& matchingCells,
|
|
const std::vector<double>& gridCellValues )
|
|
{
|
|
double maxValue = -std::numeric_limits<double>::infinity();
|
|
for ( auto [cellIdx, weight] : matchingCells )
|
|
{
|
|
double cellValue = gridCellValues[contourMapProjection.gridResultIndex( cellIdx )];
|
|
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
|
|
{
|
|
maxValue = std::max( maxValue, cellValue );
|
|
}
|
|
}
|
|
if ( maxValue == -std::numeric_limits<double>::infinity() )
|
|
{
|
|
maxValue = std::numeric_limits<double>::infinity();
|
|
}
|
|
return maxValue;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapCalculator::calculateMinValue( const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<std::pair<size_t, double>>& matchingCells,
|
|
const std::vector<double>& gridCellValues )
|
|
{
|
|
double minValue = std::numeric_limits<double>::infinity();
|
|
for ( auto [cellIdx, weight] : matchingCells )
|
|
{
|
|
double cellValue = gridCellValues[contourMapProjection.gridResultIndex( cellIdx )];
|
|
minValue = std::min( minValue, cellValue );
|
|
}
|
|
return minValue;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapCalculator::calculateSum( const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<std::pair<size_t, double>>& matchingCells,
|
|
const std::vector<double>& gridCellValues )
|
|
{
|
|
double sum = 0.0;
|
|
for ( auto [cellIdx, weight] : matchingCells )
|
|
{
|
|
double cellValue = gridCellValues[contourMapProjection.gridResultIndex( cellIdx )];
|
|
if ( std::abs( cellValue ) != std::numeric_limits<double>::infinity() )
|
|
{
|
|
sum += cellValue * weight;
|
|
}
|
|
}
|
|
return sum;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<std::vector<std::pair<size_t, double>>>
|
|
RigContourMapCalculator::generateGridMapping( RigContourMapProjection& contourMapProjection,
|
|
const RigContourMapGrid& contourMapGrid,
|
|
ResultAggregationType resultAggregation,
|
|
const std::vector<double>& weightingResultValues )
|
|
{
|
|
int nCells = contourMapGrid.numberOfCells();
|
|
std::vector<std::vector<std::pair<size_t, double>>> projected3dGridIndices( nCells );
|
|
|
|
if ( RigContourMapCalculator::isStraightSummationResult( resultAggregation ) )
|
|
{
|
|
#pragma omp parallel for
|
|
for ( int index = 0; index < nCells; ++index )
|
|
{
|
|
cvf::Vec2ui ij = contourMapGrid.ijFromCellIndex( index );
|
|
|
|
cvf::Vec2d globalPos = contourMapGrid.cellCenterPosition( ij.x(), ij.y() ) + contourMapGrid.origin2d();
|
|
projected3dGridIndices[index] =
|
|
cellRayIntersectionAndResults( contourMapProjection, contourMapGrid, globalPos, weightingResultValues );
|
|
}
|
|
}
|
|
else
|
|
{
|
|
#pragma omp parallel for
|
|
for ( int index = 0; index < nCells; ++index )
|
|
{
|
|
cvf::Vec2ui ij = contourMapGrid.ijFromCellIndex( index );
|
|
|
|
cvf::Vec2d globalPos = contourMapGrid.cellCenterPosition( ij.x(), ij.y() ) + contourMapGrid.origin2d();
|
|
projected3dGridIndices[index] =
|
|
cellOverlapVolumesAndResults( contourMapProjection, contourMapGrid, globalPos, weightingResultValues );
|
|
}
|
|
}
|
|
|
|
return projected3dGridIndices;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<RigContourMapCalculator::CellIndexAndResult>
|
|
RigContourMapCalculator::cellOverlapVolumesAndResults( const RigContourMapProjection& contourMapProjection,
|
|
const RigContourMapGrid& contourMapGrid,
|
|
const cvf::Vec2d& globalPos2d,
|
|
const std::vector<double>& weightingResultValues )
|
|
{
|
|
const cvf::BoundingBox& expandedBoundingBox = contourMapGrid.expandedBoundingBox();
|
|
cvf::Vec3d top2dElementCentroid( globalPos2d, expandedBoundingBox.max().z() );
|
|
cvf::Vec3d bottom2dElementCentroid( globalPos2d, expandedBoundingBox.min().z() );
|
|
cvf::Vec3d planarDiagonalVector( 0.5 * contourMapGrid.sampleSpacing(), 0.5 * contourMapGrid.sampleSpacing(), 0.0 );
|
|
cvf::Vec3d topNECorner = top2dElementCentroid + planarDiagonalVector;
|
|
cvf::Vec3d bottomSWCorner = bottom2dElementCentroid - planarDiagonalVector;
|
|
|
|
cvf::BoundingBox bbox2dElement( bottomSWCorner, topNECorner );
|
|
|
|
std::vector<std::pair<size_t, double>> matchingVisibleCellsAndWeight;
|
|
|
|
// Bounding box has been expanded, so 2d element may be outside actual 3d grid
|
|
if ( !bbox2dElement.intersects( contourMapGrid.originalBoundingBox() ) )
|
|
{
|
|
return matchingVisibleCellsAndWeight;
|
|
}
|
|
|
|
std::vector<size_t> allCellIndices = contourMapProjection.findIntersectingCells( bbox2dElement );
|
|
|
|
std::vector<std::vector<size_t>> kLayerCellIndexVector;
|
|
kLayerCellIndexVector.resize( contourMapProjection.kLayers() );
|
|
|
|
if ( kLayerCellIndexVector.empty() )
|
|
{
|
|
return matchingVisibleCellsAndWeight;
|
|
}
|
|
|
|
auto cellGridIdxVisibility = contourMapProjection.getCellVisibility();
|
|
for ( size_t globalCellIdx : allCellIndices )
|
|
{
|
|
if ( cellGridIdxVisibility.isNull() || ( *cellGridIdxVisibility )[globalCellIdx] )
|
|
{
|
|
kLayerCellIndexVector[contourMapProjection.kLayer( globalCellIdx )].push_back( globalCellIdx );
|
|
}
|
|
}
|
|
|
|
for ( const auto& kLayerIndices : kLayerCellIndexVector )
|
|
{
|
|
for ( size_t globalCellIdx : kLayerIndices )
|
|
{
|
|
double overlapVolume = contourMapProjection.calculateOverlapVolume( globalCellIdx, bbox2dElement );
|
|
if ( overlapVolume > 0.0 )
|
|
{
|
|
double weight =
|
|
overlapVolume * contourMapProjection.getParameterWeightForCell( contourMapProjection.gridResultIndex( globalCellIdx ),
|
|
weightingResultValues );
|
|
if ( weight > 0.0 )
|
|
{
|
|
matchingVisibleCellsAndWeight.push_back( std::make_pair( globalCellIdx, weight ) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return matchingVisibleCellsAndWeight;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<RigContourMapCalculator::CellIndexAndResult>
|
|
RigContourMapCalculator::cellRayIntersectionAndResults( const RigContourMapProjection& contourMapProjection,
|
|
const RigContourMapGrid& contourMapGrid,
|
|
const cvf::Vec2d& globalPos2d,
|
|
const std::vector<double>& weightingResultValues )
|
|
{
|
|
std::vector<std::pair<size_t, double>> matchingVisibleCellsAndWeight;
|
|
|
|
const cvf::BoundingBox& expandedBoundingBox = contourMapGrid.expandedBoundingBox();
|
|
|
|
cvf::Vec3d highestPoint( globalPos2d, expandedBoundingBox.max().z() );
|
|
cvf::Vec3d lowestPoint( globalPos2d, expandedBoundingBox.min().z() );
|
|
|
|
// Bounding box has been expanded, so ray may be outside actual 3d grid
|
|
if ( !contourMapGrid.originalBoundingBox().contains( highestPoint ) )
|
|
{
|
|
return matchingVisibleCellsAndWeight;
|
|
}
|
|
|
|
cvf::BoundingBox rayBBox;
|
|
rayBBox.add( highestPoint );
|
|
rayBBox.add( lowestPoint );
|
|
|
|
std::vector<size_t> allCellIndices = contourMapProjection.findIntersectingCells( rayBBox );
|
|
|
|
std::map<size_t, std::vector<size_t>> kLayerIndexMap;
|
|
|
|
auto cellGridIdxVisibility = contourMapProjection.getCellVisibility();
|
|
for ( size_t globalCellIdx : allCellIndices )
|
|
{
|
|
if ( cellGridIdxVisibility.isNull() || ( *cellGridIdxVisibility )[globalCellIdx] )
|
|
{
|
|
kLayerIndexMap[contourMapProjection.kLayer( globalCellIdx )].push_back( globalCellIdx );
|
|
}
|
|
}
|
|
|
|
for ( const auto& kLayerIndexPair : kLayerIndexMap )
|
|
{
|
|
double weightSumThisKLayer = 0.0;
|
|
std::vector<std::pair<size_t, double>> cellsAndWeightsThisLayer;
|
|
for ( size_t globalCellIdx : kLayerIndexPair.second )
|
|
{
|
|
double lengthInCell = contourMapProjection.calculateRayLengthInCell( globalCellIdx, highestPoint, lowestPoint );
|
|
if ( lengthInCell > 0.0 )
|
|
{
|
|
cellsAndWeightsThisLayer.push_back( std::make_pair( globalCellIdx, lengthInCell ) );
|
|
weightSumThisKLayer += lengthInCell;
|
|
}
|
|
}
|
|
for ( auto& cellWeightPair : cellsAndWeightsThisLayer )
|
|
{
|
|
cellWeightPair.second /= weightSumThisKLayer;
|
|
matchingVisibleCellsAndWeight.push_back( cellWeightPair );
|
|
}
|
|
}
|
|
|
|
return matchingVisibleCellsAndWeight;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RigContourMapCalculator::isColumnResult( ResultAggregationType aggregationType )
|
|
{
|
|
return aggregationType == OIL_COLUMN || aggregationType == GAS_COLUMN || aggregationType == HYDROCARBON_COLUMN;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RigContourMapCalculator::isMeanResult( ResultAggregationType aggregationType )
|
|
{
|
|
return aggregationType == MEAN || aggregationType == HARMONIC_MEAN || aggregationType == GEOMETRIC_MEAN;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RigContourMapCalculator::isStraightSummationResult( ResultAggregationType aggregationType )
|
|
{
|
|
return aggregationType == OIL_COLUMN || aggregationType == GAS_COLUMN || aggregationType == HYDROCARBON_COLUMN || aggregationType == SUM;
|
|
}
|