Python bmesh: Using selected edges on a curved object to create a face

Question

I am creating a parameterised model of a segmental tunnel-lining ring using Python and I have managed to create all the appropriate vertices and edges for the 'num_segs' number of segments in the ring. However, I have been having trouble creating the faces of the segments. I have written code that selects each combination of edges surrounding the 6 faces I want to create and then I run bpy.ops.mesh.edge_face_add(). But the results are not what I am looking for. The face that is created does not follow the curves and instead drops down on one side. I have tried with all edges selected in a clockwise and anti-clockwise manner and the issue persists. Here are my functions:

def create_faces(curve_res):
    # Object has been selected already
    obj = bpy.context.active_object
    bpy.context.view_layer.objects.active = obj
    bpy.ops.object.mode_set(mode='EDIT')
# Get the mesh data
mesh = bpy.context.object.data
bm = bmesh.from_edit_mesh(mesh)

all_edge_indices = generate_all_edge_indices(curve_res)

for i in range(6):
    edge_indices_to_select = all_edge_indices[i]

    for edge_index in edge_indices_to_select:
        bm.edges.ensure_lookup_table()
        bm.edges[edge_index].select_set(True)

    bmesh.update_edit_mesh(mesh)

    bpy.ops.mesh.edge_face_add()

    bmesh.update_edit_mesh(mesh)

    # Deselect edges
    for edge_index in edge_indices_to_select:
        bm.edges.ensure_lookup_table()
        bm.edges[edge_index].select_set(False)

    bpy.ops.object.mode_set(mode='OBJECT')


def generate_all_edge_indices(curve_res):
all_edge_indices = []

# To simplify
indices_01 = [i for i in range(curve_res)]
indices_23 = [curve_res + i for i in range(curve_res)]
indices_45 = [2*curve_res + i for i in range(curve_res)]
indices_67 = [3*curve_res + i for i in range(curve_res)]
index_02 = [4*curve_res]
index_13 = [4*curve_res + 1]
index_46 = [4*curve_res + 2]
index_57 = [4*curve_res + 3]
index_04 = [4*curve_res + 4]
index_15 = [4*curve_res + 5]
index_26 = [4*curve_res + 6]
index_37 = [4*curve_res + 7]

# 6 faces on all the shapes
for i in range(6):
    # Outer curve
    if i == 0:
        edge_indices_to_select = indices_01[::-1] + index_02 + indices_23 + index_13

    # Inner curve
    elif i == 1:
        edge_indices_to_select = indices_45 + index_57 + indices_67[::-1] + index_46

    # Front
    elif i == 2:
        edge_indices_to_select = indices_01 + index_15 + indices_45[::-1] + index_04

    # Back
    elif i == 3:
        edge_indices_to_select = indices_23[::-1] + index_26 + indices_67 + index_37

    # Left  
    elif i == 4:
        edge_indices_to_select = index_02 + index_04 + index_46 + index_26

    # Right
    else:
        edge_indices_to_select = index_13 + index_37 + index_57 + index_15

    all_edge_indices.append(edge_indices_to_select)

return all_edge_indices


def angle_of_vector(vector):
# Convert input vectors to NumPy arrays
vector = np.array(vector)

theta_radians = 0
vector[1] = 0

# Calculate the dot product
dot_product = np.dot((0, 0, 1), vector)

# Calculate the magnitudes
magnitude1 = np.linalg.norm((0, 0, 1))
magnitude2 = np.linalg.norm(vector)

# Calculate the cosine of the angle
cosine_theta = dot_product / (magnitude1 * magnitude2)

# Calculate the angle in radians using arccosine
theta_radians = np.arccos(cosine_theta) 

if vector[0] < 0:
    theta_radians = 2 * np.pi - theta_radians

return theta_radians


def add_circular_edge(bm, vertex_indices, radius, center_location):
    # Get the coordinates of the specified vertices
    vertex_vectors = [bm.verts[index].co for index in vertex_indices]
    vertex_coords = [(vertex_vector.x, vertex_vector.y, vertex_vector.z) for vertex_vector in vertex_vectors]
# Calculate start angle
start_angle = angle_of_vector(vertex_coords[0])
end_angle = angle_of_vector(vertex_coords[1])

if end_angle < start_angle:
    end_angle += 2 * np.pi

angle_between = end_angle - start_angle

num_steps = 16  # Adjust for smoothness

angles = np.linspace(start_angle, start_angle + angle_between, num_steps + 1)
# Add the curved edge vertices to the mesh
for ang in angles:
    x = center_location[0] + radius * math.sin(ang)
    y = center_location[1]
    z = center_location[2] + radius * math.cos(ang)

    # Create a new vertex at the calculated position
    new_vert = bm.verts.new((x, y, z))

    bm.verts.ensure_lookup_table()

    # Connect the new vertex to the previous one (except for the first iteration)
    if ang > start_angle:
        bm.edges.new([bm.verts[-2], new_vert])


def add_straight_edge(bm, vertex_indices):
    bm.edges.new((bm.verts[vertex_indices[0]], bm.verts[vertex_indices[1]]))

For the images below I put a breakpoint above bpy.ops.mesh.edge_face_add() and then ran this function manually.