diff --git a/PRIVATE b/PRIVATE
index d96bfb329..3a2f89547 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit d96bfb32920c96a8a43958f76a209d34c6bd841a
+Subproject commit 3a2f89547c264044a7bfab9d33aee78eec495a76
diff --git a/VERSION b/VERSION
index 2d8c83361..d961b0f73 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-v2.0.3-304-g7b14263c
+v2.0.3-332-g5abcca50
diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py
index caf09d536..fa47805bb 100755
--- a/processing/post/DADF5_postResults.py
+++ b/processing/post/DADF5_postResults.py
@@ -22,10 +22,15 @@ parser.add_argument('filenames', nargs='+',
                     help='DADF5 files')
 parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string',
                     help='name of subdirectory to hold output')
+parser.add_argument('--mat', nargs='+',
+                    help='labels for materialpoint/homogenization',dest='mat')
+parser.add_argument('--con', nargs='+',
+                    help='labels for constituent/crystallite/constitutive',dest='con')
 
 options = parser.parse_args()
 
-options.labels = ['Fe','Fp','xi_sl']
+if options.mat is None: options.mat=[]
+if options.con is None: options.con=[]
 
 # --- loop over input files ------------------------------------------------------------------------
 
@@ -48,16 +53,13 @@ for filename in options.filenames:
     
     data = np.array([inc['inc'] for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1])
     header+= 'inc'
-    
-    data = np.concatenate((data,np.array([j+1 for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1])),1)
-    header+=' node'
 
     coords = coords.reshape([np.product(results.grid),3])
     data = np.concatenate((data,coords),1)
     header+=' 1_pos 2_pos 3_pos'
         
     results.active['increments'] = [inc]
-    for label in options.labels:
+    for label in options.con:
       for o in results.c_output_types:
         results.active['c_output_types'] = [o]
         for c in results.constituents:
@@ -67,12 +69,33 @@ for filename in options.filenames:
             continue
           label = x[0].split('/')[-1]
           array = results.read_dataset(x,0)
-          d = np.product(np.shape(array)[1:])
+          d = int(np.product(np.shape(array)[1:]))
           array = np.reshape(array,[np.product(results.grid),d])
           data = np.concatenate((data,array),1)
           
-          header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)])
-
+          if d>1:
+            header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)])
+          else:
+            header+=' '+label
+            
+    for label in options.mat:
+      for o in results.m_output_types:
+        results.active['m_output_types'] = [o]
+        for m in results.materialpoints:
+          results.active['materialpoints'] = [m]
+          x = results.get_dataset_location(label)
+          if len(x) == 0:
+            continue
+          label = x[0].split('/')[-1]
+          array = results.read_dataset(x,0)
+          d = int(np.product(np.shape(array)[1:]))
+          array = np.reshape(array,[np.product(results.grid),d])
+          data = np.concatenate((data,array),1)
+          
+          if d>1:
+            header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)])
+          else:
+            header+=' '+label
 
     dirname  = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
     try:
diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py
index dc1177488..ef27a973c 100755
--- a/processing/post/DADF5_vtk_cells.py
+++ b/processing/post/DADF5_vtk_cells.py
@@ -23,10 +23,15 @@ parser.add_argument('filenames', nargs='+',
                     help='DADF5 files')
 parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string',
                     help='name of subdirectory to hold output')
+parser.add_argument('--mat', nargs='+',
+                    help='labels for materialpoint/homogenization',dest='mat')
+parser.add_argument('--con', nargs='+',
+                    help='labels for constituent/crystallite/constitutive',dest='con')
 
 options = parser.parse_args()
 
-options.labels = ['Fe','Fp','xi_sl']
+if options.mat is None: options.mat=[]
+if options.con is None: options.con=[]
 
 # --- loop over input files ------------------------------------------------------------------------
 
@@ -54,7 +59,9 @@ for filename in options.filenames:
     print('Output step {}/{}'.format(i+1,len(results.increments)))
     vtk_data = []
     results.active['increments'] = [inc]
-    for label in options.labels:
+    
+    for label in options.con:
+      
       for o in results.c_output_types:
         results.active['c_output_types'] = [o]
         if o != 'generic':
diff --git a/processing/pre/geom_addPrimitive.py b/processing/pre/geom_addPrimitive.py
index 7fcfdbc5c..0dfd06732 100755
--- a/processing/pre/geom_addPrimitive.py
+++ b/processing/pre/geom_addPrimitive.py
@@ -43,7 +43,7 @@ parser.add_option('-e', '--exponent',  dest='exponent',
                   1 gives a sphere (|x|^(2^1) + |y|^(2^1) + |z|^(2^1) < 1), \
                   large values produce boxes, negative turns concave.')
 parser.add_option('-f', '--fill',       dest='fill',
-                  type='int', metavar = 'int',
+                  type='float', metavar = 'float',
                   help='grain index to fill primitive. "0" selects maximum microstructure index + 1 [%default]')
 parser.add_option('-q', '--quaternion', dest='quaternion',
                   type='float', nargs = 4, metavar=' '.join(['float']*4),
@@ -60,15 +60,24 @@ parser.add_option(     '--nonperiodic', dest='periodic',
 parser.add_option(     '--realspace',  dest='realspace',
                   action='store_true',
                   help = '-c and -d span [origin,origin+size] instead of [0,grid] coordinates')
+parser.add_option(     '--invert',  dest='inside',
+                  action='store_false',
+                  help = 'invert the volume filled by the primitive (inside/outside)')
+parser.add_option('--float',        dest = 'float',
+                  action = 'store_true',
+                  help = 'use float input')
 parser.set_defaults(center = (.0,.0,.0),
-                    fill = 0,
+                    fill = 0.0,
                     degrees = False,
                     exponent = (20,20,20), # box shape by default
                     periodic = True,
                     realspace = False,
+                    inside = True,
+                    float = False,
                    )
 
 (options, filenames) = parser.parse_args()
+
 if options.dimension is None:
   parser.error('no dimension specified.')   
 if options.angleaxis is not None:
@@ -78,6 +87,8 @@ elif options.quaternion is not None:
 else:
   rotation = damask.Rotation()
 
+datatype = 'f' if options.float else 'i'
+
 options.center = np.array(options.center)
 options.dimension = np.array(options.dimension)
 # undo logarithmic sense of exponent and generate ellipsoids for negative dimensions (backward compatibility)
@@ -97,13 +108,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
-                     'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
-                     'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
-                     'homogenization:  %i'%info['homogenization'],
-                     'microstructures: %i'%info['microstructures'],
-              ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
@@ -115,7 +120,7 @@ for name in filenames:
 
 #--- read data ------------------------------------------------------------------------------------
 
-  microstructure = table.microstructure_read(info['grid'])                                          # read microstructure
+  microstructure = table.microstructure_read(info['grid'],datatype)                          # read microstructure
 
 # --- do work ------------------------------------------------------------------------------------
 
@@ -123,7 +128,7 @@ for name in filenames:
              'microstructures': 0,
             }
 
-  options.fill = microstructure.max()+1 if options.fill == 0 else options.fill
+  options.fill = np.nanmax(microstructure)+1 if options.fill == 0 else options.fill
   
   microstructure = microstructure.reshape(info['grid'],order='F')
   
@@ -193,19 +198,23 @@ for name in filenames:
                          grid[1] * j : grid[1] * (j+1),
                          grid[2] * k : grid[2] * (k+1)])**options.exponent[2] <= 1.0)
     
-    microstructure = np.where(inside, options.fill, microstructure)   
+    microstructure = np.where(inside,
+                              options.fill   if options.inside else microstructure,
+                              microstructure if options.inside else options.fill)   
 
   else: # nonperiodic, much lighter on resources
     microstructure = np.where(np.abs(X)**options.exponent[0] +
                               np.abs(Y)**options.exponent[1] +
-                              np.abs(Z)**options.exponent[2] <= 1.0, options.fill, microstructure)
+                              np.abs(Z)**options.exponent[2] <= 1.0,
+                              options.fill   if options.inside else microstructure,
+                              microstructure if options.inside else options.fill)   
 
   np.seterr(**old_settings) # Reset warnings to old state
-  newInfo['microstructures'] = microstructure.max()
+  newInfo['microstructures'] = len(np.unique(microstructure))
 
 # --- report ---------------------------------------------------------------------------------------
   if (newInfo['microstructures'] != info['microstructures']):
-    damask.util.croak('--> microstructures: %i'%newInfo['microstructures'])
+    damask.util.croak('--> microstructures: {}'.format(newInfo['microstructures']))
 
 
 #--- write header ---------------------------------------------------------------------------------
@@ -225,9 +234,9 @@ for name in filenames:
 
 # --- write microstructure information ------------------------------------------------------------
 
-  formatwidth = int(math.floor(math.log10(microstructure.max())+1))
+  format = '%g' if options.float else '%{}i'.format(int(math.floor(math.log10(np.nanmax(microstructure))+1)))
   table.data = microstructure.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose()
-  table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ')
+  table.data_writeArray(format,delimiter = ' ')
 
 #--- output finalization --------------------------------------------------------------------------
 
diff --git a/processing/pre/geom_canvas.py b/processing/pre/geom_canvas.py
index d7fd1614a..01682deb8 100755
--- a/processing/pre/geom_canvas.py
+++ b/processing/pre/geom_canvas.py
@@ -35,7 +35,7 @@ parser.add_option('-f',
                   type = 'float', metavar = 'float',
                   help = '(background) canvas grain index. "0" selects maximum microstructure index + 1 [%default]')
 parser.add_option('--float',
-                  dest = 'real',
+                  dest = 'float',
                   action = 'store_true',
                   help = 'use float input')
 parser.add_option('--blank',
@@ -45,13 +45,13 @@ parser.add_option('--blank',
 
 parser.set_defaults(grid = ['0','0','0'],
                     offset = (0,0,0),
-                    fill = 0,
-                    real = False,
+                    fill = 0.0,
+                    float = False,
                    )
 
 (options, filenames) = parser.parse_args()
 
-datatype = 'f' if options.real else 'i'
+datatype = 'f' if options.float else 'i'
 options.grid = ['1','1','1'] if options.blank and options.grid == ['0','0','0'] else options.grid
 options.fill = 1 if options.blank and options.fill == 0 else options.fill
 
@@ -107,7 +107,7 @@ for name in filenames:
   newInfo['grid'] = np.where(newInfo['grid'] > 0, newInfo['grid'],info['grid'])
 
   microstructure_cropped = np.zeros(newInfo['grid'],datatype)
-  microstructure_cropped.fill(options.fill if options.real or options.fill > 0 else microstructure.max()+1)
+  microstructure_cropped.fill(options.fill if options.float or options.fill > 0 else np.nanmax(microstructure)+1)
   
   if not options.blank:
     xindex = np.arange(max(options.offset[0],0),min(options.offset[0]+newInfo['grid'][0],info['grid'][0]))
@@ -130,7 +130,7 @@ for name in filenames:
   newInfo['size']   = info['size']/info['grid']*newInfo['grid'] if np.all(info['grid'] > 0) else newInfo['grid']
   newInfo['origin'] = info['origin']+(info['size']/info['grid'] if np.all(info['grid'] > 0) \
                               else newInfo['size']/newInfo['grid'])*options.offset
-  newInfo['microstructures'] = microstructure_cropped.max()
+  newInfo['microstructures'] = len(np.unique(microstructure_cropped))
 
 # --- report ---------------------------------------------------------------------------------------
 
@@ -172,7 +172,7 @@ for name in filenames:
 
 # --- write microstructure information ------------------------------------------------------------
 
-  format = '%g' if options.real else '%{}i'.format(int(math.floor(math.log10(microstructure_cropped.max())+1)))
+  format = '%g' if options.float else '%{}i'.format(int(math.floor(math.log10(np.nanmax(microstructure_cropped))+1)))
   table.data = microstructure_cropped.reshape((newInfo['grid'][0],newInfo['grid'][1]*newInfo['grid'][2]),order='F').transpose()
   table.data_writeArray(format,delimiter=' ')
     
diff --git a/processing/pre/geom_clean.py b/processing/pre/geom_clean.py
index 907431146..1d0769ab3 100755
--- a/processing/pre/geom_clean.py
+++ b/processing/pre/geom_clean.py
@@ -50,13 +50,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-
-  damask.util.croak(['grid     a b c:  {}'.format(' x '.join(map(str,info['grid']))),
-                     'size     x y z:  {}'.format(' x '.join(map(str,info['size']))),
-                     'origin   x y z:  {}'.format(' : '.join(map(str,info['origin']))),
-                     'homogenization:  {}'.format(info['homogenization']),
-                     'microstructures: {}'.format(info['microstructures']),
-                    ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
@@ -73,7 +67,7 @@ for name in filenames:
 # --- do work ------------------------------------------------------------------------------------
 
   microstructure = ndimage.filters.generic_filter(microstructure,mostFrequent,size=(options.stencil,)*3).astype('int_')
-  newInfo = {'microstructures': microstructure.max()}
+  newInfo = {'microstructures': len(np.unique(microstructure))}
 
 # --- report ---------------------------------------------------------------------------------------
   if (    newInfo['microstructures'] != info['microstructures']):
@@ -91,9 +85,9 @@ for name in filenames:
 
 # --- write microstructure information ------------------------------------------------------------
 
-  formatwidth = int(math.floor(math.log10(microstructure.max())+1))
+  formatwidth = int(math.floor(math.log10(np.nanmax(microstructure))+1))
   table.data = microstructure.reshape((info['grid'][0],np.prod(info['grid'][1:])),order='F').transpose()
-  table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ')
+  table.data_writeArray('%{}i'.format(formatwidth),delimiter = ' ')
 
 # --- output finalization --------------------------------------------------------------------------
 
diff --git a/processing/pre/geom_fromMinimalSurface.py b/processing/pre/geom_fromMinimalSurface.py
index 002b4800b..e0023e7ec 100755
--- a/processing/pre/geom_fromMinimalSurface.py
+++ b/processing/pre/geom_fromMinimalSurface.py
@@ -90,12 +90,7 @@ for name in filenames:
 
 #--- report ---------------------------------------------------------------------------------------
 
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
-               'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
-               'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
-               'homogenization:  %i'%info['homogenization'],
-               'microstructures: %i'%info['microstructures'],
-              ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py
index 8eb1ed8bf..7a905cd26 100755
--- a/processing/pre/geom_fromTable.py
+++ b/processing/pre/geom_fromTable.py
@@ -192,12 +192,7 @@ for name in filenames:
           'homogenization':  options.homogenization,
          }
 
-  damask.util.croak(['grid     a b c:  {}'.format(' x '.join(map(str,info['grid']))),
-                     'size     x y z:  {}'.format(' x '.join(map(str,info['size']))),
-                     'origin   x y z:  {}'.format(' : '.join(map(str,info['origin']))),
-                     'homogenization:  {}'.format(info['homogenization']),
-                     'microstructures: {}'.format(info['microstructures']),
-                    ])
+  damask.util.report_geom(info)
     
 # --- write header ---------------------------------------------------------------------------------
 
@@ -230,7 +225,7 @@ for name in filenames:
 # --- write microstructure information ------------------------------------------------------------
 
   table.data = grain.reshape(info['grid'][1]*info['grid'][2],info['grid'][0])
-  table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ')
+  table.data_writeArray('%{}i'.format(formatwidth),delimiter=' ')
   
 #--- output finalization --------------------------------------------------------------------------
 
diff --git a/processing/pre/geom_grainGrowth.py b/processing/pre/geom_grainGrowth.py
index 1afb02715..f7c50c2e5 100755
--- a/processing/pre/geom_grainGrowth.py
+++ b/processing/pre/geom_grainGrowth.py
@@ -69,13 +69,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-  
-  damask.util.croak(['grid     a b c:  {}'.format(' x '.join(list(map(str,info['grid'])))),
-                     'size     x y z:  {}'.format(' x '.join(list(map(str,info['size'])))),
-                     'origin   x y z:  {}'.format(' : '.join(list(map(str,info['origin'])))),
-                     'homogenization:  {}'.format(info['homogenization']),
-                     'microstructures: {}'.format(info['microstructures']),
-                    ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
@@ -200,8 +194,7 @@ for name in filenames:
         newID += 1
         microstructure = np.where(microstructure == microstructureID, newID, microstructure)
 
-  newInfo = {'microstructures': 0,}
-  newInfo['microstructures'] = microstructure.max()
+  newInfo = {'microstructures': len(np.unique(microstructure)),}
 
 # --- report --------------------------------------------------------------------------------------
 
@@ -226,7 +219,7 @@ for name in filenames:
 
 # --- write microstructure information ------------------------------------------------------------
 
-  formatwidth = int(math.floor(math.log10(microstructure.max())+1))
+  formatwidth = int(math.floor(math.log10(np.nanmax(microstructure))+1))
   table.data = microstructure[::1 if info['grid'][0]>1 else 2,
                               ::1 if info['grid'][1]>1 else 2,
                               ::1 if info['grid'][2]>1 else 2,].\
diff --git a/processing/pre/geom_mirror.py b/processing/pre/geom_mirror.py
index 951fb0842..853b99632 100755
--- a/processing/pre/geom_mirror.py
+++ b/processing/pre/geom_mirror.py
@@ -23,6 +23,13 @@ parser.add_option('-d','--direction',
                   dest = 'directions',
                   action = 'extend', metavar = '<string LIST>',
                   help = "directions in which to mirror {'x','y','z'}")
+parser.add_option('--float',
+                  dest = 'float',
+                  action = 'store_true',
+                  help = 'use float input')
+
+parser.set_defaults(float = False,
+                   )
 
 (options, filenames) = parser.parse_args()
 
@@ -32,6 +39,8 @@ if not set(options.directions).issubset(validDirections):
   invalidDirections = [str(e) for e in set(options.directions).difference(validDirections)]
   parser.error('invalid directions {}. '.format(*invalidDirections))
 
+datatype = 'f' if options.float else 'i'
+
 # --- loop over input files -------------------------------------------------------------------------
 
 if filenames == []: filenames = [None]
@@ -39,7 +48,8 @@ if filenames == []: filenames = [None]
 for name in filenames:
   try:
     table = damask.ASCIItable(name = name,
-                              buffered = False, labeled = False)
+                              buffered = False,
+                              labeled = False)
   except: continue
   damask.util.report(scriptName,name)
 
@@ -47,13 +57,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
-                     'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
-                     'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
-                     'homogenization:  %i'%info['homogenization'],
-                     'microstructures: %i'%info['microstructures'],
-                    ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
@@ -65,7 +69,7 @@ for name in filenames:
 
 # --- read data ------------------------------------------------------------------------------------
 
-  microstructure = table.microstructure_read(info['grid']).reshape(info['grid'],order='F')          # read microstructure
+  microstructure = table.microstructure_read(info['grid'],datatype).reshape(info['grid'],order='F')  # read microstructure
 
   if 'z' in options.directions:
     microstructure = np.concatenate([microstructure,microstructure[:,:,::-1]],2)
@@ -107,9 +111,9 @@ for name in filenames:
 
 # --- write microstructure information ------------------------------------------------------------
 
-  formatwidth = int(math.floor(math.log10(microstructure.max())+1))
+  formatwidth = int(math.floor(math.log10(np.nanmax(microstructure))+1))
   table.data = microstructure.reshape((newInfo['grid'][0],np.prod(newInfo['grid'][1:])),order='F').transpose()
-  table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ')
+  table.data_writeArray('%{}i'.format(formatwidth),delimiter = ' ')
 
 # --- output finalization --------------------------------------------------------------------------
 
diff --git a/processing/pre/geom_pack.py b/processing/pre/geom_pack.py
index 0d864bbf5..2e6080a6b 100755
--- a/processing/pre/geom_pack.py
+++ b/processing/pre/geom_pack.py
@@ -35,14 +35,8 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-  
-  damask.util.croak(['grid     a b c:  {}'.format(' x '.join(map(str,info['grid']))),
-                     'size     x y z:  {}'.format(' x '.join(map(str,info['size']))),
-                     'origin   x y z:  {}'.format(' : '.join(map(str,info['origin']))),
-                     'homogenization:  {}'.format(info['homogenization']),
-                     'microstructures: {}'.format(info['microstructures']),
-                    ])
-
+  damask.util.report_geom(info)
+ 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
   if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.')
diff --git a/processing/pre/geom_renumber.py b/processing/pre/geom_renumber.py
index 033b4a566..3faa7f449 100755
--- a/processing/pre/geom_renumber.py
+++ b/processing/pre/geom_renumber.py
@@ -35,13 +35,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
-               'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
-               'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
-               'homogenization:  %i'%info['homogenization'],
-               'microstructures: %i'%info['microstructures'],
-              ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
@@ -93,7 +87,7 @@ for name in filenames:
 
 # --- write microstructure information -----------------------------------------------------------
 
-  format = '%{}i'.format(int(math.floor(math.log10(newInfo['microstructures'])+1)))
+  format = '%{}i'.format(int(math.floor(math.log10(np.nanmax(renumbered))+1)))
   table.data = renumbered.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose()
   table.data_writeArray(format,delimiter = ' ')
 
diff --git a/processing/pre/geom_rescale.py b/processing/pre/geom_rescale.py
index b3716bd62..4a14c0050 100755
--- a/processing/pre/geom_rescale.py
+++ b/processing/pre/geom_rescale.py
@@ -31,14 +31,21 @@ parser.add_option('-r', '--renumber',
                   dest = 'renumber',
                   action = 'store_true',
                   help = 'renumber microstructure indices from 1..N [%default]')
+parser.add_option('--float',
+                  dest = 'float',
+                  action = 'store_true',
+                  help = 'use float input')
 
 parser.set_defaults(renumber = False,
                     grid = ['0','0','0'],
                     size  = ['0.0','0.0','0.0'],
+                    float = False,
                    )
 
 (options, filenames) = parser.parse_args()
 
+datatype = 'f' if options.float else 'i'
+
 # --- loop over input files -------------------------------------------------------------------------
 
 if filenames == []: filenames = [None]
@@ -46,7 +53,8 @@ if filenames == []: filenames = [None]
 for name in filenames:
   try:
     table = damask.ASCIItable(name = name,
-                              buffered = False, labeled = False)
+                              buffered = False,
+                              labeled = False)
   except: continue
   damask.util.report(scriptName,name)
 
@@ -54,13 +62,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-  
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
-               'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
-               'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
-               'homogenization:  %i'%info['homogenization'],
-               'microstructures: %i'%info['microstructures'],
-              ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
@@ -72,7 +74,7 @@ for name in filenames:
 
 # --- read data ------------------------------------------------------------------------------------
 
-  microstructure = table.microstructure_read(info['grid'])                                          # read microstructure
+  microstructure = table.microstructure_read(info['grid'],datatype)                                   # read microstructure
 
 # --- do work ------------------------------------------------------------------------------------
  
@@ -113,7 +115,7 @@ for name in filenames:
         newID += 1
         microstructure = np.where(microstructure == microstructureID, newID,microstructure).reshape(microstructure.shape)
 
-  newInfo['microstructures'] = microstructure.max()
+  newInfo['microstructures'] = len(np.unique(microstructure))
 
 # --- report ---------------------------------------------------------------------------------------
 
@@ -152,9 +154,9 @@ for name in filenames:
   
 # --- write microstructure information ------------------------------------------------------------
 
-  formatwidth = int(math.floor(math.log10(microstructure.max())+1))
+  format = '%g' if options.float else '%{}i'.format(int(math.floor(math.log10(np.nanmax(microstructure))+1)))
   table.data = microstructure.reshape((newInfo['grid'][0],newInfo['grid'][1]*newInfo['grid'][2]),order='F').transpose()
-  table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ')
+  table.data_writeArray(format,delimiter=' ')
     
 # --- output finalization --------------------------------------------------------------------------
 
diff --git a/processing/pre/geom_rotate.py b/processing/pre/geom_rotate.py
index 4da59cddf..7cce5800d 100755
--- a/processing/pre/geom_rotate.py
+++ b/processing/pre/geom_rotate.py
@@ -43,9 +43,15 @@ parser.add_option('-f', '--fill',
                   dest = 'fill',
                   type = 'int', metavar = 'int',
                   help = 'background grain index. "0" selects maximum microstructure index + 1 [%default]')
+parser.add_option('--float',
+                  dest = 'float',
+                  action = 'store_true',
+                  help = 'use float input')
 
 parser.set_defaults(degrees = False,
-                    fill = 0)
+                    fill = 0,
+                    float = False,
+                   )
 
 (options, filenames) = parser.parse_args()
 
@@ -61,6 +67,8 @@ if options.matrix is not None:
 if options.eulers is not None:
   eulers = damask.Rotation.fromEulers(np.array(options.eulers),degrees=True).asEulers(degrees=True)
 
+datatype = 'f' if options.float else 'i'
+
 # --- loop over input files -------------------------------------------------------------------------
 
 if filenames == []: filenames = [None]
@@ -77,13 +85,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-
-  damask.util.croak(['grid     a b c:  {}'.format(' x '.join(map(str,info['grid']))),
-                     'size     x y z:  {}'.format(' x '.join(map(str,info['size']))),
-                     'origin   x y z:  {}'.format(' : '.join(map(str,info['origin']))),
-                     'homogenization:  {}'.format(info['homogenization']),
-                     'microstructures: {}'.format(info['microstructures']),
-                    ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
@@ -95,9 +97,9 @@ for name in filenames:
 
 # --- read data ------------------------------------------------------------------------------------
 
-  microstructure = table.microstructure_read(info['grid']).reshape(info['grid'],order='F')                    # read microstructure
+  microstructure = table.microstructure_read(info['grid'],datatype).reshape(info['grid'],order='F')                  # read microstructure
 
-  newGrainID = options.fill if options.fill != 0 else microstructure.max()+1
+  newGrainID = options.fill if options.fill != 0 else np.nanmax(microstructure)+1
   microstructure = ndimage.rotate(microstructure,eulers[2],(0,1),order=0,prefilter=False,output=int,cval=newGrainID) # rotation around Z
   microstructure = ndimage.rotate(microstructure,eulers[1],(1,2),order=0,prefilter=False,output=int,cval=newGrainID) # rotation around X
   microstructure = ndimage.rotate(microstructure,eulers[0],(0,1),order=0,prefilter=False,output=int,cval=newGrainID) # rotation around Z
@@ -107,19 +109,18 @@ for name in filenames:
   newInfo = {
              'size':   microstructure.shape*info['size']/info['grid'],
              'grid':   microstructure.shape,
-             'microstructures': microstructure.max(),
+             'microstructures': len(np.unique(microstructure)),
             }
 
-
 # --- report ---------------------------------------------------------------------------------------
 
   remarks = []
   if (any(newInfo['grid']            != info['grid'])):
-    remarks.append('--> grid    a b c:  %s'%(' x '.join(map(str,newInfo['grid']))))
+    remarks.append('--> grid    a b c:  {}'.format(' x '.join(map(str,newInfo['grid']))))
   if (any(newInfo['size']            != info['size'])):
-    remarks.append('--> size    x y z:  %s'%(' x '.join(map(str,newInfo['size']))))
+    remarks.append('--> size    x y z:  {}'.format(' x '.join(map(str,newInfo['size']))))
   if (    newInfo['microstructures'] != info['microstructures']): 
-    remarks.append('--> microstructures: %i'%newInfo['microstructures'])
+    remarks.append('--> microstructures: {}'.format(newInfo['microstructures']))
   if remarks != []: damask.util.croak(remarks)
 
 # --- write header ---------------------------------------------------------------------------------
@@ -138,9 +139,9 @@ for name in filenames:
 
 # --- write microstructure information ------------------------------------------------------------
 
-  formatwidth = int(math.floor(math.log10(microstructure.max())+1))
+  format = '%g' if options.float else '%{}i'.format(int(math.floor(math.log10(np.nanmax(microstructure))+1)))
   table.data = microstructure.reshape((newInfo['grid'][0],np.prod(newInfo['grid'][1:])),order='F').transpose()
-  table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ')
+  table.data_writeArray(format,delimiter=' ')
 
 # --- output finalization --------------------------------------------------------------------------
 
diff --git a/processing/pre/geom_toTable.py b/processing/pre/geom_toTable.py
index 73e4888d1..0a71b335e 100755
--- a/processing/pre/geom_toTable.py
+++ b/processing/pre/geom_toTable.py
@@ -20,15 +20,15 @@ Translate geom description into ASCIItable containing position and microstructur
 """, version = scriptID)
 
 parser.add_option('--float',
-                  dest = 'real',
+                  dest = 'float',
                   action = 'store_true',
                   help = 'use float input')
 
-parser.set_defaults(real = False,
+parser.set_defaults(float = False,
                    )
 (options, filenames) = parser.parse_args()
 
-datatype = 'f' if options.real else 'i'
+datatype = 'f' if options.float else 'i'
 
 # --- loop over input files -------------------------------------------------------------------------
 
@@ -47,13 +47,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-  
-  damask.util.croak(['grid     a b c:  {}'.format(' x '.join(list(map(str,info['grid'])))),
-                     'size     x y z:  {}'.format(' x '.join(list(map(str,info['size'])))),
-                     'origin   x y z:  {}'.format(' : '.join(list(map(str,info['origin'])))),
-                     'homogenization:  {}'.format(info['homogenization']),
-                     'microstructures: {}'.format(info['microstructures']),
-                    ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
diff --git a/processing/pre/geom_translate.py b/processing/pre/geom_translate.py
index 59aaac5d5..072c270ea 100755
--- a/processing/pre/geom_translate.py
+++ b/processing/pre/geom_translate.py
@@ -31,19 +31,19 @@ parser.add_option('-s', '--substitute',
                   action = 'extend', metavar = '<string LIST>',
                   help = 'substitutions of microstructure indices from,to,from,to,...')
 parser.add_option('--float',
-                  dest = 'real',
+                  dest = 'float',
                   action = 'store_true',
                   help = 'use float input')
 
 parser.set_defaults(origin = (0.0,0.0,0.0),
                     microstructure = 0,
                     substitute = [],
-                    real = False,
+                    float = False,
                    )
 
 (options, filenames) = parser.parse_args()
 
-datatype = 'f' if options.real else 'i'
+datatype = 'f' if options.float else 'i'
 
 sub = {}
 for i in range(len(options.substitute)//2):                                                         # split substitution list into "from" -> "to"
@@ -64,13 +64,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
-               'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
-               'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
-               'homogenization:  %i'%info['homogenization'],
-               'microstructures: %i'%info['microstructures'],
-              ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
@@ -92,7 +86,7 @@ for name in filenames:
             }
 
   substituted = np.copy(microstructure)
-  for k, v in sub.items(): substituted[microstructure==k] = v                                   # substitute microstructure indices
+  for k, v in sub.items(): substituted[microstructure==k] = v                                       # substitute microstructure indices
 
   substituted += options.microstructure                                                             # shift microstructure indices
 
@@ -103,9 +97,9 @@ for name in filenames:
 
   remarks = []
   if (any(newInfo['origin']          != info['origin'])):
-    remarks.append('--> origin   x y z:  %s'%(' : '.join(map(str,newInfo['origin']))))
+    remarks.append('--> origin   x y z:  {}'.format(' : '.join(map(str,newInfo['origin']))))
   if (    newInfo['microstructures'] != info['microstructures']):
-    remarks.append('--> microstructures: %i'%newInfo['microstructures'])
+    remarks.append('--> microstructures: {}'.format(newInfo['microstructures']))
   if remarks != []: damask.util.croak(remarks)
 
 # --- write header -------------------------------------------------------------------------------
@@ -124,7 +118,7 @@ for name in filenames:
 
 # --- write microstructure information -----------------------------------------------------------
 
-  format = '%g' if options.real else '%{}i'.format(int(math.floor(math.log10(microstructure.max())+1)))
+  format = '%g' if options.float else '%{}i'.format(int(math.floor(math.log10(np.nanmax(substituted))+1)))
   table.data = substituted.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose()
   table.data_writeArray(format,delimiter = ' ')
 
diff --git a/processing/pre/geom_unpack.py b/processing/pre/geom_unpack.py
index 726e4ef04..4cac76c5f 100755
--- a/processing/pre/geom_unpack.py
+++ b/processing/pre/geom_unpack.py
@@ -43,13 +43,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-  
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
-               'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
-               'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
-               'homogenization:  %i'%info['homogenization'],
-               'microstructures: %i'%info['microstructures'],
-              ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
diff --git a/processing/pre/geom_vicinityOffset.py b/processing/pre/geom_vicinityOffset.py
index 9fce7201a..733276d01 100755
--- a/processing/pre/geom_vicinityOffset.py
+++ b/processing/pre/geom_vicinityOffset.py
@@ -73,13 +73,7 @@ for name in filenames:
 
   table.head_read()
   info,extra_header = table.head_getGeom()
-  
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
-               'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
-               'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
-               'homogenization:  %i'%info['homogenization'],
-               'microstructures: %i'%info['microstructures'],
-              ])
+  damask.util.report_geom(info)
 
   errors = []
   if np.any(info['grid'] < 1):    errors.append('invalid grid a b c.')
@@ -108,7 +102,7 @@ for name in filenames:
                                                            extra_keywords={"trigger":options.trigger,"size":1+2*options.vicinity}),
                             microstructure + options.offset,microstructure)
 
-  newInfo['microstructures'] = microstructure.max()
+  newInfo['microstructures'] = len(np.unique(microstructure))
 
 # --- report ---------------------------------------------------------------------------------------
 
@@ -131,9 +125,9 @@ for name in filenames:
   
 # --- write microstructure information ------------------------------------------------------------
 
-  formatwidth = int(math.floor(math.log10(microstructure.max())+1))
+  formatwidth = int(math.floor(math.log10(np.nanmax(microstructure))+1))
   table.data = microstructure.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose()
-  table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ')
+  table.data_writeArray('%{}i'.format(formatwidth),delimiter = ' ')
     
 # --- output finalization --------------------------------------------------------------------------
 
diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py
index 841029af8..f5fc81b16 100644
--- a/python/damask/dadf5.py
+++ b/python/damask/dadf5.py
@@ -48,16 +48,24 @@ class DADF5():
         for o in f['inc{:05}/constituent/{}'.format(self.increments[0]['inc'],c)].keys():
           self.c_output_types.append(o)
       self.c_output_types = list(set(self.c_output_types))                                          # make unique
+
+      self.m_output_types = []
+      for m in self.materialpoints:
+        for o in f['inc{:05}/materialpoint/{}'.format(self.increments[0]['inc'],m)].keys():
+          self.m_output_types.append(o)
+      self.m_output_types = list(set(self.m_output_types))                                          # make unique
       
     self.active= {'increments':     self.increments,
                   'constituents':   self.constituents,
                   'materialpoints': self.materialpoints,
                   'constituent':    self.Nconstituents,
-                  'c_output_types': self.c_output_types}
+                  'c_output_types': self.c_output_types,
+                  'm_output_types': self.m_output_types}
 
     self.filename   = filename
     self.mode       = mode
-    
+
+
   def list_data(self):
     """Shows information on all datasets in the file"""
     with h5py.File(self.filename,'r') as f:
@@ -73,6 +81,16 @@ class DADF5():
               print('    {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode()))
           except:
             pass
+      for m in self.active['materialpoints']:
+        group_materialpoint = group_inc+'/materialpoint/'+m
+        for t in self.active['m_output_types']:
+          print('  {}'.format(t))
+          group_output_types = group_materialpoint+'/'+t
+          try:
+            for x in f[group_output_types].keys():
+              print('    {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode()))
+          except:
+            pass
     
 
   def get_dataset_location(self,label):
@@ -81,14 +99,25 @@ class DADF5():
     with h5py.File(self.filename,'r') as f:
       for i in self.active['increments']:
         group_inc = 'inc{:05}'.format(i['inc'])
+        
         for c in self.active['constituents']:
           group_constituent = group_inc+'/constituent/'+c
           for t in self.active['c_output_types']:
             try:
               f[group_constituent+'/'+t+'/'+label]
               path.append(group_constituent+'/'+t+'/'+label)
-            except:
-              pass
+            except Exception as e:
+              print('unable to locate constituents dataset: '+ str(e))
+       
+        for m in self.active['materialpoints']:
+          group_materialpoint = group_inc+'/materialpoint/'+m
+          for t in self.active['m_output_types']:
+            try:
+              f[group_materialpoint+'/'+t+'/'+label]
+              path.append(group_materialpoint+'/'+t+'/'+label)
+            except Exception as e:
+              print('unable to locate materialpoints dataset: '+ str(e))
+              
     return path
     
     
@@ -100,13 +129,29 @@ class DADF5():
     """
     with h5py.File(self.filename,'r') as f:
       shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:]
+      if len(shape) == 1: shape = shape +(1,)
       dataset = np.full(shape,np.nan)
       for pa in path:
         label   = pa.split('/')[2]
-        p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0]
-        u = (f['mapping/cellResults/constituent'][p,c]['Position'])
-        dataset[p,:] = f[pa][u,:]                                                    
-    
+        try:
+          p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0]
+          u = (f['mapping/cellResults/constituent'][p,c]['Position'])
+          a = np.array(f[pa])
+          if len(a.shape) == 1:
+            a=a.reshape([a.shape[0],1])
+          dataset[p,:] = a[u,:]
+        except Exception as e:
+          print('unable to read constituent: '+ str(e))
+        try:
+          p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0]
+          u = (f['mapping/cellResults/materialpoint'][p.tolist()]['Position'])
+          a = np.array(f[pa])
+          if len(a.shape) == 1:
+            a=a.reshape([a.shape[0],1])
+          dataset[p,:] = a[u,:]
+        except Exception as e:
+          print('unable to read materialpoint: '+ str(e))
+
     return dataset
         
       
diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90
index 0ae1323f0..5e1cd71eb 100644
--- a/src/commercialFEM_fileList.f90
+++ b/src/commercialFEM_fileList.f90
@@ -46,7 +46,6 @@
 #include "plastic_nonlocal.f90"
 #include "constitutive.f90"
 #include "crystallite.f90"
-#include "homogenization_mech_RGC.f90"
 #include "thermal_isothermal.f90"
 #include "thermal_adiabatic.f90"
 #include "thermal_conduction.f90"
@@ -56,4 +55,5 @@
 #include "homogenization.f90"
 #include "homogenization_mech_none.f90"
 #include "homogenization_mech_isostrain.f90"
+#include "homogenization_mech_RGC.f90"
 #include "CPFEM.f90"
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 605cf5094..3210f02d4 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -6,7 +6,20 @@
 !--------------------------------------------------------------------------------------------------
 module homogenization
  use prec
+ use IO
+ use config
+ use debug
+ use math
  use material
+ use numerics
+ use constitutive
+ use crystallite
+ use mesh
+ use FEsolving
+#if defined(PETSc) || defined(DAMASK_HDF5)
+ use results
+ use HDF5_utilities
+#endif
 
 !--------------------------------------------------------------------------------------------------
 ! General variables for the homogenization at a  material point
@@ -22,7 +35,6 @@ module homogenization
    materialpoint_results                                                                            !< results array of material point
  integer,                                            public, protected  :: &
    materialpoint_sizeResults, &
-   homogenization_maxSizePostResults, &
    thermal_maxSizePostResults, &
    damage_maxSizePostResults
 
@@ -47,11 +59,24 @@ module homogenization
    module subroutine mech_isostrain_init
    end subroutine mech_isostrain_init
    
+   module subroutine mech_RGC_init
+   end subroutine mech_RGC_init
+   
+   
    module subroutine mech_isostrain_partitionDeformation(F,avgF)
      real(pReal),   dimension (:,:,:), intent(out) :: F                                             !< partitioned deformation gradient
      real(pReal),   dimension (3,3),   intent(in)  :: avgF                                          !< average deformation gradient at material point
    end subroutine mech_isostrain_partitionDeformation
    
+   module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
+     real(pReal),   dimension (:,:,:), intent(out) :: F                                             !< partitioned deformation gradient
+     real(pReal),   dimension (3,3),   intent(in)  :: avgF                                          !< average deformation gradient at material point
+     integer,                          intent(in)  :: &
+       instance, &
+       of
+   end subroutine mech_RGC_partitionDeformation
+   
+   
    module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
      real(pReal),   dimension (3,3),       intent(out) :: avgP                                      !< average stress at material point
      real(pReal),   dimension (3,3,3,3),   intent(out) :: dAvgPdAvgF                                !< average stiffness at material point
@@ -60,7 +85,37 @@ module homogenization
      real(pReal),   dimension (:,:,:,:,:), intent(in)  :: dPdF                                      !< partitioned stiffnesses
      integer,                              intent(in)  :: instance 
    end subroutine mech_isostrain_averageStressAndItsTangent
+   
+   module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
+     real(pReal),   dimension (3,3),       intent(out) :: avgP                                      !< average stress at material point
+     real(pReal),   dimension (3,3,3,3),   intent(out) :: dAvgPdAvgF                                !< average stiffness at material point
+ 
+     real(pReal),   dimension (:,:,:),     intent(in)  :: P                                         !< partitioned stresses
+     real(pReal),   dimension (:,:,:,:,:), intent(in)  :: dPdF                                      !< partitioned stiffnesses
+     integer,                              intent(in)  :: instance 
+   end subroutine mech_RGC_averageStressAndItsTangent
+   
+   
+   module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
+     logical, dimension(2) :: mech_RGC_updateState
+     real(pReal), dimension(:,:,:),     intent(in)    :: & 
+       P,&                                                                                          !< partitioned stresses
+       F,&                                                                                          !< partitioned deformation gradients
+       F0                                                                                           !< partitioned initial deformation gradients
+    real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF                                           !< partitioned stiffnesses
+    real(pReal), dimension(3,3),       intent(in) :: avgF                                           !< average F
+    real(pReal),                       intent(in) :: dt                                             !< time increment
+    integer,                           intent(in) :: &
+     ip, &                                                                                          !< integration point number
+     el                                                                                             !< element number
+   end function mech_RGC_updateState
 
+
+   module subroutine mech_RGC_results(instance,group)
+     integer,          intent(in) :: instance                                                       !< homogenization instance
+     character(len=*), intent(in) :: group                                                          !< group name in HDF5 file
+   end subroutine mech_RGC_results
+  
  end interface
 
  public ::  &
@@ -68,11 +123,6 @@ module homogenization
    materialpoint_stressAndItsTangent, &
    materialpoint_postResults, &
    homogenization_results
- private :: &
-   partitionDeformation, &
-   updateState, &
-   averageStressAndItsTangent, &
-   postResults
 
 contains
 
@@ -81,36 +131,12 @@ contains
 !> @brief module initialization
 !--------------------------------------------------------------------------------------------------
 subroutine homogenization_init
- use math, only: &
-   math_I3
- use debug, only: &
-   debug_level, &
-   debug_homogenization, &
-   debug_levelBasic, &
-   debug_e, &
-   debug_g
- use mesh, only: &
-   theMesh, &
-   mesh_element
- use constitutive, only: &
-   constitutive_plasticity_maxSizePostResults, &
-   constitutive_source_maxSizePostResults
- use crystallite, only: &
-   crystallite_maxSizePostResults
- use config, only: &
-  config_deallocate, &
-  config_homogenization, &
-  homogenization_name
- use homogenization_mech_RGC
  use thermal_isothermal
  use thermal_adiabatic
  use thermal_conduction
  use damage_none
  use damage_local
  use damage_nonlocal
- use IO
- use numerics, only: &
-   worldrank
 
  integer, parameter :: FILEUNIT = 200
  integer :: e,i,p
@@ -120,10 +146,9 @@ subroutine homogenization_init
  character(len=32) :: outputName                                                                    !< name of output, intermediate fix until HDF5 output is ready
  logical :: valid
 
-
  if (any(homogenization_type == HOMOGENIZATION_NONE_ID))      call mech_none_init
  if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init
- if (any(homogenization_type == HOMOGENIZATION_RGC_ID))       call homogenization_RGC_init
+ if (any(homogenization_type == HOMOGENIZATION_RGC_ID))       call mech_RGC_init
 
  if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init
  if (any(thermal_type == THERMAL_adiabatic_ID))  call thermal_adiabatic_init
@@ -135,39 +160,14 @@ subroutine homogenization_init
 
 !--------------------------------------------------------------------------------------------------
 ! write description file for homogenization output
- mainProcess2: if (worldrank == 0) then
+ mainProcess: if (worldrank == 0) then
    call IO_write_jobFile(FILEUNIT,'outputHomogenization')
    do p = 1,size(config_homogenization)
      if (any(material_homogenizationAt == p)) then
-       i = homogenization_typeInstance(p)                                                               ! which instance of this homogenization type
-       valid = .true.                                                                                   ! assume valid
-       select case(homogenization_type(p))                                                              ! split per homogenization type
-         case (HOMOGENIZATION_NONE_ID)
-           outputName = HOMOGENIZATION_NONE_label
-           thisOutput => null()
-           thisSize   => null()
-         case (HOMOGENIZATION_ISOSTRAIN_ID)
-           outputName = HOMOGENIZATION_ISOSTRAIN_label
-           thisOutput => null()
-           thisSize   => null()
-         case (HOMOGENIZATION_RGC_ID)
-           outputName = HOMOGENIZATION_RGC_label
-           thisOutput => homogenization_RGC_output
-           thisSize   => homogenization_RGC_sizePostResult
-         case default
-           valid = .false.
-       end select
        write(FILEUNIT,'(/,a,/)')  '['//trim(homogenization_name(p))//']'
-       if (valid) then
-         write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName)
-         write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p)
-         if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID .and. &
-             homogenization_type(p) /= HOMOGENIZATION_ISOSTRAIN_ID) then
-           do e = 1,size(thisOutput(:,i))
-             write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
-           enddo
-         endif
-       endif
+       write(FILEUNIT,'(a)') '(type) n/a'
+       write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p)
+       
        i = thermal_typeInstance(p)                                                                      ! which instance of this thermal type
        valid = .true.                                                                                   ! assume valid
        select case(thermal_type(p))                                                                     ! split per thermal type
@@ -197,6 +197,7 @@ subroutine homogenization_init
            enddo
          endif
        endif
+       
        i = damage_typeInstance(p)                                                                       ! which instance of this damage type
        valid = .true.                                                                                   ! assume valid
        select case(damage_type(p))                                                                      ! split per damage type
@@ -229,7 +230,7 @@ subroutine homogenization_init
      endif
    enddo
    close(FILEUNIT)
- endif mainProcess2
+ endif mainProcess
 
  call config_deallocate('material.config/homogenization')
 
@@ -237,7 +238,7 @@ subroutine homogenization_init
 ! allocate and initialize global variables
  allocate(materialpoint_dPdF(3,3,3,3,theMesh%elem%nIPs,theMesh%nElems),       source=0.0_pReal)
  allocate(materialpoint_F0(3,3,theMesh%elem%nIPs,theMesh%nElems),             source=0.0_pReal)
- materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems)                          ! initialize to identity
+ materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems)                    ! initialize to identity
  allocate(materialpoint_F(3,3,theMesh%elem%nIPs,theMesh%nElems),              source=0.0_pReal)
  materialpoint_F = materialpoint_F0                                                                 ! initialize to identity
  allocate(materialpoint_subF0(3,3,theMesh%elem%nIPs,theMesh%nElems),          source=0.0_pReal)
@@ -252,18 +253,15 @@ subroutine homogenization_init
 
 !--------------------------------------------------------------------------------------------------
 ! allocate and initialize global state and postresutls variables
- homogenization_maxSizePostResults = 0
  thermal_maxSizePostResults        = 0
  damage_maxSizePostResults         = 0
  do p = 1,size(config_homogenization)
-   homogenization_maxSizePostResults = max(homogenization_maxSizePostResults,homogState       (p)%sizePostResults)
    thermal_maxSizePostResults        = max(thermal_maxSizePostResults,       thermalState     (p)%sizePostResults)
    damage_maxSizePostResults         = max(damage_maxSizePostResults        ,damageState      (p)%sizePostResults)
  enddo
 
  materialpoint_sizeResults = 1 &                                                                    ! grain count
-                           + 1 + homogenization_maxSizePostResults &                                ! homogSize & homogResult
-                               + thermal_maxSizePostResults        &
+                           + 1 + thermal_maxSizePostResults        &
                                + damage_maxSizePostResults         &
                            + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults &      ! crystallite size & crystallite results
                                                         + 1 + constitutive_plasticity_maxSizePostResults &     ! constitutive size & constitutive results
@@ -273,11 +271,6 @@ subroutine homogenization_init
  write(6,'(/,a)')   ' <<<+-  homogenization init  -+>>>'
 
  if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
-#ifdef TODO
-   write(6,'(a32,1x,7(i8,1x))')   'homogenization_state0:          ', shape(homogenization_state0)
-   write(6,'(a32,1x,7(i8,1x))')   'homogenization_subState0:       ', shape(homogenization_subState0)
-   write(6,'(a32,1x,7(i8,1x))')   'homogenization_state:           ', shape(homogenization_state)
-#endif
    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_dPdF:             ', shape(materialpoint_dPdF)
    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_F0:               ', shape(materialpoint_F0)
    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_F:                ', shape(materialpoint_F)
@@ -303,51 +296,6 @@ end subroutine homogenization_init
 !> @brief  parallelized calculation of stress and corresponding tangent at material points
 !--------------------------------------------------------------------------------------------------
 subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
- use numerics, only: &
-   subStepMinHomog, &
-   subStepSizeHomog, &
-   stepIncreaseHomog, &
-   nMPstate
- use FEsolving, only: &
-   FEsolving_execElem, &
-   FEsolving_execIP, &
-   terminallyIll
- use mesh, only: &
-   mesh_element
- use crystallite, only: &
-   crystallite_F0, &
-   crystallite_Fp0, &
-   crystallite_Fp, &
-   crystallite_Fi0, &
-   crystallite_Fi, &
-   crystallite_Lp0, &
-   crystallite_Lp, &
-   crystallite_Li0, &
-   crystallite_Li, &
-   crystallite_S0, &
-   crystallite_S, &
-   crystallite_partionedF0, &
-   crystallite_partionedF, &
-   crystallite_partionedFp0, &
-   crystallite_partionedLp0, &
-   crystallite_partionedFi0, &
-   crystallite_partionedLi0, &
-   crystallite_partionedS0, &
-   crystallite_dt, &
-   crystallite_requested, &
-   crystallite_stress, &
-   crystallite_stressTangent, &
-   crystallite_orientations
-#ifdef DEBUG
- use debug, only: &
-   debug_level, &
-   debug_homogenization, &
-   debug_levelBasic, &
-   debug_levelExtensive, &
-   debug_levelSelective, &
-   debug_e, &
-   debug_i
-#endif
 
  real(pReal), intent(in) :: dt                                                                      !< time increment
  logical,     intent(in) :: updateJaco                                                              !< initiating Jacobian update
@@ -640,14 +588,6 @@ end subroutine materialpoint_stressAndItsTangent
 !> @brief parallelized calculation of result array at material points
 !--------------------------------------------------------------------------------------------------
 subroutine materialpoint_postResults
- use FEsolving, only: &
-   FEsolving_execElem, &
-   FEsolving_execIP
- use mesh, only: &
-   mesh_element
- use crystallite, only: &
-   crystallite_sizePostResults, &
-   crystallite_postResults
 
  integer :: &
    thePos, &
@@ -697,12 +637,6 @@ end subroutine materialpoint_postResults
 !> @brief  partition material point def grad onto constituents
 !--------------------------------------------------------------------------------------------------
 subroutine partitionDeformation(ip,el)
- use mesh, only: &
-   mesh_element
- use crystallite, only: &
-   crystallite_partionedF
- use homogenization_mech_RGC, only: &
-   homogenization_RGC_partitionDeformation
 
  integer, intent(in) :: &
    ip, &                                                                                            !< integration point
@@ -719,7 +653,7 @@ subroutine partitionDeformation(ip,el)
                           materialpoint_subF(1:3,1:3,ip,el))
 
    case (HOMOGENIZATION_RGC_ID) chosenHomogenization
-     call homogenization_RGC_partitionDeformation(&
+     call mech_RGC_partitionDeformation(&
                          crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
                          materialpoint_subF(1:3,1:3,ip,el),&
                          ip, &
@@ -734,15 +668,6 @@ end subroutine partitionDeformation
 !> "happy" with result
 !--------------------------------------------------------------------------------------------------
 function updateState(ip,el)
- use mesh, only: &
-   mesh_element
- use crystallite, only: &
-   crystallite_P, &
-   crystallite_dPdF, &
-   crystallite_partionedF,&
-   crystallite_partionedF0
- use homogenization_mech_RGC, only: &
-   homogenization_RGC_updateState
  use thermal_adiabatic, only: &
    thermal_adiabatic_updateState
  use damage_local, only: &
@@ -758,14 +683,14 @@ function updateState(ip,el)
    case (HOMOGENIZATION_RGC_ID) chosenHomogenization
      updateState = &
        updateState .and. &
-        homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
-                                       crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
-                                       crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),&
-                                       materialpoint_subF(1:3,1:3,ip,el),&
-                                       materialpoint_subdt(ip,el), &
-                                       crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
-                                       ip, &
-                                       el)
+        mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
+                             crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
+                             crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),&
+                             materialpoint_subF(1:3,1:3,ip,el),&
+                             materialpoint_subdt(ip,el), &
+                             crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
+                             ip, &
+                             el)
  end select chosenHomogenization
 
  chosenThermal: select case (thermal_type(mesh_element(3,el)))
@@ -793,12 +718,6 @@ end function updateState
 !> @brief derive average stress and stiffness from constituent quantities
 !--------------------------------------------------------------------------------------------------
 subroutine averageStressAndItsTangent(ip,el)
- use mesh, only: &
-   mesh_element
- use crystallite, only: &
-   crystallite_P,crystallite_dPdF
- use homogenization_mech_RGC, only: &
-   homogenization_RGC_averageStressAndItsTangent
 
  integer, intent(in) :: &
    ip, &                                                                                            !< integration point
@@ -818,7 +737,7 @@ subroutine averageStressAndItsTangent(ip,el)
        homogenization_typeInstance(mesh_element(3,el)))
 
    case (HOMOGENIZATION_RGC_ID) chosenHomogenization
-     call homogenization_RGC_averageStressAndItsTangent(&
+     call mech_RGC_averageStressAndItsTangent(&
        materialpoint_P(1:3,1:3,ip,el), &
        materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
        crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), &
@@ -834,10 +753,6 @@ end subroutine averageStressAndItsTangent
 !> if homogenization_sizePostResults(i,e) > 0 !!
 !--------------------------------------------------------------------------------------------------
 function postResults(ip,el)
- use mesh, only: &
-   mesh_element
- use homogenization_mech_RGC, only: &
-   homogenization_RGC_postResults
  use thermal_adiabatic, only: &
    thermal_adiabatic_postResults
  use thermal_conduction, only: &
@@ -856,23 +771,12 @@ function postResults(ip,el)
    postResults
  integer :: &
    startPos, endPos ,&
-   of, instance, homog
+   homog
 
 
  postResults = 0.0_pReal
  startPos = 1
- endPos   = homogState(material_homogenizationAt(el))%sizePostResults
- chosenHomogenization: select case (homogenization_type(mesh_element(3,el)))
-
-   case (HOMOGENIZATION_RGC_ID) chosenHomogenization
-     instance = homogenization_typeInstance(material_homogenizationAt(el))
-     of = mappingHomogenization(1,ip,el)
-     postResults(startPos:endPos) = homogenization_RGC_postResults(instance,of)
-     
- end select chosenHomogenization
-
- startPos = endPos + 1
- endPos   = endPos + thermalState(material_homogenizationAt(el))%sizePostResults
+ endPos   = thermalState(material_homogenizationAt(el))%sizePostResults
  chosenThermal: select case (thermal_type(mesh_element(3,el)))
 
    case (THERMAL_adiabatic_ID) chosenThermal
@@ -905,14 +809,9 @@ end function postResults
 !--------------------------------------------------------------------------------------------------
 subroutine homogenization_results
 #if defined(PETSc) || defined(DAMASK_HDF5)
-  use results
-  use homogenization_mech_RGC
-  use HDF5_utilities
   use config, only: &
     config_name_homogenization => homogenization_name                                               ! anticipate logical name
-   
   use material, only: &
-    homogenization_typeInstance, &
     material_homogenization_type => homogenization_type
     
   integer :: p
diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90
index d3feac1eb..1cbe837d2 100644
--- a/src/homogenization_mech_RGC.f90
+++ b/src/homogenization_mech_RGC.f90
@@ -6,16 +6,7 @@
 !> @brief Relaxed grain cluster (RGC) homogenization scheme
 !> Nconstituents is defined as p x q x r (cluster)
 !--------------------------------------------------------------------------------------------------
-module homogenization_mech_RGC
- use prec
- use material
-
- implicit none
- private
- integer,           dimension(:,:),     allocatable,target, public :: &
-   homogenization_RGC_sizePostResult
- character(len=64), dimension(:,:),     allocatable,target, public :: &
-   homogenization_RGC_output                                                                        ! name of each post result output
+submodule(homogenization) homogenization_mech_RGC
 
  enum, bind(c) 
    enumerator :: &
@@ -28,7 +19,7 @@ module homogenization_mech_RGC
      magnitudemismatch_ID
  end enum
  
- type, private :: tParameters
+ type :: tParameters
    integer, dimension(:), allocatable :: &
      Nconstituents
    real(pReal) :: &
@@ -39,11 +30,11 @@ module homogenization_mech_RGC
      angles
    integer :: &
      of_debug = 0
-   integer(kind(undefined_ID)),         dimension(:),   allocatable   :: &
+   integer(kind(undefined_ID)),  dimension(:),   allocatable   :: &
      outputID
  end type tParameters
 
- type, private :: tRGCstate
+ type :: tRGCstate
    real(pReal), pointer,     dimension(:) :: &
      work, &
      penaltyEnergy
@@ -51,7 +42,7 @@ module homogenization_mech_RGC
      relaxationVector
  end type tRGCstate
 
- type, private :: tRGCdependentState
+ type :: tRGCdependentState
    real(pReal), allocatable,     dimension(:) :: &
      volumeDiscrepancy, &
      relaxationRate_avg, &
@@ -62,56 +53,25 @@ module homogenization_mech_RGC
      orientation
  end type tRGCdependentState
 
- type(tparameters),          dimension(:), allocatable, private :: &
+ type(tparameters),          dimension(:), allocatable :: &
    param
- type(tRGCstate),            dimension(:), allocatable, private :: &
+ type(tRGCstate),            dimension(:), allocatable :: &
    state, &
    state0
- type(tRGCdependentState),   dimension(:), allocatable, private :: &
+ type(tRGCdependentState),   dimension(:), allocatable :: &
    dependentState
 
- public :: &
-   homogenization_RGC_init, &
-   homogenization_RGC_partitionDeformation, &
-   homogenization_RGC_averageStressAndItsTangent, &
-   homogenization_RGC_updateState, &
-   homogenization_RGC_postResults, &
-   mech_RGC_results ! name suited for planned submodule situation
- private :: &
-   relaxationVector, &
-   interfaceNormal, &
-   getInterface, &
-   grain1to3, &
-   grain3to1, &
-   interface4to1, &
-   interface1to4
-
 contains
 
 !--------------------------------------------------------------------------------------------------
 !> @brief allocates all necessary fields, reads information from material configuration file
 !--------------------------------------------------------------------------------------------------
-subroutine homogenization_RGC_init()
- use debug, only: &
-#ifdef DEBUG
-   debug_i, &
-   debug_e, &
-#endif
-   debug_level, &
-   debug_homogenization, &
-   debug_levelBasic
- use math, only: &
-   math_EulerToR, &
-   INRAD
- use IO, only: &
-   IO_error
- use config, only: &
-   config_homogenization
+module subroutine mech_RGC_init
 
  integer :: &
    Ninstance, &
    h, i, &
-   NofMyHomog, outputSize, &
+   NofMyHomog, &
    sizeState, nIntFaceTot
 
  character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
@@ -139,9 +99,6 @@ subroutine homogenization_RGC_init()
  allocate(state0(Ninstance))
  allocate(dependentState(Ninstance))
  
- allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),Ninstance),source=0)
- allocate(homogenization_RGC_output(maxval(homogenization_Noutput),Ninstance))
-          homogenization_RGC_output=''
 
  do h = 1, size(homogenization_type)
    if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle
@@ -176,28 +133,20 @@ subroutine homogenization_RGC_init()
      
        case('constitutivework')
          outputID = constitutivework_ID
-         outputSize = 1
        case('penaltyenergy')
          outputID = penaltyenergy_ID
-         outputSize = 1
        case('volumediscrepancy')
          outputID = volumediscrepancy_ID
-         outputSize = 1
        case('averagerelaxrate')
          outputID = averagerelaxrate_ID
-         outputSize = 1
        case('maximumrelaxrate')
          outputID = maximumrelaxrate_ID
-         outputSize = 1
        case('magnitudemismatch')
          outputID = magnitudemismatch_ID
-         outputSize = 3
 
      end select
      
      if (outputID /= undefined_ID) then
-       homogenization_RGC_output(i,homogenization_typeInstance(h)) = outputs(i)
-       homogenization_RGC_sizePostResult(i,homogenization_typeInstance(h)) = outputSize
        prm%outputID = [prm%outputID , outputID]
      endif
      
@@ -211,7 +160,7 @@ subroutine homogenization_RGC_init()
              + size(['avg constitutive work ','average penalty energy'])
 
    homogState(h)%sizeState = sizeState
-   homogState(h)%sizePostResults = sum(homogenization_RGC_sizePostResult(:,homogenization_typeInstance(h)))
+   homogState(h)%sizePostResults = 0
    allocate(homogState(h)%state0   (sizeState,NofMyHomog), source=0.0_pReal)
    allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal)
    allocate(homogState(h)%state    (sizeState,NofMyHomog), source=0.0_pReal)
@@ -235,23 +184,17 @@ subroutine homogenization_RGC_init()
    
  enddo  
  
-end subroutine homogenization_RGC_init
+end subroutine mech_RGC_init
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief partitions the deformation gradient onto the constituents
 !--------------------------------------------------------------------------------------------------
-subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of)
-#ifdef DEBUG
- use debug, only: &
-   debug_level, &
-   debug_homogenization, &
-   debug_levelExtensive
-#endif
+module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
 
  real(pReal),   dimension (:,:,:), intent(out) :: F                                                 !< partioned F  per grain
  
- real(pReal),   dimension (:,:),   intent(in)  :: avgF                                              !< averaged F
+ real(pReal),   dimension (3,3),   intent(in)  :: avgF                                              !< averaged F
  integer,                          intent(in)  :: &
    instance, &
    of
@@ -291,45 +234,14 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of)
  
  end associate
 
-end subroutine homogenization_RGC_partitionDeformation
+end subroutine mech_RGC_partitionDeformation
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief update the internal state of the homogenization scheme and tell whether "done" and 
 ! "happy" with result
 !--------------------------------------------------------------------------------------------------
-function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
-#ifdef DEBUG
- use debug, only: &
-   debug_level, &
-   debug_homogenization,&
-   debug_levelExtensive
-#endif
- use math, only: &
-   math_invert2
- use numerics, only: &
-   absTol_RGC, &
-   relTol_RGC, &
-   absMax_RGC, &
-   relMax_RGC, &
-   pPert_RGC, &
-   maxdRelax_RGC, &
-   viscPower_RGC, &
-   viscModus_RGC, &
-   refRelaxRate_RGC
-
- real(pReal), dimension(:,:,:),     intent(in)    :: & 
-   P,&                                                                                              !< array of P
-   F,&                                                                                              !< array of F
-   F0                                                                                               !< array of initial F
- real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF                                              !< array of current grain stiffness
- real(pReal), dimension(3,3),       intent(in) :: avgF                                              !< average F
- real(pReal),                       intent(in) :: dt                                                !< time increment
- integer,                           intent(in) :: &
-   ip, &                                                                                            !< integration point number
-   el                                                                                               !< element number
-
- logical, dimension(2) :: homogenization_RGC_updateState
+module procedure mech_RGC_updateState
 
  integer, dimension(4) :: intFaceN,intFaceP,faceID
  integer, dimension(3) :: nGDim,iGr3N,iGr3P
@@ -347,7 +259,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
 #endif
  
  zeroTimeStep: if(dEq0(dt)) then
-   homogenization_RGC_updateState = .true.                                                          ! pretend everything is fine and return
+   mech_RGC_updateState = .true.                                                                    ! pretend everything is fine and return
    return                                                                    
  endif zeroTimeStep
 
@@ -466,12 +378,12 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  endif
 #endif
  
- homogenization_RGC_updateState = .false.
+ mech_RGC_updateState = .false.
  
 !--------------------------------------------------------------------------------------------------
 !  If convergence reached => done and happy
  if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then 
-   homogenization_RGC_updateState = .true.
+   mech_RGC_updateState = .true.
 #ifdef DEBUG
     if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 &
         .and. prm%of_debug == of) write(6,'(1x,a55,/)')'... done and happy'
@@ -513,7 +425,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
 !--------------------------------------------------------------------------------------------------
 ! if residual blows-up => done but unhappy
  elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then                            ! try to restart when residual blows up exceeding maximum bound
-   homogenization_RGC_updateState = [.true.,.false.]                                                ! with direct cut-back
+   mech_RGC_updateState = [.true.,.false.]                                                          ! with direct cut-back
    
 #ifdef DEBUG
    if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 &
@@ -656,9 +568,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
 !--------------------------------------------------------------------------------------------------
 ! ... of the numerical viscosity traction "rmatrix"
  allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
- forall (i=1:3*nIntFaceTot) &
+ do i=1,3*nIntFaceTot
    rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* &                              ! tangent due to numerical viscosity traction appears
                   (abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal)               ! only in the main diagonal term 
+ enddo
                                                                         
 #ifdef DEBUG
  if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then
@@ -710,7 +623,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  enddo; enddo
  stt%relaxationVector(:,of) = relax + drelax                                                        ! Updateing the state variable for the next iteration
  if (any(abs(drelax) > maxdRelax_RGC)) then                                                         ! Forcing cutback when the incremental change of relaxation vector becomes too large
-   homogenization_RGC_updateState = [.true.,.false.]
+   mech_RGC_updateState = [.true.,.false.]
    !$OMP CRITICAL (write2out)
    write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
    write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax))
@@ -736,10 +649,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  !> @brief calculate stress-like penalty due to deformation mismatch
  !--------------------------------------------------------------------------------------------------
  subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of)
-  use math, only: &
-    math_civita
-  use numerics, only: &
-    xSmoo_RGC
 
   real(pReal),   dimension (:,:,:), intent(out) :: rPen                                             !< stress-like penalty
   real(pReal),   dimension (:,:),   intent(out) :: nMis                                             !< total amount of mismatch
@@ -852,13 +761,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  !> @brief calculate stress-like penalty due to volume discrepancy 
  !--------------------------------------------------------------------------------------------------
  subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of)
-  use math, only: &
-    math_det33, &
-    math_inv33
-  use numerics, only: &
-    maxVolDiscr_RGC,&
-    volDiscrMod_RGC,&
-    volDiscrPow_RGC
  
   real(pReal), dimension (:,:,:), intent(out) :: vPen                                               ! stress-like penalty due to volume
   real(pReal),                    intent(out) :: vDiscrep                                           ! total volume discrepancy
@@ -907,8 +809,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  ! deformation
  !--------------------------------------------------------------------------------------------------
  function surfaceCorrection(avgF,instance,of)
-  use math, only: &
-    math_invert33
   
   real(pReal), dimension(3)               :: surfaceCorrection
 
@@ -940,8 +840,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
  !--------------------------------------------------------------------------------------------------
  function equivalentModuli(grainID,ip,el)
-  use constitutive, only: &
-    constitutive_homogenizedC
  
   real(pReal), dimension(2)    :: equivalentModuli
 
@@ -1012,13 +910,13 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
   
  end subroutine grainDeformation
  
-end function homogenization_RGC_updateState
+end procedure mech_RGC_updateState
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief derive average stress and stiffness from constituent quantities 
 !--------------------------------------------------------------------------------------------------
-subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
+module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
 
  real(pReal), dimension (3,3),        intent(out) :: avgP                                           !< average stress at material point
  real(pReal), dimension (3,3,3,3),    intent(out) :: dAvgPdAvgF                                     !< average stiffness at material point
@@ -1030,68 +928,19 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,
  avgP       = sum(P,3)   /real(product(param(instance)%Nconstituents),pReal)
  dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal)
 
-end subroutine homogenization_RGC_averageStressAndItsTangent
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of homogenization results for post file inclusion 
-!--------------------------------------------------------------------------------------------------
-pure function homogenization_RGC_postResults(instance,of) result(postResults)
-
- integer, intent(in) :: &
-   instance, &
-   of
-   
- integer :: &
-   o,c
- real(pReal), dimension(sum(homogenization_RGC_sizePostResult(:,instance))) :: &
-   postResults
-
- associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
-
- c = 0
-
- outputsLoop: do o = 1,size(prm%outputID)
-   select case(prm%outputID(o))
-   
-     case (constitutivework_ID)
-       postResults(c+1) = stt%work(of)
-       c = c + 1
-     case (magnitudemismatch_ID)
-       postResults(c+1:c+3) = dst%mismatch(1:3,of)
-       c = c + 3
-     case (penaltyenergy_ID)
-       postResults(c+1) = stt%penaltyEnergy(of)
-       c = c + 1
-     case (volumediscrepancy_ID)
-       postResults(c+1) = dst%volumeDiscrepancy(of)
-       c = c + 1
-     case (averagerelaxrate_ID)
-       postResults(c+1) = dst%relaxationrate_avg(of)
-       c = c + 1
-     case (maximumrelaxrate_ID)
-       postResults(c+1) = dst%relaxationrate_max(of)
-       c = c + 1
-   end select
-   
- enddo outputsLoop
- 
- end associate
- 
-end function homogenization_RGC_postResults
+end subroutine mech_RGC_averageStressAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 ! ToDo: check wheter units are correct
 !--------------------------------------------------------------------------------------------------
-subroutine mech_RGC_results(instance,group)
+module subroutine mech_RGC_results(instance,group)
 #if defined(PETSc) || defined(DAMASK_HDF5)
-  use results, only: &
-    results_writeDataset
 
-  integer, intent(in) :: instance
-  character(len=*) :: group
+  integer,          intent(in) :: instance
+  character(len=*), intent(in) :: group
+  
   integer :: o
   
   associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
@@ -1122,8 +971,8 @@ subroutine mech_RGC_results(instance,group)
   end associate
   
 #else
-  integer, intent(in) :: instance
-  character(len=*) :: group
+  integer,          intent(in) :: instance
+  character(len=*), intent(in) :: group
 #endif
 
 end subroutine mech_RGC_results
@@ -1262,7 +1111,7 @@ integer pure function interface4to1(iFace4D, nGDim)
      else
        interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) &
                      + nGDim(3)*nGDim(1)*(iFace4D(3)-1) &
-                     + (nGDim(1)-1)*nGDim(2)*nGDim(3)                                               ! total number of interfaces normal //e1
+                     + (nGDim(1)-1)*nGDim(2)*nGDim(3)                                               ! total # of interfaces normal || e1
      endif
 
    case(3)
@@ -1271,8 +1120,8 @@ integer pure function interface4to1(iFace4D, nGDim)
      else
        interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) &
                      + nGDim(1)*nGDim(2)*(iFace4D(4)-1) &
-                     + (nGDim(1)-1)*nGDim(2)*nGDim(3) &                                             ! total number of interfaces normal //e1
-                     + nGDim(1)*(nGDim(2)-1)*nGDim(3)                                               ! total number of interfaces normal //e2
+                     + (nGDim(1)-1)*nGDim(2)*nGDim(3) &                                             ! total # of interfaces normal || e1
+                     + nGDim(1)*(nGDim(2)-1)*nGDim(3)                                               ! total # of interfaces normal || e2
      endif
 
    case default
@@ -1296,23 +1145,23 @@ pure function interface1to4(iFace1D, nGDim)
 
 !--------------------------------------------------------------------------------------------------
 ! compute the total number of interfaces, which ...
- nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), &                                                      ! ... normal //e1
-             nGDim(1)*(nGDim(2)-1)*nGDim(3), &                                                      ! ... normal //e2
-             nGDim(1)*nGDim(2)*(nGDim(3)-1)]                                                        ! ... normal //e3
+ nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), &                                                      ! ... normal || e1
+             nGDim(1)*(nGDim(2)-1)*nGDim(3), &                                                      ! ... normal || e2
+             nGDim(1)*nGDim(2)*(nGDim(3)-1)]                                                        ! ... normal || e3
 
 !--------------------------------------------------------------------------------------------------
 ! get the corresponding interface ID in 4D (normal and local position)
- if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then                                                 ! interface with normal //e1
+ if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then                                                 ! interface with normal || e1
    interface1to4(1) = 1
    interface1to4(3) = mod((iFace1D-1),nGDim(2))+1
    interface1to4(4) = mod(int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)),nGDim(3))+1
    interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal))+1
- elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then                   ! interface with normal //e2
+ elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then                   ! interface with normal || e2
    interface1to4(1) = 2
    interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1
    interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1
    interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1
- elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal //e3
+ elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal || e3
    interface1to4(1) = 3
    interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1
    interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1
@@ -1322,4 +1171,4 @@ pure function interface1to4(iFace1D, nGDim)
 end function interface1to4
 
 
-end module homogenization_mech_RGC
+end submodule homogenization_mech_RGC
diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90
index 7dd7bad7d..c1eba821c 100644
--- a/src/homogenization_mech_isostrain.f90
+++ b/src/homogenization_mech_isostrain.f90
@@ -6,8 +6,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(homogenization) homogenization_mech_isostrain
 
-  implicit none
-
   enum, bind(c) 
     enumerator :: &
       parallel_ID, &
@@ -30,14 +28,6 @@ contains
 !> @brief allocates all neccessary fields, reads information from material configuration file
 !--------------------------------------------------------------------------------------------------
 module subroutine mech_isostrain_init
-  use debug, only: &
-    debug_HOMOGENIZATION, &
-    debug_level, &
-    debug_levelBasic
-  use IO, only: &
-    IO_error
-  use config, only: &
-    config_homogenization
 
   integer :: &
     Ninstance, &
diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90
index e7a5a12e6..d5b24e5d1 100644
--- a/src/homogenization_mech_none.f90
+++ b/src/homogenization_mech_none.f90
@@ -6,20 +6,12 @@
 !--------------------------------------------------------------------------------------------------
 submodule(homogenization) homogenization_mech_none
 
-  implicit none
-
 contains
 
 !--------------------------------------------------------------------------------------------------
 !> @brief allocates all neccessary fields, reads information from material configuration file
 !--------------------------------------------------------------------------------------------------
 module subroutine mech_none_init
-  use debug, only: &
-    debug_HOMOGENIZATION, &
-    debug_level, &
-    debug_levelBasic
-  use config, only: &
-    config_homogenization
 
   integer :: &
     Ninstance, &
@@ -27,7 +19,7 @@ module subroutine mech_none_init
     NofMyHomog
  
   write(6,'(/,a)')   ' <<<+-  homogenization_'//HOMOGENIZATION_NONE_label//' init  -+>>>'
- 
+
   Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID)
   if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) &
     write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90
index 7db1f5f7f..cb4b3cbae 100644
--- a/src/plastic_dislotwin.f90
+++ b/src/plastic_dislotwin.f90
@@ -1067,8 +1067,8 @@ subroutine plastic_dislotwin_results(instance,group)
       case (rho_dip_ID)
         call results_writeDataset(group,stt%rho_dip,'rho_dip',&
                                   'dislocation dipole density''1/m²')
-      case (dot_gamma_sl_ID)
-        call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',&
+      case (gamma_sl_ID)
+        call results_writeDataset(group,stt%gamma_sl,'gamma_sl',&
                                   'plastic shear','1')
       case (Lambda_sl_ID)
         call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
diff --git a/src/results.f90 b/src/results.f90
index c5582b927..05db831f7 100644
--- a/src/results.f90
+++ b/src/results.f90
@@ -299,18 +299,26 @@ end subroutine results_writeVectorDataset_real
 !--------------------------------------------------------------------------------------------------
 subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit)
 
-  character(len=*), intent(in)                      :: label,group,description
-  character(len=*), intent(in), optional            :: SIunit
-  real(pReal),      intent(inout), dimension(:,:,:) :: dataset
+  character(len=*), intent(in)                   :: label,group,description
+  character(len=*), intent(in), optional         :: SIunit
+  real(pReal),      intent(in), dimension(:,:,:) :: dataset
   
+  integer :: i 
   integer(HID_T) :: groupHandle
- 
+  real(pReal), dimension(:,:,:), allocatable :: dataset_transposed
+
+
+  allocate(dataset_transposed,mold=dataset)
+  do i=1,size(dataset,3)
+    dataset_transposed(1:3,1:3,i) = transpose(dataset(1:3,1:3,i))
+  enddo
+
   groupHandle = results_openGroup(group)
   
 #ifdef PETSc
-  call HDF5_write(groupHandle,dataset,label,.true.)
+  call HDF5_write(groupHandle,dataset_transposed,label,.true.)
 #else
-  call HDF5_write(groupHandle,dataset,label,.false.)
+  call HDF5_write(groupHandle,dataset_transposed,label,.false.)
 #endif
   
   if (HDF5_objectExists(groupHandle,label)) &