pressurecell diffuse internal bounds

.vscode/settings.json

@@ -55,6 +55,8 @@
         "tracking.h": "c",
         "multigrid_navier_stokes.h": "c",
         "navier_stokes.h": "c",
-        "electrosphere_client_fluid_cache_fluidchunkdata.h": "c"
+        "electrosphere_client_fluid_cache_fluidchunkdata.h": "c",
+        "normalization.h": "c",
+        "bounds.h": "c"
     }
 }

src/main/c/includes/fluid/queue/chunk.h

@@ -69,6 +69,21 @@ typedef struct {
      * Ratio to normalize the density by
      */
     double normalizationRatio;
+
+    /**
+     * The amount of density that was recaptured
+     */
+    double recaptureDensity;
+
+    /**
+     * The density outgoing to neighbors
+     */
+    double outgoingDensity[9];
+
+    /**
+     * The pressure outgoing to neighbors
+     */
+    double outgoingPressure[9];
 } PressureCellData;
 
 /**

src/main/c/includes/fluid/sim/pressurecell/bounds.h

@@ -36,4 +36,9 @@ LIBRARY_API void pressurecell_update_bounds(Environment * environment, Chunk * c
  */
 LIBRARY_API void pressurecell_update_interest(Environment * environment, Chunk * chunk);
 
+/**
+ * Enforces velocity not running into bounds
+*/
+LIBRARY_API void pressurecell_enforce_bounds(Environment * environment, Chunk * chunk);
+
 #endif

src/main/c/src/fluid/sim/pressurecell/bounds.c

@@ -152,6 +152,42 @@ LIBRARY_API void pressurecell_update_bounds(Environment * environment, Chunk * c
     // fluid_pressurecell_set_bounds_legacy(environment,FLUID_PRESSURECELL_DIRECTION_W,chunk->w0[CENTER_LOC]);
 }
 
+
+/**
+ * Enforces velocity not running into bounds
+*/
+LIBRARY_API void pressurecell_enforce_bounds(Environment * environment, Chunk * chunk){
+    int x, y, z;
+    float * u = chunk->uTempCache;
+    float * v = chunk->vTempCache;
+    float * w = chunk->wTempCache;
+    float * b = chunk->bounds[CENTER_LOC];
+    for(x = 1; x < DIM-1; x++){
+        for(y = 1; y < DIM-1; y++){
+            for(z = 1; z < DIM-1; z++){
+                if(u[IX(x,y,z)] > 0 && b[IX(x+1,y,z)] > 0){
+                    u[IX(x,y,z)] = 0;
+                }
+                if(u[IX(x,y,z)] < 0 && b[IX(x-1,y,z)] > 0){
+                    u[IX(x,y,z)] = 0;
+                }
+                if(v[IX(x,y,z)] > 0 && b[IX(x,y+1,z)] > 0){
+                    v[IX(x,y,z)] = 0;
+                }
+                if(v[IX(x,y,z)] < 0 && b[IX(x,y-1,z)] > 0){
+                    v[IX(x,y,z)] = 0;
+                }
+                if(w[IX(x,y,z)] > 0 && b[IX(x,y,z+1)] > 0){
+                    w[IX(x,y,z)] = 0;
+                }
+                if(w[IX(x,y,z)] < 0 && b[IX(x,y,z-1)] > 0){
+                    w[IX(x,y,z)] = 0;
+                }
+            }
+        }
+    }
+}
+
 /**
  * Updates the interest tree for this chunk
  */

src/main/c/src/fluid/sim/pressurecell/density.c

@@ -27,6 +27,7 @@ LIBRARY_API void pressurecell_diffuse_density(Environment * environment, Chunk *
     int x, y, z;
     float * densityArr = chunk->d[CENTER_LOC];
     float * densityTemp = chunk->dTempCache;
+    float * bounds = chunk->bounds[CENTER_LOC];
     float a = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt / (FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING);
     float c = 1+6*a;
     // float residual = 1;
@@ -50,16 +51,26 @@ LIBRARY_API void pressurecell_diffuse_density(Environment * environment, Chunk *
     for(z = 1; z < DIM-1; z++){
         for(y = 1; y < DIM-1; y++){
             for(x = 1; x < DIM-1; x++){
+                if(bounds[IX(x,y,z)] > 0){
+                    continue;
+                }
                 densityTemp[IX(x,y,z)] = densityArr[IX(x,y,z)] +
                 (
-                    densityArr[IX(x,y,z)] * -6 +
+                    densityArr[IX(x,y,z)] * -(
+                        (1.0f - bounds[IX(x-1,y,z)]) +
+                        (1.0f - bounds[IX(x+1,y,z)]) +
+                        (1.0f - bounds[IX(x,y-1,z)]) +
+                        (1.0f - bounds[IX(x,y+1,z)]) +
+                        (1.0f - bounds[IX(x,y,z-1)]) +
+                        (1.0f - bounds[IX(x,y,z+1)])
+                    ) +
                     (
-                        densityArr[IX(x-1,y,z)] +
+                        densityArr[IX(x-1,y,z)] * (1.0f - bounds[IX(x-1,y,z)]) +
-                        densityArr[IX(x+1,y,z)] +
+                        densityArr[IX(x+1,y,z)] * (1.0f - bounds[IX(x+1,y,z)]) +
-                        densityArr[IX(x,y-1,z)] +
+                        densityArr[IX(x,y-1,z)] * (1.0f - bounds[IX(x,y-1,z)]) +
-                        densityArr[IX(x,y+1,z)] +
+                        densityArr[IX(x,y+1,z)] * (1.0f - bounds[IX(x,y+1,z)]) +
-                        densityArr[IX(x,y,z-1)] +
+                        densityArr[IX(x,y,z-1)] * (1.0f - bounds[IX(x,y,z-1)]) +
-                        densityArr[IX(x,y,z+1)]
+                        densityArr[IX(x,y,z+1)] * (1.0f - bounds[IX(x,y,z+1)])
                     )
                 ) * a
                 ;
@@ -86,9 +97,9 @@ LIBRARY_API void pressurecell_advect_density(Environment * environment, Chunk *
     int x, y, z;
     float * densityArr = chunk->d[CENTER_LOC];
     float * densityTemp = chunk->dTempCache;
-    float * u = chunk->u[CENTER_LOC];
+    float * u = chunk->uTempCache;
-    float * v = chunk->v[CENTER_LOC];
+    float * v = chunk->vTempCache;
-    float * w = chunk->w[CENTER_LOC];
+    float * w = chunk->wTempCache;
     int x0, x1, y0, y1, z0, z1;
     float xp, yp, zp;
     float s0, s1, t0, t1, u0, u1;

src/main/c/src/fluid/sim/pressurecell/normalization.c

@@ -1,4 +1,4 @@
-
+#include <math.h>
 
 #include "fluid/sim/pressurecell/normalization.h"
 #include "fluid/queue/chunkmask.h"
@@ -69,6 +69,194 @@ LIBRARY_API void fluid_pressurecell_normalize_chunk(Environment * env, Chunk * c
  * Recaptures density that has flowed into boundaries
  */
 LIBRARY_API void fluid_pressurecell_recapture_density(Environment * env, Chunk * chunk){
-    
+    int x, y;
+    float * dArr = chunk->d[CENTER_LOC];
+    float * uArr = chunk->u[CENTER_LOC];
+    float * vArr = chunk->v[CENTER_LOC];
+    float * wArr = chunk->w[CENTER_LOC];
+    int ghostIndex, adjacentIndex;
+    float overdraw, estimatedLoss, invertedForce;
+    //clear neighbor outgoing values
+    for(int i = 0; i < 9; i++){
+        chunk->pressureCellData.outgoingDensity[i] = 0;
+        chunk->pressureCellData.outgoingPressure[i] = 0;
+    }
+
+    //check +x plane
+    if(chunk->d[CK(2,1,1)] == NULL){
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(DIM-1,x,y);
+                adjacentIndex = IX(DIM-2,x,y);
+                if(uArr[adjacentIndex] > 0 && dArr[adjacentIndex] > 0){
+                    invertedForce = 1.0f - fabs(uArr[adjacentIndex]);
+                    estimatedLoss = dArr[adjacentIndex] / invertedForce;
+                    dArr[adjacentIndex] = dArr[adjacentIndex] + estimatedLoss;
+                    overdraw = dArr[adjacentIndex] - MAX_FLUID_VALUE;
+                    if(overdraw > 0){
+                        chunk->pressureCellData.recaptureDensity += overdraw;
+                        dArr[adjacentIndex] = MAX_FLUID_VALUE;
+                    }
+                }
+            }
+        }
+    } else {
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(x,y,0);
+                adjacentIndex = IX(x,y,1);
+                if(dArr[ghostIndex] > 0){
+                    chunk->pressureCellData.outgoingDensity[CK(2,1,1)] += dArr[ghostIndex];
+                }
+            }
+        }
+    }
+    //check -x plane
+    if(chunk->d[CK(0,1,1)] == NULL){
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(0,x,y);
+                adjacentIndex = IX(1,x,y);
+                if(uArr[adjacentIndex] < 0 && dArr[adjacentIndex] > 0){
+                    invertedForce = 1.0f - fabs(uArr[adjacentIndex]);
+                    estimatedLoss = dArr[adjacentIndex] / invertedForce;
+                    dArr[adjacentIndex] = dArr[adjacentIndex] + estimatedLoss;
+                    overdraw = dArr[adjacentIndex] - MAX_FLUID_VALUE;
+                    if(overdraw > 0){
+                        chunk->pressureCellData.recaptureDensity += overdraw;
+                        dArr[adjacentIndex] = MAX_FLUID_VALUE;
+                    }
+                }
+            }
+        }
+    } else {
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(0,x,y);
+                adjacentIndex = IX(1,x,y);
+                if(dArr[ghostIndex] > 0){
+                    chunk->pressureCellData.outgoingDensity[CK(0,1,1)] += dArr[ghostIndex];
+                }
+            }
+        }
+    }
+
+    //check +y plane
+    if(chunk->d[CK(1,2,1)] == NULL){
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(x,DIM-1,y);
+                adjacentIndex = IX(x,DIM-2,y);
+                if(vArr[adjacentIndex] > 0 && dArr[adjacentIndex] > 0){
+                    invertedForce = 1.0f - fabs(vArr[adjacentIndex]);
+                    estimatedLoss = dArr[adjacentIndex] / invertedForce;
+                    dArr[adjacentIndex] = dArr[adjacentIndex] + estimatedLoss;
+                    overdraw = dArr[adjacentIndex] - MAX_FLUID_VALUE;
+                    if(overdraw > 0){
+                        chunk->pressureCellData.recaptureDensity += overdraw;
+                        dArr[adjacentIndex] = MAX_FLUID_VALUE;
+                    }
+                }
+            }
+        }
+    } else {
+        for(x = 1; x < DIM-1; x++){
+            ghostIndex = IX(x,DIM-1,y);
+            adjacentIndex = IX(x,DIM-2,y);
+            for(y = 1; y < DIM-1; y++){
+                if(dArr[ghostIndex] > 0){
+                    chunk->pressureCellData.outgoingDensity[CK(1,2,1)] += dArr[ghostIndex];
+                }
+            }
+        }
+    }
+    //check -y plane
+    if(chunk->d[CK(1,0,1)] == NULL){
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(x,0,y);
+                adjacentIndex = IX(x,1,y);
+                if(vArr[adjacentIndex] < 0 && dArr[adjacentIndex] > 0){
+                    invertedForce = 1.0f - fabs(vArr[adjacentIndex]);
+                    estimatedLoss = dArr[adjacentIndex] / invertedForce;
+                    dArr[adjacentIndex] = dArr[adjacentIndex] + estimatedLoss;
+                    overdraw = dArr[adjacentIndex] - MAX_FLUID_VALUE;
+                    if(overdraw > 0){
+                        chunk->pressureCellData.recaptureDensity += overdraw;
+                        dArr[adjacentIndex] = MAX_FLUID_VALUE;
+                    }
+                }
+            }
+        }
+    } else {
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(x,0,y);
+                adjacentIndex = IX(x,1,y);
+                if(dArr[ghostIndex] > 0){
+                    chunk->pressureCellData.outgoingDensity[CK(1,0,1)] += dArr[ghostIndex];
+                }
+            }
+        }
+    }
+
+    //check +z plane
+    if(chunk->d[CK(1,1,2)] == NULL){
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(x,y,DIM-1);
+                adjacentIndex = IX(x,y,DIM-2);
+                if(wArr[adjacentIndex] > 0 && dArr[adjacentIndex] > 0){
+                    invertedForce = 1.0f - fabs(wArr[adjacentIndex]);
+                    estimatedLoss = dArr[adjacentIndex] / invertedForce;
+                    dArr[adjacentIndex] = dArr[adjacentIndex] + estimatedLoss;
+                    overdraw = dArr[adjacentIndex] - MAX_FLUID_VALUE;
+                    if(overdraw > 0){
+                        chunk->pressureCellData.recaptureDensity += overdraw;
+                        dArr[adjacentIndex] = MAX_FLUID_VALUE;
+                    }
+                }
+            }
+        }
+    } else {
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(x,y,DIM-1);
+                adjacentIndex = IX(x,y,DIM-2);
+                if(dArr[ghostIndex] > 0){
+                    chunk->pressureCellData.outgoingDensity[CK(1,1,2)] += dArr[ghostIndex];
+                }
+            }
+        }
+    }
+    //check -z plane
+    if(chunk->d[CK(1,1,0)] == NULL){
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(x,y,0);
+                adjacentIndex = IX(x,y,1);
+                if(wArr[adjacentIndex] < 0 && dArr[adjacentIndex] > 0){
+                    invertedForce = 1.0f - fabs(wArr[adjacentIndex]);
+                    estimatedLoss = dArr[adjacentIndex] / invertedForce;
+                    dArr[adjacentIndex] = dArr[adjacentIndex] + estimatedLoss;
+                    overdraw = dArr[adjacentIndex] - MAX_FLUID_VALUE;
+                    if(overdraw > 0){
+                        chunk->pressureCellData.recaptureDensity += overdraw;
+                        dArr[adjacentIndex] = MAX_FLUID_VALUE;
+                    }
+                }
+            }
+        }
+    } else {
+        for(x = 1; x < DIM-1; x++){
+            for(y = 1; y < DIM-1; y++){
+                ghostIndex = IX(x,y,0);
+                adjacentIndex = IX(x,y,1);
+                if(dArr[ghostIndex] > 0){
+                    chunk->pressureCellData.outgoingDensity[CK(1,1,0)] += dArr[ghostIndex];
+                }
+            }
+        }
+    }
 }
 

src/main/c/src/fluid/sim/pressurecell/pressurecell.c

@@ -107,6 +107,12 @@ LIBRARY_API void fluid_pressurecell_simulate(
         pressurecell_advect_velocity(environment,currentChunk);
     }
 
+    for(int i = 0; i < numChunks; i++){
+        Chunk * currentChunk = chunks[i];
+        //uTemp->div
+        pressurecell_enforce_bounds(environment,currentChunk);
+    }
+
     // printf("approx div\n");
     // fflush(stdout);
     for(int i = 0; i < numChunks; i++){
@@ -153,6 +159,7 @@ LIBRARY_API void fluid_pressurecell_simulate(
     for(int i = 0; i < numChunks; i++){
         Chunk * currentChunk = chunks[i];
         fluid_pressurecell_normalize_chunk(environment,currentChunk);
+        fluid_pressurecell_recapture_density(environment,currentChunk);
         pressurecell_copy_for_next_frame(environment,currentChunk);
         fluid_pressurecell_clearArr(currentChunk->d0[CENTER_LOC]);
         fluid_pressurecell_clearArr(currentChunk->u0[CENTER_LOC]);

src/main/c/src/fluid/sim/pressurecell/velocity.c

@@ -70,6 +70,7 @@ LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk
     float * uTemp = chunk->uTempCache;
     float * vTemp = chunk->vTempCache;
     float * wTemp = chunk->wTempCache;
+    float * bounds = chunk->bounds[CENTER_LOC];
     float a = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * environment->consts.dt / (FLUID_PRESSURECELL_SPACING * FLUID_PRESSURECELL_SPACING);
     float c = 1+6*a;
     // float residual = 1;
@@ -125,16 +126,26 @@ LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk
     for(z = 1; z < DIM-1; z++){
         for(y = 1; y < DIM-1; y++){
             for(x = 1; x < DIM-1; x++){
+                if(bounds[IX(x,y,z)] > 0){
+                    continue;
+                }
                 uArr[IX(x,y,z)] = uTemp[IX(x,y,z)] +
                 (
-                    uTemp[IX(x,y,z)] * -6 +
+                    uTemp[IX(x,y,z)] * -(
+                        (1.0f - bounds[IX(x-1,y,z)]) +
+                        (1.0f - bounds[IX(x+1,y,z)]) +
+                        (1.0f - bounds[IX(x,y-1,z)]) +
+                        (1.0f - bounds[IX(x,y+1,z)]) +
+                        (1.0f - bounds[IX(x,y,z-1)]) +
+                        (1.0f - bounds[IX(x,y,z+1)])
+                    ) +
                     (
-                        uTemp[IX(x-1,y,z)] +
+                        uTemp[IX(x-1,y,z)] * (1.0f - bounds[IX(x-1,y,z)]) +
-                        uTemp[IX(x+1,y,z)] +
+                        uTemp[IX(x+1,y,z)] * (1.0f - bounds[IX(x+1,y,z)]) +
-                        uTemp[IX(x,y-1,z)] +
+                        uTemp[IX(x,y-1,z)] * (1.0f - bounds[IX(x,y-1,z)]) +
-                        uTemp[IX(x,y+1,z)] +
+                        uTemp[IX(x,y+1,z)] * (1.0f - bounds[IX(x,y+1,z)]) +
-                        uTemp[IX(x,y,z-1)] +
+                        uTemp[IX(x,y,z-1)] * (1.0f - bounds[IX(x,y,z-1)]) +
-                        uTemp[IX(x,y,z+1)]
+                        uTemp[IX(x,y,z+1)] * (1.0f - bounds[IX(x,y,z+1)])
                     )
                 ) * a
                 ;
@@ -145,16 +156,26 @@ LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk
     for(z = 1; z < DIM-1; z++){
         for(y = 1; y < DIM-1; y++){
             for(x = 1; x < DIM-1; x++){
+                if(bounds[IX(x,y,z)] > 0){
+                    continue;
+                }
                 vArr[IX(x,y,z)] = vTemp[IX(x,y,z)] +
                 (
-                    vTemp[IX(x,y,z)] * -6 +
+                    vTemp[IX(x,y,z)] * -(
+                        (1.0f - bounds[IX(x-1,y,z)]) +
+                        (1.0f - bounds[IX(x+1,y,z)]) +
+                        (1.0f - bounds[IX(x,y-1,z)]) +
+                        (1.0f - bounds[IX(x,y+1,z)]) +
+                        (1.0f - bounds[IX(x,y,z-1)]) +
+                        (1.0f - bounds[IX(x,y,z+1)])
+                    ) +
                     (
-                        vTemp[IX(x-1,y,z)] +
+                        vTemp[IX(x-1,y,z)] * (1.0f - bounds[IX(x-1,y,z)]) +
-                        vTemp[IX(x+1,y,z)] +
+                        vTemp[IX(x+1,y,z)] * (1.0f - bounds[IX(x+1,y,z)]) +
-                        vTemp[IX(x,y-1,z)] +
+                        vTemp[IX(x,y-1,z)] * (1.0f - bounds[IX(x,y-1,z)]) +
-                        vTemp[IX(x,y+1,z)] +
+                        vTemp[IX(x,y+1,z)] * (1.0f - bounds[IX(x,y+1,z)]) +
-                        vTemp[IX(x,y,z-1)] +
+                        vTemp[IX(x,y,z-1)] * (1.0f - bounds[IX(x,y,z-1)]) +
-                        vTemp[IX(x,y,z+1)]
+                        vTemp[IX(x,y,z+1)] * (1.0f - bounds[IX(x,y,z+1)])
                     )
                 ) * a
                 ;
@@ -178,16 +199,26 @@ LIBRARY_API void pressurecell_diffuse_velocity(Environment * environment, Chunk
     for(z = 1; z < DIM-1; z++){
         for(y = 1; y < DIM-1; y++){
             for(x = 1; x < DIM-1; x++){
+                if(bounds[IX(x,y,z)] > 0){
+                    continue;
+                }
                 wArr[IX(x,y,z)] = wTemp[IX(x,y,z)] +
                 (
-                    wTemp[IX(x,y,z)] * -6 +
+                    wTemp[IX(x,y,z)] * -(
+                        (1.0f - bounds[IX(x-1,y,z)]) +
+                        (1.0f - bounds[IX(x+1,y,z)]) +
+                        (1.0f - bounds[IX(x,y-1,z)]) +
+                        (1.0f - bounds[IX(x,y+1,z)]) +
+                        (1.0f - bounds[IX(x,y,z-1)]) +
+                        (1.0f - bounds[IX(x,y,z+1)])
+                    ) +
                     (
-                    wTemp[IX(x-1,y,z)] +
+                    wTemp[IX(x-1,y,z)] * (1.0f - bounds[IX(x-1,y,z)]) +
-                    wTemp[IX(x+1,y,z)] +
+                    wTemp[IX(x+1,y,z)] * (1.0f - bounds[IX(x+1,y,z)]) +
-                    wTemp[IX(x,y-1,z)] +
+                    wTemp[IX(x,y-1,z)] * (1.0f - bounds[IX(x,y-1,z)]) +
-                    wTemp[IX(x,y+1,z)] +
+                    wTemp[IX(x,y+1,z)] * (1.0f - bounds[IX(x,y+1,z)]) +
-                    wTemp[IX(x,y,z-1)] +
+                    wTemp[IX(x,y,z-1)] * (1.0f - bounds[IX(x,y,z-1)]) +
-                    wTemp[IX(x,y,z+1)]
+                    wTemp[IX(x,y,z+1)] * (1.0f - bounds[IX(x,y,z+1)])
                     )
                 ) * a
                 ;

src/test/c/fluid/sim/pressurecell/bounds_tests.c

@@ -0,0 +1,152 @@
+#include <math.h>
+
+#include "stb/stb_ds.h"
+
+#include "fluid/queue/boundsolver.h"
+#include "fluid/queue/chunkmask.h"
+#include "fluid/queue/chunk.h"
+#include "fluid/env/environment.h"
+#include "fluid/env/utilities.h"
+#include "fluid/sim/pressurecell/bounds.h"
+#include "fluid/sim/pressurecell/density.h"
+#include "fluid/sim/pressurecell/solver_consts.h"
+#include "math/ode/multigrid.h"
+#include "../../../util/chunk_test_utils.h"
+#include "../../../util/test.h"
+
+/**
+ * Error margin for tests
+ */
+#define FLUID_PRESSURE_CELL_ERROR_MARGIN 0.01f
+
+/**
+ * Number of chunks
+ */
+#define CHUNK_DIM 4
+
+/**
+ * 
+ */
+int fluid_sim_pressurecell_bounds_test1(){
+    printf("fluid_sim_pressurecell_bounds_test1\n");
+    int rVal = 0;
+    Environment * env = fluid_environment_create();
+    Chunk ** queue = NULL;
+    queue = createChunkGrid(env,1,1,1);
+    int chunkCount = arrlen(queue);
+
+
+
+    //setup chunk values
+    float deltaDensity = 0.01f;
+    Chunk * currentChunk = queue[0];
+    currentChunk->vTempCache[IX(1,1,1)] = -0.8;
+    currentChunk->bounds[CENTER_LOC][IX(1,0,1)] = 1.0f;
+
+    //actually simulate
+    pressurecell_enforce_bounds(env,currentChunk);
+
+    //test the result
+    float expected, actual;
+
+    //
+    // cell that originall had values
+    //
+    expected = 0;
+    actual = currentChunk->vTempCache[IX(1,1,1)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to enforce advection bounds!  expected: %f    actual: %f  \n");
+    }
+
+    return rVal;
+}
+
+/**
+ * 
+ */
+int fluid_sim_pressurecell_bounds_test2(){
+    printf("fluid_sim_pressurecell_bounds_test2\n");
+    int rVal = 0;
+    Environment * env = fluid_environment_create();
+    Chunk ** queue = NULL;
+    queue = createChunkGrid(env,1,1,1);
+    int chunkCount = arrlen(queue);
+
+
+
+    //setup chunk values
+    float deltaDensity = 0.01f;
+    Chunk * currentChunk = queue[0];
+    currentChunk->vTempCache[IX(1,1,1)] = -0.8;
+    currentChunk->dTempCache[IX(1,1,1)] = MAX_FLUID_VALUE;
+    currentChunk->bounds[CENTER_LOC][IX(1,0,1)] = 1.0f;
+
+    //actually simulate
+    pressurecell_enforce_bounds(env,currentChunk);
+    pressurecell_advect_density(env,currentChunk);
+
+    //test the result
+    float expected, actual;
+
+    //
+    // cell that originall had values
+    //
+    expected = MAX_FLUID_VALUE;
+    actual = chunk_sum_density(currentChunk);
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to enforce advection bounds!  expected: %f    actual: %f  \n");
+    }
+
+    return rVal;
+}
+
+/**
+ * 
+ */
+int fluid_sim_pressurecell_bounds_test3(){
+    printf("fluid_sim_pressurecell_bounds_test3\n");
+    int rVal = 0;
+    Environment * env = fluid_environment_create();
+    Chunk ** queue = NULL;
+    queue = createChunkGrid(env,1,1,1);
+    int chunkCount = arrlen(queue);
+
+
+
+    //setup chunk values
+    float deltaDensity = 0.01f;
+    Chunk * currentChunk = queue[0];
+    currentChunk->vTempCache[IX(1,1,1)] = -0.8;
+    currentChunk->dTempCache[IX(1,1,1)] = MAX_FLUID_VALUE;
+    currentChunk->bounds[CENTER_LOC][IX(1,0,1)] = 1.0f;
+
+    //actually simulate
+    pressurecell_advect_density(env,currentChunk);
+
+    //test the result
+    float expected, actual;
+
+    //
+    // cell that originall had values
+    //
+    expected = MAX_FLUID_VALUE - MAX_FLUID_VALUE * env->consts.dt;
+    actual = chunk_sum_density(currentChunk);
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to enforce advection bounds!  expected: %f    actual: %f  \n");
+    }
+
+    return rVal;
+}
+
+/**
+ * Testing bounds logic
+ */
+int fluid_sim_pressurecell_bounds_tests(int argc, char **argv){
+    int rVal = 0;
+
+    rVal += fluid_sim_pressurecell_bounds_test1();
+    rVal += fluid_sim_pressurecell_bounds_test2();
+    rVal += fluid_sim_pressurecell_bounds_test3();
+
+    return rVal;
+}

src/test/c/fluid/sim/pressurecell/diffuse_tests.c

@@ -114,9 +114,13 @@ int fluid_sim_pressurecell_diffuse_test2(){
 
     //setup chunk values
     Chunk * currentChunk = queue[0];
-    currentChunk->u[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    float * uTemp = currentChunk->uTempCache;
-    currentChunk->v[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    float * vTemp = currentChunk->vTempCache;
-    currentChunk->w[CENTER_LOC][IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    float * wTemp = currentChunk->wTempCache;
+    float * uArr = currentChunk->u0[CENTER_LOC];
+    uTemp[IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    vTemp[IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    wTemp[IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
 
     //actually simulate
     pressurecell_diffuse_velocity(env,currentChunk);
@@ -127,7 +131,7 @@ int fluid_sim_pressurecell_diffuse_test2(){
     // cell that originall had values
     //
     expected = MAX_FLUID_VALUE - FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 6 * MAX_FLUID_VALUE * env->consts.dt;
-    actual = currentChunk->uTempCache[IX(4,4,4)];
+    actual = uArr[IX(4,4,4)];
     if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
         rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,4)!  expected: %f    actual: %f  \n");
     }
@@ -137,37 +141,370 @@ int fluid_sim_pressurecell_diffuse_test2(){
     // neighbors
     //
     expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
-    actual = currentChunk->uTempCache[IX(3,4,4)];
+    actual = uArr[IX(3,4,4)];
     if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
         rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (3,4,4)!  expected: %f    actual: %f  \n");
     }
 
     expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
-    actual = currentChunk->uTempCache[IX(5,4,4)];
+    actual = uArr[IX(5,4,4)];
     if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
         rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (5,4,4)!  expected: %f    actual: %f  \n");
     }
 
     expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
-    actual = currentChunk->uTempCache[IX(4,3,4)];
+    actual = uArr[IX(4,3,4)];
     if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
         rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,3,4)!  expected: %f    actual: %f  \n");
     }
 
     expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
-    actual = currentChunk->uTempCache[IX(4,5,4)];
+    actual = uArr[IX(4,5,4)];
     if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
         rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,5,4)!  expected: %f    actual: %f  \n");
     }
 
     expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
-    actual = currentChunk->uTempCache[IX(4,4,3)];
+    actual = uArr[IX(4,4,3)];
     if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
         rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,3)!  expected: %f    actual: %f  \n");
     }
 
     expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
-    actual = currentChunk->uTempCache[IX(4,4,5)];
+    actual = uArr[IX(4,4,5)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,5)!  expected: %f    actual: %f  \n");
+    }
+
+    return rVal;
+}
+
+/**
+ * Testing diffusing values with bounds
+ */
+int fluid_sim_pressurecell_diffuse_test3(){
+    printf("fluid_sim_pressurecell_diffuse_test3\n");
+    int rVal = 0;
+    Environment * env = fluid_environment_create();
+    Chunk ** queue = NULL;
+    queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM);
+    int chunkCount = arrlen(queue);
+
+
+
+    //setup chunk values
+    Chunk * currentChunk = queue[0];
+    float * dArr = currentChunk->d[CENTER_LOC];
+    float * dTemp = currentChunk->dTempCache;
+    float * bounds = currentChunk->bounds[CENTER_LOC];
+    dArr[IX(4,4,4)] = MAX_FLUID_VALUE;
+    bounds[IX(4,5,4)] = 1.0f;
+
+    //actually simulate
+    pressurecell_diffuse_density(env,currentChunk);
+
+    //test the result
+    float expected, actual;
+
+    //
+    // cell that originall had values
+    //
+    expected = MAX_FLUID_VALUE - FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 5 * MAX_FLUID_VALUE * env->consts.dt;
+    actual = dTemp[IX(4,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+
+    //
+    // neighbors
+    //
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
+    actual = dTemp[IX(3,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (3,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
+    actual = dTemp[IX(5,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN * env->consts.dt){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (5,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
+    actual = dTemp[IX(4,3,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,3,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = dTemp[IX(4,5,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,5,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
+    actual = dTemp[IX(4,4,3)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,3)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * MAX_FLUID_VALUE * env->consts.dt;
+    actual = dTemp[IX(4,4,5)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,5)!  expected: %f    actual: %f  \n");
+    }
+
+    return rVal;
+}
+
+/**
+ * Testing diffusing values with bounds
+ */
+int fluid_sim_pressurecell_diffuse_test4(){
+    printf("fluid_sim_pressurecell_diffuse_test4\n");
+    int rVal = 0;
+    Environment * env = fluid_environment_create();
+    Chunk ** queue = NULL;
+    queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM);
+    int chunkCount = arrlen(queue);
+
+
+
+    //setup chunk values
+    Chunk * currentChunk = queue[0];
+    float * dArr = currentChunk->d[CENTER_LOC];
+    float * dTemp = currentChunk->dTempCache;
+    float * bounds = currentChunk->bounds[CENTER_LOC];
+    dArr[IX(4,4,4)] = MAX_FLUID_VALUE;
+    bounds[IX(3,4,4)] = 1.0f;
+    bounds[IX(5,4,4)] = 1.0f;
+    bounds[IX(4,3,4)] = 1.0f;
+    bounds[IX(4,5,4)] = 1.0f;
+    bounds[IX(4,4,3)] = 1.0f;
+    bounds[IX(4,4,5)] = 1.0f;
+
+    //actually simulate
+    pressurecell_diffuse_density(env,currentChunk);
+
+    //test the result
+    float expected, actual;
+
+    //
+    // cell that originall had values
+    //
+    expected = MAX_FLUID_VALUE;
+    actual = dTemp[IX(4,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+
+    //
+    // neighbors
+    //
+    expected = 0;
+    actual = dTemp[IX(3,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (3,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = dTemp[IX(5,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN * env->consts.dt){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (5,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = dTemp[IX(4,3,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,3,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = dTemp[IX(4,5,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,5,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = dTemp[IX(4,4,3)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,3)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = dTemp[IX(4,4,5)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse density correctly (4,4,5)!  expected: %f    actual: %f  \n");
+    }
+
+    return rVal;
+}
+
+
+/**
+ * Testing diffusing values
+ */
+int fluid_sim_pressurecell_diffuse_test5(){
+    printf("fluid_sim_pressurecell_diffuse_test5\n");
+    int rVal = 0;
+    Environment * env = fluid_environment_create();
+    Chunk ** queue = NULL;
+    queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM);
+    int chunkCount = arrlen(queue);
+
+
+
+    //setup chunk values
+    Chunk * currentChunk = queue[0];
+    float * uTemp = currentChunk->uTempCache;
+    float * vTemp = currentChunk->vTempCache;
+    float * wTemp = currentChunk->wTempCache;
+    float * uArr = currentChunk->u0[CENTER_LOC];
+    float * bounds = currentChunk->bounds[CENTER_LOC];
+    uTemp[IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    vTemp[IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    wTemp[IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    bounds[IX(5,4,4)] = 1.0f;
+
+    //actually simulate
+    pressurecell_diffuse_velocity(env,currentChunk);
+
+    //test the result
+    float expected, actual;
+    //
+    // cell that originall had values
+    //
+    expected = FLUID_PRESSURECELL_MAX_VELOCITY - FLUID_PRESSURECELL_DIFFUSION_CONSTANT * 5 * FLUID_PRESSURECELL_MAX_VELOCITY * env->consts.dt;
+    actual = uArr[IX(4,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+
+    //
+    // neighbors
+    //
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * FLUID_PRESSURECELL_MAX_VELOCITY * env->consts.dt;
+    actual = uArr[IX(3,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (3,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = uArr[IX(5,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (5,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * FLUID_PRESSURECELL_MAX_VELOCITY * env->consts.dt;
+    actual = uArr[IX(4,3,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,3,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * FLUID_PRESSURECELL_MAX_VELOCITY * env->consts.dt;
+    actual = uArr[IX(4,5,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,5,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * FLUID_PRESSURECELL_MAX_VELOCITY * env->consts.dt;
+    actual = uArr[IX(4,4,3)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,3)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = FLUID_PRESSURECELL_DIFFUSION_CONSTANT * FLUID_PRESSURECELL_MAX_VELOCITY * env->consts.dt;
+    actual = uArr[IX(4,4,5)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,5)!  expected: %f    actual: %f  \n");
+    }
+
+    return rVal;
+}
+
+/**
+ * Testing diffusing values
+ */
+int fluid_sim_pressurecell_diffuse_test6(){
+    printf("fluid_sim_pressurecell_diffuse_test6\n");
+    int rVal = 0;
+    Environment * env = fluid_environment_create();
+    Chunk ** queue = NULL;
+    queue = createChunkGrid(env,CHUNK_DIM,CHUNK_DIM,CHUNK_DIM);
+    int chunkCount = arrlen(queue);
+
+
+
+    //setup chunk values
+    Chunk * currentChunk = queue[0];
+    float * uTemp = currentChunk->uTempCache;
+    float * vTemp = currentChunk->vTempCache;
+    float * wTemp = currentChunk->wTempCache;
+    float * uArr = currentChunk->u0[CENTER_LOC];
+    float * bounds = currentChunk->bounds[CENTER_LOC];
+    uTemp[IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    vTemp[IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    wTemp[IX(4,4,4)] = FLUID_PRESSURECELL_MAX_VELOCITY;
+    bounds[IX(3,4,4)] = 1.0f;
+    bounds[IX(5,4,4)] = 1.0f;
+    bounds[IX(4,3,4)] = 1.0f;
+    bounds[IX(4,5,4)] = 1.0f;
+    bounds[IX(4,4,3)] = 1.0f;
+    bounds[IX(4,4,5)] = 1.0f;
+
+    //actually simulate
+    pressurecell_diffuse_velocity(env,currentChunk);
+
+    //test the result
+    float expected, actual;
+    //
+    // cell that originall had values
+    //
+    expected = FLUID_PRESSURECELL_MAX_VELOCITY;
+    actual = uArr[IX(4,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+
+    //
+    // neighbors
+    //
+    expected = 0;
+    actual = uArr[IX(3,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (3,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = uArr[IX(5,4,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (5,4,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = uArr[IX(4,3,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,3,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = uArr[IX(4,5,4)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,5,4)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = uArr[IX(4,4,3)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,3)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = 0;
+    actual = uArr[IX(4,4,5)];
     if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
         rVal += assertEqualsFloat(expected,actual,"Failed to diffuse velocity correctly (4,4,5)!  expected: %f    actual: %f  \n");
     }
@@ -181,8 +518,12 @@ int fluid_sim_pressurecell_diffuse_test2(){
 int fluid_sim_pressurecell_diffuse_tests(int argc, char **argv){
     int rVal = 0;
 
-    // rVal += fluid_sim_pressurecell_diffuse_test1();
+    rVal += fluid_sim_pressurecell_diffuse_test1();
-    // rVal += fluid_sim_pressurecell_diffuse_test2();
+    rVal += fluid_sim_pressurecell_diffuse_test2();
+    rVal += fluid_sim_pressurecell_diffuse_test3();
+    rVal += fluid_sim_pressurecell_diffuse_test4();
+    rVal += fluid_sim_pressurecell_diffuse_test5();
+    rVal += fluid_sim_pressurecell_diffuse_test6();
 
     return rVal;
 }

src/test/c/fluid/sim/pressurecell/normalization_tests.c

@@ -0,0 +1,81 @@
+#include <math.h>
+
+#include "stb/stb_ds.h"
+
+#include "fluid/queue/boundsolver.h"
+#include "fluid/queue/chunkmask.h"
+#include "fluid/queue/chunk.h"
+#include "fluid/env/environment.h"
+#include "fluid/env/utilities.h"
+#include "fluid/sim/pressurecell/normalization.h"
+#include "fluid/sim/pressurecell/solver_consts.h"
+#include "math/ode/multigrid.h"
+#include "../../../util/chunk_test_utils.h"
+#include "../../../util/test.h"
+
+/**
+ * Error margin for tests
+ */
+#define FLUID_PRESSURE_CELL_ERROR_MARGIN 0.01f
+
+/**
+ * Number of chunks
+ */
+#define CHUNK_DIM 4
+
+/**
+ * Testing normalizing values
+ */
+int fluid_sim_pressurecell_normalization_test1(){
+    printf("fluid_sim_pressurecell_normalization_test1\n");
+    int rVal = 0;
+    Environment * env = fluid_environment_create();
+    Chunk ** queue = NULL;
+    queue = createChunkGrid(env,1,1,1);
+    int chunkCount = arrlen(queue);
+
+
+
+    //setup chunk values
+    float deltaDensity = 0.01f;
+    Chunk * currentChunk = queue[0];
+    currentChunk->v[CENTER_LOC][IX(1,1,1)] = -0.8;
+    currentChunk->d[CENTER_LOC][IX(1,1,1)] = 0.2;
+
+    currentChunk->v[CENTER_LOC][IX(DIM-2,DIM-2,DIM-2)] = 0.8;
+    currentChunk->d[CENTER_LOC][IX(DIM-2,DIM-2,DIM-2)] = 0.2;
+
+    //actually simulate
+    fluid_pressurecell_recapture_density(env,currentChunk);
+
+    //test the result
+    float expected, actual;
+
+    //
+    // cell that originall had values
+    //
+    expected = MAX_FLUID_VALUE;
+    actual = currentChunk->d[CENTER_LOC][IX(1,1,1)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to recapture density into (1,1,1)!  expected: %f    actual: %f  \n");
+    }
+
+    expected = MAX_FLUID_VALUE;
+    actual = currentChunk->d[CENTER_LOC][IX(DIM-2,DIM-2,DIM-2)];
+    if(fabs(expected - actual) > FLUID_PRESSURE_CELL_ERROR_MARGIN){
+        rVal += assertEqualsFloat(expected,actual,"Failed to recapture density into (DIM-2,DIM-2,DIM-2)!  expected: %f    actual: %f  \n");
+    }
+
+    return rVal;
+}
+
+/**
+ * Testing normalization
+ */
+int fluid_sim_pressurecell_normalization_tests(int argc, char **argv){
+    int rVal = 0;
+
+    rVal += fluid_sim_pressurecell_normalization_test1();
+
+    return rVal;
+}