From 9f86845c255bde15fd719f6f8324c954e8fd2fa0 Mon Sep 17 00:00:00 2001 From: stroblme Date: Tue, 11 Aug 2026 16:26:20 +0200 Subject: [PATCH] Add bulk UV sampler, TKHelix and MakeSphere (8.0.1.2) - n3xd_ocp.sample.face_grid: a face's UV grid of points and outward normals in one GIL-free call, sampling exactly where np.linspace does so a caller's fitted surface does not move. - n3xd_ocp.helix: OCCT 8.0's TKHelix, which upstream does not bind at all. Takes Python lists rather than NCollection_Array1. Two things the header does not say, both found by probing: SetParameters wants N+1 diameters for N segments (one per boundary, so a taper interpolates), and the builder is right-hand only -- a negative pitch is error status 12, not a mirrored helix. - Bind BRepPrimAPI_MakeSphere and give inventory.py an EXTRA_SYMBOLS addendum for symbols no app source imports. assay's gen_flow_fixtures has been unrunnable since the 10C cutover for want of it; the gap was wider than --check, since sigdiff is inventory-driven too. Gates: 105 tests, 139/139 symbols, sigdiff clean, ASAN clean, wheel self-contained with no libGL/libX11 DT_NEEDED. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01HbTQ2HYWQwdtGmGJwypt6Z --- CMakeLists.txt | 2 + README.md | 16 ++++ docs/design.md | 25 +++++- inventory.json | 7 +- pyproject.toml | 2 +- python/n3xd_ocp/__init__.py | 4 +- src/ext/ext_core.cpp | 4 + src/ext/ext_helix.cpp | 125 ++++++++++++++++++++++++++++++ src/ext/ext_sample.cpp | 122 +++++++++++++++++++++++++++++ src/modules/mod_BRepPrimAPI.cpp | 10 +++ tests/test_ext_helix.py | 118 ++++++++++++++++++++++++++++ tests/test_ext_sample.py | 133 ++++++++++++++++++++++++++++++++ tools/inventory.py | 13 ++++ 13 files changed, 573 insertions(+), 8 deletions(-) create mode 100644 src/ext/ext_helix.cpp create mode 100644 src/ext/ext_sample.cpp create mode 100644 tests/test_ext_helix.py create mode 100644 tests/test_ext_sample.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 172db9c..2edd054 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -77,6 +77,8 @@ target_link_libraries(_OCP PRIVATE # image nothing. TKService TKV3d + # 10E — OCCT 8.0's helix toolkit, reached only through n3xd_ocp.helix. + TKHelix ) target_compile_definitions(_OCP PRIVATE diff --git a/README.md b/README.md index 7358cae..b013323 100644 --- a/README.md +++ b/README.md @@ -51,6 +51,22 @@ import n3xd_ocp areas, centroids = n3xd_ocp.measure.face_surface_props(result) # one call for every face meshes = n3xd_ocp.tess.extract_meshes(result) # triangulated faces, ready to render data = n3xd_ocp.bintools.write_bytes(result) # BREP bytes, no temp file needed + +points, normals, uv_bounds = n3xd_ocp.sample.face_grid(face, 33) # a 33x33 UV grid on one face +``` + +`n3xd_ocp.helix` reaches OCCT 8.0's TKHelix, which upstream has no binding for. For N segments it wants N pitches, N turn counts and **N+1 diameters** — one per segment boundary, so consecutive values that differ taper across that segment: + +```python +from OCP.gp import gp_Ax3, gp_Dir, gp_Pnt + +axis = gp_Ax3(gp_Pnt(0, 0, 0), gp_Dir(0, 0, 1), gp_Dir(1, 0, 0)) +wire, tolerance_reached = n3xd_ocp.helix.pure_helix(axis, 8.0, [1.25], [12.0]) + +builder = n3xd_ocp.helix.BuilderHelix() # tapered, e.g. an NPT thread +builder.set_parameters(axis, [10.0, 8.0], [2.0], [4.0]) +builder.set_approx_parameters(1.0e-4) +builder.perform() ``` ## Build diff --git a/docs/design.md b/docs/design.md index 5ca5735..40e0a1a 100644 --- a/docs/design.md +++ b/docs/design.md @@ -162,9 +162,28 @@ would put a bare module in `sys.modules`, and a later `import n3xd_ocp` would then skip its `__init__.py`. Shipped today: `bintools` (shape ↔ `bytes`, byte-identical to -`OCP.BinTools`), `measure` (batched per-face area and centroid) and `tess` -(triangulated meshes and edge polylines) — see their `.pyi` stubs for full -signatures, or the [usage examples](../README.md#usage) for a quick start. +`OCP.BinTools`), `measure` (batched per-face area and centroid), `tess` +(triangulated meshes and edge polylines), `sample` (a face's UV grid of +points and normals) and `helix` (OCCT 8.0's TKHelix builder) — see their +`.pyi` stubs for full signatures, or the +[usage examples](../README.md#usage) for a quick start. + +The array extractors all reproduce the Python loop they replace exactly, +including its quirks, because their callers key on them: `tess` keeps the id +gap left by an untriangulated face, and `sample` samples the way +`numpy.linspace` does — endpoint forced onto the bound rather than +`start + (n-1)*step` — so the surface a caller fits does not move. + +`helix` is the one module with no upstream counterpart to match: TKHelix is +new in OCCT 8.0 and `cadquery-ocp` is still on 7.9.3, so a 1:1 surface would +be invented rather than reproduced. It takes Python lists and builds the +`NCollection_Array1` internally, which also avoids binding that template for +one caller. Its parameter shape is not in the OCCT header and cost a probe to +find: for **N segments** `SetParameters` wants N pitches, N turn counts and +**N+1 diameters**, one per segment *boundary*, so consecutive diameters that +differ taper across that segment. Anything else is a +`Standard_ConstructionError` reading only "wrong array dimension", so the +binding checks the shape first and says what it wanted. ## Matching upstream, and how that's checked diff --git a/inventory.json b/inventory.json index ebc908d..fcfe858 100644 --- a/inventory.json +++ b/inventory.json @@ -1,6 +1,6 @@ { - "app_root": "/home/stroblme/n3xd/app/backend", - "files_importing_ocp": 115, + "app_root": "../app/backend", + "files_importing_ocp": 116, "modules": { "BOPAlgo": { "BOPAlgo_GlueShift": [] @@ -86,7 +86,8 @@ "BRepPrimAPI_MakeCylinder": [], "BRepPrimAPI_MakeHalfSpace": [], "BRepPrimAPI_MakePrism": [], - "BRepPrimAPI_MakeRevol": [] + "BRepPrimAPI_MakeRevol": [], + "BRepPrimAPI_MakeSphere": [] }, "BRepTools": { "BRepTools": [ diff --git a/pyproject.toml b/pyproject.toml index b9d3a6e..6207f16 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,7 +9,7 @@ name = "n3xd-ocp" # .devN suffix: the registry never allows republishing a version, and its # cleanup rule collects only the dev ones. CMake asserts the prefix matches # the OCCT it found. -version = "8.0.1.1" +version = "8.0.1.2" description = "nanobind bindings for the OpenCASCADE geometry kernel (drop-in OCP)" readme = "README.md" requires-python = ">=3.12" diff --git a/python/n3xd_ocp/__init__.py b/python/n3xd_ocp/__init__.py index 5a24113..2e3ccae 100644 --- a/python/n3xd_ocp/__init__.py +++ b/python/n3xd_ocp/__init__.py @@ -16,8 +16,10 @@ import sys import OCP as _OCP_pkg # noqa: F401 (loads the extension, registering us) bintools = sys.modules["n3xd_ocp.bintools"] +helix = sys.modules["n3xd_ocp.helix"] measure = sys.modules["n3xd_ocp.measure"] +sample = sys.modules["n3xd_ocp.sample"] tess = sys.modules["n3xd_ocp.tess"] _debug = sys.modules["n3xd_ocp._debug"] -__all__ = ["bintools", "measure", "tess"] +__all__ = ["bintools", "helix", "measure", "sample", "tess"] diff --git a/src/ext/ext_core.cpp b/src/ext/ext_core.cpp index e2a98f8..7d36f8f 100644 --- a/src/ext/ext_core.cpp +++ b/src/ext/ext_core.cpp @@ -18,12 +18,16 @@ namespace nb = nanobind; void register_ext_bintools(); void register_ext_debug(); +void register_ext_helix(); void register_ext_measure(); +void register_ext_sample(); void register_ext_tess(); void register_ext(nb::module_ &) { register_ext_bintools(); register_ext_debug(); + register_ext_helix(); register_ext_measure(); + register_ext_sample(); register_ext_tess(); } diff --git a/src/ext/ext_helix.cpp b/src/ext/ext_helix.cpp new file mode 100644 index 0000000..46a6358 --- /dev/null +++ b/src/ext/ext_helix.cpp @@ -0,0 +1,125 @@ +/* + n3xd_ocp.helix — OCCT 8.0's helix builder, over Python lists. + + Lives here rather than in OCP.HelixBRep because there is nothing upstream + to match: cadquery-ocp is on 7.9.3, which has no TKHelix at all, so a 1:1 + surface would be invented rather than reproduced. When upstream does bind + it, OCP.HelixBRep can be added beside this and the two can disagree + harmlessly. + + HelixBRep_BuilderHelix::SetParameters takes NCollection_Array1, which is + not bound; the vectors are converted here instead. The array-valued form + is the one bound because it is the one that carries variable pitch and + variable diameter — tapered threads and variable-pitch springs — which the + scalar overloads cannot express. +*/ + +#include "../common/occt_module.h" +#include "../common/occt_policies.h" + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace { + +/// Copy a std::vector into the Array1 the builder wants. +template +NCollection_Array1 to_array1(const std::vector &values) { + NCollection_Array1 array(1, (int) values.size()); + for (size_t i = 0; i < values.size(); ++i) + array.SetValue((int) (i + 1), values[i]); + return array; +} + +} // namespace + +void register_ext_helix() { + nb::module_ m = ocp_named_module("n3xd_ocp.helix"); + + nb::class_(m, "BuilderHelix", + "OCCT's composite-helix wire builder.") + .def(nb::init<>()) + .def( + "set_parameters", + [](HelixBRep_BuilderHelix &self, const gp_Ax3 &axis, + const std::vector &diameters, + const std::vector &pitches, + const std::vector &nb_turns) { + if (pitches.empty() || pitches.size() != nb_turns.size() || + diameters.size() != pitches.size() + 1) + throw std::invalid_argument( + "set_parameters: for N segments OCCT wants N pitches, " + "N nb_turns and N+1 diameters (one per segment " + "boundary); got " + + std::to_string(diameters.size()) + ", " + + std::to_string(pitches.size()) + ", " + + std::to_string(nb_turns.size())); + self.SetParameters(axis, to_array1(diameters), + to_array1(pitches), to_array1(nb_turns)); + }, + "axis"_a, "diameters"_a, "pitches"_a, "nb_turns"_a, + R"doc(Parameters of an N-segment composite helix. + +``pitches`` and ``nb_turns`` carry one entry per segment; ``diameters`` carries +N+1, one per segment *boundary*, so consecutive values that differ produce a +taper across that segment. OCCT raises ``Standard_ConstructionError`` for any +other combination, so the shape is checked here first.)doc") + .def( + "set_approx_parameters", + [](HelixBRep_BuilderHelix &self, double tolerance, int max_degree, + GeomAbs_Shape continuity) { + self.SetApproxParameters(tolerance, max_degree, continuity); + }, + "tolerance"_a, "max_degree"_a = 8, "continuity"_a = GeomAbs_C1) + .def("perform", &HelixBRep_BuilderHelix::Perform, OCP_NOGIL) + .def("tolerance_reached", &HelixBRep_BuilderHelix::ToleranceReached) + .def("error_status", &HelixBRep_BuilderHelix::ErrorStatus) + .def("warning_status", &HelixBRep_BuilderHelix::WarningStatus) + .def("shape", &HelixBRep_BuilderHelix::Shape, OCP_RETURN_COPY); + + m.def( + "pure_helix", + [](const gp_Ax3 &axis, double diameter, + const std::vector &pitches, + const std::vector &nb_turns, double tolerance, + GeomAbs_Shape continuity) { + if (pitches.size() != nb_turns.size() || pitches.empty()) + throw std::invalid_argument( + "pure_helix: pitches and nb_turns must be the same " + "non-empty length"); + // N+1 boundary diameters, all equal: a constant-diameter helix. + const std::vector diameters(pitches.size() + 1, diameter); + + HelixBRep_BuilderHelix builder; + builder.SetParameters(axis, to_array1(diameters), + to_array1(pitches), to_array1(nb_turns)); + builder.SetApproxParameters(tolerance, 8, continuity); + { + nb::gil_scoped_release nogil; + builder.Perform(); + } + const int status = builder.ErrorStatus(); + if (status != 0) + throw std::runtime_error("pure_helix: HelixBRep_BuilderHelix " + "failed with error status " + + std::to_string(status)); + return nb::make_tuple(builder.Shape(), builder.ToleranceReached()); + }, + "axis"_a, "diameter"_a, "pitches"_a, "nb_turns"_a, nb::kw_only(), + "tolerance"_a = 1.0e-4, "continuity"_a = GeomAbs_C1, + R"doc(A constant-diameter helix wire, as ``(wire, tolerance_reached)``. + +One ``pitches``/``nb_turns`` pair per segment, so a single-segment call is the +ordinary thread spine. Raises ``RuntimeError`` when the approximation fails — +the returned tolerance says how close it actually got, which is what decides +whether the result is usable as a sweep spine.)doc"); +} diff --git a/src/ext/ext_sample.cpp b/src/ext/ext_sample.cpp new file mode 100644 index 0000000..c31f41f --- /dev/null +++ b/src/ext/ext_sample.cpp @@ -0,0 +1,122 @@ +/* + n3xd_ocp.sample — a face's UV grid, without the per-sample round trip. + + cad/operations/surface_pattern.py walks an (res+1)² UV grid in nested + Python loops: one BRepAdaptor_Surface::D1 plus roughly ten accessor calls + per sample, and the cross product in Python. The kernel can fill the whole + grid in one call with the GIL released. + + The output reproduces what the Python loop produces: + + * the grid is BRepTools::UVBounds, sampled the way numpy's linspace + samples it — start + i*step with the endpoint forced exactly, so the + sample locations are bit-identical and the fitted surface does not + move. The bounds come back with the arrays for the same reason: the + caller derives its normalised u/v from these numbers rather than + recomputing them. + * normals are du × dv normalised, negated for a reversed face. + * a sample whose cross product is degenerate (a pole or apex) gets the + zero vector, not a unit vector in an arbitrary direction. Callers + test it — `_inward_thickness` skips a sample on `nrm.any()`. +*/ + +#include "../common/occt_module.h" + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace { + +template nb::capsule owner_of(T *data) { + return nb::capsule(data, [](void *p) noexcept { delete[] (T *) p; }); +} + +/// numpy's linspace endpoint handling: the last sample is the bound itself, +/// not start + (n-1)*step, which can miss it by an ulp. +inline double linspace_at(double lo, double hi, int i, int n) { + if (n <= 1) + return lo; + return i == n - 1 ? hi : lo + (double) i * (hi - lo) / (double) (n - 1); +} + +constexpr double MIN_NORMAL = 1e-12; + +} // namespace + +void register_ext_sample() { + nb::module_ m = ocp_named_module("n3xd_ocp.sample"); + + m.def( + "face_grid", + [](const TopoDS_Face &face, int n) { + if (n < 1) + throw std::invalid_argument("face_grid: n must be >= 1"); + + double umin = 0.0, umax = 0.0, vmin = 0.0, vmax = 0.0; + BRepTools::UVBounds(face, umin, umax, vmin, vmax); + const double sign = + face.Orientation() == TopAbs_REVERSED ? -1.0 : 1.0; + + const size_t count = (size_t) n * (size_t) n; + double *points = new double[3 * count]; + double *normals = new double[3 * count]; + + { + nb::gil_scoped_release nogil; + BRepAdaptor_Surface adaptor(face); + gp_Pnt point; + gp_Vec du, dv; + for (int i = 0; i < n; ++i) { + const double u = linspace_at(umin, umax, i, n); + for (int j = 0; j < n; ++j) { + const double v = linspace_at(vmin, vmax, j, n); + adaptor.D1(u, v, point, du, dv); + + const size_t at = 3 * ((size_t) i * (size_t) n + j); + points[at + 0] = point.X(); + points[at + 1] = point.Y(); + points[at + 2] = point.Z(); + + const double nx = du.Y() * dv.Z() - du.Z() * dv.Y(); + const double ny = du.Z() * dv.X() - du.X() * dv.Z(); + const double nz = du.X() * dv.Y() - du.Y() * dv.X(); + const double mag = std::sqrt(nx * nx + ny * ny + nz * nz); + if (mag > MIN_NORMAL) { + normals[at + 0] = sign * nx / mag; + normals[at + 1] = sign * ny / mag; + normals[at + 2] = sign * nz / mag; + } else { + normals[at + 0] = 0.0; + normals[at + 1] = 0.0; + normals[at + 2] = 0.0; + } + } + } + } + + size_t shape[3] = {(size_t) n, (size_t) n, 3}; + return nb::make_tuple( + nb::ndarray>(points, 3, shape, + owner_of(points)), + nb::ndarray>(normals, 3, shape, + owner_of(normals)), + nb::make_tuple(umin, umax, vmin, vmax)); + }, + "face"_a, "n"_a, + R"doc(Sample *face* on an n x n grid over its UV bounds. + +Returns ``(points[n, n, 3], normals[n, n, 3], (umin, umax, vmin, vmax))`` as +float64 arrays, u along axis 0. Normals are outward (negated for a reversed +face) and unit length, or the zero vector where the surface derivatives are +degenerate. Sample locations match ``numpy.linspace`` over the returned +bounds exactly.)doc"); +} diff --git a/src/modules/mod_BRepPrimAPI.cpp b/src/modules/mod_BRepPrimAPI.cpp index 8ee354c..4cebbae 100644 --- a/src/modules/mod_BRepPrimAPI.cpp +++ b/src/modules/mod_BRepPrimAPI.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -69,6 +70,15 @@ void register_BRepPrimAPI(nb::module_ &root) { const Standard_Real>(), "Axes"_a, "R1"_a, "R2"_a, "H"_a); + // Not reached by the app, so inventory.py never turned it up. assay's + // gen_flow_fixtures.py builds spheres, and has been unrunnable against the + // app's venv since the 10C cutover without this. + nb::class_( + m, "BRepPrimAPI_MakeSphere") + .def(nb::init(), "R"_a) + .def(nb::init(), "Center"_a, "R"_a) + .def(nb::init(), "Axis"_a, "R"_a); + nb::class_( m, "BRepPrimAPI_MakeHalfSpace") .def(nb::init(), "Face"_a, diff --git a/tests/test_ext_helix.py b/tests/test_ext_helix.py new file mode 100644 index 0000000..59e9484 --- /dev/null +++ b/tests/test_ext_helix.py @@ -0,0 +1,118 @@ +"""n3xd_ocp.helix — OCCT 8.0's TKHelix builder. + +The toolkit is new in 8.0 and has no upstream binding to compare against, so +these are property checks rather than a parity gate: the wire has to be a real +helix of the requested pitch and turn count, and the failure modes have to +surface as Python exceptions rather than a silently empty shape. +""" + +from __future__ import annotations + +import math + +import pytest + +import n3xd_ocp +from OCP.BRepGProp import BRepGProp +from OCP.GeomAbs import GeomAbs_C1, GeomAbs_C2 +from OCP.gp import gp_Ax3, gp_Dir, gp_Pnt +from OCP.GProp import GProp_GProps +from OCP.TopAbs import TopAbs_EDGE +from OCP.TopExp import TopExp +from OCP.TopTools import TopTools_IndexedMapOfShape + +AXIS = gp_Ax3(gp_Pnt(0, 0, 0), gp_Dir(0, 0, 1), gp_Dir(1, 0, 0)) + + +def _length(shape) -> float: + props = GProp_GProps() + BRepGProp.LinearProperties_s(shape, props) + return props.Mass() + + +def _edge_count(shape) -> int: + edges = TopTools_IndexedMapOfShape() + TopExp.MapShapes_s(shape, TopAbs_EDGE, edges) + return edges.Extent() + + +def test_single_segment_has_the_right_arc_length(): + """One turn of a helix is sqrt(circumference^2 + pitch^2) long.""" + diameter, pitch, turns = 8.0, 1.25, 1.0 + wire, reached = n3xd_ocp.helix.pure_helix(AXIS, diameter, [pitch], [turns]) + expected = turns * math.hypot(math.pi * diameter, pitch) + assert _length(wire) == pytest.approx(expected, rel=1e-4) + assert reached > 0.0 + assert _edge_count(wire) >= 1 + + +@pytest.mark.parametrize("turns", [0.5, 3.0, 24.0]) +def test_arc_length_scales_with_the_turn_count(turns): + diameter, pitch = 8.0, 1.25 + wire, _reached = n3xd_ocp.helix.pure_helix(AXIS, diameter, [pitch], [turns]) + expected = turns * math.hypot(math.pi * diameter, pitch) + assert _length(wire) == pytest.approx(expected, rel=1e-4) + + +def test_variable_pitch_segments_compose(): + """The array form is the reason this is bound at all — a spring whose + pitch changes partway is one wire, not two.""" + diameter = 10.0 + wire, _reached = n3xd_ocp.helix.pure_helix( + AXIS, diameter, [2.0, 5.0], [3.0, 2.0] + ) + expected = 3.0 * math.hypot(math.pi * diameter, 2.0) + 2.0 * math.hypot( + math.pi * diameter, 5.0 + ) + assert _length(wire) == pytest.approx(expected, rel=1e-3) + + +def test_tapered_diameters_via_the_builder(): + """Two segments, so three boundary diameters: 10 -> 8 tapering, then 8 -> 8 + straight. A tapered thread (NPT) is the reason the array form is bound.""" + builder = n3xd_ocp.helix.BuilderHelix() + builder.set_parameters(AXIS, [10.0, 8.0, 8.0], [2.0, 2.0], [2.0, 2.0]) + builder.set_approx_parameters(1.0e-4, 8, GeomAbs_C1) + builder.perform() + assert builder.error_status() == 0 + # The straight half is exact; the tapered half is longer than a cylinder of + # its smaller diameter and shorter than one of its larger. + straight = 2.0 * math.hypot(math.pi * 8.0, 2.0) + largest = 2.0 * math.hypot(math.pi * 10.0, 2.0) + total = _length(builder.shape()) + assert straight * 2 < total < straight + largest + + +def test_the_n_plus_one_diameter_rule_is_enforced(): + """OCCT wants N+1 diameters for N segments and raises + Standard_ConstructionError otherwise; the binding checks it first so the + message names the shape it wanted.""" + builder = n3xd_ocp.helix.BuilderHelix() + with pytest.raises(ValueError, match="N\\+1 diameters"): + builder.set_parameters(AXIS, [8.0], [1.25], [1.0]) + builder.set_parameters(AXIS, [8.0, 8.0], [1.25], [1.0]) # accepted + + +def test_continuity_is_selectable(): + wire_c1, _ = n3xd_ocp.helix.pure_helix( + AXIS, 8.0, [1.25], [2.0], continuity=GeomAbs_C1 + ) + wire_c2, _ = n3xd_ocp.helix.pure_helix( + AXIS, 8.0, [1.25], [2.0], continuity=GeomAbs_C2 + ) + assert _length(wire_c1) == pytest.approx(_length(wire_c2), rel=1e-3) + + +def test_mismatched_segment_arrays_are_rejected(): + with pytest.raises(ValueError): + n3xd_ocp.helix.pure_helix(AXIS, 8.0, [1.0, 2.0], [1.0]) + with pytest.raises(ValueError): + n3xd_ocp.helix.pure_helix(AXIS, 8.0, [], []) + + +def test_builder_rejects_ragged_parameters(): + builder = n3xd_ocp.helix.BuilderHelix() + with pytest.raises(ValueError): + builder.set_parameters(AXIS, [8.0, 8.0], [1.0], [1.0, 2.0]) + with pytest.raises(ValueError): + builder.set_parameters(AXIS, [8.0, 8.0], [], []) diff --git a/tests/test_ext_sample.py b/tests/test_ext_sample.py new file mode 100644 index 0000000..fc5d3a2 --- /dev/null +++ b/tests/test_ext_sample.py @@ -0,0 +1,133 @@ +"""n3xd_ocp.sample — the bulk UV-grid sampler. + +The reference is the nested Python loop in cad/operations/surface_pattern.py's +``_sample_grid``, reproduced below. It has to match exactly rather than +closely: the samples feed a B-spline fit, so a last-ulp difference in a sample +location moves the patterned surface and would make the corpus disagree with +itself for no reason. +""" + +from __future__ import annotations + +import pytest + +np = pytest.importorskip("numpy") + +import n3xd_ocp +from OCP.BRepAdaptor import BRepAdaptor_Surface +from OCP.BRepPrimAPI import BRepPrimAPI_MakeBox, BRepPrimAPI_MakeCylinder +from OCP.BRepTools import BRepTools +from OCP.gp import gp_Ax2, gp_Dir, gp_Pnt, gp_Vec +from OCP.TopAbs import TopAbs_FACE, TopAbs_REVERSED +from OCP.TopoDS import TopoDS + +from .test_inc1_modeling import sub_shapes + +_MIN_NORMAL = 1e-12 + + +def _faces(shape): + return [TopoDS.Face_s(f) for f in sub_shapes(shape, TopAbs_FACE)] + + +def _python_grid(face, n): + """surface_pattern.py::_sample_grid, reduced to points and normals.""" + adaptor = BRepAdaptor_Surface(face) + umin, umax, vmin, vmax = BRepTools.UVBounds_s(face) + sign = -1.0 if face.Orientation() == TopAbs_REVERSED else 1.0 + us = np.linspace(umin, umax, n) + vs = np.linspace(vmin, vmax, n) + point = np.zeros((n, n, 3)) + normal = np.zeros((n, n, 3)) + pnt, du, dv = gp_Pnt(), gp_Vec(), gp_Vec() + for i in range(n): + up = float(us[i]) + for j in range(n): + adaptor.D1(up, float(vs[j]), pnt, du, dv) + point[i, j] = (pnt.X(), pnt.Y(), pnt.Z()) + nx = du.Y() * dv.Z() - du.Z() * dv.Y() + ny = du.Z() * dv.X() - du.X() * dv.Z() + nz = du.X() * dv.Y() - du.Y() * dv.X() + mag = (nx * nx + ny * ny + nz * nz) ** 0.5 + if mag > _MIN_NORMAL: + normal[i, j] = (sign * nx / mag, sign * ny / mag, sign * nz / mag) + return point, normal, (umin, umax, vmin, vmax) + + +@pytest.fixture(scope="module") +def shapes(): + return { + "box": BRepPrimAPI_MakeBox(10.0, 20.0, 30.0).Shape(), + "cylinder": BRepPrimAPI_MakeCylinder( + gp_Ax2(gp_Pnt(0, 0, 0), gp_Dir(0, 0, 1)), 5.0, 12.0 + ).Shape(), + } + + +@pytest.mark.parametrize("name", ["box", "cylinder"]) +@pytest.mark.parametrize("n", [5, 33, 65]) +def test_matches_the_python_loop_exactly(shapes, name, n): + for face in _faces(shapes[name]): + points, normals, bounds = n3xd_ocp.sample.face_grid(face, n) + ref_points, ref_normals, ref_bounds = _python_grid(face, n) + assert bounds == ref_bounds + assert np.array_equal(points, ref_points) + assert np.array_equal(normals, ref_normals) + + +def test_reversed_face_flips_the_normal(shapes): + """A box has both orientations; the sign has to follow each face.""" + seen = set() + for face in _faces(shapes["box"]): + _points, normals, _bounds = n3xd_ocp.sample.face_grid(face, 5) + _rp, ref_normals, _rb = _python_grid(face, 5) + seen.add(face.Orientation() == TopAbs_REVERSED) + assert np.array_equal(normals, ref_normals) + assert seen == {True, False}, "fixture must cover both orientations" + + +def test_normals_are_unit_length(shapes): + for face in _faces(shapes["cylinder"]): + _points, normals, _bounds = n3xd_ocp.sample.face_grid(face, 17) + lengths = np.linalg.norm(normals.reshape(-1, 3), axis=1) + assert np.allclose(lengths, 1.0, atol=1e-12) + + +def test_endpoints_hit_the_bounds_exactly(shapes): + """linspace forces the last sample onto the bound; start + (n-1)*step can + miss it by an ulp, which is enough to move a fitted surface.""" + face = _faces(shapes["cylinder"])[0] + n = 33 + points, _normals, (umin, umax, vmin, vmax) = n3xd_ocp.sample.face_grid(face, n) + adaptor = BRepAdaptor_Surface(face) + for i, u in ((0, umin), (n - 1, umax)): + for j, v in ((0, vmin), (n - 1, vmax)): + expected = adaptor.Value(u, v) + assert points[i, j, 0] == expected.X() + assert points[i, j, 1] == expected.Y() + assert points[i, j, 2] == expected.Z() + + +def test_single_sample_is_the_lower_bound(shapes): + face = _faces(shapes["box"])[0] + points, _normals, (umin, _umax, vmin, _vmax) = n3xd_ocp.sample.face_grid(face, 1) + assert points.shape == (1, 1, 3) + expected = BRepAdaptor_Surface(face).Value(umin, vmin) + assert points[0, 0, 0] == expected.X() + + +def test_rejects_a_non_positive_grid(shapes): + with pytest.raises(ValueError): + n3xd_ocp.sample.face_grid(_faces(shapes["box"])[0], 0) + + +def test_arrays_own_their_memory(shapes): + import gc + + points, normals, _bounds = n3xd_ocp.sample.face_grid( + _faces(BRepPrimAPI_MakeBox(1.0, 2.0, 3.0).Shape())[0], 9 + ) + expected = float(points.sum()) + gc.collect() + assert float(points.sum()) == expected + assert normals.base is not None # the capsule, not the face diff --git a/tools/inventory.py b/tools/inventory.py index b9a8111..bc6a3d3 100644 --- a/tools/inventory.py +++ b/tools/inventory.py @@ -32,6 +32,15 @@ HERE = pathlib.Path(__file__).resolve().parent DEFAULT_INVENTORY = HERE.parent / "inventory.json" DEFAULT_APP = HERE.parent.parent / "app" / "backend" +#: Symbols no app source imports, so the scan below cannot see them, but which +#: another consumer of this wheel needs. Hand-maintained: add the symbol and +#: name who wants it, so a later reader can tell a real dependency from a +#: leftover. Merged into the scan, so --emit keeps them. +EXTRA_SYMBOLS: dict[str, dict[str, list[str]]] = { + # assay/scripts/gen_flow_fixtures.py builds spheres for its CFD fixtures. + "BRepPrimAPI": {"BRepPrimAPI_MakeSphere": []}, +} + def scan(app_root: pathlib.Path) -> dict: """Collect {module: {symbol: [attributes reached on it]}} from the app.""" @@ -70,6 +79,10 @@ def scan(app_root: pathlib.Path) -> dict: mod, name = local[node.value.id] modules[mod][name].add(node.attr) + for mod, symbols in EXTRA_SYMBOLS.items(): + for name, attrs in symbols.items(): + modules[mod][name].update(attrs) + return { "app_root": str(app_root), "files_importing_ocp": files,