From 575da581a9403a0a50e5b5d8d382374ecce02f71 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 11 Mar 2020 22:35:58 +0100
Subject: [PATCH] support for point cloud

---
 python/damask/ktv.py    | 22 +++++++++++++++++-----
 python/damask/result.py | 25 +++++++------------------
 2 files changed, 24 insertions(+), 23 deletions(-)

diff --git a/python/damask/ktv.py b/python/damask/ktv.py
index 5e1919300..67425fd73 100644
--- a/python/damask/ktv.py
+++ b/python/damask/ktv.py
@@ -41,13 +41,12 @@ class VTK: # capitals needed/preferred?
         """
         Create an unstructured grid (mesh).
 
-        connectivity: 1 based at the moment
+        connectivity: 0 based at the moment, shape Ncell x N nodes
         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],
@@ -62,8 +61,21 @@ class VTK: # capitals needed/preferred?
 
 
     @staticmethod
-    def from_points(nodes,connectivity,cell_type):
-        pass
+    def from_polyData(points):
+        vtk_points= vtk.vtkPoints()
+        vtk_points.SetData(nps.numpy_to_vtk(points))
+
+        vertices = vtk.vtkCellArray()
+        vertices.SetNumberOfCells(points.shape[0])
+        T = np.concatenate((np.ones((points.shape[0],1),dtype=np.int64),
+                           np.arange(points.shape[0],dtype=np.int64).reshape(-1,1)),axis=1).ravel()
+        vertices.SetCells(points.shape[0],nps.numpy_to_vtk(T, deep=True, array_type=vtk.VTK_ID_TYPE))
+
+        geom = vtk.vtkPolyData()
+        geom.SetPoints(vtk_points)
+        geom.SetVerts(vertices)
+
+        return VTK(geom)
 
 
     def write(self,fname):                                              #ToDo: Discuss how to handle consistently filename extensions
@@ -95,7 +107,7 @@ class VTK: # capitals needed/preferred?
     def __repr__(self):
         """ASCII representation of the VTK data."""
         writer = vtk.vtkDataSetWriter()
-        writer.SetHeader('DAMASK.VTK v{}'.format(version))
+        writer.SetHeader('# DAMASK.VTK v{}'.format(version))
         writer.WriteToOutputStringOn()
         writer.SetInputData(self.geom)
         writer.Write()
diff --git a/python/damask/result.py b/python/damask/result.py
index 3ac442f5c..9f4c21cf6 100644
--- a/python/damask/result.py
+++ b/python/damask/result.py
@@ -975,11 +975,11 @@ class Result:
             datasets_in = {}
             lock.acquire()
             with h5py.File(self.fname,'r') as f:
-              for arg,label in datasets.items():
-                loc  = f[group+'/'+label]
-                datasets_in[arg]={'data' :loc[()],
-                                  'label':label,
-                                  'meta': {k:v.decode() for k,v in loc.attrs.items()}}
+                for arg,label in datasets.items():
+                    loc  = f[group+'/'+label]
+                    datasets_in[arg]={'data' :loc[()],
+                                      'label':label,
+                                      'meta': {k:v.decode() for k,v in loc.attrs.items()}}
             lock.release()
             r = func(**datasets_in,**args)
             return [group,r]
@@ -1042,26 +1042,15 @@ class Result:
 
           if self.structured:
               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'][()]-1,
                                                 f['/geometry/T_c'].attrs['VTK_TYPE'].decode())
-                  vtk_geom = v.geom
-
         elif mode.lower()=='point':
-          Points   = vtk.vtkPoints()
-          Vertices = vtk.vtkCellArray()
-          for c in self.cell_coordinates():
-            pointID = Points.InsertNextPoint(c)
-            Vertices.InsertNextCell(1)
-            Vertices.InsertCellPoint(pointID)
+            v = VTK.from_polyData(self.cell_coordinates())
 
-          vtk_geom = vtk.vtkPolyData()
-          vtk_geom.SetPoints(Points)
-          vtk_geom.SetVerts(Vertices)
-          vtk_geom.Modified()
+        vtk_geom = v.geom
 
         N_digits = int(np.floor(np.log10(int(self.increments[-1][3:]))))+1