`vtkProbeFilter`: invalid point on a simple 2D sqare mesh with Lagrange cell of order 6
I'm playing around with high-order elements (see e.g. https://www.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/).
In order to check that I'm using the correct ordering for the Lagrange
cell connectivity, I've tried to represent a polynomial of order 6 using a basis of order 6 on a squared mesh with 2 T3
(triangle with 3 nodes) elements. Then, I probe my field over a line and check that the values returned by VTK
match (exactly) the values of my polynomial.
I found a very strange behavior (see code below for reproducing - I'm on VTK 9.2.6
).
If my probe includes the point (0.25, 0.25), this point is "invalid" and
VTK
reports 0.
When opening my unstructured mesh with ParaView
(5.10), the same effect is seen (the line is "broken" at (0.25, 0.25), see below screenshot).
Is this a bug, or something wrong on my side ?
Thanks !
import typing
import numpy
import scipy
import vtk
from vtk.util.numpy_support import vtk_to_numpy
floating_type = typing.Union[float, numpy.ndarray]
def polynomial(order : int, x : floating_type, y : floating_type) -> floating_type:
"""
Use the Pascal triangle to create a 2D polynomial of a given `order`, and evaluate it.
"""
pascal = numpy.fliplr(scipy.sparse.triu(numpy.fliplr(scipy.linalg.pascal(n = order))).todense())
return numpy.polynomial.polynomial.polyval2d(x = x, y = y, c = pascal)
def make_cut(*,
source,
start : typing.Sequence[float],
end : typing.Sequence[float],
resolution : int,
) -> vtk.vtkProbeFilter:
"""
Probe a `source` over a line from `start` to `end`.
"""
# Create the line
line = vtk.vtkLineSource()
line.SetOutputPointsPrecision(vtk.vtkAlgorithm.DOUBLE_PRECISION)
line.SetPoint1(start)
line.SetPoint2(end)
line.SetResolution(resolution)
# Create the probe
probe_filter = vtk.vtkProbeFilter()
probe_filter.SetInputConnection(line.GetOutputPort())
probe_filter.SetSourceData(source)
probe_filter.Update()
# Check for any missing point. See also:
# * https://vtk.org/Wiki/Demystifying_the_vtkProbeFilter
if not probe_filter.GetOutput().GetNumberOfPoints() == probe_filter.GetValidPoints().GetNumberOfTuples():
raise RuntimeError(f"Error(s) occured for the probe {probe_filter}: {probe_filter.GetOutput().GetNumberOfPoints()} != {probe_filter.GetValidPoints().GetNumberOfTuples()} (some points are invalid).")
return probe_filter
def extract_coords_and_values(
dataset : typing.Any,
as_numpy : bool = False,
arrays : typing.Optional[typing.Sequence[str]] = None,
) -> typing.Tuple[
typing.Union[vtk.vtkPoints, numpy.ndarray],
typing.Union[vtk.vtkPointData, typing.Dict[str, typing.Union[numpy.ndarray, vtk.vtkDataArray]]]
]:
"""
Extract point coordinates and values from a dataset.
"""
assert isinstance(dataset, vtk.vtkDataSet)
# Extract points coordinates
coords = dataset.GetPoints()
# Extract values
values = dataset.GetPointData()
# Convert coordinates to NumPy array if requested.
if as_numpy:
coords = vtk_to_numpy(coords.GetData())
# If any of 'arrays' or 'as_numpy' is true, create the list of array names that are needed, then iterate over the names
# and get the arrays, possibly as NumPy arrays.
if arrays or as_numpy:
gen = arrays if arrays is not None else [values.GetArrayName(i) for i in range(values.GetNumberOfArrays())]
values = {
k: values.GetArray(k) if not as_numpy else vtk_to_numpy(values.GetArray(k))
for k in gen
}
return coords, values
def probe_and_check(*,
ugrid,
evaluator : typing.Callable[[floating_type, floating_type], floating_type],
start : typing.Sequence[float] = [0., 0., 0.],
end : typing.Sequence[float] = [1., 1., 0.],
probe_resolution : int = 1000,
):
"""
Check that the `ugrid` and the `evaluator` return the same values over the probing line that goes from `start` to `end`.
"""
probe = make_cut(source = ugrid, start = start, end = end, resolution = probe_resolution)
coords, values = extract_coords_and_values(dataset = probe.GetOutput(), as_numpy = True, arrays = ["my-issue-data"])
expected_values = evaluator(x = coords[:, 0], y = coords[:, 1])
assert coords.shape == (probe_resolution + 1, 3)
assert values["my-issue-data"].shape == (probe_resolution + 1,)
assert expected_values.shape == (probe_resolution + 1,)
assert numpy.allclose(a = values["my-issue-data"], b = expected_values, rtol = 1.e-14, atol = 1.e-14, equal_nan = False)
if __name__ == "__main__":
print("> VTK version:", vtk.__version__)
# Order of the polynomial that should be exactly represented by the tested basis order.
POLYNOMIAL_ORDER = 6
# Initialize unstructed grid.
ugrid = vtk.vtkUnstructuredGrid()
# Points of the grid are the "unique dofs locations".
# The domain is a mesh with 2 triangle elements, with a basis of order 6.
# The first element has 28 unique dofs, the second has 21 unique dofs (and shares dofs on the diagonal with the first element).
points = vtk.vtkPoints()
points.SetDataTypeToDouble()
points.SetNumberOfPoints(28 + 21)
h = 1. / 6.
points_raw = {
# Points of the first cell.
0 : (0. , 0. , 0.), # vertex 0
1 : (6. * h, 0. , 0.), # vertex 1
2 : (0. , 6. * h, 0.), # vertex 2
3 : (1. * h, 0. , 0.), # edge 0
4 : (2. * h, 0. , 0.), # edge 0
5 : (3. * h, 0. , 0.), # ...
6 : (4. * h, 0. , 0.),
7 : (5. * h, 0. , 0.),
8 : (5. * h, 1. * h, 0.),
9 : (4. * h, 2. * h, 0.),
10 : (3. * h, 3. * h, 0.),
11 : (2. * h, 4. * h, 0.),
12 : (1. * h, 5. * h, 0.),
13 : (0. , 5. * h, 0.),
14 : (0. , 4. * h, 0.),
15 : (0. , 3. * h, 0.),
16 : (0. , 2. * h, 0.),
17 : (0. , 1. * h, 0.),
18 : (1. * h, 1. * h, 0.),
19 : (2. * h, 1. * h, 0.),
20 : (3. * h, 1. * h, 0.),
21 : (4. * h, 1. * h, 0.),
22 : (3. * h, 2. * h, 0.),
23 : (2. * h, 3. * h, 0.),
24 : (1. * h, 4. * h, 0.),
25 : (1. * h, 3. * h, 0.),
26 : (1. * h, 2. * h, 0.),
27 : (2. * h, 2. * h, 0.),
# Points of the second cell.
28 : (6. * h, 6. * h, 0.),
29 : (6. * h, 1. * h, 0.),
30 : (6. * h, 2. * h, 0.),
31 : (6. * h, 3. * h, 0.),
32 : (6. * h, 4. * h, 0.),
33 : (6. * h, 5. * h, 0.),
34 : (5. * h, 6. * h, 0.),
35 : (4. * h, 6. * h, 0.),
36 : (3. * h, 6. * h, 0.),
37 : (2. * h, 6. * h, 0.),
38 : (1. * h, 6. * h, 0.),
39 : (5. * h, 5. * h, 0.),
40 : (5. * h, 4. * h, 0.),
41 : (5. * h, 3. * h, 0.),
42 : (5. * h, 2. * h, 0.),
43 : (4. * h, 3. * h, 0.),
44 : (4. * h, 4. * h, 0.),
45 : (4. * h, 5. * h, 0.),
46 : (3. * h, 5. * h, 0.),
47 : (3. * h, 4. * h, 0.),
48 : (2. * h, 5. * h, 0.),
}
# Ensure that all the points are unique.
assert len(points_raw) == len(set(points_raw.values()))
# Set grid points.
for k, v in points_raw.items():
points.SetPoint(k, v)
ugrid.SetPoints(points)
# Cell array placeholder.
cells = vtk.vtkCellArray()
cells.SetNumberOfCells(2)
# Cell of order 6, connectivity ordered according to
# https://www.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/..
cells.InsertNextCell(28, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 21, 24, 19, 20, 22, 23, 25, 26, 27])
cells.InsertNextCell(28, [1, 28, 2, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 12, 11, 10, 9, 8, 42, 39, 48, 41, 40, 45, 46, 47, 43, 44])
# Set grid cells
ugrid.SetCells(vtk.VTK_LAGRANGE_TRIANGLE, cells)
# Evaluator.
evaluator = lambda x, y: polynomial(order = POLYNOMIAL_ORDER + 1, x = x, y = y)
# Set grid values.
data = vtk.vtkDoubleArray()
data.SetName('my-issue-data')
data.SetNumberOfValues(points.GetNumberOfPoints())
for i in range(points.GetNumberOfPoints()):
x, y, _ = points.GetPoint(i)
value = evaluator(x = x, y = y)
data.SetValue(i, value)
ugrid.GetPointData().AddArray(data)
# Perform checks
assert ugrid.GetNumberOfCells() == 2
probe_and_check(ugrid = ugrid, evaluator = evaluator, probe_resolution = 101)
# Fails ! The point (0.25,0.25) is invalid. Note that this point is exactly at the center of the first cell.
probe_and_check(ugrid = ugrid, evaluator = evaluator, probe_resolution = 100)