如何在pyvista中以填充单元格/网格的形式填充两个网格之间的距离？

Question

我看一下这个sample，我想创建一个网格/曲面网格来填充这个距离。如何在PyVista中执行这样的操作

我所尝试的，似乎无法与Andras Deak beautiful answer搭桥：

import pyvista as pv
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import KDTree
import PVGeo
from PVGeo import interface
from PVGeo.filters import BuildSurfaceFromPoints

b = pv.read('./top.vtk') # PolyData
t = pv.read('./bottom.vtk') # PolyData

dim = (int(b.bounds[1]-b.bounds[0]), int(b.bounds[3]-b.bounds[2]), 1)
z_range =  np.arange(b.bounds[4], b.bounds[5] )
bottom = BuildSurfaceFromPoints().apply(b)
top = BuildSurfaceFromPoints().apply(t)

grid_2d = top.points.reshape(dim[:-1] + (3,), order='F')[..., :-1]

不幸的是，这在网格2d上失败了

ValueError: cannot reshape array of size 1942464 into shape (30150,26750,3)

Answer 1

我不知道是否有一种内置的方法可以在两个曲面之间进行插值，但仅使用numpy就不难做到这一点

下面是一个使用Perlin noise在同一网格上以两种不同高度生成两张数据表的示例。答案的实际代码在后面

import numpy as np
import pyvista as pv

# generate two sheets of input data
noise = pv.perlin_noise(2, (0.2, 0.2, 0.2), (0, 0, 0))
bounds_2d = (-10, 10, -10, 10)
dim = (40, 50, 1)
bottom, top = [
    pv.sample_function(noise, dim=dim, bounds=bounds_2d + (z, z)).warp_by_scalar()
    for z in [-5, 5]
]

# actual answer starts here
# the top and bottom sheets are named `top` and `bottom`
# and they share the same 2d grid

# rebuild grid points
grid_2d = top.points.reshape(dim[:-1] + (3,), order='F')[..., :-1]
values_x = grid_2d[:, 0, 0]
values_y = grid_2d[0, :, 1]

# generate full grid with equidistant interpolation in each (x, y)
nz = 10
scale = np.linspace(0, 1, nz)
scale_z = scale[:, None] * [0, 0, 1]  # shape (nz, 3)
scale_z_inv = (1 - scale[:, None]) * [0, 0, 1]  # shape (nz, 3)
z_bottom = bottom.points.reshape(dim[:-1] + (3,), order='F')[..., -1]  # shape (nx, ny)
z_top = top.points.reshape(dim[:-1] + (3,), order='F')[..., -1]  # shape (nx, ny)
interpolated_z = scale * z_bottom[..., None] + (1 - scale) * z_top[..., None]  # shape (nx, ny, nz)

grid_2d_in_3d = np.pad(grid_2d, [(0, 0), (0, 0), (0, 1)])  # shape (nx, ny, 3)
final_grid = grid_2d_in_3d[..., None, :] + interpolated_z[..., None] * [0, 0, 1]  # shape (nx, ny, nz, 3)
mesh = pv.StructuredGrid(*final_grid.transpose())

# plot the two sheets and the interpolated grid
pv.set_plot_theme('document')
plotter = pv.Plotter()
plotter.add_mesh(bottom, show_scalar_bar=False)
plotter.add_mesh(top, show_scalar_bar=False)
plotter.add_mesh(mesh, style='wireframe')
plotter.show()