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 |
---|