[974] | 1 | import numpy as np |
---|
| 2 | from scvt.timer import Timer |
---|
| 3 | |
---|
| 4 | ######################################################## |
---|
| 5 | |
---|
| 6 | def simplex(ii,jj,kk): |
---|
| 7 | if ii>jj: ii,jj = jj,ii |
---|
| 8 | if jj>kk: jj,kk = kk,jj |
---|
| 9 | if ii>jj: ii,jj = jj,ii |
---|
| 10 | # now ii<jj<kk |
---|
| 11 | return (ii,jj,kk) |
---|
| 12 | |
---|
| 13 | def enum_edges(simplices): |
---|
| 14 | for ii,jj,kk in simplices: |
---|
| 15 | # ii<jj<kk |
---|
| 16 | yield (ii,jj) |
---|
| 17 | yield (jj,kk) |
---|
| 18 | yield (ii,kk) |
---|
| 19 | |
---|
| 20 | def insert_edge(edges, regions, edge,kk,n): |
---|
| 21 | if edge in edges: |
---|
| 22 | edges[edge]=(edges[edge], (n,kk)) |
---|
| 23 | else: |
---|
| 24 | edges[edge]=(n,kk) |
---|
| 25 | ii,jj = edge |
---|
| 26 | if ii in regions: |
---|
| 27 | regions[ii].add(edge) |
---|
| 28 | else: |
---|
| 29 | regions[ii]={edge} |
---|
| 30 | if jj in regions: |
---|
| 31 | regions[jj].add(edge) |
---|
| 32 | else: |
---|
| 33 | regions[jj]={edge} |
---|
| 34 | |
---|
| 35 | def all_edges_regions(simplices): |
---|
| 36 | # returns edges, regions |
---|
| 37 | # edge = ordered pair of generators |
---|
| 38 | # edges is a dictionary edge -> ((n1,point1, (n2,point2)) |
---|
| 39 | # where n is the index of a simplex containing edge |
---|
| 40 | # and point is the index of the third point of that simplex |
---|
| 41 | # regions is a dictionary : point -> dict of edges |
---|
| 42 | edges,regions = {},{} |
---|
| 43 | for n,(ii,jj,kk) in enumerate(simplices): |
---|
| 44 | # note that ii<jj<kk |
---|
| 45 | # edges (a,b) must have a<b |
---|
| 46 | insert_edge(edges,regions, (ii,jj),kk,n) |
---|
| 47 | insert_edge(edges,regions, (jj,kk),ii,n) |
---|
| 48 | insert_edge(edges,regions, (ii,kk),jj,n) |
---|
| 49 | return edges, [ regions[i] for i in range(len(regions)) ] |
---|
| 50 | |
---|
| 51 | def dual_edges(edges): |
---|
| 52 | for (ii,jj),((n1,kk1),(n2,kk2)) in edges.items(): |
---|
| 53 | yield (n1,n2) |
---|
| 54 | |
---|
| 55 | def regions_vertices(regions, edges): |
---|
| 56 | for region in regions: |
---|
| 57 | for edge in region: |
---|
| 58 | (n1,p1), (n2,p2) = edges[edge] |
---|
| 59 | yield n1 |
---|
| 60 | yield n2 |
---|