From 46f0ad052e717ad016ce2a52d30742e054c94b91 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 30 May 2019 15:35:45 +0200
Subject: [PATCH] direct support for vtk output

- geom_check can now handle multiple files
- microstructure index is stored as integer in vtk file
---
 processing/pre/geom_check.py | 39 ++++++++++++++++++++++++
 processing/pre/geom_check.sh | 25 ---------------
 python/damask/geom.py        | 59 ++++++++++++++++++++++++++++++++++++
 3 files changed, 98 insertions(+), 25 deletions(-)
 create mode 100755 processing/pre/geom_check.py
 delete mode 100755 processing/pre/geom_check.sh

diff --git a/processing/pre/geom_check.py b/processing/pre/geom_check.py
new file mode 100755
index 000000000..a79c7e1a4
--- /dev/null
+++ b/processing/pre/geom_check.py
@@ -0,0 +1,39 @@
+#!/usr/bin/env python3
+
+import os
+import sys
+from io import StringIO
+from optparse import OptionParser
+
+import damask
+
+
+scriptName = os.path.splitext(os.path.basename(__file__))[0]
+scriptID   = ' '.join([scriptName,damask.version])
+
+
+#--------------------------------------------------------------------------------------------------
+#                                MAIN
+#--------------------------------------------------------------------------------------------------
+
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog [geomfile(s)]', description = """
+Writes vtk file for visualization.
+
+""", version = scriptID)
+
+(options, filenames) = parser.parse_args()
+
+
+if filenames == []: filenames = [None]
+
+for name in filenames:
+  damask.util.report(scriptName,name)
+
+  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
+
+  damask.util.croak(geom)
+  
+  if name is None:
+    sys.stdout.write(geom.to_vtk())
+  else:
+    geom.to_vtk(os.path.splitext(name)[0])
diff --git a/processing/pre/geom_check.sh b/processing/pre/geom_check.sh
deleted file mode 100755
index fbae91043..000000000
--- a/processing/pre/geom_check.sh
+++ /dev/null
@@ -1,25 +0,0 @@
-#!/usr/bin/env bash
-
-if [[ "$1" == "-f" || "$1" == "--float" ]]; then
-  shift
-  arg='--float'
-else
-  arg=''
-fi
-
-for geom in "$@"
-do
-  geom_toTable $arg \
-  < $geom \
-  | \
-  vtk_rectilinearGrid > ${geom%.*}.vtk
-
-  geom_toTable $arg \
-  < $geom \
-  | \
-  vtk_addRectilinearGridData \
-    --vtk ${geom%.*}.vtk \
-    --data microstructure \
-
-  rm ${geom%.*}.vtk
-done
diff --git a/python/damask/geom.py b/python/damask/geom.py
index e422516a7..365e54ade 100644
--- a/python/damask/geom.py
+++ b/python/damask/geom.py
@@ -1,8 +1,12 @@
+import os
 from io import StringIO
 
 import numpy as np
+import vtk
+from vtk.util import numpy_support
 
 from . import util
+from . import version
 
 
 class Geom():
@@ -197,6 +201,61 @@ class Geom():
     np.savetxt(fname,
                self.microstructure.reshape([grid[0],np.prod(grid[1:])],order='F').T,
                header='\n'.join(header), fmt=format_string, comments='')
+
+
+               
+  def to_vtk(self,fname=None):
+    """Generates vtk file. If file name is given, stored in file otherwise returned as string"""
+    grid = self.get_grid() + np.ones(3,dtype=int)
+    size = self.get_size()
+    origin = self.get_origin()
+
+    coords = [
+              np.linspace(0,size[0],grid[0]) - origin[0],
+              np.linspace(0,size[1],grid[1]) - origin[1],
+              np.linspace(0,size[2],grid[2]) - origin[2]
+             ]
+    
+    rGrid = vtk.vtkRectilinearGrid()
+    coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()]
+
+    rGrid.SetDimensions(*grid)
+    for d,coord in enumerate(coords):
+      for c in coord:
+        coordArray[d].InsertNextValue(c)
+
+    rGrid.SetXCoordinates(coordArray[0])
+    rGrid.SetYCoordinates(coordArray[1])
+    rGrid.SetZCoordinates(coordArray[2])
+    
+    ms = numpy_support.numpy_to_vtk(num_array=self.microstructure.flatten(order='F'),array_type=vtk.VTK_INT)
+    ms.SetName('microstructure')
+    rGrid.GetCellData().AddArray(ms)
+
+
+    if fname is not None:
+      writer = vtk.vtkXMLRectilinearGridWriter()
+      writer.SetCompressorTypeToZLib()
+      writer.SetDataModeToBinary()
+      
+      ext = os.path.splitext(fname)[1]
+      if ext == '':
+        name = fname + '.' + writer.GetDefaultFileExtension()
+      elif ext == writer.GetDefaultFileExtension():
+        name = fname
+      else:
+        raise ValueError("unknown extension {}".format(ext))
+      writer.SetFileName(name)
+    else:
+      writer = vtk.vtkDataSetWriter()
+      writer.SetHeader('damask.Geom '+version)
+      writer.WriteToOutputStringOn()
+  
+    writer.SetInputData(rGrid)
+    writer.Write()
+
+    if fname is None: return writer.GetOutputString()
+
                
   def show(self):
     """Show raw content (as in file)"""