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

Whitespace and rename

Browse Source
This commit is contained in:
Magne Sjaastad 2018-08-13 12:54:21 +02:00
parent e195eed246
commit 45f0f6b485
1 changed files with 86 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

162
ApplicationCode/ReservoirDataModel/Completions/RigWellPathStimplanIntersector.cpp
View File

@ -1,17 +1,17 @@
/////////////////////////////////////////////////////////////////////////////////
//
// 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>
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
@ -34,34 +34,44 @@
#include <cmath>
//--------------------------------------------------------------------------------------------------
///
///
//--------------------------------------------------------------------------------------------------
RigWellPathStimplanIntersector::RigWellPathStimplanIntersector(const RigWellPath* wellpathGeom, const RimFracture* rimFracture)
RigWellPathStimplanIntersector::RigWellPathStimplanIntersector(const RigWellPath* wellPathGeom, const RimFracture* rimFracture)
{
std::vector<cvf::Vec3d> wellPathPoints = wellpathGeom->wellPathPointsIncludingInterpolatedIntersectionPoint(rimFracture->fractureMD());
cvf::Mat4d fractureXf = rimFracture->transformMatrix();
double wellRadius = rimFracture->wellRadius();
std::vector<std::vector<cvf::Vec3d> > stpCellPolygons;
std::vector<cvf::Vec3d> wellPathPoints =
wellPathGeom->wellPathPointsIncludingInterpolatedIntersectionPoint(rimFracture->fractureMD());
cvf::Mat4d fractureXf = rimFracture->transformMatrix();
double wellRadius = rimFracture->wellRadius();
std::vector<std::vector<cvf::Vec3d>> fractureGridCellPolygons;
{
RimFractureTemplate* fractureTemplate = rimFracture->fractureTemplate();
if(fractureTemplate)
if (fractureTemplate && fractureTemplate->fractureGrid())
{
const std::vector<RigFractureCell>& stpCells = fractureTemplate->fractureGrid()->fractureCells();
for ( const auto& stpCell: stpCells ) stpCellPolygons.push_back(stpCell.getPolygon());
for (const auto& stpCell : stpCells)
{
fractureGridCellPolygons.push_back(stpCell.getPolygon());
}
}
}
double perforationLength = rimFracture->perforationLength();
calculate(fractureXf, wellPathPoints, wellRadius, perforationLength, stpCellPolygons, m_stimPlanCellIdxToIntersectionInfoMap);
calculate(fractureXf,
wellPathPoints,
wellRadius,
perforationLength,
fractureGridCellPolygons,
m_stimPlanCellIdxToIntersectionInfoMap);
}
//--------------------------------------------------------------------------------------------------
///
///
//--------------------------------------------------------------------------------------------------
const std::map<size_t, RigWellPathStimplanIntersector::RigWellPathStimplanIntersector::WellCellIntersection >& RigWellPathStimplanIntersector::intersections() const
const std::map<size_t, RigWellPathStimplanIntersector::RigWellPathStimplanIntersector::WellCellIntersection>&
RigWellPathStimplanIntersector::intersections() const
{
return m_stimPlanCellIdxToIntersectionInfoMap;
}
@ -69,74 +79,77 @@ const std::map<size_t, RigWellPathStimplanIntersector::RigWellPathStimplanInters
//--------------------------------------------------------------------------------------------------
/// Todo: Use only the perforated parts of the well path
//--------------------------------------------------------------------------------------------------
void RigWellPathStimplanIntersector::calculate(const cvf::Mat4d &fractureXf,
const std::vector<cvf::Vec3d>& wellPathPointsOrg,
double wellRadius,
double perforationLength,
const std::vector<std::vector<cvf::Vec3d> >& stpCellPolygons,
std::map<size_t, WellCellIntersection>& m_stimPlanCellIdxToIntersectionInfoMap)
void RigWellPathStimplanIntersector::calculate(const cvf::Mat4d& fractureXf,
const std::vector<cvf::Vec3d>& wellPathPointsDomainCoords,
double wellRadius,
double perforationLength,
const std::vector<std::vector<cvf::Vec3d>>& fractureGridCellPolygons,
std::map<size_t, WellCellIntersection>& m_stimPlanCellIdxToIntersectionInfoMap)
{
cvf::Mat4d toFractureXf = fractureXf.getInverted();
std::vector<cvf::Vec3d> perforationLengthBoundingBoxPolygon;
double cicleRadius = perforationLength / 2;
int pointsInCirclePolygon = 20;
for (int i = 0; i < pointsInCirclePolygon; i++)
{
double x = cicleRadius * cvf::Math::cos(i * (2 * cvf::PI_D / pointsInCirclePolygon));
double y = cicleRadius * cvf::Math::sin(i * (2 * cvf::PI_D / pointsInCirclePolygon));
perforationLengthBoundingBoxPolygon.push_back(cvf::Vec3d(x, y, 0));
double cicleRadius = perforationLength / 2;
int pointsInCirclePolygon = 20;
for (int i = 0; i < pointsInCirclePolygon; i++)
{
double x = cicleRadius * cvf::Math::cos(i * (2 * cvf::PI_D / pointsInCirclePolygon));
double y = cicleRadius * cvf::Math::sin(i * (2 * cvf::PI_D / pointsInCirclePolygon));
perforationLengthBoundingBoxPolygon.push_back(cvf::Vec3d(x, y, 0));
}
}
// Convert well path to fracture template system
std::vector<cvf::Vec3d> fractureRelativeWellPathPoints;
for ( auto & wellPPoint : wellPathPointsOrg ) fractureRelativeWellPathPoints.push_back(wellPPoint.getTransformedPoint( toFractureXf));
for (const auto& wellPPoint : wellPathPointsDomainCoords)
{
fractureRelativeWellPathPoints.push_back(wellPPoint.getTransformedPoint(toFractureXf));
}
// Clip well path to fracture domain
// Clip well path to fracture domain
std::vector<std::vector<cvf::Vec3d> > wellPathPartsWithinFracture =
RigCellGeometryTools::clipPolylineByPolygon(fractureRelativeWellPathPoints,
perforationLengthBoundingBoxPolygon,
RigCellGeometryTools::INTERPOLATE_LINE_Z);
std::vector<std::vector<cvf::Vec3d>> wellPathPartsWithinFracture = RigCellGeometryTools::clipPolylineByPolygon(
fractureRelativeWellPathPoints, perforationLengthBoundingBoxPolygon, 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;
std::vector<std::vector<cvf::Vec3d>> intersectingWellPathParts;
for ( const auto& part : wellPathPartsWithinFracture )
for (const auto& part : wellPathPartsWithinFracture)
{
std::vector< cvf::Vec3d > currentIntersectingWpPart;
for ( size_t vxIdx = 0; vxIdx < part.size() -1; ++vxIdx )
std::vector<cvf::Vec3d> currentIntersectingWpPart;
for (size_t vxIdx = 0; vxIdx < part.size() - 1; ++vxIdx)
{
double thisAbsZ = fabs(part[vxIdx].z());
double nextAbsZ = fabs(part[vxIdx + 1].z());
double thisZ = part[vxIdx].z();
double nextZ = part[vxIdx + 1].z();
double thisZ = part[vxIdx].z();
double nextZ = part[vxIdx + 1].z();
if ( thisAbsZ >= wellRadius && nextAbsZ >= wellRadius )
if (thisAbsZ >= wellRadius && nextAbsZ >= wellRadius)
{
if ( (thisZ >= 0 && nextZ >= 0)
|| (thisZ <= 0 && nextZ <= 0 ) )
if ((thisZ >= 0 && nextZ >= 0) || (thisZ <= 0 && nextZ <= 0))
{
continue; // Outside
}
else // In and out
{
{
double wellRadiusDistFromPlane = thisZ > 0 ? wellRadius: -wellRadius;
double wellRadiusDistFromPlane = thisZ > 0 ? wellRadius : -wellRadius;
double fraction = (wellRadiusDistFromPlane - thisZ)/ (nextZ - thisZ);
double fraction = (wellRadiusDistFromPlane - thisZ) / (nextZ - thisZ);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx+1] - part[vxIdx]);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx + 1] - part[vxIdx]);
currentIntersectingWpPart.push_back(intersectPoint);
}
{
double wellRadiusDistFromPlane = nextZ > 0 ? wellRadius: -wellRadius;
double wellRadiusDistFromPlane = nextZ > 0 ? wellRadius : -wellRadius;
double fraction = (wellRadiusDistFromPlane - thisZ)/ (nextZ - thisZ);
double fraction = (wellRadiusDistFromPlane - thisZ) / (nextZ - thisZ);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx+1] - part[vxIdx]);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx + 1] - part[vxIdx]);
currentIntersectingWpPart.push_back(intersectPoint);
intersectingWellPathParts.push_back(currentIntersectingWpPart);
@ -145,21 +158,21 @@ void RigWellPathStimplanIntersector::calculate(const cvf::Mat4d &fractureXf,
continue;
}
}
if ( thisAbsZ < wellRadius && nextAbsZ < wellRadius ) // Inside
if (thisAbsZ < wellRadius && nextAbsZ < wellRadius) // Inside
{
currentIntersectingWpPart.push_back(part[vxIdx]);
continue;
}
if ( thisAbsZ < wellRadius && nextAbsZ >= wellRadius ) // Going out
if (thisAbsZ < wellRadius && nextAbsZ >= wellRadius) // Going out
{
currentIntersectingWpPart.push_back(part[vxIdx]);
double wellRadiusDistFromPlane = nextZ > 0 ? wellRadius: -wellRadius;
double wellRadiusDistFromPlane = nextZ > 0 ? wellRadius : -wellRadius;
double fraction = (wellRadiusDistFromPlane - thisZ)/ (nextZ - thisZ);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx+1] - part[vxIdx]);
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);
@ -167,54 +180,51 @@ void RigWellPathStimplanIntersector::calculate(const cvf::Mat4d &fractureXf,
continue;
}
if ( thisAbsZ >= wellRadius && nextAbsZ < wellRadius ) // Going in
if (thisAbsZ >= wellRadius && nextAbsZ < wellRadius) // Going in
{
double wellRadiusDistFromPlane = thisZ > 0 ? wellRadius: -wellRadius;
double wellRadiusDistFromPlane = thisZ > 0 ? wellRadius : -wellRadius;
double fraction = (wellRadiusDistFromPlane - thisZ)/ (nextZ - thisZ);
double fraction = (wellRadiusDistFromPlane - thisZ) / (nextZ - thisZ);
cvf::Vec3d intersectPoint = part[vxIdx] + fraction * (part[vxIdx+1] - part[vxIdx]);
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)
if (part.size() > 1 && fabs(part.back().z()) < wellRadius)
{
currentIntersectingWpPart.push_back(part.back());
}
if ( !currentIntersectingWpPart.empty() )
if (!currentIntersectingWpPart.empty())
{
intersectingWellPathParts.push_back(currentIntersectingWpPart);
}
}
// Find the StimPlan cells touched by the intersecting well path parts
// Find the StimPlan cells touched by the intersecting well path parts
for ( size_t cIdx = 0; cIdx < stpCellPolygons.size(); ++ cIdx )
for (size_t cIdx = 0; cIdx < fractureGridCellPolygons.size(); ++cIdx)
{
const std::vector<cvf::Vec3d>& cellPolygon = stpCellPolygons[cIdx];
for ( const auto& wellpathPart :intersectingWellPathParts )
const std::vector<cvf::Vec3d>& cellPolygon = fractureGridCellPolygons[cIdx];
for (const auto& wellpathPart : intersectingWellPathParts)
{
std::vector<std::vector<cvf::Vec3d> > wellPathPartsInPolygon =
RigCellGeometryTools::clipPolylineByPolygon(wellpathPart,
cellPolygon,
RigCellGeometryTools::USE_HUGEVAL);
for ( const auto& wellPathPartInCell: wellPathPartsInPolygon )
std::vector<std::vector<cvf::Vec3d>> wellPathPartsInPolygon =
RigCellGeometryTools::clipPolylineByPolygon(wellpathPart, cellPolygon, RigCellGeometryTools::USE_HUGEVAL);
for (const auto& wellPathPartInCell : wellPathPartsInPolygon)
{
if ( !wellPathPartInCell.empty() )
if (!wellPathPartInCell.empty())
{
int endpointCount = 0;
if ( wellPathPartInCell.front().z() != HUGE_VAL ) ++endpointCount;
if ( wellPathPartInCell.back().z() != HUGE_VAL ) ++endpointCount;
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());
double xLengthInCell = fabs(intersectionLength.x());
double yLengthInCell = fabs(intersectionLength.y());
m_stimPlanCellIdxToIntersectionInfoMap[cIdx].endpointCount += endpointCount;
m_stimPlanCellIdxToIntersectionInfoMap[cIdx].hlength += xLengthInCell;