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

问题描述

所以我看这个sample，我想创建一个网格/曲面网格来填充这个距离。如何在PyVista中做这样的事情？

我尝试但似乎无法与安德拉斯·迪克美人沟通的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]

遗憾的是，在GRID_2d与

ValueError：无法将大小为1942464的数组整形 (30150,26750，3)

---

解决方案

我不知道两个曲面之间是否有内置的插补方法，但仅使用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()