Hoisted step funcs working with single chunk

src/main/c/compile.sh

@@ -66,7 +66,7 @@ gcc $COMPILE_FLAGS -I"$BASE_INCLUDE_DIR" -I"$OS_INCLUDE_DIR" $INPUT_FILES -o $OU
 #compile shared object file
 OUTPUT_FILE="libfluidsim$LIB_ENDING"
 COMPILE_FLAGS="-shared"
-INPUT_FILES="fluidsim.o densitystep.o velocitystep.o chunkmask.o"
+INPUT_FILES="densitystep.o velocitystep.o chunkmask.o"
 gcc $COMPILE_FLAGS $INPUT_FILES -o $OUTPUT_FILE
 
 #move to resources

src/main/c/densitystep.c

@@ -6,24 +6,28 @@
 #include "includes/utilities.h"
 #include "includes/chunkmask.h"
 
+void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt);
+
 /*
  * Class:     electrosphere_FluidSim
  * Method:    addDensity
  * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;F)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_addDensity
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jfloat){
+  (JNIEnv * env,
-
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jobjectArray jrx,
+  jobject x0,
+  jfloat dt){
+    int i;
+    int size=N*N*N;
+    float * x = GET_ARR(env,jrx,CENTER_LOC);
+    float * s = (*env)->GetDirectBufferAddress(env,x0);
+    for(i=0; i<size; i++){
+        x[i] += dt*s[i];
     }
-
-/*
- * Class:     electrosphere_FluidSim
- * Method:    setupDiffuseDensity
- * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
- */
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupDiffuseDensity
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat){
-
 }
 
 /*
@@ -32,8 +36,52 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupDiffuseDensity
  * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat){
+  (JNIEnv * env,
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jobjectArray jrx,
+  jobject jrx0,
+  jobjectArray jru,
+  jobjectArray jrv,
+  jobjectArray jrw,
+  jfloat DIFFUSION_CONST,
+  jfloat VISCOSITY_CONST,
+  jfloat dt){
+    float a=dt*DIFFUSION_CONST*N*N*N;
+    float c=1+6*a;
+    int i, j, k, l, m;
+    float * x = GET_ARR(env,jrx,CENTER_LOC);
+    float * x0 = (*env)->GetDirectBufferAddress(env,jrx0);
     
+    __m256 aScalar = _mm256_set1_ps(a);
+    __m256 cScalar = _mm256_set1_ps(c);
+
+    //transform u direction
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&x[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x0[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&x[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    x[IX(i,j,k)] = (x0[IX(i,j,k)] + a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]+x[IX(i,j,k-1)]+x[IX(i,j,k+1)]))/c;
+                }
+            }
+        }
+    }
 }
 
 /*
@@ -42,6 +90,139 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
  * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectDensity
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat){
+  (JNIEnv * env,
-
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jobjectArray jrx,
+  jobject jrx0,
+  jobjectArray jru,
+  jobjectArray jrv,
+  jobjectArray jrw,
+  jfloat DIFFUSION_CONST,
+  jfloat VISCOSITY_CONST,
+  jfloat dt){
+    float * x0 = (*env)->GetDirectBufferAddress(env,jrx0);
+    advectDensity(env,chunk_mask,N,3,jrx,x0,GET_ARR(env,jru,CENTER_LOC),GET_ARR(env,jrv,CENTER_LOC),GET_ARR(env,jrw,CENTER_LOC),dt);
+}
+
+void advectDensity(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt){
+    int i, j, k, i0, j0, k0, i1, j1, k1;
+    int m,n,o;
+    float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz;
+    
+    dtx=dty=dtz=dt*N;
+
+    float * d = GET_ARR(env,jrd,CENTER_LOC);
+
+    float * sampleArr = d0;
+
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            for(i=1; i<N-1; i++){
+                //calculate location to pull from
+                x = i-dtx*u[IX(i,j,k)];
+                y = j-dty*v[IX(i,j,k)];
+                z = k-dtz*w[IX(i,j,k)];
+
+                m = n = o = 1;
+
+                if(x < 1){ m -= 1; }
+                if(x >= N-1){ m += 1; }
+                if(y < 1){ n -= 1; }
+                if(y >= N-1){ n += 1; }
+                if(z < 1){ o -= 1; }
+                if(z >= N-1){ o += 1; }
+
+                //If the out of bounds coordinate is in bounds for a neighbor chunk, use that chunk as source instead
+                // if(CK(m,n,o) != CENTER_LOC){
+                //     printf("Looking in border chunk\n");
+                // }
+                // if(x > 16){
+                //     printf("%f %d %d %d\n",m,n,o);
+                // }
+                // if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
+                //     // printf("Hit other chunk\n");
+                //     sampleArr = GET_ARR(env,jrd,CK(m,n,o));
+                //     x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * N;
+                //     y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * N;
+                //     z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * N;
+                // }
+
+                //clamp location within chunk
+                if (x<0.5f) x=0.5f;
+                if (x>N+0.5f) x=N+0.5f;
+                if (y<0.5f) y=0.5f;
+                if (y>N+0.5f) y=N+0.5f;
+                if (z<0.5f) z=0.5f;
+                if (z>N+0.5f) z=N+0.5f;
+
+                //get actual indices
+                i0=(int)x;
+                i1=i0+1;
+                j0=(int)y;
+                j1=j0+1;
+                k0=(int)z;
+                k1=k0+1;
+
+                //calculate percentage of each index
+                s1 = x-i0;
+                s0 = 1-s1;
+                t1 = y-j0;
+                t0 = 1-t1;
+                u1 = z-k0;
+                u0 = 1-u1;
+
+                if(i0 >= N){
+                    i0 = N - 1;
+                }
+                // if(i0 < 0){
+                //     i0 = 0;
+                // }
+                if(j0 >= N){
+                    j0 = N - 1;
+                }
+                // if(j0 < 0){
+                //     j0 = 0;
+                // }
+                if(k0 >= N){
+                    k0 = N - 1;
+                }
+                // if(k0 < 0){
+                //     k0 = 0;
+                // }
+                if(i1 >= N){
+                    i1 = N - 1;
+                }
+                // if(i1 < 0){
+                //     i1 = 0;
+                // }
+                if(j1 >= N){
+                    j1 = N - 1;
+                }
+                // if(j1 < 0){
+                //     j1 = 0;
+                // }
+                if(k1 >= N){
+                    k1 = N - 1;
+                }
+                // if(k1 < 0){
+                //     k1 = 0;
+                // }
+                d[IX(i,j,k)] = 
+                s0*(
+                    t0*u0*sampleArr[IX(i0,j0,k0)]+
+                    t1*u0*sampleArr[IX(i0,j1,k0)]+
+                    t0*u1*sampleArr[IX(i0,j0,k1)]+
+                    t1*u1*sampleArr[IX(i0,j1,k1)]
+                )+
+                s1*(
+                    t0*u0*sampleArr[IX(i1,j0,k0)]+
+                    t1*u0*sampleArr[IX(i1,j1,k0)]+
+                    t0*u1*sampleArr[IX(i1,j0,k1)]+
+                    t1*u1*sampleArr[IX(i1,j1,k1)]
+                );
+            }
+        }
+    }
 }

src/main/c/includes/electrosphere_FluidSim.h

@@ -20,10 +20,10 @@ extern "C" {
 /*
  * Class:     electrosphere_FluidSim
  * Method:    simulate
- * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
+ * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_simulate
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
@@ -36,58 +36,50 @@ JNIEXPORT jint JNICALL Java_electrosphere_FluidSim_calculateChunkMask
 /*
  * Class:     electrosphere_FluidSim
  * Method:    addSourceToVectors
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_addSourceToVectors
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
-
-/*
- * Class:     electrosphere_FluidSim
- * Method:    setupVectorDiffuse
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
- */
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupVectorDiffuse
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
  * Method:    solveVectorDiffuse
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveVectorDiffuse
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
  * Method:    setupProjection
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
  * Method:    solveProjection
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveProjection
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
  * Method:    finalizeProjection
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
  * Method:    advectVectors
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
+ * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectVectors
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat);
+  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
 /*
  * Class:     electrosphere_FluidSim
@@ -97,14 +89,6 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectVectors
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_addDensity
   (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jfloat);
 
-/*
- * Class:     electrosphere_FluidSim
- * Method:    setupDiffuseDensity
- * Signature: (II[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;FFF)V
- */
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupDiffuseDensity
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
-
 /*
  * Class:     electrosphere_FluidSim
  * Method:    solveDiffuseDensity
@@ -121,6 +105,14 @@ JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveDiffuseDensity
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectDensity
   (JNIEnv *, jobject, jint, jint, jobjectArray, jobject, jobjectArray, jobjectArray, jobjectArray, jfloat, jfloat, jfloat);
 
+/*
+ * Class:     electrosphere_FluidSim
+ * Method:    setBoundsToNeighbors
+ * Signature: (III[Ljava/nio/ByteBuffer;)V
+ */
+JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setBoundsToNeighbors
+  (JNIEnv *, jobject, jint, jint, jint, jobjectArray);
+
 #ifdef __cplusplus
 }
 #endif

src/main/c/velocitystep.c

@@ -7,73 +7,619 @@
 #include "includes/chunkmask.h"
 
 
+#define BOUND_NO_DIR 0
+#define BOUND_DIR_U 1
+#define BOUND_DIR_V 2
+#define BOUND_DIR_W 3
+
+#define SET_BOUND_IGNORE 0
+#define SET_BOUND_USE_NEIGHBOR 1
+
+void add_source(int N, float * x, float * s, float dt);
+void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt);
+
 
 /*
- * Class:     electrosphere_FluidSim
+ * Adds force to all vectors
- * Method:    addSourceToVectors
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_addSourceToVectors
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat){
+  (JNIEnv * env,
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jobjectArray jru,
+  jobjectArray jrv,
+  jobjectArray jrw,
+  jobjectArray jru0,
+  jobjectArray jrv0,
+  jobjectArray jrw0,
+  jfloat DIFFUSION_CONST,
+  jfloat VISCOSITY_CONST,
+  jfloat dt){
+    add_source(N,GET_ARR(env,jru,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),dt);
+    add_source(N,GET_ARR(env,jrv,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),dt);
+    add_source(N,GET_ARR(env,jrw,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
+}
 
+void add_source(int N, float * x, float * s, float dt){
+	int i;
+    int size=N*N*N;
+	for(i=0; i<size; i++){
+        if(s[i] > 0){
+            printf("%f\n",s[i]);
+        }
+        x[i] += dt*s[i];
+    }
 }
 
 /*
- * Class:     electrosphere_FluidSim
+ * Solves vector diffusion along all axis
- * Method:    setupVectorDiffuse
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
- */
-JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupVectorDiffuse
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat){
-
-}
-
-/*
- * Class:     electrosphere_FluidSim
- * Method:    solveVectorDiffuse
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveVectorDiffuse
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat){
+  (JNIEnv * env,
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jobjectArray jru,
+  jobjectArray jrv,
+  jobjectArray jrw,
+  jobjectArray jru0,
+  jobjectArray jrv0,
+  jobjectArray jrw0,
+  jfloat DIFFUSION_CONST,
+  jfloat VISCOSITY_CONST,
+  jfloat dt){
+    float a=dt*DIFFUSION_CONST*N*N*N;
+    float c=1+6*a;
+    int i, j, k, l, m;
+    float * u = GET_ARR(env,jru,CENTER_LOC);
+    float * v = GET_ARR(env,jrv,CENTER_LOC);
+    float * w = GET_ARR(env,jrw,CENTER_LOC);
+    float * u0 = GET_ARR(env,jru0,CENTER_LOC);
+    float * v0 = GET_ARR(env,jrv0,CENTER_LOC);
+    float * w0 = GET_ARR(env,jrw0,CENTER_LOC);
     
+    __m256 aScalar = _mm256_set1_ps(a);
+    __m256 cScalar = _mm256_set1_ps(c);
+
+    //transform u direction
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&u[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&u0[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&u[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    u[IX(i,j,k)] = (u0[IX(i,j,k)] + a*(u[IX(i-1,j,k)]+u[IX(i+1,j,k)]+u[IX(i,j-1,k)]+u[IX(i,j+1,k)]+u[IX(i,j,k-1)]+u[IX(i,j,k+1)]))/c;
+                }
+            }
+        }
+    }
+
+    //transform v direction
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&v[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&v0[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&v[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    v[IX(i,j,k)] = (v0[IX(i,j,k)] + a*(v[IX(i-1,j,k)]+v[IX(i+1,j,k)]+v[IX(i,j-1,k)]+v[IX(i,j+1,k)]+v[IX(i,j,k-1)]+v[IX(i,j,k+1)]))/c;
+                }
+            }
+        }
+    }
+
+    //transform w direction
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&w[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&w0[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&w[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    w[IX(i,j,k)] = (w0[IX(i,j,k)] + a*(w[IX(i-1,j,k)]+w[IX(i+1,j,k)]+w[IX(i,j-1,k)]+w[IX(i,j+1,k)]+w[IX(i,j,k-1)]+w[IX(i,j,k+1)]))/c;
+                }
+            }
+        }
+    }
 }
 
 /*
- * Class:     electrosphere_FluidSim
+ * Sets up a projection system of equations
- * Method:    setupProjection
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setupProjection
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat){
+  (JNIEnv * env,
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jobjectArray jru,
+  jobjectArray jrv,
+  jobjectArray jrw,
+  jobjectArray jru0,
+  jobjectArray jrv0,
+  jobjectArray jrw0,
+  jfloat DIFFUSION_CONST,
+  jfloat VISCOSITY_CONST,
+  jfloat dt){
+    int i, j, k;
 
+    __m256 nVector = _mm256_set1_ps(N);
+    __m256 constScalar = _mm256_set1_ps(-1.0/3.0);
+    __m256 zeroVec = _mm256_set1_ps(0);
+    __m256 vector, vector2, vector3;
+
+    float * u = GET_ARR(env,jru,CENTER_LOC);
+    float * v = GET_ARR(env,jrv,CENTER_LOC);
+    float * w = GET_ARR(env,jrw,CENTER_LOC);
+
+    float * p = GET_ARR(env,jru0,CENTER_LOC);
+    float * div = GET_ARR(env,jrv0,CENTER_LOC);
+
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            i = 1;
+            //
+            //lower
+            //
+            //first part
+            vector = _mm256_loadu_ps(&u[IX(i+1,j,k)]);
+            vector = _mm256_sub_ps(vector,_mm256_loadu_ps(&u[IX(i-1,j,k)]));
+            vector = _mm256_div_ps(vector,nVector);
+            //second part
+            vector2 = _mm256_loadu_ps(&v[IX(i,j+1,k)]);
+            vector2 = _mm256_sub_ps(vector2,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
+            vector2 = _mm256_div_ps(vector2,nVector);
+            //third part
+            vector3 = _mm256_loadu_ps(&w[IX(i,j,k+1)]);
+            vector3 = _mm256_sub_ps(vector3,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
+            vector3 = _mm256_div_ps(vector3,nVector);
+            //multiply and finalize
+            vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
+            vector = _mm256_mul_ps(vector,constScalar);
+            //store
+            _mm256_storeu_ps(&div[IX(i,j,k)],vector);
+            _mm256_storeu_ps(&p[IX(i,j,k)],zeroVec);
+            i = 9;
+            //
+            //upper
+            //
+            //first part
+            vector = _mm256_loadu_ps(&u[IX(i+1,j,k)]);
+            vector = _mm256_sub_ps(vector,_mm256_loadu_ps(&u[IX(i-1,j,k)]));
+            vector = _mm256_div_ps(vector,nVector);
+            //second part
+            vector2 = _mm256_loadu_ps(&v[IX(i,j+1,k)]);
+            vector2 = _mm256_sub_ps(vector2,_mm256_loadu_ps(&v[IX(i,j-1,k)]));
+            vector2 = _mm256_div_ps(vector2,nVector);
+            //third part
+            vector3 = _mm256_loadu_ps(&w[IX(i,j,k+1)]);
+            vector3 = _mm256_sub_ps(vector3,_mm256_loadu_ps(&w[IX(i,j,k-1)]));
+            vector3 = _mm256_div_ps(vector3,nVector);
+            //multiply and finalize
+            vector = _mm256_add_ps(vector,_mm256_add_ps(vector2,vector3));
+            vector = _mm256_mul_ps(vector,constScalar);
+            //store
+            _mm256_storeu_ps(&div[IX(i,j,k)],vector);
+            _mm256_storeu_ps(&p[IX(i,j,k)],zeroVec);
+        }
+    }
+    
+    // set_bnd(env, chunk_mask, N, BOUND_NO_DIR, div);
+    // set_bnd(env, chunk_mask, N, BOUND_NO_DIR, p);
 }
 
 /*
- * Class:     electrosphere_FluidSim
+ * Solves a projection system of equations
- * Method:    solveProjection
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_solveProjection
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat){
+  (JNIEnv * env,
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jobjectArray jru,
+  jobjectArray jrv,
+  jobjectArray jrw,
+  jobjectArray jru0,
+  jobjectArray jrv0,
+  jobjectArray jrw0,
+  jfloat DIFFUSION_CONST,
+  jfloat VISCOSITY_CONST,
+  jfloat dt){
+    int a = 1;
+    int c = 6;
+    int i, j, k, l, m;
+    __m256 aScalar = _mm256_set1_ps(a);
+    __m256 cScalar = _mm256_set1_ps(c);
 
+    float * x = GET_ARR(env,jru0,CENTER_LOC);
+    float * x0 = GET_ARR(env,jrv0,CENTER_LOC);
+    // update for each cell
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            int n = 0;
+            //solve as much as possible vectorized
+            for(i = 1; i < N-1; i=i+8){
+                __m256 vector = _mm256_loadu_ps(&x[IX(i-1,j,k)]);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i+1,j,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j-1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j+1,k)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k-1)]));
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x[IX(i,j,k+1)]));
+                vector = _mm256_mul_ps(vector,aScalar);
+                vector = _mm256_add_ps(vector,_mm256_loadu_ps(&x0[IX(i,j,k)]));
+                vector = _mm256_div_ps(vector,cScalar);
+                _mm256_storeu_ps(&x[IX(i,j,k)],vector);
+            }
+            //If there is any leftover, perform manual solving
+            if(i>N-1){
+                for(i=i-8; i < N-1; i++){
+                    x[IX(i,j,k)] = (x0[IX(i,j,k)] + a*(x[IX(i-1,j,k)]+x[IX(i+1,j,k)]+x[IX(i,j-1,k)]+x[IX(i,j+1,k)]+x[IX(i,j,k-1)]+x[IX(i,j,k+1)]))/c;
+                }
+            }
+        }
+    }
+    // set_bnd(env, chunk_mask, SET_BOUND_USE_NEIGHBOR, N, b, d, jrd);
+    // set_bnd(env, chunk_mask, N, BOUND_NO_DIR, x);
 }
 
 /*
- * Class:     electrosphere_FluidSim
+ * Finalizes a projection (subtract curl, set bounds, etc)
- * Method:    finalizeProjection
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_finalizeProjection
-  (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat){
+  (JNIEnv * env,
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jobjectArray jru,
+  jobjectArray jrv,
+  jobjectArray jrw,
+  jobjectArray jru0,
+  jobjectArray jrv0,
+  jobjectArray jrw0,
+  jfloat DIFFUSION_CONST,
+  jfloat VISCOSITY_CONST,
+  jfloat dt){
+    int i, j, k;
+    __m256 nVector = _mm256_set1_ps(N);
+    __m256 constScalar = _mm256_set1_ps(0.5f*N);
+    __m256 zeroVec = _mm256_set1_ps(0);
+    __m256 vector, vector2, vector3;
 
+    float * u = GET_ARR(env,jru,CENTER_LOC);
+    float * v = GET_ARR(env,jrv,CENTER_LOC);
+    float * w = GET_ARR(env,jrw,CENTER_LOC);
+
+    float * p = GET_ARR(env,jru0,CENTER_LOC);
+    float * div = GET_ARR(env,jrv0,CENTER_LOC);
+
+    for ( k=1 ; k<N-1 ; k++ ) {
+        for ( j=1 ; j<N-1 ; j++ ) {
+            //
+            //v
+            //
+            //lower
+            vector =  _mm256_loadu_ps(&p[IX(1+1,j,k)]);
+            vector2 = _mm256_loadu_ps(&p[IX(1-1,j,k)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(1,j,k)]),vector);
+            _mm256_storeu_ps(&u[IX(1,j,k)],vector);
+            //upper
+            vector =  _mm256_loadu_ps(&p[IX(9+1,j,k)]);
+            vector2 = _mm256_loadu_ps(&p[IX(9-1,j,k)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&u[IX(9,j,k)]),vector);
+            _mm256_storeu_ps(&u[IX(9,j,k)],vector);
+            //
+            //v
+            //
+            //lower
+            vector =  _mm256_loadu_ps(&p[IX(1,j+1,k)]);
+            vector2 = _mm256_loadu_ps(&p[IX(1,j-1,k)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(1,j,k)]),vector);
+            _mm256_storeu_ps(&v[IX(1,j,k)],vector);
+            //upper
+            vector =  _mm256_loadu_ps(&p[IX(9,j+1,k)]);
+            vector2 = _mm256_loadu_ps(&p[IX(9,j-1,k)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&v[IX(9,j,k)]),vector);
+            _mm256_storeu_ps(&v[IX(9,j,k)],vector);
+            //
+            //w
+            //
+            //lower
+            vector =  _mm256_loadu_ps(&p[IX(1,j,k+1)]);
+            vector2 = _mm256_loadu_ps(&p[IX(1,j,k-1)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(1,j,k)]),vector);
+            _mm256_storeu_ps(&w[IX(1,j,k)],vector);
+            //upper
+            vector =  _mm256_loadu_ps(&p[IX(9,j,k+1)]);
+            vector2 = _mm256_loadu_ps(&p[IX(9,j,k-1)]);
+            vector =  _mm256_sub_ps(vector,vector2);
+            vector = _mm256_mul_ps(vector,constScalar);
+            vector = _mm256_sub_ps(_mm256_loadu_ps(&w[IX(9,j,k)]),vector);
+            _mm256_storeu_ps(&w[IX(9,j,k)],vector);
+        }
+    }
+    
+    // set_bnd_old(env, chunk_mask, N, 1, u);
+    // set_bnd_old(env, chunk_mask, N, 2, v);
+    // set_bnd_old(env, chunk_mask, N, 3, w);
 }
 
 /*
- * Class:     electrosphere_FluidSim
+ * Advects u, v, and w
- * Method:    advectVectors
- * Signature: (II[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;[Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;Ljava/nio/ByteBuffer;FFF)V
  */
 JNIEXPORT void JNICALL Java_electrosphere_FluidSim_advectVectors
-(JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jobjectArray, jobject, jobject, jobject, jfloat, jfloat, jfloat){
+(JNIEnv * env,
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jobjectArray jru,
+  jobjectArray jrv,
+  jobjectArray jrw,
+  jobjectArray jru0,
+  jobjectArray jrv0,
+  jobjectArray jrw0,
+  jfloat DIFFUSION_CONST,
+  jfloat VISCOSITY_CONST,
+  jfloat dt){
+    advect(env,chunk_mask,N,1,jru,GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
+    advect(env,chunk_mask,N,2,jrv,GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
+    advect(env,chunk_mask,N,3,jrw,GET_ARR(env,jrw0,CENTER_LOC),GET_ARR(env,jru0,CENTER_LOC),GET_ARR(env,jrv0,CENTER_LOC),GET_ARR(env,jrw0,CENTER_LOC),dt);
+}
+
+
+void advect(JNIEnv * env, uint32_t chunk_mask, int N, int b, jobjectArray jrd, float * d0, float * u, float * v, float * w, float dt){
+    int i, j, k, i0, j0, k0, i1, j1, k1;
+    int m,n,o;
+    float x, y, z, s0, t0, s1, t1, u1, u0, dtx,dty,dtz;
+    
+    dtx=dty=dtz=dt*N;
+
+    float * d = GET_ARR(env,jrd,CENTER_LOC);
+
+    float * sampleArr = d0;
+
+    for(k=1; k<N-1; k++){
+        for(j=1; j<N-1; j++){
+            for(i=1; i<N-1; i++){
+                //calculate location to pull from
+                x = i-dtx*u[IX(i,j,k)];
+                y = j-dty*v[IX(i,j,k)];
+                z = k-dtz*w[IX(i,j,k)];
+
+                m = n = o = 1;
+
+                if(x < 1){ m -= 1; }
+                if(x >= N-1){ m += 1; }
+                if(y < 1){ n -= 1; }
+                if(y >= N-1){ n += 1; }
+                if(z < 1){ o -= 1; }
+                if(z >= N-1){ o += 1; }
+
+                //If the out of bounds coordinate is in bounds for a neighbor chunk, use that chunk as source instead
+                // if(CK(m,n,o) != CENTER_LOC){
+                //     printf("Looking in border chunk\n");
+                // }
+                // if(x > 16){
+                //     printf("%f %d %d %d\n",m,n,o);
+                // }
+                // if(CK(m,n,o) != CENTER_LOC && ARR_EXISTS(chunk_mask,m,n,o)){
+                //     // printf("Hit other chunk\n");
+                //     sampleArr = GET_ARR(env,jrd,CK(m,n,o));
+                //     x = x + CHUNK_NORMALIZE_U[CK(m,n,o)] * N;
+                //     y = y + CHUNK_NORMALIZE_V[CK(m,n,o)] * N;
+                //     z = z + CHUNK_NORMALIZE_W[CK(m,n,o)] * N;
+                // }
+
+                //clamp location within chunk
+                if (x<0.5f) x=0.5f;
+                if (x>N+0.5f) x=N+0.5f;
+                if (y<0.5f) y=0.5f;
+                if (y>N+0.5f) y=N+0.5f;
+                if (z<0.5f) z=0.5f;
+                if (z>N+0.5f) z=N+0.5f;
+
+                //get actual indices
+                i0=(int)x;
+                i1=i0+1;
+                j0=(int)y;
+                j1=j0+1;
+                k0=(int)z;
+                k1=k0+1;
+
+                //calculate percentage of each index
+                s1 = x-i0;
+                s0 = 1-s1;
+                t1 = y-j0;
+                t0 = 1-t1;
+                u1 = z-k0;
+                u0 = 1-u1;
+
+                if(i0 >= N){
+                    i0 = N - 1;
+                }
+                // if(i0 < 0){
+                //     i0 = 0;
+                // }
+                if(j0 >= N){
+                    j0 = N - 1;
+                }
+                // if(j0 < 0){
+                //     j0 = 0;
+                // }
+                if(k0 >= N){
+                    k0 = N - 1;
+                }
+                // if(k0 < 0){
+                //     k0 = 0;
+                // }
+                if(i1 >= N){
+                    i1 = N - 1;
+                }
+                // if(i1 < 0){
+                //     i1 = 0;
+                // }
+                if(j1 >= N){
+                    j1 = N - 1;
+                }
+                // if(j1 < 0){
+                //     j1 = 0;
+                // }
+                if(k1 >= N){
+                    k1 = N - 1;
+                }
+                // if(k1 < 0){
+                //     k1 = 0;
+                // }
+                d[IX(i,j,k)] = 
+                s0*(
+                    t0*u0*sampleArr[IX(i0,j0,k0)]+
+                    t1*u0*sampleArr[IX(i0,j1,k0)]+
+                    t0*u1*sampleArr[IX(i0,j0,k1)]+
+                    t1*u1*sampleArr[IX(i0,j1,k1)]
+                )+
+                s1*(
+                    t0*u0*sampleArr[IX(i1,j0,k0)]+
+                    t1*u0*sampleArr[IX(i1,j1,k0)]+
+                    t0*u1*sampleArr[IX(i1,j0,k1)]+
+                    t1*u1*sampleArr[IX(i1,j1,k1)]
+                );
+            }
+        }
+    }
+    // set_bnd(env, chunk_mask, SET_BOUND_USE_NEIGHBOR, N, b, d, jrd);
+}
+
+JNIEXPORT void JNICALL Java_electrosphere_FluidSim_setBoundsToNeighbors
+  (JNIEnv * env,
+  jobject this,
+  jint N,
+  jint chunk_mask,
+  jint vector_dir,
+  jobjectArray neighborArray){
+    int DIM = N;
+    float * target = GET_ARR(env,neighborArray,CENTER_LOC);
+    for(int x=1; x < DIM-1; x++){
+        for(int y = 1; y < DIM-1; y++){
+            //((x)+(DIM)*(y) + (DIM)*(DIM)*(z))
+            target[IX(0,x,y)] =         vector_dir==BOUND_DIR_U ? -target[IX(1,x,y)] : target[IX(1,x,y)];
+            target[IX(DIM-1,x,y)] =     vector_dir==BOUND_DIR_U ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)];
+            target[IX(x,0,y)] =         vector_dir==BOUND_DIR_V ? -target[IX(x,1,y)] : target[IX(x,1,y)];
+            target[IX(x,DIM-1,y)] =     vector_dir==BOUND_DIR_V ? -target[IX(x,DIM-2,y)] : target[IX(x,DIM-2,y)];
+            target[IX(x,y,0)] =         vector_dir==BOUND_DIR_W ? -target[IX(x,y,1)] : target[IX(x,y,1)];
+            target[IX(x,y,DIM-1)] =     vector_dir==BOUND_DIR_W ? -target[IX(x,y,DIM-2)] : target[IX(x,y,DIM-2)];
+        }
+    }
+    for(int x = 1; x < DIM-1; x++){
+        target[IX(x,0,0)]         = (float)(0.5f * (target[IX(x,1,0)] + target[IX(x,0,1)]));
+        target[IX(x,DIM-1,0)]     = (float)(0.5f * (target[IX(x,DIM-2,0)] + target[IX(x,DIM-1,1)]));
+        target[IX(x,0,DIM-1)]     = (float)(0.5f * (target[IX(x,1,DIM-1)] + target[IX(x,0,DIM-2)]));
+        target[IX(x,DIM-1,DIM-1)] = (float)(0.5f * (target[IX(x,DIM-2,DIM-1)] + target[IX(x,DIM-1,DIM-2)]));
+
+        target[IX(0,x,0)] = (float)(0.5f * (target[IX(1,x,0)] + target[IX(0,x,1)]));
+        target[IX(DIM-1,x,0)] = (float)(0.5f * (target[IX(DIM-2,x,0)] + target[IX(DIM-1,x,1)]));
+        target[IX(0,x,DIM-1)] = (float)(0.5f * (target[IX(1,x,DIM-1)] + target[IX(0,x,DIM-2)]));
+        target[IX(DIM-1,x,DIM-1)] = (float)(0.5f * (target[IX(DIM-2,x,DIM-1)] + target[IX(DIM-1,x,DIM-2)]));
+
+
+        target[IX(0,0,x)] = (float)(0.5f * (target[IX(1,0,x)] + target[IX(0,1,x)]));
+        target[IX(DIM-1,0,x)] = (float)(0.5f * (target[IX(DIM-2,0,x)] + target[IX(DIM-1,1,x)]));
+        target[IX(0,DIM-1,x)] = (float)(0.5f * (target[IX(1,DIM-1,x)] + target[IX(0,DIM-2,x)]));
+        target[IX(DIM-1,DIM-1,x)] = (float)(0.5f * (target[IX(DIM-2,DIM-1,x)] + target[IX(DIM-1,DIM-2,x)]));
 
     }
+    target[IX(0,0,0)]             = (float)((target[IX(1,0,0)]+target[IX(0,1,0)]+target[IX(0,0,1)])/3.0);
+    target[IX(DIM-1,0,0)]         = (float)((target[IX(DIM-2,0,0)]+target[IX(DIM-1,1,0)]+target[IX(DIM-1,0,1)])/3.0);
+    target[IX(0,DIM-1,0)]         = (float)((target[IX(1,DIM-1,0)]+target[IX(0,DIM-2,0)]+target[IX(0,DIM-1,1)])/3.0);
+    target[IX(0,0,DIM-1)]         = (float)((target[IX(0,0,DIM-2)]+target[IX(1,0,DIM-1)]+target[IX(0,1,DIM-1)])/3.0);
+    target[IX(DIM-1,DIM-1,0)]     = (float)((target[IX(DIM-2,DIM-1,0)]+target[IX(DIM-1,DIM-2,0)]+target[IX(DIM-1,DIM-1,1)])/3.0);
+    target[IX(0,DIM-1,DIM-1)]     = (float)((target[IX(1,DIM-1,DIM-1)]+target[IX(0,DIM-2,DIM-1)]+target[IX(0,DIM-1,DIM-2)])/3.0);
+    target[IX(DIM-1,0,DIM-1)]     = (float)((target[IX(DIM-1,0,DIM-2)]+target[IX(DIM-2,0,DIM-1)]+target[IX(DIM-1,1,DIM-1)])/3.0);
+    target[IX(DIM-1,DIM-1,DIM-1)] = (float)((target[IX(DIM-1,DIM-1,DIM-2)]+target[IX(DIM-1,DIM-2,DIM-1)]+target[IX(DIM-1,DIM-1,DIM-2)])/3.0);
+}
+
+void set_bnd(JNIEnv * env, uint32_t chunk_mask, int ignore_neighbor, int N, int b, float * target, jobjectArray neighborMap){
+    int DIM = N;
+    for(int x=1; x < DIM-1; x++){
+        for(int y = 1; y < DIM-1; y++){
+            //((x)+(DIM)*(y) + (DIM)*(DIM)*(z))
+            target[0 + DIM * x + DIM * DIM * y] =         b==1 ? -target[1 + DIM * x + DIM * DIM * y] : target[1 + DIM * x + DIM * DIM * y];
+            target[IX(DIM-1,x,y)] =     b==1 ? -target[IX(DIM-2,x,y)] : target[IX(DIM-2,x,y)];
+            target[IX(x,0,y)] =         b==2 ? -target[IX(x,1,y)] : target[IX(x,1,y)];
+            target[IX(x,DIM-1,y)] =     b==2 ? -target[IX(x,DIM-2,y)] : target[IX(x,DIM-2,y)];
+            target[IX(x,y,0)] =         b==3 ? -target[IX(x,y,1)] : target[IX(x,y,1)];
+            target[IX(x,y,DIM-1)] =     b==3 ? -target[IX(x,y,DIM-2)] : target[IX(x,y,DIM-2)];
+        }
+    }
+    for(int x = 1; x < DIM-1; x++){
+        target[IX(x,0,0)]         = (float)(0.5f * (target[IX(x,1,0)] + target[IX(x,0,1)]));
+        target[IX(x,DIM-1,0)]     = (float)(0.5f * (target[IX(x,DIM-2,0)] + target[IX(x,DIM-1,1)]));
+        target[IX(x,0,DIM-1)]     = (float)(0.5f * (target[IX(x,1,DIM-1)] + target[IX(x,0,DIM-2)]));
+        target[IX(x,DIM-1,DIM-1)] = (float)(0.5f * (target[IX(x,DIM-2,DIM-1)] + target[IX(x,DIM-1,DIM-2)]));
+
+        target[IX(0,x,0)] = (float)(0.5f * (target[IX(1,x,0)] + target[IX(0,x,1)]));
+        target[IX(DIM-1,x,0)] = (float)(0.5f * (target[IX(DIM-2,x,0)] + target[IX(DIM-1,x,1)]));
+        target[IX(0,x,DIM-1)] = (float)(0.5f * (target[IX(1,x,DIM-1)] + target[IX(0,x,DIM-2)]));
+        target[IX(DIM-1,x,DIM-1)] = (float)(0.5f * (target[IX(DIM-2,x,DIM-1)] + target[IX(DIM-1,x,DIM-2)]));
+
+
+        target[IX(0,0,x)] = (float)(0.5f * (target[IX(1,0,x)] + target[IX(0,1,x)]));
+        target[IX(DIM-1,0,x)] = (float)(0.5f * (target[IX(DIM-2,0,x)] + target[IX(DIM-1,1,x)]));
+        target[IX(0,DIM-1,x)] = (float)(0.5f * (target[IX(1,DIM-1,x)] + target[IX(0,DIM-2,x)]));
+        target[IX(DIM-1,DIM-1,x)] = (float)(0.5f * (target[IX(DIM-2,DIM-1,x)] + target[IX(DIM-1,DIM-2,x)]));
+
+    }
+    target[IX(0,0,0)]             = (float)((target[IX(1,0,0)]+target[IX(0,1,0)]+target[IX(0,0,1)])/3.0);
+    target[IX(DIM-1,0,0)]         = (float)((target[IX(DIM-2,0,0)]+target[IX(DIM-1,1,0)]+target[IX(DIM-1,0,1)])/3.0);
+    target[IX(0,DIM-1,0)]         = (float)((target[IX(1,DIM-1,0)]+target[IX(0,DIM-2,0)]+target[IX(0,DIM-1,1)])/3.0);
+    target[IX(0,0,DIM-1)]         = (float)((target[IX(0,0,DIM-2)]+target[IX(1,0,DIM-1)]+target[IX(0,1,DIM-1)])/3.0);
+    target[IX(DIM-1,DIM-1,0)]     = (float)((target[IX(DIM-2,DIM-1,0)]+target[IX(DIM-1,DIM-2,0)]+target[IX(DIM-1,DIM-1,1)])/3.0);
+    target[IX(0,DIM-1,DIM-1)]     = (float)((target[IX(1,DIM-1,DIM-1)]+target[IX(0,DIM-2,DIM-1)]+target[IX(0,DIM-1,DIM-2)])/3.0);
+    target[IX(DIM-1,0,DIM-1)]     = (float)((target[IX(DIM-1,0,DIM-2)]+target[IX(DIM-2,0,DIM-1)]+target[IX(DIM-1,1,DIM-1)])/3.0);
+    target[IX(DIM-1,DIM-1,DIM-1)] = (float)((target[IX(DIM-1,DIM-1,DIM-2)]+target[IX(DIM-1,DIM-2,DIM-1)]+target[IX(DIM-1,DIM-1,DIM-2)])/3.0);
+}

src/main/java/electrosphere/FluidSim.java

@@ -50,11 +50,11 @@ public class FluidSim {
     //Buffers that contain w vector directions
     ByteBuffer[] wVector = new ByteBuffer[27];
     //Buffers that contain u vector directions to add to the simulation
-    ByteBuffer uAdditionVector;
+    ByteBuffer[] uAdditionVector = new ByteBuffer[27];
     //Buffers that contain v vector directions to add to the simulation
-    ByteBuffer vAdditionVector;
+    ByteBuffer[] vAdditionVector = new ByteBuffer[27];
     //Buffers that contain w vector directions to add to the simulation
-    ByteBuffer wAdditionVector;
+    ByteBuffer[] wAdditionVector = new ByteBuffer[27];
 
     //The densities for every voxel for the current frame
     float[] densityArrayView = new float[DIM * DIM * DIM];
@@ -86,18 +86,18 @@ public class FluidSim {
         uVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
         vVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
         wVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-        uAdditionVector = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        uAdditionVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-        vAdditionVector = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        vAdditionVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
-        wAdditionVector = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
+        wAdditionVector[13] = ByteBuffer.allocateDirect(DIM * DIM * DIM * 4);
         //order endian-ness
         density[13].order(ByteOrder.LITTLE_ENDIAN);
         densityAddition.order(ByteOrder.LITTLE_ENDIAN);
         uVector[13].order(ByteOrder.LITTLE_ENDIAN);
         vVector[13].order(ByteOrder.LITTLE_ENDIAN);
         wVector[13].order(ByteOrder.LITTLE_ENDIAN);
-        uAdditionVector.order(ByteOrder.LITTLE_ENDIAN);
+        uAdditionVector[13].order(ByteOrder.LITTLE_ENDIAN);
-        vAdditionVector.order(ByteOrder.LITTLE_ENDIAN);
+        vAdditionVector[13].order(ByteOrder.LITTLE_ENDIAN);
-        wAdditionVector.order(ByteOrder.LITTLE_ENDIAN);
+        wAdditionVector[13].order(ByteOrder.LITTLE_ENDIAN);
 
         //set initial values for chunk
         FloatBuffer xf = density[13].asFloatBuffer();
@@ -148,24 +148,23 @@ public class FluidSim {
             }
         }
 
-        for(int x = 0; x < simArray.length; x++){
+        //
-            for(int y = 0; y < simArray[0].length; y++){
+        //Vector stage
-                for(int z = 0; z < simArray[0][0].length; z++){
-                    simArray[x][y][z].simulate(step,timestep);
-                }
-            }
-        }
-
         solveChunkMask(simArray);
         addVectorSources(simArray, timestep);
+        swapAllVectorFields(simArray, timestep);
         solveVectorDiffusion(simArray, timestep);
-        // setupProjection(simArray, timestep);
         solveProjection(simArray, timestep);
         swapAllVectorFields(simArray, timestep);
         advectVectorsAcrossBoundaries(simArray, timestep);
-        // setupProjection(simArray, timestep);
         solveProjection(simArray, timestep);
 
+        //
+        //Density stage
+        addDensity(simArray, timestep);
+        diffuseDensity(simArray, timestep);
+        advectDensity(simArray, timestep);
+        
 
 
         //
@@ -223,17 +222,12 @@ public class FluidSim {
                     // add_source(N, v, v0, dt);
                     // add_source(N, w, w0, dt);
                     simArray[x][y][z].addSourceToVectorsWrapper(timestep);
-
+                }
+            }
+        }
         //swap
         //u <=> u0 etc for u, v, and w
-                    simArray[x][y][z].swapVectorFields();
         
-                    //setup diffuse for all 3 vectors
-                    simArray[x][y][z].setupVectorDiffuseWrapper(timestep);
-
-                }
-            }
-        }
     }
 
     private static void solveVectorDiffusion(FluidSim[][][] simArray, float timestep){
@@ -249,17 +243,13 @@ public class FluidSim {
                     }
                 }
             }
-        }
-    }
-
-    private static void setupProjection(FluidSim[][][] simArray, float timestep){
             for(int x = 0; x < simArray.length; x++){
                 for(int y = 0; y < simArray[0].length; y++){
                     for(int z = 0; z < simArray[0][0].length; z++){
-                    //then
+                        simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uVector);
-                    //setup vector projection
+                        simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vVector);
-                    //...
+                        simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wVector);
-                    simArray[x][y][z].setupProjectionWrapper(timestep);
+                    }
                 }
             }
         }
@@ -277,29 +267,54 @@ public class FluidSim {
                 }
             }
         }
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].uAdditionVector);
+                    simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].vAdditionVector);
+                }
+            }
+        }
+        //
+        //Perform main projection solver
         for(int l = 0; l < LINEARSOLVERTIMES; l++){
             for(int x = 0; x < simArray.length; x++){
                 for(int y = 0; y < simArray[0].length; y++){
                     for(int z = 0; z < simArray[0][0].length; z++){
                         //lin_solve(env, chunk_mask, N, b, x, x0, a, 1+6*a);
-                        //+
-                        //set_bnd(env, chunk_mask, N, b, x);
                         //for u, v, and w all in 1 shot
                         simArray[x][y][z].solveProjectionWrapper(timestep);
+                    }
+                }
+            }
             //be sure to set boundaries to neighbor chunk values where appropriate
+            for(int x = 0; x < simArray.length; x++){
+                for(int y = 0; y < simArray[0].length; y++){
+                    for(int z = 0; z < simArray[0][0].length; z++){
+                        simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].uAdditionVector);
+                        // simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].vAdditionVector);
                     }
                 }
             }
         }
+        //Finalize projection
         for(int x = 0; x < simArray.length; x++){
             for(int y = 0; y < simArray[0].length; y++){
                 for(int z = 0; z < simArray[0][0].length; z++){
                     //Subtract curl field from current vector field
                     //...
-
+                    simArray[x][y][z].finalizeProjectionWrapper(timestep);
+                }
+            }
+        }
         //set boundaries a final time for u,v,w
         //...
-                    simArray[x][y][z].finalizeProjectionWrapper(timestep);
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uVector);
+                    simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vVector);
+                    simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wVector);
                 }
             }
         }
@@ -316,9 +331,23 @@ public class FluidSim {
     }
 
     private void swapVectorFields(){
+        ByteBuffer tmp;
+        //swap x0 <-> x
+        tmp = densityAddition;
+        densityAddition = density[13];
+        density[13] = tmp;
         //swap u0 <-> u
+        tmp = uAdditionVector[13];
+        uAdditionVector[13] = uVector[13];
+        uVector[13] = tmp;
         //swap v0 <-> v
+        tmp = vAdditionVector[13];
+        vAdditionVector[13] = vVector[13];
+        vVector[13] = tmp;
         //swap w0 <-> w
+        tmp = wAdditionVector[13];
+        wAdditionVector[13] = wVector[13];
+        wVector[13] = tmp;
         //...
     }
 
@@ -334,6 +363,16 @@ public class FluidSim {
                 }
             }
         }
+        //mirror neighbor data
+        for(int x = 0; x < simArray.length; x++){
+            for(int y = 0; y < simArray[0].length; y++){
+                for(int z = 0; z < simArray[0][0].length; z++){
+                    simArray[x][y][z].setBoundsToNeighborsWrapper(1, simArray[x][y][z].uVector);
+                    simArray[x][y][z].setBoundsToNeighborsWrapper(2, simArray[x][y][z].vVector);
+                    simArray[x][y][z].setBoundsToNeighborsWrapper(3, simArray[x][y][z].wVector);
+                }
+            }
+        }
     }
 
     private static void addDensity(FluidSim[][][] simArray, float timestep){
@@ -345,8 +384,6 @@ public class FluidSim {
                     //swap x <=> x0
                     //swap arrays in java side...
                     simArray[x][y][z].swapDensityArrays();
-                    //setup density diffusion
-                    simArray[x][y][z].setupDiffuseDensityWrapper(timestep);
                 }
             }
         }
@@ -359,6 +396,7 @@ public class FluidSim {
     }
 
     private static void diffuseDensity(FluidSim[][][] simArray, float timestep){
+        for(int l = 0; l < LINEARSOLVERTIMES; l++){
             for(int x = 0; x < simArray.length; x++){
                 for(int y = 0; y < simArray[0].length; y++){
                     for(int z = 0; z < simArray[0][0].length; z++){
@@ -369,6 +407,14 @@ public class FluidSim {
                     }
                 }
             }
+            for(int x = 0; x < simArray.length; x++){
+                for(int y = 0; y < simArray[0].length; y++){
+                    for(int z = 0; z < simArray[0][0].length; z++){
+                        simArray[x][y][z].setBoundsToNeighborsWrapper(0, simArray[x][y][z].density);
+                    }
+                }
+            }
+        }
     }
 
     private static void advectDensity(FluidSim[][][] simArray, float timestep){
@@ -444,9 +490,9 @@ public class FluidSim {
         ByteBuffer[] u,
         ByteBuffer[] v,
         ByteBuffer[] w,
-        ByteBuffer u0,
+        ByteBuffer[] u0,
-        ByteBuffer v0,
+        ByteBuffer[] v0,
-        ByteBuffer w0,
+        ByteBuffer[] w0,
         float DIFFUSION_CONSTANT,
         float VISCOSITY_CONSTANT,
         float timestep
@@ -468,15 +514,7 @@ public class FluidSim {
     private void addSourceToVectorsWrapper(float timestep){
         addSourceToVectors(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
     }
-    private native void addSourceToVectors(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0, ByteBuffer v0, ByteBuffer w0, float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+    private native void addSourceToVectors(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
-
-    /**
-     * Setup u, v, and w to diffuse
-     */
-    private void setupVectorDiffuseWrapper(float timestep){
-        setupVectorDiffuse(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
-    }
-    private native void setupVectorDiffuse(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0, ByteBuffer v0, ByteBuffer w0, float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
 
     /**
      * Solves u, v, and w diffusion systems of equations
@@ -484,7 +522,7 @@ public class FluidSim {
     private void solveVectorDiffuseWrapper(float timestep){
         solveVectorDiffuse(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
     }
-    private native void solveVectorDiffuse(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0, ByteBuffer v0, ByteBuffer w0, float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+    private native void solveVectorDiffuse(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
 
     /**
      * Setup projection system
@@ -492,7 +530,7 @@ public class FluidSim {
     private void setupProjectionWrapper(float timestep){
         setupProjection(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
     }
-    private native void setupProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0, ByteBuffer v0, ByteBuffer w0, float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+    private native void setupProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
 
     /**
      * Solve projection system
@@ -500,7 +538,15 @@ public class FluidSim {
     private void solveProjectionWrapper(float timestep){
         solveProjection(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
     }
-    private native void solveProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0, ByteBuffer v0, ByteBuffer w0, float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+    private native void solveProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+
+    /**
+     * Solve projection system
+     */
+    // private void setProjectionBordersWrapper(float timestep){
+    //     setProjectionBorders(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
+    // }
+    // private native void setProjectionBorders(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
 
     /**
      * Does work like subtracting curl from vector field, setting boundaries, etc
@@ -508,7 +554,7 @@ public class FluidSim {
     private void finalizeProjectionWrapper(float timestep){
         finalizeProjection(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
     }
-    private native void finalizeProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0, ByteBuffer v0, ByteBuffer w0, float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+    private native void finalizeProjection(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
 
     /**
      * Advects vectors
@@ -516,7 +562,7 @@ public class FluidSim {
     private void advectVectorsWrapper(float timestep){
         advectVectors(DIM, chunkMask, uVector, vVector, wVector, uAdditionVector, vAdditionVector, wAdditionVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
     }
-    private native void advectVectors(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0, ByteBuffer v0, ByteBuffer w0, float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
+    private native void advectVectors(int DIM_X, int chunkMask, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], ByteBuffer u0[], ByteBuffer v0[], ByteBuffer w0[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
 
     /**
      * Adds density to the simulation
@@ -526,14 +572,6 @@ public class FluidSim {
     }
     private native void addDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer x0, float timestep);
 
-    /**
-     * Solve density diffusion
-     */
-    private void setupDiffuseDensityWrapper(float timestep){
-        setupDiffuseDensity(DIM, chunkMask, density, densityAddition, uVector, vVector, wVector, DIFFUSION_CONSTANT, VISCOSITY_CONSTANT, timestep);
-    }
-    private native void setupDiffuseDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer x0, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
-
     /**
      * Solve density diffusion
      */
@@ -551,6 +589,15 @@ public class FluidSim {
     private native void advectDensity(int DIM_X, int chunkMask, ByteBuffer[] x, ByteBuffer x0, ByteBuffer[] u, ByteBuffer v[], ByteBuffer w[], float DIFFUSION_CONSTANT, float VISCOSITY_CONSTANT, float timestep);
 
 
+    /**
+     * Sets the bounds of the neighbormap to neighbor values if available
+     */
+    private void setBoundsToNeighborsWrapper(int vectorDir, ByteBuffer[] neighborMap){
+        setBoundsToNeighbors(DIM, chunkMask, vectorDir, neighborMap);
+    }
+    private native void setBoundsToNeighbors(int DIM_X, int chunkMask, int vectorDir, ByteBuffer[] neighborMap);
+
+
 
 
 
@@ -644,19 +691,19 @@ public class FluidSim {
         if(densityAddition.position() > 0){
             densityAddition.position(0);
         }
-        if(uAdditionVector.position() > 0){
+        if(uAdditionVector[13].position() > 0){
-            uAdditionVector.position(0);
+            uAdditionVector[13].position(0);
         }
-        if(vAdditionVector.position() > 0){
+        if(vAdditionVector[13].position() > 0){
-            vAdditionVector.position(0);
+            vAdditionVector[13].position(0);
         }
-        if(wAdditionVector.position() > 0){
+        if(wAdditionVector[13].position() > 0){
-            wAdditionVector.position(0);
+            wAdditionVector[13].position(0);
         }
         FloatBuffer x0FloatView = densityAddition.asFloatBuffer();
-        FloatBuffer u0FloatView = uAdditionVector.asFloatBuffer();
+        FloatBuffer u0FloatView = uAdditionVector[13].asFloatBuffer();
-        FloatBuffer v0FloatView = vAdditionVector.asFloatBuffer();
+        FloatBuffer v0FloatView = vAdditionVector[13].asFloatBuffer();
-        FloatBuffer w0FloatView = wAdditionVector.asFloatBuffer();
+        FloatBuffer w0FloatView = wAdditionVector[13].asFloatBuffer();
         for(int i = 0; i < DIM; i++){
             for(int j = 0; j < DIM; j++){
                 for(int k = 0; k < DIM; k++){
@@ -753,11 +800,26 @@ public class FluidSim {
         return wVector[getNeighborIndex(1,1,1)];
     }
 
+    public ByteBuffer getUAddBuffer(){
+        return uAdditionVector[getNeighborIndex(1,1,1)];
+    }
+
+    public ByteBuffer getVAddBuffer(){
+        return vAdditionVector[getNeighborIndex(1,1,1)];
+    }
+
+    public ByteBuffer getWAddBuffer(){
+        return wAdditionVector[getNeighborIndex(1,1,1)];
+    }
+
     public void setNeighbor(int x, int y, int z, FluidSim neighbor){
         density[getNeighborIndex(x,y,z)] = neighbor.getDensityBuffer();
         uVector[getNeighborIndex(x,y,z)] = neighbor.getUBuffer();
         vVector[getNeighborIndex(x,y,z)] = neighbor.getVBuffer();
         wVector[getNeighborIndex(x,y,z)] = neighbor.getWBuffer();
+        uAdditionVector[getNeighborIndex(x,y,z)] = neighbor.getUAddBuffer();
+        vAdditionVector[getNeighborIndex(x,y,z)] = neighbor.getVAddBuffer();
+        wAdditionVector[getNeighborIndex(x,y,z)] = neighbor.getWAddBuffer();
     }
 
 }