3D Surface in PyVista is not generating correctly. I am trying to avoid closing the center opening but can't figure out how to do so

Question

Issue with PyVista filling opening that shouldn't be filled

import math
import numpy as np
import matplotlib as mpl
import pyvista as pv


mpl.use("Qt5Agg")
mpl.rcParams["toolbar"] = "None"  # Get rid of toolbar


def xy_waveguide_contour(throat, x_waveguide, ellipse_x):
    x_initial = (throat + (x_waveguide * (ellipse_x - throat))) / 2

    return x_initial


def xy_waveguide_contour_x2(throat, x_waveguide, ellipse_x):
    x_initial = (throat + (x_waveguide * (ellipse_x - throat)))

    return x_initial


def z_waveguide_contour(x_array, depth_factor, angle_factor, throat):
    angle_factor = angle_factor / 10000

    x_prime = x_array - (throat / 2)

    z = (x_prime / angle_factor) ** (1 / depth_factor)

    return z


def ellipse_contour(a, b):
    a = a / 2
    b = b / 2
    ellipse_steps = np.linspace(0, 0.5 * math.pi, 100)
    x_ellipse_array = np.array([])
    y_ellipse_array = np.array([])
    for h in range(100):
        x_ellipse_array = np.append(x_ellipse_array, a * np.cos(ellipse_steps[h]))
        y_ellipse_array = np.append(y_ellipse_array, b * np.sin(ellipse_steps[h]))
    return x_ellipse_array, y_ellipse_array


def circle_contour(throat):
    throat = (throat / 2)
    circle_steps = np.linspace(0, 0.5 * math.pi, 100)
    x_circle_array = np.array([])
    y_circle_array = np.array([])
    for j in range(100):
        x_circle_array = np.append(x_circle_array, throat * np.cos(circle_steps[j]))
        y_circle_array = np.append(y_circle_array, throat * np.sin(circle_steps[j]))
    return x_circle_array, y_circle_array


waveguide_throat = 30

ellipse_x = 250
ellipse_y = 150
depth_fact = 4
angle_fact = 40

# Total steps = 100
array_length = 100

xy_steps = np.linspace(0, 1, array_length)

# initialize x, y, z, and zero array
x_array = np.array([])
y_array = np.array([])
z_array = np.array([])
xsub_array = np.array([])
ysub_array = np.array([])
zero_array = np.zeros([array_length])

# calculate hor(x), ver(y), and height(z) contour data and add into array
for i in range(array_length):
    x_array = np.append(x_array, xy_waveguide_contour(waveguide_throat, xy_steps[i], ellipse_x))
    y_array = np.append(y_array, xy_waveguide_contour(waveguide_throat, xy_steps[i], ellipse_y))

    z_array = np.append(z_array, z_waveguide_contour(x_array[i], depth_fact, angle_fact, waveguide_throat))

# calculate data for ellipse
x_ellipse_data, y_ellipse_data = ellipse_contour(ellipse_x, ellipse_y)

# grab last point from z_array and make entire array same value to define height of ellipse/waveguide
ellipse_height = z_array[array_length - 1]

ellipse_z = np.full(shape=array_length, fill_value=ellipse_height)

# Calculate data for throat
circle_x, circle_y = circle_contour(waveguide_throat)


for j in range(0, array_length, 20):
    for i in range(array_length):

        xsub_array = np.append(xsub_array, xy_waveguide_contour_x2(circle_x[j], xy_steps[i], x_ellipse_data[j]))
        ysub_array = np.append(ysub_array, xy_waveguide_contour_x2(circle_y[j], xy_steps[i], y_ellipse_data[j]))

# X = np.concatenate((circle_x, x_ellipse_data, x_array, zero_array))
# Y = np.concatenate((circle_y, y_ellipse_data, zero_array, y_array))
# Z = np.concatenate((zero_array, ellipse_z, z_array, z_array))

# Reshape arrays into 1 column, multiple rows
x_array = x_array.reshape(-1, 1)
y_array = y_array.reshape(-1, 1)
z_array = z_array.reshape(-1, 1)
ellipse_z = ellipse_z.reshape(-1, 1)
x_ellipse_data = x_ellipse_data.reshape(-1, 1)
y_ellipse_data = y_ellipse_data.reshape(-1, 1)
circle_x = circle_x.reshape(-1, 1)
circle_y = circle_y.reshape(-1, 1)
zero_array = zero_array.reshape(-1, 1)
xsub_array = xsub_array.reshape(-1, 1)
ysub_array = ysub_array.reshape(-1, 1)

# save arrays to text

X = np.concatenate((circle_x, x_ellipse_data, x_array, zero_array, xsub_array), axis=0)
Y = np.concatenate((circle_y, y_ellipse_data, zero_array, y_array, ysub_array), axis=0)
Z = np.concatenate((zero_array, ellipse_z, z_array, z_array, z_array, z_array, z_array, z_array, z_array), axis=0)


xyz = np.concatenate((X, Y, Z), axis=1)

cloud = pv.PolyData(xyz)
surf = cloud.delaunay_2d()
surf.plot()

I am trying to create a 3D surface using pyvista (Have also tried Mayavi), whenever I perform a delaunay_2D mesh to create the surface, it closes the "mouth" opening of this surface that should still be open. I have included an image showing what part needs to be fixed and also a copy of my code that generates the data and reproduces the current issue. I would really appreciate anyone's help regarding this issue.

Answer 1

As I noted in a comment under a previous version of your question, I can't find an easy way to fix your current approach.

1. You could pass alpha=cloud.length/10 to delaunay_2d() which would be close enough, but this would still leave spurious fringes at the edge of your waveguide.
2. It would be nice to be able to create your waveguide using extrusion, but I don't see how that would be applicable here.

So the only thing I can think of is to change your whole approach: create your waveguide as a 2d structured grid, by parametrising your surface. This is actually possible, and not too difficult:

1. In 2d you have a circle and an ellipse, and linear interpolation between the two. This means taking points in the form of (r*cos(phi), r*sin(phi)) and interpolating to points at (a*cos(phi), b*sin(phi)). This is straightforward.
2. In the z direction you have a root function, which depends on the "radial" coordinate of your 2d grid.

Here's how you can do 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 = 150
# shape parameters for the z profile
depth_factor = 4
angle_factor = 40
# number of grid points in radial and angular direction
array_length = 100

# now create the actual structured grid
# 2d circular grid
r, phi = np.mgrid[0:1:array_length*1j, 0:np.pi/2:array_length*1j]

# transform to ellipse on the outside, circle on the inside
x = (ellipse_x/2 * r + waveguide_throat/2 * (1 - r))*np.cos(phi)
y = (ellipse_y/2 * 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)

waveguide = pv.StructuredGrid(x, y, z)
waveguide.plot(show_edges=True)

This gives you a dense 2d grid forming your waveguide. Whatever you need to do with this surface, you can do with the structured grid.

Here's the structured grid, with your point cloud laid over it to show they are the same: