## 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
**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`](https://nikolas-claussen.github.io/triangulax/src/halfedge_datastructure.html#hemesh)
class) can be used without modification.
2. A displacement function that takes into account the periodicity.
For a rectangular domain with lengths **L** the displacement vector
between two points **r**1, **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. **L**
dynamically during your simulation (for instance, to simulate growing
tissues). You can even take gradients with respect to **L**.
For sheared domains (e.g. to impose simple shear), we also provide a
displacement function implementing [Lees-Edwards boundary
conditions](https://doi.org/10.1088/0022-3719/5/15/006), 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
``` python
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
``` python
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.*
``` python
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.
``` python
# 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))
``` python
hemesh.bdry_loops # the mesh has no boundary
```
[]
``` python
# 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))
``` python
# 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
``` python
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
``` python
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
``` python
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`](https://nikolas-claussen.github.io/triangulax/src/geometric_quantities_periodic_bcs.html#get_periodic_face_corner_angles)).
See also
:func:[`get_periodic_cotan_weights_per_he`](https://nikolas-claussen.github.io/triangulax/src/geometric_quantities_periodic_bcs.html#get_periodic_cotan_weights_per_he)
for per-half-edge indexing.
------------------------------------------------------------------------
source
### get_periodic_corner_angles
``` python
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
``` python
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`](https://nikolas-claussen.github.io/triangulax/src/geometric_quantities_periodic_bcs.html#get_periodic_corner_angles)
for per-half-edge indexing.
------------------------------------------------------------------------
source
### get_periodic_face_centroids
``` python
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
``` python
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
``` python
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
``` python
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
``` python
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.*
``` python
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
``` python
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
``` python
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
``` python
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
``` python
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.*
``` python
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