Custom particle species

A particle species is a shape together with a set of binding sites on it. Every species is a subtype of ParticleSpecies{D,B}, where D is the spatial dimension and B is the concrete BindingSite type.

Roly ships three general implementations: PolygonParticleSpecies for regular polygons, PolyhedronParticleSpecies for any convex Polyhedron with one site per face, and PatchyParticleSpecies for disks and spheres with sites on the surface. An arbitrary convex solid needs only its corners:

PolyhedronParticleSpecies(Polyhedron([SVector(x, y, z) for x in (-1.0, 1.0)
                                                       for y in (-2.0, 2.0)
                                                       for z in (-3.0, 3.0)]))

Define a new subtype for a shape none of these covers.

The interface

A species must define these methods:

These have defaults you may override:

The graph representation

Each species carries a directed graph recording how its sites relate, ignoring geometry. Roly uses it to detect isomorphic polyforms with nauty. Every automorphism of the graph must be a symmetry of the particle. The usual choice is a directed cycle with one vertex per site, site i on vertex i. cycleencoding builds this for any number of sites. Each BindingSite records its vertices in the vertices field, so BindingSite(pose, color, i:i, ...) puts one site on vertex i.

A site may instead span a contiguous range of vertices, which is how 3D species record the twist of a face: dartencoding gives each face its own directed cycle. A graph built that way must keep each site's vertices together, so that no automorphism can carry part of one site onto part of another.

Every finite rotation group has a dart encoding

Whatever your species looks like, if its rotation group is finite you can encode it as the dart encoding of one of Roly's bodies. Three facts chain together.

1. The finite subgroups of SO(3) are exactly C_n, D_n, T, O and I (the classification of finite rotation groups). RotationGroup names these five families and nothing else.
2. Each is realized as the proper rotation group of a convex polyhedron Roly can build: C_n by Pyramid, D_n by Prism and Antiprism, and T, O, I by the Platonic solids.
3. The automorphism group of a body's dart encoding is that body's rotation group. A rotation is determined by the image of a single dart, and no rotation other than the identity fixes a dart, so the group acts freely on the darts and the graph can record it exactly.

Labeling then carves out any subgroup you need: dartencoding inherits one label per face, and the automorphisms that survive are those preserving the labeling. Taking all labels distinct leaves the trivial group, taking them all equal leaves the full rotation group, and faceorbits gives the labeling for anything in between. So picking an encoding for a custom species is a matter of finding the body with your species' rotation group, not of inventing a graph.

In 2D the question does not arise: the finite subgroups of SO(2) are the C_n, and cycleencoding covers them all.

What a binding site records

Besides its pose and color, a BindingSite carries three numbers that decide how a partner attaches:

Roly reads these to decide how many distinct bonds a pair of sites has.

In 2D both sitesym and stab are 1, so the five-argument BindingSite(pose, color, vertices, touching_tol, alignment_tol) is right. In 3D sitesym is still 1 for a site on a single vertex, but stab need not be, since a rotation about a patch's axis can carry the particle onto itself. Compute it with stabilizerorders and pass BindingSite(pose, color, vertices, tol, tol, sitesym, stab). See Orientation and twists.

A worked example: rectangle

A non-square rectangle with a site at each edge midpoint, exercising the whole interface without special geometry.

using Roly
using Roly: BindingSite, ParticleSpecies, SpeciesAndPose
import Roly: graphrep, nsites, bindingsite,
             bounding_radius, isconvex, overlap
using NautyGraphs
using StaticArrays, LinearAlgebra, Rotations

The struct

The struct holds what Roly needs plus your own geometry, here the two side lengths and the four corners in the particle's frame.

struct Rectangle{F,B<:BindingSite} <: ParticleSpecies{2,B}
    g::NautyDiGraph
    sites::Vector{B}
    corners::Vector{SVector{2,F}}
    width::F
    height::F
    skin::F
end

The parameters {2,B} fix the dimension to 2D and let Roly infer the pose type from the binding sites.

The constructor

It builds the graph, puts one site at each edge midpoint facing outward, and computes the corners.

function Rectangle(width::Real, height::Real; colors=1:4)
    F = float(promote_type(typeof(width), typeof(height)))
    w, h = F(width), F(height)
    tol = sqrt(eps(F)) * max(w, h)

    # Site poses first, so the labeling can be derived from them.
    poses = [
        Pose(SVector{2,F}( w/2, 0),  Angle2d{F}(0)),
        Pose(SVector{2,F}(0,  -h/2), Angle2d{F}(-F(π)/2)),
        Pose(SVector{2,F}(-w/2, 0),  Angle2d{F}(F(π))),
        Pose(SVector{2,F}(0,   h/2), Angle2d{F}(F(π)/2)),
    ]
    # A 2D site has no turn about its in-plane normal, so every `sitesym` is 1.
    labels = siteorbits(poses, ones(Int, 4), collect(colors))

    # Directed 4-cycle, one vertex per binding site; `ranges[i]` is `i:i` here.
    g, ranges = cycleencoding(4; labels)

    # Sites at the edge midpoints, ordered clockwise: right, bottom, left, top.
    sites = [BindingSite(poses[i], colors[i], ranges[i], tol, tol) for i in 1:4]

    corners = [
        SVector{2,F}( w/2,  h/2),
        SVector{2,F}( w/2, -h/2),
        SVector{2,F}(-w/2, -h/2),
        SVector{2,F}(-w/2,  h/2),
    ]

    return Rectangle{F,eltype(sites)}(g, sites, corners, w, h, tol)
end

skin is a small tolerance for distance comparisons, so sites that should touch count as touching despite floating-point noise.

Derive the graph labels with siteorbits rather than writing them out, as every built-in species does. It puts two sites in one orbit exactly when a rotation carries one onto the other and they have the same color, so it reads the symmetry off the geometry you already gave. For this rectangle it returns [1, 2, 1, 2] even when all four colors are equal, because opposite edges are interchangeable and adjacent ones are not. The same call on a square with one color returns [1, 1, 1, 1].

If you do write labels by hand, call check_encoding in your constructor to confirm they match the shape. It requires the labeling to be exactly the symmetry orbits, not merely fine enough to give the right symmetry number, because everything downstream reads a shared label as "a rotation carries one of these sites onto the other".

Interface methods

Four of them just read the struct:

graphrep(ps::Rectangle) = ps.g
nsites(ps::Rectangle) = length(ps.sites)
bindingsite(ps::Rectangle, i::Integer) = ps.sites[i]
bounding_radius(ps::Rectangle) = sqrt(ps.width^2 + ps.height^2) / 2
isconvex(::Rectangle) = true

Base.copy(ps::Rectangle) =
    typeof(ps)(copy(ps.g), copy(ps.sites), copy(ps.corners), ps.width, ps.height, ps.skin)

setcolors! needs no definition. Recoloring re-derives the labeling and the stabilizers along with the colors, and the generic method does all of it, finding the sites in the sites field. A species that stores them elsewhere defines its own method.

Overlap

For convex shapes use the Separating Axis Theorem: two convex bodies are disjoint exactly when their projections onto some axis do not overlap, so a finite set of candidate axes suffices. Only that set differs between dimensions, so sat_overlap takes it from you and does the rest. In 2D it is the edge normals of both polygons; in 3D both solids' face normals plus the cross products of their edge directions.

function overlap(p1::SpeciesAndPose{<:Rectangle}, p2::SpeciesAndPose{<:Rectangle}; kwargs...)
    s1, pose1 = p1
    s2, pose2 = p2
    axes = Iterators.flatten((edgenormals(s1.corners, pose1), edgenormals(s2.corners, pose2)))
    return sat_overlap(axes, s1.corners, pose1, s2.corners, pose2, s1.skin + s2.skin)
end

Axes need not be normalized or even nonzero, since sat_overlap scales each one and skips degenerate ones.

Using the new species

It now works like any built-in one:

rect = Rectangle(2.0, 1.0)
bonds = [1 1 1 3;  # right edge binds to left edge
         1 2 1 4]  # bottom edge binds to top edge
rules = BindingRules(bonds, rect)

polys = polygen(rules; maxsize=6)