OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-10 15:36:09 -06:00
Code Issues Packages Projects Releases Wiki Activity

#1514 Add Wellpath StimPlanGrid intersector. Fix errors in polyline polygon intersection. Add tests for Polyline polygon intersection

Browse Source
This commit is contained in:
Jacob Støren 2017-05-22 18:11:03 +02:00
parent 50c182aa26
commit 6cbbed15bf
12 changed files with 298 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
ApplicationCode/FileInterface/RifFractureExportTools.cpp
View File

@ -47,6 +47,7 @@
#include <QTextStream>
#include "RigStimPlanUpscalingCalc.h"
#include "RigTransmissibilityCondenser.h"
#include "RigWellPathStimplanIntersector.h"
@ -309,7 +310,14 @@ void RifFractureExportTools::exportWellPathFracturesToEclipseDataInputFile(const
// Use LinT + 1/2 radialT on each end if endpoints are inside cell
// Add to condenser
//Adding well cell to condenser
// Adding well cells to condenser
//RigWellPathStimplanIntersector wellFractureIntersector(wellPath->wellPathGeometry(), fracture);
// Todo: ..
std::pair<size_t, size_t> centerCellIJ = fracTemplateStimPlan->getStimPlanCellAtWellCenter();
size_t wellCellIndex = fracTemplateStimPlan->getGlobalIndexFromIJ(centerCellIJ.first, centerCellIJ.second);
const RigStimPlanFracTemplateCell stimPlanWellCell = fracTemplateStimPlan->stimPlanCellFromIndex(wellCellIndex);
@ -325,12 +333,8 @@ void RifFractureExportTools::exportWellPathFracturesToEclipseDataInputFile(const
{false, RigTransmissibilityCondenser::CellAddress::STIMPLAN, wellCellIndex},
wellTrans);
// std::map <size_t, std::map<RimFracture*, double>> transmissibilitiesMap; //eclipseCellIndex (size_t), fracture (RimFracture) and trans value
// Add eclipse cell transmissibilities to a map <eclipsecell index, map<fractureptr, transmissibility> >
std::set<RigTransmissibilityCondenser::CellAddress> externalCells = condenser.externalCells();
for (RigTransmissibilityCondenser::CellAddress externalCell : externalCells)
{
if (externalCell.m_cellIndexSpace == RigTransmissibilityCondenser::CellAddress::ECLIPSE)

4
ApplicationCode/ProjectDataModel/RimFracture.cpp
View File

@ -245,7 +245,7 @@ void RimFracture::computeGeometry()
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RimFracture::anchorPosition()
cvf::Vec3d RimFracture::anchorPosition() const
{
return m_anchorPosition();
}
@ -253,7 +253,7 @@ cvf::Vec3d RimFracture::anchorPosition()
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Mat4f RimFracture::transformMatrix()
cvf::Mat4f RimFracture::transformMatrix() const
{
cvf::Vec3d center = anchorPosition();

4
ApplicationCode/ProjectDataModel/RimFracture.h
View File

@ -68,10 +68,10 @@ public:
caf::PdmField< caf::AppEnum< RimDefines::UnitSystem > > fractureUnit;
cvf::Vec3d anchorPosition();
cvf::Vec3d anchorPosition() const ;
void setAnchorPosition(const cvf::Vec3d& pos);
cvf::Mat4f transformMatrix();
cvf::Mat4f transformMatrix() const;
const RigFracture* attachedRigFracture() const;

8
ApplicationCode/ProjectDataModel/RimWellPath.cpp
View File

@ -170,6 +170,14 @@ RigWellPath* RimWellPath::wellPathGeometry()
return m_wellPath.p();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const RigWellPath* RimWellPath::wellPathGeometry() const
{
return m_wellPath.p();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

2
ApplicationCode/ProjectDataModel/RimWellPath.h
View File

@ -72,6 +72,8 @@ public:
caf::PdmChildField<RimWellLogFile*> m_wellLogFile;
RigWellPath* wellPathGeometry();
const RigWellPath* wellPathGeometry() const;
caf::PdmChildField<RimWellPathFractureCollection*> fractureCollection;

2
ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake
View File

@ -56,6 +56,7 @@ ${CEE_CURRENT_LIST_DIR}RigSimulationWellCoordsAndMD.h
${CEE_CURRENT_LIST_DIR}RigFracture.h
${CEE_CURRENT_LIST_DIR}RigFractureTransCalc.h
${CEE_CURRENT_LIST_DIR}RigFractureTransmissibilityEquations.h
${CEE_CURRENT_LIST_DIR}RigWellPathStimplanIntersector.h
${CEE_CURRENT_LIST_DIR}RigTesselatorTools.h
${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.h
${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.h
@ -114,6 +115,7 @@ ${CEE_CURRENT_LIST_DIR}RigSimulationWellCoordsAndMD.cpp
${CEE_CURRENT_LIST_DIR}RigFracture.cpp
${CEE_CURRENT_LIST_DIR}RigFractureTransCalc.cpp
${CEE_CURRENT_LIST_DIR}RigFractureTransmissibilityEquations.cpp
${CEE_CURRENT_LIST_DIR}RigWellPathStimplanIntersector.cpp
${CEE_CURRENT_LIST_DIR}RigTesselatorTools.cpp
${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.cpp
${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.cpp

33
ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp
View File

@ -276,7 +276,18 @@ ClipperLib::IntPoint toClipperPoint(const cvf::Vec3d& cvfPoint)
cvf::Vec3d fromClipperPoint(const ClipperLib::IntPoint& clipPoint)
{
return cvf::Vec3d (clipPoint.X, clipPoint.Y, clipPoint.Z ) /clipperConversionFactor;
double zDValue;
if (clipPoint.Z == std::numeric_limits<int>::max())
{
zDValue = HUGE_VAL;
}
else
{
zDValue = clipPoint.Z;
}
return cvf::Vec3d (clipPoint.X, clipPoint.Y, zDValue ) /clipperConversionFactor;
}
//--------------------------------------------------------------------------------------------------
@ -348,7 +359,7 @@ void fillInterpolatedSubjectZ(ClipperLib::IntPoint& e1bot,
double fraction = e2BotPtLength/e2Length;
pt.Z = e2bot.Z + fraction*(e2top.Z - e2bot.Z);
pt.Z = std::nearbyint( e2bot.Z + fraction*(e2top.Z - e2bot.Z) );
}
//--------------------------------------------------------------------------------------------------
@ -360,27 +371,27 @@ void fillUndefinedZ(ClipperLib::IntPoint& e1bot,
ClipperLib::IntPoint& e2top,
ClipperLib::IntPoint& pt)
{
pt.Z = HUGE_VAL;
pt.Z = std::numeric_limits<int>::max();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolylineByPolygon(std::vector<cvf::Vec3d> polyLine,
std::vector<cvf::Vec3d> polygon,
std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolylineByPolygon(const std::vector<cvf::Vec3d>& polyLine,
const std::vector<cvf::Vec3d>& polygon,
ZInterpolationType interpolType)
{
int polygonScaleFactor = 10000; //For transform to clipper int
//Convert to int for clipper library and store as clipper "path"
ClipperLib::Path polyLinePath;
for (cvf::Vec3d& v : polyLine)
for (const cvf::Vec3d& v : polyLine)
{
polyLinePath.push_back(toClipperPoint(v));
}
ClipperLib::Path polygonPath;
for (cvf::Vec3d& v : polygon)
for (const cvf::Vec3d& v : polygon)
{
polygonPath.push_back(toClipperPoint(v));
}
@ -398,12 +409,16 @@ std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolylineByPolygo
clpr.ZFillFunction(&fillUndefinedZ);
}
ClipperLib::Paths solution;
ClipperLib::PolyTree solution;
clpr.Execute(ClipperLib::ctIntersection, solution, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);
// We only expect open paths from this method (unless the polyline is self intersecting, a condition we do not handle)
ClipperLib::Paths solutionPaths;
ClipperLib::OpenPathsFromPolyTree(solution, solutionPaths);
//Convert back to std::vector<std::vector<cvf::Vec3d> >
std::vector<std::vector<cvf::Vec3d> > clippedPolyline;
for (ClipperLib::Path pathInSol : solution)
for (ClipperLib::Path pathInSol : solutionPaths)
{
std::vector<cvf::Vec3d> clippedPolygon;
for (ClipperLib::IntPoint IntPosition : pathInSol)

4
ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h
View File

@ -46,7 +46,9 @@ public:
static std::vector<std::vector<cvf::Vec3d> > intersectPolygons(std::vector<cvf::Vec3d> polygon1, std::vector<cvf::Vec3d> polygon2);
enum ZInterpolationType { INTERPOLATE_LINE_Z, USE_HUGEVAL, USE_ZERO};
static std::vector<std::vector<cvf::Vec3d> > clipPolylineByPolygon(std::vector<cvf::Vec3d> polyLine, std::vector<cvf::Vec3d> polygon, ZInterpolationType interpolType = USE_ZERO);
static std::vector<std::vector<cvf::Vec3d> > clipPolylineByPolygon(const std::vector<cvf::Vec3d>& polyLine,
const std::vector<cvf::Vec3d>& polygon,
ZInterpolationType interpolType = USE_ZERO);
static std::pair<cvf::Vec3d, cvf::Vec3d> getLineThroughBoundingBox(cvf::Vec3d lineDirection, cvf::BoundingBox polygonBBox, cvf::Vec3d pointOnLine);

119
ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.cpp Normal file
View File

@ -0,0 +1,119 @@
#include "RigWellPathStimplanIntersector.h"
#include "RigWellPath.h"
#include "RimFracture.h"
#include "RigCellGeometryTools.h"
#include "RimFractureTemplate.h"
#include "cvfBase.h"
#include "cvfMatrix4.h"
#include "RigStimPlanFracTemplateCell.h"
#include "RimStimPlanFractureTemplate.h"
//--------------------------------------------------------------------------------------------------
/// Todo: Use only the perforated parts of the well path
//--------------------------------------------------------------------------------------------------
RigWellPathStimplanIntersector::RigWellPathStimplanIntersector(const RigWellPath* wellpathGeom, const RimFracture * rimFracture)
{
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;
{
std::vector<cvf::Vec3f> fracturePolygonf = stimPlanFractureTemplate->fracturePolygon(rimFracture->fractureUnit());
for ( auto fpv: fracturePolygonf ) fracturePolygon.push_back(cvf::Vec3d(fpv));
}
// Convert well path to fracture template system
for ( auto & wellPPoint : wellPathPoints ) wellPPoint.transformPoint(toFractureXf);
// Clip well path to fracture domain
std::vector<std::vector<cvf::Vec3d> > wellPathPartsWithinFracture =
RigCellGeometryTools::clipPolylineByPolygon(wellPathPoints, fracturePolygon, RigCellGeometryTools::INTERPOLATE_LINE_Z);
// Remove the part of the well path that is more than well radius away from the fracture plane
std::vector< std::vector< cvf::Vec3d > > intersectingWellPathParts;
for ( const auto& part : wellPathPartsWithinFracture )
{
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());
if ( thisZ >= wellRadius && nextZ >= wellRadius ) continue;
if ( thisZ < wellRadius && nextZ < wellRadius )
{
currentIntersectingWpPart.push_back(wellPathPoints[vxIdx]);
continue;
}
if ( thisZ < wellRadius && nextZ >= wellRadius )
{
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(intersectPoint);
intersectingWellPathParts.push_back(currentIntersectingWpPart);
currentIntersectingWpPart.clear();
continue;
}
if ( thisZ >= wellRadius && nextZ < wellRadius )
{
double fraction = (wellRadius - thisZ)/ (nextZ - thisZ);
cvf::Vec3d intersectPoint = wellPathPoints[vxIdx] + fraction * (wellPathPoints[vxIdx+1] - wellPathPoints[vxIdx]);
currentIntersectingWpPart.push_back(intersectPoint);
continue;
}
}
if ( currentIntersectingWpPart.size() )
{
intersectingWellPathParts.push_back(currentIntersectingWpPart);
}
}
// 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 )
{
std::vector<cvf::Vec3d> cellPolygon = stpCells[cIdx].getPolygon();
for ( const auto& wellpathPart :intersectingWellPathParts )
{
std::vector<std::vector<cvf::Vec3d> > wellPathPartsInPolygon =
RigCellGeometryTools::clipPolylineByPolygon(wellpathPart,
cellPolygon,
RigCellGeometryTools::USE_HUGEVAL);
for ( const auto& wellPathPartInCell: wellPathPartsInPolygon )
{
if ( wellPathPartInCell.size() )
{
int endpointCount = 0;
if ( wellPathPartInCell.front().z() != HUGE_VAL ) ++endpointCount;
if ( wellPathPartInCell.back().z() != HUGE_VAL ) ++endpointCount;
cvf::Vec3d intersectionLength = (wellPathPartInCell.back() - wellPathPartInCell.front());
double xLengthInCell = fabs(intersectionLength.x());
double yLengthInCell = fabs(intersectionLength.y());
m_stimPlanCellIdxToIntersectionInfoMap[cIdx].endpointCount += endpointCount;
m_stimPlanCellIdxToIntersectionInfoMap[cIdx].hlength += xLengthInCell;
m_stimPlanCellIdxToIntersectionInfoMap[cIdx].vlength += yLengthInCell;
}
}
}
}
}

45
ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.h Normal file
View File

@ -0,0 +1,45 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2017- 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 <map>
class RigWellPath;
class RimFracture;
class RigWellPathStimplanIntersector
{
public:
struct WellCellIntersection
{
WellCellIntersection():hlength(0.0), vlength(0.0), endpointCount(0) {}
double hlength;
double vlength;
int endpointCount;
};
RigWellPathStimplanIntersector(const RigWellPath* wellpathGeom, const RimFracture * rimFracture);
const std::map<size_t, WellCellIntersection >& intersections() { return m_stimPlanCellIdxToIntersectionInfoMap; }
private:
std::map<size_t, WellCellIntersection > m_stimPlanCellIdxToIntersectionInfoMap;
};

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

@ -221,3 +221,84 @@ TEST(RigCellGeometryTools, lengthCalcTest)
length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, polygonExample);
EXPECT_DOUBLE_EQ(length, 2.5);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST(RigCellGeometryTools, polylinePolygonIntersectionTest)
{
std::vector<cvf::Vec3d> polygonExample;
polygonExample.push_back(cvf::Vec3d(0.00, 0.00, 0.0));
polygonExample.push_back(cvf::Vec3d(0.00, 2.50, 0.0));
polygonExample.push_back(cvf::Vec3d(1.50, 2.50, 0.0));
polygonExample.push_back(cvf::Vec3d(1.50, 0.00, 0.0));
std::vector<cvf::Vec3d> polyLine;
polyLine.push_back(cvf::Vec3d(-1.00, 0.00, 1.0));
polyLine.push_back(cvf::Vec3d(1.00, 2.0, 2.0));
polyLine.push_back(cvf::Vec3d(1.00, 3.00, 3.0));
{
std::vector< std::vector<cvf::Vec3d> > clippedLines = RigCellGeometryTools::clipPolylineByPolygon(polyLine,
polygonExample,
RigCellGeometryTools::INTERPOLATE_LINE_Z);
EXPECT_EQ(1, clippedLines.size());
EXPECT_EQ(3, clippedLines.front().size());
EXPECT_EQ(0.0, clippedLines.front()[0].x());
EXPECT_EQ(1.0, clippedLines.front()[0].y());
EXPECT_EQ(1.5, clippedLines.front()[0].z());
EXPECT_EQ(1.0, clippedLines.front()[1].x());
EXPECT_EQ(2.0, clippedLines.front()[1].y());
EXPECT_EQ(2.0, clippedLines.front()[1].z());
EXPECT_EQ(1.0, clippedLines.front()[2].x());
EXPECT_EQ(2.5, clippedLines.front()[2].y());
EXPECT_EQ(2.5, clippedLines.front()[2].z());
}
{
std::vector< std::vector<cvf::Vec3d> > clippedLines = RigCellGeometryTools::clipPolylineByPolygon(polyLine,
polygonExample,
RigCellGeometryTools::USE_HUGEVAL);
EXPECT_EQ(1, clippedLines.size());
EXPECT_EQ(3, clippedLines.front().size());
EXPECT_EQ(0.0, clippedLines.front()[0].x());
EXPECT_EQ(1.0, clippedLines.front()[0].y());
EXPECT_EQ(HUGE_VAL, clippedLines.front()[0].z());
EXPECT_EQ(1.0, clippedLines.front()[1].x());
EXPECT_EQ(2.0, clippedLines.front()[1].y());
EXPECT_EQ(2.0, clippedLines.front()[1].z());
EXPECT_EQ(1.0, clippedLines.front()[2].x());
EXPECT_EQ(2.5, clippedLines.front()[2].y());
EXPECT_EQ(HUGE_VAL, clippedLines.front()[2].z());
}
polyLine.push_back({-0.5, 1.5, 0.0});
{
std::vector< std::vector<cvf::Vec3d> > clippedLines = RigCellGeometryTools::clipPolylineByPolygon(polyLine,
polygonExample,
RigCellGeometryTools::USE_HUGEVAL);
EXPECT_EQ(2, clippedLines.size());
EXPECT_EQ(2, clippedLines.front().size());
EXPECT_EQ(3, clippedLines.back().size());
EXPECT_EQ(0.5, clippedLines.front()[0].x());
EXPECT_EQ(2.5, clippedLines.front()[0].y());
EXPECT_EQ(HUGE_VAL, clippedLines.front()[0].z());
EXPECT_EQ(0.0, clippedLines.front()[1].x());
EXPECT_EQ(2.0, clippedLines.front()[1].y());
EXPECT_EQ(HUGE_VAL, clippedLines.front()[1].z());
}
}

2
ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp
View File

@ -36,7 +36,7 @@ TEST(opm_flowdiagnostics_test, basic_construction)
std::cout.precision(16);
for ( double t : tof )
{
std::cout << t << '\n';
// std::cout << t << '\n';
}
}
catch ( const std::exception& e )