415 lines
12 KiB
Rust
415 lines
12 KiB
Rust
//! Tetrahedral mesh data structures.
|
||
|
||
use nalgebra::Point3;
|
||
use serde::{Deserialize, Serialize};
|
||
use std::collections::{HashMap, HashSet};
|
||
|
||
/// A vertex in the tetrahedral mesh.
|
||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||
pub struct Vertex {
|
||
/// Vertex position.
|
||
pub position: Point3<f64>,
|
||
/// Label/region ID (for multi-material meshes).
|
||
pub label: i64,
|
||
/// Whether this is a boundary vertex.
|
||
pub is_boundary: bool,
|
||
}
|
||
|
||
impl Vertex {
|
||
/// Create a new vertex.
|
||
pub fn new(x: f64, y: f64, z: f64) -> Self {
|
||
Self {
|
||
position: Point3::new(x, y, z),
|
||
label: 0,
|
||
is_boundary: false,
|
||
}
|
||
}
|
||
|
||
/// Create a vertex with a label.
|
||
pub fn with_label(x: f64, y: f64, z: f64, label: i64) -> Self {
|
||
Self {
|
||
position: Point3::new(x, y, z),
|
||
label,
|
||
is_boundary: false,
|
||
}
|
||
}
|
||
}
|
||
|
||
/// A tetrahedron defined by four vertex indices.
|
||
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
|
||
pub struct Tetrahedron {
|
||
/// Vertex indices (0-based).
|
||
pub vertices: [usize; 4],
|
||
/// Material/region label.
|
||
pub label: i64,
|
||
}
|
||
|
||
impl Tetrahedron {
|
||
/// Create a new tetrahedron.
|
||
pub fn new(v0: usize, v1: usize, v2: usize, v3: usize) -> Self {
|
||
Self {
|
||
vertices: [v0, v1, v2, v3],
|
||
label: 0,
|
||
}
|
||
}
|
||
|
||
/// Create a tetrahedron with a label.
|
||
pub fn with_label(v0: usize, v1: usize, v2: usize, v3: usize, label: i64) -> Self {
|
||
Self {
|
||
vertices: [v0, v1, v2, v3],
|
||
label,
|
||
}
|
||
}
|
||
|
||
/// Get the four triangular faces as vertex index tuples.
|
||
pub fn faces(&self) -> [[usize; 3]; 4] {
|
||
let [v0, v1, v2, v3] = self.vertices;
|
||
[[v1, v2, v3], [v0, v3, v2], [v0, v1, v3], [v0, v2, v1]]
|
||
}
|
||
|
||
/// Get the six edges as vertex index pairs.
|
||
pub fn edges(&self) -> [[usize; 2]; 6] {
|
||
let [v0, v1, v2, v3] = self.vertices;
|
||
[[v0, v1], [v0, v2], [v0, v3], [v1, v2], [v1, v3], [v2, v3]]
|
||
}
|
||
}
|
||
|
||
/// A triangular face on the surface.
|
||
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
|
||
pub struct Triangle {
|
||
/// Vertex indices (0-based).
|
||
pub vertices: [usize; 3],
|
||
/// Face normal direction.
|
||
pub normal_outward: bool,
|
||
}
|
||
|
||
impl Triangle {
|
||
/// Create a new triangle.
|
||
pub fn new(v0: usize, v1: usize, v2: usize) -> Self {
|
||
Self {
|
||
vertices: [v0, v1, v2],
|
||
normal_outward: true,
|
||
}
|
||
}
|
||
|
||
/// Get the sorted vertex tuple (for comparison).
|
||
pub fn sorted(&self) -> [usize; 3] {
|
||
let mut v = self.vertices;
|
||
v.sort_unstable();
|
||
v
|
||
}
|
||
}
|
||
|
||
/// A tetrahedral mesh.
|
||
#[derive(Debug, Clone, Serialize, Deserialize)]
|
||
pub struct TetrahedralMesh {
|
||
/// Vertices of the mesh.
|
||
pub vertices: Vec<Vertex>,
|
||
/// Tetrahedra of the mesh.
|
||
pub tetrahedra: Vec<Tetrahedron>,
|
||
/// Surface triangles (computed on demand).
|
||
surface_triangles: Option<Vec<Triangle>>,
|
||
/// Vertex-to-tetrahedron adjacency (computed on demand).
|
||
vertex_tets: Option<Vec<Vec<usize>>>,
|
||
}
|
||
|
||
impl TetrahedralMesh {
|
||
/// Create a new empty mesh.
|
||
pub fn new() -> Self {
|
||
Self {
|
||
vertices: Vec::new(),
|
||
tetrahedra: Vec::new(),
|
||
surface_triangles: None,
|
||
vertex_tets: None,
|
||
}
|
||
}
|
||
|
||
/// Create a mesh with preallocated capacity.
|
||
pub fn with_capacity(num_vertices: usize, num_tetrahedra: usize) -> Self {
|
||
Self {
|
||
vertices: Vec::with_capacity(num_vertices),
|
||
tetrahedra: Vec::with_capacity(num_tetrahedra),
|
||
surface_triangles: None,
|
||
vertex_tets: None,
|
||
}
|
||
}
|
||
|
||
/// Add a vertex and return its index.
|
||
pub fn add_vertex(&mut self, vertex: Vertex) -> usize {
|
||
let idx = self.vertices.len();
|
||
self.vertices.push(vertex);
|
||
self.invalidate_cache();
|
||
idx
|
||
}
|
||
|
||
/// Add a tetrahedron and return its index.
|
||
pub fn add_tetrahedron(&mut self, tet: Tetrahedron) -> usize {
|
||
let idx = self.tetrahedra.len();
|
||
self.tetrahedra.push(tet);
|
||
self.invalidate_cache();
|
||
idx
|
||
}
|
||
|
||
/// Get the number of vertices.
|
||
pub fn num_vertices(&self) -> usize {
|
||
self.vertices.len()
|
||
}
|
||
|
||
/// Get the number of tetrahedra.
|
||
pub fn num_tetrahedra(&self) -> usize {
|
||
self.tetrahedra.len()
|
||
}
|
||
|
||
/// Calculate the volume of a tetrahedron.
|
||
pub fn tetrahedron_volume(&self, tet_idx: usize) -> f64 {
|
||
let tet = &self.tetrahedra[tet_idx];
|
||
let p0 = &self.vertices[tet.vertices[0]].position;
|
||
let p1 = &self.vertices[tet.vertices[1]].position;
|
||
let p2 = &self.vertices[tet.vertices[2]].position;
|
||
let p3 = &self.vertices[tet.vertices[3]].position;
|
||
|
||
let v01 = p1 - p0;
|
||
let v02 = p2 - p0;
|
||
let v03 = p3 - p0;
|
||
|
||
(v01.cross(&v02).dot(&v03) / 6.0).abs()
|
||
}
|
||
|
||
/// Calculate the total volume of the mesh.
|
||
pub fn total_volume(&self) -> f64 {
|
||
(0..self.tetrahedra.len())
|
||
.map(|i| self.tetrahedron_volume(i))
|
||
.sum()
|
||
}
|
||
|
||
/// Get the centroid of a tetrahedron.
|
||
pub fn tetrahedron_centroid(&self, tet_idx: usize) -> Point3<f64> {
|
||
let tet = &self.tetrahedra[tet_idx];
|
||
let mut centroid = Point3::origin();
|
||
for &vi in &tet.vertices {
|
||
centroid += self.vertices[vi].position.coords;
|
||
}
|
||
centroid / 4.0
|
||
}
|
||
|
||
/// Calculate the circumradius of a tetrahedron.
|
||
pub fn tetrahedron_circumradius(&self, tet_idx: usize) -> f64 {
|
||
let tet = &self.tetrahedra[tet_idx];
|
||
let p0 = &self.vertices[tet.vertices[0]].position;
|
||
let p1 = &self.vertices[tet.vertices[1]].position;
|
||
let p2 = &self.vertices[tet.vertices[2]].position;
|
||
let p3 = &self.vertices[tet.vertices[3]].position;
|
||
|
||
// Use the formula: R = |e1 × e2 · e3| / (6V)
|
||
// where e1, e2, e3 are edge vectors from p0
|
||
let e1 = p1 - p0;
|
||
let e2 = p2 - p0;
|
||
let e3 = p3 - p0;
|
||
|
||
let vol = (e1.cross(&e2).dot(&e3) / 6.0).abs();
|
||
if vol < 1e-15 {
|
||
return f64::INFINITY;
|
||
}
|
||
|
||
// Calculate edge lengths
|
||
let a = (p1 - p2).norm();
|
||
let b = (p0 - p2).norm();
|
||
let c = (p0 - p1).norm();
|
||
let d = (p0 - p3).norm();
|
||
let e = (p1 - p3).norm();
|
||
let f = (p2 - p3).norm();
|
||
|
||
// Circumradius formula
|
||
let p = (a * d) * (b * e + c * f);
|
||
let q = (b * e) * (a * d + c * f);
|
||
let r = (c * f) * (a * d + b * e);
|
||
|
||
let numerator = ((p + q + r) * (-p + q + r) * (p - q + r) * (p + q - r)).sqrt();
|
||
numerator / (24.0 * vol)
|
||
}
|
||
|
||
/// Get the shortest edge length of a tetrahedron.
|
||
pub fn tetrahedron_min_edge(&self, tet_idx: usize) -> f64 {
|
||
let tet = &self.tetrahedra[tet_idx];
|
||
let mut min_edge = f64::INFINITY;
|
||
|
||
for edge in tet.edges() {
|
||
let p0 = &self.vertices[edge[0]].position;
|
||
let p1 = &self.vertices[edge[1]].position;
|
||
let len = (p1 - p0).norm();
|
||
min_edge = min_edge.min(len);
|
||
}
|
||
|
||
min_edge
|
||
}
|
||
|
||
/// Calculate the radius-edge ratio of a tetrahedron.
|
||
pub fn tetrahedron_radius_edge_ratio(&self, tet_idx: usize) -> f64 {
|
||
let circumradius = self.tetrahedron_circumradius(tet_idx);
|
||
let min_edge = self.tetrahedron_min_edge(tet_idx);
|
||
|
||
if min_edge < 1e-15 {
|
||
return f64::INFINITY;
|
||
}
|
||
|
||
circumradius / min_edge
|
||
}
|
||
|
||
/// Get the bounding box of the mesh.
|
||
pub fn bounding_box(&self) -> Option<(Point3<f64>, Point3<f64>)> {
|
||
if self.vertices.is_empty() {
|
||
return None;
|
||
}
|
||
|
||
let mut min = self.vertices[0].position;
|
||
let mut max = self.vertices[0].position;
|
||
|
||
for v in &self.vertices {
|
||
min.x = min.x.min(v.position.x);
|
||
min.y = min.y.min(v.position.y);
|
||
min.z = min.z.min(v.position.z);
|
||
max.x = max.x.max(v.position.x);
|
||
max.y = max.y.max(v.position.y);
|
||
max.z = max.z.max(v.position.z);
|
||
}
|
||
|
||
Some((min, max))
|
||
}
|
||
|
||
/// Extract surface triangles.
|
||
pub fn surface_triangles(&mut self) -> &[Triangle] {
|
||
if self.surface_triangles.is_none() {
|
||
self.compute_surface();
|
||
}
|
||
self.surface_triangles.as_ref().unwrap()
|
||
}
|
||
|
||
/// Compute surface triangles by finding faces shared by only one tetrahedron.
|
||
fn compute_surface(&mut self) {
|
||
let mut face_count: HashMap<[usize; 3], (usize, Triangle)> = HashMap::new();
|
||
|
||
for tet in &self.tetrahedra {
|
||
for face in tet.faces() {
|
||
let mut sorted = face;
|
||
sorted.sort_unstable();
|
||
|
||
face_count
|
||
.entry(sorted)
|
||
.and_modify(|(count, _)| *count += 1)
|
||
.or_insert((1, Triangle::new(face[0], face[1], face[2])));
|
||
}
|
||
}
|
||
|
||
// Surface faces are those that appear exactly once
|
||
let surface: Vec<Triangle> = face_count
|
||
.into_iter()
|
||
.filter(|(_, (count, _))| *count == 1)
|
||
.map(|(_, (_, tri))| tri)
|
||
.collect();
|
||
|
||
self.surface_triangles = Some(surface);
|
||
}
|
||
|
||
/// Build vertex-to-tetrahedron adjacency.
|
||
pub fn vertex_tetrahedra(&mut self) -> &[Vec<usize>] {
|
||
if self.vertex_tets.is_none() {
|
||
self.compute_vertex_tets();
|
||
}
|
||
self.vertex_tets.as_ref().unwrap()
|
||
}
|
||
|
||
fn compute_vertex_tets(&mut self) {
|
||
let mut v2t: Vec<Vec<usize>> = vec![Vec::new(); self.vertices.len()];
|
||
|
||
for (ti, tet) in self.tetrahedra.iter().enumerate() {
|
||
for &vi in &tet.vertices {
|
||
v2t[vi].push(ti);
|
||
}
|
||
}
|
||
|
||
self.vertex_tets = Some(v2t);
|
||
}
|
||
|
||
/// Get unique labels in the mesh.
|
||
pub fn labels(&self) -> Vec<i64> {
|
||
let mut labels: HashSet<i64> = HashSet::new();
|
||
for tet in &self.tetrahedra {
|
||
labels.insert(tet.label);
|
||
}
|
||
let mut result: Vec<i64> = labels.into_iter().collect();
|
||
result.sort_unstable();
|
||
result
|
||
}
|
||
|
||
/// Get tetrahedra with a specific label.
|
||
pub fn tetrahedra_by_label(&self, label: i64) -> Vec<usize> {
|
||
self.tetrahedra
|
||
.iter()
|
||
.enumerate()
|
||
.filter(|(_, t)| t.label == label)
|
||
.map(|(i, _)| i)
|
||
.collect()
|
||
}
|
||
|
||
/// Invalidate cached computations.
|
||
fn invalidate_cache(&mut self) {
|
||
self.surface_triangles = None;
|
||
self.vertex_tets = None;
|
||
}
|
||
}
|
||
|
||
impl Default for TetrahedralMesh {
|
||
fn default() -> Self {
|
||
Self::new()
|
||
}
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
|
||
fn create_single_tet() -> TetrahedralMesh {
|
||
let mut mesh = TetrahedralMesh::new();
|
||
|
||
mesh.add_vertex(Vertex::new(0.0, 0.0, 0.0));
|
||
mesh.add_vertex(Vertex::new(1.0, 0.0, 0.0));
|
||
mesh.add_vertex(Vertex::new(0.5, 1.0, 0.0));
|
||
mesh.add_vertex(Vertex::new(0.5, 0.5, 1.0));
|
||
|
||
mesh.add_tetrahedron(Tetrahedron::new(0, 1, 2, 3));
|
||
|
||
mesh
|
||
}
|
||
|
||
#[test]
|
||
fn test_mesh_creation() {
|
||
let mesh = create_single_tet();
|
||
assert_eq!(mesh.num_vertices(), 4);
|
||
assert_eq!(mesh.num_tetrahedra(), 1);
|
||
}
|
||
|
||
#[test]
|
||
fn test_tetrahedron_volume() {
|
||
let mesh = create_single_tet();
|
||
let vol = mesh.tetrahedron_volume(0);
|
||
// Volume of this tet should be approximately 1/6
|
||
assert!((vol - 1.0 / 6.0).abs() < 0.01);
|
||
}
|
||
|
||
#[test]
|
||
fn test_surface_triangles() {
|
||
let mut mesh = create_single_tet();
|
||
let surface = mesh.surface_triangles();
|
||
// A single tetrahedron has 4 surface triangles
|
||
assert_eq!(surface.len(), 4);
|
||
}
|
||
|
||
#[test]
|
||
fn test_bounding_box() {
|
||
let mesh = create_single_tet();
|
||
let (min, max) = mesh.bounding_box().unwrap();
|
||
assert_eq!(min.x, 0.0);
|
||
assert_eq!(max.x, 1.0);
|
||
}
|
||
}
|