diff --git a/ApplicationLibCode/ProjectDataModel/RimEclipseCase.cpp b/ApplicationLibCode/ProjectDataModel/RimEclipseCase.cpp index 6ecbd61334..b8e456e07f 100644 --- a/ApplicationLibCode/ProjectDataModel/RimEclipseCase.cpp +++ b/ApplicationLibCode/ProjectDataModel/RimEclipseCase.cpp @@ -979,6 +979,27 @@ bool RimEclipseCase::openReserviorCase() return false; } + if ( eclipseCaseData() && eclipseCaseData()->mainGrid() && + !eclipseCaseData()->mainGrid()->hasValidCharacteristicCellSizes() ) + { + RigMainGrid* mainGrid = eclipseCaseData()->mainGrid(); + + auto activeCellInfo = eclipseCaseData()->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL ); + if ( activeCellInfo ) + { + std::vector reservoirCellIndices; + for ( size_t i = 0; i < mainGrid->cellCount(); i++ ) + { + if ( activeCellInfo->isActive( i ) ) + { + reservoirCellIndices.push_back( i ); + } + } + mainGrid->computeCharacteristicCellSize( reservoirCellIndices ); + mainGrid->computeFaceNormalsDirection( reservoirCellIndices ); + } + } + bool createPlaceholderEntries = true; if ( dynamic_cast( this ) ) { diff --git a/ApplicationLibCode/ReservoirDataModel/RigMainGrid.cpp b/ApplicationLibCode/ReservoirDataModel/RigMainGrid.cpp index f7c16d60fa..9a75358653 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigMainGrid.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigMainGrid.cpp @@ -46,6 +46,9 @@ RigMainGrid::RigMainGrid() m_useMapAxes = false; m_mapAxes = defaultMapAxes(); m_dualPorosity = false; + + m_isFaceNormalsOutwards = true; + m_isFaceNormalsOutwardsComputed = false; } RigMainGrid::~RigMainGrid() @@ -650,6 +653,23 @@ void RigMainGrid::distributeNNCsToFaults() /// but if (only) one of the flipX/Y is done, the cell is back to normal //-------------------------------------------------------------------------------------------------- bool RigMainGrid::isFaceNormalsOutwards() const +{ + if ( !m_isFaceNormalsOutwardsComputed ) + { + std::vector reservoirCellIndices; + reservoirCellIndices.resize( cellCount() ); + std::iota( reservoirCellIndices.begin(), reservoirCellIndices.end(), 0 ); + + computeFaceNormalsDirection( reservoirCellIndices ); + } + + return m_isFaceNormalsOutwards; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigMainGrid::computeFaceNormalsDirection( const std::vector& reservoirCellIndices ) const { auto isValidAndFaceNormalDir = []( const double ijSize, const double kSize, @@ -679,8 +699,9 @@ bool RigMainGrid::isFaceNormalsOutwards() const characteristicCellSizes( &iSize, &jSize, &kSize ); const double characteristicVolume = iSize * jSize * kSize; - for ( const auto& cell : m_cells ) + for ( const auto& index : reservoirCellIndices ) { + const auto& cell = m_cells[index]; if ( !cell.isInvalid() ) { // Some cells can be very twisted and distorted. Use a volume criteria to find a reasonably regular cell. @@ -701,21 +722,25 @@ bool RigMainGrid::isFaceNormalsOutwards() const if ( direction1 && direction2 && direction3 && direction4 ) { // All face normals pointing outwards - return true; + m_isFaceNormalsOutwards = true; + m_isFaceNormalsOutwardsComputed = true; + return; } if ( !direction1 && !direction2 && !direction3 && !direction4 ) { // All cell face normals pointing inwards - return false; + m_isFaceNormalsOutwards = false; + m_isFaceNormalsOutwardsComputed = true; + return; } - - // This is a mixed case, where some faces are outwards and some are inwards - continue; } } - return false; + // If this code is reached, it was not possible to get a consistent answer on the direction of a cell surface + // normal. Set a default direction for face normals. + m_isFaceNormalsOutwards = true; + m_isFaceNormalsOutwardsComputed = true; } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ReservoirDataModel/RigMainGrid.h b/ApplicationLibCode/ReservoirDataModel/RigMainGrid.h index 2f29de4920..33b7abc61d 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigMainGrid.h +++ b/ApplicationLibCode/ReservoirDataModel/RigMainGrid.h @@ -87,6 +87,7 @@ public: const RigFault* findFaultFromCellIndexAndCellFace( size_t reservoirCellIndex, cvf::StructGridInterface::FaceType face ) const; bool isFaceNormalsOutwards() const; + void computeFaceNormalsDirection( const std::vector& reservoirCellIndices ) const; void computeCachedData( std::string* aabbTreeInfo = nullptr ); void initAllSubGridsParentGridPointer(); @@ -142,4 +143,7 @@ private: std::array m_mapAxes; bool m_dualPorosity; + + mutable bool m_isFaceNormalsOutwards; + mutable bool m_isFaceNormalsOutwardsComputed; }; diff --git a/Fwk/AppFwk/CommonCode/cvfStructGrid.cpp b/Fwk/AppFwk/CommonCode/cvfStructGrid.cpp index fb7d30d48b..bbd125ee35 100644 --- a/Fwk/AppFwk/CommonCode/cvfStructGrid.cpp +++ b/Fwk/AppFwk/CommonCode/cvfStructGrid.cpp @@ -38,6 +38,8 @@ #include "cvfBase.h" #include "cvfBoundingBox.h" +#include + namespace caf { template <> @@ -319,99 +321,13 @@ void StructGridInterface::characteristicCellSizes( double* iSize, double* jSize, { CVF_ASSERT( iSize && jSize && kSize ); - if ( m_characteristicCellSizeI == cvf::UNDEFINED_DOUBLE || m_characteristicCellSizeJ == cvf::UNDEFINED_DOUBLE || - m_characteristicCellSizeK == cvf::UNDEFINED_DOUBLE ) + if ( !hasValidCharacteristicCellSizes() ) { - ubyte faceConnPosI[4]; - cellFaceVertexIndices( StructGridInterface::POS_I, faceConnPosI ); + std::vector reservoirCellIndices; + reservoirCellIndices.resize( cellCountI() * cellCountJ() * cellCountK() ); + std::iota( reservoirCellIndices.begin(), reservoirCellIndices.end(), 0 ); - ubyte faceConnNegI[4]; - cellFaceVertexIndices( StructGridInterface::NEG_I, faceConnNegI ); - - ubyte faceConnPosJ[4]; - cellFaceVertexIndices( StructGridInterface::POS_J, faceConnPosJ ); - - ubyte faceConnNegJ[4]; - cellFaceVertexIndices( StructGridInterface::NEG_J, faceConnNegJ ); - - ubyte faceConnPosK[4]; - cellFaceVertexIndices( StructGridInterface::POS_K, faceConnPosK ); - - ubyte faceConnNegK[4]; - cellFaceVertexIndices( StructGridInterface::NEG_K, faceConnNegK ); - - double iLengthAccumulated = 0.0; - double jLengthAccumulated = 0.0; - double kLengthAccumulated = 0.0; - - cvf::Vec3d cornerVerts[8]; - size_t cellCount = 0; - - size_t k; - for ( k = 0; k < cellCountK(); k++ ) - { - size_t j; - for ( j = 0; j < cellCountJ(); j++ ) - { - size_t i; - for ( i = 0; i < cellCountI(); i += 10 ) // NB! Evaluate every n-th cell - { - if ( isCellValid( i, j, k ) ) - { - size_t cellIndex = cellIndexFromIJK( i, j, k ); - cellCornerVertices( cellIndex, cornerVerts ); - - cvf::BoundingBox bb; - for ( const auto& v : cornerVerts ) - { - bb.add( v ); - } - - // Exclude cells with very small volumes - const double tolerance = 0.2; - if ( bb.extent().z() < tolerance ) continue; - - iLengthAccumulated += - ( cornerVerts[faceConnPosI[0]] - cornerVerts[faceConnNegI[0]] ).lengthSquared(); - iLengthAccumulated += - ( cornerVerts[faceConnPosI[1]] - cornerVerts[faceConnNegI[3]] ).lengthSquared(); - iLengthAccumulated += - ( cornerVerts[faceConnPosI[2]] - cornerVerts[faceConnNegI[2]] ).lengthSquared(); - iLengthAccumulated += - ( cornerVerts[faceConnPosI[3]] - cornerVerts[faceConnNegI[1]] ).lengthSquared(); - - jLengthAccumulated += - ( cornerVerts[faceConnPosJ[0]] - cornerVerts[faceConnNegJ[0]] ).lengthSquared(); - jLengthAccumulated += - ( cornerVerts[faceConnPosJ[1]] - cornerVerts[faceConnNegJ[3]] ).lengthSquared(); - jLengthAccumulated += - ( cornerVerts[faceConnPosJ[2]] - cornerVerts[faceConnNegJ[2]] ).lengthSquared(); - jLengthAccumulated += - ( cornerVerts[faceConnPosJ[3]] - cornerVerts[faceConnNegJ[1]] ).lengthSquared(); - - kLengthAccumulated += - ( cornerVerts[faceConnPosK[0]] - cornerVerts[faceConnNegK[0]] ).lengthSquared(); - kLengthAccumulated += - ( cornerVerts[faceConnPosK[1]] - cornerVerts[faceConnNegK[3]] ).lengthSquared(); - kLengthAccumulated += - ( cornerVerts[faceConnPosK[2]] - cornerVerts[faceConnNegK[2]] ).lengthSquared(); - kLengthAccumulated += - ( cornerVerts[faceConnPosK[3]] - cornerVerts[faceConnNegK[1]] ).lengthSquared(); - - cellCount++; - } - } - } - } - - double divisor = cellCount * 4.0; - - if ( divisor > 0.0 ) - { - m_characteristicCellSizeI = cvf::Math::sqrt( iLengthAccumulated / divisor ); - m_characteristicCellSizeJ = cvf::Math::sqrt( jLengthAccumulated / divisor ); - m_characteristicCellSizeK = cvf::Math::sqrt( kLengthAccumulated / divisor ); - } + computeCharacteristicCellSize( reservoirCellIndices ); } *iSize = m_characteristicCellSizeI; @@ -419,4 +335,100 @@ void StructGridInterface::characteristicCellSizes( double* iSize, double* jSize, *kSize = m_characteristicCellSizeK; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool StructGridInterface::hasValidCharacteristicCellSizes() const +{ + if ( m_characteristicCellSizeI == cvf::UNDEFINED_DOUBLE || m_characteristicCellSizeJ == cvf::UNDEFINED_DOUBLE || + m_characteristicCellSizeK == cvf::UNDEFINED_DOUBLE ) + return false; + + return true; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void StructGridInterface::computeCharacteristicCellSize( const std::vector& globalCellIndices ) const +{ + ubyte faceConnPosI[4]; + cellFaceVertexIndices( StructGridInterface::POS_I, faceConnPosI ); + + ubyte faceConnNegI[4]; + cellFaceVertexIndices( StructGridInterface::NEG_I, faceConnNegI ); + + ubyte faceConnPosJ[4]; + cellFaceVertexIndices( StructGridInterface::POS_J, faceConnPosJ ); + + ubyte faceConnNegJ[4]; + cellFaceVertexIndices( StructGridInterface::NEG_J, faceConnNegJ ); + + ubyte faceConnPosK[4]; + cellFaceVertexIndices( StructGridInterface::POS_K, faceConnPosK ); + + ubyte faceConnNegK[4]; + cellFaceVertexIndices( StructGridInterface::NEG_K, faceConnNegK ); + + double iLengthAccumulated = 0.0; + double jLengthAccumulated = 0.0; + double kLengthAccumulated = 0.0; + + cvf::Vec3d cornerVerts[8]; + size_t evaluatedCellCount = 0; + + // Evaluate N-th cells, compute the stride between each index + size_t stride = std::max( size_t( 1 ), globalCellIndices.size() / 100 ); + + size_t i, j, k = 0; + size_t index = 0; + while ( index < globalCellIndices.size() - 1 ) + { + size_t cellIndex = globalCellIndices[index]; + ijkFromCellIndex( cellIndex, &i, &j, &k ); + if ( isCellValid( i, j, k ) ) + { + cellCornerVertices( cellIndex, cornerVerts ); + + cvf::BoundingBox bb; + for ( const auto& v : cornerVerts ) + { + bb.add( v ); + } + + // Exclude cells with very small volumes + const double tolerance = 0.2; + if ( bb.extent().z() < tolerance ) continue; + + iLengthAccumulated += ( cornerVerts[faceConnPosI[0]] - cornerVerts[faceConnNegI[0]] ).lengthSquared(); + iLengthAccumulated += ( cornerVerts[faceConnPosI[1]] - cornerVerts[faceConnNegI[3]] ).lengthSquared(); + iLengthAccumulated += ( cornerVerts[faceConnPosI[2]] - cornerVerts[faceConnNegI[2]] ).lengthSquared(); + iLengthAccumulated += ( cornerVerts[faceConnPosI[3]] - cornerVerts[faceConnNegI[1]] ).lengthSquared(); + + jLengthAccumulated += ( cornerVerts[faceConnPosJ[0]] - cornerVerts[faceConnNegJ[0]] ).lengthSquared(); + jLengthAccumulated += ( cornerVerts[faceConnPosJ[1]] - cornerVerts[faceConnNegJ[3]] ).lengthSquared(); + jLengthAccumulated += ( cornerVerts[faceConnPosJ[2]] - cornerVerts[faceConnNegJ[2]] ).lengthSquared(); + jLengthAccumulated += ( cornerVerts[faceConnPosJ[3]] - cornerVerts[faceConnNegJ[1]] ).lengthSquared(); + + kLengthAccumulated += ( cornerVerts[faceConnPosK[0]] - cornerVerts[faceConnNegK[0]] ).lengthSquared(); + kLengthAccumulated += ( cornerVerts[faceConnPosK[1]] - cornerVerts[faceConnNegK[3]] ).lengthSquared(); + kLengthAccumulated += ( cornerVerts[faceConnPosK[2]] - cornerVerts[faceConnNegK[2]] ).lengthSquared(); + kLengthAccumulated += ( cornerVerts[faceConnPosK[3]] - cornerVerts[faceConnNegK[1]] ).lengthSquared(); + + evaluatedCellCount++; + } + + index += stride; + } + + double divisor = evaluatedCellCount * 4.0; + + if ( divisor > 0.0 ) + { + m_characteristicCellSizeI = cvf::Math::sqrt( iLengthAccumulated / divisor ); + m_characteristicCellSizeJ = cvf::Math::sqrt( jLengthAccumulated / divisor ); + m_characteristicCellSizeK = cvf::Math::sqrt( kLengthAccumulated / divisor ); + } +} + } // namespace cvf diff --git a/Fwk/AppFwk/CommonCode/cvfStructGrid.h b/Fwk/AppFwk/CommonCode/cvfStructGrid.h index afd21a189b..e82be4f043 100644 --- a/Fwk/AppFwk/CommonCode/cvfStructGrid.h +++ b/Fwk/AppFwk/CommonCode/cvfStructGrid.h @@ -91,6 +91,9 @@ public: virtual cvf::Vec3d maxCoordinate() const = 0; void characteristicCellSizes( double* iSize, double* jSize, double* kSize ) const; + bool hasValidCharacteristicCellSizes() const; + void computeCharacteristicCellSize( const std::vector& globalCellIndices ) const; + virtual cvf::Vec3d displayModelOffset() const; virtual bool cellIJKNeighbor( size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex ) const = 0;