Python script fails after 2 iterations

Question

I wrote a script that takes a specified octahedron and refines it. For context, I have pictures of the output with iterations set to 0,1, and 2.

When I set iterations=3, however, I get the following error:

Traceback (most recent call last):
  File "/Users/amaurydeburgos/Documents/WinterBreak2023.py", line 78, in <module>
    ax.add_collection3d(Poly3DCollection([ [ Vertices[i] for i in f] for f in Faces],edgecolors='k',facecolors='w'))
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/mpl_toolkits/mplot3d/art3d.py", line 701, in __init__
    super().__init__(verts, *args, **kwargs)
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/matplotlib/_api/deprecation.py", line 454, in wrapper
    return func(*args, **kwargs)
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/matplotlib/collections.py", line 1176, in __init__
    self.set_verts(verts, closed)
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/mpl_toolkits/mplot3d/art3d.py", line 745, in set_verts
    self.get_vector(verts)
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/mpl_toolkits/mplot3d/art3d.py", line 734, in get_vector
    xs, ys, zs = np.row_stack(segments3d).T
  File "<__array_function__ internals>", line 200, in vstack
  File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/numpy/core/shape_base.py", line 296, in vstack
    return _nx.concatenate(arrs, 0, dtype=dtype, casting=casting)
  File "<__array_function__ internals>", line 200, in concatenate
ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 3 and the array at index 689 has size 0

I don't know what's going wrong. Any help is greatly appreciated. The script I wrote is below:

import pylab as plt
import numpy as np
from numpy import array
from numpy.linalg import norm
from operator import add
from itertools import combinations
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

# Initial octahedron

Vertices = [ (0,0,1), (1,0,0), (0,1,0), (-1,0,0), (0,-1,0), (0,0,-1) ]
Edges = { frozenset({0,1}), frozenset({0,2}), frozenset({0,3}), frozenset({0,4}),
          frozenset({1,2}), frozenset({2,3}), frozenset({3,4}), frozenset({1,4}),
          frozenset({1,5}), frozenset({2,5}), frozenset({3,5}), frozenset({4,5}) }
Faces = { frozenset({0,1,2}), frozenset({0,2,3}), frozenset({0,3,4}), frozenset({0,1,4}),
          frozenset({1,2,5}), frozenset({2,3,5}), frozenset({3,4,5}), frozenset({1,4,5}) }

iterations = 3

for i in range(iterations):

    # Initializing set of new vertices, new edges, and new faces

    counter = len(Vertices)-1

    newVertices = []
    newEdges = set()
    newFaces = set()    

    # Adding elements of newVertices and trivial elements of newEdges

    for edge in Edges:
        counter = counter+1
        newVertex = np.array([0,0,0])
        for vertex in edge:
            newVertex = np.add(newVertex,np.array(Vertices[vertex]))
        
            newEdge = frozenset({vertex,counter})
            newEdges.add(newEdge)
        newVertex = np.divide(newVertex,norm(newVertex,2))
        newVertices.append(tuple(newVertex))

    # Adding non-trivial elements of newEdges and elements of newFaces

    for face in Faces:
        middleFace=set()
        SpecialEdges=set()
        pairsOfSpecialEdges={(a,b) for a,b in combinations({edge for edge in newEdges if len(edge.intersection(face))==1},2) if (a & b) and len((a-face).intersection(b-face))==1}
        for pair in pairsOfSpecialEdges:
            for edge in pair:
                SpecialEdges.add(edge)
        for vertex in face:
            incidentEdges={edge for edge in SpecialEdges if vertex in edge}
            newEdge=set()
            for edge in incidentEdges:
                for v in edge:
                    if v==vertex:
                        continue
                    else:
                        newEdge.add(v)
                        middleFace.add(v)
            newEdges.add(frozenset(newEdge))
            newEdge.add(vertex)
            newFaces.add(frozenset(newEdge))
            newFaces.add(frozenset(middleFace))
        
    Vertices = Vertices+newVertices
    Edges.clear()
    Edges.update(newEdges)
    Faces.clear()
    Faces.update(newFaces)

### Plotting ###

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set(xlim=(-1,1), ylim=(-1,1), zlim=(-1,1))
ax.add_collection3d(Poly3DCollection([ [ Vertices[i] for i in f] for f in Faces],edgecolors='k',facecolors='w'))
plt.axis('off')
plt.show()

Answer 1

I believe the problem to be on line 59, Please see my inline comments and images below:

        for vertex in face:
            incidentEdges = {edge for edge in SpecialEdges if vertex in edge}
            newEdge = set() # this wasn't checked if incidentEdges was empty
            for edge in incidentEdges:
                for v in edge:
                    if v == vertex:
                        continue
                    else:
                        newEdge.add(v)
                        middleFace.add(v)
            if not newEdge: continue # added this to prevent an empty frozenset being added later on
            newEdges.add(frozenset(newEdge)) # this was adding the empty edge
            newEdge.add(vertex)
            newFaces.add(frozenset(newEdge)) # this was creating an empty frozen set which causes it to fail later
            newFaces.add(frozenset(middleFace))

My change produces this for iteration 3 and 4.