Write a VTK file from adaptive mesh refinement with Python

Question

My concern is the following one : I have a 3D magnetic field, B, from an adaptive mesh refinement, saved and organized under a HDF format. Now I would like to create a VTK file using a Python routine. My goal is to represent it with Paraview.

Can you help me to create this routine ?

Using :

gridToVTK("./B", np.array(x), np.array(y), np.array(z), pointData = {"B" : B})

I succeed in turning HDF-5 non-uniform data to VTK but adaptive mesh refinement is a bit different in its organisation.

Thanks,

Answer 1

If you only want to visualize your mesh you could use a function similar to the following to create a vtk unstructured grid, in order to load the file in paraview:

def save_vtk_unstructured_grid(path, points, cells, point_data):
    """
    Create a vtu grid file containing quads from a list of points,
    a list of cells and additional point data.
    The list of cells references the points inside the point list via the row index.

    N Points: [[x_0, y_0, z_0],
               [x_1, y_1, z_1],
               ...
               [x_(n-1), y_(n-1), z_(n-1)]]

    M Cells: [[i_00, i_01, i_02, i_03],
              [i_10, i_11, i_12, i_13],
              ...
              [i_(m-1)0, i_(m-1)1, i_(m-1)2, i_(m-1)3]]

    E.g.:
    Cell: p0 x------x p1    =>      Cell indices inside the cell array:
             |      |               [0, 1, 2, 3]
             |      |
             |      |
          p2 x------x p3

    :param path: Save path as string
    :param points: Nx3 numpy array of point coordinates
    :param cells: Mx4 numpy array of point indices for a mesh consisting of quads.
    :param point_data: Nx1 numpy array of containing data at the point coordinates.
    """
    points = vtk.vtkPoints()
    cells = vtk.vtkCellArray()

    # insert points
    for p in points:
        points.InsertNextPoint(p[0], p[1], p[2])

    # insert the quad cells
    for idx in cells:
        # create a new quad cell
        quad = vtk.vtkQuad()
        quad.GetPointIds().SetId(0, idx[0])
        quad.GetPointIds().SetId(1, idx[1])
        quad.GetPointIds().SetId(2, idx[2])
        quad.GetPointIds().SetId(3, idx[3])
        cells.InsertNextCell(quad)

    # create the point data
    data = vtk.vtkDoubleArray()
    data.SetNumberOfComponents(1)
    data.SetNumberOfValues(len(point_data))
    for i, d in enumerate(point_data):
        data.SetTuple(i, d)

    # create the grid from the points and cells
    grid = vtk.vtkUnstructuredGrid()
    grid.SetPoints(points)
    grid.SetCells(vtk.VTK_QUAD, cells)
    grid.GetPointData().SetScalars(data)

    # write the grid
    writer = vtk.vtkXMLUnstructuredGridWriter()
    writer.SetFileName(path)
    writer.SetInputData(grid)
    writer.Write()

This will create a unstructured grid (mesh) consisting out of quads, which are often used for adaptive mesh refinement. All you need to provide is the list of points:

points = np.array([x, y, z])

and the mesh connections as index list cells. Your point data B should be a scalar value. If this has more than one component per point you have to change the number of components for the vtkDoubleArray inside the function.

Please note, that if you input your refined mesh containing hanging nodes, this might affect the connectivity of the mesh and lead to wrong results for algorithms depending on connectivity information. As far as I know vtk does not support hanging node connectivity.