program triangle c find triangulations given #vertices and Euler characteristic. implicit none integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge, . Vertices, Nvertices, V,E,F, Nb_edges, orient, . first_edge, set_connect logical dupl_edges, want_orientable common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), . from(100), NF, NV, NE, vert, edge, dupl_edges, . Vertices (100,100), Nvertices (100), V,E,F, Nb_edges, . want_orientable, orient (100), first_edge (100,2), . set_connect character*1 choice integer EC, I(100), J(100), K, L integer cnt c function definitions logical checkface, trim_tree, orientable, sym_check vert = 1 !tags that keep track of your path in the tree. edge = 2 write (6,*) 'Number of vertices?' read (5,*) V write (6,*) 'Euler characteristic?' read (5,*) EC write (6,*) 'Orientable? Y or N' read (5,'(a)') choice if (choice.eq.'Y'.or.choice.eq.'y') then want_orientable = .true. write (6,*) 'Will look for an orientable manifold' write (6,*) else want_orientable = .false. write (6,*) 'Will look for a non-orientable manifold' write (6,*) end if write (6,*) 'Fixed # vertices connected to a given vertex?' write (6,*) ' 0 means no check for this.' read (5,*) set_connect write (6,*) if (V.le.EC) stop 'No go' if (V.lt.4) stop 'No go' F = (V-EC)*2 E = (V-EC)*3 NF = 1 !also is level of recursion NV = 3 simp (1,1) = 1 simp (1,2) = 2 simp (1,3) = 3 10 continue c write (7,*) 'Simplex ', cnt+1 c if (cnt+1.eq.1999) then c write (7,*) ((simp(K,L), L = 1,3), K = 1, NF-1) c pause c end if call boundary !sets Edges c !array & finds #edges for c !Euler char. cnt = cnt + 1 if (mod(cnt, 1 000 000).eq.0) write (6,*) 'cnt is ', cnt c write (7,*) 'nf,ne,nv,nedges(nf), nvertices(nf), nb_edges' c write (7,*) nf,ne,nv,nedges(nf), nvertices(nf), nb_edges c write (7,'(a, 2i3)') 'first_edge ', (first_edge (NF, K), K = 1, 2) c write (7,'(a,20(i3))') ' Simplex', (simp(nf,K), K=1,3) c write (7,'(a,20(2i3,2x))') c . ' Edges', ((edges(nf, K,L), L=1,2), K=1, nedges(nf)) c write (7,'(a,20(i3))') c . ' Vertices', (Vertices(nf, K), K=1, nvertices(nf)) c write (7,*) c Does the structure violate symmetry conditions? if (.not.sym_check()) then NF = NF - 1 if (from(NF+1).eq.Vert) NV = NV - 1 if (NF.eq.0) then write (6,*) cnt, ' nodes of the tree were looked at' stop 'No triangulation' end if goto 20 !back in recursion end if if ((.not.dupl_edges).and. . NF.eq.F.and.NE.eq.E.and.NV.eq.V.and.Nedges(NF).eq.0) then if (.not.want_orientable) then orient(1) = 1 do K = 2, NF orient (K) = 0 end do if (orientable()) then NF = NF - 1 if (from(NF+1).eq.Vert) NV = NV - 1 goto 20 !back in recursion end if end if write (6,*) 'Triangulation:' do K = 1, NF write (6,'(3i3)') (simp (K, L), L= 1, 3) end do write ( 6,*) cnt, ' nodes of the tree were looked at' stop end if c if (trim_tree()) then NF = NF - 1 if (from(NF+1).eq.Vert) NV = NV - 1 if (NF.eq.0) then write (6,*) cnt, ' nodes of the tree were looked at' stop 'No triangulation' end if goto 20 !back one recursion level end if c c Loop through edges list and vertices list, make new simplexes. J(NF) = 0 !loop index on vertices to connect to 21 continue !top of loop J(NF) = J(NF) + 1 if (J(NF).gt.Nvertices(NF)) goto 29 !out of loop if (checkface (first_Edge(NF, 1), . first_EDGE (NF, 2), . Vertices(NF, J(NF)))) goto 10 c !down one level of recursion c checkface makes new simplex if everything is OK. 20 continue !returning from recursion goto 21 !bottom of J(NF) loop 29 continue !exit from J(NF) loop if (NF.gt.1) then !go back to node of tree you came from NF = NF - 1 if (from(NF+1).eq.Vert) NV = NV - 1 goto 20 !back one recursion level end if c write (6,*) 'No triangulation' write (6,*) cnt, ' nodes of the tree were looked at' stop end c subroutine boundary implicit none integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge, . Vertices, Nvertices, V,E,F, Nb_edges, orient, . first_edge, set_connect logical dupl_edges, want_orientable common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), . from(100), NF, NV, NE, vert, edge, dupl_edges, . Vertices (100,100), Nvertices (100), V,E,F, Nb_edges, . want_orientable, orient (100), first_edge (100,2), . set_connect integer t1, t2, t3, I, J, code_ij integer code_edge (0:175*175) dupl_edges = .false. NE = 0 Nb_edges = 0 do I = 1, NV*NV code_edge(I) = 0 end do do I = 1, NF t1 = simp(I,1) -1 + NV*(simp(I,2)-1) !store edges in !base NV code. code_edge ( t1) = code_edge (t1) + 1 !simplexes have vertices c !in ascending order t2 = simp(I,1) -1 + NV*(simp(I,3)-1) code_edge (t2) = code_edge (t2) + 1 t3 = simp(I,2) -1 + NV*(simp(I,3)-1) code_edge (t3) = code_edge (t3) + 1 end do Nedges(NF) = 0 do I = 1, NV*NV if (code_edge(I).gt.2) then dupl_edges = .true. return end if if (code_edge(I).gt.0) NE = NE + 1 if (mod (code_edge (I), 2).ne.0) then Nedges(NF) = Nedges(NF) + 1 Edges (NF, Nedges(NF), 1) = mod (I, NV) + 1 Edges (NF, Nedges(NF), 2) = I - Edges (NF, Nedges(NF), 1) + 1 Edges (NF, Nedges(NF), 2) = Edges (NF, Nedges(NF), 2)/ NV + 1 end if end do Nvertices (NF) = 0 !list of vertices on the edge do I = 1, 100 Vertices(NF, I) = 0 end do do I = 1, Nedges (NF) do J = 1,2 Vertices(NF, Edges(NF, I,J)) = . Vertices (NF, Edges(NF, I,J)) + 1 end do end do do I = 1, NV if (Vertices (NF, I).eq.1) . write (6,*) 'Sheesh! Vertex is at only one edge, nf is', nf if (Vertices (NF,I).gt.0) then Nvertices (NF) = Nvertices (NF) + 1 Vertices (NF,Nvertices(NF)) = I end if end do first_edge (NF, 1) = Vertices (NF, 1) first_edge (NF, 2) = V do I = 1, Nedges (NF) if (Edges (NF, I, 1).eq.Vertices (NF, 1)) . first_edge (NF, 2) = min (first_edge (NF, 2), Edges (NF, I, 2)) !find the "smallest" edge end do c Now check to see how many edges are between vertices on the c boundary.for trimming the tree later. do I = 1, Nvertices (NF) do J = I+1, Nvertices (NF) code_ij = Vertices(NF,I)-1 + . NV * (Vertices(NF,J)-1) if (Code_edge(code_ij).ne.0) nb_edges = nb_edges + 1 end do end do if (NV.lt.V) then Nvertices (NF) = Nvertices (NF) + 1 Vertices (NF,Nvertices(NF)) = NV + 1 end if return end logical function Checkface (v1,v2, v3) implicit none integer v1, v2, v3 integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge, . Vertices, Nvertices, V,E,F, Nb_edges, orient, . first_edge, set_connect logical dupl_edges, want_orientable common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), . from(100), NF, NV, NE, vert, edge, dupl_edges, . Vertices (100,100), Nvertices (100), V,E,F, Nb_edges, . want_orientable, orient (100), first_edge (100,2), . set_connect integer s1, s2, s3, I, J, t1, t2 c a function. logical pinch checkface = .false. if (v3.eq.v1.or.v3.eq.v2) return c if (v3.lt.v1) then !sort vertices s1 = v3 s2 = v1 s3 = v2 elseif (v3.lt.v2) then s1 = v1 s2 = v3 s3 = v2 else s1 = v1 s2 = v2 s3 = v3 end if c Is there already a face for these three vertices? do I = 1, NF if (simp(I,1).eq.s1.and.simp(I,2).eq.s2.and.simp(I,3).eq.s3) . return end do c c check to see if we'd be making a pinchpoint if (pinch(v1,v2,v3)) then c write (7,*) 'Pinchpoint' c write (7,'(a,100(2i3,2x))') c . 'Edges are', ((Edges(NF, I, J), J = 1, 2), I = 1, Nedges(NF)) c write (7,*) 'Pinchpoint at', v2 c write (7,*) 'Edge & vertex that''re being tossed out', v1,v2,v3 c write (7,*) return end if if (pinch(v2,v1,v3)) then c write (7,*) 'Pinchpoint' c write (7,'(a,100(2i3,2x))') c . 'Edges are', ((Edges(NF, I, J), J = 1, 2), I = 1, Nedges(NF)) c write (7,*) 'Pinchpoint at', v1 c write (7,*) 'Edge & vertex that''re being tossed out', v1,v2,v3 c write (7,*) return end if if (pinch(v1,v3,v2)) then c write (7,*) 'Pinchpoint' c write (7,'(a,100(2i3,2x))') c . 'Edges are', ((Edges(NF, I, J), J = 1, 2), I = 1, Nedges(NF)) c write (7,*) 'Pinchpoint at', v3 c write (7,*) 'Edge & vertex that''re being tossed out', v1,v2,v3 c write (7,*) return end if c make a new face. checkface = .true. NF = NF + 1 simp (NF, 1) = s1 simp (NF, 2) = s2 simp (NF, 3) = s3 if (v3.eq.NV+1) then NV = NV+1 from (NF) = vert else from (NF) = edge end if return end logical function trim_tree () implicit none integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge, . Vertices, Nvertices, V,E,F, Nb_edges, orient, . first_edge, set_connect logical dupl_edges, want_orientable common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), . from(100), NF, NV, NE, vert, edge, dupl_edges, . Vertices (100,100), Nvertices (100), V,E,F, Nb_edges, . want_orientable, orient (100), first_edge (100,2), . set_connect integer S1, S2, S3, S4, DV, DE, DF, Dedges c function definition. logical orientable trim_tree = .true. if (NF.eq.F) return if (dupl_edges) return if (nedges(nf).eq.0) return c Now. There are 4 basic triangling operatings: c S1: add vertex: NF -> NF+1, NE -> NE+2, NV -> NV+1, Nedges -> c Nedges+1 c S2: add edge: NF -> NF+1, NE -> NE+2, Nedges -> Nedges+1 c S3: fill in boundary: NF -> NF+1, NE -> NE+1, Nedges -> Nedges-1 c S4: fill in a triangle: NF -> NF+1, Nedges -> Nedges-3 c c This means that we can decide whether the current list of simplexes c can lead to the desired triangulation. DV = V-NV DE = E-NE DF = F-NF if (abs(DV-DE+DF).gt.DF) return !Euler characteristic too much off Dedges = - Nedges(NF) S1 = DV !# of S1 operations do S2 = 0, max (FLOAT(DF-S1), (DE-2*S1)/2.0)+1 !try values for S2 S3 = DE - 2*S1 - 2*S2 if (S3.lt.0) goto 9 S4 = DF - S1 - S2 - S3 if (S4.lt.1) goto 9 !must have at least one of these. goto 10 9 continue end do return !trim the tree 10 continue c Now, are there enough exposed vertices left over to give you the c necessary extra edges & faces? if (Vertices(NF, Nvertices(NF)).eq.NV+1) then DV = DV + Nvertices(NF) - 1 !# of vertices on the boundary else !plus unused ones. DV = DV + Nvertices(NF) end if if ((DV*(DV-1)/2) + NE - Nb_edges.lt.E) return !can only make !new edges from !these vertices if (want_orientable.and..not.orientable()) return trim_tree = .false. return end logical function Pinch (v1,v2, v3) implicit none integer v1, v2, v3 integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge, . Vertices, Nvertices, V,E,F, Nb_edges, orient, . first_edge, set_connect logical dupl_edges, want_orientable common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), . from(100), NF, NV, NE, vert, edge, dupl_edges, . Vertices (100,100), Nvertices (100), V,E,F, Nb_edges, . want_orientable, orient (100), first_edge (100,2), . set_connect integer s1,s2,s3, I, match, edge_cnt integer edge_list(175,2) logical edge_tag(175) c This function tells if v2 would be a pinch point with this new c triangulation. pinch = .false. c is (v2,v3) an edge on the boundary? s2 = min (v2,v3) s3 = max (v2,v3) do I = 1, Nedges (NF) if (Edges(NF,I,1).eq.s2.and.Edges(NF,I,2).eq.s3) goto 10 end do return ! (v2,v3) is not on the boundary 10 continue c do I = 1, Nedges (NF) if (Edges (NF,I,1).eq.v2.and. . (Edges(NF,I,2).ne.v1.and.Edges(NF,I,2).ne.v3)) goto 15 !keep on !looking. if (Edges (NF,I,2).eq.v2.and. . (Edges(NF,I,1).ne.v1.and.Edges(NF,I,1).ne.v3)) goto 15 end do 14 return !no pinch point 15 continue !keep on looking c is (v1,v2) an edge on the boundary? s1 = min (v1,v2) s2 = max (v1,v2) do I = 1, Nedges (NF) if (Edges(NF,I,1).eq.s1.and.Edges(NF,I,2).eq.s2) goto 20 end do return ! (v1,v2) is not on the boundary 20 continue c c Now are we forming a cycle of simplexes around v2? edge_cnt = 0 do I = 1, NF if (simp(I,1).eq.v2) then edge_cnt = edge_cnt + 1 edge_list (edge_cnt,1) = simp(I,2) edge_list (edge_cnt,2) = simp(I,3) elseif (simp(I,2).eq.v2) then edge_cnt = edge_cnt + 1 edge_list (edge_cnt,1) = simp(I,1) edge_list (edge_cnt,2) = simp(I,3) elseif (simp(I,3).eq.v2) then edge_cnt = edge_cnt + 1 edge_list (edge_cnt,1) = simp(I,1) edge_list (edge_cnt,2) = simp(I,2) end if end do do I = 1, edge_cnt edge_tag(I) = .true. end do c Now see if there's a cycle around v2. Starting with v1, can you get to c v3? match = v1 30 continue do I = 1, edge_cnt if (edge_tag(I)) then if (edge_list(I,1).eq.match) then match = edge_list (I,2) edge_tag(I) = .false. if (match.eq.v3) then !has to be a pinch point, we !know some single edges with !v2 remain. pinch = .true. return end if goto 30 elseif (edge_list(I,2).eq.match) then match = edge_list (I,1) edge_tag(I) = .false. if (match.eq.v3) then pinch = .true. return end if goto 30 end if end if end do return !no pinch point is being made end logical function orientable () implicit none integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge, . Vertices, Nvertices, V,E,F, Nb_edges, orient, . first_edge, set_connect logical dupl_edges, want_orientable common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), . from(100), NF, NV, NE, vert, edge, dupl_edges, . Vertices (100,100), Nvertices (100), V,E,F, Nb_edges, . want_orientable, orient (100), first_edge (100,2), . set_connect integer bdry (100*100) integer t1, t2, t3, I, J orientable = .true. if (NF.eq.1) then orient (NF) = 1 return end if orient (NF) = 0 !unset this to reset in this function. c first sum up the simplexes for which we already have orientation. do I = 1, NV*NV bdry (I) = 0 end do do I = 1, NF if (orient(I).eq.0) goto 9 !end of simplexes for c !which orientation was c !found t1 = simp(I,1) -1 + NV*(simp(I,2)-1) !store edges in !base NV code. bdry ( t1) = bdry (t1) + orient(I) !simplexes have vertices c !in ascending order t2 = simp(I,1) -1 + NV*(simp(I,3)-1) bdry (t2) = bdry (t2) - orient(I) t3 = simp(I,2) -1 + NV*(simp(I,3)-1) bdry (t3) = bdry (t3) + orient(I) end do 9 continue do J = I, NF if (orient(J).eq.0) then c Now add on the J'th simplex t1 = simp(J,1) -1 + NV*(simp(J,2)-1) t2 = simp(J,1) -1 + NV*(simp(J,3)-1) t3 = simp(J,2) -1 + NV*(simp(J,3)-1) if (bdry(t1).eq.1) then orient (J) = -1 elseif(bdry(t1).eq.-1) then orient (J) = 1 elseif (bdry (t2).eq.1) then orient (J) = 1 elseif (bdry (t2).eq.-1) then orient (J) = -1 elseif (bdry (t3).eq.1) then orient (J) = -1 elseif (bdry (t3).eq.-1) then orient (J) = 1 end if if (orient(J).ne.0) then bdry (t1) = bdry (t1) + orient (J) bdry (t2) = bdry (t2) - orient (J) bdry (t3) = bdry (t3) + orient (J) if (bdry (t1).eq.2.or.bdry(t1).eq.-2) then orientable = .false. return end if if (bdry (t2).eq.2.or.bdry(t2).eq.-2) then orientable = .false. return end if if (bdry (t3).eq.2.or.bdry(t3).eq.-2) then orientable = .false. return end if else write (6,*) 'Orientation wasn''t found ...' end if end if end do return end logical function sym_check () implicit none integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge, . Vertices, Nvertices, V,E,F, Nb_edges, orient, . first_edge, set_connect logical dupl_edges, want_orientable common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), . from(100), NF, NV, NE, vert, edge, dupl_edges, . Vertices (100,100), Nvertices (100), V,E,F, Nb_edges, . want_orientable, orient (100), first_edge (100,2), . set_connect integer connect (100), I, J logical on_edge (100) c Do symmetry checking on the triangulation. sym_check = .true. if (set_connect.ne.0) then !check to see each vertex isn't !connected to too many others. do I = 1, 100 connect (I) = 0 end do do I = 1, NF do J = 1,3 connect (Simp (I, J)) = connect (Simp (I, J)) + 1 end do end do c Set on_edge tag true if vertex is on the boundary do I = 1, NV on_edge (I) = .false. end do do I = 1, NVertices (NF) on_edge (Vertices (NF, I)) = .true. end do c If vertex is on the boundary, don't want set_connect faces already c attached to it. do I = 1, NV if (on_edge(I)) then if (connect(I).ge.set_connect) then sym_check = .false. return end if else c if vertex is not on the boundary, want set_connect faces attached. if (connect(I).ne.set_connect) then sym_check = .false. return end if end if end do end if return end