diff --git a/processing/post/marc_to_vtk.py b/processing/post/marc_to_vtk.py
new file mode 100755
index 000000000..a232d1219
--- /dev/null
+++ b/processing/post/marc_to_vtk.py
@@ -0,0 +1,167 @@
+#!/usr/bin/env python2.7
+# -*- coding: UTF-8 no BOM -*-
+
+import os,re
+import argparse
+import damask
+import vtk, numpy as np
+
+scriptName = os.path.splitext(os.path.basename(__file__))[0]
+scriptID   = ' '.join([scriptName, damask.version])
+
+parser = argparse.ArgumentParser(description='Convert from Marc input file format to VTK', version = scriptID)
+parser.add_argument('filename', type=str, nargs='+', help='files to convert')
+
+args = parser.parse_args()
+files = args.filename
+if type(files) is str:
+  files = [files]
+
+    
+for f in files:
+  with open(f, 'r') as marcfile:
+    marctext = marcfile.read();
+  # Extract connectivity chunk from file...
+  connectivity_text = re.findall(r'connectivity[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
+  connectivity_lines = re.split(r'[\n\r]+', connectivity_text, flags=(re.MULTILINE | re.DOTALL))
+  connectivity_header = connectivity_lines[0]
+  connectivity_lines = connectivity_lines[1:]
+  # Construct element map
+  elements = dict(map(lambda line: 
+    (
+      int(line[0:10]), # index
+      { 
+        'type':  int(line[10:20]), 
+        'verts': list(map(int, re.split(r' +', line[20:].strip())))
+      }
+    ), connectivity_lines)) 
+  # Extract coordinate chunk from file
+  coordinates_text = re.findall(r'coordinates[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
+  coordinates_lines = re.split(r'[\n\r]+', coordinates_text, flags=(re.MULTILINE | re.DOTALL))
+  coordinates_header = coordinates_lines[0]
+  coordinates_lines = coordinates_lines[1:]
+  # marc input file does not use "e" in scientific notation, this adds it and converts
+  fl_format = lambda string: float(re.sub(r'(\d)([\+\-])', r'\1e\2', string))
+  # Construct coordinate map
+  coordinates = dict(map(lambda line:
+    (
+      int(line[0:10]),
+      np.array([
+        fl_format(line[10:30]), 
+        fl_format(line[30:50]), 
+        fl_format(line[50:70])
+      ])
+    ), coordinates_lines))
+  
+  # Subdivide volumes
+  grid = vtk.vtkUnstructuredGrid()
+  vertex_count = len(coordinates)
+  edge_to_vert = dict() # when edges are subdivided, a new vertex in the middle is produced and placed in here
+  ordered_pair = lambda a, b: (a, b) if a < b else (b, a) # edges are bidirectional
+
+  def subdivide_edge(vert1, vert2):
+    edge = ordered_pair(vert1, vert2)
+    
+    if edge in edge_to_vert:
+      return edge_to_vert[edge]
+    
+    newvert = len(coordinates) + 1
+    coordinates[newvert] = 0.5 * (coordinates[vert1] + coordinates[vert2]) # Average
+    edge_to_vert[edge] = newvert;
+    return newvert;
+  
+  
+  
+  for el_id in range(1, len(elements) + 1):
+    el = elements[el_id]
+    if el['type'] == 7:
+      # Hexahedron, subdivided
+      
+      # There may be a better way to iterate over these, but this is consistent
+      # with the ordering scheme provided at https://damask.mpie.de/pub/Documentation/ElementType
+      
+      subverts = np.zeros((3,3,3), dtype=int)
+      # Get corners
+      subverts[0, 0, 0] = el['verts'][0]
+      subverts[2, 0, 0] = el['verts'][1]
+      subverts[2, 2, 0] = el['verts'][2]
+      subverts[0, 2, 0] = el['verts'][3]
+      subverts[0, 0, 2] = el['verts'][4]
+      subverts[2, 0, 2] = el['verts'][5]
+      subverts[2, 2, 2] = el['verts'][6]
+      subverts[0, 2, 2] = el['verts'][7]
+      
+      # lower edges
+      subverts[1, 0, 0] = subdivide_edge(subverts[0, 0, 0], subverts[2, 0, 0])
+      subverts[2, 1, 0] = subdivide_edge(subverts[2, 0, 0], subverts[2, 2, 0])
+      subverts[1, 2, 0] = subdivide_edge(subverts[2, 2, 0], subverts[0, 2, 0])
+      subverts[0, 1, 0] = subdivide_edge(subverts[0, 2, 0], subverts[0, 0, 0])
+      
+      # middle edges
+      subverts[0, 0, 1] = subdivide_edge(subverts[0, 0, 0], subverts[0, 0, 2])
+      subverts[2, 0, 1] = subdivide_edge(subverts[2, 0, 0], subverts[2, 0, 2])
+      subverts[2, 2, 1] = subdivide_edge(subverts[2, 2, 0], subverts[2, 2, 2])
+      subverts[0, 2, 1] = subdivide_edge(subverts[0, 2, 0], subverts[0, 2, 2])
+      
+      # top edges
+      subverts[1, 0, 2] = subdivide_edge(subverts[0, 0, 2], subverts[2, 0, 2])
+      subverts[2, 1, 2] = subdivide_edge(subverts[2, 0, 2], subverts[2, 2, 2])
+      subverts[1, 2, 2] = subdivide_edge(subverts[2, 2, 2], subverts[0, 2, 2])
+      subverts[0, 1, 2] = subdivide_edge(subverts[0, 2, 2], subverts[0, 0, 2])
+      
+      # then faces...  The edge_to_vert addition is due to there being two ways
+      # to calculate a face, depending which opposite vertices are used to subdivide
+      subverts[1, 1, 0] = subdivide_edge(subverts[1, 0, 0], subverts[1, 2, 0])
+      edge_to_vert[ordered_pair(subverts[0, 1, 0], subverts[2, 1, 0])] = subverts[1, 1, 0]
+      
+      subverts[1, 0, 1] = subdivide_edge(subverts[1, 0, 0], subverts[1, 0, 2])
+      edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[2, 0, 1])] = subverts[1, 0, 1]
+      
+      subverts[2, 1, 1] = subdivide_edge(subverts[2, 1, 0], subverts[2, 1, 2])
+      edge_to_vert[ordered_pair(subverts[2, 0, 1], subverts[2, 2, 1])] = subverts[2, 1, 1]
+      
+      subverts[1, 2, 1] = subdivide_edge(subverts[1, 2, 0], subverts[1, 2, 2])
+      edge_to_vert[ordered_pair(subverts[0, 2, 1], subverts[2, 2, 1])] = subverts[1, 2, 1]
+      
+      subverts[0, 1, 1] = subdivide_edge(subverts[0, 1, 0], subverts[0, 1, 2])
+      edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[0, 2, 1])] = subverts[0, 1, 1]
+      
+      subverts[1, 1, 2] = subdivide_edge(subverts[1, 0, 2], subverts[1, 2, 2])
+      edge_to_vert[ordered_pair(subverts[0, 1, 2], subverts[2, 1, 2])] = subverts[1, 1, 2]
+      
+      # and finally the center.  There are three ways to calculate, but elements should
+      # not intersect, so the edge_to_vert part isn't needed here.
+      subverts[1, 1, 1] = subdivide_edge(subverts[1, 1, 0], subverts[1, 1, 2])
+      
+      
+      # Now make the hexahedron subelements
+      # order in which vtk expects vertices for a hexahedron
+      order = np.array([(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1),(1,0,1),(1,1,1),(0,1,1)])
+      for z in range(2):
+        for y in range(2):
+          for x in range(2):
+            hex_ = vtk.vtkHexahedron()
+            for vert_id in range(8):
+              coord = order[vert_id] + (x, y, z)
+              hex_.GetPointIds().SetId(vert_id, subverts[coord[0], coord[1], coord[2]] - 1) # minus one, since vtk starts at zero but marc starts at one
+            grid.InsertNextCell(hex_.GetCellType(), hex_.GetPointIds())
+            
+      
+    else:
+      damask.util.croak('Unsupported Marc element type: {} (skipping)'.format(el['type']))
+  
+  # Load all points
+  points = vtk.vtkPoints()
+  for point in range(1, len(coordinates) + 1): # marc indices start at 1
+    points.InsertNextPoint(coordinates[point].tolist())
+
+  grid.SetPoints(points)
+  
+  # grid now contains the elements from the given marc file
+  writer = vtk.vtkXMLUnstructuredGridWriter()
+  writer.SetFileName(re.sub(r'\..+', ".vtu", f))  # *.vtk extension does not work in paraview
+  #writer.SetCompressorTypeToZLib()
+
+  if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(grid)
+  else:                          writer.SetInputData(grid)
+  writer.Write()
diff --git a/processing/post/vtk_addGridData.py b/processing/post/vtk_addGridData.py
new file mode 100755
index 000000000..e0c274dc7
--- /dev/null
+++ b/processing/post/vtk_addGridData.py
@@ -0,0 +1,206 @@
+#!/usr/bin/env python2.7
+# -*- coding: UTF-8 no BOM -*-
+
+import os,vtk
+import damask
+from vtk.util import numpy_support
+from collections import defaultdict
+from optparse import OptionParser
+
+scriptName = os.path.splitext(os.path.basename(__file__))[0]
+scriptID   = ' '.join([scriptName,damask.version])
+
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+
+msg = "Add scalars, vectors, and/or an RGB tuple from"
+msg += "an ASCIItable to existing VTK grid (.vtr/.vtk/.vtu)."
+parser = OptionParser(option_class=damask.extendableOption,
+                      usage='%prog options [file[s]]',
+                      description = msg,
+                      version = scriptID)
+
+parser.add_option(      '--vtk',
+                  dest = 'vtk',
+                  type = 'string', metavar = 'string',
+                  help = 'VTK file name')
+parser.add_option(      '--inplace',
+                  dest = 'inplace',
+                  action = 'store_true',
+                  help = 'modify VTK file in-place')
+parser.add_option('-r', '--render',
+                  dest = 'render',
+                  action = 'store_true',
+                  help = 'open output in VTK render window')
+parser.add_option('-d', '--data',
+                  dest = 'data',
+                  action = 'extend', metavar = '<string LIST>',
+                  help = 'scalar/vector value(s) label(s)')
+parser.add_option('-t', '--tensor',
+                  dest = 'tensor',
+                  action = 'extend', metavar = '<string LIST>',
+                  help = 'tensor (3x3) value label(s)')
+parser.add_option('-c', '--color',
+                  dest = 'color',
+                  action = 'extend', metavar = '<string LIST>',
+                  help = 'RGB color tuple label')
+
+parser.set_defaults(data = [],
+                    tensor = [],
+                    color = [],
+                    inplace = False,
+                    render = False,
+)
+
+(options, filenames) = parser.parse_args()
+
+if not options.vtk:                 parser.error('No VTK file specified.')
+if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
+
+if os.path.splitext(options.vtk)[1] == '.vtr':
+  reader = vtk.vtkXMLRectilinearGridReader()
+  reader.SetFileName(options.vtk)
+  reader.Update()
+  rGrid = reader.GetOutput()
+  writer = vtk.vtkXMLRectilinearGridWriter()
+  writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtr' if options.inplace else '_added.vtr'))
+elif os.path.splitext(options.vtk)[1] == '.vtk':
+  reader = vtk.vtkGenericDataObjectReader()
+  reader.SetFileName(options.vtk)
+  reader.Update()
+  rGrid = reader.GetRectilinearGridOutput()
+  writer = vtk.vtkXMLRectilinearGridWriter()
+  writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtr' if options.inplace else '_added.vtr'))
+elif os.path.splitext(options.vtk)[1] == '.vtu':
+  reader = vtk.vtkXMLUnstructuredGridReader()
+  reader.SetFileName(options.vtk)
+  reader.Update()
+  rGrid = reader.GetOutput()
+  writer = vtk.vtkXMLUnstructuredGridWriter()
+  writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtu' if options.inplace else '_added.vtu'))
+else:
+  parser.error('Unsupported VTK file type extension.')
+
+Npoints = rGrid.GetNumberOfPoints()
+Ncells  = rGrid.GetNumberOfCells()
+
+damask.util.croak('{}: {} points and {} cells...'.format(options.vtk,Npoints,Ncells))
+
+# --- loop over input files -------------------------------------------------------------------------
+
+if filenames == []: filenames = [None]
+
+for name in filenames:
+  try:    table = damask.ASCIItable(name = name,
+                                    buffered = False,
+                                    readonly = True)
+  except: continue
+  damask.util.report(scriptName, name)
+
+# --- interpret header ----------------------------------------------------------------------------
+
+  table.head_read()
+
+  remarks = []
+  errors  = []
+  VTKarray = {}
+  active = defaultdict(list)
+
+  for datatype,dimension,label in [['data',99,options.data],
+                                   ['tensor',9,options.tensor],
+                                   ['color' ,3,options.color],
+                                   ]:
+    for i,dim in enumerate(table.label_dimension(label)):
+      me = label[i]
+      if dim == -1:         remarks.append('{} "{}" not found...'.format(datatype,me))
+      elif dim > dimension: remarks.append('"{}" not of dimension {}...'.format(me,dimension))
+      else:
+        remarks.append('adding {} "{}"...'.format(datatype,me))
+        active[datatype].append(me)
+
+  if remarks != []: damask.util.croak(remarks)
+  if errors  != []:
+    damask.util.croak(errors)
+    table.close(dismiss = True)
+    continue
+
+# ------------------------------------------ process data ---------------------------------------
+
+  table.data_readArray([item for sublist in active.values() for item in sublist])                 # read all requested data
+
+  for datatype,labels in active.items():                                                          # loop over scalar,color
+    for me in labels:                                                                             # loop over all requested items
+      VTKtype = vtk.VTK_DOUBLE
+      VTKdata = table.data[:, table.label_indexrange(me)].copy()                                  # copy to force contiguous layout
+
+      if datatype == 'color':
+        VTKtype = vtk.VTK_UNSIGNED_CHAR
+        VTKdata = (VTKdata*255).astype(int)                                                       # translate to 0..255 UCHAR
+      elif datatype == 'tensor':
+        VTKdata[:,1] = VTKdata[:,3] = 0.5*(VTKdata[:,1]+VTKdata[:,3])
+        VTKdata[:,2] = VTKdata[:,6] = 0.5*(VTKdata[:,2]+VTKdata[:,6])
+        VTKdata[:,5] = VTKdata[:,7] = 0.5*(VTKdata[:,5]+VTKdata[:,7])
+
+      VTKarray[me] = numpy_support.numpy_to_vtk(num_array=VTKdata,deep=True,array_type=VTKtype)
+      VTKarray[me].SetName(me)
+
+  table.close()                                                                                     # close input ASCII table
+
+# ------------------------------------------ add data ---------------------------------------
+
+  if   len(table.data) == Npoints:  mode = 'point'
+  elif len(table.data) == Ncells:   mode = 'cell'
+  else:
+    damask.util.croak('Data count is incompatible with grid...')
+    continue
+
+  damask.util.croak('{} mode...'.format(mode))
+
+  for datatype,labels in active.items():                                                            # loop over scalar,color
+    if datatype == 'color':
+      if   mode == 'cell':  rGrid.GetCellData().SetScalars(VTKarray[active['color'][0]])
+      elif mode == 'point': rGrid.GetPointData().SetScalars(VTKarray[active['color'][0]])
+    for me in labels:                                                                               # loop over all requested items
+      if   mode == 'cell':  rGrid.GetCellData().AddArray(VTKarray[me])
+      elif mode == 'point': rGrid.GetPointData().AddArray(VTKarray[me])
+
+  rGrid.Modified()
+  if vtk.VTK_MAJOR_VERSION <= 5: rGrid.Update()
+
+# ------------------------------------------ output result ---------------------------------------
+
+  writer.SetDataModeToBinary()
+  writer.SetCompressorTypeToZLib()
+  if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid)
+  else:                          writer.SetInputData(rGrid)
+  writer.Write()
+
+# ------------------------------------------ render result ---------------------------------------
+
+if options.render:
+  mapper = vtk.vtkDataSetMapper()
+  mapper.SetInputData(rGrid)
+  actor = vtk.vtkActor()
+  actor.SetMapper(mapper)
+
+# Create the graphics structure. The renderer renders into the
+# render window. The render window interactor captures mouse events
+# and will perform appropriate camera or actor manipulation
+# depending on the nature of the events.
+
+  ren = vtk.vtkRenderer()
+
+  renWin = vtk.vtkRenderWindow()
+  renWin.AddRenderer(ren)
+
+  ren.AddActor(actor)
+  ren.SetBackground(1, 1, 1)
+  renWin.SetSize(200, 200)
+
+  iren = vtk.vtkRenderWindowInteractor()
+  iren.SetRenderWindow(renWin)
+
+  iren.Initialize()
+  renWin.Render()
+  iren.Start()