Source code for geos.processing.post_processing.GeomechanicsCalculator

# SPDX-License-Identifier: Apache-2.0
# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
# SPDX-FileContributor: Martin Lemay, Romain Baville
# ruff: noqa: E402 # disable Module level import not at top of file
from typing import Union
from typing_extensions import Self
from dataclasses import dataclass

import logging
import numpy as np
import numpy.typing as npt

import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts

from geos.mesh.utils.arrayModifiers import createAttribute
from geos.mesh.utils.arrayHelpers import (
    getArrayInObject,
    isAttributeInObject,
)

from geos.utils.Logger import (
    Logger,
    getLogger,
)
from geos.utils.GeosOutputsConstants import (
    AttributeEnum,
    ComponentNameEnum,
    GeosMeshOutputsEnum,
    PostProcessingOutputsEnum,
)
from geos.utils.PhysicalConstants import (
    DEFAULT_FRICTION_ANGLE_RAD,
    DEFAULT_GRAIN_BULK_MODULUS,
    DEFAULT_ROCK_COHESION,
    GRAVITY,
    WATER_DENSITY,
)

from vtkmodules.vtkCommonDataModel import vtkPointSet, vtkUnstructuredGrid
from vtkmodules.vtkFiltersCore import vtkCellCenters

__doc__ = """
GeomechanicsCalculator module is a vtk filter that allows to compute additional geomechanical properties from existing ones:
    - The Young modulus and the Poisson's ratio named "youngModulus" and "poissonRatio" or bulk and shear moduli named "bulkModulus" and "shearModulus"
    - The initial Young modulus and Poisson's ratio named "youngModulusInitial" and "poissonRatioInitial" or the initial bulk modulus named "bulkModulusInitial"
    - The porosity named "porosity"
    - The initial porosity named "porosityInitial"
    - The delta of pressure named "deltaPressure"
    - The density named "density"
    - The effective stress named "stressEffective"
    - The initial effective stress named "stressEffectiveInitial"
    - The pressure named "pressure"

The basic geomechanics outputs are:
    - The elastic moduli not present on the mesh
    - Biot coefficient
    - Compressibility, oedometric compressibility and real compressibility coefficient
    - Specific gravity
    - Real effective stress ratio
    - Total initial stress, total current stress and total stress ratio
    - Lithostatic stress (physic to update)
    - Elastic stain
    - Real reservoir stress path and reservoir stress path in oedometric condition

The advanced geomechanics outputs are:
    - Fracture index and threshold
    - Critical pore pressure and pressure index

GeomechanicsCalculator filter input mesh is either vtkPointSet or vtkUnstructuredGrid
and returned mesh is of same type as input.

.. Important::
     Please refer to the GeosExtractMergeBlockVolume* filters to provide the correct input.

.. Note::
    - The default physical constants used by the filter are:
        - grainBulkModulus = 38e9 Pa ( quartz value )
        - specificDensity = 1000.0 kg/m³ ( water value )
        - rockCohesion = 0.0 Pa ( fractured case )
        - frictionAngle = 10.0°

To use the filter:

.. code-block:: python

    import numpy as np
    from geos.mesh.processing.GeomechanicsCalculator import GeomechanicsCalculator

    # Define filter inputs
    mesh: Union[ vtkPointSet, vtkUnstructuredGrid ]
    computeAdvancedOutputs: bool # optional, defaults to False
    speHandler: bool # optional, defaults to False

    # Instantiate the filter
    filter: GeomechanicsCalculator = GeomechanicsCalculator( mesh, computeAdvancedOutputs, speHandler )

    # Use your own handler (if speHandler is True)
    yourHandler: logging.Handler
    filter.setLoggerHandler( yourHandler )

    # Change the physical constants if needed
    ## Basic outputs
    grainBulkModulus: float
    specificDensity: float
    filter.physicalConstants.grainBulkModulus = grainBulkModulus
    filter.physicalConstants.specificDensity = specificDensity

    ## Advanced outputs
    rockCohesion: float
    frictionAngle: float
    filter.physicalConstants.rockCohesion = rockCohesion
    filter.physicalConstants.frictionAngle = frictionAngle

    # Do calculations
    filter.applyFilter()

    # Get the mesh with the geomechanical output as attribute
    output: Union[vtkPointSet, vtkUnstructuredGrid]
    output = filter.getOutput()
"""

loggerTitle: str = "Geomechanical Calculator Filter"

# Elastic Moduli:
BULK_MODULUS: AttributeEnum = GeosMeshOutputsEnum.BULK_MODULUS
SHEAR_MODULUS: AttributeEnum = GeosMeshOutputsEnum.SHEAR_MODULUS
YOUNG_MODULUS: AttributeEnum = PostProcessingOutputsEnum.YOUNG_MODULUS
POISSON_RATIO: AttributeEnum = PostProcessingOutputsEnum.POISSON_RATIO
BULK_MODULUS_T0: AttributeEnum = PostProcessingOutputsEnum.BULK_MODULUS_INITIAL
YOUNG_MODULUS_T0: AttributeEnum = PostProcessingOutputsEnum.YOUNG_MODULUS_INITIAL
POISSON_RATIO_T0: AttributeEnum = PostProcessingOutputsEnum.POISSON_RATIO_INITIAL
ELASTIC_MODULI: tuple[ AttributeEnum, ...] = ( BULK_MODULUS, SHEAR_MODULUS, YOUNG_MODULUS, POISSON_RATIO,
                                               BULK_MODULUS_T0, YOUNG_MODULUS_T0, POISSON_RATIO_T0 )

# Mandatory attributes:
POROSITY: AttributeEnum = GeosMeshOutputsEnum.POROSITY
POROSITY_T0: AttributeEnum = GeosMeshOutputsEnum.POROSITY_INI
PRESSURE: AttributeEnum = GeosMeshOutputsEnum.PRESSURE
DELTA_PRESSURE: AttributeEnum = GeosMeshOutputsEnum.DELTA_PRESSURE
DENSITY: AttributeEnum = GeosMeshOutputsEnum.ROCK_DENSITY
STRESS_EFFECTIVE: AttributeEnum = GeosMeshOutputsEnum.STRESS_EFFECTIVE
STRESS_EFFECTIVE_T0: AttributeEnum = PostProcessingOutputsEnum.STRESS_EFFECTIVE_INITIAL
MANDATORY_ATTRIBUTES: tuple[ AttributeEnum, ...] = ( POROSITY, POROSITY_T0, PRESSURE, DELTA_PRESSURE, DENSITY,
                                                     STRESS_EFFECTIVE, STRESS_EFFECTIVE_T0 )

# Basic outputs:
BIOT_COEFFICIENT: AttributeEnum = PostProcessingOutputsEnum.BIOT_COEFFICIENT
COMPRESSIBILITY: AttributeEnum = PostProcessingOutputsEnum.COMPRESSIBILITY
COMPRESSIBILITY_OED: AttributeEnum = PostProcessingOutputsEnum.COMPRESSIBILITY_OED
COMPRESSIBILITY_REAL: AttributeEnum = PostProcessingOutputsEnum.COMPRESSIBILITY_REAL
SPECIFIC_GRAVITY: AttributeEnum = PostProcessingOutputsEnum.SPECIFIC_GRAVITY
STRESS_EFFECTIVE_RATIO_REAL: AttributeEnum = PostProcessingOutputsEnum.STRESS_EFFECTIVE_RATIO_REAL
STRESS_TOTAL: AttributeEnum = PostProcessingOutputsEnum.STRESS_TOTAL
STRESS_TOTAL_T0: AttributeEnum = PostProcessingOutputsEnum.STRESS_TOTAL_INITIAL
STRESS_TOTAL_RATIO_REAL: AttributeEnum = PostProcessingOutputsEnum.STRESS_TOTAL_RATIO_REAL
LITHOSTATIC_STRESS: AttributeEnum = PostProcessingOutputsEnum.LITHOSTATIC_STRESS
STRAIN_ELASTIC: AttributeEnum = PostProcessingOutputsEnum.STRAIN_ELASTIC
STRESS_TOTAL_DELTA: AttributeEnum = PostProcessingOutputsEnum.STRESS_TOTAL_DELTA
RSP_REAL: AttributeEnum = PostProcessingOutputsEnum.RSP_REAL
RSP_OED: AttributeEnum = PostProcessingOutputsEnum.RSP_OED
STRESS_EFFECTIVE_RATIO_OED: AttributeEnum = PostProcessingOutputsEnum.STRESS_EFFECTIVE_RATIO_OED
BASIC_OUTPUTS: tuple[ AttributeEnum,
                      ...] = ( BIOT_COEFFICIENT, COMPRESSIBILITY, COMPRESSIBILITY_OED, COMPRESSIBILITY_REAL,
                               SPECIFIC_GRAVITY, STRESS_EFFECTIVE_RATIO_REAL, STRESS_TOTAL, STRESS_TOTAL_T0,
                               STRESS_TOTAL_RATIO_REAL, LITHOSTATIC_STRESS, STRAIN_ELASTIC, STRESS_TOTAL_DELTA,
                               RSP_REAL, RSP_OED, STRESS_EFFECTIVE_RATIO_OED )

# Advanced outputs:
CRITICAL_TOTAL_STRESS_RATIO: AttributeEnum = PostProcessingOutputsEnum.CRITICAL_TOTAL_STRESS_RATIO
TOTAL_STRESS_RATIO_THRESHOLD: AttributeEnum = PostProcessingOutputsEnum.TOTAL_STRESS_RATIO_THRESHOLD
CRITICAL_PORE_PRESSURE: AttributeEnum = PostProcessingOutputsEnum.CRITICAL_PORE_PRESSURE
CRITICAL_PORE_PRESSURE_THRESHOLD: AttributeEnum = PostProcessingOutputsEnum.CRITICAL_PORE_PRESSURE_THRESHOLD
ADVANCED_OUTPUTS: tuple[ AttributeEnum, ...] = ( CRITICAL_TOTAL_STRESS_RATIO, TOTAL_STRESS_RATIO_THRESHOLD,
                                                 CRITICAL_PORE_PRESSURE, CRITICAL_PORE_PRESSURE_THRESHOLD )


[docs] class GeomechanicsCalculator:
[docs] @dataclass class PhysicalConstants: """The dataclass with the value of the physical constant used to compute geomechanics properties.""" ## Basic outputs _grainBulkModulus: float = DEFAULT_GRAIN_BULK_MODULUS _specificDensity: float = WATER_DENSITY ## Advanced outputs _rockCohesion: float = DEFAULT_ROCK_COHESION _frictionAngle: float = DEFAULT_FRICTION_ANGLE_RAD @property def grainBulkModulus( self: Self ) -> float: """Get the grain bulk modulus value.""" return self._grainBulkModulus @grainBulkModulus.setter def grainBulkModulus( self: Self, value: float ) -> None: self._grainBulkModulus = value @property def specificDensity( self: Self ) -> float: """Get the specific density value.""" return self._specificDensity @specificDensity.setter def specificDensity( self: Self, value: float ) -> None: self._specificDensity = value @property def rockCohesion( self: Self ) -> float: """Get the rock cohesion value.""" return self._rockCohesion @rockCohesion.setter def rockCohesion( self: Self, value: float ) -> None: self._rockCohesion = value @property def frictionAngle( self: Self ) -> float: """Get the friction angle value.""" return self._frictionAngle @frictionAngle.setter def frictionAngle( self: Self, value: float ) -> None: self._frictionAngle = value
[docs] @dataclass class ElasticModuliValue: """The dataclass with the value of the elastic moduli.""" _bulkModulus: npt.NDArray[ np.float64 ] | None = None _shearModulus: npt.NDArray[ np.float64 ] | None = None _youngModulus: npt.NDArray[ np.float64 ] | None = None _poissonRatio: npt.NDArray[ np.float64 ] | None = None _bulkModulusT0: npt.NDArray[ np.float64 ] | None = None _youngModulusT0: npt.NDArray[ np.float64 ] | None = None _poissonRatioT0: npt.NDArray[ np.float64 ] | None = None @property def bulkModulus( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the bulk modulus value.""" return self._bulkModulus @bulkModulus.setter def bulkModulus( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._bulkModulus = value @property def shearModulus( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the shear modulus value.""" return self._shearModulus @shearModulus.setter def shearModulus( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._shearModulus = value @property def youngModulus( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the young modulus value.""" return self._youngModulus @youngModulus.setter def youngModulus( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._youngModulus = value @property def poissonRatio( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the poisson ratio value.""" return self._poissonRatio @poissonRatio.setter def poissonRatio( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._poissonRatio = value @property def bulkModulusT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the bulk modulus at the initial time value.""" return self._bulkModulusT0 @bulkModulusT0.setter def bulkModulusT0( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._bulkModulusT0 = value @property def youngModulusT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the young modulus at the initial time value.""" return self._youngModulusT0 @youngModulusT0.setter def youngModulusT0( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._youngModulusT0 = value @property def poissonRatioT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the poisson ration at the initial time value.""" return self._poissonRatioT0 @poissonRatioT0.setter def poissonRatioT0( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._poissonRatioT0 = value
[docs] def setElasticModulusValue( self: Self, name: str, value: npt.NDArray[ np.float64 ] ) -> None: """Set the elastic modulus value wanted. Args: name (str): The name of the elastic modulus. value (npt.NDArray[np.float64]): The value to set. """ if name == BULK_MODULUS.attributeName: self.bulkModulus = value elif name == BULK_MODULUS_T0.attributeName: self.bulkModulusT0 = value elif name == SHEAR_MODULUS.attributeName: self.shearModulus = value elif name == YOUNG_MODULUS.attributeName: self.youngModulus = value elif name == YOUNG_MODULUS_T0.attributeName: self.youngModulusT0 = value elif name == POISSON_RATIO.attributeName: self.poissonRatio = value elif name == POISSON_RATIO_T0.attributeName: self.poissonRatioT0 = value
[docs] def getElasticModulusValue( self: Self, name: str ) -> npt.NDArray[ np.float64 ] | None: """Get the wanted elastic modulus value. Args: name (str): The name of the wanted elastic modulus. Returns: npt.NDArray[np.float64]: The value of the elastic modulus. """ if name == BULK_MODULUS.attributeName: return self.bulkModulus elif name == BULK_MODULUS_T0.attributeName: return self.bulkModulusT0 elif name == SHEAR_MODULUS.attributeName: return self.shearModulus elif name == YOUNG_MODULUS.attributeName: return self.youngModulus elif name == YOUNG_MODULUS_T0.attributeName: return self.youngModulusT0 elif name == POISSON_RATIO.attributeName: return self.poissonRatio elif name == POISSON_RATIO_T0.attributeName: return self.poissonRatioT0 else: raise NameError
[docs] @dataclass class MandatoryAttributesValue: """The dataclass with the value of mandatory properties to have to compute other geomechanics properties.""" _porosity: npt.NDArray[ np.float64 ] | None = None _porosityInitial: npt.NDArray[ np.float64 ] | None = None _pressure: npt.NDArray[ np.float64 ] | None = None _deltaPressure: npt.NDArray[ np.float64 ] | None = None _density: npt.NDArray[ np.float64 ] | None = None _effectiveStress: npt.NDArray[ np.float64 ] | None = None _effectiveStressT0: npt.NDArray[ np.float64 ] | None = None @property def porosity( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the porosity value.""" return self._porosity @property def porosityInitial( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the initial porosity value.""" return self._porosityInitial @property def pressure( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the pressure value.""" return self._pressure @property def deltaPressure( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the delta pressure value.""" return self._deltaPressure @property def density( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the density value.""" return self._density @property def effectiveStress( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the effective stress value.""" return self._effectiveStress @property def effectiveStressT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the initial effective stress value.""" return self._effectiveStressT0
[docs] def setMandatoryAttributeValue( self: Self, name: str, value: npt.NDArray[ np.float64 ] ) -> None: """Set the value of a mandatory property. Args: name (str): The name of the property. value (npt.NDArray[np.float64]): The value to set. """ if name == POROSITY.attributeName: self._porosity = value elif name == POROSITY_T0.attributeName: self._porosityInitial = value elif name == PRESSURE.attributeName: self._pressure = value elif name == DELTA_PRESSURE.attributeName: self._deltaPressure = value elif name == DENSITY.attributeName: self._density = value elif name == STRESS_EFFECTIVE.attributeName: self._effectiveStress = value elif name == STRESS_EFFECTIVE_T0.attributeName: self._effectiveStressT0 = value
[docs] @dataclass class BasicOutputValue: """The dataclass with the value of the basic geomechanics outputs.""" _biotCoefficient: npt.NDArray[ np.float64 ] | None = None _compressibility: npt.NDArray[ np.float64 ] | None = None _compressibilityOed: npt.NDArray[ np.float64 ] | None = None _compressibilityReal: npt.NDArray[ np.float64 ] | None = None _specificGravity: npt.NDArray[ np.float64 ] | None = None _effectiveStressRatioReal: npt.NDArray[ np.float64 ] | None = None _totalStress: npt.NDArray[ np.float64 ] | None = None _totalStressT0: npt.NDArray[ np.float64 ] | None = None _totalStressRatioReal: npt.NDArray[ np.float64 ] | None = None # _lithostaticStress: npt.NDArray[ np.float64 ] | None = None _elasticStrain: npt.NDArray[ np.float64 ] | None = None _deltaTotalStress: npt.NDArray[ np.float64 ] | None = None _rspReal: npt.NDArray[ np.float64 ] | None = None _rspOed: npt.NDArray[ np.float64 ] | None = None _effectiveStressRatioOed: npt.NDArray[ np.float64 ] | None = None @property def biotCoefficient( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the biot coefficient value.""" return self._biotCoefficient @biotCoefficient.setter def biotCoefficient( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._biotCoefficient = value @property def compressibility( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the compressibility value.""" return self._compressibility @compressibility.setter def compressibility( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._compressibility = value @property def compressibilityOed( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the compressibility in oedometric condition value.""" return self._compressibilityOed @compressibilityOed.setter def compressibilityOed( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._compressibilityOed = value @property def compressibilityReal( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the real compressibility value.""" return self._compressibilityReal @compressibilityReal.setter def compressibilityReal( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._compressibilityReal = value @property def specificGravity( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the specific gravity value.""" return self._specificGravity @specificGravity.setter def specificGravity( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._specificGravity = value @property def effectiveStressRatioReal( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the real effective stress ratio value.""" return self._effectiveStressRatioReal @effectiveStressRatioReal.setter def effectiveStressRatioReal( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._effectiveStressRatioReal = value @property def totalStress( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the total stress value.""" return self._totalStress @totalStress.setter def totalStress( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._totalStress = value @property def totalStressT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the initial total stress value.""" return self._totalStressT0 @totalStressT0.setter def totalStressT0( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._totalStressT0 = value @property def totalStressRatioReal( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the total real stress ratio value.""" return self._totalStressRatioReal @totalStressRatioReal.setter def totalStressRatioReal( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._totalStressRatioReal = value @property def lithostaticStress( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the lithostatic stress value.""" return self._lithostaticStress @lithostaticStress.setter def lithostaticStress( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._lithostaticStress = value @property def elasticStrain( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the elastic strain value.""" return self._elasticStrain @elasticStrain.setter def elasticStrain( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._elasticStrain = value @property def deltaTotalStress( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the total delta stress value.""" return self._deltaTotalStress @deltaTotalStress.setter def deltaTotalStress( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._deltaTotalStress = value @property def rspReal( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the real reservoir stress path value.""" return self._rspReal @rspReal.setter def rspReal( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._rspReal = value @property def rspOed( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the reservoir stress path in oedometric condition value.""" return self._rspOed @rspOed.setter def rspOed( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._rspOed = value @property def effectiveStressRatioOed( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the effective stress ratio oedometric value.""" return self._effectiveStressRatioOed @effectiveStressRatioOed.setter def effectiveStressRatioOed( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._effectiveStressRatioOed = value
[docs] def getBasicOutputValue( self: Self, name: str ) -> npt.NDArray[ np.float64 ] | None: """Get the value of the basic output wanted. Args: name (str): The name of the basic output. Returns: npt.NDArray[ np.float64 ]: the value of the basic output. """ if name == BIOT_COEFFICIENT.attributeName: return self.biotCoefficient elif name == COMPRESSIBILITY.attributeName: return self.compressibility elif name == COMPRESSIBILITY_OED.attributeName: return self.compressibilityOed elif name == COMPRESSIBILITY_REAL.attributeName: return self.compressibilityReal elif name == SPECIFIC_GRAVITY.attributeName: return self.specificGravity elif name == STRESS_EFFECTIVE_RATIO_REAL.attributeName: return self.effectiveStressRatioReal elif name == STRESS_TOTAL.attributeName: return self.totalStress elif name == STRESS_TOTAL_T0.attributeName: return self.totalStressT0 elif name == STRESS_TOTAL_RATIO_REAL.attributeName: return self.totalStressRatioReal elif name == LITHOSTATIC_STRESS.attributeName: return self.lithostaticStress elif name == STRAIN_ELASTIC.attributeName: return self.elasticStrain elif name == STRESS_TOTAL_DELTA.attributeName: return self.deltaTotalStress elif name == RSP_REAL.attributeName: return self.rspReal elif name == RSP_OED.attributeName: return self.rspOed elif name == STRESS_EFFECTIVE_RATIO_OED.attributeName: return self.effectiveStressRatioOed else: raise NameError
[docs] @dataclass class AdvancedOutputValue: """The dataclass with the value of the advanced geomechanics outputs.""" _criticalTotalStressRatio: npt.NDArray[ np.float64 ] | None = None _stressRatioThreshold: npt.NDArray[ np.float64 ] | None = None _criticalPorePressure: npt.NDArray[ np.float64 ] | None = None _criticalPorePressureIndex: npt.NDArray[ np.float64 ] | None = None @property def criticalTotalStressRatio( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the critical total stress ratio value.""" return self._criticalTotalStressRatio @criticalTotalStressRatio.setter def criticalTotalStressRatio( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._criticalTotalStressRatio = value @property def stressRatioThreshold( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the stress ratio threshold value.""" return self._stressRatioThreshold @stressRatioThreshold.setter def stressRatioThreshold( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._stressRatioThreshold = value @property def criticalPorePressure( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the critical pore pressure value.""" return self._criticalPorePressure @criticalPorePressure.setter def criticalPorePressure( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._criticalPorePressure = value @property def criticalPorePressureIndex( self: Self ) -> npt.NDArray[ np.float64 ] | None: """Get the critical pore pressure index value.""" return self._criticalPorePressureIndex @criticalPorePressureIndex.setter def criticalPorePressureIndex( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: self._criticalPorePressureIndex = value
[docs] def getAdvancedOutputValue( self: Self, name: str ) -> npt.NDArray[ np.float64 ] | None: """Get the value of the advanced output wanted. Args: name (str): The name of the advanced output. Returns: npt.NDArray[ np.float64 ]: the value of the advanced output. """ if name == CRITICAL_TOTAL_STRESS_RATIO.attributeName: return self.criticalTotalStressRatio elif name == TOTAL_STRESS_RATIO_THRESHOLD.attributeName: return self.stressRatioThreshold elif name == CRITICAL_PORE_PRESSURE.attributeName: return self.criticalPorePressure elif name == CRITICAL_PORE_PRESSURE_THRESHOLD.attributeName: return self.criticalPorePressureIndex else: raise NameError
physicalConstants: PhysicalConstants _elasticModuli: ElasticModuliValue _mandatoryAttributes: MandatoryAttributesValue _basicOutput: BasicOutputValue _advancedOutput: AdvancedOutputValue def __init__( self: Self, mesh: Union[ vtkPointSet, vtkUnstructuredGrid ], computeAdvancedOutputs: bool = False, speHandler: bool = False, ) -> None: """VTK Filter to perform Geomechanical output computation. Args: mesh (Union[vtkPointsSet, vtkUnstructuredGrid]): Input mesh. computeAdvancedOutputs (bool, optional): True to compute advanced geomechanical parameters, False otherwise. Defaults to False. speHandler (bool, optional): True to use a specific handler, False to use the internal handler. Defaults to False. """ self.output: Union[ vtkPointSet, vtkUnstructuredGrid ] = mesh.NewInstance() self.output.DeepCopy( mesh ) self.doComputeAdvancedOutputs: bool = computeAdvancedOutputs self.physicalConstants = self.PhysicalConstants() self._elasticModuli = self.ElasticModuliValue() self._mandatoryAttributes = self.MandatoryAttributesValue() self._basicOutput = self.BasicOutputValue() self._advancedOutput = self.AdvancedOutputValue() self._attributesToCreate: list[ AttributeEnum ] = [] # Logger. self.logger: Logger if not speHandler: self.logger = getLogger( loggerTitle, True ) else: self.logger = logging.getLogger( loggerTitle ) self.logger.setLevel( logging.INFO )
[docs] def applyFilter( self: Self ) -> bool: """Compute the geomechanical outputs of the mesh. Returns: Boolean (bool): True if calculation successfully ended, False otherwise. """ if not self._checkMandatoryAttributes(): return False if not self.computeBasicOutputs(): return False if self.doComputeAdvancedOutputs and not self.computeAdvancedOutputs(): return False # Create an attribute on the mesh for each geomechanics outputs computed: for attribute in self._attributesToCreate: attributeName: str = attribute.attributeName onPoints: bool = attribute.isOnPoints array: npt.NDArray[ np.float64 ] | None if attribute in ELASTIC_MODULI: array = self._elasticModuli.getElasticModulusValue( attributeName ) elif attribute in BASIC_OUTPUTS: array = self._basicOutput.getBasicOutputValue( attributeName ) elif attribute in ADVANCED_OUTPUTS: array = self._advancedOutput.getAdvancedOutputValue( attributeName ) componentNames: tuple[ str, ...] = () if attribute.nbComponent == 6: componentNames = ComponentNameEnum.XYZ.value if not createAttribute( self.output, array, attributeName, componentNames=componentNames, onPoints=onPoints, logger=self.logger ): return False self.logger.info( "All the geomechanics properties have been added to the mesh." ) self.logger.info( "The filter succeeded." ) return True
[docs] def getOutput( self: Self ) -> Union[ vtkPointSet, vtkUnstructuredGrid ]: """Get the mesh with the geomechanical outputs as attributes. Returns: Union[vtkPointSet, vtkUnstructuredGrid]: The mesh with the geomechanical outputs as attributes. """ return self.output
[docs] def setLoggerHandler( self: Self, handler: logging.Handler ) -> None: """Set a specific handler for the filter logger. In this filter 4 log levels are use, .info, .error, .warning and .critical, be sure to have at least the same 4 levels. Args: handler (logging.Handler): The handler to add. """ if not self.logger.hasHandlers(): self.logger.addHandler( handler ) else: self.logger.warning( "The logger already has an handler, to use yours set the argument 'speHandler' to True during the filter initialization." )
[docs] def getOutputType( self: Self ) -> str: """Get output object type. Returns: str: Type of output object. """ return self.output.GetClassName()
def _checkMandatoryAttributes( self: Self ) -> bool: """Check that mandatory attributes are present in the mesh. The mesh must contains: - The Young modulus and the Poisson's ratio named "youngModulus" and "poissonRatio" or bulk and shear moduli named "bulkModulus" and "shearModulus" - The initial Young modulus and Poisson's ratio named "youngModulusInitial" and "poissonRatioInitial" or the initial bulk modulus named "bulkModulusInitial" - The porosity named "porosity" - The initial porosity named "porosityInitial" - The pressure named "pressure" - The delta of pressure named "deltaPressure" - The density named "density" - The effective stress named "stressEffective" - The initial effective stress named "stressEffectiveInitial" Returns: bool: True if all needed attributes are present, False otherwise """ for elasticModulus in ELASTIC_MODULI: elasticModulusName: str = elasticModulus.attributeName elasticModulusOnPoints: bool = elasticModulus.isOnPoints if isAttributeInObject( self.output, elasticModulusName, elasticModulusOnPoints ): self._elasticModuli.setElasticModulusValue( elasticModulus.attributeName, getArrayInObject( self.output, elasticModulusName, elasticModulusOnPoints ) ) # Check the presence of the elastic moduli at the current time. self.computeYoungPoisson: bool if self._elasticModuli.youngModulus is None and self._elasticModuli.poissonRatio is None: if self._elasticModuli.bulkModulus is not None and self._elasticModuli.shearModulus is not None: self._elasticModuli.youngModulus = fcts.youngModulus( self._elasticModuli.bulkModulus, self._elasticModuli.shearModulus ) self._attributesToCreate.append( YOUNG_MODULUS ) self._elasticModuli.poissonRatio = fcts.poissonRatio( self._elasticModuli.bulkModulus, self._elasticModuli.shearModulus ) self._attributesToCreate.append( POISSON_RATIO ) self.computeYoungPoisson = True else: self.logger.error( f"{ BULK_MODULUS.attributeName } or { SHEAR_MODULUS.attributeName } are missing to compute geomechanical outputs." ) return False elif self._elasticModuli.bulkModulus is None and self._elasticModuli.shearModulus is None: if self._elasticModuli.youngModulus is not None and self._elasticModuli.poissonRatio is not None: self._elasticModuli.bulkModulus = fcts.bulkModulus( self._elasticModuli.youngModulus, self._elasticModuli.poissonRatio ) self._attributesToCreate.append( BULK_MODULUS ) self._elasticModuli.shearModulus = fcts.shearModulus( self._elasticModuli.youngModulus, self._elasticModuli.poissonRatio ) self._attributesToCreate.append( SHEAR_MODULUS ) self.computeYoungPoisson = False else: self.logger.error( f"{ YOUNG_MODULUS.attributeName } or { POISSON_RATIO.attributeName } are missing to compute geomechanical outputs." ) return False else: self.logger.error( f"{ BULK_MODULUS.attributeName } and { SHEAR_MODULUS.attributeName } or { YOUNG_MODULUS.attributeName } and { POISSON_RATIO.attributeName } are mandatory to compute geomechanical outputs." ) return False # Check the presence of the elastic moduli at the initial time. if self._elasticModuli.bulkModulusT0 is None: if self._elasticModuli.youngModulusT0 is not None and self._elasticModuli.poissonRatioT0 is not None: self._elasticModuli.bulkModulusT0 = fcts.bulkModulus( self._elasticModuli.youngModulusT0, self._elasticModuli.poissonRatioT0 ) self._attributesToCreate.append( BULK_MODULUS_T0 ) else: self.logger.error( f"{ BULK_MODULUS_T0.attributeName } or { YOUNG_MODULUS_T0.attributeName } and { POISSON_RATIO_T0.attributeName } are mandatory to compute geomechanical outputs." ) return False # Check the presence of the other mandatory attributes for mandatoryAttribute in MANDATORY_ATTRIBUTES: mandatoryAttributeName: str = mandatoryAttribute.attributeName mandatoryAttributeOnPoints: bool = mandatoryAttribute.isOnPoints if not isAttributeInObject( self.output, mandatoryAttributeName, mandatoryAttributeOnPoints ): self.logger.error( f"The mandatory attribute { mandatoryAttributeName } is missing to compute geomechanical outputs." ) return False else: self._mandatoryAttributes.setMandatoryAttributeValue( mandatoryAttributeName, getArrayInObject( self.output, mandatoryAttributeName, mandatoryAttributeOnPoints ) ) return True
[docs] def computeBasicOutputs( self: Self ) -> bool: """Compute basic geomechanical outputs. Returns: bool: return True if calculation successfully ended, False otherwise. """ if not self._computeBiotCoefficient(): self.logger.error( "Biot coefficient computation failed." ) return False if not self._computeCompressibilityCoefficient(): return False if not self._computeRealEffectiveStressRatio(): self.logger.error( "Effective stress ratio computation failed." ) return False if not self._computeSpecificGravity(): self.logger.error( "Specific gravity computation failed." ) return False # TODO: deactivate lithostatic stress calculation until right formula # if not self.computeLithostaticStress(): # self.logger.error( "Lithostatic stress computation failed." ) # return False if not self._computeTotalStresses(): return False if not self._computeElasticStrain(): self.logger.error( "Elastic strain computation failed." ) return False if not self._computeEffectiveStressRatioOed(): self.logger.error( "Effective stress ration in oedometric condition computation failed." ) return False if not self._computeReservoirStressPathOed(): self.logger.error( "Reservoir stress path in oedometric condition computation failed." ) return False if not self._computeReservoirStressPathReal(): return False self.logger.info( "All geomechanical basic outputs were successfully computed." ) return True
[docs] def computeAdvancedOutputs( self: Self ) -> bool: """Compute advanced geomechanical outputs. Returns: bool: return True if calculation successfully ended, False otherwise. """ if not self._computeCriticalTotalStressRatio(): return False if not self._computeCriticalPorePressure(): return False self.logger.info( "All geomechanical advanced outputs were successfully computed." ) return True
def _computeBiotCoefficient( self: Self ) -> bool: """Compute Biot coefficient from default and grain bulk modulus. Returns: bool: True if calculation successfully ended, False otherwise. """ if not isAttributeInObject( self.output, BIOT_COEFFICIENT.attributeName, BIOT_COEFFICIENT.isOnPoints ): self._basicOutput.biotCoefficient = fcts.biotCoefficient( self.physicalConstants.grainBulkModulus, self._elasticModuli.bulkModulus ) self._attributesToCreate.append( BIOT_COEFFICIENT ) else: self._basicOutput.biotCoefficient = getArrayInObject( self.output, BIOT_COEFFICIENT.attributeName, BIOT_COEFFICIENT.isOnPoints ) self.logger.warning( f"{ BIOT_COEFFICIENT.attributeName } is already on the mesh, it has not been computed by the filter." ) return True def _computeCompressibilityCoefficient( self: Self ) -> bool: """Compute compressibility coefficient from simulation outputs. Compressibility coefficient is computed from Poisson's ratio, bulk modulus, Biot coefficient and Porosity. Returns: bool: True if the attribute is correctly created, False otherwise. """ if not isAttributeInObject( self.output, COMPRESSIBILITY.attributeName, COMPRESSIBILITY.isOnPoints ): self._basicOutput.compressibility = fcts.compressibility( self._elasticModuli.poissonRatio, self._elasticModuli.bulkModulus, self._basicOutput.biotCoefficient, self._mandatoryAttributes.porosity ) self._attributesToCreate.append( COMPRESSIBILITY ) else: self._basicOutput.compressibility = getArrayInObject( self.output, COMPRESSIBILITY.attributeName, COMPRESSIBILITY.isOnPoints ) self.logger.warning( f"{ COMPRESSIBILITY.attributeName } is already on the mesh, it has not been computed by the filter." ) # oedometric compressibility if not isAttributeInObject( self.output, COMPRESSIBILITY_OED.attributeName, COMPRESSIBILITY_OED.isOnPoints ): self._basicOutput.compressibilityOed = fcts.compressibilityOed( self._elasticModuli.shearModulus, self._elasticModuli.bulkModulus, self._mandatoryAttributes.porosity ) self._attributesToCreate.append( COMPRESSIBILITY_OED ) else: self._basicOutput.compressibilityOed = getArrayInObject( self.output, COMPRESSIBILITY_OED.attributeName, COMPRESSIBILITY_OED.isOnPoints ) self.logger.warning( f"{ COMPRESSIBILITY_OED.attributeName } is already on the mesh, it has not been computed by the filter." ) # real compressibility if not isAttributeInObject( self.output, COMPRESSIBILITY_REAL.attributeName, COMPRESSIBILITY_REAL.isOnPoints ): self._basicOutput.compressibilityReal = fcts.compressibilityReal( self._mandatoryAttributes.deltaPressure, self._mandatoryAttributes.porosity, self._mandatoryAttributes.porosityInitial ) self._attributesToCreate.append( COMPRESSIBILITY_REAL ) else: self._basicOutput.compressibilityReal = getArrayInObject( self.output, COMPRESSIBILITY_REAL.attributeName, COMPRESSIBILITY_REAL.isOnPoints ) self.logger.warning( f"{ COMPRESSIBILITY_REAL.attributeName } is already on the mesh, it has not been computed by the filter." ) return True def _computeSpecificGravity( self: Self ) -> bool: """Create Specific gravity attribute. Specific gravity is computed from rock density attribute and specific density input. Returns: bool: True if the attribute is correctly created, False otherwise. """ if not isAttributeInObject( self.output, SPECIFIC_GRAVITY.attributeName, SPECIFIC_GRAVITY.isOnPoints ): self._basicOutput.specificGravity = fcts.specificGravity( self._mandatoryAttributes.density, self.physicalConstants.specificDensity ) self._attributesToCreate.append( SPECIFIC_GRAVITY ) else: self._basicOutput.specificGravity = getArrayInObject( self.output, SPECIFIC_GRAVITY.attributeName, SPECIFIC_GRAVITY.isOnPoints ) self.logger.warning( f"{ SPECIFIC_GRAVITY.attributeName } is already on the mesh, it has not been computed by the filter." ) return True def _doComputeStressRatioReal( self: Self, stress: npt.NDArray[ np.float64 ], basicOutput: AttributeEnum ) -> npt.NDArray[ np.float64 ]: """Compute the ratio between horizontal and vertical effective stress. Returns: bool: return True if calculation successfully ended, False otherwise. """ verticalStress: npt.NDArray[ np.float64 ] = stress[ :, 2 ] # keep the minimum of the 2 horizontal components horizontalStress: npt.NDArray[ np.float64 ] = np.min( stress[ :, :2 ], axis=1 ) stressRatioReal: npt.NDArray[ np.float64 ] if not isAttributeInObject( self.output, basicOutput.attributeName, basicOutput.isOnPoints ): stressRatioReal = fcts.stressRatio( horizontalStress, verticalStress ) self._attributesToCreate.append( basicOutput ) else: stressRatioReal = getArrayInObject( self.output, basicOutput.attributeName, basicOutput.isOnPoints ) self.logger.warning( f"{ basicOutput.attributeName } is already on the mesh, it has not been computed by the filter." ) return stressRatioReal def _computeRealEffectiveStressRatio( self: Self ) -> bool: """Compute effective stress ratio. Returns: bool: True if calculation successfully ended, False otherwise. """ if self._mandatoryAttributes.effectiveStress is not None: self._basicOutput.effectiveStressRatioReal = self._doComputeStressRatioReal( self._mandatoryAttributes.effectiveStress, STRESS_EFFECTIVE_RATIO_REAL ) return True else: self.logger.error( f"{ STRESS_EFFECTIVE_RATIO_REAL.attributeName } has not been computed, mandatory attribute { STRESS_EFFECTIVE.attributeName } is missing." ) return False def _doComputeTotalStress( self: Self, effectiveStress: npt.NDArray[ np.float64 ], pressure: Union[ npt.NDArray[ np.float64 ], None ], biotCoefficient: npt.NDArray[ np.float64 ], ) -> npt.NDArray[ np.float64 ]: """Compute total stress from effective stress and pressure. Args: effectiveStress (npt.NDArray[np.float64]): Effective stress. pressure (npt.NDArray[np.float64] | None): Pressure. biotCoefficient (npt.NDArray[np.float64]): Biot coefficient. Returns: npt.NDArray[ np.float64 ]: TotalStress. """ totalStress: npt.NDArray[ np.float64 ] if pressure is None: totalStress = np.copy( effectiveStress ) self.logger.warning( "Pressure attribute is undefined, total stress will be equal to effective stress." ) else: if np.isnan( pressure ).any(): self.logger.warning( "Some cells do not have pressure data, for those cells, pressure is set to 0." ) pressure[ np.isnan( pressure ) ] = 0.0 totalStress = fcts.totalStress( effectiveStress, biotCoefficient, pressure ) return totalStress def _computeTotalStressInitial( self: Self ) -> bool: """Compute total stress at initial time step. Returns: bool: True if calculation successfully ended, False otherwise. """ # Compute Biot at initial time step. biotCoefficientT0: npt.NDArray[ np.float64 ] = fcts.biotCoefficient( self.physicalConstants.grainBulkModulus, self._elasticModuli.bulkModulusT0 ) pressureT0: npt.NDArray[ np.float64 ] | None = None # Case pressureT0 is None, total stress = effective stress # (managed by doComputeTotalStress function) if self._mandatoryAttributes.pressure is not None: if self._mandatoryAttributes.deltaPressure is not None: pressureT0 = self._mandatoryAttributes.pressure - self._mandatoryAttributes.deltaPressure else: self.logger.error( f"Mandatory attribute { DELTA_PRESSURE.attributeName } is missing." ) return False if not isAttributeInObject( self.output, STRESS_TOTAL_T0.attributeName, STRESS_TOTAL_T0.isOnPoints ): if self._mandatoryAttributes.effectiveStressT0 is not None: self._basicOutput.totalStressT0 = self._doComputeTotalStress( self._mandatoryAttributes.effectiveStressT0, pressureT0, biotCoefficientT0 ) self._attributesToCreate.append( STRESS_TOTAL_T0 ) else: self.logger.error( f"{ STRESS_TOTAL_T0.attributeName } has not been computed, mandatory attribute { STRESS_EFFECTIVE_T0.attributeName } is missing." ) return False else: self._basicOutput.totalStressT0 = getArrayInObject( self.output, STRESS_TOTAL_T0.attributeName, STRESS_TOTAL_T0.isOnPoints ) self.logger.warning( f"{ STRESS_TOTAL_T0.attributeName } is already on the mesh, it has not been computed by the filter." ) return True def _computeTotalStresses( self: Self ) -> bool: """Compute total stress total stress ratio. Total stress is computed at the initial and current time steps. Total stress ratio is computed at current time step only. Returns: bool: True if calculation successfully ended, False otherwise. """ # Compute total stress at initial time step. if not self._computeTotalStressInitial(): self.logger.error( "Total stress at initial time step computation failed." ) return False # Compute total stress at current time step. if not isAttributeInObject( self.output, STRESS_TOTAL.attributeName, STRESS_TOTAL.isOnPoints ): if self._mandatoryAttributes.effectiveStress is not None and self._basicOutput.biotCoefficient is not None: self._basicOutput.totalStress = self._doComputeTotalStress( self._mandatoryAttributes.effectiveStress, self._mandatoryAttributes.pressure, self._basicOutput.biotCoefficient ) self._attributesToCreate.append( STRESS_TOTAL ) else: self.logger.error( f"{ STRESS_TOTAL.attributeName } has not been computed, mandatory attributes { STRESS_EFFECTIVE.attributeName } or { BIOT_COEFFICIENT.attributeName } are missing." ) return False else: self._basicOutput.totalStress = getArrayInObject( self.output, STRESS_TOTAL.attributeName, STRESS_TOTAL.isOnPoints ) self.logger.warning( f"{ STRESS_TOTAL.attributeName } is already on the mesh, it has not been computed by the filter." ) # Compute total stress ratio. if self._basicOutput.totalStress is not None: self._basicOutput.totalStressRatioReal = self._doComputeStressRatioReal( self._basicOutput.totalStress, STRESS_TOTAL_RATIO_REAL ) else: self.logger.error( f"{ STRESS_TOTAL_RATIO_REAL.attributeName } has not been computed, Mandatory attribute { BIOT_COEFFICIENT.attributeName } is missing." ) return False return True
[docs] def computeLithostaticStress( self: Self ) -> bool: """Compute lithostatic stress. Returns: bool: True if calculation successfully ended, False otherwise. """ if not isAttributeInObject( self.output, LITHOSTATIC_STRESS.attributeName, LITHOSTATIC_STRESS.isOnPoints ): depth: npt.NDArray[ np.float64 ] = self._doComputeDepthAlongLine( ) if LITHOSTATIC_STRESS.isOnPoints else self._doComputeDepthInMesh() self._basicOutput.lithostaticStress = fcts.lithostaticStress( depth, self._mandatoryAttributes.density, GRAVITY ) self._attributesToCreate.append( LITHOSTATIC_STRESS ) else: self._basicOutput.lithostaticStress = getArrayInObject( self.output, LITHOSTATIC_STRESS.attributeName, LITHOSTATIC_STRESS.isOnPoints ) self.logger.warning( f"{ LITHOSTATIC_STRESS.attributeName } is already on the mesh, it has not been computed by the filter." ) return True
def _doComputeDepthAlongLine( self: Self ) -> npt.NDArray[ np.float64 ]: """Compute depth along a line. Returns: npt.NDArray[np.float64]: 1D array with depth property """ # get z coordinate zCoord: npt.NDArray[ np.float64 ] = self._getZcoordinates( True ) assert zCoord is not None, "Depth coordinates cannot be computed." # TODO: to find how to compute depth in a general case # GEOS z axis is upward depth: npt.NDArray[ np.float64 ] = -1.0 * zCoord return depth def _doComputeDepthInMesh( self: Self ) -> npt.NDArray[ np.float64 ]: """Compute depth of each cell in a mesh. Returns: npt.NDArray[np.float64]: array with depth property """ # get z coordinate zCoord: npt.NDArray[ np.float64 ] = self._getZcoordinates( False ) assert zCoord is not None, "Depth coordinates cannot be computed." # TODO: to find how to compute depth in a general case depth: npt.NDArray[ np.float64 ] = -1.0 * zCoord return depth def _getZcoordinates( self: Self, onPoints: bool ) -> npt.NDArray[ np.float64 ]: """Get z coordinates from self.output. Args: onPoints (bool): True if the attribute is on points, False if it is on cells. Returns: npt.NDArray[np.float64]: 1D array with depth property """ # get z coordinate zCoord: npt.NDArray[ np.float64 ] pointCoords: npt.NDArray[ np.float64 ] = self._getPointCoordinates( onPoints ) assert pointCoords is not None, "Point coordinates are undefined." assert pointCoords.shape[ 1 ] == 2, "Point coordinates are undefined." zCoord = pointCoords[ :, 2 ] return zCoord def _getPointCoordinates( self: Self, onPoints: bool ) -> npt.NDArray[ np.float64 ]: """Get the coordinates of Points or Cell center. Args: onPoints (bool): True if the attribute is on points, False if it is on cells. Returns: npt.NDArray[np.float64]: points/cell center coordinates """ if onPoints: return self.output.GetPoints() # type: ignore[no-any-return] else: # Find cell centers filter = vtkCellCenters() filter.SetInputDataObject( self.output ) filter.Update() return filter.GetOutput().GetPoints() # type: ignore[no-any-return] def _computeElasticStrain( self: Self ) -> bool: """Compute elastic strain from effective stress and elastic modulus. Returns: bool: return True if calculation successfully ended, False otherwise. """ if self._mandatoryAttributes.effectiveStress is not None and self._mandatoryAttributes.effectiveStressT0 is not None: deltaEffectiveStress = self._mandatoryAttributes.effectiveStress - self._mandatoryAttributes.effectiveStressT0 if not isAttributeInObject( self.output, STRAIN_ELASTIC.attributeName, STRAIN_ELASTIC.isOnPoints ): if self.computeYoungPoisson: self._basicOutput.elasticStrain = fcts.elasticStrainFromBulkShear( deltaEffectiveStress, self._elasticModuli.bulkModulus, self._elasticModuli.shearModulus ) else: self._basicOutput.elasticStrain = fcts.elasticStrainFromYoungPoisson( deltaEffectiveStress, self._elasticModuli.youngModulus, self._elasticModuli.poissonRatio ) self._attributesToCreate.append( STRAIN_ELASTIC ) else: self._basicOutput.totalStressT0 = getArrayInObject( self.output, STRAIN_ELASTIC.attributeName, STRAIN_ELASTIC.isOnPoints ) self.logger.warning( f"{ STRAIN_ELASTIC.attributeName } is already on the mesh, it has not been computed by the filter." ) return True else: self.logger.error( f"{ STRAIN_ELASTIC.attributeName } has not been computed, mandatory attributes { STRESS_EFFECTIVE.attributeName } or { STRESS_EFFECTIVE_T0.attributeName } are missing." ) return False def _computeReservoirStressPathReal( self: Self ) -> bool: """Compute reservoir stress paths. Returns: bool: True if calculation successfully ended, False otherwise. """ # create delta stress attribute for QC if not isAttributeInObject( self.output, STRESS_TOTAL_DELTA.attributeName, STRESS_TOTAL_DELTA.isOnPoints ): if self._basicOutput.totalStress is not None and self._basicOutput.totalStressT0 is not None: self._basicOutput.deltaTotalStress = self._basicOutput.totalStress - self._basicOutput.totalStressT0 self._attributesToCreate.append( STRESS_TOTAL_DELTA ) else: self.logger.error( f"{ STRESS_TOTAL_DELTA.attributeName } has not been computed, mandatory attributes { STRESS_TOTAL.attributeName } or { STRESS_TOTAL_T0.attributeName } are missing." ) return False else: self._basicOutput.deltaTotalStress = getArrayInObject( self.output, STRESS_TOTAL_DELTA.attributeName, STRESS_TOTAL_DELTA.isOnPoints ) self.logger.warning( f"{ STRESS_TOTAL_DELTA.attributeName } is already on the mesh, it has not been computed by the filter." ) if not isAttributeInObject( self.output, RSP_REAL.attributeName, RSP_REAL.isOnPoints ): self._basicOutput.rspReal = fcts.reservoirStressPathReal( self._basicOutput.deltaTotalStress, self._mandatoryAttributes.deltaPressure ) self._attributesToCreate.append( RSP_REAL ) else: self._basicOutput.rspReal = getArrayInObject( self.output, RSP_REAL.attributeName, RSP_REAL.isOnPoints ) self.logger.warning( f"{ RSP_REAL.attributeName } is already on the mesh, it has not been computed by the filter." ) return True def _computeReservoirStressPathOed( self: Self ) -> bool: """Compute Reservoir Stress Path in oedometric conditions. Returns: bool: return True if calculation successfully ended, False otherwise. """ if not isAttributeInObject( self.output, RSP_OED.attributeName, RSP_OED.isOnPoints ): self._basicOutput.rspOed = fcts.reservoirStressPathOed( self._basicOutput.biotCoefficient, self._elasticModuli.poissonRatio ) self._attributesToCreate.append( RSP_OED ) else: self._basicOutput.rspOed = getArrayInObject( self.output, RSP_OED.attributeName, RSP_OED.isOnPoints ) self.logger.warning( f"{ RSP_OED.attributeName } is already on the mesh, it has not been computed by the filter." ) return True def _computeEffectiveStressRatioOed( self: Self ) -> bool: """Compute the effective stress ratio in oedometric conditions. Returns: bool: True if calculation successfully ended, False otherwise. """ if not isAttributeInObject( self.output, STRESS_EFFECTIVE_RATIO_OED.attributeName, STRESS_EFFECTIVE_RATIO_OED.isOnPoints ): self._basicOutput.effectiveStressRatioOed = fcts.deviatoricStressPathOed( self._elasticModuli.poissonRatio ) self._attributesToCreate.append( STRESS_EFFECTIVE_RATIO_OED ) else: self._basicOutput.effectiveStressRatioOed = getArrayInObject( self.output, STRESS_EFFECTIVE_RATIO_OED.attributeName, STRESS_EFFECTIVE_RATIO_OED.isOnPoints ) self.logger.warning( f"{ STRESS_EFFECTIVE_RATIO_OED.attributeName } is already on the mesh, it has not been computed by the filter." ) return True def _computeCriticalTotalStressRatio( self: Self ) -> bool: """Compute fracture index and fracture threshold. Returns: bool: return True if calculation successfully ended, False otherwise. """ if not isAttributeInObject( self.output, CRITICAL_TOTAL_STRESS_RATIO.attributeName, CRITICAL_TOTAL_STRESS_RATIO.isOnPoints ): if self._basicOutput.totalStress is not None: verticalStress: npt.NDArray[ np.float64 ] = self._basicOutput.totalStress[ :, 2 ] self._advancedOutput.criticalTotalStressRatio = fcts.criticalTotalStressRatio( self._mandatoryAttributes.pressure, verticalStress ) self._attributesToCreate.append( CRITICAL_TOTAL_STRESS_RATIO ) else: self.logger.error( f"{ CRITICAL_TOTAL_STRESS_RATIO.attributeName } has not been computed, mandatory attribute { STRESS_TOTAL.attributeName } is missing." ) return False else: self._advancedOutput.criticalTotalStressRatio = getArrayInObject( self.output, CRITICAL_TOTAL_STRESS_RATIO.attributeName, CRITICAL_TOTAL_STRESS_RATIO.isOnPoints ) self.logger.warning( f"{ CRITICAL_TOTAL_STRESS_RATIO.attributeName } is already on the mesh, it has not been computed by the filter." ) if not isAttributeInObject( self.output, TOTAL_STRESS_RATIO_THRESHOLD.attributeName, TOTAL_STRESS_RATIO_THRESHOLD.isOnPoints ): if self._basicOutput.totalStress is not None: mask: npt.NDArray[ np.bool_ ] = np.argmin( np.abs( self._basicOutput.totalStress[ :, :2 ] ), axis=1 ) horizontalStress: npt.NDArray[ np.float64 ] = self._basicOutput.totalStress[ :, :2 ][ np.arange( self._basicOutput.totalStress[ :, :2 ].shape[ 0 ] ), mask ] self._advancedOutput.stressRatioThreshold = fcts.totalStressRatioThreshold( self._mandatoryAttributes.pressure, horizontalStress ) self._attributesToCreate.append( TOTAL_STRESS_RATIO_THRESHOLD ) else: self.logger.error( f"{ TOTAL_STRESS_RATIO_THRESHOLD.attributeName } has not been computed, mandatory attribute { STRESS_TOTAL.attributeName } is missing." ) return False else: self._advancedOutput.stressRatioThreshold = getArrayInObject( self.output, TOTAL_STRESS_RATIO_THRESHOLD.attributeName, TOTAL_STRESS_RATIO_THRESHOLD.isOnPoints ) self.logger.warning( f"{ TOTAL_STRESS_RATIO_THRESHOLD.attributeName } is already on the mesh, it has not been computed by the filter." ) return True def _computeCriticalPorePressure( self: Self ) -> bool: """Compute the critical pore pressure and the pressure index. Returns: bool: return True if calculation successfully ended, False otherwise. """ if not isAttributeInObject( self.output, CRITICAL_PORE_PRESSURE.attributeName, CRITICAL_PORE_PRESSURE.isOnPoints ): if self._basicOutput.totalStress is not None: self._advancedOutput.criticalPorePressure = fcts.criticalPorePressure( -1.0 * self._basicOutput.totalStress, self.physicalConstants.rockCohesion, self.physicalConstants.frictionAngle ) self._attributesToCreate.append( CRITICAL_PORE_PRESSURE ) else: self.logger.error( f"{ CRITICAL_PORE_PRESSURE.attributeName } has not been computed, mandatory attribute { STRESS_TOTAL.attributeName } is missing." ) return False else: self._advancedOutput.criticalPorePressure = getArrayInObject( self.output, CRITICAL_PORE_PRESSURE.attributeName, CRITICAL_PORE_PRESSURE.isOnPoints ) self.logger.warning( f"{ CRITICAL_PORE_PRESSURE.attributeName } is already on the mesh, it has not been computed by the filter." ) # Add critical pore pressure index (i.e., ratio between pressure and criticalPorePressure) if not isAttributeInObject( self.output, CRITICAL_PORE_PRESSURE_THRESHOLD.attributeName, CRITICAL_PORE_PRESSURE_THRESHOLD.isOnPoints ): self._advancedOutput.criticalPorePressureIndex = fcts.criticalPorePressureThreshold( self._mandatoryAttributes.pressure, self._advancedOutput.criticalPorePressure ) self._attributesToCreate.append( CRITICAL_PORE_PRESSURE_THRESHOLD ) else: self._advancedOutput.criticalPorePressureIndex = getArrayInObject( self.output, CRITICAL_PORE_PRESSURE_THRESHOLD.attributeName, CRITICAL_PORE_PRESSURE_THRESHOLD.isOnPoints ) self.logger.warning( f"{ CRITICAL_PORE_PRESSURE_THRESHOLD.attributeName } is already on the mesh, it has not been computed by the filter." ) return True