"""Periodic tessellations and the rules that turn their cells into strut graphs.
A tessellation fills space with polyhedral cells. A rule turns the cells into struts:
============================= ===========================================================
Rule Struts
============================= ===========================================================
``"edges"`` every edge of every cell
``"cell_vertex_centroid"`` each cell's centre joined to each of its vertices
``"facet_vertex_centroid"`` each face's centre joined to each vertex of that face
``"dual"`` the centres of every two cells that share a face
``"central_axes"`` each cell's centre joined to the centre of each of its faces
============================= ===========================================================
:class:`PeriodicTessellation` holds cells that repeat with period one along U, V and W, in
normalized cell coordinates. :meth:`PeriodicTessellation.edges` returns a graph unit cell in the
catalog's edge-list format (origin-centred, spanning ``[-0.5, 0.5]``), ready for
:func:`graph_lattice`; :meth:`PeriodicTessellation.face_cell` returns the cell walls as a
``pv.FaceUnitCell`` for :func:`face_lattice`. :func:`periodic_voronoi` builds a tessellation from
seed points and :func:`tessellation` returns the built-in ones the catalog is made from. The same
rules turn the cells of a mesh into a beam lattice (:func:`mesh_lattice`).
"""
# STD
import itertools
import math
import numbers
# Third-party
import numpy as np
# OpenVCAD
import pyvcad as pv
GRAPH_RULES = (
"edges",
"cell_vertex_centroid",
"facet_vertex_centroid",
"dual",
"central_axes",
)
# Normalized coordinates are compared and returned on a grid of 10^-_DIGITS cell widths.
_DIGITS = 10
_SNAP = 1e-8
_OFFSETS = np.array(list(itertools.product((-1, 0, 1), repeat=3)), dtype=np.float64)
[docs]
class Polyhedron:
"""One cell of a tessellation: vertices, polygonal faces, a centre and optional tags.
Args:
vertices (array-like): ``(n, 3)`` vertex positions. In a
:class:`PeriodicTessellation` they are normalized cell coordinates (one repeat spans
``[0, 1]`` on each axis, and a cell may stick out of it).
faces (list[list[int]] | None): Faces as vertex-index polygons in boundary order. ``None``
(default) uses the convex hull of ``vertices``, with coplanar triangles merged into
one polygon; vertices that are not hull corners are dropped.
kind (str | None): Cell type, for selecting cells with ``kinds=...``.
face_tags (list | None): One tag per face, for selecting faces with ``face_tags=...``.
center (array-like | None): The cell centre used by the centroid, dual and axis rules.
Defaults to the average of the vertices.
Example:
>>> cube = mm.Polyhedron([(x, y, z) for x in (0, 1) for y in (0, 1) for z in (0, 1)])
>>> len(cube.faces)
6
"""
def __init__(self, vertices, faces=None, kind=None, face_tags=None, center=None):
points = np.asarray(vertices, dtype=np.float64).reshape(-1, 3)
if not np.all(np.isfinite(points)):
raise ValueError("vertices must be finite")
if faces is None:
if len(points) < 4:
raise ValueError("a convex hull needs at least four vertices")
points, faces = _hull_faces(points)
faces = [tuple(int(index) for index in face) for face in faces]
for face in faces:
if len(face) < 3 or len(set(face)) != len(face):
raise ValueError("every face needs at least three distinct vertices")
if min(face) < 0 or max(face) >= len(points):
raise ValueError("face index out of range")
self._set(points, faces, kind, face_tags, center)
[docs]
@classmethod
def prism(cls, polygon, bottom=0.0, top=1.0, kind=None):
"""Return the prism over a 2D polygon, extruded along W.
The two caps are tagged ``"cap"`` and the side walls ``"side"``, so rules and face cells
can use them separately (for example a honeycomb's walls are its ``"side"`` faces).
Args:
polygon (array-like): ``(n, 2)`` U, V corners in boundary order. It may be concave.
bottom (float): W of the bottom cap, default ``0``.
top (float): W of the top cap, default ``1``.
kind (str | None): Cell type.
Returns:
Polyhedron: The prism.
Example:
>>> square = mm.Polyhedron.prism([(0, 0), (1, 0), (1, 1), (0, 1)])
"""
corners = np.asarray(polygon, dtype=np.float64).reshape(-1, 2)
count = len(corners)
if count < 3:
raise ValueError("a prism needs a polygon with at least three corners")
if not float(top) > float(bottom):
raise ValueError("top must be above bottom")
vertices = np.vstack([
np.column_stack([corners, np.full(count, float(bottom))]),
np.column_stack([corners, np.full(count, float(top))]),
])
faces = [tuple(range(count))[::-1], tuple(range(count, 2 * count))]
faces += [(index, (index + 1) % count, count + (index + 1) % count, count + index)
for index in range(count)]
tags = ["cap", "cap"] + ["side"] * count
return cls(vertices, faces, kind=kind, face_tags=tags)
@classmethod
def _from_parts(cls, vertices, faces, kind=None, face_tags=None, center=None, vertex_ids=None):
# Unvalidated cell for meshes: faces may be 2-vertex edges (the facets of a triangle).
cell = cls.__new__(cls)
cell._set(np.asarray(vertices, dtype=np.float64).reshape(-1, 3),
[tuple(face) for face in faces], kind, face_tags, center)
cell.vertex_ids = vertex_ids
return cell
def _set(self, points, faces, kind, face_tags, center):
self.vertices = points
self.faces = faces
self.kind = kind
if face_tags is None:
face_tags = [None] * len(faces)
self.face_tags = list(face_tags)
if len(self.face_tags) != len(faces):
raise ValueError("face_tags needs one tag per face")
if center is None:
used = sorted({index for face in faces for index in face})
center = points[used].mean(axis=0)
self.center = np.asarray(center, dtype=np.float64).reshape(3)
self.vertex_ids = None
def __repr__(self):
return "Polyhedron(vertices={}, faces={}, kind={!r})".format(
len(self.vertices), len(self.faces), self.kind)
[docs]
class PeriodicTessellation:
"""Space-filling cells that repeat with period one along U, V and W.
Hold one representative of each periodic class of cells, in normalized cell coordinates.
The cells must fill space face to face once repeated. :meth:`edges` turns them into a graph
unit cell with any of the :data:`GRAPH_RULES`, and :meth:`face_cell` turns their faces into
walls.
Args:
cells (list[Polyhedron]): One cell per periodic class.
cell_size (float | tuple): The U, V, W proportions of the physical cell the tessellation
was designed for, default ``(1, 1, 1)``. ``max_length`` is measured in these units.
Give the map cells the same proportions to keep the intended shape.
Example:
>>> tess = mm.tessellation("tet_oct")
>>> edges = tess.edges("cell_vertex_centroid", kinds="tet") # the fluorite net
>>> node = mm.graph_lattice(cell_map, edges, beam_radius=0.5)
"""
def __init__(self, cells, cell_size=(1.0, 1.0, 1.0)):
self.cells = list(cells)
if not self.cells:
raise ValueError("a tessellation needs at least one cell")
if not all(isinstance(cell, Polyhedron) for cell in self.cells):
raise TypeError("cells must be Polyhedron objects")
self.cell_size = tuple(float(value) for value in _cell_size(cell_size))
def __len__(self):
return len(self.cells)
def __iter__(self):
return iter(self.cells)
def __repr__(self):
return "PeriodicTessellation(cells={}, cell_size={})".format(len(self.cells), self.cell_size)
[docs]
def edges(self, rules="edges", kinds=None, face_tags=None, max_length=None):
"""Return the graph unit cell that ``rules`` make from the cells.
Struts are repeated periodically and clipped to the cell, every strut is split where
another strut ends on it, and duplicates are removed. Struts on the cell's faces are
kept on both opposite faces, like the other catalog cells; a mapped lattice counts the
copies that neighbouring cells share as one strut.
Args:
rules (str | list[str]): One of :data:`GRAPH_RULES`, or several for their union.
Default ``"edges"``.
kinds (str | list[str] | None): Only use cells of these kinds. ``None`` uses all.
face_tags (str | list | None): Only use faces with these tags. ``"edges"`` then takes
the edges of those faces, and ``"cell_vertex_centroid"`` their vertices.
max_length (float | None): Drop struts longer than this, measured in ``cell_size``
units before clipping.
Returns:
list: Edges ``((x0, y0, z0), (x1, y1, z1))`` in origin-centred coordinates spanning
``[-0.5, 0.5]``, sorted.
Example:
>>> bcc = mm.tessellation("cube").edges("cell_vertex_centroid")
"""
names = _rule_names(rules)
kind_filter = _selection(kinds)
tag_filter = _selection(face_tags)
adjacency = _periodic_adjacency(self.cells) if "dual" in names else []
segments = []
for name in names:
segments.extend(_rule_segments(self.cells, name, kind_filter, tag_filter, adjacency))
if max_length is not None:
limit = float(max_length)
if not math.isfinite(limit) or limit <= 0.0:
raise ValueError("max_length must be positive and finite")
scale = np.asarray(self.cell_size)
segments = [
segment for segment in segments
if np.linalg.norm((segment[1] - segment[0]) * scale) <= limit * (1.0 + 1e-9)
]
return _periodic_graph(segments)
[docs]
def face_cell(self, kinds=None, face_tags=None):
"""Return the cells' faces as walls: a ``pv.FaceUnitCell`` for :func:`face_lattice`.
Each shared face is used once. Faces are repeated periodically and clipped to the cell,
and polygons with more than four corners are triangulated. Every piece of one source
face gets that face's index as its face tag, and faces with the same plane direction
share an orientation tag.
Args:
kinds (str | list[str] | None): Only use faces of cells of these kinds.
face_tags (str | list | None): Only use faces with these tags, for example
``"side"`` for a prism honeycomb's walls.
Returns:
pv.FaceUnitCell: The wall cell, vertices in ``[0, 1]``.
Example:
>>> walls = mm.tessellation("triangular_prism").face_cell(face_tags="side")
>>> node = mm.face_lattice(cell_map, walls, wall_thickness=0.6)
"""
kind_filter = _selection(kinds)
tag_filter = _selection(face_tags)
polygons = []
seen = set()
for cell in self.cells:
if not _selected(cell.kind, kind_filter):
continue
for face, tag in zip(cell.faces, cell.face_tags):
if not _selected(tag, tag_filter):
continue
points = cell.vertices[list(face)]
shifted = points - np.floor(np.round(points.mean(axis=0), 7))
key = _point_set_key(shifted)
if key in seen:
continue
seen.add(key)
polygons.append(shifted)
if not polygons:
raise ValueError("no faces selected")
return _face_unit_cell(polygons)
[docs]
def periodic_voronoi(seeds, cell_size=(1.0, 1.0, 1.0)):
"""Return the Voronoi tessellation of a periodic seed set (A1 builder).
Any set of seeds inside one cell, repeated along U, V and W, divides space into Voronoi
cells: the region closer to each seed than to any other. Their edges are a strut network
that tiles without seams, so ``periodic_voronoi(seeds).edges()`` is a graph unit cell for
:func:`graph_lattice`. Regular seeds give the classic foams; random seeds give a random foam
that still tiles, so it can follow a curved CellMap.
================================ =====================================================
Seeds (fractional coordinates) Voronoi cells
================================ =====================================================
``(0.5, 0.5, 0.5)`` cubes: the ``"cubic"`` cell
BCC: corner + centre truncated octahedra: the ``"kelvin"`` cell
FCC: corner + face centres rhombic dodecahedra: ``"rhombic_dodecahedron"``
A15 (eight sites) the Weaire-Phelan foam: ``"weaire_phelan"``
================================ =====================================================
The tessellation's cell centres are the seeds, so ``.edges("dual")`` joins seeds whose
Voronoi cells share a face (the periodic Delaunay graph, where it is unique). For the bond
network of a crystal use :func:`periodic_bonds`, which joins seeds by distance.
Args:
seeds (array-like): ``(n, 3)`` distinct seed positions in fractional cell coordinates;
they are wrapped into ``[0, 1)``.
cell_size (float | tuple): U, V, W proportions of the physical cell. Distances are
measured in this metric, so a non-cubic cell gives the right Voronoi cells when the
map has the same proportions. Default ``(1, 1, 1)``.
Returns:
PeriodicTessellation: One Voronoi polyhedron per seed.
Raises:
ValueError: If a seed is not finite or two seeds coincide.
Example:
>>> fcc = [(0, 0, 0), (0.5, 0.5, 0), (0.5, 0, 0.5), (0, 0.5, 0.5)]
>>> edges = mm.periodic_voronoi(fcc).edges() # rhombic dodecahedra
>>> node = mm.graph_lattice(cell_map, edges, beam_radius=0.5)
"""
from scipy.spatial import Voronoi
points = np.asarray(seeds, dtype=np.float64).reshape(-1, 3)
if len(points) == 0:
raise ValueError("at least one seed is required")
if not np.all(np.isfinite(points)):
raise ValueError("seeds must be finite")
points = np.mod(np.round(points, 12), 1.0)
if len({_key(point) for point in points}) != len(points):
raise ValueError("seeds must be distinct after wrapping into the cell")
size = _cell_size(cell_size)
count = len(points)
tiled = (points[None, :, :] + _OFFSETS[:, None, :]).reshape(-1, 3)
central = int(np.flatnonzero(np.all(_OFFSETS == 0.0, axis=1))[0])
diagram = Voronoi(tiled * size)
# The copies of one Voronoi vertex in different repeats come out of different Delaunay
# simplices, so they differ by rounding; give every copy the same fractional position so the
# cell is exactly periodic. Co-spherical seeds give coincident vertices; merge them too.
vertices = diagram.vertices / size
rounded = np.round(vertices, 7)
repeat = np.floor(rounded)
fractions = {}
merged = {}
vertex_ids = np.empty(len(vertices), dtype=np.int64)
unique = []
for index, vertex in enumerate(vertices):
fraction_key = tuple(rounded[index] - repeat[index])
fraction = fractions.setdefault(fraction_key, vertex - repeat[index])
key = (tuple(repeat[index]), fraction_key)
if key not in merged:
merged[key] = len(unique)
unique.append(repeat[index] + fraction)
vertex_ids[index] = merged[key]
unique = np.asarray(unique)
faces = [[] for _ in range(count)]
for (first, second), ridge in zip(diagram.ridge_points, diagram.ridge_vertices):
if -1 in ridge:
continue
for seed, other in ((first, second), (second, first)):
if seed // count != central:
continue
polygon = []
for index in ridge:
merged_index = int(vertex_ids[index])
if not polygon or polygon[-1] != merged_index:
polygon.append(merged_index)
while len(polygon) > 1 and polygon[0] == polygon[-1]:
polygon.pop()
if len(set(polygon)) >= 3:
faces[seed % count].append(polygon)
cells = []
for seed in range(count):
if len(faces[seed]) < 4:
raise ValueError("could not build the Voronoi cell of seed {}".format(seed))
used = sorted({index for face in faces[seed] for index in face})
local = {index: position for position, index in enumerate(used)}
cells.append(Polyhedron(
unique[used],
[[local[index] for index in face] for face in faces[seed]],
kind="voronoi",
center=points[seed],
))
return PeriodicTessellation(cells, cell_size=size)
[docs]
def periodic_bonds(sites, max_length, cell_size=(1.0, 1.0, 1.0), min_length=0.0):
"""Return the graph cell that joins every two periodic sites closer than ``max_length``.
This is the bond network of a crystal: give its sites in one cell and the bond
length. It needs no tessellation, so it also works where many sites are equally
far apart (the Voronoi faces of such neighbours can shrink to an edge or a point).
Args:
sites (array-like): ``(n, 3)`` distinct site positions in fractional cell
coordinates; they are wrapped into ``[0, 1)``.
max_length (float): Longest bond, in ``cell_size`` units.
cell_size (float | tuple): U, V, W proportions of the physical cell, default
``(1, 1, 1)``.
min_length (float): Shortest bond kept, default ``0``.
Returns:
list: Origin-centred edges, as :meth:`PeriodicTessellation.edges` returns them.
Example:
>>> fcc = [(0, 0, 0), (0.5, 0.5, 0), (0.5, 0, 0.5), (0, 0.5, 0.5)]
>>> holes = [(x + 0.25, y + 0.25, z + 0.25) for x, y, z in fcc]
>>> diamond = mm.periodic_bonds(fcc + holes, max_length=0.44)
"""
points = np.asarray(sites, dtype=np.float64).reshape(-1, 3)
if len(points) == 0 or not np.all(np.isfinite(points)):
raise ValueError("sites must be finite and not empty")
points = np.mod(np.round(points, 12), 1.0)
if len({_key(point) for point in points}) != len(points):
raise ValueError("sites must be distinct after wrapping into the cell")
size = _cell_size(cell_size)
longest = float(max_length)
shortest = float(min_length)
if not math.isfinite(longest) or longest <= 0.0 or not 0.0 <= shortest < longest:
raise ValueError("need 0 <= min_length < max_length")
reach = [range(-int(math.ceil(longest / size[axis])) - 1, int(math.ceil(longest / size[axis])) + 2)
for axis in range(3)]
shifts = np.array(list(itertools.product(*reach)), dtype=np.float64)
segments = []
for index, start in enumerate(points):
ends = (points[None, :, :] + shifts[:, None, :]).reshape(-1, 3)
lengths = np.linalg.norm((ends - start) * size, axis=1)
keep = (lengths > max(shortest, 1e-12)) & (lengths <= longest * (1.0 + 1e-9))
segments.extend((start, end) for end in ends[keep])
return _periodic_graph(segments)
[docs]
def tessellation(name, **parameters):
"""Return a built-in periodic tessellation by name.
=============================== ==========================================================
Name Cells (prism caps tagged ``"cap"``, walls ``"side"``)
=============================== ==========================================================
``"cube"`` one cube (a square prism along W)
``"hexagonal_prism"`` four regular hexagonal prisms in the rectangular supercell
of ``"hex_prism_edge"`` (U:V = 2:sqrt(3))
``"reentrant_hexagonal_prism"`` the re-entrant version; takes ``angle_degrees`` and
``rib_ratio`` like :func:`reentrant_hex_prism`
``"triangular_prism"`` four triangular prisms, one edge family along U
(U:V = 1:sqrt(3))
``"rotated_triangular_prism"`` the same turned 90 degrees, edges along V (U:V = sqrt(3):1)
``"rotated_square_prism"`` two square prisms turned 45 degrees about W
``"tet_oct"`` the tetrahedral-octahedral honeycomb on FCC sites: eight
regular tetrahedra (kind ``"tet"``) and four octahedra
(kind ``"oct"``)
``"truncated_cube"`` a truncated cube and an octahedron (kinds
``"truncated_cube"`` and ``"oct"``)
=============================== ==========================================================
Args:
name (str): A name from :data:`TESSELLATION_NAMES`.
**parameters: Shape parameters of a parametric tessellation.
Returns:
PeriodicTessellation: The tessellation, with its design ``cell_size``.
Example:
>>> walls = mm.tessellation("hexagonal_prism").face_cell(face_tags="side")
"""
if name not in _TESSELLATIONS:
raise ValueError("Unknown tessellation '{}'. Available: {}".format(
name, ", ".join(TESSELLATION_NAMES)))
return _TESSELLATIONS[name](**parameters)
def _cube():
return PeriodicTessellation([Polyhedron.prism([(0, 0), (1, 0), (1, 1), (0, 1)], kind="cube")])
def _honeycomb_polygons(reentrant=False, angle_degrees=-30.0, rib_ratio=2.0):
# The hexagons of the rectangular honeycomb supercell used by the prism cells; the levels
# come from the same parameterization as unit_cells._honeycomb_points.
from .unit_cells import _honeycomb_points
points = _honeycomb_points(reentrant, angle_degrees, rib_ratio)
inner_high, inner_low = points[0][1], points[1][1]
boundary_low, boundary_high = points[5][1], points[6][1]
polygons = []
for x in (0.5, 0.0):
polygons.append([(x - 0.25, inner_low), (x, boundary_low), (x + 0.25, inner_low),
(x + 0.25, inner_high), (x, boundary_high), (x - 0.25, inner_high)])
for x in (0.25, 0.75):
polygons.append([(x - 0.25, boundary_high - 1.0), (x, inner_high - 1.0),
(x + 0.25, boundary_high - 1.0), (x + 0.25, boundary_low),
(x, inner_low), (x - 0.25, boundary_low)])
return polygons
def _hexagonal_prism():
return PeriodicTessellation(
[Polyhedron.prism(polygon, kind="hexagonal_prism") for polygon in _honeycomb_polygons()],
cell_size=(2.0 / math.sqrt(3.0), 1.0, 1.0),
)
def _reentrant_hexagonal_prism(angle_degrees=-30.0, rib_ratio=2.0):
from .unit_cells import reentrant_hex_prism_reference_aspect_ratio
ratio = reentrant_hex_prism_reference_aspect_ratio(angle_degrees, rib_ratio)
polygons = _honeycomb_polygons(True, angle_degrees, rib_ratio)
return PeriodicTessellation(
[Polyhedron.prism(polygon, kind="reentrant_hexagonal_prism") for polygon in polygons],
cell_size=(ratio, 1.0, 1.0),
)
_TRIANGLES = (
((0.0, 0.0), (1.0, 0.0), (0.5, 0.5)),
((0.5, 0.5), (1.0, 0.0), (1.5, 0.5)),
((0.0, 1.0), (0.5, 0.5), (1.0, 1.0)),
((0.5, 0.5), (1.5, 0.5), (1.0, 1.0)),
)
def _triangular_prism():
return PeriodicTessellation(
[Polyhedron.prism(triangle, kind="triangular_prism") for triangle in _TRIANGLES],
cell_size=(1.0, math.sqrt(3.0), 1.0),
)
def _rotated_triangular_prism():
turned = [[(v, u) for u, v in triangle] for triangle in _TRIANGLES]
return PeriodicTessellation(
[Polyhedron.prism(triangle, kind="triangular_prism") for triangle in turned],
cell_size=(math.sqrt(3.0), 1.0, 1.0),
)
def _rotated_square_prism():
squares = (
((0.0, 0.0), (0.5, -0.5), (1.0, 0.0), (0.5, 0.5)),
((0.0, 0.0), (0.5, 0.5), (0.0, 1.0), (-0.5, 0.5)),
)
return PeriodicTessellation([Polyhedron.prism(square, kind="square_prism") for square in squares])
def _tet_oct():
cells = []
for corner in itertools.product((0.0, 1.0), repeat=3):
face_centers = []
for axis in range(3):
center = [0.5, 0.5, 0.5]
center[axis] = corner[axis]
face_centers.append(tuple(center))
cells.append(Polyhedron([corner] + face_centers, kind="tet"))
for center in ((0.5, 0.5, 0.5), (0.5, 0.0, 0.0), (0.0, 0.5, 0.0), (0.0, 0.0, 0.5)):
corners = []
for axis in range(3):
for sign in (-0.5, 0.5):
point = list(center)
point[axis] += sign
corners.append(tuple(point))
cells.append(Polyhedron(corners, kind="oct"))
return PeriodicTessellation(cells)
def _truncated_cube():
short = 0.5 * (math.sqrt(2.0) - 1.0) # half the octagon edge
corner = 0.5 - short # octahedron vertex distance from the cube corner
truncated = set()
for axis in range(3):
for first in (-short, short):
for second, third in itertools.product((-0.5, 0.5), repeat=2):
point = [0.0, 0.0, 0.0]
point[axis] = first
point[(axis + 1) % 3] = second
point[(axis + 2) % 3] = third
truncated.add(tuple(0.5 + value for value in point))
octahedron = []
for axis in range(3):
for sign in (-corner, corner):
point = [0.0, 0.0, 0.0]
point[axis] = sign
octahedron.append(tuple(point))
return PeriodicTessellation([
Polyhedron(sorted(truncated), kind="truncated_cube"),
Polyhedron(octahedron, kind="oct"),
])
_TESSELLATIONS = {
"cube": _cube,
"hexagonal_prism": _hexagonal_prism,
"reentrant_hexagonal_prism": _reentrant_hexagonal_prism,
"triangular_prism": _triangular_prism,
"rotated_triangular_prism": _rotated_triangular_prism,
"rotated_square_prism": _rotated_square_prism,
"tet_oct": _tet_oct,
"truncated_cube": _truncated_cube,
}
#: Names accepted by :func:`tessellation`.
TESSELLATION_NAMES = tuple(_TESSELLATIONS)
# ---- rules ------------------------------------------------------------------------------------
def _rule_names(rules):
names = (rules,) if isinstance(rules, str) else tuple(rules)
if not names:
raise ValueError("at least one rule is required")
for name in names:
if name not in GRAPH_RULES:
raise ValueError("Unknown rule '{}'. Available: {}".format(name, ", ".join(GRAPH_RULES)))
return tuple(dict.fromkeys(names))
def _selection(value):
if value is None:
return None
if isinstance(value, (str, numbers.Number)):
return {value}
return set(value)
def _selected(value, selection):
return selection is None or value in selection
def _face_center(cell, face):
return cell.vertices[list(face)].mean(axis=0)
def _rule_segments(cells, rule, kinds, face_tags, adjacency):
"""Return the struts of one rule as ``(start, end)`` point pairs."""
segments = []
if rule == "dual":
for first, first_face, second, _, shift in adjacency:
if not (_selected(cells[first].kind, kinds) and _selected(cells[second].kind, kinds)):
continue
if not _selected(cells[first].face_tags[first_face], face_tags):
continue
segments.append((cells[first].center, cells[second].center + shift))
return segments
for cell in cells:
if not _selected(cell.kind, kinds):
continue
faces = [face for face, tag in zip(cell.faces, cell.face_tags) if _selected(tag, face_tags)]
if rule == "edges":
seen = set()
for face in faces:
pairs = [face] if len(face) == 2 else zip(face, face[1:] + face[:1])
for start, end in pairs:
key = (min(start, end), max(start, end))
if key not in seen:
seen.add(key)
segments.append((cell.vertices[start], cell.vertices[end]))
elif rule == "cell_vertex_centroid":
for index in sorted({index for face in faces for index in face}):
segments.append((cell.center, cell.vertices[index]))
elif rule == "facet_vertex_centroid":
for face in faces:
center = _face_center(cell, face)
for index in face:
segments.append((center, cell.vertices[index]))
elif rule == "central_axes":
for face in faces:
segments.append((cell.center, _face_center(cell, face)))
return segments
def _periodic_adjacency(cells):
"""Return ``(cell, face, other_cell, other_face, shift)`` for every shared face.
``other_cell``'s face plus ``shift`` (a whole number of repeats) is ``cell``'s face.
"""
groups = {}
for cell_index, cell in enumerate(cells):
for face_index, face in enumerate(cell.faces):
points = cell.vertices[list(face)]
center = points.mean(axis=0)
key = tuple(np.mod(np.round(center, 7), 1.0).round(7))
groups.setdefault(key, []).append((cell_index, face_index, center, points))
pairs = []
for entries in groups.values():
for first, second in itertools.combinations(entries, 2):
shift = np.round(first[2] - second[2])
if _point_set_key(second[3] + shift) == _point_set_key(first[3]):
pairs.append((first[0], first[1], second[0], second[1], shift))
return pairs
# ---- periodic finishing -------------------------------------------------------------------------
def _key(point):
return tuple(float(value) + 0.0 for value in np.round(np.asarray(point, dtype=np.float64), _DIGITS))
def _point_set_key(points):
return tuple(sorted(_key(point) for point in np.asarray(points).reshape(-1, 3)))
def _cell_size(value):
if isinstance(value, numbers.Real):
size = np.full(3, float(value))
else:
size = np.asarray([float(item) for item in value], dtype=np.float64)
if size.shape != (3,) or not np.all(np.isfinite(size)) or np.any(size <= 0.0):
raise ValueError("cell_size must be a positive number or three positive numbers")
return size
def _snap(points):
points = np.where(np.abs(points) < _SNAP, 0.0, points)
points = np.where(np.abs(points - 1.0) < _SNAP, 1.0, points)
return np.clip(points, 0.0, 1.0)
def _clip_to_cell(segments):
"""Repeat segments periodically and clip every copy to the closed unit cell."""
starts = segments[:, 0, :]
directions = segments[:, 1, :] - starts
ranges = []
for axis in range(3):
low = int(math.floor(-float(np.max(segments[:, :, axis])) - _SNAP))
high = int(math.ceil(1.0 - float(np.min(segments[:, :, axis])) + _SNAP))
ranges.append(range(low, high + 1))
shifts = np.array(list(itertools.product(*ranges)), dtype=np.float64)
origin = starts[None, :, :] + shifts[:, None, :]
direction = np.broadcast_to(directions[None, :, :], origin.shape)
first = np.zeros(origin.shape[:2])
last = np.ones(origin.shape[:2])
# Crossings are computed on the exact faces, so the copies of a strut that leave through one
# face and enter through the opposite one give the same point; the tolerance only keeps
# struts that lie in a face.
with np.errstate(divide="ignore", invalid="ignore"):
for axis in range(3):
start = origin[..., axis]
step = direction[..., axis]
flat = np.abs(step) < 1e-12
inside = (start >= -_SNAP) & (start <= 1.0 + _SNAP)
enter = -start / step
leave = (1.0 - start) / step
low = np.where(flat, np.where(inside, -np.inf, np.inf), np.minimum(enter, leave))
high = np.where(flat, np.where(inside, np.inf, -np.inf), np.maximum(enter, leave))
first = np.maximum(first, low)
last = np.minimum(last, high)
keep = first < last - 1e-12
begin = _snap(origin[keep] + first[keep][:, None] * direction[keep])
end = _snap(origin[keep] + last[keep][:, None] * direction[keep])
return begin, end
def _weld(pairs, tolerance=1e-7):
"""Split each segment wherever another segment's end lies on it, and drop duplicates."""
ends = sorted({point for pair in pairs for point in pair})
if not ends:
return []
points = np.asarray(ends, dtype=np.float64)
result = set()
for start, end in pairs:
first = np.asarray(start)
step = np.asarray(end) - first
length_squared = float(step @ step)
if length_squared <= tolerance * tolerance:
continue
along = ((points - first) @ step) / length_squared
offset = points - (first + along[:, None] * step)
margin = tolerance / math.sqrt(length_squared)
on = (np.einsum("ij,ij->i", offset, offset) < tolerance * tolerance) \
& (along > margin) & (along < 1.0 - margin)
chain = [start] + [ends[index] for index in np.flatnonzero(on)[np.argsort(along[on])]] + [end]
for first_point, second_point in zip(chain, chain[1:]):
if first_point != second_point:
result.add((min(first_point, second_point), max(first_point, second_point)))
return sorted(result)
def _periodic_graph(segments):
if not segments:
return []
array = np.asarray([(np.asarray(start), np.asarray(end)) for start, end in segments],
dtype=np.float64)
begin, end = _clip_to_cell(array)
pairs = []
for first, second in zip(begin, end):
first_key = _key(first)
second_key = _key(second)
if first_key != second_key and np.linalg.norm(second - first) > 1e-7:
pairs.append((first_key, second_key))
return [
(tuple(_key(np.asarray(start) - 0.5)), tuple(_key(np.asarray(end) - 0.5)))
for start, end in _weld(pairs)
]
# ---- faces ------------------------------------------------------------------------------------
def _hull_faces(points):
from scipy.spatial import ConvexHull
hull = ConvexHull(points)
groups = {}
for simplex, equation in zip(hull.simplices, hull.equations):
groups.setdefault(tuple(np.round(equation, 7)), set()).update(int(index) for index in simplex)
used = sorted(int(index) for index in hull.vertices)
local = {index: position for position, index in enumerate(used)}
faces = []
for key, indices in groups.items():
indices = sorted(index for index in indices if index in local)
faces.append([local[index] for index in _ordered(points, indices, np.asarray(key[:3]))])
return points[used], faces
def _ordered(points, indices, normal):
"""Order coplanar points around their centre."""
polygon = points[indices]
center = polygon.mean(axis=0)
reference = np.eye(3)[int(np.argmin(np.abs(normal)))]
tangent = np.cross(normal, reference)
tangent /= np.linalg.norm(tangent)
bitangent = np.cross(normal, tangent)
angles = np.arctan2((polygon - center) @ bitangent, (polygon - center) @ tangent)
return [indices[position] for position in np.argsort(angles)]
def _polygon_normal(polygon):
# Newell's method: robust for concave and nearly degenerate polygons.
normal = np.zeros(3)
for current, following in zip(polygon, np.roll(polygon, -1, axis=0)):
normal += np.cross(current, following)
return normal
def _clip_polygon(polygon):
"""Sutherland-Hodgman clip of a planar polygon to the closed unit cell."""
points = list(polygon)
for axis in range(3):
for bound, sign in ((0.0, 1.0), (1.0, -1.0)):
if not points:
return []
clipped = []
for current, following in zip(points, points[1:] + points[:1]):
current_inside = sign * (current[axis] - bound) >= -_SNAP
following_inside = sign * (following[axis] - bound) >= -_SNAP
if current_inside:
clipped.append(current)
if current_inside != following_inside:
fraction = (bound - current[axis]) / (following[axis] - current[axis])
clipped.append(current + fraction * (following - current))
points = clipped
result = []
for point in points:
snapped = _snap(point)
if not result or np.linalg.norm(snapped - result[-1]) > 1e-9:
result.append(snapped)
while len(result) > 1 and np.linalg.norm(result[0] - result[-1]) <= 1e-9:
result.pop()
return result
def _triangulate(polygon):
"""Split a planar polygon into triangles (ear clipping), or keep it if it has 3-4 corners."""
count = len(polygon)
normal = _polygon_normal(polygon)
if count <= 4 and _is_convex(polygon, normal):
return [list(range(count))]
# Drop the dominant axis and keep the polygon counter-clockwise in the remaining two.
axis = int(np.argmax(np.abs(normal)))
flat = np.delete(np.asarray(polygon), axis, axis=1)
if (normal[axis] if axis != 1 else -normal[axis]) < 0.0:
flat = flat[:, ::-1]
remaining = list(range(count))
triangles = []
guard = 0
while len(remaining) > 3 and guard < 10 * count:
guard += 1
for position in range(len(remaining)):
previous = remaining[position - 1]
current = remaining[position]
following = remaining[(position + 1) % len(remaining)]
if _cross_2d(flat[previous], flat[current], flat[following]) <= 1e-14:
continue
if any(_in_triangle(flat[other], flat[previous], flat[current], flat[following])
for other in remaining if other not in (previous, current, following)):
continue
triangles.append([previous, current, following])
remaining.pop(position)
break
else:
break
if len(remaining) == 3:
triangles.append(remaining)
return triangles
def _cross_2d(first, second, third):
return (second[0] - first[0]) * (third[1] - first[1]) - (second[1] - first[1]) * (third[0] - first[0])
def _in_triangle(point, first, second, third):
return (_cross_2d(first, second, point) >= 0.0 and _cross_2d(second, third, point) >= 0.0
and _cross_2d(third, first, point) >= 0.0)
def _is_convex(polygon, normal):
count = len(polygon)
for index in range(count):
edge = polygon[(index + 1) % count] - polygon[index]
following = polygon[(index + 2) % count] - polygon[(index + 1) % count]
if np.dot(np.cross(edge, following), normal) < -1e-14:
return False
return True
def _face_unit_cell(polygons):
vertices = {}
faces = []
face_tags = []
orientation_tags = []
directions = {}
seen = set()
def vertex(point):
key = _key(point)
if key not in vertices:
vertices[key] = len(vertices)
return vertices[key]
for source, polygon in enumerate(polygons):
normal = _polygon_normal(polygon)
length = np.linalg.norm(normal)
if length <= 1e-12:
continue
direction = normal / length
if next((value for value in direction if abs(value) > 1e-9), 0.0) < 0.0:
direction = -direction
orientation = directions.setdefault(_key(np.round(direction, 6)), len(directions))
minimum = np.min(polygon, axis=0)
maximum = np.max(polygon, axis=0)
ranges = [range(int(math.floor(-maximum[axis] - _SNAP)), int(math.ceil(1.0 - minimum[axis] + _SNAP)) + 1)
for axis in range(3)]
for shift in itertools.product(*ranges):
piece = _clip_polygon([point + np.asarray(shift, dtype=np.float64) for point in polygon])
if len(piece) < 3 or np.linalg.norm(_polygon_normal(piece)) <= 1e-10:
continue
key = _point_set_key(piece)
if key in seen:
continue
seen.add(key)
for triangle in _triangulate(piece):
faces.append([vertex(piece[index]) for index in triangle])
face_tags.append(source)
orientation_tags.append(orientation)
if not faces:
raise ValueError("no face lies inside the cell")
ordered = sorted(vertices, key=vertices.get)
return pv.FaceUnitCell(
[pv.Vec3(*point) for point in ordered],
faces,
face_tags=face_tags,
orientation_tags=orientation_tags,
)
