refactor: Refactor geos-mesh test and test meshes

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

@@ -5,16 +5,17 @@
 import numpy as np
 import numpy.typing as npt
 import pandas as pd  # type: ignore[import-untyped]
+from typing import Union, Any
+
 import vtkmodules.util.numpy_support as vnp
-from typing import Optional, Union, Any
 from vtkmodules.util.numpy_support import vtk_to_numpy
 from vtkmodules.vtkCommonCore import vtkDataArray, vtkPoints
 from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkFieldData, vtkMultiBlockDataSet, vtkDataSet,
                                             vtkCompositeDataSet, vtkDataObject, vtkPointData, vtkCellData, vtkPolyData,
                                             vtkCell )
 from vtkmodules.vtkFiltersCore import vtkCellCenters
-from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten
 
+from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten
 from geos.utils.pieceEnum import Piece
 
 __doc__ = """
@@ -46,7 +47,7 @@
 """
 
 
-def getCellDimension( mesh: Union[ vtkMultiBlockDataSet, vtkDataSet ] ) -> set[ int ]:
+def getCellDimension( mesh: Union[ vtkMultiBlockDataSet, vtkDataSet ], ) -> set[ int ]:
     """Get the set of the different cells dimension of a mesh.
 
     Args:
@@ -183,7 +184,10 @@ def computeElementMapping(
     return elementMap
 
 
-def hasArray( mesh: vtkUnstructuredGrid, arrayNames: list[ str ] ) -> bool:
+def hasArray(
+    mesh: vtkUnstructuredGrid,
+    arrayNames: list[ str ],
+) -> bool:
     """Checks if input mesh contains at least one of input data arrays.
 
     Args:
@@ -282,9 +286,11 @@ def checkValidValuesInObject(
     return ( validValues, invalidValues )
 
 
-def getNumpyGlobalIdsArray( data: Union[ vtkCellData, vtkPointData ] ) -> npt.NDArray:
+def getNumpyGlobalIdsArray( data: Union[ vtkCellData, vtkPointData ], ) -> npt.NDArray:
     """Get a numpy array of the GlobalIds if it exist.
 
+    Note that for some cases (GEOS simulations), the attribute "localToGlobalMap" is the GlobalIds.
+
     Args:
         data (Union[ vtkCellData, vtkPointData ]): Cell or point array.
 
@@ -298,14 +304,19 @@ def getNumpyGlobalIdsArray( data: Union[ vtkCellData, vtkPointData ] ) -> npt.ND
     if not isinstance( data, vtkFieldData ):
         raise TypeError( f"data '{ data }' entered is not a vtkFieldData object." )
 
-    global_ids: Optional[ vtkDataArray ] = data.GetGlobalIds()
-    if global_ids is None:
+    globalIds: vtkDataArray = data.GetGlobalIds() if data.GetGlobalIds() is not None else data.GetArray(
+        "localToGlobalMap" )
+    if globalIds is None:
         raise AttributeError( "There is no GlobalIds in the given fieldData." )
 
-    return vtk_to_numpy( global_ids )
+    return vtk_to_numpy( globalIds )
 
 
-def getNumpyArrayByName( data: Union[ vtkCellData, vtkPointData ], name: str, sorted: bool = False ) -> npt.NDArray:
+def getNumpyArrayByName(
+    data: Union[ vtkCellData, vtkPointData ],
+    name: str,
+    sorted: bool = False,
+) -> npt.NDArray:
     """Get the numpy array of a given vtkDataArray found by its name.
 
     If sorted is selected, this allows the option to reorder the values wrt GlobalIds. If not GlobalIds was found,
@@ -333,7 +344,10 @@ def getNumpyArrayByName( data: Union[ vtkCellData, vtkPointData ], name: str, so
     return npArray
 
 
-def getAttributeSet( mesh: Union[ vtkMultiBlockDataSet, vtkDataSet ], piece: Piece ) -> set[ str ]:
+def getAttributeSet(
+    mesh: Union[ vtkMultiBlockDataSet, vtkDataSet ],
+    piece: Piece,
+) -> set[ str ]:
     """Get the set of all attributes from an mesh on points or on cells.
 
     Args:
@@ -400,7 +414,11 @@ def getAttributesWithNumberOfComponents(
     return attributes
 
 
-def isAttributeInObject( mesh: Union[ vtkMultiBlockDataSet, vtkDataSet ], attributeName: str, piece: Piece ) -> bool:
+def isAttributeInObject(
+    mesh: Union[ vtkMultiBlockDataSet, vtkDataSet ],
+    attributeName: str,
+    piece: Piece,
+) -> bool:
     """Check if an attribute is in the input mesh for the given piece.
 
     Args:
@@ -437,7 +455,11 @@ def isAttributeInObject( mesh: Union[ vtkMultiBlockDataSet, vtkDataSet ], attrib
     return False
 
 
-def isAttributeGlobal( multiBlockDataSet: vtkMultiBlockDataSet, attributeName: str, piece: Piece ) -> bool:
+def isAttributeGlobal(
+    multiBlockDataSet: vtkMultiBlockDataSet,
+    attributeName: str,
+    piece: Piece,
+) -> bool:
     """Check if an attribute is global in the input multiBlockDataSet.
 
     Args:
@@ -457,7 +479,11 @@ def isAttributeGlobal( multiBlockDataSet: vtkMultiBlockDataSet, attributeName: s
     return True
 
 
-def getArrayInObject( dataSet: vtkDataSet, attributeName: str, piece: Piece ) -> npt.NDArray[ Any ]:
+def getArrayInObject(
+    dataSet: vtkDataSet,
+    attributeName: str,
+    piece: Piece,
+) -> npt.NDArray[ Any ]:
     """Return the numpy array corresponding to input attribute name in the mesh.
 
     Args:
@@ -474,7 +500,11 @@ def getArrayInObject( dataSet: vtkDataSet, attributeName: str, piece: Piece ) ->
     return npArray
 
 
-def getVtkArrayTypeInObject( mesh: Union[ vtkDataSet, vtkMultiBlockDataSet ], attributeName: str, piece: Piece ) -> int:
+def getVtkArrayTypeInObject(
+    mesh: Union[ vtkDataSet, vtkMultiBlockDataSet ],
+    attributeName: str,
+    piece: Piece,
+) -> int:
     """Return VTK type of requested array from input mesh.
 
     Args:
@@ -506,7 +536,11 @@ def getVtkArrayTypeInObject( mesh: Union[ vtkDataSet, vtkMultiBlockDataSet ], at
         raise TypeError( "Input object must be a vtkDataSet or vtkMultiBlockDataSet." )
 
 
-def getVtkArrayInObject( dataSet: vtkDataSet, attributeName: str, piece: Piece ) -> vtkDataArray:
+def getVtkArrayInObject(
+    dataSet: vtkDataSet,
+    attributeName: str,
+    piece: Piece,
+) -> vtkDataArray:
     """Return the array corresponding to input attribute name in the mesh.
 
     Args:
@@ -621,9 +655,11 @@ def getComponentNames(
     return tuple( componentNames )
 
 
-def getAttributeValuesAsDF( surface: vtkPolyData,
-                            attributeNames: tuple[ str, ...],
-                            piece: Piece = Piece.CELLS ) -> pd.DataFrame:
+def getAttributeValuesAsDF(
+    surface: vtkPolyData,
+    attributeNames: tuple[ str, ...],
+    piece: Piece = Piece.CELLS,
+) -> pd.DataFrame:
     """Get attribute values from input surface.
 
     Args:
@@ -660,7 +696,7 @@ def getAttributeValuesAsDF( surface: vtkPolyData,
     return data
 
 
-def computeCellCenterCoordinates( mesh: vtkDataSet ) -> vtkDataArray:
+def computeCellCenterCoordinates( mesh: vtkDataSet, ) -> vtkDataArray:
     """Get the coordinates of Cell center.
 
     Args:

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

@@ -4,47 +4,26 @@
 import logging
 import numpy as np
 import numpy.typing as npt
-import vtkmodules.util.numpy_support as vnp
+
 from typing import Union, Any
-from geos.utils.Logger import ( getLogger, Logger, VTKCaptureLog, RegexExceptionFilter )
+from random import randint
 
-from vtk import (  # type: ignore[import-untyped]
-    VTK_BIT, VTK_UNSIGNED_CHAR, VTK_UNSIGNED_SHORT, VTK_UNSIGNED_LONG, VTK_UNSIGNED_INT, VTK_UNSIGNED_LONG_LONG,
-    VTK_CHAR, VTK_SIGNED_CHAR, VTK_SHORT, VTK_LONG, VTK_INT, VTK_LONG_LONG, VTK_ID_TYPE, VTK_FLOAT, VTK_DOUBLE,
-)
-from vtkmodules.vtkCommonDataModel import (
-    vtkMultiBlockDataSet,
-    vtkDataSet,
-    vtkPointSet,
-    vtkCompositeDataSet,
-    vtkDataObject,
-    vtkPointData,
-    vtkCellData,
-    vtkCell,
-)
-from vtkmodules.vtkFiltersCore import (
-    vtkArrayRename,
-    vtkCellCenters,
-    vtkPointDataToCellData,
-)
-from vtkmodules.vtkCommonCore import (
-    vtkDataArray,
-    vtkPoints,
-    vtkLogger,
-)
-from geos.mesh.utils.arrayHelpers import (
-    getComponentNames,
-    getAttributesWithNumberOfComponents,
-    getArrayInObject,
-    isAttributeInObject,
-    isAttributeGlobal,
-    getVtkArrayTypeInObject,
-    getNumberOfComponents,
-)
-from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten
-from geos.utils.Errors import VTKError
+from vtkmodules.vtkFiltersCore import vtkArrayRename, vtkCellCenters, vtkPointDataToCellData
+from vtkmodules.util.numpy_support import vtk_to_numpy, numpy_to_vtk, get_vtk_to_numpy_typemap
+from vtkmodules.vtkCommonCore import ( vtkDataArray, vtkPoints, vtkLogger, VTK_BIT, VTK_UNSIGNED_CHAR,
+                                       VTK_UNSIGNED_SHORT, VTK_UNSIGNED_LONG, VTK_UNSIGNED_INT, VTK_UNSIGNED_LONG_LONG,
+                                       VTK_CHAR, VTK_SIGNED_CHAR, VTK_SHORT, VTK_LONG, VTK_INT, VTK_LONG_LONG,
+                                       VTK_ID_TYPE, VTK_FLOAT, VTK_DOUBLE )
+from vtkmodules.vtkCommonDataModel import ( vtkMultiBlockDataSet, vtkDataSet, vtkPointSet, vtkCompositeDataSet,
+                                            vtkDataObject, vtkPointData, vtkCellData, vtkCell )
 
+from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten
+from geos.mesh.utils.arrayHelpers import ( getComponentNames, getAttributesWithNumberOfComponents, getArrayInObject,
+                                           isAttributeInObject, isAttributeGlobal, getVtkArrayTypeInObject,
+                                           getNumberOfComponents )
 from geos.utils.pieceEnum import Piece
+from geos.utils.Errors import VTKError
+from geos.utils.Logger import ( getLogger, Logger, VTKCaptureLog, RegexExceptionFilter )
 
 __doc__ = """
 ArrayModifiers contains utilities to process VTK Arrays objects.
@@ -111,7 +90,7 @@ def fillPartialAttributes(
     if nbComponents > 1:
         componentNames = getComponentNames( multiBlockDataSet, attributeName, piece )
 
-    typeMapping: dict[ int, type ] = vnp.get_vtk_to_numpy_typemap()
+    typeMapping: dict[ int, type ] = get_vtk_to_numpy_typemap()
     valueType: type = typeMapping[ vtkDataType ]
     # Set the default value depending of the type of the attribute to fill
     if listValues is None:
@@ -220,7 +199,7 @@ def createEmptyAttribute(
         vtkDataArray: The empty attribute.
     """
     # Check if the vtk data type is correct.
-    vtkNumpyTypeMap: dict[ int, type ] = vnp.get_vtk_to_numpy_typemap()
+    vtkNumpyTypeMap: dict[ int, type ] = get_vtk_to_numpy_typemap()
     if vtkDataType not in vtkNumpyTypeMap:
         raise ValueError( f"Attribute type { vtkDataType } is unknown." )
 
@@ -312,7 +291,7 @@ def createConstantAttribute(
         # Check the consistency between the given value type and the vtk array type if it exists.
         valueType = valueType().dtype
         if vtkDataType is not None:
-            vtkNumpyTypeMap: dict[ int, type ] = vnp.get_vtk_to_numpy_typemap()
+            vtkNumpyTypeMap: dict[ int, type ] = get_vtk_to_numpy_typemap()
             if vtkDataType not in vtkNumpyTypeMap:
                 raise ValueError( f"The vtk data type { vtkDataType } is unknown." )
 
@@ -394,7 +373,7 @@ def createAttribute(
 
     # Check the coherency between the given array type and the vtk array type if it exist.
     if vtkDataType is not None:
-        vtkNumpyTypeMap: dict[ int, type ] = vnp.get_vtk_to_numpy_typemap()
+        vtkNumpyTypeMap: dict[ int, type ] = get_vtk_to_numpy_typemap()
         if vtkDataType not in vtkNumpyTypeMap:
             raise ValueError( f"The vtk data type { vtkDataType } is unknown." )
 
@@ -417,7 +396,7 @@ def createAttribute(
         raise ValueError( f"The npArray must have { nbElements } elements, not { len( npArray ) }." )
 
     # Convert the numpy array int a vtkDataArray.
-    createdAttribute: vtkDataArray = vnp.numpy_to_vtk( npArray, deep=True, array_type=vtkDataType )
+    createdAttribute: vtkDataArray = numpy_to_vtk( npArray, deep=True, array_type=vtkDataType )
     createdAttribute.SetName( attributeName )
 
     nbComponents: int = createdAttribute.GetNumberOfComponents()
@@ -459,7 +438,7 @@ def copyAttribute(
     The similarity of two meshes means that the two mesh have the same number of elements (cells and points) located in the same coordinates and with the same indexation. Testing this similarity is time consuming therefore, only few metric are compared:
         - the block indexation for multiblock dataset
         - the number of the element where the attribute is located, for multiblock dataset it is done for each block
-        - the coordinates of the first element, for multiblock dataset it is done for each block
+        - the coordinates of a random element, for multiblock dataset it is done for each block
 
     Args:
         meshFrom (vtkMultiBlockDataSet | vtkDataSet): mesh from which to copy the attribute.
@@ -484,12 +463,27 @@ def copyAttribute(
         # Small check to check if the two meshes are similar.
         coordElementTo: set[ tuple[ float, ...] ] = set()
         coordElementFrom: set[ tuple[ float, ...] ] = set()
+        nbElementTo: int
+        nbElementFrom: int
+        idElemToCompare: int
         if piece == Piece.POINTS:
-            coordElementTo.add( meshTo.GetPoint( 0 ) )
-            coordElementFrom.add( meshFrom.GetPoint( 0 ) )
+            nbElementTo = meshTo.GetNumberOfPoints()
+            nbElementFrom = meshFrom.GetNumberOfPoints()
+            if nbElementFrom != nbElementTo:
+                raise ValueError( "The two meshes have not the same number of points." )
+
+            idElemToCompare = randint( 0, nbElementTo )
+            coordElementTo.add( meshTo.GetPoint( idElemToCompare ) )
+            coordElementFrom.add( meshFrom.GetPoint( idElemToCompare ) )
         elif piece == Piece.CELLS:
-            cellTo: vtkCell = meshTo.GetCell( 0 )
-            cellFrom: vtkCell = meshFrom.GetCell( 0 )
+            nbElementTo = meshTo.GetNumberOfCells()
+            nbElementFrom = meshFrom.GetNumberOfCells()
+            if nbElementFrom != nbElementTo:
+                raise ValueError( "The two meshes have not the same number of cells." )
+
+            idElemToCompare = randint( 0, nbElementTo )
+            cellTo: vtkCell = meshTo.GetCell( idElemToCompare )
+            cellFrom: vtkCell = meshFrom.GetCell( idElemToCompare )
             # Get the coordinates of each points of the cell.
             nbPointsTo: int = cellTo.GetNumberOfPoints()
             nbPointsFrom: int = cellTo.GetNumberOfPoints()
@@ -503,6 +497,7 @@ def copyAttribute(
                 coordElementFrom.add( cellPointsFrom.GetPoint( idPoint ) )
         else:
             raise ValueError( "The piece of the attribute to copy must be cells or points." )
+
         if coordElementTo != coordElementFrom:
             raise ValueError( "The two meshes have not the same element indexation." )
 
@@ -596,6 +591,9 @@ def transferAttributeWithElementMap(
     ## it is global
     if isinstance( meshFrom, vtkMultiBlockDataSet ) and not isAttributeGlobal( meshFrom, attributeName, piece ):
         raise AttributeError( f"The attribute { attributeName } must be global in the source mesh." )
+    ## it is not in the meshTo
+    if isAttributeInObject( meshTo, attributeName, piece ):
+        raise AttributeError( f"The attribute { attributeName } is already in the final mesh." )
 
     # Transfer the attribute
     if isinstance( meshTo, vtkDataSet ):
@@ -624,7 +622,7 @@ def transferAttributeWithElementMap(
         else:
             raise AttributeError( f"The attribute { attributeName } has an unknown type." )
 
-        typeMapping: dict[ int, type ] = vnp.get_vtk_to_numpy_typemap()
+        typeMapping: dict[ int, type ] = get_vtk_to_numpy_typemap()
         valueType: type = typeMapping[ vtkDataType ]
 
         arrayTo: npt.NDArray[ Any ]
@@ -648,7 +646,7 @@ def transferAttributeWithElementMap(
                 else:
                     raise TypeError( "The source mesh has to be inherited from vtkDataSet or vtkMultiBlockDataSet." )
 
-                arrayFrom: npt.NDArray[ Any ] = vnp.vtk_to_numpy( dataFrom.GetArray( attributeName ) )
+                arrayFrom: npt.NDArray[ Any ] = vtk_to_numpy( dataFrom.GetArray( attributeName ) )
                 valueToTransfer = arrayFrom[ idElementFrom ]
 
             arrayTo[ idElementTo ] = valueToTransfer
@@ -661,9 +659,6 @@ def transferAttributeWithElementMap(
                          vtkDataType=vtkDataType,
                          logger=logger )
     elif isinstance( meshTo, vtkMultiBlockDataSet ):
-        if isAttributeInObject( meshTo, attributeName, piece ):
-            raise AttributeError( f"The attribute { attributeName } is already in the final mesh." )
-
         listFlatIdDataSetTo: list[ int ] = getBlockElementIndexesFlatten( meshTo )
         for flatIdDataSetTo in listFlatIdDataSetTo:
             dataSetTo: vtkDataSet = vtkDataSet.SafeDownCast( meshTo.GetDataSet( flatIdDataSetTo ) )