From 6c152c2b7c2591f5eb0506953420fccf9d92cce9 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Thu, 27 Mar 2025 11:35:53 +0100
Subject: [PATCH 01/14] Move geomechanics processes to dedicated package

---
 .../geos/geomechanics/model}/MohrCircle.py    |  238 +-
 .../geos/geomechanics/model}/MohrCoulomb.py   |  198 +-
 .../geomechanicsCalculatorFunctions.py        | 2022 ++++++++---------
 3 files changed, 1229 insertions(+), 1229 deletions(-)
 rename {geos-posp/src/geos_posp/processing => geos-geomechanics/src/geos/geomechanics/model}/MohrCircle.py (95%)
 rename {geos-posp/src/geos_posp/processing => geos-geomechanics/src/geos/geomechanics/model}/MohrCoulomb.py (97%)
 rename {geos-posp/src/geos_posp => geos-geomechanics/src/geos/geomechanics}/processing/geomechanicsCalculatorFunctions.py (97%)

diff --git a/geos-posp/src/geos_posp/processing/MohrCircle.py b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
similarity index 95%
rename from geos-posp/src/geos_posp/processing/MohrCircle.py
rename to geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
index 5ec2cef9..f32ccf44 100644
--- a/geos-posp/src/geos_posp/processing/MohrCircle.py
+++ b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
@@ -1,119 +1,119 @@
-# SPDX-License-Identifier: Apache-2.0
-# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Alexandre Benedicto, Martin Lemay
-import numpy as np
-import numpy.typing as npt
-from typing_extensions import Self
-
-from geos_posp.processing.geomechanicsCalculatorFunctions import (
-    computeStressPrincipalComponentsFromStressVector,
-)
-
-__doc__ = """
-MohrCircle module define the Mohr's circle parameters.
-
-Inputs are a 6 component stress vector, a circle id, and the mechanical
-convention used for compression.
-The class computes principal components from stress vector during initialization.
-Accessors get access to these 3 principal components as well as circle center
-and radius.
-
-To use the object:
-
-.. code-block:: python
-
-    from processing.MohrCircle import MohrCircle
-
-    # create the object
-    stressVector :npt.NDArray[np.float64]
-    circleId :str
-    mohrCircle :MohrCircle = MohrCircle(circleId)
-
-    # either directly set principal components (p3 <= p2 <= p1)
-    mohrCircle.SetPrincipalComponents(p3, p2, p1)
-    # or compute them from stress vector
-    mohrCircle.computePrincipalComponents(stressVector)
-
-    # access to members
-    id :str = mohrCircle.getCircleId()
-    p1, p2, p3 :float = mohrCircle.getPrincipalComponents()
-    radius :float = mohrCircle.getCircleRadius()
-    center :float = mohrCircle.getCircleCenter()
-"""
-
-
-class MohrCircle:
-    def __init__(self: Self, circleId: str) -> None:
-        """Compute Mohr's Circle from input stress.
-
-        Args:
-            circleId (str): Mohr's circle id.
-        """
-        self.m_circleId: str = circleId
-
-        self.m_p1: float = 0.0
-        self.m_p2: float = 0.0
-        self.m_p3: float = 0.0
-
-    def __str__(self: Self) -> str:
-        """Overload of __str__ method."""
-        return self.m_circleId
-
-    def __repr__(self: Self) -> str:
-        """Overload of __repr__ method."""
-        return self.m_circleId
-
-    def setCircleId(self: Self, circleId: str) -> None:
-        """Set circle Id variable.
-
-        Args:
-            circleId (str): circle Id.
-        """
-        self.m_circleId = circleId
-
-    def getCircleId(self: Self) -> str:
-        """Access the Id of the Mohr circle.
-
-        Returns:
-            str: Id of the Mohr circle
-        """
-        return self.m_circleId
-
-    def getCircleRadius(self: Self) -> float:
-        """Compute and return Mohr's circle radius from principal components."""
-        return (self.m_p1 - self.m_p3) / 2.0
-
-    def getCircleCenter(self: Self) -> float:
-        """Compute and return Mohr's circle center from principal components."""
-        return (self.m_p1 + self.m_p3) / 2.0
-
-    def getPrincipalComponents(self: Self) -> tuple[float, float, float]:
-        """Get Moh's circle principal components."""
-        return (self.m_p3, self.m_p2, self.m_p1)
-
-    def setPrincipalComponents(self: Self, p3: float, p2: float, p1: float) -> None:
-        """Set principal components.
-
-        Args:
-            p3 (float): first component. Must be the lowest.
-            p2 (float): second component.
-            p1 (float): third component. Must be the greatest.
-        """
-        assert (p3 <= p2) and (p2 <= p1), "Component order is wrong."
-        self.m_p3 = p3
-        self.m_p2 = p2
-        self.m_p1 = p1
-
-    def computePrincipalComponents(
-        self: Self, stressVector: npt.NDArray[np.float64]
-    ) -> None:
-        """Calculate principal components.
-
-        Args:
-            stressVector (npt.NDArray[np.float64]): 6 components stress vector
-                Stress vector must follow GEOS convention (XX, YY, ZZ, YZ, XZ, XY)
-        """
-        # get stress principal components
-        self.m_p3, self.m_p2, self.m_p1 = (
-            computeStressPrincipalComponentsFromStressVector(stressVector)
-        )
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Alexandre Benedicto, Martin Lemay
+import numpy as np
+import numpy.typing as npt
+from typing_extensions import Self
+
+from geos_geomechanics.processing.geomechanicsCalculatorFunctions import (
+    computeStressPrincipalComponentsFromStressVector,
+)
+
+__doc__ = """
+MohrCircle module define the Mohr's circle parameters.
+
+Inputs are a 6 component stress vector, a circle id, and the mechanical
+convention used for compression.
+The class computes principal components from stress vector during initialization.
+Accessors get access to these 3 principal components as well as circle center
+and radius.
+
+To use the object:
+
+.. code-block:: python
+
+    from processing.MohrCircle import MohrCircle
+
+    # create the object
+    stressVector :npt.NDArray[np.float64]
+    circleId :str
+    mohrCircle :MohrCircle = MohrCircle(circleId)
+
+    # either directly set principal components (p3 <= p2 <= p1)
+    mohrCircle.SetPrincipalComponents(p3, p2, p1)
+    # or compute them from stress vector
+    mohrCircle.computePrincipalComponents(stressVector)
+
+    # access to members
+    id :str = mohrCircle.getCircleId()
+    p1, p2, p3 :float = mohrCircle.getPrincipalComponents()
+    radius :float = mohrCircle.getCircleRadius()
+    center :float = mohrCircle.getCircleCenter()
+"""
+
+
+class MohrCircle:
+    def __init__(self: Self, circleId: str) -> None:
+        """Compute Mohr's Circle from input stress.
+
+        Args:
+            circleId (str): Mohr's circle id.
+        """
+        self.m_circleId: str = circleId
+
+        self.m_p1: float = 0.0
+        self.m_p2: float = 0.0
+        self.m_p3: float = 0.0
+
+    def __str__(self: Self) -> str:
+        """Overload of __str__ method."""
+        return self.m_circleId
+
+    def __repr__(self: Self) -> str:
+        """Overload of __repr__ method."""
+        return self.m_circleId
+
+    def setCircleId(self: Self, circleId: str) -> None:
+        """Set circle Id variable.
+
+        Args:
+            circleId (str): circle Id.
+        """
+        self.m_circleId = circleId
+
+    def getCircleId(self: Self) -> str:
+        """Access the Id of the Mohr circle.
+
+        Returns:
+            str: Id of the Mohr circle
+        """
+        return self.m_circleId
+
+    def getCircleRadius(self: Self) -> float:
+        """Compute and return Mohr's circle radius from principal components."""
+        return (self.m_p1 - self.m_p3) / 2.0
+
+    def getCircleCenter(self: Self) -> float:
+        """Compute and return Mohr's circle center from principal components."""
+        return (self.m_p1 + self.m_p3) / 2.0
+
+    def getPrincipalComponents(self: Self) -> tuple[float, float, float]:
+        """Get Moh's circle principal components."""
+        return (self.m_p3, self.m_p2, self.m_p1)
+
+    def setPrincipalComponents(self: Self, p3: float, p2: float, p1: float) -> None:
+        """Set principal components.
+
+        Args:
+            p3 (float): first component. Must be the lowest.
+            p2 (float): second component.
+            p1 (float): third component. Must be the greatest.
+        """
+        assert (p3 <= p2) and (p2 <= p1), "Component order is wrong."
+        self.m_p3 = p3
+        self.m_p2 = p2
+        self.m_p1 = p1
+
+    def computePrincipalComponents(
+        self: Self, stressVector: npt.NDArray[np.float64]
+    ) -> None:
+        """Calculate principal components.
+
+        Args:
+            stressVector (npt.NDArray[np.float64]): 6 components stress vector
+                Stress vector must follow GEOS convention (XX, YY, ZZ, YZ, XZ, XY)
+        """
+        # get stress principal components
+        self.m_p3, self.m_p2, self.m_p1 = (
+            computeStressPrincipalComponentsFromStressVector(stressVector)
+        )
\ No newline at end of file
diff --git a/geos-posp/src/geos_posp/processing/MohrCoulomb.py b/geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py
similarity index 97%
rename from geos-posp/src/geos_posp/processing/MohrCoulomb.py
rename to geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py
index 4072191b..4ed0efcf 100644
--- a/geos-posp/src/geos_posp/processing/MohrCoulomb.py
+++ b/geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py
@@ -1,99 +1,99 @@
-# SPDX-License-Identifier: Apache-2.0
-# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Alexandre Benedicto, Martin Lemay
-import numpy as np
-import numpy.typing as npt
-from typing_extensions import Self, Union
-
-__doc__ = """
-MohrCoulomb module define the Mohr-Coulomb failure envelop class.
-
-Inputs are the rock cohesion (Pa) and the friction angle (°).
-2 methods allow to compute either shear stress values according to normal stress
-values, or the failure envelope including normal stress and corresponding shear
-stress values.
-
-To use the object:
-
-.. code-block:: python
-
-    from processing.MohrCoulomb import MohrCoulomb
-
-    # create the object
-    rockCohesion :float = 1.5e9 # Pa
-    frictionAngle :float = 10 # degree
-    mohrCoulomb = MohrCoulomb(rockCohesion, frictionAngle)
-
-    # compute shear stress values according to normal stress values
-    normalStress :npt.NDArray[np.float64] = np.linspace(1e9, 1.5e9)
-    shearStress :npt.NDArray[np.float64] = mohrCoulomb.computeShearStress(normalStress)
-
-    # compute the failure envelope including normal stress and corresponding shear
-    # stress values
-    # ones may also define minimum normal stress and the number of points
-    normalStressMax :float = 1.5e9
-    normalStress, shearStress = mohrCoulomb.computeFailureEnvelop(normalStressMax)
-"""
-
-
-class MohrCoulomb:
-    def __init__(self: Self, rockCohesion: float, frictionAngle: float) -> None:
-        """Define Mohr-Coulomb failure envelop.
-
-        Args:
-            rockCohesion (float): rock cohesion (Pa).
-            frictionAngle (float): friction angle (rad).
-        """
-        # rock cohesion
-        self.m_rockCohesion = rockCohesion
-        # failure envelop slope
-        self.m_slope: float = np.tan(frictionAngle)
-        # intersection of failure envelop and abscissa axis
-        self.m_sigmaMin: float = -rockCohesion / self.m_slope
-
-    def computeShearStress(
-        self: Self, stressNormal0: Union[float, npt.NDArray[np.float64]]
-    ) -> Union[float, npt.NDArray[np.float64]]:
-        """Compute shear stress from normal stress.
-
-        Args:
-            stressNormal0 (float | npt.NDArray[np.float64]): normal stress
-                value (Pa)
-
-        Returns:
-            float | npt.NDArray[np.float64]): shear stress value.
-        """
-        stressNormal: npt.NDArray[np.float64] = np.array(stressNormal0)
-        tau: npt.NDArray[np.float64] = self.m_slope * stressNormal + self.m_rockCohesion
-        return tau
-
-    def computeFailureEnvelop(
-        self: Self,
-        stressNormalMax: float,
-        stressNormalMin: Union[float, None] = None,
-        n: int = 100,
-    ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
-        """Compute the envelope failure between min and max normal stress.
-
-        Args:
-            stressNormalMax (float): Maximum normal stress (Pa)
-            stressNormalMin (float | None, optional): Minimum normal stress.
-                If it is None, the envelop is computed from the its intersection
-                with the abscissa axis.
-
-                Defaults to None.
-            n (int, optional): Number of points to define the envelop.
-            Defaults to 100.
-
-        Returns:
-            tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: failure
-            envelop coordinates where first element is the abscissa and second
-            element is the ordinates.
-        """
-        sigmaMin: float = (
-            self.m_sigmaMin if stressNormalMin is None else stressNormalMin
-        )
-        stressNormal: npt.NDArray[np.float64] = np.linspace(
-            sigmaMin, stressNormalMax, n
-        )
-        return (stressNormal, np.array(self.computeShearStress(stressNormal)))
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Alexandre Benedicto, Martin Lemay
+import numpy as np
+import numpy.typing as npt
+from typing_extensions import Self, Union
+
+__doc__ = """
+MohrCoulomb module define the Mohr-Coulomb failure envelop class.
+
+Inputs are the rock cohesion (Pa) and the friction angle (°).
+2 methods allow to compute either shear stress values according to normal stress
+values, or the failure envelope including normal stress and corresponding shear
+stress values.
+
+To use the object:
+
+.. code-block:: python
+
+    from processing.MohrCoulomb import MohrCoulomb
+
+    # create the object
+    rockCohesion :float = 1.5e9 # Pa
+    frictionAngle :float = 10 # degree
+    mohrCoulomb = MohrCoulomb(rockCohesion, frictionAngle)
+
+    # compute shear stress values according to normal stress values
+    normalStress :npt.NDArray[np.float64] = np.linspace(1e9, 1.5e9)
+    shearStress :npt.NDArray[np.float64] = mohrCoulomb.computeShearStress(normalStress)
+
+    # compute the failure envelope including normal stress and corresponding shear
+    # stress values
+    # ones may also define minimum normal stress and the number of points
+    normalStressMax :float = 1.5e9
+    normalStress, shearStress = mohrCoulomb.computeFailureEnvelop(normalStressMax)
+"""
+
+
+class MohrCoulomb:
+    def __init__(self: Self, rockCohesion: float, frictionAngle: float) -> None:
+        """Define Mohr-Coulomb failure envelop.
+
+        Args:
+            rockCohesion (float): rock cohesion (Pa).
+            frictionAngle (float): friction angle (rad).
+        """
+        # rock cohesion
+        self.m_rockCohesion = rockCohesion
+        # failure envelop slope
+        self.m_slope: float = np.tan(frictionAngle)
+        # intersection of failure envelop and abscissa axis
+        self.m_sigmaMin: float = -rockCohesion / self.m_slope
+
+    def computeShearStress(
+        self: Self, stressNormal0: Union[float, npt.NDArray[np.float64]]
+    ) -> Union[float, npt.NDArray[np.float64]]:
+        """Compute shear stress from normal stress.
+
+        Args:
+            stressNormal0 (float | npt.NDArray[np.float64]): normal stress
+                value (Pa)
+
+        Returns:
+            float | npt.NDArray[np.float64]): shear stress value.
+        """
+        stressNormal: npt.NDArray[np.float64] = np.array(stressNormal0)
+        tau: npt.NDArray[np.float64] = self.m_slope * stressNormal + self.m_rockCohesion
+        return tau
+
+    def computeFailureEnvelop(
+        self: Self,
+        stressNormalMax: float,
+        stressNormalMin: Union[float, None] = None,
+        n: int = 100,
+    ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
+        """Compute the envelope failure between min and max normal stress.
+
+        Args:
+            stressNormalMax (float): Maximum normal stress (Pa)
+            stressNormalMin (float | None, optional): Minimum normal stress.
+                If it is None, the envelop is computed from the its intersection
+                with the abscissa axis.
+
+                Defaults to None.
+            n (int, optional): Number of points to define the envelop.
+            Defaults to 100.
+
+        Returns:
+            tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: failure
+            envelop coordinates where first element is the abscissa and second
+            element is the ordinates.
+        """
+        sigmaMin: float = (
+            self.m_sigmaMin if stressNormalMin is None else stressNormalMin
+        )
+        stressNormal: npt.NDArray[np.float64] = np.linspace(
+            sigmaMin, stressNormalMax, n
+        )
+        return (stressNormal, np.array(self.computeShearStress(stressNormal)))
\ No newline at end of file
diff --git a/geos-posp/src/geos_posp/processing/geomechanicsCalculatorFunctions.py b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
similarity index 97%
rename from geos-posp/src/geos_posp/processing/geomechanicsCalculatorFunctions.py
rename to geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
index e50cadd8..568174b9 100644
--- a/geos-posp/src/geos_posp/processing/geomechanicsCalculatorFunctions.py
+++ b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
@@ -1,1011 +1,1011 @@
-# SPDX-License-Identifier: Apache-2.0
-# # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Alexandre Benedicto, Martin Lemay
-import numpy as np
-import numpy.typing as npt
-
-from geos_posp.processing.MohrCoulomb import MohrCoulomb
-from geos_posp.utils.geosUtils import getAttributeMatrixFromVector
-from geos_posp.utils.PhysicalConstants import (
-    EPSILON,
-)
-
-__doc__ = """Functions to compute additional geomechanical properties."""
-
-
-def specificGravity(
-    density: npt.NDArray[np.float64], specificDensity: float
-) -> npt.NDArray[np.float64]:
-    r"""Compute the specific gravity.
-
-    .. math::
-        SG = \frac{\rho}{\rho_f}
-
-    Args:
-        density (npt.NDArray[np.float64]): density (:math:`\rho` - kg/m³)
-        specificDensity (float): fluid density (:math:`\rho_f` - kg/m³)
-
-    Returns:
-        npt.NDArray[np.float64]: specific gravity (*SG* - no unit)
-
-    """
-    assert density is not None, "Density data must be defined"
-
-    if abs(specificDensity) < EPSILON:
-        return np.full_like(density, np.nan)
-    return density / specificDensity
-
-
-# https://en.wikipedia.org/wiki/Elastic_modulus
-def youngModulus(
-    bulkModulus: npt.NDArray[np.float64], shearModulus: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute Young modulus.
-
-    .. math::
-        E = \frac{9K.G}{3K+G}
-
-    Args:
-        bulkModulus (npt.NDArray[np.float64]): Bulk modulus (*K* - Pa)
-        shearModulus (npt.NDArray[np.float64]): Shear modulus (*G* - Pa)
-
-    Returns:
-        npt.NDArray[np.float64]: Young modulus (*E* - Pa)
-
-    """
-    assert bulkModulus is not None, "Bulk modulus must be defined"
-    assert shearModulus is not None, "Shear modulus must be defined"
-    # manage division by 0 by replacing with nan
-    assert bulkModulus.size == shearModulus.size, (
-        "Bulk modulus array and Shear modulus array"
-        + " sizes (i.e., number of cells) must be equal."
-    )
-
-    den: npt.NDArray[np.float64] = 3.0 * bulkModulus + shearModulus
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    young: npt.NDArray[np.float64] = 9.0 * bulkModulus * shearModulus / den
-    young[mask] = np.nan
-    return young
-
-
-def poissonRatio(
-    bulkModulus: npt.NDArray[np.float64], shearModulus: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute Poisson's ratio.
-
-    .. math::
-        \nu = \frac{3K-2G}{2(3K+G)}
-
-    Args:
-        bulkModulus (npt.NDArray[np.float64]): Bulk modulus (*K* - Pa)
-        shearModulus (npt.NDArray[np.float64]): Shear modulus (*G* - Pa)
-
-    Returns:
-        npt.NDArray[np.float64]: Poisson's ratio (:math:`\nu`)
-
-    """
-    assert bulkModulus is not None, "Bulk modulus must be defined"
-    assert shearModulus is not None, "Shear modulus must be defined"
-    assert bulkModulus.size == shearModulus.size, (
-        "Bulk modulus array and Shear modulus array"
-        + " sizes (i.e., number of cells) must be equal."
-    )
-    # manage division by 0 by replacing with nan
-    den: npt.NDArray[np.float64] = 2.0 * (3.0 * bulkModulus + shearModulus)
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    poisson: npt.NDArray[np.float64] = (3.0 * bulkModulus - 2.0 * shearModulus) / den
-    poisson[mask] = np.nan
-    return poisson
-
-
-def bulkModulus(
-    youngModulus: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute bulk Modulus from young modulus and poisson ratio.
-
-    .. math::
-        K = \frac{E}{3(1-2\nu)}
-
-    Args:
-        youngModulus (npt.NDArray[np.float64]): Young modulus (*E* - Pa)
-        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
-
-    Returns:
-        npt.NDArray[np.float64]: Bulk modulus (*K* - Pa)
-
-    """
-    assert poissonRatio is not None, "Poisson's ratio must be defined"
-    assert youngModulus is not None, "Young modulus must be defined"
-
-    den: npt.NDArray[np.float64] = 3.0 * (1.0 - 2.0 * poissonRatio)
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    bulkMod: npt.NDArray[np.float64] = youngModulus / den
-    bulkMod[mask] = np.nan
-    return bulkMod
-
-
-def shearModulus(
-    youngModulus: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute shear Modulus from young modulus and poisson ratio.
-
-    .. math::
-        G = \frac{E}{2(1+\nu)}
-
-    Args:
-        youngModulus (npt.NDArray[np.float64]): Young modulus (*E* - Pa)
-        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
-
-    Returns:
-        npt.NDArray[np.float64]: Shear modulus (*G* - Pa)
-
-    """
-    assert poissonRatio is not None, "Poisson's ratio must be defined"
-    assert youngModulus is not None, "Young modulus must be defined"
-
-    den: npt.NDArray[np.float64] = 2.0 * (1.0 + poissonRatio)
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    shearMod: npt.NDArray[np.float64] = youngModulus / den
-    shearMod[mask] = np.nan
-    return shearMod
-
-
-def lambdaCoefficient(
-    youngModulus: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute lambda coefficient from young modulus and Poisson ratio.
-
-    .. math::
-        \lambda = \frac{E*\nu}{(1+\nu)(1-2\nu)}
-
-    Args:
-        youngModulus (npt.NDArray[np.float64]): Young modulus (*E* - Pa)
-        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
-
-    Returns:
-        npt.NDArray[np.float64]: lambda coefficient (:math:`\lambda`)
-
-    """
-    lambdaCoeff: npt.NDArray[np.float64] = poissonRatio * youngModulus
-    den: npt.NDArray[np.float64] = (1.0 + poissonRatio) * (1.0 - 2.0 * poissonRatio)
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    lambdaCoeff /= den
-    lambdaCoeff[mask] = np.nan
-    return lambdaCoeff
-
-
-def oedometricModulus(
-    Edef: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute Oedometric modulus.
-
-    .. math::
-        M_{oed} = \frac{E_{def}}{1-2\frac{\nu^2}{1-\nu}}
-
-    Args:
-        Edef (npt.NDArray[np.float64]): Deformation modulus (:math:`E_{def}` - Pa)
-        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
-
-    Returns:
-        npt.NDArray[np.float64]: Oedometric modulus (:math:`M_{oed}` - Pa)
-
-    """
-    assert poissonRatio is not None, "Poisson's ratio must be defined"
-    assert Edef is not None, "Deformation modulus must be defined"
-
-    assert Edef.size == poissonRatio.size, (
-        "Deformation modulus array and Poisson's"
-        + " ratio array sizes (i.e., number of cells) must be equal."
-    )
-    den: npt.NDArray[np.float64] = 1.0 - (2.0 * poissonRatio * poissonRatio) / (
-        1.0 - poissonRatio
-    )
-
-    # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    EodMod: npt.NDArray[np.float64] = Edef / den
-    EodMod[mask] = np.nan
-    return EodMod
-
-
-def biotCoefficient(
-    Kgrain: float, bulkModulus: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute Biot coefficient.
-
-    .. math::
-        b = 1-\frac{K}{K_{grain}}
-
-    Args:
-        Kgrain (float): grain bulk modulus (:math:`K_{grain}` - Pa)
-        bulkModulus (npt.NDArray[np.float64]): default bulk modulus (*K* - Pa)
-
-    Returns:
-        npt.NDArray[np.float64]: biot coefficient (*b*)
-
-    """
-    assert bulkModulus is not None, "Bulk modulus must be defined"
-
-    # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(Kgrain) < EPSILON
-    den: npt.NDArray[np.float64] = np.copy(Kgrain)
-    den[mask] = 1.0
-    biot: npt.NDArray[np.float64] = 1.0 - bulkModulus / den
-    biot[mask] = np.nan
-    return biot
-
-
-def totalStress(
-    effectiveStress: npt.NDArray[np.float64],
-    biot: npt.NDArray[np.float64],
-    pressure: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
-    r"""Compute total stress from effective stress, pressure, and Biot coeff.
-
-    .. math::
-        \sigma_{tot} = \sigma_{eff}-bP
-
-    Args:
-        effectiveStress (npt.NDArray[np.float64]): effective stress
-            (:math:`\sigma_{eff}` - Pa) using Geos convention
-        biot (npt.NDArray[np.float64]): Biot coefficient (*b*)
-        pressure (npt.NDArray[np.float64]): Pore pressure (*P* - Pa)
-
-    Returns:
-        npt.NDArray[np.float64]: total stress (:math:`\sigma_{tot}` - Pa)
-
-    """
-    assert effectiveStress is not None, "Effective stress must be defined"
-    assert biot is not None, "Biot coefficient must be defined"
-    assert pressure is not None, "Pressure must be defined"
-
-    assert effectiveStress.shape[0] == biot.size, (
-        "Biot coefficient array and "
-        + "effective stress sizes (i.e., number of cells) must be equal."
-    )
-    assert biot.size == pressure.size, (
-        "Biot coefficient array and pressure array"
-        + "sizes (i.e., number of cells) must be equal."
-    )
-
-    totalStress: npt.NDArray[np.float64] = np.copy(effectiveStress)
-    # pore pressure has an effect on normal stresses only
-    # (cf. https://dnicolasespinoza.github.io/node5.html)
-    nb: int = totalStress.shape[1] if totalStress.shape[1] < 4 else 3
-    for j in range(nb):
-        totalStress[:, j] = effectiveStress[:, j] - biot * pressure
-    return totalStress
-
-
-def stressRatio(
-    horizontalStress: npt.NDArray[np.float64], verticalStress: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute horizontal to vertical stress ratio.
-
-    .. math::
-        r = \frac{\sigma_h}{\sigma_v}
-
-    Args:
-        horizontalStress (npt.NDArray[np.float64]): horizontal stress
-            (:math:`\sigma_h` - Pa)
-        verticalStress (npt.NDArray[np.float64]): vertical stress
-            (:math:`\sigma_v` - Pa)
-
-    Returns:
-        npt.NDArray[np.float64]: stress ratio (:math:`\sigma` - Pa)
-
-    """
-    assert horizontalStress is not None, "Horizontal stress must be defined"
-    assert verticalStress is not None, "Vertical stress must be defined"
-
-    assert horizontalStress.size == verticalStress.size, (
-        "Horizontal stress array "
-        + "and vertical stress array sizes (i.e., number of cells) must be equal."
-    )
-
-    # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(verticalStress) < EPSILON
-    verticalStress2: npt.NDArray[np.float64] = np.copy(verticalStress)
-    verticalStress2[mask] = 1.0
-    ratio: npt.NDArray[np.float64] = horizontalStress / verticalStress2
-    ratio[mask] = np.nan
-    return ratio
-
-
-def lithostaticStress(
-    depth: npt.NDArray[np.float64], density: npt.NDArray[np.float64], gravity: float
-) -> npt.NDArray[np.float64]:
-    """Compute the lithostatic stress.
-
-    Args:
-        depth (npt.NDArray[np.float64]): depth from surface - m)
-        density (npt.NDArray[np.float64]): density of the overburden (kg/m³)
-        gravity (float): gravity (m²/s)
-
-    Returns:
-        npt.NDArray[np.float64]: lithostatic stress (Pa)
-
-    """
-    # compute average density of the overburden of each point (replacing nan value by 0)
-
-    # TODO: the formula should not depends on the density of the cell but the average
-    # density of the overburden.
-    # But how to compute the average density of the overburden in an unstructured mesh?
-
-    # density2 = np.copy(density)
-    # density2[np.isnan(density)] = 0
-    # overburdenMeanDensity = np.cumsum(density) / np.arange(1, density.size+1, 1)
-
-    # TODO: which convention? + or -?
-
-    # TODO: Warning z coordinate may be + or - if 0 is sea level. Need to take dz from surface.
-    # Is the surface always on top of the model?
-    assert depth is not None, "Depth must be defined"
-    assert density is not None, "Density must be defined"
-
-    assert depth.size == density.size, (
-        "Depth array "
-        + "and density array sizes (i.e., number of cells) must be equal."
-    )
-    # use -1*depth to agree with Geos convention (i.e., compression with negative stress)
-    return -depth * density * gravity
-
-
-def elasticStrainFromBulkShear(
-    deltaEffectiveStress: npt.NDArray[np.float64],
-    bulkModulus: npt.NDArray[np.float64],
-    shearModulus: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
-    r"""Compute elastic strain from Bulk and Shear moduli.
-
-    See documentation on https://dnicolasespinoza.github.io/node5.html.
-
-    .. math::
-        \epsilon=\Delta\sigma_{eff}.C^{-1}
-
-
-        C=\begin{pmatrix}
-          K+\frac{4}{3}G & K-\frac{2}{3}G & K-\frac{2}{3}G &   0 &   0 &   0\\
-          K-\frac{2}{3}G & K+\frac{4}{3}G & K-\frac{2}{3}G &   0 &   0 &   0\\
-          K-\frac{2}{3}G & K-\frac{2}{3}G & K+\frac{4}{3}G &   0 &   0 &   0\\
-                       0 &              0 &              0 & \nu &   0 &   0\\
-                       0 &              0 &              0 &   0 & \nu &   0\\
-                       0 &              0 &              0 &   0 &   0 & \nu\\
-          \end{pmatrix}
-
-    where *C* is stiffness tensor.
-
-
-    Args:
-        deltaEffectiveStress (npt.NDArray[np.float64]): effective stress
-            variation (:math:`\Delta\sigma_{eff}` - Pa) [S11, S22, S33, S23, S13, S12]
-        bulkModulus (npt.NDArray[np.float64]): Bulk modulus (*K* - Pa)
-        shearModulus (npt.NDArray[np.float64]): Shear modulus (*G* - Pa)
-
-    Returns:
-        npt.NDArray[np.float64]: elastic strain (:math:`\epsilon`)
-
-    """
-    assert (
-        deltaEffectiveStress is not None
-    ), "Effective stress variation must be defined"
-    assert bulkModulus is not None, "Bulk modulus must be defined"
-    assert shearModulus is not None, "Shear modulus must be defined"
-
-    assert deltaEffectiveStress.shape[0] == bulkModulus.size, (
-        "Effective stress variation "
-        + " and bulk modulus array sizes (i.e., number of cells) must be equal."
-    )
-    assert shearModulus.size == bulkModulus.size, (
-        "Shear modulus "
-        + "and bulk modulus array sizes (i.e., number of cells) must be equal."
-    )
-    assert deltaEffectiveStress.shape[1] == 6, (
-        "Effective stress variation " + "number of components must be equal to 6."
-    )
-
-    elasticStrain: npt.NDArray[np.float64] = np.full_like(deltaEffectiveStress, np.nan)
-    for i, stressVector in enumerate(deltaEffectiveStress):
-        stiffnessTensor: npt.NDArray[np.float64] = shearModulus[i] * np.identity(
-            6, dtype=float
-        )
-        for k in range(3):
-            for m in range(3):
-                val: float = (
-                    (bulkModulus[i] + 4.0 / 3.0 * shearModulus[i])
-                    if k == m
-                    else (bulkModulus[i] - 2.0 / 3.0 * shearModulus[i])
-                )
-                stiffnessTensor[k, m] = val
-        elasticStrain[i] = stressVector @ np.linalg.inv(stiffnessTensor)
-    return elasticStrain
-
-
-def elasticStrainFromYoungPoisson(
-    deltaEffectiveStress: npt.NDArray[np.float64],
-    youngModulus: npt.NDArray[np.float64],
-    poissonRatio: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
-    r"""Compute elastic strain from Young modulus and Poisson ratio.
-
-    See documentation on https://dnicolasespinoza.github.io/node5.html.
-
-    .. math::
-        \epsilon=\Delta\sigma_{eff}.C^{-1}
-
-
-        C=\begin{pmatrix}
-          \lambda+2G &    \lambda &    \lambda &   0 &   0 &   0\\
-             \lambda & \lambda+2G &    \lambda &   0 &   0 &   0\\
-             \lambda &    \lambda & \lambda+2G &   0 &   0 &   0\\
-                   0 &          0 &          0 & \nu &   0 &   0\\
-                   0 &          0 &          0 &   0 & \nu &   0\\
-                   0 &          0 &          0 &   0 &   0 & \nu\\
-          \end{pmatrix}
-
-    where *C* is stiffness tensor, :math:`\nu` is shear modulus,
-    :math:`\lambda` is lambda coefficient.
-
-    Args:
-        deltaEffectiveStress (npt.NDArray[np.float64]): effective stress
-            variation (:math:`\Delta\sigma_{eff}` - Pa) [S11, S22, S33, S23, S13, S12]
-        youngModulus (npt.NDArray[np.float64]): Young modulus (*E* - Pa)
-        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
-
-    Returns:
-        npt.NDArray[np.float64]: elastic strain (:math:`\epsilon`)
-
-    """
-    assert (
-        deltaEffectiveStress is not None
-    ), "Effective stress variation must be defined"
-    assert youngModulus is not None, "Bulk modulus must be defined"
-    assert poissonRatio is not None, "Shear modulus must be defined"
-
-    assert deltaEffectiveStress.shape[0] == youngModulus.size, (
-        "Effective stress variation "
-        + " and bulk modulus array sizes (i.e., number of cells) must be equal."
-    )
-    assert poissonRatio.size == youngModulus.size, (
-        "Shear modulus "
-        + "and bulk modulus array sizes (i.e., number of cells) must be equal."
-    )
-    assert deltaEffectiveStress.shape[1] == 6, (
-        "Effective stress variation " + "number of components must be equal to 6."
-    )
-
-    # use of Lamé's equations
-    lambdaCoeff: npt.NDArray[np.float64] = lambdaCoefficient(youngModulus, poissonRatio)
-    shearMod: npt.NDArray[np.float64] = shearModulus(youngModulus, poissonRatio)
-
-    elasticStrain: npt.NDArray[np.float64] = np.full_like(deltaEffectiveStress, np.nan)
-    for i, deltaStressVector in enumerate(deltaEffectiveStress):
-        stiffnessTensor: npt.NDArray[np.float64] = shearMod[i] * np.identity(
-            6, dtype=float
-        )
-        for k in range(3):
-            for m in range(3):
-                val: float = (
-                    (lambdaCoeff[i] + 2.0 * shearMod[i]) if k == m else (lambdaCoeff[i])
-                )
-                stiffnessTensor[k, m] = val
-
-        elasticStrain[i] = deltaStressVector @ np.linalg.inv(stiffnessTensor)
-    return elasticStrain
-
-
-def deviatoricStressPathOed(
-    poissonRatio: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
-    r"""Compute the Deviatoric Stress Path parameter in oedometric conditions.
-
-    This parameter corresponds to the ratio between horizontal and vertical
-    stresses in oedometric conditions.
-
-    .. math::
-        DSP=\frac{\nu}{1-\nu}
-
-    Args:
-        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
-
-    Returns:
-        npt.NDArray[np.float64]: Deviatoric Stress Path parameter in
-            oedometric conditions (*DSP*)
-
-    """
-    assert poissonRatio is not None, "Poisson's ratio must be defined"
-
-    # manage division by 0 by replacing with nan
-    den: npt.NDArray[np.float64] = 1 - poissonRatio
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    ratio: npt.NDArray[np.float64] = poissonRatio / den
-    ratio[mask] = np.nan
-    return ratio
-
-
-def reservoirStressPathReal(
-    deltaStress: npt.NDArray[np.float64], deltaPressure: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute real reservoir stress path.
-
-    .. math::
-        RSP_{real}=\frac{\Delta\sigma}{\Delta P}
-
-    Args:
-        deltaStress (npt.NDArray[np.float64]): stress difference from start
-            (:math:`\Delta\sigma` - Pa)
-        deltaPressure (npt.NDArray[np.float64]): pressure difference from start
-            (:math:`\Delta P` - Pa)
-
-    Returns:
-        npt.NDArray[np.float64]: reservoir stress path (:math:`RSP_{real}`)
-
-    """
-    assert deltaPressure is not None, "Pressure deviation must be defined"
-    assert deltaStress is not None, "Stress deviation must be defined"
-
-    assert deltaStress.shape[0] == deltaPressure.size, (
-        "Total stress array and pressure variation "
-        + "array sizes (i.e., number of cells) must be equal."
-    )
-
-    # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(deltaPressure) < EPSILON
-    den: npt.NDArray[np.float64] = np.copy(deltaPressure)
-    den[mask] = 1.0
-    # use -1 to agree with Geos convention (i.e., compression with negative stress)
-    # take the xx, yy, and zz components only
-    rsp: npt.NDArray[np.float64] = np.copy(deltaStress[:, :3])
-    for j in range(rsp.shape[1]):
-        rsp[:, j] /= den
-        rsp[mask, j] = np.nan
-    return rsp
-
-
-def reservoirStressPathOed(
-    biotCoefficient: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute reservoir stress path in oedometric conditions.
-
-    .. math::
-        RSP_{oed}=b\frac{1-2\nu}{1-\nu}
-
-    Args:
-        biotCoefficient (npt.NDArray[np.float64]): biot coefficient (*b*)
-        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
-
-    Returns:
-        npt.NDArray[np.float64]: reservoir stress path (:math:`RSP_{oed}`)
-
-    """
-    assert biotCoefficient is not None, "Biot coefficient must be defined"
-    assert poissonRatio is not None, "Poisson's ratio must be defined"
-
-    assert biotCoefficient.size == poissonRatio.size, (
-        "Biot coefficient array and "
-        + "Poisson's ratio array sizes (i.e., number of cells) must be equal."
-    )
-
-    # manage division by 0 by replacing with nan
-    den: npt.NDArray[np.float64] = 1.0 - poissonRatio
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    rsp: npt.NDArray[np.float64] = biotCoefficient * (1.0 - 2.0 * poissonRatio) / den
-    rsp[mask] = np.nan
-    return rsp
-
-
-def criticalTotalStressRatio(
-    pressure: npt.NDArray[np.float64], verticalStress: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute critical total stress ratio.
-
-    Corresponds to the fracture index from Lemgruber-Traby et al (2024).
-    Fracturing can occur in areas where K > Total stress ratio.
-    (see Lemgruber-Traby, A., Cacas, M. C., Bonte, D., Rudkiewicz, J. L., Gout, C.,
-    & Cornu, T. (2024). Basin modelling workflow applied to the screening of deep
-    aquifers for potential CO2 storage. Geoenergy, geoenergy2024-010.
-    https://doi.org/10.1144/geoenergy2024-010)
-
-    .. math::
-        \sigma_{Ĉr}=\frac{P}{\sigma_v}
-
-    Args:
-        pressure (npt.NDArray[np.float64]): Pressure (*P* - Pa)
-        verticalStress (npt.NDArray[np.float64]): Vertical total stress
-            (:math:`\sigma_v` - Pa) using Geos convention
-
-    Returns:
-        npt.NDArray[np.float64]: critical total stress ratio
-            (:math:`\sigma_{Cr}`)
-
-    """
-    assert pressure is not None, "Pressure must be defined"
-    assert verticalStress is not None, "Vertical stress must be defined"
-
-    assert pressure.size == verticalStress.size, (
-        "pressure array and "
-        + "vertical stress array sizes (i.e., number of cells) must be equal."
-    )
-    # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(verticalStress) < EPSILON
-    # use -1 to agree with Geos convention (i.e., compression with negative stress)
-    verticalStress2: npt.NDArray[np.float64] = -1.0 * np.copy(verticalStress)
-    verticalStress2[mask] = 1.0
-    fi: npt.NDArray[np.float64] = pressure / verticalStress2
-    fi[mask] = np.nan
-    return fi
-
-
-def totalStressRatioThreshold(
-    pressure: npt.NDArray[np.float64], horizontalStress: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute total stress ratio threshold.
-
-    Corresponds to the fracture threshold from Lemgruber-Traby et al (2024).
-    Fracturing can occur in areas where FT > 1.
-    Equals FractureIndex / totalStressRatio.
-    (see Lemgruber-Traby, A., Cacas, M. C., Bonte, D., Rudkiewicz, J. L., Gout, C.,
-    & Cornu, T. (2024). Basin modelling workflow applied to the screening of deep
-    aquifers for potential CO2 storage. Geoenergy, geoenergy2024-010.
-    https://doi.org/10.1144/geoenergy2024-010)
-
-    .. math::
-        \sigma_{Th}=\frac{P}{\sigma_h}
-
-    Args:
-        pressure (npt.NDArray[np.float64]): Pressure (*P* - Pa)
-        horizontalStress (npt.NDArray[np.float64]): minimal horizontal total
-            stress (:math:`\sigma_h` - Pa) using Geos convention
-
-    Returns:
-        npt.NDArray[np.float64]: fracture threshold (:math:`\sigma_{Th}`)
-
-    """
-    assert pressure is not None, "Pressure must be defined"
-    assert horizontalStress is not None, "Horizontal stress must be defined"
-
-    assert pressure.size == horizontalStress.size, (
-        "pressure array and "
-        + "horizontal stress array sizes (i.e., number of cells) must be equal."
-    )
-
-    # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(horizontalStress) < EPSILON
-    # use -1 to agree with Geos convention (i.e., compression with negative stress)
-    horizontalStress2: npt.NDArray[np.float64] = -1.0 * np.copy(horizontalStress)
-    horizontalStress2[mask] = 1.0
-    ft: npt.NDArray[np.float64] = pressure / horizontalStress2
-    ft[mask] = np.nan
-    return ft
-
-
-def criticalPorePressure(
-    stressVector: npt.NDArray[np.float64],
-    rockCohesion: float,
-    frictionAngle: float = 0.0,
-) -> npt.NDArray[np.float64]:
-    r"""Compute the critical pore pressure.
-
-    Fracturing can occur in areas where Critical pore pressure is greater than
-    the pressure.
-    (see Khan, S., Khulief, Y., Juanes, R., Bashmal, S., Usman, M., & Al-Shuhail, A.
-    (2024). Geomechanical Modeling of CO2 Sequestration: A Review Focused on CO2
-    Injection and Monitoring. Journal of Environmental Chemical Engineering, 112847.
-    https://doi.org/10.1016/j.jece.2024.112847)
-
-    .. math::
-        P_{Cr}=\frac{c.\cos(\alpha)}{1-\sin(\alpha)} + \frac{3\sigma_3-\sigma_1}{2}
-
-    Args:
-        stressVector (npt.NDArray[np.float64]): stress vector
-            (:math:`\sigma` - Pa).
-        rockCohesion (float): rock cohesion (*c* - Pa)
-        frictionAngle (float, optional): friction angle (:math:`\alpha` - rad).
-
-            Defaults to 0 rad.
-
-    Returns:
-        npt.NDArray[np.float64]: critical pore pressure (:math:`P_{Cr}` - Pa)
-
-    """
-    assert stressVector is not None, "Stress vector must be defined"
-    assert stressVector.shape[1] == 6, "Stress vector must be of size 6."
-
-    assert (frictionAngle >= 0.0) and (frictionAngle < np.pi / 2.0), (
-        "Fristion angle " + "must range between 0 and pi/2."
-    )
-
-    minimumPrincipalStress: npt.NDArray[np.float64] = np.full(
-        stressVector.shape[0], np.nan
-    )
-    maximumPrincipalStress: npt.NDArray[np.float64] = np.copy(minimumPrincipalStress)
-    for i in range(minimumPrincipalStress.shape[0]):
-        p3, p2, p1 = computeStressPrincipalComponentsFromStressVector(stressVector[i])
-        minimumPrincipalStress[i] = p3
-        maximumPrincipalStress[i] = p1
-
-    # assertion frictionAngle < np.pi/2., so sin(frictionAngle) != 1
-    cohesiveTerm: npt.NDArray[np.float64] = (
-        rockCohesion * np.cos(frictionAngle) / (1 - np.sin(frictionAngle))
-    )
-    residualTerm: npt.NDArray[np.float64] = (
-        3 * minimumPrincipalStress - maximumPrincipalStress
-    ) / 2.0
-    return cohesiveTerm + residualTerm
-
-
-def criticalPorePressureThreshold(
-    pressure: npt.NDArray[np.float64], criticalPorePressure: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    r"""Compute the critical pore pressure threshold.
-
-    Defined as the ratio between pressure and critical pore pressure.
-    Fracturing can occur in areas where critical pore pressure threshold is >1.
-
-    .. math::
-        P_{Th}=\frac{P}{P_{Cr}}
-
-    Args:
-        pressure (npt.NDArray[np.float64]): pressure (*P* - Pa)
-        criticalPorePressure (npt.NDArray[np.float64]): Critical pore pressure
-            (:math:`P_{Cr}` - Pa)
-
-    Returns:
-        npt.NDArray[np.float64]: Critical pore pressure threshold (:math:`P_{Th}`)
-
-    """
-    assert pressure is not None, "Pressure must be defined"
-    assert criticalPorePressure is not None, "Critical pore pressure must be defined"
-
-    assert pressure.size == criticalPorePressure.size, (
-        "Pressure array and critical"
-        + " pore pressure array sizes (i.e., number of cells) must be equal."
-    )
-
-    # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(criticalPorePressure) < EPSILON
-    den = np.copy(criticalPorePressure)
-    den[mask] = 1.0
-    index: npt.NDArray[np.float64] = pressure / den
-    index[mask] = np.nan
-    return index
-
-
-def compressibilityOed(
-    shearModulus: npt.NDArray[np.float64],
-    bulkModulus: npt.NDArray[np.float64],
-    porosity: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
-    r"""Compute compressibility from elastic moduli and porosity.
-
-    Compressibility formula is:
-
-    .. math::
-        C = \frac{1}{\phi}.\frac{3}{3K+4G}
-
-    Args:
-        shearModulus (npt.NDArray[np.float64]): shear modulus (*G* - Pa).
-        bulkModulus (npt.NDArray[np.float64]): Bulk Modulus (*K* - Pa).
-        porosity (npt.NDArray[np.float64]): Rock porosity (:math:`\phi`).
-
-    Returns:
-        npt.NDArray[np.float64]: Oedometric Compressibility (*C* - Pa^-1).
-
-    """
-    assert bulkModulus is not None, "Bulk modulus must be defined"
-    assert shearModulus is not None, "Shear modulus must be defined"
-    assert porosity is not None, "Porosity must be defined"
-    mask1: npt.NDArray[np.bool_] = np.abs(porosity) < EPSILON
-    den1: npt.NDArray[np.float64] = np.copy(porosity)
-    den1[mask1] = 1.0
-
-    den2: npt.NDArray[np.float64] = 3.0 * bulkModulus + 4.0 * shearModulus
-    mask2: npt.NDArray[np.bool_] = np.abs(den2) < EPSILON
-    den2[mask2] = 1.0
-
-    comprOed: npt.NDArray[np.float64] = 1.0 / den1 * 3.0 / den2
-    comprOed[mask1] = np.nan
-    comprOed[mask2] = np.nan
-    return comprOed
-
-
-def compressibilityReal(
-    deltaPressure: npt.NDArray[np.float64],
-    porosity: npt.NDArray[np.float64],
-    porosityInitial: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
-    r"""Compute compressibility from elastic moduli and porosity.
-
-    Compressibility formula is:
-
-    .. math::
-        C = \frac{\phi-\phi_0}{\Delta P.\phi_0}
-
-    Args:
-        deltaPressure (npt.NDArray[np.float64]): Pressure deviation
-            (:math:`\Delta P` - Pa).
-        porosity (npt.NDArray[np.float64]): Rock porosity (:math:`\phi`).
-        porosityInitial (npt.NDArray[np.float64]): initial porosity
-            (:math:`\phi_0`).
-
-    Returns:
-        npt.NDArray[np.float64]: Real compressibility (*C* - Pa^-1).
-
-    """
-    assert deltaPressure is not None, "Pressure deviation must be defined"
-    assert porosity is not None, "Porosity must be defined"
-    assert porosityInitial is not None, "Initial porosity must be defined"
-
-    den: npt.NDArray[np.float64] = deltaPressure * porosityInitial
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-
-    comprReal: npt.NDArray[np.float64] = (porosity - porosityInitial) / den
-    comprReal[mask] = np.nan
-    return comprReal
-
-
-def compressibility(
-    poissonRatio: npt.NDArray[np.float64],
-    bulkModulus: npt.NDArray[np.float64],
-    biotCoefficient: npt.NDArray[np.float64],
-    porosity: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
-    r"""Compute compressibility from elastic moduli, biot coefficient and porosity.
-
-    Compressibility formula is:
-
-    .. math::
-        C = \frac{1-2\nu}{\phi K}\left(\frac{b²(1+\nu)}{1-\nu} + 3(b-\phi)(1-b)\right)
-
-    Args:
-        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`).
-        bulkModulus (npt.NDArray[np.float64]): Bulk Modulus (*K* - Pa)
-        biotCoefficient (npt.NDArray[np.float64]): Biot coefficient (*b*).
-        porosity (npt.NDArray[np.float64]): Rock porosity (:math:`\phi`).
-
-    Returns:
-        npt.NDArray[np.float64]: Compressibility array (*C* - Pa^-1).
-
-    """
-    assert poissonRatio is not None, "Poisson's ratio must be defined"
-    assert bulkModulus is not None, "Bulk modulus must be defined"
-    assert biotCoefficient is not None, "Biot coefficient must be defined"
-    assert porosity is not None, "Porosity must be defined"
-
-    term1: npt.NDArray[np.float64] = 1.0 - 2.0 * poissonRatio
-
-    mask: npt.NDArray[np.bool_] = (np.abs(bulkModulus) < EPSILON) * (
-        np.abs(porosity) < EPSILON
-    )
-    denFac1: npt.NDArray[np.float64] = porosity * bulkModulus
-    term1[mask] = 1.0
-    term1 /= denFac1
-
-    term2M1: npt.NDArray[np.float64] = (
-        biotCoefficient * biotCoefficient * (1 + poissonRatio)
-    )
-    denTerm2M1: npt.NDArray[np.float64] = 1 - poissonRatio
-    mask2: npt.NDArray[np.bool_] = np.abs(denTerm2M1) < EPSILON
-    denTerm2M1[mask2] = 1.0
-    term2M1 /= denTerm2M1
-    term2M1[mask2] = np.nan
-
-    term2M2: npt.NDArray[np.float64] = (
-        3.0 * (biotCoefficient - porosity) * (1 - biotCoefficient)
-    )
-    term2: npt.NDArray[np.float64] = term2M1 + term2M2
-    return term1 * term2
-
-
-def shearCapacityUtilization(
-    traction: npt.NDArray[np.float64], rockCohesion: float, frictionAngle: float
-) -> npt.NDArray[np.float64]:
-    r"""Compute shear capacity utilization (SCU).
-
-    .. math::
-        SCU = \frac{abs(\tau_1)}{\tau_{max}}
-
-    where \tau_{max} is the Mohr-Coulomb failure threshold.
-
-    Args:
-        traction (npt.NDArray[np.float64]): traction vector
-            (:math:`(\sigma, \tau_1, \tau2)` - Pa)
-        rockCohesion (float): rock cohesion (*c* - Pa).
-        frictionAngle (float): friction angle (:math:`\alpha` - rad).
-
-    Returns:
-        npt.NDArray[np.float64]: *SCU*
-
-    """
-    assert traction is not None, "Traction must be defined"
-    assert traction.shape[1] == 3, "Traction vector must have 3 components."
-
-    scu: npt.NDArray[np.float64] = np.full(traction.shape[0], np.nan)
-    for i, tractionVec in enumerate(traction):
-        # use -1 to agree with Geos convention (i.e., compression with negative stress)
-        stressNormal: npt.NDArray[np.float64] = -1.0 * tractionVec[0]
-
-        # compute failure envelope
-        mohrCoulomb: MohrCoulomb = MohrCoulomb(rockCohesion, frictionAngle)
-        tauFailure: float = float(mohrCoulomb.computeShearStress(stressNormal))
-        scu_i: float = np.nan
-        if tauFailure > 0:
-            scu_i = np.abs(tractionVec[1]) / tauFailure
-        # compute SCU
-        scu[i] = scu_i
-    return scu
-
-
-def computeStressPrincipalComponentsFromStressVector(
-    stressVector: npt.NDArray[np.float64],
-) -> tuple[float, float, float]:
-    """Compute stress principal components from stress vector.
-
-    Args:
-        stressVector (npt.NDArray[np.float64]): stress vector.
-
-    Returns:
-        tuple[float, float, float]: Principal components sorted in ascending
-            order.
-
-    """
-    assert stressVector.size == 6, "Stress vector dimension is wrong."
-    stressTensor: npt.NDArray[np.float64] = getAttributeMatrixFromVector(stressVector)
-    return computeStressPrincipalComponents(stressTensor)
-
-
-def computeStressPrincipalComponents(
-    stressTensor: npt.NDArray[np.float64],
-) -> tuple[float, float, float]:
-    """Compute stress principal components.
-
-    Args:
-        stressTensor (npt.NDArray[np.float64]): stress tensor.
-
-    Returns:
-        tuple[float, float, float]: Principal components sorted in ascending
-            order.
-
-    """
-    # get eigen values
-    e_val, e_vec = np.linalg.eig(stressTensor)
-    # sort principal stresses from smallest to largest
-    p3, p2, p1 = np.sort(e_val)
-    return (p3, p2, p1)
-
-
-def computeNormalShearStress(
-    stressTensor: npt.NDArray[np.float64], directionVector: npt.NDArray[np.float64]
-) -> tuple[float, float]:
-    """Compute normal and shear stress according to stress tensor and direction.
-
-    Args:
-        stressTensor (npt.NDArray[np.float64]): 3x3 stress tensor
-        directionVector (npt.NDArray[np.float64]): direction vector
-
-    Returns:
-        tuple[float, float]: normal and shear stresses.
-
-    """
-    assert stressTensor.shape == (3, 3), "Stress tensor must be 3x3 matrix."
-    assert directionVector.size == 3, "Direction vector must have 3 components"
-
-    # normalization of direction vector
-    directionVector = directionVector / np.linalg.norm(directionVector)
-    # stress vector
-    T: npt.NDArray[np.float64] = np.dot(stressTensor, directionVector)
-    # normal stress
-    sigmaN: float = np.dot(T, directionVector)
-    # shear stress
-    tauVec: npt.NDArray[np.float64] = T - np.dot(sigmaN, directionVector)
-    tau: float = float(np.linalg.norm(tauVec))
-    return (sigmaN, tau)
+# SPDX-License-Identifier: Apache-2.0
+# # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Alexandre Benedicto, Martin Lemay
+import numpy as np
+import numpy.typing as npt
+
+from geos-geomechanics.processing.MohrCoulomb import MohrCoulomb
+from geos_posp.utils.geosUtils import getAttributeMatrixFromVector
+from geos_posp.utils.PhysicalConstants import (
+    EPSILON,
+)
+
+__doc__ = """Functions to compute additional geomechanical properties."""
+
+
+def specificGravity(
+    density: npt.NDArray[np.float64], specificDensity: float
+) -> npt.NDArray[np.float64]:
+    r"""Compute the specific gravity.
+
+    .. math::
+        SG = \frac{\rho}{\rho_f}
+
+    Args:
+        density (npt.NDArray[np.float64]): density (:math:`\rho` - kg/m³)
+        specificDensity (float): fluid density (:math:`\rho_f` - kg/m³)
+
+    Returns:
+        npt.NDArray[np.float64]: specific gravity (*SG* - no unit)
+
+    """
+    assert density is not None, "Density data must be defined"
+
+    if abs(specificDensity) < EPSILON:
+        return np.full_like(density, np.nan)
+    return density / specificDensity
+
+
+# https://en.wikipedia.org/wiki/Elastic_modulus
+def youngModulus(
+    bulkModulus: npt.NDArray[np.float64], shearModulus: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute Young modulus.
+
+    .. math::
+        E = \frac{9K.G}{3K+G}
+
+    Args:
+        bulkModulus (npt.NDArray[np.float64]): Bulk modulus (*K* - Pa)
+        shearModulus (npt.NDArray[np.float64]): Shear modulus (*G* - Pa)
+
+    Returns:
+        npt.NDArray[np.float64]: Young modulus (*E* - Pa)
+
+    """
+    assert bulkModulus is not None, "Bulk modulus must be defined"
+    assert shearModulus is not None, "Shear modulus must be defined"
+    # manage division by 0 by replacing with nan
+    assert bulkModulus.size == shearModulus.size, (
+        "Bulk modulus array and Shear modulus array"
+        + " sizes (i.e., number of cells) must be equal."
+    )
+
+    den: npt.NDArray[np.float64] = 3.0 * bulkModulus + shearModulus
+    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
+    den[mask] = 1.0
+    young: npt.NDArray[np.float64] = 9.0 * bulkModulus * shearModulus / den
+    young[mask] = np.nan
+    return young
+
+
+def poissonRatio(
+    bulkModulus: npt.NDArray[np.float64], shearModulus: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute Poisson's ratio.
+
+    .. math::
+        \nu = \frac{3K-2G}{2(3K+G)}
+
+    Args:
+        bulkModulus (npt.NDArray[np.float64]): Bulk modulus (*K* - Pa)
+        shearModulus (npt.NDArray[np.float64]): Shear modulus (*G* - Pa)
+
+    Returns:
+        npt.NDArray[np.float64]: Poisson's ratio (:math:`\nu`)
+
+    """
+    assert bulkModulus is not None, "Bulk modulus must be defined"
+    assert shearModulus is not None, "Shear modulus must be defined"
+    assert bulkModulus.size == shearModulus.size, (
+        "Bulk modulus array and Shear modulus array"
+        + " sizes (i.e., number of cells) must be equal."
+    )
+    # manage division by 0 by replacing with nan
+    den: npt.NDArray[np.float64] = 2.0 * (3.0 * bulkModulus + shearModulus)
+    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
+    den[mask] = 1.0
+    poisson: npt.NDArray[np.float64] = (3.0 * bulkModulus - 2.0 * shearModulus) / den
+    poisson[mask] = np.nan
+    return poisson
+
+
+def bulkModulus(
+    youngModulus: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute bulk Modulus from young modulus and poisson ratio.
+
+    .. math::
+        K = \frac{E}{3(1-2\nu)}
+
+    Args:
+        youngModulus (npt.NDArray[np.float64]): Young modulus (*E* - Pa)
+        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
+
+    Returns:
+        npt.NDArray[np.float64]: Bulk modulus (*K* - Pa)
+
+    """
+    assert poissonRatio is not None, "Poisson's ratio must be defined"
+    assert youngModulus is not None, "Young modulus must be defined"
+
+    den: npt.NDArray[np.float64] = 3.0 * (1.0 - 2.0 * poissonRatio)
+    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
+    den[mask] = 1.0
+    bulkMod: npt.NDArray[np.float64] = youngModulus / den
+    bulkMod[mask] = np.nan
+    return bulkMod
+
+
+def shearModulus(
+    youngModulus: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute shear Modulus from young modulus and poisson ratio.
+
+    .. math::
+        G = \frac{E}{2(1+\nu)}
+
+    Args:
+        youngModulus (npt.NDArray[np.float64]): Young modulus (*E* - Pa)
+        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
+
+    Returns:
+        npt.NDArray[np.float64]: Shear modulus (*G* - Pa)
+
+    """
+    assert poissonRatio is not None, "Poisson's ratio must be defined"
+    assert youngModulus is not None, "Young modulus must be defined"
+
+    den: npt.NDArray[np.float64] = 2.0 * (1.0 + poissonRatio)
+    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
+    den[mask] = 1.0
+    shearMod: npt.NDArray[np.float64] = youngModulus / den
+    shearMod[mask] = np.nan
+    return shearMod
+
+
+def lambdaCoefficient(
+    youngModulus: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute lambda coefficient from young modulus and Poisson ratio.
+
+    .. math::
+        \lambda = \frac{E*\nu}{(1+\nu)(1-2\nu)}
+
+    Args:
+        youngModulus (npt.NDArray[np.float64]): Young modulus (*E* - Pa)
+        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
+
+    Returns:
+        npt.NDArray[np.float64]: lambda coefficient (:math:`\lambda`)
+
+    """
+    lambdaCoeff: npt.NDArray[np.float64] = poissonRatio * youngModulus
+    den: npt.NDArray[np.float64] = (1.0 + poissonRatio) * (1.0 - 2.0 * poissonRatio)
+    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
+    den[mask] = 1.0
+    lambdaCoeff /= den
+    lambdaCoeff[mask] = np.nan
+    return lambdaCoeff
+
+
+def oedometricModulus(
+    Edef: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute Oedometric modulus.
+
+    .. math::
+        M_{oed} = \frac{E_{def}}{1-2\frac{\nu^2}{1-\nu}}
+
+    Args:
+        Edef (npt.NDArray[np.float64]): Deformation modulus (:math:`E_{def}` - Pa)
+        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
+
+    Returns:
+        npt.NDArray[np.float64]: Oedometric modulus (:math:`M_{oed}` - Pa)
+
+    """
+    assert poissonRatio is not None, "Poisson's ratio must be defined"
+    assert Edef is not None, "Deformation modulus must be defined"
+
+    assert Edef.size == poissonRatio.size, (
+        "Deformation modulus array and Poisson's"
+        + " ratio array sizes (i.e., number of cells) must be equal."
+    )
+    den: npt.NDArray[np.float64] = 1.0 - (2.0 * poissonRatio * poissonRatio) / (
+        1.0 - poissonRatio
+    )
+
+    # manage division by 0 by replacing with nan
+    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
+    den[mask] = 1.0
+    EodMod: npt.NDArray[np.float64] = Edef / den
+    EodMod[mask] = np.nan
+    return EodMod
+
+
+def biotCoefficient(
+    Kgrain: float, bulkModulus: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute Biot coefficient.
+
+    .. math::
+        b = 1-\frac{K}{K_{grain}}
+
+    Args:
+        Kgrain (float): grain bulk modulus (:math:`K_{grain}` - Pa)
+        bulkModulus (npt.NDArray[np.float64]): default bulk modulus (*K* - Pa)
+
+    Returns:
+        npt.NDArray[np.float64]: biot coefficient (*b*)
+
+    """
+    assert bulkModulus is not None, "Bulk modulus must be defined"
+
+    # manage division by 0 by replacing with nan
+    mask: npt.NDArray[np.bool_] = np.abs(Kgrain) < EPSILON
+    den: npt.NDArray[np.float64] = np.copy(Kgrain)
+    den[mask] = 1.0
+    biot: npt.NDArray[np.float64] = 1.0 - bulkModulus / den
+    biot[mask] = np.nan
+    return biot
+
+
+def totalStress(
+    effectiveStress: npt.NDArray[np.float64],
+    biot: npt.NDArray[np.float64],
+    pressure: npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
+    r"""Compute total stress from effective stress, pressure, and Biot coeff.
+
+    .. math::
+        \sigma_{tot} = \sigma_{eff}-bP
+
+    Args:
+        effectiveStress (npt.NDArray[np.float64]): effective stress
+            (:math:`\sigma_{eff}` - Pa) using Geos convention
+        biot (npt.NDArray[np.float64]): Biot coefficient (*b*)
+        pressure (npt.NDArray[np.float64]): Pore pressure (*P* - Pa)
+
+    Returns:
+        npt.NDArray[np.float64]: total stress (:math:`\sigma_{tot}` - Pa)
+
+    """
+    assert effectiveStress is not None, "Effective stress must be defined"
+    assert biot is not None, "Biot coefficient must be defined"
+    assert pressure is not None, "Pressure must be defined"
+
+    assert effectiveStress.shape[0] == biot.size, (
+        "Biot coefficient array and "
+        + "effective stress sizes (i.e., number of cells) must be equal."
+    )
+    assert biot.size == pressure.size, (
+        "Biot coefficient array and pressure array"
+        + "sizes (i.e., number of cells) must be equal."
+    )
+
+    totalStress: npt.NDArray[np.float64] = np.copy(effectiveStress)
+    # pore pressure has an effect on normal stresses only
+    # (cf. https://dnicolasespinoza.github.io/node5.html)
+    nb: int = totalStress.shape[1] if totalStress.shape[1] < 4 else 3
+    for j in range(nb):
+        totalStress[:, j] = effectiveStress[:, j] - biot * pressure
+    return totalStress
+
+
+def stressRatio(
+    horizontalStress: npt.NDArray[np.float64], verticalStress: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute horizontal to vertical stress ratio.
+
+    .. math::
+        r = \frac{\sigma_h}{\sigma_v}
+
+    Args:
+        horizontalStress (npt.NDArray[np.float64]): horizontal stress
+            (:math:`\sigma_h` - Pa)
+        verticalStress (npt.NDArray[np.float64]): vertical stress
+            (:math:`\sigma_v` - Pa)
+
+    Returns:
+        npt.NDArray[np.float64]: stress ratio (:math:`\sigma` - Pa)
+
+    """
+    assert horizontalStress is not None, "Horizontal stress must be defined"
+    assert verticalStress is not None, "Vertical stress must be defined"
+
+    assert horizontalStress.size == verticalStress.size, (
+        "Horizontal stress array "
+        + "and vertical stress array sizes (i.e., number of cells) must be equal."
+    )
+
+    # manage division by 0 by replacing with nan
+    mask: npt.NDArray[np.bool_] = np.abs(verticalStress) < EPSILON
+    verticalStress2: npt.NDArray[np.float64] = np.copy(verticalStress)
+    verticalStress2[mask] = 1.0
+    ratio: npt.NDArray[np.float64] = horizontalStress / verticalStress2
+    ratio[mask] = np.nan
+    return ratio
+
+
+def lithostaticStress(
+    depth: npt.NDArray[np.float64], density: npt.NDArray[np.float64], gravity: float
+) -> npt.NDArray[np.float64]:
+    """Compute the lithostatic stress.
+
+    Args:
+        depth (npt.NDArray[np.float64]): depth from surface - m)
+        density (npt.NDArray[np.float64]): density of the overburden (kg/m³)
+        gravity (float): gravity (m²/s)
+
+    Returns:
+        npt.NDArray[np.float64]: lithostatic stress (Pa)
+
+    """
+    # compute average density of the overburden of each point (replacing nan value by 0)
+
+    # TODO: the formula should not depends on the density of the cell but the average
+    # density of the overburden.
+    # But how to compute the average density of the overburden in an unstructured mesh?
+
+    # density2 = np.copy(density)
+    # density2[np.isnan(density)] = 0
+    # overburdenMeanDensity = np.cumsum(density) / np.arange(1, density.size+1, 1)
+
+    # TODO: which convention? + or -?
+
+    # TODO: Warning z coordinate may be + or - if 0 is sea level. Need to take dz from surface.
+    # Is the surface always on top of the model?
+    assert depth is not None, "Depth must be defined"
+    assert density is not None, "Density must be defined"
+
+    assert depth.size == density.size, (
+        "Depth array "
+        + "and density array sizes (i.e., number of cells) must be equal."
+    )
+    # use -1*depth to agree with Geos convention (i.e., compression with negative stress)
+    return -depth * density * gravity
+
+
+def elasticStrainFromBulkShear(
+    deltaEffectiveStress: npt.NDArray[np.float64],
+    bulkModulus: npt.NDArray[np.float64],
+    shearModulus: npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
+    r"""Compute elastic strain from Bulk and Shear moduli.
+
+    See documentation on https://dnicolasespinoza.github.io/node5.html.
+
+    .. math::
+        \epsilon=\Delta\sigma_{eff}.C^{-1}
+
+
+        C=\begin{pmatrix}
+          K+\frac{4}{3}G & K-\frac{2}{3}G & K-\frac{2}{3}G &   0 &   0 &   0\\
+          K-\frac{2}{3}G & K+\frac{4}{3}G & K-\frac{2}{3}G &   0 &   0 &   0\\
+          K-\frac{2}{3}G & K-\frac{2}{3}G & K+\frac{4}{3}G &   0 &   0 &   0\\
+                       0 &              0 &              0 & \nu &   0 &   0\\
+                       0 &              0 &              0 &   0 & \nu &   0\\
+                       0 &              0 &              0 &   0 &   0 & \nu\\
+          \end{pmatrix}
+
+    where *C* is stiffness tensor.
+
+
+    Args:
+        deltaEffectiveStress (npt.NDArray[np.float64]): effective stress
+            variation (:math:`\Delta\sigma_{eff}` - Pa) [S11, S22, S33, S23, S13, S12]
+        bulkModulus (npt.NDArray[np.float64]): Bulk modulus (*K* - Pa)
+        shearModulus (npt.NDArray[np.float64]): Shear modulus (*G* - Pa)
+
+    Returns:
+        npt.NDArray[np.float64]: elastic strain (:math:`\epsilon`)
+
+    """
+    assert (
+        deltaEffectiveStress is not None
+    ), "Effective stress variation must be defined"
+    assert bulkModulus is not None, "Bulk modulus must be defined"
+    assert shearModulus is not None, "Shear modulus must be defined"
+
+    assert deltaEffectiveStress.shape[0] == bulkModulus.size, (
+        "Effective stress variation "
+        + " and bulk modulus array sizes (i.e., number of cells) must be equal."
+    )
+    assert shearModulus.size == bulkModulus.size, (
+        "Shear modulus "
+        + "and bulk modulus array sizes (i.e., number of cells) must be equal."
+    )
+    assert deltaEffectiveStress.shape[1] == 6, (
+        "Effective stress variation " + "number of components must be equal to 6."
+    )
+
+    elasticStrain: npt.NDArray[np.float64] = np.full_like(deltaEffectiveStress, np.nan)
+    for i, stressVector in enumerate(deltaEffectiveStress):
+        stiffnessTensor: npt.NDArray[np.float64] = shearModulus[i] * np.identity(
+            6, dtype=float
+        )
+        for k in range(3):
+            for m in range(3):
+                val: float = (
+                    (bulkModulus[i] + 4.0 / 3.0 * shearModulus[i])
+                    if k == m
+                    else (bulkModulus[i] - 2.0 / 3.0 * shearModulus[i])
+                )
+                stiffnessTensor[k, m] = val
+        elasticStrain[i] = stressVector @ np.linalg.inv(stiffnessTensor)
+    return elasticStrain
+
+
+def elasticStrainFromYoungPoisson(
+    deltaEffectiveStress: npt.NDArray[np.float64],
+    youngModulus: npt.NDArray[np.float64],
+    poissonRatio: npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
+    r"""Compute elastic strain from Young modulus and Poisson ratio.
+
+    See documentation on https://dnicolasespinoza.github.io/node5.html.
+
+    .. math::
+        \epsilon=\Delta\sigma_{eff}.C^{-1}
+
+
+        C=\begin{pmatrix}
+          \lambda+2G &    \lambda &    \lambda &   0 &   0 &   0\\
+             \lambda & \lambda+2G &    \lambda &   0 &   0 &   0\\
+             \lambda &    \lambda & \lambda+2G &   0 &   0 &   0\\
+                   0 &          0 &          0 & \nu &   0 &   0\\
+                   0 &          0 &          0 &   0 & \nu &   0\\
+                   0 &          0 &          0 &   0 &   0 & \nu\\
+          \end{pmatrix}
+
+    where *C* is stiffness tensor, :math:`\nu` is shear modulus,
+    :math:`\lambda` is lambda coefficient.
+
+    Args:
+        deltaEffectiveStress (npt.NDArray[np.float64]): effective stress
+            variation (:math:`\Delta\sigma_{eff}` - Pa) [S11, S22, S33, S23, S13, S12]
+        youngModulus (npt.NDArray[np.float64]): Young modulus (*E* - Pa)
+        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
+
+    Returns:
+        npt.NDArray[np.float64]: elastic strain (:math:`\epsilon`)
+
+    """
+    assert (
+        deltaEffectiveStress is not None
+    ), "Effective stress variation must be defined"
+    assert youngModulus is not None, "Bulk modulus must be defined"
+    assert poissonRatio is not None, "Shear modulus must be defined"
+
+    assert deltaEffectiveStress.shape[0] == youngModulus.size, (
+        "Effective stress variation "
+        + " and bulk modulus array sizes (i.e., number of cells) must be equal."
+    )
+    assert poissonRatio.size == youngModulus.size, (
+        "Shear modulus "
+        + "and bulk modulus array sizes (i.e., number of cells) must be equal."
+    )
+    assert deltaEffectiveStress.shape[1] == 6, (
+        "Effective stress variation " + "number of components must be equal to 6."
+    )
+
+    # use of Lamé's equations
+    lambdaCoeff: npt.NDArray[np.float64] = lambdaCoefficient(youngModulus, poissonRatio)
+    shearMod: npt.NDArray[np.float64] = shearModulus(youngModulus, poissonRatio)
+
+    elasticStrain: npt.NDArray[np.float64] = np.full_like(deltaEffectiveStress, np.nan)
+    for i, deltaStressVector in enumerate(deltaEffectiveStress):
+        stiffnessTensor: npt.NDArray[np.float64] = shearMod[i] * np.identity(
+            6, dtype=float
+        )
+        for k in range(3):
+            for m in range(3):
+                val: float = (
+                    (lambdaCoeff[i] + 2.0 * shearMod[i]) if k == m else (lambdaCoeff[i])
+                )
+                stiffnessTensor[k, m] = val
+
+        elasticStrain[i] = deltaStressVector @ np.linalg.inv(stiffnessTensor)
+    return elasticStrain
+
+
+def deviatoricStressPathOed(
+    poissonRatio: npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
+    r"""Compute the Deviatoric Stress Path parameter in oedometric conditions.
+
+    This parameter corresponds to the ratio between horizontal and vertical
+    stresses in oedometric conditions.
+
+    .. math::
+        DSP=\frac{\nu}{1-\nu}
+
+    Args:
+        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
+
+    Returns:
+        npt.NDArray[np.float64]: Deviatoric Stress Path parameter in
+            oedometric conditions (*DSP*)
+
+    """
+    assert poissonRatio is not None, "Poisson's ratio must be defined"
+
+    # manage division by 0 by replacing with nan
+    den: npt.NDArray[np.float64] = 1 - poissonRatio
+    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
+    den[mask] = 1.0
+    ratio: npt.NDArray[np.float64] = poissonRatio / den
+    ratio[mask] = np.nan
+    return ratio
+
+
+def reservoirStressPathReal(
+    deltaStress: npt.NDArray[np.float64], deltaPressure: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute real reservoir stress path.
+
+    .. math::
+        RSP_{real}=\frac{\Delta\sigma}{\Delta P}
+
+    Args:
+        deltaStress (npt.NDArray[np.float64]): stress difference from start
+            (:math:`\Delta\sigma` - Pa)
+        deltaPressure (npt.NDArray[np.float64]): pressure difference from start
+            (:math:`\Delta P` - Pa)
+
+    Returns:
+        npt.NDArray[np.float64]: reservoir stress path (:math:`RSP_{real}`)
+
+    """
+    assert deltaPressure is not None, "Pressure deviation must be defined"
+    assert deltaStress is not None, "Stress deviation must be defined"
+
+    assert deltaStress.shape[0] == deltaPressure.size, (
+        "Total stress array and pressure variation "
+        + "array sizes (i.e., number of cells) must be equal."
+    )
+
+    # manage division by 0 by replacing with nan
+    mask: npt.NDArray[np.bool_] = np.abs(deltaPressure) < EPSILON
+    den: npt.NDArray[np.float64] = np.copy(deltaPressure)
+    den[mask] = 1.0
+    # use -1 to agree with Geos convention (i.e., compression with negative stress)
+    # take the xx, yy, and zz components only
+    rsp: npt.NDArray[np.float64] = np.copy(deltaStress[:, :3])
+    for j in range(rsp.shape[1]):
+        rsp[:, j] /= den
+        rsp[mask, j] = np.nan
+    return rsp
+
+
+def reservoirStressPathOed(
+    biotCoefficient: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute reservoir stress path in oedometric conditions.
+
+    .. math::
+        RSP_{oed}=b\frac{1-2\nu}{1-\nu}
+
+    Args:
+        biotCoefficient (npt.NDArray[np.float64]): biot coefficient (*b*)
+        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`)
+
+    Returns:
+        npt.NDArray[np.float64]: reservoir stress path (:math:`RSP_{oed}`)
+
+    """
+    assert biotCoefficient is not None, "Biot coefficient must be defined"
+    assert poissonRatio is not None, "Poisson's ratio must be defined"
+
+    assert biotCoefficient.size == poissonRatio.size, (
+        "Biot coefficient array and "
+        + "Poisson's ratio array sizes (i.e., number of cells) must be equal."
+    )
+
+    # manage division by 0 by replacing with nan
+    den: npt.NDArray[np.float64] = 1.0 - poissonRatio
+    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
+    den[mask] = 1.0
+    rsp: npt.NDArray[np.float64] = biotCoefficient * (1.0 - 2.0 * poissonRatio) / den
+    rsp[mask] = np.nan
+    return rsp
+
+
+def criticalTotalStressRatio(
+    pressure: npt.NDArray[np.float64], verticalStress: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute critical total stress ratio.
+
+    Corresponds to the fracture index from Lemgruber-Traby et al (2024).
+    Fracturing can occur in areas where K > Total stress ratio.
+    (see Lemgruber-Traby, A., Cacas, M. C., Bonte, D., Rudkiewicz, J. L., Gout, C.,
+    & Cornu, T. (2024). Basin modelling workflow applied to the screening of deep
+    aquifers for potential CO2 storage. Geoenergy, geoenergy2024-010.
+    https://doi.org/10.1144/geoenergy2024-010)
+
+    .. math::
+        \sigma_{Ĉr}=\frac{P}{\sigma_v}
+
+    Args:
+        pressure (npt.NDArray[np.float64]): Pressure (*P* - Pa)
+        verticalStress (npt.NDArray[np.float64]): Vertical total stress
+            (:math:`\sigma_v` - Pa) using Geos convention
+
+    Returns:
+        npt.NDArray[np.float64]: critical total stress ratio
+            (:math:`\sigma_{Cr}`)
+
+    """
+    assert pressure is not None, "Pressure must be defined"
+    assert verticalStress is not None, "Vertical stress must be defined"
+
+    assert pressure.size == verticalStress.size, (
+        "pressure array and "
+        + "vertical stress array sizes (i.e., number of cells) must be equal."
+    )
+    # manage division by 0 by replacing with nan
+    mask: npt.NDArray[np.bool_] = np.abs(verticalStress) < EPSILON
+    # use -1 to agree with Geos convention (i.e., compression with negative stress)
+    verticalStress2: npt.NDArray[np.float64] = -1.0 * np.copy(verticalStress)
+    verticalStress2[mask] = 1.0
+    fi: npt.NDArray[np.float64] = pressure / verticalStress2
+    fi[mask] = np.nan
+    return fi
+
+
+def totalStressRatioThreshold(
+    pressure: npt.NDArray[np.float64], horizontalStress: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute total stress ratio threshold.
+
+    Corresponds to the fracture threshold from Lemgruber-Traby et al (2024).
+    Fracturing can occur in areas where FT > 1.
+    Equals FractureIndex / totalStressRatio.
+    (see Lemgruber-Traby, A., Cacas, M. C., Bonte, D., Rudkiewicz, J. L., Gout, C.,
+    & Cornu, T. (2024). Basin modelling workflow applied to the screening of deep
+    aquifers for potential CO2 storage. Geoenergy, geoenergy2024-010.
+    https://doi.org/10.1144/geoenergy2024-010)
+
+    .. math::
+        \sigma_{Th}=\frac{P}{\sigma_h}
+
+    Args:
+        pressure (npt.NDArray[np.float64]): Pressure (*P* - Pa)
+        horizontalStress (npt.NDArray[np.float64]): minimal horizontal total
+            stress (:math:`\sigma_h` - Pa) using Geos convention
+
+    Returns:
+        npt.NDArray[np.float64]: fracture threshold (:math:`\sigma_{Th}`)
+
+    """
+    assert pressure is not None, "Pressure must be defined"
+    assert horizontalStress is not None, "Horizontal stress must be defined"
+
+    assert pressure.size == horizontalStress.size, (
+        "pressure array and "
+        + "horizontal stress array sizes (i.e., number of cells) must be equal."
+    )
+
+    # manage division by 0 by replacing with nan
+    mask: npt.NDArray[np.bool_] = np.abs(horizontalStress) < EPSILON
+    # use -1 to agree with Geos convention (i.e., compression with negative stress)
+    horizontalStress2: npt.NDArray[np.float64] = -1.0 * np.copy(horizontalStress)
+    horizontalStress2[mask] = 1.0
+    ft: npt.NDArray[np.float64] = pressure / horizontalStress2
+    ft[mask] = np.nan
+    return ft
+
+
+def criticalPorePressure(
+    stressVector: npt.NDArray[np.float64],
+    rockCohesion: float,
+    frictionAngle: float = 0.0,
+) -> npt.NDArray[np.float64]:
+    r"""Compute the critical pore pressure.
+
+    Fracturing can occur in areas where Critical pore pressure is greater than
+    the pressure.
+    (see Khan, S., Khulief, Y., Juanes, R., Bashmal, S., Usman, M., & Al-Shuhail, A.
+    (2024). Geomechanical Modeling of CO2 Sequestration: A Review Focused on CO2
+    Injection and Monitoring. Journal of Environmental Chemical Engineering, 112847.
+    https://doi.org/10.1016/j.jece.2024.112847)
+
+    .. math::
+        P_{Cr}=\frac{c.\cos(\alpha)}{1-\sin(\alpha)} + \frac{3\sigma_3-\sigma_1}{2}
+
+    Args:
+        stressVector (npt.NDArray[np.float64]): stress vector
+            (:math:`\sigma` - Pa).
+        rockCohesion (float): rock cohesion (*c* - Pa)
+        frictionAngle (float, optional): friction angle (:math:`\alpha` - rad).
+
+            Defaults to 0 rad.
+
+    Returns:
+        npt.NDArray[np.float64]: critical pore pressure (:math:`P_{Cr}` - Pa)
+
+    """
+    assert stressVector is not None, "Stress vector must be defined"
+    assert stressVector.shape[1] == 6, "Stress vector must be of size 6."
+
+    assert (frictionAngle >= 0.0) and (frictionAngle < np.pi / 2.0), (
+        "Fristion angle " + "must range between 0 and pi/2."
+    )
+
+    minimumPrincipalStress: npt.NDArray[np.float64] = np.full(
+        stressVector.shape[0], np.nan
+    )
+    maximumPrincipalStress: npt.NDArray[np.float64] = np.copy(minimumPrincipalStress)
+    for i in range(minimumPrincipalStress.shape[0]):
+        p3, p2, p1 = computeStressPrincipalComponentsFromStressVector(stressVector[i])
+        minimumPrincipalStress[i] = p3
+        maximumPrincipalStress[i] = p1
+
+    # assertion frictionAngle < np.pi/2., so sin(frictionAngle) != 1
+    cohesiveTerm: npt.NDArray[np.float64] = (
+        rockCohesion * np.cos(frictionAngle) / (1 - np.sin(frictionAngle))
+    )
+    residualTerm: npt.NDArray[np.float64] = (
+        3 * minimumPrincipalStress - maximumPrincipalStress
+    ) / 2.0
+    return cohesiveTerm + residualTerm
+
+
+def criticalPorePressureThreshold(
+    pressure: npt.NDArray[np.float64], criticalPorePressure: npt.NDArray[np.float64]
+) -> npt.NDArray[np.float64]:
+    r"""Compute the critical pore pressure threshold.
+
+    Defined as the ratio between pressure and critical pore pressure.
+    Fracturing can occur in areas where critical pore pressure threshold is >1.
+
+    .. math::
+        P_{Th}=\frac{P}{P_{Cr}}
+
+    Args:
+        pressure (npt.NDArray[np.float64]): pressure (*P* - Pa)
+        criticalPorePressure (npt.NDArray[np.float64]): Critical pore pressure
+            (:math:`P_{Cr}` - Pa)
+
+    Returns:
+        npt.NDArray[np.float64]: Critical pore pressure threshold (:math:`P_{Th}`)
+
+    """
+    assert pressure is not None, "Pressure must be defined"
+    assert criticalPorePressure is not None, "Critical pore pressure must be defined"
+
+    assert pressure.size == criticalPorePressure.size, (
+        "Pressure array and critical"
+        + " pore pressure array sizes (i.e., number of cells) must be equal."
+    )
+
+    # manage division by 0 by replacing with nan
+    mask: npt.NDArray[np.bool_] = np.abs(criticalPorePressure) < EPSILON
+    den = np.copy(criticalPorePressure)
+    den[mask] = 1.0
+    index: npt.NDArray[np.float64] = pressure / den
+    index[mask] = np.nan
+    return index
+
+
+def compressibilityOed(
+    shearModulus: npt.NDArray[np.float64],
+    bulkModulus: npt.NDArray[np.float64],
+    porosity: npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
+    r"""Compute compressibility from elastic moduli and porosity.
+
+    Compressibility formula is:
+
+    .. math::
+        C = \frac{1}{\phi}.\frac{3}{3K+4G}
+
+    Args:
+        shearModulus (npt.NDArray[np.float64]): shear modulus (*G* - Pa).
+        bulkModulus (npt.NDArray[np.float64]): Bulk Modulus (*K* - Pa).
+        porosity (npt.NDArray[np.float64]): Rock porosity (:math:`\phi`).
+
+    Returns:
+        npt.NDArray[np.float64]: Oedometric Compressibility (*C* - Pa^-1).
+
+    """
+    assert bulkModulus is not None, "Bulk modulus must be defined"
+    assert shearModulus is not None, "Shear modulus must be defined"
+    assert porosity is not None, "Porosity must be defined"
+    mask1: npt.NDArray[np.bool_] = np.abs(porosity) < EPSILON
+    den1: npt.NDArray[np.float64] = np.copy(porosity)
+    den1[mask1] = 1.0
+
+    den2: npt.NDArray[np.float64] = 3.0 * bulkModulus + 4.0 * shearModulus
+    mask2: npt.NDArray[np.bool_] = np.abs(den2) < EPSILON
+    den2[mask2] = 1.0
+
+    comprOed: npt.NDArray[np.float64] = 1.0 / den1 * 3.0 / den2
+    comprOed[mask1] = np.nan
+    comprOed[mask2] = np.nan
+    return comprOed
+
+
+def compressibilityReal(
+    deltaPressure: npt.NDArray[np.float64],
+    porosity: npt.NDArray[np.float64],
+    porosityInitial: npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
+    r"""Compute compressibility from elastic moduli and porosity.
+
+    Compressibility formula is:
+
+    .. math::
+        C = \frac{\phi-\phi_0}{\Delta P.\phi_0}
+
+    Args:
+        deltaPressure (npt.NDArray[np.float64]): Pressure deviation
+            (:math:`\Delta P` - Pa).
+        porosity (npt.NDArray[np.float64]): Rock porosity (:math:`\phi`).
+        porosityInitial (npt.NDArray[np.float64]): initial porosity
+            (:math:`\phi_0`).
+
+    Returns:
+        npt.NDArray[np.float64]: Real compressibility (*C* - Pa^-1).
+
+    """
+    assert deltaPressure is not None, "Pressure deviation must be defined"
+    assert porosity is not None, "Porosity must be defined"
+    assert porosityInitial is not None, "Initial porosity must be defined"
+
+    den: npt.NDArray[np.float64] = deltaPressure * porosityInitial
+    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
+    den[mask] = 1.0
+
+    comprReal: npt.NDArray[np.float64] = (porosity - porosityInitial) / den
+    comprReal[mask] = np.nan
+    return comprReal
+
+
+def compressibility(
+    poissonRatio: npt.NDArray[np.float64],
+    bulkModulus: npt.NDArray[np.float64],
+    biotCoefficient: npt.NDArray[np.float64],
+    porosity: npt.NDArray[np.float64],
+) -> npt.NDArray[np.float64]:
+    r"""Compute compressibility from elastic moduli, biot coefficient and porosity.
+
+    Compressibility formula is:
+
+    .. math::
+        C = \frac{1-2\nu}{\phi K}\left(\frac{b²(1+\nu)}{1-\nu} + 3(b-\phi)(1-b)\right)
+
+    Args:
+        poissonRatio (npt.NDArray[np.float64]): Poisson's ratio (:math:`\nu`).
+        bulkModulus (npt.NDArray[np.float64]): Bulk Modulus (*K* - Pa)
+        biotCoefficient (npt.NDArray[np.float64]): Biot coefficient (*b*).
+        porosity (npt.NDArray[np.float64]): Rock porosity (:math:`\phi`).
+
+    Returns:
+        npt.NDArray[np.float64]: Compressibility array (*C* - Pa^-1).
+
+    """
+    assert poissonRatio is not None, "Poisson's ratio must be defined"
+    assert bulkModulus is not None, "Bulk modulus must be defined"
+    assert biotCoefficient is not None, "Biot coefficient must be defined"
+    assert porosity is not None, "Porosity must be defined"
+
+    term1: npt.NDArray[np.float64] = 1.0 - 2.0 * poissonRatio
+
+    mask: npt.NDArray[np.bool_] = (np.abs(bulkModulus) < EPSILON) * (
+        np.abs(porosity) < EPSILON
+    )
+    denFac1: npt.NDArray[np.float64] = porosity * bulkModulus
+    term1[mask] = 1.0
+    term1 /= denFac1
+
+    term2M1: npt.NDArray[np.float64] = (
+        biotCoefficient * biotCoefficient * (1 + poissonRatio)
+    )
+    denTerm2M1: npt.NDArray[np.float64] = 1 - poissonRatio
+    mask2: npt.NDArray[np.bool_] = np.abs(denTerm2M1) < EPSILON
+    denTerm2M1[mask2] = 1.0
+    term2M1 /= denTerm2M1
+    term2M1[mask2] = np.nan
+
+    term2M2: npt.NDArray[np.float64] = (
+        3.0 * (biotCoefficient - porosity) * (1 - biotCoefficient)
+    )
+    term2: npt.NDArray[np.float64] = term2M1 + term2M2
+    return term1 * term2
+
+
+def shearCapacityUtilization(
+    traction: npt.NDArray[np.float64], rockCohesion: float, frictionAngle: float
+) -> npt.NDArray[np.float64]:
+    r"""Compute shear capacity utilization (SCU).
+
+    .. math::
+        SCU = \frac{abs(\tau_1)}{\tau_{max}}
+
+    where \tau_{max} is the Mohr-Coulomb failure threshold.
+
+    Args:
+        traction (npt.NDArray[np.float64]): traction vector
+            (:math:`(\sigma, \tau_1, \tau2)` - Pa)
+        rockCohesion (float): rock cohesion (*c* - Pa).
+        frictionAngle (float): friction angle (:math:`\alpha` - rad).
+
+    Returns:
+        npt.NDArray[np.float64]: *SCU*
+
+    """
+    assert traction is not None, "Traction must be defined"
+    assert traction.shape[1] == 3, "Traction vector must have 3 components."
+
+    scu: npt.NDArray[np.float64] = np.full(traction.shape[0], np.nan)
+    for i, tractionVec in enumerate(traction):
+        # use -1 to agree with Geos convention (i.e., compression with negative stress)
+        stressNormal: npt.NDArray[np.float64] = -1.0 * tractionVec[0]
+
+        # compute failure envelope
+        mohrCoulomb: MohrCoulomb = MohrCoulomb(rockCohesion, frictionAngle)
+        tauFailure: float = float(mohrCoulomb.computeShearStress(stressNormal))
+        scu_i: float = np.nan
+        if tauFailure > 0:
+            scu_i = np.abs(tractionVec[1]) / tauFailure
+        # compute SCU
+        scu[i] = scu_i
+    return scu
+
+
+def computeStressPrincipalComponentsFromStressVector(
+    stressVector: npt.NDArray[np.float64],
+) -> tuple[float, float, float]:
+    """Compute stress principal components from stress vector.
+
+    Args:
+        stressVector (npt.NDArray[np.float64]): stress vector.
+
+    Returns:
+        tuple[float, float, float]: Principal components sorted in ascending
+            order.
+
+    """
+    assert stressVector.size == 6, "Stress vector dimension is wrong."
+    stressTensor: npt.NDArray[np.float64] = getAttributeMatrixFromVector(stressVector)
+    return computeStressPrincipalComponents(stressTensor)
+
+
+def computeStressPrincipalComponents(
+    stressTensor: npt.NDArray[np.float64],
+) -> tuple[float, float, float]:
+    """Compute stress principal components.
+
+    Args:
+        stressTensor (npt.NDArray[np.float64]): stress tensor.
+
+    Returns:
+        tuple[float, float, float]: Principal components sorted in ascending
+            order.
+
+    """
+    # get eigen values
+    e_val, e_vec = np.linalg.eig(stressTensor)
+    # sort principal stresses from smallest to largest
+    p3, p2, p1 = np.sort(e_val)
+    return (p3, p2, p1)
+
+
+def computeNormalShearStress(
+    stressTensor: npt.NDArray[np.float64], directionVector: npt.NDArray[np.float64]
+) -> tuple[float, float]:
+    """Compute normal and shear stress according to stress tensor and direction.
+
+    Args:
+        stressTensor (npt.NDArray[np.float64]): 3x3 stress tensor
+        directionVector (npt.NDArray[np.float64]): direction vector
+
+    Returns:
+        tuple[float, float]: normal and shear stresses.
+
+    """
+    assert stressTensor.shape == (3, 3), "Stress tensor must be 3x3 matrix."
+    assert directionVector.size == 3, "Direction vector must have 3 components"
+
+    # normalization of direction vector
+    directionVector = directionVector / np.linalg.norm(directionVector)
+    # stress vector
+    T: npt.NDArray[np.float64] = np.dot(stressTensor, directionVector)
+    # normal stress
+    sigmaN: float = np.dot(T, directionVector)
+    # shear stress
+    tauVec: npt.NDArray[np.float64] = T - np.dot(sigmaN, directionVector)
+    tau: float = float(np.linalg.norm(tauVec))
+    return (sigmaN, tau)
\ No newline at end of file

From 5dbbf94d381e440be461c4ae4287fed971e1e33a Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Thu, 27 Mar 2025 12:00:44 +0100
Subject: [PATCH 02/14] Update dependencies

---
 .../processing/geomechanicsCalculatorFunctions.py             | 2 +-
 geos-posp/src/PVplugins/PVMohrCirclePlot.py                   | 2 +-
 geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py     | 2 +-
 geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py        | 2 +-
 .../src/geos_posp/visu/mohrCircles/functionsMohrCircle.py     | 4 ++--
 geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py   | 4 ++--
 6 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
index 568174b9..39e43974 100644
--- a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
+++ b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
@@ -4,7 +4,7 @@
 import numpy as np
 import numpy.typing as npt
 
-from geos-geomechanics.processing.MohrCoulomb import MohrCoulomb
+from geos_geomechanics.model.MohrCoulomb import MohrCoulomb
 from geos_posp.utils.geosUtils import getAttributeMatrixFromVector
 from geos_posp.utils.PhysicalConstants import (
     EPSILON,
diff --git a/geos-posp/src/PVplugins/PVMohrCirclePlot.py b/geos-posp/src/PVplugins/PVMohrCirclePlot.py
index ca5a438f..57cf7300 100644
--- a/geos-posp/src/PVplugins/PVMohrCirclePlot.py
+++ b/geos-posp/src/PVplugins/PVMohrCirclePlot.py
@@ -34,7 +34,7 @@
 
 import geos_posp.visu.mohrCircles.functionsMohrCircle as mcf
 import geos_posp.visu.PVUtils.paraviewTreatments as pvt
-from geos_posp.processing.MohrCircle import MohrCircle
+from geos_geomechanics.model.MohrCircle import MohrCircle
 from geos_posp.processing.vtkUtils import getArrayInObject, mergeBlocks
 from geos_posp.utils.enumUnits import Pressure, enumerationDomainUnit
 from geos_posp.utils.GeosOutputsConstants import (
diff --git a/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py b/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
index 0f8e9ba1..81a22117 100644
--- a/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
+++ b/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
@@ -16,7 +16,7 @@
 )
 from vtkmodules.vtkFiltersCore import vtkCellCenters
 
-import geos_posp.processing.geomechanicsCalculatorFunctions as fcts
+import geos_geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 from geos_posp.processing.vtkUtils import (
     createAttribute,
     getArrayInObject,
diff --git a/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py b/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
index d1e2d8a9..fa178c9a 100644
--- a/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
+++ b/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
@@ -17,7 +17,7 @@
     vtkPolyData,
 )
 
-import geos_posp.processing.geomechanicsCalculatorFunctions as fcts
+import geos_geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 import geos_posp.processing.geometryFunctions as geom
 from geos_posp.processing.vtkUtils import (
     createAttribute,
diff --git a/geos-posp/src/geos_posp/visu/mohrCircles/functionsMohrCircle.py b/geos-posp/src/geos_posp/visu/mohrCircles/functionsMohrCircle.py
index 832d8ed5..bfe194cd 100644
--- a/geos-posp/src/geos_posp/visu/mohrCircles/functionsMohrCircle.py
+++ b/geos-posp/src/geos_posp/visu/mohrCircles/functionsMohrCircle.py
@@ -7,8 +7,8 @@
 import numpy as np
 import numpy.typing as npt
 
-from geos_posp.processing.MohrCircle import MohrCircle
-from geos_posp.processing.MohrCoulomb import MohrCoulomb
+from geos_geomechanics.model.MohrCircle import MohrCircle
+from geos_geomechanics.model.MohrCoulomb import MohrCoulomb
 from geos_posp.visu.mohrCircles import (
     MOHR_CIRCLE_ANALYSIS_MAIN,
     MOHR_CIRCLE_PATH,
diff --git a/geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py b/geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py
index c3a19d11..f0c21a77 100644
--- a/geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py
+++ b/geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py
@@ -12,8 +12,8 @@
 from matplotlib.lines import Line2D  # type: ignore[import-untyped]
 
 import geos_posp.visu.mohrCircles.functionsMohrCircle as mcf
-from geos_posp.processing.MohrCircle import MohrCircle
-from geos_posp.processing.MohrCoulomb import MohrCoulomb
+from geos_geomechanics.model.MohrCircle import MohrCircle
+from geos_geomechanics.model.MohrCoulomb import MohrCoulomb
 from geos_posp.utils.enumUnits import Pressure, Unit, convert
 from geos_posp.utils.GeosOutputsConstants import FAILURE_ENVELOPE
 from geos_posp.visu.PVUtils.matplotlibOptions import (

From 708f5a7eb09d878a6e7797927b3e9d8de84c9359 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Thu, 27 Mar 2025 12:03:03 +0100
Subject: [PATCH 03/14] Move tests to geomechanics package

---
 .../tests/testsFunctionsGeomechanicsCalculator.py             | 3 +--
 {geos-posp => geos-geomechanics}/tests/testsMohrCircle.py     | 4 ++--
 2 files changed, 3 insertions(+), 4 deletions(-)
 rename {geos-posp => geos-geomechanics}/tests/testsFunctionsGeomechanicsCalculator.py (96%)
 rename {geos-posp => geos-geomechanics}/tests/testsMohrCircle.py (98%)

diff --git a/geos-posp/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
similarity index 96%
rename from geos-posp/tests/testsFunctionsGeomechanicsCalculator.py
rename to geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
index 9a8b2ead..0a3eb0a1 100644
--- a/geos-posp/tests/testsFunctionsGeomechanicsCalculator.py
+++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
@@ -17,10 +17,9 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
-import geos_posp.processing.geomechanicsCalculatorFunctions as fcts
+import geos_geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 from geos_posp.utils import PhysicalConstants
 from geos_posp.utils.geosUtils import getAttributeMatrixFromVector
-from geos_posp.utils.UnitRepository import Unit, UnitRepository
 
 # geomechanical outputs - Testing variables - Unit is GPa
 bulkModulus: npt.NDArray[np.float64] = np.array([9.0, 50.0, 65.0, np.nan, 150.0])
diff --git a/geos-posp/tests/testsMohrCircle.py b/geos-geomechanics/tests/testsMohrCircle.py
similarity index 98%
rename from geos-posp/tests/testsMohrCircle.py
rename to geos-geomechanics/tests/testsMohrCircle.py
index 83bda4e9..5cf1bf97 100644
--- a/geos-posp/tests/testsMohrCircle.py
+++ b/geos-geomechanics/tests/testsMohrCircle.py
@@ -15,8 +15,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
-from geos_posp.processing.MohrCircle import MohrCircle
-from geos_posp.processing.MohrCoulomb import MohrCoulomb
+from geos_geomechanics.model.MohrCircle import MohrCircle
+from geos_geomechanics.model.MohrCoulomb import MohrCoulomb
 
 circleId = "12453"
 stressVector: npt.NDArray[np.float64] = np.array([1.8, 2.5, 5.2, 0.3, 0.4, 0.1])

From 9f306fdb2ddff3ded9e9259128456c12551f16c3 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Fri, 28 Mar 2025 14:43:21 +0100
Subject: [PATCH 04/14] Update build process

---
 geos-geomechanics/pyproject.toml | 71 ++++++++++++++++++++++++++++++++
 geos-geomechanics/setup.py       | 15 +++++++
 geos-posp/pyproject.toml         |  1 +
 geos-posp/setup.py               | 15 +++++++
 4 files changed, 102 insertions(+)
 create mode 100644 geos-geomechanics/pyproject.toml
 create mode 100644 geos-geomechanics/setup.py
 create mode 100644 geos-posp/setup.py

diff --git a/geos-geomechanics/pyproject.toml b/geos-geomechanics/pyproject.toml
new file mode 100644
index 00000000..72f020f9
--- /dev/null
+++ b/geos-geomechanics/pyproject.toml
@@ -0,0 +1,71 @@
+[build-system]
+requires = ["setuptools>=61.2"]
+build-backend = "setuptools.build_meta"
+
+[tool.setuptools]
+include-package-data = true
+
+[tool.setuptools.packages.find]
+where = ["src"]
+include = ["geos_geomechanics*"]
+exclude = ['tests*']
+
+[project]
+name = "geos-geomechanics"
+version = "0.1.0"
+description = "Geomechanics models and processing tools"
+authors = [{name = "GEOS contributors" }]
+maintainers = [{name = "Martin Lemay", email = "martin.lemay@external.totalenergies.com"}]
+license = {text = "Apache-2.0"}
+classifiers = [
+    "Development Status :: 4 - Beta",
+    "Programming Language :: Python"
+]
+dependencies = [
+    "numpy",
+    "typing_extensions",
+]
+dynamic = ["dependencies"]
+requires-python = ">= 3.9"
+
+[project.optional-dependencies]
+build = [
+  "build ~= 1.2"
+]
+dev = [
+    "pylint",
+    "mypy",
+    "black",
+    "sphinx", 
+    "sphinx-rtd-theme", 
+    "sphinx-autodoc-typehints"
+]
+test = [
+   "pytest",
+   "pytest-cov"
+]
+
+[project.scripts]
+
+
+[tool.pytest.ini_options]
+addopts = "--import-mode=importlib"
+console_output_style = "count"
+pythonpath = [".", "src"]
+python_classes = "Test"
+python_files = "test*.py"
+python_functions = "test*"
+testpaths = ["tests"]
+norecursedirs = "bin"
+filterwarnings = []
+
+
+[tool.coverage.run]
+branch = true
+source = ["geos_geomechanics"]
+
+
+[tool.mypy]
+python_version = "3.9"
+warn_return_any = true
+warn_unused_configs = true
diff --git a/geos-geomechanics/setup.py b/geos-geomechanics/setup.py
new file mode 100644
index 00000000..afc06c24
--- /dev/null
+++ b/geos-geomechanics/setup.py
@@ -0,0 +1,15 @@
+from pathlib import Path
+from setuptools import setup
+
+# This is where you add any fancy path resolution to the local lib:
+geos_utils_path: str = (Path(__file__).parent.parent / "geos-utils").as_uri()
+
+setup(
+    install_requires=[
+        "vtk >= 9.3",
+        "numpy >= 1.26",
+        "pandas >= 2.2",
+        "typing_extensions >= 4.12",
+        f"geos-utils @ {geos_utils_path}",
+    ]
+)
\ No newline at end of file
diff --git a/geos-posp/pyproject.toml b/geos-posp/pyproject.toml
index 3bdbbd26..c40646de 100644
--- a/geos-posp/pyproject.toml
+++ b/geos-posp/pyproject.toml
@@ -35,6 +35,7 @@ dependencies = [
     "pandas >= 2.2",
     "typing_extensions >= 4.12",
 ]
+dynamic = ["dependencies"]
 requires-python = ">= 3.9"
 
 [project.urls]
diff --git a/geos-posp/setup.py b/geos-posp/setup.py
new file mode 100644
index 00000000..08ac268e
--- /dev/null
+++ b/geos-posp/setup.py
@@ -0,0 +1,15 @@
+from pathlib import Path
+from setuptools import setup
+
+# This is where you add any fancy path resolution to the local lib:
+geos_geomecha_path: str = (Path(__file__).parent.parent / "geos-geomechanics").as_uri()
+
+setup(
+    install_requires=[
+        "vtk >= 9.3",
+        "numpy >= 1.26",
+        "pandas >= 2.2",
+        "typing_extensions >= 4.12",
+        f"geos-geomechanics @ {geos_geomecha_path}",
+    ]
+)
\ No newline at end of file

From cd64f590aad3a3cfd7bfc83cb315804077f29b57 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Mon, 31 Mar 2025 09:49:39 +0200
Subject: [PATCH 05/14] Correction to build process

---
 docs/conf.py                     | 4 ++--
 geos-geomechanics/pyproject.toml | 6 +-----
 2 files changed, 3 insertions(+), 7 deletions(-)

diff --git a/docs/conf.py b/docs/conf.py
index 8e22d041..54e31a69 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -17,8 +17,8 @@
 
 # Add python modules to be documented
 python_root = '..'
-python_modules = ( 'geos-ats', 'geos-mesh', 'geos-posp', 'geos-timehistory', 'geos-utils', 'geos-xml-tools', 
-                   'hdf5-wrapper', 'pygeos-tools' )
+python_modules = ( 'geos-ats', 'geos-mesh', 'geos-posp', 'geos-timehistory', 'geos-utils', 'geos-xml-tools', 'hdf5-wrapper',
+                   'pygeos-tools', 'geos-geomechanics' )
 for m in python_modules:
     sys.path.insert( 0, os.path.abspath( os.path.join( python_root, m, 'src' ) ) )
 
diff --git a/geos-geomechanics/pyproject.toml b/geos-geomechanics/pyproject.toml
index 72f020f9..0ba2c2d6 100644
--- a/geos-geomechanics/pyproject.toml
+++ b/geos-geomechanics/pyproject.toml
@@ -21,10 +21,6 @@ classifiers = [
     "Development Status :: 4 - Beta",
     "Programming Language :: Python"
 ]
-dependencies = [
-    "numpy",
-    "typing_extensions",
-]
 dynamic = ["dependencies"]
 requires-python = ">= 3.9"
 
@@ -62,7 +58,7 @@ filterwarnings = []
 
 [tool.coverage.run]
 branch = true
-source = ["geos_geomechanics"]
+source = ["geos/geomechanics"]
 
 
 [tool.mypy]

From 9a7aad5af1d49bc517278a2e6f8acbcc2d36d8d8 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Mon, 31 Mar 2025 11:27:51 +0200
Subject: [PATCH 06/14] Update documentation

---
 docs/geos-geomechanics.rst                 | 10 ++++++++++
 docs/geos_geomechanics_docs/home.rst       |  4 ++++
 docs/geos_geomechanics_docs/model.rst      | 19 +++++++++++++++++++
 docs/geos_geomechanics_docs/modules.rst    | 10 ++++++++++
 docs/geos_geomechanics_docs/processing.rst | 11 +++++++++++
 docs/index.rst                             |  2 ++
 6 files changed, 56 insertions(+)
 create mode 100644 docs/geos-geomechanics.rst
 create mode 100644 docs/geos_geomechanics_docs/home.rst
 create mode 100644 docs/geos_geomechanics_docs/model.rst
 create mode 100644 docs/geos_geomechanics_docs/modules.rst
 create mode 100644 docs/geos_geomechanics_docs/processing.rst

diff --git a/docs/geos-geomechanics.rst b/docs/geos-geomechanics.rst
new file mode 100644
index 00000000..d3d12109
--- /dev/null
+++ b/docs/geos-geomechanics.rst
@@ -0,0 +1,10 @@
+GEOS Post-Processing tools
+=============================
+
+.. toctree::
+   :maxdepth: 5
+   :caption: Contents:
+
+   ./geos_geomechanics_docs/home.rst
+
+   ./geos_geomechanics_docs/modules.rst
\ No newline at end of file
diff --git a/docs/geos_geomechanics_docs/home.rst b/docs/geos_geomechanics_docs/home.rst
new file mode 100644
index 00000000..794b98a5
--- /dev/null
+++ b/docs/geos_geomechanics_docs/home.rst
@@ -0,0 +1,4 @@
+GEOS Geomechanics Tools
+=========================
+
+This package contains geomechanics functions and data models.
\ No newline at end of file
diff --git a/docs/geos_geomechanics_docs/model.rst b/docs/geos_geomechanics_docs/model.rst
new file mode 100644
index 00000000..a7429572
--- /dev/null
+++ b/docs/geos_geomechanics_docs/model.rst
@@ -0,0 +1,19 @@
+Geomechanics models
+======================
+This 
+
+geos_geomechanics.model.MohrCircle module
+------------------------------------------
+
+.. automodule:: geos.geomechanics.model.MohrCircle
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+geos_geomechanics.model.MohrCoulomb module
+-------------------------------------------
+
+.. automodule:: geos.geomechanics.model.MohrCoulomb
+   :members:
+   :undoc-members:
+   :show-inheritance:
\ No newline at end of file
diff --git a/docs/geos_geomechanics_docs/modules.rst b/docs/geos_geomechanics_docs/modules.rst
new file mode 100644
index 00000000..b23b857a
--- /dev/null
+++ b/docs/geos_geomechanics_docs/modules.rst
@@ -0,0 +1,10 @@
+GEOS geomechanics tools
+=========================
+
+.. toctree::
+   :maxdepth: 5
+   
+   model
+   
+   processing
+
diff --git a/docs/geos_geomechanics_docs/processing.rst b/docs/geos_geomechanics_docs/processing.rst
new file mode 100644
index 00000000..9314f6f8
--- /dev/null
+++ b/docs/geos_geomechanics_docs/processing.rst
@@ -0,0 +1,11 @@
+Processing
+============
+The processing part of `geos-geomechanics` package contains functions tools to compute geomechanical properties.
+
+geos_geomechanics.processing.geomechanicsCalculatorFunctions module
+---------------------------------------------------------------
+
+.. automodule:: geos_geomechanics.processing.geomechanicsCalculatorFunctions
+   :members:
+   :undoc-members:
+   :show-inheritance:
\ No newline at end of file
diff --git a/docs/index.rst b/docs/index.rst
index 363634c9..f05ca861 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -84,6 +84,8 @@ Packages
 
     geos-posp
 
+    geos-geomechanics
+
     geos-timehistory
 
     geos-utils

From 37364c10f951f1ab47d1fe164792dee1175b6f27 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Tue, 1 Apr 2025 10:28:56 +0200
Subject: [PATCH 07/14] Build typo and missing files

---
 geos-geomechanics/src/geos/geomechanics/__init__.py             | 0
 geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py     | 2 +-
 geos-geomechanics/src/geos/geomechanics/model/__init__.py       | 0
 geos-geomechanics/src/geos/geomechanics/processing/__init__.py  | 0
 .../geomechanics/processing/geomechanicsCalculatorFunctions.py  | 2 +-
 geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py          | 2 +-
 6 files changed, 3 insertions(+), 3 deletions(-)
 create mode 100644 geos-geomechanics/src/geos/geomechanics/__init__.py
 create mode 100644 geos-geomechanics/src/geos/geomechanics/model/__init__.py
 create mode 100644 geos-geomechanics/src/geos/geomechanics/processing/__init__.py

diff --git a/geos-geomechanics/src/geos/geomechanics/__init__.py b/geos-geomechanics/src/geos/geomechanics/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
index f32ccf44..94857a6c 100644
--- a/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
+++ b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 from typing_extensions import Self
 
-from geos_geomechanics.processing.geomechanicsCalculatorFunctions import (
+from geos.geomechanics.processing.geomechanicsCalculatorFunctions import (
     computeStressPrincipalComponentsFromStressVector,
 )
 
diff --git a/geos-geomechanics/src/geos/geomechanics/model/__init__.py b/geos-geomechanics/src/geos/geomechanics/model/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/geos-geomechanics/src/geos/geomechanics/processing/__init__.py b/geos-geomechanics/src/geos/geomechanics/processing/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
index c2bcf271..2327d574 100644
--- a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
+++ b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
@@ -4,7 +4,7 @@
 import numpy as np
 import numpy.typing as npt
 
-from geos_geomechanics.model.MohrCoulomb import MohrCoulomb
+from geos.geomechanics.model.MohrCoulomb import MohrCoulomb
 from geos.utils.algebraFunctions import getAttributeMatrixFromVector
 from geos.utils.PhysicalConstants import (
     EPSILON,
diff --git a/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py b/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
index 36ac8144..a1c6f094 100644
--- a/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
+++ b/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
@@ -17,7 +17,7 @@
     vtkPolyData,
 )
 
-import geos_geomechanics.processing.geomechanicsCalculatorFunctions as fcts
+import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 
 from geos_posp.processing.vtkUtils import (
     createAttribute,

From 4c2cf4b8fdad4b5ff108fb9d81580bb077f771e6 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Tue, 1 Apr 2025 10:37:57 +0200
Subject: [PATCH 08/14] Doc automodule typo

---
 docs/geos_geomechanics_docs/model.rst      | 4 ++--
 docs/geos_geomechanics_docs/processing.rst | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/geos_geomechanics_docs/model.rst b/docs/geos_geomechanics_docs/model.rst
index a7429572..d1b98d41 100644
--- a/docs/geos_geomechanics_docs/model.rst
+++ b/docs/geos_geomechanics_docs/model.rst
@@ -2,7 +2,7 @@ Geomechanics models
 ======================
 This 
 
-geos_geomechanics.model.MohrCircle module
+geos.geomechanics.model.MohrCircle module
 ------------------------------------------
 
 .. automodule:: geos.geomechanics.model.MohrCircle
@@ -10,7 +10,7 @@ geos_geomechanics.model.MohrCircle module
    :undoc-members:
    :show-inheritance:
 
-geos_geomechanics.model.MohrCoulomb module
+geos.geomechanics.model.MohrCoulomb module
 -------------------------------------------
 
 .. automodule:: geos.geomechanics.model.MohrCoulomb
diff --git a/docs/geos_geomechanics_docs/processing.rst b/docs/geos_geomechanics_docs/processing.rst
index 9314f6f8..35e416a7 100644
--- a/docs/geos_geomechanics_docs/processing.rst
+++ b/docs/geos_geomechanics_docs/processing.rst
@@ -2,10 +2,10 @@ Processing
 ============
 The processing part of `geos-geomechanics` package contains functions tools to compute geomechanical properties.
 
-geos_geomechanics.processing.geomechanicsCalculatorFunctions module
+geos.geomechanics.processing.geomechanicsCalculatorFunctions module
 ---------------------------------------------------------------
 
-.. automodule:: geos_geomechanics.processing.geomechanicsCalculatorFunctions
+.. automodule:: geos.geomechanics.processing.geomechanicsCalculatorFunctions
    :members:
    :undoc-members:
    :show-inheritance:
\ No newline at end of file

From 1cbda7b4577505a20ddf09bbe373b56caf66b678 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Tue, 1 Apr 2025 11:27:33 +0200
Subject: [PATCH 09/14] Add module path to sys path for PV plugins

---
 geos-posp/src/PVplugins/PVAttributeMapping.py      |  2 ++
 .../PVCreateConstantAttributePerRegion.py          |  2 ++
 .../src/PVplugins/PVExtractMergeBlocksVolume.py    |  2 ++
 .../PVplugins/PVExtractMergeBlocksVolumeSurface.py |  2 ++
 .../PVExtractMergeBlocksVolumeSurfaceWell.py       |  2 ++
 .../PVplugins/PVExtractMergeBlocksVolumeWell.py    |  2 ++
 geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py  |  2 ++
 .../src/PVplugins/PVGeomechanicsWorkflowVolume.py  |  2 ++
 .../PVGeomechanicsWorkflowVolumeSurface.py         |  1 +
 .../PVGeomechanicsWorkflowVolumeSurfaceWell.py     |  1 +
 .../PVplugins/PVGeomechanicsWorkflowVolumeWell.py  |  2 ++
 geos-posp/src/PVplugins/PVGeosLogReader.py         |  2 ++
 geos-posp/src/PVplugins/PVMergeBlocksEnhanced.py   |  2 ++
 geos-posp/src/PVplugins/PVMohrCirclePlot.py        |  4 +++-
 .../src/PVplugins/PVPythonViewConfigurator.py      |  2 ++
 geos-posp/src/PVplugins/PVSurfaceGeomechanics.py   |  2 ++
 .../PVplugins/PVTransferAttributesVolumeSurface.py |  2 ++
 geos-posp/src/PVplugins/__init__.py                | 14 ++++++++++++++
 18 files changed, 47 insertions(+), 1 deletion(-)

diff --git a/geos-posp/src/PVplugins/PVAttributeMapping.py b/geos-posp/src/PVplugins/PVAttributeMapping.py
index d8234423..f1901de1 100644
--- a/geos-posp/src/PVplugins/PVAttributeMapping.py
+++ b/geos-posp/src/PVplugins/PVAttributeMapping.py
@@ -12,6 +12,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
     smdomain,
diff --git a/geos-posp/src/PVplugins/PVCreateConstantAttributePerRegion.py b/geos-posp/src/PVplugins/PVCreateConstantAttributePerRegion.py
index bb8fbbdf..ff8f5520 100644
--- a/geos-posp/src/PVplugins/PVCreateConstantAttributePerRegion.py
+++ b/geos-posp/src/PVplugins/PVCreateConstantAttributePerRegion.py
@@ -15,6 +15,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 import vtkmodules.util.numpy_support as vnp
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
diff --git a/geos-posp/src/PVplugins/PVExtractMergeBlocksVolume.py b/geos-posp/src/PVplugins/PVExtractMergeBlocksVolume.py
index 5423bbe9..8b5ac542 100644
--- a/geos-posp/src/PVplugins/PVExtractMergeBlocksVolume.py
+++ b/geos-posp/src/PVplugins/PVExtractMergeBlocksVolume.py
@@ -16,6 +16,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
     smdomain,
diff --git a/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeSurface.py b/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeSurface.py
index 62d35ea4..77f9d005 100644
--- a/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeSurface.py
+++ b/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeSurface.py
@@ -16,6 +16,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
     smdomain,
diff --git a/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeSurfaceWell.py b/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeSurfaceWell.py
index 4ab7010a..20fce2f1 100644
--- a/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeSurfaceWell.py
+++ b/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeSurfaceWell.py
@@ -16,6 +16,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
     smdomain,
diff --git a/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeWell.py b/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeWell.py
index 2c1c5652..cc6860f8 100644
--- a/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeWell.py
+++ b/geos-posp/src/PVplugins/PVExtractMergeBlocksVolumeWell.py
@@ -23,6 +23,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from geos_posp.filters.GeosBlockExtractor import GeosBlockExtractor
 from geos_posp.filters.GeosBlockMerge import GeosBlockMerge
 from geos_posp.processing.vtkUtils import (
diff --git a/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py b/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py
index 25b6f9d5..c6e90add 100644
--- a/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py
+++ b/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py
@@ -23,6 +23,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from geos.utils.Logger import Logger, getLogger
 from geos.utils.PhysicalConstants import (
     DEFAULT_FRICTION_ANGLE_DEG,
diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolume.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolume.py
index 92c0142a..6cb49c94 100644
--- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolume.py
+++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolume.py
@@ -17,6 +17,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
     smdomain,
diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py
index d498d767..a5451efb 100644
--- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py
+++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py
@@ -17,6 +17,7 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
 
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py
index bbfdd053..a5f798b2 100644
--- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py
+++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py
@@ -17,6 +17,7 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
 
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeWell.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeWell.py
index f2c01541..29ac5e48 100644
--- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeWell.py
+++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeWell.py
@@ -17,6 +17,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
     smdomain,
diff --git a/geos-posp/src/PVplugins/PVGeosLogReader.py b/geos-posp/src/PVplugins/PVGeosLogReader.py
index a5dd35dd..7a5b65d2 100644
--- a/geos-posp/src/PVplugins/PVGeosLogReader.py
+++ b/geos-posp/src/PVplugins/PVGeosLogReader.py
@@ -17,6 +17,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 import vtkmodules.util.numpy_support as vnp
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
diff --git a/geos-posp/src/PVplugins/PVMergeBlocksEnhanced.py b/geos-posp/src/PVplugins/PVMergeBlocksEnhanced.py
index 4b8a7c86..7c7c155d 100644
--- a/geos-posp/src/PVplugins/PVMergeBlocksEnhanced.py
+++ b/geos-posp/src/PVplugins/PVMergeBlocksEnhanced.py
@@ -12,6 +12,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
     smdomain,
diff --git a/geos-posp/src/PVplugins/PVMohrCirclePlot.py b/geos-posp/src/PVplugins/PVMohrCirclePlot.py
index 9527eafb..eb4c73b4 100644
--- a/geos-posp/src/PVplugins/PVMohrCirclePlot.py
+++ b/geos-posp/src/PVplugins/PVMohrCirclePlot.py
@@ -32,9 +32,11 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 import geos_posp.visu.mohrCircles.functionsMohrCircle as mcf
 import geos_posp.visu.PVUtils.paraviewTreatments as pvt
-from geos_geomechanics.model.MohrCircle import MohrCircle
+from geos.geomechanics.model.MohrCircle import MohrCircle
 from geos_posp.processing.vtkUtils import getArrayInObject, mergeBlocks
 from geos.utils.enumUnits import Pressure, enumerationDomainUnit
 from geos.utils.GeosOutputsConstants import (
diff --git a/geos-posp/src/PVplugins/PVPythonViewConfigurator.py b/geos-posp/src/PVplugins/PVPythonViewConfigurator.py
index 9d63e707..e8c850ce 100644
--- a/geos-posp/src/PVplugins/PVPythonViewConfigurator.py
+++ b/geos-posp/src/PVplugins/PVPythonViewConfigurator.py
@@ -14,6 +14,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.simple import (  # type: ignore[import-not-found]
     GetActiveSource,
     GetActiveView,
diff --git a/geos-posp/src/PVplugins/PVSurfaceGeomechanics.py b/geos-posp/src/PVplugins/PVSurfaceGeomechanics.py
index 2ff1b4f5..cc91a209 100644
--- a/geos-posp/src/PVplugins/PVSurfaceGeomechanics.py
+++ b/geos-posp/src/PVplugins/PVSurfaceGeomechanics.py
@@ -13,6 +13,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.util.vtkAlgorithm import (  # type: ignore[import-not-found]
     VTKPythonAlgorithmBase,
     smdomain,
diff --git a/geos-posp/src/PVplugins/PVTransferAttributesVolumeSurface.py b/geos-posp/src/PVplugins/PVTransferAttributesVolumeSurface.py
index 2bc6e1ca..26af9867 100644
--- a/geos-posp/src/PVplugins/PVTransferAttributesVolumeSurface.py
+++ b/geos-posp/src/PVplugins/PVTransferAttributesVolumeSurface.py
@@ -12,6 +12,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
+import PVplugins #required to update sys path
+
 from paraview.simple import (  # type: ignore[import-not-found]
     FindSource,
     GetActiveSource,
diff --git a/geos-posp/src/PVplugins/__init__.py b/geos-posp/src/PVplugins/__init__.py
index e69de29b..7cd2f2f3 100644
--- a/geos-posp/src/PVplugins/__init__.py
+++ b/geos-posp/src/PVplugins/__init__.py
@@ -0,0 +1,14 @@
+import os
+import sys
+
+
+# Add other packages path to sys path
+dir_path = os.path.dirname(os.path.realpath(__file__))
+python_root = '../../..'
+
+python_modules = ( 'geos-posp', 'geos-utils', 'geos-geomechanics' )
+
+for m in python_modules:
+    m_path = os.path.abspath( os.path.join( dir_path, python_root, m, 'src' ) )
+    if m_path not in sys.path:
+        sys.path.insert( 0, m_path)
\ No newline at end of file

From 101b7d7e80f91adfc37c1703022dd0c55de1eec6 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Tue, 1 Apr 2025 11:29:35 +0200
Subject: [PATCH 10/14] Typo in module import

---
 .../tests/testsFunctionsGeomechanicsCalculator.py             | 2 +-
 geos-geomechanics/tests/testsMohrCircle.py                    | 4 ++--
 geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py     | 2 +-
 .../src/geos_posp/visu/mohrCircles/functionsMohrCircle.py     | 4 ++--
 geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py   | 4 ++--
 5 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
index a18972fb..827cc815 100644
--- a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
+++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
@@ -17,7 +17,7 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
-import geos_geomechanics.processing.geomechanicsCalculatorFunctions as fcts
+import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 from geos.utils import PhysicalConstants
 from geos.utils.algebraFunctions import getAttributeMatrixFromVector
 
diff --git a/geos-geomechanics/tests/testsMohrCircle.py b/geos-geomechanics/tests/testsMohrCircle.py
index 5cf1bf97..133b30b3 100644
--- a/geos-geomechanics/tests/testsMohrCircle.py
+++ b/geos-geomechanics/tests/testsMohrCircle.py
@@ -15,8 +15,8 @@
 if parent_dir_path not in sys.path:
     sys.path.append(parent_dir_path)
 
-from geos_geomechanics.model.MohrCircle import MohrCircle
-from geos_geomechanics.model.MohrCoulomb import MohrCoulomb
+from geos.geomechanics.model.MohrCircle import MohrCircle
+from geos.geomechanics.model.MohrCoulomb import MohrCoulomb
 
 circleId = "12453"
 stressVector: npt.NDArray[np.float64] = np.array([1.8, 2.5, 5.2, 0.3, 0.4, 0.1])
diff --git a/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py b/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
index 71d22d57..720771f7 100644
--- a/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
+++ b/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
@@ -16,7 +16,7 @@
 )
 from vtkmodules.vtkFiltersCore import vtkCellCenters
 
-import geos_geomechanics.processing.geomechanicsCalculatorFunctions as fcts
+import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 from geos_posp.processing.vtkUtils import (
     createAttribute,
     getArrayInObject,
diff --git a/geos-posp/src/geos_posp/visu/mohrCircles/functionsMohrCircle.py b/geos-posp/src/geos_posp/visu/mohrCircles/functionsMohrCircle.py
index bfe194cd..f41a43ce 100644
--- a/geos-posp/src/geos_posp/visu/mohrCircles/functionsMohrCircle.py
+++ b/geos-posp/src/geos_posp/visu/mohrCircles/functionsMohrCircle.py
@@ -7,8 +7,8 @@
 import numpy as np
 import numpy.typing as npt
 
-from geos_geomechanics.model.MohrCircle import MohrCircle
-from geos_geomechanics.model.MohrCoulomb import MohrCoulomb
+from geos.geomechanics.model.MohrCircle import MohrCircle
+from geos.geomechanics.model.MohrCoulomb import MohrCoulomb
 from geos_posp.visu.mohrCircles import (
     MOHR_CIRCLE_ANALYSIS_MAIN,
     MOHR_CIRCLE_PATH,
diff --git a/geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py b/geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py
index 70962ce3..8e22a18e 100644
--- a/geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py
+++ b/geos-posp/src/geos_posp/visu/mohrCircles/plotMohrCircles.py
@@ -12,8 +12,8 @@
 from matplotlib.lines import Line2D  # type: ignore[import-untyped]
 
 import geos_posp.visu.mohrCircles.functionsMohrCircle as mcf
-from geos_geomechanics.model.MohrCircle import MohrCircle
-from geos_geomechanics.model.MohrCoulomb import MohrCoulomb
+from geos.geomechanics.model.MohrCircle import MohrCircle
+from geos.geomechanics.model.MohrCoulomb import MohrCoulomb
 from geos.utils.enumUnits import Pressure, Unit, convert
 from geos.utils.GeosOutputsConstants import FAILURE_ENVELOPE
 from geos_posp.visu.PVUtils.matplotlibOptions import (

From 0750d9618adaccc889442f83e4556ad9fcc47e83 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Wed, 2 Apr 2025 11:25:56 +0200
Subject: [PATCH 11/14] Fix import path and move line to top of file

---
 geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py         | 3 ++-
 geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py b/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py
index c6e90add..e5cd865a 100644
--- a/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py
+++ b/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py
@@ -18,6 +18,8 @@
 from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector
 from vtkmodules.vtkCommonDataModel import vtkPointSet, vtkUnstructuredGrid
 
+from geos_posp.filters.GeomechanicsCalculator import GeomechanicsCalculator
+
 dir_path = os.path.dirname(os.path.realpath(__file__))
 parent_dir_path = os.path.dirname(dir_path)
 if parent_dir_path not in sys.path:
@@ -342,7 +344,6 @@ def RequestData(
         """
         try:
             self.m_logger.info(f"Apply filter {__name__}")
-            from filters.GeomechanicsCalculator import GeomechanicsCalculator
 
             inData = self.GetInputData(inInfoVec, 0, 0)
             assert inData is not None
diff --git a/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py b/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
index 720771f7..b8afd5cb 100644
--- a/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
+++ b/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
@@ -50,7 +50,7 @@
 .. code-block:: python
 
     import numpy as np
-    from filters.GeomechanicsCalculator import GeomechanicsCalculator
+    from geos_posp.filters.GeomechanicsCalculator import GeomechanicsCalculator
 
     # filter inputs
     logger :Logger

From 4424031227a4b2f61f6aed0b3a4aa458672ff472 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Wed, 2 Apr 2025 11:46:58 +0200
Subject: [PATCH 12/14] linting

---
 geos-geomechanics/setup.py                    |  18 +-
 .../src/geos/geomechanics/model/MohrCircle.py |  38 +-
 .../geos/geomechanics/model/MohrCoulomb.py    |  27 +-
 .../geomechanicsCalculatorFunctions.py        | 590 ++++++++----------
 .../testsFunctionsGeomechanicsCalculator.py   | 500 ++++++---------
 geos-geomechanics/tests/testsMohrCircle.py    | 420 ++++++-------
 6 files changed, 684 insertions(+), 909 deletions(-)

diff --git a/geos-geomechanics/setup.py b/geos-geomechanics/setup.py
index afc06c24..3c542ff6 100644
--- a/geos-geomechanics/setup.py
+++ b/geos-geomechanics/setup.py
@@ -2,14 +2,12 @@
 from setuptools import setup
 
 # This is where you add any fancy path resolution to the local lib:
-geos_utils_path: str = (Path(__file__).parent.parent / "geos-utils").as_uri()
+geos_utils_path: str = ( Path( __file__ ).parent.parent / "geos-utils" ).as_uri()
 
-setup(
-    install_requires=[
-        "vtk >= 9.3",
-        "numpy >= 1.26",
-        "pandas >= 2.2",
-        "typing_extensions >= 4.12",
-        f"geos-utils @ {geos_utils_path}",
-    ]
-)
\ No newline at end of file
+setup( install_requires=[
+    "vtk >= 9.3",
+    "numpy >= 1.26",
+    "pandas >= 2.2",
+    "typing_extensions >= 4.12",
+    f"geos-utils @ {geos_utils_path}",
+] )
diff --git a/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
index 94857a6c..a917c3de 100644
--- a/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
+++ b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
@@ -6,8 +6,7 @@
 from typing_extensions import Self
 
 from geos.geomechanics.processing.geomechanicsCalculatorFunctions import (
-    computeStressPrincipalComponentsFromStressVector,
-)
+    computeStressPrincipalComponentsFromStressVector, )
 
 __doc__ = """
 MohrCircle module define the Mohr's circle parameters.
@@ -43,7 +42,8 @@
 
 
 class MohrCircle:
-    def __init__(self: Self, circleId: str) -> None:
+
+    def __init__( self: Self, circleId: str ) -> None:
         """Compute Mohr's Circle from input stress.
 
         Args:
@@ -55,15 +55,15 @@ def __init__(self: Self, circleId: str) -> None:
         self.m_p2: float = 0.0
         self.m_p3: float = 0.0
 
-    def __str__(self: Self) -> str:
+    def __str__( self: Self ) -> str:
         """Overload of __str__ method."""
         return self.m_circleId
 
-    def __repr__(self: Self) -> str:
+    def __repr__( self: Self ) -> str:
         """Overload of __repr__ method."""
         return self.m_circleId
 
-    def setCircleId(self: Self, circleId: str) -> None:
+    def setCircleId( self: Self, circleId: str ) -> None:
         """Set circle Id variable.
 
         Args:
@@ -71,7 +71,7 @@ def setCircleId(self: Self, circleId: str) -> None:
         """
         self.m_circleId = circleId
 
-    def getCircleId(self: Self) -> str:
+    def getCircleId( self: Self ) -> str:
         """Access the Id of the Mohr circle.
 
         Returns:
@@ -79,19 +79,19 @@ def getCircleId(self: Self) -> str:
         """
         return self.m_circleId
 
-    def getCircleRadius(self: Self) -> float:
+    def getCircleRadius( self: Self ) -> float:
         """Compute and return Mohr's circle radius from principal components."""
-        return (self.m_p1 - self.m_p3) / 2.0
+        return ( self.m_p1 - self.m_p3 ) / 2.0
 
-    def getCircleCenter(self: Self) -> float:
+    def getCircleCenter( self: Self ) -> float:
         """Compute and return Mohr's circle center from principal components."""
-        return (self.m_p1 + self.m_p3) / 2.0
+        return ( self.m_p1 + self.m_p3 ) / 2.0
 
-    def getPrincipalComponents(self: Self) -> tuple[float, float, float]:
+    def getPrincipalComponents( self: Self ) -> tuple[ float, float, float ]:
         """Get Moh's circle principal components."""
-        return (self.m_p3, self.m_p2, self.m_p1)
+        return ( self.m_p3, self.m_p2, self.m_p1 )
 
-    def setPrincipalComponents(self: Self, p3: float, p2: float, p1: float) -> None:
+    def setPrincipalComponents( self: Self, p3: float, p2: float, p1: float ) -> None:
         """Set principal components.
 
         Args:
@@ -99,14 +99,12 @@ def setPrincipalComponents(self: Self, p3: float, p2: float, p1: float) -> None:
             p2 (float): second component.
             p1 (float): third component. Must be the greatest.
         """
-        assert (p3 <= p2) and (p2 <= p1), "Component order is wrong."
+        assert ( p3 <= p2 ) and ( p2 <= p1 ), "Component order is wrong."
         self.m_p3 = p3
         self.m_p2 = p2
         self.m_p1 = p1
 
-    def computePrincipalComponents(
-        self: Self, stressVector: npt.NDArray[np.float64]
-    ) -> None:
+    def computePrincipalComponents( self: Self, stressVector: npt.NDArray[ np.float64 ] ) -> None:
         """Calculate principal components.
 
         Args:
@@ -114,6 +112,4 @@ def computePrincipalComponents(
                 Stress vector must follow GEOS convention (XX, YY, ZZ, YZ, XZ, XY)
         """
         # get stress principal components
-        self.m_p3, self.m_p2, self.m_p1 = (
-            computeStressPrincipalComponentsFromStressVector(stressVector)
-        )
\ No newline at end of file
+        self.m_p3, self.m_p2, self.m_p1 = ( computeStressPrincipalComponentsFromStressVector( stressVector ) )
diff --git a/geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py b/geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py
index 4ed0efcf..93409149 100644
--- a/geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py
+++ b/geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py
@@ -37,7 +37,8 @@
 
 
 class MohrCoulomb:
-    def __init__(self: Self, rockCohesion: float, frictionAngle: float) -> None:
+
+    def __init__( self: Self, rockCohesion: float, frictionAngle: float ) -> None:
         """Define Mohr-Coulomb failure envelop.
 
         Args:
@@ -47,13 +48,13 @@ def __init__(self: Self, rockCohesion: float, frictionAngle: float) -> None:
         # rock cohesion
         self.m_rockCohesion = rockCohesion
         # failure envelop slope
-        self.m_slope: float = np.tan(frictionAngle)
+        self.m_slope: float = np.tan( frictionAngle )
         # intersection of failure envelop and abscissa axis
         self.m_sigmaMin: float = -rockCohesion / self.m_slope
 
     def computeShearStress(
-        self: Self, stressNormal0: Union[float, npt.NDArray[np.float64]]
-    ) -> Union[float, npt.NDArray[np.float64]]:
+            self: Self,
+            stressNormal0: Union[ float, npt.NDArray[ np.float64 ] ] ) -> Union[ float, npt.NDArray[ np.float64 ] ]:
         """Compute shear stress from normal stress.
 
         Args:
@@ -63,16 +64,16 @@ def computeShearStress(
         Returns:
             float | npt.NDArray[np.float64]): shear stress value.
         """
-        stressNormal: npt.NDArray[np.float64] = np.array(stressNormal0)
-        tau: npt.NDArray[np.float64] = self.m_slope * stressNormal + self.m_rockCohesion
+        stressNormal: npt.NDArray[ np.float64 ] = np.array( stressNormal0 )
+        tau: npt.NDArray[ np.float64 ] = self.m_slope * stressNormal + self.m_rockCohesion
         return tau
 
     def computeFailureEnvelop(
         self: Self,
         stressNormalMax: float,
-        stressNormalMin: Union[float, None] = None,
+        stressNormalMin: Union[ float, None ] = None,
         n: int = 100,
-    ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
+    ) -> tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]:
         """Compute the envelope failure between min and max normal stress.
 
         Args:
@@ -90,10 +91,6 @@ def computeFailureEnvelop(
             envelop coordinates where first element is the abscissa and second
             element is the ordinates.
         """
-        sigmaMin: float = (
-            self.m_sigmaMin if stressNormalMin is None else stressNormalMin
-        )
-        stressNormal: npt.NDArray[np.float64] = np.linspace(
-            sigmaMin, stressNormalMax, n
-        )
-        return (stressNormal, np.array(self.computeShearStress(stressNormal)))
\ No newline at end of file
+        sigmaMin: float = ( self.m_sigmaMin if stressNormalMin is None else stressNormalMin )
+        stressNormal: npt.NDArray[ np.float64 ] = np.linspace( sigmaMin, stressNormalMax, n )
+        return ( stressNormal, np.array( self.computeShearStress( stressNormal ) ) )
diff --git a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
index 2327d574..48317c34 100644
--- a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
+++ b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
@@ -7,15 +7,12 @@
 from geos.geomechanics.model.MohrCoulomb import MohrCoulomb
 from geos.utils.algebraFunctions import getAttributeMatrixFromVector
 from geos.utils.PhysicalConstants import (
-    EPSILON,
-)
+    EPSILON, )
 
 __doc__ = """Functions to compute additional geomechanical properties."""
 
 
-def specificGravity(
-    density: npt.NDArray[np.float64], specificDensity: float
-) -> npt.NDArray[np.float64]:
+def specificGravity( density: npt.NDArray[ np.float64 ], specificDensity: float ) -> npt.NDArray[ np.float64 ]:
     r"""Compute the specific gravity.
 
     .. math::
@@ -31,15 +28,14 @@ def specificGravity(
     """
     assert density is not None, "Density data must be defined"
 
-    if abs(specificDensity) < EPSILON:
-        return np.full_like(density, np.nan)
+    if abs( specificDensity ) < EPSILON:
+        return np.full_like( density, np.nan )
     return density / specificDensity
 
 
 # https://en.wikipedia.org/wiki/Elastic_modulus
-def youngModulus(
-    bulkModulus: npt.NDArray[np.float64], shearModulus: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def youngModulus( bulkModulus: npt.NDArray[ np.float64 ],
+                  shearModulus: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute Young modulus.
 
     .. math::
@@ -56,22 +52,19 @@ def youngModulus(
     assert bulkModulus is not None, "Bulk modulus must be defined"
     assert shearModulus is not None, "Shear modulus must be defined"
     # manage division by 0 by replacing with nan
-    assert bulkModulus.size == shearModulus.size, (
-        "Bulk modulus array and Shear modulus array"
-        + " sizes (i.e., number of cells) must be equal."
-    )
-
-    den: npt.NDArray[np.float64] = 3.0 * bulkModulus + shearModulus
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    young: npt.NDArray[np.float64] = 9.0 * bulkModulus * shearModulus / den
-    young[mask] = np.nan
+    assert bulkModulus.size == shearModulus.size, ( "Bulk modulus array and Shear modulus array" +
+                                                    " sizes (i.e., number of cells) must be equal." )
+
+    den: npt.NDArray[ np.float64 ] = 3.0 * bulkModulus + shearModulus
+    mask: npt.NDArray[ np.bool_ ] = np.abs( den ) < EPSILON
+    den[ mask ] = 1.0
+    young: npt.NDArray[ np.float64 ] = 9.0 * bulkModulus * shearModulus / den
+    young[ mask ] = np.nan
     return young
 
 
-def poissonRatio(
-    bulkModulus: npt.NDArray[np.float64], shearModulus: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def poissonRatio( bulkModulus: npt.NDArray[ np.float64 ],
+                  shearModulus: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute Poisson's ratio.
 
     .. math::
@@ -87,22 +80,19 @@ def poissonRatio(
     """
     assert bulkModulus is not None, "Bulk modulus must be defined"
     assert shearModulus is not None, "Shear modulus must be defined"
-    assert bulkModulus.size == shearModulus.size, (
-        "Bulk modulus array and Shear modulus array"
-        + " sizes (i.e., number of cells) must be equal."
-    )
+    assert bulkModulus.size == shearModulus.size, ( "Bulk modulus array and Shear modulus array" +
+                                                    " sizes (i.e., number of cells) must be equal." )
     # manage division by 0 by replacing with nan
-    den: npt.NDArray[np.float64] = 2.0 * (3.0 * bulkModulus + shearModulus)
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    poisson: npt.NDArray[np.float64] = (3.0 * bulkModulus - 2.0 * shearModulus) / den
-    poisson[mask] = np.nan
+    den: npt.NDArray[ np.float64 ] = 2.0 * ( 3.0 * bulkModulus + shearModulus )
+    mask: npt.NDArray[ np.bool_ ] = np.abs( den ) < EPSILON
+    den[ mask ] = 1.0
+    poisson: npt.NDArray[ np.float64 ] = ( 3.0 * bulkModulus - 2.0 * shearModulus ) / den
+    poisson[ mask ] = np.nan
     return poisson
 
 
-def bulkModulus(
-    youngModulus: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def bulkModulus( youngModulus: npt.NDArray[ np.float64 ],
+                 poissonRatio: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute bulk Modulus from young modulus and poisson ratio.
 
     .. math::
@@ -119,17 +109,16 @@ def bulkModulus(
     assert poissonRatio is not None, "Poisson's ratio must be defined"
     assert youngModulus is not None, "Young modulus must be defined"
 
-    den: npt.NDArray[np.float64] = 3.0 * (1.0 - 2.0 * poissonRatio)
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    bulkMod: npt.NDArray[np.float64] = youngModulus / den
-    bulkMod[mask] = np.nan
+    den: npt.NDArray[ np.float64 ] = 3.0 * ( 1.0 - 2.0 * poissonRatio )
+    mask: npt.NDArray[ np.bool_ ] = np.abs( den ) < EPSILON
+    den[ mask ] = 1.0
+    bulkMod: npt.NDArray[ np.float64 ] = youngModulus / den
+    bulkMod[ mask ] = np.nan
     return bulkMod
 
 
-def shearModulus(
-    youngModulus: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def shearModulus( youngModulus: npt.NDArray[ np.float64 ],
+                  poissonRatio: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute shear Modulus from young modulus and poisson ratio.
 
     .. math::
@@ -146,17 +135,16 @@ def shearModulus(
     assert poissonRatio is not None, "Poisson's ratio must be defined"
     assert youngModulus is not None, "Young modulus must be defined"
 
-    den: npt.NDArray[np.float64] = 2.0 * (1.0 + poissonRatio)
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    shearMod: npt.NDArray[np.float64] = youngModulus / den
-    shearMod[mask] = np.nan
+    den: npt.NDArray[ np.float64 ] = 2.0 * ( 1.0 + poissonRatio )
+    mask: npt.NDArray[ np.bool_ ] = np.abs( den ) < EPSILON
+    den[ mask ] = 1.0
+    shearMod: npt.NDArray[ np.float64 ] = youngModulus / den
+    shearMod[ mask ] = np.nan
     return shearMod
 
 
-def lambdaCoefficient(
-    youngModulus: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def lambdaCoefficient( youngModulus: npt.NDArray[ np.float64 ],
+                       poissonRatio: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute lambda coefficient from young modulus and Poisson ratio.
 
     .. math::
@@ -170,18 +158,17 @@ def lambdaCoefficient(
         npt.NDArray[np.float64]: lambda coefficient (:math:`\lambda`)
 
     """
-    lambdaCoeff: npt.NDArray[np.float64] = poissonRatio * youngModulus
-    den: npt.NDArray[np.float64] = (1.0 + poissonRatio) * (1.0 - 2.0 * poissonRatio)
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
+    lambdaCoeff: npt.NDArray[ np.float64 ] = poissonRatio * youngModulus
+    den: npt.NDArray[ np.float64 ] = ( 1.0 + poissonRatio ) * ( 1.0 - 2.0 * poissonRatio )
+    mask: npt.NDArray[ np.bool_ ] = np.abs( den ) < EPSILON
+    den[ mask ] = 1.0
     lambdaCoeff /= den
-    lambdaCoeff[mask] = np.nan
+    lambdaCoeff[ mask ] = np.nan
     return lambdaCoeff
 
 
-def oedometricModulus(
-    Edef: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def oedometricModulus( Edef: npt.NDArray[ np.float64 ],
+                       poissonRatio: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute Oedometric modulus.
 
     .. math::
@@ -198,25 +185,19 @@ def oedometricModulus(
     assert poissonRatio is not None, "Poisson's ratio must be defined"
     assert Edef is not None, "Deformation modulus must be defined"
 
-    assert Edef.size == poissonRatio.size, (
-        "Deformation modulus array and Poisson's"
-        + " ratio array sizes (i.e., number of cells) must be equal."
-    )
-    den: npt.NDArray[np.float64] = 1.0 - (2.0 * poissonRatio * poissonRatio) / (
-        1.0 - poissonRatio
-    )
+    assert Edef.size == poissonRatio.size, ( "Deformation modulus array and Poisson's" +
+                                             " ratio array sizes (i.e., number of cells) must be equal." )
+    den: npt.NDArray[ np.float64 ] = 1.0 - ( 2.0 * poissonRatio * poissonRatio ) / ( 1.0 - poissonRatio )
 
     # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    EodMod: npt.NDArray[np.float64] = Edef / den
-    EodMod[mask] = np.nan
+    mask: npt.NDArray[ np.bool_ ] = np.abs( den ) < EPSILON
+    den[ mask ] = 1.0
+    EodMod: npt.NDArray[ np.float64 ] = Edef / den
+    EodMod[ mask ] = np.nan
     return EodMod
 
 
-def biotCoefficient(
-    Kgrain: float, bulkModulus: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def biotCoefficient( Kgrain: float, bulkModulus: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute Biot coefficient.
 
     .. math::
@@ -233,19 +214,19 @@ def biotCoefficient(
     assert bulkModulus is not None, "Bulk modulus must be defined"
 
     # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(Kgrain) < EPSILON
-    den: npt.NDArray[np.float64] = np.copy(Kgrain)
-    den[mask] = 1.0
-    biot: npt.NDArray[np.float64] = 1.0 - bulkModulus / den
-    biot[mask] = np.nan
+    mask: npt.NDArray[ np.bool_ ] = np.abs( Kgrain ) < EPSILON
+    den: npt.NDArray[ np.float64 ] = np.copy( Kgrain )
+    den[ mask ] = 1.0
+    biot: npt.NDArray[ np.float64 ] = 1.0 - bulkModulus / den
+    biot[ mask ] = np.nan
     return biot
 
 
 def totalStress(
-    effectiveStress: npt.NDArray[np.float64],
-    biot: npt.NDArray[np.float64],
-    pressure: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
+    effectiveStress: npt.NDArray[ np.float64 ],
+    biot: npt.NDArray[ np.float64 ],
+    pressure: npt.NDArray[ np.float64 ],
+) -> npt.NDArray[ np.float64 ]:
     r"""Compute total stress from effective stress, pressure, and Biot coeff.
 
     .. math::
@@ -265,27 +246,22 @@ def totalStress(
     assert biot is not None, "Biot coefficient must be defined"
     assert pressure is not None, "Pressure must be defined"
 
-    assert effectiveStress.shape[0] == biot.size, (
-        "Biot coefficient array and "
-        + "effective stress sizes (i.e., number of cells) must be equal."
-    )
-    assert biot.size == pressure.size, (
-        "Biot coefficient array and pressure array"
-        + "sizes (i.e., number of cells) must be equal."
-    )
+    assert effectiveStress.shape[ 0 ] == biot.size, ( "Biot coefficient array and " +
+                                                      "effective stress sizes (i.e., number of cells) must be equal." )
+    assert biot.size == pressure.size, ( "Biot coefficient array and pressure array" +
+                                         "sizes (i.e., number of cells) must be equal." )
 
-    totalStress: npt.NDArray[np.float64] = np.copy(effectiveStress)
+    totalStress: npt.NDArray[ np.float64 ] = np.copy( effectiveStress )
     # pore pressure has an effect on normal stresses only
     # (cf. https://dnicolasespinoza.github.io/node5.html)
-    nb: int = totalStress.shape[1] if totalStress.shape[1] < 4 else 3
-    for j in range(nb):
-        totalStress[:, j] = effectiveStress[:, j] - biot * pressure
+    nb: int = totalStress.shape[ 1 ] if totalStress.shape[ 1 ] < 4 else 3
+    for j in range( nb ):
+        totalStress[ :, j ] = effectiveStress[ :, j ] - biot * pressure
     return totalStress
 
 
-def stressRatio(
-    horizontalStress: npt.NDArray[np.float64], verticalStress: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def stressRatio( horizontalStress: npt.NDArray[ np.float64 ],
+                 verticalStress: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute horizontal to vertical stress ratio.
 
     .. math::
@@ -305,22 +281,19 @@ def stressRatio(
     assert verticalStress is not None, "Vertical stress must be defined"
 
     assert horizontalStress.size == verticalStress.size, (
-        "Horizontal stress array "
-        + "and vertical stress array sizes (i.e., number of cells) must be equal."
-    )
+        "Horizontal stress array " + "and vertical stress array sizes (i.e., number of cells) must be equal." )
 
     # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(verticalStress) < EPSILON
-    verticalStress2: npt.NDArray[np.float64] = np.copy(verticalStress)
-    verticalStress2[mask] = 1.0
-    ratio: npt.NDArray[np.float64] = horizontalStress / verticalStress2
-    ratio[mask] = np.nan
+    mask: npt.NDArray[ np.bool_ ] = np.abs( verticalStress ) < EPSILON
+    verticalStress2: npt.NDArray[ np.float64 ] = np.copy( verticalStress )
+    verticalStress2[ mask ] = 1.0
+    ratio: npt.NDArray[ np.float64 ] = horizontalStress / verticalStress2
+    ratio[ mask ] = np.nan
     return ratio
 
 
-def lithostaticStress(
-    depth: npt.NDArray[np.float64], density: npt.NDArray[np.float64], gravity: float
-) -> npt.NDArray[np.float64]:
+def lithostaticStress( depth: npt.NDArray[ np.float64 ], density: npt.NDArray[ np.float64 ],
+                       gravity: float ) -> npt.NDArray[ np.float64 ]:
     """Compute the lithostatic stress.
 
     Args:
@@ -349,19 +322,17 @@ def lithostaticStress(
     assert depth is not None, "Depth must be defined"
     assert density is not None, "Density must be defined"
 
-    assert depth.size == density.size, (
-        "Depth array "
-        + "and density array sizes (i.e., number of cells) must be equal."
-    )
+    assert depth.size == density.size, ( "Depth array " +
+                                         "and density array sizes (i.e., number of cells) must be equal." )
     # use -1*depth to agree with Geos convention (i.e., compression with negative stress)
     return -depth * density * gravity
 
 
 def elasticStrainFromBulkShear(
-    deltaEffectiveStress: npt.NDArray[np.float64],
-    bulkModulus: npt.NDArray[np.float64],
-    shearModulus: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
+    deltaEffectiveStress: npt.NDArray[ np.float64 ],
+    bulkModulus: npt.NDArray[ np.float64 ],
+    shearModulus: npt.NDArray[ np.float64 ],
+) -> npt.NDArray[ np.float64 ]:
     r"""Compute elastic strain from Bulk and Shear moduli.
 
     See documentation on https://dnicolasespinoza.github.io/node5.html.
@@ -392,46 +363,34 @@ def elasticStrainFromBulkShear(
         npt.NDArray[np.float64]: elastic strain (:math:`\epsilon`)
 
     """
-    assert (
-        deltaEffectiveStress is not None
-    ), "Effective stress variation must be defined"
+    assert ( deltaEffectiveStress is not None ), "Effective stress variation must be defined"
     assert bulkModulus is not None, "Bulk modulus must be defined"
     assert shearModulus is not None, "Shear modulus must be defined"
 
-    assert deltaEffectiveStress.shape[0] == bulkModulus.size, (
-        "Effective stress variation "
-        + " and bulk modulus array sizes (i.e., number of cells) must be equal."
-    )
+    assert deltaEffectiveStress.shape[ 0 ] == bulkModulus.size, (
+        "Effective stress variation " + " and bulk modulus array sizes (i.e., number of cells) must be equal." )
     assert shearModulus.size == bulkModulus.size, (
-        "Shear modulus "
-        + "and bulk modulus array sizes (i.e., number of cells) must be equal."
-    )
-    assert deltaEffectiveStress.shape[1] == 6, (
-        "Effective stress variation " + "number of components must be equal to 6."
-    )
-
-    elasticStrain: npt.NDArray[np.float64] = np.full_like(deltaEffectiveStress, np.nan)
-    for i, stressVector in enumerate(deltaEffectiveStress):
-        stiffnessTensor: npt.NDArray[np.float64] = shearModulus[i] * np.identity(
-            6, dtype=float
-        )
-        for k in range(3):
-            for m in range(3):
-                val: float = (
-                    (bulkModulus[i] + 4.0 / 3.0 * shearModulus[i])
-                    if k == m
-                    else (bulkModulus[i] - 2.0 / 3.0 * shearModulus[i])
-                )
-                stiffnessTensor[k, m] = val
-        elasticStrain[i] = stressVector @ np.linalg.inv(stiffnessTensor)
+        "Shear modulus " + "and bulk modulus array sizes (i.e., number of cells) must be equal." )
+    assert deltaEffectiveStress.shape[ 1 ] == 6, ( "Effective stress variation " +
+                                                   "number of components must be equal to 6." )
+
+    elasticStrain: npt.NDArray[ np.float64 ] = np.full_like( deltaEffectiveStress, np.nan )
+    for i, stressVector in enumerate( deltaEffectiveStress ):
+        stiffnessTensor: npt.NDArray[ np.float64 ] = shearModulus[ i ] * np.identity( 6, dtype=float )
+        for k in range( 3 ):
+            for m in range( 3 ):
+                val: float = ( ( bulkModulus[ i ] + 4.0 / 3.0 * shearModulus[ i ] ) if k == m else
+                               ( bulkModulus[ i ] - 2.0 / 3.0 * shearModulus[ i ] ) )
+                stiffnessTensor[ k, m ] = val
+        elasticStrain[ i ] = stressVector @ np.linalg.inv( stiffnessTensor )
     return elasticStrain
 
 
 def elasticStrainFromYoungPoisson(
-    deltaEffectiveStress: npt.NDArray[np.float64],
-    youngModulus: npt.NDArray[np.float64],
-    poissonRatio: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
+    deltaEffectiveStress: npt.NDArray[ np.float64 ],
+    youngModulus: npt.NDArray[ np.float64 ],
+    poissonRatio: npt.NDArray[ np.float64 ],
+) -> npt.NDArray[ np.float64 ]:
     r"""Compute elastic strain from Young modulus and Poisson ratio.
 
     See documentation on https://dnicolasespinoza.github.io/node5.html.
@@ -462,47 +421,34 @@ def elasticStrainFromYoungPoisson(
         npt.NDArray[np.float64]: elastic strain (:math:`\epsilon`)
 
     """
-    assert (
-        deltaEffectiveStress is not None
-    ), "Effective stress variation must be defined"
+    assert ( deltaEffectiveStress is not None ), "Effective stress variation must be defined"
     assert youngModulus is not None, "Bulk modulus must be defined"
     assert poissonRatio is not None, "Shear modulus must be defined"
 
-    assert deltaEffectiveStress.shape[0] == youngModulus.size, (
-        "Effective stress variation "
-        + " and bulk modulus array sizes (i.e., number of cells) must be equal."
-    )
+    assert deltaEffectiveStress.shape[ 0 ] == youngModulus.size, (
+        "Effective stress variation " + " and bulk modulus array sizes (i.e., number of cells) must be equal." )
     assert poissonRatio.size == youngModulus.size, (
-        "Shear modulus "
-        + "and bulk modulus array sizes (i.e., number of cells) must be equal."
-    )
-    assert deltaEffectiveStress.shape[1] == 6, (
-        "Effective stress variation " + "number of components must be equal to 6."
-    )
+        "Shear modulus " + "and bulk modulus array sizes (i.e., number of cells) must be equal." )
+    assert deltaEffectiveStress.shape[ 1 ] == 6, ( "Effective stress variation " +
+                                                   "number of components must be equal to 6." )
 
     # use of Lamé's equations
-    lambdaCoeff: npt.NDArray[np.float64] = lambdaCoefficient(youngModulus, poissonRatio)
-    shearMod: npt.NDArray[np.float64] = shearModulus(youngModulus, poissonRatio)
-
-    elasticStrain: npt.NDArray[np.float64] = np.full_like(deltaEffectiveStress, np.nan)
-    for i, deltaStressVector in enumerate(deltaEffectiveStress):
-        stiffnessTensor: npt.NDArray[np.float64] = shearMod[i] * np.identity(
-            6, dtype=float
-        )
-        for k in range(3):
-            for m in range(3):
-                val: float = (
-                    (lambdaCoeff[i] + 2.0 * shearMod[i]) if k == m else (lambdaCoeff[i])
-                )
-                stiffnessTensor[k, m] = val
-
-        elasticStrain[i] = deltaStressVector @ np.linalg.inv(stiffnessTensor)
+    lambdaCoeff: npt.NDArray[ np.float64 ] = lambdaCoefficient( youngModulus, poissonRatio )
+    shearMod: npt.NDArray[ np.float64 ] = shearModulus( youngModulus, poissonRatio )
+
+    elasticStrain: npt.NDArray[ np.float64 ] = np.full_like( deltaEffectiveStress, np.nan )
+    for i, deltaStressVector in enumerate( deltaEffectiveStress ):
+        stiffnessTensor: npt.NDArray[ np.float64 ] = shearMod[ i ] * np.identity( 6, dtype=float )
+        for k in range( 3 ):
+            for m in range( 3 ):
+                val: float = ( ( lambdaCoeff[ i ] + 2.0 * shearMod[ i ] ) if k == m else ( lambdaCoeff[ i ] ) )
+                stiffnessTensor[ k, m ] = val
+
+        elasticStrain[ i ] = deltaStressVector @ np.linalg.inv( stiffnessTensor )
     return elasticStrain
 
 
-def deviatoricStressPathOed(
-    poissonRatio: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
+def deviatoricStressPathOed( poissonRatio: npt.NDArray[ np.float64 ], ) -> npt.NDArray[ np.float64 ]:
     r"""Compute the Deviatoric Stress Path parameter in oedometric conditions.
 
     This parameter corresponds to the ratio between horizontal and vertical
@@ -522,17 +468,16 @@ def deviatoricStressPathOed(
     assert poissonRatio is not None, "Poisson's ratio must be defined"
 
     # manage division by 0 by replacing with nan
-    den: npt.NDArray[np.float64] = 1 - poissonRatio
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    ratio: npt.NDArray[np.float64] = poissonRatio / den
-    ratio[mask] = np.nan
+    den: npt.NDArray[ np.float64 ] = 1 - poissonRatio
+    mask: npt.NDArray[ np.bool_ ] = np.abs( den ) < EPSILON
+    den[ mask ] = 1.0
+    ratio: npt.NDArray[ np.float64 ] = poissonRatio / den
+    ratio[ mask ] = np.nan
     return ratio
 
 
-def reservoirStressPathReal(
-    deltaStress: npt.NDArray[np.float64], deltaPressure: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def reservoirStressPathReal( deltaStress: npt.NDArray[ np.float64 ],
+                             deltaPressure: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute real reservoir stress path.
 
     .. math::
@@ -551,27 +496,24 @@ def reservoirStressPathReal(
     assert deltaPressure is not None, "Pressure deviation must be defined"
     assert deltaStress is not None, "Stress deviation must be defined"
 
-    assert deltaStress.shape[0] == deltaPressure.size, (
-        "Total stress array and pressure variation "
-        + "array sizes (i.e., number of cells) must be equal."
-    )
+    assert deltaStress.shape[ 0 ] == deltaPressure.size, ( "Total stress array and pressure variation " +
+                                                           "array sizes (i.e., number of cells) must be equal." )
 
     # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(deltaPressure) < EPSILON
-    den: npt.NDArray[np.float64] = np.copy(deltaPressure)
-    den[mask] = 1.0
+    mask: npt.NDArray[ np.bool_ ] = np.abs( deltaPressure ) < EPSILON
+    den: npt.NDArray[ np.float64 ] = np.copy( deltaPressure )
+    den[ mask ] = 1.0
     # use -1 to agree with Geos convention (i.e., compression with negative stress)
     # take the xx, yy, and zz components only
-    rsp: npt.NDArray[np.float64] = np.copy(deltaStress[:, :3])
-    for j in range(rsp.shape[1]):
-        rsp[:, j] /= den
-        rsp[mask, j] = np.nan
+    rsp: npt.NDArray[ np.float64 ] = np.copy( deltaStress[ :, :3 ] )
+    for j in range( rsp.shape[ 1 ] ):
+        rsp[ :, j ] /= den
+        rsp[ mask, j ] = np.nan
     return rsp
 
 
-def reservoirStressPathOed(
-    biotCoefficient: npt.NDArray[np.float64], poissonRatio: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def reservoirStressPathOed( biotCoefficient: npt.NDArray[ np.float64 ],
+                            poissonRatio: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute reservoir stress path in oedometric conditions.
 
     .. math::
@@ -589,22 +531,19 @@ def reservoirStressPathOed(
     assert poissonRatio is not None, "Poisson's ratio must be defined"
 
     assert biotCoefficient.size == poissonRatio.size, (
-        "Biot coefficient array and "
-        + "Poisson's ratio array sizes (i.e., number of cells) must be equal."
-    )
+        "Biot coefficient array and " + "Poisson's ratio array sizes (i.e., number of cells) must be equal." )
 
     # manage division by 0 by replacing with nan
-    den: npt.NDArray[np.float64] = 1.0 - poissonRatio
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
-    rsp: npt.NDArray[np.float64] = biotCoefficient * (1.0 - 2.0 * poissonRatio) / den
-    rsp[mask] = np.nan
+    den: npt.NDArray[ np.float64 ] = 1.0 - poissonRatio
+    mask: npt.NDArray[ np.bool_ ] = np.abs( den ) < EPSILON
+    den[ mask ] = 1.0
+    rsp: npt.NDArray[ np.float64 ] = biotCoefficient * ( 1.0 - 2.0 * poissonRatio ) / den
+    rsp[ mask ] = np.nan
     return rsp
 
 
-def criticalTotalStressRatio(
-    pressure: npt.NDArray[np.float64], verticalStress: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def criticalTotalStressRatio( pressure: npt.NDArray[ np.float64 ],
+                              verticalStress: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute critical total stress ratio.
 
     Corresponds to the fracture index from Lemgruber-Traby et al (2024).
@@ -631,22 +570,19 @@ def criticalTotalStressRatio(
     assert verticalStress is not None, "Vertical stress must be defined"
 
     assert pressure.size == verticalStress.size, (
-        "pressure array and "
-        + "vertical stress array sizes (i.e., number of cells) must be equal."
-    )
+        "pressure array and " + "vertical stress array sizes (i.e., number of cells) must be equal." )
     # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(verticalStress) < EPSILON
+    mask: npt.NDArray[ np.bool_ ] = np.abs( verticalStress ) < EPSILON
     # use -1 to agree with Geos convention (i.e., compression with negative stress)
-    verticalStress2: npt.NDArray[np.float64] = -1.0 * np.copy(verticalStress)
-    verticalStress2[mask] = 1.0
-    fi: npt.NDArray[np.float64] = pressure / verticalStress2
-    fi[mask] = np.nan
+    verticalStress2: npt.NDArray[ np.float64 ] = -1.0 * np.copy( verticalStress )
+    verticalStress2[ mask ] = 1.0
+    fi: npt.NDArray[ np.float64 ] = pressure / verticalStress2
+    fi[ mask ] = np.nan
     return fi
 
 
-def totalStressRatioThreshold(
-    pressure: npt.NDArray[np.float64], horizontalStress: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def totalStressRatioThreshold( pressure: npt.NDArray[ np.float64 ],
+                               horizontalStress: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute total stress ratio threshold.
 
     Corresponds to the fracture threshold from Lemgruber-Traby et al (2024).
@@ -673,25 +609,23 @@ def totalStressRatioThreshold(
     assert horizontalStress is not None, "Horizontal stress must be defined"
 
     assert pressure.size == horizontalStress.size, (
-        "pressure array and "
-        + "horizontal stress array sizes (i.e., number of cells) must be equal."
-    )
+        "pressure array and " + "horizontal stress array sizes (i.e., number of cells) must be equal." )
 
     # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(horizontalStress) < EPSILON
+    mask: npt.NDArray[ np.bool_ ] = np.abs( horizontalStress ) < EPSILON
     # use -1 to agree with Geos convention (i.e., compression with negative stress)
-    horizontalStress2: npt.NDArray[np.float64] = -1.0 * np.copy(horizontalStress)
-    horizontalStress2[mask] = 1.0
-    ft: npt.NDArray[np.float64] = pressure / horizontalStress2
-    ft[mask] = np.nan
+    horizontalStress2: npt.NDArray[ np.float64 ] = -1.0 * np.copy( horizontalStress )
+    horizontalStress2[ mask ] = 1.0
+    ft: npt.NDArray[ np.float64 ] = pressure / horizontalStress2
+    ft[ mask ] = np.nan
     return ft
 
 
 def criticalPorePressure(
-    stressVector: npt.NDArray[np.float64],
+    stressVector: npt.NDArray[ np.float64 ],
     rockCohesion: float,
     frictionAngle: float = 0.0,
-) -> npt.NDArray[np.float64]:
+) -> npt.NDArray[ np.float64 ]:
     r"""Compute the critical pore pressure.
 
     Fracturing can occur in areas where Critical pore pressure is greater than
@@ -717,34 +651,27 @@ def criticalPorePressure(
 
     """
     assert stressVector is not None, "Stress vector must be defined"
-    assert stressVector.shape[1] == 6, "Stress vector must be of size 6."
+    assert stressVector.shape[ 1 ] == 6, "Stress vector must be of size 6."
 
-    assert (frictionAngle >= 0.0) and (frictionAngle < np.pi / 2.0), (
-        "Fristion angle " + "must range between 0 and pi/2."
-    )
+    assert ( frictionAngle >= 0.0 ) and ( frictionAngle < np.pi / 2.0 ), ( "Fristion angle " +
+                                                                           "must range between 0 and pi/2." )
 
-    minimumPrincipalStress: npt.NDArray[np.float64] = np.full(
-        stressVector.shape[0], np.nan
-    )
-    maximumPrincipalStress: npt.NDArray[np.float64] = np.copy(minimumPrincipalStress)
-    for i in range(minimumPrincipalStress.shape[0]):
-        p3, p2, p1 = computeStressPrincipalComponentsFromStressVector(stressVector[i])
-        minimumPrincipalStress[i] = p3
-        maximumPrincipalStress[i] = p1
+    minimumPrincipalStress: npt.NDArray[ np.float64 ] = np.full( stressVector.shape[ 0 ], np.nan )
+    maximumPrincipalStress: npt.NDArray[ np.float64 ] = np.copy( minimumPrincipalStress )
+    for i in range( minimumPrincipalStress.shape[ 0 ] ):
+        p3, p2, p1 = computeStressPrincipalComponentsFromStressVector( stressVector[ i ] )
+        minimumPrincipalStress[ i ] = p3
+        maximumPrincipalStress[ i ] = p1
 
     # assertion frictionAngle < np.pi/2., so sin(frictionAngle) != 1
-    cohesiveTerm: npt.NDArray[np.float64] = (
-        rockCohesion * np.cos(frictionAngle) / (1 - np.sin(frictionAngle))
-    )
-    residualTerm: npt.NDArray[np.float64] = (
-        3 * minimumPrincipalStress - maximumPrincipalStress
-    ) / 2.0
+    cohesiveTerm: npt.NDArray[ np.float64 ] = ( rockCohesion * np.cos( frictionAngle ) /
+                                                ( 1 - np.sin( frictionAngle ) ) )
+    residualTerm: npt.NDArray[ np.float64 ] = ( 3 * minimumPrincipalStress - maximumPrincipalStress ) / 2.0
     return cohesiveTerm + residualTerm
 
 
-def criticalPorePressureThreshold(
-    pressure: npt.NDArray[np.float64], criticalPorePressure: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
+def criticalPorePressureThreshold( pressure: npt.NDArray[ np.float64 ],
+                                   criticalPorePressure: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute the critical pore pressure threshold.
 
     Defined as the ratio between pressure and critical pore pressure.
@@ -766,24 +693,22 @@ def criticalPorePressureThreshold(
     assert criticalPorePressure is not None, "Critical pore pressure must be defined"
 
     assert pressure.size == criticalPorePressure.size, (
-        "Pressure array and critical"
-        + " pore pressure array sizes (i.e., number of cells) must be equal."
-    )
+        "Pressure array and critical" + " pore pressure array sizes (i.e., number of cells) must be equal." )
 
     # manage division by 0 by replacing with nan
-    mask: npt.NDArray[np.bool_] = np.abs(criticalPorePressure) < EPSILON
-    den = np.copy(criticalPorePressure)
-    den[mask] = 1.0
-    index: npt.NDArray[np.float64] = pressure / den
-    index[mask] = np.nan
+    mask: npt.NDArray[ np.bool_ ] = np.abs( criticalPorePressure ) < EPSILON
+    den = np.copy( criticalPorePressure )
+    den[ mask ] = 1.0
+    index: npt.NDArray[ np.float64 ] = pressure / den
+    index[ mask ] = np.nan
     return index
 
 
 def compressibilityOed(
-    shearModulus: npt.NDArray[np.float64],
-    bulkModulus: npt.NDArray[np.float64],
-    porosity: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
+    shearModulus: npt.NDArray[ np.float64 ],
+    bulkModulus: npt.NDArray[ np.float64 ],
+    porosity: npt.NDArray[ np.float64 ],
+) -> npt.NDArray[ np.float64 ]:
     r"""Compute compressibility from elastic moduli and porosity.
 
     Compressibility formula is:
@@ -803,25 +728,25 @@ def compressibilityOed(
     assert bulkModulus is not None, "Bulk modulus must be defined"
     assert shearModulus is not None, "Shear modulus must be defined"
     assert porosity is not None, "Porosity must be defined"
-    mask1: npt.NDArray[np.bool_] = np.abs(porosity) < EPSILON
-    den1: npt.NDArray[np.float64] = np.copy(porosity)
-    den1[mask1] = 1.0
+    mask1: npt.NDArray[ np.bool_ ] = np.abs( porosity ) < EPSILON
+    den1: npt.NDArray[ np.float64 ] = np.copy( porosity )
+    den1[ mask1 ] = 1.0
 
-    den2: npt.NDArray[np.float64] = 3.0 * bulkModulus + 4.0 * shearModulus
-    mask2: npt.NDArray[np.bool_] = np.abs(den2) < EPSILON
-    den2[mask2] = 1.0
+    den2: npt.NDArray[ np.float64 ] = 3.0 * bulkModulus + 4.0 * shearModulus
+    mask2: npt.NDArray[ np.bool_ ] = np.abs( den2 ) < EPSILON
+    den2[ mask2 ] = 1.0
 
-    comprOed: npt.NDArray[np.float64] = 1.0 / den1 * 3.0 / den2
-    comprOed[mask1] = np.nan
-    comprOed[mask2] = np.nan
+    comprOed: npt.NDArray[ np.float64 ] = 1.0 / den1 * 3.0 / den2
+    comprOed[ mask1 ] = np.nan
+    comprOed[ mask2 ] = np.nan
     return comprOed
 
 
 def compressibilityReal(
-    deltaPressure: npt.NDArray[np.float64],
-    porosity: npt.NDArray[np.float64],
-    porosityInitial: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
+    deltaPressure: npt.NDArray[ np.float64 ],
+    porosity: npt.NDArray[ np.float64 ],
+    porosityInitial: npt.NDArray[ np.float64 ],
+) -> npt.NDArray[ np.float64 ]:
     r"""Compute compressibility from elastic moduli and porosity.
 
     Compressibility formula is:
@@ -844,21 +769,21 @@ def compressibilityReal(
     assert porosity is not None, "Porosity must be defined"
     assert porosityInitial is not None, "Initial porosity must be defined"
 
-    den: npt.NDArray[np.float64] = deltaPressure * porosityInitial
-    mask: npt.NDArray[np.bool_] = np.abs(den) < EPSILON
-    den[mask] = 1.0
+    den: npt.NDArray[ np.float64 ] = deltaPressure * porosityInitial
+    mask: npt.NDArray[ np.bool_ ] = np.abs( den ) < EPSILON
+    den[ mask ] = 1.0
 
-    comprReal: npt.NDArray[np.float64] = (porosity - porosityInitial) / den
-    comprReal[mask] = np.nan
+    comprReal: npt.NDArray[ np.float64 ] = ( porosity - porosityInitial ) / den
+    comprReal[ mask ] = np.nan
     return comprReal
 
 
 def compressibility(
-    poissonRatio: npt.NDArray[np.float64],
-    bulkModulus: npt.NDArray[np.float64],
-    biotCoefficient: npt.NDArray[np.float64],
-    porosity: npt.NDArray[np.float64],
-) -> npt.NDArray[np.float64]:
+    poissonRatio: npt.NDArray[ np.float64 ],
+    bulkModulus: npt.NDArray[ np.float64 ],
+    biotCoefficient: npt.NDArray[ np.float64 ],
+    porosity: npt.NDArray[ np.float64 ],
+) -> npt.NDArray[ np.float64 ]:
     r"""Compute compressibility from elastic moduli, biot coefficient and porosity.
 
     Compressibility formula is:
@@ -881,34 +806,27 @@ def compressibility(
     assert biotCoefficient is not None, "Biot coefficient must be defined"
     assert porosity is not None, "Porosity must be defined"
 
-    term1: npt.NDArray[np.float64] = 1.0 - 2.0 * poissonRatio
+    term1: npt.NDArray[ np.float64 ] = 1.0 - 2.0 * poissonRatio
 
-    mask: npt.NDArray[np.bool_] = (np.abs(bulkModulus) < EPSILON) * (
-        np.abs(porosity) < EPSILON
-    )
-    denFac1: npt.NDArray[np.float64] = porosity * bulkModulus
-    term1[mask] = 1.0
+    mask: npt.NDArray[ np.bool_ ] = ( np.abs( bulkModulus ) < EPSILON ) * ( np.abs( porosity ) < EPSILON )
+    denFac1: npt.NDArray[ np.float64 ] = porosity * bulkModulus
+    term1[ mask ] = 1.0
     term1 /= denFac1
 
-    term2M1: npt.NDArray[np.float64] = (
-        biotCoefficient * biotCoefficient * (1 + poissonRatio)
-    )
-    denTerm2M1: npt.NDArray[np.float64] = 1 - poissonRatio
-    mask2: npt.NDArray[np.bool_] = np.abs(denTerm2M1) < EPSILON
-    denTerm2M1[mask2] = 1.0
+    term2M1: npt.NDArray[ np.float64 ] = ( biotCoefficient * biotCoefficient * ( 1 + poissonRatio ) )
+    denTerm2M1: npt.NDArray[ np.float64 ] = 1 - poissonRatio
+    mask2: npt.NDArray[ np.bool_ ] = np.abs( denTerm2M1 ) < EPSILON
+    denTerm2M1[ mask2 ] = 1.0
     term2M1 /= denTerm2M1
-    term2M1[mask2] = np.nan
+    term2M1[ mask2 ] = np.nan
 
-    term2M2: npt.NDArray[np.float64] = (
-        3.0 * (biotCoefficient - porosity) * (1 - biotCoefficient)
-    )
-    term2: npt.NDArray[np.float64] = term2M1 + term2M2
+    term2M2: npt.NDArray[ np.float64 ] = ( 3.0 * ( biotCoefficient - porosity ) * ( 1 - biotCoefficient ) )
+    term2: npt.NDArray[ np.float64 ] = term2M1 + term2M2
     return term1 * term2
 
 
-def shearCapacityUtilization(
-    traction: npt.NDArray[np.float64], rockCohesion: float, frictionAngle: float
-) -> npt.NDArray[np.float64]:
+def shearCapacityUtilization( traction: npt.NDArray[ np.float64 ], rockCohesion: float,
+                              frictionAngle: float ) -> npt.NDArray[ np.float64 ]:
     r"""Compute shear capacity utilization (SCU).
 
     .. math::
@@ -927,27 +845,26 @@ def shearCapacityUtilization(
 
     """
     assert traction is not None, "Traction must be defined"
-    assert traction.shape[1] == 3, "Traction vector must have 3 components."
+    assert traction.shape[ 1 ] == 3, "Traction vector must have 3 components."
 
-    scu: npt.NDArray[np.float64] = np.full(traction.shape[0], np.nan)
-    for i, tractionVec in enumerate(traction):
+    scu: npt.NDArray[ np.float64 ] = np.full( traction.shape[ 0 ], np.nan )
+    for i, tractionVec in enumerate( traction ):
         # use -1 to agree with Geos convention (i.e., compression with negative stress)
-        stressNormal: npt.NDArray[np.float64] = -1.0 * tractionVec[0]
+        stressNormal: npt.NDArray[ np.float64 ] = -1.0 * tractionVec[ 0 ]
 
         # compute failure envelope
-        mohrCoulomb: MohrCoulomb = MohrCoulomb(rockCohesion, frictionAngle)
-        tauFailure: float = float(mohrCoulomb.computeShearStress(stressNormal))
+        mohrCoulomb: MohrCoulomb = MohrCoulomb( rockCohesion, frictionAngle )
+        tauFailure: float = float( mohrCoulomb.computeShearStress( stressNormal ) )
         scu_i: float = np.nan
         if tauFailure > 0:
-            scu_i = np.abs(tractionVec[1]) / tauFailure
+            scu_i = np.abs( tractionVec[ 1 ] ) / tauFailure
         # compute SCU
-        scu[i] = scu_i
+        scu[ i ] = scu_i
     return scu
 
 
 def computeStressPrincipalComponentsFromStressVector(
-    stressVector: npt.NDArray[np.float64],
-) -> tuple[float, float, float]:
+    stressVector: npt.NDArray[ np.float64 ], ) -> tuple[ float, float, float ]:
     """Compute stress principal components from stress vector.
 
     Args:
@@ -959,13 +876,11 @@ def computeStressPrincipalComponentsFromStressVector(
 
     """
     assert stressVector.size == 6, "Stress vector dimension is wrong."
-    stressTensor: npt.NDArray[np.float64] = getAttributeMatrixFromVector(stressVector)
-    return computeStressPrincipalComponents(stressTensor)
+    stressTensor: npt.NDArray[ np.float64 ] = getAttributeMatrixFromVector( stressVector )
+    return computeStressPrincipalComponents( stressTensor )
 
 
-def computeStressPrincipalComponents(
-    stressTensor: npt.NDArray[np.float64],
-) -> tuple[float, float, float]:
+def computeStressPrincipalComponents( stressTensor: npt.NDArray[ np.float64 ], ) -> tuple[ float, float, float ]:
     """Compute stress principal components.
 
     Args:
@@ -977,15 +892,14 @@ def computeStressPrincipalComponents(
 
     """
     # get eigen values
-    e_val, e_vec = np.linalg.eig(stressTensor)
+    e_val, e_vec = np.linalg.eig( stressTensor )
     # sort principal stresses from smallest to largest
-    p3, p2, p1 = np.sort(e_val)
-    return (p3, p2, p1)
+    p3, p2, p1 = np.sort( e_val )
+    return ( p3, p2, p1 )
 
 
-def computeNormalShearStress(
-    stressTensor: npt.NDArray[np.float64], directionVector: npt.NDArray[np.float64]
-) -> tuple[float, float]:
+def computeNormalShearStress( stressTensor: npt.NDArray[ np.float64 ],
+                              directionVector: npt.NDArray[ np.float64 ] ) -> tuple[ float, float ]:
     """Compute normal and shear stress according to stress tensor and direction.
 
     Args:
@@ -996,16 +910,16 @@ def computeNormalShearStress(
         tuple[float, float]: normal and shear stresses.
 
     """
-    assert stressTensor.shape == (3, 3), "Stress tensor must be 3x3 matrix."
+    assert stressTensor.shape == ( 3, 3 ), "Stress tensor must be 3x3 matrix."
     assert directionVector.size == 3, "Direction vector must have 3 components"
 
     # normalization of direction vector
-    directionVector = directionVector / np.linalg.norm(directionVector)
+    directionVector = directionVector / np.linalg.norm( directionVector )
     # stress vector
-    T: npt.NDArray[np.float64] = np.dot(stressTensor, directionVector)
+    T: npt.NDArray[ np.float64 ] = np.dot( stressTensor, directionVector )
     # normal stress
-    sigmaN: float = np.dot(T, directionVector)
+    sigmaN: float = np.dot( T, directionVector )
     # shear stress
-    tauVec: npt.NDArray[np.float64] = T - np.dot(sigmaN, directionVector)
-    tau: float = float(np.linalg.norm(tauVec))
-    return (sigmaN, tau)
\ No newline at end of file
+    tauVec: npt.NDArray[ np.float64 ] = T - np.dot( sigmaN, directionVector )
+    tau: float = float( np.linalg.norm( tauVec ) )
+    return ( sigmaN, tau )
diff --git a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
index 827cc815..82b90bff 100644
--- a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
+++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
@@ -12,369 +12,267 @@
 import pandas as pd  # type: ignore[import-untyped]
 from typing_extensions import Self
 
-dir_path = os.path.dirname(os.path.realpath(__file__))
-parent_dir_path = os.path.join(os.path.dirname(dir_path), "src")
+dir_path = os.path.dirname( os.path.realpath( __file__ ) )
+parent_dir_path = os.path.join( os.path.dirname( dir_path ), "src" )
 if parent_dir_path not in sys.path:
-    sys.path.append(parent_dir_path)
+    sys.path.append( parent_dir_path )
 
 import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 from geos.utils import PhysicalConstants
 from geos.utils.algebraFunctions import getAttributeMatrixFromVector
 
 # geomechanical outputs - Testing variables - Unit is GPa
-bulkModulus: npt.NDArray[np.float64] = np.array([9.0, 50.0, 65.0, np.nan, 150.0])
-shearModulus: npt.NDArray[np.float64] = np.array([3.2, 24.0, np.nan, 70.0, 100.0])
-poissonRatio: npt.NDArray[np.float64] = np.array([0.34, 0.29, np.nan, np.nan, 0.23])
-density: npt.NDArray[np.float64] = np.array([2600.0, 3100.0, 2800.0, np.nan, 2900.0])
-biotCoefficient: npt.NDArray[np.float64] = np.array([0.1, 0.3, np.nan, 0.5, 0.6])
-pressure: npt.NDArray[np.float64] = np.array([10.0, 280.0, np.nan, 80.0, 100.0])
-effectiveStress: npt.NDArray[np.float64] = -1.0 * np.array(
-    [
-        [160.0, 250.0, 180.0, np.nan, 200.0, 190.0],
-        [180.0, 260.0, 190.0, np.nan, 220.0, 210.0],
-        [200.0, 280.0, 210.0, np.nan, 230.0, 230.0],
-        [220.0, 290.0, 220.0, np.nan, 250.0, 250.0],
-        [260.0, 310.0, 230.0, np.nan, 290.0, 270.0],
-    ]
-)
-verticalStress: npt.NDArray[np.float64] = effectiveStress[:, 2]
-horizontalStress: npt.NDArray[np.float64] = np.min(effectiveStress[:, :2], axis=1)
-
-stressVector: npt.NDArray[np.float64] = np.array([1.8, 2.5, 5.2, 0.3, 0.4, 0.1])
-stressTensor: npt.NDArray[np.float64] = getAttributeMatrixFromVector(stressVector)
-
-
-class TestsFunctionsGeomechanicsCalculator(unittest.TestCase):
-    def test_SpecificGravity(self: Self) -> None:
+bulkModulus: npt.NDArray[ np.float64 ] = np.array( [ 9.0, 50.0, 65.0, np.nan, 150.0 ] )
+shearModulus: npt.NDArray[ np.float64 ] = np.array( [ 3.2, 24.0, np.nan, 70.0, 100.0 ] )
+poissonRatio: npt.NDArray[ np.float64 ] = np.array( [ 0.34, 0.29, np.nan, np.nan, 0.23 ] )
+density: npt.NDArray[ np.float64 ] = np.array( [ 2600.0, 3100.0, 2800.0, np.nan, 2900.0 ] )
+biotCoefficient: npt.NDArray[ np.float64 ] = np.array( [ 0.1, 0.3, np.nan, 0.5, 0.6 ] )
+pressure: npt.NDArray[ np.float64 ] = np.array( [ 10.0, 280.0, np.nan, 80.0, 100.0 ] )
+effectiveStress: npt.NDArray[ np.float64 ] = -1.0 * np.array( [
+    [ 160.0, 250.0, 180.0, np.nan, 200.0, 190.0 ],
+    [ 180.0, 260.0, 190.0, np.nan, 220.0, 210.0 ],
+    [ 200.0, 280.0, 210.0, np.nan, 230.0, 230.0 ],
+    [ 220.0, 290.0, 220.0, np.nan, 250.0, 250.0 ],
+    [ 260.0, 310.0, 230.0, np.nan, 290.0, 270.0 ],
+] )
+verticalStress: npt.NDArray[ np.float64 ] = effectiveStress[ :, 2 ]
+horizontalStress: npt.NDArray[ np.float64 ] = np.min( effectiveStress[ :, :2 ], axis=1 )
+
+stressVector: npt.NDArray[ np.float64 ] = np.array( [ 1.8, 2.5, 5.2, 0.3, 0.4, 0.1 ] )
+stressTensor: npt.NDArray[ np.float64 ] = getAttributeMatrixFromVector( stressVector )
+
+
+class TestsFunctionsGeomechanicsCalculator( unittest.TestCase ):
+
+    def test_SpecificGravity( self: Self ) -> None:
         """Test calculation of specific Gravity."""
         specificDensity: float = 1000.0
 
-        obtained: npt.NDArray[np.float64] = fcts.specificGravity(
-            density, specificDensity
-        )
-        expected: npt.NDArray[np.float64] = np.array([2.6, 3.1, 2.8, np.nan, 2.9])
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
+        obtained: npt.NDArray[ np.float64 ] = fcts.specificGravity( density, specificDensity )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 2.6, 3.1, 2.8, np.nan, 2.9 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
 
-    def test_YoungModulus(self: Self) -> None:
+    def test_YoungModulus( self: Self ) -> None:
         """Test calculation of young Modulus."""
-        obtained: npt.NDArray[np.float64] = fcts.youngModulus(bulkModulus, shearModulus)
-        expected: npt.NDArray[np.float64] = np.array(
-            [8.58, 62.07, np.nan, np.nan, 245.45]
-        )
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_PoissonRatio(self: Self) -> None:
+        obtained: npt.NDArray[ np.float64 ] = fcts.youngModulus( bulkModulus, shearModulus )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 8.58, 62.07, np.nan, np.nan, 245.45 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_PoissonRatio( self: Self ) -> None:
         """Test calculation of poisson Ratio."""
-        obtained: npt.NDArray[np.float64] = fcts.poissonRatio(bulkModulus, shearModulus)
-        expected: npt.NDArray[np.float64] = np.array([0.34, 0.29, np.nan, np.nan, 0.23])
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
+        obtained: npt.NDArray[ np.float64 ] = fcts.poissonRatio( bulkModulus, shearModulus )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 0.34, 0.29, np.nan, np.nan, 0.23 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
 
-    def test_BulkModulus(self: Self) -> None:
+    def test_BulkModulus( self: Self ) -> None:
         """Test calculation of bulk Modulus."""
-        E: npt.NDArray[np.float64] = fcts.youngModulus(bulkModulus, shearModulus)
-        u: npt.NDArray[np.float64] = fcts.poissonRatio(bulkModulus, shearModulus)
+        E: npt.NDArray[ np.float64 ] = fcts.youngModulus( bulkModulus, shearModulus )
+        u: npt.NDArray[ np.float64 ] = fcts.poissonRatio( bulkModulus, shearModulus )
 
-        obtained: npt.NDArray[np.float64] = fcts.bulkModulus(E, u)
-        expected: npt.NDArray[np.float64] = bulkModulus
-        expected[2] = np.nan
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
+        obtained: npt.NDArray[ np.float64 ] = fcts.bulkModulus( E, u )
+        expected: npt.NDArray[ np.float64 ] = bulkModulus
+        expected[ 2 ] = np.nan
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
 
-    def test_ShearModulus(self: Self) -> None:
+    def test_ShearModulus( self: Self ) -> None:
         """Test calculation of shear Modulus."""
-        E: npt.NDArray[np.float64] = fcts.youngModulus(bulkModulus, shearModulus)
-        u: npt.NDArray[np.float64] = fcts.poissonRatio(bulkModulus, shearModulus)
+        E: npt.NDArray[ np.float64 ] = fcts.youngModulus( bulkModulus, shearModulus )
+        u: npt.NDArray[ np.float64 ] = fcts.poissonRatio( bulkModulus, shearModulus )
 
-        obtained: npt.NDArray[np.float64] = fcts.shearModulus(E, u)
-        expected: npt.NDArray[np.float64] = shearModulus
-        expected[3] = np.nan
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
+        obtained: npt.NDArray[ np.float64 ] = fcts.shearModulus( E, u )
+        expected: npt.NDArray[ np.float64 ] = shearModulus
+        expected[ 3 ] = np.nan
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
 
-    def test_lambdaCoefficient(self: Self) -> None:
+    def test_lambdaCoefficient( self: Self ) -> None:
         """Test calculation of shear Modulus."""
-        E: npt.NDArray[np.float64] = fcts.youngModulus(bulkModulus, shearModulus)
-        u: npt.NDArray[np.float64] = fcts.poissonRatio(bulkModulus, shearModulus)
-
-        obtained: npt.NDArray[np.float64] = fcts.lambdaCoefficient(E, u)
-        expected: npt.NDArray[np.float64] = np.array(
-            [6.87, 34.0, np.nan, np.nan, 83.33]
-        )
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_OedometricModulus(self: Self) -> None:
+        E: npt.NDArray[ np.float64 ] = fcts.youngModulus( bulkModulus, shearModulus )
+        u: npt.NDArray[ np.float64 ] = fcts.poissonRatio( bulkModulus, shearModulus )
+
+        obtained: npt.NDArray[ np.float64 ] = fcts.lambdaCoefficient( E, u )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 6.87, 34.0, np.nan, np.nan, 83.33 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_OedometricModulus( self: Self ) -> None:
         """Test calculation of oedometric Modulus."""
-        Edef: npt.NDArray[np.float64] = np.array([9.0, 50.0, 65.0, np.nan, 150.0])
+        Edef: npt.NDArray[ np.float64 ] = np.array( [ 9.0, 50.0, 65.0, np.nan, 150.0 ] )
 
-        obtained: npt.NDArray[np.float64] = fcts.oedometricModulus(Edef, poissonRatio)
-        expected: npt.NDArray[np.float64] = np.array(
-            [13.85, 65.52, np.nan, np.nan, 173.89]
-        )
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
+        obtained: npt.NDArray[ np.float64 ] = fcts.oedometricModulus( Edef, poissonRatio )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 13.85, 65.52, np.nan, np.nan, 173.89 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
 
-    def test_BiotCoefficient(self: Self) -> None:
+    def test_BiotCoefficient( self: Self ) -> None:
         """Test calculation of biot Coefficient."""
         Kgrain: float = 150.0
-        obtained: npt.NDArray[np.float64] = fcts.biotCoefficient(Kgrain, bulkModulus)
-        expected: npt.NDArray[np.float64] = np.array([0.94, 0.67, 0.57, np.nan, 0.0])
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
+        obtained: npt.NDArray[ np.float64 ] = fcts.biotCoefficient( Kgrain, bulkModulus )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 0.94, 0.67, 0.57, np.nan, 0.0 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
 
-    def test_TotalStress(self: Self) -> None:
+    def test_TotalStress( self: Self ) -> None:
         """Test calculation of Total Stress."""
-        obtained: npt.NDArray[np.float64] = fcts.totalStress(
-            effectiveStress, biotCoefficient, pressure
-        )
-        expected: npt.NDArray[np.float64] = np.array(
-            [
-                [-161.0, -251.0, -181.0, np.nan, -200.0, -190.0],
-                [-264.0, -344.0, -274.0, np.nan, -220.0, -210.0],
-                [np.nan, np.nan, np.nan, np.nan, -230.0, -230.0],
-                [-260.0, -330.0, -260.0, np.nan, -250.0, -250.0],
-                [-320.0, -370.0, -290.0, np.nan, -290.0, -270.0],
-            ]
-        )
-
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_StressRatio(self: Self) -> None:
+        obtained: npt.NDArray[ np.float64 ] = fcts.totalStress( effectiveStress, biotCoefficient, pressure )
+        expected: npt.NDArray[ np.float64 ] = np.array( [
+            [ -161.0, -251.0, -181.0, np.nan, -200.0, -190.0 ],
+            [ -264.0, -344.0, -274.0, np.nan, -220.0, -210.0 ],
+            [ np.nan, np.nan, np.nan, np.nan, -230.0, -230.0 ],
+            [ -260.0, -330.0, -260.0, np.nan, -250.0, -250.0 ],
+            [ -320.0, -370.0, -290.0, np.nan, -290.0, -270.0 ],
+        ] )
+
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_StressRatio( self: Self ) -> None:
         """Test calculation of Stress Ratio."""
-        obtained: npt.NDArray[np.float64] = fcts.stressRatio(
-            horizontalStress, verticalStress
-        )
-        expected: npt.NDArray[np.float64] = np.array([1.39, 1.37, 1.33, 1.32, 1.35])
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_LithostaticStress(self: Self) -> None:
+        obtained: npt.NDArray[ np.float64 ] = fcts.stressRatio( horizontalStress, verticalStress )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 1.39, 1.37, 1.33, 1.32, 1.35 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_LithostaticStress( self: Self ) -> None:
         """Test calculation of lithostatic Stress."""
-        zCoords: npt.NDArray[np.float64] = np.array(
-            [-2.0, -100.0, -500.0, np.nan, -2500.0]
-        )
+        zCoords: npt.NDArray[ np.float64 ] = np.array( [ -2.0, -100.0, -500.0, np.nan, -2500.0 ] )
         gravity: float = PhysicalConstants.GRAVITY
 
-        obtained: npt.NDArray[np.float64] = fcts.lithostaticStress(
-            zCoords, density, gravity
-        )
-        expected: npt.NDArray[np.float64] = np.array(
-            [51012.0, 3041100.0, 13734000.0, np.nan, 71122500.0]
-        )
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_ElasticStrainFromBulkShear(self: Self) -> None:
+        obtained: npt.NDArray[ np.float64 ] = fcts.lithostaticStress( zCoords, density, gravity )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 51012.0, 3041100.0, 13734000.0, np.nan, 71122500.0 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_ElasticStrainFromBulkShear( self: Self ) -> None:
         """Test calculation of Elastic Strain."""
-        deltaEffectiveStress: npt.NDArray[np.float64] = -1.3 * effectiveStress
-        deltaEffectiveStress[np.isnan(deltaEffectiveStress)] = 150.0
-        obtained: npt.NDArray[np.float64] = fcts.elasticStrainFromBulkShear(
-            deltaEffectiveStress, bulkModulus, shearModulus
-        )
-        expected: npt.NDArray[np.float64] = np.array(
-            [
-                [2.02, 20.3, 6.08, 46.88, 81.25, 77.19],
-                [1.01, 3.17, 1.28, 6.25, 11.92, 11.38],
-                [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-                [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-                [0.73, 1.05, 0.53, 1.5, 3.77, 3.51],
-            ]
-        )
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_ElasticStrainFromYoungPoisson(self: Self) -> None:
+        deltaEffectiveStress: npt.NDArray[ np.float64 ] = -1.3 * effectiveStress
+        deltaEffectiveStress[ np.isnan( deltaEffectiveStress ) ] = 150.0
+        obtained: npt.NDArray[ np.float64 ] = fcts.elasticStrainFromBulkShear( deltaEffectiveStress, bulkModulus,
+                                                                               shearModulus )
+        expected: npt.NDArray[ np.float64 ] = np.array( [
+            [ 2.02, 20.3, 6.08, 46.88, 81.25, 77.19 ],
+            [ 1.01, 3.17, 1.28, 6.25, 11.92, 11.38 ],
+            [ np.nan, np.nan, np.nan, np.nan, np.nan, np.nan ],
+            [ np.nan, np.nan, np.nan, np.nan, np.nan, np.nan ],
+            [ 0.73, 1.05, 0.53, 1.5, 3.77, 3.51 ],
+        ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_ElasticStrainFromYoungPoisson( self: Self ) -> None:
         """Test calculation of Elastic Strain."""
-        deltaEffectiveStress: npt.NDArray[np.float64] = -1.3 * effectiveStress
-        deltaEffectiveStress[np.isnan(deltaEffectiveStress)] = 150.0
-        young: npt.NDArray[np.float64] = fcts.youngModulus(bulkModulus, shearModulus)
-        obtained: npt.NDArray[np.float64] = fcts.elasticStrainFromYoungPoisson(
-            deltaEffectiveStress, young, poissonRatio
-        )
-        print(np.round(obtained, 2).tolist())
-
-        expected: npt.NDArray[np.float64] = np.array(
-            [
-                [2.09, 20.36, 6.15, 46.84, 81.19, 77.13],
-                [1.04, 3.2, 1.31, 6.24, 11.89, 11.35],
-                [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-                [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan],
-                [0.72, 1.04, 0.52, 1.5, 3.78, 3.52],
-            ]
-        )
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_DeviatoricStressPathOed(self: Self) -> None:
+        deltaEffectiveStress: npt.NDArray[ np.float64 ] = -1.3 * effectiveStress
+        deltaEffectiveStress[ np.isnan( deltaEffectiveStress ) ] = 150.0
+        young: npt.NDArray[ np.float64 ] = fcts.youngModulus( bulkModulus, shearModulus )
+        obtained: npt.NDArray[ np.float64 ] = fcts.elasticStrainFromYoungPoisson( deltaEffectiveStress, young,
+                                                                                  poissonRatio )
+        print( np.round( obtained, 2 ).tolist() )
+
+        expected: npt.NDArray[ np.float64 ] = np.array( [
+            [ 2.09, 20.36, 6.15, 46.84, 81.19, 77.13 ],
+            [ 1.04, 3.2, 1.31, 6.24, 11.89, 11.35 ],
+            [ np.nan, np.nan, np.nan, np.nan, np.nan, np.nan ],
+            [ np.nan, np.nan, np.nan, np.nan, np.nan, np.nan ],
+            [ 0.72, 1.04, 0.52, 1.5, 3.78, 3.52 ],
+        ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_DeviatoricStressPathOed( self: Self ) -> None:
         """Test calculation of deviatoric Stress Path Oedometric conditions."""
-        obtained: npt.NDArray[np.float64] = fcts.deviatoricStressPathOed(poissonRatio)
-        expected: npt.NDArray[np.float64] = np.array([0.52, 0.41, np.nan, np.nan, 0.3])
+        obtained: npt.NDArray[ np.float64 ] = fcts.deviatoricStressPathOed( poissonRatio )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 0.52, 0.41, np.nan, np.nan, 0.3 ] )
 
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
 
-    def test_ReservoirStressPathReal(self: Self) -> None:
+    def test_ReservoirStressPathReal( self: Self ) -> None:
         """Test calculation of reservoir Stress Path Real conditions."""
-        deltaStress: npt.NDArray[np.float64] = -0.1 * effectiveStress
-        deltaPressure: npt.NDArray[np.float64] = 0.2 * pressure
-        obtained: npt.NDArray[np.float64] = fcts.reservoirStressPathReal(
-            deltaStress, deltaPressure
-        )
-        expected: npt.NDArray[np.float64] = np.array(
-            [
-                [8.0, 12.5, 9.0],
-                [0.32, 0.46, 0.34],
-                [np.nan, np.nan, np.nan],
-                [1.38, 1.81, 1.38],
-                [1.30, 1.55, 1.15],
-            ]
-        )
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_ReservoirStressPathOed(self: Self) -> None:
+        deltaStress: npt.NDArray[ np.float64 ] = -0.1 * effectiveStress
+        deltaPressure: npt.NDArray[ np.float64 ] = 0.2 * pressure
+        obtained: npt.NDArray[ np.float64 ] = fcts.reservoirStressPathReal( deltaStress, deltaPressure )
+        expected: npt.NDArray[ np.float64 ] = np.array( [
+            [ 8.0, 12.5, 9.0 ],
+            [ 0.32, 0.46, 0.34 ],
+            [ np.nan, np.nan, np.nan ],
+            [ 1.38, 1.81, 1.38 ],
+            [ 1.30, 1.55, 1.15 ],
+        ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_ReservoirStressPathOed( self: Self ) -> None:
         """Test calculation of reservoir Stress Path Oedometric conditions."""
-        obtained: npt.NDArray[np.float64] = fcts.reservoirStressPathOed(
-            biotCoefficient, poissonRatio
-        )
-        expected: npt.NDArray[np.float64] = np.array([0.05, 0.18, np.nan, np.nan, 0.42])
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_criticalTotalStressRatio(self: Self) -> None:
+        obtained: npt.NDArray[ np.float64 ] = fcts.reservoirStressPathOed( biotCoefficient, poissonRatio )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 0.05, 0.18, np.nan, np.nan, 0.42 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_criticalTotalStressRatio( self: Self ) -> None:
         """Test calculation of fracture Index."""
-        obtained: npt.NDArray[np.float64] = fcts.criticalTotalStressRatio(
-            pressure, verticalStress
-        )
-        expected: npt.NDArray[np.float64] = np.array([0.06, 1.47, np.nan, 0.36, 0.43])
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_totalStressRatioThreshold(self: Self) -> None:
+        obtained: npt.NDArray[ np.float64 ] = fcts.criticalTotalStressRatio( pressure, verticalStress )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 0.06, 1.47, np.nan, 0.36, 0.43 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_totalStressRatioThreshold( self: Self ) -> None:
         """Test calculation of fracture Threshold."""
-        obtained: npt.NDArray[np.float64] = fcts.totalStressRatioThreshold(
-            pressure, horizontalStress
-        )
-        expected: npt.NDArray[np.float64] = np.array([0.04, 1.08, np.nan, 0.28, 0.32])
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_CriticalPorePressure(self: Self) -> None:
+        obtained: npt.NDArray[ np.float64 ] = fcts.totalStressRatioThreshold( pressure, horizontalStress )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 0.04, 1.08, np.nan, 0.28, 0.32 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_CriticalPorePressure( self: Self ) -> None:
         """Test calculation of critical Pore Pressure."""
         rockCohesion: float = 20  # GPa
         frictionAngle: float = 30 * np.pi / 180.0  # friction angle in radian
-        effectiveStress2: npt.NDArray[np.float64] = -1.0 * effectiveStress
-        effectiveStress2[np.isnan(effectiveStress)] = 180.0
-        obtained: npt.NDArray[np.float64] = fcts.criticalPorePressure(
-            effectiveStress2, rockCohesion, frictionAngle
-        )
-        expected: npt.NDArray[np.float64] = np.array(
-            [-302.43, -333.72, -348.1, -383.01, -438.03]
-        )
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_CriticalPorePressureThreshold(self: Self) -> None:
+        effectiveStress2: npt.NDArray[ np.float64 ] = -1.0 * effectiveStress
+        effectiveStress2[ np.isnan( effectiveStress ) ] = 180.0
+        obtained: npt.NDArray[ np.float64 ] = fcts.criticalPorePressure( effectiveStress2, rockCohesion, frictionAngle )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ -302.43, -333.72, -348.1, -383.01, -438.03 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_CriticalPorePressureThreshold( self: Self ) -> None:
         """Test calculation of critical Pore Pressure Index."""
-        criticalPorePressure: npt.NDArray[np.float64] = np.array(
-            [149.64, 174.64, 194.64, 219.64, 224.64]
-        )
-
-        obtained: npt.NDArray[np.float64] = fcts.criticalPorePressureThreshold(
-            pressure, criticalPorePressure
-        )
-        expected: npt.NDArray[np.float64] = np.array([0.07, 1.6, np.nan, 0.36, 0.45])
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    def test_compressibility(self: Self) -> None:
+        criticalPorePressure: npt.NDArray[ np.float64 ] = np.array( [ 149.64, 174.64, 194.64, 219.64, 224.64 ] )
+
+        obtained: npt.NDArray[ np.float64 ] = fcts.criticalPorePressureThreshold( pressure, criticalPorePressure )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 0.07, 1.6, np.nan, 0.36, 0.45 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    def test_compressibility( self: Self ) -> None:
         """Test calculation of compressibility from elastic moduli."""
-        porosity: npt.NDArray[np.float64] = np.array([0.05, 0.1, 0.15, 0.30, 0.20])
+        porosity: npt.NDArray[ np.float64 ] = np.array( [ 0.05, 0.1, 0.15, 0.30, 0.20 ] )
 
         # multiply by 1000 to prevent from to small values
-        obtained: npt.NDArray[np.float64] = 1000.0 * fcts.compressibility(
-            poissonRatio, bulkModulus, biotCoefficient, porosity
-        )
-        expected: npt.NDArray[np.float64] = np.array(
-            [110.438, 49.016, np.nan, np.nan, 18.991]
-        )
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 3), np.round(expected, 3), equal_nan=True)
-        )
-
-    def test_shearCapacityUtilization(self: Self) -> None:
+        obtained: npt.NDArray[ np.float64 ] = 1000.0 * fcts.compressibility( poissonRatio, bulkModulus, biotCoefficient,
+                                                                             porosity )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 110.438, 49.016, np.nan, np.nan, 18.991 ] )
+        self.assertTrue( np.array_equal( np.round( obtained, 3 ), np.round( expected, 3 ), equal_nan=True ) )
+
+    def test_shearCapacityUtilization( self: Self ) -> None:
         """Test calculation of shear capacity utilization."""
         # inputs
-        traction: npt.NDArray[np.float64] = np.copy(effectiveStress[:, :3])
+        traction: npt.NDArray[ np.float64 ] = np.copy( effectiveStress[ :, :3 ] )
         # replace nan values for calculation
-        traction[np.isnan(traction)] = 150.0
+        traction[ np.isnan( traction ) ] = 150.0
         rockCohesion = 250.0
         frictionAngle = 10.0 * np.pi / 180.0
 
         # calculation
-        obtained: npt.NDArray[np.float64] = fcts.shearCapacityUtilization(
-            traction, rockCohesion, frictionAngle
-        )
-        expected: npt.NDArray[np.float64] = np.array(
-            [0.899, 0.923, 0.982, 1.004, 1.048]
-        )
-
-        self.assertSequenceEqual(
-            np.round(obtained, 3).flatten().tolist(), expected.tolist()
-        )
-
-    def test_computeStressPrincipalComponents(self: Self) -> None:
+        obtained: npt.NDArray[ np.float64 ] = fcts.shearCapacityUtilization( traction, rockCohesion, frictionAngle )
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 0.899, 0.923, 0.982, 1.004, 1.048 ] )
+
+        self.assertSequenceEqual( np.round( obtained, 3 ).flatten().tolist(), expected.tolist() )
+
+    def test_computeStressPrincipalComponents( self: Self ) -> None:
         """Test calculation of stress principal components from stress tensor."""
-        obtained: list[float] = [
-            round(val, 3) for val in fcts.computeStressPrincipalComponents(stressTensor)
-        ]
-        expected: tuple[float, float, float] = (1.748, 2.471, 5.281)
-        self.assertSequenceEqual(tuple(obtained), expected)
+        obtained: list[ float ] = [ round( val, 3 ) for val in fcts.computeStressPrincipalComponents( stressTensor ) ]
+        expected: tuple[ float, float, float ] = ( 1.748, 2.471, 5.281 )
+        self.assertSequenceEqual( tuple( obtained ), expected )
 
-    def test_computeStressPrincipalComponentsFromStressVector(self: Self) -> None:
+    def test_computeStressPrincipalComponentsFromStressVector( self: Self ) -> None:
         """Test calculation of stress principal components from stress vector."""
-        obtained: list[float] = [
-            round(val, 3)
-            for val in fcts.computeStressPrincipalComponentsFromStressVector(
-                stressVector
-            )
+        obtained: list[ float ] = [
+            round( val, 3 ) for val in fcts.computeStressPrincipalComponentsFromStressVector( stressVector )
         ]
-        expected: tuple[float, float, float] = (1.748, 2.471, 5.281)
-        self.assertSequenceEqual(tuple(obtained), expected)
+        expected: tuple[ float, float, float ] = ( 1.748, 2.471, 5.281 )
+        self.assertSequenceEqual( tuple( obtained ), expected )
 
-    def test_computeNormalShearStress(self: Self) -> None:
+    def test_computeNormalShearStress( self: Self ) -> None:
         """Test calculation of normal and shear stress."""
-        directionVector = np.array([1.0, 1.0, 0.0])
-        obtained: list[float] = [
-            round(val, 3)
-            for val in fcts.computeNormalShearStress(stressTensor, directionVector)
+        directionVector = np.array( [ 1.0, 1.0, 0.0 ] )
+        obtained: list[ float ] = [
+            round( val, 3 ) for val in fcts.computeNormalShearStress( stressTensor, directionVector )
         ]
-        expected: tuple[float, float] = (2.250, 0.606)
-        self.assertSequenceEqual(tuple(obtained), expected)
+        expected: tuple[ float, float ] = ( 2.250, 0.606 )
+        self.assertSequenceEqual( tuple( obtained ), expected )
 
 
 if __name__ == "__main__":
diff --git a/geos-geomechanics/tests/testsMohrCircle.py b/geos-geomechanics/tests/testsMohrCircle.py
index 133b30b3..cd3b2e55 100644
--- a/geos-geomechanics/tests/testsMohrCircle.py
+++ b/geos-geomechanics/tests/testsMohrCircle.py
@@ -10,285 +10,257 @@
 import numpy.typing as npt
 from typing_extensions import Self
 
-dir_path = os.path.dirname(os.path.realpath(__file__))
-parent_dir_path = os.path.join(os.path.dirname(dir_path), "src")
+dir_path = os.path.dirname( os.path.realpath( __file__ ) )
+parent_dir_path = os.path.join( os.path.dirname( dir_path ), "src" )
 if parent_dir_path not in sys.path:
-    sys.path.append(parent_dir_path)
+    sys.path.append( parent_dir_path )
 
 from geos.geomechanics.model.MohrCircle import MohrCircle
 from geos.geomechanics.model.MohrCoulomb import MohrCoulomb
 
 circleId = "12453"
-stressVector: npt.NDArray[np.float64] = np.array([1.8, 2.5, 5.2, 0.3, 0.4, 0.1])
+stressVector: npt.NDArray[ np.float64 ] = np.array( [ 1.8, 2.5, 5.2, 0.3, 0.4, 0.1 ] )
 
 rockCohesion: float = 8.0e8
 frictionAngle: float = 10 * np.pi / 180.0  # in radian
 
-principalComponentsExpected: tuple[float, float, float] = (1.748, 2.471, 5.281)
+principalComponentsExpected: tuple[ float, float, float ] = ( 1.748, 2.471, 5.281 )
 
 
-class TestsMohrCircle(unittest.TestCase):
+class TestsMohrCircle( unittest.TestCase ):
 
-    def test_MohrCircleInit(self: Self) -> None:
+    def test_MohrCircleInit( self: Self ) -> None:
         """Test instanciation of MohrCircle objects."""
-        mohrCircle: MohrCircle = MohrCircle(circleId)
-        mohrCircle.setPrincipalComponents(*principalComponentsExpected)
+        mohrCircle: MohrCircle = MohrCircle( circleId )
+        mohrCircle.setPrincipalComponents( *principalComponentsExpected )
 
-        self.assertSequenceEqual(circleId, mohrCircle.getCircleId())
+        self.assertSequenceEqual( circleId, mohrCircle.getCircleId() )
 
-        obtained: list[float] = [
-            round(val, 3) for val in mohrCircle.getPrincipalComponents()
-        ]
-        self.assertSequenceEqual(obtained, principalComponentsExpected)
+        obtained: list[ float ] = [ round( val, 3 ) for val in mohrCircle.getPrincipalComponents() ]
+        self.assertSequenceEqual( obtained, principalComponentsExpected )
 
-    def test_MohrCircleRadius(self: Self) -> None:
+    def test_MohrCircleRadius( self: Self ) -> None:
         """Test the calculation of MohrCircle radius."""
-        mohrCircle: MohrCircle = MohrCircle(circleId)
-        mohrCircle.computePrincipalComponents(stressVector)
+        mohrCircle: MohrCircle = MohrCircle( circleId )
+        mohrCircle.computePrincipalComponents( stressVector )
 
         obtained: float = mohrCircle.getCircleRadius()
-        expected: float = (
-            principalComponentsExpected[2] - principalComponentsExpected[0]
-        ) / 2.0
-        self.assertAlmostEqual(obtained, expected, 3)
+        expected: float = ( principalComponentsExpected[ 2 ] - principalComponentsExpected[ 0 ] ) / 2.0
+        self.assertAlmostEqual( obtained, expected, 3 )
 
-    def test_MohrCircleCenter(self: Self) -> None:
+    def test_MohrCircleCenter( self: Self ) -> None:
         """Test the calculation of MohrCircle center."""
-        mohrCircle: MohrCircle = MohrCircle(circleId)
-        mohrCircle.computePrincipalComponents(stressVector)
+        mohrCircle: MohrCircle = MohrCircle( circleId )
+        mohrCircle.computePrincipalComponents( stressVector )
 
         obtained: float = mohrCircle.getCircleCenter()
-        expected: float = (
-            principalComponentsExpected[2] + principalComponentsExpected[0]
-        ) / 2.0
-        self.assertAlmostEqual(obtained, expected, 3)
+        expected: float = ( principalComponentsExpected[ 2 ] + principalComponentsExpected[ 0 ] ) / 2.0
+        self.assertAlmostEqual( obtained, expected, 3 )
 
 
-class TestsMohrCoulomb(unittest.TestCase):
+class TestsMohrCoulomb( unittest.TestCase ):
 
-    def test_MohrCoulombInit(self: Self) -> None:
+    def test_MohrCoulombInit( self: Self ) -> None:
         """Test instanciation of MohrCoulomb objects."""
         # expected values
-        expectedSlope: float = np.tan(frictionAngle)
+        expectedSlope: float = np.tan( frictionAngle )
         expectedSigmaMin: float = -rockCohesion / expectedSlope
         # instantiate object
-        mohrCoulomb: MohrCoulomb = MohrCoulomb(rockCohesion, frictionAngle)
+        mohrCoulomb: MohrCoulomb = MohrCoulomb( rockCohesion, frictionAngle )
         # test results
-        self.assertEqual(mohrCoulomb.m_rockCohesion, rockCohesion)
-        self.assertEqual(mohrCoulomb.m_slope, expectedSlope)
-        self.assertEqual(mohrCoulomb.m_sigmaMin, expectedSigmaMin)
+        self.assertEqual( mohrCoulomb.m_rockCohesion, rockCohesion )
+        self.assertEqual( mohrCoulomb.m_slope, expectedSlope )
+        self.assertEqual( mohrCoulomb.m_sigmaMin, expectedSigmaMin )
 
-    def test_computeShearStress(self: Self) -> None:
+    def test_computeShearStress( self: Self ) -> None:
         """Test calculation of shear stress from normal stress."""
         # inputs
-        stressNormal: npt.NDArray[np.float64] = np.linspace(5.0e8, 1.0e9, 100)
+        stressNormal: npt.NDArray[ np.float64 ] = np.linspace( 5.0e8, 1.0e9, 100 )
         # expected values
-        expectedValues: npt.NDArray[np.float64] = np.array(
-            [
-                888163490.0,
-                889054031.0,
-                889944571.0,
-                890835111.0,
-                891725652.0,
-                892616192.0,
-                893506732.0,
-                894397273.0,
-                895287813.0,
-                896178353.0,
-                897068893.0,
-                897959434.0,
-                898849974.0,
-                899740514.0,
-                900631055.0,
-                901521595.0,
-                902412135.0,
-                903302676.0,
-                904193216.0,
-                905083756.0,
-                905974296.0,
-                906864837.0,
-                907755377.0,
-                908645917.0,
-                909536458.0,
-                910426998.0,
-                911317538.0,
-                912208079.0,
-                913098619.0,
-                913989159.0,
-                914879700.0,
-                915770240.0,
-                916660780.0,
-                917551320.0,
-                918441861.0,
-                919332401.0,
-                920222941.0,
-                921113482.0,
-                922004022.0,
-                922894562.0,
-                923785103.0,
-                924675643.0,
-                925566183.0,
-                926456724.0,
-                927347264.0,
-                928237804.0,
-                929128344.0,
-                930018885.0,
-                930909425.0,
-                931799965.0,
-                932690506.0,
-                933581046.0,
-                934471586.0,
-                935362127.0,
-                936252667.0,
-                937143207.0,
-                938033748.0,
-                938924288.0,
-                939814828.0,
-                940705368.0,
-                941595909.0,
-                942486449.0,
-                943376989.0,
-                944267530.0,
-                945158070.0,
-                946048610.0,
-                946939151.0,
-                947829691.0,
-                948720231.0,
-                949610772.0,
-                950501312.0,
-                951391852.0,
-                952282392.0,
-                953172933.0,
-                954063473.0,
-                954954013.0,
-                955844554.0,
-                956735094.0,
-                957625634.0,
-                958516175.0,
-                959406715.0,
-                960297255.0,
-                961187795.0,
-                962078336.0,
-                962968876.0,
-                963859416.0,
-                964749957.0,
-                965640497.0,
-                966531037.0,
-                967421578.0,
-                968312118.0,
-                969202658.0,
-                970093199.0,
-                970983739.0,
-                971874279.0,
-                972764819.0,
-                973655360.0,
-                974545900.0,
-                975436440.0,
-                976326981.0,
-            ]
-        )
+        expectedValues: npt.NDArray[ np.float64 ] = np.array( [
+            888163490.0,
+            889054031.0,
+            889944571.0,
+            890835111.0,
+            891725652.0,
+            892616192.0,
+            893506732.0,
+            894397273.0,
+            895287813.0,
+            896178353.0,
+            897068893.0,
+            897959434.0,
+            898849974.0,
+            899740514.0,
+            900631055.0,
+            901521595.0,
+            902412135.0,
+            903302676.0,
+            904193216.0,
+            905083756.0,
+            905974296.0,
+            906864837.0,
+            907755377.0,
+            908645917.0,
+            909536458.0,
+            910426998.0,
+            911317538.0,
+            912208079.0,
+            913098619.0,
+            913989159.0,
+            914879700.0,
+            915770240.0,
+            916660780.0,
+            917551320.0,
+            918441861.0,
+            919332401.0,
+            920222941.0,
+            921113482.0,
+            922004022.0,
+            922894562.0,
+            923785103.0,
+            924675643.0,
+            925566183.0,
+            926456724.0,
+            927347264.0,
+            928237804.0,
+            929128344.0,
+            930018885.0,
+            930909425.0,
+            931799965.0,
+            932690506.0,
+            933581046.0,
+            934471586.0,
+            935362127.0,
+            936252667.0,
+            937143207.0,
+            938033748.0,
+            938924288.0,
+            939814828.0,
+            940705368.0,
+            941595909.0,
+            942486449.0,
+            943376989.0,
+            944267530.0,
+            945158070.0,
+            946048610.0,
+            946939151.0,
+            947829691.0,
+            948720231.0,
+            949610772.0,
+            950501312.0,
+            951391852.0,
+            952282392.0,
+            953172933.0,
+            954063473.0,
+            954954013.0,
+            955844554.0,
+            956735094.0,
+            957625634.0,
+            958516175.0,
+            959406715.0,
+            960297255.0,
+            961187795.0,
+            962078336.0,
+            962968876.0,
+            963859416.0,
+            964749957.0,
+            965640497.0,
+            966531037.0,
+            967421578.0,
+            968312118.0,
+            969202658.0,
+            970093199.0,
+            970983739.0,
+            971874279.0,
+            972764819.0,
+            973655360.0,
+            974545900.0,
+            975436440.0,
+            976326981.0,
+        ] )
         # instantiate object
-        mohrCoulomb: MohrCoulomb = MohrCoulomb(rockCohesion, frictionAngle)
-        obtained: npt.NDArray[np.float64] = np.array(
-            mohrCoulomb.computeShearStress(stressNormal)
-        )
+        mohrCoulomb: MohrCoulomb = MohrCoulomb( rockCohesion, frictionAngle )
+        obtained: npt.NDArray[ np.float64 ] = np.array( mohrCoulomb.computeShearStress( stressNormal ) )
 
         # test results
-        self.assertSequenceEqual(expectedValues.tolist(), np.round(obtained).tolist())
+        self.assertSequenceEqual( expectedValues.tolist(), np.round( obtained ).tolist() )
 
-    def test_computeFailureEnvelop1(self: Self) -> None:
+    def test_computeFailureEnvelop1( self: Self ) -> None:
         """Test calculation of failure envelop from minimum normal stress."""
         # inputs
         stressNormalMax: float = 1.0e9
         # expected values
-        normalStressExpected: npt.NDArray[np.float64] = np.array(
-            [
-                -4537025456.0,
-                -3921800405.0,
-                -3306575354.0,
-                -2691350304.0,
-                -2076125253.0,
-                -1460900203.0,
-                -845675152.0,
-                -230450101.0,
-                384774949.0,
-                1000000000.0,
-            ]
-        )
-        shearStressExpected: npt.NDArray[np.float64] = np.array(
-            [
-                0.0,
-                108480776.0,
-                216961551.0,
-                325442327.0,
-                433923103.0,
-                542403878.0,
-                650884654.0,
-                759365429.0,
-                867846205.0,
-                976326981.0,
-            ]
-        )
+        normalStressExpected: npt.NDArray[ np.float64 ] = np.array( [
+            -4537025456.0,
+            -3921800405.0,
+            -3306575354.0,
+            -2691350304.0,
+            -2076125253.0,
+            -1460900203.0,
+            -845675152.0,
+            -230450101.0,
+            384774949.0,
+            1000000000.0,
+        ] )
+        shearStressExpected: npt.NDArray[ np.float64 ] = np.array( [
+            0.0,
+            108480776.0,
+            216961551.0,
+            325442327.0,
+            433923103.0,
+            542403878.0,
+            650884654.0,
+            759365429.0,
+            867846205.0,
+            976326981.0,
+        ] )
 
         # instantiate object
-        mohrCoulomb: MohrCoulomb = MohrCoulomb(rockCohesion, frictionAngle)
-        normalStressObtained, shearStressObtained = mohrCoulomb.computeFailureEnvelop(
-            stressNormalMax, n=10
-        )
+        mohrCoulomb: MohrCoulomb = MohrCoulomb( rockCohesion, frictionAngle )
+        normalStressObtained, shearStressObtained = mohrCoulomb.computeFailureEnvelop( stressNormalMax, n=10 )
 
         # test results
-        self.assertSequenceEqual(
-            normalStressExpected.tolist(), np.round(normalStressObtained).tolist()
-        )
-        self.assertSequenceEqual(
-            shearStressExpected.tolist(), np.round(shearStressObtained).tolist()
-        )
+        self.assertSequenceEqual( normalStressExpected.tolist(), np.round( normalStressObtained ).tolist() )
+        self.assertSequenceEqual( shearStressExpected.tolist(), np.round( shearStressObtained ).tolist() )
 
-    def test_computeFailureEnvelop2(self: Self) -> None:
+    def test_computeFailureEnvelop2( self: Self ) -> None:
         """Test calculation of failure envelop in user-defined range."""
         # inputs
         stressNormalMin: float = 8.0e8
         stressNormalMax: float = 1.0e9
         # expected values
-        normalStressExpected: npt.NDArray[np.float64] = np.array(
-            [
-                800000000.0,
-                822222222.0,
-                844444444.0,
-                866666667.0,
-                888888889.0,
-                911111111.0,
-                933333333.0,
-                955555556.0,
-                977777778.0,
-                1000000000.0,
-            ]
-        )
-        shearStressExpected: npt.NDArray[np.float64] = np.array(
-            [
-                941061585.0,
-                944979962.0,
-                948898339.0,
-                952816717.0,
-                956735094.0,
-                960653471.0,
-                964571849.0,
-                968490226.0,
-                972408603.0,
-                976326981.0,
-            ]
-        )
+        normalStressExpected: npt.NDArray[ np.float64 ] = np.array( [
+            800000000.0,
+            822222222.0,
+            844444444.0,
+            866666667.0,
+            888888889.0,
+            911111111.0,
+            933333333.0,
+            955555556.0,
+            977777778.0,
+            1000000000.0,
+        ] )
+        shearStressExpected: npt.NDArray[ np.float64 ] = np.array( [
+            941061585.0,
+            944979962.0,
+            948898339.0,
+            952816717.0,
+            956735094.0,
+            960653471.0,
+            964571849.0,
+            968490226.0,
+            972408603.0,
+            976326981.0,
+        ] )
 
         # instantiate object
-        mohrCoulomb: MohrCoulomb = MohrCoulomb(rockCohesion, frictionAngle)
-        normalStressObtained, shearStressObtained = mohrCoulomb.computeFailureEnvelop(
-            stressNormalMax, stressNormalMin, n=10
-        )
+        mohrCoulomb: MohrCoulomb = MohrCoulomb( rockCohesion, frictionAngle )
+        normalStressObtained, shearStressObtained = mohrCoulomb.computeFailureEnvelop( stressNormalMax,
+                                                                                       stressNormalMin,
+                                                                                       n=10 )
 
         # test results
-        self.assertSequenceEqual(
-            normalStressExpected.tolist(), np.round(normalStressObtained).tolist()
-        )
-        self.assertSequenceEqual(
-            shearStressExpected.tolist(), np.round(shearStressObtained).tolist()
-        )
+        self.assertSequenceEqual( normalStressExpected.tolist(), np.round( normalStressObtained ).tolist() )
+        self.assertSequenceEqual( shearStressExpected.tolist(), np.round( shearStressObtained ).tolist() )

From 86c892096646d8d7dc06a522e952a1b6dfa68d68 Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Wed, 2 Apr 2025 13:35:20 +0200
Subject: [PATCH 13/14] Fix doc

---
 docs/geos_geomechanics_docs/processing.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/geos_geomechanics_docs/processing.rst b/docs/geos_geomechanics_docs/processing.rst
index 35e416a7..78c9e462 100644
--- a/docs/geos_geomechanics_docs/processing.rst
+++ b/docs/geos_geomechanics_docs/processing.rst
@@ -3,7 +3,7 @@ Processing
 The processing part of `geos-geomechanics` package contains functions tools to compute geomechanical properties.
 
 geos.geomechanics.processing.geomechanicsCalculatorFunctions module
----------------------------------------------------------------
+---------------------------------------------------------------------
 
 .. automodule:: geos.geomechanics.processing.geomechanicsCalculatorFunctions
    :members:

From 7fdd8c38893f27cdf98cadc835d8a2b141b65c4b Mon Sep 17 00:00:00 2001
From: Paloma Martinez <104762252+paloma-martinez@users.noreply.github.com>
Date: Wed, 2 Apr 2025 15:12:11 +0200
Subject: [PATCH 14/14] Consistency with other packages pyproject.toml

---
 geos-geomechanics/pyproject.toml | 19 ++++---------------
 1 file changed, 4 insertions(+), 15 deletions(-)

diff --git a/geos-geomechanics/pyproject.toml b/geos-geomechanics/pyproject.toml
index 0ba2c2d6..6ec7b661 100644
--- a/geos-geomechanics/pyproject.toml
+++ b/geos-geomechanics/pyproject.toml
@@ -22,23 +22,18 @@ classifiers = [
     "Programming Language :: Python"
 ]
 dynamic = ["dependencies"]
-requires-python = ">= 3.9"
+requires-python = ">= 3.10"
 
 [project.optional-dependencies]
 build = [
   "build ~= 1.2"
 ]
 dev = [
-    "pylint",
+    "yapf",
     "mypy",
-    "black",
-    "sphinx", 
-    "sphinx-rtd-theme", 
-    "sphinx-autodoc-typehints"
 ]
 test = [
    "pytest",
-   "pytest-cov"
 ]
 
 [project.scripts]
@@ -47,7 +42,7 @@ test = [
 [tool.pytest.ini_options]
 addopts = "--import-mode=importlib"
 console_output_style = "count"
-pythonpath = [".", "src"]
+pythonpath = ["src"]
 python_classes = "Test"
 python_files = "test*.py"
 python_functions = "test*"
@@ -58,10 +53,4 @@ filterwarnings = []
 
 [tool.coverage.run]
 branch = true
-source = ["geos/geomechanics"]
-
-
-[tool.mypy]
-python_version = "3.9"
-warn_return_any = true
-warn_unused_configs = true
+source = ["src/geos"]
\ No newline at end of file