r""" Edgewise subdivision of a simplicial complex, implementing the algorithm of H. Edelsbrunner and D.R. Grayson (2000). Edgewise subdivision of a simplex. Discrete & Computational Geometry, 24(4), 707-719. Authors: - Fedor Manin: initial version (6/1/2018) """ def edgewise_subdivision(X,k): faces_list=X.faces() # start with the vertices of the original complex new_faces_list=[[((v[0],k),)] for v in faces_list[0]] # now generate every subdivided face, disjointly for d in faces_list: if d<1: continue the_tuples=[tuple(tup) for tup in Tuples([0 .. k],d+1).list() if sum(tup)==k] the_partitions=UnorderedTuples([1 .. k],d).list() subdivided_faces={} for dface in faces_list[d]: # Generate the vertices. # These are points in the faces whose barycentric coordinates are # integer multiples of 1/k. We represent them by said integers. face_vertices={} for tup in the_tuples: vertex_draft=zip(dface,tup) vertex_draft=[(v,i) for (v,i) in vertex_draft if i!=0] vertex_draft.sort() # the _values_ represent the vertex as a thing living in X # the _keys_ represent the position of the vertex in the face face_vertices[tup]=tuple(vertex_draft) # simplices correspond to pairs: # (list of breakpoints, permutation of [1 .. d]) # (not all of them, see later) for part in the_partitions: for perm in Permutations(d): a_flag=0 the_face=[] for i in [0 .. d]: # generate the ith vertex of the face as a list of # breakpoints which partition k vertex_part=[] for j in [1 .. d]: # curse 0- vs. 1-indexing if perm[j-1] <= i: vertex_part.append(part[j-1]-1) else: vertex_part.append(part[j-1]) # convert the list into barycentric coordinates; # if any of the coordinates are negative, # it's not a real vertex... vertex_part=[0]+vertex_part+[k] new_vertex=[] for j in [0 .. d]: the_diff=vertex_part[j+1]-vertex_part[j] if the_diff<0: a_flag=1 continue new_vertex.append(the_diff) if a_flag==1: continue the_face.append(face_vertices[tuple(new_vertex)]) #... and so the face is not a real face if a_flag==0: # print "subdivided face "+str(the_face) new_faces_list.append(the_face) return SimplicialComplex(new_faces_list)