#1776 Updating calculation of length of stimplan-eclipse-overlap polygon for sl/pi-term in matrix-to-fracture transmissibility to avoid possible NaN values in export. Length is now calculated in x direction in the fracture coordinate system.

This commit is contained in:
astridkbjorke 2017-08-23 09:46:00 +02:00
parent 0bb892541a
commit 5134a461d9
5 changed files with 41 additions and 59 deletions

View File

@ -153,11 +153,8 @@ void RigCellGeometryTools::findCellLocalXYZ(const std::array<cvf::Vec3d, 8>& hex
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
//TODO: Suggested rename: polygonLengthWeightedByArea. polygonAverageLengthWeightedByArea
double RigCellGeometryTools::polygonAreaWeightedLength(cvf::Vec3d directionOfLength, std::vector<cvf::Vec3d> polygonToCalcLengthOf)
double RigCellGeometryTools::polygonLengthInLocalXdirWeightedByArea(std::vector<cvf::Vec3d> polygonToCalcLengthOf)
{
//TODO: Check that polygon is in xy plane
//Find bounding box
cvf::BoundingBox polygonBBox;
for (cvf::Vec3d nodeCoord : polygonToCalcLengthOf) polygonBBox.add(nodeCoord);
@ -167,16 +164,25 @@ double RigCellGeometryTools::polygonAreaWeightedLength(cvf::Vec3d directionOfLen
//Split bounding box in multiple polygons (2D)
int resolutionOfLengthCalc = 20;
cvf::Vec3d widthOfPolygon = polygonBBox.extent() / resolutionOfLengthCalc;
double widthOfPolygon = polygonBBox.extent().y() / resolutionOfLengthCalc;
std::vector<double> areasOfPolygonContributions;
std::vector<double> lengthOfPolygonContributions;
cvf::Vec3d directionOfLength(1, 0, 0);
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));
cvf::Vec3d pointOnLine1(bboxCorners[0].x(), bboxCorners[0].y() + i*widthOfPolygon, 0);
cvf::Vec3d pointOnLine2(bboxCorners[0].x(), bboxCorners[0].y() + (i + 1)*widthOfPolygon, 0);
std::pair<cvf::Vec3d, cvf::Vec3d> line1 = getLineThroughBoundingBox(directionOfLength,
polygonBBox,
pointOnLine1);
std::pair<cvf::Vec3d, cvf::Vec3d> line2 = getLineThroughBoundingBox(directionOfLength,
polygonBBox,
pointOnLine2);
std::vector<cvf::Vec3d> polygon;
polygon.push_back(line1.first);
polygon.push_back(line1.second);
@ -195,7 +201,8 @@ double RigCellGeometryTools::polygonAreaWeightedLength(cvf::Vec3d directionOfLen
{
areaVector = cvf::GeometryTools::polygonAreaNormal3D(clippedPolygon);
area += areaVector.length();
length += getLengthOfPolygonAlongLine(line1, clippedPolygon); //For testing that parts of code fit together...
length += (getLengthOfPolygonAlongLine(line1, clippedPolygon)
+ getLengthOfPolygonAlongLine(line2, clippedPolygon)) / 2;
}
areasOfPolygonContributions.push_back(area);
lengthOfPolygonContributions.push_back(length);

View File

@ -40,7 +40,7 @@ public:
static void findCellLocalXYZ(const std::array<cvf::Vec3d, 8>& hexCorners, cvf::Vec3d &localXdirection, cvf::Vec3d &localYdirection, cvf::Vec3d &localZdirection);
static double polygonAreaWeightedLength(cvf::Vec3d directionOfLength, std::vector<cvf::Vec3d> polygon2d);
static double polygonLengthInLocalXdirWeightedByArea(std::vector<cvf::Vec3d> polygon2d);
static std::vector<std::vector<cvf::Vec3d> > intersectPolygons(std::vector<cvf::Vec3d> polygon1, std::vector<cvf::Vec3d> polygon2);

View File

@ -146,17 +146,6 @@ void RigEclipseToStimPlanCellTransmissibilityCalculator::calculateStimPlanCellsM
}
}
cvf::Vec3d localZinFracPlane;
localZinFracPlane = localZ;
localZinFracPlane.transformVector(invertedTransMatrix);
cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO;
directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1));
directionOfLength.normalize();
//Fracture plane is XY - so localZinFracPlane is in this plane.
//Crossing this vector with a vector normal to this plane (0,0,1) gives a vector for in the ij-direction in the frac plane
//This is the direction in which we calculate the length of the fracture element in the cell,
//to use in the skinfactor contribution (S*l/pi) to the transmissibility.
std::vector<std::vector<cvf::Vec3d> > polygonsForStimPlanCellInEclipseCell;
cvf::Vec3d areaVector;
std::vector<cvf::Vec3d> stimPlanPolygon = m_stimPlanCell.getPolygon();
@ -184,8 +173,7 @@ void RigEclipseToStimPlanCellTransmissibilityCalculator::calculateStimPlanCellsM
area = areaVector.length();
areaOfFractureParts.push_back(area);
//TODO: should the l in the sl/pi term in the denominator of the Tmj expression be the length of the full Eclipse cell or fracture?
length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, fracturePartPolygon);
length = RigCellGeometryTools::polygonLengthInLocalXdirWeightedByArea(fracturePartPolygon);
lengthXareaOfFractureParts.push_back(length * area);
cvf::Plane fracturePlane;

View File

@ -113,8 +113,14 @@ double RigFractureTransmissibilityEquations::matrixToFractureTrans(double perm,
{
double transmissibility;
double slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D;
double slDivPi = 0.0;
if (abs(skinfactor) > 1e-9)
{
slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D;
}
transmissibility = 8 * cDarcy * (perm * NTG) * A / (cellSizeLength + slDivPi);
CVF_ASSERT(transmissibility == transmissibility);
return transmissibility;
}

View File

@ -149,51 +149,32 @@ 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);
double length = RigCellGeometryTools::polygonLengthInLocalXdirWeightedByArea(polygonExample);
EXPECT_DOUBLE_EQ(length, 1.5);
directionOfLength = cvf::Vec3d(0, 1, 0);
length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, polygonExample);
EXPECT_DOUBLE_EQ(length, 2.5);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST(RigCellGeometryTools, lengthCalcTestTriangle)
{
std::vector<cvf::Vec3d> trianglePolygonExample;
trianglePolygonExample.push_back(cvf::Vec3d(0.00, 0.00, 0.0));
trianglePolygonExample.push_back(cvf::Vec3d(2.50, 2.50, 0.0));
trianglePolygonExample.push_back(cvf::Vec3d(2.50, 0.00, 0.0));
double length = RigCellGeometryTools::polygonLengthInLocalXdirWeightedByArea(trianglePolygonExample);
EXPECT_GT(length, 1.7);
EXPECT_LT(length, 1.8);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------