City-Lifeline-Cockpit/tools/city-massing/verify-roof-geometry.py

"""Exercise the actual Blender building generator, without running a city export.
blender --background --factory-startup --python-exit-code 1 --python tools/city-massing/verify-roof-geometry.py
"""
import ast
import math
import sys
from pathlib import Path

import bpy
import bmesh
from mathutils import Vector
from mathutils.geometry import tessellate_polygon

source = Path(__file__).with_name('assemble-cad-bim.py')
module = ast.parse(source.read_text(encoding='utf-8'))
# Load the real function definitions; city intake and export are deliberately not executed.
functions = {'ring_area', 'as_ccw', 'prism_triangles', 'facade_window_count', 'massing_triangles',
             'poly_prism', 'prism_volume_sign'}
namespace = {'bpy': bpy, 'bmesh': bmesh, 'math': math, 'Vector': Vector,
             'tessellate_polygon': tessellate_polygon, 'ARCHITECTURE': '--legacy' not in sys.argv}
exec(compile(ast.Module(body=[node for node in module.body
                             if isinstance(node, ast.FunctionDef) and node.name in functions],
                        type_ignores=[]), str(source), 'exec'), namespace)
materials = [bpy.data.materials.new(name) for name in ('Wall', 'Roof', 'Window', 'Door')]
fixtures = [
    ('rectangle', [(0, 0), (12, 0), (12, 8), (0, 8)]),
    ('concave-L', [(0, 0), (10, 0), (10, 3), (4, 3), (4, 10), (0, 10)]),
    ('clockwise', [(0, 0), (0, 8), (12, 8), (12, 0)]),
    ('thin', [(0, 0), (12, 0), (12, 0.5), (0, 0.5)]),
    ('centimetre-backtrack', [(-881.93, 844.99), (-881.94, 844.99), (-874.09, 842.72),
                             (-876.51, 834.71), (-884.25, 837.04)]),
]
failures = []
for name, ring in fixtures:
    for level in ('lod1', 'lod0'):
        outline = namespace['as_ccw'](ring)
        bottom, top = 140.0, 156.4
        obj = namespace['poly_prism'](outline, bottom, top, materials[0], materials[1],
                                     (6, 4, 150), 2, level, materials[2], materials[3])
        mesh = obj.data
        mesh.calc_loop_triangles()
        top_faces = [face for face in mesh.polygons if all(
            abs(mesh.vertices[i].co.z + obj.location.z - top) < 0.001 for i in face.vertices)]
        bottom_faces = [face for face in mesh.polygons if all(
            abs(mesh.vertices[i].co.z + obj.location.z - bottom) < 0.001 for i in face.vertices)]
        expected = abs(namespace['ring_area'](outline))
        area = sum(face.area for face in top_faces if face.normal.z > 0.99 and face.material_index == 1)
        bottom_area = sum(face.area for face in bottom_faces if face.normal.z < -0.99)
        print(f'ROOF_GEOMETRY {name}/{level}: upwardRoofArea={area:.3f}, expected={expected:.3f}, bottomArea={bottom_area:.3f}', flush=True)
        if abs(area - expected) > 0.001 or abs(bottom_area - expected) > 0.001:
            failures.append(f'{name}/{level}: roof or bottom cap missing/reversed')
        if namespace['ARCHITECTURE']:
            door_faces = [face for face in mesh.polygons if face.material_index == 3 and len(face.vertices) == 4
                          and min(mesh.vertices[i].co.z + obj.location.z for i in face.vertices) < 151]
            door_vertices = {i for face in door_faces for i in face.vertices}
            if not door_vertices or min(mesh.vertices[i].co.z + obj.location.z for i in door_vertices) < 150:
                failures.append(f'{name}/{level}: entrance is buried below seated floor')
            longest = max(range(len(outline)), key=lambda i: math.dist(outline[i], outline[(i + 1) % len(outline)]))
            a, b = outline[longest], outline[(longest + 1) % len(outline)]
            length = math.dist(a, b)
            projected = [(mesh.vertices[i].co.x * (b[0]-a[0]) + mesh.vertices[i].co.y * (b[1]-a[1])) / length
                         for i in door_vertices]
            if projected and max(projected) - min(projected) > 1.6:
                failures.append(f'{name}/{level}: entrance spans an implausible facade fraction')
        # Boundary edges on the solid shell must be closed; facade quads can have their own edges.
        edit = bmesh.new()
        edit.from_mesh(mesh)
        edit.verts.ensure_lookup_table()
        shell_open = sum(not edge.is_manifold for edge in edit.edges
                         if all(vertex.index < 2 * len(outline) for vertex in edge.verts))
        edit.free()
        if shell_open:
            failures.append(f'{name}/{level}: {shell_open} open shell edges')
        bpy.data.objects.remove(obj, do_unlink=True)
        bpy.data.meshes.remove(mesh)
if failures:
    raise AssertionError('; '.join(failures))
print(f'ROOF_GEOMETRY passed: {len(fixtures) * 2} cases', flush=True)