Point based vtk file in DADF5 class

python/damask/dadf5.py

@@ -846,40 +846,57 @@ class DADF5():
     pool.wait_completion()
 
 
-  def to_vtk(self,labels):
+  def to_vtk(self,labels,mode='Cell'):
     """
-    Export to vtk cell data.
+    Export to vtk cell/point data.
     
     Parameters
     ----------
     labels : list of str
       Labels of the datasets to be exported.
+    mode   : str
+      Export in cell format or point format.
+      Default value is 'Cell'
 
     """
-    if self.structured:
+    if mode=='Cell':
 
-      coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()]
+      if self.structured:
-      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()
+        coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()]
-      grid.SetDimensions(*(self.grid+1))
+        for dim in [0,1,2]:
-      grid.SetXCoordinates(coordArray[0])
+          for c in np.linspace(0,self.size[dim],1+self.grid[dim]):
-      grid.SetYCoordinates(coordArray[1])
+            coordArray[dim].InsertNextValue(c)
-      grid.SetZCoordinates(coordArray[2])
   
+        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!
     else:
+      Points   = vtk.vtkPoints()
+      Vertices = vtk.vtkCellArray()  
+      for c in self.cell_coordinates():
+        pointID = Points.InsertNextPoint(c)
+        Vertices.InsertNextCell(1)
+        Vertices.InsertCellPoint(pointID)
       
-      nodes = vtk.vtkPoints()
+      Polydata = vtk.vtkPolyData()
-      with h5py.File(self.fname) as f:
+      Polydata.SetPoints(Points)
-        nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True))
+      Polydata.SetVerts(Vertices)
-        
+      Polydata.Modified()
-        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
     
@@ -900,7 +917,10 @@ class DADF5():
               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])
+              if mode=='Cell':
+                grid.GetCellData().AddArray(vtk_data[-1])
+              else:
+                Polydata.GetCellData().AddArray(vtk_data[-1])
           else:
             x = self.get_dataset_location(label)
             if len(x) == 0:
@@ -912,7 +932,10 @@ class DADF5():
             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])
+            if mode=='Cell':
+              grid.GetCellData().AddArray(vtk_data[-1])
+            else:
+              Polydata.GetCellData().AddArray(vtk_data[-1])
       self.set_visible('materialpoints',materialpoints_backup)
   
       constituents_backup = self.visible['constituents'].copy()
@@ -929,7 +952,10 @@ class DADF5():
               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])
+              if mode=='Cell':
+                grid.GetCellData().AddArray(vtk_data[-1])
+              else:
+                Polydata.GetCellData().AddArray(vtk_data[-1])
           else:
             x = self.get_dataset_location(label)
             if len(x) == 0:
@@ -939,17 +965,23 @@ class DADF5():
             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])
+            if mode=='Cell':
+              grid.GetCellData().AddArray(vtk_data[-1])
+            else:
+              Polydata.GetCellData().AddArray(vtk_data[-1])
       self.set_visible('constituents',constituents_backup)
      
-      writer = vtk.vtkXMLRectilinearGridWriter() if self.structured else \
+      if mode=='Cell': 
-               vtk.vtkXMLUnstructuredGridWriter()
+        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])
+      else:
+        writer = vtk.vtkXMLPolyDataWriter()
 
-      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),
@@ -958,6 +990,9 @@ class DADF5():
       writer.SetCompressorTypeToZLib()
       writer.SetDataModeToBinary()
       writer.SetFileName(file_out)
-      writer.SetInputData(grid)
+      if mode=='Cell':
+        writer.SetInputData(grid)
+      else:
+        writer.SetInputData(Polydata)
 
       writer.Write()