@i common.w @* Navigating Simplicial Complexes. Only two types of simplices have explicit structures, vertices and tetrahedra. Tetrahedra have many fields giving characteristics of their edges and faces. Simplices are referred to by a |simplex_id|. Each vertex has a unique id. Each tetrahedron has a unique id. Each triangle has two possible ids, one for each containing tetrahedron. Each edge has many ids (on average about 6 ids), one for each containing tetrahedron. Although edges and triangles have more than one id, but there is one optimal id with which information about the simplex can be retrieved most quickly. Consequently, many functions replace a given id with the optimal id. @c #include "Morse.h" #include "globals.h" #include "simplex.wh" @@; @@; @ @(simplex.h@>+= typedef SInt32 simplex_id; @*1 Vertex structure. First we have the vertex structure. The field |h| is the the value of the vertex. I expect 16 bit resolution of the function is good enough, since there is deterministic tie breaking. The field |tetra_id| is the id of a tetrahedron containing the simplex. If the vertex is in $K$ and not critical, this tetrahedron contains the edge with which the vertex is paired. The field |type| contains various bit flags with the following meanings: \item{$\cdot$} Bit 0 is set if the vertex is in $K$. \item{$\cdot$} Bit 1 is set if the vertex is critical. \item{$\cdot$} If bit 1 is not set then bits 3,4,5 give the index of its paired edge in the tetrahedron |tetra_id|. @(simplex.h@>+= typedef struct vertexstruct { UInt8 type; SInt16 h; simplex_id tetra_id; GLfloat coords[3]; } vertex; @*1 Tetrahedron structure. Now we have the tetrahedron structure. The field |vertex_ids| gives the ids of the four vertices. These are in order of decreasing value. The field |face_ids| gives the alternate ids of the four faces of the tetrahedron, i.e., the ids based on the four neighboring tetrahedra. These are in an order consistent with |vertex_ids|, so |face_ids[i]| is the face which doen't contain |vertex_ids[i]|. Note however that when when first read in the fields |face_ids| have a very different meaning which we explain later in the section `Reading in a Complex'. \smallskip The field |flag| is used when reading in the complex. If |optional_link_values| is set then the lower order 2 bits give the index |i| of the vertex such that the quotient |(value(vertex_ids[0])-value(vertex_ids[i]))/|length of edge between |vertex_ids[0]| and |vertex_ids[i]| is as small as possible. The next 2 bits give the index with second smallest quotient and the third two give the index with the largest quotient. If |optional_link_values| is not set then flag is set to |0x39| which gives the indices defining the order of the values of |vertex_ids[i]|, which are always in decreasing order. Note that in either case, |vertex_ids[0]| is always the vertex with highest value and |vertex_ids[1]| the vertex with next highest value, and so on. \smallskip The field |types[0]| has bit fields with various info about the tetrahedron: \item{$\cdot$} Bit 0 is set if the tetrahedron is in $K$. \item{$\cdot$} Bit 1 is set if the tetrahedron is critical. \item{$\cdot$} Bits 3 and 4 give the index of the paired face if the tetrahedron is not critical. \item{$\cdot$} Bit 6 is set if it is known that no gradient paths from the tetrahedron go to critical triangles. \item{$\cdot$} Bit 7 is set if this tetrahedron has been processed. If not set many fields are untrustworthy. ? obsolete? \smallskip The bits in character |types[i]| have the folowing meaning for $1\le i\le 4$, refering to the triangle face opposite the $i-1$-th vertex: \item{$\cdot$} Bit 0 is set if the triangle is in $K$. \item{$\cdot$} Bit 1 is set if the triangle is critical. \item{$\cdot$} Bit 2 is set if the triangle is paired with an edge. \item{$\cdot$} Bits 3-4 give the relative index of that edge if bit 2 is set. \item{$\cdot$} Bit 5 is set if the triangle is paired down and it is known that all gradient paths to the triangle do not come from critical triangles. \item{$\cdot$} Bit 6 is set if the triangle is paired up and it is known that all gradient paths from the triangle do not lead to critical triangles. Bit 6 is also set if the triangle is paired down and it is known that all gradient paths from the triangle do not lead to critical edges. \item{$\cdot$} Bit 7 is set if this type character is valid. A given triangle in $K$ is valid in exactly one tetrahedron, and if it is paired to a tetrahedron it is valid in that tetrahedron. Bits 0 and 7 are correct in all two tetrahedra containing the triangle, the others only if bit 7 is set. Perhaps they should be valid in both tetrahedra. \smallskip The bits in character |types[i]| have the folowing meaning for $5\le i\le 10$, refering to the $i-5$-th edge: \item{$\cdot$} Bit 0 is set if the edge is in $K$. \item{$\cdot$} Bit 1 is set if the edge is critical. \item{$\cdot$} Bit 2 is set if the edge is paired with a vertex. \item{$\cdot$} Bit 3 gives the relative index of that vertex if bit 2 is set. (the i so the vertex is |get_vertex(e,i)|.) Bits 3 and 4 give the index of the paired triangle if bits 1 and 2 are not set. \item{$\cdot$} Bit 5 is set if the edge is paired up and it is known that all gradient paths to the edge do not come from critical triangles. \item{$\cdot$} Bit 6 is set if the edge is paired up and it is known that all gradient paths from the edge do not lead to critical edges. \item{$\cdot$} Bit 7 is set if this type character is valid. A given edge in $K$ is valid in exactly one tetrahedron. Bits 0 and 7 are correct in all tetrahedra containing the edge, the others only if bit 7 is set. However, bit 6, if set in some tetrahedron containing the edge, is valid. @(simplex.h@>+= typedef struct tetrahedronstruct { UInt8 types[11]; UInt8 flag; simplex_id vertex_ids[4]; simplex_id face_ids[4]; } tetrahedron; @*1 Converting to and from a simplex id. @(globals.c@>+= #include "Morse.h" int vtab[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}}; int etab[4][3] = {{10,9,8}, {10,7,6}, {9,7,5}, {8,6,5}}; int etab2[4][4] = {{0,5,6,7}, {5,0,8,9}, {6,8,0,10}, {7,9,10,0}}; int id_dim_tab[16] = {3,2,2,2,2,1,1,1,1,1,1,0,3,-2,-2,-1}; int id_val_tab[11] = {0,1,0,0,0,0,0,0,1,1,2}; int best_edge_tab[6] = {3,3,2,3,2,1}; UInt32 number_of_vertices, number_of_tetrahedra; vertex *vertexlist; tetrahedron *tlist; list *crit[4]; @ @(globals.h@>+= extern int vtab[6][2]; extern int etab[4][3]; extern int etab2[4][4]; extern int id_dim_tab[16]; extern int id_val_tab[11]; extern int best_edge_tab[6]; extern UInt32 number_of_vertices, number_of_tetrahedra; extern vertex *vertexlist; extern tetrahedron *tlist; extern list *crit[4]; @ @(simplex.h@>+= #define id2vertex(id) (vertexlist+((id)>>4)) #define vertex_id(n) ((simplex_id)(((n)<<4) + 11)) #define id2tetrahedron(id) (tlist +((id)>>4)) #define id2tetra_id(id) ((id)&0xfffffff0) #define tetrahedron_id(n) ((simplex_id)((n)<<4)) #define id_is_null(id) (((id)&15)==15) #define NULL_ID 15 #define id_dimension(id) (id_dim_tab[(id)&15]) #define edge_index(id) (((id)&15)-5) #define face_index(id) (((id)&15)-1) void get_edge_vertices(simplex_id e,simplex_id *vl); void get_triangle_vertices(simplex_id f, simplex_id *vlist); int triangle_in_tetrahedron(simplex_id f,simplex_id t); list *vertex_star(simplex_id v); @*1 Getting faces of a simplex. @(simplex.h@>+= #define get_vertex(e,i) (id2tetrahedron(e)->vertex_ids[vtab[ edge_index(e)][i]]) #define get_tetrahedron_vertices(t,vlist)@/ \ memcpy(vlist,id2tetrahedron(t)->vertex_ids,4*sizeof(simplex_id)) simplex_id other_vertex_in_edge(simplex_id u,simplex_id e); @ @+= simplex_id other_vertex_in_edge(simplex_id u,simplex_id e) { tetrahedron *t; simplex_id v; int j = edge_index(e); t = id2tetrahedron(e); v = t->vertex_ids[vtab[j][0]]; if (u==v) v = t->vertex_ids[vtab[j][1]]; return v; } void get_edge_vertices(simplex_id e,simplex_id *vl) { tetrahedron *tp; int j = edge_index(e); tp = id2tetrahedron(e); vl[0] = tp->vertex_ids[vtab[j][0]]; vl[1] = tp->vertex_ids[vtab[j][1]]; } void get_triangle_vertices(simplex_id f,simplex_id *vlist) // return a list of the vertices of |f| { int k = face_index(f); simplex_id *q = vlist; simplex_id *p = id2tetrahedron(f)->vertex_ids; if(k!=0) *q++ = *p; p++; if(k!=1) *q++ = *p; p++; if(k!=2) *q++ = *p; p++; if(k!=3) *q++ = *p; } @ |get_edge(f,i)| returns an id of the |i|-th edge of the triangle |f|. |get_face(t,i)| returns an id of the i-th face of the terahedron |t|. |get_edge_of_tetrahedron(s,i,j)| returns the edge whose vertices are the i-th and j-th vertices of the tetrahedron containing |s|. Usually |s| is a tetrahedron, but it could be an edge or triangle. |coface(f,i)| returns a tetrahedron containing the triangle |f|. |coface0(f) = coface(f,0)| and |coface1(f) = coface(f,1)|. @(simplex.h@>+= #define get_edge(f,i) (id2tetra_id(f)+etab[face_index(f)][i]) #define get_face(t,i) (id2tetra_id(t)+(i)+1) #define get_edge_of_tetrahedron(s,i,j) (id2tetra_id(s)+etab2[i][j]) #define coface0(f) id2tetra_id(f) #define coface1(f) id2tetra_id(id2tetrahedron(f)->face_ids[face_index(f)]) #define coface(f,i) (((i)==0)? coface0(f): coface1(f)) #define other_face_id(f) id2tetrahedron(f)->face_ids[face_index(f)] simplex_id other_coface(simplex_id s,simplex_id t); simplex_id other_edge(simplex_id v,simplex_id e,simplex_id t); @ @+= simplex_id other_coface(simplex_id s,simplex_id t) { if(id_dimension(s)!=2 || id_dimension(t)!=3) abort_message("internal coface"); if(id2tetra_id(s) != t) return id2tetra_id(s); return id2tetra_id(id2tetrahedron(t)->face_ids[face_index(s)]); } @ Given a tetrahedron |t| and a list |fvl| of three vertices of |t|, |vertices2triangle| and |vertices2trianglex| return the id of the corresponding triangle. Use |vertices2trianglex| if the list of vertices |fvl| might not be in order. Likewise, |vertices2edge| returns the id of the edge corresponding to a list of two vertices (in order). @+= simplex_id vertices2triangle(simplex_id *fvl,simplex_id t) { simplex_id *tvl = id2tetrahedron(t)->vertex_ids; if(fvl[0]==tvl[0] && fvl[1]==tvl[1] && fvl[2]==tvl[2]) return get_face(t,3); if(fvl[0]==tvl[0] && fvl[1]==tvl[1] && fvl[2]==tvl[3]) return get_face(t,2); if(fvl[0]==tvl[0] && fvl[1]==tvl[2] && fvl[2]==tvl[3]) return get_face(t,1); if(fvl[0]==tvl[1] && fvl[1]==tvl[2] && fvl[2]==tvl[3]) return get_face(t,0); return NULL_ID; } simplex_id vertices2trianglex(simplex_id *vl,simplex_id t) { simplex_id *tvl = id2tetrahedron(t)->vertex_ids; simplex_id fvl[3]; if(vl[0]>vl[1]) { fvl[0]=vl[0]; fvl[1]=vl[1]; } else { fvl[0]=vl[1]; fvl[1]=vl[0]; } if(vl[2]vertex_ids; int i,j; for(i=0;i<3;i++) { if(fvl[0]==tvl[i]) { for(j=i+1;j<4;j++) { if(fvl[1]==tvl[j]) return get_edge_of_tetrahedron(t,i,j); } } } return NULL_ID; } @ These routines take a codimension one face of a simplex and find the unique other codimension one face which also contains the given codimension 2 face. @+= simplex_id other_edge(simplex_id v,simplex_id e,simplex_id t) // return the edge of the triangle t which contains v but is not e { simplex_id evl[2],tvl[3]; get_edge_vertices(e,evl); get_triangle_vertices(t,tvl); if(evl[0]==tvl[0] && evl[1]==tvl[1]) return (v==evl[0])? get_edge(t,1): get_edge(t,0); if(evl[0]==tvl[0] && evl[1]==tvl[2]) return (v==evl[0])? get_edge(t,2): get_edge(t,0); if(evl[0]==tvl[1] && evl[1]==tvl[2]) return (v==evl[0])? get_edge(t,2): get_edge(t,1); abort_message("edge not in triangle"); } simplex_id other_face(simplex_id e,simplex_id f,simplex_id t) // Find the face of t which contains the edge e but is not f { simplex_id evl[2],fvl[3]; simplex_id *tvl = id2tetrahedron(t)->vertex_ids; int i,j; get_edge_vertices(e,evl); get_triangle_vertices(f,fvl); if(evl[0]==fvl[0] && evl[1]==fvl[1]) i=2; else if(evl[0]==fvl[0] && evl[1]==fvl[2]) i=1; else if(evl[0]==fvl[1] && evl[1]==fvl[2]) i=0; else abort_message("edge not in triangle"); for(j=0;j<4;j++) { if(tvl[j]==fvl[i]) return get_face(t,j); } abort_message("triangle not in tetrahedron"); } @ Convert the dege id |e| to an id based on the tetrhedron |t|. @+= simplex_id edge_id_in_tetrahedron(simplex_id e,simplex_id t) { simplex_id evl[2]; simplex_id *tvl = id2tetrahedron(t)->vertex_ids; get_edge_vertices(e,evl); if(evl[0]==tvl[0]) { if(evl[1]==tvl[1]) return t+5; if(evl[1]==tvl[2]) return t+6; if(evl[1]==tvl[3]) return t+7; abort_message("internal error 23"); } else if(evl[0]==tvl[1]) { if(evl[1]==tvl[2]) return t+8; if(evl[1]==tvl[3]) return t+9; abort_message("internal error 24"); } return t + 10; } @ On input, |vert| is three vertices of the tetrahedron |t|. (or |t| could be a face of some tetrahedron which is implied). It must have |vert[0] > vert[1]|. On output, |vert[2]| is changed to the vertex of |t| which is not in |vert|. The return value is the unique neighboring tetrahedron |tp| which contains the vertices |vert[0]|, |vert[1]|, and the new |vert[2]|. This allows you to successively travel the tetrahedra containing a given edge. This is sort of equivalent to |other_coface(other_face(e,f,t), t)| where |e| is the edge with vertices |vert[0]| and |vert[1]|, and |f| is the triangle with vertices |vert[i]|. But it is a bit faster and returns a face id rather than a tetrahedron id. @+= simplex_id next_face(simplex_id *vert,simplex_id t) { simplex_id *v = id2tetrahedron(t)->vertex_ids; int i,j; @@; vert[2] = v[i]; return id2tetrahedron(t)->face_ids[j]; } @ On input, |vert| is three vertices of a tetrahedron |t| implied by the edge id |e|. It must have |vert[0] > vert[1]| and the two vertices of |e| are |vert[0]| and |vert[1]| On output, |vert[2]| is changed to the vertex of |t| which is not in |vert|. The return value is edge id of |e| based in the unique neighboring tetrahedron |tp| which contains |e|, and the new vertex |vert[2]|. This allows you to successively travel the tetrahedra containing a given edge. @+= simplex_id next_edge(simplex_id *vert,simplex_id e) { tetrahedron *tt = id2tetrahedron(e); simplex_id *v = tt->vertex_ids; simplex_id f; int i,j; @@; vert[2] = v[i]; f = tt->face_ids[j]; tt = id2tetrahedron(tt->face_ids[j]); v = tt->vertex_ids; if(v[0] == vert[0]) { if(v[1] == vert[1]) return get_edge_of_tetrahedron(f,0,1); if(v[2] == vert[1]) return get_edge_of_tetrahedron(f,0,2); if(v[3] == vert[1]) return get_edge_of_tetrahedron(f,0,3); abort_message("next tetrahedron error 5"); } if(v[1] == vert[0]) { if(v[2] == vert[1]) return get_edge_of_tetrahedron(f,1,2); if(v[3] == vert[1]) return get_edge_of_tetrahedron(f,1,3); abort_message("next tetrahedron error 6"); } if(v[2] != vert[0]) abort_message("next tetrahedron error 7"); if(v[3] != vert[1]) abort_message("next tetrahedron error 8"); return get_edge_of_tetrahedron(f,2,3); } @ @= if(v[0] == vert[0]) { if(v[1] == vert[2]) { j=1; if(v[2] == vert[1]) i=3; else i=2; } else if(v[2] == vert[2]) { j=2; if(v[1] == vert[1]) i=3; else i=1; } else if(v[3] == vert[2]) { j=3; if(v[2] == vert[1]) i=1; else i=2; } else abort_message("next tetrahedron error 1"); } else if(v[1] == vert[0]) { if(v[0] == vert[2]) { j=0; if(v[2] == vert[1]) i=3; else i=2; } else if(v[2] == vert[2]) { j=2; i = 0; } else if(v[3] == vert[2]) { j=3; i=0; } else abort_message("next tetrahedron error 2"); } else if(v[2] == vert[0]) { if(v[0] == vert[2]) { i=1; j=0; } else if(v[1] == vert[2]) { i=0; j=1; } else abort_message("next tetrahedron error 3"); } else abort_message("next tetrahedron error 4"); @ @+= #ifdef debug void print_simp(simplex_id s) { simplex_id vl[4]; switch(id_dimension(s)) { case 0: printf("vertex %ld",s>>4); break; case 1: get_edge_vertices(s,vl); printf("edge [%ld,%ld]",vl[0]>>4,vl[1]>>4); break; case 2: get_triangle_vertices(s,vl); printf("triangle [%ld,%ld,%ld]",vl[0]>>4,vl[1]>>4,vl[2]>>4); break; case 3: get_tetrahedron_vertices(s,vl); printf("tetrahedron [%ld,%ld,%ld,%ld]",vl[0]>>4,vl[1]>>4,vl[2]>>4,vl[3]>>4); break; default: abort_message("Unknown simplex dimension"); } } #endif @*1 Containment in the complex. @(simplex.h@>+= #define vertex_is_in_K(v) ((id2vertex(v)->type & 1)!=0) #define edge_is_in_K(e) ((id2tetrahedron(e)->types[(e)&15] & 1)!=0) #define triangle_is_in_K(f) ((id2tetrahedron(f)->types[(f)&15] & 1)!=0) #define tetrahedron_is_in_K(t) ((id2tetrahedron(t)->types[0] & 1)!=0) @*1 Simplex values. The simplices are stored in order so you can compare values just by comparing their ids. @(simplex.h@>+= #define vertex_value(v) ((SInt32)id2vertex(v)->h) #define vertex_compare(v,w) (v-w) #define smallest_vertex_in_edge(e) get_vertex(e,1) #define largest_vertex_in_edge(e) get_vertex(e,0) @ @+= SInt32 value(simplex_id t) { if((t&15)<11) return vertex_value(id2tetrahedron(t)->vertex_ids[id_val_tab[t&15]]); else if ((t&15)>11) abort_message("invalid id"); return vertex_value(t); } @ Find the smallest or largest valued vertex in a list of vertices of length n @+= int smallest_vertex(simplex_id *v,int n) { int i,j; for(i=0,j=1;j0) i=j; } return i; } @ This routine tests whether the tetrahedron |t| is in the lower Star of the vertex |v|. @(simplex.h@>+= #define is_in_lower_Star(t,v) ((v) == id2tetrahedron(t)->vertex_ids[0]) Boolean edges_equal(simplex_id e0,simplex_id e1); SInt32 min_value(simplex_id s,int d); @ find the maximum value of the vertices of the simplex |s| @+= SInt32 max_value(simplex_id s,int d) { SInt32 m; tetrahedron *tp; if(id_is_null(s)) return -0x8000; // smallest value switch(d) { case 3: tp = id2tetrahedron(s); m = value(tp->vertex_ids[0]); break; case 2: tp = id2tetrahedron(s); m = (face_index(s)==0)? value(tp->vertex_ids[1]): value(tp->vertex_ids[0]); break; case 1: m = value(get_vertex(s,0)); break; default: abort_message("max value of simplex of unknown dimension"); } return m; } @ find the minimum value of the vertices of the simplex |s| @+= SInt32 min_value(simplex_id s,int d) { SInt32 m; tetrahedron *tp; if(id_is_null(s)) return 0x7fff; // largest value switch(d) { case 3: tp = id2tetrahedron(s); m = value(tp->vertex_ids[3]); break; case 2: tp = id2tetrahedron(s); m = (face_index(s)==3)? value(tp->vertex_ids[2]): value(tp->vertex_ids[3]); break; case 1: m = value(get_vertex(s,1)); break; default: abort_message("max value of simplex of unknown dimension"); } return m; } @ Put the distance between vertices |v0| and |v1| in |length| and return |(value(v1)-value(v0))/length|. @+= double alternate_vertex_value(simplex_id v0,simplex_id v1,double *length) { GLfloat *u0,*u1,d,h; u0 = id2vertex(v0)->coords; u1 = id2vertex(v1)->coords; d = (u0[0]-u1[0])*(u0[0]-u1[0])+(u0[1]-u1[1])*(u0[1]-u1[1])+(u0[2]-u1[2])*(u0[2]-u1[2]); d = sqrt(d); if(length!=NULL) *length=d; h = value(v1)-value(v0); return h/d; } @ return the length of the cross product |(v1-v0)\times (v2-v0)|. This is the sine of the angle formed by |v1,v0,v2| times the distance from |v0| to |v1| times the distance from |v0| to |v2|. @+= double cross_length(simplex_id v0,simplex_id v1,simplex_id v2) { GLfloat *u2,*u0,*u1,c[3]; u2 = id2vertex(v2)->coords; u0 = id2vertex(v0)->coords; u1 = id2vertex(v1)->coords; c[0] = (u1[1]-u0[1])*(u2[2]-u0[2])-(u1[2]-u0[2])*(u2[1]-u0[1]); c[1] = (u1[2]-u0[2])*(u2[0]-u0[0])-(u1[0]-u0[0])*(u2[2]-u0[2]); c[2] = (u1[0]-u0[0])*(u2[1]-u0[1])-(u1[1]-u0[1])*(u2[0]-u0[0]); return sqrt(c[0]*c[0]+c[1]*c[1]+c[2]*c[2]); } @ return the alternate value of edge |e1| relative to edge |e0|. @+= double alternate_edge_value_basic(simplex_id e0,simplex_id e1) { double l0,l1,h0,h1,c; simplex_id v,v0,v1; v = get_vertex(e0,0); if(v!=get_vertex(e1,0)) abort_message("alternate_edge_value error"); v0 = get_vertex(e0,1); v1 = get_vertex(e1,1); h0 = alternate_vertex_value(v,v0,&l0); h1 = alternate_vertex_value(v,v1,&l1); c = cross_length(v,v0,v1); return (h1-h0)*l0*l1/c; } @ A speedier version of the above. Here |v1| is the lower vertex of the edge |e1| and the alternate vale and length of |e0| are a;ready known so they do not need to be recalculated. @+= double alternate_edge_value(simplex_id e0,simplex_id v1,double e0_value,double e0_length) { double l1,h1,c; simplex_id v,v0; v = get_vertex(e0,0); v0 = get_vertex(e0,1); h1 = alternate_vertex_value(v,v1,&l1); c = cross_length(v,v0,v1); return (h1-e0_value)*e0_length*l1/c; } @*1 Deadend simplices. @(simplex.h@>+= Boolean edge_is_deadend(simplex_id *e); Boolean edge_is_deadendx(simplex_id *e); Boolean triangle_is_deadend(simplex_id *t); Boolean tetrahedron_is_deadend(simplex_id t); void make_edge_deadend(simplex_id *e); void make_triangle_deadend(simplex_id *t); void make_tetrahedron_deadend(simplex_id t); #define edge_is_deadend(e) ((id2tetrahedron(*e)->types[(*e)&15] & 64)!=0) #define edge_is_deadendx(e) edge_is_deadend(e) @ @+= Boolean triangle_is_deadend(simplex_id *t) { tetrahedron *tp; *t = best_face_id(*t); tp = id2tetrahedron(*t); return (tp->types[(*t)&15] & 64)!=0; } Boolean tetrahedron_is_deadend(simplex_id t) { tetrahedron *tp; tp = id2tetrahedron(t); return (tp->types[0] & 64)!=0; } void make_edge_deadend(simplex_id *e) { simplex_id vlist[3]; simplex_id t,ee; t = id2tetra_id(*e); get_edge_vertices(*e,vlist); vlist[2] = id2tetrahedron(t)->vertex_ids[best_edge_tab[edge_index(*e)]]; ee = *e; do { id2tetrahedron(ee)->types[ee&15] |= 64; ee = next_edge(vlist,ee); } while(ee!=*e); } void make_triangle_deadend(simplex_id *t) { tetrahedron *tp; *t = best_face_id(*t); tp = id2tetrahedron(*t); tp->types[(*t)&15] |= 64; } void make_tetrahedron_deadend(simplex_id t) { tetrahedron *tp; tp = id2tetrahedron(t); tp->types[0] |= 64; } Boolean edge_is_xdeadend(simplex_id *e) { tetrahedron *tp; *e = best_edge_id(*e); tp = id2tetrahedron(*e); return (tp->types[(*e)&15] & 32)!=0; } Boolean triangle_is_xdeadend(simplex_id *t) { tetrahedron *tp; *t = best_face_id(*t); tp = id2tetrahedron(*t); return (tp->types[(*t)&15] & 32)!=0; } void make_edge_xdeadend(simplex_id *e) { tetrahedron *tp; tp = id2tetrahedron(*e); tp->types[(*e)&15] |= 32; *e = best_edge_id(*e); tp = id2tetrahedron(*e); tp->types[(*e)&15] |= 32; } void make_triangle_xdeadend(simplex_id *t) { tetrahedron *tp; *t = best_face_id(*t); tp = id2tetrahedron(*t); tp->types[(*t)&15] |= 32; } @*1 face containment tests. These routines test whether a given simplex is a codimension one face of another, and if it is, returns the index of that face. If it is not, then -1 is returned. @(simplex.h@>+= #define vertex_in_edge(v,e) (((v) == get_vertex(e,0))?0:(((v) == get_vertex(e,1))?1:-1)) @ @+= int edge_in_triangle(simplex_id e,simplex_id t) { simplex_id evl[2],tvl[3]; get_edge_vertices(e,evl); get_triangle_vertices(t,tvl); if(evl[0]==tvl[0] && evl[1]==tvl[1]) return 2; if(evl[0]==tvl[0] && evl[1]==tvl[2]) return 1; if(evl[0]==tvl[1] && evl[1]==tvl[2]) return 0; return -1; } int edge_verts_in_triangle(simplex_id *evl,simplex_id f) { simplex_id *p = id2tetrahedron(f)->vertex_ids; switch(face_index(f)) { case 0: if(evl[0]==p[1] && evl[1]==p[2]) return 2; if(evl[0]==p[1] && evl[1]==p[3]) return 1; if(evl[0]==p[2] && evl[1]==p[3]) return 0; break; case 1: if(evl[0]==p[0] && evl[1]==p[2]) return 2; if(evl[0]==p[0] && evl[1]==p[3]) return 1; if(evl[0]==p[2] && evl[1]==p[3]) return 0; break; case 2: if(evl[0]==p[0] && evl[1]==p[1]) return 2; if(evl[0]==p[0] && evl[1]==p[3]) return 1; if(evl[0]==p[1] && evl[1]==p[3]) return 0; break; case 3: if(evl[0]==p[0] && evl[1]==p[1]) return 2; if(evl[0]==p[0] && evl[1]==p[2]) return 1; if(evl[0]==p[1] && evl[1]==p[2]) return 0; break; } return -1; } int triangle_in_tetrahedron(simplex_id f,simplex_id t) { simplex_id fvl[3]; simplex_id *tvl = id2tetrahedron(t)->vertex_ids; if(t == id2tetra_id(f)) return(face_index(f)); get_triangle_vertices(f,fvl); if(fvl[0]==tvl[0] && fvl[1]==tvl[1] && fvl[2]==tvl[2]) return 3; if(fvl[0]==tvl[0] && fvl[1]==tvl[1] && fvl[2]==tvl[3]) return 2; if(fvl[0]==tvl[0] && fvl[1]==tvl[2] && fvl[2]==tvl[3]) return 1; if(fvl[0]==tvl[1] && fvl[1]==tvl[2] && fvl[2]==tvl[3]) return 0; return -1; } @*1 Simplex equality tests. These test whether two edges or triangles are equal. They assume no two different edges have the same vertices. If this were used on a CW complex where this might not be the case, these would need to be changed to determine the optimum ids and compare them. @+= Boolean edges_equal(simplex_id e0,simplex_id e1) { simplex_id v0[2],v1[2]; if(e0==e1) return true; if(id_is_null(e1) || id_is_null(e0)) return false; get_edge_vertices(e0,v0); get_edge_vertices(e1,v1); return (v0[0]==v1[0] && v0[1]==v1[1]); } Boolean triangles_equal(simplex_id f0,simplex_id f1) { tetrahedron *tp; simplex_id f; int n; if (f0==f1) return true; if(id_is_null(f1) || id_is_null(f0)) return false; tp = id2tetrahedron(f1); f = tp->face_ids[face_index(f1)]; return f0 == f; } @*1 Lists of simplices. Often there are lists of simplex ids. So the following are useful. @(simplex.h@>+= #define id_list_push(l,t) { simplex_id id_list_push_simplex = t; \ (void)list_push(l,&id_list_push_simplex); } simplex_id id_list_pop(list *l); simplex_id id_list_read(list *l); @ @+= simplex_id id_list_pop(list *l) { simplex_id t; list_pop(l,&t); return t; } simplex_id id_list_read(list *l) { simplex_id *r; r = (simplex_id *)list_read(l); if(r!=NULL) return *r; else return 0xffffffff; } @*1 Critical simplices. Here are various routines for critical simplices. In |unmake_critical|, the formerly critical simplex is not actually removed from the list of critical simplices. So when running through |crit| you may need to check if the entries are in fact critical. The routine |clean_crit| will delete all the noncritical entries from |crit|. @(simplex.h@>+= #define vertex_is_critical(v) ( (id2vertex(v)->type & 2)!=0 ) #define unmake_vertex_critical(v) ( id2vertex(v)->type &= 0x81 ) #define tetrahedron_is_critical(t) ( (id2tetrahedron(t)->types[0] & 2)!=0 ) #define unmake_tetrahedron_critical(t) ( id2tetrahedron(t)->types[0] &= 0x81 ) Boolean edge_is_critical(simplex_id *e); Boolean triangle_is_critical(simplex_id *t); void make_vertex_critical(simplex_id v); void make_edge_critical(simplex_id *e); void make_triangle_critical(simplex_id *f); void make_tetrahedron_critical(simplex_id t); void unmake_triangle_critical(simplex_id *f); void unmake_edge_critical(simplex_id *e); void clean_crit(void); @ @+= Boolean edge_is_critical(simplex_id *e) { tetrahedron *t; *e = best_edge_id(*e); t = id2tetrahedron(*e); return (id2tetrahedron(*e)->types[(*e)&15] & 2)!=0; } Boolean triangle_is_critical(simplex_id *t) { tetrahedron *tt; *t = best_face_id(*t); tt = id2tetrahedron(*t); return (id2tetrahedron(*t)->types[(*t)&15] & 2)!=0; } @ @+= void make_vertex_critical(simplex_id v) { @q+ printf("Making 0 ");@> @q+ print_simp(v);@> @q+ printf(" critical\n");@> id2vertex(v)->type |= 2; id_list_push(crit[0],v); } void make_edge_critical(simplex_id *e) { @q+ printf("Making 1 ");@> @q+ print_simp(*e);@> @q+ printf(" critical\n");@> *e = best_edge_id(*e); id2tetrahedron(*e)->types[(*e)&15] |= 2; list_push(crit[1],e); } void make_triangle_critical(simplex_id *f) { @q+ printf("Making 2 ");@> @q+ print_simp(*f);@> @q+ printf(" critical\n");@> *f = best_face_id(*f); id2tetrahedron(*f)->types[(*f)&15] |= 2; list_push(crit[2],f); } void make_tetrahedron_critical(simplex_id t) { @q+ printf("Making 3 ");@> @q+ print_simp(t);@> @q+ printf(" critical\n");@> id2tetrahedron(t)->types[(t)&15] |= 2; id_list_push(crit[3],t); } @ @+= void unmake_edge_critical(simplex_id *e) { *e = best_edge_id(*e); id2tetrahedron(*e)->types[(*e)&15] &= 0xe1; } void unmake_triangle_critical(simplex_id *f) { *f = best_face_id(*f); id2tetrahedron(*f)->types[(*f)&15] &= 0x81; } void clean_crit(void) { simplex_id *s; list *l; l = crit[0]; list_read_init(l); while ((s = list_read(l))!= NULL) { if(!vertex_is_critical(*s)) list_read_delete(l); } l = crit[1]; list_read_init(l); while ((s = list_read(l))!= NULL) { if(!edge_is_critical(s)) list_read_delete(l); } l = crit[2]; list_read_init(l); while ((s = list_read(l))!= NULL) { if(!triangle_is_critical(s)) list_read_delete(l); } l = crit[3]; list_read_init(l); while ((s = list_read(l))!= NULL) { if(!tetrahedron_is_critical(*s)) list_read_delete(l); } } void clean_crit_at(UInt32 at[4]) { simplex_id *s; list *l; l = crit[0]; list_read_init_at(l,at[0]); while ((s = list_read(l))!= NULL) { if(!vertex_is_critical(*s)) list_read_delete(l); } l = crit[1]; list_read_init_at(l,at[1]); while ((s = list_read(l))!= NULL) { if(!edge_is_critical(s)) list_read_delete(l); } l = crit[2]; list_read_init_at(l,at[2]); while ((s = list_read(l))!= NULL) { if(!triangle_is_critical(s)) list_read_delete(l); } l = crit[3]; list_read_init_at(l,at[3]); while ((s = list_read(l))!= NULL) { if(!tetrahedron_is_critical(*s)) list_read_delete(l); } } @*1 Making simplices valid. A given triangle has two possible ids, and a given edge has many possible ids, one for each tetrahedron which contains it. In each case we pick one of these ids to be the preferred one. We mark the corresponding tetrahedron valid to indicate this. The tetrahedron with the valid id will correctly indicate whether the edge or triangle is critical and if not, its paired simplex. @(simplex.h@>+= #define make_edge_valid(e) (id2tetrahedron(e)->types[(e)&15] |= 128) #define make_triangle_valid(f) (id2tetrahedron(f)->types[(f)&15] |= 128) #define edge_is_valid(e) ((id2tetrahedron(e)->types[(e)&15] & 128)!=0) #define triangle_is_valid(f) ((id2tetrahedron(f)->types[(f)&15] & 128)!=0) @ @+= void unmake_edge_valid(simplex_id *e) { *e = best_edge_id(*e); id2tetrahedron(*e)->types[(*e)&15] &= 0x7f; } void unmake_triangle_valid(simplex_id *f) { *f = best_face_id(*f); id2tetrahedron(*f)->types[(*f)&15] &= 0x7f; } // the following may be obsolete void mark_tetrahedron_done(simplex_id v,simplex_id t) { tetrahedron *tp; int i; tp = id2tetrahedron(t); for(i=0;i<4;i++) { if (v==tp->vertex_ids[i]) { tp->flag |= (1<vertex_ids[i]) { j = (tp->flag & (1<+= simplex_id best_edge_idx(simplex_id e); simplex_id best_face_id(simplex_id f); #define best_edge_id(e) (edge_is_valid(e)?e:best_edge_idx(e)) #define best_face_id(f) (triangle_is_valid(f)?f:other_face_id(f)) @ @+= simplex_id best_edge_idx(simplex_id e) { simplex_id vlist[3]; simplex_id ee; get_edge_vertices(e,vlist); vlist[2] = id2tetrahedron(e)->vertex_ids[best_edge_tab[edge_index(e)]]; for(ee = e;!edge_is_valid(ee);) { ee = next_edge(vlist,ee); if (ee==e) abort_message("internal error 25"); } return ee; } @*1 Pairing simplices. The routines |..._is_paired_up| find whether the simplex |t| is paired with a simplex of one higher dimension. Likewise |..._is_paired_down| find whether the simplex |t| is paired with a simplex of one lower dimension. The routines |rij| return the dimension |j| simplex paired with the dimension |i| simplex |t|. The routines |pairij| pair up the i simplex |s0| and the j simplex |s1|. @(simplex.h@>+= #define vertex_is_paired_up(v) ( (id2vertex(v)->type & 2)==0 ) #define tetrahedron_is_paired_down(t) ( (id2tetrahedron(t)->types[0] & 2)==0 ) #define r32(t) get_face(t,(id2tetrahedron(t)->types[0]>>3)&3) #define r01(v) ( id2vertex(v)->tetra_id + ((id2vertex(v)->type>>3)&7)+5 ) #define edge_is_paired_downx(e) (edge_is_valid(e)?(id2tetrahedron(e)->types[(e)&15]&6)==4:eipdx(e)) Boolean edge_is_paired_up(simplex_id *e); Boolean triangle_is_paired_up(simplex_id *f); Boolean edge_is_paired_down(simplex_id *e); Boolean triangle_is_paired_down(simplex_id *f); simplex_id r21(simplex_id *f); simplex_id r10(simplex_id *e); simplex_id r23(simplex_id *f); simplex_id r12(simplex_id *e); void pair01(simplex_id s0,simplex_id s1,Boolean flag); void pair12(simplex_id *s0,simplex_id s1,Boolean flag); void pair23(simplex_id *s0,simplex_id s1,Boolean flag); @ @+= Boolean edge_is_paired_up(simplex_id *e) { *e = best_edge_id(*e); return (id2tetrahedron(*e)->types[(*e)&15] & 6)==0; } Boolean triangle_is_paired_up(simplex_id *f) { *f = best_face_id(*f); return (id2tetrahedron(*f)->types[(*f)&15] & 6)==0; } Boolean edge_is_paired_down(simplex_id *e) { *e = best_edge_id(*e); return (id2tetrahedron(*e)->types[(*e)&15] & 6)==4; } Boolean eipdx(simplex_id e) { simplex_id vl[2],e1; get_edge_vertices(e,vl); if(vertex_is_paired_up(vl[0])) { e1 = r01(vl[0]); if(vl[1] == get_vertex(e1,1)) return true; } if(vertex_is_paired_up(vl[1])) { e1 = r01(vl[1]); if(vl[0] == get_vertex(e1,0)) return true; } return false; } Boolean triangle_is_paired_down(simplex_id *f) { *f = best_face_id(*f); return (id2tetrahedron(*f)->types[(*f)&15] & 6)==4; } simplex_id r21(simplex_id *f) { *f = best_face_id(*f); return get_edge(*f,(id2tetrahedron(*f)->types[(*f)&15]>>3)&3); } simplex_id r10(simplex_id *e) { *e = best_edge_id(*e); return get_vertex(*e,(id2tetrahedron(*e)->types[(*e)&15]>>3)&1); } simplex_id r23(simplex_id *f) { *f = best_face_id(*f); return id2tetra_id(*f); } simplex_id r12(simplex_id *e) { *e = best_edge_id(*e); return get_face(id2tetra_id(*e),(id2tetrahedron(*e)->types[(*e)&15]>>3)&3); } @ @+= void pair01(simplex_id s0,simplex_id s1,Boolean flag) { simplex_id s; @q+ printf("01 Pairing ");@> @q+ print_simp(s0);@> @q+ printf(" with ");@> @q+ print_simp(s1);@> @q+ printf("\n");@> if(flag) { s = s1; unmake_edge_valid(&s); } id2vertex(s0)->tetra_id = id2tetra_id(s1); id2vertex(s0)->type = (UInt8)(1+(edge_index(s1)<<3)); id2tetrahedron(s1)->types[s1&15] = (UInt8)((s0==get_vertex(s1,0))? 5+128: 5+128+8); } void pair12(simplex_id *s0,simplex_id s1,Boolean flag) { int n; tetrahedron *tp; simplex_id s; @q+ printf("12 Pairing ");@> @q+ print_simp(*s0);@> @q+ printf(" with ");@> @q+ print_simp(s1);@> @q+ printf("\n");@> if(flag) { s = s1; unmake_triangle_valid(&s); s = *s0; unmake_edge_valid(&s); } n = edge_in_triangle(*s0,s1); tp = id2tetrahedron(s1); tp->types[s1&15] = (UInt8)(5 + 128 + (n<<3)); *s0 = get_edge(s1,n); tp->types[(*s0)&15] = (UInt8)(1 + 128 + (face_index(s1)<<3)); } void pair23(simplex_id *s0,simplex_id s1,Boolean flag) { int n; tetrahedron *tp; simplex_id s; @q+ printf("23 Pairing ");@> @q+ print_simp(*s0);@> @q+ printf(" with ");@> @q+ print_simp(s1);@> @q+ printf("\n");@> if(flag) { s = *s0; unmake_triangle_valid(&s); } n = triangle_in_tetrahedron(*s0,s1); tp = id2tetrahedron(s1); *s0 = get_face(s1,n); tp->types[0] = (UInt8)(1 + 128 + (n<<3)); tp->types[(*s0)&15] = 1 + 128; } @*1 Links of simplices. @ This creates a list of tetrahedra |simplex_id|s in the star of a vertex |v|. The list has a vague ordering, first on the list are tetrahedra in the lower star of |v|, next come tetrahedra for which |v| is the second highest vertex (so some face(s) are in the lower star), next come tetrahedra for which |v| is the third highest vertex (so some edge(s) are in the lower star), and last come tetrahedra for which |v| is the lowest vertex. @+= list *vertex_star(simplex_id v) { hlist *l; list *tl; int i; vertex *vv; hlist_initialize(&l,sizeof(tetrahedron *),4,2,0); list_initialize(&tl,sizeof(simplex_id)); vv = id2vertex(v); if(id_is_null(vv->tetra_id)) @@; @@; for(i=0;i<4;i++) @@; hlist_abandon(&l); return tl; } @ @= { tetrahedron *tt; simplex_id t; SInt32 m[2]; Boolean flag; list_push(tl,&(vv->tetra_id)); while(!list_is_empty(tl)) { list_pop(tl,m+1); tt = id2tetrahedron(m[1]); m[1] = id2tetra_id(m[1]); for(i=0;i<4;i++) { if(tt->vertex_ids[i]==v) {m[0] = i; break;} } hlist_find_add(l,m,&tt,&flag); if(!flag) { for(i=0;i<4;i++) { if(tt->vertex_ids[i]!=v) { list_push(tl,tt->face_ids+i); } } } } } @ @= { SInt32 m[2]; m[0]=i; while(hlist_sub_match(l,m,NULL,1)) { list_push(tl,m+1); } } @ When the complex was read in, critical vertices did not have their |tetra_id| field filled in. This only becomes a problem now, so we remedy it. We could also do this more efficiently- when a tetrahedron is read in, set the |tetra_id| of its lowest vertex to 15 + the id of the tetrahedron. This will be interpreted as a |NULL_ID|, but is easily converted to a valid id. @= { UInt32 j; tetrahedron *tt; for(j=0;jvertex_ids[3]==v) { vv->tetra_id = tetrahedron_id(j); break; } } if(j==number_of_tetrahedra) abort_message("internal error: tetra_id"); } @* Checking Complexes. Check the discrete Morse function in the parts of $K$ in the lower stars of vertices $\ge w$ for consistency. Form lists of all their edges and triangles and check that each simplex is either critical or paired with a simplex which is paired to it. @+= void complex_check(simplex_id w) { hlist *edges,*faces; SInt32 i,j,k; simplex_id v,e,f,t,vlist[4]; hlist_initialize(&edges,sizeof(simplex_id),12713,2,0); hlist_initialize(&faces,sizeof(simplex_id),12713,3,0); for(i=(w>>4);i@; } for(i=0;i0||vlist[1]>=w)) @@; for(k=j+1;k<4;k++) { e = get_edge_of_tetrahedron(t,j,k); if(edge_is_in_K(e)&&vlist[j]>=w) @@; } } } hlist_abandon(&edges); hlist_abandon(&faces); } @ @= { if(!vertex_is_critical(v)) { e = r01(v); @@; @@; } } @ @= { e = best_edge_id(e); @@; if(edge_is_paired_up(&e)) { f = r12(&e); @@; @@; @@; } else if(edge_is_paired_down(&e)) { v = r10(&e); @@; } else if(!edge_is_critical(&e)) abort_message("nothing edge"); } @ @= { if(!edge_is_paired_down(&e)) abort_message("should be paired-down edge"); if(v!=r10(&e)) abort_message("nonsymmetric pairing"); } @ @= { if(!edge_is_paired_up(&e)) abort_message("should be paired-up edge"); if(f!=r12(&e)) abort_message("nonsymmetric pairing"); } @ @= { Boolean flag; SInt32 index; if(id_dimension(e)!=1) abort_message("edge not an edge"); if(!edge_is_valid(e)) abort_message("non-valid edge"); get_edge_vertices(e,vlist); index = hlist_find_add(edges,vlist,&e,&flag); if(flag) { if(e!= *((simplex_id *)hlist_entry(edges,index))) abort_message("edge has two valid ids"); } } @ @= { simplex_id t1; tetrahedron *ttt,*tttt; f = best_face_id(f); ttt = id2tetrahedron(f); tttt = id2tetrahedron(t); @@; if(triangle_is_paired_up(&f)) { t1 = r23(&f); @@; @@; } else if(triangle_is_paired_down(&f)) { e = r21(&f); @@; @@; @@; } else if(!triangle_is_critical(&f)) abort_message("nothing triangle"); } @ @= { if(!triangle_is_paired_down(&f)) abort_message("should be paired-down face"); if(e!=r21(&f)) abort_message("nonsymmetric pairing"); } @ @= { if(!triangle_is_paired_up(&f)) abort_message("should be paired-up face"); if(t1!=r23(&f)) abort_message("nonsymmetric pairing"); } @ @= { Boolean flag; SInt32 index; if(id_dimension(f)!=2) abort_message("face not a face"); if(!triangle_is_valid(f)) abort_message("non-valid face"); get_triangle_vertices(f,vlist); index = hlist_find_add(faces,vlist,&f,&flag); if(flag) { if(f!= *((simplex_id *)hlist_entry(faces,index))) abort_message("face has two valid ids"); } } @ @= { if(!tetrahedron_is_paired_down(t1)) abort_message("should be paired-down tetrahedron"); if(f!=r32(t1)) abort_message("nonsymmetric pairing"); } @* Constructing Complexes. @*1 Allocating vertices and tetrahedra. The following routines allocate a new simplex. The last parameter |flag| is $=0$ if the simplex is not in |K|, and is $!=0$ if it is a vertex in |K| or if it is a tetrahedron and you want it to be in |K| if and only if all its vertices are in |K|. @+= vertex *new_vertex(UInt32 n,SInt16 val,float x,float y,float z,Boolean flag) { vertex *v; v = id2vertex( vertex_id(n)); v->type=(UInt8)((flag)?1:0); // indicate in $K$ or not v->tetra_id=NULL_ID; v->h = val; v->coords[0]=x; v->coords[1]=y; v->coords[2]=z; return v; } @ @+= tetrahedron *new_tetrahedron(UInt32 n,UInt32 *vs,Boolean flag) { tetrahedron *t; int i; Boolean flagp; t = id2tetrahedron( tetrahedron_id(n)); t->flag = 0; flagp = flag; for(i=0;i<4;i++) { t->vertex_ids[i] = vertex_id(vs[i]); flagp = flagp && vertex_is_in_K(t->vertex_ids[i]); t->face_ids[i] = NULL_ID; } t->types[0]=(UInt8)((flagp)?1:0); @@; @@; @@; return t; } @ @= { simplex_id temp; int j; for(j=0;j<3;j++) { for(i=j+1;i<4;i++) { if(vertex_compare(t->vertex_ids[j],t->vertex_ids[i]) < 0) { temp = t->vertex_ids[j]; t->vertex_ids[j] = t->vertex_ids[i]; t->vertex_ids[i] = temp; } } } } @ @= { int j; for(i=0;i<4;i++) { flagp = true; for(j=0;j<4;j++) { if(i==j) continue; flagp = flagp && vertex_is_in_K(t->vertex_ids[j]); } t->types[1+i] = (UInt8)((flagp)?1:0); } } @ @= { int j; for(i=0;i<3;i++) { for(j=i+1;j<4;j++) { flagp = vertex_is_in_K(t->vertex_ids[i]) && vertex_is_in_K(t->vertex_ids[j]); t->types[etab2[i][j]] = (UInt8)((flagp)?1:0); } } } @*1 Reading in a complex. For now we suppose that the ambient complex is a manifold, later versions could add in tetrahedra as needed to make the complex a pseudo manifold if it is not already. We also assume that the complex $K$ of interest is almost full, in the sense that an edge or triangle is in $K$ if and only if all its vertices are in $K$. It would be a trivial matter to read in a file of exceptions to this. @(simplex.h@>+= void read_in_complex(FILE *dfv,FILE *dft,double min,double max,Boolean qdelaunay, double scalex,double scaley,double scalez, Boolean optional_link_values); @ @+= @@; void read_in_complex(FILE *dfv,FILE *dft,double min,double max,Boolean qdelaunay, double scalex,double scaley,double scalez, Boolean optional_link_values) { UInt32 i; UInt32 *translation; // permutation which sorts the vertices if(max<=min) abort_message("max is less than min"); @@; @@; free(translation); @@; } @ Read in and sort the vertices by value. @+= list *vertex_blocks; // blocks of 256 vertices int vertex_sort(const void *nn, const void *mm); @ @= { float actual_min,actual_max; UInt32 *translation1; list_initialize(&vertex_blocks,sizeof(float *)); number_of_vertices=0; @@; vertexlist = (vertex *)malloc( (number_of_vertices) * sizeof(vertex) ); if (vertexlist == NULL ) abort_message("Out of memory"); translation1 = (UInt32 *)malloc( (number_of_vertices) * sizeof(UInt32) ); if (translation1 == NULL ) abort_message("Out of memory"); for(i=0;i@; for(i=0;i@; list_abandon(&vertex_blocks); translation = (UInt32 *)malloc( (number_of_vertices) * sizeof(UInt32) ); if (translation == NULL ) abort_message("Out of memory"); // invert the permutation |translation1| to |translation| for(i=0;i= { float x,y,z,h; int n; int c; float *p; do { c = getc(dfv); if(c==EOF) break; if(c!='x') { (void)ungetc(c,dfv); n = fscanf(dfv,"%f %f %f %f\n",&x,&y,&z,&h); if(n!=4) abort_message("bad vertex file"); if(number_of_vertices==0) actual_min=actual_max=h; else if(hactual_max) actual_max=h; x/=scalex; y/=scaley; z/=scalez; } else { @@; x = y = z = h = PLUS_INFINITY; } @@; number_of_vertices++; } while(true); } @ @= { do { c = getc(dfv); } while(c!=EOF && c!= '\n' && c!= '\r'); } @ @= { if((number_of_vertices & 0xff)==0) { p = (float *)malloc(256*4*sizeof(float)); if(p==NULL) abort_message("out of memory"); list_push(vertex_blocks,&p); } *p++ = h; *p++ = x; *p++ = y; *p++ = z; } @ @= qsort(translation1, number_of_vertices, sizeof(UInt32),vertex_sort); @ @+= int vertex_sort(const void *nn, const void *mm) { float *p,*q; UInt32 n,m; n = *(UInt32 *)nn; m = *(UInt32 *)mm; p = *(float **)list_entry(vertex_blocks,n>>8); q = *(float **)list_entry(vertex_blocks,m>>8); p+=((n&0xff)*4); q+=((m&0xff)*4); if(*p<*q) return -1; if(*p>*q) return 1; return 0; } @ @= { float *p; SInt16 val; p = *(float **)list_entry(vertex_blocks,translation1[i]>>8); p+=((translation1[i]&0xff)*4); if(*pmax) val = 32767; else val = (SInt32)(-32768.0+65535.0*(*p - min)/(max-min)); new_vertex(i, val ,p[1],p[2],p[3], *p!=PLUS_INFINITY); } @ @= { @@; tlist = (tetrahedron *)malloc( number_of_tetrahedra * sizeof(tetrahedron) ); if (tlist == NULL ) abort_message("Out of memory"); @@; } @ @= { long vi[4]; int c; int n; if(qdelaunay) { fscanf(dft,"%li\n",vi); number_of_tetrahedra = vi[0]+1; } else { number_of_tetrahedra = 0; do { @@; if(n!=4) break; number_of_tetrahedra++; } while (true); rewind(dft); } } @ @= { long vi[4]; int c; UInt32 vn[4]; int j,n; tetrahedron *t; for(i=0;i@; for(j=0;j<4;j++) vn[j] = translation[vi[j]]; t = new_tetrahedron(i,vn,c!='x'); @@; @q+1 {@> @q+1 int j;@> @q+1 print_simp(tetrahedron_id(i));@> @q+1 printf(" ");@> @q+1 for(j=0;j<11;j++) printf("%d",t->types[j]);@> @q+1 printf("\n");@> @q+1 }@> } } @ @= { if(qdelaunay && i==0) { // this is a tetrahedron at infinity, not in the input file vi[0] = number_of_vertices - 4; vi[1] = number_of_vertices - 3; vi[2] = number_of_vertices - 2; vi[3] = number_of_vertices - 1; c = 'x'; n = 4; } else { c = getc(dft); if(c!='x') (void)ungetc(c,dft); n = fscanf(dft,"%li %li %li %li\n",vi,vi+1,vi+2,vi+3); } } @ When first read in, a tetrahedron's |face_ids| fields are set to very different values from their eventual use, with very different meanings as follows: For each edge with a nonempty lower star, there is one tetrahedron in its lower star whose fields are as follows: \item{} |face_ids[0]| is |NULL_ID| or the id of a tetrahedron with the same first vertex but with a second vertex of lower value. \item{} |face_ids[1]| is |NULL_ID| or the id of a tetrahedron with the same first vertex but with a second vertex of higher value. \item{} |face_ids[2]| is |NULL_ID| or the id of a tetrahedron with the same first and second vertex, i.e., another tetrahedron in the lower star of the edge. \item{} |face_ids[3]| is ignored. The first two give a binary tree of the locally three dimensional edges in the lower star of the first vertex. The third starts a linked list of the tetrahedra in the lower star of that edge. All the other tetrahedra in the lower star of the edge have fields as follows: \item{} |face_ids[i]| is ignored for $i=0,1,3$. \item{} |face_ids[2]| is |NULL_ID| or the id of a tetrahedron with the same first and second vertex. @= { simplex_id *tp = &(id2vertex(t->vertex_ids[0])->tetra_id); tetrahedron *tpt; SInt32 res; double best_val,tree_val; if(optional_link_values) { if(t->types[0]&1) @flag| to optional vertex order@>@; else { // not in K t->flag = 0x39; best_val = PLUS_INFINITY; } } else t->flag = 0x39; while(*tp != NULL_ID) { tpt = id2tetrahedron(*tp); if(tpt->vertex_ids[tpt->flag&3] == t->vertex_ids[t->flag&3]) break; if(optional_link_values) { tree_val = alternate_vertex_value(t->vertex_ids[0], tpt->vertex_ids[tpt->flag&3],NULL); if (tree_val < best_val) res = 1; else if (tree_val > best_val) res = -1; else res = t->vertex_ids[t->flag&3] - tpt->vertex_ids[tpt->flag&3]; } else res = vertex_compare(t->vertex_ids[1],tpt->vertex_ids[1]); if(res<0) tp = tpt->face_ids; else if(res>0) tp = tpt->face_ids+1; else abort_message("Huh?"); } if(*tp != NULL_ID) { tp = tpt->face_ids+2; t->face_ids[2] = *tp; } *tp = tetrahedron_id(i); } @ @flag| to optional vertex order@>= { double vals[3]; vals[0] = alternate_vertex_value(t->vertex_ids[0],t->vertex_ids[1],NULL); vals[1] = alternate_vertex_value(t->vertex_ids[0],t->vertex_ids[2],NULL); vals[2] = alternate_vertex_value(t->vertex_ids[0],t->vertex_ids[3],NULL); if(vals[0]>=vals[1]) { if(vals[0]>=vals[2]) { if(vals[1]>=vals[2]) t->flag = 0x39; else t->flag = 0x2d; } else t->flag = 0x27; } else { if(vals[1]>=vals[2]) { if(vals[0]>=vals[2]) t->flag = 0x36; else t->flag = 0x1e; } else t->flag = 0x1b; } best_val = vals[(t->flag&3)-1]; } @ This bit of code does two things at once. It links up all the tetrahedra with neighboring tetrahedra which share a face with it, and it puts a discrete Morse function on the complex (a procedure called |ExtractRaw| in \cite{KKM}). In particular it runs through all the vertices |v| in the complex, first those not in |K| and then those in |K| in order of value, largest value first. For each vertex |v|, it links up each tetrahedron in the lower star of |v| with any of its neighboring tetrahedra which were previously encontered. (The lower Star of |v| is all simplices for which |v| is the vertex of maximal value.) For any neighboring tetrahedron not yet encountered, the common face is put on either the list |unglued_faces| (if the face does not contain |v|) or |local_unglued_faces| (if the face contains |v| but not the second highest vertex of the tetrahedron) or on |more_local_unglued_faces| (if the face contains both |v| and the second highest vertex of the tetrahedron). For each vertex |v|, if |v| is in |K| it also finds an optimal discrete Morse function on the lower Star of |v|. It then takes one critical edge |beste| in the lower Star of |v|, makes it noncritical, and pairs it with |v|. The choice of |beste| dependes on |optional_link_values|. If |optional_link_values| is false, it chooses |beste| so that its other vertex |bestv| has value as small as possible. If |optional_link_values| is true, it chooses |beste| so that the quantity (|value|(|v|) - |value|(|bestv|))/|length|(|beste|) is as negative as possible. Thus the discrete gradient flow from |v| follows the path of steepest descent in some sense. It is the genuine steepest edge if |optional_link_values| is true. But if |optional_link_values| is false, it follows the edge which decreases the function the most. This may not be the steepest edge if it is relatively long. The next cancelling step cancels all pairs of critical simplices it can in the lower star so that in the end you get as few of them as possible. Note that if the lower Star of |v| is empty, you make the vertex |v| critical, since |v| is a local minimum. This also sets up lists |lone_faces|, |possible_lone_edges|, |local_done_edges|, |local_lone_edges|, and |components| which will be explained later. @= hlist *unglued_faces; // triangles we have only seen in one tetrahedron hlist *local_unglued_faces; // and those in the lower star of v hlist *lone_faces; // triangles isolated in the lower star of their maximal vertex hlist *local_lone_faces; // and those in the lower star of v hlist *possible_lone_edges; // edges perhaps isolated in the lower star of their maximal vertex struct star_component { simplex_id end_face[2]; // faces at boundary of lower link component simplex_id end_vert[2]; // lowest vertex in end face SInt32 index[2]; // indices of |end_vert[2]| in |more_local_unglued_faces| simplex_id bestv; // minimal vertex in lower link component double bestv_value; // value of mininal vertex }; @ @= { simplex_id v; simplex_id beste; // this will be the steepest edge from |v| simplex_id bestv; // the other vertex of |bestv| hlist *more_local_unglued_faces; // and those in the lower star of an edge hlist *local_done_edges; // edges we have already processed hlist *local_lone_edges; // edges perhaps isolated in the lower star of their maximal vertex list *components; // list of components of the deleted lower star of an edge. UInt32 critsave[4]; // save top of lists of critical points for(i=0;i<4;i++) list_initialize(crit+i,sizeof(simplex_id)); hlist_initialize(&unglued_faces,sizeof(simplex_id),12713,3,0); //12713 prime hlist_initialize(&local_unglued_faces,sizeof(simplex_id),103,2,0); hlist_initialize(&more_local_unglued_faces,sizeof(SInt32),13,1,0); hlist_initialize(&lone_faces,sizeof(simplex_id),2713,3,0); //2713 prime hlist_initialize(&local_lone_faces,sizeof(simplex_id),103,2,0); hlist_initialize(&local_done_edges,sizeof(simplex_id),103,1,0); hlist_initialize(&local_lone_edges,sizeof(simplex_id),103,1,0); hlist_initialize(&possible_lone_edges,sizeof(simplex_id),2713,2,0); //2713 prime list_initialize(&components,sizeof(struct star_component)); for(i = number_of_vertices; i>0 ;) { i--; v = vertex_id(i); @q+ printf("\n\nExtracting LowLk of vertex %ld at (%f,%f,%f) with value %d\n",@> @q+ i,id2vertex(v)->coords[0],id2vertex(v)->coords[1],id2vertex(v)->coords[2],@> @q+ id2vertex(v)->h);@> beste = NULL_ID; bestv = NULL_ID; hlist_clear(local_unglued_faces); hlist_clear(local_lone_faces); hlist_clear(local_done_edges); hlist_clear(local_lone_edges); critsave[0] = list_count(crit[0]); critsave[1] = list_count(crit[1]); critsave[2] = list_count(crit[2]); critsave[3] = list_count(crit[3]); @@; if(vertex_is_in_K(v)) { if(id_is_null(beste)) { make_vertex_critical(v); } else { pair01(v,beste,false); // pair |v| with the steepest edge down LocalCancel(v,critsave[2]); // cancel critical simplices in the lower Star of |v| clean_crit_at(critsave); // delete cancelled critical points } } } hlist_abandon(&unglued_faces); hlist_abandon(&local_unglued_faces); hlist_abandon(&more_local_unglued_faces); hlist_abandon(&lone_faces); hlist_abandon(&local_lone_faces); hlist_abandon(&local_done_edges); hlist_abandon(&local_lone_edges); hlist_abandon(&possible_lone_edges); list_abandon(&components); } @ Here we link up the tetrahedra in the lower star of |v| and if $v\in K$ we find a discrete Morse function on the lower star of |v|. We do this by taking each edge |e| in the lower Star of |v|, linking up the tetrahedra in its lower star, and (if |e|$\in K$) putting a discrete Morse function on its lower Star. The lower Star of |e| is all simplices $\sigma$ of $K$ containing |e| so that the maximal and next--to--maximal vertices of $\sigma$ are in |e|. On entry the |face_ids| give a tree structure on the lower star of |v|, whose root is in |vv->tetra_id|. On exit each |face_ids[i]| will be the id of a face shared with a neighboring tetrahedron (if the neighboring tetrahedron was already encountered) or |NULL_ID| (if not). There are three types of edges in the lower star of |v| which we must account for. The first type is those edges whose lower star is three dimensional. We process those edges first, in order (highest value of the other vertex first). The next type of edge has two dimensional lower star. We process these next in |@|. As we'll see later the triangles in the lower star of these edges are all in the lists |lone_faces| or |local_lone_faces|. Finally we have the edges with empty lower star. These are on the lists |local_lone_edges| or |possible_lone_edges| and are processed last. @= /* |(simplex_id v,*simplex_id beste,*simplex_id bestv, hlist *local_lone_edges,hlist *local_done_edges, Boolean optional_link_values, list *components,hlist *more_local_unglued_faces)| */ { vertex *vv=id2vertex(v); simplex_id edge; tetrahedron *tt; simplex_id t; double beste_value; while(!id_is_null(vv->tetra_id)) { @@; @@; } @@; @@; } @ @= { simplex_id *last; last = &vv->tetra_id; do { t = *last; tt = id2tetrahedron(t); if (id_is_null(tt->face_ids[1])) break; last = tt->face_ids + 1; } while(true); *last = tt->face_ids[0]; edge = get_edge_of_tetrahedron(t,0,tt->flag&3); } @ We have an edge named |edge| whose lower star is three dimensional and will link up all tetrahedra in its lower star. If |edge| is in |K| we will also put an optimum discrete Morse function on the lower star of |edge|. |face_with_lowest_vertex| is the triangle in the lower star of |edge| whose third vertex |lowest_vertex| has the lowest value. This triangle will be paired with |edge|. There are two possiblities for the lower star of |edge|. It could be that every tetrahedron containing |edge| is in the lower star. So the lower link is a circle. In this case we make one of the tetrahedra containing |face_with_highest_vertex| critical. There will be no other critical simplices in the lower star of |edge|. The other possibility is that the lower link consists of a number of intervals and points. Each of the intervals will be represented by an entry in the list |components|. Recall that an entry of |components| is a structure with fields |simplex_id end_face[2]|, |simplex_id end_vert[2]|, |SInt32 index[2]|, |simplex_id bestv|, and |double bestv_value|. The field |end_vert[2]| is the two vertices at the ends of the interval. The field |end_face[2]| is the corresponding triangles in the lower star, i.e., the joins of |end_vert[2]| with |edge|. The field |index[2]| is the indices in the list |more_local_unglued_faces| of the vertices |end_vert[2]|. The field |bestv| is the minimal vertex in the lower link and |bestv_value| is either the simplex_id itself, or in the alternate case, the alternate value |value(bestv)-value(v)/length(v,bestv)|. The isolated points in the lower link correspond to triangles in |local_unglued_faces| and |unglued_faces| whose highest edge is |edge|. @= /* |(tetrahedron *tt, simplex_id t, hlist *more_local_unglued_faces, list *components, simplex_id edge, Boolean optional_link_values, *simplex_id beste, *simplex_id bestv, *double beste_value, hlist *local_done_edges, hlist *local_lone_edges, simplex_id v)| */ { simplex_id next; simplex_id lowest_vertex = NULL_ID; simplex_id face_with_lowest_vertex = NULL_ID; double lowest_vertex_value; simplex_id highest_vertex = NULL_ID; simplex_id face_with_highest_vertex = NULL_ID; double highest_vertex_value; double edge_value,edge_length; if(optional_link_values) edge_value = alternate_vertex_value(v,get_vertex(edge,1),&edge_length); @q+ printf("\nExtracting LowLk of ");@> @q+ print_simp(edge);@> @q+ printf("\n");@> hlist_clear(more_local_unglued_faces); list_clear(components); do { next = id2tetrahedron(t)->face_ids[2]; @q+ printf("Processing ");@> @q+ print_simp(t);@> if(tetrahedron_is_in_K(t)) @@; else @@; t = next; } while(!id_is_null(t)); if(edge_is_in_K(edge)) @@; } @ Here we take a tetrahedron in the lower star of |edge| and add it to our representation of the lower star. Perhaps it is easier to visualize what this does to the lower link. Then we are adding an edge to the portion of the lower link we have already discovered, which is a collection of intervals. One of three things could happen: \item{1.} The new edge could intersect the portion of the link in both vertices. Thus it either joins two intervals together or makes an interval into a circle. \item{2.} The new edge could intersect in a single vertex, thus lengthening one of the intervals. \item{3.} The new edge could be disjoint from the portion of the link, thus creating a new interval. @= /* |(tetrahedron *tt, simplex_id t, hlist *more_local_unglued_faces, list *components, simplex_id next, Boolean optional_link_values, simplex_id edge, double edge_value, double edge_length, *double lowest_vertex_value, *double highest_vertex_value, *simplex_id lowest_vertex, *simplex_id face_with_lowest_vertex, *simplex_id highest_vertex, *simplex_id face_with_highest_vertex)| */ { SInt32 more_local_index[2]; SInt32 comp_index[2]; int vertex_index[2]; simplex_id f; int j; tt = id2tetrahedron(t); vertex_index[0] = ((tt->flag)&12)>>2; vertex_index[1] = ((tt->flag)&48)>>4; tt->face_ids[0] = NULL_ID; tt->face_ids[1] = NULL_ID; tt->face_ids[2] = NULL_ID; tt->face_ids[3] = NULL_ID; more_local_index[0] = hlist_add_delete(more_local_unglued_faces, tt->vertex_ids+vertex_index[0],comp_index); more_local_index[1] = hlist_add_delete(more_local_unglued_faces, tt->vertex_ids+vertex_index[1],comp_index+1); if(more_local_index[0]<0 && more_local_index[1]<0) @@; else if(more_local_index[0]>=0 && more_local_index[1]>=0) @@; else { if(more_local_index[0]>=0) j=1; else j=0; @@; } } @ This code either joins two intervals together or makes a single interval into a circle. The last possibility could only occur if |t| is the very last tetrahedron in the lower star and every tetrahedron containing |edge| is in the lower star of |edge|. In the first possibility, we combine the two components in the struct for the first component and abandon the struct for the last component, since I doubt whether the list of components would ever be long enough in practice to worry about garbage collection. The maximum number of items needed in the list |components| is the number of tetrahedra in the lower link of |edge| which would typically be small. @= /* |(list *components, SInt32 *comp_index, tetrahedron *tt, simplex_id t, int *vertex_index, Boolean optional_link_values, hlist *more_local_unglued_faces, simplex_id next)| */ { struct star_component *comp[2]; int j0,j1; simplex_id f0,f1; @q+ printf(" -joining\n");@> comp[0] = (struct star_component *)list_entry(components,comp_index[0]); comp[1] = (struct star_component *)list_entry(components,comp_index[1]); if(comp[0]->end_vert[0] == tt->vertex_ids[vertex_index[0]]) j0 = 0; else j0 = 1; if(comp[1]->end_vert[0] == tt->vertex_ids[vertex_index[1]]) j1 = 0; else j1 = 1; f0 = get_face(t,vertex_index[0]); f1 = comp[1]->end_face[j1]; glue_faces(f0,f1); f0 = get_face(t,vertex_index[1]); f1 = comp[0]->end_face[j0]; glue_faces(f0,f1); if(comp_index[0]==comp_index[1]) { // we are joining an interval into a circle. // so there should be no other simplices to add // and no other components in the lower star. if(!id_is_null(next) || !hlist_is_empty(more_local_unglued_faces)) abort_message("link of an edge is disconnected"); } else { comp[0]->end_vert[j0] = comp[1]->end_vert[1-j1]; comp[0]->end_face[j0] = comp[1]->end_face[1-j1]; comp[0]->index[j0] = comp[1]->index[1-j1]; if(comp[0]->bestv_value>comp[1]->bestv_value) { comp[0]->bestv = comp[1]->bestv; comp[0]->bestv_value = comp[1]->bestv_value; } *((SInt32 *)hlist_entry(more_local_unglued_faces,comp[0]->index[j0])) = comp_index[0]; } } @ This adds in a new component to the lower link of |edge|. @= /* |(tetrahedron *tt, simplex_id t, SInt32 *more_local_index, *simplex_id lowest_vertex, *simplex_id face_with_lowest_vertex, *simplex_id highest_vertex, *simplex_id face_with_highest_vertex, *double lowest_vertex_value, *double highest_vertex_value, Boolean optional_link_values, int *vertex_index, int j, simplex_id edge, double edge_value, double edge_length, hlist *more_local_unglued_faces, list *components)| */ { struct star_component comp; double val[2]; @q+ printf(" -isolated\n");@> if(optional_link_values) { val[0] = alternate_edge_value(edge,tt->vertex_ids[vertex_index[0]], edge_value,edge_length); val[1] = alternate_edge_value(edge,tt->vertex_ids[vertex_index[1]], edge_value,edge_length); } else { val[0] = tt->vertex_ids[vertex_index[0]]; val[1] = tt->vertex_ids[vertex_index[1]]; } if(val[0] > val[1]) j=1; else j=0; comp.end_vert[0] = tt->vertex_ids[vertex_index[0]]; comp.end_vert[1] = tt->vertex_ids[vertex_index[1]]; comp.end_face[0] = get_face(t,vertex_index[1]); comp.end_face[1] = get_face(t,vertex_index[0]); comp.index[0] = more_local_index[0]; comp.index[1] = more_local_index[1]; comp.bestv = tt->vertex_ids[vertex_index[j]]; comp.bestv_value = val[j]; if(id_is_null(lowest_vertex) || lowest_vertex_value > val[j]) { lowest_vertex = tt->vertex_ids[vertex_index[j]]; face_with_lowest_vertex = get_face(t,vertex_index[1-j]); lowest_vertex_value = val[j]; } if(id_is_null(highest_vertex) || highest_vertex_value < val[1-j]) { highest_vertex_value = val[1-j]; highest_vertex = tt->vertex_ids[vertex_index[1-j]]; face_with_highest_vertex = get_face(t,vertex_index[j]); } *((SInt32 *)hlist_entry(more_local_unglued_faces,more_local_index[0])) = (SInt32)list_count(components); *((SInt32 *)hlist_entry(more_local_unglued_faces,more_local_index[1])) = (SInt32)list_count(components); list_push(components,&comp); } @ Here the tetrahedron |t| intersects the already discovered portion of the link in just a single face, |get_face(t,3-j)|. @= /* |(list *components, SInt32 *comp_index, tetrahedron *tt, int j, simplex_id t, Boolean optional_link_values, int *vertex_index, *simplex_id face_with_lowest_vertex, *simplex_id lowest_vertex, *simplex_id highest_vertex, *simplex_id face_with_highest_vertex, *double lowest_vertex_value, *double highest_vertex_value, simplex_id edge, double edge_value, double edge_length, hlist *more_local_unglued_faces, SInt32 *more_local_index)| */ { struct star_component *comp; int jj; simplex_id f0,f1; simplex_id new_vertex; simplex_id glued_vertex; double new_vertex_value; @q+ printf(" -merging\n");@> glued_vertex = tt->vertex_ids[vertex_index[j]]; new_vertex = tt->vertex_ids[vertex_index[1-j]]; if(optional_link_values) { new_vertex_value = alternate_edge_value(edge,new_vertex,edge_value,edge_length); } else new_vertex_value = new_vertex; comp = (struct star_component *)list_entry(components,comp_index[j]); if(comp->end_vert[0] == glued_vertex) jj = 0; else jj = 1; f0 = get_face(t,vertex_index[1-j]); f1 = comp->end_face[jj]; glue_faces(f0,f1); comp->end_vert[jj] = new_vertex; comp->end_face[jj] = get_face(t,vertex_index[j]); comp->index[jj] = more_local_index[1-j]; if(comp->bestv_value > new_vertex_value) { comp->bestv = new_vertex; comp->bestv_value = new_vertex_value; if(lowest_vertex_value > new_vertex_value) { lowest_vertex = new_vertex; lowest_vertex_value = new_vertex_value; face_with_lowest_vertex = get_face(t,vertex_index[j]); } } if(highest_vertex_value < new_vertex_value) { highest_vertex = new_vertex; highest_vertex_value = new_vertex_value; face_with_highest_vertex = get_face(t,vertex_index[j]); } *((SInt32 *)hlist_entry(more_local_unglued_faces,comp->index[jj])) = comp_index[j]; } @ At this point we have linked up all the tetrahedra in the lower star of |edge| to themselves. We have also linked up all the tetrahedra not in |K| which have highest edge equal to |edge|. @= /* |(simplex_id t, *simplex_id lowest_vertex, hlist *more_local_unglued_faces, list *components, *simplex_id face_with_highest_vertex, tetrahedron *tt, *simplex_id beste, *simplex_id bestv, *double beste_value, *double lowest_vertex_value, simplex_id edge, double edge_value, double edge_length, simplex_id face_with_lowest_vertex, hlist *local_done_edges, simplex_id v, Boolean optional_link_values, hlist *local_lone_edges)| */ { simplex_id w; SInt32 index; struct star_component *comp; simplex_id f,f0,f1,local_lowest_vertex; if(hlist_is_empty(more_local_unglued_faces) && !id_is_null(face_with_lowest_vertex)) @@; else { while(hlist_get(more_local_unglued_faces,&w,&index)) { comp = (struct star_component *)list_entry(components,index); if(comp->end_vert[0] == w) w = comp->end_vert[1]; else w = comp->end_vert[0]; if(!hlist_find_delete(more_local_unglued_faces,&w,NULL)) abort_message("internal error 46"); f = comp->end_face[0]; @@; f = comp->end_face[1]; @@; local_lowest_vertex = comp->bestv; f0 = comp->end_face[0]; f1 = comp->end_face[1]; @@; } @@; } @@; } @ @= /* |(*simplex_id beste,*simplex_id bestv,simplex_id edge,simplex_id w, Boolean optional_link_values, simplex_id v, *double beste_value, simplex_id face_with_lowest_vertex,hlist *local_done_edges)| */ { w = get_vertex(edge,1); if(id_is_null(face_with_lowest_vertex)) @@; else { pair12(&edge,face_with_lowest_vertex,false); } // indicate this edge has been looked at hlist_find_add(local_done_edges,&w,&edge,NULL); } @ @= /* |(*simplex_id beste,*simplex_id bestv,simplex_id edge, *double beste_value, Boolean optional_link_values, simplex_id v, simplex_id w)| */ { double w_value; if(optional_link_values) w_value = alternate_vertex_value(v,w,NULL); else w_value = w; if(id_is_null(bestv)) { bestv = w; beste = edge; beste_value = w_value; } else { if(beste_value>w_value) { make_edge_valid(beste); make_edge_critical(&beste); bestv = w; beste = edge; beste_value = w_value; } else { make_edge_valid(edge); make_edge_critical(&edge); } } } @ @= /* |(simplex_id f, simplex_id t, simplex_id lowest_vertex, simplex_id f0, simplex_id f1, simplex_id local_lowest_vertex, simplex_id edge, Boolean optional_link_values, simplex_id face_with_highest_vertex)| */ { int flag; // We should really look at the other coface and pick the optimal one. t = id2tetra_id(face_with_highest_vertex); f = get_face(t,0); @@; flag = id2tetrahedron(t)->flag; f = get_face(t,flag&3); @@; make_tetrahedron_critical(t); f0 = other_face_id(get_face(t,(flag>>2)&3)); f1 = other_face_id(get_face(t,(flag>>4)&3)); local_lowest_vertex = lowest_vertex; @@; } @ @= /* |(simplex_id f,*simplex_id lowest_vertex,simplex_id edge, *double lowest_vertex_value, double edge_value, double edge_length, *simplex_id face_with_lowest_vertex, Boolean optional_link_values, tetrahedron *tt, simplex_id v, hlist *local_lone_edges, hlist *local_done_edges)| */ { simplex_id vlist[3]; get_edge_vertices(edge,vlist); @q+ printf("Isolated triangles:\n");@> while(hlist_sub_match(lone_faces,vlist,&f,2)) @@; while(hlist_sub_match(local_lone_faces,vlist+1,&f,1)) @@; } @ @= /* |(simplex_id f,simplex_id *vlist,*simplex_id lowest_vertex, *double lowest_vertex_value, *simplex_id face_with_lowest_vertex, Boolean optional_link_values, simplex_id edge, double edge_value, double edge_length, tetrahedron *tt, simplex_id v, hlist *local_lone_edges, hlist *local_done_edges)| */ { simplex_id e; double new_vertex_value; if(triangle_is_in_K(f)) { @q+ printf("\t");@> @q+ print_simp(f);@> @q+ printf("\n");@> if(optional_link_values) { new_vertex_value = alternate_edge_value(edge,vlist[2],edge_value,edge_length); } else new_vertex_value = vlist[2]; if(id_is_null(lowest_vertex)) { lowest_vertex = vlist[2]; face_with_lowest_vertex = f; lowest_vertex_value = new_vertex_value; } else { if(lowest_vertex_value > new_vertex_value) { make_triangle_valid(face_with_lowest_vertex); make_triangle_critical(&face_with_lowest_vertex); lowest_vertex = vlist[2]; lowest_vertex_value = new_vertex_value; face_with_lowest_vertex = f; } else { make_triangle_valid(f); make_triangle_critical(&f); } } } e = get_edge(f,1); if(get_vertex(e,1)!=vlist[2]) e = get_edge(f,2); if(get_vertex(e,1)!=vlist[2]) abort_message("internal 59"); @@; e = get_edge(f,0); @@; } @ @= /* |( tetrahedron *tt, simplex_id e, simplex_id v, Boolean optional_link_values, hlist *local_lone_edges, hlist *local_done_edges)| */ { simplex_id vl[2]; if(edge_is_in_K(e)) { get_edge_vertices(e,vl); @q+ printf("\tadding to local_lone_edges ");@> @q+ print_simp(e);@> @q+ printf("\n");@> if(hlist_find(local_done_edges,vl+1)<0) hlist_find_add(local_lone_edges,vl+1,&e,NULL); } } @ @= /* |( tetrahedron *tt, Boolean optional_link_values, simplex_id e)| */ { simplex_id vl[2], *last; if(edge_is_in_K(e)) { get_edge_vertices(e,vl); @q+ printf("\tadding to possible_lone_edges ");@> @q+ print_simp(e);@> @q+ printf("\n");@> hlist_find_add(possible_lone_edges,vl,&e,NULL); } } @ @= /* |(simplex_id v,simplex_id edge,hlist *local_done_edges, Boolean optional_link_values, *double beste_value, tetrahedron *tt,*simplex_id beste,*simplex_id bestv, hlist *local_lone_edges)| */ { simplex_id vlist[3]; simplex_id lowest_vertex; double lowest_vertex_value; simplex_id face_with_lowest_vertex; simplex_id f,w; double edge_value,edge_length; vlist[0] = v; if(!hlist_is_empty(local_unglued_faces)) abort_message("internal 34"); if(hlist_sub_match(unglued_faces,vlist,NULL,1)) abort_message("internal 35"); do { if(hlist_sub_match(local_lone_faces,vlist+1,&f,0)) ; else if(!hlist_sub_match(lone_faces,vlist,&f,1)) break; if(hlist_find(local_done_edges,vlist+1)>=0) continue; edge = vertices2edge(vlist,id2tetra_id(f)); if(optional_link_values) edge_value = alternate_vertex_value(v,get_vertex(edge,1),&edge_length); @q+ printf("Extracting 2D Low Link of ");@> @q+ print_simp(edge);@> @q+ printf("\n");@> lowest_vertex = NULL_ID; face_with_lowest_vertex = NULL_ID; do { @@; if(!hlist_sub_match(lone_faces,vlist,&f,2)) { if(!hlist_sub_match(local_lone_faces,vlist+1,&f,1)) break; } } while(true); if(edge_is_in_K(edge)) @@; } while(true); } @ Here we add any edges whose lower star is empty. They will all be on either of the lists |local_lone_edges| or |possible_lone_edges|. But those lists may contain duplicates or edges whose lower star is not in fact empty, so we test them against the list |local_done_edges| to make sure. Let me see why these edges are on |local_lone_edges| or |possible_lone_edges|. An edge is added to these lists in the code segment |@|, when |getedge(f,1)| is added to |local_lone_edges| and |getedge(f,0)| is added to |possible_lone_edges|. The code segment |@| is executed for faces |f| which appear on the lists |lone_faces| and |local_lone_faces|. @= /* |(simplex_id v,*simplex_id beste,*simplex_id bestv,simplex_id edge, Boolean optional_link_values, *double beste_value, hlist *local_lone_edges,hlist *local_done_edges)| */ { simplex_id w,vlist[2]; vlist[0] = v; do { if(hlist_sub_match(local_lone_edges,vlist+1,&edge,0)) ; else if(!hlist_sub_match(possible_lone_edges,vlist,&edge,1)) break; if(!edge_is_in_K(edge) || hlist_find(local_done_edges,vlist+1)>=0) continue; @q+ printf("Extracting 1D Low Link of ");@> @q+ print_simp(edge);@> @q+ printf("\n");@> w = vlist[1]; @@; // indicate this edge has been looked at hlist_find_add(local_done_edges,&w,&edge,NULL); } while(true); } @ @= /* |(Boolean optional_link_values, simplex_id t)| */ { simplex_id f; @q+ printf(" -outside K\n");@> f = get_face(t,0); @@; f = get_face(t,1); @@; f = get_face(t,2); @@; f = get_face(t,3); @@; } @ @+= void glue_faces(simplex_id f0,simplex_id f1) { simplex_id v0[3]; simplex_id v1[3]; if(id2tetra_id(f0)==id2tetra_id(f1)) abort_message("err101"); get_triangle_vertices(f0,v0); get_triangle_vertices(f1,v1); if(v0[0]!=v1[0]) abort_message("err102"); if(v0[1]!=v1[1]) abort_message("err102"); if(v0[2]!=v1[2]) abort_message("err102"); glue_face(f0,f1); glue_face(f1,f0); } @ @d glue_face(f0,f1) (id2tetrahedron(f1)->face_ids[face_index(f1)] = f0) @d glue_faces_old(f0,f1) (glue_face(f0,f1), glue_face(f1,f0)) @= // |(Boolean optional_link_values, simplex_id f)| { simplex_id vl[3]; simplex_id f2 = f; get_triangle_vertices(f,vl); if(hlist_add_delete(unglued_faces,vl,&f2) < 0) { if(face_index(f2)!=0) abort_message("internal error 68"); glue_faces(f,f2); @@; hlist_add(lone_faces,vl,&f2); } } @ @= { double val0,val1; simplex_id temp; if(optional_link_values) { val0 = alternate_vertex_value(vl[0],vl[1],NULL); val1 = alternate_vertex_value(vl[0],vl[2],NULL); if(val1 > val0) { temp = vl[2]; vl[2] = vl[1]; vl[1] = temp; } } } @ @= /* |(Boolean optional_link_values, simplex_id f)| */ { simplex_id vl[3]; simplex_id f2 = f; get_triangle_vertices(f,vl); if(hlist_find_delete(unglued_faces,vl,&f2)) { if(face_index(f2)!=0) abort_message("internal error 68"); glue_faces(f,f2); @ hlist_add(local_lone_faces,vl+1,&f2); } else if(hlist_add_delete(local_unglued_faces,vl+1,&f2) < 0) { glue_faces(f,f2); @ hlist_add(local_lone_faces,vl+1,&f2); } } @ @= // |(Boolean optional_link_values, simplex_id f)| { simplex_id vl[3]; simplex_id f2 = f; get_triangle_vertices(f,vl); if(hlist_find_delete(unglued_faces,vl,&f2)) { if(face_index(f2)!=0) abort_message("internal error 68"); glue_faces(f,f2); } else if(hlist_add_delete(local_unglued_faces,vl+1,&f2) < 0) { glue_faces(f,f2); } else abort_message("link faces: I don't think this should happen"); } @ This takes an interval in the lower link of an edge whose minimum vertex is |local_min| and whose ends are the triangles |f0| and |f1| (after joining with the edge). It puts a discrete Morse function on the join of this interval with the edge. The vertex |global_min| is the minimum vertex in the whole lower link of the edge. @= { simplex_id vlist[3]; get_triangle_vertices(f0,vlist); @@; t = id2tetra_id(f0); f=f0; @@; get_triangle_vertices(f1,vlist); @@; t = id2tetra_id(f1); f=f1; @@; if(local_lowest_vertex != lowest_vertex) { make_triangle_valid(f); make_triangle_critical(&f); } } @ @= if(vlist[2] == get_vertex(edge,1)) { vlist[2] = vlist[1]; vlist[1] = get_vertex(edge,1); } @ @= while(vlist[2]!=local_lowest_vertex) { pair23(&f,t,false); f = get_face(t,0); @@; f = get_face(t,id2tetrahedron(t)->flag&3); @@; f = next_face(vlist,t); t = id2tetra_id(f); }