/*
* IndexGrid.java
*
* Copyright 2011 Chris Martin.
*
* This file is part of Flowy.
*
* Flowy is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Flowy is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Flowy. If not, see <http://www.gnu.org/licenses/>.
*
*
*/
package org.cmj.flowy.simulation;
import java.util.List;
import java.util.Random;
import org.cmj.flowy.simulation.math.RectangleF;
import org.cmj.flowy.simulation.math.Vector2;
import org.cmj.flowy.simulation.particles.FluidParticle;
import org.cmj.flowy.simulation.smothingkernel.SKPoly6;
import org.cmj.flowy.simulation.smothingkernel.SKSpiky;
import org.cmj.flowy.simulation.smothingkernel.SKViscosity;
import org.cmj.flowy.simulation.smothingkernel.SmoothingKernel;
/**
* Class representing a smoothed particle hydrodynamic simulation (see
* http://en.wikipedia.org/wiki/Smoothed-particle_hydrodynamics).
* @author Chris Martin
*/
public class SPHSimulation {
/**
* Random number generator. Used to generate random numbers for adding
* Brownian motion
*/
private static final Random RAND = new Random();
/** The Grid to use for particle positions */
private IndexGrid theGrid;
/** The spacing between grid points */
private float theCellSpacing;
/** The maximum x position (in domain space) that a particle may have*/
private final float theMaxWidth;
/** The minimum x position (in domain space) that a particle may have*/
private final float theMinWidth;
/** The maximum y position (in domain space) that a particle may have*/
private final float theMaxHeight;
/** The minimum y position (in domain space) that a particle may have*/
private final float theMinHeight;
/** Smoothing Kernel Used to spread out the mass of each particle */
private final SmoothingKernel theSKGeneral;
/** Smoothing Kernel used to spread out the pressure of each particle */
private final SmoothingKernel theSKPressure;
/** Smoothing Kernel used to spread out the viscosity of each particle */
private final SmoothingKernel theSKViscosity;
/** The viscosity of the fluid to be simulated */
private final float theViscosity;
/**
* Creates a new Simulation. The simulation will take place in continuous
* space defined by the supplied domain and will be mapped into discrete
* space using the supplied cell spacing.
* @param cellSpace The spacing between
* @param domain The continuous space over which the simulation will occur.
* @param viscosity The viscosity of the fluid to be simulated.
*/
public SPHSimulation(final float cellSpace, final RectangleF domain, final float viscosity) {
theCellSpacing = cellSpace;
theViscosity = viscosity;
theGrid = new IndexGrid(cellSpace, domain);
theSKGeneral = new SKPoly6(cellSpace);
theSKPressure = new SKSpiky(cellSpace);
theSKViscosity = new SKViscosity(cellSpace);
// Set the bounds
theMaxWidth = domain.x + domain.width;
theMinWidth = domain.x;
theMaxHeight = domain.y + domain.height;
theMinHeight = domain.y;
}
/**
* Updates the sate of all particles in the simulation.
* @param particles The particles to update.
* @param globalForce The force that acts on all particles
* @param timeStep The amount of simulation time that has passed since the
* simulation was last updated
*/
public void update(final List<FluidParticle> particles, final Vector2 globalForce, final float timeStep) {
theGrid.refresh(particles);
updatePressureAndDensities(particles);
updateForces(particles, globalForce);
updateParticles(particles, timeStep);
checkParticleDistance(particles);
}
/**
* Updates the pressures and densities of all particles in the supplied list.
* @param particles The particles whose densities and pressures are to be updated.
*/
private void updatePressureAndDensities(final List<FluidParticle> particles) {
final Vector2 distance = new Vector2();
for (int i = particles.size() - 1; i >= 0; --i) {
final FluidParticle particle = particles.get(i);
particle.density = 0.0f;
final List<Integer> neighbours = theGrid.GetNeighbourIndex(particle);
for (int j = neighbours.size()-1; j>=0; --j ) {
final int neighbourIndex = neighbours.get(j);
final FluidParticle neighbour = particles.get(neighbourIndex);
distance.x = particle.position.x - neighbour.position.x;
distance.y = particle.position.y - neighbour.position.y;
particle.density += particle.mass * this.theSKGeneral.Calculate(distance);
}
particle.updatePressure();
}
}
/**
* Updates the forces which are acting on all particles in the supplied list.
* @param particles The particles whose forces are to be updated
* @param globalForce A force which acts on all particles.
*/
private void updateForces(final List<FluidParticle> particles, final Vector2 globalForce) {
// Optimise slightly by precomputing a few things needed in our calculations.
final float halfMass = Constants.PARTICLE_MASS / 2.0f;
final float massViscosityProduct = Constants.PARTICLE_MASS * theViscosity;
final Vector2 distance = new Vector2();
for (int i = particles.size()-1; i >= 0 ; --i) {
final FluidParticle particle = particles.get(i);
// Add global force to every particle
particle.force.add(globalForce);
// Get neighbours for each of the particles
final List<Integer> neighbours = theGrid.GetNeighbourIndex(particle);
for (final int neighbourIndex : neighbours) {
// Prevent double tests
if (neighbourIndex > i) {
final FluidParticle neighbour = particles.get(neighbourIndex);
if (neighbour.density > Constants.FLOAT_EPSILON) {
distance.x = particle.position.x - neighbour.position.x;
distance.y = particle.position.y - neighbour.position.y;
// Pressure
float scalar = halfMass * (particle.pressure + neighbour.pressure) / neighbour.density;
Vector2 force = theSKPressure.CalculateGradient(distance);
force.multiply(scalar);
particle.force.subtract(force);
neighbour.force.add(force);
// Viscosity
scalar = massViscosityProduct * theSKViscosity.CalculateLaplacian(distance)
/ neighbour.density;
force = Vector2.Subtract(neighbour.velocity, particle.velocity);
force.multiply(scalar);
particle.force.add(force);
neighbour.force.subtract(force);
}
}
}
// Add some Brownian motion (i.e a small random force).
particle.force.x += 5.0f*(-0.5f + RAND.nextFloat());
particle.force.y += 5.0f*(-0.5f + RAND.nextFloat());
}
}
/**
* Updates the positions and velocities of all particles in the supplied
* list. Any particles which leave the domain space will be clamped back
* into position.
* @param particles The particles whose positions and velocities are to be updated.
* @param timestep The time over which hte positions and velocities should be updated.
*/
private void updateParticles(final List<FluidParticle> particles, final float timestep) {
for (int i=particles.size()-1; i >= 0; --i) {
final FluidParticle particle = particles.get(i);
// Clip positions to domain space
if (particle.position.x < theMinWidth) {
particle.position.x = (theMinWidth + Constants.FLOAT_EPSILON);
}
else if (particle.position.x > theMaxWidth) {
particle.position.x = (theMaxWidth - Constants.FLOAT_EPSILON);
}
if (particle.position.y < theMinHeight) {
particle.position.y = (theMinHeight + Constants.FLOAT_EPSILON);
}
else if (particle.position.y > theMaxHeight) {
particle.position.y = (theMaxHeight - Constants.FLOAT_EPSILON);
}
// Update velocity + position using forces
particle.updatePositionAndVelocity(timestep);
// Reset force
particle.force.x = 0;
particle.force.y = 0;
}
}
/**
* Checks particles to see if they are too close together. Any aprticles
* that are found to be too close together will be moved to a safe distance.
* @param particles Particles to check.
*/
private void checkParticleDistance(final List<FluidParticle> particles) {
final float minDist = 0.5f * theCellSpacing;
final float minDistSq = minDist * minDist;
final Vector2 particleSeparation = new Vector2();
for (int i = particles.size() -1; i >=0; --i) {
final FluidParticle particle = particles.get(i);
final List<Integer> neighbours = theGrid.GetNeighbourIndex(particle);
for (int j=neighbours.size()-1; j>=0; --j) {
final int neighbourIndex = neighbours.get(j);
// Prevent double tests
if (neighbourIndex > i) {
final FluidParticle neighbour = particles.get(neighbourIndex);
particleSeparation.x = neighbour.position.x - particle.position.x;
particleSeparation.y = neighbour.position.y - particle.position.y;
final float particleSeparationSquared = particleSeparation.getLengthSqaured();
if (particleSeparationSquared < minDistSq) {
if (particleSeparationSquared > Constants.FLOAT_EPSILON) {
final float distLen = particleSeparation.getLength();
particleSeparation.multiply(0.5f * (distLen - minDist) / distLen);
neighbour.position.subtract(particleSeparation);
neighbour.oldPosition.subtract(particleSeparation);
particle.position.add(particleSeparation);
particle.oldPosition.add(particleSeparation);
}
else {
final float diff = 0.5f * minDist;
neighbour.position.y -= diff;
neighbour.oldPosition.y -= diff;
particle.position.y += diff;
particle.oldPosition.y += diff;
}
}
}
}
}
}
}