triangulax

`adjacency`: Adjacency-based operators on half-edge meshes

The HeMesh data structure enables efficient mesh “traversal”. Using such traversals, one can express many adjacency-based operators, for example:

Sum over all half-edges “incoming” to a vertex (special case: count the incoming edges, i.e., compute the coordination number)
Compute the finite-element gradient of a function defined on vertices
Find the minimum of a per-vertex function across the neighbors of a vertex

These operations can be done efficiently using a “gather/scatter” approach, see jax.numpy.ndarray.at. There is no need to explicitly instantiate a matrix for the operators.

All operators defined in this notebook depend only on the mesh topology, not the geometry (vertex/face positions). All functions are compatible with jax.jit and jax.vmap.

Note on boundary half-edges: On meshes with boundary, some half-edges are boundary half-edges (heface == -1). These are included in scatter-add operations to vertices (e.g. sum_he_to_vertex_incoming), which is correct for most use cases (e.g., coordination number counts all neighbors). If you need to exclude boundary contributions, mask with hemesh.is_bdry_he before calling.

from triangulax.triangular import TriMesh

# load test data

mesh = TriMesh.read_obj("../test_meshes/disk.obj", dim=2)
hemesh = msh.HeMesh.from_triangles(mesh.vertices.shape[0], mesh.faces)
geommesh = msh.GeomMesh(mesh.vertices, mesh.face_positions)

mesh_3d = TriMesh.read_obj("../test_meshes/disk.obj", dim=3)
geommesh_3d = msh.GeomMesh(mesh_3d.vertices, mesh_3d.face_positions)

Warning: readOBJ() ignored non-comment line 3:
  o flat_tri_ecmc
Warning: readOBJ() ignored non-comment line 3:
  o flat_tri_ecmc

Discrete (exterior) derivative

On a triangular mesh, there are two natural “derivatives”: for a per-vertex field, the difference across half-edges, and for a per-half-edge field, the circulation around a face (this is the basis of discrete exterior calculus).

source

get_exterior_circulation


def get_exterior_circulation(
    hemesh:HeMesh, he_field:Float[Array, 'n_hes ...']
)->Float[Array, 'n_faces ...']:

Discrete exterior derivative: circulation of a half-edge field around each face.

source

get_exterior_gradient


def get_exterior_gradient(
    hemesh:HeMesh, v_field:Float[Array, 'n_vertices ...']
)->Float[Array, 'n_hes ...']:

Discrete exterior derivative: difference of a vertex field across each half-edge.

# define a random scalar field on vertices and compute its gradient on halfedges
v_field = jax.random.uniform(jax.random.PRNGKey(0), (hemesh.n_vertices,))
he_gradient = get_exterior_gradient(hemesh, v_field)
f_circulation = get_exterior_circulation(hemesh, he_gradient)

assert he_gradient.shape == (hemesh.n_hes,)
assert f_circulation.shape == (hemesh.n_faces,)
# d^2 = 0: the circulation of a gradient vanishes
assert jnp.allclose(f_circulation, 0, atol=1e-12)
# the exterior gradient is antisymmetric under the twin map
assert jnp.allclose(he_gradient + he_gradient[hemesh.twin], 0, atol=1e-12)

Summing over adjacent mesh elements

A second important class of operation is summing over adjacent mesh elements. For example, to get the coordination number of a vertex, one sums the value \(1\) over all incoming half-edges. For computing things like cell areas, it’s also useful to sum over half-edges opposite to a vertex.

source

sum_he_to_vertex_opposite


def sum_he_to_vertex_opposite(
    hemesh:HeMesh, he_field:Float[Array, 'n_hes ...']
)->Float[Array, 'n_vertices ...']:

Sum a half-edge field onto opposite vertices (the vertex across the face from the half-edge).

Warning: includes boundary half-edges, whose “opposite vertex” may not be meaningful.

source

sum_he_to_vertex_outgoing


def sum_he_to_vertex_outgoing(
    hemesh:HeMesh, he_field:Float[Array, 'n_hes ...']
)->Float[Array, 'n_vertices ...']:

Sum a half-edge field onto origin vertices (scatter-add over outgoing half-edges).

Includes boundary half-edges. To exclude them, zero out the field at boundary half-edges before calling.

source

sum_he_to_vertex_incoming


def sum_he_to_vertex_incoming(
    hemesh:HeMesh, he_field:Float[Array, 'n_hes ...']
)->Float[Array, 'n_vertices ...']:

Sum a half-edge field onto destination vertices (scatter-add over incoming half-edges).

Includes boundary half-edges. To exclude them, zero out the field at boundary half-edges before calling.

# test: sum_he_to_vertex_outgoing should equal sum_he_to_vertex_incoming of the twin field
he_field = jax.random.normal(jax.random.PRNGKey(7), (hemesh.n_hes,))
assert jnp.allclose(sum_he_to_vertex_outgoing(hemesh, he_field),
                     sum_he_to_vertex_incoming(hemesh, he_field[hemesh.twin]))
print("sum_he_to_vertex_outgoing test passed")

sum_he_to_vertex_outgoing test passed

source

sum_face_to_he


def sum_face_to_he(
    hemesh:HeMesh, f_field:Float[Array, 'n_faces ...']
)->Float[Array, 'n_hes ...']:

Sum face-field to half-edges. Each half-edge receives the sum of its face and its twin’s face.

Boundary half-edges (heface == -1) contribute zero from the boundary side.

source

sum_he_to_face


def sum_he_to_face(
    hemesh:HeMesh, he_field:Float[Array, 'n_hes ...']
)->Float[Array, 'n_faces ...']:

Sum over the three half-edges of each face. Same as get_exterior_circulation.

# test sum_face_to_he: each interior half-edge should get f[face] + f[twin's face]
f_field = jax.random.normal(jax.random.PRNGKey(42), (hemesh.n_faces,))
he_summed = sum_face_to_he(hemesh, f_field)

# verify on an interior half-edge
interior_he = jnp.where(~hemesh.is_bdry_he, size=1)[0][0]
expected = f_field[hemesh.heface[interior_he]] + f_field[hemesh.heface[hemesh.twin[interior_he]]]
assert jnp.allclose(he_summed[interior_he], expected), "sum_face_to_he mismatch on interior he"

# boundary half-edges: the boundary side contributes zero
if hemesh.is_bdry_he.any():
    bdry_he = jnp.where(hemesh.is_bdry_he, size=1)[0][0]
    # boundary he has heface == -1, so its contribution should be 0
    twin_face = hemesh.heface[hemesh.twin[bdry_he]]
    expected_bdry = f_field[twin_face]  # only the twin's face contributes
    assert jnp.allclose(he_summed[bdry_he], expected_bdry), "sum_face_to_he mismatch on boundary he"
    print("boundary he test passed")

print("sum_face_to_he test passed, shape:", he_summed.shape)

boundary he test passed
sum_face_to_he test passed, shape: (708,)

source

average_face_to_vertex


def average_face_to_vertex(
    hemesh:HeMesh, f_field:Float[Array, 'n_faces ...']
)->Float[Array, 'n_vertices ...']:

Average face-field to vertices (uniform weights, i.e., divided by number of incident faces).

source

sum_face_to_vertex


def sum_face_to_vertex(
    hemesh:HeMesh, f_field:Float[Array, 'n_faces ...']
)->Float[Array, 'n_vertices ...']:

Sum face-field to vertices. Each vertex receives the sum over its incident faces.

source

average_vertex_to_face


def average_vertex_to_face(
    hemesh:HeMesh, v_field:Float[Array, 'n_vertices ...']
)->Float[Array, 'n_faces ...']:

Average vertex-field to faces (uniform 1/3 weights).

source

sum_vertex_to_face


def sum_vertex_to_face(
    hemesh:HeMesh, v_field:Float[Array, 'n_vertices ...']
)->Float[Array, 'n_faces ...']:

Sum vertex-field to faces. Sums over the three vertices of each face.

# tests vs libigl
key = jax.random.PRNGKey(123)

u_v = jax.random.normal(key, (hemesh.n_vertices,))
faces_avg_jax = average_vertex_to_face(hemesh, u_v)
faces_avg_igl = igl.average_onto_faces(np.asarray(hemesh.faces), np.asarray(u_v))

rel_err_faces = jnp.linalg.norm(faces_avg_jax - faces_avg_igl) / jnp.linalg.norm(faces_avg_igl)
print("vertex->face rel. error:", rel_err_faces)

vertex->face rel. error: 0.0

u_f = jax.random.normal(key, (hemesh.n_faces,))
verts_avg_jax = average_face_to_vertex(hemesh, u_f)
verts_avg_igl = igl.average_onto_vertices(mesh.vertices, np.asarray(hemesh.faces), np.asarray(u_f))
rel_err_verts = jnp.linalg.norm(verts_avg_jax-verts_avg_igl) / jnp.linalg.norm(verts_avg_igl)

print("face->vertex rel. error:", rel_err_verts)

face->vertex rel. error: 8.339340577730768e-17

# also works for vector fields
u_f = jax.random.normal(key, (hemesh.n_faces, 10))
verts_avg_jax = average_face_to_vertex(hemesh, u_f)
verts_avg_jax.shape

(131, 10)

source

get_coordination_number


def get_coordination_number(
    hemesh:HeMesh
)->Float[Array, 'n_vertices']:

Coordination number (vertex degree): number of neighboring vertices.

Includes boundary half-edges, so boundary vertices count all their neighbors.

get_coordination_number(hemesh).mean()

Array(5.40458015, dtype=float64)

Uniform/graph Laplacian

The uniform (or graph) Laplacian is \(L = D - A\), where \(D\) is the diagonal degree matrix and \(A\) is the adjacency matrix. It is positive semi-definite, with a zero eigenvalue for constant fields. This purely topological Laplacian ignores edge lengths and angles; for the geometry-aware cotangent Laplacian, see linops.

source

get_uniform_laplacian


def get_uniform_laplacian(
    hemesh:HeMesh
)->BCOO:

Assemble the uniform Laplacian as a sparse BCOO matrix. Non-normalized, positive semi-definite.

source

compute_uniform_laplacian


def compute_uniform_laplacian(
    hemesh:HeMesh, v_field:Float[Array, 'n_vertices ...']
)->Float[Array, 'n_vertices ...']:

Apply the uniform Laplacian to a vertex field. Non-normalized, positive semi-definite.

# test that the matrix and function versions are equivalent

laplace_mat = get_uniform_laplacian(hemesh)

assert jnp.allclose(laplace_mat @ v_field, compute_uniform_laplacian(hemesh, v_field), atol=1e-12)
# the uniform Laplacian is symmetric, positive semi-definite, with vanishing row sums
dense = laplace_mat.todense()
assert jnp.allclose(dense, dense.T, atol=1e-12)
assert jnp.allclose(dense.sum(axis=1), 0, atol=1e-12)
assert jnp.linalg.eigvalsh(dense).min() > -1e-10
assert jnp.dot(laplace_mat @ v_field, v_field) > 0
# and equals D - A from igl
deg = jnp.diag(get_coordination_number(hemesh))
A_igl = jnp.array(igl.adjacency_matrix(np.asarray(hemesh.faces, dtype=np.int64)).todense())
assert jnp.allclose(dense, deg - A_igl, atol=1e-12)

adjacency: Adjacency-based operators on half-edge meshes

Discrete (exterior) derivative

get_exterior_circulation

get_exterior_gradient

Summing over adjacent mesh elements

sum_he_to_vertex_opposite

sum_he_to_vertex_outgoing

sum_he_to_vertex_incoming

sum_face_to_he

sum_he_to_face

average_face_to_vertex

sum_face_to_vertex

average_vertex_to_face

sum_vertex_to_face

get_coordination_number

Uniform/graph Laplacian

get_uniform_laplacian

compute_uniform_laplacian

`adjacency`: Adjacency-based operators on half-edge meshes