"""Beam lattices built from the cells of a mesh, with the tessellation rules.
The rules of :data:`GRAPH_RULES` (``tessellations.py``) apply to any cells: here they run on
the tetrahedra or hexahedra of a volume mesh, or the triangles of a surface mesh, and the
struts become a ``pv.BeamLattice``. An adaptive tetrahedral mesh gives a graded cell size for
free.
"""
# STD
import math
import numbers
# Third-party
import numpy as np
# OpenVCAD
import pyvcad as pv
# Local
from .fields import _coerce_attribute
from .tessellations import Polyhedron, _rule_names, _rule_segments
_TRIANGLE_FACETS = ((0, 1), (1, 2), (2, 0))
_TET_FACES = ((0, 2, 1), (0, 1, 3), (0, 3, 2), (1, 2, 3))
# VTK node order, as HexahedralMesh writes it.
_HEX_FACES = ((0, 3, 2, 1), (4, 5, 6, 7), (0, 1, 5, 4), (1, 2, 6, 5), (2, 3, 7, 6), (3, 0, 4, 7))
_FACES = {3: ("triangle", _TRIANGLE_FACETS), 4: ("tet", _TET_FACES), 8: ("hex", _HEX_FACES)}
def _points(values):
return np.asarray([
(value.x, value.y, value.z) if isinstance(value, pv.Vec3) else tuple(value)
for value in values
], dtype=np.float64).reshape(-1, 3)
def _mesh_cells(mesh):
"""Return the mesh's vertex positions and its cells as vertex-index tuples."""
if isinstance(mesh, (pv.TetrahedralMesh, pv.HexahedralMesh)):
return _points(mesh.nodes), [tuple(int(index) for index in cell) for cell in mesh.elements]
if isinstance(mesh, pv.SurfaceMesh):
return _points(mesh.vertices), [tuple(int(index) for index in cell) for cell in mesh.triangles]
if isinstance(mesh, (tuple, list)) and len(mesh) == 2:
vertices, cells = mesh
return np.asarray(vertices, dtype=np.float64).reshape(-1, 3), \
[tuple(int(index) for index in cell) for cell in cells]
raise TypeError("mesh must be a pv.TetrahedralMesh, pv.HexahedralMesh, pv.SurfaceMesh or "
"(vertices, cells)")
def _complex(points, cells):
"""Build the rule engine's cells and their face adjacency (shared vertex sets)."""
polyhedra = []
for ids in cells:
if len(ids) not in _FACES:
raise ValueError("cells must have 3 (triangle), 4 (tetrahedron) or 8 (hexahedron) vertices")
if min(ids) < 0 or max(ids) >= len(points):
raise ValueError("cell vertex index out of range")
kind, faces = _FACES[len(ids)]
polyhedra.append(Polyhedron._from_parts(points[list(ids)], faces, kind=kind, vertex_ids=ids))
owners = {}
for cell_index, cell in enumerate(polyhedra):
for face_index, face in enumerate(cell.faces):
key = frozenset(cell.vertex_ids[index] for index in face)
owners.setdefault(key, []).append((cell_index, face_index))
zero = np.zeros(3)
adjacency = [
(first[0], first[1], second[0], second[1], zero)
for entries in owners.values() if len(entries) == 2
for first, second in [entries]
]
return polyhedra, adjacency, owners
def _radius_source(value, name):
"""A number, a per-vertex array, or a field; returns (kind, data)."""
if isinstance(value, numbers.Real):
if not math.isfinite(float(value)) or float(value) <= 0.0:
raise ValueError(name + " must be positive")
return "constant", float(value)
if isinstance(value, (list, tuple, np.ndarray)):
return "vertices", np.asarray(value, dtype=np.float64).reshape(-1)
return "field", _coerce_attribute(value)
def _vertex_values(source, name, lattice_points, mesh_points, cells, owners, tolerance):
"""Evaluate a radius source at every lattice vertex."""
kind, data = source
if kind == "constant":
return np.full(len(lattice_points), data)
if kind == "field":
values = np.asarray([
float(data.sample(float(x), float(y), float(z), 0.0)) for x, y, z in lattice_points
])
else:
if len(data) != len(mesh_points):
raise ValueError("{} has {} values; the mesh has {} vertices".format(
name, len(data), len(mesh_points)))
# Mesh vertices take their own value; cell and face centres the mean of their vertices.
known = {}
for ids in [tuple(sorted(key)) for key in owners] + cells:
centre = mesh_points[list(ids)].mean(axis=0)
known.setdefault(tuple(np.round(centre / tolerance).astype(np.int64)), float(data[list(ids)].mean()))
for index, point in enumerate(mesh_points):
known[tuple(np.round(point / tolerance).astype(np.int64))] = float(data[index])
values = []
for point in lattice_points:
key = tuple(np.round(point / tolerance).astype(np.int64))
if key not in known:
raise ValueError("could not find the mesh vertices behind a lattice point")
values.append(known[key])
values = np.asarray(values)
if not np.all(np.isfinite(values)) or np.any(values <= 0.0):
raise ValueError(name + " must be positive and finite")
return values
[docs]
def mesh_lattice(mesh, rules="edges", beam_radius=0.5, node_radius=None, min_length=0.0):
"""Build a beam lattice from the cells of a volume or surface mesh.
The tessellation rules of :data:`GRAPH_RULES` run on the mesh cells: ``"edges"`` keeps
every mesh edge, ``"cell_vertex_centroid"`` joins each cell's centre to its corners,
``"facet_vertex_centroid"`` joins each face's centre to its corners, ``"dual"`` joins the
centres of cells that share a face (on a triangle surface, a strut across every interior
edge), and ``"central_axes"`` joins each cell's centre to its face centres. Struts that
neighbouring cells share are kept once. A list of rules gives their union.
Struts are centred on the mesh, so they bulge out of the part by their radius where they
run along its surface; intersect the result with the part to trim them.
Args:
mesh (pv.TetrahedralMesh | pv.HexahedralMesh | pv.SurfaceMesh | tuple): The cells. A
``(vertices, cells)`` tuple takes an ``(n, 3)`` array and cells of 3 (triangle),
4 (tetrahedron) or 8 (hexahedron, VTK order) vertex indices.
rules (str | list[str]): One or more of :data:`GRAPH_RULES`, default ``"edges"``.
beam_radius (float | array-like | pv.FloatAttribute | callable): Strut radius in mm: a
number, one value per mesh vertex (cell and face centres take the mean of their
vertices), or a field sampled at both ends of every strut.
node_radius (float | array-like | pv.FloatAttribute | callable | None): Radius in mm of
a ball at every lattice vertex, given the same ways; ``None`` (default) adds none.
min_length (float): Struts shorter than this (mm) are dropped, default ``0``.
Returns:
pv.BeamLattice: The lattice (new geometry; its distance is exact outside the struts and
a bound inside, like any union of struts).
Example:
>>> settings = pv.TetrahedralMeshFixedSettings()
>>> settings.cell_size = 6.0
>>> tets = pv.TetrahedralMesh.from_tree(pv.Sphere(pv.Vec3(0, 0, 0), 20.0), settings)
>>> lattice = mm.mesh_lattice(tets, "edges", beam_radius=0.6)
"""
names = _rule_names(rules)
points, cells = _mesh_cells(mesh)
if not cells:
raise ValueError("the mesh has no cells")
polyhedra, adjacency, owners = _complex(points, cells)
segments = []
for name in names:
segments.extend(_rule_segments(polyhedra, name, None, None, adjacency))
if not segments:
raise ValueError("the rules made no struts on this mesh")
size = float(np.linalg.norm(points.max(axis=0) - points.min(axis=0)))
tolerance = max(1e-12, 1e-9 * size)
ends = np.asarray(segments, dtype=np.float64).reshape(-1, 3)
keys = np.round(ends / tolerance).astype(np.int64)
_, first, inverse = np.unique(keys, axis=0, return_index=True, return_inverse=True)
lattice_points = ends[first]
pairs = np.sort(inverse.reshape(-1, 2), axis=1)
pairs = np.unique(pairs[pairs[:, 0] != pairs[:, 1]], axis=0)
radii = _vertex_values(_radius_source(beam_radius, "beam_radius"), "beam_radius",
lattice_points, points, cells, owners, tolerance)
vertices = [pv.Vec3(float(x), float(y), float(z)) for x, y, z in lattice_points]
beams = [pv.Beam(int(a), int(b), r1=float(radii[a]), r2=float(radii[b])) for a, b in pairs]
if node_radius is None:
return pv.BeamLattice(vertices, beams, float(radii.max()), min_length=float(min_length))
joints = _vertex_values(_radius_source(node_radius, "node_radius"), "node_radius",
lattice_points, points, cells, owners, tolerance)
balls = [pv.Ball(index, radius=float(joints[index])) for index in range(len(vertices))]
return pv.BeamLattice(vertices, beams, float(radii.max()), min_length=float(min_length),
ball_mode=pv.BallMode.Mixed, ball_radius=float(joints.max()), balls=balls)
