# Geometry helpers adapted from the read-only Breton atom builder.
import bpy,bmesh,math
from mathutils import Vector
MATS=[]
TILES=[]
class Mesh:
 def __init__(self,name):self.name=name;self.v=[];self.f=[];self.mi=[];self.uv=[];self.smooth=[];self.colors=[];self.color=(1,1,1,1)
 def add(self,verts,faces,mat,smooth=False,uv=None):
  vs=[Vector(v) for v in verts];off=len(self.v);self.v.extend(vs)
  for face in faces:
   self.f.append(tuple(off+j for j in face));self.mi.append(mat);self.smooth.append(smooth);self.colors.append(self.color)
   if uv:self.uv.append([uv[j] for j in face]);continue
   normal=(vs[face[1]]-vs[face[0]]).cross(vs[face[2]]-vs[face[0]]).normalized()
   # Orthonormal surface basis preserves metres even on bevels and pitched roofs.
   if abs(normal.y)<.95:
    up=Vector((0,1,0));u=up.cross(normal).normalized();v=normal.cross(u).normalized()
   else:
    u=(Vector((1,0,0))-normal*normal.x).normalized();v=normal.cross(u)
    if v.z<0:v=-v
   self.uv.append([(p.dot(u)/TILES[mat],p.dot(v)/TILES[mat]) for p in [vs[j] for j in face]])
 def finish(self):
  me=bpy.data.meshes.new(self.name);me.from_pydata(self.v,[],self.f);me.update();ob=bpy.data.objects.new(self.name,me);bpy.context.collection.objects.link(ob)
  for m in MATS:me.materials.append(m)
  uv=me.uv_layers.new(name='UVMap')
  for p,mi,coords,sm in zip(me.polygons,self.mi,self.uv,self.smooth):
   p.material_index=mi;p.use_smooth=sm
   for li,co in zip(p.loop_indices,coords):uv.data[li].uv=co
  color=me.color_attributes.new(name='Color',type='FLOAT_COLOR',domain='CORNER')
  for p,c in zip(me.polygons,self.colors):
   for li in p.loop_indices:color.data[li].color=c
  bm=bmesh.new();bm.from_mesh(me);bmesh.ops.recalc_face_normals(bm,faces=bm.faces);bmesh.ops.triangulate(bm,faces=list(bm.faces));bm.to_mesh(me);bm.free();me.update()
  # Flat faces with bevel normals preserve the dressed silhouette without a modifier.
  return ob

def box(g,c,s,m,b=.01,rot=None):
 bm=bmesh.new();bmesh.ops.create_cube(bm,size=1)
 for v in bm.verts:v.co=Vector((v.co.x*s[0],v.co.y*s[1],v.co.z*s[2]))
 if b:bmesh.ops.bevel(bm,geom=list(bm.edges),offset=min(b,min(s)*.24),segments=1,affect='EDGES')
 bm.verts.ensure_lookup_table();bm.verts.index_update();vs=[v.co.copy() for v in bm.verts];fs=[tuple(v.index for v in f.verts) for f in bm.faces];bm.free()
 if rot:vs=[rot@v for v in vs]
 g.add([v+Vector(c) for v in vs],fs,m)

def prism(g,poly,depth,m,axis='z'):
 # Polygon in XY, extrusion along Z; optional rigid axis mapping by caller.
 vs=[(x,y,z) for z in depth for x,y in poly];n=len(poly)
 fs=[tuple(range(n-1,-1,-1)),tuple(range(n,2*n))]+[(i,(i+1)%n,(i+1)%n+n,i+n) for i in range(n)]
 if axis=='y':vs=[(x,z,y) for x,y,z in vs]
 g.add(vs,fs,m)

def cyl(g,a,b,r,m,n=8,r2=None):
 a=Vector(a);b=Vector(b);L=(b-a).length;q=(b-a).to_track_quat('Z','Y');r2=r if r2 is None else r2
 vs=[p+q@Vector((rr*math.cos(i*math.tau/n),rr*math.sin(i*math.tau/n),0)) for p,rr in [(a,r),(b,r2)] for i in range(n)]
 fs=[tuple(range(n-1,-1,-1)),tuple(range(n,2*n))]+[(i,(i+1)%n,(i+1)%n+n,i+n) for i in range(n)]
 g.add(vs,fs,m,True)

def tube(g,points,r,m,n=8):
 for a,b in zip(points,points[1:]):cyl(g,a,b,r,m,n)