Search code examples
pythonvtkparaviewopenfoam

Python script for plotting the evolution of charts such as Paraview does


I want to code a python script that generates a plot like the one shown to the right in the next screenshot from Paraview:

enter image description here

I have a series of files generated with the command foamToVTK:

enter image description here

Is there in VTK any function similar to PlotOverLine method of Paraview?


Solution

  • Inspired by @Nico Vuaille and https://stackoverflow.com/a/21990296/4286662 answers, I came across another different solution based on probe filters. This solution remains faithfully to the Paraview results. This is the working code:

    from vtk.util import numpy_support as VN
    from matplotlib import pyplot as plt
    import vtk
    import numpy as np
    from os import walk, path, system
    import pandas as pd
    
    def getFilenames():
        ### Get the file names of each step of the simulation
        (dirpath, dirnames, filenames) = next(walk('VTK'))
        ids=[]
        for dir in dirnames:
            ids.append(int(dir.split("_")[1]))
        ids = sorted(ids)
        basename = dirnames[0].split("_")[0]
    
        return ids, basename, dirpath
    
    def createLine(p1, p2):
        # Create the line along which you want to sample
        line = vtk.vtkLineSource()
        line.SetResolution(numPoints)
        line.SetPoint1(p1)
        line.SetPoint2(p2)
        line.Update()
        return line
    
    def probeOverLine(line,reader):
        ### Interpolate the data from the VTK-file on the created line.
        probe = vtk.vtkProbeFilter()
        probe.SetInputConnection(line.GetOutputPort())
        probe.SetSourceConnection(reader.GetOutputPort())
        probe.Update()
        ### Get the velocity from the probe
        return VN.vtk_to_numpy(probe.GetOutput().GetPointData().GetArray('U'))
    
    
    ### Initialization of variables
    cnt=1
    fig= plt.figure()
    numPoints=100
    
    ids, basename, dirpath = getFilenames()
    
    ### Iteration of time steps
    for id in ids[1:]:
        ### Read values from the file of this time step
        filename = "%s/%s_%d/internal.vtu" % (dirpath, basename, id)
        reader = vtk.vtkXMLUnstructuredGridReader()
        reader.SetFileName(filename)
        reader.Update()
    
        if cnt==1:
            ### Get points for a line in the center of the object
            bounds = reader.GetOutput().GetBounds()
            p1 = [(bounds[1] - bounds[0]) / 2.0 + bounds[0], bounds[2], 0]
            p2 = [(bounds[3] - bounds[2]) / 2.0 + bounds[2], bounds[3], 0]
            line = createLine(p1, p2)
    
        plotOverLine = probeOverLine(line, reader)
    
        df = pd.DataFrame(plotOverLine, columns=['X', 'Y', 'Z'])
        plt.clf()
        plt.title('Frame %d' % cnt)
        plt.plot(df)
        plt.legend(df.columns)
        axes = plt.gca()
        plt.draw()
        plt.pause(.05)
    
        cnt+=1
    
    plt.show()
    

    With the following result:

    enter image description here