runmat-meshing-core 0.6.0

Deterministic meshing preparation contracts for RunMat
Documentation
use super::{
    types::{Point3, Tetrahedron3, Triangle3},
    vector::{add, cross, distance, dot, norm, scale, sub},
};

pub fn triangle_centroid(points: Triangle3) -> Point3 {
    scale(add(add(points[0], points[1]), points[2]), 1.0 / 3.0)
}

pub fn tetrahedron_centroid(points: Tetrahedron3) -> Point3 {
    scale(
        add(add(points[0], points[1]), add(points[2], points[3])),
        0.25,
    )
}

pub fn triangle_area(points: Triangle3) -> f64 {
    0.5 * norm(cross(sub(points[1], points[0]), sub(points[2], points[0])))
}

pub fn tetrahedron_signed_volume(points: Tetrahedron3) -> f64 {
    dot(
        sub(points[1], points[0]),
        cross(sub(points[2], points[0]), sub(points[3], points[0])),
    ) / 6.0
}

pub fn tetrahedron_volume(points: Tetrahedron3) -> f64 {
    tetrahedron_signed_volume(points).abs()
}

pub fn orient_tetrahedron_node_ids(
    mut node_ids: [u32; 4],
    points: Tetrahedron3,
) -> ([u32; 4], f64) {
    let mut signed_volume = tetrahedron_signed_volume(points);
    if signed_volume < 0.0 {
        node_ids.swap(1, 2);
        signed_volume = -signed_volume;
    }
    (node_ids, signed_volume)
}

pub fn tetrahedron_edge_aspect_ratio(points: Tetrahedron3) -> f64 {
    let mut min_edge = f64::INFINITY;
    let mut max_edge = 0.0_f64;
    for left_index in 0..4 {
        for right_index in (left_index + 1)..4 {
            let length = distance(points[left_index], points[right_index]);
            min_edge = min_edge.min(length);
            max_edge = max_edge.max(length);
        }
    }
    max_edge / min_edge.max(f64::EPSILON)
}

pub fn tetrahedron_scaled_jacobian(points: Tetrahedron3) -> f64 {
    let corners = [
        (points[0], points[1], points[2], points[3]),
        (points[1], points[0], points[3], points[2]),
        (points[2], points[0], points[1], points[3]),
        (points[3], points[0], points[2], points[1]),
    ];
    corners
        .into_iter()
        .map(|(origin, first, second, third)| {
            let first = sub(first, origin);
            let second = sub(second, origin);
            let third = sub(third, origin);
            let denominator = norm(first) * norm(second) * norm(third);
            if denominator <= f64::EPSILON {
                return 0.0;
            }
            (2.0_f64.sqrt() * dot(first, cross(second, third)) / denominator).abs()
        })
        .fold(f64::INFINITY, f64::min)
}