From a78026288156c4f6766e7ff7848a3d3ae9ee2eff Mon Sep 17 00:00:00 2001 From: Griffin Long <151978260+DrunkOnJava@users.noreply.github.com> Date: Tue, 3 Nov 2026 11:44:39 -0500 Subject: [PATCH 1/3] ci(measure): mesh_bodies_vs_ifc and saved_meshes_vs_ifc mesh across the CPUs and measure surface distance by exact pruning, output unchanged (#1237) mesh_bodies_vs_ifc opens each file once, meshes in a process pool with one create_shape per part instead of one after another, and prunes surface_distance with block, box and plane lower bounds, evaluating the survivors with the arithmetic the exhaustive search uses so the result is the same float. saved_meshes_vs_ifc meshes through the same pool in place of the geometry iterator, which took twice the CPU per element, and lists the bodies in the iterator's one-thread (ascending instance id) order, which the tables' ties follow. --- tools/re/ifc_common.py | 47 ++++++++ tools/re/mesh_bodies_vs_ifc.py | 183 +++++++++++++++++++++++++------- tools/re/saved_meshes_vs_ifc.py | 23 ++-- 3 files changed, 203 insertions(+), 50 deletions(-) diff --git a/tools/re/ifc_common.py b/tools/re/ifc_common.py index 62c686e9..cc55ec32 100644 --- a/tools/re/ifc_common.py +++ b/tools/re/ifc_common.py @@ -1,7 +1,13 @@ """Helpers the tools/re scorers share, one copy of each so they cannot drift.""" +import concurrent.futures +import concurrent.futures.process +import multiprocessing +import os import re +import sys +import ifcopenshell.geom import ifcopenshell.util.placement import numpy as np @@ -36,3 +42,44 @@ def report_mesh_error(error): if message not in MESH_ERRORS: MESH_ERRORS.add(message) print(f"IfcOpenShell cannot mesh: {message}") + + +_MESHING = {} + + +def _mesh_one(part_id): + try: + shape = ifcopenshell.geom.create_shape(_MESHING["settings"], _MESHING["model"].by_id(part_id)) + except RuntimeError: + return None + return ( + np.array(shape.geometry.verts, dtype=float).reshape(-1, 3), + np.array(shape.geometry.faces, dtype=int).reshape(-1, 3), + ) + + +def mesh_parts(model, settings, parts): + """Instance id -> (vertices (N, 3), faces (M, 3)) for each of `parts` that + IfcOpenShell meshes, one `create_shape` per part spread over the CPUs. The + triangles are those `create_shape` gives, in the same order (the + iterator's are the same, but it took twice the CPU time per element in + a measurement on the BIM4LCA office, and no less wall time with 12 threads + than with 4); a part it cannot mesh is left out.""" + ids = list(dict.fromkeys(part.id() for part in parts)) + cpus = len(os.sched_getaffinity(0)) if hasattr(os, "sched_getaffinity") else os.cpu_count() or 1 + workers = min(cpus, len(ids)) + _MESHING.update(model=model, settings=settings) + meshes = None + if workers > 1: + # A forked worker flushes the buffered output it inherits when it ends. + sys.stdout.flush() + sys.stderr.flush() + try: + with concurrent.futures.ProcessPoolExecutor(workers, mp_context=multiprocessing.get_context("fork")) as pool: + meshes = list(pool.map(_mesh_one, ids, chunksize=1)) + except concurrent.futures.process.BrokenProcessPool: + meshes = None + if meshes is None: + meshes = [_mesh_one(i) for i in ids] + _MESHING.clear() + return {i: mesh for i, mesh in zip(ids, meshes) if mesh is not None} diff --git a/tools/re/mesh_bodies_vs_ifc.py b/tools/re/mesh_bodies_vs_ifc.py index eeaffca4..96dfdda6 100644 --- a/tools/re/mesh_bodies_vs_ifc.py +++ b/tools/re/mesh_bodies_vs_ifc.py @@ -36,12 +36,17 @@ import numpy as np sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) -from opening_boxes_vs_ifc import site_frame # noqa: E402 +from ifc_common import mesh_parts, site_frame # noqa: E402 from roof_bodies_vs_ifc import shapes_of # noqa: E402 from saved_meshes_vs_ifc import body_sources # noqa: E402 SOURCE = "saved_graphics_mesh" PAIR_BUDGET = 1_000_000 +# surface_distance measures every pair up to this many, and beyond it prunes. +EXHAUSTIVE = 200_000 +BLOCK = 8 +SEEDS = 4 +SLIVER = 1e-2 # Most elements scored per class (an even stride through them by Tag), and # most vertices of a body tested against the other's surface, so a file # with thousands of mullions is scored in the time Measure allows. @@ -49,13 +54,14 @@ VERTEX_CAP = 400 -def meshes(path, tags=None, source=None): +def meshes(path, tags=None, source=None, model=None): """Tag -> (class, ObjectType, triangles (T, 3, 3) in internal feet).""" - model = ifcopenshell.open(path) + if model is None: + model = ifcopenshell.open(path) frame = site_frame(model) settings = ifcopenshell.geom.settings() settings.set("use-world-coords", True) - out = {} + chosen = [] for element in model.by_type("IfcElement"): if element.is_a("IfcOpeningElement") or not element.Tag: continue @@ -63,19 +69,20 @@ def meshes(path, tags=None, source=None): continue if source is not None and body_sources(element) != source: continue - parts = [] - for part in shapes_of(element): - try: - shape = ifcopenshell.geom.create_shape(settings, part) - except RuntimeError: + chosen.append((element, shapes_of(element))) + shapes = mesh_parts(model, settings, [part for _, parts in chosen for part in parts]) + out = {} + for element, parts in chosen: + triangles = [] + for part in parts: + if part.id() not in shapes: continue - verts = np.array(shape.geometry.verts, dtype=float).reshape(-1, 3) - faces = np.array(shape.geometry.faces, dtype=int).reshape(-1, 3) + verts, faces = shapes[part.id()] if len(verts) and len(faces): - parts.append(frame(verts)[faces]) - if not parts: + triangles.append(frame(verts)[faces]) + if not triangles: continue - triangles = np.concatenate(parts) + triangles = np.concatenate(triangles) if element.Tag in out: triangles = np.concatenate([out[element.Tag][2], triangles]) out[element.Tag] = (element.is_a(), element.ObjectType or "", triangles) @@ -119,42 +126,140 @@ def find(k): return sorted(shells.values(), key=lambda v: -abs(v)) +def dot(x, y): + return np.einsum("...k,...k->...", x, y) + + def segment_distances(p, a, b): ab = b - a - t = np.einsum("ntk,tk->nt", p - a, ab) / np.maximum(np.einsum("tk,tk->t", ab, ab), 1e-24) + t = dot(p - a, ab) / np.maximum(dot(ab, ab), 1e-24) q = a + np.clip(t, 0.0, 1.0)[..., None] * ab return np.linalg.norm(p - q, axis=-1) -def surface_distance(points, triangles): - """The largest distance from any of `points` to the nearest triangle.""" +def triangle_terms(triangles): a, b, c = triangles[:, 0], triangles[:, 1], triangles[:, 2] ab, ac = b - a, c - a normal = np.cross(ab, ac) length = np.linalg.norm(normal, axis=1) unit = normal / np.maximum(length, 1e-24)[:, None] - d00 = np.einsum("tk,tk->t", ab, ab) - d01 = np.einsum("tk,tk->t", ab, ac) - d11 = np.einsum("tk,tk->t", ac, ac) + d00, d01, d11 = dot(ab, ab), dot(ab, ac), dot(ac, ac) denom = np.where(length > 1e-12, d00 * d11 - d01 * d01, 1.0) - worst = 0.0 - step = max(1, PAIR_BUDGET // len(triangles)) - for start in range(0, len(points), step): - p = points[start : start + step][:, None, :] - ap = p - a - height = np.einsum("ntk,tk->nt", ap, unit) - q = ap - height[..., None] * unit - d20 = np.einsum("ntk,tk->nt", q, ab) - d21 = np.einsum("ntk,tk->nt", q, ac) - v = (d11 * d20 - d01 * d21) / denom - w = (d00 * d21 - d01 * d20) / denom - inside = (v >= -1e-9) & (w >= -1e-9) & (v + w <= 1 + 1e-9) & (length > 1e-12) - plane = np.where(inside, np.abs(height), np.inf) - nearest = np.minimum.reduce( - [plane, segment_distances(p, a, b), segment_distances(p, b, c), segment_distances(p, c, a)] - ).min(axis=1) - worst = max(worst, float(nearest.max())) - return worst + return a, b, c, unit, length, d00, d01, d11, denom + + +def distances(p, terms): + """Distance from points `p` (..., 1, 3) to triangles whose terms broadcast with them.""" + a, b, c, unit, length, d00, d01, d11, denom = terms + ab, ac = b - a, c - a + ap = p - a + height = dot(ap, unit) + q = ap - height[..., None] * unit + d20 = dot(q, ab) + d21 = dot(q, ac) + v = (d11 * d20 - d01 * d21) / denom + w = (d00 * d21 - d01 * d20) / denom + inside = (v >= -1e-9) & (w >= -1e-9) & (v + w <= 1 + 1e-9) & (length > 1e-12) + plane = np.where(inside, np.abs(height), np.inf) + return np.minimum.reduce([plane, segment_distances(p, a, b), segment_distances(p, b, c), segment_distances(p, c, a)]) + + +def morton(centres): + """The order of points along a Z-order curve, so triangles near each other are near in the order.""" + lo, hi = centres.min(axis=0), centres.max(axis=0) + cells = ((centres - lo) / np.maximum(hi - lo, 1e-12) * 1023).astype(np.uint64) + code = np.zeros(len(centres), dtype=np.uint64) + for bit in range(10): + for axis in range(3): + code |= ((cells[:, axis] >> np.uint64(bit)) & np.uint64(1)) << np.uint64(3 * bit + axis) + return np.argsort(code, kind="stable") + + +def surface_distance(points, triangles): + """The largest distance from any of `points` to the nearest triangle. + + Exact, not approximate: for each point the nearest distance is the minimum of + `distances` over a set of triangles that holds every one that could be nearest, + and each of them is measured by the very arithmetic that measures all of them in + a small job (`EXHAUSTIVE` pairs or fewer, and any input that is not finite), so + the result is the same float. A triangle is left out for a point only when a lower + bound on its distance exceeds the nearest distance already found for that point + plus a slack (1e-6 ft and 1e-9 of the coordinates' extent) that is far above the + rounding of these distances and the 1e-9 barycentric tolerance of `inside`. The + bounds hold for every triangle but the slivers (the sine of the angle between two + edges under SLIVER, whose computed plane is not reliable) and the zero-area ones (which + have no plane): the bounds are the distance to a block of eight neighbouring + triangles' box, then to a triangle's box, and for a triangle with a plane the + distance to its plane and to the side planes of its edges. Slivers are measured + against every point. + """ + terms = triangle_terms(triangles) + n = len(points) + if n * len(triangles) <= EXHAUSTIVE or not (np.isfinite(points).all() and np.isfinite(triangles).all()): + worst = 0.0 + step = max(1, PAIR_BUDGET // len(triangles)) + for start in range(0, n, step): + nearest = distances(points[start : start + step][:, None, :], terms).min(axis=1) + worst = max(worst, float(nearest.max())) + return worst + a, b, c, unit, length, d00, d01, d11, denom = terms + sliver = (length > 1e-12) & (length < SLIVER * np.sqrt(d00 * d11)) + flat = (length > 1e-12) & ~sliver + lo = np.minimum(np.minimum(a, b), c) + hi = np.maximum(np.maximum(a, b), c) + extent = max(float(np.abs(triangles).max()), float(np.abs(points).max()), float(np.sqrt(max(d00.max(), d11.max())))) + slack = 1e-6 + 1e-9 * extent + sides = [] + for start, end in ((a, b), (b, c), (c, a)): + out = np.cross(end - start, unit) + sides.append((start, out / np.maximum(np.linalg.norm(out, axis=1), 1e-24)[:, None])) + order = morton((lo + hi) / 2) + order = np.concatenate([order, np.full(-len(order) % BLOCK, order[-1])]) + blocks = order.reshape(-1, BLOCK) + block_lo, block_hi = lo[blocks].min(axis=1), hi[blocks].max(axis=1) + seeds = min(SEEDS, len(blocks)) + nearest = np.empty(n) + step = max(1, PAIR_BUDGET // len(blocks)) + for start in range(0, n, step): + p = points[start : start + step] + gap = np.maximum(np.maximum(block_lo[None] - p[:, None], p[:, None] - block_hi[None]), 0.0) + bound = np.sqrt(dot(gap, gap)) + rows = np.arange(len(p)) + seed = np.argpartition(bound, seeds - 1, axis=1)[:, :seeds] if seeds < len(blocks) else np.tile(np.arange(seeds), (len(p), 1)) + best = np.full(len(p), np.inf) + for column in range(seeds): + best = np.minimum(best, distances(p[:, None, :], tuple(x[blocks[seed[:, column]]] for x in terms)).min(axis=1)) + bound[rows, seed[:, column]] = np.inf + pair_rows, pair_blocks = np.nonzero(bound <= (best + slack)[:, None]) + done = 0 + while done < len(pair_rows): + stop = int(np.searchsorted(pair_rows, pair_rows[min(len(pair_rows) - 1, done + PAIR_BUDGET // BLOCK)], side="right")) + row = np.repeat(pair_rows[done:stop], BLOCK) + tri = blocks[pair_blocks[done:stop]].reshape(-1) + done = stop + q = p[row] + gap = np.maximum(np.maximum(lo[tri] - q, q - hi[tri]), 0.0) + floor = np.sqrt(dot(gap, gap)) + height = dot(q - a[tri], unit[tri]) + outside = np.zeros(len(tri)) + for corner, side in sides: + outside = np.maximum(outside, dot(q - corner[tri], side[tri])) + floor = np.where(flat[tri], np.maximum(floor, np.sqrt(height * height + outside * outside)), floor) + first = np.flatnonzero(np.r_[True, row[1:] != row[:-1]]) + lowest = np.minimum.reduceat(floor, first) + group = np.repeat(np.arange(len(first)), np.diff(np.r_[first, len(row)])) + pick = np.flatnonzero(floor == lowest[group]) + pick = pick[np.unique(group[pick], return_index=True)[1]] + best[row[pick]] = np.minimum(best[row[pick]], distances(q[pick][:, None, :], tuple(x[tri[pick]][:, None] for x in terms))[:, 0]) + keep = floor <= best[row] + slack + row, tri = row[keep], tri[keep] + np.minimum.at(best, row, distances(p[row][:, None, :], tuple(x[tri][:, None] for x in terms))[:, 0]) + nearest[start : start + step] = best + slivers = np.flatnonzero(sliver) + for start in range(0, len(slivers), BLOCK): + part = slivers[start : start + BLOCK] + nearest = np.minimum(nearest, distances(points[:, None, :], tuple(x[part] for x in terms)).min(axis=1)) + return float(nearest.max()) def sample(triangles): @@ -179,8 +284,8 @@ def main(): chosen.update(tags[::step]) if step > 1: print(f"{cls}: {len(tags)} bodies, every {step}th scored") + ours = meshes(sys.argv[1], tags=chosen, source=SOURCE, model=model) del model - ours = meshes(sys.argv[1], tags=chosen, source=SOURCE) ref = meshes(sys.argv[2], tags=set(ours)) print(f"{SOURCE} bodies scored: {len(ours)}; in Revit's export {len(ref)}") rows = collections.defaultdict(lambda: collections.defaultdict(collections.Counter)) diff --git a/tools/re/saved_meshes_vs_ifc.py b/tools/re/saved_meshes_vs_ifc.py index b20193b1..e98db67e 100644 --- a/tools/re/saved_meshes_vs_ifc.py +++ b/tools/re/saved_meshes_vs_ifc.py @@ -35,14 +35,16 @@ import argparse import collections import math -import multiprocessing +import os import statistics +import sys import ifcopenshell import ifcopenshell.geom import numpy as np -from ifc_common import site_frame +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +from ifc_common import mesh_parts, site_frame # noqa: E402 FT = 0.3048 @@ -80,12 +82,14 @@ def world_boxes(path, tags=None, transform=None, wanted=None, frame=True): if not elements: return out to_internal = site_frame(model) if frame else None - it = ifcopenshell.geom.iterator(settings, model, multiprocessing.cpu_count(), include=elements) - if not it.initialize(): - return out - while True: - shape = it.get() - verts = np.array(shape.geometry.verts, dtype=float).reshape(-1, 3) + # The order the geometry iterator returned them in, one thread, which the + # tables' ties (equal counts) follow. + elements.sort(key=lambda e: e.id()) + shapes = mesh_parts(model, settings, elements) + for element in elements: + if element.id() not in shapes: + continue + verts = shapes[element.id()][0] if len(verts): if to_internal is not None: points = to_internal(verts) @@ -94,10 +98,7 @@ def world_boxes(path, tags=None, transform=None, wanted=None, frame=True): if transform: points = np.array([transform(p) for p in points]) box = [float(v) for v in points.min(axis=0)] + [float(v) for v in points.max(axis=0)] - element = model.by_id(shape.id) out[element.Tag] = (element.is_a(), box, sources.get(element.Tag)) - if not it.next(): - break return out From 3ff51de3e856976527e7d10feb863c879a12b741 Mon Sep 17 00:00:00 2001 From: Griffin Long <151978260+DrunkOnJava@users.noreply.github.com> Date: Tue, 3 Nov 2026 12:35:47 -0500 Subject: [PATCH 2/3] ci(measure): the two mesh scorers run as one task per model and share what they meshed, through a cache Measure deletes after the model (#1237) mesh_bodies_vs_ifc and saved_meshes_vs_ifc mesh the same bodies of the same two files, so Revit's side was meshed twice. With RVT_MESH_CACHE set, mesh_parts keeps what it meshed in npz shards keyed by the file's path, size and mtime, and reads what an earlier scorer kept; the script chains the two scorers as one task per model, sets it to the model's folder and removes it with the model's other bulky outputs. Fewer than 16 parts are meshed without a pool. --- tools/ci/measure-reference-models.sh | 17 +++++-- tools/re/ifc_common.py | 74 +++++++++++++++++++++++++--- tools/re/mesh_bodies_vs_ifc.py | 2 +- tools/re/saved_meshes_vs_ifc.py | 2 +- 4 files changed, 84 insertions(+), 11 deletions(-) diff --git a/tools/ci/measure-reference-models.sh b/tools/ci/measure-reference-models.sh index c11ebf85..51fad878 100755 --- a/tools/ci/measure-reference-models.sh +++ b/tools/ci/measure-reference-models.sh @@ -256,7 +256,7 @@ finished() { # dir name started # longest ones do not start last (the BIM4LCA office's mesh_bodies_vs_ifc took # 477 s of its 1,953 run one after another, run 41589642386). model_tasks() { # name rvt ref - local name="$1" rvt="$2" ref="$3" s + local name="$1" rvt="$2" ref="$3" s meshes= echo "1 base" # ports_vs_ifc reads the findings the schema step's validation wrote, so it # runs after it, in the same task. @@ -268,8 +268,10 @@ model_tasks() { # name rvt ref wanted "$s" || continue [ "$s" = ports_vs_ifc ] && wanted schema && continue case "$s" in - mesh_bodies_vs_ifc) echo "10 ifc:$s" ;; - saved_meshes_vs_ifc) echo "8 ifc:$s" ;; + # These two mesh the same bodies of the same two files, so they run + # one after the other as one task and the second reads what the + # first meshed. + mesh_bodies_vs_ifc | saved_meshes_vs_ifc) meshes=1 ;; layer_sets_vs_ifc | wall_heights_vs_ifc | door_bodies_vs_ifc) echo "2 ifc:$s" ;; *) echo "1 ifc:$s" ;; esac @@ -280,6 +282,7 @@ model_tasks() { # name rvt ref echo "1 ifc:$s" done fi + [ -n "$meshes" ] && echo "18 meshes" [ -f "$MODELS/${rvt%.rvt}.revit-sets.ifc" ] && wanted revit_property_sets_vs_ifc && echo "2 sets" return 0 } @@ -324,6 +327,13 @@ model_task() { # name rvt ref task s=${task#ifc:} score "$dir/$s.txt" "$SCORERS/$s.py" "$dir/model.ifc" "$MODELS/$ref" ;; + meshes) + for s in saved_meshes_vs_ifc mesh_bodies_vs_ifc; do + wanted "$s" || continue + RVT_MESH_CACHE="$dir/mesh-cache" score "$dir/$s.txt" "$SCORERS/$s.py" "$dir/model.ifc" "$MODELS/$ref" + done + rm -rf "$dir/mesh-cache" + ;; sets) # RE-210: where Revit's export with its own property sets is beside the # model (.revit-sets.ifc, the R2025 models), an export with @@ -348,6 +358,7 @@ model_task() { # name rvt ref task model_done() { # name local dir ms dir=$(model_dir "$1") + rm -rf "$dir/mesh-cache" rm -f "$dir/model.glb" "$dir/validation.jsonl" [ -n "${KEEP_IFC:-}" ] || rm -f "$dir/model.ifc" ms=$(awk '{ms += $2} END {print ms + 0}' "$dir/timings.log" 2>/dev/null || echo 0) diff --git a/tools/re/ifc_common.py b/tools/re/ifc_common.py index cc55ec32..8b5e51c6 100644 --- a/tools/re/ifc_common.py +++ b/tools/re/ifc_common.py @@ -2,6 +2,7 @@ import concurrent.futures import concurrent.futures.process +import hashlib import multiprocessing import os import re @@ -45,6 +46,9 @@ def report_mesh_error(error): _MESHING = {} +# Fewer parts than this are meshed in this process: a forked worker parses +# again the parts of the file it meshes, which a handful do not repay. +POOL_FROM = 16 def _mesh_one(part_id): @@ -58,16 +62,69 @@ def _mesh_one(part_id): ) -def mesh_parts(model, settings, parts): +def _mesh_cache(path, settings): + root = os.environ.get("RVT_MESH_CACHE") + if not root or not path: + return None + stat = os.stat(path) + key = f"{os.path.realpath(path)} {stat.st_size} {stat.st_mtime_ns} {settings.get('use-world-coords')}" + return os.path.join(root, hashlib.sha256(key.encode()).hexdigest()[:16]) + + +def _read_meshes(path): + """What `_write_meshes` kept: instance id -> (vertices, faces), or None for a part IfcOpenShell could not mesh.""" + out = {} + with np.load(path, allow_pickle=False) as kept: + ids, meshed, nverts, nfaces = kept["ids"], kept["meshed"], kept["nverts"], kept["nfaces"] + verts, faces = kept["verts"], kept["faces"] + vstart = np.concatenate([[0], np.cumsum(nverts)]) + fstart = np.concatenate([[0], np.cumsum(nfaces)]) + for k, part_id in enumerate(ids): + out[int(part_id)] = (verts[vstart[k] : vstart[k + 1]], faces[fstart[k] : fstart[k + 1]]) if meshed[k] else None + return out + + +def _write_meshes(path, meshes): + ids = list(meshes) + parts = [meshes[i] if meshes[i] is not None else (np.empty((0, 3)), np.empty((0, 3), dtype=int)) for i in ids] + scratch = f"{path}.tmp" + with open(scratch, "wb") as f: + np.savez( + f, + ids=np.array(ids, dtype=np.int64), + meshed=np.array([meshes[i] is not None for i in ids], dtype=bool), + nverts=np.array([len(v) for v, _ in parts], dtype=np.int64), + nfaces=np.array([len(t) for _, t in parts], dtype=np.int64), + verts=np.concatenate([v for v, _ in parts]), + faces=np.concatenate([t for _, t in parts]), + ) + os.replace(scratch, path) + + +def mesh_parts(model, settings, parts, path=None): """Instance id -> (vertices (N, 3), faces (M, 3)) for each of `parts` that IfcOpenShell meshes, one `create_shape` per part spread over the CPUs. The triangles are those `create_shape` gives, in the same order (the iterator's are the same, but it took twice the CPU time per element in a measurement on the BIM4LCA office, and no less wall time with 12 threads - than with 4); a part it cannot mesh is left out.""" + than with 4); a part it cannot mesh is left out. + + With `RVT_MESH_CACHE` set to a directory and the model's `path`, what is + meshed is kept there for the next scorer to read, which Measure sets for + the scorers it runs one after another on one pair of files and deletes + afterwards.""" ids = list(dict.fromkeys(part.id() for part in parts)) + cache = _mesh_cache(path, settings) + known = {} + if cache and os.path.isdir(cache): + for name in sorted(os.listdir(cache)): + if name.endswith(".npz"): + known.update(_read_meshes(os.path.join(cache, name))) + todo = [i for i in ids if i not in known] cpus = len(os.sched_getaffinity(0)) if hasattr(os, "sched_getaffinity") else os.cpu_count() or 1 - workers = min(cpus, len(ids)) + workers = min(int(os.environ.get("RVT_MESH_WORKERS") or cpus), len(todo)) + if len(todo) < POOL_FROM: + workers = 1 _MESHING.update(model=model, settings=settings) meshes = None if workers > 1: @@ -76,10 +133,15 @@ def mesh_parts(model, settings, parts): sys.stderr.flush() try: with concurrent.futures.ProcessPoolExecutor(workers, mp_context=multiprocessing.get_context("fork")) as pool: - meshes = list(pool.map(_mesh_one, ids, chunksize=1)) + meshes = list(pool.map(_mesh_one, todo, chunksize=1)) except concurrent.futures.process.BrokenProcessPool: meshes = None if meshes is None: - meshes = [_mesh_one(i) for i in ids] + meshes = [_mesh_one(i) for i in todo] _MESHING.clear() - return {i: mesh for i, mesh in zip(ids, meshes) if mesh is not None} + fresh = dict(zip(todo, meshes)) + if cache and fresh: + os.makedirs(cache, exist_ok=True) + _write_meshes(os.path.join(cache, f"{os.getpid()}.npz"), fresh) + known.update(fresh) + return {i: known[i] for i in ids if known[i] is not None} diff --git a/tools/re/mesh_bodies_vs_ifc.py b/tools/re/mesh_bodies_vs_ifc.py index 96dfdda6..0bca0fd4 100644 --- a/tools/re/mesh_bodies_vs_ifc.py +++ b/tools/re/mesh_bodies_vs_ifc.py @@ -70,7 +70,7 @@ def meshes(path, tags=None, source=None, model=None): if source is not None and body_sources(element) != source: continue chosen.append((element, shapes_of(element))) - shapes = mesh_parts(model, settings, [part for _, parts in chosen for part in parts]) + shapes = mesh_parts(model, settings, [part for _, parts in chosen for part in parts], path) out = {} for element, parts in chosen: triangles = [] diff --git a/tools/re/saved_meshes_vs_ifc.py b/tools/re/saved_meshes_vs_ifc.py index e98db67e..f5ba01d4 100644 --- a/tools/re/saved_meshes_vs_ifc.py +++ b/tools/re/saved_meshes_vs_ifc.py @@ -85,7 +85,7 @@ def world_boxes(path, tags=None, transform=None, wanted=None, frame=True): # The order the geometry iterator returned them in, one thread, which the # tables' ties (equal counts) follow. elements.sort(key=lambda e: e.id()) - shapes = mesh_parts(model, settings, elements) + shapes = mesh_parts(model, settings, elements, path) for element in elements: if element.id() not in shapes: continue From 9851ae6c259ca7ebbe56ac12c56e935bef366c98 Mon Sep 17 00:00:00 2001 From: Griffin Long <151978260+DrunkOnJava@users.noreply.github.com> Date: Tue, 3 Nov 2026 13:42:14 -0500 Subject: [PATCH 3/3] ci(measure): the surface distance's pruning slack covers the barycentric tolerance of a very large triangle, and the mesher's workers, shards and cache key are firmer (#1237) The slack was 1e-9 of the extent, less than the 1e-9 tolerance of inside can move a point (up to 2e-9 of two edges), so a point 0.99e-9 outside a 1e7 ft triangle measured 0.003 where the exhaustive search says 0; it is 4e-9 now, with the argument in the docstring. A pool worker stops at its next part when its parent has died, a shard is named by process and time, the cache key holds every setting, a pool that breaks says so on stderr, a part IfcOpenShell cannot mesh keeps its reason in MESH_FAILURES (and in the cache), and the script logs the cache's size before deleting it. --- tools/ci/measure-reference-models.sh | 1 + tools/re/ifc_common.py | 56 +++++++++++++++++++++------- tools/re/mesh_bodies_vs_ifc.py | 30 +++++++++------ 3 files changed, 63 insertions(+), 24 deletions(-) diff --git a/tools/ci/measure-reference-models.sh b/tools/ci/measure-reference-models.sh index 84b8d910..c51071e0 100755 --- a/tools/ci/measure-reference-models.sh +++ b/tools/ci/measure-reference-models.sh @@ -334,6 +334,7 @@ model_task() { # name rvt ref task wanted "$s" || continue RVT_MESH_CACHE="$dir/mesh-cache" score "$dir/$s.txt" "$SCORERS/$s.py" "$dir/model.ifc" "$MODELS/$ref" done + echo "$name: mesh cache $(du -sh "$dir/mesh-cache" 2>/dev/null | cut -f1), removed" rm -rf "$dir/mesh-cache" ;; sets) diff --git a/tools/re/ifc_common.py b/tools/re/ifc_common.py index 8b5e51c6..06a6aca8 100644 --- a/tools/re/ifc_common.py +++ b/tools/re/ifc_common.py @@ -7,6 +7,7 @@ import os import re import sys +import time import ifcopenshell.geom import ifcopenshell.util.placement @@ -46,53 +47,75 @@ def report_mesh_error(error): _MESHING = {} +# Instance id -> the first line of IfcOpenShell's reason, for each part +# `mesh_parts` could not mesh (a part kept in the cache keeps its reason too). +# Nothing prints them: a scorer's output does not change when a part fails. +MESH_FAILURES = {} # Fewer parts than this are meshed in this process: a forked worker parses # again the parts of the file it meshes, which a handful do not repay. POOL_FROM = 16 +def _in_worker(): + _MESHING["parent"] = os.getppid() + + def _mesh_one(part_id): + """(vertices, faces) of the part, or the first line of why IfcOpenShell would not mesh it.""" + # An OOM-killed parent leaves its workers with the queued parts: stop at the next one. + if "parent" in _MESHING and os.getppid() != _MESHING["parent"]: + os._exit(1) try: shape = ifcopenshell.geom.create_shape(_MESHING["settings"], _MESHING["model"].by_id(part_id)) - except RuntimeError: - return None + except RuntimeError as error: + return (str(error).strip().splitlines() or [""])[0] or type(error).__name__ return ( np.array(shape.geometry.verts, dtype=float).reshape(-1, 3), np.array(shape.geometry.faces, dtype=int).reshape(-1, 3), ) +def _settings_key(settings): + items = [] + for name in sorted(settings.setting_names()): + try: + items.append(f"{name}={settings.get(name)!r}") + except RuntimeError: + items.append(f"{name}=unset") + return ";".join(items) + + def _mesh_cache(path, settings): root = os.environ.get("RVT_MESH_CACHE") if not root or not path: return None stat = os.stat(path) - key = f"{os.path.realpath(path)} {stat.st_size} {stat.st_mtime_ns} {settings.get('use-world-coords')}" + key = f"{os.path.realpath(path)} {stat.st_size} {stat.st_mtime_ns} {_settings_key(settings)}" return os.path.join(root, hashlib.sha256(key.encode()).hexdigest()[:16]) def _read_meshes(path): - """What `_write_meshes` kept: instance id -> (vertices, faces), or None for a part IfcOpenShell could not mesh.""" + """What `_write_meshes` kept: instance id -> (vertices, faces), or the reason IfcOpenShell gave for a part it could not mesh.""" out = {} with np.load(path, allow_pickle=False) as kept: - ids, meshed, nverts, nfaces = kept["ids"], kept["meshed"], kept["nverts"], kept["nfaces"] + ids, reasons, nverts, nfaces = kept["ids"], kept["reasons"], kept["nverts"], kept["nfaces"] verts, faces = kept["verts"], kept["faces"] vstart = np.concatenate([[0], np.cumsum(nverts)]) fstart = np.concatenate([[0], np.cumsum(nfaces)]) for k, part_id in enumerate(ids): - out[int(part_id)] = (verts[vstart[k] : vstart[k + 1]], faces[fstart[k] : fstart[k + 1]]) if meshed[k] else None + out[int(part_id)] = str(reasons[k]) or (verts[vstart[k] : vstart[k + 1]], faces[fstart[k] : fstart[k + 1]]) return out def _write_meshes(path, meshes): ids = list(meshes) - parts = [meshes[i] if meshes[i] is not None else (np.empty((0, 3)), np.empty((0, 3), dtype=int)) for i in ids] + parts = [(np.empty((0, 3)), np.empty((0, 3), dtype=int)) if isinstance(meshes[i], str) else meshes[i] for i in ids] scratch = f"{path}.tmp" with open(scratch, "wb") as f: np.savez( f, ids=np.array(ids, dtype=np.int64), - meshed=np.array([meshes[i] is not None for i in ids], dtype=bool), + reasons=np.array([meshes[i] if isinstance(meshes[i], str) else "" for i in ids]), nverts=np.array([len(v) for v, _ in parts], dtype=np.int64), nfaces=np.array([len(t) for _, t in parts], dtype=np.int64), verts=np.concatenate([v for v, _ in parts]), @@ -107,12 +130,15 @@ def mesh_parts(model, settings, parts, path=None): triangles are those `create_shape` gives, in the same order (the iterator's are the same, but it took twice the CPU time per element in a measurement on the BIM4LCA office, and no less wall time with 12 threads - than with 4); a part it cannot mesh is left out. + than with 4); a part it cannot mesh is left out, with its reason in + `MESH_FAILURES`. With `RVT_MESH_CACHE` set to a directory and the model's `path`, what is meshed is kept there for the next scorer to read, which Measure sets for the scorers it runs one after another on one pair of files and deletes - afterwards.""" + afterwards. Each call keeps its own shard, named by process and time, and + the cache is keyed by the file's path, size and modification time and by + every setting.""" ids = list(dict.fromkeys(part.id() for part in parts)) cache = _mesh_cache(path, settings) known = {} @@ -132,9 +158,12 @@ def mesh_parts(model, settings, parts, path=None): sys.stdout.flush() sys.stderr.flush() try: - with concurrent.futures.ProcessPoolExecutor(workers, mp_context=multiprocessing.get_context("fork")) as pool: + with concurrent.futures.ProcessPoolExecutor( + workers, mp_context=multiprocessing.get_context("fork"), initializer=_in_worker + ) as pool: meshes = list(pool.map(_mesh_one, todo, chunksize=1)) except concurrent.futures.process.BrokenProcessPool: + print("mesh_parts: a worker died, so the parts are meshed one after another", file=sys.stderr) meshes = None if meshes is None: meshes = [_mesh_one(i) for i in todo] @@ -142,6 +171,7 @@ def mesh_parts(model, settings, parts, path=None): fresh = dict(zip(todo, meshes)) if cache and fresh: os.makedirs(cache, exist_ok=True) - _write_meshes(os.path.join(cache, f"{os.getpid()}.npz"), fresh) + _write_meshes(os.path.join(cache, f"{os.getpid()}-{time.time_ns()}.npz"), fresh) known.update(fresh) - return {i: known[i] for i in ids if known[i] is not None} + MESH_FAILURES.update({i: known[i] for i in ids if isinstance(known[i], str)}) + return {i: known[i] for i in ids if not isinstance(known[i], str)} diff --git a/tools/re/mesh_bodies_vs_ifc.py b/tools/re/mesh_bodies_vs_ifc.py index 0bca0fd4..5fcece3d 100644 --- a/tools/re/mesh_bodies_vs_ifc.py +++ b/tools/re/mesh_bodies_vs_ifc.py @@ -182,16 +182,24 @@ def surface_distance(points, triangles): `distances` over a set of triangles that holds every one that could be nearest, and each of them is measured by the very arithmetic that measures all of them in a small job (`EXHAUSTIVE` pairs or fewer, and any input that is not finite), so - the result is the same float. A triangle is left out for a point only when a lower - bound on its distance exceeds the nearest distance already found for that point - plus a slack (1e-6 ft and 1e-9 of the coordinates' extent) that is far above the - rounding of these distances and the 1e-9 barycentric tolerance of `inside`. The - bounds hold for every triangle but the slivers (the sine of the angle between two - edges under SLIVER, whose computed plane is not reliable) and the zero-area ones (which - have no plane): the bounds are the distance to a block of eight neighbouring - triangles' box, then to a triangle's box, and for a triangle with a plane the - distance to its plane and to the side planes of its edges. Slivers are measured - against every point. + the result is the same float. + + A triangle is left out for a point only when a lower bound on its distance + exceeds the nearest distance already found for that point plus a slack. The + slack is 1e-6 ft and 4e-9 of the extent, the largest of any coordinate and any + triangle edge. It has to cover two things. `inside` takes a point up to 1e-9 of + two edges outside the triangle in its plane (or 1e-9 past the third) as inside + and measures it by its plane alone, so the distance computed for it falls short + of the distance to the triangle, which the bounds are on, by up to 1e-9 times + the sum of two edges: under 2e-9 of the extent. And the rounding of these + distances, which is many orders of magnitude smaller. + + The bounds hold for every triangle but the slivers (the sine of the angle between + two edges under SLIVER, whose computed plane is not reliable) and the zero-area + ones (which have no plane, so only their box bounds them): the distance to a + block of eight neighbouring triangles' box, then to a triangle's box, and for a + triangle with a plane the distance to its plane combined with that to the side + planes of its edges. Slivers are measured against every point. """ terms = triangle_terms(triangles) n = len(points) @@ -208,7 +216,7 @@ def surface_distance(points, triangles): lo = np.minimum(np.minimum(a, b), c) hi = np.maximum(np.maximum(a, b), c) extent = max(float(np.abs(triangles).max()), float(np.abs(points).max()), float(np.sqrt(max(d00.max(), d11.max())))) - slack = 1e-6 + 1e-9 * extent + slack = 1e-6 + 4e-9 * extent sides = [] for start, end in ((a, b), (b, c), (c, a)): out = np.cross(end - start, unit)