#1320 - pre-proto - Prototyping calculation of upscaled parameter values from stimPlan. Bug in calculation so output values are wrong.

ApplicationCode/FileInterface/RifEclipseExportTools.cpp

@@ -83,7 +83,6 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName,  c
     }
 
     caf::ProgressInfo pi(fractures.size(), QString("Writing data to file %1").arg(fileName));
-
     RimEclipseWell* simWell = nullptr;
     RimWellPath* wellPath = nullptr;
 
@@ -99,6 +98,45 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName,  c
     if (caseUnit == RigEclipseCaseData::UNITS_FIELD) out << "-- Using field unit system" << "\n";
     out << "\n";
 
+
+    for (RimFracture* fracture : fractures) //For testing upscaling...
+    {
+        fracture->computeUpscaledPropertyFromStimPlan(caseToApply);
+        std::vector<RigFractureData> 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().left(15) + " ";
+
+
+            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 << qSetFieldWidth(10);
+            out << QString::number(fracData.upscaledStimPlanValue, 'f', 3);
+
+            out << "\n";
+        }
+    }
+
+
+
+
     printBackgroundDataHeaderLine(out);
 
 

ApplicationCode/ModelVisualization/RivWellFracturePartMgr.h

@@ -54,6 +54,8 @@ public:
     void appendGeometryPartsToModel(cvf::ModelBasicList* model, caf::DisplayCoordTransform* displayCoordTransform);
     void clearGeometryCache();
 
+    static std::vector<double>                 mirrorDataAtSingleDepth(std::vector<double> depthData);
+
 private:
     void updatePartGeometry(caf::DisplayCoordTransform* displayCoordTransform);
     void updatePartGeometryTexture(caf::DisplayCoordTransform* displayCoordTransform);
@@ -69,7 +71,6 @@ private:
     std::vector<cvf::Vec3f>             transfromToFractureDisplayCoords(std::vector<cvf::Vec3f> polygon, cvf::Mat4f m, caf::DisplayCoordTransform* displayCoordTransform);
     bool                                stimPlanCellTouchesPolygon(const std::vector<cvf::Vec3f>& polygon, double xMin, double xMax, double yMin, double yMax, float polygonXmin, float polygonXmax, float polygonYmin, float polygonYmax);
 
-    std::vector<double>                 mirrorDataAtSingleDepth(std::vector<double> depthData);
     static cvf::ref<cvf::DrawableGeo>   createGeo(const std::vector<cvf::uint>& triangleIndices, const std::vector<cvf::Vec3f>& nodeCoords);
      
 

ApplicationCode/ProjectDataModel/RimFracture.cpp

@@ -497,6 +497,153 @@ void RimFracture::computeTransmissibility(RimEclipseCase* caseToApply)
 }
 
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+void RimFracture::computeUpscaledPropertyFromStimPlan(RimEclipseCase* caseToApply)
+{
+    if (!attachedFractureDefinition()) return;
+    
+    RimStimPlanFractureTemplate* fracTemplateStimPlan;
+    if (dynamic_cast<RimStimPlanFractureTemplate*>(attachedFractureDefinition()))
+    {
+        fracTemplateStimPlan = dynamic_cast<RimStimPlanFractureTemplate*>(attachedFractureDefinition());
+    }
+    else return;
+
+    //Get correct unit system: 
+    RigEclipseCaseData::UnitsType caseUnit = caseToApply->eclipseCaseData()->unitsType();
+    RimDefines::UnitSystem unitForExport;
+
+    if (caseUnit == RigEclipseCaseData::UNITS_METRIC)
+    {
+        RiaLogging::debug(QString("Calculating transmissibilities for %1 in metric units").arg(name()));
+        unitForExport = RimDefines::UNITS_METRIC;
+    }
+    else if (caseUnit == RigEclipseCaseData::UNITS_FIELD)
+    {
+        RiaLogging::debug(QString("Calculating transmissibilities for %1 in field units").arg(name()));
+        unitForExport = RimDefines::UNITS_FIELD;
+    }
+    else
+    {
+        RiaLogging::error(QString("Unit system for case not supported for fracture export."));
+        return;
+    }
+
+
+    //TODO: Get these more generally: 
+    QString resultName = "CONDUCTIVITY";
+    QString resultUnit = "md-m";
+    size_t timeStepIndex = 0;
+
+
+    std::vector<std::vector<cvf::Vec3d> > stimPlanCellsAsPolygons;
+    std::vector<double> stimPlanParameterValues;
+    fracTemplateStimPlan->getStimPlanDataAsPolygonsAndValues(stimPlanCellsAsPolygons, stimPlanParameterValues, resultName, resultUnit, timeStepIndex);
+
+
+
+    //TODO: A lot of common code with function above... Can be cleaned up...?
+    std::vector<size_t> fracCells = getPotentiallyFracturedCells();
+ 
+
+    RigEclipseCaseData* eclipseCaseData = caseToApply->eclipseCaseData();
+
+    RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS;
+    RimReservoirCellResultsStorage* gridCellResults = caseToApply->results(porosityModel);
+    RigActiveCellInfo* activeCellInfo = eclipseCaseData->activeCellInfo(porosityModel);
+
+
+    std::vector<RigFractureData> fracDataVec;
+
+    for (size_t fracCell : fracCells)
+    {
+        bool cellIsActive = activeCellInfo->isActive(fracCell);
+
+        cvf::Vec3d localX;
+        cvf::Vec3d localY;
+        cvf::Vec3d localZ;
+        std::vector<std::vector<cvf::Vec3d> > 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/stimPlan mesh already is located)
+        cvf::Mat4f invertedTransMatrix = transformMatrix().getInverted();
+        for (std::vector<cvf::Vec3d> & planeCellPolygon : planeCellPolygons)
+        {
+            for (cvf::Vec3d& v : planeCellPolygon)
+            {
+                v.transformPoint(static_cast<cvf::Mat4d>(invertedTransMatrix));
+            }
+        }
+
+        cvf::Vec3d localZinFracPlane;
+        localZinFracPlane = localZ;
+        localZinFracPlane.transformVector(static_cast<cvf::Mat4d>(invertedTransMatrix));
+        cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO;
+        directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1));
+        directionOfLength.normalize();
+
+        RigFractureData fracData;
+        fracData.reservoirCellIndex = fracCell;
+
+
+
+        std::vector<cvf::Vec3f> fracPolygon = attachedFractureDefinition()->fracturePolygon(unitForExport);
+        std::vector<std::vector<cvf::Vec3d> > polygonsDescribingFractureInCell;
+
+        double area;
+        std::vector<double> areaOfFractureParts;
+        std::vector<double> areaXvalueOfFractureParts;
+
+
+        for (std::vector<cvf::Vec3d> planeCellPolygon : planeCellPolygons)
+        {
+
+            for (int i = 0; i < stimPlanParameterValues.size(); i++)
+            {
+                double stimPlanParameterValue = stimPlanParameterValues[i];
+                if (stimPlanParameterValue != 0)
+                {
+                    std::vector<cvf::Vec3d> stimPlanCell = stimPlanCellsAsPolygons[i];
+                    std::vector<std::vector<cvf::Vec3d> >clippedStimPlanPolygons = RigCellGeometryTools::clipPolygons(stimPlanCell, planeCellPolygon);
+                    if (clippedStimPlanPolygons.size() > 0)
+                    {
+                        for (auto clippedStimPlanPolygon : clippedStimPlanPolygons)
+                        {
+                            area = cvf::GeometryTools::polygonAreaNormal3D(clippedStimPlanPolygon).length();
+                            areaOfFractureParts.push_back(area);
+                            areaXvalueOfFractureParts.push_back(area*stimPlanParameterValue);
+                        }
+                    }
+                }
+            }
+
+
+        }
+
+        if (areaXvalueOfFractureParts.size() > 0)
+        {
+            double fractureCellArea = 0.0;
+            for (double area : areaOfFractureParts) fractureCellArea += area;
+            double fractureCellAreaXvalue = 0.0;
+            for (double areaXvalue : areaXvalueOfFractureParts) fractureCellAreaXvalue += area;
+            double upscaledValue = fractureCellAreaXvalue / fractureCellArea;
+
+            fracData.upscaledStimPlanValue = upscaledValue;
+            fracDataVec.push_back(fracData);
+        }
+
+        
+    }
+
+    m_rigFracture->setFractureData(fracDataVec);
+
+
+
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimFracture.h

@@ -85,6 +85,8 @@ public:
 
     std::vector<size_t>             getPotentiallyFracturedCells();
     void                            computeTransmissibility(RimEclipseCase* caseToApply);
+    void                            computeUpscaledPropertyFromStimPlan(RimEclipseCase* caseToApply);
+
 
     virtual void                    fieldChangedByUi(const caf::PdmFieldHandle* changedField, const QVariant& oldValue, const QVariant& newValue) override;
     cvf::Vec3d                      fracturePosition() const;

ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.cpp

@@ -39,6 +39,8 @@
 #include <QMessageBox>
 
 #include <algorithm>
+#include <vector>
+#include "RivWellFracturePartMgr.h"
 
 
 
@@ -311,7 +313,6 @@ void RimStimPlanFractureTemplate::loadDataAndUpdate()
     updateConnectedEditors();
 }
 
-
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
@@ -320,6 +321,24 @@ std::vector<std::vector<double>> RimStimPlanFractureTemplate::getDataAtTimeIndex
     return m_stimPlanFractureDefinitionData->getDataAtTimeIndex(resultName, unitName, timeStepIndex);
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+std::vector<std::vector<double>> RimStimPlanFractureTemplate::getMirroredDataAtTimeIndex(const QString& resultName, const QString& unitName, size_t timeStepIndex) const
+{
+    std::vector<std::vector<double>> notMirrordedData = m_stimPlanFractureDefinitionData->getDataAtTimeIndex(resultName, unitName, timeStepIndex);
+    std::vector<std::vector<double>> mirroredData;
+
+    for (std::vector<double> depthData : notMirrordedData)
+    {
+        std::vector<double> mirrordDepthData = RivWellFracturePartMgr::mirrorDataAtSingleDepth(depthData);
+        mirroredData.push_back(mirrordDepthData);
+    }
+
+
+    return mirroredData;
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
@@ -709,6 +728,46 @@ void RimStimPlanFractureTemplate::computeMinMax(const QString& resultName, const
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+void RimStimPlanFractureTemplate::getStimPlanDataAsPolygonsAndValues(std::vector<std::vector<cvf::Vec3d> > &cellsAsPolygons, std::vector<double> &parameterValues, const QString& resultName, const QString& unitName, size_t timeStepIndex)
+{
+
+    std::vector<std::vector<double>> propertyValuesAtTimeStep = getMirroredDataAtTimeIndex(resultName, unitName, timeStepIndex);
+
+    cellsAsPolygons.clear();
+    parameterValues.clear();
+
+    //TODO: Code partly copied from RivWellFracturePartMgr - can this be combined in some function?
+    std::vector<double> depthCoordsAtNodes = adjustedDepthCoordsAroundWellPathPosition();
+    std::vector<double> xCoordsAtNodes = getNegAndPosXcoords();
+
+    //Cells are around nodes instead of between nodes
+    std::vector<double> xCoords;
+    for (int i = 0; i < xCoordsAtNodes.size() - 1; i++) xCoords.push_back((xCoordsAtNodes[i] + xCoordsAtNodes[i + 1]) / 1);
+    std::vector<double> depthCoords;
+    for (int i = 0; i < depthCoordsAtNodes.size() - 1; i++) depthCoords.push_back((depthCoordsAtNodes[i] + depthCoordsAtNodes[i + 1]) / 1);
+
+    for (int i = 0; i < xCoords.size() - 1; i++)
+    {
+        for (int j = 0; j < depthCoords.size() - 1; j++)
+        {
+            std::vector<cvf::Vec3d> cellAsPolygon;
+            cellAsPolygon.push_back(cvf::Vec3d(static_cast<float>(xCoords[i]), static_cast<float>(depthCoords[j]), 0.0));
+            cellAsPolygon.push_back(cvf::Vec3d(static_cast<float>(xCoords[i + 1]), static_cast<float>(depthCoords[j]), 0.0));
+            cellAsPolygon.push_back(cvf::Vec3d(static_cast<float>(xCoords[i + 1]), static_cast<float>(depthCoords[j + 1]), 0.0));
+            cellAsPolygon.push_back(cvf::Vec3d(static_cast<float>(xCoords[i]), static_cast<float>(depthCoords[j + 1]), 0.0));
+            cellsAsPolygons.push_back(cellAsPolygon);
+            //TODO: Values for both neg and pos x values...
+            parameterValues.push_back(propertyValuesAtTimeStep[j+1][i+1]); //TODO test that this value exsist...
+
+        }
+    }
+
+
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimStimPlanFractureTemplate.h

@@ -62,7 +62,6 @@ public:
 
     void                                    fractureGeometry(std::vector<cvf::Vec3f>* nodeCoords, std::vector<cvf::uint>* triangleIndices, caf::AppEnum< RimDefines::UnitSystem > fractureUnit) override;
     std::vector<cvf::Vec3f>                 fracturePolygon(caf::AppEnum< RimDefines::UnitSystem > fractureUnit);
-
     void sortPolygon(std::vector<cvf::Vec3f> &polygon);
 
     std::vector<double>                     getNegAndPosXcoords();
@@ -70,12 +69,12 @@ public:
     std::vector<double>                     getStimPlanTimeValues();
     std::vector<std::pair<QString, QString> > getStimPlanPropertyNamesUnits() const;
     void                                    computeMinMax(const QString& resultName, const QString& unitName, double* minValue, double* maxValue) const;
+    void                                    getStimPlanDataAsPolygonsAndValues(std::vector<std::vector<cvf::Vec3d> > &cellsAsPolygons, std::vector<double> &parameterValue, const QString& resultName, const QString& unitName, size_t timeStepIndex);
 
     void                                    loadDataAndUpdate();
     void                                    setDefaultsBasedOnXMLfile();
-
     std::vector<std::vector<double>>        getDataAtTimeIndex(const QString& resultName, const QString& unitName, size_t timeStepIndex) const;
-
+    std::vector<std::vector<double>>        getMirroredDataAtTimeIndex(const QString& resultName, const QString& unitName, size_t timeStepIndex) const;
     
     virtual QList<caf::PdmOptionItemInfo>   calculateValueOptions(const caf::PdmFieldHandle* fieldNeedingOptions, bool* useOptionsOnly) override;
 
@@ -88,21 +87,19 @@ private:
     void                                    updateUiTreeName();
 
     void                                    readStimPlanXMLFile(QString * errorMessage);
-
-
     size_t                                  readStimplanGridAndTimesteps(QXmlStreamReader &xmlStream);
-
     static double                           getAttributeValueDouble(QXmlStreamReader &xmlStream, QString parameterName);
     static QString                          getAttributeValueString(QXmlStreamReader &xmlStream, QString parameterName);
     void                                    getGriddingValues(QXmlStreamReader &xmlStream, std::vector<double>& gridValues, size_t& startNegValues);
     std::vector<std::vector<double>>        getAllDepthDataAtTimeStep(QXmlStreamReader &xmlStream, size_t startingNegValuesXs);
 
-    caf::PdmField<QString>                      m_stimPlanFileName;
-    cvf::ref<RigStimPlanFractureDefinition>     m_stimPlanFractureDefinitionData;
-
     bool                                    numberOfParameterValuesOK(std::vector<std::vector<double>> propertyValuesAtTimestep);
     bool                                    setPropertyForPolygonDefault();
     void                                    setDepthOfWellPathAtFracture();
     QString                                 getUnitForStimPlanParameter(QString parameterName);
 
+    caf::PdmField<QString>                      m_stimPlanFileName;
+    cvf::ref<RigStimPlanFractureDefinition>     m_stimPlanFractureDefinitionData;
+
+
 };

ApplicationCode/ReservoirDataModel/RigFracture.h

@@ -46,6 +46,7 @@ public:
 
     bool cellIsActive;
     
+    double upscaledStimPlanValue;
 
 
 };