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

#1514, #1506 Fix several wellpath stimplan intersection bugs.

Browse Source 
First numbers to file achieved, but they are wrong and too few
This commit is contained in:
Jacob Støren 2017-05-26 17:34:31 +02:00
parent cedd8c42d0
commit e4abc737cd
7 changed files with 280 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
ApplicationCode/Commands/RicExportSelectedWellPathFractureWellCompletionFeature.cpp
View File

@ -52,21 +52,9 @@ CAF_CMD_SOURCE_INIT(RicExportSelectedWellPathFractureWellCompletionFeature, "Ric
//--------------------------------------------------------------------------------------------------
void RicExportSelectedWellPathFractureWellCompletionFeature::onActionTriggered(bool isChecked)
{
std::vector<RimWellPath*> selection;
caf::SelectionManager::instance()->objectsByType(&selection);
std::vector<RimFracture*> fractures;
for (RimWellPath* well : selection)
{
std::vector<RimFracture*> fracListForWell;
well->descendantsIncludingThisOfType(fracListForWell);
for (RimFracture* fracture : fracListForWell)
{
fractures.push_back(fracture);
}
}
if (!selection.size()) return;
RimFractureExportSettings exportSettings;
@ -74,20 +62,15 @@ void RicExportSelectedWellPathFractureWellCompletionFeature::onActionTriggered(b
QString projectFolder = app->currentProjectPath();
RimView* view = app->activeReservoirView();
caf::PdmObjectHandle* objHandle = dynamic_cast<caf::PdmObjectHandle*>(view);
if (!objHandle) return;
RimEclipseCase* caseToApply;
objHandle->firstAncestorOrThisOfType(caseToApply);
exportSettings.caseToApply = caseToApply;
view->firstAncestorOrThisOfType(caseToApply);
if (!caseToApply) return;
exportSettings.caseToApply = caseToApply;
if (projectFolder.isEmpty())
{
RimView* activeView = RiaApplication::instance()->activeReservoirView();
if (!activeView) return;
RimEclipseView * activeRiv = dynamic_cast<RimEclipseView*>(activeView);
if (!activeRiv) return;
projectFolder = activeRiv->eclipseCase()->locationOnDisc();
projectFolder = caseToApply->locationOnDisc();
}
QString defaultDir = RiaApplication::instance()->lastUsedDialogDirectoryWithFallback("FRACTURE_EXPORT_DIR", projectFolder);
@ -100,17 +83,10 @@ void RicExportSelectedWellPathFractureWellCompletionFeature::onActionTriggered(b
{
RiaApplication::instance()->setLastUsedDialogDirectory("FRACTURE_EXPORT_DIR", QFileInfo(exportSettings.fileName).absolutePath());
bool isOk = RifFractureExportTools::exportFracturesToEclipseDataInputFile(exportSettings.fileName, fractures, exportSettings.caseToApply);
if (!isOk)
{
QMessageBox::critical(NULL, "File export", "Failed to exported current result to " + exportSettings.fileName);
RifFractureExportTools::exportWellPathFracturesToEclipseDataInputFile(exportSettings.fileName, selection[0], exportSettings.caseToApply);
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

2
ApplicationCode/FileInterface/RifFractureExportTools.cpp
View File

@ -338,7 +338,7 @@ void RifFractureExportTools::exportWellPathFracturesToEclipseDataInputFile(const
double totalWellTrans = 0.5 * intersection.endpointCount * radialTrans + linearTrans;
transCondenser.addNeighborTransmissibility( { true, RigTransmissibilityCondenser::CellAddress::WELL, stpWellCellIdx},
transCondenser.addNeighborTransmissibility( { true, RigTransmissibilityCondenser::CellAddress::WELL, 1},
{ false, RigTransmissibilityCondenser::CellAddress::STIMPLAN, stpWellCellIdx },
totalWellTrans);
}

42
ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp
View File

@ -341,25 +341,39 @@ void fillInterpolatedSubjectZ(ClipperLib::IntPoint& e1bot,
ClipperLib::IntPoint& e2top,
ClipperLib::IntPoint& pt)
{
int e2XRange = (e2top.X - e2bot.X);
int e2YRange = (e2top.Y - e2bot.Y);
ClipperLib::IntPoint ePLbot;
ClipperLib::IntPoint ePLtop;
double e2Length = sqrt(e2XRange*e2XRange + e2YRange*e2YRange);
if (e2Length <= 1)
if (e1top.Z == std::numeric_limits<int>::max())
{
pt.Z = e2bot.Z;
ePLtop = e2top;
ePLbot = e2bot;
}
else
{
ePLtop = e1top;
ePLbot = e1bot;
}
double ePLXRange = (ePLtop.X - ePLbot.X);
double ePLYRange = (ePLtop.Y - ePLbot.Y);
double ePLLength = sqrt(ePLXRange*ePLXRange + ePLYRange*ePLYRange);
if (ePLLength <= 1)
{
pt.Z = ePLbot.Z;
return;
}
int e2BotPtXRange = pt.X - e2bot.X;
int e2BotPtYRange = pt.Y - e2bot.Y;
double ePLBotPtXRange = pt.X - ePLbot.X;
double ePLBotPtYRange = pt.Y - ePLbot.Y;
double e2BotPtLength = sqrt(e2BotPtXRange*e2BotPtXRange + e2BotPtYRange*e2BotPtYRange);
double ePLBotPtLength = sqrt(ePLBotPtXRange*ePLBotPtXRange + ePLBotPtYRange*ePLBotPtYRange);
double fraction = e2BotPtLength/e2Length;
double fraction = ePLBotPtLength/ePLLength;
pt.Z = std::nearbyint( e2bot.Z + fraction*(e2top.Z - e2bot.Z) );
pt.Z = std::nearbyint( ePLbot.Z + fraction*(ePLtop.Z - ePLbot.Z) );
}
//--------------------------------------------------------------------------------------------------
@ -375,7 +389,7 @@ void fillUndefinedZ(ClipperLib::IntPoint& e1bot,
}
//--------------------------------------------------------------------------------------------------
///
/// Assumes x.y plane polygon. Polyline might have a Z, and the returned Z is the polyline Z, interpolated if it is clipped.
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolylineByPolygon(const std::vector<cvf::Vec3d>& polyLine,
const std::vector<cvf::Vec3d>& polygon,
@ -393,7 +407,9 @@ std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolylineByPolygo
ClipperLib::Path polygonPath;
for (const cvf::Vec3d& v : polygon)
{
polygonPath.push_back(toClipperPoint(v));
ClipperLib::IntPoint intp = toClipperPoint(v);
intp.Z = std::numeric_limits<int>::max();
polygonPath.push_back(intp);
}
ClipperLib::Clipper clpr;

3
ApplicationCode/ReservoirDataModel/RigTransmissibilityCondenser.cpp
View File

@ -27,9 +27,10 @@
//--------------------------------------------------------------------------------------------------
void RigTransmissibilityCondenser::addNeighborTransmissibility(CellAddress cell1, CellAddress cell2, double transmissibility)
{
if (transmissibility < 1e-9) return;
m_condensedTransmissibilities.clear();
m_externalCellAddrSet.clear();
if ( cell1 < cell2 )
m_neighborTransmissibilities[cell1][cell2] = transmissibility;
else

125
ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.cpp
View File

@ -10,31 +10,55 @@
#include <cmath>
//--------------------------------------------------------------------------------------------------
/// Todo: Use only the perforated parts of the well path
//--------------------------------------------------------------------------------------------------
RigWellPathStimplanIntersector::RigWellPathStimplanIntersector(const RigWellPath* wellpathGeom, const RimFracture * rimFracture)
{
std::vector<cvf::Vec3d> wellPathPoints = wellpathGeom->m_wellPathPoints;
cvf::Mat4f fractureXf = rimFracture->transformMatrix();
double wellRadius = rimFracture->wellRadius();
std::vector<cvf::Vec3f> fracturePolygonf ;
std::vector<std::vector<cvf::Vec3d> > stpCellPolygons;
{
auto stimPlanFractureTemplate = dynamic_cast<RimStimPlanFractureTemplate*> (rimFracture->attachedFractureDefinition());
CVF_ASSERT(stimPlanFractureTemplate);
std::vector<cvf::Vec3d> wellPathPoints = wellpathGeom->m_wellPathPoints;
cvf::Mat4d toFractureXf = cvf::Mat4d (rimFracture->transformMatrix().getInverted());
double wellRadius = rimFracture->wellRadius();
std::vector<cvf::Vec3d> fracturePolygon;
fracturePolygonf = stimPlanFractureTemplate->fracturePolygon(rimFracture->fractureUnit());
{
std::vector<cvf::Vec3f> fracturePolygonf = stimPlanFractureTemplate->fracturePolygon(rimFracture->fractureUnit());
for ( auto fpv: fracturePolygonf ) fracturePolygon.push_back(cvf::Vec3d(fpv));
const std::vector<RigStimPlanFracTemplateCell>& stpCells = stimPlanFractureTemplate->getStimPlanCells();
for ( const auto& stpCell: stpCells ) stpCellPolygons.push_back(stpCell.getPolygon());
}
}
calculate(fractureXf, fracturePolygonf, wellPathPoints, wellRadius, stpCellPolygons, m_stimPlanCellIdxToIntersectionInfoMap);
}
//--------------------------------------------------------------------------------------------------
/// Todo: Use only the perforated parts of the well path
//--------------------------------------------------------------------------------------------------
void RigWellPathStimplanIntersector::calculate(const cvf::Mat4f &fractureXf,
const std::vector<cvf::Vec3f>& fracturePolygonf,
const std::vector<cvf::Vec3d>& wellPathPointsOrg,
double wellRadius,
const std::vector<std::vector<cvf::Vec3d> >& stpCellPolygons,
std::map<size_t, WellCellIntersection>& m_stimPlanCellIdxToIntersectionInfoMap)
{
cvf::Mat4d toFractureXf = cvf::Mat4d(fractureXf.getInverted());
std::vector<cvf::Vec3d> fracturePolygon;
for ( auto fpv: fracturePolygonf ) fracturePolygon.push_back(cvf::Vec3d(fpv));
// Convert well path to fracture template system
for ( auto & wellPPoint : wellPathPoints ) wellPPoint.transformPoint(toFractureXf);
std::vector<cvf::Vec3d> fractureRelativeWellPathPoints;
for ( auto & wellPPoint : wellPathPointsOrg ) fractureRelativeWellPathPoints.push_back(wellPPoint.getTransformedPoint( toFractureXf));
// Clip well path to fracture domain
std::vector<std::vector<cvf::Vec3d> > wellPathPartsWithinFracture =
RigCellGeometryTools::clipPolylineByPolygon(wellPathPoints, fracturePolygon, RigCellGeometryTools::INTERPOLATE_LINE_Z);
RigCellGeometryTools::clipPolylineByPolygon(fractureRelativeWellPathPoints, fracturePolygon, RigCellGeometryTools::INTERPOLATE_LINE_Z);
// Remove the part of the well path that is more than well radius away from the fracture plane
@ -45,22 +69,57 @@ RigWellPathStimplanIntersector::RigWellPathStimplanIntersector(const RigWellPath
std::vector< cvf::Vec3d > currentIntersectingWpPart;
for ( size_t vxIdx = 0; vxIdx < part.size() -1; ++vxIdx )
{
double thisZ = fabs(wellPathPoints[vxIdx].z());
double nextZ = fabs(wellPathPoints[vxIdx + 1].z());
double thisAbsZ = fabs(part[vxIdx].z());
double nextAbsZ = fabs(part[vxIdx + 1].z());
double thisZ = part[vxIdx].z();
double nextZ = part[vxIdx + 1].z();
if ( thisZ >= wellRadius && nextZ >= wellRadius ) continue;
if ( thisZ < wellRadius && nextZ < wellRadius )
if ( thisAbsZ >= wellRadius && nextAbsZ >= wellRadius )
{
currentIntersectingWpPart.push_back(wellPathPoints[vxIdx]);
if ( (thisZ >= 0 && nextZ >= 0)
|| (thisZ <= 0 && nextZ <= 0 ) )
{
continue; // Outside
}
else // In and out
{
{
double wellRadiusDistFromPlane = thisZ > 0 ? wellRadius: -wellRadius;
double fraction = (wellRadiusDistFromPlane - thisZ)/ (nextZ - thisZ);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx+1] - part[vxIdx]);
currentIntersectingWpPart.push_back(intersectPoint);
}
{
double wellRadiusDistFromPlane = nextZ > 0 ? wellRadius: -wellRadius;
double fraction = (wellRadiusDistFromPlane - thisZ)/ (nextZ - thisZ);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx+1] - part[vxIdx]);
currentIntersectingWpPart.push_back(intersectPoint);
intersectingWellPathParts.push_back(currentIntersectingWpPart);
currentIntersectingWpPart.clear();
}
continue;
}
}
if ( thisAbsZ < wellRadius && nextAbsZ < wellRadius ) // Inside
{
currentIntersectingWpPart.push_back(part[vxIdx]);
continue;
}
if ( thisZ < wellRadius && nextZ >= wellRadius )
if ( thisAbsZ < wellRadius && nextAbsZ >= wellRadius ) // Going out
{
currentIntersectingWpPart.push_back(wellPathPoints[vxIdx]);
double fraction = (wellRadius - thisZ)/ (nextZ - thisZ);
cvf::Vec3d intersectPoint = wellPathPoints[vxIdx] + fraction * (wellPathPoints[vxIdx+1] - wellPathPoints[vxIdx]);
currentIntersectingWpPart.push_back(part[vxIdx]);
double wellRadiusDistFromPlane = nextZ > 0 ? wellRadius: -wellRadius;
double fraction = (wellRadiusDistFromPlane - thisZ)/ (nextZ - thisZ);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx+1] - part[vxIdx]);
currentIntersectingWpPart.push_back(intersectPoint);
intersectingWellPathParts.push_back(currentIntersectingWpPart);
@ -68,13 +127,24 @@ RigWellPathStimplanIntersector::RigWellPathStimplanIntersector(const RigWellPath
continue;
}
if ( thisZ >= wellRadius && nextZ < wellRadius )
if ( thisAbsZ >= wellRadius && nextAbsZ < wellRadius ) // Going in
{
double fraction = (wellRadius - thisZ)/ (nextZ - thisZ);
cvf::Vec3d intersectPoint = wellPathPoints[vxIdx] + fraction * (wellPathPoints[vxIdx+1] - wellPathPoints[vxIdx]);
double wellRadiusDistFromPlane = thisZ > 0 ? wellRadius: -wellRadius;
double fraction = (wellRadiusDistFromPlane - thisZ)/ (nextZ - thisZ);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx+1] - part[vxIdx]);
currentIntersectingWpPart.push_back(intersectPoint);
continue;
}
}
// Add last point if it is within the radius
if (part.size() > 1 && fabs(part.back().z()) < wellRadius)
{
currentIntersectingWpPart.push_back(part.back());
}
if ( currentIntersectingWpPart.size() )
@ -85,12 +155,9 @@ RigWellPathStimplanIntersector::RigWellPathStimplanIntersector(const RigWellPath
// Find the StimPlan cells touched by the intersecting well path parts
const std::vector<RigStimPlanFracTemplateCell>& stpCells = stimPlanFractureTemplate->getStimPlanCells();
for ( size_t cIdx = 0; cIdx < stpCells.size(); ++ cIdx )
for ( size_t cIdx = 0; cIdx < stpCellPolygons.size(); ++ cIdx )
{
std::vector<cvf::Vec3d> cellPolygon = stpCells[cIdx].getPolygon();
const std::vector<cvf::Vec3d>& cellPolygon = stpCellPolygons[cIdx];
for ( const auto& wellpathPart :intersectingWellPathParts )
{
std::vector<std::vector<cvf::Vec3d> > wellPathPartsInPolygon =

26
ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.h
View File

@ -18,9 +18,14 @@
#pragma once
#include <map>
#include <vector>
#include "cvfBase.h"
#include "cvfMatrix4.h"
class RigWellPath;
class RimFracture;
class RigWellPathStimplanIntersectorTester;
class RigWellPathStimplanIntersector
{
@ -38,8 +43,29 @@ public:
const std::map<size_t, WellCellIntersection >& intersections() { return m_stimPlanCellIdxToIntersectionInfoMap; }
private:
friend class RigWellPathStimplanIntersectorTester;
static void calculate(const cvf::Mat4f& fractureXf,
const std::vector<cvf::Vec3f>& fracturePolygon,
const std::vector<cvf::Vec3d>& wellPathPoints,
double wellRadius,
const std::vector<std::vector<cvf::Vec3d> >& stpCellPolygons,
std::map<size_t, WellCellIntersection>& stimPlanCellIdxToIntersectionInfoMap);
std::map<size_t, WellCellIntersection > m_stimPlanCellIdxToIntersectionInfoMap;
};
class RigWellPathStimplanIntersectorTester
{
public:
static void testCalculate(const cvf::Mat4f& fractureXf,
const std::vector<cvf::Vec3f>& fracturePolygon,
const std::vector<cvf::Vec3d>& wellPathPoints,
double wellRadius,
const std::vector<std::vector<cvf::Vec3d> >& stpCellPolygons,
std::map<size_t, RigWellPathStimplanIntersector::WellCellIntersection>& stimPlanCellIdxToIntersectionInfoMap)
{
RigWellPathStimplanIntersector::calculate(fractureXf, fracturePolygon, wellPathPoints, wellRadius, stpCellPolygons, stimPlanCellIdxToIntersectionInfoMap);
}
};

117
ApplicationCode/UnitTests/RigCellGeometryTools-Test.cpp
View File

@ -302,3 +302,120 @@ TEST(RigCellGeometryTools, polylinePolygonIntersectionTest)
}
}
#include "RigWellPathStimplanIntersector.h"
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST(RigWellPathStimplanIntersector, intersection)
{
{
cvf::Mat4f fractureXf = cvf::Mat4f::IDENTITY;
fractureXf.setTranslation({ 50.0f, 0.0f, 0.0f });
std::vector<cvf::Vec3f> fracturePolygon ={ {0.0f, 0.0f, 0.0f}, {5.0f, 10.0f, 0.0f}, {10.0f, 0.0f, 0.0f} };
std::vector<cvf::Vec3d> wellPathPoints ={ {50.0f-4.0f, 6.0f, 10.0f}, {50.0f+6.0f, 6.0f, 0.0f}, {50.0f+10.0f, 10.0f, -100.0f} };
double wellRadius = 1.5;
std::vector<std::vector<cvf::Vec3d> > stpCellPolygons =
{
{ { 0.0f, 0.0f, 0.0f }, { 0.0f, 5.0f, 0.0f }, { 5.0f, 5.0f, 0.0f }, { 5.0f, 0.0f, 0.0f } },
{ { 0.5f, 0.0f, 0.0f }, { 0.5f, 5.0f, 0.0f }, {10.0f, 5.0f, 0.0f }, {10.0f, 0.0f, 0.0f } },
{ { 0.0f, 5.0f, 0.0f }, { 0.0f,10.0f, 0.0f }, { 5.0f,10.0f, 0.0f }, { 5.0f, 5.0f, 0.0f } },
{ { 5.0f, 5.0f, 0.0f }, { 5.0f,10.0f, 0.0f }, {10.0f,10.0f, 0.0f }, {10.0f, 5.0f, 0.0f } },
};
std::map<size_t, RigWellPathStimplanIntersector::WellCellIntersection> stimPlanCellIdxToIntersectionInfoMap;
RigWellPathStimplanIntersectorTester::testCalculate(fractureXf,
fracturePolygon,
wellPathPoints,
wellRadius,
stpCellPolygons,
stimPlanCellIdxToIntersectionInfoMap);
EXPECT_EQ(2, stimPlanCellIdxToIntersectionInfoMap.size());
auto it = stimPlanCellIdxToIntersectionInfoMap.begin();
EXPECT_EQ(2, it->first);
EXPECT_EQ(1, it->second.endpointCount);
++it;
EXPECT_EQ(3, it->first);
EXPECT_EQ(1, it->second.endpointCount);
}
{
cvf::Mat4f fractureXf = cvf::Mat4f::IDENTITY;
std::vector<cvf::Vec3f> fracturePolygon ={ {0.0f, 0.0f, 0.0f}, {5.0f, 10.0f, 0.0f}, {10.0f, 0.0f, 0.0f} };
double wellRadius = 1.5;
std::vector<std::vector<cvf::Vec3d> > stpCellPolygons =
{
{ { 0.0f, 0.0f, 0.0f }, { 0.0f, 5.0f, 0.0f }, { 5.0f, 5.0f, 0.0f }, { 5.0f, 0.0f, 0.0f } },
{ { 5.0f, 0.0f, 0.0f }, { 5.0f, 5.0f, 0.0f }, {10.0f, 5.0f, 0.0f }, {10.0f, 0.0f, 0.0f } },
{ { 0.0f, 5.0f, 0.0f }, { 0.0f,10.0f, 0.0f }, { 5.0f,10.0f, 0.0f }, { 5.0f, 5.0f, 0.0f } },
{ { 5.0f, 5.0f, 0.0f }, { 5.0f,10.0f, 0.0f }, {10.0f,10.0f, 0.0f }, {10.0f, 5.0f, 0.0f } },
};
{
std::map<size_t, RigWellPathStimplanIntersector::WellCellIntersection> stimPlanCellIdxToIntersectionInfoMap;
std::vector<cvf::Vec3d> wellPathPoints ={ {1.0f, 0.5f, 10.0f}, {1.0f, 1.5f, -10.0f} };
RigWellPathStimplanIntersectorTester::testCalculate(fractureXf,
fracturePolygon,
wellPathPoints,
wellRadius,
stpCellPolygons,
stimPlanCellIdxToIntersectionInfoMap);
EXPECT_EQ(1, stimPlanCellIdxToIntersectionInfoMap.size());
auto it = stimPlanCellIdxToIntersectionInfoMap.begin();
EXPECT_EQ(0, it->first);
EXPECT_EQ(2, it->second.endpointCount);
}
{
std::map<size_t, RigWellPathStimplanIntersector::WellCellIntersection> stimPlanCellIdxToIntersectionInfoMap;
std::vector<cvf::Vec3d> wellPathPoints ={ {1.0f, 0.5f, 10.0f}, {1.0f, 1.0f, 0.5f} };
RigWellPathStimplanIntersectorTester::testCalculate(fractureXf,
fracturePolygon,
wellPathPoints,
wellRadius,
stpCellPolygons,
stimPlanCellIdxToIntersectionInfoMap);
EXPECT_EQ(1, stimPlanCellIdxToIntersectionInfoMap.size());
auto it = stimPlanCellIdxToIntersectionInfoMap.begin();
EXPECT_EQ(0, it->first);
EXPECT_EQ(2, it->second.endpointCount);
}
{
std::map<size_t, RigWellPathStimplanIntersector::WellCellIntersection> stimPlanCellIdxToIntersectionInfoMap;
std::vector<cvf::Vec3d> wellPathPoints ={ {1.0f, 0.5f, 10.0f}, {1.0f, 1.0f, 0.5f}, {1.0f, 1.5f, -0.5f}, {1.0f, 2.0f, -10.0f}};
RigWellPathStimplanIntersectorTester::testCalculate(fractureXf,
fracturePolygon,
wellPathPoints,
wellRadius,
stpCellPolygons,
stimPlanCellIdxToIntersectionInfoMap);
EXPECT_EQ(1, stimPlanCellIdxToIntersectionInfoMap.size());
auto it = stimPlanCellIdxToIntersectionInfoMap.begin();
EXPECT_EQ(0, it->first);
EXPECT_EQ(2, it->second.endpointCount);
}
}
}