From 744e3bb50bc309e53c59c62ce9ebc813c98d4270 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 11 Mar 2020 20:56:58 +0100
Subject: [PATCH] vectorized cell assignment + bugfix for writing out

---
 python/damask/ktv.py    | 51 ++++++++++++++++++++++-------------------
 python/damask/result.py |  4 ++--
 2 files changed, 30 insertions(+), 25 deletions(-)

diff --git a/python/damask/ktv.py b/python/damask/ktv.py
index 7afbffcc6..5e1919300 100644
--- a/python/damask/ktv.py
+++ b/python/damask/ktv.py
@@ -3,7 +3,7 @@ import os
 import pandas as pd
 import numpy as np
 import vtk
-from vtk.util import numpy_support
+from vtk.util import numpy_support as nps
 
 from . import table
 from . import version
@@ -24,11 +24,11 @@ class VTK: # capitals needed/preferred?
         """Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data."""
         coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()]
         for dim in [0,1,2]:
-            for c in np.linspace(0,size[dim],1+grid[dim]):
+            for c in np.linspace(origin[dim],origin[dim]+size[dim],grid[dim]):
                 coordArray[dim].InsertNextValue(c)
 
         geom = vtk.vtkRectilinearGrid()
-        geom.SetDimensions(*(grid+1))
+        geom.SetDimensions(*grid)
         geom.SetXCoordinates(coordArray[0])
         geom.SetYCoordinates(coordArray[1])
         geom.SetZCoordinates(coordArray[2])
@@ -37,35 +37,40 @@ class VTK: # capitals needed/preferred?
 
 
     @staticmethod
-    def from_unstructuredGrid(nodes,connectivity,elem):
+    def from_unstructuredGrid(nodes,connectivity,cell_type):
+        """
+        Create an unstructured grid (mesh).
+
+        connectivity: 1 based at the moment
+        cell_type: TRIANGLE, 'QUAD', 'TETRA','HEXAHEDRON'
+
+        """
+        vtk_nodes = vtk.vtkPoints()
+        vtk_nodes.SetData(nps.numpy_to_vtk(nodes))
+
+        cells = vtk.vtkCellArray()
+        cells.SetNumberOfCells(connectivity.shape[0])
+        T = np.concatenate((np.ones((connectivity.shape[0],1),dtype=np.int64)*connectivity.shape[1],
+                           connectivity),axis=1).ravel()
+        cells.SetCells(connectivity.shape[0],nps.numpy_to_vtk(T, deep=True, array_type=vtk.VTK_ID_TYPE))
+
         geom = vtk.vtkUnstructuredGrid()
-        geom.SetPoints(numpy_support.numpy_to_vtk(nodes)) #,deep=True)
-        geom.Allocate(connectivity.shape[0])
-
-        if   elem == 'TRIANGLE':
-            vtk_type = vtk.VTK_TRIANGLE
-            n_nodes = 3
-        elif elem == 'QUAD':
-            vtk_type = vtk.VTK_QUAD
-            n_nodes = 4
-        elif elem == 'TETRA':
-            vtk_type = vtk.VTK_TETRA
-            n_nodes = 4
-        elif elem == 'HEXAHEDRON':
-            vtk_type = vtk.VTK_HEXAHEDRON
-            n_nodes = 8
-
-        for i in connectivity:
-            geom.InsertNextCell(vtk_type,n_nodes,i-1)
+        geom.SetPoints(vtk_nodes)
+        geom.SetCells(eval('vtk.VTK_{}'.format(cell_type.upper())),cells)
 
         return VTK(geom)
 
 
+    @staticmethod
+    def from_points(nodes,connectivity,cell_type):
+        pass
+
+
     def write(self,fname):                                              #ToDo: Discuss how to handle consistently filename extensions
         if  (isinstance(self.geom,vtk.vtkRectilinearGrid)):
             writer = vtk.vtkXMLRectilinearGridWriter()
         elif(isinstance(self.geom,vtk.vtkUnstructuredGrid)):
-            writer = vtk.vtkUnstructuredGrid()
+            writer = vtk.vtkXMLUnstructuredGridWriter()
         elif(isinstance(self.geom,vtk.vtkPolyData)):
             writer = vtk.vtkXMLPolyDataWriter()
 
diff --git a/python/damask/result.py b/python/damask/result.py
index 3c1cb15a1..3ac442f5c 100644
--- a/python/damask/result.py
+++ b/python/damask/result.py
@@ -1041,12 +1041,12 @@ class Result:
         if mode.lower()=='cell':
 
           if self.structured:
-              v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin)
+              v = VTK.from_rectilinearGrid(self.grid+1,self.size,self.origin)
               vtk_geom = v.geom
           else:
               with h5py.File(self.fname,'r') as f:
                   v = VTK.from_unstructuredGrid(f['/geometry/x_n'][()],
-                                                f['/geometry/T_c'],
+                                                f['/geometry/T_c'][()]-1,
                                                 f['/geometry/T_c'].attrs['VTK_TYPE'].decode())
                   vtk_geom = v.geom