Constified some methods in FemPart. WIP range filter mapping first untested code.

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPart.cpp

@@ -71,7 +71,7 @@ void RigFemPart::appendElement(RigElementType elmType, int id, const int* connec
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-const RigFemPartGrid* RigFemPart::structGrid()
+const RigFemPartGrid* RigFemPart::structGrid() const
 {
     if (m_structGrid.isNull())
     {
@@ -250,7 +250,7 @@ void RigFemPart::calculateElmNeighbors()
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-cvf::Vec3f RigFemPart::faceNormal(int elmIdx, int faceIdx)
+cvf::Vec3f RigFemPart::faceNormal(int elmIdx, int faceIdx) const
 {
     const std::vector<cvf::Vec3f>& nodeCoordinates = this->nodes().coordinates;
 

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPart.h

@@ -82,9 +82,9 @@ public:
     const std::vector<int>&     possibleGridCornerElements() const { return m_possibleGridCornerElements; }
     void                        findIntersectingCells(const cvf::BoundingBox& inputBB, std::vector<size_t>* elementIndices) const;
 
-    cvf::Vec3f                  faceNormal(int elmentIndex, int faceIndex);
+    cvf::Vec3f                  faceNormal(int elmentIndex, int faceIndex) const;
 
-    const RigFemPartGrid*       structGrid();   
+    const RigFemPartGrid*       structGrid() const;   
 
 private:
     int                         m_elementPartId;
@@ -96,7 +96,7 @@ private:
 
     RigFemPartNodes             m_nodes;
 
-    cvf::ref<RigFemPartGrid>    m_structGrid;
+    mutable cvf::ref<RigFemPartGrid>    m_structGrid;
 
     void calculateNodeToElmRefs();
     std::vector<std::vector<int> > m_nodeToElmRefs; // Needs a more memory friendly structure

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartGrid.cpp

@@ -26,7 +26,7 @@
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-RigFemPartGrid::RigFemPartGrid(RigFemPart* femPart)
+RigFemPartGrid::RigFemPartGrid(const RigFemPart* femPart)
 {
     m_femPart = femPart;
     generateStructGridData();
@@ -172,6 +172,16 @@ void RigFemPartGrid::generateStructGridData()
 
         if (kCoord > static_cast<int>(m_elmentIJKCounts[2])) m_elmentIJKCounts[2] = kCoord;
     }
+
+
+    
+    m_elmIdxPrIJK.resize(m_elmentIJKCounts[0], m_elmentIJKCounts[1],m_elmentIJKCounts[2]);
+    
+    for (int elmIdx = 0; elmIdx < m_femPart->elementCount(); ++elmIdx)
+    {
+        cvf::Vec3i ijk = m_ijkPrElement[elmIdx];
+        m_elmIdxPrIJK.at(ijk[0], ijk[1],  ijk[2]) = elmIdx;
+    }
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -380,8 +390,7 @@ bool RigFemPartGrid::cellIJKNeighbor(size_t i, size_t j, size_t k, FaceType face
 //--------------------------------------------------------------------------------------------------
 size_t RigFemPartGrid::cellIndexFromIJK(size_t i, size_t j, size_t k) const
 {
-    CVF_ASSERT(false);
+    return m_elmIdxPrIJK.at(i,j,k);
-    return cvf::UNDEFINED_SIZE_T;
 }
 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartGrid.h

@@ -27,10 +27,11 @@ class RigFemPart;
 class RigFemPartGrid : public cvf::StructGridInterface
 {
 public:
-    RigFemPartGrid(RigFemPart* femPart);
+    RigFemPartGrid(const RigFemPart* femPart);
     virtual ~RigFemPartGrid();
 
     virtual bool        ijkFromCellIndex(size_t cellIndex, size_t* i, size_t* j, size_t* k) const;
+    virtual size_t      cellIndexFromIJK(size_t i, size_t j, size_t k) const;
 
     virtual size_t      gridPointCountI() const;
     virtual size_t      gridPointCountJ() const;
@@ -43,7 +44,7 @@ public:
     int                 findElmIdxForIJK000();
     int                 perpendicularFaceInDirection(cvf::Vec3f direction, int perpFaceIdx, int elmIdx);
 
-    RigFemPart*         m_femPart;
+    const RigFemPart*  m_femPart;
 
     std::vector<cvf::Vec3i> m_ijkPrElement;
     cvf::Vec3st         m_elmentIJKCounts;
@@ -53,7 +54,6 @@ private: // Unused, Not implemented
     virtual cvf::Vec3d  minCoordinate() const;
     virtual cvf::Vec3d  maxCoordinate() const;
     virtual bool        cellIJKNeighbor(size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex) const;
-    virtual size_t      cellIndexFromIJK(size_t i, size_t j, size_t k) const;
 
 
     virtual bool        cellIJKFromCoordinate(const cvf::Vec3d& coord, size_t* i, size_t* j, size_t* k) const;
@@ -63,6 +63,37 @@ private: // Unused, Not implemented
     virtual size_t      gridPointIndexFromIJK(size_t i, size_t j, size_t k) const;
     virtual cvf::Vec3d  gridPointCoordinate(size_t i, size_t j, size_t k) const;
 
+    class IJKArray
+    {
+    public:
+        IJKArray(): m_iCount(0), m_jCount(0){}
+
+        void resize(size_t iCount, size_t jCount, size_t kCount)
+        {
+            data.resize(iCount*jCount*kCount, cvf::UNDEFINED_SIZE_T);
+            m_iCount = iCount;
+            m_jCount = jCount;
+        }
+
+        size_t& at(size_t i, size_t j, size_t k)
+        {
+            return data[i + j* m_iCount + k * m_iCount*m_jCount];
+        }
+        
+        size_t at(size_t i, size_t j, size_t k) const
+        {
+            return data[i + j* m_iCount + k * m_iCount*m_jCount];
+        }
+
+    private:
+        size_t m_iCount;
+        size_t m_jCount;
+
+        std::vector< size_t > data;
+    };
+
+    IJKArray m_elmIdxPrIJK; 
+
 };
 
 

ApplicationCode/ReservoirDataModel/RigCaseToCaseCellMapper.cpp

@@ -382,7 +382,7 @@ int RigCaseToCaseMapperTools::quadVxClosestToXYOfPoint( const cvf::Vec3d point,
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-bool RigCaseToCaseMapperTools::elementCorners(RigFemPart* femPart, int elmIdx, cvf::Vec3d elmCorners[8])
+bool RigCaseToCaseMapperTools::elementCorners(const RigFemPart* femPart, int elmIdx, cvf::Vec3d elmCorners[8])
 {
     RigElementType elmType = femPart->elementType(elmIdx);
     if (!(elmType == HEX8 || elmType == HEX8P)) return false;
@@ -556,8 +556,6 @@ void RigCaseToCaseCellMapper::calculateEclToGeomCellMapping(RigMainGrid* masterE
         std::vector<size_t> closeElements;
         dependentFemPart->findIntersectingCells(elmBBox, &closeElements);
 
-        std::vector<int> matchingCells;
-
         for (size_t ccIdx = 0; ccIdx < closeElements.size(); ++ccIdx)
         {
             int elmIdx = static_cast<int>(closeElements[ccIdx]);
@@ -651,11 +649,80 @@ for (size_t ij = 0; ij < minIJCount; ++ij )
 }
 #endif
 
-class RigCaseToCaseSingleCellMapper
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+cvf::Vec3d RigCaseToCaseMapperTools::calculateCellCenter(cvf::Vec3d elmCorners[8])
+{
+    cvf::Vec3d avg(cvf::Vec3d::ZERO);
+
+    size_t i;
+    for (i = 0; i < 8; i++)
+    {
+        avg += elmCorners[i];
+    }
+
+    avg /= 8.0;
+
+    return avg;
+}
+
+#include "RimCellRangeFilter.h"
+
+class RigCaseToCaseRangeFilterMapper
 {
 public:
-    static void findBestFemCellFromEclCell(RigMainGrid* masterEclGrid, int ei, int ej, int ek,
+   
-                                           RigFemPart* dependentFemPart, int* fi, int * fj, int* fk)
+
+    static void convertRangeFilter(RimCellRangeFilter* srcFilter, const RigMainGrid* srcEclGrid,
+                                   RimCellRangeFilter* dstFilter, const RigFemPart*  dstFemPart)
+    {
+        CVF_ASSERT(srcFilter && srcEclGrid && dstFilter && dstFemPart);
+        CVF_ASSERT(srcFilter->gridIndex() == 0); // LGR not supported yet
+        
+
+        size_t startFi, startFj, startFk;
+        bool isStartExactMatch  = findBestFemCellFromEclCell(srcEclGrid, 
+                                   srcFilter->startIndexI() -1, 
+                                   srcFilter->startIndexJ() -1, 
+                                   srcFilter->startIndexK() -1 ,
+                                   dstFemPart, &startFi, &startFj, &startFk);
+
+ 
+        size_t endFi, endFj, endFk;
+        bool isEndExactMatch  = findBestFemCellFromEclCell(srcEclGrid,
+                                   srcFilter->startIndexI() -1 + srcFilter->cellCountI(),
+                                   srcFilter->startIndexJ() -1 + srcFilter->cellCountJ(),
+                                   srcFilter->startIndexK() -1 + srcFilter->cellCountK(),
+                                   dstFemPart, &endFi, &endFj, &endFk);
+
+        if (   (startFi != cvf::UNDEFINED_SIZE_T && startFj != cvf::UNDEFINED_SIZE_T && startFk != cvf::UNDEFINED_SIZE_T)
+            && (endFi   != cvf::UNDEFINED_SIZE_T && endFj   != cvf::UNDEFINED_SIZE_T && endFk   != cvf::UNDEFINED_SIZE_T))
+        {
+            dstFilter->startIndexJ = static_cast<int>(startFj + 1);
+            dstFilter->startIndexK = static_cast<int>(startFk + 1);
+            dstFilter->startIndexI = static_cast<int>(startFi + 1);
+
+            dstFilter->cellCountI = static_cast<int>(endFi - (startFi-1));
+            dstFilter->cellCountJ = static_cast<int>(endFj - (startFj-1));
+            dstFilter->cellCountK = static_cast<int>(endFk - (startFk-1));
+            dstFilter->computeAndSetValidValues();
+        }
+        else
+        {
+            dstFilter->startIndexI = 1;
+            dstFilter->startIndexJ = 1;
+            dstFilter->startIndexK = 1;
+
+            dstFilter->cellCountI = 0;
+            dstFilter->cellCountJ = 0;
+            dstFilter->cellCountK = 0;
+            dstFilter->computeAndSetValidValues();
+        }
+    }
+
+    static bool findBestFemCellFromEclCell(const RigMainGrid* masterEclGrid, size_t ei, size_t ej, size_t ek,
+                                           const RigFemPart* dependentFemPart, size_t* fi, size_t * fj, size_t* fk)
     {
         // Find tolerance
 
@@ -667,8 +734,6 @@ public:
 
         bool isEclFaceNormalsOutwards = masterEclGrid->isFaceNormalsOutwards();
 
-        cvf::Vec3d elmCorners[8];
-
         size_t cellIdx =  masterEclGrid->cellIndexFromIJK( ei,  ej,  ek);
 
         cvf::Vec3d geoMechConvertedEclCell[8];
@@ -680,25 +745,124 @@ public:
         std::vector<size_t> closeElements;
         dependentFemPart->findIntersectingCells(elmBBox, &closeElements);
 
-        std::vector<int> matchingCells;
+        cvf::Vec3d elmCorners[8];
+        int elmIdxToBestMatch = -1;
+        double sqDistToClosestElmCenter = HUGE_VAL;
+        cvf::Vec3d convEclCellCenter = RigCaseToCaseMapperTools::calculateCellCenter(geoMechConvertedEclCell);
+
+        bool foundExactMatch = false;
 
         for (size_t ccIdx = 0; ccIdx < closeElements.size(); ++ccIdx)
         {
             int elmIdx = static_cast<int>(closeElements[ccIdx]);
 
             RigCaseToCaseMapperTools::elementCorners(dependentFemPart, elmIdx, elmCorners);
+            
+            cvf::Vec3d cellCenter = RigCaseToCaseMapperTools::calculateCellCenter(elmCorners);
+            double sqDist = (cellCenter - convEclCellCenter).lengthSquared();
+            if (sqDist < sqDistToClosestElmCenter) 
+            {
+                elmIdxToBestMatch = elmIdx;
+                sqDistToClosestElmCenter = sqDist;
+            }
 
             RigCaseToCaseMapperTools::rotateCellTopologicallyToMatchBaseCell(geoMechConvertedEclCell, isEclFaceNormalsOutwards, elmCorners);
 
-            bool isMatching = RigCaseToCaseMapperTools::isEclFemCellsMatching(geoMechConvertedEclCell, elmCorners,
+            foundExactMatch = RigCaseToCaseMapperTools::isEclFemCellsMatching(geoMechConvertedEclCell, elmCorners,
                                                                               xyTolerance, zTolerance);
 
-            if (isMatching)
+            if (foundExactMatch)
             {
-               
+                elmIdxToBestMatch = elmIdx;
+                break;
             }
         }
 
+        if (elmIdxToBestMatch != -1)
+        {
+            dependentFemPart->structGrid()->ijkFromCellIndex(elmIdxToBestMatch, fi, fj, fk);
+        }
+        else
+        {
+            (*fi) = cvf::UNDEFINED_SIZE_T;
+            (*fj) = cvf::UNDEFINED_SIZE_T;
+            (*fk) = cvf::UNDEFINED_SIZE_T;
+        }
+
+        return foundExactMatch;
+
+    }
+
+
+    static bool findBestEclCellFromFemCell(const RigFemPart* dependentFemPart, size_t fi, size_t fj, size_t fk,
+                                           const RigMainGrid* masterEclGrid, size_t* ei, size_t* ej, size_t* ek)
+    {
+        // Find tolerance
+
+        double cellSizeI, cellSizeJ, cellSizeK;
+        masterEclGrid->characteristicCellSizes(&cellSizeI, &cellSizeJ, &cellSizeK);
+
+        double xyTolerance  = cellSizeI* 0.4;
+        double zTolerance   = cellSizeK* 0.4;
+
+        bool isEclFaceNormalsOutwards = masterEclGrid->isFaceNormalsOutwards();
+
+        int elementIdx =  static_cast<int>(dependentFemPart->structGrid()->cellIndexFromIJK( fi,  fj,  fk));
+
+        cvf::Vec3d elmCorners[8];
+        RigCaseToCaseMapperTools::elementCorners(dependentFemPart, elementIdx, elmCorners);
+
+        cvf::BoundingBox elmBBox;
+        for (int i = 0; i < 8 ; ++i) elmBBox.add(elmCorners[i]);
+
+        std::vector<size_t> closeCells;
+        masterEclGrid->findIntersectingCells(elmBBox, &closeCells); // This might actually miss the exact one, but we have no other alternative yet.
+       
+        size_t globCellIdxToBestMatch = cvf::UNDEFINED_SIZE_T;
+        double sqDistToClosestCellCenter = HUGE_VAL;
+        cvf::Vec3d elmCenter = RigCaseToCaseMapperTools::calculateCellCenter(elmCorners);
+
+        bool foundExactMatch = false;
+
+        for (size_t ccIdx = 0; ccIdx < closeCells.size(); ++ccIdx)
+        {
+            size_t cellIdx = closeCells[ccIdx];
+            cvf::Vec3d geoMechConvertedEclCell[8];
+            RigCaseToCaseMapperTools::estimatedFemCellFromEclCell(masterEclGrid, cellIdx, geoMechConvertedEclCell);
+            
+            cvf::Vec3d cellCenter = RigCaseToCaseMapperTools::calculateCellCenter(geoMechConvertedEclCell);
+            double sqDist = (cellCenter - elmCenter).lengthSquared();
+            if (sqDist < sqDistToClosestCellCenter) 
+            {
+                globCellIdxToBestMatch = cellIdx;
+                sqDistToClosestCellCenter = sqDist;
+            }
+
+            RigCaseToCaseMapperTools::rotateCellTopologicallyToMatchBaseCell(geoMechConvertedEclCell, isEclFaceNormalsOutwards, elmCorners);
+
+            foundExactMatch = RigCaseToCaseMapperTools::isEclFemCellsMatching(geoMechConvertedEclCell, elmCorners,
+                                                                              xyTolerance, zTolerance);
+
+            if (foundExactMatch)
+            {
+                globCellIdxToBestMatch = cellIdx;
+                break;
+            }
+        }
+
+        if (globCellIdxToBestMatch != cvf::UNDEFINED_SIZE_T)
+        {
+            masterEclGrid->ijkFromCellIndex(globCellIdxToBestMatch, ei, ej, ek);
+        }
+        else
+        {
+            (*ei) = cvf::UNDEFINED_SIZE_T;
+            (*ej) = cvf::UNDEFINED_SIZE_T;
+            (*ek) = cvf::UNDEFINED_SIZE_T;
+        }
+
+        return foundExactMatch;
+       
     }
 
 };

ApplicationCode/ReservoirDataModel/RigCaseToCaseCellMapper.h

@@ -53,6 +53,7 @@ private:
     void storeMapping(int depCaseCellIdx, const std::vector<int>& masterCaseMatchingCells);
     void addMapping(int depCaseCellIdx, int masterCaseMatchingCell);
     void calculateEclToGeomCellMapping(RigMainGrid* masterEclGrid, RigFemPart* dependentFemPart, bool eclipseIsMaster);
+
     std::vector<int>    m_masterCellOrIntervalIndex;
     std::vector<std::vector<int> > m_masterCellIndexSeries;
 
@@ -75,8 +76,9 @@ public:
     static void rotateQuad(cvf::Vec3d quad[4], int idxToNewStart);
     static void flipQuadWinding(cvf::Vec3d quad[4]);
     static int  quadVxClosestToXYOfPoint(const cvf::Vec3d point, const cvf::Vec3d quad[4]);
-    static bool elementCorners(RigFemPart* femPart, int elmIdx, cvf::Vec3d elmCorners[8]);
+    static bool elementCorners(const RigFemPart* femPart, int elmIdx, cvf::Vec3d elmCorners[8]);
     static int  findMatchingPOSKFaceIdx(const cvf::Vec3d baseCell[8], bool isBaseCellNormalsOutwards, const cvf::Vec3d c2[8]);
     static bool isEclFemCellsMatching(const cvf::Vec3d baseCell[8], cvf::Vec3d cell[8], double xyTolerance, double zTolerance);
     static void rotateCellTopologicallyToMatchBaseCell(const cvf::Vec3d * baseCell, bool baseCellFaceNormalsIsOutwards, cvf::Vec3d * cell);
+    static cvf::Vec3d calculateCellCenter(cvf::Vec3d elmCorners[8]);
 };