"""Unit-cell edge generators for strut-based lattices.
Each generator returns a list of edges for a single unit cell as
``[((x0, y0, z0), (x1, y1, z1)), ...]`` in **normalized, origin-centered**
coordinates that span ``[-0.5, 0.5]`` on every axis. ``lattices.py`` shifts them
into the native ``GraphUnitCell`` interval ``[0, 1]``; the mapped C++ node then
tiles and deforms them through the supplied ``CellMap``.
The four classic cells (``cubic``, ``bcc``, ``fcc``, ``kelvin``) mirror the edge
definitions used by the C++ ``GraphLattice`` presets in
``core-lib/src/tree/leaf/lattice/graph_lattice.cpp`` so results are consistent
with the built-in node. Additional crystallographic and engineering presets are
defined with explicit junctions so graph metadata agrees with the geometry. Most
of the newer cells are a rule applied to a tessellation (``tessellations.py``):
the docstring of each says which.
"""
# STD
import math
# Local
from .tessellations import _periodic_graph, periodic_bonds, periodic_voronoi, tessellation
_H = 0.5 # half cell (normalized)
_Q = 0.25 # quarter cell (normalized)
# Eight cube corners, indexed to match the C++ GraphLattice vertex numbering.
_CORNERS = [
(-_H, -_H, -_H), # v0
(_H, -_H, -_H), # v1
(-_H, _H, -_H), # v2
(_H, _H, -_H), # v3
(-_H, -_H, _H), # v4
(_H, -_H, _H), # v5
(-_H, _H, _H), # v6
(_H, _H, _H), # v7
]
def cubic():
"""Simple cubic cell: the twelve edges of the cube."""
v = _CORNERS
idx = [
(0, 1), (1, 3), (3, 2), (2, 0), # bottom face
(4, 5), (5, 7), (7, 6), (6, 4), # top face
(0, 4), (1, 5), (2, 6), (3, 7), # vertical connectors
]
return [(v[a], v[b]) for a, b in idx]
def bcc():
"""Body-centered cubic cell: the four body diagonals (corner to opposite
corner through the center)."""
v = _CORNERS
idx = [(0, 7), (1, 6), (2, 5), (3, 4)]
return [(v[a], v[b]) for a, b in idx]
def fcc():
"""Face-centered cubic cell: the twelve face diagonals."""
v = _CORNERS
idx = [
(0, 6), (0, 5), (0, 3),
(1, 7), (1, 4), (1, 2),
(2, 7), (2, 4),
(3, 5), (3, 6),
(5, 6), (7, 4),
]
return [(v[a], v[b]) for a, b in idx]
def kelvin():
"""Kelvin cell (truncated octahedron) edges."""
# Vertices grouped by face, in origin-centered normalized coordinates.
v = [
(_Q, 0.0, _H), (0.0, _Q, _H), (-_Q, 0.0, _H), (0.0, -_Q, _H), # top
(_Q, 0.0, -_H), (0.0, _Q, -_H), (-_Q, 0.0, -_H), (0.0, -_Q, -_H), # bottom
(_H, _Q, 0.0), (_H, 0.0, _Q), (_H, -_Q, 0.0), (_H, 0.0, -_Q), # right
(-_H, _Q, 0.0), (-_H, 0.0, _Q), (-_H, -_Q, 0.0), (-_H, 0.0, -_Q), # left
(_Q, _H, 0.0), (0.0, _H, _Q), (-_Q, _H, 0.0), (0.0, _H, -_Q), # front
(_Q, -_H, 0.0), (0.0, -_H, _Q), (-_Q, -_H, 0.0), (0.0, -_H, -_Q), # back
]
idx = [
(0, 1), (1, 2), (2, 3), (3, 0), # top square
(4, 5), (5, 6), (6, 7), (7, 4), # bottom square
(8, 9), (9, 10), (10, 11), (11, 8), # right square
(12, 13), (13, 14), (14, 15), (15, 12), # left square
(16, 17), (17, 18), (18, 19), (19, 16), # front square
(20, 21), (21, 22), (22, 23), (23, 20), # back square
(0, 9), (3, 21), (2, 13), (1, 17), # top connectors
(4, 11), (5, 19), (6, 15), (7, 23), # bottom connectors
(14, 22), (12, 18), (8, 16), (10, 20), # side connectors
]
return [(v[a], v[b]) for a, b in idx]
def octet():
"""Octet truss cell: face-center to corner struts (tetrahedra) plus the
octahedron formed by the six face centers."""
corners = _CORNERS
faces = {
"x+": (_H, 0.0, 0.0),
"x-": (-_H, 0.0, 0.0),
"y+": (0.0, _H, 0.0),
"y-": (0.0, -_H, 0.0),
"z+": (0.0, 0.0, _H),
"z-": (0.0, 0.0, -_H),
}
edges = []
# Each face center to the four corners lying on that face.
axis_of = {"x+": (0, _H), "x-": (0, -_H), "y+": (1, _H),
"y-": (1, -_H), "z+": (2, _H), "z-": (2, -_H)}
for key, center in faces.items():
axis, value = axis_of[key]
for corner in corners:
if abs(corner[axis] - value) < 1e-9:
edges.append((center, corner))
# Octahedron: every face center connects to every non-opposite face center.
opposite = {"x+": "x-", "x-": "x+", "y+": "y-",
"y-": "y+", "z+": "z-", "z-": "z+"}
keys = list(faces)
for i in range(len(keys)):
for j in range(i + 1, len(keys)):
if opposite[keys[i]] != keys[j]:
edges.append((faces[keys[i]], faces[keys[j]]))
return edges
def diamond():
"""Diamond cubic cell: tetrahedral bonds between the two interpenetrating
FCC sub-lattices."""
# Defined in the [0, 1] cube, then shifted to be origin-centered.
tetra_bonds = [
((0.25, 0.25, 0.25), [(0, 0, 0), (0.5, 0.5, 0), (0.5, 0, 0.5), (0, 0.5, 0.5)]),
((0.75, 0.75, 0.25), [(0.5, 0.5, 0), (1, 1, 0), (1, 0.5, 0.5), (0.5, 1, 0.5)]),
((0.75, 0.25, 0.75), [(0.5, 0, 0.5), (1, 0.5, 0.5), (1, 0, 1), (0.5, 0.5, 1)]),
((0.25, 0.75, 0.75), [(0, 0.5, 0.5), (0.5, 1, 0.5), (0.5, 0.5, 1), (0, 1, 1)]),
]
def shift(p):
return (p[0] - 0.5, p[1] - 0.5, p[2] - 0.5)
edges = []
for node, neighbors in tetra_bonds:
for neighbor in neighbors:
edges.append((shift(node), shift(neighbor)))
return edges
def isotruss():
"""SC+BCC+SCC Isotruss topology with every center junction explicit.
The topology combines the twelve cube edges, eight center-to-corner struts,
and six center-to-face-center struts. The name identifies the graph
topology; isotropic mechanical behavior still depends on member radii.
"""
center = (0.0, 0.0, 0.0)
face_centers = [
(-_H, 0.0, 0.0), (_H, 0.0, 0.0),
(0.0, -_H, 0.0), (0.0, _H, 0.0),
(0.0, 0.0, -_H), (0.0, 0.0, _H),
]
return cubic() + [
(center, corner) for corner in _CORNERS
] + [
(center, face_center) for face_center in face_centers
]
def fluorite():
"""Fluorite crystal net: FCC calcium sites joined to tetrahedral sites."""
calcium = list(_CORNERS) + [
(-_H, 0.0, 0.0), (_H, 0.0, 0.0),
(0.0, -_H, 0.0), (0.0, _H, 0.0),
(0.0, 0.0, -_H), (0.0, 0.0, _H),
]
tetrahedral = [
(x, y, z)
for x in (-_Q, _Q)
for y in (-_Q, _Q)
for z in (-_Q, _Q)
]
edges = []
target_squared_length = 3.0 / 16.0
for site in tetrahedral:
for neighbor in calcium:
squared_length = sum(
(site[axis] - neighbor[axis]) ** 2 for axis in range(3)
)
if abs(squared_length - target_squared_length) < 1e-12:
edges.append((site, neighbor))
return edges
_HONEYCOMB_EDGE_INDICES = [
(3, 2), (3, 5), (3, 9),
(1, 0), (1, 9), (1, 13),
(10, 0), (10, 2), (10, 11),
(6, 7), (14, 15), (4, 5), (8, 9), (13, 12),
(0, 14), (2, 6),
]
def _honeycomb_points(reentrant=False, angle_degrees=-30.0, rib_ratio=2.0):
if reentrant:
angle = math.radians(float(angle_degrees))
rib_ratio = float(rib_ratio)
if not -90.0 < float(angle_degrees) < 0.0:
raise ValueError("angle_degrees must lie strictly between -90 and 0")
if not math.isfinite(rib_ratio) or rib_ratio <= 0.0:
raise ValueError("rib_ratio must be positive and finite")
inset = -math.sin(angle)
height = rib_ratio + 2.0 * inset
inner_low = inset / height
inner_high = (rib_ratio + inset) / height
boundary_low = 2.0 * inset / height
boundary_high = rib_ratio / height
else:
inner_low = 1.0 / 3.0
inner_high = 2.0 / 3.0
boundary_low = 1.0 / 6.0
boundary_high = 5.0 / 6.0
return [
(0.75, inner_high), (0.75, inner_low),
(0.25, inner_high), (0.25, inner_low),
(0.0, 0.0), (0.0, boundary_low),
(0.0, boundary_high), (0.0, 1.0),
(0.5, 0.0), (0.5, boundary_low),
(0.5, boundary_high), (0.5, 1.0),
(1.0, 0.0), (1.0, boundary_low),
(1.0, boundary_high), (1.0, 1.0),
]
def hex_prism_edge():
"""Regular hexagonal-prism edge network: the edges of ``tessellation("hexagonal_prism")``.
The U/V unit-cell pitch should use the reference aspect ratio returned by
:func:`hex_prism_reference_aspect_ratio` to preserve regular hexagons. W
struts stand only at hexagon corners.
"""
return tessellation("hexagonal_prism").edges()
def reentrant_hex_prism(angle_degrees=-30.0, rib_ratio=2.0):
"""Extruded re-entrant hexagonal graph: the edges of the re-entrant prism tessellation.
``rib_ratio`` is the vertical-to-inclined member length ratio. The catalog
defaults follow the foldable equal-height geometry h/l=2 at -30 degrees. W
struts stand only at hexagon corners.
"""
return tessellation(
"reentrant_hexagonal_prism",
angle_degrees=angle_degrees,
rib_ratio=rib_ratio,
).edges()
# Struts of the "column" rule on a prism: the axis through the centre and spokes from the
# centre of each cap to its corners.
_COLUMN_RULES = ("central_axes", "facet_vertex_centroid")
_FCC_SITES = ((0.0, 0.0, 0.0), (0.5, 0.5, 0.0), (0.5, 0.0, 0.5), (0.0, 0.5, 0.5))
_A15_SITES = (
(0.0, 0.0, 0.0), (0.5, 0.5, 0.5),
(0.25, 0.5, 0.0), (0.75, 0.5, 0.0),
(0.0, 0.25, 0.5), (0.0, 0.75, 0.5),
(0.5, 0.0, 0.25), (0.5, 0.0, 0.75),
)
# Hexagonal crystal cells use the rectangular supercell of the hexagonal prism cells: U = 2a,
# V = sqrt(3) a, W = c, with the ideal c / a = sqrt(8 / 3) of close packing.
_C_OVER_A = math.sqrt(8.0 / 3.0)
_HEXAGONAL_CELL_SIZE = (2.0 / math.sqrt(3.0), 1.0, _C_OVER_A / math.sqrt(3.0))
def _hexagonal_sites(fractional):
"""Map sites given in a hexagonal cell (a1, a2 at 120 degrees, c) into the supercell."""
sites = set()
for f1, f2, f3 in fractional:
for n1 in range(-3, 4):
for n2 in range(-3, 4):
# Cartesian in units of a, shifted half a cell along U so the atom columns of
# the hexagonal diamond sit on the corners of hex_prism_edge's hexagons.
x = (f1 + n1) - 0.5 * (f2 + n2) + 0.5
y = 0.5 * math.sqrt(3.0) * (f2 + n2)
u = x / 2.0
v = y / math.sqrt(3.0)
if -1e-9 <= u < 1.0 - 1e-9 and -1e-9 <= v < 1.0 - 1e-9:
sites.add((round(u, 12) + 0.0, round(v, 12) + 0.0, round(f3 % 1.0, 12) + 0.0))
return sorted(sites)
def _bond_net(fractional, bond_length):
"""Nearest-neighbour bonds of hexagonal sites, ``bond_length`` in units of a."""
unit = _HEXAGONAL_CELL_SIZE[1] / math.sqrt(3.0) # a in the cell_size units (V = 1)
return periodic_bonds(_hexagonal_sites(fractional), 1.02 * bond_length * unit,
cell_size=_HEXAGONAL_CELL_SIZE)
def rhombic_dodecahedron():
"""Rhombic-dodecahedron foam: the Voronoi edges of FCC sites (:func:`periodic_voronoi`)."""
return periodic_voronoi(_FCC_SITES).edges()
def weaire_phelan():
"""Weaire-Phelan foam: the Voronoi edges of the eight A15 sites.
Two pentagonal dodecahedra and six tetrakaidecahedra fill each cell; every
interior node joins four struts.
"""
return periodic_voronoi(_A15_SITES).edges()
def truncated_cube():
"""Edges of the truncated cubic honeycomb: a truncated cube centred in the cell, with
octahedra at the cell corners. Every strut lies in a cell face plane."""
return tessellation("truncated_cube").edges()
def column():
"""One W column per cell, tied at the top and bottom by X struts in the cell faces.
The column rule on ``tessellation("cube")``: the axis through the centre and
spokes from each cap's centre to its corners.
"""
return tessellation("cube").edges(_COLUMN_RULES, face_tags="cap")
def columns():
"""Two W columns per cell, tied by a half-pitch square grid in every layer.
The column rule of :func:`column` on ``tessellation("rotated_square_prism")``,
so the columns stand at the centres of the cell's U and V edges.
"""
return tessellation("rotated_square_prism").edges(_COLUMN_RULES, face_tags="cap")
def reentrant_3d(angle_degrees=-30.0, rib_ratio=2.0):
"""Three-dimensional re-entrant (auxetic) strut cell.
Two orthogonal families of the planar re-entrant strut pattern of
:func:`reentrant_hex_prism`, turned so that the re-entrant direction is W:
U-W walls at V = 1/4 and 3/4 and V-W walls at U = 1/4 and 3/4. The walls
share their full W ribs, so compression along W pulls the walls in along
both U and V. Give the map the U:V:W proportions from
:func:`reference_cell_size` with the same parameters.
Args:
angle_degrees (float): Inclined-member angle in degrees, strictly
between ``-90`` and ``0``, default ``-30``.
rib_ratio (float): W-rib to inclined member length ratio, positive,
default ``2.0``.
Returns:
list: Origin-centred edges.
"""
import numpy as np
points = _honeycomb_points(reentrant=True, angle_degrees=angle_degrees, rib_ratio=rib_ratio)
segments = []
for offset in (0.25, 0.75):
for start, end in _HONEYCOMB_EDGE_INDICES:
first, second = points[start], points[end]
segments.append((np.array((first[0], offset, first[1])), np.array((second[0], offset, second[1]))))
segments.append((np.array((offset, first[0], first[1])), np.array((offset, second[0], second[1]))))
return _periodic_graph(segments)
def triangular_prism_edge():
"""Triangular-prism edge network (a strut triangular honeycomb), one edge family along U.
Edges of ``tessellation("triangular_prism")``; keep U:V = 1:sqrt(3) for
equilateral triangles.
"""
return tessellation("triangular_prism").edges()
def rotated_triangular_prism_edge():
"""Triangular-prism edge network turned 90 degrees: one edge family along V.
Keep U:V = sqrt(3):1 for equilateral triangles.
"""
return tessellation("rotated_triangular_prism").edges()
def rotated_square_prism_edge():
"""Square-prism edge network turned 45 degrees about W: X struts in every
layer and W struts at the cell's vertical edges and centre."""
return tessellation("rotated_square_prism").edges()
def hex_prism_vertex_centroid():
"""Hexagonal prisms with each prism's centre joined to its twelve corners
(``"cell_vertex_centroid"`` on ``tessellation("hexagonal_prism")``)."""
return tessellation("hexagonal_prism").edges("cell_vertex_centroid")
def hex_prism_central_axis():
"""Hexagonal-prism edges plus a central W column in every prism, whose ends
are joined to the cap corners by spokes (the column rule of :func:`column`)."""
hexagonal = tessellation("hexagonal_prism")
return _periodic_union(hexagonal.edges(), hexagonal.edges(_COLUMN_RULES, face_tags="cap"))
def hex_prism_diamond():
"""Hexagonal diamond (lonsdaleite): tetrahedral bonds, with puckered
hexagons in every layer joined by W bonds.
The atom columns sit on the corners of :func:`hex_prism_edge`'s hexagons.
Ideal bond angles need the U:V:W proportions from
:func:`reference_cell_size`.
"""
sites = ((1.0 / 3.0, 2.0 / 3.0, 0.0), (1.0 / 3.0, 2.0 / 3.0, 3.0 / 8.0),
(2.0 / 3.0, 1.0 / 3.0, 0.5), (2.0 / 3.0, 1.0 / 3.0, 7.0 / 8.0))
return _bond_net(sites, 3.0 / 8.0 * _C_OVER_A)
def hex_prism_laves():
"""Hexagonal Laves-phase net: the bonds between the smaller atoms of the C14
structure, a network of corner- and face-sharing tetrahedra (kagome layers
joined through triangles). Uses the proportions from :func:`reference_cell_size`.
"""
x = -1.0 / 6.0
sites = [(0.0, 0.0, 0.0), (0.0, 0.0, 0.5),
(x, 2.0 * x, 0.25), (-2.0 * x, -x, 0.25), (x, -x, 0.25),
(-x, -2.0 * x, 0.75), (2.0 * x, x, 0.75), (-x, x, 0.75)]
return _bond_net(sites, 0.5)
def tet_oct_vertex_centroid():
"""The tetrahedral-octahedral honeycomb with each tetrahedron's and
octahedron's centre joined to its corners (``"cell_vertex_centroid"``)."""
return tessellation("tet_oct").edges("cell_vertex_centroid")
def oct_vertex_centroid():
"""The face centres of the octahedra of the tetrahedral-octahedral
honeycomb joined to their three corners (``"facet_vertex_centroid"`` on the
octahedra). The octahedra's own centre-to-corner struts would be
:func:`cubic` at half the pitch."""
return tessellation("tet_oct").edges("facet_vertex_centroid", kinds="oct")
def _periodic_union(*edge_lists):
import numpy as np
segments = [
(np.asarray(start) + 0.5, np.asarray(end) + 0.5)
for edges in edge_lists for start, end in edges
]
return _periodic_graph(segments)
[docs]
def hex_prism_reference_aspect_ratio():
"""Return the U/V cell pitch ratio for regular hexagons, ``2 / sqrt(3)``.
Use it to size the map cells for :func:`hex_prism_edge`, e.g. a U pitch of
``ratio * v_pitch``.
Returns:
float: U pitch divided by V pitch.
"""
return 2.0 / math.sqrt(3.0)
[docs]
def reentrant_hex_prism_reference_aspect_ratio(
angle_degrees=-30.0,
rib_ratio=2.0,
):
"""Return the U/V cell pitch ratio preserving a re-entrant geometry.
Pass the same values given to :func:`reentrant_hex_prism`.
Args:
angle_degrees (float): Inclined-member angle in degrees, strictly
between ``-90`` and ``0``, default ``-30``.
rib_ratio (float): Vertical-to-inclined member length ratio, positive,
default ``2.0``.
Returns:
float: U pitch divided by V pitch.
Raises:
ValueError: If either parameter is out of range.
"""
angle = math.radians(float(angle_degrees))
rib_ratio = float(rib_ratio)
if not -90.0 < float(angle_degrees) < 0.0:
raise ValueError("angle_degrees must lie strictly between -90 and 0")
if not math.isfinite(rib_ratio) or rib_ratio <= 0.0:
raise ValueError("rib_ratio must be positive and finite")
return 4.0 * math.cos(angle) / (rib_ratio - 2.0 * math.sin(angle))
#: Registry mapping a friendly cell name to its edge generator.
UNIT_CELLS = {
"cubic": cubic,
"bcc": bcc,
"fcc": fcc,
"kelvin": kelvin,
"octet": octet,
"diamond": diamond,
"isotruss": isotruss,
"fluorite": fluorite,
"hex_prism_edge": hex_prism_edge,
"reentrant_hex_prism": reentrant_hex_prism,
"rhombic_dodecahedron": rhombic_dodecahedron,
"weaire_phelan": weaire_phelan,
"truncated_cube": truncated_cube,
"column": column,
"columns": columns,
"reentrant_3d": reentrant_3d,
"triangular_prism_edge": triangular_prism_edge,
"rotated_triangular_prism_edge": rotated_triangular_prism_edge,
"rotated_square_prism_edge": rotated_square_prism_edge,
"hex_prism_vertex_centroid": hex_prism_vertex_centroid,
"hex_prism_central_axis": hex_prism_central_axis,
"hex_prism_diamond": hex_prism_diamond,
"hex_prism_laves": hex_prism_laves,
"tet_oct_vertex_centroid": tet_oct_vertex_centroid,
"oct_vertex_centroid": oct_vertex_centroid,
}
_HEX_PRISM_SIZE = (2.0 / math.sqrt(3.0), 1.0, 1.0)
_REFERENCE_CELL_SIZES = {
"hex_prism_edge": _HEX_PRISM_SIZE,
"hex_prism_vertex_centroid": _HEX_PRISM_SIZE,
"hex_prism_central_axis": _HEX_PRISM_SIZE,
"hexagonal_honeycomb": _HEX_PRISM_SIZE,
"hex_prism_diamond": _HEXAGONAL_CELL_SIZE,
"hex_prism_laves": _HEXAGONAL_CELL_SIZE,
"triangular_prism_edge": (1.0, math.sqrt(3.0), 1.0),
"triangular_honeycomb": (1.0, math.sqrt(3.0), 1.0),
"rotated_triangular_prism_edge": (math.sqrt(3.0), 1.0, 1.0),
"rotated_triangular_honeycomb": (math.sqrt(3.0), 1.0, 1.0),
}
[docs]
def reference_cell_size(name, angle_degrees=-30.0, rib_ratio=2.0):
"""Return the U:V:W cell proportions that keep a catalog cell's intended shape.
Give a map cells with these proportions (any common scale) to get regular
hexagons and triangles, ideal bond angles, or the requested re-entrant angle.
Cells without a preferred shape return ``(1, 1, 1)``; for extruded prism and
honeycomb cells the W value is a free choice.
Args:
name (str): A graph cell name from ``LATTICE_NAMES`` or a face cell name
from ``FACE_LATTICE_NAMES``.
angle_degrees (float): For ``reentrant_hex_prism``, ``reentrant_3d`` and
``reentrant_honeycomb``: the inclined-member angle, default ``-30``.
rib_ratio (float): For the same cells: the rib ratio, default ``2.0``.
Returns:
tuple[float, float, float]: Relative U, V, W pitches.
Example:
>>> u, v, w = mm.reference_cell_size("hex_prism_diamond")
>>> pitch = 8.0 / v # 8 mm along V
>>> cell_map = mm.rectangular_cell_map(bounds, cell_size=(u * pitch, v * pitch, w * pitch))
"""
if name in ("reentrant_hex_prism", "reentrant_honeycomb"):
return (reentrant_hex_prism_reference_aspect_ratio(angle_degrees, rib_ratio), 1.0, 1.0)
if name == "reentrant_3d":
ratio = reentrant_hex_prism_reference_aspect_ratio(angle_degrees, rib_ratio)
return (ratio, ratio, 1.0)
if name in _REFERENCE_CELL_SIZES:
return tuple(float(value) for value in _REFERENCE_CELL_SIZES[name])
from .faces import FACE_LATTICE_NAMES
if name in UNIT_CELLS or name in FACE_LATTICE_NAMES:
return (1.0, 1.0, 1.0)
raise ValueError("Unknown unit cell '{}'".format(name))
def unit_cell_edges(name):
"""Return the normalized edge list for a named unit cell."""
if name not in UNIT_CELLS:
raise ValueError("Unknown unit cell '{}'. Available: {}".format(
name, ", ".join(sorted(UNIT_CELLS))))
return UNIT_CELLS[name]()
