add 3D example of cutting of vtk meshes

Author: sclaus2

python/demo/cut_mesh_2D/cut_mesh_2D.py python/demo/cut_vtk_mesh/cut_mesh_2D.py

@@ -26,7 +26,7 @@ def create_rectangle_mesh(x0,y0,x1,y1,Nx,Ny):
 def mark_cells(mesh):
   cells = mesh.cells
   points = mesh.points
-  ls_values = np.zeros(3)
+  ls_values = np.zeros(len(points))
   inside_cells = np.empty(0, dtype=int)
   intersected_cells = np.empty(0, dtype=int)
 
@@ -53,7 +53,7 @@ def mark_cells(mesh):
 def create_cut_mesh(mesh,inner_mesh, intersected_cells):
   cell_type = cutcells.CellType.triangle
   triangulate = True
-  gdim = 3
+  gdim = 3 #vtk assumes 3 dimensions and for 2d sets third dimension to 0
 
   for cell_id in intersected_cells:
     c = mesh.get_cell(cell_id)

python/demo/cut_vtk_mesh/cut_mesh_3D.py

@@ -0,0 +1,86 @@
+import cutcells
+import numpy as np
+import pyvista as pv
+
+def level_set(x):
+  r = 0.6
+  c = np.array([0,0,0])
+  value = 0
+  for i in range(0,3):
+    value = value  + (x[i]-c[i])**2
+  value = np.sqrt(value)-r
+  return value
+
+def create_box_mesh(x0,y0,z0,x1,y1,z1,Nx,Ny,Nz):
+  x = np.linspace(x0, x1, num=Nx)
+  y = np.linspace(y0,y1, num=Ny)
+  z = np.linspace(z0,z1, num=Nz)
+  #produce grid of points by tensor product
+  xx, yy, zz = np.meshgrid(x, y, z)
+  points = np.c_[xx.reshape(-1), yy.reshape(-1), zz.reshape(-1)]
+  poly_points = pv.PolyData(points)
+  poly_mesh = poly_points.delaunay_3d()
+  grid = pv.UnstructuredGrid(poly_mesh)
+
+  return grid
+
+def mark_cells(mesh):
+  cells = mesh.cells
+  points = mesh.points
+  ls_values = np.zeros(len(points))
+  inside_cells = np.empty(0, dtype=int)
+  intersected_cells = np.empty(0, dtype=int)
+
+  id = 0
+
+  for i in range(mesh.n_cells):
+    num_points = cells[id]
+    id = id + 1
+
+    for j in range(num_points):
+       point = points[cells[id]]
+       ls_values[j] = level_set(point)
+       id = id + 1
+
+    domain_str = cutcells.classify_cell_domain(ls_values)
+    if domain_str=="inside":
+      inside_cells=np.append(inside_cells,i)
+
+    if domain_str=="intersected":
+      intersected_cells=np.append(intersected_cells,i)
+
+  return inside_cells, intersected_cells
+
+def create_cut_mesh(mesh,inner_mesh, intersected_cells):
+  cell_type = cutcells.CellType.tetrahedron
+  triangulate = True
+  gdim = 3
+
+  for cell_id in intersected_cells:
+    c = mesh.get_cell(cell_id)
+    vertex_coordinates = c.points.flatten()
+    points = c.points
+    j = 0
+    ls_values = np.zeros(len(points))
+    for point in points:
+      ls_values[j] = level_set(point)
+      j = j+1
+    cut_cell_int = cutcells.cut(cell_type, vertex_coordinates,  gdim, ls_values, "phi<0", triangulate)
+    grid_cell = pv.UnstructuredGrid(cut_cell_int.connectivity, cut_cell_int.types, cut_cell_int.vertex_coords)
+    inner_mesh = inner_mesh.merge(grid_cell)
+
+  return inner_mesh
+
+
+N = 10
+grid = create_box_mesh(-1,-1,-1, 1,1,1, N,N,N)
+
+inside_cells, intersected_cells = mark_cells(grid)
+
+extract = grid.extract_cells(inside_cells)
+
+mesh = create_cut_mesh(grid,extract,intersected_cells)
+
+mesh.plot(cpos="xy", show_edges=True)
+
+