OilfieldToolsNet/ResInsight
4
0 Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2026-10-24 22:27:55 -05:00
Code Issues Packages Projects Releases Wiki Activity

VTK GeoMech: handle partial stress/strain tensor rows

Browse Source 
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.
This commit is contained in:
Kristian Bendiksen authored and Magne Sjaastad committed 2026-10-24 13:59:38 +02:00
1 parent d7b3e60520
commit 5fa48d6e60
4 files changed
+166 -17
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

ApplicationLibCode/GeoMech/GeoMechFileInterface/RifVtkReader.cpp
+42 -17
View File
@@ -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<std::string> stressComponents = { S11, S22, S33, S12, S13, S23 };
const std::vector<std::string> strainComponents = { E11, E22, E33, E12, E13, E23 };
} // namespace RifVtkFieldNames
namespace
{
//--------------------------------------------------------------------------------------------------
/// The components of a tensor field that have data, stored as "<field>__<component>"
//--------------------------------------------------------------------------------------------------
std::vector<std::string> availableComponents( const std::map<std::string, std::map<int, std::map<int, std::vector<double>>>>& data,
const std::string& fieldName,
const std::vector<std::string>& components )
{
std::vector<std::string> 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<std::string, std::vector<std::string>> 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<std::string, std::vector<std::string>> RifVtkReader::scalarElementNodeF
using namespace RifVtkFieldNames;
std::map<std::string, std::vector<std::string>> 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<std::vector<float>*>* 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<std::string> sComponents = { S11, S22, S33, S12, S13, S23 };
static const std::vector<std::string> 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<double> compValues( dataSize, 0.0 );
for ( size_t i = 0; i < dataSize; i++ )
{
ApplicationLibCode/UnitTests/RifVtkReader-Test.cpp
+69
View File
@@ -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<RigFemPartCollection> femParts = new RigFemPartCollection();
reader.readFemParts( femParts.p() );
auto elementFields = reader.scalarElementFieldAndComponentNames();
ASSERT_TRUE( elementFields.count( "E_VTK" ) > 0 );
const std::vector<std::string> expectedComponents = { "E11", "E22", "E33", "E12", "E13", "E23" };
EXPECT_EQ( elementFields["E_VTK"], expectedComponents );
std::vector<std::vector<float>> values( 6 );
std::vector<std::vector<float>*> resultValues;
for ( auto& v : values )
resultValues.push_back( &v );
reader.readElementField( "E_VTK", 0, 1, 0, &resultValues );
const std::vector<float> expectedCell0 = { 0.1f, 0.4f, 0.6f, 0.2f, -0.3f, -0.5f };
const std::vector<float> 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<RigFemPartCollection> femParts = new RigFemPartCollection();
reader.readFemParts( femParts.p() );
auto elementFields = reader.scalarElementFieldAndComponentNames();
const std::vector<std::string> expectedS = { "S11", "S12", "S13" };
const std::vector<std::string> expectedE = { "E22", "E33", "E23" };
EXPECT_EQ( elementFields["S_VTK"], expectedS );
EXPECT_EQ( elementFields["E_VTK"], expectedE );
EXPECT_TRUE( reader.scalarElementNodeFieldAndComponentNames().empty() );
std::vector<std::vector<float>> s( 3 ), e( 3 );
std::vector<std::vector<float>*> sValues = { &s[0], &s[1], &s[2] };
std::vector<std::vector<float>*> 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<float> expectedSCell1 = { 111, 112, -113 };
const std::vector<float> 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] );
}
}
ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial-00000.vtu
+49
View File
@@ -0,0 +1,49 @@
<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
<UnstructuredGrid>
<Piece NumberOfCells="2" NumberOfPoints="12">
<PointData>
<DataArray type="Float32" Name="disp" NumberOfComponents="3" format="ascii">
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
</DataArray>
</PointData>
<CellData>
<DataArray type="Float32" Name="Stress" NumberOfComponents="1" format="ascii">
100.0 200.0
</DataArray>
<DataArray type="Float32" Name="stress[0]" NumberOfComponents="3" format="ascii">
11 12 13 111 112 113
</DataArray>
<DataArray type="Float32" Name="strain[1]" NumberOfComponents="3" format="ascii">
0.12 0.22 0.23 1.12 1.22 1.23
</DataArray>
<DataArray type="Float32" Name="strain[2]" NumberOfComponents="3" format="ascii">
0.13 0.23 0.33 1.13 1.23 1.33
</DataArray>
</CellData>
<Points>
<DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
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
</DataArray>
</Points>
<Cells>
<DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
0 1 4 3 6 7 10 9
1 2 5 4 7 8 11 10
</DataArray>
<DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
8 16
</DataArray>
<DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
12 12
</DataArray>
</Cells>
</Piece>
</UnstructuredGrid>
</VTKFile>
ApplicationLibCode/UnitTests/TestData/RifVtkReader/partial.pvd
+6
View File
@@ -0,0 +1,6 @@
<?xml version="1.0"?>
<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">
<Collection>
<DataSet timestep="1.0" file="partial-00000.vtu"/>
</Collection>
</VTKFile>