From 9979eb58f44174e1cf3c1333f5654dd31a896115 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 12:41:03 +0200
Subject: [PATCH 01/19] consistent layout for grid data

has now always the shape ([x,y,z,...]) with x fastest. For conversion from or
to linear layout ([x*y*z,...]), e.g. storage in ASCII table, reshape needs
to have the 'F' option.
Credits to Vitesh and Fran for pointing this out.
---
 processing/post/addCurl.py         |  7 ++-
 processing/post/addDisplacement.py | 14 ++---
 processing/post/addDivergence.py   |  7 ++-
 processing/post/addGradient.py     |  7 ++-
 python/damask/_geom.py             |  7 +--
 python/damask/grid_filters.py      | 88 ++++++++++++++++--------------
 python/tests/test_grid_filters.py  | 14 ++---
 7 files changed, 75 insertions(+), 69 deletions(-)

diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py
index 87d1ab2f6..17459a2df 100755
--- a/processing/post/addCurl.py
+++ b/processing/post/addCurl.py
@@ -49,9 +49,10 @@ for name in filenames:
     for label in options.labels:
         field = table.get(label)
         shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3)                         # vector or tensor
-        field = field.reshape(np.append(grid[::-1],shape))
+        field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape)
+        curl  = damask.grid_filters.curl(size,field)
         table.add('curlFFT({})'.format(label),
-                  damask.grid_filters.curl(size[::-1],field).reshape(-1,np.prod(shape)),
+                  curl.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape),order='F'),
                   scriptID+' '+' '.join(sys.argv[1:]))
-    
+
     table.to_ASCII(sys.stdout if name is None else name)
diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py
index f74d876bc..a9424f5e1 100755
--- a/processing/post/addDisplacement.py
+++ b/processing/post/addDisplacement.py
@@ -51,23 +51,23 @@ for name in filenames:
 
     table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
     grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
-    
-    F = table.get(options.f).reshape(np.append(grid[::-1],(3,3)))
+
+    F = table.get(options.f).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3))
     if options.nodal:
-        table = damask.Table(damask.grid_filters.node_coord0(grid[::-1],size[::-1]).reshape(-1,3),
+        table = damask.Table(damask.grid_filters.node_coord0(grid,size).reshape(-1,3,order='F'),
                              {'pos':(3,)})
         table.add('avg({}).{}'.format(options.f,options.pos),
-                  damask.grid_filters.node_displacement_avg(size[::-1],F).reshape(-1,3),
+                  damask.grid_filters.node_displacement_avg(size,F).reshape(-1,3,order='F'),
                   scriptID+' '+' '.join(sys.argv[1:]))
         table.add('fluct({}).{}'.format(options.f,options.pos),
-                  damask.grid_filters.node_displacement_fluct(size[::-1],F).reshape(-1,3),
+                  damask.grid_filters.node_displacement_fluct(size,F).reshape(-1,3,order='F'),
                   scriptID+' '+' '.join(sys.argv[1:]))
         table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt')
     else:
         table.add('avg({}).{}'.format(options.f,options.pos),
-                  damask.grid_filters.cell_displacement_avg(size[::-1],F).reshape(-1,3),
+                  damask.grid_filters.cell_displacement_avg(size,F).reshape(-1,3,order='F'),
                   scriptID+' '+' '.join(sys.argv[1:]))
         table.add('fluct({}).{}'.format(options.f,options.pos),
-                  damask.grid_filters.cell_displacement_fluct(size[::-1],F).reshape(-1,3),
+                  damask.grid_filters.cell_displacement_fluct(size,F).reshape(-1,3,order='F'),
                   scriptID+' '+' '.join(sys.argv[1:]))
         table.to_ASCII(sys.stdout if name is None else name)
diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py
index 2619bc499..50048b44e 100755
--- a/processing/post/addDivergence.py
+++ b/processing/post/addDivergence.py
@@ -49,9 +49,10 @@ for name in filenames:
     for label in options.labels:
         field = table.get(label)
         shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3)                         # vector or tensor
-        field = field.reshape(np.append(grid[::-1],shape))
+        field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape)
+        div   = damask.grid_filters.divergence(size,field)
         table.add('divFFT({})'.format(label),
-                  damask.grid_filters.divergence(size[::-1],field).reshape(-1,np.prod(shape)//3),
+                  div.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)//3,order='F'),
                   scriptID+' '+' '.join(sys.argv[1:]))
-    
+
     table.to_ASCII(sys.stdout if name is None else name)
diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py
index 409b2ce6d..c6df0eacf 100755
--- a/processing/post/addGradient.py
+++ b/processing/post/addGradient.py
@@ -49,9 +49,10 @@ for name in filenames:
     for label in options.labels:
         field = table.get(label)
         shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,)                          # scalar or vector
-        field = field.reshape(np.append(grid[::-1],shape))
+        field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape)
+        grad  = damask.grid_filters.gradient(size,field)
         table.add('gradFFT({})'.format(label),
-                  damask.grid_filters.gradient(size[::-1],field).reshape(-1,np.prod(shape)*3),
+                  grad.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)*3,order='F'),
                   scriptID+' '+' '.join(sys.argv[1:]))
-    
+
     table.to_ASCII(sys.stdout if name is None else name)
diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index e2c2428fe..a33e855d8 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -357,7 +357,6 @@ class Geom:
             seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
             seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
             coords  = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F')
-
         else:
             weights_p = weights.flatten()
             seeds_p   = seeds
@@ -370,10 +369,10 @@ class Geom:
         microstructure = np.array(result.get())
 
         if periodic:
-            microstructure = microstructure.reshape(grid*3)
+            microstructure = microstructure.reshape(grid*3,order='F')
             microstructure = microstructure[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
         else:
-            microstructure = microstructure.reshape(grid)
+            microstructure = microstructure.reshape(grid,order='F')
 
         #comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version)
         return Geom(microstructure+1,size,homogenization=1)
@@ -401,7 +400,7 @@ class Geom:
         devNull,microstructure = KDTree.query(coords)
 
         #comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version)
-        return Geom(microstructure.reshape(grid)+1,size,homogenization=1)
+        return Geom(microstructure.reshape(grid,order='F')+1,size,homogenization=1)
 
 
     def to_file(self,fname,pack=None):
diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py
index d8b136a6b..0c1c3867a 100644
--- a/python/damask/grid_filters.py
+++ b/python/damask/grid_filters.py
@@ -7,7 +7,7 @@ def _ks(size,grid,first_order=False):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
 
     """
@@ -19,8 +19,7 @@ def _ks(size,grid,first_order=False):
 
     k_si = _np.arange(grid[2]//2+1)/size[2]
 
-    kk, kj, ki = _np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
-    return _np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3)
+    return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
 
 
 def curl(size,field):
@@ -29,7 +28,7 @@ def curl(size,field):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
 
     """
@@ -53,7 +52,7 @@ def divergence(size,field):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
 
     """
@@ -73,7 +72,7 @@ def gradient(size,field):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
 
     """
@@ -93,9 +92,9 @@ def cell_coord0(grid,size,origin=_np.zeros(3)):
 
     Parameters
     ----------
-    grid : numpy.ndarray
+    grid : numpy.ndarray of shape (3)
         number of grid points.
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
     origin : numpy.ndarray, optional
         physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
@@ -103,7 +102,11 @@ def cell_coord0(grid,size,origin=_np.zeros(3)):
     """
     start = origin        + size/grid*.5
     end   = origin + size - size/grid*.5
-    return _np.mgrid[start[0]:end[0]:grid[0]*1j,start[1]:end[1]:grid[1]*1j,start[2]:end[2]:grid[2]*1j].T
+
+    return _np.stack(_np.meshgrid(_np.linspace(start[0],end[0],grid[0]),
+                                  _np.linspace(start[1],end[1],grid[1]),
+                                  _np.linspace(start[2],end[2],grid[2]),indexing = 'ij'),
+                     axis = -1)
 
 
 def cell_displacement_fluct(size,F):
@@ -112,7 +115,7 @@ def cell_displacement_fluct(size,F):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
     F : numpy.ndarray
         deformation gradient field.
@@ -139,14 +142,14 @@ def cell_displacement_avg(size,F):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
     F : numpy.ndarray
         deformation gradient field.
 
     """
     F_avg = _np.average(F,axis=(0,1,2))
-    return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),cell_coord0(F.shape[:3][::-1],size))
+    return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),cell_coord0(F.shape[:3],size))
 
 
 def cell_displacement(size,F):
@@ -155,7 +158,7 @@ def cell_displacement(size,F):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
     F : numpy.ndarray
         deformation gradient field.
@@ -170,30 +173,30 @@ def cell_coord(size,F,origin=_np.zeros(3)):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
     F : numpy.ndarray
         deformation gradient field.
-    origin : numpy.ndarray, optional
+    origin : numpy.ndarray of shape (3), optional
         physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
 
     """
-    return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F)
+    return cell_coord0(F.shape[:3],size,origin) + cell_displacement(size,F)
 
 
 def cell_coord0_gridSizeOrigin(coord0,ordered=True):
     """
-    Return grid 'DNA', i.e. grid, size, and origin from array of cell positions.
+    Return grid 'DNA', i.e. grid, size, and origin from 1D array of cell positions.
 
     Parameters
     ----------
-    coord0 : numpy.ndarray
-        array of undeformed cell coordinates.
+    coord0 : numpy.ndarray of shape (:,3)
+        undeformed cell coordinates.
     ordered : bool, optional
         expect coord0 data to be ordered (x fast, z slow).
 
     """
-    coords    = [_np.unique(coord0[:,i]) for i in range(3)]
+    coords    = [_np.unique(coord0[:,i]) for i in range(3)] # _np.unique(coord0, axis=1)
     mincorner = _np.array(list(map(min,coords)))
     maxcorner = _np.array(list(map(max,coords)))
     grid      = _np.array(list(map(len,coords)),'i')
@@ -216,7 +219,7 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
            _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2])):
         raise ValueError('Regular grid spacing violated.')
 
-    if ordered and not _np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin)):
+    if ordered and not _np.allclose(coord0.reshape(tuple(grid)+(3,),order='F'),cell_coord0(grid,size,origin)):
         raise ValueError('Input data is not a regular grid.')
 
     return (grid,size,origin)
@@ -241,17 +244,18 @@ def node_coord0(grid,size,origin=_np.zeros(3)):
 
     Parameters
     ----------
-    grid : numpy.ndarray
+    grid : numpy.ndarray of shape (3)
         number of grid points.
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
-    origin : numpy.ndarray, optional
+    origin : numpy.ndarray of shape (3), optional
         physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
 
     """
-    return _np.mgrid[origin[0]:size[0]+origin[0]:(grid[0]+1)*1j,
-                     origin[1]:size[1]+origin[1]:(grid[1]+1)*1j,
-                     origin[2]:size[2]+origin[2]:(grid[2]+1)*1j].T
+    return _np.stack(_np.meshgrid(_np.linspace(origin[0],size[0]+origin[0],grid[0]+1),
+                                  _np.linspace(origin[1],size[1]+origin[1],grid[1]+1),
+                                  _np.linspace(origin[2],size[2]+origin[2],grid[2]+1),indexing = 'ij'),
+                     axis = -1)
 
 
 def node_displacement_fluct(size,F):
@@ -260,7 +264,7 @@ def node_displacement_fluct(size,F):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
     F : numpy.ndarray
         deformation gradient field.
@@ -275,14 +279,14 @@ def node_displacement_avg(size,F):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
     F : numpy.ndarray
         deformation gradient field.
 
     """
     F_avg = _np.average(F,axis=(0,1,2))
-    return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),node_coord0(F.shape[:3][::-1],size))
+    return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),node_coord0(F.shape[:3],size))
 
 
 def node_displacement(size,F):
@@ -291,7 +295,7 @@ def node_displacement(size,F):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
     F : numpy.ndarray
         deformation gradient field.
@@ -306,15 +310,15 @@ def node_coord(size,F,origin=_np.zeros(3)):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size of the periodic field.
     F : numpy.ndarray
         deformation gradient field.
-    origin : numpy.ndarray, optional
+    origin : numpy.ndarray of shape (3), optional
         physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
 
     """
-    return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F)
+    return node_coord0(F.shape[:3],size,origin) + node_displacement(size,F)
 
 
 def cell_2_node(cell_data):
@@ -335,19 +339,19 @@ def node_2_cell(node_data):
     return c[:-1,:-1,:-1]
 
 
-def node_coord0_gridSizeOrigin(coord0,ordered=False):
+def node_coord0_gridSizeOrigin(coord0,ordered=True):
     """
-    Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions.
+    Return grid 'DNA', i.e. grid, size, and origin from 1D array of nodal positions.
 
     Parameters
     ----------
-    coord0 : numpy.ndarray
-        array of undeformed nodal coordinates.
+    coord0 : numpy.ndarray of shape (:,3)
+        undeformed nodal coordinates.
     ordered : bool, optional
         expect coord0 data to be ordered (x fast, z slow).
 
     """
-    coords    = [_np.unique(coord0[:,i]) for i in range(3)]
+    coords    = [_np.unique(coord0[:,i]) for i in range(3)] # _np.unique(coord0, axis=1)
     mincorner = _np.array(list(map(min,coords)))
     maxcorner = _np.array(list(map(max,coords)))
     grid      = _np.array(list(map(len,coords)),'i') - 1
@@ -362,7 +366,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False):
            _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1)):
         raise ValueError('Regular grid spacing violated.')
 
-    if ordered and not _np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)):
+    if ordered and not _np.allclose(coord0.reshape(tuple(grid+1)+(3,),order='F'),node_coord0(grid,size,origin)):
         raise ValueError('Input data is not a regular grid.')
 
     return (grid,size,origin)
@@ -374,7 +378,7 @@ def regrid(size,F,new_grid):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    size : numpy.ndarray of shape (3)
         physical size
     F : numpy.ndarray
         deformation gradient field
@@ -382,7 +386,7 @@ def regrid(size,F,new_grid):
         new grid for undeformed coordinates
 
     """
-    c = cell_coord0(F.shape[:3][::-1],size) \
+    c = cell_coord0(F.shape[:3],size) \
       + cell_displacement_avg(size,F) \
       + cell_displacement_fluct(size,F)
 
diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py
index acbdbf688..3152072b7 100644
--- a/python/tests/test_grid_filters.py
+++ b/python/tests/test_grid_filters.py
@@ -4,18 +4,18 @@ import numpy as np
 from damask import grid_filters
 
 class TestGridFilters:
-   
+
     def test_cell_coord0(self):
          size = np.random.random(3)
          grid = np.random.randint(8,32,(3))
          coord = grid_filters.cell_coord0(grid,size)
-         assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid[::-1]) + (3,)
+         assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid) + (3,)
 
     def test_node_coord0(self):
          size = np.random.random(3)
          grid = np.random.randint(8,32,(3))
          coord = grid_filters.node_coord0(grid,size)
-         assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid[::-1]+1) + (3,)
+         assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid+1) + (3,)
 
     def test_coord0(self):
          size = np.random.random(3)
@@ -31,7 +31,7 @@ class TestGridFilters:
          size   = np.random.random(3)
          origin = np.random.random(3)
          coord0 = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode))                     # noqa
-         _grid,_size,_origin = eval('grid_filters.{}_coord0_gridSizeOrigin(coord0.reshape(-1,3))'.format(mode))
+         _grid,_size,_origin = eval('grid_filters.{}_coord0_gridSizeOrigin(coord0.reshape(-1,3,order="F"))'.format(mode))
          assert np.allclose(grid,_grid) and np.allclose(size,_size) and np.allclose(origin,_origin)
 
     def test_displacement_fluct_equivalence(self):
@@ -57,9 +57,9 @@ class TestGridFilters:
          shifted   = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode))
          unshifted = eval('grid_filters.{}_coord0(grid,size)'.format(mode))
          if   mode == 'cell':
-            assert  np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1])  +(3,)))
+            assert  np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid)  +(3,)))
          elif mode == 'node':
-            assert  np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]+1)+(3,)))
+            assert  np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid+1)+(3,)))
 
     @pytest.mark.parametrize('function',[grid_filters.cell_displacement_avg,
                                          grid_filters.node_displacement_avg])
@@ -83,5 +83,5 @@ class TestGridFilters:
     def test_regrid(self):
          size = np.random.random(3)
          grid = np.random.randint(8,32,(3))
-         F    = np.broadcast_to(np.eye(3), tuple(grid[::-1])+(3,3))
+         F    = np.broadcast_to(np.eye(3), tuple(grid)+(3,3))
          assert all(grid_filters.regrid(size,F,grid) == np.arange(grid.prod()))

From 02dde3c2556a4e8c7ccf83a4b21670e9cc713083 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 13:10:13 +0200
Subject: [PATCH 02/19] test also invalid operations

---
 python/tests/test_grid_filters.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py
index 3152072b7..9eba75298 100644
--- a/python/tests/test_grid_filters.py
+++ b/python/tests/test_grid_filters.py
@@ -80,6 +80,11 @@ class TestGridFilters:
          F    = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3))
          assert np.allclose(function(size,F),0.0)
 
+    def test_invalid_coordinates(self):
+        invalid_coordinates = np.random.random((np.random.randint(12,52),3))
+        with pytest.raises(ValueError):
+            grid_filters.coord0_check(invalid_coordinates)
+
     def test_regrid(self):
          size = np.random.random(3)
          grid = np.random.randint(8,32,(3))

From f519e62cd5713d2ba5b80ec574e695ec358969cf Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 13:56:33 +0200
Subject: [PATCH 03/19] testing data layout checks

---
 python/damask/grid_filters.py     | 20 +++++++++---------
 python/tests/test_grid_filters.py | 34 +++++++++++++++++++++++++++++--
 2 files changed, 42 insertions(+), 12 deletions(-)

diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py
index 0c1c3867a..2ee99d276 100644
--- a/python/damask/grid_filters.py
+++ b/python/damask/grid_filters.py
@@ -196,7 +196,7 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
         expect coord0 data to be ordered (x fast, z slow).
 
     """
-    coords    = [_np.unique(coord0[:,i]) for i in range(3)] # _np.unique(coord0, axis=1)
+    coords    = [_np.unique(coord0[:,i]) for i in range(3)]
     mincorner = _np.array(list(map(min,coords)))
     maxcorner = _np.array(list(map(max,coords)))
     grid      = _np.array(list(map(len,coords)),'i')
@@ -214,13 +214,13 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
     start = origin + delta*.5
     end   = origin - delta*.5 + size
 
-    if not _np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0])) and \
-           _np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1])) and \
-           _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2])):
+    if not (_np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0])) and \
+            _np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1])) and \
+            _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2]))):
         raise ValueError('Regular grid spacing violated.')
 
     if ordered and not _np.allclose(coord0.reshape(tuple(grid)+(3,),order='F'),cell_coord0(grid,size,origin)):
-        raise ValueError('Input data is not a regular grid.')
+        raise ValueError('Input data is not ordered (x fast, z slow).')
 
     return (grid,size,origin)
 
@@ -351,7 +351,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=True):
         expect coord0 data to be ordered (x fast, z slow).
 
     """
-    coords    = [_np.unique(coord0[:,i]) for i in range(3)] # _np.unique(coord0, axis=1)
+    coords    = [_np.unique(coord0[:,i]) for i in range(3)]
     mincorner = _np.array(list(map(min,coords)))
     maxcorner = _np.array(list(map(max,coords)))
     grid      = _np.array(list(map(len,coords)),'i') - 1
@@ -361,13 +361,13 @@ def node_coord0_gridSizeOrigin(coord0,ordered=True):
     if (grid+1).prod() != len(coord0):
         raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid))
 
-    if not _np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \
-           _np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \
-           _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1)):
+    if not (_np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \
+            _np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \
+            _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1))):
         raise ValueError('Regular grid spacing violated.')
 
     if ordered and not _np.allclose(coord0.reshape(tuple(grid+1)+(3,),order='F'),node_coord0(grid,size,origin)):
-        raise ValueError('Input data is not a regular grid.')
+        raise ValueError('Input data is not ordered (x fast, z slow).')
 
     return (grid,size,origin)
 
diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py
index 9eba75298..eb359006a 100644
--- a/python/tests/test_grid_filters.py
+++ b/python/tests/test_grid_filters.py
@@ -80,10 +80,40 @@ class TestGridFilters:
          F    = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3))
          assert np.allclose(function(size,F),0.0)
 
-    def test_invalid_coordinates(self):
+    @pytest.mark.parametrize('function',[grid_filters.coord0_check,
+                                         grid_filters.node_coord0_gridSizeOrigin,
+                                         grid_filters.cell_coord0_gridSizeOrigin])
+    def test_invalid_coordinates(self,function):
         invalid_coordinates = np.random.random((np.random.randint(12,52),3))
         with pytest.raises(ValueError):
-            grid_filters.coord0_check(invalid_coordinates)
+            function(invalid_coordinates)
+
+    @pytest.mark.parametrize('function',[grid_filters.node_coord0_gridSizeOrigin,
+                                         grid_filters.cell_coord0_gridSizeOrigin])
+    def test_uneven_spaced_coordinates(self,function):
+        start = np.random.random(3)
+        end   = np.random.random(3)*10. + start
+        grid  = np.random.randint(8,32,(3))
+        uneven = np.stack(np.meshgrid(np.logspace(start[0],end[0],grid[0]),
+                                      np.logspace(start[1],end[1],grid[1]),
+                                      np.logspace(start[2],end[2],grid[2]),indexing = 'ij'),
+                           axis = -1).reshape((grid.prod(),3),order='F')
+        with pytest.raises(ValueError):
+            function(uneven)
+
+    @pytest.mark.parametrize('mode',[True,False])
+    @pytest.mark.parametrize('function',[grid_filters.node_coord0_gridSizeOrigin,
+                                         grid_filters.cell_coord0_gridSizeOrigin])
+    def test_unordered_coordinates(self,function,mode):
+        origin = np.random.random(3)
+        size   = np.random.random(3)*10.+origin
+        grid  = np.random.randint(8,32,(3))
+        unordered = grid_filters.node_coord0(grid,size,origin).reshape(-1,3)
+        if mode:
+            with pytest.raises(ValueError):
+                function(unordered,mode)
+        else:
+            function(unordered,mode)
 
     def test_regrid(self):
          size = np.random.random(3)

From a367128436fb099c3fe4f8b920f61b31e169f1f4 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 13:57:33 +0200
Subject: [PATCH 04/19] not used anymore

---
 processing/post/addDisplacement.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py
index a9424f5e1..6f4a60192 100755
--- a/processing/post/addDisplacement.py
+++ b/processing/post/addDisplacement.py
@@ -5,8 +5,6 @@ import sys
 from io import StringIO
 from optparse import OptionParser
 
-import numpy as np
-
 import damask
 
 

From 8465c3ab1b7b19a84b4aa00c56b7aada29131ea1 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 19:13:45 +0200
Subject: [PATCH 05/19] adjusting to new data layout

---
 processing/pre/geom_fromVoronoiTessellation.py | 2 +-
 processing/pre/seeds_fromGeom.py               | 8 ++++----
 processing/pre/seeds_fromRandom.py             | 2 +-
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py
index df40176d8..2b4da133b 100755
--- a/processing/pre/geom_fromVoronoiTessellation.py
+++ b/processing/pre/geom_fromVoronoiTessellation.py
@@ -224,7 +224,7 @@ for name in filenames:
 
     header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
            + config_header
-    geom = damask.Geom(indices.reshape(grid),size,origin,
+    geom = damask.Geom(indices.reshape(grid,order='F'),size,origin,
                        homogenization=options.homogenization,comments=header)
     damask.util.croak(geom)
 
diff --git a/processing/pre/seeds_fromGeom.py b/processing/pre/seeds_fromGeom.py
index 80dc0d6f5..0b741a077 100755
--- a/processing/pre/seeds_fromGeom.py
+++ b/processing/pre/seeds_fromGeom.py
@@ -45,7 +45,7 @@ options.blacklist = [int(i) for i in options.blacklist]
 
 for name in filenames:
     damask.util.report(scriptName,name)
-    
+
     geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
     microstructure = geom.get_microstructure().reshape((-1,1),order='F')
 
@@ -53,9 +53,9 @@ for name in filenames:
                           np.full(geom.grid.prod(),True,dtype=bool),
                           np.in1d(microstructure,options.blacklist,invert=True)  if options.blacklist else \
                           np.full(geom.grid.prod(),True,dtype=bool))
-   
-    seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3)
-    
+
+    seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
+
     comments = geom.comments \
              + [scriptID + ' ' + ' '.join(sys.argv[1:]),
                 'grid\ta {}\tb {}\tc {}'.format(*geom.grid),
diff --git a/processing/pre/seeds_fromRandom.py b/processing/pre/seeds_fromRandom.py
index 2de513c2c..be690713c 100755
--- a/processing/pre/seeds_fromRandom.py
+++ b/processing/pre/seeds_fromRandom.py
@@ -128,7 +128,7 @@ for name in filenames:
 
 
     if not options.selective:
-        coords = damask.grid_filters.cell_coord0(grid,size).reshape(-1,3)
+        coords = damask.grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
         seeds = coords[np.random.choice(np.prod(grid), options.N, replace=False)] \
               + np.broadcast_to(size/grid,(options.N,3))*(np.random.rand(options.N,3)*.5-.25)       # wobble without leaving grid
     else:

From 4db91fff03ab68f68cf26ac2c5311fd10c82b8c8 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 20:14:27 +0200
Subject: [PATCH 06/19] simplified

---
 python/damask/_geom.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index a33e855d8..f01277f95 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -322,11 +322,10 @@ class Geom:
         if i != grid.prod():
             raise TypeError('Invalid file: expected {} entries, found {}'.format(grid.prod(),i))
 
-        microstructure = microstructure.reshape(grid,order='F')
-        if not np.any(np.mod(microstructure.flatten(),1) != 0.0):                                   # no float present
+        if not np.any(np.mod(microstructure,1) != 0.0):                                             # no float present
             microstructure = microstructure.astype('int')
 
-        return Geom(microstructure.reshape(grid),size,origin,homogenization,comments)
+        return Geom(microstructure.reshape(grid,order='F'),size,origin,homogenization,comments)
 
 
     @staticmethod

From 1ccda22ae354de44c78da730543ccd7b747cb906 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 20:16:25 +0200
Subject: [PATCH 07/19] nothing to flatten here

---
 python/damask/_geom.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index f01277f95..8010cd167 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -351,13 +351,13 @@ class Geom:
 
         """
         if periodic:
-            weights_p = np.tile(weights,27).flatten(order='F')                                      # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
+            weights_p = np.tile(weights,27)                                                         # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
             seeds_p = np.vstack((seeds  -np.array([size[0],0.,0.]),seeds,  seeds  +np.array([size[0],0.,0.])))
             seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
             seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
             coords  = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F')
         else:
-            weights_p = weights.flatten()
+            weights_p = weights
             seeds_p   = seeds
             coords    = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
 

From 25d29bf438d00390f105bbdff29c6c85ddc3472d Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 20:24:55 +0200
Subject: [PATCH 08/19] no need for Fortran order forward/backward

---
 python/damask/_geom.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 8010cd167..892000b7c 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -355,11 +355,11 @@ class Geom:
             seeds_p = np.vstack((seeds  -np.array([size[0],0.,0.]),seeds,  seeds  +np.array([size[0],0.,0.])))
             seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
             seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
-            coords  = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F')
+            coords  = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3)
         else:
             weights_p = weights
             seeds_p   = seeds
-            coords    = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
+            coords    = grid_filters.cell_coord0(grid,size).reshape(-1,3)
 
         pool = multiprocessing.Pool(processes = int(Environment().options['DAMASK_NUM_THREADS']))
         result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
@@ -368,10 +368,10 @@ class Geom:
         microstructure = np.array(result.get())
 
         if periodic:
-            microstructure = microstructure.reshape(grid*3,order='F')
+            microstructure = microstructure.reshape(grid*3)
             microstructure = microstructure[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
         else:
-            microstructure = microstructure.reshape(grid,order='F')
+            microstructure = microstructure.reshape(grid)
 
         #comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version)
         return Geom(microstructure+1,size,homogenization=1)
@@ -394,12 +394,12 @@ class Geom:
             perform a periodic tessellation. Defaults to True.
 
         """
-        coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
+        coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
         KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
         devNull,microstructure = KDTree.query(coords)
 
         #comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version)
-        return Geom(microstructure.reshape(grid,order='F')+1,size,homogenization=1)
+        return Geom(microstructure.reshape(grid)+1,size,homogenization=1)
 
 
     def to_file(self,fname,pack=None):

From a8b75a23c5f746fec18fd771db2c02dbe30fb320 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 20:40:47 +0200
Subject: [PATCH 09/19] simplified

---
 processing/pre/geom_fromVoronoiTessellation.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py
index 2b4da133b..e6755c28e 100755
--- a/processing/pre/geom_fromVoronoiTessellation.py
+++ b/processing/pre/geom_fromVoronoiTessellation.py
@@ -24,22 +24,22 @@ def findClosestSeed(seeds, weights, point):
 def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), periodic = True, cpus = 2):
 
     if periodic:
-        weights_p = np.tile(weights,27).flatten(order='F')                                          # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
+        weights_p = np.tile(weights.squeeze(),27)                                                   # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
         seeds_p = np.vstack((seeds  -np.array([size[0],0.,0.]),seeds,  seeds  +np.array([size[0],0.,0.])))
         seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
         seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
-        coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3,order='F')
+        coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3)
     else:
-        weights_p = weights.flatten()
+        weights_p = weights.squeeze()
         seeds_p   = seeds
-        coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
+        coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3)
 
     if cpus > 1:
         pool = multiprocessing.Pool(processes = cpus)
         result = pool.map_async(partial(findClosestSeed,seeds_p,weights_p), [coord for coord in coords])
         pool.close()
         pool.join()
-        closest_seed = np.array(result.get())
+        closest_seed = np.array(result.get()).reshape(-1,3)
     else:
         closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords])
 
@@ -52,7 +52,7 @@ def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), peri
 
 def Voronoi_tessellation(grid, size, seeds, origin = np.zeros(3), periodic = True):
 
-    coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
+    coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3)
     KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
     devNull,closest_seed = KDTree.query(coords)
 
@@ -224,7 +224,7 @@ for name in filenames:
 
     header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
            + config_header
-    geom = damask.Geom(indices.reshape(grid,order='F'),size,origin,
+    geom = damask.Geom(indices.reshape(grid),size,origin,
                        homogenization=options.homogenization,comments=header)
     damask.util.croak(geom)
 

From 9a54c326e23a5e876cd18febb154452591f7904a Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 21:28:23 +0200
Subject: [PATCH 10/19] adopted to new data layout

---
 processing/post/averageDown.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/processing/post/averageDown.py b/processing/post/averageDown.py
index 39d925d19..341cc748d 100755
--- a/processing/post/averageDown.py
+++ b/processing/post/averageDown.py
@@ -61,7 +61,7 @@ if any(shift != 0): prefix += 'shift{:+}{:+}{:+}_'.format(*shift)
 
 for name in filenames:
   damask.util.report(scriptName,name)
-  
+
   table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
 
   if (options.grid is None or options.size is None):
@@ -87,11 +87,11 @@ for name in filenames:
                   origin = list(-(packing//2)) + [0])\
                   [::packing[0],::packing[1],::packing[2],:].reshape((packedGrid.prod(),-1),order = 'F')
 
-  
+
   table = damask.Table(averagedDown,table.shapes,table.comments)
 
   coords = damask.grid_filters.cell_coord0(packedGrid,size,shift/packedGrid*size+origin)
-  table.set(options.pos, coords.reshape(-1,3))
+  table.set(options.pos, coords.reshape(-1,3,order='F'))
 
 
   outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name))

From 4278ba32ae2ddd71a4522bfe390e0f7efec518a8 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 21:44:38 +0200
Subject: [PATCH 11/19] adopting to standard data layout

---
 processing/post/addCompatibilityMismatch.py | 55 ++++++++++-----------
 1 file changed, 27 insertions(+), 28 deletions(-)

diff --git a/processing/post/addCompatibilityMismatch.py b/processing/post/addCompatibilityMismatch.py
index 0678e78ba..7cab9c2ca 100755
--- a/processing/post/addCompatibilityMismatch.py
+++ b/processing/post/addCompatibilityMismatch.py
@@ -16,8 +16,8 @@ scriptID   = ' '.join([scriptName,damask.version])
 def volTetrahedron(coords):
   """
   Return the volume of the tetrahedron with given vertices or sides.
-  
-  Ifvertices are given they must be in a NumPy array with shape (4,3): the
+
+  If vertices are given they must be in a NumPy array with shape (4,3): the
   position vectors of the 4 vertices in 3 dimensions; if the six sides are
   given, they must be an array of length 6. If both are given, the sides
   will be used in the calculation.
@@ -62,19 +62,18 @@ def volTetrahedron(coords):
 def volumeMismatch(size,F,nodes):
   """
   Calculates the volume mismatch.
-  
-  volume mismatch is defined as the difference between volume of reconstructed 
+
+  volume mismatch is defined as the difference between volume of reconstructed
   (compatible) cube and determinant of deformation gradient at Fourier point.
   """
   coords = np.empty([8,3])
-  vMismatch = np.empty(grid[::-1])
-  volInitial = size.prod()/grid.prod()
- 
+  vMismatch = np.empty(F.shape[:3])
+
 #--------------------------------------------------------------------------------------------------
 # calculate actual volume and volume resulting from deformation gradient
-  for k in range(grid[2]):
+  for k in range(grid[0]):
     for j in range(grid[1]):
-      for i in range(grid[0]):
+      for i in range(grid[2]):
         coords[0,0:3] = nodes[k,  j,  i  ,0:3]
         coords[1,0:3] = nodes[k  ,j,  i+1,0:3]
         coords[2,0:3] = nodes[k  ,j+1,i+1,0:3]
@@ -91,21 +90,21 @@ def volumeMismatch(size,F,nodes):
          + abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[5,0:3]])) \
          + abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[0,0:3]]))) \
         /np.linalg.det(F[k,j,i,0:3,0:3])
-  return vMismatch/volInitial
+  return vMismatch/(size.prod()/grid.prod())
 
 
 
 def shapeMismatch(size,F,nodes,centres):
   """
   Routine to calculate the shape mismatch.
-  
+
   shape mismatch is defined as difference between the vectors from the central point to
   the corners of reconstructed (combatible) volume element and the vectors calculated by deforming
   the initial volume element with the  current deformation gradient.
   """
   coordsInitial = np.empty([8,3])
-  sMismatch    = np.empty(grid[::-1])
-   
+  sMismatch     = np.empty(F.shape[:3])
+
 #--------------------------------------------------------------------------------------------------
 # initial positions
   coordsInitial[0,0:3] = [-size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
@@ -117,21 +116,21 @@ def shapeMismatch(size,F,nodes,centres):
   coordsInitial[6,0:3] = [+size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
   coordsInitial[7,0:3] = [-size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
   coordsInitial = coordsInitial/2.0
- 
+
 #--------------------------------------------------------------------------------------------------
 # compare deformed original and deformed positions to actual positions
-  for k in range(grid[2]):
+  for k in range(grid[0]):
     for j in range(grid[1]):
-      for i in range(grid[0]):
+      for i in range(grid[2]):
        sMismatch[k,j,i] = \
          + np.linalg.norm(nodes[k,  j,  i  ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[0,0:3]))\
-         + np.linalg.norm(nodes[k,  j,  i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\
-         + np.linalg.norm(nodes[k,  j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\
+         + np.linalg.norm(nodes[k+1,j,  i  ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\
+         + np.linalg.norm(nodes[k+1,j+1,i  ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\
          + np.linalg.norm(nodes[k,  j+1,i  ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[3,0:3]))\
-         + np.linalg.norm(nodes[k+1,j,  i  ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\
+         + np.linalg.norm(nodes[k,  j,  i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\
          + np.linalg.norm(nodes[k+1,j,  i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[5,0:3]))\
          + np.linalg.norm(nodes[k+1,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[6,0:3]))\
-         + np.linalg.norm(nodes[k+1,j+1,i  ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3]))
+         + np.linalg.norm(nodes[k  ,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3]))
   return sMismatch
 
 
@@ -174,24 +173,24 @@ if filenames == []: filenames = [None]
 
 for name in filenames:
   damask.util.report(scriptName,name)
-  
+
   table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
   grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
 
-  F = table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3)
+  F = table.get(options.defgrad).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3))
   nodes = damask.grid_filters.node_coord(size,F)
-  
+
   if options.shape:
     centers = damask.grid_filters.cell_coord(size,F)
-    shapeMismatch = shapeMismatch(  size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes,centers)
+    shapeMismatch = shapeMismatch(size,F,nodes,centers)
     table.add('shapeMismatch(({}))'.format(options.defgrad),
-              shapeMismatch.reshape(-1,1),
+              shapeMismatch.reshape(-1,1,order='F'),
               scriptID+' '+' '.join(sys.argv[1:]))
-    
+
   if options.volume:
-    volumeMismatch = volumeMismatch(size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes)
+    volumeMismatch = volumeMismatch(size,F,nodes)
     table.add('volMismatch(({}))'.format(options.defgrad),
-              volumeMismatch.reshape(-1,1),
+              volumeMismatch.reshape(-1,1,order='F'),
               scriptID+' '+' '.join(sys.argv[1:]))
 
   table.to_ASCII(sys.stdout if name is None else name)

From 85f831d5b44594db2c6460e693715a0d92d77164 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 22:14:57 +0200
Subject: [PATCH 12/19] simplified

---
 processing/post/addCompatibilityMismatch.py | 20 +++++++++-----------
 1 file changed, 9 insertions(+), 11 deletions(-)

diff --git a/processing/post/addCompatibilityMismatch.py b/processing/post/addCompatibilityMismatch.py
index 7cab9c2ca..59eafc077 100755
--- a/processing/post/addCompatibilityMismatch.py
+++ b/processing/post/addCompatibilityMismatch.py
@@ -93,7 +93,6 @@ def volumeMismatch(size,F,nodes):
   return vMismatch/(size.prod()/grid.prod())
 
 
-
 def shapeMismatch(size,F,nodes,centres):
   """
   Routine to calculate the shape mismatch.
@@ -102,20 +101,19 @@ def shapeMismatch(size,F,nodes,centres):
   the corners of reconstructed (combatible) volume element and the vectors calculated by deforming
   the initial volume element with the  current deformation gradient.
   """
-  coordsInitial = np.empty([8,3])
   sMismatch     = np.empty(F.shape[:3])
 
 #--------------------------------------------------------------------------------------------------
 # initial positions
-  coordsInitial[0,0:3] = [-size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
-  coordsInitial[1,0:3] = [+size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
-  coordsInitial[2,0:3] = [+size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
-  coordsInitial[3,0:3] = [-size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
-  coordsInitial[4,0:3] = [-size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
-  coordsInitial[5,0:3] = [+size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
-  coordsInitial[6,0:3] = [+size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
-  coordsInitial[7,0:3] = [-size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
-  coordsInitial = coordsInitial/2.0
+  delta = size/grid*.5
+  coordsInitial = np.vstack((delta * np.array((-1,-1,-1)),
+                             delta * np.array((+1,-1,-1)),
+                             delta * np.array((+1,+1,-1)),
+                             delta * np.array((-1,+1,-1)),
+                             delta * np.array((-1,-1,+1)),
+                             delta * np.array((+1,-1,+1)),
+                             delta * np.array((+1,+1,+1)),
+                             delta * np.array((-1,+1,+1))))
 
 #--------------------------------------------------------------------------------------------------
 # compare deformed original and deformed positions to actual positions

From 93db1b37cacc4ccf6b0ac3d3fbaa7660f19d76c5 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 22:26:29 +0200
Subject: [PATCH 13/19] some useful information

---
 python/damask/grid_filters.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py
index 2ee99d276..3925c8db3 100644
--- a/python/damask/grid_filters.py
+++ b/python/damask/grid_filters.py
@@ -1,3 +1,17 @@
+"""
+Filters for operations on regular grids.
+
+Notes
+-----
+The grids are defined as (x,y,z,...) where x is fastest and z is slowest.
+This convention is consistent with the geom file format.
+When converting to/from a plain list (e.g. storage in ASCII table),
+the following operations are required for tensorial:
+
+D3 = D1.reshape(grid+(-1,),order='F').reshape(grid+(3,3))
+D1 = D3.reshape(grid+(-1,)).reshape(-1,9,order='F')
+
+"""
 from scipy import spatial as _spatial
 import numpy as _np
 

From 00cb17ccfaf867a6f7e4f321a7153f89eec34098 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 22:47:19 +0200
Subject: [PATCH 14/19] commenting

---
 python/damask/grid_filters.py | 18 ++++++++++++++----
 1 file changed, 14 insertions(+), 4 deletions(-)

diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py
index 3925c8db3..f74bc92ac 100644
--- a/python/damask/grid_filters.py
+++ b/python/damask/grid_filters.py
@@ -6,7 +6,7 @@ Notes
 The grids are defined as (x,y,z,...) where x is fastest and z is slowest.
 This convention is consistent with the geom file format.
 When converting to/from a plain list (e.g. storage in ASCII table),
-the following operations are required for tensorial:
+the following operations are required for tensorial data:
 
 D3 = D1.reshape(grid+(-1,),order='F').reshape(grid+(3,3))
 D1 = D3.reshape(grid+(-1,)).reshape(-1,9,order='F')
@@ -23,6 +23,10 @@ def _ks(size,grid,first_order=False):
     ----------
     size : numpy.ndarray of shape (3)
         physical size of the periodic field.
+    grid : numpy.ndarray of shape (3)
+        number of grid points.
+    first_order : bool, optional
+        correction for first order derivatives, defaults to False.
 
     """
     k_sk = _np.where(_np.arange(grid[0])>grid[0]//2,_np.arange(grid[0])-grid[0],_np.arange(grid[0]))/size[0]
@@ -44,6 +48,8 @@ def curl(size,field):
     ----------
     size : numpy.ndarray of shape (3)
         physical size of the periodic field.
+    field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
+        periodic field of which the curl is calculated.
 
     """
     n = _np.prod(field.shape[3:])
@@ -68,6 +74,8 @@ def divergence(size,field):
     ----------
     size : numpy.ndarray of shape (3)
         physical size of the periodic field.
+    field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
+        periodic field of which the divergence is calculated.
 
     """
     n = _np.prod(field.shape[3:])
@@ -82,12 +90,14 @@ def divergence(size,field):
 
 def gradient(size,field):
     """
-    Calculate gradient of a vector or scalar field in Fourier space.
+    Calculate gradient of a scalar or tensor field in Fourier space.
 
     Parameters
     ----------
     size : numpy.ndarray of shape (3)
         physical size of the periodic field.
+    field : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3)
+        periodic field of which the gradient is calculated.
 
     """
     n = _np.prod(field.shape[3:])
@@ -394,9 +404,9 @@ def regrid(size,F,new_grid):
     ----------
     size : numpy.ndarray of shape (3)
         physical size
-    F : numpy.ndarray
+    F : numpy.ndarray of shape (:,:,:,3,3)
         deformation gradient field
-    new_grid : numpy.ndarray
+    new_grid : numpy.ndarray of shape (3)
         new grid for undeformed coordinates
 
     """

From 528378eba22536d7b377db1de9f47d68923f51ce Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 22:51:51 +0200
Subject: [PATCH 15/19] polishing/bugfix

* HDF5 data has x fast, therefore reshape needs to be aware of 'F' order
* Datatype for True/False is 'bool' in python, not 'Boolean'
---
 python/damask/_result.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/python/damask/_result.py b/python/damask/_result.py
index ee7d696a9..65ac7a01b 100644
--- a/python/damask/_result.py
+++ b/python/damask/_result.py
@@ -105,7 +105,7 @@ class Result:
             select from 'set', 'add', and 'del'
         what : str
             attribute to change (must be from self.selection)
-        datasets : list of str or Boolean
+        datasets : list of str or bool
            name of datasets as list, supports ? and * wildcards.
             True is equivalent to [*], False is equivalent to []
 
@@ -197,7 +197,7 @@ class Result:
         ----------
         what : str
             attribute to change (must be from self.selection)
-        datasets : list of str or Boolean
+        datasets : list of str or bool
             name of datasets as list, supports ? and * wildcards.
             True is equivalent to [*], False is equivalent to []
 
@@ -213,7 +213,7 @@ class Result:
         ----------
         what : str
             attribute to change (must be from self.selection)
-        datasets : list of str or Boolean
+        datasets : list of str or bool
             name of datasets as list, supports ? and * wildcards.
             True is equivalent to [*], False is equivalent to []
 
@@ -229,7 +229,7 @@ class Result:
         ----------
         what : str
             attribute to change (must be from self.selection)
-        datasets : list of str or Boolean
+        datasets : list of str or bool
             name of datasets as list, supports ? and * wildcards.
             True is equivalent to [*], False is equivalent to []
 
@@ -256,10 +256,10 @@ class Result:
           datasets : iterable or str
           component : int
               homogenization component to consider for constituent data
-          tagged : Boolean
+          tagged : bool
               tag Table.column name with '#component'
               defaults to False
-          split : Boolean
+          split : bool
               split Table by increment and return dictionary of Tables
               defaults to True
 
@@ -320,7 +320,7 @@ class Result:
 
         Parameters
         ----------
-            datasets : iterable or str or Boolean
+            datasets : iterable or str or bool
 
         Examples
         --------
@@ -454,7 +454,7 @@ class Result:
     def cell_coordinates(self):
         """Return initial coordinates of the cell centers."""
         if self.structured:
-            return grid_filters.cell_coord0(self.grid,self.size,self.origin).reshape(-1,3)
+            return grid_filters.cell_coord0(self.grid,self.size,self.origin).reshape(-1,3,order='F')
         else:
             with h5py.File(self.fname,'r') as f:
                 return f['geometry/x_c'][()]

From 7084fc9616d5833fc554faf471c2d2d3bb3754d4 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 20 Apr 2020 23:19:38 +0200
Subject: [PATCH 16/19] adjusting to consistent data layout

---
 processing/post/blowUp.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/processing/post/blowUp.py b/processing/post/blowUp.py
index 316b74753..23c6b2ef2 100755
--- a/processing/post/blowUp.py
+++ b/processing/post/blowUp.py
@@ -59,13 +59,13 @@ for name in filenames:
   packing = np.array(options.packing,'i')
   outSize = grid*packing
   
-  data = table.data.values.reshape(tuple(grid)+(-1,))
-  blownUp = ndimage.interpolation.zoom(data,tuple(packing)+(1,),order=0,mode='nearest').reshape(outSize.prod(),-1)
+  data = table.data.values.reshape(tuple(grid)+(-1,),order='F')
+  blownUp = ndimage.interpolation.zoom(data,tuple(packing)+(1,),order=0,mode='nearest').reshape(outSize.prod(),-1,order='F')
 
   table = damask.Table(blownUp,table.shapes,table.comments)
 
   coords = damask.grid_filters.cell_coord0(outSize,size,origin)
-  table.set(options.pos,coords.reshape(-1,3))
+  table.set(options.pos,coords.reshape(-1,3,order='F'))
   table.set('elem',np.arange(1,outSize.prod()+1))
   
   outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name))

From 4a05cdbbbb0b6fa977ecf083c6ca8c0bddff1e93 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 21 Apr 2020 03:26:26 +0200
Subject: [PATCH 17/19] copy and paste error

---
 python/damask/grid_filters.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py
index f74bc92ac..eaee42924 100644
--- a/python/damask/grid_filters.py
+++ b/python/damask/grid_filters.py
@@ -90,7 +90,7 @@ def divergence(size,field):
 
 def gradient(size,field):
     """
-    Calculate gradient of a scalar or tensor field in Fourier space.
+    Calculate gradient of a scalar or vector field in Fourier space.
 
     Parameters
     ----------

From b26bc78a44f42e9507e017c9f1d8aa7cb5a61c74 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 21 Apr 2020 18:49:50 +0200
Subject: [PATCH 18/19] adopting to new grid layout

---
 processing/post/DADF5_postResults.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py
index 70c883aeb..2230ba563 100755
--- a/processing/post/DADF5_postResults.py
+++ b/processing/post/DADF5_postResults.py
@@ -33,7 +33,7 @@ for filename in options.filenames:
     results = damask.Result(filename)
 
     if not results.structured: continue
-    coords = damask.grid_filters.cell_coord0(results.grid,results.size,results.origin)
+    coords = damask.grid_filters.cell_coord0(results.grid,results.size,results.origin).reshape(-1,3,order='F')
 
     N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1
     N_digits = 5 # hack to keep test intact

From 936ce6a160020e2b038805e8091757665362ce8e Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 22 Apr 2020 18:26:10 +0200
Subject: [PATCH 19/19] reshape for scalar/vector not needed

---
 processing/post/addGradient.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py
index c6df0eacf..f63f24789 100755
--- a/processing/post/addGradient.py
+++ b/processing/post/addGradient.py
@@ -49,7 +49,7 @@ for name in filenames:
     for label in options.labels:
         field = table.get(label)
         shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,)                          # scalar or vector
-        field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape)
+        field = field.reshape(tuple(grid)+(-1,),order='F')
         grad  = damask.grid_filters.gradient(size,field)
         table.add('gradFFT({})'.format(label),
                   grad.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)*3,order='F'),