feat: add principal val and direction output for GM calc

geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py

@@ -137,4 +137,7 @@ def computePrincipalComponents( self: Self, stressVector: npt.NDArray[ np.float6
                 Stress vector must follow GEOS convention (XX, YY, ZZ, YZ, XZ, XY)
         """
         # Get stress principal components
-        self.p3, self.p2, self.p1 = ( computeStressPrincipalComponentsFromStressVector( stressVector ) )
+        vp1, vp2, vp3 = np.unstack(
+            np.swapaxes( computeStressPrincipalComponentsFromStressVector( stressVector.reshape( -1, 6 ) )[ 0 ], 0,
+                         1 ) )
+        self.p1, self.p2, self.p3 = ( vp1[ 0 ], vp2[ 0 ], vp3[ 0 ] )

geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py

@@ -672,21 +672,41 @@ def criticalPorePressure(
     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
+    pca, _ = computeStressPrincipalComponentsFromStressVector( stressVector )
 
     # assertion frictionAngle < np.pi/2., so sin(frictionAngle) != 1
     cohesiveTerm: npt.NDArray[ np.floating[ Any ] ] = ( rockCohesion * np.cos( frictionAngle ) /
                                                         ( 1.0 - np.sin( frictionAngle ) ) )
-    residualTerm: npt.NDArray[ np.floating[ Any ] ] = ( 3.0 * minimumPrincipalStress - maximumPrincipalStress ) / 2.0
+    residualTerm: npt.NDArray[ np.floating[ Any ] ] = ( 3.0 * np.min( pca, axis=1 ) - np.max( pca, axis=1 ) ) / 2.0
 
     return cohesiveTerm + residualTerm
 
 
+def principalAxesAndDirections(
+    stressVector: npt.NDArray[ np.float64 ], ) -> tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]:
+    r"""Getting the principal axes and directions separately by fast eigh decomposition.
+
+    Beware they are ordered in reverse natural ordre so that p1>=p2>=p3.
+
+    Args:
+        stressVector (npt.NDArray[np.float64]): stress vector of Nx6 dims
+            (:math:`\sigma` - Pa).
+
+    Returns:
+        tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ]: tuple containing principal stresses values and directions
+            (p1, p2, p3) and (d1, d2, d3) where p1>=p2>=p3 and d1, d2, d3 are the corresponding eigenvectors.
+
+    """
+    assert stressVector is not None, "Stress vector must be defined"
+    assert stressVector.shape[ 1 ] == 6, "Stress vector must be of size 6."
+
+    principalStresses: npt.NDArray[ np.float64 ] = np.empty( shape=( stressVector.shape[ 0 ], 3 ) )
+    principalDirs: npt.NDArray[ np.float64 ] = np.empty( shape=( stressVector.shape[ 0 ], 3, 3 ) )
+    principalStresses, principalDirs = computeStressPrincipalComponentsFromStressVector( stressVector )
+
+    return principalStresses, principalDirs
+
+
 def criticalPorePressureThreshold( pressure: npt.NDArray[ np.float64 ],
                                    criticalPorePressure: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     r"""Compute the critical pore pressure threshold.
@@ -882,62 +902,49 @@ def shearCapacityUtilization( traction: npt.NDArray[ np.float64 ], rockCohesion:
 
 
 def computeStressPrincipalComponentsFromStressVector(
-    stressVector: npt.NDArray[ np.float64 ], ) -> tuple[ float, float, float ]:
+    stressVector: npt.NDArray[ np.float64 ], ) -> tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]:
     """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
+        tuple[float, float, float]: Principal components sorted in descending
             order.
 
     """
-    assert stressVector.size == 6, "Stress vector dimension is wrong."
+    assert stressVector.shape[ 1 ] == 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 )
+    eVal, eVec = np.linalg.eigh( stressTensor )
+
+    return eVal[ :, ::-1 ], eVec[ :, :, ::-1 ]
 
 
-def computeNormalShearStress( stressTensor: npt.NDArray[ np.float64 ],
-                              directionVector: npt.NDArray[ np.float64 ] ) -> tuple[ float, float ]:
+def computeNormalShearStress(
+        stressTensors: npt.NDArray[ np.float64 ],
+        directionVectors: npt.NDArray[ np.float64 ] ) -> tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]:
     """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
+        stressTensors (npt.NDArray[np.float64]): stress vector.
+        directionVectors: npt.NDArray[ np.float64 ]: the direction vectors of the local frame.
 
     Returns:
-        tuple[float, float]: normal and shear stresses.
+        tuple[npt.NDArray[np.float64] , npt.NDArray[np.float64]]: 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"
+    assert stressTensors.shape[ 0 ] == directionVectors.shape[ 0 ], "First dim should be number of points"
+    assert stressTensors.shape[ 1: ] == ( 3, 3 ), "Stress tensor must be nx3x3 matrix."
+    assert directionVectors.shape[ 1 ] == 3, "Direction vector must have 3 components"
 
     # normalization of direction vector
-    directionVector = directionVector / np.linalg.norm( directionVector )
+    directionVectors = directionVectors / np.linalg.norm( directionVectors, axis=1, keepdims=True )
     # stress vector
-    T: npt.NDArray[ np.float64 ] = np.dot( stressTensor, directionVector )
+    T: npt.NDArray[ np.float64 ] = np.einsum( 'nij,nj->ni', stressTensors, directionVectors )
     # normal stress
-    sigmaN: float = np.dot( T, directionVector )
+    sigmaN: npt.NDArray[ np.float64 ] = np.einsum( 'ni,ni->n', T, directionVectors )
     # shear stress
-    tauVec: npt.NDArray[ np.float64 ] = T - np.dot( sigmaN, directionVector )
-    tau: float = float( np.linalg.norm( tauVec ) )
+    tauVec: npt.NDArray[ np.float64 ] = T - np.einsum( 'n,nj->nj', sigmaN, directionVectors )
+    tau: npt.NDArray[ np.float64 ] = np.linalg.norm( tauVec, axis=1 )
     return ( sigmaN, tau )

geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py

@@ -37,7 +37,7 @@
 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 ] )
+stressVector: npt.NDArray[ np.float64 ] = np.array( [ 1.8, 2.5, 5.2, 0.3, 0.4, 0.1 ] ).reshape( 1, 6 )
 stressTensor: npt.NDArray[ np.float64 ] = getAttributeMatrixFromVector( stressVector )
 
 
@@ -250,28 +250,21 @@ def test_shearCapacityUtilization( self: Self ) -> None:
 
         self.assertSequenceEqual( np.round( obtained, 3 ).flatten().tolist(), expected )
 
-    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 )
+        obtained: npt.NDArray[ np.float64 ] = fcts.computeStressPrincipalComponentsFromStressVector( stressVector )[ 0 ]
+        expected: npt.NDArray[ np.float64 ] = np.array( [ 5.281, 2.471, 1.748 ] ).reshape( 1, -1 )
+        assert np.linalg.norm( obtained - expected ) < 1e-3
 
     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 )
+        directionVector = np.array( [ 1.0, 1.0, 0.0 ] ).reshape( 1, 3 )
+        dv = np.vstack( [ directionVector, directionVector ] )
+        sv = np.vstack( [ stressTensor, stressTensor ] )
+        obtained: tuple[ npt.NDArray[ np.float64 ],
+                         npt.NDArray[ np.float64 ] ] = fcts.computeNormalShearStress( sv, dv )
+        expected: tuple[ list[ float ], list[ float ] ] = ( [ 2.250, 2.250 ], [ 0.606, 0.606 ] )
+        assert np.linalg.norm( [ obt - exp for obt, exp in zip( obtained, expected ) ] ) < 1e-3
 
 
 if __name__ == "__main__":

geos-mesh/src/geos/mesh/utils/genericHelpers.py

@@ -7,12 +7,11 @@
 from typing import ( Iterator, List, Sequence, Any, Union, Tuple )
 
 from vtkmodules.util.numpy_support import ( numpy_to_vtk, vtk_to_numpy )
-from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints, reference, vtkDataArray, vtkLogger, vtkFloatArray
+from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints, reference, vtkLogger
 from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkMultiBlockDataSet, vtkPolyData, vtkDataSet,
                                             vtkDataObject, vtkPlane, vtkCellTypes, vtkIncrementalOctreePointLocator,
                                             VTK_TRIANGLE, vtkSelection, vtkSelectionNode )
-from vtkmodules.vtkFiltersCore import ( vtk3DLinearGridPlaneCutter, vtkPolyDataNormals, vtkPolyDataTangents )
-from vtkmodules.vtkFiltersTexture import vtkTextureMapToPlane
+from vtkmodules.vtkFiltersCore import ( vtk3DLinearGridPlaneCutter, vtkPolyDataNormals )
 from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter, vtkGeometryFilter
 from vtkmodules.vtkFiltersGeneral import vtkDataSetTriangleFilter
 from vtkmodules.util.numpy_support import numpy_to_vtkIdTypeArray
@@ -22,7 +21,7 @@
 from geos.mesh.utils.multiblockHelpers import ( getBlockElementIndexesFlatten, getBlockFromFlatIndex )
 
 from geos.utils.algebraFunctions import ( getAttributeMatrixFromVector, getAttributeVectorFromMatrix )
-from geos.utils.geometryFunctions import ( getChangeOfBasisMatrix, CANONICAL_BASIS_3D )
+from geos.utils.geometryFunctions import ( getChangeOfBasisMatrix, CANONICAL_BASIS_3D, _normalize, _cross )
 from geos.utils.Logger import ( getLogger, Logger, VTKCaptureLog, RegexExceptionFilter )
 from geos.utils.Errors import ( VTKError )
 
@@ -397,7 +396,7 @@ def createVertices( cellPtsCoord: list[ npt.NDArray[ np.float64 ] ],
     return points, cellVertexMapAll
 
 
-def convertAttributeFromLocalToXYZForOneCell(
+def convertAttributeFromLocalToXYZ(
     vector: npt.NDArray[ np.float64 ], localBasisVectors: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ],
                                                                  npt.NDArray[ np.float64 ] ]
 ) -> npt.NDArray[ np.float64 ]:
@@ -421,10 +420,11 @@ def convertAttributeFromLocalToXYZForOneCell(
     transformMatrix: npt.NDArray[ np.float64 ] = getChangeOfBasisMatrix( localBasisVectors, CANONICAL_BASIS_3D )
 
     # Apply transformation
-    arrayXYZ: npt.NDArray[ np.float64 ] = transformMatrix @ matrix3x3 @ transformMatrix.T
+    # arrayXYZ: npt.NDArray[ np.float64 ] = np.einsum( 'nij,njk,nlk->nil', transformMatrix, matrix3x3, transformMatrix )
+    arrayXYZ: npt.NDArray[ np.float64 ] = np.einsum( 'nki,nkl,nlj->nij', transformMatrix, matrix3x3, transformMatrix )
 
     # Convert back to GEOS type attribute and return
-    return getAttributeVectorFromMatrix( arrayXYZ, vector.size )
+    return getAttributeVectorFromMatrix( arrayXYZ, vector.shape )
 
 
 def getNormalVectors( surface: vtkPolyData ) -> npt.NDArray[ np.float64 ]:
@@ -460,20 +460,15 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
     Returns:
         Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: The tangents vectors of the input surface.
     """
-    # Get first tangential component
-    vtkTangents: Union[ vtkFloatArray, None ] = surface.GetCellData().GetTangents()  # type: ignore[no-untyped-call]
-    tangents1: npt.NDArray[ np.float64 ]
-
     try:
-        tangents1 = vtk_to_numpy( vtkTangents )
+        tangents1: npt.NDArray[ np.float64 ] = _normalize( vtk_to_numpy( surface.GetCellData().GetTangents() ) )
+        # Compute second tangential component
+        normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface )
+        tangents2: npt.NDArray[ np.float64 ] = _cross( normals, tangents1 )
+        tangents2 = _normalize( tangents2 )
     except AttributeError as err:
         context: str = f"No tangential attribute found in the mesh. Use the computeTangents function beforehand.\n{ err }"
         raise VTKError( context ) from err
-    else:
-        # Compute second tangential component
-        normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface )
-
-        tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 )
 
     return ( tangents1, tangents2 )
 
@@ -513,7 +508,7 @@ def getLocalBasisVectors(
         surfaceWithTangents: vtkPolyData = computeTangents( surfaceWithNormals, logger )
         tangents = getTangentsVectors( surfaceWithTangents )
 
-    return np.array( ( normals, *tangents ) )
+    return np.stack( [ normals, *tangents ], axis=2 )
 
 
 def computeNormals(
@@ -581,94 +576,27 @@ def computeTangents(
         vtkPolyData: The surface with tangent attribute
     """
     # Need to compute texture coordinates required for tangent calculation
-    surface1: vtkPolyData = computeSurfaceTextureCoordinates( triangulatedSurface, logger )
+    # surface1: vtkPolyData = computeSurfaceTextureCoordinates( triangulatedSurface, logger )
 
     # TODO: triangulate the surface before computation of the tangents if needed.
 
     # Compute tangents
-    vtkErrorLogger: Logger
-    if logger is None:
-        vtkErrorLogger = getLogger( "Compute Surface Tangents vtkError Logger", True )
-    else:
-        vtkErrorLogger = logging.getLogger( f"{ logger.name } vtkError Logger" )
-        vtkErrorLogger.setLevel( logging.INFO )
-        vtkErrorLogger.addHandler( logger.handlers[ 0 ] )
-        vtkErrorLogger.propagate = False
-
-    vtkLogger.SetStderrVerbosity( vtkLogger.VERBOSITY_ERROR )
-    vtkErrorLogger.addFilter( RegexExceptionFilter() )  # will raise VTKError if captured VTK Error
-
-    with VTKCaptureLog() as capturedLog:
-
-        tangentFilter: vtkPolyDataTangents = vtkPolyDataTangents()
-        tangentFilter.SetInputData( surface1 )
-        tangentFilter.ComputeCellTangentsOn()
-        tangentFilter.ComputePointTangentsOff()
-        tangentFilter.Update()
-        surfaceOut: vtkPolyData = tangentFilter.GetOutput()
-
-        capturedLog.seek( 0 )
-        captured = capturedLog.read().decode()
-
-    if captured != "":
-        vtkErrorLogger.error( captured.strip() )
-
-    if surfaceOut is None:
-        raise VTKError( "Something went wrong in VTK calculation." )
+    ptArray: npt.NDArray[ np.float64 ] = vtk_to_numpy( triangulatedSurface.GetPoints().GetData() )
+    # Cooking the indexes
+    ids: list[ list[ int ] ] = [ [
+        triangulatedSurface.GetCell( i ).GetPointIds().GetId( j )
+        for j in range( triangulatedSurface.GetCell( i ).GetNumberOfPoints() )
+    ] for i in range( triangulatedSurface.GetNumberOfCells() ) ]
+    # to honor orientation t1 is always p1-p0 normalized !! Beware to pinched cells
+    array = _normalize( ptArray[ ids ][ :, 1, : ] - ptArray[ ids ][ :, 0, : ] )
 
     # copy tangents attributes into filter input surface because surface attributes
     # are not transferred to the output (bug that should be corrected by Kitware)
-    array: vtkDataArray = surfaceOut.GetCellData().GetTangents()
-    if array is None:
-        raise VTKError( "Attribute Tangents is not in the mesh." )
-
-    surface1.GetCellData().SetTangents( array )
-    surface1.GetCellData().Modified()
-    surface1.Modified()
-    return surface1
-
-
-def computeSurfaceTextureCoordinates(
-    surface: vtkPolyData,
-    logger: Union[ Logger, None ] = None,
-) -> vtkPolyData:
-    """Generate the 2D texture coordinates required for tangent computation. The dataset points are mapped onto a plane generated automatically.
-
-    Args:
-        surface (vtkPolyData): The input surface.
-        logger (Union[Logger, None], optional): A logger to manage the output messages.
-            Defaults to None, an internal logger is used.
-
-    Returns:
-        vtkPolyData: The input surface with generated texture map.
-    """
-    # Need to compute texture coordinates required for tangent calculation
-    vtkErrorLogger: Logger
-    if logger is None:
-        vtkErrorLogger = getLogger( "Compute Surface Texture Coordinates vtkError Logger", True )
-    else:
-        vtkErrorLogger = logging.getLogger( f"{ logger.name } vtkError Logger" )
-        vtkErrorLogger.setLevel( logging.INFO )
-        vtkErrorLogger.addHandler( logger.handlers[ 0 ] )
-        vtkErrorLogger.propagate = False
-
-    vtkLogger.SetStderrVerbosity( vtkLogger.VERBOSITY_ERROR )
-    vtkErrorLogger.addFilter( RegexExceptionFilter() )  # will raise VTKError if captured VTK Error
-
-    with VTKCaptureLog() as capturedLog:
-
-        textureFilter: vtkTextureMapToPlane = vtkTextureMapToPlane()
-        textureFilter.SetInputData( surface )
-        textureFilter.AutomaticPlaneGenerationOn()
-        textureFilter.Update()
-
-        capturedLog.seek( 0 )
-        captured = capturedLog.read().decode()
-
-    if captured != "":
-        vtkErrorLogger.error( captured.strip() )
 
-    return textureFilter.GetOutput()
+    triangulatedSurface.GetCellData().SetTangents( numpy_to_vtk( array ) )
+    triangulatedSurface.GetCellData().Modified()
+    triangulatedSurface.Modified()
+    return triangulatedSurface
 
 
 def extractCellSelection( mesh: vtkUnstructuredGrid, ids: list[ int ] ) -> vtkUnstructuredGrid: