OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-01-26 00:06:49 -06:00
Code Issues Packages Projects Releases Wiki Activity
c99fb71c36
ResInsight/ApplicationCode/ProjectDataModel/RimContourMapProjection.cpp

1167 lines
45 KiB
C++
Raw Blame History

#include "RimContourMapProjection.h"
#include "RiaWeightedGeometricMeanCalculator.h"
#include "RiaWeightedHarmonicMeanCalculator.h"
#include "RiaWeightedMeanCalculator.h"
#include "RigActiveCellInfo.h"
#include "RigCaseCellResultsData.h"
#include "RigCell.h"
#include "RigCellGeometryTools.h"
#include "RigEclipseCaseData.h"
#include "RigHexIntersectionTools.h"
#include "RigMainGrid.h"
#include "RigResultAccessor.h"
#include "RigResultAccessorFactory.h"
#include "RivReservoirViewPartMgr.h"
#include "RimCellRangeFilterCollection.h"
#include "RimContourMapView.h"
#include "RimEclipseCellColors.h"
#include "RimEclipseView.h"
#include "RimEclipseResultCase.h"
#include "RimEclipseResultDefinition.h"
#include "RimProject.h"
#include "RimRegularLegendConfig.h"
#include "cafContourLines.h"
#include "cafPdmUiDoubleSliderEditor.h"
#include "cafPdmUiTreeOrdering.h"
#include "cvfArray.h"
#include "cvfCellRange.h"
#include "cvfScalarMapper.h"
#include "cvfStructGridGeometryGenerator.h"
#include <QDebug>
#include <algorithm>
namespace caf
{
template<>
void RimContourMapProjection::ResultAggregation::setUp()
{
addItem(RimContourMapProjection::RESULTS_TOP_VALUE, "TOP_VALUE", "Top Value");
addItem(RimContourMapProjection::RESULTS_MEAN_VALUE, "MEAN_VALUE", "Arithmetic Mean");
addItem(RimContourMapProjection::RESULTS_HARM_VALUE, "HARM_VALUE", "Harmonic Mean");
addItem(RimContourMapProjection::RESULTS_GEOM_VALUE, "GEOM_VALUE", "Geometric Mean");
addItem(RimContourMapProjection::RESULTS_MIN_VALUE, "MIN_VALUE", "Min Value");
addItem(RimContourMapProjection::RESULTS_MAX_VALUE, "MAX_VALUE", "Max Value");
addItem(RimContourMapProjection::RESULTS_VOLUME_SUM, "VOLUME_SUM", "Volume Weighted Sum");
addItem(RimContourMapProjection::RESULTS_SUM, "SUM", "Sum");
addItem(RimContourMapProjection::RESULTS_OIL_COLUMN, "OIL_COLUMN", "Oil Column");
addItem(RimContourMapProjection::RESULTS_GAS_COLUMN, "GAS_COLUMN", "Gas Column");
addItem(RimContourMapProjection::RESULTS_HC_COLUMN, "HC_COLUMN", "Hydrocarbon Column");
setDefault(RimContourMapProjection::RESULTS_MEAN_VALUE);
}
}
CAF_PDM_SOURCE_INIT(RimContourMapProjection, "RimContourMapProjection");
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimContourMapProjection::RimContourMapProjection()
{
CAF_PDM_InitObject("RimContourMapProjection", ":/draw_style_meshlines_24x24.png", "", "");
CAF_PDM_InitField(&m_relativeSampleSpacing, "SampleSpacing", 0.75, "Sample Spacing Factor", "", "", "");
m_relativeSampleSpacing.uiCapability()->setUiEditorTypeName(caf::PdmUiDoubleSliderEditor::uiEditorTypeName());
CAF_PDM_InitFieldNoDefault(&m_resultAggregation, "ResultAggregation", "Result Aggregation", "", "", "");
CAF_PDM_InitField(&m_showContourLines, "ContourLines", true, "Show Contour Lines", "", "", "");
CAF_PDM_InitField(&m_weightByParameter, "WeightByParameter", false, "Weight by Result Parameter", "", "", "");
CAF_PDM_InitFieldNoDefault(&m_weightingResult, "WeightingResult", "", "", "", "");
m_weightingResult.uiCapability()->setUiHidden(true);
m_weightingResult.uiCapability()->setUiTreeChildrenHidden(true);
m_weightingResult = new RimEclipseResultDefinition;
m_weightingResult->findField("MResultType")->uiCapability()->setUiName("Result Type");
setName("Map Projection");
nameField()->uiCapability()->setUiReadOnly(true);
m_resultAccessor = new RigHugeValResultAccessor;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimContourMapProjection::~RimContourMapProjection()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::BoundingBox RimContourMapProjection::expandedBoundingBox() const
{
cvf::BoundingBox boundingBox = eclipseCase()->activeCellsBoundingBox();
boundingBox.expand(sampleSpacing() * 0.5);
return boundingBox;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::generateGridMapping()
{
calculateTotalCellVisibility();
cvf::Vec3d gridExtent = expandedBoundingBox().extent();
cvf::Vec2ui gridSize2d = surfaceGridSize();
RimEclipseResultCase* eclipseCase = nullptr;
firstAncestorOrThisOfTypeAsserted(eclipseCase);
m_projected3dGridIndices.resize(vertexCount());
int nVertices = vertexCount();
const std::vector<double>* weightingResultValues = nullptr;
if (m_weightByParameter())
{
size_t gridScalarResultIdx = m_weightingResult->scalarResultIndex();
if (gridScalarResultIdx != cvf::UNDEFINED_SIZE_T)
{
m_weightingResult->loadResult();
int timeStep = 0;
if (m_weightingResult->hasDynamicResult())
{
timeStep = view()->currentTimeStep();
}
weightingResultValues = &(m_weightingResult->currentGridCellResults()->cellScalarResults(gridScalarResultIdx)[timeStep]);
}
}
if (isStraightSummationResult())
{
for (int index = 0; index < nVertices; ++index)
{
cvf::Vec2ui ij = ijFromGridIndex(index);
cvf::Vec2d globalPos = globalPos2d(ij.x(), ij.y());
m_projected3dGridIndices[index] = visibleCellsAndLengthInCellFrom2dPoint(globalPos, weightingResultValues);
}
}
else
{
#pragma omp parallel for
for (int index = 0; index < nVertices; ++index)
{
cvf::Vec2ui ij = ijFromGridIndex(index);
cvf::Vec2d globalPos = globalPos2d(ij.x(), ij.y());
m_projected3dGridIndices[index] = visibleCellsAndOverlapVolumeFrom2dPoint(globalPos, weightingResultValues);
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::generateVertices(cvf::Vec3fArray* vertices, const caf::DisplayCoordTransform* displayCoordTransform)
{
CVF_ASSERT(vertices);
vertices->resize(vertexCount());
cvf::Vec2ui gridSize2d = surfaceGridSize();
cvf::BoundingBox boundingBox = expandedBoundingBox();
int nVertices = vertexCount();
#pragma omp parallel for
for (int index = 0; index < nVertices; ++index)
{
cvf::Vec2ui ij = ijFromGridIndex(index);
cvf::Vec2d globalPos = globalPos2d(ij.x(), ij.y());
cvf::Vec3d globalVertexPos(globalPos, boundingBox.min().z() - 1.0);
cvf::Vec3f displayVertexPos(displayCoordTransform->transformToDisplayCoord(globalVertexPos));
(*vertices)[index] = displayVertexPos;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimContourMapProjection::ContourPolygons RimContourMapProjection::generateContourPolygons(const caf::DisplayCoordTransform* displayCoordTransform)
{
std::vector<cvf::ref<cvf::Vec3fArray>> contourPolygons;
if (minValue() != std::numeric_limits<double>::infinity() &&
maxValue() != -std::numeric_limits<double>::infinity() &&
std::fabs(maxValue() - minValue()) > 1.0e-8)
{
cvf::BoundingBox boundingBox = expandedBoundingBox();
std::vector<double> contourLevels;
legendConfig()->scalarMapper()->majorTickValues(&contourLevels);
int nContourLevels = static_cast<int>(contourLevels.size());
if (nContourLevels > 2)
{
contourLevels[0] += (contourLevels[1] - contourLevels[0]) * 0.01;
contourLevels[nContourLevels - 1] -= (contourLevels[nContourLevels - 1] - contourLevels[nContourLevels - 2]) * 0.01;
std::vector<std::vector<cvf::Vec2d>> contourLines;
caf::ContourLines::create(m_aggregatedResults, xPositions(), yPositions(), contourLevels, &contourLines);
contourPolygons.reserve(contourLines.size());
for (size_t i = 0; i < contourLines.size(); ++i)
{
if (!contourLines[i].empty())
{
cvf::ref<cvf::Vec3fArray> contourPolygon = new cvf::Vec3fArray(contourLines[i].size());
for (size_t j = 0; j < contourLines[i].size(); ++j)
{
cvf::Vec3d contourPoint3d = cvf::Vec3d(contourLines[i][j], boundingBox.min().z());
cvf::Vec3d displayPoint3d = displayCoordTransform->transformToDisplayCoord(contourPoint3d);
(*contourPolygon)[j] = cvf::Vec3f(displayPoint3d);
}
contourPolygons.push_back(contourPolygon);
}
}
}
}
return contourPolygons;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::generateResults()
{
generateGridMapping();
int nVertices = vertexCount();
m_aggregatedResults = std::vector<double>(nVertices, std::numeric_limits<double>::infinity());
int timeStep = view()->currentTimeStep();
RimEclipseCellColors* cellColors = view()->cellResult();
RimEclipseResultCase* eclipseCase = nullptr;
firstAncestorOrThisOfTypeAsserted(eclipseCase);
{
if (!cellColors->isTernarySaturationSelected())
{
RigCaseCellResultsData* resultData = eclipseCase->results(RiaDefines::MATRIX_MODEL);
if (isColumnResult())
{
resultData->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "PORO");
resultData->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "NTG");
resultData->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "DZ");
if (m_resultAggregation == RESULTS_OIL_COLUMN || m_resultAggregation == RESULTS_HC_COLUMN)
{
resultData->findOrLoadScalarResultForTimeStep(RiaDefines::DYNAMIC_NATIVE, "SOIL", timeStep);
}
if (m_resultAggregation == RESULTS_GAS_COLUMN || m_resultAggregation == RESULTS_HC_COLUMN)
{
resultData->findOrLoadScalarResultForTimeStep(RiaDefines::DYNAMIC_NATIVE, "SGAS", timeStep);
}
}
else
{
m_resultAccessor = RigResultAccessorFactory::createFromResultDefinition(eclipseCase->eclipseCaseData(), 0, timeStep, cellColors);
if (m_resultAccessor.isNull())
{
m_resultAccessor = new RigHugeValResultAccessor;
}
}
#pragma omp parallel for
for (int index = 0; index < nVertices; ++index)
{
cvf::Vec2ui ij = ijFromGridIndex(index);
m_aggregatedResults[index] = calculateValue(ij.x(), ij.y());
}
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::maxValue() const
{
double maxV = -std::numeric_limits<double>::infinity();
int nVertices = vertexCount();
for (int index = 0; index < nVertices; ++index)
{
if (m_aggregatedResults[index] != std::numeric_limits<double>::infinity())
{
maxV = std::max(maxV, m_aggregatedResults[index]);
}
}
return maxV;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::minValue() const
{
double minV = std::numeric_limits<double>::infinity();
int nVertices = vertexCount();
for (int index = 0; index < nVertices; ++index)
{
if (m_aggregatedResults[index] != std::numeric_limits<double>::infinity())
{
minV = std::min(minV, m_aggregatedResults[index]);
}
}
return minV;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::meanValue() const
{
return sumAllValues() / validVertexCount();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::sumAllValues() const
{
double sum = 0.0;
int nVertices = vertexCount();
for (int index = 0; index < nVertices; ++index)
{
if (m_aggregatedResults[index] != std::numeric_limits<double>::infinity())
{
sum += m_aggregatedResults[index];
}
}
return sum;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::sampleSpacing() const
{
return m_relativeSampleSpacing * mainGrid()->characteristicIJCellSize();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::sampleSpacingFactor() const
{
return m_relativeSampleSpacing();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::showContourLines() const
{
return m_showContourLines();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RimContourMapProjection::aggregatedResults() const
{
return m_aggregatedResults;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RimContourMapProjection::weightingParameter() const
{
QString parameter = "None";
if (m_weightByParameter() && !m_weightingResult->isTernarySaturationSelected())
{
parameter = m_weightingResult->resultVariableUiShortName();
}
return parameter;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::isMeanResult() const
{
return m_resultAggregation() == RESULTS_MEAN_VALUE ||
m_resultAggregation() == RESULTS_HARM_VALUE ||
m_resultAggregation() == RESULTS_GEOM_VALUE;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::isSummationResult() const
{
return isStraightSummationResult() || m_resultAggregation() == RESULTS_VOLUME_SUM;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::isStraightSummationResult() const
{
return isColumnResult() || m_resultAggregation() == RESULTS_SUM;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::isColumnResult() const
{
return m_resultAggregation() == RESULTS_OIL_COLUMN ||
m_resultAggregation() == RESULTS_GAS_COLUMN ||
m_resultAggregation() == RESULTS_HC_COLUMN;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::value(uint i, uint j) const
{
size_t index = gridIndex(i, j);
if (index < vertexCount())
{
return m_aggregatedResults.at(gridIndex(i, j));
}
return std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::calculateValue(uint i, uint j) const
{
if (!isColumnResult() && view()->cellResult()->scalarResultIndex() == cvf::UNDEFINED_SIZE_T)
{
return 0.0; // Special case of NONE-result. Show 0 all over to ensure we see something.
}
const std::vector<std::pair<size_t, double>>& matchingCells = cellsAtPos2d(i, j);
if (!matchingCells.empty())
{
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<double> calculator;
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = m_resultAccessor->cellScalarGlobIdx(cellIdx);
calculator.addValueAndWeight(cellValue, cellIdxAndWeight.second);
}
if (calculator.validAggregatedWeight())
{
return calculator.weightedMean();
}
return std::numeric_limits<double>::infinity();
}
case RESULTS_GEOM_VALUE:
{
RiaWeightedGeometricMeanCalculator calculator;
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = m_resultAccessor->cellScalarGlobIdx(cellIdx);
if (cellValue < 1.0e-8)
{
return 0.0;
}
calculator.addValueAndWeight(cellValue, cellIdxAndWeight.second);
}
if (calculator.validAggregatedWeight())
{
qDebug() << calculator.weightedMean();
return calculator.weightedMean();
}
return std::numeric_limits<double>::infinity();
}
case RESULTS_HARM_VALUE:
{
RiaWeightedHarmonicMeanCalculator calculator;
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = m_resultAccessor->cellScalarGlobIdx(cellIdx);
if (std::fabs(cellValue) < 1.0e-8)
{
return 0.0;
}
calculator.addValueAndWeight(cellValue, cellIdxAndWeight.second);
}
if (calculator.validAggregatedWeight())
{
return calculator.weightedMean();
}
return std::numeric_limits<double>::infinity();
}
case RESULTS_MAX_VALUE:
{
double maxValue = -std::numeric_limits<double>::infinity();
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = m_resultAccessor->cellScalarGlobIdx(cellIdx);
maxValue = std::max(maxValue, cellValue);
}
return maxValue;
}
case RESULTS_MIN_VALUE:
{
double minValue = std::numeric_limits<double>::infinity();
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = m_resultAccessor->cellScalarGlobIdx(cellIdx);
minValue = std::min(minValue, cellValue);
}
return minValue;
}
case RESULTS_VOLUME_SUM:
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;
}
case RESULTS_OIL_COLUMN:
case RESULTS_GAS_COLUMN:
case RESULTS_HC_COLUMN:
{
double sum = 0.0;
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = findColumnResult(m_resultAggregation(), cellIdx);
sum += cellValue * cellIdxAndWeight.second;
}
return sum;
}
default:
CVF_TIGHT_ASSERT(false);
}
}
return std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::hasResultAt(uint i, uint j) const
{
RimEclipseCellColors* cellColors = view()->cellResult();
if (cellColors->isTernarySaturationSelected())
{
return false;
}
return !cellsAtPos2d(i, j).empty();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2ui RimContourMapProjection::surfaceGridSize() const
{
cvf::BoundingBox boundingBox = expandedBoundingBox();
cvf::Vec3d gridExtent = boundingBox.extent();
uint projectionSizeX = static_cast<uint>(std::ceil(gridExtent.x() / sampleSpacing())) + 1u;
uint projectionSizeY = static_cast<uint>(std::ceil(gridExtent.y() / sampleSpacing())) + 1u;
return cvf::Vec2ui(projectionSizeX, projectionSizeY);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
uint RimContourMapProjection::vertexCount() const
{
cvf::Vec2ui gridSize2d = surfaceGridSize();
return gridSize2d.x() * gridSize2d.y();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
uint RimContourMapProjection::validVertexCount() const
{
uint validCount = 0u;
for (uint i = 0; i < vertexCount(); ++i)
{
cvf::Vec2ui ij = ijFromGridIndex(i);
if (hasResultAt(ij.x(), ij.y()))
{
validCount++;
}
}
return validCount;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimRegularLegendConfig* RimContourMapProjection::legendConfig() const
{
return view()->cellResult()->legendConfig();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::calculateTotalCellVisibility()
{
m_cellGridIdxVisibility = view()->currentTotalCellVisibility();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2d RimContourMapProjection::globalPos2d(uint i, uint j) const
{
cvf::BoundingBox boundingBox = expandedBoundingBox();
cvf::Vec3d gridExtent = boundingBox.extent();
cvf::Vec2d origin(boundingBox.min().x(), boundingBox.min().y());
cvf::Vec2ui gridSize2d = surfaceGridSize();
return origin + cvf::Vec2d((i * gridExtent.x()) / (gridSize2d.x() - 1),
(j * gridExtent.y()) / (gridSize2d.y() - 1));
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2ui RimContourMapProjection::ijFromLocalPos(const cvf::Vec2d& localPos2d) const
{
cvf::Vec2ui ijCoords(localPos2d.x() / sampleSpacing(), localPos2d.y() / sampleSpacing());
return ijCoords;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<std::pair<size_t, double>>& RimContourMapProjection::cellsAtPos2d(uint i, uint j) const
{
return m_projected3dGridIndices[gridIndex(i, j)];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<size_t, double>> RimContourMapProjection::visibleCellsAndOverlapVolumeFrom2dPoint(const cvf::Vec2d& globalPos2d, const std::vector<double>* weightingResultValues) const
{
cvf::BoundingBox gridBoundingBox = expandedBoundingBox();
cvf::Vec3d top2dElementCentroid(globalPos2d, gridBoundingBox.max().z());
cvf::Vec3d bottom2dElementCentroid(globalPos2d, gridBoundingBox.min().z());
cvf::Vec3d planarDiagonalVector(0.5 * sampleSpacing(), 0.5 * sampleSpacing(), 0.0);
cvf::Vec3d topNECorner = top2dElementCentroid + planarDiagonalVector;
cvf::Vec3d bottomSWCorner = bottom2dElementCentroid - planarDiagonalVector;
cvf::BoundingBox bbox2dElement(bottomSWCorner, topNECorner);
std::vector<size_t> allCellIndices;
mainGrid()->findIntersectingCells(bbox2dElement, &allCellIndices);
typedef std::map<size_t, std::vector<std::pair<size_t, double>>> KLayerCellWeightMap;
KLayerCellWeightMap matchingVisibleCellsWeightPerKLayer;
std::array<cvf::Vec3d, 8> hexCorners;
for (size_t globalCellIdx : allCellIndices)
{
if ((*m_cellGridIdxVisibility)[globalCellIdx])
{
RigCell cell = mainGrid()->globalCellArray()[globalCellIdx];
size_t mainGridCellIdx = cell.mainGridCellIndex();
size_t i, j, k;
mainGrid()->ijkFromCellIndex(mainGridCellIdx, &i, &j, &k);
size_t localCellIdx = cell.gridLocalCellIndex();
RigGridBase* localGrid = cell.hostGrid();
localGrid->cellCornerVertices(localCellIdx, hexCorners.data());
cvf::BoundingBox overlapBBox;
std::array<cvf::Vec3d, 8> overlapCorners = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(hexCorners, bbox2dElement, &overlapBBox);
double overlapVolume = RigCellGeometryTools::calculateCellVolume(overlapCorners);
if (overlapVolume > 0.0)
{
double weight = overlapVolume;
if (weightingResultValues)
{
const RigActiveCellInfo* activeCellInfo = eclipseCase()->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL);
size_t cellResultIdx = activeCellInfo->cellResultIndex(globalCellIdx);
double result = std::max((*weightingResultValues)[cellResultIdx], 0.0);
if (result < 1.0e-6)
{
result = 0.0;
}
weight *= result;
}
if (weight > 0.0)
{
matchingVisibleCellsWeightPerKLayer[k].push_back(std::make_pair(globalCellIdx, weight));
}
}
}
}
std::vector<std::pair<size_t, double>> matchingVisibleCellsAndWeight;
for (auto kLayerCellWeight : matchingVisibleCellsWeightPerKLayer)
{
for (auto cellWeight : kLayerCellWeight.second)
{
matchingVisibleCellsAndWeight.push_back(std::make_pair(cellWeight.first, cellWeight.second));
}
}
return matchingVisibleCellsAndWeight;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<size_t, double>> RimContourMapProjection::visibleCellsAndLengthInCellFrom2dPoint(const cvf::Vec2d& globalPos2d, const std::vector<double>* weightingResultValues /*= nullptr*/) const
{
cvf::BoundingBox boundingBox = expandedBoundingBox();
cvf::Vec3d highestPoint(globalPos2d, boundingBox.max().z());
cvf::Vec3d lowestPoint(globalPos2d, boundingBox.min().z());
cvf::BoundingBox rayBBox;
rayBBox.add(highestPoint);
rayBBox.add(lowestPoint);
std::vector<size_t> allCellIndices;
mainGrid()->findIntersectingCells(rayBBox, &allCellIndices);
std::map<size_t, std::vector<std::pair<size_t, double>>> matchingVisibleCellsAndWeightPerKLayer;
cvf::Vec3d hexCorners[8];
for (size_t globalCellIdx : allCellIndices)
{
if ((*m_cellGridIdxVisibility)[globalCellIdx])
{
RigCell cell = mainGrid()->globalCellArray()[globalCellIdx];
size_t mainGridCellIdx = cell.mainGridCellIndex();
size_t i, j, k;
mainGrid()->ijkFromCellIndex(mainGridCellIdx, &i, &j, &k);
size_t localCellIdx = cell.gridLocalCellIndex();
RigGridBase* localGrid = cell.hostGrid();
localGrid->cellCornerVertices(localCellIdx, hexCorners);
std::vector<HexIntersectionInfo> intersections;
if (RigHexIntersectionTools::lineHexCellIntersection(highestPoint, lowestPoint, hexCorners, 0, &intersections))
{
double lengthInCell = (intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint).length();
matchingVisibleCellsAndWeightPerKLayer[k].push_back(std::make_pair(globalCellIdx, lengthInCell));
}
}
}
std::vector<std::pair<size_t, double>> matchingVisibleCellsAndWeight;
for (auto kLayerCellWeight : matchingVisibleCellsAndWeightPerKLayer)
{
// Make sure the sum of all weights in the same K-layer is 1.
double weightSumThisKLayer = 0.0;
for (auto cellWeight : kLayerCellWeight.second)
{
weightSumThisKLayer += cellWeight.second;
}
for (auto cellWeight : kLayerCellWeight.second)
{
matchingVisibleCellsAndWeight.push_back(std::make_pair(cellWeight.first, cellWeight.second / weightSumThisKLayer));
}
}
return matchingVisibleCellsAndWeight;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::findColumnResult(ResultAggregation resultAggregation, size_t cellGlobalIdx) const
{
const RigCaseCellResultsData* resultData = eclipseCase()->results(RiaDefines::MATRIX_MODEL);
size_t poroResultIndex = resultData->findScalarResultIndex(RiaDefines::STATIC_NATIVE, "PORO");
size_t ntgResultIndex = resultData->findScalarResultIndex(RiaDefines::STATIC_NATIVE, "NTG");
size_t dzResultIndex = resultData->findScalarResultIndex(RiaDefines::STATIC_NATIVE, "DZ");
if (poroResultIndex == cvf::UNDEFINED_SIZE_T || ntgResultIndex == cvf::UNDEFINED_SIZE_T)
{
return std::numeric_limits<double>::infinity();
}
const std::vector<double>& poroResults = resultData->cellScalarResults(poroResultIndex)[0];
const std::vector<double>& ntgResults = resultData->cellScalarResults(ntgResultIndex)[0];
const std::vector<double>& dzResults = resultData->cellScalarResults(dzResultIndex)[0];
const RigActiveCellInfo* activeCellInfo = eclipseCase()->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL);
size_t cellResultIdx = activeCellInfo->cellResultIndex(cellGlobalIdx);
if (cellResultIdx >= poroResults.size() || cellResultIdx >= ntgResults.size())
{
return std::numeric_limits<double>::infinity();
}
double poro = poroResults.at(cellResultIdx);
double ntg = ntgResults.at(cellResultIdx);
double dz = dzResults.at(cellResultIdx);
int timeStep = view()->currentTimeStep();
double resultValue = 0.0;
if (resultAggregation == RESULTS_OIL_COLUMN || resultAggregation == RESULTS_HC_COLUMN)
{
size_t soilResultIndex = resultData->findScalarResultIndex(RiaDefines::DYNAMIC_NATIVE, "SOIL");
const std::vector<double>& soilResults = resultData->cellScalarResults(soilResultIndex)[timeStep];
if (cellResultIdx < soilResults.size())
{
resultValue = soilResults.at(cellResultIdx);
}
}
if (resultAggregation == RESULTS_GAS_COLUMN || resultAggregation == RESULTS_HC_COLUMN)
{
size_t sgasResultIndex = resultData->findScalarResultIndex(RiaDefines::DYNAMIC_NATIVE, "SGAS");
const std::vector<double>& sgasResults = resultData->cellScalarResults(sgasResultIndex)[timeStep];
if (cellResultIdx < sgasResults.size())
{
resultValue += sgasResults.at(cellResultIdx);
}
}
return resultValue * poro * ntg * dz;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const RimEclipseResultCase* RimContourMapProjection::eclipseCase() const
{
const RimEclipseResultCase* eclipseCase = nullptr;
firstAncestorOrThisOfType(eclipseCase);
return eclipseCase;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimEclipseResultCase* RimContourMapProjection::eclipseCase()
{
RimEclipseResultCase* eclipseCase = nullptr;
firstAncestorOrThisOfType(eclipseCase);
return eclipseCase;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimContourMapView* RimContourMapProjection::view() const
{
RimContourMapView* view = nullptr;
firstAncestorOrThisOfTypeAsserted(view);
return view;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RimContourMapProjection::gridIndex(uint i, uint j) const
{
cvf::Vec2ui gridSize2d = surfaceGridSize();
CVF_ASSERT(i < gridSize2d.x());
CVF_ASSERT(j < gridSize2d.y());
return i + j * gridSize2d.x();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec2ui RimContourMapProjection::ijFromGridIndex(size_t index) const
{
CVF_TIGHT_ASSERT(index < vertexCount());
cvf::Vec2ui gridSize2d = surfaceGridSize();
uint quotientX = static_cast<uint>(index) / gridSize2d.x();
uint remainderX = static_cast<uint>(index) % gridSize2d.x();
return cvf::Vec2ui(remainderX, quotientX);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::updateLegend()
{
RimEclipseCellColors* cellColors = view()->cellResult();
if (getLegendRangeFrom3dGrid())
{
cellColors->updateLegendData(view()->currentTimeStep(), legendConfig());
}
else
{
double minVal = minValue();
double maxVal = maxValue();
legendConfig()->setAutomaticRanges(minVal, maxVal, minVal, maxVal);
}
if (m_resultAggregation() == RESULTS_OIL_COLUMN ||
m_resultAggregation() == RESULTS_GAS_COLUMN ||
m_resultAggregation() == RESULTS_HC_COLUMN)
{
legendConfig()->setTitle(QString("Map Projection\n%1").arg(m_resultAggregation().uiText()));
}
else
{
QString projectionLegendText = QString("Map Projection\n%1").arg(m_resultAggregation().uiText());
if (weightingParameter() != "None")
{
projectionLegendText += QString("(W: %1)").arg(weightingParameter());
}
projectionLegendText += QString("\nResult: %1").arg(cellColors->resultVariableUiShortName());
legendConfig()->setTitle(projectionLegendText);
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimContourMapProjection::ResultAggregation RimContourMapProjection::resultAggregation() const
{
return m_resultAggregation();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RimContourMapProjection::resultAggregationText() const
{
return m_resultAggregation().uiText();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RimContourMapProjection::resultDescriptionText() const
{
QString resultText = resultAggregationText();
if (!isColumnResult())
{
resultText += QString(", %1").arg(view()->cellResult()->resultVariable());
}
return resultText;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::updatedWeightingResult()
{
this->updateConnectedEditors();
this->generateResults();
this->updateLegend();
RimProject* proj;
this->firstAncestorOrThisOfTypeAsserted(proj);
proj->scheduleCreateDisplayModelAndRedrawAllViews();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::checkForMapIntersection(const cvf::Vec3d& localPoint3d, cvf::Vec2ui* contourMapIJ, double* valueAtPoint) const
{
CVF_TIGHT_ASSERT(contourMapIJ);
CVF_TIGHT_ASSERT(valueAtPoint);
cvf::Vec2d localPos2d(localPoint3d.x(), localPoint3d.y());
*contourMapIJ = ijFromLocalPos(localPos2d);
*valueAtPoint = value(contourMapIJ->x(), contourMapIJ->y());
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RimContourMapProjection::xPositions() const
{
cvf::BoundingBox boundingBox = expandedBoundingBox();
cvf::Vec3d gridExtent = boundingBox.extent();
double origin = boundingBox.min().x();
cvf::Vec2ui gridSize2d = surfaceGridSize();
std::vector<double> positions;
positions.reserve(gridSize2d.x());
for (uint i = 0; i < gridSize2d.x(); ++i)
{
positions.push_back(origin + (i * gridExtent.x()) / (gridSize2d.x() - 1));
}
return positions;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RimContourMapProjection::yPositions() const
{
cvf::BoundingBox boundingBox = expandedBoundingBox();
cvf::Vec3d gridExtent = boundingBox.extent();
double origin = boundingBox.min().y();
cvf::Vec2ui gridSize2d = surfaceGridSize();
std::vector<double> positions;
positions.reserve(gridSize2d.y());
for (uint j = 0; j < gridSize2d.y(); ++j)
{
positions.push_back(origin + (j * gridExtent.y()) / (gridSize2d.y() - 1));
}
return positions;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigMainGrid* RimContourMapProjection::mainGrid() const
{
RimEclipseResultCase* eclipseCase = nullptr;
firstAncestorOrThisOfTypeAsserted(eclipseCase);
return eclipseCase->eclipseCaseData()->mainGrid();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::fieldChangedByUi(const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue)
{
legendConfig()->disableAllTimeStepsRange(!getLegendRangeFrom3dGrid());
m_weightingResult->loadResult();
view()->updateConnectedEditors();
RimProject* proj;
firstAncestorOrThisOfTypeAsserted(proj);
proj->scheduleCreateDisplayModelAndRedrawAllViews();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::defineEditorAttribute(const caf::PdmFieldHandle* field, QString uiConfigName, caf::PdmUiEditorAttribute* attribute)
{
if (&m_relativeSampleSpacing == field)
{
caf::PdmUiDoubleSliderEditorAttribute* myAttr = dynamic_cast<caf::PdmUiDoubleSliderEditorAttribute*>(attribute);
if (myAttr)
{
myAttr->m_minimum = 0.25;
myAttr->m_maximum = 2.0;
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::defineUiOrdering(QString uiConfigName, caf::PdmUiOrdering& uiOrdering)
{
caf::PdmUiGroup* mainGroup = uiOrdering.addNewGroup("Projection Settings");
mainGroup->add(&m_relativeSampleSpacing);
mainGroup->add(&m_resultAggregation);
mainGroup->add(&m_showContourLines);
caf::PdmUiGroup* weightingGroup = uiOrdering.addNewGroup("Mean Weighting Options");
weightingGroup->add(&m_weightByParameter);
weightingGroup->setCollapsedByDefault(true);
m_weightByParameter.uiCapability()->setUiReadOnly(!isMeanResult());
if (!isMeanResult())
{
m_weightByParameter = false;
}
if (m_weightByParameter())
{
m_weightingResult->uiOrdering(uiConfigName, *weightingGroup);
}
uiOrdering.skipRemainingFields(true);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::defineUiTreeOrdering(caf::PdmUiTreeOrdering& uiTreeOrdering, QString uiConfigName /*= ""*/)
{
uiTreeOrdering.skipRemainingChildren(true);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::initAfterRead()
{
legendConfig()->disableAllTimeStepsRange(!getLegendRangeFrom3dGrid());
if (eclipseCase())
{
m_weightingResult->setEclipseCase(eclipseCase());
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimContourMapProjection::getLegendRangeFrom3dGrid() const
{
if (isMeanResult())
{
return true;
}
else if (m_resultAggregation == RESULTS_TOP_VALUE)
{
return true;
}
return false;
}
Reference in New Issue View Git Blame Copy Permalink