OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-11 07:56:08 -06:00
Code Issues Packages Projects Releases Wiki Activity

#3597 2d Maps: Implement Top, Sum, Harmonic and Geometric mean aggregation

Browse Source
This commit is contained in:
Gaute Lindkvist 2018-11-01 15:21:35 +01:00
parent 0291dd2237
commit 9730048319
2 changed files with 134 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

164
ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp
View File

@ -1,5 +1,7 @@
#include "Rim2dGridProjection.h"
#include "RiaWeightedGeometricMeanCalculator.h"
#include "RiaWeightedHarmonicMeanCalculator.h"
#include "RiaWeightedMeanCalculator.h"
#include "RigActiveCellInfo.h"
@ -37,11 +39,15 @@ namespace caf
template<>
void Rim2dGridProjection::ResultAggregation::setUp()
{
addItem(Rim2dGridProjection::RESULTS_MEAN_VALUE, "MEAN_VALUE", "Mean Value");
addItem(Rim2dGridProjection::RESULTS_TOP_VALUE, "TOP_VALUE", "Top Value");
addItem(Rim2dGridProjection::RESULTS_MEAN_VALUE, "MEAN_VALUE", "Arithmetic Mean");
addItem(Rim2dGridProjection::RESULTS_HARM_VALUE, "HARM_VALUE", "Harmonic Mean");
addItem(Rim2dGridProjection::RESULTS_GEOM_VALUE, "GEOM_VALUE", "Geometric Mean");
addItem(Rim2dGridProjection::RESULTS_MIN_VALUE, "MIN_VALUE", "Min Value");
addItem(Rim2dGridProjection::RESULTS_MAX_VALUE, "MAX_VALUE", "Max Value");
addItem(Rim2dGridProjection::RESULTS_SUM, "SUM", "Sum");
setDefault(Rim2dGridProjection::RESULTS_MEAN_VALUE);
setDefault(Rim2dGridProjection::RESULTS_TOP_VALUE);
}
}
CAF_PDM_SOURCE_INIT(Rim2dGridProjection, "Rim2dGridProjection");
@ -283,6 +289,12 @@ double Rim2dGridProjection::value(uint i, uint j) const
{
switch (m_resultAggregation())
{
case RESULTS_TOP_VALUE:
{
size_t cellIdx = matchingCells.front().first;
double cellValue = m_resultAccessor->cellScalarGlobIdx(cellIdx);
return cellValue;
}
case RESULTS_MEAN_VALUE:
{
RiaWeightedMeanCalculator<float> calculator;
@ -294,6 +306,28 @@ double Rim2dGridProjection::value(uint i, uint j) const
}
return calculator.weightedMean();
}
case RESULTS_GEOM_VALUE:
{
RiaWeightedGeometricMeanCalculator calculator;
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = m_resultAccessor->cellScalarGlobIdx(cellIdx);
calculator.addValueAndWeight(cellValue, cellIdxAndWeight.second);
}
return calculator.weightedMean();
}
case RESULTS_HARM_VALUE:
{
RiaWeightedHarmonicMeanCalculator calculator;
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = m_resultAccessor->cellScalarGlobIdx(cellIdx);
calculator.addValueAndWeight(cellValue, cellIdxAndWeight.second);
}
return calculator.weightedMean();
}
case RESULTS_MAX_VALUE:
{
double maxValue = -std::numeric_limits<double>::infinity();
@ -316,6 +350,17 @@ double Rim2dGridProjection::value(uint i, uint j) const
}
return minValue;
}
case RESULTS_SUM:
{
double sum = 0.0;
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = m_resultAccessor->cellScalarGlobIdx(cellIdx);
sum += cellValue * cellIdxAndWeight.second;
}
return sum;
}
default:
CVF_TIGHT_ASSERT(false);
}
@ -377,45 +422,51 @@ RimRegularLegendConfig* Rim2dGridProjection::legendConfig() const
//--------------------------------------------------------------------------------------------------
void Rim2dGridProjection::calculateCellRangeVisibility()
{
const RigMainGrid* grid = mainGrid();
m_cellVisibility = new cvf::UByteArray(grid->cellCount());
for (size_t gridIndex = 0u; gridIndex < mainGrid()->gridCount(); ++gridIndex)
{
const RigGridBase* grid = mainGrid()->gridByIndex(gridIndex);
m_cellGridIdxVisibilityMap[gridIndex] = new cvf::UByteArray(grid->cellCount());
RimEclipseView* view = nullptr;
firstAncestorOrThisOfTypeAsserted(view);
RimEclipseView* view = nullptr;
firstAncestorOrThisOfTypeAsserted(view);
RimCellRangeFilterCollection* rangeFilterCollection = view->rangeFilterCollection();
RimCellRangeFilterCollection* rangeFilterCollection = view->rangeFilterCollection();
const RigActiveCellInfo* activeCellInfo = eclipseCase()->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL);
const RigActiveCellInfo* activeCellInfo = eclipseCase()->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL);
#pragma omp parallel for
for (int cellIndex = 0; cellIndex < static_cast<int>(grid->cellCount()); ++cellIndex)
{
(*m_cellVisibility)[cellIndex] = activeCellInfo->isActive(cellIndex);
}
if (rangeFilterCollection && rangeFilterCollection->isActive())
{
cvf::CellRangeFilter cellRangeFilter;
rangeFilterCollection->compoundCellRangeFilter(&cellRangeFilter, grid->gridIndex());
if (cellRangeFilter.hasIncludeRanges())
for (int cellIndex = 0; cellIndex < static_cast<int>(grid->cellCount()); ++cellIndex)
{
#pragma omp parallel for
for (int cellIndex = 0; cellIndex < static_cast<int>(grid->cellCount()); ++cellIndex)
{
size_t i, j, k;
grid->ijkFromCellIndex(cellIndex, &i, &j, &k);
(*m_cellVisibility)[cellIndex] = (*m_cellVisibility)[cellIndex] && cellRangeFilter.isCellVisible(i, j, k, false);
}
size_t globalCellIdx = grid->reservoirCellIndex(cellIndex);
(*m_cellGridIdxVisibilityMap[gridIndex])[cellIndex] = activeCellInfo->isActive(globalCellIdx);
}
else
if (rangeFilterCollection && rangeFilterCollection->isActive())
{
#pragma omp parallel for
for (int cellIndex = 0; cellIndex < static_cast<int>(grid->cellCount()); ++cellIndex)
cvf::CellRangeFilter cellRangeFilter;
rangeFilterCollection->compoundCellRangeFilter(&cellRangeFilter, gridIndex);
if (cellRangeFilter.hasIncludeRanges())
{
size_t i, j, k;
grid->ijkFromCellIndex(cellIndex, &i, &j, &k);
(*m_cellVisibility)[cellIndex] = (*m_cellVisibility)[cellIndex] && !cellRangeFilter.isCellExcluded(i, j, k, false);
#pragma omp parallel for
for (int cellIndex = 0; cellIndex < static_cast<int>(grid->cellCount()); ++cellIndex)
{
size_t i, j, k;
grid->ijkFromCellIndex(cellIndex, &i, &j, &k);
(*m_cellGridIdxVisibilityMap[gridIndex])[cellIndex] =
(*m_cellGridIdxVisibilityMap[gridIndex])[cellIndex] && cellRangeFilter.isCellVisible(i, j, k, false);
}
}
else
{
#pragma omp parallel for
for (int cellIndex = 0; cellIndex < static_cast<int>(grid->cellCount()); ++cellIndex)
{
size_t i, j, k;
grid->ijkFromCellIndex(cellIndex, &i, &j, &k);
(*m_cellGridIdxVisibilityMap[gridIndex])[cellIndex] =
(*m_cellGridIdxVisibilityMap[gridIndex])[cellIndex] && !cellRangeFilter.isCellExcluded(i, j, k, false);
}
}
}
}
@ -430,7 +481,11 @@ void Rim2dGridProjection::calculatePropertyFilterVisibility()
firstAncestorOrThisOfTypeAsserted(view);
int timeStep = view->currentTimeStep();
RivReservoirViewPartMgr::computePropertyVisibility(m_cellVisibility.p(), mainGrid(), timeStep, m_cellVisibility.p(), view->eclipsePropertyFilterCollection());
for (size_t gridIndex = 0u; gridIndex < mainGrid()->gridCount(); ++gridIndex)
{
const RigGridBase* grid = mainGrid()->gridByIndex(gridIndex);
RivReservoirViewPartMgr::computePropertyVisibility(m_cellGridIdxVisibilityMap[gridIndex].p(), grid, timeStep, m_cellGridIdxVisibilityMap[gridIndex].p(), view->eclipsePropertyFilterCollection());
}
}
//--------------------------------------------------------------------------------------------------
@ -471,29 +526,50 @@ std::vector<std::pair<size_t, float>> Rim2dGridProjection::visibleCellsAndWeight
std::vector<size_t> allCellIndices;
mainGrid()->findIntersectingCells(rayBBox, &allCellIndices);
std::vector<std::pair<size_t, float>> matchingVisibleCellsAndWeight;
std::vector<std::tuple<size_t, float, float>> matchingVisibleCellsWeightAndHeight;
cvf::Vec3d hexCorners[8];
for (size_t globalCellIdx : allCellIndices)
{
size_t localCellIdx = 0u;
RigGridBase* localGrid = mainGrid()->gridAndGridLocalIdxFromGlobalCellIdx(globalCellIdx, &localCellIdx);
if (localGrid == mainGrid())
if ((*m_cellGridIdxVisibilityMap.at(localGrid->gridIndex()))[localCellIdx])
{
if ((*m_cellVisibility)[globalCellIdx])
localGrid->cellCornerVertices(localCellIdx, hexCorners);
std::vector<HexIntersectionInfo> intersections;
double height = 0.0;
double weight = 1.0;
if (RigHexIntersectionTools::lineHexCellIntersection(highestPoint, lowestPoint, hexCorners, 0, &intersections))
{
localGrid->cellCornerVertices(localCellIdx, hexCorners);
std::vector<HexIntersectionInfo> intersections;
float weight = 1.0f;
if (RigHexIntersectionTools::lineHexCellIntersection(highestPoint, lowestPoint, hexCorners, 0, &intersections))
{
weight = std::max(1.0, (intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint).length());
}
matchingVisibleCellsAndWeight.push_back(std::make_pair(globalCellIdx, weight));
height = intersections.front().m_intersectionPoint.z();
weight = (intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint).length();
}
else
{
for (cvf::Vec3d corner : hexCorners)
{
height += corner.z();
}
height /= 8;
}
matchingVisibleCellsWeightAndHeight.push_back(std::make_tuple(globalCellIdx, weight, height));
}
}
std::sort(matchingVisibleCellsWeightAndHeight.begin(), matchingVisibleCellsWeightAndHeight.end(), [](const auto& lhs, const auto& rhs)
{
return std::get<2>(lhs) > std::get<2>(rhs);
});
std::vector<std::pair<size_t, float>> matchingVisibleCellsAndWeight;
for (const auto& visWeightHeight : matchingVisibleCellsWeightAndHeight)
{
matchingVisibleCellsAndWeight.push_back(std::make_pair(std::get<0>(visWeightHeight), std::get<1>(visWeightHeight)));
}
return matchingVisibleCellsAndWeight;
}

10
ApplicationCode/ProjectDataModel/Rim2dGridProjection.h
View File

@ -45,9 +45,13 @@ class Rim2dGridProjection : public RimCheckableNamedObject
public:
enum ResultAggregationEnum
{
RESULTS_MEAN_VALUE = 0,
RESULTS_TOP_VALUE,
RESULTS_MEAN_VALUE,
RESULTS_GEOM_VALUE,
RESULTS_HARM_VALUE,
RESULTS_MIN_VALUE,
RESULTS_MAX_VALUE
RESULTS_MAX_VALUE,
RESULTS_SUM
};
typedef caf::AppEnum<ResultAggregationEnum> ResultAggregation;
typedef std::vector<cvf::ref<cvf::Vec3fArray>> ContourPolygons;
@ -98,7 +102,7 @@ protected:
protected:
caf::PdmField<double> m_sampleSpacing;
caf::PdmField<ResultAggregation> m_resultAggregation;
cvf::ref<cvf::UByteArray> m_cellVisibility;
std::map<size_t, cvf::ref<cvf::UByteArray>> m_cellGridIdxVisibilityMap;
std::vector<double> m_aggregatedResults;
std::vector<std::vector<std::pair<size_t, float>>> m_projected3dGridIndices;