search

pyvista / pyvista-support

[moved] Please head over to the Discussions tab of the PyVista repository
https://github.com/pyvista/pyvista/discussions
59 stars 4 forks source link

crash when loading a vtk unstructured grid #277

Open Boorhin opened 3 years ago

Boorhin commented 3 years ago

Description

I have received a mesh generated with SMS in a netcdf format and I have written a little script to convert it in vtk legacy format. Hopefully this was done correctly. Although in the original mesh, it seems that nodes are cells and points are element (correct me if I am wrong). The resulting unstructured grid is strange as cells(nodes) have between 1 and 9 elements. What I did is filter any cell with less than 3 elements and generated a vtk with it. However, when I try to read the created vtk, bott = pv.read("ClydeMesh.vtk") the kernel dies and I get this output.

$ jupyter qtconsole Fatal Python error: Segmentation fault

Thread 0x00007fe617fff700 (most recent call first): File "/usr/lib/python3/dist-packages/zmq/utils/garbage.py", line 47 in run File "/usr/lib/python3.8/threading.py", line 932 in _bootstrap_inner File "/usr/lib/python3.8/threading.py", line 890 in _bootstrap

Thread 0x00007fe62cc33700 (most recent call first): File "/usr/lib/python3/dist-packages/ipykernel/parentpoller.py", line 39 in run File "/usr/lib/python3.8/threading.py", line 932 in _bootstrap_inner File "/usr/lib/python3.8/threading.py", line 890 in _bootstrap

Thread 0x00007fe62d474700 (most recent call first): File "/usr/lib/python3.8/threading.py", line 302 in wait File "/usr/lib/python3.8/threading.py", line 558 in wait File "/usr/lib/python3/dist-packages/IPython/core/history.py", line 829 in run File "/usr/lib/python3/dist-packages/IPython/core/history.py", line 58 in needs_sqlite File "", line 2 in run File "/usr/lib/python3.8/threading.py", line 932 in _bootstrap_inner File "/usr/lib/python3.8/threading.py", line 890 in _bootstrap

Thread 0x00007fe62ecb7700 (most recent call first): File "/usr/lib/python3/dist-packages/ipykernel/heartbeat.py", line 100 in run File "/usr/lib/python3.8/threading.py", line 932 in _bootstrap_inner File "/usr/lib/python3.8/threading.py", line 890 in _bootstrap

Thread 0x00007fe62f4b8700 (most recent call first): File "/usr/lib/python3.8/selectors.py", line 468 in select File "/usr/lib/python3.8/asyncio/base_events.py", line 1823 in _run_once File "/usr/lib/python3.8/asyncio/base_events.py", line 570 in run_forever File "/usr/lib/python3/dist-packages/tornado/platform/asyncio.py", line 132 in start File "/usr/lib/python3/dist-packages/ipykernel/iostream.py", line 78 in _thread_main File "/usr/lib/python3.8/threading.py", line 870 in run File "/usr/lib/python3.8/threading.py", line 932 in _bootstrap_inner File "/usr/lib/python3.8/threading.py", line 890 in _bootstrap

Current thread 0x00007fe632e66740 (most recent call first): File "/home/julien/.local/lib/python3.8/site-packages/pyvista/utilities/fileio.py", line 156 in read_legacy File "/home/julien/.local/lib/python3.8/site-packages/pyvista/utilities/fileio.py", line 239 in read File "", line 1 in File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3331 in run_code File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3254 in run_ast_nodes File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3062 in run_cell_async File "/usr/lib/python3/dist-packages/IPython/core/async_helpers.py", line 68 in _pseudo_sync_runner File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2886 in _run_cell File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2857 in run_cell File "/usr/lib/python3/dist-packages/ipykernel/zmqshell.py", line 536 in run_cell File "/usr/lib/python3/dist-packages/ipykernel/ipkernel.py", line 300 in do_execute File "/usr/lib/python3/dist-packages/tornado/gen.py", line 326 in wrapper File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 543 in execute_request File "/usr/lib/python3/dist-packages/tornado/gen.py", line 326 in wrapper File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 268 in dispatch_shell File "/usr/lib/python3/dist-packages/tornado/gen.py", line 326 in wrapper File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 365 in process_one File "/usr/lib/python3/dist-packages/tornado/gen.py", line 1162 in run File "/usr/lib/python3/dist-packages/tornado/gen.py", line 1248 in inner File "/usr/lib/python3/dist-packages/tornado/stack_context.py", line 300 in null_wrapper File "/usr/lib/python3/dist-packages/tornado/ioloop.py", line 758 in _run_callback File "/usr/lib/python3.8/asyncio/events.py", line 81 in _run File "/usr/lib/python3.8/asyncio/base_events.py", line 1859 in _run_once File "/usr/lib/python3.8/asyncio/base_events.py", line 570 in run_forever File "/usr/lib/python3/dist-packages/tornado/platform/asyncio.py", line 132 in start File "/usr/lib/python3/dist-packages/ipykernel/kernelapp.py", line 583 in start File "/usr/lib/python3/dist-packages/traitlets/config/application.py", line 664 in launch_instance File "/usr/lib/python3/dist-packages/ipykernel_launcher.py", line 16 in File "/usr/lib/python3.8/runpy.py", line 87 in _run_code File "/usr/lib/python3.8/runpy.py", line 194 in _run_module_as_main [JupyterQtConsoleApp] KernelRestarter: restarting kernel (1/5), keep random ports [JupyterQtConsoleApp] WARNING | kernel restarted

Example Data

VTK zipped: https://drive.google.com/file/d/1WLvEdzNCEkvx14eqpdLjt-NJwkoYlyb9/view?usp=sharing

akaszynski commented 3 years ago

This fails using pure vtk with just:

import vtk
reader = vtk.vtkDataSetReader()
reader.SetFileName('ClydeMesh.vtk')
reader.Update()
# segfault follows...

My guess is the exported file isn't correctly formatted. I suggest that instead of converting your netcdf to a legacy VTK file, you instead convert it to a pyvista.UnstructuredGrid. Might be easier to catch bugs. If you need to save it as a file, you can always use my_grid.save('my_cool_mesh.vtk') for later reference.

banesullivan commented 3 years ago

yeah, this isn't really something we can resolve on our end. The VTK file was improperly generated. If it helps, here is an error message I captured:

vtkUnstructuredGrid (0x7f82fbd5d5d0): Cell array bathymetry with 1 components, has 29729 tuples but there are only 0 cells
Boorhin commented 3 years ago

I will work on it, this is from a FVCOM model which is a standard in ocean modelling so if I manage I will post the results here

Boorhin commented 3 years ago

I simply made a polydata and not a mesh and it worked fine I have now generated a correct legacy ascii and I can load it as a mesh but have been unsuccessful making one from scratch in pyvista so far. If you want to play with the data I attach a zip with the arrays extracted from the netcdf4 file arrays.zip

import netCDF4 as nc
import numpy as np
import pyvista as pv
from pyvistaqt import BackgroundPlotter
def loadarr():
    Nx= np.load("Nx.npy")
    Ny=np.load("Ny.npy")
    bath=np.load("bath.npy")
    cells=np.load("cells.npy")
    return Nx, Ny, cells, bath

def Make_legacy(Nx,Ny,bath,C):
    headerp = "# vtk DataFile Version 2.0\nUnstructured Grid Example\nASCII\n\nDATASET UNSTRUCTURED_GRID\nPOINTS {} float\n".format(len(Nx))
    with open("ClydeMesh2.vtk","w") as f:
        f.write(headerp)
        for i in range(len(Nx)):
            f.write("{:.3f} {:.3f} {:.3f}\n".format(Nx[i],Ny[i],bath[i]*-1))
        f.write("\nCELLS {} {}\n".format(len(C), len(C)*4))
        for j in range(len(C)):
            f.write("3 {:d} {:d} {:d}\n".format(C[j,1],C[j,2],C[j,3]))
        CT = ""
        for k in range(len(C)):
            CT += "5 "
            if k%20 == 0:
                 CT+="\n"
        f.write("\nCELL_TYPES {:d}\n{}".format(len(C),CT))
        dP =""
        for h in range(len(Nx)):
            dP +="{:.1f} ".format(bath[h])
            if h%20 == 0:
                 dP+="\n"
        f.write("\n POINT_DATA {}\nSCALARS Bathymetry float 1\nLOOKUP_TABLE default\n{}".format(len(Nx), dP))

Nx, Ny, cells, bath = loadarr()
C = np.ones((len(cells[0]), 4), dtype= np.int)*3
C[:,1:] = cells[:].T-1
 #Make_legacy(Nx,Ny,bath,C)
points = np.zeros((len(Nx), 3))
points[:,0],points[:,1],points[:,2]= Nx.T, Ny.T, bath.T
#CellType = np.ones(len(C),dtype=np.int)*vtk.VTK_TRIANGLE
PD= pv.PolyData(points,C)
plotter = BackgroundPlotter(show = True)
plotter.add_mesh(PD, scalars=bath,flip_scalars=True)
akaszynski commented 3 years ago

This creates an unstructured grid by extending your script :

import vtk
cells = C.ravel()
celltypes = np.ones(len(C), dtype=np.int)*vtk.VTK_TRIANGLE
offset = np.arange(0, cells.size, 4)
grid = pv.UnstructuredGrid(offset, cells, celltypes, points)

You can also convert your polydata to an unstructured grid with cast_to_unstructured_grid. Since your data is only a triangular mesh, you may want to consider just keeping it a PolyData