6 files changed, 176 insertions, 118 deletions

diff --git a/Smoke/Week 2.ncb b/Smoke/Week 2.ncb
index 30727ff..a2ce9f8 100644
--- a/Smoke/Week 2.ncb
+++ b/Smoke/Week 2.ncb
diff --git a/Smoke/Week 2.suo b/Smoke/Week 2.suo
index 0f9fc0c..3f0e660 100644
--- a/Smoke/Week 2.suo
+++ b/Smoke/Week 2.suo
diff --git a/Smoke/fluids.c b/Smoke/fluids.c
index 9ae7d7f..4c96636 100644
--- a/Smoke/fluids.c
+++ b/Smoke/fluids.c
@@ -27,7 +27,8 @@ fftw_real *vx, *vy;             //(vx,vy)   = velocity field at the current mome
 fftw_real *vx0, *vy0;           //(vx0,vy0) = velocity field at the previous moment
 fftw_real *fx, *fy;	            //(fx,fy)   = user-controlled simulation forces, steered with the mouse 
 fftw_real *rho, *rho0;		    	//smoke density at the current (rho) and previous (rho0) moment 
-fftw_real *height_array;         //used for height plot
+fftw_real *height_array;        //used for height plot
+struct point *normal_array;     //used for normal vectors
 int *frame_hist;
 int frame_index = 0;
 rfftwnd_plan plan_rc, plan_cr;  //simulation domain discretization
@@ -81,7 +82,8 @@ fftw_real *init_simulation(int n)
 	plan_cr = rfftw2d_create_plan(n, n, FFTW_COMPLEX_TO_REAL, FFTW_IN_PLACE);
 
   height_array = (fftw_real*) malloc(dim);
-  frame_hist  = (fftw_real*) malloc(dim * n);
+  normal_array = (struct point*) malloc(dim);
+  frame_hist   = (fftw_real*) malloc(dim * n);
 
   for (i = 0; i <= n; i++)
   {
@@ -212,7 +214,107 @@ void set_forces(void)
 }
 
 
-void calculate_hight_plot(void)
+
+void vector_normal(struct point vert1, struct point vert2, struct point vert3, struct point *normal)
+{
+  struct point v1, v2;
+  float length;
+
+  // vector v1
+  v1.x = vert1.x - vert2.x;
+  v1.y = vert1.y - vert2.y;
+  v1.z = vert1.z - vert2.z;
+
+  // vector v2
+  v2.x = vert2.x - vert3.x;
+	v2.y = vert2.y - vert3.y;
+	v2.z = vert2.z - vert3.z;
+
+  // calculate cross produkt
+  normal->x = v1.y * v2.z - v1.z * v2.y;
+  normal->y = v1.z * v2.x - v1.x * v2.z;
+  normal->z = v1.x * v2.y - v1.y * v2.x;
+
+  length = vec_len3f(normal->x, normal->y, normal->z);
+
+  if(length == 0.0f)
+  {
+    length = 1.0f;
+  }
+
+  // normalize
+  normal->x /= length;
+  normal->y /= length;
+  normal->z /= length;
+}
+
+
+void calculate_normal_vectors(void)
+{
+  int i;
+  struct point p1, p2, p3;
+
+	for (i = 0; i < DIM * DIM; i++) 
+	{
+    switch (vis_dataset)
+    {
+		  case DATASET_FORCE:
+        p1.x = fx[i];
+        p1.y = fy[i];
+        p1.z = height_array[i];
+
+        p2.x = fx[i + 1];
+        p2.y = fy[i + 1];
+        p2.z = height_array[i + 1];
+
+        p3.x = fx[i + DIM];
+        p3.y = fy[i + DIM];
+        p3.z = height_array[i + DIM];
+			  break;
+
+		  case DATASET_VEL:
+        p1.x = vx[i];
+        p1.y = vy[i];
+        p1.z = height_array[i];
+
+        p2.x = vx[i + 1];
+        p2.y = vy[i + 1];
+        p2.z = height_array[i + 1];
+
+        p3.x = vx[i + DIM];
+        p3.y = vy[i + DIM];
+        p3.z = height_array[i + DIM];
+			  break;
+
+		  case DATASET_RHO:
+        p1.x = rho[i];
+        p1.y = rho[i];
+        p1.z = height_array[i];
+
+        p2.x = rho[i + 1];
+        p2.y = rho[i + 1];
+        p2.z = height_array[i + 1];
+
+        p3.x = rho[i + DIM];
+        p3.y = rho[i + DIM];
+        p3.z = height_array[i + DIM];
+			  break;
+
+      default:
+		  case DATASET_DIVV:
+		  case DATASET_DIVF:
+        p1.x = 0.0f; p1.y = 0.0f; p1.z = 0.0f;
+        p2.x = 0.0f; p2.y = 0.0f; p2.z = 0.0f;
+        p3.x = 0.0f; p3.y = 0.0f; p3.z = 0.0f;
+			  break;
+	  }
+
+    vector_normal(p1, p2, p3, &normal_array[i]);
+  }
+}
+
+
+void calculate_height_plot(void)
 {
   int i;
 	for (i = 0; i < DIM * DIM; i++) 
@@ -240,23 +342,6 @@ void calculate_hight_plot(void)
 void copy_frame(fftw_real *field)
 {
   memcpy(field, rho, DIM * 2 * (DIM / 2 + 1) * sizeof(fftw_real));
-#if 0
-  switch (vis_dataset)
-  {
-	  case DATASET_FORCE:
-		  break;
-	  case DATASET_VEL:
-		  break;
-	  case DATASET_RHO:
-      //memcpy(frame_hist[frame_index], rho, sizeof(DIM * 2 * (DIM / 2 + 1) * sizeof(fftw_real)));
-      memcpy(frame_hist[frame_index], rho, DIM * 2 * (DIM / 2 + 1) * sizeof(fftw_real));
-      break;
-    default:
-	  case DATASET_DIVV:
-	  case DATASET_DIVF:
-		  break;
-  }
-#endif
 }
 
 //do_one_simulation_step: Do one complete cycle of the simulation:
@@ -272,7 +357,8 @@ void calculate_one_simulation_step(fftw_real *field)
       set_forces();
       solve(DIM, vx, vy, vx0, vy0, visc, dt);
       diffuse_matter(DIM, vx, vy, rho, rho0, dt);
-      calculate_hight_plot();
+      calculate_height_plot();
+      calculate_normal_vectors();
       copy_frame(field);
   }
 }
diff --git a/Smoke/fluids.h b/Smoke/fluids.h
index db487d4..9814e16 100644
--- a/Smoke/fluids.h
+++ b/Smoke/fluids.h
@@ -3,15 +3,16 @@
 
 //--- SIMULATION PARAMETERS ------------------------------------------------------------------------
 
-extern const int DIM;				//size of simulation grid
-extern double dt;				//simulation time step
-extern float visc;				//fluid viscosity
-extern fftw_real *vx, *vy;             //(vx,vy)   = velocity field at the current moment
-extern fftw_real *vx0, *vy0;           //(vx0,vy0) = velocity field at the previous moment
+extern const int DIM;			             	//size of simulation grid
+extern double dt;		                		//simulation time step
+extern float visc;              				//fluid viscosity
+extern fftw_real *vx, *vy;              //(vx,vy)   = velocity field at the current moment
+extern fftw_real *vx0, *vy0;            //(vx0,vy0) = velocity field at the previous moment
 extern fftw_real *fx, *fy;	            //(fx,fy)   = user-controlled simulation forces, steered with the mouse 
-extern fftw_real *rho, *rho0;			//smoke density at the current (rho) and previous (rho0) moment 
-extern rfftwnd_plan plan_rc, plan_cr;  //simulation domain discretization
+extern fftw_real *rho, *rho0;			      //smoke density at the current (rho) and previous (rho0) moment 
+extern rfftwnd_plan plan_rc, plan_cr;   //simulation domain discretization
 extern fftw_real *height_array;         //used for hight plot
+struct point *normal_array;             //used for normal vectors
 extern int *frame_hist;
 extern int frame_index;
 
diff --git a/Smoke/isolines.c b/Smoke/isolines.c
index f29345b..883b4b6 100644
--- a/Smoke/isolines.c
+++ b/Smoke/isolines.c
@@ -12,7 +12,7 @@
 #include "isolines.h"
 
 
-static int isolines_render = FALSE;
+static int isolines_render = TRUE;
 
 static int isolines_num_colors = PALETTE_MAXCOLORS;
 
diff --git a/Smoke/renderer_gl.c b/Smoke/renderer_gl.c
index bcdded0..0728063 100644
--- a/Smoke/renderer_gl.c
+++ b/Smoke/renderer_gl.c
@@ -35,40 +35,6 @@ float xPos  = -298.0f;
 float yPos  = -295.0f;
 float zPos  = -725.0f;
 
-void vector_normal(struct point vert1, struct point vert2, struct point vert3, struct point *normal)
-{
-  struct point v1, v2;
-  float length;
-
-  // vector v1
-  v1.x = vert1.x - vert2.x;
-  v1.y = vert1.y - vert2.y;
-  v1.z = vert1.z - vert2.z;
-
-  // vector v2
-  v2.x = vert2.x - vert3.x;
-	v2.y = vert2.y - vert3.y;
-	v2.z = vert2.z - vert3.z;
-
-  // calculate cross produkt
-  normal->x = v1.y * v2.z - v1.z * v2.y;
-  normal->y = v1.z * v2.x - v1.x * v2.z;
-  normal->z = v1.x * v2.y - v1.y * v2.x;
-
-  //normal->x = fabs(normal->x);
-  //normal->y = fabs(normal->y);
-  //normal->z = fabs(normal->z);
-
-  if(length == 0.0f) {
-    length = 1.0f;
-  }
-
-  // normalize
-  normal->x /= length;
-  normal->y /= length;
-  normal->z /= length;
-}
-
 
 static void render_legend(void)
 {
@@ -155,9 +121,6 @@ static void render_smoke(void)
 	fftw_real  wn = (fftw_real)winWidth / (fftw_real)(DIM + 1);   // Grid cell width
 	fftw_real  hn = (fftw_real)winHeight / (fftw_real)(DIM + 1);  // Grid cell height
   struct color4f color;
-  struct point vert[3];
-  struct point normal;
-  int vert_nr = 0;
 
   glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
 
@@ -172,7 +135,9 @@ static void render_smoke(void)
 		
     color = colormap_get_color(get_dataset(idx));
     glColor4f(color.r, color.g, color.b, color.a);
-	  glVertex3f(px,py, height_array[idx]); // vertex 1
+	  glVertex3f(px, py, height_array[idx]); // vertex 1
+
+    glNormal3f(normal_array[idx].x, normal_array[idx].y, normal_array[idx].z);
 						
 	  for (i = 0; i < DIM - 1; i++) 
 	  {
@@ -183,19 +148,8 @@ static void render_smoke(void)
       color = colormap_get_color(get_dataset(idx));
       glColor4f(color.r, color.g, color.b, color.a);
       glVertex3f(px, py, height_array[idx]); // vertex 2
-      
-      if (vert_nr > 2)
-      {
-        vert[vert_nr % 3].x = px;
-        vert[vert_nr % 3].y = px;
-        vert[vert_nr % 3].z = height_array[idx];
-      } vert_nr++;
 
-      if (vert_nr % 3 == 0)
-      {
-        vector_normal(vert[0], vert[1], vert[2], &normal);
-        glNormal3f(normal.x, normal.y, normal.z);
-      }
+      glNormal3f(normal_array[idx].x, normal_array[idx].y, normal_array[idx].z);
 
 		  px = wn + (fftw_real)(i + 1) * wn;
 		  py = hn + (fftw_real)j * hn;
@@ -204,15 +158,7 @@ static void render_smoke(void)
       glColor4f(color.r, color.g, color.b, color.a);
 		  glVertex3f(px, py, height_array[idx]); // vertex 3
 
-      vert[vert_nr % 3].x = px;
-      vert[vert_nr % 3].y = px;
-      vert[vert_nr % 3].z = height_array[idx];
-
-      if (vert_nr % 3 == 0)
-      {
-        vector_normal(vert[0], vert[1], vert[2], &normal);
-        glNormal3f(normal.x, normal.y, normal.z);
-      } vert_nr++;
+      glNormal3f(normal_array[idx].x, normal_array[idx].y, normal_array[idx].z);
 	  }
 
 	  px = wn + (fftw_real)(DIM - 1) * wn;
@@ -221,13 +167,8 @@ static void render_smoke(void)
 	  color = colormap_get_color(get_dataset(idx));
     glColor4f(color.r, color.g, color.b, color.a);
 	  glVertex3f(px, py, height_array[idx]); // vertex 4
-    vert_nr++;
 
-    if (vert_nr % 3 == 0)
-    {
-      vector_normal(vert[0], vert[1], vert[2], &normal);
-      glNormal3f(normal.x, normal.y, normal.z);
-    }
+    glNormal3f(normal_array[idx].x, normal_array[idx].y, normal_array[idx].z);
 
 	  glEnd();
   }
@@ -382,7 +323,7 @@ static void render_isolines(void)
     float iso_scale;
 
     if (isolines_nr) {
-      iso_scale = fabs(threshold1 - threshold2) / isolines_nr;
+      iso_scale = (float)(fabs(threshold1 - threshold2) / (float)isolines_nr);
     } else {
       iso_scale = 0.0f;
     }
@@ -403,8 +344,9 @@ static void render_isolines(void)
       
       threshold = min(threshold1, threshold2) + count * iso_scale;
 
-      if (isolines_nr == 1) glLineWidth(2.0f);
-      else glLineWidth(2.0f);
+      glDisable(GL_LIGHTING);
+
+      glLineWidth(2.0f);
       glBegin(GL_LINES);
 
       for (i = 0; i < DIM - 1; i++)
@@ -414,7 +356,8 @@ static void render_isolines(void)
           state = 0;
           idx = (j * DIM) + i;
           color = isolines_get_color(get_dataset(idx));
-          glColor4f(color.r, color.g, color.b, color.a);
+          //glColor4f(color.r, color.g, color.b, color.a);
+          glColor3f(1.0f, 0.0f, 0.0f);
 
           v0 = get_dataset(idx + DIM);
           v1 = get_dataset(idx + 1 + DIM);
@@ -434,54 +377,82 @@ static void render_isolines(void)
             case 5:
             case 10:
               x0 = 0;
-              y0 = percentage(fabs(v3 - v0), threshold) * hn;
-              x1 = percentage(fabs(v1 - v0), threshold) * wn;
+              //y0 = percentage(fabs(v3 - v0), threshold) * hn;
+              //x1 = percentage(fabs(v1 - v0), threshold) * wn;
+              y0 = 0.5f * hn;
+              x1 = 0.5f * wn;
+              //y0 = percentage(max(v3, v0) - min(v3, v0), threshold) * hn;
+              //x1 = percentage(max(v1, v0) - min(v1, v0), threshold) * wn;
               y1 = hn;
               if (x0 != y0 != x1 != y1 != 0.0f) {
-                glVertex3i(x_offset + x0, y_offset + y0, 0.0f);
-                glVertex3i(x_offset + x1, y_offset + y1, 0.0f);
+                glVertex3i(x_offset + x0, y_offset + y0, 5.0f);
+                glVertex3i(x_offset + x1, y_offset + y1, 5.0f);
               }
               // no break !!
             case 4:
             case 11:
-              x0 = percentage(fabs(v3 - v2), threshold) * wn;
+              //x0 = percentage(fabs(v3 - v2), threshold) * wn;
+              //x0 = percentage(max(v3, v0) - min(v3, v2), threshold) * wn;
+              x0 = 0.5f * wn;
               y0 = 0;
               x1 = wn;
-              y1 = percentage(fabs(v2 - v1), threshold) * hn;
+              //y1 = percentage(fabs(v2 - v1), threshold) * hn;
+              //y1 = percentage(max(v2, v1) - min(v2, v1), threshold) * hn;
+              y1 = 0.5f * hn;
               break;
             case 6:
             case 9:
-              x0 = percentage(fabs(v3 - v2), threshold) * wn;
+              //x0 = percentage(fabs(v3 - v2), threshold) * wn;
+              //x0 = percentage(max(v3, v2) - min(v3, v2), threshold) * wn;
+              x0 = 0.5f * wn;
               y0 = 0;
-              x1 = percentage(fabs(v1 - v0), threshold) * wn;
+              //x1 = percentage(fabs(v1 - v0), threshold) * wn;
+              //x1 = percentage(max(v1, v0) - min(v1, v0), threshold) * wn;
+              x1 = 0.5f * wn;
               y1 = hn;
               break;
             case 7:
             case 8:
-              x0 = percentage(fabs(v3 - v2), threshold) * wn;
+              //x0 = percentage(fabs(v3 - v2), threshold) * wn;
+              //x0 = percentage(max(v3, v2) - min(v3, v2), threshold) * wn;
+              x0 = 0.5f * wn;
               y0 = 0;
               x1 = 0;
-              y1 = percentage(fabs(v3 - v0), threshold) * hn;
+              //y1 = percentage(fabs(v3 - v0), threshold) * hn;
+              //y1 = percentage(max(v3, v0) - min(v3, v0), threshold) * hn;
+              y1 = 0.5f * hn;
               break;
             case 3:
             case 12:
               x0 = 0;
-              y0 = percentage(fabs(v3 - v0), threshold) * hn;
+              //y0 = percentage(fabs(v3 - v0), threshold) * hn;
+              //y0 = percentage(max(v3, v0) - min(v3, v0), threshold) * hn;
+              y0 = 0.5f * hn;
               x1 = wn;
-              y1 = percentage(fabs(v2 - v1), threshold) * hn;
+              //y1 = percentage(fabs(v2 - v1), threshold) * hn;
+              //y1 = percentage(max(v2, v1) - min(v2, v1), threshold) * hn;
+              y1 = 0.5f * hn;
               break;
             case 2:
             case 13:
-              x0 = percentage(fabs(v1 - v0), threshold) * wn;
+              //x0 = percentage(fabs(v1 - v0), threshold) * wn;
+              //x0 = percentage(max(v1, v0) - min(v1, v0), threshold) * wn;
+              x0 = 0.5f * wn;
               y0 = hn;
               x1 = wn;
-              y1 = percentage(fabs(v2 - v1), threshold) * hn;
+              //y1 = percentage(fabs(v2 - v1), threshold) * hn;
+              //y1 = percentage(max(v2, v1) - min(v2, v1), threshold) * hn;
+              y1 = 0.5f * hn;
               break;
             case 1:
             case 14:
               x0 = 0;
-              y0 = percentage(fabs(v3 - v0), threshold) * hn;
-              x1 = percentage(fabs(v1 - v0), threshold) * wn;
+              //y0 = percentage(fabs(v3 - v0), threshold) * hn;
+              //x1 = percentage(fabs(v1 - v0), threshold) * wn;
+              //y0 = percentage(max(v3, v0) - min(v3, v0), threshold) * hn;
+              //x1 = percentage(max(v1, v0) - min(v1, v0), threshold) * wn;
+              y0 = 0.5f * hn;
+              x1 = 0.5f * wn;
               y1 = hn;
               break;
             case 0: case 15: default:
@@ -492,8 +463,8 @@ static void render_isolines(void)
           // draw line
           if (x0 != y0 != x1 != y1 != 0.0f)
           {
-            glVertex3i(x_offset + x0, y_offset + y0, 0.0f);
-            glVertex3i(x_offset + x1, y_offset + y1, 0.0f);
+            glVertex3i(x_offset + x0, y_offset + y0, 5.0f);
+            glVertex3i(x_offset + x1, y_offset + y1, 5.0f);
           }
 
         }
@@ -502,6 +473,7 @@ static void render_isolines(void)
       glEnd();
     }
   }
+  glEnable(GL_LIGHTING);
 }
 
 
@@ -588,8 +560,7 @@ void render_flowvis(fftw_real *field)
 
 void init_gl(void)
 {
-  //float LightAmbient[]  = { 0.25f, 0.25f,  0.25f, 1.00f }; // Ambient light values
-  float LightAmbient[]  = { 0.0f, 0.0f,  0.0f, 1.00f }; // Ambient light values
+  float LightAmbient[]  = { 0.25f, 0.25f,  0.25f, 1.00f }; // Ambient light values
   float LightDiffuse[]  = { 0.50f, 0.50f,  0.50f, 1.00f }; // Diffuse light values
   float LightPosition[] = { 0.0f,  400.0f, -700,  1.00f }; // Position of the light source