diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp index ecaf785128..94ff59cae5 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp @@ -1380,8 +1380,9 @@ RigFemClosestResultIndexCalculator::RigFemClosestResultIndexCalculator(RigFemPar { m_resultIndexToClosestResult = -1; m_closestNodeId = -1; + m_closestElementNodeResIdx = -1; - if ( resultPosition != RIG_ELEMENT_NODAL_FACE ) + if ( resultPosition != RIG_ELEMENT_NODAL_FACE || m_face == -1 ) { RigElementType elmType = femPart->elementType(elementIndex); const int* elmentConn = femPart->connectivities(elementIndex); @@ -1406,13 +1407,19 @@ RigFemClosestResultIndexCalculator::RigFemClosestResultIndexCalculator(RigFemPar { float scalarValue = std::numeric_limits::infinity(); int nodeIdx = elmentConn[closestLocalNode]; + m_closestElementNodeResIdx = static_cast(femPart->elementNodeResultIdx(elementIndex, closestLocalNode)); + if ( resultPosition == RIG_NODAL ) { m_resultIndexToClosestResult = nodeIdx; } + else if (resultPosition == RIG_ELEMENT_NODAL_FACE) + { + m_resultIndexToClosestResult = -1; + } else { - m_resultIndexToClosestResult = static_cast(femPart->elementNodeResultIdx(elementIndex, closestLocalNode)); + m_resultIndexToClosestResult = m_closestElementNodeResIdx; } m_closestNodeId = femPart->nodes().nodeIds[nodeIdx]; @@ -1424,6 +1431,7 @@ RigFemClosestResultIndexCalculator::RigFemClosestResultIndexCalculator(RigFemPar int closestNodeIdx = -1; { int closestLocFaceNode = -1; + int closestLocalElmNode = -1; { RigElementType elmType = femPart->elementType(elementIndex); int faceCount = RigFemTypes::elmentFaceCount(elmType); @@ -1441,6 +1449,7 @@ RigFemClosestResultIndexCalculator::RigFemClosestResultIndexCalculator(RigFemPar { closestLocFaceNode = faceNodIdx; closestNodeIdx = nodeIdx; + closestLocalElmNode = localElmNodeIndicesForFace[faceNodIdx]; minDist = dist; } } @@ -1449,7 +1458,11 @@ RigFemClosestResultIndexCalculator::RigFemClosestResultIndexCalculator(RigFemPar int elmNodFaceResIdxElmStart = elementIndex * 24; // HACK should get from part int elmNodFaceResIdxFaceStart = elmNodFaceResIdxElmStart + 4*m_face; - if ( closestLocFaceNode >= 0 ) elmNodFaceResIdx = elmNodFaceResIdxFaceStart + closestLocFaceNode; + if ( closestLocFaceNode >= 0 ) + { + elmNodFaceResIdx = elmNodFaceResIdxFaceStart + closestLocFaceNode; + m_closestElementNodeResIdx = static_cast(femPart->elementNodeResultIdx(elementIndex, closestLocalElmNode)); + } } m_resultIndexToClosestResult = elmNodFaceResIdx; diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h index 2eb50de9c7..3b9e656cee 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h @@ -138,8 +138,10 @@ public: int resultIndexToClosestResult() { return m_resultIndexToClosestResult; } int closestNodeId() { return m_closestNodeId; } + int closestElementNodeResIdx () { return m_closestElementNodeResIdx; } private: int m_resultIndexToClosestResult; int m_closestNodeId; + int m_closestElementNodeResIdx; }; diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.cpp index ef9e28c1e9..6f8009def6 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.cpp @@ -25,6 +25,7 @@ #include "RigFemPartResultsCollection.h" #include "RigFemTypes.h" #include "RigGeoMechCaseData.h" +#include "RiuGeoMechXfTensorResultAccessor.h" #include // Needed for HUGE_VAL on Linux @@ -34,15 +35,38 @@ RigFemTimeHistoryResultAccessor::RigFemTimeHistoryResultAccessor(RigGeoMechCaseData* geomData, RigFemResultAddress femResultAddress, size_t gridIndex, - size_t cellIndex, + int elementIndex, int face, const cvf::Vec3d& intersectionPoint) : m_geoMechCaseData(geomData), m_femResultAddress(femResultAddress), m_gridIndex(gridIndex), - m_cellIndex(cellIndex), + m_elementIndex(elementIndex), m_face(face), - m_intersectionPoint(intersectionPoint) + m_intersectionPoint(intersectionPoint), + m_hasIntersectionTriangle(false) +{ + computeTimeHistoryData(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemTimeHistoryResultAccessor::RigFemTimeHistoryResultAccessor(RigGeoMechCaseData* geomData, + RigFemResultAddress femResultAddress, + size_t gridIndex, + int elementIndex, + int face, + const cvf::Vec3d& intersectionPoint, + const std::array& intersectionTriangle) + : m_geoMechCaseData(geomData), + m_femResultAddress(femResultAddress), + m_gridIndex(gridIndex), + m_elementIndex(elementIndex), + m_face(face), + m_intersectionPoint(intersectionPoint), + m_hasIntersectionTriangle(true), + m_intersectionTriangle(intersectionTriangle) { computeTimeHistoryData(); } @@ -57,13 +81,13 @@ QString RigFemTimeHistoryResultAccessor::topologyText() const if (m_geoMechCaseData) { RigFemPart* femPart = m_geoMechCaseData->femParts()->part(m_gridIndex); - int elementId = femPart->elmId(m_cellIndex); + int elementId = femPart->elmId(m_elementIndex); text += QString("Element : Id[%1]").arg(elementId); size_t i = 0; size_t j = 0; size_t k = 0; - if (m_geoMechCaseData->femParts()->part(m_gridIndex)->structGrid()->ijkFromCellIndex(m_cellIndex, &i, &j, &k)) + if (m_geoMechCaseData->femParts()->part(m_gridIndex)->structGrid()->ijkFromCellIndex(m_elementIndex, &i, &j, &k)) { // Adjust to 1-based Eclipse indexing i++; @@ -100,23 +124,40 @@ void RigFemTimeHistoryResultAccessor::computeTimeHistoryData() RigFemClosestResultIndexCalculator closestCalc(m_geoMechCaseData->femParts()->part(m_gridIndex), m_femResultAddress.resultPosType, - m_cellIndex, + m_elementIndex, m_face, m_intersectionPoint ); int scalarResultIndex = closestCalc.resultIndexToClosestResult(); + m_closestNodeId = closestCalc.closestNodeId(); - if (scalarResultIndex < 0) return; - RigFemPartResultsCollection* femPartResultsColl = m_geoMechCaseData->femPartResults(); - for (int frameIdx = 0; frameIdx < femPartResultsColl->frameCount(); frameIdx++) - { - const std::vector& scalarResults = m_geoMechCaseData->femPartResults()->resultValues(m_femResultAddress, static_cast(m_gridIndex), frameIdx); - if (scalarResults.size()) - { - float scalarValue = scalarResults[scalarResultIndex]; + if (m_femResultAddress.resultPosType == RIG_ELEMENT_NODAL_FACE && m_hasIntersectionTriangle) + { + int closestElmNodeResIndex = closestCalc.closestElementNodeResIdx(); + + for ( int frameIdx = 0; frameIdx < femPartResultsColl->frameCount(); frameIdx++ ) + { + RiuGeoMechXfTensorResultAccessor stressXfAccessor(femPartResultsColl, m_femResultAddress, frameIdx); + float scalarValue = stressXfAccessor.calculateElmNodeValue(m_intersectionTriangle, closestElmNodeResIndex); m_timeHistoryValues.push_back(scalarValue); } } + else + { + if ( scalarResultIndex < 0 ) return; + + for ( int frameIdx = 0; frameIdx < femPartResultsColl->frameCount(); frameIdx++ ) + { + const std::vector& scalarResults = m_geoMechCaseData->femPartResults()->resultValues(m_femResultAddress, static_cast(m_gridIndex), frameIdx); + if ( scalarResults.size() ) + { + float scalarValue = scalarResults[scalarResultIndex]; + + m_timeHistoryValues.push_back(scalarValue); + } + } + } + } diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.h index 7cd7090490..b1e5e78311 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.h @@ -23,6 +23,7 @@ #include "cvfStructGrid.h" #include "cvfVector3.h" +#include class RigGeoMechCaseData; @@ -33,12 +34,21 @@ public: RigFemTimeHistoryResultAccessor(RigGeoMechCaseData* geomData, RigFemResultAddress femResultAddress, size_t gridIndex, - size_t cellIndex, + int elementIndex, int face, const cvf::Vec3d& intersectionPoint); + RigFemTimeHistoryResultAccessor(RigGeoMechCaseData* geomData, + RigFemResultAddress femResultAddress, + size_t gridIndex, + int elementIndex, + int face, + const cvf::Vec3d& intersectionPoint, + const std::array& m_intersectionTriangle); + QString topologyText() const; std::vector timeHistoryValues() const; + int closestNodeId() const { return m_closestNodeId; } private: void computeTimeHistoryData(); @@ -48,12 +58,16 @@ private: RigFemResultAddress m_femResultAddress; size_t m_gridIndex; - size_t m_cellIndex; + int m_elementIndex; size_t m_scalarResultIndex; int m_face; + int m_closestNodeId; cvf::Vec3d m_intersectionPoint; + bool m_hasIntersectionTriangle; + std::array m_intersectionTriangle; + std::vector m_timeHistoryValues; }; diff --git a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.cpp b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.cpp index 9c5202243c..0dbb75a0fc 100644 --- a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.cpp +++ b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.cpp @@ -40,6 +40,19 @@ const std::vector& RivIntersectionBoxSourceInfo::triangleToCellIndex() c return m_intersectionBoxGeometryGenerator->triangleToCellIndex(); } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::array RivIntersectionBoxSourceInfo::triangle(int triangleIdx) const +{ + std::array tri; + tri[0] = (*m_intersectionBoxGeometryGenerator->triangleVxes())[triangleIdx*3]; + tri[1] = (*m_intersectionBoxGeometryGenerator->triangleVxes())[triangleIdx*3+1]; + tri[2] = (*m_intersectionBoxGeometryGenerator->triangleVxes())[triangleIdx*3+2]; + + return tri; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.h b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.h index 7718c35893..b8a9901dda 100644 --- a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.h +++ b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.h @@ -21,6 +21,7 @@ #include "cvfBase.h" #include "cvfObject.h" #include "cvfArray.h" +#include class RivIntersectionBoxGeometryGenerator; class RimIntersectionBox; @@ -32,6 +33,7 @@ public: const std::vector& triangleToCellIndex() const; + std::array triangle(int triangleIdx) const; const RimIntersectionBox* intersectionBox() const; private: diff --git a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp index 7b24a9c1d4..ca80711b39 100644 --- a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp +++ b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp @@ -52,6 +52,7 @@ #include "cvfRenderStatePoint.h" #include "cafTensor3.h" #include "cvfGeometryTools.h" +#include "RiuGeoMechXfTensorResultAccessor.h" //-------------------------------------------------------------------------------------------------- @@ -252,8 +253,6 @@ void RivIntersectionPartMgr::calculateGeoMechTextureCoords(cvf::Vec2fArray* text } } - - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -264,98 +263,28 @@ void RivIntersectionPartMgr::calculateGeoMechTensorXfTextureCoords(cvf::Vec2fArr const RigFemResultAddress& resVarAddress, int timeStepIdx, const cvf::ScalarMapper* mapper) -{ - RigFemResultAddress tensComp = resVarAddress; - tensComp.resultPosType = RIG_ELEMENT_NODAL; +{ - tensComp.componentName = "S11"; - const std::vector& tens11 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx); - tensComp.componentName = "S22"; - const std::vector& tens22 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx); - tensComp.componentName = "S33"; - const std::vector& tens33 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx); - tensComp.componentName = "S12"; - const std::vector& tens12 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx); - tensComp.componentName = "S23"; - const std::vector& tens23 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx); - tensComp.componentName = "S13"; - const std::vector& tens13 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx); + RiuGeoMechXfTensorResultAccessor accessor(caseData->femPartResults(), resVarAddress, timeStepIdx); textureCoords->resize(vertexWeights.size()); - - caf::Ten3f::TensorComponentEnum resultComponent = caf::Ten3f::SZZ; + cvf::Vec2f* rawPtr = textureCoords->ptr(); + int vxCount = static_cast(vertexWeights.size()); + int triCount = vxCount/3; - if ( resVarAddress.componentName == "SN" ) resultComponent = caf::Ten3f::SZZ; - if ( resVarAddress.componentName == "STH" ) resultComponent = caf::Ten3f::SXX; - if ( resVarAddress.componentName == "STQV" ) resultComponent = caf::Ten3f::SYY; - if ( resVarAddress.componentName == "TNH" ) resultComponent = caf::Ten3f::SZX; - if ( resVarAddress.componentName == "TNQV" ) resultComponent = caf::Ten3f::SYZ; - if ( resVarAddress.componentName == "THQV" ) resultComponent = caf::Ten3f::SXY; - - if(tens11.size() == 0) + #pragma omp parallel for schedule(dynamic) + for ( int triangleIdx = 0; triangleIdx < triCount; ++triangleIdx ) { - textureCoords->setAll(cvf::Vec2f(0.0, 1.0f)); + int triangleVxStartIdx = triangleIdx*3; + float values[3]; + + accessor.calculateInterpolatedValue(&((*triangelVertices)[triangleVxStartIdx]), &(vertexWeights[triangleVxStartIdx]), values ); + + rawPtr[triangleVxStartIdx + 0] = (values[0] != std::numeric_limits::infinity()) ? mapper->mapToTextureCoord(values[0]) : cvf::Vec2f(0.0f, 1.0f); + rawPtr[triangleVxStartIdx + 1] = (values[1] != std::numeric_limits::infinity()) ? mapper->mapToTextureCoord(values[1]) : cvf::Vec2f(0.0f, 1.0f); + rawPtr[triangleVxStartIdx + 2] = (values[2] != std::numeric_limits::infinity()) ? mapper->mapToTextureCoord(values[2]) : cvf::Vec2f(0.0f, 1.0f); } - else - { - cvf::Vec2f* rawPtr = textureCoords->ptr(); - int vxCount = static_cast(vertexWeights.size()); - int triCount = vxCount/3; - - #pragma omp parallel for schedule(dynamic) - for (int triangleIdx = 0; triangleIdx < triCount; ++triangleIdx) - { - int triangleVxStartIdx = triangleIdx*3; - cvf::Vec3f p0 = triangelVertices->get(triangleVxStartIdx); - cvf::Vec3f p1 = triangelVertices->get(triangleVxStartIdx + 1); - cvf::Vec3f p2 = triangelVertices->get(triangleVxStartIdx + 2); - - cvf::Mat3f triangleXf = cvf::GeometryTools::computePlaneHorizontalRotationMx(p1 - p0, p2 - p0); - - for(int triangleVxIdx = triangleVxStartIdx; triangleVxIdx < triangleVxStartIdx+3; ++triangleVxIdx) - { - float ipT11 = 0; - float ipT22 = 0; - float ipT33 = 0; - float ipT12 = 0; - float ipT23 = 0; - float ipT13 = 0; - - int weightCount = vertexWeights[triangleVxIdx].size(); - for(int wIdx = 0; wIdx < weightCount; ++wIdx) - { - size_t resIdx = vertexWeights[triangleVxIdx].vxId(wIdx) ; - float interpolationWeight = vertexWeights[triangleVxIdx].weight(wIdx); - ipT11 += tens11[resIdx] * interpolationWeight; - ipT22 += tens22[resIdx] * interpolationWeight; - ipT33 += tens33[resIdx] * interpolationWeight; - ipT12 += tens12[resIdx] * interpolationWeight; - ipT23 += tens23[resIdx] * interpolationWeight; - ipT13 += tens13[resIdx] * interpolationWeight; - } - - if ( ipT11 == HUGE_VAL || ipT11 != ipT11 - || ipT22 == HUGE_VAL || ipT22 != ipT22 - || ipT33 == HUGE_VAL || ipT33 != ipT33 - || ipT12 == HUGE_VAL || ipT12 != ipT12 - || ipT23 == HUGE_VAL || ipT23 != ipT23 - || ipT13 == HUGE_VAL || ipT13 != ipT13 ) // a != a is true for NAN's - { - rawPtr[triangleVxIdx][1] = 1.0f; - } - else - { - caf::Ten3f tensor(ipT11, ipT22, ipT33, - ipT12, ipT23, ipT13); - caf::Ten3f xfTen = tensor.rotated(triangleXf); - - - rawPtr[triangleVxIdx] = mapper->mapToTextureCoord(xfTen[resultComponent]); - } - } - } - } } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.cpp b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.cpp index 8ba3c18b99..470c5a63db 100644 --- a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.cpp +++ b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.cpp @@ -41,6 +41,19 @@ const std::vector& RivIntersectionSourceInfo::triangleToCellIndex() cons return m_crossSectionGeometryGenerator->triangleToCellIndex(); } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::array RivIntersectionSourceInfo::triangle(int triangleIdx) const +{ + std::array tri; + tri[0] = (*m_crossSectionGeometryGenerator->triangleVxes())[triangleIdx*3]; + tri[1] = (*m_crossSectionGeometryGenerator->triangleVxes())[triangleIdx*3+1]; + tri[2] = (*m_crossSectionGeometryGenerator->triangleVxes())[triangleIdx*3+2]; + + return tri; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.h b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.h index 5cd5aafeb3..436a0a6c75 100644 --- a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.h +++ b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.h @@ -22,6 +22,7 @@ #include "cvfBase.h" #include "cvfObject.h" #include "cvfArray.h" +#include class RivIntersectionGeometryGenerator; class RimIntersection; @@ -32,7 +33,7 @@ public: RivIntersectionSourceInfo(RivIntersectionGeometryGenerator* geometryGenerator); const std::vector& triangleToCellIndex() const; - + std::array triangle(int triangleIdx) const; const RimIntersection* crossSection() const; private: diff --git a/ApplicationCode/UserInterface/CMakeLists_files.cmake b/ApplicationCode/UserInterface/CMakeLists_files.cmake index bc6b0877b0..1167e10147 100644 --- a/ApplicationCode/UserInterface/CMakeLists_files.cmake +++ b/ApplicationCode/UserInterface/CMakeLists_files.cmake @@ -36,6 +36,7 @@ ${CEE_CURRENT_LIST_DIR}RiuViewer.h ${CEE_CURRENT_LIST_DIR}RiuViewerCommands.h ${CEE_CURRENT_LIST_DIR}RiuWellLogPlot.h ${CEE_CURRENT_LIST_DIR}RiuWellLogTrack.h +${CEE_CURRENT_LIST_DIR}RiuGeoMechXfTensorResultAccessor.h ) set (SOURCE_GROUP_SOURCE_FILES @@ -70,6 +71,8 @@ ${CEE_CURRENT_LIST_DIR}RiuViewer.cpp ${CEE_CURRENT_LIST_DIR}RiuViewerCommands.cpp ${CEE_CURRENT_LIST_DIR}RiuWellLogPlot.cpp ${CEE_CURRENT_LIST_DIR}RiuWellLogTrack.cpp +${CEE_CURRENT_LIST_DIR}RiuGeoMechXfTensorResultAccessor.cpp + ) list(APPEND CODE_HEADER_FILES diff --git a/ApplicationCode/UserInterface/RiuFemResultTextBuilder.cpp b/ApplicationCode/UserInterface/RiuFemResultTextBuilder.cpp index e3d5967de7..3e5f8ee9f8 100644 --- a/ApplicationCode/UserInterface/RiuFemResultTextBuilder.cpp +++ b/ApplicationCode/UserInterface/RiuFemResultTextBuilder.cpp @@ -29,13 +29,18 @@ #include "RimGeoMechCase.h" #include "RimGeoMechResultDefinition.h" #include "RimGeoMechView.h" +#include "RiuGeoMechXfTensorResultAccessor.h" //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -RiuFemResultTextBuilder::RiuFemResultTextBuilder(RimGeoMechView* reservoirView, int gridIndex, int cellIndex, int timeStepIndex) +RiuFemResultTextBuilder::RiuFemResultTextBuilder(RimGeoMechView* reservoirView, + int gridIndex, + int cellIndex, + int timeStepIndex) + : m_isIntersectionTriangleSet(false) { CVF_ASSERT(reservoirView); @@ -56,6 +61,15 @@ void RiuFemResultTextBuilder::setIntersectionPoint(cvf::Vec3d intersectionPoint) m_intersectionPoint = intersectionPoint; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RiuFemResultTextBuilder::setIntersectionTriangle(const std::array& triangle) +{ + m_intersectionTriangle = triangle; + m_isIntersectionTriangleSet = true; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -341,7 +355,8 @@ QString RiuFemResultTextBuilder::closestNodeResultText(RimGeoMechResultDefinitio m_intersectionPoint); int resultIndex = closestIndexCalc.resultIndexToClosestResult(); int closestNodeId = closestIndexCalc.closestNodeId(); - + int closestElmNodResIdx = closestIndexCalc.closestElementNodeResIdx(); + float scalarValue = (resultIndex >= 0) ? scalarResults[resultIndex]: std::numeric_limits::infinity(); @@ -356,6 +371,16 @@ QString RiuFemResultTextBuilder::closestNodeResultText(RimGeoMechResultDefinitio .arg(caf::AppEnum::textFromIndex(m_face)) .arg(scalarValue)); } + else if (m_isIntersectionTriangleSet && activeResultPosition == RIG_ELEMENT_NODAL_FACE) + { + RiuGeoMechXfTensorResultAccessor tensAccessor(geomData->femPartResults(), resultColors->resultAddress(), m_timeStepIndex); + float tensValue = tensAccessor.calculateElmNodeValue(m_intersectionTriangle, closestElmNodResIdx); + + text.append(QString("Closest result: N[%1], in Element[%2] transformed onto intersection: %3 \n") + .arg(closestNodeId) + .arg(femPart->elmId(m_cellIndex)) + .arg(tensValue)); + } } } diff --git a/ApplicationCode/UserInterface/RiuFemResultTextBuilder.h b/ApplicationCode/UserInterface/RiuFemResultTextBuilder.h index 13a9a6d926..2b3e3d45b7 100644 --- a/ApplicationCode/UserInterface/RiuFemResultTextBuilder.h +++ b/ApplicationCode/UserInterface/RiuFemResultTextBuilder.h @@ -24,6 +24,7 @@ #include "cvfStructGrid.h" #include +#include class RigGeoMechCaseData; class RimEclipseCellColors; @@ -44,6 +45,7 @@ public: RiuFemResultTextBuilder(RimGeoMechView* reservoirView, int gridIndex, int cellIndex, int timeStepIndex); void setFace(int face); void setIntersectionPoint(cvf::Vec3d intersectionPoint); + void setIntersectionTriangle(const std::array& triangle); QString mainResultText(); @@ -67,6 +69,8 @@ private: int m_timeStepIndex; int m_face; + bool m_isIntersectionTriangleSet; + std::array m_intersectionTriangle; cvf::Vec3d m_intersectionPoint; }; diff --git a/ApplicationCode/UserInterface/RiuGeoMechXfTensorResultAccessor.cpp b/ApplicationCode/UserInterface/RiuGeoMechXfTensorResultAccessor.cpp new file mode 100644 index 0000000000..0f4f14db12 --- /dev/null +++ b/ApplicationCode/UserInterface/RiuGeoMechXfTensorResultAccessor.cpp @@ -0,0 +1,147 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) Statoil 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 "RiuGeoMechXfTensorResultAccessor.h" +#include "RigFemPartResultsCollection.h" +#include "cvfGeometryTools.h" +#include "RivHexGridIntersectionTools.h" + + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RiuGeoMechXfTensorResultAccessor::RiuGeoMechXfTensorResultAccessor(RigFemPartResultsCollection * femResCollection, const RigFemResultAddress& resVarAddress, int timeStepIdx) +{ + RigFemResultAddress tensComp = resVarAddress; + tensComp.resultPosType = RIG_ELEMENT_NODAL; + + tensComp.componentName = "S11"; + tens11 = &femResCollection->resultValues(tensComp, 0, timeStepIdx); + tensComp.componentName = "S22"; + tens22 = &femResCollection->resultValues(tensComp, 0, timeStepIdx); + tensComp.componentName = "S33"; + tens33 = &femResCollection->resultValues(tensComp, 0, timeStepIdx); + tensComp.componentName = "S12"; + tens12 = &femResCollection->resultValues(tensComp, 0, timeStepIdx); + tensComp.componentName = "S23"; + tens23 = &femResCollection->resultValues(tensComp, 0, timeStepIdx); + tensComp.componentName = "S13"; + tens13 = &femResCollection->resultValues(tensComp, 0, timeStepIdx); + + resultComponent = caf::Ten3f::SZZ; + + if ( resVarAddress.componentName == "SN" ) resultComponent = caf::Ten3f::SZZ; + if ( resVarAddress.componentName == "STH" ) resultComponent = caf::Ten3f::SXX; + if ( resVarAddress.componentName == "STQV" ) resultComponent = caf::Ten3f::SYY; + if ( resVarAddress.componentName == "TNH" ) resultComponent = caf::Ten3f::SZX; + if ( resVarAddress.componentName == "TNQV" ) resultComponent = caf::Ten3f::SYZ; + if ( resVarAddress.componentName == "THQV" ) resultComponent = caf::Ten3f::SXY; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RiuGeoMechXfTensorResultAccessor::calculateInterpolatedValue(const cvf::Vec3f triangle[3], const RivIntersectionVertexWeights vertexWeights[3], float returnValues[3]) +{ + if ( tens11->size() == 0 ) + { + returnValues[0] = returnValues[1] = returnValues[2] = std::numeric_limits::infinity(); + return; + } + + cvf::Mat3f triangleXf = cvf::GeometryTools::computePlaneHorizontalRotationMx(triangle[1] - triangle[0], triangle[2] - triangle[0]); + + for ( int triangleVxIdx = 0; triangleVxIdx < 3; ++triangleVxIdx ) + { + float ipT11 = 0; + float ipT22 = 0; + float ipT33 = 0; + float ipT12 = 0; + float ipT23 = 0; + float ipT13 = 0; + + int weightCount = vertexWeights[triangleVxIdx].size(); + for ( int wIdx = 0; wIdx < weightCount; ++wIdx ) + { + size_t resIdx = vertexWeights[triangleVxIdx].vxId(wIdx) ; + float interpolationWeight = vertexWeights[triangleVxIdx].weight(wIdx); + ipT11 += (*tens11)[resIdx] * interpolationWeight; + ipT22 += (*tens22)[resIdx] * interpolationWeight; + ipT33 += (*tens33)[resIdx] * interpolationWeight; + ipT12 += (*tens12)[resIdx] * interpolationWeight; + ipT23 += (*tens23)[resIdx] * interpolationWeight; + ipT13 += (*tens13)[resIdx] * interpolationWeight; + } + + if ( ipT11 == HUGE_VAL || ipT11 != ipT11 + || ipT22 == HUGE_VAL || ipT22 != ipT22 + || ipT33 == HUGE_VAL || ipT33 != ipT33 + || ipT12 == HUGE_VAL || ipT12 != ipT12 + || ipT23 == HUGE_VAL || ipT23 != ipT23 + || ipT13 == HUGE_VAL || ipT13 != ipT13 ) // a != a is true for NAN's + { + returnValues[triangleVxIdx] = std::numeric_limits::infinity(); + } + else + { + caf::Ten3f tensor(ipT11, ipT22, ipT33, + ipT12, ipT23, ipT13); + caf::Ten3f xfTen = tensor.rotated(triangleXf); + + returnValues[triangleVxIdx] = xfTen[resultComponent]; + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +float RiuGeoMechXfTensorResultAccessor::calculateElmNodeValue(const std::array & triangle, int globalElmNodeResIndex) +{ + if ( tens11->size() == 0 ) return std::numeric_limits::infinity(); + + cvf::Mat3f triangleXf = cvf::GeometryTools::computePlaneHorizontalRotationMx(triangle[1] - triangle[0], triangle[2] - triangle[0]); + + float ipT11 = (*tens11)[globalElmNodeResIndex]; + float ipT22 = (*tens22)[globalElmNodeResIndex]; + float ipT33 = (*tens33)[globalElmNodeResIndex]; + float ipT12 = (*tens12)[globalElmNodeResIndex]; + float ipT23 = (*tens23)[globalElmNodeResIndex]; + float ipT13 = (*tens13)[globalElmNodeResIndex]; + + if ( ipT11 == HUGE_VAL || ipT11 != ipT11 + || ipT22 == HUGE_VAL || ipT22 != ipT22 + || ipT33 == HUGE_VAL || ipT33 != ipT33 + || ipT12 == HUGE_VAL || ipT12 != ipT12 + || ipT23 == HUGE_VAL || ipT23 != ipT23 + || ipT13 == HUGE_VAL || ipT13 != ipT13 ) // a != a is true for NAN's + { + return std::numeric_limits::infinity(); + } + else + { + caf::Ten3f tensor(ipT11, ipT22, ipT33, + ipT12, ipT23, ipT13); + caf::Ten3f xfTen = tensor.rotated(triangleXf); + + float scalarValue = xfTen[resultComponent]; + + return scalarValue; + } +} diff --git a/ApplicationCode/UserInterface/RiuGeoMechXfTensorResultAccessor.h b/ApplicationCode/UserInterface/RiuGeoMechXfTensorResultAccessor.h new file mode 100644 index 0000000000..69b052821c --- /dev/null +++ b/ApplicationCode/UserInterface/RiuGeoMechXfTensorResultAccessor.h @@ -0,0 +1,51 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) Statoil 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. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once +#include +#include + +#include "cafTensor3.h" + +class RigFemPartResultsCollection; +class RigFemResultAddress; +class RivIntersectionVertexWeights; + +class RiuGeoMechXfTensorResultAccessor +{ +public: + RiuGeoMechXfTensorResultAccessor(RigFemPartResultsCollection * femResCollection, const RigFemResultAddress& resVarAddress, int timeStepIdx); + + void calculateInterpolatedValue(const cvf::Vec3f triangle[3], const RivIntersectionVertexWeights vertexWeights[3], float returnValues[3]); + + float calculateElmNodeValue(const std::array & triangle, int globalElmNodeResIndex); + + +private: + + const std::vector* tens11; + const std::vector* tens22; + const std::vector* tens33; + const std::vector* tens12; + const std::vector* tens23; + const std::vector* tens13; + + caf::Ten3f::TensorComponentEnum resultComponent; +}; + + diff --git a/ApplicationCode/UserInterface/RiuSelectionChangedHandler.cpp b/ApplicationCode/UserInterface/RiuSelectionChangedHandler.cpp index 79400988c0..7c203d1ecc 100644 --- a/ApplicationCode/UserInterface/RiuSelectionChangedHandler.cpp +++ b/ApplicationCode/UserInterface/RiuSelectionChangedHandler.cpp @@ -44,6 +44,8 @@ #include #include +#include "RigFemPartResultsCollection.h" +#include "RigFemPartCollection.h" //-------------------------------------------------------------------------------------------------- /// @@ -143,12 +145,29 @@ void RiuSelectionChangedHandler::addCurveFromSelectionItem(const RiuGeoMechSelec geoMechView->geoMechCase() && geoMechView->geoMechCase()->geoMechData()) { - RigFemTimeHistoryResultAccessor timeHistResultAccessor(geoMechView->geoMechCase()->geoMechData(), - geoMechView->cellResultResultDefinition()->resultAddress(), - geomSelectionItem->m_gridIndex, - geomSelectionItem->m_cellIndex, - geomSelectionItem->m_elementFace, - geomSelectionItem->m_localIntersectionPoint); + std::unique_ptr timeHistResultAccessor; + + if ( geomSelectionItem->m_hasIntersectionTriangle ) + { + timeHistResultAccessor = std::unique_ptr( + new RigFemTimeHistoryResultAccessor(geoMechView->geoMechCase()->geoMechData(), + geoMechView->cellResultResultDefinition()->resultAddress(), + geomSelectionItem->m_gridIndex, + static_cast(geomSelectionItem->m_cellIndex), + geomSelectionItem->m_elementFace, + geomSelectionItem->m_localIntersectionPoint, + geomSelectionItem->m_intersectionTriangle)); + } + else + { + timeHistResultAccessor = std::unique_ptr( + new RigFemTimeHistoryResultAccessor(geoMechView->geoMechCase()->geoMechData(), + geoMechView->cellResultResultDefinition()->resultAddress(), + geomSelectionItem->m_gridIndex, + static_cast(geomSelectionItem->m_cellIndex), + geomSelectionItem->m_elementFace, + geomSelectionItem->m_localIntersectionPoint)); + } QString curveName; curveName.append(geoMechView->geoMechCase()->caseUserDescription() + ", "); @@ -158,16 +177,22 @@ void RiuSelectionChangedHandler::addCurveFromSelectionItem(const RiuGeoMechSelec curveName.append(geoMechView->cellResultResultDefinition()->resultFieldUiName()+ ", ") ; curveName.append(geoMechView->cellResultResultDefinition()->resultComponentUiName() + " "); - if ( resPosAppEnum == RIG_ELEMENT_NODAL_FACE && geomSelectionItem->m_elementFace >= 0 ) - { - curveName.append(", " + caf::AppEnum::textFromIndex(geomSelectionItem->m_elementFace)); + if ( resPosAppEnum == RIG_ELEMENT_NODAL_FACE ) + { + if ( geomSelectionItem->m_elementFace >= 0 ) + { + curveName.append(", " + caf::AppEnum::textFromIndex(geomSelectionItem->m_elementFace)); + } + else + { + curveName.append(", from N[" + QString::number(timeHistResultAccessor->closestNodeId()) + "] transformed onto intersection"); + } } + curveName.append("\n"); - curveName.append(":\n"); + curveName.append(timeHistResultAccessor->topologyText()); - curveName.append(timeHistResultAccessor.topologyText()); - - std::vector timeHistoryValues = timeHistResultAccessor.timeHistoryValues(); + std::vector timeHistoryValues = timeHistResultAccessor->timeHistoryValues(); QStringList stepNames = geoMechView->geoMechCase()->timeStepStrings(); std::vector dates = RimGeoMechCase::dateTimeVectorFromTimeStepStrings(stepNames); @@ -258,6 +283,8 @@ void RiuSelectionChangedHandler::updateResultInfo(const RiuSelectionItem* itemAd RiuFemResultTextBuilder textBuilder(geomView, (int)geomSelectionItem->m_gridIndex, (int)geomSelectionItem->m_cellIndex, geomView->currentTimeStep()); textBuilder.setIntersectionPoint(geomSelectionItem->m_localIntersectionPoint); textBuilder.setFace(geomSelectionItem->m_elementFace); + if (geomSelectionItem->m_hasIntersectionTriangle) textBuilder.setIntersectionTriangle(geomSelectionItem->m_intersectionTriangle); + resultInfo = textBuilder.mainResultText(); pickInfo = textBuilder.topologyText(", "); diff --git a/ApplicationCode/UserInterface/RiuSelectionManager.cpp b/ApplicationCode/UserInterface/RiuSelectionManager.cpp index e89ca92e5a..379fa14277 100644 --- a/ApplicationCode/UserInterface/RiuSelectionManager.cpp +++ b/ApplicationCode/UserInterface/RiuSelectionManager.cpp @@ -140,6 +140,29 @@ RiuGeoMechSelectionItem::RiuGeoMechSelectionItem(RimGeoMechView* view, m_cellIndex(cellIndex), m_color(color), m_elementFace(elementFace), - m_localIntersectionPoint(localIntersectionPoint) + m_localIntersectionPoint(localIntersectionPoint), + m_hasIntersectionTriangle(false) { } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RiuGeoMechSelectionItem::RiuGeoMechSelectionItem(RimGeoMechView* view, + size_t gridIndex, + size_t cellIndex, + cvf::Color3f color, + int elementFace, + const cvf::Vec3d& localIntersectionPoint, + const std::array& intersectionTriangle) + : m_view(view), + m_gridIndex(gridIndex), + m_cellIndex(cellIndex), + m_color(color), + m_elementFace(elementFace), + m_localIntersectionPoint(localIntersectionPoint), + m_hasIntersectionTriangle(true), + m_intersectionTriangle(m_intersectionTriangle) +{ + +} diff --git a/ApplicationCode/UserInterface/RiuSelectionManager.h b/ApplicationCode/UserInterface/RiuSelectionManager.h index 769af31fff..4b91fc8608 100644 --- a/ApplicationCode/UserInterface/RiuSelectionManager.h +++ b/ApplicationCode/UserInterface/RiuSelectionManager.h @@ -27,6 +27,7 @@ #include #include +#include class RimEclipseView; class RiuSelectionChangedHandler; @@ -138,6 +139,14 @@ public: cvf::Color3f color, int elementFace, const cvf::Vec3d& localIntersectionPoint); + + explicit RiuGeoMechSelectionItem(RimGeoMechView* view, + size_t gridIndex, + size_t cellIndex, + cvf::Color3f color, + int elementFace, + const cvf::Vec3d& localIntersectionPoint, + const std::array& m_intersectionTriangle ); virtual ~RiuGeoMechSelectionItem() {}; virtual RiuSelectionType type() const @@ -151,6 +160,8 @@ public: size_t m_cellIndex; cvf::Color3f m_color; int m_elementFace; + bool m_hasIntersectionTriangle; + std::array m_intersectionTriangle; cvf::Vec3d m_localIntersectionPoint; }; diff --git a/ApplicationCode/UserInterface/RiuViewerCommands.cpp b/ApplicationCode/UserInterface/RiuViewerCommands.cpp index d157bc626a..6b30f75ebf 100644 --- a/ApplicationCode/UserInterface/RiuViewerCommands.cpp +++ b/ApplicationCode/UserInterface/RiuViewerCommands.cpp @@ -80,6 +80,7 @@ #include #include #include +#include @@ -470,6 +471,8 @@ void RiuViewerCommands::handlePickAction(int winPosX, int winPosY, Qt::KeyboardM size_t nncIndex = cvf::UNDEFINED_SIZE_T; cvf::StructGridInterface::FaceType face = cvf::StructGridInterface::NO_FACE; int gmFace = -1; + bool intersectionHit = false; + std::array intersectionTriangleHit; cvf::Vec3d localIntersectionPoint(cvf::Vec3d::ZERO); @@ -525,14 +528,20 @@ void RiuViewerCommands::handlePickAction(int winPosX, int winPosY, Qt::KeyboardM else if (crossSectionSourceInfo) { findCellAndGridIndex(crossSectionSourceInfo, firstPartTriangleIndex, &cellIndex, &gridIndex); + intersectionHit = true; + intersectionTriangleHit = crossSectionSourceInfo->triangle(firstPartTriangleIndex); RiuMainWindow::instance()->selectAsCurrentItem(const_cast(crossSectionSourceInfo->crossSection())); + } else if (intersectionBoxSourceInfo) { findCellAndGridIndex(intersectionBoxSourceInfo, firstPartTriangleIndex, &cellIndex, &gridIndex); + intersectionHit = true; + intersectionTriangleHit = intersectionBoxSourceInfo->triangle(firstPartTriangleIndex); RiuMainWindow::instance()->selectAsCurrentItem(const_cast(intersectionBoxSourceInfo->intersectionBox())); + } else if (eclipseWellSourceInfo) { @@ -580,9 +589,10 @@ void RiuViewerCommands::handlePickAction(int winPosX, int winPosY, Qt::KeyboardM } RimGeoMechView* geomView = dynamic_cast(m_reservoirView.p()); - if (geomView) + if (geomView ) { - selItem = new RiuGeoMechSelectionItem(geomView, gridIndex, cellIndex, curveColor, gmFace, localIntersectionPoint); + if(intersectionHit) selItem = new RiuGeoMechSelectionItem(geomView, gridIndex, cellIndex, curveColor, gmFace, localIntersectionPoint, intersectionTriangleHit); + else selItem = new RiuGeoMechSelectionItem(geomView, gridIndex, cellIndex, curveColor, gmFace, localIntersectionPoint); } }