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/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..d1b98d41
--- /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..78c9e462
--- /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
diff --git a/geos-geomechanics/pyproject.toml b/geos-geomechanics/pyproject.toml
new file mode 100644
index 00000000..6ec7b661
--- /dev/null
+++ b/geos-geomechanics/pyproject.toml
@@ -0,0 +1,56 @@
+[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"
+]
+dynamic = ["dependencies"]
+requires-python = ">= 3.10"
+
+[project.optional-dependencies]
+build = [
+  "build ~= 1.2"
+]
+dev = [
+    "yapf",
+    "mypy",
+]
+test = [
+   "pytest",
+]
+
+[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 = ["src/geos"]
\ No newline at end of file
diff --git a/geos-geomechanics/setup.py b/geos-geomechanics/setup.py
new file mode 100644
index 00000000..3c542ff6
--- /dev/null
+++ b/geos-geomechanics/setup.py
@@ -0,0 +1,13 @@
+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}",
+] )
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-posp/src/geos_posp/processing/MohrCircle.py b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
similarity index 70%
rename from geos-posp/src/geos_posp/processing/MohrCircle.py
rename to geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
index 5ec2cef9..a917c3de 100644
--- a/geos-posp/src/geos_posp/processing/MohrCircle.py
+++ b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
@@ -1,119 +1,115 @@
-# 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 ) )
diff --git a/geos-posp/src/geos_posp/processing/MohrCoulomb.py b/geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py
similarity index 75%
rename from geos-posp/src/geos_posp/processing/MohrCoulomb.py
rename to geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py
index 4072191b..93409149 100644
--- a/geos-posp/src/geos_posp/processing/MohrCoulomb.py
+++ b/geos-geomechanics/src/geos/geomechanics/model/MohrCoulomb.py
@@ -1,99 +1,96 @@
-# 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 ) ) )
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-posp/src/geos_posp/processing/geomechanicsCalculatorFunctions.py b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
similarity index 56%
rename from geos-posp/src/geos_posp/processing/geomechanicsCalculatorFunctions.py
rename to geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
index c1f700d0..48317c34 100644
--- a/geos-posp/src/geos_posp/processing/geomechanicsCalculatorFunctions.py
+++ b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
@@ -1,1011 +1,925 @@
-# 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.utils.algebraFunctions import getAttributeMatrixFromVector
-from geos.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.model.MohrCoulomb import MohrCoulomb
+from geos.utils.algebraFunctions import getAttributeMatrixFromVector
+from geos.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 )
diff --git a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
new file mode 100644
index 00000000..82b90bff
--- /dev/null
+++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
@@ -0,0 +1,279 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Alexandre Benedicto, Martin Lemay
+# SPDX-License-Identifier: Apache 2.0
+# ruff: noqa: E402 # disable Module level import not at top of file
+import os
+import sys
+import unittest
+
+import numpy as np
+import numpy.typing as npt
+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" )
+if parent_dir_path not in sys.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:
+        """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 ) )
+
+    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:
+        """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 ) )
+
+    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 )
+
+        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:
+        """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.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:
+        """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:
+        """Test calculation of oedometric Modulus."""
+        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 ) )
+
+    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 ) )
+
+    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:
+        """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:
+        """Test calculation of lithostatic Stress."""
+        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:
+        """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:
+        """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:
+        """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 ] )
+
+        self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) )
+
+    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:
+        """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:
+        """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:
+        """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:
+        """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:
+        """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:
+        """Test calculation of compressibility from elastic moduli."""
+        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:
+        """Test calculation of shear capacity utilization."""
+        # inputs
+        traction: npt.NDArray[ np.float64 ] = np.copy( effectiveStress[ :, :3 ] )
+        # replace nan values for calculation
+        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:
+        """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 )
+
+    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 )
+        ]
+        expected: tuple[ float, float, float ] = ( 1.748, 2.471, 5.281 )
+        self.assertSequenceEqual( tuple( obtained ), expected )
+
+    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 )
+        ]
+        expected: tuple[ float, float ] = ( 2.250, 0.606 )
+        self.assertSequenceEqual( tuple( obtained ), expected )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/geos-geomechanics/tests/testsMohrCircle.py b/geos-geomechanics/tests/testsMohrCircle.py
new file mode 100644
index 00000000..cd3b2e55
--- /dev/null
+++ b/geos-geomechanics/tests/testsMohrCircle.py
@@ -0,0 +1,266 @@
+# SPDX-License-Identifier: Apache-2.0
+# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
+# SPDX-FileContributor: Martin Lemay
+# ruff: noqa: E402 # disable Module level import not at top of file
+import os
+import sys
+import unittest
+
+import numpy as np
+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" )
+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
+
+circleId = "12453"
+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 )
+
+
+class TestsMohrCircle( unittest.TestCase ):
+
+    def test_MohrCircleInit( self: Self ) -> None:
+        """Test instanciation of MohrCircle objects."""
+        mohrCircle: MohrCircle = MohrCircle( circleId )
+        mohrCircle.setPrincipalComponents( *principalComponentsExpected )
+
+        self.assertSequenceEqual( circleId, mohrCircle.getCircleId() )
+
+        obtained: list[ float ] = [ round( val, 3 ) for val in mohrCircle.getPrincipalComponents() ]
+        self.assertSequenceEqual( obtained, principalComponentsExpected )
+
+    def test_MohrCircleRadius( self: Self ) -> None:
+        """Test the calculation of MohrCircle radius."""
+        mohrCircle: MohrCircle = MohrCircle( circleId )
+        mohrCircle.computePrincipalComponents( stressVector )
+
+        obtained: float = mohrCircle.getCircleRadius()
+        expected: float = ( principalComponentsExpected[ 2 ] - principalComponentsExpected[ 0 ] ) / 2.0
+        self.assertAlmostEqual( obtained, expected, 3 )
+
+    def test_MohrCircleCenter( self: Self ) -> None:
+        """Test the calculation of MohrCircle center."""
+        mohrCircle: MohrCircle = MohrCircle( circleId )
+        mohrCircle.computePrincipalComponents( stressVector )
+
+        obtained: float = mohrCircle.getCircleCenter()
+        expected: float = ( principalComponentsExpected[ 2 ] + principalComponentsExpected[ 0 ] ) / 2.0
+        self.assertAlmostEqual( obtained, expected, 3 )
+
+
+class TestsMohrCoulomb( unittest.TestCase ):
+
+    def test_MohrCoulombInit( self: Self ) -> None:
+        """Test instanciation of MohrCoulomb objects."""
+        # expected values
+        expectedSlope: float = np.tan( frictionAngle )
+        expectedSigmaMin: float = -rockCohesion / expectedSlope
+        # instantiate object
+        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 )
+
+    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 )
+        # 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,
+        ] )
+        # instantiate object
+        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() )
+
+    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,
+        ] )
+
+        # instantiate object
+        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() )
+
+    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,
+        ] )
+
+        # instantiate object
+        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() )
diff --git a/geos-posp/setup.py b/geos-posp/setup.py
index d933e5f6..21bd8e95 100644
--- a/geos-posp/setup.py
+++ b/geos-posp/setup.py
@@ -2,8 +2,8 @@
 from setuptools import setup
 
 # This is where you add any fancy path resolution to the local lib:
-package_name = "geos-utils"
-geos_utils_path: str = (Path(__file__).parent.parent / package_name).as_uri()
+geos_geomecha_path: str = (Path(__file__).parent.parent / "geos-geomechanics").as_uri()
+geos_utils_path: str = (Path(__file__).parent.parent / "geos-utils").as_uri()
 
 setup(
     install_requires=[
@@ -11,6 +11,7 @@
         "numpy >= 1.26",
         "pandas >= 2.2",
         "typing_extensions >= 4.12",
-        f"{package_name} @ {geos_utils_path}",
+        f"geos-geomechanics @ {geos_geomecha_path}",
+        f"geos-utils @ {geos_utils_path}",
     ]
 )
\ No newline at end of file
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/PVplugins/PVMohrCirclePlot.py b/geos-posp/src/PVplugins/PVMohrCirclePlot.py
index c464e951..eb4c73b4 100644
--- a/geos-posp/src/PVplugins/PVMohrCirclePlot.py
+++ b/geos-posp/src/PVplugins/PVMohrCirclePlot.py
@@ -36,7 +36,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.utils.enumUnits import Pressure, enumerationDomainUnit
 from geos.utils.GeosOutputsConstants import (
diff --git a/geos-posp/src/PVplugins/__init__.py b/geos-posp/src/PVplugins/__init__.py
index 1fface1c..7cd2f2f3 100644
--- a/geos-posp/src/PVplugins/__init__.py
+++ b/geos-posp/src/PVplugins/__init__.py
@@ -11,4 +11,4 @@
 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)
+        sys.path.insert( 0, m_path)
\ No newline at end of file
diff --git a/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py b/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py
index f376a661..b8afd5cb 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,
@@ -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
diff --git a/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py b/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
index 2c541469..a1c6f094 100644
--- a/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
+++ b/geos-posp/src/geos_posp/filters/SurfaceGeomechanics.py
@@ -17,7 +17,8 @@
     vtkPolyData,
 )
 
-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/visu/mohrCircles/functionsMohrCircle.py b/geos-posp/src/geos_posp/visu/mohrCircles/functionsMohrCircle.py
index 832d8ed5..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_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 29257aa8..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_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.utils.enumUnits import Pressure, Unit, convert
 from geos.utils.GeosOutputsConstants import FAILURE_ENVELOPE
 from geos_posp.visu.PVUtils.matplotlibOptions import (
diff --git a/geos-posp/tests/testsFunctionsGeomechanicsCalculator.py b/geos-posp/tests/testsFunctionsGeomechanicsCalculator.py
deleted file mode 100644
index fcad9e0d..00000000
--- a/geos-posp/tests/testsFunctionsGeomechanicsCalculator.py
+++ /dev/null
@@ -1,382 +0,0 @@
-# SPDX-License-Identifier: Apache-2.0
-# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Alexandre Benedicto, Martin Lemay
-# SPDX-License-Identifier: Apache 2.0
-# ruff: noqa: E402 # disable Module level import not at top of file
-import os
-import sys
-import unittest
-
-import numpy as np
-import numpy.typing as npt
-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")
-if parent_dir_path not in sys.path:
-    sys.path.append(parent_dir_path)
-
-import geos_posp.processing.geomechanicsCalculatorFunctions as fcts
-from geos.utils import PhysicalConstants
-from geos.utils.algebraFunctions import getAttributeMatrixFromVector
-from geos.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])
-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)
-        )
-
-    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:
-        """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)
-        )
-
-    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)
-
-        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:
-        """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.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:
-        """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:
-        """Test calculation of oedometric Modulus."""
-        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)
-        )
-
-    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)
-        )
-
-    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:
-        """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:
-        """Test calculation of lithostatic Stress."""
-        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:
-        """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:
-        """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:
-        """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])
-
-        self.assertTrue(
-            np.array_equal(np.round(obtained, 2), np.round(expected, 2), equal_nan=True)
-        )
-
-    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:
-        """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:
-        """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:
-        """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:
-        """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:
-        """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:
-        """Test calculation of compressibility from elastic moduli."""
-        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:
-        """Test calculation of shear capacity utilization."""
-        # inputs
-        traction: npt.NDArray[np.float64] = np.copy(effectiveStress[:, :3])
-        # replace nan values for calculation
-        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:
-        """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)
-
-    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
-            )
-        ]
-        expected: tuple[float, float, float] = (1.748, 2.471, 5.281)
-        self.assertSequenceEqual(tuple(obtained), expected)
-
-    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)
-        ]
-        expected: tuple[float, float] = (2.250, 0.606)
-        self.assertSequenceEqual(tuple(obtained), expected)
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/geos-posp/tests/testsMohrCircle.py b/geos-posp/tests/testsMohrCircle.py
deleted file mode 100644
index 83bda4e9..00000000
--- a/geos-posp/tests/testsMohrCircle.py
+++ /dev/null
@@ -1,294 +0,0 @@
-# SPDX-License-Identifier: Apache-2.0
-# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Martin Lemay
-# ruff: noqa: E402 # disable Module level import not at top of file
-import os
-import sys
-import unittest
-
-import numpy as np
-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")
-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
-
-circleId = "12453"
-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)
-
-
-class TestsMohrCircle(unittest.TestCase):
-
-    def test_MohrCircleInit(self: Self) -> None:
-        """Test instanciation of MohrCircle objects."""
-        mohrCircle: MohrCircle = MohrCircle(circleId)
-        mohrCircle.setPrincipalComponents(*principalComponentsExpected)
-
-        self.assertSequenceEqual(circleId, mohrCircle.getCircleId())
-
-        obtained: list[float] = [
-            round(val, 3) for val in mohrCircle.getPrincipalComponents()
-        ]
-        self.assertSequenceEqual(obtained, principalComponentsExpected)
-
-    def test_MohrCircleRadius(self: Self) -> None:
-        """Test the calculation of MohrCircle radius."""
-        mohrCircle: MohrCircle = MohrCircle(circleId)
-        mohrCircle.computePrincipalComponents(stressVector)
-
-        obtained: float = mohrCircle.getCircleRadius()
-        expected: float = (
-            principalComponentsExpected[2] - principalComponentsExpected[0]
-        ) / 2.0
-        self.assertAlmostEqual(obtained, expected, 3)
-
-    def test_MohrCircleCenter(self: Self) -> None:
-        """Test the calculation of MohrCircle center."""
-        mohrCircle: MohrCircle = MohrCircle(circleId)
-        mohrCircle.computePrincipalComponents(stressVector)
-
-        obtained: float = mohrCircle.getCircleCenter()
-        expected: float = (
-            principalComponentsExpected[2] + principalComponentsExpected[0]
-        ) / 2.0
-        self.assertAlmostEqual(obtained, expected, 3)
-
-
-class TestsMohrCoulomb(unittest.TestCase):
-
-    def test_MohrCoulombInit(self: Self) -> None:
-        """Test instanciation of MohrCoulomb objects."""
-        # expected values
-        expectedSlope: float = np.tan(frictionAngle)
-        expectedSigmaMin: float = -rockCohesion / expectedSlope
-        # instantiate object
-        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)
-
-    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)
-        # 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,
-            ]
-        )
-        # instantiate object
-        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())
-
-    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,
-            ]
-        )
-
-        # instantiate object
-        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()
-        )
-
-    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,
-            ]
-        )
-
-        # instantiate object
-        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()
-        )