diff --git a/ApplicationLibCode/FileInterface/RifVtkImportUtil.cpp b/ApplicationLibCode/FileInterface/RifVtkImportUtil.cpp index 266fcc768b..0afc6070b2 100644 --- a/ApplicationLibCode/FileInterface/RifVtkImportUtil.cpp +++ b/ApplicationLibCode/FileInterface/RifVtkImportUtil.cpp @@ -71,8 +71,8 @@ std::vector RifVtkImportUtil::readDisplacements( const pugi::xml_nod auto pointData = piece.child( "PointData" ); if ( !pointData ) return {}; - auto dataArray = pointData.child( "DataArray" ); - if ( !dataArray || std::string( dataArray.attribute( "Name" ).value() ) != "disp" ) return {}; + auto dataArray = pointData.find_child_by_attribute( "DataArray", "Name", "disp" ); + if ( !dataArray ) return {}; std::string_view text = dataArray.text().get(); std::vector displacements = parseVec3fs( text ); @@ -119,12 +119,27 @@ std::vector RifVtkImportUtil::readConnectivity( const pugi::xml_node& //-------------------------------------------------------------------------------------------------- std::map> RifVtkImportUtil::readProperties( const pugi::xml_node& piece ) { - auto cellData = piece.child( "CellData" ); - if ( !cellData ) return {}; + return readDataArrays( piece.child( "CellData" ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::map> RifVtkImportUtil::readPointProperties( const pugi::xml_node& piece ) +{ + return readDataArrays( piece.child( "PointData" ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::map> RifVtkImportUtil::readDataArrays( const pugi::xml_node& data ) +{ + if ( !data ) return {}; std::map> properties; - for ( const auto& dataArray : cellData.children( "DataArray" ) ) + for ( const auto& dataArray : data.children( "DataArray" ) ) { const char* name = dataArray.attribute( "Name" ).value(); if ( name && *name ) diff --git a/ApplicationLibCode/FileInterface/RifVtkImportUtil.h b/ApplicationLibCode/FileInterface/RifVtkImportUtil.h index 10ebac0dbb..46204729a0 100644 --- a/ApplicationLibCode/FileInterface/RifVtkImportUtil.h +++ b/ApplicationLibCode/FileInterface/RifVtkImportUtil.h @@ -45,8 +45,10 @@ public: static std::vector readDisplacements( const pugi::xml_node& piece ); static std::vector readConnectivity( const pugi::xml_node& piece ); static std::map> readProperties( const pugi::xml_node& piece ); + static std::map> readPointProperties( const pugi::xml_node& piece ); private: - static std::vector parseVec3ds( std::string_view text ); - static std::vector parseVec3fs( std::string_view text ); + static std::map> readDataArrays( const pugi::xml_node& data ); + static std::vector parseVec3ds( std::string_view text ); + static std::vector parseVec3fs( std::string_view text ); }; diff --git a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp index caa37d96a6..e6565d54bf 100644 --- a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp +++ b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp @@ -55,6 +55,11 @@ constexpr const char* E33 = "E33"; constexpr const char* E12 = "E12"; constexpr const char* E13 = "E13"; constexpr const char* E23 = "E23"; + +// Fields in the convention of the derived results +constexpr const char* derivedS = "S"; +constexpr const char* derivedE = "E"; +constexpr const char* por = "POR"; } // namespace RifVtkFieldNames //-------------------------------------------------------------------------------------------------- @@ -83,7 +88,8 @@ void RifVtkReader::close() //-------------------------------------------------------------------------------------------------- bool RifVtkReader::populateDerivedResultNames() const { - return false; + using namespace RifVtkFieldNames; + return m_propertyPartDataElements.contains( std::string( S ) + separator + S11 ) && m_propertyPartDataNodes.contains( por ); } //-------------------------------------------------------------------------------------------------- @@ -129,8 +135,9 @@ bool RifVtkReader::readFemParts( RigFemPartCollection* femParts ) std::map>>> elements; std::map>>> elementSets; std::map>>> properties; + std::map>>> pointProperties; - auto result = read( m_inputPath, parts, nodes, elements, elementSets, properties, m_displacements, m_stepNames ); + auto result = read( m_inputPath, parts, nodes, elements, elementSets, properties, pointProperties, m_displacements, m_stepNames ); if ( !result ) { RiaLogging::error( result.error() ); @@ -223,12 +230,20 @@ bool RifVtkReader::readFemParts( RigFemPartCollection* femParts ) femPart->setElementPartId( femParts->partCount() ); femParts->addFemPart( femPart ); + std::vector& elementNodes = m_elementNodes[femPart->elementPartId()]; + for ( int elmIdx = 0; elmIdx < elmCount; ++elmIdx ) + { + const int* conn = femPart->connectivities( elmIdx ); + elementNodes.insert( elementNodes.end(), conn, conn + RigFemTypes::elementNodeCount( femPart->elementType( elmIdx ) ) ); + } + modelProgress.incrementProgress(); } size_t elementCount = femParts->part( 0 )->elementCount(); readScalarData( femParts, parts, properties, elementCount ); + readNodalData( femParts, pointProperties ); return true; } @@ -243,6 +258,7 @@ std::expected std::map>>>& elements, std::map>>>& elementSets, std::map>>>& properties, + std::map>>>& pointProperties, std::map>>& displacements, std::vector& stepNames ) { @@ -338,6 +354,7 @@ std::expected .toStdString() ); properties[partId].push_back( partProperties ); + pointProperties[partId].push_back( RifVtkImportUtil::readPointProperties( piece ) ); displacements[partId].push_back( partDisplacements ); isFirst = false; @@ -455,7 +472,18 @@ std::map> RifVtkReader::scalarElementField //-------------------------------------------------------------------------------------------------- std::map> RifVtkReader::scalarElementNodeFieldAndComponentNames() { - return {}; + using namespace RifVtkFieldNames; + + std::map> retVal; + if ( m_propertyPartDataElements.contains( std::string( S ) + separator + S11 ) && m_propertyPartDataNodes.contains( por ) ) + { + retVal[derivedS] = { S11, S22, S33, S12, S13, S23 }; + } + if ( m_propertyPartDataElements.contains( std::string( E ) + separator + E11 ) ) + { + retVal[derivedE] = { E11, E22, E33, E12, E13, E23 }; + } + return retVal; } //-------------------------------------------------------------------------------------------------- @@ -463,7 +491,7 @@ std::map> RifVtkReader::scalarElementNodeF //-------------------------------------------------------------------------------------------------- std::map> RifVtkReader::scalarIntegrationPointFieldAndComponentNames() { - return {}; + return scalarElementNodeFieldAndComponentNames(); } //-------------------------------------------------------------------------------------------------- @@ -585,7 +613,43 @@ void RifVtkReader::readElementNodeField( const std::string& field int frameIndex, std::vector*>* resultValues ) { - readField( RigFemResultPosEnum::RIG_ELEMENT_NODAL, fieldName, partIndex, stepIndex, resultValues ); + CAF_ASSERT( resultValues ); + + // 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. + 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 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; + + 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 bool isNormal = c < 3; + values->resize( elementNodes.size() ); + for ( size_t i = 0; i < elementNodes.size(); i++ ) + { + 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 ); + } + } } //-------------------------------------------------------------------------------------------------- @@ -597,7 +661,7 @@ void RifVtkReader::readIntegrationPointField( const std::string& int frameIndex, std::vector*>* resultValues ) { - readField( RigFemResultPosEnum::RIG_INTEGRATION_POINT, fieldName, partIndex, stepIndex, resultValues ); + readElementNodeField( fieldName, partIndex, stepIndex, frameIndex, resultValues ); } //-------------------------------------------------------------------------------------------------- @@ -716,6 +780,30 @@ void RifVtkReader::readScalarData( RigFemPartCollection* } } +//-------------------------------------------------------------------------------------------------- +/// Scalar node data, e.g. from OPM Flow with --vtk-write-geo-mech-nodal. "pressure" becomes the pore pressure POR. +//-------------------------------------------------------------------------------------------------- +void RifVtkReader::readNodalData( RigFemPartCollection* femParts, + const std::map>>>& pointProperties ) +{ + for ( const auto& [partId, propertiesPerTimeStep] : pointProperties ) + { + if ( partId >= femParts->partCount() ) continue; + const size_t nodeCount = femParts->part( partId )->nodeCount(); + + for ( int stepId = 0; stepId < static_cast( propertiesPerTimeStep.size() ); stepId++ ) + { + for ( const auto& [propertyName, values] : propertiesPerTimeStep[stepId] ) + { + if ( values.size() != nodeCount ) continue; + + const std::string fieldName = ( propertyName == "pressure" ) ? RifVtkFieldNames::por : propertyName; + m_propertyPartDataNodes[fieldName][stepId][partId].assign( values.begin(), values.end() ); + } + } + } +} + //-------------------------------------------------------------------------------------------------- /// Map keys: result name / time step / part id //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.h b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.h index 94e45f61f1..43e9bb779b 100644 --- a/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.h +++ b/ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.h @@ -95,16 +95,21 @@ private: const std::map>>>& properties, size_t numElements ); + void readNodalData( RigFemPartCollection* femParts, + const std::map>>>& pointProperties ); + std::map>>>* propertyDataMap( RigFemResultPosEnum resultType ); - static std::expected read( const std::filesystem::path& filePath, - std::map& parts, - std::map>>& nodes, - std::map>>>& elements, - std::map>>>& elementSets, - std::map>>>& properties, - std::map>>& displacements, - std::vector& stepNames ); + static std::expected + read( const std::filesystem::path& filePath, + std::map& parts, + std::map>>& nodes, + std::map>>>& elements, + std::map>>>& elementSets, + std::map>>>& properties, + std::map>>>& pointProperties, + std::map>>& displacements, + std::vector& stepNames ); private: std::map> m_partElementSetNames; @@ -113,4 +118,5 @@ private: std::map>>> m_propertyPartDataNodes; std::map>>> m_propertyPartDataElements; std::map>> m_displacements; + std::map> m_elementNodes; // Node index per element corner, by part }; diff --git a/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp b/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp index 3eb41672aa..f9e2704ad7 100644 --- a/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp +++ b/ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp @@ -198,3 +198,48 @@ TEST( RifVtkReaderTest, ReadElementFieldStressTensorRows ) EXPECT_FLOAT_EQ( values[c][1], expectedCell1[c] ); } } + +//-------------------------------------------------------------------------------------------------- +/// Point data "pressure" is POR; S and E are the cell values at every element corner, S = POR * I - stress +//-------------------------------------------------------------------------------------------------- +TEST( RifVtkReaderTest, ReadElementNodalStressFromCellsAndPor ) +{ + RifVtkReader reader; + reader.openFile( ( testDataDir() / "model.pvd" ).string(), nullptr ); + + cvf::ref femParts = new RigFemPartCollection(); + reader.readFemParts( femParts.p() ); + + EXPECT_TRUE( reader.populateDerivedResultNames() ); + EXPECT_TRUE( reader.scalarNodeFieldAndComponentNames().count( "POR" ) > 0 ); + EXPECT_TRUE( reader.scalarElementNodeFieldAndComponentNames().count( "S" ) > 0 ); + EXPECT_TRUE( reader.scalarElementNodeFieldAndComponentNames().count( "E" ) > 0 ); + + std::vector por; + std::vector*> porValues = { &por }; + reader.readNodeField( "POR", 0, 1, 0, &porValues ); + ASSERT_EQ( por.size(), 12u ); + EXPECT_FLOAT_EQ( por[5], 50.0f ); + + std::vector> s( 6 ), e( 6 ); + std::vector*> sValues, eValues; + for ( int c = 0; c < 6; c++ ) + { + sValues.push_back( &s[c] ); + eValues.push_back( &e[c] ); + } + reader.readElementNodeField( "S", 0, 1, 0, &sValues ); + reader.readElementNodeField( "E", 0, 1, 0, &eValues ); + + // Element 1, local node 2 is node 5 (POR 50). Cell stress rows (111 112 113), (112 122 123), (113 123 133); + // z is negated, so xz and yz change sign. + const size_t idx = 8 + 2; + const std::vector expectedS = { 50 - 111, 50 - 122, 50 - 133, -112, 113, 123 }; + const std::vector expectedE = { 1.1f, 1.4f, 1.6f, 1.2f, -1.3f, -1.5f }; + for ( int c = 0; c < 6; c++ ) + { + ASSERT_EQ( s[c].size(), 16u ); + EXPECT_FLOAT_EQ( s[c][idx], expectedS[c] ); + EXPECT_FLOAT_EQ( e[c][idx], expectedE[c] ); + } +} diff --git a/ApplicationLibCode/UnitTests/TestData/RifVtkReader/model-00001.vtu b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/model-00001.vtu index 0855e9c731..5f7acfca7c 100644 --- a/ApplicationLibCode/UnitTests/TestData/RifVtkReader/model-00001.vtu +++ b/ApplicationLibCode/UnitTests/TestData/RifVtkReader/model-00001.vtu @@ -3,6 +3,9 @@ + + 0 10 20 30 40 50 60 70 80 90 100 110 + 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 @@ -23,6 +26,15 @@ 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 +