diff --git a/ApplicationLibCode/Application/Tools/RiaEnsembleNameTools.cpp b/ApplicationLibCode/Application/Tools/RiaEnsembleNameTools.cpp index b8b8fa1390..d8af83d049 100644 --- a/ApplicationLibCode/Application/Tools/RiaEnsembleNameTools.cpp +++ b/ApplicationLibCode/Application/Tools/RiaEnsembleNameTools.cpp @@ -120,12 +120,23 @@ QString RiaEnsembleNameTools::uniqueShortName( const QString& sourceFileName const QStringList& allFileNames, const QString& ensembleCaseName ) { - QRegularExpression trimRe( "^[^a-zA-Z0-9]+" ); - std::map keyFileComponentsForAllFiles = RiaFilePathTools::keyPathComponentsForEachFilePath( allFileNames ); - QStringList keyFileComponents = keyFileComponentsForAllFiles[sourceFileName]; + return uniqueShortNameFromComponents( sourceFileName, keyFileComponentsForAllFiles, ensembleCaseName ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +QString RiaEnsembleNameTools::uniqueShortNameFromComponents( const QString& sourceFileName, + const std::map& keyFileComponentsForAllFiles, + const QString& ensembleCaseName ) +{ + QRegularExpression trimRe( "^[^a-zA-Z0-9]+" ); + + auto modifyableMap( keyFileComponentsForAllFiles ); + QStringList keyFileComponents = modifyableMap[sourceFileName]; if ( keyFileComponents.empty() ) return "Unnamed"; if ( !ensembleCaseName.isEmpty() ) diff --git a/ApplicationLibCode/Application/Tools/RiaEnsembleNameTools.h b/ApplicationLibCode/Application/Tools/RiaEnsembleNameTools.h index 93b288a517..58474656c6 100644 --- a/ApplicationLibCode/Application/Tools/RiaEnsembleNameTools.h +++ b/ApplicationLibCode/Application/Tools/RiaEnsembleNameTools.h @@ -19,6 +19,7 @@ #pragma once #include +#include //================================================================================================== // @@ -32,4 +33,8 @@ public: static QString uniqueShortName( const QString& sourceFileName, const QStringList& allFileNames, const QString& ensembleCaseName = QString() ); + + static QString uniqueShortNameFromComponents( const QString& sourceFileName, + const std::map& keyFileComponentsForAllFiles, + const QString& ensembleCaseName ); }; diff --git a/ApplicationLibCode/Commands/RicCreateEnsembleSurfaceFeature.cpp b/ApplicationLibCode/Commands/RicCreateEnsembleSurfaceFeature.cpp index ff4d97e9c3..d7831c3061 100644 --- a/ApplicationLibCode/Commands/RicCreateEnsembleSurfaceFeature.cpp +++ b/ApplicationLibCode/Commands/RicCreateEnsembleSurfaceFeature.cpp @@ -84,7 +84,7 @@ void RicCreateEnsembleSurfaceFeature::openDialogAndExecuteCommand() if ( propertyDialog.exec() == QDialog::Accepted ) { - executeCommand( *ui, result.files.toStdList() ); + executeCommand( *ui, result.files.toVector().toStdVector() ); } } @@ -92,19 +92,32 @@ void RicCreateEnsembleSurfaceFeature::openDialogAndExecuteCommand() /// //-------------------------------------------------------------------------------------------------- void RicCreateEnsembleSurfaceFeature::executeCommand( const RicCreateEnsembleSurfaceUi& ui, - const std::list& fileNames ) + const std::vector& fileNames ) { std::vector layers = ui.layers(); caf::ProgressInfo progress( fileNames.size(), "Generating ensemble surface" ); QStringList allSurfaceFileNames; - for ( auto fileName : fileNames ) + + auto fileCount = static_cast( fileNames.size() ); +#pragma omp parallel for + for ( int i = 0; i < fileCount; i++ ) { - auto task = progress.task( QString( "Extracting surfaces for %1" ).arg( fileName ) ); - auto [isOk, surfaceFileNames] = RimcCommandRouter_extractSurfaces::extractSurfaces( fileName, layers ); - if ( isOk ) allSurfaceFileNames << surfaceFileNames; + auto fileName = fileNames[i]; + + // Not possible to use structured bindings here due to a bug in clang + auto surfaceResult = RimcCommandRouter_extractSurfaces::extractSurfaces( fileName, layers ); + auto isOk = surfaceResult.first; + auto surfaceFileNames = surfaceResult.second; + +#pragma omp critical( RicCreateEnsembleSurfaceFeature ) + { + auto task = progress.task( QString( "Extracting surfaces for %1" ).arg( fileName ) ); + if ( isOk ) allSurfaceFileNames << surfaceFileNames; + } } + progress.setProgress( fileNames.size() ); if ( ui.autoCreateEnsembleSurfaces() ) RicImportEnsembleSurfaceFeature::importEnsembleSurfaceFromFiles( allSurfaceFileNames ); diff --git a/ApplicationLibCode/Commands/RicCreateEnsembleSurfaceFeature.h b/ApplicationLibCode/Commands/RicCreateEnsembleSurfaceFeature.h index 4ef62c1889..d95ee48397 100644 --- a/ApplicationLibCode/Commands/RicCreateEnsembleSurfaceFeature.h +++ b/ApplicationLibCode/Commands/RicCreateEnsembleSurfaceFeature.h @@ -33,7 +33,7 @@ class RicCreateEnsembleSurfaceFeature : public caf::CmdFeature public: static void openDialogAndExecuteCommand(); - static void executeCommand( const RicCreateEnsembleSurfaceUi& ui, const std::list& fileNames ); + static void executeCommand( const RicCreateEnsembleSurfaceUi& ui, const std::vector& fileNames ); protected: bool isCommandEnabled() override; diff --git a/ApplicationLibCode/Commands/RicImportEnsembleSurfaceFeature.cpp b/ApplicationLibCode/Commands/RicImportEnsembleSurfaceFeature.cpp index ceb910972b..1eb179438c 100644 --- a/ApplicationLibCode/Commands/RicImportEnsembleSurfaceFeature.cpp +++ b/ApplicationLibCode/Commands/RicImportEnsembleSurfaceFeature.cpp @@ -20,6 +20,7 @@ #include "RiaApplication.h" #include "RiaEnsembleNameTools.h" +#include "RiaFilePathTools.h" #include "RiaLogging.h" #include "RiaSummaryTools.h" @@ -75,42 +76,86 @@ void RicImportEnsembleSurfaceFeature::importEnsembleSurfaceFromFiles( const QStr { if ( fileNames.isEmpty() ) return; - QString ensembleName = RiaEnsembleNameTools::findSuitableEnsembleName( fileNames ); - QString layerName = RiaEnsembleNameTools::findCommonBaseName( fileNames ); - if ( !layerName.isEmpty() ) + // Create a list of file names for each layer + std::map fileNamesForEachLayer; + for ( const auto& name : fileNames ) { - ensembleName += QString( " : %1" ).arg( layerName ); + QFileInfo fi( name ); + + auto layerName = fi.baseName(); + fileNamesForEachLayer[layerName].push_back( name ); } - if ( ensembleName.isEmpty() ) ensembleName = "Ensemble Surface"; - - std::vector surfaces; - for ( const auto& fileName : fileNames ) + RimEnsembleSurface* ensembleToSelect = nullptr; + for ( const auto& fileNamesForLayer : fileNamesForEachLayer ) { - RimFileSurface* fileSurface = new RimFileSurface; - fileSurface->setSurfaceFilePath( fileName ); + auto filenames = fileNamesForLayer.second; + QString ensembleName = RiaEnsembleNameTools::findSuitableEnsembleName( fileNames ); - auto shortName = RiaEnsembleNameTools::uniqueShortName( fileName, fileNames ); - fileSurface->setUserDescription( shortName ); - - if ( fileSurface->onLoadData() ) + QString layerName = fileNamesForLayer.first; + if ( !layerName.isEmpty() ) { - surfaces.push_back( fileSurface ); + ensembleName += QString( " : %1" ).arg( layerName ); } + + if ( ensembleName.isEmpty() ) ensembleName = "Ensemble Surface"; + + std::map keyFileComponentsForAllFiles = + RiaFilePathTools::keyPathComponentsForEachFilePath( fileNames ); + + std::vector surfaces; + for ( size_t i = 0; i < fileNames.size(); i++ ) + { + surfaces.push_back( new RimFileSurface ); + } + + auto fileCount = static_cast( fileNames.size() ); +#pragma omp parallel for + for ( int i = 0; i < fileCount; i++ ) + { + auto fileName = fileNames[i]; + + auto fileSurface = surfaces[i]; + fileSurface->setSurfaceFilePath( fileName ); + + auto shortName = + RiaEnsembleNameTools::uniqueShortNameFromComponents( fileName, keyFileComponentsForAllFiles, ensembleName ); + fileSurface->setUserDescription( shortName ); + + auto isOk = fileSurface->onLoadData(); + if ( !isOk ) + { + delete fileSurface; + surfaces[i] = nullptr; + } + } + + { + // Remove null pointers from vector of surfaces + std::vector tmp; + for ( auto s : surfaces ) + { + if ( s != nullptr ) tmp.push_back( s ); + } + + surfaces.swap( tmp ); + } + + if ( surfaces.empty() ) return; + + auto ensemble = new RimEnsembleSurface; + ensemble->setCollectionName( ensembleName ); + for ( auto surface : surfaces ) + ensemble->addFileSurface( surface ); + + ensemble->loadDataAndUpdate(); + RimProject::current()->activeOilField()->surfaceCollection->addEnsembleSurface( ensemble ); + + ensembleToSelect = ensemble; } - if ( surfaces.empty() ) return; - - RimEnsembleSurface* ensemble = new RimEnsembleSurface; - ensemble->setCollectionName( ensembleName ); - for ( auto surface : surfaces ) - ensemble->addFileSurface( surface ); - - ensemble->loadDataAndUpdate(); - RimProject::current()->activeOilField()->surfaceCollection->addEnsembleSurface( ensemble ); - RimProject::current()->activeOilField()->surfaceCollection->updateConnectedEditors(); - Riu3DMainWindowTools::selectAsCurrentItem( ensemble ); + Riu3DMainWindowTools::selectAsCurrentItem( ensembleToSelect ); } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/FileInterface/RifSurfaceImporter.cpp b/ApplicationLibCode/FileInterface/RifSurfaceImporter.cpp index 975ab44152..1d4533ad6a 100644 --- a/ApplicationLibCode/FileInterface/RifSurfaceImporter.cpp +++ b/ApplicationLibCode/FileInterface/RifSurfaceImporter.cpp @@ -17,6 +17,7 @@ ///////////////////////////////////////////////////////////////////////////////// #include "RifSurfaceImporter.h" +#include "RiaStdStringTools.h" #include "RigGocadData.h" #include "cafProgressInfo.h" @@ -61,7 +62,12 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g std::vector> propertyValues; { - std::ifstream stream( filename.toLatin1().data() ); + std::stringstream stream; + { + // Read the file content into a stringstream to avoid expensive file operations + std::ifstream t( filename.toLatin1().data() ); + stream << t.rdbuf(); + } bool isInTfaceSection = false; GocadZPositive zDir = GocadZPositive::Unknown; @@ -71,75 +77,77 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g std::string line; std::getline( stream, line ); - std::transform( line.begin(), line.end(), line.begin(), ::toupper ); - - std::istringstream lineStream( line ); + auto tokens = RiaStdStringTools::splitString( line, ' ' ); std::string firstToken; - lineStream >> firstToken; + if ( !tokens.empty() ) firstToken = tokens.front(); if ( isInTfaceSection ) { if ( firstToken.compare( "VRTX" ) == 0 ) { - int vertexId = -1; - double x{ std::numeric_limits::infinity() }; - double y{ std::numeric_limits::infinity() }; - double z{ std::numeric_limits::infinity() }; - std::string endVertex; - - lineStream >> vertexId >> x >> y >> z >> endVertex; - - if ( vertexId > -1 ) + if ( tokens.size() > 4 ) { - if ( zDir == GocadZPositive::Depth ) - { - z = -z; - } + int vertexId = RiaStdStringTools::toInt( tokens[1] ); + double x = RiaStdStringTools::toDouble( tokens[2] ); + double y = RiaStdStringTools::toDouble( tokens[3] ); + double z = RiaStdStringTools::toDouble( tokens[4] ); - vertices.emplace_back( cvf::Vec3d( x, y, z ) ); - vertexIdToIndex[vertexId] = static_cast( vertices.size() - 1 ); + if ( vertexId > -1 ) + { + if ( zDir == GocadZPositive::Depth ) + { + z = -z; + } + + vertices.emplace_back( cvf::Vec3d( x, y, z ) ); + vertexIdToIndex[vertexId] = static_cast( vertices.size() - 1 ); + } } } else if ( firstToken.compare( "PVRTX" ) == 0 ) { - int vertexId = -1; - double x{ std::numeric_limits::infinity() }; - double y{ std::numeric_limits::infinity() }; - double z{ std::numeric_limits::infinity() }; - - lineStream >> vertexId >> x >> y >> z; - - if ( vertexId > -1 ) + if ( tokens.size() > 4 ) { - if ( zDir == GocadZPositive::Depth ) z = -z; + int vertexId = RiaStdStringTools::toInt( tokens[1] ); + double x = RiaStdStringTools::toDouble( tokens[2] ); + double y = RiaStdStringTools::toDouble( tokens[3] ); + double z = RiaStdStringTools::toDouble( tokens[4] ); - vertices.emplace_back( cvf::Vec3d( x, y, z ) ); - vertexIdToIndex[vertexId] = static_cast( vertices.size() - 1 ); - } + if ( vertexId > -1 ) + { + if ( zDir == GocadZPositive::Depth ) z = -z; - for ( size_t i = 0; i < propertyNames.size(); i++ ) - { - float value = std::numeric_limits::infinity(); + vertices.emplace_back( cvf::Vec3d( x, y, z ) ); + vertexIdToIndex[vertexId] = static_cast( vertices.size() - 1 ); - lineStream >> value; + for ( size_t i = 0; i < propertyNames.size(); i++ ) + { + float value = std::numeric_limits::infinity(); - propertyValues[i].push_back( value ); + auto tokenIndex = 5 + i; + if ( tokenIndex < tokens.size() ) + value = RiaStdStringTools::toDouble( tokens[tokenIndex] ); + + propertyValues[i].push_back( value ); + } + } } } else if ( firstToken.compare( "TRGL" ) == 0 ) { - int id1{ -1 }; - int id2{ -1 }; - int id3{ -1 }; - - lineStream >> id1 >> id2 >> id3; - - if ( id1 >= 0 && id2 >= 0 && id3 >= 0 ) + if ( tokens.size() > 3 ) { - trianglesByIds.emplace_back( static_cast( id1 ) ); - trianglesByIds.emplace_back( static_cast( id2 ) ); - trianglesByIds.emplace_back( static_cast( id3 ) ); + auto id1 = RiaStdStringTools::toInt( tokens[1] ); + auto id2 = RiaStdStringTools::toInt( tokens[2] ); + auto id3 = RiaStdStringTools::toInt( tokens[3] ); + + if ( id1 >= 0 && id2 >= 0 && id3 >= 0 ) + { + trianglesByIds.emplace_back( static_cast( id1 ) ); + trianglesByIds.emplace_back( static_cast( id2 ) ); + trianglesByIds.emplace_back( static_cast( id3 ) ); + } } } else if ( firstToken.compare( "END" ) == 0 ) @@ -153,19 +161,9 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g } else if ( firstToken.compare( "PROPERTIES" ) == 0 ) { - QString qstringLine = QString::fromStdString( line ); - - qstringLine.remove( "PROPERTIES" ); - -#if QT_VERSION >= QT_VERSION_CHECK( 5, 15, 0 ) - QStringList words = qstringLine.split( " ", Qt::SkipEmptyParts ); -#else - QStringList words = qstringLine.split( " ", QString::SkipEmptyParts ); -#endif - - for ( auto w : words ) + for ( size_t i = 1; i < tokens.size(); i++ ) { - propertyNames.push_back( w ); + propertyNames.push_back( QString::fromStdString( tokens[i] ) ); } propertyValues.resize( propertyNames.size() ); @@ -173,8 +171,8 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g else if ( firstToken.compare( "ZPOSITIVE" ) == 0 ) { std::string secondToken; - lineStream >> secondToken; + if ( tokens.size() > 1 ) secondToken = RiaStdStringTools::toUpper( tokens[1] ); if ( secondToken == "DEPTH" ) { zDir = GocadZPositive::Depth; @@ -368,8 +366,8 @@ std::pair, std::vector> auto to2d = []( const cvf::Vec3d vector ) -> cvf::Vec2d { return cvf::Vec2d( vector.x(), vector.y() ); }; auto to3d = []( const cvf::Vec2d vector ) -> cvf::Vec3d { return cvf::Vec3d( vector.x(), vector.y(), 0.0 ); }; - // Checks if the given vector is a possible new candidate for an axis vector and adds it to the given list of - // axesVectorCandidates. Also increases the number of occurrences of vector candidates. + // Checks if the given vector is a possible new candidate for an axis vector and adds it to the given list + // of axesVectorCandidates. Also increases the number of occurrences of vector candidates. auto maybeInsertAxisVectorCandidate = [epsilon]( const cvf::Vec2d vector, std::map& axesVectorCandidates, diff --git a/ApplicationLibCode/ProjectDataModel/Surfaces/RimEnsembleSurface.cpp b/ApplicationLibCode/ProjectDataModel/Surfaces/RimEnsembleSurface.cpp index 529a7b4347..2e430cc6ff 100644 --- a/ApplicationLibCode/ProjectDataModel/Surfaces/RimEnsembleSurface.cpp +++ b/ApplicationLibCode/ProjectDataModel/Surfaces/RimEnsembleSurface.cpp @@ -103,11 +103,14 @@ std::vector RimEnsembleSurface::sourceFileSurfaces() const //-------------------------------------------------------------------------------------------------- void RimEnsembleSurface::loadDataAndUpdate() { - for ( auto& w : sourceFileSurfaces() ) { - if ( !w->onLoadData() ) + auto fileSurfaces = sourceFileSurfaces(); + auto surfaceCount = static_cast( fileSurfaces.size() ); +#pragma omp parallel for + for ( int i = 0; i < surfaceCount; i++ ) { - RiaLogging::warning( QString( "Failed to load surface: %1" ).arg( w->surfaceFilePath() ) ); + auto surf = fileSurfaces[i]; + surf->onLoadData(); } } @@ -121,9 +124,16 @@ void RimEnsembleSurface::loadDataAndUpdate() { cvf::ref firstSurface = sourceSurfaceForStatistics[0]->surfaceData(); - std::vector> sourceSurfaces; - for ( auto& w : sourceSurfaceForStatistics ) - sourceSurfaces.push_back( RigSurfaceResampler::resampleSurface( firstSurface, w->surfaceData() ) ); + auto surfaceCount = static_cast( sourceSurfaceForStatistics.size() ); + std::vector> sourceSurfaces( surfaceCount ); + +#pragma omp parallel for + for ( int i = 0; i < surfaceCount; i++ ) + { + auto surf = sourceSurfaceForStatistics[i]; + auto resampledSurface = RigSurfaceResampler::resampleSurface( firstSurface, surf->surfaceData() ); + sourceSurfaces[i] = resampledSurface; + } m_statisticsSurface = RigSurfaceStatisticsCalculator::computeStatistics( sourceSurfaces ); if ( !m_statisticsSurface.isNull() ) diff --git a/ApplicationLibCode/ReservoirDataModel/RigSurface.cpp b/ApplicationLibCode/ReservoirDataModel/RigSurface.cpp index 790765ccd9..1ca3c64758 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigSurface.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigSurface.cpp @@ -104,6 +104,22 @@ void RigSurface::findIntersectingTriangles( const cvf::BoundingBox& inputBB, std m_surfaceBoundingBoxTree->findIntersections( inputBB, triangleStartIndices ); } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigSurface::maxExtentTriangleInXDirection() const +{ + return m_maxExtentTriangleInXDirection; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigSurface::maxExtentTriangleInYDirection() const +{ + return m_maxExtentTriangleInYDirection; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -118,6 +134,9 @@ void RigSurface::ensureIntersectionSearchTreeIsBuilt() cellBoundingBoxes.resize( itemCount ); boundingBoxIds.resize( itemCount ); + double maxX = -1.0; + double maxY = -1.0; + for ( size_t triangleIdx = 0; triangleIdx < itemCount; ++triangleIdx ) { cvf::BoundingBox& cellBB = cellBoundingBoxes[triangleIdx]; @@ -126,8 +145,14 @@ void RigSurface::ensureIntersectionSearchTreeIsBuilt() cellBB.add( m_vertices[m_triangleIndices[triangleIdx * 3 + 2]] ); boundingBoxIds[triangleIdx] = triangleIdx * 3; + + if ( cellBB.extent().x() > maxX ) maxX = cellBB.extent().x(); + if ( cellBB.extent().y() > maxY ) maxY = cellBB.extent().y(); } + m_maxExtentTriangleInXDirection = maxX; + m_maxExtentTriangleInYDirection = maxY; + m_surfaceBoundingBoxTree = new cvf::BoundingBoxTree; m_surfaceBoundingBoxTree->buildTreeFromBoundingBoxes( cellBoundingBoxes, &boundingBoxIds ); } diff --git a/ApplicationLibCode/ReservoirDataModel/RigSurface.h b/ApplicationLibCode/ReservoirDataModel/RigSurface.h index e9d0b0d44b..3d53081dcd 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigSurface.h +++ b/ApplicationLibCode/ReservoirDataModel/RigSurface.h @@ -49,6 +49,9 @@ public: void ensureIntersectionSearchTreeIsBuilt(); void findIntersectingTriangles( const cvf::BoundingBox& inputBB, std::vector* triangleStartIndices ) const; + double maxExtentTriangleInXDirection() const; + double maxExtentTriangleInYDirection() const; + private: size_t triangleCount() const; @@ -58,4 +61,7 @@ private: std::map> m_verticeResults; cvf::ref m_surfaceBoundingBoxTree; + + double m_maxExtentTriangleInXDirection; + double m_maxExtentTriangleInYDirection; }; diff --git a/ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.cpp b/ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.cpp index 7f3e3c79ac..fd8a180227 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.cpp @@ -91,9 +91,8 @@ bool RigSurfaceResampler::resamplePoint( RigSurface* surface, } } - // Handle cases where no match is found due to floating point imprecision, - // or when falling off resulting grid slightly. - double maxDistance = 10.0; + double maxDistance = computeMaxDistance( surface ); + return findClosestPointXY( pointAbove, vertices, maxDistance, intersectionPoint ); } @@ -106,24 +105,26 @@ bool RigSurfaceResampler::findClosestPointXY( const cvf::Vec3d& tar double maxDistance, cvf::Vec3d& intersectionPoint ) { + double maxDistanceSquared = maxDistance * maxDistance; + // Find closest vertices - double shortestDistance = std::numeric_limits::max(); - double closestZ = std::numeric_limits::infinity(); + double shortestDistanceSquared = std::numeric_limits::max(); + double closestZ = std::numeric_limits::infinity(); for ( auto v : vertices ) { // Ignore height (z) component when finding closest by // moving point to same height as target point above cvf::Vec3d p( v.x(), v.y(), targetPoint.z() ); - double distance = p.pointDistance( targetPoint ); - if ( distance < shortestDistance ) + double distanceSquared = p.pointDistanceSquared( targetPoint ); + if ( distanceSquared < shortestDistanceSquared ) { - shortestDistance = distance; - closestZ = v.z(); + shortestDistanceSquared = distanceSquared; + closestZ = v.z(); } } // Check if the closest point is not to far away to be valid - if ( shortestDistance < maxDistance ) + if ( shortestDistanceSquared < maxDistanceSquared ) { intersectionPoint = cvf::Vec3d( targetPoint.x(), targetPoint.y(), closestZ ); return true; @@ -133,3 +134,26 @@ bool RigSurfaceResampler::findClosestPointXY( const cvf::Vec3d& tar return false; } } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigSurfaceResampler::computeMaxDistance( RigSurface* surface ) +{ + // Handle cases where no match is found due to floating point imprecision, + // or when falling off resulting grid slightly. + // Use the XY extent of a triangle to define a suitable search distance + + const double minimumDistance = 10.0; + + if ( !surface ) return minimumDistance; + + auto maxX = surface->maxExtentTriangleInXDirection() / 2.0; + auto maxY = surface->maxExtentTriangleInYDirection() / 2.0; + + auto candidate = std::min( maxX, maxY ); + + double distance = std::max( minimumDistance, candidate ); + + return distance; +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.h b/ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.h index afe5391229..ceb1af445d 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.h +++ b/ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.h @@ -39,4 +39,6 @@ private: const std::vector& vertices, double maxDistance, cvf::Vec3d& intersectionPoint ); + + static double computeMaxDistance( RigSurface* surface ); };