diff --git a/ApplicationCode/FileInterface/RifEclipseExportTools.cpp b/ApplicationCode/FileInterface/RifEclipseExportTools.cpp index 6f7634fe45..838e0154c4 100644 --- a/ApplicationCode/FileInterface/RifEclipseExportTools.cpp +++ b/ApplicationCode/FileInterface/RifEclipseExportTools.cpp @@ -34,6 +34,7 @@ #include #include +#include @@ -59,6 +60,15 @@ RifEclipseExportTools::~RifEclipseExportTools() //-------------------------------------------------------------------------------------------------- bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, const std::vector< RimFracture*>& fractures) { + RiaApplication* app = RiaApplication::instance(); + RimView* activeView = RiaApplication::instance()->activeReservoirView(); + if (!activeView) return false; + RimEclipseView* activeRiv = dynamic_cast(activeView); + if (!activeRiv) return false; + + const RigMainGrid* mainGrid = activeRiv->mainGrid(); + if (!mainGrid) return false; + QFile file(fileName); if (!file.open(QIODevice::WriteOnly | QIODevice::Text)) @@ -66,11 +76,6 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c return false; } - QTextStream out(&file); - out << "\n"; - out << "-- Exported from ResInsight" << "\n\n"; - out << "COMPDAT" << "\n" << right << qSetFieldWidth(8); - caf::ProgressInfo pi(fractures.size(), QString("Writing data to file %1").arg(fileName)); RimEclipseWell* simWell = nullptr; @@ -79,6 +84,89 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c size_t progress =0; std::vector ijk; + QTextStream out(&file); + out << "\n"; + out << "-- Exported from ResInsight" << "\n\n"; + + out << "-- Background data for calculation" << "\n\n"; + + + //Write header line + out << qSetFieldWidth(4); + out << "--"; + out << qSetFieldWidth(12); + out << "Well name"; + + out << qSetFieldWidth(16); + out << "Fracture Name"; + + out << qSetFieldWidth(5); + out << "i"; + out << "j"; + out << "k"; + + out << qSetFieldWidth(8); + out << "Ax"; + out << " "; + out << "Ay"; + out << " "; + out << "Az"; + out << " "; + out << "TotalArea"; + + out << "\n"; + + + + //Write background numbers (with -- first to comment out line) + for (RimFracture* fracture : fractures) + { + fracture->computeTransmissibility(); + std::vector fracDataVector = fracture->attachedRigFracture()->fractureData(); + + for (RigFractureData fracData : fracDataVector) + { + out << qSetFieldWidth(4); + out << "--"; + out << qSetFieldWidth(12); + + wellPath, simWell = nullptr; + fracture->firstAncestorOrThisOfType(simWell); + if (simWell) out << simWell->name; // 1. Well name + fracture->firstAncestorOrThisOfType(wellPath); + if (wellPath) out << wellPath->name; // 1. Well name + + out << qSetFieldWidth(16); + out << fracture->name(); + + + out << qSetFieldWidth(5); + size_t i, j, k; + mainGrid->ijkFromCellIndex(fracData.reservoirCellIndex, &i, &j, &k); + out << i + 1; // 2. I location grid block, adding 1 to go to eclipse 1-based grid definition + out << j + 1; // 3. J location grid block, adding 1 to go to eclipse 1-based grid definition + out << k + 1; // 4. K location of upper connecting grid block, adding 1 to go to eclipse 1-based grid definition + out << " "; + + + + out << fracData.projectedAreas.x(); + out << " "; + out << fracData.projectedAreas.y(); + out << " "; + out << fracData.projectedAreas.z(); + out << " "; + out << fracData.totalArea; + + out << "\n"; + + } + } + + out << "\n"; + + //Write COMPDAT keyword and entries + out << "COMPDAT" << "\n" << right << qSetFieldWidth(8); for (RimFracture* fracture : fractures) { fracture->computeTransmissibility(); @@ -97,15 +185,6 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c out << qSetFieldWidth(5); - RiaApplication* app = RiaApplication::instance(); - RimView* activeView = RiaApplication::instance()->activeReservoirView(); - if (!activeView) return false; - RimEclipseView* activeRiv = dynamic_cast(activeView); - if (!activeRiv) return false; - - const RigMainGrid* mainGrid = activeRiv->mainGrid(); - if (!mainGrid) return false; - size_t i, j, k; mainGrid->ijkFromCellIndex(fracData.reservoirCellIndex, &i, &j, &k); out << i+1; // 2. I location grid block, adding 1 to go to eclipse 1-based grid definition @@ -113,8 +192,9 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c out << k+1; // 4. K location of upper connecting grid block, adding 1 to go to eclipse 1-based grid definition out << k+1; // 5. K location of lower connecting grid block, adding 1 to go to eclipse 1-based grid definition - out << "OPEN"; // 6. Open / Shut flag of connection - out << "1* "; // 7. Saturation table number for connection rel perm. Default value + out << "2*"; // Default value for + //6. Open / Shut flag of connection + // 7. Saturation table number for connection rel perm. Default value out << qSetFieldWidth(8); // 8. Transmissibility @@ -122,29 +202,27 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c else out << "UNDEF"; out << qSetFieldWidth(4); - out << "1* "; // 9. Well bore diameter. Set to default - - out << qSetFieldWidth(8); + out << "2*"; // Default value for + // 9. Well bore diameter. Set to default + // 10. Effective Kh (perm times width) + if (fracture->attachedFractureDefinition()) { - out << fracture->attachedFractureDefinition()->effectiveKh(); // 10. Effective Kh (perm times width) - out << qSetFieldWidth(4); out << fracture->attachedFractureDefinition()->skinFactor; // 11. Skin factor } else //If no attached fracture definition these parameters are set to UNDEF { - out << "UNDEF"; - out << qSetFieldWidth(4); out << "UNDEF"; } - out << "1*"; // 12. D-factor for handling non-Darcy flow of free gas. Default value. - out << "Z"; // 13. Direction well is penetrating the grid block. Z is default. - out << "1*"; // 14. Pressure equivalent radius, Default - - out << "/" << "\n"; + out << "/"; + out << " " << fracture->name(); //Fracture name as comment + out << "\n"; // Terminating entry } + + //TODO: If same cell is used for multiple fractures, the sum of contributions should be added to table. + progress++; pi.setProgress(progress); } diff --git a/ApplicationCode/ProjectDataModel/RimFracture.cpp b/ApplicationCode/ProjectDataModel/RimFracture.cpp index ddd788a129..055d382b51 100644 --- a/ApplicationCode/ProjectDataModel/RimFracture.cpp +++ b/ApplicationCode/ProjectDataModel/RimFracture.cpp @@ -253,8 +253,8 @@ void RimFracture::computeTransmissibility() std::vector fracDataVec; std::vector fracCells = getPotentiallyFracturedCells(); - - for (auto fracCell : fracCells) + + for (size_t fracCell : fracCells) { //TODO: Remove - only for simplifying debugging... RiaApplication* app = RiaApplication::instance(); @@ -274,6 +274,7 @@ void RimFracture::computeTransmissibility() //Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon already is located) cvf::Mat4f invertedTransMatrix = transformMatrix().getInverted(); + for (std::vector & planeCellPolygon : planeCellPolygons) { for (cvf::Vec3d& v : planeCellPolygon) @@ -287,11 +288,18 @@ void RimFracture::computeTransmissibility() localZinFracPlane.transformVector(static_cast(invertedTransMatrix)); cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO; directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1)); + directionOfLength.normalize(); RigFractureData fracData; fracData.reservoirCellIndex = fracCell; double transmissibility; + double fractureArea = 0.0; + double fractureAreaWeightedlength = 0.0; + double Ax = 0.0; + double Ay = 0.0; + double Az = 0.0; + if (attachedFractureDefinition()) { @@ -326,34 +334,26 @@ void RimFracture::computeTransmissibility() areaOfFractureParts.push_back(area); length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, fracturePartPolygon); lengthXareaOfFractureParts.push_back(length * area); - double AreaX = calculateProjectedArea(fracturePartPolygon, localX); - - - //Calculating area in x, y and z direction - cvf::Vec3d planeNormal = cvf::Vec3d::ZERO; - planeNormal.cross(localX, -localZ); - double Ax = calculateProjectedArea(fracturePartPolygon, planeNormal); - - planeNormal.cross(localY, localZ); - double Ay = calculateProjectedArea(fracturePartPolygon, planeNormal); - - planeNormal.cross(localX, localY); - double Az = calculateProjectedArea(fracturePartPolygon, planeNormal); + + cvf::Plane fracturePlane; + cvf::Mat4f m = transformMatrix(); + bool isCellIntersected = false; + fracturePlane.setFromPointAndNormal(static_cast(m.translation()), + static_cast(m.col(2))); + Ax += abs(area*(fracturePlane.normal().dot(localX))); + Ay += abs(area*(fracturePlane.normal().dot(localY))); + Az += abs(area*(fracturePlane.normal().dot(localZ))); } - double totalArea = 0.0; - for (double area : areaOfFractureParts) totalArea += area; + + fractureArea = 0.0; + for (double area : areaOfFractureParts) fractureArea += area; double totalAreaXLength = 0.0; for (double lengtXarea : lengthXareaOfFractureParts) totalAreaXLength += lengtXarea; - double fractureAreaWeightedlength = totalAreaXLength / totalArea; - - - - - + double fractureAreaWeightedlength = totalAreaXLength / fractureArea; //TODO: FInd direction for length calculation (normal to z, in fracture plane) double flowLength = 2.718281828; @@ -361,13 +361,17 @@ void RimFracture::computeTransmissibility() //TODO: read permeability from file (should use matrix permeability and not fracture perm)... - transmissibility = totalArea; - //transmissibility = 8 * c * attachedFractureDefinition()->permeability * areaOfCellPlaneFractureOverlap / + transmissibility = 8 * c * attachedFractureDefinition()->permeability * fractureArea / ( flowLength + (attachedFractureDefinition()->skinFactor * fractureAreaWeightedlength) / cvf::PI_D); } - else transmissibility = cvf::UNDEFINED_DOUBLE; - + else transmissibility = cvf::UNDEFINED_DOUBLE; + fracData.transmissibility = transmissibility; + fracData.fractureLenght = fractureAreaWeightedlength; + fracData.totalArea = fractureArea; + fracData.projectedAreas = cvf::Vec3d(Ax, Ay, Az); + + //only keep fracData if transmissibility is non-zero if (transmissibility > 0) { @@ -438,24 +442,6 @@ bool RimFracture::planeCellIntersectionPolygons(size_t cellindex, std::vector polygon, cvf::Vec3d planeNormal) -{ - //Set up plane - cvf::Plane plane; - plane.setFromPointAndNormal(polygon[0], planeNormal); - - //Project points - for (cvf::Vec3d v : polygon) plane.projectPoint(v); - - //calculate Area - cvf::Vec3d areaVector = cvf::GeometryTools::polygonAreaNormal3D(polygon); - double AreaInPlane = areaVector.length(); - return AreaInPlane; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ProjectDataModel/RimFracture.h b/ApplicationCode/ProjectDataModel/RimFracture.h index c36b705162..200d0cc409 100644 --- a/ApplicationCode/ProjectDataModel/RimFracture.h +++ b/ApplicationCode/ProjectDataModel/RimFracture.h @@ -93,8 +93,6 @@ private: //Functions for area calculations - should these be in separate class bool planeCellIntersectionPolygons(size_t cellindex, std::vector > & polygons, cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ); - - double calculateProjectedArea(std::vector fracturePartPolygon, cvf::Vec3d localX); protected: caf::PdmPtrField m_fractureTemplate; diff --git a/ApplicationCode/ReservoirDataModel/RigFracture.h b/ApplicationCode/ReservoirDataModel/RigFracture.h index 2a1e2d3b90..e789b0affc 100644 --- a/ApplicationCode/ReservoirDataModel/RigFracture.h +++ b/ApplicationCode/ReservoirDataModel/RigFracture.h @@ -33,6 +33,11 @@ public: size_t reservoirCellIndex; double transmissibility; + + double totalArea; + double fractureLenght; + cvf::Vec3d projectedAreas; + };