mirror of
https://github.com/OPM/ResInsight.git
synced 2026-10-24 14:17:52 -05:00
#14128 Use KD-tree (nanoflann) for fault distance nearest-face search
computeFaultDistances faked a nearest-point query with cvf::BoundingBoxTree, which only supports box overlap: it grew a cell-sized box until it overlapped a fault face. For cells far from any fault this did many full tree traversals and then brute-forced the distance to a large, over-collected candidate set, and the answer was only approximate. This was slow on grids with ~10 million cells. The metric is the distance from each cell center to the nearest fault face center, which is exactly a nearest-point search. Replace the bounding-box tree and the grow-the-box loop with a nanoflann KD-tree built over the fault face centers and a single nearest-neighbor query per active cell. A zero-copy adaptor reads the existing std::vector<cvf::Vec3d> directly. The OpenMP loop is retained since nanoflann queries are read-only. The result is O(F log F) build plus O(N log F) queries, the full grid-node array copy is gone, and distances are now exact rather than approximate. Add nanoflann as a vcpkg dependency, wired in like Clipper2 as an imported target. Refs: https://github.com/OPM/ResInsight/issues/14128
This commit is contained in:
1 parent
5dae9a0a19
commit
5ee62c61e1
4 files changed
+47
-59
No files matched your search
@@ -307,6 +307,7 @@ endif()
|
||||
set(LINK_LIBRARIES
|
||||
Eigen3::Eigen
|
||||
Clipper2::Clipper2Z
|
||||
nanoflann::nanoflann
|
||||
${THIRD_PARTY_LIBRARIES}
|
||||
${OPENGL_LIBRARIES}
|
||||
${QT_LIBRARIES}
|
||||
|
||||
+39
-58
@@ -23,10 +23,34 @@
|
||||
#include "RigFault.h"
|
||||
#include "RigMainGrid.h"
|
||||
|
||||
#include "cvfBoundingBoxTree.h"
|
||||
#include <nanoflann.hpp>
|
||||
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
|
||||
namespace
|
||||
{
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Adaptor exposing the collected fault-face-center points to nanoflann without copying them.
|
||||
/// cvf::Vec3d stores its three doubles contiguously, so ptr()[dim] is a direct coordinate lookup.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
struct FaceCenterCloud
|
||||
{
|
||||
const std::vector<cvf::Vec3d>& points;
|
||||
|
||||
inline size_t kdtree_get_point_count() const { return points.size(); }
|
||||
inline double kdtree_get_pt( size_t idx, size_t dim ) const { return points[idx].ptr()[dim]; }
|
||||
|
||||
template <class BBOX>
|
||||
bool kdtree_get_bbox( BBOX& ) const
|
||||
{
|
||||
return false;
|
||||
}
|
||||
};
|
||||
|
||||
using FaceCenterKdTree = nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, FaceCenterCloud>, FaceCenterCloud, 3>;
|
||||
} // namespace
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
@@ -61,25 +85,12 @@ void RigFaultDistanceCalculator::computeFaultDistances( const RigMainGrid*
|
||||
|
||||
if ( faultFaceCenters.empty() ) return;
|
||||
|
||||
// Create bounding box tree for all face centers
|
||||
cvf::BoundingBoxTree searchTree;
|
||||
{
|
||||
std::vector<size_t> faceIndicesForBoundingBoxes;
|
||||
std::vector<cvf::BoundingBox> faceBBs;
|
||||
|
||||
size_t faceCenterIndex = 0;
|
||||
for ( const auto& faultFaceCenter : faultFaceCenters )
|
||||
{
|
||||
cvf::BoundingBox bb;
|
||||
bb.add( faultFaceCenter );
|
||||
faceBBs.push_back( bb );
|
||||
faceIndicesForBoundingBoxes.push_back( faceCenterIndex++ );
|
||||
}
|
||||
searchTree.buildTreeFromBoundingBoxes( faceBBs, &faceIndicesForBoundingBoxes );
|
||||
}
|
||||
|
||||
const auto nodes = mainGrid->nodes();
|
||||
const auto mainGridBB = mainGrid->boundingBox();
|
||||
// Build a KD-tree over the fault face centers and query the nearest one for every active cell.
|
||||
// This is an exact nearest-point search in O(log N) per cell, replacing the previous expanding
|
||||
// bounding-box heuristic.
|
||||
FaceCenterCloud cloud{ faultFaceCenters };
|
||||
FaceCenterKdTree kdTree( 3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams( 10 /* max leaf size */ ) );
|
||||
kdTree.buildIndex();
|
||||
|
||||
#pragma omp parallel for
|
||||
for ( int activeIndex = 0; activeIndex < static_cast<int>( activeCells.size() ); activeIndex++ )
|
||||
@@ -90,45 +101,15 @@ void RigFaultDistanceCalculator::computeFaultDistances( const RigMainGrid*
|
||||
const RigCell& cell = mainGrid->cell( cellIdx.value() );
|
||||
if ( cell.isInvalid() ) continue;
|
||||
|
||||
std::vector<size_t> candidateFaceIndices;
|
||||
{
|
||||
cvf::BoundingBox bb;
|
||||
const auto& cellIndices = cell.cornerIndices();
|
||||
for ( const auto& i : cellIndices )
|
||||
{
|
||||
bb.add( nodes[i] );
|
||||
}
|
||||
const cvf::Vec3d cellCenter = cell.center();
|
||||
const double queryPoint[] = { cellCenter.x(), cellCenter.y(), cellCenter.z() };
|
||||
|
||||
searchTree.findIntersections( bb, &candidateFaceIndices );
|
||||
size_t nearestIndex = 0;
|
||||
double nearestDistanceSquared = std::numeric_limits<double>::infinity();
|
||||
nanoflann::KNNResultSet<double> resultSet( 1 );
|
||||
resultSet.init( &nearestIndex, &nearestDistanceSquared );
|
||||
kdTree.findNeighbors( resultSet, queryPoint );
|
||||
|
||||
bool bbIsBelowThreshold = true;
|
||||
while ( candidateFaceIndices.empty() && bbIsBelowThreshold )
|
||||
{
|
||||
if ( bb.extent().x() > mainGridBB.extent().x() * 2 )
|
||||
{
|
||||
bbIsBelowThreshold = false;
|
||||
break;
|
||||
}
|
||||
if ( bb.extent().y() > mainGridBB.extent().y() * 2 )
|
||||
{
|
||||
bbIsBelowThreshold = false;
|
||||
break;
|
||||
}
|
||||
|
||||
bb.expand( bb.extent().x() );
|
||||
searchTree.findIntersections( bb, &candidateFaceIndices );
|
||||
}
|
||||
}
|
||||
|
||||
// Find closest fault face
|
||||
double shortestDistance = std::numeric_limits<double>::infinity();
|
||||
|
||||
for ( const auto& faultFaceIndex : candidateFaceIndices )
|
||||
{
|
||||
const cvf::Vec3d& faultFaceCenter = faultFaceCenters[faultFaceIndex];
|
||||
shortestDistance = std::min( cell.center().pointDistance( faultFaceCenter ), shortestDistance );
|
||||
}
|
||||
|
||||
resultValues[activeIndex] = shortestDistance;
|
||||
resultValues[activeIndex] = std::sqrt( nearestDistanceSquared );
|
||||
}
|
||||
}
|
||||
@@ -758,6 +758,11 @@ find_package(Clipper2 CONFIG REQUIRED)
|
||||
# THIRD_PARTY_LIBRARIES which is used for setting folder properties on targets
|
||||
# we build
|
||||
|
||||
find_package(nanoflann CONFIG REQUIRED)
|
||||
|
||||
# nanoflann::nanoflann is a header-only imported target, so like Clipper2 it is
|
||||
# not added to THIRD_PARTY_LIBRARIES
|
||||
|
||||
# ##############################################################################
|
||||
# roffcpp
|
||||
# ##############################################################################
|
||||
|
||||
+2
-1
@@ -13,7 +13,8 @@
|
||||
"type-lite",
|
||||
"fast-float",
|
||||
"spdlog",
|
||||
"pugixml"
|
||||
"pugixml",
|
||||
"nanoflann"
|
||||
],
|
||||
"overrides": [
|
||||
{
|
||||
|
||||
Reference in new issue
Block a user