mod extended;
mod surface;
pub use extended::*;
pub use surface::*;
use crate::mesh::attrib::*;
use crate::mesh::topology::*;
use crate::mesh::vertex_positions::*;
use crate::prim::Tetrahedron;
use crate::Real;
use math::{Matrix3, SquareMatrix};
use std::slice::{Iter, IterMut};
#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};
#[derive(Clone, Debug, PartialEq, Attrib, Intrinsic)]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
pub struct TetMesh<T: Real> {
#[intrinsic(VertexPositions)]
pub vertex_positions: IntrinsicAttribute<[T; 3], VertexIndex>,
pub indices: IntrinsicAttribute<[usize; 4], CellIndex>,
pub vertex_attributes: AttribDict<VertexIndex>,
pub cell_attributes: AttribDict<CellIndex>,
pub cell_vertex_attributes: AttribDict<CellVertexIndex>,
pub cell_face_attributes: AttribDict<CellFaceIndex>,
}
impl<T: Real> TetMesh<T> {
pub const TET_FACES: [[usize; 3]; 4] = [[1, 3, 2], [0, 2, 3], [0, 3, 1], [0, 1, 2]];
pub fn new(verts: Vec<[T; 3]>, indices: Vec<usize>) -> TetMesh<T> {
assert_eq!(indices.len() % 4, 0);
TetMesh {
vertex_positions: IntrinsicAttribute::from_vec(verts),
indices: IntrinsicAttribute::from_vec(crate::into_chunked_vec4(indices)),
vertex_attributes: AttribDict::new(),
cell_attributes: AttribDict::new(),
cell_vertex_attributes: AttribDict::new(),
cell_face_attributes: AttribDict::new(),
}
}
#[inline]
pub fn cell_iter(&self) -> Iter<[usize; 4]> {
self.indices.iter()
}
#[cfg(feature = "rayon")]
#[inline]
pub fn cell_par_iter(&self) -> rayon::slice::Iter<[usize; 4]> {
self.indices.par_iter()
}
#[inline]
pub fn cell_iter_mut(&mut self) -> IterMut<[usize; 4]> {
self.indices.iter_mut()
}
#[cfg(feature = "rayon")]
#[inline]
pub fn cell_par_iter_mut(&mut self) -> rayon::slice::IterMut<[usize; 4]> {
self.indices.par_iter_mut()
}
#[inline]
pub fn cell<CI: Into<CellIndex>>(&self, cidx: CI) -> &[usize; 4] {
&self.indices[cidx.into()]
}
#[inline]
pub fn cells(&self) -> &[[usize; 4]] {
self.indices.as_slice()
}
#[inline]
pub fn tet_iter<'a>(&'a self) -> impl Iterator<Item = Tetrahedron<T>> + 'a {
self.cell_iter().map(move |tet| self.tet_from_indices(tet))
}
#[inline]
pub fn tet_from_indices(&self, indices: &[usize; 4]) -> Tetrahedron<T> {
Tetrahedron::from_indexed_slice(indices, self.vertex_positions.as_slice())
}
#[inline]
pub fn tet<CI: Into<CellIndex>>(&self, cidx: CI) -> Tetrahedron<T> {
self.tet_from_indices(self.cell(cidx))
}
#[inline]
pub fn inverted(mut self) -> TetMesh<T> {
self.invert();
self
}
#[inline]
fn invert_tet_cell(cell: &mut [usize; 4]) {
cell.swap(2, 3);
}
#[inline]
pub fn invert(&mut self) {
for cell in self.indices.iter_mut() {
Self::invert_tet_cell(cell);
}
}
#[inline]
pub fn canonicalized(mut self) -> TetMesh<T> {
self.canonicalize();
self
}
#[inline]
pub fn canonicalize(&mut self) {
use crate::ops::ShapeMatrix;
let TetMesh {
ref vertex_positions,
ref mut indices,
..
} = *self;
for cell in indices.iter_mut() {
let tet = Tetrahedron::from_indexed_slice(cell, vertex_positions.as_slice());
if Matrix3::from(tet.shape_matrix()).determinant() < T::zero() {
Self::invert_tet_cell(cell);
}
}
}
}
impl<T: Real> Default for TetMesh<T> {
fn default() -> Self {
TetMesh::new(vec![], vec![])
}
}
impl<T: Real> NumVertices for TetMesh<T> {
#[inline]
fn num_vertices(&self) -> usize {
self.vertex_positions.len()
}
}
impl<T: Real> NumCells for TetMesh<T> {
#[inline]
fn num_cells(&self) -> usize {
self.indices.len()
}
}
impl<T: Real> CellVertex for TetMesh<T> {
#[inline]
fn vertex<CVI>(&self, cv_idx: CVI) -> VertexIndex
where
CVI: Copy + Into<CellVertexIndex>,
{
let cv_idx = usize::from(cv_idx.into());
debug_assert!(cv_idx < self.num_cell_vertices());
self.indices[cv_idx / 4][cv_idx % 4].into()
}
#[inline]
fn cell_vertex<CI>(&self, cidx: CI, which: usize) -> Option<CellVertexIndex>
where
CI: Copy + Into<CellIndex>,
{
if which >= 4 {
None
} else {
Some((4 * usize::from(cidx.into()) + which).into())
}
}
#[inline]
fn num_cell_vertices(&self) -> usize {
self.indices.len() * 4
}
#[inline]
fn num_vertices_at_cell<CI>(&self, _: CI) -> usize
where
CI: Copy + Into<CellIndex>,
{
4
}
}
impl<T: Real> CellFace for TetMesh<T> {
#[inline]
fn face<CFI>(&self, cf_idx: CFI) -> FaceIndex
where
CFI: Copy + Into<CellFaceIndex>,
{
let cf_idx = usize::from(cf_idx.into());
debug_assert!(cf_idx < self.num_cell_faces());
self.indices[cf_idx / 4][cf_idx % 4].into()
}
#[inline]
fn cell_face<CI>(&self, cidx: CI, which: usize) -> Option<CellFaceIndex>
where
CI: Copy + Into<CellIndex>,
{
if which >= 4 {
None
} else {
Some((4 * usize::from(cidx.into()) + which).into())
}
}
#[inline]
fn num_cell_faces(&self) -> usize {
self.indices.len() * 4
}
#[inline]
fn num_faces_at_cell<CI>(&self, _: CI) -> usize
where
CI: Copy + Into<CellIndex>,
{
4
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::index::Index;
fn simple_tetmesh() -> TetMesh<f64> {
let points = vec![
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[1.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
[1.0, 0.0, 1.0],
];
let indices = vec![0, 4, 2, 5, 0, 5, 2, 3, 5, 3, 0, 1];
TetMesh::new(points, indices)
}
#[test]
fn tetmesh_test() {
let mesh = simple_tetmesh();
assert_eq!(mesh.num_vertices(), 6);
assert_eq!(mesh.num_cells(), 3);
assert_eq!(mesh.num_cell_vertices(), 12);
assert_eq!(mesh.num_cell_faces(), 12);
assert_eq!(Index::from(mesh.cell_vertex(1, 1)), 5);
assert_eq!(Index::from(mesh.cell_vertex(0, 2)), 2);
assert_eq!(Index::from(mesh.cell_face(2, 3)), 11);
}
#[test]
fn tet_iter_test() {
use math::Vector3;
let mesh = simple_tetmesh();
let points = mesh.vertex_positions().to_vec();
let pt = |i| Vector3::from(points[i]);
let tets = vec![
Tetrahedron(pt(0), pt(4), pt(2), pt(5)),
Tetrahedron(pt(0), pt(5), pt(2), pt(3)),
Tetrahedron(pt(5), pt(3), pt(0), pt(1)),
];
for (tet, exptet) in mesh.tet_iter().zip(tets.into_iter()) {
assert_eq!(tet, exptet);
}
}
#[test]
fn invert_test() {
use crate::ops::Volume;
let mut mesh = simple_tetmesh();
let mut vols = Vec::new();
for ref tet in mesh.tet_iter() {
vols.push(tet.volume());
assert!(tet.signed_volume() > 0.0);
}
mesh.invert();
for (tet, vol) in mesh.tet_iter().zip(vols) {
assert_eq!(tet.signed_volume(), -vol);
}
}
#[test]
fn canonicalize_test() {
use crate::ops::Volume;
let points = vec![
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[1.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
[1.0, 0.0, 1.0],
];
let indices = vec![0, 4, 2, 5, 0, 5, 3, 2, 5, 3, 1, 0];
let mut mesh = TetMesh::new(points.clone(), indices);
let vols: Vec<_> = mesh.tet_iter().map(|t| t.volume()).collect();
assert!(mesh.tet(0).signed_volume() > 0.0);
assert!(mesh.tet(1).signed_volume() < 0.0);
assert!(mesh.tet(2).signed_volume() < 0.0);
mesh.canonicalize();
for (tet, vol) in mesh.tet_iter().zip(vols) {
assert_eq!(tet.signed_volume(), vol);
}
}
}