Merge branch '351-damask-grid-load_spparks' into 'development'

load data from SPPARKS

Closes #351

See merge request damask/DAMASK!861

python/damask/_geomgrid.py

@@ -197,10 +197,41 @@ class GeomGrid:
 
 
     @staticmethod
-    def load(fname: Union[str, Path]) -> 'GeomGrid':
+    def _load(fname: Union[str, Path], label: str) -> 'GeomGrid':
         """
         Load from VTK ImageData file.
 
+        Parameters
+        ----------
+        fname : str or pathlib.Path
+            VTK ImageData file to read.
+            Valid extension is .vti, which will be appended if not given.
+        label : str
+            Label of the dataset containing the material IDs.
+
+        Returns
+        -------
+        loaded : damask.GeomGrid
+            GeomGrid-based geometry from file.
+
+        """
+        v = VTK.load(fname if str(fname).endswith('.vti') else str(fname)+'.vti')
+        cells = np.array(v.vtk_data.GetDimensions())-1
+        bbox  = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
+        ic = {l:v.get(l).reshape(cells,order='F') for l in set(v.labels['Cell Data']) - {label}}
+
+        return GeomGrid(material = v.get(label).reshape(cells,order='F'),
+                        size     = bbox[1] - bbox[0],
+                        origin   = bbox[0],
+                        initial_conditions = ic,
+                        comments = v.comments,
+                       )
+
+    @staticmethod
+    def load(fname: Union[str, Path]) -> 'GeomGrid':
+        """
+        Load from VTK ImageData file with material IDs stored as 'material'.
+
         Parameters
         ----------
         fname : str or pathlib.Path
@@ -213,17 +244,32 @@ class GeomGrid:
             GeomGrid-based geometry from file.
 
         """
-        v = VTK.load(fname if str(fname).endswith('.vti') else str(fname)+'.vti')
-        cells = np.array(v.vtk_data.GetDimensions())-1
-        bbox  = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
-        ic = {label:v.get(label).reshape(cells,order='F') for label in set(v.labels['Cell Data']) - {'material'}}
+        return GeomGrid._load(fname,'material')
 
-        return GeomGrid(material = v.get('material').reshape(cells,order='F'),
-                        size     = bbox[1] - bbox[0],
-                        origin   = bbox[0],
-                        initial_conditions = ic,
-                        comments = v.comments,
-                       )
+
+    @staticmethod
+    def load_SPPARKS(fname: Union[str, Path]) -> 'GeomGrid':
+        """
+        Load from SPPARKS VTK dump.
+
+        Parameters
+        ----------
+        fname : str or pathlib.Path
+            SPPARKS VTK dump file to read.
+            Valid extension is .vti, which will be appended if not given.
+
+        Returns
+        -------
+        loaded : damask.GeomGrid
+            GeomGrid-based geometry from file.
+
+        Notes
+        -----
+        A SPPARKS VTI dump is equivalent to a DAMASK VTI file,
+        but stores the materialID information as 'Spin' rather than 'material'.
+
+        """
+        return GeomGrid._load(fname,'Spin')
 
 
     @typing.no_type_check

python/damask/_vtk.py

@@ -457,7 +457,7 @@ class VTK:
             data: Union[None, np.ndarray, np.ma.MaskedArray] = None,
             info: Optional[str] = None,
             *,
-            table: Optional['Table'] = None):
+            table: Optional['Table'] = None) -> 'VTK':
         """
         Add new or replace existing point or cell data.
 
@@ -534,7 +534,6 @@ class VTK:
         else:
             raise TypeError
 
-
         return dup
 
 
@@ -559,20 +558,18 @@ class VTK:
 
         """
         cell_data = self.vtk_data.GetCellData()
-        for a in range(cell_data.GetNumberOfArrays()):
-            if cell_data.GetArrayName(a) == label:
-                try:
-                    return vtk_to_numpy(cell_data.GetArray(a))
-                except AttributeError:
-                    vtk_array = cell_data.GetAbstractArray(a)                                       # string array
+        if label in [cell_data.GetArrayName(a) for a in range(cell_data.GetNumberOfArrays())]:
+            try:
+                return vtk_to_numpy(cell_data.GetArray(label))
+            except AttributeError:
+                vtk_array = cell_data.GetAbstractArray(label)                                       # string array
 
         point_data = self.vtk_data.GetPointData()
-        for a in range(point_data.GetNumberOfArrays()):
-            if point_data.GetArrayName(a) == label:
-                try:
-                    return vtk_to_numpy(point_data.GetArray(a))
-                except AttributeError:
-                    vtk_array = point_data.GetAbstractArray(a)                                      # string array
+        if label in [point_data.GetArrayName(a) for a in range(point_data.GetNumberOfArrays())]:
+            try:
+                return vtk_to_numpy(point_data.GetArray(label))
+            except AttributeError:
+                vtk_array = point_data.GetAbstractArray(label)                                      # string array
 
         try:
             # string array
@@ -581,6 +578,41 @@ class VTK:
             raise KeyError(f'array "{label}" not found')
 
 
+    def delete(self,
+               label: str) -> 'VTK':
+        """
+        Delete either cell or point data.
+
+        Cell data takes precedence over point data, i.e. this
+        function assumes that labels are unique among cell and
+        point data.
+
+        Parameters
+        ----------
+        label : str
+            Data label.
+
+        Returns
+        -------
+        updated : damask.VTK
+            Updated VTK-based geometry.
+
+        """
+        dup = self.copy()
+
+        cell_data = dup.vtk_data.GetCellData()
+        if label in [cell_data.GetArrayName(a) for a in range(cell_data.GetNumberOfArrays())]:
+            dup.vtk_data.GetCellData().RemoveArray(label)
+            return dup
+
+        point_data = self.vtk_data.GetPointData()
+        if label in [point_data.GetArrayName(a) for a in range(point_data.GetNumberOfArrays())]:
+            dup.vtk_data.GetPointData().RemoveArray(label)
+            return dup
+
+        raise KeyError(f'array "{label}" not found')
+
+
     def show(self,
              label: Optional[str] = None,
              colormap: Union[Colormap, str] = 'cividis'):

python/tests/test_GeomGrid.py

@@ -92,6 +92,12 @@ class TestGeomGrid:
         default.material -= 1
         assert GeomGrid.load_ASCII(tmp_path/'ASCII') == default
 
+    def test_save_load_SPPARKS(self,res_path,tmp_path):
+        v = VTK.load(res_path/'SPPARKS_dump.vti')
+        v.set('material',v.get('Spin')).delete('Spin').save(tmp_path/'SPPARKS_dump.vti',parallel=False)
+        assert GeomGrid.load_SPPARKS(res_path/'SPPARKS_dump.vti') == \
+               GeomGrid.load(tmp_path/'SPPARKS_dump.vti')
+
     def test_invalid_origin(self,default):
         with pytest.raises(ValueError):
             GeomGrid(default.material[1:,1:,1:],

python/tests/test_VTK.py

@@ -199,6 +199,12 @@ class TestVTK:
         mask_manual = default.set('D',np.where(masked.mask,masked.fill_value,masked))
         assert mask_manual == mask_auto
 
+    @pytest.mark.parametrize('mode',['cells','points'])
+    def test_delete(self,default,mode):
+        data = np.random.rand(default.N_cells if mode == 'cells' else default.N_points).astype(np.float32)
+        v = default.set('D',data)
+        assert (v.get('D') == data).all()
+        assert v.delete('D') == default
 
     @pytest.mark.parametrize('data_type,shape',[(float,(3,)),
                                                 (float,(3,3)),