diff --git a/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h b/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h index 452cefe716..5a7c9e940d 100644 --- a/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h +++ b/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h @@ -127,6 +127,7 @@ public: const RigFemPartCollection* parts() const; int partCount() const; int timeStepCount() const; + std::vector filteredTimeStepNames() const; int frameCount( int timeStepIndex ) const; int totalSteps(); @@ -198,7 +199,6 @@ public: private: RigFemScalarResultFrames* calculateDerivedResult( int partIndex, const RigFemResultAddress& resVarAddr ); - std::vector filteredTimeStepNames() const; private: cvf::Collection m_femPartResults; diff --git a/ApplicationLibCode/ProjectDataModelCommands/RimcGeoMechCase.cpp b/ApplicationLibCode/ProjectDataModelCommands/RimcGeoMechCase.cpp index 979836768c..ed09da76bf 100644 --- a/ApplicationLibCode/ProjectDataModelCommands/RimcGeoMechCase.cpp +++ b/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 + +namespace +{ +struct RegisterScriptEnumNames +{ + RegisterScriptEnumNames() { caf::PdmScriptEnumNameRegistry::registerName( "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 RimGeoMechCase_timeStepStrings::execute() +{ + auto* geoMechCase = self(); + 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 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( RIG_ELEMENT_NODAL ), + "Result Position", + "", + "", + "Result position, for example NODAL or ELEMENT_NODAL" ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::expected RimGeoMechCase_resultNames::execute() +{ + auto* geoMechCase = self(); + 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 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( 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 RimGeoMechCase_readResultValuesInternal::execute() +{ + auto* geoMechCase = self(); + 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 values; + results->globalResultValues( address, m_timeStep(), frameIndex, values ); + + RiaApplication::instance()->keyValueStore()->set( m_key().toStdString(), RiaKeyValueStoreUtil::convertToByteVector( values ) ); + + return nullptr; +} diff --git a/ApplicationLibCode/ProjectDataModelCommands/RimcGeoMechCase.h b/ApplicationLibCode/ProjectDataModelCommands/RimcGeoMechCase.h index 54c6f034ac..103bdd7ed0 100644 --- a/ApplicationLibCode/ProjectDataModelCommands/RimcGeoMechCase.h +++ b/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 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 execute() override; + QString classKeywordReturnedType() const override; + +private: + caf::PdmField> 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 execute() override; + +private: + caf::PdmField> m_resultPosition; + caf::PdmField m_fieldName; + caf::PdmField m_componentName; + caf::PdmField m_timeStep; + caf::PdmField m_frameIndex; + caf::PdmField m_key; +}; diff --git a/GrpcInterface/Python/rips/PythonExamples/export_and_plotting/alter_wbs_plot.py b/GrpcInterface/Python/rips/PythonExamples/geomech/alter_wbs_plot.py similarity index 100% rename from GrpcInterface/Python/rips/PythonExamples/export_and_plotting/alter_wbs_plot.py rename to GrpcInterface/Python/rips/PythonExamples/geomech/alter_wbs_plot.py diff --git a/GrpcInterface/Python/rips/PythonExamples/export_and_plotting/create_wbs_plot.py b/GrpcInterface/Python/rips/PythonExamples/geomech/create_wbs_plot.py similarity index 100% rename from GrpcInterface/Python/rips/PythonExamples/export_and_plotting/create_wbs_plot.py rename to GrpcInterface/Python/rips/PythonExamples/geomech/create_wbs_plot.py diff --git a/GrpcInterface/Python/rips/PythonExamples/geomech/read_geomech_results.py b/GrpcInterface/Python/rips/PythonExamples/geomech/read_geomech_results.py new file mode 100644 index 0000000000..9afb15252a --- /dev/null +++ b/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 +###################################################################### +import math +import sys + +import rips + +if len(sys.argv) < 2: + print("Usage: read_geomech_results.py ") + 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}" + ) diff --git a/GrpcInterface/Python/rips/PythonExamples/case_and_grid_operations/set_geomech_cell_result.py b/GrpcInterface/Python/rips/PythonExamples/geomech/select_geomech_cell_result.py similarity index 100% rename from GrpcInterface/Python/rips/PythonExamples/case_and_grid_operations/set_geomech_cell_result.py rename to GrpcInterface/Python/rips/PythonExamples/geomech/select_geomech_cell_result.py diff --git a/GrpcInterface/Python/rips/case.py b/GrpcInterface/Python/rips/case.py index cfa9340262..909b0c6b31 100644 --- a/GrpcInterface/Python/rips/case.py +++ b/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 diff --git a/GrpcInterface/Python/rips/tests/test_cases.py b/GrpcInterface/Python/rips/tests/test_cases.py index e367d79996..a77a185ffe 100644 --- a/GrpcInterface/Python/rips/tests/test_cases.py +++ b/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""" + + + + + +{coordinates} + + + + {" ".join(map(str, connectivity))} + {offsets} + {types} + + + +{displacements} + + + + {" ".join(map(str, values))} + + + + +""" + with open(os.path.join(folder, f"step_{step}.vtu"), "w") as f: + f.write(vtu) + + pvd = """ + + + + + + +""" + 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 + ) diff --git a/docs/agents/python-script-writer.md b/docs/agents/python-script-writer.md index 950ba92f60..cdf7fe8436 100644 --- a/docs/agents/python-script-writer.md +++ b/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 |