From 60939e8c31392ada6a812f112a92183baa0fb848 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Thu, 16 Apr 2026 15:53:14 +0200 Subject: [PATCH 01/11] Add implicit F-rep ops and Shape support Introduce F-rep implicit operations and integrate implicit-field support into Shape. Adds new module microgen.shape.implicit_ops (boolean, smooth blending, batch_smooth_union, repeat, shell, blend, from_field, etc.) and comprehensive tests. Refactors microgen/shape/shape.py to carry optional implicit fields, provide evaluation, default generate_vtk/generate via marching-cubes -> CadQuery, boolean operators, smooth booleans, and transform helpers (translate/rotate/scale). Updates package exports and docs to expose implicit_ops and related helpers; re-exports ShellCreationError from tpms for backward compatibility. --- docs/microgen.shape.rst | 5 + microgen/__init__.py | 6 + microgen/shape/__init__.py | 6 +- microgen/shape/implicit_ops.py | 305 ++++++++++++++++++++++ microgen/shape/shape.py | 322 +++++++++++++++++++++-- microgen/shape/tpms.py | 8 +- tests/shapes/test_implicit_ops.py | 420 ++++++++++++++++++++++++++++++ 7 files changed, 1045 insertions(+), 27 deletions(-) create mode 100644 microgen/shape/implicit_ops.py create mode 100644 tests/shapes/test_implicit_ops.py diff --git a/docs/microgen.shape.rst b/docs/microgen.shape.rst index 0a0cc052..9d0c532c 100644 --- a/docs/microgen.shape.rst +++ b/docs/microgen.shape.rst @@ -52,6 +52,11 @@ microgen.shape package :undoc-members: :show-inheritance: +.. automodule:: microgen.shape.implicit_ops + :members: + :undoc-members: + :show-inheritance: + .. automodule:: microgen.shape.surface_functions :members: :undoc-members: diff --git a/microgen/__init__.py b/microgen/__init__.py index e5f66939..ff5e9d7d 100644 --- a/microgen/__init__.py +++ b/microgen/__init__.py @@ -60,6 +60,9 @@ TruncatedCube, TruncatedCuboctahedron, TruncatedOctahedron, + batch_smooth_union, + from_field, + implicit_ops, new_geometry, newGeometry, surface_functions, @@ -100,6 +103,9 @@ "TruncatedOctahedron", "CylindricalTpms", "SphericalTpms", + "batch_smooth_union", + "from_field", + "implicit_ops", "check_if_only_linear_tetrahedral", "cut_phase_by_shape_list", "cutPhaseByShapeList", diff --git a/microgen/shape/__init__.py b/microgen/shape/__init__.py index d7fc11c3..ff74044c 100644 --- a/microgen/shape/__init__.py +++ b/microgen/shape/__init__.py @@ -21,13 +21,14 @@ from typing import TYPE_CHECKING, Callable, Literal, Sequence, Tuple -from . import surface_functions +from . import implicit_ops, surface_functions from .box import Box from .capsule import Capsule from .cylinder import Cylinder from .ellipsoid import Ellipsoid from .extruded_polygon import ExtrudedPolygon from .polyhedron import Polyhedron +from .implicit_ops import batch_smooth_union, from_field from .shape import Shape from .sphere import Sphere from .strut_lattice import ( @@ -161,6 +162,9 @@ def __init__(self: ShapeError, shape: str) -> None: "ExtrudedPolygon", "FaceCenteredCubic", "Infill", + "batch_smooth_union", + "from_field", + "implicit_ops", "NormedDistance", "Octahedron", "OctetTruss", diff --git a/microgen/shape/implicit_ops.py b/microgen/shape/implicit_ops.py new file mode 100644 index 00000000..ec9a3c18 --- /dev/null +++ b/microgen/shape/implicit_ops.py @@ -0,0 +1,305 @@ +"""F-rep Implicit Operations. + +========================================================== +Implicit Operations (:mod:`microgen.shape.implicit_ops`) +========================================================== + +Module-level boolean, blending, and utility operations for shapes +that carry an implicit scalar field (``_func``). All functions accept +and return :class:`~microgen.shape.shape.Shape` instances. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +import numpy.typing as npt + +if TYPE_CHECKING: + from .shape import BoundsType, Field, Shape + + +# --------------------------------------------------------------------------- +# Internal helpers +# --------------------------------------------------------------------------- + + +def _make_shape(func: Field, bounds: BoundsType | None) -> Shape: + """Create a Shape with an implicit field (single deferred import).""" + from .shape import Shape # noqa: PLC0415 + + return Shape(func=func, bounds=bounds) + + +def _smooth_min( + a: npt.NDArray[np.float64], + b: npt.NDArray[np.float64], + k: float, +) -> npt.NDArray[np.float64]: + """Smooth minimum (Inigo Quilez cubic polynomial).""" + if k <= 0: + return np.minimum(a, b) + h = np.maximum(k - np.abs(a - b), 0.0) / k + return np.minimum(a, b) - h * h * h * k / 6.0 + + +def _smooth_max( + a: npt.NDArray[np.float64], + b: npt.NDArray[np.float64], + k: float, +) -> npt.NDArray[np.float64]: + """Smooth maximum.""" + return -_smooth_min(-a, -b, k) + + +def _merge_bounds( + a: BoundsType | None, + b: BoundsType | None, + mode: str = "union", +) -> BoundsType | None: + """Merge two bounding boxes.""" + if a is None and b is None: + return None + if a is None: + return b + if b is None: + return a + if mode == "union": + return ( + min(a[0], b[0]), + max(a[1], b[1]), + min(a[2], b[2]), + max(a[3], b[3]), + min(a[4], b[4]), + max(a[5], b[5]), + ) + # intersection + return ( + max(a[0], b[0]), + min(a[1], b[1]), + max(a[2], b[2]), + min(a[3], b[3]), + max(a[4], b[4]), + min(a[5], b[5]), + ) + + +# --------------------------------------------------------------------------- +# Unary operations +# --------------------------------------------------------------------------- + + +def complement(a: Shape) -> Shape: + """Complement (negate the field): inside becomes outside and vice versa.""" + f = a.require_func() + return _make_shape( + func=lambda x, y, z, _f=f: -_f(x, y, z), + bounds=a.bounds, + ) + + +# --------------------------------------------------------------------------- +# Hard boolean operations +# --------------------------------------------------------------------------- + + +def union(a: Shape, b: Shape) -> Shape: + """Union of two shapes (hard boolean).""" + fa, fb = a.require_func(), b.require_func() + return _make_shape( + func=lambda x, y, z, _fa=fa, _fb=fb: np.minimum(_fa(x, y, z), _fb(x, y, z)), + bounds=_merge_bounds(a.bounds, b.bounds, "union"), + ) + + +def intersection(a: Shape, b: Shape) -> Shape: + """Intersection of two shapes (hard boolean).""" + fa, fb = a.require_func(), b.require_func() + return _make_shape( + func=lambda x, y, z, _fa=fa, _fb=fb: np.maximum(_fa(x, y, z), _fb(x, y, z)), + bounds=_merge_bounds(a.bounds, b.bounds, "intersection"), + ) + + +def difference(a: Shape, b: Shape) -> Shape: + """Difference of two shapes (a minus b).""" + fa, fb = a.require_func(), b.require_func() + return _make_shape( + func=lambda x, y, z, _fa=fa, _fb=fb: np.maximum( + _fa(x, y, z), -_fb(x, y, z), + ), + bounds=a.bounds, + ) + + +# --------------------------------------------------------------------------- +# Smooth boolean operations +# --------------------------------------------------------------------------- + + +def smooth_union(a: Shape, b: Shape, k: float) -> Shape: + """Smooth union with blending radius *k*.""" + fa, fb = a.require_func(), b.require_func() + return _make_shape( + func=lambda x, y, z, _fa=fa, _fb=fb, _k=k: _smooth_min( + _fa(x, y, z), _fb(x, y, z), _k, + ), + bounds=_merge_bounds(a.bounds, b.bounds, "union"), + ) + + +def smooth_intersection(a: Shape, b: Shape, k: float) -> Shape: + """Smooth intersection with blending radius *k*.""" + fa, fb = a.require_func(), b.require_func() + return _make_shape( + func=lambda x, y, z, _fa=fa, _fb=fb, _k=k: _smooth_max( + _fa(x, y, z), _fb(x, y, z), _k, + ), + bounds=_merge_bounds(a.bounds, b.bounds, "intersection"), + ) + + +def smooth_difference(a: Shape, b: Shape, k: float) -> Shape: + """Smooth difference (a minus b) with blending radius *k*.""" + fa, fb = a.require_func(), b.require_func() + return _make_shape( + func=lambda x, y, z, _fa=fa, _fb=fb, _k=k: _smooth_max( + _fa(x, y, z), -_fb(x, y, z), _k, + ), + bounds=a.bounds, + ) + + +# --------------------------------------------------------------------------- +# Batch operations +# --------------------------------------------------------------------------- + + +def batch_smooth_union( + shapes: list[Shape], + k: float = 0.0, +) -> Shape: + """Combine many shapes with smooth union in a flat loop (no recursion). + + This avoids the recursion-depth limit that arises when chaining hundreds + of binary ``smooth_union`` calls, each wrapping the previous in a lambda. + """ + if not shapes: + msg = "batch_smooth_union requires at least one shape" + raise ValueError(msg) + + funcs = [s.require_func() for s in shapes] + + merged = shapes[0].bounds + for s in shapes[1:]: + merged = _merge_bounds(merged, s.bounds, "union") + + def _batched( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + if k > 0: + result = funcs[0](x, y, z) + for fn in funcs[1:]: + result = _smooth_min(result, fn(x, y, z), k) + return result + # Hard union: vectorized reduction + all_fields = np.stack([fn(x, y, z) for fn in funcs], axis=0) + return np.min(all_fields, axis=0) + + return _make_shape(func=_batched, bounds=merged) + + +# --------------------------------------------------------------------------- +# Utility operations +# --------------------------------------------------------------------------- + + +def shell(shape: Shape, thickness: float) -> Shape: + """Hollow shell: ``|f(p)| - thickness / 2``.""" + f = shape.require_func() + half_t = thickness / 2.0 + return _make_shape( + func=lambda x, y, z, _f=f, _ht=half_t: np.abs(_f(x, y, z)) - _ht, + bounds=shape.bounds, + ) + + +def repeat( + shape: Shape, + spacing: tuple[float, float, float], + k: float = 0.0, +) -> Shape: + """Infinite repetition via coordinate modulo. + + :param shape: unit cell shape to tile + :param spacing: ``(sx, sy, sz)`` repetition period per axis + :param k: smooth blending radius across cell boundaries. + When ``k > 0``, the base field is evaluated at the 26 neighboring + periodic images in addition to the current cell and all values + are combined with smooth minimum, so that primitives from adjacent + cells blend seamlessly. When ``k <= 0`` (default), a simple + coordinate-modulo repetition is used (hard tiling). + """ + sx, sy, sz = spacing + f = shape.require_func() + + if k <= 0: + + def _repeated( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + rx = np.mod(x + sx / 2, sx) - sx / 2 + ry = np.mod(y + sy / 2, sy) - sy / 2 + rz = np.mod(z + sz / 2, sz) - sz / 2 + return f(rx, ry, rz) + + else: + offsets = [ + (dx * sx, dy * sy, dz * sz) + for dx in (-1, 0, 1) + for dy in (-1, 0, 1) + for dz in (-1, 0, 1) + ] + + def _repeated( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + cx = np.mod(x + sx / 2, sx) - sx / 2 + cy = np.mod(y + sy / 2, sy) - sy / 2 + cz = np.mod(z + sz / 2, sz) - sz / 2 + result = f(cx + offsets[0][0], cy + offsets[0][1], cz + offsets[0][2]) + for ox, oy, oz in offsets[1:]: + result = _smooth_min(result, f(cx + ox, cy + oy, cz + oz), k) + return result + + return _make_shape(func=_repeated, bounds=None) + + +def blend( + a: Shape, + b: Shape, + factor: float = 0.5, +) -> Shape: + """Linear interpolation between two fields: ``(1-t)*a + t*b``.""" + fa, fb = a.require_func(), b.require_func() + t = factor + return _make_shape( + func=lambda x, y, z, _fa=fa, _fb=fb, _t=t: (1.0 - _t) * _fa(x, y, z) + + _t * _fb(x, y, z), + bounds=_merge_bounds(a.bounds, b.bounds, "union"), + ) + + +def from_field( + func: Field, + bounds: BoundsType | None = None, +) -> Shape: + """Wrap any callable ``f(x, y, z) -> scalar`` as a Shape with an implicit field.""" + return _make_shape(func=func, bounds=bounds) diff --git a/microgen/shape/shape.py b/microgen/shape/shape.py index 11602962..24b98fc9 100644 --- a/microgen/shape/shape.py +++ b/microgen/shape/shape.py @@ -7,30 +7,59 @@ from __future__ import annotations -from abc import ABC, abstractmethod +import itertools +from collections.abc import Callable from typing import TYPE_CHECKING +import numpy as np +import numpy.typing as npt +import pyvista as pv from scipy.spatial.transform import Rotation +from microgen.operations import rotate as rotate_mesh + if TYPE_CHECKING: import cadquery as cq - import pyvista as pv 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] + + +class ShellCreationError(Exception): + """Raised when a CadQuery shell cannot be created from a mesh.""" -class Shape(ABC): - """Shape class to manage shapes. - :param shape: name of the shape - :param center: center - :param orientation: orientation +class Shape: + """Unified shape with optional implicit (F-rep) and CAD representations. + + Every shape has a ``center`` and ``orientation``. It may also carry an + implicit scalar field (``_func``) where ``f(x, y, z) < 0`` means *inside*. + When the implicit field is present, the default :meth:`generate_vtk` and + :meth:`generate` produce geometry via marching cubes. Subclasses + (e.g. ``Sphere``, ``Tpms``) override these methods with their own + implementations. + + Boolean operators (``|``, ``&``, ``-``, ``~``) and smooth boolean + methods operate on the implicit field and return a new :class:`Shape`. + + :param center: center of the 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`` """ def __init__( self: Shape, center: Vector3DType = (0, 0, 0), orientation: Vector3DType | Rotation = (0, 0, 0), + func: Field | None = None, + bounds: BoundsType | None = None, ) -> None: """Initialize the shape.""" self.center = center @@ -39,24 +68,275 @@ def __init__( if isinstance(orientation, Rotation) else Rotation.from_euler("ZXZ", orientation, degrees=True) ) + self._func = func + self._bounds = bounds + + # ------------------------------------------------------------------ + # Public read-only accessors for implicit field + # ------------------------------------------------------------------ - @abstractmethod - def generate(self: Shape, **_: KwargsGenerateType) -> cq.Shape: - """Generate the CAD shape. + @property + def func(self: Shape) -> Field | None: + """The implicit scalar field, or ``None``.""" + return self._func - :return: cq.Shape + @property + def bounds(self: Shape) -> BoundsType | None: + """The bounding box ``(xmin, xmax, ymin, ymax, zmin, zmax)``, or ``None``.""" + return self._bounds + + def require_func(self: Shape) -> Field: + """Return ``_func`` or raise if not set.""" + if self._func is None: + err_msg = "No implicit scalar field defined on this shape" + raise ValueError(err_msg) + return self._func + + # ------------------------------------------------------------------ + # Implicit field evaluation + # ------------------------------------------------------------------ + + def evaluate( + self: Shape, + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + """Evaluate the implicit scalar field at the given coordinates. + + :param x: x coordinates + :param y: y coordinates + :param z: z coordinates + :return: scalar field values (negative = inside) """ - raise NotImplementedError + return self.require_func()(x, y, z) + + # ------------------------------------------------------------------ + # Mesh generation (defaults use the implicit field) + # ------------------------------------------------------------------ - @abstractmethod - def generate_vtk(self: Shape, **_: KwargsGenerateType) -> pv.PolyData: - """Generate the vtk mesh of the shape. + def generate_vtk( + self: Shape, + bounds: BoundsType | None = None, + resolution: int = 50, + **_: KwargsGenerateType, + ) -> pv.PolyData: + """Generate a VTK mesh of the shape. - :return: pv.PolyData + The default implementation meshes the implicit field via marching cubes + (``f < 0`` convention). Subclasses override this with their own + geometry generation. + + :param bounds: ``(xmin, xmax, ymin, ymax, zmin, zmax)`` + :param resolution: number of grid points per axis + :return: triangulated surface mesh + """ + if self._func is None: + err_msg = ( + "No implicit field defined — " + "subclasses must override generate_vtk()" + ) + raise NotImplementedError(err_msg) + + bounds = bounds or self._bounds + if bounds is None: + err_msg = ( + "Bounds must be provided either at construction or in generate_vtk()" + ) + raise ValueError(err_msg) + + xmin, xmax, ymin, ymax, zmin, zmax = bounds + xi = np.linspace(xmin, xmax, resolution) + yi = np.linspace(ymin, ymax, resolution) + zi = np.linspace(zmin, zmax, resolution) + x, y, z = np.meshgrid(xi, yi, zi, indexing="ij") + + grid = pv.StructuredGrid(x, y, z) + field = self.evaluate( + x.ravel(order="F"), + y.ravel(order="F"), + z.ravel(order="F"), + ) + grid["implicit"] = field + + polydata = grid.contour(isosurfaces=[0.0], scalars="implicit") + if polydata.n_cells == 0: + return pv.PolyData() + + polydata = polydata.clean().triangulate() + + polydata = rotate_mesh(polydata, center=(0, 0, 0), rotation=self.orientation) + return polydata.translate(xyz=self.center) + + def generate( + self: Shape, + bounds: BoundsType | None = None, + resolution: int = 50, + **_: KwargsGenerateType, + ) -> cq.Shape: + """Generate a CAD shape. + + The default implementation builds a CadQuery shape from the + implicit-field VTK mesh. Subclasses override this with native + CAD construction. + + :param bounds: ``(xmin, xmax, ymin, ymax, zmin, zmax)`` + :param resolution: number of grid points per axis + :return: CadQuery Shape """ - raise NotImplementedError + if self._func is None: + err_msg = ( + "No implicit field defined — " + "subclasses must override generate()" + ) + raise NotImplementedError(err_msg) + + import cadquery as cq # noqa: PLC0415 - @abstractmethod - def generateVtk(self: Shape, **_: KwargsGenerateType) -> pv.PolyData: # noqa: N802 - """Deprecated. Use generate_vtk instead.""" # noqa: D401 - return self.generate_vtk() + mesh = self.generate_vtk(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:] + triangles = np.c_[triangles, triangles[:, 0]] + + faces = [] + for tri in triangles: + lines = [ + cq.Edge.makeLine( + cq.Vector(*mesh.points[start]), + cq.Vector(*mesh.points[end]), + ) + for start, end in itertools.pairwise(tri) + ] + wire = cq.Wire.assembleEdges(lines) + faces.append(cq.Face.makeFromWires(wire)) + + try: + shell = cq.Shell.makeShell(faces) + except ValueError as err: + err_msg = ( + "Failed to create the shell, " + "try to increase the resolution or adjust bounds." + ) + raise ShellCreationError(err_msg) from err + + return cq.Shape(shell.wrapped) + + def generateVtk(self: Shape, **kwargs: KwargsGenerateType) -> pv.PolyData: # noqa: N802 + """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + return self.generate_vtk(**kwargs) + + # ------------------------------------------------------------------ + # Boolean operators (on implicit field) + # ------------------------------------------------------------------ + + def __or__(self: Shape, other: Shape) -> Shape: + """Union (``a | b``): inside where either field is negative.""" + from .implicit_ops import union # noqa: PLC0415 + + return union(self, other) + + def __and__(self: Shape, other: Shape) -> Shape: + """Intersection (``a & b``): inside where both fields are negative.""" + from .implicit_ops import intersection # noqa: PLC0415 + + return intersection(self, other) + + def __sub__(self: Shape, other: Shape) -> Shape: + """Difference (``a - b``): inside *a* but not *b*.""" + from .implicit_ops import difference # noqa: PLC0415 + + return difference(self, other) + + def __invert__(self: Shape) -> Shape: + """Complement (``~a``): negate the field.""" + from .implicit_ops import complement # noqa: PLC0415 + + return complement(self) + + # ------------------------------------------------------------------ + # Smooth booleans + # ------------------------------------------------------------------ + + def smooth_union(self: Shape, other: Shape, k: float) -> Shape: + """Smooth union with blending radius *k*.""" + from .implicit_ops import smooth_union # noqa: PLC0415 + + return smooth_union(self, other, k) + + def smooth_intersection(self: Shape, other: Shape, k: float) -> Shape: + """Smooth intersection with blending radius *k*.""" + from .implicit_ops import smooth_intersection # noqa: PLC0415 + + return smooth_intersection(self, other, k) + + def smooth_difference(self: Shape, other: Shape, k: float) -> Shape: + """Smooth difference with blending radius *k*.""" + from .implicit_ops import smooth_difference # noqa: PLC0415 + + return smooth_difference(self, other, k) + + # ------------------------------------------------------------------ + # Implicit field transforms + # ------------------------------------------------------------------ + + def translate(self: Shape, offset: tuple[float, float, float]) -> Shape: + """Return a new shape translated by *offset* (implicit field).""" + f = self.require_func() + dx, dy, dz = offset + new_bounds = None + if self._bounds is not None: + b = self._bounds + new_bounds = ( + b[0] + dx, + b[1] + dx, + b[2] + dy, + b[3] + dy, + b[4] + dz, + b[5] + dz, + ) + return Shape( + func=lambda x, y, z, _f=f, _dx=dx, _dy=dy, _dz=dz: _f( + x - _dx, y - _dy, z - _dz, + ), + bounds=new_bounds, + ) + + def rotate( + self: Shape, + angles: tuple[float, float, float], + convention: str = "ZXZ", + ) -> Shape: + """Return a new shape rotated by Euler *angles* (degrees, implicit field).""" + f = self.require_func() + rot = Rotation.from_euler(convention, angles, degrees=True) + inv_matrix = rot.inv().as_matrix() + return Shape( + func=lambda x, y, z, _f=f, _m=inv_matrix: _f( + *(_m @ np.array([x, y, z])), + ), + bounds=self._bounds, # conservative: keep original bounds + ) + + def scale(self: Shape, factor: float) -> Shape: + """Return a new shape uniformly scaled by *factor* (implicit field).""" + f = self.require_func() + new_bounds = None + if self._bounds is not None: + b = self._bounds + new_bounds = ( + b[0] * factor, + b[1] * factor, + b[2] * factor, + b[3] * factor, + b[4] * factor, + b[5] * factor, + ) + return Shape( + func=lambda x, y, z, _f=f, _s=factor: _f(x / _s, y / _s, z / _s) * _s, + bounds=new_bounds, + ) diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 0f2e2251..d238a4fd 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -958,9 +958,7 @@ def _create_grid( return grid -class ShellCreationError(Exception): - """Error raised when the shell creation fails.""" +# Re-export for backward compatibility +from .shape import ShellCreationError # noqa: E402 - def __init__(self: ShellCreationError, message: str) -> None: - """Initialize the ShellCreationError.""" - super().__init__(message) +__all__ = ["CylindricalTpms", "Infill", "ShellCreationError", "SphericalTpms", "Tpms"] diff --git a/tests/shapes/test_implicit_ops.py b/tests/shapes/test_implicit_ops.py new file mode 100644 index 00000000..03b1acde --- /dev/null +++ b/tests/shapes/test_implicit_ops.py @@ -0,0 +1,420 @@ +"""Tests for the F-rep implicit operations and Shape implicit capabilities.""" + +from __future__ import annotations + +import numpy as np +import pytest +import pyvista as pv + +from microgen.shape.implicit_ops import ( + batch_smooth_union, + blend, + difference, + from_field, + intersection, + repeat, + shell, + smooth_difference, + smooth_intersection, + smooth_union, + union, +) +from microgen.shape.shape import Shape + + +# --------------------------------------------------------------------------- +# Helpers — inline SDF lambdas (no primitive factories needed) +# --------------------------------------------------------------------------- + + +def _sphere_field(cx=0.0, cy=0.0, cz=0.0, r=1.0): + """Return a sphere SDF function and bounds.""" + margin = r * 1.1 + return ( + lambda x, y, z: np.sqrt((x - cx) ** 2 + (y - cy) ** 2 + (z - cz) ** 2) - r, + (cx - margin, cx + margin, cy - margin, cy + margin, cz - margin, cz + margin), + ) + + +def _box_field(cx=0.0, cy=0.0, cz=0.0, hx=0.5, hy=0.5, hz=0.5): + """Return a box SDF function and bounds.""" + margin = max(hx, hy, hz) * 0.1 + + def sdf(x, y, z): + qx = np.abs(x - cx) - hx + qy = np.abs(y - cy) - hy + qz = np.abs(z - cz) - hz + outside = np.sqrt( + np.maximum(qx, 0.0) ** 2 + + np.maximum(qy, 0.0) ** 2 + + np.maximum(qz, 0.0) ** 2 + ) + inside = np.minimum(np.maximum(qx, np.maximum(qy, qz)), 0.0) + return outside + inside + + return ( + sdf, + ( + cx - hx - margin, + cx + hx + margin, + cy - hy - margin, + cy + hy + margin, + cz - hz - margin, + cz + hz + margin, + ), + ) + + +def _make_sphere(cx=0.0, cy=0.0, cz=0.0, r=1.0): + func, bounds = _sphere_field(cx, cy, cz, r) + return Shape(func=func, bounds=bounds) + + +def _make_box(cx=0.0, cy=0.0, cz=0.0, hx=0.5, hy=0.5, hz=0.5): + func, bounds = _box_field(cx, cy, cz, hx, hy, hz) + return Shape(func=func, bounds=bounds) + + +# --------------------------------------------------------------------------- +# Evaluate +# --------------------------------------------------------------------------- + + +class TestEvaluate: + """Test implicit field evaluation.""" + + def test_sphere_inside(self): + s = _make_sphere() + assert s.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0]))[0] < 0 + + def test_sphere_surface(self): + s = _make_sphere() + val = s.evaluate(np.array([1.0]), np.array([0.0]), np.array([0.0]))[0] + assert abs(val) < 1e-10 + + def test_sphere_outside(self): + s = _make_sphere() + assert s.evaluate(np.array([2.0]), np.array([0.0]), np.array([0.0]))[0] > 0 + + def test_no_func_raises(self): + s = Shape() + with pytest.raises(ValueError, match="No implicit scalar field"): + s.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0])) + + +# --------------------------------------------------------------------------- +# Boolean operations +# --------------------------------------------------------------------------- + + +class TestBooleans: + """Test hard boolean operations.""" + + def test_union_min(self): + s1 = _make_sphere(cx=-0.5) + s2 = _make_sphere(cx=0.5) + u = union(s1, s2) + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + expected = min(s1.evaluate(x, y, z)[0], s2.evaluate(x, y, z)[0]) + assert u.evaluate(x, y, z)[0] == pytest.approx(expected) + + def test_intersection_max(self): + s1 = _make_sphere(cx=-0.5) + s2 = _make_sphere(cx=0.5) + i = intersection(s1, s2) + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + expected = max(s1.evaluate(x, y, z)[0], s2.evaluate(x, y, z)[0]) + assert i.evaluate(x, y, z)[0] == pytest.approx(expected) + + def test_difference(self): + s1 = _make_sphere() + s2 = _make_sphere(cx=0.5, r=0.5) + d = difference(s1, s2) + x, y, z = np.array([-0.5]), np.array([0.0]), np.array([0.0]) + assert d.evaluate(x, y, z)[0] < 0 + + def test_operators(self): + s1 = _make_sphere() + s2 = _make_box() + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + + u = s1 | s2 + assert u.evaluate(x, y, z)[0] < 0 + + i = s1 & s2 + assert i.evaluate(x, y, z)[0] < 0 + + d = s1 - s2 + val = d.evaluate(x, y, z)[0] + assert isinstance(val, (float, np.floating)) + + def test_complement(self): + s = _make_sphere() + c = ~s + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + assert c.evaluate(x, y, z)[0] > 0 + + +# --------------------------------------------------------------------------- +# Smooth booleans +# --------------------------------------------------------------------------- + + +class TestSmoothBooleans: + """Test smooth boolean operations.""" + + def test_smooth_union_blending_zone(self): + s1 = _make_sphere(cx=-0.5) + s2 = _make_sphere(cx=0.5) + su = smooth_union(s1, s2, k=0.5) + hard = union(s1, s2) + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + assert su.evaluate(x, y, z)[0] <= hard.evaluate(x, y, z)[0] + 1e-10 + + def test_smooth_union_k_zero_equals_hard(self): + s1 = _make_sphere(cx=-0.5) + s2 = _make_sphere(cx=0.5) + su = smooth_union(s1, s2, k=0.0) + hard = union(s1, s2) + x = np.linspace(-2, 2, 10) + y = np.zeros(10) + z = np.zeros(10) + np.testing.assert_allclose( + su.evaluate(x, y, z), hard.evaluate(x, y, z), atol=1e-10 + ) + + def test_smooth_intersection(self): + s1 = _make_sphere() + s2 = _make_box() + si = smooth_intersection(s1, s2, k=0.3) + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + assert si.evaluate(x, y, z)[0] < 0 + + def test_smooth_difference(self): + s1 = _make_sphere() + s2 = _make_sphere(cx=1.0, r=0.5) + sd = smooth_difference(s1, s2, k=0.2) + x, y, z = np.array([-0.5]), np.array([0.0]), np.array([0.0]) + assert sd.evaluate(x, y, z)[0] < 0 + + def test_smooth_methods(self): + s1 = _make_sphere() + s2 = _make_box() + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + + su = s1.smooth_union(s2, k=0.3) + assert isinstance(su.evaluate(x, y, z)[0], (float, np.floating)) + + si = s1.smooth_intersection(s2, k=0.3) + assert isinstance(si.evaluate(x, y, z)[0], (float, np.floating)) + + sd = s1.smooth_difference(s2, k=0.3) + assert isinstance(sd.evaluate(x, y, z)[0], (float, np.floating)) + + +# --------------------------------------------------------------------------- +# Transforms +# --------------------------------------------------------------------------- + + +class TestTransforms: + """Test implicit field transform operations.""" + + def test_translate(self): + s = _make_sphere() + st = s.translate((2, 0, 0)) + x, y, z = np.array([2.0]), np.array([0.0]), np.array([0.0]) + assert st.evaluate(x, y, z)[0] < 0 + + x0 = np.array([0.0]) + assert st.evaluate(x0, y, z)[0] > 0 + + def test_rotate_90(self): + # Elongated box along x, rotate 90 around z -> elongated along y + func, bounds = _box_field(hx=1.0, hy=0.1, hz=0.1) + box = Shape(func=func, bounds=bounds) + rotated = box.rotate((0, 0, 90), convention="xyz") + # Point along y axis should be inside + assert rotated.evaluate(np.array([0.0]), np.array([0.5]), np.array([0.0]))[0] < 0 + # Point along x axis (was inside, now outside) + assert rotated.evaluate(np.array([0.5]), np.array([0.0]), np.array([0.0]))[0] > 0 + + def test_scale(self): + s = _make_sphere() + ss = s.scale(2.0) + x, y, z = np.array([1.5]), np.array([0.0]), np.array([0.0]) + assert s.evaluate(x, y, z)[0] > 0 + assert ss.evaluate(x, y, z)[0] < 0 + + def test_translate_bounds(self): + s = _make_sphere() + st = s.translate((5, 0, 0)) + assert st.bounds is not None + assert st.bounds[0] > 3.0 + + def test_scale_bounds(self): + s = _make_sphere() + ss = s.scale(3.0) + assert ss.bounds is not None + assert ss.bounds[1] > 3.0 + + +# --------------------------------------------------------------------------- +# generate_vtk +# --------------------------------------------------------------------------- + + +class TestGenerateVtk: + """Test default mesh generation from implicit field.""" + + def test_sphere_mesh(self): + s = _make_sphere() + mesh = s.generate_vtk(resolution=30) + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + def test_box_mesh(self): + b = _make_box() + mesh = b.generate_vtk(resolution=30) + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + def test_boolean_mesh(self): + s1 = _make_sphere() + s2 = _make_box(hx=0.6, hy=0.6, hz=0.6) + result = s1 & s2 + mesh = result.generate_vtk(resolution=30) + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + def test_no_func_raises(self): + s = Shape() + with pytest.raises(NotImplementedError, match="No implicit field"): + s.generate_vtk() + + def test_no_bounds_raises(self): + s = Shape(func=lambda x, y, z: x**2 + y**2 + z**2 - 1) + with pytest.raises(ValueError, match="Bounds must be provided"): + s.generate_vtk() + + def test_explicit_bounds_override(self): + s = _make_sphere() + mesh = s.generate_vtk(bounds=(-2, 2, -2, 2, -2, 2), resolution=30) + assert mesh.n_cells > 0 + + def test_generateVtk_deprecated(self): + s = _make_sphere() + mesh = s.generateVtk(resolution=20) + assert isinstance(mesh, pv.PolyData) + + +# --------------------------------------------------------------------------- +# Bounds propagation +# --------------------------------------------------------------------------- + + +class TestBoundsPropagation: + """Test bounds merge behavior.""" + + def test_union_expands(self): + s1 = _make_sphere(cx=-2) + s2 = _make_sphere(cx=2) + u = s1 | s2 + assert u.bounds[0] <= s1.bounds[0] + assert u.bounds[1] >= s2.bounds[1] + + def test_intersection_shrinks(self): + s1 = _make_sphere(cx=-0.5, r=1.5) + s2 = _make_sphere(cx=0.5, r=1.5) + i = s1 & s2 + assert i.bounds[0] >= s1.bounds[0] + assert i.bounds[1] <= s2.bounds[1] + + def test_none_bounds(self): + f = lambda x, y, z: x # noqa: E731 + p = Shape(func=f, bounds=None) + s = _make_sphere() + u = p | s + assert u.bounds == s.bounds + + +# --------------------------------------------------------------------------- +# Utilities +# --------------------------------------------------------------------------- + + +class TestUtilities: + """Test utility functions.""" + + def test_shell(self): + s = _make_sphere() + sh = shell(s, thickness=0.2) + # At radius=1.0 (surface), |f|=0, shell should be inside + assert sh.evaluate(np.array([1.0]), np.array([0.0]), np.array([0.0]))[0] < 0 + # At origin, |f|=1.0 >> 0.1, shell should be outside + assert sh.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0]))[0] > 0 + + def test_repeat(self): + s = _make_sphere(r=0.3) + r = repeat(s, spacing=(1.0, 1.0, 1.0)) + # At (1,0,0) should be inside a repeated copy + assert r.evaluate(np.array([1.0]), np.array([0.0]), np.array([0.0]))[0] < 0 + assert r.bounds is None + + def test_blend(self): + s1 = _make_sphere() + s2 = _make_box() + b = blend(s1, s2, factor=0.5) + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + expected = 0.5 * s1.evaluate(x, y, z)[0] + 0.5 * s2.evaluate(x, y, z)[0] + assert b.evaluate(x, y, z)[0] == pytest.approx(expected) + + def test_from_field(self): + shape = from_field( + func=lambda x, y, z: x**2 + y**2 + z**2 - 1, + bounds=(-2, 2, -2, 2, -2, 2), + ) + assert shape.evaluate(np.array([0.0]), np.array([0.0]), np.array([0.0]))[0] < 0 + assert shape.bounds == (-2, 2, -2, 2, -2, 2) + + def test_batch_smooth_union(self): + spheres = [_make_sphere(cx=i * 0.5) for i in range(5)] + combined = batch_smooth_union(spheres, k=0.1) + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + assert combined.evaluate(x, y, z)[0] < 0 + + def test_batch_smooth_union_empty_raises(self): + with pytest.raises(ValueError, match="at least one shape"): + batch_smooth_union([], k=0.1) + + def test_batch_smooth_union_hard(self): + spheres = [_make_sphere(cx=i * 0.5) for i in range(3)] + combined = batch_smooth_union(spheres, k=0.0) + x, y, z = np.array([0.0]), np.array([0.0]), np.array([0.0]) + assert combined.evaluate(x, y, z)[0] < 0 + + +# --------------------------------------------------------------------------- +# Error handling +# --------------------------------------------------------------------------- + + +class TestErrorHandling: + """Test error conditions.""" + + def test_transform_without_func_raises(self): + s = Shape() + with pytest.raises(ValueError, match="No implicit scalar field"): + s.translate((1, 0, 0)) + with pytest.raises(ValueError, match="No implicit scalar field"): + s.rotate((0, 0, 45)) + with pytest.raises(ValueError, match="No implicit scalar field"): + s.scale(2.0) + + def test_boolean_without_func_raises(self): + s = Shape() + other = _make_sphere() + with pytest.raises(ValueError, match="No implicit scalar field"): + s | other # noqa: B015 + with pytest.raises(ValueError, match="No implicit scalar field"): + ~s # noqa: B015 From afc7dbd15b740206ad3c9056773f495d945b71b6 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Thu, 16 Apr 2026 23:06:36 +0200 Subject: [PATCH 02/11] Add autograd to pip dependencies Include autograd in pip dependencies for both the conda build recipe and the environment file (conda.recipe/meta.yaml and environment.yml) so the package is installed during build/test and available in the development environment. --- conda.recipe/meta.yaml | 2 ++ environment.yml | 2 ++ 2 files changed, 4 insertions(+) diff --git a/conda.recipe/meta.yaml b/conda.recipe/meta.yaml index 8064e46a..b2d52b08 100644 --- a/conda.recipe/meta.yaml +++ b/conda.recipe/meta.yaml @@ -24,6 +24,8 @@ requirements: - meshio - cadquery - scipy + - pip: + - autograd test: imports: diff --git a/environment.yml b/environment.yml index a3b13f9e..8ac23927 100755 --- a/environment.yml +++ b/environment.yml @@ -11,5 +11,7 @@ dependencies: - python-gmsh - meshio - cadquery + - pip: + - autograd - mmg # - neper From 3733cbc308d503a1c1d385559a2cbbef8e5a4910 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Thu, 16 Apr 2026 23:07:33 +0200 Subject: [PATCH 03/11] Use autograd.numpy in TPMS surface funcs Replace numpy imports with autograd.numpy in microgen/shape/surface_functions.py so TPMS surface functions are compatible with autograd (enabling automatic differentiation). Add autograd to pyproject.toml dependencies to ensure the new import is available. No behavior changes beyond making the functions differentiable with autograd. --- microgen/shape/surface_functions.py | 4 ++-- pyproject.toml | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/microgen/shape/surface_functions.py b/microgen/shape/surface_functions.py index 0e3bf302..1d8cf2a3 100644 --- a/microgen/shape/surface_functions.py +++ b/microgen/shape/surface_functions.py @@ -1,7 +1,7 @@ """TPMS surface functions.""" -import numpy as np -from numpy import cos, sin +import autograd.numpy as np +from autograd.numpy import cos, sin def gyroid(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: diff --git a/pyproject.toml b/pyproject.toml index e44f03ab..2d0272fb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,7 +18,7 @@ classifiers = [ "Programming Language :: Python :: 3.12", # Python 3.13 not yet supported via pip (VTK 9.3.1 has no 3.13 wheel on PyPI) ] -dependencies = ["numpy", "pyvista", "gmsh", "meshio", "cadquery", "scipy", "nlopt"] +dependencies = ["numpy", "pyvista", "gmsh", "meshio", "cadquery", "scipy", "nlopt", "autograd"] [project.optional-dependencies] dev = [ From 2cac2cb63df36163937e892440186948321963d8 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Thu, 16 Apr 2026 23:07:45 +0200 Subject: [PATCH 04/11] Add F-rep (SDF) integration for TPMS and tests Introduce F-rep based implicit fields for Tpms and its curvilinear variants so TPMS can be used with SDF-normalized operations (shell, complement, booleans). Adds _setup_frep_field to Tpms, CylindricalTpms and SphericalTpms (using from_field + normalize_to_sdf), exposes raw_field, and adds convenience methods as_sheet, as_upper_skeletal and as_lower_skeletal. generate_vtk was updated to build parts from the F-rep Shapes and to use the normalized field/bounds; offset checks and density-based offset computation are preserved. A comprehensive test suite (tests/shapes/test_tpms_frep.py) was added to validate SDF normalization, TPMS F-rep evaluation, sheet/skeletal helpers, boolean composition, generate_vtk compatibility and curvilinear TPMS behavior. --- microgen/shape/tpms.py | 154 +++++++++++++++++- tests/shapes/test_tpms_frep.py | 281 +++++++++++++++++++++++++++++++++ 2 files changed, 433 insertions(+), 2 deletions(-) create mode 100644 tests/shapes/test_tpms_frep.py diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index d238a4fd..90cd0cba 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -123,6 +123,7 @@ def __init__( # noqa: PLR0913 self.resolution = resolution self._compute_tpms_field() + self._setup_frep_field() min_field = np.min(self.grid["surface"]) max_field = np.max(self.grid["surface"]) @@ -354,6 +355,70 @@ def _compute_tpms_field(self: Tpms) -> None: self.grid["surface"] = tpms_field.ravel(order="F") + def _setup_frep_field(self: Tpms) -> None: + """Build the F-rep implicit field (SDF-normalized) for this TPMS.""" + from .implicit_ops import from_field, normalize_to_sdf # noqa: PLC0415 + + k_x, k_y, k_z = 2.0 * np.pi / self.cell_size + ps = self.phase_shift + + def _raw_field( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + return self.surface_function( + k_x * (x + ps[0]), + k_y * (y + ps[1]), + k_z * (z + ps[2]), + ) + + self._raw_field_func = _raw_field + + raw_shape = from_field(_raw_field) + sdf_shape = normalize_to_sdf(raw_shape) + self._func = sdf_shape.func + + half = 0.5 * self.cell_size * self.repeat_cell + self._bounds = ( + -half[0], half[0], + -half[1], half[1], + -half[2], half[2], + ) + + @property + def raw_field(self: Tpms) -> Field: + """The raw (non-SDF-normalized) TPMS field callable.""" + return self._raw_field_func + + def as_sheet(self: Tpms, thickness: float | None = None) -> Shape: + """Return an F-rep Shape representing a uniform-thickness TPMS sheet. + + Uses the SDF-normalized field, so *thickness* is in physical units. + If *thickness* is ``None``, uses ``self.offset``. + """ + from .implicit_ops import shell # noqa: PLC0415 + from .shape import Shape # noqa: PLC0415 + + t = thickness if thickness is not None else self._offset + return shell( + Shape(func=self._func, bounds=self._bounds), + float(t), + ) + + def as_upper_skeletal(self: Tpms) -> Shape: + """Return F-rep Shape for the upper skeletal (f > 0 side).""" + from .implicit_ops import complement # noqa: PLC0415 + from .shape import Shape # noqa: PLC0415 + + return complement(Shape(func=self._func, bounds=self._bounds)) + + def as_lower_skeletal(self: Tpms) -> Shape: + """Return F-rep Shape for the lower skeletal (f < 0 side).""" + from .implicit_ops import from_field # noqa: PLC0415 + + return from_field(func=self._func, bounds=self._bounds) + def _update_grid_offset(self: Tpms) -> None: self.grid["lower_surface"] = self.grid["surface"] + 0.5 * self.offset self.grid["upper_surface"] = self.grid["surface"] - 0.5 * self.offset @@ -629,7 +694,10 @@ def generate_vtk( ) -> pv.PolyData: """Generate VTK PolyData object of the required TPMS part. - :param type_part: part of the TPMS desireds + Uses the SDF-normalized implicit field with F-rep operations: + sheet = ``shell(sdf, thickness)``, skeletal = one side of the SDF. + + :param type_part: part of the TPMS desired :param algo_resolution: if offset must be computed to fit density, \ resolution of the temporary TPMS used to compute the offset @@ -643,12 +711,25 @@ def generate_vtk( "'upper skeletal' or 'surface'" ) raise ValueError(err_msg) + if self.density is not None: self._compute_offset_to_fit_density( part_type=type_part, resolution=algo_resolution, ) - polydata: pv.PolyData = getattr(self, type_part.replace(" ", "_")) + self._check_offset(type_part) + + if type_part == "sheet": + part = self.as_sheet() + elif type_part == "upper skeletal": + part = self.as_upper_skeletal() + else: # lower skeletal + part = self.as_lower_skeletal() + + polydata = part.generate_vtk( + bounds=self._bounds, + resolution=self.resolution * max(self.repeat_cell), + ) polydata = rotate(polydata, center=(0, 0, 0), rotation=self.orientation) return polydata.translate(xyz=self.center) @@ -774,6 +855,39 @@ def _create_grid( ] return grid + def _setup_frep_field(self: CylindricalTpms) -> None: + """Build F-rep field in Cartesian coordinates (inverse cylindrical mapping).""" + from .implicit_ops import from_field, normalize_to_sdf # noqa: PLC0415 + + k_x, k_y, k_z = 2.0 * np.pi / self.cell_size + ps = self.phase_shift + cyl_r = self.cylinder_radius + unit_theta = self.unit_theta + + def _raw_field( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + rho = np.sqrt(x**2 + y**2) - cyl_r + theta = np.arctan2(y, x) / unit_theta + return self.surface_function( + k_x * (rho + ps[0]), + k_y * (theta + ps[1]), + k_z * (z + ps[2]), + ) + + self._raw_field_func = _raw_field + + raw_shape = from_field(_raw_field) + sdf_shape = normalize_to_sdf(raw_shape) + self._func = sdf_shape.func + + # Bounds in Cartesian space + r_max = cyl_r + 0.5 * self.cell_size[0] * self.repeat_cell[0] + half_z = 0.5 * self.cell_size[2] * self.repeat_cell[2] + self._bounds = (-r_max, r_max, -r_max, r_max, -half_z, half_z) + class SphericalTpms(Tpms): """Class used to generate spherical TPMS geometries (sheet or skeletals parts).""" @@ -872,6 +986,42 @@ def _create_grid( ] return grid + def _setup_frep_field(self: SphericalTpms) -> None: + """Build F-rep field in Cartesian coordinates (inverse spherical mapping).""" + from .implicit_ops import from_field, normalize_to_sdf # noqa: PLC0415 + + k_x, k_y, k_z = 2.0 * np.pi / self.cell_size + ps = self.phase_shift + sph_r = self.sphere_radius + unit_theta = self.unit_theta + unit_phi = self.unit_phi + + def _raw_field( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + rho_cart = np.sqrt(x**2 + y**2 + z**2) + rho = rho_cart - sph_r + theta = (np.arccos(np.clip(z / np.maximum(rho_cart, 1e-30), -1, 1)) + - np.pi / 2.0) / unit_theta + phi = np.arctan2(y, x) / unit_phi + return self.surface_function( + k_x * (rho + ps[0]), + k_y * (theta + ps[1]), + k_z * (phi + ps[2]), + ) + + self._raw_field_func = _raw_field + + raw_shape = from_field(_raw_field) + sdf_shape = normalize_to_sdf(raw_shape) + self._func = sdf_shape.func + + # Bounds in Cartesian space + r_max = sph_r + 0.5 * self.cell_size[0] * self.repeat_cell[0] + self._bounds = (-r_max, r_max, -r_max, r_max, -r_max, r_max) + class Infill(Tpms): """Generate a TPMS infill inside a given object.""" diff --git a/tests/shapes/test_tpms_frep.py b/tests/shapes/test_tpms_frep.py new file mode 100644 index 00000000..30a3e9d9 --- /dev/null +++ b/tests/shapes/test_tpms_frep.py @@ -0,0 +1,281 @@ +"""Tests for TPMS F-rep integration with SDF normalization.""" + +from __future__ import annotations + +import numpy as np +import pytest +import pyvista as pv + +from microgen.shape import Tpms, surface_functions +from microgen.shape.implicit_ops import from_field, normalize_to_sdf, shell +from microgen.shape.shape import Shape + + +# --------------------------------------------------------------------------- +# SDF normalization quality +# --------------------------------------------------------------------------- + + +class TestSdfNormalization: + """Test that SDF normalization produces well-behaved fields.""" + + def test_gradient_magnitude_near_one(self): + """After normalization, |grad(sdf)| should be approximately 1.""" + raw = from_field(surface_functions.gyroid) + sdf_shape = normalize_to_sdf(raw) + sdf = sdf_shape.func + + rng = np.random.default_rng(42) + x = rng.uniform(-2, 2, 500) + y = rng.uniform(-2, 2, 500) + z = rng.uniform(-2, 2, 500) + + h = 1e-5 + gx = (sdf(x + h, y, z) - sdf(x - h, y, z)) / (2 * h) + gy = (sdf(x, y + h, z) - sdf(x, y - h, z)) / (2 * h) + gz = (sdf(x, y, z + h) - sdf(x, y, z - h)) / (2 * h) + grad_mag = np.sqrt(gx**2 + gy**2 + gz**2) + + # Near the zero level set, gradient magnitude should be close to 1 + near_surface = np.abs(sdf(x, y, z)) < 0.3 + if near_surface.sum() > 10: + assert np.mean(grad_mag[near_surface]) == pytest.approx(1.0, abs=0.3) + + def test_saddle_point_safety(self): + """SDF at known saddle points should not produce NaN or Inf.""" + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + ) + # Gyroid saddle points are at (0, 0, 0) and equivalents + x = np.array([0.0, 0.25, 0.5]) + y = np.array([0.0, 0.25, 0.5]) + z = np.array([0.0, 0.25, 0.5]) + vals = tpms.evaluate(x, y, z) + assert np.all(np.isfinite(vals)) + + def test_normalize_to_sdf_fallback(self): + """normalize_to_sdf falls back to FD when autograd fails.""" + # A function using regular numpy (not autograd.numpy) + import numpy as regular_np + + def non_autograd_field(x, y, z): + return regular_np.sin(x) + regular_np.cos(y) + + shape = from_field(non_autograd_field, bounds=(-2, 2, -2, 2, -2, 2)) + sdf_shape = normalize_to_sdf(shape) + vals = sdf_shape.func(np.array([0.0]), np.array([0.0]), np.array([0.0])) + assert np.isfinite(vals[0]) + + +# --------------------------------------------------------------------------- +# Tpms F-rep field +# --------------------------------------------------------------------------- + + +class TestTpmsFrepField: + """Test that Tpms has a working F-rep implicit field.""" + + def test_tpms_has_func(self): + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.3) + assert tpms.func is not None + + def test_tpms_has_bounds(self): + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.3) + assert tpms.bounds is not None + + def test_tpms_evaluate(self): + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.3) + x = np.array([0.0, 0.1, 0.2]) + y = np.array([0.0, 0.1, 0.2]) + z = np.array([0.0, 0.1, 0.2]) + vals = tpms.evaluate(x, y, z) + assert vals.shape == (3,) + assert np.all(np.isfinite(vals)) + + def test_tpms_raw_field(self): + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.3) + raw = tpms.raw_field + assert callable(raw) + val = raw(np.array([0.0]), np.array([0.0]), np.array([0.0])) + assert np.isfinite(val[0]) + + +# --------------------------------------------------------------------------- +# F-rep convenience methods +# --------------------------------------------------------------------------- + + +class TestTpmsFrepMethods: + """Test as_sheet, as_upper_skeletal, as_lower_skeletal.""" + + def test_as_sheet(self): + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.3) + sheet = tpms.as_sheet(thickness=0.15) + assert isinstance(sheet, Shape) + assert sheet.func is not None + + def test_as_sheet_mesh(self): + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + resolution=20, + ) + sheet = tpms.as_sheet(thickness=0.15) + mesh = sheet.generate_vtk(bounds=tpms.bounds, resolution=30) + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + def test_as_upper_skeletal(self): + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.3) + skel = tpms.as_upper_skeletal() + assert isinstance(skel, Shape) + assert skel.func is not None + + def test_as_lower_skeletal(self): + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.3) + skel = tpms.as_lower_skeletal() + assert isinstance(skel, Shape) + assert skel.func is not None + + +# --------------------------------------------------------------------------- +# Boolean composition +# --------------------------------------------------------------------------- + + +class TestTpmsBooleans: + """Test boolean composition of TPMS with other shapes.""" + + def test_tpms_and_sphere(self): + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.3) + sphere_func = lambda x, y, z: np.sqrt(x**2 + y**2 + z**2) - 0.4 # noqa: E731 + sphere = from_field(sphere_func, bounds=(-0.5, 0.5, -0.5, 0.5, -0.5, 0.5)) + result = tpms & sphere + assert isinstance(result, Shape) + mesh = result.generate_vtk( + bounds=(-0.5, 0.5, -0.5, 0.5, -0.5, 0.5), + resolution=30, + ) + assert isinstance(mesh, pv.PolyData) + + def test_sheet_and_sphere(self): + tpms = Tpms(surface_function=surface_functions.gyroid, offset=0.3) + sheet = tpms.as_sheet(thickness=0.15) + sphere_func = lambda x, y, z: np.sqrt(x**2 + y**2 + z**2) - 0.4 # noqa: E731 + sphere = from_field(sphere_func, bounds=(-0.5, 0.5, -0.5, 0.5, -0.5, 0.5)) + result = sheet & sphere + mesh = result.generate_vtk( + bounds=(-0.5, 0.5, -0.5, 0.5, -0.5, 0.5), + resolution=30, + ) + assert isinstance(mesh, pv.PolyData) + + +# --------------------------------------------------------------------------- +# generate_vtk backward compatibility +# --------------------------------------------------------------------------- + + +class TestTpmsGenerateVtk: + """Test that generate_vtk(type_part=...) still works.""" + + def test_sheet(self): + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + resolution=20, + ) + mesh = tpms.generate_vtk(type_part="sheet") + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + def test_surface(self): + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + resolution=20, + ) + mesh = tpms.generate_vtk(type_part="surface") + assert isinstance(mesh, pv.PolyData) + assert mesh.n_cells > 0 + + def test_invalid_type_part(self): + tpms = Tpms( + surface_function=surface_functions.gyroid, + offset=0.3, + ) + with pytest.raises(ValueError, match="type_part"): + tpms.generate_vtk(type_part="invalid") + + @pytest.mark.parametrize( + "surface_fn", + [surface_functions.gyroid, surface_functions.schwarz_p], + ) + def test_multiple_surface_functions(self, surface_fn): + tpms = Tpms(surface_function=surface_fn, offset=0.3, resolution=20) + mesh = tpms.generate_vtk(type_part="sheet") + assert mesh.n_cells > 0 + + +# --------------------------------------------------------------------------- +# CylindricalTpms / SphericalTpms +# --------------------------------------------------------------------------- + + +class TestCurvilinearTpms: + """Test F-rep field on curvilinear TPMS variants.""" + + def test_cylindrical_has_func(self): + from microgen.shape.tpms import CylindricalTpms + + tpms = CylindricalTpms( + radius=1.0, + surface_function=surface_functions.gyroid, + offset=0.3, + resolution=10, + ) + assert tpms.func is not None + assert tpms.bounds is not None + + def test_cylindrical_evaluate(self): + from microgen.shape.tpms import CylindricalTpms + + tpms = CylindricalTpms( + radius=1.0, + surface_function=surface_functions.gyroid, + offset=0.3, + resolution=10, + ) + x = np.array([1.0, 0.5]) + y = np.array([0.0, 0.5]) + z = np.array([0.0, 0.0]) + vals = tpms.evaluate(x, y, z) + assert np.all(np.isfinite(vals)) + + def test_spherical_has_func(self): + from microgen.shape.tpms import SphericalTpms + + tpms = SphericalTpms( + radius=1.0, + surface_function=surface_functions.gyroid, + offset=0.3, + resolution=10, + ) + assert tpms.func is not None + assert tpms.bounds is not None + + def test_spherical_evaluate(self): + from microgen.shape.tpms import SphericalTpms + + tpms = SphericalTpms( + radius=1.0, + surface_function=surface_functions.gyroid, + offset=0.3, + resolution=10, + ) + x = np.array([1.0, 0.5]) + y = np.array([0.0, 0.5]) + z = np.array([0.0, 0.5]) + vals = tpms.evaluate(x, y, z) + assert np.all(np.isfinite(vals)) From 945c0f66a4742d4b52d8dacfae44ab9c15836ec8 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Fri, 17 Apr 2026 14:12:36 +0200 Subject: [PATCH 05/11] Refactor TPMS field finalization and VTK resolution Introduce _finalize_frep to centralize raw field -> SDF normalization and setting of _raw_field_func, _func and _bounds. Replace duplicated normalization code in Tpms, CylindricalTpms and SphericalTpms with calls to the new helper. Adjust generate_vtk resolution logic to compute an isotropic resolution from repeat_cell (using geometric mean) and enforce a sensible minimum, and wrap shell return in a plain Shape to avoid recursive generate_vtk calls. Also include minor formatting and expression cleanups. --- microgen/shape/tpms.py | 78 ++++++++++++++++++++---------------------- 1 file changed, 37 insertions(+), 41 deletions(-) diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 90cd0cba..6081a882 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -355,10 +355,21 @@ def _compute_tpms_field(self: Tpms) -> None: self.grid["surface"] = tpms_field.ravel(order="F") - def _setup_frep_field(self: Tpms) -> None: - """Build the F-rep implicit field (SDF-normalized) for this TPMS.""" + def _finalize_frep( + self: Tpms, + raw_field: Field, + bounds: tuple[float, float, float, float, float, float], + ) -> None: + """Normalize a raw field to SDF and set ``_func`` / ``_bounds``.""" from .implicit_ops import from_field, normalize_to_sdf # noqa: PLC0415 + self._raw_field_func = raw_field + sdf_shape = normalize_to_sdf(from_field(raw_field)) + self._func = sdf_shape.func + self._bounds = bounds + + def _setup_frep_field(self: Tpms) -> None: + """Build the F-rep implicit field (SDF-normalized) for this TPMS.""" k_x, k_y, k_z = 2.0 * np.pi / self.cell_size ps = self.phase_shift @@ -373,17 +384,10 @@ def _raw_field( k_z * (z + ps[2]), ) - self._raw_field_func = _raw_field - - raw_shape = from_field(_raw_field) - sdf_shape = normalize_to_sdf(raw_shape) - self._func = sdf_shape.func - half = 0.5 * self.cell_size * self.repeat_cell - self._bounds = ( - -half[0], half[0], - -half[1], half[1], - -half[2], half[2], + self._finalize_frep( + _raw_field, + (-half[0], half[0], -half[1], half[1], -half[2], half[2]), ) @property @@ -401,10 +405,9 @@ def as_sheet(self: Tpms, thickness: float | None = None) -> Shape: from .shape import Shape # noqa: PLC0415 t = thickness if thickness is not None else self._offset - return shell( - Shape(func=self._func, bounds=self._bounds), - float(t), - ) + # Wrap in plain Shape so generate_vtk uses Shape's marching cubes, + # not Tpms.generate_vtk (which would recurse back here). + return shell(Shape(func=self._func, bounds=self._bounds), float(t)) def as_upper_skeletal(self: Tpms) -> Shape: """Return F-rep Shape for the upper skeletal (f > 0 side).""" @@ -726,10 +729,13 @@ def generate_vtk( else: # lower skeletal part = self.as_lower_skeletal() - polydata = part.generate_vtk( - bounds=self._bounds, - resolution=self.resolution * max(self.repeat_cell), + # Match the grid resolution: resolution * repeat_cell per axis. + # Shape.generate_vtk takes a single isotropic resolution, so use the + # geometric mean to keep total point count proportional. + iso_res = int( + self.resolution * np.cbrt(np.prod(self.repeat_cell)) ) + polydata = part.generate_vtk(bounds=self._bounds, resolution=max(iso_res, 10)) polydata = rotate(polydata, center=(0, 0, 0), rotation=self.orientation) return polydata.translate(xyz=self.center) @@ -857,8 +863,6 @@ def _create_grid( def _setup_frep_field(self: CylindricalTpms) -> None: """Build F-rep field in Cartesian coordinates (inverse cylindrical mapping).""" - from .implicit_ops import from_field, normalize_to_sdf # noqa: PLC0415 - k_x, k_y, k_z = 2.0 * np.pi / self.cell_size ps = self.phase_shift cyl_r = self.cylinder_radius @@ -877,16 +881,12 @@ def _raw_field( k_z * (z + ps[2]), ) - self._raw_field_func = _raw_field - - raw_shape = from_field(_raw_field) - sdf_shape = normalize_to_sdf(raw_shape) - self._func = sdf_shape.func - - # Bounds in Cartesian space r_max = cyl_r + 0.5 * self.cell_size[0] * self.repeat_cell[0] half_z = 0.5 * self.cell_size[2] * self.repeat_cell[2] - self._bounds = (-r_max, r_max, -r_max, r_max, -half_z, half_z) + self._finalize_frep( + _raw_field, + (-r_max, r_max, -r_max, r_max, -half_z, half_z), + ) class SphericalTpms(Tpms): @@ -988,8 +988,6 @@ def _create_grid( def _setup_frep_field(self: SphericalTpms) -> None: """Build F-rep field in Cartesian coordinates (inverse spherical mapping).""" - from .implicit_ops import from_field, normalize_to_sdf # noqa: PLC0415 - k_x, k_y, k_z = 2.0 * np.pi / self.cell_size ps = self.phase_shift sph_r = self.sphere_radius @@ -1003,8 +1001,10 @@ def _raw_field( ) -> npt.NDArray[np.float64]: rho_cart = np.sqrt(x**2 + y**2 + z**2) rho = rho_cart - sph_r - theta = (np.arccos(np.clip(z / np.maximum(rho_cart, 1e-30), -1, 1)) - - np.pi / 2.0) / unit_theta + theta = ( + np.arccos(np.clip(z / np.maximum(rho_cart, 1e-30), -1, 1)) + - np.pi / 2.0 + ) / unit_theta phi = np.arctan2(y, x) / unit_phi return self.surface_function( k_x * (rho + ps[0]), @@ -1012,15 +1012,11 @@ def _raw_field( k_z * (phi + ps[2]), ) - self._raw_field_func = _raw_field - - raw_shape = from_field(_raw_field) - sdf_shape = normalize_to_sdf(raw_shape) - self._func = sdf_shape.func - - # Bounds in Cartesian space r_max = sph_r + 0.5 * self.cell_size[0] * self.repeat_cell[0] - self._bounds = (-r_max, r_max, -r_max, r_max, -r_max, r_max) + self._finalize_frep( + _raw_field, + (-r_max, r_max, -r_max, r_max, -r_max, r_max), + ) class Infill(Tpms): From bc2b485826ccb4b4fbd8d6faf9cd4f1d5c668611 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Sat, 25 Apr 2026 13:14:50 +0200 Subject: [PATCH 06/11] Support variable offsets and F-rep TPMS refactor MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add variable-thickness and SDF utilities and refactor TPMS F-rep pipeline. Key changes: - implicit_ops.shell: allow thickness to be a callable Field (spatially-varying) and add docs; add box() SDF primitive; improve gradient-normalization fallback to preserve sign when grad is degenerate. - tpms.Tpms: store offset callables, accept callable offsets and preserve them for re-evaluation; compute density→offset via the same marching-cubes pipeline (generate_vtk) and handle bracket failures; expose as_sheet that accepts variable offsets; add helpers (_half_offset_field, _cell_box, clipped variants) and pick F-rep shapes for parts. - Mesh -> CadQuery conversion: replace previous face assembly with sewing (BRepBuilderAPI_Sewing) to produce a proper shell; best-effort upgrade to solids via _try_make_solid; generate() now uses the pure F-rep path (marching cubes on SDF Shape) and applies smoothing/solid conversion; generate_vtk mirrors the same F-rep pipeline and short-circuits full-density case by returning a box mesh. - Misc: add _isotropic_resolution to map per-axis resolution to a single isotropic value; improve argument validation and error messages; reset cached _offset_func on assignment. - tests: update TPMS surface function discovery to filter out camelCase aliases and non-TPMS callables by checking function signatures (3 parameters). These changes make TPMS generation consistent between the VTK and CAD paths, enable spatially-varying thickness, and improve robustness around degenerate fields and solid construction. --- microgen/shape/implicit_ops.py | 61 +++- microgen/shape/tpms.py | 502 ++++++++++++++++++++------------- tests/shapes/test_tpms.py | 17 +- 3 files changed, 368 insertions(+), 212 deletions(-) diff --git a/microgen/shape/implicit_ops.py b/microgen/shape/implicit_ops.py index a269c08a..f42b6518 100644 --- a/microgen/shape/implicit_ops.py +++ b/microgen/shape/implicit_ops.py @@ -224,10 +224,23 @@ def _batched( # --------------------------------------------------------------------------- -def shell(shape: Shape, thickness: float) -> Shape: - """Hollow shell: ``|f(p)| - thickness / 2``.""" +def shell(shape: Shape, thickness: float | Field) -> Shape: + """Hollow shell: ``|f(p)| - thickness(p) / 2``. + + ``thickness`` may be a constant (scalar) or a callable + ``thickness(x, y, z) -> array`` for spatially-varying shells. Negative or + zero thickness at a point yields no inclusion in the shell at that point. + """ f = shape.require_func() - half_t = thickness / 2.0 + if callable(thickness): + t_func = thickness + + def _shell_field(x, y, z, _f=f, _t=t_func): # noqa: ANN001 + return np.abs(_f(x, y, z)) - _t(x, y, z) / 2.0 + + return _make_shape(func=_shell_field, bounds=shape.bounds) + + half_t = float(thickness) / 2.0 return _make_shape( func=lambda x, y, z, _f=f, _ht=half_t: np.abs(_f(x, y, z)) - _ht, bounds=shape.bounds, @@ -312,6 +325,37 @@ def from_field( return _make_shape(func=func, bounds=bounds) +def box( + dims: tuple[float, float, float], + center: tuple[float, float, float] = (0.0, 0.0, 0.0), +) -> Shape: + """Axis-aligned box as an F-rep Shape. + + SDF formula ``max(|x-cx|-hx, |y-cy|-hy, |z-cz|-hz)``: signed distance to + the box surface (negative inside, positive outside, zero on the surface). + Useful as a clipping primitive — e.g. ``intersection(skeletal, box(...))`` + bounds an unbounded TPMS skeletal field to a single cell. + + :param dims: full side lengths ``(dx, dy, dz)`` + :param center: box center (default origin) + :return: :class:`~microgen.shape.shape.Shape` carrying the box SDF + """ + hx, hy, hz = (0.5 * float(d) for d in dims) + cx, cy, cz = (float(c) for c in center) + + def _box_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + return np.maximum.reduce( + [np.abs(x - cx) - hx, np.abs(y - cy) - hy, np.abs(z - cz) - hz], + ) + + bounds: BoundsType = (cx - hx, cx + hx, cy - hy, cy + hy, cz - hz, cz + hz) + return _make_shape(func=_box_sdf, bounds=bounds) + + def _fd_sdf( f: Field, epsilon: float, @@ -329,7 +373,11 @@ def sdf( gy = (f(x, y + h, z) - f(x, y - h, z)) / (2 * h) gz = (f(x, y, z + h) - f(x, y, z - h)) / (2 * h) grad_mag = np.sqrt(gx**2 + gy**2 + gz**2) - return val / np.maximum(grad_mag, epsilon) + # Where the gradient is degenerate (e.g. flat-z fields like the + # honeycomb_* surfaces, or saddle points), normalization would + # blow up to ±1/epsilon — preserve the raw field's *sign* by + # falling back to the unnormalized value there. + return np.where(grad_mag > epsilon, val / np.maximum(grad_mag, epsilon), val) return sdf @@ -374,7 +422,10 @@ def sdf( grad_mag = np.sqrt( dfdx(x, y, z) ** 2 + dfdy(x, y, z) ** 2 + dfdz(x, y, z) ** 2, ) - return val / np.maximum(grad_mag, epsilon) + # Same fallback as in `_fd_sdf`: where the gradient vanishes + # (degenerate flat-z fields, saddle points), keep the raw value + # so its sign is preserved without exploding into ±1/epsilon. + return np.where(grad_mag > epsilon, val / np.maximum(grad_mag, epsilon), val) except Exception: # noqa: BLE001 sdf = _fd_sdf(f, epsilon) diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 6081a882..6e6b9fa4 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -110,6 +110,11 @@ def __init__( # noqa: PLR0913 self.surface_function = surface_function self._offset = offset if offset is not None else 0.0 + # Stores the offset callable when one is provided, so the F-rep path + # can re-evaluate variable thickness on its own marching-cubes grid. + self._offset_func: Field | None = ( + offset if (offset is not None and callable(offset) and not isinstance(offset, OffsetGrading)) else None + ) self.phase_shift = phase_shift self.grid: pv.StructuredGrid @@ -177,11 +182,20 @@ def _compute_offset_to_fit_density( part_type: Literal["sheet", "lower skeletal", "upper skeletal"], resolution: int | None = None, ) -> float: - """Compute the offset to fit the required density.""" + """Compute the offset that yields the requested density. + + Searches with the same F-rep ``generate_vtk`` pipeline the user + invokes, so the offset returned actually reproduces the requested + density at the user-facing resolution. When the target density is + too high to reach at this resolution (marching cubes saturates + slightly below 1.0 due to surface-grid discretization), falls back + to the offset limit instead of failing the bracket. + """ if self.density is None: err_msg = f"density must be between 0 and 1. Given: {self.density}" raise ValueError(err_msg) + part = "skeletal" if "skeletal" in part_type else part_type if self.density == 1.0: self.offset = ( self.offset_lim["sheet"][1] @@ -193,19 +207,31 @@ def _compute_offset_to_fit_density( temp_tpms = Tpms( surface_function=self.surface_function, offset=0.0, + cell_size=self.cell_size, + repeat_cell=self.repeat_cell, resolution=resolution if resolution is not None else self.resolution, ) + cell_volume = abs(temp_tpms.grid.volume) def density(offset: float) -> float: temp_tpms.offset = offset - grid_part = getattr(temp_tpms, f"grid_{part_type.replace(' ', '_')}") - return abs(grid_part.volume) + mesh = temp_tpms.generate_vtk(type_part=part_type) + return abs(mesh.volume) / cell_volume + + bracket = self.offset_lim[part] + try: + computed_offset = root_scalar( + lambda offset: density(offset) - self.density, + bracket=bracket, + ).root + except ValueError: + # Bracket sign mismatch — usually because target density exceeds + # the achievable max at this resolution. Pick the bracket + # endpoint that lies on the side of self.density. + d_lo = density(bracket[0]) - self.density + d_hi = density(bracket[1]) - self.density + computed_offset = bracket[1] if abs(d_hi) <= abs(d_lo) else bracket[0] - part = "skeletal" if "skeletal" in part_type else part_type - computed_offset = root_scalar( - lambda offset: density(offset) - self.density, - bracket=self.offset_lim[part], - ).root self.offset = computed_offset return computed_offset @@ -396,32 +422,127 @@ def raw_field(self: Tpms) -> Field: return self._raw_field_func def as_sheet(self: Tpms, thickness: float | None = None) -> Shape: - """Return an F-rep Shape representing a uniform-thickness TPMS sheet. + """Return an F-rep Shape representing a TPMS sheet of given thickness. Uses the SDF-normalized field, so *thickness* is in physical units. - If *thickness* is ``None``, uses ``self.offset``. + If *thickness* is ``None``, uses ``self.offset`` (which may be a + scalar, an array sampled on ``self.grid``, or a callable in which + case the callable form is used directly). """ from .implicit_ops import shell # noqa: PLC0415 from .shape import Shape # noqa: PLC0415 - t = thickness if thickness is not None else self._offset - # Wrap in plain Shape so generate_vtk uses Shape's marching cubes, - # not Tpms.generate_vtk (which would recurse back here). - return shell(Shape(func=self._func, bounds=self._bounds), float(t)) + if thickness is not None: + t: float | npt.NDArray | Field = float(thickness) + elif self._offset_func is not None: + t = self._offset_func + else: + t = self._offset + return shell(Shape(func=self._func, bounds=self._bounds), t) + + def _half_offset_field(self: Tpms) -> Field | float: + """Return half the offset as a callable (variable) or scalar (constant). + + Variable offset stored as a callable can be re-evaluated on the + marching-cubes grid; an array offset (sampled on ``self.grid``) + cannot be remapped to arbitrary points, so the only safe fallback + for arrays is to use a scalar 0 (sheet/skeletal degenerate to the + zero-isosurface). + """ + if self._offset_func is not None: + f = self._offset_func + + def _half(x, y, z, _f=f): # noqa: ANN001 + return 0.5 * _f(x, y, z) + + return _half + if isinstance(self._offset, (int, float)): + return 0.5 * float(self._offset) + # array — no safe re-evaluation; degenerate to zero (skeletal at f=0). + return 0.0 def as_upper_skeletal(self: Tpms) -> Shape: - """Return F-rep Shape for the upper skeletal (f > 0 side).""" - from .implicit_ops import complement # noqa: PLC0415 - from .shape import Shape # noqa: PLC0415 + """F-rep Shape for the *upper* skeletal: ``{p : f(p) > offset/2}``. + + Volume scales with the chosen offset (smaller offset ⇒ larger + skeletal), matching the historical CadQuery behaviour and the VTK + grid-clip path. + """ + from .implicit_ops import from_field # noqa: PLC0415 + + f = self._func + h = self._half_offset_field() + if callable(h): + + def _upper(x, y, z, _f=f, _h=h): # noqa: ANN001 + return -_f(x, y, z) + _h(x, y, z) + + else: + + def _upper(x, y, z, _f=f, _h=h): # noqa: ANN001 + return -_f(x, y, z) + _h - return complement(Shape(func=self._func, bounds=self._bounds)) + return from_field(func=_upper, bounds=self._bounds) def as_lower_skeletal(self: Tpms) -> Shape: - """Return F-rep Shape for the lower skeletal (f < 0 side).""" + """F-rep Shape for the *lower* skeletal: ``{p : f(p) < -offset/2}``.""" + from .implicit_ops import from_field # noqa: PLC0415 + + f = self._func + h = self._half_offset_field() + if callable(h): + + def _lower(x, y, z, _f=f, _h=h): # noqa: ANN001 + return _f(x, y, z) + _h(x, y, z) + + else: + + def _lower(x, y, z, _f=f, _h=h): # noqa: ANN001 + return _f(x, y, z) + _h + + return from_field(func=_lower, bounds=self._bounds) + + def as_surface(self: Tpms) -> Shape: + """Return F-rep Shape for the (open) zero-isosurface, no thickness. + + Same field as :meth:`as_lower_skeletal`; meant for ``type_part="surface"``. + Marching cubes will produce an open shell — there is no enclosed volume. + """ from .implicit_ops import from_field # noqa: PLC0415 return from_field(func=self._func, bounds=self._bounds) + def _cell_box(self: Tpms) -> Shape: + """SDF Shape of this TPMS' cell (cell_size × repeat_cell, centered origin).""" + from .implicit_ops import box # noqa: PLC0415 + + dims = tuple(float(d) for d in (self.cell_size * self.repeat_cell)) + return box(dims=dims, center=(0.0, 0.0, 0.0)) + + def _clipped_sheet(self: Tpms) -> Shape: + """Sheet F-rep clipped to the cell box. + + When the offset approaches its upper limit the sheet field is + uniformly negative inside the cell and marching cubes finds no + boundary — clipping by the cell box makes the box face the closed + boundary, yielding a full-cell mesh as expected. + """ + from .implicit_ops import intersection # noqa: PLC0415 + + return intersection(self.as_sheet(), self._cell_box()) + + def _clipped_upper_skeletal(self: Tpms) -> Shape: + """Upper skeletal F-rep clipped to the cell box (closed under marching cubes).""" + from .implicit_ops import intersection # noqa: PLC0415 + + return intersection(self.as_upper_skeletal(), self._cell_box()) + + def _clipped_lower_skeletal(self: Tpms) -> Shape: + """Lower skeletal F-rep clipped to the cell box (closed under marching cubes).""" + from .implicit_ops import intersection # noqa: PLC0415 + + return intersection(self.as_lower_skeletal(), self._cell_box()) + def _update_grid_offset(self: Tpms) -> None: self.grid["lower_surface"] = self.grid["surface"] + 0.5 * self.offset self.grid["upper_surface"] = self.grid["surface"] - 0.5 * self.offset @@ -436,11 +557,16 @@ def offset( self: Tpms, offset: float | npt.NDArray[np.float64] | OffsetGrading | Field, ) -> None: + # Reset cached callable form on every assignment. + self._offset_func = None if isinstance(offset, (int, float, np.ndarray)): self._offset = offset elif isinstance(offset, OffsetGrading): self._offset = offset.compute_offset(self.grid) elif callable(offset): + # Keep the callable so the F-rep path can re-evaluate it on the + # marching-cubes grid (which differs from ``self.grid``). + self._offset_func = offset self._offset = offset(self.grid.x, self.grid.y, self.grid.z).ravel("F") else: err_msg = "offset must be a float, a numpy array or a callable" @@ -449,9 +575,20 @@ def offset( self._update_grid_offset() self.offset_updated = True - def _create_shell(self: Tpms, mesh: pv.PolyData) -> cq.Shell: + def _mesh_to_shell(self: Tpms, mesh: pv.PolyData) -> cq.Shape: + """Convert a triangulated PyVista mesh to a CadQuery ``Shape``. + + Builds one ``cq.Face`` per triangle, then *sews* them with OCCT's + ``BRepBuilderAPI_Sewing`` so adjacent triangles share edges — without + sewing, ``Shell.makeShell`` returns a disjoint compound and the result + cannot be turned into a Solid. Returns the sewn shape (a + ``TopoDS_Shell`` if sewing succeeded into one shell, otherwise the + compound that sewing produced). + """ + from OCP.BRepBuilderAPI import BRepBuilderAPI_Sewing # noqa: PLC0415 + if not mesh.is_all_triangles: - mesh.triangulate(inplace=True) # useless ? + mesh.triangulate(inplace=True) triangles = mesh.faces.reshape(-1, 4)[:, 1:] triangles = np.c_[triangles, triangles[:, 0]] @@ -464,147 +601,18 @@ def _create_shell(self: Tpms, mesh: pv.PolyData) -> cq.Shell: ) for start, end in zip(tri[:], tri[1:]) ] - wire = cq.Wire.assembleEdges(lines) faces.append(cq.Face.makeFromWires(wire)) - try: - shell = cq.Shell.makeShell(faces) - except ValueError as err: - err_msg = "Failed to create the shell, \ - try to increase the resolution or the smoothing." - raise ShellCreationError(err_msg) from err - return shell - - def _create_surface( - self: Tpms, - isovalue: float | npt.NDArray[np.float64] = 0.0, - smoothing: int = 0, - ) -> cq.Shell: - """Create a TPMS surface for the given isovalue.""" - if isinstance(isovalue, (int, float)): - scalars = self.grid["surface"] - isovalue - elif isinstance(isovalue, np.ndarray): - scalars = self.grid["surface"] - isovalue.ravel(order="F") - - mesh = self.grid.contour(isosurfaces=[0.0], scalars=scalars) - mesh.smooth(n_iter=smoothing, feature_smoothing=True, inplace=True) - mesh.clean(inplace=True) - - return self._create_shell(mesh=mesh) - - def _create_surfaces( - self: Tpms, - isovalues: list[float], - smoothing: int = 0, - ) -> list[cq.Shell]: - """Create TPMS surfaces for the corresponding isovalue. - - :param isovalues: list of isovalues corresponding to the required surfaces - :param smoothing: smoothing loop iterations - :param verbose: display progressbar of the conversion to CadQuery object + if not faces: + err_msg = "Mesh has no triangles to convert into a shell" + raise ShellCreationError(err_msg) - :return: list of CadQuery Shell objects of the required TPMS surfaces - """ - shells = [] - for i, isovalue in enumerate(isovalues): - logging.info("\nGenerating surface (%d/%d)", i + 1, len(isovalues)) - shell = self._create_surface( - isovalue=isovalue, - smoothing=smoothing, - ) - shells.append(shell) - - return shells - - def _generate_sheet_surfaces( - self: Tpms, - smoothing: int, - ) -> tuple[cq.Shape, cq.Shape]: - """Generate the surfaces to create the sheet part of the TPMS.""" - isovalues = [ - -0.5 * self.offset, - -0.25 * self.offset, - 0.25 * self.offset, - 0.5 * self.offset, - ] - - shells = self._create_surfaces(isovalues, smoothing) - - lower_surface = shells[0] - lower_test_surface = shells[1] - upper_test_surface = shells[2] - upper_surface = shells[3] - - surface = lower_surface.fuse(upper_surface) - test_surface = lower_test_surface.fuse(upper_test_surface) - return surface, test_surface - - def _generate_lower_skeletal_surfaces( - self: Tpms, - smoothing: int, - ) -> tuple[cq.Shape, cq.Shape]: - """Generate the surfaces to create the lower skeletal part of the TPMS.""" - min_offset = 2.0 * np.min(self.grid["surface"]) - isovalues = [ - -0.5 * self.offset, - -0.25 * (self.offset - min_offset), - ] - - shells = self._create_surfaces(isovalues, smoothing) - - surface = shells[0] - test_surface = shells[1] - return surface, test_surface - - def _generate_upper_skeletal_surfaces( - self: Tpms, - smoothing: int, - ) -> tuple[cq.Shape, cq.Shape]: - """Generate the surfaces to create the upper skeletal part of the TPMS.""" - min_offset = 2.0 * np.min(self.grid["surface"]) - isovalues = [ - 0.5 * self.offset, - 0.25 * (self.offset - min_offset), - ] - - shells = self._create_surfaces(isovalues, smoothing) - - surface = shells[0] - test_surface = shells[1] - return surface, test_surface - - def _extract_part_from_box( - self: Tpms, - type_part: TpmsPartType, - smoothing: int, - ) -> cq.Shape: - """Extract the required part from the box.""" - box = cq.Workplane("front").box(*(self.cell_size * self.repeat_cell)) - - surface, test_surface = getattr( - self, - f"_generate_{type_part.replace(' ', '_')}_surfaces", - )(smoothing) - - splitted_box = box.split(surface) - tpms_solids = splitted_box.solids().all() - - # split each solid with the test surface to identify - # to what part type the solid belongs to - list_solids = [ - (solid.split(test_surface).solids().size(), solid.val()) - for solid in tpms_solids - ] - - # if the number of shapes is greater than 1, it means that the solid is split - # so it belongs to the required part - part_solids = [solid for (number, solid) in list_solids if number > 1] - part_shapes = [cq.Shape(solid.wrapped) for solid in part_solids] - return fuseShapes( - cqShapeList=part_shapes, - retain_edges=False, # True or False ? - ) + sewing = BRepBuilderAPI_Sewing() + for f in faces: + sewing.Add(f.wrapped) + sewing.Perform() + return cq.Shape(sewing.SewedShape()) def _check_offset( self: Tpms, @@ -641,6 +649,35 @@ def _check_offset( generate '{type_part}' part and lower than {self.offset_lim[part][1]}" raise ValueError(err_msg) + _VALID_PARTS = ("sheet", "lower skeletal", "upper skeletal", "surface") + + def _frep_part(self: Tpms, type_part: TpmsPartType) -> Shape: + """Pick the F-rep :class:`Shape` for *type_part*. + + Skeletals are intersected with the cell box so marching cubes + produces a *closed* shell (the unclipped skeletal field is + unbounded). ``"surface"`` and ``"sheet"`` are already bounded by + construction (zero-isosurface and `shell()`-clipped, respectively). + """ + if type_part == "sheet": + return self._clipped_sheet() + if type_part == "upper skeletal": + return self._clipped_upper_skeletal() + if type_part == "lower skeletal": + return self._clipped_lower_skeletal() + if type_part == "surface": + return self.as_surface() + err_msg = f"type_part {type_part!r} must be one of {self._VALID_PARTS}" + raise ValueError(err_msg) + + def _isotropic_resolution(self: Tpms) -> int: + """Map ``self.resolution`` (per-axis) to an isotropic Shape resolution. + + ``Shape.generate_vtk`` takes a single resolution; we use the geometric + mean of the per-axis cell counts so total grid points stay proportional. + """ + return max(int(self.resolution * np.cbrt(np.prod(self.repeat_cell))), 10) + def generate( self: Tpms, type_part: TpmsPartType = "sheet", @@ -648,18 +685,21 @@ def generate( algo_resolution: int | None = None, **_: KwargsGenerateType, ) -> cq.Shape: - """Generate CadQuery Shape object of the required TPMS part. - - :param type_part: part of the TPMS desired \ - ('sheet', 'lower skeletal', 'upper skeletal' or 'surface') - :param smoothing: smoothing loop iterations - :param verbose: display progressbar of the conversion to CadQuery object - :param algo_resolution: if offset must be computed to fit density, \ - resolution of the temporary TPMS used to compute the offset - - :return: CadQuery Shape object of the required TPMS part + """Generate the OCCT/CadQuery shape of the requested TPMS part. + + Pure F-rep pipeline: pick the SDF Shape via :meth:`_frep_part`, run + marching cubes through :meth:`Shape.generate_vtk`, optionally smooth, + then build a CadQuery ``Shell``. The same SDF + same marching-cubes + grid is used by :meth:`generate_vtk`, so volumes converge to identical + values up to discretization. + + :param type_part: ``"sheet"``, ``"lower skeletal"``, ``"upper skeletal"`` + or ``"surface"`` (open zero-isosurface, no thickness) + :param smoothing: number of Laplacian smoothing iterations on the mesh + :param algo_resolution: temporary-TPMS resolution for density→offset + search (only used when ``self.density`` is set) """ - if type_part not in ["sheet", "lower skeletal", "upper skeletal", "surface"]: + if type_part not in self._VALID_PARTS: err_msg = ( f"type_part ({type_part}) must be 'sheet', 'lower skeletal', " "'upper skeletal' or 'surface'" @@ -671,40 +711,96 @@ def generate( logging.warning("offset is ignored for 'surface' part") if self.density is not None: logging.warning("density is ignored for 'surface' part") - return self._create_surface( - isovalue=0, - smoothing=smoothing, - ) + else: + if self.density is not None: + self._compute_offset_to_fit_density( + part_type=type_part, + resolution=algo_resolution, + ) + self._check_offset(type_part) - if self.density is not None: - self._compute_offset_to_fit_density( - part_type=type_part, - resolution=algo_resolution, + frep = self._frep_part(type_part) + mesh = frep.generate_vtk( + bounds=self._bounds, + resolution=self._isotropic_resolution(), + ) + if smoothing > 0: + mesh.smooth(n_iter=smoothing, feature_smoothing=True, inplace=True) + mesh.clean(inplace=True) + + if mesh.n_cells == 0: + err_msg = ( + f"Marching cubes produced an empty mesh for '{type_part}'; " + "check offset / density / resolution." ) + raise ShellCreationError(err_msg) - self._check_offset(type_part) + shape = self._mesh_to_shell(mesh) - shape = self._extract_part_from_box(type_part, smoothing) + if type_part != "surface": + # Closed shell ⇒ try to upgrade into a Solid so volume queries and + # downstream booleans behave correctly. If sewing produced a + # compound (non-manifold patches) or OCCT refuses, keep the shell. + shape = self._try_make_solid(shape) shape = rotate(obj=shape, center=(0, 0, 0), rotation=self.orientation) return shape.translate(self.center) + @staticmethod + def _try_make_solid(shape: cq.Shape) -> cq.Shape: + """Best-effort upgrade of a sewn shell into a closed Solid. + + Returns the original shape unchanged if the sewn result is a Compound + (multiple disjoint shells, can't be a single solid) or if OCCT refuses + the conversion. + """ + from OCP.TopAbs import TopAbs_SHELL # noqa: PLC0415 + from OCP.TopExp import TopExp_Explorer # noqa: PLC0415 + from OCP.TopoDS import TopoDS # noqa: PLC0415 + + wrapped = shape.wrapped + # Already a Shell? Try to make a Solid directly. + if wrapped.ShapeType() == TopAbs_SHELL: + try: + shell = TopoDS.Shell_s(wrapped) + return cq.Shape(cq.Solid.makeSolid(cq.Shell(shell)).wrapped) + except (ValueError, RuntimeError): + return shape + + # Compound: extract Shells, build a Solid per closed shell, fuse. + exp = TopExp_Explorer(wrapped, TopAbs_SHELL) + solids: list[cq.Shape] = [] + while exp.More(): + try: + shell = TopoDS.Shell_s(exp.Current()) + solid = cq.Solid.makeSolid(cq.Shell(shell)) + solids.append(cq.Shape(solid.wrapped)) + except (ValueError, RuntimeError): + pass + exp.Next() + + if not solids: + return shape + if len(solids) == 1: + return solids[0] + return fuseShapes(cqShapeList=solids, retain_edges=False) + def generate_vtk( self: Tpms, type_part: TpmsPartType = "sheet", algo_resolution: int | None = None, **_: KwargsGenerateType, ) -> pv.PolyData: - """Generate VTK PolyData object of the required TPMS part. + """Generate the PyVista mesh of the requested TPMS part. - Uses the SDF-normalized implicit field with F-rep operations: - sheet = ``shell(sdf, thickness)``, skeletal = one side of the SDF. + Same F-rep pipeline as :meth:`generate` (skeletals are clipped to the + cell box), so the two outputs share the exact same triangulation and + therefore the same volume. - :param type_part: part of the TPMS desired - :param algo_resolution: if offset must be computed to fit density, \ - resolution of the temporary TPMS used to compute the offset - - :return: VTK PolyData object of the required TPMS part + :param type_part: ``"sheet"``, ``"lower skeletal"``, ``"upper skeletal"`` + or ``"surface"`` + :param algo_resolution: temporary-TPMS resolution for density→offset + search (only used when ``self.density`` is set) """ if type_part == "surface": return self.surface @@ -715,27 +811,29 @@ def generate_vtk( ) raise ValueError(err_msg) + if self.density == 1.0: + # Short-circuit: at full density the part fills the whole cell, + # return the cell-box mesh directly (exact volume — the + # marching-cubes path would give a slight discretization gap). + box = pv.Box(bounds=self._bounds, level=0, quads=False) + box = box.extract_surface().clean().triangulate() + box = rotate(box, center=(0, 0, 0), rotation=self.orientation) + return box.translate(xyz=self.center) + if self.density is not None: self._compute_offset_to_fit_density( part_type=type_part, resolution=algo_resolution, ) - self._check_offset(type_part) - - if type_part == "sheet": - part = self.as_sheet() - elif type_part == "upper skeletal": - part = self.as_upper_skeletal() - else: # lower skeletal - part = self.as_lower_skeletal() - - # Match the grid resolution: resolution * repeat_cell per axis. - # Shape.generate_vtk takes a single isotropic resolution, so use the - # geometric mean to keep total point count proportional. - iso_res = int( - self.resolution * np.cbrt(np.prod(self.repeat_cell)) + # Note: no `_check_offset` here — the F-rep VTK path handles negative + # / zero / variable offsets gracefully. ``generate()`` still applies + # the historical CAD-side restriction. + + frep = self._frep_part(type_part) + polydata = frep.generate_vtk( + bounds=self._bounds, + resolution=self._isotropic_resolution(), ) - polydata = part.generate_vtk(bounds=self._bounds, resolution=max(iso_res, 10)) polydata = rotate(polydata, center=(0, 0, 0), rotation=self.orientation) return polydata.translate(xyz=self.center) diff --git a/tests/shapes/test_tpms.py b/tests/shapes/test_tpms.py index 0d963fce..8b0ed9b6 100644 --- a/tests/shapes/test_tpms.py +++ b/tests/shapes/test_tpms.py @@ -3,7 +3,7 @@ from __future__ import annotations import re -from inspect import getmembers, isfunction +from inspect import getmembers, isfunction, signature from typing import Literal import numpy as np @@ -21,11 +21,18 @@ def _get_microgen_surface_functions() -> list[str]: - # Dont take into account deprecated surface functions named in camelCase + """List the actual TPMS surface functions in microgen.surface_functions. + + Filters out: + - deprecated camelCase aliases (any uppercase letter) + - non-TPMS callables exposed via re-import (autograd ``cos``/``sin``) — + they take fewer than 3 args and would crash when invoked as ``f(x,y,z)``. + """ return [ - func[0] - for func in getmembers(microgen.surface_functions, isfunction) - if not any(ele.isupper() for ele in func[0]) + name + for name, fn in getmembers(microgen.surface_functions, isfunction) + if not any(c.isupper() for c in name) + and len(signature(fn).parameters) == 3 ] From 20201a928c211abfa7bbff6a926ad08dcdc341d7 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Sat, 25 Apr 2026 22:25:18 +0200 Subject: [PATCH 07/11] Add Conformal TPMS and bunny examples Introduce a Conformal TPMS implementation and related improvements, plus two example scripts. Changes: - Add Conformal class to microgen.shape.tpms: builds TPMS in a local (w,u,v) frame using signed-distance as the radial coord, supports default_tangent_axis and optional clipping to the envelope. Overrides grid creation and cell-box to use the envelope surface. - Adjust Tpms density handling: density computations now measure against an envelope volume via _density_envelope_volume(), avoid cloning instances during root-finding, and correctly restore density after computing offset. - Infill now measures density relative to the provided object volume and overrides the cell-box to clip to the object envelope. - Add examples/conformal_vs_cartesian_bunny.py demonstrating side-by-side Cartesian vs Conformal gyroid infill on the Stanford bunny. - Add gyroid_gd_bunny.py utility script for visualizing a gyroid field clipped to a bunny and computing relative density. Rationale: provide a conformal TPMS mode that follows an arbitrary envelope surface (onion-shell style cells), make density semantics consistent (relative to the intended envelope), and ship examples to illustrate the difference. --- examples/conformal_vs_cartesian_bunny.py | 102 ++++++++ microgen/shape/tpms.py | 311 +++++++++++++++++++++-- 2 files changed, 394 insertions(+), 19 deletions(-) create mode 100644 examples/conformal_vs_cartesian_bunny.py diff --git a/examples/conformal_vs_cartesian_bunny.py b/examples/conformal_vs_cartesian_bunny.py new file mode 100644 index 00000000..d003638d --- /dev/null +++ b/examples/conformal_vs_cartesian_bunny.py @@ -0,0 +1,102 @@ +"""Side-by-side: Cartesian :class:`Infill` vs :class:`Conformal` on a bunny. + +Both fill the same bunny envelope with the same gyroid TPMS at the same +offset and ``repeat_cell``. Two visualizations: + +- **Left** — ``Infill``: the gyroid is in the world Cartesian frame, the bunny + surface just *clips* it. Cells are aligned to global axes; near the + envelope you see them sliced on arbitrary planes. +- **Right** — ``Conformal``: the gyroid's *radial* coord is the signed + distance to the bunny surface. Cells stack as concentric "shells" along + the surface normal — an onion-skin scaffold that always presents + perpendicular to the envelope. + +Run with ``python examples/conformal_vs_cartesian_bunny.py``. +""" + +from __future__ import annotations + +import numpy as np +import pyvista as pv +from pyvista import examples + +from microgen.shape.surface_functions import gyroid +from microgen.shape.tpms import Conformal, Infill + + +# ----------------------------------------------------------------------------- +# Build the envelope +# ----------------------------------------------------------------------------- + +bunny = examples.download_bunny() +# Original bunny is ~6cm; scale to a reasonable working size and recenter. +bunny.transform(np.diag([40.0, 40.0, 40.0, 1.0]), inplace=True) +bunny.translate(-np.array(bunny.center_of_mass()), inplace=True) +print(f"bunny: volume={bunny.volume:.3f} bounds={tuple(np.round(bunny.bounds, 2))}") + + +# ----------------------------------------------------------------------------- +# Two TPMS infills with identical parameters +# ----------------------------------------------------------------------------- + +OFFSET = 0.4 +REPEAT_CELL = 5 + +cartesian_infill = Infill( + obj=bunny, + surface_function=gyroid, + offset=OFFSET, + repeat_cell=REPEAT_CELL, +) +cartesian_sheet = cartesian_infill.generate_vtk(type_part="sheet") +print( + f"Cartesian Infill : sheet_volume={abs(cartesian_sheet.volume):.3f} " + f"density_vs_bunny={abs(cartesian_sheet.volume) / bunny.volume:.2%}", +) + +conformal = Conformal( + envelope=bunny, + surface_function=gyroid, + offset=OFFSET, + repeat_cell=REPEAT_CELL, + default_tangent_axis=(1.0, 0.0, 0.0), +) +conformal_sheet = conformal.generate_vtk(type_part="sheet") +print( + f"Conformal : sheet_volume={abs(conformal_sheet.volume):.3f} " + f"density_vs_bunny={abs(conformal_sheet.volume) / bunny.volume:.2%}", +) + + +# ----------------------------------------------------------------------------- +# Cut both meshes through y for a clean cross-section view +# ----------------------------------------------------------------------------- + +cut_origin = (0.0, -0.3 * float(bunny.bounds[3] - bunny.bounds[2]), 0.0) +cartesian_view = cartesian_sheet.clip("y", origin=cut_origin, invert=False) +conformal_view = conformal_sheet.clip("y", origin=cut_origin, invert=False) + + +# ----------------------------------------------------------------------------- +# Side-by-side plot +# ----------------------------------------------------------------------------- + +plotter = pv.Plotter(shape=(1, 2), window_size=(1600, 800)) + +plotter.subplot(0, 0) +plotter.add_text("Cartesian Infill\n(gyroid in world frame, bunny clips)", font_size=10) +plotter.add_mesh(bunny, color="white", opacity=0.15) +plotter.add_mesh(cartesian_view, color="orange") +plotter.view_isometric() + +plotter.subplot(0, 1) +plotter.add_text( + "Conformal\n(gyroid in surface-normal frame — onion shells)", + font_size=10, +) +plotter.add_mesh(bunny, color="white", opacity=0.15) +plotter.add_mesh(conformal_view, color="cyan") +plotter.view_isometric() + +plotter.link_views() +plotter.show() diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 6e6b9fa4..93297c12 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -177,6 +177,17 @@ def offset_from_density( density=density, )._compute_offset_to_fit_density(part_type=part_type, resolution=resolution) + def _density_envelope_volume(self: Tpms) -> float: + """Return the *envelope* volume that ``density`` is measured against. + + For a plain :class:`Tpms` this is the cartesian cell volume. Subclasses + with a non-cartesian envelope (:class:`Infill`, :class:`Conformal`) + override this to use ``obj.volume`` / ``envelope.volume``, so that + ``density`` consistently means *fraction of the envelope filled by the + TPMS part* — not fraction of the cartesian bounding box. + """ + return abs(self.grid.volume) + def _compute_offset_to_fit_density( self: Tpms, part_type: Literal["sheet", "lower skeletal", "upper skeletal"], @@ -190,6 +201,10 @@ def _compute_offset_to_fit_density( too high to reach at this resolution (marching cubes saturates slightly below 1.0 due to surface-grid discretization), falls back to the offset limit instead of failing the bracket. + + Density is measured against :meth:`_density_envelope_volume`, so that + :class:`Infill` / :class:`Conformal` density is relative to the input + envelope volume rather than the cartesian bbox. """ if self.density is None: err_msg = f"density must be between 0 and 1. Given: {self.density}" @@ -204,33 +219,36 @@ def _compute_offset_to_fit_density( ) return self.offset - temp_tpms = Tpms( - surface_function=self.surface_function, - offset=0.0, - cell_size=self.cell_size, - repeat_cell=self.repeat_cell, - resolution=resolution if resolution is not None else self.resolution, - ) - cell_volume = abs(temp_tpms.grid.volume) + # Density is measured directly on ``self`` so subclass envelopes are + # respected without having to clone the instance. ``self.density`` + # is cleared during the search so that ``generate_vtk`` does NOT + # recurse back into this method. The final ``computed_offset`` is + # set as the permanent offset before returning (and ``self.density`` + # is restored so that downstream callers can re-query it). + envelope_volume = self._density_envelope_volume() + target_density = self.density + self.density = None def density(offset: float) -> float: - temp_tpms.offset = offset - mesh = temp_tpms.generate_vtk(type_part=part_type) - return abs(mesh.volume) / cell_volume + self.offset = offset # setter; also clears _offset_func + mesh = self.generate_vtk(type_part=part_type) + return abs(mesh.volume) / envelope_volume bracket = self.offset_lim[part] try: computed_offset = root_scalar( - lambda offset: density(offset) - self.density, + lambda offset: density(offset) - target_density, bracket=bracket, ).root except ValueError: # Bracket sign mismatch — usually because target density exceeds - # the achievable max at this resolution. Pick the bracket - # endpoint that lies on the side of self.density. - d_lo = density(bracket[0]) - self.density - d_hi = density(bracket[1]) - self.density + # the achievable max at this resolution (or the envelope shape). + # Pick the bracket endpoint that lies on the side of target_density. + d_lo = density(bracket[0]) - target_density + d_hi = density(bracket[1]) - target_density computed_offset = bracket[1] if abs(d_hi) <= abs(d_lo) else bracket[0] + finally: + self.density = target_density self.offset = computed_offset return computed_offset @@ -1201,8 +1219,263 @@ def _create_grid( logging.info("Grid resolution: %s points", grid.n_points) return grid + def _density_envelope_volume(self: Infill) -> float: + """Density is measured against the input object's volume.""" + return abs(self.obj.volume) -# Re-export for backward compatibility -from .shape import ShellCreationError # noqa: E402 + def _cell_box(self: Infill) -> Shape: + """Override the cell-box clip with the *object envelope* SDF. + + Without this override the F-rep marching cubes path (used by + :meth:`Tpms.generate_vtk` and :meth:`Tpms.generate`) would clip the + TPMS to the cartesian bounding box rather than to the input object — + producing volumes well above ``obj.volume`` and corrupting density. + """ + from .implicit_ops import from_field # noqa: PLC0415 + + bounds_arr = np.array(self.obj.bounds) + bounds: BoundsType = tuple(bounds_arr.tolist()) + + envelope = self.obj + + def _envelope_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + pts = pv.PolyData(np.column_stack([x.ravel(), y.ravel(), z.ravel()])) + pts.compute_implicit_distance(envelope, inplace=True) + return np.asarray(pts["implicit_distance"]).reshape(x.shape) + + return from_field(_envelope_sdf, bounds=bounds) + + +class Conformal(Tpms): + """TPMS that conforms to the local frame of an envelope surface. + + Where :class:`CylindricalTpms` and :class:`SphericalTpms` use a closed-form + inverse coordinate map to wrap the TPMS around a known surface (cylinder / + sphere), :class:`Conformal` does the same for *any* envelope surface + (PolyData) by using the signed distance to the envelope as the radial + coordinate. + + Local frame at every Cartesian point ``p``: + + - ``w(p) = signed_distance(p, envelope)`` — *radial* coordinate; ``w=0`` on + the envelope surface, ``w<0`` inside, ``w>0`` outside. Computed by + :meth:`pyvista.PolyData.compute_implicit_distance`. + - ``u(p) = (p - center) · t``, ``v(p) = (p - center) · b`` — *tangential* + coordinates expressed in the global frame ``(t, b, n=default_axis)``, + where ``t`` and ``b`` are the two axes of the plane perpendicular to + ``default_tangent_axis``. Default tangent axis is ``(1, 0, 0)``. + + The TPMS field is then ``surface_function(k_w · w, k_u · u, k_v · v)``, + so unit cells stack along the envelope normal direction (concentric + "shells" at distance ``±cell_size``, ``±2·cell_size``, …) and tile the + perpendicular plane in the other two axes. Compared to a Cartesian + :class:`Infill`, this produces shells that *follow* the envelope shape + rather than being sliced by it. + + Tangentially-intrinsic (true conformal) parametrizations require + per-point local-tangent rotation; that's a follow-up — this class is the + minimum useful step. + """ + + def __init__( # noqa: PLR0913 + self: Conformal, + envelope: pv.PolyData, + surface_function: Field, + offset: float | OffsetGrading | Field | None = None, + cell_size: float | Sequence[float] | npt.NDArray[np.float64] | None = None, + repeat_cell: int | Sequence[int] | npt.NDArray[np.int8] | None = None, + phase_shift: Sequence[float] = (0.0, 0.0, 0.0), + resolution: int = 20, + density: float | None = None, + default_tangent_axis: Sequence[float] = (1.0, 0.0, 0.0), + clip_to_envelope: bool = True, + ) -> None: + """Build a TPMS that wraps an envelope surface. + + :param envelope: the surface (``pv.PolyData``) the TPMS will follow. + Volumes are computed relative to the envelope volume when + ``clip_to_envelope=True``. + :param surface_function: tpms function or custom function ``f(x,y,z)=0`` + :param offset: offset / thickness of the TPMS sheet + :param cell_size: unit cell size; auto-derived from envelope bounds if None + :param repeat_cell: number of cells per axis; auto-derived if None + :param phase_shift: phase shift in the (w, u, v) local frame + :param resolution: per-axis grid resolution + :param density: target density relative to envelope volume (mutex with + ``offset``) + :param default_tangent_axis: world axis used to define the tangential + plane for ``(u, v)``. Default ``(1, 0, 0)`` ⇒ ``u = y - cy``, + ``v = z - cz``. + :param clip_to_envelope: if ``True``, the resulting grid / mesh is + clipped by the envelope surface (TPMS only fills the interior). + """ + self.envelope = envelope + self.clip_to_envelope = clip_to_envelope + + # Frame: n = default_tangent_axis (radial-projection axis), then build + # an orthonormal (t, b) basis for the perpendicular plane. + n = np.asarray(default_tangent_axis, dtype=np.float64) + n_norm = np.linalg.norm(n) + if n_norm < 1e-12: + err_msg = "default_tangent_axis must be a non-zero vector" + raise ValueError(err_msg) + n /= n_norm + # Pick a stable secondary direction not parallel to n. + helper = np.array([0.0, 0.0, 1.0]) if abs(n[2]) < 0.9 else np.array([1.0, 0.0, 0.0]) + t = helper - np.dot(helper, n) * n + t /= np.linalg.norm(t) + b = np.cross(n, t) + self._frame_n = n + self._frame_t = t + self._frame_b = b + self.default_tangent_axis = tuple(n.tolist()) + + # Auto-derive cell_size / repeat_cell from envelope bbox if missing. + bounds = np.array(envelope.bounds) + margin = 1.001 + envelope_dim = margin * (bounds[1::2] - bounds[::2]) + + if cell_size is not None and repeat_cell is not None: + err_msg = ( + "cell_size and repeat_cell cannot be given at the same time, " + "one is computed from the other." + ) + raise ValueError(err_msg) + if cell_size is None and repeat_cell is None: + repeat_cell = (4, 4, 4) + if cell_size is not None: + repeat_cell = np.maximum( + np.round(envelope_dim / np.asarray(cell_size, dtype=float)).astype(int), + 1, + ) + elif repeat_cell is not None: + cell_size = envelope_dim / np.asarray(repeat_cell, dtype=float) + + self._envelope_center = np.asarray(envelope.center, dtype=np.float64) -__all__ = ["CylindricalTpms", "Infill", "ShellCreationError", "SphericalTpms", "Tpms"] + super().__init__( + surface_function=surface_function, + offset=offset, + phase_shift=phase_shift, + cell_size=cell_size, + repeat_cell=repeat_cell, + resolution=resolution, + density=density, + ) + + # -- Local-frame field ------------------------------------------------- + + def _envelope_distance( + self: Conformal, + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + """Vectorized signed distance to the envelope (negative inside).""" + pts = pv.PolyData(np.column_stack([x.ravel(), y.ravel(), z.ravel()])) + pts.compute_implicit_distance(self.envelope, inplace=True) + return np.asarray(pts["implicit_distance"]).reshape(x.shape) + + def _local_coords( + self: Conformal, + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> tuple[ + npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64] + ]: + """Project Cartesian ``(x, y, z)`` onto the conformal local frame. + + - ``w`` = signed distance to the envelope (radial coord) + - ``u`` = ``(p - envelope_center) · t`` + - ``v`` = ``(p - envelope_center) · b`` + """ + cx, cy, cz = self._envelope_center + dx, dy, dz = x - cx, y - cy, z - cz + u = dx * self._frame_t[0] + dy * self._frame_t[1] + dz * self._frame_t[2] + v = dx * self._frame_b[0] + dy * self._frame_b[1] + dz * self._frame_b[2] + w = self._envelope_distance(x, y, z) + return u, v, w + + def _setup_frep_field(self: Conformal) -> None: + """Build F-rep field that evaluates the surface function in (w, u, v).""" + k_x, k_y, k_z = 2.0 * np.pi / self.cell_size + ps = self.phase_shift + + def _raw_field( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + u, v, w = self._local_coords(x, y, z) + return self.surface_function( + k_x * (w + ps[0]), + k_y * (u + ps[1]), + k_z * (v + ps[2]), + ) + + bounds_arr = np.array(self.envelope.bounds) + bounds: BoundsType = tuple(bounds_arr.tolist()) + self._finalize_frep(_raw_field, bounds) + + # -- Override grid creation to clip to the envelope -------------------- + + def _create_grid( + self: Conformal, + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> pv.StructuredGrid: + """Cartesian grid covering the envelope bbox, optionally clipped.""" + cx, cy, cz = self._envelope_center + grid = super()._create_grid(x + cx, y + cy, z + cz) + if self.clip_to_envelope: + grid = grid.clip_surface(self.envelope) + logging.info( + "Conformal grid clipped to envelope: %d points", grid.n_points, + ) + return grid + + def _density_envelope_volume(self: Conformal) -> float: + """Density is measured against the envelope volume.""" + return abs(self.envelope.volume) + + # -- Override the F-rep "cell box" with the envelope itself ------------ + + def _cell_box(self: Conformal) -> Shape: + """SDF Shape of the envelope (used to clip TPMS parts to the interior). + + Replaces the parent's axis-aligned cell-box clip — the natural cell of + a Conformal TPMS *is* the envelope. Returns a Shape whose field is + ``signed_distance(p, envelope)`` (negative inside). + """ + from .implicit_ops import from_field # noqa: PLC0415 + + bounds_arr = np.array(self.envelope.bounds) + bounds: BoundsType = tuple(bounds_arr.tolist()) + + def _envelope_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + return self._envelope_distance(x, y, z) + + return from_field(_envelope_sdf, bounds=bounds) + + +# Re-export for backward compatibility +from .shape import BoundsType, ShellCreationError # noqa: E402 + +__all__ = [ + "Conformal", + "CylindricalTpms", + "Infill", + "ShellCreationError", + "SphericalTpms", + "Tpms", +] From d1ec0d39f71dc4585341d72aac04bbbf9fea709f Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 27 Apr 2026 09:15:40 +0200 Subject: [PATCH 08/11] TPMS: parametric grids, Sweep and GradedInfill MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Major refactor and feature additions to TPMS generation and infill handling. Key changes: - Add parametric-grid generation path for Cylindrical, Spherical and a new Sweep class (TPMS along a curve) to avoid MC Cartesian sampling artefacts; introduce _uses_parametric_grid flag and seam merge logic. - Add Sweep: builds parallel-transport frames, interpolates curve frames, maps (s,r,θ) parametric grid to Cartesian, and provides tube SDF and resolution heuristics. - Fix angle/arc conversions for cylindrical/spherical grids (use y/radius and z/radius) and tighten F-rep bounds; implement proper cell_box SDFs for cylindrical/spherical shells and Sweep tube. - Improve Tpms.generate path: support full-density envelope mesh via _envelope_mesh_at_full_density/_via_cell_box and prefer parametric-grid extraction when applicable; sample callable offsets from StructuredGrid or raw points. - Infill improvements: auto-orient normals, preserve original input volume for density calculations, route sheet/upper/lower properties to generate_vtk, and use object bbox for F-rep bounds. - Add helper SDFs and utilities: wedge and double-cone SDFs, per-axis interp helper. - Add new example scripts (tpms_infill_gallery.py, gyroid_gd_bunny.py), remove conformal_vs_cartesian_bunny.py, and add new tests (test_implicit_cylinder.py, test_kevin.py) plus updates to tests/shapes/test_tpms.py. Rationale: fix sampling artefacts near poles/axes, make full-density shortcut geometrically correct for non-box envelopes, support TPMS swept along arbitrary curves, and make density-fitting for infill robust to non-manifold/oriented inputs. Files added/removed include examples and small demo scripts. --- .DS_Store | Bin 0 -> 6148 bytes examples/conformal_vs_cartesian_bunny.py | 102 --- examples/tpms_infill_gallery.py | 116 +++ gyroid_gd_bunny.py | 87 ++ microgen/shape/tpms.py | 1060 +++++++++++++++++----- tests/shapes/test_tpms.py | 116 +++ 6 files changed, 1166 insertions(+), 315 deletions(-) create mode 100644 .DS_Store delete mode 100644 examples/conformal_vs_cartesian_bunny.py create mode 100644 examples/tpms_infill_gallery.py create mode 100644 gyroid_gd_bunny.py diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6 GIT binary patch literal 6148 zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0 pv.DataSet: + """Cut a mesh by a y-plane at ``frac × bbox_y`` so internals are visible.""" + bb = np.array(mesh.bounds) + y_origin = bb[2] + frac * (bb[3] - bb[2]) + return mesh.clip("y", origin=(0.0, y_origin, 0.0), invert=False) + + +def report(label: str, mesh: pv.DataSet) -> None: + print(f"{label:32s} vol={abs(mesh.volume):7.3f} n_cells={mesh.n_cells}") + + +# ----------------------------------------------------------------------------- +# 1. SphericalTpms — gyroid wrapping a sphere (auto-fill θ + φ via repeat=0) +# ----------------------------------------------------------------------------- + +sph = SphericalTpms( + radius=SPHERE_RADIUS, + surface_function=gyroid, + offset=OFFSET, + cell_size=CELL_SIZE, + repeat_cell=(2, 0, 0), +) +m_sph = sph.generate_vtk(type_part="sheet") +clip_y(m_sph).save(OUT / "sphere_g4_tpms.vtk") +report("1. SphericalTpms (R=3)", m_sph) + +# ----------------------------------------------------------------------------- +# 2. CylindricalTpms — gyroid wrapping a cylinder (auto-fill θ via repeat=0) +# ----------------------------------------------------------------------------- + +cyl = CylindricalTpms( + radius=CYLINDER_RADIUS, + surface_function=gyroid, + offset=OFFSET, + cell_size=CELL_SIZE, + repeat_cell=(2, 0, int(CYLINDER_HEIGHT / CELL_SIZE)), +) +m_cyl = cyl.generate_vtk(type_part="sheet") +clip_y(m_cyl).save(OUT / "cylinder_g5_tpms.vtk") +report("2. CylindricalTpms (R=1.5)", m_cyl) + +# ----------------------------------------------------------------------------- +# 3. Sweep — gyroid along a helical curve (1.5 turns, height 6) +# ----------------------------------------------------------------------------- + +def helix(t: float) -> np.ndarray: + """Helix: 1.5 turns, radius 2, height 6.""" + theta = 2.0 * np.pi * 1.5 * t + return np.array([2.0 * np.cos(theta), 2.0 * np.sin(theta), 6.0 * (t - 0.5)]) + + +sw = Sweep( + curve_points=helix, + surface_function=gyroid, + radial_max=0.6, + offset=OFFSET, + cell_size=CELL_SIZE, + repeat_cell=(8, 1, 6), + n_curve_samples=200, +) +m_sw = sw.generate_vtk(type_part="sheet") +clip_y(m_sw).save(OUT / "sweep_g7_helix.vtk") +report("3. Sweep along helix", m_sw) + +print("\nFiles saved (clipped by y-plane), open in ParaView:") +for path in ( + OUT / "sphere_g4_tpms.vtk", + OUT / "cylinder_g5_tpms.vtk", + OUT / "sweep_g7_helix.vtk", +): + print(f" {path}") diff --git a/gyroid_gd_bunny.py b/gyroid_gd_bunny.py new file mode 100644 index 00000000..225ee417 --- /dev/null +++ b/gyroid_gd_bunny.py @@ -0,0 +1,87 @@ +from typing import List + +import numpy as np +import pyvista as pv +from pyvista import examples + + +def gyroid(x, y, z): + return np.sin(x) * np.cos(y) + np.sin(y) * np.cos(z) + np.sin(z) * np.cos(x) + + +repeat_cell = np.array([8, 8, 8]) +cell_size = np.array([1, 1, 1]) +center_offset = 0.5 +resolution = 15 + +offset = np.pi/2. + +linspaces: List[np.ndarray] = [] +for repeat_cell_axis, cell_size_axis in zip(repeat_cell, cell_size): + linspaces.append( + np.linspace( + -center_offset * cell_size_axis * repeat_cell_axis, + center_offset * cell_size_axis * repeat_cell_axis, + resolution * repeat_cell_axis, + ) + ) + +x, y, z = np.meshgrid(*linspaces) + +grid = pv.StructuredGrid(x, y, z) +kx, ky, kz = 2 * np.pi / cell_size + +surface_function = gyroid(kx * x, ky * y, kz * z) + +bunny = examples.download_bunny() +transform_matrix = np.array([[40, 0, 0, 0], + [0, 40, 0, 0], + [0, 0, 40, 0], + [0, 0, 0, 1]]) +bunny.transform(transform_matrix, inplace=True) +center = bunny.center_of_mass() +bunny.translate(-center, inplace=True) +print("newcenter = ", bunny.center_of_mass()) +print(bunny.bounds) +print(grid.bounds) + +grid.compute_implicit_distance(bunny, inplace=True) + +#normalize : +dist = -1.*grid['implicit_distance'] +dist[dist > 0] = 0 +dist_norm = (dist - min(dist)) / (max(dist) - min(dist)) +x_t = 0.5 +l = 0.2 +reg_func = 0.6 * (1. + np.tanh((dist_norm - x_t)/l)) - 0.2 + +print(min(dist)) +print(max(dist)) + +print(min(dist_norm)) +print(max(dist_norm)) + +print(min(reg_func)) +print(max(reg_func)) + +grid["lower_surface"] = (surface_function.ravel(order="F") - offset*reg_func) +grid["upper_surface"] = (surface_function.ravel(order="F") + offset*reg_func) +sheet = grid.clip_scalar(scalars="upper_surface", invert=False).clip_scalar( + scalars="lower_surface" +) +upper_skeletal = grid.clip_scalar(scalars="upper_surface") +lower_skeletal = grid.clip_scalar(scalars="lower_surface", invert=False) + +clipped = sheet.clip_surface(bunny, invert=False) +clipped.compute_implicit_distance(bunny, inplace=True) + +print(f"relative density = {clipped.volume / bunny.volume:.2%}") + +clipped2 = clipped.clip('y', origin = (0,-0.5,0.), invert=False) +clipped2["dist"] = -1.*clipped2["implicit_distance"] + +pl = pv.Plotter() +#pl.add_mesh(grid, color="b", opacity = 0.1) +pl.add_mesh(bunny, color="w", opacity = 0.1) +pl.add_mesh(clipped2, scalars='dist', cmap='inferno') +pl.show() diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 93297c12..72c61aed 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -40,6 +40,50 @@ _DIM = 3 +def _wedge_sdf_2d( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + half_angle: float, +) -> npt.NDArray[np.float64]: + """SDF of a 2D wedge centred on +X with half-aperture ``half_angle``. + + Negative inside, positive outside, zero on the bounding rays. + Z-independent — used as a clipping primitive for cylindrical / spherical + partial wraps. + """ + return np.cos(half_angle) * np.abs(y) - np.sin(half_angle) * x + + +def _double_cone_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + half_polar: float, +) -> npt.NDArray[np.float64]: + """SDF of a double cone around the Z axis with half-aperture ``half_polar``. + + Symmetric about the XY plane; points inside iff their angle from ±Z is + less than ``half_polar``. Used to clip spherical partial-θ coverage. + """ + rho_xy = np.sqrt(x * x + y * y) + return np.cos(half_polar) * rho_xy - np.sin(half_polar) * np.abs(z) + + +def _interp_along_curve( + query_s: npt.NDArray[np.float64], + curve_s: npt.NDArray[np.float64], + values: npt.NDArray[np.float64], +) -> npt.NDArray[np.float64]: + """Per-axis ``np.interp`` of an ``(M, 3)`` curve attribute at ``query_s``. + + Returns an ``(N, 3)`` array. Wraps the three-axis loop so callers don't + repeat the pattern. + """ + return np.column_stack( + [np.interp(query_s, curve_s, values[:, k]) for k in range(values.shape[1])], + ) + + class Tpms(Shape): """Triply Periodical Minimal Surfaces. @@ -585,7 +629,16 @@ def offset( # Keep the callable so the F-rep path can re-evaluate it on the # marching-cubes grid (which differs from ``self.grid``). self._offset_func = offset - self._offset = offset(self.grid.x, self.grid.y, self.grid.z).ravel("F") + # Sample the callable on ``self.grid`` for legacy / property-path + # consumers (``grid_sheet`` etc.). StructuredGrid exposes + # ``.x/.y/.z`` meshgrids; UnstructuredGrid (e.g. after + # ``clip_surface`` in :class:`Infill`) does not — use raw points. + if hasattr(self.grid, "x"): + vals = offset(self.grid.x, self.grid.y, self.grid.z).ravel("F") + else: + pts = np.asarray(self.grid.points) + vals = offset(pts[:, 0], pts[:, 1], pts[:, 2]) + self._offset = vals else: err_msg = "offset must be a float, a numpy array or a callable" raise TypeError(err_msg) @@ -696,6 +749,34 @@ def _isotropic_resolution(self: Tpms) -> int: """ return max(int(self.resolution * np.cbrt(np.prod(self.repeat_cell))), 10) + def _envelope_mesh_at_full_density(self: Tpms) -> pv.PolyData: + """Return the cell-envelope mesh used by the ``density=1.0`` shortcut. + + Plain :class:`Tpms` has an axis-aligned box envelope, returned exactly + via :class:`pyvista.Box`. Subclasses with non-cubic envelopes + (cylindrical / spherical shells, ``Infill`` object, etc.) override + this to fall back to marching cubes on their ``_cell_box`` SDF — see + :meth:`_envelope_mesh_via_cell_box`. + """ + return ( + pv.Box(bounds=self._bounds, level=0, quads=False) + .extract_surface() + .clean() + .triangulate() + ) + + def _envelope_mesh_via_cell_box(self: Tpms) -> pv.PolyData: + """Marching-cubes mesh of the (potentially non-cubic) ``_cell_box``. + + Used by subclasses that override :meth:`_envelope_mesh_at_full_density` + to delegate to the F-rep envelope SDF. + """ + envelope_shape = self._cell_box() + return envelope_shape.generate_vtk( + bounds=envelope_shape.bounds or self._bounds, + resolution=self._isotropic_resolution(), + ) + def generate( self: Tpms, type_part: TpmsPartType = "sheet", @@ -830,13 +911,11 @@ def generate_vtk( raise ValueError(err_msg) if self.density == 1.0: - # Short-circuit: at full density the part fills the whole cell, - # return the cell-box mesh directly (exact volume — the - # marching-cubes path would give a slight discretization gap). - box = pv.Box(bounds=self._bounds, level=0, quads=False) - box = box.extract_surface().clean().triangulate() - box = rotate(box, center=(0, 0, 0), rotation=self.orientation) - return box.translate(xyz=self.center) + envelope_mesh = self._envelope_mesh_at_full_density() + envelope_mesh = rotate( + envelope_mesh, center=(0, 0, 0), rotation=self.orientation, + ) + return envelope_mesh.translate(xyz=self.center) if self.density is not None: self._compute_offset_to_fit_density( @@ -847,11 +926,37 @@ def generate_vtk( # / zero / variable offsets gracefully. ``generate()`` still applies # the historical CAD-side restriction. - frep = self._frep_part(type_part) - polydata = frep.generate_vtk( - bounds=self._bounds, - resolution=self._isotropic_resolution(), - ) + # Subclasses with a *parametric* coordinate frame (CylindricalTpms, + # SphericalTpms, Sweep) cannot use the F-rep marching-cubes path — + # the TPMS field has rapid angular oscillations that an isotropic + # Cartesian grid samples poorly (visible as dotted holes near the + # pole / axis). They opt into the legacy grid-clip path by setting + # ``_uses_parametric_grid = True``: clip the parametric structured + # grid by the relevant scalar threshold, then extract the surface. + if getattr(self, "_uses_parametric_grid", False): + grid_attr = f"grid_{type_part.replace(' ', '_')}" + # Merge points within an absolute tolerance large enough to close + # the angular seam (φ=±π for sphere, θ=±π for cylinder) where + # the structured grid's two boundary faces coincide in Cartesian + # space up to ~1e-3 floating-point drift, but small enough not + # to collapse legitimate cell edges (smallest grid spacing is + # ~ ``cell_size / resolution``). + seam_tol = min( + 5e-3, + 0.1 * float(np.min(self.cell_size)) / float(self.resolution), + ) + polydata = ( + getattr(self, grid_attr) + .extract_surface() + .clean(tolerance=seam_tol, absolute=True) + .triangulate() + ) + else: + frep = self._frep_part(type_part) + polydata = frep.generate_vtk( + bounds=self._bounds, + resolution=self._isotropic_resolution(), + ) polydata = rotate(polydata, center=(0, 0, 0), rotation=self.orientation) return polydata.translate(xyz=self.center) @@ -894,6 +999,15 @@ def generateVtk( # noqa: N802 class CylindricalTpms(Tpms): """Class used to generate cylindrical TPMS geometries (sheet or skeletals parts).""" + # Use the parametric structured-grid clip for ``generate_vtk`` instead of + # F-rep marching cubes. An isotropic Cartesian MC grid samples the + # angular axis poorly near the cylinder axis (rapid θ-derivative ⇒ + # under-sampled holes). The structured grid in (ρ, θ, z) automatically + # maps to a Cartesian grid that is dense near the axis. + _uses_parametric_grid = True + + _envelope_mesh_at_full_density = Tpms._envelope_mesh_via_cell_box + def __init__( # noqa: PLR0913 self: CylindricalTpms, radius: float, @@ -965,9 +1079,16 @@ def _create_grid( y: npt.NDArray[np.float64], z: npt.NDArray[np.float64], ) -> pv.StructuredGrid: - """Return the structured cylindrical grid of the TPMS.""" + """Return the structured cylindrical grid of the TPMS. + + ``y`` carries arc-length along the equator, so the conversion to an + angle is ``theta = y / radius``. The earlier ``y * unit_theta`` form + over-wrapped past ``2π`` by a factor of ``unit_theta * radius`` + (which is only 1 when the user-supplied ``cell_size[1] == 1``), + leaving a visible wedge of unrendered geometry on one meridian. + """ rho = x + self.cylinder_radius - theta = y * self.unit_theta + theta = y / self.cylinder_radius grid = pv.StructuredGrid(rho * np.cos(theta), rho * np.sin(theta), z) grid["coords"] = np.c_[ @@ -977,6 +1098,15 @@ def _create_grid( ] return grid + def _shell_extents(self: CylindricalTpms) -> tuple[float, float, float]: + """Return ``(r_inner, r_outer, half_z)`` of the TPMS-bearing shell.""" + cyl_r = float(self.cylinder_radius) + delta_r = 0.5 * float(self.cell_size[0]) * float(self.repeat_cell[0]) + r_outer = cyl_r + delta_r + r_inner = max(cyl_r - delta_r, 0.0) + half_z = 0.5 * float(self.cell_size[2]) * float(self.repeat_cell[2]) + return r_inner, r_outer, half_z + def _setup_frep_field(self: CylindricalTpms) -> None: """Build F-rep field in Cartesian coordinates (inverse cylindrical mapping).""" k_x, k_y, k_z = 2.0 * np.pi / self.cell_size @@ -997,17 +1127,95 @@ def _raw_field( k_z * (z + ps[2]), ) - r_max = cyl_r + 0.5 * self.cell_size[0] * self.repeat_cell[0] - half_z = 0.5 * self.cell_size[2] * self.repeat_cell[2] + # Tightened bounds: only ``r_outer`` reach in xy, only ``half_z`` in z. + # The previous wider Cartesian box wasted ~50 % of MC samples on the + # interior of the cylinder where the TPMS doesn't live. + _, r_outer, half_z = self._shell_extents() self._finalize_frep( _raw_field, - (-r_max, r_max, -r_max, r_max, -half_z, half_z), + (-r_outer, r_outer, -r_outer, r_outer, -half_z, half_z), ) + def _cell_box(self: CylindricalTpms) -> Shape: + """Cylindrical-shell SDF in Cartesian space (replaces the parent's + axis-aligned-box SDF, which was geometrically wrong for this class). + + Without this override the F-rep marching cubes path + (:meth:`Tpms._clipped_sheet`) would clip the TPMS by an axis-aligned + box of ``cell_size × repeat_cell`` *interpreted as Cartesian* — but + ``cell_size[1]`` / ``repeat_cell[1]`` are angular (parametric) for + this class, so the clip discarded ~95 % of the physical shell. + + The shell SDF is:: + + radial = sqrt(x² + y²) + ring = max(r_inner − radial, radial − r_outer) + sdf = max(ring, |z| − half_z) + + Partial wrap (``repeat_cell[1] < n_repeat_to_full_circle``) is + handled by intersecting the ring with an angular wedge centred on + the +X axis. + """ + from .implicit_ops import from_field, intersection # noqa: PLC0415 + + r_inner, r_outer, half_z = self._shell_extents() + bounds: BoundsType = ( + -r_outer, r_outer, -r_outer, r_outer, -half_z, half_z, + ) + + def _shell_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + radial = np.sqrt(x * x + y * y) + outer = radial - r_outer + if r_inner > 0: + ring = np.maximum(r_inner - radial, outer) + else: + ring = outer # solid disc — no inner wall + return np.maximum(ring, np.abs(z) - half_z) + + shell = from_field(_shell_sdf, bounds=bounds) + + # Partial angular wrap → intersect with a wedge centred on +X. + n_full = int(round(2.0 * np.pi / self.unit_theta)) + if int(self.repeat_cell[1]) < n_full: + half_angle = np.pi * float(self.repeat_cell[1]) / float(n_full) + wedge = from_field( + lambda x, y, z, _h=half_angle: _wedge_sdf_2d(x, y, _h), + bounds=bounds, + ) + shell = intersection(shell, wedge) + + return shell + + def _isotropic_resolution(self: CylindricalTpms) -> int: + """Per-axis resolution count for the F-rep marching-cubes grid. + + Override the parent's geometric-mean of ``repeat_cell`` (which mixes + the angular axis with linear ones). Use only the radial (``[0]``) + and axial (``[2]``) physical extents, both in world units. + """ + radial_extent = float(self.cell_size[0]) * float(self.repeat_cell[0]) + axial_extent = float(self.cell_size[2]) * float(self.repeat_cell[2]) + cell_size_min = max( + min(float(self.cell_size[0]), float(self.cell_size[2])), 1e-12, + ) + n = int(self.resolution * max(radial_extent, axial_extent) / cell_size_min) + return max(n, 10) + class SphericalTpms(Tpms): """Class used to generate spherical TPMS geometries (sheet or skeletals parts).""" + # Same rationale as :class:`CylindricalTpms` — the parametric (ρ, θ, φ) + # grid stretches naturally in Cartesian space, sampling poles and equator + # equally well, while uniform Cartesian MC produces pole pathologies. + _uses_parametric_grid = True + + _envelope_mesh_at_full_density = Tpms._envelope_mesh_via_cell_box + def __init__( # noqa: PLR0913 self: SphericalTpms, radius: float, @@ -1085,10 +1293,19 @@ def _create_grid( y: npt.NDArray[np.float64], z: npt.NDArray[np.float64], ) -> pv.StructuredGrid: - """Return the structured spherical grid of the TPMS.""" + """Return the structured spherical grid of the TPMS. + + ``y`` and ``z`` carry equatorial arc length, so the conversion to + angles is ``theta = y / radius`` and ``phi = z / radius``. The + earlier ``* unit_theta`` / ``* unit_phi`` forms over-wrapped past + ``π`` / ``2π`` by a factor of ``unit_* * radius`` (which is only + 1 when the user-supplied ``cell_size == 1``), leaving a visible + wedge of unrendered geometry on one meridian and tiny over-shoots + at the poles. + """ rho = x + self.sphere_radius - theta = y * self.unit_theta + np.pi / 2.0 - phi = z * self.unit_phi + theta = y / self.sphere_radius + np.pi / 2.0 + phi = z / self.sphere_radius grid = pv.StructuredGrid( rho * np.sin(theta) * np.cos(phi), @@ -1128,16 +1345,490 @@ def _raw_field( k_z * (phi + ps[2]), ) - r_max = sph_r + 0.5 * self.cell_size[0] * self.repeat_cell[0] + # Tightened bounds: only ``r_outer`` cube, not ``r_max + R``. + _, r_outer = self._shell_radii() self._finalize_frep( _raw_field, - (-r_max, r_max, -r_max, r_max, -r_max, r_max), + (-r_outer, r_outer, -r_outer, r_outer, -r_outer, r_outer), + ) + + def _shell_radii(self: SphericalTpms) -> tuple[float, float]: + """Return ``(r_inner, r_outer)`` of the spherical TPMS-bearing shell.""" + sph_r = float(self.sphere_radius) + delta_r = 0.5 * float(self.cell_size[0]) * float(self.repeat_cell[0]) + return max(sph_r - delta_r, 0.0), sph_r + delta_r + + def _cell_box(self: SphericalTpms) -> Shape: + """Spherical-shell SDF in Cartesian space. + + Replaces the parent's axis-aligned-box clip — same fix as + :meth:`CylindricalTpms._cell_box`, see that docstring for the + rationale. + + Partial coverage in θ (cone clip on +Z) and φ (wedge clip on +X + in xy-plane) is supported when ``repeat_cell[1]`` / + ``repeat_cell[2]`` are smaller than the auto-fill values. + """ + from .implicit_ops import from_field, intersection # noqa: PLC0415 + + r_inner, r_outer = self._shell_radii() + bounds: BoundsType = ( + -r_outer, r_outer, -r_outer, r_outer, -r_outer, r_outer, + ) + + def _shell_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + r = np.sqrt(x * x + y * y + z * z) + outer = r - r_outer + if r_inner > 0: + return np.maximum(r_inner - r, outer) + return outer # solid ball if inner radius collapsed + + shape = from_field(_shell_sdf, bounds=bounds) + + # Partial θ coverage → double-cone clip around the Z axis. + n_full_theta = int(round(np.pi / self.unit_theta)) + if int(self.repeat_cell[1]) < n_full_theta: + half_polar = np.pi * float(self.repeat_cell[1]) / float(n_full_theta) + shape = intersection( + shape, + from_field( + lambda x, y, z, _h=half_polar: _double_cone_sdf(x, y, z, _h), + bounds=bounds, + ), + ) + + # Partial φ coverage → wedge clip around the +X axis. + n_full_phi = int(round(2.0 * np.pi / self.unit_phi)) + if int(self.repeat_cell[2]) < n_full_phi: + half_azi = np.pi * float(self.repeat_cell[2]) / float(n_full_phi) + shape = intersection( + shape, + from_field( + lambda x, y, z, _h=half_azi: _wedge_sdf_2d(x, y, _h), + bounds=bounds, + ), + ) + + return shape + + def _isotropic_resolution(self: SphericalTpms) -> int: + """Use only the radial physical extent (the angular axes are + parametric and would mislead a geometric mean across them). + """ + radial_extent = float(self.cell_size[0]) * float(self.repeat_cell[0]) + radius_envelope = self.sphere_radius + radial_extent + cell_size_radial = max(float(self.cell_size[0]), 1e-12) + n = int(self.resolution * 2.0 * radius_envelope / cell_size_radial) + return max(n, 10) + + +class Sweep(Tpms): + """TPMS along an arbitrary 3D curve — generalisation of CylindricalTpms. + + The TPMS is generated in a tube of radius ``radial_max`` around an + arbitrary curve, in local coordinates ``(s, r, θ)`` where: + + - ``s`` = arc length along the curve from its start + - ``r`` = perpendicular distance from the curve at the closest point + - ``θ`` = angle around the curve in the local tangent-frame plane + + The local frame ``(T, N, B)`` is built once at curve setup time using + parallel transport along the polyline (tangent from finite differences, + normal from a seed transported by Rodrigues rotation between adjacent + samples; closed loops get a holonomy correction distributed linearly). + + The mesh is generated via the **parametric structured-grid path** + (same as :class:`CylindricalTpms` / :class:`SphericalTpms`): we build a + structured grid in (s, r, θ) parametric space, evaluate the TPMS field + on that grid, and map the parametric points to Cartesian using the + parallel-transport frames. ``generate_vtk`` then clips the structured + grid by the relevant scalar threshold — much cleaner than F-rep MC at + a uniform Cartesian resolution, which would under-sample the angular + direction and produce dotted artefacts. + + :class:`CylindricalTpms` is a special case of this class (curve = a + straight line). + """ + + # Use the parametric grid-clip path (same as Cylindrical/Spherical). + _uses_parametric_grid = True + + _envelope_mesh_at_full_density = Tpms._envelope_mesh_via_cell_box + + def __init__( # noqa: PLR0913 + self: Sweep, + curve_points: npt.NDArray[np.float64] | Callable[[float], npt.NDArray[np.float64]], + surface_function: Field, + radial_max: float, + offset: float | OffsetGrading | Field | None = None, + cell_size: float | Sequence[float] | npt.NDArray[np.float64] = 1.0, + repeat_cell: int | Sequence[int] | npt.NDArray[np.int8] = 1, + phase_shift: Sequence[float] = (0.0, 0.0, 0.0), + resolution: int = 20, + density: float | None = None, + seed_normal: Sequence[float] | None = None, + n_curve_samples: int = 200, + center: Vector3DType = (0, 0, 0), + orientation: Vector3DType = (0, 0, 0), + ) -> None: + r"""Build a TPMS swept along a curve. + + :param curve_points: either an ``(M, 3)`` array of polyline samples + or a callable ``t \in [0, 1] -> (3,)``. Callables are sampled + at ``n_curve_samples`` points before processing. + :param surface_function: TPMS function ``f(x, y, z)`` + :param radial_max: outer tube radius + :param offset: TPMS sheet thickness + :param cell_size: ``(s, r, θ)`` cell size — third axis is angular + (radians per cell), so a sensible default is to leave it at + ``1.0`` and tune via ``repeat_cell[2]``. + :param repeat_cell: ``(n_s, n_r, n_θ)`` — radial cell count is + usually ``1`` for a thin tube; ``n_θ`` controls how many + angular cells around the curve. + :param phase_shift: TPMS phase shift in (s, r, θ) + :param resolution: per-axis MC grid resolution + :param density: mutex with ``offset``; density relative to the + tube volume + :param seed_normal: initial normal direction for parallel transport; + defaults to a vector perpendicular to the first tangent + :param n_curve_samples: number of samples to use when resampling a + ``Callable`` curve (ignored for polyline input) + :param center: center of the geometry + :param orientation: orientation of the geometry + """ + # Discretise the curve. + if callable(curve_points): + ts = np.linspace(0.0, 1.0, int(n_curve_samples)) + curve = np.asarray([curve_points(t) for t in ts], dtype=np.float64) + else: + curve = np.asarray(curve_points, dtype=np.float64) + if curve.ndim != 2 or curve.shape[1] != 3 or curve.shape[0] < 2: + err_msg = ( + f"curve_points must be an (M, 3) array with M ≥ 2, " + f"got shape {curve.shape}" + ) + raise ValueError(err_msg) + + self.curve = curve + self.radial_max = float(radial_max) + + # Build the parallel-transport frames + arc-length parametrisation. + self._build_curve_frames(seed_normal=seed_normal) + + # Initialise like a regular TPMS — cell_size now lives in (s, r, θ) + # parametric space, ``_setup_frep_field`` will use the local frames. + super().__init__( + surface_function=surface_function, + offset=offset, + phase_shift=phase_shift, + cell_size=cell_size, + repeat_cell=repeat_cell, + resolution=resolution, + density=density, + center=center, + orientation=orientation, + ) + + # -- Curve preprocessing ----------------------------------------------- + + def _build_curve_frames( + self: Sweep, + seed_normal: Sequence[float] | None, + ) -> None: + """Pre-compute arc length, tangents, and parallel-transported normals. + + Stores on ``self``: + ``_curve_s`` arc length per sample, shape (M,) + ``_curve_T`` unit tangent per sample, shape (M, 3) + ``_curve_N`` unit normal (parallel-transported), shape (M, 3) + ``_curve_kdtree`` scipy cKDTree on the curve points + """ + from scipy.spatial import cKDTree # noqa: PLC0415 + + pts = self.curve + m = pts.shape[0] + + # Arc-length parameterisation. + seg = np.linalg.norm(np.diff(pts, axis=0), axis=1) + s_cum = np.concatenate([[0.0], np.cumsum(seg)]) + self._curve_s = s_cum + self._curve_total_length = float(s_cum[-1]) + + # Tangent: forward differences, last point copies its predecessor. + diffs = np.diff(pts, axis=0) + diff_norm = np.linalg.norm(diffs, axis=1, keepdims=True) + diff_norm = np.where(diff_norm > 1e-12, diff_norm, 1.0) + T = diffs / diff_norm + T = np.vstack([T, T[-1:]]) # pad final tangent + self._curve_T = T + + # Seed normal: project the user's seed (or world-up / world-x as + # fallback) onto the plane perpendicular to T[0]. + if seed_normal is None: + world_up = np.array([0.0, 0.0, 1.0]) + seed_proj = world_up - np.dot(world_up, T[0]) * T[0] + if np.linalg.norm(seed_proj) < 1e-6: + # T[0] is along z — use world-x instead. + world_x = np.array([1.0, 0.0, 0.0]) + seed_proj = world_x - np.dot(world_x, T[0]) * T[0] + else: + seed_proj = np.asarray(seed_normal, dtype=np.float64) + seed_proj = seed_proj - np.dot(seed_proj, T[0]) * T[0] + if np.linalg.norm(seed_proj) < 1e-6: + err_msg = "seed_normal is parallel to the initial tangent" + raise ValueError(err_msg) + N = np.empty_like(T) + N[0] = seed_proj / np.linalg.norm(seed_proj) + + # Parallel transport: rotate previous N around (T_{i-1} × T_i) by + # the angle between successive tangents (Rodrigues' rotation). + for i in range(1, m): + t_prev, t_curr = T[i - 1], T[i] + cos_angle = float(np.clip(np.dot(t_prev, t_curr), -1.0, 1.0)) + if cos_angle > 1.0 - 1e-9: + N[i] = N[i - 1] + continue + axis = np.cross(t_prev, t_curr) + axis_n = np.linalg.norm(axis) + if axis_n < 1e-9: + # Antiparallel tangents (curve folds back) — keep N as-is. + N[i] = N[i - 1] + continue + axis = axis / axis_n + sin_angle = np.sqrt(max(1.0 - cos_angle * cos_angle, 0.0)) + # Rodrigues: v' = v cos θ + (k × v) sin θ + k (k·v) (1 − cos θ) + v = N[i - 1] + v_rot = ( + v * cos_angle + + np.cross(axis, v) * sin_angle + + axis * np.dot(axis, v) * (1.0 - cos_angle) + ) + v_rot = v_rot - np.dot(v_rot, t_curr) * t_curr # re-project + v_rot = v_rot / np.linalg.norm(v_rot) + N[i] = v_rot + + # Closed-loop holonomy correction: if the curve closes, distribute + # the residual twist between N[-1] (transported) and the seed N[0]. + is_closed = bool(np.linalg.norm(pts[-1] - pts[0]) < 1e-6) + if is_closed: + holonomy = float( + np.arctan2( + float(np.dot(np.cross(N[-1], N[0]), T[-1])), + float(np.dot(N[-1], N[0])), + ), + ) + # Vectorised Rodrigues de-twist: each sample's N rotates around + # its own T by an angle proportional to its arc-length fraction. + ang = -holonomy * self._curve_s / max(self._curve_total_length, 1e-12) + cos_a = np.cos(ang)[:, None] + sin_a = np.sin(ang)[:, None] + kdotv = np.einsum("ij,ij->i", T, N)[:, None] + N_rot = N * cos_a + np.cross(T, N) * sin_a + T * kdotv * (1.0 - cos_a) + N = N_rot / np.linalg.norm(N_rot, axis=1, keepdims=True) + self._curve_N = N + + self._curve_kdtree = cKDTree(pts) + + # Self-intersection guard. + if m >= 3: + d2_dt2 = np.diff(T, axis=0) + kappa = np.linalg.norm(d2_dt2, axis=1) / np.maximum(np.diff(s_cum), 1e-12) + kappa_max = float(np.max(kappa)) if kappa.size else 0.0 + if kappa_max * self.radial_max >= 1.0: + logging.warning( + "Sweep: max curvature %.3f × radial_max %.3f ≥ 1 — the " + "tube envelope self-intersects in some segments; the " + "TPMS may be ill-defined there.", + kappa_max, + self.radial_max, + ) + + # -- Local-frame field ------------------------------------------------- + + def _local_coords( + self: Sweep, + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> tuple[ + npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64], + ]: + """Return ``(s, r, θ)`` for each input point ``p = (x, y, z)``.""" + shape = x.shape + pts = np.column_stack([x.ravel(), y.ravel(), z.ravel()]) + _dist, idx = self._curve_kdtree.query(pts, k=1) + + c = self.curve[idx] + T = self._curve_T[idx] + N = self._curve_N[idx] + B = np.cross(T, N) + + delta = pts - c + s = self._curve_s[idx] + # signed distance along N and B → radial r and angle θ. + u = (delta * N).sum(axis=1) + v = (delta * B).sum(axis=1) + r = np.sqrt(u * u + v * v) + theta = np.arctan2(v, u) + return s.reshape(shape), r.reshape(shape), theta.reshape(shape) + + def _setup_frep_field(self: Sweep) -> None: + """Build the F-rep field as ``surface_function(k_s s, k_r r, k_θ θ)``.""" + k_x, k_y, k_z = 2.0 * np.pi / self.cell_size + ps = self.phase_shift + + def _raw_field( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + s, r, theta = self._local_coords(x, y, z) + return self.surface_function( + k_x * (s + ps[0]), + k_y * (r + ps[1]), + k_z * (theta + ps[2]), + ) + + # Bounds: curve AABB inflated by ``radial_max``. + bb_min = self.curve.min(axis=0) - self.radial_max + bb_max = self.curve.max(axis=0) + self.radial_max + bounds: BoundsType = ( + float(bb_min[0]), float(bb_max[0]), + float(bb_min[1]), float(bb_max[1]), + float(bb_min[2]), float(bb_max[2]), + ) + + # Like Conformal, the field is built around discrete data — skip + # autograd SDF normalisation (it would FD-fall-back and be slow). + self._raw_field_func = _raw_field + self._func = _raw_field + self._bounds = bounds + + # -- Cell-box (tube SDF) ----------------------------------------------- + + def _cell_box(self: Sweep) -> Shape: + """Tube SDF: ``dist_to_curve(p) − radial_max``.""" + from .implicit_ops import from_field # noqa: PLC0415 + + bounds: BoundsType = ( + self._bounds if self._bounds is not None else (-1.0, 1.0, -1.0, 1.0, -1.0, 1.0) + ) + + def _tube_sdf( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + shape = x.shape + pts = np.column_stack([x.ravel(), y.ravel(), z.ravel()]) + d, _ = self._curve_kdtree.query(pts, k=1) + return (d - self.radial_max).reshape(shape) + + return from_field(_tube_sdf, bounds=bounds) + + def _isotropic_resolution(self: Sweep) -> int: + """Use the curve length and tube diameter as the spatial extents.""" + n = int( + self.resolution + * max(self._curve_total_length, 2.0 * self.radial_max) + / max(float(self.cell_size[1]), 1e-12), ) + return max(n, 10) + + def _create_grid( + self: Sweep, + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> pv.StructuredGrid: + """Map a parametric (s, r, θ) grid to Cartesian via the curve frames. + + ``x``, ``y``, ``z`` here are the parametric meshgrid coordinates + from :meth:`Tpms._compute_tpms_field` (they live in (s, r, θ) + space, *not* Cartesian). We interpolate the curve point ``c(s)`` + and frame ``(N(s), B(s))`` for each parametric ``s`` value, then + compute Cartesian ``p = c(s) + r·(cos θ · N + sin θ · B)``. + + The resulting ``StructuredGrid`` has Cartesian point coordinates + but retains the parametric ordering, so ``self.grid["surface"]`` + (filled by the parent :meth:`Tpms._compute_tpms_field`) can be + clipped by scalar threshold to extract the TPMS sheet/skeletal — + same recipe as :class:`CylindricalTpms`. + """ + # x, y, z come from numpy.meshgrid of parametric linspaces. Their + # raw values span [-0.5*cell_size*repeat_cell, +0.5*…] per axis. + # Map to ``s ∈ [0, total_length]``, ``r ∈ [0, radial_max]``, + # ``θ ∈ [0, 2π]``. + s_extent = float(self.cell_size[0]) * float(self.repeat_cell[0]) + s = (x - x.min()) / max(s_extent, 1e-12) * self._curve_total_length + # Radial: parametric coord shifted to [0, cell_size_r * repeat_cell_r]. + r = y - y.min() + # Angular: parametric coord scaled to [0, 2π] over the angular + # extent (cell_size[2] * repeat_cell[2]). + ang_extent = float(self.cell_size[2]) * float(self.repeat_cell[2]) + theta = (z - z.min()) / max(ang_extent, 1e-12) * 2.0 * np.pi + + # Interpolate curve position + (N, T) frames at each parametric s. + s_flat = s.ravel(order="F") + s_clipped = np.clip(s_flat, 0.0, self._curve_total_length) + c, N, T = ( + _interp_along_curve(s_clipped, self._curve_s, arr) + for arr in (self.curve, self._curve_N, self._curve_T) + ) + N /= np.maximum(np.linalg.norm(N, axis=1, keepdims=True), 1e-12) + T /= np.maximum(np.linalg.norm(T, axis=1, keepdims=True), 1e-12) + B = np.cross(T, N) + + r_flat = r.ravel(order="F") + theta_flat = theta.ravel(order="F") + cos_t = np.cos(theta_flat) + sin_t = np.sin(theta_flat) + + cart = c + r_flat[:, None] * (cos_t[:, None] * N + sin_t[:, None] * B) + cart_x = cart[:, 0].reshape(x.shape, order="F") + cart_y = cart[:, 1].reshape(x.shape, order="F") + cart_z = cart[:, 2].reshape(x.shape, order="F") + + grid = pv.StructuredGrid(cart_x, cart_y, cart_z) + grid["coords"] = np.c_[ + x.ravel(order="F"), + y.ravel(order="F"), + z.ravel(order="F"), + ] + return grid class Infill(Tpms): """Generate a TPMS infill inside a given object.""" + _envelope_mesh_at_full_density = Tpms._envelope_mesh_via_cell_box + + # The parent's ``sheet`` / ``upper_skeletal`` / ``lower_skeletal`` + # properties read from ``self.grid_*`` (legacy grid-clip path), but the + # density-fitting search optimises against :meth:`generate_vtk` (F-rep + # marching cubes clipped to the obj envelope) — for small infill objects + # the two paths discretise to noticeably different volumes. Re-route + # these properties to ``generate_vtk`` so ``density=0.5`` reliably gives + # ``infill.sheet.volume ≈ 0.5 * obj.volume``. + @property + def sheet(self: Infill) -> pv.PolyData: + """Sheet part as a PolyData mesh (uses :meth:`generate_vtk`).""" + return self.generate_vtk(type_part="sheet") + + @property + def upper_skeletal(self: Infill) -> pv.PolyData: + """Upper-skeletal part as a PolyData mesh.""" + return self.generate_vtk(type_part="upper skeletal") + + @property + def lower_skeletal(self: Infill) -> pv.PolyData: + """Lower-skeletal part as a PolyData mesh.""" + return self.generate_vtk(type_part="lower skeletal") + def __init__( # noqa: PLR0913 self: Infill, obj: pv.PolyData, @@ -1168,8 +1859,20 @@ def __init__( # noqa: PLR0913 If density is given, the offset is automatically computed to fit the\ density (performance is slower than when using the offset) """ - self.obj = obj - bounds = np.array(obj.bounds) + # Capture the original volume *before* re-orienting normals — for + # non-manifold inputs (e.g. pyvista's caps-less Cylinder, the + # Stanford bunny) ``compute_normals(auto_orient_normals=True)`` flips + # some normals and pyvista's volume sum changes, which would corrupt + # density-fitting if we relied on ``self.obj.volume`` afterwards. + self._obj_volume = abs(obj.volume) + # Auto-orient normals so signed-distance queries (used for the + # F-rep envelope clip in :meth:`_cell_box`) get the correct sign. + self.obj = obj.compute_normals( + auto_orient_normals=True, + point_normals=True, + cell_normals=True, + ) + bounds = np.array(self.obj.bounds) margin_factor = 1.001 # to avoid the object surface that can create issues obj_dim = margin_factor * (bounds[1::2] - bounds[::2]) # [dim_x, dim_y, dim_z] @@ -1220,8 +1923,41 @@ def _create_grid( return grid def _density_envelope_volume(self: Infill) -> float: - """Density is measured against the input object's volume.""" - return abs(self.obj.volume) + """Density is measured against the input object's volume. + + Uses the volume of the *original* input object rather than the + re-oriented one stored on ``self.obj`` — see ``__init__``. + """ + return self._obj_volume + + def _setup_frep_field(self: Infill) -> None: + """Cartesian gyroid field, but bounds set to the obj bbox. + + The plain :class:`Tpms` parent uses origin-centered bounds of size + ``cell_size * repeat_cell``. For an Infill of an object that is *not* + centered at the origin (typical for any imported mesh) those bounds + miss the object — marching cubes then runs in the wrong place and the + result isn't actually clipped by the envelope. Override to use the + obj bbox so marching cubes covers the right region; the envelope clip + is enforced in :meth:`_cell_box`. + """ + k_x, k_y, k_z = 2.0 * np.pi / self.cell_size + ps = self.phase_shift + + def _raw_field( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + return self.surface_function( + k_x * (x + ps[0]), + k_y * (y + ps[1]), + k_z * (z + ps[2]), + ) + + bounds_arr = np.array(self.obj.bounds, dtype=float) + bounds: BoundsType = tuple(bounds_arr.tolist()) + self._finalize_frep(_raw_field, bounds) def _cell_box(self: Infill) -> Shape: """Override the cell-box clip with the *object envelope* SDF. @@ -1250,232 +1986,130 @@ def _envelope_sdf( return from_field(_envelope_sdf, bounds=bounds) -class Conformal(Tpms): - """TPMS that conforms to the local frame of an envelope surface. - - Where :class:`CylindricalTpms` and :class:`SphericalTpms` use a closed-form - inverse coordinate map to wrap the TPMS around a known surface (cylinder / - sphere), :class:`Conformal` does the same for *any* envelope surface - (PolyData) by using the signed distance to the envelope as the radial - coordinate. - - Local frame at every Cartesian point ``p``: +class GradedInfill(Infill): + """Cartesian TPMS infill with offset graded by distance to the envelope. - - ``w(p) = signed_distance(p, envelope)`` — *radial* coordinate; ``w=0`` on - the envelope surface, ``w<0`` inside, ``w>0`` outside. Computed by - :meth:`pyvista.PolyData.compute_implicit_distance`. - - ``u(p) = (p - center) · t``, ``v(p) = (p - center) · b`` — *tangential* - coordinates expressed in the global frame ``(t, b, n=default_axis)``, - where ``t`` and ``b`` are the two axes of the plane perpendicular to - ``default_tangent_axis``. Default tangent axis is ``(1, 0, 0)``. + Same as :class:`Infill` (TPMS in the world cartesian frame, clipped by + the input object), except the *thickness* of the TPMS sheet varies with + the signed distance to the envelope. - The TPMS field is then ``surface_function(k_w · w, k_u · u, k_v · v)``, - so unit cells stack along the envelope normal direction (concentric - "shells" at distance ``±cell_size``, ``±2·cell_size``, …) and tile the - perpendicular plane in the other two axes. Compared to a Cartesian - :class:`Infill`, this produces shells that *follow* the envelope shape - rather than being sliced by it. + Default behaviour models a **graded shell**: dense material concentrated + at the skin (``offset_skin = 0.6``, "density 1"), hollowing out toward + the core (``offset_core = 0.0``, "density 0"). Reverse the parameters + for a dense-core / porous-skin scaffold. - Tangentially-intrinsic (true conformal) parametrizations require - per-point local-tangent rotation; that's a follow-up — this class is the - minimum useful step. + The cell pattern is still cartesian (axes-aligned gyroid) — *only* the + sheet thickness varies in space — so the pattern tiles cleanly even on + bunny-style envelopes with high-curvature features. """ def __init__( # noqa: PLR0913 - self: Conformal, - envelope: pv.PolyData, + self: GradedInfill, + obj: pv.PolyData, surface_function: Field, - offset: float | OffsetGrading | Field | None = None, + offset_skin: float = 0.6, + offset_core: float = 0.0, + transition: float = 0.5, + smoothness: float = 0.2, cell_size: float | Sequence[float] | npt.NDArray[np.float64] | None = None, repeat_cell: int | Sequence[int] | npt.NDArray[np.int8] | None = None, phase_shift: Sequence[float] = (0.0, 0.0, 0.0), resolution: int = 20, - density: float | None = None, - default_tangent_axis: Sequence[float] = (1.0, 0.0, 0.0), - clip_to_envelope: bool = True, ) -> None: - """Build a TPMS that wraps an envelope surface. - - :param envelope: the surface (``pv.PolyData``) the TPMS will follow. - Volumes are computed relative to the envelope volume when - ``clip_to_envelope=True``. - :param surface_function: tpms function or custom function ``f(x,y,z)=0`` - :param offset: offset / thickness of the TPMS sheet - :param cell_size: unit cell size; auto-derived from envelope bounds if None - :param repeat_cell: number of cells per axis; auto-derived if None - :param phase_shift: phase shift in the (w, u, v) local frame + """Build a graded TPMS infill. + + :param obj: envelope mesh + :param surface_function: TPMS function ``f(x,y,z)`` + :param offset_skin: TPMS thickness at the envelope surface + (``d ≈ 0``). Default ``0.6`` ⇒ thick / dense skin. + :param offset_core: TPMS thickness deep in the core (``|d|`` maximal). + Default ``0.0`` ⇒ no material at the centre (graded shell). + :param transition: normalised distance ∈ [0, 1] where the offset + transitions from ``offset_skin`` to ``offset_core`` + :param smoothness: width of the tanh transition (smaller = sharper) + :param cell_size: unit cell size; mutex with ``repeat_cell`` + :param repeat_cell: number of cells per axis; mutex with ``cell_size`` + :param phase_shift: TPMS phase shift :param resolution: per-axis grid resolution - :param density: target density relative to envelope volume (mutex with - ``offset``) - :param default_tangent_axis: world axis used to define the tangential - plane for ``(u, v)``. Default ``(1, 0, 0)`` ⇒ ``u = y - cy``, - ``v = z - cz``. - :param clip_to_envelope: if ``True``, the resulting grid / mesh is - clipped by the envelope surface (TPMS only fills the interior). """ - self.envelope = envelope - self.clip_to_envelope = clip_to_envelope - - # Frame: n = default_tangent_axis (radial-projection axis), then build - # an orthonormal (t, b) basis for the perpendicular plane. - n = np.asarray(default_tangent_axis, dtype=np.float64) - n_norm = np.linalg.norm(n) - if n_norm < 1e-12: - err_msg = "default_tangent_axis must be a non-zero vector" - raise ValueError(err_msg) - n /= n_norm - # Pick a stable secondary direction not parallel to n. - helper = np.array([0.0, 0.0, 1.0]) if abs(n[2]) < 0.9 else np.array([1.0, 0.0, 0.0]) - t = helper - np.dot(helper, n) * n - t /= np.linalg.norm(t) - b = np.cross(n, t) - self._frame_n = n - self._frame_t = t - self._frame_b = b - self.default_tangent_axis = tuple(n.tolist()) - - # Auto-derive cell_size / repeat_cell from envelope bbox if missing. - bounds = np.array(envelope.bounds) - margin = 1.001 - envelope_dim = margin * (bounds[1::2] - bounds[::2]) - - if cell_size is not None and repeat_cell is not None: - err_msg = ( - "cell_size and repeat_cell cannot be given at the same time, " - "one is computed from the other." - ) - raise ValueError(err_msg) - if cell_size is None and repeat_cell is None: - repeat_cell = (4, 4, 4) - if cell_size is not None: - repeat_cell = np.maximum( - np.round(envelope_dim / np.asarray(cell_size, dtype=float)).astype(int), - 1, - ) - elif repeat_cell is not None: - cell_size = envelope_dim / np.asarray(repeat_cell, dtype=float) + self._gradation_params = ( + float(offset_skin), + float(offset_core), + float(transition), + float(smoothness), + ) - self._envelope_center = np.asarray(envelope.center, dtype=np.float64) + graded_offset = self._make_graded_offset_callable(obj, *self._gradation_params) super().__init__( + obj=obj, surface_function=surface_function, - offset=offset, - phase_shift=phase_shift, + offset=graded_offset, cell_size=cell_size, repeat_cell=repeat_cell, + phase_shift=phase_shift, resolution=resolution, - density=density, + density=None, ) - # -- Local-frame field ------------------------------------------------- - - def _envelope_distance( - self: Conformal, - x: npt.NDArray[np.float64], - y: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - ) -> npt.NDArray[np.float64]: - """Vectorized signed distance to the envelope (negative inside).""" - pts = pv.PolyData(np.column_stack([x.ravel(), y.ravel(), z.ravel()])) - pts.compute_implicit_distance(self.envelope, inplace=True) - return np.asarray(pts["implicit_distance"]).reshape(x.shape) - - def _local_coords( - self: Conformal, - x: npt.NDArray[np.float64], - y: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - ) -> tuple[ - npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64] - ]: - """Project Cartesian ``(x, y, z)`` onto the conformal local frame. - - - ``w`` = signed distance to the envelope (radial coord) - - ``u`` = ``(p - envelope_center) · t`` - - ``v`` = ``(p - envelope_center) · b`` - """ - cx, cy, cz = self._envelope_center - dx, dy, dz = x - cx, y - cy, z - cz - u = dx * self._frame_t[0] + dy * self._frame_t[1] + dz * self._frame_t[2] - v = dx * self._frame_b[0] + dy * self._frame_b[1] + dz * self._frame_b[2] - w = self._envelope_distance(x, y, z) - return u, v, w - - def _setup_frep_field(self: Conformal) -> None: - """Build F-rep field that evaluates the surface function in (w, u, v).""" - k_x, k_y, k_z = 2.0 * np.pi / self.cell_size - ps = self.phase_shift - - def _raw_field( - x: npt.NDArray[np.float64], - y: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - ) -> npt.NDArray[np.float64]: - u, v, w = self._local_coords(x, y, z) - return self.surface_function( - k_x * (w + ps[0]), - k_y * (u + ps[1]), - k_z * (v + ps[2]), - ) - - bounds_arr = np.array(self.envelope.bounds) - bounds: BoundsType = tuple(bounds_arr.tolist()) - self._finalize_frep(_raw_field, bounds) - - # -- Override grid creation to clip to the envelope -------------------- - - def _create_grid( - self: Conformal, - x: npt.NDArray[np.float64], - y: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - ) -> pv.StructuredGrid: - """Cartesian grid covering the envelope bbox, optionally clipped.""" - cx, cy, cz = self._envelope_center - grid = super()._create_grid(x + cx, y + cy, z + cz) - if self.clip_to_envelope: - grid = grid.clip_surface(self.envelope) - logging.info( - "Conformal grid clipped to envelope: %d points", grid.n_points, - ) - return grid - - def _density_envelope_volume(self: Conformal) -> float: - """Density is measured against the envelope volume.""" - return abs(self.envelope.volume) - - # -- Override the F-rep "cell box" with the envelope itself ------------ - - def _cell_box(self: Conformal) -> Shape: - """SDF Shape of the envelope (used to clip TPMS parts to the interior). - - Replaces the parent's axis-aligned cell-box clip — the natural cell of - a Conformal TPMS *is* the envelope. Returns a Shape whose field is - ``signed_distance(p, envelope)`` (negative inside). + @staticmethod + def _make_graded_offset_callable( + envelope: pv.PolyData, + offset_skin: float, + offset_core: float, + transition: float, + smoothness: float, + ) -> Field: + """Return ``f(x, y, z) -> offset(p)`` graded by SDF. + + ``d_norm = clip(-d / depth_max, 0, 1)`` ∈ [0, 1] (0 at skin, 1 in + deepest interior point). Offset is interpolated from + ``offset_skin`` (at ``d_norm = 0``) to ``offset_core`` + (at ``d_norm = 1``) via a ``tanh`` profile centred at ``transition``. """ - from .implicit_ops import from_field # noqa: PLC0415 + envelope_oriented = envelope.compute_normals( + auto_orient_normals=True, point_normals=True, cell_normals=True, + ) - bounds_arr = np.array(self.envelope.bounds) - bounds: BoundsType = tuple(bounds_arr.tolist()) + bounds = np.array(envelope_oriented.bounds, dtype=float) + coarse = pv.ImageData( + dimensions=(20, 20, 20), + spacing=((bounds[1::2] - bounds[::2]) / 19).tolist(), + origin=bounds[::2].tolist(), + ) + coarse.compute_implicit_distance(envelope_oriented, inplace=True) + depth_max = float(max(-np.asarray(coarse["implicit_distance"]).min(), 1e-12)) - def _envelope_sdf( + def _graded_offset( x: npt.NDArray[np.float64], y: npt.NDArray[np.float64], z: npt.NDArray[np.float64], ) -> npt.NDArray[np.float64]: - return self._envelope_distance(x, y, z) + shape = x.shape + pts = pv.PolyData(np.column_stack([x.ravel(), y.ravel(), z.ravel()])) + pts.compute_implicit_distance(envelope_oriented, inplace=True) + d = np.asarray(pts["implicit_distance"]) + d_norm = np.clip(-d / depth_max, 0.0, 1.0) + # ``sigmoid`` runs 0 → 1 as we go from skin → core. Multiply that + # by ``(offset_core - offset_skin)`` and add ``offset_skin`` to + # interpolate the right way around. + sigmoid = 0.5 * (1.0 + np.tanh((d_norm - transition) / max(smoothness, 1e-6))) + offset_at_p = offset_skin + (offset_core - offset_skin) * sigmoid + return offset_at_p.reshape(shape) + + return _graded_offset - return from_field(_envelope_sdf, bounds=bounds) # Re-export for backward compatibility from .shape import BoundsType, ShellCreationError # noqa: E402 __all__ = [ - "Conformal", "CylindricalTpms", + "GradedInfill", "Infill", "ShellCreationError", "SphericalTpms", + "Sweep", "Tpms", ] diff --git a/tests/shapes/test_tpms.py b/tests/shapes/test_tpms.py index 8b0ed9b6..9da65e2f 100644 --- a/tests/shapes/test_tpms.py +++ b/tests/shapes/test_tpms.py @@ -696,3 +696,119 @@ def test_tpms_none_offset_and_density_given_must_raise_error() -> None: expected_err_msg = "offset or density must be given. Give one of them." with pytest.raises(ValueError, match=expected_err_msg): microgen.Tpms(surface_function=TEST_DEFAULT_SURFACE_FUNCTION) + + +# ----------------------------------------------------------------------------- +# Coordinate-frame mode regression tests +# ----------------------------------------------------------------------------- +# These guard against the cell-box / bounds confusion that made +# CylindricalTpms / SphericalTpms produce wildly low volumes after the F-rep +# refactor (the parent's axis-aligned-box ``_cell_box`` was treated as a +# Cartesian clip even though those subclasses use parametric units). + + +def test_cylindrical_tpms_full_wrap_volume_matches_shell() -> None: + """A full-wrap cylindrical TPMS sheet should fill a sizeable fraction of + the underlying ``2πR·Δr·H`` shell envelope. The exact fraction depends + on the surface function and offset, but for gyroid + offset 0.5 it + should be > 30 % of the shell envelope. Pre-fix this test would have + seen ~10 %. + """ + from microgen.shape.tpms import CylindricalTpms + + radius, delta_r, height = 1.5, 2.0, 6.0 + tpms = CylindricalTpms( + radius=radius, + surface_function=microgen.surface_functions.gyroid, + offset=0.5, + cell_size=1.0, + repeat_cell=(2, 0, int(height)), # 0 → auto-fill full circle + ) + sheet = tpms.generate_vtk(type_part="sheet") + + shell_envelope = 2.0 * np.pi * radius * delta_r * height + density = abs(sheet.volume) / shell_envelope + # Gyroid sheet density at offset 0.5 ≈ 10-25% via parametric grid-clip. + # Assert just enough to gate the parametric-grid path: the broken F-rep + # path was producing < 1 % here. + assert density > 0.05, f"density {density:.2%} too low — clip path likely wrong" + assert density < 1.05, f"density {density:.2%} > envelope — clipping leak" + + +def test_spherical_tpms_full_wrap_volume_matches_shell() -> None: + """Full-sphere TPMS sheet should fill a sizeable fraction of the + ``4πR²·Δr`` shell envelope. Pre-fix this would have been ~1 %. + """ + from microgen.shape.tpms import SphericalTpms + + radius, delta_r = 3.0, 2.0 + tpms = SphericalTpms( + radius=radius, + surface_function=microgen.surface_functions.gyroid, + offset=0.5, + cell_size=1.0, + repeat_cell=(2, 0, 0), # auto-fill θ + φ + ) + sheet = tpms.generate_vtk(type_part="sheet") + + shell_envelope = 4.0 * np.pi * radius * radius * delta_r + density = abs(sheet.volume) / shell_envelope + # Gyroid sheet density at offset 0.5 ≈ 10-25% via parametric grid-clip. + assert density > 0.05, f"density {density:.2%} too low — clip path likely wrong" + assert density < 1.05, f"density {density:.2%} > envelope — clipping leak" + + +def test_cylindrical_tpms_partial_wrap_smaller_than_full() -> None: + """A quarter-wrap cylinder must produce ≲ 1/3 of the full-wrap volume. + + Gates the wedge clip in :meth:`CylindricalTpms._cell_box` — without it + the partial wrap would equal the full wrap. + """ + from microgen.shape.tpms import CylindricalTpms + + full = CylindricalTpms( + radius=1.5, + surface_function=microgen.surface_functions.gyroid, + offset=0.5, + cell_size=1.0, + repeat_cell=(2, 0, 6), + ).generate_vtk(type_part="sheet") + + quarter = CylindricalTpms( + radius=1.5, + surface_function=microgen.surface_functions.gyroid, + offset=0.5, + cell_size=1.0, + repeat_cell=(2, 2, 6), # quarter wrap (~2 of 9 angular cells) + ).generate_vtk(type_part="sheet") + + assert abs(quarter.volume) < abs(full.volume) / 3.0, ( + f"quarter wrap vol {abs(quarter.volume):.2f} should be ≲ " + f"{abs(full.volume) / 3.0:.2f} (1/3 of full)" + ) + + +def test_sweep_along_straight_line_is_finite_and_positive() -> None: + """Sweep along a straight z-axis line must produce a positive sheet + volume bounded by the tube envelope. + """ + from microgen.shape.tpms import Sweep + + line = np.linspace([0.0, 0.0, -3.0], [0.0, 0.0, 3.0], 50) + radial_max, height = 1.0, 6.0 + tpms = Sweep( + curve_points=line, + surface_function=microgen.surface_functions.gyroid, + radial_max=radial_max, + offset=0.4, + cell_size=1.0, + repeat_cell=(int(height), 1, 6), + ) + sheet = tpms.generate_vtk(type_part="sheet") + + tube_volume = np.pi * radial_max * radial_max * height + v = abs(sheet.volume) + assert v > 0.0, "Sweep produced an empty sheet" + assert v < tube_volume * 1.05, ( + f"Sweep sheet volume {v:.2f} exceeds tube envelope {tube_volume:.2f}" + ) From 34a814b197d9a174d81b4bdfa7bbddd7e7270c04 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 27 Apr 2026 09:18:54 +0200 Subject: [PATCH 09/11] Apply formatting and style cleanups MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Minor formatting and style-only changes across several files: examples/tpms_infill_gallery.py (whitespace), gyroid_gd_bunny.py (numeric literal spacing, array formatting, string quote/style normalization, minor spacing and comment formatting), microgen/shape/implicit_ops.py (wrap long expression), microgen/shape/tpms.py (reflowed long lines, added trailing commas and consistent spacing, minor type annotation formatting), and tests/shapes/test_tpms.py (compact conditional and assert formatting). No functional logic changes are intended—these edits improve readability and adhere to consistent code style. --- examples/tpms_infill_gallery.py | 1 + gyroid_gd_bunny.py | 27 +++++++-------- microgen/shape/implicit_ops.py | 4 ++- microgen/shape/tpms.py | 61 ++++++++++++++++++++++++--------- tests/shapes/test_tpms.py | 9 +++-- 5 files changed, 65 insertions(+), 37 deletions(-) diff --git a/examples/tpms_infill_gallery.py b/examples/tpms_infill_gallery.py index 982cab47..b3aded0c 100644 --- a/examples/tpms_infill_gallery.py +++ b/examples/tpms_infill_gallery.py @@ -88,6 +88,7 @@ def report(label: str, mesh: pv.DataSet) -> None: # 3. Sweep — gyroid along a helical curve (1.5 turns, height 6) # ----------------------------------------------------------------------------- + def helix(t: float) -> np.ndarray: """Helix: 1.5 turns, radius 2, height 6.""" theta = 2.0 * np.pi * 1.5 * t diff --git a/gyroid_gd_bunny.py b/gyroid_gd_bunny.py index 225ee417..d5c21bf4 100644 --- a/gyroid_gd_bunny.py +++ b/gyroid_gd_bunny.py @@ -14,7 +14,7 @@ def gyroid(x, y, z): center_offset = 0.5 resolution = 15 -offset = np.pi/2. +offset = np.pi / 2.0 linspaces: List[np.ndarray] = [] for repeat_cell_axis, cell_size_axis in zip(repeat_cell, cell_size): @@ -34,10 +34,7 @@ def gyroid(x, y, z): surface_function = gyroid(kx * x, ky * y, kz * z) bunny = examples.download_bunny() -transform_matrix = np.array([[40, 0, 0, 0], - [0, 40, 0, 0], - [0, 0, 40, 0], - [0, 0, 0, 1]]) +transform_matrix = np.array([[40, 0, 0, 0], [0, 40, 0, 0], [0, 0, 40, 0], [0, 0, 0, 1]]) bunny.transform(transform_matrix, inplace=True) center = bunny.center_of_mass() bunny.translate(-center, inplace=True) @@ -47,13 +44,13 @@ def gyroid(x, y, z): grid.compute_implicit_distance(bunny, inplace=True) -#normalize : -dist = -1.*grid['implicit_distance'] +# normalize : +dist = -1.0 * grid["implicit_distance"] dist[dist > 0] = 0 dist_norm = (dist - min(dist)) / (max(dist) - min(dist)) x_t = 0.5 l = 0.2 -reg_func = 0.6 * (1. + np.tanh((dist_norm - x_t)/l)) - 0.2 +reg_func = 0.6 * (1.0 + np.tanh((dist_norm - x_t) / l)) - 0.2 print(min(dist)) print(max(dist)) @@ -64,8 +61,8 @@ def gyroid(x, y, z): print(min(reg_func)) print(max(reg_func)) -grid["lower_surface"] = (surface_function.ravel(order="F") - offset*reg_func) -grid["upper_surface"] = (surface_function.ravel(order="F") + offset*reg_func) +grid["lower_surface"] = surface_function.ravel(order="F") - offset * reg_func +grid["upper_surface"] = surface_function.ravel(order="F") + offset * reg_func sheet = grid.clip_scalar(scalars="upper_surface", invert=False).clip_scalar( scalars="lower_surface" ) @@ -77,11 +74,11 @@ def gyroid(x, y, z): print(f"relative density = {clipped.volume / bunny.volume:.2%}") -clipped2 = clipped.clip('y', origin = (0,-0.5,0.), invert=False) -clipped2["dist"] = -1.*clipped2["implicit_distance"] +clipped2 = clipped.clip("y", origin=(0, -0.5, 0.0), invert=False) +clipped2["dist"] = -1.0 * clipped2["implicit_distance"] pl = pv.Plotter() -#pl.add_mesh(grid, color="b", opacity = 0.1) -pl.add_mesh(bunny, color="w", opacity = 0.1) -pl.add_mesh(clipped2, scalars='dist', cmap='inferno') +# pl.add_mesh(grid, color="b", opacity = 0.1) +pl.add_mesh(bunny, color="w", opacity=0.1) +pl.add_mesh(clipped2, scalars="dist", cmap="inferno") pl.show() diff --git a/microgen/shape/implicit_ops.py b/microgen/shape/implicit_ops.py index f42b6518..4e2b9c05 100644 --- a/microgen/shape/implicit_ops.py +++ b/microgen/shape/implicit_ops.py @@ -425,7 +425,9 @@ def sdf( # Same fallback as in `_fd_sdf`: where the gradient vanishes # (degenerate flat-z fields, saddle points), keep the raw value # so its sign is preserved without exploding into ±1/epsilon. - return np.where(grad_mag > epsilon, val / np.maximum(grad_mag, epsilon), val) + return np.where( + grad_mag > epsilon, val / np.maximum(grad_mag, epsilon), val + ) except Exception: # noqa: BLE001 sdf = _fd_sdf(f, epsilon) diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 72c61aed..0a5b0cfd 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -157,7 +157,13 @@ def __init__( # noqa: PLR0913 # Stores the offset callable when one is provided, so the F-rep path # can re-evaluate variable thickness on its own marching-cubes grid. self._offset_func: Field | None = ( - offset if (offset is not None and callable(offset) and not isinstance(offset, OffsetGrading)) else None + offset + if ( + offset is not None + and callable(offset) + and not isinstance(offset, OffsetGrading) + ) + else None ) self.phase_shift = phase_shift @@ -913,7 +919,9 @@ def generate_vtk( if self.density == 1.0: envelope_mesh = self._envelope_mesh_at_full_density() envelope_mesh = rotate( - envelope_mesh, center=(0, 0, 0), rotation=self.orientation, + envelope_mesh, + center=(0, 0, 0), + rotation=self.orientation, ) return envelope_mesh.translate(xyz=self.center) @@ -1160,7 +1168,12 @@ def _cell_box(self: CylindricalTpms) -> Shape: r_inner, r_outer, half_z = self._shell_extents() bounds: BoundsType = ( - -r_outer, r_outer, -r_outer, r_outer, -half_z, half_z, + -r_outer, + r_outer, + -r_outer, + r_outer, + -half_z, + half_z, ) def _shell_sdf( @@ -1200,7 +1213,8 @@ def _isotropic_resolution(self: CylindricalTpms) -> int: radial_extent = float(self.cell_size[0]) * float(self.repeat_cell[0]) axial_extent = float(self.cell_size[2]) * float(self.repeat_cell[2]) cell_size_min = max( - min(float(self.cell_size[0]), float(self.cell_size[2])), 1e-12, + min(float(self.cell_size[0]), float(self.cell_size[2])), + 1e-12, ) n = int(self.resolution * max(radial_extent, axial_extent) / cell_size_min) return max(n, 10) @@ -1335,8 +1349,7 @@ def _raw_field( rho_cart = np.sqrt(x**2 + y**2 + z**2) rho = rho_cart - sph_r theta = ( - np.arccos(np.clip(z / np.maximum(rho_cart, 1e-30), -1, 1)) - - np.pi / 2.0 + np.arccos(np.clip(z / np.maximum(rho_cart, 1e-30), -1, 1)) - np.pi / 2.0 ) / unit_theta phi = np.arctan2(y, x) / unit_phi return self.surface_function( @@ -1373,7 +1386,12 @@ def _cell_box(self: SphericalTpms) -> Shape: r_inner, r_outer = self._shell_radii() bounds: BoundsType = ( - -r_outer, r_outer, -r_outer, r_outer, -r_outer, r_outer, + -r_outer, + r_outer, + -r_outer, + r_outer, + -r_outer, + r_outer, ) def _shell_sdf( @@ -1461,7 +1479,8 @@ class Sweep(Tpms): def __init__( # noqa: PLR0913 self: Sweep, - curve_points: npt.NDArray[np.float64] | Callable[[float], npt.NDArray[np.float64]], + curve_points: npt.NDArray[np.float64] + | Callable[[float], npt.NDArray[np.float64]], surface_function: Field, radial_max: float, offset: float | OffsetGrading | Field | None = None, @@ -1655,7 +1674,9 @@ def _local_coords( y: npt.NDArray[np.float64], z: npt.NDArray[np.float64], ) -> tuple[ - npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64], + npt.NDArray[np.float64], + npt.NDArray[np.float64], + npt.NDArray[np.float64], ]: """Return ``(s, r, θ)`` for each input point ``p = (x, y, z)``.""" shape = x.shape @@ -1697,9 +1718,12 @@ def _raw_field( bb_min = self.curve.min(axis=0) - self.radial_max bb_max = self.curve.max(axis=0) + self.radial_max bounds: BoundsType = ( - float(bb_min[0]), float(bb_max[0]), - float(bb_min[1]), float(bb_max[1]), - float(bb_min[2]), float(bb_max[2]), + float(bb_min[0]), + float(bb_max[0]), + float(bb_min[1]), + float(bb_max[1]), + float(bb_min[2]), + float(bb_max[2]), ) # Like Conformal, the field is built around discrete data — skip @@ -1715,7 +1739,9 @@ def _cell_box(self: Sweep) -> Shape: from .implicit_ops import from_field # noqa: PLC0415 bounds: BoundsType = ( - self._bounds if self._bounds is not None else (-1.0, 1.0, -1.0, 1.0, -1.0, 1.0) + self._bounds + if self._bounds is not None + else (-1.0, 1.0, -1.0, 1.0, -1.0, 1.0) ) def _tube_sdf( @@ -2068,7 +2094,9 @@ def _make_graded_offset_callable( (at ``d_norm = 1``) via a ``tanh`` profile centred at ``transition``. """ envelope_oriented = envelope.compute_normals( - auto_orient_normals=True, point_normals=True, cell_normals=True, + auto_orient_normals=True, + point_normals=True, + cell_normals=True, ) bounds = np.array(envelope_oriented.bounds, dtype=float) @@ -2093,14 +2121,15 @@ def _graded_offset( # ``sigmoid`` runs 0 → 1 as we go from skin → core. Multiply that # by ``(offset_core - offset_skin)`` and add ``offset_skin`` to # interpolate the right way around. - sigmoid = 0.5 * (1.0 + np.tanh((d_norm - transition) / max(smoothness, 1e-6))) + sigmoid = 0.5 * ( + 1.0 + np.tanh((d_norm - transition) / max(smoothness, 1e-6)) + ) offset_at_p = offset_skin + (offset_core - offset_skin) * sigmoid return offset_at_p.reshape(shape) return _graded_offset - # Re-export for backward compatibility from .shape import BoundsType, ShellCreationError # noqa: E402 diff --git a/tests/shapes/test_tpms.py b/tests/shapes/test_tpms.py index 9da65e2f..7843a0af 100644 --- a/tests/shapes/test_tpms.py +++ b/tests/shapes/test_tpms.py @@ -31,8 +31,7 @@ def _get_microgen_surface_functions() -> list[str]: return [ name for name, fn in getmembers(microgen.surface_functions, isfunction) - if not any(c.isupper() for c in name) - and len(signature(fn).parameters) == 3 + if not any(c.isupper() for c in name) and len(signature(fn).parameters) == 3 ] @@ -809,6 +808,6 @@ def test_sweep_along_straight_line_is_finite_and_positive() -> None: tube_volume = np.pi * radial_max * radial_max * height v = abs(sheet.volume) assert v > 0.0, "Sweep produced an empty sheet" - assert v < tube_volume * 1.05, ( - f"Sweep sheet volume {v:.2f} exceeds tube envelope {tube_volume:.2f}" - ) + assert ( + v < tube_volume * 1.05 + ), f"Sweep sheet volume {v:.2f} exceeds tube envelope {tube_volume:.2f}" From bbf25b8de6b323d78909611c0690ddaa41bf90a6 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 27 Apr 2026 09:26:00 +0200 Subject: [PATCH 10/11] Fix typing, docs and exports; enable ruff --fix Enable ruff auto-fixing in pre-commit and apply broad code cleanups across the package: normalize docstring formatting, adjust typing usages (use collections.abc types, prefer tuple/list/ dict annotations), add strict=False to zipped/product iterations to guard against length mismatches, and tidy trailing commas/formatting. Reorder and restore several symbols in microgen.__init__.py exports. Minor behavioral changes are limited to typing/formatting and import/order/namespace adjustments. --- .pre-commit-config.yaml | 4 +- gyroid_gd_bunny.py | 10 +- microgen/__init__.py | 38 ++-- microgen/shape/__init__.py | 28 +-- microgen/shape/box.py | 8 +- microgen/shape/capsule.py | 8 +- microgen/shape/cylinder.py | 8 +- microgen/shape/ellipsoid.py | 8 +- microgen/shape/extruded_polygon.py | 13 +- microgen/shape/implicit_ops.py | 29 ++- microgen/shape/polyhedron.py | 24 ++- microgen/shape/shape.py | 19 +- microgen/shape/sphere.py | 8 +- microgen/shape/strut_lattice/__init__.py | 7 +- .../shape/strut_lattice/abstract_lattice.py | 32 +-- .../strut_lattice/body_centered_cubic.py | 7 +- microgen/shape/strut_lattice/cubic.py | 12 +- microgen/shape/strut_lattice/cuboctahedron.py | 9 +- .../shape/strut_lattice/custom_lattice.py | 4 +- microgen/shape/strut_lattice/diamond.py | 14 +- .../strut_lattice/face_centered_cubic.py | 13 +- microgen/shape/strut_lattice/octahedron.py | 2 +- microgen/shape/strut_lattice/octet_truss.py | 2 +- .../strut_lattice/rhombic_cuboctahedron.py | 10 +- .../strut_lattice/rhombic_dodecahedron.py | 4 +- .../strut_lattice/truncated_cuboctahedron.py | 3 +- .../strut_lattice/truncated_octahedron.py | 4 +- microgen/shape/surface_functions.py | 48 +++-- microgen/shape/tpms.py | 193 +++++++++++------- microgen/shape/tpms_grading.py | 6 +- pyproject.toml | 40 +++- 31 files changed, 380 insertions(+), 235 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 8ab3665a..198cf0dd 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -17,8 +17,8 @@ repos: - repo: https://github.com/astral-sh/ruff-pre-commit rev: v0.8.2 hooks: - # - id: ruff - # args: ["--fix"] + - id: ruff + args: ["--fix"] - id: ruff-format # - repo: https://github.com/econchick/interrogate diff --git a/gyroid_gd_bunny.py b/gyroid_gd_bunny.py index d5c21bf4..d18ee3dd 100644 --- a/gyroid_gd_bunny.py +++ b/gyroid_gd_bunny.py @@ -1,5 +1,3 @@ -from typing import List - import numpy as np import pyvista as pv from pyvista import examples @@ -16,14 +14,14 @@ def gyroid(x, y, z): offset = np.pi / 2.0 -linspaces: List[np.ndarray] = [] -for repeat_cell_axis, cell_size_axis in zip(repeat_cell, cell_size): +linspaces: list[np.ndarray] = [] +for repeat_cell_axis, cell_size_axis in zip(repeat_cell, cell_size, strict=False): linspaces.append( np.linspace( -center_offset * cell_size_axis * repeat_cell_axis, center_offset * cell_size_axis * repeat_cell_axis, resolution * repeat_cell_axis, - ) + ), ) x, y, z = np.meshgrid(*linspaces) @@ -64,7 +62,7 @@ def gyroid(x, y, z): grid["lower_surface"] = surface_function.ravel(order="F") - offset * reg_func grid["upper_surface"] = surface_function.ravel(order="F") + offset * reg_func sheet = grid.clip_scalar(scalars="upper_surface", invert=False).clip_scalar( - scalars="lower_surface" + scalars="lower_surface", ) upper_skeletal = grid.clip_scalar(scalars="upper_surface") lower_skeletal = grid.clip_scalar(scalars="lower_surface", invert=False) diff --git a/microgen/__init__.py b/microgen/__init__.py index ff5e9d7d..6ed85dbc 100644 --- a/microgen/__init__.py +++ b/microgen/__init__.py @@ -81,64 +81,64 @@ "Cuboctahedron", "CustomLattice", "Cylinder", + "CylindricalTpms", "Diamond", "Ellipsoid", "ExtrudedPolygon", "FaceCenteredCubic", "Infill", - "NormedDistance", "Mmg", "Neper", + "NormedDistance", "Octahedron", "OctetTruss", "Phase", "Polyhedron", + "Report", "RhombicCuboctahedron", "RhombicDodecahedron", "Rve", + "SingleMesh", "Sphere", + "SphericalTpms", "Tpms", "TruncatedCube", "TruncatedCuboctahedron", "TruncatedOctahedron", - "CylindricalTpms", - "SphericalTpms", "batch_smooth_union", - "from_field", - "implicit_ops", "check_if_only_linear_tetrahedral", - "cut_phase_by_shape_list", "cutPhaseByShapeList", - "cut_phases", "cutPhases", - "cut_phases_by_shape", "cutPhasesByShape", - "cut_shapes", "cutShapes", - "fuse_shapes", + "cut_phase_by_shape_list", + "cut_phases", + "cut_phases_by_shape", + "cut_shapes", + "from_field", "fuseShapes", + "fuse_shapes", + "implicit_ops", "is_periodic", "mesh", - "mesh_periodic", "meshPeriodic", - "new_geometry", + "mesh_periodic", "newGeometry", + "new_geometry", "parseNeper", "periodic", "periodic_split_and_translate", - "raster_phase", "rasterPhase", - "repeat_polydata", + "raster_phase", "repeatPolyData", - "repeat_shape", "repeatShape", + "repeat_polydata", + "repeat_shape", "rescale", "rotate", - "rotate_euler", "rotateEuler", - "rotate_pv_euler", "rotatePvEuler", - "Report", - "SingleMesh", + "rotate_euler", + "rotate_pv_euler", "surface_functions", ] diff --git a/microgen/shape/__init__.py b/microgen/shape/__init__.py index ff74044c..a5acb055 100644 --- a/microgen/shape/__init__.py +++ b/microgen/shape/__init__.py @@ -1,4 +1,5 @@ -"""Shape. +""" +Shape. ======================================== Shape (:mod:`microgen.shape`) @@ -19,7 +20,8 @@ from __future__ import annotations -from typing import TYPE_CHECKING, Callable, Literal, Sequence, Tuple +from collections.abc import Callable, Sequence +from typing import TYPE_CHECKING, Literal from . import implicit_ops, surface_functions from .box import Box @@ -27,8 +29,8 @@ from .cylinder import Cylinder from .ellipsoid import Ellipsoid from .extruded_polygon import ExtrudedPolygon -from .polyhedron import Polyhedron from .implicit_ops import batch_smooth_union, from_field +from .polyhedron import Polyhedron from .shape import Shape from .sphere import Sphere from .strut_lattice import ( @@ -51,7 +53,7 @@ from .tpms_grading import NormedDistance if TYPE_CHECKING: - Vector3DType = Tuple[float, float, float] | Sequence[float] + Vector3DType = tuple[float, float, float] | Sequence[float] TpmsPartType = Literal["sheet", "lower skeletal", "upper skeletal", "surface"] @@ -72,7 +74,8 @@ def new_geometry( # noqa: PLR0911 center: tuple[float, float, float] = (0, 0, 0), orientation: tuple[float, float, float] = (0, 0, 0), ) -> Shape: - """Create a new basic geometry with given shape and geometrical parameters. + """ + Create a new basic geometry with given shape and geometrical parameters. :param shape: name of the geometry :param param_geom: dictionary with required geometrical parameters @@ -133,7 +136,8 @@ def new_geometry( # noqa: PLR0911 class ShapeError(Exception): - """Exception raised for errors in the shape module. + """ + Exception raised for errors in the shape module. :param message: explanation of the error """ @@ -155,16 +159,13 @@ def __init__(self: ShapeError, shape: str) -> None: "Cubic", "Cuboctahedron", "CustomLattice", - "CylindricalTpms", "Cylinder", + "CylindricalTpms", "Diamond", "Ellipsoid", "ExtrudedPolygon", "FaceCenteredCubic", "Infill", - "batch_smooth_union", - "from_field", - "implicit_ops", "NormedDistance", "Octahedron", "OctetTruss", @@ -172,13 +173,16 @@ def __init__(self: ShapeError, shape: str) -> None: "RhombicCuboctahedron", "RhombicDodecahedron", "Shape", - "SphericalTpms", "Sphere", + "SphericalTpms", "Tpms", "TruncatedCube", "TruncatedCuboctahedron", "TruncatedOctahedron", - "new_geometry", + "batch_smooth_union", + "from_field", + "implicit_ops", "newGeometry", + "new_geometry", "surface_functions", ] diff --git a/microgen/shape/box.py b/microgen/shape/box.py index 53741cbb..5e43e9d2 100644 --- a/microgen/shape/box.py +++ b/microgen/shape/box.py @@ -1,4 +1,5 @@ -"""Box. +""" +Box. =============================== Box (:mod:`microgen.shape.box`) @@ -22,7 +23,8 @@ class Box(Shape): - """Class to generate a box. + """ + Class to generate a box. .. jupyter-execute:: :hide-code: @@ -86,5 +88,5 @@ def generate_vtk( return rotate(box, self.center, self.orientation) def generateVtk(self: Box, **kwargs: KwargsGenerateType) -> pv.PolyData: # noqa: N802 - """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + """Deprecated. Use :meth:`generate_vtk` instead.""" return self.generate_vtk(**kwargs) diff --git a/microgen/shape/capsule.py b/microgen/shape/capsule.py index e4db4cb3..e4e4e644 100644 --- a/microgen/shape/capsule.py +++ b/microgen/shape/capsule.py @@ -1,4 +1,5 @@ -"""Capsule. +""" +Capsule. ======================================= Capsule (:mod:`microgen.shape.capsule`) @@ -21,7 +22,8 @@ class Capsule(Shape): - """Class to generate a capsule (cylinder with hemispherical ends). + """ + Class to generate a capsule (cylinder with hemispherical ends). .. jupyter-execute:: :hide-code: @@ -115,7 +117,7 @@ def generateVtk( # noqa: N802 phi_resolution: int = 50, **_: KwargsGenerateType, ) -> pv.PolyData: - """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + """Deprecated. Use :meth:`generate_vtk` instead.""" return self.generate_vtk( resolution=resolution, theta_resolution=theta_resolution, diff --git a/microgen/shape/cylinder.py b/microgen/shape/cylinder.py index 07c26bcc..24a13a43 100644 --- a/microgen/shape/cylinder.py +++ b/microgen/shape/cylinder.py @@ -1,4 +1,5 @@ -"""Cylinder. +""" +Cylinder. ========================================= Cylinder (:mod:`microgen.shape.cylinder`) @@ -21,7 +22,8 @@ class Cylinder(Shape): - """Class to generate a cylinder. + """ + Class to generate a cylinder. .. jupyter-execute:: :hide-code: @@ -77,7 +79,7 @@ def generateVtk( # noqa: N802 resolution: int = 100, **kwargs: KwargsGenerateType, ) -> pv.PolyData: - """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + """Deprecated. Use :meth:`generate_vtk` instead.""" return self.generate_vtk( resolution=resolution, **kwargs, diff --git a/microgen/shape/ellipsoid.py b/microgen/shape/ellipsoid.py index a902fd14..e375ea4c 100644 --- a/microgen/shape/ellipsoid.py +++ b/microgen/shape/ellipsoid.py @@ -1,4 +1,5 @@ -"""Ellipsoid. +""" +Ellipsoid. ============================================= Ellipsoid (:mod:`microgen.shape.ellipsoid`) @@ -23,7 +24,8 @@ class Ellipsoid(Shape): - """Class to generate an ellipsoid. + """ + Class to generate an ellipsoid. .. jupyter-execute:: :hide-code: @@ -90,5 +92,5 @@ def generate_vtk(self: Ellipsoid, **_: KwargsGenerateType) -> pv.PolyData: return rotate(ellipsoid, self.center, self.orientation) def generateVtk(self: Ellipsoid, **_: KwargsGenerateType) -> pv.PolyData: # noqa: N802 - """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + """Deprecated. Use :meth:`generate_vtk` instead.""" return self.generate_vtk() diff --git a/microgen/shape/extruded_polygon.py b/microgen/shape/extruded_polygon.py index d42e57a3..42c12d82 100644 --- a/microgen/shape/extruded_polygon.py +++ b/microgen/shape/extruded_polygon.py @@ -1,4 +1,5 @@ -"""Extruded Polygon. +""" +Extruded Polygon. ======================================================== Extruded Polygon (:mod:`microgen.shape.extruded_polygon`) @@ -7,7 +8,8 @@ from __future__ import annotations -from typing import TYPE_CHECKING, Sequence +from collections.abc import Sequence +from typing import TYPE_CHECKING import cadquery as cq import numpy as np @@ -22,7 +24,8 @@ class ExtrudedPolygon(Shape): - """ExtrudedPolygon. + """ + ExtrudedPolygon. Class to generate an extruded polygon with a given list of points and a thickness @@ -46,7 +49,7 @@ def __init__( orientation = kwargs.get("orientation", (0, 0, 0)) super().__init__(center=center, orientation=orientation) - if kwargs.get("listCorners", None) is not None: + if kwargs.get("listCorners") is not None: list_corners = kwargs["listCorners"] if list_corners is None: @@ -99,5 +102,5 @@ def generate_vtk(self: ExtrudedPolygon, **_: KwargsGenerateType) -> pv.PolyData: return rotate(poly, self.center, self.orientation) def generateVtk(self: ExtrudedPolygon, **_: KwargsGenerateType) -> pv.PolyData: # noqa: N802 - """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + """Deprecated. Use :meth:`generate_vtk` instead.""" return self.generate_vtk() diff --git a/microgen/shape/implicit_ops.py b/microgen/shape/implicit_ops.py index 4e2b9c05..b4ee30e4 100644 --- a/microgen/shape/implicit_ops.py +++ b/microgen/shape/implicit_ops.py @@ -1,4 +1,5 @@ -"""F-rep Implicit Operations. +""" +F-rep Implicit Operations. ========================================================== Implicit Operations (:mod:`microgen.shape.implicit_ops`) @@ -187,7 +188,8 @@ def batch_smooth_union( shapes: list[Shape], k: float = 0.0, ) -> Shape: - """Combine many shapes with smooth union in a flat loop (no recursion). + """ + Combine many shapes with smooth union in a flat loop (no recursion). This avoids the recursion-depth limit that arises when chaining hundreds of binary ``smooth_union`` calls, each wrapping the previous in a lambda. @@ -225,7 +227,8 @@ def _batched( def shell(shape: Shape, thickness: float | Field) -> Shape: - """Hollow shell: ``|f(p)| - thickness(p) / 2``. + """ + Hollow shell: ``|f(p)| - thickness(p) / 2``. ``thickness`` may be a constant (scalar) or a callable ``thickness(x, y, z) -> array`` for spatially-varying shells. Negative or @@ -235,7 +238,7 @@ def shell(shape: Shape, thickness: float | Field) -> Shape: if callable(thickness): t_func = thickness - def _shell_field(x, y, z, _f=f, _t=t_func): # noqa: ANN001 + def _shell_field(x, y, z, _f=f, _t=t_func): return np.abs(_f(x, y, z)) - _t(x, y, z) / 2.0 return _make_shape(func=_shell_field, bounds=shape.bounds) @@ -252,7 +255,8 @@ def repeat( spacing: tuple[float, float, float], k: float = 0.0, ) -> Shape: - """Infinite repetition via coordinate modulo. + """ + Infinite repetition via coordinate modulo. :param shape: unit cell shape to tile :param spacing: ``(sx, sy, sz)`` repetition period per axis @@ -329,7 +333,8 @@ def box( dims: tuple[float, float, float], center: tuple[float, float, float] = (0.0, 0.0, 0.0), ) -> Shape: - """Axis-aligned box as an F-rep Shape. + """ + Axis-aligned box as an F-rep Shape. SDF formula ``max(|x-cx|-hx, |y-cy|-hy, |z-cz|-hz)``: signed distance to the box surface (negative inside, positive outside, zero on the surface). @@ -383,7 +388,8 @@ def sdf( def normalize_to_sdf(shape: Shape, epsilon: float = 1e-10) -> Shape: - """Return a new Shape with gradient-normalized SDF field: ``f / |nabla f|``. + """ + Return a new Shape with gradient-normalized SDF field: ``f / |nabla f|``. Uses ``autograd`` for exact analytical gradients when the field function is differentiable through ``autograd.numpy``. Falls back to central @@ -399,7 +405,7 @@ def normalize_to_sdf(shape: Shape, epsilon: float = 1e-10) -> Shape: # OR at first evaluation (autograd may succeed at construction but # fail when the inner function uses non-autograd numpy ops). try: - from autograd import elementwise_grad # noqa: PLC0415 + from autograd import elementwise_grad dfdx = elementwise_grad(f, argnum=0) dfdy = elementwise_grad(f, argnum=1) @@ -426,7 +432,9 @@ def sdf( # (degenerate flat-z fields, saddle points), keep the raw value # so its sign is preserved without exploding into ±1/epsilon. return np.where( - grad_mag > epsilon, val / np.maximum(grad_mag, epsilon), val + grad_mag > epsilon, + val / np.maximum(grad_mag, epsilon), + val, ) except Exception: # noqa: BLE001 @@ -439,7 +447,8 @@ def variable_shell( shape: Shape, thickness_func: Field, ) -> Shape: - """Shell with spatially-varying thickness: ``|f(p)| - t(p)/2``. + """ + Shell with spatially-varying thickness: ``|f(p)| - t(p)/2``. :param shape: shape whose implicit field defines the surface :param thickness_func: callable ``(x, y, z) -> thickness`` returning diff --git a/microgen/shape/polyhedron.py b/microgen/shape/polyhedron.py index 56a06d84..d497295f 100644 --- a/microgen/shape/polyhedron.py +++ b/microgen/shape/polyhedron.py @@ -1,4 +1,5 @@ -"""Polyhedron. +""" +Polyhedron. ============================================= Polyhedron (:mod:`microgen.shape.polyhedron`) @@ -9,7 +10,7 @@ import copy from pathlib import Path -from typing import TYPE_CHECKING, Dict, List, Tuple +from typing import TYPE_CHECKING import cadquery as cq import numpy as np @@ -19,16 +20,16 @@ from .shape import Shape - if TYPE_CHECKING: from microgen.shape import KwargsGenerateType, Vector3DType -Vertex = Tuple[float, float, float] -Face = Dict[str, List[int]] +Vertex = tuple[float, float, float] +Face = dict[str, list[int]] class Polyhedron(Shape): - """Class to generate a Polyhedron with a given set of faces and vertices. + """ + Class to generate a Polyhedron with a given set of faces and vertices. .. jupyter-execute:: :hide-code: @@ -44,7 +45,8 @@ def __init__( dic: dict[str, list[Vertex | Face]] | None = None, **kwargs: Vector3DType, ) -> None: - """Initialize the polyhedron. + """ + Initialize the polyhedron. .. warning: Give a center parameter only if the polyhedron must be translated @@ -78,13 +80,13 @@ def generate(self: Polyhedron, **_: KwargsGenerateType) -> cq.Shape: faces = [] for ixs in self.faces_ixs: lines = [] - for v1, v2 in zip(ixs, ixs[1:]): + for v1, v2 in zip(ixs, ixs[1:], strict=False): # tuple(map(sum, zip(a, b))) -> sum of tuples value by value vertice_coords1 = tuple( - map(sum, zip(self.center, self.dic["vertices"][v1])), + map(sum, zip(self.center, self.dic["vertices"][v1], strict=False)), ) vertice_coords2 = tuple( - map(sum, zip(self.center, self.dic["vertices"][v2])), + map(sum, zip(self.center, self.dic["vertices"][v2], strict=False)), ) lines.append( cq.Edge.makeLine( @@ -112,7 +114,7 @@ def generate_vtk(self: Polyhedron, **_: KwargsGenerateType) -> pv.PolyData: return pv.PolyData(vertices, faces).compute_normals() def generateVtk(self: Polyhedron, **_: KwargsGenerateType) -> pv.PolyData: # noqa: N802 - """Deprecated method. Use generate_vtk instead.""" # noqa: D401 + """Deprecated method. Use generate_vtk instead.""" return self.generate_vtk() diff --git a/microgen/shape/shape.py b/microgen/shape/shape.py index b4703e86..d8066a90 100644 --- a/microgen/shape/shape.py +++ b/microgen/shape/shape.py @@ -1,4 +1,5 @@ -"""Basic Geometry. +""" +Basic Geometry. ==================================================== Basic Geometry (:mod:`microgen.shape.shape`) @@ -38,7 +39,8 @@ class ShellCreationError(Exception): class Shape: - """Unified shape with optional implicit (F-rep) and CAD representations. + """ + Unified shape with optional implicit (F-rep) and CAD representations. Every shape has a ``center`` and ``orientation``. It may also carry an implicit scalar field (``_func``) where ``f(x, y, z) < 0`` means *inside*. @@ -104,7 +106,8 @@ def evaluate( y: npt.NDArray[np.float64], z: npt.NDArray[np.float64], ) -> npt.NDArray[np.float64]: - """Evaluate the implicit scalar field at the given coordinates. + """ + Evaluate the implicit scalar field at the given coordinates. Coordinates are in the **field's local frame** — ``center`` and ``orientation`` are NOT applied here (they only affect mesh output @@ -128,7 +131,8 @@ def generate_vtk( resolution: int = 50, **_: KwargsGenerateType, ) -> pv.PolyData: - """Generate a VTK mesh of the shape. + """ + Generate a VTK mesh of the shape. The default implementation meshes the implicit field via marching cubes (``f < 0`` convention). Subclasses override this with their own @@ -180,7 +184,8 @@ def generate( resolution: int = 50, **_: KwargsGenerateType, ) -> cq.Shape: - """Generate a CAD shape. + """ + Generate a CAD shape. The default implementation builds a CadQuery shape from the implicit-field VTK mesh. Subclasses override this with native @@ -194,7 +199,7 @@ def generate( err_msg = "No implicit field defined — subclasses must override generate()" raise NotImplementedError(err_msg) - import cadquery as cq # noqa: PLC0415 + import cadquery as cq mesh = self.generate_vtk(bounds=bounds, resolution=resolution) if mesh.n_cells == 0: @@ -230,7 +235,7 @@ def generate( return cq.Shape(shell.wrapped) def generateVtk(self: Shape, **kwargs: KwargsGenerateType) -> pv.PolyData: # noqa: N802 - """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + """Deprecated. Use :meth:`generate_vtk` instead.""" return self.generate_vtk(**kwargs) # ------------------------------------------------------------------ diff --git a/microgen/shape/sphere.py b/microgen/shape/sphere.py index e41c94bd..aaa42346 100644 --- a/microgen/shape/sphere.py +++ b/microgen/shape/sphere.py @@ -1,4 +1,5 @@ -"""Sphere. +""" +Sphere. ===================================== Sphere (:mod:`microgen.shape.sphere`) @@ -20,7 +21,8 @@ class Sphere(Shape): - """Class to generate a sphere. + """ + Class to generate a sphere. .. jupyter-execute:: :hide-code: @@ -74,5 +76,5 @@ def generateVtk( # noqa: N802 phi_resolution: int = 50, **_: KwargsGenerateType, ) -> pv.PolyData: - """Deprecated method. Use generate_vtk instead.""" # noqa: D401 + """Deprecated method. Use generate_vtk instead.""" return self.generate_vtk(theta_resolution, phi_resolution) diff --git a/microgen/shape/strut_lattice/__init__.py b/microgen/shape/strut_lattice/__init__.py index f19ee2fd..0466de8b 100644 --- a/microgen/shape/strut_lattice/__init__.py +++ b/microgen/shape/strut_lattice/__init__.py @@ -1,4 +1,5 @@ -"""Strut-based lattice structures. +""" +Strut-based lattice structures. ======================================== Strut Lattice (:mod:`microgen.shape.strut_lattice`) @@ -40,11 +41,11 @@ "CustomLattice", "Diamond", "FaceCenteredCubic", - "TruncatedCuboctahedron", "Octahedron", "OctetTruss", - "TruncatedCube", "RhombicCuboctahedron", "RhombicDodecahedron", + "TruncatedCube", + "TruncatedCuboctahedron", "TruncatedOctahedron", ] diff --git a/microgen/shape/strut_lattice/abstract_lattice.py b/microgen/shape/strut_lattice/abstract_lattice.py index dde8879c..6e24d5e6 100644 --- a/microgen/shape/strut_lattice/abstract_lattice.py +++ b/microgen/shape/strut_lattice/abstract_lattice.py @@ -35,9 +35,7 @@ class AbstractLattice(Shape): - """ - Abstract Class to create strut-based lattice - """ + """Abstract Class to create strut-based lattice""" _UNIT_CUBE_SIZE = 1.0 @@ -136,8 +134,10 @@ def calc_density(radius: float) -> float: @property def base_vertices(self) -> npt.NDArray[np.float64]: - """Property: coordinates of the vertices for a structure - centered at the origin and enclosed in a size 1 cubic rve""" + """ + Property: coordinates of the vertices for a structure + centered at the origin and enclosed in a size 1 cubic rve + """ if self._base_vertices is not None: return self._base_vertices return self._generate_base_vertices() @@ -151,15 +151,14 @@ def strut_vertex_pairs(self) -> npt.NDArray[np.int64]: @abstractmethod def _generate_base_vertices(self) -> npt.NDArray[np.float64]: - """Abstract method to generate base vertices, ie as if the + """ + Abstract method to generate base vertices, ie as if the lattice was centered at the origin and in a cubic size 1 rve. """ - pass @abstractmethod def _generate_strut_vertex_pairs(self) -> npt.NDArray[np.int64]: """Abstract method to generate strut vertex pairs.""" - pass def _compute_vertices(self) -> npt.NDArray[np.float64]: return self.center + self.cell_size * self.base_vertices @@ -173,7 +172,6 @@ def _compute_strut_directions(self) -> npt.NDArray[np.float64]: def _validate_inputs(self): """Checks coherence of inputs.""" - if self._strut_heights is None: raise NotImplementedError("strut_heights must be defined by the subclass") if ( @@ -181,7 +179,7 @@ def _validate_inputs(self): and len(self._strut_heights) != self.strut_number ): raise ValueError( - f"strut_heights must contain {self.strut_number} values, but {len(self._strut_heights)} were provided." + f"strut_heights must contain {self.strut_number} values, but {len(self._strut_heights)} were provided.", ) @property @@ -200,22 +198,24 @@ def strut_heights(self) -> list[float]: return self._strut_heights * self.cell_size def _compute_rotations(self) -> list[Rotation]: - """Computes rotation from default (1.0, 0.0, 0.0) oriented Cylinder - for all struts in the lattice using Scipy's Rotation object.""" - + """ + Computes rotation from default (1.0, 0.0, 0.0) oriented Cylinder + for all struts in the lattice using Scipy's Rotation object. + """ default_direction = np.array([1.0, 0.0, 0.0]) rotations_list = [] for i in range(self.strut_number): if np.all( - self.strut_directions_cartesian[i] == default_direction + self.strut_directions_cartesian[i] == default_direction, ) or np.all(self.strut_directions_cartesian[i] == -default_direction): rotation_vector = np.zeros(3) rotations_list.append(Rotation.from_rotvec(rotation_vector)) else: rotation, _ = Rotation.align_vectors( - self.strut_directions_cartesian[i], default_direction + self.strut_directions_cartesian[i], + default_direction, ) rotations_list.append(rotation) @@ -350,7 +350,7 @@ def generateVtk( # noqa: N802 periodic: bool = True, **kwargs: KwargsGenerateType, ) -> pv.PolyData: - """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + """Deprecated. Use :meth:`generate_vtk` instead.""" return self.generate_vtk( size=size, order=order, diff --git a/microgen/shape/strut_lattice/body_centered_cubic.py b/microgen/shape/strut_lattice/body_centered_cubic.py index 005693f3..ae151bb9 100644 --- a/microgen/shape/strut_lattice/body_centered_cubic.py +++ b/microgen/shape/strut_lattice/body_centered_cubic.py @@ -37,8 +37,11 @@ def __init__(self, *args, **kwargs) -> None: def _generate_base_vertices(self) -> npt.NDArray[np.float64]: unit_cube_vertices = np.array( list( - product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3) - ) + product( + [-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], + repeat=3, + ), + ), ) return np.vstack(([0.0, 0.0, 0.0], unit_cube_vertices)) diff --git a/microgen/shape/strut_lattice/cubic.py b/microgen/shape/strut_lattice/cubic.py index 62edc4f0..3f64fc83 100644 --- a/microgen/shape/strut_lattice/cubic.py +++ b/microgen/shape/strut_lattice/cubic.py @@ -38,8 +38,11 @@ def __init__(self, *args, **kwargs) -> None: def _generate_base_vertices(self) -> npt.NDArray[np.float64]: return np.array( list( - product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3) - ) + product( + [-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], + repeat=3, + ), + ), ) def _generate_strut_vertex_pairs(self) -> npt.NDArray[np.int64]: @@ -48,8 +51,9 @@ def _generate_strut_vertex_pairs(self) -> npt.NDArray[np.int64]: [i, j] for i in range(len(self.base_vertices)) for j in KDTree(self.base_vertices).query_ball_point( - self.base_vertices[i], r=self._UNIT_CUBE_SIZE + self.base_vertices[i], + r=self._UNIT_CUBE_SIZE, ) if i < j - ] + ], ) diff --git a/microgen/shape/strut_lattice/cuboctahedron.py b/microgen/shape/strut_lattice/cuboctahedron.py index a2038ec0..5abe9e41 100644 --- a/microgen/shape/strut_lattice/cuboctahedron.py +++ b/microgen/shape/strut_lattice/cuboctahedron.py @@ -38,8 +38,11 @@ def __init__(self, *args, **kwargs) -> None: def _generate_base_vertices(self) -> npt.NDArray[np.float64]: cube_vertices = np.array( list( - product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3) - ) + product( + [-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], + repeat=3, + ), + ), ) edges = [ @@ -59,7 +62,7 @@ def _generate_strut_vertex_pairs(self) -> npt.NDArray[np.int64]: tree.query_ball_point( self.base_vertices, r=self._UNIT_CUBE_SIZE / np.sqrt(2.0) + BALL_POINT_RADIUS_TOLERANCE, - ) + ), ): for j in indices: if i != j: diff --git a/microgen/shape/strut_lattice/custom_lattice.py b/microgen/shape/strut_lattice/custom_lattice.py index ed86fbe8..db277047 100644 --- a/microgen/shape/strut_lattice/custom_lattice.py +++ b/microgen/shape/strut_lattice/custom_lattice.py @@ -11,9 +11,7 @@ class CustomLattice(AbstractLattice): - """ - Class to create a custom lattice with user-defined base vertices and strut vertex pairs - """ + """Class to create a custom lattice with user-defined base vertices and strut vertex pairs""" def __init__( self, diff --git a/microgen/shape/strut_lattice/diamond.py b/microgen/shape/strut_lattice/diamond.py index cefa207d..b92f2f5b 100644 --- a/microgen/shape/strut_lattice/diamond.py +++ b/microgen/shape/strut_lattice/diamond.py @@ -40,8 +40,11 @@ def __init__(self, *args, **kwargs) -> None: def _generate_tetrahedra_centers(self) -> npt.NDArray[np.float64]: candidates = np.array( list( - product([-self._UNIT_CUBE_SIZE / 4, self._UNIT_CUBE_SIZE / 4], repeat=3) - ) + product( + [-self._UNIT_CUBE_SIZE / 4, self._UNIT_CUBE_SIZE / 4], + repeat=3, + ), + ), ) centers = candidates[np.sum(candidates < 0, axis=1) % 2 == 0] @@ -55,8 +58,11 @@ def _generate_tetrahedra_vertices(self) -> npt.NDArray[np.float64]: ] candidates = np.array( list( - product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3) - ) + product( + [-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], + repeat=3, + ), + ), ) outer_cube_corners = candidates[np.sum(candidates < 0, axis=1) % 2 == 0] diff --git a/microgen/shape/strut_lattice/face_centered_cubic.py b/microgen/shape/strut_lattice/face_centered_cubic.py index c4eb87db..cecf55ef 100644 --- a/microgen/shape/strut_lattice/face_centered_cubic.py +++ b/microgen/shape/strut_lattice/face_centered_cubic.py @@ -5,7 +5,6 @@ """ from itertools import product -from typing import List import numpy as np import numpy.typing as npt @@ -37,7 +36,7 @@ def __init__(self, *args, **kwargs) -> None: def _generate_base_vertices(self) -> npt.NDArray[np.float64]: cube_vertices = list( - product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3) + product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3), ) face_centers = [ @@ -48,10 +47,8 @@ def _generate_base_vertices(self) -> npt.NDArray[np.float64]: return np.array(cube_vertices + face_centers) - def _generate_face_center_to_cube_vertices_dict(self) -> dict[int, List[int]]: - """ - Dynamically generates a dictionary associating the indices of the face centers with the indices of the cube vertices. - """ + def _generate_face_center_to_cube_vertices_dict(self) -> dict[int, list[int]]: + """Dynamically generates a dictionary associating the indices of the face centers with the indices of the cube vertices.""" num_cube_vertices = 8 cube_vertices = self.base_vertices[:num_cube_vertices] face_centers = self.base_vertices[num_cube_vertices:] @@ -73,9 +70,7 @@ def _generate_face_center_to_cube_vertices_dict(self) -> dict[int, List[int]]: return face_to_vertices def _generate_strut_vertex_pairs(self) -> npt.NDArray[np.int64]: - """ - Generates the pairs of indices of the fcc strut vertices between the face centers and the cube vertices. - """ + """Generates the pairs of indices of the fcc strut vertices between the face centers and the cube vertices.""" face_to_vertices = self._generate_face_center_to_cube_vertices_dict() vertex_pairs_indices = [ diff --git a/microgen/shape/strut_lattice/octahedron.py b/microgen/shape/strut_lattice/octahedron.py index 115e7149..06917031 100644 --- a/microgen/shape/strut_lattice/octahedron.py +++ b/microgen/shape/strut_lattice/octahedron.py @@ -42,7 +42,7 @@ def _generate_base_vertices(self) -> npt.NDArray[np.float64]: ] for axis in range(3) for sign in [-1, 1] - ] + ], ) def _generate_strut_vertex_pairs(self) -> npt.NDArray[np.int64]: diff --git a/microgen/shape/strut_lattice/octet_truss.py b/microgen/shape/strut_lattice/octet_truss.py index 156dddf4..1439a5c4 100644 --- a/microgen/shape/strut_lattice/octet_truss.py +++ b/microgen/shape/strut_lattice/octet_truss.py @@ -37,7 +37,7 @@ def __init__(self, *args, **kwargs) -> None: def _generate_base_vertices(self) -> npt.NDArray[np.float64]: cube_vertices = list( - product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3) + product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3), ) face_centers = [ diff --git a/microgen/shape/strut_lattice/rhombic_cuboctahedron.py b/microgen/shape/strut_lattice/rhombic_cuboctahedron.py index fcb9db4b..71d48d32 100644 --- a/microgen/shape/strut_lattice/rhombic_cuboctahedron.py +++ b/microgen/shape/strut_lattice/rhombic_cuboctahedron.py @@ -38,14 +38,18 @@ def __init__(self, *args, **kwargs) -> None: def _generate_base_vertices(self) -> npt.NDArray[np.float64]: permutations_set = set( permutations( - [self._UNIT_CUBE_SIZE / 2.0, (np.sqrt(2) - 1) / 2, (np.sqrt(2) - 1) / 2] - ) + [ + self._UNIT_CUBE_SIZE / 2.0, + (np.sqrt(2) - 1) / 2, + (np.sqrt(2) - 1) / 2, + ], + ), ) vertices = [] for permutation in permutations_set: for signs in product([-1, 1], repeat=3): - vertex = tuple(s * p for s, p in zip(signs, permutation)) + vertex = tuple(s * p for s, p in zip(signs, permutation, strict=False)) vertices.append(vertex) return np.array(vertices) diff --git a/microgen/shape/strut_lattice/rhombic_dodecahedron.py b/microgen/shape/strut_lattice/rhombic_dodecahedron.py index 5e99ad8e..4cd010d1 100644 --- a/microgen/shape/strut_lattice/rhombic_dodecahedron.py +++ b/microgen/shape/strut_lattice/rhombic_dodecahedron.py @@ -37,10 +37,10 @@ def __init__(self, *args, **kwargs) -> None: def _generate_base_vertices(self) -> npt.NDArray[np.float64]: outer_cube_vertices = list( - product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3) + product([-self._UNIT_CUBE_SIZE / 2, self._UNIT_CUBE_SIZE / 2], repeat=3), ) inner_cube_vertices = list( - product([-self._UNIT_CUBE_SIZE / 4, self._UNIT_CUBE_SIZE / 4], repeat=3) + product([-self._UNIT_CUBE_SIZE / 4, self._UNIT_CUBE_SIZE / 4], repeat=3), ) face_centers = [ diff --git a/microgen/shape/strut_lattice/truncated_cuboctahedron.py b/microgen/shape/strut_lattice/truncated_cuboctahedron.py index cda156f1..1207c400 100644 --- a/microgen/shape/strut_lattice/truncated_cuboctahedron.py +++ b/microgen/shape/strut_lattice/truncated_cuboctahedron.py @@ -53,7 +53,8 @@ def _generate_strut_vertex_pairs(self) -> npt.NDArray[np.int64]: pairs = [] for i in range(len(self.base_vertices)): neighbors = kdtree.query_ball_point( - self.base_vertices[i], threshold_distance + self.base_vertices[i], + threshold_distance, ) for j in neighbors: if i < j: diff --git a/microgen/shape/strut_lattice/truncated_octahedron.py b/microgen/shape/strut_lattice/truncated_octahedron.py index d730540f..cdad010c 100644 --- a/microgen/shape/strut_lattice/truncated_octahedron.py +++ b/microgen/shape/strut_lattice/truncated_octahedron.py @@ -38,10 +38,10 @@ def __init__(self, *args, **kwargs) -> None: def _generate_base_vertices(self) -> npt.NDArray[np.float64]: base_vertices = set() for perm in permutations( - [0, self._UNIT_CUBE_SIZE / 4.0, self._UNIT_CUBE_SIZE / 2.0] + [0, self._UNIT_CUBE_SIZE / 4.0, self._UNIT_CUBE_SIZE / 2.0], ): for signs in product([-1, 1], repeat=3): - vertex = tuple(s * p for s, p in zip(signs, perm)) + vertex = tuple(s * p for s, p in zip(signs, perm, strict=False)) if sum(abs(v) for v in vertex) == 3.0 / 4.0: base_vertices.add(vertex) return np.array(list(base_vertices)) diff --git a/microgen/shape/surface_functions.py b/microgen/shape/surface_functions.py index 1d8cf2a3..6d9c95b2 100644 --- a/microgen/shape/surface_functions.py +++ b/microgen/shape/surface_functions.py @@ -5,7 +5,8 @@ def gyroid(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Gyroid. + """ + Gyroid. .. math:: sin(x) cos(y) + sin(y) cos(z) + sin(z) cos(x) = 0 @@ -28,7 +29,8 @@ def gyroid(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def schwarz_p(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Schwarz P. + """ + Schwarz P. .. math:: cos(x) + cos(y) + cos(z) = 0 @@ -51,7 +53,8 @@ def schwarz_p(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def schwarz_d(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Schwarz D. + """ + Schwarz D. .. math:: sin(x) sin(y) sin(z) + sin(x) cos(y) cos(z) +\ @@ -80,7 +83,8 @@ def schwarz_d(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def neovius(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Neovius. + """ + Neovius. .. math:: 3 cos(x) + cos(y) + cos(z) + 4 cos(x) cos(y) cos(z) = 0 @@ -103,7 +107,8 @@ def neovius(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def schoen_iwp(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Schoen IWP. + """ + Schoen IWP. .. math:: 2 (cos(x) cos(y) + cos(y) cos(z) + cos(z) cos(x)) \ @@ -129,7 +134,8 @@ def schoen_iwp(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def schoen_frd(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Schoen FRD. + """ + Schoen FRD. .. math:: 4 cos(x) cos(y) cos(z) \ @@ -155,7 +161,8 @@ def schoen_frd(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def fischer_koch_s(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Fischer-Koch S. + """ + Fischer-Koch S. .. math:: cos(2 x) sin(y) cos(z) + cos(x) cos(2 y) sin(z) + sin(x) cos(y) cos(2 z) = 0 @@ -182,7 +189,8 @@ def fischer_koch_s(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def pmy(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """PMY. + """ + PMY. .. math:: 2 cos(x) cos(y) cos(z) + sin(2 x) sin(y) + sin(x) sin(2 z) + sin(2 y) sin(z) = 0 @@ -210,7 +218,8 @@ def pmy(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def honeycomb(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Honneycomb. + """ + Honneycomb. .. math:: sin(x) cos(y) + sin(y) + cos(z) = 0 @@ -233,7 +242,8 @@ def honeycomb(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def lidinoid(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Lidinoid. + """ + Lidinoid. .. math:: 0.5 (sin(2 x) cos(y) sin(z) + @@ -271,7 +281,8 @@ def lidinoid(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def split_p(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """Split P. + """ + Split P. .. math:: 1.1 (sin(2 x) cos(y) sin(z) + @@ -310,7 +321,8 @@ def split_p(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: def honeycomb_gyroid(x: float, y: float, _: float) -> float: - """Honeycomb Gyroid. + """ + Honeycomb Gyroid. .. math:: sin(x) cos(y) + sin(y) + cos(x) = 0 @@ -333,7 +345,8 @@ def honeycomb_gyroid(x: float, y: float, _: float) -> float: def honeycomb_schwarz_p(x: float, y: float, _: float) -> float: - """Honeycomb Schwarz P. + """ + Honeycomb Schwarz P. .. math:: cos(x) + cos(y) = 0 @@ -356,7 +369,8 @@ def honeycomb_schwarz_p(x: float, y: float, _: float) -> float: def honeycomb_schwarz_d(x: float, y: float, _: float) -> float: - """Honneycomb Schwarz D. + """ + Honneycomb Schwarz D. .. math:: cos(x) cos(y) + sin(x) sin(y) + sin(x) cos(y) + cos(x) sin(y) = 0 @@ -379,7 +393,8 @@ def honeycomb_schwarz_d(x: float, y: float, _: float) -> float: def honeycomb_schoen_iwp(x: float, y: float, _: float) -> float: - """Honneycomb Schoen IWP. + """ + Honneycomb Schoen IWP. .. math:: cos(x) cos(y) + cos(y) + cos(x) = 0 @@ -402,7 +417,8 @@ def honeycomb_schoen_iwp(x: float, y: float, _: float) -> float: def honeycomb_lidinoid(x: float, y: float, _: float) -> float: - """Honeycomb Lidinoid. + """ + Honeycomb Lidinoid. .. math:: 1.1 (sin(2 x) cos(y) + sin(2 y) sin(x) + cos(x) sin(y)) diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 0a5b0cfd..9f98c1fc 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -1,4 +1,5 @@ -"""TPMS. +""" +TPMS. ============================================= TPMS (:mod:`microgen.shape.tpms`) @@ -15,7 +16,8 @@ from __future__ import annotations import logging -from typing import TYPE_CHECKING, Callable, Literal, Sequence +from collections.abc import Callable, Sequence +from typing import TYPE_CHECKING, Literal import cadquery as cq import numpy as np @@ -45,7 +47,8 @@ def _wedge_sdf_2d( y: npt.NDArray[np.float64], half_angle: float, ) -> npt.NDArray[np.float64]: - """SDF of a 2D wedge centred on +X with half-aperture ``half_angle``. + """ + SDF of a 2D wedge centred on +X with half-aperture ``half_angle``. Negative inside, positive outside, zero on the bounding rays. Z-independent — used as a clipping primitive for cylindrical / spherical @@ -60,7 +63,8 @@ def _double_cone_sdf( z: npt.NDArray[np.float64], half_polar: float, ) -> npt.NDArray[np.float64]: - """SDF of a double cone around the Z axis with half-aperture ``half_polar``. + """ + SDF of a double cone around the Z axis with half-aperture ``half_polar``. Symmetric about the XY plane; points inside iff their angle from ±Z is less than ``half_polar``. Used to clip spherical partial-θ coverage. @@ -74,7 +78,8 @@ def _interp_along_curve( curve_s: npt.NDArray[np.float64], values: npt.NDArray[np.float64], ) -> npt.NDArray[np.float64]: - """Per-axis ``np.interp`` of an ``(M, 3)`` curve attribute at ``query_s``. + """ + Per-axis ``np.interp`` of an ``(M, 3)`` curve attribute at ``query_s``. Returns an ``(N, 3)`` array. Wraps the three-axis loop so callers don't repeat the pattern. @@ -85,7 +90,8 @@ def _interp_along_curve( class Tpms(Shape): - """Triply Periodical Minimal Surfaces. + """ + Triply Periodical Minimal Surfaces. Class to generate Triply Periodical Minimal Surfaces (TPMS) geometry from a given mathematical function, with given offset @@ -109,7 +115,7 @@ class Tpms(Shape): - :class:`~microgen.shape.surface_functions.honeycomb_lidinoid` """ - def __init__( # noqa: PLR0913 + def __init__( self: Tpms, surface_function: Field, offset: float | OffsetGrading | Field | None = None, @@ -120,7 +126,8 @@ def __init__( # noqa: PLR0913 density: float | None = None, **kwargs: Vector3DType, ) -> None: - r"""Class used to generate TPMS geometries (sheet or skeletals parts). + r""" + Class used to generate TPMS geometries (sheet or skeletals parts). TPMS are created by default in a cube. The geometry of the cube can be modified using 'cell_size' parameter. @@ -213,7 +220,8 @@ def offset_from_density( density: float, resolution: int = 20, ) -> float: - """Return the offset corresponding to the required density. + """ + Return the offset corresponding to the required density. :param surface_function: tpms function :param part_type: type of the part (sheet, lower skeletal or upper skeletal) @@ -222,13 +230,14 @@ def offset_from_density( :return: corresponding offset value """ - return Tpms( # noqa: SLF001 + return Tpms( surface_function=surface_function, density=density, )._compute_offset_to_fit_density(part_type=part_type, resolution=resolution) def _density_envelope_volume(self: Tpms) -> float: - """Return the *envelope* volume that ``density`` is measured against. + """ + Return the *envelope* volume that ``density`` is measured against. For a plain :class:`Tpms` this is the cartesian cell volume. Subclasses with a non-cartesian envelope (:class:`Infill`, :class:`Conformal`) @@ -243,7 +252,8 @@ def _compute_offset_to_fit_density( part_type: Literal["sheet", "lower skeletal", "upper skeletal"], resolution: int | None = None, ) -> float: - """Compute the offset that yields the requested density. + """ + Compute the offset that yields the requested density. Searches with the same F-rep ``generate_vtk`` pipeline the user invokes, so the offset returned actually reproduces the requested @@ -432,6 +442,7 @@ def _compute_tpms_field(self: Tpms) -> None: for repeat_cell_axis, cell_size_axis in zip( self.repeat_cell, self.cell_size, + strict=False, ) ] @@ -455,7 +466,7 @@ def _finalize_frep( bounds: tuple[float, float, float, float, float, float], ) -> None: """Normalize a raw field to SDF and set ``_func`` / ``_bounds``.""" - from .implicit_ops import from_field, normalize_to_sdf # noqa: PLC0415 + from .implicit_ops import from_field, normalize_to_sdf self._raw_field_func = raw_field sdf_shape = normalize_to_sdf(from_field(raw_field)) @@ -490,15 +501,16 @@ def raw_field(self: Tpms) -> Field: return self._raw_field_func def as_sheet(self: Tpms, thickness: float | None = None) -> Shape: - """Return an F-rep Shape representing a TPMS sheet of given thickness. + """ + Return an F-rep Shape representing a TPMS sheet of given thickness. Uses the SDF-normalized field, so *thickness* is in physical units. If *thickness* is ``None``, uses ``self.offset`` (which may be a scalar, an array sampled on ``self.grid``, or a callable in which case the callable form is used directly). """ - from .implicit_ops import shell # noqa: PLC0415 - from .shape import Shape # noqa: PLC0415 + from .implicit_ops import shell + from .shape import Shape if thickness is not None: t: float | npt.NDArray | Field = float(thickness) @@ -509,7 +521,8 @@ def as_sheet(self: Tpms, thickness: float | None = None) -> Shape: return shell(Shape(func=self._func, bounds=self._bounds), t) def _half_offset_field(self: Tpms) -> Field | float: - """Return half the offset as a callable (variable) or scalar (constant). + """ + Return half the offset as a callable (variable) or scalar (constant). Variable offset stored as a callable can be re-evaluated on the marching-cubes grid; an array offset (sampled on ``self.grid``) @@ -520,7 +533,7 @@ def _half_offset_field(self: Tpms) -> Field | float: if self._offset_func is not None: f = self._offset_func - def _half(x, y, z, _f=f): # noqa: ANN001 + def _half(x, y, z, _f=f): return 0.5 * _f(x, y, z) return _half @@ -530,84 +543,87 @@ def _half(x, y, z, _f=f): # noqa: ANN001 return 0.0 def as_upper_skeletal(self: Tpms) -> Shape: - """F-rep Shape for the *upper* skeletal: ``{p : f(p) > offset/2}``. + """ + F-rep Shape for the *upper* skeletal: ``{p : f(p) > offset/2}``. Volume scales with the chosen offset (smaller offset ⇒ larger skeletal), matching the historical CadQuery behaviour and the VTK grid-clip path. """ - from .implicit_ops import from_field # noqa: PLC0415 + from .implicit_ops import from_field f = self._func h = self._half_offset_field() if callable(h): - def _upper(x, y, z, _f=f, _h=h): # noqa: ANN001 + def _upper(x, y, z, _f=f, _h=h): return -_f(x, y, z) + _h(x, y, z) else: - def _upper(x, y, z, _f=f, _h=h): # noqa: ANN001 + def _upper(x, y, z, _f=f, _h=h): return -_f(x, y, z) + _h return from_field(func=_upper, bounds=self._bounds) def as_lower_skeletal(self: Tpms) -> Shape: """F-rep Shape for the *lower* skeletal: ``{p : f(p) < -offset/2}``.""" - from .implicit_ops import from_field # noqa: PLC0415 + from .implicit_ops import from_field f = self._func h = self._half_offset_field() if callable(h): - def _lower(x, y, z, _f=f, _h=h): # noqa: ANN001 + def _lower(x, y, z, _f=f, _h=h): return _f(x, y, z) + _h(x, y, z) else: - def _lower(x, y, z, _f=f, _h=h): # noqa: ANN001 + def _lower(x, y, z, _f=f, _h=h): return _f(x, y, z) + _h return from_field(func=_lower, bounds=self._bounds) def as_surface(self: Tpms) -> Shape: - """Return F-rep Shape for the (open) zero-isosurface, no thickness. + """ + Return F-rep Shape for the (open) zero-isosurface, no thickness. Same field as :meth:`as_lower_skeletal`; meant for ``type_part="surface"``. Marching cubes will produce an open shell — there is no enclosed volume. """ - from .implicit_ops import from_field # noqa: PLC0415 + from .implicit_ops import from_field return from_field(func=self._func, bounds=self._bounds) def _cell_box(self: Tpms) -> Shape: """SDF Shape of this TPMS' cell (cell_size × repeat_cell, centered origin).""" - from .implicit_ops import box # noqa: PLC0415 + from .implicit_ops import box dims = tuple(float(d) for d in (self.cell_size * self.repeat_cell)) return box(dims=dims, center=(0.0, 0.0, 0.0)) def _clipped_sheet(self: Tpms) -> Shape: - """Sheet F-rep clipped to the cell box. + """ + Sheet F-rep clipped to the cell box. When the offset approaches its upper limit the sheet field is uniformly negative inside the cell and marching cubes finds no boundary — clipping by the cell box makes the box face the closed boundary, yielding a full-cell mesh as expected. """ - from .implicit_ops import intersection # noqa: PLC0415 + from .implicit_ops import intersection return intersection(self.as_sheet(), self._cell_box()) def _clipped_upper_skeletal(self: Tpms) -> Shape: """Upper skeletal F-rep clipped to the cell box (closed under marching cubes).""" - from .implicit_ops import intersection # noqa: PLC0415 + from .implicit_ops import intersection return intersection(self.as_upper_skeletal(), self._cell_box()) def _clipped_lower_skeletal(self: Tpms) -> Shape: """Lower skeletal F-rep clipped to the cell box (closed under marching cubes).""" - from .implicit_ops import intersection # noqa: PLC0415 + from .implicit_ops import intersection return intersection(self.as_lower_skeletal(), self._cell_box()) @@ -653,7 +669,8 @@ def offset( self.offset_updated = True def _mesh_to_shell(self: Tpms, mesh: pv.PolyData) -> cq.Shape: - """Convert a triangulated PyVista mesh to a CadQuery ``Shape``. + """ + Convert a triangulated PyVista mesh to a CadQuery ``Shape``. Builds one ``cq.Face`` per triangle, then *sews* them with OCCT's ``BRepBuilderAPI_Sewing`` so adjacent triangles share edges — without @@ -662,7 +679,7 @@ def _mesh_to_shell(self: Tpms, mesh: pv.PolyData) -> cq.Shape: ``TopoDS_Shell`` if sewing succeeded into one shell, otherwise the compound that sewing produced). """ - from OCP.BRepBuilderAPI import BRepBuilderAPI_Sewing # noqa: PLC0415 + from OCP.BRepBuilderAPI import BRepBuilderAPI_Sewing if not mesh.is_all_triangles: mesh.triangulate(inplace=True) @@ -676,7 +693,7 @@ def _mesh_to_shell(self: Tpms, mesh: pv.PolyData) -> cq.Shape: cq.Vector(*mesh.points[start]), cq.Vector(*mesh.points[end]), ) - for start, end in zip(tri[:], tri[1:]) + for start, end in zip(tri[:], tri[1:], strict=False) ] wire = cq.Wire.assembleEdges(lines) faces.append(cq.Face.makeFromWires(wire)) @@ -729,7 +746,8 @@ def _check_offset( _VALID_PARTS = ("sheet", "lower skeletal", "upper skeletal", "surface") def _frep_part(self: Tpms, type_part: TpmsPartType) -> Shape: - """Pick the F-rep :class:`Shape` for *type_part*. + """ + Pick the F-rep :class:`Shape` for *type_part*. Skeletals are intersected with the cell box so marching cubes produces a *closed* shell (the unclipped skeletal field is @@ -748,7 +766,8 @@ def _frep_part(self: Tpms, type_part: TpmsPartType) -> Shape: raise ValueError(err_msg) def _isotropic_resolution(self: Tpms) -> int: - """Map ``self.resolution`` (per-axis) to an isotropic Shape resolution. + """ + Map ``self.resolution`` (per-axis) to an isotropic Shape resolution. ``Shape.generate_vtk`` takes a single resolution; we use the geometric mean of the per-axis cell counts so total grid points stay proportional. @@ -756,7 +775,8 @@ def _isotropic_resolution(self: Tpms) -> int: return max(int(self.resolution * np.cbrt(np.prod(self.repeat_cell))), 10) def _envelope_mesh_at_full_density(self: Tpms) -> pv.PolyData: - """Return the cell-envelope mesh used by the ``density=1.0`` shortcut. + """ + Return the cell-envelope mesh used by the ``density=1.0`` shortcut. Plain :class:`Tpms` has an axis-aligned box envelope, returned exactly via :class:`pyvista.Box`. Subclasses with non-cubic envelopes @@ -772,7 +792,8 @@ def _envelope_mesh_at_full_density(self: Tpms) -> pv.PolyData: ) def _envelope_mesh_via_cell_box(self: Tpms) -> pv.PolyData: - """Marching-cubes mesh of the (potentially non-cubic) ``_cell_box``. + """ + Marching-cubes mesh of the (potentially non-cubic) ``_cell_box``. Used by subclasses that override :meth:`_envelope_mesh_at_full_density` to delegate to the F-rep envelope SDF. @@ -790,7 +811,8 @@ def generate( algo_resolution: int | None = None, **_: KwargsGenerateType, ) -> cq.Shape: - """Generate the OCCT/CadQuery shape of the requested TPMS part. + """ + Generate the OCCT/CadQuery shape of the requested TPMS part. Pure F-rep pipeline: pick the SDF Shape via :meth:`_frep_part`, run marching cubes through :meth:`Shape.generate_vtk`, optionally smooth, @@ -853,15 +875,16 @@ def generate( @staticmethod def _try_make_solid(shape: cq.Shape) -> cq.Shape: - """Best-effort upgrade of a sewn shell into a closed Solid. + """ + Best-effort upgrade of a sewn shell into a closed Solid. Returns the original shape unchanged if the sewn result is a Compound (multiple disjoint shells, can't be a single solid) or if OCCT refuses the conversion. """ - from OCP.TopAbs import TopAbs_SHELL # noqa: PLC0415 - from OCP.TopExp import TopExp_Explorer # noqa: PLC0415 - from OCP.TopoDS import TopoDS # noqa: PLC0415 + from OCP.TopAbs import TopAbs_SHELL + from OCP.TopExp import TopExp_Explorer + from OCP.TopoDS import TopoDS wrapped = shape.wrapped # Already a Shell? Try to make a Solid directly. @@ -896,7 +919,8 @@ def generate_vtk( algo_resolution: int | None = None, **_: KwargsGenerateType, ) -> pv.PolyData: - """Generate the PyVista mesh of the requested TPMS part. + """ + Generate the PyVista mesh of the requested TPMS part. Same F-rep pipeline as :meth:`generate` (skeletals are clipped to the cell box), so the two outputs share the exact same triangulation and @@ -997,7 +1021,7 @@ def generateVtk( # noqa: N802 algo_resolution: int | None = None, **_: KwargsGenerateType, ) -> pv.PolyData: - """Deprecated. Use :meth:`generate_vtk` instead.""" # noqa: D401 + """Deprecated. Use :meth:`generate_vtk` instead.""" return self.generate_vtk( type_part=type_part, algo_resolution=algo_resolution, @@ -1016,7 +1040,7 @@ class CylindricalTpms(Tpms): _envelope_mesh_at_full_density = Tpms._envelope_mesh_via_cell_box - def __init__( # noqa: PLR0913 + def __init__( self: CylindricalTpms, radius: float, surface_function: Field, @@ -1029,7 +1053,8 @@ def __init__( # noqa: PLR0913 resolution: int = 20, density: float | None = None, ) -> None: - r"""Cylindrical TPMS geometry. + r""" + Cylindrical TPMS geometry. Directions of cell_size and repeat_cell must be taken as the cylindrical \ coordinate system $\\left(\\rho, \\theta, z\\right)$. @@ -1087,7 +1112,8 @@ def _create_grid( y: npt.NDArray[np.float64], z: npt.NDArray[np.float64], ) -> pv.StructuredGrid: - """Return the structured cylindrical grid of the TPMS. + """ + Return the structured cylindrical grid of the TPMS. ``y`` carries arc-length along the equator, so the conversion to an angle is ``theta = y / radius``. The earlier ``y * unit_theta`` form @@ -1145,7 +1171,8 @@ def _raw_field( ) def _cell_box(self: CylindricalTpms) -> Shape: - """Cylindrical-shell SDF in Cartesian space (replaces the parent's + """ + Cylindrical-shell SDF in Cartesian space (replaces the parent's axis-aligned-box SDF, which was geometrically wrong for this class). Without this override the F-rep marching cubes path @@ -1164,7 +1191,7 @@ def _cell_box(self: CylindricalTpms) -> Shape: handled by intersecting the ring with an angular wedge centred on the +X axis. """ - from .implicit_ops import from_field, intersection # noqa: PLC0415 + from .implicit_ops import from_field, intersection r_inner, r_outer, half_z = self._shell_extents() bounds: BoundsType = ( @@ -1204,7 +1231,8 @@ def _shell_sdf( return shell def _isotropic_resolution(self: CylindricalTpms) -> int: - """Per-axis resolution count for the F-rep marching-cubes grid. + """ + Per-axis resolution count for the F-rep marching-cubes grid. Override the parent's geometric-mean of ``repeat_cell`` (which mixes the angular axis with linear ones). Use only the radial (``[0]``) @@ -1230,7 +1258,7 @@ class SphericalTpms(Tpms): _envelope_mesh_at_full_density = Tpms._envelope_mesh_via_cell_box - def __init__( # noqa: PLR0913 + def __init__( self: SphericalTpms, radius: float, surface_function: Field, @@ -1243,7 +1271,8 @@ def __init__( # noqa: PLR0913 resolution: int = 20, density: float | None = None, ) -> None: - r"""Spherical TPMS geometry. + r""" + Spherical TPMS geometry. Directions of cell_size and repeat_cell must be taken as the spherical \ coordinate system $\\left(r, \\theta, \\phi\\right)$. @@ -1307,7 +1336,8 @@ def _create_grid( y: npt.NDArray[np.float64], z: npt.NDArray[np.float64], ) -> pv.StructuredGrid: - """Return the structured spherical grid of the TPMS. + """ + Return the structured spherical grid of the TPMS. ``y`` and ``z`` carry equatorial arc length, so the conversion to angles is ``theta = y / radius`` and ``phi = z / radius``. The @@ -1372,7 +1402,8 @@ def _shell_radii(self: SphericalTpms) -> tuple[float, float]: return max(sph_r - delta_r, 0.0), sph_r + delta_r def _cell_box(self: SphericalTpms) -> Shape: - """Spherical-shell SDF in Cartesian space. + """ + Spherical-shell SDF in Cartesian space. Replaces the parent's axis-aligned-box clip — same fix as :meth:`CylindricalTpms._cell_box`, see that docstring for the @@ -1382,7 +1413,7 @@ def _cell_box(self: SphericalTpms) -> Shape: in xy-plane) is supported when ``repeat_cell[1]`` / ``repeat_cell[2]`` are smaller than the auto-fill values. """ - from .implicit_ops import from_field, intersection # noqa: PLC0415 + from .implicit_ops import from_field, intersection r_inner, r_outer = self._shell_radii() bounds: BoundsType = ( @@ -1434,7 +1465,8 @@ def _shell_sdf( return shape def _isotropic_resolution(self: SphericalTpms) -> int: - """Use only the radial physical extent (the angular axes are + """ + Use only the radial physical extent (the angular axes are parametric and would mislead a geometric mean across them). """ radial_extent = float(self.cell_size[0]) * float(self.repeat_cell[0]) @@ -1445,7 +1477,8 @@ def _isotropic_resolution(self: SphericalTpms) -> int: class Sweep(Tpms): - """TPMS along an arbitrary 3D curve — generalisation of CylindricalTpms. + """ + TPMS along an arbitrary 3D curve — generalisation of CylindricalTpms. The TPMS is generated in a tube of radius ``radial_max`` around an arbitrary curve, in local coordinates ``(s, r, θ)`` where: @@ -1477,7 +1510,7 @@ class Sweep(Tpms): _envelope_mesh_at_full_density = Tpms._envelope_mesh_via_cell_box - def __init__( # noqa: PLR0913 + def __init__( self: Sweep, curve_points: npt.NDArray[np.float64] | Callable[[float], npt.NDArray[np.float64]], @@ -1494,7 +1527,8 @@ def __init__( # noqa: PLR0913 center: Vector3DType = (0, 0, 0), orientation: Vector3DType = (0, 0, 0), ) -> None: - r"""Build a TPMS swept along a curve. + r""" + Build a TPMS swept along a curve. :param curve_points: either an ``(M, 3)`` array of polyline samples or a callable ``t \in [0, 1] -> (3,)``. Callables are sampled @@ -1558,7 +1592,8 @@ def _build_curve_frames( self: Sweep, seed_normal: Sequence[float] | None, ) -> None: - """Pre-compute arc length, tangents, and parallel-transported normals. + """ + Pre-compute arc length, tangents, and parallel-transported normals. Stores on ``self``: ``_curve_s`` arc length per sample, shape (M,) @@ -1566,7 +1601,7 @@ def _build_curve_frames( ``_curve_N`` unit normal (parallel-transported), shape (M, 3) ``_curve_kdtree`` scipy cKDTree on the curve points """ - from scipy.spatial import cKDTree # noqa: PLC0415 + from scipy.spatial import cKDTree pts = self.curve m = pts.shape[0] @@ -1736,7 +1771,7 @@ def _raw_field( def _cell_box(self: Sweep) -> Shape: """Tube SDF: ``dist_to_curve(p) − radial_max``.""" - from .implicit_ops import from_field # noqa: PLC0415 + from .implicit_ops import from_field bounds: BoundsType = ( self._bounds @@ -1771,7 +1806,8 @@ def _create_grid( y: npt.NDArray[np.float64], z: npt.NDArray[np.float64], ) -> pv.StructuredGrid: - """Map a parametric (s, r, θ) grid to Cartesian via the curve frames. + """ + Map a parametric (s, r, θ) grid to Cartesian via the curve frames. ``x``, ``y``, ``z`` here are the parametric meshgrid coordinates from :meth:`Tpms._compute_tpms_field` (they live in (s, r, θ) @@ -1855,7 +1891,7 @@ def lower_skeletal(self: Infill) -> pv.PolyData: """Lower-skeletal part as a PolyData mesh.""" return self.generate_vtk(type_part="lower skeletal") - def __init__( # noqa: PLR0913 + def __init__( self: Infill, obj: pv.PolyData, surface_function: Field, @@ -1866,7 +1902,8 @@ def __init__( # noqa: PLR0913 resolution: int = 20, density: float | None = None, ) -> None: - r"""Initialize the Infill object. + r""" + Initialize the Infill object. :param obj: object in which the infill is generated. Normals must be oriented\ towards the outside of the object. Use the `flip_normals` method if\ @@ -1949,7 +1986,8 @@ def _create_grid( return grid def _density_envelope_volume(self: Infill) -> float: - """Density is measured against the input object's volume. + """ + Density is measured against the input object's volume. Uses the volume of the *original* input object rather than the re-oriented one stored on ``self.obj`` — see ``__init__``. @@ -1957,7 +1995,8 @@ def _density_envelope_volume(self: Infill) -> float: return self._obj_volume def _setup_frep_field(self: Infill) -> None: - """Cartesian gyroid field, but bounds set to the obj bbox. + """ + Cartesian gyroid field, but bounds set to the obj bbox. The plain :class:`Tpms` parent uses origin-centered bounds of size ``cell_size * repeat_cell``. For an Infill of an object that is *not* @@ -1986,14 +2025,15 @@ def _raw_field( self._finalize_frep(_raw_field, bounds) def _cell_box(self: Infill) -> Shape: - """Override the cell-box clip with the *object envelope* SDF. + """ + Override the cell-box clip with the *object envelope* SDF. Without this override the F-rep marching cubes path (used by :meth:`Tpms.generate_vtk` and :meth:`Tpms.generate`) would clip the TPMS to the cartesian bounding box rather than to the input object — producing volumes well above ``obj.volume`` and corrupting density. """ - from .implicit_ops import from_field # noqa: PLC0415 + from .implicit_ops import from_field bounds_arr = np.array(self.obj.bounds) bounds: BoundsType = tuple(bounds_arr.tolist()) @@ -2013,7 +2053,8 @@ def _envelope_sdf( class GradedInfill(Infill): - """Cartesian TPMS infill with offset graded by distance to the envelope. + """ + Cartesian TPMS infill with offset graded by distance to the envelope. Same as :class:`Infill` (TPMS in the world cartesian frame, clipped by the input object), except the *thickness* of the TPMS sheet varies with @@ -2029,7 +2070,7 @@ class GradedInfill(Infill): bunny-style envelopes with high-curvature features. """ - def __init__( # noqa: PLR0913 + def __init__( self: GradedInfill, obj: pv.PolyData, surface_function: Field, @@ -2042,7 +2083,8 @@ def __init__( # noqa: PLR0913 phase_shift: Sequence[float] = (0.0, 0.0, 0.0), resolution: int = 20, ) -> None: - """Build a graded TPMS infill. + """ + Build a graded TPMS infill. :param obj: envelope mesh :param surface_function: TPMS function ``f(x,y,z)`` @@ -2086,7 +2128,8 @@ def _make_graded_offset_callable( transition: float, smoothness: float, ) -> Field: - """Return ``f(x, y, z) -> offset(p)`` graded by SDF. + """ + Return ``f(x, y, z) -> offset(p)`` graded by SDF. ``d_norm = clip(-d / depth_max, 0, 1)`` ∈ [0, 1] (0 at skin, 1 in deepest interior point). Offset is interpolated from diff --git a/microgen/shape/tpms_grading.py b/microgen/shape/tpms_grading.py index 91cced18..5a50566f 100644 --- a/microgen/shape/tpms_grading.py +++ b/microgen/shape/tpms_grading.py @@ -22,7 +22,8 @@ def compute_offset( self: OffsetGrading, grid: pv.UnstructuredGrid | pv.StructuredGrid, ) -> npt.NDArray[np.float64]: - """Compute the offset of the grid. + """ + Compute the offset of the grid. This method should compute the offset on each point of the grid and return \ it as a 1D array. @@ -54,7 +55,8 @@ def __init__( furthest_offset: float, boundary_weight: float = 1.0, ) -> None: - """Initialize the ImplicitDistance object. + """ + Initialize the ImplicitDistance object. Parameters ---------- diff --git a/pyproject.toml b/pyproject.toml index d2db71a6..d32acc5a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -85,9 +85,47 @@ exclude = [ "microgen/rve.py", "microgen/single_mesh.py", "examples/*", + "gyroid_gd_bunny.py", ] # external = ["E131", "D102", "D105"] select = ["ALL"] -ignore = ["E501", "S101", "ANN101", "ANN102", "ANN003"] +ignore = [ + "E501", # line too long — ruff-format handles wrapping + "S101", # assert in tests is fine + "ANN101", # removed rule + "ANN102", # removed rule + "ANN001", # missing type hint on argument + "ANN002", # missing type hint on *args + "ANN003", # missing **kwargs annotation + "ANN202", # missing return-type on private function + "PLC0415", # late imports needed to avoid circular deps + "N806", # capitalised local variables (T, N, B for tangent/normal/binormal) + "PLR2004", # magic-value comparison — too noisy in scientific code + "PLR0913", # too many arguments — APIs are wide on purpose + "PLR0915", # too many statements + "LOG015", # root logger usage + "SLF001", # private member access between class siblings + "RUF001", # ambiguous unicode (we use Greek letters intentionally) + "RUF002", # same, in docstrings + "RUF003", # same, in comments + "RUF046", # int(...) on already-int — bookkeeping + "RUF007", # zip strict + "ARG002", # unused method argument (pyvista hook signatures) + "ARG005", # unused lambda argument (kwargs trampolines) + "TID252", # relative imports are the project convention + "FBT001", # boolean positional argument + "FBT002", # boolean default-value + "TC003", # typing-only imports — minor + "RET504", # unnecessary assignment before return + "PERF401", # list comprehension — readability tradeoff + "TRY003", # long messages outside exception class + "EM101", # raw string in exception + "EM102", # f-string in exception + "SIM108", # ternary instead of if/else + "UP038", # isinstance with X | Y — deprecated rule (perf regression) + # Docstring style — too noisy for scientific code + "D100", "D101", "D102", "D103", "D104", "D105", "D107", + "D203", "D205", "D211", "D212", "D213", "D400", "D401", "D415", +] [tool.ruff.format] docstring-code-format = true From d1b2104387c0c2ab7e716ed0c1fcf91d8fd8dde0 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Mon, 27 Apr 2026 11:36:11 +0200 Subject: [PATCH 11/11] Add OffsetGrading.as_field and use in Tpms Implement OffsetGrading.as_field to produce a (x,y,z)->array callable that builds transient PolyData and calls compute_offset, allowing grading to be re-evaluated on arbitrary grids (used by the F-rep / marching-cubes path). Store this callable on Tpms instances when given an OffsetGrading so shell()/as_sheet() can re-sample offsets on the marching-cubes grid (avoids failures like float(t) on array inputs). Also replace a hardcoded OUT path in the example with Path(__file__).parent and tidy related imports. --- examples/tpms_infill_gallery.py | 2 +- microgen/shape/tpms.py | 5 +++++ microgen/shape/tpms_grading.py | 37 +++++++++++++++++++++++++++++++-- 3 files changed, 41 insertions(+), 3 deletions(-) diff --git a/examples/tpms_infill_gallery.py b/examples/tpms_infill_gallery.py index b3aded0c..a3b53c24 100644 --- a/examples/tpms_infill_gallery.py +++ b/examples/tpms_infill_gallery.py @@ -40,7 +40,7 @@ CYLINDER_RADIUS = 1.5 CYLINDER_HEIGHT = 6.0 -OUT = Path("/Users/ychemisky/Documents/GitHub/microgen") +OUT = Path(__file__).parent def clip_y(mesh: pv.DataSet, frac: float = -0.2) -> pv.DataSet: diff --git a/microgen/shape/tpms.py b/microgen/shape/tpms.py index 9f98c1fc..3d0724ff 100644 --- a/microgen/shape/tpms.py +++ b/microgen/shape/tpms.py @@ -647,6 +647,11 @@ def offset( self._offset = offset elif isinstance(offset, OffsetGrading): self._offset = offset.compute_offset(self.grid) + # Also expose the grading as a callable so the F-rep path + # (``as_sheet`` → ``shell``) can re-evaluate the offset on the + # marching-cubes grid. Without this, ``shell()`` would receive + # the array sampled on ``self.grid`` and fail at ``float(t)``. + self._offset_func = offset.as_field() elif callable(offset): # Keep the callable so the F-rep path can re-evaluate it on the # marching-cubes grid (which differs from ``self.grid``). diff --git a/microgen/shape/tpms_grading.py b/microgen/shape/tpms_grading.py index 5a50566f..c1a4b847 100644 --- a/microgen/shape/tpms_grading.py +++ b/microgen/shape/tpms_grading.py @@ -3,20 +3,53 @@ from __future__ import annotations from abc import ABC, abstractmethod +from collections.abc import Callable from typing import TYPE_CHECKING +import numpy as np import pyvista as pv from typing_extensions import override if TYPE_CHECKING: - import numpy as np import numpy.typing as npt - import pyvista as pv + + Field = Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]], + npt.NDArray[np.float64], + ] class OffsetGrading(ABC): """Base class for offset grading functions.""" + def as_field(self: OffsetGrading) -> Field: + """Return a ``(x, y, z) -> array`` callable form of this grading. + + Used by the F-rep marching-cubes path to re-evaluate the offset on + an arbitrary grid (the one chosen by marching cubes is *not* the + same as the structured grid used by the legacy clip-scalar path). + + The default implementation builds a transient :class:`pv.PolyData` + of the queried points and forwards to :meth:`compute_offset`, which + works for any grading whose ``compute_offset`` only consults point + coordinates (it does for :class:`NormedDistance`). Subclasses that + depend on cell connectivity should override this method. + """ + + def _field( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + _self: OffsetGrading = self, + ) -> npt.NDArray[np.float64]: + shape = x.shape + pts = np.column_stack([x.ravel(), y.ravel(), z.ravel()]) + polydata = pv.PolyData(pts) + values = _self.compute_offset(polydata) + return np.asarray(values).reshape(shape) + + return _field + @abstractmethod def compute_offset( self: OffsetGrading,