1

我正在尝试创建一个 3D 表面,其外部为 1/4 矩形,内部为 1/4 圆形。之前我得到了帮助,以椭圆作为外部创建 3D 表面,但由于某种原因,我无法为矩形执行此操作。我已经手工完成了有意义的数学计算,但我的代码没有。我将不胜感激任何帮助。

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

phase_plug = 0
phase_plug_dia = 20
plug_offset = 5
dome_dia = 28
# theta is angle where x and y intersect
theta = np.arctan(ellipse_x / ellipse_y)
# chi is for x direction and lhi is for y direction
chi = np.linspace(0, theta, 100)
lhi = np.linspace(theta, np.pi/2, 100)

# mgrid to create structured grid
r, phi = np.mgrid[0:1:array_length*1j, 0:np.pi/2:array_length*1j]

# Rectangle exterior, circle interior 
x = (ellipse_y * np.tan(chi)) * r + ((waveguide_throat / 2 * (1 - r)) * np.cos(phi))
y = (ellipse_x / np.tan(lhi)) * r + ((waveguide_throat / 2 * (1 - r)) * np.sin(phi))


# compute z profile
angle_factor = angle_factor / 10000
z = (ellipse_x / 2 * r / angle_factor) ** (1 / depth_factor)

plotter = pv.Plotter()

waveguide_mesh = pv.StructuredGrid(x, y, z)
plotter.add_mesh(waveguide_mesh)
plotter.show()

电流输出

4

1 回答 1

1

您尝试使用的线性插值是一种应该工作的通用工具(有一个小警告)。所以问题首先在于你的矩形边缘。

这是绘制内部和外部线条的完整性检查:

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

曲线看起来很合理: 固定边缘的屏幕截图

与插值曲面一样: 带有插值曲面线框的屏幕截图; 看起来不错

于 2022-01-26T01:07:28.390 回答