(#557, #560, #561) WLP: Started to return a cleaner intersection list from the extractor.

Created a common baseclass to be able to use same intersection cleanup code.
This commit is contained in:
Jacob Støren 2015-10-14 08:08:13 +02:00
parent 8698291e7d
commit 6014312056
8 changed files with 347 additions and 79 deletions

View File

@ -19,6 +19,7 @@ ${CEE_CURRENT_LIST_DIR}RigCombTransResultAccessor.h
${CEE_CURRENT_LIST_DIR}RigCombMultResultAccessor.h
${CEE_CURRENT_LIST_DIR}RigResultModifier.h
${CEE_CURRENT_LIST_DIR}RigResultModifierFactory.h
${CEE_CURRENT_LIST_DIR}RigWellLogExtractor.h
${CEE_CURRENT_LIST_DIR}RigEclipseWellLogExtractor.h
${CEE_CURRENT_LIST_DIR}RigLocalGrid.h
${CEE_CURRENT_LIST_DIR}RigMainGrid.h
@ -52,6 +53,7 @@ ${CEE_CURRENT_LIST_DIR}RigCellEdgeResultAccessor.cpp
${CEE_CURRENT_LIST_DIR}RigCombTransResultAccessor.cpp
${CEE_CURRENT_LIST_DIR}RigCombMultResultAccessor.cpp
${CEE_CURRENT_LIST_DIR}RigResultModifierFactory.cpp
${CEE_CURRENT_LIST_DIR}RigWellLogExtractor.cpp
${CEE_CURRENT_LIST_DIR}RigEclipseWellLogExtractor.cpp
${CEE_CURRENT_LIST_DIR}RigLocalGrid.cpp
${CEE_CURRENT_LIST_DIR}RigMainGrid.cpp

View File

@ -33,7 +33,7 @@
//==================================================================================================
RigEclipseWellLogExtractor::RigEclipseWellLogExtractor(const RigCaseData* aCase, const RigWellPath* wellpath)
: m_caseData(aCase), m_wellPath(wellpath)
: m_caseData(aCase), RigWellLogExtractor(wellpath)
{
calculateIntersection();
}
@ -50,6 +50,8 @@ void RigEclipseWellLogExtractor::calculateIntersection()
if (!m_wellPath->m_wellPathPoints.size()) return ;
std::map<RigMDCellIdxEnterLeaveIntersectionSorterKey, HexIntersectionInfo > uniqueIntersections;
for (size_t wpp = 0; wpp < m_wellPath->m_wellPathPoints.size() - 1; ++wpp)
{
cvf::BoundingBox bb;
@ -83,8 +85,8 @@ void RigEclipseWellLogExtractor::calculateIntersection()
// Now, with all the intersections of this piece of line, we need to
// sort them in order, and set the measured depth and corresponding cell index
// map <WellPathDepthPoint, (CellIdx, intersectionPoint)>
std::map<WellPathDepthPoint, HexIntersectionInfo > sortedIntersections;
// Inserting the intersections in this map will remove identical intersections
// and sort them according to MD, CellIdx, Leave/enter
double md1 = m_wellPath->m_measuredDepths[wpp];
double md2 = m_wellPath->m_measuredDepths[wpp+1];
@ -92,11 +94,13 @@ void RigEclipseWellLogExtractor::calculateIntersection()
for (size_t intIdx = 0; intIdx < intersections.size(); ++intIdx)
{
if (!isCellFaceNormalsOut) intersections[intIdx].m_isIntersectionEntering = !intersections[intIdx].m_isIntersectionEntering ;
double lenghtAlongLineSegment1 = (intersections[intIdx].m_intersectionPoint - p1).length();
double lenghtAlongLineSegment2 = (p2 - intersections[intIdx].m_intersectionPoint).length();
double measuredDepthDiff = md2 - md1;
double lineLength = lenghtAlongLineSegment1 + lenghtAlongLineSegment2;
double measuredDepthOfPoint = 0.0;
double measuredDepthDiff = md2 - md1;
double lineLength = lenghtAlongLineSegment1 + lenghtAlongLineSegment2;
double measuredDepthOfPoint = 0.0;
if (lineLength > 0.00001)
{
measuredDepthOfPoint = md1 + measuredDepthDiff*lenghtAlongLineSegment1/(lineLength);
@ -105,14 +109,125 @@ void RigEclipseWellLogExtractor::calculateIntersection()
{
measuredDepthOfPoint = md1;
}
sortedIntersections.insert(std::make_pair(WellPathDepthPoint(measuredDepthOfPoint, intersections[intIdx].m_isIntersectionEntering), intersections[intIdx]));
uniqueIntersections.insert(std::make_pair(RigMDCellIdxEnterLeaveIntersectionSorterKey(measuredDepthOfPoint,
intersections[intIdx].m_hexIndex,
intersections[intIdx].m_isIntersectionEntering),
intersections[intIdx]));
}
}
{
// For same MD and same cell, remove enter/leave pairs, as they only touches the wellpath, and should not contribute.
std::map<RigMDCellIdxEnterLeaveIntersectionSorterKey, HexIntersectionInfo >::iterator it1 = uniqueIntersections.begin();
std::map<RigMDCellIdxEnterLeaveIntersectionSorterKey, HexIntersectionInfo >::iterator it2 = uniqueIntersections.begin();
std::vector<std::map<RigMDCellIdxEnterLeaveIntersectionSorterKey, HexIntersectionInfo >::iterator> iteratorsToIntersectonsToErase;
while (it2 != uniqueIntersections.end())
{
++it2;
if (it2 != uniqueIntersections.end())
{
if (RigHexIntersector::isEqualDepth(it1->first.measuredDepth, it2->first.measuredDepth))
{
if (it1->first.hexIndex == it2->first.hexIndex)
{
// Remove the two from the map, as they just are a touch of the cell surface
CVF_TIGHT_ASSERT(!it1->first.isEnteringCell && it2->first.isEnteringCell);
iteratorsToIntersectonsToErase.push_back(it1);
iteratorsToIntersectonsToErase.push_back(it2);
}
}
}
++it1;
}
// Now populate the return arrays
// Erase all the intersections that is not needed
for (size_t erItIdx = 0; erItIdx < iteratorsToIntersectonsToErase.size(); ++erItIdx)
{
uniqueIntersections.erase(iteratorsToIntersectonsToErase[erItIdx]);
}
}
std::map<WellPathDepthPoint, HexIntersectionInfo >::iterator it;
it = sortedIntersections.begin();
while (it != sortedIntersections.end())
// Copy the map into a different sorting regime, with enter leave more significant than cell index
std::map<RigMDEnterLeaveCellIdxIntersectionSorterKey, HexIntersectionInfo > sortedUniqueIntersections;
{
std::map<RigMDCellIdxEnterLeaveIntersectionSorterKey, HexIntersectionInfo >::iterator it = uniqueIntersections.begin();
while (it != uniqueIntersections.end())
{
sortedUniqueIntersections.insert(std::make_pair(RigMDEnterLeaveCellIdxIntersectionSorterKey(it->first.measuredDepth, it->first.isEnteringCell, it->first.hexIndex),
it->second));
++it;
}
}
// Make sure we have sensible pairs of intersections. One pair for each in/out of a cell
{
// Add an intersection for the well startpoint that is inside the first cell
std::map<RigMDEnterLeaveCellIdxIntersectionSorterKey, HexIntersectionInfo >::iterator it = sortedUniqueIntersections.begin();
if (it != sortedUniqueIntersections.end() && !it->first.isEnteringCell) // Leaving a cell as first intersection. Well starts inside a cell.
{
// Needs wellpath start point in front
HexIntersectionInfo firstLeavingPoint = it->second;
firstLeavingPoint.m_intersectionPoint = m_wellPath->m_wellPathPoints[0];
sortedUniqueIntersections.insert(std::make_pair(RigMDEnterLeaveCellIdxIntersectionSorterKey(m_wellPath->m_measuredDepths[0], firstLeavingPoint.m_hexIndex, true),
firstLeavingPoint));
}
// Add an intersection for the well endpoint possibly inside the last cell.
std::map<RigMDEnterLeaveCellIdxIntersectionSorterKey, HexIntersectionInfo >::reverse_iterator rit = sortedUniqueIntersections.rbegin();
if (rit != sortedUniqueIntersections.rend() && rit->first.isEnteringCell) // Entering a cell as last intersection. Well ends inside a cell.
{
// Needs wellpath end point at end
HexIntersectionInfo lastEnterPoint = rit->second;
lastEnterPoint.m_intersectionPoint = m_wellPath->m_wellPathPoints.back();
sortedUniqueIntersections.insert(std::make_pair(RigMDEnterLeaveCellIdxIntersectionSorterKey(m_wellPath->m_measuredDepths.back(), lastEnterPoint.m_hexIndex, false),
lastEnterPoint));
}
}
std::map<RigMDEnterLeaveIntersectionSorterKey, HexIntersectionInfo > filteredSortedIntersections;
{
std::map<RigMDEnterLeaveCellIdxIntersectionSorterKey, HexIntersectionInfo >::iterator it1 = sortedUniqueIntersections.begin();
std::map<RigMDEnterLeaveCellIdxIntersectionSorterKey, HexIntersectionInfo >::iterator it2;
while (it1 != sortedUniqueIntersections.end())
{
it2 = it1;
++it2;
if (it2 == sortedUniqueIntersections.end()) break;
// If we have a proper pair, insert in the filtered list and continue
if ( !RigMDEnterLeaveCellIdxIntersectionSorterKey::isProperPair( it1->first, it2->first))
{
//CVF_ASSERT(false);
cvf::Trace::show(cvf::String("Well log curve is inaccurate around MD: ") + cvf::String::number((double)(it1->first.measuredDepth)) + ", " + cvf::String::number((double)(it2->first.measuredDepth)));
}
{
filteredSortedIntersections.insert(std::make_pair(RigMDEnterLeaveIntersectionSorterKey(it1->first.measuredDepth, it1->first.isEnteringCell),
it1->second));
++it1;
filteredSortedIntersections.insert(std::make_pair(RigMDEnterLeaveIntersectionSorterKey(it1->first.measuredDepth, it1->first.isEnteringCell),
it1->second));
++it1;
}
}
}
{
// Now populate the return arrays
std::map<RigMDEnterLeaveIntersectionSorterKey, HexIntersectionInfo >::iterator it;
it = filteredSortedIntersections.begin();
while (it != filteredSortedIntersections.end())
{
m_measuredDepth.push_back(it->first.measuredDepth);
m_trueVerticalDepth.push_back(abs(it->second.m_intersectionPoint[2]));
@ -130,7 +245,7 @@ void RigEclipseWellLogExtractor::calculateIntersection()
void RigEclipseWellLogExtractor::curveData(const RigResultAccessor* resultAccessor, std::vector<double>* values)
{
CVF_TIGHT_ASSERT(values);
values->resize(m_intersections.size());// + 1); // Plus one for the end of the wellpath stopping inside a cell
values->resize(m_intersections.size());
for (size_t cpIdx = 0; cpIdx < m_intersections.size(); ++cpIdx)
{
@ -138,9 +253,7 @@ void RigEclipseWellLogExtractor::curveData(const RigResultAccessor* resultAccess
cvf::StructGridInterface::FaceType cellFace = m_intersectedCellFaces[cpIdx];
(*values)[cpIdx] = resultAccessor->cellFaceScalarGlobIdx(cellIdx, cellFace);
}
// What do we do with the endpoint of the wellpath ?
// Ignore it for now ...
}
//--------------------------------------------------------------------------------------------------

View File

@ -19,6 +19,8 @@
#pragma once
#include "RigWellLogExtractor.h"
#include "cvfBase.h"
#include "cvfObject.h"
#include "cvfMath.h"
@ -38,17 +40,12 @@ namespace cvf {
//==================================================================================================
///
//==================================================================================================
class RigEclipseWellLogExtractor : public cvf::Object
class RigEclipseWellLogExtractor : public RigWellLogExtractor
{
public:
RigEclipseWellLogExtractor(const RigCaseData* aCase, const RigWellPath* wellpath);
const std::vector<double>& measuredDepth(){ return m_measuredDepth; }
const std::vector<double>& trueVerticalDepth(){return m_trueVerticalDepth;}
const RigWellPath* wellPathData() { return m_wellPath.p();}
void curveData(const RigResultAccessor* resultAccessor, std::vector<double>* values );
const RigCaseData* caseData() { return m_caseData.p();}
protected:
@ -56,14 +53,4 @@ protected:
std::vector<size_t> findCloseCells(const cvf::BoundingBox& bb);
cvf::cref<RigCaseData> m_caseData;
std::vector<double> m_measuredDepth;
std::vector<double> m_trueVerticalDepth;
std::vector<cvf::Vec3d> m_intersections;
std::vector<size_t> m_intersectedCells;
std::vector<cvf::StructGridInterface::FaceType>
m_intersectedCellFaces;
cvf::cref<RigWellPath> m_wellPath;
};

View File

@ -32,7 +32,7 @@
///
//--------------------------------------------------------------------------------------------------
RigGeoMechWellLogExtractor::RigGeoMechWellLogExtractor(RigGeoMechCaseData* aCase, const RigWellPath* wellpath)
:m_caseData(aCase), m_wellPath(wellpath)
:m_caseData(aCase), RigWellLogExtractor(wellpath)
{
calculateIntersection();
}
@ -155,7 +155,7 @@ void RigGeoMechWellLogExtractor::calculateIntersection()
// sort them in order, and set the measured depth and corresponding cell index
// map <WellPathDepthPoint, (CellIdx, intersectionPoint)>
std::map<WellPathDepthPoint, HexIntersectionInfo > sortedIntersections;
std::map<RigMDEnterLeaveIntersectionSorterKey, HexIntersectionInfo > sortedIntersections;
double md1 = m_wellPath->m_measuredDepths[wpp];
double md2 = m_wellPath->m_measuredDepths[wpp+1];
@ -176,12 +176,12 @@ void RigGeoMechWellLogExtractor::calculateIntersection()
measuredDepthOfPoint = md1;
}
sortedIntersections.insert(std::make_pair(WellPathDepthPoint(measuredDepthOfPoint, intersections[intIdx].m_isIntersectionEntering), intersections[intIdx]));
sortedIntersections.insert(std::make_pair(RigMDEnterLeaveIntersectionSorterKey(measuredDepthOfPoint, intersections[intIdx].m_isIntersectionEntering), intersections[intIdx]));
}
// Now populate the return arrays
std::map<WellPathDepthPoint, HexIntersectionInfo >::iterator it;
std::map<RigMDEnterLeaveIntersectionSorterKey, HexIntersectionInfo >::iterator it;
it = sortedIntersections.begin();
while (it != sortedIntersections.end())
{

View File

@ -19,6 +19,8 @@
#pragma once
#include "RigWellLogExtractor.h"
#include "cvfBase.h"
#include "cvfObject.h"
#include "cvfMath.h"
@ -38,18 +40,12 @@ namespace cvf {
//==================================================================================================
///
//==================================================================================================
class RigGeoMechWellLogExtractor : public cvf::Object
class RigGeoMechWellLogExtractor : public RigWellLogExtractor
{
public:
RigGeoMechWellLogExtractor(RigGeoMechCaseData* aCase, const RigWellPath* wellpath);
const std::vector<double>& measuredDepth(){ return m_measuredDepth; }
const std::vector<double>& trueVerticalDepth(){return m_trueVerticalDepth;}
const RigWellPath* wellPathData() { return m_wellPath.p();}
void curveData(const RigFemResultAddress& resAddr, int frameIndex, std::vector<double>* values );
const RigGeoMechCaseData* caseData() { return m_caseData.p();}
private:
@ -57,16 +53,6 @@ private:
std::vector<size_t> findCloseCells(const cvf::BoundingBox& bb);
cvf::ref<RigGeoMechCaseData> m_caseData;
std::vector<double> m_measuredDepth;
std::vector<double> m_trueVerticalDepth;
std::vector<cvf::Vec3d> m_intersections;
std::vector<size_t> m_intersectedCells;
std::vector<cvf::StructGridInterface::FaceType>
m_intersectedCellFaces;
cvf::cref<RigWellPath> m_wellPath;
};

View File

@ -72,39 +72,53 @@ struct RigHexIntersector
int next = i < 3 ? i+1 : 0;
int intsStatus = cvf::GeometryTools::intersectLineSegmentTriangle(p1, p2,
hexCorners[faceVertexIndices[i]], hexCorners[faceVertexIndices[next]], faceCenter,
hexCorners[faceVertexIndices[i]],
hexCorners[faceVertexIndices[next]],
faceCenter,
&intersection,
&isEntering);
if (intsStatus == 1)
{
intersectionCount++;
intersections->push_back(HexIntersectionInfo(intersection, isEntering, static_cast<cvf::StructGridInterface::FaceType>(face), hexIndex));
intersections->push_back(HexIntersectionInfo(intersection,
isEntering,
static_cast<cvf::StructGridInterface::FaceType>(face), hexIndex));
}
}
}
return intersectionCount;
}
static bool isEqualDepth(double d1, double d2)
{
double depthDiff = d1 - d2;
const double tolerance = 0.1;// Meters To handle inaccuracies across faults
return (fabs(depthDiff) < tolerance); // Equal depth
}
};
//==================================================================================================
/// Class used to sort the intersections along the wellpath
//==================================================================================================
struct WellPathDepthPoint
struct RigMDEnterLeaveIntersectionSorterKey
{
WellPathDepthPoint(double md, bool entering): measuredDepth(md), isEnteringCell(entering){}
RigMDEnterLeaveIntersectionSorterKey(double md, bool entering): measuredDepth(md), isEnteringCell(entering){}
double measuredDepth;
bool isEnteringCell; // As opposed to leaving.
bool isLeavingCell() { return !isEnteringCell;}
bool operator < (const WellPathDepthPoint& other) const
// Sorting according to MD first, then Leaving before entering cell
bool operator < (const RigMDEnterLeaveIntersectionSorterKey& other) const
{
double depthDiff = other.measuredDepth - measuredDepth;
const double tolerance = 1e-6;
if (fabs(depthDiff) < tolerance) // Equal depth
if (RigHexIntersector::isEqualDepth(measuredDepth, other.measuredDepth) )
{
if (isEnteringCell == other.isEnteringCell)
{
@ -112,29 +126,99 @@ struct WellPathDepthPoint
return false;
}
if(!isEnteringCell) // Leaving this cell
{
return true;
}
else
{
return false;
}
return !isEnteringCell; // Sort Leaving cell before (less than) entering cell
}
// The depths are not equal
if (measuredDepth < other.measuredDepth)
{
return true;
}
else if (measuredDepth > other.measuredDepth)
{
return false;
}
CVF_ASSERT(false);
return false;
return (measuredDepth < other.measuredDepth);
}
};
//==================================================================================================
/// Class used to sort the intersections along the wellpath,
/// Use as key in a map
/// Sorting according to MD first, then Cell Idx, then Leaving before entering cell
//==================================================================================================
struct RigMDCellIdxEnterLeaveIntersectionSorterKey
{
RigMDCellIdxEnterLeaveIntersectionSorterKey(double md, size_t cellIdx, bool entering): measuredDepth(md), hexIndex(cellIdx), isEnteringCell(entering){}
double measuredDepth;
size_t hexIndex;
bool isEnteringCell; // As opposed to leaving.
bool isLeavingCell() const { return !isEnteringCell;}
bool operator < (const RigMDCellIdxEnterLeaveIntersectionSorterKey& other) const
{
if (RigHexIntersector::isEqualDepth(measuredDepth, other.measuredDepth))
{
if (hexIndex == other.hexIndex)
{
if (isEnteringCell == other.isEnteringCell)
{
// Completely equal: intersections at cell edges or corners or edges of the face triangles
return false;
}
return !isEnteringCell; // Sort Leaving cell before (less than) entering cell
}
return (hexIndex < other.hexIndex);
}
// The depths are not equal
return (measuredDepth < other.measuredDepth);
}
};
//==================================================================================================
/// Class used to sort the intersections along the wellpath,
/// Use as key in a map
/// Sorting according to MD first,then Leaving before entering cell, then Cell Idx,
//==================================================================================================
struct RigMDEnterLeaveCellIdxIntersectionSorterKey
{
RigMDEnterLeaveCellIdxIntersectionSorterKey(double md, bool entering, size_t cellIdx ): measuredDepth(md), hexIndex(cellIdx), isEnteringCell(entering){}
double measuredDepth;
bool isEnteringCell; // As opposed to leaving.
bool isLeavingCell() const { return !isEnteringCell;}
size_t hexIndex;
bool operator < (const RigMDEnterLeaveCellIdxIntersectionSorterKey& other) const
{
if (RigHexIntersector::isEqualDepth(measuredDepth, other.measuredDepth))
{
if (isEnteringCell == other.isEnteringCell)
{
if (hexIndex == other.hexIndex)
{
// Completely equal: intersections at cell edges or corners or edges of the face triangles
return false;
}
return (hexIndex < other.hexIndex);
}
return isLeavingCell(); // Sort Leaving cell before (less than) entering cell
}
// The depths are not equal
return (measuredDepth < other.measuredDepth);
}
static bool isProperPair(const RigMDEnterLeaveCellIdxIntersectionSorterKey& key1, const RigMDEnterLeaveCellIdxIntersectionSorterKey& key2 )
{
return ( key1.hexIndex == key2.hexIndex
&& key1.isEnteringCell && key2.isLeavingCell()
&& !RigHexIntersector::isEqualDepth(key1.measuredDepth, key2.measuredDepth));
}
};

View File

@ -0,0 +1,37 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) Statoil ASA
// Copyright (C) Ceetron Solutions AS
//
// 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 "RigWellLogExtractor.h"
#include "RigWellPath.h"
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigWellLogExtractor::RigWellLogExtractor(const RigWellPath* wellpath) : m_wellPath(wellpath)
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigWellLogExtractor::~RigWellLogExtractor()
{
}

View File

@ -0,0 +1,59 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) Statoil ASA
// Copyright (C) Ceetron Solutions AS
//
// 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 "cvfBase.h"
#include "cvfObject.h"
#include "cvfMath.h"
#include "cvfVector3.h"
#include <vector>
#include "cvfStructGrid.h"
class RigWellPath;
//==================================================================================================
///
//==================================================================================================
class RigWellLogExtractor : public cvf::Object
{
public:
RigWellLogExtractor(const RigWellPath* wellpath);
virtual ~RigWellLogExtractor();
const std::vector<double>& measuredDepth() { return m_measuredDepth; }
const std::vector<double>& trueVerticalDepth() {return m_trueVerticalDepth;}
const RigWellPath* wellPathData() { return m_wellPath.p();}
protected:
std::vector<double> m_measuredDepth;
std::vector<double> m_trueVerticalDepth;
std::vector<cvf::Vec3d> m_intersections;
std::vector<size_t> m_intersectedCells;
std::vector<cvf::StructGridInterface::FaceType>
m_intersectedCellFaces;
cvf::cref<RigWellPath> m_wellPath;
};