Rendering cutting planes created with OpenFOAM

Support
1

Hello everyone,

I am working on a Python script to process an OpenFOAM simulation automatically. I’d like to do in Python with vtk, because the rest of my post processing is also in Python. My goal is to render all the cutting planes I created in OpenFOAM in the vtk-format (shown in the following section) and to safe images of these planes in png-format or something similar. Currently my problem is that I don’t get the geometry coloured with the field data. The only thing I get is a grey representation of the planes geometry.
I am using version ‘8.2.0’ of vtk for Python.

The file I’d like to process looks like this:
#vtk DataFile Version 2.0
sampleSurface
ASCII
DATASET POLYDATA
POINTS 3540 float
-0.0314316 0.00495 0.000881417 -0.031924 0.00495 0.000968241 -0.031924 0.00495 0.00102257 -0.0314316 0.00495
#many lines following with coordinates
0.02412 0.00495 -7.30676e-05 0.02484 0.00495 -0.0001 0.02556 0.00495 -0.0001
POLYGONS 3382 16910
4 0 1 2 3 4 4 0 3 5 4
#many lines following
POINT_DATA 3540
FIELD attributes 1
H2O 1 3540 float
0.042 0.042 0.042 0.042 0.042 0.042 0.042 0.042 0.042 0.042 0.042
#many lines with data following

My rendering code looks like this:

pathToPlane =‘pathToData’
#print(pathToPlane)

reader = vtk.vtkDataSetReader()
reader.SetFileName(pathToPlane)
reader.ReadAllScalarsOn()
reader.ReadAllColorScalarsOn()
reader.ReadAllNormalsOn()
reader.ReadAllTCoordsOn()
reader.ReadAllVectorsOn()
reader.Update() # Needed because of GetScalarRange
polydata = reader.GetOutput()

#mapper
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(reader.GetOutputPort())
mapper.SetScalarModeToUseCellFieldData()
mapper.SelectColorArray(“H2O”)
mapper.Update()

actor = vtk.vtkActor()
actor.SetMapper(mapper)

#A renderer and render window
renderer = vtk.vtkRenderer()

#render window
renwin = vtk.vtkRenderWindow()
renwin.SetSize(800,600)
renwin.AddRenderer(renderer)
renwin.Render()

#add the actor
renderer.AddActor(actor)

#An interactor
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(renwin)

#Start
interactor.Initialize()
interactor.Start()

Hope you can help me. Thanks in advance!

2

Could take a look at the in-situ processing (runTimePostProcessing)

1 Like
3

Hi, thanks for the link. Looks quite interesting.
I managed to find my error: My cutting plane only contains point data. Therefore the mapper must be set to mapper.SetScalarModeToUsePointFieldData()

The whole script with export of the screenshot looks like:
pathToPlane = casePath + ‘postProcessing/surfaces/’ + time + ‘/’ + fieldName + ‘_’ + cuttingPlaneName + ‘.vtk’
#print(pathToPlane)

reader = vtk.vtkDataSetReader()
reader.SetFileName(pathToPlane)
reader.ReadAllScalarsOn()
reader.ReadAllColorScalarsOn()
reader.ReadAllNormalsOn()
reader.ReadAllTCoordsOn()
reader.ReadAllVectorsOn()
reader.Update() # Needed because of GetScalarRange
polydata = reader.GetOutput()   
scalar_range = polydata.GetScalarRange()

points = polydata.GetPoints()
array = points.GetData()
numpy_coordinates = vtk_to_numpy(array)

fieldData_vtk_array = reader.GetOutput().GetPointData().GetArray(0)
fieldData_np_array = np.array(vtk_to_numpy(fieldData_vtk_array))

x_y_Coordinates = numpy_coordinates[:,[0,2]]

z_Coordinate = numpy_coordinates[0,2]

#column_names = ['x', 'y', 'z', 'H2O']
column_names = ['x', 'z', fieldName] #information only

dataArray = np.ones([len(fieldData_np_array), len(column_names)])
dataArray[:,0:2] = x_y_Coordinates
dataArray[:,2] = fieldData_np_array


zMin = np.min(fieldData_np_array)
zMax = np.max(fieldData_np_array)


#mapper
mapper = vtk.vtkPolyDataMapper()

mapper.SetInputConnection(reader.GetOutputPort())
mapper.SetDebug(True)
mapper.SetArrayName('H2O')
mapper.SetScalarRange(zMin, zMax)
mapper.ScalarVisibilityOn()
mapper.SetInterpolateScalarsBeforeMapping(1)
mapper.SetScalarModeToUsePointFieldData()
mapper.SelectColorArray("H2O")
mapper.Update()

actor = vtk.vtkActor()
actor.SetMapper(mapper)

#A renderer and render window
renderer = vtk.vtkRenderer()
#add the actor
renderer.AddActor(actor)


camera = renderer.GetActiveCamera()
camera.SetPosition([0, 0.1, 0])
camera.SetFocalPoint([0.0, 0.027000000700354576, 0.000881417])
camera.SetViewUp([0, 0, -1])
renderer.SetActiveCamera(camera)

#render window
renwin = vtk.vtkRenderWindow()
renwin.SetSize(800,600)
renwin.AddRenderer(renderer)
renwin.Render()

#An interactor
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(renwin)

w2if = vtk.vtkWindowToImageFilter()
w2if.SetInput(renwin)
w2if.Update()

writer = vtk.vtkPNGWriter()
writer.SetFileName("screenshot.png")
writer.SetInputConnection(w2if.GetOutputPort())
writer.Write()

#Start
interactor.Initialize()
interactor.Start()