vtkCellSizeFilter() generates non plausible values for cell types 13 and 41
Describe the bug, what's wrong, and what you expected.
I created an unstructured grid with different types of prisms. Then I calculated the cell volumes with vtkCellSizeFilter().
For cell type vtk.VTK_WEDGE (13) I find negativ volumes.
For cell type vtk.VTK_CONVEX_POINT_SET (41) the volumes do not seem to be plausible when there are more than one cell in the grid.
Summary of reproduction
- When the prism is of type 13, the volume is negative. When the same prism is of type 41, the volume is positive and as expected.
- The given order of points results in positive volumes for all types except for type 13. (-> I assume, the calculation of the volumes is not consistent for different cell types)
- The volume of prisms of type 41 is as expected, when this prism is the first one in the grid. (-> I assume, the algorithm to calculate the volume is correct)
- The volume of prisms of type 41 are not as expected, when there is more than one prism in the grid and the prism is not the first element in the grid. (-> I assume, the algorithm picks wrong points to calculate the volume)
Steps to reproduce the bug.
Here is a volume comparision for prisms with different number of edges.
import numpy as np
import matplotlib.pyplot as plt
import vtk
import pyvista as pv
def area_polygon(p):
"""
Calculates the area on a polygon given by an array of points.
https://en.wikipedia.org/wiki/Shoelace_formula#Generalization
"""
npoints = np.shape(p)[0]
A = 0
for i in np.arange(npoints):
A = A + np.cross(p[i, :], p[(i+1) % npoints, :])
return np.linalg.norm(A)/2
def area(n=100, radius=1):
"""
Calculates the area of a polygon with n edges given by polygon_edges().
For radius 1 the value tends towards pi.
"""
p = polygon_edges(n, radius)
return area_polygon(p)
def polygon_edges(n, radius=1, plot=False):
"""
Returns n evenly distributed 3D points on a circle with the given radius
"""
i = np.linspace(0, 2*np.pi, n+1)
i = np.delete(i, n)
p = np.column_stack((np.sin(i)*radius, np.cos(i)*radius, np.ones(n)))
if plot:
plt.plot(p[:, 0], p[:, 1], 'o')
plt.show()
return p
def UnstructuredGrid(coords, ind, cell_types):
# Erstelle ein VTK-UnstructuredGrid-Objekt
grid = vtk.vtkUnstructuredGrid()
# Erstelle Punkte
points = vtk.vtkPoints()
for coord in coords:
points.InsertNextPoint(coord)
grid.SetPoints(points)
# Erstelle Zellen
cell_array = vtk.vtkCellArray()
for i, cell_type_id in enumerate(cell_types):
cell = vtk.vtkIdList()
for j in range(len(ind[i])):
cell.InsertNextId(ind[i][j])
cell_array.InsertNextCell(cell)
grid.SetCells(cell_types, cell_array)
return grid
def get_element_volumes(mesh):
# Create a vtkCellSizeFilter
sizeFilter = vtk.vtkCellSizeFilter()
sizeFilter.SetInputData(mesh)
sizeFilter.Update()
# Get the output data set
output = sizeFilter.GetOutput()
# Get the volume array
volumes = output.GetCellData().GetArray("Volume")
# Convert the volume array to a list
volumes_list = [volumes.GetValue(i) for i in range(volumes.GetNumberOfTuples())]
return volumes_list
# When there are many prism in the grid, cell type 41 gives wrong volumes
edge_min, edge_max = 3, 9
# When there is only one cell in the grid, cell type 41 gives correct volumes
#edge_min, edge_max = 7, 7
test = 1 # Check cell types from ctype (see below)
# test = 2 # Check 41 -> vtk.VTK_CONVEX_POINT_SET only
shift_factor = 3
edges = np.arange(edge_min, edge_max+1)
# Type of geometry
ctype = {3: vtk.VTK_WEDGE, # -> 13
4: vtk.VTK_HEXAHEDRON, # -> 12
5: vtk.VTK_PENTAGONAL_PRISM, # -> 15
6: vtk.VTK_HEXAGONAL_PRISM, # -> 16
7: vtk.VTK_CONVEX_POINT_SET} # -> 41
# Init arrays
ind = list()
coords = np.empty((0, 3))
cell_type = list()
index = 0
# Create prism with different number of edges for unstructured grid
for n in edges:
p = polygon_edges(n)
p[:, 0] = p[:, 0] + shift_factor * n
p = np.vstack([p, p])
p[n:, 2] = 0
coords = np.vstack((coords, p))
#ind.append(2*n)
#ind = ind + list(range(index, index+2*n))
ind.append(list(range(index, index+2*n)))
index = index + 2*n
if test == 1:
cell_type.append(ctype.get(n, vtk.VTK_CONVEX_POINT_SET))
else:
cell_type.append(vtk.VTK_CONVEX_POINT_SET)
grid = UnstructuredGrid(coords, ind, cell_type)
pv.UnstructuredGrid(grid).plot(show_edges=True)
# Calculate actual volume of cells
volumes = get_element_volumes(grid)
# Calculate reference values for volumes of cells
# Height of cells is 1 -> volume = area * 1
# Reference values are plausible, because with increasing number of edges n
# the value tends towards pi for edges on a circle with radius 1.
reference_volume = [area(n) for n in np.arange(edge_min, edge_max+1)]
# Compare actual and reference volumes
print("Cell type | Actual value | Reference value -> match?")
for [vol_ac, vol_des, ct] in zip(volumes, reference_volume, cell_type):
match = ":)" if abs(vol_ac - vol_des) < 1e-4 else "X"
print(f'{ct} | {vol_ac: .8f} | {vol_des: .8f} -> {match}')
- This test gives for "test = 1"
Cell type | Actual value | Reference value -> match?
13 | -1.29903746 | 1.29903811 -> X
12 | 2.00000000 | 2.00000000 -> :)
15 | 2.37764177 | 2.37764129 -> :)
16 | 2.59807777 | 2.59807621 -> :)
41 | -1.04465834 | 2.73641019 -> X
41 | 0.76634582 | 2.82842712 -> X
41 | 0.66068331 | 2.89254424 -> X
Here you can see, that the volume calculated for prism 13 is negativ, but needs to be positiv. The volumes calculated for prisms of type 41 is always not matching the reference values. This grid has seven cells.
- This test gives for "test = 2" and "edge_min, edge_max = 3, 3":
Cell type | Actual value | Reference value -> match?
41 | 1.29903746 | 1.29903811 -> :)
Here is only one cell in the grid and the volume for a prism with 3 edges is calculated as expected.
- And finally this test gives for "test = 2" and "edge_min, edge_max = 7, 7"
Cell type | Actual value | Reference value -> match?
41 | 2.73641041 | 2.73641019 -> :)
Again, here is only one cell in the grid and the volume for a prism with 7 edges is calculated as expected.