riTRANXYZ: Added to build system. First implementation.

ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake

@@ -17,6 +17,7 @@ ${CEE_CURRENT_LIST_DIR}RigActiveCellsResultAccessor.h
 ${CEE_CURRENT_LIST_DIR}RigCellEdgeResultAccessor.h
 ${CEE_CURRENT_LIST_DIR}RigCombTransResultAccessor.h
 ${CEE_CURRENT_LIST_DIR}RigCombMultResultAccessor.h
+${CEE_CURRENT_LIST_DIR}RigCombRiTransResultAccessor.h
 ${CEE_CURRENT_LIST_DIR}RigResultModifier.h
 ${CEE_CURRENT_LIST_DIR}RigResultModifierFactory.h
 ${CEE_CURRENT_LIST_DIR}RigLocalGrid.h
@@ -50,6 +51,7 @@ ${CEE_CURRENT_LIST_DIR}RigActiveCellsResultAccessor.cpp
 ${CEE_CURRENT_LIST_DIR}RigCellEdgeResultAccessor.cpp
 ${CEE_CURRENT_LIST_DIR}RigCombTransResultAccessor.cpp
 ${CEE_CURRENT_LIST_DIR}RigCombMultResultAccessor.cpp
+${CEE_CURRENT_LIST_DIR}RigCombRiTransResultAccessor.cpp
 ${CEE_CURRENT_LIST_DIR}RigResultModifierFactory.cpp
 ${CEE_CURRENT_LIST_DIR}RigLocalGrid.cpp
 ${CEE_CURRENT_LIST_DIR}RigMainGrid.cpp

ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.cpp

@@ -18,10 +18,12 @@
 
 #include "RigCombRiTransResultAccessor.h"
 
+#include "RigMainGrid.h"
 #include "RigGridBase.h"
+#include "RigCell.h"
 
 #include <cmath>
-
+#include "cvfGeometryTools.h"
 
 //--------------------------------------------------------------------------------------------------
 /// 
@@ -65,7 +67,7 @@ double RigCombRiTransResultAccessor::cellScalar(size_t gridLocalCellIndex) const
     return HUGE_VAL;
 }
 
-double RigCombRiTransResultAccessor::getPermValue(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId)
+double RigCombRiTransResultAccessor::getPermValue(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId) const
 {
     switch (faceId)
     {
@@ -84,19 +86,19 @@ double RigCombRiTransResultAccessor::getPermValue(size_t gridLocalCellIndex, cvf
     case cvf::StructGridInterface::POS_K:
     case cvf::StructGridInterface::NEG_K:
         {
-            return zPermAccessor->cellScalar(gridLocalCellIndex);
+            return m_zPermAccessor->cellScalar(gridLocalCellIndex);
         }
         break;
     }
     return HUGE_VAL;
 }
 
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-double RigCombRiTransResultAccessor::calculateHalfCellTrans( size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId)
+double RigCombRiTransResultAccessor::getNtgValue(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId ) const
 {
-    double perm = getPermValue(gridLocalCellIndex, faceId);
     double ntg = 1.0;
 
     if (faceId != cvf::StructGridInterface::POS_K && faceId != cvf::StructGridInterface::NEG_K)
@@ -104,40 +106,141 @@ double RigCombRiTransResultAccessor::calculateHalfCellTrans( size_t gridLocalCel
         m_ntgAccessor->cellScalar(gridLocalCellIndex);
     }
 
-    RigCell& cell = m_grid->cell(gridLocalCellIndex);
-    
-    cvf::Vec3d centerToFace = cell.faceCenter(faceId) - cell.center();
-    cvf::Vec3d faceAreaVec = cell.faceNormalWithAreaLenght(faceId);
-
-    double halfCellTrans = perm*ntg*(faceAreaVec*centerToFace) / (centerToFace*centerToFace);
-
-    return halfCellTrans;
+    return ntg;
 }
 
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-double RigCombRiTransResultAccessor::cellFaceScalar(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId) const
+double halfCellTransmissibility( double perm , double ntg, const cvf::Vec3d& centerToFace, const cvf::Vec3d& faceAreaVec) 
 {
-    size_t reservoirCellIndex = m_grid->reservoirCellIndex(gridLocalCellIndex);
-    RigFault* fault = m_grid->mainGrid()->findFaultFromCellIndexAndCellFace(reservoirCellIndex, faceId);
+    return perm*ntg*(faceAreaVec*centerToFace) / (centerToFace*centerToFace);
+}
 
-    if (fault)
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+double newtran(double cdarchy, double mult, double halfCellTrans, double neighborHalfCellTrans)
+{
+    return cdarchy * mult / ( ( 1 / halfCellTrans) + (1 / neighborHalfCellTrans) );
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+double RigCombRiTransResultAccessor::calculateHalfCellTrans(size_t gridLocalCellIndex, size_t neighborGridCellIdx, cvf::StructGridInterface::FaceType faceId, bool isFaultFace) const
+{
+    const RigCell& cell = m_grid->cell(gridLocalCellIndex);
+
+  
+    cvf::Vec3d faceAreaVec;
+    cvf::Vec3d faceCenter;
+
+    if (isFaultFace)
     {
         
+        calculateConnectionGeometry( gridLocalCellIndex, neighborGridCellIdx, faceId, &faceCenter, &faceAreaVec);
     }
     else
     {
+        const RigCell& cell = m_grid->cell(gridLocalCellIndex);
+        faceCenter = cell.faceCenter(faceId);
+        
+        faceAreaVec = cell.faceNormalWithAreaLenght(faceId);
+    }
+
+    cvf::Vec3d centerToFace = faceCenter - cell.center();
+
+    double perm = getPermValue(gridLocalCellIndex, faceId);
+    double ntg = getNtgValue(gridLocalCellIndex, faceId );
+
+    return halfCellTransmissibility(perm, ntg, centerToFace, faceAreaVec);
+}
+
+void RigCombRiTransResultAccessor::calculateConnectionGeometry( size_t gridLocalCellIndex, size_t neighborGridCellIdx, 
+                                                                cvf::StructGridInterface::FaceType faceId,
+                                                                cvf::Vec3d* faceCenter, cvf::Vec3d* faceAreaVec) const
+{
+    CVF_TIGHT_ASSERT(faceCenter && faceAreaVec);
+
+    *faceCenter = cvf::Vec3d::ZERO;
+    *faceAreaVec = cvf::Vec3d::ZERO; 
+    const RigMainGrid* mainGrid = m_grid->mainGrid();
+
+    const RigCell& c1 = m_grid->cell(gridLocalCellIndex);
+    const RigCell& c2 = m_grid->cell(neighborGridCellIdx);
+
+    std::vector<size_t> polygon;
+    std::vector<cvf::Vec3d> intersections;
+    caf::SizeTArray4 face1;
+    caf::SizeTArray4 face2;
+    c1.faceIndices(faceId, &face1);
+    c2.faceIndices(cvf::StructGridInterface::oppositeFace(faceId), &face2);
+
+    bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads(
+        &polygon, 
+        &intersections, 
+        (cvf::EdgeIntersectStorage<size_t>*)NULL, 
+        cvf::wrapArrayConst(&(mainGrid->nodes())), 
+        face1.data(), 
+        face2.data(), 
+        1e-6);
+
+
+    if (foundOverlap)
+    {
+        std::vector<cvf::Vec3d> realPolygon;
+
+        for (size_t pIdx = 0; pIdx < polygon.size(); ++pIdx)
+        {
+            if (polygon[pIdx] < m_grid->mainGrid()->nodes().size())
+                realPolygon.push_back(m_grid->mainGrid()->nodes()[polygon[pIdx]]);
+            else
+                realPolygon.push_back(intersections[polygon[pIdx] - m_grid->mainGrid()->nodes().size()]);
+        }
+
+        // Polygon center
+        for (size_t pIdx = 0; pIdx < realPolygon.size(); ++pIdx)
+        {
+            *faceCenter += realPolygon[pIdx];
+        }
+
+        *faceCenter *= 1.0/realPolygon.size();
+
+        // Polygon area vector
+
+        *faceAreaVec = cvf::GeometryTools::polygonAreaNormal3D(realPolygon);
+
+    }
+
+}
+//--------------------------------------------------------------------------------------------------
+/// 
+/// Neighbour cell transmisibilities only is calculated here
+/// Not NNC transmisibilities. Thart has to be done separatly elsewhere 
+///
+/// Todo: What about Grid to Grid connections ?
+/// Todo: needs optimization. Things are done several times. Caching of the results should be considered. etc.
+///
+//--------------------------------------------------------------------------------------------------
+double RigCombRiTransResultAccessor::cellFaceScalar(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId) const
+{
     size_t i, j, k, neighborGridCellIdx;
     m_grid->ijkFromCellIndex(gridLocalCellIndex, &i, &j, &k);
 
     if (m_grid->cellIJKNeighbor(i, j, k, faceId, &neighborGridCellIdx))
     {
-            double halfCellTrans = calculateHalfCellTrans(gridLocalCellIndex, faceId);
-            double neighborHalfCellTrans = calculateHalfCellTrans(neighborGridCellIdx, cvf::StructGridInterface::oppositeFace(faceId));
+        size_t reservoirCellIndex = m_grid->reservoirCellIndex(gridLocalCellIndex);
+        const RigFault* fault = m_grid->mainGrid()->findFaultFromCellIndexAndCellFace(reservoirCellIndex, faceId);
+        bool isOnFault = fault;
 
-            return m_cdarchy/ ( ( 1 / halfCellTrans) + (1 / neighborHalfCellTrans) );
-        }
+        double halfCellTrans = 0;
+        double neighborHalfCellTrans = 0;
+
+        halfCellTrans = calculateHalfCellTrans(gridLocalCellIndex, neighborGridCellIdx, faceId, isOnFault);
+        neighborHalfCellTrans = calculateHalfCellTrans(neighborGridCellIdx, gridLocalCellIndex, cvf::StructGridInterface::oppositeFace(faceId), isOnFault);
+
+        return newtran(m_cdarchy, 1.0, halfCellTrans, neighborHalfCellTrans);
     }
 
     return HUGE_VAL;

ApplicationCode/ReservoirDataModel/RigCombRiTransResultAccessor.h

@@ -49,6 +49,12 @@ public:
     virtual double  cellFaceScalar(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId) const;
 
 private:
+    double getPermValue(size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId) const;
+    double calculateHalfCellTrans( size_t gridLocalCellIndex, size_t neighborGridCellIdx, cvf::StructGridInterface::FaceType faceId, bool isFaultFace) const;
+
+    double getNtgValue( size_t gridLocalCellIndex, cvf::StructGridInterface::FaceType faceId ) const;
+    void   calculateConnectionGeometry( size_t gridLocalCellIndex, size_t neighborGridCellIdx, cvf::StructGridInterface::FaceType faceId, cvf::Vec3d* centerToFace, cvf::Vec3d* faceAreaVec) const;
+
 
     cvf::ref<RigResultAccessor> m_xPermAccessor;
     cvf::ref<RigResultAccessor> m_yPermAccessor;

ApplicationCode/ReservoirDataModel/cvfGeometryTools.cpp

@@ -517,6 +517,60 @@ void GeometryTools::addMidEdgeNodes(std::list<std::pair<cvf::uint, bool> >* poly
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+/// Based on http://geomalgorithms.com/a01-_area.html
+/// This method returns the polygon normal with length equal to the polygon area.
+/// The components of the normal is thus the size of projected area into each of the main axis planes
+//--------------------------------------------------------------------------------------------------
+cvf::Vec3d GeometryTools::polygonAreaNormal3D(const std::vector<cvf::Vec3d>& polygon)
+{
+    size_t pSize = polygon.size();
+    switch (pSize)
+    {
+    case 0:
+    case 1:
+    case 2:
+        {
+            return cvf::Vec3d::ZERO;
+        }
+        break;
+    case 3:
+        {
+            return 0.5 * ((polygon[1]-polygon[0]) ^ (polygon[2] - polygon[0]));
+        }
+        break;
+    case 4:
+        {
+            // Cross product of diagonal = 2*A
+            return 0.5* ((polygon[2]-polygon[0]) ^ (polygon[3] - polygon[1]));
+        }
+        break;
+    default:
+        {
+            /// JJS:
+            // This is possibly not the fastest approach with large polygons, where a scaled projections approach would be better, 
+            // but I suspect this (simpler) approach is faster for small polygons, as long as we do not have the polygon normal up front.
+            //
+            cvf::Vec3d areaNormal(cvf::Vec3d::ZERO);
+            size_t h = (pSize - 1)/2;
+            size_t k = (pSize % 2) ? 0 : pSize - 1;
+
+            // First quads
+            for (size_t i = 1; i < h; ++i)
+            {
+                areaNormal += ( (polygon[2*i] - polygon[0]) ^ (polygon[2*i + 1] - polygon[2*i-1] ) );
+            }
+
+            // Last triangle or quad
+            areaNormal += ( (polygon[2*h] - polygon[0]) ^ (polygon[k] - polygon[2*h-1] ) );
+
+            areaNormal *= 0.5;
+
+            return areaNormal;
+        }
+    }
+}
+
 
 
 

ApplicationCode/ReservoirDataModel/cvfGeometryTools.h

@@ -27,6 +27,8 @@ public:
     static double     getAngle(const cvf::Vec3d& positiveNormalAxis, const cvf::Vec3d& v1, const cvf::Vec3d& v2);
     static double     getAngle(const cvf::Vec3d& v1, const cvf::Vec3d& v2);
 
+    static cvf::Vec3d polygonAreaNormal3D(const std::vector<cvf::Vec3d>& polygon);
+
     enum IntersectionStatus
     {
         NO_INTERSECTION,