How to compute force resultant using VTK

I'm postprocessing some CFD data using python VTK. After reading the points coordinates and field values from the .vtp I want to integrate the pressure on the body to obtain the 3 components of the force resultant. I'm having troubles in calculating the surface normals and areas for each cell.

This is my code so far:

import vtk
import numpy as np

# Get file_path
file_path = "/path/test.vtp"

# Reading file
if file_path.endswith(".vtk"):
    reader = vtk.vtkDataSetReader() 
elif file_path.endswith(".vtu"):
    reader = vtk.vtkXMLUnstructuredGridReader()
elif file_path.endswith(".vtp"):
    reader = vtk.vtkXMLPolyDataReader()
else:
    print("ERROR: File extension is not supported")
    exit()

reader.SetFileName(file_path)
reader.Update()

# Extracting points coordinates
VTKoutput = reader.GetOutput()
npoints = VTKoutput.GetNumberOfPoints()
coords = np.array([VTKoutput.GetPoint(i) for i in range(npoints)])

# Extracting fields in .vtk
nfields = VTKoutput.GetPointData().GetNumberOfArrays()
field_names = []
for i in range(nfields):
    field_names.append(VTKoutput.GetPointData().GetArrayName(i))

# Extract fields
p = VTKoutput.GetPointData().GetArray('p')
p_array = np.array(p)

I've tried some of the functions that I've found on vtk.org that I thought would get my the normals and the areas, but I don't fully understand the documentation since I'm used to programming in MATLAB. As a last resort I tried doing something like what I did in the first part of the script in order to obtain the vertices of the cells with the idea of manually compute the normals and the area. Something like this:

ncells = VTKoutput.GetNumberOfCells()
F = [0,0,0]
for cell in ncells
    cell_verts = np.array([VTKoutput.GetCell(cell)])
    cell_normal = function_of_verts
    cell_area = function_of_verts
    F = F + p*cell_area*cell_normal

Doing this I also found something I didn't really expect ncells and npoints are not the same, why?

Thanks for helping.

Solution

I am trying to understand your problem, why you expect npoints = ncells in the first place. Cells can be made up of 'n' points.

To calculate cell normal, you can use `vtkPolyDataNormals'

# Create the vtkPolyDataNormals object
normals = vtk.vtkPolyDataNormals()
normals.SetInputData(reader.GetOutput())
normals.ComputeCellNormalsOn()
normals.ComputePointNormalsOff()
normals.Update()

# Access the output of vtkPolyDataNormals
output = normals.GetOutput()
cellNormals = output.GetCellData().GetNormals()

# Print the cell normals
numCells = output.GetNumberOfCells()
for i in range(numCells):
    normal = cellNormals.GetTuple(i)
    print(f"Cell {i}: Normal = {normal}")

Second thing, in the code above where you are trying to calculate Force, I think you should keep the `F = [0,0,0]' inside the for loop.

Edit: Calculate cell area: Based on your cell type, you can calculate the cell area as follows: I am giving example for triangle and quadrilateral cell types, you can extend this for polygonal cell type if required. For polygonal cell, you can triangulate each polygon using any triangulation method (I may prefer Delaunay2D).

for cellId in range(numCells):
    cell = polydata.GetCell(cellId)

    # Check if the cell is a triangle
    if isinstance(cell, vtk.vtkTriangle):
        # Calculate the cell area
        area = cell.TriangleArea()
        print(f"Cell {cellId}: Area = {area}")

    # Check if the cell is a quadrilateral
    elif isinstance(cell, vtk.vtkQuad):
        # Calculate the cell area for quadrilaterals
        pointIds = cell.GetPointIds()
        p0 = polydata.GetPoint(pointIds.GetId(0))
        p1 = polydata.GetPoint(pointIds.GetId(1))
        p2 = polydata.GetPoint(pointIds.GetId(2))
        p3 = polydata.GetPoint(pointIds.GetId(3))

        # Calculate the area by splitting the quad into two triangles
        triangle1 = vtk.vtkTriangle()
        triangle1.GetPointIds().SetId(0, pointIds.GetId(0))
        triangle1.GetPointIds().SetId(1, pointIds.GetId(1))
        triangle1.GetPointIds().SetId(2, pointIds.GetId(2))
        area1 = triangle1.TriangleArea()

        triangle2 = vtk.vtkTriangle()
        triangle2.GetPointIds().SetId(0, pointIds.GetId(0))
        triangle2.GetPointIds().SetId(1, pointIds.GetId(2))
        triangle2.GetPointIds().SetId(2, pointIds.GetId(3))
        area2 = triangle2.TriangleArea()

        area = area1 + area2
        print(f"Cell {cellId}: Area = {area}")