use math::{BaseNum, BaseFloat, Vector3, Matrix3, SquareMatrix};
use crate::ops::{Centroid, ShapeMatrix, Volume};
use std::ops::{Add, Sub, Mul};
use crate::Pod;
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct Tetrahedron<T>(
pub Vector3<T>,
pub Vector3<T>,
pub Vector3<T>,
pub Vector3<T>,
);
impl<T: Clone> Tetrahedron<T> {
#[inline]
pub fn from_indexed_slice<V: Into<[T; 3]> + Clone>(
indices: &[usize; 4],
slice: &[V],
) -> Tetrahedron<T> {
Tetrahedron(
slice[indices[0]].clone().into().into(),
slice[indices[1]].clone().into().into(),
slice[indices[2]].clone().into().into(),
slice[indices[3]].clone().into().into(),
)
}
#[inline]
pub fn new([a,b,c,d]: [[T; 3]; 4]) -> Self {
Tetrahedron(
a.into(),
b.into(),
c.into(),
d.into(),
)
}
}
impl<T: Pod> Tetrahedron<T> {
#[inline]
pub fn as_array(&self) -> &[[T; 3]; 4] {
debug_assert_eq!(
std::mem::size_of::<[[T; 3]; 4]>(),
std::mem::size_of::<Tetrahedron<T>>(),
);
unsafe { std::mem::transmute_copy(&self) }
}
}
impl<T> Tetrahedron<T> {
#[inline]
pub fn into_array(self) -> [[T; 3]; 4] {
[self.0.into(), self.1.into(), self.2.into(), self.3.into()]
}
}
impl<T: BaseNum> Add for Tetrahedron<T> {
type Output = Tetrahedron<T>;
fn add(self, other: Tetrahedron<T>) -> Tetrahedron<T> {
Tetrahedron(
self.0 + other.0,
self.1 + other.1,
self.2 + other.2,
self.3 + other.3,
)
}
}
impl<'a, T: BaseNum> Add for &'a Tetrahedron<T> {
type Output = Tetrahedron<T>;
fn add(self, other: &Tetrahedron<T>) -> Tetrahedron<T> {
Tetrahedron(
self.0 + other.0,
self.1 + other.1,
self.2 + other.2,
self.3 + other.3,
)
}
}
impl<T: BaseNum> Sub for Tetrahedron<T> {
type Output = Tetrahedron<T>;
fn sub(self, other: Tetrahedron<T>) -> Tetrahedron<T> {
Tetrahedron(
self.0 - other.0,
self.1 - other.1,
self.2 - other.2,
self.3 - other.3,
)
}
}
impl<'a, T: BaseNum> Sub for &'a Tetrahedron<T> {
type Output = Tetrahedron<T>;
fn sub(self, other: &Tetrahedron<T>) -> Tetrahedron<T> {
Tetrahedron(
self.0 - other.0,
self.1 - other.1,
self.2 - other.2,
self.3 - other.3,
)
}
}
impl<T: BaseNum> Mul<T> for Tetrahedron<T> {
type Output = Self;
fn mul(self, rhs: T) -> Self {
Tetrahedron(self.0 * rhs, self.1 * rhs, self.2 * rhs, self.3 * rhs)
}
}
impl<'a, T: BaseNum> std::ops::Mul<T> for &'a Tetrahedron<T> {
type Output = Tetrahedron<T>;
fn mul(self, rhs: T) -> Tetrahedron<T> {
Tetrahedron(self.0 * rhs, self.1 * rhs, self.2 * rhs, self.3 * rhs)
}
}
impl<'a, T: BaseNum + num_traits::FromPrimitive> Centroid<[T; 3]> for &'a Tetrahedron<T> {
#[inline]
fn centroid(self) -> [T; 3] {
((self.0 + self.1 + self.2 + self.3) * T::from(0.25).unwrap()).into()
}
}
impl<'a, T: BaseNum> ShapeMatrix<[[T; 3]; 3]> for &'a Tetrahedron<T> {
#[inline]
fn shape_matrix(self) -> [[T; 3]; 3] {
[
(self.0 - self.3).into(),
(self.1 - self.3).into(),
(self.2 - self.3).into(),
]
}
}
impl<T: BaseNum> ShapeMatrix<[[T; 3]; 3]> for Tetrahedron<T> {
#[inline]
fn shape_matrix(self) -> [[T; 3]; 3] {
[
(self.0 - self.3).into(),
(self.1 - self.3).into(),
(self.2 - self.3).into(),
]
}
}
impl<'a, T: BaseFloat> Volume<T> for &'a Tetrahedron<T> {
#[inline]
fn volume(self) -> T {
self.signed_volume().abs()
}
#[inline]
fn signed_volume(self) -> T {
Matrix3::from(self.shape_matrix()).determinant() / T::from(6.0).unwrap()
}
}
impl<T: BaseFloat> Volume<T> for Tetrahedron<T> {
#[inline]
fn volume(self) -> T {
self.signed_volume().abs()
}
#[inline]
fn signed_volume(self) -> T {
Matrix3::from(self.shape_matrix()).determinant() / T::from(6.0).unwrap()
}
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
fn tet() -> Tetrahedron<f64> {
Tetrahedron(
Vector3::from([0.0, 1.0, 0.0]),
Vector3::from([-0.94281, -0.33333, 0.0]),
Vector3::from([0.471405, -0.33333, 0.816498]),
Vector3::from([0.471405, -0.33333, -0.816498]),
)
}
#[test]
fn volume_test() {
assert_relative_eq!(tet().volume(), 0.5132, epsilon = 1e-4);
assert_relative_eq!(tet().signed_volume(), 0.5132, epsilon = 1e-4);
}
#[test]
fn shape_matrix_test() {
let mtx = tet().shape_matrix();
assert_relative_eq!(mtx[0][0], -0.471405, epsilon = 1e-4);
assert_relative_eq!(mtx[0][1], 1.33333, epsilon = 1e-4);
assert_relative_eq!(mtx[0][2], 0.816498, epsilon = 1e-4);
assert_relative_eq!(mtx[1][0], -1.41421, epsilon = 1e-4);
assert_relative_eq!(mtx[1][1], 0.0, epsilon = 1e-4);
assert_relative_eq!(mtx[1][2], 0.816498, epsilon = 1e-4);
assert_relative_eq!(mtx[2][0], 0.0, epsilon = 1e-4);
assert_relative_eq!(mtx[2][1], 0.0, epsilon = 1e-4);
assert_relative_eq!(mtx[2][2], 1.633, epsilon = 1e-4);
}
#[test]
fn centroid_test() {
let c = tet().centroid();
assert_relative_eq!(c[0], 0.0, epsilon = 1e-4);
assert_relative_eq!(c[1], 0.0, epsilon = 1e-4);
assert_relative_eq!(c[2], 0.0, epsilon = 1e-4);
}
}