#[cfg(feature = "parallel")]
use rayon::prelude::*;
use approx::AbsDiffEq;
use na::{self, RealField, Unit};
use crate::geometry::ParticlesContacts;
use crate::math::{Real, Vector};
use crate::object::{Boundary, Fluid};
use crate::solver::NonPressureForce;
use crate::TimestepManager;
#[derive(Clone)]
pub struct Akinci2013SurfaceTension {
fluid_tension_coefficient: Real,
boundary_adhesion_coefficient: Real,
normals: Vec<Vector<Real>>,
}
impl Akinci2013SurfaceTension {
pub fn new(fluid_tension_coefficient: Real, boundary_adhesion_coefficient: Real) -> Self {
Self {
fluid_tension_coefficient,
boundary_adhesion_coefficient,
normals: Vec::new(),
}
}
fn init(&mut self, fluid: &Fluid) {
if self.normals.len() != fluid.num_particles() {
self.normals.resize(fluid.num_particles(), Vector::zeros());
}
}
fn compute_normals(
&mut self,
kernel_radius: Real,
fluid_fluid_contacts: &ParticlesContacts,
fluid: &Fluid,
densities: &[Real],
) {
par_iter_mut!(self.normals)
.enumerate()
.for_each(|(i, normal_i)| {
let mut normal = Vector::zeros();
for c in fluid_fluid_contacts
.particle_contacts(i)
.read()
.unwrap()
.iter()
{
if c.i_model == c.j_model {
normal += c.gradient * (fluid.particle_mass(c.j) / densities[c.j]);
}
}
*normal_i = normal * kernel_radius;
})
}
}
fn cohesion_kernel(r: Real, h: Real) -> Real {
#[cfg(feature = "dim3")]
let normalizer = na::convert::<_, Real>(32.0f64) / (Real::pi() * h.powi(9));
#[cfg(feature = "dim2")]
let normalizer = na::convert::<_, Real>(32.0f64) / (Real::pi() * h.powi(8));
let _2: Real = na::convert::<_, Real>(2.0f64);
let coeff = if r <= h / _2 {
_2 * (h - r).powi(3) * r.powi(3) - h.powi(6) / na::convert::<_, Real>(64.0f64)
} else if r <= h {
(h - r).powi(3) * r.powi(3)
} else {
na::zero::<Real>()
};
normalizer * coeff
}
fn adhesion_kernel(r: Real, h: Real) -> Real {
let _2: Real = na::convert::<_, Real>(2.0f64);
if r > h / _2 && r <= h {
let _4: Real = na::convert::<_, Real>(4.0f64);
let _6: Real = na::convert::<_, Real>(6.0f64);
#[cfg(feature = "dim3")]
let normalizer = na::convert::<_, Real>(0.007) / h.powf(na::convert::<_, Real>(3.25f64));
#[cfg(feature = "dim2")]
let normalizer = na::convert::<_, Real>(0.007) / h.powf(na::convert::<_, Real>(2.25f64));
let coeff = ((-_4 * r * r / h + _6 * r - _2 * h).max(na::zero::<Real>()))
.powf(na::convert::<_, Real>(0.25));
normalizer * coeff
} else {
na::zero::<Real>()
}
}
impl NonPressureForce for Akinci2013SurfaceTension {
fn solve(
&mut self,
_timestep: &TimestepManager,
kernel_radius: Real,
fluid_fluid_contacts: &ParticlesContacts,
fluid_boundaries_contacts: &ParticlesContacts,
fluid: &mut Fluid,
boundaries: &[Boundary],
densities: &[Real],
) {
self.init(fluid);
let _2: Real = na::convert::<_, Real>(2.0f64);
self.compute_normals(kernel_radius, fluid_fluid_contacts, fluid, densities);
let normals = &self.normals;
let fluid_tension_coefficient = self.fluid_tension_coefficient;
let boundary_adhesion_coefficient = self.boundary_adhesion_coefficient;
let volumes = &mut fluid.volumes;
let density0 = fluid.density0;
let positions = &fluid.positions;
par_iter_mut!(fluid.accelerations)
.enumerate()
.for_each(|(i, acceleration_i)| {
if self.fluid_tension_coefficient != na::zero::<Real>() {
for c in fluid_fluid_contacts
.particle_contacts(i)
.read()
.unwrap()
.iter()
{
if c.i_model == c.j_model {
let dpos = positions[c.i] - positions[c.j];
let cohesion_vec = if let Some((dir, dist)) =
Unit::try_new_and_get(dpos, Real::default_epsilon())
{
*dir * cohesion_kernel(dist, kernel_radius)
} else {
Vector::zeros()
};
let cohesion_acc = cohesion_vec
* (-fluid_tension_coefficient * volumes[c.j] * density0);
let curvature_acc =
(normals[c.i] - normals[c.j]) * -fluid_tension_coefficient;
let kij = _2 * density0 / (densities[c.i] + densities[c.j]);
*acceleration_i += (curvature_acc + cohesion_acc) * kij;
}
}
}
if boundary_adhesion_coefficient != na::zero::<Real>() {
for c in fluid_boundaries_contacts
.particle_contacts(i)
.read()
.unwrap()
.iter()
{
let dpos = positions[c.i] - boundaries[c.j_model].positions[c.j];
let adhesion_vec = if let Some((dir, dist)) =
Unit::try_new_and_get(dpos, Real::default_epsilon())
{
*dir * adhesion_kernel(dist, kernel_radius)
} else {
Vector::zeros()
};
let mi = volumes[c.i] * density0;
let mj = boundaries[c.j_model].volumes[c.j] * density0;
let adhesion_acc = adhesion_vec * (boundary_adhesion_coefficient * mj);
*acceleration_i -= adhesion_acc;
boundaries[c.j_model].apply_force(c.j, adhesion_acc * mi);
}
}
})
}
fn apply_permutation(&mut self, _: &[usize]) {}
}