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 @@ + + + + + +