diff --git a/microgen/__init__.py b/microgen/__init__.py index 1623109e..dac1514e 100644 --- a/microgen/__init__.py +++ b/microgen/__init__.py @@ -17,7 +17,6 @@ rescale, rotate, rotate_euler, - rotate_pv_euler, ) from .periodic import periodic_split_and_translate from .phase import Phase @@ -37,6 +36,7 @@ Ellipsoid, ExtrudedPolygon, FaceCenteredCubic, + GradedInfill, Infill, NormedDistance, Octahedron, @@ -48,6 +48,7 @@ Sphere, SphericalTpms, Spinodoid, + Sweep, Tpms, TruncatedCube, TruncatedCuboctahedron, @@ -76,6 +77,7 @@ "Ellipsoid", "ExtrudedPolygon", "FaceCenteredCubic", + "GradedInfill", "Infill", "NonBoxMeshError", "NormedDistance", @@ -91,6 +93,7 @@ "Shape", "Sphere", "Spinodoid", + "Sweep", "Tpms", "TruncatedCube", "TruncatedCuboctahedron", @@ -117,7 +120,6 @@ "rescale", "rotate", "rotate_euler", - "rotate_pv_euler", "Report", "SingleMesh", "surface_functions", diff --git a/microgen/cad.py b/microgen/cad.py deleted file mode 100644 index e6dd3668..00000000 --- a/microgen/cad.py +++ /dev/null @@ -1,1282 +0,0 @@ -"""CAD backend — direct OCCT (via OCP) replacement for CadQuery. - -========================================================= -CAD backend (:mod:`microgen.cad`) -========================================================= - -All CadQuery calls in microgen have been replaced by direct OCP -(``cadquery-ocp-novtk``) calls housed in this module. OCP is an *optional* -dependency — install via ``pip install microgen[cad]``. - -The module's top-level body does not import OCP, so ``import microgen.cad`` -always succeeds. The OCP-dependent functions import it lazily and raise a -helpful ``ImportError`` with install instructions if OCP is missing. - -Return type ------------ - -CAD-producing functions return a :class:`CadShape` — a thin wrapper around -an OCCT ``TopoDS_Shape`` exposing ``.wrapped`` for downstream OCP calls, plus -convenience methods (``translate``, ``rotate``, ``fuse``, ``cut``, -``export_stl``, ``export_step``, ``export_brep``). ``.wrapped`` matches the -attribute name CadQuery's ``Shape`` exposed, so most legacy call sites keep -working unchanged. -""" - -from __future__ import annotations - -from collections.abc import Iterable, Sequence -from typing import TYPE_CHECKING, Any - -import numpy as np -import numpy.typing as npt - -if TYPE_CHECKING: - from pathlib import Path - - from OCP.TopoDS import TopoDS_Shape - - -_INSTALL_HINT = ( - "microgen's CAD backend requires the OCP (OCCT) Python bindings. " - "Install with: pip install 'microgen[cad]' " - "(this pulls cadquery-ocp-novtk; on conda-forge use `ocp` instead)." -) - - -class _Centre(tuple): - """Tuple-like 3D point exposing ``.x``, ``.y``, ``.z`` and ``.to_tuple()``. - - Returned by :meth:`CadShape.center`. Mimics just enough of - ``cadquery.Vector`` to drop-in where existing code uses - ``shape.center().to_tuple()`` or ``shape.center().x``. - """ - - __slots__ = () - - def __new__(cls, x: float, y: float, z: float) -> _Centre: - """Create a 3-tuple ``(x, y, z)``.""" - return super().__new__(cls, (float(x), float(y), float(z))) - - @property - def x(self) -> float: - """X coordinate.""" - return self[0] - - @property - def y(self) -> float: - """Y coordinate.""" - return self[1] - - @property - def z(self) -> float: - """Z coordinate.""" - return self[2] - - def to_tuple(self) -> tuple[float, float, float]: - """Return ``(x, y, z)`` as a plain tuple.""" - return (self[0], self[1], self[2]) - - -class _BBox: - """Axis-aligned bounding box exposing ``xmin`` / ``xmax`` / …. - - Returned by :meth:`CadShape.bounding_box`. Also indexable as a 6-tuple - ``(xmin, ymin, zmin, xmax, ymax, zmax)`` matching OCCT's ``Bnd_Box.Get``. - """ - - __slots__ = ("xmax", "xmin", "ymax", "ymin", "zmax", "zmin") - - def __init__( - self, - xmin: float, - ymin: float, - zmin: float, - xmax: float, - ymax: float, - zmax: float, - ) -> None: - """Initialize from the 6 axis-aligned extents.""" - self.xmin = float(xmin) - self.ymin = float(ymin) - self.zmin = float(zmin) - self.xmax = float(xmax) - self.ymax = float(ymax) - self.zmax = float(zmax) - - @property - def diagonal_length(self) -> float: - """Length of the box's space diagonal.""" - dx = self.xmax - self.xmin - dy = self.ymax - self.ymin - dz = self.zmax - self.zmin - return float((dx * dx + dy * dy + dz * dz) ** 0.5) - - -def require_cad() -> None: - """Raise :class:`ImportError` if the CAD backend (OCP) is not importable.""" - try: - import OCP # noqa: F401 - except ImportError as err: - raise ImportError(_INSTALL_HINT) from err - - -def _run_boolean(op_cls: Any, a: CadShape, b: CadShape, label: str) -> TopoDS_Shape: - """Run an OCCT boolean op and raise on failure. - - Older OCP releases don't expose ``HasErrors()`` / ``IsDone()`` on the - ``BRepAlgoAPI_*`` classes; we probe via ``getattr`` and skip the check - when the API isn't available. - """ - op = op_cls(a.wrapped, b.wrapped) - has_errors = getattr(op, "HasErrors", None) - if callable(has_errors) and has_errors(): - err_msg = f"BRepAlgoAPI_{label} failed" - raise RuntimeError(err_msg) - return op.Shape() - - -def _topods_cast(name: str) -> Any: - """Return ``TopoDS.`` cast helper, tolerant of OCP version drift. - - Older OCP releases expose the static cast as ``TopoDS.Shell_s`` (pybind11 - ``_s`` convention); newer releases expose the unsuffixed ``TopoDS.Shell``. - Try the suffixed form first, fall back to unsuffixed. - """ - from OCP.TopoDS import TopoDS - - return getattr(TopoDS, f"{name}_s", None) or getattr(TopoDS, name) - - -# --------------------------------------------------------------------------- -# CadShape wrapper -# --------------------------------------------------------------------------- - - -class CadShape: - """Thin wrapper around an OCCT ``TopoDS_Shape``. - - Preserves the ``.wrapped`` attribute name used by CadQuery so downstream - OCP calls (``BRepAlgoAPI_Fuse(a.wrapped, b.wrapped)``) keep working. - - ``_mesh_volume`` (optional) is a trusted volume in the source mesh's - units, set by mesh-derived constructors (e.g. the TPMS periodic shell) - where OCCT's surface-integral volume is unreliable on invalid topology. - :meth:`volume` prefers it over the OCCT integral when present. - """ - - __slots__ = ("_mesh_volume", "wrapped") - - def __init__(self, shape: TopoDS_Shape) -> None: - """Wrap an OCCT ``TopoDS_Shape``.""" - self.wrapped = shape - self._mesh_volume: float | None = None - - # -- transforms -------------------------------------------------------- - - def translate(self, offset: Sequence[float]) -> CadShape: - """Return a translated copy.""" - from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform - from OCP.gp import gp_Trsf, gp_Vec - - trsf = gp_Trsf() - trsf.SetTranslation( - gp_Vec(float(offset[0]), float(offset[1]), float(offset[2])) - ) - transformed = BRepBuilderAPI_Transform(self.wrapped, trsf, True).Shape() - return CadShape(transformed) - - def rotate( - self, - center: Sequence[float], - axis: Sequence[float], - angle_degrees: float, - ) -> CadShape: - """Return a rotated copy (angle in degrees, axis is a unit vector).""" - from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform - from OCP.gp import gp_Ax1, gp_Dir, gp_Pnt, gp_Trsf - - trsf = gp_Trsf() - ax = gp_Ax1( - gp_Pnt(float(center[0]), float(center[1]), float(center[2])), - gp_Dir(float(axis[0]), float(axis[1]), float(axis[2])), - ) - trsf.SetRotation(ax, float(np.deg2rad(angle_degrees))) - transformed = BRepBuilderAPI_Transform(self.wrapped, trsf, True).Shape() - return CadShape(transformed) - - def copy(self) -> CadShape: - """Return an independent copy (deep topology copy).""" - from OCP.BRepBuilderAPI import BRepBuilderAPI_Copy - - return CadShape(BRepBuilderAPI_Copy(self.wrapped).Shape()) - - # -- boolean ops ------------------------------------------------------- - - def fuse(self, other: CadShape) -> CadShape: - """Boolean fusion: ``self ∪ other``.""" - from OCP.BRepAlgoAPI import BRepAlgoAPI_Fuse - - return CadShape(_run_boolean(BRepAlgoAPI_Fuse, self, other, "Fuse")) - - def cut(self, other: CadShape) -> CadShape: - """Boolean difference: ``self \\ other``.""" - from OCP.BRepAlgoAPI import BRepAlgoAPI_Cut - - return CadShape(_run_boolean(BRepAlgoAPI_Cut, self, other, "Cut")) - - def intersect(self, other: CadShape) -> CadShape: - """Boolean intersection: ``self ∩ other``.""" - from OCP.BRepAlgoAPI import BRepAlgoAPI_Common - - return CadShape(_run_boolean(BRepAlgoAPI_Common, self, other, "Common")) - - # -- topology queries -------------------------------------------------- - - def solids(self) -> list[CadShape]: - """Enumerate contained solids.""" - from OCP.TopAbs import TopAbs_SOLID - from OCP.TopExp import TopExp_Explorer - - cast_solid = _topods_cast("Solid") - out: list[CadShape] = [] - exp = TopExp_Explorer(self.wrapped, TopAbs_SOLID) - while exp.More(): - out.append(CadShape(cast_solid(exp.Current()))) - exp.Next() - return out - - def vertices(self) -> list[tuple[float, float, float]]: - """Enumerate the vertex coordinates of the shape. - - Callers use this to check that a generated mesh has any vertices at - all (``assert np.any(shape.vertices())``). - """ - from OCP.BRep import BRep_Tool - from OCP.TopAbs import TopAbs_VERTEX - from OCP.TopExp import TopExp_Explorer - - cast_vertex = _topods_cast("Vertex") - pnt = getattr(BRep_Tool, "Pnt_s", None) or BRep_Tool.Pnt - out: list[tuple[float, float, float]] = [] - exp = TopExp_Explorer(self.wrapped, TopAbs_VERTEX) - while exp.More(): - v = cast_vertex(exp.Current()) - p = pnt(v) - out.append((float(p.X()), float(p.Y()), float(p.Z()))) - exp.Next() - return out - - def faces(self) -> list[CadShape]: - """Enumerate the faces of the shape.""" - from OCP.TopAbs import TopAbs_FACE - from OCP.TopExp import TopExp_Explorer - - cast_face = _topods_cast("Face") - out: list[CadShape] = [] - exp = TopExp_Explorer(self.wrapped, TopAbs_FACE) - while exp.More(): - out.append(CadShape(cast_face(exp.Current()))) - exp.Next() - return out - - def is_closed(self) -> bool: - """Whether the shape is topologically closed. - - Solids are always closed; for shells/compounds we read OCCT's - per-shape ``Closed`` flag (set by ``BRep_Builder::IsClosed`` when the - shell was built from a watertight set of faces). - """ - from OCP.TopAbs import TopAbs_SOLID - - if self.wrapped.ShapeType() == TopAbs_SOLID: - return True - return bool(self.wrapped.Closed()) - - def volume(self) -> float: - """Return the (unsigned) volume of the shape. - - OCCT's ``BRepGProp::VolumeProperties`` returns a *signed* volume that - depends on face orientation; mesh-built shells from - :func:`mesh_to_shell_brep` can carry inverted orientation and yield a - negative value. We return ``abs(...)`` so volumes are non-negative. - - If a mesh-derived volume was stashed on ``_mesh_volume`` AND the OCCT - solid is not valid (BRepCheck_Analyzer flags self-intersection / - non-manifold edges, common on raw marching-cubes input), we trust the - mesh volume — the OCCT surface integral on an invalid topology is - meaningless. - """ - from OCP.BRepCheck import BRepCheck_Analyzer - from OCP.BRepGProp import BRepGProp - from OCP.GProp import GProp_GProps - - if ( - self._mesh_volume is not None - and not BRepCheck_Analyzer( - self.wrapped, - ).IsValid() - ): - return float(abs(self._mesh_volume)) - - props = GProp_GProps() - BRepGProp.VolumeProperties_s(self.wrapped, props) - return float(abs(props.Mass())) - - def center(self) -> _Centre: - """Return the volumetric center of mass. - - The result is a :class:`_Centre` — exposes ``.x``, ``.y``, ``.z``, - ``.to_tuple()``, and unpacks like a tuple. - """ - from OCP.BRepGProp import BRepGProp - from OCP.GProp import GProp_GProps - - props = GProp_GProps() - BRepGProp.VolumeProperties_s(self.wrapped, props) - c = props.CentreOfMass() - return _Centre(float(c.X()), float(c.Y()), float(c.Z())) - - def bounding_box(self) -> _BBox: - """Return the axis-aligned bounding box. - - The result exposes ``xmin`` / ``xmax`` / … attributes (see - :class:`_BBox`). - """ - from OCP.Bnd import Bnd_Box - from OCP.BRepBndLib import BRepBndLib - - box = Bnd_Box() - # AddOptimal uses exact geometric bounds (not cached triangulation). - BRepBndLib.AddOptimal_s(self.wrapped, box, True, True) - xmin, ymin, zmin, xmax, ymax, zmax = box.Get() - return _BBox(xmin, ymin, zmin, xmax, ymax, zmax) - - # -- exports ----------------------------------------------------------- - - def export_stl( - self, - path: str | Path, - linear_deflection: float = 0.01, - angular_deflection: float = 0.5, - *, - ascii_mode: bool = False, - ) -> None: - """Export to STL. Mesh is regenerated at the given deflection.""" - from OCP.BRepMesh import BRepMesh_IncrementalMesh - from OCP.StlAPI import StlAPI_Writer - - BRepMesh_IncrementalMesh( - self.wrapped, - float(linear_deflection), - False, - float(angular_deflection), - True, - ) - writer = StlAPI_Writer() - writer.ASCIIMode = bool(ascii_mode) - if not writer.Write(self.wrapped, str(path)): - err_msg = f"STL write failed for {path!r}" - raise RuntimeError(err_msg) - - def export_step(self, path: str | Path) -> None: - """Export to STEP (AP214).""" - from OCP.IFSelect import IFSelect_RetDone - from OCP.STEPControl import ( - STEPControl_AsIs, - STEPControl_Writer, - ) - - writer = STEPControl_Writer() - status = writer.Transfer(self.wrapped, STEPControl_AsIs) - if status != IFSelect_RetDone: - err_msg = f"STEP transfer failed with status {status!r}" - raise RuntimeError(err_msg) - status = writer.Write(str(path)) - if status != IFSelect_RetDone: - err_msg = f"STEP write failed with status {status!r}" - raise RuntimeError(err_msg) - - def export_brep(self, path: str | Path) -> None: - """Export to OCCT native BREP.""" - from OCP.BRepTools import BRepTools - - ok = BRepTools.Write_s(self.wrapped, str(path)) - if not ok: - err_msg = f"BREP write failed for {path!r}" - raise RuntimeError(err_msg) - - -def import_step(path: str | Path) -> CadShape: - """Import a STEP file and return the resulting :class:`CadShape`. - - Multi-root STEP files are merged into a single ``TopoDS_Compound``. - """ - require_cad() - from OCP.IFSelect import IFSelect_RetDone - from OCP.STEPControl import STEPControl_Reader - - reader = STEPControl_Reader() - status = reader.ReadFile(str(path)) - if status != IFSelect_RetDone: - err_msg = f"STEP read failed for {path!r} with status {status!r}" - raise RuntimeError(err_msg) - reader.TransferRoots() - return CadShape(reader.OneShape()) - - -# --------------------------------------------------------------------------- -# Mesh → Shell (SOTA path: one TopoDS_Face with attached Poly_Triangulation) -# --------------------------------------------------------------------------- - - -class ShellCreationError(RuntimeError): - """Raised when a mesh cannot be converted into an OCCT shell.""" - - -def mesh_to_shape( - points: npt.NDArray[np.float64], - triangles: npt.NDArray[np.int64], -) -> CadShape: - """Convert a triangle mesh to a :class:`CadShape` via ``Poly_Triangulation``. - - One ``TopoDS_Face`` carries the full triangulation (OCCT native tessellated - BREP representation). This is the SOTA fast path: O(N) in pure OCCT C++ - with no Python-per-triangle overhead, and exports cleanly to STEP AP242 - (tessellated) and STL. - - :param points: ``(N, 3)`` array of vertex coordinates - :param triangles: ``(M, 3)`` array of 0-indexed triangle vertex indices - :return: wrapped ``TopoDS_Shell`` containing one tessellated face - :raises ShellCreationError: if the triangulation cannot be built - """ - require_cad() - from OCP.BRep import BRep_Builder - from OCP.gp import gp_Pnt - from OCP.Poly import Poly_Triangle, Poly_Triangulation - from OCP.TopoDS import TopoDS_Face, TopoDS_Shell - - pts = np.asarray(points, dtype=np.float64) - tris = np.asarray(triangles, dtype=np.int64) - if pts.ndim != 2 or pts.shape[1] != 3: - err_msg = f"points must be (N, 3), got {pts.shape}" - raise ValueError(err_msg) - if tris.ndim != 2 or tris.shape[1] != 3: - err_msg = f"triangles must be (M, 3), got {tris.shape}" - raise ValueError(err_msg) - if tris.size == 0: - err_msg = "Cannot build a shell from an empty triangle list" - raise ShellCreationError(err_msg) - - nb_nodes = int(pts.shape[0]) - nb_tri = int(tris.shape[0]) - triangulation = Poly_Triangulation(nb_nodes, nb_tri, False) - for i in range(nb_nodes): - triangulation.SetNode( - i + 1, gp_Pnt(float(pts[i, 0]), float(pts[i, 1]), float(pts[i, 2])) - ) - for i in range(nb_tri): - a, b, c = int(tris[i, 0]), int(tris[i, 1]), int(tris[i, 2]) - triangulation.SetTriangle(i + 1, Poly_Triangle(a + 1, b + 1, c + 1)) - - builder = BRep_Builder() - face = TopoDS_Face() - try: - builder.MakeFace(face, triangulation) - except Exception as err: - err_msg = "OCCT refused the triangulation — check bounds and field." - raise ShellCreationError(err_msg) from err - - shell = TopoDS_Shell() - builder.MakeShell(shell) - builder.Add(shell, face) - return CadShape(shell) - - -def _triangle_components( - triangles: npt.NDArray[np.int64], -) -> npt.NDArray[np.int64]: - """Label connected components of a triangle mesh by edge adjacency.""" - from collections import defaultdict - - n_tri = int(triangles.shape[0]) - edges_sorted = np.sort( - np.vstack( - [triangles[:, [0, 1]], triangles[:, [1, 2]], triangles[:, [2, 0]]], - ), - axis=1, - ) - _keys, inv = np.unique(edges_sorted, axis=0, return_inverse=True) - owner = np.tile(np.arange(n_tri), 3) - - edge_to_tris: dict[int, list[int]] = defaultdict(list) - for global_i, key_i in enumerate(inv): - edge_to_tris[int(key_i)].append(int(owner[global_i])) - - adj: list[list[int]] = [[] for _ in range(n_tri)] - for tlist in edge_to_tris.values(): - if len(tlist) == 2: - adj[tlist[0]].append(tlist[1]) - adj[tlist[1]].append(tlist[0]) - - component = -np.ones(n_tri, dtype=np.int64) - n_comp = 0 - for start in range(n_tri): - if component[start] >= 0: - continue - component[start] = n_comp - stack = [start] - while stack: - t = stack.pop() - for nb in adj[t]: - if component[nb] < 0: - component[nb] = n_comp - stack.append(nb) - n_comp += 1 - return component - - -def _walk_boundary_loops( - comp_tris: npt.NDArray[np.int64], -) -> list[list[int]]: - """Extract directed closed loops from the boundary edges of a triangle group. - - A boundary edge appears in exactly one triangle of the group; the loops - are the closed chains of those edges in their original triangle direction. - """ - from collections import defaultdict - - ce = np.vstack( - [comp_tris[:, [0, 1]], comp_tris[:, [1, 2]], comp_tris[:, [2, 0]]], - ) - cs = np.sort(ce, axis=1) - _keys, inv, counts = np.unique( - cs, - axis=0, - return_inverse=True, - return_counts=True, - ) - bedges = ce[counts[inv] == 1] - if len(bedges) == 0: - return [] - - used = np.zeros(len(bedges), dtype=bool) - start_map: dict[int, list[int]] = defaultdict(list) - for i, (a, _b) in enumerate(bedges): - start_map[int(a)].append(i) - - loops: list[list[int]] = [] - for start_i in range(len(bedges)): - if used[start_i]: - continue - cur = start_i - used[cur] = True - loop = [int(bedges[cur, 0])] - while True: - b = int(bedges[cur, 1]) - loop.append(b) - if b == loop[0]: - break - nxt = next( - (cand for cand in start_map[b] if not used[cand]), - None, - ) - if nxt is None: - break - used[nxt] = True - cur = nxt - if len(loop) >= 4 and loop[-1] == loop[0]: - loops.append(loop) - return loops - - -def mesh_to_planar_face( # noqa: C901 - points: npt.NDArray[np.float64], - triangles: npt.NDArray[np.int64], - plane_origin: Sequence[float], - plane_normal: Sequence[float], -) -> CadShape: - """Build planar BREP face(s) whose wires trace the triangle-group boundary. - - Each connected component of the triangle group becomes a - :class:`TopoDS_Face` whose underlying surface is ``Geom_Plane`` and - whose outer/inner wires follow the boundary edges of the group on the - plane. Suitable both for STEP export (the BRep represents the right - region, not a bounding rectangle) and for gmsh ``setPeriodic`` (the - plane equation is recognised and the trimmed wires define the slave/ - master mesh region identically on opposite cell sides). - - All triangle vertices must lie on the plane within OCCT tolerance. - - :param points: ``(N, 3)`` array of vertex coordinates - :param triangles: ``(M, 3)`` array of 0-indexed triangle vertex indices - :param plane_origin: a point on the plane - :param plane_normal: the plane's outward unit normal - :return: wrapped face (single component) or shell (multiple components) - :raises ShellCreationError: if no usable boundary loop is found - """ - require_cad() - from OCP.BRep import BRep_Builder - from OCP.BRepBuilderAPI import ( - BRepBuilderAPI_MakeEdge, - BRepBuilderAPI_MakeFace, - BRepBuilderAPI_MakeWire, - ) - from OCP.gp import gp_Dir, gp_Pln, gp_Pnt - from OCP.TopoDS import TopoDS_Shell - - cast_wire = _topods_cast("Wire") - - pts = np.asarray(points, dtype=np.float64) - tris = np.asarray(triangles, dtype=np.int64) - if tris.size == 0: - err_msg = "Cannot build a planar face from an empty triangle list" - raise ShellCreationError(err_msg) - - plane = gp_Pln( - gp_Pnt(float(plane_origin[0]), float(plane_origin[1]), float(plane_origin[2])), - gp_Dir(float(plane_normal[0]), float(plane_normal[1]), float(plane_normal[2])), - ) - - n = np.asarray(plane_normal, dtype=np.float64) - n = n / (float(np.linalg.norm(n)) or 1.0) - helper = np.array([1.0, 0.0, 0.0]) if abs(n[0]) < 0.9 else np.array([0.0, 1.0, 0.0]) - u_ax = np.cross(n, helper) - u_ax /= float(np.linalg.norm(u_ax)) or 1.0 - v_ax = np.cross(n, u_ax) - o_arr = np.asarray(plane_origin, dtype=np.float64) - - component = _triangle_components(tris) - n_comp = int(component.max()) + 1 if component.size else 0 - - pnt_cache: dict[int, gp_Pnt] = {} - - def _occ_pnt(idx: int) -> gp_Pnt: - cached = pnt_cache.get(idx) - if cached is None: - v = pts[idx] - cached = gp_Pnt(float(v[0]), float(v[1]), float(v[2])) - pnt_cache[idx] = cached - return cached - - def _signed_area(loop: list[int]) -> float: - rel = pts[loop[:-1]] - o_arr - u = rel @ u_ax - v = rel @ v_ax - return 0.5 * float(np.sum(u * np.roll(v, -1) - np.roll(u, -1) * v)) - - def _build_wire(loop: list[int]): - wb = BRepBuilderAPI_MakeWire() - for k in range(len(loop) - 1): - edge_builder = BRepBuilderAPI_MakeEdge( - _occ_pnt(loop[k]), _occ_pnt(loop[k + 1]) - ) - if not edge_builder.IsDone(): - err_msg = "BRepBuilderAPI_MakeEdge failed for boundary segment" - raise ShellCreationError(err_msg) - wb.Add(edge_builder.Edge()) - if not wb.IsDone(): - err_msg = "BRepBuilderAPI_MakeWire failed for boundary loop" - raise ShellCreationError(err_msg) - return wb.Wire() - - def _build_component_face(loops: list[list[int]]): - areas = [_signed_area(L) for L in loops] - outer_local = int(np.argmax([abs(a) for a in areas])) - if areas[outer_local] < 0: - loops = [list(reversed(L)) for L in loops] - wires = [_build_wire(L) for L in loops] - face_builder = BRepBuilderAPI_MakeFace(plane, wires[outer_local]) - for k, w in enumerate(wires): - if k == outer_local: - continue - face_builder.Add(cast_wire(w.Reversed())) - if not face_builder.IsDone(): - err_msg = "BRepBuilderAPI_MakeFace failed building a planar face" - raise ShellCreationError(err_msg) - return face_builder.Face() - - component_faces = [] - for ci in range(n_comp): - comp_tris = tris[np.where(component == ci)[0]] - loops = _walk_boundary_loops(comp_tris) - if loops: - component_faces.append(_build_component_face(loops)) - - if not component_faces: - err_msg = "No planar face could be built from the triangle group" - raise ShellCreationError(err_msg) - - if len(component_faces) == 1: - return CadShape(component_faces[0]) - - builder = BRep_Builder() - shell = TopoDS_Shell() - builder.MakeShell(shell) - for f in component_faces: - builder.Add(shell, f) - return CadShape(shell) - - -def mesh_to_shell_brep( - points: npt.NDArray[np.float64], - triangles: npt.NDArray[np.int64], -) -> CadShape: - """Convert a triangle mesh to a shell with one planar BREP face per triangle. - - Slower than :func:`mesh_to_shape` (which uses a single tessellated face) - but produces a shell with *real* geometric surfaces — required whenever - the resulting shell is used as a cutting tool in boolean ops - (``BRepAlgoAPI_Cut``, ``Workplane.split()``, etc.), which refuse a - tessellated-only face. - - :param points: ``(N, 3)`` array of vertex coordinates - :param triangles: ``(M, 3)`` array of 0-indexed triangle vertex indices - :return: wrapped ``TopoDS_Shell`` with one planar face per triangle - :raises ShellCreationError: if any triangle cannot be built into a face - """ - require_cad() - from OCP.BRep import BRep_Builder - from OCP.BRepBuilderAPI import ( - BRepBuilderAPI_MakeEdge, - BRepBuilderAPI_MakeFace, - BRepBuilderAPI_MakeWire, - ) - from OCP.gp import gp_Pnt - from OCP.TopoDS import TopoDS_Shell - - pts = np.asarray(points, dtype=np.float64) - tris = np.asarray(triangles, dtype=np.int64) - if tris.size == 0: - err_msg = "Cannot build a shell from an empty triangle list" - raise ShellCreationError(err_msg) - - occ_points = [gp_Pnt(float(p[0]), float(p[1]), float(p[2])) for p in pts] - - builder = BRep_Builder() - shell = TopoDS_Shell() - builder.MakeShell(shell) - - try: - for a, b, c in tris: - e1 = BRepBuilderAPI_MakeEdge(occ_points[int(a)], occ_points[int(b)]).Edge() - e2 = BRepBuilderAPI_MakeEdge(occ_points[int(b)], occ_points[int(c)]).Edge() - e3 = BRepBuilderAPI_MakeEdge(occ_points[int(c)], occ_points[int(a)]).Edge() - wire = BRepBuilderAPI_MakeWire(e1, e2, e3).Wire() - face = BRepBuilderAPI_MakeFace(wire).Face() - builder.Add(shell, face) - except Exception as err: - err_msg = ( - "Failed to build the OCCT shell from the mesh; " - "try to increase the resolution or adjust bounds." - ) - raise ShellCreationError(err_msg) from err - - return CadShape(shell) - - -def mesh_to_sewn_shell( - points: npt.NDArray[np.float64], - triangles: npt.NDArray[np.int64], - tolerance: float | None = None, -) -> CadShape: - """Convert a triangle mesh to a sewn shell with shared edges. - - Builds one planar BREP face per triangle then sews them via - :class:`BRepBuilderAPI_Sewing` so coincident edges/vertices are merged - into shared topology. The resulting shell has valid topology, which - makes any subsequent boolean op (``Cut``, ``Common``, ``Fuse``) tractable - — unlike :func:`mesh_to_shell_brep` whose triangle-soup output is - pathologically slow for booleans. - - :param points: ``(N, 3)`` array of vertex coordinates - :param triangles: ``(M, 3)`` array of 0-indexed triangle vertex indices - :param tolerance: sewing tolerance; defaults to ``1e-6 * bbox_diag`` - :return: wrapped sewn ``TopoDS_Shape`` (Shell, or Compound of shells if - the input is disconnected) - :raises ShellCreationError: if any triangle cannot be built into a face - """ - require_cad() - from OCP.BRepBuilderAPI import ( - BRepBuilderAPI_MakeEdge, - BRepBuilderAPI_MakeFace, - BRepBuilderAPI_MakeWire, - BRepBuilderAPI_Sewing, - ) - from OCP.gp import gp_Pnt - - pts = np.asarray(points, dtype=np.float64) - tris = np.asarray(triangles, dtype=np.int64) - if tris.size == 0: - err_msg = "Cannot build a shell from an empty triangle list" - raise ShellCreationError(err_msg) - - if tolerance is None: - bbox_diag = float(np.linalg.norm(pts.max(axis=0) - pts.min(axis=0))) - tolerance = max(1e-9, 1e-6 * bbox_diag) - - occ_points = [gp_Pnt(float(p[0]), float(p[1]), float(p[2])) for p in pts] - - sewing = BRepBuilderAPI_Sewing(tolerance) - try: - for a, b, c in tris: - e1 = BRepBuilderAPI_MakeEdge(occ_points[int(a)], occ_points[int(b)]).Edge() - e2 = BRepBuilderAPI_MakeEdge(occ_points[int(b)], occ_points[int(c)]).Edge() - e3 = BRepBuilderAPI_MakeEdge(occ_points[int(c)], occ_points[int(a)]).Edge() - wire = BRepBuilderAPI_MakeWire(e1, e2, e3).Wire() - face = BRepBuilderAPI_MakeFace(wire).Face() - sewing.Add(face) - except Exception as err: - err_msg = ( - "Failed to build the OCCT shell from the mesh; " - "try to increase the resolution or adjust bounds." - ) - raise ShellCreationError(err_msg) from err - - sewing.Perform() - return CadShape(sewing.SewedShape()) - - -# --------------------------------------------------------------------------- -# Primitive builders -# --------------------------------------------------------------------------- - - -def make_box(dim: Sequence[float], center: Sequence[float]) -> CadShape: - """Axis-aligned box of size ``dim`` centered at ``center``.""" - require_cad() - from OCP.BRepPrimAPI import BRepPrimAPI_MakeBox - from OCP.gp import gp_Pnt - - dx, dy, dz = (float(d) for d in dim) - cx, cy, cz = (float(c) for c in center) - corner = gp_Pnt(cx - dx / 2.0, cy - dy / 2.0, cz - dz / 2.0) - return CadShape(BRepPrimAPI_MakeBox(corner, dx, dy, dz).Shape()) - - -def make_sphere(radius: float, center: Sequence[float]) -> CadShape: - """Sphere of given radius at ``center``.""" - require_cad() - from OCP.BRepPrimAPI import BRepPrimAPI_MakeSphere - from OCP.gp import gp_Pnt - - pnt = gp_Pnt(float(center[0]), float(center[1]), float(center[2])) - return CadShape(BRepPrimAPI_MakeSphere(pnt, float(radius)).Shape()) - - -def make_cylinder( - radius: float, - height: float, - center: Sequence[float], - axis: Sequence[float] = (1.0, 0.0, 0.0), -) -> CadShape: - """Cylinder of given radius and height centered at ``center`` along ``axis``.""" - require_cad() - from OCP.BRepPrimAPI import BRepPrimAPI_MakeCylinder - from OCP.gp import gp_Ax2, gp_Dir, gp_Pnt - - h = float(height) - ax_vec = np.asarray(axis, dtype=np.float64) - ax_vec = ax_vec / np.linalg.norm(ax_vec) - base = ( - float(center[0]) - h / 2.0 * ax_vec[0], - float(center[1]) - h / 2.0 * ax_vec[1], - float(center[2]) - h / 2.0 * ax_vec[2], - ) - ax = gp_Ax2( - gp_Pnt(*base), - gp_Dir(float(ax_vec[0]), float(ax_vec[1]), float(ax_vec[2])), - ) - return CadShape(BRepPrimAPI_MakeCylinder(ax, float(radius), h).Shape()) - - -def make_capsule( - radius: float, - height: float, - center: Sequence[float], -) -> CadShape: - """Capsule (cylinder along X with hemispherical caps).""" - require_cad() - from OCP.BRepPrimAPI import ( - BRepPrimAPI_MakeCylinder, - BRepPrimAPI_MakeSphere, - ) - from OCP.gp import gp_Ax2, gp_Dir, gp_Pnt - - cx, cy, cz = (float(c) for c in center) - h = float(height) - r = float(radius) - base_axis = gp_Ax2(gp_Pnt(cx - h / 2.0, cy, cz), gp_Dir(1.0, 0.0, 0.0)) - cyl = BRepPrimAPI_MakeCylinder(base_axis, r, h).Shape() - left = BRepPrimAPI_MakeSphere(gp_Pnt(cx - h / 2.0, cy, cz), r).Shape() - right = BRepPrimAPI_MakeSphere(gp_Pnt(cx + h / 2.0, cy, cz), r).Shape() - fused = CadShape(cyl).fuse(CadShape(left)).fuse(CadShape(right)) - return fused - - -def make_ellipsoid(radii: Sequence[float], center: Sequence[float]) -> CadShape: - """Ellipsoid of the given axis-aligned radii at ``center``. - - Built as a unit sphere transformed by a non-uniform scaling matrix. - """ - require_cad() - from OCP.BRepBuilderAPI import BRepBuilderAPI_GTransform - from OCP.BRepPrimAPI import BRepPrimAPI_MakeSphere - from OCP.gp import gp_GTrsf, gp_Mat, gp_Pnt, gp_XYZ - - rx, ry, rz = (float(r) for r in radii) - cx, cy, cz = (float(c) for c in center) - sphere = BRepPrimAPI_MakeSphere(gp_Pnt(0.0, 0.0, 0.0), 1.0).Shape() - - gtrsf = gp_GTrsf() - gtrsf.SetVectorialPart( - gp_Mat( - rx, - 0.0, - 0.0, - 0.0, - ry, - 0.0, - 0.0, - 0.0, - rz, - ), - ) - gtrsf.SetTranslationPart(gp_XYZ(cx, cy, cz)) - return CadShape(BRepBuilderAPI_GTransform(sphere, gtrsf, True).Shape()) - - -def make_polyhedron( - vertices: Sequence[Sequence[float]], - faces_ixs: Sequence[Sequence[int]], - center: Sequence[float] = (0.0, 0.0, 0.0), -) -> CadShape: - """Polyhedron from vertex list and face→vertex index list. - - Each face's vertex index list must be closed (last == first). - - Uses ``BRepBuilderAPI_Sewing`` to stitch the independently-constructed - faces into a shared-edge shell with consistent outward orientation, - then ``BRepBuilderAPI_MakeSolid`` to close it. (Raw ``BRep_Builder.Add`` - does not orient; the resulting solid would have mixed-sign volume.) - """ - require_cad() - from OCP.BRepBuilderAPI import ( - BRepBuilderAPI_MakeEdge, - BRepBuilderAPI_MakeFace, - BRepBuilderAPI_MakeSolid, - BRepBuilderAPI_MakeWire, - BRepBuilderAPI_Sewing, - ) - from OCP.gp import gp_Pnt - from OCP.ShapeFix import ShapeFix_Solid - from OCP.TopAbs import TopAbs_SHELL - from OCP.TopExp import TopExp_Explorer - - cx, cy, cz = (float(c) for c in center) - points = [ - gp_Pnt(float(v[0]) + cx, float(v[1]) + cy, float(v[2]) + cz) for v in vertices - ] - - sewing = BRepBuilderAPI_Sewing() - for ixs in faces_ixs: - wire_builder = BRepBuilderAPI_MakeWire() - for i1, i2 in zip(ixs, ixs[1:], strict=False): - edge_builder = BRepBuilderAPI_MakeEdge(points[i1], points[i2]) - if not edge_builder.IsDone(): - err_msg = "BRepBuilderAPI_MakeEdge failed for polyhedron face" - raise ShellCreationError(err_msg) - wire_builder.Add(edge_builder.Edge()) - if not wire_builder.IsDone(): - err_msg = "BRepBuilderAPI_MakeWire failed for polyhedron face" - raise ShellCreationError(err_msg) - face_builder = BRepBuilderAPI_MakeFace(wire_builder.Wire()) - if not face_builder.IsDone(): - err_msg = "BRepBuilderAPI_MakeFace failed for polyhedron face" - raise ShellCreationError(err_msg) - sewing.Add(face_builder.Face()) - sewing.Perform() - sewn = sewing.SewedShape() - - # Extract the shell (sewing may return it directly or wrapped in a compound). - exp = TopExp_Explorer(sewn, TopAbs_SHELL) - if not exp.More(): - err_msg = "Sewing did not produce a shell — check face connectivity" - raise ShellCreationError(err_msg) - shell = _topods_cast("Shell")(exp.Current()) - - solid_builder = BRepBuilderAPI_MakeSolid(shell) - if not solid_builder.IsDone(): - err_msg = "BRepBuilderAPI_MakeSolid failed; sewn shell is not closed" - raise ShellCreationError(err_msg) - fixer = ShapeFix_Solid(solid_builder.Solid()) - fixer.Perform() - return CadShape(fixer.Solid()) - - -def make_extruded_polygon( - list_corners: Sequence[tuple[float, float]], - height: float, - center: Sequence[float], -) -> CadShape: - """Extrude a 2D polygon (in the YZ plane) along the X axis.""" - require_cad() - from OCP.BRepBuilderAPI import ( - BRepBuilderAPI_MakeEdge, - BRepBuilderAPI_MakeFace, - BRepBuilderAPI_MakeWire, - ) - from OCP.BRepPrimAPI import BRepPrimAPI_MakePrism - from OCP.gp import gp_Pnt, gp_Vec - - cx, cy, cz = (float(c) for c in center) - h = float(height) - x_base = cx - h / 2.0 - pts = [gp_Pnt(x_base, cy + float(y), cz + float(z)) for (y, z) in list_corners] - if (list_corners[0][0], list_corners[0][1]) != ( - list_corners[-1][0], - list_corners[-1][1], - ): - pts.append(pts[0]) - - wire_builder = BRepBuilderAPI_MakeWire() - for i in range(len(pts) - 1): - edge_builder = BRepBuilderAPI_MakeEdge(pts[i], pts[i + 1]) - if not edge_builder.IsDone(): - err_msg = "BRepBuilderAPI_MakeEdge failed for extruded polygon" - raise ShellCreationError(err_msg) - wire_builder.Add(edge_builder.Edge()) - if not wire_builder.IsDone(): - err_msg = "BRepBuilderAPI_MakeWire failed for extruded polygon" - raise ShellCreationError(err_msg) - face_builder = BRepBuilderAPI_MakeFace(wire_builder.Wire()) - if not face_builder.IsDone(): - err_msg = "BRepBuilderAPI_MakeFace failed for extruded polygon" - raise ShellCreationError(err_msg) - extruded = BRepPrimAPI_MakePrism(face_builder.Face(), gp_Vec(h, 0.0, 0.0)).Shape() - return CadShape(extruded) - - -# --------------------------------------------------------------------------- -# Compound assembly (for lattice) -# --------------------------------------------------------------------------- - - -def make_compound(shapes: Iterable[CadShape]) -> CadShape: - """Assemble shapes into a single OCCT ``TopoDS_Compound``.""" - require_cad() - from OCP.BRep import BRep_Builder - from OCP.TopoDS import TopoDS_Compound - - builder = BRep_Builder() - compound = TopoDS_Compound() - builder.MakeCompound(compound) - for s in shapes: - builder.Add(compound, s.wrapped) - return CadShape(compound) - - -def make_compound_from_solids(solids: Iterable[Any]) -> CadShape: - """Assemble raw OCCT ``TopoDS_Shape`` solids (not ``CadShape``) into a compound.""" - require_cad() - from OCP.BRep import BRep_Builder - from OCP.TopoDS import TopoDS_Compound - - builder = BRep_Builder() - compound = TopoDS_Compound() - builder.MakeCompound(compound) - for s in solids: - shape = s.wrapped if hasattr(s, "wrapped") else s - builder.Add(compound, shape) - return CadShape(compound) - - -# --------------------------------------------------------------------------- -# Topology exploration and splitting -# --------------------------------------------------------------------------- - - -def enumerate_solids(shape: CadShape) -> list[Any]: - """Return the list of ``TopoDS_Solid`` inside a shape (empty if none).""" - require_cad() - from OCP.TopAbs import TopAbs_SOLID - from OCP.TopExp import TopExp_Explorer - - cast_solid = _topods_cast("Solid") - out: list[Any] = [] - exp = TopExp_Explorer(shape.wrapped, TopAbs_SOLID) - while exp.More(): - out.append(cast_solid(exp.Current())) - exp.Next() - return out - - -def split_shape( - shape: CadShape, - tool: CadShape | Iterable[CadShape], - *, - fuzzy_value: float = 1e-4, -) -> CadShape: - """Split *shape* by *tool* using OCCT's ``BRepAlgoAPI_Splitter``. - - The result is a :class:`CadShape` wrapping a ``TopoDS_Compound`` that - contains the sub-shapes produced by the split. Use - :func:`enumerate_solids` to iterate over the resulting solids. - - :param tool: a single :class:`CadShape` or an iterable of them. Pass - multiple tools when each tool is a separate ``TopoDS_Shell`` and - you want to keep them topologically distinct — the OCCT splitter - treats each list entry as one cutting tool. This matters because - ``BRepAlgoAPI_Fuse`` on two shells *decomposes* them into a - compound of per-face shells, which then defeats the splitter. - :param fuzzy_value: tolerance forwarded to ``SetFuzzyValue`` so OCCT - recognises near-coincident geometry as touching. Necessary when - a tool is a tessellated shell whose boundary edges lie on - ``shape``'s planar faces only up to a few microns of drift. - """ - require_cad() - from OCP.BRepAlgoAPI import BRepAlgoAPI_Splitter - from OCP.TopTools import TopTools_ListOfShape - - args = TopTools_ListOfShape() - args.Append(shape.wrapped) - tools = TopTools_ListOfShape() - if isinstance(tool, CadShape): - tools.Append(tool.wrapped) - else: - for t in tool: - tools.Append(t.wrapped) - - splitter = BRepAlgoAPI_Splitter() - splitter.SetArguments(args) - splitter.SetTools(tools) - if fuzzy_value > 0: - splitter.SetFuzzyValue(float(fuzzy_value)) - splitter.Build() - return CadShape(splitter.Shape()) - - -def make_plane_face( - base_pnt: Sequence[float], - direction: Sequence[float], - half_size: float = 1.0e6, -) -> CadShape: - """Build a large planar face used as a cutting tool. - - :param base_pnt: a point on the plane - :param direction: plane normal - :param half_size: half-edge of the square face (default large so the plane - reaches far outside any realistic shape) - """ - require_cad() - from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeFace - from OCP.gp import gp_Ax3, gp_Dir, gp_Pln, gp_Pnt - - pnt = gp_Pnt(float(base_pnt[0]), float(base_pnt[1]), float(base_pnt[2])) - nrm = gp_Dir(float(direction[0]), float(direction[1]), float(direction[2])) - plane = gp_Pln(gp_Ax3(pnt, nrm)) - face = BRepBuilderAPI_MakeFace( - plane, - -float(half_size), - float(half_size), - -float(half_size), - float(half_size), - ).Face() - return CadShape(face) - - -def transform_geometry(shape: CadShape, matrix: npt.NDArray[np.float64]) -> CadShape: - """Apply a 3x4 affine matrix (linear + translation). - - Wraps OCCT ``BRepBuilderAPI_GTransform``. - - :param matrix: ``(3, 4)`` array; rows are ``[a b c tx; d e f ty; g h i tz]``. - """ - require_cad() - from OCP.BRepBuilderAPI import BRepBuilderAPI_GTransform - from OCP.gp import gp_GTrsf, gp_Mat, gp_XYZ - - m = np.asarray(matrix, dtype=np.float64) - gtrsf = gp_GTrsf() - gtrsf.SetVectorialPart( - gp_Mat( - float(m[0, 0]), - float(m[0, 1]), - float(m[0, 2]), - float(m[1, 0]), - float(m[1, 1]), - float(m[1, 2]), - float(m[2, 0]), - float(m[2, 1]), - float(m[2, 2]), - ), - ) - gtrsf.SetTranslationPart(gp_XYZ(float(m[0, 3]), float(m[1, 3]), float(m[2, 3]))) - return CadShape(BRepBuilderAPI_GTransform(shape.wrapped, gtrsf, True).Shape()) - - -def translate_solid(solid: Any, offset: Sequence[float]) -> Any: - """Translate a raw OCCT solid/shape by ``offset`` and return the same type.""" - require_cad() - from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform - from OCP.gp import gp_Trsf, gp_Vec - - shape = solid.wrapped if hasattr(solid, "wrapped") else solid - trsf = gp_Trsf() - trsf.SetTranslation( - gp_Vec(float(offset[0]), float(offset[1]), float(offset[2])), - ) - return BRepBuilderAPI_Transform(shape, trsf, True).Shape() - - -def solid_center(shape: Any) -> tuple[float, float, float]: - """Return the volumetric center of mass of a raw OCCT shape/solid.""" - require_cad() - from OCP.BRepGProp import BRepGProp - from OCP.GProp import GProp_GProps - - s = shape.wrapped if hasattr(shape, "wrapped") else shape - props = GProp_GProps() - BRepGProp.VolumeProperties_s(s, props) - com = props.CentreOfMass() - return (float(com.X()), float(com.Y()), float(com.Z())) - - -def select_solids_on_side( - shape: CadShape, - base_pnt: Sequence[float], - side_direction: Sequence[float], -) -> list[Any]: - """Enumerate a compound's solids, keep those whose centroid is on the - positive side of the plane through ``base_pnt`` normal to ``side_direction``. - - Matches CadQuery's ``.solids(">X")`` / ``.solids(" 0: - selected.append(solid) - return selected - - -def intersect_solids_with_box(solids: Iterable[Any], box: CadShape) -> CadShape: - """Intersect each solid with *box* and fuse the results into a ``CadShape``. - - Returns a :class:`CadShape` wrapping a compound (possibly empty). - """ - require_cad() - from OCP.BRepAlgoAPI import BRepAlgoAPI_Common - - parts: list[Any] = [] - for solid in solids: - s = solid.wrapped if hasattr(solid, "wrapped") else solid - common = BRepAlgoAPI_Common(s, box.wrapped).Shape() - # Keep if it contains at least one solid. - if enumerate_solids(CadShape(common)): - parts.append(common) - return make_compound_from_solids(parts) diff --git a/microgen/cad/__init__.py b/microgen/cad/__init__.py new file mode 100644 index 00000000..e532525c --- /dev/null +++ b/microgen/cad/__init__.py @@ -0,0 +1,124 @@ +"""CAD backend — direct OCCT (via OCP) replacement for CadQuery. + +========================================================= +CAD backend (:mod:`microgen.cad`) +========================================================= + +All CadQuery calls in microgen have been replaced by direct OCP +(``cadquery-ocp-novtk``) calls housed in this subpackage. OCP is an +*optional* dependency — install via ``pip install microgen[cad]``. + +The subpackage's submodules don't import OCP at top level, so +``import microgen.cad`` always succeeds. OCP-dependent functions import +it lazily and raise a helpful ``ImportError`` with install instructions if +OCP is missing. + +Return type +----------- + +CAD-producing functions return a :class:`CadShape` — a thin wrapper around +an OCCT ``TopoDS_Shape`` exposing ``.wrapped`` for downstream OCP calls, plus +convenience methods (``translate``, ``rotate``, ``fuse``, ``cut``, +``export_stl``, ``export_step``, ``export_brep``). ``.wrapped`` matches the +attribute name CadQuery's ``Shape`` exposed, so most legacy call sites keep +working unchanged. + +Layout +------ + +- ``cad._install`` — ``require_cad`` / ``_INSTALL_HINT`` gate +- ``cad.shape`` — ``CadShape`` class + ``_Centre`` / ``_BBox`` / + ``_run_boolean`` / ``_topods_cast`` / ``ShellCreationError`` +- ``cad.io`` — ``import_step`` (export methods live on ``CadShape``) +- ``cad.meshbridge`` — ``mesh_to_shape`` / ``shape_to_cad`` / + ``mesh_to_planar_face`` / ``mesh_to_shell_brep`` / ``mesh_to_sewn_shell`` +- ``cad.primitives`` — ``make_box`` / ``make_sphere`` / ``make_cylinder`` / + ``make_capsule`` / ``make_ellipsoid`` / ``make_polyhedron`` / + ``make_extruded_polygon`` +- ``cad.topo`` — ``make_compound`` / ``make_compound_from_solids`` / + ``enumerate_solids`` / ``split_shape`` / ``make_plane_face`` / + ``transform_geometry`` / ``translate_solid`` / ``solid_center`` / + ``select_solids_on_side`` / ``intersect_solids_with_box`` + +All public symbols are re-exported here so existing +``from microgen.cad import …`` imports keep working unchanged. +""" + +from __future__ import annotations + +from ._install import _INSTALL_HINT, require_cad +from .io import import_step +from .meshbridge import ( + _triangle_components, + _walk_boundary_loops, + mesh_to_planar_face, + mesh_to_sewn_shell, + mesh_to_shape, + mesh_to_shell_brep, + shape_to_cad, +) +from .primitives import ( + make_box, + make_capsule, + make_cylinder, + make_ellipsoid, + make_extruded_polygon, + make_polyhedron, + make_sphere, +) +from .shape import ( + CadShape, + ShellCreationError, + _BBox, + _Centre, + _run_boolean, + _topods_cast, +) +from .topo import ( + enumerate_solids, + intersect_solids_with_box, + make_compound, + make_compound_from_solids, + make_plane_face, + select_solids_on_side, + solid_center, + split_shape, + transform_geometry, + translate_solid, +) + +__all__ = [ + "CadShape", + "ShellCreationError", + "_BBox", + "_Centre", + "_INSTALL_HINT", + "_run_boolean", + "_topods_cast", + "_triangle_components", + "_walk_boundary_loops", + "enumerate_solids", + "import_step", + "intersect_solids_with_box", + "make_box", + "make_capsule", + "make_compound", + "make_compound_from_solids", + "make_cylinder", + "make_ellipsoid", + "make_extruded_polygon", + "make_plane_face", + "make_polyhedron", + "make_sphere", + "mesh_to_planar_face", + "mesh_to_sewn_shell", + "mesh_to_shape", + "mesh_to_shell_brep", + "require_cad", + "select_solids_on_side", + "shape_to_cad", + "solid_center", + "split_shape", + "transform_geometry", + "translate_solid", +] diff --git a/microgen/cad/_install.py b/microgen/cad/_install.py new file mode 100644 index 00000000..a46f54b6 --- /dev/null +++ b/microgen/cad/_install.py @@ -0,0 +1,22 @@ +"""CAD backend install gate. + +``require_cad()`` is the canonical way for lazy CAD code paths to verify +that the optional ``[cad]`` extra is available before touching OCP, and to +raise a user-friendly ``ImportError`` with install instructions otherwise. +""" + +from __future__ import annotations + +_INSTALL_HINT = ( + "microgen's CAD backend requires the OCP (OCCT) Python bindings. " + "Install with: pip install 'microgen[cad]' " + "(this pulls cadquery-ocp-novtk; on conda-forge use `ocp` instead)." +) + + +def require_cad() -> None: + """Raise :class:`ImportError` if the CAD backend (OCP) is not importable.""" + try: + import OCP # noqa: F401, PLC0415 + except ImportError as err: + raise ImportError(_INSTALL_HINT) from err diff --git a/microgen/cad/io.py b/microgen/cad/io.py new file mode 100644 index 00000000..22599c8b --- /dev/null +++ b/microgen/cad/io.py @@ -0,0 +1,34 @@ +"""STEP file I/O (read). + +STEP/BREP/STL *export* is exposed as methods on :class:`CadShape` itself +(``export_step``, ``export_brep``, ``export_stl``) — see ``cad/shape.py``. +This module holds the standalone ``import_step`` reader. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from ._install import require_cad +from .shape import CadShape + +if TYPE_CHECKING: + from pathlib import Path + + +def import_step(path: str | Path) -> CadShape: + """Import a STEP file and return the resulting :class:`CadShape`. + + Multi-root STEP files are merged into a single ``TopoDS_Compound``. + """ + require_cad() + from OCP.IFSelect import IFSelect_RetDone # noqa: PLC0415 + from OCP.STEPControl import STEPControl_Reader # noqa: PLC0415 + + reader = STEPControl_Reader() + status = reader.ReadFile(str(path)) + if status != IFSelect_RetDone: + err_msg = f"STEP read failed for {path!r} with status {status!r}" + raise RuntimeError(err_msg) + reader.TransferRoots() + return CadShape(reader.OneShape()) diff --git a/microgen/cad/meshbridge.py b/microgen/cad/meshbridge.py new file mode 100644 index 00000000..ceab2d43 --- /dev/null +++ b/microgen/cad/meshbridge.py @@ -0,0 +1,478 @@ +"""Mesh ↔ BREP bridges. + +Functions in this module convert triangle meshes (numpy point/index arrays, +or pyvista ``PolyData``) into OCCT ``CadShape`` representations, or vice +versa (via the implicit ``Shape`` → BREP bridge in :func:`shape_to_cad`). +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +import numpy.typing as npt + +from ._install import require_cad +from .shape import CadShape, ShellCreationError, _topods_cast + +if TYPE_CHECKING: + from collections.abc import Sequence + + from ..shape._types import BoundsType + from ..shape.shape import Shape + + +def mesh_to_shape( + points: npt.NDArray[np.float64], + triangles: npt.NDArray[np.int64], +) -> CadShape: + """Convert a triangle mesh to a :class:`CadShape` via ``Poly_Triangulation``. + + One ``TopoDS_Face`` carries the full triangulation (OCCT native tessellated + BREP representation). This is the SOTA fast path: O(N) in pure OCCT C++ + with no Python-per-triangle overhead, and exports cleanly to STEP AP242 + (tessellated) and STL. + + :param points: ``(N, 3)`` array of vertex coordinates + :param triangles: ``(M, 3)`` array of 0-indexed triangle vertex indices + :return: wrapped ``TopoDS_Shell`` containing one tessellated face + :raises ShellCreationError: if the triangulation cannot be built + """ + require_cad() + from OCP.BRep import BRep_Builder # noqa: PLC0415 + from OCP.gp import gp_Pnt # noqa: PLC0415 + from OCP.Poly import Poly_Triangle, Poly_Triangulation # noqa: PLC0415 + from OCP.TopoDS import TopoDS_Face, TopoDS_Shell # noqa: PLC0415 + + pts = np.asarray(points, dtype=np.float64) + tris = np.asarray(triangles, dtype=np.int64) + if pts.ndim != 2 or pts.shape[1] != 3: + err_msg = f"points must be (N, 3), got {pts.shape}" + raise ValueError(err_msg) + if tris.ndim != 2 or tris.shape[1] != 3: + err_msg = f"triangles must be (M, 3), got {tris.shape}" + raise ValueError(err_msg) + if tris.size == 0: + err_msg = "Cannot build a shell from an empty triangle list" + raise ShellCreationError(err_msg) + + nb_nodes = int(pts.shape[0]) + nb_tri = int(tris.shape[0]) + triangulation = Poly_Triangulation(nb_nodes, nb_tri, False) + for i in range(nb_nodes): + triangulation.SetNode( + i + 1, gp_Pnt(float(pts[i, 0]), float(pts[i, 1]), float(pts[i, 2])) + ) + for i in range(nb_tri): + a, b, c = int(tris[i, 0]), int(tris[i, 1]), int(tris[i, 2]) + triangulation.SetTriangle(i + 1, Poly_Triangle(a + 1, b + 1, c + 1)) + + builder = BRep_Builder() + face = TopoDS_Face() + try: + builder.MakeFace(face, triangulation) + except Exception as err: + err_msg = "OCCT refused the triangulation — check bounds and field." + raise ShellCreationError(err_msg) from err + + shell = TopoDS_Shell() + builder.MakeShell(shell) + builder.Add(shell, face) + return CadShape(shell) + + +def shape_to_cad( + shape: Shape, + bounds: BoundsType | None = None, + resolution: int = 50, +) -> CadShape: + """Bridge an implicit :class:`~microgen.shape.shape.Shape` to a CAD BREP. + + Runs the shape's :meth:`generate_surface_mesh` (marching cubes on the + SDF for free Shapes; native renderer for concrete subclasses) then + wraps the resulting triangle mesh into a single tessellated BREP face + via :func:`mesh_to_shape`. + + Concrete subclasses with native primitive paths (``Box``, ``Sphere``, + ``Tpms``, ``Spinodoid`` …) usually skip this bridge and call their + own ``generate_cad``. This function is the generic fallback for + bare ``Shape`` instances built via :func:`microgen.shape.implicit_ops.from_field` + or boolean composition (``a | b``, ``a - b``, ...). + + Requires the optional ``[cad]`` install extra. + + :param shape: the implicit shape to materialise + :param bounds: ``(xmin, xmax, ymin, ymax, zmin, zmax)``; defaults to + ``shape.bounds`` if set, else raises ``ValueError`` + :param resolution: marching-cubes grid resolution per axis + :return: :class:`CadShape` wrapping the tessellated ``TopoDS_Shell`` + """ + require_cad() + if shape.func is None: + err_msg = "No implicit field defined — cannot build BREP from an empty Shape" + raise NotImplementedError(err_msg) + + mesh = shape.generate_surface_mesh(bounds=bounds, resolution=resolution) + if mesh.n_cells == 0: + err_msg = "Generated mesh is empty — check bounds and field function" + raise ValueError(err_msg) + + if not mesh.is_all_triangles: + mesh.triangulate(inplace=True) + triangles = mesh.faces.reshape(-1, 4)[:, 1:] + points = np.asarray(mesh.points, dtype=np.float64) + + from ..shape.shape import ShellCreationError as _ShapeShellError # noqa: PLC0415 + + try: + return mesh_to_shape(points, triangles) + except Exception as err: + err_msg = ( + "Failed to build the OCCT shell from the mesh; " + "try to increase the resolution or adjust bounds." + ) + raise _ShapeShellError(err_msg) from err + + +def _triangle_components( + triangles: npt.NDArray[np.int64], +) -> npt.NDArray[np.int64]: + """Label connected components of a triangle mesh by edge adjacency.""" + from collections import defaultdict # noqa: PLC0415 + + n_tri = int(triangles.shape[0]) + edges_sorted = np.sort( + np.vstack( + [triangles[:, [0, 1]], triangles[:, [1, 2]], triangles[:, [2, 0]]], + ), + axis=1, + ) + _keys, inv = np.unique(edges_sorted, axis=0, return_inverse=True) + owner = np.tile(np.arange(n_tri), 3) + + edge_to_tris: dict[int, list[int]] = defaultdict(list) + for global_i, key_i in enumerate(inv): + edge_to_tris[int(key_i)].append(int(owner[global_i])) + + adj: list[list[int]] = [[] for _ in range(n_tri)] + for tlist in edge_to_tris.values(): + if len(tlist) == 2: + adj[tlist[0]].append(tlist[1]) + adj[tlist[1]].append(tlist[0]) + + component = -np.ones(n_tri, dtype=np.int64) + n_comp = 0 + for start in range(n_tri): + if component[start] >= 0: + continue + component[start] = n_comp + stack = [start] + while stack: + t = stack.pop() + for nb in adj[t]: + if component[nb] < 0: + component[nb] = n_comp + stack.append(nb) + n_comp += 1 + return component + + +def _walk_boundary_loops( + comp_tris: npt.NDArray[np.int64], +) -> list[list[int]]: + """Extract directed closed loops from the boundary edges of a triangle group. + + A boundary edge appears in exactly one triangle of the group; the loops + are the closed chains of those edges in their original triangle direction. + """ + from collections import defaultdict # noqa: PLC0415 + + ce = np.vstack( + [comp_tris[:, [0, 1]], comp_tris[:, [1, 2]], comp_tris[:, [2, 0]]], + ) + cs = np.sort(ce, axis=1) + _keys, inv, counts = np.unique( + cs, + axis=0, + return_inverse=True, + return_counts=True, + ) + bedges = ce[counts[inv] == 1] + if len(bedges) == 0: + return [] + + used = np.zeros(len(bedges), dtype=bool) + start_map: dict[int, list[int]] = defaultdict(list) + for i, (a, _b) in enumerate(bedges): + start_map[int(a)].append(i) + + loops: list[list[int]] = [] + for start_i in range(len(bedges)): + if used[start_i]: + continue + cur = start_i + used[cur] = True + loop = [int(bedges[cur, 0])] + while True: + b = int(bedges[cur, 1]) + loop.append(b) + if b == loop[0]: + break + nxt = next( + (cand for cand in start_map[b] if not used[cand]), + None, + ) + if nxt is None: + break + used[nxt] = True + cur = nxt + if len(loop) >= 4 and loop[-1] == loop[0]: + loops.append(loop) + return loops + + +def mesh_to_planar_face( # noqa: C901 + points: npt.NDArray[np.float64], + triangles: npt.NDArray[np.int64], + plane_origin: Sequence[float], + plane_normal: Sequence[float], +) -> CadShape: + """Build planar BREP face(s) whose wires trace the triangle-group boundary. + + Each connected component of the triangle group becomes a + :class:`TopoDS_Face` whose underlying surface is ``Geom_Plane`` and + whose outer/inner wires follow the boundary edges of the group on the + plane. Suitable both for STEP export (the BRep represents the right + region, not a bounding rectangle) and for gmsh ``setPeriodic`` (the + plane equation is recognised and the trimmed wires define the slave/ + master mesh region identically on opposite cell sides). + + All triangle vertices must lie on the plane within OCCT tolerance. + + :param points: ``(N, 3)`` array of vertex coordinates + :param triangles: ``(M, 3)`` array of 0-indexed triangle vertex indices + :param plane_origin: a point on the plane + :param plane_normal: the plane's outward unit normal + :return: wrapped face (single component) or shell (multiple components) + :raises ShellCreationError: if no usable boundary loop is found + """ + require_cad() + from OCP.BRep import BRep_Builder # noqa: PLC0415 + from OCP.BRepBuilderAPI import ( # noqa: PLC0415 + BRepBuilderAPI_MakeEdge, + BRepBuilderAPI_MakeFace, + BRepBuilderAPI_MakeWire, + ) + from OCP.gp import gp_Dir, gp_Pln, gp_Pnt # noqa: PLC0415 + from OCP.TopoDS import TopoDS_Shell # noqa: PLC0415 + + cast_wire = _topods_cast("Wire") + + pts = np.asarray(points, dtype=np.float64) + tris = np.asarray(triangles, dtype=np.int64) + if tris.size == 0: + err_msg = "Cannot build a planar face from an empty triangle list" + raise ShellCreationError(err_msg) + + plane = gp_Pln( + gp_Pnt(float(plane_origin[0]), float(plane_origin[1]), float(plane_origin[2])), + gp_Dir(float(plane_normal[0]), float(plane_normal[1]), float(plane_normal[2])), + ) + + n = np.asarray(plane_normal, dtype=np.float64) + n = n / (float(np.linalg.norm(n)) or 1.0) + helper = np.array([1.0, 0.0, 0.0]) if abs(n[0]) < 0.9 else np.array([0.0, 1.0, 0.0]) + u_ax = np.cross(n, helper) + u_ax /= float(np.linalg.norm(u_ax)) or 1.0 + v_ax = np.cross(n, u_ax) + o_arr = np.asarray(plane_origin, dtype=np.float64) + + component = _triangle_components(tris) + n_comp = int(component.max()) + 1 if component.size else 0 + + pnt_cache: dict[int, gp_Pnt] = {} + + def _occ_pnt(idx: int) -> gp_Pnt: + cached = pnt_cache.get(idx) + if cached is None: + v = pts[idx] + cached = gp_Pnt(float(v[0]), float(v[1]), float(v[2])) + pnt_cache[idx] = cached + return cached + + def _signed_area(loop: list[int]) -> float: + rel = pts[loop[:-1]] - o_arr + u = rel @ u_ax + v = rel @ v_ax + return 0.5 * float(np.sum(u * np.roll(v, -1) - np.roll(u, -1) * v)) + + def _build_wire(loop: list[int]): + wb = BRepBuilderAPI_MakeWire() + for k in range(len(loop) - 1): + edge_builder = BRepBuilderAPI_MakeEdge( + _occ_pnt(loop[k]), _occ_pnt(loop[k + 1]) + ) + if not edge_builder.IsDone(): + err_msg = "BRepBuilderAPI_MakeEdge failed for boundary segment" + raise ShellCreationError(err_msg) + wb.Add(edge_builder.Edge()) + if not wb.IsDone(): + err_msg = "BRepBuilderAPI_MakeWire failed for boundary loop" + raise ShellCreationError(err_msg) + return wb.Wire() + + def _build_component_face(loops: list[list[int]]): + areas = [_signed_area(L) for L in loops] + outer_local = int(np.argmax([abs(a) for a in areas])) + if areas[outer_local] < 0: + loops = [list(reversed(L)) for L in loops] + wires = [_build_wire(L) for L in loops] + face_builder = BRepBuilderAPI_MakeFace(plane, wires[outer_local]) + for k, w in enumerate(wires): + if k == outer_local: + continue + face_builder.Add(cast_wire(w.Reversed())) + if not face_builder.IsDone(): + err_msg = "BRepBuilderAPI_MakeFace failed building a planar face" + raise ShellCreationError(err_msg) + return face_builder.Face() + + component_faces = [] + for ci in range(n_comp): + comp_tris = tris[np.where(component == ci)[0]] + loops = _walk_boundary_loops(comp_tris) + if loops: + component_faces.append(_build_component_face(loops)) + + if not component_faces: + err_msg = "No planar face could be built from the triangle group" + raise ShellCreationError(err_msg) + + if len(component_faces) == 1: + return CadShape(component_faces[0]) + + builder = BRep_Builder() + shell = TopoDS_Shell() + builder.MakeShell(shell) + for f in component_faces: + builder.Add(shell, f) + return CadShape(shell) + + +def mesh_to_shell_brep( + points: npt.NDArray[np.float64], + triangles: npt.NDArray[np.int64], +) -> CadShape: + """Convert a triangle mesh to a shell with one planar BREP face per triangle. + + Slower than :func:`mesh_to_shape` (which uses a single tessellated face) + but produces a shell with *real* geometric surfaces — required whenever + the resulting shell is used as a cutting tool in boolean ops + (``BRepAlgoAPI_Cut``, ``Workplane.split()``, etc.), which refuse a + tessellated-only face. + + :param points: ``(N, 3)`` array of vertex coordinates + :param triangles: ``(M, 3)`` array of 0-indexed triangle vertex indices + :return: wrapped ``TopoDS_Shell`` with one planar face per triangle + :raises ShellCreationError: if any triangle cannot be built into a face + """ + require_cad() + from OCP.BRep import BRep_Builder # noqa: PLC0415 + from OCP.BRepBuilderAPI import ( # noqa: PLC0415 + BRepBuilderAPI_MakeEdge, + BRepBuilderAPI_MakeFace, + BRepBuilderAPI_MakeWire, + ) + from OCP.gp import gp_Pnt # noqa: PLC0415 + from OCP.TopoDS import TopoDS_Shell # noqa: PLC0415 + + pts = np.asarray(points, dtype=np.float64) + tris = np.asarray(triangles, dtype=np.int64) + if tris.size == 0: + err_msg = "Cannot build a shell from an empty triangle list" + raise ShellCreationError(err_msg) + + occ_points = [gp_Pnt(float(p[0]), float(p[1]), float(p[2])) for p in pts] + + builder = BRep_Builder() + shell = TopoDS_Shell() + builder.MakeShell(shell) + + try: + for a, b, c in tris: + e1 = BRepBuilderAPI_MakeEdge(occ_points[int(a)], occ_points[int(b)]).Edge() + e2 = BRepBuilderAPI_MakeEdge(occ_points[int(b)], occ_points[int(c)]).Edge() + e3 = BRepBuilderAPI_MakeEdge(occ_points[int(c)], occ_points[int(a)]).Edge() + wire = BRepBuilderAPI_MakeWire(e1, e2, e3).Wire() + face = BRepBuilderAPI_MakeFace(wire).Face() + builder.Add(shell, face) + except Exception as err: + err_msg = ( + "Failed to build the OCCT shell from the mesh; " + "try to increase the resolution or adjust bounds." + ) + raise ShellCreationError(err_msg) from err + + return CadShape(shell) + + +def mesh_to_sewn_shell( + points: npt.NDArray[np.float64], + triangles: npt.NDArray[np.int64], + tolerance: float | None = None, +) -> CadShape: + """Convert a triangle mesh to a sewn shell with shared edges. + + Builds one planar BREP face per triangle then sews them via + :class:`BRepBuilderAPI_Sewing` so coincident edges/vertices are merged + into shared topology. The resulting shell has valid topology, which + makes any subsequent boolean op (``Cut``, ``Common``, ``Fuse``) tractable + — unlike :func:`mesh_to_shell_brep` whose triangle-soup output is + pathologically slow for booleans. + + :param points: ``(N, 3)`` array of vertex coordinates + :param triangles: ``(M, 3)`` array of 0-indexed triangle vertex indices + :param tolerance: sewing tolerance; defaults to ``1e-6 * bbox_diag`` + :return: wrapped sewn ``TopoDS_Shape`` (Shell, or Compound of shells if + the input is disconnected) + :raises ShellCreationError: if any triangle cannot be built into a face + """ + require_cad() + from OCP.BRepBuilderAPI import ( # noqa: PLC0415 + BRepBuilderAPI_MakeEdge, + BRepBuilderAPI_MakeFace, + BRepBuilderAPI_MakeWire, + BRepBuilderAPI_Sewing, + ) + from OCP.gp import gp_Pnt # noqa: PLC0415 + + pts = np.asarray(points, dtype=np.float64) + tris = np.asarray(triangles, dtype=np.int64) + if tris.size == 0: + err_msg = "Cannot build a shell from an empty triangle list" + raise ShellCreationError(err_msg) + + if tolerance is None: + bbox_diag = float(np.linalg.norm(pts.max(axis=0) - pts.min(axis=0))) + tolerance = max(1e-9, 1e-6 * bbox_diag) + + occ_points = [gp_Pnt(float(p[0]), float(p[1]), float(p[2])) for p in pts] + + sewing = BRepBuilderAPI_Sewing(tolerance) + try: + for a, b, c in tris: + e1 = BRepBuilderAPI_MakeEdge(occ_points[int(a)], occ_points[int(b)]).Edge() + e2 = BRepBuilderAPI_MakeEdge(occ_points[int(b)], occ_points[int(c)]).Edge() + e3 = BRepBuilderAPI_MakeEdge(occ_points[int(c)], occ_points[int(a)]).Edge() + wire = BRepBuilderAPI_MakeWire(e1, e2, e3).Wire() + face = BRepBuilderAPI_MakeFace(wire).Face() + sewing.Add(face) + except Exception as err: + err_msg = ( + "Failed to build the OCCT shell from the mesh; " + "try to increase the resolution or adjust bounds." + ) + raise ShellCreationError(err_msg) from err + + sewing.Perform() + return CadShape(sewing.SewedShape()) diff --git a/microgen/cad/primitives.py b/microgen/cad/primitives.py new file mode 100644 index 00000000..567b9c36 --- /dev/null +++ b/microgen/cad/primitives.py @@ -0,0 +1,233 @@ +"""OCCT primitive builders. + +Free factory functions producing :class:`CadShape` instances for the core +parametric primitives (box, sphere, cylinder, capsule, ellipsoid, +polyhedron, extruded polygon). These are the back-ends the implicit shape +classes (``Box``, ``Sphere``, …) call from their ``generate_cad`` methods. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np + +from ._install import require_cad +from .shape import CadShape, ShellCreationError, _topods_cast + +if TYPE_CHECKING: + from collections.abc import Sequence + + +def make_box(dim: Sequence[float], center: Sequence[float]) -> CadShape: + """Axis-aligned box of size ``dim`` centered at ``center``.""" + require_cad() + from OCP.BRepPrimAPI import BRepPrimAPI_MakeBox # noqa: PLC0415 + from OCP.gp import gp_Pnt # noqa: PLC0415 + + dx, dy, dz = (float(d) for d in dim) + cx, cy, cz = (float(c) for c in center) + corner = gp_Pnt(cx - dx / 2.0, cy - dy / 2.0, cz - dz / 2.0) + return CadShape(BRepPrimAPI_MakeBox(corner, dx, dy, dz).Shape()) + + +def make_sphere(radius: float, center: Sequence[float]) -> CadShape: + """Sphere of given radius at ``center``.""" + require_cad() + from OCP.BRepPrimAPI import BRepPrimAPI_MakeSphere # noqa: PLC0415 + from OCP.gp import gp_Pnt # noqa: PLC0415 + + pnt = gp_Pnt(float(center[0]), float(center[1]), float(center[2])) + return CadShape(BRepPrimAPI_MakeSphere(pnt, float(radius)).Shape()) + + +def make_cylinder( + radius: float, + height: float, + center: Sequence[float], + axis: Sequence[float] = (1.0, 0.0, 0.0), +) -> CadShape: + """Cylinder of given radius and height centered at ``center`` along ``axis``.""" + require_cad() + from OCP.BRepPrimAPI import BRepPrimAPI_MakeCylinder # noqa: PLC0415 + from OCP.gp import gp_Ax2, gp_Dir, gp_Pnt # noqa: PLC0415 + + h = float(height) + ax_vec = np.asarray(axis, dtype=np.float64) + ax_vec = ax_vec / np.linalg.norm(ax_vec) + base = ( + float(center[0]) - h / 2.0 * ax_vec[0], + float(center[1]) - h / 2.0 * ax_vec[1], + float(center[2]) - h / 2.0 * ax_vec[2], + ) + ax = gp_Ax2( + gp_Pnt(*base), + gp_Dir(float(ax_vec[0]), float(ax_vec[1]), float(ax_vec[2])), + ) + return CadShape(BRepPrimAPI_MakeCylinder(ax, float(radius), h).Shape()) + + +def make_capsule( + radius: float, + height: float, + center: Sequence[float], +) -> CadShape: + """Capsule (cylinder along X with hemispherical caps).""" + require_cad() + from OCP.BRepPrimAPI import ( # noqa: PLC0415 + BRepPrimAPI_MakeCylinder, + BRepPrimAPI_MakeSphere, + ) + from OCP.gp import gp_Ax2, gp_Dir, gp_Pnt # noqa: PLC0415 + + cx, cy, cz = (float(c) for c in center) + h = float(height) + r = float(radius) + base_axis = gp_Ax2(gp_Pnt(cx - h / 2.0, cy, cz), gp_Dir(1.0, 0.0, 0.0)) + cyl = BRepPrimAPI_MakeCylinder(base_axis, r, h).Shape() + left = BRepPrimAPI_MakeSphere(gp_Pnt(cx - h / 2.0, cy, cz), r).Shape() + right = BRepPrimAPI_MakeSphere(gp_Pnt(cx + h / 2.0, cy, cz), r).Shape() + return CadShape(cyl).fuse(CadShape(left)).fuse(CadShape(right)) + + +def make_ellipsoid(radii: Sequence[float], center: Sequence[float]) -> CadShape: + """Ellipsoid of the given axis-aligned radii at ``center``. + + Built as a unit sphere transformed by a non-uniform scaling matrix. + """ + require_cad() + from OCP.BRepBuilderAPI import BRepBuilderAPI_GTransform # noqa: PLC0415 + from OCP.BRepPrimAPI import BRepPrimAPI_MakeSphere # noqa: PLC0415 + from OCP.gp import gp_GTrsf, gp_Mat, gp_Pnt, gp_XYZ # noqa: PLC0415 + + rx, ry, rz = (float(r) for r in radii) + cx, cy, cz = (float(c) for c in center) + sphere = BRepPrimAPI_MakeSphere(gp_Pnt(0.0, 0.0, 0.0), 1.0).Shape() + + gtrsf = gp_GTrsf() + gtrsf.SetVectorialPart( + gp_Mat( + rx, + 0.0, + 0.0, + 0.0, + ry, + 0.0, + 0.0, + 0.0, + rz, + ), + ) + gtrsf.SetTranslationPart(gp_XYZ(cx, cy, cz)) + return CadShape(BRepBuilderAPI_GTransform(sphere, gtrsf, True).Shape()) + + +def make_polyhedron( + vertices: Sequence[Sequence[float]], + faces_ixs: Sequence[Sequence[int]], + center: Sequence[float] = (0.0, 0.0, 0.0), +) -> CadShape: + """Polyhedron from vertex list and face→vertex index list. + + Each face's vertex index list must be closed (last == first). + + Uses ``BRepBuilderAPI_Sewing`` to stitch the independently-constructed + faces into a shared-edge shell with consistent outward orientation, + then ``BRepBuilderAPI_MakeSolid`` to close it. (Raw ``BRep_Builder.Add`` + does not orient; the resulting solid would have mixed-sign volume.) + """ + require_cad() + from OCP.BRepBuilderAPI import ( # noqa: PLC0415 + BRepBuilderAPI_MakeEdge, + BRepBuilderAPI_MakeFace, + BRepBuilderAPI_MakeSolid, + BRepBuilderAPI_MakeWire, + BRepBuilderAPI_Sewing, + ) + from OCP.gp import gp_Pnt # noqa: PLC0415 + from OCP.ShapeFix import ShapeFix_Solid # noqa: PLC0415 + from OCP.TopAbs import TopAbs_SHELL # noqa: PLC0415 + from OCP.TopExp import TopExp_Explorer # noqa: PLC0415 + + cx, cy, cz = (float(c) for c in center) + points = [ + gp_Pnt(float(v[0]) + cx, float(v[1]) + cy, float(v[2]) + cz) for v in vertices + ] + + sewing = BRepBuilderAPI_Sewing() + for ixs in faces_ixs: + wire_builder = BRepBuilderAPI_MakeWire() + for i1, i2 in zip(ixs, ixs[1:], strict=False): + edge_builder = BRepBuilderAPI_MakeEdge(points[i1], points[i2]) + if not edge_builder.IsDone(): + err_msg = "BRepBuilderAPI_MakeEdge failed for polyhedron face" + raise ShellCreationError(err_msg) + wire_builder.Add(edge_builder.Edge()) + if not wire_builder.IsDone(): + err_msg = "BRepBuilderAPI_MakeWire failed for polyhedron face" + raise ShellCreationError(err_msg) + face_builder = BRepBuilderAPI_MakeFace(wire_builder.Wire()) + if not face_builder.IsDone(): + err_msg = "BRepBuilderAPI_MakeFace failed for polyhedron face" + raise ShellCreationError(err_msg) + sewing.Add(face_builder.Face()) + sewing.Perform() + sewn = sewing.SewedShape() + + # Extract the shell (sewing may return it directly or wrapped in a compound). + exp = TopExp_Explorer(sewn, TopAbs_SHELL) + if not exp.More(): + err_msg = "Sewing did not produce a shell — check face connectivity" + raise ShellCreationError(err_msg) + shell = _topods_cast("Shell")(exp.Current()) + + solid_builder = BRepBuilderAPI_MakeSolid(shell) + if not solid_builder.IsDone(): + err_msg = "BRepBuilderAPI_MakeSolid failed; sewn shell is not closed" + raise ShellCreationError(err_msg) + fixer = ShapeFix_Solid(solid_builder.Solid()) + fixer.Perform() + return CadShape(fixer.Solid()) + + +def make_extruded_polygon( + list_corners: Sequence[tuple[float, float]], + height: float, + center: Sequence[float], +) -> CadShape: + """Extrude a 2D polygon (in the YZ plane) along the X axis.""" + require_cad() + from OCP.BRepBuilderAPI import ( # noqa: PLC0415 + BRepBuilderAPI_MakeEdge, + BRepBuilderAPI_MakeFace, + BRepBuilderAPI_MakeWire, + ) + from OCP.BRepPrimAPI import BRepPrimAPI_MakePrism # noqa: PLC0415 + from OCP.gp import gp_Pnt, gp_Vec # noqa: PLC0415 + + cx, cy, cz = (float(c) for c in center) + h = float(height) + x_base = cx - h / 2.0 + pts = [gp_Pnt(x_base, cy + float(y), cz + float(z)) for (y, z) in list_corners] + if (list_corners[0][0], list_corners[0][1]) != ( + list_corners[-1][0], + list_corners[-1][1], + ): + pts.append(pts[0]) + + wire_builder = BRepBuilderAPI_MakeWire() + for i in range(len(pts) - 1): + edge_builder = BRepBuilderAPI_MakeEdge(pts[i], pts[i + 1]) + if not edge_builder.IsDone(): + err_msg = "BRepBuilderAPI_MakeEdge failed for extruded polygon" + raise ShellCreationError(err_msg) + wire_builder.Add(edge_builder.Edge()) + if not wire_builder.IsDone(): + err_msg = "BRepBuilderAPI_MakeWire failed for extruded polygon" + raise ShellCreationError(err_msg) + face_builder = BRepBuilderAPI_MakeFace(wire_builder.Wire()) + if not face_builder.IsDone(): + err_msg = "BRepBuilderAPI_MakeFace failed for extruded polygon" + raise ShellCreationError(err_msg) + extruded = BRepPrimAPI_MakePrism(face_builder.Face(), gp_Vec(h, 0.0, 0.0)).Shape() + return CadShape(extruded) diff --git a/microgen/cad/shape.py b/microgen/cad/shape.py new file mode 100644 index 00000000..108be59c --- /dev/null +++ b/microgen/cad/shape.py @@ -0,0 +1,376 @@ +"""``CadShape`` wrapper and OCCT topology helpers. + +This module holds the ``CadShape`` thin wrapper around ``TopoDS_Shape`` plus +the low-level OCP utilities used by the rest of the CAD subpackage +(``_run_boolean``, ``_topods_cast``, ``_Centre``, ``_BBox``, +``ShellCreationError``). +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING, Any + +import numpy as np + +from ._install import require_cad # noqa: F401 (re-exported for back-compat) + +if TYPE_CHECKING: + from collections.abc import Sequence + from pathlib import Path + + from OCP.TopoDS import TopoDS_Shape + + +class _Centre(tuple): + """Tuple-like 3D point exposing ``.x``, ``.y``, ``.z`` and ``.to_tuple()``. + + Returned by :meth:`CadShape.center`. Mimics just enough of + ``cadquery.Vector`` to drop-in where existing code uses + ``shape.center().to_tuple()`` or ``shape.center().x``. + """ + + __slots__ = () + + def __new__(cls, x: float, y: float, z: float) -> _Centre: + """Create a 3-tuple ``(x, y, z)``.""" + return super().__new__(cls, (float(x), float(y), float(z))) + + @property + def x(self) -> float: + """X coordinate.""" + return self[0] + + @property + def y(self) -> float: + """Y coordinate.""" + return self[1] + + @property + def z(self) -> float: + """Z coordinate.""" + return self[2] + + def to_tuple(self) -> tuple[float, float, float]: + """Return ``(x, y, z)`` as a plain tuple.""" + return (self[0], self[1], self[2]) + + +class _BBox: + """Axis-aligned bounding box exposing ``xmin`` / ``xmax`` / …. + + Returned by :meth:`CadShape.bounding_box`. Also indexable as a 6-tuple + ``(xmin, ymin, zmin, xmax, ymax, zmax)`` matching OCCT's ``Bnd_Box.Get``. + """ + + __slots__ = ("xmax", "xmin", "ymax", "ymin", "zmax", "zmin") + + def __init__( + self, + xmin: float, + ymin: float, + zmin: float, + xmax: float, + ymax: float, + zmax: float, + ) -> None: + """Initialize from the 6 axis-aligned extents.""" + self.xmin = float(xmin) + self.ymin = float(ymin) + self.zmin = float(zmin) + self.xmax = float(xmax) + self.ymax = float(ymax) + self.zmax = float(zmax) + + @property + def diagonal_length(self) -> float: + """Length of the box's space diagonal.""" + dx = self.xmax - self.xmin + dy = self.ymax - self.ymin + dz = self.zmax - self.zmin + return float((dx * dx + dy * dy + dz * dz) ** 0.5) + + +class ShellCreationError(RuntimeError): + """Raised when a mesh cannot be converted into an OCCT shell.""" + + +def _run_boolean(op_cls: Any, a: CadShape, b: CadShape, label: str) -> TopoDS_Shape: + """Run an OCCT boolean op and raise on failure. + + Older OCP releases don't expose ``HasErrors()`` / ``IsDone()`` on the + ``BRepAlgoAPI_*`` classes; we probe via ``getattr`` and skip the check + when the API isn't available. + """ + op = op_cls(a.wrapped, b.wrapped) + has_errors = getattr(op, "HasErrors", None) + if callable(has_errors) and has_errors(): + err_msg = f"BRepAlgoAPI_{label} failed" + raise RuntimeError(err_msg) + return op.Shape() + + +def _topods_cast(name: str) -> Any: + """Return ``TopoDS.`` cast helper, tolerant of OCP version drift. + + Older OCP releases expose the static cast as ``TopoDS.Shell_s`` (pybind11 + ``_s`` convention); newer releases expose the unsuffixed ``TopoDS.Shell``. + Try the suffixed form first, fall back to unsuffixed. + """ + from OCP.TopoDS import TopoDS # noqa: PLC0415 + + return getattr(TopoDS, f"{name}_s", None) or getattr(TopoDS, name) + + +class CadShape: + """Thin wrapper around an OCCT ``TopoDS_Shape``. + + Preserves the ``.wrapped`` attribute name used by CadQuery so downstream + OCP calls (``BRepAlgoAPI_Fuse(a.wrapped, b.wrapped)``) keep working. + + ``_mesh_volume`` (optional) is a trusted volume in the source mesh's + units, set by mesh-derived constructors (e.g. the TPMS periodic shell) + where OCCT's surface-integral volume is unreliable on invalid topology. + :meth:`volume` prefers it over the OCCT integral when present. + """ + + __slots__ = ("_mesh_volume", "wrapped") + + def __init__(self, shape: TopoDS_Shape) -> None: + """Wrap an OCCT ``TopoDS_Shape``.""" + self.wrapped = shape + self._mesh_volume: float | None = None + + # -- transforms -------------------------------------------------------- + + def translate(self, offset: Sequence[float]) -> CadShape: + """Return a translated copy.""" + from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform # noqa: PLC0415 + from OCP.gp import gp_Trsf, gp_Vec # noqa: PLC0415 + + trsf = gp_Trsf() + trsf.SetTranslation( + gp_Vec(float(offset[0]), float(offset[1]), float(offset[2])) + ) + transformed = BRepBuilderAPI_Transform(self.wrapped, trsf, True).Shape() + return CadShape(transformed) + + def rotate( + self, + center: Sequence[float], + axis: Sequence[float], + angle_degrees: float, + ) -> CadShape: + """Return a rotated copy (angle in degrees, axis is a unit vector).""" + from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform # noqa: PLC0415 + from OCP.gp import gp_Ax1, gp_Dir, gp_Pnt, gp_Trsf # noqa: PLC0415 + + trsf = gp_Trsf() + ax = gp_Ax1( + gp_Pnt(float(center[0]), float(center[1]), float(center[2])), + gp_Dir(float(axis[0]), float(axis[1]), float(axis[2])), + ) + trsf.SetRotation(ax, float(np.deg2rad(angle_degrees))) + transformed = BRepBuilderAPI_Transform(self.wrapped, trsf, True).Shape() + return CadShape(transformed) + + def copy(self) -> CadShape: + """Return an independent copy (deep topology copy).""" + from OCP.BRepBuilderAPI import BRepBuilderAPI_Copy # noqa: PLC0415 + + return CadShape(BRepBuilderAPI_Copy(self.wrapped).Shape()) + + # -- boolean ops ------------------------------------------------------- + + def fuse(self, other: CadShape) -> CadShape: + """Boolean fusion: ``self ∪ other``.""" + from OCP.BRepAlgoAPI import BRepAlgoAPI_Fuse # noqa: PLC0415 + + return CadShape(_run_boolean(BRepAlgoAPI_Fuse, self, other, "Fuse")) + + def cut(self, other: CadShape) -> CadShape: + """Boolean difference: ``self \\ other``.""" + from OCP.BRepAlgoAPI import BRepAlgoAPI_Cut # noqa: PLC0415 + + return CadShape(_run_boolean(BRepAlgoAPI_Cut, self, other, "Cut")) + + def intersect(self, other: CadShape) -> CadShape: + """Boolean intersection: ``self ∩ other``.""" + from OCP.BRepAlgoAPI import BRepAlgoAPI_Common # noqa: PLC0415 + + return CadShape(_run_boolean(BRepAlgoAPI_Common, self, other, "Common")) + + # -- topology queries -------------------------------------------------- + + def solids(self) -> list[CadShape]: + """Enumerate contained solids.""" + from OCP.TopAbs import TopAbs_SOLID # noqa: PLC0415 + from OCP.TopExp import TopExp_Explorer # noqa: PLC0415 + + cast_solid = _topods_cast("Solid") + out: list[CadShape] = [] + exp = TopExp_Explorer(self.wrapped, TopAbs_SOLID) + while exp.More(): + out.append(CadShape(cast_solid(exp.Current()))) + exp.Next() + return out + + def vertices(self) -> list[tuple[float, float, float]]: + """Enumerate the vertex coordinates of the shape. + + Callers use this to check that a generated mesh has any vertices at + all (``assert np.any(shape.vertices())``). + """ + from OCP.BRep import BRep_Tool # noqa: PLC0415 + from OCP.TopAbs import TopAbs_VERTEX # noqa: PLC0415 + from OCP.TopExp import TopExp_Explorer # noqa: PLC0415 + + cast_vertex = _topods_cast("Vertex") + pnt = getattr(BRep_Tool, "Pnt_s", None) or BRep_Tool.Pnt + out: list[tuple[float, float, float]] = [] + exp = TopExp_Explorer(self.wrapped, TopAbs_VERTEX) + while exp.More(): + v = cast_vertex(exp.Current()) + p = pnt(v) + out.append((float(p.X()), float(p.Y()), float(p.Z()))) + exp.Next() + return out + + def faces(self) -> list[CadShape]: + """Enumerate the faces of the shape.""" + from OCP.TopAbs import TopAbs_FACE # noqa: PLC0415 + from OCP.TopExp import TopExp_Explorer # noqa: PLC0415 + + cast_face = _topods_cast("Face") + out: list[CadShape] = [] + exp = TopExp_Explorer(self.wrapped, TopAbs_FACE) + while exp.More(): + out.append(CadShape(cast_face(exp.Current()))) + exp.Next() + return out + + def is_closed(self) -> bool: + """Whether the shape is topologically closed. + + Solids are always closed; for shells/compounds we read OCCT's + per-shape ``Closed`` flag (set by ``BRep_Builder::IsClosed`` when the + shell was built from a watertight set of faces). + """ + from OCP.TopAbs import TopAbs_SOLID # noqa: PLC0415 + + if self.wrapped.ShapeType() == TopAbs_SOLID: + return True + return bool(self.wrapped.Closed()) + + def volume(self) -> float: + """Return the (unsigned) volume of the shape. + + OCCT's ``BRepGProp::VolumeProperties`` returns a *signed* volume that + depends on face orientation; mesh-built shells from + :func:`microgen.cad.mesh_to_shell_brep` can carry inverted orientation + and yield a negative value. We return ``abs(...)`` so volumes are + non-negative. + + If a mesh-derived volume was stashed on ``_mesh_volume`` AND the OCCT + solid is not valid (BRepCheck_Analyzer flags self-intersection / + non-manifold edges, common on raw marching-cubes input), we trust the + mesh volume — the OCCT surface integral on an invalid topology is + meaningless. + """ + from OCP.BRepCheck import BRepCheck_Analyzer # noqa: PLC0415 + from OCP.BRepGProp import BRepGProp # noqa: PLC0415 + from OCP.GProp import GProp_GProps # noqa: PLC0415 + + if ( + self._mesh_volume is not None + and not BRepCheck_Analyzer( + self.wrapped, + ).IsValid() + ): + return float(abs(self._mesh_volume)) + + props = GProp_GProps() + BRepGProp.VolumeProperties_s(self.wrapped, props) + return float(abs(props.Mass())) + + def center(self) -> _Centre: + """Return the volumetric center of mass. + + The result is a :class:`_Centre` — exposes ``.x``, ``.y``, ``.z``, + ``.to_tuple()``, and unpacks like a tuple. + """ + from OCP.BRepGProp import BRepGProp # noqa: PLC0415 + from OCP.GProp import GProp_GProps # noqa: PLC0415 + + props = GProp_GProps() + BRepGProp.VolumeProperties_s(self.wrapped, props) + c = props.CentreOfMass() + return _Centre(float(c.X()), float(c.Y()), float(c.Z())) + + def bounding_box(self) -> _BBox: + """Return the axis-aligned bounding box. + + The result exposes ``xmin`` / ``xmax`` / … attributes (see + :class:`_BBox`). + """ + from OCP.Bnd import Bnd_Box # noqa: PLC0415 + from OCP.BRepBndLib import BRepBndLib # noqa: PLC0415 + + box = Bnd_Box() + # AddOptimal uses exact geometric bounds (not cached triangulation). + BRepBndLib.AddOptimal_s(self.wrapped, box, True, True) + xmin, ymin, zmin, xmax, ymax, zmax = box.Get() + return _BBox(xmin, ymin, zmin, xmax, ymax, zmax) + + # -- exports ----------------------------------------------------------- + + def export_stl( + self, + path: str | Path, + linear_deflection: float = 0.01, + angular_deflection: float = 0.5, + *, + ascii_mode: bool = False, + ) -> None: + """Export to STL. Mesh is regenerated at the given deflection.""" + from OCP.BRepMesh import BRepMesh_IncrementalMesh # noqa: PLC0415 + from OCP.StlAPI import StlAPI_Writer # noqa: PLC0415 + + BRepMesh_IncrementalMesh( + self.wrapped, + float(linear_deflection), + False, + float(angular_deflection), + True, + ) + writer = StlAPI_Writer() + writer.ASCIIMode = bool(ascii_mode) + if not writer.Write(self.wrapped, str(path)): + err_msg = f"STL write failed for {path!r}" + raise RuntimeError(err_msg) + + def export_step(self, path: str | Path) -> None: + """Export to STEP (AP214).""" + from OCP.IFSelect import IFSelect_RetDone # noqa: PLC0415 + from OCP.STEPControl import ( # noqa: PLC0415 + STEPControl_AsIs, + STEPControl_Writer, + ) + + writer = STEPControl_Writer() + status = writer.Transfer(self.wrapped, STEPControl_AsIs) + if status != IFSelect_RetDone: + err_msg = f"STEP transfer failed with status {status!r}" + raise RuntimeError(err_msg) + status = writer.Write(str(path)) + if status != IFSelect_RetDone: + err_msg = f"STEP write failed with status {status!r}" + raise RuntimeError(err_msg) + + def export_brep(self, path: str | Path) -> None: + """Export to OCCT native BREP.""" + from OCP.BRepTools import BRepTools # noqa: PLC0415 + + ok = BRepTools.Write_s(self.wrapped, str(path)) + if not ok: + err_msg = f"BREP write failed for {path!r}" + raise RuntimeError(err_msg) diff --git a/microgen/cad/topo.py b/microgen/cad/topo.py new file mode 100644 index 00000000..293bad30 --- /dev/null +++ b/microgen/cad/topo.py @@ -0,0 +1,238 @@ +"""Topology assembly, exploration, and splitting helpers. + +Free functions for working with OCCT ``TopoDS_Solid`` / ``TopoDS_Compound`` +collections: building compounds, enumerating solids, splitting shapes by +tool, plane-face cutting tools, affine transforms, side selection, and +box-clipping a list of solids. + +These power the periodic split-and-translate algorithm in +``microgen.periodic`` and the rasterisation pipeline on ``Phase``. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING, Any + +import numpy as np +import numpy.typing as npt + +from ._install import require_cad +from .shape import CadShape, _topods_cast + +if TYPE_CHECKING: + from collections.abc import Iterable, Sequence + + +def make_compound(shapes: Iterable[CadShape]) -> CadShape: + """Assemble shapes into a single OCCT ``TopoDS_Compound``.""" + require_cad() + from OCP.BRep import BRep_Builder # noqa: PLC0415 + from OCP.TopoDS import TopoDS_Compound # noqa: PLC0415 + + builder = BRep_Builder() + compound = TopoDS_Compound() + builder.MakeCompound(compound) + for s in shapes: + builder.Add(compound, s.wrapped) + return CadShape(compound) + + +def make_compound_from_solids(solids: Iterable[Any]) -> CadShape: + """Assemble raw OCCT ``TopoDS_Shape`` solids (not ``CadShape``) into a compound.""" + require_cad() + from OCP.BRep import BRep_Builder # noqa: PLC0415 + from OCP.TopoDS import TopoDS_Compound # noqa: PLC0415 + + builder = BRep_Builder() + compound = TopoDS_Compound() + builder.MakeCompound(compound) + for s in solids: + shape = s.wrapped if hasattr(s, "wrapped") else s + builder.Add(compound, shape) + return CadShape(compound) + + +def enumerate_solids(shape: CadShape) -> list[Any]: + """Return the list of ``TopoDS_Solid`` inside a shape (empty if none).""" + require_cad() + from OCP.TopAbs import TopAbs_SOLID # noqa: PLC0415 + from OCP.TopExp import TopExp_Explorer # noqa: PLC0415 + + cast_solid = _topods_cast("Solid") + out: list[Any] = [] + exp = TopExp_Explorer(shape.wrapped, TopAbs_SOLID) + while exp.More(): + out.append(cast_solid(exp.Current())) + exp.Next() + return out + + +def split_shape( + shape: CadShape, + tool: CadShape | Iterable[CadShape], + *, + fuzzy_value: float = 1e-4, +) -> CadShape: + """Split *shape* by *tool* using OCCT's ``BRepAlgoAPI_Splitter``. + + The result is a :class:`CadShape` wrapping a ``TopoDS_Compound`` that + contains the sub-shapes produced by the split. Use + :func:`enumerate_solids` to iterate over the resulting solids. + + :param tool: a single :class:`CadShape` or an iterable of them. Pass + multiple tools when each tool is a separate ``TopoDS_Shell`` and + you want to keep them topologically distinct — the OCCT splitter + treats each list entry as one cutting tool. This matters because + ``BRepAlgoAPI_Fuse`` on two shells *decomposes* them into a + compound of per-face shells, which then defeats the splitter. + :param fuzzy_value: tolerance forwarded to ``SetFuzzyValue`` so OCCT + recognises near-coincident geometry as touching. Necessary when + a tool is a tessellated shell whose boundary edges lie on + ``shape``'s planar faces only up to a few microns of drift. + """ + require_cad() + from OCP.BRepAlgoAPI import BRepAlgoAPI_Splitter # noqa: PLC0415 + from OCP.TopTools import TopTools_ListOfShape # noqa: PLC0415 + + args = TopTools_ListOfShape() + args.Append(shape.wrapped) + tools = TopTools_ListOfShape() + if isinstance(tool, CadShape): + tools.Append(tool.wrapped) + else: + for t in tool: + tools.Append(t.wrapped) + + splitter = BRepAlgoAPI_Splitter() + splitter.SetArguments(args) + splitter.SetTools(tools) + if fuzzy_value > 0: + splitter.SetFuzzyValue(float(fuzzy_value)) + splitter.Build() + return CadShape(splitter.Shape()) + + +def make_plane_face( + base_pnt: Sequence[float], + direction: Sequence[float], + half_size: float = 1.0e6, +) -> CadShape: + """Build a large planar face used as a cutting tool. + + :param base_pnt: a point on the plane + :param direction: plane normal + :param half_size: half-edge of the square face (default large so the plane + reaches far outside any realistic shape) + """ + require_cad() + from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeFace # noqa: PLC0415 + from OCP.gp import gp_Ax3, gp_Dir, gp_Pln, gp_Pnt # noqa: PLC0415 + + pnt = gp_Pnt(float(base_pnt[0]), float(base_pnt[1]), float(base_pnt[2])) + nrm = gp_Dir(float(direction[0]), float(direction[1]), float(direction[2])) + plane = gp_Pln(gp_Ax3(pnt, nrm)) + face = BRepBuilderAPI_MakeFace( + plane, + -float(half_size), + float(half_size), + -float(half_size), + float(half_size), + ).Face() + return CadShape(face) + + +def transform_geometry(shape: CadShape, matrix: npt.NDArray[np.float64]) -> CadShape: + """Apply a 3x4 affine matrix (linear + translation). + + Wraps OCCT ``BRepBuilderAPI_GTransform``. + + :param matrix: ``(3, 4)`` array; rows are ``[a b c tx; d e f ty; g h i tz]``. + """ + require_cad() + from OCP.BRepBuilderAPI import BRepBuilderAPI_GTransform # noqa: PLC0415 + from OCP.gp import gp_GTrsf, gp_Mat, gp_XYZ # noqa: PLC0415 + + m = np.asarray(matrix, dtype=np.float64) + gtrsf = gp_GTrsf() + gtrsf.SetVectorialPart( + gp_Mat( + float(m[0, 0]), + float(m[0, 1]), + float(m[0, 2]), + float(m[1, 0]), + float(m[1, 1]), + float(m[1, 2]), + float(m[2, 0]), + float(m[2, 1]), + float(m[2, 2]), + ), + ) + gtrsf.SetTranslationPart(gp_XYZ(float(m[0, 3]), float(m[1, 3]), float(m[2, 3]))) + return CadShape(BRepBuilderAPI_GTransform(shape.wrapped, gtrsf, True).Shape()) + + +def translate_solid(solid: Any, offset: Sequence[float]) -> Any: + """Translate a raw OCCT solid/shape by ``offset`` and return the same type.""" + require_cad() + from OCP.BRepBuilderAPI import BRepBuilderAPI_Transform # noqa: PLC0415 + from OCP.gp import gp_Trsf, gp_Vec # noqa: PLC0415 + + shape = solid.wrapped if hasattr(solid, "wrapped") else solid + trsf = gp_Trsf() + trsf.SetTranslation( + gp_Vec(float(offset[0]), float(offset[1]), float(offset[2])), + ) + return BRepBuilderAPI_Transform(shape, trsf, True).Shape() + + +def solid_center(shape: Any) -> tuple[float, float, float]: + """Return the volumetric center of mass of a raw OCCT shape/solid.""" + require_cad() + from OCP.BRepGProp import BRepGProp # noqa: PLC0415 + from OCP.GProp import GProp_GProps # noqa: PLC0415 + + s = shape.wrapped if hasattr(shape, "wrapped") else shape + props = GProp_GProps() + BRepGProp.VolumeProperties_s(s, props) + com = props.CentreOfMass() + return (float(com.X()), float(com.Y()), float(com.Z())) + + +def select_solids_on_side( + shape: CadShape, + base_pnt: Sequence[float], + side_direction: Sequence[float], +) -> list[Any]: + """Enumerate a compound's solids; keep those on the positive side of a plane. + + The plane passes through ``base_pnt`` normal to ``side_direction``. + Matches CadQuery's ``.solids(">X")`` / ``.solids(" 0: + selected.append(solid) + return selected + + +def intersect_solids_with_box(solids: Iterable[Any], box: CadShape) -> CadShape: + """Intersect each solid with *box* and fuse the results into a ``CadShape``. + + Returns a :class:`CadShape` wrapping a compound (possibly empty). + """ + require_cad() + from OCP.BRepAlgoAPI import BRepAlgoAPI_Common # noqa: PLC0415 + + parts: list[Any] = [] + for solid in solids: + s = solid.wrapped if hasattr(solid, "wrapped") else solid + common = BRepAlgoAPI_Common(s, box.wrapped).Shape() + # Keep if it contains at least one solid. + if enumerate_solids(CadShape(common)): + parts.append(common) + return make_compound_from_solids(parts) diff --git a/microgen/operations.py b/microgen/operations.py index 1fa167d9..1a3982bc 100644 --- a/microgen/operations.py +++ b/microgen/operations.py @@ -126,39 +126,6 @@ def rotate_euler( return rotate(obj, center, rotation) -def rotate_pv_euler( - obj: pv.PolyData, - center: Sequence[float], - angles_or_rotation: Sequence[float] | Rotation, -) -> pv.PolyData: - """Rotate object according to ZXZ Euler angle convention. - - :param obj: Object to rotate - :param center: numpy array (x, y, z) - :param angles_or_rotation: list of Euler angles (psi, theta, phi) in - degrees, or a scipy ``Rotation`` object - - :return: Rotated object - """ - if isinstance(angles_or_rotation, Rotation): - rotation = angles_or_rotation - else: - rotation = Rotation.from_euler("ZXZ", angles_or_rotation, degrees=True) - - rotvec = rotation.as_rotvec(degrees=True) - angle = np.linalg.norm(rotvec) - if angle == 0: - return obj - axis = rotvec / angle - - return obj.rotate_vector( - vector=axis, - angle=angle, - point=tuple(center), - inplace=False, - ) - - def rescale(shape: CadShape, scale: float | tuple[float, float, float]) -> CadShape: """Rescale given object according to scale parameters [dim_x, dim_y, dim_z].""" from .phase import Phase # noqa: PLC0415 diff --git a/microgen/rve.py b/microgen/rve.py index 1ce5355d..a9c2f02a 100644 --- a/microgen/rve.py +++ b/microgen/rve.py @@ -1,12 +1,17 @@ """Representative Volume Element (RVE). -The ``Rve.box`` attribute is a :class:`~microgen.cad.CadShape` wrapping an -OCCT box; it is built lazily on first access and requires the ``[cad]`` -extra (``cadquery-ocp-novtk``). +Frozen, immutable container for the RVE bounding box and its periodicity flags. + +``Rve.box`` is a :class:`~microgen.cad.CadShape` wrapping an OCCT box; it is +built lazily on first access and requires the ``[cad]`` extra. ``Rve.grid`` +returns a :class:`pyvista.StructuredGrid` aligned with the cell — used by +implicit shapes to sample SDF/level-set fields. """ from __future__ import annotations +from dataclasses import dataclass +from functools import cached_property from typing import TYPE_CHECKING import numpy as np @@ -16,81 +21,144 @@ if TYPE_CHECKING: from collections.abc import Sequence - from .cad import CadShape - - Vector3DType = tuple[float, float, float] | Sequence[float] + import numpy.typing as npt + import pyvista as pv + from .cad import CadShape + Vector3DType = ( + tuple[float, float, float] | Sequence[float] | npt.NDArray[np.float64] + ) + ResolutionType = int | tuple[int, int, int] | Sequence[int] + + +def _validate_center(center: object) -> np.ndarray: + if isinstance(center, (tuple, list)) and len(center) == _DIM: + return np.asarray(center, dtype=float) + if isinstance(center, np.ndarray) and center.shape == (_DIM,): + return center.astype(float, copy=True) + err_msg = f"center must be an array or Sequence of length {_DIM}" + raise ValueError(err_msg) + + +def _validate_dim(dim: object) -> np.ndarray: + if isinstance(dim, (int, float)) and not isinstance(dim, bool): + arr = np.array([dim] * _DIM, dtype=float) + elif isinstance(dim, (tuple, list)) and len(dim) == _DIM: + arr = np.asarray(dim, dtype=float) + elif isinstance(dim, np.ndarray) and dim.shape == (_DIM,): + arr = dim.astype(float, copy=True) + else: + err_msg = f"dim must be an array or Sequence of length {_DIM}" + raise ValueError(err_msg) + if np.any(arr <= 0): + err_msg = f"dimensions of the RVE must be greater than 0, got {arr.tolist()}" + raise ValueError(err_msg) + return arr + + +def _validate_pbc(pbc: object) -> tuple[bool, bool, bool]: + if isinstance(pbc, bool): + return (pbc, pbc, pbc) + if isinstance(pbc, (tuple, list)) and len(pbc) == _DIM: + return (bool(pbc[0]), bool(pbc[1]), bool(pbc[2])) + err_msg = f"pbc must be a bool or Sequence[bool] of length {_DIM}" + raise ValueError(err_msg) + + +@dataclass(frozen=True, init=False, eq=False, repr=False) class Rve: - """Representative Volume Element (RVE). + """Representative Volume Element (RVE) — frozen. :param center: center of the RVE - :param dim: dimensions of the RVE + :param dim: dimensions of the RVE (scalar → cube) + :param pbc: periodic-boundary-condition flags per axis ``(x, y, z)``; + defaults to fully periodic. A single ``bool`` is broadcast to all axes. """ + center: np.ndarray + dim: np.ndarray + pbc: tuple[bool, bool, bool] + def __init__( self: Rve, center: Vector3DType = (0, 0, 0), dim: float | Vector3DType = 1, + pbc: bool | tuple[bool, bool, bool] | Sequence[bool] = (True, True, True), ) -> None: """Initialize the RVE.""" - if isinstance(center, (tuple, list)) and len(center) == _DIM: - self.center = np.array(center) - elif isinstance(center, np.ndarray) and center.shape == (_DIM,): - self.center = center - else: - err_msg = f"center must be an array or Sequence of length {_DIM}" - raise ValueError(err_msg) - - if isinstance(dim, (int, float)): - self.dim = np.array([dim for _ in range(_DIM)]) - elif isinstance(dim, (tuple, list)) and len(dim) == _DIM: - self.dim = np.array(dim) - elif isinstance(dim, np.ndarray) and dim.shape == (_DIM,): - self.dim = dim - else: - err_msg = f"dim must be an array or Sequence of length {_DIM}" - raise ValueError(err_msg) - - if np.any(self.dim <= 0): - err_msg = f"dimensions of the RVE must be greater than 0, got {self.dim}" - raise ValueError(err_msg) - - self.min_point = self.center - 0.5 * self.dim - self.max_point = self.center + 0.5 * self.dim + object.__setattr__(self, "center", _validate_center(center)) + object.__setattr__(self, "dim", _validate_dim(dim)) + object.__setattr__(self, "pbc", _validate_pbc(pbc)) + + @cached_property + def min_point(self: Rve) -> np.ndarray: + """Min corner ``center - 0.5 * dim``.""" + return self.center - 0.5 * self.dim + + @cached_property + def max_point(self: Rve) -> np.ndarray: + """Max corner ``center + 0.5 * dim``.""" + return self.center + 0.5 * self.dim + + @cached_property + def box(self: Rve) -> CadShape: + """Return a :class:`~microgen.cad.CadShape` box of the RVE (cached). - self.is_matrix = False - self.matrix_number = 0 + Requires the ``[cad]`` extra. + """ + from .cad import make_box # noqa: PLC0415 - self._cached_box: CadShape | None = None + return make_box(tuple(self.dim), tuple(self.center)) - @property - def box(self) -> CadShape: - """Return a :class:`~microgen.cad.CadShape` box of the RVE (cached). + def grid(self: Rve, resolution: ResolutionType) -> pv.StructuredGrid: + """Return a structured grid aligned with the RVE. - Requires the ``[cad]`` extra. + :param resolution: points per axis — int (broadcast to all 3 axes) or + length-3 sequence. The grid spans ``min_point`` → ``max_point`` + with endpoints included on every axis. """ - if self._cached_box is None: - from .cad import make_box # noqa: PLC0415 + import pyvista as pv # noqa: PLC0415 + + if isinstance(resolution, int): + nx, ny, nz = resolution, resolution, resolution + elif isinstance(resolution, (tuple, list)) and len(resolution) == _DIM: + nx, ny, nz = (int(r) for r in resolution) + elif isinstance(resolution, np.ndarray) and resolution.shape == (_DIM,): + nx, ny, nz = (int(r) for r in resolution.tolist()) + else: + err_msg = f"resolution must be an int or Sequence of length {_DIM}" + raise ValueError(err_msg) - self._cached_box = make_box(tuple(self.dim), tuple(self.center)) - return self._cached_box + xi = np.linspace(self.min_point[0], self.max_point[0], nx) + yi = np.linspace(self.min_point[1], self.max_point[1], ny) + zi = np.linspace(self.min_point[2], self.max_point[2], nz) + x, y, z = np.meshgrid(xi, yi, zi, indexing="ij") + return pv.StructuredGrid(x, y, z) - @box.setter - def box(self, value: CadShape) -> None: - self._cached_box = value + def __repr__(self: Rve) -> str: + return ( + f"Rve(center={tuple(self.center.tolist())}, " + f"dim={tuple(self.dim.tolist())}, pbc={self.pbc})" + ) @classmethod def from_min_max( cls: type[Rve], min_point: Vector3DType = (-0.5, -0.5, -0.5), max_point: Vector3DType = (0.5, 0.5, 0.5), + pbc: bool | tuple[bool, bool, bool] | Sequence[bool] = (True, True, True), ) -> Rve: """Generate a Rve from min and max corner points. :param min_point: ``(x_min, y_min, z_min)`` corner of the RVE :param max_point: ``(x_max, y_max, z_max)`` corner of the RVE + :param pbc: periodic-boundary-condition flags (see ``__init__``) """ lo = np.asarray(min_point, dtype=float) hi = np.asarray(max_point, dtype=float) - return cls(center=tuple(0.5 * (lo + hi)), dim=tuple(np.abs(hi - lo))) + return cls( + center=tuple(0.5 * (lo + hi)), + dim=tuple(np.abs(hi - lo)), + pbc=pbc, + ) diff --git a/microgen/shape/__init__.py b/microgen/shape/__init__.py index 2176527f..1f9f5218 100644 --- a/microgen/shape/__init__.py +++ b/microgen/shape/__init__.py @@ -48,7 +48,7 @@ TruncatedOctahedron, ) from .spinodoid import Spinodoid -from .tpms import CylindricalTpms, Infill, SphericalTpms, Tpms +from .tpms import CylindricalTpms, GradedInfill, Infill, SphericalTpms, Sweep, Tpms from .tpms_grading import NormedDistance if TYPE_CHECKING: @@ -162,6 +162,7 @@ def __init__(self: ShapeError, shape: str) -> None: "Ellipsoid", "ExtrudedPolygon", "FaceCenteredCubic", + "GradedInfill", "Infill", "NormedDistance", "Octahedron", @@ -173,6 +174,7 @@ def __init__(self: ShapeError, shape: str) -> None: "Sphere", "SphericalTpms", "Spinodoid", + "Sweep", "Tpms", "TruncatedCube", "TruncatedCuboctahedron", diff --git a/microgen/shape/_types.py b/microgen/shape/_types.py new file mode 100644 index 00000000..d6959833 --- /dev/null +++ b/microgen/shape/_types.py @@ -0,0 +1,29 @@ +"""Shared type aliases for the shape package. + +Centralising these here avoids duplicate definitions across ``shape.py``, +``tpms.py`` and downstream modules. All implicit shapes use the same +``Field`` callable and ``BoundsType`` AABB representation. +""" + +from __future__ import annotations + +from collections.abc import Callable + +import numpy as np +import numpy.typing as npt + +# Implicit scalar field: ``(x, y, z) -> array``, with negative values inside. +Field = Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]], + npt.NDArray[np.float64], +] + +# Axis-aligned bounding box: ``(xmin, xmax, ymin, ymax, zmin, zmax)``. +BoundsType = tuple[float, float, float, float, float, float] + +# Period of an intrinsically-periodic shape: ``(Lx, Ly, Lz)``. +# When set, ``field(x + Lx, y, z) == field(x, y, z)`` (and analogously for y, z). +# A ``None`` period means the shape is not intrinsically periodic. +PeriodType = tuple[float, float, float] + +__all__ = ["BoundsType", "Field", "PeriodType"] diff --git a/microgen/shape/box.py b/microgen/shape/box.py index 64fc46da..863447ef 100644 --- a/microgen/shape/box.py +++ b/microgen/shape/box.py @@ -21,6 +21,7 @@ if TYPE_CHECKING: from microgen.cad import CadShape from microgen.shape import KwargsGenerateType, Vector3DType + from microgen.shape.shape import BoundsType class Box(Shape): @@ -105,10 +106,24 @@ def generate_cad(self: Box, **_: KwargsGenerateType) -> CadShape: def generate_surface_mesh( self: Box, + bounds: BoundsType | None = None, + resolution: int | None = None, level: int = 0, **_: KwargsGenerateType, ) -> pv.PolyData: - """Generate a box VTK shape using the given parameters.""" + """Generate a box VTK shape using the given parameters. + + When ``bounds`` or ``resolution`` is explicitly provided, fall back to + the implicit (marching-cubes-on-SDF) base implementation so polymorphic + callers that ask for a specific sampling get what they asked for. + Otherwise the native :class:`pyvista.Box` (fixed 6-face quad mesh + controlled by ``level``) is used. + """ + if bounds is not None or resolution is not None: + return super().generate_surface_mesh( + bounds=bounds, + resolution=resolution if resolution is not None else 50, + ) box = pv.Box( bounds=( self.center[0] - 0.5 * self.dim[0], diff --git a/microgen/shape/periodic_shell.py b/microgen/shape/periodic_shell.py new file mode 100644 index 00000000..3c1a4e7c --- /dev/null +++ b/microgen/shape/periodic_shell.py @@ -0,0 +1,124 @@ +"""Sew a triangulated periodic surface into a closed OCCT shell. + +Given a marching-cubes (or otherwise periodic-aligned) triangle mesh of an +implicit field whose iso-surface meets the cell boundary on cap planes, +:func:`mesh_to_periodic_shell` groups boundary triangles per cap plane, +builds **one planar BREP face per cap** (via :func:`microgen.cad.mesh_to_planar_face`) +carrying the actual cap-wire trace, and sews them together with the +interior triangles into a closed shell suitable for boolean ops, STEP +export, and gmsh ``setPeriodic`` constraints. + +Shared by :class:`microgen.shape.tpms.Tpms` and +:class:`microgen.shape.spinodoid.Spinodoid` — the two TPMS/F-rep shapes +that produce periodic meshes that need this cap-aware sewing. + +Requires the optional ``[cad]`` install extra. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +import numpy.typing as npt + +if TYPE_CHECKING: + from collections.abc import Sequence + + from microgen.cad import CadShape + + from ._types import BoundsType + + +def mesh_to_periodic_shell( + points: npt.NDArray[np.float64], + triangles: npt.NDArray[np.int64], + bounds: BoundsType | Sequence[float], +) -> CadShape: + """Build a sewn OCCT shell from a triangulated surface in an AABB cell. + + Triangles whose three vertices share the mesh's exact extremum on an + AABB cap plane are grouped per face and converted to a single planar + BREP face via :func:`microgen.cad.mesh_to_planar_face` (so STEP shows + the actual gyroid/spinodoid cuts, not a bounding cube, and gmsh + ``setPeriodic`` can pair opposite cap faces by plane equation). + Remaining triangles become one planar BREP face per triangle (raw, + not pre-sewn — sewing must stitch them to the cap wires along the + seam). Everything is sewn into a closed shell via + :class:`OCP.BRepBuilderAPI.BRepBuilderAPI_Sewing`. + + :param points: ``(N, 3)`` vertex coordinates + :param triangles: ``(M, 3)`` vertex-index triplets + :param bounds: ``(xmin, xmax, ymin, ymax, zmin, zmax)`` of the cell + :return: :class:`microgen.cad.CadShape` wrapping the sewn ``TopoDS_Shell`` + """ + from OCP.BRepBuilderAPI import ( # noqa: PLC0415 + BRepBuilderAPI_MakeEdge, + BRepBuilderAPI_MakeFace, + BRepBuilderAPI_MakeWire, + BRepBuilderAPI_Sewing, + ) + from OCP.gp import gp_Pnt # noqa: PLC0415 + from OCP.TopAbs import TopAbs_FACE # noqa: PLC0415 + from OCP.TopExp import TopExp_Explorer # noqa: PLC0415 + + from microgen.cad import CadShape as _CadShape # noqa: PLC0415 + from microgen.cad import mesh_to_planar_face # noqa: PLC0415 + + pts = np.asarray(points, dtype=np.float64) + tris = np.asarray(triangles, dtype=np.int64).reshape(-1, 3) + + drift_tol = 1e-9 * float(max(abs(b) for b in bounds) or 1.0) + + consumed = np.zeros(tris.shape[0], dtype=bool) + on_plane: list[tuple[int, int, float, npt.NDArray[np.int64]]] = [] + for axis in range(3): + for sign in (-1, +1): + extremum = ( + float(pts[:, axis].max()) if sign > 0 else float(pts[:, axis].min()) + ) + expected = bounds[2 * axis + (1 if sign > 0 else 0)] + if abs(extremum - expected) > drift_tol: + continue + vert_on = pts[:, axis] == extremum + tri_on = np.all(vert_on[tris], axis=1) & ~consumed + if tri_on.any(): + on_plane.append((axis, sign, extremum, np.where(tri_on)[0])) + consumed |= tri_on + interior_idx = np.where(~consumed)[0] + + bbox_diag = float(np.linalg.norm(pts.max(axis=0) - pts.min(axis=0))) + sew_tol = max(1e-9, 1e-6 * bbox_diag) + sewing = BRepBuilderAPI_Sewing(sew_tol) + + for axis, sign, extremum, tri_idx in on_plane: + origin = [0.0, 0.0, 0.0] + origin[axis] = extremum + normal = [0.0, 0.0, 0.0] + normal[axis] = float(sign) + planar = mesh_to_planar_face(pts, tris[tri_idx], origin, normal) + exp = TopExp_Explorer(planar.wrapped, TopAbs_FACE) + while exp.More(): + sewing.Add(exp.Current()) + exp.Next() + + interior_tris = tris[interior_idx] + used_vertices = np.unique(interior_tris) + pnt_cache = { + int(i): gp_Pnt(float(pts[i, 0]), float(pts[i, 1]), float(pts[i, 2])) + for i in used_vertices + } + for a, b, c in interior_tris: + ga, gb, gc = pnt_cache[int(a)], pnt_cache[int(b)], pnt_cache[int(c)] + e1 = BRepBuilderAPI_MakeEdge(ga, gb).Edge() + e2 = BRepBuilderAPI_MakeEdge(gb, gc).Edge() + e3 = BRepBuilderAPI_MakeEdge(gc, ga).Edge() + wire = BRepBuilderAPI_MakeWire(e1, e2, e3).Wire() + face = BRepBuilderAPI_MakeFace(wire).Face() + sewing.Add(face) + + sewing.Perform() + return _CadShape(sewing.SewedShape()) + + +__all__ = ["mesh_to_periodic_shell"] diff --git a/microgen/shape/shape.py b/microgen/shape/shape.py index 7e153e74..fa074a39 100644 --- a/microgen/shape/shape.py +++ b/microgen/shape/shape.py @@ -8,7 +8,6 @@ from __future__ import annotations import itertools -from collections.abc import Callable from typing import TYPE_CHECKING import numpy as np @@ -17,18 +16,12 @@ from scipy.spatial.transform import Rotation from . import implicit_ops as _ops +from ._types import BoundsType, Field, PeriodType if TYPE_CHECKING: from microgen.cad import CadShape from microgen.shape import KwargsGenerateType, Vector3DType -Field = Callable[ - [npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]], - npt.NDArray[np.float64], -] - -BoundsType = tuple[float, float, float, float, float, float] - # Scalar-field name on the StructuredGrid used by the default mesh generators. _IMPLICIT_SCALAR = "implicit" @@ -91,6 +84,9 @@ class Shape: :param orientation: orientation of the shape :param func: implicit scalar field ``(x, y, z) -> array``, or ``None`` :param bounds: ``(xmin, xmax, ymin, ymax, zmin, zmax)`` or ``None`` + :param period: ``(Lx, Ly, Lz)`` if the field is intrinsically periodic + (``func(p + L) == func(p)`` along each axis), or ``None``. + Set by ``Tpms`` and ``Spinodoid`` from ``cell_size * repeat_cell``. """ def __init__( @@ -99,6 +95,7 @@ def __init__( orientation: Vector3DType | Rotation = (0, 0, 0), func: Field | None = None, bounds: BoundsType | None = None, + period: PeriodType | None = None, ) -> None: """Initialize the shape.""" self._center = center @@ -109,6 +106,7 @@ def __init__( ) self._func = func self._bounds = bounds + self._period: PeriodType | None = period # Cache of sampled structured grids keyed on (bounds, resolution). # Shared between generate_surface_mesh and generate_volume_mesh so # users calling both on the same instance only pay one N^3 field @@ -152,6 +150,17 @@ def bounds(self: Shape) -> BoundsType | None: """The bounding box ``(xmin, xmax, ymin, ymax, zmin, zmax)``, or ``None``.""" return self._bounds + @property + def period(self: Shape) -> PeriodType | None: + """The intrinsic period ``(Lx, Ly, Lz)`` if the field is periodic, else ``None``. + + When non-``None``, ``self.evaluate(x + Lx, y, z) == self.evaluate(x, y, z)`` + (and analogously for y, z) — i.e. periodicity is a data-structure + invariant of the field, not a runtime flag. ``Tpms`` and + ``Spinodoid`` set this from ``cell_size * repeat_cell``. + """ + return self._period + def require_func(self: Shape) -> Field: """Return ``_func`` or raise if not set.""" if self._func is None: @@ -297,10 +306,12 @@ def generate_cad( ) -> CadShape: """Generate a CAD shape. - The default implementation builds an OCCT tessellated BREP from the - implicit-field VTK mesh, via :func:`microgen.cad.mesh_to_shape` - (single ``TopoDS_Face`` carrying a ``Poly_Triangulation``). Subclasses - override this with native primitive construction. + The default implementation delegates to + :func:`microgen.cad.shape_to_cad`, which builds a tessellated OCCT BREP + from the implicit-field marching-cubes mesh. Concrete subclasses + with native primitive paths (``Box``, ``Sphere``, ``Cylinder``, + ``Capsule``, ``Ellipsoid``, ``Tpms``, ``Spinodoid`` …) override + this method with native OCCT construction. Requires the optional ``[cad]`` install extra (``cadquery-ocp-novtk``). @@ -308,32 +319,9 @@ def generate_cad( :param resolution: number of grid points per axis :return: :class:`microgen.cad.CadShape` wrapping an OCCT ``TopoDS_Shell`` """ - if self._func is None: - err_msg = ( - "No implicit field defined — subclasses must override generate_cad()" - ) - raise NotImplementedError(err_msg) - - from microgen.cad import mesh_to_shape # noqa: PLC0415 + from microgen.cad import shape_to_cad # noqa: PLC0415 - mesh = self.generate_surface_mesh(bounds=bounds, resolution=resolution) - if mesh.n_cells == 0: - err_msg = "Generated mesh is empty — check bounds and field function" - raise ValueError(err_msg) - - if not mesh.is_all_triangles: - mesh.triangulate(inplace=True) - triangles = mesh.faces.reshape(-1, 4)[:, 1:] - points = np.asarray(mesh.points, dtype=np.float64) - - try: - return mesh_to_shape(points, triangles) - except Exception as err: - err_msg = ( - "Failed to build the OCCT shell from the mesh; " - "try to increase the resolution or adjust bounds." - ) - raise ShellCreationError(err_msg) from err + return shape_to_cad(self, bounds=bounds, resolution=resolution) # ------------------------------------------------------------------ # Boolean operators (on implicit field) diff --git a/microgen/shape/sphere.py b/microgen/shape/sphere.py index 32e0010d..afb1cc85 100644 --- a/microgen/shape/sphere.py +++ b/microgen/shape/sphere.py @@ -18,6 +18,7 @@ if TYPE_CHECKING: from microgen.cad import CadShape from microgen.shape import KwargsGenerateType, Vector3DType + from microgen.shape.shape import BoundsType class Sphere(Shape): @@ -78,11 +79,25 @@ def generate_cad(self: Sphere, **_: KwargsGenerateType) -> CadShape: def generate_surface_mesh( self: Sphere, + bounds: BoundsType | None = None, + resolution: int | None = None, theta_resolution: int = 50, phi_resolution: int = 50, **_: KwargsGenerateType, ) -> pv.PolyData: - """Generate a sphere VTK shape using the given parameters.""" + """Generate a sphere VTK shape using the given parameters. + + When ``bounds`` or ``resolution`` is explicitly provided, fall back to + the implicit (marching-cubes-on-SDF) base implementation so polymorphic + callers that ask for a specific sampling get what they asked for. + Otherwise the native :class:`pyvista.Sphere` (parametric, controlled + by ``theta_resolution``/``phi_resolution``) is used. + """ + if bounds is not None or resolution is not None: + return super().generate_surface_mesh( + bounds=bounds, + resolution=resolution if resolution is not None else 50, + ) return pv.Sphere( radius=self.radius, center=tuple(self.center), diff --git a/microgen/shape/spinodoid.py b/microgen/shape/spinodoid.py index f184f075..685a98b3 100644 --- a/microgen/shape/spinodoid.py +++ b/microgen/shape/spinodoid.py @@ -29,12 +29,15 @@ from ..operations import rotate from ._frep_grf import _FrepGRF, _normalize_cell_size, compute_threshold_for_porosity +from .periodic_shell import mesh_to_periodic_shell from .shape import Shape if TYPE_CHECKING: from microgen.cad import CadShape from microgen.shape import KwargsGenerateType, Vector3DType + from ._types import BoundsType + class Spinodoid(Shape): """Class to generate a periodic Spinodoid microstructure. @@ -180,6 +183,11 @@ def _signed_field( self._func = _signed_field lx, ly, lz = (self.cell_size * self.repeat_cell).tolist() self._bounds = (0.0, float(lx), 0.0, float(ly), 0.0, float(lz)) + # Spinodoid's field is *bit-exact* periodic on ``cell_size`` along each + # axis — every kept Fourier mode lives on the reciprocal lattice, so + # ``frep.evaluate(p + cell_size) == frep.evaluate(p)`` (see _frep_grf.py). + cs = np.asarray(self.cell_size, dtype=float) + self._period = (float(cs[0]), float(cs[1]), float(cs[2])) @cached_property def grid(self: Spinodoid) -> pv.StructuredGrid: @@ -215,9 +223,23 @@ def generate_volume_mesh( def generate_surface_mesh( self: Spinodoid, + bounds: BoundsType | None = None, + resolution: int | None = None, **_: KwargsGenerateType, ) -> pv.PolyData: - """Generate the iso-surface as a PolyData with center+rotation applied.""" + """Generate the iso-surface as a PolyData with center+rotation applied. + + The native path uses the cached structured-grid surface whose resolution + is fixed at construction time (``Spinodoid(resolution=…)``). When + ``bounds`` or ``resolution`` is explicitly provided here, fall back to + the implicit (marching-cubes-on-SDF) base implementation so polymorphic + callers can re-sample the field at an arbitrary resolution / bbox. + """ + if bounds is not None or resolution is not None: + return super().generate_surface_mesh( + bounds=bounds, + resolution=resolution if resolution is not None else 50, + ) polydata = self.surface.copy() polydata = rotate(polydata, center=(0, 0, 0), rotation=self.orientation) return polydata.translate(xyz=self.center) @@ -241,7 +263,7 @@ def generate_cad( pts = np.asarray(mesh.points, dtype=np.float64) tris = mesh.faces.reshape(-1, 4)[:, 1:].astype(np.int64) - shape = _try_make_solid(_mesh_to_periodic_shell(pts, tris, self._bounds)) + shape = _try_make_solid(mesh_to_periodic_shell(pts, tris, self._bounds)) shape = rotate(obj=shape, center=(0, 0, 0), rotation=self.orientation) shape = shape.translate(self.center) # Fallback for OCCT Volume() on solids it flags invalid (rigid transforms preserve volume). @@ -250,93 +272,6 @@ def generate_cad( return shape -def _mesh_to_periodic_shell( - points: npt.NDArray[np.float64], - triangles: npt.NDArray[np.int64], - bounds: Sequence[float], -) -> CadShape: - """Build a sewn OCCT shell from a triangulated surface inside an axis-aligned box. - - Triangles whose three vertices share the mesh's exact extremum on a cube - plane are grouped per face and converted to a single planar BREP face via - :func:`microgen.cad.mesh_to_planar_face`. Remaining triangles become one - planar BREP face per triangle (raw, not pre-sewn — sewing must stitch - them to cap wires along the seam). Everything is sewn into a closed shell. - - Mirrors :meth:`microgen.Tpms._mesh_to_periodic_shell` (tpms.py:683) but - parameterised by ``bounds`` rather than an instance's - ``cell_size × repeat_cell`` (Spinodoid lives at ``[0, L]^3``, Tpms at - ``[-L/2, +L/2]^3``). - """ - from OCP.BRepBuilderAPI import ( - BRepBuilderAPI_MakeEdge, - BRepBuilderAPI_MakeFace, - BRepBuilderAPI_MakeWire, - BRepBuilderAPI_Sewing, - ) - from OCP.gp import gp_Pnt - from OCP.TopAbs import TopAbs_FACE - from OCP.TopExp import TopExp_Explorer - - from microgen.cad import CadShape as _CadShape - from microgen.cad import mesh_to_planar_face - - pts = np.asarray(points, dtype=np.float64) - tris = np.asarray(triangles, dtype=np.int64).reshape(-1, 3) - - drift_tol = 1e-9 * float(max(abs(b) for b in bounds) or 1.0) - - consumed = np.zeros(tris.shape[0], dtype=bool) - on_plane: list[tuple[int, int, float, npt.NDArray[np.int64]]] = [] - for axis in range(3): - for sign in (-1, +1): - extremum = ( - float(pts[:, axis].max()) if sign > 0 else float(pts[:, axis].min()) - ) - expected = bounds[2 * axis + (1 if sign > 0 else 0)] - if abs(extremum - expected) > drift_tol: - continue - vert_on = pts[:, axis] == extremum - tri_on = np.all(vert_on[tris], axis=1) & ~consumed - if tri_on.any(): - on_plane.append((axis, sign, extremum, np.where(tri_on)[0])) - consumed |= tri_on - interior_idx = np.where(~consumed)[0] - - bbox_diag = float(np.linalg.norm(pts.max(axis=0) - pts.min(axis=0))) - sew_tol = max(1e-9, 1e-6 * bbox_diag) - sewing = BRepBuilderAPI_Sewing(sew_tol) - - for axis, sign, extremum, tri_idx in on_plane: - origin = [0.0, 0.0, 0.0] - origin[axis] = extremum - normal = [0.0, 0.0, 0.0] - normal[axis] = float(sign) - planar = mesh_to_planar_face(pts, tris[tri_idx], origin, normal) - exp = TopExp_Explorer(planar.wrapped, TopAbs_FACE) - while exp.More(): - sewing.Add(exp.Current()) - exp.Next() - - interior_tris = tris[interior_idx] - used_vertices = np.unique(interior_tris) - pnt_cache = { - int(i): gp_Pnt(float(pts[i, 0]), float(pts[i, 1]), float(pts[i, 2])) - for i in used_vertices - } - for a, b, c in interior_tris: - ga, gb, gc = pnt_cache[int(a)], pnt_cache[int(b)], pnt_cache[int(c)] - e1 = BRepBuilderAPI_MakeEdge(ga, gb).Edge() - e2 = BRepBuilderAPI_MakeEdge(gb, gc).Edge() - e3 = BRepBuilderAPI_MakeEdge(gc, ga).Edge() - wire = BRepBuilderAPI_MakeWire(e1, e2, e3).Wire() - face = BRepBuilderAPI_MakeFace(wire).Face() - sewing.Add(face) - - sewing.Perform() - return _CadShape(sewing.SewedShape()) - - def _try_make_solid(shape: CadShape) -> CadShape: """Best-effort upgrade of a sewn shell (or compound of shells) to a Solid. diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 66a37723..8d1f77d5 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -26,7 +26,8 @@ from microgen.operations import fuse_shapes, rotate -from .shape import BoundsType, Shape, ShellCreationError +from ._types import BoundsType, Field +from .shape import Shape, ShellCreationError if TYPE_CHECKING: from microgen.cad import CadShape @@ -36,10 +37,6 @@ from .tpms_grading import OffsetGrading logging.basicConfig(level=logging.INFO) -Field = Callable[ - [npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]], - npt.NDArray[np.float64], -] _DIM = 3 @@ -480,13 +477,18 @@ def _finalize_frep( raw_field: Field, bounds: tuple[float, float, float, float, float, float], ) -> None: - """Normalize a raw field to SDF and set ``_func`` / ``_bounds``.""" + """Normalize a raw field to SDF and set ``_func`` / ``_bounds`` / ``_period``.""" from .implicit_ops import from_field, normalize_to_sdf self._raw_field_func = raw_field sdf_shape = normalize_to_sdf(from_field(raw_field)) self._func = sdf_shape.func self._bounds = bounds + # The TPMS field is intrinsically periodic on ``cell_size`` along each + # axis (the 2π/cell_size wavenumbers in ``_setup_frep_field`` make + # this a data-structure invariant, not a flag). + cs = np.asarray(self.cell_size, dtype=float) + self._period = (float(cs[0]), float(cs[1]), float(cs[2])) def _setup_frep_field(self: Tpms) -> None: """Build the F-rep implicit field (SDF-normalized) for this TPMS.""" @@ -691,107 +693,17 @@ def offset( def _mesh_to_periodic_shell(self: Tpms, mesh: pv.PolyData) -> CadShape: """Build a closed, sewn OCCT shell with planar BREP faces on cell sides. - Splits the triangle mesh into seven groups: one per unit-cell - boundary plane (x=±half, y=±half, z=±half) plus the TPMS interior. - Each cell-side group goes through :func:`mesh_to_planar_face`, - producing one or more :class:`TopoDS_Face` whose underlying surface - is a ``Geom_Plane`` and whose wires trace the actual cell-side - outline (so STEP shows the gyroid cuts, not a bounding cube; gmsh - ``setPeriodic`` matches opposite cell sides by their plane equation). - The TPMS interior contributes one planar face per triangle. - - All faces — caps + interior — are sewn together via - :class:`BRepBuilderAPI_Sewing` so the cap/interior seam shares edges - and face orientations are reconciled. The result is a closed shell - whose ``BRepGProp::VolumeProperties`` matches the underlying VTK - volume, and which is a valid input to boolean ops. + Delegates to :func:`microgen.shape.periodic_shell.mesh_to_periodic_shell`, + which sews TPMS-interior triangles with per-face cap planes into a + single closed shell (see that function's docstring for details). """ - from OCP.BRepBuilderAPI import ( - BRepBuilderAPI_MakeEdge, - BRepBuilderAPI_MakeFace, - BRepBuilderAPI_MakeWire, - BRepBuilderAPI_Sewing, - ) - from OCP.gp import gp_Pnt - from OCP.TopAbs import TopAbs_FACE - from OCP.TopExp import TopExp_Explorer - - from microgen.cad import ( - CadShape as _CadShape, - ) - from microgen.cad import ( - mesh_to_planar_face, - ) + from .periodic_shell import mesh_to_periodic_shell # noqa: PLC0415 if not mesh.is_all_triangles: mesh.triangulate(inplace=True) pts = np.asarray(mesh.points, dtype=np.float64) tris = mesh.faces.reshape(-1, 4)[:, 1:].astype(np.int64) - - half = 0.5 * np.asarray(self.cell_size) * np.asarray(self.repeat_cell) - # A vertex is "on" a cell-side plane iff its coordinate on that axis - # equals the mesh's exact extremum on that axis. Marching cubes on - # the periodic-aligned grid places cap vertices on in-plane grid - # edges, so they share the linspace endpoint exactly — no tolerance - # ball, just exact equality on the actual data. Slanted surface - # triangles near the boundary land on perpendicular grid edges and - # never produce that exact value. - # Sanity-check the extremum is the expected ±half (skip otherwise so - # we never misclassify a surface-only mesh whose extremum is not on - # a cube face). - drift_tol = 1e-9 * float(np.max(np.abs(half)) or 1.0) - - consumed = np.zeros(tris.shape[0], dtype=bool) - on_plane: list[tuple[int, int, npt.NDArray[np.int64]]] = [] - for axis in range(3): - for sign in (-1, +1): - extremum = ( - float(pts[:, axis].max()) if sign > 0 else float(pts[:, axis].min()) - ) - if abs(extremum - sign * float(half[axis])) > drift_tol: - continue - vert_on = pts[:, axis] == extremum - tri_on = np.all(vert_on[tris], axis=1) & ~consumed - if tri_on.any(): - on_plane.append((axis, sign, np.where(tri_on)[0])) - consumed |= tri_on - interior_idx = np.where(~consumed)[0] - - # Sewing tolerance: a small fraction of the bbox diagonal absorbs - # numerical drift between cap-wire vertices and interior-triangle - # vertices that should match exactly along the cell-side seam. - bbox_diag = float(np.linalg.norm(pts.max(axis=0) - pts.min(axis=0))) - sew_tol = max(1e-9, 1e-6 * bbox_diag) - sewing = BRepBuilderAPI_Sewing(sew_tol) - - for axis, sign, tri_idx in on_plane: - origin = [0.0, 0.0, 0.0] - origin[axis] = sign * float(half[axis]) - normal = [0.0, 0.0, 0.0] - normal[axis] = float(sign) - planar = mesh_to_planar_face(pts, tris[tri_idx], origin, normal) - exp = TopExp_Explorer(planar.wrapped, TopAbs_FACE) - while exp.More(): - sewing.Add(exp.Current()) - exp.Next() - - # Interior TPMS triangles: contribute as raw per-triangle faces so - # sewing can stitch them to the cap wires (a pre-sewn interior would - # have shared edges already, blocking the seam stitch). - for a, b, c in tris[interior_idx]: - pa, pb, pc = pts[int(a)], pts[int(b)], pts[int(c)] - ga = gp_Pnt(float(pa[0]), float(pa[1]), float(pa[2])) - gb = gp_Pnt(float(pb[0]), float(pb[1]), float(pb[2])) - gc = gp_Pnt(float(pc[0]), float(pc[1]), float(pc[2])) - e1 = BRepBuilderAPI_MakeEdge(ga, gb).Edge() - e2 = BRepBuilderAPI_MakeEdge(gb, gc).Edge() - e3 = BRepBuilderAPI_MakeEdge(gc, ga).Edge() - wire = BRepBuilderAPI_MakeWire(e1, e2, e3).Wire() - face = BRepBuilderAPI_MakeFace(wire).Face() - sewing.Add(face) - - sewing.Perform() - return _CadShape(sewing.SewedShape()) + return mesh_to_periodic_shell(pts, tris, self._bounds) def _mesh_to_shell(self: Tpms, mesh: pv.PolyData) -> CadShape: """Convert a triangulated PyVista mesh to an OCCT ``CadShape``. diff --git a/tests/shapes/test_shape_period.py b/tests/shapes/test_shape_period.py new file mode 100644 index 00000000..72ad52ef --- /dev/null +++ b/tests/shapes/test_shape_period.py @@ -0,0 +1,100 @@ +"""Tests for ``Shape.period``: the intrinsic-periodicity attribute. + +When set, ``shape.period == (Lx, Ly, Lz)`` is a data-structure invariant +guaranteeing ``shape.evaluate(p + L) == shape.evaluate(p)`` along each axis. +``Tpms`` / ``Spinodoid`` populate it from ``cell_size * repeat_cell`` / +``cell_size`` respectively; free shapes built via ``from_field`` or boolean +composition leave it ``None``. +""" + +# ruff: noqa: S101 + +from __future__ import annotations + +import numpy as np + +from microgen import Box, Sphere, Spinodoid, Tpms, surface_functions +from microgen.shape.implicit_ops import from_field +from microgen.shape.shape import Shape + + +def test_bare_shape_period_is_none() -> None: + """A free ``Shape(func=...)`` has no intrinsic period.""" + s = Shape(func=lambda x, y, z: x * x + y * y + z * z - 1.0) + assert s.period is None + + +def test_box_sphere_have_no_period() -> None: + """Non-periodic primitives expose ``period is None``.""" + assert Box().period is None + assert Sphere().period is None + + +def test_tpms_period_matches_cell_size() -> None: + """``Tpms.period`` equals ``cell_size`` (per-axis).""" + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.5, cell_size=2.0) + assert tpms.period == (2.0, 2.0, 2.0) + + +def test_tpms_anisotropic_cell_size_period() -> None: + """``Tpms.period`` is per-axis when ``cell_size`` is a tuple.""" + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.5, + cell_size=(1.0, 2.0, 3.0), + ) + assert tpms.period == (1.0, 2.0, 3.0) + + +def test_spinodoid_period_matches_cell_size() -> None: + """``Spinodoid.period`` equals ``cell_size``.""" + sp = Spinodoid(offset=0.0, cell_size=1.5, resolution=8, seed=42) + assert sp.period == (1.5, 1.5, 1.5) + + +def test_from_field_propagates_no_period() -> None: + """``from_field`` produces a free shape with ``period is None``.""" + s = from_field(lambda x, y, z: x + y + z) + assert s.period is None + + +def test_boolean_composition_does_not_inherit_period() -> None: + """``a | b`` produces a free shape; intrinsic periodicity is not preserved + by composition (the union of two periodic fields may or may not be + periodic — we don't claim it is). + """ + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.5) + sphere = Sphere(radius=0.5) + combined = tpms | sphere + assert combined.period is None + + +def test_tpms_field_actually_periodic_at_declared_period() -> None: + """If ``shape.period == (Lx, Ly, Lz)`` then ``f(p + L) == f(p)``.""" + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.5, cell_size=1.0) + lx, ly, lz = tpms.period + rng = np.random.default_rng(seed=0) + pts = rng.uniform(-0.5, 0.5, size=(20, 3)) + x, y, z = pts[:, 0], pts[:, 1], pts[:, 2] + f0 = tpms.evaluate(x, y, z) + fp = tpms.evaluate(x + lx, y + ly, z + lz) + # TPMS field is approximately periodic after SDF normalization (not bit-exact); + # tight but not exact tolerance. + assert np.allclose(f0, fp, atol=1e-6) + + +def test_spinodoid_field_periodic_at_declared_period() -> None: + """Spinodoid's reciprocal-lattice modes give (float-precision) periodicity. + + Bit-exactness is asserted by ``test_spinodoid_field_is_bit_exact_periodic`` + in tests/test_spinodoid.py on the raw ``_frep.evaluate``; via ``Shape.evaluate`` + the sign-flip + shifted-coord path picks up ULP-level drift (~1e-15). + """ + sp = Spinodoid(offset=0.0, cell_size=1.0, resolution=8, seed=42) + lx, ly, lz = sp.period + rng = np.random.default_rng(seed=1) + pts = rng.uniform(0.0, 1.0, size=(20, 3)) + x, y, z = pts[:, 0], pts[:, 1], pts[:, 2] + f0 = sp.evaluate(x, y, z) + fp = sp.evaluate(x + lx, y + ly, z + lz) + assert np.allclose(f0, fp, atol=1e-12, rtol=0) diff --git a/tests/test_no_top_level_ocp_imports.py b/tests/test_no_top_level_ocp_imports.py new file mode 100644 index 00000000..bc060c0d --- /dev/null +++ b/tests/test_no_top_level_ocp_imports.py @@ -0,0 +1,86 @@ +"""Lint-style test: no top-level OCP imports outside the CAD boundary. + +The ``[cad]`` extra (``cadquery-ocp-novtk``) is strictly optional. ``import +microgen`` must succeed without it, which means no module reachable from +``microgen.__init__`` may execute ``import OCP`` / ``from OCP …`` at +module-load time. Every OCP import must be either: + + * inside a function/method body (lazy); + * guarded by ``if TYPE_CHECKING:`` (only consumed by type checkers); or + * located inside the CAD boundary (``microgen/cad.py`` today; the + ``microgen/cad/`` subpackage after PR 4). + +This test walks every ``.py`` under ``microgen/`` and AST-checks the +top level for stray OCP imports. Failing here is what catches a future +refactor that accidentally drops an eager ``from OCP.… import …`` at +module scope and silently breaks ``pip install microgen`` (no extras). +""" +# ruff: noqa: S101 + +from __future__ import annotations + +import ast +from pathlib import Path + +MICROGEN_ROOT = Path(__file__).resolve().parent.parent / "microgen" + +# The CAD boundary: every file under ``microgen/cad/`` is allowed to have +# top-level OCP imports. Everything outside it must keep OCP imports lazy +# (inside a function body) or guarded by ``if TYPE_CHECKING:``. +_CAD_BOUNDARY_PREFIX = "cad/" + + +def _is_type_checking_block(node: ast.stmt) -> bool: + """True if ``node`` is ``if TYPE_CHECKING:`` (any common spelling).""" + if not isinstance(node, ast.If): + return False + test = node.test + if isinstance(test, ast.Name) and test.id == "TYPE_CHECKING": + return True + if ( + isinstance(test, ast.Attribute) + and isinstance(test.value, ast.Name) + and test.value.id == "typing" + and test.attr == "TYPE_CHECKING" + ): + return True + return False + + +def _top_level_imports(tree: ast.Module) -> list[ast.stmt]: + """Return module-level import statements, skipping ``if TYPE_CHECKING:``.""" + out: list[ast.stmt] = [] + for node in tree.body: + if isinstance(node, (ast.Import, ast.ImportFrom)): + out.append(node) + elif _is_type_checking_block(node): + # Skip — these imports never execute at runtime. + continue + return out + + +def _imports_ocp(node: ast.stmt) -> bool: + if isinstance(node, ast.ImportFrom): + return (node.module or "").split(".", 1)[0] == "OCP" + if isinstance(node, ast.Import): + return any(alias.name.split(".", 1)[0] == "OCP" for alias in node.names) + return False + + +def test_no_top_level_ocp_imports_outside_cad_boundary() -> None: + """Walk microgen/*.py and assert no eager OCP import escapes the CAD boundary.""" + offenders: list[str] = [] + for py in MICROGEN_ROOT.rglob("*.py"): + rel = py.relative_to(MICROGEN_ROOT).as_posix() + if rel.startswith(_CAD_BOUNDARY_PREFIX): + continue + tree = ast.parse(py.read_text(encoding="utf-8"), filename=str(py)) + for node in _top_level_imports(tree): + if _imports_ocp(node): + offenders.append(f"{rel}:{node.lineno}") + assert not offenders, ( + "Top-level OCP imports outside the CAD boundary " + f"(microgen/{_CAD_BOUNDARY_PREFIX}*): {offenders}. " + "Move them inside a function body or under " + "``if TYPE_CHECKING:``." + ) diff --git a/tests/test_rve.py b/tests/test_rve.py index 63dbdbf6..21397342 100644 --- a/tests/test_rve.py +++ b/tests/test_rve.py @@ -1,7 +1,10 @@ """Tests for Rve class.""" +import dataclasses + import numpy as np import pytest +import pyvista as pv from microgen import Rve @@ -80,3 +83,84 @@ def test_rve_from_min_max_must_return_expected_rve() -> None: rve = Rve.from_min_max(min_point=(-1, 0, -3), max_point=(0, 2, 3)) assert np.all(rve.center == [-0.5, 1, 0]) assert np.all(rve.dim == [1, 2, 6]) + + +def test_rve_pbc_default_is_fully_periodic() -> None: + """``pbc`` defaults to ``(True, True, True)``.""" + assert Rve().pbc == (True, True, True) + + +def test_rve_pbc_scalar_broadcasts_to_all_axes() -> None: + """A single bool is broadcast to all three axes.""" + assert Rve(pbc=False).pbc == (False, False, False) + assert Rve(pbc=True).pbc == (True, True, True) + + +def test_rve_pbc_accepts_per_axis_tuple() -> None: + """A length-3 sequence sets each axis independently.""" + assert Rve(pbc=(True, False, True)).pbc == (True, False, True) + assert Rve(pbc=[False, False, True]).pbc == (False, False, True) + + +def test_rve_invalid_pbc_raises() -> None: + """``pbc`` of wrong length or type raises ``ValueError``.""" + with pytest.raises(ValueError, match="pbc must be a bool"): + Rve(pbc=(True, False)) # type: ignore[arg-type] + with pytest.raises(ValueError, match="pbc must be a bool"): + Rve(pbc=1.5) # type: ignore[arg-type] + + +def test_rve_is_frozen() -> None: + """Reassigning ``center``/``dim``/``pbc`` after construction raises FrozenInstanceError.""" + rve = Rve(dim=2) + with pytest.raises(dataclasses.FrozenInstanceError): + rve.center = np.array([1.0, 2.0, 3.0]) # type: ignore[misc] + with pytest.raises(dataclasses.FrozenInstanceError): + rve.dim = np.array([3.0, 3.0, 3.0]) # type: ignore[misc] + with pytest.raises(dataclasses.FrozenInstanceError): + rve.pbc = (False, False, False) # type: ignore[misc] + + +def test_rve_min_max_points_are_consistent() -> None: + """``min_point`` / ``max_point`` are derived from center ± 0.5·dim.""" + rve = Rve(center=(1.0, 2.0, 3.0), dim=(2.0, 4.0, 6.0)) + assert np.allclose(rve.min_point, [0.0, 0.0, 0.0]) + assert np.allclose(rve.max_point, [2.0, 4.0, 6.0]) + + +def test_rve_grid_with_scalar_resolution() -> None: + """``Rve.grid(n)`` returns a StructuredGrid spanning the cell with n^3 points.""" + rve = Rve(center=(0.0, 0.0, 0.0), dim=(2.0, 4.0, 6.0)) + grid = rve.grid(5) + assert isinstance(grid, pv.StructuredGrid) + assert grid.dimensions == (5, 5, 5) + pts = np.asarray(grid.points) + assert np.isclose(pts[:, 0].min(), -1.0) + assert np.isclose(pts[:, 0].max(), 1.0) + assert np.isclose(pts[:, 1].min(), -2.0) + assert np.isclose(pts[:, 1].max(), 2.0) + assert np.isclose(pts[:, 2].min(), -3.0) + assert np.isclose(pts[:, 2].max(), 3.0) + + +def test_rve_grid_with_per_axis_resolution() -> None: + """``Rve.grid((nx, ny, nz))`` honors per-axis resolution.""" + rve = Rve(center=(0.0, 0.0, 0.0), dim=1.0) + grid = rve.grid((3, 4, 5)) + assert grid.dimensions == (3, 4, 5) + + +def test_rve_grid_invalid_resolution_raises() -> None: + """Wrong-length resolution raises ``ValueError``.""" + rve = Rve() + with pytest.raises(ValueError, match="resolution must be an int"): + rve.grid((3, 4)) # type: ignore[arg-type] + + +def test_rve_repr_round_trips_inputs() -> None: + """``repr(rve)`` contains the canonical ``center``/``dim``/``pbc`` triple.""" + rve = Rve(center=(1.0, 2.0, 3.0), dim=(2.0, 2.0, 2.0), pbc=(True, False, True)) + text = repr(rve) + assert "center=(1.0, 2.0, 3.0)" in text + assert "dim=(2.0, 2.0, 2.0)" in text + assert "pbc=(True, False, True)" in text