From ab48c808be564a0347e17778c5925e435b64fcc7 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Fri, 3 Apr 2026 18:00:06 +0200
Subject: [PATCH 01/22] adding Nicolvectorizing GMCalculator

---
 .../geomechanicsCalculatorFunctions.py        | 68 ++++++++-----------
 .../testsFunctionsGeomechanicsCalculator.py   | 25 +++----
 .../post_processing/GeomechanicsCalculator.py | 25 ++++++-
 3 files changed, 65 insertions(+), 53 deletions(-)

diff --git a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
index 84f21f0fe..205ca6afb 100644
--- a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
+++ b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
@@ -672,20 +672,27 @@ 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 ]]:
+    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))
+    for i in range( principalStresses.shape[ 0 ] ):
+        principalStresses[i,:], principalDirs[i,:,:] = computeStressPrincipalComponentsFromStressVector( stressVector[ i ] )
+
+    return principalStresses, principalDirs
 
 def criticalPorePressureThreshold( pressure: npt.NDArray[ np.float64 ],
                                    criticalPorePressure: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
@@ -882,7 +889,7 @@ 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:
@@ -893,51 +900,36 @@ def computeStressPrincipalComponentsFromStressVector(
             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,eVec
 
-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
+        stressTensor (npt.NDArray[np.float64]): nx3x3 stress tensor
+        directionVector (npt.NDArray[np.float64]): nx3 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"
+    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: np.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: float = np.linalg.norm( tauVec, axis=1 )
     return ( sigmaN, tau )
diff --git a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
index 7fc8b4b5a..365bebaf9 100644
--- a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
+++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
@@ -11,6 +11,9 @@
 import numpy.typing as npt
 from typing_extensions import Self
 
+
+sys.path.insert(0,"/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-geomechanics/src")
+sys.path.insert(0,"/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-utils/src")
 dir_path = os.path.dirname( os.path.realpath( __file__ ) )
 parent_dir_path = os.path.join( os.path.dirname( dir_path ), "src" )
 if parent_dir_path not in sys.path:
@@ -37,7 +40,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 +253,22 @@ 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 )
+            round(val,3) for val in list( fcts.computeStressPrincipalComponentsFromStressVector( stressVector )[0] )
         ]
         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 )
+        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__":
diff --git a/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py b/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py
index da65b869c..07777b1a9 100644
--- a/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py
+++ b/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py
@@ -10,6 +10,9 @@
 import numpy as np
 import numpy.typing as npt
 
+import sys
+sys.path.insert(0,"/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-processing/src")
+
 import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 
 from geos.mesh.utils.arrayModifiers import createAttribute
@@ -47,6 +50,7 @@
     - Total initial stress, total current stress and total stress ratio
     - Elastic stain
     - Real reservoir stress path and reservoir stress path in oedometric condition
+    - Average stress Principal values and direction
 
 The advanced geomechanics properties computed on the mesh are:
     - Fracture index and threshold
@@ -146,11 +150,16 @@
 RSP_REAL: AttributeEnum = PostProcessingOutputsEnum.RSP_REAL
 RSP_OED: AttributeEnum = PostProcessingOutputsEnum.RSP_OED
 STRESS_EFFECTIVE_RATIO_OED: AttributeEnum = PostProcessingOutputsEnum.STRESS_EFFECTIVE_RATIO_OED
+PRINCIPAL_AXIS_VAL: AttributeEnum = PostProcessingOutputsEnum.PRINCIPAL_AXIS_VAL
+PRINCIPAL_AXIS_DIR_1: AttributeEnum = PostProcessingOutputsEnum.PRINCIPAL_AXIS_DIR_1
+PRINCIPAL_AXIS_DIR_2: AttributeEnum = PostProcessingOutputsEnum.PRINCIPAL_AXIS_DIR_2
+PRINCIPAL_AXIS_DIR_3: AttributeEnum = PostProcessingOutputsEnum.PRINCIPAL_AXIS_DIR_3
 BASIC_PROPERTIES: tuple[ AttributeEnum,
                          ...] = ( BIOT_COEFFICIENT, COMPRESSIBILITY, COMPRESSIBILITY_OED, COMPRESSIBILITY_REAL,
                                   SPECIFIC_GRAVITY, STRESS_EFFECTIVE_RATIO_REAL, STRESS_TOTAL, STRESS_TOTAL_T0,
                                   STRESS_TOTAL_RATIO_REAL, LITHOSTATIC_STRESS, AVERAGE_STRAIN, STRESS_TOTAL_DELTA,
-                                  RSP_REAL, RSP_OED, STRESS_EFFECTIVE_RATIO_OED )
+                                  RSP_REAL, RSP_OED, STRESS_EFFECTIVE_RATIO_OED, PRINCIPAL_AXIS_VAL, 
+                                  PRINCIPAL_AXIS_DIR_1, PRINCIPAL_AXIS_DIR_2, PRINCIPAL_AXIS_DIR_3 )
 
 # Advanced properties:
 CRITICAL_TOTAL_STRESS_RATIO: AttributeEnum = PostProcessingOutputsEnum.CRITICAL_TOTAL_STRESS_RATIO
@@ -423,6 +432,10 @@ class BasicProperties:
         _rspReal: npt.NDArray[ np.float64 ] | None = None
         _rspOed: npt.NDArray[ np.float64 ] | None = None
         _effectiveStressRatioOed: npt.NDArray[ np.float64 ] | None = None
+        _principalAxesVal: npt.NDArray[ np.float64 ] | None = None
+        _principalAxesDir1: npt.NDArray[ np.float64 ] | None = None
+        _principalAxesDir2: npt.NDArray[ np.float64 ] | None = None
+        _principalAxesDir3: npt.NDArray[ np.float64 ] | None = None
 
         @property
         def biotCoefficient( self: Self ) -> npt.NDArray[ np.float64 ] | None:
@@ -896,9 +909,19 @@ def _computeBasicProperties( self: Self ) -> None:
         self._computeEffectiveStressRatioOed()
         self._computeReservoirStressPathOed()
         self._computeReservoirStressPathReal()
+        self._computePrincipalAxesAndDirections()
 
         self.logger.info( "All geomechanics basic properties have been successfully computed." )
         return
+    
+    def _computePrincipalAxesAndDirections( self: Self )->None:
+        """Compute the Principal axis and directions."""
+        self._basicProperties._principalAxesVal, dirs = fcts.computeStressPrincipalComponentsFromStressVector(self._mandatoryProperties._effectiveStress)
+        for i in range(3):
+            self._basicProperties._principalAxesDir1[:,:,i] = dirs[:,i]
+        self.logger.info( "All geomechanics basic properties have been successfully computed." )
+        return
+        
 
     def _computeAdvancedProperties( self: Self ) -> None:
         """Compute the advanced geomechanics properties."""

From 66c36a8d78fa34a6fe607e1758251d4d28fb81dc Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Tue, 7 Apr 2026 21:22:23 +0200
Subject: [PATCH 02/22] vectorizing + Nicola s script

---
 .../geomechanicsCalculatorFunctions.py        | 63 ++++++++++-------
 .../testsFunctionsGeomechanicsCalculator.py   | 23 ++++---
 .../src/geos/mesh/utils/genericHelpers.py     |  9 +--
 .../post_processing/GeomechanicsCalculator.py | 27 +++++---
 .../post_processing/SurfaceGeomechanics.py    | 15 ++--
 .../tests/test_GeomechanicsCalculator.py      |  8 +++
 .../tests/test_MeshToMeshInterpolator.py      |  4 --
 .../src/geos/utils/GeosOutputsConstants.py    |  4 ++
 geos-utils/src/geos/utils/algebraFunctions.py | 69 ++++++++++---------
 .../src/geos/utils/geometryFunctions.py       | 25 +++----
 10 files changed, 143 insertions(+), 104 deletions(-)

diff --git a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
index 205ca6afb..745f98b9f 100644
--- a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
+++ b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py
@@ -677,23 +677,36 @@ def criticalPorePressure(
     # 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 * np.min(pca,axis=1) - np.max(pca,axis=1) ) / 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 ]]:
+    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))
-    for i in range( principalStresses.shape[ 0 ] ):
-        principalStresses[i,:], principalDirs[i,:,:] = computeStressPrincipalComponentsFromStressVector( stressVector[ i ] )
+    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.
@@ -889,47 +902,49 @@ def shearCapacityUtilization( traction: npt.NDArray[ np.float64 ], rockCohesion:
 
 
 def computeStressPrincipalComponentsFromStressVector(
-    stressVector: npt.NDArray[ np.float64 ], ) -> tuple[ npt.NDArray[np.float64], npt.NDArray[np.float64] ]:
+    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.shape[1] == 6, "Stress vector dimension is wrong."
+    assert stressVector.shape[ 1 ] == 6, "Stress vector dimension is wrong."
     stressTensor: npt.NDArray[ np.float64 ] = getAttributeMatrixFromVector( stressVector )
     # sort principal stresses from smallest to largest
-    eVal, eVec =np.linalg.eigh( stressTensor )
+    eVal, eVec = np.linalg.eigh( stressTensor )
+
+    return eVal[ :, ::-1 ], eVec[ :, :, ::-1 ]
 
-    return eVal,eVec
 
-def computeNormalShearStress( stressTensors: npt.NDArray[ np.float64 ],
-                              directionVectors: npt.NDArray[ np.float64 ] ) -> tuple[ npt.NDArray[np.float64], npt.NDArray[np.float64] ]:
+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]): nx3x3 stress tensor
-        directionVector (npt.NDArray[np.float64]): nx3 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 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"
+    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
     directionVectors = directionVectors / np.linalg.norm( directionVectors, axis=1, keepdims=True )
     # stress vector
-    T: npt.NDArray[ np.float64 ] = np.einsum('nij,nj->ni', stressTensors, directionVectors )
+    T: npt.NDArray[ np.float64 ] = np.einsum( 'nij,nj->ni', stressTensors, directionVectors )
     # normal stress
-    sigmaN: np.NDArray[ np.float64 ] = np.einsum('ni,ni->n', T, directionVectors )
+    sigmaN: npt.NDArray[ np.float64 ] = np.einsum( 'ni,ni->n', T, directionVectors )
     # shear stress
-    tauVec: npt.NDArray[ np.float64 ] = T - np.einsum('n,nj->nj', sigmaN ,directionVectors )
-    tau: float = np.linalg.norm( tauVec, axis=1 )
+    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 )
diff --git a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
index 365bebaf9..608f71a6b 100644
--- a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
+++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
@@ -11,9 +11,8 @@
 import numpy.typing as npt
 from typing_extensions import Self
 
-
-sys.path.insert(0,"/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-geomechanics/src")
-sys.path.insert(0,"/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-utils/src")
+# sys.path.insert( 0, "/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-geomechanics/src" )
+# sys.path.insert( 0, "/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-utils/src" )
 dir_path = os.path.dirname( os.path.realpath( __file__ ) )
 parent_dir_path = os.path.join( os.path.dirname( dir_path ), "src" )
 if parent_dir_path not in sys.path:
@@ -40,7 +39,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 ] ).reshape(1,6)
+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 )
 
 
@@ -256,19 +255,21 @@ def test_shearCapacityUtilization( self: Self ) -> None:
     def test_computeStressPrincipalComponentsFromStressVector( self: Self ) -> None:
         """Test calculation of stress principal components from stress vector."""
         obtained: list[ float ] = [
-            round(val,3) for val in list( fcts.computeStressPrincipalComponentsFromStressVector( stressVector )[0] )
+            round( val, 3 )
+            for val in list( fcts.computeStressPrincipalComponentsFromStressVector( stressVector )[ 0 ].flatten() )
         ]
         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 ] ).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
+        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__":
diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index d6c2321ad..b30993bdc 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -421,10 +421,10 @@ 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 )
 
     # 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 ]:
@@ -465,7 +465,7 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
     tangents1: npt.NDArray[ np.float64 ]
 
     try:
-        tangents1 = vtk_to_numpy( vtkTangents )
+        tangents1 = np.einsum('ni,n->ni',vtk_to_numpy( vtkTangents ),1/np.linalg.norm(vtk_to_numpy( vtkTangents ), axis=1))
     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
@@ -474,6 +474,7 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
         normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface )
 
         tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 )
+        tangents2 = np.einsum('ni,n->ni',tangents2, 1/np.linalg.norm(tangents2, axis=1))
 
     return ( tangents1, tangents2 )
 
@@ -513,7 +514,7 @@ def getLocalBasisVectors(
         surfaceWithTangents: vtkPolyData = computeTangents( surfaceWithNormals, logger )
         tangents = getTangentsVectors( surfaceWithTangents )
 
-    return np.array( ( normals, *tangents ) )
+    return np.swapaxes(np.array( ( normals, *tangents ) ), 0, 1)
 
 
 def computeNormals(
diff --git a/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py b/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py
index 07777b1a9..960ddfd58 100644
--- a/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py
+++ b/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py
@@ -11,7 +11,8 @@
 import numpy.typing as npt
 
 import sys
-sys.path.insert(0,"/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-processing/src")
+
+sys.path.insert( 0, "/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-processing/src" )
 
 import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 
@@ -158,7 +159,7 @@
                          ...] = ( BIOT_COEFFICIENT, COMPRESSIBILITY, COMPRESSIBILITY_OED, COMPRESSIBILITY_REAL,
                                   SPECIFIC_GRAVITY, STRESS_EFFECTIVE_RATIO_REAL, STRESS_TOTAL, STRESS_TOTAL_T0,
                                   STRESS_TOTAL_RATIO_REAL, LITHOSTATIC_STRESS, AVERAGE_STRAIN, STRESS_TOTAL_DELTA,
-                                  RSP_REAL, RSP_OED, STRESS_EFFECTIVE_RATIO_OED, PRINCIPAL_AXIS_VAL, 
+                                  RSP_REAL, RSP_OED, STRESS_EFFECTIVE_RATIO_OED, PRINCIPAL_AXIS_VAL,
                                   PRINCIPAL_AXIS_DIR_1, PRINCIPAL_AXIS_DIR_2, PRINCIPAL_AXIS_DIR_3 )
 
 # Advanced properties:
@@ -613,6 +614,14 @@ def getBasicPropertyValue( self: Self, name: str ) -> npt.NDArray[ np.float64 ]
                 return self.rspOed
             elif name == STRESS_EFFECTIVE_RATIO_OED.attributeName:
                 return self.effectiveStressRatioOed
+            elif name == PRINCIPAL_AXIS_VAL.attributeName:
+                return self._principalAxesVal
+            elif name == PRINCIPAL_AXIS_DIR_1.attributeName:
+                return self._principalAxesDir1
+            elif name == PRINCIPAL_AXIS_DIR_2.attributeName:
+                return self._principalAxesDir2
+            elif name == PRINCIPAL_AXIS_DIR_3.attributeName:
+                return self._principalAxesDir3
             else:
                 raise NameError( f"The property { name } is not a basic property." )
 
@@ -913,15 +922,17 @@ def _computeBasicProperties( self: Self ) -> None:
 
         self.logger.info( "All geomechanics basic properties have been successfully computed." )
         return
-    
-    def _computePrincipalAxesAndDirections( self: Self )->None:
+
+    def _computePrincipalAxesAndDirections( self: Self ) -> None:
         """Compute the Principal axis and directions."""
-        self._basicProperties._principalAxesVal, dirs = fcts.computeStressPrincipalComponentsFromStressVector(self._mandatoryProperties._effectiveStress)
-        for i in range(3):
-            self._basicProperties._principalAxesDir1[:,:,i] = dirs[:,i]
+        self._basicProperties._principalAxesVal, dir = fcts.computeStressPrincipalComponentsFromStressVector(
+            self._mandatoryProperties._effectiveStress )
+        self._basicProperties._principalAxesDir1, self._basicProperties._principalAxesDir2, self._basicProperties._principalAxesDir3 = np.unstack(
+            dir, axis=2 )
+        self._attributesToCreate.extend(
+            [ PRINCIPAL_AXIS_VAL, PRINCIPAL_AXIS_DIR_1, PRINCIPAL_AXIS_DIR_2, PRINCIPAL_AXIS_DIR_3 ] )
         self.logger.info( "All geomechanics basic properties have been successfully computed." )
         return
-        
 
     def _computeAdvancedProperties( self: Self ) -> None:
         """Compute the advanced geomechanics properties."""
diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
index 843ac9222..0f69ab39b 100644
--- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
+++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
@@ -359,20 +359,19 @@ def __computeXYZCoordinates(
         Raises:
             ValueError: Error with the shape of attrArray or the computation of the attribute coordinates.
         """
-        attrXYZ: npt.NDArray[ np.float64 ] = np.full_like( attrArray, np.nan )
+        attrXYZ: npt.NDArray[ np.float64 ] = np.zeros_like( attrArray )
 
         # Get all local basis vectors
         localBasis: npt.NDArray[ np.float64 ] = getLocalBasisVectors( self.outputMesh, self.logger )
 
-        for i, cellAttribute in enumerate( attrArray ):
-            if len( cellAttribute ) not in ( 3, 6, 9 ):
-                raise ValueError(
-                    f"Inconsistent number of components for attribute. Expected 3, 6 or 9 but got { len( cellAttribute.shape ) }."
-                )
+        if len( attrArray[0] ) not in ( 3, 6, 9 ):
+            raise ValueError(
+                f"Inconsistent number of components for attribute. Expected 3, 6 or 9 but got { len( attrArray[0].shape ) }."
+            )
 
             # Compute attribute XYZ components
-            cellLocalBasis: npt.NDArray[ np.float64 ] = localBasis[ :, i, : ]
-            attrXYZ[ i ] = convertAttributeFromLocalToXYZForOneCell( cellAttribute, cellLocalBasis )
+        # cellLocalBasis: npt.NDArray[ np.float64 ] = localBasis[ :, i, : ]
+        attrXYZ = convertAttributeFromLocalToXYZForOneCell( attrArray, localBasis )
 
         if not np.any( np.isfinite( attrXYZ ) ):
             raise ValueError( "Attribute new coordinates calculation failed." )
diff --git a/geos-processing/tests/test_GeomechanicsCalculator.py b/geos-processing/tests/test_GeomechanicsCalculator.py
index 78e319185..2f73baf18 100644
--- a/geos-processing/tests/test_GeomechanicsCalculator.py
+++ b/geos-processing/tests/test_GeomechanicsCalculator.py
@@ -8,6 +8,8 @@
 
 from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid
 
+from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridWriter
+
 from geos.mesh.utils.arrayHelpers import isAttributeInObject
 from geos.processing.post_processing.GeomechanicsCalculator import ( GeomechanicsCalculator, BASIC_PROPERTIES,
                                                                      ADVANCED_PROPERTIES, LITHOSTATIC_STRESS )
@@ -36,3 +38,9 @@ def test_GeomechanicsCalculator(
     if computeAdvancedProperties:
         for attribute in ADVANCED_PROPERTIES:
             assert isAttributeInObject( output, attribute.attributeName, attribute.piece )
+
+    w = vtkXMLUnstructuredGridWriter()
+    w.SetInputData( output )
+    w.SetFileName( "/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/tmp/testNicola.vtu" )
+    w.Update()
+    w.Write()
diff --git a/geos-processing/tests/test_MeshToMeshInterpolator.py b/geos-processing/tests/test_MeshToMeshInterpolator.py
index 9a9eb937a..1a63de2ff 100644
--- a/geos-processing/tests/test_MeshToMeshInterpolator.py
+++ b/geos-processing/tests/test_MeshToMeshInterpolator.py
@@ -8,10 +8,6 @@
 from collections import Counter
 import numpy as np
 
-import sys
-
-sys.path.insert( 0, "/data/pau901/SIM_CS/users/jfranc/geosPythonPackages/geos-processing" )
-
 from geos.processing.generic_processing_tools.MeshToMeshInterpolator import MeshToMeshInterpolator
 from vtkmodules.vtkCommonDataModel import vtkDataSet
 from vtkmodules.util.numpy_support import vtk_to_numpy
diff --git a/geos-utils/src/geos/utils/GeosOutputsConstants.py b/geos-utils/src/geos/utils/GeosOutputsConstants.py
index 2eb89314e..c484c04ee 100644
--- a/geos-utils/src/geos/utils/GeosOutputsConstants.py
+++ b/geos-utils/src/geos/utils/GeosOutputsConstants.py
@@ -203,6 +203,10 @@ class PostProcessingOutputsEnum( AttributeEnum ):
     AVERAGE_STRAIN = ( "averageStrain", 6, Piece.CELLS )
     RSP_OED = ( "rsp_oed", 1, Piece.CELLS )
     RSP_REAL = ( "rsp_real", 6, Piece.CELLS )
+    PRINCIPAL_AXIS_VAL = ( "principalAxisVal", 3, Piece.CELLS )
+    PRINCIPAL_AXIS_DIR_1 = ( "principalAxisDir1", 3, Piece.CELLS )
+    PRINCIPAL_AXIS_DIR_2 = ( "principalAxisDir2", 3, Piece.CELLS )
+    PRINCIPAL_AXIS_DIR_3 = ( "principalAxisDir3", 3, Piece.CELLS )
 
     # advanced geomechanical outputs
     CRITICAL_TOTAL_STRESS_RATIO = ( "totalStressRatioCritical_real", 1, Piece.CELLS )
diff --git a/geos-utils/src/geos/utils/algebraFunctions.py b/geos-utils/src/geos/utils/algebraFunctions.py
index a7027f9e1..dbb9e4ac9 100644
--- a/geos-utils/src/geos/utils/algebraFunctions.py
+++ b/geos-utils/src/geos/utils/algebraFunctions.py
@@ -54,37 +54,42 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt
     Returns:
         npt.NDArray[np.float64]: The matrix of attribute values.
     """
-    if attrArray.size not in ( 3, 6, 9 ):
+    if attrArray.shape[ 1 ] not in ( 3, 6, 9 ):
         raise ValueError( "Vectorial attribute must contain 3, 6 or 9 components." )
 
     # diagonal terms
-    matrix: npt.NDArray[ np.float64 ] = np.diagflat( attrArray[ :3 ] )
+    matrices: npt.NDArray[ np.float64 ] = np.empty( shape=( attrArray.shape[ 0 ], 3, 3 ) )
+
+    #diag
+    matrices[ :, 0, 0 ] = attrArray[ :, 0 ]  #
+    matrices[ :, 1, 1 ] = attrArray[ :, 1 ]  #
+    matrices[ :, 2, 2 ] = attrArray[ :, 2 ]  #
 
     # shear stress components
-    if attrArray.size == 6:
-        matrix[ 0, 1 ] = attrArray[ 5 ]  # v1
-        matrix[ 1, 0 ] = attrArray[ 5 ]
+    if attrArray.shape[ 1 ] == 6:
+        matrices[ :, 0, 1 ] = attrArray[ :, 5 ]  # v1
+        matrices[ :, 1, 0 ] = attrArray[ :, 5 ]
 
-        matrix[ 0, 2 ] = attrArray[ 4 ]  # v5
-        matrix[ 2, 0 ] = attrArray[ 4 ]
+        matrices[ :, 0, 2 ] = attrArray[ :, 4 ]  # v5
+        matrices[ :, 2, 0 ] = attrArray[ :, 4 ]
 
-        matrix[ 1, 2 ] = attrArray[ 3 ]  # v4
-        matrix[ 2, 1 ] = attrArray[ 3 ]
+        matrices[ :, 1, 2 ] = attrArray[ :, 3 ]  # v4
+        matrices[ :, 2, 1 ] = attrArray[ :, 3 ]
 
-    elif attrArray.size == 9:
-        matrix[ 0, 1 ] = attrArray[ 5 ]  # v1
-        matrix[ 1, 0 ] = attrArray[ 8 ]  # v9
+    elif attrArray.shape[ 1 ] == 9:
+        matrices[ :, 0, 1 ] = attrArray[ :, 5 ]  # v1
+        matrices[ :, 1, 0 ] = attrArray[ :, 8 ]  # v9
 
-        matrix[ 0, 2 ] = attrArray[ 4 ]  # v5
-        matrix[ 2, 0 ] = attrArray[ 7 ]  # v8
+        matrices[ :, 0, 2 ] = attrArray[ :, 4 ]  # v5
+        matrices[ :, 2, 0 ] = attrArray[ :, 7 ]  # v8
 
-        matrix[ 1, 2 ] = attrArray[ 3 ]  # v4
-        matrix[ 2, 1 ] = attrArray[ 6 ]  # v7
+        matrices[ :, 1, 2 ] = attrArray[ :, 3 ]  # v4
+        matrices[ :, 2, 1 ] = attrArray[ :, 6 ]  # v7
 
-    return matrix
+    return matrices
 
 
-def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], size: int ) -> npt.NDArray[ np.float64 ]:
+def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], shape: tuple[int,int,int] ) -> npt.NDArray[ np.float64 ]:
     r"""Get the vector of attribute values from the matrix.
 
     Matrix to vector conversion is the following:
@@ -133,27 +138,29 @@ def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], size: i
     Returns:
         npt.NDArray[np.float64]: Vector of attribute values.
     """
-    attrArray: npt.NDArray[ np.float64 ] = np.full( size, np.nan )
-    if size not in ( 3, 6, 9 ):
+    attrArray: npt.NDArray[ np.float64 ] = np.zeros( shape )
+    if shape[1] not in ( 3, 6, 9 ):
         raise ValueError( "Requested output size can only be 3, 9 or 6 (symmetrical case)." )
-    if attrMatrix.shape != ( 3, 3 ):
+    if attrMatrix.shape[1:] != ( 3, 3 ):
         raise ValueError( "Input matrix shape must be (3,3)." )
 
     # Diagonal terms
-    attrArray[ :3 ] = np.diag( attrMatrix )
+    attrArray[:, 0 ] = attrMatrix[:,0,0]
+    attrArray[:, 1 ] = attrMatrix[:,1,1]
+    attrArray[:, 2 ] = attrMatrix[:,2,2]
 
     # shear stress components
     if attrArray.size == 6:
-        attrArray[ 3 ] = attrMatrix[ 1, 2 ]
-        attrArray[ 4 ] = attrMatrix[ 0, 2 ]
-        attrArray[ 5 ] = attrMatrix[ 0, 1 ]
+        attrArray[:, 3 ] = attrMatrix[:, 1, 2 ]
+        attrArray[:, 4 ] = attrMatrix[:, 0, 2 ]
+        attrArray[:, 5 ] = attrMatrix[:, 0, 1 ]
 
     elif attrArray.size == 9:
-        attrArray[ 3 ] = attrMatrix[ 1, 2 ]
-        attrArray[ 4 ] = attrMatrix[ 0, 2 ]
-        attrArray[ 5 ] = attrMatrix[ 0, 1 ]
-        attrArray[ 6 ] = attrMatrix[ 2, 1 ]
-        attrArray[ 7 ] = attrMatrix[ 2, 0 ]
-        attrArray[ 8 ] = attrMatrix[ 1, 0 ]
+        attrArray[:, 3 ] = attrMatrix[:, 1, 2 ]
+        attrArray[:, 4 ] = attrMatrix[:, 0, 2 ]
+        attrArray[:, 5 ] = attrMatrix[:, 0, 1 ]
+        attrArray[:, 6 ] = attrMatrix[:, 2, 1 ]
+        attrArray[:, 7 ] = attrMatrix[:, 2, 0 ]
+        attrArray[:, 8 ] = attrMatrix[:, 1, 0 ]
 
     return attrArray
diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index eb1843758..f3ccc4419 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -9,16 +9,12 @@
 
 __doc__ = """Functions to permform geometry calculations."""
 
-CANONICAL_BASIS_3D: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ],
-                           npt.NDArray[ np.float64 ] ] = ( np.array( [ 1.0, 0.0, 0.0 ] ), np.array( [ 0.0, 1.0, 0.0 ] ),
-                                                           np.array( [ 0.0, 0.0, 1.0 ] ) )
+CANONICAL_BASIS_3D: npt.NDArray[ np.float64 ] = np.array( [ [ 1.0, 0.0, 0.0 ],[ 0.0, 1.0, 0.0 ],[ 0.0, 0.0, 1.0 ] ] ) 
 
 
 def getChangeOfBasisMatrix(
-    basisFrom: tuple[ npt.NDArray[ np.floating[ Any ] ], npt.NDArray[ np.floating[ Any ] ],
-                      npt.NDArray[ np.floating[ Any ] ] ],
-    basisTo: tuple[ npt.NDArray[ np.floating[ Any ] ], npt.NDArray[ np.floating[ Any ] ],
-                    npt.NDArray[ np.floating[ Any ] ] ],
+    basisFrom:  npt.NDArray[ np.floating[ Any ] ],
+    basisTo: npt.NDArray[ np.floating[ Any ] ],
 ) -> Any:
     """Get the change of basis matrix from basis basisFrom to basis basisTo.
 
@@ -36,19 +32,20 @@ def getChangeOfBasisMatrix(
     Returns:
         npt.NDArray[np.floating[Any]]: Change of basis matrix.
     """
-    assert ( basisFrom[ 0 ].size == basisFrom[ 1 ].size ) and ( basisFrom[ 0 ].size == basisFrom[ 2 ].size ), (
+    assert ( basisFrom.shape[1] == basisFrom.shape[2]) , (
         "Origin space vectors must have the same size." )
 
-    assert ( basisTo[ 0 ].size == basisTo[ 1 ].size ) and ( basisTo[ 0 ].size == basisTo[ 2 ].size ), (
+    basisTo = np.repeat(basisTo[None,:],basisFrom.shape[0], axis=0) 
+    assert ( basisTo.shape[1] == basisTo.shape[2]) , (
         "Destination space vectors must have the same size." )
 
     # build the matrices where columns are the vectors of the bases
-    B = np.transpose( np.array( basisFrom ) )
-    C = np.transpose( np.array( basisTo ) )
-    # compute the inverse of C
-    C1: npt.NDArray[ np.floating[ Any ] ] = np.linalg.inv( C )
+    # B = np.transpose( np.array( basisFrom ) )
+    # C = np.transpose( np.array( basisTo ) )
+    # no need to compute the inverse of C as it is orthonormal checked - transpose is enough 
+    assert np.linalg.norm( np.einsum('nji,nij->nij',basisTo,basisTo) - np.repeat(np.eye(3)[None,:],basisFrom.shape[0], axis=0) ) < 1e-6
     # get the change of basis matrix
-    return np.dot( C1, B )
+    return np.einsum('nji,nij->nij',basisTo,basisFrom)
 
 
 def computeCoordinatesInNewBasis( vec: npt.NDArray[ np.floating[ Any ] ],

From d9f9830e19878a736e03ec3a54c594934a347acc Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Wed, 8 Apr 2026 13:15:30 +0200
Subject: [PATCH 03/22] fixes

---
 geos-mesh/src/geos/mesh/utils/genericHelpers.py               | 2 +-
 geos-mesh/tests/test_genericHelpers_normalAndTangents.py      | 4 ++--
 .../geos/processing/post_processing/SurfaceGeomechanics.py    | 4 ++--
 geos-utils/src/geos/utils/algebraFunctions.py                 | 2 +-
 geos-utils/src/geos/utils/geometryFunctions.py                | 2 +-
 5 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index b30993bdc..2e183a8bd 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -397,7 +397,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 ]:
diff --git a/geos-mesh/tests/test_genericHelpers_normalAndTangents.py b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
index 5cfc3f276..0e41adda9 100644
--- a/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
+++ b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
@@ -14,7 +14,7 @@
 
 from geos.mesh.utils.genericHelpers import ( computeSurfaceTextureCoordinates, computeTangents, computeNormals,
                                              getLocalBasisVectors, getTangentsVectors, getNormalVectors,
-                                             convertAttributeFromLocalToXYZForOneCell )
+                                             convertAttributeFromLocalToXYZ )
 
 from geos.utils.Errors import VTKError
 
@@ -297,6 +297,6 @@ def test_convertAttributesToXYZ( testcase: AttributeConversionTestCase ) -> None
     """Test the conversion of one cell attribute from local to canonic basis."""
     localAttr: npt.NDArray[ np.float64 ] = testcase.vector
 
-    attrXYZ: npt.NDArray[ np.float64 ] = convertAttributeFromLocalToXYZForOneCell(
+    attrXYZ: npt.NDArray[ np.float64 ] = convertAttributeFromLocalToXYZ(
         localAttr, ( testcase.normal, testcase.tangent1, testcase.tangent2 ) )
     assert np.allclose( attrXYZ, testcase.expectedVectorXYZ, rtol=1e-6 )
diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
index 0f69ab39b..a48d32ace 100644
--- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
+++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
@@ -13,7 +13,7 @@
 
 from geos.mesh.utils.arrayModifiers import createAttribute
 from geos.mesh.utils.arrayHelpers import ( getArrayInObject, getAttributeSet, isAttributeInObject )
-from geos.mesh.utils.genericHelpers import ( getLocalBasisVectors, convertAttributeFromLocalToXYZForOneCell )
+from geos.mesh.utils.genericHelpers import ( getLocalBasisVectors, convertAttributeFromLocalToXYZ )
 import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts
 from geos.utils.pieceEnum import Piece
 from geos.utils.Logger import ( getLogger, Logger, CountVerbosityHandler, isHandlerInLogger, getLoggerHandlerType )
@@ -371,7 +371,7 @@ def __computeXYZCoordinates(
 
             # Compute attribute XYZ components
         # cellLocalBasis: npt.NDArray[ np.float64 ] = localBasis[ :, i, : ]
-        attrXYZ = convertAttributeFromLocalToXYZForOneCell( attrArray, localBasis )
+        attrXYZ = convertAttributeFromLocalToXYZ( attrArray, localBasis )
 
         if not np.any( np.isfinite( attrXYZ ) ):
             raise ValueError( "Attribute new coordinates calculation failed." )
diff --git a/geos-utils/src/geos/utils/algebraFunctions.py b/geos-utils/src/geos/utils/algebraFunctions.py
index dbb9e4ac9..54f397bdb 100644
--- a/geos-utils/src/geos/utils/algebraFunctions.py
+++ b/geos-utils/src/geos/utils/algebraFunctions.py
@@ -58,7 +58,7 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt
         raise ValueError( "Vectorial attribute must contain 3, 6 or 9 components." )
 
     # diagonal terms
-    matrices: npt.NDArray[ np.float64 ] = np.empty( shape=( attrArray.shape[ 0 ], 3, 3 ) )
+    matrices: npt.NDArray[ np.float64 ] = np.zeros( shape=( attrArray.shape[ 0 ], 3, 3 ) )
 
     #diag
     matrices[ :, 0, 0 ] = attrArray[ :, 0 ]  #
diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index f3ccc4419..0095dcc3c 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -33,7 +33,7 @@ def getChangeOfBasisMatrix(
         npt.NDArray[np.floating[Any]]: Change of basis matrix.
     """
     assert ( basisFrom.shape[1] == basisFrom.shape[2]) , (
-        "Origin space vectors must have the same size." )
+        f"Origin space vectors must have the same size. shape: {basisFrom.shape}" )
 
     basisTo = np.repeat(basisTo[None,:],basisFrom.shape[0], axis=0) 
     assert ( basisTo.shape[1] == basisTo.shape[2]) , (

From 76233f789d361334d460ed30c398ea3d5e2ac68a Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Wed, 8 Apr 2026 13:22:41 +0200
Subject: [PATCH 04/22] linting

---
 .../src/geos/mesh/utils/genericHelpers.py     |  9 ++---
 .../post_processing/SurfaceGeomechanics.py    |  2 +-
 geos-utils/src/geos/utils/algebraFunctions.py | 33 ++++++++++---------
 .../src/geos/utils/geometryFunctions.py       | 19 ++++++-----
 4 files changed, 33 insertions(+), 30 deletions(-)

diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index 2e183a8bd..a379d4f7c 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -421,7 +421,7 @@ def convertAttributeFromLocalToXYZ(
     transformMatrix: npt.NDArray[ np.float64 ] = getChangeOfBasisMatrix( localBasisVectors, CANONICAL_BASIS_3D )
 
     # Apply transformation
-    arrayXYZ: npt.NDArray[ np.float64 ] = np.einsum('nij,njk,nlk->nil', transformMatrix,matrix3x3,transformMatrix )
+    arrayXYZ: npt.NDArray[ np.float64 ] = np.einsum( 'nij,njk,nlk->nil', transformMatrix, matrix3x3, transformMatrix )
 
     # Convert back to GEOS type attribute and return
     return getAttributeVectorFromMatrix( arrayXYZ, vector.shape )
@@ -465,7 +465,8 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
     tangents1: npt.NDArray[ np.float64 ]
 
     try:
-        tangents1 = np.einsum('ni,n->ni',vtk_to_numpy( vtkTangents ),1/np.linalg.norm(vtk_to_numpy( vtkTangents ), axis=1))
+        tangents1 = np.einsum( 'ni,n->ni', vtk_to_numpy( vtkTangents ),
+                               1 / np.linalg.norm( vtk_to_numpy( vtkTangents ), axis=1 ) )
     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
@@ -474,7 +475,7 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
         normals: npt.NDArray[ np.float64 ] = getNormalVectors( surface )
 
         tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 )
-        tangents2 = np.einsum('ni,n->ni',tangents2, 1/np.linalg.norm(tangents2, axis=1))
+        tangents2 = np.einsum( 'ni,n->ni', tangents2, 1 / np.linalg.norm( tangents2, axis=1 ) )
 
     return ( tangents1, tangents2 )
 
@@ -514,7 +515,7 @@ def getLocalBasisVectors(
         surfaceWithTangents: vtkPolyData = computeTangents( surfaceWithNormals, logger )
         tangents = getTangentsVectors( surfaceWithTangents )
 
-    return np.swapaxes(np.array( ( normals, *tangents ) ), 0, 1)
+    return np.swapaxes( np.array( ( normals, *tangents ) ), 0, 1 )
 
 
 def computeNormals(
diff --git a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
index a48d32ace..335ceeab9 100644
--- a/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
+++ b/geos-processing/src/geos/processing/post_processing/SurfaceGeomechanics.py
@@ -364,7 +364,7 @@ def __computeXYZCoordinates(
         # Get all local basis vectors
         localBasis: npt.NDArray[ np.float64 ] = getLocalBasisVectors( self.outputMesh, self.logger )
 
-        if len( attrArray[0] ) not in ( 3, 6, 9 ):
+        if len( attrArray[ 0 ] ) not in ( 3, 6, 9 ):
             raise ValueError(
                 f"Inconsistent number of components for attribute. Expected 3, 6 or 9 but got { len( attrArray[0].shape ) }."
             )
diff --git a/geos-utils/src/geos/utils/algebraFunctions.py b/geos-utils/src/geos/utils/algebraFunctions.py
index 54f397bdb..988346f14 100644
--- a/geos-utils/src/geos/utils/algebraFunctions.py
+++ b/geos-utils/src/geos/utils/algebraFunctions.py
@@ -89,7 +89,8 @@ def getAttributeMatrixFromVector( attrArray: npt.NDArray[ np.float64 ], ) -> npt
     return matrices
 
 
-def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], shape: tuple[int,int,int] ) -> npt.NDArray[ np.float64 ]:
+def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ],
+                                  shape: tuple[ int, int, int ] ) -> npt.NDArray[ np.float64 ]:
     r"""Get the vector of attribute values from the matrix.
 
     Matrix to vector conversion is the following:
@@ -129,7 +130,7 @@ def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], shape:
 
     Args:
         attrMatrix (npt.NDArray[np.float64]): Matrix of attribute values.
-        size (int): Size of the final vector. Values accepted are 3, 9 or 6.
+        shape (tuple[int, int, int]): Shape of the output vector (n, 3), (n, 6) or (n, 9).
 
     Raises:
         ValueError: The output vector size requested can only be 3, 9 or 6.
@@ -139,28 +140,28 @@ def getAttributeVectorFromMatrix( attrMatrix: npt.NDArray[ np.float64 ], shape:
         npt.NDArray[np.float64]: Vector of attribute values.
     """
     attrArray: npt.NDArray[ np.float64 ] = np.zeros( shape )
-    if shape[1] not in ( 3, 6, 9 ):
+    if shape[ 1 ] not in ( 3, 6, 9 ):
         raise ValueError( "Requested output size can only be 3, 9 or 6 (symmetrical case)." )
-    if attrMatrix.shape[1:] != ( 3, 3 ):
+    if attrMatrix.shape[ 1: ] != ( 3, 3 ):
         raise ValueError( "Input matrix shape must be (3,3)." )
 
     # Diagonal terms
-    attrArray[:, 0 ] = attrMatrix[:,0,0]
-    attrArray[:, 1 ] = attrMatrix[:,1,1]
-    attrArray[:, 2 ] = attrMatrix[:,2,2]
+    attrArray[ :, 0 ] = attrMatrix[ :, 0, 0 ]
+    attrArray[ :, 1 ] = attrMatrix[ :, 1, 1 ]
+    attrArray[ :, 2 ] = attrMatrix[ :, 2, 2 ]
 
     # shear stress components
     if attrArray.size == 6:
-        attrArray[:, 3 ] = attrMatrix[:, 1, 2 ]
-        attrArray[:, 4 ] = attrMatrix[:, 0, 2 ]
-        attrArray[:, 5 ] = attrMatrix[:, 0, 1 ]
+        attrArray[ :, 3 ] = attrMatrix[ :, 1, 2 ]
+        attrArray[ :, 4 ] = attrMatrix[ :, 0, 2 ]
+        attrArray[ :, 5 ] = attrMatrix[ :, 0, 1 ]
 
     elif attrArray.size == 9:
-        attrArray[:, 3 ] = attrMatrix[:, 1, 2 ]
-        attrArray[:, 4 ] = attrMatrix[:, 0, 2 ]
-        attrArray[:, 5 ] = attrMatrix[:, 0, 1 ]
-        attrArray[:, 6 ] = attrMatrix[:, 2, 1 ]
-        attrArray[:, 7 ] = attrMatrix[:, 2, 0 ]
-        attrArray[:, 8 ] = attrMatrix[:, 1, 0 ]
+        attrArray[ :, 3 ] = attrMatrix[ :, 1, 2 ]
+        attrArray[ :, 4 ] = attrMatrix[ :, 0, 2 ]
+        attrArray[ :, 5 ] = attrMatrix[ :, 0, 1 ]
+        attrArray[ :, 6 ] = attrMatrix[ :, 2, 1 ]
+        attrArray[ :, 7 ] = attrMatrix[ :, 2, 0 ]
+        attrArray[ :, 8 ] = attrMatrix[ :, 1, 0 ]
 
     return attrArray
diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index 0095dcc3c..a456c5b7f 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -9,11 +9,11 @@
 
 __doc__ = """Functions to permform geometry calculations."""
 
-CANONICAL_BASIS_3D: npt.NDArray[ np.float64 ] = np.array( [ [ 1.0, 0.0, 0.0 ],[ 0.0, 1.0, 0.0 ],[ 0.0, 0.0, 1.0 ] ] ) 
+CANONICAL_BASIS_3D: npt.NDArray[ np.float64 ] = np.array( [ [ 1.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 1.0 ] ] )
 
 
 def getChangeOfBasisMatrix(
-    basisFrom:  npt.NDArray[ np.floating[ Any ] ],
+    basisFrom: npt.NDArray[ np.floating[ Any ] ],
     basisTo: npt.NDArray[ np.floating[ Any ] ],
 ) -> Any:
     """Get the change of basis matrix from basis basisFrom to basis basisTo.
@@ -32,20 +32,21 @@ def getChangeOfBasisMatrix(
     Returns:
         npt.NDArray[np.floating[Any]]: Change of basis matrix.
     """
-    assert ( basisFrom.shape[1] == basisFrom.shape[2]) , (
+    assert ( basisFrom.shape[ 1 ] == basisFrom.shape[ 2 ] ), (
         f"Origin space vectors must have the same size. shape: {basisFrom.shape}" )
 
-    basisTo = np.repeat(basisTo[None,:],basisFrom.shape[0], axis=0) 
-    assert ( basisTo.shape[1] == basisTo.shape[2]) , (
-        "Destination space vectors must have the same size." )
+    basisTo = np.repeat( basisTo[ None, : ], basisFrom.shape[ 0 ], axis=0 )
+    assert ( basisTo.shape[ 1 ] == basisTo.shape[ 2 ] ), ( "Destination space vectors must have the same size." )
 
     # build the matrices where columns are the vectors of the bases
     # B = np.transpose( np.array( basisFrom ) )
     # C = np.transpose( np.array( basisTo ) )
-    # no need to compute the inverse of C as it is orthonormal checked - transpose is enough 
-    assert np.linalg.norm( np.einsum('nji,nij->nij',basisTo,basisTo) - np.repeat(np.eye(3)[None,:],basisFrom.shape[0], axis=0) ) < 1e-6
+    # no need to compute the inverse of C as it is orthonormal checked - transpose is enough
+    assert np.linalg.norm(
+        np.einsum( 'nji,nij->nij', basisTo, basisTo ) -
+        np.repeat( np.eye( 3 )[ None, : ], basisFrom.shape[ 0 ], axis=0 ) ) < 1e-6
     # get the change of basis matrix
-    return np.einsum('nji,nij->nij',basisTo,basisFrom)
+    return np.einsum( 'nji,nij->nij', basisTo, basisFrom )
 
 
 def computeCoordinatesInNewBasis( vec: npt.NDArray[ np.floating[ Any ] ],

From 20758f49a96ccff38f5f63caf272c71028107ca5 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Thu, 9 Apr 2026 13:08:43 +0200
Subject: [PATCH 05/22] refactor unused - vectorized others

---
 .../src/geos/mesh/utils/genericHelpers.py     |   4 +-
 .../pre_processing/MeshQualityEnhanced.py     |   3 +
 .../src/geos/utils/geometryFunctions.py       | 403 +++++++++---------
 geos-utils/tests/testsGeometryFunctions.py    | 370 ++++++++--------
 4 files changed, 391 insertions(+), 389 deletions(-)

diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index a379d4f7c..468f3ecde 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -466,7 +466,7 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
 
     try:
         tangents1 = np.einsum( 'ni,n->ni', vtk_to_numpy( vtkTangents ),
-                               1 / np.linalg.norm( vtk_to_numpy( vtkTangents ), axis=1 ) )
+                               1 / np.linalg.norm( vtk_to_numpy( vtkTangents ), axis=1) )
     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
@@ -474,7 +474,7 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
         # 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 )
+        # tangents2: npt.NDArray[ np.float64 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 )
         tangents2 = np.einsum( 'ni,n->ni', tangents2, 1 / np.linalg.norm( tangents2, axis=1 ) )
 
     return ( tangents1, tangents2 )
diff --git a/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py b/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py
index 84ad960f4..7a23acff5 100644
--- a/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py
+++ b/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py
@@ -778,6 +778,7 @@ def _computeSquishIndex( self: Self ) -> None:
         copyData: vtkUnstructuredGrid = vtkUnstructuredGrid()
         copyData.ShallowCopy( self._outputMesh )
         points: vtkPoints = copyData.GetPoints()
+
         for c in range( copyData.GetNumberOfCells() ):
             cell: vtkCell = copyData.GetCell( c )
             # Applies only to polyhedra
@@ -796,6 +797,7 @@ def _computeSquishIndex( self: Self ) -> None:
                 for i in range( ptsIdsList.GetNumberOfIds() ):
                     ptsIds.InsertNextValue( ptsIdsList.GetId( i ) )
                 faceCenter: npt.NDArray[ np.float64 ] = self._getCellCenter( face, ptsIds, points )
+                # TODO use vtkFilter instead of reinvening the wheel a 1000th time !!!
                 faceNormal: npt.NDArray[ np.float64 ] = self._getNormalVector( points, face )
 
                 vec: npt.NDArray[ np.float64 ] = cellCenter - faceCenter
@@ -811,6 +813,7 @@ def _computeSquishIndex( self: Self ) -> None:
         cellArrays.Modified()
         self._outputMesh.Modified()
 
+
     def _getCellCenter( self: Self,
                         cell: vtkCell,
                         ptsIds: Optional[ vtkIdTypeArray ] = None,
diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index a456c5b7f..5984d25a4 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -9,12 +9,25 @@
 
 __doc__ = """Functions to permform geometry calculations."""
 
-CANONICAL_BASIS_3D: npt.NDArray[ np.float64 ] = np.array( [ [ 1.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 1.0 ] ] )
+CANONICAL_BASIS_3D: npt.NDArray[ np.float64 ] = np.array( [ [ 1.0, 0.0, 0.0 ],[ 0.0, 1.0, 0.0 ],[ 0.0, 0.0, 1.0 ] ] ) 
 
+# for batch cross product
+EPS = np.zeros((3, 3, 3), dtype=int)
+EPS[0,1,2] = EPS[1,2,0] = EPS[2,0,1] =  1
+EPS[0,2,1] = EPS[2,1,0] = EPS[1,0,2] = -1
+
+def _normBasis( basis: npt.NDArray[np.float64] ):
+    """Norm and orthonormalize basis vector wise"""
+    Q, _ = np.linalg.qr( basis )
+    det = np.linalg.det( Q )
+    Q[ det < 0 ] *= -1
+
+    return Q
+    
 
 def getChangeOfBasisMatrix(
-    basisFrom: npt.NDArray[ np.floating[ Any ] ],
-    basisTo: npt.NDArray[ np.floating[ Any ] ],
+    basisFrom:  npt.NDArray[ np.float64 ],
+    basisTo: npt.NDArray[ np.float64 ],
 ) -> Any:
     """Get the change of basis matrix from basis basisFrom to basis basisTo.
 
@@ -26,202 +39,209 @@ def getChangeOfBasisMatrix(
     C is then Vc = P.Vb
 
     Args:
-        basisFrom (tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]): Origin basis
-        basisTo (tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]): Destination basis
+        basisFrom (tuple[npt.NDArray[np.float64 npt.NDArray[np.float64 npt.NDArray[np.float64): Origin basis
+        basisTo (tuple[npt.NDArray[np.float64 npt.NDArray[np.float64 npt.NDArray[np.float64): Destination basis
 
     Returns:
-        npt.NDArray[np.floating[Any]]: Change of basis matrix.
+        npt.NDArray[np.float64 Change of basis matrix.
     """
-    assert ( basisFrom.shape[ 1 ] == basisFrom.shape[ 2 ] ), (
+    basisFrom = _normBasis( basisFrom )
+    basisTo = _normBasis( basisTo )
+
+    assert ( basisFrom.shape[1] == basisFrom.shape[2]) , (
         f"Origin space vectors must have the same size. shape: {basisFrom.shape}" )
 
-    basisTo = np.repeat( basisTo[ None, : ], basisFrom.shape[ 0 ], axis=0 )
-    assert ( basisTo.shape[ 1 ] == basisTo.shape[ 2 ] ), ( "Destination space vectors must have the same size." )
+    if len(basisTo.shape) < len(basisFrom.shape):
+        basisTo = np.repeat(basisTo[None,:],basisFrom.shape[0], axis=0) 
+    assert ( basisTo.shape[1] == basisTo.shape[2]) , (
+        "Destination space vectors must have the same size." )
 
     # build the matrices where columns are the vectors of the bases
     # B = np.transpose( np.array( basisFrom ) )
     # C = np.transpose( np.array( basisTo ) )
     # no need to compute the inverse of C as it is orthonormal checked - transpose is enough
     assert np.linalg.norm(
-        np.einsum( 'nji,nij->nij', basisTo, basisTo ) -
+        np.einsum( 'nlj,nkj->nlk', basisTo, basisTo ) -
         np.repeat( np.eye( 3 )[ None, : ], basisFrom.shape[ 0 ], axis=0 ) ) < 1e-6
     # get the change of basis matrix
-    return np.einsum( 'nji,nij->nij', basisTo, basisFrom )
-
-
-def computeCoordinatesInNewBasis( vec: npt.NDArray[ np.floating[ Any ] ],
-                                  changeOfBasisMatrix: npt.NDArray[ np.floating[ Any ] ] ) -> Any:
-    """Computes the coordinates of a matrix from a basis B in the new basis B'.
-
-    Args:
-        vec (npt.NDArray[np.floating[Any]]): Vector to compute the new coordinates
-        changeOfBasisMatrix (npt.NDArray[np.floating[Any]]): Change of basis matrix
-
-    Returns:
-        npt.NDArray[np.floating[Any]]: The new coordinates of the matrix in the basis
-        B'.
-    """
-    assert ( vec.size == changeOfBasisMatrix.shape[ 1 ] ), """The size of the vector
-    must be equal to the number of columns of and change of basis matrix."""
-    return np.dot( changeOfBasisMatrix, vec )
-
-
-def computePlaneFrom3Points(
-    pt1: npt.NDArray[ np.floating[ Any ] ],
-    pt2: npt.NDArray[ np.floating[ Any ] ],
-    pt3: npt.NDArray[ np.floating[ Any ] ],
-) -> tuple[ float, float, float, float ]:
-    """Compute the 4 coefficients of a plane equation.
-
-    Let's defined a plane from the following equation: ax + by + cz + d = 0.
-    This function determines a, b, c, d from 3 points of the plane.
-
-    Args:
-        pt1 (npt.NDArray[np.floating[Any]]): 1st point of the plane.
-        pt2 (npt.NDArray[np.floating[Any]]): 2nd point of the plane.
-        pt3 (npt.NDArray[np.floating[Any]]): 3rd point of the plane.
-
-    Returns:
-        tuple[float, float, float, float]: Tuple of the 4 coefficients.
-    """
-    # plane vectors
-    v1: npt.NDArray[ np.floating[ Any ] ] = pt2 - pt1
-    v2: npt.NDArray[ np.floating[ Any ] ] = pt3 - pt1
-    # normal vector
-    normal: npt.NDArray[ np.floating[ Any ] ] = np.cross( v1, v2 )
-    assert np.linalg.norm( normal ) != 0, "Vectors of the plane must not be colinear."
-    # first 3 coefficients of the plane equation
-    a, b, c = normal
-    # last coefficient of the plane equation
-    d: float = -np.dot( normal, pt1 )
-    return a, b, c, d
-
-
-def getCellSideAgainstPlane(
-    cellPtsCoords: npt.NDArray[ np.floating[ Any ] ],
-    planePt: npt.NDArray[ np.floating[ Any ] ],
-    planeNormal: npt.NDArray[ np.floating[ Any ] ],
-) -> bool:
-    """Get the side of input cell against input plane.
-
-    Input plane is defined by a point on it and its normal vector.
-
-    Args:
-        cellPtsCoords (npt.NDArray[np.floating[Any]]): Coordinates of the vertices of
-            the cell to get the side.
-        planePt (npt.NDArray[np.floating[Any]]): Point on the plane.
-        planeNormal (npt.NDArray[np.floating[Any]]): Normal vector to the plane.
-
-    Returns:
-        bool: True if the cell is in the direction of the normal vector,
-        False otherwise.
-    """
-    assert ( len( cellPtsCoords ) > 1 ), "The list of points must contains more than one element"
-    ptCenter: npt.NDArray[ np.floating[ Any ] ] = np.mean( cellPtsCoords, axis=0 )
-    return getPointSideAgainstPlane( ptCenter, planePt, planeNormal )
-
-
-def getPointSideAgainstPlane(
-    ptCoords: npt.NDArray[ np.floating[ Any ] ],
-    planePt: npt.NDArray[ np.floating[ Any ] ],
-    planeNormal: npt.NDArray[ np.floating[ Any ] ],
-) -> bool:
-    """Get the side of input point against input plane.
-
-    Input plane is defined by a point on it and its normal vector.
-
-    Args:
-        ptCoords (npt.NDArray[np.floating[Any]]): Coordinates of the point to get
-            the side.
-        planePt (npt.NDArray[np.floating[Any]]): Point on the plane.
-        planeNormal (npt.NDArray[np.floating[Any]]): Normal vector to the plane.
-
-    Returns:
-        bool: True if the point is in the direction of the normal vector,
-        False otherwise.
-    """
-    assert ptCoords.size == 3, "Point coordinates must have 3 components"
-    assert planeNormal.size == 3, "Plane normal vector must have 3 components"
-    assert planePt.size == 3, "Plane point must have 3 components"
-    vec: npt.NDArray[ np.floating[ Any ] ] = ptCoords - planePt
-    dot: float = np.dot( planeNormal, vec )
-    assert abs( dot ) > EPSILON, "The point is on the plane."
-    return dot > 0
-
-
-def computeAngleFromPoints( pt1: npt.NDArray[ np.float64 ], pt2: npt.NDArray[ np.float64 ],
-                            pt3: npt.NDArray[ np.float64 ] ) -> float:
-    """Compute angle from 3 points.
-
-    Args:
-        pt1 (npt.NDArray[np.float64]): First point
-        pt2 (npt.NDArray[np.float64]): Second point
-        pt3 (npt.NDArray[np.float64]): Third point
-
-    Returns:
-        float: Angle
-    """
-    # compute vectors
-    vec1: npt.NDArray[ np.float64 ] = pt1 - pt2
-    vec2: npt.NDArray[ np.float64 ] = pt3 - pt2
-    return computeAngleFromVectors( vec1, vec2 )
-
-
-def computeAngleFromVectors(
-    vec1: npt.NDArray[ np.float64 ],
-    vec2: npt.NDArray[ np.float64 ],
-) -> float:
-    """Compute angle from 2 vectors.
-
-    Args:
-        vec1 (npt.NDArray[np.float64]): First vector
-        vec2 (npt.NDArray[np.float64]): Second vector
-
-    Returns:
-        float: angle
-    """
-    assert abs( np.linalg.norm( vec1 ) ) > 0., "First vector cannot be null"
-    assert abs( np.linalg.norm( vec2 ) ) > 0., "Second vector cannot be null"
-    # normalization
-    vec1_norm: npt.NDArray[ np.float64 ] = vec1 / np.linalg.norm( vec1 )
-    vec2_norm: npt.NDArray[ np.float64 ] = vec2 / np.linalg.norm( vec2 )
-
-    # Use normalized vectors for dot product
-    dot: float = np.dot( vec1_norm, vec2_norm )
-
-    # Clamp to valid range for arccos
-    dot = np.clip( dot, -1.0, 1.0 )
-
-    # Handle 2D vs 3D cases
-    if vec1.size == 2 and vec2.size == 2:
-        # For 2D, use cross product to determine orientation
-        cross: float = vec1_norm[ 0 ] * vec2_norm[ 1 ] - vec1_norm[ 1 ] * vec2_norm[ 0 ]
-        angle: float = np.arccos( dot )
-        return angle if cross >= 0 else 2.0 * np.pi - angle
-    else:
-        # For 3D, return angle in [0, π]
-        return np.arccos( dot )
-
-
-def computeCosineFromVectors(
-    vec1: npt.NDArray[ np.float64 ],
-    vec2: npt.NDArray[ np.float64 ],
-) -> float:
-    """Compute cosine from 2 vectors.
-
-    Args:
-        vec1 (npt.NDArray[np.float64]): First vector
-        vec2 (npt.NDArray[np.float64]): Second vector
-
-    Returns:
-        float: Cosine
-    """
-    assert abs( np.linalg.norm( vec1 ) ) > 0., "First vector cannot be null"
-    assert abs( np.linalg.norm( vec2 ) ) > 0., "Second vector cannot be null"
-    # normalization
-    vec1_norm: npt.NDArray[ np.float64 ] = vec1 / np.linalg.norm( vec1 )
-    vec2_norm: npt.NDArray[ np.float64 ] = vec2 / np.linalg.norm( vec2 )
-    return np.dot( vec1_norm, vec2_norm )
+    return np.einsum( 'nlj,nkj->nlk', basisTo, basisFrom )
+
+
+# def computeCoordinatesInNewBasis( vec: npt.NDArray[ np.float64 ],
+#                                   changeOfBasisMatrix: npt.NDArray[ np.float64 ] ) -> Any:
+#     """Computes the coordinates of a matrix from a basis B in the new basis B'.
+
+#     Args:
+#         vec (npt.NDArray[np.float64: Vector to compute the new coordinates
+#         changeOfBasisMatrix (npt.NDArray[np.float64: Change of basis matrix
+
+#     Returns:
+#         npt.NDArray[np.float64 The new coordinates of the matrix in the basis
+#         B'.
+#     """
+#     assert ( vec.shape[ 1 ] == changeOfBasisMatrix.shape[ 1 ] ), """The size of the vector
+#     must be equal to the number of columns of and change of basis matrix."""
+#     return np.einsum('nij,ni->ni' ,changeOfBasisMatrix, vec )
+
+
+# def computePlaneFrom3Points(
+#     pt1: npt.NDArray[ np.float64 ],
+#     pt2: npt.NDArray[ np.float64 ],
+#     pt3: npt.NDArray[ np.float64 ],
+# ) -> tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]:
+#     """Compute the 4 coefficients of a plane equation.
+
+#     Let's defined a plane from the following equation: ax + by + cz + d = 0.
+#     This function determines a, b, c, d from 3 points of the plane.
+
+#     Args:
+#         pt1 (npt.NDArray[np.float64: 1st point of the plane.
+#         pt2 (npt.NDArray[np.float64: 2nd point of the plane.
+#         pt3 (npt.NDArray[np.float64: 3rd point of the plane.
+
+#     Returns:
+#         tuple[float, float, float, float]: Tuple of the 4 coefficients.
+#     """
+#     assert pt1.shape[0] == pt2.shape[0] and pt2.shape[0] == pt3.shape[0]
+
+#     # plane vectors
+#     v1: npt.NDArray[ np.float64 ] = pt2 - pt1
+#     v2: npt.NDArray[ np.float64 ] = pt3 - pt1
+#     # normal vector
+#     normal: npt.NDArray[ np.float64 ] = np.einsum('ijk,nj,nk->ni', EPS, v1, v2)
+#     assert np.linalg.norm( normal ) != 0, "Vectors of the plane must not be colinear."
+#     # first 3 coefficients of the plane equation
+#     a, b, c = np.unstack( normal, axis=1 )
+#     # last coefficient of the plane equation
+#     d: float = -np.einsum('ni,ni->n', normal, pt1 )
+#     return a, b, c, d
+
+
+# def getCellSideAgainstPlane(
+#     cellPtsCoords: npt.NDArray[ np.float64 ],
+#     planePt: npt.NDArray[ np.float64 ],
+#     planeNormal: npt.NDArray[ np.float64 ],
+# ) -> bool:
+#     """Get the side of input cell against input plane.
+
+#     Input plane is defined by a point on it and its normal vector.
+
+#     Args:
+#         cellPtsCoords (npt.NDArray[np.float64: Coordinates of the vertices of
+#             the cell to get the side.
+#         planePt (npt.NDArray[np.float64: Point on the plane.
+#         planeNormal (npt.NDArray[np.float64: Normal vector to the plane.
+
+#     Returns:
+#         bool: True if the cell is in the direction of the normal vector,
+#         False otherwise.
+#     """
+#     assert ( len( cellPtsCoords.shape[0] ) > 1 ), "The list of points must contains more than one element"
+#     ptCenter: npt.NDArray[ np.float64 ] = np.mean( cellPtsCoords, axis=0 )
+#     return getPointSideAgainstPlane( ptCenter, planePt, planeNormal )
+
+
+# def getPointSideAgainstPlane(
+#     ptCoords: npt.NDArray[ np.float64 ],
+#     planePt: npt.NDArray[ np.float64 ],
+#     planeNormal: npt.NDArray[ np.float64 ],
+# ) -> bool:
+#     """Get the side of input point against input plane.
+
+#     Input plane is defined by a point on it and its normal vector.
+
+#     Args:
+#         ptCoords (npt.NDArray[np.float64: Coordinates of the point to get
+#             the side.
+#         planePt (npt.NDArray[np.float64: Point on the plane.
+#         planeNormal (npt.NDArray[np.float64: Normal vector to the plane.
+
+#     Returns:
+#         bool: True if the point is in the direction of the normal vector,
+#         False otherwise.
+#     """
+#     assert ptCoords.shape[1] == 3, "Point coordinates must have 3 components"
+#     assert planeNormal.shape[1] == 3, "Plane normal vector must have 3 components"
+#     assert planePt.shape[1] == 3, "Plane point must have 3 components"
+#     vec: npt.NDArray[ np.float64 ] = ptCoords - planePt
+#     dot: float = np.einsum('ni,ni->n', planeNormal, vec )
+#     assert np.linalg.norm( dot ) > EPSILON, "The point is on the plane."
+#     return dot > 0
+
+
+# def computeAngleFromPoints( pt1: npt.NDArray[ np.float64 ], pt2: npt.NDArray[ np.float64 ],
+#                             pt3: npt.NDArray[ np.float64 ] ) -> float:
+#     """Compute angle from 3 points.
+
+#     Args:
+#         pt1 (npt.NDArray[np.float64]): First point
+#         pt2 (npt.NDArray[np.float64]): Second point
+#         pt3 (npt.NDArray[np.float64]): Third point
+
+#     Returns:
+#         float: Angle
+#     """
+#     # compute vectors
+#     vec1: npt.NDArray[ np.float64 ] = pt1 - pt2
+#     vec2: npt.NDArray[ np.float64 ] = pt3 - pt2
+#     return computeAngleFromVectors( vec1, vec2 )
+
+
+# def computeAngleFromVectors(
+#     vec1: npt.NDArray[ np.float64 ],
+#     vec2: npt.NDArray[ np.float64 ],
+# ) -> float:
+#     """Compute angle from 2 vectors.
+
+#     Args:
+#         vec1 (npt.NDArray[np.float64]): First vector
+#         vec2 (npt.NDArray[np.float64]): Second vector
+
+#     Returns:
+#         float: angle
+#     """
+#     assert np.where( np.abs( np.linalg.norm( vec1, axis=0 ) ) ) == 0., "First vector cannot be null"
+#     assert np.where( np.abs( np.linalg.norm( vec2, axis=0 ) ) ) == 0., "Second vector cannot be null"
+#     # normalization
+#     vec1_norm: npt.NDArray[ np.float64 ] = np.einsum( 'ni,nj->ni', vec1, 1. / np.linalg.norm( vec1, axis=1, keepdims=True ) )
+#     vec2_norm: npt.NDArray[ np.float64 ] = np.einsum( 'ni,nj->ni' , vec2, 1. / np.linalg.norm( vec2, axis=1, keepdims=True ) )
+
+#     # Use normalized vectors for dot product
+#     dot: float = np.einsum('ni,nj->n', vec1_norm, vec2_norm)
+
+#     # Clamp to valid range for arccos
+#     dot = np.clip( dot, -1.0, 1.0 )
+
+#     # Handle 2D vs 3D cases
+#     if vec1.shape[1] == 2 and vec2.shape[1] == 2:
+#         # For 2D, use cross product to determine orientation
+#         cross: float = np.einsum('ijk,nj,nk->ni', EPS, vec1_norm, vec2_norm)
+#         angle: float = np.arccos( dot )
+#         return angle if cross >= 0 else 2.0 * np.pi - angle
+#     else:
+#         # For 3D, return angle in [0, π]
+#         return np.arccos( dot )
+
+
+# def computeCosineFromVectors(
+#     vec1: npt.NDArray[ np.float64 ],
+#     vec2: npt.NDArray[ np.float64 ],
+# ) -> float:
+#     """Compute cosine from 2 vectors.
+
+#     Args:
+#         vec1 (npt.NDArray[np.float64]): First vector
+#         vec2 (npt.NDArray[np.float64]): Second vector
+
+#     Returns:
+#         float: Cosine
+#     """
+#     assert abs( np.linalg.norm( vec1 ) ) > 0., "First vector cannot be null"
+#     assert abs( np.linalg.norm( vec2 ) ) > 0., "Second vector cannot be null"
+#     # normalization
+#     vec1_norm: npt.NDArray[ np.float64 ] = vec1 / np.linalg.norm( vec1 )
+#     vec2_norm: npt.NDArray[ np.float64 ] = vec2 / np.linalg.norm( vec2 )
+#     return np.dot( vec1_norm, vec2_norm )
 
 
 def computeNormalFromPoints( pt1: npt.NDArray[ np.float64 ], pt2: npt.NDArray[ np.float64 ],
@@ -236,9 +256,12 @@ def computeNormalFromPoints( pt1: npt.NDArray[ np.float64 ], pt2: npt.NDArray[ n
     Returns:
         npt.NDArray[np.float64]: Normal vector coordinates
     """
+    assert pt1.shape[0] == pt2.shape[0] and pt2.shape[0] == pt3.shape[0]
     # compute vectors
     vec1: npt.NDArray[ np.float64 ] = pt1 - pt2
     vec2: npt.NDArray[ np.float64 ] = pt3 - pt2
+    assert np.all( np.abs( np.linalg.norm( vec1, axis=1 ) ) > 0.), "First and second points must be different"
+    assert np.all( np.abs( np.linalg.norm( vec2, axis=1 ) ) > 0.), "First and second points must be different"
     return computeNormalFromVectors( vec1, vec2 )
 
 
@@ -255,9 +278,7 @@ def computeNormalFromVectors(
     Returns:
         npt.NDArray[np.float64]: Normal vector coordinates
     """
-    assert abs( np.linalg.norm( vec1 ) ) > 0., "First and second points must be different"
-    assert abs( np.linalg.norm( vec2 ) ) > 0., "Second and third points must be different"
     # normalization
-    vec1_norm = vec1 / np.linalg.norm( vec1 )
-    vec2_norm = vec2 / np.linalg.norm( vec2 )
-    return np.cross( vec1_norm, vec2_norm )
+    vec1_norm = np.einsum('ni,n->ni',vec1, 1./ np.linalg.norm( vec1, axis=1))
+    vec2_norm = np.einsum('ni,n->ni', vec2, 1./ np.linalg.norm( vec2, axis=1))
+    return np.einsum('ijk,nj,nk->ni', EPS , vec1_norm, vec2_norm )
diff --git a/geos-utils/tests/testsGeometryFunctions.py b/geos-utils/tests/testsGeometryFunctions.py
index dbb83717f..d91fdbe1f 100644
--- a/geos-utils/tests/testsGeometryFunctions.py
+++ b/geos-utils/tests/testsGeometryFunctions.py
@@ -1,6 +1,6 @@
 # SPDX-License-Identifier: Apache-2.0
 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
-# SPDX-FileContributor: Martin Lemay
+# SPDX-FileContributor: Martin Lemay, Jacques Franc
 # ruff: noqa: E402 # disable Module level import not at top of file
 import pytest
 import numpy as np
@@ -10,233 +10,211 @@
 from itertools import combinations
 import geos.utils.geometryFunctions as fcts
 
-basisCanon: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ],
-                   npt.NDArray[ np.float64 ] ] = ( np.array( [ 1.0, 0.0, 0.0 ] ), np.array( [ 0.0, 1.0, 0.0 ] ),
-                                                   np.array( [ 0.0, 0.0, 1.0 ] ) )
+basisCanon: npt.NDArray[ np.float64 ] = np.array( [[[ 1.0, 0.0, 0.0 ] ,[ 0.0, 1.0, 0.0 ], 
+                                            [ 0.0, 0.0, 1.0 ]]] )
 # destination basis according to canonic coordinates
-basisTo0: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ],
-                 npt.NDArray[ np.float64 ] ] = ( np.array( [ 2.0, 3.0, 0.0 ] ), np.array( [ 4.0, 5.0, 0.0 ] ),
-                                                 np.array( [ 0.0, 0.0, 1.0 ] ) )
-basisTo1: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ],
-                 npt.NDArray[ np.float64 ] ] = ( np.array( [ 1.0, 1.0, 0.0 ] ), np.array( [ 0.0, 1.0, 1.0 ] ),
-                                                 np.array( [ 1.0, 0.0, 1.0 ] ) )
-basisTo2: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ],
-                 npt.NDArray[ np.float64 ] ] = ( np.array( [ 0.0, 2.0, 0.0 ] ), np.array( [ -2.0, 0.0, 0.0 ] ),
-                                                 np.array( [ 0.0, 0.0, 2.0 ] ) )
+basisTo0: npt.NDArray[ np.float64 ] = np.array( [[ [1.0, 1.0, 0.0 ] ,[ -1.0, 1.0, 0.0 ], 
+                                                  [ 0.0, 0.0, 1.0 ] ]] )
+basisTo1: npt.NDArray[ np.float64 ] = np.array( [[ [1.0, -1.0, 0.0 ] ,[ 1.0, 1.0, 0.0 ], 
+                                                 [ 0.0, 0.0, 2.0 ] ]])
+# basisTo2: npt.NDArray[ np.float64 ] = np.array( [[ [0.0, 2.0, 0.0 ] ,[ -2.0, 0.0, 0.0 ],
+#                                                  [0.0, 0.0, 2.0 ] ]])
 
 
 def test_getChangeOfBasisMatrixToCanonic() -> None:
     """Test change of basis matrix using canonic basis."""
     obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisCanon )
     # matrix where the columns are the vectors
-    expected: npt.NDArray[ np.float64 ] = np.transpose( np.array( basisTo0 ) )
-    assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ),
-                           equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}"
-
-
-def test_getChangeOfBasisMatrix() -> None:
-    """Test change of basis matrix format from basis vectors."""
+    expected: npt.NDArray[ np.float64 ] = fcts._normBasis(basisTo0).transpose(0, 2, 1)
+    assert np.linalg.norm( obtained-expected ) < 10e-12 , f"Expected array is {np.round( expected, 2 ).tolist()}"
+    # 
     obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 )
-    expected: npt.NDArray[ np.float64 ] = np.array( [ [ 2.5, 4.5, -0.5 ], [ 0.5, 0.5, 0.5 ], [ -0.5, -0.5, 0.5 ] ] )
-    assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ),
-                           equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}"
-
-
-def test_computeCoordinatesInNewBasisIdentity() -> None:
-    """Test calculation of coordinates of a vector in the same basis."""
-    vec: npt.NDArray[ np.float64 ] = np.array( [ 2.0, 3.0, 4.0 ] )
-
-    # get change of basis matrix
-    changeOfBasisMatrix: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisCanon, basisCanon )
-    obtained: npt.NDArray[ np.float64 ] = fcts.computeCoordinatesInNewBasis( vec, changeOfBasisMatrix )
-    expected: npt.NDArray[ np.float64 ] = vec
-    assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ),
-                           equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}"
-
-
-def test_computeCoordinatesInNewBasis() -> None:
-    """Test calculation of coordinates of a vector in another basis."""
-    vec: npt.NDArray[ np.float64 ] = np.array( [ 2.0, 3.0, 4.0 ] )
-
-    # get change of basis matrix and the inverse
-    changeOfBasisMatrix: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 )
-    obtained = fcts.computeCoordinatesInNewBasis( vec, changeOfBasisMatrix )
-    expected: npt.NDArray[ np.float64 ] = np.array( [ 16.5, 4.5, -0.5 ] )
-    assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ),
-                           equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}"
-
-
-def test_computePlaneFrom3Points() -> None:
-    """Test calculation of plane coefficients from 3 points."""
-    pt1: npt.NDArray[ np.float64 ] = np.array( [ 1.0, 2.0, 1.0 ] )
-    pt2: npt.NDArray[ np.float64 ] = np.array( [ 1.0, 1.0, 2.0 ] )
-    pt3: npt.NDArray[ np.float64 ] = np.array( [ 3.0, 2.0, 2.0 ] )
-    obtained: tuple[ float, float, float, float ] = fcts.computePlaneFrom3Points( pt1, pt2, pt3 )
-    expected: tuple[ float, float, float, float ] = ( -1.0, 2.0, 2.0, -5.0 )
-    assert obtained == expected, f"Expected tuple is {expected}"
-
-
-def test_getPointSideAgainstPlaneAssertion() -> None:
-    """Test get point side against a plane."""
-    planePt: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 0.0 ] )
-
-    # assertion error - Point on the plane
-    planeNormal: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 1.0 ] )
-    with pytest.raises( AssertionError ):
-        fcts.getPointSideAgainstPlane( planePt, planePt, planeNormal )
-
-
-# yapf: disable
-listPtsCoords_all = ( [ np.array( [ 0.5, 0.5, 0.5 ] ), np.array( [ 0.5, 0.5, -0.5 ] ) ],
-                      [ np.array( [ 0.5, 0.5, 0.5 ] ), np.array( [ -0.5, 0.5, 0.5 ] ) ] )
-planePt_all = ( np.array( [ 0.0, 0.0, 0.0 ] ), np.array( [ 0.0, 0.0, 0.0 ] ) )
-planeNormal_all = ( np.array( [ 0.0, 0.0, 1.0 ] ), np.array( [ 1.0, 0.0, 0.0 ] ) )
-expected_all = ( ( True, False ), ( True, False ) )
-# yapf: enable
-
-
-@dataclass( frozen=True )
-class TestCasePointSideAgainstPlane:
-    """Test case."""
-    __test__ = False
-    #: list of points
-    listPtsCoords: list[ npt.NDArray[ np.float64 ] ]
-    #: plane point
-    planePt: npt.NDArray[ np.float64 ]
-    #: plane normal
-    planeNormal: npt.NDArray[ np.float64 ]
-    #: expected result
-    expected: list[ bool ]
-
-
-def __generate_PointSideAgainstPlane_test_data() -> Iterator[ TestCasePointSideAgainstPlane ]:
-    """Generate test cases.
-
-    Yields:
-        Iterator[ TestCase ]: Iterator on test cases
-    """
-    for listPtsCoords, planePt, planeNormal, expected in zip( listPtsCoords_all,
-                                                              planePt_all,
-                                                              planeNormal_all,
-                                                              expected_all,
-                                                              strict=True ):
-        yield TestCasePointSideAgainstPlane( listPtsCoords, planePt, planeNormal, list( expected ) )
-
-
-ids: tuple[ str, str ] = ( "Horizontal plane", "Vertical plane" )
-
-
-@pytest.mark.parametrize( "test_case", __generate_PointSideAgainstPlane_test_data(), ids=ids )
-def test_getPointSideAgainstPlane( test_case: TestCasePointSideAgainstPlane ) -> None:
-    """Test get point side against a plane."""
-    obtained: list[ bool ] = []
-    for ptCoords in test_case.listPtsCoords:
-        side: bool = fcts.getPointSideAgainstPlane( ptCoords, test_case.planePt, test_case.planeNormal )
-        obtained += [ side ]
-    assert obtained == test_case.expected, f"Expected tuple is {test_case.expected}"
-
-
-def test_getCellSideAgainstPlaneRandom() -> None:
-    """Test get cell side against a plane."""
-    # random plane
-    planePt: npt.NDArray[ np.float64 ] = np.array( [ 125.58337, 1386.0465, -2782.502 ] )
-    listCellPtsCoords: list[ npt.NDArray[ np.float64 ] ] = [
-        np.array( [
-            [ 135.49551, 1374.7644, -2786.884 ],
-            [ 125.58337, 1376.7441, -2782.502 ],
-            [ 132.19525, 1382.2516, -2771.0508 ],
-            [ 125.58337, 1386.0465, -2782.502 ],
-        ] ),
-        np.array( [
-            [ 111.9148, 1377.0265, -2764.875 ],
-            [ 132.19525, 1382.2516, -2771.0508 ],
-            [ 125.58337, 1376.7441, -2782.502 ],
-            [ 125.58337, 1386.0465, -2782.502 ],
-        ] ),
-    ]
-    planeNormal: npt.NDArray[ np.float64 ] = np.array( [ 0.8660075, 0.0, -0.5000311 ] )
-    obtained: list[ bool ] = []
-    for cellPtsCoords in listCellPtsCoords:
-        side: bool = fcts.getCellSideAgainstPlane( cellPtsCoords, planePt, planeNormal )
-        obtained += [ side ]
-    expected: list[ bool ] = [ True, False ]
-    assert obtained == expected, f"Expected tuple is {expected}"
-
-
-pts_all: tuple[ npt.NDArray[ np.float64 ], ...] = (
-    np.array( [ 0.0, 0.0, 0.0 ] ),
-    np.array( [ 1.0, 0.0, 0.0 ] ),
-    np.array( [ 0.5, 0.5, 0.0 ] ),
-    np.array( [ -0.5, 0.5, 0.0 ] ),
-    np.array( [ -0.5, -0.5, 0.0 ] ),
-    np.array( [ 0.5, -0.5, -1.0 ] ),
-)
+    # matrix where the columns are the vectors
+    expected: npt.NDArray[ np.float64 ] = np.array([ [0., -1., 0.],[1., 0., 0.],[0., 0., 1.] ]) 
+    assert np.linalg.norm( obtained-expected ) < 10e-12 , f"Expected array is {np.round( expected, 2 ).tolist()}"
+
+
+# def test_computeCoordinatesInNewBasis() -> None:
+#     """Test calculation of coordinates of a vector in another basis."""
+#     vec: npt.NDArray[ np.float64 ] = np.array( [[ 3.0, -3.0, 0.0 ], [0.0, 0.0, 50.0]] )
+
+#     # get change of basis matrix and the inverse
+#     changeOfBasisMatrix: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 )
+#     obtained = fcts.computeCoordinatesInNewBasis( vec, changeOfBasisMatrix )
+#     expected: npt.NDArray[ np.float64 ] = np.array( [[ -3.0, -3.0, 0. ], [0.0, 0.0, 50.0]] )
+#     assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ),
+#                            equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}"
+
+
+# def test_computePlaneFrom3Points() -> None:
+#     """Test calculation of plane coefficients from 3 points."""
+#     pt1: npt.NDArray[ np.float64 ] = np.array( [[ 1.0, 2.0, 1.0 ]] )
+#     pt2: npt.NDArray[ np.float64 ] = np.array( [[ 1.0, 1.0, 2.0 ]] )
+#     pt3: npt.NDArray[ np.float64 ] = np.array( [[ 3.0, 2.0, 2.0 ]] )
+#     obtained: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]   = fcts.computePlaneFrom3Points( pt1, pt2, pt3 )
+#     expected: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]  = ( np.array([-1.0]), np.array([2.0]), np.array([2.0]), np.array([-5.0]) )
+#     assert obtained == expected, f"Expected tuple is {expected}"
+
+
+# def test_getPointSideAgainstPlaneAssertion() -> None:
+#     """Test get point side against a plane."""
+#     planePt: npt.NDArray[ np.float64 ] = np.array( [[ 0.0, 0.0, 0.0 ]] )
+
+#     # assertion error - Point on the plane
+#     planeNormal: npt.NDArray[ np.float64 ] = np.array( [[ 0.0, 0.0, 1.0 ]] )
+#     with pytest.raises( AssertionError ):
+#         fcts.getPointSideAgainstPlane( planePt, planePt, planeNormal )
+
+
+# # yapf: disable
+# listPtsCoords_all = ( [ np.array( [ 0.5, 0.5, 0.5 ] ), np.array( [ 0.5, 0.5, -0.5 ] ) ],
+#                       [ np.array( [ 0.5, 0.5, 0.5 ] ), np.array( [ -0.5, 0.5, 0.5 ] ) ] )
+# planePt_all = ( np.array( [ 0.0, 0.0, 0.0 ] ), np.array( [ 0.0, 0.0, 0.0 ] ) )
+# planeNormal_all = ( np.array( [ 0.0, 0.0, 1.0 ] ), np.array( [ 1.0, 0.0, 0.0 ] ) )
+# expected_all = ( ( True, False ), ( True, False ) )
+# # yapf: enable
+
+
+# @dataclass( frozen=True )
+# class TestCasePointSideAgainstPlane:
+#     """Test case."""
+#     __test__ = False
+#     #: list of points
+#     listPtsCoords: list[ npt.NDArray[ np.float64 ] ]
+#     #: plane point
+#     planePt: npt.NDArray[ np.float64 ]
+#     #: plane normal
+#     planeNormal: npt.NDArray[ np.float64 ]
+#     #: expected result
+#     expected: list[ bool ]
+
+
+# def __generate_PointSideAgainstPlane_test_data() -> Iterator[ TestCasePointSideAgainstPlane ]:
+#     """Generate test cases.
+
+#     Yields:
+#         Iterator[ TestCase ]: Iterator on test cases
+#     """
+#     for listPtsCoords, planePt, planeNormal, expected in zip( listPtsCoords_all,
+#                                                               planePt_all,
+#                                                               planeNormal_all,
+#                                                               expected_all,
+#                                                               strict=True ):
+#         yield TestCasePointSideAgainstPlane( listPtsCoords, planePt, planeNormal, list( expected ) )
+
+
+# ids: tuple[ str, str ] = ( "Horizontal plane", "Vertical plane" )
+
+
+# @pytest.mark.parametrize( "test_case", __generate_PointSideAgainstPlane_test_data(), ids=ids )
+# def test_getPointSideAgainstPlane( test_case: TestCasePointSideAgainstPlane ) -> None:
+#     """Test get point side against a plane."""
+#     obtained: list[ bool ] = []
+#     for ptCoords in test_case.listPtsCoords:
+#         side: bool = fcts.getPointSideAgainstPlane( ptCoords, test_case.planePt, test_case.planeNormal )
+#         obtained += [ side ]
+#     assert obtained == test_case.expected, f"Expected tuple is {test_case.expected}"
+
+
+# def test_getCellSideAgainstPlaneRandom() -> None:
+#     """Test get cell side against a plane."""
+#     # random plane
+#     planePt: npt.NDArray[ np.float64 ] = np.array( [ 125.58337, 1386.0465, -2782.502 ] )
+#     listCellPtsCoords: list[ npt.NDArray[ np.float64 ] ] = [
+#         np.array( [
+#             [ 135.49551, 1374.7644, -2786.884 ],
+#             [ 125.58337, 1376.7441, -2782.502 ],
+#             [ 132.19525, 1382.2516, -2771.0508 ],
+#             [ 125.58337, 1386.0465, -2782.502 ],
+#         ] ),
+#         np.array( [
+#             [ 111.9148, 1377.0265, -2764.875 ],
+#             [ 132.19525, 1382.2516, -2771.0508 ],
+#             [ 125.58337, 1376.7441, -2782.502 ],
+#             [ 125.58337, 1386.0465, -2782.502 ],
+#         ] ),
+#     ]
+#     planeNormal: npt.NDArray[ np.float64 ] = np.array( [ 0.8660075, 0.0, -0.5000311 ] )
+#     obtained: list[ bool ] = []
+#     for cellPtsCoords in listCellPtsCoords:
+#         side: bool = fcts.getCellSideAgainstPlane( cellPtsCoords, planePt, planeNormal )
+#         obtained += [ side ]
+#     expected: list[ bool ] = [ True, False ]
+#     assert obtained == expected, f"Expected tuple is {expected}"
+
+
+# pts_all: tuple[ npt.NDArray[ np.float64 ], ...] = (
+#     np.array( [ 0.0, 0.0, 0.0 ] ),
+#     np.array( [ 1.0, 0.0, 0.0 ] ),
+#     np.array( [ 0.5, 0.5, 0.0 ] ),
+#     np.array( [ -0.5, 0.5, 0.0 ] ),
+#     np.array( [ -0.5, -0.5, 0.0 ] ),
+#     np.array( [ 0.5, -0.5, -1.0 ] ),
+# )
 # yapf: disable
-angleExp_all: tuple[ float, ...] = (
-    0.78539816, 0.32175055, 0.32175055, 1.15026200, 0.78539816, 2.0e-8, 1.04719755,
-    0.78539816, 0.61547971, 1.04719755, 2.35619449, 1.57079633, 1.04719755, 1.24904577,
-    0.75204009, 0.83548187, 1.57079633, 0.95531662, 1.04719755, 1.57079633
-)
-# yapf: enable
+# angleExp_all: tuple[ float, ...] = (
+#     0.78539816, 0.32175055, 0.32175055, 1.15026200, 0.78539816, 2.0e-8, 1.04719755,
+#     0.78539816, 0.61547971, 1.04719755, 2.35619449, 1.57079633, 1.04719755, 1.24904577,
+#     0.75204009, 0.83548187, 1.57079633, 0.95531662, 1.04719755, 1.57079633
+# )
+# # yapf: enable
 
 
-@dataclass( frozen=True )
-class TestCaseAngle:
-    """Test case."""
-    __test__ = False
-    pt1: npt.NDArray[ np.float64 ]
-    pt2: npt.NDArray[ np.float64 ]
-    pt3: npt.NDArray[ np.float64 ]
-    angleExp: float
+# @dataclass( frozen=True )
+# class TestCaseAngle:
+#     """Test case."""
+#     __test__ = False
+#     pt1: npt.NDArray[ np.float64 ]
+#     pt2: npt.NDArray[ np.float64 ]
+#     pt3: npt.NDArray[ np.float64 ]
+#     angleExp: float
 
 
-def __generate_Angle_test_data() -> Iterator[ TestCaseAngle ]:
-    """Generate test cases.
+# def __generate_Angle_test_data() -> Iterator[ TestCaseAngle ]:
+#     """Generate test cases.
 
-    Yields:
-        Iterator[ TestCase ]: Iterator on test cases
-    """
-    for pts, angle in zip( list( combinations( pts_all, 3 ) ), angleExp_all, strict=True ):
-        yield TestCaseAngle( pts[ 0 ], pts[ 1 ], pts[ 2 ], angle )
+#     Yields:
+#         Iterator[ TestCase ]: Iterator on test cases
+#     """
+#     for pts, angle in zip( list( combinations( pts_all, 3 ) ), angleExp_all, strict=True ):
+#         yield TestCaseAngle( pts[ 0 ], pts[ 1 ], pts[ 2 ], angle )
 
 
-@pytest.mark.parametrize( "test_case", __generate_Angle_test_data() )
-def test_computeAngleFromPoints( test_case: TestCaseAngle ) -> None:
-    """Test computeAngleFromPoints method."""
-    obs: float = fcts.computeAngleFromPoints( test_case.pt1, test_case.pt2, test_case.pt3 )
-    assert np.isclose( obs, test_case.angleExp ), f"Expected angle {test_case.angleExp}, but got {obs}."
+# @pytest.mark.parametrize( "test_case", __generate_Angle_test_data() )
+# def test_computeAngleFromPoints( test_case: TestCaseAngle ) -> None:
+#     """Test computeAngleFromPoints method."""
+#     obs: float = fcts.computeAngleFromPoints( test_case.pt1, test_case.pt2, test_case.pt3 )
+#     assert np.isclose( obs, test_case.angleExp ), f"Expected angle {test_case.angleExp}, but got {obs}."
 
 
-@pytest.mark.parametrize( "test_case", __generate_Angle_test_data() )
-def test_computeAngleFromVectors( test_case: TestCaseAngle ) -> None:
-    """Test computeAngleFromVectors method."""
-    vec1: npt.NDArray[ np.float64 ] = test_case.pt1 - test_case.pt2
-    vec2: npt.NDArray[ np.float64 ] = test_case.pt3 - test_case.pt2
-    obs: float = fcts.computeAngleFromVectors( vec1, vec2 )
-    assert np.isclose( obs, test_case.angleExp ), f"Expected angle {test_case.angleExp}, but got {obs}."
+# @pytest.mark.parametrize( "test_case", __generate_Angle_test_data() )
+# def test_computeAngleFromVectors( test_case: TestCaseAngle ) -> None:
+#     """Test computeAngleFromVectors method."""
+#     vec1: npt.NDArray[ np.float64 ] = test_case.pt1 - test_case.pt2
+#     vec2: npt.NDArray[ np.float64 ] = test_case.pt3 - test_case.pt2
+#     obs: float = fcts.computeAngleFromVectors( vec1, vec2 )
+#     assert np.isclose( obs, test_case.angleExp ), f"Expected angle {test_case.angleExp}, but got {obs}."
 
 
 def test_computeNormalFromPoints() -> None:
     """Test computeNormalFromPoints method."""
-    pt1 = np.array( [ 1.0, 0.0, 0.0 ] )
-    pt2 = np.array( [ 0.0, 0.0, 0.0 ] )
-    pt3 = np.array( [ 0.0, 1.0, 0.0 ] )
-    # vec1 = pt1 - pt2 = (1, 0, 0)
-    # vec2 = pt3 - pt2 = (0, 1, 0)
+    pt1 = np.array( [[ 1.0, 0.0, 0.0 ]] )
+    pt2 = np.array( [[ 0.0, 0.0, 0.0 ]] )
+    pt3 = np.array( [[ 0.0, 1.0, 0.0 ]] )
     # The function normalizes these vectors before computing the cross product.
     # The cross product of (1,0,0) and (0,1,0) is (0,0,1).
     obtained: npt.NDArray[ np.float64 ] = fcts.computeNormalFromPoints( pt1, pt2, pt3 )
-    expected: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 1.0 ] )
+    expected: npt.NDArray[ np.float64 ] = np.array( [[ 0.0, 0.0, 1.0 ]] )
     assert np.allclose( obtained, expected ), f"Expected normal {expected}, but got {obtained}."
 
 
 def test_computeNormalFromVectors() -> None:
     """Test computeNormalFromVectors method."""
     # Use non-unit vectors to test normalization step
-    vec1 = np.array( [ 5.0, 0.0, 0.0 ] )
-    vec2 = np.array( [ 0.0, 2.0, 0.0 ] )
+    vec1 = np.array( [[ 5.0, 0.0, 0.0 ]] )
+    vec2 = np.array( [[ 0.0, 2.0, 0.0 ]] )
     # The function should normalize vec1 to (1,0,0) and vec2 to (0,1,0)
     # The resulting cross product should be (0,0,1)
     obtained: npt.NDArray[ np.float64 ] = fcts.computeNormalFromVectors( vec1, vec2 )
-    expected: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 0.0, 1.0 ] )
+    expected: npt.NDArray[ np.float64 ] = np.array( [[ 0.0, 0.0, 1.0 ]] )
     assert np.allclose( obtained, expected ), f"Expected normal {expected}, but got {obtained}."

From 00d5747c0807b410690f26133dd104a6dff7c66c Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Thu, 9 Apr 2026 14:27:54 +0200
Subject: [PATCH 06/22] quick fix

---
 .../pre_processing/MeshQualityEnhanced.py     |  6 ++--
 .../tests/test_MeshQualityEnhanced.py         |  1 +
 geos-utils/tests/test_algebraFunctions.py     | 34 +++++++++----------
 3 files changed, 22 insertions(+), 19 deletions(-)

diff --git a/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py b/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py
index 7a23acff5..bfeced258 100644
--- a/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py
+++ b/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py
@@ -801,7 +801,8 @@ def _computeSquishIndex( self: Self ) -> None:
                 faceNormal: npt.NDArray[ np.float64 ] = self._getNormalVector( points, face )
 
                 vec: npt.NDArray[ np.float64 ] = cellCenter - faceCenter
-                angle: float = vtkMath.AngleBetweenVectors( vec, faceNormal )  # type: ignore[arg-type]
+                # TODO vtk Batch ??
+                angle: float = vtkMath.AngleBetweenVectors( vec, faceNormal[0] )  # type: ignore[arg-type]
                 squishIndex[ f ] = np.sin( angle )
             newArray.InsertValue( c, np.nanmax( squishIndex ) )
 
@@ -864,4 +865,5 @@ def _getNormalVector( self: Self, points: vtkPoints, face: vtkCell ) -> npt.NDAr
         ptsCoords: npt.NDArray[ np.float64 ] = np.zeros( ( 3, 3 ), dtype=float )
         for i in range( 3 ):
             points.GetPoint( facePtsIds.GetId( i ), ptsCoords[ i ] )
-        return geom.computeNormalFromPoints( ptsCoords[ 0 ], ptsCoords[ 1 ], ptsCoords[ 2 ] )
+        # TODO vectorize !!! 
+        return geom.computeNormalFromPoints( np.array([ptsCoords[ 0 ] ]), np.array([ptsCoords[ 1 ]]), np.array([ptsCoords[ 2 ]]) )
diff --git a/geos-processing/tests/test_MeshQualityEnhanced.py b/geos-processing/tests/test_MeshQualityEnhanced.py
index e60c16eb1..7d048b1ec 100644
--- a/geos-processing/tests/test_MeshQualityEnhanced.py
+++ b/geos-processing/tests/test_MeshQualityEnhanced.py
@@ -81,6 +81,7 @@ class TestCase:
 ] )
 def test_MeshQualityEnhanced( dataSetTest: Any, case: int ) -> None:
     """Test of MeshQualityEnhanced filter."""
+    print(f"test [{meshName_all[ case ]}]")
     test_case: TestCase = TestCase( dataSetTest( meshName_all[ case ] ), cellTypes_all[ case ],
                                     qualityMetrics_all[ case ], cellTypeCounts_all[ case ], metricsSummary_all[ case ] )
     mesh: vtkUnstructuredGrid = test_case.mesh
diff --git a/geos-utils/tests/test_algebraFunctions.py b/geos-utils/tests/test_algebraFunctions.py
index 104a8136b..f90e957e7 100644
--- a/geos-utils/tests/test_algebraFunctions.py
+++ b/geos-utils/tests/test_algebraFunctions.py
@@ -18,30 +18,30 @@ class TestAttributeMatrixFromVector( TestCase ):
 
     def test_wrongInputVectorSize( self: Self ) -> None:
         """Test failure on incorrect input vector size."""
-        emptyVector: npt.NDArray[ np.float64 ] = np.array( [] )
+        emptyVector: npt.NDArray[ np.float64 ] = np.array( [[]] )
         with self.assertRaises( ValueError ):
             getAttributeMatrixFromVector( emptyVector )
 
     def test_vector3size( self: Self ) -> None:
         """Test for an input vector size of 3."""
         vector3 = np.arange( 1, 4 )
-        expectedMatrix: npt.NDArray[ np.float64 ] = np.array( [ [ 1, 0, 0 ], [ 0, 2, 0 ], [ 0, 0, 3 ] ] )
+        expectedMatrix: npt.NDArray[ np.float64 ] = np.array( [ [ [ 1, 0, 0 ], [ 0, 2, 0 ], [ 0, 0, 3 ] ] ] )
 
-        self.assertTrue( np.array_equal( expectedMatrix, getAttributeMatrixFromVector( vector3 ) ) )
+        self.assertTrue( np.array_equal( expectedMatrix, getAttributeMatrixFromVector( np.array( [ vector3 ] ) ) ) )
 
     def test_vector6( self: Self ) -> None:
         """Test for an input vector size of 6."""
         vector6 = np.arange( 1, 7 )
-        expectedMatrix = np.array( [ [ 1, 6, 5 ], [ 6, 2, 4 ], [ 5, 4, 3 ] ] )
+        expectedMatrix = np.array( [ [ [ 1, 6, 5 ], [ 6, 2, 4 ], [ 5, 4, 3 ] ] ] )
 
-        self.assertTrue( np.array_equal( expectedMatrix, getAttributeMatrixFromVector( vector6 ) ) )
+        self.assertTrue( np.array_equal( expectedMatrix, getAttributeMatrixFromVector( np.array( [ vector6 ] ) ) ) )
 
     def test_vector9( self: Self ) -> None:
         """Test for an input vector size of 9."""
         vector9 = np.arange( 1, 10 )
-        expectedMatrix = np.array( [ [ 1, 6, 5 ], [ 9, 2, 4 ], [ 8, 7, 3 ] ] )
+        expectedMatrix = np.array( [ [ [ 1, 6, 5 ], [ 9, 2, 4 ], [ 8, 7, 3 ] ] ] )
 
-        self.assertTrue( np.array_equal( expectedMatrix, getAttributeMatrixFromVector( vector9 ) ) )
+        self.assertTrue( np.array_equal( expectedMatrix, getAttributeMatrixFromVector( np.array( [ vector9 ] ) ) ) )
 
 
 class TestAttributeVectorFromMatrix( TestCase ):
@@ -49,25 +49,25 @@ class TestAttributeVectorFromMatrix( TestCase ):
 
     def setUp( self: Self ) -> None:
         """Set up parameters."""
-        self.rdMatrix = np.arange( 1, 10 ).reshape( 3, 3 )
-        self.expected: npt.NDArray[ np.float64 ] = np.array( [ 1, 5, 9, 6, 3, 2, 8, 7, 4 ] )
+        self.rdMatrix = np.arange( 1, 10 ).reshape(1, 3, 3 )
+        self.expected: npt.NDArray[ np.float64 ] = np.array( [[ 1, 5, 9, 6, 3, 2, 8, 7, 4 ]] )
 
     def test_wrongInputMatrixShape( self ) -> None:
         """Test failure on empty input matrix."""
-        emptyMatrix = np.array( [] )
+        emptyMatrix = np.array( [[]] )
         with self.assertRaises( ValueError ):
-            getAttributeVectorFromMatrix( emptyMatrix, 3 )
+            getAttributeVectorFromMatrix( emptyMatrix, (1,3) )
 
-    def test_wrongOutputSize( self: Self ) -> None:
-        """Test failure on incorrect output size requested."""
-        size = 4
+    def test_wrongOutputShape( self: Self ) -> None:
+        """Test failure on incorrect output shape requested."""
+        shape = (1, 4)
         with self.assertRaises( ValueError ):
-            getAttributeVectorFromMatrix( self.rdMatrix, size )
+            getAttributeVectorFromMatrix( self.rdMatrix, shape )
 
     def test_vecOutput( self: Self ) -> None:
         """Test correct output for requested size."""
         listSize = ( 3, 6, 9 )
         for size in listSize:
             with self.subTest( size ):
-                expectedVector = np.array( self.expected[ :size ] )
-                self.assertTrue( np.array_equal( expectedVector, getAttributeVectorFromMatrix( self.rdMatrix, size ) ) )
+                expectedVector = np.array( self.expected[:, :size ] )
+                self.assertTrue( np.array_equal( expectedVector, getAttributeVectorFromMatrix( self.rdMatrix, (1,size) ) ) )

From f8153c97e7169b8881a7a80b167ed12f332f4022 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Thu, 9 Apr 2026 15:07:52 +0200
Subject: [PATCH 07/22] refactor einsums

---
 .../pre_processing/MeshQualityEnhanced.py     |  8 +--
 .../tests/test_MeshQualityEnhanced.py         |  2 +-
 .../src/geos/utils/geometryFunctions.py       | 61 ++++++++++---------
 geos-utils/tests/test_algebraFunctions.py     | 17 +++---
 geos-utils/tests/testsGeometryFunctions.py    | 44 +++++--------
 5 files changed, 62 insertions(+), 70 deletions(-)

diff --git a/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py b/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py
index bfeced258..6690fbb30 100644
--- a/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py
+++ b/geos-processing/src/geos/processing/pre_processing/MeshQualityEnhanced.py
@@ -802,7 +802,7 @@ def _computeSquishIndex( self: Self ) -> None:
 
                 vec: npt.NDArray[ np.float64 ] = cellCenter - faceCenter
                 # TODO vtk Batch ??
-                angle: float = vtkMath.AngleBetweenVectors( vec, faceNormal[0] )  # type: ignore[arg-type]
+                angle: float = vtkMath.AngleBetweenVectors( vec, faceNormal[ 0 ] )  # type: ignore[arg-type]
                 squishIndex[ f ] = np.sin( angle )
             newArray.InsertValue( c, np.nanmax( squishIndex ) )
 
@@ -814,7 +814,6 @@ def _computeSquishIndex( self: Self ) -> None:
         cellArrays.Modified()
         self._outputMesh.Modified()
 
-
     def _getCellCenter( self: Self,
                         cell: vtkCell,
                         ptsIds: Optional[ vtkIdTypeArray ] = None,
@@ -865,5 +864,6 @@ def _getNormalVector( self: Self, points: vtkPoints, face: vtkCell ) -> npt.NDAr
         ptsCoords: npt.NDArray[ np.float64 ] = np.zeros( ( 3, 3 ), dtype=float )
         for i in range( 3 ):
             points.GetPoint( facePtsIds.GetId( i ), ptsCoords[ i ] )
-        # TODO vectorize !!! 
-        return geom.computeNormalFromPoints( np.array([ptsCoords[ 0 ] ]), np.array([ptsCoords[ 1 ]]), np.array([ptsCoords[ 2 ]]) )
+        # TODO vectorize !!!
+        return geom.computeNormalFromPoints( np.array( [ ptsCoords[ 0 ] ] ), np.array( [ ptsCoords[ 1 ] ] ),
+                                             np.array( [ ptsCoords[ 2 ] ] ) )
diff --git a/geos-processing/tests/test_MeshQualityEnhanced.py b/geos-processing/tests/test_MeshQualityEnhanced.py
index 7d048b1ec..6984c42d2 100644
--- a/geos-processing/tests/test_MeshQualityEnhanced.py
+++ b/geos-processing/tests/test_MeshQualityEnhanced.py
@@ -81,7 +81,7 @@ class TestCase:
 ] )
 def test_MeshQualityEnhanced( dataSetTest: Any, case: int ) -> None:
     """Test of MeshQualityEnhanced filter."""
-    print(f"test [{meshName_all[ case ]}]")
+    print( f"test [{meshName_all[ case ]}]" )
     test_case: TestCase = TestCase( dataSetTest( meshName_all[ case ] ), cellTypes_all[ case ],
                                     qualityMetrics_all[ case ], cellTypeCounts_all[ case ], metricsSummary_all[ case ] )
     mesh: vtkUnstructuredGrid = test_case.mesh
diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index 5984d25a4..9e4df7f58 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -5,28 +5,40 @@
 
 import numpy as np
 import numpy.typing as npt
-from geos.utils.PhysicalConstants import EPSILON
 
 __doc__ = """Functions to permform geometry calculations."""
 
-CANONICAL_BASIS_3D: npt.NDArray[ np.float64 ] = np.array( [ [ 1.0, 0.0, 0.0 ],[ 0.0, 1.0, 0.0 ],[ 0.0, 0.0, 1.0 ] ] ) 
+CANONICAL_BASIS_3D: npt.NDArray[ np.float64 ] = np.array( [ [ 1.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 1.0 ] ] )
 
 # for batch cross product
-EPS = np.zeros((3, 3, 3), dtype=int)
-EPS[0,1,2] = EPS[1,2,0] = EPS[2,0,1] =  1
-EPS[0,2,1] = EPS[2,1,0] = EPS[1,0,2] = -1
+EPS = np.zeros( ( 3, 3, 3 ), dtype=int )
+EPS[ 0, 1, 2 ] = EPS[ 1, 2, 0 ] = EPS[ 2, 0, 1 ] = 1
+EPS[ 0, 2, 1 ] = EPS[ 2, 1, 0 ] = EPS[ 1, 0, 2 ] = -1
 
-def _normBasis( basis: npt.NDArray[np.float64] ):
-    """Norm and orthonormalize basis vector wise"""
+
+# (n,x,x) opertors
+def _normalize( arr : npt.NDArray[ np.float64 ]) -> npt.NDArray[ np.float64 ]:
+    """N generatlization of normalization."""
+    return np.einsum( 'ni,n->ni', arr, 1 / np.linalg.norm( arr, axis=1 ) )
+
+def _transposeProd( basisTo: npt.NDArray[ np.float64 ], basisFrom : npt.NDArray[ np.float64 ])-> npt.NDArray[ np.float64 ]:
+    """N generalization of transpose product."""
+    return np.einsum( 'nlj,nkj->nlk', basisTo, basisFrom )
+
+def _cross( vec1 : npt.NDArray[ np.float64 ],  vec2 : npt.NDArray[ np.float64 ]) -> npt.NDArray[ np.float64 ]:
+    """N generatlization of cross product."""
+    return np.einsum( 'ijk,nj,nk->ni', EPS, vec1, vec2 )
+
+def _normBasis( basis: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
+    """Norm and orthonormalize basis vector wise."""
     Q, _ = np.linalg.qr( basis )
     det = np.linalg.det( Q )
     Q[ det < 0 ] *= -1
 
     return Q
-    
 
 def getChangeOfBasisMatrix(
-    basisFrom:  npt.NDArray[ np.float64 ],
+    basisFrom: npt.NDArray[ np.float64 ],
     basisTo: npt.NDArray[ np.float64 ],
 ) -> Any:
     """Get the change of basis matrix from basis basisFrom to basis basisTo.
@@ -48,23 +60,22 @@ def getChangeOfBasisMatrix(
     basisFrom = _normBasis( basisFrom )
     basisTo = _normBasis( basisTo )
 
-    assert ( basisFrom.shape[1] == basisFrom.shape[2]) , (
+    assert ( basisFrom.shape[ 1 ] == basisFrom.shape[ 2 ] ), (
         f"Origin space vectors must have the same size. shape: {basisFrom.shape}" )
 
-    if len(basisTo.shape) < len(basisFrom.shape):
-        basisTo = np.repeat(basisTo[None,:],basisFrom.shape[0], axis=0) 
-    assert ( basisTo.shape[1] == basisTo.shape[2]) , (
-        "Destination space vectors must have the same size." )
+    if len( basisTo.shape ) < len( basisFrom.shape ):
+        basisTo = np.repeat( basisTo[ None, : ], basisFrom.shape[ 0 ], axis=0 )
+    assert ( basisTo.shape[ 1 ] == basisTo.shape[ 2 ] ), ( "Destination space vectors must have the same size." )
 
     # build the matrices where columns are the vectors of the bases
     # B = np.transpose( np.array( basisFrom ) )
     # C = np.transpose( np.array( basisTo ) )
     # no need to compute the inverse of C as it is orthonormal checked - transpose is enough
     assert np.linalg.norm(
-        np.einsum( 'nlj,nkj->nlk', basisTo, basisTo ) -
+        _transposeProd(basisTo,basisTo)
         np.repeat( np.eye( 3 )[ None, : ], basisFrom.shape[ 0 ], axis=0 ) ) < 1e-6
     # get the change of basis matrix
-    return np.einsum( 'nlj,nkj->nlk', basisTo, basisFrom )
+    return _transposeProd(basisTo,basisFrom)
 
 
 # def computeCoordinatesInNewBasis( vec: npt.NDArray[ np.float64 ],
@@ -83,7 +94,6 @@ def getChangeOfBasisMatrix(
 #     must be equal to the number of columns of and change of basis matrix."""
 #     return np.einsum('nij,ni->ni' ,changeOfBasisMatrix, vec )
 
-
 # def computePlaneFrom3Points(
 #     pt1: npt.NDArray[ np.float64 ],
 #     pt2: npt.NDArray[ np.float64 ],
@@ -116,7 +126,6 @@ def getChangeOfBasisMatrix(
 #     d: float = -np.einsum('ni,ni->n', normal, pt1 )
 #     return a, b, c, d
 
-
 # def getCellSideAgainstPlane(
 #     cellPtsCoords: npt.NDArray[ np.float64 ],
 #     planePt: npt.NDArray[ np.float64 ],
@@ -140,7 +149,6 @@ def getChangeOfBasisMatrix(
 #     ptCenter: npt.NDArray[ np.float64 ] = np.mean( cellPtsCoords, axis=0 )
 #     return getPointSideAgainstPlane( ptCenter, planePt, planeNormal )
 
-
 # def getPointSideAgainstPlane(
 #     ptCoords: npt.NDArray[ np.float64 ],
 #     planePt: npt.NDArray[ np.float64 ],
@@ -168,7 +176,6 @@ def getChangeOfBasisMatrix(
 #     assert np.linalg.norm( dot ) > EPSILON, "The point is on the plane."
 #     return dot > 0
 
-
 # def computeAngleFromPoints( pt1: npt.NDArray[ np.float64 ], pt2: npt.NDArray[ np.float64 ],
 #                             pt3: npt.NDArray[ np.float64 ] ) -> float:
 #     """Compute angle from 3 points.
@@ -186,7 +193,6 @@ def getChangeOfBasisMatrix(
 #     vec2: npt.NDArray[ np.float64 ] = pt3 - pt2
 #     return computeAngleFromVectors( vec1, vec2 )
 
-
 # def computeAngleFromVectors(
 #     vec1: npt.NDArray[ np.float64 ],
 #     vec2: npt.NDArray[ np.float64 ],
@@ -222,7 +228,6 @@ def getChangeOfBasisMatrix(
 #         # For 3D, return angle in [0, π]
 #         return np.arccos( dot )
 
-
 # def computeCosineFromVectors(
 #     vec1: npt.NDArray[ np.float64 ],
 #     vec2: npt.NDArray[ np.float64 ],
@@ -256,12 +261,12 @@ def computeNormalFromPoints( pt1: npt.NDArray[ np.float64 ], pt2: npt.NDArray[ n
     Returns:
         npt.NDArray[np.float64]: Normal vector coordinates
     """
-    assert pt1.shape[0] == pt2.shape[0] and pt2.shape[0] == pt3.shape[0]
+    assert pt1.shape[ 0 ] == pt2.shape[ 0 ] and pt2.shape[ 0 ] == pt3.shape[ 0 ]
     # compute vectors
     vec1: npt.NDArray[ np.float64 ] = pt1 - pt2
     vec2: npt.NDArray[ np.float64 ] = pt3 - pt2
-    assert np.all( np.abs( np.linalg.norm( vec1, axis=1 ) ) > 0.), "First and second points must be different"
-    assert np.all( np.abs( np.linalg.norm( vec2, axis=1 ) ) > 0.), "First and second points must be different"
+    assert np.all( np.abs( np.linalg.norm( vec1, axis=1 ) ) > 0. ), "First and second points must be different"
+    assert np.all( np.abs( np.linalg.norm( vec2, axis=1 ) ) > 0. ), "First and second points must be different"
     return computeNormalFromVectors( vec1, vec2 )
 
 
@@ -279,6 +284,6 @@ def computeNormalFromVectors(
         npt.NDArray[np.float64]: Normal vector coordinates
     """
     # normalization
-    vec1_norm = np.einsum('ni,n->ni',vec1, 1./ np.linalg.norm( vec1, axis=1))
-    vec2_norm = np.einsum('ni,n->ni', vec2, 1./ np.linalg.norm( vec2, axis=1))
-    return np.einsum('ijk,nj,nk->ni', EPS , vec1_norm, vec2_norm )
+    vec1_norm = _normalize( vec1 )
+    vec2_norm = _normalize( vec2 ) 
+    return _cross( vec1_norm, vec2_norm) 
diff --git a/geos-utils/tests/test_algebraFunctions.py b/geos-utils/tests/test_algebraFunctions.py
index f90e957e7..a3bba38ef 100644
--- a/geos-utils/tests/test_algebraFunctions.py
+++ b/geos-utils/tests/test_algebraFunctions.py
@@ -18,7 +18,7 @@ class TestAttributeMatrixFromVector( TestCase ):
 
     def test_wrongInputVectorSize( self: Self ) -> None:
         """Test failure on incorrect input vector size."""
-        emptyVector: npt.NDArray[ np.float64 ] = np.array( [[]] )
+        emptyVector: npt.NDArray[ np.float64 ] = np.array( [ [] ] )
         with self.assertRaises( ValueError ):
             getAttributeMatrixFromVector( emptyVector )
 
@@ -49,18 +49,18 @@ class TestAttributeVectorFromMatrix( TestCase ):
 
     def setUp( self: Self ) -> None:
         """Set up parameters."""
-        self.rdMatrix = np.arange( 1, 10 ).reshape(1, 3, 3 )
-        self.expected: npt.NDArray[ np.float64 ] = np.array( [[ 1, 5, 9, 6, 3, 2, 8, 7, 4 ]] )
+        self.rdMatrix = np.arange( 1, 10 ).reshape( 1, 3, 3 )
+        self.expected: npt.NDArray[ np.float64 ] = np.array( [ [ 1, 5, 9, 6, 3, 2, 8, 7, 4 ] ] )
 
     def test_wrongInputMatrixShape( self ) -> None:
         """Test failure on empty input matrix."""
-        emptyMatrix = np.array( [[]] )
+        emptyMatrix = np.array( [ [] ] )
         with self.assertRaises( ValueError ):
-            getAttributeVectorFromMatrix( emptyMatrix, (1,3) )
+            getAttributeVectorFromMatrix( emptyMatrix, ( 1, 3 ) )
 
     def test_wrongOutputShape( self: Self ) -> None:
         """Test failure on incorrect output shape requested."""
-        shape = (1, 4)
+        shape = ( 1, 4 )
         with self.assertRaises( ValueError ):
             getAttributeVectorFromMatrix( self.rdMatrix, shape )
 
@@ -69,5 +69,6 @@ def test_vecOutput( self: Self ) -> None:
         listSize = ( 3, 6, 9 )
         for size in listSize:
             with self.subTest( size ):
-                expectedVector = np.array( self.expected[:, :size ] )
-                self.assertTrue( np.array_equal( expectedVector, getAttributeVectorFromMatrix( self.rdMatrix, (1,size) ) ) )
+                expectedVector = np.array( self.expected[ :, :size ] )
+                self.assertTrue(
+                    np.array_equal( expectedVector, getAttributeVectorFromMatrix( self.rdMatrix, ( 1, size ) ) ) )
diff --git a/geos-utils/tests/testsGeometryFunctions.py b/geos-utils/tests/testsGeometryFunctions.py
index d91fdbe1f..9b4ff0d26 100644
--- a/geos-utils/tests/testsGeometryFunctions.py
+++ b/geos-utils/tests/testsGeometryFunctions.py
@@ -2,21 +2,14 @@
 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies.
 # SPDX-FileContributor: Martin Lemay, Jacques Franc
 # ruff: noqa: E402 # disable Module level import not at top of file
-import pytest
 import numpy as np
 import numpy.typing as npt
-from typing import Iterator
-from dataclasses import dataclass
-from itertools import combinations
 import geos.utils.geometryFunctions as fcts
 
-basisCanon: npt.NDArray[ np.float64 ] = np.array( [[[ 1.0, 0.0, 0.0 ] ,[ 0.0, 1.0, 0.0 ], 
-                                            [ 0.0, 0.0, 1.0 ]]] )
+basisCanon: npt.NDArray[ np.float64 ] = np.array( [ [ [ 1.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 1.0 ] ] ] )
 # destination basis according to canonic coordinates
-basisTo0: npt.NDArray[ np.float64 ] = np.array( [[ [1.0, 1.0, 0.0 ] ,[ -1.0, 1.0, 0.0 ], 
-                                                  [ 0.0, 0.0, 1.0 ] ]] )
-basisTo1: npt.NDArray[ np.float64 ] = np.array( [[ [1.0, -1.0, 0.0 ] ,[ 1.0, 1.0, 0.0 ], 
-                                                 [ 0.0, 0.0, 2.0 ] ]])
+basisTo0: npt.NDArray[ np.float64 ] = np.array( [ [ [ 1.0, 1.0, 0.0 ], [ -1.0, 1.0, 0.0 ], [ 0.0, 0.0, 1.0 ] ] ] )
+basisTo1: npt.NDArray[ np.float64 ] = np.array( [ [ [ 1.0, -1.0, 0.0 ], [ 1.0, 1.0, 0.0 ], [ 0.0, 0.0, 2.0 ] ] ] )
 # basisTo2: npt.NDArray[ np.float64 ] = np.array( [[ [0.0, 2.0, 0.0 ] ,[ -2.0, 0.0, 0.0 ],
 #                                                  [0.0, 0.0, 2.0 ] ]])
 
@@ -25,13 +18,13 @@ def test_getChangeOfBasisMatrixToCanonic() -> None:
     """Test change of basis matrix using canonic basis."""
     obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisCanon )
     # matrix where the columns are the vectors
-    expected: npt.NDArray[ np.float64 ] = fcts._normBasis(basisTo0).transpose(0, 2, 1)
-    assert np.linalg.norm( obtained-expected ) < 10e-12 , f"Expected array is {np.round( expected, 2 ).tolist()}"
-    # 
+    expected: npt.NDArray[ np.float64 ] = fcts._normBasis( basisTo0 ).transpose( 0, 2, 1 )
+    assert np.linalg.norm( obtained - expected ) < 10e-12, f"Expected array is {np.round( expected, 2 ).tolist()}"
+    #
     obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 )
     # matrix where the columns are the vectors
-    expected: npt.NDArray[ np.float64 ] = np.array([ [0., -1., 0.],[1., 0., 0.],[0., 0., 1.] ]) 
-    assert np.linalg.norm( obtained-expected ) < 10e-12 , f"Expected array is {np.round( expected, 2 ).tolist()}"
+    expected: npt.NDArray[ np.float64 ] = np.array( [ [ 0., -1., 0. ], [ 1., 0., 0. ], [ 0., 0., 1. ] ] )
+    assert np.linalg.norm( obtained - expected ) < 10e-12, f"Expected array is {np.round( expected, 2 ).tolist()}"
 
 
 # def test_computeCoordinatesInNewBasis() -> None:
@@ -45,7 +38,6 @@ def test_getChangeOfBasisMatrixToCanonic() -> None:
 #     assert np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ),
 #                            equal_nan=True ), f"Expected array is {np.round( expected, 2 ).tolist()}"
 
-
 # def test_computePlaneFrom3Points() -> None:
 #     """Test calculation of plane coefficients from 3 points."""
 #     pt1: npt.NDArray[ np.float64 ] = np.array( [[ 1.0, 2.0, 1.0 ]] )
@@ -55,7 +47,6 @@ def test_getChangeOfBasisMatrixToCanonic() -> None:
 #     expected: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ]  = ( np.array([-1.0]), np.array([2.0]), np.array([2.0]), np.array([-5.0]) )
 #     assert obtained == expected, f"Expected tuple is {expected}"
 
-
 # def test_getPointSideAgainstPlaneAssertion() -> None:
 #     """Test get point side against a plane."""
 #     planePt: npt.NDArray[ np.float64 ] = np.array( [[ 0.0, 0.0, 0.0 ]] )
@@ -65,7 +56,6 @@ def test_getChangeOfBasisMatrixToCanonic() -> None:
 #     with pytest.raises( AssertionError ):
 #         fcts.getPointSideAgainstPlane( planePt, planePt, planeNormal )
 
-
 # # yapf: disable
 # listPtsCoords_all = ( [ np.array( [ 0.5, 0.5, 0.5 ] ), np.array( [ 0.5, 0.5, -0.5 ] ) ],
 #                       [ np.array( [ 0.5, 0.5, 0.5 ] ), np.array( [ -0.5, 0.5, 0.5 ] ) ] )
@@ -159,7 +149,6 @@ def test_getChangeOfBasisMatrixToCanonic() -> None:
 # )
 # # yapf: enable
 
-
 # @dataclass( frozen=True )
 # class TestCaseAngle:
 #     """Test case."""
@@ -169,7 +158,6 @@ def test_getChangeOfBasisMatrixToCanonic() -> None:
 #     pt3: npt.NDArray[ np.float64 ]
 #     angleExp: float
 
-
 # def __generate_Angle_test_data() -> Iterator[ TestCaseAngle ]:
 #     """Generate test cases.
 
@@ -179,14 +167,12 @@ def test_getChangeOfBasisMatrixToCanonic() -> None:
 #     for pts, angle in zip( list( combinations( pts_all, 3 ) ), angleExp_all, strict=True ):
 #         yield TestCaseAngle( pts[ 0 ], pts[ 1 ], pts[ 2 ], angle )
 
-
 # @pytest.mark.parametrize( "test_case", __generate_Angle_test_data() )
 # def test_computeAngleFromPoints( test_case: TestCaseAngle ) -> None:
 #     """Test computeAngleFromPoints method."""
 #     obs: float = fcts.computeAngleFromPoints( test_case.pt1, test_case.pt2, test_case.pt3 )
 #     assert np.isclose( obs, test_case.angleExp ), f"Expected angle {test_case.angleExp}, but got {obs}."
 
-
 # @pytest.mark.parametrize( "test_case", __generate_Angle_test_data() )
 # def test_computeAngleFromVectors( test_case: TestCaseAngle ) -> None:
 #     """Test computeAngleFromVectors method."""
@@ -198,23 +184,23 @@ def test_getChangeOfBasisMatrixToCanonic() -> None:
 
 def test_computeNormalFromPoints() -> None:
     """Test computeNormalFromPoints method."""
-    pt1 = np.array( [[ 1.0, 0.0, 0.0 ]] )
-    pt2 = np.array( [[ 0.0, 0.0, 0.0 ]] )
-    pt3 = np.array( [[ 0.0, 1.0, 0.0 ]] )
+    pt1 = np.array( [ [ 1.0, 0.0, 0.0 ] ] )
+    pt2 = np.array( [ [ 0.0, 0.0, 0.0 ] ] )
+    pt3 = np.array( [ [ 0.0, 1.0, 0.0 ] ] )
     # The function normalizes these vectors before computing the cross product.
     # The cross product of (1,0,0) and (0,1,0) is (0,0,1).
     obtained: npt.NDArray[ np.float64 ] = fcts.computeNormalFromPoints( pt1, pt2, pt3 )
-    expected: npt.NDArray[ np.float64 ] = np.array( [[ 0.0, 0.0, 1.0 ]] )
+    expected: npt.NDArray[ np.float64 ] = np.array( [ [ 0.0, 0.0, 1.0 ] ] )
     assert np.allclose( obtained, expected ), f"Expected normal {expected}, but got {obtained}."
 
 
 def test_computeNormalFromVectors() -> None:
     """Test computeNormalFromVectors method."""
     # Use non-unit vectors to test normalization step
-    vec1 = np.array( [[ 5.0, 0.0, 0.0 ]] )
-    vec2 = np.array( [[ 0.0, 2.0, 0.0 ]] )
+    vec1 = np.array( [ [ 5.0, 0.0, 0.0 ] ] )
+    vec2 = np.array( [ [ 0.0, 2.0, 0.0 ] ] )
     # The function should normalize vec1 to (1,0,0) and vec2 to (0,1,0)
     # The resulting cross product should be (0,0,1)
     obtained: npt.NDArray[ np.float64 ] = fcts.computeNormalFromVectors( vec1, vec2 )
-    expected: npt.NDArray[ np.float64 ] = np.array( [[ 0.0, 0.0, 1.0 ]] )
+    expected: npt.NDArray[ np.float64 ] = np.array( [ [ 0.0, 0.0, 1.0 ] ] )
     assert np.allclose( obtained, expected ), f"Expected normal {expected}, but got {obtained}."

From c161a8e0e7ad883df8dfce08e24f1baef50b791f Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Thu, 9 Apr 2026 15:13:03 +0200
Subject: [PATCH 08/22] some clean up

---
 geos-mesh/src/geos/mesh/utils/genericHelpers.py | 11 +++++------
 geos-mesh/tests/test_QualityMetricSummary.py    |  2 +-
 geos-mesh/tests/test_meshQualityHelpers.py      |  2 +-
 3 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index 468f3ecde..0800ab494 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -22,7 +22,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 )
 
@@ -238,7 +238,7 @@ def getMonoBlockBounds( input: vtkUnstructuredGrid, ) -> tuple[ float, float, fl
     Args:
         input (vtkMultiBlockDataSet): input single block mesh
 
-    Returns:
+     Returns:
         tuple[float, float, float, float, float, float]: tuple containing
             bounds (xmin, xmax, ymin, ymax, zmin, zmax)
 
@@ -465,17 +465,16 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
     tangents1: npt.NDArray[ np.float64 ]
 
     try:
-        tangents1 = np.einsum( 'ni,n->ni', vtk_to_numpy( vtkTangents ),
-                               1 / np.linalg.norm( vtk_to_numpy( vtkTangents ), axis=1) )
+        tangents1 = _normalize( tangents1 )
     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 )
-        tangents2 = np.einsum( 'ni,n->ni', tangents2, 1 / np.linalg.norm( tangents2, axis=1 ) )
+        tangents2 : npt.NDArray[ np.float64 ] = _cross( normals, tangents1 )
+        tangents2 = _normalize( tangents2 )
 
     return ( tangents1, tangents2 )
 
diff --git a/geos-mesh/tests/test_QualityMetricSummary.py b/geos-mesh/tests/test_QualityMetricSummary.py
index 91335f8d5..4e2b513d2 100644
--- a/geos-mesh/tests/test_QualityMetricSummary.py
+++ b/geos-mesh/tests/test_QualityMetricSummary.py
@@ -225,4 +225,4 @@ def test_QualityMetricSummary_plotSummaryFigure( test_case: TestCase ) -> None:
     assert fig is not None, "Figure is undefined"
     # metrics + global counts + cell count + legend
     nbAxesExp: int = len( test_case.cellMetricIndexes ) + len( test_case.otherMetricIndexes ) + 3
-    assert len( fig.get_axes() ) == nbAxesExp, f"Number of Axes is expected to be {nbAxesExp}."
\ No newline at end of file
+    assert len( fig.get_axes() ) == nbAxesExp, f"Number of Axes is expected to be {nbAxesExp}."
diff --git a/geos-mesh/tests/test_meshQualityHelpers.py b/geos-mesh/tests/test_meshQualityHelpers.py
index 2b077d96f..67a0e6add 100644
--- a/geos-mesh/tests/test_meshQualityHelpers.py
+++ b/geos-mesh/tests/test_meshQualityHelpers.py
@@ -361,4 +361,4 @@ def test_getAllCellTypesExtended() -> None:
 def test_getAllCellTypes() -> None:
     """Test of getAllCellTypes method."""
     cellTypesExp: list[ int ] = [ VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON ]
-    assert cellTypesExp == getAllCellTypes(), "Cell types differ."
\ No newline at end of file
+    assert cellTypesExp == getAllCellTypes(), "Cell types differ."

From 334c30c97a1217239acc30eedf37e8d76c7be923 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Thu, 9 Apr 2026 16:35:22 +0200
Subject: [PATCH 09/22] few fixes

---
 .../src/geos/mesh/utils/genericHelpers.py     | 22 ++++++++-----------
 geos-mesh/tests/test_QualityMetricSummary.py  | 10 ++++-----
 .../src/geos/utils/geometryFunctions.py       | 20 ++++++++++-------
 geos-utils/tests/testsGeometryFunctions.py    | 14 ++++++------
 4 files changed, 33 insertions(+), 33 deletions(-)

diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index 0800ab494..6c63fd5dc 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -7,7 +7,7 @@
 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, vtkDataArray, vtkLogger
 from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkMultiBlockDataSet, vtkPolyData, vtkDataSet,
                                             vtkDataObject, vtkPlane, vtkCellTypes, vtkIncrementalOctreePointLocator,
                                             VTK_TRIANGLE, vtkSelection, vtkSelectionNode )
@@ -22,7 +22,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 , _normalize, _cross)
+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 )
 
@@ -238,7 +238,7 @@ def getMonoBlockBounds( input: vtkUnstructuredGrid, ) -> tuple[ float, float, fl
     Args:
         input (vtkMultiBlockDataSet): input single block mesh
 
-     Returns:
+    Returns:
         tuple[float, float, float, float, float, float]: tuple containing
             bounds (xmin, xmax, ymin, ymax, zmin, zmax)
 
@@ -460,21 +460,17 @@ 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 = _normalize( tangents1 )
-    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:
+        tangents1: npt.NDArray[ np.float64 ] = _normalize(
+            np.array( [ vtk_to_numpy( surface.GetCellData().GetTangents() ) ] ) )
         # 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 )
-        tangents2 : npt.NDArray[ np.float64 ] = _cross( normals, tangents1 )
+        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
 
     return ( tangents1, tangents2 )
 
diff --git a/geos-mesh/tests/test_QualityMetricSummary.py b/geos-mesh/tests/test_QualityMetricSummary.py
index 4e2b513d2..ccc8166a2 100644
--- a/geos-mesh/tests/test_QualityMetricSummary.py
+++ b/geos-mesh/tests/test_QualityMetricSummary.py
@@ -7,7 +7,7 @@
 import pandas as pd
 import random as rd
 import pytest
-from typing import Iterator
+from typing import Iterator, Any
 from matplotlib.figure import Figure
 
 from geos.mesh.model.QualityMetricSummary import QualityMetricSummary, StatTypes
@@ -24,7 +24,7 @@
 
 cellType_all: tuple[ int, ...] = ( VTK_TRIANGLE, VTK_QUAD, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_HEXAHEDRON,
                                    VTK_POLYGON, VTK_POLYHEDRON )
-cellMetricIndexes_all: tuple[ int, ...] = (
+cellMetricIndexes_all: tuple[ Any, ...] = (
     ( vtkMeshQuality.QualityMeasureTypes.EDGE_RATIO, vtkMeshQuality.QualityMeasureTypes.AREA,
       vtkMeshQuality.QualityMeasureTypes.CONDITION ),
     ( vtkMeshQuality.QualityMeasureTypes.MED_ASPECT_FROBENIUS, vtkMeshQuality.QualityMeasureTypes.AREA,
@@ -41,7 +41,7 @@
     ( vtkMeshQuality.QualityMeasureTypes.SHAPE, vtkMeshQuality.QualityMeasureTypes.VOLUME ),
 )
 
-metricIndexesFail_all: tuple[ int, ...] = (
+metricIndexesFail_all: tuple[ Any, ...] = (
     ( vtkMeshQuality.QualityMeasureTypes.MED_ASPECT_FROBENIUS, vtkMeshQuality.QualityMeasureTypes.VOLUME, 200 ),
     ( vtkMeshQuality.QualityMeasureTypes.NORMALIZED_INRADIUS, vtkMeshQuality.QualityMeasureTypes.VOLUME, 200 ),
     ( vtkMeshQuality.QualityMeasureTypes.ODDY, vtkMeshQuality.QualityMeasureTypes.AREA, 200 ),
@@ -58,7 +58,7 @@ class TestCase_StatsType:
     """Test case."""
     __test__ = False
     statType: StatTypes
-    values: tuple[ float, ...]
+    values: npt.NDArray[ Any ]
     expected: int | float
 
 
@@ -135,7 +135,7 @@ def __generate_failed_test_data() -> Iterator[ TestCase ]:
     """
     for metricIndexes, cellType in zip( metricIndexesFail_all, cellType_all, strict=True ):
         otherMetricIndexes: tuple = ()
-        values: tuple[ float, float, float, float, int ] = ( 0., 0., 0., 0., 0., 0 )
+        values: tuple[ float, float, float, float, int ] = ( 0., 0., 0., 0., 0 )
         yield TestCase( metricIndexes, otherMetricIndexes, cellType, statTypes, values )
 
 
diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index 9e4df7f58..1cc23fad4 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -17,18 +17,22 @@
 
 
 # (n,x,x) opertors
-def _normalize( arr : npt.NDArray[ np.float64 ]) -> npt.NDArray[ np.float64 ]:
+def _normalize( arr: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     """N generatlization of normalization."""
     return np.einsum( 'ni,n->ni', arr, 1 / np.linalg.norm( arr, axis=1 ) )
 
-def _transposeProd( basisTo: npt.NDArray[ np.float64 ], basisFrom : npt.NDArray[ np.float64 ])-> npt.NDArray[ np.float64 ]:
+
+def _transposeProd( basisTo: npt.NDArray[ np.float64 ],
+                    basisFrom: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     """N generalization of transpose product."""
     return np.einsum( 'nlj,nkj->nlk', basisTo, basisFrom )
 
-def _cross( vec1 : npt.NDArray[ np.float64 ],  vec2 : npt.NDArray[ np.float64 ]) -> npt.NDArray[ np.float64 ]:
+
+def _cross( vec1: npt.NDArray[ np.float64 ], vec2: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     """N generatlization of cross product."""
     return np.einsum( 'ijk,nj,nk->ni', EPS, vec1, vec2 )
 
+
 def _normBasis( basis: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     """Norm and orthonormalize basis vector wise."""
     Q, _ = np.linalg.qr( basis )
@@ -37,6 +41,7 @@ def _normBasis( basis: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
 
     return Q
 
+
 def getChangeOfBasisMatrix(
     basisFrom: npt.NDArray[ np.float64 ],
     basisTo: npt.NDArray[ np.float64 ],
@@ -72,10 +77,9 @@ def getChangeOfBasisMatrix(
     # C = np.transpose( np.array( basisTo ) )
     # no need to compute the inverse of C as it is orthonormal checked - transpose is enough
     assert np.linalg.norm(
-        _transposeProd(basisTo,basisTo)
-        np.repeat( np.eye( 3 )[ None, : ], basisFrom.shape[ 0 ], axis=0 ) ) < 1e-6
+        _transposeProd( basisTo, basisTo ) - np.repeat( np.eye( 3 )[ None, : ], basisFrom.shape[ 0 ], axis=0 ) ) < 1e-6
     # get the change of basis matrix
-    return _transposeProd(basisTo,basisFrom)
+    return _transposeProd( basisTo, basisFrom )
 
 
 # def computeCoordinatesInNewBasis( vec: npt.NDArray[ np.float64 ],
@@ -285,5 +289,5 @@ def computeNormalFromVectors(
     """
     # normalization
     vec1_norm = _normalize( vec1 )
-    vec2_norm = _normalize( vec2 ) 
-    return _cross( vec1_norm, vec2_norm) 
+    vec2_norm = _normalize( vec2 )
+    return _cross( vec1_norm, vec2_norm )
diff --git a/geos-utils/tests/testsGeometryFunctions.py b/geos-utils/tests/testsGeometryFunctions.py
index 9b4ff0d26..cb8b9b2a3 100644
--- a/geos-utils/tests/testsGeometryFunctions.py
+++ b/geos-utils/tests/testsGeometryFunctions.py
@@ -16,14 +16,14 @@
 
 def test_getChangeOfBasisMatrixToCanonic() -> None:
     """Test change of basis matrix using canonic basis."""
-    obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisCanon )
+    obtained: npt.NDArray[ np.float64 ] = np.vstack(
+        [ fcts.getChangeOfBasisMatrix( basisTo0, basisCanon ),
+          fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 ) ] )
     # matrix where the columns are the vectors
-    expected: npt.NDArray[ np.float64 ] = fcts._normBasis( basisTo0 ).transpose( 0, 2, 1 )
-    assert np.linalg.norm( obtained - expected ) < 10e-12, f"Expected array is {np.round( expected, 2 ).tolist()}"
-    #
-    obtained: npt.NDArray[ np.float64 ] = fcts.getChangeOfBasisMatrix( basisTo0, basisTo1 )
-    # matrix where the columns are the vectors
-    expected: npt.NDArray[ np.float64 ] = np.array( [ [ 0., -1., 0. ], [ 1., 0., 0. ], [ 0., 0., 1. ] ] )
+    expected: npt.NDArray[ np.float64 ] = np.vstack( [
+        fcts._normBasis( basisTo0 ).transpose( 0, 2, 1 ),
+        np.array( [ [ 0., -1., 0. ], [ 1., 0., 0. ], [ 0., 0., 1. ] ] )
+    ] )
     assert np.linalg.norm( obtained - expected ) < 10e-12, f"Expected array is {np.round( expected, 2 ).tolist()}"
 
 

From e6cc3d8677815efc3b431707bfe972d01a8c0c59 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Thu, 9 Apr 2026 17:32:55 +0200
Subject: [PATCH 10/22] fix lost in dimensions

---
 geos-mesh/src/geos/mesh/utils/genericHelpers.py | 2 +-
 geos-utils/tests/testsGeometryFunctions.py      | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index 6c63fd5dc..4760fd87d 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -462,7 +462,7 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
     """
     try:
         tangents1: npt.NDArray[ np.float64 ] = _normalize(
-            np.array( [ vtk_to_numpy( surface.GetCellData().GetTangents() ) ] ) )
+           vtk_to_numpy( surface.GetCellData().GetTangents() ) )
         # 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 )
diff --git a/geos-utils/tests/testsGeometryFunctions.py b/geos-utils/tests/testsGeometryFunctions.py
index cb8b9b2a3..1a1990e6f 100644
--- a/geos-utils/tests/testsGeometryFunctions.py
+++ b/geos-utils/tests/testsGeometryFunctions.py
@@ -22,7 +22,7 @@ def test_getChangeOfBasisMatrixToCanonic() -> None:
     # matrix where the columns are the vectors
     expected: npt.NDArray[ np.float64 ] = np.vstack( [
         fcts._normBasis( basisTo0 ).transpose( 0, 2, 1 ),
-        np.array( [ [ 0., -1., 0. ], [ 1., 0., 0. ], [ 0., 0., 1. ] ] )
+        np.array( [ [ [ 0., -1., 0. ], [ 1., 0., 0. ], [ 0., 0., 1. ] ] ] )
     ] )
     assert np.linalg.norm( obtained - expected ) < 10e-12, f"Expected array is {np.round( expected, 2 ).tolist()}"
 

From d27b42d470db5e034086556bb7f2c5798dc26543 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Thu, 9 Apr 2026 17:53:33 +0200
Subject: [PATCH 11/22] more fixes

---
 .../tests/testsFunctionsGeomechanicsCalculator.py        | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
index 608f71a6b..0d13a1a0c 100644
--- a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
+++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
@@ -254,12 +254,9 @@ def test_shearCapacityUtilization( self: Self ) -> None:
 
     def test_computeStressPrincipalComponentsFromStressVector( self: Self ) -> None:
         """Test calculation of stress principal components from stress vector."""
-        obtained: list[ float ] = [
-            round( val, 3 )
-            for val in list( fcts.computeStressPrincipalComponentsFromStressVector( stressVector )[ 0 ].flatten() )
-        ]
-        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."""

From 5edb31bd23aba0b50fc79e265eeea0c3e3b04874 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Thu, 9 Apr 2026 18:44:11 +0200
Subject: [PATCH 12/22] last fixes

---
 geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py   | 4 +++-
 .../tests/testsFunctionsGeomechanicsCalculator.py             | 2 +-
 geos-mesh/src/geos/mesh/utils/genericHelpers.py               | 3 +--
 3 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
index b0916da13..9d7b582f9 100644
--- a/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
+++ b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
@@ -137,4 +137,6 @@ 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])
+
diff --git a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
index 0d13a1a0c..4fda7ad46 100644
--- a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
+++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
@@ -255,7 +255,7 @@ def test_shearCapacityUtilization( self: Self ) -> None:
     def test_computeStressPrincipalComponentsFromStressVector( self: Self ) -> None:
         """Test calculation of stress principal components from stress vector."""
         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)
+        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:
diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index 4760fd87d..d58b24752 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -461,8 +461,7 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
         Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: The tangents vectors of the input surface.
     """
     try:
-        tangents1: npt.NDArray[ np.float64 ] = _normalize(
-           vtk_to_numpy( surface.GetCellData().GetTangents() ) )
+        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 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 )

From d7abbf95049348859c0d2c63fb9d26441199d726 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Sat, 11 Apr 2026 00:19:25 +0200
Subject: [PATCH 13/22] test fixes

---
 .../src/geos/geomechanics/model/MohrCircle.py |   7 +-
 .../src/geos/mesh/utils/genericHelpers.py     | 147 +++++-------
 .../test_genericHelpers_normalAndTangents.py  | 216 +++++++++---------
 3 files changed, 171 insertions(+), 199 deletions(-)

diff --git a/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
index 9d7b582f9..228722ab6 100644
--- a/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
+++ b/geos-geomechanics/src/geos/geomechanics/model/MohrCircle.py
@@ -137,6 +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
-        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])
-
+        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 ] )
diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index d58b24752..eb14b7bcb 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/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
+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
@@ -421,7 +420,8 @@ def convertAttributeFromLocalToXYZ(
     transformMatrix: npt.NDArray[ np.float64 ] = getChangeOfBasisMatrix( localBasisVectors, CANONICAL_BASIS_3D )
 
     # Apply transformation
-    arrayXYZ: npt.NDArray[ np.float64 ] = np.einsum( 'nij,njk,nlk->nil', transformMatrix, matrix3x3, transformMatrix )
+    # 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.shape )
@@ -464,7 +464,6 @@ def getTangentsVectors( surface: vtkPolyData ) -> Tuple[ npt.NDArray[ np.float64
         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 ] = np.cross( normals, tangents1, axis=1 ).astype( np.float64 )
         tangents2: npt.NDArray[ np.float64 ] = _cross( normals, tangents1 )
         tangents2 = _normalize( tangents2 )
     except AttributeError as err:
@@ -509,7 +508,7 @@ def getLocalBasisVectors(
         surfaceWithTangents: vtkPolyData = computeTangents( surfaceWithNormals, logger )
         tangents = getTangentsVectors( surfaceWithTangents )
 
-    return np.swapaxes( np.array( ( normals, *tangents ) ), 0, 1 )
+    return np.stack( [ normals, *tangents ], axis=2 )
 
 
 def computeNormals(
@@ -577,94 +576,70 @@ 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 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()
 
 
 def extractCellSelection( mesh: vtkUnstructuredGrid, ids: list[ int ] ) -> vtkUnstructuredGrid:
diff --git a/geos-mesh/tests/test_genericHelpers_normalAndTangents.py b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
index 0e41adda9..b6651b8a8 100644
--- a/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
+++ b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
@@ -12,68 +12,109 @@
 from vtkmodules.vtkCommonCore import vtkPoints
 from vtkmodules.util.numpy_support import vtk_to_numpy
 
-from geos.mesh.utils.genericHelpers import ( computeSurfaceTextureCoordinates, computeTangents, computeNormals,
-                                             getLocalBasisVectors, getTangentsVectors, getNormalVectors,
-                                             convertAttributeFromLocalToXYZ )
+from geos.mesh.utils.genericHelpers import ( computeTangents, computeNormals, getLocalBasisVectors, getTangentsVectors,
+                                             getNormalVectors, convertAttributeFromLocalToXYZ )
 
 from geos.utils.Errors import VTKError
 
 # yapf: disable
 pointsCoordsAll: list[ list[ list[ float ] ] ] = [
     [ [ 0., 0., 0. ], [ 0., 1., 0. ], [ 0., 2., 0. ], [ 0., 2., 1. ], [ 0., 1., 1. ], [ 0., 0., 1. ], ],
-    [ [ 0., 0., 0. ], [ 0., 1., 0. ], [ 0., 2., 0. ], [ 0., 2., 1. ], [ 1., 1., 1.5 ], [ 0., 0., 1. ], ],
+    [ [ 0., 0., 0. ], [ 0., 1., 0. ], [ 0., 2., 0. ], [ 0., 2., 1. ], [ 1., 1., 1. ], [ 0., 0., 1. ], ],
 ]
 trianglesAll: list[ list[ tuple[ int, int, int ] ] ] = [
-    [ ( 0, 1, 5 ), ( 1, 4, 5 ), ( 1, 2, 3 ), ( 1, 4, 3 ) ],
-    [ ( 0, 1, 5 ), ( 1, 4, 5 ), ( 1, 2, 3 ), ( 1, 4, 3 ) ],
+    [ ( 0, 1, 5 ), ( 1, 4, 5 ), ( 1, 2, 3 ), ( 3, 4, 1 ) ],
+    [ ( 0, 1, 5 ), ( 1, 4, 5 ), ( 1, 2, 3 ), ( 3, 4, 1 ) ],
 ]
-expectedNormalsAll: list[ npt.NDArray[ np.float64 ] ] = [
-    np.array( [ [ 1., 0., 0. ], [ 1., 0., 0. ], [ 1., 0., 0. ], [ 1., 0., 0. ], ] ),
-    np.array( [ [ 1., 0., 0. ], [ 0.7276069, -0.48507124, -0.48507124 ], [ 1., 0., 0. ],
-                [ 0.7276069, 0.48507124, -0.48507124 ] ] ),
+
+B_015: npt.NDArray[ np.float64 ] = np.array([
+    [1.0, 0.0, 0.0],   # n
+    [0.0, 1.0, 0.0],   # t1
+    [0.0, 0.0, 1.0],   # t2
+])
+
+B_145: list[ npt.NDArray[ np.float64 ] ] = [
+np.array([
+    [1.0, 0.0, 0.0],   # n
+    [0.0, 0.0, 1.0],   # t1
+    [0.0, -1.0, 0.0],  # t2
+]),
+np.array([
+    [ 0.577350, -0.577350, -0.577350],  # n
+    [ 0.707107,  0.000000,  0.707107],  # t1
+    [-0.408248, -0.816497,  0.408248],  # t2
+])
 ]
-expectedTangentsAll: list[ npt.NDArray[ np.float64 ] ] = [
-    np.array( [ [ [ 0., 2., 0. ], [ -0., 2., 0. ], [ 0., 2., 0. ], [ 0., 2., 0. ] ],
-                [ [ 0., 0., 2. ], [ 0., -0., 2. ], [ 0., 0., 2. ], [ 0., 0., 2. ],
-                ] ] ),
-    np.array( [ [ [ 0., 2., 0. ], [ 0.8301887, 2., -0.754717 ], [ 0., 2., 0. ], [ -0.8301887, 2., 0.754717 ] ],
-                [ [ 0., 0., 2. ], [ 1.33623397, 0.14643663, 1.85791445 ], [ 0., 0., 2. ],
-                  [ 1.33623397, -0.14643663, 1.85791445 ] ] ] ),
+
+B_123: npt.NDArray[ np.float64 ]  = np.array([
+    [1.0, 0.0, 0.0],   # n
+    [0.0, 1.0, 0.0],   # t1
+    [0.0, 0.0, 1.0],   # t2
+])
+
+B_341: list[ npt.NDArray[ np.float64 ] ] = [
+np.array([
+    [1.0, 0.0, 0.0],   # n
+    [0.0, -1.0, 0.0],   # t1
+    [0.0, 0.0, -1.0],   # t2
+]),
+np.array([
+    [ 0.577350,  0.577350, -0.577350],  # n
+    [ 0.707107, -0.707107,  0.000000],  # t1
+    [-0.408248, -0.408248, -0.816497],  # t2
+])
 ]
-expectedTexturedCoordsAll: list[ npt.NDArray[ np.float64 ] ] = [
-    np.array( [ [ 0., 0. ], [ 0.5, 0. ], [ 1., 0. ], [ 1., 1. ], [ 0.5, 1. ], [ 0., 1. ],
-    ] ),
-    np.array( [ [ [ 0., 0. ], [ 0.5, 0. ], [ 1., 0. ], [ 1., 0.41509435 ], [ 0.5, 1. ], [ 0., 0.41509435 ] ],
-    ] ),
+
+expectedBasisAll: list[ npt.NDArray[ np.float64 ] ] = [
+   np.array([B_015.T, B_145[0].T, B_123.T, B_341[0].T]),
+   np.array([B_015.T, B_145[1].T, B_123.T, B_341[1].T]),
 ]
+
 # Same attribute for all cells, test for vector size 3, 6, 9
-attributes: list[ npt.NDArray[ np.float64 ] ] = [ np.arange( 1., 4., dtype=np.float64 ), np.arange( 1., 7., dtype=np.float64 ), np.arange( 1., 10., dtype=np.float64 )
+attributes: list[ npt.NDArray[ np.float64 ] ] = [ np.arange( 1., 4., dtype=np.float64 ).reshape(1,3),
+np.arange( 1., 7., dtype=np.float64 ).reshape(1,6), np.arange( 1., 10., dtype=np.float64 ).reshape(1,9)
 ]
-expectedAttributesXYZ: list[ list[ npt.NDArray[ np.float64 ] ] ] = [
-    [
-        np.array( [ [ 1., 8., 12. ], [ 1., 8., 12. ], [ 1., 8., 12. ], [ 1., 8., 12. ],
-        ] ),
-        np.array( [ [ 1., 8., 12., 16., 10., 12. ], [ 1., 8., 12., 16., 10., 12. ], [ 1., 8., 12., 16., 10., 12. ],
-                    [ 1., 8., 12., 16., 10., 12. ] ] ),
-        np.array( [ [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ], [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ],
-                    [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ], [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ],
-                    [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ], [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ] ] )
+
+# -2.309401, 2.828427,-0.816497
+expectedAttributesXYZ: list[ list[ list[ npt.NDArray[ np.float64 ] ] ] ] = [[
+        [
+        np.array(  [ 1., 2., 3. ]),
+        np.array( [ 1., 3., 2. ]),
+        np.array( [ 1., 2., 3. ]),
+        np.array( [ 1., 2., 3. ]), ],
+        [
+           np.array( [ 1., 2., 3., 4., 5., 6.]),
+           np.array( [ 1., 3., 2., -4., 6. , -5.]),
+           np.array( [ 1., 2., 3., 4., 5., 6.]),
+           np.array( [ 1., 2., 3., 4., -5., -6.]) ],
+        [
+
+          np.array(       [1.,2.,3.,4.,5.,6.,7.,8.,9.]),
+          np.array(  [1., 3., 2., -7., 6., -5., -4., 9., -8. ] ),
+          np.array(   [1.,2.,3.,4.,5.,6.,7.,8.,9.] ),
+          np.array(   [1.,2.,3., 4., -5., -6., 7., -8., -9.] )
+        ]
     ],  # case 1
     [
-        np.array( [ [ 1., 8., 12. ], [ 7.26440201, 8.29962517, 11.73002787 ], [ 1., 8., 12. ],
-                    [ 7.26440201, 8.29962517, 11.73002787 ] ] ),
-        np.array( [ [ 1., 8., 12., 16., 10., 12. ],
-                    [ 33.11015535, -1.70942051, -4.10667954, 3.96829788, 5.78481908, 18.33796323 ],
-                    [ 1., 8., 12., 16., 10., 12. ],
-                    [ 0.86370966, 16.88802688, 9.54231792, 17.62557587, 12.93534579, 16.64449814 ] ] ),
-        np.array( [ [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ],
-                    [ 41.16704654, -3.95432477, -9.91866643, 0.51321919, -0.39322437, 23.20275907, 13.51039649,
-                        12.82015999, 23.3879596
-                    ], [ 1., 8., 12., 16., 10., 12., 28., 16., 18. ],
-                    [ -1.35966323, 18.70673795, 9.9469796, 14.59669037, 15.22437721, 25.39830602, 32.57499969,
-                        14.01099303, 21.0552047
-                    ] ] )
-    ],  # case 2
+        [
+        np.array( [ 1., 2., 3. ] ),
+         np.array( [1.99999814e+00, 2.00000124e+00, 2.00000042e+00, 5.77350037e-01, 1.15470000e-06,-8.16496453e-01] ),
+           np.array( [ 1., 2. , 3. ] ),
+           np.array([ 1.99999814, 1.50000093,2.50000134, 0.28867502, 0.70710768,-0.40824823 ])
+           ],
+           [
+        np.array( [1.,2.,3.,4.,5.,6.] ),
+        np.array(     [-2.66666418,7.00000433,1.66666919,-5.1961574,1.88561586,-4.89897872]),
+        np.array(     [1.,2.,3.,4.,5.,6.]),
+        np.array(     [2.54962629e-17,-4.50000279e+00,1.04999973e+01,-2.88675726e-01,-4.24263915e+00,-8.16496453e-01])
+            ],
+            [
+        np.array([1.,2.,3.,4.,5.,6.,7.,8.,9.]),
+        np.array([-3.66666325,8.50000526,1.16666991,-7.79423469,4.71404139,-7.34846808,-6.06218458,0.47140223,-4.89897872]),
+        np.array([1.,2.,3.,4.,5.,6.,7.,8.,9.]),
+        np.array([-0.99999907,-6.00000371,12.9999962,0.57734933,-3.53553321,0.40824823,-1.15470078,-7.77817236,-2.04124113])
+    ],# case 2
+]
 ]
 # yapf: enable
 
@@ -83,13 +124,11 @@ class Options( IntEnum ):
 
     - RAW : bare mesh
     - NORMALS : only normals are calculated
-    - NORMALS_TEXTURED : normals and textured coordinates are computed
     - TANGENTS : normals and tangents are present in the mesh
     """
     RAW = 0
     NORMALS = 1
-    NORMALS_TEXTURED = 2
-    TANGENTS = 3
+    TANGENTS = 2
 
 
 @dataclass
@@ -98,10 +137,10 @@ class TriangulatedSurfaceTestCase:
     pointsCoords: list[ list[ float ] ]
     triangles: list[ tuple[ int, int, int ] ]
     attribute: list[ npt.NDArray ]
-    expectedNormals: npt.NDArray
-    expectedTangents: npt.NDArray
-    expectedTexturedCoords: npt.NDArray
-    expectedAttributeInXYZ: list[ npt.NDArray ]
+    # expectedNormals: npt.NDArray
+    # expectedTangents: npt.NDArray
+    expectedBasis: npt.NDArray
+    expectedAttributeInXYZ: list[ list[ npt.NDArray ] ]
     options: Options
     mesh: vtkPolyData = field( init=False )
 
@@ -127,16 +166,9 @@ def __post_init__( self: Self ) -> None:
         if self.options == Options.NORMALS:
             self.mesh = computeNormals( polydata )
 
-        elif self.options == Options.NORMALS_TEXTURED:
-            mesh: vtkPolyData = computeSurfaceTextureCoordinates( polydata )
-            mesh2: vtkPolyData = computeNormals( mesh )
-            self.mesh = mesh2
-
         elif self.options == Options.TANGENTS:
-            mesh = computeSurfaceTextureCoordinates( polydata )
-            mesh2 = computeNormals( mesh )
-            mesh3: vtkPolyData = computeTangents( mesh2 )
-            self.mesh = mesh3
+            mesh2 = computeNormals( polydata )
+            self.mesh = computeTangents( mesh2 )
 
         else:
             # Unknown cases and case 0
@@ -155,25 +187,7 @@ def __generateSurfacicTestCase( options: tuple[ Options, ...] ) -> Iterator[ Tri
     for opt in options:
         for i in range( len( pointsCoordsAll ) ):
             yield TriangulatedSurfaceTestCase( pointsCoordsAll[ i ], trianglesAll[ i ], attributes,
-                                               expectedNormalsAll[ i ], expectedTangentsAll[ i ],
-                                               expectedTexturedCoordsAll[ i ], expectedAttributesXYZ[ i ], opt )
-
-
-@pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.RAW, ) ) )
-def test_computeTextureCoords( case: TriangulatedSurfaceTestCase ) -> None:
-    """Test the computation of texture coordinates."""
-    stc: vtkPolyData = computeSurfaceTextureCoordinates( case.mesh )
-    texturedMap: npt.NDArray[ np.float64 ] = vtk_to_numpy( stc.GetPointData().GetArray( "Texture Coordinates" ) )
-    assert np.allclose( texturedMap, case.expectedTexturedCoords )
-
-
-@pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.NORMALS_TEXTURED, ) ) )
-def test_computeTangents( case: TriangulatedSurfaceTestCase ) -> None:
-    """Test the computation of tangents."""
-    surfaceWithTangents: vtkPolyData = computeTangents( case.mesh )
-    tangents: npt.NDArray[ np.float64 ] = vtk_to_numpy( surfaceWithTangents.GetCellData().GetTangents() )
-
-    assert np.allclose( tangents, case.expectedTangents[ 0 ] )
+                                               expectedBasisAll[ i ], expectedAttributesXYZ[ i ], opt )
 
 
 @pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.TANGENTS, ) ) )
@@ -184,8 +198,8 @@ def test_getTangents( case: TriangulatedSurfaceTestCase ) -> None:
 
     tangents1, tangents2 = getTangentsVectors( case.mesh )
 
-    assert np.allclose( tangents1, case.expectedTangents[ 0 ], rtol=1e-3 )
-    assert np.allclose( tangents2, case.expectedTangents[ 1 ], rtol=1e-3 )
+    assert np.allclose( tangents1, case.expectedBasis[ :, :, 1 ], rtol=1e-3 )
+    assert np.allclose( tangents2, case.expectedBasis[ :, :, 2 ], rtol=1e-3 )
 
 
 @pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.RAW, ) ) )
@@ -195,7 +209,7 @@ def test_computeNormals( case: TriangulatedSurfaceTestCase ) -> None:
 
     normals = vtk_to_numpy( surfaceWithNormals.GetCellData().GetNormals() )
 
-    assert np.allclose( normals, case.expectedNormals, rtol=1e-3 )
+    assert np.allclose( normals, case.expectedBasis[ :, :, 0 ], rtol=1e-3 )
 
 
 @pytest.mark.parametrize( "case", __generateSurfacicTestCase( options=( Options.NORMALS, ) ) )
@@ -203,7 +217,7 @@ def test_getNormals( case: TriangulatedSurfaceTestCase ) -> None:
     """Test normals getter."""
     normals = getNormalVectors( case.mesh )
 
-    assert np.allclose( normals, case.expectedNormals, rtol=1e-3 )
+    assert np.allclose( normals, case.expectedBasis[ :, :, 0 ], rtol=1e-3 )
 
 
 @pytest.mark.parametrize( "case",
@@ -214,11 +228,9 @@ def test_getNormals( case: TriangulatedSurfaceTestCase ) -> None:
                           ) ) )
 def test_getLocalBasis( case: TriangulatedSurfaceTestCase ) -> None:
     """Test local basis getter."""
-    normals, tangents1, tangents2 = getLocalBasisVectors( case.mesh )
+    basis = getLocalBasisVectors( case.mesh )
 
-    assert np.allclose( normals, case.expectedNormals, rtol=1e-3 )
-    assert np.allclose( tangents1, case.expectedTangents[ 0 ], rtol=1e-3 )
-    assert np.allclose( tangents2, case.expectedTangents[ 1 ], rtol=1e-3 )
+    assert np.linalg.norm( basis - case.expectedBasis ) < 1e-3
 
 
 #########################################################################################
@@ -234,18 +246,6 @@ def emptySurface() -> vtkPolyData:
     return vtkPolyData()
 
 
-def test_failingComputeSurfaceTextureCoords( emptySurface: vtkPolyData ) -> None:
-    """Test VTK error raising of texture coordinate calculation."""
-    with pytest.raises( VTKError ):
-        computeSurfaceTextureCoordinates( emptySurface )
-
-
-def test_failingComputeTangents( emptySurface: vtkPolyData ) -> None:
-    """Test VTK error raising during tangent calculation."""
-    with pytest.raises( VTKError ):
-        computeTangents( emptySurface )
-
-
 def test_failingComputeNormals( emptySurface: vtkPolyData ) -> None:
     """Test VTK error raising of normals calculation."""
     with pytest.raises( VTKError ):
@@ -271,9 +271,7 @@ def test_failingGetNormals( emptySurface: vtkPolyData ) -> None:
 class AttributeConversionTestCase:
     """Test case for attribute conversion from local basis to XYZ basis for one cell."""
     vector: npt.NDArray[ np.float64 ]
-    normal: npt.NDArray[ np.float64 ]
-    tangent1: npt.NDArray[ np.float64 ]
-    tangent2: npt.NDArray[ np.float64 ]
+    basis: npt.NDArray[ np.float64 ]
     expectedVectorXYZ: npt.NDArray[ np.float64 ]
 
 
@@ -286,10 +284,9 @@ def __generateAttributeConversionTestCase() -> Iterator[ AttributeConversionTest
             for nvec, localAttribute in enumerate( attributes ):
                 for ncell in range( testcase.mesh.GetNumberOfCells() ):
 
-                    yield AttributeConversionTestCase( localAttribute, testcase.expectedNormals[ ncell ],
-                                                       testcase.expectedTangents[ 0 ][ ncell ],
-                                                       testcase.expectedTangents[ 1 ][ ncell ],
-                                                       testcase.expectedAttributeInXYZ[ nvec ][ ncell ] )
+                    yield AttributeConversionTestCase(
+                        localAttribute, testcase.expectedBasis[ ncell, :, : ].reshape( 1, 3, 3 ),
+                        testcase.expectedAttributeInXYZ[ nvec ][ ncell ].reshape( 1, -1 ) )
 
 
 @pytest.mark.parametrize( "testcase", __generateAttributeConversionTestCase() )
@@ -297,6 +294,5 @@ def test_convertAttributesToXYZ( testcase: AttributeConversionTestCase ) -> None
     """Test the conversion of one cell attribute from local to canonic basis."""
     localAttr: npt.NDArray[ np.float64 ] = testcase.vector
 
-    attrXYZ: npt.NDArray[ np.float64 ] = convertAttributeFromLocalToXYZ(
-        localAttr, ( testcase.normal, testcase.tangent1, testcase.tangent2 ) )
+    attrXYZ: npt.NDArray[ np.float64 ] = convertAttributeFromLocalToXYZ( localAttr, testcase.basis )
     assert np.allclose( attrXYZ, testcase.expectedVectorXYZ, rtol=1e-6 )

From bf55ddc508fac280b1e7d40753fbfc1df0bfb320 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Mon, 13 Apr 2026 08:54:18 +0200
Subject: [PATCH 14/22] ref val changes

---
 .../test_genericHelpers_normalAndTangents.py  | 24 +++++++++----------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/geos-mesh/tests/test_genericHelpers_normalAndTangents.py b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
index b6651b8a8..abbc2235c 100644
--- a/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
+++ b/geos-mesh/tests/test_genericHelpers_normalAndTangents.py
@@ -78,21 +78,21 @@
 # -2.309401, 2.828427,-0.816497
 expectedAttributesXYZ: list[ list[ list[ npt.NDArray[ np.float64 ] ] ] ] = [[
         [
-        np.array(  [ 1., 2., 3. ]),
-        np.array( [ 1., 3., 2. ]),
-        np.array( [ 1., 2., 3. ]),
-        np.array( [ 1., 2., 3. ]), ],
+        np.array( [ 1., 2., 3. ] ),
+        np.array( [ 1., 3., 2. ] ),
+        np.array( [ 1., 2., 3. ] ),
+        np.array( [ 1., 2., 3. ] ), ],
         [
-           np.array( [ 1., 2., 3., 4., 5., 6.]),
-           np.array( [ 1., 3., 2., -4., 6. , -5.]),
-           np.array( [ 1., 2., 3., 4., 5., 6.]),
-           np.array( [ 1., 2., 3., 4., -5., -6.]) ],
+           np.array( [ 1., 2., 3., 4., 5., 6.] ),
+           np.array( [ 1., 3., 2., -4., 6. , -5.] ),
+           np.array( [ 1., 2., 3., 4., 5., 6.] ),
+           np.array( [ 1., 2., 3., 4., -5., -6.] ) ],
         [
 
-          np.array(       [1.,2.,3.,4.,5.,6.,7.,8.,9.]),
-          np.array(  [1., 3., 2., -7., 6., -5., -4., 9., -8. ] ),
-          np.array(   [1.,2.,3.,4.,5.,6.,7.,8.,9.] ),
-          np.array(   [1.,2.,3., 4., -5., -6., 7., -8., -9.] )
+          np.array( [1.,2.,3.,4.,5.,6.,7.,8.,9.] ),
+          np.array( [1., 3., 2., -7., 6., -5., -4., 9., -8. ] ),
+          np.array( [1.,2.,3.,4.,5.,6.,7.,8.,9.] ),
+          np.array( [1.,2.,3., 4., -5., -6., 7., -8., -9.] )
         ]
     ],  # case 1
     [

From 84a5b7fc648259a0c92e1c09f0d8dcac17f57ff5 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Mon, 13 Apr 2026 09:16:37 +0200
Subject: [PATCH 15/22] using polar for sign consistency

---
 geos-utils/src/geos/utils/geometryFunctions.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index 1cc23fad4..56ee78fe6 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -4,6 +4,7 @@
 from typing import Any
 
 import numpy as np
+from scipy.linalg import polar
 import numpy.typing as npt
 
 __doc__ = """Functions to permform geometry calculations."""
@@ -35,11 +36,12 @@ def _cross( vec1: npt.NDArray[ np.float64 ], vec2: npt.NDArray[ np.float64 ] ) -
 
 def _normBasis( basis: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     """Norm and orthonormalize basis vector wise."""
-    Q, _ = np.linalg.qr( basis )
-    det = np.linalg.det( Q )
-    Q[ det < 0 ] *= -1
+    # Q, R = np.linalg.qr( basis )
+    U, P = polar( basis )
+    # det = np.linalg.det( Q )
+    # Q[ det < 0 ] *= -1
 
-    return Q
+    return U
 
 
 def getChangeOfBasisMatrix(

From 98ed643e5a93c6fc245c65be766a89f6e3a55d47 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Mon, 13 Apr 2026 10:24:54 +0200
Subject: [PATCH 16/22] add scipy

---
 geos-mesh/pyproject.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/geos-mesh/pyproject.toml b/geos-mesh/pyproject.toml
index e0138b0c8..eddfcfecb 100644
--- a/geos-mesh/pyproject.toml
+++ b/geos-mesh/pyproject.toml
@@ -31,6 +31,7 @@ dependencies = [
     "vtk >= 9.3, < 9.6",
     "networkx >= 2.4",
     "tqdm >= 4.67",
+    "scipy >= 1.16.3 ",
     "numpy >= 2.2",
     "pandas >= 2.2",
     "meshio >= 5.3",

From 3f8572aad2179c371d89017cb4c80533b18ae059 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Mon, 13 Apr 2026 10:50:43 +0200
Subject: [PATCH 17/22] relocate scipy

---
 geos-mesh/pyproject.toml  | 1 -
 geos-utils/pyproject.toml | 1 +
 2 files changed, 1 insertion(+), 1 deletion(-)

diff --git a/geos-mesh/pyproject.toml b/geos-mesh/pyproject.toml
index eddfcfecb..e0138b0c8 100644
--- a/geos-mesh/pyproject.toml
+++ b/geos-mesh/pyproject.toml
@@ -31,7 +31,6 @@ dependencies = [
     "vtk >= 9.3, < 9.6",
     "networkx >= 2.4",
     "tqdm >= 4.67",
-    "scipy >= 1.16.3 ",
     "numpy >= 2.2",
     "pandas >= 2.2",
     "meshio >= 5.3",
diff --git a/geos-utils/pyproject.toml b/geos-utils/pyproject.toml
index bc11980d8..ff8d46946 100644
--- a/geos-utils/pyproject.toml
+++ b/geos-utils/pyproject.toml
@@ -28,6 +28,7 @@ requires-python = ">=3.10"
 
 dependencies = [
     "numpy >= 2.2",
+    "scipy >= 1.16.3 ",
     "typing_extensions >= 4.12",
 ]
 

From e4f10aeaee21a73c5e5fe84fe2f55261c000e016 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Mon, 13 Apr 2026 12:52:39 +0200
Subject: [PATCH 18/22] lower version for 3.10

---
 geos-utils/pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/geos-utils/pyproject.toml b/geos-utils/pyproject.toml
index ff8d46946..4ce1f94f4 100644
--- a/geos-utils/pyproject.toml
+++ b/geos-utils/pyproject.toml
@@ -28,7 +28,7 @@ requires-python = ">=3.10"
 
 dependencies = [
     "numpy >= 2.2",
-    "scipy >= 1.16.3 ",
+    "scipy >= 1.15.3 ",
     "typing_extensions >= 4.12",
 ]
 

From 4cdfe5488eddd3e1bfae2c58236a123c398c4c54 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Mon, 13 Apr 2026 13:53:39 +0200
Subject: [PATCH 19/22] version disjunction

---
 geos-utils/src/geos/utils/geometryFunctions.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index 56ee78fe6..2261fa979 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -5,6 +5,8 @@
 
 import numpy as np
 from scipy.linalg import polar
+from packaging.version import Version
+from scipy import __version__ as __scipy_version__
 import numpy.typing as npt
 
 __doc__ = """Functions to permform geometry calculations."""
@@ -37,7 +39,13 @@ def _cross( vec1: npt.NDArray[ np.float64 ], vec2: npt.NDArray[ np.float64 ] ) -
 def _normBasis( basis: npt.NDArray[ np.float64 ] ) -> npt.NDArray[ np.float64 ]:
     """Norm and orthonormalize basis vector wise."""
     # Q, R = np.linalg.qr( basis )
-    U, P = polar( basis )
+    if Version( __scipy_version__ ) > Version( "1.15.3" ):
+        U, _ = polar( basis )
+    else:
+        U = np.zeros_like( basis )
+        for i in range( basis.shape[ 0 ] ):
+            U[ i, :, : ], _ = polar( basis[ i, :, : ] )
+
     # det = np.linalg.det( Q )
     # Q[ det < 0 ] *= -1
 

From 00757d206e3e957c9bd3b0ec736630744dfa241f Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Mon, 13 Apr 2026 14:08:07 +0200
Subject: [PATCH 20/22] reordered ops

---
 geos-utils/src/geos/utils/geometryFunctions.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/geos-utils/src/geos/utils/geometryFunctions.py b/geos-utils/src/geos/utils/geometryFunctions.py
index 2261fa979..ed2a0164d 100644
--- a/geos-utils/src/geos/utils/geometryFunctions.py
+++ b/geos-utils/src/geos/utils/geometryFunctions.py
@@ -73,13 +73,13 @@ def getChangeOfBasisMatrix(
         npt.NDArray[np.float64 Change of basis matrix.
     """
     basisFrom = _normBasis( basisFrom )
+    if len( basisTo.shape ) < len( basisFrom.shape ):
+        basisTo = np.repeat( basisTo[ None, : ], basisFrom.shape[ 0 ], axis=0 )
     basisTo = _normBasis( basisTo )
 
     assert ( basisFrom.shape[ 1 ] == basisFrom.shape[ 2 ] ), (
         f"Origin space vectors must have the same size. shape: {basisFrom.shape}" )
 
-    if len( basisTo.shape ) < len( basisFrom.shape ):
-        basisTo = np.repeat( basisTo[ None, : ], basisFrom.shape[ 0 ], axis=0 )
     assert ( basisTo.shape[ 1 ] == basisTo.shape[ 2 ] ), ( "Destination space vectors must have the same size." )
 
     # build the matrices where columns are the vectors of the bases

From ab0c4b293ef2d3bc7d20aa1c093250807c660aac Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Mon, 13 Apr 2026 18:23:42 +0200
Subject: [PATCH 21/22] clean up

---
 .../testsFunctionsGeomechanicsCalculator.py   |  2 -
 .../src/geos/mesh/utils/genericHelpers.py     | 44 -------------------
 2 files changed, 46 deletions(-)

diff --git a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
index 4fda7ad46..3f9875bb8 100644
--- a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
+++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py
@@ -11,8 +11,6 @@
 import numpy.typing as npt
 from typing_extensions import Self
 
-# sys.path.insert( 0, "/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-geomechanics/src" )
-# sys.path.insert( 0, "/data/pau901/SIM_CS/04_WORKSPACE/USERS/jfranc/geosPythonPackages/geos-utils/src" )
 dir_path = os.path.dirname( os.path.realpath( __file__ ) )
 parent_dir_path = os.path.join( os.path.dirname( dir_path ), "src" )
 if parent_dir_path not in sys.path:
diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index eb14b7bcb..e24fdf733 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -598,50 +598,6 @@ def computeTangents(
     triangulatedSurface.Modified()
     return triangulatedSurface
 
-
-# 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()
-
-
 def extractCellSelection( mesh: vtkUnstructuredGrid, ids: list[ int ] ) -> vtkUnstructuredGrid:
     """Extract cell selection from list of cell Ids.
 

From f68dc9f2f4d1679c39f5cbfc90cada9d2eca01f1 Mon Sep 17 00:00:00 2001
From: jacques franc <jacquesfrancdev@gmail.com>
Date: Tue, 14 Apr 2026 09:28:22 +0200
Subject: [PATCH 22/22] yapf

---
 geos-mesh/src/geos/mesh/utils/genericHelpers.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
index e24fdf733..2e1031ebd 100644
--- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py
+++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py
@@ -598,6 +598,7 @@ def computeTangents(
     triangulatedSurface.Modified()
     return triangulatedSurface
 
+
 def extractCellSelection( mesh: vtkUnstructuredGrid, ids: list[ int ] ) -> vtkUnstructuredGrid:
     """Extract cell selection from list of cell Ids.