From b14c15fd9ee6720d79677ed30d4558e3928269a5 Mon Sep 17 00:00:00 2001
From: Martin Diehl <m.diehl@mpie.de>
Date: Fri, 13 Dec 2019 12:15:45 +0100
Subject: [PATCH] directly output DADF5 to vtk from python

---
 python/damask/dadf5.py | 120 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 120 insertions(+)

diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py
index beced188d..1142d02b4 100644
--- a/python/damask/dadf5.py
+++ b/python/damask/dadf5.py
@@ -1,7 +1,10 @@
 from queue import Queue
 import re
 import glob
+import os
 
+import vtk
+from vtk.util import numpy_support
 import h5py
 import numpy as np
 
@@ -841,3 +844,120 @@ class DADF5():
         N_added +=1
 
     pool.wait_completion()
+
+
+  def to_vtk(self,labels):
+    """
+    Export to vtk cell data.
+    
+    Parameters
+    ----------
+    labels : list of str
+      Labels of the datasets to be exported.
+
+    """
+    if self.structured:
+
+      coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()]
+      for dim in [0,1,2]:
+        for c in np.linspace(0,self.size[dim],1+self.grid[dim]):
+          coordArray[dim].InsertNextValue(c)
+
+      grid = vtk.vtkRectilinearGrid()
+      grid.SetDimensions(*(self.grid+1))
+      grid.SetXCoordinates(coordArray[0])
+      grid.SetYCoordinates(coordArray[1])
+      grid.SetZCoordinates(coordArray[2])
+
+    else:
+      
+      nodes = vtk.vtkPoints()
+      with h5py.File(self.fname) as f:
+        nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True))
+        
+        grid = vtk.vtkUnstructuredGrid()
+        grid.SetPoints(nodes)
+        grid.Allocate(f['/geometry/T_c'].shape[0])
+        for i in f['/geometry/T_c']:
+          grid.InsertNextCell(vtk.VTK_HEXAHEDRON,8,i-1) # not for all elements!
+
+    N_digits = int(np.floor(np.log10(int(self.increments[-1][3:]))))+1
+    
+    for i,inc in enumerate(self.iter_visible('increments')):
+      vtk_data = []
+      
+      materialpoints_backup = self.visible['materialpoints'].copy()
+      self.set_visible('materialpoints',False)
+      for label in labels:
+        for p in self.iter_visible('con_physics'):
+          if p != 'generic':
+            for c in self.iter_visible('constituents'):
+              x = self.get_dataset_location(label)
+              if len(x) == 0:
+                continue
+              array = self.read_dataset(x,0)
+              shape = [array.shape[0],np.product(array.shape[1:])]
+              vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),
+                              deep=True,array_type= vtk.VTK_DOUBLE))
+              vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: hard coded 1!
+              grid.GetCellData().AddArray(vtk_data[-1])
+          else:
+            x = self.get_dataset_location(label)
+            if len(x) == 0:
+              continue
+            array = self.read_dataset(x,0)
+            shape = [array.shape[0],np.product(array.shape[1:])]
+            vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),
+                            deep=True,array_type= vtk.VTK_DOUBLE))
+            ph_name   = re.compile(r'(\/[1-9])_([A-Z][a-z]*)_(([a-z]*)|([A-Z]*))')  #looking for phase name in dataset name
+            dset_name = '1_' + re.sub(ph_name,r'',x[0].split('/',1)[1])             #removing phase name from generic dataset
+            vtk_data[-1].SetName(dset_name)
+            grid.GetCellData().AddArray(vtk_data[-1])
+      self.set_visible('materialpoints',materialpoints_backup)
+
+      constituents_backup = self.visible['constituents'].copy()
+      self.set_visible('constituents',False)
+      for label in labels:
+        for p in self.iter_visible('mat_physics'):
+          if p != 'generic':
+            for m in self.iter_visible('materialpoints'):
+              x = self.get_dataset_location(label)
+              if len(x) == 0:
+                continue
+              array = self.read_dataset(x,0)
+              shape = [array.shape[0],np.product(array.shape[1:])]
+              vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),
+                              deep=True,array_type= vtk.VTK_DOUBLE))
+              vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: why 1_?
+              grid.GetCellData().AddArray(vtk_data[-1])
+          else:
+            x = self.get_dataset_location(label)
+            if len(x) == 0:
+              continue
+            array = self.read_dataset(x,0)
+            shape = [array.shape[0],np.product(array.shape[1:])]
+            vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),
+                            deep=True,array_type= vtk.VTK_DOUBLE))
+            vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
+            grid.GetCellData().AddArray(vtk_data[-1])
+      self.set_visible('constituents',constituents_backup)
+      
+      writer = vtk.vtkXMLRectilinearGridWriter() if self.structured else \
+               vtk.vtkXMLUnstructuredGridWriter()
+
+      x = self.get_dataset_location('u_n')
+      vtk_data.append(numpy_support.numpy_to_vtk(num_array=self.read_dataset(x,0),
+                      deep=True,array_type=vtk.VTK_DOUBLE))
+      vtk_data[-1].SetName('u')
+      grid.GetPointData().AddArray(vtk_data[-1])
+      
+      file_out = '{}_inc{}.{}'.format(os.path.splitext(os.path.basename(self.fname))[0],
+                                      inc[3:].zfill(N_digits),
+                                      writer.GetDefaultFileExtension())
+      
+      writer.SetCompressorTypeToZLib()
+      writer.SetDataModeToBinary()
+      writer.SetFileName(file_out)
+      writer.SetInputData(grid)
+      
+      writer.Write()