diff --git a/ApplicationCode/ProjectDataModel/RimFracture.cpp b/ApplicationCode/ProjectDataModel/RimFracture.cpp index 47c5eab00d..6037b56e17 100644 --- a/ApplicationCode/ProjectDataModel/RimFracture.cpp +++ b/ApplicationCode/ProjectDataModel/RimFracture.cpp @@ -243,6 +243,16 @@ void RimFracture::computeTransmissibility() for (auto fracCell : fracCells) { + //TODO: Remove - only for simplifying debugging... + RiaApplication* app = RiaApplication::instance(); + RimView* activeView = RiaApplication::instance()->activeReservoirView(); + RimEclipseView* activeRiv = dynamic_cast(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 > planeCellPolygons; @@ -258,6 +268,13 @@ void RimFracture::computeTransmissibility() { double areaOfCellPlaneFractureOverlap = 0.0; std::vector fracPolygon = attachedFractureDefinition()->fracturePolygon(); + + //TODO: remove if - only intended for adding break point for debug! + if (fracCell == 186234) + { + double c = 0.008527; // TODO: Get value with units, is defined in RimReservoirCellResultsStorage + } + calculateFracturePlaneCellPolygonOverlap(planeCellPolygons, fracPolygon, areaOfCellPlaneFractureOverlap); //TODO: get correct input values... @@ -267,6 +284,7 @@ void RimFracture::computeTransmissibility() double c = 0.008527; // TODO: Get value with units, is defined in RimReservoirCellResultsStorage + transmissibility = 8 * c * attachedFractureDefinition()->permeability * areaOfCellPlaneFractureOverlap / ( flowLength + (attachedFractureDefinition()->skinFactor * fractureLength) / cvf::PI_D); } @@ -274,7 +292,11 @@ void RimFracture::computeTransmissibility() fracData.transmissibility = transmissibility; //only keep fracData if transmissibility is non-zero - if (transmissibility > 0) fracDataVec.push_back(fracData); + if (transmissibility > 0) + { + fracDataVec.push_back(fracData); + } + } m_rigFracture->setFractureData(fracDataVec); @@ -307,6 +329,11 @@ bool RimFracture::planeCellIntersection(size_t cellindex, std::vectorglobalCellArray()[cellindex]; if (cell.isInvalid()) return isCellIntersected; + if (cellindex == 186234) + { + cvf::Vec3d cellcenter = cell.center(); + } + //Copied (and adapted) from RigEclipseWellLogExtractor cvf::Vec3d hexCorners[8]; const std::vector& nodeCoords = mainGrid->nodes(); @@ -374,9 +401,14 @@ bool RimFracture::planeCellIntersection(size_t cellindex, std::vector >::iterator lIt = intersectionLineSegments.begin(); lIt != intersectionLineSegments.end(); lIt++) { - if (lIt->first.equals(polygon.back())) + cvf::Vec3d lineSegmentStart = lIt->first; + cvf::Vec3d lineSegmentEnd = lIt->second; + cvf::Vec3d polygonEnd = polygon.back(); + if(((lineSegmentStart - polygonEnd).lengthSquared() < tolerance*tolerance )) { polygon.push_back(lIt->second); intersectionLineSegments.erase(lIt); @@ -384,7 +416,7 @@ bool RimFracture::planeCellIntersection(size_t cellindex, std::vectorsecond.equals(polygon.back())) + if (((lineSegmentEnd - polygonEnd).lengthSquared() < tolerance*tolerance)) { polygon.push_back(lIt->first); intersectionLineSegments.erase(lIt); @@ -399,7 +431,6 @@ bool RimFracture::planeCellIntersection(size_t cellindex, std::vector clippedPolygon; for (ClipperLib::IntPoint IntPosition : pathInSol) { - cvf::Vec3d v; - v.x() = IntPosition.X / polygonScaleFactor; - v.y() = IntPosition.Y / polygonScaleFactor; + 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); } - //We should only have one solution polygon / path!?!? //calculate area for (std::vector areaPolygon : clippedPolygons)