"""Cell-map construction helpers for mapped metamaterials."""
# STD
import math
import numbers
# Third-party
import numpy as np
# OpenVCAD
import pyvcad as pv
# Local
from ._domain import as_vec3
from .fields import _coerce_attribute
def _as_cells(value):
values = list(value)
if len(values) != 3:
raise ValueError("cells must contain exactly three U/V/W counts")
cells = tuple(int(item) for item in values)
if any(item <= 0 for item in cells):
raise ValueError("cells must contain positive U/V/W counts")
if any(float(values[index]) != cells[index] for index in range(3)):
raise ValueError("cells must contain integer U/V/W counts")
return cells
def _is_cell_triplet(value):
try:
values = list(value)
except TypeError:
return False
return len(values) == 3 and all(isinstance(item, numbers.Real) for item in values)
def _broadcast_cells(cells, count):
if _is_cell_triplet(cells):
return [_as_cells(cells)] * count
values = list(cells)
if len(values) != count:
raise ValueError("per-face cells must contain one U/V/W triplet for every face")
return [_as_cells(value) for value in values]
def _as_maps(cell_map):
if isinstance(cell_map, pv.CellMap):
return [cell_map]
if isinstance(cell_map, pv.CellMapCollection):
return list(cell_map.maps)
raise TypeError("expected a pv.CellMap or pv.CellMapCollection")
[docs]
def rectangular_cell_map(bounds, cell_size=None, cells=None, linear=True):
"""Create an axis-aligned finite cell map.
Exactly one of ``cell_size`` and ``cells`` is required. Requested cell sizes
are converted to integer counts with ceiling so the supplied bounds remain
exact and opposite boundaries never drift.
"""
if (cell_size is None) == (cells is None):
raise ValueError("provide exactly one of cell_size or cells")
minimum, maximum = (as_vec3(value) for value in bounds)
if cells is None:
size = as_vec3(cell_size)
if size.x <= 0.0 or size.y <= 0.0 or size.z <= 0.0:
raise ValueError("cell_size must be positive on every axis")
cells = (
int(math.ceil((maximum.x - minimum.x) / size.x)),
int(math.ceil((maximum.y - minimum.y) / size.y)),
int(math.ceil((maximum.z - minimum.z) / size.z)),
)
return pv.CellMap.rectangular(minimum, maximum, _as_cells(cells), linear=linear)
def _positive_size(value, name):
try:
size = float(value)
except (TypeError, ValueError):
raise ValueError("{} must be a positive finite number".format(name)) from None
if not math.isfinite(size) or size <= 0.0:
raise ValueError("{} must be a positive finite number".format(name))
return size
def _positive_count(value, name):
if not isinstance(value, numbers.Real):
raise ValueError("{} must be a positive integer".format(name))
count = int(value)
if count <= 0 or float(value) != count:
raise ValueError("{} must be a positive integer".format(name))
return count
[docs]
def cylindrical_cell_map(
center,
inner_radius,
outer_radius,
height,
cells=None,
axis=pv.Vec3(0.0, 0.0, 1.0),
reference_direction=pv.Vec3(1.0, 0.0, 0.0),
start_angle_degrees=0.0,
sweep_degrees=360.0,
linear=False,
arc_count=None,
axial_cell_size=None,
radial_cell_size=None,
):
"""Create an analytic annular-cylinder cell map.
Provide either explicit ``cells=(circumferential, axial, radial)`` or the
complete sizing group ``arc_count``, ``axial_cell_size``, and
``radial_cell_size``. Target sizes are rounded up to counts so the requested
radii and height remain exact.
"""
sizing = (arc_count, axial_cell_size, radial_cell_size)
if cells is not None and any(value is not None for value in sizing):
raise ValueError(
"provide either cells or the complete cylindrical sizing group, not both"
)
if cells is None:
if not all(value is not None for value in sizing):
raise ValueError(
"provide either cells or the complete cylindrical sizing group: "
"arc_count, axial_cell_size, and radial_cell_size"
)
axial_size = _positive_size(axial_cell_size, "axial_cell_size")
radial_size = _positive_size(radial_cell_size, "radial_cell_size")
height_value = _positive_size(height, "height")
radial_thickness = float(outer_radius) - float(inner_radius)
if not math.isfinite(radial_thickness) or radial_thickness <= 0.0:
raise ValueError("outer_radius must be greater than inner_radius")
cells = (
_positive_count(arc_count, "arc_count"),
int(math.ceil(height_value / axial_size)),
int(math.ceil(radial_thickness / radial_size)),
)
return pv.CellMap.cylindrical(
as_vec3(center),
float(inner_radius),
float(outer_radius),
float(height),
_as_cells(cells),
axis=as_vec3(axis),
reference_direction=as_vec3(reference_direction),
start_angle_degrees=float(start_angle_degrees),
sweep_degrees=float(sweep_degrees),
linear=linear,
)
[docs]
def spherical_cell_map(
center,
inner_radius,
outer_radius,
latitude_bounds_degrees,
cells=None,
north_axis=pv.Vec3(0.0, 0.0, 1.0),
reference_direction=pv.Vec3(1.0, 0.0, 0.0),
start_longitude_degrees=0.0,
longitude_sweep_degrees=360.0,
linear=False,
longitude_count=None,
latitude_count=None,
radial_cell_size=None,
):
"""Create an analytic spherical-shell band cell map.
Provide either explicit ``cells=(longitude, latitude, radial)`` or the
complete sizing group ``longitude_count``, ``latitude_count``, and
``radial_cell_size``. Radial target size is rounded up to a count so the
requested radii remain exact.
"""
sizing = (longitude_count, latitude_count, radial_cell_size)
if cells is not None and any(value is not None for value in sizing):
raise ValueError(
"provide either cells or the complete spherical sizing group, not both"
)
if cells is None:
if not all(value is not None for value in sizing):
raise ValueError(
"provide either cells or the complete spherical sizing group: "
"longitude_count, latitude_count, and radial_cell_size"
)
radial_size = _positive_size(radial_cell_size, "radial_cell_size")
radial_thickness = float(outer_radius) - float(inner_radius)
if not math.isfinite(radial_thickness) or radial_thickness <= 0.0:
raise ValueError("outer_radius must be greater than inner_radius")
cells = (
_positive_count(longitude_count, "longitude_count"),
_positive_count(latitude_count, "latitude_count"),
int(math.ceil(radial_thickness / radial_size)),
)
latitude_bounds = tuple(float(value) for value in latitude_bounds_degrees)
if len(latitude_bounds) != 2:
raise ValueError("latitude_bounds_degrees must contain exactly two values")
return pv.CellMap.spherical(
as_vec3(center),
float(inner_radius),
float(outer_radius),
latitude_bounds,
_as_cells(cells),
north_axis=as_vec3(north_axis),
reference_direction=as_vec3(reference_direction),
start_longitude_degrees=float(start_longitude_degrees),
longitude_sweep_degrees=float(longitude_sweep_degrees),
linear=linear,
)
def _as_trim_mode(trim):
if isinstance(trim, pv.TrimMode):
return trim
if isinstance(trim, str):
value = getattr(pv.TrimMode, trim, None)
if isinstance(value, pv.TrimMode):
return value
raise ValueError("trim must be a pv.TrimMode or one of 'none', 'cells', 'boundary'")
def _surface_height(surface, height):
"""Coerce one surface's map height into a float or a ``pv.FloatAttribute``.
Accepts a number, a ``pv.FloatAttribute``, a numba callable, or a per-vertex array (list,
tuple, or NumPy array) for a ``pv.TriangleMeshSurface``.
"""
if isinstance(height, pv.FloatAttribute):
return height
if isinstance(height, numbers.Real):
return float(height)
if isinstance(height, (list, tuple, np.ndarray)):
if not isinstance(surface, pv.TriangleMeshSurface):
raise TypeError("a per-vertex height array requires a pv.TriangleMeshSurface")
values = np.asarray(height, dtype=float).ravel()
if values.size != surface.vertex_count:
raise ValueError(
f"per-vertex height has {values.size} values but the surface has "
f"{surface.vertex_count} vertices"
)
if not np.all(np.isfinite(values)):
raise ValueError("per-vertex height values must be finite")
return surface.vertex_field(values.tolist())
if callable(height):
return _coerce_attribute(height)
raise TypeError(
"height must be a number, a pv.FloatAttribute, a numba callable, or a per-vertex array"
)
[docs]
def cell_map_from_surface(surface_or_surfaces, cells, height, linear=False, trim=pv.TrimMode.boundary, standoff=0.0):
"""Create one signed-height map per parametric surface.
Accepts a single ``ParametricSurface`` (a CADFace or TriangleMeshSurface), or a
non-empty iterable that mixes any surface types. A single surface returns a
``pv.CellMap``; an iterable returns a ``pv.CellMapCollection`` with per-surface
cell counts broadcast from ``cells``.
``height`` is the signed thickness of the map in mm along the surface normal. It may be:
- a number, for a constant height (the original behaviour);
- a ``pv.FloatAttribute`` (or a numba callable), sampled once per cell-map column at the
column's base point on the surface, for a variable height;
- a per-vertex array (list or NumPy array with one value per vertex) for a
``pv.TriangleMeshSurface``, interpolated across the patch. ``surface_gap`` computes one
by ray casting to a target such as the inside of a shell, so the cells end on it.
A variable height must stay finite, non-zero, and of one sign over the whole map. With
several surfaces, ``height`` is either broadcast to all of them or given as a list with one
entry per surface.
``trim`` selects how the map is restricted to the surface's true domain: a
``pv.TrimMode`` or one of the strings ``"none"``, ``"cells"``, ``"boundary"``.
``standoff`` is a signed stand-off distance along the surface normal at which
the map begins, so it spans ``[standoff, standoff + height]`` off the surface.
"""
trim = _as_trim_mode(trim)
if isinstance(surface_or_surfaces, pv.ParametricSurface):
return pv.CellMap.from_surface(
surface_or_surfaces,
_as_cells(cells),
_surface_height(surface_or_surfaces, height),
linear=linear,
trim=trim,
standoff=float(standoff),
)
surfaces = list(surface_or_surfaces)
if not surfaces or not all(isinstance(item, pv.ParametricSurface) for item in surfaces):
raise TypeError(
"surface_or_surfaces must be a ParametricSurface or a non-empty iterable of them"
)
counts = _broadcast_cells(cells, len(surfaces))
if isinstance(height, (list, tuple)):
if len(height) != len(surfaces):
raise ValueError("per-surface height must contain one entry for every surface")
heights = list(height)
else:
heights = [height] * len(surfaces)
maps = [
pv.CellMap.from_surface(
surface, count, _surface_height(surface, item), linear=linear, trim=trim, standoff=float(standoff)
)
for surface, count, item in zip(surfaces, counts, heights)
]
return pv.CellMapCollection(maps)
[docs]
def surface_gap(surface, target, max_distance, offset=0.0, missing="raise", tolerance=1e-3, voxel_size=0.25):
"""Measure the gap from every vertex of a mesh patch along its normal to a target.
This is the head-to-shell gap that drives an exact-fit variable-height map: with
``gap = mm.surface_gap(head, shell, 60.0)``, ``mm.cell_map_from_surface(head, cells,
height=gap)`` grows cells from the head patch that end on the shell.
Arguments:
surface (pv.TriangleMeshSurface): Patch whose vertices and vertex normals cast the rays.
target: What the rays hit. A ``pv.TriangleMeshSurface`` or ``pv.SurfaceMesh`` is hit by
exact ray-triangle intersection, so open or non-watertight meshes (a shell's inner
surface exported from CAD) work. Any other ``pyvcad`` node is sphere-traced against its
signed distance field (a private clone is prepared with ``voxel_size``).
max_distance (float): Longest ray in mm.
offset (float): Subtracted from every distance, for example a stand-off or a clearance.
missing: What to do for rays that hit nothing within ``max_distance``: ``"raise"``
(default) raises ``ValueError``, ``"nearest"`` copies the value of the nearest vertex
that did hit, and a number is used as is.
tolerance (float): Hit tolerance in mm for sphere tracing.
voxel_size (float): Voxel size in mm used to prepare a node target.
Returns:
numpy.ndarray: One distance in mm per surface vertex, minus ``offset``.
"""
if not isinstance(surface, pv.TriangleMeshSurface):
raise TypeError("surface_gap requires a pv.TriangleMeshSurface")
distances = np.asarray(
pv.ray_distances(
surface.vertices,
surface.vertex_normals,
target,
float(max_distance),
tolerance=float(tolerance),
voxel_size=float(voxel_size),
),
dtype=float,
)
misses = ~np.isfinite(distances)
if np.any(misses):
if isinstance(missing, numbers.Real):
distances[misses] = float(missing)
elif missing == "nearest":
if np.all(misses):
raise ValueError("surface_gap: no ray hit the target within max_distance")
points = np.array([[p.x, p.y, p.z] for p in surface.vertices], dtype=float)
hit_points = points[~misses]
hit_values = distances[~misses]
for index in np.flatnonzero(misses):
nearest = np.argmin(np.sum((hit_points - points[index]) ** 2, axis=1))
distances[index] = hit_values[nearest]
elif missing == "raise":
raise ValueError(
f"surface_gap: {int(np.count_nonzero(misses))} of {distances.size} rays hit nothing "
"within max_distance; increase max_distance or pass missing='nearest' or a number"
)
else:
raise ValueError("missing must be 'raise', 'nearest', or a number")
return distances - float(offset)
[docs]
def cell_map_from_cad_face(face_or_faces, cells, height, linear=False, trim=pv.TrimMode.boundary, standoff=0.0):
"""Create one native-UV signed-height map per CAD face."""
if isinstance(face_or_faces, pv.CADFace):
faces_ok = True
else:
try:
faces = list(face_or_faces)
except TypeError:
faces_ok = False
else:
faces_ok = bool(faces) and all(isinstance(face, pv.CADFace) for face in faces)
if not faces_ok:
raise TypeError("face_or_faces must be a CADFace or a non-empty iterable of CADFace objects")
return cell_map_from_surface(face_or_faces, cells, height, linear=linear, trim=trim, standoff=standoff)
[docs]
def cell_map_between_surfaces(
first_surface,
second_surface,
cells,
linear=False,
trim=pv.TrimMode.boundary,
auto_align=True,
first_origin=None,
flip_first=False,
second_origin=None,
flip_second=False,
standoff=0.0,
):
"""Create a curvilinear volume between two parametric surfaces of any type.
``first_surface`` and ``second_surface`` may be any mix of CADFace and
TriangleMeshSurface. In advanced mode, ``flip_first`` and ``flip_second``
request opposing UV orientation changes. Origins are normalized ``(u, v)``
seam shifts and non-zero values are valid only for periodic directions.
``standoff`` shifts the first surface along its oriented normal (in mm) before the two
surfaces are interpolated, so the volume starts a clearance away from it. Zero (the
default) leaves the map unchanged; the second surface never moves.
"""
if not isinstance(first_surface, pv.ParametricSurface) or not isinstance(
second_surface, pv.ParametricSurface
):
raise TypeError("first_surface and second_surface must be pv.ParametricSurface objects")
relative_flip = bool(flip_first) != bool(flip_second)
return pv.CellMap.between_surfaces(
first_surface,
second_surface,
_as_cells(cells),
linear=linear,
trim=_as_trim_mode(trim),
auto_align=auto_align,
flip_u=relative_flip,
flip_v=relative_flip,
exchange_uv=False,
first_origin=first_origin,
second_origin=second_origin,
standoff=float(standoff),
)
[docs]
def cell_map_between_cad_faces(
first_face,
second_face,
cells,
linear=False,
trim=pv.TrimMode.boundary,
auto_align=True,
first_origin=None,
flip_first=False,
second_origin=None,
flip_second=False,
standoff=0.0,
):
"""Create a curvilinear volume between two compatible CAD faces.
In advanced mode, ``flip_first`` and ``flip_second`` request opposing UV
orientation changes. Origins are normalized ``(u, v)`` seam shifts and
non-zero values are valid only for periodic CAD directions. ``standoff``
shifts the first face along its oriented normal before interpolating.
"""
if not isinstance(first_face, pv.CADFace) or not isinstance(second_face, pv.CADFace):
raise TypeError("first_face and second_face must be pv.CADFace objects")
return cell_map_between_surfaces(
first_face,
second_face,
cells,
linear=linear,
trim=trim,
auto_align=auto_align,
first_origin=first_origin,
flip_first=flip_first,
second_origin=second_origin,
flip_second=flip_second,
standoff=standoff,
)
[docs]
def rim(cell_map, width=1.0):
"""Build a solid rim along the edges of each surface map.
The rim is a signed distance field from the map's edge wall (the surface's
cropped boundary swept through the map height for boundary-trimmed maps, and
the exposed lateral cell faces otherwise). It extends ``width`` inward from the edge,
so unioning it with a lattice built on the same map seals the lattice's open edge cells.
``width`` may be a constant or a ``pv.FloatAttribute`` sampled through the map.
"""
maps = _as_maps(cell_map)
nodes = [pv.MappedRim(item, _coerce_attribute(width)) for item in maps]
return nodes[0] if isinstance(cell_map, pv.CellMap) else pv.MappedStructureCollection(maps, nodes)
[docs]
def domain(cell_map):
"""Return the active volume of a cell map as a solid node (``pv.CellMapDomain``).
The node's field is the map's signed domain distance in mm: negative inside the active,
trimmed map volume and positive outside, including trimmed margins and holes. Use it to clip
another node to a liner volume (``pv.Intersection``), as a keep-out (``pv.Difference``), or as
the host of a stochastic lattice. A ``pv.CellMapCollection`` returns the union of its maps'
domains.
"""
maps = _as_maps(cell_map)
nodes = [pv.CellMapDomain(item) for item in maps]
if isinstance(cell_map, pv.CellMap):
return nodes[0]
return pv.Union(0.0, nodes)
[docs]
def map_node(cell_map, node, space="physical", scale=None, clip_to_domain=True):
"""Wrap any node through a cell map so geometry designed flat follows the map.
The node is authored in a flat child box and evaluated at the map's logical coordinates
scaled to child units. With ``space="physical"`` (default) the scale is the map's average
millimetres per cell along U, V, and W, so the child box ``pv.MappedNode.child_bounds``
covers the whole map in approximate millimetres; ``space="logical"`` uses cell units, and
``scale`` (a ``pv.Vec3`` or 3-tuple) sets explicit child units per cell. Typical uses are pads,
bosses, labels, or cut-outs placed on a curved liner. Distances are first-order Jacobian
estimates: accurate near the surface (safe for meshing) but approximate far from it.
By default the wrapped node ends where the map ends (its field is clipped with the map's
domain distance). Pass ``clip_to_domain=False`` when the node marks a zone of the map, for
example a region for :func:`zones` or :func:`overlap`: inside the domain it then keeps its own
depth, so growing or shrinking the region by a band width works as expected.
"""
if not isinstance(cell_map, pv.CellMap):
raise TypeError("map_node requires a pv.CellMap")
if scale is not None:
return pv.MappedNode(cell_map, node, space=space, scale=as_vec3(scale),
clip_to_domain=clip_to_domain)
return pv.MappedNode(cell_map, node, space=space, clip_to_domain=clip_to_domain)
[docs]
def warp_cell_map(
cell_map,
u_scale=0.0,
v_scale=0.0,
w_scale=0.0,
lock_u_bounds=True,
lock_v_bounds=True,
lock_w_bounds=True,
):
"""Warp relative U/V/W spacing with bounded positive exponential weights.
Each cell's share of the map along an axis is weighted by
``exp(clamp(scale, -4, 4))``, sampled on the unwarped base map. A scale of
``0`` keeps the original spacing; larger values give that region more room.
The base map is not modified.
Args:
cell_map (pv.CellMap | pv.CellMapCollection): Base map or maps.
u_scale, v_scale, w_scale (float | pv.FloatAttribute | callable):
Relative spacing controls per axis, default ``0.0``. These are
weights, not cell sizes in mm. ``cell_map.logical_position`` is a
convenient input.
lock_u_bounds, lock_v_bounds, lock_w_bounds (bool): Keep the two outer
faces on that axis fixed and redistribute the cells between them
(default ``True``).
Returns:
pv.CellMap | pv.CellMapCollection: The warped map, same kind as the input.
Raises:
RuntimeError: If a map was built with ``TrimMode.boundary`` (rebuild it
with ``trim="none"`` before warping), or if the warp would fold
the domain.
Example:
>>> warped = mm.warp_cell_map(cell_map, u_scale=cell_map.logical_position.x * 0.2)
"""
maps = _as_maps(cell_map)
warped = [
pv.CellMap.warp(
item,
_coerce_attribute(u_scale),
_coerce_attribute(v_scale),
_coerce_attribute(w_scale),
lock_u_bounds=lock_u_bounds,
lock_v_bounds=lock_v_bounds,
lock_w_bounds=lock_w_bounds,
)
for item in maps
]
return warped[0] if isinstance(cell_map, pv.CellMap) else pv.CellMapCollection(warped)
