#3980 Optimizing and expanding geomech boundaries to outside Pore pressure region

This commit is contained in:
Gaute Lindkvist 2019-01-21 15:45:35 +01:00
parent 810a991ba4
commit 9a3a61811f
8 changed files with 487 additions and 307 deletions

View File

@ -462,7 +462,13 @@ bool RigFemPartGrid::cellIJKFromCoordinate(const cvf::Vec3d& coord, size_t* i, s
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::cellCornerVertices(size_t cellIndex, cvf::Vec3d vertices[8]) const
{
CVF_ASSERT(false);
const std::vector<cvf::Vec3f>& nodeCoords = m_femPart->nodes().coordinates;
const int* cornerIndices = m_femPart->connectivities(cellIndex);
for (size_t i = 0; i < 8; ++i)
{
vertices[i] = cvf::Vec3d(nodeCoords[cornerIndices[i]]);
}
}
//--------------------------------------------------------------------------------------------------

View File

@ -40,6 +40,7 @@ public:
cvf::Vec3i findMainIJKFaces(int elementIndex) const;
std::pair<cvf::Vec3st, cvf::Vec3st> reservoirIJKBoundingBox() const;
void cellCornerVertices(size_t cellIndex, cvf::Vec3d vertices[8]) const override;
private:
void generateStructGridData();
@ -59,7 +60,6 @@ private: // Unused, Not implemented
bool cellIJKFromCoordinate(const cvf::Vec3d& coord, size_t* i, size_t* j, size_t* k) const override;
void cellCornerVertices(size_t cellIndex, cvf::Vec3d vertices[8]) const override;
cvf::Vec3d cellCentroid(size_t cellIndex) const override;
void cellMinMaxCordinates(size_t cellIndex, cvf::Vec3d* minCoordinate, cvf::Vec3d* maxCoordinate) const override;
size_t gridPointIndexFromIJK(size_t i, size_t j, size_t k) const override;

View File

@ -109,9 +109,11 @@ void RimContourMapProjection::generateResultsIfNecessary(int timeStep)
{
updateGridInformation();
if (gridMappingNeedsUpdating())
if (gridMappingNeedsUpdating() || mapCellVisibilityNeedsUpdating())
{
generateGridMapping();
clearResults();
m_projected3dGridIndices = generateGridMapping();
m_mapCellVisibility = getMapCellVisibility();
}
if (resultVariableChanged())
@ -417,7 +419,15 @@ cvf::Vec3d RimContourMapProjection::origin3d() const
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimContourMapProjection::calculateValueInMapCell(uint i, uint j) const
size_t RimContourMapProjection::gridResultIndex(size_t globalCellIdx) const
{
return globalCellIdx;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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);
if (!matchingCells.empty())
@ -426,9 +436,16 @@ double RimContourMapProjection::calculateValueInMapCell(uint i, uint j) const
{
case RESULTS_TOP_VALUE:
{
size_t cellIdx = matchingCells.front().first;
double cellValue = gridCellValue(cellIdx);
return cellValue;
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = gridCellValues[cellIdx];
if (cellValue != std::numeric_limits<double>::infinity())
{
return cellValue;
}
}
return std::numeric_limits<double>::infinity();
}
case RESULTS_MEAN_VALUE:
{
@ -436,7 +453,7 @@ double RimContourMapProjection::calculateValueInMapCell(uint i, uint j) const
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = gridCellValue(cellIdx);
double cellValue = gridCellValues[cellIdx];
if (cellValue != std::numeric_limits<double>::infinity())
{
calculator.addValueAndWeight(cellValue, cellIdxAndWeight.second);
@ -454,7 +471,7 @@ double RimContourMapProjection::calculateValueInMapCell(uint i, uint j) const
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = gridCellValue(cellIdx);
double cellValue = gridCellValues[cellIdx];
if (cellValue < 1.0e-8)
{
return 0.0;
@ -476,7 +493,7 @@ double RimContourMapProjection::calculateValueInMapCell(uint i, uint j) const
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = gridCellValue(cellIdx);
double cellValue = gridCellValues[cellIdx];
if (std::fabs(cellValue) < 1.0e-8)
{
return 0.0;
@ -498,7 +515,7 @@ double RimContourMapProjection::calculateValueInMapCell(uint i, uint j) const
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = gridCellValue(cellIdx);
double cellValue = gridCellValues[cellIdx];
if (cellValue != std::numeric_limits<double>::infinity())
{
maxValue = std::max(maxValue, cellValue);
@ -512,19 +529,22 @@ double RimContourMapProjection::calculateValueInMapCell(uint i, uint j) const
for (auto cellIdxAndWeight : matchingCells)
{
size_t cellIdx = cellIdxAndWeight.first;
double cellValue = gridCellValue(cellIdx);
double cellValue = gridCellValues[cellIdx];
minValue = std::min(minValue, cellValue);
}
return minValue;
}
case RESULTS_VOLUME_SUM:
case RESULTS_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 = gridCellValue(cellIdx);
double cellValue = gridCellValues[cellIdx];
if (cellValue != std::numeric_limits<double>::infinity())
{
sum += cellValue * cellIdxAndWeight.second;
@ -532,7 +552,6 @@ double RimContourMapProjection::calculateValueInMapCell(uint i, uint j) const
}
return sum;
}
default:
CVF_TIGHT_ASSERT(false);
}
@ -561,6 +580,7 @@ bool RimContourMapProjection::gridMappingNeedsUpdating() const
{
if ((*currentVisibility)[i] != (*m_cellGridIdxVisibility)[i]) return true;
}
return false;
}
@ -615,6 +635,7 @@ void RimContourMapProjection::clearGridMapping()
clearResults();
clearTimeStepRange();
m_projected3dGridIndices.clear();
m_mapCellVisibility.clear();
}
//--------------------------------------------------------------------------------------------------
@ -677,20 +698,22 @@ double RimContourMapProjection::minValue(const std::vector<double>& aggregatedRe
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<double, double> RimContourMapProjection::minmaxValuesAllTimeSteps(int skipSteps /*= 1*/)
std::pair<double, double> RimContourMapProjection::minmaxValuesAllTimeSteps()
{
if (timestepRangeNeedsUpdating())
{
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 += skipSteps)
for (int i = 0; i < (int)baseView()->ownerCase()->timeStepStrings().size() - 1; ++i)
{
std::vector<double> aggregatedResults = generateResults(i);
if (i == m_currentResultTimestep)
{
std::vector<double> aggregatedResults = generateResults(i, true);
m_minResultAllTimeSteps = std::min(m_minResultAllTimeSteps, minValue(m_aggregatedResults));
m_maxResultAllTimeSteps = std::max(m_maxResultAllTimeSteps, maxValue(m_aggregatedResults));
}
else
{
std::vector<double> aggregatedResults = generateResults(i);
m_minResultAllTimeSteps = std::min(m_minResultAllTimeSteps, minValue(aggregatedResults));
m_maxResultAllTimeSteps = std::max(m_maxResultAllTimeSteps, maxValue(aggregatedResults));
}
@ -707,17 +730,32 @@ cvf::ref<cvf::UByteArray> RimContourMapProjection::getCellVisibility() const
return baseView()->currentTotalCellVisibility();
}
//--------------------------------------------------------------------------------------------------
/// Empty default implementation
//--------------------------------------------------------------------------------------------------
std::vector<bool> RimContourMapProjection::getMapCellVisibility()
{
return std::vector<bool>(numberOfCells(), true);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimContourMapProjection::generateGridMapping()
bool RimContourMapProjection::mapCellVisibilityNeedsUpdating()
{
clearResults();
std::vector<bool> mapCellVisiblity = getMapCellVisibility();
return !(mapCellVisiblity == m_mapCellVisibility);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<std::pair<size_t, double>>> RimContourMapProjection::generateGridMapping()
{
m_cellGridIdxVisibility = getCellVisibility();
int nCells = numberOfCells();
m_projected3dGridIndices.resize(nCells);
std::vector<std::vector<std::pair<size_t, double>>> projected3dGridIndices(nCells);
std::vector<double> weightingResultValues = retrieveParameterWeights();
@ -728,8 +766,8 @@ void RimContourMapProjection::generateGridMapping()
{
cvf::Vec2ui ij = ijFromCellIndex(index);
cvf::Vec2d globalPos = cellCenterPosition(ij.x(), ij.y()) + origin2d();
m_projected3dGridIndices[index] = cellRayIntersectionAndResults(globalPos, weightingResultValues);
cvf::Vec2d globalPos = cellCenterPosition(ij.x(), ij.y()) + origin2d();
projected3dGridIndices[index] = cellRayIntersectionAndResults(globalPos, weightingResultValues);
}
}
else
@ -739,10 +777,12 @@ void RimContourMapProjection::generateGridMapping()
{
cvf::Vec2ui ij = ijFromCellIndex(index);
cvf::Vec2d globalPos = cellCenterPosition(ij.x(), ij.y()) + origin2d();
m_projected3dGridIndices[index] = cellOverlapVolumesAndResults(globalPos, weightingResultValues);
cvf::Vec2d globalPos = cellCenterPosition(ij.x(), ij.y()) + origin2d();
projected3dGridIndices[index] = cellOverlapVolumesAndResults(globalPos, weightingResultValues);
}
}
return projected3dGridIndices;
}
//--------------------------------------------------------------------------------------------------
@ -973,7 +1013,7 @@ void RimContourMapProjection::generateTrianglesWithVertexValues()
}
double triangleAreasThisLevel = sumTriangleAreas(allTrianglesThisLevel);
if (triangleAreasThisLevel > 1.0e-3 * m_contourLevelCumulativeAreas[c])
if (c >= m_contourLevelCumulativeAreas.size() || triangleAreasThisLevel > 1.0e-3 * m_contourLevelCumulativeAreas[c])
{
trianglesPerLevel[c] = allTrianglesThisLevel;
}
@ -1206,9 +1246,9 @@ double RimContourMapProjection::sumTriangleAreas(const std::vector<cvf::Vec4d>&
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<size_t, double>>
std::vector<RimContourMapProjection::CellIndexAndResult>
RimContourMapProjection::cellOverlapVolumesAndResults(const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues) const
const std::vector<double>& weightingResultValues) const
{
cvf::Vec3d top2dElementCentroid(globalPos2d, m_expandedBoundingBox.max().z());
cvf::Vec3d bottom2dElementCentroid(globalPos2d, m_expandedBoundingBox.min().z());
@ -1228,33 +1268,30 @@ std::vector<std::pair<size_t, double>>
std::vector<size_t> allCellIndices = findIntersectingCells(bbox2dElement);
typedef std::map<size_t, std::vector<std::pair<size_t, double>>> KLayerCellWeightMap;
KLayerCellWeightMap matchingVisibleCellsWeightPerKLayer;
std::map<size_t, std::vector<size_t>> kLayerIndexMap;
for (size_t globalCellIdx : allCellIndices)
{
if ((*m_cellGridIdxVisibility)[globalCellIdx])
{
size_t cellKLayer = 0u;
double overlapVolume = calculateOverlapVolume(globalCellIdx, bbox2dElement, &cellKLayer);
if (overlapVolume > 0.0)
{
double weight = overlapVolume * getParameterWeightForCell(globalCellIdx, weightingResultValues);
if (weight > 0.0)
{
matchingVisibleCellsWeightPerKLayer[cellKLayer].push_back(std::make_pair(globalCellIdx, weight));
}
}
kLayerIndexMap[kLayer(globalCellIdx)].push_back(globalCellIdx);
}
}
for (auto kLayerCellWeight : matchingVisibleCellsWeightPerKLayer)
for (const auto& kLayerIndexPair : kLayerIndexMap)
{
for (auto cellWeight : kLayerCellWeight.second)
for (size_t globalCellIdx : kLayerIndexPair.second)
{
matchingVisibleCellsAndWeight.push_back(std::make_pair(cellWeight.first, cellWeight.second));
double overlapVolume = calculateOverlapVolume(globalCellIdx, bbox2dElement);
if (overlapVolume > 0.0)
{
size_t resultIndex = gridResultIndex(globalCellIdx);
double weight = overlapVolume * getParameterWeightForCell(resultIndex, weightingResultValues);
if (weight > 0.0)
{
matchingVisibleCellsAndWeight.push_back(std::make_pair(resultIndex, weight));
}
}
}
}
@ -1264,16 +1301,15 @@ std::vector<std::pair<size_t, double>>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<size_t, double>>
RimContourMapProjection::cellRayIntersectionAndResults(const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues) const
std::vector<RimContourMapProjection::CellIndexAndResult> RimContourMapProjection::cellRayIntersectionAndResults(const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues) const
{
std::vector<std::pair<size_t, double>> matchingVisibleCellsAndWeight;
cvf::Vec3d highestPoint(globalPos2d, m_expandedBoundingBox.max().z());
cvf::Vec3d lowestPoint(globalPos2d, m_expandedBoundingBox.min().z());
// Bounding box has been expanded, so ray may be outside actual grid
// Bounding box has been expanded, so ray may be outside actual 3d grid
if (!m_gridBoundingBox.contains(highestPoint))
{
return matchingVisibleCellsAndWeight;
@ -1285,33 +1321,33 @@ std::vector<std::pair<size_t, double>>
std::vector<size_t> allCellIndices = findIntersectingCells(rayBBox);
std::map<size_t, std::vector<std::pair<size_t, double>>> matchingVisibleCellsAndWeightPerKLayer;
std::map<size_t, std::vector<size_t>> kLayerIndexMap;
for (size_t globalCellIdx : allCellIndices)
{
if ((*m_cellGridIdxVisibility)[globalCellIdx])
{
size_t cellKLayer = 0u;
double lengthInCell = calculateRayLengthInCell(globalCellIdx, highestPoint, lowestPoint, &cellKLayer);
if (lengthInCell > 0.0)
{
matchingVisibleCellsAndWeightPerKLayer[cellKLayer].push_back(std::make_pair(globalCellIdx, lengthInCell));
}
kLayerIndexMap[kLayer(globalCellIdx)].push_back(globalCellIdx);
}
}
for (auto kLayerCellWeight : matchingVisibleCellsAndWeightPerKLayer)
for (const auto& kLayerIndexPair : kLayerIndexMap)
{
// Make sure the sum of all weights in the same K-layer is 1.
double weightSumThisKLayer = 0.0;
for (auto cellWeight : kLayerCellWeight.second)
std::vector<std::pair<size_t, double>> cellsAndWeightsThisLayer;
for (size_t globalCellIdx : kLayerIndexPair.second)
{
weightSumThisKLayer += cellWeight.second;
double lengthInCell = calculateRayLengthInCell(globalCellIdx, highestPoint, lowestPoint);
if (lengthInCell > 0.0)
{
cellsAndWeightsThisLayer.push_back(std::make_pair(gridResultIndex(globalCellIdx), lengthInCell));
weightSumThisKLayer += lengthInCell;
}
}
for (auto cellWeight : kLayerCellWeight.second)
for (auto& cellWeightPair : cellsAndWeightsThisLayer)
{
matchingVisibleCellsAndWeight.push_back(std::make_pair(cellWeight.first, cellWeight.second / weightSumThisKLayer));
cellWeightPair.second /= weightSumThisKLayer;
matchingVisibleCellsAndWeight.push_back(cellWeightPair);
}
}

View File

@ -125,13 +125,15 @@ protected:
virtual bool resultVariableChanged() const = 0;
virtual void clearResultVariable() = 0;
virtual RimGridView* baseView() const = 0;
virtual size_t kLayer(size_t globalCellIdx) const = 0;
virtual std::vector<size_t> findIntersectingCells(const cvf::BoundingBox& bbox) const = 0;
virtual double calculateOverlapVolume(size_t globalCellIdx, const cvf::BoundingBox& bbox, size_t* cellKLayerOut) const = 0;
virtual double calculateRayLengthInCell(size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint, size_t* cellKLayerOut) 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 double gridCellValue(size_t globalCellIdx) const = 0;
virtual double calculateValueInMapCell(uint i, uint j) const;
virtual size_t gridResultIndex(size_t globalCellIdx) const;
double calculateValueInMapCell(uint i, uint j, const std::vector<double>& gridCellValues) const;
protected:
// Keep track of whether cached data needs updating
@ -145,10 +147,13 @@ protected:
double maxValue(const std::vector<double>& aggregatedResults) const;
double minValue(const std::vector<double>& aggregatedResults) const;
std::pair<double, double> minmaxValuesAllTimeSteps(int skipSteps = 1);
std::pair<double, double> minmaxValuesAllTimeSteps();
virtual cvf::ref<cvf::UByteArray> getCellVisibility() const;
virtual std::vector<bool> getMapCellVisibility();
bool mapCellVisibilityNeedsUpdating();
std::vector<std::vector<std::pair<size_t, double>>> generateGridMapping();
virtual cvf::ref<cvf::UByteArray> getCellVisibility() const;
void generateGridMapping();
void generateVertexResults();
void generateTrianglesWithVertexValues();
std::vector<cvf::Vec3d> generateVertices() const;
@ -157,9 +162,9 @@ protected:
static double sumPolygonArea(const ContourPolygons& contourPolygons);
static double sumTriangleAreas(const std::vector<cvf::Vec4d>& triangles);
std::vector<CellIndexAndResult> cellOverlapVolumesAndResults(const cvf::Vec2d& globalPos2d,
std::vector<CellIndexAndResult> cellOverlapVolumesAndResults(const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues) const;
std::vector<CellIndexAndResult> cellRayIntersectionAndResults(const cvf::Vec2d& globalPos2d,
std::vector<CellIndexAndResult> cellRayIntersectionAndResults(const cvf::Vec2d& globalPos2d,
const std::vector<double>& weightingResultValues) const;
bool isMeanResult() const;
@ -219,6 +224,7 @@ protected:
std::vector<double> m_contourLevelCumulativeAreas;
std::vector<cvf::Vec4d> m_trianglesWithVertexValues;
int m_currentResultTimestep;
std::vector<bool> m_mapCellVisibility;
double m_minResultAllTimeSteps;
double m_maxResultAllTimeSteps;

View File

@ -29,8 +29,6 @@
#include "RigEclipseCaseData.h"
#include "RigHexIntersectionTools.h"
#include "RigMainGrid.h"
#include "RigResultAccessor.h"
#include "RigResultAccessorFactory.h"
#include "RimCellRangeFilterCollection.h"
#include "RimEclipseContourMapView.h"
@ -75,8 +73,6 @@ RimEclipseContourMapProjection::RimEclipseContourMapProjection()
setName("Map Projection");
nameField()->uiCapability()->setUiReadOnly(true);
m_resultAccessor = new RigHugeValResultAccessor;
}
//--------------------------------------------------------------------------------------------------
@ -186,7 +182,7 @@ std::vector<double> RimEclipseContourMapProjection::generateResults(int timeStep
if (!cellColors->isTernarySaturationSelected())
{
RigCaseCellResultsData* resultData = eclipseCase->results(RiaDefines::MATRIX_MODEL);
std::vector<double> gridResultValues;
if (isColumnResult())
{
m_currentResultName = "";
@ -201,26 +197,23 @@ std::vector<double> RimEclipseContourMapProjection::generateResults(int timeStep
{
resultData->findOrLoadScalarResultForTimeStep(RiaDefines::DYNAMIC_NATIVE, "SGAS", timeStep);
}
gridResultValues = calculateColumnResult(m_resultAggregation());
}
else
{
m_currentResultName = cellColors->resultVariable();
m_resultAccessor =
RigResultAccessorFactory::createFromResultDefinition(eclipseCase->eclipseCaseData(), 0, timeStep, cellColors);
if (m_resultAccessor.isNull())
{
m_resultAccessor = new RigHugeValResultAccessor;
}
size_t resultIndex = resultData->findScalarResultIndex(cellColors->resultType(), cellColors->resultVariable());
gridResultValues = resultData->cellScalarResults(resultIndex)[timeStep];
}
int everyNCells = temporaryResult ? 10 : 1;
#pragma omp parallel for
for (int index = 0; index < static_cast<int>(nCells); ++index)
if (!gridResultValues.empty())
{
cvf::Vec2ui ij = ijFromCellIndex(index);
aggregatedResults[index] = calculateValueInMapCell(ij.x(), ij.y());
#pragma omp parallel for
for (int index = 0; index < static_cast<int>(nCells); ++index)
{
cvf::Vec2ui ij = ijFromCellIndex(index);
aggregatedResults[index] = calculateValueInMapCell(ij.x(), ij.y(), gridResultValues);
}
}
}
}
@ -254,99 +247,51 @@ void RimEclipseContourMapProjection::clearResultVariable()
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimEclipseContourMapProjection::calculateValueInMapCell(uint i, uint j) const
{
if (isColumnResult())
{
const std::vector<std::pair<size_t, double>>& matchingCells = cellsAtIJ(i, j);
if (!matchingCells.empty())
{
switch (m_resultAggregation())
{
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 = calculateColumnResult(m_resultAggregation(), cellIdx);
sum += cellValue * cellIdxAndWeight.second;
}
return sum;
}
default:
CVF_TIGHT_ASSERT(false);
}
}
}
else
{
if (!view()->cellResult()->isFlowDiagOrInjectionFlooding() &&
view()->cellResult()->scalarResultIndex() == cvf::UNDEFINED_SIZE_T)
{
return 0.0; // Special case of NONE-result. Show 0 all over to ensure we see something.
}
return RimContourMapProjection::calculateValueInMapCell(i, j);
}
return std::numeric_limits<double>::infinity();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimEclipseContourMapProjection::calculateColumnResult(ResultAggregation resultAggregation, size_t cellGlobalIdx) const
std::vector<double> RimEclipseContourMapProjection::calculateColumnResult(ResultAggregation resultAggregation) 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)
if (poroResultIndex == cvf::UNDEFINED_SIZE_T || ntgResultIndex == cvf::UNDEFINED_SIZE_T || dzResultIndex == cvf::UNDEFINED_SIZE_T)
{
return std::numeric_limits<double>::infinity();
return std::vector<double>();
}
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);
CVF_ASSERT(poroResults.size() == ntgResults.size() && ntgResults.size() == dzResults.size());
int timeStep = view()->currentTimeStep();
double resultValue = 0.0;
std::vector<double> resultValues(poroResults.size(), 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())
for (size_t cellResultIdx = 0; cellResultIdx < resultValues.size(); ++cellResultIdx)
{
resultValue = soilResults.at(cellResultIdx);
resultValues[cellResultIdx] = soilResults[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())
for (size_t cellResultIdx = 0; cellResultIdx < resultValues.size(); ++cellResultIdx)
{
resultValue += sgasResults.at(cellResultIdx);
resultValues[cellResultIdx] += sgasResults[cellResultIdx];
}
}
return resultValue * poro * ntg * dz;
for (size_t cellResultIdx = 0; cellResultIdx < resultValues.size(); ++cellResultIdx)
{
resultValues[cellResultIdx] *= poroResults[cellResultIdx] * ntgResults[cellResultIdx] * dzResults[cellResultIdx];
}
return resultValues;
}
//--------------------------------------------------------------------------------------------------
@ -368,7 +313,7 @@ void RimEclipseContourMapProjection::updateGridInformation()
cvf::Vec3d maxPoint = m_expandedBoundingBox.max();
maxPoint.x() = minPoint.x() + m_mapSize.x() * m_sampleSpacing;
maxPoint.y() = minPoint.y() + m_mapSize.y() * m_sampleSpacing;
m_expandedBoundingBox = cvf::BoundingBox(minPoint, maxPoint);
m_expandedBoundingBox = cvf::BoundingBox(minPoint, maxPoint);
}
//--------------------------------------------------------------------------------------------------
@ -425,17 +370,23 @@ std::vector<size_t> RimEclipseContourMapProjection::findIntersectingCells(const
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimEclipseContourMapProjection::calculateOverlapVolume(size_t globalCellIdx, const cvf::BoundingBox& bbox, size_t* cellKLayerOut) const
size_t RimEclipseContourMapProjection::kLayer(size_t globalCellIdx) const
{
CVF_ASSERT(cellKLayerOut != nullptr);
std::array<cvf::Vec3d, 8> hexCorners;
RigCell cell = m_mainGrid->globalCellArray()[globalCellIdx];
const RigCell& cell = m_mainGrid->globalCellArray()[globalCellIdx];
size_t mainGridCellIdx = cell.mainGridCellIndex();
size_t i, j, k;
m_mainGrid->ijkFromCellIndex(mainGridCellIdx, &i, &j, &k);
return k;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimEclipseContourMapProjection::calculateOverlapVolume(size_t globalCellIdx, const cvf::BoundingBox& bbox) const
{
std::array<cvf::Vec3d, 8> hexCorners;
const RigCell& cell = m_mainGrid->globalCellArray()[globalCellIdx];
size_t localCellIdx = cell.gridLocalCellIndex();
RigGridBase* localGrid = cell.hostGrid();
@ -447,9 +398,6 @@ double RimEclipseContourMapProjection::calculateOverlapVolume(size_t globalCellI
RigCellGeometryTools::estimateHexOverlapWithBoundingBox(hexCorners, bbox, &overlapBBox);
double overlapVolume = RigCellGeometryTools::calculateCellVolume(overlapCorners);
*cellKLayerOut = k;
return overlapVolume;
}
@ -458,19 +406,12 @@ double RimEclipseContourMapProjection::calculateOverlapVolume(size_t globalCellI
//--------------------------------------------------------------------------------------------------
double RimEclipseContourMapProjection::calculateRayLengthInCell(size_t globalCellIdx,
const cvf::Vec3d& highestPoint,
const cvf::Vec3d& lowestPoint,
size_t* cellKLayerOut) const
const cvf::Vec3d& lowestPoint) const
{
CVF_ASSERT(cellKLayerOut);
std::array<cvf::Vec3d, 8> hexCorners;
RigCell cell = m_mainGrid->globalCellArray()[globalCellIdx];
size_t mainGridCellIdx = cell.mainGridCellIndex();
size_t i, j, k;
m_mainGrid->ijkFromCellIndex(mainGridCellIdx, &i, &j, &k);
size_t localCellIdx = cell.gridLocalCellIndex();
RigGridBase* localGrid = cell.hostGrid();
@ -480,7 +421,6 @@ double RimEclipseContourMapProjection::calculateRayLengthInCell(size_t
if (RigHexIntersectionTools::lineHexCellIntersection(highestPoint, lowestPoint, hexCorners.data(), 0, &intersections))
{
double lengthInCell = (intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint).length();
*cellKLayerOut = k;
return lengthInCell;
}
return 0.0;
@ -489,12 +429,10 @@ double RimEclipseContourMapProjection::calculateRayLengthInCell(size_t
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimEclipseContourMapProjection::getParameterWeightForCell(size_t globalCellIdx, const std::vector<double>& cellWeights) const
double RimEclipseContourMapProjection::getParameterWeightForCell(size_t cellResultIdx, const std::vector<double>& cellWeights) const
{
if (cellWeights.empty()) return 1.0;
const RigActiveCellInfo* activeCellInfo = eclipseCase()->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL);
size_t cellResultIdx = activeCellInfo->cellResultIndex(globalCellIdx);
double result = std::max(cellWeights[cellResultIdx], 0.0);
if (result < 1.0e-6)
{
@ -506,9 +444,10 @@ double RimEclipseContourMapProjection::getParameterWeightForCell(size_t globalCe
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimEclipseContourMapProjection::gridCellValue(size_t globalCellIdx) const
size_t RimEclipseContourMapProjection::gridResultIndex(size_t globalCellIdx) const
{
return m_resultAccessor->cellScalarGlobIdx(globalCellIdx);
const RigActiveCellInfo* activeCellInfo = eclipseCase()->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL);
return activeCellInfo->cellResultIndex(globalCellIdx);
}
//--------------------------------------------------------------------------------------------------

View File

@ -68,14 +68,14 @@ protected:
void clearResultVariable() override;
RimGridView* baseView() const override;
std::vector<size_t> findIntersectingCells(const cvf::BoundingBox& bbox) const override;
double calculateOverlapVolume(size_t globalCellIdx, const cvf::BoundingBox& bbox, size_t* cellKLayerOut) const override;
double calculateRayLengthInCell(size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint, size_t* cellKLayerOut) const override;
double getParameterWeightForCell(size_t globalCellIdx, const std::vector<double>& parameterWeights) const override;
double gridCellValue(size_t globalCellIdx) const override;
size_t kLayer(size_t globalCellIdx) 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
double calculateValueInMapCell(uint i, uint j) const;
double calculateColumnResult(ResultAggregation resultAggregation, size_t cellGlobalIdx) const;
std::vector<double> calculateColumnResult(ResultAggregation resultAggregation) const;
RimEclipseResultCase* eclipseCase() const;
RimEclipseContourMapView* view() const;
@ -91,8 +91,6 @@ protected:
caf::PdmField<bool> m_weightByParameter;
caf::PdmChildField<RimEclipseResultDefinition*> m_weightingResult;
cvf::ref<RigResultAccessor> m_resultAccessor;
cvf::ref<RigMainGrid> m_mainGrid;
QString m_currentResultName;
};

View File

@ -36,6 +36,8 @@
#include "RivFemElmVisibilityCalculator.h"
#include "cafPdmUiDoubleSliderEditor.h"
#include "cvfArray.h"
#include "cvfCellRange.h"
#include "cvfGeometryTools.h"
@ -58,7 +60,8 @@ RimGeoMechContourMapProjection::RimGeoMechContourMapProjection()
{
CAF_PDM_InitObject("RimContourMapProjection", ":/2DMapProjection16x16.png", "", "");
CAF_PDM_InitField(&m_limitToPorePressureRegions, "LimitToPorRegion", true, "Limit to Pore Pressure regions", "", "", "");
CAF_PDM_InitField(&m_includePaddingAroundPorePressureRegion, "PaddingAroundPorRegion", true, "Include Padding around Pore Pressure regions", "", "", "");
CAF_PDM_InitField(&m_paddingAroundPorePressureRegion, "PaddingAroundPorRegion", 0.0, "Padding around PP regions", "", "", "");
m_paddingAroundPorePressureRegion.uiCapability()->setUiEditorTypeName(caf::PdmUiDoubleSliderEditor::uiEditorTypeName());
setName("Map Projection");
nameField()->uiCapability()->setUiReadOnly(true);
}
@ -95,7 +98,7 @@ void RimGeoMechContourMapProjection::updateLegend()
double minVal = minValue(m_aggregatedResults);
double maxVal = maxValue(m_aggregatedResults);
std::pair<double, double> minmaxValAllTimeSteps = minmaxValuesAllTimeSteps(2);
std::pair<double, double> minmaxValAllTimeSteps = minmaxValuesAllTimeSteps();
legendConfig()->setAutomaticRanges(minmaxValAllTimeSteps.first, minmaxValAllTimeSteps.second, minVal, maxVal);
@ -122,32 +125,157 @@ cvf::ref<cvf::UByteArray> RimGeoMechContourMapProjection::getCellVisibility() co
RivFemElmVisibilityCalculator::computeRangeVisibility(cellGridIdxVisibility.p(), m_femPart.p(), cellRangeFilter);
RivFemElmVisibilityCalculator::computePropertyVisibility(cellGridIdxVisibility.p(), m_femPart.p(), view()->currentTimeStep(), cellGridIdxVisibility.p(), view()->geoMechPropertyFilterCollection());
if (m_limitToPorePressureRegions)
{
ensureOnlyValidPorBarVisible(cellGridIdxVisibility.p(), view()->currentTimeStep());
}
return cellGridIdxVisibility;
}
//--------------------------------------------------------------------------------------------------
/// Meijster, Roerdink, Hesselink
/// A GENERAL ALGORITHM FOR COMPUTING DISTANCE TRANSFORMS IN LINEAR TIME
/// http://fab.cba.mit.edu/classes/S62.12/docs/Meijster_distance.pdf
//--------------------------------------------------------------------------------------------------
void distanceTransform(std::vector<std::vector<uint>>& image)
{
if (image.empty())
{
return;
}
if (image.front().empty())
{
return;
}
const int64_t M = (int64_t) image.size();
const int64_t N = (int64_t) image.front().size();
unsigned int uinf = M + N;
// First phase
std::vector<std::vector<uint>> g(M);
#pragma omp parallel for
for (int64_t x = 0; x < M; ++x)
{
g[x].resize(N, uinf);
if (image[x][0])
{
g[x][0] = 0;
}
for (int64_t y = 1; y < N - 1; ++y)
{
if (image[x][y])
{
g[x][y] = 0;
}
else
{
g[x][y] = 1 + g[x][y - 1];
}
}
for (int64_t y = N - 2; y > 0; --y)
{
if (g[x][y + 1] < g[x][y])
{
g[x][y] = 1 + g[x][y + 1];
}
}
}
auto f = [](int64_t x, int64_t i, const std::vector<std::vector<uint>>& g, int64_t y)
{
return (x - i) * (x - i) + g[i][y] * g[i][y];
};
auto sep = [](int64_t i, int64_t u, const std::vector<std::vector<uint>>& g, int64_t y)
{
if (i == u) return (int64_t) 0;
int64_t numerator = u * u - i * i + g[u][y] * g[u][y] - g[i][y] * g[i][y];
int64_t divisor = 2 * (u - i);
return numerator / divisor;
};
// Second phase
#pragma omp parallel for
for (int64_t y = 0; y < N; ++y)
{
int64_t q = 0;
std::vector<uint> s(std::max(N, M), 0u);
std::vector<uint> t(std::max(N, M), 0u);
for (int64_t u = 1; u < M - 1; ++u)
{
while (q >= 0 && f(t[q], s[q], g, y) > f(t[q], u, g, y))
{
q--;
}
if (q < 0)
{
q = 0;
s[0] = u;
}
else
{
int64_t w = 1 + sep((int64_t) s[q], u, g, y);
if (w < M)
{
q++;
s[q] = u;
t[q] = w;
}
}
}
for (int64_t u = M - 1; u > 0; --u)
{
image[u][y] = f(u, s[q], g, y);
if (u == t[q])
{
q = q - 1;
}
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimGeoMechContourMapProjection::ensureOnlyValidPorBarVisible(cvf::UByteArray* visibility, int timeStep) const
cvf::BoundingBox RimGeoMechContourMapProjection::calculateExpandedPorBarBBox(int timeStep) const
{
RigFemResultAddress porBarAddr(RigFemResultPosEnum::RIG_ELEMENT_NODAL, "POR-Bar", view()->cellResult()->resultComponentName().toStdString());
RigFemResultAddress porBarAddr(
RigFemResultPosEnum::RIG_ELEMENT_NODAL, "POR-Bar", view()->cellResult()->resultComponentName().toStdString());
RigGeoMechCaseData* caseData = geoMechCase()->geoMechData();
RigFemPartResultsCollection* resultCollection = caseData->femPartResults();
const std::vector<float>& resultValues = resultCollection->resultValues(porBarAddr, 0, timeStep);
for (int i = 0; i < static_cast<int>(visibility->size()); ++i)
cvf::BoundingBox boundingBox;
for (int i = 0; i < m_femPart->elementCount(); ++i)
{
size_t resValueIdx = m_femPart->elementNodeResultIdx((int) i, 0);
size_t resValueIdx = m_femPart->elementNodeResultIdx((int)i, 0);
CVF_ASSERT(resValueIdx < resultValues.size());
double scalarValue = resultValues[resValueIdx];
(*visibility)[i] &= scalarValue != std::numeric_limits<double>::infinity();
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;
boxMin.y() -= boxExtent.y() * 0.5;
boxMin.z() -= geoMechCase()->characteristicCellSize();
boxMax.x() += boxExtent.x() * 0.5;
boxMax.y() += boxExtent.y() * 0.5;
boxMax.z() += geoMechCase()->characteristicCellSize();
return cvf::BoundingBox(boxMin, boxMax);
}
//--------------------------------------------------------------------------------------------------
@ -159,12 +287,19 @@ void RimGeoMechContourMapProjection::updateGridInformation()
m_femPart = geoMechCase->geoMechData()->femParts()->part(0);
m_femPartGrid = m_femPart->getOrCreateStructGrid();
m_sampleSpacing = m_relativeSampleSpacing * geoMechCase->characteristicCellSize();
m_gridBoundingBox = geoMechCase->activeCellsBoundingBox();
m_femPart->ensureIntersectionSearchTreeIsBuilt();
if (m_limitToPorePressureRegions)
{
m_gridBoundingBox = calculateExpandedPorBarBBox(view()->currentTimeStep());
}
else
{
m_gridBoundingBox = geoMechCase->activeCellsBoundingBox();
}
cvf::Vec3d minExpandedPoint = m_gridBoundingBox.min() - cvf::Vec3d(gridEdgeOffset(), gridEdgeOffset(), 0.0);
cvf::Vec3d maxExpandedPoint = m_gridBoundingBox.max() + cvf::Vec3d(gridEdgeOffset(), gridEdgeOffset(), 0.0);
m_expandedBoundingBox = cvf::BoundingBox(minExpandedPoint, maxExpandedPoint);
m_expandedBoundingBox = cvf::BoundingBox(minExpandedPoint, maxExpandedPoint);
m_mapSize = calculateMapSize();
@ -176,6 +311,62 @@ void RimGeoMechContourMapProjection::updateGridInformation()
m_expandedBoundingBox = cvf::BoundingBox(minPoint, maxPoint);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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 porBarAddr(RigFemResultPosEnum::RIG_ELEMENT_NODAL, "POR-Bar", "");
std::vector<double> cellResults = generateResultsFromAddress(porBarAddr, mapCellVisibility, view()->currentTimeStep());
mapCellVisibility.resize(numberOfCells(), true);
CVF_ASSERT(mapCellVisibility.size() == cellResults.size());
if (m_limitToPorePressureRegions)
{
cvf::BoundingBox validPorBoundingBox;
for (size_t cellIndex = 0; cellIndex < cellResults.size(); ++cellIndex)
{
cvf::Vec2ui ij = ijFromCellIndex(cellIndex);
validPorBoundingBox.add(cvf::Vec3d(cellCenterPosition(ij.x(), ij.y()), 0.0));
if (cellResults[cellIndex] != std::numeric_limits<double>::infinity())
{
distanceImage[ij.x()][ij.y()] = 1u;
}
else
{
mapCellVisibility[cellIndex] = false;
}
}
if (m_paddingAroundPorePressureRegion > 0.0)
{
distanceTransform(distanceImage);
cvf::Vec3d porExtent = validPorBoundingBox.extent();
double radius = std::max(porExtent.x(), porExtent.y()) * 0.25;
double expansion = m_paddingAroundPorePressureRegion * radius;
size_t cellPadding = std::ceil(expansion / m_sampleSpacing);
for (size_t cellIndex = 0; cellIndex < cellResults.size(); ++cellIndex)
{
if (!mapCellVisibility[cellIndex])
{
cvf::Vec2ui ij = ijFromCellIndex(cellIndex);
if (distanceImage[ij.x()][ij.y()] < cellPadding * cellPadding)
{
mapCellVisibility[cellIndex] = true;
}
}
}
}
}
return mapCellVisibility;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@ -188,42 +379,54 @@ std::vector<double> RimGeoMechContourMapProjection::retrieveParameterWeights()
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RimGeoMechContourMapProjection::generateResults(int timeStep)
{
RimGeoMechCellColors* cellColors = view()->cellResult();
RigFemResultAddress resultAddress = cellColors->resultAddress();
std::vector<double> aggregatedResults = generateResultsFromAddress(resultAddress, m_mapCellVisibility, timeStep);
return aggregatedResults;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RimGeoMechContourMapProjection::generateResultsFromAddress(RigFemResultAddress resultAddress, const std::vector<bool>& mapCellVisibility, int timeStep)
{
RigGeoMechCaseData* caseData = geoMechCase()->geoMechData();
RigFemPartResultsCollection* resultCollection = caseData->femPartResults();
RimGeoMechCellColors* cellColors = view()->cellResult();
RigFemResultAddress resAddr = cellColors->resultAddress();
size_t nCells = numberOfCells();
std::vector<double> aggregatedResults = std::vector<double>(nCells, std::numeric_limits<double>::infinity());
if (!resAddr.isValid())
if (!resultAddress.isValid())
return aggregatedResults;
if (resAddr.fieldName == "PP")
if (resultAddress.fieldName == "PP")
{
resAddr.fieldName = "POR-Bar"; // More likely to be in memory than POR
resultAddress.fieldName = "POR-Bar"; // More likely to be in memory than POR
}
if (resultAddress.fieldName == "POR-Bar")
{
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.
}
if (resAddr.fieldName == "POR-Bar") resAddr.resultPosType = RIG_ELEMENT_NODAL;
m_resultValues = resultCollection->resultValues(resAddr, 0, timeStep);
int everyNCells = temporaryResult ? 10 : 1;
std::vector<float> resultValuesF = resultCollection->resultValues(resultAddress, 0, timeStep);
std::vector<double> resultValues = gridCellValues(resultAddress, resultValuesF);
#pragma omp parallel for
for (int index = 0; index < static_cast<int>(nCells); ++index)
{
cvf::Vec2ui ij = ijFromCellIndex(index);
aggregatedResults[index] = calculateValueInMapCell(ij.x(), ij.y());
if (mapCellVisibility.empty() || mapCellVisibility[index])
{
cvf::Vec2ui ij = ijFromCellIndex(index);
aggregatedResults[index] = calculateValueInMapCell(ij.x(), ij.y(), resultValues);
}
}
m_currentResultAddr = resAddr;
if (temporaryResult)
{
resultCollection->deleteResultFrame(resAddr, 0, timeStep);
}
return aggregatedResults;
}
@ -273,40 +476,28 @@ std::vector<size_t> RimGeoMechContourMapProjection::findIntersectingCells(const
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimGeoMechContourMapProjection::calculateOverlapVolume(size_t globalCellIdx,
const cvf::BoundingBox& bbox,
size_t* cellKLayerOut) const
size_t RimGeoMechContourMapProjection::kLayer(size_t globalCellIdx) const
{
CVF_ASSERT(cellKLayerOut != nullptr);
size_t i, j, k;
m_femPartGrid->ijkFromCellIndex(globalCellIdx, &i, &j, &k);
return k;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimGeoMechContourMapProjection::calculateOverlapVolume(size_t globalCellIdx,
const cvf::BoundingBox& bbox) 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]]);
m_femPartGrid->cellCornerVertices(globalCellIdx, hexCorners.data());
cvf::BoundingBox overlapBBox;
std::array<cvf::Vec3d, 8> overlapCorners =
RigCellGeometryTools::estimateHexOverlapWithBoundingBox(hexCorners, bbox, &overlapBBox);
double overlapVolume = RigCellGeometryTools::calculateCellVolume(overlapCorners);
if (overlapVolume > 0.0)
{
size_t i, j, k;
m_femPartGrid->ijkFromCellIndex(globalCellIdx, &i, &j, &k);
*cellKLayerOut = k;
return overlapVolume;
}
return 0.0;
return overlapVolume;
}
//--------------------------------------------------------------------------------------------------
@ -314,11 +505,8 @@ double RimGeoMechContourMapProjection::calculateOverlapVolume(size_t
//--------------------------------------------------------------------------------------------------
double RimGeoMechContourMapProjection::calculateRayLengthInCell(size_t globalCellIdx,
const cvf::Vec3d& highestPoint,
const cvf::Vec3d& lowestPoint,
size_t* cellKLayerOut) const
const cvf::Vec3d& lowestPoint) const
{
CVF_ASSERT(cellKLayerOut != nullptr);
std::array<cvf::Vec3d, 8> hexCorners;
const std::vector<cvf::Vec3f>& nodeCoords = m_femPart->nodes().coordinates;
@ -338,10 +526,6 @@ double RimGeoMechContourMapProjection::calculateRayLengthInCell(size_t
if (RigHexIntersectionTools::lineHexCellIntersection(highestPoint, lowestPoint, hexCorners.data(), 0, &intersections))
{
double lengthInCell = (intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint).length();
size_t i, j, k;
m_femPartGrid->ijkFromCellIndex(globalCellIdx, &i, &j, &k);
*cellKLayerOut = k;
return lengthInCell;
}
return 0.0;
@ -361,53 +545,42 @@ double RimGeoMechContourMapProjection::getParameterWeightForCell(size_t
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimGeoMechContourMapProjection::gridCellValue(size_t globalCellIdx) const
std::vector<double> RimGeoMechContourMapProjection::gridCellValues(RigFemResultAddress resAddr, std::vector<float>& resultValues) const
{
RimGeoMechCellColors* cellColors = view()->cellResult();
RigFemResultAddress resAddr = cellColors->resultAddress();
if (resAddr.fieldName == "PP")
std::vector<double> gridCellValues(m_femPart->elementCount(), std::numeric_limits<double>::infinity());
for (size_t globalCellIdx = 0; globalCellIdx < m_femPart->elementCount(); ++globalCellIdx)
{
resAddr.fieldName = "POR-Bar"; // More likely to be in memory than POR
}
if (resAddr.fieldName == "POR-Bar") resAddr.resultPosType = RIG_ELEMENT_NODAL;
RigElementType elmType = m_femPart->elementType(globalCellIdx);
if (!(elmType == HEX8 || elmType == HEX8P)) continue;
RigFemResultPosEnum resultPosType = resAddr.resultPosType;
RigElementType elmType = m_femPart->elementType(globalCellIdx);
if (!(elmType == HEX8 || elmType == HEX8P)) return 0.0;
if (resultPosType == RIG_FORMATION_NAMES)
{
resultPosType = RIG_ELEMENT_NODAL; // formation indices are stored per element node result.
}
if (resultPosType == RIG_ELEMENT)
{
return m_resultValues[globalCellIdx];
}
else if (resultPosType == RIG_ELEMENT_NODAL)
{
RiaWeightedMeanCalculator<float> cellAverage;
for (int i = 0; i < 8; ++i)
if (resAddr.resultPosType == RIG_ELEMENT)
{
size_t gridResultValueIdx = m_femPart->resultValueIdxFromResultPosType(resultPosType, static_cast<int>(globalCellIdx), i);
cellAverage.addValueAndWeight(m_resultValues[gridResultValueIdx], 1.0);
gridCellValues[globalCellIdx] = static_cast<double>(resultValues[globalCellIdx]);
}
return cellAverage.weightedMean();
}
else
{
RiaWeightedMeanCalculator<float> cellAverage;
const int* elmNodeIndices = m_femPart->connectivities(globalCellIdx);
for (int i = 0; i < 8; ++i)
else if (resAddr.resultPosType == RIG_ELEMENT_NODAL)
{
cellAverage.addValueAndWeight(m_resultValues[elmNodeIndices[i]], 1.0);
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 cellAverage.weightedMean();
}
return gridCellValues;
}
//--------------------------------------------------------------------------------------------------
@ -438,7 +611,7 @@ void RimGeoMechContourMapProjection::fieldChangedByUi(const caf::PdmFieldHandle*
const QVariant& newValue)
{
RimContourMapProjection::fieldChangedByUi(changedField, oldValue, newValue);
if (changedField == &m_limitToPorePressureRegions)
if (changedField == &m_limitToPorePressureRegions || changedField == &m_paddingAroundPorePressureRegion)
{
clearGridMapping();
}
@ -478,7 +651,27 @@ void RimGeoMechContourMapProjection::defineUiOrdering(QString uiConfigName, caf:
RimContourMapProjection::defineUiOrdering(uiConfigName, uiOrdering);
caf::PdmUiGroup* group = uiOrdering.addNewGroup("Grid Boundaries");
group->add(&m_limitToPorePressureRegions);
group->add(&m_paddingAroundPorePressureRegion);
m_paddingAroundPorePressureRegion.uiCapability()->setUiReadOnly(!m_limitToPorePressureRegions());
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimGeoMechContourMapProjection::defineEditorAttribute(const caf::PdmFieldHandle* field,
QString uiConfigName,
caf::PdmUiEditorAttribute* attribute)
{
RimContourMapProjection::defineEditorAttribute(field, uiConfigName, attribute);
if (field == &m_paddingAroundPorePressureRegion)
{
caf::PdmUiDoubleSliderEditorAttribute* myAttr = dynamic_cast<caf::PdmUiDoubleSliderEditorAttribute*>(attribute);
if (myAttr)
{
myAttr->m_minimum = 0.0;
myAttr->m_maximum = 1.0;
myAttr->m_sliderTickCount = 4;
}
}
}

View File

@ -59,22 +59,23 @@ public:
protected:
typedef RimContourMapProjection::CellIndexAndResult CellIndexAndResult;
// GeoMech implementation specific data generation methods
cvf::ref<cvf::UByteArray> getCellVisibility() const override;
void ensureOnlyValidPorBarVisible(cvf::UByteArray* visibility, int timeStep) const;
cvf::BoundingBox calculateExpandedPorBarBBox(int timeStep) const;
void updateGridInformation() override;
std::vector<bool> getMapCellVisibility() override;
std::vector<double> retrieveParameterWeights() override;
std::vector<double> generateResults(int timeStep) override;
std::vector<double> generateResultsFromAddress(RigFemResultAddress resultAddress, const std::vector<bool>& mapCellVisibility, int timeStep);
bool resultVariableChanged() const override;
void clearResultVariable() override;
RimGridView* baseView() const override;
std::vector<size_t> findIntersectingCells(const cvf::BoundingBox& bbox) const override;
double calculateOverlapVolume(size_t globalCellIdx, const cvf::BoundingBox& bbox, size_t* cellKLayerOut) const override;
double calculateRayLengthInCell(size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint, size_t* cellKLayerOut) const override;
size_t kLayer(size_t globalCellIdx) 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;
// GeoMech implementation specific data generation methods
double gridCellValue(size_t globalCellIdx) const override;
std::vector<double> gridCellValues(RigFemResultAddress resAddr, std::vector<float>& resultValues) const;
RimGeoMechCase* geoMechCase() const;
RimGeoMechContourMapView* view() const;
@ -84,15 +85,16 @@ protected:
QList<caf::PdmOptionItemInfo> calculateValueOptions(const caf::PdmFieldHandle* fieldNeedingOptions,
bool* useOptionsOnly) override;
void defineUiOrdering(QString uiConfigName, caf::PdmUiOrdering& uiOrdering) override;
void defineEditorAttribute(const caf::PdmFieldHandle* field,
QString uiConfigName,
caf::PdmUiEditorAttribute* attribute) override;
protected:
caf::PdmField<bool> m_limitToPorePressureRegions;
caf::PdmField<bool> m_includePaddingAroundPorePressureRegion;
caf::PdmField<double> m_paddingAroundPorePressureRegion;
cvf::ref<RigFemPart> m_femPart;
cvf::cref<RigFemPartGrid> m_femPartGrid;
RigFemResultAddress m_currentResultAddr;
std::vector<float> m_resultValues;
};