use super::SoftBody;
use crate::math::{DIM, Matrix, Real, Vector};
#[cfg(not(feature = "std"))]
#[allow(unused_imports)]
use simba::scalar::{ComplexField as _, RealField as _};
impl SoftBody {
pub(crate) fn update_orientation(&mut self) {
let (particles, boundary) = (&self.particles, &self.boundary);
for piece in &mut self.volume_pieces {
if piece.rest_volume != 0.0 {
let volume = Self::elements_volume(piece.boundary_elements(boundary), |i| {
particles[i as usize].position
});
let ratio = volume / piece.rest_volume;
if ratio < -0.25 {
piece.inverted = true;
} else if ratio > 0.25 {
piece.inverted = false;
}
}
}
self.for_each_mesh_mut(|body, mesh| mesh.update_orientation(body));
}
pub(crate) fn boundary_volume(
boundary: &[[u32; DIM]],
position: impl Fn(u32) -> Vector,
) -> Real {
Self::elements_volume(boundary, position)
}
pub(crate) fn elements_volume<'a>(
elements: impl IntoIterator<Item = &'a [u32; DIM]>,
position: impl Fn(u32) -> Vector,
) -> Real {
let mut vol = 0.0;
for element in elements {
#[cfg(feature = "dim2")]
{
let a = position(element[0]);
let b = position(element[1]);
vol += a.perp_dot(b) * 0.5;
}
#[cfg(feature = "dim3")]
{
let a = position(element[0]);
let b = position(element[1]);
let c = position(element[2]);
vol += a.dot(b.cross(c)) / 6.0;
}
}
vol
}
pub(crate) fn boundary_volume_gradients<'a>(
boundary: impl IntoIterator<Item = &'a [u32; DIM]>,
position: impl Fn(u32) -> Vector,
gradients: &mut [Vector],
) {
for element in boundary {
#[cfg(feature = "dim2")]
{
let a = position(element[0]);
let b = position(element[1]);
gradients[element[0] as usize] += Vector::new(b.y, -b.x) * 0.5;
gradients[element[1] as usize] += Vector::new(-a.y, a.x) * 0.5;
}
#[cfg(feature = "dim3")]
{
let a = position(element[0]);
let b = position(element[1]);
let c = position(element[2]);
gradients[element[0] as usize] += b.cross(c) / 6.0;
gradients[element[1] as usize] += c.cross(a) / 6.0;
gradients[element[2] as usize] += a.cross(b) / 6.0;
}
}
}
pub(crate) fn cell_volume(x: [Vector; DIM + 1]) -> Real {
#[cfg(feature = "dim2")]
{
(x[1] - x[0]).perp_dot(x[2] - x[0]) * 0.5
}
#[cfg(feature = "dim3")]
{
(x[1] - x[0]).dot((x[2] - x[0]).cross(x[3] - x[0])) / 6.0
}
}
pub(crate) fn cell_volume_gradients(x: [Vector; DIM + 1]) -> [Vector; DIM + 1] {
let mut grad = [Vector::ZERO; DIM + 1];
#[cfg(feature = "dim2")]
{
let a = x[1] - x[0];
let b = x[2] - x[0];
grad[1] = Vector::new(b.y, -b.x) * 0.5;
grad[2] = Vector::new(-a.y, a.x) * 0.5;
grad[0] = -(grad[1] + grad[2]);
}
#[cfg(feature = "dim3")]
{
let a = x[1] - x[0];
let b = x[2] - x[0];
let c = x[3] - x[0];
grad[1] = b.cross(c) / 6.0;
grad[2] = c.cross(a) / 6.0;
grad[3] = a.cross(b) / 6.0;
grad[0] = -(grad[1] + grad[2] + grad[3]);
}
grad
}
pub(crate) fn cell_edge_matrix(x: [Vector; DIM + 1]) -> Matrix {
#[cfg(feature = "dim2")]
{
Matrix::from_cols(x[1] - x[0], x[2] - x[0])
}
#[cfg(feature = "dim3")]
{
Matrix::from_cols(x[1] - x[0], x[2] - x[0], x[3] - x[0])
}
}
}