From 917e8552df6dd74fc5e20e56c8bd3ce3c895b9ea Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Fri, 23 Oct 2026 21:03:24 +0200 Subject: [PATCH] VTK GeoMech: read element nodal results without changing reader state Look up POR and the tensor components with find(), and POR only for S. Reading E on a model without node pressure no longer made S and the derived results available. Use per-element node offsets instead of assuming eight nodes per element, clear the element nodes when the parts are read again, require the complete stress tensor for derived results, and document the Biot coefficient of 1. --- .../GeoMechFileInterface/RifVtkReader.cpp | 66 ++++++++++++---- .../GeoMechFileInterface/RifVtkReader.h | 9 ++- .../UnitTests/RifVtkReader-Test.cpp | 77 +++++++++++++++++++ .../TestData/RifVtkReader/nopor-00000.vtu | 58 ++++++++++++++ .../UnitTests/TestData/RifVtkReader/nopor.pvd | 6 ++ 5 files changed, 198 insertions(+), 18 deletions(-) create mode 100644 ApplicationLibCode/UnitTests/TestData/RifVtkReader/nopor-00000.vtu create mode 100644 ApplicationLibCode/UnitTests/TestData/RifVtkReader/nopor.pvd diff --git a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp index 56395e8042..0d627ffed0 100644 --- a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp +++ b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp @@ -81,6 +81,23 @@ std::vector availableComponents( const std::map* findValues( const std::map>>>& data, + const std::string& fieldName, + int stepIndex, + int partIndex ) +{ + auto fieldIt = data.find( fieldName ); + if ( fieldIt == data.end() ) return nullptr; + auto stepIt = fieldIt->second.find( stepIndex ); + if ( stepIt == fieldIt->second.end() ) return nullptr; + auto partIt = stepIt->second.find( partIndex ); + if ( partIt == stepIt->second.end() || partIt->second.empty() ) return nullptr; + return &partIt->second; +} } // namespace //-------------------------------------------------------------------------------------------------- @@ -110,7 +127,8 @@ void RifVtkReader::close() bool RifVtkReader::populateDerivedResultNames() const { using namespace RifVtkFieldNames; - return m_propertyPartDataElements.contains( std::string( S ) + separator + S11 ) && m_propertyPartDataNodes.contains( por ); + return availableComponents( m_propertyPartDataElements, S, stressComponents ).size() == stressComponents.size() && + m_propertyPartDataNodes.contains( por ); } //-------------------------------------------------------------------------------------------------- @@ -165,6 +183,8 @@ bool RifVtkReader::readFemParts( RigFemPartCollection* femParts ) return false; } + m_elementNodes.clear(); + auto elementType = result.value(); if ( m_stepNames.empty() ) m_stepNames.push_back( "Geostatic" ); @@ -251,11 +271,14 @@ bool RifVtkReader::readFemParts( RigFemPartCollection* femParts ) femPart->setElementPartId( femParts->partCount() ); femParts->addFemPart( femPart ); - std::vector& elementNodes = m_elementNodes[femPart->elementPartId()]; + ElementNodes& elementNodes = m_elementNodes[femPart->elementPartId()]; + elementNodes.offsets.reserve( elmCount + 1 ); + elementNodes.offsets.push_back( 0 ); for ( int elmIdx = 0; elmIdx < elmCount; ++elmIdx ) { const int* conn = femPart->connectivities( elmIdx ); - elementNodes.insert( elementNodes.end(), conn, conn + RigFemTypes::elementNodeCount( femPart->elementType( elmIdx ) ) ); + elementNodes.nodes.insert( elementNodes.nodes.end(), conn, conn + RigFemTypes::elementNodeCount( femPart->elementType( elmIdx ) ) ); + elementNodes.offsets.push_back( elementNodes.nodes.size() ); } modelProgress.incrementProgress(); @@ -639,37 +662,46 @@ void RifVtkReader::readElementNodeField( const std::string& field // Stress and strain are cell values; every corner of an element gets its cell's value. S is effective // stress, tension positive, from the simulator's total stress (compression positive) and the node pore - // pressure: S = POR * I - S_VTK. Then ST = -S + POR gives the simulator's stress back. + // pressure: S = POR * I - S_VTK. Then ST = -S + POR gives the simulator's stress back. This assumes a Biot + // coefficient of 1. A step without POR or tensor data gives no values. using namespace RifVtkFieldNames; const auto components = scalarElementNodeFieldAndComponentNames()[fieldName]; const auto nodesIt = m_elementNodes.find( partIndex ); if ( components.empty() || nodesIt == m_elementNodes.end() ) return; - const std::vector& elementNodes = nodesIt->second; + const ElementNodes& elementNodes = nodesIt->second; + const size_t elementCount = elementNodes.offsets.size() - 1; const bool isStress = fieldName == derivedS; const std::string cellField = isStress ? S : E; - const std::vector& porValues = m_propertyPartDataNodes[por][stepIndex][partIndex]; - if ( isStress && porValues.empty() ) return; + const std::vector* porValues = nullptr; + if ( isStress ) + { + porValues = findValues( m_propertyPartDataNodes, por, stepIndex, partIndex ); + if ( porValues == nullptr ) return; + } - constexpr int nodesPerElement = 8; for ( size_t c = 0; c < components.size() && c < resultValues->size(); c++ ) { std::vector* values = ( *resultValues )[c]; if ( values == nullptr ) continue; - const std::vector& cellValues = m_propertyPartDataElements[cellField + separator + components[c]][stepIndex][partIndex]; - if ( cellValues.empty() ) continue; + const std::vector* cellValues = + findValues( m_propertyPartDataElements, cellField + separator + components[c], stepIndex, partIndex ); + if ( cellValues == nullptr || cellValues->size() < elementCount ) continue; const bool isNormal = c < 3; - values->resize( elementNodes.size() ); - for ( size_t i = 0; i < elementNodes.size(); i++ ) + values->resize( elementNodes.nodes.size() ); + for ( size_t elmIdx = 0; elmIdx < elementCount; elmIdx++ ) { - const double cellValue = cellValues[i / nodesPerElement]; - if ( isStress ) - ( *values )[i] = static_cast( ( isNormal ? porValues[elementNodes[i]] : 0.0 ) - cellValue ); - else - ( *values )[i] = static_cast( cellValue ); + const double cellValue = ( *cellValues )[elmIdx]; + for ( size_t i = elementNodes.offsets[elmIdx]; i < elementNodes.offsets[elmIdx + 1]; i++ ) + { + if ( isStress ) + ( *values )[i] = static_cast( ( isNormal ? ( *porValues )[elementNodes.nodes[i]] : 0.0 ) - cellValue ); + else + ( *values )[i] = static_cast( cellValue ); + } } } } diff --git a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.h b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.h index 43e9bb779b..5dab4f1221 100644 --- a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.h +++ b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.h @@ -118,5 +118,12 @@ private: std::map>>> m_propertyPartDataNodes; std::map>>> m_propertyPartDataElements; std::map>> m_displacements; - std::map> m_elementNodes; // Node index per element corner, by part + + // Node index per element corner, by part. Element i has the corners [offsets[i], offsets[i + 1]). + struct ElementNodes + { + std::vector nodes; + std::vector offsets; + }; + std::map m_elementNodes; }; diff --git a/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp b/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp index a09bbf8e28..c3d9335205 100644 --- a/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp +++ b/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp @@ -312,3 +312,80 @@ TEST( RifVtkReaderTest, ReadElementFieldPartialTensorRows ) EXPECT_FLOAT_EQ( e[c][1], expectedECell1[c] ); } } + +//-------------------------------------------------------------------------------------------------- +/// Complete stress and strain but no node "pressure": E is available, S and the derived results are not, +/// also after reading E +//-------------------------------------------------------------------------------------------------- +TEST( RifVtkReaderTest, ReadElementNodalStrainWithoutPor ) +{ + RifVtkReader reader; + reader.openFile( ( testDataDir() / "nopor.pvd" ).string(), nullptr ); + + cvf::ref femParts = new RigFemPartCollection(); + reader.readFemParts( femParts.p() ); + + EXPECT_FALSE( reader.populateDerivedResultNames() ); + EXPECT_EQ( reader.scalarElementNodeFieldAndComponentNames().count( "S" ), 0u ); + ASSERT_EQ( reader.scalarElementNodeFieldAndComponentNames().count( "E" ), 1u ); + + std::vector> e( 6 ); + std::vector*> eValues; + for ( auto& v : e ) + eValues.push_back( &v ); + reader.readElementNodeField( "E", 0, 0, 0, &eValues ); + ASSERT_EQ( e[0].size(), 16u ); + EXPECT_FLOAT_EQ( e[0][8], 1.1f ); + + EXPECT_FALSE( reader.populateDerivedResultNames() ); + EXPECT_EQ( reader.scalarElementNodeFieldAndComponentNames().count( "S" ), 0u ); + EXPECT_EQ( reader.scalarNodeFieldAndComponentNames().count( "POR" ), 0u ); +} + +//-------------------------------------------------------------------------------------------------- +/// Step 0 has neither tensor nor POR data: S and E give no values, and reading them adds no fields +//-------------------------------------------------------------------------------------------------- +TEST( RifVtkReaderTest, ReadElementNodalStressAtStepWithoutData ) +{ + RifVtkReader reader; + reader.openFile( ( testDataDir() / "model.pvd" ).string(), nullptr ); + + cvf::ref femParts = new RigFemPartCollection(); + reader.readFemParts( femParts.p() ); + + const auto elementFieldsBefore = reader.scalarElementFieldAndComponentNames(); + const auto nodeFieldsBefore = reader.scalarNodeFieldAndComponentNames(); + + for ( const std::string fieldName : { "S", "E" } ) + { + std::vector> values( 6 ); + std::vector*> resultValues; + for ( auto& v : values ) + resultValues.push_back( &v ); + reader.readElementNodeField( fieldName, 0, 0, 0, &resultValues ); + for ( const auto& v : values ) + EXPECT_TRUE( v.empty() ); + } + + EXPECT_EQ( reader.scalarElementFieldAndComponentNames(), elementFieldsBefore ); + EXPECT_EQ( reader.scalarNodeFieldAndComponentNames(), nodeFieldsBefore ); +} + +//-------------------------------------------------------------------------------------------------- +/// Reading the parts twice gives the same element node results +//-------------------------------------------------------------------------------------------------- +TEST( RifVtkReaderTest, ReadFemPartsTwice ) +{ + RifVtkReader reader; + reader.openFile( ( testDataDir() / "model.pvd" ).string(), nullptr ); + + cvf::ref femParts = new RigFemPartCollection(); + reader.readFemParts( femParts.p() ); + cvf::ref femPartsAgain = new RigFemPartCollection(); + reader.readFemParts( femPartsAgain.p() ); + + std::vector e11; + std::vector*> eValues = { &e11 }; + reader.readElementNodeField( "E", 0, 1, 0, &eValues ); + EXPECT_EQ( e11.size(), 16u ); +} diff --git a/ApplicationLibCode/UnitTests/TestData/RifVtkReader/nopor-00000.vtu b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/nopor-00000.vtu new file mode 100644 index 0000000000..f4903a4283 --- /dev/null +++ b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/nopor-00000.vtu @@ -0,0 +1,58 @@ + + + + + + + 0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 + 0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 + 0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 + 0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 + + + + + 150.0 250.0 + + + 11 12 13 111 112 113 + + + 12 22 23 112 122 123 + + + 13 23 33 113 123 133 + + + 0.1 0.2 0.3 1.1 1.2 1.3 + + + 0.2 0.4 0.5 1.2 1.4 1.5 + + + 0.3 0.5 0.6 1.3 1.5 1.6 + + + + + 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/nopor.pvd b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/nopor.pvd new file mode 100644 index 0000000000..2ee484c7b9 --- /dev/null +++ b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/nopor.pvd @@ -0,0 +1,6 @@ + + + + + +