swh:1:snp:317d7c0ae76920777e7a0fe55b5c4d239132d019
Tip revision: 26e88518c0695c213d0013bf42c1ae4fce6ed0a9 authored by lhf on 20 March 2024, 11:25:09 UTC
Create LICENSE
Create LICENSE
Tip revision: 26e8851
kf.py
# A vertex-centric representation for kite fractals
# based on dk.py and the paper
# Fractal tilings based on kite- and dart-shaped prototiles by Fathauer (2001)
# https://doi.org/10.1016/S0097-8493(00)00134-5
from __future__ import print_function
from sys import argv
# depth of refinement
R= int(argv[1]) if len(argv)>1 else 4
# make array
def array(n):
return list(range(n))
# normalize 3-adic lattice coordinates [a,b,m]=(a+b*w^2)/(3^m)
def normalize(a,b,m):
a,b=int(a),int(b)
while m>0 and a%3==0 and b%3==0:
a=a/3
b=b/3
m=m-1
return [a,b,m]
# add 3-adic lattice points with scaling
def add(v1,v2,n):
a1,b1,m1=v1
a2,b2,m2=v2
m2=m2+n//2
d1=3**m1
d2=3**m2
a=d2*a1+d1*a2
b=d2*b1+d1*b2
m=m1+m2
return normalize(a,b,m)
# multiply 3-adic lattice points
# PolynomialRemainder[(a_1 + b_1 z)(a_2 + b_2 z),z^2-z+1,z]
# z (a_2 b_1 + a_1 b_2 + b_2 b_1) + a_1 a_2 - b_1 b_2
def multiply(v1,v2):
a1,b1,m1=v1
a2,b2,m2=v2
a=a1*a2-b1*b2
b=a2*b1+a1*b2+b2*b1
m=m1+m2
return normalize(a,b,m)
# rotate star
def rotate(a,z):
n=len(a)
b=array(n)
for i in range(n):
b[i]=multiply(a[i],z)
return b
# 3-adic lattice coordinates of basic directions; odd directions are shorter
W=array(12)
W[0]=[1,0,0] # 1
W[1]=[1,1,1] # (1+w^2)/3
W[2]=[0,1,0] # w^2
for i in range(3,12):
W[i]=multiply(W[i-2],W[2])
# stars
TYPE=[20,31,32,41,42,43,50,60]
STAR={}
for s in TYPE:
d=s//10
STAR[s]=array(12)
STAR[s][0]=array(d)
# 20
STAR[20][0]=[W[4],W[8]]
# 31
STAR[31][0]=[W[3],W[6],W[9]]
# 41
STAR[41][0]=[W[3],W[6],W[9],W[0]]
STAR[41][0][3]=[1,0,1]
# 32
STAR[32][0]=[W[2],W[6],W[10]]
# 42
STAR[42][0]=[W[2],W[6],W[10],W[0]]
# 43
STAR[43][0]=[W[3],W[6],W[9],W[0]]
# 50
STAR[50][0]=[W[4],W[8],W[10],W[0],W[2]]
# 60
STAR[60][0]=[W[0],W[2],W[4],W[6],W[8],W[10]]
# rotate basic stars
for s in TYPE:
STAR[s][1]=rotate(STAR[s][0],W[1])
for k in range(2,12):
STAR[s][k]=rotate(STAR[s][k-2],W[2])
# opposites
OPP={}
for s in TYPE:
d=s//10
if d<4: d=d-1
OPP[s]=array(12)
OPP[s][0]=array(d)
# 20
OPP[20][0][0]=[-2,0,0]
# 31
OPP[31][0][0]=[-1,1,0]
OPP[31][0][1]=[0,-1,0]
# 41
OPP[41][0][0]=[-1,1,0]
OPP[41][0][1]=[0,-1,0]
OPP[41][0][2]=[2,-1,1]
OPP[41][0][3]=[1,1,1]
# 32
OPP[32][0][0]=[-2,2,0]
OPP[32][0][1]=[0,-2,0]
# 42
OPP[42][0][0]=[-2,2,0]
OPP[42][0][1]=[0,-2,0]
OPP[42][0][2]=[4,-2,1]
OPP[42][0][3]=[2,2,1]
# 43
OPP[43][0][0]=[-1,1,0]
OPP[43][0][1]=[0,-1,0]
OPP[43][0][2]=[1,-1,0]
OPP[43][0][3]=[0,1,0]
# 50
OPP[50][0][0]=[-2,0,0]
OPP[50][0][3]=[2,2,1]
OPP[50][0][4]=multiply(OPP[50][0][3],W[2])
OPP[50][0][2]=multiply(OPP[50][0][3],W[10])
OPP[50][0][1]=multiply(OPP[50][0][3],W[8])
# 60
OPP[60][0][0]=[2,2,1]
for i in range(1,6):
OPP[60][0][i]=multiply(OPP[60][0][i-1],W[2])
# rotate basic opposites
for s in TYPE:
OPP[s][1]=rotate(OPP[s][0],W[1])
for k in range(2,12):
OPP[s][k]=rotate(OPP[s][k-2],W[2])
# Cartesian coordinates
w2=complex(0.5,0.866025403784438646763723170752936183471402626905190314027)
# vertex cloud
V={}
# topological data
U={}
E={}
F={}
nU=0
# use dot notation for dict -- https://stackoverflow.com/a/74214556/107090
class DOTTED(dict):
__getattr__ = dict.get
__setattr__ = dict.__setitem__
__delattr__ = dict.__delitem__
# add vertex to cloud
def addvertex(a,b,m,s,k,n):
global nU
a,b,m=normalize(a,b,m)
if not (a,b,m) in V:
z=(a+b*w2)/(3**m)
t=(a,b,m)
d=s//10
U[nU]=V[a,b,m]=DOTTED({'t':t,'d':d,'k':k,'n':n,'z':z,'id':nU,'s':s})
nU=nU+1
return V[a,b,m]
# base mesh
def basemesh():
v=addvertex(0,0,0,60,0,0)
s=thestar(v)
for i in range(6):
w=s[i]
w.s=31
w.d=3
w.k=2*i
w.n=0
s=theopposites(v)
for i in range(6):
w=s[i]
w.s=20
w.d=2
w.k=2*i+1
w.n=1
def loadmesh(filename):
N=0
for line in open(filename):
if N>0:
a,b,m,d,k,n=list(map(int,line.split(",")))
addvertex(a,b,m,d,k,n)
N=N+1
# stars
def adda(v,i,A):
w=A[v.s][v.k][i]
return add(v.t,w,v.n)
def adj(v,i):
return adda(v,i,STAR)
def adjacent(v,i):
a,b,m=adj(v,i)
return V[a,b,m]
def star(v):
s=array(v.d)
for i in range(v.d):
s[i]=adjacent(v,i)
return s
def thestar(v):
s=array(v.d)
for i in range(v.d):
a,b,m=adj(v,i)
s[i]=addvertex(a,b,m,0,0,0)
return s
# opposites
def opp(v,i):
return adda(v,i,OPP)
def opposite(v,i):
a,b,m=opp(v,i)
return V[a,b,m]
def opposites(v):
s=array(v.d)
for i in range(v.d):
s[i]=opposite(v,i)
return s
def theopposites(v):
s=array(v.d)
for i in range(v.d):
a,b,m=opp(v,i)
s[i]=addvertex(a,b,m,0,0,0)
return s
# faces and edges
def addedge(v1,v2):
a=[v1,v2]
j=a.index(min(a))
a=(a[j],a[(j+1)%2])
E[a]=True
def addface(v1,v2,v3,v4):
a=[v1.id,v2.id,v3.id,v4.id]
j=a.index(min(a))
a=(a[j],a[(j+1)%4],a[(j+2)%4],a[(j+3)%4])
F[a]=True
addedge(a[0],a[1])
addedge(a[1],a[2])
addedge(a[2],a[3])
addedge(a[3],a[0])
def addfaces(v):
s=star(v)
o=opposites(v)
for i in range(v.d):
addface(v,s[i],o[i],s[(i+1)%v.d])
def facetype(k):
v1,v2,v3,v4=k
n=max(U[v1].n,U[v2].n,U[v3].n,U[v4].n)
return str(n)+" qt"
# refinement
def refine(v):
if v.s==0:
return
elif v.s==20:
v.s=50
v.d=5
elif v.s==31:
v.s=41
v.d=4
elif v.s==32:
v.s=42
v.d=4
thestar(v)
theopposites(v)
if v.s==41:
w=adjacent(v,3)
w.s=32
w.d=3
w.k=v.k
w.n=v.n+2
elif v.s==42:
w=adjacent(v,3)
w.s=31
w.d=3
w.k=v.k
w.n=v.n
elif v.s==50:
for j in [2,3,4]:
w=adjacent(v,j)
old= w.s!=0
w.s=31
w.d=3
w.k=(v.k+2*j-6)%12
w.n=v.n
if old:
w.s=43
w.d=4
for j in [2,3]:
w=opposite(v,j)
k1=w.k
old= w.s!=0
w.s=20
w.d=2
w.k=(v.k+2*j-5)%12
w.n=v.n+1
if old:
k2=w.k
if (k1+4)%12==k2:
w.k=(k1+2)%12
else:
w.k=(k2+2)%12
w.s=32
w.d=3
# main -----------------------------------------------------------------------
# initial mesh
basemesh()
# refine mesh
for r in range(R):
#print "%=U",r,R,len(U)
for k in range(len(U)):
v=U[k]
if v.s<40:
refine(v)
#print "%=U",r,R,len(U)
# reconstruct mesh
for k in V:
v=V[k]
if v.s>40:
addfaces(v)
# output PostScript
print('''\
%!PS-Adobe-2.0 EPSF-2.0
%%BoundingBox: 0 0 500 500
%%BoundingBox: 400 200 460 300
250 dup translate
90 dup scale
0 setlinewidth
1 setlinejoin
1 setlinecap
/c { 0.02 0 360 arc fill } bind def
/p { 0.005 0 360 arc fill } bind def
/p { 0.04 0 360 arc fill } bind def
/p { 0.01 0 360 arc fill } bind def
/l { moveto lineto stroke } bind def
/q { moveto lineto lineto lineto closepath fill } bind def
/p0 { 0.5 0.5 0.5 setrgbcolor p } bind def
/p20 { 0.5 0.5 1 setrgbcolor p } bind def
/p31 { 0 1 0 setrgbcolor p } bind def
/p32 { 0 1 1 setrgbcolor p } bind def
/p41 { 1 0 0 setrgbcolor p } bind def
/p42 { 1 0 1 setrgbcolor p } bind def
/p43 { 1.0 0.65 0 setrgbcolor p } bind def
/p50 { 1 1 0 setrgbcolor p } bind def
/p60 { 1 1 1 setrgbcolor p } bind def
/qt { 0.2 mul 0.1 add 1 sub neg dup 1 setrgbcolor q } bind def
/qt { 0.8 6 div mul 0.1 add 1 sub neg dup 1 setrgbcolor q } bind def
''')
print("%=ARG","R",R)
print("%=VEF",R,len(V),len(U),len(E),len(F))
print("")
print("% faces")
print("1 0.8 0.8 setrgbcolor")
for k in F:
z1=U[k[0]].z
z2=U[k[1]].z
z3=U[k[2]].z
z4=U[k[3]].z
print(z1.real,z1.imag, z2.real,z2.imag, z3.real,z3.imag, z4.real,z4.imag,facetype(k))
print("")
print("% edges")
print("0 0 0 setrgbcolor")
print("0.01 setlinewidth")
for k in E:
z1=U[k[0]].z
z2=U[k[1]].z
print(z1.real,z1.imag, z2.real,z2.imag, "l")
"""
print ""
print "% vertices"
print "0 0 0 setrgbcolor"
for k in U:
v=U[k]
if v.s<90:
print v.z.real, v.z.imag, "p" + str(v.s)
"""
print("")
print("showpage")
print("%%EOF")
exit() # if extra files not needed
# output CSV
print("")
PREFIX="%=CSV"
print(PREFIX, "a,b,m,s,k,n")
for k in range(len(U)):
v=U[k]
print(PREFIX, ','.join(map(str,[v.t[0],v.t[1],v.t[2],v.s,v.k,v.n])))
# done
print("")
print("%=EOF")