您尝试使用的线性插值是一种应该工作的通用工具(有一个小警告)。所以问题首先在于你的矩形边缘。
这是绘制内部和外部线条的完整性检查:
# debugging: plot interior and exterior
exterior_points = np.array([
ellipse_y * np.tan(chi),
ellipse_x / np.tan(lhi),
np.zeros_like(chi)
]).T
phi_aux = np.linspace(0, np.pi/2, array_length)
interior_points = np.array([
waveguide_throat / 2 * np.cos(phi_aux),
waveguide_throat / 2 * np.sin(phi_aux),
np.zeros_like(phi_aux)
]).T
plotter = pv.Plotter()
plotter.add_mesh(pv.wrap(exterior_points))
plotter.add_mesh(pv.wrap(interior_points))
plotter.show()
左下角是你的内圈,看起来不错。右上角应该是矩形,但不是。
要了解为什么您的原始曲面看起来如此,我们还必须注意一件事(这是我提到的小警告):曲线的方向也是相反的。这意味着您将内部曲线的“顶部”(在屏幕截图中)点与外部曲线的“底部”点进行插值。这就解释了奇怪的扇形。
所以你需要修复外部曲线,并确保两个边缘的方向相同。请注意,您可以为两条边创建两个一维数组,然后对它们进行插值。您不必想出插入插值步骤的符号公式。如果您有相同形状的一维数组,x_interior, y_interior, x_exterior, y_exterior
那么您可以x_exterior * r + x_interior * (1 - r)
对y
. 这意味着删除mgrid
调用,仅使用r
shape数组(n, 1)
,并利用数组广播进行插值。这意味着做r = np.linspace(0, 1, array_length)[:, None]
.
所以问题是如何定义你的矩形。矩形曲线上的点数需要与圆上的点数相同(我强烈建议在array_length
任何地方使用该参数来确保这一点!)。由于您想跨越整个矩形,我相信您必须选择一个数组索引(即圆弧中的某个角度),它将映射到矩形的角。y
然后,仅针对该索引之前的点以及x
其余点(反之亦然)进行更改是一件简单的事情。
这就是我的意思:你知道矩形的角在你的代码中是成角度的(尽管如果我们假设“x”、“y”和角度的切线之间的传统关系,theta
我认为你有x
并且混淆了)。y
由于theta
从 0 到 pi/2,并且您的phi
值也从 0 到 pi/2,您应该选择(array_length * (2*theta/np.pi)).round().astype(int) - 1
将映射到矩形角的索引(或类似的东西)。如果你有一个正方形,这会给你theta = pi/4
,因此(array_length / 2).round().astype(int) - 1
。为此array_length = 3
是 index (2 - 1) == 1
,它是 3 长度数组的中间索引。(你在边缘的点越多,如果你在这里犯一个错误的问题就越不重要。)
剩下的唯一复杂性是我们必须将一维z
数组显式广播到公共形状。我们可以使用相同的数学方法来获得角度等距的矩形边。
您的代码已使用此建议修复(请注意,我已将 1 添加到角索引,因为我将其用作右排他范围索引):
import numpy as np
import pyvista as pv
# parameters for the waveguide
# diameter of the inner circle
waveguide_throat = 30
# axes of the outer ellipse
ellipse_x = 250
ellipse_y = 170
# shape parameters for the z profile
depth_factor = 4
angle_factor = 40
# number of grid points in radial and angular direction
array_length = 100
# quarter circle interior line
phi = np.linspace(0, np.pi/2, array_length)
x_interior = waveguide_throat / 2 * np.cos(phi)
y_interior = waveguide_throat / 2 * np.sin(phi)
# theta is angle where x and y intersect
theta = np.arctan2(ellipse_y, ellipse_x)
# find array index which maps to the corner of the rectangle
corner_index = (array_length * (2*theta/np.pi)).round().astype(int)
# construct rectangular coordinates manually
x_exterior = np.zeros_like(x_interior)
y_exterior = x_exterior.copy()
phi_aux = np.linspace(0, theta, corner_index)
x_exterior[:corner_index] = ellipse_x
y_exterior[:corner_index] = ellipse_x * np.tan(phi_aux)
phi_aux = np.linspace(np.pi/2, theta, array_length - corner_index, endpoint=False)[::-1] # mind the reverse!
x_exterior[corner_index:] = ellipse_y / np.tan(phi_aux)
y_exterior[corner_index:] = ellipse_y
# interpolate between two curves
r = np.linspace(0, 1, array_length)[:, None] # shape (array_length, 1) for broadcasting
x = x_exterior * r + x_interior * (1 - r)
y = y_exterior * r + y_interior * (1 - r)
# debugging: plot interior and exterior
exterior_points = np.array([
x_exterior,
y_exterior,
np.zeros_like(x_exterior),
]).T
interior_points = np.array([
x_interior,
y_interior,
np.zeros_like(x_interior),
]).T
plotter = pv.Plotter()
plotter.add_mesh(pv.wrap(exterior_points))
plotter.add_mesh(pv.wrap(interior_points))
plotter.show()
# compute z profile
angle_factor = angle_factor / 10000
z = (ellipse_x / 2 * r / angle_factor) ** (1 / depth_factor)
# explicitly broadcast to the shape of x and y
z = np.broadcast_to(z, x.shape)
plotter = pv.Plotter()
waveguide_mesh = pv.StructuredGrid(x, y, z)
plotter.add_mesh(waveguide_mesh, style='wireframe')
plotter.show()
曲线看起来很合理:
与插值曲面一样: