I wrote this code as I am auditing a course "Physically-based and biomechanical modeling of natural phenomena" by Dr. Prusinkiewicz at the University of Calgary (Jan-April 2015).
I would like to implement fluid flow in 3D using the Smoothed Particle Hydro-Dyanmics. You can see the animation and the code for the 2D version (will compile only with L-studio) below:
//Code written by Jey for 2D fluid simulation using the smoothed particle hydrodynamics
//The code will compile only in the recent Version of L-studio
#include <vve.h>
#include <util/parms.h>
#include <util/palette.h>
#include <algorithms/insert.h>
#include <geometry/geometry.h>
#include <vector>
#include <list>
#include <storage/storage.h>
#include <storage/graph.h>
#include <storage/vector.h>
#include <vector>
#include <math.h>
using geometry::Point3d;
using util::norm;
using algorithms::insert;
using std::cout;
using std::endl;
struct particleContent
{
Point3d pos; // Position of the particle
Point3d acc; // acc applied on the particle
Point3d velocity; // Current velocity of the particle
float density; //Density of the particles
float press; //Pressure
};
// Class defining your model
// If you change the name, don't forget to change the last line of the file too
class MyModel : public Model
{
public:
util::Palette palette;
double dt; //Time step
double TwoPi;
double N; //Total number of particles
double M; //Number of columns of particles
double sigma2; //Controls smoothing distance
int mechSteps;
int color1, color2, color3, colorFail, background;
double lineWidth;
std::vector<particleContent> particles;
void readParms()
{
util::Parms parms("view.v");
parms("Main", "Dt", dt);
parms("Main", "TwoPi", TwoPi);
parms("Main", "NumParticles", N);
parms("Main", "NumColumns", M);
parms("Main", "StepsPerIteration", mechSteps);
parms("Main", "sigma2", sigma2);
parms("View", "Color1", color1);
parms("View", "Color2", color2);
parms("View", "Color3", color3);
parms("View", "ColorFail", colorFail);
parms("View", "BackGround", background);
parms("View", "LineWidth", lineWidth);
}
// Here, reread the files when they are modified
void modifiedFiles( const std::set<std::string>& filenames )
{
forall(const std::string& fn, filenames)
{
if(fn == "pal.map")
palette.reread();
else if(fn == "view.v")
readParms();
}
}
MyModel(QObject *parent)
: Model(parent)
, palette("pal.map")
{
readParms();
// Registering the configuration files
registerFile("pal.map");
registerFile("view.v");
cout<<"Value of N "<<N<<endl;
cout<<"Value of M "<<M<<endl;
cout<<"Value of N/M "<<N/M<<endl;
for(int i = 0; i<N/M; i++)
{
for(int j=0; j<M; j++)
{
particleContent p;
p.pos = Point3d(i+(float)-1, 1.0 + (float) j, 0); // Point3d(0.0+ (float) i, 0.1 + (float) j, 0.0);
particles.push_back(p);
}
}
cout<<"Particles size "<<particles.size()<<endl;
}
float kernel(const Point3d &p)
{
float normalized_distance2 = (p.x()*p.x() + p.y()*p.y())/(2.0*sigma2);
if(normalized_distance2>50)
return 0;
else
return (1.0/(TwoPi*sigma2))*exp(-normalized_distance2);
}
Point3d kernel_gradient(const Point3d &p)
{
Point3d temp(p.x(), p.y(),0.0);
return (-temp/sigma2) * kernel(p);
}
float density(int i)
{
float result = 0;
for(int j = 0; j<N; j++)
{
result+=kernel(particles[i].pos-particles[j].pos);
}
return result;
}
float press(int i)
{
return 0.01*pow(particles[i].density - 0.05, 9.0);
}
Point3d accDueToPressure(int i)
{
Point3d result(0.0, 0.0, 0.0);
for(int j=0; j<N; j++)
{
result+=kernel_gradient(particles[j].pos - particles[i].pos) *
1.0*(particles[j].press+particles[i].press)/particles[j].density;
}
return result/particles[i].density;
}
Point3d accDueToViscosity(int i)
{
Point3d result(0.0, 0.0, 0.0);
for(int j=0; j<N; j++)
{
result+=0.1*kernel(particles[j].pos-particles[i].pos)*
(particles[j].velocity - particles[i].velocity)/particles[j].density;
}
return result/particles[i].density;
}
Point3d gravity()
{
return Point3d(0.0, -0.003, 0.0);
}
void step()
{
for(int i = 0; i<mechSteps; i++)
{
// First, mass-sprint simulation
for(int i = 0; i<N; i++)
{
particles[i].density = density(i);
particles[i].press = press(i);
}
for(int i = 0; i<N; i++)
{
particles[i].acc = accDueToPressure(i) + accDueToViscosity(i) + gravity();
}
for(int i = 0; i<N; i++)
{
particles[i].velocity+=particles[i].acc*dt;
if(particles[i].pos.y()<=0)
{
particles[i].velocity.y()*=-0.5;
}
if(particles[i].pos.x()<=0 || particles[i].pos.x()>=40)
{
particles[i].velocity.x()*=-0.5;
}
particles[i].pos+=particles[i].velocity*dt;
}
}
}
void preDraw()
{
util::Palette::Color c = palette.getColor(background);
//glClearColor(c.r()/255.0, c.g()/255.0, c.b()/255.0, 1);
glClearColor(1.0, 1.0, 1.0, 1.0);
glLineWidth(lineWidth);
}
void draw(Viewer* viewer)
{
Point3d pmin, pmax;
particleContent p = particles.front();
pmax = pmin = p.pos;
for(std::vector<particleContent>::iterator it = particles.begin() ; it != particles.end(); ++it)
{
Point3d& pos = it->pos;
if(pos.x() < pmin.x())
pmin.x() = pos.x();
if(pos.x() > pmax.x())
pmax.x() = pos.x();
if(pos.y() < pmin.y())
pmin.y() = pos.y();
if(pos.y() > pmax.y())
pmax.y() = pos.y();
}
viewer->setSceneBoundingBox(Vec(pmin), Vec(pmax));
//for(std::vector<particleContent>::iterator it = particles.begin() ; it != particles.end(); ++it)
for(int i = 0; i<N; i++)
{
palette.useColor(color1);
//drawCircle(it->pos);
drawCircle(particles[i].pos);
}
}
void drawCircle(const Point3d& pos)
{
glBegin(GL_TRIANGLE_FAN);
glVertex3f(pos.x(),pos.y(),pos.z()+0.2);
for(float i=0;i<2*M_PI;i+=0.1){
glVertex3f(0.25*cos(i)+pos.x(),0.25*sin(i)+pos.y(),pos.z()+0.2);
}
glVertex3f(0.25*cos(0)+pos.x(),0.25*sin(0)+pos.y(),pos.z()+0.2);
glEnd();
}
};
// Define 'MyModel' as the model to be used
DEFINE_MODEL(MyModel);
DEFINE_VIEWER(Viewer);