import numpy as np
import igl

Adjacency-based operators on half-edge meshes

Using the HeMesh data structure, we can efficiently “traverse” our mesh. 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")
hemesh = msh.HeMesh.from_triangles(mesh.vertices.shape[0], mesh.faces)
geommesh = msh.GeomMesh(*hemesh.n_items, mesh.vertices, mesh.face_positions)

mesh_3d = TriMesh.read_obj("../test_meshes/disk.obj", dim=3)
geommesh_3d = msh.GeomMesh(*hemesh.n_items, 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)

hemesh, he_gradient.shape, f_circulation.shape, jnp.allclose(f_circulation, 0 )

(HeMesh(N_V=131, N_HE=708, N_F=224), (708,), (224,), Array(True, dtype=bool))

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, you want to sum 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)

jnp.allclose(laplace_mat @ v_field, compute_uniform_laplacian(hemesh, v_field)), jnp.dot(laplace_mat@v_field, v_field) > 0

(Array(True, dtype=bool), Array(True, dtype=bool))