1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
use crate::coupling::CouplingManager;
use crate::geometry::{HGrid, HGridEntry};
use crate::object::{BoundaryHandle, BoundarySet, Fluid};
use crate::solver::DFSPHSolver;
use crate::LiquidWorld;
use crate::TimestepManager;
use approx::AbsDiffEq;
use na::Unit;
use rapier::dynamics::RigidBodySet;
use rapier::geometry::{ColliderHandle, ColliderSet};
use rapier::math::{Point, Vector};
use rapier::parry::bounding_volume::BoundingVolume;
use rapier::parry::shape::FeatureId;
use std::collections::HashMap;
use std::sync::RwLock;
/// Pipeline for particle-based fluid simulation.
pub struct FluidsPipeline {
// FIXME: keep these public?
/// The liquid world responsible for the actual fluid simulation.
pub liquid_world: LiquidWorld,
/// Set of couplings between fluids and rigid-bodies.
pub coupling: ColliderCouplingSet,
}
impl FluidsPipeline {
/// Initialize a new pipeline for fluids simulation.
///
/// # Parameters
///
/// - `particle_radius`: the radius of every particle for the fluid simulation.
/// - `smoothing_factor`: the smoothing factor used to compute the SPH kernel radius.
/// The kernel radius will be computed as `particle_radius * smoothing_factor * 2.0.
pub fn new(particle_radius: f32, smoothing_factor: f32) -> Self {
let dfsph: DFSPHSolver = DFSPHSolver::new();
Self {
liquid_world: LiquidWorld::new(dfsph, particle_radius, smoothing_factor),
coupling: ColliderCouplingSet::new(),
}
}
/// Advances the fluid simulation by `dt` milliseconds.
///
/// All the fluid particles will be affected by an acceleration equal to `gravity`.
/// This `step` function may apply forces to some rigid-bodies that interact with fluids.
/// However, it will not integrate these forces. Use the `PhysicsPipeline` for this integration.
pub fn step(
&mut self,
gravity: &Vector<f32>,
dt: f32,
colliders: &ColliderSet,
bodies: &mut RigidBodySet,
) {
self.liquid_world.step_with_coupling(
dt,
gravity,
&mut self.coupling.as_manager_mut(colliders, bodies),
)
}
}
/// The way a collider is coupled to a boundary object.
pub enum ColliderSampling {
/// The collider shape is approximated with the given sample points in local-space.
///
/// It is recommended that those points are separated by a distance smaller or equal to twice
/// the particle radius used to initialize the LiquidWorld.
StaticSampling(Vec<Point<f32>>),
/// The collider shape is approximated by a dynamic set of points automatically computed based on contacts with fluid particles.
DynamicContactSampling,
}
pub struct ColliderCouplingEntry {
pub sampling_method: ColliderSampling,
pub boundary: BoundaryHandle,
features: Vec<FeatureId>,
}
/// Structure managing all the coupling between colliders from rapier with boundaries and fluids from salva.
pub struct ColliderCouplingSet {
/// The hashmap containing mappings from ColliderHandle to ColliderCouplingEntry
pub entries: HashMap<ColliderHandle, ColliderCouplingEntry>,
}
impl ColliderCouplingSet {
/// Create a new collider coupling manager.
pub fn new() -> Self {
Self {
entries: HashMap::new(),
}
}
/// Register a coupling between a boundary and a collider.
/// There can be only up to one coupling between a collider and a boundary object. If a coupling
/// already exists for this collider when calling this function, the handle of the previously coupled
/// boundary is returned.
pub fn register_coupling(
&mut self,
boundary: BoundaryHandle,
collider: ColliderHandle,
sampling_method: ColliderSampling,
) -> Option<BoundaryHandle> {
let old = self.entries.insert(
collider,
ColliderCouplingEntry {
sampling_method,
boundary,
features: Vec::new(),
},
);
old.map(|e| e.boundary)
}
/// Unregister a coupling between a boundary and a collider.
/// Note that this does not remove the boundary itself from the liquid world.
/// Returns the handle of the boundary this collider was coupled with.
pub fn unregister_coupling(&mut self, collider: ColliderHandle) -> Option<BoundaryHandle> {
let deleted = self.entries.remove(&collider);
deleted.map(|e| e.boundary)
}
/// Use this collider coupling set as a coupling manager.
pub fn as_manager_mut<'a>(
&'a mut self,
colliders: &'a ColliderSet,
bodies: &'a mut RigidBodySet,
) -> ColliderCouplingManager {
ColliderCouplingManager {
coupling: self,
colliders,
bodies,
}
}
}
/// A manager for coupling colliders from rapier2d/rapier3D with the boundary
/// objects from salva.
pub struct ColliderCouplingManager<'a> {
coupling: &'a mut ColliderCouplingSet,
colliders: &'a ColliderSet,
bodies: &'a mut RigidBodySet,
}
impl<'a> CouplingManager for ColliderCouplingManager<'a> {
fn update_boundaries(
&mut self,
timestep: &TimestepManager,
h: f32,
particle_radius: f32,
hgrid: &HGrid<HGridEntry>,
fluids: &mut [Fluid],
boundaries: &mut BoundarySet,
) {
let bodies = &self.bodies;
for (collider, coupling) in &mut self.coupling.entries {
if let (Some(collider), Some(boundary)) = (
self.colliders.get(*collider),
boundaries.get_mut(coupling.boundary),
) {
// Update the boundary's ability to receive forces.
let body = collider.parent().and_then(|p| bodies.get(p));
if let Some(body) = body {
if !body.is_dynamic() {
boundary.forces = None;
} else {
boundary.forces = Some(RwLock::new(Vec::new()));
boundary.clear_forces(true);
}
}
// Update positions and velocities.
boundary.positions.clear();
boundary.velocities.clear();
boundary.volumes.clear();
coupling.features.clear();
match &coupling.sampling_method {
ColliderSampling::StaticSampling(points) => {
for pt in points {
boundary.positions.push(collider.position() * pt);
let velocity = body.map(|b| b.velocity_at_point(pt));
boundary
.velocities
.push(velocity.unwrap_or(Vector::zeros()));
}
boundary.volumes.resize(points.len(), na::zero::<f32>());
}
ColliderSampling::DynamicContactSampling => {
let prediction = h * na::convert::<_, f32>(0.5);
let margin = particle_radius * na::convert::<_, f32>(0.1);
let collider_pos = collider.position();
let aabb = collider
.shape()
.compute_aabb(&collider_pos)
.loosened(h + prediction);
for particle in hgrid
.cells_intersecting_aabb(&aabb.mins, &aabb.maxs)
.flat_map(|e| e.1)
{
match particle {
HGridEntry::FluidParticle(fluid_id, particle_id) => {
let fluid = &mut fluids[*fluid_id];
let particle_pos = fluid.positions[*particle_id]
+ fluid.velocities[*particle_id] * timestep.dt();
if aabb.contains_local_point(&particle_pos) {
let (proj, feature) =
collider.shape().project_point_and_get_feature(
&collider_pos,
&particle_pos,
);
let dpt = particle_pos - proj.point;
if let Some((normal, depth)) =
Unit::try_new_and_get(dpt, f32::default_epsilon())
{
if proj.is_inside {
fluid.positions[*particle_id] -=
*normal * (depth + margin);
let vel_err =
normal.dot(&fluid.velocities[*particle_id]);
if vel_err > na::zero::<f32>() {
fluid.velocities[*particle_id] -=
*normal * vel_err;
}
} else if depth > h + prediction {
continue;
}
}
let velocity =
body.map(|b| b.velocity_at_point(&proj.point));
boundary
.velocities
.push(velocity.unwrap_or(Vector::zeros()));
boundary.positions.push(proj.point);
boundary.volumes.push(na::zero::<f32>());
coupling.features.push(feature);
}
}
HGridEntry::BoundaryParticle(..) => {
// Not yet implemented.
}
}
}
}
}
boundary.clear_forces(true);
}
}
}
fn transmit_forces(&mut self, timestep: &TimestepManager, boundaries: &BoundarySet) {
for (collider, coupling) in &self.coupling.entries {
if let (Some(collider), Some(boundary)) = (
self.colliders.get(*collider),
boundaries.get(coupling.boundary),
) {
if boundary.positions.is_empty() {
continue;
}
if let Some(forces) = &boundary.forces {
let forces = forces.read().unwrap();
if let Some(parent) = collider.parent() {
if let Some(body) = self.bodies.get_mut(parent) {
for (pos, force) in
boundary.positions.iter().zip(forces.iter().cloned())
{
body.apply_impulse_at_point(force * timestep.dt(), *pos, true)
}
}
}
}
}
}
}
}