VTK writer improvements

[BinWang0213]

Is your feature request related to a problem? Please describe.
The vtk writer for PolyMesh was actually a facets writer instead of cell wrtier. It also doesn't include the attributes such as volume, phase number and seed numbers.

I have wrote a vtk writer for polymesh. Now you can do following analyze:

This is a example of checking node ind and cell ind for any particular polyhedron. It will be very useful for debuging purpose.

Describe the solution you'd like
The working code is attached, please feel free to integrate it into the MicroStructPy.

seed_basename = 'seeds.vtk'
seed_filename = os.path.join(directory, seed_basename)
#seeds2VTK(seed_filename,seeds)
seeds.write(seed_filename)

poly_basename = 'polymesh.vtk'
poly_filename = os.path.join(directory, poly_basename)
#polymesh2VTK(poly_filename,pmesh)
pmesh.write(poly_filename)

import vtk
import vtk.util.numpy_support as vtk_np

def seeds2VTK(fname, seeds):
    numPts=len(seeds)
    geom_name = type(seeds.seeds[0].geometry).__name__.lower().strip()
    centers=np.array([s.position for s in seeds.seeds])
    
    if(centers.shape[1]==2):#extend to 3D pts
        temp = np.zeros((numPts,3))
        temp[:,:-1] = centers
        centers=temp
    
    vtkPoly=vtk.vtkPolyData()
    if(geom_name in ['circle','sphere']):
        verts = vtk.vtkPoints()
        cells = vtk.vtkCellArray()
        
        verts.SetNumberOfPoints(numPts)
        verts.SetData(vtk_np.numpy_to_vtk(centers))
        vtkPoly.SetPoints(verts)

        cells_npy = np.vstack([np.ones(numPts,dtype=np.int64),
                               np.arange(numPts,dtype=np.int64)]).T.flatten()
        cells.SetCells(numPts,vtk_np.numpy_to_vtkIdTypeArray(cells_npy))
        vtkPoly.SetVerts(cells)

        diameters=np.array([s.geometry.diameter for s in seeds.seeds])
        #print(diameters)
        diameters = vtk_np.numpy_to_vtk(diameters)
        diameters.SetName("Diameters")
        vtkPoly.GetCellData().SetScalars(diameters)
        
        
        MaterialID=np.array([s.phase for s in seeds.seeds])
        MaterialID = vtk_np.numpy_to_vtk(MaterialID)
        MaterialID.SetName("MaterialID")
        vtkPoly.GetCellData().AddArray(MaterialID)

        writer = vtk.vtkPolyDataWriter()
        writer.SetFileName(fname)
        writer.SetInputData(vtkPoly)
        writer.SetFileTypeToASCII() #We can switch to binary for large dataset
        writer.Write() 
    else:
        print('Other shape is not supported to write to VTK yet!')

def polymesh2VTK(fname,pmesh):
    ugrid = vtk.vtkUnstructuredGrid()

    # Create real polyhedron
    points = vtk.vtkPoints()

    points.SetNumberOfPoints(len(pmesh.points))
    pmesh.points=np.array(pmesh.points)
    for i in range(len(pmesh.points)):
        points.SetPoint(i, pmesh.points[i,0],pmesh.points[i,1],pmesh.points[i,2])

    ugrid.SetPoints(points)

    # List of pointIds that make up the face
    for cell in pmesh.regions:
        faceId = vtk.vtkIdList()
        faceId.InsertNextId(len(cell)) # Number faces that make up the cell.
        for face in cell:
            face_pts=pmesh.facets[face]
            faceId.InsertNextId(len(face_pts)) # Number of points in face
            [faceId.InsertNextId(i) for i in face_pts] # Insert the pointIds for the face
        ugrid.InsertNextCell(vtk.VTK_POLYHEDRON, faceId)

    MaterialID = vtk_np.numpy_to_vtk(num_array=np.array(pmesh.phase_numbers))
    MaterialID.SetName("MaterialID")
    ugrid.GetCellData().SetScalars(MaterialID)
    
    SeedsID = vtk_np.numpy_to_vtk(num_array=np.array(pmesh.seed_numbers))
    SeedsID.SetName("SeedsID")
    ugrid.GetCellData().AddArray(SeedsID)

    Vol=vtk_np.numpy_to_vtk(num_array=np.array(pmesh.volumes)) 
    Vol.SetName("Volume")
    ugrid.GetCellData().AddArray(Vol)

    #writer = vtk.vtkXMLUnstructuredGridWriter()
    writer = vtk.vtkUnstructuredGridWriter()
    writer.SetInputData(ugrid)
    writer.SetFileName(fname)
    writer.SetFileTypeToASCII()
    writer.Write()

[kip-hart]

Hi Bin,

Thanks so much for this. I was unfamiliar with the vtk package at the time I wrote file I/O code. I'll take a look at adding this in - I definitely don't have a problem with adding functionality!

Best,
Kip

[kip-hart]

Hi Bin,

SeedList and PolyMesh now output proper VTK files. I ended up not using the vtk package and instead followed the format of the output files. Calling seeds.write(seed_filename, format='vtk') will write the seeds to VTK, and similarly pmesh.write(poly_filename, format='vtk'). Feel free to re-open if there are any problems. Once I add gmsh support, I'll create a new release with these features. Thanks for the feature recommendation!

Best,
Kip