如何使用mgrid在矩形和圓形之間進行插值

我正在嘗試創建一個 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()

uj5u.com熱心網友回復：

您嘗試使用的線性插值是一種應該作業的通用工具（有一個小警告）。所以問題首先在于你的矩形邊緣。

這是繪制內部和外部線條的完整性檢查：

# 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其余點（反之亦然）進行更改是一件簡單的事情。

Here's what I mean: you know that the rectangle's corner is at angle theta in your code (although I think you have x and y mixed up if we assume the conventional relationship between "x", "y" and the tangent of the angle). Since theta goes from 0 to pi/2, and your phi values also go from 0 to pi/2, you should choose index (array_length * (2*theta/np.pi)).round().astype(int) - 1 (or something similar) that will map to the rectangle's corner. If you have a square, this gives you theta = pi/4, and consequently (array_length / 2).round().astype(int) - 1. For array_length = 3 this is index (2 - 1) == 1, which is the middle index for 3-length arrays. (The more points you have along the edge, the less it will matter if you commit an off-by-one error here.)

The only remaining complication then is that we have to explicitly broadcast the 1d z array to the common shape. And we can use the same math you used to get a rectangular edge that is equidistant in angles.

Your code fixed with this suggestion (note that I've added 1 to the corner index because I'm using it as a right-exclusive range index):

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

The curves look reasonable:

As does the interpolated surface:

標籤：

上一篇：如何從WordPress主題檔案或插件重命名管理選單？

下一篇：從只有數字運算、沒有串列或字串的數字中洗掉偶數位