#3186 Fix crash when selecting element nodal (and some other types) for some GeoMech well log curves.

This commit is contained in:
Gaute Lindkvist 2018-08-07 12:38:48 +02:00
parent 8425ea96ad
commit c869ce4430
4 changed files with 64 additions and 20 deletions

View File

@ -386,3 +386,32 @@ size_t RigFemPart::elementNodeResultCount() const
return lastElmResultIdx + 1;
}
//--------------------------------------------------------------------------------------------------
/// Generate a sensible index into the result vector based on which result position type is used.
//--------------------------------------------------------------------------------------------------
size_t RigFemPart::resultValueIdxFromResultPosType(RigFemResultPosEnum resultPosType, int elementIdx, int elmLocalNodeIdx) const
{
if (resultPosType == RIG_ELEMENT || resultPosType == RIG_FORMATION_NAMES)
{
CVF_ASSERT(elementIdx < m_elementId.size());
return elementIdx;
}
size_t elementNodeResultIdx = this->elementNodeResultIdx(static_cast<int>(elementIdx), elmLocalNodeIdx);
CVF_ASSERT(elementNodeResultIdx < elementNodeResultCount());
if (resultPosType == RIG_ELEMENT_NODAL || resultPosType == RIG_INTEGRATION_POINT)
{
return elementNodeResultIdx;
}
else if (resultPosType == RIG_NODAL)
{
size_t nodeIdx = nodeIdxFromElementNodeResultIdx(elementNodeResultIdx);
CVF_ASSERT(nodeIdx < m_nodes.nodeIds.size());
return nodeIdx;
}
CVF_ASSERT(false);
return 0u;
}

View File

@ -22,6 +22,7 @@
#include <stdlib.h>
#include "RigFemTypes.h"
#include "RigFemResultPosEnum.h"
#include "cvfObject.h"
#include "cvfAssert.h"
#include "cvfBoundingBox.h"
@ -63,7 +64,7 @@ public:
size_t elementNodeResultIdx(int elementIdx, int elmLocalNodeIdx) const { return m_elementConnectivityStartIndices[elementIdx] + elmLocalNodeIdx;}
size_t elementNodeResultCount() const;
int nodeIdxFromElementNodeResultIdx(size_t elmNodeResultIdx) const { return m_allElementConnectivities[elmNodeResultIdx]; }
size_t resultValueIdxFromResultPosType(RigFemResultPosEnum resultPosType, int elementIdx, int elmLocalNodeIdx) const;
RigFemPartNodes& nodes() {return m_nodes;}
const RigFemPartNodes& nodes() const {return m_nodes;}

View File

@ -352,39 +352,42 @@ T RigGeoMechWellLogExtractor::interpolateGridResultValue(RigFemResultPosEnum r
if (cellFace == cvf::StructGridInterface::NO_FACE)
{
if (resultPosType == RIG_ELEMENT_NODAL_FACE)
{
return std::numeric_limits<T>::infinity(); // undefined value. ELEMENT_NODAL_FACE values are only defined on a face.
}
// TODO: Should interpolate within the whole hexahedron. This requires converting to locals coordinates.
// For now just pick the average value for the cell.
T sumOfVertexValues = gridResultValues[femPart->elementNodeResultIdx(static_cast<int>(elmIdx), 0)];
size_t gridResultValueIdx = femPart->resultValueIdxFromResultPosType(resultPosType, static_cast<int>(elmIdx), 0);
T sumOfVertexValues = gridResultValues[gridResultValueIdx];
for (int i = 1; i < 8; ++i)
{
sumOfVertexValues = sumOfVertexValues + gridResultValues[femPart->elementNodeResultIdx(static_cast<int>(elmIdx), i)];
gridResultValueIdx = femPart->resultValueIdxFromResultPosType(resultPosType, static_cast<int>(elmIdx), i);
sumOfVertexValues = sumOfVertexValues + gridResultValues[gridResultValueIdx];
}
return sumOfVertexValues * (1.0 / 8.0);
}
int faceNodeCount = 0;
const int* faceLocalIndices = RigFemTypes::localElmNodeIndicesForFace(elmType, cellFace, &faceNodeCount);
const int* elementLocalIndicesForFace = RigFemTypes::localElmNodeIndicesForFace(elmType, cellFace, &faceNodeCount);
const int* elmNodeIndices = femPart->connectivities(elmIdx);
cvf::Vec3d v0(nodeCoords[elmNodeIndices[faceLocalIndices[0]]]);
cvf::Vec3d v1(nodeCoords[elmNodeIndices[faceLocalIndices[1]]]);
cvf::Vec3d v2(nodeCoords[elmNodeIndices[faceLocalIndices[2]]]);
cvf::Vec3d v3(nodeCoords[elmNodeIndices[faceLocalIndices[3]]]);
cvf::Vec3d v0(nodeCoords[elmNodeIndices[elementLocalIndicesForFace[0]]]);
cvf::Vec3d v1(nodeCoords[elmNodeIndices[elementLocalIndicesForFace[1]]]);
cvf::Vec3d v2(nodeCoords[elmNodeIndices[elementLocalIndicesForFace[2]]]);
cvf::Vec3d v3(nodeCoords[elmNodeIndices[elementLocalIndicesForFace[3]]]);
std::vector<size_t> nodeResIdx(4, cvf::UNDEFINED_SIZE_T);
if (resultPosType == RIG_NODAL)
for (size_t i = 0; i < nodeResIdx.size(); ++i)
{
for (size_t i = 0; i < nodeResIdx.size(); ++i)
if (resultPosType == RIG_ELEMENT_NODAL_FACE)
{
nodeResIdx[i] = elmNodeIndices[faceLocalIndices[i]];
nodeResIdx[i] = gridResultIndexFace(elmIdx, cellFace, static_cast<int>(i));
}
}
else
{
for (size_t i = 0; i < nodeResIdx.size(); ++i)
else
{
nodeResIdx[i] = (size_t)femPart->elementNodeResultIdx((int)elmIdx, faceLocalIndices[i]);
nodeResIdx[i] = femPart->resultValueIdxFromResultPosType(resultPosType, static_cast<int>(elmIdx), elementLocalIndicesForFace[i]);
}
}
@ -398,12 +401,12 @@ T RigGeoMechWellLogExtractor::interpolateGridResultValue(RigFemResultPosEnum r
int nodeIndex = femPart->nodeIdxFromElementNodeResultIdx(nodeResIdx[i]);
const std::vector<int>& elements = femPart->elementsUsingNode(nodeIndex);
const std::vector<unsigned char>& localIndices = femPart->elementLocalIndicesForNode(nodeIndex);
size_t otherNodeResIdx = femPart->elementNodeResultIdx(elements[0], static_cast<int>(localIndices[0]));
T nodeResultValue = gridResultValues[otherNodeResIdx];
size_t otherGridResultValueIdx = femPart->resultValueIdxFromResultPosType(resultPosType, elements[0], static_cast<int>(localIndices[0]));
T nodeResultValue = gridResultValues[otherGridResultValueIdx];
for (size_t j = 1; j < elements.size(); ++j)
{
otherNodeResIdx = femPart->elementNodeResultIdx(elements[j], static_cast<int>(localIndices[j]));
nodeResultValue = nodeResultValue + gridResultValues[otherNodeResIdx];
otherGridResultValueIdx = femPart->resultValueIdxFromResultPosType(resultPosType, elements[j], static_cast<int>(localIndices[j]));
nodeResultValue = nodeResultValue + gridResultValues[otherGridResultValueIdx];
}
nodeResultValue = nodeResultValue * (1.0 / elements.size());
nodeResultValues.push_back(nodeResultValue);
@ -427,6 +430,15 @@ T RigGeoMechWellLogExtractor::interpolateGridResultValue(RigFemResultPosEnum r
return interpolatedValue;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigGeoMechWellLogExtractor::gridResultIndexFace(size_t elementIdx, cvf::StructGridInterface::FaceType cellFace, int faceLocalNodeIdx) const
{
CVF_ASSERT(cellFace != cvf::StructGridInterface::NO_FACE && faceLocalNodeIdx < 4);
return elementIdx * 24 + static_cast<int>(cellFace) * 4 + faceLocalNodeIdx;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

View File

@ -33,6 +33,7 @@
#include <vector>
class RigFemPart;
class RigFemResultAddress;
class RigGeoMechCaseData;
class RigWellPath;
@ -60,6 +61,7 @@ private:
template<typename T>
T interpolateGridResultValue(RigFemResultPosEnum resultPosType, const std::vector<T>& gridResultValues, int64_t intersectionIdx, bool averageNodeElementResults) const;
size_t gridResultIndexFace(size_t elementIdx, cvf::StructGridInterface::FaceType cellFace, int faceLocalNodeIdx) const;
void calculateIntersection();
std::vector<size_t> findCloseCells(const cvf::BoundingBox& bb);
virtual cvf::Vec3d calculateLengthInCell(size_t cellIndex,