Merge branch 'numpy-2.0-compatible' into 'development'

compatible with NumPy 2.0

See merge request damask/DAMASK!953

python/damask/_geomgrid.py

@@ -132,12 +132,12 @@ def material(self,
                  material: np.ndarray):
         if len(material.shape) != 3:
             raise ValueError(f'invalid material shape {material.shape}')
-        if material.dtype not in np.sctypes['float'] and material.dtype not in np.sctypes['int']:
+        if material.dtype not in [np.float32,np.float64, np.int32,np.int64]:
             raise TypeError(f'invalid material data type "{material.dtype}"')
 
         self._material = np.copy(material)
 
-        if self.material.dtype in np.sctypes['float'] and \
+        if self.material.dtype in [np.float32,np.float64] and \
            np.all(self.material == self.material.astype(np.int64).astype(float)):
             self._material = self.material.astype(np.int64)
 
@@ -153,7 +153,7 @@ def size(self,
         if len(size) != 3 or any(np.array(size) < 0):
             raise ValueError(f'invalid size {size}')
 
-        self._size = np.array(size)
+        self._size = np.array(size,np.float64)
 
     @property
     def origin(self) -> np.ndarray:
@@ -166,7 +166,7 @@ def origin(self,
         if len(origin) != 3:
             raise ValueError(f'invalid origin {origin}')
 
-        self._origin = np.array(origin)
+        self._origin = np.array(origin,np.float64)
 
     @property
     def initial_conditions(self) -> Dict[str,np.ndarray]:
@@ -471,7 +471,7 @@ def load_DREAM3D(fname: Union[str, Path],
                 O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa
                 unique,unique_inverse = np.unique(np.hstack([O,phase]),return_inverse=True,axis=0)
                 ma = np.arange(cells.prod()) if len(unique) == cells.prod() else \
-                     np.arange(unique.size)[np.argsort(pd.unique(unique_inverse))][unique_inverse]
+                     np.arange(unique.size)[np.argsort(pd.unique(unique_inverse.squeeze()))][unique_inverse]
             else:
                 ma = f['/'.join([b,c,feature_IDs])][()].flatten()
 
@@ -512,7 +512,7 @@ def from_table(table: Table,
         unique,unique_inverse = np.unique(np.hstack([table.get(l) for l in labels_]),return_inverse=True,axis=0)
 
         ma = np.arange(cells.prod()) if len(unique) == cells.prod() else \
-             np.arange(unique.size)[np.argsort(pd.unique(unique_inverse))][unique_inverse]
+             np.arange(unique.size)[np.argsort(pd.unique(unique_inverse.squeeze()))][unique_inverse]
 
         return GeomGrid(material = ma.reshape(cells,order='F'),
                         size     = size,
@@ -812,7 +812,7 @@ def save_ASCII(self,
                    'homogenization 1',
                   ]
 
-        format_string = '%g' if self.material.dtype in np.sctypes['float'] else \
+        format_string = '%g' if self.material.dtype in [np.float32,np.float64] else \
                         '%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.material)))))
         np.savetxt(fname,
                    self.material.reshape([self.cells[0],np.prod(self.cells[1:])],order='F').T,
@@ -1347,14 +1347,14 @@ def add_primitive(self,
 
         """
         # radius and center
-        r = np.array(dimension)/2.0*self.size/self.cells if np.array(dimension).dtype in np.sctypes['int'] else \
+        r = np.array(dimension)/2.0*self.size/self.cells if np.issubdtype(np.array(dimension).dtype,np.integer) else \
             np.array(dimension)/2.0
-        c = (np.array(center) + .5)*self.size/self.cells if np.array(center).dtype    in np.sctypes['int'] else \
+        c = (np.array(center) + .5)*self.size/self.cells if np.issubdtype(np.array(center).dtype,   np.integer) else \
             (np.array(center) - self.origin)
 
         coords = grid_filters.coordinates0_point(self.cells,self.size,
                                           -(0.5*(self.size + (self.size/self.cells
-                                                              if np.array(center).dtype in np.sctypes['int'] else
+                                                              if np.issubdtype(np.array(center).dtype,np.integer) else
                                                               0)) if periodic else c))
         coords_rot = R.broadcast_to(tuple(self.cells))@coords
 

python/damask/_result.py

@@ -25,11 +25,14 @@
 from . import util
 from ._typehints import FloatSequence, IntSequence, DADF5Dataset
 
+
 h5py3 = h5py.__version__[0] == '3'
 
 chunk_size = 1024**2//8                                                                             # for compression in HDF5
+
 prefix_inc = 'increment_'
 
+
 def _read(dataset: h5py._hl.dataset.Dataset) -> np.ndarray:
     """Read a dataset and its metadata into a numpy.ndarray."""
     metadata = {k:(v.decode() if not h5py3 and type(v) is bytes else v) for k,v in dataset.attrs.items()}
@@ -59,7 +62,7 @@ def _empty_like(dataset: np.ma.core.MaskedArray,
                 fill_int: int) -> np.ma.core.MaskedArray:
     """Create empty numpy.ma.MaskedArray."""
     return ma.array(np.empty((N_materialpoints,)+dataset.shape[1:],dataset.dtype),
-                    fill_value = fill_float if dataset.dtype in np.sctypes['float'] else fill_int,
+                    fill_value = fill_float if np.issubdtype(dataset.dtype,np.floating) else fill_int,
                     mask = True)
 
 
@@ -1783,9 +1786,10 @@ def export_XDMF(self,
         attribute_type_map = defaultdict(lambda:'Matrix', ( ((),'Scalar'), ((3,),'Vector'), ((3,3),'Tensor')) )
 
         def number_type_map(dtype):
-            if dtype in np.sctypes['int']:   return 'Int'
-            if dtype in np.sctypes['uint']:  return 'UInt'
-            if dtype in np.sctypes['float']: return 'Float'
+            if np.issubdtype(dtype,np.signedinteger):   return 'Int'
+            if np.issubdtype(dtype,np.unsignedinteger): return 'UInt'
+            if np.issubdtype(dtype,np.floating):        return 'Float'
+            raise TypeError(f'invalid type "{dtype}"')
 
 
         xdmf = ET.Element('Xdmf')

python/damask/util.py

@@ -23,6 +23,7 @@
 from . import version as _version
 from ._typehints import FloatSequence as _FloatSequence, NumpyRngSeed as _NumpyRngSeed, FileHandle as _FileHandle
 
+
 # https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py
 # https://stackoverflow.com/questions/287871
 _colors = {
@@ -323,8 +324,9 @@ def lcm(a,b):
 
     max_denominator = int(10**(N_significant-1))
 
-    if (v_ := _np.asarray(v)).dtype in _np.sctypes['float']:
-        v_ = _np.round(_np.asarray(v,'float64')/_np.max(_np.abs(v)),N_significant)
+    v_ = _np.asarray(v)
+    if _np.issubdtype(v_.dtype,_np.inexact):
+        v_ = _np.round(_np.asarray(v,_np.float64)/_np.max(_np.abs(v)),N_significant)
     m = (v_ * _reduce(lcm, map(lambda x: int(get_square_denominator(x,max_denominator)),v_))**0.5).astype(_np.int64)
     m = m//_reduce(_np.gcd,m)
 

python/tests/test_GeomGrid.py

@@ -80,8 +80,8 @@ def test_invalid_material_type(self):
             GeomGrid(np.zeros((3,3,3),dtype='complex'),np.ones(3))
 
     def test_cast_to_int(self):
-        g = GeomGrid(np.zeros((3,3,3)),np.ones(3))
-        assert g.material.dtype in np.sctypes['int']
+        g = GeomGrid(np.zeros((3,3,3),dtype=np.float64),np.ones(3))
+        assert g.material.dtype in [np.int32,np.int64]
 
     def test_invalid_size(self,default):
         with pytest.raises(ValueError):

python/tests/test_grid_filters.py

@@ -394,9 +394,9 @@ def test_div(self,field_def,div_def):
     def test_ravel_index(self):
         cells = np.random.randint(8,32,(3))
 
-        indices = np.block(np.meshgrid(np.arange(cells[0]),
-                                       np.arange(cells[1]),
-                                       np.arange(cells[2]),indexing='ij')).reshape(tuple(cells)+(3,),order='F')
+        indices = np.block(list(np.meshgrid(np.arange(cells[0]),
+                                            np.arange(cells[1]),
+                                            np.arange(cells[2]),indexing='ij'))).reshape(tuple(cells)+(3,),order='F')
         x,y,z = map(np.random.randint,cells)
         assert grid_filters.ravel_index(indices)[x,y,z] == np.arange(0,np.prod(cells)).reshape(cells,order='F')[x,y,z]
 

src/grid/VTI.f90

@@ -186,8 +186,8 @@ subroutine VTI_readCellsSizeOrigin(cells,geomSize,origin, &
 
   end do outer
 
-  if (any(geomSize<=0))                 call IO_error(error_ID = 844, ext_msg='size')
-  if (any(cells<1))                     call IO_error(error_ID = 844, ext_msg='cells')
+  if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
+  if (any(cells<1))     call IO_error(error_ID = 844, ext_msg='cells')
 
 end subroutine VTI_readCellsSizeOrigin