mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
#1091 - pre-proto - Adding function for calculating fracture length.
This commit is contained in:
astridkbjorke
parent
aedf184d14
commit
9c9592ccce
ApplicationCode
ProjectDataModel
ReservoirDataModel
UnitTests
@ -261,20 +261,15 @@ void RimFracture::computeTransmissibility()
|
||||
RimView* activeView = RiaApplication::instance()->activeReservoirView();
|
||||
RimEclipseView* activeRiv = dynamic_cast<RimEclipseView*>(activeView);
|
||||
const RigMainGrid* mainGrid = activeRiv->mainGrid();
|
||||
|
||||
size_t i, j, k;
|
||||
mainGrid->ijkFromCellIndex(fracCell, &i, &j, &k);
|
||||
|
||||
//End of code only for debugging...
|
||||
|
||||
|
||||
std::vector<std::vector<cvf::Vec3d> > planeCellPolygons;
|
||||
cvf::Vec3d localX;
|
||||
cvf::Vec3d localY;
|
||||
cvf::Vec3d localZ;
|
||||
|
||||
std::vector<std::vector<cvf::Vec3d> > planeCellPolygons;
|
||||
bool isPlanIntersected = planeCellIntersectionPolygons(fracCell, planeCellPolygons, localX, localY, localZ);
|
||||
|
||||
if (!isPlanIntersected || planeCellPolygons.size()==0) continue;
|
||||
|
||||
//Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon already is located)
|
||||
@ -287,9 +282,11 @@ void RimFracture::computeTransmissibility()
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
//TODO: Make copy of z dir vector, we need it both in fracture coords and domain cords
|
||||
localZ.transformVector(static_cast<cvf::Mat4d>(invertedTransMatrix));
|
||||
cvf::Vec3d localZinFracPlane;
|
||||
localZinFracPlane = localZ;
|
||||
localZinFracPlane.transformVector(static_cast<cvf::Mat4d>(invertedTransMatrix));
|
||||
cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO;
|
||||
directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1));
|
||||
|
||||
RigFractureData fracData;
|
||||
fracData.reservoirCellIndex = fracCell;
|
||||
@ -301,13 +298,26 @@ void RimFracture::computeTransmissibility()
|
||||
double areaOfCellPlaneFractureOverlap = 0.0;
|
||||
std::vector<cvf::Vec3f> fracPolygon = attachedFractureDefinition()->fracturePolygon();
|
||||
|
||||
calculateFracturePlaneCellPolygonOverlapArea(planeCellPolygons, fracPolygon, areaOfCellPlaneFractureOverlap);
|
||||
|
||||
std::vector<std::vector<cvf::Vec3d> > clippedPolygons;
|
||||
calculateFracturePlaneCellPolygonOverlapArea(planeCellPolygons, fracPolygon, clippedPolygons, areaOfCellPlaneFractureOverlap);
|
||||
|
||||
//TODO: get correct input values...
|
||||
double fractureLength = 1.2345;
|
||||
//TODO: FInd direction for length calculation (normal to z, in fracture plane)
|
||||
double flowLength = 2.718281828;
|
||||
|
||||
double c = 0.008527; // TODO: Get value with units, is defined in RimReservoirCellResultsStorage
|
||||
|
||||
double fractureLength = 0.0;
|
||||
for (std::vector<cvf::Vec3d> clippedPolygon : clippedPolygons)
|
||||
//TODO: Should probably do some something more accurate here...
|
||||
{
|
||||
fractureLength += RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, clippedPolygon);
|
||||
}
|
||||
|
||||
|
||||
|
||||
double c = 0.008527; // TODO: Get value with units, is defined in RimReservoirCellResultsStorage
|
||||
|
||||
|
||||
//TODO: read permeability from file (should use matrix permeability and not fracture perm)...
|
||||
transmissibility = 8 * c * attachedFractureDefinition()->permeability * areaOfCellPlaneFractureOverlap /
|
||||
@ -390,7 +400,7 @@ bool RimFracture::planeCellIntersectionPolygons(size_t cellindex, std::vector<st
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RimFracture::calculateFracturePlaneCellPolygonOverlapArea(std::vector<std::vector<cvf::Vec3d> > planeCellPolygons,
|
||||
std::vector<cvf::Vec3f> fracturePolygon, double & areaOfCellPlaneFractureOverlap)
|
||||
std::vector<cvf::Vec3f> fracturePolygon, std::vector<std::vector<cvf::Vec3d> > & clippedPolygons, double & areaOfCellPlaneFractureOverlap)
|
||||
{
|
||||
int polygonScaleFactor = 10000; //For transform to clipper int
|
||||
int xInt, yInt;
|
||||
@ -423,7 +433,7 @@ void RimFracture::calculateFracturePlaneCellPolygonOverlapArea(std::vector<std::
|
||||
clpr.Execute(ClipperLib::ctIntersection, solution, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);
|
||||
|
||||
//Convert back to std::vector<std::vector<cvf::Vec3d> >
|
||||
std::vector<std::vector<cvf::Vec3d> > clippedPolygons;
|
||||
//std::vector<std::vector<cvf::Vec3d> > clippedPolygons;
|
||||
for (ClipperLib::Path pathInSol : solution)
|
||||
{
|
||||
std::vector<cvf::Vec3d> clippedPolygon;
|
||||
|
||||
@ -95,7 +95,7 @@ private:
|
||||
bool planeCellIntersectionPolygons(size_t cellindex, std::vector<std::vector<cvf::Vec3d> > & polygons, cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ);
|
||||
|
||||
void calculateFracturePlaneCellPolygonOverlapArea(std::vector<std::vector<cvf::Vec3d> > planeCellPolygons,
|
||||
std::vector<cvf::Vec3f> fracturePolygon, double & area);
|
||||
std::vector<cvf::Vec3f> fracturePolygon, std::vector<std::vector<cvf::Vec3d> > & clippedPolygons, double & area);
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
@ -23,6 +23,9 @@
|
||||
|
||||
#include "cafHexGridIntersectionTools/cafHexGridIntersectionTools.h"
|
||||
#include "cvfBoundingBox.h"
|
||||
#include "custom-clipper/clipper/clipper.hpp"
|
||||
|
||||
#include <vector>
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
@ -186,26 +189,165 @@ void RigCellGeometryTools::findCellLocalXYZ(cvf::Vec3d * hexCorners, cvf::Vec3d
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RigCellGeometryTools::polygonAreaWeightedLength(cvf::Vec3d directionOfLength, std::vector<cvf::Vec3d> polygon2d)
|
||||
double RigCellGeometryTools::polygonAreaWeightedLength(cvf::Vec3d directionOfLength, std::vector<cvf::Vec3d> polygonToCalcLengthOf)
|
||||
{
|
||||
//TODO: Check that polygon is in xy plane
|
||||
|
||||
//Find bounding box
|
||||
cvf::BoundingBox polygonBBox;
|
||||
for (cvf::Vec3d nodeCoord : polygon2d) polygonBBox.add(nodeCoord);
|
||||
for (cvf::Vec3d nodeCoord : polygonToCalcLengthOf) polygonBBox.add(nodeCoord);
|
||||
cvf::Vec3d bboxCorners[8];
|
||||
polygonBBox.cornerVertices(bboxCorners);
|
||||
|
||||
|
||||
//Split bounding box in multplie polygons (2D)
|
||||
//Split bounding box in multiple polygons (2D)
|
||||
int resolutionOfLengthCalc = 20;
|
||||
cvf::Vec3d widthOfPolygon = polygonBBox.extent() / resolutionOfLengthCalc;
|
||||
|
||||
std::vector<double> areasOfPolygonContributions;
|
||||
std::vector<double> lengthOfPolygonContributions;
|
||||
|
||||
//Use clipper to find overlap between bbpolygon and fracture
|
||||
for (int i = 0; i < resolutionOfLengthCalc; i++)
|
||||
{
|
||||
//Doing the same thing twice, this part can be optimized...
|
||||
std::pair<cvf::Vec3d, cvf::Vec3d> line1 = getLineThroughBoundingBox(directionOfLength, polygonBBox, bboxCorners[0] + widthOfPolygon*i);
|
||||
std::pair<cvf::Vec3d, cvf::Vec3d> line2 = getLineThroughBoundingBox(directionOfLength, polygonBBox, bboxCorners[0] + widthOfPolygon*(i+1));
|
||||
std::vector<cvf::Vec3d> polygon;
|
||||
polygon.push_back(line1.first);
|
||||
polygon.push_back(line1.second);
|
||||
polygon.push_back(line2.second);
|
||||
polygon.push_back(line2.first);
|
||||
|
||||
//Calculate length (max-min) and area
|
||||
//Use clipper to find overlap between bbpolygon and fracture
|
||||
std::vector<std::vector<cvf::Vec3d> > clippedPolygons = clipPolygons(polygonToCalcLengthOf, polygon);
|
||||
|
||||
//Calculate area-weighted average of above vectors.
|
||||
double area = 0;
|
||||
double length = 0;
|
||||
cvf::Vec3d areaVector = cvf::Vec3d::ZERO;
|
||||
|
||||
//Calculate length (max-min) and area
|
||||
for (std::vector<cvf::Vec3d> clippedPolygon : clippedPolygons)
|
||||
{
|
||||
areaVector = cvf::GeometryTools::polygonAreaNormal3D(clippedPolygon);
|
||||
area += areaVector.length();
|
||||
length += getLengthOfPolygonAlongLine(line1, clippedPolygon); //For testing that parts of code fit together...
|
||||
}
|
||||
areasOfPolygonContributions.push_back(area);
|
||||
lengthOfPolygonContributions.push_back(length);
|
||||
}
|
||||
|
||||
//Calculate area-weighted length average.
|
||||
double totalArea = 0.0;
|
||||
double totalAreaXlength = 0.0;
|
||||
|
||||
for (int i = 0; i < areasOfPolygonContributions.size(); i++)
|
||||
{
|
||||
totalArea += areasOfPolygonContributions[i];
|
||||
totalAreaXlength += (areasOfPolygonContributions[i] * lengthOfPolygonContributions[i]);
|
||||
}
|
||||
|
||||
double areaWeightedLength = totalAreaXlength / totalArea;
|
||||
return areaWeightedLength;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<std::vector<cvf::Vec3d> > RigCellGeometryTools::clipPolygons(std::vector<cvf::Vec3d> polygon1, std::vector<cvf::Vec3d> polygon2)
|
||||
{
|
||||
int polygonScaleFactor = 10000; //For transform to clipper int
|
||||
int xInt, yInt;
|
||||
|
||||
//Convert to int for clipper library and store as clipper "path"
|
||||
ClipperLib::Path polygon1path;
|
||||
for (cvf::Vec3d& v : polygon1)
|
||||
{
|
||||
xInt = v.x()*polygonScaleFactor;
|
||||
yInt = v.y()*polygonScaleFactor;
|
||||
polygon1path.push_back(ClipperLib::IntPoint(xInt, yInt));
|
||||
}
|
||||
|
||||
ClipperLib::Path polygon2path;
|
||||
for (cvf::Vec3d& v : polygon2)
|
||||
{
|
||||
xInt = v.x()*polygonScaleFactor;
|
||||
yInt = v.y()*polygonScaleFactor;
|
||||
polygon2path.push_back(ClipperLib::IntPoint(xInt, yInt));
|
||||
}
|
||||
|
||||
ClipperLib::Clipper clpr;
|
||||
clpr.AddPath(polygon1path, ClipperLib::ptSubject, true);
|
||||
clpr.AddPath(polygon2path, ClipperLib::ptClip, true);
|
||||
|
||||
ClipperLib::Paths solution;
|
||||
clpr.Execute(ClipperLib::ctIntersection, solution, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);
|
||||
|
||||
//Convert back to std::vector<std::vector<cvf::Vec3d> >
|
||||
std::vector<std::vector<cvf::Vec3d> > clippedPolygons;
|
||||
for (ClipperLib::Path pathInSol : solution)
|
||||
{
|
||||
std::vector<cvf::Vec3d> clippedPolygon;
|
||||
for (ClipperLib::IntPoint IntPosition : pathInSol)
|
||||
{
|
||||
cvf::Vec3d v = cvf::Vec3d::ZERO;
|
||||
v.x() = (float)IntPosition.X / (float)polygonScaleFactor;
|
||||
v.y() = (float)IntPosition.Y / (float)polygonScaleFactor;
|
||||
clippedPolygon.push_back(v);
|
||||
}
|
||||
clippedPolygons.push_back(clippedPolygon);
|
||||
}
|
||||
|
||||
return clippedPolygons;
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::pair<cvf::Vec3d, cvf::Vec3d> RigCellGeometryTools::getLineThroughBoundingBox(cvf::Vec3d lineDirection, cvf::BoundingBox polygonBBox, cvf::Vec3d pointOnLine)
|
||||
{
|
||||
cvf::Vec3d bboxCorners[8];
|
||||
polygonBBox.cornerVertices(bboxCorners);
|
||||
|
||||
cvf::Vec3d startPoint = pointOnLine;
|
||||
cvf::Vec3d endPoint = pointOnLine;
|
||||
|
||||
|
||||
//To avoid doing many iterations in loops below linedirection should be quite large.
|
||||
lineDirection.normalize();
|
||||
lineDirection = lineDirection * polygonBBox.extent().length() / 5;
|
||||
|
||||
//Extend line in positive direction
|
||||
while (polygonBBox.contains(startPoint))
|
||||
{
|
||||
startPoint = startPoint + lineDirection;
|
||||
}
|
||||
//Extend line in negative direction
|
||||
while (polygonBBox.contains(endPoint))
|
||||
{
|
||||
endPoint = endPoint - lineDirection;
|
||||
}
|
||||
|
||||
std::pair<cvf::Vec3d, cvf::Vec3d> line;
|
||||
line = { startPoint, endPoint };
|
||||
return line;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigCellGeometryTools::getLengthOfPolygonAlongLine(std::pair<cvf::Vec3d, cvf::Vec3d> line, std::vector<cvf::Vec3d> polygon)
|
||||
{
|
||||
cvf::BoundingBox lineBoundingBox;
|
||||
|
||||
std::vector<cvf::Vec3d> pointsOnLine;
|
||||
|
||||
for (cvf::Vec3d polygonPoint : polygon)
|
||||
{
|
||||
cvf::Vec3d pointOnLine = cvf::GeometryTools::projectPointOnLine(line.first, line.second, polygonPoint, nullptr);
|
||||
lineBoundingBox.add(pointOnLine);
|
||||
}
|
||||
|
||||
double length = lineBoundingBox.extent().length();
|
||||
return length;
|
||||
}
|
||||
|
||||
@ -21,9 +21,11 @@
|
||||
#include "cvfBase.h"
|
||||
#include "cvfVector3.h"
|
||||
|
||||
#include "cvfBoundingBox.h"
|
||||
#include "cvfPlane.h"
|
||||
|
||||
#include <list>
|
||||
#include <vector>
|
||||
|
||||
|
||||
|
||||
@ -39,5 +41,11 @@ public:
|
||||
|
||||
static void findCellLocalXYZ(cvf::Vec3d * hexCorners, cvf::Vec3d &localXdirection, cvf::Vec3d &localYdirection, cvf::Vec3d &localZdirection);
|
||||
|
||||
static void polygonAreaWeightedLength(cvf::Vec3d directionOfLength, std::vector<cvf::Vec3d> polygon2d);
|
||||
static double polygonAreaWeightedLength(cvf::Vec3d directionOfLength, std::vector<cvf::Vec3d> polygon2d);
|
||||
|
||||
static std::vector<std::vector<cvf::Vec3d> > clipPolygons(std::vector<cvf::Vec3d> polygon1, std::vector<cvf::Vec3d> polygon2);
|
||||
|
||||
static std::pair<cvf::Vec3d, cvf::Vec3d> getLineThroughBoundingBox(cvf::Vec3d lineDirection, cvf::BoundingBox polygonBBox, cvf::Vec3d pointOnLine);
|
||||
|
||||
static double getLengthOfPolygonAlongLine(std::pair<cvf::Vec3d, cvf::Vec3d> line, std::vector<cvf::Vec3d> polygon);
|
||||
};
|
||||
|
||||
@ -143,3 +143,51 @@ TEST(RigCellGeometryTools, findCellAverageZTest)
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST(RigCellGeometryTools, lengthCalcTest)
|
||||
{
|
||||
|
||||
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));
|
||||
|
||||
|
||||
// polygonExample.push_back(cvf::Vec3d(0.00, 0.50, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(-7.73, 0.48, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(-14.69, 0.40, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(-20.23, 0.29, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(-23.78, 0.15, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(-25.00, 0.00, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(-23.78, -0.15, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(-20.23, -0.29, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(-14.69, -0.40, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(-7.73, -0.48, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(0.00, -0.50, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(7.73, -0.48, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(14.69, -0.40, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(20.23, -0.29, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(23.78, -0.15, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(25.00, 0.00, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(23.78, 0.15, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(20.23, 0.29, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(14.69, 0.40, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(7.73, 0.48, 0.0));
|
||||
// polygonExample.push_back(cvf::Vec3d(0.00, 0.50, 0.0));
|
||||
|
||||
double length = 0.0;
|
||||
cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO;
|
||||
|
||||
|
||||
directionOfLength = cvf::Vec3d(1, 0, 0);
|
||||
length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, polygonExample);
|
||||
EXPECT_DOUBLE_EQ(length, 1.5);
|
||||
|
||||
directionOfLength = cvf::Vec3d(0, 1, 0);
|
||||
length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, polygonExample);
|
||||
EXPECT_DOUBLE_EQ(length, 2.5);
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user