Merge branch 'development' into Marc2020

.gitlab-ci.yml

@@ -141,13 +141,6 @@ Pre_General:
     - release
 
 ###################################################################################################
-Post_AverageDown:
-  stage: postprocessing
-  script: averageDown/test.py
-  except:
-    - master
-    - release
-
 Post_ASCIItable:
   stage: postprocessing
   script: ASCIItable/test.py

PRIVATE

@@ -1 +1 @@
-Subproject commit 3fc9d58a35614fd8ffa1179e634431eb457d0150
+Subproject commit c6dc3f116fc6007caf076772d5c29f0a0523544d

README

@@ -8,6 +8,6 @@ Max-Planck-Str. 1
 40237 Düsseldorf
 Germany
 
-Email: DAMASK@mpie.de
+damask@mpie.de
 https://damask.mpie.de
 https://magit1.mpie.de

VERSION

@@ -1 +1 @@
-v2.0.3-2945-gb7e03364
+v3.0.0-alpha-89-g23d61511

examples/SpectralMethod/Polycrystal/material.config

@@ -1,126 +0,0 @@
-#-------------------#
-<homogenization>
-#-------------------#
-
-[SX]
-mech	none
-
-#-------------------#
-<phase>
-#-------------------#
-[Aluminum_phenopowerlaw]
-elasticity              hooke
-plasticity              phenopowerlaw
-
-(output) resistance_slip
-(output) orientation    # quaternion
-(output) F              # deformation gradient tensor
-(output) Fe             # elastic deformation gradient tensor
-(output) Fp             # plastic deformation gradient tensor
-(output) P              # first Piola-Kichhoff stress tensor
-(output) Lp             # plastic velocity gradient tensor
-
-
-lattice_structure       fcc
-Nslip                   12   # per family
-
-c11                     106.75e9
-c12                     60.41e9
-c44                     28.34e9
-
-gdot0_slip              0.001
-n_slip                  20
-tau0_slip                 31e6 # per family
-tausat_slip               63e6 # per family
-a_slip                  2.25
-h0_slipslip             75e6
-interaction_slipslip    1 1 1.4 1.4 1.4 1.4
-
-#-------------------#
-<microstructure>
-#-------------------#
-[Grain01]
-(constituent)  phase 1   texture 01   fraction 1.0
-[Grain02]
-(constituent)  phase 1   texture 02   fraction 1.0
-[Grain03]
-(constituent)  phase 1   texture 03   fraction 1.0
-[Grain04]
-(constituent)  phase 1   texture 04   fraction 1.0
-[Grain05]
-(constituent)  phase 1   texture 05   fraction 1.0
-[Grain06]
-(constituent)  phase 1   texture 06   fraction 1.0
-[Grain07]
-(constituent)  phase 1   texture 07   fraction 1.0
-[Grain08]
-(constituent)  phase 1   texture 08   fraction 1.0
-[Grain09]
-(constituent)  phase 1   texture 09   fraction 1.0
-[Grain10]
-(constituent)  phase 1   texture 10   fraction 1.0
-[Grain11]
-(constituent)  phase 1   texture 11   fraction 1.0
-[Grain12]
-(constituent)  phase 1   texture 12   fraction 1.0
-[Grain13]
-(constituent)  phase 1   texture 13   fraction 1.0
-[Grain14]
-(constituent)  phase 1   texture 14   fraction 1.0
-[Grain15]
-(constituent)  phase 1   texture 15   fraction 1.0
-[Grain16]
-(constituent)  phase 1   texture 16   fraction 1.0
-[Grain17]
-(constituent)  phase 1   texture 17   fraction 1.0
-[Grain18]
-(constituent)  phase 1   texture 18   fraction 1.0
-[Grain19]
-(constituent)  phase 1   texture 19   fraction 1.0
-[Grain20]
-(constituent)  phase 1   texture 20   fraction 1.0
-
-
-#-------------------#
-<texture>
-#-------------------#
-[Grain01]
-(gauss)  phi1 0.0    Phi 0.0    phi2 0.0
-[Grain02]
-(gauss)  phi1 257.468172    Phi 53.250534    phi2 157.331503
-[Grain03]
-(gauss)  phi1 216.994815    Phi 94.418518    phi2 251.147231
-[Grain04]
-(gauss)  phi1 196.157946    Phi 55.870978    phi2 21.68117
-[Grain05]
-(gauss)  phi1 152.515728    Phi 139.769395    phi2 240.036018
-[Grain06]
-(gauss)  phi1 232.521881    Phi 73.749222    phi2 241.429633
-[Grain07]
-(gauss)  phi1 157.531396    Phi 135.503513    phi2 75.737722
-[Grain08]
-(gauss)  phi1 321.03828    Phi 27.209843    phi2 46.413467
-[Grain09]
-(gauss)  phi1 346.918594    Phi 87.495569    phi2 113.554206
-[Grain10]
-(gauss)  phi1 138.038947    Phi 99.827132    phi2 130.935878
-[Grain11]
-(gauss)  phi1 285.021014    Phi 118.092004    phi2 205.270837
-[Grain12]
-(gauss)  phi1 190.402171    Phi 56.738068    phi2 157.896545
-[Grain13]
-(gauss)  phi1 204.496042    Phi 95.031265    phi2 355.814582
-[Grain14]
-(gauss)  phi1 333.21479    Phi 82.133355    phi2 36.736132
-[Grain15]
-(gauss)  phi1 25.572981    Phi 164.242648    phi2 75.195632
-[Grain16]
-(gauss)  phi1 31.366548    Phi 76.392403    phi2 58.071426
-[Grain17]
-(gauss)  phi1 7.278623    Phi 77.044663    phi2 235.118997
-[Grain18]
-(gauss)  phi1 299.743144    Phi 76.475096    phi2 91.184977
-[Grain19]
-(gauss)  phi1 280.13643    Phi 27.439718    phi2 167.871878
-[Grain20]
-(gauss)  phi1 313.204373    Phi 68.676053    phi2 87.993213

examples/SpectralMethod/Polycrystal/material.yaml

@@ -0,0 +1,123 @@
+homogenization:
+  SX:
+    mech: {type: none}
+microstructure:
+- constituents:
+  - fraction: 1.0
+    orientation: [1.0, 0.0, 0.0, 0.0]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.33012772942625784, -0.6781865350268957, 0.6494525351030648, 0.09638521992649676]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.43596817439583935, -0.5982537129781701, 0.046599032277502436, 0.6707106499919265]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.169734823419553, -0.699615227367322, -0.6059581215838098, -0.33844257746495854]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.9698864809294915, 0.1729052643205874, -0.15948307917616958, 0.06315956884687175]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.46205660912967883, 0.3105054068891252, -0.617849551030653, 0.555294529545738]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.4512443497461787, -0.7636045534540555, -0.04739348426715133, -0.45939142396805815]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.2161856212656443, -0.6581450184826598, -0.5498086209601588, 0.4667112513346289]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.8753220715350803, -0.4561599367657419, -0.13298279533852678, -0.08969369719975541]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.11908260752431069, 0.18266024809834172, -0.7144822594012615, -0.664807992845101]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.751104669484278, 0.5585633382623958, -0.34579336397009175, 0.06538900566860861]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.08740438971703973, 0.8991264096610437, -0.4156704205935976, 0.10559485570696363]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.5584325870096193, 0.6016408353068798, -0.14280340445801173, 0.5529814994483859]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.4052725440888093, 0.25253073423599154, 0.5693263597910454, -0.669215876471182]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.7570164606888676, 0.15265448024694664, -0.5998021466848317, 0.20942796551297105]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.6987659297138081, -0.132172211261028, -0.19693254724422338, 0.6748883269678543]
+    phase: Aluminum
+  homogenization: SX
+- constituents:
+  - fraction: 1.0
+    orientation: [0.7729330445886478, 0.21682179052722322, -0.5207379472917645, 0.2905078484066341]
+    phase: Aluminum
+  homogenization: SX
+phase:
+  Aluminum:
+    elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}
+    generic:
+      output: [F, P, Fe, Fp, Lp]
+    lattice: fcc
+    plasticity:
+      N_sl: [12]
+      a_sl: 2.25
+      atol_xi: 1.0
+      dot_gamma_0_sl: 0.001
+      h_0_sl_sl: 75e6
+      h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
+      n_sl: 20
+      output: [xi_sl]
+      xi_0_sl: [31e6]
+      xi_inf_sl: [63e6]
+      type: phenopowerlaw
+

processing/legacy/addAPS34IDEstrainCoords.py

@@ -49,4 +49,4 @@ for name in filenames:
 
     table.add('coord',rot_to_TSL.broadcast_to(coord.shape[0]) @ coord,scriptID+' '+' '.join(sys.argv[1:]))
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

processing/legacy/addCumulative.py

@@ -45,4 +45,4 @@ for name in filenames:
                   np.cumprod(table.get(label),0) if options.product else np.cumsum(table.get(label),0),
                   scriptID+' '+' '.join(sys.argv[1:]))
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

processing/legacy/imageData.py

@@ -1,184 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from optparse import OptionParser
-
-import numpy as np
-from PIL import Image
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
-Generate PNG image from data in given column (or 2D data of overall table).
-
-""", version = scriptID)
-
-parser.add_option('-l','--label',
-                  dest = 'label',
-                  type = 'string', metavar = 'string',
-                  help = 'column containing data [all]')
-parser.add_option('-r','--range',
-                  dest = 'range',
-                  type = 'float', nargs = 2, metavar = 'float float',
-                  help = 'data range (min max) [auto]')
-parser.add_option('--gap', '--transparent',
-                  dest = 'gap',
-                  type = 'float', metavar = 'float',
-                  help = 'value to treat as transparent [%default]')
-parser.add_option('-d','--dimension',
-                  dest = 'dimension',
-                  type = 'int', nargs = 2, metavar = 'int int',
-                  help = 'data dimension (width height) [native]')
-parser.add_option('--color',
-                  dest = 'color',
-                  type = 'string', metavar = 'string',
-                  help = 'color scheme [%default]')
-parser.add_option('--invert',
-                  dest = 'invert',
-                  action = 'store_true',
-                  help = 'invert color scheme')
-parser.add_option('--abs',
-                  dest = 'abs',
-                  action = 'store_true',
-                  help = 'magnitude of values')
-parser.add_option('--log',
-                  dest = 'log',
-                  action = 'store_true',
-                  help = 'log10 of values')
-parser.add_option('--fliplr',
-                  dest = 'flipLR',
-                  action = 'store_true',
-                  help = 'flip around vertical axis')
-parser.add_option('--flipud',
-                  dest = 'flipUD',
-                  action = 'store_true',
-                  help = 'flip around horizontal axis')
-parser.add_option('--crop',
-                  dest = 'crop',
-                  type = 'int', nargs = 4, metavar = 'int int int int',
-                  help = 'pixels cropped on left, right, top, bottom')
-parser.add_option('-N','--pixelsize',
-                  dest = 'pixelsize',
-                  type = 'int', metavar = 'int',
-                  help = 'pixel per data point')
-parser.add_option('-x','--pixelsizex',
-                  dest = 'pixelsizex',
-                  type = 'int', metavar = 'int',
-                  help = 'pixel per data point along x')
-parser.add_option('-y','--pixelsizey',
-                  dest = 'pixelsizey',
-                  type = 'int', metavar = 'int',
-                  help = 'pixel per data point along y')
-parser.add_option('--show',
-                  dest = 'show',
-                  action = 'store_true',
-                  help = 'show resulting image')
-
-parser.set_defaults(label = None,
-                    range = [0.0,0.0],
-                    gap = None,
-                    dimension = [],
-                    abs = False,
-                    log = False,
-                    flipLR = False,
-                    flipUD = False,
-                    color = "gray",
-                    invert = False,
-                    crop = [0,0,0,0],
-                    pixelsize  = 1,
-                    pixelsizex = 1,
-                    pixelsizey = 1,
-                    show = False,
-                   )
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if options.pixelsize > 1: (options.pixelsizex,options.pixelsizey) = [options.pixelsize]*2
-
-# --- color palette ---------------------------------------------------------------------------------
-
-theMap = damask.Colormap(predefined = options.color)
-if options.invert: theMap = theMap.invert()
-theColors = np.uint8(np.array(theMap.export(format = 'list',steps = 256))*255)
-
-# --- loop over input files -------------------------------------------------------------------------
-for name in filenames:
-  try:
-    table = damask.ASCIItable(name = name, labeled = options.label is not None, readonly = True)
-  except IOError:
-    continue
-  damask.util.report(scriptName,name)
-
-# ------------------------------------------ read header ------------------------------------------
-
-  table.head_read()
-
-# ------------------------------------------ process data ------------------------------------------
-
-  missing_labels = table.data_readArray(options.label)
-  if len(missing_labels) > 0:
-    damask.util.croak('column {} not found.'.format(options.label))
-    table.close(dismiss = True)                                                                     # close ASCIItable and remove empty file
-    continue
-# convert data to values between 0 and 1 and arrange according to given options
-  if options.dimension != []: table.data = table.data.reshape(options.dimension[1],options.dimension[0])
-  if options.abs:             table.data = np.abs(table.data)
-  if options.log:             table.data = np.log10(table.data);options.range = np.log10(options.range)
-  if options.flipLR:          table.data = np.fliplr(table.data)
-  if options.flipUD:          table.data = np.flipud(table.data)
-
-  mask = np.logical_or(table.data == options.gap, np.isnan(table.data))\
-                                          if options.gap else np.logical_not(np.isnan(table.data))  # mask gap and NaN (if gap present)
-  if np.all(np.array(options.range) == 0.0):
-    options.range = [table.data[mask].min(),
-                     table.data[mask].max()]
-    damask.util.croak('data range: {0} – {1}'.format(*options.range))
-
-  delta =      max(options.range) - min(options.range)
-  avg   = 0.5*(max(options.range) + min(options.range))
-
-  if delta * 1e8 <= avg:                                                                           # delta around numerical noise
-    options.range = [min(options.range) - 0.5*avg, max(options.range) + 0.5*avg]                   # extend range to have actual data centered within
-
-  table.data =         (table.data - min(options.range)) / \
-               (max(options.range) - min(options.range))
-  
-  table.data = np.clip(table.data,0.0,1.0).\
-                  repeat(options.pixelsizex,axis = 1).\
-                  repeat(options.pixelsizey,axis = 0)
-
-  mask =       mask.\
-                  repeat(options.pixelsizex,axis = 1).\
-                  repeat(options.pixelsizey,axis = 0)
-
-  (height,width) = table.data.shape
-  damask.util.croak('image dimension: {0} x {1}'.format(width,height))
-
-  im = Image.fromarray(np.dstack((theColors[np.array(255*table.data,dtype = np.uint8)],
-                                  255*mask.astype(np.uint8))), 'RGBA').\
-             crop((       options.crop[0],
-                          options.crop[2],
-                   width -options.crop[1],
-                   height-options.crop[3]))
-
-# ------------------------------------------ output result -----------------------------------------
-
-  im.save(sys.stdout if not name else
-          os.path.splitext(name)[0]+ \
-          ('' if options.label is None else '_'+options.label)+ \
-          '.png',
-          format = "PNG")
-
-  table.close()                                                                                     # close ASCII table
-  if options.show: im.show()

processing/legacy/imageDataRGB.py

@@ -1,134 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from optparse import OptionParser
-
-import numpy as np
-from PIL import Image
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
-Generate PNG image from data in given column vector containing RGB tuples.
-
-""", version = scriptID)
-
-parser.add_option('-l','--label',
-                 dest = 'label',
-                 type = 'string', metavar = 'string',
-                 help = 'column containing RGB triplet')
-parser.add_option('-d','--dimension',
-                  dest = 'dimension',
-                  type = 'int', nargs = 2, metavar = 'int int',
-                  help = 'data dimension (width height)')
-parser.add_option('--fliplr',
-                  dest = 'flipLR',
-                  action = 'store_true',
-                  help = 'flip around vertical axis')
-parser.add_option('--flipud',
-                  dest = 'flipUD',
-                  action = 'store_true',
-                  help = 'flip around horizontal axis')
-parser.add_option('--crop',
-                  dest = 'crop',
-                  type = 'int', nargs = 4, metavar = ' '.join(['int']*4),
-                  help = 'pixels cropped on left, right, top, bottom')
-parser.add_option('-N','--pixelsize',
-                  dest = 'pixelsize',
-                  type = 'int', metavar = 'int',
-                  help = 'pixels per data point')
-parser.add_option('-x','--pixelsizex',
-                  dest = 'pixelsizex',
-                  type = 'int', metavar = 'int',
-                  help = 'pixels per data point along x')
-parser.add_option('-y','--pixelsizey',
-                  dest = 'pixelsizey',
-                  type = 'int', metavar = 'int',
-                  help = 'pixels per data point along y')
-parser.add_option('--show',
-                  dest = 'show',
-                  action = 'store_true',
-                  help = 'show resulting image')
-
-parser.set_defaults(label = None,
-                    dimension = [],
-                    flipLR = False,
-                    flipUD = False,
-                    crop = [0,0,0,0],
-                    pixelsize  = 1,
-                    pixelsizex = 1,
-                    pixelsizey = 1,
-                    show = False,
-                   )
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if options.dimension == []: parser.error('dimension of data array missing')
-if options.pixelsize > 1: (options.pixelsizex,options.pixelsizey) = [options.pixelsize]*2
-
-# --- loop over input files -------------------------------------------------------------------------
-for name in filenames:
-  try:
-    table = damask.ASCIItable(name = name, labeled = options.label is not None, readonly = True)
-  except IOError:
-    continue
-  damask.util.report(scriptName,name)
-
-# ------------------------------------------ read header ------------------------------------------
-
-  table.head_read()
-
-# ------------------------------------------ process data ------------------------------------------
-
-  errors = []
-  
-  missing_labels = table.data_readArray(options.label)
-  if len(missing_labels) > 0:
-    errors.append('column{} {} not found'.format('s' if len(missing_labels) > 1 else '',
-                                                 ', '.join(missing_labels)))
-  if table.label_dimension(options.label) != 3:
-    errors.append('column {} does not have dimension'.format(options.label))
-
-  if errors != []:
-    damask.util.croak(errors)
-    table.close(dismiss = True)                                                                     # close ASCII table file handles and delete output file
-    continue
-# convert data to shape and arrange according to given options
-  if options.dimension != []: table.data = table.data.reshape(options.dimension[1],options.dimension[0],3)
-  if options.flipLR:          table.data = np.fliplr(table.data)
-  if options.flipUD:          table.data = np.flipud(table.data)
-
-  table.data = table.data.repeat(options.pixelsizex,axis=1).\
-                          repeat(options.pixelsizey,axis=0)
-
-  table.data *= 1. if np.any(table.data > 1.0) else 255.0                                          # ensure 8 bit data range
-
-  (height,width,bands) = table.data.shape
-  damask.util.croak('image dimension: {0} x {1}'.format(width,height))
-
-  im = Image.fromarray(table.data.astype('uint8'), 'RGB').\
-             crop((       options.crop[0],
-                          options.crop[2],
-                   width -options.crop[1],
-                   height-options.crop[3]))
-
-# ------------------------------------------ output result -----------------------------------------
-
-  im.save(os.path.splitext(name)[0]+ \
-          ('_'+options.label if options.label else '')+ \
-          '.png' if name else sys.stdout,
-          format = "PNG")
-
-  table.close()                                                                                     # close ASCII table
-  if options.show: im.show()

processing/post/DADF5_postResults.py

@@ -60,4 +60,4 @@ for filename in options.filenames:
             os.mkdir(dirname,0o755)
         file_out = '{}_inc{}.txt'.format(os.path.splitext(os.path.split(filename)[-1])[0],
                                          inc[3:].zfill(N_digits))
-        table.to_ASCII(os.path.join(dirname,file_out))
+        table.to_file(os.path.join(dirname,file_out))

processing/post/addCompatibilityMismatch.py

@@ -191,4 +191,4 @@ for name in filenames:
               volumeMismatch.reshape(-1,1,order='F'),
               scriptID+' '+' '.join(sys.argv[1:]))
 
-  table.to_ASCII(sys.stdout if name is None else name)
+  table.to_file(sys.stdout if name is None else name)

processing/post/addCurl.py

@@ -55,4 +55,4 @@ for name in filenames:
                   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)
+    table.to_file(sys.stdout if name is None else name)

processing/post/addDerivative.py

@@ -71,4 +71,4 @@ for name in filenames:
                   derivative(table.get(options.coordinates),table.get(label)),
                   scriptID+' '+' '.join(sys.argv[1:]))
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

processing/post/addDisplacement.py

@@ -60,7 +60,7 @@ for name in filenames:
         table.add('fluct({}).{}'.format(options.f,options.pos),
                   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')
+        table.to_file(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,F).reshape(-1,3,order='F'),
@@ -68,4 +68,4 @@ for name in filenames:
         table.add('fluct({}).{}'.format(options.f,options.pos),
                   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)
+        table.to_file(sys.stdout if name is None else name)

processing/post/addDivergence.py

@@ -55,4 +55,4 @@ for name in filenames:
                   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)
+    table.to_file(sys.stdout if name is None else name)

processing/post/addEuclideanDistance.py

@@ -184,4 +184,4 @@ for name in filenames:
                   distance[i,:],
                   scriptID+' '+' '.join(sys.argv[1:]))
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

processing/post/addGaussian.py

@@ -73,4 +73,4 @@ for name in filenames:
                                                   mode = 'wrap' if options.periodic else 'nearest'),
                   scriptID+' '+' '.join(sys.argv[1:]))
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

processing/post/addGradient.py

@@ -55,4 +55,4 @@ for name in filenames:
                   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)
+    table.to_file(sys.stdout if name is None else name)

processing/post/addIndexed.py

@@ -57,8 +57,8 @@ for name in filenames:
     table        = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
     indexedTable = damask.Table.from_ASCII(options.asciitable)
     idx = np.reshape(table.get(options.index).astype(int) + options.offset,(-1))-1
-    
+
     for data in options.label:
         table.add(data+'_addIndexed',indexedTable.get(data)[idx],scriptID+' '+' '.join(sys.argv[1:]))
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

processing/post/addOrientations.py

@@ -147,4 +147,4 @@ for name in filenames:
     if 'axisangle' in options.output:
         table.add('om({})'.format(label),o.as_axisangle(options.degrees), scriptID+' '+' '.join(sys.argv[1:]))
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

processing/post/addSchmidfactors.py

@@ -189,4 +189,4 @@ for name in filenames:
     for i,label in enumerate(labels):
         table.add(label,S[:,i],scriptID+' '+' '.join(sys.argv[1:]))
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

processing/post/averageDown.py

@@ -1,98 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import numpy as np
-import scipy.ndimage
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Average each data block of size 'packing' into single values thus reducing the former grid to grid/packing.
-
-""", version = scriptID)
-
-parser.add_option('-c','--coordinates',
-                  dest = 'pos',
-                  type = 'string', metavar = 'string',
-                  help = 'column label of coordinates [%default]')
-parser.add_option('-p','--packing',
-                  dest = 'packing',
-                  type = 'int', nargs = 3, metavar = 'int int int',
-                  help = 'size of packed group [%default]')
-parser.add_option('--shift',
-                  dest = 'shift',
-                  type = 'int', nargs = 3, metavar = 'int int int',
-                  help = 'shift vector of packing stencil [%default]')
-parser.add_option('-g', '--grid',
-                  dest = 'grid',
-                  type = 'int', nargs = 3, metavar = 'int int int',
-                  help = 'grid in x,y,z (optional)')
-parser.add_option('-s', '--size',
-                  dest = 'size',
-                  type = 'float', nargs = 3, metavar = 'float float float',
-                  help = 'size in x,y,z (optional)')
-parser.set_defaults(pos     = 'pos',
-                    packing = (2,2,2),
-                    shift   = (0,0,0),
-                   )
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-packing = np.array(options.packing,dtype = int)
-shift   = np.array(options.shift,  dtype = int)
-
-prefix = 'averagedDown{}x{}x{}_'.format(*packing)
-if any(shift != 0): prefix += 'shift{:+}{:+}{:+}_'.format(*shift)
-
-
-for name in filenames:
-  damask.util.report(scriptName,name)
-
-  table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-
-  if (options.grid is None or options.size is None):
-    grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
-  else:
-    grid   = np.array(options.grid,'i')
-    size   = np.array(options.size,'d')
-
-  packing = np.where(grid == 1,1,packing)                                                           # reset packing to 1 where grid==1
-  shift   = np.where(grid == 1,0,shift)                                                             # reset   shift to 0 where grid==1
-  packedGrid = np.maximum(np.ones(3,'i'),grid//packing)
-
-  data = table.data.values.reshape(tuple(grid)+(-1,),order = 'F')
-  averagedDown = scipy.ndimage.filters.uniform_filter( \
-                  np.roll(
-                  np.roll(
-                  np.roll(data,
-                          -shift[0],axis = 0),
-                          -shift[1],axis = 1),
-                          -shift[2],axis = 2),
-                  size = list(packing) + [1],
-                  mode = 'wrap',
-                  origin = list(-(packing//2)) + [0])\
-                  [::packing[0],::packing[1],::packing[2],:].reshape((packedGrid.prod(),-1),order = 'F')
-
-
-  table = damask.Table(averagedDown,table.shapes,table.comments)
-
-  coords = damask.grid_filters.cell_coord0(packedGrid,size,shift/packedGrid*size+origin)
-  table.set(options.pos, coords.reshape(-1,3,order='F'))
-
-
-  outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name))
-  table.to_ASCII(sys.stdout if name is None else outname)

processing/post/blowUp.py

@@ -1,72 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-from scipy import ndimage
-import numpy as np
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Blows up each value to a surrounding data block of size 'packing' thus increasing the former resolution
-to resolution*packing.
-
-""", version = scriptID)
-
-parser.add_option('-c','--coordinates',
-                  dest = 'pos', metavar = 'string',
-                  help = 'column label of coordinates [%default]')
-parser.add_option('-p','--packing',
-                  dest = 'packing', type = 'int', nargs = 3, metavar = 'int int int',
-                  help = 'dimension of packed group [%default]')
-parser.add_option('-g','--grid',
-                  dest = 'resolution', type = 'int', nargs = 3, metavar = 'int int int',
-                  help = 'grid in x,y,z (optional)')
-parser.add_option('-s','--size',
-                  dest = 'dimension', type = 'float', nargs = 3, metavar = 'int int int',
-                  help = 'size in x,y,z (optional)')
-parser.set_defaults(pos  = 'pos',
-                    packing = (2,2,2),
-                    grid    = (0,0,0),
-                    size    = (0.0,0.0,0.0),
-                   )
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-options.packing = np.array(options.packing)
-prefix = 'blowUp{}x{}x{}_'.format(*options.packing)
-
-
-for name in filenames:
-  damask.util.report(scriptName,name)
-  
-  table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-  grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos))
-  
-  packing = np.array(options.packing,'i')
-  outSize = grid*packing
-  
-  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,order='F'))
-  table.set('elem',np.arange(1,outSize.prod()+1))
-  
-  outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name))
-  table.to_ASCII(sys.stdout if name is None else outname)

processing/post/permuteData.py

@@ -58,4 +58,4 @@ for name in filenames:
         rng.shuffle(uniques)
         table.set(label,uniques[inverse], scriptID+' '+' '.join(sys.argv[1:]))
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

processing/post/vtk2ang.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python3
 
 import os
 import sys
@@ -11,54 +11,50 @@ import vtk
 
 import damask
 
- 
+
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
 scriptID   = ' '.join([scriptName,damask.version])
 
 # -----------------------------
 def getHeader(filename,sizeFastIndex,sizeSlowIndex,stepsize):
-  """Returns header for ang file step size in micrometer"""
+  """Returns header for ang file step size in micrometer."""
-  return '\n'.join([ \
+  return '\n'.join([
-    '# TEM_PIXperUM          1.000000', \
+    '# TEM_PIXperUM          1.000000',
-    '# x-star                1.000000', \
+    '# x-star                1.000000',
-    '# y-star                1.000000', \
+    '# y-star                1.000000',
-    '# z-star                1.000000', \
+    '# z-star                1.000000',
-    '# WorkingDistance       18.000000', \
+    '# WorkingDistance       18.000000',
-    '#', \
+    '#',
-    '# Phase                 1', \
+    '# Phase                 1',
-    '# MaterialName          XX', \
+    '# MaterialName          XX',
-    '# Formula               XX', \
+    '# Formula               XX',
-    '# Info', \
+    '# Info',
-    '# Symmetry              43', \
+    '# Symmetry              43',
-    '# LatticeConstants      2.870 2.870 2.870  90.000  90.000  90.000', \
+    '# LatticeConstants      2.870 2.870 2.870  90.000  90.000  90.000',
-    '# NumberFamilies        1', \
+    '# NumberFamilies        1',
-    '# hklFamilies           1  1  0 1 0.000000 1', \
+    '# hklFamilies           1  1  0 1 0.000000 1',
-    '# Categories            0 0 0 0 0 ', \
+    '# Categories            0 0 0 0 0 ',
-    '#', \
+    '#',
-    '# GRID: SqrGrid', \
+    '# GRID: SqrGrid',
-    '# XSTEP: ' + str(stepsize*1e6), \
+    '# XSTEP: ' + str(stepsize*1e6),
-    '# YSTEP: ' + str(stepsize*1e6), \
+    '# YSTEP: ' + str(stepsize*1e6),
-    '# NCOLS_ODD: ' + str(sizeFastIndex), \
+    '# NCOLS_ODD: ' + str(sizeFastIndex),
-    '# NCOLS_EVEN: ' + str(sizeFastIndex), \
+    '# NCOLS_EVEN: ' + str(sizeFastIndex),
-    '# NROWS: ' + str(sizeSlowIndex), \
+    '# NROWS: ' + str(sizeSlowIndex),
-    '#', \
+    '#',
-    '# OPERATOR: ' + string.replace('$Id$','\n','\\n'), \
+    '# OPERATOR: ' + string.replace('$Id$','\n','\\n'),
-    '#', \
+    '#',
-    '# SAMPLEID: %s'%filename, \
+    '# SAMPLEID: {}'.format(filename),
-    '#', \
+    '#',
-    '# SCANID: ', \
+    '# SCANID: ',
-    '#', \
+    '#',
     ]) + '\n'
 
 
 # -----------------------------
 def positiveRadians(angle):
-  """Returns positive angle in radians from angle in degrees"""
+  """Returns positive angle in radians from angle in degrees."""
-  angle = math.radians(float(angle))
+  return math.radians(float(angle)) % (2.*math.pi)
-  while angle < 0.0:
-    angle += 2.0 * math.pi
-
-  return angle
 
 
 # -----------------------------
@@ -67,7 +63,7 @@ def getDataLine(angles,x,y,validData=True):
   Returns string of one line in ang file.
 
   Convention in ang file: y coordinate comes first and is fastest index
-  positions in micrometer
+  positions in micrometer.
   """
   info = {True:  (9999.9, 1.0, 0,99999,0.0),
           False: (  -1.0,-1.0,-1,   -1,1.0)}
@@ -75,18 +71,16 @@ def getDataLine(angles,x,y,validData=True):
                 %(tuple(map(positiveRadians,angles))+(y*1e6,x*1e6)+info[validData])
 
 
-
 # --------------------------------------------------------------------
 # MAIN FUNCTION STARTS HERE
 # --------------------------------------------------------------------
-
 parser = OptionParser(usage='%prog options [file[s]]', description = """
 Builds a ang files from a vtk file.
 
 """, version = scriptID)
 
 
-parser.add_option('--disp','--displacement',dest='dispLabel', 
+parser.add_option('--disp','--displacement',dest='dispLabel',
                                             metavar ='string',
                                             help='label of displacements [%default]')
 parser.add_option('--euler',                dest='eulerLabel', nargs=3,
@@ -110,11 +104,6 @@ parser.add_option('-s','--scale',           dest='scale', type='float',
 parser.add_option('-r','--resolution',      dest='resolution', type='float',
                                             metavar ='float',
                                             help='scaling factor for resolution [%default]')
-parser.add_option('--hex','--hexagonal',    dest='hexagonal', action='store_true',
-                                            help='use in plane hexagonal grid')
-parser.add_option('--interpolation',        dest='interpolation', type='int',
-                                            metavar='float',
-                                            help='number of points for linear interpolation [%default]')
 parser.add_option('--verbose',              dest='verbose', action='store_true',
                                             help='verbose mode')
 parser.add_option('--visualize',            dest='visualize', action='store_true',
@@ -122,7 +111,6 @@ parser.add_option('--visualize',            dest='visualize', action='store_true
 
 parser.set_defaults(dispLabel = 'displacement')
 parser.set_defaults(eulerLabel = ['1_1_eulerangles','1_2_eulerangles','1_3_eulerangles'])
-parser.set_defaults(hexagonal = False)
 parser.set_defaults(normal = [0.0,0.0,-1.0])
 parser.set_defaults(up = [0.0,1.0,0.0])
 parser.set_defaults(Nslices = 1)
@@ -130,7 +118,6 @@ parser.set_defaults(distance = 0.0)
 parser.set_defaults(scale = 1.0)
 parser.set_defaults(resolution = 1.0)
 parser.set_defaults(dispScaling = 1.0)
-parser.set_defaults(interpolation = 1)
 parser.set_defaults(verbose = False)
 parser.set_defaults(visualize = False)
 (options,filenames) = parser.parse_args()
@@ -153,35 +140,26 @@ if np.dot(np.array(options.normal),np.array(options.up)) > 1e-3:
   parser.error('normal vector and up vector have to be orthogonal')
 
 
-# check for options that are not yet implemented
-
-if options.interpolation > 1:
-  parser.error('interpolation not yet supported')
-if options.hexagonal:
-  parser.error('hexagonal grid not yet supported')
-
-
-
 #--- ITERATE OVER FILES AND PROCESS THEM
 
 for filename in filenames:
-  
+
   if options.verbose: sys.stdout.write("\nREADING VTK FILE\n")
 # Read the source file
   reader = vtk.vtkUnstructuredGridReader()
   reader.SetFileName(filename)
   reader.ReadAllScalarsOn()
   reader.ReadAllVectorsOn()
-  reader.Update() 
+  reader.Update()
   undeformedMesh = reader.GetOutput()
 
-  
+
 # Get euler angles from cell data
-  
+
   if options.verbose: sys.stdout.write("\nGETTING EULER ANGLES\n")
   angles = {}
-  for i in range(reader.GetNumberOfScalarsInFile()):
+  for i in range(undeformedMesh.GetPointData().GetNumberOfArrays()):
-    scalarName = reader.GetScalarsNameInFile(i)
+    scalarName = undeformedMesh.GetPointData().GetArrayName(i)
     if scalarName in options.eulerLabel:
       angles[scalarName] = undeformedMesh.GetCellData().GetScalars(scalarName)
       if options.verbose: sys.stdout.write("  found scalar with name %s\n"%scalarName)
@@ -189,14 +167,14 @@ for filename in filenames:
     for label in options.eulerLabel:
       if label not in angles.keys():
         parser.error('Could not find scalar data with name %s'%label)
-  
+
-  
+
 # Get deformed mesh
-  
+
   if options.verbose: sys.stdout.write("\nDEFORM MESH\n")
   warpVector = vtk.vtkWarpVector()
   undeformedMesh.GetPointData().SetActiveVectors(options.dispLabel)
-  warpVector.SetInput(undeformedMesh)
+  warpVector.SetInputData(undeformedMesh)
   warpVector.Update()
   deformedMesh = warpVector.GetOutput()
   box = deformedMesh.GetBounds()              # bounding box in mesh system
@@ -208,24 +186,24 @@ for filename in filenames:
 
 
 # Get cell centers of deformed mesh (position of ips)
-  
+
   if options.verbose: sys.stdout.write("\nGETTING CELL CENTERS OF DEFORMED MESH\n")
   cellCenter = vtk.vtkCellCenters()
   cellCenter.SetVertexCells(0)   # do not generate vertex cells, just points
-  cellCenter.SetInput(deformedMesh)
+  cellCenter.SetInputData(deformedMesh)
   cellCenter.Update()
   meshIPs = cellCenter.GetOutput()
 
 
 # Get outer surface of deformed mesh
-  
+
   if options.verbose: sys.stdout.write("\nGETTING OUTER SURFACE OF DEFORMED MESH\n")
   surfaceFilter = vtk.vtkDataSetSurfaceFilter()
-  surfaceFilter.SetInput(deformedMesh)
+  surfaceFilter.SetInputData(deformedMesh)
   surfaceFilter.Update()
   surface = surfaceFilter.GetOutput()
-  
+
-  
+
 # Get coordinate system for ang files
 # z-vector is normal to slices
 # x-vector corresponds to the up-direction
@@ -248,11 +226,11 @@ for filename in filenames:
 # Get bounding box in rotated system (x,y,z)
 
   if options.verbose: sys.stdout.write("\nGETTING BOUNDING BOX IN ROTATED SYSTEM\n")
-  rotatedbox = [[np.inf,-np.inf] for i in range(3)]     # bounding box in rotated TSL system
+  rotatedbox = [[np.inf,-np.inf] for i in range(3)]        # bounding box in rotated TSL system
-  for n in range(8):                                          # loop over eight vertices of mesh bounding box 
+  for n in range(8):                                       # loop over eight vertices of mesh bounding box
-    vert = np.array([box[0+(n/1)%2], 
+    vert = np.array([box[0+(n//1)%2],
-                        box[2+(n/2)%2], 
+                     box[2+(n//2)%2],
-                        box[4+(n/4)%2]])                      # vertex in mesh system
+                     box[4+(n//4)%2]])                     # vertex in mesh system
     rotatedvert = np.dot(R,vert)                           # vertex in rotated system
     for i in range(3):
       rotatedbox[i][0] = min(rotatedbox[i][0],rotatedvert[i])
@@ -274,7 +252,7 @@ for filename in filenames:
   for i in range(2):
     Npoints.extend([int(math.ceil(extent[i] / options.resolution))])
     correction.extend([float(Npoints[i]) * options.resolution - extent[i]])
-  if options.distance > 0.0: 
+  if options.distance > 0.0:
     Npoints.extend([int(math.ceil(extent[2] / options.distance))])
     correction.extend([float(Npoints[2]) * options.distance - extent[2]])
   else:
@@ -282,8 +260,8 @@ for filename in filenames:
     correction.extend([0.0])
     options.distance = extent[2] / float(options.Nslices)
   for i in range(3):
-    rotatedbox[i][0] = rotatedbox[i][0] - 0.5 * correction[i]
+    rotatedbox[i][0] -= 0.5 * correction[i]
-    rotatedbox[i][1] = rotatedbox[i][1] + 0.5 * correction[i]
+    rotatedbox[i][1] += 0.5 * correction[i]
     extent[i] = rotatedbox[i][1] - rotatedbox[i][0]
   NpointsPerSlice = Npoints[0] * Npoints[1]
   totalNpoints = NpointsPerSlice * Npoints[2]
@@ -301,20 +279,20 @@ for filename in filenames:
   if options.verbose: sys.stdout.write("\nGENERATING POINTS FOR POINT GRID")
   points = vtk.vtkPoints()
   for k in range(Npoints[2]):
-    for j in range(Npoints[0]):  
+    for j in range(Npoints[0]):
       for i in range(Npoints[1]):    # y is fastest index
         rotatedpoint = np.array([rotatedbox[0][0] + (float(j) + 0.5) * options.resolution,
-                                    rotatedbox[1][0] + (float(i) + 0.5) * options.resolution,
+                                 rotatedbox[1][0] + (float(i) + 0.5) * options.resolution,
-                                    rotatedbox[2][0] + (float(k) + 0.5) * options.distance ])       # point in rotated system
+                                 rotatedbox[2][0] + (float(k) + 0.5) * options.distance ])       # point in rotated system
         point = np.dot(R.T,rotatedpoint)                                                            # point in mesh system
         points.InsertNextPoint(list(point))
-        if options.verbose: 
+        if options.verbose:
           sys.stdout.write("\rGENERATING POINTS FOR POINT GRID %d%%" %(100*(Npoints[1]*(k*Npoints[0]+j)+i+1)/totalNpoints))
           sys.stdout.flush()
-  if options.verbose: 
+  if options.verbose:
     sys.stdout.write("\n  number of slices: %i\n"%Npoints[2])
     sys.stdout.write("  slice spacing: %.8f\n"%options.distance)
-    if Npoints[2] > 1: 
+    if Npoints[2] > 1:
       sys.stdout.write("  number of points per slice: %i = %i rows * %i points in row\n"%(NpointsPerSlice,Npoints[0],Npoints[1]))
     sys.stdout.write("  grid resolution: %.8f\n"%options.resolution)
 
@@ -324,7 +302,7 @@ for filename in filenames:
     vertex = vtk.vtkVertex()
     vertex.GetPointIds().SetId(0,i)  # each vertex consists of exactly one (index 0) point with ID "i"
     vertices.InsertNextCell(vertex)
-    if options.verbose: 
+    if options.verbose:
       sys.stdout.write("\rGENERATING VERTICES FOR POINT GRID %d%%" %(100*(i+1)/totalNpoints))
       sys.stdout.flush()
 
@@ -357,34 +335,34 @@ for filename in filenames:
     if enclosedPoints.IsInside(i):
       NenclosedPoints += 1
 # here one could use faster(?) "FindClosestPoint" if only first nearest neighbor required
-      kdTree.FindClosestNPoints(options.interpolation,pointgrid.GetPoint(i),ids) 
+      kdTree.FindClosestNPoints(1,pointgrid.GetPoint(i),ids)
-      for j in range(ids.GetNumberOfIds()): 
+      for j in range(ids.GetNumberOfIds()):
         gridToMesh[-1].extend([ids.GetId(j)])
-    if options.verbose: 
+    if options.verbose:
       sys.stdout.write("\rBUILDING MAPPING OF GRID POINTS %d%%" %(100*(i+1)/totalNpoints))
       sys.stdout.flush()
   if options.verbose:
     sys.stdout.write("\n  Number of points inside mesh geometry %i\n"%NenclosedPoints)
     sys.stdout.write("  Number of points outside mesh geometry %i\n"%(totalNpoints - NenclosedPoints))
-  
+
 
 
 # ITERATE OVER SLICES AND CREATE ANG FILE
-  
+
-  if options.verbose: 
+  if options.verbose:
     sys.stdout.write("\nWRITING OUT ANG FILES\n")
     sys.stdout.write("  scaling all length with %f\n"%options.scale)
   x0,y0,z0 = np.dot(R,pointgrid.GetPoint(0))                    # first point on slice defines origin
   for sliceN in range(Npoints[2]):
-    
+
     # Open file and write header
-    
+
     angfilename = eval('"'+eval("'%%s_slice%%0%ii.ang'%(math.log10(Npoints[2])+1)")+'"%(os.path.splitext(filename)[0],sliceN+1)')
     with open(angfilename,'w') as angfile:
       if options.verbose: sys.stdout.write("  %s\n"%angfilename)
       angfile.write(getHeader(filename,Npoints[1],Npoints[0],options.resolution*options.scale))
       for i in range(sliceN*NpointsPerSlice,(sliceN+1)*NpointsPerSlice):   # Iterate over points on slice
-        
+
 
         # Get euler angles of closest IDs
 
@@ -397,26 +375,26 @@ for filename in filenames:
               phi[-1].extend([angles[options.eulerLabel[k]].GetValue(IP)])
         else:
           phi = [[720,720,720]] # fake angles
-        
+
-        
+
         # Interpolate Euler angle
         # NOT YET IMPLEMENTED, simply take the nearest neighbors values
 
         interpolatedPhi = phi[0]
-        
+
-        
+
         # write data to ang file
 
         x,y,z = np.dot(R,pointgrid.GetPoint(i))                  # point in rotated TSL system
-        x -= x0                                                     # first point on slice defines origin 
+        x -= x0                                                     # first point on slice defines origin
-        y -= y0                                                     # first point on slice defines origin 
+        y -= y0                                                     # first point on slice defines origin
         x *= options.scale
         y *= options.scale
         angfile.write(getDataLine(interpolatedPhi,x,y,enclosedPoints.IsInside(i)))
-      
+
 
 # Visualize slices
-  
+
   if options.visualize:
     meshMapper = vtk.vtkDataSetMapper()
     meshMapper.SetInput(surface)
@@ -426,11 +404,11 @@ for filename in filenames:
     meshActor.GetProperty().SetOpacity(0.2)
     meshActor.GetProperty().SetColor(1.0,1.0,0)
     meshActor.GetProperty().BackfaceCullingOn()
-    
+
     boxpoints = vtk.vtkPoints()
     for n in range(8):
       P = [rotatedbox[0][(n/1)%2],
-           rotatedbox[1][(n/2)%2], 
+           rotatedbox[1][(n/2)%2],
            rotatedbox[2][(n/4)%2]]
       boxpoints.InsertNextPoint(list(np.dot(R.T,np.array(P))))
     box = vtk.vtkHexahedron()
@@ -469,7 +447,7 @@ for filename in filenames:
     renderer.AddActor(boxActor)
     renderer.AddActor(gridActor)
     renderer.SetBackground(1,1,1)
-     
+
     renderWindow.Render()
     renderWindowInteractor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
     renderWindowInteractor.Start()

processing/pre/geom_addPrimitive.py

@@ -1,135 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import numpy as np
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [geomfile(s)]', description = """
-Inserts a primitive geometric object at a given position.
-These objects can be boxes, cylinders, or ellipsoids.
-
-""", version = scriptID)
-
-parser.add_option('-c', '--center',
-                  dest='center',
-                  type='float', nargs = 3, metavar=' '.join(['float']*3),
-                  help='a,b,c origin of primitive %default')
-parser.add_option('-d', '--dimension',
-                  dest='dimension',
-                  type='float', nargs = 3, metavar=' '.join(['float']*3),
-                  help='a,b,c extension of hexahedral box')
-parser.add_option('-e', '--exponent',
-                  dest='exponent',
-                  type='float', nargs = 3, metavar=' '.join(['float']*3),
-                  help='i,j,k exponents for axes: '+
-                  '0 gives octahedron (|x|^(2^0) + |y|^(2^0) + |z|^(2^0) < 1), '+
-                  '1 gives a sphere (|x|^(2^1) + |y|^(2^1) + |z|^(2^1) < 1), '+
-                  'large values produce boxes, negative turn concave.')
-parser.add_option('-f', '--fill',
-                  dest='fill',
-                  type='float', metavar = 'int',
-                  help='microstructure index to fill primitive, defaults to max microstructure index + 1')
-parser.add_option('-q', '--quaternion',
-                  dest='quaternion',
-                  type='float', nargs = 4, metavar=' '.join(['float']*4),
-                  help = 'rotation of primitive as quaternion')
-parser.add_option('-a', '--angleaxis',
-                  dest='angleaxis',
-                  type=float, nargs = 4, metavar=' '.join(['float']*4),
-                  help = 'axis and angle to rotate primitive')
-parser.add_option(     '--degrees',
-                  dest='degrees',
-                  action='store_true',
-                  help = 'angle is given in degrees')
-parser.add_option(     '--nonperiodic',
-                  dest='periodic',
-                  action='store_false',
-                  help = 'wrap around edges')
-parser.add_option(     '--realspace',
-                  dest='realspace',
-                  action='store_true',
-                  help = '-c and -d span [origin,origin+size] instead of [0,grid] coordinates')
-parser.add_option(     '--invert',
-                  dest='inside',
-                  action='store_false',
-                  help = 'invert the volume filled by the primitive (inside/outside)')
-
-parser.set_defaults(center = (.0,.0,.0),
-                    degrees = False,
-                    exponent = (20,20,20), # box shape by default
-                    periodic = True,
-                    realspace = False,
-                    inside = True,
-                   )
-
-(options, filenames) = parser.parse_args()
-
-if options.dimension is None:
-  parser.error('no dimension specified.')
-if [options.angleaxis,options.quaternion].count(None) == 0:
-  parser.error('more than one rotation specified.')
-
-if options.angleaxis is not None:
-  rotation = damask.Rotation.from_axis_angle(np.array(options.angleaxis),options.degrees,normalise=True)
-elif options.quaternion is not None:
-  rotation = damask.Rotation.from_quaternion(options.quaternion)
-else:
-  rotation = damask.Rotation()
-
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-  damask.util.report(scriptName,name)
-
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-  grid = geom.get_grid()
-  size = geom.get_size()
-
-  # scale to box of size [1.0,1.0,1.0]
-  if options.realspace:
-    center = (np.array(options.center) - geom.get_origin())/size
-    r = np.array(options.dimension)/size/2.0
-  else:
-    center = (np.array(options.center) + 0.5)/grid
-    r = np.array(options.dimension)/grid/2.0
-
-  if np.any(center<0.0) or np.any(center>=1.0): print('error')
-
-  offset = np.ones(3)*0.5 if options.periodic else center
-  mask = np.full(grid,False)
-  # High exponents can cause underflow & overflow - okay here, just compare to 1, so +infinity and 0 are fine
-  np.seterr(over='ignore', under='ignore')
-
-  e = 2.0**np.array(options.exponent)
-  for x in range(grid[0]):
-    for y in range(grid[1]):
-      for z in range(grid[2]):
-        coords = np.array([x+0.5,y+0.5,z+0.5])/grid
-        mask[x,y,z] = np.sum(np.abs((rotation*(coords-offset))/r)**e) < 1
-
-  if options.periodic:
-    shift = ((offset-center)*grid).astype(int)
-    mask = np.roll(mask,shift,(0,1,2))
-
-  if options.inside: mask = np.logical_not(mask)
-  fill = np.nanmax(geom.microstructure)+1 if options.fill is None else options.fill
-
-  damask.util.croak(geom.update(np.where(mask,geom.microstructure,fill)))
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-
-  geom.to_file(sys.stdout if name is None else name,pack=False)

processing/pre/geom_fromDREAM3D.py

@@ -154,4 +154,4 @@ for name in filenames:
                      homogenization=options.homogenization,comments=header)
   damask.util.croak(geom)
 
-  geom.to_file(os.path.splitext(name)[0]+'.geom',pack=False)
+  geom.to_file(os.path.splitext(name)[0]+'.geom',format='ASCII',pack=False)

processing/pre/geom_fromMinimalSurface.py

@@ -89,4 +89,4 @@ geom=damask.Geom(microstructure,options.size,
                  comments=[scriptID + ' ' + ' '.join(sys.argv[1:])])
 damask.util.croak(geom)
 
-geom.to_file(sys.stdout if name is None else name,pack=False)
+geom.to_file(sys.stdout if name is None else name,format='ASCII',pack=False)

processing/pre/geom_fromOsteonGeometry.py

@@ -142,4 +142,4 @@ geom = damask.Geom(microstructure.reshape(grid),
                    homogenization=options.homogenization,comments=header)
 damask.util.croak(geom)
 
-geom.to_file(sys.stdout if name is None else name,pack=False)
+geom.to_file(sys.stdout if name is None else name,format='ASCII',pack=False)

processing/pre/geom_fromTable.py

@@ -105,4 +105,5 @@ for name in filenames:
                        homogenization=options.homogenization,comments=header)
     damask.util.croak(geom)
 
-    geom.to_file(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',pack=False)
+    geom.to_file(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',
+                 format='ASCII',pack=False)

processing/pre/geom_fromVoronoiTessellation.py

@@ -228,4 +228,5 @@ for name in filenames:
                        homogenization=options.homogenization,comments=header)
     damask.util.croak(geom)
 
-    geom.to_file(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',pack=False)
+    geom.to_file(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',
+                 format='ASCII',pack=False)

processing/pre/geom_grainGrowth.py

@@ -169,7 +169,7 @@ for name in filenames:
     # undo any changes involving immutable microstructures
     microstructure = np.where(immutable, microstructure_original,microstructure)
 
-  damask.util.croak(geom.update(microstructure[0:grid_original[0],0:grid_original[1],0:grid_original[2]]))
+  geom=geom.duplicate(microstructure[0:grid_original[0],0:grid_original[1],0:grid_original[2]])
   geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
 
-  geom.to_file(sys.stdout if name is None else name,pack=False)
+  geom.to_file(sys.stdout if name is None else name,format='ASCII',pack=False)

processing/pre/geom_vicinityOffset.py

@@ -1,85 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-from scipy import ndimage
-import numpy as np
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-def taintedNeighborhood(stencil,trigger=[],size=1):
-
-  me = stencil[stencil.shape[0]//2]
-  if len(trigger) == 0:
-    return np.any(stencil != me)
-  if me in trigger:
-    trigger = set(trigger)
-    trigger.remove(me)
-    trigger = list(trigger)
-  return np.any(np.in1d(stencil,np.array(trigger)))
-
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [geomfile(s)]', description = """
-Offset microstructure index for points which see a microstructure different from themselves
-(or listed as triggers) within a given (cubic) vicinity, i.e. within the region close to a grain/phase boundary.
-
-""", version = scriptID)
-
-parser.add_option('-v', '--vicinity',
-                  dest = 'vicinity',
-                  type = 'int', metavar = 'int',
-                  help = 'voxel distance checked for presence of other microstructure [%default]')
-parser.add_option('-o', '--offset',
-                  dest='offset',
-                  type = 'int', metavar = 'int',
-                  help='offset (positive or negative) to tag microstructure indices, defaults to max microstructure index')
-parser.add_option('-t', '--trigger',
-                  dest = 'trigger',
-                  action = 'extend', metavar = '<int LIST>',
-                  help = 'list of microstructure indices triggering a change')
-parser.add_option('-n', '--nonperiodic',
-                  dest = 'mode',
-                  action = 'store_const', const = 'nearest',
-                  help = 'assume geometry to be non-periodic')
-
-parser.set_defaults(vicinity = 1,
-                    trigger  = [],
-                    mode     = 'wrap',
-                   )
-
-(options, filenames) = parser.parse_args()
-
-options.trigger = np.array(options.trigger, dtype=int)
-
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-  damask.util.report(scriptName,name)
-
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-
-  offset = np.nanmax(geom.microstructure) if options.offset is None else options.offset
-
-  damask.util.croak(geom.update(np.where(ndimage.filters.generic_filter(
-                                           geom.microstructure,
-                                           taintedNeighborhood,
-                                           size=1+2*options.vicinity,mode=options.mode,
-                                           extra_arguments=(),
-                                           extra_keywords={"trigger":options.trigger,"size":1+2*options.vicinity}),
-                                           geom.microstructure + offset,geom.microstructure)))
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-
-  geom.to_file(sys.stdout if name is None else name,pack=False)

processing/pre/mentat_pbcOnBoxMesh.py

@@ -1,9 +1,14 @@
 #!/usr/bin/env python3
-# -*- coding: UTF-8 no BOM -*-
 
-import sys,os,re,time,tempfile
+import sys
-import numpy as np
+import os
+import re
+import time
+import tempfile
 from optparse import OptionParser
+
+import numpy as np
+
 import damask
 
 sys.path.append(str(damask.solver.Marc().library_path))

processing/pre/mentat_spectralBox.py

@@ -20,7 +20,6 @@ def outMentat(cmd,locals):
     py_mentat.py_send(cmd)
   else:
     py_mentat.py_send(cmd)
-  return
 
 #-------------------------------------------------------------------------------------------------
 def outFile(cmd,locals,dest):
@@ -31,7 +30,6 @@ def outFile(cmd,locals,dest):
     dest.write(cmd+'\n')
   else:
     dest.write(cmd+'\n')
-  return
 
 #-------------------------------------------------------------------------------------------------
 def output(cmds,locals,dest):
@@ -43,8 +41,6 @@ def output(cmds,locals,dest):
         outMentat(str(cmd),locals)
       else:
         outFile(str(cmd),locals,dest)
-  return
-
 
   
 #-------------------------------------------------------------------------------------------------
@@ -105,7 +101,7 @@ def mesh(r,d):
 
 #-------------------------------------------------------------------------------------------------
 def material():
-  cmds = [\
+  return [\
   "*new_mater standard",
   "*mater_option general:state:solid",
   "*mater_option structural:type:hypo_elast",
@@ -119,12 +115,10 @@ def material():
   "all_existing",
   ]
   
-  return cmds
-  
 
 #-------------------------------------------------------------------------------------------------
 def geometry():
-  cmds = [\
+  return [\
   "*geometry_type mech_three_solid",
 #  "*geometry_option red_integ_capacity:on",
   "*add_geometry_elements",
@@ -134,8 +128,6 @@ def geometry():
   "all_existing",
   ]
   
-  return cmds
-  
 
 #-------------------------------------------------------------------------------------------------
 def initial_conditions(homogenization,microstructures):

processing/pre/seeds_fromDistribution.py

@@ -90,7 +90,7 @@ class myThread (threading.Thread):
           direction[i]*=2.
           i+= 1
       perturbedSeedsTable.set('pos',coords)
-      perturbedSeedsTable.to_ASCII(perturbedSeedsVFile)
+      perturbedSeedsTable.to_file(perturbedSeedsVFile)
 
 #--- do tesselation with perturbed seed file ------------------------------------------------------
       perturbedGeomVFile.close()

processing/pre/seeds_fromGeom.py

@@ -65,5 +65,5 @@ for name in filenames:
 
     table = damask.Table(seeds[mask],{'pos':(3,)},comments)
     table.add('microstructure',microstructure[mask])
-    table.to_ASCII(sys.stdout if name is None else \
+    table.to_file(sys.stdout if name is None else \
-                   os.path.splitext(name)[0]+'.seeds')
+                  os.path.splitext(name)[0]+'.seeds')

processing/pre/seeds_fromPokes.py

@@ -92,5 +92,5 @@ for name in filenames:
 
     table = damask.Table(seeds,{'pos':(3,),'microstructure':(1,)},comments)
     table.set('microstructure',table.get('microstructure').astype(np.int))
-    table.to_ASCII(sys.stdout if name is None else \
+    table.to_file(sys.stdout if name is None else \
-                   os.path.splitext(name)[0]+'_poked_{}.seeds'.format(options.N))
+                  os.path.splitext(name)[0]+f'_poked_{options.N}.seeds')

processing/pre/seeds_fromRandom.py

@@ -162,4 +162,4 @@ for name in filenames:
              else np.random.normal(loc = options.mean, scale = options.sigma, size = options.N)
         table.add('weight',weights)
 
-    table.to_ASCII(sys.stdout if name is None else name)
+    table.to_file(sys.stdout if name is None else name)

python/.coveragerc

@@ -2,4 +2,3 @@
 omit = tests/*
        damask/_asciitable.py
        damask/_test.py
-       damask/config/*

python/damask/__init__.py

@@ -24,5 +24,4 @@ from .             import solver           # noqa
 Environment = _
 from ._asciitable  import ASCIItable       # noqa
 from ._test        import Test             # noqa
-from .config       import Material         # noqa
 from .util         import extendableOption # noqa

python/damask/_colormap.py

@@ -8,8 +8,9 @@ if os.name == 'posix' and 'DISPLAY' not in os.environ:
     mpl.use('Agg')
 import matplotlib.pyplot as plt
 from matplotlib import cm
+from PIL import Image
 
-import damask
+from . import util
 from . import Table
 
 _eps   = 216./24389.
@@ -50,7 +51,7 @@ class Colormap(mpl.colors.ListedColormap):
             Color definition for minimum value.
         high : numpy.ndarray of shape (3)
             Color definition for maximum value.
-        N : integer, optional
+        N : int, optional
             The number of color quantization levels. Defaults to 256.
         name : str, optional
             The name of the colormap. Defaults to `DAMASK colormap`.
@@ -161,6 +162,47 @@ class Colormap(mpl.colors.ListedColormap):
             print('  '+', '.join(cat[1]))
 
 
+    def shade(self,field,bounds=None,gap=None):
+        """
+        Generate PIL image of 2D field using colormap.
+
+        Parameters
+        ----------
+        field : numpy.array of shape(:,:)
+            Data to be shaded.
+        bounds : iterable of len(2), optional
+            Colormap value range (low,high).
+        gap : field.dtype, optional
+            Transparent value. NaN will always be rendered transparent.
+
+        Returns
+        -------
+        PIL.Image
+            RGBA image of shaded data.
+
+        """
+        N = len(self.colors)
+        mask = np.logical_not(np.isnan(field) if gap is None else \
+               np.logical_or (np.isnan(field), field == gap))                                       # mask NaN (and gap if present)
+
+        lo,hi = (field[mask].min(),field[mask].max()) if bounds is None else \
+                (min(bounds[:2]),max(bounds[:2]))
+
+        delta,avg = hi-lo,0.5*(hi+lo)
+
+        if delta * 1e8 <= avg:                                                                      # delta is similar to numerical noise
+            hi,lo = hi+0.5*avg,lo-0.5*avg                                                           # extend range to have actual data centered within
+
+        return Image.fromarray(
+            (np.dstack((
+                        self.colors[(np.round(np.clip((field-lo)/(hi-lo),0.0,1.0)*(N-1))).astype(np.uint16),:3],
+                        mask.astype(float)
+                       )
+                      )*255
+            ).astype(np.uint8),
+            mode='RGBA')
+
+
     def show(self,aspect=10,vertical=False):
         """Show colormap as matplotlib figure."""
         fig = plt.figure(figsize=(5/aspect,5) if vertical else (5,5/aspect))
@@ -238,7 +280,7 @@ class Colormap(mpl.colors.ListedColormap):
             colors+=[i]+c
 
         out = [{
-                'Creator':f'damask.Colormap v{damask.version}',
+                'Creator':util.execution_stamp('Colormap'),
                 'ColorSpace':'RGB',
                 'Name':colormap.name,
                 'DefaultMap':True,
@@ -254,13 +296,13 @@ class Colormap(mpl.colors.ListedColormap):
     def _export_ASCII(colormap,fhandle=None):
         """Write colormap to ASCII table."""
         labels = {'RGBA':4} if colormap.colors.shape[1] == 4 else {'RGB': 3}
-        t = Table(colormap.colors,labels,f'Creator: damask.Colormap v{damask.version}')
+        t = Table(colormap.colors,labels,f'Creator: {util.execution_stamp("Colormap")}')
 
         if fhandle is None:
             with open(colormap.name.replace(' ','_')+'.txt', 'w') as f:
-                t.to_ASCII(f,True)
+                t.to_file(f,new_style=True)
         else:
-            t.to_ASCII(fhandle,True)
+            t.to_file(fhandle,new_style=True)
 
     @staticmethod
     def _export_GOM(colormap,fhandle=None):

python/damask/_geom.py

@@ -1,12 +1,14 @@
 import sys
-from io import StringIO
+import copy
 import multiprocessing
+from io import StringIO
 from functools import partial
 
 import numpy as np
 from scipy import ndimage,spatial
 
-import damask
+from . import environment
+from . import Rotation
 from . import VTK
 from . import util
 from . import grid_filters
@@ -22,15 +24,15 @@ class Geom:
         Parameters
         ----------
         microstructure : numpy.ndarray
-            microstructure array (3D)
+            Microstructure array (3D)
         size : list or numpy.ndarray
-            physical size of the microstructure in meter.
+            Physical size of the microstructure in meter.
         origin : list or numpy.ndarray, optional
-            physical origin of the microstructure in meter.
+            Physical origin of the microstructure in meter.
-        homogenization : integer, optional
+        homogenization : int, optional
-            homogenization index.
+            Homogenization index.
         comments : list of str, optional
-            comments lines.
+            Comment lines.
 
         """
         self.set_microstructure(microstructure)
@@ -43,74 +45,91 @@ class Geom:
     def __repr__(self):
         """Basic information on geometry definition."""
         return util.srepr([
-               f'grid     a b c:      {util.srepr(self.get_grid  ()," x ")}',
+               f'grid     a b c:     {util.srepr(self.get_grid  ()," x ")}',
-               f'size     x y z:      {util.srepr(self.get_size  ()," x ")}',
+               f'size     x y z:     {util.srepr(self.get_size  ()," x ")}',
-               f'origin   x y z:      {util.srepr(self.get_origin(),"   ")}',
+               f'origin   x y z:     {util.srepr(self.get_origin(),"   ")}',
-               f'homogenization:      {self.get_homogenization()}',
+               f'# materialpoints:   {self.N_microstructure}',
-               f'# microstructures:   {self.N_microstructure}',
+               f'max materialpoint:  {np.nanmax(self.microstructure)}',
-               f'max microstructure:  {np.nanmax(self.microstructure)}',
               ])
 
 
-    def update(self,microstructure=None,size=None,origin=None,rescale=False):
+    def __copy__(self):
+        """Copy geometry."""
+        return copy.deepcopy(self)
+
+
+    def copy(self):
+        """Copy geometry."""
+        return self.__copy__()
+
+
+    def duplicate(self,microstructure=None,size=None,origin=None,comments=None,autosize=False):
         """
-        Update microstructure and size.
+        Create a duplicate having updated microstructure, size, and origin.
 
         Parameters
         ----------
         microstructure : numpy.ndarray, optional
-            microstructure array (3D).
+            Microstructure array (3D).
         size : list or numpy.ndarray, optional
-            physical size of the microstructure in meter.
+            Physical size of the microstructure in meter.
         origin : list or numpy.ndarray, optional
-            physical origin of the microstructure in meter.
+            Physical origin of the microstructure in meter.
-        rescale : bool, optional
+        comments : list of str, optional
-            ignore size parameter and rescale according to change of grid points.
+            Comment lines.
+        autosize : bool, optional
+            Ignore size parameter and rescale according to change of grid points.
 
         """
-        grid_old    = self.get_grid()
+        if size is not None and autosize:
-        size_old    = self.get_size()
+            raise ValueError('Auto-sizing conflicts with explicit size parameter.')
-        origin_old  = self.get_origin()
-        unique_old  = self.N_microstructure
-        max_old     = np.nanmax(self.microstructure)
 
-        if size is not None and rescale:
+        grid_old = self.get_grid()
-            raise ValueError('Either set size explicitly or rescale automatically')
+        dup = self.copy()
+        dup.set_microstructure(microstructure)
+        dup.set_origin(origin)
 
-        self.set_microstructure(microstructure)
+        if comments is not None:
-        self.set_origin(origin)
+            dup.set_comments(comments)
 
         if size is not None:
-            self.set_size(size)
+            dup.set_size(size)
-        elif rescale:
+        elif autosize:
-            self.set_size(self.get_grid()/grid_old*self.size)
+            dup.set_size(dup.get_grid()/grid_old*self.get_size())
 
-        message = [f'grid     a b c:      {util.srepr(grid_old," x ")}']
+        return dup
-        if np.any(grid_old != self.get_grid()):
-            message[-1] = util.delete(message[-1])
-            message.append(util.emph(f'grid     a b c:      {util.srepr(self.get_grid()," x ")}'))
 
-        message.append(f'size     x y z:      {util.srepr(size_old," x ")}')
-        if np.any(size_old != self.get_size()):
-            message[-1] = util.delete(message[-1])
-            message.append(util.emph(f'size     x y z:      {util.srepr(self.get_size()," x ")}'))
 
-        message.append(f'origin   x y z:      {util.srepr(origin_old,"   ")}')
+    def diff(self,other):
-        if np.any(origin_old != self.get_origin()):
+        """
-            message[-1] = util.delete(message[-1])
+        Report property differences of self relative to other.
-            message.append(util.emph(f'origin   x y z:      {util.srepr(self.get_origin(),"   ")}'))
 
-        message.append(f'homogenization:      {self.get_homogenization()}')
+        Parameters
+        ----------
+        other : Geom
+            Geometry to compare self against.
 
-        message.append(f'# microstructures:   {unique_old}')
+        """
-        if unique_old != self.N_microstructure:
+        message = []
-            message[-1] = util.delete(message[-1])
+        if np.any(other.get_grid() != self.get_grid()):
-            message.append(util.emph(f'# microstructures:   {self.N_microstructure}'))
+            message.append(util.delete(f'grid     a b c:     {util.srepr(other.get_grid()," x ")}'))
+            message.append(util.emph(  f'grid     a b c:     {util.srepr( self.get_grid()," x ")}'))
 
-        message.append(f'max microstructure:  {max_old}')
+        if np.any(other.get_size() != self.get_size()):
-        if max_old != np.nanmax(self.microstructure):
+            message.append(util.delete(f'size     x y z:     {util.srepr(other.get_size()," x ")}'))
-            message[-1] = util.delete(message[-1])
+            message.append(util.emph(  f'size     x y z:     {util.srepr( self.get_size()," x ")}'))
-            message.append(util.emph(f'max microstructure:  {np.nanmax(self.microstructure)}'))
+
+        if np.any(other.get_origin() != self.get_origin()):
+            message.append(util.delete(f'origin   x y z:     {util.srepr(other.get_origin(),"   ")}'))
+            message.append(util.emph(  f'origin   x y z:     {util.srepr( self.get_origin(),"   ")}'))
+
+        if other.N_microstructure != self.N_microstructure:
+            message.append(util.delete(f'# materialpoints:   {other.N_microstructure}'))
+            message.append(util.emph(  f'# materialpoints:   { self.N_microstructure}'))
+
+        if np.nanmax(other.microstructure) != np.nanmax(self.microstructure):
+            message.append(util.delete(f'max materialpoint:  {np.nanmax(other.microstructure)}'))
+            message.append(util.emph(  f'max materialpoint:  {np.nanmax( self.microstructure)}'))
 
         return util.return_message(message)
 
@@ -122,7 +141,7 @@ class Geom:
         Parameters
         ----------
         comments : list of str
-            new comments.
+            All comments.
 
         """
         self.comments = []
@@ -136,7 +155,7 @@ class Geom:
         Parameters
         ----------
         comments : list of str
-            new comments.
+            New comments.
 
         """
         self.comments += [str(c) for c in comments] if isinstance(comments,list) else [str(comments)]
@@ -146,20 +165,32 @@ class Geom:
         """
         Replace the existing microstructure representation.
 
+        The complete microstructure is replaced (indcluding grid definition),
+        unless a masked array is provided in which case the grid dimensions
+        need to match and masked entries are not replaced.
+
         Parameters
         ----------
-        microstructure : numpy.ndarray
+        microstructure : numpy.ndarray or numpy.ma.core.MaskedArray of shape (:,:,:)
-            microstructure array (3D).
+            Microstructure indices.
 
         """
         if microstructure is not None:
-            if len(microstructure.shape) != 3:
+            if isinstance(microstructure,np.ma.core.MaskedArray):
-                raise ValueError(f'Invalid microstructure shape {microstructure.shape}')
+                self.microstructure = np.where(microstructure.mask,
-            elif microstructure.dtype not in np.sctypes['float'] + np.sctypes['int']:
+                                               self.microstructure,microstructure.data)
-                raise TypeError(f'Invalid microstructue data type {microstructure.dtype}')
             else:
                 self.microstructure = np.copy(microstructure)
 
+        if         self.microstructure.dtype in np.sctypes['float'] and \
+            np.all(self.microstructure == self.microstructure.astype(int).astype(float)):
+            self.microstructure = self.microstructure.astype(int)
+
+        if len(self.microstructure.shape) != 3:
+            raise ValueError(f'Invalid microstructure shape {microstructure.shape}')
+        elif self.microstructure.dtype not in np.sctypes['float'] + np.sctypes['int']:
+            raise TypeError(f'Invalid microstructure data type {microstructure.dtype}')
+
 
     def set_size(self,size):
         """
@@ -168,13 +199,10 @@ class Geom:
         Parameters
         ----------
         size : list or numpy.ndarray
-            physical size of the microstructure in meter.
+            Physical size of the microstructure in meter.
 
         """
-        if size is None:
+        if size is not None:
-            grid = np.asarray(self.microstructure.shape)
-            self.size = grid/np.max(grid)
-        else:
             if len(size) != 3 or any(np.array(size) <= 0):
                 raise ValueError(f'Invalid size {size}')
             else:
@@ -188,7 +216,7 @@ class Geom:
         Parameters
         ----------
         origin : list or numpy.ndarray
-            physical origin of the microstructure in meter
+            Physical origin of the microstructure in meter.
 
         """
         if origin is not None:
@@ -204,13 +232,13 @@ class Geom:
 
         Parameters
         ----------
-        homogenization : integer
+        homogenization : int
-            homogenization index
+            Homogenization index.
 
         """
         if homogenization is not None:
             if not isinstance(homogenization,int) or homogenization < 1:
-                raise TypeError(f'Invalid homogenization {homogenization}')
+                raise TypeError(f'Invalid homogenization {homogenization}.')
             else:
                 self.homogenization = homogenization
 
@@ -255,16 +283,6 @@ class Geom:
         return self.comments[:]
 
 
-    def get_header(self):
-        """Return the full header (grid, size, origin, homogenization, comments)."""
-        header =  [f'{len(self.comments)+4} header'] + self.comments
-        header.append('grid   a {} b {} c {}'.format(*self.get_grid()))
-        header.append('size   x {} y {} z {}'.format(*self.get_size()))
-        header.append('origin x {} y {} z {}'.format(*self.get_origin()))
-        header.append(f'homogenization {self.get_homogenization()}')
-        return header
-
-
     @staticmethod
     def from_file(fname):
         """
@@ -273,7 +291,7 @@ class Geom:
         Parameters
         ----------
         fname : str or file handle
-            geometry file to read.
+            Geometry file to read.
 
         """
         try:
@@ -282,13 +300,16 @@ class Geom:
             f = fname
 
         f.seek(0)
-        header_length,keyword = f.readline().split()[:2]
+        try:
-        header_length = int(header_length)
+            header_length,keyword = f.readline().split()[:2]
-        content = f.readlines()
+            header_length = int(header_length)
-
+        except ValueError:
+            header_length,keyword = (-1, 'invalid')
         if not keyword.startswith('head') or header_length < 3:
             raise TypeError('Header length information missing or invalid')
 
+        content = f.readlines()
+
         comments = []
         for i,line in enumerate(content[:header_length]):
             items = line.split('#')[0].lower().strip().split()
@@ -328,6 +349,27 @@ class Geom:
         return Geom(microstructure.reshape(grid,order='F'),size,origin,homogenization,comments)
 
 
+    @staticmethod
+    def from_vtr(fname):
+        """
+        Read a VTK rectilinear grid.
+
+        Parameters
+        ----------
+        fname : str or or pathlib.Path
+            Geometry file to read.
+            Valid extension is .vtr, it will be appended if not given.
+
+        """
+        v = VTK.from_file(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
+        comments = v.get_comments()
+        grid = np.array(v.vtk_data.GetDimensions())-1
+        bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
+        size = bbox[1] - bbox[0]
+
+        return Geom(v.get('materialpoint').reshape(grid,order='F'),size,bbox[0],comments=comments)
+
+
     @staticmethod
     def _find_closest_seed(seeds, weights, point):
         return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
@@ -339,16 +381,16 @@ class Geom:
 
         Parameters
         ----------
-        grid : numpy.ndarray of shape (3)
+        grid : int numpy.ndarray of shape (3)
-            number of grid points in x,y,z direction.
+            Number of grid points in x,y,z direction.
         size : list or numpy.ndarray of shape (3)
-            physical size of the microstructure in meter.
+            Physical size of the microstructure in meter.
         seeds : numpy.ndarray of shape (:,3)
-            position of the seed points in meter. All points need to lay within the box.
+            Position of the seed points in meter. All points need to lay within the box.
         weights : numpy.ndarray of shape (seeds.shape[0])
-            weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
+            Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
         periodic : Boolean, optional
-            perform a periodic tessellation. Defaults to True.
+            Perform a periodic tessellation. Defaults to True.
 
         """
         if periodic:
@@ -362,7 +404,7 @@ class Geom:
             seeds_p   = seeds
             coords    = grid_filters.cell_coord0(grid,size).reshape(-1,3)
 
-        pool = multiprocessing.Pool(processes = int(damask.environment.options['DAMASK_NUM_THREADS']))
+        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])
         pool.close()
         pool.join()
@@ -374,8 +416,9 @@ class Geom:
         else:
             microstructure = microstructure.reshape(grid)
 
-        #ToDo: comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version)
+        return Geom(microstructure+1,size,homogenization=1,
-        return Geom(microstructure+1,size,homogenization=1)
+                    comments=util.execution_stamp('Geom','from_Laguerre_tessellation'),
+                   )
 
 
     @staticmethod
@@ -385,115 +428,203 @@ class Geom:
 
         Parameters
         ----------
-        grid : numpy.ndarray of shape (3)
+        grid : int numpy.ndarray of shape (3)
-            number of grid points in x,y,z direction.
+            Number of grid points in x,y,z direction.
         size : list or numpy.ndarray of shape (3)
-            physical size of the microstructure in meter.
+            Physical size of the microstructure in meter.
         seeds : numpy.ndarray of shape (:,3)
-            position of the seed points in meter. All points need to lay within the box.
+            Position of the seed points in meter. All points need to lay within the box.
         periodic : Boolean, optional
-            perform a periodic tessellation. Defaults to True.
+            Perform a periodic tessellation. Defaults to True.
 
         """
         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)
 
-        #ToDo: comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version)
+        return Geom(microstructure.reshape(grid)+1,size,homogenization=1,
-        return Geom(microstructure.reshape(grid)+1,size,homogenization=1)
+                    comments=util.execution_stamp('Geom','from_Voronoi_tessellation'),
+                   )
 
 
-    def to_file(self,fname,pack=None):
+    def to_file(self,fname,format='vtr',pack=None):
         """
         Writes a geom file.
 
         Parameters
         ----------
         fname : str or file handle
-            geometry file to write.
+            Geometry file to write.
+        format : {'vtr', 'ASCII'}, optional
+            File format, defaults to 'vtr'. Available formats are:
+            - vtr: VTK rectilinear grid file, extension '.vtr'.
+            - ASCII: Plain text file, extension '.geom'.
         pack : bool, optional
-            compress geometry with 'x of y' and 'a to b'.
+            Compress ASCII geometry with 'x of y' and 'a to b'.
 
         """
-        header = self.get_header()
+        def _to_ASCII(geom,fname,pack=None):
-        grid   = self.get_grid()
+            """
+            Writes a geom file.
 
-        if pack is None:
+            Parameters
-            plain = grid.prod()/self.N_microstructure < 250
+            ----------
+            geom : Geom object
+                Geometry to write.
+            fname : str or file handle
+                Geometry file to write.
+            pack : bool, optional
+                Compress geometry with 'x of y' and 'a to b'.
+
+            """
+            header =  [f'{len(geom.comments)+4} header'] + geom.comments
+            header.append('grid   a {} b {} c {}'.format(*geom.get_grid()))
+            header.append('size   x {} y {} z {}'.format(*geom.get_size()))
+            header.append('origin x {} y {} z {}'.format(*geom.get_origin()))
+            header.append(f'homogenization {geom.get_homogenization()}')
+
+            grid = geom.get_grid()
+
+            if pack is None:
+                plain = grid.prod()/geom.N_microstructure < 250
+            else:
+                plain = not pack
+
+            if plain:
+                format_string = '%g' if geom.microstructure.dtype in np.sctypes['float'] else \
+                                '%{}i'.format(1+int(np.floor(np.log10(np.nanmax(geom.microstructure)))))
+                np.savetxt(fname,
+                           geom.microstructure.reshape([grid[0],np.prod(grid[1:])],order='F').T,
+                           header='\n'.join(header), fmt=format_string, comments='')
+            else:
+                try:
+                    f = open(fname,'w')
+                except TypeError:
+                    f = fname
+
+                compressType = None
+                former = start = -1
+                reps = 0
+                for current in geom.microstructure.flatten('F'):
+                    if abs(current - former) == 1 and (start - current) == reps*(former - current):
+                        compressType = 'to'
+                        reps += 1
+                    elif current == former and start == former:
+                        compressType = 'of'
+                        reps += 1
+                    else:
+                        if   compressType is None:
+                            f.write('\n'.join(header)+'\n')
+                        elif compressType == '.':
+                            f.write(f'{former}\n')
+                        elif compressType == 'to':
+                            f.write(f'{start} to {former}\n')
+                        elif compressType == 'of':
+                            f.write(f'{reps} of {former}\n')
+
+                        compressType = '.'
+                        start = current
+                        reps = 1
+
+                    former = current
+
+                if compressType == '.':
+                    f.write(f'{former}\n')
+                elif compressType == 'to':
+                    f.write(f'{start} to {former}\n')
+                elif compressType == 'of':
+                    f.write(f'{reps} of {former}\n')
+
+
+        def _to_vtr(geom,fname=None):
+            """
+            Generates vtk rectilinear grid.
+
+            Parameters
+            ----------
+            geom : Geom object
+                Geometry to write.
+            fname : str, optional
+                Filename to write. If no file is given, a string is returned.
+                Valid extension is .vtr, it will be appended if not given.
+
+            """
+            v = VTK.from_rectilinearGrid(geom.grid,geom.size,geom.origin)
+            v.add(geom.microstructure.flatten(order='F'),'materialpoint')
+            v.add_comments(geom.comments)
+
+            if fname:
+                v.to_file(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
+            else:
+                sys.stdout.write(v.__repr__())
+
+        if format.lower() == 'ascii':
+            return _to_ASCII(self,fname,pack)
+        elif format.lower() == 'vtr':
+            return _to_vtr(self,fname)
         else:
-            plain = not pack
+            raise TypeError(f'Unknown format {format}.')
 
-        if plain:
+    def as_ASCII(self,pack=False):
-            format_string = '%g' if self.microstructure.dtype in np.sctypes['float'] else \
+        """Format geometry as human-readable ASCII."""
-                            '%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.microstructure)))))
+        f = StringIO()
-            np.savetxt(fname,
+        self.to_file(f,'ASCII',pack)
-                       self.microstructure.reshape([grid[0],np.prod(grid[1:])],order='F').T,
+        f.seek(0)
-                       header='\n'.join(header), fmt=format_string, comments='')
+        return ''.join(f.readlines())
-        else:
-            try:
-                f = open(fname,'w')
-            except TypeError:
-                f = fname
-
-            compressType = None
-            former = start = -1
-            reps = 0
-            for current in self.microstructure.flatten('F'):
-                if abs(current - former) == 1 and (start - current) == reps*(former - current):
-                    compressType = 'to'
-                    reps += 1
-                elif current == former and start == former:
-                    compressType = 'of'
-                    reps += 1
-                else:
-                    if   compressType is None:
-                        f.write('\n'.join(self.get_header())+'\n')
-                    elif compressType == '.':
-                        f.write(f'{former}\n')
-                    elif compressType == 'to':
-                        f.write(f'{start} to {former}\n')
-                    elif compressType == 'of':
-                        f.write(f'{reps} of {former}\n')
-
-                    compressType = '.'
-                    start = current
-                    reps = 1
-
-                former = current
-
-            if compressType == '.':
-                f.write(f'{former}\n')
-            elif compressType == 'to':
-                f.write(f'{start} to {former}\n')
-            elif compressType == 'of':
-                f.write(f'{reps} of {former}\n')
 
 
-    def to_vtk(self,fname=None):
+    def add_primitive(self,dimension,center,exponent,
+                      fill=None,R=Rotation(),inverse=False,periodic=True):
         """
-        Generates vtk file.
+        Inserts a primitive geometric object at a given position.
 
         Parameters
         ----------
-        fname : str, optional
+        dimension : int or float numpy.ndarray of shape(3)
-            vtk file to write. If no file is given, a string is returned.
+            Dimension (diameter/side length) of the primitive. If given as
+            integers, grid point locations (cell centers) are addressed.
+            If given as floats, coordinates are addressed.
+        center : int or float numpy.ndarray of shape(3)
+            Center of the primitive. If given as integers, grid point
+            locations (cell centers) are addressed.
+            If given as floats, coordinates are addressed.
+        exponent : numpy.ndarray of shape(3) or float
+            Exponents for the three axis.
+            0 gives octahedron (|x|^(2^0) + |y|^(2^0) + |z|^(2^0) < 1)
+            1 gives a sphere (|x|^(2^1) + |y|^(2^1) + |z|^(2^1) < 1)
+        fill : int, optional
+            Fill value for primitive. Defaults to microstructure.max() + 1.
+        R : damask.Rotation, optional
+            Rotation of primitive. Defaults to no rotation.
+        inverse : Boolean, optional
+            Retain original microstructure within primitive and fill outside.
+            Defaults to False.
+        periodic : Boolean, optional
+            Repeat primitive over boundaries. Defaults to True.
 
         """
-        v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
+        # normalized 'radius' and center
-        v.add(self.microstructure.flatten(order='F'),'microstructure')
+        r = np.array(dimension)/self.grid/2.0 if np.array(dimension).dtype in np.sctypes['int'] else \
+            np.array(dimension)/self.size/2.0
+        c = (np.array(center) + .5)/self.grid if np.array(center).dtype in np.sctypes['int'] else \
+            (np.array(center) - self.origin)/self.size
 
-        if fname:
+        coords = grid_filters.cell_coord0(self.grid,np.ones(3)) \
-            v.write(fname)
+               - ((np.ones(3)-(1./self.grid if np.array(center).dtype in np.sctypes['int'] else 0))*0.5 if periodic else c)  # periodic center is always at CoG
-        else:
+        coords_rot = R.broadcast_to(tuple(self.grid))@coords
-            sys.stdout.write(v.__repr__())
 
+        with np.errstate(all='ignore'):
+            mask = np.where(np.sum(np.power(coords_rot/r,2.0**exponent),axis=-1) > 1.0,True,False)
 
-    def show(self):
+        if periodic:                                                                                # translate back to center
-        """Show raw content (as in file)."""
+            mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))
-        f=StringIO()
+
-        self.to_file(f)
+        fill_ = np.full_like(self.microstructure,np.nanmax(self.microstructure)+1 if fill is None else fill)
-        f.seek(0)
+        ms = np.ma.MaskedArray(fill_,np.logical_not(mask) if inverse else mask)
-        return ''.join(f.readlines())
+
+        return self.duplicate(ms,
+                              comments=self.get_comments()+[util.execution_stamp('Geom','add_primitive')],
+                             )
 
 
     def mirror(self,directions,reflect=False):
@@ -503,15 +634,14 @@ class Geom:
         Parameters
         ----------
         directions : iterable containing str
-            direction(s) along which the microstructure is mirrored. Valid entries are 'x', 'y', 'z'.
+            Direction(s) along which the microstructure is mirrored.
+            Valid entries are 'x', 'y', 'z'.
         reflect : bool, optional
-            reflect (include) outermost layers.
+            Reflect (include) outermost layers. Defaults to False.
 
         """
-        valid = {'x','y','z'}
+        valid = ['x','y','z']
-        if not all(isinstance(d, str) for d in directions):
+        if not set(directions).issubset(valid):
-            raise TypeError('Directions are not of type str.')
-        elif not set(directions).issubset(valid):
             raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
 
         limits = [None,None] if reflect else [-2,0]
@@ -524,54 +654,88 @@ class Geom:
         if 'x' in directions:
             ms = np.concatenate([ms,ms[limits[0]:limits[1]:-1,:,:]],0)
 
-        #ToDo: self.add_comments('geom.py:mirror v{}'.format(version)
+        return self.duplicate(ms,
-        return self.update(ms,rescale=True)
+                              comments=self.get_comments()+[util.execution_stamp('Geom','mirror')],
+                              autosize=True)
 
 
-    def scale(self,grid):
+    def flip(self,directions):
+        """
+        Flip microstructure along given directions.
+
+        Parameters
+        ----------
+        directions : iterable containing str
+            Direction(s) along which the microstructure is flipped.
+            Valid entries are 'x', 'y', 'z'.
+
+        """
+        valid = ['x','y','z']
+        if not set(directions).issubset(valid):
+            raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
+
+        ms = np.flip(self.microstructure, (valid.index(d) for d in directions if d in valid))
+
+        return self.duplicate(ms,
+                              comments=self.get_comments()+[util.execution_stamp('Geom','flip')],
+                             )
+
+
+    def scale(self,grid,periodic=True):
         """
         Scale microstructure to new grid.
 
         Parameters
         ----------
         grid : numpy.ndarray of shape (3)
-            number of grid points in x,y,z direction.
+            Number of grid points in x,y,z direction.
+        periodic : Boolean, optional
+            Assume geometry to be periodic. Defaults to True.
 
         """
-        #ToDo: self.add_comments('geom.py:scale v{}'.format(version)
+        return self.duplicate(ndimage.interpolation.zoom(
-        return self.update(
+                                                         self.microstructure,
-                           ndimage.interpolation.zoom(
+                                                         grid/self.get_grid(),
-                                                      self.microstructure,
+                                                         output=self.microstructure.dtype,
-                                                      grid/self.get_grid(),
+                                                         order=0,
-                                                      output=self.microstructure.dtype,
+                                                         mode=('wrap' if periodic else 'nearest'),
-                                                      order=0,
+                                                         prefilter=False
-                                                      mode='nearest',
+                                                        ),
-                                                      prefilter=False
+                              comments=self.get_comments()+[util.execution_stamp('Geom','scale')],
-                                                     )
+                             )
-                          )
 
 
-    def clean(self,stencil=3):
+    def clean(self,stencil=3,selection=None,periodic=True):
         """
         Smooth microstructure by selecting most frequent index within given stencil at each location.
 
         Parameters
         ----------
         stencil : int, optional
-            size of smoothing stencil.
+            Size of smoothing stencil.
+        selection : list, optional
+            Field values that can be altered. Defaults to all.
+        periodic : Boolean, optional
+            Assume geometry to be periodic. Defaults to True.
 
         """
-        def mostFrequent(arr):
+        def mostFrequent(arr,selection=None):
-            unique, inverse = np.unique(arr, return_inverse=True)
+            me = arr[arr.size//2]
-            return unique[np.argmax(np.bincount(inverse))]
+            if selection is None or me in selection:
+                unique, inverse = np.unique(arr, return_inverse=True)
+                return unique[np.argmax(np.bincount(inverse))]
+            else:
+                return me
 
-        #ToDo: self.add_comments('geom.py:clean v{}'.format(version)
+        return self.duplicate(ndimage.filters.generic_filter(
-        return self.update(ndimage.filters.generic_filter(
+                                                             self.microstructure,
-                                                          self.microstructure,
+                                                             mostFrequent,
-                                                          mostFrequent,
+                                                             size=(stencil if selection is None else stencil//2*2+1,)*3,
-                                                          size=(stencil,)*3
+                                                             mode=('wrap' if periodic else 'nearest'),
-                                                         ).astype(self.microstructure.dtype)
+                                                             extra_keywords=dict(selection=selection),
-                          )
+                                                            ).astype(self.microstructure.dtype),
+                              comments=self.get_comments()+[util.execution_stamp('Geom','clean')],
+                             )
 
 
     def renumber(self):
@@ -580,8 +744,9 @@ class Geom:
         for i, oldID in enumerate(np.unique(self.microstructure)):
             renumbered = np.where(self.microstructure == oldID, i+1, renumbered)
 
-        #ToDo: self.add_comments('geom.py:renumber v{}'.format(version)
+        return self.duplicate(renumbered,
-        return self.update(renumbered)
+                              comments=self.get_comments()+[util.execution_stamp('Geom','renumber')],
+                             )
 
 
     def rotate(self,R,fill=None):
@@ -591,9 +756,9 @@ class Geom:
         Parameters
         ----------
         R : damask.Rotation
-            rotation to apply to the microstructure.
+            Rotation to apply to the microstructure.
         fill : int or float, optional
-            microstructure index to fill the corners. Defaults to microstructure.max() + 1.
+            Microstructure index to fill the corners. Defaults to microstructure.max() + 1.
 
         """
         if fill is None: fill = np.nanmax(self.microstructure) + 1
@@ -615,8 +780,11 @@ class Geom:
 
         origin = self.origin-(np.asarray(microstructure_in.shape)-self.grid)*.5 * self.size/self.grid
 
-        #ToDo: self.add_comments('geom.py:rotate v{}'.format(version)
+        return self.duplicate(microstructure_in,
-        return self.update(microstructure_in,origin=origin,rescale=True)
+                              origin=origin,
+                              comments=self.get_comments()+[util.execution_stamp('Geom','rotate')],
+                              autosize=True,
+                             )
 
 
     def canvas(self,grid=None,offset=None,fill=None):
@@ -626,46 +794,96 @@ class Geom:
         Parameters
         ----------
         grid : numpy.ndarray of shape (3)
-            number of grid points in x,y,z direction.
+            Number of grid points in x,y,z direction.
         offset : numpy.ndarray of shape (3)
-            offset (measured in grid points) from old to new microstructue[0,0,0].
+            Offset (measured in grid points) from old to new microstructure[0,0,0].
         fill : int or float, optional
-            microstructure index to fill the corners. Defaults to microstructure.max() + 1.
+            Microstructure index to fill the background. Defaults to microstructure.max() + 1.
 
         """
-        if fill is None: fill = np.nanmax(self.microstructure) + 1
         if offset is None: offset = 0
-        dtype = float if int(fill) != fill or self.microstructure.dtype==np.float else int
+        if fill is None: fill = np.nanmax(self.microstructure) + 1
+        dtype = float if int(fill) != fill or self.microstructure.dtype in np.sctypes['float'] else int
 
-        canvas = np.full(self.grid if grid is None else grid,
+        canvas = np.full(self.grid if grid is None else grid,fill,dtype)
-                         fill if fill is not None else np.nanmax(self.microstructure)+1,dtype)
 
-        l = np.clip( offset,          0,np.minimum(self.grid  +offset,grid))                        # noqa
+        LL = np.clip( offset,          0,np.minimum(self.grid,     grid+offset))
-        r = np.clip( offset+self.grid,0,np.minimum(self.grid*2+offset,grid))
+        UR = np.clip( offset+grid,     0,np.minimum(self.grid,     grid+offset))
-        L = np.clip(-offset,          0,np.minimum(grid       -offset,self.grid))
+        ll = np.clip(-offset,          0,np.minimum(     grid,self.grid-offset))
-        R = np.clip(-offset+grid,     0,np.minimum(grid*2     -offset,self.grid))
+        ur = np.clip(-offset+self.grid,0,np.minimum(     grid,self.grid-offset))
 
-        canvas[l[0]:r[0],l[1]:r[1],l[2]:r[2]] = self.microstructure[L[0]:R[0],L[1]:R[1],L[2]:R[2]]
+        canvas[ll[0]:ur[0],ll[1]:ur[1],ll[2]:ur[2]] = self.microstructure[LL[0]:UR[0],LL[1]:UR[1],LL[2]:UR[2]]
 
-        #ToDo: self.add_comments('geom.py:canvas v{}'.format(version)
+        return self.duplicate(canvas,
-        return self.update(canvas,origin=self.origin+offset*self.size/self.grid,rescale=True)
+                              origin=self.origin+offset*self.size/self.grid,
+                              comments=self.get_comments()+[util.execution_stamp('Geom','canvas')],
+                              autosize=True,
+                             )
 
 
     def substitute(self,from_microstructure,to_microstructure):
         """
-        Substitude microstructure indices.
+        Substitute microstructure indices.
 
         Parameters
         ----------
         from_microstructure : iterable of ints
-            microstructure indices to be substituted.
+            Microstructure indices to be substituted.
         to_microstructure : iterable of ints
-            new microstructure indices.
+            New microstructure indices.
 
         """
         substituted = self.get_microstructure()
         for from_ms,to_ms in zip(from_microstructure,to_microstructure):
             substituted[self.microstructure==from_ms] = to_ms
 
-        #ToDo: self.add_comments('geom.py:substitute v{}'.format(version)
+        return self.duplicate(substituted,
-        return self.update(substituted)
+                              comments=self.get_comments()+[util.execution_stamp('Geom','substitute')],
+                             )
+
+
+    def vicinity_offset(self,vicinity=1,offset=None,trigger=[],periodic=True):
+        """
+        Offset microstructure index of points in the vicinity of xxx.
+
+        Different from themselves (or listed as triggers) within a given (cubic) vicinity,
+        i.e. within the region close to a grain/phase boundary.
+        ToDo: use include/exclude as in seeds.from_geom
+
+        Parameters
+        ----------
+        vicinity : int, optional
+            Voxel distance checked for presence of other microstructure.
+            Defaults to 1.
+        offset : int, optional
+            Offset (positive or negative) to tag microstructure indices,
+            defaults to microstructure.max() + 1.
+        trigger : list of ints, optional
+            List of microstructure indices triggering a change.
+            Defaults to [], meaning that different neigboors trigger a change.
+        periodic : Boolean, optional
+            Assume geometry to be periodic. Defaults to True.
+
+        """
+        def tainted_neighborhood(stencil,trigger):
+
+            me = stencil[stencil.shape[0]//2]
+            if len(trigger) == 0:
+                return np.any(stencil != me)
+            if me in trigger:
+                trigger = set(trigger)
+                trigger.remove(me)
+                trigger = list(trigger)
+            return np.any(np.in1d(stencil,np.array(trigger)))
+
+        offset_ = np.nanmax(self.microstructure) if offset is None else offset
+        mask = ndimage.filters.generic_filter(self.microstructure,
+                                              tainted_neighborhood,
+                                              size=1+2*vicinity,
+                                              mode='wrap' if periodic else 'nearest',
+                                              extra_keywords={'trigger':trigger})
+        microstructure = np.ma.MaskedArray(self.microstructure + offset_, np.logical_not(mask))
+
+        return self.duplicate(microstructure,
+                              comments=self.get_comments()+[util.execution_stamp('Geom','vicinity_offset')],
+                             )

python/damask/_lattice.py

@@ -164,16 +164,16 @@ class Symmetry:
         with np.errstate(invalid='ignore'):
             # using '*'/prod for 'and'
             if self.system == 'cubic':
-                return np.where(np.prod(np.sqrt(2)-1. >= rho_abs,axis=-1) * \
+                return np.where(np.prod(np.sqrt(2)-1. >= rho_abs,axis=-1) *
                                 (1. >= np.sum(rho_abs,axis=-1)),True,False)
             elif self.system == 'hexagonal':
-                return np.where(np.prod(1.             >= rho_abs,axis=-1) * \
+                return np.where(np.prod(1.             >= rho_abs,axis=-1) *
-                                (2. >= np.sqrt(3)*rho_abs[...,0] + rho_abs[...,1]) * \
+                                (2. >= np.sqrt(3)*rho_abs[...,0] + rho_abs[...,1]) *
-                                (2. >= np.sqrt(3)*rho_abs[...,1] + rho_abs[...,0]) * \
+                                (2. >= np.sqrt(3)*rho_abs[...,1] + rho_abs[...,0]) *
                                 (2. >= np.sqrt(3) + rho_abs[...,2]),True,False)
             elif self.system == 'tetragonal':
-                return np.where(np.prod(1.             >= rho_abs[...,:2],axis=-1) * \
+                return np.where(np.prod(1.             >= rho_abs[...,:2],axis=-1) *
-                                (np.sqrt(2) >= rho_abs[...,0] + rho_abs[...,1]) * \
+                                (np.sqrt(2) >= rho_abs[...,0] + rho_abs[...,1]) *
                                 (np.sqrt(2) >= rho_abs[...,2] + 1.),True,False)
             elif self.system == 'orthorhombic':
                 return np.where(np.prod(1.             >= rho_abs,axis=-1),True,False)
@@ -321,6 +321,7 @@ class Lattice: # ToDo: Make a subclass of Symmetry!
     """
 
     lattices = {
+                'iso':      {'system':None},
                 'triclinic':{'system':None},
                 'bct':      {'system':'tetragonal'},
                 'hex':      {'system':'hexagonal'},

python/damask/_result.py

@@ -1,6 +1,5 @@
 import multiprocessing as mp
 import re
-import inspect
 import glob
 import os
 import datetime
@@ -49,7 +48,7 @@ class Result:
                 self.version_major = f.attrs['DADF5-major']
                 self.version_minor = f.attrs['DADF5-minor']
 
-            if self.version_major != 0 or not 2 <= self.version_minor <= 6:
+            if self.version_major != 0 or not 2 <= self.version_minor <= 7:
                 raise TypeError(f'Unsupported DADF5 version {self.version_major}.{self.version_minor}')
 
             self.structured = 'grid' in f['geometry'].attrs.keys()
@@ -401,8 +400,9 @@ class Result:
                             if sets is True:
                                 groups.append(group)
                             else:
-                                match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_]
+                                if group in f.keys():
-                                if len(set(match)) == len(sets): groups.append(group)
+                                    match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_]
+                                    if len(set(match)) == len(sets): groups.append(group)
         return groups
 
 
@@ -466,6 +466,11 @@ class Result:
             return f[self.get_dataset_location('orientation')[0]].attrs['Lattice'].astype('str')    # np.bytes_ to string
 
 
+    def enable_user_function(self,func):
+        globals()[func.__name__]=func
+        print(f'Function {func.__name__} enabled in add_calculation.')
+
+
     def read_dataset(self,path,c=0,plain=False):
         """
         Dataset for all points/cells.
@@ -504,7 +509,7 @@ class Result:
         else:
             return dataset
 
-
+    @property
     def cell_coordinates(self):
         """Return initial coordinates of the cell centers."""
         if self.structured:
@@ -513,6 +518,7 @@ class Result:
             with h5py.File(self.fname,'r') as f:
                 return f['geometry/x_c'][()]
 
+    @property
     def node_coordinates(self):
         """Return initial coordinates of the cell centers."""
         if self.structured:
@@ -530,7 +536,7 @@ class Result:
                 'meta':  {
                           'Unit':        x['meta']['Unit'],
                           'Description': f"Absolute value of {x['label']} ({x['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_absolute'
                           }
                  }
     def add_absolute(self,x):
@@ -558,10 +564,10 @@ class Result:
                 'meta':  {
                           'Unit':        kwargs['unit'],
                           'Description': f"{kwargs['description']} (formula: {kwargs['formula']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_calculation'
                           }
                  }
-    def add_calculation(self,label,formula,unit='n/a',description=None,vectorized=True):
+    def add_calculation(self,label,formula,unit='n/a',description=None):
         """
         Add result of a general formula.
 
@@ -575,13 +581,8 @@ class Result:
             Physical unit of the result.
         description : str, optional
             Human-readable description of the result.
-        vectorized : bool, optional
-            Indicate whether the formula can be used in vectorized form. Defaults to ‘True’.
 
         """
-        if not vectorized:
-            raise NotImplementedError
-
         dataset_mapping  = {d:d for d in set(re.findall(r'#(.*?)#',formula))}                       # datasets used in the formula
         args             = {'formula':formula,'label':label,'unit':unit,'description':description}
         self._add_generic_pointwise(self._add_calculation,dataset_mapping,args)
@@ -597,7 +598,7 @@ class Result:
                           'Description': "Cauchy stress calculated "
                                          f"from {P['label']} ({P['meta']['Description']})"
                                          f" and {F['label']} ({F['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_Cauchy'
                           }
                 }
     def add_Cauchy(self,P='P',F='F'):
@@ -623,7 +624,7 @@ class Result:
                 'meta':  {
                           'Unit':        T['meta']['Unit'],
                           'Description': f"Determinant of tensor {T['label']} ({T['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_determinant'
                           }
                 }
     def add_determinant(self,T):
@@ -647,7 +648,7 @@ class Result:
                 'meta':  {
                           'Unit':        T['meta']['Unit'],
                           'Description': f"Deviator of tensor {T['label']} ({T['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_deviator'
                           }
                  }
     def add_deviator(self,T):
@@ -678,7 +679,7 @@ class Result:
                 'meta' : {
                           'Unit':         T_sym['meta']['Unit'],
                           'Description': f"{label} eigenvalue of {T_sym['label']} ({T_sym['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_eigenvalue'
                          }
                 }
     def add_eigenvalue(self,T_sym,eigenvalue='max'):
@@ -711,7 +712,7 @@ class Result:
                           'Unit':        '1',
                           'Description': f"Eigenvector corresponding to {label} eigenvalue"
                                          f" of {T_sym['label']} ({T_sym['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_eigenvector'
                          }
                }
     def add_eigenvector(self,T_sym,eigenvalue='max'):
@@ -744,7 +745,7 @@ class Result:
                           'Unit':        '8-bit RGB',
                           'Lattice':     q['meta']['Lattice'],
                           'Description': 'Inverse Pole Figure (IPF) colors along sample direction [{} {} {}]'.format(*m),
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_IPF_color'
                          }
                }
     def add_IPF_color(self,q,l):
@@ -770,7 +771,7 @@ class Result:
                 'meta':  {
                           'Unit':        T_sym['meta']['Unit'],
                           'Description': f"Maximum shear component of {T_sym['label']} ({T_sym['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_maximum_shear'
                           }
                  }
     def add_maximum_shear(self,T_sym):
@@ -797,7 +798,7 @@ class Result:
                 'meta':  {
                           'Unit':        T_sym['meta']['Unit'],
                           'Description': f"Mises equivalent {t} of {T_sym['label']} ({T_sym['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_Mises'
                           }
                 }
     def add_Mises(self,T_sym):
@@ -833,7 +834,7 @@ class Result:
                 'meta':  {
                           'Unit':        x['meta']['Unit'],
                           'Description': f"{o}-norm of {t} {x['label']} ({x['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_norm'
                           }
                  }
     def add_norm(self,x,ord=None):
@@ -861,7 +862,7 @@ class Result:
                           'Description': "2. Piola-Kirchhoff stress calculated "
                                          f"from {P['label']} ({P['meta']['Description']})"
                                          f" and {F['label']} ({F['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_PK2'
                           }
                 }
     def add_PK2(self,P='P',F='F'):
@@ -897,7 +898,7 @@ class Result:
                           'Unit':        '1',
                           'Description': '{} coordinates of stereographic projection of pole (direction/plane) in crystal frame'\
                                          .format('Polar' if polar else 'Cartesian'),
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_pole'
                          }
                }
     def add_pole(self,q,p,polar=False):
@@ -925,7 +926,7 @@ class Result:
                 'meta':  {
                           'Unit':        F['meta']['Unit'],
                           'Description': f"Rotational part of {F['label']} ({F['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_rotational_part'
                           }
                  }
     def add_rotational_part(self,F):
@@ -949,7 +950,7 @@ class Result:
                 'meta':  {
                           'Unit':        T['meta']['Unit'],
                           'Description': f"Spherical component of tensor {T['label']} ({T['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_spherical'
                           }
                  }
     def add_spherical(self,T):
@@ -973,7 +974,7 @@ class Result:
                 'meta':  {
                           'Unit':        F['meta']['Unit'],
                           'Description': f"Strain tensor of {F['label']} ({F['meta']['Description']})",
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_strain_tensor'
                           }
                  }
     def add_strain_tensor(self,F='F',t='V',m=0.0):
@@ -1005,7 +1006,7 @@ class Result:
                           'Unit':        F['meta']['Unit'],
                           'Description': '{} stretch tensor of {} ({})'.format('Left' if t.upper() == 'V' else 'Right',
                                                                                F['label'],F['meta']['Description']),
-                          'Creator':     inspect.stack()[0][3][1:]
+                          'Creator':     'add_stretch_tensor'
                           }
                  }
     def add_stretch_tensor(self,F='F',t='V'):
@@ -1064,6 +1065,10 @@ class Result:
         lock = mp.Manager().Lock()
 
         groups = self.groups_with_datasets(datasets.values())
+        if len(groups) == 0:
+            print('No matching dataset found, no data was added.')
+            return
+
         default_arg = partial(self._job,func=func,datasets=datasets,args=args,lock=lock)
 
         for result in util.show_progress(pool.imap_unordered(default_arg,groups),len(groups)):
@@ -1263,4 +1268,4 @@ class Result:
             u = self.read_dataset(self.get_dataset_location('u_n' if mode.lower() == 'cell' else 'u_p'))
             v.add(u,'u')
 
-            v.write(f'{self.fname.stem}_inc{inc[3:].zfill(N_digits)}')
+            v.to_file(f'{self.fname.stem}_inc{inc[3:].zfill(N_digits)}')

python/damask/_rotation.py

@@ -357,7 +357,7 @@ class Rotation:
         accept_homomorph : boolean, optional
             Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere).
             Defaults to False.
-        P : integer ∈ {-1,1}, optional
+        P : int ∈ {-1,1}, optional
             Convention used. Defaults to -1.
 
         """
@@ -422,7 +422,7 @@ class Rotation:
             Angle ω is given in degrees. Defaults to False.
         normalize: boolean, optional
             Allow |n| ≠ 1. Defaults to False.
-        P : integer ∈ {-1,1}, optional
+        P : int ∈ {-1,1}, optional
             Convention used. Defaults to -1.
 
         """
@@ -505,7 +505,7 @@ class Rotation:
             (n_1, n_2, n_3, tan(ω/2)), |n| = 1  and ω ∈ [0,π].
         normalize : boolean, optional
             Allow |n| ≠ 1. Defaults to False.
-        P : integer ∈ {-1,1}, optional
+        P : int ∈ {-1,1}, optional
             Convention used. Defaults to -1.
 
         """
@@ -534,7 +534,7 @@ class Rotation:
         ----------
         h : numpy.ndarray of shape (...,3)
             Homochoric vector: (h_1, h_2, h_3), |h| < (3/4*π)^(1/3).
-        P : integer ∈ {-1,1}, optional
+        P : int ∈ {-1,1}, optional
             Convention used. Defaults to -1.
 
         """
@@ -561,7 +561,7 @@ class Rotation:
         ----------
         c : numpy.ndarray of shape (...,3)
             Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3).
-        P : integer ∈ {-1,1}, optional
+        P : int ∈ {-1,1}, optional
             Convention used. Defaults to -1.
 
         """

python/damask/_table.py

@@ -3,7 +3,6 @@ import re
 import pandas as pd
 import numpy as np
 
-from . import version
 from . import util
 
 class Table:
@@ -49,7 +48,9 @@ class Table:
 
     def _add_comment(self,label,shape,info):
         if info is not None:
-            self.comments.append(f'{label}{" "+str(shape) if np.prod(shape,dtype=int) > 1 else ""}: {info}')
+            specific = f'{label}{" "+str(shape) if np.prod(shape,dtype=int) > 1 else ""}: {info}'
+            general  = util.execution_stamp('Table')
+            self.comments.append(f'{specific} / {general}')
 
 
     @staticmethod
@@ -135,7 +136,7 @@ class Table:
 
         content = f.readlines()
 
-        comments = [f'table.py:from_ang v{version}']
+        comments = [util.execution_stamp('Table','from_ang')]
         for line in content:
             if line.startswith('#'):
                 comments.append(line.strip())
@@ -312,7 +313,7 @@ class Table:
                 self.shapes[key] = other.shapes[key]
 
 
-    def to_ASCII(self,fname,new_style=False):
+    def to_file(self,fname,format='ASCII',new_style=False):
         """
         Store as plain text file.
 
@@ -320,32 +321,51 @@ class Table:
         ----------
         fname : file, str, or pathlib.Path
             Filename or file for writing.
+        format : {ASCII'}, optional
+            File format, defaults to 'ASCII'. Available formats are:
+            - ASCII: Plain text file, extension '.txt'.
         new_style : Boolean, optional
             Write table in new style, indicating header lines by comment sign ('#') only.
 
         """
-        seen = set()
+        def _to_ASCII(table,fname,new_style=False):
-        labels = []
+            """
-        for l in [x for x in self.data.columns if not (x in seen or seen.add(x))]:
+            Store as plain text file.
-            if self.shapes[l] == (1,):
-                labels.append(f'{l}')
-            elif len(self.shapes[l]) == 1:
-                labels += [f'{i+1}_{l}' \
-                          for i in range(self.shapes[l][0])]
-            else:
-                labels += [f'{util.srepr(self.shapes[l],"x")}:{i+1}_{l}' \
-                          for i in range(np.prod(self.shapes[l]))]
 
-        if new_style:
+            Parameters
-            header = [f'# {comment}' for comment in self.comments]
+            ----------
+            table : Table object
+                Table to write.
+            fname : file, str, or pathlib.Path
+                Filename or file for writing.
+            new_style : Boolean, optional
+                Write table in new style, indicating header lines by comment sign ('#') only.
+
+            """
+            seen = set()
+            labels = []
+            for l in [x for x in table.data.columns if not (x in seen or seen.add(x))]:
+                if table.shapes[l] == (1,):
+                    labels.append(f'{l}')
+                elif len(table.shapes[l]) == 1:
+                    labels += [f'{i+1}_{l}' \
+                              for i in range(table.shapes[l][0])]
+                else:
+                    labels += [f'{util.srepr(table.shapes[l],"x")}:{i+1}_{l}' \
+                              for i in range(np.prod(table.shapes[l]))]
+
+            header = [f'# {comment}' for comment in table.comments] if new_style else \
+                     [f'{len(table.comments)+1} header'] + table.comments
+
+            try:
+                f = open(fname,'w')
+            except TypeError:
+                f = fname
+
+            for line in header + [' '.join(labels)]: f.write(line+'\n')
+            table.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False)
+
+        if format.lower() == 'ascii':
+            return _to_ASCII(self,fname,new_style)
         else:
-            header = [f'{len(self.comments)+1} header'] \
+            raise TypeError(f'Unknown format {format}.')
-                   + self.comments \
-
-        try:
-            f = open(fname,'w')
-        except TypeError:
-            f = fname
-
-        for line in header + [' '.join(labels)]: f.write(line+'\n')
-        self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False)

python/damask/_vtk.py

@@ -6,8 +6,10 @@ import numpy as np
 import vtk
 from vtk.util.numpy_support import numpy_to_vtk            as np_to_vtk
 from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray
+from vtk.util.numpy_support import vtk_to_numpy            as vtk_to_np
 
-import damask
+from . import util
+from . import environment
 from . import Table
 
 
@@ -18,18 +20,19 @@ class VTK:
     High-level interface to VTK.
     """
 
-    def __init__(self,geom):
+    def __init__(self,vtk_data):
         """
-        Set geometry and topology.
+        Initialize from vtk dataset.
 
         Parameters
         ----------
-        geom : subclass of vtk.vtkDataSet
+        vtk_data : subclass of vtk.vtkDataSet
-            Description of geometry and topology. Valid types are vtk.vtkRectilinearGrid,
+            Description of geometry and topology, optionally with attached data.
-            vtk.vtkUnstructuredGrid, or vtk.vtkPolyData.
+            Valid types are vtk.vtkRectilinearGrid, vtk.vtkUnstructuredGrid,
+            or vtk.vtkPolyData.
 
         """
-        self.geom = geom
+        self.vtk_data = vtk_data
 
 
     @staticmethod
@@ -49,15 +52,15 @@ class VTK:
             Spatial origin.
 
         """
-        geom = vtk.vtkRectilinearGrid()
+        vtk_data = vtk.vtkRectilinearGrid()
-        geom.SetDimensions(*(grid+1))
+        vtk_data.SetDimensions(*(grid+1))
         coord = [np_to_vtk(np.linspace(origin[i],origin[i]+size[i],grid[i]+1),deep=True) for i in [0,1,2]]
         [coord[i].SetName(n) for i,n in enumerate(['x','y','z'])]
-        geom.SetXCoordinates(coord[0])
+        vtk_data.SetXCoordinates(coord[0])
-        geom.SetYCoordinates(coord[1])
+        vtk_data.SetYCoordinates(coord[1])
-        geom.SetZCoordinates(coord[2])
+        vtk_data.SetZCoordinates(coord[2])
 
-        return VTK(geom)
+        return VTK(vtk_data)
 
 
     @staticmethod
@@ -85,11 +88,11 @@ class VTK:
                             connectivity),axis=1).ravel()
         cells.SetCells(connectivity.shape[0],np_to_vtkIdTypeArray(T,deep=True))
 
-        geom = vtk.vtkUnstructuredGrid()
+        vtk_data = vtk.vtkUnstructuredGrid()
-        geom.SetPoints(vtk_nodes)
+        vtk_data.SetPoints(vtk_nodes)
-        geom.SetCells(eval(f'vtk.VTK_{cell_type.split("_",1)[-1].upper()}'),cells)
+        vtk_data.SetCells(eval(f'vtk.VTK_{cell_type.split("_",1)[-1].upper()}'),cells)
 
-        return VTK(geom)
+        return VTK(vtk_data)
 
 
     @staticmethod
@@ -108,10 +111,10 @@ class VTK:
         vtk_points = vtk.vtkPoints()
         vtk_points.SetData(np_to_vtk(points))
 
-        geom = vtk.vtkPolyData()
+        vtk_data = vtk.vtkPolyData()
-        geom.SetPoints(vtk_points)
+        vtk_data.SetPoints(vtk_points)
 
-        return VTK(geom)
+        return VTK(vtk_data)
 
 
     @staticmethod
@@ -129,18 +132,20 @@ class VTK:
 
         """
         ext = Path(fname).suffix
-        if ext == '.vtk' or dataset_type:
+        if ext == '.vtk' or dataset_type is not None:
             reader = vtk.vtkGenericDataObjectReader()
             reader.SetFileName(str(fname))
-            reader.Update()
             if dataset_type is None:
                 raise TypeError('Dataset type for *.vtk file not given.')
             elif dataset_type.lower().endswith('rectilineargrid'):
-                geom = reader.GetRectilinearGridOutput()
+                reader.Update()
+                vtk_data = reader.GetRectilinearGridOutput()
             elif dataset_type.lower().endswith('unstructuredgrid'):
-                geom = reader.GetUnstructuredGridOutput()
+                reader.Update()
+                vtk_data = reader.GetUnstructuredGridOutput()
             elif dataset_type.lower().endswith('polydata'):
-                geom = reader.GetPolyDataOutput()
+                reader.Update()
+                vtk_data = reader.GetPolyDataOutput()
             else:
                 raise TypeError(f'Unknown dataset type {dataset_type} for vtk file')
         else:
@@ -155,15 +160,15 @@ class VTK:
 
             reader.SetFileName(str(fname))
             reader.Update()
-            geom = reader.GetOutput()
+            vtk_data = reader.GetOutput()
 
-        return VTK(geom)
+        return VTK(vtk_data)
 
     @staticmethod
     def _write(writer):
         """Wrapper for parallel writing."""
         writer.Write()
-    def write(self,fname,parallel=True):
+    def to_file(self,fname,parallel=True):
         """
         Write to file.
 
@@ -175,11 +180,11 @@ class VTK:
             Write data in parallel background process. Defaults to True.
 
         """
-        if   isinstance(self.geom,vtk.vtkRectilinearGrid):
+        if   isinstance(self.vtk_data,vtk.vtkRectilinearGrid):
             writer = vtk.vtkXMLRectilinearGridWriter()
-        elif isinstance(self.geom,vtk.vtkUnstructuredGrid):
+        elif isinstance(self.vtk_data,vtk.vtkUnstructuredGrid):
             writer = vtk.vtkXMLUnstructuredGridWriter()
-        elif isinstance(self.geom,vtk.vtkPolyData):
+        elif isinstance(self.vtk_data,vtk.vtkPolyData):
             writer = vtk.vtkXMLPolyDataWriter()
 
         default_ext = writer.GetDefaultFileExtension()
@@ -189,10 +194,14 @@ class VTK:
         writer.SetFileName(str(Path(fname).with_suffix('.'+default_ext)))
         writer.SetCompressorTypeToZLib()
         writer.SetDataModeToBinary()
-        writer.SetInputData(self.geom)
+        writer.SetInputData(self.vtk_data)
+
         if parallel:
-            mp_writer = mp.Process(target=self._write,args=(writer,))
+            try:
-            mp_writer.start()
+                mp_writer = mp.Process(target=self._write,args=(writer,))
+                mp_writer.start()
+            except TypeError:
+                writer.Write()
         else:
             writer.Write()
 
@@ -200,9 +209,20 @@ class VTK:
     # Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data
     # Needs support for pd.DataFrame and/or table
     def add(self,data,label=None):
-        """Add data to either cells or points."""
+        """
-        N_points = self.geom.GetNumberOfPoints()
+        Add data to either cells or points.
-        N_cells  = self.geom.GetNumberOfCells()
+
+        Parameters
+        ----------
+        data : numpy.ndarray
+            Data to add. First dimension need to match either
+            number of cells or number of points
+        label : str
+            Data label.
+
+        """
+        N_points = self.vtk_data.GetNumberOfPoints()
+        N_cells  = self.vtk_data.GetNumberOfCells()
 
         if   isinstance(data,np.ndarray):
             if label is None:
@@ -215,9 +235,9 @@ class VTK:
             d.SetName(label)
 
             if   data.shape[0] == N_cells:
-                self.geom.GetCellData().AddArray(d)
+                self.vtk_data.GetCellData().AddArray(d)
             elif data.shape[0] == N_points:
-                self.geom.GetPointData().AddArray(d)
+                self.vtk_data.GetPointData().AddArray(d)
             else:
                 raise ValueError(f'Invalid shape {data.shape[0]}')
         elif isinstance(data,pd.DataFrame):
@@ -228,12 +248,79 @@ class VTK:
             raise TypeError
 
 
+    def get(self,label):
+        """
+        Get either cell or point data.
+
+        Cell data takes precedence over point data, i.e. this
+        function assumes that labels are unique among cell and
+        point data.
+
+        Parameters
+        ----------
+        label : str
+            Data label.
+
+        """
+        cell_data = self.vtk_data.GetCellData()
+        for a in range(cell_data.GetNumberOfArrays()):
+            if cell_data.GetArrayName(a) == label:
+                return vtk_to_np(cell_data.GetArray(a))
+
+        point_data = self.vtk_data.GetPointData()
+        for a in range(point_data.GetNumberOfArrays()):
+            if point_data.GetArrayName(a) == label:
+                return vtk_to_np(point_data.GetArray(a))
+
+        raise ValueError(f'array "{label}" not found')
+
+
+    def get_comments(self):
+        """Return the comments."""
+        fielddata = self.vtk_data.GetFieldData()
+        for a in range(fielddata.GetNumberOfArrays()):
+            if fielddata.GetArrayName(a) == 'comments':
+                comments = fielddata.GetAbstractArray(a)
+                return [comments.GetValue(i) for i in range(comments.GetNumberOfValues())]
+        return []
+
+
+    def set_comments(self,comments):
+        """
+        Set comments.
+
+        Parameters
+        ----------
+        comments : str or list of str
+            Comments.
+
+        """
+        s = vtk.vtkStringArray()
+        s.SetName('comments')
+        for c in [comments] if isinstance(comments,str) else comments:
+            s.InsertNextValue(c)
+        self.vtk_data.GetFieldData().AddArray(s)
+
+
+    def add_comments(self,comments):
+        """
+        Add comments.
+
+        Parameters
+        ----------
+        comments : str or list of str
+            Comments to add.
+
+        """
+        self.set_comments(self.get_comments() + ([comments] if isinstance(comments,str) else comments))
+
+
     def __repr__(self):
         """ASCII representation of the VTK data."""
         writer = vtk.vtkDataSetWriter()
-        writer.SetHeader(f'# damask.VTK v{damask.version}')
+        writer.SetHeader(f'# {util.execution_stamp("VTK")}')
         writer.WriteToOutputStringOn()
-        writer.SetInputData(self.geom)
+        writer.SetInputData(self.vtk_data)
         writer.Write()
         return writer.GetOutputString()
 
@@ -245,7 +332,7 @@ class VTK:
         See http://compilatrix.com/article/vtk-1 for further ideas.
         """
         mapper = vtk.vtkDataSetMapper()
-        mapper.SetInputData(self.geom)
+        mapper.SetInputData(self.vtk_data)
         actor = vtk.vtkActor()
         actor.SetMapper(mapper)
 
@@ -257,7 +344,7 @@ class VTK:
         ren.AddActor(actor)
         ren.SetBackground(0.2,0.2,0.2)
 
-        window.SetSize(damask.environment.screen_size[0],damask.environment.screen_size[1])
+        window.SetSize(environment.screen_size[0],environment.screen_size[1])
 
         iren = vtk.vtkRenderWindowInteractor()
         iren.SetRenderWindow(window)

python/damask/config/__init__.py

@@ -1,3 +0,0 @@
-"""Aggregator for configuration file handling."""
-
-from .material import Material    # noqa

python/damask/config/material.py

@@ -1,282 +0,0 @@
-import re
-import os
-
-from damask import util
-
-class Section():
-  def __init__(self,data = {'__order__':[]},part = ''):
-    """New material.config section."""
-    classes = {
-                'homogenization':Homogenization,
-                'microstructure':Microstructure,
-                'crystallite':Crystallite,
-                'phase':Phase,
-                'texture':Texture,
-              }
-    self.parameters = {}
-    for key in data:
-      self.parameters[key] = data[key] if isinstance(data[key], list) else [data[key]]
-
-    if '__order__' not in self.parameters:
-      self.parameters['__order__'] = list(self.parameters.keys())
-    if part.lower() in classes:
-      self.__class__ = classes[part.lower()]
-      self.__init__(data)
-
-  def add_multiKey(self,key,data):
-    multiKey = '(%s)'%key
-    if multiKey not in self.parameters: self.parameters[multiKey] = []
-    if multiKey not in self.parameters['__order__']: self.parameters['__order__'] += [multiKey]
-    self.parameters[multiKey] += [[item] for item in data] if isinstance(data, list) else [[data]]
-
-  def data(self):
-    return self.parameters
-
-
-class Homogenization(Section):
-  def __init__(self,data = {'__order__':[]}):
-    """New material.config <homogenization> section."""
-    Section.__init__(self,data)
-
-
-class Crystallite(Section):
-  def __init__(self,data = {'__order__':[]}):
-    """New material.config <crystallite> section."""
-    Section.__init__(self,data)
-
-
-class Phase(Section):
-  def __init__(self,data = {'__order__':[]}):
-    """New material.config <Phase> section."""
-    Section.__init__(self,data)
-
-
-class Microstructure(Section):
-  def __init__(self,data = {'__order__':[]}):
-    """New material.config <microstructure> section."""
-    Section.__init__(self,data)
-
-
-class Texture(Section):
-  def __init__(self,data = {'__order__':[]}):
-    """New material.config <texture> section."""
-    Section.__init__(self,data)
-
-  def add_component(self,theType,properties):
-
-    scatter  = properties['scatter']  if 'scatter'  in list(map(str.lower,list(properties.keys()))) else 0.0
-    fraction = properties['fraction'] if 'fraction' in list(map(str.lower,list(properties.keys()))) else 1.0
-
-    try:
-      multiKey = theType.lower()
-    except AttributeError:
-      pass
-
-    if multiKey == 'gauss':
-      self.add_multiKey(multiKey,'phi1 %g\tPhi %g\tphi2 %g\tscatter %g\tfraction %g'%(
-                                        properties['eulers'][0],
-                                        properties['eulers'][1],
-                                        properties['eulers'][2],
-                                        scatter,
-                                        fraction,
-                                        )
-                        )
-
-
-class Material():
-  """Read, manipulate, and write material.config files."""
-
-  def __init__(self,verbose=True):
-    """Generates ordered list of parts."""
-    self.parts = [
-            'homogenization',
-            'crystallite',
-            'phase',
-            'texture',
-            'microstructure',
-            ]
-    self.data = {
-            'homogenization': {'__order__': []},
-            'microstructure': {'__order__': []},
-            'crystallite':    {'__order__': []},
-            'phase':          {'__order__': []},
-            'texture':        {'__order__': []},
-            }
-    self.verbose = verbose
-
-  def __repr__(self):
-    """Returns current data structure in material.config format."""
-    me = []
-    for part in self.parts:
-      if self.verbose: print(f'processing <{part}>')
-      me += ['',
-             '#'*100,
-             f'<{part}>',
-             '#'*100,
-            ]
-      for section in self.data[part]['__order__']:
-        me += [f'[{section}] {"#"+"-"*max(0,96-len(section))}']
-        for key in self.data[part][section]['__order__']:
-          if key.startswith('(') and key.endswith(')'):                       # multiple (key)
-            me += [f'{key}\t{" ".join(values)}' for values in self.data[part][section][key]]
-          else:                                                               # plain key
-            me += [f'{key}\t{util.srepr(self.data[part][section][key]," ")}']
-    return '\n'.join(me) + '\n'
-
-  def parse(self, part=None, sections=[], content=None):
-
-    re_part = re.compile(r'^<(.+)>$')                     # pattern for part
-    re_sec  = re.compile(r'^\[(.+)\]$')                   # pattern for section
-
-    name_section = ''
-    active = False
-
-    for line in content:
-      line = line.split('#')[0].strip()                   # kill comments and extra whitespace
-      line = line.split('/echo/')[0].strip()              # remove '/echo/' tags
-      line = line.lower()                                 # be case insensitive
-
-      if line:                                            # content survives...
-        match_part = re_part.match(line)
-        if match_part:                                    # found <...> separator
-          active = (match_part.group(1) == part)          # only active in <part>
-          continue
-        if active:
-          match_sec  = re_sec.match(line)
-          if match_sec:                                   # found [section]
-            name_section = match_sec.group(1)             # remember name ...
-            if '__order__' not in self.data[part]: self.data[part]['__order__'] = []
-            self.data[part]['__order__'].append(name_section)  # ... and position
-            self.data[part][name_section] = {'__order__':[]}
-            continue
-
-          if sections == [] or name_section in sections:  # possibly restrict to subset
-            items = line.split()
-            if items[0] not in self.data[part][name_section]:         # first encounter of key?
-              self.data[part][name_section][items[0]] = []            # create item
-              self.data[part][name_section]['__order__'].append(items[0])
-            if items[0].startswith('(') and items[0].endswith(')'):   # multiple "(key)"
-              self.data[part][name_section][items[0]].append(items[1:])
-            else:                                                     # plain key
-              self.data[part][name_section][items[0]] = items[1:]
-
-
-
-  def read(self,filename=None):
-    """Read material.config file."""
-    def recursiveRead(filename):
-      """Takes care of include statements like '{}'."""
-      result = []
-      re_include  = re.compile(r'^{(.+)}$')
-      with open(filename) as f: lines = f.readlines()
-      for line in lines:
-        match = re_include.match(line.split()[0]) if line.strip() else False
-        result += [line] if not match else \
-                  recursiveRead(match.group(1) if match.group(1).startswith('/') else
-                                os.path.normpath(os.path.join(os.path.dirname(filename),match.group(1))))
-      return result
-
-    c = recursiveRead(filename)
-    for p in self.parts:
-      self.parse(part=p, content=c)
-
-  def write(self,filename='material.config', overwrite=False):
-    """Write to material.config."""
-    i = 0
-    outname = filename
-    while os.path.exists(outname) and not overwrite:
-     i += 1
-     outname = f'{filename}_{i}'
-
-    if self.verbose: print(f'Writing material data to {outname}')
-    with open(outname,'w') as f:
-      f.write(str(self))
-    return outname
-
-  def add_section(self, part=None, section=None, initialData=None, merge=False):
-    """Add Update."""
-    part    = part.lower()
-    section = section.lower()
-    if part not in self.parts: raise Exception(f'invalid part {part}')
-
-    if not isinstance(initialData, dict):
-      initialData = initialData.data()
-
-    if section not in self.data[part]: self.data[part]['__order__'] += [section]
-    if section in self.data[part] and merge:
-      for existing in self.data[part][section]['__order__']:                               # replace existing
-        if existing in initialData['__order__']:
-          if existing.startswith('(') and existing.endswith(')'):                          # multiple (key)
-            self.data[part][section][existing] += initialData[existing]                    # add new multiple entries to existing ones
-          else:                                                                            # regular key
-            self.data[part][section][existing] = initialData[existing]                     # plain replice
-      for new in initialData['__order__']:                                                 # merge new content
-        if new not in self.data[part][section]['__order__']:
-          self.data[part][section][new] = initialData[new]
-          self.data[part][section]['__order__'] += [new]
-    else:
-      self.data[part][section] = initialData
-
-
-
-
-  def add_microstructure(self, section='',
-                               components={},      # dict of phase,texture, and fraction lists
-                        ):
-    """Experimental! Needs expansion to multi-constituent microstructures..."""
-    microstructure = Microstructure()
-    # make keys lower case (http://stackoverflow.com/questions/764235/dictionary-to-lowercase-in-python)
-    components=dict((k.lower(), v) for k,v in components.items())
-
-    for key in ['phase','texture','fraction','crystallite']:
-      if isinstance(components[key], list):
-        for i, x in enumerate(components[key]):
-          try:
-            components[key][i] = x.lower()
-          except AttributeError:
-            pass
-      else:
-        try:
-          components[key] = [components[key].lower()]
-        except AttributeError:
-          components[key] = [components[key]]
-
-    for (phase,texture,fraction,crystallite) in zip(components['phase'],components['texture'],
-                                                    components['fraction'],components['crystallite']):
-      microstructure.add_multiKey('constituent','phase %i\ttexture %i\tfraction %g\ncrystallite %i'%(
-                                    self.data['phase']['__order__'].index(phase)+1,
-                                    self.data['texture']['__order__'].index(texture)+1,
-                                    fraction,
-                                    self.data['crystallite']['__order__'].index(crystallite)+1))
-
-    self.add_section('microstructure',section,microstructure)
-
-
-  def change_value(self, part=None,
-                         section=None,
-                         key=None,
-                         value=None):
-    if not isinstance(value,list):
-      if not isinstance(value,str):
-        value = '%s'%value
-      value = [value]
-    newlen = len(value)
-    oldval = self.data[part.lower()][section.lower()][key.lower()]
-    oldlen = len(oldval)
-    print('changing %s:%s:%s from %s to %s '%(part.lower(),section.lower(),key.lower(),oldval,value))
-    self.data[part.lower()][section.lower()][key.lower()] = value
-    if newlen is not oldlen:
-      print('Length of value was changed from %i to %i!'%(oldlen,newlen))
-
-
-  def add_value(self, part=None,
-                  section=None,
-                  key=None,
-                  value=None):
-    if not isinstance(value,list):
-      if not isinstance(value,str):
-        value = '%s'%value
-      value = [value]
-    print('adding %s:%s:%s with value %s '%(part.lower(),section.lower(),key.lower(),value))
-    self.data[part.lower()][section.lower()][key.lower()] = value
-    self.data[part.lower()][section.lower()]['__order__'] += [key.lower()]

python/damask/solver/_marc.py

@@ -3,12 +3,12 @@ import shlex
 import string
 from pathlib import Path
 
-import damask
+from .. import environment
 
 class Marc:
     """Wrapper to run DAMASK with MSCMarc."""
 
-    def __init__(self,version=damask.environment.options['MARC_VERSION']):
+    def __init__(self,version=environment.options['MARC_VERSION']):
         """
         Create a Marc solver object.
 
@@ -24,7 +24,7 @@ class Marc:
     @property
     def library_path(self):
 
-        path_MSC = damask.environment.options['MSC_ROOT']
+        path_MSC = environment.options['MSC_ROOT']
         path_lib = Path(f'{path_MSC}/mentat{self.version}/shlib/linux64')
 
         return path_lib if path_lib.is_dir() else None
@@ -33,7 +33,7 @@ class Marc:
     @property
     def tools_path(self):
 
-        path_MSC   = damask.environment.options['MSC_ROOT']
+        path_MSC   = environment.options['MSC_ROOT']
         path_tools = Path(f'{path_MSC}/marc{self.version}/tools')
 
         return path_tools if path_tools.is_dir() else None
@@ -49,7 +49,7 @@ class Marc:
                   ):
 
 
-        usersub = damask.environment.root_dir/'src/DAMASK_marc'
+        usersub = environment.root_dir/'src/DAMASK_marc'
         usersub = usersub.parent/(usersub.name + ('.f90' if compile else '.marc'))
         if not usersub.is_file():
             raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),usersub))

python/damask/util.py

@@ -9,6 +9,8 @@ from optparse import Option
 
 import numpy as np
 
+from . import version
+
 # limit visibility
 __all__=[
          'srepr',
@@ -20,6 +22,7 @@ __all__=[
          'scale_to_coprime',
          'return_message',
          'extendableOption',
+         'execution_stamp'
         ]
 
 ####################################################################################################
@@ -175,6 +178,13 @@ def scale_to_coprime(v):
     return m
 
 
+def execution_stamp(class_name,function_name=None):
+    """Timestamp the execution of a (function within a) class."""
+    now = datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
+    _function_name = '' if function_name is None else f'.{function_name}'
+    return f'damask.{class_name}{_function_name} v{version} ({now})'
+
+
 ####################################################################################################
 # Classes
 ####################################################################################################

python/setup.py

@@ -16,10 +16,11 @@ setuptools.setup(
     packages=setuptools.find_packages(),
     include_package_data=True,
     install_requires = [
-        "pandas",
+        "pandas",                                                                                   # requires numpy
         "scipy",
-        "h5py",
+        "h5py",                                                                                     # requires numpy
         "vtk",
+        "matplotlib",                                                                               # requires numpy, pillow
     ],
     classifiers = [
         "Intended Audience :: Science/Research",

python/tests/conftest.py

@@ -1,11 +1,38 @@
 from pathlib import Path
-import numpy as np
+import datetime
 
+import numpy as np
 import pytest
 
-# Use to monkeypatch damask.version (for comparsion to reference files that contain version information)
+import damask
-def pytest_configure():
+
-    pytest.dummy_version = '99.99.99-9999-pytest'
+
+patched_version = '99.99.99-9999-pytest'
+@pytest.fixture
+def patch_damask_version(monkeypatch):
+    """Set damask.version for reproducible tests results."""
+    monkeypatch.setattr(damask, 'version', patched_version)
+
+
+patched_date = datetime.datetime(2019, 11, 2, 11, 58, 0)
+@pytest.fixture
+def patch_datetime_now(monkeypatch):
+    """Set datetime.datetime.now for reproducible tests results."""
+    class mydatetime:
+        @classmethod
+        def now(cls):
+            return patched_date
+
+    monkeypatch.setattr(datetime, 'datetime', mydatetime)
+
+@pytest.fixture
+def execution_stamp(monkeypatch):
+    """Set damask.util.execution_stamp for reproducible tests results."""
+    def execution_stamp(class_name,function_name=None):
+        _function_name = '' if function_name is None else f'.{function_name}'
+        return f'damask.{class_name}{_function_name} v{patched_version} ({patched_date})'
+
+    monkeypatch.setattr(damask.util, 'execution_stamp', execution_stamp)
 
 
 def pytest_addoption(parser):

python/tests/reference/Colormap/binary.json

@@ -1,6 +1,6 @@
 [
     {
-        "Creator": "damask.Colormap v99.99.99-9999-pytest",
+        "Creator": "damask.Colormap v99.99.99-9999-pytest (2019-11-02 11:58:00)",
         "ColorSpace": "RGB",
         "Name": "binary",
         "DefaultMap": true,

python/tests/reference/Colormap/binary.txt

@@ -1,4 +1,4 @@
-# Creator: damask.Colormap v99.99.99-9999-pytest
+# Creator: damask.Colormap v99.99.99-9999-pytest (2019-11-02 11:58:00)
 1_RGBA 2_RGBA 3_RGBA 4_RGBA
 1.0 1.0 1.0 1.0
 0.996078431372549 0.996078431372549 0.996078431372549 1.0

python/tests/reference/Geom/clean_1_1+2+3_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_1_1+2+3_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_1_1_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_1_1_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_1_None_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_1_None_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_2_1+2+3_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAZwAAAA==eF7t0rcOgmAAhVEgNmyo2AuoWN//BR04EwsJcfzvcvabL47qxcFOJg177HPAIUdMOeaEU844Z8YFl1wx55obbrnjngceeeKZFxYseeWNd1Z88MkX3/zwy+Z/wf8YOqzX1uEPlgwHCA==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_2_1+2+3_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAagAAAA==eF7t0rkOglAARFExLrgCKuKuqLj8/w9acCoSY7B+05x+cqNOvSj4l92GPfY54JAxRxxzwilnnDNhyowLLrlizjULbrjljnseeOSJZ15Y8sob76z44JMvvtn8L9jObz2GDuv96vADk5QHBg==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_2_1_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAZAAAAA==eF7t0scRglAAQEEBAyZUMCuomPtv0ANbgMNw/O+yDbyo1xQFWxkzYZ8DDjliyjEnnHLGOTMuuOSKOQuuueGWO+554JEnnlmy4oVX3ljzzgeffPHND7+Mg50aPmz698MfmvQHCg==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_2_1_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_2_2_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 20 0 18 0 8">
+  <Piece Extent="0 20 0 18 0 8">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int32" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAALQAAOgAAAA==eF7t1CEOACAMBMHS/z8aXwNJSagYMe6y8jIislgNtTW9d9oD/PL6r6b3AAAAAAAAAAAA4MYGlRYLYA==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACoAAAAVQAAAA==eF5jYICAWTNBYKU9hN5pb2IMAoeh/JP2EFUXoOKX7dPTQOAaVP6m/dkzIHAHqu4BVPwhVP1jqPwTqL5nUHUvoOpeQtW9hqp7A1X3Dqrugz0ASSZF3Q==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACYAAAAVgAAAA==eF5jYIAAmeOFQLTGHkLvsQ8Fg6NQ/hn7IjDjIlT8qr1F32MgugGVv2MPMeUBVN1D+8cQBVD1T+3BymSeQ/W9sF8FBq+g+t/Yg4Ut3kHN+WAPAAVdQE4=
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAABIAAAAIgAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9Bkp/sAcAAU8I6Q==
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_2_2_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 20 0 18 0 8">
+  <Piece Extent="0 20 0 18 0 8">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int32" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAALQAANwAAAA==eF7t1KERADAMA7Gk+w9dWpYCswiI+R66q6rDzmPa/kj3ALZK/2m6BwAAAAAAAAAAAJMLZrELTQ==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACoAAAAVQAAAA==eF5jYICAWTNBYKU9hN5pb2IMAoeh/JP2EFUXoOKX7dPTQOAaVP6m/dkzIHAHqu4BVPwhVP1jqPwTqL5nUHUvoOpeQtW9hqp7A1X3Dqrugz0ASSZF3Q==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACYAAAAVgAAAA==eF5jYIAAmeOFQLTGHkLvsQ8Fg6NQ/hn7IjDjIlT8qr1F32MgugGVv2MPMeUBVN1D+8cQBVD1T+3BymSeQ/W9sF8FBq+g+t/Yg4Ut3kHN+WAPAAVdQE4=
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAABIAAAAIgAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9Bkp/sAcAAU8I6Q==
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_2_None_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAABQAAGwAAAA==eF5jZIAAxlF6lB4AmmmUpogeDUfKaAD7jwDw
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_2_None_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAABQAAGQAAAA==eF5jZIAAxlF6lB4AmmmUHqUHkAYA/M8A8Q==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_3_1+2+3_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAZwAAAA==eF7t0rcOgmAAhVEgNmyo2AuoWN//BR04EwsJcfzvcvabL47qxcFOJg177HPAIUdMOeaEU844Z8YFl1wx55obbrnjngceeeKZFxYseeWNd1Z88MkX3/zwy+Z/wf8YOqzX1uEPlgwHCA==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_3_1+2+3_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAagAAAA==eF7t0rkOglAARFExLrgCKuKuqLj8/w9acCoSY7B+05x+cqNOvSj4l92GPfY54JAxRxxzwilnnDNhyowLLrlizjULbrjljnseeOSJZ15Y8sob76z44JMvvtn8L9jObz2GDuv96vADk5QHBg==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_3_1_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAZAAAAA==eF7t0scRglAAQEEBAyZUMCuomPtv0ANbgMNw/O+yDbyo1xQFWxkzYZ8DDjliyjEnnHLGOTMuuOSKOQuuueGWO+554JEnnlmy4oVX3ljzzgeffPHND7+Mg50aPmz698MfmvQHCg==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_3_1_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAYwAAAA==eF7t0scBgkAAAEHBgBEwgDmBsf8GfTANCN/bzzSwUa8pCrYyZp8DDjliwjEnnHLGORdMmTHnkiuuuWHBklvuuOeBR5545oVX3nhnxZoPPvnimx9+GQc7GT5sqvjvhz+ZtAcJ
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_3_2_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 20 0 18 0 8">
+  <Piece Extent="0 20 0 18 0 8">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int32" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAALQAAOgAAAA==eF7t1CEOACAMBMHS/z8aXwNJSagYMe6y8jIislgNtTW9d9oD/PL6r6b3AAAAAAAAAAAA4MYGlRYLYA==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACoAAAAVQAAAA==eF5jYICAWTNBYKU9hN5pb2IMAoeh/JP2EFUXoOKX7dPTQOAaVP6m/dkzIHAHqu4BVPwhVP1jqPwTqL5nUHUvoOpeQtW9hqp7A1X3Dqrugz0ASSZF3Q==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACYAAAAVgAAAA==eF5jYIAAmeOFQLTGHkLvsQ8Fg6NQ/hn7IjDjIlT8qr1F32MgugGVv2MPMeUBVN1D+8cQBVD1T+3BymSeQ/W9sF8FBq+g+t/Yg4Ut3kHN+WAPAAVdQE4=
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAABIAAAAIgAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9Bkp/sAcAAU8I6Q==
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_3_2_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 20 0 18 0 8">
+  <Piece Extent="0 20 0 18 0 8">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int32" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAALQAANwAAAA==eF7t1KERADAMA7Gk+w9dWpYCswiI+R66q6rDzmPa/kj3ALZK/2m6BwAAAAAAAAAAAJMLZrELTQ==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACoAAAAVQAAAA==eF5jYICAWTNBYKU9hN5pb2IMAoeh/JP2EFUXoOKX7dPTQOAaVP6m/dkzIHAHqu4BVPwhVP1jqPwTqL5nUHUvoOpeQtW9hqp7A1X3Dqrugz0ASSZF3Q==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACYAAAAVgAAAA==eF5jYIAAmeOFQLTGHkLvsQ8Fg6NQ/hn7IjDjIlT8qr1F32MgugGVv2MPMeUBVN1D+8cQBVD1T+3BymSeQ/W9sF8FBq+g+t/Yg4Ut3kHN+WAPAAVdQE4=
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAABIAAAAIgAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9Bkp/sAcAAU8I6Q==
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_3_None_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAABQAAIgAAAA==eF5jZIAAxlGaLJoJjSakntr6hzqN7v9RepSmJw0AC04A9Q==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_3_None_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAABQAALwAAAA==eF5jZIAAxlGaLJoJjSakHpc+cvUTUkdrmlL3j9KU0dROF5TqH2iaVPcDAALOANU=
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_4_1+2+3_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAcQAAAA==eF7t0rkOglAUBFAxKu6igvsKrv//gxYcm9fQGEPBNKe6yc1kolaZqPEndthljzH7HHDIEceccMoZE8654JIpM6645oZb7rjngUeeeOaFV+YseOOdDz754pthf+3Aqr7rdv9vw3+/NjssU7XDD0/8BuQ=
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_4_1+2+3_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAYQAAAA==eF7t0scVglAAAEHgqZgBA2ZExdR/gx6YCpDj38s0sEnUlgR7ccAhR0w55oRTzjjngktmzFlwxTU33LLkjnseeOSJZ15Y8cqaN975YMMnX3zzwy/j4F+GD9u6fvgD+gwHCA==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_4_1_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="41">
+        AQAAAACAAAAABQAAZAAAAA==eF7t0scRglAAQEEBAyZUMCuomPtv0ANbgMNw/O+yDbyo1xQFWxkzYZ8DDjliyjEnnHLGOTMuuOSKOQuuueGWO+554JEnnlmy4oVX3ljzzgeffPHND7+Mg50aPmz698MfmvQHCg==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_4_1_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="2" RangeMax="41">
+        AQAAAACAAAAABQAAZAAAAA==eF7t0rcSglAARFEHE0bAgBkE8///oAWnF8b2bXP6nRv1mkXBv+xzwCFHHDPmhFPOOOeCSyZMmXHFNTfcMueOex545IlnXliw5JUVa95454NPvvjmh79+DXYzdNisbYdfSqMHMg==
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_4_2_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 20 0 18 0 8">
+  <Piece Extent="0 20 0 18 0 8">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int32" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAALQAANQAAAA==eF7t1CEOACAMBMGj/380Fk+TQjJi7MqtJHVYlypv9wB+0f2+7h4AAAAAAAAAAABM2HWwC1M=
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACoAAAAVQAAAA==eF5jYICAWTNBYKU9hN5pb2IMAoeh/JP2EFUXoOKX7dPTQOAaVP6m/dkzIHAHqu4BVPwhVP1jqPwTqL5nUHUvoOpeQtW9hqp7A1X3Dqrugz0ASSZF3Q==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACYAAAAVgAAAA==eF5jYIAAmeOFQLTGHkLvsQ8Fg6NQ/hn7IjDjIlT8qr1F32MgugGVv2MPMeUBVN1D+8cQBVD1T+3BymSeQ/W9sF8FBq+g+t/Yg4Ut3kHN+WAPAAVdQE4=
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAABIAAAAIgAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9Bkp/sAcAAU8I6Q==
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_4_2_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 20 0 18 0 8">
+  <Piece Extent="0 20 0 18 0 8">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int32" Name="materialpoint" format="binary" RangeMin="1" RangeMax="1">
+        AQAAAACAAAAALQAAJAAAAA==eF7twwEJAAAMBKH7/qWXY6DgqqmqqqqqqqqqqqqqPnhyUwtB
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACoAAAAVQAAAA==eF5jYICAWTNBYKU9hN5pb2IMAoeh/JP2EFUXoOKX7dPTQOAaVP6m/dkzIHAHqu4BVPwhVP1jqPwTqL5nUHUvoOpeQtW9hqp7A1X3Dqrugz0ASSZF3Q==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAACYAAAAVgAAAA==eF5jYIAAmeOFQLTGHkLvsQ8Fg6NQ/hn7IjDjIlT8qr1F32MgugGVv2MPMeUBVN1D+8cQBVD1T+3BymSeQ/W9sF8FBq+g+t/Yg4Ut3kHN+WAPAAVdQE4=
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="1">
+        AQAAAACAAABIAAAAIgAAAA==eF5jYEAGB+wh9AUofQNKP4DST6D0Cyj9Bkp/sAcAAU8I6Q==
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_4_None_False.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAABQAAIAAAAA==eF5jZIAAxlF6lB4AmokAPdj1DzRNyP2jNH4aAMufANU=
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_4_None_True.vtr

@@ -0,0 +1,25 @@
+<?xml version="1.0"?>
+<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
+  <RectilinearGrid WholeExtent="0 8 0 5 0 4">
+  <Piece Extent="0 8 0 5 0 4">
+    <PointData>
+    </PointData>
+    <CellData>
+      <DataArray type="Int64" Name="materialpoint" format="binary" RangeMin="1" RangeMax="2">
+        AQAAAACAAAAABQAAMAAAAA==eF5jYoAAJhw0IwEalz566aeUptT+oa6fUppS+4e6fkppSu0f6voppSm1HwBAngDh
+      </DataArray>
+    </CellData>
+    <Coordinates>
+      <DataArray type="Float64" Name="x" format="binary" RangeMin="0" RangeMax="0.000008">
+        AQAAAACAAABIAAAAOgAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2D3IeNxx9YfV6DiN+22x9tFGsbchco/sAMA/fQl6g==
+      </DataArray>
+      <DataArray type="Float64" Name="y" format="binary" RangeMin="0" RangeMax="0.000005">
+        AQAAAACAAAAwAAAAKwAAAA==eF5jYICAvrdbF3w/tsEOQh+wC30iUFisdRLKv2D3MeNxx9YfV+wAD5wZgw==
+      </DataArray>
+      <DataArray type="Float64" Name="z" format="binary" RangeMin="0" RangeMax="0.000004">
+        AQAAAACAAAAoAAAAIwAAAA==eF5jYICA3rdbF3w/tsEOQh+wC3kiUFisdRLKv2AHAFVBE/w=
+      </DataArray>
+    </Coordinates>
+  </Piece>
+  </RectilinearGrid>
+</VTKFile>

python/tests/reference/Geom/clean_stencil=1.geom

@@ -1,25 +0,0 @@
-4 header
-grid   a 8 b 5 c 4
-size   x 8e-06 y 5e-06 z 4e-06
-origin x 0.0 y 0.0 z 0.0
-homogenization 1
- 1  1  2 22  2  2  1 21
- 1  1  6 26  2  2  5 25
- 1  1 10 30  2  2  9 29
- 1  1 14 34  2  2 13 33
- 1  1 18 38  2  2 17 37
- 1  1  3 23  2  2  2 22
- 1  1  7 27  2  2  6 26
- 1  1 11 31  2  2 10 30
- 1  1 15 35  2  2 14 34
- 1  1 19 39  2  2 18 38
- 1  1  4 24  2  2  3 23
- 1  1  8 28  2  2  7 27
- 1  1 12 32  2  2 11 31
- 1  1 16 36  2  2 15 35
- 1  1 20 40  2  2 19 39
- 1  1  5 25  2  2  4 24
- 1  1  9 29  2  2  8 28
- 1  1 13 33  2  2 12 32
- 1  1 17 37  2  2 16 36
- 1  1 21 41  2  2 20 40

python/tests/reference/Geom/clean_stencil=2.geom

@@ -1,25 +0,0 @@
-4 header
-grid   a 8 b 5 c 4
-size   x 8e-06 y 5e-06 z 4e-06
-origin x 0.0 y 0.0 z 0.0
-homogenization 1
- 1  1  1  2  2  2  1  1
- 1  1  1  2  2  2  2  1
- 1  1  1  6  2  2  2  5
- 1  1  1 10  2  2  2  9
- 1  1  1 14  2  2  2 13
- 1  1  1  2  2  2  2  1
- 1  1  1  2  2  2  2  1
- 1  1  1  6  2  2  2  5
- 1  1  1 10  2  2  2  9
- 1  1  1 14  2  2  2 13
- 1  1  1  3  2  2  2  2
- 1  1  1  3  2  2  2  2
- 1  1  1  7  2  2  2  6
- 1  1  1 11  2  2  2 10
- 1  1  1 15  2  2  2 14
- 1  1  1  4  2  2  2  3
- 1  1  1  4  2  2  2  3
- 1  1  1  8  2  2  2  7
- 1  1  1 12  2  2  2 11
- 1  1  1 16  2  2  2 15

python/tests/reference/Geom/clean_stencil=3.geom

@@ -1,25 +0,0 @@
-4 header
-grid   a 8 b 5 c 4
-size   x 8e-06 y 5e-06 z 4e-06
-origin x 0.0 y 0.0 z 0.0
-homogenization 1
- 1  1  1  2  2  2  2 21
- 1  1  1  2  2  2  2 21
- 1  1  1  2  2  2  2 25
- 1  1  1  2  2  2  2 29
- 1  1  1  2  2  2  2 37
- 1  1  1  2  2  2  2 21
- 1  1  1  2  2  2  2 21
- 1  1  1  2  2  2  2 25
- 1  1  1  2  2  2  2 29
- 1  1  1  2  2  2  2 37
- 1  1  1  2  2  2  2 22
- 1  1  1  2  2  2  2 22
- 1  1  1  2  2  2  2 26
- 1  1  1  2  2  2  2 30
- 1  1  1  2  2  2  2 38
- 1  1  1  2  2  2  2 24
- 1  1  1  2  2  2  2 24
- 1  1  1  2  2  2  2 28
- 1  1  1  2  2  2  2 32
- 1  1  1  2  2  2  2 40

python/tests/reference/Geom/clean_stencil=4.geom

@@ -1,25 +0,0 @@
-4 header
-grid   a 8 b 5 c 4
-size   x 8e-06 y 5e-06 z 4e-06
-origin x 0.0 y 0.0 z 0.0
-homogenization 1
-1 1 1 2 2 2 2 2
-1 1 1 2 2 2 2 2
-1 1 1 2 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 2 2 2 2 2
-1 1 1 2 2 2 2 2
-1 1 1 2 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 2 2 2 2 2
-1 1 1 2 2 2 2 2
-1 1 1 2 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 1 2 2 2 2
-1 1 1 1 2 2 2 2

python/tests/reference/Geom/flip_directions=x-y-z.geom

@@ -0,0 +1,25 @@
+4 header
+grid   a 8 b 5 c 4
+size   x 8e-06 y 5e-06 z 4e-06
+origin x 0.0 y 0.0 z 0.0
+homogenization 1
+40 39 38 37 36 35 34 33
+32 31 30 29 28 27 26 25
+24 23 22 21 20 19 18 17
+16 15 14 13 12 11 10  9
+ 8  7  6  5  4  3  2  1
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+41 40 39 38 37 36 35 34
+33 32 31 30 29 28 27 26
+25 24 23 22 21 20 19 18
+17 16 15 14 13 12 11 10
+ 9  8  7  6  5  4  3  2
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1

python/tests/reference/Geom/flip_directions=x.geom

@@ -0,0 +1,25 @@
+4 header
+grid   a 8 b 5 c 4
+size   x 8e-06 y 5e-06 z 4e-06
+origin x 0.0 y 0.0 z 0.0
+homogenization 1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 9  8  7  6  5  4  3  2
+17 16 15 14 13 12 11 10
+25 24 23 22 21 20 19 18
+33 32 31 30 29 28 27 26
+41 40 39 38 37 36 35 34
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 8  7  6  5  4  3  2  1
+16 15 14 13 12 11 10  9
+24 23 22 21 20 19 18 17
+32 31 30 29 28 27 26 25
+40 39 38 37 36 35 34 33

python/tests/reference/Geom/flip_directions=y-z.geom

@@ -0,0 +1,25 @@
+4 header
+grid   a 8 b 5 c 4
+size   x 8e-06 y 5e-06 z 4e-06
+origin x 0.0 y 0.0 z 0.0
+homogenization 1
+33 34 35 36 37 38 39 40
+25 26 27 28 29 30 31 32
+17 18 19 20 21 22 23 24
+ 9 10 11 12 13 14 15 16
+ 1  2  3  4  5  6  7  8
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+34 35 36 37 38 39 40 41
+26 27 28 29 30 31 32 33
+18 19 20 21 22 23 24 25
+10 11 12 13 14 15 16 17
+ 2  3  4  5  6  7  8  9
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1

python/tests/reference/Geom/flip_directions=z-x-y.geom

@@ -0,0 +1,25 @@
+4 header
+grid   a 8 b 5 c 4
+size   x 8e-06 y 5e-06 z 4e-06
+origin x 0.0 y 0.0 z 0.0
+homogenization 1
+40 39 38 37 36 35 34 33
+32 31 30 29 28 27 26 25
+24 23 22 21 20 19 18 17
+16 15 14 13 12 11 10  9
+ 8  7  6  5  4  3  2  1
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+ 2  2  2  2  2  2  2  2
+41 40 39 38 37 36 35 34
+33 32 31 30 29 28 27 26
+25 24 23 22 21 20 19 18
+17 16 15 14 13 12 11 10
+ 9  8  7  6  5  4  3  2
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1
+ 1  1  1  1  1  1  1  1

python/tests/reference/Geom/mirror_directions=x-y-z_reflect=True.geom

@@ -3,83 +3,83 @@ grid   a 16 b 10 c 8
 size   x 1.6e-05 y 1e-05 z 8e-06
 origin x 0.0 y 0.0 z 0.0
 homogenization 1
- 1  1  2 22  2  2  1 21 21  1  2  2 22  2  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1  6 26  2  2  5 25 25  5  2  2 26  6  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 10 30  2  2  9 29 29  9  2  2 30 10  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 14 34  2  2 13 33 33 13  2  2 34 14  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 18 38  2  2 17 37 37 17  2  2 38 18  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 18 38  2  2 17 37 37 17  2  2 38 18  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 14 34  2  2 13 33 33 13  2  2 34 14  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 10 30  2  2  9 29 29  9  2  2 30 10  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1  6 26  2  2  5 25 25  5  2  2 26  6  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1  2 22  2  2  1 21 21  1  2  2 22  2  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1  3 23  2  2  2 22 22  2  2  2 23  3  1  1
+ 2  3  4  5  6  7  8  9  9  8  7  6  5  4  3  2
- 1  1  7 27  2  2  6 26 26  6  2  2 27  7  1  1
+10 11 12 13 14 15 16 17 17 16 15 14 13 12 11 10
- 1  1 11 31  2  2 10 30 30 10  2  2 31 11  1  1
+18 19 20 21 22 23 24 25 25 24 23 22 21 20 19 18
- 1  1 15 35  2  2 14 34 34 14  2  2 35 15  1  1
+26 27 28 29 30 31 32 33 33 32 31 30 29 28 27 26
- 1  1 19 39  2  2 18 38 38 18  2  2 39 19  1  1
+34 35 36 37 38 39 40 41 41 40 39 38 37 36 35 34
- 1  1 19 39  2  2 18 38 38 18  2  2 39 19  1  1
+34 35 36 37 38 39 40 41 41 40 39 38 37 36 35 34
- 1  1 15 35  2  2 14 34 34 14  2  2 35 15  1  1
+26 27 28 29 30 31 32 33 33 32 31 30 29 28 27 26
- 1  1 11 31  2  2 10 30 30 10  2  2 31 11  1  1
+18 19 20 21 22 23 24 25 25 24 23 22 21 20 19 18
- 1  1  7 27  2  2  6 26 26  6  2  2 27  7  1  1
+10 11 12 13 14 15 16 17 17 16 15 14 13 12 11 10
- 1  1  3 23  2  2  2 22 22  2  2  2 23  3  1  1
+ 2  3  4  5  6  7  8  9  9  8  7  6  5  4  3  2
- 1  1  4 24  2  2  3 23 23  3  2  2 24  4  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  8 28  2  2  7 27 27  7  2  2 28  8  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 12 32  2  2 11 31 31 11  2  2 32 12  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 16 36  2  2 15 35 35 15  2  2 36 16  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 20 40  2  2 19 39 39 19  2  2 40 20  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 20 40  2  2 19 39 39 19  2  2 40 20  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 16 36  2  2 15 35 35 15  2  2 36 16  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 12 32  2  2 11 31 31 11  2  2 32 12  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  8 28  2  2  7 27 27  7  2  2 28  8  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  4 24  2  2  3 23 23  3  2  2 24  4  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  5 25  2  2  4 24 24  4  2  2 25  5  1  1
+ 1  2  3  4  5  6  7  8  8  7  6  5  4  3  2  1
- 1  1  9 29  2  2  8 28 28  8  2  2 29  9  1  1
+ 9 10 11 12 13 14 15 16 16 15 14 13 12 11 10  9
- 1  1 13 33  2  2 12 32 32 12  2  2 33 13  1  1
+17 18 19 20 21 22 23 24 24 23 22 21 20 19 18 17
- 1  1 17 37  2  2 16 36 36 16  2  2 37 17  1  1
+25 26 27 28 29 30 31 32 32 31 30 29 28 27 26 25
- 1  1 21 41  2  2 20 40 40 20  2  2 41 21  1  1
+33 34 35 36 37 38 39 40 40 39 38 37 36 35 34 33
- 1  1 21 41  2  2 20 40 40 20  2  2 41 21  1  1
+33 34 35 36 37 38 39 40 40 39 38 37 36 35 34 33
- 1  1 17 37  2  2 16 36 36 16  2  2 37 17  1  1
+25 26 27 28 29 30 31 32 32 31 30 29 28 27 26 25
- 1  1 13 33  2  2 12 32 32 12  2  2 33 13  1  1
+17 18 19 20 21 22 23 24 24 23 22 21 20 19 18 17
- 1  1  9 29  2  2  8 28 28  8  2  2 29  9  1  1
+ 9 10 11 12 13 14 15 16 16 15 14 13 12 11 10  9
- 1  1  5 25  2  2  4 24 24  4  2  2 25  5  1  1
+ 1  2  3  4  5  6  7  8  8  7  6  5  4  3  2  1
- 1  1  5 25  2  2  4 24 24  4  2  2 25  5  1  1
+ 1  2  3  4  5  6  7  8  8  7  6  5  4  3  2  1
- 1  1  9 29  2  2  8 28 28  8  2  2 29  9  1  1
+ 9 10 11 12 13 14 15 16 16 15 14 13 12 11 10  9
- 1  1 13 33  2  2 12 32 32 12  2  2 33 13  1  1
+17 18 19 20 21 22 23 24 24 23 22 21 20 19 18 17
- 1  1 17 37  2  2 16 36 36 16  2  2 37 17  1  1
+25 26 27 28 29 30 31 32 32 31 30 29 28 27 26 25
- 1  1 21 41  2  2 20 40 40 20  2  2 41 21  1  1
+33 34 35 36 37 38 39 40 40 39 38 37 36 35 34 33
- 1  1 21 41  2  2 20 40 40 20  2  2 41 21  1  1
+33 34 35 36 37 38 39 40 40 39 38 37 36 35 34 33
- 1  1 17 37  2  2 16 36 36 16  2  2 37 17  1  1
+25 26 27 28 29 30 31 32 32 31 30 29 28 27 26 25
- 1  1 13 33  2  2 12 32 32 12  2  2 33 13  1  1
+17 18 19 20 21 22 23 24 24 23 22 21 20 19 18 17
- 1  1  9 29  2  2  8 28 28  8  2  2 29  9  1  1
+ 9 10 11 12 13 14 15 16 16 15 14 13 12 11 10  9
- 1  1  5 25  2  2  4 24 24  4  2  2 25  5  1  1
+ 1  2  3  4  5  6  7  8  8  7  6  5  4  3  2  1
- 1  1  4 24  2  2  3 23 23  3  2  2 24  4  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  8 28  2  2  7 27 27  7  2  2 28  8  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 12 32  2  2 11 31 31 11  2  2 32 12  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 16 36  2  2 15 35 35 15  2  2 36 16  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 20 40  2  2 19 39 39 19  2  2 40 20  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 20 40  2  2 19 39 39 19  2  2 40 20  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 16 36  2  2 15 35 35 15  2  2 36 16  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 12 32  2  2 11 31 31 11  2  2 32 12  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  8 28  2  2  7 27 27  7  2  2 28  8  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  4 24  2  2  3 23 23  3  2  2 24  4  1  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  3 23  2  2  2 22 22  2  2  2 23  3  1  1
+ 2  3  4  5  6  7  8  9  9  8  7  6  5  4  3  2
- 1  1  7 27  2  2  6 26 26  6  2  2 27  7  1  1
+10 11 12 13 14 15 16 17 17 16 15 14 13 12 11 10
- 1  1 11 31  2  2 10 30 30 10  2  2 31 11  1  1
+18 19 20 21 22 23 24 25 25 24 23 22 21 20 19 18
- 1  1 15 35  2  2 14 34 34 14  2  2 35 15  1  1
+26 27 28 29 30 31 32 33 33 32 31 30 29 28 27 26
- 1  1 19 39  2  2 18 38 38 18  2  2 39 19  1  1
+34 35 36 37 38 39 40 41 41 40 39 38 37 36 35 34
- 1  1 19 39  2  2 18 38 38 18  2  2 39 19  1  1
+34 35 36 37 38 39 40 41 41 40 39 38 37 36 35 34
- 1  1 15 35  2  2 14 34 34 14  2  2 35 15  1  1
+26 27 28 29 30 31 32 33 33 32 31 30 29 28 27 26
- 1  1 11 31  2  2 10 30 30 10  2  2 31 11  1  1
+18 19 20 21 22 23 24 25 25 24 23 22 21 20 19 18
- 1  1  7 27  2  2  6 26 26  6  2  2 27  7  1  1
+10 11 12 13 14 15 16 17 17 16 15 14 13 12 11 10
- 1  1  3 23  2  2  2 22 22  2  2  2 23  3  1  1
+ 2  3  4  5  6  7  8  9  9  8  7  6  5  4  3  2
- 1  1  2 22  2  2  1 21 21  1  2  2 22  2  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1  6 26  2  2  5 25 25  5  2  2 26  6  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 10 30  2  2  9 29 29  9  2  2 30 10  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 14 34  2  2 13 33 33 13  2  2 34 14  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 18 38  2  2 17 37 37 17  2  2 38 18  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 18 38  2  2 17 37 37 17  2  2 38 18  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 14 34  2  2 13 33 33 13  2  2 34 14  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 10 30  2  2  9 29 29  9  2  2 30 10  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1  6 26  2  2  5 25 25  5  2  2 26  6  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1  2 22  2  2  1 21 21  1  2  2 22  2  1  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1

python/tests/reference/Geom/mirror_directions=x_reflect=False.geom

@@ -3,23 +3,23 @@ grid   a 14 b 5 c 4
 size   x 1.4e-05 y 5e-06 z 4e-06
 origin x 0.0 y 0.0 z 0.0
 homogenization 1
- 1  1  2 22  2  2  1 21  1  2  2 22  2  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1  6 26  2  2  5 25  5  2  2 26  6  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 10 30  2  2  9 29  9  2  2 30 10  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 14 34  2  2 13 33 13  2  2 34 14  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1 18 38  2  2 17 37 17  2  2 38 18  1
+ 1  1  1  1  1  1  1  1  1  1  1  1  1  1
- 1  1  3 23  2  2  2 22  2  2  2 23  3  1
+ 2  3  4  5  6  7  8  9  8  7  6  5  4  3
- 1  1  7 27  2  2  6 26  6  2  2 27  7  1
+10 11 12 13 14 15 16 17 16 15 14 13 12 11
- 1  1 11 31  2  2 10 30 10  2  2 31 11  1
+18 19 20 21 22 23 24 25 24 23 22 21 20 19
- 1  1 15 35  2  2 14 34 14  2  2 35 15  1
+26 27 28 29 30 31 32 33 32 31 30 29 28 27
- 1  1 19 39  2  2 18 38 18  2  2 39 19  1
+34 35 36 37 38 39 40 41 40 39 38 37 36 35
- 1  1  4 24  2  2  3 23  3  2  2 24  4  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  8 28  2  2  7 27  7  2  2 28  8  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 12 32  2  2 11 31 11  2  2 32 12  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 16 36  2  2 15 35 15  2  2 36 16  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1 20 40  2  2 19 39 19  2  2 40 20  1
+ 2  2  2  2  2  2  2  2  2  2  2  2  2  2
- 1  1  5 25  2  2  4 24  4  2  2 25  5  1
+ 1  2  3  4  5  6  7  8  7  6  5  4  3  2
- 1  1  9 29  2  2  8 28  8  2  2 29  9  1
+ 9 10 11 12 13 14 15 16 15 14 13 12 11 10
- 1  1 13 33  2  2 12 32 12  2  2 33 13  1
+17 18 19 20 21 22 23 24 23 22 21 20 19 18
- 1  1 17 37  2  2 16 36 16  2  2 37 17  1
+25 26 27 28 29 30 31 32 31 30 29 28 27 26
- 1  1 21 41  2  2 20 40 20  2  2 41 21  1
+33 34 35 36 37 38 39 40 39 38 37 36 35 34

python/tests/reference/Geom/mirror_directions=y-z_reflect=False.geom

@@ -3,51 +3,51 @@ grid   a 8 b 8 c 6
 size   x 8e-06 y 8.000000000000001e-06 z 6e-06
 origin x 0.0 y 0.0 z 0.0
 homogenization 1
- 1  1  2 22  2  2  1 21
+ 1  1  1  1  1  1  1  1
- 1  1  6 26  2  2  5 25
+ 1  1  1  1  1  1  1  1
- 1  1 10 30  2  2  9 29
+ 1  1  1  1  1  1  1  1
- 1  1 14 34  2  2 13 33
+ 1  1  1  1  1  1  1  1
- 1  1 18 38  2  2 17 37
+ 1  1  1  1  1  1  1  1
- 1  1 14 34  2  2 13 33
+ 1  1  1  1  1  1  1  1
- 1  1 10 30  2  2  9 29
+ 1  1  1  1  1  1  1  1
- 1  1  6 26  2  2  5 25
+ 1  1  1  1  1  1  1  1
- 1  1  3 23  2  2  2 22
+ 2  3  4  5  6  7  8  9
- 1  1  7 27  2  2  6 26
+10 11 12 13 14 15 16 17
- 1  1 11 31  2  2 10 30
+18 19 20 21 22 23 24 25
- 1  1 15 35  2  2 14 34
+26 27 28 29 30 31 32 33
- 1  1 19 39  2  2 18 38
+34 35 36 37 38 39 40 41
- 1  1 15 35  2  2 14 34
+26 27 28 29 30 31 32 33
- 1  1 11 31  2  2 10 30
+18 19 20 21 22 23 24 25
- 1  1  7 27  2  2  6 26
+10 11 12 13 14 15 16 17
- 1  1  4 24  2  2  3 23
+ 2  2  2  2  2  2  2  2
- 1  1  8 28  2  2  7 27
+ 2  2  2  2  2  2  2  2
- 1  1 12 32  2  2 11 31
+ 2  2  2  2  2  2  2  2
- 1  1 16 36  2  2 15 35
+ 2  2  2  2  2  2  2  2
- 1  1 20 40  2  2 19 39
+ 2  2  2  2  2  2  2  2
- 1  1 16 36  2  2 15 35
+ 2  2  2  2  2  2  2  2
- 1  1 12 32  2  2 11 31
+ 2  2  2  2  2  2  2  2
- 1  1  8 28  2  2  7 27
+ 2  2  2  2  2  2  2  2
- 1  1  5 25  2  2  4 24
+ 1  2  3  4  5  6  7  8
- 1  1  9 29  2  2  8 28
+ 9 10 11 12 13 14 15 16
- 1  1 13 33  2  2 12 32
+17 18 19 20 21 22 23 24
- 1  1 17 37  2  2 16 36
+25 26 27 28 29 30 31 32
- 1  1 21 41  2  2 20 40
+33 34 35 36 37 38 39 40
- 1  1 17 37  2  2 16 36
+25 26 27 28 29 30 31 32
- 1  1 13 33  2  2 12 32
+17 18 19 20 21 22 23 24
- 1  1  9 29  2  2  8 28
+ 9 10 11 12 13 14 15 16
- 1  1  4 24  2  2  3 23
+ 2  2  2  2  2  2  2  2
- 1  1  8 28  2  2  7 27
+ 2  2  2  2  2  2  2  2
- 1  1 12 32  2  2 11 31
+ 2  2  2  2  2  2  2  2
- 1  1 16 36  2  2 15 35
+ 2  2  2  2  2  2  2  2
- 1  1 20 40  2  2 19 39
+ 2  2  2  2  2  2  2  2
- 1  1 16 36  2  2 15 35
+ 2  2  2  2  2  2  2  2
- 1  1 12 32  2  2 11 31
+ 2  2  2  2  2  2  2  2
- 1  1  8 28  2  2  7 27
+ 2  2  2  2  2  2  2  2
- 1  1  3 23  2  2  2 22
+ 2  3  4  5  6  7  8  9
- 1  1  7 27  2  2  6 26
+10 11 12 13 14 15 16 17
- 1  1 11 31  2  2 10 30
+18 19 20 21 22 23 24 25
- 1  1 15 35  2  2 14 34
+26 27 28 29 30 31 32 33
- 1  1 19 39  2  2 18 38
+34 35 36 37 38 39 40 41
- 1  1 15 35  2  2 14 34
+26 27 28 29 30 31 32 33
- 1  1 11 31  2  2 10 30
+18 19 20 21 22 23 24 25
- 1  1  7 27  2  2  6 26
+10 11 12 13 14 15 16 17