OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-01-26 16:26:48 -06:00
Code Issues Packages Projects Releases Wiki Activity
796c6a0e0a
ResInsight/ApplicationCode/ProjectDataModel/RimGeoMechContourMapProjection.cpp
Gaute Lindkvist 796c6a0e0a #3987 Fix picking on Geo Mech contour maps
2019-01-23 09:23:22 +01:00

579 lines
25 KiB
C++
Raw Blame History

/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2018- Equinor ASA
//
// ResInsight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "RimGeoMechContourMapProjection.h"
#include "RiaImageTools.h"
#include "RiaWeightedGeometricMeanCalculator.h"
#include "RiaWeightedHarmonicMeanCalculator.h"
#include "RiaWeightedMeanCalculator.h"
#include "RigCellGeometryTools.h"
#include "RigGeoMechCaseData.h"
#include "RigFemPart.h"
#include "RigFemPartGrid.h"
#include "RigFemPartCollection.h"
#include "RigFemPartResultsCollection.h"
#include "RigHexIntersectionTools.h"
#include "RimCellRangeFilterCollection.h"
#include "RimGeoMechContourMapView.h"
#include "RimGeoMechCellColors.h"
#include "RimGeoMechPropertyFilterCollection.h"
#include "RivFemElmVisibilityCalculator.h"
#include "cafPdmUiDoubleSliderEditor.h"
#include "cvfArray.h"
#include "cvfCellRange.h"
#include "cvfGeometryTools.h"
#include "cvfGeometryUtils.h"
#include "cvfScalarMapper.h"
#include "cvfStructGridGeometryGenerator.h"
#include "cvfVector3.h"
#include <algorithm>
#include <array>
#include <omp.h>
#include <QDebug>
CAF_PDM_SOURCE_INIT(RimGeoMechContourMapProjection, "RimGeoMechContourMapProjection");
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimGeoMechContourMapProjection::RimGeoMechContourMapProjection()
{
CAF_PDM_InitObject("RimContourMapProjection", ":/2DMapProjection16x16.png", "", "");
CAF_PDM_InitField(&m_limitToPorePressureRegions, "LimitToPorRegion", true, "Limit to 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);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimGeoMechContourMapProjection::~RimGeoMechContourMapProjection() {}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QString RimGeoMechContourMapProjection::resultDescriptionText() const
{
QString resultText = QString("%1, %2").arg(resultAggregationText()).arg(view()->cellResult()->resultFieldUiName());
return resultText;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimRegularLegendConfig* RimGeoMechContourMapProjection::legendConfig() const
{
return view()->cellResult()->legendConfig();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimGeoMechContourMapProjection::updateLegend()
{
RimGeoMechCellColors* cellColors = view()->cellResult();
double minVal = minValue(m_aggregatedResults);
double maxVal = maxValue(m_aggregatedResults);
std::pair<double, double> minmaxValAllTimeSteps = minmaxValuesAllTimeSteps();
legendConfig()->setAutomaticRanges(minmaxValAllTimeSteps.first, minmaxValAllTimeSteps.second, minVal, maxVal);
QString projectionLegendText = QString("Map Projection\n%1").arg(m_resultAggregation().uiText());
projectionLegendText += QString("\nResult: %1").arg(cellColors->resultFieldUiName());
if (!cellColors->resultComponentUiName().isEmpty())
{
projectionLegendText += QString(", %1").arg(cellColors->resultComponentUiName());
}
legendConfig()->setTitle(projectionLegendText);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::ref<cvf::UByteArray> RimGeoMechContourMapProjection::getCellVisibility() const
{
cvf::ref<cvf::UByteArray> cellGridIdxVisibility = new cvf::UByteArray(m_femPart->elementCount());
RivFemElmVisibilityCalculator::computeAllVisible(cellGridIdxVisibility.p(), m_femPart.p());
cvf::CellRangeFilter cellRangeFilter;
view()->rangeFilterCollection()->compoundCellRangeFilter(&cellRangeFilter, 0);
RivFemElmVisibilityCalculator::computeRangeVisibility(cellGridIdxVisibility.p(), m_femPart.p(), cellRangeFilter);
RivFemElmVisibilityCalculator::computePropertyVisibility(cellGridIdxVisibility.p(), m_femPart.p(), view()->currentTimeStep(), cellGridIdxVisibility.p(), view()->geoMechPropertyFilterCollection());
return cellGridIdxVisibility;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::BoundingBox RimGeoMechContourMapProjection::calculateExpandedPorBarBBox(int timeStep) const
{
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);
cvf::BoundingBox boundingBox;
for (int i = 0; i < m_femPart->elementCount(); ++i)
{
size_t resValueIdx = m_femPart->elementNodeResultIdx((int)i, 0);
CVF_ASSERT(resValueIdx < resultValues.size());
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);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimGeoMechContourMapProjection::updateGridInformation()
{
RimGeoMechCase* geoMechCase = this->geoMechCase();
m_femPart = geoMechCase->geoMechData()->femParts()->part(0);
m_femPartGrid = m_femPart->getOrCreateStructGrid();
m_sampleSpacing = m_relativeSampleSpacing * geoMechCase->characteristicCellSize();
m_femPart->ensureIntersectionSearchTreeIsBuilt();
m_gridBoundingBox = geoMechCase->activeCellsBoundingBox();
if (m_limitToPorePressureRegions)
{
m_gridBoundingBox = calculateExpandedPorBarBBox(view()->currentTimeStep());
}
else
{
m_expandedBoundingBox = m_gridBoundingBox;
}
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_mapSize = calculateMapSize();
// Re-jig max point to be an exact multiple of cell size
cvf::Vec3d minPoint = m_expandedBoundingBox.min();
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);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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 resAddr = view()->cellResult()->resultAddress();
if (m_limitToPorePressureRegions)
{
resAddr = RigFemResultAddress(RigFemResultPosEnum::RIG_ELEMENT_NODAL, "POR-Bar", "");
}
std::vector<double> cellResults = generateResultsFromAddress(resAddr, mapCellVisibility, view()->currentTimeStep());
mapCellVisibility.resize(numberOfCells(), true);
CVF_ASSERT(mapCellVisibility.size() == cellResults.size());
{
cvf::BoundingBox validResBoundingBox;
for (size_t cellIndex = 0; cellIndex < cellResults.size(); ++cellIndex)
{
cvf::Vec2ui ij = ijFromCellIndex(cellIndex);
if (cellResults[cellIndex] != std::numeric_limits<double>::infinity())
{
distanceImage[ij.x()][ij.y()] = 1u;
validResBoundingBox.add(cvf::Vec3d(cellCenterPosition(ij.x(), ij.y()), 0.0));
}
else
{
mapCellVisibility[cellIndex] = false;
}
}
if (m_limitToPorePressureRegions && m_paddingAroundPorePressureRegion > 0.0)
{
RiaImageTools::distanceTransform2d(distanceImage);
cvf::Vec3d porExtent = validResBoundingBox.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;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RimGeoMechContourMapProjection::retrieveParameterWeights()
{
return std::vector<double>();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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();
size_t nCells = numberOfCells();
std::vector<double> aggregatedResults = std::vector<double>(nCells, std::numeric_limits<double>::infinity());
if (!resultAddress.isValid())
return aggregatedResults;
if (resultAddress.fieldName == "PP")
{
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.
}
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)
{
if (mapCellVisibility.empty() || mapCellVisibility[index])
{
cvf::Vec2ui ij = ijFromCellIndex(index);
aggregatedResults[index] = calculateValueInMapCell(ij.x(), ij.y(), resultValues);
}
}
return aggregatedResults;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimGeoMechContourMapProjection::resultVariableChanged() const
{
RimGeoMechCellColors* cellColors = view()->cellResult();
RigFemResultAddress resAddr = cellColors->resultAddress();
if (resAddr.fieldName == "PP")
{
resAddr.fieldName = "POR-Bar"; // More likely to be in memory than POR
}
if (resAddr.fieldName == "POR-Bar") resAddr.resultPosType = RIG_ELEMENT_NODAL;
return !(m_currentResultAddr == resAddr);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimGeoMechContourMapProjection::clearResultVariable()
{
m_currentResultAddr = RigFemResultAddress();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimGridView* RimGeoMechContourMapProjection::baseView() const
{
return view();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<size_t> RimGeoMechContourMapProjection::findIntersectingCells(const cvf::BoundingBox& bbox) const
{
std::vector<size_t> allCellIndices;
m_femPart->findIntersectingCellsWithExistingSearchTree(bbox, &allCellIndices);
return allCellIndices;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RimGeoMechContourMapProjection::kLayer(size_t globalCellIdx) const
{
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;
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);
return overlapVolume;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimGeoMechContourMapProjection::calculateRayLengthInCell(size_t globalCellIdx,
const cvf::Vec3d& highestPoint,
const cvf::Vec3d& lowestPoint) 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]]);
std::vector<HexIntersectionInfo> intersections;
if (RigHexIntersectionTools::lineHexCellIntersection(highestPoint, lowestPoint, hexCorners.data(), 0, &intersections))
{
double lengthInCell = (intersections.back().m_intersectionPoint - intersections.front().m_intersectionPoint).length();
return lengthInCell;
}
return 0.0;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimGeoMechContourMapProjection::getParameterWeightForCell(size_t globalCellIdx,
const std::vector<double>& parameterWeights) const
{
if (parameterWeights.empty()) return 1.0;
return parameterWeights[globalCellIdx];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<double> RimGeoMechContourMapProjection::gridCellValues(RigFemResultAddress resAddr, std::vector<float>& resultValues) const
{
std::vector<double> gridCellValues(m_femPart->elementCount(), std::numeric_limits<double>::infinity());
for (size_t globalCellIdx = 0; globalCellIdx < m_femPart->elementCount(); ++globalCellIdx)
{
RigElementType elmType = m_femPart->elementType(globalCellIdx);
if (!(elmType == HEX8 || elmType == HEX8P)) continue;
if (resAddr.resultPosType == RIG_ELEMENT)
{
gridCellValues[globalCellIdx] = static_cast<double>(resultValues[globalCellIdx]);
}
else if (resAddr.resultPosType == RIG_ELEMENT_NODAL)
{
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 gridCellValues;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimGeoMechCase* RimGeoMechContourMapProjection::geoMechCase() const
{
RimGeoMechCase* geoMechCase = nullptr;
firstAncestorOrThisOfType(geoMechCase);
return geoMechCase;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimGeoMechContourMapView* RimGeoMechContourMapProjection::view() const
{
RimGeoMechContourMapView* view = nullptr;
firstAncestorOrThisOfTypeAsserted(view);
return view;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimGeoMechContourMapProjection::fieldChangedByUi(const caf::PdmFieldHandle* changedField,
const QVariant& oldValue,
const QVariant& newValue)
{
RimContourMapProjection::fieldChangedByUi(changedField, oldValue, newValue);
if (changedField == &m_limitToPorePressureRegions || changedField == &m_paddingAroundPorePressureRegion)
{
clearGridMapping();
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
QList<caf::PdmOptionItemInfo>
RimGeoMechContourMapProjection::calculateValueOptions(const caf::PdmFieldHandle* fieldNeedingOptions, bool* useOptionsOnly)
{
QList<caf::PdmOptionItemInfo> options;
if (fieldNeedingOptions == &m_resultAggregation)
{
std::vector<ResultAggregationEnum> validOptions = {RESULTS_TOP_VALUE,
RESULTS_MEAN_VALUE,
RESULTS_GEOM_VALUE,
RESULTS_HARM_VALUE,
RESULTS_MIN_VALUE,
RESULTS_MAX_VALUE,
RESULTS_SUM};
for (ResultAggregationEnum option : validOptions)
{
options.push_back(caf::PdmOptionItemInfo(ResultAggregation::uiText(option), option));
}
}
return options;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimGeoMechContourMapProjection::defineUiOrdering(QString uiConfigName, caf::PdmUiOrdering& uiOrdering)
{
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;
myAttr->m_delaySliderUpdateUntilRelease = true;
}
}
}
Reference in New Issue View Git Blame Copy Permalink