How to use mgrid to interpolate between a rectangle and a circle

Question

I am trying to create a 3D surface that has a 1/4 rectangle for the exterior and 1/4 circle for the interior. I had help before to create the 3D surface with an ellipse as an exterior but I cannot do this for a rectangle for some reason. I have done the math by hand which makes sense, but my code does not. I would greatly appreciate any help with this.

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

Answer 1

The linear interpolation you're trying to use is a general tool that should work (with one small caveat). So the issue is first with your rectangular edge.

Here's a sanity check which plots your interior and exterior lines:

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

The bottom left is your interior circle, looks good. The top right is what's supposed to be a rectangle, but isn't.

To see why your original surface looks the way it does, we have to note one more thing (this is the small caveat I mentioned): the orientation of your curves is also the opposite. This implies that you interpolate the "top" (in the screenshot) point of your interior curve with the "bottom" point of the exterior curve. This explains the weird fan shape.

So you need to fix the exterior curve, and make sure the orientation of the two edges is the same. Note that you can just create the two 1d arrays for the two edges, and then interpolate them. You don't have to come up with a symbolic formula that you plug into the interpolation step. If you have 1d arrays of the same shape x_interior, y_interior, x_exterior, y_exterior then you can then do x_exterior * r + x_interior * (1 - r) and the same for y. This means removing the mgrid call, only using an array r of shape (n, 1), and making use of array broadcasting to do the interpolation. This means doing r = np.linspace(0, 1, array_length)[:, None].

So the question is how to define your rectangle. You need to have the same number of points on the rectangular curve than what you have on the circle (I would strongly recommend using the array_length parameter everywhere to ensure this!). Since you want to span the whole rectangle, I believe you have to choose an array index (i.e. a certain angle in the circular arc) which will map to the corner of the rectangle. Then it's a simple matter of varying only y for the points until that index, and x for the rest (or vice versa).

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: