Geometric quantities in periodic boundary conditions

In many 2D biophysical simulations (see tutorial 3 on “vertex models”), it is convenient to work with periodic boundary conditions, i.e. simulate cells in a box with lengths \(\mathbf{L} = [L_x, L_y]\) where opposite sides are identified. This module contains tools for mesh geometry in periodic boundary conditions.

In triangulax, mesh connectivity and geometry are decoupled, so periodic boundary conditions are easy to implement. One needs two ingredients:

  1. A triangulation whose connectivity has the desired periodicity (e.g. a triangulation of a torus).

There are different ways to generate a triangulation of the torus - one example is included in test_meshes/torus_2d.obj. Note: edge flips/T1 on this mesh will preserve its topology, no further bookkeeping is needed. In particular, one does not need to keep track of any boundary vertices since there is no boundary. All triangulax tools related to the mesh connectivity (like the HeMesh class) can be used without modification.

  1. A displacement function that takes into account the periodicity.

For a rectangular domain with lengths \(\mathbf{L}\) the displacement vector between two points \(\mathbf{r}_1, \mathbf{r}_2\) can be computed as \[\mathbf{d} = \mathbf{r}_2 - \mathbf{r}_1 - \mathrm{round}\left(\frac{\mathbf{r}_2 - \mathbf{r}_1}{L} \right)\mathbf{L}\]

This always gives the shortest displacement vector between the two points (minimum-image convention). With this displacement function, one can compute edge lengths, angles, Voronoi duals, etc. Note: nothing prevents you from changing the shape of your domain, i.e. \(\mathbf{L}\) dynamically during your simulation (for instance, to simulate growing tissues). You can even take gradients with respect to \(\mathbf{L}\).

For sheared domains (e.g. to impose simple shear), we also provide a displacement function implementing Lees-Edwards boundary conditions, which add an x-shift proportional to a shear factor \(s\) when wrapping in the y-direction.

Note: throughout this module, the displacement_fn argument is a callable (r_1, r_2) -> d that returns a 2D displacement vector, not a scalar distance.

source

displacement_periodic_twisted


def displacement_periodic_twisted(
    r_1:Float[Array, '2'], r_2:Float[Array, '2'], L:Float[Array, '2'], # Box lengths [L_x, L_y].
    s:float, # Shear factor: wrapping in y shifts x by ``s * L_x``.
)->Float[Array, '2']:

Return the minimum-image displacement on a sheared periodic torus (Lees-Edwards BCs).

source

displacement_periodic


def displacement_periodic(
    r_1:Float[Array, '2'], r_2:Float[Array, '2'], L:Float[Array, '2'], # Box lengths [L_x, L_y].
)->Float[Array, '2']:

Return the minimum-image displacement on a rectangular torus.

L = jnp.array([2.0, 3.0])

assert jnp.allclose(
    displacement_periodic(
        jnp.array([0.1, 0.2]),
        jnp.array([1.9, 2.8]),
        L,
    ),
    jnp.array([-0.2, -0.4]),
)

assert jnp.allclose(
    displacement_periodic(
        jnp.array([1.9, 2.8]),
        jnp.array([0.1, 0.2]),
        L,
    ),
    jnp.array([0.2, 0.4]),
)

assert jnp.allclose(
    displacement_periodic_twisted(
        jnp.array([0.0, 0.0]),
        jnp.array([0.9, 1.1]),
        jnp.array([1.0, 1.0]),
        0.25,
    ),
    jnp.array([0.15, 0.1]),
)

assert jnp.allclose(
    displacement_periodic_twisted(
        jnp.array([0.9, 0.9]),
        jnp.array([0.1, -0.1]),
        jnp.array([1.0, 1.0]),
        0.25,
    ),
    jnp.array([-0.05, 0.0]),
)

Edge lengths, angles, and areas

Geometric quantities like angles and triangle areas can be computed using the displacement functions defined above. To do this, we can use the “intrinsic” geometry functions from the trigonometry module to compute everything in terms of edge lengths.

# this is an example of a 2d mesh with "periodic boundary conditions" - topologically, it is just a torus.
# Thus, the faces that wrap around the boundary.

trimesh = TriMesh.read_obj("../test_meshes/torus_2d.obj", dim=2)
hemesh = HeMesh.from_triangles(trimesh.vertices.shape[0], trimesh.faces)

vertices = trimesh.vertices

plt.triplot(vertices[:,0], vertices[:,1], trimesh.faces, lw=0.5, color="k")
plt.axis("equal")
(np.float64(-0.028125350000000004),
 np.float64(1.04895835),
 np.float64(0.03749965),
 np.float64(1.04583335))

hemesh.bdry_loops  # the mesh has no boundary
[]
# to compute angles, we use the intrinsic geometry of the mesh

L = jnp.array([1.0, 1.0]) # box lengths
my_displacement_fn = lambda r_1, r_2: displacement_periodic(r_1, r_2, L)

edge_lengths = jax.vmap(my_displacement_fn)(vertices[hemesh.orig], vertices[hemesh.dest])
edge_lengths_nonperiodic = jnp.linalg.norm(vertices[hemesh.orig] - vertices[hemesh.dest], axis=1)

edge_lengths.max(), edge_lengths_nonperiodic.max()
(Array(0.083334, dtype=float64), Array(1.34128535, dtype=float64))
# let's group together the edge lengths for each face, then compute the area

la, lb, lc = (edge_lengths[hemesh.prv[hemesh.face_incident]],
              edge_lengths[hemesh.face_incident],
              edge_lengths[hemesh.nxt[hemesh.face_incident]],)

angles = jax.vmap(trig.get_angles_from_lengths)(la, lb, lc)
areas = jax.vmap(trig.get_triangle_area_from_lengths)(la, lb, lc)

source

get_periodic_cotan_weights_per_edge


def get_periodic_cotan_weights_per_edge(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_hes']:

Average cotangent weight per edge, using a periodic displacement function.

Returns (cot_he + cot_twin) / 2 (same convention as :func:geometry.get_cotan_weights_per_edge).

source

get_periodic_cotan_weights_per_he


def get_periodic_cotan_weights_per_he(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_hes']:

Cotangent of the angle opposite to each half-edge, using a periodic displacement function.

Returns a per-half-edge array (same convention as :func:geometry.get_cotan_weights_per_he).

source

get_periodic_face_corner_cotangents


def get_periodic_face_corner_cotangents(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_faces 3']:

Compute corner cotangents for each face using a periodic displacement function.

Returns an (n_faces, 3) array (same layout as :func:get_periodic_face_corner_angles). See also :func:get_periodic_cotan_weights_per_he for per-half-edge indexing.

source

get_periodic_corner_angles


def get_periodic_corner_angles(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_hes']:

Compute the angle opposite to each half-edge, using a periodic displacement function.

Returns a per-half-edge array (same convention as :func:geometry.get_corner_angles).

source

get_periodic_face_corner_angles


def get_periodic_face_corner_angles(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_faces 3']:

Compute corner angles for each face using a periodic displacement function.

Returns an (n_faces, 3) array where columns correspond to the three corners of each face in half-edge order (prv, face_incident, nxt). See also :func:get_periodic_corner_angles for per-half-edge indexing.

source

get_periodic_face_centroids


def get_periodic_face_centroids(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_faces 2']:

Compute face centroids (barycenters) using a periodic displacement function.

Returns positions relative to vertex a of each face, not wrapped to the periodic box.

source

get_periodic_barycentric_cell_areas


def get_periodic_barycentric_cell_areas(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_vertices']:

Get area of barycentric dual cell around each vertex, using a periodic displacement function.

Defined as 1/3 * sum of adjacent triangle areas.

source

get_periodic_area


def get_periodic_area(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, '']:

Total mesh area using a periodic displacement function.

source

get_periodic_triangle_areas


def get_periodic_triangle_areas(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_faces']:

Compute triangle areas using a periodic displacement function.

source

get_periodic_he_lengths


def get_periodic_he_lengths(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_hes']:

Get lengths of half-edges using a periodic displacement function.

periodic_displacement = lambda r_1, r_2: displacement_periodic(r_1, r_2, L)

periodic_areas = get_periodic_triangle_areas(vertices, hemesh, periodic_displacement)
periodic_face_angles = get_periodic_face_corner_angles(vertices, hemesh, periodic_displacement)
periodic_face_cotangents = get_periodic_face_corner_cotangents(vertices, hemesh, periodic_displacement)
periodic_he_angles = get_periodic_corner_angles(vertices, hemesh, periodic_displacement)
periodic_he_cotangents = get_periodic_cotan_weights_per_he(vertices, hemesh, periodic_displacement)

manual_edge_lengths = jnp.linalg.norm(
    jax.vmap(periodic_displacement)(vertices[hemesh.orig], vertices[hemesh.dest]),
    axis=-1,)
la_test, lb_test, lc_test = (manual_edge_lengths[hemesh.prv[hemesh.face_incident]],
                             manual_edge_lengths[hemesh.face_incident],
                             manual_edge_lengths[hemesh.nxt[hemesh.face_incident]],)

# shape checks
assert periodic_areas.shape == (hemesh.n_faces,)
assert periodic_face_angles.shape == (hemesh.n_faces, 3)
assert periodic_face_cotangents.shape == (hemesh.n_faces, 3)
assert periodic_he_angles.shape == (hemesh.n_hes,)
assert periodic_he_cotangents.shape == (hemesh.n_hes,)

# face-indexed consistency
assert jnp.allclose(
    periodic_areas,
    jax.vmap(trig.get_triangle_area_from_lengths)(la_test, lb_test, lc_test),)
assert jnp.allclose(
    periodic_face_angles,
    jax.vmap(trig.get_angles_from_lengths)(la_test, lb_test, lc_test),)
assert jnp.allclose(
    periodic_face_cotangents,
    jax.vmap(trig.get_cotangents_from_lengths)(la_test, lb_test, lc_test),)

# per-halfedge vs per-face consistency
assert jnp.allclose(periodic_he_angles[hemesh.face_incident], periodic_face_angles[:, 1])
assert jnp.allclose(periodic_he_cotangents[hemesh.face_incident], periodic_face_cotangents[:, 1])

# total area
assert jnp.allclose(get_periodic_area(vertices, hemesh, periodic_displacement), periodic_areas.sum())

# barycentric cell areas
bary_areas = get_periodic_barycentric_cell_areas(vertices, hemesh, periodic_displacement)
assert bary_areas.shape == (hemesh.n_vertices,)
assert jnp.allclose(bary_areas.sum(), periodic_areas.sum())

# he lengths
he_lens = get_periodic_he_lengths(vertices, hemesh, periodic_displacement)
assert he_lens.shape == (hemesh.n_hes,)
assert jnp.allclose(he_lens, manual_edge_lengths)

# face centroids
centroids = get_periodic_face_centroids(vertices, hemesh, periodic_displacement)
assert centroids.shape == (hemesh.n_faces, 2)

# zero-twist equivalence
zero_twist = lambda r_1, r_2: displacement_periodic_twisted(r_1, r_2, L, 0.0)
assert jnp.allclose(periodic_areas, get_periodic_triangle_areas(vertices, hemesh, zero_twist))
assert jnp.allclose(periodic_face_angles, get_periodic_face_corner_angles(vertices, hemesh, zero_twist))
assert jnp.allclose(periodic_face_cotangents, get_periodic_face_corner_cotangents(vertices, hemesh, zero_twist))
assert jnp.allclose(periodic_he_angles, get_periodic_corner_angles(vertices, hemesh, zero_twist))

print("periodic geometry tests passed")
periodic geometry tests passed

Voronoi dual

We can also compute the Voronoi dual for periodic boundary conditions. Since vertex coordinates cannot be subtracted directly across the periodic boundary, we compute the circumcenter in intrinsic (edge-length) barycentric coordinates via trig.get_circumcenter_from_lengths, then map back to Euclidean positions using displacement vectors that respect the periodicity.

source

get_periodic_voronoi_perimeters


def get_periodic_voronoi_perimeters(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_vertices']:

Compute Voronoi cell perimeters by summing dual edge lengths per vertex.

source

get_periodic_dual_he_length


def get_periodic_dual_he_length(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_hes']:

Voronoi dual edge lengths computed from cotangent weights.

Equivalent to cotan_weights_per_edge * he_length.

source

get_periodic_voronoi_face_positions


def get_periodic_voronoi_face_positions(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_faces 2']:

Compute periodic Voronoi dual positions (circumcenters) from intrinsic barycentric coordinates.

source

get_periodic_voronoi_areas


def get_periodic_voronoi_areas(
    vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh, displacement_fn:Callable
)->Float[Array, 'n_vertices']:

Compute Voronoi cell areas from periodic edge lengths and cotangent weights.

periodic_voronoi_areas = get_periodic_voronoi_areas(vertices, hemesh, periodic_displacement)
periodic_face_positions = get_periodic_voronoi_face_positions(vertices, hemesh, periodic_displacement)
periodic_he_lengths = get_periodic_he_lengths(vertices, hemesh, periodic_displacement)
periodic_cotan_weights = get_periodic_cotan_weights_per_edge(vertices, hemesh, periodic_displacement)
periodic_dual_lengths = get_periodic_dual_he_length(vertices, hemesh, periodic_displacement)
periodic_perimeters = get_periodic_voronoi_perimeters(vertices, hemesh, periodic_displacement)

assert periodic_voronoi_areas.shape == (hemesh.n_vertices,)
assert periodic_face_positions.shape == (hemesh.n_faces, 2)
assert periodic_dual_lengths.shape == (hemesh.n_hes,)
assert periodic_perimeters.shape == (hemesh.n_vertices,)

# voronoi areas from cotan weights
assert jnp.allclose(
    periodic_voronoi_areas,
    adj.sum_he_to_vertex_incoming(
        hemesh,
        periodic_cotan_weights * periodic_he_lengths**2 / 4,),)
assert jnp.allclose(periodic_voronoi_areas.sum(), periodic_areas.sum())

# dual edge lengths = cotan_weights_per_edge * he_length
assert jnp.allclose(periodic_dual_lengths, periodic_cotan_weights * periodic_he_lengths)

# perimeters = sum of dual edge lengths
assert jnp.allclose(periodic_perimeters,
                     adj.sum_he_to_vertex_incoming(hemesh, periodic_dual_lengths))

# circumradius check: circumcenter equidistant from all three vertices
face_hes = hemesh.face_incident
a = vertices[hemesh.orig[face_hes]]
ab = jax.vmap(periodic_displacement)(vertices[hemesh.orig[face_hes]], vertices[hemesh.dest[face_hes]])
bc = jax.vmap(periodic_displacement)(vertices[hemesh.orig[hemesh.nxt[face_hes]]],
                                     vertices[hemesh.dest[hemesh.nxt[face_hes]]])
b = a + ab
c = b + bc
radii_a = jnp.linalg.norm(periodic_face_positions - a, axis=1)
radii_b = jnp.linalg.norm(periodic_face_positions - b, axis=1)
radii_c = jnp.linalg.norm(periodic_face_positions - c, axis=1)
assert jnp.allclose(radii_a, radii_b)
assert jnp.allclose(radii_b, radii_c)

# zero-twist equivalence
zero_twist = lambda r_1, r_2: displacement_periodic_twisted(r_1, r_2, L, 0.0)
assert jnp.allclose(periodic_voronoi_areas, get_periodic_voronoi_areas(vertices, hemesh, zero_twist))
assert jnp.allclose(periodic_face_positions, get_periodic_voronoi_face_positions(vertices, hemesh, zero_twist))
assert jnp.allclose(periodic_perimeters, get_periodic_voronoi_perimeters(vertices, hemesh, zero_twist))

print("periodic Voronoi tests passed")
periodic Voronoi tests passed