From 5fa48d6e60f69ea8d2f7bcaa5efe16d5bf85fa0d Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Fri, 23 Oct 2026 20:57:31 +0200 Subject: [PATCH] VTK GeoMech: handle partial stress/strain tensor rows List only the S_VTK/E_VTK components that have data, and require the complete tensor for the derived S/E results. Guard against too few values in the tensor rows. Add tests for strain and partial data. --- .../GeoMechFileInterface/RifVtkReader.cpp | 59 +++++++++++----- .../UnitTests/RifVtkReader-Test.cpp | 69 +++++++++++++++++++ .../TestData/RifVtkReader/partial-00000.vtu | 49 +++++++++++++ .../TestData/RifVtkReader/partial.pvd | 6 ++ 4 files changed, 166 insertions(+), 17 deletions(-) create mode 100644 ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial-00000.vtu create mode 100644 ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial.pvd diff --git a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp index e6565d54bf..56395e8042 100644 --- a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp +++ b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp @@ -60,8 +60,29 @@ constexpr const char* E23 = "E23"; constexpr const char* derivedS = "S"; constexpr const char* derivedE = "E"; constexpr const char* por = "POR"; + +const std::vector stressComponents = { S11, S22, S33, S12, S13, S23 }; +const std::vector strainComponents = { E11, E22, E33, E12, E13, E23 }; } // namespace RifVtkFieldNames +namespace +{ +//-------------------------------------------------------------------------------------------------- +/// The components of a tensor field that have data, stored as "__" +//-------------------------------------------------------------------------------------------------- +std::vector availableComponents( const std::map>>>& data, + const std::string& fieldName, + const std::vector& components ) +{ + std::vector available; + for ( const auto& component : components ) + { + if ( data.contains( fieldName + RifVtkFieldNames::separator + component ) ) available.push_back( component ); + } + return available; +} +} // namespace + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -452,16 +473,15 @@ std::map> RifVtkReader::scalarElementField retVal[key] = {}; } - // Add S field with components if any S component data was loaded - if ( m_propertyPartDataElements.count( std::string( S ) + separator + S11 ) > 0 ) + // A file may contain only some of the tensor rows, so list only the components that have data + if ( auto components = availableComponents( m_propertyPartDataElements, S, stressComponents ); !components.empty() ) { - retVal[S] = { S11, S22, S33, S12, S13, S23 }; + retVal[S] = components; } - // Add E field with components if any E component data was loaded - if ( m_propertyPartDataElements.count( std::string( E ) + separator + E11 ) > 0 ) + if ( auto components = availableComponents( m_propertyPartDataElements, E, strainComponents ); !components.empty() ) { - retVal[E] = { E11, E22, E33, E12, E13, E23 }; + retVal[E] = components; } return retVal; @@ -475,13 +495,15 @@ std::map> RifVtkReader::scalarElementNodeF using namespace RifVtkFieldNames; std::map> retVal; - if ( m_propertyPartDataElements.contains( std::string( S ) + separator + S11 ) && m_propertyPartDataNodes.contains( por ) ) + // The derived results need the complete tensor + if ( availableComponents( m_propertyPartDataElements, S, stressComponents ).size() == stressComponents.size() && + m_propertyPartDataNodes.contains( por ) ) { - retVal[derivedS] = { S11, S22, S33, S12, S13, S23 }; + retVal[derivedS] = stressComponents; } - if ( m_propertyPartDataElements.contains( std::string( E ) + separator + E11 ) ) + if ( availableComponents( m_propertyPartDataElements, E, strainComponents ).size() == strainComponents.size() ) { - retVal[derivedE] = { E11, E22, E33, E12, E13, E23 }; + retVal[derivedE] = strainComponents; } return retVal; } @@ -584,14 +606,14 @@ void RifVtkReader::readElementField( const std::string& fieldName int frameIndex, std::vector*>* resultValues ) { - // S and E fields have 6 components stored as individual internal keys with a separator + // S and E fields have up to 6 components stored as individual internal keys with a separator. The result + // values follow the component list from scalarElementFieldAndComponentNames(), which has only those with data. if ( fieldName == RifVtkFieldNames::S || fieldName == RifVtkFieldNames::E ) { using namespace RifVtkFieldNames; - static const std::vector sComponents = { S11, S22, S33, S12, S13, S23 }; - static const std::vector eComponents = { E11, E22, E33, E12, E13, E23 }; - const auto& components = ( fieldName == RifVtkFieldNames::S ) ? sComponents : eComponents; - size_t count = std::min( components.size(), resultValues->size() ); + const auto components = + availableComponents( m_propertyPartDataElements, fieldName, fieldName == S ? stressComponents : strainComponents ); + size_t count = std::min( components.size(), resultValues->size() ); for ( size_t i = 0; i < count; i++ ) { std::string internalKey = fieldName + separator + components[i]; @@ -729,7 +751,8 @@ void RifVtkReader::readScalarData( RigFemPartCollection* if ( map == nullptr ) continue; // OPM Flow writes a symmetric 3x3 tensor as three rows "stress[0..2]" of 3 values each. - // The upper triangle is read. z is negated on import, so the xz and yz components change sign. + // The upper triangle is read. z is negated on import (RifVtkImportUtil::readPoints and + // readDisplacements), so the xz and yz components change sign. struct StressStrainMapping { std::string vtkName; @@ -763,7 +786,9 @@ void RifVtkReader::readScalarData( RigFemPartCollection* for ( int partIdx = 0; partIdx < femParts->partCount(); partIdx++ ) { - size_t dataSize = femParts->part( partIdx )->elementCount(); + size_t dataSize = femParts->part( partIdx )->elementCount(); + if ( values.size() < 3 * dataSize ) continue; + std::vector compValues( dataSize, 0.0 ); for ( size_t i = 0; i < dataSize; i++ ) { diff --git a/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp b/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp index f9e2704ad7..a09bbf8e28 100644 --- a/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp +++ b/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp @@ -243,3 +243,72 @@ TEST( RifVtkReaderTest, ReadElementNodalStressFromCellsAndPor ) EXPECT_FLOAT_EQ( e[c][idx], expectedE[c] ); } } + +//-------------------------------------------------------------------------------------------------- +/// strain[0..2] are the rows of a symmetric tensor; xz and yz change sign with the negated z +//-------------------------------------------------------------------------------------------------- +TEST( RifVtkReaderTest, ReadElementFieldStrainTensorRows ) +{ + RifVtkReader reader; + reader.openFile( ( testDataDir() / "model.pvd" ).string(), nullptr ); + + cvf::ref femParts = new RigFemPartCollection(); + reader.readFemParts( femParts.p() ); + + auto elementFields = reader.scalarElementFieldAndComponentNames(); + ASSERT_TRUE( elementFields.count( "E_VTK" ) > 0 ); + const std::vector expectedComponents = { "E11", "E22", "E33", "E12", "E13", "E23" }; + EXPECT_EQ( elementFields["E_VTK"], expectedComponents ); + + std::vector> values( 6 ); + std::vector*> resultValues; + for ( auto& v : values ) + resultValues.push_back( &v ); + + reader.readElementField( "E_VTK", 0, 1, 0, &resultValues ); + + const std::vector expectedCell0 = { 0.1f, 0.4f, 0.6f, 0.2f, -0.3f, -0.5f }; + const std::vector expectedCell1 = { 1.1f, 1.4f, 1.6f, 1.2f, -1.3f, -1.5f }; + for ( size_t c = 0; c < 6; c++ ) + { + ASSERT_EQ( values[c].size(), 2u ); + EXPECT_FLOAT_EQ( values[c][0], expectedCell0[c] ); + EXPECT_FLOAT_EQ( values[c][1], expectedCell1[c] ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// Only stress[0] and strain[1..2]: list only the components with data, and no derived results +//-------------------------------------------------------------------------------------------------- +TEST( RifVtkReaderTest, ReadElementFieldPartialTensorRows ) +{ + RifVtkReader reader; + reader.openFile( ( testDataDir() / "partial.pvd" ).string(), nullptr ); + + cvf::ref femParts = new RigFemPartCollection(); + reader.readFemParts( femParts.p() ); + + auto elementFields = reader.scalarElementFieldAndComponentNames(); + const std::vector expectedS = { "S11", "S12", "S13" }; + const std::vector expectedE = { "E22", "E33", "E23" }; + EXPECT_EQ( elementFields["S_VTK"], expectedS ); + EXPECT_EQ( elementFields["E_VTK"], expectedE ); + + EXPECT_TRUE( reader.scalarElementNodeFieldAndComponentNames().empty() ); + + std::vector> s( 3 ), e( 3 ); + std::vector*> sValues = { &s[0], &s[1], &s[2] }; + std::vector*> eValues = { &e[0], &e[1], &e[2] }; + reader.readElementField( "S_VTK", 0, 0, 0, &sValues ); + reader.readElementField( "E_VTK", 0, 0, 0, &eValues ); + + const std::vector expectedSCell1 = { 111, 112, -113 }; + const std::vector expectedECell1 = { 1.22f, 1.33f, -1.23f }; + for ( size_t c = 0; c < 3; c++ ) + { + ASSERT_EQ( s[c].size(), 2u ); + ASSERT_EQ( e[c].size(), 2u ); + EXPECT_FLOAT_EQ( s[c][1], expectedSCell1[c] ); + EXPECT_FLOAT_EQ( e[c][1], expectedECell1[c] ); + } +} diff --git a/ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial-00000.vtu b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial-00000.vtu new file mode 100644 index 0000000000..75639e9d78 --- /dev/null +++ b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial-00000.vtu @@ -0,0 +1,49 @@ + + + + + + + 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 + + + + + 100.0 200.0 + + + 11 12 13 111 112 113 + + + 0.12 0.22 0.23 1.12 1.22 1.23 + + + 0.13 0.23 0.33 1.13 1.23 1.33 + + + + + 0 0 0 1 0 0 2 0 0 + 0 1 0 1 1 0 2 1 0 + 0 0 1 1 0 1 2 0 1 + 0 1 1 1 1 1 2 1 1 + + + + + 0 1 4 3 6 7 10 9 + 1 2 5 4 7 8 11 10 + + + 8 16 + + + 12 12 + + + + + diff --git a/ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial.pvd b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial.pvd new file mode 100644 index 0000000000..fc44613c64 --- /dev/null +++ b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial.pvd @@ -0,0 +1,6 @@ + + + + + +