🐛[BUG]: XAeronet-V Preprocessor.py unstructured2voxel function

[LorenzoMattia]

Version

1.0.1

On which installation method(s) does this occur?

Pip

Describe the issue

I was trying to run the XAeronet Volume example with the AhmedML Dataset and having some troubles in the Preprocessor.py script. Theorically the dataset should be structured in the same way as the DrivAerML that is suggested in the repo.

The function unstrctured2voxel takes as input bounds variable containing a tuple of three elements. However, according to the VTK library the function SetSamplingBounds of the vtk.vtkResampleToImage() takes as input a tuple of six elements representing the min and max of the bound of each dimension. Indeed, in the following is the error I am getting.

Error executing job with overrides: []
multiprocessing.pool.RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/usr/lib64/python3.12/multiprocessing/pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
                    ^^^^^^^^^^^^^^^^^^^
  File "/physicsnemo/examples/cfd/external_aerodynamics/xaeronet/volume/preprocessor.py", line 147, in process_file
    voxel_grid = unstructured2voxel(vtu_mesh, grid_size_expanded, bounds)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/physicsnemo/examples/cfd/external_aerodynamics/xaeronet/volume/preprocessor.py", line 62, in unstructured2voxel
    resampler.SetSamplingBounds(bounds)
TypeError: SetSamplingBounds argument 1

In order to avoid this error, I have also tried to not pass the bounds to the functions (since the code also provides theirs auto computing), and in that case I have problems matching dimensions in line 173.

volume_vertices = volume_vertices.reshape(3, *grid_size, order="F")

I don't know exactly the correct size of the AhmedML Dataset but once the grid size is defined for each axis, volume vertices generated should be rearranged to match dimensions, that rearrangement should also take into account the spacing value defined in the config.yaml file but it is never used in the provided code.

I hope to receive from you also some clarifications about how this part of the code works.

Minimum reproducible example

Relevant log output

Environment details

[Alexey-Kamenev]

Thank you for reporting this issue. Indeed, SetSamplingBounds takes min-max values for coordinates. We'll take a closer look and update this thread.

[LorenzoMattia]

Hi @Alexey-Kamenev, thank you for your answer. I hope to receive feedback early! By the way, I have temporarily skipped the problem leaving passing bounds = None to the process_directory function and leaving in this way the bounds computation to the code itself exploting the vtk GetBounds() function. Indeed, this function returns a tuple of six elements, which is perfectly suitable as inpu for SetSamplingBounds .

However, I still got troubles because of this line in the process_file function:

volume_vertices = volume_vertices.reshape(3, *grid_size, order="F")

The problem here was that volume_vertices is computed with get_cell_centers that exploits the vtk.vtkCellCenters() function which, apparently, skips cells that are blank or empty due to lack of data for some interpolation/computation.
Thus, I attach a proposal for a new version of get_cell_centers, still based on the vtk library, that should be more robust to this problem and guarantee dimensions matching between the requested voxel_grid and the computed cell_centers.

def get_cell_centers(voxel_grid):
    """Extracts the cell centers from a voxel grid.
    This version generates cell centers for ALL cells, including blank ones,
    by creating a regular grid based on the voxel grid's origin, spacing, and dimensions.
    """
    # Get grid properties
    dims = voxel_grid.GetDimensions()  # Number of points (nodes)
    origin = voxel_grid.GetOrigin()
    spacing = voxel_grid.GetSpacing()
    # Number of cells in each direction
    nx, ny, nz = dims[0] - 1, dims[1] - 1, dims[2] - 1
    # Create a regular grid of cell centers
    # Cell centers are offset by half spacing from origin
    x = origin[0] + spacing[0] * (0.5 + np.arange(nx))
    y = origin[1] + spacing[1] * (0.5 + np.arange(ny))
    z = origin[2] + spacing[2] * (0.5 + np.arange(nz))
    # Create meshgrid and flatten to get all cell centers
    X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
    centers = np.stack([X.ravel(order='F'), Y.ravel(order='F'), Z.ravel(order='F')], axis=1)
    return centers

This new implementation of get_cell_centers() generates cell centers mathematically rather than through the VTK filter.
A regular grid of coordinates is created for each axis (x, y, z), where centers are positioned halfway between grid nodes. Using np.meshgrid with indexing='ij', the 3D coordinates of all centers are generated and then flattened into an array of shape (N, 3) where N is the total number of cells.
This approach guarantees that the number of returned centers always matches exactly the product of the grid dimensions (nx × ny × nz), eliminating reshape issues caused by missing cells. The function maintains the same interface and output format as the original version, making it a drop-in replacement.

In this way I was able to run all the preprocessor script.

Hope this is helpful for you!