python - 如何在pyvista中填充两个网格作为填充单元格/网格之间的距离？

Question

所以我查看了这个示例，我想创建一个填充该距离的网格/表面网格。如何在 PyVista 中做这样的事情？

我尝试但似乎无法与 Andras Deak 美丽的答案桥接：

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]

遗憾的是，这与 grid_2d 失败了

ValueError：无法将大小为 1942464 的数组重塑为形状 (30150,26750,3)

Answer 1

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

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

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()