diff --git a/ApplicationLibCode/ProjectDataModel/Completions/RimEnsembleFractureStatistics.cpp b/ApplicationLibCode/ProjectDataModel/Completions/RimEnsembleFractureStatistics.cpp
index 82c98c97bb..f97303a47d 100644
--- a/ApplicationLibCode/ProjectDataModel/Completions/RimEnsembleFractureStatistics.cpp
+++ b/ApplicationLibCode/ProjectDataModel/Completions/RimEnsembleFractureStatistics.cpp
@@ -899,17 +899,6 @@ double RimEnsembleFractureStatistics::linearSampling( double               minVa
     return sampleDistance;
 }
 
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-bool RimEnsembleFractureStatistics::isCoordinateInsideFractureCell( double x, double y, const RigFractureCell& cell )
-{
-    const cvf::Vec3d& minPoint = cell.getPolygon()[0];
-    const cvf::Vec3d& maxPoint = cell.getPolygon()[2];
-    // TODO: Investigate strange ordering for y coords.
-    return ( x > minPoint.x() && x <= maxPoint.x() && y <= minPoint.y() && y > maxPoint.y() );
-}
-
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -939,24 +928,24 @@ void RimEnsembleFractureStatistics::sampleAllGrids( const std::vector<cvf::cref<
                                                     const std::vector<double>&                     samplesY,
                                                     std::vector<std::vector<double>>&              samples )
 {
-    for ( size_t y = 0; y < samplesY.size(); y++ )
+    const int ny = samplesY.size();
+#pragma omp parallel for
+    for ( int y = 0; y < ny; y++ )
     {
         for ( size_t x = 0; x < samplesX.size(); x++ )
         {
-            double posX = samplesX[x];
-            double posY = samplesY[y];
+            double     posX = samplesX[x];
+            double     posY = samplesY[y];
+            cvf::Vec3d pos( posX, posY, 0.0 );
 
             for ( auto fractureGrid : fractureGrids )
             {
-                for ( auto fractureCell : fractureGrid->fractureCells() )
+                const RigFractureCell* fractureCell = fractureGrid->getCellFromPosition( pos );
+                if ( fractureCell )
                 {
-                    if ( isCoordinateInsideFractureCell( posX, posY, fractureCell ) )
-                    {
-                        size_t idx   = y * samplesX.size() + x;
-                        double value = fractureCell.getConductivityValue();
-                        if ( !std::isinf( value ) ) samples[idx].push_back( value );
-                        break;
-                    }
+                    size_t idx   = y * samplesX.size() + x;
+                    double value = fractureCell->getConductivityValue();
+                    if ( !std::isinf( value ) ) samples[idx].push_back( value );
                 }
             }
         }
@@ -1032,7 +1021,9 @@ void RimEnsembleFractureStatistics::generateStatisticsGrids(
 
     RigSlice2D occurrenceGrid( numSamplesX, numSamplesY );
 
-    for ( size_t y = 0; y < numSamplesY; y++ )
+    const int ny = static_cast<int>( numSamplesY );
+#pragma omp parallel for
+    for ( int y = 0; y < ny; y++ )
     {
         for ( size_t x = 0; x < numSamplesX; x++ )
         {
diff --git a/ApplicationLibCode/ReservoirDataModel/RigFractureGrid.cpp b/ApplicationLibCode/ReservoirDataModel/RigFractureGrid.cpp
index f3fa79d41d..cf6b1362a0 100644
--- a/ApplicationLibCode/ReservoirDataModel/RigFractureGrid.cpp
+++ b/ApplicationLibCode/ReservoirDataModel/RigFractureGrid.cpp
@@ -20,6 +20,9 @@
 
 #include "RiaLogging.h"
 
+#include "cvfBoundingBox.h"
+#include "cvfBoundingBoxTree.h"
+
 #include <QString>
 
 //--------------------------------------------------------------------------------------------------
@@ -122,3 +125,54 @@ std::pair<size_t, size_t> RigFractureGrid::fractureCellAtWellCenter() const
 {
     return m_wellCenterFractureCellIJ;
 }
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RigFractureGrid::ensureCellSearchTreeIsBuilt()
+{
+    if ( m_cellBoundingBoxTree.isNull() )
+    {
+        size_t itemCount = m_fractureCells.size();
+
+        std::vector<cvf::BoundingBox> cellBoundingBoxes( itemCount );
+        std::vector<size_t>           boundingBoxIds( itemCount );
+
+        for ( size_t idx = 0; idx < itemCount; ++idx )
+        {
+            const RigFractureCell&         cell    = m_fractureCells[idx];
+            cvf::BoundingBox&              cellBB  = cellBoundingBoxes[idx];
+            const std::vector<cvf::Vec3d>& corners = cell.getPolygon();
+            for ( auto c : corners )
+                cellBB.add( c );
+
+            boundingBoxIds[idx] = idx;
+        }
+
+        m_cellBoundingBoxTree = new cvf::BoundingBoxTree;
+        m_cellBoundingBoxTree->buildTreeFromBoundingBoxes( cellBoundingBoxes, &boundingBoxIds );
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+const RigFractureCell* RigFractureGrid::getCellFromPosition( const cvf::Vec3d& position ) const
+{
+    cvf::BoundingBox inputBB;
+    inputBB.add( position );
+
+    std::vector<size_t> indexes;
+    m_cellBoundingBoxTree->findIntersections( inputBB, &indexes );
+
+    if ( !indexes.empty() )
+    {
+        // Hit: should only one cell since they have no overlap
+        return &m_fractureCells[indexes[0]];
+    }
+    else
+    {
+        // No hit
+        return nullptr;
+    }
+}
diff --git a/ApplicationLibCode/ReservoirDataModel/RigFractureGrid.h b/ApplicationLibCode/ReservoirDataModel/RigFractureGrid.h
index 526ec80957..655ac88f1b 100644
--- a/ApplicationLibCode/ReservoirDataModel/RigFractureGrid.h
+++ b/ApplicationLibCode/ReservoirDataModel/RigFractureGrid.h
@@ -23,6 +23,12 @@
 
 #include <vector>
 
+namespace cvf
+{
+class BoundingBox;
+class BoundingBoxTree;
+} // namespace cvf
+
 class RigFractureCell;
 
 //==================================================================================================
@@ -47,9 +53,14 @@ public:
 
     std::pair<size_t, size_t> fractureCellAtWellCenter() const;
 
+    const RigFractureCell* getCellFromPosition( const cvf::Vec3d& position ) const;
+    void                   ensureCellSearchTreeIsBuilt();
+
 private:
     std::vector<RigFractureCell> m_fractureCells;
     std::pair<size_t, size_t>    m_wellCenterFractureCellIJ;
     size_t                       m_iCellCount;
     size_t                       m_jCellCount;
+
+    cvf::ref<cvf::BoundingBoxTree> m_cellBoundingBoxTree;
 };
diff --git a/ApplicationLibCode/ReservoirDataModel/RigStimPlanFractureDefinition.cpp b/ApplicationLibCode/ReservoirDataModel/RigStimPlanFractureDefinition.cpp
index ccef694bd5..d441c737e7 100644
--- a/ApplicationLibCode/ReservoirDataModel/RigStimPlanFractureDefinition.cpp
+++ b/ApplicationLibCode/ReservoirDataModel/RigStimPlanFractureDefinition.cpp
@@ -487,6 +487,7 @@ cvf::cref<RigFractureGrid>
     fractureGrid->setWellCenterFractureCellIJ( wellCenterStimPlanCellIJ );
     fractureGrid->setICellCount( this->m_Xs.size() - 2 );
     fractureGrid->setJCellCount( this->m_Ys.size() - 2 );
+    fractureGrid->ensureCellSearchTreeIsBuilt();
 
     return cvf::cref<RigFractureGrid>( fractureGrid.p() );
 }