From 34a08bd8b8939d786d0ce1fda1015e0eff0e2243 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 22 Aug 2020 17:04:21 +0200
Subject: [PATCH 01/41] pip package is pillow, not PIL

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

diff --git a/python/setup.py b/python/setup.py
index 360918b38..c13ddbae0 100644
--- a/python/setup.py
+++ b/python/setup.py
@@ -21,7 +21,7 @@ setuptools.setup(
         "h5py",
         "vtk",
         "matplotlib",
-        "PIL",
+        "pillow",
     ],
     classifiers = [
         "Intended Audience :: Science/Research",

From 9ca0e409661661d09a97413b0cd17c2ef956d97a Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 22 Aug 2020 20:22:34 +0200
Subject: [PATCH 02/41] simplified

---
 README | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README b/README
index 7fc372881..81ea47599 100644
--- a/README
+++ b/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

From bd4fb562bcbc6b983abd17f89d25d6e07559458c Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 22 Aug 2020 20:57:42 +0200
Subject: [PATCH 03/41] make use of inhereted dependencies

---
 python/setup.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/python/setup.py b/python/setup.py
index c13ddbae0..ac2363160 100644
--- a/python/setup.py
+++ b/python/setup.py
@@ -16,12 +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",
-        "pillow",
+        "matplotlib",                                                                               # requires numpy, pillow
     ],
     classifiers = [
         "Intended Audience :: Science/Research",

From e007aed407af44a3af2ca617540df7080c68f652 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 23 Aug 2020 13:28:15 +0200
Subject: [PATCH 04/41] updates in PRIVATE

- lightweigth tarball: no examples, no Makefile (explicit cmake)
- improved build instructions for ubuntu/debian package
---
 PRIVATE | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/PRIVATE b/PRIVATE
index a16d1e45a..8d7f2b665 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit a16d1e45a2ed925e12244b0879b9d7e5a58d973b
+Subproject commit 8d7f2b665cfee7fbafb387201558bd27c54b2abb

From 7b50a3b36436acb9e442e8bdb6af58f51f7d5af2 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 23 Aug 2020 15:20:11 +0200
Subject: [PATCH 05/41] material.config is not used anymore

---
 python/.coveragerc               |   1 -
 python/damask/__init__.py        |   1 -
 python/damask/config/__init__.py |   3 -
 python/damask/config/material.py | 282 -------------------------------
 4 files changed, 287 deletions(-)
 delete mode 100644 python/damask/config/__init__.py
 delete mode 100644 python/damask/config/material.py

diff --git a/python/.coveragerc b/python/.coveragerc
index 5daa25bb2..97114fb82 100644
--- a/python/.coveragerc
+++ b/python/.coveragerc
@@ -2,4 +2,3 @@
 omit = tests/*
        damask/_asciitable.py
        damask/_test.py
-       damask/config/*
diff --git a/python/damask/__init__.py b/python/damask/__init__.py
index b0feee4ac..1404e88d1 100644
--- a/python/damask/__init__.py
+++ b/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
diff --git a/python/damask/config/__init__.py b/python/damask/config/__init__.py
deleted file mode 100644
index f5c57a879..000000000
--- a/python/damask/config/__init__.py
+++ /dev/null
@@ -1,3 +0,0 @@
-"""Aggregator for configuration file handling."""
-
-from .material import Material    # noqa
diff --git a/python/damask/config/material.py b/python/damask/config/material.py
deleted file mode 100644
index 433c70d03..000000000
--- a/python/damask/config/material.py
+++ /dev/null
@@ -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()]

From a814db5b643c016c5936e5497c6180913591bd3c Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 23 Aug 2020 19:31:30 +0200
Subject: [PATCH 06/41] PRIVATE: Master includes yaml changes

---
 PRIVATE | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/PRIVATE b/PRIVATE
index 8d7f2b665..4715b68a1 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 8d7f2b665cfee7fbafb387201558bd27c54b2abb
+Subproject commit 4715b68a1b55d8f18fe6b405e312785489d9c9d0

From 446ac03b072117ec8979973027df0667206084ed Mon Sep 17 00:00:00 2001
From: Philip Eisenlohr <eisenlohr@egr.msu.edu>
Date: Sun, 23 Aug 2020 18:46:01 -0400
Subject: [PATCH 07/41] All geom methods are now out-of-place, i.e. return an
 updated duplicate (to allow for daisy chaining).

* Added comments when methods acted.
* Added diff method
* Added flip method
* Fixed add_primitive inversion bug (again...)
* Fixed cell centering bug in add_primitive

* Added missing tests
---
 python/damask/_geom.py    | 255 ++++++++++++++++++++++++++++----------
 python/tests/test_Geom.py | 124 ++++++++++--------
 2 files changed, 263 insertions(+), 116 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 0f1c35e42..c4bbbbdef 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -7,6 +7,7 @@ from functools import partial
 import numpy as np
 from scipy import ndimage,spatial
 
+from . import version
 from . import environment
 from . import Rotation
 from . import VTK
@@ -24,15 +25,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 : int, optional
-            homogenization index.
+            Homogenization index.
         comments : list of str, optional
-            comments lines.
+            Comment lines.
 
         """
         self.set_microstructure(microstructure)
@@ -50,7 +51,7 @@ class Geom:
                f'origin   x y z:      {util.srepr(self.get_origin(),"   ")}',
                f'# microstructures:   {self.N_microstructure}',
                f'max microstructure:  {np.nanmax(self.microstructure)}',
-              ])
+              ]+self.get_comments())
 
 
     def __copy__(self):
@@ -63,20 +64,57 @@ class Geom:
         return self.__copy__()
 
 
-    def update(self,microstructure=None,size=None,origin=None,rescale=False):
+    def duplicate(self,microstructure=None,size=None,origin=None,comments=None,autosize=False):
+        """
+        Create a duplicate having updated microstructure, size, and origin.
+
+        Parameters
+        ----------
+        microstructure : numpy.ndarray, optional
+            Microstructure array (3D).
+        size : list or numpy.ndarray, optional
+            Physical size of the microstructure in meter.
+        origin : list or numpy.ndarray, optional
+            Physical origin of the microstructure in meter.
+        comments : list of str, optional
+            Comment lines.
+        autosize : bool, optional
+            Ignore size parameter and rescale according to change of grid points.
+
+        """
+        if size is not None and autosize:
+            raise ValueError('Auto sizing conflicts with explicit size parameter.')
+
+        grid_old = self.get_grid()
+        dup = self.copy()
+        dup.set_microstructure(microstructure)
+        dup.set_origin(origin)
+
+        if comments is not None:
+            dup.set_comments(comments)
+
+        if size is not None:
+            dup.set_size(size)
+        elif autosize:
+            dup.set_size(dup.get_grid()/grid_old*self.get_size())
+
+        return dup
+
+
+    def update(self,microstructure=None,size=None,origin=None,autosize=False):
         """
         Update microstructure and size.
 
         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.
-        rescale : bool, optional
-            ignore size parameter and rescale according to change of grid points.
+            Physical origin of the microstructure in meter.
+        autosize : bool, optional
+            Ignore size parameter and rescale according to change of grid points.
 
         """
         grid_old    = self.get_grid()
@@ -85,15 +123,15 @@ class Geom:
         unique_old  = self.N_microstructure
         max_old     = np.nanmax(self.microstructure)
 
-        if size is not None and rescale:
-            raise ValueError('Either set size explicitly or rescale automatically')
+        if size is not None and autosize:
+            raise ValueError('Auto sizing conflicts with explicit size parameter.')
 
         self.set_microstructure(microstructure)
         self.set_origin(origin)
 
         if size is not None:
             self.set_size(size)
-        elif rescale:
+        elif autosize:
             self.set_size(self.get_grid()/grid_old*self.size)
 
         message = [f'grid     a b c:      {util.srepr(grid_old," x ")}']
@@ -124,6 +162,40 @@ class Geom:
         return util.return_message(message)
 
 
+    def diff(self,other):
+        """
+        Report property differences of self relative to other.
+
+        Parameters
+        ----------
+        other : Geom
+            Geometry to compare self against.
+
+        """
+        message = []
+        if np.any(other.get_grid() != self.get_grid()):
+            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 ")}'))
+
+        if np.any(other.get_size() != self.get_size()):
+            message.append(util.delete(f'size     x y z:      {util.srepr(other.get_size()," x ")}'))
+            message.append(util.emph(  f'size     x y z:      {util.srepr( self.get_size()," x ")}'))
+
+        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'# microstructures:   {other.N_microstructure}'))
+            message.append(util.emph(  f'# microstructures:   { self.N_microstructure}'))
+
+        if np.nanmax(other.microstructure) != np.nanmax(self.microstructure):
+            message.append(util.delete(f'max microstructure:  {np.nanmax(other.microstructure)}'))
+            message.append(util.emph(  f'max microstructure:  {np.nanmax( self.microstructure)}'))
+
+        return util.return_message(message)
+
+
     def set_comments(self,comments):
         """
         Replace all existing comments.
@@ -131,7 +203,7 @@ class Geom:
         Parameters
         ----------
         comments : list of str
-            new comments.
+            All comments.
 
         """
         self.comments = []
@@ -145,7 +217,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)]
@@ -172,6 +244,10 @@ class Geom:
             else:
                 self.microstructure = np.copy(microstructure)
 
+        if         self.microstructure.dtype == 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']:
@@ -185,7 +261,7 @@ 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:
@@ -205,7 +281,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:
@@ -222,12 +298,12 @@ class Geom:
         Parameters
         ----------
         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
 
@@ -419,8 +495,11 @@ 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=f'geom.py:from_Laguerre_tessellation v{version}',
+                   )
 
 
     @staticmethod
@@ -444,8 +523,11 @@ class Geom:
         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=f'geom.py:from_Voronoi_tessellation v{version}',
+                   )
 
 
     def to_file(self,fname,pack=None):
@@ -535,7 +617,7 @@ class Geom:
 
     def show(self):
         """Show raw content (as in file)."""
-        f=StringIO()
+        f = StringIO()
         self.to_file(f)
         f.seek(0)
         return ''.join(f.readlines())
@@ -565,10 +647,10 @@ class Geom:
         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.
+            Retain original microstructure within primitive and fill outside.
+            Defaults to False.
         periodic : Boolean, optional
-            Repeat primitive over boundaries. Defaults to False.
+            Repeat primitive over boundaries. Defaults to True.
 
         """
         # normalized 'radius' and center
@@ -578,11 +660,11 @@ class Geom:
             (np.array(center) - self.origin)/self.size
 
         coords = grid_filters.cell_coord0(self.grid,np.ones(3)) \
-               - (np.ones(3)*0.5 if periodic else c)                                                # center if periodic
+               - ((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
         coords_rot = R.broadcast_to(tuple(self.grid))@coords
 
         with np.errstate(over='ignore',under='ignore'):
-            mask = np.where(np.sum(np.abs(coords_rot/r)**(2.0**exponent),axis=-1) < 1,True,False)
+            mask = np.where(np.sum(np.abs(coords_rot/r)**(2.0**exponent),axis=-1) <= 1.0,False,True)
 
         if periodic:                                                                                # translate back to center
             mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))
@@ -590,7 +672,9 @@ class Geom:
         fill_ = np.full_like(self.microstructure,np.nanmax(self.microstructure)+1 if fill is None else fill)
         ms = np.ma.MaskedArray(fill_,np.logical_not(mask) if inverse else mask)
 
-        return self.update(ms)
+        return self.duplicate(ms,
+                              comments=self.get_comments()+[f'geom.py:add_primitive v{version}'],
+                             )
 
 
     def mirror(self,directions,reflect=False):
@@ -622,8 +706,41 @@ 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.update(ms,rescale=True)
+        return self.duplicate(ms,
+                              comments=self.get_comments()+[f'geom.py:mirror {set(directions)} v{version}'],
+                              autosize=True)
+
+
+    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 all(isinstance(d, str) for d in directions):
+            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]
+        ms = self.get_microstructure()
+
+        if 'z' in directions:
+            ms = ms[:,:,limits[0]:limits[1]:-1]
+        if 'y' in directions:
+            ms = ms[:,limits[0]:limits[1]:-1,:]
+        if 'x' in directions:
+            ms = ms[limits[0]:limits[1]:-1,:,:]
+
+        return self.duplicate(ms,
+                              comments=self.get_comments()+[f'geom.py:flip {set(directions)} v{version}'],
+                             )
 
 
     def scale(self,grid):
@@ -636,17 +753,16 @@ class Geom:
             Number of grid points in x,y,z direction.
 
         """
-        #ToDo: self.add_comments('geom.py:scale v{}'.format(version)
-        return self.update(
-                           ndimage.interpolation.zoom(
-                                                      self.microstructure,
-                                                      grid/self.get_grid(),
-                                                      output=self.microstructure.dtype,
-                                                      order=0,
-                                                      mode='nearest',
-                                                      prefilter=False
-                                                     )
-                          )
+        return self.duplicate(ndimage.interpolation.zoom(
+                                                         self.microstructure,
+                                                         grid/self.get_grid(),
+                                                         output=self.microstructure.dtype,
+                                                         order=0,
+                                                         mode='nearest',
+                                                         prefilter=False
+                                                        ),
+                              comments=self.get_comments()+[f'geom.py:scale {grid} v{version}'],
+                             )
 
 
     def clean(self,stencil=3,mode='nearest',selection=None):
@@ -672,15 +788,15 @@ class Geom:
             else:
                 return me
 
-        #ToDo: self.add_comments('geom.py:clean v{}'.format(version)
-        return self.update(ndimage.filters.generic_filter(
-                                                          self.microstructure,
-                                                          mostFrequent,
-                                                          size=(stencil if selection is None else stencil//2*2+1,)*3,
-                                                          mode=mode,
-                                                          extra_keywords=dict(selection=selection),
-                                                         ).astype(self.microstructure.dtype)
-                          )
+        return self.duplicate(ndimage.filters.generic_filter(
+                                                             self.microstructure,
+                                                             mostFrequent,
+                                                             size=(stencil if selection is None else stencil//2*2+1,)*3,
+                                                             mode=mode,
+                                                             extra_keywords=dict(selection=selection),
+                                                            ).astype(self.microstructure.dtype),
+                              comments=self.get_comments()+[f'geom.py:clean {stencil} {mode} v{version}'],
+                             )
 
 
     def renumber(self):
@@ -689,8 +805,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.update(renumbered)
+        return self.duplicate(renumbered,
+                              comments=self.get_comments()+[f'geom.py:renumber v{version}'],
+                             )
 
 
     def rotate(self,R,fill=None):
@@ -724,8 +841,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.update(microstructure_in,origin=origin,rescale=True)
+        return self.duplicate(microstructure_in,
+                              origin=origin,
+                              comments=self.get_comments()+[f'geom.py:rotate {R.as_quaternion()} v{version}'],
+                              autosize=True,
+                             )
 
 
     def canvas(self,grid=None,offset=None,fill=None):
@@ -757,8 +877,11 @@ class Geom:
 
         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.update(canvas,origin=self.origin+offset*self.size/self.grid,rescale=True)
+        return self.duplicate(canvas,
+                              origin=self.origin+offset*self.size/self.grid,
+                              comments=self.get_comments()+[f'geom.py:canvas {grid} {offset} {fill} v{version}'],
+                              autosize=True,
+                             )
 
 
     def substitute(self,from_microstructure,to_microstructure):
@@ -777,8 +900,9 @@ class Geom:
         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.update(substituted)
+        return self.duplicate(substituted,
+                              comments=self.get_comments()+[f'geom.py:substitute v{version}'],
+                             )
 
 
     def vicinity_offset(self,vicinity=1,offset=None,trigger=[],periodic=True):
@@ -819,9 +943,10 @@ class Geom:
         mask = ndimage.filters.generic_filter(self.microstructure,
                                               tainted_neighborhood,
                                               size=1+2*vicinity,
-                                              mode=('wrap' if periodic else 'nearest'),
+                                              mode='wrap' if periodic else 'nearest',
                                               extra_keywords={'trigger':trigger})
         microstructure = np.ma.MaskedArray(self.microstructure + offset_, np.logical_not(mask))
 
-        #ToDo: self.add_comments('geom.py:vicinity_offset v{}'.format(version)
-        return self.update(microstructure)
+        return self.duplicate(microstructure,
+                              comments=self.get_comments()+[f'geom.py:vicinity_offset {vicinity} v{version}'],
+                             )
diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 4fb1ae00a..e568b4a53 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -28,6 +28,19 @@ def reference_dir(reference_dir_base):
 
 class TestGeom:
 
+    @pytest.mark.parametrize('flavor',['plain','explicit'])
+    def test_duplicate(self,default,flavor):
+        if flavor == 'plain':
+            modified = default.duplicate()
+        elif flavor == 'explicit':
+            modified = default.duplicate(
+                            default.get_microstructure(),
+                            default.get_size(),
+                            default.get_origin()
+                           )
+        print(modified)
+        assert geom_equal(default,modified)
+
     def test_update(self,default):
         modified = default.copy()
         modified.update(
@@ -36,7 +49,7 @@ class TestGeom:
                         default.get_origin()
                        )
         print(modified)
-        assert geom_equal(modified,default)
+        assert geom_equal(default,modified)
 
     @pytest.mark.parametrize('masked',[True,False])
     def test_set_microstructure(self,default,masked):
@@ -47,36 +60,38 @@ class TestGeom:
 
 
     def test_write_read_str(self,default,tmpdir):
-        default.to_file(str(tmpdir.join('default.geom')))
-        new = Geom.from_file(str(tmpdir.join('default.geom')))
-        assert geom_equal(new,default)
+        default.to_file(tmpdir/'default.geom')
+        new = Geom.from_file(tmpdir/'default.geom')
+        assert geom_equal(default,new)
 
     def test_write_read_file(self,default,tmpdir):
-        with open(tmpdir.join('default.geom'),'w') as f:
+        with open(tmpdir/'default.geom','w') as f:
             default.to_file(f)
-        with open(tmpdir.join('default.geom')) as f:
+        with open(tmpdir/'default.geom') as f:
             new = Geom.from_file(f)
-        assert geom_equal(new,default)
+        assert geom_equal(default,new)
 
     def test_write_show(self,default,tmpdir):
-        with open(tmpdir.join('str.geom'),'w') as f:
+        with open(tmpdir/'str.geom','w') as f:
             f.write(default.show())
-        with open(tmpdir.join('str.geom')) as f:
+        with open(tmpdir/'str.geom') as f:
             new = Geom.from_file(f)
-        assert geom_equal(new,default)
+        assert geom_equal(default,new)
+
+    def test_export_import_vtk(self,default,tmpdir):
+        default.to_vtk(tmpdir/'default')
+        assert geom_equal(default,Geom.from_vtk(tmpdir/'default.vtr'))
 
-    def test_export_vtk(self,default,tmpdir):
-        default.to_vtk(str(tmpdir.join('default')))
 
     @pytest.mark.parametrize('pack',[True,False])
     def test_pack(self,default,tmpdir,pack):
-        default.to_file(tmpdir.join('default.geom'),pack=pack)
-        new = Geom.from_file(tmpdir.join('default.geom'))
+        default.to_file(tmpdir/'default.geom',pack=pack)
+        new = Geom.from_file(tmpdir/'default.geom')
         assert geom_equal(new,default)
 
     def test_invalid_combination(self,default):
         with pytest.raises(ValueError):
-            default.update(default.microstructure[1:,1:,1:],size=np.ones(3), rescale=True)
+            default.update(default.microstructure[1:,1:,1:],size=np.ones(3), autosize=True)
 
     def test_invalid_size(self,default):
         with pytest.raises(ValueError):
@@ -108,21 +123,36 @@ class TestGeom:
                                                   ]
                             )
     def test_mirror(self,default,update,reference_dir,directions,reflect):
-        modified = default.copy()
-        modified.mirror(directions,reflect)
+        modified = default.mirror(directions,reflect)
         tag = f'directions={"-".join(directions)}_reflect={reflect}'
         reference = reference_dir/f'mirror_{tag}.geom'
         if update: modified.to_file(reference)
-        assert geom_equal(modified,Geom.from_file(reference))
+        assert geom_equal(Geom.from_file(reference),
+                          modified)
+
+    @pytest.mark.parametrize('directions',[
+                                           ['x'],
+                                           ['x','y','z'],
+                                           ['z','x','y'],
+                                           ['y','z'],
+                                          ]
+                            )
+    def test_flip(self,default,update,reference_dir,directions):
+        modified = default.flip(directions)
+        tag = f'directions={"-".join(directions)}'
+        reference = reference_dir/f'flip_{tag}.geom'
+        if update: modified.to_file(reference)
+        assert geom_equal(Geom.from_file(reference),
+                          modified)
 
     @pytest.mark.parametrize('stencil',[1,2,3,4])
     def test_clean(self,default,update,reference_dir,stencil):
-        modified = default.copy()
-        modified.clean(stencil)
+        modified = default.clean(stencil)
         tag = f'stencil={stencil}'
         reference = reference_dir/f'clean_{tag}.geom'
         if update: modified.to_file(reference)
-        assert geom_equal(modified,Geom.from_file(reference))
+        assert geom_equal(Geom.from_file(reference),
+                          modified)
 
     @pytest.mark.parametrize('grid',[
                                      (10,11,10),
@@ -134,33 +164,29 @@ class TestGeom:
                                     ]
                             )
     def test_scale(self,default,update,reference_dir,grid):
-        modified = default.copy()
-        modified.scale(grid)
+        modified = default.scale(grid)
         tag = f'grid={util.srepr(grid,"-")}'
         reference = reference_dir/f'scale_{tag}.geom'
         if update: modified.to_file(reference)
-        assert geom_equal(modified,Geom.from_file(reference))
+        assert geom_equal(Geom.from_file(reference),
+                          modified)
 
     def test_renumber(self,default):
-        modified = default.copy()
-        microstructure = modified.get_microstructure()
+        microstructure = default.get_microstructure()
         for m in np.unique(microstructure):
             microstructure[microstructure==m] = microstructure.max() + np.random.randint(1,30)
-        modified.update(microstructure)
+        modified = default.duplicate(microstructure)
         assert not geom_equal(modified,default)
-        modified.renumber()
-        assert geom_equal(modified,default)
+        assert geom_equal(default,
+                          modified.renumber())
 
     def test_substitute(self,default):
-        modified = default.copy()
-        microstructure = modified.get_microstructure()
         offset = np.random.randint(1,500)
-        microstructure += offset
-        modified.update(microstructure)
+        modified = default.duplicate(default.get_microstructure() + offset)
         assert not geom_equal(modified,default)
-        modified.substitute(np.arange(default.microstructure.max())+1+offset,
-                            np.arange(default.microstructure.max())+1)
-        assert geom_equal(modified,default)
+        assert geom_equal(default,
+                          modified.substitute(np.arange(default.microstructure.max())+1+offset,
+                                              np.arange(default.microstructure.max())+1))
 
     @pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]),
                                            np.array([1,0,0,175]),np.array([0,-1,0,178]),np.array([0,0,1,180])])
@@ -168,38 +194,35 @@ class TestGeom:
         modified = default.copy()
         for i in range(np.rint(360/axis_angle[3]).astype(int)):
             modified.rotate(Rotation.from_axis_angle(axis_angle,degrees=True))
-        assert geom_equal(modified,default)
+        assert geom_equal(default,modified)
 
     @pytest.mark.parametrize('Eulers',[[32.0,68.0,21.0],
                                        [0.0,32.0,240.0]])
     def test_rotate(self,default,update,reference_dir,Eulers):
-        modified = default.copy()
-        modified.rotate(Rotation.from_Eulers(Eulers,degrees=True))
+        modified = default.rotate(Rotation.from_Eulers(Eulers,degrees=True))
         tag = f'Eulers={util.srepr(Eulers,"-")}'
         reference = reference_dir/f'rotate_{tag}.geom'
         if update: modified.to_file(reference)
-        assert geom_equal(modified,Geom.from_file(reference))
+        assert geom_equal(Geom.from_file(reference),
+                          modified)
 
     def test_canvas(self,default):
+        grid = default.grid
         grid_add = np.random.randint(0,30,(3))
-        modified = default.copy()
-        modified.canvas(modified.grid + grid_add)
-        e = default.grid
-        assert np.all(modified.microstructure[:e[0],:e[1],:e[2]] == default.microstructure)
+        modified = default.canvas(grid + grid_add)
+        assert np.all(modified.microstructure[:grid[0],:grid[1],:grid[2]] == default.microstructure)
 
     @pytest.mark.parametrize('center1,center2',[(np.random.random(3)*.5,np.random.random(3)),
                                                 (np.random.randint(4,8,(3)),np.random.randint(9,12,(3)))])
     @pytest.mark.parametrize('diameter',[np.random.random(3)*.5,
-                                        np.random.randint(4,10,(3))])
+                                         np.random.randint(4,10,(3))])
     def test_add_primitive(self,diameter,center1,center2):
         """Same volume fraction for periodic microstructures and different center."""
         o = np.random.random(3)-.5
         g = np.random.randint(8,32,(3))
         s = np.random.random(3)+.5
-        G_1 = Geom(np.ones(g,'i'),s,o)
-        G_2 = Geom(np.ones(g,'i'),s,o)
-        G_1.add_primitive(diameter,center1,1)
-        G_2.add_primitive(diameter,center2,1)
+        G_1 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center1,1)
+        G_2 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center2,1)
         assert np.count_nonzero(G_1.microstructure!=2) == np.count_nonzero(G_2.microstructure!=2)
 
     @pytest.mark.parametrize('trigger',[[1],[]])
@@ -218,8 +241,7 @@ class TestGeom:
         if len(trigger) > 0:
             m2[m==1] = 1
 
-        geom = Geom(m,np.random.rand(3))
-        geom.vicinity_offset(vicinity,offset,trigger=trigger)
+        geom = Geom(m,np.random.rand(3)).vicinity_offset(vicinity,offset,trigger=trigger)
 
         assert np.all(m2==geom.microstructure)
 

From 3738002e4a77ad9c99f99d603e3bc871de231a79 Mon Sep 17 00:00:00 2001
From: Philip Eisenlohr <eisenlohr@egr.msu.edu>
Date: Sun, 23 Aug 2020 22:27:53 -0400
Subject: [PATCH 08/41] forgot the new reference files...

---
 .../reference/Geom/flip_directions=x-y-z.geom | 25 +++++++++++++++++++
 .../reference/Geom/flip_directions=x.geom     | 25 +++++++++++++++++++
 .../reference/Geom/flip_directions=y-z.geom   | 25 +++++++++++++++++++
 .../reference/Geom/flip_directions=z-x-y.geom | 25 +++++++++++++++++++
 4 files changed, 100 insertions(+)
 create mode 100644 python/tests/reference/Geom/flip_directions=x-y-z.geom
 create mode 100644 python/tests/reference/Geom/flip_directions=x.geom
 create mode 100644 python/tests/reference/Geom/flip_directions=y-z.geom
 create mode 100644 python/tests/reference/Geom/flip_directions=z-x-y.geom

diff --git a/python/tests/reference/Geom/flip_directions=x-y-z.geom b/python/tests/reference/Geom/flip_directions=x-y-z.geom
new file mode 100644
index 000000000..99e55ad7f
--- /dev/null
+++ b/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
diff --git a/python/tests/reference/Geom/flip_directions=x.geom b/python/tests/reference/Geom/flip_directions=x.geom
new file mode 100644
index 000000000..9d4ee74a9
--- /dev/null
+++ b/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
diff --git a/python/tests/reference/Geom/flip_directions=y-z.geom b/python/tests/reference/Geom/flip_directions=y-z.geom
new file mode 100644
index 000000000..ecd22f902
--- /dev/null
+++ b/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
diff --git a/python/tests/reference/Geom/flip_directions=z-x-y.geom b/python/tests/reference/Geom/flip_directions=z-x-y.geom
new file mode 100644
index 000000000..99e55ad7f
--- /dev/null
+++ b/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

From 56157af7db80f85774605901b4e6c9441143aab7 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 24 Aug 2020 19:03:27 +0200
Subject: [PATCH 09/41] updated script for building fedora (rpm) package

---
 PRIVATE | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/PRIVATE b/PRIVATE
index 4715b68a1..1ca2223c6 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 4715b68a1b55d8f18fe6b405e312785489d9c9d0
+Subproject commit 1ca2223c68475bbcb9da633353dbe4a98c18db0d

From 71e08ea66aa1d5c425430acf12c25f8197afcf99 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 24 Aug 2020 23:23:47 +0200
Subject: [PATCH 10/41] small fixes:

- numpy has multiple float variants
- start renaming microstructure -> materialpoint where it does not hurt
- no need for type check if comparing against set of valid directions
  (or I miss some corner cases)
---
 python/damask/_geom.py | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 083e0ec0d..8d98700c4 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -45,12 +45,12 @@ 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'size     x y z:      {util.srepr(self.get_size  ()," x ")}',
-               f'origin   x y z:      {util.srepr(self.get_origin(),"   ")}',
-               f'# microstructures:   {self.N_microstructure}',
-               f'max microstructure:  {np.nanmax(self.microstructure)}',
-              ]+self.get_comments())
+               f'grid     a b c:     {util.srepr(self.get_grid  ()," x ")}',
+               f'size     x y z:     {util.srepr(self.get_size  ()," x ")}',
+               f'origin   x y z:     {util.srepr(self.get_origin(),"   ")}',
+               f'# materialpoints:   {self.N_microstructure}',
+               f'max materialpoint:  {np.nanmax(self.microstructure)}',
+              ])
 
 
     def __copy__(self):
@@ -182,7 +182,7 @@ class Geom:
             else:
                 self.microstructure = np.copy(microstructure)
 
-        if         self.microstructure.dtype == float and \
+        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)
 
@@ -640,9 +640,7 @@ class Geom:
 
         """
         valid = {'x','y','z'}
-        if not all(isinstance(d, str) for d in directions):
-            raise TypeError('Directions are not of type str.')
-        elif not set(directions).issubset(valid):
+        if not set(directions).issubset(valid):
             raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
 
         limits = [None,None]

From 0ceba2a6d31e1653e636c4ed16640232ef31828c Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 24 Aug 2020 23:28:26 +0200
Subject: [PATCH 11/41] use central functionality

---
 python/damask/_table.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/python/damask/_table.py b/python/damask/_table.py
index e8d35c545..fd8cf0fe3 100644
--- a/python/damask/_table.py
+++ b/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())

From 9a2ac3154550731752f145de277d9f205cf1b897 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Mon, 24 Aug 2020 23:35:46 +0200
Subject: [PATCH 12/41] explicit is better then implicit

---
 python/damask/_result.py | 35 +++++++++++++++++------------------
 1 file changed, 17 insertions(+), 18 deletions(-)

diff --git a/python/damask/_result.py b/python/damask/_result.py
index bf3ee67b9..5d0da804d 100644
--- a/python/damask/_result.py
+++ b/python/damask/_result.py
@@ -1,6 +1,5 @@
 import multiprocessing as mp
 import re
-import inspect
 import glob
 import os
 import datetime
@@ -536,7 +535,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):
@@ -564,7 +563,7 @@ 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):
@@ -598,7 +597,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'):
@@ -624,7 +623,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):
@@ -648,7 +647,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):
@@ -679,7 +678,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'):
@@ -712,7 +711,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'):
@@ -745,7 +744,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):
@@ -771,7 +770,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):
@@ -798,7 +797,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):
@@ -834,7 +833,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):
@@ -862,7 +861,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'):
@@ -898,7 +897,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):
@@ -926,7 +925,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):
@@ -950,7 +949,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):
@@ -974,7 +973,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):
@@ -1006,7 +1005,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'):

From 2751cdb6afb4c0c9adcdd8308682d65214448bb4 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 25 Aug 2020 00:40:14 +0200
Subject: [PATCH 13/41] support comments also in vtr

---
 python/damask/_geom.py | 4 +++-
 python/damask/_vtk.py  | 6 +++---
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 8d98700c4..7185f0963 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -359,11 +359,12 @@ class Geom:
 
         """
         v = VTK.from_file(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
+        comments = v.get_comments()
         grid = np.array(v.geom.GetDimensions())-1
         bbox = np.array(v.geom.GetBounds()).reshape(3,2).T
         size = bbox[1] - bbox[0]
 
-        return Geom(v.get('materialpoint').reshape(grid,order='F'),size,bbox[0])
+        return Geom(v.get('materialpoint').reshape(grid,order='F'),size,bbox[0],comments=comments)
 
 
     @staticmethod
@@ -527,6 +528,7 @@ class Geom:
         """
         v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
         v.add(self.microstructure.flatten(order='F'),'materialpoint')
+        v.add_comments(self.comments)
 
         if fname:
             v.write(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py
index 15a571f4f..cccfc5524 100644
--- a/python/damask/_vtk.py
+++ b/python/damask/_vtk.py
@@ -284,7 +284,7 @@ class VTK:
 
     def set_comments(self,comments):
         """
-        Set Comments.
+        Set comments.
 
         Parameters
         ----------
@@ -301,7 +301,7 @@ class VTK:
 
     def add_comments(self,comments):
         """
-        Add Comments.
+        Add comments.
 
         Parameters
         ----------
@@ -309,7 +309,7 @@ class VTK:
             Comments to add.
 
         """
-        self.set_comments(self.get_comments + ([comments] if isinstance(comments,str) else comments))
+        self.set_comments(self.get_comments() + ([comments] if isinstance(comments,str) else comments))
 
 
     def __repr__(self):

From 453f5a14c539d5211c1ee6530a977e312400dd32 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 25 Aug 2020 00:53:05 +0200
Subject: [PATCH 14/41] duplicated test (w/o wait for parallel out)

---
 python/tests/test_Geom.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index d1657f1df..4fe07fe54 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -73,10 +73,6 @@ class TestGeom:
             new = Geom.from_file(f)
         assert geom_equal(default,new)
 
-    def test_export_import_vtk(self,default,tmpdir):
-        default.to_vtr(tmpdir/'default')
-        assert geom_equal(default,Geom.from_vtr(tmpdir/'default.vtr'))
-
     def test_read_write_vtr(self,default,tmpdir):
         default.to_vtr(tmpdir/'default')
         for _ in range(10):

From 2d98325fa4de0dfc741c8ad880c59c7be07d8fb2 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 25 Aug 2020 00:59:41 +0200
Subject: [PATCH 15/41] increasing test coverage

---
 python/damask/_geom.py    | 20 ++++++++++----------
 python/tests/test_Geom.py | 25 ++++++++++++++++++++-----
 2 files changed, 30 insertions(+), 15 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 7185f0963..07d3a73a9 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -112,24 +112,24 @@ class Geom:
         """
         message = []
         if np.any(other.get_grid() != self.get_grid()):
-            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(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 ")}'))
 
         if np.any(other.get_size() != self.get_size()):
-            message.append(util.delete(f'size     x y z:      {util.srepr(other.get_size()," x ")}'))
-            message.append(util.emph(  f'size     x y z:      {util.srepr( self.get_size()," x ")}'))
+            message.append(util.delete(f'size     x y z:     {util.srepr(other.get_size()," x ")}'))
+            message.append(util.emph(  f'size     x y z:     {util.srepr( self.get_size()," x ")}'))
 
         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(),"   ")}'))
+            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'# microstructures:   {other.N_microstructure}'))
-            message.append(util.emph(  f'# microstructures:   { 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 microstructure:  {np.nanmax(other.microstructure)}'))
-            message.append(util.emph(  f'max 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)
 
diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 4fe07fe54..4b7f05f91 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -13,7 +13,8 @@ from damask import util
 def geom_equal(a,b):
     return np.all(a.get_microstructure() == b.get_microstructure()) and \
            np.all(a.get_grid()           == b.get_grid()) and \
-           np.allclose(a.get_size(), b.get_size())
+           np.allclose(a.get_size(), b.get_size()) and \
+           str(a.diff(b)) == str(b.diff(a))
 
 @pytest.fixture
 def default():
@@ -45,6 +46,13 @@ class TestGeom:
         print(modified)
         assert geom_equal(default,modified)
 
+    def test_diff_equal(self,default):
+        assert str(default.diff(default)) == ''
+
+    def test_diff_not_equal(self,default):
+        new = Geom(default.microstructure[1:,1:,1:]+1,default.size*.9,np.ones(3)-default.origin,comments=['modified'])
+        assert str(default.diff(new)) != ''
+
 
     @pytest.mark.parametrize('masked',[True,False])
     def test_set_microstructure(self,default,masked):
@@ -55,8 +63,8 @@ class TestGeom:
 
 
     def test_write_read_str(self,default,tmpdir):
-        default.to_file(tmpdir/'default.geom')
-        new = Geom.from_file(tmpdir/'default.geom')
+        default.to_file(str(tmpdir/'default.geom'))
+        new = Geom.from_file(str(tmpdir/'default.geom'))
         assert geom_equal(default,new)
 
     def test_write_read_file(self,default,tmpdir):
@@ -139,6 +147,12 @@ class TestGeom:
         assert geom_equal(Geom.from_file(reference),
                           modified)
 
+    @pytest.mark.parametrize('directions',[(1,2,'y'),('a','b','x'),[1]])
+    def test_mirror_invalid(self,default,directions):
+        with pytest.raises(ValueError):
+            default.mirror(directions)
+
+
     @pytest.mark.parametrize('directions',[
                                            ['x'],
                                            ['x','y','z'],
@@ -155,9 +169,10 @@ class TestGeom:
                           modified)
 
     @pytest.mark.parametrize('directions',[(1,2,'y'),('a','b','x'),[1]])
-    def test_mirror_invalid(self,default,directions):
+    def test_flip_invalid(self,default,directions):
         with pytest.raises(ValueError):
-            default.mirror(directions)
+            default.flip(directions)
+
 
     @pytest.mark.parametrize('stencil',[1,2,3,4])
     @pytest.mark.parametrize('selection',[None,[1],[1,2,3]])

From 51da63210876708198495fb8de79a3e53232a619 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 25 Aug 2020 07:49:56 +0200
Subject: [PATCH 16/41] aiming at testing each individual statement

---
 python/damask/_vtk.py    | 22 ++++++++++---------
 python/tests/test_VTK.py | 46 +++++++++++++++++++++++++++++++++++-----
 2 files changed, 53 insertions(+), 15 deletions(-)

diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py
index cccfc5524..60e556f84 100644
--- a/python/damask/_vtk.py
+++ b/python/damask/_vtk.py
@@ -131,17 +131,19 @@ 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'):
+                reader.Update()
                 geom = reader.GetRectilinearGridOutput()
             elif dataset_type.lower().endswith('unstructuredgrid'):
+                reader.Update()
                 geom = reader.GetUnstructuredGridOutput()
             elif dataset_type.lower().endswith('polydata'):
+                reader.Update()
                 geom = reader.GetPolyDataOutput()
             else:
                 raise TypeError(f'Unknown dataset type {dataset_type} for vtk file')
@@ -259,15 +261,15 @@ class VTK:
             Data label.
 
         """
-        celldata = self.geom.GetCellData()
-        for a in range(celldata.GetNumberOfArrays()):
-            if celldata.GetArrayName(a) == label:
-                return vtk_to_np(celldata.GetArray(a))
+        cell_data = self.geom.GetCellData()
+        for a in range(cell_data.GetNumberOfArrays()):
+            if cell_data.GetArrayName(a) == label:
+                return vtk_to_np(cell_data.GetArray(a))
 
-        pointdata = self.geom.GetPointData()
-        for a in range(celldata.GetNumberOfArrays()):
-            if pointdata.GetArrayName(a) == label:
-                return vtk_to_np(pointdata.GetArray(a))
+        point_data = self.geom.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')
 
diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py
index 91ff4033c..35155caa2 100644
--- a/python/tests/test_VTK.py
+++ b/python/tests/test_VTK.py
@@ -13,8 +13,19 @@ def reference_dir(reference_dir_base):
     """Directory containing reference results."""
     return reference_dir_base/'VTK'
 
+@pytest.fixture
+def default():
+    """Simple VTK."""
+    grid = np.array([5,6,7],int)
+    size = np.array([.6,1.,.5])
+    return VTK.from_rectilinearGrid(grid,size)
+
 class TestVTK:
 
+    @pytest.fixture(autouse=True)
+    def _execution_stamp(self, execution_stamp):
+        print('patched damask.util.execution_stamp')
+
     def test_rectilinearGrid(self,tmp_path):
         grid   = np.random.randint(5,10,3)*2
         size   = np.random.random(3) + 1.0
@@ -77,9 +88,33 @@ class TestVTK:
     @pytest.mark.parametrize('name,dataset_type',[('this_file_does_not_exist.vtk', None),
                                                   ('this_file_does_not_exist.vtk','vtk'),
                                                   ('this_file_does_not_exist.vtx', None)])
-    def test_invalid_dataset_type(self,dataset_type,name):
+    def test_invalid_dataset_type(self,name,dataset_type):
         with pytest.raises(TypeError):
-            VTK.from_file('this_file_does_not_exist.vtk',dataset_type)
+            VTK.from_file(name,dataset_type)
+
+    def test_invalid_extension_write(self,default):
+        with pytest.raises(ValueError):
+            default.write('default.txt')
+
+    def test_invalid_get(self,default):
+        with pytest.raises(ValueError):
+            default.get('does_not_exist')
+
+    @pytest.mark.parametrize('data,label',[(np.ones(3),'valid'),
+                                           (np.ones(3),None)])
+    def test_invalid_add(self,default,data,label):
+        with pytest.raises(ValueError):
+            default.add(np.ones(3),label)
+
+    def test_invalid_add_type(self,default):
+        with pytest.raises(TypeError):
+            default.add('invalid_type','label')
+
+    def test_comments(self,tmp_path,default):
+        default.add_comments(['this is a comment'])
+        default.write(tmp_path/'with_comments',parallel=False)
+        new = VTK.from_file(tmp_path/'with_comments.vtr')
+        assert new.get_comments() == ['this is a comment']
 
 
     def test_compare_reference_polyData(self,update,reference_dir,tmp_path):
@@ -90,7 +125,8 @@ class TestVTK:
              polyData.write(reference_dir/'polyData')
         else:
              reference = VTK.from_file(reference_dir/'polyData.vtp')
-             assert polyData.__repr__() == reference.__repr__()
+             assert polyData.__repr__() == reference.__repr__() and \
+                    np.allclose(polyData.get('coordinates'),points)
 
     def test_compare_reference_rectilinearGrid(self,update,reference_dir,tmp_path):
         grid = np.array([5,6,7],int)
@@ -104,5 +140,5 @@ class TestVTK:
              rectilinearGrid.write(reference_dir/'rectilinearGrid')
         else:
              reference = VTK.from_file(reference_dir/'rectilinearGrid.vtr')
-             assert rectilinearGrid.__repr__() == reference.__repr__()
-
+             assert rectilinearGrid.__repr__() == reference.__repr__() and \
+                    np.allclose(rectilinearGrid.get('cell'),c)

From 96d4975fc4e0349cc7f82261f5fcdb0361514dc7 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 25 Aug 2020 08:34:04 +0200
Subject: [PATCH 17/41] using numpy functionality

---
 python/damask/_geom.py    | 16 ++++------------
 python/tests/test_Geom.py |  3 +++
 2 files changed, 7 insertions(+), 12 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 07d3a73a9..1dac0fa84 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -608,10 +608,10 @@ class Geom:
             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 set(directions).issubset(valid):
             raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
 
@@ -641,19 +641,11 @@ class Geom:
             Valid entries are 'x', 'y', 'z'.
 
         """
-        valid = {'x','y','z'}
+        valid = ['x','y','z']
         if not set(directions).issubset(valid):
             raise ValueError(f'Invalid direction {set(directions).difference(valid)} specified.')
 
-        limits = [None,None]
-        ms = self.get_microstructure()
-
-        if 'z' in directions:
-            ms = ms[:,:,limits[0]:limits[1]:-1]
-        if 'y' in directions:
-            ms = ms[:,limits[0]:limits[1]:-1,:]
-        if 'x' in directions:
-            ms = ms[limits[0]:limits[1]:-1,:,:]
+        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')],
diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 4b7f05f91..68d727f11 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -168,6 +168,9 @@ class TestGeom:
         assert geom_equal(Geom.from_file(reference),
                           modified)
 
+    def test_flip_invariant(self,default):
+        assert geom_equal(default,default.flip([]))
+
     @pytest.mark.parametrize('directions',[(1,2,'y'),('a','b','x'),[1]])
     def test_flip_invalid(self,default,directions):
         with pytest.raises(ValueError):

From 800dac5d0171fb57bd92b75e6bb251a0e08d38d1 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 25 Aug 2020 08:46:08 +0200
Subject: [PATCH 18/41] correct type handling

no reason to calculate fill twice
---
 python/damask/_geom.py | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 1dac0fa84..8daf5f476 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -769,16 +769,14 @@ class Geom:
         offset : numpy.ndarray of shape (3)
             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,
-                         np.nanmax(self.microstructure)+1 if fill is None else fill,
-                         dtype)
+        canvas = np.full(self.grid if grid is None else grid,fill,dtype)
 
         LL = np.clip( offset,          0,np.minimum(self.grid,     grid+offset))
         UR = np.clip( offset+grid,     0,np.minimum(self.grid,     grid+offset))

From bf401e56cfc9c4fe51d71f64496d8207cf7caf8b Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 25 Aug 2020 09:56:24 +0200
Subject: [PATCH 19/41] avoid name clash with damask.Geom

---
 python/damask/_geom.py |  4 +--
 python/damask/_vtk.py  | 77 +++++++++++++++++++++---------------------
 2 files changed, 41 insertions(+), 40 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 8daf5f476..c9848e641 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -360,8 +360,8 @@ class Geom:
         """
         v = VTK.from_file(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
         comments = v.get_comments()
-        grid = np.array(v.geom.GetDimensions())-1
-        bbox = np.array(v.geom.GetBounds()).reshape(3,2).T
+        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)
diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py
index 60e556f84..bbc2d2e69 100644
--- a/python/damask/_vtk.py
+++ b/python/damask/_vtk.py
@@ -20,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
-            Description of geometry and topology. Valid types are vtk.vtkRectilinearGrid,
-            vtk.vtkUnstructuredGrid, or vtk.vtkPolyData.
+        vtk_data : subclass of vtk.vtkDataSet
+            Description of geometry and topology, optionally with attached data.
+            Valid types are vtk.vtkRectilinearGrid, vtk.vtkUnstructuredGrid,
+            or vtk.vtkPolyData.
 
         """
-        self.geom = geom
+        self.vtk_data = vtk_data
 
 
     @staticmethod
@@ -51,15 +52,15 @@ class VTK:
             Spatial origin.
 
         """
-        geom = vtk.vtkRectilinearGrid()
-        geom.SetDimensions(*(grid+1))
+        vtk_data = vtk.vtkRectilinearGrid()
+        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])
-        geom.SetYCoordinates(coord[1])
-        geom.SetZCoordinates(coord[2])
+        vtk_data.SetXCoordinates(coord[0])
+        vtk_data.SetYCoordinates(coord[1])
+        vtk_data.SetZCoordinates(coord[2])
 
-        return VTK(geom)
+        return VTK(vtk_data)
 
 
     @staticmethod
@@ -87,11 +88,11 @@ class VTK:
                             connectivity),axis=1).ravel()
         cells.SetCells(connectivity.shape[0],np_to_vtkIdTypeArray(T,deep=True))
 
-        geom = vtk.vtkUnstructuredGrid()
-        geom.SetPoints(vtk_nodes)
-        geom.SetCells(eval(f'vtk.VTK_{cell_type.split("_",1)[-1].upper()}'),cells)
+        vtk_data = vtk.vtkUnstructuredGrid()
+        vtk_data.SetPoints(vtk_nodes)
+        vtk_data.SetCells(eval(f'vtk.VTK_{cell_type.split("_",1)[-1].upper()}'),cells)
 
-        return VTK(geom)
+        return VTK(vtk_data)
 
 
     @staticmethod
@@ -110,10 +111,10 @@ class VTK:
         vtk_points = vtk.vtkPoints()
         vtk_points.SetData(np_to_vtk(points))
 
-        geom = vtk.vtkPolyData()
-        geom.SetPoints(vtk_points)
+        vtk_data = vtk.vtkPolyData()
+        vtk_data.SetPoints(vtk_points)
 
-        return VTK(geom)
+        return VTK(vtk_data)
 
 
     @staticmethod
@@ -138,13 +139,13 @@ class VTK:
                 raise TypeError('Dataset type for *.vtk file not given.')
             elif dataset_type.lower().endswith('rectilineargrid'):
                 reader.Update()
-                geom = reader.GetRectilinearGridOutput()
+                vtk_data = reader.GetRectilinearGridOutput()
             elif dataset_type.lower().endswith('unstructuredgrid'):
                 reader.Update()
-                geom = reader.GetUnstructuredGridOutput()
+                vtk_data = reader.GetUnstructuredGridOutput()
             elif dataset_type.lower().endswith('polydata'):
                 reader.Update()
-                geom = reader.GetPolyDataOutput()
+                vtk_data = reader.GetPolyDataOutput()
             else:
                 raise TypeError(f'Unknown dataset type {dataset_type} for vtk file')
         else:
@@ -159,9 +160,9 @@ 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):
@@ -179,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()
@@ -193,7 +194,7 @@ 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:
             try:
@@ -220,8 +221,8 @@ class VTK:
             Data label.
 
         """
-        N_points = self.geom.GetNumberOfPoints()
-        N_cells  = self.geom.GetNumberOfCells()
+        N_points = self.vtk_data.GetNumberOfPoints()
+        N_cells  = self.vtk_data.GetNumberOfCells()
 
         if   isinstance(data,np.ndarray):
             if label is None:
@@ -234,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):
@@ -261,12 +262,12 @@ class VTK:
             Data label.
 
         """
-        cell_data = self.geom.GetCellData()
+        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.geom.GetPointData()
+        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))
@@ -276,7 +277,7 @@ class VTK:
 
     def get_comments(self):
         """Return the comments."""
-        fielddata = self.geom.GetFieldData()
+        fielddata = self.vtk_data.GetFieldData()
         for a in range(fielddata.GetNumberOfArrays()):
             if fielddata.GetArrayName(a) == 'comments':
                 comments = fielddata.GetAbstractArray(a)
@@ -298,7 +299,7 @@ class VTK:
         s.SetName('comments')
         for c in [comments] if isinstance(comments,str) else comments:
             s.InsertNextValue(c)
-        self.geom.GetFieldData().AddArray(s)
+        self.vtk_data.GetFieldData().AddArray(s)
 
 
     def add_comments(self,comments):
@@ -319,7 +320,7 @@ class VTK:
         writer = vtk.vtkDataSetWriter()
         writer.SetHeader(f'# {util.execution_stamp("VTK")}')
         writer.WriteToOutputStringOn()
-        writer.SetInputData(self.geom)
+        writer.SetInputData(self.vtk_data)
         writer.Write()
         return writer.GetOutputString()
 
@@ -331,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)
 

From fee21cbd9c0c56f34d8204fc88a474b3ba93575f Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 25 Aug 2020 17:17:49 +0200
Subject: [PATCH 20/41] more sensible tests

---
 python/tests/test_Geom.py |  4 ++++
 python/tests/test_VTK.py  | 14 ++++++++------
 2 files changed, 12 insertions(+), 6 deletions(-)

diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 68d727f11..1639d9819 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -171,6 +171,10 @@ class TestGeom:
     def test_flip_invariant(self,default):
         assert geom_equal(default,default.flip([]))
 
+    @pytest.mark.parametrize('direction',[['x'],['x','y']])
+    def test_flip_double(self,default,direction):
+        assert geom_equal(default,default.flip(direction).flip(direction))
+
     @pytest.mark.parametrize('directions',[(1,2,'y'),('a','b','x'),[1]])
     def test_flip_invalid(self,default,directions):
         with pytest.raises(ValueError):
diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py
index 35155caa2..ab9c4fa8b 100644
--- a/python/tests/test_VTK.py
+++ b/python/tests/test_VTK.py
@@ -100,15 +100,18 @@ class TestVTK:
         with pytest.raises(ValueError):
             default.get('does_not_exist')
 
-    @pytest.mark.parametrize('data,label',[(np.ones(3),'valid'),
-                                           (np.ones(3),None)])
-    def test_invalid_add(self,default,data,label):
+    def test_invalid_add_shape(self,default):
         with pytest.raises(ValueError):
-            default.add(np.ones(3),label)
+            default.add(np.ones(3),'valid')
+
+    def test_invalid_add_missing_label(self,default):
+        data = np.random.randint(9,size=np.prod(np.array(default.vtk_data.GetDimensions())-1))
+        with pytest.raises(ValueError):
+            default.add(data)
 
     def test_invalid_add_type(self,default):
         with pytest.raises(TypeError):
-            default.add('invalid_type','label')
+            default.add('invalid_type','valid')
 
     def test_comments(self,tmp_path,default):
         default.add_comments(['this is a comment'])
@@ -116,7 +119,6 @@ class TestVTK:
         new = VTK.from_file(tmp_path/'with_comments.vtr')
         assert new.get_comments() == ['this is a comment']
 
-
     def test_compare_reference_polyData(self,update,reference_dir,tmp_path):
         points=np.dstack((np.linspace(0.,1.,10),np.linspace(0.,2.,10),np.linspace(-1.,1.,10))).squeeze()
         polyData = VTK.from_polyData(points)

From 7dba9518cc81f1986362b072a9b07eaf7c87c36e Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 25 Aug 2020 18:08:23 +0200
Subject: [PATCH 21/41] adjust to new Geom behavior

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

diff --git a/processing/pre/geom_grainGrowth.py b/processing/pre/geom_grainGrowth.py
index bdf8d8efe..6d1985ee3 100755
--- a/processing/pre/geom_grainGrowth.py
+++ b/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)

From 06b524d13e267aeb88715e7d8ee7f4670f6c3c60 Mon Sep 17 00:00:00 2001
From: Philip Eisenlohr <eisenlohr@egr.msu.edu>
Date: Wed, 26 Aug 2020 21:39:20 +0200
Subject: [PATCH 22/41] added 'iso' as possible lattice (with no symmetry)

---
 python/damask/_lattice.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py
index 9df451c69..da8f76b64 100644
--- a/python/damask/_lattice.py
+++ b/python/damask/_lattice.py
@@ -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'},

From da46e5ea9a79ef922deb39350991d14391d3b590 Mon Sep 17 00:00:00 2001
From: Philip Eisenlohr <eisenlohr@egr.msu.edu>
Date: Wed, 26 Aug 2020 15:32:48 -0400
Subject: [PATCH 23/41] cannot use np.linalg.norm when dealing with exponent
 triple

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

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index c9848e641..08b0b8757 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -585,7 +585,7 @@ class Geom:
         coords_rot = R.broadcast_to(tuple(self.grid))@coords
 
         with np.errstate(over='ignore',under='ignore'):
-            mask = np.where(np.linalg.norm(coords_rot/r,2.0**exponent,axis=-1) <= 1.0,False,True)
+            mask = np.where(np.sum(np.power(coords_rot/r,2.0**exponent),axis=-1) <= 1.0,False,True)
 
         if periodic:                                                                                # translate back to center
             mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))

From b9f1421c6bd9c8edd06d542e041fbd7f1f8a63dd Mon Sep 17 00:00:00 2001
From: Philip Eisenlohr <eisenlohr@egr.msu.edu>
Date: Wed, 26 Aug 2020 17:27:08 -0400
Subject: [PATCH 24/41] [skip ci] removed unnecessary line continuations "\"

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

diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py
index da8f76b64..143fa50f1 100644
--- a/python/damask/_lattice.py
+++ b/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) * \
-                                (2. >= np.sqrt(3)*rho_abs[...,0] + rho_abs[...,1]) * \
-                                (2. >= np.sqrt(3)*rho_abs[...,1] + rho_abs[...,0]) * \
+                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[...,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) * \
-                                (np.sqrt(2) >= rho_abs[...,0] + rho_abs[...,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[...,2] + 1.),True,False)
             elif self.system == 'orthorhombic':
                 return np.where(np.prod(1.             >= rho_abs,axis=-1),True,False)

From 229f6139c837d79321b13bc11a8039dcfafa9a5b Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 26 Aug 2020 23:54:56 +0200
Subject: [PATCH 25/41] better user experience

- no meaningless warnings regarding floating point precision
- meaningful error message for invalid header
---
 python/damask/_geom.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 08b0b8757..1570c3a9f 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -300,13 +300,16 @@ class Geom:
             f = fname
 
         f.seek(0)
-        header_length,keyword = f.readline().split()[:2]
-        header_length = int(header_length)
-        content = f.readlines()
-
+        try:
+            header_length,keyword = f.readline().split()[:2]
+            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()
@@ -584,7 +587,7 @@ class Geom:
                - ((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
         coords_rot = R.broadcast_to(tuple(self.grid))@coords
 
-        with np.errstate(over='ignore',under='ignore'):
+        with np.errstate(over='ignore',under='ignore',invalid='ignore'):
             mask = np.where(np.sum(np.power(coords_rot/r,2.0**exponent),axis=-1) <= 1.0,False,True)
 
         if periodic:                                                                                # translate back to center

From 248ef8ef97f23ba7075ba026d554a3cedf83d75c Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 26 Aug 2020 23:56:20 +0200
Subject: [PATCH 26/41] more tests

- invalid files
- scalar/array size,center,exponent for add_primitive
- rotation invariance for add primitive
---
 python/tests/test_Geom.py | 40 +++++++++++++++++++++++++++++++++------
 1 file changed, 34 insertions(+), 6 deletions(-)

diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 1639d9819..665aa5d3b 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -69,7 +69,7 @@ class TestGeom:
 
     def test_write_read_file(self,default,tmpdir):
         with open(tmpdir/'default.geom','w') as f:
-            default.to_file(f)
+            default.to_file(f,pack=True)
         with open(tmpdir/'default.geom') as f:
             new = Geom.from_file(f)
         assert geom_equal(default,new)
@@ -83,6 +83,7 @@ class TestGeom:
 
     def test_read_write_vtr(self,default,tmpdir):
         default.to_vtr(tmpdir/'default')
+        print(default.to_vtr())
         for _ in range(10):
             time.sleep(.2)
             if os.path.exists(tmpdir/'default.vtr'): break
@@ -90,6 +91,13 @@ class TestGeom:
         new = Geom.from_vtr(tmpdir/'default.vtr')
         assert geom_equal(new,default)
 
+    def test_invalid_geom(self,tmpdir):
+        with open('invalid_file','w') as f:
+            f.write('this is not a valid header')
+        with open('invalid_file','r') as f:
+            with pytest.raises(TypeError):
+                Geom.from_file(f)
+
     def test_invalid_vtr(self,tmpdir):
         v = VTK.from_rectilinearGrid(np.random.randint(5,10,3)*2,np.random.random(3) + 1.0)
         v.write(tmpdir/'no_materialpoint.vtr')
@@ -254,19 +262,39 @@ class TestGeom:
         modified = default.canvas(grid + grid_add)
         assert np.all(modified.microstructure[:grid[0],:grid[1],:grid[2]] == default.microstructure)
 
-    @pytest.mark.parametrize('center1,center2',[(np.random.random(3)*.5,np.random.random(3)),
+    @pytest.mark.parametrize('center1,center2',[(np.random.random(3)*.5,np.random.random()*8),
                                                 (np.random.randint(4,8,(3)),np.random.randint(9,12,(3)))])
     @pytest.mark.parametrize('diameter',[np.random.random(3)*.5,
-                                         np.random.randint(4,10,(3))])
-    def test_add_primitive(self,diameter,center1,center2):
+                                         np.random.randint(4,10,(3)),
+                                         np.random.rand(),
+                                         np.random.randint(30)])
+    @pytest.mark.parametrize('exponent',[np.random.random(3)*.5,
+                                         np.random.randint(4,10,(3)),
+                                         np.random.rand()*4,
+                                         np.random.randint(20)])
+    def test_add_primitive(self,center1,center2,diameter,exponent):
         """Same volume fraction for periodic microstructures and different center."""
         o = np.random.random(3)-.5
         g = np.random.randint(8,32,(3))
         s = np.random.random(3)+.5
-        G_1 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center1,1)
-        G_2 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center2,1)
+        G_1 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center1,exponent)
+        G_2 = Geom(np.ones(g,'i'),s,o).add_primitive(diameter,center2,exponent)
         assert np.count_nonzero(G_1.microstructure!=2) == np.count_nonzero(G_2.microstructure!=2)
 
+    @pytest.mark.parametrize('center',[np.random.randint(4,10,(3)),
+                                       np.random.randint(2,10),
+                                       np.random.rand()*4,
+                                       np.random.rand(3)*10])
+    @pytest.mark.parametrize('inverse',[True,False])
+    @pytest.mark.parametrize('periodic',[True,False])
+    def test_add_primitive_rotation(self,center,inverse,periodic):
+        g = np.random.randint(8,32,(3))
+        s = np.random.random()+.5
+        fill = np.random.randint(10)+2
+        G_1 = Geom(np.ones(g,'i'),[s,s,s]).add_primitive(s*.3,center,1,fill,inverse=inverse,periodic=periodic)
+        G_2 = Geom(np.ones(g,'i'),[s,s,s]).add_primitive(s*.3,center,1,fill,Rotation.from_random(),inverse,periodic=periodic)
+        assert geom_equal(G_1,G_2)
+
     @pytest.mark.parametrize('trigger',[[1],[]])
     def test_vicinity_offset(self,trigger):
         offset = np.random.randint(2,4)

From c2191fd930bd8db39b542be977831ff4b1539d01 Mon Sep 17 00:00:00 2001
From: Philip Eisenlohr <eisenlohr@egr.msu.edu>
Date: Thu, 27 Aug 2020 00:06:52 +0200
Subject: [PATCH 27/41] Make logic to mask (outside of) primitive more apparent
 (hopefully).

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

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index 1570c3a9f..c5833af6c 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -588,7 +588,7 @@ class Geom:
         coords_rot = R.broadcast_to(tuple(self.grid))@coords
 
         with np.errstate(over='ignore',under='ignore',invalid='ignore'):
-            mask = np.where(np.sum(np.power(coords_rot/r,2.0**exponent),axis=-1) <= 1.0,False,True)
+            mask = np.where(np.sum(np.power(coords_rot/r,2.0**exponent),axis=-1) > 1.0,True,False)
 
         if periodic:                                                                                # translate back to center
             mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))

From c5761831e23d17ac7cabed63ef646371bff951ec Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 27 Aug 2020 00:14:37 +0200
Subject: [PATCH 28/41] more meaningful message

if add_xxx (Results) does not find matching datasets, inform the user
about this fact instead of saying
TypeError: object of type 'IMapUnorderedIterator' has no len()
---
 python/damask/_result.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/python/damask/_result.py b/python/damask/_result.py
index 8e5bb2e8b..2ada4cd01 100644
--- a/python/damask/_result.py
+++ b/python/damask/_result.py
@@ -1065,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)):

From 94797f832580d55e2e5adc3a214bffee381eecb4 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 27 Aug 2020 09:32:49 +0200
Subject: [PATCH 29/41] ignore all kinds of FPE

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

diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index c5833af6c..8a5dfcdc9 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -587,7 +587,7 @@ class Geom:
                - ((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
         coords_rot = R.broadcast_to(tuple(self.grid))@coords
 
-        with np.errstate(over='ignore',under='ignore',invalid='ignore'):
+        with np.errstate(all='ignore'):
             mask = np.where(np.sum(np.power(coords_rot/r,2.0**exponent),axis=-1) > 1.0,True,False)
 
         if periodic:                                                                                # translate back to center

From 0d5279863db1682159ae69492b3372e5d2d97f85 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 27 Aug 2020 09:33:09 +0200
Subject: [PATCH 30/41] avoid rounding errors related to rotation

---
 python/tests/test_Geom.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 665aa5d3b..423f1b961 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -272,7 +272,7 @@ class TestGeom:
                                          np.random.randint(4,10,(3)),
                                          np.random.rand()*4,
                                          np.random.randint(20)])
-    def test_add_primitive(self,center1,center2,diameter,exponent):
+    def test_add_primitive_shift(self,center1,center2,diameter,exponent):
         """Same volume fraction for periodic microstructures and different center."""
         o = np.random.random(3)-.5
         g = np.random.randint(8,32,(3))
@@ -288,11 +288,12 @@ class TestGeom:
     @pytest.mark.parametrize('inverse',[True,False])
     @pytest.mark.parametrize('periodic',[True,False])
     def test_add_primitive_rotation(self,center,inverse,periodic):
-        g = np.random.randint(8,32,(3))
-        s = np.random.random()+.5
+        """Rotation should not chage result for sphere (avoid discretization errors."""
+        g = np.array([32,32,32])
         fill = np.random.randint(10)+2
-        G_1 = Geom(np.ones(g,'i'),[s,s,s]).add_primitive(s*.3,center,1,fill,inverse=inverse,periodic=periodic)
-        G_2 = Geom(np.ones(g,'i'),[s,s,s]).add_primitive(s*.3,center,1,fill,Rotation.from_random(),inverse,periodic=periodic)
+        eu=np.array([np.random.randint(4),np.random.randint(2),np.random.randint(4)])*.5*np.pi
+        G_1 = Geom(np.ones(g,'i'),[1.,1.,1.]).add_primitive(.3,center,1,fill,inverse=inverse,periodic=periodic)
+        G_2 = Geom(np.ones(g,'i'),[1.,1.,1.]).add_primitive(.3,center,1,fill,Rotation.from_Eulers(eu),inverse,periodic=periodic)
         assert geom_equal(G_1,G_2)
 
     @pytest.mark.parametrize('trigger',[[1],[]])

From 4090fdb1dd53cc7b147528ff33e7a10dd52f0894 Mon Sep 17 00:00:00 2001
From: Test User <damask_user@mpie.de>
Date: Thu, 27 Aug 2020 15:51:27 +0200
Subject: [PATCH 31/41] [skip ci] updated version information after successful
 test of v3.0.0-alpha-32-g792db297

---
 VERSION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/VERSION b/VERSION
index 9aba6acf5..7fd44eedf 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-v3.0.0-alpha-27-g68c2908b
+v3.0.0-alpha-32-g792db297

From 338281e0abf5d6e3f70acd0f3df93d92ef04c70f Mon Sep 17 00:00:00 2001
From: Test User <damask_user@mpie.de>
Date: Thu, 27 Aug 2020 19:32:28 +0200
Subject: [PATCH 32/41] [skip ci] updated version information after successful
 test of v3.0.0-alpha-41-g94574356

---
 VERSION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/VERSION b/VERSION
index 7fd44eedf..73dd90e38 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-v3.0.0-alpha-32-g792db297
+v3.0.0-alpha-41-g94574356

From c58693328e7fdd2b71c9b5caf830c78b16b64ecb Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Fri, 28 Aug 2020 12:45:41 +0200
Subject: [PATCH 33/41] don't try to access folders that don't exist

this happened in the case when different types of output exist in
different phases, e.g.

Aluminum: generic, plastic, sources
Steel: generic, plastic

We are a little bit inconsistent because 'generic' and 'plastic' are
always created (even if empty) but 'sources' will only exist if it
contains output. In future, we should have only folders that actually
contain data
---
 python/damask/_result.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/python/damask/_result.py b/python/damask/_result.py
index 2ada4cd01..d739db8f5 100644
--- a/python/damask/_result.py
+++ b/python/damask/_result.py
@@ -401,8 +401,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 len(set(match)) == len(sets): groups.append(group)
+                                if group in f.keys():
+                                    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
 
 

From 03b02ad641c428729b9e6334919b155fe05960cc Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 29 Aug 2020 16:06:10 +0200
Subject: [PATCH 34/41] typos/sentence

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

diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 423f1b961..91cb61591 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -288,7 +288,7 @@ class TestGeom:
     @pytest.mark.parametrize('inverse',[True,False])
     @pytest.mark.parametrize('periodic',[True,False])
     def test_add_primitive_rotation(self,center,inverse,periodic):
-        """Rotation should not chage result for sphere (avoid discretization errors."""
+        """Rotation should not change result for sphere (except for discretization errors)."""
         g = np.array([32,32,32])
         fill = np.random.randint(10)+2
         eu=np.array([np.random.randint(4),np.random.randint(2),np.random.randint(4)])*.5*np.pi

From 44bb99c57a3048723fb52222e21683aaa6506ecc Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 29 Aug 2020 16:34:19 +0200
Subject: [PATCH 35/41] not needed

---
 src/material.f90 | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/src/material.f90 b/src/material.f90
index c89113311..efdeeafcb 100644
--- a/src/material.f90
+++ b/src/material.f90
@@ -241,7 +241,6 @@ subroutine material_parseHomogenization
     homogDamage
 
   integer :: h
-  logical, dimension(:), allocatable :: homogenization_active
 
   material_homogenization => material_root%get('homogenization')
   material_Nhomogenization = material_homogenization%length
@@ -253,13 +252,9 @@ subroutine material_parseHomogenization
   allocate(thermal_typeInstance(material_Nhomogenization),          source=0)
   allocate(damage_typeInstance(material_Nhomogenization),           source=0)
   allocate(homogenization_Ngrains(material_Nhomogenization),        source=0)
-  allocate(homogenization_active(material_Nhomogenization),         source=.false.)  !!!!!!!!!!!!!!!
   allocate(thermal_initialT(material_Nhomogenization),              source=300.0_pReal)
   allocate(damage_initialPhi(material_Nhomogenization),             source=1.0_pReal)
 
-  forall (h = 1:material_Nhomogenization) &
-    homogenization_active(h) = any(discretization_homogenizationAt == h)  !ToDo: SR: needed??
-
   do h=1, material_Nhomogenization
     homog => material_homogenization%get(h)
     homogMech => homog%get('mech')
@@ -317,7 +312,7 @@ subroutine material_parseHomogenization
     damage_typeInstance(h)          = count(damage_type        (1:h) == damage_type        (h))
   enddo
 
-  homogenization_maxNgrains = maxval(homogenization_Ngrains,homogenization_active)
+  homogenization_maxNgrains = maxval(homogenization_Ngrains)
 
 
 end subroutine material_parseHomogenization

From 23bf51cca0176664abdfe12912ef1f8e4c242cbc Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 29 Aug 2020 16:45:18 +0200
Subject: [PATCH 36/41] homogenizationAt from discretization not needed anymore

---
 src/discretization.f90           |  5 +----
 src/grid/discretization_grid.f90 | 20 +++++---------------
 src/marc/discretization_marc.f90 | 29 ++++++++++++-----------------
 src/material.f90                 |  2 +-
 src/mesh/discretization_mesh.f90 |  4 +---
 5 files changed, 20 insertions(+), 40 deletions(-)

diff --git a/src/discretization.f90 b/src/discretization.f90
index e52784c6d..4520bf2ca 100644
--- a/src/discretization.f90
+++ b/src/discretization.f90
@@ -15,7 +15,6 @@ module discretization
     discretization_nElem
     
   integer,     public, protected, dimension(:),   allocatable :: &
-    discretization_homogenizationAt, &
     discretization_microstructureAt   
 
   real(pReal), public, protected, dimension(:,:), allocatable :: & 
@@ -38,12 +37,11 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief stores the relevant information in globally accesible variables
 !--------------------------------------------------------------------------------------------------
-subroutine discretization_init(homogenizationAt,microstructureAt,&
+subroutine discretization_init(microstructureAt,&
                                IPcoords0,NodeCoords0,&
                                sharedNodesBegin)
 
   integer,     dimension(:),   intent(in) :: &
-    homogenizationAt, &
     microstructureAt
   real(pReal), dimension(:,:), intent(in) :: &
     IPcoords0, &
@@ -56,7 +54,6 @@ subroutine discretization_init(homogenizationAt,microstructureAt,&
   discretization_nElem = size(microstructureAt,1)
   discretization_nIP   = size(IPcoords0,2)/discretization_nElem
 
-  discretization_homogenizationAt = homogenizationAt
   discretization_microstructureAt = microstructureAt  
 
   discretization_IPcoords0   = IPcoords0
diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90
index 6d128eb0b..995f9e104 100644
--- a/src/grid/discretization_grid.f90
+++ b/src/grid/discretization_grid.f90
@@ -53,8 +53,7 @@ subroutine discretization_grid_init(restart)
     myGrid                                                                                          !< domain grid of this process
 
   integer,     dimension(:),   allocatable :: &
-    microstructureAt, &
-    homogenizationAt
+    microstructureAt
 
   integer :: &
     j, &
@@ -65,7 +64,7 @@ subroutine discretization_grid_init(restart)
 
   write(6,'(/,a)') ' <<<+-  discretization_grid init  -+>>>'; flush(6)
 
-  call readGeom(grid,geomSize,origin,microstructureAt,homogenizationAt)
+  call readGeom(grid,geomSize,origin,microstructureAt)
 
 !--------------------------------------------------------------------------------------------------
 ! grid solver specific quantities
@@ -94,10 +93,8 @@ subroutine discretization_grid_init(restart)
 ! general discretization
   microstructureAt = microstructureAt(product(grid(1:2))*grid3Offset+1: &
                                       product(grid(1:2))*(grid3Offset+grid3))                       ! reallocate/shrink in case of MPI
-  homogenizationAt = homogenizationAt(product(grid(1:2))*grid3Offset+1: &
-                                      product(grid(1:2))*(grid3Offset+grid3))                       ! reallocate/shrink in case of MPI
 
-  call discretization_init(homogenizationAt,microstructureAt, &
+  call discretization_init(microstructureAt, &
                            IPcoordinates0(myGrid,mySize,grid3Offset), &
                            Nodes0(myGrid,mySize,grid3Offset),&
                            merge((grid(1)+1) * (grid(2)+1) * (grid3+1),&                            ! write bottom layer
@@ -139,7 +136,7 @@ end subroutine discretization_grid_init
 !> @details important variables have an implicit "save" attribute. Therefore, this function is
 ! supposed to be called only once!
 !--------------------------------------------------------------------------------------------------
-subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
+subroutine readGeom(grid,geomSize,origin,microstructure)
 
   integer,     dimension(3), intent(out) :: &
     grid                                                                                            ! grid (for all processes!)
@@ -147,8 +144,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
     geomSize, &                                                                                     ! size (for all processes!)
     origin                                                                                          ! origin (for all processes!)
   integer,     dimension(:), intent(out), allocatable :: &
-    microstructure, &
-    homogenization
+    microstructure
 
   character(len=:),      allocatable :: rawData
   character(len=65536)               :: line
@@ -249,24 +245,18 @@ subroutine readGeom(grid,geomSize,origin,microstructure,homogenization)
           enddo
         endif
 
-      case ('homogenization')
-        if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2)
-
     end select
 
   enddo
 
 !--------------------------------------------------------------------------------------------------
 ! sanity checks
-  if(h < 1) &
-    call IO_error(error_ID = 842, ext_msg='homogenization (readGeom)')
   if(any(grid < 1)) &
     call IO_error(error_ID = 842, ext_msg='grid (readGeom)')
   if(any(geomSize < 0.0_pReal)) &
     call IO_error(error_ID = 842, ext_msg='size (readGeom)')
 
   allocate(microstructure(product(grid)), source = -1)                                              ! too large in case of MPI (shrink later, not very elegant)
-  allocate(homogenization(product(grid)), source = h)                                               ! too large in case of MPI (shrink later, not very elegant)
 
 !--------------------------------------------------------------------------------------------------
 ! read and interpret content
diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90
index e5c382fe1..618e3195b 100644
--- a/src/marc/discretization_marc.f90
+++ b/src/marc/discretization_marc.f90
@@ -52,8 +52,7 @@ subroutine discretization_marc_init
   type(tElement) :: elem
 
   integer,     dimension(:),   allocatable :: &
-    microstructureAt, &
-    homogenizationAt
+    microstructureAt
   integer:: &
     Nnodes, &                                                                                       !< total number of nodes in the mesh
     Nelems, &                                                                                       !< total number of elements in the mesh
@@ -84,7 +83,7 @@ subroutine discretization_marc_init
   mesh_unitlength = num_commercialFEM%get_asFloat('unitlength',defaultVal=1.0_pReal)                ! set physical extent of a length unit in mesh
   if (mesh_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength')
 
-  call inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt)
+  call inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
   nElems = size(connectivity_elem,2)
 
   if (debug_e < 1 .or. debug_e > nElems)    call IO_error(602,ext_msg='element')
@@ -104,7 +103,7 @@ subroutine discretization_marc_init
   call buildIPcoordinates(IP_reshaped,reshape(connectivity_cell,[elem%NcellNodesPerCell,&
                           elem%nIPs*nElems]),node0_cell)
 
-  call discretization_init(microstructureAt,homogenizationAt,&
+  call discretization_init(microstructureAt,&
                            IP_reshaped,&
                            node0_cell)
 
@@ -173,7 +172,7 @@ end subroutine writeGeometry
 !--------------------------------------------------------------------------------------------------
 !> @brief Read mesh from marc input file
 !--------------------------------------------------------------------------------------------------
-subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogenizationAt)
+subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
 
   type(tElement), intent(out) :: elem
   real(pReal), dimension(:,:), allocatable, intent(out) :: &
@@ -181,8 +180,7 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
   integer, dimension(:,:),     allocatable, intent(out) :: &
     connectivity_elem
   integer,     dimension(:),   allocatable, intent(out) :: &
-    microstructureAt, &
-    homogenizationAt
+    microstructureAt
 
   integer :: &
     fileFormatVersion, &
@@ -228,9 +226,9 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
 
   connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,inputFile)
 
-  call inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, &
-                                                 nElems,elem%nNodes,nameElemSet,mapElemSet,&
-                                                 initialcondTableStyle,inputFile)
+  call inputRead_microstructureAnd(microstructureAt, &
+                                   nElems,elem%nNodes,nameElemSet,mapElemSet,&
+                                   initialcondTableStyle,inputFile)
 end subroutine inputRead
 
 
@@ -677,14 +675,13 @@ end function inputRead_connectivityElem
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Stores homogenization and microstructure ID
+!> @brief Store microstructure ID
 !--------------------------------------------------------------------------------------------------
-subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogenizationAt, &
+subroutine inputRead_microstructure(microstructureAt,&
                                     nElem,nNodes,nameElemSet,mapElemSet,initialcondTableStyle,fileContent)
 
   integer, dimension(:), allocatable, intent(out) :: &
-    microstructureAt, &
-    homogenizationAt
+    microstructureAt
   integer, intent(in) :: &
     nElem, &
     nNodes, &                                                                                       !< number of nodes per element
@@ -700,7 +697,6 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
 
 
   allocate(microstructureAt(nElem),source=0)
-  allocate(homogenizationAt(nElem),source=0)
 
   do l = 1, size(fileContent)
     chunkPos = IO_stringPos(fileContent(l))
@@ -720,7 +716,6 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
           do i = 1,contInts(1)
             e = mesh_FEM2DAMASK_elem(contInts(1+i))
             if (sv == 2) microstructureAt(e) = myVal
-            if (sv == 3) homogenizationAt(e) = myVal
           enddo
           if (initialcondTableStyle == 0) m = m + 1
         enddo
@@ -728,7 +723,7 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
     endif
   enddo
 
-end subroutine inputRead_microstructureAndHomogenization
+end subroutine inputRead_microstructure
 
 
 !--------------------------------------------------------------------------------------------------
diff --git a/src/material.f90 b/src/material.f90
index efdeeafcb..864760c76 100644
--- a/src/material.f90
+++ b/src/material.f90
@@ -80,7 +80,7 @@ module material
     damage_initialPhi                                                                               !< initial damage per each homogenization
 
   integer, dimension(:),     allocatable, public, protected :: &                                    ! (elem)
-    material_homogenizationAt                                                                       !< homogenization ID of each element (copy of discretization_homogenizationAt)
+    material_homogenizationAt                                                                       !< homogenization ID of each element
   integer, dimension(:,:),   allocatable, public, target :: &                                       ! (ip,elem) ToDo: ugly target for mapping hack
     material_homogenizationMemberAt                                                                 !< position of the element within its homogenization instance
   integer, dimension(:,:), allocatable, public, protected :: &                                      ! (constituent,elem)
diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90
index 7964e1220..ac6ccf528 100644
--- a/src/mesh/discretization_mesh.f90
+++ b/src/mesh/discretization_mesh.f90
@@ -77,7 +77,6 @@ subroutine discretization_mesh_init(restart)
   IS :: faceSetIS
   PetscErrorCode :: ierr
   integer, dimension(:), allocatable :: &
-    homogenizationAt, &
     microstructureAt
   class(tNode), pointer :: &
     num_mesh
@@ -165,7 +164,6 @@ subroutine discretization_mesh_init(restart)
   call mesh_FEM_build_ipVolumes(dimPlex)
 
   allocate(microstructureAt(mesh_NcpElems))
-  allocate(homogenizationAt(mesh_NcpElems),source=1)
   do j = 1, mesh_NcpElems
     call DMGetLabelValue(geomMesh,'material',j-1,microstructureAt(j),ierr)
     CHKERRQ(ierr)
@@ -179,7 +177,7 @@ subroutine discretization_mesh_init(restart)
 
   allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
 
-  call discretization_init(microstructureAt,homogenizationAt,&
+  call discretization_init(microstructureAt,&
                            reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), &
                            mesh_node0)
 

From 80b84cf76a2bb9bc75bcb8f99eafca06197ad7fe Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 29 Aug 2020 23:37:59 +0200
Subject: [PATCH 37/41] fixed hickup of state variables

two times the wrong location resulted in the correct behavior before
---
 src/marc/discretization_marc.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90
index 618e3195b..2930e02b1 100644
--- a/src/marc/discretization_marc.f90
+++ b/src/marc/discretization_marc.f90
@@ -715,7 +715,7 @@ subroutine inputRead_microstructure(microstructureAt,&
           contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements
           do i = 1,contInts(1)
             e = mesh_FEM2DAMASK_elem(contInts(1+i))
-            if (sv == 2) microstructureAt(e) = myVal
+            if (sv == 3) microstructureAt(e) = myVal
           enddo
           if (initialcondTableStyle == 0) m = m + 1
         enddo

From d66cdc8324298abc0f3613a18a71a2d42171c788 Mon Sep 17 00:00:00 2001
From: Martin Diehl <m.diehl@mpie.de>
Date: Sun, 30 Aug 2020 06:52:15 +0200
Subject: [PATCH 38/41] fixed incomplete rename

---
 src/marc/discretization_marc.f90 | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90
index 2930e02b1..fd0852257 100644
--- a/src/marc/discretization_marc.f90
+++ b/src/marc/discretization_marc.f90
@@ -226,9 +226,9 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt)
 
   connectivity_elem = inputRead_connectivityElem(nElems,elem%nNodes,inputFile)
 
-  call inputRead_microstructureAnd(microstructureAt, &
-                                   nElems,elem%nNodes,nameElemSet,mapElemSet,&
-                                   initialcondTableStyle,inputFile)
+  call inputRead_microstructure(microstructureAt, &
+                                nElems,elem%nNodes,nameElemSet,mapElemSet,&
+                                initialcondTableStyle,inputFile)
 end subroutine inputRead
 
 

From ffd45aa5fcf5b6f5ea11e3c8998cf3e8ffcad440 Mon Sep 17 00:00:00 2001
From: Test User <damask_user@mpie.de>
Date: Sun, 30 Aug 2020 20:54:25 +0200
Subject: [PATCH 39/41] [skip ci] updated version information after successful
 test of v3.0.0-alpha-47-g3ee8c471

---
 VERSION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/VERSION b/VERSION
index 73dd90e38..bcb442f66 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-v3.0.0-alpha-41-g94574356
+v3.0.0-alpha-47-g3ee8c471

From 33527c126315b6f731bdcf7b597a1ea2ed65005a Mon Sep 17 00:00:00 2001
From: Test User <damask_user@mpie.de>
Date: Mon, 31 Aug 2020 15:00:07 +0200
Subject: [PATCH 40/41] [skip ci] updated version information after successful
 test of v3.0.0-alpha-51-g31282973

---
 VERSION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/VERSION b/VERSION
index bcb442f66..70c328358 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-v3.0.0-alpha-47-g3ee8c471
+v3.0.0-alpha-51-g31282973

From 41f7c6597ac10a6667f9594d913377a4efbccc06 Mon Sep 17 00:00:00 2001
From: Test User <damask_user@mpie.de>
Date: Mon, 31 Aug 2020 20:03:24 +0200
Subject: [PATCH 41/41] [skip ci] updated version information after successful
 test of v3.0.0-alpha-75-g6986d20b

---
 VERSION | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/VERSION b/VERSION
index bcb442f66..812818d60 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-v3.0.0-alpha-47-g3ee8c471
+v3.0.0-alpha-75-g6986d20b