您尝试使用的线性插值是一种应该工作的通用工具（有一个小警告）。所以问题首先在于你的矩形边缘。
这是绘制内部和外部线条的完整性检查：
# 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调用，仅使用rshape数组(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()
曲线看起来很合理：
与插值曲面一样：
