Merge branch 'development' into less-shell-scripts

.gitlab-ci.yml

@@ -203,7 +203,6 @@ Post_OrientationConversion:
   stage: postprocessing
   script:
     - OrientationConversion/test.py
-    - OrientationConversion/test2.py
   except:
     - master
     - release

VERSION

@@ -1 +1 @@
-v2.0.3-2303-g2a6132b7
+v2.0.3-2390-g524706d3

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

processing/post/addCompatibilityMismatch.py

@@ -17,7 +17,7 @@ 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.
@@ -67,14 +67,13 @@ def volumeMismatch(size,F,nodes):
   (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,8 +90,7 @@ 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):
@@ -103,35 +101,34 @@ 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(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]]
-  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
-  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
 
 
@@ -178,20 +175,20 @@ 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.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)

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)

processing/post/addDisplacement.py

@@ -5,8 +5,6 @@ import sys
 from io import StringIO
 from optparse import OptionParser
 
-import numpy as np
-
 import damask
 
 
@@ -52,22 +50,22 @@ 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)

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)

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')
+        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)

processing/post/addOrientations.py

@@ -172,7 +172,7 @@ for name in filenames:
       
     elif inputtype == 'matrix':
       d = representations['matrix'][1]
-      o = damask.Rotation.fromMatrix(list(map(float,table.data[column:column+d])))
+      o = damask.Rotation.fromMatrix(np.array(list(map(float,table.data[column:column+d]))).reshape(3,3))
 
     elif inputtype == 'frame':
       M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \

processing/post/addSchmidfactors.py

@@ -214,7 +214,7 @@ for name in filenames:
 
   outputAlive = True
   while outputAlive and table.data_read():                                                          # read next data line of ASCII table
-    o = damask.Rotation(list(map(float,table.data[column:column+4])))
+    o = damask.Rotation(np.array(list(map(float,table.data[column:column+4]))))
 
     table.data_append(  np.abs(  np.sum(slip_direction * (o * force) ,axis=1) \
                                * np.sum(slip_normal    * (o * normal),axis=1)))

processing/post/averageDown.py

@@ -91,7 +91,7 @@ for name in filenames:
   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))

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))

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)
 

processing/pre/seeds_fromGeom.py

@@ -54,7 +54,7 @@ for name in filenames:
                           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:]),

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:

python/.coveragerc

@@ -1,2 +1,5 @@
 [run]
 omit = tests/*
+       damask/_asciitable.py
+       damask/_test.py
+       damask/config/*

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
@@ -352,16 +351,15 @@ 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')
-
+            coords  = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3)
         else:
-            weights_p = weights.flatten()
+            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])
@@ -396,7 +394,7 @@ 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)
 

python/damask/_orientation.py

@@ -38,6 +38,9 @@ class Orientation:
         else:
             self.rotation = Rotation.fromQuaternion(rotation)                                       # assume quaternion
 
+        if self.rotation.quaternion.shape != (4,):
+            raise NotImplementedError('Support for multiple rotations missing')
+
     def disorientation(self,
                        other,
                        SST = True,

python/damask/_result.py

@@ -111,7 +111,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 []
 
@@ -203,7 +203,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 []
 
@@ -219,7 +219,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 []
 
@@ -235,7 +235,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 []
 
@@ -262,10 +262,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
 
@@ -326,7 +326,7 @@ class Result:
 
         Parameters
         ----------
-            datasets : iterable or str or Boolean
+            datasets : iterable or str or bool
 
         Examples
         --------
@@ -460,7 +460,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'][()]

python/damask/_rotation.py

@@ -1,6 +1,7 @@
 import numpy as np
 
 from ._Lambert import ball_to_cube, cube_to_ball
+from . import mechanics
 
 _P = -1
 
@@ -61,6 +62,8 @@ class Rotation:
 
     def __repr__(self):
         """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
+        if self.quaternion.shape != (4,):
+            raise NotImplementedError('Support for multiple rotations missing')
         return '\n'.join([
                'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
                'Matrix:\n{}'.format(self.asMatrix()),
@@ -83,6 +86,8 @@ class Rotation:
         considere rotation of (3,3,3,3)-matrix
 
         """
+        if self.quaternion.shape != (4,):
+            raise NotImplementedError('Support for multiple rotations missing')
         if isinstance(other, Rotation):                                                             # rotate a rotation
             self_q  = self.quaternion[0]
             self_p  = self.quaternion[1:]
@@ -107,7 +112,7 @@ class Rotation:
             elif other.shape == (3,3,):                                                             # rotate a single (3x3)-matrix
                 return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
             elif other.shape == (3,3,3,3,):
-                raise NotImplementedError
+                raise NotImplementedError('Support for rotation of 4th order tensors missing')
             else:
                 return NotImplemented
         else:
@@ -116,7 +121,7 @@ class Rotation:
 
     def inverse(self):
         """In-place inverse rotation/backward rotation."""
-        self.quaternion[1:] *= -1
+        self.quaternion[...,1:] *= -1
         return self
 
     def inversed(self):
@@ -125,12 +130,12 @@ class Rotation:
 
 
     def standardize(self):
-        """In-place quaternion representation with positive q."""
-        if self.quaternion[0] < 0.0: self.quaternion*=-1
+        """In-place quaternion representation with positive real part."""
+        self.quaternion[self.quaternion[...,0] < 0.0] *= -1
         return self
 
     def standardized(self):
-        """Quaternion representation with positive q."""
+        """Quaternion representation with positive real part."""
         return self.copy().standardize()
 
 
@@ -157,15 +162,17 @@ class Rotation:
             Rotation from which the average is rotated.
 
         """
+        if self.quaternion.shape != (4,) or other.quaternion.shape != (4,):
+            raise NotImplementedError('Support for multiple rotations missing')
         return Rotation.fromAverage([self,other])
 
 
     ################################################################################################
     # convert to different orientation representations (numpy arrays)
 
-    def asQuaternion(self):
+    def as_quaternion(self):
         """
-        Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object.
+        Unit quaternion [q, p_1, p_2, p_3].
 
         Parameters
         ----------
@@ -175,8 +182,8 @@ class Rotation:
         """
         return self.quaternion
 
-    def asEulers(self,
-                 degrees = False):
+    def as_Eulers(self,
+                  degrees = False):
         """
         Bunge-Euler angles: (φ_1, ϕ, φ_2).
 
@@ -190,9 +197,9 @@ class Rotation:
         if degrees: eu = np.degrees(eu)
         return eu
 
-    def asAxisAngle(self,
-                    degrees = False,
-                    pair = False):
+    def as_axis_angle(self,
+                      degrees = False,
+                      pair = False):
         """
         Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω).
 
@@ -205,15 +212,15 @@ class Rotation:
 
         """
         ax = Rotation.qu2ax(self.quaternion)
-        if degrees: ax[3] = np.degrees(ax[3])
-        return (ax[:3],ax[3]) if pair else ax
+        if degrees: ax[...,3] = np.degrees(ax[...,3])
+        return (ax[...,:3],ax[...,3]) if pair else ax
 
-    def asMatrix(self):
+    def as_matrix(self):
         """Rotation matrix."""
         return Rotation.qu2om(self.quaternion)
 
-    def asRodrigues(self,
-                    vector = False):
+    def as_Rodrigues(self,
+                     vector = False):
         """
         Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2).
 
@@ -224,9 +231,9 @@ class Rotation:
 
         """
         ro = Rotation.qu2ro(self.quaternion)
-        return ro[:3]*ro[3] if vector else ro
+        return ro[...,:3]*ro[...,3] if vector else ro
 
-    def asHomochoric(self):
+    def as_homochoric(self):
         """Homochoric vector: (h_1, h_2, h_3)."""
         return Rotation.qu2ho(self.quaternion)
 
@@ -234,7 +241,7 @@ class Rotation:
         """Cubochoric vector: (c_1, c_2, c_3)."""
         return Rotation.qu2cu(self.quaternion)
 
-    def asM(self):
+    def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M
         """
         Intermediate representation supporting quaternion averaging.
 
@@ -244,114 +251,133 @@ class Rotation:
         https://doi.org/10.2514/1.28949
 
         """
-        return np.outer(self.quaternion,self.quaternion)
+        return np.einsum('...i,...j',self.quaternion,self.quaternion)
 
+    # for compatibility (old names do not follow convention)
+    asM          = M
+    asQuaternion = as_quaternion
+    asEulers     = as_Eulers
+    asAxisAngle  = as_axis_angle
+    asMatrix     = as_matrix
+    asRodrigues  = as_Rodrigues
+    asHomochoric = as_homochoric
 
     ################################################################################################
-    # static constructors. The input data needs to follow the convention, options allow to
-    # relax these convections
+    # Static constructors. The input data needs to follow the conventions, options allow to
+    # relax the conventions.
     @staticmethod
-    def fromQuaternion(quaternion,
-                       acceptHomomorph = False,
-                       P = -1):
+    def from_quaternion(quaternion,
+                        acceptHomomorph = False,
+                        P = -1):
 
-        qu =   quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \
-                          else np.array(quaternion,dtype=float)
-        if P > 0: qu[1:4] *= -1                                                                     # convert from P=1 to P=-1
-        if qu[0] < 0.0:
-            if acceptHomomorph:
-                qu *= -1.
-            else:
-                raise ValueError('Quaternion has negative first component: {}.'.format(qu[0]))
-        if not np.isclose(np.linalg.norm(qu), 1.0):
-            raise ValueError('Quaternion is not of unit length: {} {} {} {}.'.format(*qu))
+        qu = np.array(quaternion,dtype=float)
+        if qu.shape[:-2:-1] != (4,):
+            raise ValueError('Invalid shape.')
+
+        if P > 0: qu[...,1:4] *= -1                                                                 # convert from P=1 to P=-1
+        if acceptHomomorph:
+            qu[qu[...,0] < 0.0] *= -1
+        else:
+            if np.any(qu[...,0] < 0.0):
+                raise ValueError('Quaternion with negative first (real) component.')
+        if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0)):
+            raise ValueError('Quaternion is not of unit length.')
 
         return Rotation(qu)
 
     @staticmethod
-    def fromEulers(eulers,
-                   degrees = False):
+    def from_Eulers(eulers,
+                    degrees = False):
+
+        eu = np.array(eulers,dtype=float)
+        if eu.shape[:-2:-1] != (3,):
+            raise ValueError('Invalid shape.')
 
-        eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \
-                    else np.array(eulers,dtype=float)
         eu = np.radians(eu) if degrees else eu
-        if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi:
-            raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]: {} {} {}.'.format(*eu))
+        if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or np.any(eu[...,1] > np.pi):                 # ToDo: No separate check for PHI
+            raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].')
 
         return Rotation(Rotation.eu2qu(eu))
 
     @staticmethod
-    def fromAxisAngle(angleAxis,
-                      degrees = False,
-                      normalise = False,
-                      P = -1):
+    def from_axis_angle(axis_angle,
+                        degrees = False,
+                        normalise = False,
+                        P = -1):
 
-        ax = angleAxis if isinstance(angleAxis, np.ndarray) and angleAxis.dtype == np.dtype(float) \
-                       else np.array(angleAxis,dtype=float)
-        if P > 0:     ax[0:3] *= -1                                                                 # convert from P=1 to P=-1
-        if degrees:   ax[  3]  = np.radians(ax[3])
-        if normalise: ax[0:3] /= np.linalg.norm(ax[0:3])
-        if ax[3] < 0.0 or ax[3] > np.pi:
-            raise ValueError('Axis angle rotation angle outside of [0..π]: {}.'.format(ax[3]))
-        if not np.isclose(np.linalg.norm(ax[0:3]), 1.0):
-            raise ValueError('Axis angle rotation axis is not of unit length: {} {} {}.'.format(*ax[0:3]))
+        ax = np.array(axis_angle,dtype=float)
+        if ax.shape[:-2:-1] != (4,):
+            raise ValueError('Invalid shape.')
+
+        if P > 0:     ax[...,0:3] *= -1                                                             # convert from P=1 to P=-1
+        if degrees:   ax[...,  3]  = np.radians(ax[...,3])
+        if normalise: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1)
+        if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi):
+            raise ValueError('Axis angle rotation angle outside of [0..π].')
+        if not np.all(np.isclose(np.linalg.norm(ax[...,0:3],axis=-1), 1.0)):
+            raise ValueError('Axis angle rotation axis is not of unit length.')
 
         return Rotation(Rotation.ax2qu(ax))
 
     @staticmethod
-    def fromBasis(basis,
-                  orthonormal = True,
-                  reciprocal = False,
-                 ):
+    def from_basis(basis,
+                   orthonormal = True,
+                   reciprocal = False):
+
+        om = np.array(basis,dtype=float)
+        if om.shape[:-3:-1] != (3,3):
+            raise ValueError('Invalid shape.')
 
-        om = basis if isinstance(basis, np.ndarray) else np.array(basis).reshape(3,3)
         if reciprocal:
-            om = np.linalg.inv(om.T/np.pi)                                                          # transform reciprocal basis set
+            om = np.linalg.inv(mechanics.transpose(om)/np.pi)                                       # transform reciprocal basis set
             orthonormal = False                                                                     # contains stretch
         if not orthonormal:
             (U,S,Vh) = np.linalg.svd(om)                                                            # singular value decomposition
-            om = np.dot(U,Vh)
-        if not np.isclose(np.linalg.det(om),1.0):
-            raise ValueError('matrix is not a proper rotation: {}.'.format(om))
-        if    not np.isclose(np.dot(om[0],om[1]), 0.0) \
-           or not np.isclose(np.dot(om[1],om[2]), 0.0) \
-           or not np.isclose(np.dot(om[2],om[0]), 0.0):
-            raise ValueError('matrix is not orthogonal: {}.'.format(om))
+            om = np.einsum('...ij,...jl->...il',U,Vh)
+        if not np.all(np.isclose(np.linalg.det(om),1.0)):
+            raise ValueError('Orientation matrix has determinant ≠ 1.')
+        if    not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \
+           or not np.all(np.isclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0)) \
+           or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)):
+            raise ValueError('Orientation matrix is not orthogonal.')
 
         return Rotation(Rotation.om2qu(om))
 
     @staticmethod
-    def fromMatrix(om,
-                  ):
+    def from_matrix(om):
 
-        return Rotation.fromBasis(om)
+        return Rotation.from_basis(om)
 
     @staticmethod
-    def fromRodrigues(rodrigues,
-                      normalise = False,
-                      P = -1):
+    def from_Rodrigues(rodrigues,
+                       normalise = False,
+                       P = -1):
 
-        ro = rodrigues if isinstance(rodrigues, np.ndarray) and rodrigues.dtype == np.dtype(float) \
-                       else np.array(rodrigues,dtype=float)
-        if P > 0:     ro[0:3] *= -1                                                                 # convert from P=1 to P=-1
-        if normalise: ro[0:3] /= np.linalg.norm(ro[0:3])
-        if not np.isclose(np.linalg.norm(ro[0:3]), 1.0):
-            raise ValueError('Rodrigues rotation axis is not of unit length: {} {} {}.'.format(*ro[0:3]))
-        if ro[3] < 0.0:
-            raise ValueError('Rodrigues rotation angle not positive: {}.'.format(ro[3]))
+        ro = np.array(rodrigues,dtype=float)
+        if ro.shape[:-2:-1] != (4,):
+            raise ValueError('Invalid shape.')
+
+        if P > 0:     ro[...,0:3] *= -1                                                             # convert from P=1 to P=-1
+        if normalise: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1)
+        if np.any(ro[...,3] < 0.0):
+            raise ValueError('Rodrigues vector rotation angle not positive.')
+        if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)):
+            raise ValueError('Rodrigues vector rotation axis is not of unit length.')
 
         return Rotation(Rotation.ro2qu(ro))
 
     @staticmethod
-    def fromHomochoric(homochoric,
-                       P = -1):
+    def from_homochoric(homochoric,
+                        P = -1):
+
+        ho = np.array(homochoric,dtype=float)
+        if ho.shape[:-2:-1] != (3,):
+            raise ValueError('Invalid shape.')
 
-        ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
-                        else np.array(homochoric,dtype=float)
         if P > 0: ho *= -1                                                                          # convert from P=1 to P=-1
 
-        if np.linalg.norm(ho) > (3.*np.pi/4.)**(1./3.)+1e-9:
-            raise ValueError('Coordinate outside of the sphere: {} {} {}.'.format(ho))
+        if np.any(np.linalg.norm(ho,axis=-1) > (3.*np.pi/4.)**(1./3.)+1e-9):
+            raise ValueError('Homochoric coordinate outside of the sphere.')
 
         return Rotation(Rotation.ho2qu(ho))
 
@@ -359,11 +385,12 @@ class Rotation:
     def fromCubochoric(cubochoric,
                        P = -1):
 
-        cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \
-                        else np.array(cubochoric,dtype=float)
+        cu = np.array(cubochoric,dtype=float)
+        if cu.shape[:-2:-1] != (3,):
+            raise ValueError('Invalid shape.')
 
         if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9:
-            raise ValueError('Coordinate outside of the cube: {} {} {}.'.format(*cu))
+            raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu))
 
         ho = Rotation.cu2ho(cu)
         if P > 0: ho *= -1                                                                          # convert from P=1 to P=-1
@@ -403,17 +430,34 @@ class Rotation:
 
         return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
 
-
     @staticmethod
-    def fromRandom():
-        r = np.random.random(3)
-        A = np.sqrt(r[2])
-        B = np.sqrt(1.0-r[2])
-        return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A,
-                                  np.sin(2.0*np.pi*r[1])*B,
-                                  np.cos(2.0*np.pi*r[1])*B,
-                                  np.sin(2.0*np.pi*r[0])*A])).standardize()
+    def from_random(shape=None):
+        if shape is None:
+            r = np.random.random(3)
+        elif hasattr(shape, '__iter__'):
+            r = np.random.random(tuple(shape)+(3,))
+        else:
+            r = np.random.random((shape,3))
 
+        A = np.sqrt(r[...,2])
+        B = np.sqrt(1.0-r[...,2])
+        q = np.stack([np.cos(2.0*np.pi*r[...,0])*A,
+                      np.sin(2.0*np.pi*r[...,1])*B,
+                      np.cos(2.0*np.pi*r[...,1])*B,
+                      np.sin(2.0*np.pi*r[...,0])*A],axis=-1)
+
+        return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q).standardize()
+
+
+    # for compatibility (old names do not follow convention)
+    fromQuaternion = from_quaternion
+    fromEulers     = from_Eulers
+    fromAxisAngle  = from_axis_angle
+    fromBasis      = from_basis
+    fromMatrix     = from_matrix
+    fromRodrigues  = from_Rodrigues
+    fromHomochoric = from_homochoric
+    fromRandom     = from_random
 
 ####################################################################################################
 # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations
@@ -808,12 +852,11 @@ class Rotation:
                 c = np.cos(ax[3]*0.5)
                 s = np.sin(ax[3]*0.5)
                 qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ])
-            return qu
         else:
             c = np.cos(ax[...,3:4]*.5)
             s = np.sin(ax[...,3:4]*.5)
             qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s]))
-            return qu
+        return qu
 
     @staticmethod
     def ax2om(ax):
@@ -859,7 +902,7 @@ class Rotation:
                 # 180 degree case
                 ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \
                       [np.tan(ax[3]*0.5)]
-            return np.array(ro)
+            ro = np.array(ro)
         else:
             ro = np.block([ax[...,:3],
                            np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0),
@@ -867,7 +910,7 @@ class Rotation:
                                     np.tan(ax[...,3:4]*0.5))
                           ])
             ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0]
-            return ro
+        return ro
 
     @staticmethod
     def ax2ho(ax):
@@ -875,11 +918,10 @@ class Rotation:
         if len(ax.shape) == 1:
             f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0)
             ho = ax[0:3] * f
-            return ho
         else:
             f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0)
             ho = ax[...,:3] * f
-            return ho
+        return ho
 
     @staticmethod
     def ax2cu(ax):
@@ -936,7 +978,6 @@ class Rotation:
             f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi)
             ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape),
                           np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0))
-
         return ho
 
     @staticmethod
@@ -1010,7 +1051,7 @@ class Rotation:
         if len(ho.shape) == 1:
             return ball_to_cube(ho)
         else:
-            raise NotImplementedError
+            raise NotImplementedError('Support for multiple rotations missing')
 
 
     #---------- Cubochoric ----------
@@ -1045,4 +1086,4 @@ class Rotation:
         if len(cu.shape) == 1:
             return cube_to_ball(cu)
         else:
-            raise NotImplementedError
+            raise NotImplementedError('Support for multiple rotations missing')

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 data:
+
+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
 
@@ -7,8 +21,12 @@ def _ks(size,grid,first_order=False):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    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]
@@ -19,8 +37,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,8 +46,10 @@ def curl(size,field):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    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:])
@@ -53,8 +72,10 @@ def divergence(size,field):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    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:])
@@ -69,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 vector field in Fourier space.
 
     Parameters
     ----------
-    size : numpy.ndarray
+    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:])
@@ -93,9 +116,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 +126,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 +139,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 +166,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 +182,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,25 +197,25 @@ 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).
 
@@ -211,13 +238,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[::-1])+(3,)),cell_coord0(grid,size,origin)):
-        raise ValueError('Input data is not a regular grid.')
+    if ordered and not _np.allclose(coord0.reshape(tuple(grid)+(3,),order='F'),cell_coord0(grid,size,origin)):
+        raise ValueError('Input data is not ordered (x fast, z slow).')
 
     return (grid,size,origin)
 
@@ -241,17 +268,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 +288,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 +303,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 +319,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 +334,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,14 +363,14 @@ 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).
 
@@ -357,13 +385,13 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False):
     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)[::-1])+(3,)),node_coord0(grid,size,origin)):
-        raise ValueError('Input data is not a regular grid.')
+    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 ordered (x fast, z slow).')
 
     return (grid,size,origin)
 
@@ -374,15 +402,15 @@ def regrid(size,F,new_grid):
 
     Parameters
     ----------
-    size : numpy.ndarray
+    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
 
     """
-    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)
 

python/damask/mechanics.py

@@ -135,16 +135,16 @@ def PK2(P,F):
 
     Parameters
     ----------
-    P : numpy.ndarray of shape (:,3,3) or (3,3)
+    P : numpy.ndarray of shape (...,3,3) or (3,3)
         First Piola-Kirchhoff stress.
-    F : numpy.ndarray of shape (:,3,3) or (3,3)
+    F : numpy.ndarray of shape (...,3,3) or (3,3)
         Deformation gradient.
 
     """
     if _np.shape(F) == _np.shape(P) == (3,3):
         S = _np.dot(_np.linalg.inv(F),P)
     else:
-        S = _np.einsum('ijk,ikl->ijl',_np.linalg.inv(F),P)
+        S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P)
     return symmetric(S)
 
 
@@ -241,7 +241,7 @@ def symmetric(T):
 
     Parameters
     ----------
-    T : numpy.ndarray of shape (:,3,3) or (3,3)
+    T : numpy.ndarray of shape (...,3,3) or (3,3)
         Tensor of which the symmetrized values are computed.
 
     """
@@ -254,12 +254,12 @@ def transpose(T):
 
     Parameters
     ----------
-    T : numpy.ndarray of shape (:,3,3) or (3,3)
+    T : numpy.ndarray of shape (...,3,3) or (3,3)
         Tensor of which the transpose is computed.
 
     """
     return T.T if _np.shape(T) == (3,3) else \
-           _np.transpose(T,(0,2,1))
+           _np.swapaxes(T,axis2=-2,axis1=-1)
 
 
 def _polar_decomposition(T,requested):

python/tests/test_Rotation.py

@@ -157,6 +157,30 @@ class TestRotation:
             print(m,o,rot.asQuaternion())
             assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9
 
+    @pytest.mark.parametrize('function',[Rotation.from_quaternion,
+                                         Rotation.from_Eulers,
+                                         Rotation.from_axis_angle,
+                                         Rotation.from_matrix,
+                                         Rotation.from_Rodrigues,
+                                         Rotation.from_homochoric])
+    def test_invalid_shape(self,function):
+        invalid_shape = np.random.random(np.random.randint(8,32,(3)))
+        with pytest.raises(ValueError):
+            function(invalid_shape)
+
+    @pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])),
+                                                 (Rotation.from_quaternion, np.array([1,1,1,0])),
+                                                 (Rotation.from_Eulers,     np.array([1,4,0])),
+                                                 (Rotation.from_axis_angle, np.array([1,0,0,4])),
+                                                 (Rotation.from_axis_angle, np.array([1,1,0,1])),
+                                                 (Rotation.from_matrix,     np.random.rand(3,3)),
+                                                 (Rotation.from_Rodrigues,  np.array([1,0,0,-1])),
+                                                 (Rotation.from_Rodrigues,  np.array([1,1,0,1])),
+                                                 (Rotation.from_homochoric, np.array([2,2,2]))  ])
+    def test_invalid(self,function,invalid):
+        with pytest.raises(ValueError):
+            function(invalid)
+
     @pytest.mark.parametrize('conversion',[Rotation.qu2om,
                                            Rotation.qu2eu,
                                            Rotation.qu2ax,

python/tests/test_grid_filters.py

@@ -9,13 +9,13 @@ class TestGridFilters:
          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])
@@ -80,8 +80,43 @@ class TestGridFilters:
          F    = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3))
          assert np.allclose(function(size,F),0.0)
 
+    @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):
+            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)
          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()))

src/CPFEM.f90

@@ -73,28 +73,24 @@ subroutine CPFEM_initAll(el,ip)
   integer(pInt), intent(in) :: el, &                                                                !< FE el number
                                ip                                                                   !< FE integration point number
 
-  !$OMP CRITICAL(init)
-  if (.not. CPFEM_init_done) then
-    call DAMASK_interface_init
-    call prec_init
-    call IO_init
-    call numerics_init
-    call debug_init
-    call config_init
-    call math_init
-    call rotations_init
-    call HDF5_utilities_init
-    call results_init
-    call discretization_marc_init(ip, el)
-    call lattice_init
-    call material_init
-    call constitutive_init
-    call crystallite_init
-    call homogenization_init
-    call CPFEM_init
-    CPFEM_init_done = .true.
-  endif
-  !$OMP END CRITICAL(init)
+  CPFEM_init_done = .true.
+  call DAMASK_interface_init
+  call prec_init
+  call IO_init
+  call numerics_init
+  call debug_init
+  call config_init
+  call math_init
+  call rotations_init
+  call HDF5_utilities_init
+  call results_init
+  call discretization_marc_init(ip, el)
+  call lattice_init
+  call material_init
+  call constitutive_init
+  call crystallite_init
+  call homogenization_init
+  call CPFEM_init
 
 end subroutine CPFEM_initAll
 

src/DAMASK_marc.f90

@@ -261,11 +261,10 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
  endif
 
  !$ defaultNumThreadsInt = omp_get_num_threads()                                                    ! remember number of threads set by Marc
+ !$ call omp_set_num_threads(DAMASK_NumThreadsInt)                                                  ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
 
  if (.not. CPFEM_init_done) call CPFEM_initAll(m(1),nn)
 
- !$ call omp_set_num_threads(DAMASK_NumThreadsInt)                                                  ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
-
  computationMode = 0                                                                                ! save initialization value, since it does not result in any calculation
  if (lovl == 4 ) then                                                                               ! jacobian requested by marc
    if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) &                              ! first after cutback

src/constitutive.f90

@@ -327,7 +327,7 @@ module constitutive
     constitutive_initialFi, &
     constitutive_SandItsTangents, &
     constitutive_collectDotState, &
-    constitutive_collectDeltaState, &
+    constitutive_deltaState, &
     constitutive_results
 
 contains
@@ -709,12 +709,14 @@ end subroutine constitutive_hooke_SandItsTangents
 !--------------------------------------------------------------------------------------------------
 !> @brief contains the constitutive equation for calculating the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
-subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el)
+function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken)
 
   integer, intent(in) :: &
     ipc, &                                                                                          !< component-ID of integration point
     ip, &                                                                                           !< integration point
-    el                                                                                              !< element
+    el, &                                                                                              !< element
+    phase, &
+    of
   real(pReal),  intent(in) :: &
     subdt                                                                                           !< timestep
   real(pReal),  intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
@@ -730,16 +732,16 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip,
     ho, &                                                                                           !< homogenization
     tme, &                                                                                          !< thermal member position
     i, &                                                                                            !< counter in source loop
-    instance, of
+    instance
+  logical :: broken
 
   ho = material_homogenizationAt(el)
   tme = thermalMapping(ho)%p(ip,el)
-  of = material_phasememberAt(ipc,ip,el)
-  instance = phase_plasticityInstance(material_phaseAt(ipc,el))
+  instance = phase_plasticityInstance(phase)
 
   Mp = matmul(matmul(transpose(Fi),Fi),S)
 
-  plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
+  plasticityType: select case (phase_plasticity(phase))
 
     case (PLASTICITY_ISOTROPIC_ID) plasticityType
       call plastic_isotropic_dotState    (Mp,instance,of)
@@ -760,10 +762,11 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip,
       call plastic_nonlocal_dotState     (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, &
                                           instance,of,ip,el)
   end select plasticityType
+  broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of)))
 
-  SourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el))
+  SourceLoop: do i = 1, phase_Nsources(phase)
 
-    sourceType: select case (phase_source(i,material_phaseAt(ipc,el)))
+    sourceType: select case (phase_source(i,phase))
 
       case (SOURCE_damage_anisoBrittle_ID) sourceType
         call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress?
@@ -775,25 +778,29 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip,
         call source_damage_anisoDuctile_dotState (   ipc, ip, el)
 
       case (SOURCE_thermal_externalheat_ID) sourceType
-        call source_thermal_externalheat_dotState(material_phaseAt(ipc,el),of)
+        call source_thermal_externalheat_dotState(phase,of)
 
     end select sourceType
 
+    broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of)))
+
   enddo SourceLoop
 
-end subroutine constitutive_collectDotState
+end function constitutive_collectDotState
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief for constitutive models having an instantaneous change of state
 !> will return false if delta state is not needed/supported by the constitutive model
 !--------------------------------------------------------------------------------------------------
-subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
+function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broken)
 
   integer, intent(in) :: &
     ipc, &                                                                                          !< component-ID of integration point
     ip, &                                                                                           !< integration point
-    el                                                                                              !< element
+    el, &                                                                                           !< element
+    phase, &
+    of
   real(pReal),   intent(in), dimension(3,3) :: &
     S, &                                                                                            !< 2nd Piola Kirchhoff stress
     Fe, &                                                                                           !< elastic deformation gradient
@@ -802,35 +809,62 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
     Mp
   integer :: &
     i, &
-    instance, of
+    instance, &
+    myOffset, &
+    mySize
+  logical :: &
+    broken
 
   Mp  = matmul(matmul(transpose(Fi),Fi),S)
-  of = material_phasememberAt(ipc,ip,el)
-  instance = phase_plasticityInstance(material_phaseAt(ipc,el))
+  instance = phase_plasticityInstance(phase)
 
-  plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
+  plasticityType: select case (phase_plasticity(phase))
 
     case (PLASTICITY_KINEHARDENING_ID) plasticityType
       call plastic_kinehardening_deltaState(Mp,instance,of)
+      broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
 
     case (PLASTICITY_NONLOCAL_ID) plasticityType
       call plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
+      broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of)))
+
+    case default
+      broken = .false.
 
   end select plasticityType
 
-  sourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el))
+  if(.not. broken) then
+    select case(phase_plasticity(phase))
+      case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
 
-     sourceType: select case (phase_source(i,material_phaseAt(ipc,el)))
+        myOffset = plasticState(phase)%offsetDeltaState
+        mySize   = plasticState(phase)%sizeDeltaState
+        plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = &
+        plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of)
+    end select
+  endif
+
+
+  sourceLoop: do i = 1, phase_Nsources(phase)
+
+     sourceType: select case (phase_source(i,phase))
 
       case (SOURCE_damage_isoBrittle_ID) sourceType
         call source_damage_isoBrittle_deltaState  (constitutive_homogenizedC(ipc,ip,el), Fe, &
                                                    ipc, ip, el)
+        broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of)))
+        if(.not. broken) then
+          myOffset = sourceState(phase)%p(i)%offsetDeltaState
+          mySize   = sourceState(phase)%p(i)%sizeDeltaState
+          sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) = &
+          sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) + sourceState(phase)%p(i)%deltaState(1:mySize,of)
+        endif
 
     end select sourceType
 
   enddo SourceLoop
 
-end subroutine constitutive_collectDeltaState
+end function constitutive_deltaState
 
 
 !--------------------------------------------------------------------------------------------------

src/constitutive_plastic_disloUCLA.f90

@@ -209,7 +209,7 @@ module subroutine plastic_disloUCLA_init
     sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
     sizeState = sizeDotState
 
-    call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
+    call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization

src/constitutive_plastic_dislotwin.f90

@@ -399,7 +399,7 @@ module subroutine plastic_dislotwin_init
                  + size(['f_tr'])                           * prm%sum_N_tr
     sizeState = sizeDotState
 
-    call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
+    call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! locally defined state aliases and initialization of state0 and atol

src/constitutive_plastic_isotropic.f90

@@ -117,7 +117,7 @@ module subroutine plastic_isotropic_init
     sizeDotState = size(['xi               ','accumulated_shear'])
     sizeState = sizeDotState
 
-    call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
+    call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization

src/constitutive_plastic_kinehardening.f90

@@ -164,7 +164,7 @@ module subroutine plastic_kinehardening_init
     sizeDeltaState = size(['sense ',   'chi0  ',    'gamma0'   ]) * prm%sum_N_sl
     sizeState = sizeDotState + sizeDeltaState
 
-    call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
+    call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization

src/constitutive_plastic_none.f90

@@ -29,7 +29,7 @@ module subroutine plastic_none_init
     if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
 
     NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
-    call material_allocatePlasticState(p,NipcMyPhase,0,0,0)
+    call material_allocateState(plasticState(p),NipcMyPhase,0,0,0)
 
   enddo
 

src/constitutive_plastic_nonlocal.f90

@@ -320,6 +320,7 @@ module subroutine plastic_nonlocal_init
       prm%fEdgeMultiplication         = config%getFloat('edgemultiplication')
       prm%shortRangeStressCorrection  = config%keyExists('/shortrangestresscorrection/')
 
+
 !--------------------------------------------------------------------------------------------------
 !  sanity checks
       if (any(prm%burgers          <  0.0_pReal)) extmsg = trim(extmsg)//' burgers'
@@ -384,9 +385,9 @@ module subroutine plastic_nonlocal_init
                                 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl     !< other dependent state variables that are not updated by microstructure
     sizeDeltaState            = sizeDotState
 
-    call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
+    call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
 
-    plasticState(p)%nonlocal = .true.
+    plasticState(p)%nonlocal         = config%KeyExists('/nonlocal/')
     plasticState(p)%offsetDeltaState = 0                                                            ! ToDo: state structure does not follow convention
 
     st0%rho => plasticState(p)%state0                             (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)

src/constitutive_plastic_phenopowerlaw.f90

@@ -213,7 +213,7 @@ module subroutine plastic_phenopowerlaw_init
                  + size(['xi_tw   ','gamma_tw']) * prm%sum_N_tw
     sizeState = sizeDotState
 
-    call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
+    call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization

src/material.f90

@@ -11,7 +11,6 @@ module material
   use results
   use IO
   use debug
-  use numerics
   use rotations
   use discretization
 
@@ -174,8 +173,7 @@ module material
 
   public :: &
     material_init, &
-    material_allocatePlasticState, &
-    material_allocateSourceState, &
+    material_allocateState, &
     ELASTICITY_HOOKE_ID ,&
     PLASTICITY_NONE_ID, &
     PLASTICITY_ISOTROPIC_ID, &
@@ -700,63 +698,35 @@ end subroutine material_parseTexture
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief allocates the plastic state of a phase
+!> @brief Allocate the components of the state structure for a given phase
 !--------------------------------------------------------------------------------------------------
-subroutine material_allocatePlasticState(phase,NipcMyPhase,&
-                                         sizeState,sizeDotState,sizeDeltaState)
+subroutine material_allocateState(state, &
+                                  NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
 
+  class(tState), intent(out) :: &
+    state
   integer, intent(in) :: &
-    phase, &
     NipcMyPhase, &
     sizeState, &
     sizeDotState, &
     sizeDeltaState
 
-  plasticState(phase)%sizeState        = sizeState
-  plasticState(phase)%sizeDotState     = sizeDotState
-  plasticState(phase)%sizeDeltaState   = sizeDeltaState
-  plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState                                   ! deltaState occupies latter part of state by definition
+  state%sizeState        = sizeState
+  state%sizeDotState     = sizeDotState
+  state%sizeDeltaState   = sizeDeltaState
+  state%offsetDeltaState = sizeState-sizeDeltaState                                                 ! deltaState occupies latter part of state by definition
 
-  allocate(plasticState(phase)%atol                (sizeState),               source=0.0_pReal)
-  allocate(plasticState(phase)%state0              (sizeState,NipcMyPhase),   source=0.0_pReal)
-  allocate(plasticState(phase)%partionedState0     (sizeState,NipcMyPhase),   source=0.0_pReal)
-  allocate(plasticState(phase)%subState0           (sizeState,NipcMyPhase),   source=0.0_pReal)
-  allocate(plasticState(phase)%state               (sizeState,NipcMyPhase),   source=0.0_pReal)
+  allocate(state%atol           (sizeState),             source=0.0_pReal)
+  allocate(state%state0         (sizeState,NipcMyPhase), source=0.0_pReal)
+  allocate(state%partionedState0(sizeState,NipcMyPhase), source=0.0_pReal)
+  allocate(state%subState0      (sizeState,NipcMyPhase), source=0.0_pReal)
+  allocate(state%state          (sizeState,NipcMyPhase), source=0.0_pReal)
 
-  allocate(plasticState(phase)%dotState            (sizeDotState,NipcMyPhase),source=0.0_pReal)
+  allocate(state%dotState    (sizeDotState,NipcMyPhase), source=0.0_pReal)
 
-  allocate(plasticState(phase)%deltaState        (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
+  allocate(state%deltaState(sizeDeltaState,NipcMyPhase), source=0.0_pReal)
 
-end subroutine material_allocatePlasticState
+end subroutine material_allocateState
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief allocates the source state of a phase
-!--------------------------------------------------------------------------------------------------
-subroutine material_allocateSourceState(phase,of,NipcMyPhase,&
-                                        sizeState,sizeDotState,sizeDeltaState)
-
-  integer, intent(in) :: &
-    phase, &
-    of, &
-    NipcMyPhase, &
-    sizeState, sizeDotState,sizeDeltaState
-
-  sourceState(phase)%p(of)%sizeState        = sizeState
-  sourceState(phase)%p(of)%sizeDotState     = sizeDotState
-  sourceState(phase)%p(of)%sizeDeltaState   = sizeDeltaState
-  sourceState(phase)%p(of)%offsetDeltaState = sizeState-sizeDeltaState                              ! deltaState occupies latter part of state by definition
-
-  allocate(sourceState(phase)%p(of)%atol                (sizeState),               source=0.0_pReal)
-  allocate(sourceState(phase)%p(of)%state0              (sizeState,NipcMyPhase),   source=0.0_pReal)
-  allocate(sourceState(phase)%p(of)%partionedState0     (sizeState,NipcMyPhase),   source=0.0_pReal)
-  allocate(sourceState(phase)%p(of)%subState0           (sizeState,NipcMyPhase),   source=0.0_pReal)
-  allocate(sourceState(phase)%p(of)%state               (sizeState,NipcMyPhase),   source=0.0_pReal)
-
-  allocate(sourceState(phase)%p(of)%dotState            (sizeDotState,NipcMyPhase),source=0.0_pReal)
-
-  allocate(sourceState(phase)%p(of)%deltaState        (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
-
-end subroutine material_allocateSourceState
-
 end module material

src/numerics.f90

@@ -20,8 +20,7 @@ module numerics
    iJacoStiffness             =  1, &                                                               !< frequency of stiffness update
    randomSeed                 =  0, &                                                               !< fixed seeding for pseudo-random number generator, Default 0: use random seed
    worldrank                  =  0, &                                                               !< MPI worldrank (/=0 for MPI simulations only)
-   worldsize                  =  1, &                                                               !< MPI worldsize (/=1 for MPI simulations only)
-   numerics_integrator        =  1                                                                  !< method used for state integration Default 1: fix-point iteration
+   worldsize                  =  1                                                                  !< MPI worldsize (/=1 for MPI simulations only)
  integer(4), protected, public :: &
    DAMASK_NumThreadsInt       =  0                                                                  !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
  real(pReal), protected, public :: &
@@ -134,8 +133,6 @@ subroutine numerics_init
          defgradTolerance = IO_floatValue(line,chunkPos,2)
        case ('ijacostiffness')
          iJacoStiffness = IO_intValue(line,chunkPos,2)
-       case ('integrator')
-         numerics_integrator = IO_intValue(line,chunkPos,2)
        case ('usepingpong')
          usepingpong = IO_intValue(line,chunkPos,2) > 0
        case ('unitlength')
@@ -176,6 +173,11 @@ subroutine numerics_init
        case ('maxstaggerediter')
          stagItMax = IO_intValue(line,chunkPos,2)
 
+#ifdef PETSC
+       case ('petsc_options')
+         petsc_options = trim(line(chunkPos(4):))
+#endif
+
 !--------------------------------------------------------------------------------------------------
 ! spectral parameters
 #ifdef Grid
@@ -187,8 +189,6 @@ subroutine numerics_init
          err_stress_tolrel = IO_floatValue(line,chunkPos,2)
        case ('err_stress_tolabs')
          err_stress_tolabs = IO_floatValue(line,chunkPos,2)
-       case ('petsc_options')
-         petsc_options = trim(line(chunkPos(4):))
        case ('err_curl_tolabs')
          err_curl_tolAbs = IO_floatValue(line,chunkPos,2)
        case ('err_curl_tolrel')
@@ -206,8 +206,6 @@ subroutine numerics_init
          integrationorder = IO_intValue(line,chunkPos,2)
        case ('structorder')
          structorder = IO_intValue(line,chunkPos,2)
-       case ('petsc_options')
-         petsc_options = trim(line(chunkPos(4):))
        case ('bbarstabilisation')
          BBarStabilisation = IO_intValue(line,chunkPos,2) > 0
 #endif
@@ -223,7 +221,6 @@ subroutine numerics_init
 ! writing parameters to output
  write(6,'(a24,1x,es8.1)')  ' defgradTolerance:       ',defgradTolerance
  write(6,'(a24,1x,i8)')     ' iJacoStiffness:         ',iJacoStiffness
- write(6,'(a24,1x,i8)')     ' integrator:             ',numerics_integrator
  write(6,'(a24,1x,L8)')     ' use ping pong scheme:   ',usepingpong
  write(6,'(a24,1x,es8.1,/)')' unitlength:             ',numerics_unitlength
 
@@ -266,7 +263,6 @@ subroutine numerics_init
  write(6,'(a24,1x,es8.1)')   ' err_curl_tolRel:        ',err_curl_tolRel
  write(6,'(a24,1x,es8.1)')   ' polarAlpha:             ',polarAlpha
  write(6,'(a24,1x,es8.1)')   ' polarBeta:              ',polarBeta
- write(6,'(a24,1x,a)')       ' PETSc_options:          ',trim(petsc_options)
 #endif
 
 !--------------------------------------------------------------------------------------------------
@@ -274,16 +270,17 @@ subroutine numerics_init
 #ifdef FEM
  write(6,'(a24,1x,i8)')      ' integrationOrder:       ',integrationOrder
  write(6,'(a24,1x,i8)')      ' structOrder:            ',structOrder
- write(6,'(a24,1x,a)')       ' PETSc_options:          ',trim(petsc_options)
  write(6,'(a24,1x,L8)')      ' B-Bar stabilisation:    ',BBarStabilisation
 #endif
 
+#ifdef PETSC
+ write(6,'(a24,1x,a)')       ' PETSc_options:          ',trim(petsc_options)
+#endif
+
 !--------------------------------------------------------------------------------------------------
 ! sanity checks
  if (defgradTolerance <= 0.0_pReal)        call IO_error(301,ext_msg='defgradTolerance')
  if (iJacoStiffness < 1)                   call IO_error(301,ext_msg='iJacoStiffness')
- if (numerics_integrator <= 0 .or. numerics_integrator >= 6) &
-                                           call IO_error(301,ext_msg='integrator')
  if (numerics_unitlength <= 0.0_pReal)     call IO_error(301,ext_msg='unitlength')
  if (residualStiffness < 0.0_pReal)        call IO_error(301,ext_msg='residualStiffness')
  if (itmax <= 1)                           call IO_error(301,ext_msg='itmax')

src/prec.f90

@@ -53,8 +53,7 @@ module prec
     logical :: &
       nonlocal = .false.
     real(pReal), pointer,     dimension(:,:) :: &
-      slipRate, &                                                                                   !< slip rate
-      accumulatedSlip                                                                               !< accumulated plastic slip
+      slipRate                                                                                      !< slip rate
   end type
 
   type :: tSourceState

src/source_damage_anisoBrittle.f90

@@ -107,7 +107,7 @@ subroutine source_damage_anisoBrittle_init
     if (any(prm%critDisp <  0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critDisp'
 
     NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
-    call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0)
+    call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
     sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
     if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
 

src/source_damage_anisoDuctile.f90

@@ -89,7 +89,7 @@ subroutine source_damage_anisoDuctile_init
     if (any(prm%critPlasticStrain <  0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain'
 
     NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
-    call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0)
+    call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
     sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisoductile_atol',defaultVal=1.0e-3_pReal)
     if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
 

src/source_damage_isoBrittle.f90

@@ -83,7 +83,7 @@ subroutine source_damage_isoBrittle_init
     if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy'
 
     NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
-    call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,1)
+    call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,1)
     sourceState(p)%p(sourceOffset)%atol = config%getFloat('isobrittle_atol',defaultVal=1.0e-3_pReal)
     if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
 

src/source_damage_isoDuctile.f90

@@ -82,7 +82,7 @@ subroutine source_damage_isoDuctile_init
     if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain'
 
     NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
-    call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0)
+    call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
     sourceState(p)%p(sourceOffset)%atol = config%getFloat('isoductile_atol',defaultVal=1.0e-3_pReal)
     if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
 

src/source_thermal_dissipation.f90

@@ -67,7 +67,7 @@ subroutine source_thermal_dissipation_init
     prm%kappa = config%getFloat('dissipation_coldworkcoeff')
 
     NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
-    call material_allocateSourceState(p,sourceOffset,NipcMyPhase,0,0,0)
+    call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,0,0,0)
 
     end associate
   enddo

src/source_thermal_externalheat.f90

@@ -74,7 +74,7 @@ subroutine source_thermal_externalheat_init
     prm%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(prm%time))
 
     NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
-    call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0)
+    call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
 
     end associate
   enddo