Ensemble Surface improvements

* Performance : Improve surface import
* Performance: Use opm when importing files
* Surface : Use the triangle size as basis for the maximum search distance
* Performance : Resample surfaces in parallell
* Performance: Import file surfaces in parallell
* Ensemble Surface : Create one ensemble per surface

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<QString, QStringList> keyFileComponentsForAllFiles =
         RiaFilePathTools::keyPathComponentsForEachFilePath( allFileNames );
 
-    QStringList keyFileComponents = keyFileComponentsForAllFiles[sourceFileName];
+    return uniqueShortNameFromComponents( sourceFileName, keyFileComponentsForAllFiles, ensembleCaseName );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+QString RiaEnsembleNameTools::uniqueShortNameFromComponents( const QString&                        sourceFileName,
+                                                             const std::map<QString, QStringList>& 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() )

ApplicationLibCode/Application/Tools/RiaEnsembleNameTools.h

@@ -19,6 +19,7 @@
 #pragma once
 
 #include <QString>
+#include <map>
 
 //==================================================================================================
 //
@@ -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<QString, QStringList>& keyFileComponentsForAllFiles,
+                                                  const QString&                        ensembleCaseName );
 };

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<QString>&         fileNames )
+                                                      const std::vector<QString>&       fileNames )
 {
     std::vector layers = ui.layers();
 
     caf::ProgressInfo progress( fileNames.size(), "Generating ensemble surface" );
 
     QStringList allSurfaceFileNames;
-    for ( auto fileName : fileNames )
+
+    auto fileCount = static_cast<int>( fileNames.size() );
+#pragma omp parallel for
+    for ( int i = 0; i < fileCount; i++ )
+    {
+        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 ) );
-        auto [isOk, surfaceFileNames] = RimcCommandRouter_extractSurfaces::extractSurfaces( fileName, layers );
             if ( isOk ) allSurfaceFileNames << surfaceFileNames;
         }
+    }
+    progress.setProgress( fileNames.size() );
 
     if ( ui.autoCreateEnsembleSurfaces() )
         RicImportEnsembleSurfaceFeature::importEnsembleSurfaceFromFiles( allSurfaceFileNames );

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<QString>& fileNames );
+    static void executeCommand( const RicCreateEnsembleSurfaceUi& ui, const std::vector<QString>& fileNames );
 
 protected:
     bool isCommandEnabled() override;

ApplicationLibCode/Commands/RicImportEnsembleSurfaceFeature.cpp

@@ -20,6 +20,7 @@
 
 #include "RiaApplication.h"
 #include "RiaEnsembleNameTools.h"
+#include "RiaFilePathTools.h"
 #include "RiaLogging.h"
 #include "RiaSummaryTools.h"
 
@@ -75,8 +76,23 @@ void RicImportEnsembleSurfaceFeature::importEnsembleSurfaceFromFiles( const QStr
 {
     if ( fileNames.isEmpty() ) return;
 
+    // Create a list of file names for each layer
+    std::map<QString, QStringList> fileNamesForEachLayer;
+    for ( const auto& name : fileNames )
+    {
+        QFileInfo fi( name );
+
+        auto layerName = fi.baseName();
+        fileNamesForEachLayer[layerName].push_back( name );
+    }
+
+    RimEnsembleSurface* ensembleToSelect = nullptr;
+    for ( const auto& fileNamesForLayer : fileNamesForEachLayer )
+    {
+        auto    filenames    = fileNamesForLayer.second;
         QString ensembleName = RiaEnsembleNameTools::findSuitableEnsembleName( fileNames );
-    QString layerName    = RiaEnsembleNameTools::findCommonBaseName( fileNames );
+
+        QString layerName = fileNamesForLayer.first;
         if ( !layerName.isEmpty() )
         {
             ensembleName += QString( " : %1" ).arg( layerName );
@@ -84,24 +100,50 @@ void RicImportEnsembleSurfaceFeature::importEnsembleSurfaceFromFiles( const QStr
 
         if ( ensembleName.isEmpty() ) ensembleName = "Ensemble Surface";
 
+        std::map<QString, QStringList> keyFileComponentsForAllFiles =
+            RiaFilePathTools::keyPathComponentsForEachFilePath( fileNames );
+
         std::vector<RimFileSurface*> surfaces;
-    for ( const auto& fileName : fileNames )
+        for ( size_t i = 0; i < fileNames.size(); i++ )
         {
-        RimFileSurface* fileSurface = new RimFileSurface;
+            surfaces.push_back( new RimFileSurface );
+        }
+
+        auto fileCount = static_cast<int>( 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::uniqueShortName( fileName, fileNames );
+            auto shortName =
+                RiaEnsembleNameTools::uniqueShortNameFromComponents( fileName, keyFileComponentsForAllFiles, ensembleName );
             fileSurface->setUserDescription( shortName );
 
-        if ( fileSurface->onLoadData() )
+            auto isOk = fileSurface->onLoadData();
+            if ( !isOk )
             {
-            surfaces.push_back( fileSurface );
+                delete fileSurface;
+                surfaces[i] = nullptr;
             }
         }
 
+        {
+            // Remove null pointers from vector of surfaces
+            std::vector<RimFileSurface*> tmp;
+            for ( auto s : surfaces )
+            {
+                if ( s != nullptr ) tmp.push_back( s );
+            }
+
+            surfaces.swap( tmp );
+        }
+
         if ( surfaces.empty() ) return;
 
-    RimEnsembleSurface* ensemble = new RimEnsembleSurface;
+        auto ensemble = new RimEnsembleSurface;
         ensemble->setCollectionName( ensembleName );
         for ( auto surface : surfaces )
             ensemble->addFileSurface( surface );
@@ -109,8 +151,11 @@ void RicImportEnsembleSurfaceFeature::importEnsembleSurfaceFromFiles( const QStr
         ensemble->loadDataAndUpdate();
         RimProject::current()->activeOilField()->surfaceCollection->addEnsembleSurface( ensemble );
 
+        ensembleToSelect = ensemble;
+    }
+
     RimProject::current()->activeOilField()->surfaceCollection->updateConnectedEditors();
-    Riu3DMainWindowTools::selectAsCurrentItem( ensemble );
+    Riu3DMainWindowTools::selectAsCurrentItem( ensembleToSelect );
 }
 
 //--------------------------------------------------------------------------------------------------

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<std::vector<float>> 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,24 +77,21 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g
             std::string line;
             std::getline( stream, line );
 
-            std::transform( line.begin(), line.end(), line.begin(), ::toupper );
+            auto tokens = RiaStdStringTools::splitString( line, ' ' );
-
-            std::istringstream lineStream( line );
 
             std::string firstToken;
-            lineStream >> firstToken;
+            if ( !tokens.empty() ) firstToken = tokens.front();
 
             if ( isInTfaceSection )
             {
                 if ( firstToken.compare( "VRTX" ) == 0 )
                 {
-                    int         vertexId = -1;
+                    if ( tokens.size() > 4 )
-                    double      x{ std::numeric_limits<double>::infinity() };
+                    {
-                    double      y{ std::numeric_limits<double>::infinity() };
+                        int    vertexId = RiaStdStringTools::toInt( tokens[1] );
-                    double      z{ std::numeric_limits<double>::infinity() };
+                        double x        = RiaStdStringTools::toDouble( tokens[2] );
-                    std::string endVertex;
+                        double y        = RiaStdStringTools::toDouble( tokens[3] );
-
+                        double z        = RiaStdStringTools::toDouble( tokens[4] );
-                    lineStream >> vertexId >> x >> y >> z >> endVertex;
 
                         if ( vertexId > -1 )
                         {
@@ -101,14 +104,15 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g
                             vertexIdToIndex[vertexId] = static_cast<unsigned>( vertices.size() - 1 );
                         }
                     }
+                }
                 else if ( firstToken.compare( "PVRTX" ) == 0 )
                 {
-                    int    vertexId = -1;
+                    if ( tokens.size() > 4 )
-                    double x{ std::numeric_limits<double>::infinity() };
+                    {
-                    double y{ std::numeric_limits<double>::infinity() };
+                        int    vertexId = RiaStdStringTools::toInt( tokens[1] );
-                    double z{ std::numeric_limits<double>::infinity() };
+                        double x        = RiaStdStringTools::toDouble( tokens[2] );
-
+                        double y        = RiaStdStringTools::toDouble( tokens[3] );
-                    lineStream >> vertexId >> x >> y >> z;
+                        double z        = RiaStdStringTools::toDouble( tokens[4] );
 
                         if ( vertexId > -1 )
                         {
@@ -116,24 +120,27 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g
 
                             vertices.emplace_back( cvf::Vec3d( x, y, z ) );
                             vertexIdToIndex[vertexId] = static_cast<unsigned>( vertices.size() - 1 );
-                    }
 
                             for ( size_t i = 0; i < propertyNames.size(); i++ )
                             {
                                 float value = std::numeric_limits<double>::infinity();
 
-                        lineStream >> 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 };
+                    if ( tokens.size() > 3 )
-                    int id2{ -1 };
+                    {
-                    int id3{ -1 };
+                        auto id1 = RiaStdStringTools::toInt( tokens[1] );
-
+                        auto id2 = RiaStdStringTools::toInt( tokens[2] );
-                    lineStream >> id1 >> id2 >> id3;
+                        auto id3 = RiaStdStringTools::toInt( tokens[3] );
 
                         if ( id1 >= 0 && id2 >= 0 && id3 >= 0 )
                         {
@@ -142,6 +149,7 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g
                             trianglesByIds.emplace_back( static_cast<unsigned int>( id3 ) );
                         }
                     }
+                }
                 else if ( firstToken.compare( "END" ) == 0 )
                 {
                     isInTfaceSection = false;
@@ -153,19 +161,9 @@ void RifSurfaceImporter::readGocadFile( const QString& filename, RigGocadData* g
             }
             else if ( firstToken.compare( "PROPERTIES" ) == 0 )
             {
-                QString qstringLine = QString::fromStdString( line );
+                for ( size_t i = 1; i < tokens.size(); i++ )
-
-                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 )
                 {
-                    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<cvf::Vec3d>, std::vector<unsigned>>
     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
+    // Checks if the given vector is a possible new candidate for an axis vector and adds it to the given list
-    // axesVectorCandidates. Also increases the number of occurrences of vector candidates.
+    // of axesVectorCandidates. Also increases the number of occurrences of vector candidates.
     auto maybeInsertAxisVectorCandidate =
         [epsilon]( const cvf::Vec2d                              vector,
                    std::map<cvf::Vec2d, double, vec2dCompare>&   axesVectorCandidates,

ApplicationLibCode/ProjectDataModel/Surfaces/RimEnsembleSurface.cpp

@@ -103,11 +103,14 @@ std::vector<RimFileSurface*> RimEnsembleSurface::sourceFileSurfaces() const
 //--------------------------------------------------------------------------------------------------
 void RimEnsembleSurface::loadDataAndUpdate()
 {
-    for ( auto& w : sourceFileSurfaces() )
     {
-        if ( !w->onLoadData() )
+        auto fileSurfaces = sourceFileSurfaces();
+        auto surfaceCount = static_cast<int>( 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<RigSurface> firstSurface = sourceSurfaceForStatistics[0]->surfaceData();
 
-        std::vector<cvf::ref<RigSurface>> sourceSurfaces;
+        auto                              surfaceCount = static_cast<int>( sourceSurfaceForStatistics.size() );
-        for ( auto& w : sourceSurfaceForStatistics )
+        std::vector<cvf::ref<RigSurface>> sourceSurfaces( surfaceCount );
-            sourceSurfaces.push_back( RigSurfaceResampler::resampleSurface( firstSurface, w->surfaceData() ) );
+
+#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() )

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 );
     }

ApplicationLibCode/ReservoirDataModel/RigSurface.h

@@ -49,6 +49,9 @@ public:
     void ensureIntersectionSearchTreeIsBuilt();
     void findIntersectingTriangles( const cvf::BoundingBox& inputBB, std::vector<size_t>* triangleStartIndices ) const;
 
+    double maxExtentTriangleInXDirection() const;
+    double maxExtentTriangleInYDirection() const;
+
 private:
     size_t triangleCount() const;
 
@@ -58,4 +61,7 @@ private:
     std::map<QString, std::vector<float>> m_verticeResults;
 
     cvf::ref<cvf::BoundingBoxTree> m_surfaceBoundingBoxTree;
+
+    double m_maxExtentTriangleInXDirection;
+    double m_maxExtentTriangleInYDirection;
 };

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,
+    double maxDistance = computeMaxDistance( surface );
-    // or when falling off resulting grid slightly.
+
-    double maxDistance = 10.0;
     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<double>::max();
+    double shortestDistanceSquared = std::numeric_limits<double>::max();
     double closestZ                = std::numeric_limits<double>::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 );
+        double     distanceSquared = p.pointDistanceSquared( targetPoint );
-        if ( distance < shortestDistance )
+        if ( distanceSquared < shortestDistanceSquared )
         {
-            shortestDistance = distance;
+            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;
+}

ApplicationLibCode/ReservoirDataModel/RigSurfaceResampler.h

@@ -39,4 +39,6 @@ private:
                                     const std::vector<cvf::Vec3d>& vertices,
                                     double                         maxDistance,
                                     cvf::Vec3d&                    intersectionPoint );
+
+    static double computeMaxDistance( RigSurface* surface );
 };