diff --git a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp index 5032ee2cb1..a535375494 100644 --- a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp +++ b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp @@ -344,6 +344,10 @@ std::vector > RigCellGeometryTools::clipPolylineByPolygo const std::vector& polygon, ZInterpolationType interpolType) { + + //Adjusting polygon to avoid clipper issue with interpolating z-values when lines crosses though polygon vertecies + std::vector adjustedPolygon = ajustPolygonToAvoidIntersectionsAtVertex(polyLine, polygon); + int polygonScaleFactor = 10000; //For transform to clipper int //Convert to int for clipper library and store as clipper "path" @@ -354,7 +358,7 @@ std::vector > RigCellGeometryTools::clipPolylineByPolygo } ClipperLib::Path polygonPath; - for (const cvf::Vec3d& v : polygon) + for (const cvf::Vec3d& v : adjustedPolygon) { ClipperLib::IntPoint intp = toClipperPoint(v); intp.Z = std::numeric_limits::max(); @@ -446,3 +450,44 @@ double RigCellGeometryTools::getLengthOfPolygonAlongLine(std::pair RigCellGeometryTools::ajustPolygonToAvoidIntersectionsAtVertex(const std::vector& polyLine, + const std::vector& polygon) +{ + std::vector adjustedPolygon; + + double treshold = (1.0 / 10000.0) * 5; //5 times polygonScaleFactor for converting to int for clipper + + for (cvf::Vec3d polygonPoint : polygon) + { + + for (int i = 0; i < polyLine.size() - 1; i++) + { + cvf::Vec3d linePoint1(polyLine[i].x(), polyLine[i].y(), 0.0); + cvf::Vec3d linePoint2(polyLine[i + 1].x(), polyLine[i + 1].y(), 0.0); + + double pointDistanceFromLine = cvf::GeometryTools::linePointSquareDist(linePoint1, linePoint2, polygonPoint); + if (pointDistanceFromLine < treshold) + { + //calculate new polygonPoint + cvf::Vec3d directionOfLineSegment = linePoint2 - linePoint1; + //finding normal to the direction of the line segment in the XY plane (z=0) + cvf::Vec3d normalToLine(-directionOfLineSegment.y(), directionOfLineSegment.x(), 0.0); + normalToLine.normalize(); + polygonPoint = polygonPoint + normalToLine * 0.005; + } + } + adjustedPolygon.push_back(polygonPoint); + } + + return adjustedPolygon; +} + + + + + + diff --git a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h index 555a5bad79..90bcc736ab 100644 --- a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h +++ b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h @@ -52,4 +52,9 @@ public: static std::pair getLineThroughBoundingBox(cvf::Vec3d lineDirection, cvf::BoundingBox polygonBBox, cvf::Vec3d pointOnLine); static double getLengthOfPolygonAlongLine(std::pair line, std::vector polygon); + +private: + static std::vector ajustPolygonToAvoidIntersectionsAtVertex(const std::vector& polyLine, + const std::vector& polygon); + }; diff --git a/ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.cpp b/ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.cpp index f79d197c80..52814562c3 100644 --- a/ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.cpp +++ b/ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.cpp @@ -67,7 +67,9 @@ void RigWellPathStimplanIntersector::calculate(const cvf::Mat4d &fractureXf, // Clip well path to fracture domain std::vector > wellPathPartsWithinFracture = - RigCellGeometryTools::clipPolylineByPolygon(fractureRelativeWellPathPoints, perforationLengthBoundingBoxPolygon, RigCellGeometryTools::INTERPOLATE_LINE_Z); + 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