OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

#933 Fixed the closes result for intersections also. Display the tensor data from the closest element node transformed onto the plane.

Browse Source
This commit is contained in:
Jacob Støren 2016-10-26 17:05:39 +02:00
parent 1a23e6c670
commit c93c7e41df
18 changed files with 454 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp
View File

@ -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<float>::infinity();
int nodeIdx = elmentConn[closestLocalNode];
m_closestElementNodeResIdx = static_cast<int>(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<int>(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<int>(femPart->elementNodeResultIdx(elementIndex, closestLocalElmNode));
}
}
m_resultIndexToClosestResult = elmNodFaceResIdx;

2
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h
View File

@ -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;
};

69
ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.cpp
View File

@ -25,6 +25,7 @@
#include "RigFemPartResultsCollection.h"
#include "RigFemTypes.h"
#include "RigGeoMechCaseData.h"
#include "RiuGeoMechXfTensorResultAccessor.h"
#include <cmath> // 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<cvf::Vec3f, 3>& 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<float>& scalarResults = m_geoMechCaseData->femPartResults()->resultValues(m_femResultAddress, static_cast<int>(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<float>& scalarResults = m_geoMechCaseData->femPartResults()->resultValues(m_femResultAddress, static_cast<int>(m_gridIndex), frameIdx);
if ( scalarResults.size() )
{
float scalarValue = scalarResults[scalarResultIndex];
m_timeHistoryValues.push_back(scalarValue);
}
}
}
}

18
ApplicationCode/GeoMech/GeoMechDataModel/RigFemTimeHistoryResultAccessor.h
View File

@ -23,6 +23,7 @@
#include "cvfStructGrid.h"
#include "cvfVector3.h"
#include <array>
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<cvf::Vec3f, 3>& m_intersectionTriangle);
QString topologyText() const;
std::vector<double> 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<cvf::Vec3f, 3> m_intersectionTriangle;
std::vector<double> m_timeHistoryValues;
};

13
ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.cpp
View File

@ -40,6 +40,19 @@ const std::vector<size_t>& RivIntersectionBoxSourceInfo::triangleToCellIndex() c
return m_intersectionBoxGeometryGenerator->triangleToCellIndex();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::array<cvf::Vec3f, 3> RivIntersectionBoxSourceInfo::triangle(int triangleIdx) const
{
std::array<cvf::Vec3f, 3> 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;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

2
ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxSourceInfo.h
View File

@ -21,6 +21,7 @@
#include "cvfBase.h"
#include "cvfObject.h"
#include "cvfArray.h"
#include <array>
class RivIntersectionBoxGeometryGenerator;
class RimIntersectionBox;
@ -32,6 +33,7 @@ public:
const std::vector<size_t>& triangleToCellIndex() const;
std::array<cvf::Vec3f, 3> triangle(int triangleIdx) const;
const RimIntersectionBox* intersectionBox() const;
private:

103
ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp
View File

@ -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<float>& tens11 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx);
tensComp.componentName = "S22";
const std::vector<float>& tens22 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx);
tensComp.componentName = "S33";
const std::vector<float>& tens33 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx);
tensComp.componentName = "S12";
const std::vector<float>& tens12 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx);
tensComp.componentName = "S23";
const std::vector<float>& tens23 = caseData->femPartResults()->resultValues(tensComp, 0, timeStepIdx);
tensComp.componentName = "S13";
const std::vector<float>& 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<int>(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<float>::infinity()) ? mapper->mapToTextureCoord(values[0]) : cvf::Vec2f(0.0f, 1.0f);
rawPtr[triangleVxStartIdx + 1] = (values[1] != std::numeric_limits<float>::infinity()) ? mapper->mapToTextureCoord(values[1]) : cvf::Vec2f(0.0f, 1.0f);
rawPtr[triangleVxStartIdx + 2] = (values[2] != std::numeric_limits<float>::infinity()) ? mapper->mapToTextureCoord(values[2]) : cvf::Vec2f(0.0f, 1.0f);
}
else
{
cvf::Vec2f* rawPtr = textureCoords->ptr();
int vxCount = static_cast<int>(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]);
}
}
}
}
}
//--------------------------------------------------------------------------------------------------

13
ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.cpp
View File

@ -41,6 +41,19 @@ const std::vector<size_t>& RivIntersectionSourceInfo::triangleToCellIndex() cons
return m_crossSectionGeometryGenerator->triangleToCellIndex();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::array<cvf::Vec3f, 3> RivIntersectionSourceInfo::triangle(int triangleIdx) const
{
std::array<cvf::Vec3f, 3> 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;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

3
ApplicationCode/ModelVisualization/Intersections/RivIntersectionSourceInfo.h
View File

@ -22,6 +22,7 @@
#include "cvfBase.h"
#include "cvfObject.h"
#include "cvfArray.h"
#include <array>
class RivIntersectionGeometryGenerator;
class RimIntersection;
@ -32,7 +33,7 @@ public:
RivIntersectionSourceInfo(RivIntersectionGeometryGenerator* geometryGenerator);
const std::vector<size_t>& triangleToCellIndex() const;
std::array<cvf::Vec3f, 3> triangle(int triangleIdx) const;
const RimIntersection* crossSection() const;
private:

3
ApplicationCode/UserInterface/CMakeLists_files.cmake
View File

@ -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

29
ApplicationCode/UserInterface/RiuFemResultTextBuilder.cpp
View File

@ -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<cvf::Vec3f, 3>& 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<float>::infinity();
@ -356,6 +371,16 @@ QString RiuFemResultTextBuilder::closestNodeResultText(RimGeoMechResultDefinitio
.arg(caf::AppEnum<cvf::StructGridInterface::FaceType>::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));
}
}
}

4
ApplicationCode/UserInterface/RiuFemResultTextBuilder.h
View File

@ -24,6 +24,7 @@
#include "cvfStructGrid.h"
#include <QString>
#include <array>
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<cvf::Vec3f, 3>& triangle);
QString mainResultText();
@ -67,6 +69,8 @@ private:
int m_timeStepIndex;
int m_face;
bool m_isIntersectionTriangleSet;
std::array<cvf::Vec3f, 3> m_intersectionTriangle;
cvf::Vec3d m_intersectionPoint;
};

147
ApplicationCode/UserInterface/RiuGeoMechXfTensorResultAccessor.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/gpl.html>
// 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<float>::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<float>::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<cvf::Vec3f, 3> & triangle, int globalElmNodeResIndex)
{
if ( tens11->size() == 0 ) return std::numeric_limits<float>::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<float>::infinity();
}
else
{
caf::Ten3f tensor(ipT11, ipT22, ipT33,
ipT12, ipT23, ipT13);
caf::Ten3f xfTen = tensor.rotated(triangleXf);
float scalarValue = xfTen[resultComponent];
return scalarValue;
}
}

51
ApplicationCode/UserInterface/RiuGeoMechXfTensorResultAccessor.h Normal file
View File

@ -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 <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#pragma once
#include <vector>
#include <array>
#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<cvf::Vec3f, 3> & triangle, int globalElmNodeResIndex);
private:
const std::vector<float>* tens11;
const std::vector<float>* tens22;
const std::vector<float>* tens33;
const std::vector<float>* tens12;
const std::vector<float>* tens23;
const std::vector<float>* tens13;
caf::Ten3f::TensorComponentEnum resultComponent;
};

53
ApplicationCode/UserInterface/RiuSelectionChangedHandler.cpp
View File

@ -44,6 +44,8 @@
#include <QStatusBar>
#include <assert.h>
#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<RigFemTimeHistoryResultAccessor> timeHistResultAccessor;
if ( geomSelectionItem->m_hasIntersectionTriangle )
{
timeHistResultAccessor = std::unique_ptr<RigFemTimeHistoryResultAccessor>(
new RigFemTimeHistoryResultAccessor(geoMechView->geoMechCase()->geoMechData(),
geoMechView->cellResultResultDefinition()->resultAddress(),
geomSelectionItem->m_gridIndex,
static_cast<int>(geomSelectionItem->m_cellIndex),
geomSelectionItem->m_elementFace,
geomSelectionItem->m_localIntersectionPoint,
geomSelectionItem->m_intersectionTriangle));
}
else
{
timeHistResultAccessor = std::unique_ptr<RigFemTimeHistoryResultAccessor>(
new RigFemTimeHistoryResultAccessor(geoMechView->geoMechCase()->geoMechData(),
geoMechView->cellResultResultDefinition()->resultAddress(),
geomSelectionItem->m_gridIndex,
static_cast<int>(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<cvf::StructGridInterface::FaceType>::textFromIndex(geomSelectionItem->m_elementFace));
if ( resPosAppEnum == RIG_ELEMENT_NODAL_FACE )
{
if ( geomSelectionItem->m_elementFace >= 0 )
{
curveName.append(", " + caf::AppEnum<cvf::StructGridInterface::FaceType>::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<double> timeHistoryValues = timeHistResultAccessor.timeHistoryValues();
std::vector<double> timeHistoryValues = timeHistResultAccessor->timeHistoryValues();
QStringList stepNames = geoMechView->geoMechCase()->timeStepStrings();
std::vector<QDateTime> 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(", ");

25
ApplicationCode/UserInterface/RiuSelectionManager.cpp
View File

@ -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<cvf::Vec3f, 3>& 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)
{
}

11
ApplicationCode/UserInterface/RiuSelectionManager.h
View File

@ -27,6 +27,7 @@
#include <vector>
#include <assert.h>
#include <array>
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<cvf::Vec3f, 3>& 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<cvf::Vec3f, 3> m_intersectionTriangle;
cvf::Vec3d m_localIntersectionPoint;
};

14
ApplicationCode/UserInterface/RiuViewerCommands.cpp
View File

@ -80,6 +80,7 @@
#include <QMenu>
#include <QMouseEvent>
#include <QStatusBar>
#include <array>
@ -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<cvf::Vec3f, 3> 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<RimIntersection*>(crossSectionSourceInfo->crossSection()));
}
else if (intersectionBoxSourceInfo)
{
findCellAndGridIndex(intersectionBoxSourceInfo, firstPartTriangleIndex, &cellIndex, &gridIndex);
intersectionHit = true;
intersectionTriangleHit = intersectionBoxSourceInfo->triangle(firstPartTriangleIndex);
RiuMainWindow::instance()->selectAsCurrentItem(const_cast<RimIntersectionBox*>(intersectionBoxSourceInfo->intersectionBox()));
}
else if (eclipseWellSourceInfo)
{
@ -580,9 +589,10 @@ void RiuViewerCommands::handlePickAction(int winPosX, int winPosY, Qt::KeyboardM
}
RimGeoMechView* geomView = dynamic_cast<RimGeoMechView*>(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);
}
}