grid assemble + corrected grid.scale

python/damask/_grid.py

@@ -10,7 +10,7 @@ from pathlib import Path
 import numpy as np
 import pandas as pd
 import h5py
-from scipy import ndimage, spatial
+from scipy import ndimage, spatial, interpolate
 
 from . import VTK
 from . import util
@@ -41,7 +41,7 @@ class Grid:
 
         Parameters
         ----------
-        material : numpy.ndarray, shape (:,:,:)
+        material : numpy.ndarray of int, shape (:,:,:)
             Material indices. The shape of the material array defines
             the number of cells.
         size : sequence of float, len (3)
@@ -50,7 +50,7 @@ class Grid:
             Coordinates of grid origin in meter. Defaults to [0.0,0.0,0.0].
         initial_conditions : dictionary, optional
             Labels and values of the inital conditions at each material point.
-        comments : str or iterable of str, optional
+        comments : str or sequence of str, optional
             Additional, human-readable information, e.g. history of operations.
 
         """
@@ -183,7 +183,7 @@ class Grid:
     @comments.setter
     def comments(self,
                  comments: Union[str, Sequence[str]]):
-        self._comments = [str(c) for c in comments] if isinstance(comments,list) else [str(comments)]
+        self._comments = [str(c) for c in comments] if isinstance(comments,Sequence) else [str(comments)]
 
 
     @property
@@ -427,7 +427,7 @@ class Grid:
         coordinates : str
             Label of the vector column containing the spatial coordinates.
             Need to be ordered (1./x fast, 3./z slow).
-        labels : (list of) str
+        labels : str or sequence of str
             Label(s) of the columns containing the material definition.
             Each unique combination of values results in one material ID.
 
@@ -474,7 +474,7 @@ class Grid:
             Number of cells in x,y,z direction.
         size : sequence of float, len (3)
             Physical size of the grid in meter.
-        seeds : numpy.ndarray, shape (:,3)
+        seeds : numpy.ndarray of float, shape (:,3)
             Position of the seed points in meter. All points need to lay within the box.
         weights : sequence of float, len (seeds.shape[0])
             Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
@@ -531,7 +531,7 @@ class Grid:
             Number of cells in x,y,z direction.
         size : sequence of float, len (3)
             Physical size of the grid in meter.
-        seeds : numpy.ndarray, shape (:,3)
+        seeds : numpy.ndarray of float, shape (:,3)
             Position of the seed points in meter. All points need to lay within the box.
         material : sequence of int, len (seeds.shape[0]), optional
             Material ID of the seeds.
@@ -940,17 +940,14 @@ class Grid:
 
 
     def scale(self,
-              cells: IntSequence,
+              cells: IntSequence) -> 'Grid':
-              periodic: bool = True) -> 'Grid':
         """
-        Scale grid to new cells.
+        Scale grid to new cell count.
 
         Parameters
         ----------
         cells : sequence of int, len (3)
             Number of cells in x,y,z direction.
-        periodic : bool, optional
-            Assume grid to be periodic. Defaults to True.
 
         Returns
         -------
@@ -963,7 +960,11 @@ class Grid:
 
         >>> import numpy as np
         >>> import damask
-        >>> g = damask.Grid(np.zeros([32]*3,int),np.ones(3)*1e-4)
+        >>> (g := damask.Grid(np.zeros([32]*3,int),np.ones(3)*1e-4))
+        cells:  32 × 32 × 32
+        size:   0.0001 × 0.0001 × 0.0001 m³
+        origin: 0.0   0.0   0.0 m
+        # materials: 1
         >>> g.scale(g.cells*2)
         cells : 64 x 64 x 64
         size  : 0.0001 x 0.0001 x 0.0001 m³
@@ -971,20 +972,49 @@ class Grid:
         # materials: 1
 
         """
-        return Grid(material = ndimage.interpolation.zoom(
+        options = ('nearest',False,None)
-                                                          self.material,
+        orig = tuple(map(np.linspace,self.origin             + self.size/self.cells*.5,
-                                                          cells/self.cells,
+                                     self.origin + self.size - self.size/self.cells*.5,self.cells))
-                                                          output=self.material.dtype,
+        new = grid_filters.coordinates0_point(cells,self.size,self.origin)
-                                                          order=0,
+
-                                                          mode='wrap' if periodic else 'nearest',
+        return Grid(material = interpolate.RegularGridInterpolator(orig,self.material,*options)(new).astype(int),
-                                                          prefilter=False
-                                                         ),
                     size     = self.size,
                     origin   = self.origin,
+                    initial_conditions = {k: interpolate.RegularGridInterpolator(orig,v,*options)(new)
+                                          for k,v in self.initial_conditions.items()},
                     comments = self.comments+[util.execution_stamp('Grid','scale')],
                    )
 
 
+    def assemble(self,
+                 idx: np.ndarray) -> 'Grid':
+        """
+        Assemble new grid from index map.
+
+        Parameters
+        ----------
+        idx : numpy.ndarray of int, shape (:,:,:) or (:,:,:,3)
+          Grid of flat indices or coordinate indices.
+
+        Returns
+        -------
+        updated : damask.Grid
+            Updated grid-based geometry.
+            Cell count of resulting grid matches shape of index map.
+
+        """
+        cells = idx.shape[:3]
+        flat = (idx if len(idx.shape)==3 else grid_filters.ravel_index(idx)).flatten(order='F')
+        ic = {k: v.flatten(order='F')[flat].reshape(cells,order='F') for k,v in self.initial_conditions.items()}
+
+        return Grid(material = self.material.flatten(order='F')[flat].reshape(cells,order='F'),
+                    size     = self.size,
+                    origin   = self.origin,
+                    initial_conditions = ic,
+                    comments = self.comments+[util.execution_stamp('Grid','assemble')],
+                   )
+
+
     def renumber(self) -> 'Grid':
         """
         Renumber sorted material indices as 0,...,N-1.
@@ -1113,7 +1143,7 @@ class Grid:
         selection_ = None if selection is None else \
                      set(self.material.flatten()) - set(util.aslist(selection)) if invert_selection else \
                      set(self.material.flatten()) & set(util.aslist(selection))
-        material = ndimage.filters.generic_filter(
+        material = ndimage.generic_filter(
                                           self.material,
                                           most_frequent,
                                           footprint=footprint,
@@ -1269,7 +1299,7 @@ class Grid:
         selection_ = None if selection is None else \
                      set(self.material.flatten()) - set(util.aslist(selection)) if invert_selection else \
                      set(self.material.flatten()) & set(util.aslist(selection))
-        mask = ndimage.filters.generic_filter(self.material,
+        mask = ndimage.generic_filter(self.material,
                                       tainted_neighborhood,
                                       footprint=footprint,
                                       mode='wrap' if periodic else 'nearest',

python/damask/grid_filters.py

@@ -541,28 +541,118 @@ def coordinates0_valid(coordinates0: _np.ndarray) -> bool:
         return False
 
 
+def unravel_index(idx: _np.ndarray) -> _np.ndarray:
+    """
+    Convert flat indices to coordinate indices.
+
+    Parameters
+    ----------
+    idx : numpy.ndarray, shape (:,:,:)
+        Grid of flat indices.
+
+    Returns
+    -------
+    unravelled : numpy.ndarray, shape (:,:,:,3)
+        Grid of coordinate indices.
+
+    Examples
+    --------
+    Unravel a linearly increasing sequence of material indices on a 3 × 2 × 1 grid.
+
+    >>> import numpy as np
+    >>> import damask
+    >>> seq = np.arange(6).reshape((3,2,1),order='F')
+    >>> (coord_idx := damask.grid_filters.unravel_index(seq))
+    array([[[[0, 0, 0]],
+
+            [[0, 1, 0]]],
+
+
+           [[[1, 0, 0]],
+
+            [[1, 1, 0]]],
+
+
+           [[[2, 0, 0]],
+
+            [[2, 1, 0]]]])
+    >>> coord_idx[1,1,0]
+    array([1, 1, 0])
+
+    """
+    cells = idx.shape
+    idx_ = _np.expand_dims(idx,3)
+    return _np.block([  idx_ %cells[0],
+                       (idx_//cells[0]) %cells[1],
+                      ((idx_//cells[0])//cells[1])%cells[2]])
+
+def ravel_index(idx: _np.ndarray) -> _np.ndarray:
+    """
+    Convert coordinate indices to flat indices.
+
+    Parameters
+    ----------
+    idx : numpy.ndarray, shape (:,:,:,3)
+        Grid of coordinate indices.
+
+    Returns
+    -------
+    ravelled : numpy.ndarray, shape (:,:,:)
+        Grid of flat indices.
+
+    Examples
+    --------
+    Ravel a reversed sequence of coordinate indices on a 2 × 2 × 1 grid.
+
+    >>> import numpy as np
+    >>> import damask
+    >>> (rev := np.array([[1,1,0],[0,1,0],[1,0,0],[0,0,0]]).reshape((2,2,1,3)))
+    array([[[[1, 1, 0]],
+
+            [[0, 1, 0]]],
+
+
+           [[[1, 0, 0]],
+
+            [[0, 0, 0]]]])
+    >>> (flat_idx := damask.grid_filters.ravel_index(rev))
+    array([[[3],
+            [2]],
+
+           [[1],
+            [0]]])
+
+    """
+    cells = idx.shape[:3]
+    return   idx[:,:,:,0] \
+           + idx[:,:,:,1]*cells[0] \
+           + idx[:,:,:,2]*cells[0]*cells[1]
+
+
 def regrid(size: _FloatSequence,
            F: _np.ndarray,
            cells: _IntSequence) -> _np.ndarray:
     """
-    Return mapping from coordinates in deformed configuration to a regular grid.
+    Map a deformed grid A back to a rectilinear grid B.
+
+    The size of grid B is chosen as the average deformed size of grid A.
 
     Parameters
     ----------
     size : sequence of float, len (3)
-        Physical size.
+        Physical size of grid A.
-    F : numpy.ndarray, shape (:,:,:,3,3), shape (:,:,:,3,3)
+    F : numpy.ndarray, shape (:,:,:,3,3)
-        Deformation gradient field.
+        Deformation gradient field on grid A.
     cells : sequence of int, len (3)
-        Cell count along x,y,z of remapping grid.
+        Cell count along x,y,z of grid B.
+
+    Returns
+    -------
+    idx : numpy.ndarray of int, shape (cells)
+        Flat index of closest point on deformed grid A for each point on grid B.
 
     """
-    c = coordinates_point(size,F)
+    box = _np.dot(_np.average(F,axis=(0,1,2)),size)
-
+    c = coordinates_point(size,F)%box
-    outer = _np.dot(_np.average(F,axis=(0,1,2)),size)
+    tree = _spatial.cKDTree(c.reshape((-1,3),order='F'),boxsize=box)
-    for d in range(3):
+    return tree.query(coordinates0_point(cells,box))[1]
-        c[_np.where(c[:,:,:,d]<0)]        += outer[d]
-        c[_np.where(c[:,:,:,d]>outer[d])] -= outer[d]
-
-    tree = _spatial.cKDTree(c.reshape(-1,3),boxsize=outer)
-    return tree.query(coordinates0_point(cells,outer))[1].flatten()

python/tests/reference/Grid/scale_grid_10-10-10.vti

@@ -3,7 +3,7 @@
   <ImageData WholeExtent="0 10 0 10 0 10" Origin="0 0 0" Spacing="8e-7 5.000000000000001e-7 4e-7" Direction="1 0 0 0 1 0 0 0 1">
     <FieldData>
       <Array type="String" Name="comments" NumberOfTuples="1" format="binary">
-        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wc0/Nz9UrTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFYKEKs=
+        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wcy/KTNErTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFW2EKk=
       </Array>
     </FieldData>
   <Piece Extent="0 10 0 10 0 10">
@@ -11,7 +11,7 @@
     </PointData>
     <CellData>
       <DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
-        AQAAAACAAABAHwAA+QAAAA==eF7t2DcSAjEQRFF28d57t3i4/wUJaCVdpYCqH46Sl/1oklZR+70iDMMwDP+wlHXZMJuyJdtmR3Yl3evJvhyYQzmSY3Mip5LuzeRcLsylXMm1uZFbSfd2ci8P5lGe5Nms5EXSvau8ybv5kE/5Mt/yI+kefc90j75nukffM92j75nu0fdM9+h7pnv0PdM9+p7pHn3PdI++Z7pXhmEYhmEYhqi5fyzfD74bfD+k3UD3cjvE94PvBt8PaTfQvdwO8f3gu8H3Q9oNdC+3Q3w/+G7w/VBJupfbIb4ffDf4fki7ge7R90z36Hume/Q90z36nukefc907wt6sixX
+        AQAAAACAAABAHwAACgEAAA==eF7t1zcSAjEQRFEW7723i4f7X5CAVtJVCqB+OEpe9qMJWkXl+4owDMMw/MGqrMm62ZBN2TLbsiPpXlf2ZN8cyKEcmWM5kXRvKmdybi7kUq7MtdxIureVO7k3D/IoT2Ypz5LuXeRV3sy7fMin+ZJvSffoe6Z79D3TPfqe6R59z3SPvme6R98z3aPvme7R90z36Hume/Q9071qGIZhGP5h7p/u+8h3ke+jtIvoXm5n+T7yXeT7KO0iupfbWb6PfBf5Pkq7iO7ldpbvI99Fvo9KSfdyO8v3ke8i30dpF9E9+p7pHn3PdI++Z7pH3zPdo++Z7tH3TPfoe6Z79D3TPfqe6R59z3TvAxSdM5E=
       </DataArray>
     </CellData>
   </Piece>

python/tests/reference/Grid/scale_grid_10-11-10.vti

@@ -3,7 +3,7 @@
   <ImageData WholeExtent="0 10 0 11 0 10" Origin="0 0 0" Spacing="8e-7 4.5454545454545457e-7 4e-7" Direction="1 0 0 0 1 0 0 0 1">
     <FieldData>
       <Array type="String" Name="comments" NumberOfTuples="1" format="binary">
-        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wc0/Nz9UrTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFYKEKs=
+        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wcy/KTNErTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFW2EKk=
       </Array>
     </FieldData>
   <Piece Extent="0 10 0 11 0 10">
@@ -11,7 +11,7 @@
     </PointData>
     <CellData>
       <DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
-        AQAAAACAAABgIgAA/gAAAA==eF7t2CcOA0EQRFHnnHPO8f4XNHAtKWmQPxhZveSxj6ZBbbn0/cphGIZh+EdWZFXWzLpsyKbZkm1J9zqyK3tmXw7k0BzJsaR7EzmVM3MuF3JpruRa5t7byK3cmXt5kEfzJM+S7l3kVd7Mu3zIp/mSb0n36Puge/R90D36Pefeo++D7tH3Qffo+6B79H3QPfo9596j74Pu0fdB9yphGIZhGIZh+IOp/5S+v3x3+f4qdhfdS+0431++u3x/FbuL7qV2nO8l30m+l4qdlHsvteN8f/nu8v11knQvteN8f/nu8v1V7C66R98H3aPvg+7R7zn3Hn0fdI++D7r3AZqNMMY=
+        AQAAAACAAABgIgAAEAEAAA==eF7t2DcOAkEQRFG8996zeLj/BQmoTUqa7Acj1Ju87EfTUmmrld9XDcMwDMM/sibrsmE2ZUu2zY7sSrrXk305MIdyJMfmRE4l3ZvJuVyYS7mSa3MjtzL33k7u5cE8ypM8m4W8SLp3lTd5Nx/yKV/mW34k3aPvg+7R90H36Pece4++D7pH3wfdo++D7tH3Qffo95x7j74PukffB92rhWEYhuEfmvpv4nvQd6DvwXIH0r3UrvQ96DvQ92C5A+lealf6fvPd5vut3G2591K70veg70Dfg4Wke6ld6XvQd6DvwXIH0j36PugefR90j37Puffo+6B79H3QPfo+6B59H3SPfs+59+j7oHv0fdC9L3s9OLk=
       </DataArray>
     </CellData>
   </Piece>

python/tests/reference/Grid/scale_grid_10-13-10.vti

@@ -3,7 +3,7 @@
   <ImageData WholeExtent="0 10 0 13 0 10" Origin="0 0 0" Spacing="8e-7 3.8461538461538463e-7 4e-7" Direction="1 0 0 0 1 0 0 0 1">
     <FieldData>
       <Array type="String" Name="comments" NumberOfTuples="1" format="binary">
-        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wc0/Nz9UrTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFYKEKs=
+        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wcy/KTNErTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFW2EKk=
       </Array>
     </FieldData>
   <Piece Extent="0 10 0 13 0 10">
@@ -11,7 +11,7 @@
     </PointData>
     <CellData>
       <DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
-        AQAAAACAAACgKAAACQEAAA==eF7t2blWAlEURFEBlcGJQcAZZBD4/x804HRSa72sghtUJzs72Q3qde/m+vVijDHGGLGPA7wV7/Aeh+IIx+juTfABH8UnfMYXcYozrN6b4wJfxSWucC2+4TtW733gJ36J3/iDG3GLv1i9t8M9HsQj/uFJPOMF3T33/bp77nur3nPfW/We+96q99z36+6579fdc99b9Z773qr33PdWvee+X3evH2OMMcYYY4wNW/8ZdP/r7tf93+1+d6/1jqB7XXe67vVup1fvtd4RdK/rTte93u306r3WO4Ludd3pute3WL3XekfQ/a+7X/d/t/vdPff9unvue6vec99b9Z773qr33Pfr7v0DeUA5pA==
+        AQAAAACAAACgKAAAFwEAAA==eF7t2DduA1EUQ1ErWcHK0co57X+DKsRpCPyOBQtOc7rbEXh/aj/frxZjjDHGCOuwAZtkC/7CNtmBXeje68E/2CcHcAhH5BhOoLo3hTM4JxdwCVfkP1xD994GbuGO3MMDPJIneIbq3gVe4Y28wwd8ki/4hu499d7ce+r9qnvqvbn31PtV99R7c++p9+beU+9X3VPvzb2n3q+6p96be68eY4wxxkiW/ivyfc53Od/n1V3u3iu9G/je5zuf7/3qzlf3Su8Gvs/5Luf7vLrL3XuldwPf+3zn871/gupe6d3A9znf5XyfV3e5e0+9N/eeer/qnnpv7j31ftU99d7ce+q9uffU+1X31Htz76n3q+6p9+be+wA2m0MJ
       </DataArray>
     </CellData>
   </Piece>

python/tests/reference/Grid/scale_grid_10-20-2.vti

@@ -3,7 +3,7 @@
   <ImageData WholeExtent="0 10 0 20 0 2" Origin="0 0 0" Spacing="8e-7 2.5000000000000004e-7 0.000002" Direction="1 0 0 0 1 0 0 0 1">
     <FieldData>
       <Array type="String" Name="comments" NumberOfTuples="1" format="binary">
-        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wc0/Nz9UrTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFYKEKs=
+        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wcy/KTNErTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFW2EKk=
       </Array>
     </FieldData>
   <Piece Extent="0 10 0 20 0 2">
@@ -11,7 +11,7 @@
     </PointData>
     <CellData>
       <DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="40">
-        AQAAAACAAACADAAAdgAAAA==eF7tzMsWgQAYReEuKJUoRKGLS3n/JzSwJ51ZRq3l35Nvtl3nm2uapmmaP+ihLy5wiSsxwBDn/ltjhLGY4AZTcYs7tN/Yqb8Mc9yLBzxiIZ7wjP/2K7HCi3jFG9Zigy3ab+zUX4d3fIhPfGEvDvjGuf8+MEUQzQ==
+        AQAAAACAAACADAAAbQAAAA==eF7tzDcCggAQBVGiooCYQcFM8P4ntGCqX9ruTvO6CYO50HVd13X/MMJYTDDFhbjEDK39VrjGXCywxI1Y4Rat/Xa4x4N4xBOexRobtPa74BVbscMb3sUHPtHa74Vv/Ig9DjiKE37R2u8H6ycQzQ==
       </DataArray>
     </CellData>
   </Piece>

python/tests/reference/Grid/scale_grid_5-4-20.vti

@@ -3,7 +3,7 @@
   <ImageData WholeExtent="0 5 0 4 0 20" Origin="0 0 0" Spacing="0.0000016 0.00000125 2e-7" Direction="1 0 0 0 1 0 0 0 1">
     <FieldData>
       <Array type="String" Name="comments" NumberOfTuples="1" format="binary">
-        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wc0/Nz9UrTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFYKEKs=
+        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wcy/KTNErTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFW2EKk=
       </Array>
     </FieldData>
   <Piece Extent="0 5 0 4 0 20">
@@ -11,7 +11,7 @@
     </PointData>
     <CellData>
       <DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
-        AQAAAACAAACADAAAeAAAAA==eF7t1TkWggAUQ1EREZS5BFQmB9j/Bin++40b8OBJmtumSoKDJZDyBx7xhGeM8YJXzLDAEmts8YYP7HHECZ/4xg+uqH4W9bOon+W7n5Ryf/oPhxih70yCvjMp5ug7U2GDHd7Rd2ZA35kZX+g7s6D6Wf613wZ8oBHR
+        AQAAAACAAACADAAAeAAAAA==eF7t1bkWglAUQ1GfiICAQMngwKTy/z9occ+t6LVJmt2myUo4WIKUcucRTxhjghlesMASK2ywwwFv+MARJ1xwxTduqH4W9bP8q5+Ucq//SIS+4zOm6DvO0Xd8xRpb7NF3fMcn+o5n9B2/8IPqZ1E/y6/7fQE+LBGB
       </DataArray>
     </CellData>
   </Piece>

python/tests/reference/Grid/scale_grid_8-10-12.vti

@@ -3,7 +3,7 @@
   <ImageData WholeExtent="0 8 0 10 0 12" Origin="0 0 0" Spacing="0.000001 5.000000000000001e-7 3.333333333333333e-7" Direction="1 0 0 0 1 0 0 0 1">
     <FieldData>
       <Array type="String" Name="comments" NumberOfTuples="1" format="binary">
-        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wc0/Nz9UrTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFYKEKs=
+        AQAAAACAAAA+AAAAQQAAAA==eF5LScxNLM7Wcy/KTNErTk7MSVUos7TUAyNdSyDQLagsSS0uUdAwMjC01DU01DUwUjA0tDK1sDIw0GQAAFW2EKk=
       </Array>
     </FieldData>
   <Piece Extent="0 8 0 10 0 12">
@@ -11,7 +11,7 @@
     </PointData>
     <CellData>
       <DataArray type="Int64" Name="material" format="binary" RangeMin="1" RangeMax="41">
-        AQAAAACAAAAAHgAAyAAAAA==eF7t2LcBAkEUxFDu8N5776H/BglQBwpvNnmRskn+FrX/K2KMsWKWWMcGNrGFbexgF23fwz4OcIgjHOMEp2j7Gc5xgUtc4Ro3uEXb73CPBzziCc94wSva/oZ3fOATX/jGD37R9nY/trf7sb3dj+3tfmxv92N7ux/b2/3Y3u7H9nY/trf7sb3dj+3tfmxv92N7ux/b2/3YvowxxhhjjJXS/gPa3t4htrd3iO3tHWJ7e4fY3u7H9nY/trf7sb3dj+3tfmz/A1V7KtE=
+        AQAAAACAAAAAHgAAyQAAAA==eF7t17cRw0AUxFCR8t5776X+G1QgdIAZRnvJi5Bt8K+o/V8RY4wxxsossY4NbGIL29jBLtq+h30c4BBHOMYJTtH2M5zjApe4wjVucIu23+EeD3jEE57xgle0/Q3v+MAnvvCNH/yi7e1+bG/3Y3u7H9vb/dje7sf2dj+2t/uxvd2P7e1+bG/3Y/syxhhjjJVr/8G2t3eg7e0daHt7B9re3oG2t/uxvd2P7e1+bG/3Y3u7H9vb/dje7sf2dj+2t/uxvd2P7X9ILyox
       </DataArray>
     </CellData>
   </Piece>

python/tests/test_Grid.py

@@ -212,6 +212,14 @@ class TestGrid:
         assert     default == modified.renumber()
 
 
+    def test_assemble(self):
+        cells = np.random.randint(8,16,3)
+        N = cells.prod()
+        g = Grid(np.arange(N).reshape(cells),np.ones(3))
+        idx = np.random.randint(0,N,N).reshape(cells)
+        assert (idx == g.assemble(idx).material).all
+
+
     def test_substitute(self,default):
         offset = np.random.randint(1,500)
         modified = Grid(default.material + offset,

python/tests/test_grid_filters.py

@@ -144,16 +144,15 @@ class TestGridFilters:
     def test_regrid_identity(self):
          size = np.random.random(3)                                                                 # noqa
          cells = np.random.randint(8,32,(3))
-         F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
+         F = np.broadcast_to(np.eye(3), (*cells,3,3))
-         assert all(grid_filters.regrid(size,F,cells) == np.arange(cells.prod()))
+         assert (grid_filters.regrid(size,F,cells).flatten() == np.arange(cells.prod())).all
 
     def test_regrid_double_cells(self):
          size = np.random.random(3)                                                                 # noqa
          cells = np.random.randint(8,32,(3))
          g = Grid.from_Voronoi_tessellation(cells,size,seeds.from_random(size,10))
-         F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
+         F = np.broadcast_to(np.eye(3), (*cells,3,3))
-         assert all(g.scale(cells*2).material.flatten() ==
+         assert g.scale(cells*2) == g.assemble(grid_filters.regrid(size,F,cells*2))
-                    g.material.flatten()[grid_filters.regrid(size,F,cells*2)])
 
     @pytest.mark.parametrize('differential_operator',[grid_filters.curl,
                                                       grid_filters.divergence,
@@ -319,7 +318,6 @@ class TestGridFilters:
      ]
 
     @pytest.mark.parametrize('field_def,div_def',div_test_data)
-
     def test_div(self,field_def,div_def):
         size = np.random.random(3)+1.0
         cells = np.random.randint(8,32,(3))
@@ -336,3 +334,31 @@ class TestGridFilters:
             div=div.reshape(tuple(cells))
 
         assert np.allclose(div,grid_filters.divergence(size,field))
+
+
+    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')
+        x,y,z = map(np.random.randint,cells)
+        assert grid_filters.ravel_index(indices)[x,y,z] == np.arange(0,np.product(cells)).reshape(cells,order='F')[x,y,z]
+
+    def test_unravel_index(self):
+        cells = np.random.randint(8,32,(3))
+        indices = np.arange(np.prod(cells)).reshape(cells,order='F')
+        x,y,z = map(np.random.randint,cells)
+        assert np.all(grid_filters.unravel_index(indices)[x,y,z] == [x,y,z])
+
+    def test_ravel_unravel_index(self):
+        cells = np.random.randint(8,32,(3))
+        indices = np.random.randint(0,np.prod(cells),cells).reshape(cells)
+        assert np.all(indices==grid_filters.ravel_index(grid_filters.unravel_index(indices)))
+
+    def test_unravel_ravel_index(self):
+        cells = np.hstack([np.random.randint(8,32,(3)),1])
+        indices = np.block([np.random.randint(0,cells[0],cells),
+                            np.random.randint(0,cells[1],cells),
+                            np.random.randint(0,cells[2],cells)])
+        assert np.all(indices==grid_filters.unravel_index(grid_filters.ravel_index(indices)))