extract-plane-from-fluent-cff-file

#!/usr/bin/env python
import vtk
import numpy as np
import matplotlib.pyplot as plt

r = vtk.vtkFLUENTCFFReader()
r.SetFileName("/tmp/ansys/mesh_3ddp.cas.h5")
r.Update()

# The cut plane
plane = vtk.vtkPlane()
plane.SetOrigin(0, 0, 0)
plane.SetNormal(0, 0, 1)

cut0 = vtk.vtkCutter()
cut0.SetCutFunction(plane)
cut0.SetInputConnection(r.GetOutputPort())
cut0.Update()

data = cut0.GetOutputDataObject(0).GetBlock(1)

triangles = data.GetPolys().GetData()
points = data.GetPoints()

mapper = vtk.vtkCellDataToPointData()
mapper.AddInputData(data)
mapper.Update()

vels = mapper.GetOutput().GetPointData().GetArray(1)

ntri = int(triangles.GetNumberOfTuples()/4)
npts = int(points.GetNumberOfPoints())
nvls = int(vels.GetNumberOfTuples())

tri = np.zeros((ntri,3))
x = np.zeros(npts)
y = np.zeros(npts)
ux = np.zeros(nvls)
uy = np.zeros(nvls)

for i in range(0, ntri):
    tri[i, 0] = triangles.GetTuple(4*i + 1)[0]
    tri[i, 1] = triangles.GetTuple(4*i + 2)[0]
    tri[i, 2] = triangles.GetTuple(4*i + 3)[0]

for i in range(npts):
    pt = points.GetPoint(i)
    x[i] = pt[0]
    y[i] = pt[1]

for i in range(0, nvls):
    U = vels.GetTuple(i)
    ux[i] = U[0]
#    uy[i] = U[1]

# Mesh
plt.figure(figsize=(8,8))
plt.triplot(x, y, tri)
plt.gca().set_aspect('equal')
plt.savefig('fit1.png')

# Velocity x-component
plt.figure(figsize=(8,8))
plt.tricontourf(x, y, tri, ux, 16)
plt.tricontour(x, y, tri, ux, 16)
plt.savefig('fit2.png')