///////////////////////////////////////////////////////////////////////////////// // // 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 // 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 #include #include #include 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_applyPPRegionLimitVertically, "VerticalLimit", false, "Apply Limit Vertically", "", "", ""); CAF_PDM_InitField(&m_paddingAroundPorePressureRegion, "PaddingAroundPorRegion", 0.0, "Horizontal 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 minmaxValAllTimeSteps = minmaxValuesAllTimeSteps(); legendConfig()->setAutomaticRanges(minmaxValAllTimeSteps.first, minmaxValAllTimeSteps.second, minVal, maxVal); QString projectionLegendText = QString("Map Projection\n%1").arg(m_resultAggregation().uiText()); if (cellColors->resultAddress().isValid()) { projectionLegendText += QString("\nResult: %1").arg(cellColors->resultFieldUiName()); if (!cellColors->resultComponentUiName().isEmpty()) { projectionLegendText += QString(", %1").arg(cellColors->resultComponentUiName()); } } else { projectionLegendText += QString("\nNo Result Selected"); } legendConfig()->setTitle(projectionLegendText); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- cvf::ref RimGeoMechContourMapProjection::getCellVisibility() const { cvf::ref cellGridIdxVisibility = new cvf::UByteArray(m_femPart->elementCount()); RivFemElmVisibilityCalculator::computeAllVisible(cellGridIdxVisibility.p(), m_femPart.p()); if (view()->rangeFilterCollection()->isActive()) { cvf::CellRangeFilter cellRangeFilter; view()->rangeFilterCollection()->compoundCellRangeFilter(&cellRangeFilter, 0); RivFemElmVisibilityCalculator::computeRangeVisibility(cellGridIdxVisibility.p(), m_femPart.p(), cellRangeFilter); } if (view()->propertyFilterCollection()->isActive()) { 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& 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::infinity(); if (validPorValue) { std::array 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 * m_paddingAroundPorePressureRegion(); boxMin.y() -= boxExtent.y() * 0.5 * m_paddingAroundPorePressureRegion(); boxMax.x() += boxExtent.x() * 0.5 * m_paddingAroundPorePressureRegion(); boxMax.y() += boxExtent.y() * 0.5 * m_paddingAroundPorePressureRegion(); 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_expandedBoundingBox = calculateExpandedPorBarBBox(view()->currentTimeStep()); } else { m_expandedBoundingBox = m_gridBoundingBox; } cvf::Vec3d minExpandedPoint = m_expandedBoundingBox.min() - cvf::Vec3d(gridEdgeOffset(), gridEdgeOffset(), 0.0); cvf::Vec3d maxExpandedPoint = m_expandedBoundingBox.max() + cvf::Vec3d(gridEdgeOffset(), gridEdgeOffset(), 0.0); if (m_limitToPorePressureRegions && !m_applyPPRegionLimitVertically) { minExpandedPoint.z() = m_gridBoundingBox.min().z(); maxExpandedPoint.z() = m_gridBoundingBox.max().z(); } 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 RimGeoMechContourMapProjection::getMapCellVisibility() { cvf::Vec2ui nCellsIJ = numberOfElementsIJ(); std::vector> distanceImage (nCellsIJ.x(), std::vector(nCellsIJ.y(), 0u)); std::vector mapCellVisibility; RigFemResultAddress resAddr = view()->cellResult()->resultAddress(); if (m_limitToPorePressureRegions) { resAddr = RigFemResultAddress(RigFemResultPosEnum::RIG_ELEMENT_NODAL, "POR-Bar", ""); } std::vector 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::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 RimGeoMechContourMapProjection::retrieveParameterWeights() { return std::vector(); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RimGeoMechContourMapProjection::generateResults(int timeStep) { RimGeoMechCellColors* cellColors = view()->cellResult(); RigFemResultAddress resultAddress = cellColors->resultAddress(); std::vector aggregatedResults = generateResultsFromAddress(resultAddress, m_mapCellVisibility, timeStep); return aggregatedResults; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RimGeoMechContourMapProjection::generateResultsFromAddress(RigFemResultAddress resultAddress, const std::vector& mapCellVisibility, int timeStep) { RigGeoMechCaseData* caseData = geoMechCase()->geoMechData(); RigFemPartResultsCollection* resultCollection = caseData->femPartResults(); size_t nCells = numberOfCells(); std::vector aggregatedResults = std::vector(nCells, std::numeric_limits::infinity()); bool wasInvalid = false; if (!resultAddress.isValid()) { wasInvalid = true; resultAddress = RigFemResultAddress(RigFemResultPosEnum::RIG_ELEMENT_NODAL, "POR-Bar", ""); } 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 resultValuesF = resultCollection->resultValues(resultAddress, 0, timeStep); std::vector resultValues = gridCellValues(resultAddress, resultValuesF); if (wasInvalid) { // For invalid result addresses we just use the POR-Bar result to get the reservoir region // And display a dummy 0-result in the region. for (double& value : resultValues) { if (value != std::numeric_limits::infinity()) { value = 0.0; } } } #pragma omp parallel for for (int index = 0; index < static_cast(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 RimGeoMechContourMapProjection::findIntersectingCells(const cvf::BoundingBox& bbox) const { std::vector 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 hexCorners; m_femPartGrid->cellCornerVertices(globalCellIdx, hexCorners.data()); cvf::BoundingBox overlapBBox; std::array overlapCorners; if (RigCellGeometryTools::estimateHexOverlapWithBoundingBox(hexCorners, bbox, &overlapCorners, &overlapBBox)) { double overlapVolume = RigCellGeometryTools::calculateCellVolume(overlapCorners); return overlapVolume; } return 0.0; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- double RimGeoMechContourMapProjection::calculateRayLengthInCell(size_t globalCellIdx, const cvf::Vec3d& highestPoint, const cvf::Vec3d& lowestPoint) const { std::array hexCorners; const std::vector& 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 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& parameterWeights) const { if (parameterWeights.empty()) return 1.0; return parameterWeights[globalCellIdx]; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RimGeoMechContourMapProjection::gridCellValues(RigFemResultAddress resAddr, std::vector& resultValues) const { std::vector gridCellValues(m_femPart->elementCount(), std::numeric_limits::infinity()); for (size_t globalCellIdx = 0; globalCellIdx < static_cast(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(resultValues[globalCellIdx]); } else if (resAddr.resultPosType == RIG_ELEMENT_NODAL) { RiaWeightedMeanCalculator cellAverage; for (int i = 0; i < 8; ++i) { size_t gridResultValueIdx = m_femPart->resultValueIdxFromResultPosType(resAddr.resultPosType, static_cast(globalCellIdx), i); cellAverage.addValueAndWeight(resultValues[gridResultValueIdx], 1.0); } gridCellValues[globalCellIdx] = static_cast(cellAverage.weightedMean()); } else { RiaWeightedMeanCalculator 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(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_applyPPRegionLimitVertically || changedField == &m_paddingAroundPorePressureRegion) { clearGridMapping(); } } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- QList RimGeoMechContourMapProjection::calculateValueOptions(const caf::PdmFieldHandle* fieldNeedingOptions, bool* useOptionsOnly) { QList options; if (fieldNeedingOptions == &m_resultAggregation) { std::vector 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("Map Boundaries"); group->add(&m_limitToPorePressureRegions); group->add(&m_applyPPRegionLimitVertically); group->add(&m_paddingAroundPorePressureRegion); m_applyPPRegionLimitVertically.uiCapability()->setUiReadOnly(!m_limitToPorePressureRegions()); 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(attribute); if (myAttr) { myAttr->m_minimum = 0.0; myAttr->m_maximum = 2.0; myAttr->m_sliderTickCount = 4; myAttr->m_delaySliderUpdateUntilRelease = true; } } }