///////////////////////////////////////////////////////////////////////////////// // // Copyright (C) 2017 Statoil ASA // // ResInsight is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY // WARRANTY; without even the implied warranty of MERCHANTABILITY or // FITNESS FOR A PARTICULAR PURPOSE. // // See the GNU General Public License at // for more details. // ///////////////////////////////////////////////////////////////////////////////// #include "RimFracture.h" #include "RiaApplication.h" #include "RifReaderInterface.h" #include "RigCaseCellResultsData.h" #include "RigCell.h" #include "RigCellGeometryTools.h" #include "RigEclipseCaseData.h" #include "RigFracture.h" #include "RigMainGrid.h" #include "RigResultAccessor.h" #include "RigResultAccessorFactory.h" #include "RigTesselatorTools.h" #include "RimEclipseCase.h" #include "RimEclipseCellColors.h" #include "RimEclipseView.h" #include "RimEllipseFractureTemplate.h" #include "RimFractureDefinitionCollection.h" #include "RimOilField.h" #include "RimProject.h" #include "RimView.h" #include "RivWellFracturePartMgr.h" #include "cafPdmUiDoubleSliderEditor.h" #include "cafHexGridIntersectionTools/cafHexGridIntersectionTools.h" #include "cvfBoundingBox.h" #include "cvfGeometryTools.h" #include "cvfMath.h" #include "cvfMatrix4.h" #include "cvfPlane.h" #include "clipper/clipper.hpp" #include #include "RimReservoirCellResultsStorage.h" #include "RigActiveCellInfo.h" CAF_PDM_XML_ABSTRACT_SOURCE_INIT(RimFracture, "Fracture"); //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RimFracture::RimFracture(void) { CAF_PDM_InitObject("Fracture", "", "", ""); CAF_PDM_InitField(&name, "Name", QString("Fracture Name"), "Name", "", "", ""); CAF_PDM_InitFieldNoDefault(&m_fractureTemplate, "FractureDef", "Fracture Template", "", "", ""); CAF_PDM_InitFieldNoDefault(&m_anchorPosition, "anchorPosition", "Anchor Position", "", "", ""); m_anchorPosition.uiCapability()->setUiHidden(true); CAF_PDM_InitFieldNoDefault(&m_uiAnchorPosition, "ui_positionAtWellpath", "Fracture Position", "", "", ""); m_uiAnchorPosition.registerGetMethod(this, &RimFracture::fracturePositionForUi); m_uiAnchorPosition.uiCapability()->setUiReadOnly(true); CAF_PDM_InitField(&azimuth, "Azimuth", 0.0, "Azimuth", "", "", ""); CAF_PDM_InitField(&perforationLength, "PerforationLength", 0.0, "PerforationLength", "", "", ""); azimuth.uiCapability()->setUiEditorTypeName(caf::PdmUiDoubleSliderEditor::uiEditorTypeName()); CAF_PDM_InitField(&m_i, "I", 1, "Fracture location cell I", "", "", ""); m_i.uiCapability()->setUiHidden(true); CAF_PDM_InitField(&m_j, "J", 1, "Fracture location cell J", "", "", ""); m_j.uiCapability()->setUiHidden(true); CAF_PDM_InitField(&m_k, "K", 1, "Fracture location cell K", "", "", ""); m_k.uiCapability()->setUiHidden(true); CAF_PDM_InitFieldNoDefault(&m_displayIJK, "Cell_IJK", "Cell IJK", "", "", ""); m_displayIJK.registerGetMethod(this, &RimFracture::createOneBasedIJK); m_displayIJK.uiCapability()->setUiReadOnly(true); m_rigFracture = new RigFracture; m_recomputeGeometry = true; m_rivFracture = new RivWellFracturePartMgr(this); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RimFracture::~RimFracture() { } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const std::vector& RimFracture::triangleIndices() const { return m_rigFracture->triangleIndices(); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const std::vector& RimFracture::nodeCoords() const { return m_rigFracture->nodeCoords(); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector RimFracture::getPotentiallyFracturedCells() { std::vector cellindecies; RiaApplication* app = RiaApplication::instance(); RimView* activeView = RiaApplication::instance()->activeReservoirView(); if (!activeView) return cellindecies; RimEclipseView* activeRiv = dynamic_cast(activeView); if (!activeRiv) return cellindecies; const RigMainGrid* mainGrid = activeRiv->mainGrid(); if (!mainGrid) return cellindecies; const std::vector& nodeCoordVec = nodeCoords(); if (!hasValidGeometry()) computeGeometry(); cvf::BoundingBox polygonBBox; for (cvf::Vec3f nodeCoord : nodeCoordVec) polygonBBox.add(nodeCoord); mainGrid->findIntersectingCells(polygonBBox, &cellindecies); return cellindecies; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFracture::fieldChangedByUi(const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue) { if (changedField == &m_fractureTemplate) { perforationLength = m_fractureTemplate->perforationLength(); azimuth = m_fractureTemplate->azimuthAngle(); } if (changedField == &azimuth || changedField == &m_fractureTemplate) { setRecomputeGeometryFlag(); RimView* rimView = nullptr; this->firstAncestorOrThisOfType(rimView); if (rimView) { rimView->createDisplayModelAndRedraw(); } else { // Can be triggered from well path, find active view RimView* activeView = RiaApplication::instance()->activeReservoirView(); if (activeView) { activeView->createDisplayModelAndRedraw(); } } } } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- caf::PdmFieldHandle* RimFracture::userDescriptionField() { return &name; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFracture::computeGeometry() { std::vector nodeCoords; std::vector polygonIndices; RimEllipseFractureTemplate* fractureDef = attachedFractureDefinition(); if (fractureDef ) { fractureDef->fractureGeometry(&nodeCoords, &polygonIndices); } cvf::Mat4f m = transformMatrix(); for (cvf::Vec3f& v : nodeCoords) { v.transformPoint(m); } m_rigFracture->setGeometry(polygonIndices, nodeCoords); m_recomputeGeometry = false; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- cvf::Vec3d RimFracture::anchorPosition() { return m_anchorPosition(); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- cvf::Mat4f RimFracture::transformMatrix() { cvf::Vec3d center = anchorPosition(); // Ellipsis geometry is produced in XY-plane, rotate 90 deg around X to get zero azimuth along Y cvf::Mat4f rotationFromTesselator = cvf::Mat4f::fromRotation(cvf::Vec3f::X_AXIS, cvf::Math::toRadians(90.0f)); // Azimuth rotation cvf::Mat4f azimuthRotation = cvf::Mat4f::fromRotation(cvf::Vec3f::Z_AXIS, cvf::Math::toRadians(-azimuth())); cvf::Mat4f m = azimuthRotation * rotationFromTesselator; m.setTranslation(cvf::Vec3f(center)); return m; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFracture::computeTransmissibility() { RiaApplication* app = RiaApplication::instance(); RimView* activeView = RiaApplication::instance()->activeReservoirView(); RimEclipseView* activeRiv = dynamic_cast(activeView); caf::PdmObjectHandle* objHandle = dynamic_cast(this); if (!objHandle) return; RimEclipseCase* eclipseCase; objHandle->firstAncestorOrThisOfType(eclipseCase); RigEclipseCaseData* eclipseCaseData = eclipseCase->reservoirData(); RimEclipseCellColors* resultColors = activeRiv->cellResult(); RimReservoirCellResultsStorage* gridCellResults = resultColors->currentGridCellResults(); RifReaderInterface::PorosityModelResultType porosityModel = RigCaseCellResultsData::convertFromProjectModelPorosityModel(resultColors->porosityModel()); size_t scalarSetIndex; scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DX"); cvf::ref dataAccessObjectDx = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DX"); //assuming 0 time step and main grid (so grid index =0) scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DY"); cvf::ref dataAccessObjectDy = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DY"); //assuming 0 time step and main grid (so grid index =0) scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DZ"); cvf::ref dataAccessObjectDz = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DZ"); //assuming 0 time step and main grid (so grid index =0) scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMX"); cvf::ref dataAccessObjectPermX = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMX"); //assuming 0 time step and main grid (so grid index =0) scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMY"); cvf::ref dataAccessObjectPermY = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMY"); //assuming 0 time step and main grid (so grid index =0) scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMZ"); cvf::ref dataAccessObjectPermZ = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMZ"); //assuming 0 time step and main grid (so grid index =0) scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "NTG"); cvf::ref dataAccessObjectNTG = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "NTG"); //assuming 0 time step and main grid (so grid index =0) RigActiveCellInfo* activeCellInfo = eclipseCaseData->activeCellInfo(porosityModel); std::vector fracDataVec; std::vector fracCells = getPotentiallyFracturedCells(); for (size_t fracCell : fracCells) { //TODO: Remove - only for simplifying debugging... const RigMainGrid* mainGrid = activeRiv->mainGrid(); size_t i, j, k; mainGrid->ijkFromCellIndex(fracCell, &i, &j, &k); bool cellIsActive = activeCellInfo->isActive(fracCell); double permX = dataAccessObjectPermX->cellScalarGlobIdx(fracCell); double permY = dataAccessObjectPermY->cellScalarGlobIdx(fracCell); double permZ = dataAccessObjectPermZ->cellScalarGlobIdx(fracCell); double dx = dataAccessObjectDx->cellScalarGlobIdx(fracCell); double dy = dataAccessObjectDy->cellScalarGlobIdx(fracCell); double dz = dataAccessObjectDz->cellScalarGlobIdx(fracCell); double NTG = dataAccessObjectDx->cellScalarGlobIdx(fracCell); cvf::Vec3d localX; cvf::Vec3d localY; cvf::Vec3d localZ; std::vector > 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) cvf::Mat4f invertedTransMatrix = transformMatrix().getInverted(); for (std::vector & planeCellPolygon : planeCellPolygons) { for (cvf::Vec3d& v : planeCellPolygon) { v.transformPoint(static_cast(invertedTransMatrix)); } } cvf::Vec3d localZinFracPlane; localZinFracPlane = localZ; 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; double skinfactor = 0.0; double transmissibility_X = 0.0; double transmissibility_Y = 0.0; double transmissibility_Z = 0.0; if (attachedFractureDefinition()) { std::vector fracPolygon = attachedFractureDefinition()->fracturePolygon(); std::vector fracPolygonDouble; for (auto v : fracPolygon) fracPolygonDouble.push_back(static_cast(v)); std::vector > polygonsDescribingFractureInCell; cvf::Vec3d areaVector; for (std::vector planeCellPolygon : planeCellPolygons) { std::vector >clippedPolygons = RigCellGeometryTools::clipPolygons(planeCellPolygon, fracPolygonDouble); for (std::vector clippedPolygon : clippedPolygons) { polygonsDescribingFractureInCell.push_back(clippedPolygon); } } double area; std::vector areaOfFractureParts; double length; std::vector lengthXareaOfFractureParts; for (std::vector fracturePartPolygon : polygonsDescribingFractureInCell) { areaVector = cvf::GeometryTools::polygonAreaNormal3D(fracturePartPolygon); area = areaVector.length(); areaOfFractureParts.push_back(area); length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, fracturePartPolygon); lengthXareaOfFractureParts.push_back(length * area); 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(localY))); Ay += abs(area*(fracturePlane.normal().dot(localX))); Az += abs(area*(fracturePlane.normal().dot(localZ))); //TODO: resulting values have only been checked for vertical fracture... } fractureArea = 0.0; for (double area : areaOfFractureParts) fractureArea += area; double totalAreaXLength = 0.0; for (double lengtXarea : lengthXareaOfFractureParts) totalAreaXLength += lengtXarea; fractureAreaWeightedlength = totalAreaXLength / fractureArea; double c = 0.008527; // TODO: Get value with units, is defined in RimReservoirCellResultsStorage double skinfactor = attachedFractureDefinition()->skinFactor; double slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D; transmissibility_X = 8 * c * ( permY * NTG ) * Ay / (dx + slDivPi); transmissibility_Y = 8 * c * ( permX * NTG ) * Ax / (dy + slDivPi); transmissibility_Z = 8 * c * permZ * Az / (dz + slDivPi); transmissibility = sqrt(transmissibility_X * transmissibility_X + transmissibility_Y * transmissibility_Y + transmissibility_Z * transmissibility_Z); } else { transmissibility = cvf::UNDEFINED_DOUBLE; } fracData.transmissibility = transmissibility; fracData.transmissibilities = cvf::Vec3d(transmissibility_X, transmissibility_Y, transmissibility_Z); fracData.totalArea = fractureArea; fracData.projectedAreas = cvf::Vec3d(Ax, Ay, Az); fracData.fractureLenght = fractureAreaWeightedlength; fracData.cellSizes = cvf::Vec3d(dx, dy, dz); fracData.permeabilities = cvf::Vec3d(permX, permY, permZ); fracData.NTG = NTG; fracData.skinFactor = skinfactor; fracData.cellIsActive = cellIsActive; //Since we loop over all potentially fractured cells, we only keep FractureData for cells where fracture have an non-zero area. if (fractureArea > 1e-5) { fracDataVec.push_back(fracData); } } m_rigFracture->setFractureData(fracDataVec); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- bool RimFracture::planeCellIntersectionPolygons(size_t cellindex, std::vector > & polygons, cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ) { cvf::Plane fracturePlane; cvf::Mat4f m = transformMatrix(); bool isCellIntersected = false; fracturePlane.setFromPointAndNormal(static_cast(m.translation()), static_cast(m.col(2))); RiaApplication* app = RiaApplication::instance(); RimView* activeView = RiaApplication::instance()->activeReservoirView(); if (!activeView) return isCellIntersected; RimEclipseView* activeRiv = dynamic_cast(activeView); if (!activeRiv) return isCellIntersected; const RigMainGrid* mainGrid = activeRiv->mainGrid(); if (!mainGrid) return isCellIntersected; RigCell cell = mainGrid->globalCellArray()[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(); const caf::SizeTArray8& cornerIndices = cell.cornerIndices(); hexCorners[0] = nodeCoords[cornerIndices[0]]; hexCorners[1] = nodeCoords[cornerIndices[1]]; hexCorners[2] = nodeCoords[cornerIndices[2]]; hexCorners[3] = nodeCoords[cornerIndices[3]]; hexCorners[4] = nodeCoords[cornerIndices[4]]; hexCorners[5] = nodeCoords[cornerIndices[5]]; hexCorners[6] = nodeCoords[cornerIndices[6]]; hexCorners[7] = nodeCoords[cornerIndices[7]]; //Find line-segments where cell and fracture plane intersects std::list > intersectionLineSegments; isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments); RigCellGeometryTools::createPolygonFromLineSegments(intersectionLineSegments, polygons); RigCellGeometryTools::findCellLocalXYZ(hexCorners, localX, localY, localZ); return isCellIntersected; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFracture::setRecomputeGeometryFlag() { m_recomputeGeometry = true; m_rivFracture->clearGeometryCache(); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- bool RimFracture::isRecomputeGeometryFlagSet() { return m_recomputeGeometry; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- cvf::Vec3d RimFracture::fracturePositionForUi() const { cvf::Vec3d v = m_anchorPosition; v.z() = -v.z(); return v; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- QString RimFracture::createOneBasedIJK() const { return QString("Cell : [%1, %2, %3]").arg(m_i + 1).arg(m_j + 1).arg(m_k + 1); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- QList RimFracture::calculateValueOptions(const caf::PdmFieldHandle* fieldNeedingOptions, bool * useOptionsOnly) { QList options; RimProject* proj = RiaApplication::instance()->project(); CVF_ASSERT(proj); RimOilField* oilField = proj->activeOilField(); if (oilField == nullptr) return options; if (fieldNeedingOptions == &m_fractureTemplate) { RimFractureDefinitionCollection* fracDefColl = oilField->fractureDefinitionCollection(); if (fracDefColl == nullptr) return options; for (RimEllipseFractureTemplate* fracDef : fracDefColl->fractureDefinitions()) { options.push_back(caf::PdmOptionItemInfo(fracDef->name(), fracDef)); } } return options; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFracture::defineUiOrdering(QString uiConfigName, caf::PdmUiOrdering& uiOrdering) { uiOrdering.add(&name); caf::PdmUiGroup* geometryGroup = uiOrdering.addNewGroup("Properties"); geometryGroup->add(&azimuth); geometryGroup->add(&perforationLength); geometryGroup->add(&m_fractureTemplate); caf::PdmUiGroup* fractureCenterGroup = uiOrdering.addNewGroup("Fracture Center Info"); fractureCenterGroup->add(&m_uiAnchorPosition); fractureCenterGroup->add(&m_displayIJK); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFracture::defineEditorAttribute(const caf::PdmFieldHandle* field, QString uiConfigName, caf::PdmUiEditorAttribute * attribute) { if (field == &azimuth) { caf::PdmUiDoubleSliderEditorAttribute* myAttr = dynamic_cast(attribute); if (myAttr) { myAttr->m_minimum = 0; myAttr->m_maximum = 360; } } } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFracture::setAnchorPosition(const cvf::Vec3d& pos) { m_anchorPosition = pos; setRecomputeGeometryFlag(); // Update ijk { std::vector cellindecies; RiaApplication* app = RiaApplication::instance(); RimView* activeView = RiaApplication::instance()->activeReservoirView(); if (!activeView) return; RimEclipseView* activeRiv = dynamic_cast(activeView); if (!activeRiv) return; const RigMainGrid* mainGrid = activeRiv->mainGrid(); if (!mainGrid) return; cvf::BoundingBox polygonBBox; polygonBBox.add(m_anchorPosition); mainGrid->findIntersectingCells(polygonBBox, &cellindecies); if (cellindecies.size() > 0) { size_t i = cvf::UNDEFINED_SIZE_T; size_t j = cvf::UNDEFINED_SIZE_T; size_t k = cvf::UNDEFINED_SIZE_T; size_t gridCellIndex = cellindecies[0]; if (mainGrid->ijkFromCellIndex(gridCellIndex, &i, &j, &k)) { m_i = static_cast(i); m_j = static_cast(j); m_k = static_cast(k); } } } } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const RigFracture* RimFracture::attachedRigFracture() const { CVF_ASSERT(m_rigFracture.notNull()); return m_rigFracture.p(); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RimFracture::setFractureTemplate(RimEllipseFractureTemplate* fractureTemplate) { m_fractureTemplate = fractureTemplate; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RimEllipseFractureTemplate* RimFracture::attachedFractureDefinition() const { return m_fractureTemplate(); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RivWellFracturePartMgr* RimFracture::fracturePartManager() { CVF_ASSERT(m_rivFracture.notNull()); return m_rivFracture.p(); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- bool RimFracture::hasValidGeometry() const { if (m_recomputeGeometry) return false; return (nodeCoords().size() > 0 && triangleIndices().size() > 0); }