3ea43ae37f6a — Leonard Ritter tip 10 days ago
* more work on bcc lattice
2 files changed, 323 insertions(+), 34 deletions(-)

A => testing/bcc_sh_table.sage
M testing/bcc_shearing.sc
A => testing/bcc_sh_table.sage +194 -0
@@ 0,0 1,194 @@ 
+
+A = 4 * pi
+SH1_C0 = sqrt(3 / (4 * pi))
+SH1_C1 = sqrt(1 / (4 * pi))
+SH1_Basis = (SH1_C0, SH1_C0, SH1_C0, SH1_C1)
+
+# clean solid angle computation for bcc grid
+
+def vxor3 (a, b):
+    return (int(a[0]) ^^ int(b[0]), int(a[1]) ^^ int(b[1]), int(a[2]) ^^ int(b[2]))
+def vand3 (a, b):
+    return (a[0] & b[0], a[1] & b[1], a[2] & b[2])
+def vor3 (a, b):
+    return (a[0] | b[0], a[1] | b[1], a[2] | b[2])
+def vadd3 (a, b):
+    return (a[0] + b[0], a[1] + b[1], a[2] + b[2])
+def vadd3s (a, b):
+    return (a[0] + b, a[1] + b, a[2] + b)
+def vsub3 (a, b):
+    return (a[0] - b[0], a[1] - b[1], a[2] - b[2])
+def vsub3s (a, b):
+    return (a[0] - b, a[1] - b, a[2] - b)
+def vmul3 (a, b):
+    return (a[0] * b[0], a[1] * b[1], a[2] * b[2])
+def vmul4 (a, b):
+    return (a[0] * b[0], a[1] * b[1], a[2] * b[2], a[3] * b[3])
+def vmul4s (a, b):
+    return (a[0] * b, a[1] * b, a[2] * b, a[3] * b)
+def vmul3s (a, b):
+    return (a[0] * b, a[1] * b, a[2] * b)
+
+def dot3 (a, b):
+    return (a[0] * b[0] + a[1] * b[1] + a[2] * b[2])
+def length (v):
+    return sqrt(dot3(v,v))
+def normalize (v):
+    return vmul3s(v, 1 / length(v))
+
+def center3 (a,b,c):
+    return vmul3s(vadd3(vadd3(a,b),c),1/3)
+def center4 (a,b,c,d):
+    return vmul3s(vadd3(vadd3(a,b),vadd3(c,d)),1/4)
+
+def bcc (v):
+    s = (v[0] + v[1] + v[2]) / 6
+    return vsub3s(v, s)
+
+def cc (v):
+    s = (v[0] + v[1] + v[2]) / 3
+    return vadd3s(v, s)
+
+def anglebasis(a, b):
+    ab = dot3(a, b)
+    aabb = dot3(a,a) * dot3(b, b)
+    return (ab / sqrt(aabb), sqrt(1 - ab * ab / aabb))
+
+def solid_angle (u, v, w):
+    return \
+        acos ((u[0] - v[0] * w[0]) / (v[1] * w[1])) + \
+        acos ((v[0] - w[0] * u[0]) / (w[1] * u[1])) + \
+        acos ((w[0] - u[0] * v[0]) / (u[1] * v[1])) - pi
+
+def simplex_solid_angle(U, V, W):
+    u = anglebasis(V, W)
+    v = anglebasis(W, U)
+    w = anglebasis(U, V)
+    return solid_angle(u, v, w)
+
+def edge (code):
+    return (code & 1, (code >> 1) & 1, (code >> 2) & 1)
+
+def simplex (n):
+    a = (0b101100110010011001 >> (n * 3)) & 7
+    b = (0b100110010011001101 >> (n * 3)) & 7
+    return (bcc(edge(a)), bcc(edge(b)), bcc(edge(0b000)), bcc(edge(0b111)))
+
+def simplex_code (n):
+    a = (0b101100110010011001 >> (n * 3)) & 7
+    b = (0b100110010011001101 >> (n * 3)) & 7
+    return (edge(a), edge(b), edge(0b000), edge(0b111))
+
+def face_vertices (s, n):
+    return (
+        s[(0b00011011 >> (n * 2)) & 3],
+        s[(0b01001110 >> (n * 2)) & 3],
+        s[(0b10110001 >> (n * 2)) & 3])
+
+def adj (s, f):
+    o = [0,0,0]
+    if (f == 0):
+        simplex = s + 5
+    elif (f == 1):
+        simplex = s + 1
+    elif (f == 2):
+        o[s >> 1] = 1
+        simplex = s + 4 - ((s & 1) << 1)
+    elif (f == 3): # (f == 3)
+        o[((s + 3) % 6) >> 1] = -1
+        simplex = s + 2 + ((s & 1) << 1)
+    else:
+        assert False
+    return ((simplex % 6), (f ^^ 1), bcc(tuple(o)))
+
+syms = []
+symtable = {}
+symtable[0] = 0
+usedsyms = {}
+rsymtable = {}
+def symlify(x, s = None):
+    global nextsym
+    if (not (x in symtable)):
+        assert s
+        k = 1
+        while (s + str(k)) in usedsyms:
+            k += 1
+        usedsyms[s + str(k)] = True
+        d = var(s + str(k))
+        sx = x.simplify()
+        if sx.substitute(rsymtable) < 0:
+            sx = -sx
+        symtable[-sx] = -d
+        symtable[sx] = d
+        symtable[d] = d
+        symtable[-d] = -d
+        rsymtable[d] = sx
+        syms.append((d,sx))
+    return symtable[x]
+
+table = [[[] for i in range(4)] for i in range(6)]
+
+# gather pattern
+for i in range(6):
+    #print("simplex",i)
+    s = simplex(i)
+    for f in range(4):
+        #print("face",f)
+        ai,af,ofs = adj(i,f)
+        aa = simplex(ai)
+        aa = vadd3(aa[0],ofs),vadd3(aa[1],ofs),vadd3(aa[2],ofs),vadd3(aa[3],ofs)
+        c = center4(*aa)
+        for gg in range(3):
+            g = (f + 1 + gg) % 4
+            vs = face_vertices(s, g)
+            v0 = vsub3(vs[0], c)
+            v1 = vsub3(vs[1], c)
+            v2 = vsub3(vs[2], c)
+            vc = center3(vs[0],vs[1],vs[2])
+            n = vsub3(vc,c)
+            n = (*normalize(n), 1)
+            sa = symlify(simplex_solid_angle(v0,v1,v2).simplify(), 'sa')
+            shn = vmul4(n,SH1_Basis)
+            shn = ([symlify(x,'t') for x in shn])
+            table[i][f].append(vmul4s(shn, sa))
+
+    #print()
+
+for i in range(6):
+    for f in range(4):
+        m = table[i][f]
+        m.sort(key = lambda k: str(k[3]))
+
+for q in range(3):
+    for f in range(4):
+        for i in range(6):
+            W = table[i][f][q]
+            x,y,z,w = W
+            if (q < 2):
+                b = 'u'
+            else:
+                b = 'v'
+            x = symlify(x, b)
+            y = symlify(y, b)
+            z = symlify(z, b)
+            w = symlify(w, 'q')
+
+print("Tx3x4x6")
+for i in range(6):
+    print("    Tx3x4")
+    for f in range(4):
+        s = "        Tx3 "
+        for q in range(3):
+            W = table[i][f][q]
+            x,y,z,w = W
+            x = symlify(x)
+            y = symlify(y)
+            z = symlify(z)
+            w = symlify(w)
+            s += "(T " + repr(x) + " " + repr(y) + " " + repr(z) + " " + repr(w) + ") "
+        print(s)
+print()
+print("symbols:")
+for k,v in syms:
+    print(k,":=",v)
+

          
M testing/bcc_shearing.sc +129 -34
@@ 17,6 17,10 @@ 
     * plane along +YZ splits simplices 1 2 3 | 4 5 0
     * plane along +XZ splits simplices 0 1 2 | 3 4 5
     * plane along +XY splits simplices 5 0 1 | 2 3 4
+    * the solid angles of a simplices neighbors visible faces are, approximately
+        a := 2*acos(3/sqrt(14)) - acos(6/7); for the smaller one
+        and
+        b := (pi - a)/2 = pi/2 - acos(3/sqrt(14)) + acos(6/7)/2
 
 using import glm
 

          
@@ 38,14 42,14 @@ fn simplex-volume (A B C D)
 fn anglebasis (a b)
     ab := (dot a b)
     aabb := (dot a a) * (dot b b)
-    vec2
+    dvec2
         ab / (sqrt aabb)
-        sqrt (1.0 - ab * ab / aabb)
+        sqrt (1.0:f64 - ab * ab / aabb)
 
 # for when both arguments are already normalized
 fn nanglebasis (a b)
     ab := (dot a b)
-    vec2 ab (sqrt (1.0 - ab * ab))
+    dvec2 ab (sqrt (1.0:f64 - ab * ab))
 
 # solid angle of basis angles, specified as sin/cos pairs
 fn simplex-solid-angle (u v w)

          
@@ 54,8 58,8 @@ fn simplex-solid-angle (u v w)
             acos ((u.x - v.x * w.x) / (v.y * w.y))
             acos ((v.x - w.x * u.x) / (w.y * u.y))
             acos ((w.x - u.x * v.x) / (u.y * v.y))
-            -pi
-    ? (A != A) 0.0 A
+            -pi:f64
+    ? (A != A) 0.0:f64 A
 
 # solid angle of basis from normalized vectors
 fn simplex-points-solid-angle (U V W)

          
@@ 76,6 80,10 @@ fn pyramid-points-solid-angle (A B C D)
         simplex-solid-angle bc ac ab
         simplex-solid-angle da ac cd
 
+A:f64 := 4.0:f64 * pi:f64
+SH1_Basis:f64 := (sqrt ((dvec4 (dvec3 3.0:f64) 1.0:f64) / A:f64)) #* (dvec4 1 1 1 1)
+SH1_Basis := (vec4 SH1_Basis:f64)
+
 inline cc->bcc (p)
     p - (p.x + p.y + p.z) / 6.0
 

          
@@ 89,6 97,9 @@ fn cos_angle (a b)
 
 # constructive edges are either unit length (orthogonal) or (sqrt 3) / 2 (body diagonal)
 
+fn iedge (code)
+    (ivec3 code (code >> 1) (code >> 2)) & 1
+
 fn edge (code)
     vec3 ((ivec3 code (code >> 1) (code >> 2)) & 1)
 

          
@@ 178,6 189,59 @@ fn neighbor (s f)
             trap;
     _ (simplex % 6) (f ^ 1) (copy o)
 
+define sh1-constants
+    sa2 := 2 * (acos (3 / (sqrt 14.0))) - (acos (6 / 7)); sa1 := (pi - sa2) / 2
+    t4 := 1 / (2 * (sqrt pi))
+    q1 := sa1 * t4; q2 := sa2 * t4
+    t11 := 1 / (sqrt (87 * pi))
+    t1 := 7 / 2 * t11; t2 := 2 * t11; t3 := 7 * t11
+    t9 := 8 * t11; t10 := 1 / 2 * t11; t12 := 4 * t11
+    t6 := (sqrt (3 / (29 * pi)))
+    t5 := 5 / 2 * t6; t7 := 3 / 2 * t6; t8 := 2 * t6
+    u1 := sa1 * t1; u2 := sa1 * t2; u3 := sa1 * t3
+    u4 := sa1 * t5; u5 := sa1 * t6; u6 := sa1 * t8
+    u7 := sa1 * t7; u8 := sa1 * t10; u9 := sa1 * t9
+    u10 := sa1 * t11
+    v1 := sa2 * t7; v2 := sa2 * t6; v3 := sa2 * t8
+    v4 := sa2 * t3; v5 := sa2 * t10; v6 := sa2 * t12
+    T := vec4; Tx3 := (array T 3); Tx3x4 := (array Tx3 4)
+    (array Tx3x4 6)
+        Tx3x4
+            Tx3 (T u1 u2 -u3 q1) (T u4 0 -u5 q1) (T v1 -v2 -v3 q2)
+            Tx3 (T 0 -u4 u5 q1) (T -u2 -u1 u3 q1) (T v2 -v1 v3 q2)
+            Tx3 (T -u6 -u7 -u5 q1) (T -u9 -u8 u10 q1) (T -v4 v5 -v6 q2)
+            Tx3 (T u8 u9 -u10 q1) (T u7 u6 u5 q1) (T -v5 v4 v6 q2)
+        Tx3x4
+            Tx3 (T -u2 u3 -u1 q1) (T 0 u5 -u4 q1) (T v2 v3 -v1 q2)
+            Tx3 (T u4 -u5 0 q1) (T u1 -u3 u2 q1) (T v1 -v3 -v2 q2)
+            Tx3 (T -u6 -u5 -u7 q1) (T -u9 u10 -u8 q1) (T -v4 -v6 v5 q2)
+            Tx3 (T u8 -u10 u9 q1) (T u7 u5 u6 q1) (T -v5 v6 v4 q2)
+        Tx3x4
+            Tx3 (T -u3 u1 u2 q1) (T -u5 u4 0 q1) (T -v3 v1 -v2 q2)
+            Tx3 (T u5 0 -u4 q1) (T u3 -u2 -u1 q1) (T v3 v2 -v1 q2)
+            Tx3 (T -u5 -u6 -u7 q1) (T u10 -u9 -u8 q1) (T -v6 -v4 v5 q2)
+            Tx3 (T -u10 u8 u9 q1) (T u5 u7 u6 q1) (T v6 -v5 v4 q2)
+        Tx3x4
+            Tx3 (T -u1 -u2 u3 q1) (T -u4 0 u5 q1) (T -v1 v2 v3 q2)
+            Tx3 (T 0 u4 -u5 q1) (T u2 u1 -u3 q1) (T -v2 v1 -v3 q2)
+            Tx3 (T -u7 -u6 -u5 q1) (T -u8 -u9 u10 q1) (T v5 -v4 -v6 q2)
+            Tx3 (T u9 u8 -u10 q1) (T u6 u7 u5 q1) (T v4 -v5 v6 q2)
+        Tx3x4
+            Tx3 (T u2 -u3 u1 q1) (T 0 -u5 u4 q1) (T -v2 -v3 v1 q2)
+            Tx3 (T -u4 u5 0 q1) (T -u1 u3 -u2 q1) (T -v1 v3 v2 q2)
+            Tx3 (T -u7 -u5 -u6 q1) (T -u8 u10 -u9 q1) (T v5 -v6 -v4 q2)
+            Tx3 (T u9 -u10 u8 q1) (T u6 u5 u7 q1) (T v4 v6 -v5 q2)
+        Tx3x4
+            Tx3 (T u3 -u1 -u2 q1) (T u5 -u4 0 q1) (T v3 -v1 v2 q2)
+            Tx3 (T -u5 0 u4 q1) (T -u3 u2 u1 q1) (T -v3 -v2 v1 q2)
+            Tx3 (T -u5 -u7 -u6 q1) (T u10 -u8 -u9 q1) (T -v6 v5 -v4 q2)
+            Tx3 (T -u10 u9 u8 q1) (T u5 u6 u7 q1) (T v6 v4 -v5 q2)
+
+global sh1-constants = sh1-constants
+
+inline gather-solid-angle (s f)
+    unpack (sh1-constants @ s @ f)
+
 inline sort3 (pred a b c)
     let a b =
         if (pred a b) (_ b a)

          
@@ 217,6 281,23 @@ inline simplexfacecoords (s f ofs)
                         print "simplices" i j "faces" f g
                         #print va vb
 
+fn dv4xvT (a)
+    dmat4
+        a.x * a
+        a.y * a
+        a.z * a
+        a.w * a
+
+type+ dmat4x4
+    inline __+ (cls T)
+        static-if (cls == T)
+            inline (a b)
+                dmat4
+                    (a @ 0) + (b @ 0)
+                    (a @ 1) + (b @ 1)
+                    (a @ 2) + (b @ 2)
+                    (a @ 3) + (b @ 3)
+
 for i in (range 6)
     s := (simplex i)
     print "simplex" i

          
@@ 227,45 308,59 @@ for i in (range 6)
                 inline (i) (cc->bcc (edge i))
                 #inline (i) (edge i)
                 unpack s
-    center :=
-        /
-            +
-                v @ 0
-                v @ 1
-                v @ 2
-                v @ 3
-            4
-    print
-        simplex-volume (v @ 0) (v @ 1) (v @ 2) (v @ 3)
+    local sa = 0.0:f64
+    probe := (dvec4 (normalize (dvec3 1 1 1)) 1) * pi:f64 * SH1_Basis:f64
+    #probe := (dvec4 1 0 0 1) * 4 * pi:f64 * SH1_Basis:f64
+    local sh = (dvec4 0)
     for f in (range 4)
-        let v1 v2 v3 =
+        print "face" f
+        #let v1 v2 v3 =
             va-map
                 inline (i) (v @ i)
                 unpack (simplex-local-face-vertices f)
-        c := ((+ v1 v2 v3) / 3)
         let i1 f1 d1 = (neighbor i f)
         local vb =
             arrayof vec3
                 va-map
-                    inline (i) (cc->bcc (edge i))
+                    inline (i) (cc->bcc ((edge i) + (vec3 d1)))
                     #inline (i) (edge i)
                     unpack (simplex i1)
-        let vb1 vb2 vb3 =
-            va-map
-                inline (i) (vb @ i)
-                unpack (simplex-local-face-vertices f1)
-        c1 := (((+ vb1 vb2 vb3) / 3) + (cc->bcc (vec3 d1)))
-        #c1 := ((+ vb1 vb2 vb3) / 3)
-        #A := (triangle-area v1 v2 v3)
-        #print v1 v2 v3 A
-        print (triangle-normal v1 v2 v3)
-        print (c - c1)
-        print (dot (triangle-normal v1 v2 v3) (triangle-normal vb1 vb2 vb3))
-        #print
-            simplex-points-solid-angle
-                normalize (v1 - center)
-                normalize (v2 - center)
-                normalize (v3 - center)
+        center :=
+            /
+                +
+                    vb @ 0
+                    vb @ 1
+                    vb @ 2
+                    vb @ 3
+                4
+        for gg in (range 3)
+            let g = ((f + 1 + gg) % 4)
+            let vb1 vb2 vb3 =
+                va-map
+                    inline (i) (v @ i)
+                    unpack (simplex-local-face-vertices g)
+            ct := (((vb1 + vb2 + vb3) / 3) - center)
+            ct := (normalize ct)
+            gsa :=
+                simplex-points-solid-angle
+                    normalize (dvec3 (vb1 - center))
+                    normalize (dvec3 (vb2 - center))
+                    normalize (dvec3 (vb3 - center))
+            shb := (dvec4 ct 1) * SH1_Basis:f64 * gsa
+            let n0 n1 n2 = (gather-solid-angle i f)
+            let l1 l2 l3 =
+                length ((vec4 shb) - n0)
+                length ((vec4 shb) - n1)
+                length ((vec4 shb) - n2)
+            print (min l1 l2 l3)
+            assert ((min l1 l2 l3) < 0.00001)
+            sh += (dot shb probe) * shb
+            print "neighbor face" i1 g ct gsa
+            sa += gsa
+        #sh += shc * probe
+    print "total" (sa / (4 * pi:f64))
+    print
+        (dot sh sh) / (4 * pi)
 
     print;