Try subdividing boundaries to break up facets near torsion

examples.py

@@ -136,6 +136,8 @@ def gen_twisted_boundary(count=4, dx0=2, dz=0.2, ang=0.1):
         [1, 0, 1],
         [0, 0, 1],
     ], dtype=numpy.float64) - [0.5, 0, 0.5]
+    b = meshutil.subdivide_boundary(b)
+    b = meshutil.subdivide_boundary(b)
     # Generate 'seed' transformations:
     xfs = [meshutil.Transform().translate(dx0, 0, 0).rotate([0,1,0], numpy.pi * 2 * i / count)
            for i in range(count)]

meshutil.py

@@ -1,6 +1,7 @@
 import stl.mesh
 import numpy
 import quaternion
+import random
 
 import quat
 
@@ -158,7 +159,7 @@ def cube_distort(angle, open_xz=False):
         # winding order?
     return FaceVertexMesh(verts, faces)
 
-def subdivide_boundary(bound):
+def split_boundary(bound):
     # assume bound1 has shape (4,3).
     # Midpoints of every segment:
     mids = (bound + numpy.roll(bound, 1, axis=0)) / 2
@@ -173,6 +174,16 @@ def subdivide_boundary(bound):
     ]
     return bounds
 
+def subdivide_boundary(bound):
+    # assume bound1 has shape (4,3).
+    # Midpoints of every segment:
+    mids = (bound + numpy.roll(bound, -1, axis=0)) / 2
+    b2 = numpy.zeros((bound.shape[0]*2, bound.shape[1]))
+    for i,row in enumerate(bound):
+        b2[2*i,:] = bound[i,:]
+        b2[2*i+1,:] = mids[i,:]
+    return b2
+
 def join_boundary_simple(bound1, bound2):
     # bound1 & bound2 are both arrays of shape (N,3), representing
     # the points of a boundary.  This joins the two boundaries by
@@ -188,8 +199,12 @@ def join_boundary_simple(bound1, bound2):
     for i in range(n):
         v0 = i
         v1 = (i + 1) % n
-        fs[2*i]     = [n + v1, n + v0, v0]
-        fs[2*i + 1] = [v1,     n + v1, v0]
+        if random.random() < 0.5:
+            fs[2*i]     = [n + v1, n + v0, v0]
+            fs[2*i + 1] = [v1,     n + v1, v0]
+        else:
+            fs[2*i]     = [n + v1, n + v0, v1]
+            fs[2*i + 1] = [v1,     n + v0, v0]
     return FaceVertexMesh(vs, fs)
 
 def close_boundary_simple(bound):