Simple Fluid Simulation added

Makefile

@@ -1,2 +1,9 @@
-itarget:
-	g++ -Wall -std=c++0x -O4 -fopenmp -mfpmath=sse -msse4 -fexpensive-optimizations main.cpp -I. -o test
+targets=SimpleFluidSimulation
+CFLAGS=-Wall -std=c++0x -I.
+CFLAGS+=-O3 -fopenmp -msse3
+
+SimpleFluidSimulation: SimpleFluidSimulation.cpp
+	g++ $(CFLAGS) $< -o $@
+
+SimpleAlgebraTest: SimpleAlgebraTest.cpp
+	g++ $(CFLAGS) $< -o $@

SimpleAlgebraTest.cpp

@@ -0,0 +1,40 @@
+#include <iostream>
+#include <cstdio>
+#include <chrono>
+#include "construct/Construct.h"
+using namespace Construct;
+using namespace std;
+
+int main(int argc, char **argv) {
+  const Vec3 p(1,1,1);
+
+	// Create some simple fields
+  ScalarField c = 1.f;
+  VectorField x = identity();
+
+  // Do some simple algebraic operations
+  auto a1 = x / c; 
+  cout << a1.eval(p).transpose() << endl;
+
+  auto a2 = x / (constant(1.f) + c*c);
+  cout << a2.eval(p).transpose() << endl;
+
+  auto a3 = (c - dot(x, constant(Vec3(0,1,2)))) / (constant(1.f) - c*dot(x,x));
+  cout << a3.eval(p) << endl;
+
+  // And take derivatives of these expressions!
+  auto d1 = grad(a3);
+  cout << d1.eval(p) << endl;
+
+  auto d2 = grad(x);
+  cout << d2.eval(p) << endl;
+
+  // Can not take a spatial derivative of a matrix field though!
+  // (This would generate a third-order tensor, which isn't supported (yet?)
+  // auto d3 = grad(constant(Mat3(0,0,0)))
+
+  // Also, you can not take two derivatives analytically here...
+  // auto d4 = grad(grad(1.f / dot(x,x)))
+
+  return 0;
+}

SimpleFluidSimulation.cpp

@@ -0,0 +1,92 @@
+#include <iostream>
+#include <cstdio>
+#include <chrono>
+#include "construct/Construct.h"
+using namespace Construct;
+using namespace std;
+
+// Output a _very_ crude PPM down the middle slice of the volume
+void render_ppm(const char *path, ScalarField field, VectorField color, Domain domain) {
+	FILE *f = fopen(path, "wb");
+	int W = 512;//domain.res[0] * 2;
+	int H = 512;//domain.res[2] * 2;
+
+	fprintf(f,"P6\n%d %d\n255\n",W,H);
+	for(int y=H-1;y>=0;--y) {
+		for(int x=0;x<W;++x) {
+
+			Vec3 C(0,0,0); // Output color
+			float T=1; // Transmittance
+			const float ds = domain.H[2] * .5f;
+
+      // Ray march through the volume
+			for(float z=domain.bmin[2];z<=domain.bmax[2];z+=ds) {
+				Vec3 X;
+				X[0] = domain.bmin[0] + domain.extent[0] * (float)x / (float)(W-1);
+				X[1] = domain.bmin[1] + domain.extent[1] * (float)y / (float)(H-1);
+				X[2] = z;
+				float rho = field.eval(X);
+				if(rho <= 0) continue;
+				Vec3 col = color.eval(X);
+
+				const float dT = expf(rho * -ds);
+				T *= dT;
+				C += (1.f-dT) * T * col;
+			}
+			// Gamma adjust + convert to the [0,255] range
+			const float gamma = 1.f / 2.f;
+			for(int i=0;i<3;++i)
+				C[i] = C[i] > 1.f ? 255 : powf(C[i], gamma) * 255;
+			fprintf(f, "%c%c%c", (unsigned char)C[0], (unsigned char)C[1], (unsigned char)C[2]);
+		}
+	}
+
+	fclose(f);
+}
+
+// Semi-lagrangian advection by warping space along 
+// characteristic lines in the flow field
+template<typename T>
+Field<T> advect(Field<T> f, VectorField u, ScalarField dt) {
+	return warp(f, identity() - u * dt);
+}
+
+// Signed distance function for a sphere
+ScalarField sphere(Vec3 center, float radius) {
+  return constant(radius) - length(identity() - constant(center));
+}
+
+int main(int argc, char **argv) {
+	const unsigned int R = 128; // Resolution
+  Domain domain(R, R, R, Vec3(-1,-1,-1), Vec3(1,1,1));
+
+  auto density_source = mask(sphere(Vec3(0,-1,0), .3f)) * 1.f;
+	auto density = constant(0.f);
+  auto velocity = constant(Vec3(0,0,0));
+	auto dt = constant(.1f);
+
+	for(int iter=0; iter<1000; ++iter) {
+		//////////////////////////////////////////////////////////	
+		// Advect density using semi-lagrangian advection		
+		density = advect(density, velocity, dt) + dt * density_source;
+		density = writeToGrid(density, constant(0.f), domain);
+
+		// Advect velocity similarly
+		velocity = advect(velocity, velocity, dt);
+    // Add force upward, proportional to density
+    velocity = velocity + dt * density * constant(Vec3(0,1,0));
+    // Div-Free Projection
+		velocity = divFree(velocity, constant(0.f), domain, 100);
+		//////////////////////////////////////////////////////////	
+
+		// Output results
+		char path[256];
+		sprintf(path, "frame.%04d.ppm", iter);
+		ScalarField render_density = density ;
+		VectorField render_color = density * constant(Vec3(1,1,1));
+		render_ppm(path, render_density, render_color, domain);
+
+		cout << "Finished time step " << iter+1 << endl;
+	}
+	return 0;
+}

construct/ConstructField.h

@@ -70,10 +70,14 @@ struct Field {
 
   Field(ConstructFieldNode<T>* node) : node(NodePtr(node)) { }
   Field(NodePtr node) : node(node) { }
- 
+
   //! Evaluates the underlying expression tree
   T eval(const Vec3& x) const 
   { return node->eval(x); }
+  T operator()(const Vec3& x) const
+  { return eval(x); }
+
+  // Return the gradient of this expression
   typename FieldInfo<T>::GradType grad(const Vec3& x) const
   { return node->grad(x); }
 };

construct/ConstructGrid.h

@@ -168,6 +168,7 @@ template<> void ConstructGrid<Vec3>::divFree(ScalarField boundary, int iteration
       skip.set(i,j,k,1);
   }
 
+  // TODO: Generalize this for other boundaries and conditions (Dirichlet, etc)
 	// Set no flux for velocity
 #pragma omp parallel for
     for(int k=0;k<domain.res[2];++k) 
@@ -189,7 +190,8 @@ template<> void ConstructGrid<Vec3>::divFree(ScalarField boundary, int iteration
       divergence.set(i,j,k, D ); // ASSUMED CUBIC CELLS!
     }
 
-    // CG
+    // Conjugate Gradient 
+    // (http://en.wikipedia.org/wiki/Conjugate_gradient_method)
     // r = b - Ax
 #pragma omp parallel for
     for(int k=1;k<domain.res[2]-1;++k) 
@@ -207,7 +209,6 @@ template<> void ConstructGrid<Vec3>::divFree(ScalarField boundary, int iteration
       r.set(i,j,k, skip.get(i,j,k)==1 ? 0. : R);
     }
 
-
     // d = r
 #pragma omp parallel for
     for(int k=1;k<domain.res[2]-1;++k) 
@@ -227,7 +228,7 @@ template<> void ConstructGrid<Vec3>::divFree(ScalarField boundary, int iteration
 
     // delta = deltaNew
     const real eps = 1.e-4;
-    real maxR = 2. * eps;
+    real maxR = 1.f;
     int iter=0;
     while((iter<iterations) && (maxR > eps)) {
       // q = A d
@@ -298,11 +299,13 @@ template<> void ConstructGrid<Vec3>::divFree(ScalarField boundary, int iteration
 
       // Next iteration...
       ++iter;
-			if(iter%10==0)
+#if 0
+      if(iter%10==0)
       {
         using namespace std;
         cout << "Iteration " << iter << " -- Error: " << maxR << endl;
       }
+#endif
     }
 
 #pragma omp parallel for