#!/usr/bin/env python from __future__ import division from other.core import * import numpy as np from numpy import linalg as LA import sys #import other.cgal as cg octa_base_vec = { 'cube':[1,1,1], 'cuboctahedron':[1,0,1], 'octahedron':[1,0,0], 'small-rhombicuboctahedron': [1,1,sqrt(2)+1], 'truncated-cube':[1,1,sqrt(2)-1], 'truncated-octahedron':[1,0,2], } phi = (1+sqrt(5))/2 icos_base_vec = { 'dodecahedron':[1+phi,1,0], 'icosahedron':[0,1,phi], 'icosidodecahedron':[0,0,1], 'small-rhombicosidodecahedron':[-1,phi,2], 'truncated-dodecahedron':[2,phi,1+phi], 'truncated-icosahedron':[1+2*phi,phi,2] } def close_find(l,x,tol=1e-5): for i,y in enumerate(l): if relative_error(x,y)=0: break else: G.append(ab) grew = 1 return G def real_eigvec(m): val,vec = LA.eig(m) i = close_find(val,1.) return vec[i] def icosahedral_group(): P = Matrix([[0,1,0],[0,0,1],[1,0,0]],dtype=float) p = (+1+sqrt(5))/2 q = (-1+sqrt(5))/2 A = Matrix([[q,-p,1],[p,1,q],[-1,q,p]])/2 assert allclose(dot(A.T,A),eye(3)) G = close_group([P,A]) assert len(G)==60 return G def octahedral_group(): R1 = Matrix([[0,0,1],[1,0,0],[0,1,0]],dtype=float) R2 = Matrix([[0,1,0],[-1,0,0],[0,0,1]]) G = close_group([R1,R2]) assert len(G)==24 return G def orbit(G,v): v = np.array(v) X = v.reshape(1,3) close_points = np.array([]) for g in G: y = g*v if close_find(X,y)<0: X=vstack((X,y)) return X def find_char_length(G,v): o = orbit(G,v) l = [magnitude(x-v) for x in o] mins,ind = close_mins(l[1:],1e-5) return mins[0] def eig(G): return np.linalg.eig(G) def gv(name): if name in icos_base_vec.keys(): g = icosahedral_group() v = icos_base_vec[name] elif name in octa_base_vec.keys(): g = octahedral_group() v = octa_base_vec[name] else: print "polyhedron name not found" return -1 return g,v def find_edges(pts,l,tol=1e-5): edges = [] for i,x in enumerate(pts): mags = [magnitude(x-y) for y in pts] inds = close_finds(mags,l) for j in inds: if j>i: edges.append([i,j]) return edges def m_rep(p): return "{%(x)04f,%(y)04f,%(z)04f}" % {"x":p[0],"y":p[1],"z":p[2]} def e_rep(p1,p2): return "%(x1)05f,%(y1)05f,%(z1)05f,%(x2)05f,%(y2)05f,%(z2)05f" % {"x1":p1[0],"y1":p1[1],"z1":p1[2],"x2":p2[0],"y2":p2[1],"z2":p2[2]} def v_rep(p1): return "%(x1)05f,%(y1)05f,%(z1)05f" % {"x1":p1[0],"y1":p1[1],"z1":p1[2]} if __name__ == "__main__": name = 'small-rhombicosidodecahedron' g,v = gv(name) v = asarray(v)/magnitude(v) l = find_char_length(g,v) pts = orbit(g,v) edges = find_edges(pts,l) # print edges # print pts print "n_verts=",len(pts)," n_edges=",len(edges)," char_length=",l f = open(name+'.verts','w') for pt in pts: f.write(v_rep(pt)+'\n') f = open(name+'.edges','w') for e in edges: f.write(e_rep(pts[e[0]],pts[e[1]])+'\n') # for e in edges: # print "Line[{"+m_rep(pts[e[0]])+','+m_rep(pts[e[1]])+'}],' # for poly,vector in icos_base_vec.items(): # print poly # print find_tri(icosahedral_group(), vector) # for poly,vector in octa_base_vec.items(): # print poly # print find_tri(octahedral_group(), vector)