Merge branch 'development' into unified-RK4-RKCK45

.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-2364-g62f7363a

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

@@ -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])
+  vMismatch = np.empty(F.shape[:3])
-  volInitial = size.prod()/grid.prod()
 
 #--------------------------------------------------------------------------------------------------
 # 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(F.shape[:3])
-  sMismatch    = np.empty(grid[::-1])
 
 #--------------------------------------------------------------------------------------------------
 # initial positions
-  coordsInitial[0,0:3] = [-size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
+  delta = size/grid*.5
-  coordsInitial[1,0:3] = [+size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
+  coordsInitial = np.vstack((delta * np.array((-1,-1,-1)),
-  coordsInitial[2,0:3] = [+size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
+                             delta * np.array((+1,-1,-1)),
-  coordsInitial[3,0:3] = [-size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
+                             delta * np.array((+1,+1,-1)),
-  coordsInitial[4,0:3] = [-size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
+                             delta * np.array((-1,+1,-1)),
-  coordsInitial[5,0:3] = [+size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
+                             delta * np.array((-1,-1,+1)),
-  coordsInitial[6,0:3] = [+size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
+                             delta * np.array((+1,-1,+1)),
-  coordsInitial[7,0:3] = [-size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
+                             delta * np.array((+1,+1,+1)),
-  coordsInitial = coordsInitial/2.0
+                             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+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,  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+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/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,))
+  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)
+  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/post/permuteData.py

@@ -97,10 +97,10 @@ for name in filenames:
   for col,dim in zip(columns,dims):
     if options.unique:
       s = set(map(tuple,table.data[:,col:col+dim]))                                                 # generate set of (unique) values
-      uniques = np.array(map(np.array,s))                                                           # translate set to np.array
+      uniques = np.array(list(map(np.array,s)))                                                     # translate set to np.array
       shuffler = dict(zip(s,np.random.permutation(len(s))))                                         # random permutation
-      table.data[:,col:col+dim] = uniques[np.array(map(lambda x: shuffler[tuple(x)],
+      table.data[:,col:col+dim] = uniques[np.array(list(map(lambda x: shuffler[tuple(x)],
-                                                       table.data[:,col:col+dim]))]                 # fill table with mapped uniques
+                                                            table.data[:,col:col+dim])))]           # fill table with mapped uniques
     else:
       np.random.shuffle(table.data[:,col:col+dim])                                                  # independently shuffle every row
 

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/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,1) != 0.0):                                             # no float present
-        if not np.any(np.mod(microstructure.flatten(),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/_result.py

@@ -86,12 +86,18 @@ class Result:
     def __repr__(self):
         """Show selected data."""
         all_selected_increments = self.selection['increments']
+
         self.pick('increments',all_selected_increments[0:1])
         first = self.list_data()
+
         self.pick('increments',all_selected_increments[-1:])
-        last  = self.list_data()
+        last  = '' if len(all_selected_increments) < 2 else self.list_data()
+
         self.pick('increments',all_selected_increments)
-        in_between = ''.join(['\n{}\n  ...\n'.format(inc) for inc in all_selected_increments[1:-2]])
+
+        in_between = '' if len(all_selected_increments) < 3 else \
+                     ''.join(['\n{}\n  ...\n'.format(inc) for inc in all_selected_increments[1:-2]])
+
         return util.srepr(first + in_between + last)
 
 
@@ -105,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 []
 
@@ -197,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 []
 
@@ -213,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 []
 
@@ -229,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 []
 
@@ -256,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
 
@@ -320,7 +326,7 @@ class Result:
 
         Parameters
         ----------
-            datasets : iterable or str or Boolean
+            datasets : iterable or str or bool
 
         Examples
         --------
@@ -454,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'][()]
@@ -1009,7 +1015,7 @@ class Result:
                 continue
             lock.acquire()
             with h5py.File(self.fname, 'a') as f:
-                try:
+                try:                                                                                # ToDo: Replace if exists?
                     dataset = f[result[0]].create_dataset(result[1]['label'],data=result[1]['data'])
                     for l,v in result[1]['meta'].items():
                         dataset.attrs[l]=v.encode()

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.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
-    return _np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3)
 
 
 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
+    coord0 : numpy.ndarray of shape (:,3)
-        array of undeformed cell coordinates.
+        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 \
+    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])):
+            _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('I_nput data is not a regular grid.')
+        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,
+    return _np.stack(_np.meshgrid(_np.linspace(origin[0],size[0]+origin[0],grid[0]+1),
-                    origin[1]:size[1]+origin[1]:(grid[1]+1)*1j,
+                                  _np.linspace(origin[1],size[1]+origin[1],grid[1]+1),
-                    origin[2]:size[2]+origin[2]:(grid[2]+1)*1j].T
+                                  _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
+    coord0 : numpy.ndarray of shape (:,3)
-        array of undeformed nodal coordinates.
+        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 \
+    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)):
+            _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('I_nput data is not a regular grid.')
+        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/tests/test_Result.py

@@ -24,6 +24,10 @@ def reference_dir(reference_dir_base):
 
 class TestResult:
 
+    def test_self_report(self,default):
+        print(default)
+
+
     def test_time_increments(self,default):
         shape = default.read_dataset(default.get_dataset_location('F'),0).shape
         default.set_by_time(0.0,20.0)

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/homogenization.f90

@@ -184,8 +184,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
     e, &                                                                                            !< element number
     mySource, &
     myNgrains
-  real(pReal), dimension(3,3) :: &
-    subF
   real(pReal), dimension(discretization_nIP,discretization_nElem) :: &
     subFrac, &
     subStep
@@ -376,16 +374,15 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
 
 !--------------------------------------------------------------------------------------------------
 ! deformation partitioning
-! based on materialpoint_subF0,.._subF,crystallite_partionedF0, and homogenization_state,
+      !$OMP PARALLEL DO PRIVATE(myNgrains)
-! results in crystallite_partionedF
-      !$OMP PARALLEL DO PRIVATE(myNgrains,subF)
       elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
         myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
         IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2)
           if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then                                   ! requested but not yet done
-            subF = materialpoint_F0(1:3,1:3,i,e) &
+            call partitionDeformation(materialpoint_F0(1:3,1:3,i,e) &
-                 + (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e))*(subStep(i,e)+subFrac(i,e))
+                                      + (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e))&
-            call partitionDeformation(subF,i,e)                                                     ! partition deformation onto constituents
+                                         *(subStep(i,e)+subFrac(i,e)), &
+                                      i,e)
             crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e)                                       ! propagate materialpoint dt to grains
             crystallite_requested(1:myNgrains,i,e) = .true.                                         ! request calculation for constituents
           else
@@ -397,23 +394,22 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
 
 !--------------------------------------------------------------------------------------------------
 ! crystallite integration
-! based on crystallite_partionedF0,.._partionedF
-! incrementing by crystallite_dt
-
       converged = crystallite_stress() !ToDo: MD not sure if that is the best logic
 
 !--------------------------------------------------------------------------------------------------
 ! state update
-     !$OMP PARALLEL DO PRIVATE(subF)
+     !$OMP PARALLEL DO
       elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
         IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
           if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then
             if (.not. converged(i,e)) then
               doneAndHappy(1:2,i,e) = [.true.,.false.]
             else
-              subF = materialpoint_F0(1:3,1:3,i,e) &
+              doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), &
-                   + (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e))*(subStep(i,e)+subFrac(i,e))
+                                                  materialpoint_F0(1:3,1:3,i,e) &
-              doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e),subF,i,e)
+                                                  + (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e)) &
+                                                     *(subStep(i,e)+subFrac(i,e)), &
+                                                 i,e)
               converged(i,e) = all(doneAndHappy(1:2,i,e))                                           ! converged if done and happy
             endif
           endif

src/quaternions.f90

@@ -118,7 +118,7 @@ end subroutine quaternions_init
 
 
 !---------------------------------------------------------------------------------------------------
-!> construct a quaternion from a 4-vector
+!> @brief construct a quaternion from a 4-vector
 !---------------------------------------------------------------------------------------------------
 type(quaternion) pure function init__(array)
 
@@ -133,7 +133,7 @@ end function init__
 
 
 !---------------------------------------------------------------------------------------------------
-!> assign a quaternion
+!> @brief assign a quaternion
 !---------------------------------------------------------------------------------------------------
 elemental pure subroutine assign_quat__(self,other)
 
@@ -146,7 +146,7 @@ end subroutine assign_quat__
 
 
 !---------------------------------------------------------------------------------------------------
-!> assign a 4-vector
+!> @brief assign a 4-vector
 !---------------------------------------------------------------------------------------------------
 pure subroutine assign_vec__(self,other)
 
@@ -162,7 +162,7 @@ end subroutine assign_vec__
 
 
 !---------------------------------------------------------------------------------------------------
-!> add a quaternion
+!> @brief add a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function add__(self,other)
 
@@ -175,7 +175,7 @@ end function add__
 
 
 !---------------------------------------------------------------------------------------------------
-!> return (unary positive operator)
+!> @brief return (unary positive operator)
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function pos__(self)
 
@@ -187,7 +187,7 @@ end function pos__
 
 
 !---------------------------------------------------------------------------------------------------
-!> subtract a quaternion
+!> @brief subtract a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function sub__(self,other)
 
@@ -200,7 +200,7 @@ end function sub__
 
 
 !---------------------------------------------------------------------------------------------------
-!> negate (unary negative operator)
+!> @brief negate (unary negative operator)
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function neg__(self)
 
@@ -212,7 +212,7 @@ end function neg__
 
 
 !---------------------------------------------------------------------------------------------------
-!> multiply with a quaternion
+!> @brief multiply with a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function mul_quat__(self,other)
 
@@ -227,7 +227,7 @@ end function mul_quat__
 
 
 !---------------------------------------------------------------------------------------------------
-!> multiply with a scalar
+!> @brief multiply with a scalar
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function mul_scal__(self,scal)
 
@@ -240,7 +240,7 @@ end function mul_scal__
 
 
 !---------------------------------------------------------------------------------------------------
-!> divide by a quaternion
+!> @brief divide by a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function div_quat__(self,other)
 
@@ -252,7 +252,7 @@ end function div_quat__
 
 
 !---------------------------------------------------------------------------------------------------
-!> divide by a scalar
+!> @brief divide by a scalar
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function div_scal__(self,scal)
 
@@ -265,7 +265,7 @@ end function div_scal__
 
 
 !---------------------------------------------------------------------------------------------------
-!> test equality
+!> @brief test equality
 !---------------------------------------------------------------------------------------------------
 logical elemental pure function eq__(self,other)
 
@@ -278,7 +278,7 @@ end function eq__
 
 
 !---------------------------------------------------------------------------------------------------
-!> test inequality
+!> @brief test inequality
 !---------------------------------------------------------------------------------------------------
 logical elemental pure function neq__(self,other)
 
@@ -290,7 +290,7 @@ end function neq__
 
 
 !---------------------------------------------------------------------------------------------------
-!> raise to the power of a quaternion
+!> @brief raise to the power of a quaternion
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function pow_quat__(self,expon)
 
@@ -303,7 +303,7 @@ end function pow_quat__
 
 
 !---------------------------------------------------------------------------------------------------
-!> raise to the power of a scalar
+!> @brief raise to the power of a scalar
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function pow_scal__(self,expon)
 
@@ -316,7 +316,7 @@ end function pow_scal__
 
 
 !---------------------------------------------------------------------------------------------------
-!> take exponential
+!> @brief take exponential
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function exp__(a)
 
@@ -336,7 +336,7 @@ end function exp__
 
 
 !---------------------------------------------------------------------------------------------------
-!> take logarithm
+!> @brief take logarithm
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function log__(a)
 
@@ -356,7 +356,7 @@ end function log__
 
 
 !---------------------------------------------------------------------------------------------------
-!> return norm
+!> @brief return norm
 !---------------------------------------------------------------------------------------------------
 real(pReal) elemental pure function abs__(self)
 
@@ -368,7 +368,7 @@ end function abs__
 
 
 !---------------------------------------------------------------------------------------------------
-!> calculate dot product
+!> @brief calculate dot product
 !---------------------------------------------------------------------------------------------------
 real(pReal) elemental pure function dot_product__(a,b)
 
@@ -380,7 +380,7 @@ end function dot_product__
 
 
 !---------------------------------------------------------------------------------------------------
-!> take conjugate complex
+!> @brief take conjugate complex
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function conjg__(self)
 
@@ -392,7 +392,7 @@ end function conjg__
 
 
 !---------------------------------------------------------------------------------------------------
-!> homomorph
+!> @brief homomorph
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function homomorphed(self)
 
@@ -404,7 +404,7 @@ end function homomorphed
 
 
 !---------------------------------------------------------------------------------------------------
-!> return as plain array
+!> @brief return as plain array
 !---------------------------------------------------------------------------------------------------
 pure function asArray(self)
 
@@ -417,7 +417,7 @@ end function asArray
 
 
 !---------------------------------------------------------------------------------------------------
-!> real part (scalar)
+!> @brief real part (scalar)
 !---------------------------------------------------------------------------------------------------
 pure function real__(self)
 
@@ -430,7 +430,7 @@ end function real__
 
 
 !---------------------------------------------------------------------------------------------------
-!> imaginary part (3-vector)
+!> @brief imaginary part (3-vector)
 !---------------------------------------------------------------------------------------------------
 pure function aimag__(self)
 
@@ -443,7 +443,7 @@ end function aimag__
 
 
 !---------------------------------------------------------------------------------------------------
-!> inverse
+!> @brief inverse
 !---------------------------------------------------------------------------------------------------
 type(quaternion) elemental pure function inverse(self)