Scripting: read GeoMech time steps, result names and values

ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h

@@ -127,6 +127,7 @@ public:
     const RigFemPartCollection* parts() const;
     int                         partCount() const;
     int                         timeStepCount() const;
+    std::vector<std::string>    filteredTimeStepNames() const;
     int                         frameCount( int timeStepIndex ) const;
 
     int                                     totalSteps();
@@ -198,7 +199,6 @@ public:
 
 private:
     RigFemScalarResultFrames* calculateDerivedResult( int partIndex, const RigFemResultAddress& resVarAddr );
-    std::vector<std::string>  filteredTimeStepNames() const;
 
 private:
     cvf::Collection<RigFemPartResults>  m_femPartResults;

ApplicationLibCode/ProjectDataModelCommands/RimcGeoMechCase.cpp

@@ -20,12 +20,31 @@
 
 #include "WellLogCommands/RicNewWellBoreStabilityPlotFeature.h"
 
+#include "KeyValueStore/RiaKeyValueStoreUtil.h"
+#include "RiaApplication.h"
+
+#include "RigFemPartResultsCollection.h"
+#include "RigGeoMechCaseData.h"
+
 #include "RimGeoMechCase.h"
 #include "RimWbsParameters.h"
 #include "RimWellBoreStabilityPlot.h"
 #include "RimWellPath.h"
+#include "RimcDataContainerString.h"
 
 #include "cafPdmFieldScriptingCapability.h"
+#include "cafPdmScriptEnumNameRegistry.h"
+
+#include <algorithm>
+
+namespace
+{
+struct RegisterScriptEnumNames
+{
+    RegisterScriptEnumNames() { caf::PdmScriptEnumNameRegistry::registerName<RigFemResultPosEnum>( "GeoMechResultPosition" ); }
+};
+const RegisterScriptEnumNames s_registerScriptEnumNames;
+} // namespace
 
 CAF_PDM_OBJECT_METHOD_SOURCE_INIT( RimGeoMechCase, RimGeoMechCase_createWellBoreStabilityPlot, "createWellBoreStabilityPlot" );
 
@@ -97,3 +116,191 @@ QString RimGeoMechCase_createWellBoreStabilityPlot::classKeywordReturnedType() c
 {
     return RimWellBoreStabilityPlot::classKeywordStatic();
 }
+
+CAF_PDM_OBJECT_METHOD_SOURCE_INIT( RimGeoMechCase, RimGeoMechCase_timeStepStrings, "TimeStepStringsInternal" );
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RimGeoMechCase_timeStepStrings::RimGeoMechCase_timeStepStrings( caf::PdmObjectHandle* self )
+    : PdmObjectMethod( self, PdmObjectMethod::NullPointerType::NULL_IS_INVALID, PdmObjectMethod::ResultType::PERSISTENT_FALSE )
+{
+    CAF_PDM_InitObject( "Time Step Strings", "", "", "Time step strings for this case." );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::expected<caf::PdmObjectHandle*, QString> RimGeoMechCase_timeStepStrings::execute()
+{
+    auto* geoMechCase = self<RimGeoMechCase>();
+    if ( !geoMechCase ) return std::unexpected( "No GeoMech case is available." );
+
+    auto* container = new RimcDataContainerString();
+
+    if ( geoMechCase->geoMechData() && geoMechCase->geoMechData()->femPartResults() )
+    {
+        // Filtered time step names match the time step index used by result_values()
+        std::vector<QString> names;
+        for ( const auto& name : geoMechCase->geoMechData()->femPartResults()->filteredTimeStepNames() )
+        {
+            names.push_back( QString::fromStdString( name ) );
+        }
+        container->m_stringValues = names;
+    }
+
+    return container;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+QString RimGeoMechCase_timeStepStrings::classKeywordReturnedType() const
+{
+    return RimcDataContainerString::classKeywordStatic();
+}
+
+CAF_PDM_OBJECT_METHOD_SOURCE_INIT( RimGeoMechCase, RimGeoMechCase_resultNames, "ResultNamesInternal" );
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RimGeoMechCase_resultNames::RimGeoMechCase_resultNames( caf::PdmObjectHandle* self )
+    : PdmObjectMethod( self, PdmObjectMethod::NullPointerType::NULL_IS_INVALID, PdmObjectMethod::ResultType::PERSISTENT_FALSE )
+{
+    CAF_PDM_InitObject( "Result Names", "", "", "Result names for this case." );
+
+    CAF_PDM_InitScriptableField( &m_resultPosition,
+                                 "ResultPosition",
+                                 caf::AppEnum<RigFemResultPosEnum>( RIG_ELEMENT_NODAL ),
+                                 "Result Position",
+                                 "",
+                                 "",
+                                 "Result position, for example NODAL or ELEMENT_NODAL" );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::expected<caf::PdmObjectHandle*, QString> RimGeoMechCase_resultNames::execute()
+{
+    auto* geoMechCase = self<RimGeoMechCase>();
+    if ( !geoMechCase ) return std::unexpected( "No GeoMech case is available." );
+
+    if ( !geoMechCase->geoMechData() || !geoMechCase->geoMechData()->femPartResults() )
+    {
+        return std::unexpected( "The GeoMech case has no result data. Open the case first." );
+    }
+
+    const auto fieldsWithComponents = geoMechCase->geoMechData()->femPartResults()->scalarFieldAndComponentNames( m_resultPosition() );
+
+    std::vector<QString> names;
+    for ( const auto& [field, components] : fieldsWithComponents )
+    {
+        const QString fieldName = QString::fromStdString( field );
+        if ( components.empty() )
+        {
+            names.push_back( fieldName );
+            continue;
+        }
+
+        for ( const auto& component : components )
+        {
+            names.push_back( fieldName + ":" + QString::fromStdString( component ) );
+        }
+    }
+
+    auto* container           = new RimcDataContainerString();
+    container->m_stringValues = std::move( names );
+
+    return container;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+QString RimGeoMechCase_resultNames::classKeywordReturnedType() const
+{
+    return RimcDataContainerString::classKeywordStatic();
+}
+
+CAF_PDM_OBJECT_METHOD_SOURCE_INIT( RimGeoMechCase, RimGeoMechCase_readResultValuesInternal, "ReadResultValuesInternal" );
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RimGeoMechCase_readResultValuesInternal::RimGeoMechCase_readResultValuesInternal( caf::PdmObjectHandle* self )
+    : caf::PdmVoidObjectMethod( self )
+{
+    CAF_PDM_InitObject( "Read Result Values", "", "", "Read GeoMech result values" );
+
+    CAF_PDM_InitScriptableField( &m_resultPosition,
+                                 "ResultPosition",
+                                 caf::AppEnum<RigFemResultPosEnum>( RIG_ELEMENT_NODAL ),
+                                 "Result Position",
+                                 "",
+                                 "",
+                                 "Result position, for example NODAL or ELEMENT_NODAL" );
+    CAF_PDM_InitScriptableField( &m_fieldName, "FieldName", QString(), "Field Name" );
+    CAF_PDM_InitScriptableField( &m_componentName, "ComponentName", QString(), "Component Name" );
+    CAF_PDM_InitScriptableField( &m_timeStep, "TimeStep", 0, "Time Step", "", "", "Zero-based time step index" );
+    CAF_PDM_InitScriptableField( &m_frameIndex, "FrameIndex", -1, "Frame Index", "", "", "Zero-based frame index. Negative selects the last frame" );
+    CAF_PDM_InitScriptableField( &m_key, "Key", QString(), "Key" );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::expected<caf::PdmObjectHandle*, QString> RimGeoMechCase_readResultValuesInternal::execute()
+{
+    auto* geoMechCase = self<RimGeoMechCase>();
+    if ( !geoMechCase ) return std::unexpected( "No GeoMech case is available." );
+
+    if ( m_key().isEmpty() ) return std::unexpected( "No key specified." );
+    if ( m_fieldName().isEmpty() ) return std::unexpected( "No result field name specified." );
+
+    if ( !geoMechCase->geoMechData() || !geoMechCase->geoMechData()->femPartResults() )
+    {
+        return std::unexpected( "The GeoMech case has no result data. Open the case first." );
+    }
+
+    auto* results = geoMechCase->geoMechData()->femPartResults();
+
+    if ( m_timeStep() < 0 || m_timeStep() >= results->timeStepCount() )
+    {
+        return std::unexpected( QString( "Time step %1 is out of range [0, %2)." ).arg( m_timeStep() ).arg( results->timeStepCount() ) );
+    }
+
+    const int frameCount = results->frameCount( m_timeStep() );
+    const int frameIndex = m_frameIndex() < 0 ? frameCount - 1 : m_frameIndex();
+    if ( frameIndex < 0 || frameIndex >= frameCount )
+    {
+        return std::unexpected(
+            QString( "Frame %1 (resolved to %2) is out of range [0, %3)." ).arg( m_frameIndex() ).arg( frameIndex ).arg( frameCount ) );
+    }
+
+    const RigFemResultAddress address( m_resultPosition(), m_fieldName().toStdString(), m_componentName().toStdString() );
+    if ( !address.isValid() ) return std::unexpected( QString( "Invalid result %1 %2." ).arg( m_fieldName(), m_componentName() ) );
+
+    const auto fieldsWithComponents = results->scalarFieldAndComponentNames( m_resultPosition() );
+    const auto fieldIt              = fieldsWithComponents.find( address.fieldName );
+    if ( fieldIt == fieldsWithComponents.end() )
+    {
+        return std::unexpected( QString( "Unknown result field '%1'." ).arg( m_fieldName() ) );
+    }
+    const auto& components  = fieldIt->second;
+    const bool  componentOk = components.empty()
+                                  ? address.componentName.empty()
+                                  : std::find( components.begin(), components.end(), address.componentName ) != components.end();
+    if ( !componentOk )
+    {
+        return std::unexpected( QString( "Unknown component '%1' for result field '%2'." ).arg( m_componentName(), m_fieldName() ) );
+    }
+
+    std::vector<float> values;
+    results->globalResultValues( address, m_timeStep(), frameIndex, values );
+
+    RiaApplication::instance()->keyValueStore()->set( m_key().toStdString(), RiaKeyValueStoreUtil::convertToByteVector( values ) );
+
+    return nullptr;
+}

ApplicationLibCode/ProjectDataModelCommands/RimcGeoMechCase.h

@@ -18,6 +18,9 @@
 
 #pragma once
 
+#include "RigFemResultPosEnum.h"
+
+#include "cafAppEnum.h"
 #include "cafPdmField.h"
 #include "cafPdmObjectHandle.h"
 #include "cafPdmObjectMethod.h"
@@ -54,3 +57,60 @@ private:
 
     const RimWbsParameters* m_parameters = nullptr;
 };
+
+//==================================================================================================
+/// Get the time step strings of this GeoMech case.
+//==================================================================================================
+class RimGeoMechCase_timeStepStrings : public caf::PdmObjectMethod
+{
+    CAF_PDM_HEADER_INIT;
+
+public:
+    RimGeoMechCase_timeStepStrings( caf::PdmObjectHandle* self );
+
+    std::expected<caf::PdmObjectHandle*, QString> execute() override;
+    QString                                       classKeywordReturnedType() const override;
+};
+
+//==================================================================================================
+/// Get the result names of this GeoMech case for a result position.
+///
+/// Names are returned as "FIELD:COMPONENT", or "FIELD" for results without components.
+//==================================================================================================
+class RimGeoMechCase_resultNames : public caf::PdmObjectMethod
+{
+    CAF_PDM_HEADER_INIT;
+
+public:
+    RimGeoMechCase_resultNames( caf::PdmObjectHandle* self );
+
+    std::expected<caf::PdmObjectHandle*, QString> execute() override;
+    QString                                       classKeywordReturnedType() const override;
+
+private:
+    caf::PdmField<caf::AppEnum<RigFemResultPosEnum>> m_resultPosition;
+};
+
+//==================================================================================================
+/// Read the values of a GeoMech result into the key-value store.
+///
+/// A negative frame index selects the last frame of the time step. The values of all parts are
+/// stored as a float array under the given key.
+//==================================================================================================
+class RimGeoMechCase_readResultValuesInternal : public caf::PdmVoidObjectMethod
+{
+    CAF_PDM_HEADER_INIT;
+
+public:
+    RimGeoMechCase_readResultValuesInternal( caf::PdmObjectHandle* self );
+
+    std::expected<caf::PdmObjectHandle*, QString> execute() override;
+
+private:
+    caf::PdmField<caf::AppEnum<RigFemResultPosEnum>> m_resultPosition;
+    caf::PdmField<QString>                           m_fieldName;
+    caf::PdmField<QString>                           m_componentName;
+    caf::PdmField<int>                               m_timeStep;
+    caf::PdmField<int>                               m_frameIndex;
+    caf::PdmField<QString>                           m_key;
+};

GrpcInterface/Python/rips/PythonExamples/geomech/read_geomech_results.py

@@ -0,0 +1,56 @@
+######################################################################
+# This script loads a geomechanical case (.odb, .inp or VTK .pvd) and
+# reads time steps, result names and result values into Python lists.
+#
+# Usage: python read_geomech_results.py <path to geomech case (.odb, .inp or VTK .pvd)>
+######################################################################
+import math
+import sys
+
+import rips
+
+if len(sys.argv) < 2:
+    print("Usage: read_geomech_results.py <path to geomech case>")
+    sys.exit(1)
+
+try:
+    resinsight = rips.Instance.find()
+    case = resinsight.project.load_case(sys.argv[1])
+except Exception as exc:
+    print(f"Failed to load geomechanical case: {exc}")
+    sys.exit(1)
+
+if not isinstance(case, rips.GeoMechCase):
+    print("The file is not a geomechanical case")
+    sys.exit(1)
+
+time_steps = case.time_step_strings()
+print("Time steps:", time_steps)
+
+# Result names are "FIELD:COMPONENT", or "FIELD" for results without components
+for position in rips.GeoMechResultPosition:
+    if position in (
+        rips.GeoMechResultPosition.WELLPATH_DERIVED,
+        rips.GeoMechResultPosition.DIFFERENTIALS,
+    ):
+        continue
+    print(position, case.result_names(position))
+
+# Print statistics for every result with defined values at one position.
+# A negative frame_index (the default) selects the last frame of the time step.
+position = rips.GeoMechResultPosition.NODAL
+for name in case.result_names(position):
+    field_name, _, component_name = name.partition(":")
+    for time_step, label in enumerate(time_steps):
+        values = case.result_values(
+            field_name,
+            component_name,
+            time_step=time_step,
+            result_position=position,
+        )
+        finite = [v for v in values if math.isfinite(v)]
+        if finite:
+            print(
+                f"{name}, {label}: {len(values)} values, "
+                f"min={min(finite):.4g}, max={max(finite):.4g}"
+            )

GrpcInterface/Python/rips/case.py

@@ -54,6 +54,7 @@ from .resinsight_classes import (
     Case as Case,
     EclipseCase as EclipseCase,
     GeoMechCase as GeoMechCase,
+    GeoMechResultPosition as GeoMechResultPosition,
     PorosityModelType as PorosityModelType,
     PropertyDataType as PropertyDataType,
     PropertyType as PropertyType,
@@ -303,6 +304,76 @@ def time_steps(self) -> List[Case_pb2.TimeStepDate]:
     return self.__case_stub.GetTimeSteps(self.__request()).dates
 
 
+@add_method(GeoMechCase)
+def time_step_strings(self: GeoMechCase) -> List[str]:
+    """Get the time step strings of this geomechanical case
+
+    Returns:
+        A list of time step strings, one per time step. The list index
+        matches the time_step argument of result_values().
+    """
+    return self.time_step_strings_internal().values
+
+
+@add_method(GeoMechCase)
+def result_names(
+    self: GeoMechCase,
+    result_position: GeoMechResultPosition = GeoMechResultPosition.ELEMENT_NODAL,
+) -> List[str]:
+    """Get the result names of this geomechanical case
+
+    Arguments:
+        result_position(GeoMechResultPosition): Result position, for example
+            GeoMechResultPosition.NODAL or GeoMechResultPosition.ELEMENT_NODAL.
+
+    Returns:
+        A list of names on the form "FIELD:COMPONENT", or "FIELD" for results
+        without components.
+    """
+    return self.result_names_internal(result_position=result_position).values
+
+
+@add_method(GeoMechCase)
+def result_values(
+    self: GeoMechCase,
+    field_name: str,
+    component_name: str = "",
+    time_step: int = 0,
+    result_position: GeoMechResultPosition = GeoMechResultPosition.ELEMENT_NODAL,
+    frame_index: int = -1,
+) -> List[float]:
+    """Read the values of a geomechanical result
+
+    Arguments:
+        field_name(str): Result field name, e.g. "ST" or "POR-Bar".
+        component_name(str): Result component name, e.g. "S1". Empty for results without components.
+        time_step(int): Zero-based time step index.
+        result_position(GeoMechResultPosition): Result position of the result.
+        frame_index(int): Zero-based frame index. Negative selects the last frame of the time step.
+
+    Returns:
+        The values of all parts, in the order used by ResInsight for the result position.
+        Values are NaN where the result is not defined.
+    """
+    project = self.ancestor(Project)
+    if project is None:
+        raise RipsError("Could not find parent project")
+
+    key = f"{uuid.uuid4()}_geomech_result"
+    try:
+        self.read_result_values_internal(
+            result_position=result_position,
+            field_name=field_name,
+            component_name=component_name,
+            time_step=time_step,
+            frame_index=frame_index,
+            key=key,
+        )
+        return project.key_values(key)  # type: ignore[attr-defined]
+    finally:
+        project.remove_key_values(key)  # type: ignore[attr-defined]
+
+
 @add_method(Case)
 def reservoir_boundingbox(self) -> Case_pb2.BoundingBox:
     """Get the reservoir bounding box

GrpcInterface/Python/rips/tests/test_cases.py

@@ -470,3 +470,102 @@ def test_create_well_bore_stability_plot_is_geomech_only(
     case = rips_instance.project.load_case(case_path)
     assert not hasattr(case, "create_well_bore_stability_plot")
     assert hasattr(rips.GeoMechCase, "create_well_bore_stability_plot")
+
+
+def _write_small_vtk_case(folder):
+    # Two hexahedron cells sharing a face, with two time steps and one cell property
+    points = [(x, y, z) for x in range(3) for y in range(2) for z in range(2)]
+    point_index = {p: i for i, p in enumerate(points)}
+
+    def hex_connectivity(x0):
+        corners = [
+            (x0, 0, 0),
+            (x0 + 1, 0, 0),
+            (x0 + 1, 1, 0),
+            (x0, 1, 0),
+            (x0, 0, 1),
+            (x0 + 1, 0, 1),
+            (x0 + 1, 1, 1),
+            (x0, 1, 1),
+        ]
+        return [point_index[c] for c in corners]
+
+    connectivity = hex_connectivity(0) + hex_connectivity(1)
+    coordinates = "\n".join(f"{x} {y} {z}" for x, y, z in points)
+    displacements = "\n".join("0 0 0" for _ in points)
+    offsets = " ".join(str(8 * (i + 1)) for i in range(2))
+    types = " ".join("12" for _ in range(2))
+
+    property_values = [[0.1, 0.2], [0.3, 0.4]]
+    for step, values in enumerate(property_values):
+        vtu = f"""<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="{len(points)}" NumberOfCells="2">
+      <Points>
+        <DataArray type="Float64" Name="Coordinates" NumberOfComponents="3" format="ascii">
+{coordinates}
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int32" Name="connectivity" format="ascii">{" ".join(map(str, connectivity))}</DataArray>
+        <DataArray type="Int32" Name="offsets" format="ascii">{offsets}</DataArray>
+        <DataArray type="UInt8" Name="types" format="ascii">{types}</DataArray>
+      </Cells>
+      <PointData>
+        <DataArray type="Float32" Name="disp" NumberOfComponents="3" format="ascii">
+{displacements}
+        </DataArray>
+      </PointData>
+      <CellData>
+        <DataArray type="Float32" Name="porosity" format="ascii">{" ".join(map(str, values))}</DataArray>
+      </CellData>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>
+"""
+        with open(os.path.join(folder, f"step_{step}.vtu"), "w") as f:
+            f.write(vtu)
+
+    pvd = """<?xml version="1.0"?>
+<VTKFile type="Collection" version="0.1">
+  <Collection>
+    <DataSet timestep="0" file="step_0.vtu"/>
+    <DataSet timestep="1" file="step_1.vtu"/>
+  </Collection>
+</VTKFile>
+"""
+    pvd_path = os.path.join(folder, "small.pvd")
+    with open(pvd_path, "w") as f:
+        f.write(pvd)
+    return pvd_path, property_values
+
+
+def test_geomech_case_results(rips_instance, initialize_test, tmp_path):
+    pvd_path, property_values = _write_small_vtk_case(str(tmp_path))
+    case = rips_instance.project.load_case(pvd_path)
+    assert isinstance(case, rips.GeoMechCase)
+
+    assert len(case.time_step_strings()) == 2
+
+    positions = rips.GeoMechResultPosition
+    names = case.result_names(positions.ELEMENT)
+    assert "porosity" in names
+
+    for time_step, expected in enumerate(property_values):
+        values = case.result_values(
+            "porosity", time_step=time_step, result_position=positions.ELEMENT
+        )
+        assert values == pytest.approx(expected)
+
+    # The last frame is the default, and is the only frame in this case
+    assert case.result_values(
+        "porosity", time_step=1, result_position=positions.ELEMENT, frame_index=0
+    ) == pytest.approx(property_values[1])
+
+    with pytest.raises(rips.RipsError):
+        case.result_values("porosity", time_step=2, result_position=positions.ELEMENT)
+    with pytest.raises(rips.RipsError):
+        case.result_values(
+            "porosity", time_step=0, result_position=positions.ELEMENT, frame_index=1
+        )

docs/agents/python-script-writer.md

@@ -63,6 +63,7 @@ Browse `GrpcInterface/Python/rips/PythonExamples/` for reference scripts coverin
 | Directory | Topics covered |
 |-----------|---------------|
 | `case_and_grid_operations/` | Grid access, cell properties, polygon regions, bounding box |
+| `geomech/` | Geomechanical cases: time steps, result names and values, cell results |
 | `surfaces_and_visualization/` | Surfaces, polygons, view setup |
 | `well_operations/` | Well paths, well log data |
 | `simulation_and_results/` | Time-step iteration, dynamic properties |