From 9730048319b24685fa476ab93b11ce9c9c1e3553 Mon Sep 17 00:00:00 2001 From: Gaute Lindkvist Date: Thu, 1 Nov 2018 15:21:35 +0100 Subject: [PATCH] #3597 2d Maps: Implement Top, Sum, Harmonic and Geometric mean aggregation --- .../ProjectDataModel/Rim2dGridProjection.cpp | 178 +++++++++++++----- .../ProjectDataModel/Rim2dGridProjection.h | 10 +- 2 files changed, 134 insertions(+), 54 deletions(-) diff --git a/ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp b/ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp index 4ebf9e23aa..997632c056 100644 --- a/ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp +++ b/ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp @@ -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 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::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,48 +422,54 @@ 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); - - RimCellRangeFilterCollection* rangeFilterCollection = view->rangeFilterCollection(); + RimEclipseView* view = nullptr; + firstAncestorOrThisOfTypeAsserted(view); - const RigActiveCellInfo* activeCellInfo = eclipseCase()->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL); + RimCellRangeFilterCollection* rangeFilterCollection = view->rangeFilterCollection(); + + const RigActiveCellInfo* activeCellInfo = eclipseCase()->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL); #pragma omp parallel for - for (int cellIndex = 0; cellIndex < static_cast(grid->cellCount()); ++cellIndex) - { - (*m_cellVisibility)[cellIndex] = activeCellInfo->isActive(cellIndex); + for (int cellIndex = 0; cellIndex < static_cast(grid->cellCount()); ++cellIndex) + { + size_t globalCellIdx = grid->reservoirCellIndex(cellIndex); + (*m_cellGridIdxVisibilityMap[gridIndex])[cellIndex] = activeCellInfo->isActive(globalCellIdx); + } + + if (rangeFilterCollection && rangeFilterCollection->isActive()) + { + cvf::CellRangeFilter cellRangeFilter; + rangeFilterCollection->compoundCellRangeFilter(&cellRangeFilter, gridIndex); + + if (cellRangeFilter.hasIncludeRanges()) + { +#pragma omp parallel for + for (int cellIndex = 0; cellIndex < static_cast(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(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); + } + } + } } - - if (rangeFilterCollection && rangeFilterCollection->isActive()) - { - cvf::CellRangeFilter cellRangeFilter; - rangeFilterCollection->compoundCellRangeFilter(&cellRangeFilter, grid->gridIndex()); - - if (cellRangeFilter.hasIncludeRanges()) - { -#pragma omp parallel for - for (int cellIndex = 0; cellIndex < static_cast(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); - } - } - else - { -#pragma omp parallel for - for (int cellIndex = 0; cellIndex < static_cast(grid->cellCount()); ++cellIndex) - { - size_t i, j, k; - grid->ijkFromCellIndex(cellIndex, &i, &j, &k); - (*m_cellVisibility)[cellIndex] = (*m_cellVisibility)[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> Rim2dGridProjection::visibleCellsAndWeight std::vector allCellIndices; mainGrid()->findIntersectingCells(rayBBox, &allCellIndices); - std::vector> matchingVisibleCellsAndWeight; + + std::vector> 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 intersections; + + double height = 0.0; + double weight = 1.0; + if (RigHexIntersectionTools::lineHexCellIntersection(highestPoint, lowestPoint, hexCorners, 0, &intersections)) { - localGrid->cellCornerVertices(localCellIdx, hexCorners); - std::vector 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> matchingVisibleCellsAndWeight; + for (const auto& visWeightHeight : matchingVisibleCellsWeightAndHeight) + { + matchingVisibleCellsAndWeight.push_back(std::make_pair(std::get<0>(visWeightHeight), std::get<1>(visWeightHeight))); + } + + return matchingVisibleCellsAndWeight; } diff --git a/ApplicationCode/ProjectDataModel/Rim2dGridProjection.h b/ApplicationCode/ProjectDataModel/Rim2dGridProjection.h index 9e97ca393b..98e1b27e5d 100644 --- a/ApplicationCode/ProjectDataModel/Rim2dGridProjection.h +++ b/ApplicationCode/ProjectDataModel/Rim2dGridProjection.h @@ -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 ResultAggregation; typedef std::vector> ContourPolygons; @@ -98,7 +102,7 @@ protected: protected: caf::PdmField m_sampleSpacing; caf::PdmField m_resultAggregation; - cvf::ref m_cellVisibility; + std::map> m_cellGridIdxVisibilityMap; std::vector m_aggregatedResults; std::vector>> m_projected3dGridIndices;