pyvista / pyvista-support

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

Boolean operations with UnstructuredGrid #176

Closed pempmu closed 4 years ago

pempmu commented 4 years ago

Hi pyvista community,

I have recently come across this great project and I am currently using its capabilities to do some manufacturing related simulations with it. For my purpuse I need to make usage of Boolean Operations, but I have obtained a first error that I am not quite sure how to figure out. Essentially, what I am trying to achieve is perform a Boolean operation between a unstructured grid and structured grid. Here you find the full code to reproduce the issue:

import numpy as np
import math as m
import panel as pn
import pyvista as pv

pn.extension('vtk')

### The first part of the script generates the unstructured grid shown in the attached picture

A=np.array([1.0,0.0,0.0])

def normalize_vector(v):
    norm = np.linalg.norm(v)
    if norm == 0:
       return v
    return v / norm

def generate_perpendicular(v):
    if v[0]==0:
        return [1.0,0.0,0.0]
    if v[1]==0:
        return [0.0,1.0,0.0]   
    if v[2]==0:
        return [0.0,0.0,1.0]
    if v[0]!=0 and v[1]!=0 and v[2]!=0 :
        return normalize_vector([1.0,1.0,(-v[0]-v[1])/v[2]])

def point_gen(t,r,n,u,c):
    return r*m.cos(t*m.pi/180)*u+r*m.sin(t*m.pi/180)*np.cross(n,u)+c

def tool_fine(tool_raw):
    results=[]
    for i in range(len(tool_raw)+1):
        if i==0:
            rdif=tool_raw[0][0]
            res = [[(rdif/10*x+rdif/10),0.0] for x in range(10-1)]
        if i==len(tool_raw):
            rdif=tool_raw[-1][0]
            res = [[(-0.1*x+rdif),tool_raw[-1][1]] for x in range(10)]          
        if i!=len(tool_raw) and i!=0:
            rdif=tool_raw[i][0]-tool_raw[i-1][0]
            zdif=tool_raw[i][1]-tool_raw[i-1][1]
            res = [[tool_raw[i-1][0]+rdif/10*x,tool_raw[i-1][1]+zdif/10*x] for x in range(10)]     
        results+=res
    return results

def tool_3d(tool,n):
    result=[]
    u=np.array(generate_perpendicular(n))
    for i in range(len(tool)):
        r=tool[i][0]
        z=tool[i][1]
        c=np.array([n[0]*z,n[1]*z,n[2]*z])
        res=[point_gen(x*360/60,r,n,u,c) for x in range(60)]
        result+=res
    return result

def tool2d(tool):
    t_i=tool[::-1]
    side2=[[-t_i[i][0],t_i[i][1]] for i in range(len(t_i))]
    return tool+side2

n=normalize_vector(A)
u=np.array(generate_perpendicular(A))

c=np.array([0.0,0.0,0.0])
r=1.0
t=7.0

tool_raw=[[1.0,0.0],[1.2,0.5],[1.2,5.0]]
tool=tool_fine(tool_raw)
result = tool_3d(tool,n)

x=[result[i][0] for i in range(len(result))]
y=[result[i][1] for i in range(len(result))]
z=[result[i][2] for i in range(len(result))]

point_cloud = pv.PolyData(np.array(result))

volume = point_cloud.delaunay_3d()

### The second part of the script generates the structured grid 

def make_cube():
    x = np.linspace(0.4, -1.5, 25)
    y = np.linspace(-0.4, -1.5, 25)
    z = np.linspace(-0.4, -1.5, 25)
    grid = pv.StructuredGrid(*np.meshgrid(x, y, z))
    return grid.extract_surface().triangulate()

cube = make_cube()

### The third part of the script generates the plot shown in the picture in a jupyter notebook 

pl = pv.Plotter(notebook=True);
pl.add_mesh(volume,show_edges=True)
pl.add_mesh(cube, color="red", show_edges=False, opacity=0.5)

orientation_widget = True
vtkpan = pn.panel(pl.ren_win, sizing_mode='stretch_both', height=600)
vtkpan

The result of executing this code in a jupyter notebook is shown in the attached picture. The problem apperars when executing the next line of code to perform the boolean:

cut = cube - volume

The following error appears:

AttributeError: 'UnstructuredGrid' object has no attribute 'is_all_triangles'

Could anyone give me a hint how to interpret this error and find a solution?

Thanks in advance!

pic
akaszynski commented 4 years ago

AFAIK, vtk doesn't support boolean operations between unstructured grids, only PolyData. We need to add type checking to provide a better error message.

You might be able to use select_enclosed_points to subselect a part of the unstructred grid and then extract those points. voxelize uses this:

def voxelize(mesh, density):
    x_min, x_max, y_min, y_max, z_min, z_max = mesh.bounds
    x = np.arange(x_min, x_max, density)
    y = np.arange(y_min, y_max, density)
    z = np.arange(z_min, z_max, density)
    x, y, z = np.meshgrid(x, y, z)

    # Create unstructured grid from the structured grid
    grid = pyvista.StructuredGrid(x, y, z)
    ugrid = pyvista.UnstructuredGrid(grid)

    # get part of the mesh within the mesh
    selection = ugrid.select_enclosed_points(mesh, tolerance=0.0)
    mask = selection.point_arrays['SelectedPoints'].view(np.bool)

    # extract cells from point indices
    return ugrid.extract_points(mask)

Your boolean operations won't be smooth (i.e. the surface of the cut will be quite rough), but it will at least get you started.

pempmu commented 4 years ago

@akaszynski thank you for the hints. I will try to figure out some solution for my problem and get back once I am satisfied with the results.

banesullivan commented 4 years ago

My impression: there's no need for you to use UnstructuredGrids. Try switching to PolyData and using PyVista's Box helper to create your cutter:

def make_cube():
    bounds = [-1.5, 0.4, -1.5, -0.4, -1.5, -0.4]
    return pv.Box(bounds).triangulate()

cube = make_cube()

# You're point cloud is only a surface, so extract the bounding surface of the delaunay filter and clean it
volume = point_cloud.delaunay_3d().extract_surface().clean()

# Now you can do the cut
cut = cube - volume

pl = pv.Plotter(notebook=False, shape=(1,2));
pl.add_mesh(volume, show_edges=True, opacity=0.5, color='grey')
pl.add_mesh(cube, color="orange", show_edges=False, opacity=0.5)
pl.add_mesh(cut, color="red")
pl.subplot(0,1)
pl.add_mesh(cut, color="red")
pl.link_views()
pl.show()
Screen Shot 2020-06-20 at 9 49 00 PM
pempmu commented 4 years ago

Great, this will keep me working on my project. Thank you for the tip @banesullivan