diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 00000000..5008ddfc Binary files /dev/null and b/.DS_Store differ 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/conda.recipe/meta.yaml b/conda.recipe/meta.yaml index bd895dea..26295ea4 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 0d9ed25c..2c1171cc 100755 --- a/environment.yml +++ b/environment.yml @@ -11,5 +11,7 @@ dependencies: - python-gmsh - meshio - cadquery + - pip: + - autograd - mmg # - neper diff --git a/examples/tpms_infill_gallery.py b/examples/tpms_infill_gallery.py new file mode 100644 index 00000000..a3b53c24 --- /dev/null +++ b/examples/tpms_infill_gallery.py @@ -0,0 +1,117 @@ +"""TPMS coordinate-frame gallery — sphere, cylinder, sweep. + +Three reference geometries to show off the parametric-grid TPMS classes, +each clipped by a y-axis plane so the interior is visible: + +| File | Class | What it shows | +| -------------------------- | -------------------- | ---------------------------------------------- | +| ``sphere_g4_tpms.vtk`` | ``SphericalTpms`` | gyroid wrapping a sphere of radius ``R`` | +| ``cylinder_g5_tpms.vtk`` | ``CylindricalTpms`` | gyroid wrapping a cylinder | +| ``sweep_g7_helix.vtk`` | ``Sweep`` | gyroid following an arbitrary helical curve | + +Knobs: + +- ``CELL_SIZE`` — uniform cubic cell edge (scalar → no stretched cells) +- ``OFFSET`` — sheet thickness +- ``SPHERE_RADIUS`` — radius of the spherical demo +- ``CYLINDER_RADIUS``, ``CYLINDER_HEIGHT`` — for the cylindrical demo + +Run with ``python examples/tpms_infill_gallery.py``. +""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pyvista as pv + +from microgen.shape.surface_functions import gyroid +from microgen.shape.tpms import CylindricalTpms, SphericalTpms, Sweep + + +# ----------------------------------------------------------------------------- +# Knobs +# ----------------------------------------------------------------------------- + +CELL_SIZE = 1.0 +OFFSET = 0.5 +SPHERE_RADIUS = 3.0 +CYLINDER_RADIUS = 1.5 +CYLINDER_HEIGHT = 6.0 + +OUT = Path(__file__).parent + + +def clip_y(mesh: pv.DataSet, frac: float = -0.2) -> 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..d18ee3dd --- /dev/null +++ b/gyroid_gd_bunny.py @@ -0,0 +1,82 @@ +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.0 + +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) + +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.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.0 + 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.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.show() 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 a269c08a..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. @@ -224,10 +226,24 @@ 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): + 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, @@ -239,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 @@ -312,6 +329,38 @@ 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,13 +378,18 @@ 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 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 @@ -351,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) @@ -374,7 +428,14 @@ 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) @@ -386,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 0e3bf302..6d9c95b2 100644 --- a/microgen/shape/surface_functions.py +++ b/microgen/shape/surface_functions.py @@ -1,11 +1,12 @@ """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: - """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 d238a4fd..3d0724ff 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 @@ -40,8 +42,56 @@ _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. + """ + Triply Periodical Minimal Surfaces. Class to generate Triply Periodical Minimal Surfaces (TPMS) geometry from a given mathematical function, with given offset @@ -65,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, @@ -76,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. @@ -110,6 +161,17 @@ 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 @@ -123,6 +185,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"]) @@ -157,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) @@ -166,21 +230,47 @@ 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. + + 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"], 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. + + 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}" 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] @@ -189,22 +279,37 @@ def _compute_offset_to_fit_density( ) return self.offset - temp_tpms = Tpms( - surface_function=self.surface_function, - offset=0.0, - resolution=resolution if resolution is not None else self.resolution, - ) + # 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 - grid_part = getattr(temp_tpms, f"grid_{part_type.replace(' ', '_')}") - return abs(grid_part.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) - target_density, + bracket=bracket, + ).root + except ValueError: + # Bracket sign mismatch — usually because target density exceeds + # 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 - 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 @@ -337,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, ) ] @@ -354,6 +460,173 @@ def _compute_tpms_field(self: Tpms) -> None: self.grid["surface"] = tpms_field.ravel(order="F") + 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 + + 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 + + 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]), + ) + + half = 0.5 * self.cell_size * self.repeat_cell + self._finalize_frep( + _raw_field, + (-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 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 + from .shape import Shape + + 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): + 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: + """ + 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 + + f = self._func + h = self._half_offset_field() + if callable(h): + + 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): + 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 + + f = self._func + h = self._half_offset_field() + if callable(h): + + 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): + 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 + + 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 + + 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 + + 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 + + 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 + + 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 @@ -368,12 +641,31 @@ 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) + # 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): - self._offset = offset(self.grid.x, self.grid.y, self.grid.z).ravel("F") + # 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 + # 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) @@ -381,9 +673,21 @@ 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 + 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]] @@ -394,149 +698,20 @@ def _create_shell(self: Tpms, mesh: pv.PolyData) -> cq.Shell: 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)) - 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 - - :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), - ] + if not faces: + err_msg = "Mesh has no triangles to convert into a shell" + raise ShellCreationError(err_msg) - 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, @@ -573,6 +748,67 @@ 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 _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", @@ -580,18 +816,22 @@ 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 """ - if type_part not in ["sheet", "lower skeletal", "upper skeletal", "surface"]: + 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 self._VALID_PARTS: err_msg = ( f"type_part ({type_part}) must be 'sheet', 'lower skeletal', " "'upper skeletal' or 'surface'" @@ -603,37 +843,98 @@ 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 + from OCP.TopExp import TopExp_Explorer + from OCP.TopoDS import TopoDS + + 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. - :param type_part: part of the TPMS desireds - :param algo_resolution: if offset must be computed to fit density, \ - resolution of the temporary TPMS used to compute the offset + 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. - :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 @@ -643,12 +944,56 @@ def generate_vtk( "'upper skeletal' or 'surface'" ) raise ValueError(err_msg) + + 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, + ) + return envelope_mesh.translate(xyz=self.center) + 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(" ", "_")) + # Note: no `_check_offset` here — the F-rep VTK path handles negative + # / zero / variable offsets gracefully. ``generate()`` still applies + # the historical CAD-side restriction. + + # 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) @@ -681,7 +1026,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, @@ -691,7 +1036,16 @@ def generateVtk( # noqa: N802 class CylindricalTpms(Tpms): """Class used to generate cylindrical TPMS geometries (sheet or skeletals parts).""" - def __init__( # noqa: PLR0913 + # 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__( self: CylindricalTpms, radius: float, surface_function: Field, @@ -704,7 +1058,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)$. @@ -762,9 +1117,17 @@ 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_[ @@ -774,11 +1137,133 @@ 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 + 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]), + ) + + # 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_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 + + 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).""" - def __init__( # noqa: PLR0913 + # 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__( self: SphericalTpms, radius: float, surface_function: Field, @@ -791,7 +1276,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)$. @@ -855,10 +1341,20 @@ 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), @@ -872,11 +1368,535 @@ def _create_grid( ] return grid + def _setup_frep_field(self: SphericalTpms) -> None: + """Build F-rep field in Cartesian coordinates (inverse spherical mapping).""" + 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]), + ) + + # Tightened bounds: only ``r_outer`` cube, not ``r_max + R``. + _, r_outer = self._shell_radii() + self._finalize_frep( + _raw_field, + (-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 + + 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__( + 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 + + 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 + + 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.""" - def __init__( # noqa: PLR0913 + _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__( self: Infill, obj: pv.PolyData, surface_function: Field, @@ -887,7 +1907,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\ @@ -906,8 +1927,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] @@ -957,8 +1990,203 @@ 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. + + Uses the volume of the *original* input object rather than the + re-oriented one stored on ``self.obj`` — see ``__init__``. + """ + return self._obj_volume -# Re-export for backward compatibility -from .shape import ShellCreationError # noqa: E402 + 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) -__all__ = ["CylindricalTpms", "Infill", "ShellCreationError", "SphericalTpms", "Tpms"] + 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 + + 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 GradedInfill(Infill): + """ + 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 + the signed distance to the envelope. + + 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. + + 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__( + self: GradedInfill, + obj: pv.PolyData, + surface_function: Field, + 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, + ) -> None: + """ + 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 + """ + self._gradation_params = ( + float(offset_skin), + float(offset_core), + float(transition), + float(smoothness), + ) + + graded_offset = self._make_graded_offset_callable(obj, *self._gradation_params) + + super().__init__( + obj=obj, + surface_function=surface_function, + offset=graded_offset, + cell_size=cell_size, + repeat_cell=repeat_cell, + phase_shift=phase_shift, + resolution=resolution, + density=None, + ) + + @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``. + """ + envelope_oriented = envelope.compute_normals( + auto_orient_normals=True, + point_normals=True, + cell_normals=True, + ) + + 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 _graded_offset( + x: npt.NDArray[np.float64], + y: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + ) -> npt.NDArray[np.float64]: + 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 + + +# Re-export for backward compatibility +from .shape import BoundsType, ShellCreationError # noqa: E402 + +__all__ = [ + "CylindricalTpms", + "GradedInfill", + "Infill", + "ShellCreationError", + "SphericalTpms", + "Sweep", + "Tpms", +] diff --git a/microgen/shape/tpms_grading.py b/microgen/shape/tpms_grading.py index 91cced18..c1a4b847 100644 --- a/microgen/shape/tpms_grading.py +++ b/microgen/shape/tpms_grading.py @@ -3,26 +3,60 @@ 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, 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 +88,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 f5c2b01d..d32acc5a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,7 +18,7 @@ classifiers = [ "Programming Language :: Python :: 3.12", # Python 3.13 not yet supported (CadQuery does not support 3.13) ] -dependencies = ["numpy", "pyvista", "gmsh", "meshio", "cadquery", "scipy", "nlopt"] +dependencies = ["numpy", "pyvista", "gmsh", "meshio", "cadquery", "scipy", "nlopt", "autograd"] [project.optional-dependencies] dev = [ @@ -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 diff --git a/tests/shapes/test_tpms.py b/tests/shapes/test_tpms.py index 0d963fce..7843a0af 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,17 @@ 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 ] @@ -689,3 +695,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}" 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))