#1507 Setting up fractureGrid for ellipse fracture template, and adapting transmissibility calculation to be able to handle both ellipse and stimplan fractures. Missing setup of fractureCells for ellipse when loading old project.

ApplicationCode/FileInterface/RifFractureExportTools.cpp

@@ -224,10 +224,8 @@ void RifFractureExportTools::exportWellPathFracturesToEclipseDataInputFile(const
     {
         using CellIdxSpace = RigTransmissibilityCondenser::CellAddress;
 
-        RimStimPlanFractureTemplate* fracTemplateStimPlan = dynamic_cast<RimStimPlanFractureTemplate*>(fracture->attachedFractureDefinition());
+        RimFractureTemplate* fracTemplate = fracture->attachedFractureDefinition();
-        const RigFractureGrid* fractureGrid = fracTemplateStimPlan->fractureGrid();
+        const RigFractureGrid* fractureGrid = fracTemplate->fractureGrid();
-
-        if (!fracTemplateStimPlan) continue; // We do not handle Elliptical fractures yet
 
         RigTransmissibilityCondenser transCondenser;
 
@@ -330,24 +328,24 @@ void RifFractureExportTools::exportWellPathFracturesToEclipseDataInputFile(const
             if (intersection.endpointCount)
             {
                 radialTrans = RigFractureTransmissibilityEquations::fractureCellToWellRadialTrans(fractureWellCell.getConductivtyValue(),
-                                                                                                                       fractureWellCell.cellSizeX(),
+                                                                                                                   fractureWellCell.cellSizeX(),
-                                                                                                                       fractureWellCell.cellSizeZ(),
+                                                                                                                   fractureWellCell.cellSizeZ(),
-                                                                                                                       fracture->wellRadius(),
+                                                                                                                   fracture->wellRadius(),
-                                                                                                                       fracTemplateStimPlan->skinFactor(),
+                                                                                                                   fracTemplate->skinFactor(),
-                                                                                                                       cDarcyInCorrectUnit);
+                                                                                                                   cDarcyInCorrectUnit);
             }
 
             double linearTrans = 0.0;
             if (intersection.hlength > 0.0 || intersection.vlength > 0.0 )
             {
                 linearTrans = RigFractureTransmissibilityEquations::fractureCellToWellLinearTrans(fractureWellCell.getConductivtyValue(),
-                                                                                                                      fractureWellCell.cellSizeX(),
+                                                                                                                   fractureWellCell.cellSizeX(),
-                                                                                                                      fractureWellCell.cellSizeZ(),
+                                                                                                                   fractureWellCell.cellSizeZ(),
-                                                                                                                      intersection.vlength,
+                                                                                                                   intersection.vlength,
-                                                                                                                      intersection.hlength ,
+                                                                                                                   intersection.hlength ,
-                                                                                                                      fracture->perforationEfficiency, 
+                                                                                                                   fracture->perforationEfficiency, 
-                                                                                                                      fracTemplateStimPlan->skinFactor(), 
+                                                                                                                   fracTemplate->skinFactor(),
-                                                                                                                      cDarcyInCorrectUnit);
+                                                                                                                   cDarcyInCorrectUnit);
             }
             
             double totalWellTrans = 0.5 * intersection.endpointCount * radialTrans + linearTrans;
@@ -371,9 +369,9 @@ void RifFractureExportTools::exportWellPathFracturesToEclipseDataInputFile(const
             }
         }
         out << "\n" << "\n" << "\n----------- All Transimissibilities " << fracture->name() << " -------------------- \n\n";
-        out << QString::fromStdString(transCondenser.neighborTransDebugOutput(mainGrid, fracTemplateStimPlan));
+        out << QString::fromStdString(transCondenser.neighborTransDebugOutput(mainGrid, fractureGrid));
         out << "\n" << "\n" << "\n----------- Condensed Results -------------------- \n\n";
-        out << QString::fromStdString(transCondenser.condensedTransDebugOutput(mainGrid, fracTemplateStimPlan));
+        out << QString::fromStdString(transCondenser.condensedTransDebugOutput(mainGrid, fractureGrid));
         out << "\n" ;
     } 
 

ApplicationCode/ProjectDataModel/RimEllipseFractureTemplate.cpp

@@ -20,6 +20,9 @@
 
 #include "RiaLogging.h"
 
+#include "RigCellGeometryTools.h"
+#include "RigFractureCell.h"
+#include "RigFractureGrid.h"
 #include "RigTesselatorTools.h"
 
 #include "RimDefines.h"
@@ -30,6 +33,7 @@
 #include "cafPdmObject.h"
 
 #include "cvfVector3.h"
+#include "cvfGeometryTools.h"
 
 
 
@@ -47,6 +51,8 @@ RimEllipseFractureTemplate::RimEllipseFractureTemplate(void)
     CAF_PDM_InitField(&width,       "Width",            1.0f,    "Width", "", "", "");
 
     CAF_PDM_InitField(&permeability,"Permeability",     22000.f, "Permeability [mD]", "", "", "");
+
+    m_fractureGrid = new RigFractureGrid();
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -84,8 +90,13 @@ void RimEllipseFractureTemplate::fieldChangedByUi(const caf::PdmFieldHandle* cha
             }
 
             proj->createDisplayModelAndRedrawAllViews();
+            setupFractureGridCells();
         }
     }
+    if (changedField == &width || changedField == &permeability)
+    {
+        setupFractureGridCells();
+    }
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -173,6 +184,105 @@ void RimEllipseFractureTemplate::changeUnits()
     this->updateConnectedEditors();
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+void RimEllipseFractureTemplate::setupFractureGridCells()
+{
+
+    std::vector<RigFractureCell> fractureCells;
+    std::pair<size_t, size_t> wellCenterFractureCellIJ = std::make_pair(0, 0);
+
+    bool wellCenterFractureCellFound = false;
+
+    int numberOfCellsX = 5;
+    int numberOfCellsY = 5;
+    
+    double cellSizeX = (halfLength * 2) / numberOfCellsX;
+    double cellSizeZ = height / numberOfCellsY;
+
+    double cellArea = cellSizeX * cellSizeZ;
+    double areaTresholdForIncludingCell = 0.5 * cellArea;
+
+
+    for (int i = 0; i < numberOfCellsX; i++)
+    {
+        for (int j = 0; j < numberOfCellsX; j++)
+        {
+            double X1 = - halfLength +  i    * cellSizeX;
+            double X2 = - halfLength + (i+1) * cellSizeX;
+            double Y1 = - height / 2 +  j    * cellSizeZ;
+            double Y2 = - height / 2 + (j+1) * cellSizeZ;
+
+            std::vector<cvf::Vec3d> cellPolygon;
+            cellPolygon.push_back(cvf::Vec3d(X1, Y1, 0.0));
+            cellPolygon.push_back(cvf::Vec3d(X2, Y1, 0.0));
+            cellPolygon.push_back(cvf::Vec3d(X2, Y2, 0.0));
+            cellPolygon.push_back(cvf::Vec3d(X1, Y2, 0.0));
+            
+            double cond = cvf::UNDEFINED_DOUBLE;
+            if (fractureTemplateUnit == RimDefines::UNITS_METRIC)
+            {
+                //Conductivity should be md-m, width is in m
+                cond = permeability * width;
+            }
+            else if(fractureTemplateUnit == RimDefines::UNITS_FIELD)
+            {
+                //Conductivity should be md-ft, but width is in inches 
+                cond = permeability * RimDefines::inchToFeet(width);
+            }
+
+            std::vector<cvf::Vec3f> ellipseFracPolygon = fracturePolygon(fractureTemplateUnit());
+            std::vector<cvf::Vec3d> ellipseFracPolygonDouble;
+            for (auto v : ellipseFracPolygon) ellipseFracPolygonDouble.push_back(static_cast<cvf::Vec3d>(v));
+            std::vector<std::vector<cvf::Vec3d> >clippedFracturePolygons = RigCellGeometryTools::intersectPolygons(cellPolygon, ellipseFracPolygonDouble);
+            if (clippedFracturePolygons.size() > 0)
+            {
+                for (auto clippedFracturePolygon : clippedFracturePolygons)
+                {
+                    double areaCutPolygon = cvf::GeometryTools::polygonAreaNormal3D(clippedFracturePolygon).length();
+                    if (areaCutPolygon < areaTresholdForIncludingCell) cond = 0.0; //Cell is excluded from calculation, cond is set to zero. Must be included for indexing to be correct
+                }
+            }
+            else cond = 0.0;
+
+            RigFractureCell fractureCell(cellPolygon, i, j);
+            fractureCell.setConductivityValue(cond);
+
+            if (cellPolygon[0].x() < 0.0 && cellPolygon[1].x() > 0.0)
+            {
+                if (cellPolygon[1].y() > 0.0 && cellPolygon[2].y() < 0.0)
+                {
+                    wellCenterFractureCellIJ = std::make_pair(fractureCell.getI(), fractureCell.getJ());
+                    RiaLogging::debug(QString("Setting wellCenterStimPlanCell at cell %1, %2").
+                                      arg(QString::number(fractureCell.getI()), QString::number(fractureCell.getJ())));
+
+                    wellCenterFractureCellFound = true;
+                }
+            }
+
+            fractureCells.push_back(fractureCell);
+
+        }
+    }
+
+    
+    m_fractureGrid->setFractureCells(fractureCells);
+    m_fractureGrid->setWellCenterFractureCellIJ(wellCenterFractureCellIJ);
+    m_fractureGrid->setICellCount(numberOfCellsX);
+    m_fractureGrid->setJCellCount(numberOfCellsY);
+
+
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+const RigFractureGrid* RimEllipseFractureTemplate::fractureGrid() const
+{
+    return m_fractureGrid.p();
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimEllipseFractureTemplate.h

@@ -27,10 +27,12 @@
 #include "cafPdmFieldCvfVec3d.h"
 
 #include "cvfBase.h"
+#include "cvfObject.h"
 #include "cvfVector3.h"
 
 #include <vector>
 
+class RigFractureGrid;
 
 //==================================================================================================
 ///  
@@ -56,7 +58,13 @@ public:
     std::vector<cvf::Vec3f>         fracturePolygon(RimDefines::UnitSystem  fractureTemplateUnit);
     void                            changeUnits();
     
+    const RigFractureGrid*          fractureGrid() const;
+
 protected:
     virtual void defineUiOrdering(QString uiConfigName, caf::PdmUiOrdering& uiOrdering);
 
+private:
+    void                             setupFractureGridCells();
+    cvf::ref<RigFractureGrid>        m_fractureGrid;
+
 };

ApplicationCode/ProjectDataModel/RimFractureTemplate.h

@@ -32,6 +32,7 @@
 #include <vector>
 
 class RigEclipseCaseData;
+class RigFractureGrid;
 
 //==================================================================================================
 ///  
@@ -75,6 +76,9 @@ public:
     
     virtual void                    fractureGeometry(std::vector<cvf::Vec3f>* nodeCoords, std::vector<cvf::uint>* triangleIndices, RimDefines::UnitSystem fractureTemplateUnit) = 0;
     virtual std::vector<cvf::Vec3f> fracturePolygon(RimDefines::UnitSystem fractureTemplateUnit) = 0;
+
+    virtual const RigFractureGrid*  fractureGrid() const = 0;
+
 protected:
     virtual void                    defineUiOrdering(QString uiConfigName, caf::PdmUiOrdering& uiOrdering);
 

ApplicationCode/ReservoirDataModel/RigTransmissibilityCondenser.cpp

@@ -181,7 +181,7 @@ void RigTransmissibilityCondenser::calculateCondensedTransmissibilitiesIfNeeded(
 
 void printCellAddress(std::stringstream& str, 
                       const RigMainGrid* mainGrid, 
-                      const RimStimPlanFractureTemplate* fractureGrid,
+                      const RigFractureGrid* fractureGrid,
                       RigTransmissibilityCondenser::CellAddress cellAddr)
 {
     using CellAddress =  RigTransmissibilityCondenser::CellAddress;
@@ -200,7 +200,7 @@ void printCellAddress(std::stringstream& str,
         case CellAddress::STIMPLAN: 
         {
             str << "STP ";
-            const RigFractureCell& stpCell = fractureGrid->fractureGrid()->cellFromIndex(cellAddr.m_globalCellIdx);
+            const RigFractureCell& stpCell = fractureGrid->cellFromIndex(cellAddr.m_globalCellIdx);
             str << std::setw(5) << stpCell.getI()+1 << std::setw(5) << stpCell.getJ()+1  << std::setw(5) << " ";
         }
         break;
@@ -220,7 +220,7 @@ void printCellAddress(std::stringstream& str,
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-std::string RigTransmissibilityCondenser::neighborTransDebugOutput(const RigMainGrid* mainGrid, const RimStimPlanFractureTemplate* fractureGrid)
+std::string RigTransmissibilityCondenser::neighborTransDebugOutput(const RigMainGrid* mainGrid, const RigFractureGrid* fractureGrid)
 {
     std::stringstream debugText;
     for ( const auto& adrEqIdxPair : m_neighborTransmissibilities )
@@ -243,7 +243,7 @@ std::string RigTransmissibilityCondenser::neighborTransDebugOutput(const RigMain
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-std::string RigTransmissibilityCondenser::condensedTransDebugOutput(const RigMainGrid* mainGrid, const RimStimPlanFractureTemplate* fractureGrid)
+std::string RigTransmissibilityCondenser::condensedTransDebugOutput(const RigMainGrid* mainGrid, const RigFractureGrid* fractureGrid)
 {
     std::stringstream debugText;
     for ( const auto& adrEqIdxPair : m_condensedTransmissibilities )

ApplicationCode/ReservoirDataModel/RigTransmissibilityCondenser.h

@@ -26,6 +26,7 @@
 
 class RigMainGrid;
 class RimStimPlanFractureTemplate;
+class RigFractureGrid;
 
 class RigTransmissibilityCondenser
 {
@@ -74,8 +75,8 @@ public:
 
     double condensedTransmissibility( CellAddress externalCell1, CellAddress externalCell2);
 
-    std::string neighborTransDebugOutput(const RigMainGrid* mainGrid, const RimStimPlanFractureTemplate* fractureGrid);
+    std::string neighborTransDebugOutput(const RigMainGrid* mainGrid, const RigFractureGrid* fractureGrid);
-    std::string condensedTransDebugOutput(const RigMainGrid* mainGrid, const RimStimPlanFractureTemplate* fractureGrid);
+    std::string condensedTransDebugOutput(const RigMainGrid* mainGrid, const RigFractureGrid* fractureGrid);
 
 private:
     void calculateCondensedTransmissibilitiesIfNeeded();

ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.cpp

@@ -20,18 +20,15 @@ RigWellPathStimplanIntersector::RigWellPathStimplanIntersector(const RigWellPath
     std::vector<cvf::Vec3f> fracturePolygonf ; 
     std::vector<std::vector<cvf::Vec3d> > stpCellPolygons;
     {
-        auto stimPlanFractureTemplate = dynamic_cast<RimStimPlanFractureTemplate*> (rimFracture->attachedFractureDefinition());
+        RimFractureTemplate* fractureTemplate = rimFracture->attachedFractureDefinition();
 
-        CVF_ASSERT(stimPlanFractureTemplate);
-
-        fracturePolygonf = stimPlanFractureTemplate->fracturePolygon(rimFracture->fractureUnit());
         {
-            const std::vector<RigFractureCell>& stpCells = stimPlanFractureTemplate->fractureGrid()->fractureCells();
+            const std::vector<RigFractureCell>& stpCells = fractureTemplate->fractureGrid()->fractureCells();
             for ( const auto& stpCell: stpCells ) stpCellPolygons.push_back(stpCell.getPolygon());
         }
     }
 
-    calculate(fractureXf, fracturePolygonf, wellPathPoints, wellRadius, stpCellPolygons, m_stimPlanCellIdxToIntersectionInfoMap);
+    calculate(fractureXf, wellPathPoints, wellRadius, stpCellPolygons, m_stimPlanCellIdxToIntersectionInfoMap);
 
 }
 
@@ -40,7 +37,6 @@ RigWellPathStimplanIntersector::RigWellPathStimplanIntersector(const RigWellPath
 //--------------------------------------------------------------------------------------------------
 
 void RigWellPathStimplanIntersector::calculate(const cvf::Mat4f &fractureXf, 
-                                               const std::vector<cvf::Vec3f>& fracturePolygonf, 
                                                const std::vector<cvf::Vec3d>& wellPathPointsOrg, 
                                                double wellRadius, 
                                                const std::vector<std::vector<cvf::Vec3d> >& stpCellPolygons, 
@@ -48,8 +44,22 @@ void RigWellPathStimplanIntersector::calculate(const cvf::Mat4f &fractureXf,
 {
     cvf::Mat4d toFractureXf = cvf::Mat4d(fractureXf.getInverted());
 
-    std::vector<cvf::Vec3d> fracturePolygon;
+    //Find bounding box
-    for ( auto fpv: fracturePolygonf ) fracturePolygon.push_back(cvf::Vec3d(fpv));
+    cvf::BoundingBox polygonBBox;
+    for (std::vector<cvf::Vec3d> fracCellPolygon : stpCellPolygons)
+    {
+        for (cvf::Vec3d nodeCoord : fracCellPolygon) polygonBBox.add(nodeCoord);
+    }
+    cvf::Vec3d bboxCorners[8];
+    polygonBBox.cornerVertices(bboxCorners);
+
+    //put BB corners into polygon, to use clipPolylineByPolygon function below
+    //Since we are in 2D, the 4 corners of the boundingbox are repeated so only 4 are included in the polygon
+    std::vector<cvf::Vec3d> boundingBoxPolygon;
+    for (int i = 0; i < 4; i++)
+    {
+        boundingBoxPolygon.push_back(bboxCorners[i]);
+    }
 
     // Convert well path to fracture template system
 
@@ -59,7 +69,7 @@ void RigWellPathStimplanIntersector::calculate(const cvf::Mat4f &fractureXf,
     // Clip well path to fracture domain 
 
     std::vector<std::vector<cvf::Vec3d> > wellPathPartsWithinFracture =
-        RigCellGeometryTools::clipPolylineByPolygon(fractureRelativeWellPathPoints, fracturePolygon, RigCellGeometryTools::INTERPOLATE_LINE_Z);
+        RigCellGeometryTools::clipPolylineByPolygon(fractureRelativeWellPathPoints, boundingBoxPolygon, RigCellGeometryTools::INTERPOLATE_LINE_Z);
 
     // Remove the part of the well path that is more than well radius away from the fracture plane
 

ApplicationCode/ReservoirDataModel/RigWellPathStimplanIntersector.h

@@ -45,7 +45,6 @@ public:
 private:
     friend class RigWellPathStimplanIntersectorTester;
     static void calculate(const cvf::Mat4f& fractureXf,
-                          const std::vector<cvf::Vec3f>& fracturePolygon,
                           const std::vector<cvf::Vec3d>& wellPathPoints,
                           double wellRadius,
                           const std::vector<std::vector<cvf::Vec3d> >& stpCellPolygons,
@@ -65,7 +64,7 @@ public:
                               const std::vector<std::vector<cvf::Vec3d> >& stpCellPolygons,
                               std::map<size_t, RigWellPathStimplanIntersector::WellCellIntersection>& stimPlanCellIdxToIntersectionInfoMap)
     {
-        RigWellPathStimplanIntersector::calculate(fractureXf, fracturePolygon, wellPathPoints, wellRadius, stpCellPolygons, stimPlanCellIdxToIntersectionInfoMap);
+        RigWellPathStimplanIntersector::calculate(fractureXf, wellPathPoints, wellRadius, stpCellPolygons, stimPlanCellIdxToIntersectionInfoMap);
     }
 
 };