## Mesh geometry
Using the half-edge mesh and the adjacency-like operators it defines, we
can compute all sorts of geometric quantities of interest: edge lengths,
cell areas, curvature in 3d, etc.
#### Discrete exterior calculus and discrete Hodge star
Note: triangle area, cell area, edge length, and dual edge lengths are
what’s required for [discrete exterior
calculus](https://www.cs.cmu.edu/~kmcrane/Projects/DDG/paper.pdf).
``` python
from triangulax.triangular import TriMesh
```
``` python
# 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
### Basic mesh geometry
------------------------------------------------------------------------
source
### get_he_length
``` python
def get_he_length(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_hes']:
```
*Get lengths of half-edges (triangulation/primal edges).*
------------------------------------------------------------------------
source
### get_dihedral_angles
``` python
def get_dihedral_angles(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_hes ...']:
```
*Get dihedral angles (angle between adjacent face normals).*
------------------------------------------------------------------------
source
### get_vertex_normals
``` python
def get_vertex_normals(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_faces ...']:
```
*Compute per-vertex normals by averaging over faces*
------------------------------------------------------------------------
source
### get_triangle_normals
``` python
def get_triangle_normals(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_faces ...']:
```
*Compute normals. In 2d, this just returns +/-1.*
------------------------------------------------------------------------
source
### get_oriented_triangle_areas
``` python
def get_oriented_triangle_areas(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_faces ...']:
```
*Compute oriented triangle areas in a mesh. In 3d, this is a vector.*
------------------------------------------------------------------------
source
### get_triangle_areas
``` python
def get_triangle_areas(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_faces ...']:
```
*Compute triangle areas in a mesh.*
------------------------------------------------------------------------
source
### set_voronoi_face_positions
``` python
def set_voronoi_face_positions(
geommesh:GeomMesh, hemesh:HeMesh
)->GeomMesh:
```
*Set face positions of geommesh to the circumcenters of the faces
defined by hemesh.*
------------------------------------------------------------------------
source
### get_voronoi_face_positions
``` python
def get_voronoi_face_positions(
vertices:Float[Array, 'n_vertices 2'], hemesh:HeMesh
)->Float[Array, 'n_faces 2']:
```
*Get face positions of geommesh to the circumcenters of the faces
defined by hemesh.*
------------------------------------------------------------------------
source
### get_oriented_dual_he_length
``` python
def get_oriented_dual_he_length(
vertices:Float[Array, 'n_vertices 2'], face_positions:Float[Array, 'n_faces 2'], hemesh:HeMesh
)->Float[Array, 'n_hes']:
```
*Compute lengths of dual edges. Boundary dual edges get length 1.
Negative sign = flipped edge.*
------------------------------------------------------------------------
source
### get_dual_he_length
``` python
def get_dual_he_length(
face_positions:Float[Array, 'n_faces dim'], hemesh:HeMesh
)->Float[Array, 'n_hes']:
```
*Get lengths of dual/cell half-edges.*
``` python
a = get_dual_he_length(mesh.face_positions, hemesh)
b = get_oriented_dual_he_length(mesh.vertices, mesh.face_positions, hemesh)
jnp.allclose(a[~hemesh.is_bdry_edge], jnp.abs(b)[~hemesh.is_bdry_edge])
```
Array(True, dtype=bool)
``` python
# edges and dual edges should be orthogonal since we are using circumcenters
face_positions = get_voronoi_face_positions(mesh.vertices, hemesh)
edges = mesh.vertices[hemesh.orig]-mesh.vertices[hemesh.dest]
dual_edges = (face_positions[hemesh.heface]-face_positions[hemesh.heface[hemesh.twin]])
jnp.allclose(jnp.einsum('vi,vi->v', edges[~hemesh.is_bdry_edge], dual_edges[~hemesh.is_bdry_edge]), 0)
```
Array(True, dtype=bool)
``` python
# computing the signed edge length shows that there are some "flipped" edges.
dual_length = get_oriented_dual_he_length(mesh.vertices, face_positions, hemesh)
jnp.where((dual_length < -0.0) & ~hemesh.is_bdry_edge )[0]
```
Array([ 9, 185, 191, 335, 363, 539, 545, 689], dtype=int64)
------------------------------------------------------------------------
source
### get_voronoi_edge_lengths
``` python
def get_voronoi_edge_lengths(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_hes']:
```
*Computed directly from angles. Accurate in any dimension*
------------------------------------------------------------------------
source
### get_cotan_weights_per_egde
``` python
def get_cotan_weights_per_egde(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_hes']:
```
*Average of cotangent of angles opposite to edge.*
------------------------------------------------------------------------
source
### get_cotan_weights_per_he
``` python
def get_cotan_weights_per_he(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_hes']:
```
*Cotangent of angle opposite to half-edge*
------------------------------------------------------------------------
source
### get_angle_sum
``` python
def get_angle_sum(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_vertices']:
```
*Angle sum around vertices. 2*pi-angle sum measures Gaussian curvature\*
------------------------------------------------------------------------
source
### get_corner_angles
``` python
def get_corner_angles(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_hes']:
```
*Get angles in mesh corners (opposite to half-edges).*
``` python
angles = get_corner_angles(mesh.vertices, hemesh)
np.allclose(get_angle_sum(mesh.vertices, hemesh)[~hemesh.is_bdry], 2*jnp.pi) # mesh is not curved
```
True
``` python
jnp.allclose(1/jnp.tan(angles), get_cotan_weights_per_he(mesh.vertices, hemesh))
```
Array(True, dtype=bool)
``` python
# we can either compute the Voronoi-length of a dual edge directly, or from the face positions
voronoi_edge_lengths = get_voronoi_edge_lengths(mesh.vertices, hemesh)
dual_edge_length = get_oriented_dual_he_length(mesh.vertices, mesh.face_positions, hemesh)
```
``` python
jnp.allclose(voronoi_edge_lengths[~hemesh.is_bdry_edge], dual_edge_length[~hemesh.is_bdry_edge])
```
Array(True, dtype=bool)
### Computing cell areas, perimeters, etc via corners
To compute, for instance, the cell area using the shoelace formula, you
need to iterate around the faces adjacent to a vertex. This is not
straightforward to vectorize because the number of adjacent faces per
vertex can vary (there can be 5-, 6-, 7-sided cells etc.). One way to
solve this is a scheme in which the lists of adjacent faces are “padded”
in some manner, so that they are all the same length. This is
cumbersome. Instead, let us split all “cell-based” quantities into
contributions from “corners”, i.e., half-edges, like this:
Source:
[CGAL](https://doc.cgal.org/latest/Weights/group__PkgWeightsRefVoronoiRegionWeights.html)
To compute the total area, we can sum over all half-edges (*r*, *p*)
opposite to a vertex *q*.
This approach can also compute cell perimeter, …
------------------------------------------------------------------------
source
### get_voronoi_areas
``` python
def get_voronoi_areas(
vertices:Float[Array, 'n_vertices dim'], hemesh:HeMesh
)->Float[Array, 'n_vertices']:
```
*Compute Voronoi area for each vertex.*
------------------------------------------------------------------------
source
### get_cell_areas_traversal
``` python
def get_cell_areas_traversal(
geommesh:GeomMesh, hemesh:HeMesh
)->Float[Array, 'n_vertices']:
```
*Compute areas of cells by mesh traversal (don’t use for simulation,
inefficient).*
Boundary vertices get area 0.
``` python
## Let's use the adjacency matrix to compute the area of all cells. First, compute all corner areas
a, b, c = (hemesh.dest[hemesh.nxt], hemesh.dest[hemesh.prv], hemesh.dest)
corner_areas = jax.vmap(trig.get_voronoi_corner_area)(geommesh.vertices[a], geommesh.vertices[b], geommesh.vertices[c])
cell_areas_corner = adj.sum_he_to_vertex_opposite(hemesh, corner_areas)
cell_areas_corner = cell_areas_corner.at[hemesh.is_bdry].set(0)
```
``` python
# for comparison, compute the areas by mesh traversal
cell_areas_iterative = get_cell_areas_traversal(geommesh, hemesh)
np.abs(cell_areas_iterative-cell_areas_corner).max() # works!
```
np.float64(5.551115123125783e-17)