From 03eff53507cfa1b5a34caf316c052373c8999f48 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 14 Oct 2026 15:26:10 +0200 Subject: [PATCH] VTK GeoMech: read stress/strain as the rows of a 3x3 tensor OPM Flow writes each symmetric tensor as three cell arrays "stress[0]", "stress[1]", "stress[2]", the rows (Si1, Si2, Si3). The reader took "stress[0]" as (S11, S22, S33) and "stress[1]" as (S12, S13, S23), so S22, S33 and S13 held the wrong components and "stress[2]" was ignored. Read the upper triangle from the three rows. Coordinates are imported with z negated; S13/S23 and E13/E23 now change sign with it, as U3 already does. Co-Authored-By: Claude Opus 5.5 --- .../GeoMechFileInterface/RifVtkReader.cpp | 56 +++++++++---------- .../UnitTests/RifVtkReader-Test.cpp | 33 +++++++++++ .../TestData/RifVtkReader/model-00001.vtu | 9 +++ 3 files changed, 69 insertions(+), 29 deletions(-) diff --git a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp index c11f188591..caa37d96a6 100644 --- a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp +++ b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp @@ -664,50 +664,48 @@ void RifVtkReader::readScalarData( RigFemPartCollection* auto map = propertyDataMap( RigFemResultPosEnum::RIG_ELEMENT ); if ( map == nullptr ) continue; - // Multi-component element results: stress[0/1] -> S_VTK, strain[0/1] -> E_VTK - // Each VTK array contains 3 components per element. - // Mapping: - // stress[0]: S_VTK.S11 (index 0), S_VTK.S22 (index 1), S_VTK.S33 (index 2) - // stress[1]: S_VTK.S12 (index 0), S_VTK.S13 (index 1), S_VTK.S23 (index 2) - // strain[0]: E_VTK.E11 (index 0), E_VTK.E22 (index 1), E_VTK.E33 (index 2) - // strain[1]: E_VTK.E12 (index 0), E_VTK.E13 (index 1), E_VTK.E23 (index 2) - + // 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. struct StressStrainMapping { - std::string vtkName; - std::string fieldName; - std::array componentNames; + std::string vtkName; + int column; + std::string fieldName; + std::string componentName; + double sign; }; using namespace RifVtkFieldNames; static const std::vector mappings = { - { "stress[0]", S, { S11, S22, S33 } }, - { "stress[1]", S, { S12, S13, S23 } }, - { "strain[0]", E, { E11, E22, E33 } }, - { "strain[1]", E, { E12, E13, E23 } }, + { "stress[0]", 0, S, S11, 1.0 }, + { "stress[0]", 1, S, S12, 1.0 }, + { "stress[0]", 2, S, S13, -1.0 }, + { "stress[1]", 1, S, S22, 1.0 }, + { "stress[1]", 2, S, S23, -1.0 }, + { "stress[2]", 2, S, S33, 1.0 }, + { "strain[0]", 0, E, E11, 1.0 }, + { "strain[0]", 1, E, E12, 1.0 }, + { "strain[0]", 2, E, E13, -1.0 }, + { "strain[1]", 1, E, E22, 1.0 }, + { "strain[1]", 2, E, E23, -1.0 }, + { "strain[2]", 2, E, E33, 1.0 }, }; for ( const auto& mapping : mappings ) { if ( propertyName != mapping.vtkName ) continue; - for ( int compIdx = 0; compIdx < 3; compIdx++ ) + std::string internalKey = mapping.fieldName + separator + mapping.componentName; + + for ( int partIdx = 0; partIdx < femParts->partCount(); partIdx++ ) { - std::string internalKey = mapping.fieldName + separator + mapping.componentNames[compIdx]; - - if ( map->count( internalKey ) == 0 ) ( *map )[internalKey] = {}; - if ( ( *map )[internalKey].count( stepId ) == 0 ) ( *map )[internalKey][stepId] = {}; - - for ( int partIdx = 0; partIdx < femParts->partCount(); partIdx++ ) + size_t dataSize = femParts->part( partIdx )->elementCount(); + std::vector compValues( dataSize, 0.0 ); + for ( size_t i = 0; i < dataSize; i++ ) { - size_t dataSize = femParts->part( partIdx )->elementCount(); - std::vector compValues( dataSize, 0.0 ); - for ( size_t i = 0; i < dataSize; i++ ) - { - compValues[i] = values[i * 3 + compIdx]; - } - ( *map )[internalKey][stepId][partIdx] = compValues; + compValues[i] = mapping.sign * values[i * 3 + mapping.column]; } + ( *map )[internalKey][stepId][partIdx] = compValues; } } } diff --git a/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp b/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp index d9ab31d096..3eb41672aa 100644 --- a/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp +++ b/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp @@ -165,3 +165,36 @@ TEST( RifVtkReaderTest, ReadNodeFieldDisplacement ) EXPECT_NEAR( u3[i], -0.3f, 1e-5f ); } } + +//-------------------------------------------------------------------------------------------------- +/// stress[0..2] are the rows of a symmetric tensor; xz and yz change sign with the negated z +//-------------------------------------------------------------------------------------------------- +TEST( RifVtkReaderTest, ReadElementFieldStressTensorRows ) +{ + 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( "S_VTK" ) > 0 ); + const std::vector expectedComponents = { "S11", "S22", "S33", "S12", "S13", "S23" }; + EXPECT_EQ( elementFields["S_VTK"], expectedComponents ); + + std::vector> values( 6 ); + std::vector*> resultValues; + for ( auto& v : values ) + resultValues.push_back( &v ); + + reader.readElementField( "S_VTK", 0, 1, 0, &resultValues ); + + const std::vector expectedCell0 = { 11, 22, 33, 12, -13, -23 }; + const std::vector expectedCell1 = { 111, 122, 133, 112, -113, -123 }; + 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] ); + } +} diff --git a/ApplicationLibCode/UnitTests/TestData/RifVtkReader/model-00001.vtu b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/model-00001.vtu index 39228a5cd1..0855e9c731 100644 --- a/ApplicationLibCode/UnitTests/TestData/RifVtkReader/model-00001.vtu +++ b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/model-00001.vtu @@ -14,6 +14,15 @@ 150.0 250.0 + + 11 12 13 111 112 113 + + + 12 22 23 112 122 123 + + + 13 23 33 113 123 133 +