0

所以我查看了这个示例,我想创建一个填充该距离的网格/表面网格。如何在 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)

4

1 回答 1

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

绘制两个 Perlin 噪声表面和中间的密集网格

于 2021-07-01T09:49:03.303 回答